基于數(shù)字高程模型的混合流向算法

夏譽(yù)玲，李小娟，王 濤

1. 首都師范大學(xué)北京市成像技術(shù)高精尖創(chuàng)新中心，北京 100094； 2. 首都師范大學(xué)地球空間信息科學(xué)與技術(shù)國際化示范學(xué)院，北京 100048； 3. 首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京100048； 4. 首都師范大學(xué)三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048； 5. 首都師范大學(xué)水資源安全北京實(shí)驗(yàn)室，北京 100048

從數(shù)字高程模型(digital elevation model，DEM)中自動(dòng)地提取匯水網(wǎng)絡(luò)和流域范圍可以為分布式水文模型及相關(guān)應(yīng)用提供重要的基礎(chǔ)信息。在現(xiàn)有提取方法中，基于模擬地表徑流累積過程的算法得到了廣泛應(yīng)用，該類算法中，地面各點(diǎn)徑流方向和流量分配對(duì)提取結(jié)果有著直接影響。目前的流向確定算法大致可以分為單流向算法、多流向算法和混合流向算法[1]。

單流向算法[2]是基于模擬地表徑流累積過程提取匯水信息的主要算法之一，其基本思想是給定地面點(diǎn)上的積水將流向較矮的最陡坡向相鄰地面點(diǎn)，在格網(wǎng)DEM中，給定柵格單元周圍有8個(gè)相鄰柵格單元，該算法又稱為D8算法。為了模擬地表流水受地物影響，文獻(xiàn)[3]在地形計(jì)算中引入隨機(jī)量，提出了隨機(jī)八方向算法。D8算法在8個(gè)離散的方向選擇流向，隨著徑流的延伸，得到的水流路徑與真實(shí)徑流的偏差會(huì)越來越大，為克服該局限性，基于全局搜索策略的D8算法[4]和改進(jìn)后的GD8算法[5]在D8的基礎(chǔ)上擴(kuò)展了檢索范圍，并通過定義第2流向以及滿足使用第2流向的一系列條件修正了逐漸偏離的水流路徑。文獻(xiàn)[6]用柵格DEM構(gòu)建規(guī)則三角網(wǎng)，根據(jù)三角面的坡向決定流向。單流向算法適于模擬水流的匯流機(jī)制，在溝谷地區(qū)能夠提取出效果良好的河網(wǎng)，然而因其流向單一，無法表現(xiàn)自然環(huán)境下的坡面漫流，在平緩區(qū)域傾向于產(chǎn)生偽平行流徑，在后續(xù)水文分析中引入了不確定性因素。

多流向算法適于模擬平緩地形中地表漫流過程。文獻(xiàn)[7—8]對(duì)D8和Rho8進(jìn)行了改進(jìn)，分別提出了水流分配權(quán)重為1的FD8和權(quán)重值為1.1的FRho8，建立漫散流動(dòng)模型，使得水可以流向中心柵格周圍的所有低高程方向柵格。考慮到不同地形條件下水流分配的不同，文獻(xiàn)[9]按FRho8方法計(jì)算匯流面積，確定匯流面積閾值，再將水流分配權(quán)重值設(shè)定為隨匯流面積變化的函數(shù)，即隨著累積匯流面積的增加該算法模擬的水流情況越來越類似于單流向算法。文獻(xiàn)[10]提出了基于局域形態(tài)確定流向的算法，該方法分析八鄰域高程值所反映的地表凸凹幾何形態(tài)，區(qū)分鄰域中所有可能被分配水流的獨(dú)立單元，依據(jù)該單元涉及的像素?cái)?shù)目決定其分配到的流量。文獻(xiàn)[11]認(rèn)為地形變化對(duì)水流分配有著重要的影響，設(shè)計(jì)了水流分配策略隨著下坡坡度變化的多流向算法。

在地表坡面徑流過程中，坡度較大的區(qū)域徑流通常呈匯聚形式，通常表現(xiàn)為單方向流動(dòng)，因此單流向算法能夠較好地模擬這些區(qū)域的水流運(yùn)動(dòng)；然而，在坡度平緩地表上模擬徑流過程時(shí)，該算法的結(jié)果中會(huì)出現(xiàn)不合理的平行流現(xiàn)象[2,12]。在相對(duì)平坦的坡面區(qū)域水流往往會(huì)分流到多個(gè)方向，這種情況下多流向算法相比于單流向算法能準(zhǔn)確地反映出水流向周圍較低高程方向分配的思想。于是出現(xiàn)了一種混合的流向算法，該方法根據(jù)地形地貌特征決定適用的流向算法，這樣綜合了多種算法的優(yōu)點(diǎn)，避免了無差別的應(yīng)用單一算法帶來的部分問題，能夠在不同的地形條件下模擬得到更加合理的水流運(yùn)動(dòng)以及水量分配。文獻(xiàn)[13]提出了在同一個(gè)DEM同時(shí)采用不同的單流向算法的方法，該算法用3×3窗口依次對(duì)DEM進(jìn)行曲面擬合，根據(jù)該曲面計(jì)算得到的正切曲率和流線曲率給中心柵格指定流向算法。文獻(xiàn)[14]將研究區(qū)分為正地形和負(fù)地形，對(duì)正地形采用多流向算法，負(fù)地形采用D8算法。本文設(shè)計(jì)實(shí)現(xiàn)了一種混合的流向算法，該方法首先對(duì)地形進(jìn)行分類，然后根據(jù)地形特征決定適用的流向算法。

最初的地形地貌提取是依靠人工判讀勾畫，隨著數(shù)據(jù)的豐富和數(shù)字地形分析技術(shù)和方法的發(fā)展，出現(xiàn)了許多地貌形態(tài)自動(dòng)分類方法?，F(xiàn)有的地形地貌分類方法有：

(1) 基于地表局部幾何形態(tài)的方法，該方法的原理是分析局部地表相對(duì)起伏度等屬性來劃分地貌類型。如文獻(xiàn)[15]使用了3×3大小的模板掃描DEM，將高程值大于南北或東西向的中心柵格標(biāo)記為山谷。文獻(xiàn)[16]提出了對(duì)高程值的下四分位數(shù)柵格區(qū)域采用自定義直徑的圓形模板進(jìn)行掃描檢測(cè)山谷線的方法。文獻(xiàn)[17]提出的Geomorphons方法，利用地形開放度選取地形特征點(diǎn)來決定與中心柵格高程的相對(duì)關(guān)系，該算法的統(tǒng)計(jì)單元搜索半徑不同，具有一定的局域地形條件自適應(yīng)能力。

(2) 對(duì)坡度、曲率等地形因子集合進(jìn)行聚類的方法[18-20]，即對(duì)DEM所派生的地形屬性集合進(jìn)行聚類，之后將聚出的類識(shí)別為地形元素類型。

(3) 基于地表整體結(jié)構(gòu)的方法。如文獻(xiàn)[21]根據(jù)地形表面水流方向提出了基于形態(tài)學(xué)的溝沿線提取方法。文獻(xiàn)[22]通過對(duì)等高線進(jìn)行凹凸段劃分，獲得等高線特征段，對(duì)特征段構(gòu)建三角網(wǎng)，從而獲取特征點(diǎn)，匹配特征段和特征點(diǎn)，完成地形特征線的追蹤，生成地形特征線。在山谷或山脊這類陡峭區(qū)域，徑流通常指向一個(gè)方向，由于高程、坡度、剖面曲率等因素，坡面地區(qū)容易產(chǎn)生不同程度的漫流。為了簡(jiǎn)單直接地對(duì)單流向和多流向算法進(jìn)行分配，只需要提取山脊線和山谷線，與坡面地區(qū)區(qū)分開來。本文采用了基于地表局部幾何形態(tài)的分類方法，利用3×3窗口內(nèi)高程的相對(duì)差異信息分別提取山谷線和山脊線[15]，對(duì)研究區(qū)的地形進(jìn)行分類。

1 混合流向算法

本文提出的混合流向算法實(shí)現(xiàn)流程如圖1所示。基于地表匯水過程模擬的算法實(shí)現(xiàn)過程中，假定地表每一點(diǎn)均有下游點(diǎn)，因此需要對(duì)DEM進(jìn)行洼地填平、平地賦流向等預(yù)處理[24-25]。

圖1 混合流向算法工作流程Fig.1 Hybrid flow direction algorithm

1.1 地形分類

本文采用了基于局部地表幾何形態(tài)檢測(cè)的方法進(jìn)行地形分類。在3×3模板窗口中，如果南北或者東西相鄰柵格的高程值大于中心柵格，則該中心柵格被標(biāo)記為山谷點(diǎn)[15]；如果南北或者東西相鄰柵格的高程值小于中心柵格，則該中心柵格被標(biāo)記為山脊點(diǎn)；如果同時(shí)滿足以上兩個(gè)條件，則標(biāo)記為鞍部點(diǎn)。當(dāng)模板完成掃描時(shí)，這些被標(biāo)記的柵格可以連接為山脊線和山谷線。由于消除了平坦區(qū)域，經(jīng)過預(yù)處理的DEM中每個(gè)柵格單元都有一個(gè)坡度值，在模板窗口檢測(cè)后的DEM中，除了山谷線，山脊線以外的區(qū)域，其余區(qū)域定義為坡面地區(qū)，地形坡度較大區(qū)域通常表現(xiàn)為徑流匯聚特性，因此利用坡度閾值將坡面區(qū)域分為陡坡區(qū)和緩坡區(qū)。經(jīng)過計(jì)算得到的用于混合流向算法中流向分配的地形圖包含了山脊線、山谷線、鞍部、陡坡和緩坡5類。

1.2 流向算法分配

現(xiàn)有的單流向算法中應(yīng)用最多的是D8，該算法為每個(gè)柵格單元指定其較低的鄰居?xùn)鸥駟卧凶疃钙孪驗(yàn)樗鞣较?。D8在模擬山谷、山脊等陡峭區(qū)域的徑流過程中能夠獲得較為理想的效果，且相較于其他單流向算法，其算法更加簡(jiǎn)單，易于實(shí)現(xiàn)，計(jì)算效率高。這里沿用以往做法，采用2的整數(shù)次冪來標(biāo)記流向，8個(gè)流向分別使用1、2、4、8、16、32、64、128，從當(dāng)前柵格單元左側(cè)柵格開始以順時(shí)針方向的不同流向進(jìn)行標(biāo)識(shí)，如圖2(a)所示。

在緩坡以及鞍部區(qū)域，采用MFD-fg算法建立徑流模型。MFD-fg以最大下坡坡度的線性函數(shù)對(duì)水流分配進(jìn)行加權(quán)，水流分配策略隨著地形變化而變化，使得陡峭的方向分配的水流量越多，平緩的方向分配的水流量越少。通過這種自適應(yīng)的水流分配算法計(jì)算得到的匯水面積(本文用匯水區(qū)柵格數(shù)代表匯水面積)與D8算法得到的值銜接相對(duì)自然，因此選用MFD-fg算法模擬平緩地區(qū)的水流運(yùn)動(dòng)。MFD-fg的水流分配策略由式(1)定義

(1)

式中，di為第i個(gè)鄰域柵格的水流分配比例；tanβi為坡度；f(e)(f(e)>0)為流量分配權(quán)重，其中e代表的是最大坡度值；emin、emax分別表示的3×3窗口中e的最大值和最小值[11]。Li為第i個(gè)相鄰柵格的等高線長(zhǎng)度加權(quán)因子，如式(2)所示

(2)

混合流向算法采用與單流向算法一致的流向標(biāo)記方式。如圖2(b)所示：具有多個(gè)流向的柵格，其流向值由代表各個(gè)方向的數(shù)值之和表示，如：流向值為3的柵格其水流方向指向東和東南，流向值為255的柵格其水流方向指向周圍8個(gè)柵格。

圖2 混合流向值示意圖Fig.2 Schematic plot of hybrid flow direction value

2 試驗(yàn)與結(jié)果

本文采用了2000年“航天飛機(jī)雷達(dá)地形測(cè)量計(jì)劃”獲取的DEM數(shù)據(jù)，分別在山西省永和縣的黃土地貌區(qū)和重慶市巫溪縣的中山切割地貌區(qū)各選擇了一個(gè)試驗(yàn)樣區(qū)，二者均不位于城區(qū)，研究區(qū)地貌特征如圖3所示，圖3(a)樣區(qū)面積為14.9 km2，平均高程為1 109.1 m，平均坡度為18.6°。圖3(b)樣區(qū)面積為43.6 km2，平均高程為1 399.9 m，平均坡度為32.6°。使用的DEM空間分辨率約3″和1″。

圖3 研究區(qū)概況Fig.3 The general situation of study area

首先對(duì)原始SRTM數(shù)據(jù)進(jìn)行投影轉(zhuǎn)換，1″和3″數(shù)據(jù)分別被轉(zhuǎn)換為高斯投影下30 m和90 m的DEM?？紤]到SRTM高程數(shù)據(jù)為地表高程模型(digital surface model)，同時(shí)在試驗(yàn)區(qū)內(nèi)存在低矮建筑、植被等信息，為了減少噪聲影響，進(jìn)行了高斯濾波處理。DEM中的洼地和平坦地區(qū)會(huì)導(dǎo)致柵格沒有流向或者流向不明確，因此，進(jìn)行了洼地填充處理。

接著，利用幾何形態(tài)法提取了山谷線和山脊線。根據(jù)自然斷點(diǎn)分級(jí)法，確定坡度閾值。永和樣區(qū)30 m分辨率的DEM的坡度閾值為17.2°，90 m分辨率DEM的坡度閾值為11.6°；巫溪樣區(qū)30 m分辨率的DEM的坡度閾值為30.7°，90 m分辨率DEM的坡度閾值為27.2°。然后，通過特征線和坡度閾值將流域內(nèi)的地形劃分為：山脊、山谷、鞍部、陡坡和緩坡5類，對(duì)山谷、山脊以及陡坡區(qū)域采用單流向的D8算法，對(duì)鞍部和緩坡地區(qū)采用多流向的MFD-fg算法。如圖4所示，在地形分類圖中，白色柵格表示山脊，黑色柵格表示山谷，灰色系表示坡地和鞍部；在流向算法分配柵格中，灰色代表單流向算法區(qū)，黑色代表多流向算法區(qū)。

圖4 分類結(jié)果Fig.4 The result of classification

分別利用單流向算法D8、Rho8，多流向算法FD8、MFD-fg，以及混合流向算法對(duì)兩個(gè)試驗(yàn)區(qū)計(jì)算了地表各點(diǎn)匯水面積。如圖5和圖6所示，D8算法和Rho8算法得到的集水區(qū)面積在緩坡地區(qū)生成了大量與實(shí)際不符的單向徑流和平行直線。FD8算法和MFD-fg算法得到的結(jié)果在山脊、山谷和陡坡地區(qū)產(chǎn)生了不合理的漫散分流的水流運(yùn)動(dòng)，MFD-fg算法雖在水量分配上有改進(jìn)，但依然沒有克服此類水流運(yùn)動(dòng)在陡坡區(qū)域的不合理現(xiàn)象?；旌狭飨蛩惴ńY(jié)合了單流向算法能較好模擬陡峭區(qū)域水流運(yùn)動(dòng)的特點(diǎn)，以及多流向算法在平緩坡區(qū)模擬水流漫散流動(dòng)的優(yōu)勢(shì)，相比于另外4種算法得到的匯水面積柵格圖，可以看到陡峭區(qū)域水流匯聚，徑流寬度相對(duì)變窄；平緩坡區(qū)的徑流分流，根據(jù)地形特點(diǎn)水量合理分配到中心柵格周圍的低高程?hào)鸥??；旌狭飨蛩惴ǜ雍侠淼胤从沉说匦巫兓瘜?duì)水流運(yùn)動(dòng)以及水量分配的影響。兩個(gè)樣區(qū)30 m分辨率的匯水面積結(jié)果相比于90 m的結(jié)果都更加詳細(xì)地反映了復(fù)雜地形的水量分配。

圖5 永和樣區(qū)5種流向算法的匯水面積圖Fig.5 Catchment area results of five algorithms in Yonghe area

圖6 巫溪樣區(qū)5種流向算法的匯水面積結(jié)果Fig.6 Catchment area results of five algorithms in Wuxi area

3 分析與討論

為了進(jìn)一步驗(yàn)證混合流向算法在傳統(tǒng)算法上的改進(jìn)，分別采用了匯水面積差異圖、累積頻率圖、散點(diǎn)分布圖矩陣來定量化分析和描述混合流向算法與傳統(tǒng)算法的相似性和差異性。

3.1 匯水面積差異圖

用混合流向算法分別與另外4種算法得到的匯水面積進(jìn)行相減計(jì)算，將得到的匯水面積差值結(jié)果進(jìn)一步分層設(shè)色表示(見圖7、圖8)。

圖7 永和樣區(qū)匯水面積差異Fig.7 The difference of catchment area in Yonghe area

圖8 巫溪樣區(qū)匯水面積差異Fig.8 The difference of catchment area in Wuxi area

因?yàn)閰R流區(qū)域的上游柵格數(shù)量大，累積了所有上游柵格匯水面積的變化量，所以得到的變化較大的地方主要集中在溝谷區(qū)域。4種傳統(tǒng)算法與混合流向算法相減的差異從小到大依次排序是：D8、MFD-fg、Rho8、FD8。對(duì)比不同分辨率DEM，30 m分辨率DEM計(jì)算的匯水面積差值明顯大于90 m分辨率DEM，因?yàn)楦叻直鍰EM的柵格數(shù)多，相應(yīng)位置匯流區(qū)域得到的匯水面積值大于低分辨率DEM；并且高分辨率DEM中的地形分類更加復(fù)雜，混合流向算法相應(yīng)的流向分配策略更加詳細(xì)，因此在高分辨率柵格中，混合流向算法與其他算法的差異體現(xiàn)得更加明顯。

3.2 累積頻率

由于匯水面積數(shù)值區(qū)間較大，為了方便對(duì)比，首先對(duì)其進(jìn)行對(duì)數(shù)變換，在此基礎(chǔ)上統(tǒng)計(jì)累積頻率[23]。圖9為兩個(gè)樣區(qū)不同分辨率DEM計(jì)算所得的匯流面積累積頻率分布圖。累積頻率分布圖在ln(catchment area)≤6.4時(shí)趨于一致，約占總數(shù)的98%。圖中單流向算法和多流向算法可以明顯區(qū)分開來，混合流向算法由于在陡峭地區(qū)采用了單流向算法，平緩區(qū)域采用了多流向算法，從累積頻率圖上可以看出，D8和Rho8中匯水面積較小的柵格占比大于多流向柵格，其原因是單流向算法的下游區(qū)域完全接受上游匯入，而多流向算法中較高區(qū)域徑流向多個(gè)方向擴(kuò)散?；旌狭飨蚶鄯e曲線介于D8和MFD-fg之間，說明匯水面積分布介于單流向算法和多流向算法結(jié)果之間，能較好地兼顧不同區(qū)域水流匯聚和坡面散流的特性。

3.3 散點(diǎn)圖

圖10是在30 m和90 m DEM中分別使用4種傳統(tǒng)流徑算法和混合流向算法所得到的匯水面積的散點(diǎn)圖。

在同樣分辨率DEM上，混合流向算法的匯水面積值與D8和MFD-fg的分布接近，相關(guān)系數(shù)最高，其次是FD8和Rho8。4種算法與混合流向算法在較大值處都有部分一致分布。對(duì)比不同分辨率DEM，30 m分辨率的DEM上混合流向算法與其他4種算法的相似性低于90 m分辨率的DEM。因?yàn)镈EM的分辨率越高，所能刻畫的地形細(xì)節(jié)越豐富，混合流向算法的流向分配方案適用性更強(qiáng)，所以得到結(jié)果與其他算法的差異比低分辨率DEM得到的結(jié)果差異更大，與累積頻率分析一致。也就是說，在更高分辨率的DEM上，混合流向算法的改進(jìn)效果更加明顯。對(duì)比兩個(gè)地貌試驗(yàn)區(qū)的結(jié)果，中山切割地貌水流匯聚特征得到了反映，抑制了多流向算法結(jié)果的散流現(xiàn)象，而在黃土地貌坡度平緩區(qū)域單流向造成的漫流特征也得到了加強(qiáng)。

圖9 ln(catchment area)累積頻率分布圖Fig.9 Cumulative frequency distribution of ln(catchment area)

圖10 ln(catchment area)散點(diǎn)分布圖Fig.10 Scatter plot matrix for ln(catchment area)

4 結(jié) 論

流向算法直接影響了匯水面積的計(jì)算結(jié)果，是匯水網(wǎng)絡(luò)提取以及匯水單元?jiǎng)澐值年P(guān)鍵步驟，對(duì)研究分布式水文模型具有重要的意義。多年來，研究者們提出了各種模擬水流運(yùn)動(dòng)以及流量分配的方法。單流向算法雖然擅長(zhǎng)模擬溝谷等地形中的水流的匯流機(jī)制，但在模擬坡面水流時(shí)會(huì)受到單流向的局限，無法表現(xiàn)漫流情景，導(dǎo)致結(jié)果中存在大量平行流現(xiàn)象。多流向算法中的水能流向所有高程較低的方向，彌補(bǔ)了單流向算法模擬坡面漫流的缺陷，但也失去了單流向算法模擬陡峭地形水流運(yùn)動(dòng)優(yōu)勢(shì)。

本文研究提出的混合流向算法實(shí)現(xiàn)了單流向算法和多流向算法的結(jié)合，充分考慮了各種地形變化對(duì)水流運(yùn)動(dòng)以及水量分配的影響，對(duì)比結(jié)果顯示，該算法得到的匯水面積更加合理，并且在相對(duì)較高分辨率的DEM中能夠得到更加明顯的改進(jìn)。

地形類別的劃分對(duì)于流向算法的分配有著至關(guān)重要的影響，在后續(xù)的工作中，將對(duì)地形分類和坡度閾值的判斷選擇作進(jìn)一步研究，以便改進(jìn)算法。

參考文獻(xiàn)：

[1] 盧慶輝, 熊禮陽, 蔣如喬, 等. 一種融合Priority-Flood算法與D8算法特點(diǎn)的河網(wǎng)提取方法[J]. 地理與地理信息科學(xué), 2017, 33(4): 40-46.

LU Qinghui, XIONG Liyang, JIANG Ruqiao, et al. A Drainage Networks Extraction Algorithm by Integrating the Characteristics of Priority-flood and D8 Algorithm[J]. Geography and Geo-Information Science, 2017, 33(4): 40-46.

[2] 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(3): 323-344.

[3] FAIRFIELD J, LEYMARIE P. Drainage Networks from Grid Digital Elevation Models[J]. Water Resources Research, 1991, 27(5): 709-717.

[4] PAIK K. Global Search Algorithm for Nondispersive Flow Path Extraction[J]. Journal of Geophysical Research: Earth Surface, 2008, 113(F4): F04001.

[5] SHIN S, PAIK K. An Improved Method for Single Flow Direction Calculation in Grid Digital Elevation Models[J]. Hydrological Processes, 2017, 31(8): 1650-1661.

[6] ZHOU Qiming, PILESJ? P, CHEN Yumin. Estimating Surface Flow Paths on a Digital Elevation Model Using a Triangular Facet Network[J]. Water Resources Research, 2011, 47(7): W07522.

[7] FREEMAN T G. Calculating Catchment Area with Divergent Flow Based on a Regular Grid[J]. Computers & Geosciences, 1991, 17(3): 413-422.

[8] QUINN P, BEVEN K, CHEVALLIER P, et al. The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling Using Digital Terrain Models[J]. Hydrological Processes, 1991, 5(1): 59-79.

[9] QUINN P F, BEVEN K J, LAMB R. The in(a/tan/β) Index: How to Calculate It and How to Use It within the Topmodel Framework[J]. Hydrological Processes, 1995, 9(2): 161-182.

[10] PILESJ? P, ZHOU Qiming, HARRIE L. Estimating Flow Distribution over Digital Elevation Models Using a Form-based Algorithm[J]. Geographic Information Sciences: A Journal of the Association of Chinese Professionals in Geographic Information Systems, 1998, 4(1-2): 44-51.

[11] 秦承志, 李寶林, 朱阿興, 等. 水流分配策略隨下坡坡度變化的多流向算法[J]. 水科學(xué)進(jìn)展, 2006, 17(4): 450-456.

QIN Chengzhi, LI Baolin, ZHU Axing, et al. Multiple Flow Direction Algorithm with Flow Partition Scheme Based on Downslope Gradient[J]. Advances in Water Science, 2006, 17(4): 450-456.

[12] GALLANT J C, WILSON J P. Primary Topographic Attributes[M]∥WILSON J P, GALLANT J C. Terrain Analysis: Principles and Applications. New York: John Wiley & Sons, 2000: 51-85.

[13] 張占陽, 鄭紅梅, 郭新成, 等. 局部曲面正切曲率與流線曲率的水流路徑算法[J]. 測(cè)繪科學(xué), 2013, 38(2): 122-124.

ZHANG Zhanyang, ZHENG Hongmei, GUO Xincheng, et al. Flow Paths Algorithm Based on Tangential Curvature and Flow-path Curvature of Partial Surface[J]. Science of Surveying and Mapping, 2013, 38(2): 122-124.

[14] XIONG Liyang, TANG Guoan, YAN Shijiang, et al. Landform-oriented Flow-routing Algorithm for the Dual-structure Loess Terrain Based on Digital Elevation Models[J]. Hydrological Processes, 2014, 28(4): 1756-1766.

[15] JOHNSTON E G, ROSENFELD A. Digital Detection of Pits, Peaks, Ridges, and Ravines[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1975, SMC-5(4): 472-480.

[16] LINDSAY J B. Sensitivity of Channel Mapping Techniques to Uncertainty in Digital Elevation Data[J]. International Journal of Geographical Information Science, 2006, 20(6): 669-692.

[17] JASIEWICZ J, STEPINSKI T F. Geomorphons —A Pattern Recognition Approach to Classification and Mapping of Landforms[J]. Geomorphology, 2013, 182: 147-156.

[18] BURROUGH P A, VAN GAANS P F M, MACMILLAN R A. High-resolution Landform Classification Using Fuzzy K-means[J]. Fuzzy Sets and Systems, 2000, 113(1): 37-52.

[21] 閭國年, 錢亞東, 陳鐘明. 基于柵格數(shù)字高程模型自動(dòng)提取黃土地貌溝沿線技術(shù)研究[J]. 地理科學(xué), 1998, 18(6): 567-573.

LU Guonian, QIAN Yadong, CHEN Zhongming. Study of Automated Extraction of Shoulder Line of Valley from Grid Digital Elevation Data[J]. Scientia Geographica Sinica, 1998, 18(6): 567-573.

[22] 張堯, 樊紅, 李玉娥. 一種基于等高線的地形特征線提取方法[J]. 測(cè)繪學(xué)報(bào), 2013, 42(4): 574-580.

ZHANG Yao, FAN Hong, LI Yu’e. A Method of Terrain Feature Extraction Based on Contour[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(4): 574-580.

[23] 劉學(xué)軍, 晉蓓, 王彥芳. DEM流徑算法的相似性分析[J]. 地理研究, 2008, 27(6): 1347-1357.

LIU Xuejun, JIN Bei, WANG Yanfang. Similarity Analysis of Flow Route Algorithms for Extracting Drainage Network from Grid-based Terrain Model[J]. Geographical Research, 2008, 27(6): 1347-1357.

[24] JENSON S K, DOMINGUE J O. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis[J]. Photogrammetric Engineering and Remote Sensing, 1988, 54(11): 1593-1600.

[25] 李輝, 陳曉玲, 張利華, 等. 基于三方向搜索的DEM中洼地處理方法[J]. 水科學(xué)進(jìn)展, 2009, 20(4): 473-479.

LI Hui, CHEN Xiaoling, ZHANG Lihua, et al. Depression Removal Method for Grid DEM Based on Three-direction Search[J]. Advances in Water Science, 2009, 20(4): 473-479.