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

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

0 引 言

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

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

1.1 研究區(qū)域概況

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

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

1.2 實測數(shù)據(jù)

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

圖1 采樣點位置

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

表1 枯水期實測采樣點懸浮泥沙濃度

表2 豐水期實測采樣點懸浮泥沙濃度

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

本文使用的遙感數(shù)據(jù)均來自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個波段，各波段參數(shù)見表3。

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

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

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

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

2 反演模型構建

2.1 敏感波段選取

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

2.2 波段相關性分析

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

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

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

2.3 反演模型建立

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

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

表5 懸浮泥沙濃度經(jīng)驗模型(枯水期)

表6 懸浮泥沙濃度經(jīng)驗模型(豐水期)

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

(1)

f(x)=1.118lnx+0.112

(2)

2.4 模型精度驗證

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

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

(3)

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

(4)

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

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

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

表7 驗證點模擬值與實測值

2.5 反演效果圖

得出枯水期和豐水期的反演模型并進行模型精度驗證后，運用ENVI軟件制作出枯水期和豐水期懸浮泥沙濃度的反演結果見圖2。

圖2 反演效果

3 結 論

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

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