基于遙感技術(shù)對(duì)松花江懸浮泥沙濃度的分析

韓 旭，杜 崇，劉福全，李 瑞

(黑龍江大學(xué) 水利電力學(xué)院，哈爾濱 150080)

0 引 言

河流的懸浮泥沙是河流重要的水質(zhì)參數(shù)之一，泥沙輸移是河流中重要的水文現(xiàn)象。水體懸浮泥沙濃度影響一系列的水體光學(xué)特性，如河流水色、水體渾濁度、水體透明度等，進(jìn)而影響水生生物的生態(tài)環(huán)境，例如光合作用、初級(jí)生產(chǎn)力、營(yíng)養(yǎng)流動(dòng)、河流生物多樣性等。在天然江河修建水壩，水壩上游水位抬高，流速變緩，河流中的懸浮泥沙在壩底淤積，泥沙壓力可能造成水壩應(yīng)力分布不均，破壞壩體結(jié)構(gòu)和穩(wěn)定性。因此在泥沙含量偏高的河流中修建水壩，應(yīng)考慮到應(yīng)力的影響。此外，河流泥沙含量對(duì)河道平衡也有很大影響，局部泥沙淤積造成河床變形等問(wèn)題。因此對(duì)河流泥沙的動(dòng)態(tài)監(jiān)測(cè)有重要意義[1]。傳統(tǒng)測(cè)量懸浮泥沙的方法是對(duì)采樣區(qū)內(nèi)的采樣點(diǎn)進(jìn)行逐個(gè)點(diǎn)的長(zhǎng)時(shí)間監(jiān)測(cè)，這種方法不僅耗時(shí)過(guò)長(zhǎng)，而且存在著成本高和無(wú)法對(duì)大面積水域進(jìn)行長(zhǎng)時(shí)間監(jiān)測(cè)等多方面缺點(diǎn)。遙感技術(shù)的監(jiān)測(cè)范圍廣、耗時(shí)短，且遙感圖像獲取的成本較低，運(yùn)用遙感技術(shù)對(duì)河流懸浮泥沙進(jìn)行監(jiān)測(cè)成為新的發(fā)展趨勢(shì)。受制于遙感圖像分辨率等多方面因素，遙感技術(shù)對(duì)懸浮泥沙的反演多是對(duì)海洋和湖泊等面狀水體，對(duì)河流線狀水體的研究較少。韓震等[2]對(duì)遙感衛(wèi)星單波段的反射率和實(shí)測(cè)懸浮泥沙濃度構(gòu)建了數(shù)學(xué)模型對(duì)其進(jìn)行研究，但單波段的反演模型無(wú)法排除水中葉綠素等物質(zhì)的干擾。方馨蕊等[3]使用隨機(jī)森林回歸模型對(duì)懸浮泥沙濃度進(jìn)行遙感估算，但機(jī)器模型需要大量的訓(xùn)練樣本，很難在河流遙感的定量反演中應(yīng)用。林承達(dá)等[4]利用LandsatETM+遙感影像對(duì)長(zhǎng)江中游懸浮泥沙濃度進(jìn)行遙感反演，將遙感技術(shù)運(yùn)用到線狀水體懸浮泥沙的監(jiān)測(cè)上，但LandsatETM+影像分辨率為30 m，對(duì)于線狀水體，其分辨率較低。本文使用Sentinel-2衛(wèi)星數(shù)據(jù)，其遙感影像分辨率經(jīng)過(guò)SNAP軟件重采樣后分辨率為10 m，且相較于Landsat系列陸地衛(wèi)星，Sentinel-2衛(wèi)星重訪周期更短[5]。

1 數(shù)據(jù)與方法

1.1 研究區(qū)域概況

由于大頂子山航電樞紐工程，大壩上游松花江哈爾濱段水位已蓄至115 m，兩岸大面積的平地和灘地被回水覆蓋，形成從樞紐大壩至肇源縣長(zhǎng)達(dá)128 km的人工湖。大頂子山航電樞紐工程蓄水后，松花江左岸原有的濕地資源被成倍放大。位于呼蘭區(qū)的“黑龍江呼蘭河濕地自然保護(hù)區(qū)”申請(qǐng)立項(xiàng)，經(jīng)批準(zhǔn)的保護(hù)區(qū)總面積為19 262 hm2，其中林地1 680 hm2，蘆葦濕地、濕生草地、草甸等近萬(wàn)hm2，是國(guó)內(nèi)目前最大的城區(qū)濕地。大頂子山航電樞紐建成后，松花江哈爾濱段的通航期由原來(lái)的100 d左右提高至210 d左右，在枯水期時(shí)，1 000 t船隊(duì)也可以在松花江哈爾濱段暢通無(wú)阻。對(duì)于河流遙感的定量反演，河面寬度增大，衛(wèi)星傳感器足夠清晰識(shí)別河面目標(biāo)，適合作為本次實(shí)地采樣區(qū)域。

近年來(lái)，松花江大力發(fā)展航運(yùn)，大頂子山樞紐附近的水運(yùn)量在1×107t以上，研究該區(qū)域的泥沙分布，不僅為航道建設(shè)、航路規(guī)劃提供數(shù)據(jù)支持，也為防止大頂子山航電樞紐泥沙淤積提供參考。

1.2 實(shí)測(cè)數(shù)據(jù)

枯水期和豐水期野外實(shí)地采樣時(shí)間分別為5月12日，7月13日。兩次采集的采樣點(diǎn)均是14個(gè)，采樣點(diǎn)選取時(shí)，選擇江面平緩的位置，并且當(dāng)日無(wú)降水，反演模型建立點(diǎn)均勻分布在河流橫斷面上，確保反演模型的準(zhǔn)確性。14個(gè)采樣點(diǎn)中，1～8號(hào)點(diǎn)用于反演模型的建立，9～14號(hào)點(diǎn)用于模型精度的驗(yàn)證。使用手持GPS進(jìn)行定位，確保兩次采集點(diǎn)位置相同，懸浮泥沙采樣方法與實(shí)驗(yàn)室測(cè)量方法參照國(guó)家河流懸移質(zhì)泥沙監(jiān)測(cè)規(guī)范(GB/T 50159-2015)。采樣點(diǎn)位置見(jiàn)圖1。

圖1 采樣點(diǎn)位置

實(shí)地測(cè)量前先測(cè)量每個(gè)采樣點(diǎn)對(duì)應(yīng)容器和濾膜的質(zhì)量，采用容器為500 mL寬口塑料瓶，因?yàn)楦鱾€(gè)波段對(duì)水體的穿透性有限，現(xiàn)場(chǎng)采樣水深為0.5 m，每個(gè)采樣點(diǎn)采集3瓶水樣取平均值，目的是為了減少誤差和后期數(shù)據(jù)篩選。采集水樣后，即到實(shí)驗(yàn)室進(jìn)行處理?？菟诤拓S水期各采樣點(diǎn)懸浮泥沙濃度見(jiàn)表1和表2。

表1 枯水期實(shí)測(cè)采樣點(diǎn)懸浮泥沙濃度

表2 豐水期實(shí)測(cè)采樣點(diǎn)懸浮泥沙濃度

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

本文使用的遙感數(shù)據(jù)均來(lái)自Sentinel-2衛(wèi)星，該衛(wèi)星遙感數(shù)據(jù)在歐空局網(wǎng)站(https:∥scihub.copernicus.eu)下載，Sentinel-2衛(wèi)星于2015年6月23日升空，遙感圖像寬幅為290 km，重訪周期為5 d，擁有13個(gè)波段，各波段參數(shù)見(jiàn)表3。

表3 Sentinel-2衛(wèi)星各波段參數(shù)

實(shí)測(cè)時(shí)間為5月12日下午14：00左右，7月13日上午10：00左右。由于Sentinel-2衛(wèi)星圖像標(biāo)注時(shí)間為格林尼治時(shí)間，格林尼治時(shí)間與北京時(shí)間相差8 h，所以5月12日的Sentinel-2衛(wèi)星圖像時(shí)間選擇盡可能接近上午6：00，7月13日的Sentinel-2衛(wèi)星圖像時(shí)間選擇盡可能接近凌晨2：00，減小由時(shí)間差異造成的誤差。最終選擇5月12日上午6：32的圖像和7月13日上午4：17的圖像，該圖像云層覆蓋率分別為3.12%和9.22%，均小于10%，能見(jiàn)度高，適合運(yùn)用于反演模型的建立。

因Sentinel-2衛(wèi)星經(jīng)過(guò)處理的L2A數(shù)據(jù)沒(méi)有與本次實(shí)際采樣時(shí)間對(duì)應(yīng)的遙感影像，本次Sentinel-2衛(wèi)星數(shù)據(jù)下載的是L1C級(jí)產(chǎn)品，經(jīng)過(guò)幾何校正，沒(méi)有進(jìn)行輻射定標(biāo)和大氣校正，需要進(jìn)行遙感數(shù)據(jù)的預(yù)處理。大氣校正方法包括ENVI中FLASSH模塊大氣校正、6S模型大氣校正和Sen2cor等多種方法。由于本文選用Sentinel-2衛(wèi)星遙感圖像，選用由歐空局提供的Sen2cor將L1C級(jí)數(shù)據(jù)進(jìn)行處理，得到L2A級(jí)大氣底層反射率數(shù)據(jù)[6]。

Sentinel-2衛(wèi)星各個(gè)波段的空間分辨率不同，需要進(jìn)行重采樣將空間分辨率統(tǒng)一，運(yùn)用SNAP軟件將各波段的空間分辨率重采樣為10 m，重采樣方法為雙線性插值法。輸出數(shù)據(jù)為ENVI可以打開(kāi)的格式，運(yùn)用ENVI中的LayerStacking工具進(jìn)行波段融合后，進(jìn)行圖像裁剪和拼接處理，獲得各個(gè)波段實(shí)測(cè)點(diǎn)的遙感反射率數(shù)據(jù)。

2 反演模型構(gòu)建

2.1 敏感波段選取

Senlinel-2A擁有13個(gè)波段，第10波段用于大氣校正，不參與相關(guān)性研究，從第1波斷開(kāi)始，中心波長(zhǎng)不斷增加，第1波段到第12波段其波長(zhǎng)為0.433～2.190 μm，由于松花江相較于國(guó)內(nèi)其他水系是少沙河流，在較清澈水體中，河流的反射率和波長(zhǎng)大致呈反比關(guān)系，即波長(zhǎng)越大，反射率越低。通過(guò)影像分析和波譜分析，在波長(zhǎng)1.0～2.190 μm，泥沙較少的自然清澈水體反射率極低，近乎為0，可以忽略不計(jì)。此外，Band5、Band6、Band7是植被紅邊波段，對(duì)水體的敏感性較低，最終選擇Band2藍(lán)光波段、Band3綠光波段、Band4紅光波段、Band8近紅外波段作為研究波段[7]。

2.2 波段相關(guān)性分析

目前用于懸浮泥沙反演的波段組合有單波段、波段比值、多波段組合等多種方法。波段比值和多波段組合，可以減輕水體葉綠素和其他懸浮物質(zhì)對(duì)反演結(jié)果的影響。本文利用ENVI中的Band math工具計(jì)算單波段、波段比值、3波段組合、4波段組合數(shù)據(jù)，并與懸浮泥沙濃度做Pearson相關(guān)性分析[8]。各波段與懸浮泥沙的相關(guān)性系數(shù)見(jiàn)表4。

由表4可見(jiàn)，單波段的相關(guān)性都比較低，尤其Band8近紅外波段的相關(guān)性在枯水期和豐水期的相關(guān)性系數(shù)只有0.096和0.041，相關(guān)性極低。只有Band3綠光波段和Band4紅光波段的P值小于0.5呈正相關(guān)性，并且其相關(guān)性系數(shù)不大。在波段比值分析中，總體相關(guān)性高于單波段相關(guān)性，枯水期的B2/B3、B2/B4、B3/B4的相關(guān)性大于0.5，豐水期的B2/B3、B2/B4的相關(guān)性大于0.5，表現(xiàn)良好。在多波段的相關(guān)性分析中，相較于波段比值和單波段其相關(guān)性稍好，枯水期的B2/(B3+B4)、B3/(B2+B4)相關(guān)性大于0.6，豐水期的B3/(B2+B4)、(B2+B3)/B4相關(guān)性大于0.6，可以滿足反演模型建立的相關(guān)性需求。總體來(lái)說(shuō)，無(wú)論單一Band8近紅外波段，還是波段比值和多波段組合中涉及Band8近紅外波段，其與懸浮泥沙的相關(guān)性都比較低。在對(duì)黃河懸浮泥沙含量進(jìn)行監(jiān)測(cè)時(shí)，近紅外波段是一個(gè)不可或缺的波段，但對(duì)松花江懸浮泥沙監(jiān)測(cè)時(shí)相關(guān)性極低。分析原因，可能是由于松花江是少沙河流，Band8近紅外波段對(duì)其反應(yīng)不敏感，或是采樣區(qū)域葉綠素等有機(jī)質(zhì)含量過(guò)多，干擾了Band8近紅外波段對(duì)懸浮泥沙的監(jiān)測(cè)效果[9]。

表4 各波段與懸浮泥沙相關(guān)性系數(shù)

2.3 反演模型建立

根據(jù)波段相關(guān)性分析，本文選取相關(guān)性系數(shù)大于0.6的波段組合建立反演模型，反演模型選用線性模型、指數(shù)模型、二次多項(xiàng)式模型、對(duì)數(shù)模型。運(yùn)用SPSS曲線估算模塊進(jìn)行反演模型的建立，遙感反射率作為自變量，懸浮泥沙濃度作為因變量[10]。選出擬合度最高的模型作為最終的反演模型?？菟诤拓S水期具體模型關(guān)系式和擬合度分別見(jiàn)表5和表6。

表5的擬合度分析結(jié)果表明B3/(B2+B4)的二次模型擬合度最高，擬合度為0.904，表6的擬合度分析結(jié)果表明B3/(B2+B4)的對(duì)數(shù)模型擬合度最高，擬合度為0.832，基于此結(jié)果，提出松花江枯水期和豐水期懸浮泥沙濃度的反演模型分別為

表5 懸浮泥沙濃度經(jīng)驗(yàn)?zāi)Ｐ?枯水期)

表6 懸浮泥沙濃度經(jīng)驗(yàn)?zāi)Ｐ?豐水期)

f(x)=0.589x2-0.291x+0.124

(1)

f(x)=1.118lnx+0.112

(2)

2.4 模型精度驗(yàn)證

對(duì)松花江懸浮泥沙濃度反演模型的精度評(píng)價(jià)，采用平均相對(duì)誤差(MAPE)和均方根誤差(RMSE)兩種方法。

1)平均相對(duì)誤差(MAPE)：求出每個(gè)點(diǎn)相對(duì)誤差的絕對(duì)值后取平均值，平均相對(duì)誤差越小證明模型精度越高，反之平均相對(duì)誤差越大則模型精度越低。

(3)

2)均方根誤差(RMSE)：又稱標(biāo)準(zhǔn)誤差，是模擬值與實(shí)測(cè)值之差的平方與測(cè)量次數(shù)n的比值的平方根，均方根誤差越小表明模型精度越高，反之均方根誤差越大表明模型精度越低[11]。

(4)

式中：n為樣本數(shù)量；yi為實(shí)測(cè)懸浮泥沙含量；fi為模擬懸浮泥沙含量。

本次模型驗(yàn)證運(yùn)用實(shí)測(cè)采集14個(gè)點(diǎn)中的6個(gè)驗(yàn)證點(diǎn)，在枯水期，6個(gè)驗(yàn)證點(diǎn)中最大相對(duì)誤差為27.17%，最小相對(duì)誤差為12.3%，平均相對(duì)誤差為17.42%，均方根誤差為0.034 g·L-1。在豐水期，6個(gè)驗(yàn)證點(diǎn)中最大相對(duì)誤差為26.08%，最小相對(duì)誤差為9.8%，平均相對(duì)誤差為16.02%，均方根誤差為0.016 g·L-1。模型精度基本滿足要求。實(shí)測(cè)值與模擬值的對(duì)比見(jiàn)表7。

由表7可見(jiàn)，模擬值總體高于實(shí)測(cè)值，除個(gè)別點(diǎn)外，模擬值的大小與實(shí)測(cè)值的大小基本呈擬合狀態(tài)。造成誤差的原因：①遙感圖像時(shí)間與實(shí)測(cè)時(shí)間有近8 h的時(shí)差，河流流速等水文條件有所變化[12]；②遙感圖像有少量云層覆蓋，進(jìn)而影響后續(xù)大氣校正結(jié)果；③河流中含有葉綠素等物質(zhì)，影響敏感波段遙感反射率；④實(shí)驗(yàn)過(guò)程中的過(guò)濾，烘干等操作也會(huì)造成誤差[13]。

表7 驗(yàn)證點(diǎn)模擬值與實(shí)測(cè)值

2.5 反演效果圖

得出枯水期和豐水期的反演模型并進(jìn)行模型精度驗(yàn)證后，運(yùn)用ENVI軟件制作出枯水期和豐水期懸浮泥沙濃度的反演結(jié)果見(jiàn)圖2。

圖2 反演效果

3 結(jié) 論

豐水期的懸浮泥沙濃度總體大于枯水期的懸浮泥沙濃度，其主要原因是豐水期水位上升，將岸邊大量的泥沙帶入河中，致使河流中的泥沙濃度上升。

懸浮泥沙濃度在河流中的分布總體呈河中心低，岸邊高的趨勢(shì)，其主要原因是岸邊流速較低，泥沙容易淤積，且影響遙感反射率的其他雜質(zhì)較多，反演效果圖呈現(xiàn)岸邊的懸浮泥沙濃度偏高。