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

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

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

0 引言

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

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

1 研究區域與數據

1.1 研究區概況

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

圖1 研究區位置Fig.1 Location of the study area

1.2 數據來源

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

2 研究方法

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

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

2.1 AMHG方法原理

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

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為由河流斷面寬度變化決定的經驗參數。

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

2.2 RivWidth_v04工具

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

2.3 精度評價

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

(5)

(6)

(7)

3 結果分析

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

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

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

將RivWidth_v04的輸出值與使用ArcMap測量工具手動測量的河流斷面寬度進行比較(圖4)，可以看出：對于Landsat8影像，瑪曲河段二者決定系數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。對比發現，使用高分辨率水體掩膜數據可以降低河寬誤差。

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

3.2 AMHG RS Slope回歸關系求解

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

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

根據多時相Landsat8影像提取的河流斷面寬度建立RS Slope回歸關系(圖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回歸方程的斜率絕對值和截距均較接近，該結果進一步表明可以利用AMHG方法估算河流流量。

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

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

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

表1 研究區GA算法參數設置Table 1 Parameterization of the genetic algorithm for the study area

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

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

采用相同的GA算法參數設置，對比兩個河段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)，模擬結果均較好。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方法估算流量模擬效果，發現非汛期模擬效果較好，汛期估算流量低于實測流量，該結果與Durga Rao等[27]在印度4條主要河流應用AMHG方法的研究結論相同。分析原因，可能是由于當汛期河流水量超過河岸時，AHG冪律行為遭到破壞[24]，也可能是由于基于衛星影像的AMHG方法估算的河流流量反映的是衛星影像過境時的瞬時流量，與日徑流實測值存在一定偏差。

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

4 結論

本文通過自動提取不同時相Landsat8影像的河流斷面系列寬度信息，分析了AMHG-Landsat8在黃河源區的適用性，在此基礎上，結合高時空分辨率Sentinel-2影像增加河流流量的估算頻次，以反映研究河段的水文過程。結論如下：1)應用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%，模擬結果均較好；3)結合高時空分辨率Sentinel-2影像可顯著增加河流流量的估算頻次，更好地反映研究河段的水文過程；4)河流流量范圍取值對AMHG模擬結果有很大影響，應用AMHG方法時需考慮河流實際的最小和最大流量限制。

本文采用Sentinel-2影像一定程度上提高了AMHG方法在河流監測和水資源管理領域的實用性，但影像易受天氣影響，導致觀測頻次不足。下一步擬嘗試聯合Sentinel-1 SAR數據，實現不受云和天氣影響的連續河流寬度測量和長時間序列的河流流量觀測，滿足水資源管理要求，并將其應用于其他河流系統。