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

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

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

0 引 言

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

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

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

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

圖1 瀾湄流域概況

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

2 研究方法

2.1 河網(wǎng)提取方法

2.1.1 河網(wǎng)修正Agree算法

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

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

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

2.1.2 流域水系提取

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

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

2.2 確定最佳集水閾值方法

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

擬合集水面積閾值與河網(wǎng)密度關系曲線，構建集水面積閾值與河網(wǎng)密度關系方程f(x)，定義域為[x1，x2]，值域為[f(x2)，f(x1)]，令f(x)的一階導數(shù)為f′(x)。在定義域內(nèi)，曲線割線的斜率k=[f(x1)-f(x2)]/(x1-x2)，以該斜率作為擬合曲線的切線，即k=f′(x0)，求解x0，切線與擬合曲線的切點(x0，y0)代表了河網(wǎng)密度隨集水面積閾值由劇烈變化變?yōu)槠骄徸兓霓D(zhuǎn)折點。即，x0為目標閾值[16]。

2.3 河網(wǎng)精度評價

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

圖3 河網(wǎng)套合差示意

2.4 流域面雨量計算

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

3 結(jié)果與分析

3.1 集水面積閾值與水系河網(wǎng)的變化關系

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

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

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

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

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

圖4 河網(wǎng)柵格數(shù)閾值—河網(wǎng)密度關系曲線

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

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

(1)

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

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

圖5 最佳閾值的確定

k=y′=-2.867x-1.465

(2)

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

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

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

圖7 提取水系與數(shù)字化水系河網(wǎng)套河差

3.3 河網(wǎng)分級和子流域劃分

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

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

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

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

圖8 面雨量空間分布

4 結(jié)論與討論

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

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