基于精細(xì)DEM的同震滑坡隱患自動(dòng)提取方法

孫昆襄 丁雨淋 曾浩煒 謝 瀟 朱 慶

1西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都，611756

2浙江中海達(dá)空間信息技術(shù)有限公司，浙江 湖州，313299

3四川視慧智圖空間信息技術(shù)有限公司，四川 成都，610036

同震滑坡是指由地震引發(fā)，震時(shí)突發(fā)滑坡或震后一定時(shí)間發(fā)生的滑坡（地震滯后滑坡）［1］，是發(fā)生在山區(qū)地震的一種重要的次生災(zāi)害類型，其造成的人員傷亡與財(cái)產(chǎn)損失往往占地震災(zāi)害相當(dāng)大的比例［2］。同震滑坡隱患是指“裂”而未“滑”、“松”而未“動(dòng)”的震裂山體，這些受過損傷的山體在之后的余震階段或受到連續(xù)降雨影響及在長(zhǎng)期重力作用下很容易形成滯后滑坡［3］。

傳統(tǒng)的群測(cè)群防人工排查方式，能夠近距離地觀察識(shí)別滑坡隱患，但工作量大、效率低，且同震滑坡隱患多處于地質(zhì)環(huán)境惡劣的高山茂林地區(qū)，群測(cè)群防調(diào)查人員的人身安全難以保障［4，5］。光學(xué)衛(wèi)星遙感能夠獲取大尺度、多期的遙感數(shù)據(jù)，能夠顧及周圍環(huán)境充分結(jié)合背景語(yǔ)義信息識(shí)別滑坡隱患［6-8］。但由于山區(qū)滑坡背景環(huán)境復(fù)雜，對(duì)于“同譜異物”的裸地與滑坡，光學(xué)遙感影像上難以識(shí)別。且光學(xué)遙感手段難以消除植被影響，無(wú)法獲得真實(shí)地表信息，對(duì)于山區(qū)植被茂密地區(qū)無(wú)法適用。SAR（synthetic apenture radar）具有全天候、全天時(shí)的優(yōu)點(diǎn)，InSAR（interferometric SAR）技術(shù)可以進(jìn)行地表形變監(jiān)測(cè)，其形變監(jiān)測(cè)精度可以達(dá)到厘米級(jí)甚至毫米級(jí)，結(jié)合光學(xué)影像可判定滑坡隱患區(qū)［4，5］。但是SAR的植被“穿透”能力差，在植被茂密地區(qū)，無(wú)法通過干涉圖獲得真實(shí)地表形變［9］。

機(jī)載激光雷達(dá)（light laser detection and ranging，LiDAR）是一種主動(dòng)式對(duì)地觀測(cè)系統(tǒng)，能直接得到帶地表3D坐標(biāo)的點(diǎn)云數(shù)據(jù)，其時(shí)間、空間分辨率高，探測(cè)范圍廣，不受日照限制，能全天時(shí)觀測(cè)，具有一定的植被“穿透”能力，通過多次回波技術(shù)和濾波操作，能快速獲取高精度的真實(shí)地面三維信息［9-11］，為滑坡隱患三維地形特征識(shí)別提供了可靠的、高精度的DEM（digital elevation model）數(shù)據(jù)。本文從滑坡隱患地形形態(tài)特征出發(fā)，通過地形形態(tài)約束方法從高精度DEM上提取滑坡隱患，為滑坡隱患提取提供了一種新的解決思路。

1 同震滑坡隱患自動(dòng)識(shí)別原理

基于三維地形特征隱患提取有如下幾點(diǎn)優(yōu)勢(shì)［12］：①三維地形值突變可識(shí)別斜坡的凹陷和突變；②三維地形地面粗糙度可用于識(shí)別斜坡結(jié)構(gòu)形變，裂隙、張裂縫及滑痕等一些地形特征；③對(duì)于高山植被覆蓋區(qū)，三維地形識(shí)別更能發(fā)揮優(yōu)勢(shì)?；诟呔菵EM數(shù)據(jù)的地形形態(tài)特征約束的同震滑坡隱患自動(dòng)提取方法，核心關(guān)鍵技術(shù)包括：①同震滑坡隱患地形特征分析；②提煉地形特征因子，建立同震滑坡隱患地形表征指標(biāo)體系；③同震滑坡隱患地形表征識(shí)別指標(biāo)體系量化，便于計(jì)算機(jī)自動(dòng)提取；④顧及地形形態(tài)特征的同震滑坡隱患自動(dòng)提取，從圖形形態(tài)學(xué)角度出發(fā)，以地形形態(tài)指標(biāo)為約束，實(shí)現(xiàn)同震滑坡隱患區(qū)自動(dòng)提取，如圖1所示。

圖1 同震滑坡隱患自動(dòng)識(shí)別框架Fig.1 Automatic Identification Framework for Co-seismic Landslide Hazards

1.1 同震滑坡隱患地形特征分析

滑坡災(zāi)害的發(fā)生與滑體和滑床沿滑動(dòng)面的相對(duì)運(yùn)動(dòng)有關(guān)，由雙體災(zāi)變力學(xué)可知，滑坡發(fā)生是滑動(dòng)面上的牛頓力突然下降的結(jié)果［13］。同震滑坡，是邊坡對(duì)地震動(dòng)的響應(yīng)，邊坡受力發(fā)生改變，牛頓力突變，從而發(fā)育出滑坡隱患區(qū)或引發(fā)滑坡。發(fā)生地震時(shí)，斜坡所受牛頓力發(fā)生改變，隨著地震的持續(xù)，斜坡逐漸產(chǎn)生形變，斜坡坡腳受到嚴(yán)重破損，坡面中部出現(xiàn)弧形裂隙，隨著地震的不斷增強(qiáng)，沿坡面向下出現(xiàn)剪切裂隙［14］，發(fā)育出完整的“L”或雙“L”型同震滑坡隱患形態(tài)特征（見圖2中黃色虛線標(biāo)記），此時(shí)斜坡受外界因素作用發(fā)育成滑坡的可能性達(dá)到最大，呈現(xiàn)出“裂”而未“滑”，“松”而未“動(dòng)”的特點(diǎn)。本文將主要針對(duì)此種特征的隱患類型進(jìn)行提取。圖2是九寨溝地區(qū)典型同震滑坡隱患的光學(xué)影像和數(shù)字格網(wǎng)模型影像。地震進(jìn)一步增強(qiáng)時(shí)，圖2中的斜坡受力平衡將被打破，坡體中部弧形裂隙將產(chǎn)生錯(cuò)動(dòng)，滑體將沿斜面產(chǎn)生滑移，斜坡進(jìn)一步發(fā)育成滑坡。

圖2 九寨溝典型同震滑坡隱患Fig.2 Typical Co-seismic Landslide Hidden Danger of Jiuzhaigou

此時(shí)滑坡還處于引而未發(fā)的微小變形階段，茂密的植被覆蓋，使得光學(xué)影像上很難發(fā)現(xiàn)滑坡隱患。由LiDAR點(diǎn)云數(shù)據(jù)得到反映地表真實(shí)地形的高精度的DEM數(shù)據(jù)上可以明顯看出同震滑坡隱患特征，其同震滑坡隱患自動(dòng)提取地形表征如表1所示。

表1 同震滑坡隱患自動(dòng)提取地形表征Tab.1Topographical Characterizationof Co-seismic Landslides

1.2 地形表征指標(biāo)體系及指標(biāo)量化

依據(jù)同震滑坡隱患三維地形表征可構(gòu)建三維地形表征指標(biāo)體系。地形屬性特征可通過DEM高程值及其衍生出的坡度值、地表粗糙度進(jìn)行量化表達(dá)［15］。山體陰影圖是地形起伏可視化的一種增強(qiáng)表達(dá)，反映了地表相對(duì)高度情況。對(duì)圖2中的典型同震滑坡隱患形態(tài)特征進(jìn)行分析，可得出隱患區(qū)域形成的弧形裂隙及剪切裂隙均呈現(xiàn)出坡度較大且相對(duì)斜坡周圍區(qū)域發(fā)生陡變的特征。因此，可通過設(shè)置坡度閾值，對(duì)由DEM數(shù)據(jù)生成得到的坡度圖數(shù)據(jù)進(jìn)行特征提取。

地表粗糙度（R）是指一定范圍內(nèi)地形表面面積S地表與其在水平面的投影面積S水平之比，反映了地形形態(tài)的宏觀變化，是反映地形起伏變化和侵蝕程度的一個(gè)重要指標(biāo)。進(jìn)一步可得出其與地形坡度角θ之間的關(guān)系，計(jì)算公式為：

幾何形態(tài)特征是同震滑坡隱患區(qū)別于自然存在的高陡坡區(qū)的關(guān)鍵要素，滑坡隱患裂隙呈現(xiàn)出“L”或雙“L”型裂隙形態(tài)，根據(jù)圖形學(xué)知識(shí)，可通過占空比這一因子來(lái)衡量。占空比（K），即圖形自身面積與圖形外包矩形面積之比，針對(duì)同震滑坡隱患提取，即指滑坡隱患特征輪廓線所包圍的面積S輪廓與特征外包矩形面積S外包矩形之比：

“L”或雙“L”型的滑坡隱患裂隙形態(tài)具有較低的占空比。

對(duì)照組用雌激素軟膏（國(guó)藥準(zhǔn)字：J20090033）治療，清洗外陰后將藥物涂抹于陰道內(nèi)，1次/d；觀察組在對(duì)照組的基礎(chǔ)上給予保婦康栓治療，清洗外陰后將藥物置于陰道內(nèi)，1次/d。兩組均治療14 d。

1.3 地形形態(tài)約束的同震滑坡隱患自動(dòng)識(shí)別

依據(jù)已經(jīng)構(gòu)建的量化指標(biāo)體系，本文將采用“三步走”步驟完成同震滑坡隱患自動(dòng)提?。孩俣捣ㄌ崛〉匦翁卣?；②噪點(diǎn)過濾及裂隙連通性分析；③輪廓提取及地形形態(tài)約束提取滑坡隱患。

山區(qū)地形地勢(shì)起伏，多陡坡、斷崖，同震滑坡隱患裂隙坡度值大且坡度與周圍發(fā)生陡變。首先依據(jù)隱患裂隙分布坡度特征一般在60°～90°之間，選擇適當(dāng)?shù)钠露乳撝祵?duì)坡度圖進(jìn)行二值化處理。

針對(duì)高山地區(qū)背景坡度大、地勢(shì)多起伏，特征提取得到的結(jié)果圖中存在噪點(diǎn)現(xiàn)象，本文基于灰度形態(tài)學(xué)中的腐蝕操作對(duì)其進(jìn)行去噪。灰度形態(tài)學(xué)中的腐蝕是定義在圖像空間上的類似卷積的一種操作，首先選定結(jié)構(gòu)元素矩陣模板，然后將結(jié)構(gòu)元素矩陣模板在圖像矩陣上平移，用圖像矩陣P減去結(jié)構(gòu)元素矩陣B形成的小矩形，取其中最小值賦值到模板中心對(duì)應(yīng)位置的圖像矩陣元素，即可得到濾澡處理之后的圖像結(jié)果，如圖3所示。對(duì)于二值圖像，模板元素較為特殊，為單位矩陣，且模板操作由相減取最小值變?yōu)槟０寰仃嚺c圖像矩陣像元做“與”運(yùn)算，即若整個(gè)模板范圍內(nèi)的原始圖像像元值均為前景值1，則中心像元保持為圖像前景值1不變，否則像素值變?yōu)?，如此可使得靠近背景的所有像素點(diǎn)均被腐蝕（平滑）掉，從而可對(duì)大背景下的孤立異常值點(diǎn)進(jìn)行平滑。模板大小不同，腐蝕程度不同，模板越大腐蝕程度越大，可快速平滑掉噪點(diǎn)，但過大的模板也會(huì)導(dǎo)致目標(biāo)地形異常值被腐蝕掉，破壞目標(biāo)地形異常原有結(jié)構(gòu)，因此需要選擇適當(dāng)?shù)哪０宕笮　?/p>

圖3 灰度圖像腐蝕運(yùn)算示意圖Fig.3 Schematic Diagram of Gray Image Corrosion Operation

滑坡隱患裂隙發(fā)育程度不同，裂隙連接不連貫。連通性分析，即對(duì)斷斷續(xù)續(xù)相連的線段做聯(lián)通操作，使之完全相連，也即平滑對(duì)象邊緣，減少或者填充對(duì)象之間的距離。該操作與上文的腐蝕操作相對(duì)，在形態(tài)學(xué)操作中被稱為膨脹。不同于腐蝕操作的“與”運(yùn)算，膨脹操作是模板與圖像做“或”運(yùn)算。即整個(gè)模板范圍內(nèi)的原始圖像像元值若存在一個(gè)為前景值1，則中心像元即設(shè)置為前景值1。如此可使得前景對(duì)象邊緣得到擴(kuò)大，進(jìn)而可使兩個(gè)分開的前景對(duì)象連接在一起，所能連接的兩個(gè)前景對(duì)象的距離取決于模板大小，模板越大對(duì)于相隔較遠(yuǎn)的前景對(duì)象連接能力越強(qiáng)，但過大的模板會(huì)錯(cuò)誤引入背景為前景對(duì)象，因此需要選擇合適的模板。

相較于陡崖、陡坡，滑坡隱患裂隙地形形態(tài)呈“L”或“雙L”型，是提取同震滑坡隱患的關(guān)鍵要素。根據(jù)邊緣為圖像灰度值變化最顯著的部分，提取每一個(gè)連通區(qū)域的邊緣像素作為該連通區(qū)域的輪廓［16］，并對(duì)輪廓進(jìn)行形態(tài)特征分析，通過占空比指標(biāo)約束，提取出最終同震滑坡隱患區(qū)。

2 實(shí)驗(yàn)驗(yàn)證與分析

2.1 研究區(qū)及數(shù)據(jù)情況

本文研究數(shù)據(jù)為四川省自然資源廳提供的九寨溝熊貓海景區(qū)部分LiDAR點(diǎn)云處理后得到的DEM數(shù)據(jù)［17］。該區(qū)域LiDAR點(diǎn)云數(shù)據(jù)平均密度優(yōu)于50點(diǎn)/m2，光學(xué)影像分辨率優(yōu)于0.2 m，最大限度地保證了地質(zhì)災(zāi)害信息的有效獲取。其處理后制作得到的高精度DEM數(shù)據(jù)分辨率優(yōu)于0.5 m。如圖4所示。

圖4 研究區(qū)數(shù)據(jù)Fig.4 Study Area Data

從光學(xué)影像可以看到，研究區(qū)植被茂密，遮擋嚴(yán)重，無(wú)法識(shí)別出地表真實(shí)情況，DEM數(shù)據(jù)上可大致識(shí)別出3處滑坡隱患特征。由于315°方位角生成得到的山體陰影數(shù)據(jù)，最為符合人眼實(shí)際觀測(cè)體驗(yàn)，本文利用ArcGIS軟件進(jìn)行數(shù)據(jù)處理，設(shè)置方位角參數(shù)315°，高度參數(shù)45°，生成得到研究區(qū)山體陰影數(shù)據(jù)，并疊加等高線圖層，如圖5所示，可明顯觀察到3處典型同震滑坡隱患區(qū)，分別對(duì)3處研究區(qū)采用本文方法進(jìn)行滑坡隱患提取實(shí)驗(yàn)。

圖5 研究區(qū)典型滑坡Fig.5 T ypical Landslides in Study Area

2.2 同震滑坡隱患區(qū)識(shí)別

從圖像處理角度分析，坡度是高程的一階導(dǎo)，坡度特征圖可用于識(shí)別高程突變區(qū)域，進(jìn)一步可對(duì)圖像進(jìn)行二值化操作，區(qū)分目標(biāo)區(qū)和背景區(qū)。對(duì)研究區(qū)DEM數(shù)據(jù)處理，得到如圖6所示各研究區(qū)坡度圖。二值化過程的關(guān)鍵是尋找一個(gè)最優(yōu)閾值，將每一個(gè)像素的特征值與最優(yōu)閾值進(jìn)行比較，若大于閾值則使其為目標(biāo)區(qū)，否則為背景區(qū)。通過分析裂隙形成原理，根據(jù)專家經(jīng)驗(yàn)知識(shí)，裂隙處的坡度一般在60°～90°之間，由于3處研究區(qū)裂隙發(fā)育程度不同，分別設(shè)置3種不同坡度閾值，對(duì)研究區(qū)進(jìn)行地形特征提取，如圖7所示。分析提取效果，最終選定各研究區(qū)坡度閾值為75°、65°、70°。

圖6 研究區(qū)坡度圖Fig.6 Slope Map of Study Area

圖7 不同坡度閾值粗提取結(jié)果Fig.7 Rough Extraction Results for Different Slope Thresholds

2）噪點(diǎn)過濾。

研究區(qū)二值化處理后得到的粗提取結(jié)果存在裂隙特征連接不完整和噪點(diǎn)現(xiàn)象，實(shí)驗(yàn)借助openCV的dilate與opening算法對(duì)粗提取結(jié)果進(jìn)行濾澡處理，得到圖8所示處理結(jié)果。

3）輪廓提取及地形形態(tài)約束提取滑坡隱患。

滑坡隱患裂隙以其獨(dú)特的“L”或雙“L”型地形特征區(qū)別于山區(qū)自然高陡區(qū)域。對(duì)濾噪之后的二值化圖像進(jìn)行連通性分析，提取連通區(qū)域的邊緣輪廓，如圖9所示。設(shè)置占空比閾值分別為0.45、0.35、0.40，對(duì)濾澡后的數(shù)據(jù)進(jìn)行處理，最終提取出隱患區(qū)如圖10所示。紅色線為隱患區(qū)裂隙特征輪廓，綠色框?yàn)殡[患區(qū)外包矩形。對(duì)比提取結(jié)果、解譯結(jié)果及實(shí)地考察結(jié)果，本文方法整體提取效果較好，隱患主體基本圈出；2號(hào)區(qū)域在圈出隱患裂隙主體的同時(shí)，也提取出了坡體中部殘留的滑后裂隙特征。

圖9 輪廓提取結(jié)果Fig.9 Results of Contour Extraction

圖10 同震滑坡隱患識(shí)別結(jié)果Fig.10 Identification Results of Co-seismic Landslide Hidden Danger

3 結(jié)束語(yǔ)

本文提出了一種地形形態(tài)約束的同震滑坡隱患自動(dòng)提取方法，并對(duì)九寨溝熊貓海地區(qū)3處典型同震滑坡隱患區(qū)進(jìn)行了隱患提取。針對(duì)各自研究區(qū)時(shí)，由于隱患裂隙發(fā)育程度不同，坡度閾值、濾澡kernel核窗口大小及占空比閾值均需要進(jìn)行相應(yīng)地調(diào)整，因此進(jìn)一步的研究可以采用一種局部自適應(yīng)的特征提取方法，根據(jù)研究區(qū)地形特性設(shè)置自適應(yīng)的坡度閾值、kernel核大小及占空比閾值。除此之外，還可進(jìn)一步挖掘隱患區(qū)裂隙特征的其他特征表現(xiàn)形式，擴(kuò)充滑坡隱患識(shí)別地形表征，更好地排查同震滑坡隱患。