張路路，劉召芹，劉峰，李巍

（1.山東科技大學(xué)，山東青島266590；2.中國科學(xué)院遙感與數(shù)字地球研究所遙感科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京100101）

基于數(shù)字高程模型的古滑坡區(qū)提取

張路路1，劉召芹2，劉峰1，李巍2

（1.山東科技大學(xué)，山東青島266590；2.中國科學(xué)院遙感與數(shù)字地球研究所遙感科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京100101）

針對古滑坡的滑前影像無法獲得，植被、紋理信息都已恢復(fù)，無法通過對比分析滑坡滑動前后的植被、紋理等信息的變化來提取滑坡區(qū)域的問題，提出了一種新的基于數(shù)字高程模型的滑坡區(qū)域范圍提取方法。該方法基于簡化的滑坡體模型及特征分析，對滑坡區(qū)進(jìn)行水流方向、坡度、山脊山谷線提取，通過流域分析獲取滑坡區(qū)域范圍；利用坡度圖實(shí)現(xiàn)滑坡壁與滑坡體提取。實(shí)驗(yàn)利用全球30m分辨率ASTER GDEM數(shù)據(jù)，提取了四川理縣3個古滑坡體區(qū)域范圍，驗(yàn)證了該方法的有效性。

古滑坡；滑坡壁；滑坡體；數(shù)字高程模型；流域分析

0 引 言

據(jù)初步統(tǒng)計(jì)，全國至少有400多個市、縣、區(qū)、鎮(zhèn)，10000多個村莊受到滑坡災(zāi)害嚴(yán)重侵害［1］，有證可查的滑坡災(zāi)害點(diǎn)約為41萬多處，總面積為1.7352×104km2，占國土總面積的18.10%［2］（截至2000年）。上世紀(jì)八十年代以來，我國造成10人以上死亡和經(jīng)濟(jì)損失超過500萬元的滑坡就超過130起。這些大型災(zāi)難性滑坡的發(fā)生，不僅帶來重大的人員傷亡和財(cái)產(chǎn)損失，而且引發(fā)了嚴(yán)重的社會和公共安全問題。除此之外，古滑坡的潛在危害也是不容忽視的，古滑坡在強(qiáng)烈的外部作用下能夠重新被激活，并演化成巨型滑坡［3］。公路工程建設(shè)中，大量的邊坡開挖，古滑坡體的平衡與穩(wěn)定遭到嚴(yán)重改變和破壞，將會形成規(guī)模更大，破壞更嚴(yán)重，難于治理的新滑坡，嚴(yán)重影響工程建設(shè)質(zhì)量和安全［4］。而我國西部山區(qū)有很多村莊都位于古滑坡體上，并且這些滑坡現(xiàn)在仍在變動中。因此，為減少和避免滑坡災(zāi)害造成的損失，我們首先要勘查圈定滑坡邊界，以估算滑坡致災(zāi)面積，推測滑坡的規(guī)模，為救災(zāi)快速應(yīng)急響應(yīng)提供有效數(shù)據(jù)支持。

目前已有不少學(xué)者開展了滑坡邊界的提取研究，并取得了一定的成果。葉潤青等［5］以遙感影像光譜信息為基礎(chǔ)，利用遙感圖像分類方法將滑坡區(qū)分為植被、土體和水體來圈定滑坡邊界；陳瑩［6］以遙感影像為基礎(chǔ)，通過災(zāi)前和災(zāi)后水體的提取，檢測河流堵塞區(qū)域，獲得準(zhǔn)滑坡范圍，再利用簡單的數(shù)學(xué)形態(tài)學(xué)處理方法與人機(jī)交互判斷，最終識別出滑坡區(qū)域；Mondini［7］利用超高分辨率全色影像和高分辨率多光譜影像，通過對比分析滑前滑后影像的變化，圈定了滑坡區(qū)范圍；魏冬梅［8］以無人機(jī)彩色圖像為基礎(chǔ)，進(jìn)行滑坡邊緣的色度學(xué)差異研究，利用圖像分割技術(shù)對滑坡區(qū)域進(jìn)行提??；葉樹林［9］依據(jù)放射性核素鈾、鐳、氡等在地質(zhì)體中普通存在并有一定差異的基本事實(shí)提出了利用放射性勘查方法來圈定滑坡邊界；徐金明［10］針對巖體滑坡提出了利用流體包裹體來確定滑坡邊界的方法，主要就是對滑坡滑動過程中流體包裹體跡面特征參數(shù)和包裹體熱動力學(xué)參數(shù)進(jìn)行了類別分析與對應(yīng)分析，進(jìn)而確定滑坡邊界。

綜合國內(nèi)外文獻(xiàn)資料來看，勘查滑坡邊界的方法主要有兩種：①通過遙感影像來監(jiān)測滑坡區(qū)河流、植被、陰影、建筑物等信息的變化來確定滑坡邊界；②通過地質(zhì)調(diào)查鉆孔分析其巖性來確定滑坡邊界。但對于古滑坡而言，無法得到滑前遙感影像，對于一些光禿滑坡區(qū)而言，其上河流、植被、陰影、建筑物等信息較缺乏，無法利用影像分析確定邊界，而實(shí)地踏勘將耗費(fèi)大量人力物力，且作業(yè)艱苦，周期較長。針對這一問題，本文從古滑坡區(qū)地形特征出發(fā)，結(jié)合DEM水文分析，提出了一種利用DEM水文分析進(jìn)行滑坡區(qū)域范圍提取的簡便方法。

1 理論與方法

1.1 古滑坡特征分析

在自然地質(zhì)作用和人類活動等因素的影響下，斜坡上的巖土體在重力作用下沿一定的軟弱面“整體”或局部保持巖土結(jié)構(gòu)完整而向下滑動的過程和現(xiàn)象及其形成的地貌形態(tài)，稱之為滑坡［8］，其中發(fā)生在一級河流階地侵蝕時期或更早階地侵蝕或堆積期的滑坡稱為古滑坡［11］。因古滑坡發(fā)生時間較長，受人類活動影響較大，除滑坡壁、滑坡體、滑坡周界外，其他要素難以識別，為方便研究，本文將古滑坡簡化為如圖1所示模型，并對模型中的三大要素進(jìn)行提取。

圖1 古滑坡簡化模型

具有與山體分離的部分斜坡可能向前運(yùn)動的臨空面是滑坡發(fā)育的三個基本地質(zhì)環(huán)境條件之一［11］，故山體滑坡大多發(fā)生在溝谷兩側(cè)，周圍由山脊高陡的巖壁圍繞。受雨水沖刷，人類活動等影響，圖1所示的古滑坡壁坡度已不似新滑坡陡峭，且上沿模糊不清，難以識別；滑坡末端靠近溝谷區(qū)域，受人類活動或河水長期沖刷侵蝕，難以辨識。為提取滑坡范圍邊界，此處將滑坡壁上沿延伸至山脊，以山體與溝谷交界線作為滑坡前緣邊界，二者共同構(gòu)成滑坡區(qū)域范圍。因此滑坡區(qū)域提取分為兩步：滑坡壁上沿提取和滑坡前緣邊界提取。其提取流程如圖2所示。

圖2 滑坡區(qū)域提取流程圖

1.2 提取方法

（1）滑坡壁上沿提取

對于古山體滑坡而言，將滑坡壁上沿延伸至山脊線，即為分水嶺，因此，本文采用流域提取的方法進(jìn)行滑坡壁上沿的提取。目前應(yīng)用最廣泛的于DEM數(shù)據(jù)上提取水文特征的方法是基于地表徑流模擬的方法［12］。地表徑流模擬的方法主要包括三個步驟：一是模擬確定地表徑流方向；二是計(jì)算匯水面積累加矩陣；三是確定匯水面積閾值，定義匯水面積超過閾值的像素為地表河網(wǎng)。

①確定地表徑流方向。目前確定水流方向的方法有兩種：單向流法和多向流法［13］。本文采用前者，即D8算法。D8法假設(shè)單個網(wǎng)格中的水流方向只有8種可能，即流入它的8個鄰域方向。使用D8法確定地表徑流方向有很多實(shí)現(xiàn)方法，本文使用Jenson的方法［14］，通過計(jì)算中心網(wǎng)格與各相鄰網(wǎng)格間的距離權(quán)落差來確定最陡坡度流向。可簡單描述為：在3×3掃描窗口中，設(shè)中心柵格為C，鄰域8個水流可能流出的方向分別為E、SE、S、SW、W、NW、N、NE（東、東南、南、西南、西、西北、北、東北），為方便說明和計(jì)算，每個流向用一個特定編碼來表示：i（i＝1，2，4，8，16，32，64，128；東，東南，南，西南，西，西北，北，東北）［15］。

圖3 流向及其編碼

②計(jì)算匯水面積累加矩陣。確定格網(wǎng)水流方向后，便可進(jìn)行河網(wǎng)提取。匯水面積累加矩陣值直觀上表示區(qū)域地形給水點(diǎn)的水流累積量。其基本思想是：以規(guī)則格網(wǎng)表示的數(shù)字地面高程模型每點(diǎn)處有一個單位的水量，按照自然水流從高處往低處的自然規(guī)律，根據(jù)區(qū)域地形的水流方向數(shù)據(jù)計(jì)算每點(diǎn)處所流過的水量數(shù)值，便得到了該區(qū)域的匯流累積量。

③河網(wǎng)的提取。設(shè)定合適的集流閾值，從水流累積矩陣圖中提取河流柵格圖。當(dāng)柵格特征值大于此閾值時，柵格賦值為1，表示河道，小于該閾值的柵格賦值為0，表示產(chǎn)流區(qū)。各水道按有效的水流方向連接，便可生成連續(xù)的河網(wǎng)［16］。

以上步驟在ArcGIS中即可完成，由于要對5個實(shí)驗(yàn)區(qū)進(jìn)行操作，工作量較大，研究采用ArcGIS10.1中的空間分析建模功能，其模型生成器為設(shè)計(jì)和實(shí)現(xiàn)空間處理模型提供了一個圖形化的建模框架，能將一系列的工具和數(shù)據(jù)串起來以創(chuàng)建高級的功能和流程［12］。河網(wǎng)提取模型如圖4所示。

圖4 河網(wǎng)提取模型

圖5 實(shí)驗(yàn)區(qū)1河網(wǎng)提取結(jié)果

由于這里是對反地形進(jìn)行提取的，所以開始輸入的DEM應(yīng)該為進(jìn)行反地形計(jì)算后的DEM，柵格計(jì)算要根據(jù)需求設(shè)定合適的閾值。圖5是河網(wǎng)提取結(jié)果圖。

圖中水流方向以及河網(wǎng)都是基于反地形數(shù)據(jù)進(jìn)行提取的，由圖5可以看出，河網(wǎng)支流較多，滑坡壁位于河網(wǎng)干流，因此首先進(jìn)行干流提取。由圖5可知，滑坡壁由左右兩條河流主干組成，它們分別為兩邊最長的那根，即由m點(diǎn)和n點(diǎn)出發(fā)的河流，它們最后匯聚于p點(diǎn)，一起流向q點(diǎn)，二者在pq處有一小部分的重疊。以此為基礎(chǔ)，本文對各個節(jié)點(diǎn)進(jìn)行遍歷分析，首先找出一條最長的河流線mq（nq），然后找出一條與該河流重疊長度不大于其本身長度10%的另一條最長河流nq（mq）。刪除pq段，即可得到河流主干mpn。

由于ASTER GDEM（V1）數(shù)據(jù)本身存在一些缺陷，邊界堆疊導(dǎo)致數(shù)據(jù)顯示為直線、坑、隆起或其他異常的幾何形狀等［17］，導(dǎo)致提取的河網(wǎng)信息有時發(fā)生一些錯誤，因此還要結(jié)合坡度以及正地形水流方向信息，對其做進(jìn)一步的判識，對明顯錯誤區(qū)進(jìn)行編輯，從而得到滑坡壁上沿，如圖6所示。

圖6 實(shí)驗(yàn)區(qū)1滑坡壁邊界

（2）滑坡前緣邊界提取

滑坡前緣位于山體與溝谷的交界處，滑坡前緣以上的坡面，其坡度大于滑坡前緣以下坡面的坡度?；谶@一特征，本文利用坡度這一基本地形因子對滑坡前緣進(jìn)行提取，但坡度是一種微觀的地形因子，表示的僅僅是局部地面傾斜，與滑坡前緣具有的宏觀特性不相匹配，所以需要對其作用形式稍作修正［18］。考慮到實(shí)際地形復(fù)雜多變，某些地區(qū)存在局部突變，同時也為了從滑坡體與溝谷的宏觀地貌差異入手尋找滑坡前緣，需要一種表征地面坡度整體變化的統(tǒng)計(jì)指標(biāo)，因此選擇計(jì)算平均坡度，來考察不同區(qū)域范圍坡度陡緩情況。具體步驟如下：

①用圖4河網(wǎng)提取模型所示的河網(wǎng)提取模型對滑坡區(qū)DEM進(jìn)行河網(wǎng)提取，得到河網(wǎng)主干線；

②以河網(wǎng)線為基線，向坡體方向做平行線，平行線間距同DEM格網(wǎng)大小，由河網(wǎng)線到坡體方向依次標(biāo)記為a，b，c，…h(huán)；

③計(jì)算各平行線處DEM的平均坡度；

④通過分析圖7所示的各個面內(nèi)DEM數(shù)據(jù)的平均坡度圖，橫軸為分割面序號，縱軸為坡度值。比較各相鄰面內(nèi)平均坡度差異，找出變化最陡處，即差值最大者對應(yīng)的相鄰面（面d和面e），它們之間的分界線即可視為滑坡前緣邊界。

由以上步驟得到滑坡前緣邊界后，與滑坡壁上沿進(jìn)行疊加，便可得到滑坡區(qū)范圍。

（3）滑坡壁與滑坡體提取

滑坡壁是由滑坡體脫離未移動山體而形成的陡壁，坡度較高；滑坡體脫離山體向下滑動，其坡度較緩，因此在滑坡壁與滑坡體交界處，存在坡度陡變。據(jù)此，可利用滑坡坡度圖實(shí)現(xiàn)滑坡壁的提取，提取算法流程如圖8所示。

圖7 DEM平均坡度圖

圖8 滑坡壁提取流程

提取過程中，首先要對坡度圖像進(jìn)行二值化，滑坡壁提取要選擇坡度大于某一閾值的數(shù)據(jù)，滑坡體則相反。由于圖像二值化會造成邊緣損失，這里將二值化后的圖像與原圖像邊緣進(jìn)行疊加進(jìn)行邊緣填充，以保證滑坡壁邊緣完整性。而滑坡體本身就處于滑坡內(nèi)部，不涉及邊緣，故進(jìn)行滑坡體提取時不需要進(jìn)行邊緣填充。

2 實(shí)驗(yàn)與討論

2.1 實(shí)驗(yàn)區(qū)地理位置

所選的3個滑坡實(shí)驗(yàn)區(qū)位于四川省阿壩藏族羌族自治州理縣，地質(zhì)結(jié)構(gòu)屬龍門山斷裂帶中斷，地形呈蜿蜒起伏的立體單元，地表由西北向東南傾斜。地貌類型為低中山－中山－高山－極高山，是典型的中高山峽谷區(qū)。其地理位置坐標(biāo)如表1所示，圖9展示了其位置分布。

表1 滑坡實(shí)驗(yàn)區(qū)地理位置

圖9 研究區(qū)地理分布

2.2 數(shù)據(jù)

本文采用中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心國際科學(xué)數(shù)據(jù)鏡像網(wǎng)站（http：／／www.gscloud.cn）2009年發(fā)布的全球空間分辨率為30m的數(shù)字高程數(shù)據(jù)產(chǎn)品，此數(shù)據(jù)集利用ASTER GDEM第一版本（V1）的數(shù)據(jù)進(jìn)行加工得來，ASTER GDEM是使用全自動化的方法對所有從Terra衛(wèi)星發(fā)射后到2008年8月獲取的150萬景ASTER近紅外影像進(jìn)行處理得到的［19］，覆蓋范圍為83°N～83°S的所有陸地區(qū)域，達(dá)到了地球陸地表面的99%，其數(shù)據(jù)精度估計(jì)在95%誤差置信水平下，平面誤差為30m，高程誤差為15m～30m［20］。

2.3 實(shí)驗(yàn)結(jié)果討論

通過對滑坡區(qū)DEM的分析，獲得了滑坡實(shí)驗(yàn)區(qū)范圍邊界、滑坡壁以及滑坡體，顯示結(jié)果如圖10所示。

由圖7中顯示的3個實(shí)驗(yàn)區(qū)的結(jié)果圖可以清楚的看出，利用本文所介紹的邊界提取方法獲得的結(jié)果與實(shí)際邊界基本吻合，可以有效地表達(dá)滑坡區(qū)域。為了進(jìn)一步評價(jià)提取的邊界精度，利用人工目視解譯出滑坡范圍邊界，由于目視解譯存在人為誤差，本文以目視解譯邊界的3×3鄰域進(jìn)行判定，如果上述方法提取出來的邊界點(diǎn)落入目視解譯邊界的鄰域中，即可認(rèn)定該邊界點(diǎn)為正確的邊界點(diǎn)，計(jì)算出正確的百分比。經(jīng)計(jì)算得出精度為表2所示。

圖10 滑坡實(shí)驗(yàn)區(qū)提取結(jié)果

表2 提取精度評估

對邊緣定位的精度評估只是一種定性的評估，從精度的結(jié)果來看，本文方法實(shí)現(xiàn)的邊界結(jié)果與目視解譯相差不大，基本滿足解譯要求。

3 結(jié)束語

滑坡壁和滑坡體是滑坡的重要地質(zhì)參數(shù)，它們的確定對滑坡災(zāi)害的防治有著重要的意義。以上結(jié)果表明利用DEM圈定滑坡范圍邊界是有一定理論依據(jù)且方法是可行的、有效的。利用DEM水文分析方法在滑坡地區(qū)，通過河網(wǎng)提取、水流方向分析、坡度分析對滑坡范圍邊界進(jìn)行圈定；利用DEM坡度圖對滑坡壁和滑坡體進(jìn)行提取，此方法無需依靠滑前影像特征，也無需分析紋理特征，克服了傳統(tǒng)方法數(shù)據(jù)較難獲取，紋理特征不易分析等困難，為滑坡體調(diào)查和滑坡制圖提供了一種有效的方法。

［1］黃潤秋.20世紀(jì)以來中國的大型滑坡及其發(fā)生機(jī)制［J］.巖石力學(xué)與工程學(xué)報(bào)，2007，26（3）：433－454.

［2］段永侯.我國地質(zhì)災(zāi)害的基本特征與發(fā)展趨勢［J］.第四紀(jì)研究，1999（3）：208－216.

［3］LIN C W，TSENG C M，TSENG Y H，et al.Recognition of large scale deep－seated landslides in forest areas of Taiwan using high resolution topography［J］.Journal of Asian Earth Sciences，2013，62：389－400.

［4］李勇飛，嵇其偉.淺談古滑坡的判別［J］.西部探礦工程，2006（12）：286－287.

［5］葉潤清，鄧清祿，王海慶.基于圖像分類方法滑坡識別與特征提?。跩］.工程地球物理學(xué)報(bào)，2007，4（6）：574－577.

［6］陳瑩.基于遙感影像的變化檢測方法在滑坡體提取中的應(yīng)用［D］.重慶：西南大學(xué)，2011.

［7］MONDINI A，GUZZETTI F，REICHENBACH P，et al.Semi－automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images［J］.Remote Sensing of Environment，2011，115（7）：1743－1757.

［8］魏冬梅.基于遙感影像滑坡邊界自動提取方法的研究［D］.成都：西南交通大學(xué)，2013.

［9］葉樹林，黃國夫.放射性勘查方法圈定滑坡邊界的機(jī)理及應(yīng)用［J］.中國地質(zhì)災(zāi)害與防治學(xué)報(bào)，2000，11（3）：86－88.

［10］徐金明.流體包裹體法在巖體滑坡周界預(yù)測中的應(yīng)用［J］.土木工程學(xué)報(bào)，2007，40（1）：69－73.

［11］王治華.滑坡遙感［M］.北京：科學(xué)出版社，2012.

［12］李麗，郝振純.基于DEM的流域特征提取綜述［J］.地球科學(xué)進(jìn)展，2003，18（2）：251－256.

［13］黃娜娜，寧芊.基于DEM的數(shù)字河網(wǎng)提取方法及應(yīng)用研究［J］.人民長江，2011，42（24）：50－53.

［14］JENSON S，DOMINGUE J.Extracting topographic structure from digital elevation data for geographic information system analysis［J］.Photogrammetric Engineering and Remote Sensing，1988，54（11）：1593－1600.

［15］張宏鳴.流域分布式土壤侵蝕學(xué)坡長提取與分析［J］.咸陽：西北農(nóng)林科技大學(xué)，2012.

［16］沈中原，李占斌，李鵬，等.基于DEM的流域數(shù)字河網(wǎng)提取算法研究［J］.水資源與水工程學(xué)報(bào)，2009，20（1）：20－23.

［17］張朝忙.中國地區(qū)SRTM3DEM與ASTER GDEM高程精度質(zhì)量評價(jià)［D］.武漢：華中農(nóng)業(yè)大學(xué)，2013.

［18］蘇婧.地形特征線的提取方法研究［D］.西安：西安建筑科技大學(xué)，2012.

［19］趙海濤，張兵，左正立，等.中國及周邊區(qū)域ASTER GDEM與STRM DEM高程對比分析及互補(bǔ)修復(fù)［J］.測繪科學(xué)，2012，37（1）：8－11.

［20］FOURNIADIS I，LIU J，MASON P.Landslide hazard assessment in the Three Gorges area，China，using ASTER imagery：Wushan－Badong［J］.Geomorphology，2007，84（1）：126－144.

Extraction of Ancient Landslide Based on DEM

ZHANG Lu－lu1，LIU Zhao－qin2，LIU Feng1，LI Wei2
（1.Geometric College，Shandong University of Science and Technology，Qingdao 266590；2.State Key Laboratory of Remote Sensing Science，Institute of Remote Sensing and Digital Earth，Chinese Academy of Sciences，Beijing100101）

This paper presents a new method of extracting landslide area based on DEM.Flow direction，slope，ridge line and valley line are extracted based on the landslide model and feature analysis，then the landslide area is obtained through watershed analysis；finally the landslide scarp and the displaced material are acquired based on slope map.The ranges of three ancient landslides of Li county from Sichuan province are extracted from ASTER GDEM data with resolution of 30m，verifying the reliability of this method.

ancient landslides；landslide scarp；displaced material；DEM；watershed analysis

10.3969／j.issn.1000－3177.2015.06.007

TP79

A

1000－3177（2015）142－0037－05

2014－09－29

2014－11－17

973計(jì)劃課題（2013CB733202、2013CB733201）。

張路路（1991—），女，碩士研究生，主要從事遙感信息提取研究。

E－mail：skzll08＠163.com

劉召芹（1973—），男，副研究員，主要從事遙感制圖與導(dǎo)航定位研究。

E－mail：liuzq＠radi.ac.cn