基于多站水力幾何法的黃河源區(qū)河流流量估算研究

白 娟，張 亦 弛，甘 甫 平，郭 藝，閆 柏 琨，楊 勝 天，李 和 謀,邢 乃 琛，馬 燕 妮，劉 琪

(1.中國自然資源航空物探遙感中心，北京 100083；2.北京師范大學(xué)地理科學(xué)學(xué)部，北京100875；3.北京師范大學(xué)水科學(xué)研究院，北京100875)

0 引言

河流流量是水文循環(huán)的關(guān)鍵參數(shù)之一，及時掌握河流流量變化對于流域水資源管理、生態(tài)修復(fù)和氣候變化研究等具有重要意義[1,2]。目前河流流量主要通過水文觀測站獲取，受復(fù)雜地形和惡劣氣候等因素限制，許多地區(qū)水文觀測站的建設(shè)和維護(hù)運(yùn)營極為困難[3]，缺少相關(guān)的水文觀測數(shù)據(jù)。利用遙感數(shù)據(jù)估算河流流量可以彌補(bǔ)地面觀測站的不足[4,5]，研究方法主要包括基于水量平衡模型[6-8]、基于河流斷面水位—流量關(guān)系曲線[9-12]、基于河流多水力特征參數(shù)[13-15]、基于河流寬度[16-21]4類反演方法。由于遙感技術(shù)對水平方向的信息捕捉能力和分辨率優(yōu)于垂直方向，因此，基于水位信息的單變量模型和多變量模型應(yīng)用相對較少，基于河流寬度反演河流流量的方法成為研究熱點[20]。

Gleason等提出的多站水力幾何(At-Many-stations Hydraulic Geometry,AMHG)法連接時空上的河流截面，本質(zhì)上擴(kuò)展了傳統(tǒng)的單站水力幾何(At-a-station Hydraulic Geometry，AHG)法，不需要任何地面觀測或先驗信息，僅根據(jù)河流寬度的時空變化就可以計算得到河流流量[22-25]。Gleason等[23]選取全球34條不同地形和氣候條件的大型河流分析AMHG方法的適用性和精度，結(jié)果表明估算流量與實測流量的相對均方根誤差(RRMSE)為26%～41%；Hagemann等[26]將該理論發(fā)展為貝葉斯水量平衡徑流反演方法(McFLI)，即BAM(Bayesian AMHG-Manning)算法，使用AMHG方法和曼寧方程以概率的方式估計河流流量；Durga Rao等[27]應(yīng)用AMHG方法估算印度4條主要河流流量，納什效率系數(shù)(NSE)均大于0.8；Mengen等[17]采用Sentinel-1A/1B數(shù)據(jù)在湄公河開展AMHG方法應(yīng)用，將時間分辨率從16 d提高至6 d，RRMSE為19.5%。Sentinel-1A/1B數(shù)據(jù)不受云和陰影限制，且具有高空間分辨率,然而，這種C波段的SAR影像在不同水體分類時存在椒鹽噪聲及反射率差異，特別是在水體與沉積物或植被相互作用區(qū)域?？紤]到AMHG方法在國內(nèi)的適用性研究較少，本研究以典型水文資料缺失地區(qū)黃河源區(qū)為研究對象，基于GEE平臺采用Landsat8和Sentinel-2影像提取河道掩膜，結(jié)合RivWidth_v04工具自動提取不同時相河流寬度，而后利用AMHG方法估算河流流量，并結(jié)合高時空分辨率的Sentinel-2影像，增加河流流量的估算頻次，盡可能反映研究河段的水文過程，對及時掌握徑流變化和開展黃河流域水資源管理及氣候變化研究均具有重要意義。

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)概況

黃河源區(qū)指黃河流域唐乃亥水文站以上區(qū)域，位于青藏高原東北部，面積約12.2萬km2，海拔2 664～6 277 m[28]，屬典型內(nèi)陸高原氣候，年均氣溫-1.6 ℃，年均降水量457 mm，5-10月降水量占全年降水量的90.4%[29]；黃河源區(qū)約占黃河流域面積的15.3%，多年平均徑流量卻占全流域的34.1%，是黃河流域重要的水源涵養(yǎng)區(qū)和產(chǎn)水區(qū)[30]，同時也是典型的水文資料缺失地區(qū)，其水文站網(wǎng)平均密度為9 380 km2/站。達(dá)日—瑪曲—軍功一帶是黃河源區(qū)的重點產(chǎn)流區(qū)，本文選取位于干流的瑪曲水文站(33.97°N，102.08°E，集水面積86 048 km2)、軍功水文站(34.7°N，100.65°E，集水面積98 414 km2)所在河段為研究對象(圖1)，兩河段流量與河寬變化明顯相關(guān)，滿足AMHG方法的適用條件。

圖1 研究區(qū)位置Fig.1 Location of the study area

1.2 數(shù)據(jù)來源

研究數(shù)據(jù)包括：1) Landsat8 OLI數(shù)據(jù)(https://search.earthdata.nasa.gov/)和Sentinel-2 數(shù)據(jù)(https://scihub.copernicus.eu/dhus/#/home)，Landsat8數(shù)據(jù)重訪周期為16 d，剔除研究河段被云遮擋的影像后，可用影像較少，故補(bǔ)充重訪周期為5 d的Sentinel-2影像(空間分辨率為10 m)，用于提取水體，兩顆衛(wèi)星互補(bǔ)重訪周期為5 d；2)GRWL(Global River Widths from Landsat)數(shù)據(jù)，用于剔除水體提取結(jié)果中的非河道水體(https://zenodo.org/record/1297434/)[31]；3)徑流實測資料，用于驗證估算流量精度，來源于水利部公布的水文站逐日平均流量數(shù)據(jù)。由于本文收集到的部分徑流實測數(shù)據(jù)僅涵蓋2015年、2017年和2021年，故模擬時段選擇2015-2021年，其中，20210422、20210711、20211001、20211113缺少實測資料，取前后兩天徑流實測均值替代，不參與模擬精度評價。

2 研究方法

研究區(qū)河流流量估算流程(圖2)為：1)河流寬度提取。采用Zou等的水體提取方法[32]，在GEE平臺上基于Landsat8/Sentinel-2影像將MNDWI>EVI或MNDWI>NDVI的像元劃分為水體，并去除EVI<0.1的水體和植被混合像元，得到研究區(qū)水體掩膜，在此基礎(chǔ)上，結(jié)合GRWL數(shù)據(jù)剔除非河流水體信息，得到河道掩膜后代入RivWidth_v04工具，得到河流寬度信息。2)河流流量估算。首先，根據(jù)實測數(shù)據(jù)和河流寬度判斷研究河段是否滿足參數(shù)a和b呈對數(shù)—線性關(guān)系，若滿足，則采用河流寬度信息建立AMHG RS Slope回歸關(guān)系，近似得到a-b對數(shù)—線性方程，而后判斷是否可以利用AMHG方法估算河流流量，若可以則代入遺傳算法(GA)估算河流流量，否則不能用AMHG方法估算河流流量。

圖2 河流流量估算流程Fig.2 Workflow of river discharge estimation

2.1 AMHG方法原理

根據(jù)單站水力幾何(AHG)方程，河流斷面寬度w與河流斷面流量Q存在指數(shù)關(guān)系(式(1))，將等式兩邊分別取對數(shù)可得式(2)，多站水力幾何(AMHG)方法證明AHG方程中參數(shù)a和b呈對數(shù)-線性關(guān)系(式(3))，可通過對同一河段內(nèi)多個固定斷面的河寬進(jìn)行多次測量以估算河道的特定常數(shù)E，斜率1/logE可由式(4)中的經(jīng)驗回歸參數(shù)y近似表示[31]。本文利用遙感影像提取的河流斷面寬度和式(3)進(jìn)行回歸分析，據(jù)此計算得到式(3)中的截距1/logE×log(wglob)，進(jìn)而確定AHG中b和loga的線性相關(guān)關(guān)系，即RS Slope回歸關(guān)系。

w=aQb

(1)

logQ=(logw-loga)/b

(2)

(3)

式中：x為河流斷面編號；wglob為所有斷面寬度的平均值，由遙感影像獲得。

max(wx1,x2,…，xn)=p(max(wx1,x2,…，xn)2-
min(wx1,x2,…，xn)2)y

(4)

式中：x1,x2,…,xn為每個河流斷面的編號；p和y為由河流斷面寬度變化決定的經(jīng)驗參數(shù)。

根據(jù)參數(shù)a、b的取值范圍以及a-b對數(shù)—線性方程，給定一對a-b的值和河流寬度則可計算得到任意斷面的流量。對此，AMHG方法采用遺傳算法(GA)估算河流流量，以不同斷面的流量差值最小化為率定目標(biāo)，同時保證估算出的流量值位于合理范圍，最后取所有流量值的均值作為該段河流的流量。

2.2 RivWidth_v04工具

RivWidth_v04工具用于測量柵格圖像中連續(xù)河流寬度[33]，由ITT可視化信息解決方案(ITTVIS)IDL語言開發(fā)而成(http://uncglobalhydrology.org/rivwidth/)。該工具的輸入文件為ENVI格式的二值水體掩膜數(shù)據(jù)，其中水體像素值為1，非水體像素值為0，主要步驟包括：提取河流掩膜、提取河流中心線和沿中線測量河流寬度。

2.3 精度評價

基于日均流量觀測數(shù)據(jù)，選取納什效率系數(shù)NSE(式(5))、均方根誤差RMSE(式(6))和相對均方根誤差RRMSE(式(7))對AMHG方法模擬結(jié)果進(jìn)行精度評價，NSE值越大，RMSE和RRMSE值越小，模擬效果越好。

(5)

(6)

(7)

3 結(jié)果分析

3.1 河流寬度提取及精度驗證

AMHG方法要求研究河段應(yīng)選在無支流匯入或流出的地區(qū)，以滿足水量守恒條件；此外，由于該方法對于辮狀河道模擬效果較差，需要在選取斷面的過程中避開河道中間季節(jié)性出露的心灘。因此，本研究在瑪曲站附近沿河道選取17個斷面，在軍功站附近選取16個斷面(圖3)，保持?jǐn)嗝婵臻g位置不變，將利用遙感影像提取的河道掩膜信息代入RivWidth_v04工具計算各斷面寬度。

圖3 瑪曲、軍功河段斷面位置Fig.3 Location of river cross-sections in Maqu and Jungong reaches

將RivWidth_v04的輸出值與使用ArcMap測量工具手動測量的河流斷面寬度進(jìn)行比較(圖4)，可以看出：對于Landsat8影像，瑪曲河段二者決定系數(shù)R2為0.87、RMSE為32.07 m，軍功河段R2為0.62、RMSE為29.22 m；對于Sentinel-2影像，瑪曲河段R2為0.97、RMSE為19.68 m，軍功河段R2為0.87、RMSE為14.23 m。對比發(fā)現(xiàn)，使用高分辨率水體掩膜數(shù)據(jù)可以降低河寬誤差。

圖4 RivWidth_v04提取河流寬度精度Fig.4 Accuracy of river width extracted by RivWidth_v04

3.2 AMHG RS Slope回歸關(guān)系求解

利用實測流量和多時相Landsat影像提取的河流斷面寬度進(jìn)行擬合得到瑪曲河段和軍功河段AMHG RS Slope回歸關(guān)系(圖5)。由圖5可以看出，瑪曲河段AMHG的回歸方程為y=-0.26x+0.64，R2為0.92，軍功河段AMHG的回歸方程為y=-0.35x+0.75，R2為0.96，表明瑪曲和軍功河流斷面寬度滿足AMHG中參數(shù)a和b的對數(shù)—線性關(guān)系，可以運(yùn)用AMHG方法計算河流流量。

圖5 瑪曲河段和軍功河段AMHG 回歸關(guān)系Fig.5 Regression relationships of AMHG for Maqu and Jungong reaches

根據(jù)多時相Landsat8影像提取的河流斷面寬度建立RS Slope回歸關(guān)系(圖6)，可以看出，瑪曲河段RS Slope的回歸方程為y=0.29x+1.09，R2為0.35，代入式(3)計算的截距為0.66，軍功河段RS Slope的回歸方程為y=0.34x+0.74，R2為0.78，代入式(3)計算的截距為0.70，兩個回歸方程斜率的絕對值和截距與AMHG回歸方程的斜率絕對值和截距均較接近，該結(jié)果進(jìn)一步表明可以利用AMHG方法估算河流流量。

圖6 瑪曲河段和軍功河段RS Slope回歸關(guān)系Fig.6 Regression relationships of RS slope for Maqu and Jungong reaches

3.3 AMHG-Landsat8/Sentinel-2估算河道流量

利用AMHG方法估算河流流量時，GA共包括8個參數(shù)：河流流量范圍(lowfilter，hifilter)、AHG方程中參數(shù)a和b的取值范圍(amin，amax，bmin，bmax)、遺傳算法參數(shù)(numGen，numGA)。本研究中參數(shù)a和b的取值范圍以及遺傳算法參數(shù)均采用默認(rèn)值，河流流量范圍根據(jù)河道附近徑流實測資料設(shè)置，具體參數(shù)設(shè)置如表1所示。

表1 研究區(qū)GA算法參數(shù)設(shè)置Table 1 Parameterization of the genetic algorithm for the study area

2015-2021年瑪曲、軍功河段AMHG-Landsat8估算流量與實測流量如圖7所示。由于GA算法缺乏穩(wěn)定性，本研究在兩個河段分別進(jìn)行10次模擬，圖7中紅色區(qū)間表示10次模擬的估算結(jié)果取值范圍，本文對10次模擬結(jié)果取平均值作為最終估算結(jié)果與實測數(shù)據(jù)進(jìn)行對比。結(jié)果表明，瑪曲河段河流流量估算值與實測值的NSE為0.75、RMSE為155.66 m3/s、RRMSE為34.04%，軍功河段河流流量估算值與實測值的NSE值為0.67、RMSE為225.15 m3/s、RRMSE為38.86%，總體上，兩個河段的模擬結(jié)果均較好。為分析河流流量范圍參數(shù)取值對估算結(jié)果的影響，采用表1中l(wèi)owfilter，hifilter推薦取值(瑪曲河段取值分別為7.5和13 000，軍功河段取值分別為4.6和7 754)進(jìn)行模擬，其余參數(shù)保持不變。由圖7可知，河流流量范圍取值對模擬結(jié)果有很大影響，采用推薦取值的河流流量范圍在兩個河段的模擬結(jié)果均存在嚴(yán)重低估，與Gleason等[34]的結(jié)論一致，因此，利用AMHG方法估算河流流量時應(yīng)考慮河流的實際最小和最大流量限制。

圖7 2015-2021年瑪曲和軍功河段AMHG估算流量和實測流量Fig.7 Estimated discharge based on AMHG and measured discharge of Maqu and Jungong reaches from 2015 to 2021

采用相同的GA算法參數(shù)設(shè)置，對比兩個河段AMHG-Landsat8、AMHG-Sentinel-2估算流量的模擬效果。2021年瑪曲河段AMHG-Landsat8估算流量頻次為7次，與實測流量的NSE為0.86、RMSE為162.34 m3/s、RRMSE為31.06%，AMHG-Sentinel-2估算流量頻次為34次，與實測流量的NSE為0.81、RMSE為126.59 m3/s、RRMSE為21.96%(圖8a)，模擬結(jié)果均較好。2021年軍功河段AMHG-Landsat8估算流量頻次為6次，與實測流量的NSE為0.74、RMSE為250.1 m3/s、RRMSE為33.48%，AMHG-Sentinel-2估算流量頻次為14次，與實測流量的NSE為0.68、RMSE為145.9 m3/s、RRMSE為22.02%(圖8b)，表明采用Sentinel-2影像可顯著增加河流流量估算頻次和精度。對比圖8中瑪曲河段和軍功河段AMHG方法估算流量模擬效果，發(fā)現(xiàn)非汛期模擬效果較好，汛期估算流量低于實測流量，該結(jié)果與Durga Rao等[27]在印度4條主要河流應(yīng)用AMHG方法的研究結(jié)論相同。分析原因，可能是由于當(dāng)汛期河流水量超過河岸時，AHG冪律行為遭到破壞[24]，也可能是由于基于衛(wèi)星影像的AMHG方法估算的河流流量反映的是衛(wèi)星影像過境時的瞬時流量，與日徑流實測值存在一定偏差。

圖8 2021年瑪曲和軍功河段AMHG估算流量和實測流量Fig.8 Estimated discharge based on AMHG and measured discharge of Maqu and Jungong reaches in 2021

4 結(jié)論

本文通過自動提取不同時相Landsat8影像的河流斷面系列寬度信息，分析了AMHG-Landsat8在黃河源區(qū)的適用性，在此基礎(chǔ)上，結(jié)合高時空分辨率Sentinel-2影像增加河流流量的估算頻次，以反映研究河段的水文過程。結(jié)論如下：1)應(yīng)用RivWidth_v04工具提取河流斷面寬度精度較高，使用高分辨率遙感影像可降低河寬誤差；2)2015-2021年瑪曲、軍功河段AMHG-Landsat8模擬值與實測值的NSE分別為0.75和0.67，RMSE分別為155.66 m3/s和225.15 m3/s，RRMSE分別為34.04%和38.86%，模擬結(jié)果均較好；3)結(jié)合高時空分辨率Sentinel-2影像可顯著增加河流流量的估算頻次，更好地反映研究河段的水文過程；4)河流流量范圍取值對AMHG模擬結(jié)果有很大影響，應(yīng)用AMHG方法時需考慮河流實際的最小和最大流量限制。

本文采用Sentinel-2影像一定程度上提高了AMHG方法在河流監(jiān)測和水資源管理領(lǐng)域的實用性，但影像易受天氣影響，導(dǎo)致觀測頻次不足。下一步擬嘗試聯(lián)合Sentinel-1 SAR數(shù)據(jù)，實現(xiàn)不受云和天氣影響的連續(xù)河流寬度測量和長時間序列的河流流量觀測，滿足水資源管理要求，并將其應(yīng)用于其他河流系統(tǒng)。