基于QAA算法的遼河口懸浮物濃度遙感反演*

鄧智天 孫永華 邱 琦 孫 薇 倪 萍 邢 瑞

(首都師范大學資源環(huán)境與旅游學院，北京 100048)

0 引 言

懸浮物濃度(suspended sediment concentration, SSC)是海水最重要的參數(shù)之一[1].不同濃度的懸浮物直接導致光散射、折射和吸收的差異，最終導致水的光學性質(zhì)的差異.超過一定濃度的懸浮物會使海水渾濁，水的透明度降低，從而導致入射到水中的太陽輻射減弱，海水中浮游生物吸收的能量減少，使得光合作用被削弱，并最終影響海水中浮游生物的生長.懸浮物對水質(zhì)的影響也較大，原因在于懸浮物包括一定數(shù)量的膠體物質(zhì)和黏土礦物.一方面，懸浮物具有一定的吸收能力，可以吸收污染物，起到凈化水的作用；另一方面，它成為污染物的遷移和再循環(huán)的重要依附物.同時，懸浮物也是研究泥沙輸移和地貌演化的重要依據(jù)，可以幫助預測沉積物的承載能力，并在沿海水域泥沙輸移規(guī)律的研究中發(fā)揮重要作用.其研究不僅可以探索海水污染的情況，還可以進一步追蹤污染物的流動路徑和海水中污染物的方向，然后探討海水中污染物的運動路徑和變化規(guī)律，為海水污染控制提供技術(shù)支持.

遙感可以有效地克服現(xiàn)場測量的不足，具有范圍廣、效率高的優(yōu)點，成為了研究海水中懸浮物的主要技術(shù)手段.20世紀90年代初，Ritchie和Cooper[2]和Topliss等[3]使用Landsat 3 多光譜掃描儀(multispectral scanner, MSS)進行了水體懸浮物的濃度反演，指出MSS對懸浮物濃度有相關(guān)性；Landsat 5在軌穩(wěn)定運行之后，很多研究者對專題繪圖儀(thematic mapper，TM)數(shù)據(jù)進行了水體懸浮物的研究[4-12]，使用了單波段、波段組合等方法構(gòu)建模型，對水體懸浮物進行了定量研究；還有學者使用了MODIS數(shù)據(jù)對更大范圍的湖泊、海洋等進行了研究[13-19]，由于其波段較多且數(shù)據(jù)較多，包括卡爾曼濾波、機器學習等方法相繼應用于水體懸浮物的定量研究；相關(guān)學者分別對海景寬視場傳感器[20]、機載高光譜掃描儀[21]、小型機載光譜成像儀[22]、中分辨率成像光譜儀[23]、伊科諾斯[24]等水色遙感、高光譜遙感、高分遙感傳感器進行懸浮物反演，沿用了MODIS和Landsat的反演方法，并使用了更為專業(yè)的傳感器，因此在反演精度上有所提高.

QAA由Lee等[25]在2002年提出，用于推導出深水光學的固有光學特性.固有光學特性(inherent optical properties，IOP)是指僅與水體自身組成成分相關(guān)，由介質(zhì)決定，不隨光照條件改變而改變的水體光學參量，主要包含海水的吸收系數(shù)、散射系數(shù)和衰減系數(shù)等[26]，其相關(guān)研究集中在IOP與懸浮物顆粒的特性[27].國際海洋水色遙感組織在2006年提交了IOP算法報告[28]，在該報告中，表明QAA算法對于清澈海水和沿海水域都有非常高適用性.直到2018年，QAA已經(jīng)歷了5次更新，最新正式發(fā)布的版本是QAAv6[29].這些更新使QAA更加完美，性能也大大提高.在Lee之外，許多學者對算法進行了改進，如Chen和Zhang[30]基于MODIS數(shù)據(jù)，提出了一種改進的算法QAA-RGR, 直接選擇衛(wèi)星遙感數(shù)據(jù)而非實測高光譜，證明衛(wèi)星遙感數(shù)據(jù)在一定條件下能替代實測高光譜進行QAA算法的計算，并能取得較高的精度，這也開啟了近幾年使用不同的遙感數(shù)據(jù)關(guān)于QAA算法的研究.

近年來，基于QAA算法懸浮物的研究有所增加.陳莉瓊[31]使用實測高光譜進行研究，指出QAA模型比經(jīng)驗模型具有更高的精度和時相穩(wěn)定性;王建國等[32]、段化杰等[33]分別使用MODIS和Sentinel 3 OLCI數(shù)據(jù)對懸浮物濃度進行了相關(guān)研究，發(fā)現(xiàn)其精度與實測高光譜相當.

通過QAA算法，可以實現(xiàn)由遙感反射率數(shù)據(jù)推導IOP，既減少了IOP數(shù)據(jù)獲取的成本，具有時效性與經(jīng)濟性，又可以對遙感反射率數(shù)據(jù)進行深層次的數(shù)據(jù)挖掘，而非局限于反射率數(shù)值.而研究衛(wèi)星多光譜數(shù)據(jù)QAA算法適用性，探討其能否在一定條件下替代現(xiàn)場實測高光譜數(shù)據(jù)，可以進一步放大QAA算法的時效性與經(jīng)濟性.

本研究將以含沙量較高的遼河口三角洲為例，基于QAA算法，使用天宮二號寬波段成像儀數(shù)據(jù)，選取682.5 nm為參考波長，求出水體總吸收系數(shù)a(λ)，利用2018年6月7—9日實測水體懸浮物濃度數(shù)據(jù)，建立基于a(λ)的懸浮物濃度反演模型.

1 研究區(qū)和數(shù)據(jù)源

1.1 研究區(qū)及采樣位置

遼河平均徑流量95.27億m3，河道彎曲，呈不規(guī)則河型，水系發(fā)育，大小支流70余條，中下游河道寬淺，河道寬1 000～2 000 m；鐵嶺水文站多年平均含沙量3.60 kg/m3，多年平均輸沙量2 098萬噸；水流緩慢，泥沙淤積；河床質(zhì)為沙壤土.治理水土流失和開展水土保持也是遼河流域的一個主要問題，50年代初遼河流域水土流失面積約1 220萬hm2，占流域面積的55.7%；1985年全流域水土流失面積為951萬hm2，占全流域水土流失面積的45.2%；其中老哈、教來河位于冀北遼西山地和黃土丘陵區(qū)，植被覆蓋率不到30%，水土流失十分嚴重.

樣本采集于6月7—9日沿著雙臺子河和大遼河進行，如圖1所示，圖例中綠色三角形為采樣位置，底圖為全球海圖數(shù)據(jù)，兩條河流屬于遼河支流.本研究采樣方案著重在河道中間深水處，避開灘涂，即圖中綠色部分(水深≤0 m)；由于海船體積較大，吃水深，行駛時可能會帶動底層沉積物及攪渾海水，本研究采樣時采用線性采樣，沿著行駛的路線等間距采樣，并在駛?cè)氩蓸訁^(qū)時逐漸減速，以減少底層沉積物帶來的誤差；北邊采樣重點在雙臺子河的懸浮物濃度；中間部分采樣重點在河口三角洲灘涂附近的懸浮物濃度，并在大潮時(水深>1.5 m)對灘涂正上方進行了數(shù)次采樣；東邊采樣重點在大遼河的懸浮物濃度.本次樣本采集一共采集104個水樣進行測量.

1.2 懸浮物濃度測定

試驗標準參照中華人民共和國生態(tài)環(huán)境部確定懸浮物濃度標準《水質(zhì) 懸浮物的測定 重量法：GB 11901—89》[34].該標準指出水質(zhì)中的懸浮物是指水樣通過孔徑為0.45 μm的濾膜，截留在濾膜上并于103～105 ℃烘干至恒重的固體物質(zhì).

(1)用扁咀無齒鑷子夾取微孔濾膜放于事先恒重的稱量瓶里，移入烘箱中于103～105 ℃下烘烤半小時，烘干半小時后取出置干燥器內(nèi)冷卻至室溫，稱其重量.反復烘干、冷卻、稱量，直至兩次稱量的重量差<0.2 mg；

(2)量取充分混合均勻的試樣100 mL抽吸過濾；使水分全部通過濾膜；再以每次10 mL蒸餾水連續(xù)洗滌三次，繼續(xù)吸濾以除去痕量水分；

(3)停止吸濾后，仔細取出載有懸浮物的濾膜放在原恒重的稱量瓶里，移入烘箱中于103～105 ℃下烘干1小時后移入干燥器中，使冷卻到室溫，稱其重量；反復烘干、冷卻、稱量，直至2次稱量的重量差<0.4 mg為止.

圖1 研究區(qū)及采樣位置示意圖

1.3 遙感數(shù)據(jù)處理

天宮二號空間實驗室，是繼天宮一號任務(wù)完成其使命后，發(fā)射的第二個空間實驗室，安排了一批體現(xiàn)科學前沿和戰(zhàn)略高技術(shù)發(fā)展方向的科學與應用任務(wù)，開展相應的應用與試驗，對相關(guān)新技術(shù)進行體制驗證.天宮二號空間實驗室發(fā)射后，與神舟十一號載人飛船對接，進行人在太空的中期駐留實驗；與天舟一號貨運飛船進行對接，開展推進劑補給等相關(guān)試驗.

天宮二號空間實驗室搭載了全新的空間應用載荷設(shè)備，載荷數(shù)量及規(guī)模都超過了以往我國各次載人航天任務(wù).開展10余項應用與實驗，涉及對地觀測和空間地球科學、空間天文、微重力基礎(chǔ)物理、微重力流體物理及空間材料科學、空間生命科學和空間環(huán)境與空間物理等多個領(lǐng)域.其中對地觀測和空間地球科學有寬波段成像儀、三維成像微波高度計和紫外臨邊成像光譜儀，其中寬波段成像儀的波段設(shè)置如表1.

表1 天宮二號寬波段成像儀波段設(shè)置

本研究使用影像獲取時間與采樣時間間隔最小的一幅TG-2數(shù)據(jù)，2018年6月1日的數(shù)據(jù)，由于2018年6月初無強降雨臺風等天氣，故認為2018年6月1日的影像能反映懸浮物在2018年6月7—9日采樣時的實際分布.

TG-2 MWI Level-2數(shù)據(jù)通過載人航天空間應用數(shù)據(jù)推廣服務(wù)平臺申請數(shù)據(jù)產(chǎn)品(http：//www.msadc.cn/)，該Level-2數(shù)據(jù)產(chǎn)品已經(jīng)經(jīng)過幾何校正和輻射糾正.Level-2數(shù)據(jù)包含表觀反射率的校準數(shù)據(jù)，圖像的表觀反射率由ENVI中的大氣校正模塊Quick Atmospheric Correction(QUAC)生成.

2 方 法

2.1 波長λ0的獲取

為去除水面折射對反射率的影響，通過公式(1)將表觀反射率Rrs(λ)換算為水面以下反射率rrs(λ):

(1)

后向散射與吸收系數(shù)比值u(λ)和總吸收系數(shù)a(λ)及后向散射系數(shù)bb(λ)的定量關(guān)系如公式(2):

(2)

其中，g0=0.089，g1=0.1245.

表2 QAAv6算法

懸浮粒子后向散射系數(shù)bbp(λ)和波長λ滿足冪函數(shù)規(guī)律，如公式(3):

(3)

其中λ0代表參考波長.

在海洋中，純水后向散射bbw(λ)占bb(λ)的1/10[35]，如公式(4):

(4)

根據(jù)懸浮物反射特征，700 nm處的吸收系數(shù)迅速下降，我們將682.5 nm指定為λ0，即VNI Band 7的中心波長.

2.2 QAAv6算法

QAAv6算法步驟見表2，步驟①將水面的表觀反射率轉(zhuǎn)換到水面以下的反射率；步驟②計算后向散射與吸收系數(shù)比值，經(jīng)過查找天宮數(shù)據(jù)元文件，Rrs(670)>0.001 5sr-1，執(zhí)行右邊的公式；步驟③～④為參考波長的固有光學參數(shù)估算，其中純水吸收系數(shù)aw(670)根據(jù)Deng等[36]所做的室內(nèi)純水吸收系數(shù)表格查找得出；步驟⑤～⑥為任意波長的固有光學參數(shù)估算，通過計算的λ0與λ的系數(shù)，估算出任意波長的后向散射系數(shù)bbp(λ)；步驟⑦高光譜波長的總吸收系數(shù)估算，通過前面計算得出的u(λ)、bbw(λ)、bbp(λ)，估算總吸收系數(shù)a(λ).本研究使用線性回歸模型建立a(λ)和SSC之間的反演模型.

本研究構(gòu)建懸浮物濃度反演模型主要有2個部分：第一部分是輸入遙感反射率為，利用QAA估算水體中總吸收系數(shù)a(λ)；第二部分是根據(jù)QAA算法估算的水體中總吸收系數(shù)a(λ)構(gòu)建波段歸一化因子與懸浮物濃度值建立回歸模型，進而求解水體懸浮物濃度.

3 結(jié) 果

3.1 懸浮物濃度結(jié)果

懸浮物濃度測定結(jié)果見圖2，大部分樣品濃度位于18～33 mg/L的區(qū)間內(nèi)，其中雙臺子河口采樣區(qū)域43個點均值為23.91 mg/L，大遼河口采樣區(qū)域35個點均值為24.66 mg/L，河口三角洲采樣區(qū)域26個點均值為26.08 mg/L；有4個離群點，應為樣本保存或?qū)嶒炦^程中的誤差，本研究將該4個離群值去除.

圖2 懸浮物濃度測定結(jié)果

3.2 光譜分析

除去異常點，在該研究中使用100個點；隨機選擇70個點進行回歸分析，其余30個點用于驗證.如圖3所示，原始光譜曲線在700 nm附近有一個反射峰，這與懸浮物的光譜特性一致.

圖3 樣本表觀反射率Rrs(λ)

bbp(λ)與Rrs(λ)和a(λ)相比，由于計算方法的特性，僅需要一個參考點λ0便可以計算出所需的400～1 000 nm內(nèi)所有波段，如圖4所示，其光譜特征具有一致性，后向散射系數(shù)隨著波長的增加呈冪函數(shù)遞減.

圖4 樣本后向散射系數(shù)bbp(λ)

在QAA算法模型的第⑦步之后，得到了總樣本的吸收系數(shù)，結(jié)果如圖5所示.此處由于涉及到原始多光譜的相關(guān)波段，所以估算出來每一個樣品的a(λ)只有14個點，后文a(565)代表著565 nm處的a(λ).

圖5 樣本總吸收系數(shù)a(λ)

3.3 懸浮物濃度反演模型構(gòu)建

在構(gòu)建模型時，使用顯著性檢驗來輔助波段選擇，如圖6所示，選擇超過顯著性水平0.01的波段，即圖中黑線外的點.

圖6 a(λ)與SSC的相關(guān)系數(shù)

根據(jù)懸浮物的反射特性，本研究選取Band9和Band10波段的水體吸收系數(shù)a(620)和a(565)構(gòu)建波段歸一化因子，并將該歸一化因子與懸浮物濃度值建立線性回歸方程：

(5)

3.4 模型驗證

本研究將驗證組的30個點的相關(guān)數(shù)據(jù)，通過QAA模型進行演算，對比其模擬值與實測值，評估總吸收系數(shù)的反演效果并分析QAA模型準確性.本研究使用均方根誤差(average relative error，簡稱RMSE)[37]，平均相對誤差(average relative error，ARE)，如公式(6)和(7)，以及相關(guān)系數(shù)(R2)來評估.

(6)

(7)

QAA模型的準確度R2為 0.685 9，RMSE為2.814 8，ARE為9.68%.通常，QAA估算得到的a(λ)與綠色(500～550 nm)和紅色(600～650 nm)波段一致，但在其他波段，尤其是藍色波段(400～500 nm)，估算結(jié)果較差且明顯被低估，而在550～600 nm，則有過高估算[38].對于大多數(shù)內(nèi)陸渾濁水域，浮游植物色素吸收系數(shù)aφ(λ)和黃色物質(zhì)吸收系數(shù)ag(λ)是相互獨立，而海洋或沿海水樣中的水質(zhì)參數(shù)濃度或吸收系數(shù)相對低于內(nèi)陸渾水，aφ(λ)和ag(λ)之間存在相關(guān)關(guān)系.以往的研究表明，ag(λ)在水中呈現(xiàn)出指數(shù)衰減規(guī)律，導致短波段遙感反射中出現(xiàn)大量復雜的水色信息，從而增加了短波段算法的不確定性.

本研究利用2018年6月1日的TG-2 MWI數(shù)據(jù)，基于QAA算法，選擇682.5 nm作為參考波長λ0，建立了適用于了河口的懸浮物濃度反演模型，將該模型應用于TG-2 WIS數(shù)據(jù)，其結(jié)果如圖7所示.

圖7 遼河三角洲懸浮物濃度空間分布

從圖7可知，懸浮物濃度呈現(xiàn)自北向南逐漸降低的趨勢，濃度范圍11～32 mg/L，雙臺子河、大遼河河口濃度都大于 20 mg/L，與實測值相近；碼頭靜止的水的懸浮物濃度較碼頭外低了10 mg/L 左右，符合實際情況.

第一個誤差是由于大多數(shù)采樣點靠近岸邊，其濃度值高于研究區(qū)域的其他部分，導致反演結(jié)果相對集中在高值.第二個誤差在于QAA算法是根據(jù)水的輻射度逐步建立光學和生化特征之間的關(guān)系，其計算步驟容易傳輸并累積誤差.第三個誤差主要來自水面以上的遙感反射轉(zhuǎn)換過程，計算無法直接獲得的水下反射率rrs(λ)，目前，對不同的測量角度進行精確建模需要考慮風速、陰影等，而本方法使用經(jīng)驗參數(shù)來替代.此外，即使排除河床底部反射并假設(shè)其為零誤差，QAA算法的計算仍存在許多不確定性.如在QAA算法的步驟1中，g0和g1不是常數(shù)，但受到太陽天頂角和水體光學特性的影響，并且隨水體散射特性而變化，所以參考波長處的吸收系數(shù)也是QAA算法中的主要誤差源之一.

4 結(jié) 論

本研究以含沙量較高的遼河口三角洲為例，利用2018年6月7—9日實測水體懸浮物濃度數(shù)據(jù)，基于QAA算法，使用天宮二號寬波段成像儀數(shù)據(jù)，選取682.5 nm為參考波長，求出水體總吸收系數(shù)a(λ)，建立基于a(λ)的懸浮物濃度反演模型，得到2018年6月遼河三角洲懸浮物濃度空間分布圖.主要結(jié)論包括：

(1)根據(jù)懸浮物反射特征，700 nm處的吸收系數(shù)最低，將682.5 nm指定為λ0是合適的；(2)基于遙感反射率與吸收系數(shù)推出的總吸收系數(shù)a(λ)-懸浮物模型的擬合精度R2為0.685 9，能較為正確的反映遼河三角洲懸浮物濃度變化趨勢，但是精度有待提高；(3)QAA算法強大但復雜，為了減少其誤差，需要在采樣時留意距岸距離、天氣等因素，才能更好發(fā)揮其作用.

本研究探討了天宮二號寬波段成像儀數(shù)據(jù)通過QAA模型演算IOP的適用性，通過將遙感反射率轉(zhuǎn)化為總吸收系數(shù)，研究總吸收系數(shù)與懸浮物濃度的相關(guān)性，選取最佳波段，總結(jié)反演模型.受限于作者水平，本研究IOP演算止于總吸收系數(shù)，如果能從總吸收系數(shù)中分離出懸浮物單位吸收系數(shù)會有助于提高反演模型的精度.