基于DEM與河網密度的水系提取與應用
——以瀾湄流域為例

曹 爽，王 蒙，包紅軍，張雨鳳

(1.國家氣象中心，北京 100081；2.中國氣象局-河海大學水文氣象研究聯合實驗室，北京 100081；3.江蘇省水文水資源勘測局蘇州分局，江蘇 蘇州 215011)

0 引 言

流域水系和流域邊界能有效反映流域地形地貌特征，是流域面雨量監測預報和分布式水文模型構建的重要基礎數據，是數字孿生流域建設的主要依據[1]。它們也是完成面雨量實況監測、預報預測及預警服務的關鍵環節和前提條件，數據精度對相關計算和模擬結果有著直接影響[2]。隨著地理信息系統(GIS)和數字高程模型(DEM)的迅猛發展，國內外不斷有學者對提高提取流域河網精度和子流域劃分進行深入研究和改進，以期得到最接近真實河網的數字化水系。一些學者[3-6]以DEM數據源為基礎通過優化改進算法提出了移動窗口法、坡面徑流模擬法、谷線搜索法、從DEM直接提取河網與劃分子流域的方法等；還有一些學者[7-12]引入矢量河網作為DEM控制條件，通過Agree算法、Burn-in算法、Stream Burning算法等對主要河流進行高程與流向修正，從而解決平原區和受人類活動干擾的城市化的水系提取問題；另外，設置合理的網格數和閾值范圍也是控制河網精度的決定性因素之一。有大量研究表明[13-19]，可以通過河網密度法、水系分形維數法、均值變點分析法、河道平均坡降法、流域寬度分布法等方法確定最佳集水面積閾值；在子流域劃分方面，有些學者提出了利用河流分汊和拓撲關系[20]、多閾值虛擬河網融合技術[21]等劃分子流域，還提出了考慮高山流域[22]、湖泊水庫范圍[23]、山區平原地貌差異[24]、反映地表下墊面類型特征[25]的子流域劃分和編碼方法。

在前人研究成果的基礎上，本文以瀾滄江-湄公河流域(以下簡稱“瀾湄流域”)為研究對象，從遙感影像獲取矢量河網數據，采用Agree算法、應用Arcgis軟件疊加矢量河網從而修正主要河道DEM高程數據；分析河網密度與集水面積閾值的變化關系，以曲線割線斜率法定量確定最佳集水面積閾值，提取流域水系、邊界并劃分子流域；基于河網“套合差”法進行精度評價，驗證提取河網的可靠性；再以一次降雨過程為例分析流域面雨量的空間分布，將結果應用于流域面雨量監測、預報服務業務，為水文氣象預報預警工作提供技術支撐。

1 研究區概況與數據

本文以瀾湄流域作為研究區，流域范圍介于9°26′～34°2′ N，93°46′～108°52′ E，河長4 900 km左右，流域面積約81萬km2，是東南亞地區最重要也是最大的國際河流。該流域發源于中國青海省唐古拉山脈，流經中國、緬甸、老撾、泰國、柬埔寨、越南等六國，以云南省南臘河口為界，中國境內河段稱為瀾滄江，境外河段稱為湄公河，終由越南胡志明市流入南海。瀾湄流域地形地勢呈北高南低態勢，高度落差大，表現出很強的垂直地帶性，高程隨緯度減小逐漸下降，地貌由北向南依次是高山峽谷區、中低山寬谷區和沖積平原，海拔跨度從3 500～5 000 m降至1 000～3 000 m，再降至1 000 m以下。受地形地貌影響，瀾滄江-湄公河表現為典型的南北向狹長型河流，流域形狀似“帚狀”；上游支流短小且少，呈 “樹枝狀”、“羽狀”水系；中游東岸支流水系發育較好且西岸少有大支流，呈“梳狀”水系；下游多為平原和三角洲，河網特別發育，呈“辮狀”、“格狀”水系特征[26-27]，見圖1。此外，在氣候變化和人類活動的共同作用下，瀾湄流域頻發旱澇災害[28]，坡地山區易發山洪常伴有泥石流，河谷壩區、平原和三角洲的耕地城鎮被洪水侵襲，瀾滄江中游常有“焚風”效應導致干旱河谷等災害事件。

圖1 瀾湄流域概況

本文采用NASA SRTM的空間分辨率為90 m的DEM數據，以及國家基礎地理信息中心提供的全球1∶1 000 000包含行政區劃、湖泊、水系等信息的矢量數據，而降水數據采用中央氣象臺5 km×5 km分辨率的智能網格降水。

2 研究方法

2.1 河網提取方法

2.1.1 河網修正Agree算法

DEM數據是高程信息的反演產品，其本身就存在一定誤差。在提取水系研究河網分布時無法綜合考慮地形地貌、雨水沖刷、人工河道等其他因素的共同作用，會直接影響河網提取精度。尤其是在地勢平坦、受人類活動影響大的地區表現明顯[9,11]。

瀾湄流域下游地區即為地勢落差小的沖積平原，且耕地、人工開鑿水渠、城市化等人類活動影響大。以DEM作為單一要素提取的河網水系與實際分布會存在較大差異，遂引入Agree算法修正DEM。該算法原理是通過疊加矢量河網，降低實際河網所在柵格的高程，控制水流方向，增加提取河網柵格的“匯流”能力。

算法主要步驟包括：①通過衛星遙感影像解譯、Google earth清繪河道中線或國家基礎地理信息中心提供的地理信息數據等途徑獲取所需的實際主要河網數字化水系的矢量文件；②疊加矢量河網和填洼后的DEM，以矢量線要素為中心做緩沖區，緩沖區鄰域半徑的設置應不小于1/2個DEM分辨率；③運用轉換工具將矢量緩沖區范圍轉為柵格文件，進而運算降低重疊部分DEM柵格高程，獲取修正后DEM，命名為DEM-A。

2.1.2 流域水系提取

基于DEM的水系和流域邊界提取技術路線見圖2。工作原理是基于地表徑流漫流模型，關鍵步驟包括：迭代計算洼地填平；D8單流向算法確定水流流向；沿水流流向計算每個柵格單元的上游匯流能力即匯流累積量；設置匯流閾值，匯流量不小于閾值的柵格是潛在河網；根據Strahler水系分級法對河流進行分級；捕捉傾斜點和計算分水嶺，獲取自然全流域和自然子流域出口和流域邊界；水系河網和流域邊界矢量化[14]。

圖2 水系和流域邊界提取技術路線

2.2 確定最佳集水閾值方法

在流域提取過程中有一個重要的變量參數是集水面積閾值。閾值的變化控制著水系河網的疏密程度，目視判讀的方式存在一定的主觀性，計算結果難以達成統一。于是，本文采用河網密度法，運用擬合指數函數曲線的割線所對應的斜率與擬合曲線相切的數學方法來確定切點，該切點的物理意義即為河網密度變化趨于平穩時的最佳集水面積閾值。

擬合集水面積閾值與河網密度關系曲線，構建集水面積閾值與河網密度關系方程f(x)，定義域為[x1，x2]，值域為[f(x2)，f(x1)]，令f(x)的一階導數為f′(x)。在定義域內，曲線割線的斜率k=[f(x1)-f(x2)]/(x1-x2)，以該斜率作為擬合曲線的切線，即k=f′(x0)，求解x0，切線與擬合曲線的切點(x0，y0)代表了河網密度隨集水面積閾值由劇烈變化變為平緩變化的轉折點。即，x0為目標閾值[16]。

2.3 河網精度評價

本文對河網精度評價分為兩個方面：一方面，檢驗提取的水系和數字化水系主要河道河長的相對誤差，有效反映偏離真值的實際大??；另一方面，引入“河網套合差”的概念，指的是疊加提取水系和數字化水系兩部分，兩水系由于位置偏移會形成細碎的多邊形，計算這些細碎多邊形面積占流域總面積的比值即為套合差，該比值越小代表兩條水系的位置偏差越小[14,29](見圖3)。

圖3 河網套合差示意

2.4 流域面雨量計算

流域面雨量是指某一流域或區域整個面上的平均降雨量。它能客觀地反映流域的降水情況，在水文模型和水情預報等工作中應用廣泛。本文的降水資料數據采用的是中央氣象臺5 km×5 km網格的智能網格預報降水產品，以算術平均法估算流域面雨量預報結果適用可行，計算簡單。具體計算原理及公式不再詳述[30]。

3 結果與分析

3.1 集水面積閾值與水系河網的變化關系

本文選擇1 000、10 000、15 000、20 000、25 000、30 000、40 000、50 000、60 000、100 000共10個河網柵格數為閾值樣本，在提取水系河網的過程中以此判斷匯流累積量柵格數量的可變參數影響河網疏密程度。不同集水面積閾值對水系提取的影響特征見表1。

表1 受集水面積閾值影響提取水系河網特征變化

由表1可看出，隨著閾值面積由8.1 km2逐漸增加至810 km2的變化過程中，河道總數從44 429條減少為433條，河源數從22 261個減少為217個，總河長從198 854.57 km減至23 363.51 km，河道級別從八級河流減少至五級河流；在集水面積閾值擴大了100倍的情況下，河道總數、河源數相應減少為1/100，總河長則減少成了近原來的1/9，流域面積雖受到一定影響有一些減少但變化不大。另外，流域河網分級結果也存在一定變化規律。在同一集水面積閾值下，河道數和河長隨河網級別增加逐漸減少且減小幅度逐漸變緩，一級河道數占河道總數50%左右，二級河道數和河長是一級河道數和河長的1/2，一級和二級的河長和河道數在總數的占比達到75%左右。

3.2 流域最佳閾值的確定與精度評價

依據表1中河網閾值柵格數對應的河網密度繪制圖4關系曲線。由圖4可知，隨閾值增大，河網密度逐漸減小，閾值范圍為103～104，最大減小幅度呈斷崖式；閾值范圍為104～4×104，仍呈減小趨勢但下降幅度逐漸變緩；閾值范圍為4×104～105，變化逐漸趨于穩定。

圖4 河網柵格數閾值—河網密度關系曲線

本文引入多種函數關系擬合河網柵格數閾值和河網密度，最終確定采用擬合效果最好的冪函數關系，得到擬合方程

y=6.165 7x-0.465(R2=0999 3)

(1)

式中，y為河網密度，km/km2；x為河網柵格數閾值；R2為相關系數，代表擬合程度，R2=0.999 3表明擬合程度高，兩者相關性強。

如圖5所示，連接閾值樣本首尾兩端點做割線，以割線斜率為擬合曲線切線斜率做切線，求得切點即為擬合曲線由劇烈到平緩變化的轉折點。對擬合方程(1)求導，得到擬合曲線一階導數

圖5 最佳閾值的確定

k=y′=-2.867x-1.465

(2)

進一步求解方程(2)與斜率k的關系，獲得切點坐標為(15 933，0.068)，此時河網柵格數閾值15 933為提取水系的最佳閾值，即最佳集水面積閾值約為129.1 km2，以該閾值提取水系結果如圖6所示。切點縱坐標0.068其物理意義為，擬合曲線上河網柵格數閾值15 933對應的河網密度是0.068 km/km2；利用Arcgis水文分析工具以15 933為閾值提取水系，計算得到河網密度為0.069 km/km2，相對誤差為1.5%，擬合效果好。

圖6 90 m分辨率下瀾湄流域最佳集水面積閾值提取水系

為驗證提取水系精度，反映提取水系與矢量化數字水系偏移程度，對河流總長的相對誤差和河網“套合差”兩個指標進行計算分析。提取水系河流總長52 061.6 km，數字水系河流總長48 395.5 km，相對誤差為7.5%?；贏rcGIS圖層加載1∶1 000 000矢量化數字水系圖作為標準做檢驗，對比數字化水系和提取水系之間的位置偏移，統計由此偏移產生的兩水系之間的細碎多邊形面積(見圖7)。統計結果顯示，整個瀾湄流域內的細碎多邊形面積約為20 000 km2，河網套合差為2.5%，小于3%，吻合程度較好。由圖7a可見，上游地區海拔較高、地勢變化明顯且無較大湖泊，兩水系疊加形成的多邊形相對較少，多邊形面積約為2 100 km2，對應的流域面積為166 000 km2，河網套河差為1.3%；圖7b顯示為流域中下游地區，該范圍內地勢逐漸變緩，水系豐富，在柬埔寨境內有較大湖泊洞里薩湖，兩水系之間的多邊形較上游明顯增多且面積增大，計算得中下游的河網套河差為2.9%，上游和中下游的河網套河差結果均在有效范圍內，但中下游地區河網偏移程度較上游大，提取水系精度略差于上游地區?？傮w而言，本文提取水系結果與實際水系較為吻合，能有效、客觀地反映瀾湄流域水系整體特征，可為后續研究提供數據支撐。

圖7 提取水系與數字化水系河網套河差

3.3 河網分級和子流域劃分

由90 m分辨率高程數據提取最優集水面積閾值的瀾湄流域水系情況如圖6所示，河網分為1～6級，一級為河道最低級別的河源水系，共計1 342條，占總河道數的50%，一級河網發育系數約為2.9，河系不均勻系數約為1.3。

基于Arcgis自動化提取流域出口，共劃分小流域3 259個，流域面積普遍在50～900 km2范圍。為滿足業務應用需要，依據水系上下游匯流關系，建立河道水系與小流域以及各小流域間的拓撲關系；參考流域內地形特征以及湖泊范圍，滿足水系連續性、分水線完整性、出水口準確性等關鍵條件；對流域面積在5 km2以下的微小流域或中間匯流區，以匯流關系為主要依據合并到相鄰小流域內[31-32]；另外，又以流域內干流上的8個主要水文站點作為流域出口對其進行流域劃分；進而完成人工合并小流域的工作。將瀾湄流域的子流域劃分出3個等級，根據子流域面積由小到大依次定義為一級、二級、三級流域。一級流域共計9個子流域，子流域面積最小的是昌都—舊州區間，約為1.6萬km2，最大的是穆達漢—上丁區間，約為24萬km2。二級流域共計61個子流域，子流域面積小于5 000 km2的占比25%，在5 000～10 000 km2的占23%，在10 000～20 000 km2的占33%，大于20 000 km2的占19%。三級流域共計328個子流域，子流域面積普遍控制在500～5 000 km2，占比達到90%。其中，小于500 km2的為不適合與相鄰合并的中間匯流區，另存在一個大于5 000 km2的子流域，受洞里薩湖水域限制，包含湖域區域，子流域面積約6 700 km2。

3.4 不同子流域尺度流域面雨量預報分析

選取2021年7月24日～25日的一次降水過程預報為例，即為24日8時起報的24 h面雨量，分析不劃分子流域的瀾湄流域以及劃分出一級流域、二級流域的面雨量結果對此次過程在空間分布上的響應。未劃分子流域的全流域經算術平均計算面雨量為25 mm，全流域面積大僅以一個數值表達面雨量，在降水分布均勻且降雨量較小時可起到一定參考作用，但在降雨空間分布不均時無法體現降雨中心。圖8a中一級流域面雨量計算結果顯示，昌都—舊州和穆達漢—上丁面雨量量級為小雨，舊州—允景洪量級為中雨，允景洪—清盛、清盛—瑯勃拉邦和瑯勃拉邦—萬象段量級為大雨，舊州—允景洪和萬象—穆達漢量級達到暴雨。圖8b顯示，二級流域進一步縮小子流域面積，面雨量結果與智能網格預報降雨結果在暴雨中心的空間體現上更為一致，降雨較大范圍集中于中下游地區，泰國中東部和老撾中南部交界的位置，有6個子流域面雨量超過60 mm，量級為大暴雨，最大在泰國境內達到143 mm，與一級流域萬象—穆達漢段暴雨級面雨量相呼應。

圖8 面雨量空間分布

4 結論與討論

本文分析了瀾湄流域集水面積閾值與水系特征變化的數量關系，建立了集水面積閾值與河網密度之間的冪函數，通過曲線割線斜率法客服主觀性，確定了90 m分辨率高程下的最佳集水面積閾值為129.1km2。并以1∶1 000 000的矢量化數字水系圖做驗證依據，總河長相對誤差為7.5%，河網套河差為2.5%，均說明提取水系與實際情況吻合程度較高，提取水系有效，進一步劃分了瀾湄流域1～6級的水系和3個等級的子流域?；谥悄芫W格預報降水，以一次降水過程展現面雨量在不同等級子流域上的空間分布，子流域等級越大、劃分越細致對暴雨中心的體現越明顯。本文的研究成果，可在水文氣象業務上為瀾湄流域的面雨量監測、預報和水文模型的構建提供參考，為數字孿生流域建設與流域防洪減災奠定基礎。

雖然本文對DEM進行了修正工作，但是在處理較大湖泊、入?？谌侵薜忍厥獾匦紊吓c真實情況方面仍存在差異；研究區域面積大，在高精度高程數據情況下運行速度慢，選取的集水閾值數據量有限。下一步工作可以考慮按照地形、水系類型等因素先將整個流域劃分為幾個不同的小區域，以求在保證高程精度情況下提高運算效率，進而綜合模擬流域水系。