ArcGIS中基于DEM提取溝道特征

嚴建鋼，金復鑫，2，周小程，馬海洋

（1.海軍航空工程學院指揮系，山東煙臺264001；2.陸軍航空兵學院，北京101123）

從數字高程模型（Digital Elevation Model，DEM）生成的積水流域和河道溝系網絡數據是大多數地表水文分析模型的輸入數據[1]。因為DEM 中包含有各種分辨率的地形高程信息，可以提取出大量地表形態信息，如流域的坡度、坡向和流域單元格的流向等。由DEM 提取地形特征的研究在國外從20 世紀60年代就已經開始，80年代之前研究比較少，范圍也僅限于分水線和山谷的識別和提取；其高峰出現在80～90年代，相繼出現了各種提取河網、流域邊界以及劃分子流域的方法[1]。

隨著流域特征提取研究方法的不斷成熟，各種專業軟件和計算機模型相繼推出，所有這些軟件和模型中，尤以ESRI 的ArcGIS 最為成熟，應用也最為普及。因而本文以ArcGIS9.3為軟件平臺對研究區進行流域溝道特征的提取。

1 資料與方法

1.1 研究區概況

大南溝流域地處陜北黃土丘陵溝壑區，位于安塞縣城西北約7.5 km 處，屬延河一級支流，流域面積約3.6 km2，海拔1 100～1 327 m。流域屬暖溫帶半干旱氣候區，年均溫9 ℃，年均降雨量549 mm。流域內地形破碎，溝壑密度達6.9 km/km2；地勢起伏率和坡度大，溝峁頻繁相間，土壤侵蝕強烈，塑造成典型的梁狀黃土丘陵地貌和復雜多樣的土地類型；流域內部地勢相對開闊，多平緩塌地，有大面積的梯田。流域內絕大部分土壤是黃綿土，還有少量二色土和紅膠土[2]。

1.2 數據源

本文用于提取流域特征的是大南溝流域5 m分辨率的DEM數據[3]。在這里，DEM數據是通過數字化研究區1∶10 000 地形圖，再插值生成TIN 數據，然后由TIN數據生成ESRI的grid格式柵格數據而成。在X方向上有580個柵格，在Y方向上有519個柵格，Z方向上有1個柵格。

1.3 提取過程

ArcGIS9.3 中沒有可視化的水文分析模塊，其相關的功能要通過在Grid 模塊中調用一系列函數來實現。圖1的流程圖簡要地說明了分析的過程。

圖1 ArcGIS進行水文分析流程圖

1.3.1 數據預處理

預處理主要是對原始DEM中存在的洼地進行填平處理。洼地是指高程低于周圍柵格的一個柵格或空間上相聯系的柵格的集合，洼地的出現可能是由DEM在數據采集和插值時的誤差造成的，也可能是真實地形特征的反映（如喀斯特地區、采石場等）[4]。由于洼地被較高的地形包圍，不能確定其水流方向，因而是進行水文分析的一大障礙，必須將其填平[5]。

在ArcGIS 中對洼地的填平處理包括以下步驟[6]：首先，基于原始DEM 數據生成流域的流向數據；然后，結合流向數據調用SINK 函數對洼地進行標定。SINK函數的功能就是標定出流域內高程低于周圍的封閉區域。語句格式為SINK（＜flowdir＞），flowdir 為流向數據。

標定出洼地之后要確定洼地的影響區域，即有多少柵格的匯水到達這些洼地中。本步驟的實現調用Watershed函數，該函數需要流向數據和上一步中的洼地作為輸入。語句格式為watershed（＜flowdir＞，＜sinks＞），sinks為標定的洼地柵格。

下一步要結合高程數據和影響面積數據分別確定洼地的最低高程（洼地底部）和最高高程（洼地邊緣的最低高程），需要分別調用zonalmin 和zonalfill 函數。函數的語句格式分別為zonalmin（＜sink_areas＞，＜elevation＞）和zonalfill（＜sink_areas＞，＜elevation＞）。

完成上一步之后要確定需要填充洼地的最深值，即超過這一深度的洼地被認為是合理的地形特征，而不再進行填充，在此深度范圍內的洼地才填充，該值用洼地最高高程柵格減去最低高程柵格即得。

最后一步就是根據以上步驟的結果對原始數據進行填洼，調用的函數是Fill。Fill 的工作原理是：掃描各個單元格的時候，比較該單元格與相鄰的8 個單元格的高程，如果是洼地，那么該單元格的高程值將被賦予相鄰8個單元格中高程最低的那個；Fill函數的語句格式為FILL＜in_grid＞＜out_grid＞{SINK|PEAK}{z_limit}{out_dir_grid}，in_grid 為待處理DEM 數據，out_grid 為處理后DEM 數據的路徑和文件名，SINK PEAK是選擇填洼或是削峰，在此處選擇填洼，z_limit即為填充洼地的最深值，超過此深度不再填充，out_dir_grid為生成的新流向數據。

因為在填洼的同時改變了原有的柵格高程，可能在洼地的邊緣處產生新的洼地，所以以上填洼過程要反復進行，直到用SINK 函數標定不出洼地時才可以認為填洼完成，生成了無洼地DEM數據。

一般認為，公司的盈余質量、信息披露質量、會計穩健性等越好，財務報告質量越好。目前，以審計委員會中獨立董事的比例度量審計委員會獨立性的研究，發現審計委員會獨立性越高，內部控制質量越好（謝海娟等2016、陳文娟等2016），盈余質量越好（王振秀2017、李建紅2016），會計穩健性越好（謝香兵2011），信息披露質量越好（蔡衛星等2009、柯明等2011、劉彬2014），審計意見越好（何衛紅2016），財務報告質量越好（周國華2011、周國平等2013）。另外，潘珺等（2017）發現審計委員會獨立董事和非執行董事比例之和越高，公司的財務報告質量越好。

1.3.2 流向數據的生成

完成原始DEM 數據的填洼預處理之后，則可基于無洼地DEM數據生成流向數據。

在ArcGIS 中生成流向數據用的是D8 算法[7]，該方法假設單個網格中的水流只有8 種可能的流向，用最陡坡度法來確定水流的方向，每個柵格單元水流的流向為其鄰近8 個柵格單元中坡度最大的那個單元格，坡度的計算公式為θ=arctan[(hi-hj)/D]，其中hi和hj為2 個單元格的高程值，D為2 個單元格中心的距離，當相鄰柵格為共邊時，D為柵格像元的尺寸，當相鄰柵格為對角時，D為 2倍的柵格尺寸。

在ArcGIS 中調用Flowdirection 函數生成流域8種方向值的流向數據，函數的語句格式為flowdirection（＜elevation＞）。圖2 所示為流向柵格數據每個數值所代表的方向。

圖2 柵格流向示意圖

1.3.3 累積匯水面積數據生成

每個柵格單元的累積匯水面積表現了該柵格匯集上游來水的能力，匯水能力越大的柵格就越可能是河道。在ArcGIS中，這一過程被稱為求柵格的累計流量[8]，通過調用Flow accumulation 函數來實現。此功能的原理是假想在集水區的每一網格上降下一單位的水量，而后按網格的流向來向下移動，其移動經過的網格則使其累積流量值提升一個單位，因而每一網格都能計算出其所累積的上游流量值。由于投入每一網格的水量皆為一單位，故流量累積值亦代表各網格的上游集流網格數量，將之乘上網格面積便可得到每一網點的上游集水面積[9]。函數的語句格式為flow accumulation（flowdir）。

1.3.4 臨界匯水面積的確定和溝道提取

臨界匯水面積，即CSA（Critical Source Area），是區分河道柵格和流域內非河道柵格的臨界值，大于該數值的柵格被賦值為1（溝道），小于該數值的賦值為NULL。這一值在不同的氣候帶和土壤覆蓋條件下是不同的，如果取值太小，則會在提取中出現大量偽溝道，如果取值太小則會忽略真實水系。國內外很多學者對CSA的取值進行了研究，有些學者認為在地形復雜區應根據地貌特征選擇不同的值，國內一些學者采取與光照暈渲圖擬合的方法確定CSA取值[10]，本文即采用此種方法來確定臨界值。

確定出匯水臨界面積之后，要構造條件語句進行流域匯水能力數據的二值化，即流域溝道的提取。在Grid 中調用con 函數，語句格式為con（＜condition＞，＜true_expression＞，{＜condition＞，＜true_expression＞}，…{＜condition＞， ＜true_expression＞}， {false_expression}）。其中{false_expression}為可選項。可選項為空時，如果＜condition＞語句返回的值為假，那么單元格將被默認賦值為NULL。例如：con（flow_accumulation＞Const，1）表示將柵格數據flow_accumulation 中的匯流的值大于某一常數的全部賦值為1，其余賦值為NULL。

1.3.5 溝道的分級

對溝道的分級有2 種方案[10]：Strahler 水系系統分級方案和Shreve 水系系統分級方案。前者的原理是以源頭溝道開始為第1 級，相同級別的溝道相交后級別增1級，不同級別溝道相交后級別不增，保持原溝道中較高的級別，依次向下游追蹤；后者考慮河網中所有的結點，最外層溝道賦級別為1，向下游追索到的結點級別為上游級別的累加。分級要調用streamorder函 數，語句格式為streamorder（streamnet，flowdir，strahler/shreve），streamnet 為河道數據，flowdir 為流向數據。strahler/shreve為提供備選的分級方案，本文選擇Strahler方案。

2 結果與分析

在本次研究中自動統計生成的關鍵的數據有3個：填洼之后生成的無洼地DEM 數據，基于無洼地DEM數據生成的流向數據和流域累計匯水面積數據。

2.1 無洼地DEM

在預處理中通過掃描原始DEM 數據標定了5 處洼地區域，最低的洼地底部高程為1 145 m，最高的洼地底部高程為1 195 m。計算得出的填充洼地的最深值為1，以此值為限制用Fill 函數進行填充，對填充后的DEM高程數據用SINK函數查找洼地得到空值，說明流域中影響流向確定的洼地已經被移除，見圖3。

圖3 流域無洼地DEM圖

2.2 累積匯水面積分析

流域累積匯水面積數據（見圖4）是進行溝道提取的基礎，本文結合生成的地形光照暈渲圖來設定合適的CSA 值，經過多次選值與暈渲圖疊加分析，認為當CSA 取300 時，所提取的溝道與暈渲圖的地形特征吻合較好（圖5 白色線狀部分為提取的溝道），可以把二值化得出的數據作為最終的溝道提取結果。

圖4 流域累積匯水面積圖

圖5 溝道與地形暈渲圖疊加示意圖

2.3 溝道分級結果

通過對提取出的流域溝道數據進行Strahler 方案水系分級，可將流域內的溝道系統分為4級，其中從源頭開始的1 級水系有58 條，一級水系交匯而成的2 級水系有15 條，3 級水系有3 條，4 級水系有1 條。水系的分級基本符合流域溝道的實際分布特征，見圖6。

圖6 溝道分級圖

3 結論

通過與大南溝小流域1∶10 000地形圖所反映的實際流域特征對比分析，可以認為提取出的溝道特征基本符合實際，即用ArcGIS的水文分析功能可以進行流域溝道特征的提取，這在水文分析和溝道系統研究領域有著重要的意義。

但在分析中也存在一些問題有待于做進一步的研究和探討，主要有以下幾方面：

1）填洼預處理過程有可能會導致實際地形的改變。在ArcGIS 中把洼地的成因歸為是由于對地形的過低估計而致，故而采用增加高程的填平方案，而實際上洼地還可能是由于周圍高程過高估計而成，或是（有坡度地區）被高地阻擋而成（阻擋型洼地），單純的填平會導致這類洼地實際地形被改變。

2）對平地的處理效果不明顯，在地形起伏小的地區，溝道的模擬會出現誤差，尤其是在填平處理后產生的平坦地區，流向的確定仍不理想。

3）CSA值的選取，本文采用地形光照暈渲圖做參照的方法來確定臨界匯水面積值，但這種方法的主觀性太強，缺乏定量的數學基礎，且光照暈渲圖只能模糊地反映地形起伏特征，以其為參照會產生誤差，故需要進一步討論更好的方法，以使提取的溝道能正確地反映實際狀況。

[1] 李麗，郝振純.基于DEM的流域特征提取綜述[J].地球科學進展，2003，4（18）：251-252.

LI LI，HAO ZHENCHUN. The automated extraction of catchment properties from digital elevation models[J].Advance in Earth Sciences，2003，4（18）：251-252.（in Chinese）

[2] 許明祥，劉國彬，溫仲明，等.黃土丘陵區小流域土壤特性時空動態變化研究[J].水土保持通報，2000，20（1）：21-23.

XU MINGXIANG，LIU GUOBIN，WEN ZHONGMING，et al. Temporal and spatial variation of soil characters in small catchment of loess hilly areas[J]. Bulletin of Soil and Water Conservation，2000，20（1）：21-23.（in Chinese）

[3] 趙健，賈忠華，羅紈.ARCGIS 環境下基于DEM 的流域特征提取[J].水資源與水工程學報，2006，17（1）：74-76.

ZHAO JIAN，JIA ZHONGHUA，LUO WAN. Extracted catchment properties from DEM using ARCGIS[J]. Journal of Water Resources & Water Engineering，2006，17（1）：74-76.（in Chinese）

[4] 肖飛，張百平，凌峰，等.基于DEM的地貌實體單元自動提取方法[J].地理研究，2008，27（2）：459-466.

XIAO FEI，ZHANG BAIPING，LING FENG，et al. DEM based auto-extraction of geomorphic units[J]. Geographical Research，2008，27（2）：459-466.（in Chinese）

[5] JENSON S K，DOMINGUE J O. Extracting topographic structure from digital elevation data for geographical information system analysis[J]. Photogram Metric Engineering and Remote Sensing，1988，54（11）：1593-1600.

[6] 張超，鄭鈞，張尚弘.ArcGis9.0中基于DEM的水文信息提取方法[J].水利水電技術，2005，36（1）：1-4.

ZHANG CHAO，ZHENG JUN，ZHANG SHANGHONG.Extraction of hydrological information from digital elevation model with ArcGis 9.0[J]. Water Resources and Hydropower Engineering，2005，36（1）：1-4.（in Chinese）

[7] 孫慶艷，余新曉，胡淑萍，等.基于ArcGIS 環境下DEM流域特征提取及應用[J].北京林業大學學報，2008，30（2）：144-147.

SUN QINGYAN，YU XINXIAO，HU SHUPING，et al.Extraction and application of hydrological information based on DEM in ArcGIS environment[J].Journal of Beijing Forestry University，2008，30（2）：144-147.（in Chinese）

[8] 李翀，楊大文.基于柵格數字高程模型DEM的河網提取及實現[J].中國水利水電科學研究院學報，2004，2（3）：208-214.

LI CHONG，YANG DAWEN. Deriving drainage networks and catchment boundaries from Grid Digital Elevation Model[J]. Journal of China Institute of Water Resources and Hydropower Research，2004，2（3）：208-214.（in Chinese）

[9] 易紅偉，湯國安，劉詠梅，等. 河網徑流節點及其基于DEM 的自動提取[J].水土保持學報，2003，19（3）：113-118.

YI HONGWEI，TANG GUOAN，LIU YONGMEI，et al.Stream runoff nodes and their derivation based on DEM[J].Journal of Soil and Water Conservation，2003，19（3）：113-118.（in Chinese）

[10]李金朝，國慶喜，葛劍平.基于DEM 的黃土高原溝壑區溝道系統的自動提取[J]. 西北林學院學報，2009，24（6）：220-223.

LI JINZHAO，GUO QINGXI，GE JIANPING，et al.DEM based automated extraction system of gully in the loess upland gully are[J].Journal of Northwest Forestry University，2009，24（6）：220-223.（in Chinese）