基于GIS和DEM數據的流域水系提取與優化分析

王成文，李 英,2,3，黃小琴，張 勃，徐兆祥，李 陽

(1.寧夏回族自治區水文環境地質勘察院，銀川 750011；2.中國地質大學(北京)水資源與環境學院，北京 100083；3.寧夏回族自治區地質局，銀川 750021)

近年來，隨著GIS技術的不斷發展，區域地形地表水系的提取研究逐步由傳統的野外測量和地學統計方法轉向高精度的DEM(數字高程模型)數據數字化和自動化提取工作中來。在流域水系數字化提取研究中，如何利用DEM數據更精確真實地反映實際地表水系特征，是當前水文模型研究熱點和難點，同時水系的數值化提取研究也是數據信息反演解譯實際的重要組成部分[1]。目前，對于流域水系的自動提取研究方法豐富多樣，總體可分為兩種思路，一種是基于已有遙感圖像進行自動提取，另一種是利用DEM數據進行水系自動提取。由于DEM數據分辨率不斷更新提升，且能夠表達豐富的地貌信息而被廣泛研究和應用[2]。

大量研究實踐發現[3-8]，在實際地形中凹陷洼地和平緩地帶的是普遍存在的，DEM數據中對于柵格單元局部凹陷和平坦地帶細節處理不夠完善，因此在提取水系之前必須要對DEM進行預處理。此外，基于GIS技術提取水系過程中，如何合理的設定網格數和閾值范圍是關乎河流數字化成敗的關鍵因素[9-11]，而以往研究中大多是進行網格和閾值的大范圍值域的模擬、篩選，甄選出最佳值。本文在參閱國內外文獻和反復實踐的基礎上，利用ArcGIS10.4軟件中水文分析模型，進行基于DEM數據的流域水系提取與優化處理分析，探討流域水系提取的優化方法，特別是對網格數和閾值的選擇進行數學擬合分析，選擇合理的網格和閾值，為流域水系特征的數字化提取研究提供一定參考。

1 研究區概況

本文擬以寧夏海原縣境內范圍為研究區，進行山區流域水系提取。海原縣位于寧夏回族自治區中南部，六盤山西北麓，行政區劃屬中衛市管轄，地理位置介于東經105°09′～106°10′，北緯36°06′～37°04′之間，東西寬約100 km，南北長約95 km，面積4 990 km2，海拔在1 400～2 800 m之間，最高點位于南、西華山山區，最高海拔2 954 m，見圖1。

圖1 海原縣等高線示意圖Fig.1 Contour lines diagram in Haiyuan county

2 研究方法與過程

2.1 數據及預處理

本文中的DEM數據是源于“地理空間數據云”信息中心的高程數據，它是集結了ASTER GDEM第一版本(V1)的數據連接而成，數字高程分辨率為30 m。由于ASTER GDEM V1數據中有個別數據異常現象，需結合美國NASA的SRTM高程數據，分辨率為90 m。研究中以30 m分辨率的DEM數據為主，90 m輔助使用，影像圖見圖2。

圖2 研究區DEM影像圖Fig.2 DEM image of study area

由于DEM數據是由等高線生成，使得柵格數據中帶有一個或一組凹陷點，表現為四周高中間低的洼地，在進行水系提取前必須填充，對DEM數據進行預處理，稱為“填洼”處理[12]。利用ArcGIS中水文分析的Fill Sink工具進行填洼，先對每一個柵格單元進行掃描、收集，然后運用Zlimit選項確定合適的填充閾值，該值即為被填補深度的臨界值(地形標高)，只有當洼地低于該臨界值時才被識別、填充，而超過該臨界值的地形標高保留不變以達修復目的，修復結果見圖3。

圖3 DEM的預處理影像Fig.3 DEM preprocessing images

2.2 流域水系的提取

基于處理后的DEM數據，進行水系提取工作，采用目前廣泛應用D8單流向算法，該方法是通過DEM數據判別流域水系水流方向，計算流域匯流累積閾[13]。水系提取過程：①地形分析，②水流方向判別，③柵格數的確定和匯流累計量算，④設定合理集流閾值，提取水系信息，⑤柵格水系矢量化。水系提取技術路線見圖4。

圖4 水系提取實現技術路線Fig.4 Technical route for basin river extraction

2.2.1 研究區地形分析

坡度、坡向的地形特征是流域水系流向分析的基礎。坡度表示地表傾斜程度，坡向是表征某點高程值變幅的最大變化方向[14]。一般運用擬合曲線面法求解，并結合ArcGIS中3D analyst模塊對坡度、坡向進行分析、提取，公式如下：

(1)

A=Sx/Sy

(2)

式中：S為坡度；A為坡向；Sx為x方向坡度；Sy為y方向坡度。

坡度、坡向提取結果見圖5，研究區內坡度在0.01°～36.72°之間，由西南往北東坡度呈減小趨勢，到東北一帶坡度<20°，坡向取值在0～360°之間。從坡度坡向可以看出，研究區坡度較大、坡向多變，地形極其復雜，中部、南部等大部分地區為縱深很長的山地，坡向清晰可辨；東部坡度較為平緩，辨析度降低。

圖5 研究區坡向坡度影像分析Fig.5 Image analysis of slope direction and gradient in study area

2.2.2 水流方向確定

運用D8算法[11,15]來確定水流方向，首先計算DEM數據每個柵格單元與四周的坡度和坡向關系，然后選擇最陡坡度，設定最陡坡度為該單元的水流流向，每個柵格單元中的水流向四周有8個方向流出可能，并用1、2、4、8、16、32、64、128這8個流向編碼分別表示各個流出方向，通過對比坡度，選擇某一方向，確定水流方向，水流流向編碼見圖6 ，其最大優勢是在基于ArcGIS 等的環境下能夠較為快速、準確地運算。

圖6 水流流向編碼Fig.6 Flow direction code

2.2.3 匯流累積量的計算

利用ArcGIS 水文分析模塊下的Flow Accumulation函數確定流域水系的匯流量，設置柵格臨界值，識別有效柵格數，計算出有效柵格上累積的匯流柵格數，確定匯水面積，匯流累積量為柵格單元數目與柵格單元面積乘積之和[16,17]。計算識別的有效柵格數的多少可以表征流域水系的匯流能力，即某流域單元內上的匯流柵格數越大，形成的匯流面積越大，形成地表的徑流越明顯。

2.2.4 流域水系的提取

科學確定符合研究區地形條件的閾值是提取流域水系的關鍵,以往研究一般設定閾值主要是通過反復對比不同閾值提取的流網，主觀判斷最佳閾值[18]。本文水系的閾值設定是根據需要對不同級別的河流設定不同閾值，要綜合考慮研究區的基本狀況、流域地形地貌等方面，分別將匯流量閾值范圍為2 500～25 000，提取河網水系，并運用多種數學函數擬合確定合理閾值提取水系。

3 結果分析與優化

3.1 結果驗證與誤差分析

本文以野外實測數據為準值，通過野外實地勘察測量，室內數字化得到河流的中心線，并以該中心線的長度為真實河長。由實測結果可知，海原縣境內共發育清水河和祖歷河兩大水系，均屬黃河支流，按其支流集水區域的不同，自北向南又可分為西河、莧麻河和楊明河等流域區。

本文以海原縣境內河流為研究對象，結合ArcGIS中的水文分析模塊繪制其流域信息，其流域長度見表1。在基于30 m DEM數據基礎上，利用D8算法提取了4級水系，其中一級水系為11條，二級水系為15條，三級水系為8條，四級水系均為11條。提取的水系總長度為350.05 km。為了定量比較提取的結果，表1給出了基于30 m DEM數據提取值和野外實測值的數據對比，研究區流域水系提取效果見圖7。對比可得，基于GIS和DEM數據的流域水系提取計算結果與野外實測量算成果較為接近，尤其大中河道長度偏差僅在7.5%以內，中小河流河道長度偏差在9.2%以內，而小河流河道長度的偏差值可達20%以上，可見提取的低級別水系數目的誤差較大。總之，從提取的水系長度及數量看，基于GIS和DEM數據提取的流域水系有一定參考意義。

表1 水系提取與高精度DEM測值對比Tab.1 Comparison of basin river extraction and high-precision DEM measurements

圖7 研究區流域水系提取影像圖Fig.7 Image of basin river extraction in the study area

另外，運用SPSS統計軟件分析，對比研究區提取的36條河流，發現DEM數據在流域水系提取應用中，在地形起伏較大的山地、山谷區域的提取識別較好，誤差也小，一般水系提取的誤差小于5.32%；而在起伏變化不大的丘陵、平原地帶提取識別效果一般，且誤差變大，提取的誤差值一般介于5.5%～12.8%，個別河流誤差大于20%，見圖8。

圖8 不同地區DEM數據提取誤差分析Fig.8 Error analysis of DEM data extraction in different regions

3.2 網格數的確定與改進

合理識別有效DEM數據的柵格單元數是影響匯流網格數的主要因素，而網格數的確定是提取水系流域面積的決定性因素，真實流域的水系與基于GIS提取水系的距離誤差是隨著網格數的變化而變化，通過擬合網格數與距離誤差的函數，得到距離誤差最小值所對應的網格數，并通過此網格數提取流域水系。在2 000～40 000網格數范圍內，選取若干網格數分別生成水系圖。在Arcgis10.4下加載DEM數據和提取的水系圖，逐一量測提取水系河源與實際河源坐標的距離誤差，從而得到距離誤差與網格數的相關關系，如圖9所示。

圖9 距離誤差隨網格數的變化曲線Fig.9 Range error curve along with grid numbers

由圖9可知，當網格數在10 000～15 000范圍時，實際河源與提取水系河源之間的距離誤差變化呈最小趨勢，而當網格數從15 000左右繼續增加時，距離誤差又呈增大的變化趨勢。為了進一步確定最小誤差及其對應的網格數，通過Origin9.0中的函數擬合，運用多種函數及偏導數的計算，最終選取采用多項式擬合，函數最次冪為3次，相關系數R=0.976。同時，求取函數式的一階偏導數，并令其一階偏導數等于0，再根據網格數的變化區間，確定誤差最小時的網格數為12 245，即網格數為12 245時所對應的水系即為該區的水系(見表2)。

表2 距離誤差擬合曲線參數Tab.2 Range error fitting curve parameters

3.3 集流閾值格數的合理取值

利用Hydrology水文處理工具，分別設定2 500，5 000，8 000，12 000，18 000，25 000 共6個匯流累積量閾值，生成柵格河網，見圖10，并結合不同閾值的流域水系特征，對比實際水系流量，計算其河流長度、匯流面積與河源數等水文信息，運用數學函數擬合分析，綜合確定合理閾值。

圖10 不同閾值條件下的流域水系提取Fig.10 Basin river extraction under different threshold value conditions

以研究區南部水系的楊明河為例，從累積匯流量圖層中提取各閾值下的流域水系，計算不同閾值下的河流長度、匯流面積與河源支流個數等水文信息，詳見表3。隨著閾值的不斷增加，河流長度、匯流面積與支流個數均呈減小的趨勢，且不同水文參數的變化趨勢各有差異。結合各閾值流域特征擬合曲線分析(圖11)，當閾值在8 000之內時，河流長度隨閾值的增加呈驟降趨勢，由464.36 km急降到88.72 km；當閾值大于8 000時，河流長度隨閾值的增加呈緩慢減少趨勢，表明河網密度變化較平緩，總體變幅較小。當閾值在12 000內時，流域面積隨閾值的增加呈驟降趨勢，由3 545.32 km2驟降到552.52 km2；當閾值大于12 000 時，流域面積隨閾值的增加呈緩慢減少趨勢，表明河流面積變化較平緩，變幅較小。同樣，河流支流個數同樣呈相識的規律變化，在閾值12 000為河流支流個數的變化拐點，呈“驟降-緩減”的變化趨勢。可見閾值在8 000和12 000為該研究區的關節點，尤其是閾值12 000，對整個河流的影響呈顯著影響。

表3 不同閾值的流域水系特征(楊明河流域)Tab.3 Basin river characteristic with different threshold values (Yangming River)

圖11 不同流域特征因數隨閾值變化的擬合曲線Fig.11 Fitting curve of different basin characteristic factors with threshold value

另外，運用函數擬合手段，分別采用冪函數、指數函數、多項式等多種函數對河流長度、流域面積及支流個數進行趨勢擬合分析，以確定提取水系流網的最佳閾值。通過多次試驗對比，發現河流長度和支流個數采用冪函數(y=AxB)擬合效果較好，相關系數分別為0.981和0.973，曲線的閾值變化拐點范圍均為8 000～12 000；而流域面積則更適合采用指數函數擬合，其相關系數為0.998，擬合效果最好，曲線的閾值變化拐點范圍同樣是在8 000～12 000，詳見表4。

表4 不同流域特征因數隨閾值變化的擬合參數Tab.4 Fitting parameters of different basin characteristic factors with threshold value

結合研究區地形氣候特征和實際測量分析，海原地區為西北干旱山區，常年地表蒸發遠遠大于降水補給，易選取閾值較大值，即選擇水系流網較稀疏，水流面積較小值，更符合干旱山區的實際水系特征。由此可以得出，選取的閾值設為12 000時，模擬的海原地區水系河網與實際水系吻合效果最佳，同時更科學地指導基于DEM數據源的水系和實測水系之間的差異分析。

4 結 語

(1)基于ArcGIS中水文模型和DEM數據的技術分析，實現了海原縣地區流域水系的提取過程，形成了山區水系提取模型，并優化了水系提取方法。

(2)結合研究區域地貌特征，對DEM影像的坡度圖進行相應的處理，基于ArcGIS拓展模塊3D analyst中的坡度、坡向分析工具，對研究區進行坡度、坡向分析，發現山區的坡向辨識度更清晰，而丘陵、平原的坡向辨識度依次降低。

(3)結合研究區水文流域信息，運用多種數學函數分析，科學合理地選定網格數網格數為12 245，設定閾值為12 000，提取出的流域河網水系和實際情況吻合度最高。

(4)在研究網格數對區域水系提取的誤差影響中，發現3次冪多項式擬合曲線最佳，為今后山區水系提研究提供一定參考，同時其擬合效果的穩定性也需要更多的實踐驗證。

(5)本文僅提取了流域水系的河流長度、匯流面積和支流個數參數，將來可以考慮多重因素，提取更加豐富的地形與水文信息。

□