ArcGIS中基于DEM提取溝道特征

嚴(yán)建鋼，金復(fù)鑫，2，周小程，馬海洋

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

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

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

1 資料與方法

1.1 研究區(qū)概況

大南溝流域地處陜北黃土丘陵溝壑區(qū)，位于安塞縣城西北約7.5 km 處，屬延河一級(jí)支流，流域面積約3.6 km2，海拔1 100～1 327 m。流域?qū)倥瘻貛О敫珊禋夂騾^(qū)，年均溫9 ℃，年均降雨量549 mm。流域內(nèi)地形破碎，溝壑密度達(dá)6.9 km/km2；地勢(shì)起伏率和坡度大，溝峁頻繁相間，土壤侵蝕強(qiáng)烈，塑造成典型的梁狀黃土丘陵地貌和復(fù)雜多樣的土地類型；流域內(nèi)部地勢(shì)相對(duì)開闊，多平緩塌地，有大面積的梯田。流域內(nèi)絕大部分土壤是黃綿土，還有少量二色土和紅膠土[2]。

1.2 數(shù)據(jù)源

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

1.3 提取過程

ArcGIS9.3 中沒有可視化的水文分析模塊，其相關(guān)的功能要通過在Grid 模塊中調(diào)用一系列函數(shù)來實(shí)現(xiàn)。圖1的流程圖簡(jiǎn)要地說明了分析的過程。

圖1 ArcGIS進(jìn)行水文分析流程圖

1.3.1 數(shù)據(jù)預(yù)處理

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

在ArcGIS 中對(duì)洼地的填平處理包括以下步驟[6]：首先，基于原始DEM 數(shù)據(jù)生成流域的流向數(shù)據(jù)；然后，結(jié)合流向數(shù)據(jù)調(diào)用SINK 函數(shù)對(duì)洼地進(jìn)行標(biāo)定。SINK函數(shù)的功能就是標(biāo)定出流域內(nèi)高程低于周圍的封閉區(qū)域。語句格式為SINK（＜flowdir＞），flowdir 為流向數(shù)據(jù)。

標(biāo)定出洼地之后要確定洼地的影響區(qū)域，即有多少柵格的匯水到達(dá)這些洼地中。本步驟的實(shí)現(xiàn)調(diào)用Watershed函數(shù)，該函數(shù)需要流向數(shù)據(jù)和上一步中的洼地作為輸入。語句格式為watershed（＜flowdir＞，＜sinks＞），sinks為標(biāo)定的洼地柵格。

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

完成上一步之后要確定需要填充洼地的最深值，即超過這一深度的洼地被認(rèn)為是合理的地形特征，而不再進(jìn)行填充，在此深度范圍內(nèi)的洼地才填充，該值用洼地最高高程?hào)鸥駵p去最低高程?hào)鸥窦吹谩?/p>

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

因?yàn)樵谔钔莸耐瑫r(shí)改變了原有的柵格高程，可能在洼地的邊緣處產(chǎn)生新的洼地，所以以上填洼過程要反復(fù)進(jìn)行，直到用SINK 函數(shù)標(biāo)定不出洼地時(shí)才可以認(rèn)為填洼完成，生成了無洼地DEM數(shù)據(jù)。

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

1.3.2 流向數(shù)據(jù)的生成

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

在ArcGIS 中生成流向數(shù)據(jù)用的是D8 算法[7]，該方法假設(shè)單個(gè)網(wǎng)格中的水流只有8 種可能的流向，用最陡坡度法來確定水流的方向，每個(gè)柵格單元水流的流向?yàn)槠溧徑? 個(gè)柵格單元中坡度最大的那個(gè)單元格，坡度的計(jì)算公式為θ=arctan[(hi-hj)/D]，其中hi和hj為2 個(gè)單元格的高程值，D為2 個(gè)單元格中心的距離，當(dāng)相鄰柵格為共邊時(shí)，D為柵格像元的尺寸，當(dāng)相鄰柵格為對(duì)角時(shí)，D為 2倍的柵格尺寸。

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

圖2 柵格流向示意圖

1.3.3 累積匯水面積數(shù)據(jù)生成

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

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

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

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

1.3.5 溝道的分級(jí)

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

2 結(jié)果與分析

在本次研究中自動(dòng)統(tǒng)計(jì)生成的關(guān)鍵的數(shù)據(jù)有3個(gè)：填洼之后生成的無洼地DEM 數(shù)據(jù)，基于無洼地DEM數(shù)據(jù)生成的流向數(shù)據(jù)和流域累計(jì)匯水面積數(shù)據(jù)。

2.1 無洼地DEM

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

圖3 流域無洼地DEM圖

2.2 累積匯水面積分析

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

圖4 流域累積匯水面積圖

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

2.3 溝道分級(jí)結(jié)果

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

圖6 溝道分級(jí)圖

3 結(jié)論

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

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

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

2）對(duì)平地的處理效果不明顯，在地形起伏小的地區(qū)，溝道的模擬會(huì)出現(xiàn)誤差，尤其是在填平處理后產(chǎn)生的平坦地區(qū)，流向的確定仍不理想。

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

[1] 李麗，郝振純.基于DEM的流域特征提取綜述[J].地球科學(xué)進(jìn)展，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] 許明祥，劉國彬，溫仲明，等.黃土丘陵區(qū)小流域土壤特性時(shí)空動(dòng)態(tài)變化研究[J].水土保持通報(bào)，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 環(huán)境下基于DEM 的流域特征提取[J].水資源與水工程學(xué)報(bào)，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的地貌實(shí)體單元自動(dòng)提取方法[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].水利水電技術(shù)，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 環(huán)境下DEM流域特征提取及應(yīng)用[J].北京林業(yè)大學(xué)學(xué)報(bào)，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] 李翀，楊大文.基于柵格數(shù)字高程模型DEM的河網(wǎng)提取及實(shí)現(xiàn)[J].中國水利水電科學(xué)研究院學(xué)報(bào)，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] 易紅偉，湯國安，劉詠梅，等. 河網(wǎng)徑流節(jié)點(diǎn)及其基于DEM 的自動(dòng)提取[J].水土保持學(xué)報(bào)，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 的黃土高原溝壑區(qū)溝道系統(tǒng)的自動(dòng)提取[J]. 西北林學(xué)院學(xué)報(bào)，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）