基于精細DEM的同震滑坡隱患自動提取方法

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

1西南交通大學地球科學與環境工程學院，四川 成都，611756

2浙江中海達空間信息技術有限公司，浙江 湖州，313299

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

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

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

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

1 同震滑坡隱患自動識別原理

基于三維地形特征隱患提取有如下幾點優勢［12］：①三維地形值突變可識別斜坡的凹陷和突變；②三維地形地面粗糙度可用于識別斜坡結構形變，裂隙、張裂縫及滑痕等一些地形特征；③對于高山植被覆蓋區，三維地形識別更能發揮優勢?；诟呔菵EM數據的地形形態特征約束的同震滑坡隱患自動提取方法，核心關鍵技術包括：①同震滑坡隱患地形特征分析；②提煉地形特征因子，建立同震滑坡隱患地形表征指標體系；③同震滑坡隱患地形表征識別指標體系量化，便于計算機自動提??；④顧及地形形態特征的同震滑坡隱患自動提取，從圖形形態學角度出發，以地形形態指標為約束，實現同震滑坡隱患區自動提取，如圖1所示。

圖1 同震滑坡隱患自動識別框架Fig.1 Automatic Identification Framework for Co-seismic Landslide Hazards

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

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

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

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

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

1.2 地形表征指標體系及指標量化

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

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

幾何形態特征是同震滑坡隱患區別于自然存在的高陡坡區的關鍵要素，滑坡隱患裂隙呈現出“L”或雙“L”型裂隙形態，根據圖形學知識，可通過占空比這一因子來衡量。占空比（K），即圖形自身面積與圖形外包矩形面積之比，針對同震滑坡隱患提取，即指滑坡隱患特征輪廓線所包圍的面積S輪廓與特征外包矩形面積S外包矩形之比：

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

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

1.3 地形形態約束的同震滑坡隱患自動識別

依據已經構建的量化指標體系，本文將采用“三步走”步驟完成同震滑坡隱患自動提?。孩俣捣ㄌ崛〉匦翁卣?；②噪點過濾及裂隙連通性分析；③輪廓提取及地形形態約束提取滑坡隱患。

山區地形地勢起伏，多陡坡、斷崖，同震滑坡隱患裂隙坡度值大且坡度與周圍發生陡變。首先依據隱患裂隙分布坡度特征一般在60°～90°之間，選擇適當的坡度閾值對坡度圖進行二值化處理。

針對高山地區背景坡度大、地勢多起伏，特征提取得到的結果圖中存在噪點現象，本文基于灰度形態學中的腐蝕操作對其進行去噪?；叶刃螒B學中的腐蝕是定義在圖像空間上的類似卷積的一種操作，首先選定結構元素矩陣模板，然后將結構元素矩陣模板在圖像矩陣上平移，用圖像矩陣P減去結構元素矩陣B形成的小矩形，取其中最小值賦值到模板中心對應位置的圖像矩陣元素，即可得到濾澡處理之后的圖像結果，如圖3所示。對于二值圖像，模板元素較為特殊，為單位矩陣，且模板操作由相減取最小值變為模板矩陣與圖像矩陣像元做“與”運算，即若整個模板范圍內的原始圖像像元值均為前景值1，則中心像元保持為圖像前景值1不變，否則像素值變為0，如此可使得靠近背景的所有像素點均被腐蝕（平滑）掉，從而可對大背景下的孤立異常值點進行平滑。模板大小不同，腐蝕程度不同，模板越大腐蝕程度越大，可快速平滑掉噪點，但過大的模板也會導致目標地形異常值被腐蝕掉，破壞目標地形異常原有結構，因此需要選擇適當的模板大小。

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

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

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

2 實驗驗證與分析

2.1 研究區及數據情況

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

圖4 研究區數據Fig.4 Study Area Data

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

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

2.2 同震滑坡隱患區識別

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

圖6 研究區坡度圖Fig.6 Slope Map of Study Area

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

2）噪點過濾。

研究區二值化處理后得到的粗提取結果存在裂隙特征連接不完整和噪點現象，實驗借助openCV的dilate與opening算法對粗提取結果進行濾澡處理，得到圖8所示處理結果。

3）輪廓提取及地形形態約束提取滑坡隱患。

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

圖9 輪廓提取結果Fig.9 Results of Contour Extraction

圖10 同震滑坡隱患識別結果Fig.10 Identification Results of Co-seismic Landslide Hidden Danger

3 結束語

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