基于GIS和DEM數(shù)據(jù)的流域水系提取與優(yōu)化分析

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

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

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

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

1 研究區(qū)概況

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

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

2 研究方法與過程

2.1 數(shù)據(jù)及預處理

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

圖2 研究區(qū)DEM影像圖Fig.2 DEM image of study area

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

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

2.2 流域水系的提取

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

圖4 水系提取實現(xiàn)技術(shù)路線Fig.4 Technical route for basin river extraction

2.2.1 研究區(qū)地形分析

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

(1)

A=Sx/Sy

(2)

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

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

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

2.2.2 水流方向確定

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

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

2.2.3 匯流累積量的計算

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

2.2.4 流域水系的提取

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

3 結(jié)果分析與優(yōu)化

3.1 結(jié)果驗證與誤差分析

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

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

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

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

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

圖8 不同地區(qū)DEM數(shù)據(jù)提取誤差分析Fig.8 Error analysis of DEM data extraction in different regions

3.2 網(wǎng)格數(shù)的確定與改進

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

圖9 距離誤差隨網(wǎng)格數(shù)的變化曲線Fig.9 Range error curve along with grid numbers

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

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

3.3 集流閾值格數(shù)的合理取值

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

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

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

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

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

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

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

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

4 結(jié) 語

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

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

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

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

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

□