基于DEM的水系提取及集水閾值確定方法研究

羅大游, 溫興平, 沈 攀, 張皓楠, 張麗娟, 李進波, 郁 智

(1.昆明理工大學 國土資源工程學院, 云南 昆明 650093;2.賀州學院, 廣西 賀州 542899; 3.云南省礦產資源預測評價工程實驗室, 云南 昆明 650093)

基于DEM的水系提取及集水閾值確定方法研究

羅大游1,3, 溫興平1,3, 沈 攀1,3, 張皓楠1,3, 張麗娟2, 李進波1,3, 郁 智1,3

(1.昆明理工大學國土資源工程學院,云南昆明650093;2.賀州學院,廣西賀州542899; 3.云南省礦產資源預測評價工程實驗室,云南昆明650093)

[目的] 確定適宜的集水閾值，使自動提取的水系河網與實際河道相符。[方法] 以DEM數據為基礎，利用ArcGIS水文分析模塊對研究區流域河網水系進行自動提取，基于河網密度與集水閾值的相關性，在已提出的擬合函數一階導數求轉折點、二階導數求拐點法確定集水閾值的基礎上，以龍川江流域為例，提出求解河網密度變化率等于集水閾值變化率的數值方法，得到適宜的集水閾值。[結果] 通過設置不同集水閾值生成河網，發現不同集水閾值對主河道長度及地理空間位置影響較小，但對提取的河網特征影響較大，最終確定龍川江流域集水閾值設置為0.12 km2為宜。[結論] 集水閾值的確定影響著河網提取的精度，通過變化率確定集水閾值的方法主觀因素較少，避免了人為干擾，相對客觀，可對區域水土流失監測提供一定依據。

DEM; ArcGIS; 河網密度; 集水閾值

文獻參數： 羅大游, 溫興平, 沈攀, 等.基于DEM的水系提取及集水閾值確定方法研究[J].水土保持通報，2017,37(4)：189-193.DOI:10.13961/j.cnki.stbctb.2017.04.032； Luo Dayou, Wen Xingping, Shen Pan, et al. Information extraction of river networks and determination of drainage area threshold using DEM data[J]. Bulletin of Soil and Water Conservation, 2017,37(4)：189-193.DOI:10.13961/j.cnki.stbctb.2017.04.032

DEM(digital elevation model)是對地形表面形態數字化表達的一種數字高程模型，它在水文領域中有著廣泛的應用，是地形劃分、水文分析的重要基礎數據。通過利用ArcGIS水文分析模塊，結合DEM數據自動獲取包含河網在內的水文參數，對數據加以分析，最終實現地形模型可視化。目前水文特征提取與分析在水土保持中起到了重要作用，在暴雨洪水災害中，通過DEM生成的水文模型提取的河網作為研究區危險性評價指標[1-2]；在泥石流、滑坡等地質災害頻發，水土流失嚴重地區，為災后研究區的水土治理提供基礎依據[3]；在河流、湖水流域作為生態環境治理實踐數據支撐[4-5]。

以DEM數據為基礎，利用GIS空間分析方法提取水系以提取研究區水文特征。國內外在這方面已開展了大量的工作[6-8]。國內在集水閾值的確定方面也取得了長足的進步。傳統方法中孔凡哲等[9-10]觀察集水閾值與河網密度的相關性，取河網密度變化趨于平緩時對應的值作為集水閾值，此方法簡單，易于操作但人為主觀性較高；數學方法中葉章蕊等[11]通過擬合函數一階導數的求解得到轉折點。關穎慧等[12]利用擬合函數二階倒數求拐點以獲取集水閾值。該方法將閾值確定在較小范圍但可能出現擬合函數如冪函數無法求得轉折點、拐點的情況；其他方法中王林等[13]基于網格法、Horton定理的估算方法計算流域水系分維數，結合水系分維與集水閾值的相關關系，把分維值趨于平緩的點作為集水閾值，此方法同樣存在人為主觀影響。在前人的基礎上，提出基于河網密度與集水閾值的相關性，求解河網密度變化率等于集水閾值變化率的數值方法，將河網密度與集水閾值的冪函數同切線方程聯立，求解唯一值作為河網提取的閾值。此方法避免了水系特征提取的人為主觀性，具有速度快，成本低的優點，真實反映水系的空間分布規律的優點。

本文擬以龍川江流域為研究對象，DEM數據為基礎，利用ArcGIS水文分析模塊對流域進行水文特征提取分析，對了解地區地形發育特征、水土流失監測、地質災害提取等有著重要意義。

1 研究區概況

龍川江為金沙江南岸一級支流，水源豐富，發源于云南省楚雄彝族自治州西部，由西向東流經楚雄市，橫穿元謀壩區，最終在元謀北部的江邊鄉匯入金沙江。流域面積約6 109 km2，境內最高高程與入江口相對落差約2 000 m，流域平均海拔約1 850 m。該流域以山區和丘陵為主，屬低緯度高原季風氣候，受大氣環流和季風影響，干濕季分明，降水量年內分配不均勻，雨季6—10月降水量占年降水量的80%以上，多年平均降水量847.2 mm，年最大降水量為1 075.2 mm，最小降水量為634.6 mm。2016年9月17日，因持續強降雨，研究區發生近年來最大山洪泥石流災害，在龍川江內形成朱布、海洛2個堰塞湖，嚴重威脅著上下游區域人民群眾的生命財產安全，因此龍川江是云南省防汛抗災的重點區域。

2 數據來源及研究方法

2.1 數據來源

本文采用ASTER GDEM作為基礎數據。ASTER GDEM(星載熱發射和反射輻射儀全球數字高程模型)是由美國航天局(NASA)與日本經濟產業省(METI)共同推出的地球電子地形數據,根據NASA的新一代對地觀測衛星Terra的詳盡觀測結果制作完成的。其數據覆蓋范圍為北緯83°到南緯83°之間的所有陸地區域，達到了地球陸地表面的99%。目前，中國科學院計算機網絡信息中心科學數據中心已經加工產生了中國及周邊區域范圍內ASTER GDEM 30 m分辨率系列數據產品。從中國科學院國際科學數據服務平臺找到龍川江下游所在圖幅，選擇數字高程數據下載即獲得所需的DEM數據。

ASTER GDEM數據可免費獲取，為檢驗數據的可行性，對比分析了ASTER GDEM(裸露地表高程+植被和人工建筑物)與地形圖制作的DEM(裸露地表高程)。由于研究區地面多為低矮植被，結果是二者的差距不明顯，高程信息基本一致，提取的水系格局也基本一致。因此可用ASTER GDEM替代地形圖制作的DEM提取水系特征。

2.2 研究方法

以DEM數據為基礎，利用ArcGIS水文分析模塊提取地表水流徑流模型的水流方向、匯流累積量、河網提取、研究區流域的劃分。通過對基本水文因子的提取分析，可以再現地表水的流動過程，最終實現研究區流域水文特征的數字化分析。圖1簡要描述了水文要素提取的基本流程。

圖1 水文信息提取流程

2.2.1 水流方向提取與無洼地DEM生成 水流方向是指水流經過每一個柵格單元時的指向，由計算中心柵格與領域柵格的最大距離權落差來確定。ArcGIS中水流方向是利用D8算法[14](最大坡降法)，用鄰域柵格的坡降來確定水流流向，坡降最大的單元格即為水流出的方向。

一般認為DEM是較光滑的地形表面模擬，但由于真實存在的凹陷部分(如喀斯特地貌)，可導致在進行水流流向計算時在研究區內得到不合理或錯誤的水流方向。因此在進行水流方向的計算前，應先對DEM數據進行填洼處理，得到無洼地的DEM[15]。洼地填充的基本過程是先計算出DEM數據中的洼地區域，然后計算出洼地區域的洼地深度，最后以最大洼地深度為參考而設定填充閾值進行研究區洼地填充，反復填洼直至無洼地DEM生成。

2.2.2 匯流累積量及河網提取 匯流累積量的基本思想是以規格網表示的DEM每點處有一個單位水量，按照自然水流由高處流向低處的自然規律，根據區域流向柵格計算每點流過的水量值，從而得到該區域的匯流累積量。當匯流量累積到一定程度的時候，地表就有水流產生，所有匯流量超過臨界值得柵格就是潛在的水流路徑，這些流水路徑構成河網。目前常用的河網提取方法是地表徑流漫流模型。計算步驟是以水流方向及匯流累積量為基礎前提，設置合理的集水閾值提取河網。目前較為常用“試誤法”確定集水閾值。閾值設定的大小影響河網的密度見圖2。閾值越大，河道數目越少，河網密度越小；閾值越小，河道數目越多，河網密度越大。

圖2 研究區河網隨集水閾值變化

2.2.3 流域的分割 流域是由分水嶺分割而成的匯水區域，通過對水流方向數據的分析確定出所有相互連接并處于同一流域盆地的柵格。流域的分割可以減少大尺度空間上不同因素差異導致的誤差，其步驟是確定小級別流域的出水口位置，然后分析找出所有上游出水點流過出水口的柵格，結合流向數據，確定集水區柵格的位置，從而實現子流域的劃分。

3 結果與分析

設置不同的集水閾值，獲取的河網也不同(圖2)。選取0.01～1.0 km2范圍內的13個集水閾值，分別計算出對應的河網密度和子流域數量，進而得出集水閾值與河網密度、子流域數量關系圖(圖3—4)。

從圖3—4中看出隨著集水閾值的變化，河網密度、子流域數量的變化趨勢大致相同。

圖3 研究區集水閾值與子流域數量的關系

圖4 研究區集水閾值與河網密度的關系

3.1 集水閾值對河網提取的影響

從圖2可看出集水閾值取值較小時，生成的河網密集，在平地區域容易生成平行狀河道(偽河道)，且生成大量破碎的子流域。隨著集水閾值取值逐漸增大，河網密度和子流域數量呈下降趨勢，但是降低的趨勢越來越緩慢。當集水閾值由0.01 km2變化至0.1 km2，關系曲線下降明顯。河網密度由0.834 km/km2減少到0.261 km/km2，子流域數量由4 686下降到464。當集水閾值從0.1 km2增至1.0 km2時，河網密度減少到0.087 km/km2，子流域數量下降到47。

3.2 集水閾值的確定

分別采用冪函數、指數、對數、線性、多項式、移動平均6種函數進行趨勢擬合，發現采用乘函數進行趨勢擬合時兩者的擬合度最好，得到的擬合方程如下：

y=0.084 6x-0.495(R2=0.999 8)

(1)

冪函數的一階導數、二階導數曲線關系見圖5。河網密度隨集水閾值的增大無限逼近于0，即不能通過冪函數的一階、二階求導的方法求出轉折點以提取集水閾值。因此本文提出用冪函數與直線相切的方法提取集水閾值。

圖5 研究區河網密度隨集水閾值變化冪函數的一階導數、二階導數曲線

研究中采用公式(2)為切線方程：

y=-x+a

(2)

式中：y——河網密度(km/km2)；x——集水閾值(km2)；a——截距。

冪函數與直線相切的點為拐點，拐點在數學上也稱為穩定點，出現拐點之前，集水閾值變化率Δx小于河網密度變化率Δy，反之Δy<Δx；嘗試求拐點即河網密度隨集水閾值變化達到穩定時，得到合理集水閾值。當x有唯一解時，Δy=Δx，所取的值即為最佳集水閾值。經計算，當a=0.359時，x有唯一解0.12；當a<0.359時，x無解；當a>0.359時，x有雙解。a取值0.359后，x<0.12時，Δy>Δx；x=0.12時，Δy=Δx；x>0.12時，Δy<Δx。通過上述確定集水閾值的方法與傳統方法相比，避免了集水閾值取值的主觀性，運用數學求值的方法簡化了反復試驗閾值的繁瑣過程。

3.3 提取河網與真實地形對比

隨著集水閾值的不斷增大，河道逐漸向地勢平坦處縮回，平行狀河道逐漸消失。當集水閾值為0.12 km2時，地勢平坦地區未出現平行線段，即偽河道全部被消除，生成的河網曲線較為光滑符合流域的河道走向。采用部分參考文獻[9-10]中河網密度變化趨于穩定提取閾值的方法，可將研究區的集水閾值定為0.2 km2；在本文運用的數值方法中求得的集水閾值為0.12 km2。通過與傳統方法比較，發現集水閾值為0.2 km2的河網中部分常流河道缺失，而本文提出的方法保持了常流河道的完整性，生成的河網水系與元謀縣中部SPOT6遙感影像吻合。

4 結果與討論

本文以龍川江流域為例，應用DEM和ArcGIS的水文分析模塊對研究區進行水文特征提取與分析。研究區地形起伏較大，在沒有高精度的DEM數據的情況下，提取的河網也能較好的與自然水系吻合。集水閾值對主干河道的空間地理位置影響小，主河流長度基本不發生變化，對研究區流域提取的河網特征影響極大。采用求解河網密度變化率等于集水閾值變化率的數值方法，確定研究區適宜的集水閾值為0.12 km2，進而得出自動提取的河網密度為0.239 km/km2，生成的河網與河網密度見圖6。從河網隨集水閾值變化趨勢可以看出，當集水閾值較小時河網密集，存在眾多小級別河道，集水閾值超過0.3 km2時小級別河道消失，主河道長度基本不發生改變，說明龍川江存在季節性河流。當集水閾值取值偏小時，提取的水文特征可用于檢測季節性河流，進而對地質災害點的提取起到積極作用。河網密度最密集的地區(圖6b)分別是元謀縣姜驛鄉、江邊鄉、黃瓜園鎮，這與地質災害統計數據元謀縣泥石流頻發區域相吻合；當集水閾值取值偏大時，可用以識別非季節性河流的長度與地理位置。

圖6 研究區河網與河網密度

不同DEM分辨率下提取的流域水系形態有差異，分辨率越高，河網越復雜精細，但河網結構大體上是一致的,提取的流域面積、河網密度相差不大[16]。隨著分辨率的降低，在地勢起伏較平坦地區河網特征變化大[17-18],研究區地形起伏較大，用低分辨率DEM替代高分辨率DEM，仍可取得良好的水系特征[19]。生成的河網存在部分缺失的情況與DEM數據來源有一定關系，因研究區海拔差異較大，且區內多為低矮植被，最終取得的水系特征良好。

綜上所述，利用DEM和ArcGIS進行的水文特征提取具有速度快，成本低，提取效果良好的優點。通過此方法可為研究區的地形研究、水土流失防治等工作提供研究基礎和參考依據。

致謝：本文得到中國科學院計算機網絡信息中心提供的DEM數據，昆明理工大學地質過程與礦產資源省創新團隊和昆明理工大學遙感地球化學學科方向團隊聯合資助，在此表示感謝！

[1] 張游,王紹強,葛全勝,等.基于GIS的江西省洪澇災害風險評估[J].長江流域資源與環境,2011(S1):166-172.

[2] 李慧敏,張建軍,黃明,等.基于DEM的黃土高原典型流域特征參數分析[J].北京林業大學學報,2012,34(2):90-95.

[3] 曾超.GIS支持下岷江上游水文特征定量分析[D].成都:四川師范大學,2011.

[4] 李晶,張征,朱建剛,等.基于DEM的太湖流域水文特征提取[J].環境科學與管理,2009,34(5):138-142.

[5] 宋向陽,吳發啟,趙龍山,等.基于DEM的延河流域水文特征提取與分析[J].干旱地區農業研究,2012,30(4):200-206.

[6] Murphy P N C, Ogilvie J, Meng Fanrui, et al. Stream network modelling using lidar and photogrammetric digital elevation models: A comparison and field verification[J]. Hydrological Processes, 2008,22(12):1747-1754.

[7] Ariza-Villaverde A B, Jiménez-Hornero F J, Gutiérrez de Ravé E. Influence of DEM resolution on drainage network extraction: A multifractal analysis[J]. Geomorphology, 2015,241:243-254.

[8] 任立良,劉新仁.基于數字流域的水文過程模擬研究[J].自然災害學報,2000,9(4):45-52.

[9] 孔凡哲,李莉莉.利用DEM提取河網時集水面積閾值的確定[J].水電能源科學,2005(4):65-67.

[10] 朱海玲,楊曉暉,張學培,等.基于DEM的密云水庫上游流域特征提取與分析[J].中國水土保持科學,2013,11(3):66-72.

[11] 葉章蕊,盧毅敏,張永田.基于曲線割線斜率法的水文特征提取[J].人民黃河,2016,38(2):28-31.

[12] 關穎慧,鄭粉莉,王彬,等.基于DEM的黑龍江賓州河流域水系提取試驗研究[J].水土保持通報,2012,32(1):127-131.

[13] 王林,陳興偉.基于DEM的晉江流域水系分維估算方法探討[J].地球科學信息,2007,9(4):133-137.

[14] O’Callaghan J F, Mark D M. The extraction of drainage networks from digital elevation data[J]. Computer Vision Graphics & Image Processing, 1984,28(3):323-344.

[15] Martz L W, Garbrecht J. An outlet breaching algorithm for the treatment of closed depressions in a raster DEM[J]. Computers & Geosciences, 1999,25(7):835-844.

[16] 王培法.柵格DEM的尺度與水平分辨率對流域特征提取的分析:以黃土嶺流域為例[J].江西師范大學學報:自然版,2004,28(6):549-554.

[17] 吳險峰,劉昌明,王中根.柵格DEM的水平分辨率對流域特征的影響分析[J].自然資源學報,2003,18(2):148-154.

[18] 林聲盼,荊長偉, MOORE Nathan,等.數字高程模型分辨率對流域地形特征參數的影響[J].水科學進展,2012,23(4):457-463.

[19] 康婷婷,王茜,趙苗琦,等.不同比例尺水系修正DEM的流域河網提取對比[J].地理與地理信息科學,2011,27(1):111-112.

Information Extraction of River Networks and Determination of Drainage Area Threshold Using DEM Data

LUO Dayou1,3, WEN Xingping1,3, SHEN Pan1,3, ZHANG Haonan1,3, ZHANG Lijuan2, LI Jinbo1,3, YU Zhi1,3

(1.Faculty of Land Resource Engineering, Kunming University of Science and Technology, Kunming, Yunnan 650093, China; 2.Hezhou University, Hezhou, Guangxi 5428999, China； 3.Mineral Resources Prediction and Evaluation Engineering Laboratory of Yunnan Province, Kunming, Yunnan 650093, China)

[Objective] Method of determining appropriate drain area threshold was demonstrated to reflect the actual river drainage network. [Methods] Based on the DEM data of Longchuan River basin, the automatic extraction module of watershed drainage network in ArcGIS software was used to extract relevant drainage information. After that, a function between river network density and catchment threshold was fitted, which was used to determine the appropriate drainage threshold by finding the turning point of its first derivative function and /or the inflection point of the second order derivative. This was simultaneously exemplified by the Longchuan River basin. In this case, the threshold was determined by finding where the ratio change of river network density equal to the change rate of the water drainage threshold. [Results] The river networks under different catchment drainage threshold were compared and found that no obvious differences existed with regard to the catchment threshold length of the main channel and the geographic location, but other features of river network were apparently different. In the case of the Longchuan River basin, its appropriate threshold was thought to be 0.12 km2. [Conclusion] The river network extraction precision varied with the drain area threshold. The proposed method in the determination of the catchment threshold is less subjective, and is more objective, which can be used to provide some bases for the monitoring of regional soil erosion.

DEM;ArcGIS;riverdensity;drainageareathreshold

A

： 1000-288X(2017)04-0189-05

： P627

2016-12-05

：2017-02-28

國家自然科學基金項目“天空開闊度對城市熱島效應的影響研究”(41101343)

羅大游(1993—),女(漢族),四川省成都市人,碩士研究生,主要研究方向為遙感地質。E-mail:746206934@qq.com。

溫興平(1970—),男(漢族),山西省興縣人,博士,教授,博士生導師,主要從事遙感地質及應用研究。E-mail:wfxyp@qq.com。