基于Landsat 8的長(zhǎng)江口懸浮泥沙濃度遙感反演

馬 驊，鐘煌亮，羅 章，朱 磊，施 峰

（中交疏浚技術(shù)裝備國(guó)家工程研究中心有限公司，上海 201208）

相比傳統(tǒng)的數(shù)據(jù)采集方式而言，遙感數(shù)據(jù)具有空間覆蓋、時(shí)間周期覆蓋以及歷史數(shù)據(jù)再現(xiàn)等方面的優(yōu)勢(shì)，能為大面積水域監(jiān)測(cè)提供有效途徑，還能發(fā)現(xiàn)常規(guī)方法難以發(fā)現(xiàn)的物質(zhì)時(shí)空變化特征。雖然水色遙感監(jiān)測(cè)的方法還不能完全替代精確的實(shí)地測(cè)量方法，但是在宏觀監(jiān)測(cè)、評(píng)測(cè)、管理方面還是具有很大的優(yōu)勢(shì)。海水按照其懸浮物的來源和物質(zhì)組成的不同分為一類水體和二類水體，一類水體主要是大洋開闊水體，二類水體主要是近岸河口水體。從1978年第一顆水色傳感器CACS（Coastal Zone Color Scanner）到各類水色傳感器的出現(xiàn)，很多學(xué)者借助實(shí)測(cè)數(shù)據(jù)和經(jīng)驗(yàn)知識(shí)為不同區(qū)域或者特定時(shí)段，建立了不同的反演模式。這類算法在一類水體中精度較高，很多成果已經(jīng)得到業(yè)務(wù)化應(yīng)用，并取得了良好的效果[1]。一類水體的光學(xué)特性主要受到浮游植物和降解物質(zhì)色素的影響，但是在研究二類水體光學(xué)特性時(shí)，必須考慮當(dāng)?shù)厮w中葉綠素、總懸浮顆粒物和有色可溶性有機(jī)物（CDOM）三要素對(duì)離水光譜輻射的影響[2-3]。像長(zhǎng)江河口這類大河口地區(qū)，受到潮汐和徑流交匯的影響，其泥沙含量高，水體渾濁，水動(dòng)力情況復(fù)雜且存在周期性變化，因此，其水體光學(xué)特性更加復(fù)雜多變。雖然已有考慮不同含沙量水體的光譜特性[4]，也有討論水體中不同泥沙顆粒物粒徑對(duì)反射率的影響[5]，但是至今仍沒有理想的模型能夠?qū)崿F(xiàn)穩(wěn)定監(jiān)測(cè)，而且，自上游對(duì)干流徑流調(diào)蓄以后，進(jìn)入河口的水沙通量和季節(jié)分配與之前有很大不同[6]，給河口水體的光學(xué)特性可能帶來影響。在杭州灣水域利用MERIS數(shù)據(jù)曾建立了復(fù)雜反演算法，取得了良好效果[7]。Landsat 8衛(wèi)星提供了4個(gè)可見光和1個(gè)近紅外波段，本項(xiàng)研究嘗試針對(duì)Landsat 8波段設(shè)置，在長(zhǎng)江口地區(qū)建立復(fù)雜經(jīng)驗(yàn)算法。在測(cè)得表觀光學(xué)特征的基礎(chǔ)上，分析懸浮泥沙在潮周期內(nèi)的變化，并討論各波段對(duì)TSM濃度的敏感性，對(duì)今后建立本地化反演算法、提高遙感反演算法精度打下基礎(chǔ)。

1 數(shù)據(jù)與測(cè)量方法

本文研究區(qū)域?yàn)楸备酆湍喜?，見圖1，該區(qū)域在一個(gè)大小潮周期內(nèi)受到潮汐作用，三要素變化明顯，導(dǎo)致表觀光學(xué)變化。2014年5～6月南槽和北港的連續(xù)觀測(cè)數(shù)據(jù)，兩次現(xiàn)場(chǎng)測(cè)量共獲得67組水沙數(shù)據(jù)和41組葉綠素、CDOM和表面光學(xué)同步數(shù)據(jù)。由于攔門沙是涇潮流相互作用的平衡地帶，其規(guī)模和形態(tài)取決于河流徑流量、含沙量的大小和潮汐作用的強(qiáng)弱，又是河口地區(qū)的主要航行通道，一般情況下也是河口地帶懸浮泥沙含量最高的區(qū)域，因此選擇攔門沙作為研究站點(diǎn)具有重要的理論價(jià)值和實(shí)際意義。

水的表面光學(xué)量測(cè)量使用Satlantic公司的Hyper-SAS高光譜數(shù)字光譜儀，采用水面以上測(cè)量法[3]，在避開太陽(yáng)直射、忽略或者避免白帽現(xiàn)象的情況下，保持儀器的觀測(cè)平面與太陽(yáng)入射平面的夾角在之間，向下總?cè)肷涮筋^與海面法線方向保持在之間，Rrs可由(1)式計(jì)算而得[8]：

圖1 研究區(qū)域與觀測(cè)站點(diǎn)Fig.1 The research area and fixed observation site

(1)式中，Ls,dif是天空漫散射光，是干擾水體光譜的噪聲信息，必須去除；ρsky(W)Ls,dif是氣-水面對(duì)天空漫散射的反射率，它取決于太陽(yáng)高度角等諸多因素，當(dāng)風(fēng)速較?。ā?m/s），天空無云或云量較少時(shí)，可據(jù)Ruddick等于2006年提出的式(2)來確定：

本文取0.0256，由此可以由(3)式計(jì)算出離水輻亮度為：

Ed(λ)是經(jīng)過HySas自帶軟件計(jì)算而得的海面總?cè)肷漭椪斩?，因此，根?jù)遙感反射率的定義，Rrs(λ,Θ0,W)可以由上式計(jì)算得到。

現(xiàn)場(chǎng)采集600ml表層水樣，于每天8:00-17:00，每整點(diǎn)取水樣600ml，特征潮情時(shí)刻（漲憩、落憩、漲急、落急）和光譜測(cè)量時(shí)刻加采，運(yùn)回實(shí)驗(yàn)室采用過濾稱量法測(cè)量懸浮泥沙濃度。使用孔徑63μm玻璃纖維濾膜，對(duì)含沙水樣進(jìn)行抽濾，然后放入105℃恒溫箱內(nèi)烘干，再放入干燥器中冷卻10min后，放入1/10000分析天平稱重。水體中懸沙濃度值等于泥沙和濾膜的總重減去濾膜重量再除以過濾水樣體積。

遙感影像的處理流程主要分為輻射定標(biāo)和大氣校正，輻射定標(biāo)使用ENVI的輻射定標(biāo)模塊，大氣校正使用6S大氣輻射傳輸模型，由于6S模型依賴于輸入的參數(shù)，其中大氣模型、氣溶膠類型和550nm處的光學(xué)厚度都能影響計(jì)算出的糾正公式，從而影響大氣校正的精度，本文所采用的氣溶膠含量按照NASA當(dāng)天的反演AOD值取平均值。表1為2014年6月13日OLI傳感器大氣校正參數(shù)設(shè)置。

2 分析與建模

2.1 光譜特征分析

光譜分析對(duì)于遙感反演算法的波段選擇有著重要的意義，同時(shí)也能反映出研究區(qū)域內(nèi)水體的表面光學(xué)特性。本文分析2014年5月定點(diǎn)觀測(cè)遙感反射率數(shù)據(jù)共90個(gè)，如圖2所示，可以看出在590～690nm區(qū)間，大小潮遙感反射率形狀略有區(qū)別，大潮反射基臺(tái)較平，而小潮反射基臺(tái)有隨著波段增加逐漸緩慢變小的趨勢(shì)。況潤(rùn)元曾將長(zhǎng)江口遙感反射率曲線分為8種典型的形狀[9]，本航次遙感反射率曲線，可以看出從大潮到小潮，有從高泥沙濃度低葉綠素濃度向泥沙濃度降低的形狀發(fā)展的趨勢(shì)，另外在820附近都有明顯的第二峰值。長(zhǎng)江口枯季光譜特征無論在大潮還是在小潮都是典型的高泥沙濃度低葉綠素含量形式。

表1 Landsat8遙感影像6S大氣校正模型參數(shù)Table 1 Model parameters of 6S atmospheric correction of Landsat8 remote sensing

為了避免葉綠素對(duì)TSM光學(xué)特性的影響，選取了4個(gè)TSM濃度差異較大的光譜數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)TSM主要影響光譜的長(zhǎng)波波段。如圖3(a)所示，隨著泥沙濃度的增加，最大反射率波段向長(zhǎng)波方向移動(dòng)，也就是所謂的紅移現(xiàn)象，在580～690波段范圍內(nèi)存在明顯的反射峰，并且泥沙濃度越大，峰面就越平，之后隨著波段增加，遙感反射率明顯下降，在750nm附近形成一個(gè)反射谷，主要是因?yàn)樗w自身的吸收造成，另外在806nm左右波段也存在明顯的反射峰，遙感反射率迅速增加。顯然，長(zhǎng)江口地區(qū)的水體光譜特征，總懸浮物顆粒物是絕對(duì)的主導(dǎo)因素，由于其強(qiáng)烈的吸收和散射作用，其粒徑大小、顆粒物成分、圓度以及相應(yīng)的鹽度、流速、流向等水文特征和非色素顆粒物的本身特性是造成長(zhǎng)江口地區(qū)光譜特征變化的主要因素。圖3(b)所示為對(duì)光譜曲線進(jìn)行的一階微分處理，可以更客觀地反映不同泥沙濃度對(duì)光譜的響應(yīng)，其中正值表示光譜上升，負(fù)值則相反。0值的位置表示光譜曲線出現(xiàn)波峰或波谷的位置，在580～600nm附近第一次出現(xiàn)0值的位置基本一致，在806nm附近的第二反射峰位置也基本一致，第二反射峰前后光譜都有明顯的上升和下降趨勢(shì)，其變化幅度隨著TSM變大也會(huì)相應(yīng)增大。而在350～560nm之間的短波波段，遙感反射率一階微分不隨TSM變化發(fā)生明顯變化且變幅較小，說明懸浮泥沙主要影響580nm之后的中長(zhǎng)波段，這也是很多TSM反演算法選擇長(zhǎng)波波段的原因。

圖3 不同TSM濃度、相近葉綠素a濃度的遙感反射率(a)及其一階微分(b)Fig.3 The remote sensing reflectance (a) and its first order derivatives(b) with different TSM concentrations and similar chlorophyll-a concentrations

2.2 模型構(gòu)建

選擇Landsat8的五個(gè)典型可見光與近紅外波段進(jìn)行分析，即443nm、483nm、563nm、655nm、858nm。在東海綜合四種不同經(jīng)驗(yàn)和簡(jiǎn)化半經(jīng)驗(yàn)算法指數(shù)的優(yōu)勢(shì)，實(shí)現(xiàn)了CPTSM算法，在長(zhǎng)江口和杭州灣地區(qū)的高渾濁水體中也取得了良好效果[7]，借鑒其模型設(shè)計(jì)方法，本文也利用四種典型的常規(guī)經(jīng)驗(yàn)反演算法建立模型。四種典型模式分別為單波段比值算法，兩波段比值算法，可見光與近紅外波段組合算法和波段微分組合算法，表達(dá)式如下：

式(4)中C1是理論衍生值，取值不影響模型精度，采用C1=0.093的常量值，將備選的五個(gè)波段，組成所有的可能波段組合分別帶入上述四種形式進(jìn)行運(yùn)算，運(yùn)用統(tǒng)計(jì)分析方法，選擇與實(shí)測(cè)懸浮泥沙濃度相關(guān)性最好的組合建立回歸方程，相關(guān)性最好的四種組合分別為：X1=Rrs(858)/(0.093-Rrs(858))，X2=Rrs(858))/(Rrs(563)，X3=Rrs(443)+Rrs(858)，X4=Rrs(858)-Rrs(650)/(858-650)+(Rrs(650)-Rrs(563)/(650-563)+(Rrs(563)-Rrs(488))/(563-488)。

最大相關(guān)系數(shù)分別為0.928、0.894、0.964和0.832，具體分析過程在下節(jié)中討論。X1，X2，X3，X4分別表示四種不同的算法，然后將這四種算法分別取權(quán)重相加，得到復(fù)雜的經(jīng)驗(yàn)算法，如式(5)式所示：

3 討論

定點(diǎn)航次采集水樣67個(gè)經(jīng)過室內(nèi)實(shí)驗(yàn)分析，南槽大潮表層水體平均含沙量為0.1364kg/m3，小潮表層水體平均含沙量為0.0995kg/m3；北港大潮表層水體平均含沙量為0.3029 kg/m3，小潮表層水體平均含沙量為0.1056 kg/m3。表層懸浮總顆粒物濃度的變化范圍0.0365～0.5808g/l，0.1kg/m3以下TSM濃度出現(xiàn)頻率最大，多數(shù)為北潮站點(diǎn)，可能是由于該站點(diǎn)是長(zhǎng)江口徑流的主泓方向，而且位置較南槽站點(diǎn)更靠近口外，落潮時(shí)徑流攜帶大量的泥沙流經(jīng)此處，而漲潮時(shí)由于靠近口外，受到較強(qiáng)的潮流作用，泥沙含量急劇下降。葉綠素濃度總體較低，在0.4～1.4 g/m3的區(qū)間范圍內(nèi)分布較為均勻。

分析TSM在時(shí)間上的分布情況，如圖4所示可以看出TSM濃度與長(zhǎng)江口徑流與東海漲落潮周期基本一致，大潮含沙量平均值明顯高于小潮。泥沙含量最小的時(shí)候出現(xiàn)在小潮的漲急階段，由于大潮時(shí)北港站點(diǎn)天氣狀況不好，所以沒有測(cè)到完整的北港大潮數(shù)據(jù)，泥沙含量最大值5月31日北港下午13、14時(shí)，為該站點(diǎn)落急階段，春季入海懸浮泥沙的擴(kuò)散形式從冬季格局向夏季的懸沙分布格局過渡，此外風(fēng)應(yīng)力也是顯著影響渾水區(qū)的重要因素，風(fēng)浪掀沙引起河口泥沙再懸浮，可使近岸泥沙擴(kuò)散范圍增大。

利用復(fù)雜反演算法思想，分析四個(gè)指數(shù)在各波段或波段組合上與K1、K2站點(diǎn)實(shí)測(cè)TSM濃度之間的相關(guān)關(guān)系時(shí)，將OLI傳感器的五個(gè)可見光/近紅外波段和全色波段作為備選波段。分別將6個(gè)波段帶入X1，對(duì)其進(jìn)行波段敏感性測(cè)試，測(cè)試結(jié)果用X1與TSM濃度之間的相關(guān)系數(shù)表示。X1中C1為Nechad基于生物光學(xué)模型的半分析算法中的理論值，從表2可以看出C1＜0.1時(shí)會(huì)造成X1在各波段的敏感性降低，個(gè)別波段甚至?xí)霈F(xiàn)負(fù)值，而當(dāng)C1＞0.1時(shí)，不會(huì)影響X1在各波段上與TSM濃度之間的相關(guān)性，可能是TSM濃度分布的不同造成的。X1在綠色波段上敏感性較低，而在紅外波段上與TSM濃度相關(guān)性在0.9以上，選擇B6來計(jì)算X1。

表2 X1指數(shù)與TSM濃度之間的相關(guān)性分析以及C1對(duì)其相關(guān)性的影響Table 2 Correlation analysis between X1 index and TSM concentration and the effect of C1 on its correlation

分析X2與TSM濃度之間相關(guān)性，其中橫行波段為X2指數(shù)的分子，縱列波段為分母，共有30種波段組合選擇。從表3中可以看出B6波段與各波段的比值都與TSM高度相關(guān)，相關(guān)系數(shù)都接近0.9，為了降低CDOM在443nm處對(duì)TSM濃度反演的影響，X2的最佳波段選擇是B6與B2比值。

表3 X2與TSM濃度的相關(guān)性分析Table 3 Correlation Analysis of X2 and TSM Concentration

選擇了部分波段組合分析X3與TSM濃度之間的相關(guān)性，共有16種備選組合，B1-B4在第一排分別測(cè)試其單波段上的相關(guān)性，第二排為四個(gè)波段分別加上B5，最后一排為分別加上B5和B6。從表4中可以看出，敏感性最高的是后三個(gè)波段之和，可以看出在B3波段加入后的值均為各組合值之中最小的，而B1的加入使得X3與TSM濃度的相關(guān)性相對(duì)較高，主要是443nm能有效體現(xiàn)出CDOM吸收，而CDOM的物質(zhì)來源決定了其與TSM的高相關(guān)性。

表4 X3與TSM濃度的相關(guān)性分析Table 4 Correlation Analysis of X3 and TSM Concentration

X4與TSM濃度之間的相關(guān)系數(shù)如表5，大部分值都是正的，說明隨著TSM濃度的增加大部分波段上的遙感反射率都在增加，而且隨著波長(zhǎng)增加，增長(zhǎng)速度越來越快，在最后的近紅外波段增長(zhǎng)速度有所下降，在B6波段上的增長(zhǎng)速度應(yīng)為最大值，X4選用五斜率之和計(jì)算。

表5 X4與TSM濃度的相關(guān)性分析Table 5 Correlation Analysis of X4 and TSM Concentration

在2014年5月航次數(shù)據(jù)中X1的變化范圍在0.02～0.14之間，而X4變化范圍在0.0001～0.0003之間，各指數(shù)加權(quán)平均后得到復(fù)雜算法指數(shù)CP=0.19X1+0.05X2+X3+0.98X4/(0.19+0.05+1+0.98)。建立CP指數(shù)與實(shí)測(cè)TSM濃度的對(duì)數(shù)值之間的線性擬合關(guān)系，得到如(6)式的反演公式，其擬合決定系數(shù)為0.81：

通過上述公式計(jì)算出2015年3月13日長(zhǎng)江口TSM濃度分布，結(jié)果如圖4所示。

長(zhǎng)江口懸浮泥沙從徐六涇到口門再到口外，懸浮泥沙濃度不斷升高[10-12]，口外部分泥沙含量從南向北部顯示出中間高而兩邊低的形態(tài)，南部略高于北部，與河口汊道方向基本一致，最高濃度的區(qū)域出現(xiàn)在攔門沙前緣渾水峰水域。泥沙擴(kuò)散形狀呈羽毛狀和噴射狀，噴射狀主要是由于渾濁水體徑流在河槽的約束下受到落潮動(dòng)力軸慣性作用的影響，而河口沖淡水在鹽水楔之上，漂浮于海面，在風(fēng)浪作用下運(yùn)動(dòng)方向容易發(fā)生改變，在擴(kuò)散前沿形成羽毛狀分布[13]。

圖4 Landsat8遙感影像表層TSM濃度反演結(jié)果Fig.4 Landsat8 remote sensing inversion results of surface TSM concentration

4 結(jié)論

（1）長(zhǎng)江口地區(qū)表層懸浮泥沙濃度隨著潮周期的變化發(fā)生周期性變化，大潮的變差明顯大于小潮，而且平均值也較小潮大，是典型的高濁度二類水體。

（2）遙感反射率在各波段與TSM濃度的相關(guān)性不同，波段選擇與要反演的遙感影像的季節(jié)和相應(yīng)潮情有關(guān)，因此，建立適用性強(qiáng)的復(fù)雜算法十分有必要。

（3）懸浮泥沙的經(jīng)驗(yàn)反演算法對(duì)于有效利用Landsat8數(shù)據(jù)十分重要，復(fù)雜經(jīng)驗(yàn)反演算法較常規(guī)經(jīng)驗(yàn)算法更能有效反演長(zhǎng)江口地區(qū)懸浮泥沙，可為上海市懸浮泥沙濃度監(jiān)測(cè)業(yè)務(wù)化運(yùn)行提供參考。

參考文獻(xiàn)(References)

[1]Doxaran D, Froidefond J, Lavender S, et al. Spectral signature of highly turbid waters: Application with SPOT data to quantify suspended particulate matter concentrations[J].Remote Sensing of Environment, 2002,81(1):149-161.

[2]陳楚群,潘志林,施平. 海水光譜模擬及其在黃色物質(zhì)遙感反演中的應(yīng)用[J]. 熱帶海洋學(xué)報(bào),2003,22(5):33-39.Chen C Q, Pan Z L, Shi P. Simulation of sea water reflectance and its application in retrieval of yellow substance by remote sensing data[J].Journal of Tropical Oceanography, 2003,22(5):33-39.

[3]唐軍武,王曉梅,宋慶君,等. 黃 東海二類水體水色要素的統(tǒng)計(jì)反演模式[J]. 海洋科學(xué)進(jìn)展,2004,22(S1):1-7.Tang J W, Wang X M, Song Q J, et al. Statistical inversion models for case II water color elements in the Yellow Sea and East China Sea[J].Advances in Marine Science, 2004,22(S1):1-7.

[4]汪小勇,李銅基,周虹麗,等. 中國(guó)近海典型海區(qū)表觀光學(xué)特性分析與研究[J]. 光學(xué)技術(shù),2011,37(5):574-581.Wang X Y, Li T J, Zhou H L, et al. Analysis and research on apparent optical properties about typical area of Chinese offshore[J].Optical Technique, 2011,37(5):574-581.

[5]蔣雪中,王維佳,張俊儒,等. 長(zhǎng)江口最大渾濁帶表層水體懸沙粒徑對(duì)離水光譜反射率的影響[J]. 泥沙研究,2012 (5):1-7.Jiang X Z, Wang W J, Zhang J R, et al. Impact of suspended sediment particle size on water-leaving radiance reflectance in turbidity maximum zone of Changjiang estuary[J].Journal of Sediment Research, 2012(5):1-7.

[6]劉紅,何青,王亞,等. 長(zhǎng)江河口懸浮泥沙的混合過程[J]. 地理學(xué)報(bào),2012,67(9):1269-1281.Liu H, He Q, Wang Y, et al. Processes of suspended sediment mixture in the Yangtze river estuary[J].Acta Geographica Sinica,2012,67(9):1269-1281.

[7]Mao Z H, Chen J Y, Pan D L, et al. A regional remote sensing algorithm for total suspended matter in the East China Sea[J].Remote Sensing of Environment, 2012,124:819-831.

[8]Mobley C D. Estimation of the remote-sensing reflectance from above-surface measurements[J].Application Optics,1999,38(36):7442-7455.

[9]況潤(rùn)元. 長(zhǎng)江口水色遙感參數(shù)模擬研究[D]. 上海:華東師范大學(xué),2010.Kuang R Y. Remote sensing model of ocean colour parameters in Yangtze estuary[D]. Shanghai: East China Normal University, 2010.

[10]陳沈良,張國(guó)安,楊世倫,等. 長(zhǎng)江口水域懸沙濃度時(shí)空變化與泥沙再懸浮[J]. 地理學(xué)報(bào),2004,59(2):260-266.Chen S L, Zhang G A, Yang S L, et al. Temporal and spatial changes of suspended sediment concentration and resuspension in the Yangtze river estuary and its adjacent waters[J].Acta Geographica Sinica, 2004,59(2):260-266.

[11]謝火艷,王如生,張國(guó)安,等. 長(zhǎng)江口北槽沉積物的粒度特征和輸運(yùn)趨勢(shì)探討[J]. 上海國(guó)土資源,2016,37(2):84-88,95.Xie H Y, Wang R S, Zhang G A, et al. Discussion on grain-size characteristics of seafloor sediment and transport pattern in the northern passage of the Yangtze River estuary[J].Shanghai Land& Resources, 2016,37(2):84-88,95.

[12]李懿淼,李茂田,艾威,等. 長(zhǎng)江流域水庫(kù)“泥沙過濾器”效應(yīng)對(duì)溶解硅的滯留行為[J]. 上海國(guó)土資源,2017,38(2):49-53.Li Y M, Li M T, Ai W, et al. Sediment filter effect on dissolved silica retention in reservoirs of the Yangtze river basin[J].Shanghai Land & Resources, 2017,38(2):49-53.

[13]韓震,金亞秋,惲才興. 神舟三號(hào)CMODIS數(shù)據(jù)獲取長(zhǎng)江口懸浮泥沙含量的時(shí)空分布[J]. 遙感學(xué)報(bào),2006,10(3):381-386.Han Z, Jin Y Q, Yun C X. Spatial and temporal distributions of suspended sediment contents in the Yangtze river estuary using the CMODIS image data from China's SZ-3 spacecraft[J].Journal of Remote Sensing, 2006,10(3):381-386.