基于DEM不同路徑算法的溝壑密度提取

郭蘭勤，丑述仁

(1.西北農(nóng)林科技大學(xué) 資源環(huán)境學(xué)院，陜西 楊凌712100;2.西北大學(xué) 城市與環(huán)境學(xué)院，陜西，西安，710127)

溝壑是黃土侵蝕地貌中很重要的一種地貌類型，其形成過(guò)程是流域地貌的演化過(guò)程，也是土壤侵蝕的過(guò)程［1］。溝壑密度反映了一個(gè)地區(qū)被水流切割破碎的程度，同時(shí)溝壑密度也是地形發(fā)育階段、降水量或地勢(shì)高差、土壤滲透能力和地表抗侵蝕能力的重要特征值［2］。溝壑密度的調(diào)查方法有傳統(tǒng)的外業(yè)調(diào)查法、地形圖量算法和基于DEM提取法等。外業(yè)調(diào)查工作量巨大，適合范圍狹小，不能滿足大范圍水土流失程度調(diào)查、評(píng)價(jià)及預(yù)測(cè)工作。張麗萍建立了溝壑密度隨切割深度變化的理論極值模型，并進(jìn)行了溝壑密度的理論極值建模。蔣忠信對(duì)張麗萍等建立的流域溝壑密度理論極值圓形數(shù)學(xué)模式進(jìn)行了更正。由于改正的模式又與實(shí)際相矛盾，又推導(dǎo)出更符合實(shí)際的正六邊形教學(xué)模式［3］。韋中亞用數(shù)學(xué)形態(tài)學(xué)的溝整密度提取算法，從TM影像中提取溝谷線信息，由于該算法的計(jì)算結(jié)果隨窗口的大小變化較大，未形成標(biāo)準(zhǔn)［4－5］。隨著數(shù)字地形分析技術(shù)和地理信息系統(tǒng)技術(shù)的發(fā)展，基于數(shù)字高程模型(DEM)提取溝壑密度已經(jīng)得到越來(lái)越廣泛的使用，但是傳統(tǒng)的通過(guò)DEM提取溝壑密度都是基于D8的單流向算法，且匯流閾值的選取主觀性太大，沒(méi)有具體的規(guī)范［6－9］。D8算法計(jì)算結(jié)果在平緩地區(qū)會(huì)產(chǎn)生平行水流，不能模擬分流，且會(huì)把2維流路簡(jiǎn)化為1維，且匯流閾值的確定具有一定的地貌局限性，錯(cuò)誤的使用其他地貌類型的匯流閾值，將產(chǎn)生大量的偽水道，國(guó)內(nèi)學(xué)者孫友波，熊立華等對(duì)匯流閾值的確定進(jìn)行了大量的研究，但都沒(méi)有得到非常滿意的結(jié)果［10－17］。而利用多流向算法(DEMON)能給出較可靠的匯水面積值和較理想的匯水面積分布。利用多流向算法、和比較合理的匯流閾值設(shè)定方法來(lái)計(jì)算溝壑密度的方法在國(guó)內(nèi)還比較少見(jiàn)。

本文選擇陜北黃土高原縣南溝流域?yàn)檠芯繀^(qū)，分別利用坡向流量分配的多流向算法(DEMON)和D8的單流向算法計(jì)算該地區(qū)的溝壑密度，并對(duì)這兩種方法計(jì)算的結(jié)果進(jìn)行比較分析，得出較合理的溝壑密度提取方法，對(duì)以后的黃土高原溝壑密度空間分異的研究，以及水土保持規(guī)劃具有重大的意義。

1 研究區(qū)域概況

縣南溝流域位于陜北黃土高原中部，屬延河流域子流域，該流域位于東經(jīng)109°11'15″～109°22'30″，北緯 36°42'30″～36°47'30″，平均海拔為1 220 m，屬于典型的黃土丘陵溝壑地貌類型，其地形復(fù)雜，地表破碎，溝壑縱橫，地面坡度陡峭，平均坡度達(dá)28°，是典型的受人類活動(dòng)影響的水土流失嚴(yán)重區(qū)。

2 研究方法

2.1 實(shí)驗(yàn)數(shù)據(jù)和軟件

實(shí)驗(yàn)數(shù)據(jù)由國(guó)家基礎(chǔ)地理數(shù)據(jù)庫(kù)提供的1:10 000比例尺DLG數(shù)據(jù)線劃圖生成的分辨率為5 m的Hc－DEM。軟件采用ESRI公司開發(fā)的地理信息系統(tǒng)平臺(tái)ArcGIS9.2以及SagaGIS。

2.2 實(shí)驗(yàn)原理與過(guò)程

溝壑密度(gully density)也叫溝道密度或溝谷密度，指單位面積內(nèi)溝壑的總長(zhǎng)度，單位一般以km/km2表示，數(shù)學(xué)表達(dá)式為:

式中:D為溝壑密度，∑L為流域范圍內(nèi)溝壑的總長(zhǎng)度，一般以km為單位，A為流域面積，其單位為km2。

基于格網(wǎng)DEM的路徑算法主要關(guān)系到兩個(gè)問(wèn)題:(1)當(dāng)前柵格單元流向的確定，(2)決定當(dāng)前柵格單元向下游較低單元的流量分配比例。對(duì)應(yīng)得產(chǎn)生了兩種流向算法，單流向算法(Single Flow Direction Algorithm，SFD)和多流向算法(Multi Flow Direction Algorithm，MFD)。單流向算法認(rèn)為水全部流入8鄰域中處于最陡下坡方向上的那個(gè)像元，由于SFD不適于模擬實(shí)際表面水流的問(wèn)題，所以出現(xiàn)了MFD。多流向算法認(rèn)為水流具有分散性，中心象元的流量并不是流入最低的鄰域象元，而是流入比其低的所有鄰域象元。本文選擇單流向的D8算法、多流向的DEMON算法來(lái)確定流向。

由于原始DEM中數(shù)據(jù)噪音，內(nèi)插方法的影響，其數(shù)據(jù)中經(jīng)常包含一些洼地，洼地使流域中水流不通暢，不能形成完整的流域網(wǎng)絡(luò)，所以，必須對(duì)DEM中的洼地進(jìn)行填充處理。然后分別用多流向和單流向算法計(jì)算水流方向，進(jìn)而分別提取出流域的匯流面積，也就是流向該格網(wǎng)的所有的上游格網(wǎng)單元的水流累積量。提取溝谷線時(shí)，必須給匯流面積設(shè)定一個(gè)合理的匯流閾值，才能提取出符合實(shí)際的溝谷線，具體過(guò)程如圖1所示。

3 結(jié)果與分析

3.1 匯水面積的結(jié)果分析

匯水面積是流經(jīng)一段等高線上游的所有地形的投影面積，分別用D8算法和DEMON算法計(jì)算縣南溝流域的匯水面積，然后分別對(duì)這兩種方法計(jì)算的結(jié)果進(jìn)行分析比較，從而得出多流向算法計(jì)算得到的匯水面積分布更合理，能反映出更多地形變化對(duì)水流分配的影響細(xì)節(jié)。

D8算法計(jì)算的匯水面積可以總體上表現(xiàn)水流路徑，隨著高程的降低，水流會(huì)從匯水點(diǎn)開始一直流出流域的出口，越在下游，接受的流量越多，但是單流向算法只用下游的一個(gè)柵格便接受了全部的上游流量，流域最高點(diǎn)到出口處，只能看到匯流的基本軌跡，如果實(shí)際溝谷的寬度遠(yuǎn)遠(yuǎn)大于單位柵格尺寸，單流向算法就無(wú)能為力了。多流向算法提取的匯水面積可以很好的克服單流向算法的缺點(diǎn)，其提取的匯水面積，不僅可以很好的表現(xiàn)出水流的路徑，而且在溝道比象元尺寸大的地方更恰當(dāng)?shù)谋磉_(dá)了匯水面積的范圍，更加接近實(shí)際，同時(shí)也克服了單流向算法遇到分散水流時(shí)不能正確表達(dá)水流路徑的缺點(diǎn)。

圖1 溝壑密度的計(jì)算流程圖

3.2 匯流閾值的確定

在水文分析中，提取水系異常重要，然而水系提取中關(guān)鍵性的一步是匯流閾值的確定，在單流向算法中，對(duì)匯水面積設(shè)置合理的閾值，就可以提取大于匯流閾值的溝道，在ArcGIS中，用CON函數(shù)計(jì)算出大于匯流閾值的溝道，然后用水文分析工具下的Stream to feature，就可以把柵格的溝谷轉(zhuǎn)換為矢量的溝谷網(wǎng)絡(luò)。本研究用不同的匯流閾值計(jì)算出其相對(duì)應(yīng)的流域的溝谷長(zhǎng)度、最大值、最小值、標(biāo)準(zhǔn)差和溝壑密度，匯流閾值從100到1500，每間隔100取值一次，如表1所示。

表1 不同閾值所對(duì)應(yīng)的特征數(shù)據(jù)表

從表中可以看出，隨著閾值的增大，溝谷總長(zhǎng)度在減小，河流條數(shù)也在減少，最大值先增大后減小，平均值和標(biāo)準(zhǔn)方差都逐漸增大，匯流閾值在500時(shí)，溝谷長(zhǎng)度開始趨于穩(wěn)定，然后一直保持緩慢變小趨勢(shì)，河流條數(shù)也開始變的相對(duì)穩(wěn)定，最大值突然增大然后保持穩(wěn)定，平均值和標(biāo)準(zhǔn)差增加速度變慢，從匯流閾值和溝谷總長(zhǎng)度的坐標(biāo)圖上也可以看出(圖2)，當(dāng)匯流面積大于500時(shí)，曲線出現(xiàn)轉(zhuǎn)折，明顯的變緩。

所以當(dāng)匯流閾值為500時(shí)，各特征參數(shù)都趨于穩(wěn)定，所有坡地溝谷從流域中剔除，當(dāng)然閾值的確定只從該圖中唯一確定缺乏精確性，因?yàn)檫x取的閾值是每隔100取一次，所以本文確定的閾值應(yīng)該是一個(gè)范圍，位于400～600之間，可以看出該方法優(yōu)于主管臆斷得到的匯流閾值。并且若提取溝谷網(wǎng)絡(luò)時(shí)閾值太小，則提取的溝谷過(guò)于詳細(xì)，有很多的平行流偽溝谷;若閾值太大，則提取的溝谷網(wǎng)絡(luò)過(guò)于簡(jiǎn)單，不能反映溝谷實(shí)際情況，而本研究選取的閾值提取的溝谷網(wǎng)絡(luò)圖則較接近實(shí)際。

3.3 多流向算法和單流向算法結(jié)果分析

確定合理的匯流閾值后，則用不同路徑算法來(lái)提取溝谷網(wǎng)絡(luò)，并計(jì)算最終的溝壑密度。結(jié)果表明多流向算法在計(jì)算溝谷網(wǎng)絡(luò)和溝壑密度時(shí)，能準(zhǔn)確的表達(dá)地形和水流在坡面時(shí)的流向分散性，提取的溝谷網(wǎng)絡(luò)比較準(zhǔn)確，彌補(bǔ)了單流向算法在這方面的不足。然而正是由于多流向的這種特性，也造成了它的一些缺陷，在提取溝谷網(wǎng)絡(luò)時(shí)，會(huì)出現(xiàn)很多的偽溝谷和一些交叉溝谷，并且在同一河道出現(xiàn)幾條并列的溝谷，給進(jìn)一步計(jì)算溝谷長(zhǎng)度帶來(lái)困難，所以要對(duì)這些不正確的溝谷進(jìn)行手工修改，去除這些不符合實(shí)際的河谷。然后用DEMON多流向算法計(jì)算出該流域的溝谷總長(zhǎng)度和溝壑密度，這里假定單流向算法確定的匯流閾值是該地區(qū)溝谷的起始點(diǎn)，所以在多流向算法中也采用這一閾值，提取的溝谷如下圖3。

圖2 閾值與溝谷長(zhǎng)度的關(guān)系圖

計(jì)算出多流向DEMON算法提取的溝谷總長(zhǎng)度，用溝谷總長(zhǎng)度除以流域面積，便得到了該地區(qū)的溝壑密度，然后和單流向算法計(jì)算的結(jié)果進(jìn)行分析比較。表2表明多流向和單流向算法都可以提取到較合適的溝壑密度，都和縣南溝流域的實(shí)地調(diào)查溝壑密度4.4km/km較接近，但是多流向算法DEMON的提取精度更好一些。

圖3 FMFD和DEMON提取的溝谷網(wǎng)絡(luò)及光照陰影疊加圖

表2 不同路徑算法計(jì)算結(jié)果

4 結(jié)論與討論

(1)本研究根據(jù)匯流閾值大小與溝谷長(zhǎng)度的關(guān)系圖，確定了縣南溝流域合理的匯流閾值為400－600。利用SFD單流向算法和DEMON多流向算法都可以提取到較合適的溝壑密度，且多流向算法DEMON提取結(jié)果更接近實(shí)地調(diào)查值一些。

(2)是否多流向算法用在任何地方都能得出合理的結(jié)果，是值得思考和研究的問(wèn)題。且對(duì)不同分辨率DEM，不同的流域，不同路徑算法進(jìn)行溝壑密度的計(jì)算需進(jìn)一步的研究。溝壑密度的提取過(guò)程中，匯流閾值的設(shè)置不僅和當(dāng)?shù)氐牡刭|(zhì)地貌、水文、氣象等因素有關(guān)，而且和當(dāng)?shù)氐娜祟惢顒?dòng)也有很大的關(guān)系，故溝壑密度提取中匯流閾值的確定尚需通過(guò)溝壑密度的不確定性研究來(lái)解決。

［1］張麗萍，馬志正.流域地貌演化的不同階段溝壑密度與切割深度關(guān)系研究［J］，地理研究，1998，17(3):273 －278.

［2］周啟鳴，劉學(xué)軍.數(shù)字地形分析［M］.北京:科學(xué)出版社.2006.

［3］陳濤，朱金兆等.黃土高原溝壑區(qū)溝道提取時(shí)閾值范圍的研究［J］，長(zhǎng)江科學(xué)院院報(bào)，2008，53(3):28 －30.

［4］熊立華，郭練生.基于DEM的數(shù)字河網(wǎng)生成方法的探討［J］.長(zhǎng)江科學(xué)院院報(bào).2003，20(4):14 －17.

［5］韋中亞，周貴云，羅萬(wàn)勤.一種基于數(shù)學(xué)形態(tài)學(xué)的溝壑密度提取算法［J］，地理學(xué)與國(guó)土研究，2001，17(2):25－27.

［6］Yang，QK，Van Niel，TG，McVicar，TR，et.al.Developing a digital elevation model using ANUDEM for the Coarse Sandy Hilly Catchments of the Loess Plateau，China.CSIRO Land and Water Technical Report，7(05):1 － 74.

［7］蔣忠信.流域溝壑密度理論極值數(shù)學(xué)模式商討［J］.地理研究，1999，18(2):220 －223.

［8］孔凡哲，李莉莉.利用DEM提取河網(wǎng)時(shí)集水面積閾值的確定［J］.水電能源科學(xué).2005，23(4):65 －67.

［9］楊勤科，李銳，王占禮.區(qū)域水土流失監(jiān)測(cè)與評(píng)價(jià)指標(biāo)體系研究.水土保持通報(bào)，2000，20(7).

［10］孫友波，宮輝力.基于DEM的數(shù)字河網(wǎng)生成方法的淺議［M］.首都師范大學(xué)學(xué)報(bào)(自然科學(xué)版).2005，26(2)，106 －111.

［11］O’CALLAGHAN J F，MARK D M.The extraction of drainage networks from digital elevation data［J］.Computer Vision，Graphics and Image Processing，1984，28:323 ～344.

［12］鄔倫，汪大明，張毅.基于DEM的水流方向算法研究［J］.中國(guó)圖形圖像學(xué)報(bào)，2006，11(7):998 ～1003.

［13］TARBOTON D.A new method for the determination of flow directions and contributing areas in grid digital elevation models［J］.Water Resources Research，1997，33:309 －319.

［14］MARIZA C，COSTA － CABRAL，STEPHEN J.BURGES.Digital elevation model networks(DEMON):A model of flow over hillslopes for computation of contributing and dispersal areas［J］.Water Resource Research，1994，30(6):1681 ～1692.

［15］QUINN P，BEVEN K，CHEVALIER V，et al.The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models［J］.Hydrological Processes，1991，5(1):59 ～ 79.

［16］秦承志，朱阿興，李寶林等.基于柵格DEM的多流向算法述評(píng)［J］.地學(xué)前緣，2006，13(3):91 －98.

［17］鄭子彥，張萬(wàn)昌，邰慶國(guó).基于DEM與數(shù)字化河道提取流域河網(wǎng)的不同方案比較研究［J］.資源科學(xué)，2009，31(10):1730－1739.

［18］李俊，湯國(guó)安等.利用DEM提取陜北黃土高原溝谷網(wǎng)絡(luò)的匯流閾值研究［J］，水土保持通報(bào)，2007，27(2):75－78.

［19］黨福江，戈素芬.溝壑密度調(diào)查應(yīng)采取水文網(wǎng)絡(luò)法［J］.水土保持通報(bào)，1998，18(6):31 －33.

［20］湯國(guó)安，楊昕.ArcGIS地理信息系統(tǒng)空間分析試驗(yàn)教程［M］.北京:科學(xué)出版社，2006.