基于精细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］：①三维地形值突变可识别斜坡的凹陷和突变；②三维地形地面粗糙度可用于识别斜坡结构形变，裂隙、张裂缝及滑痕等一些地形特征；③对于高山植被覆盖区，三维地形识别更能发挥优势。基于高精度DEM数据的地形形态特征约束的同震滑坡隐患自动提取方法，核心关键技术包括：①同震滑坡隐患地形特征分析；②提炼地形特征因子，建立同震滑坡隐患地形表征指标体系；③同震滑坡隐患地形表征识别指标体系量化，便于计算机自动提取；④顾及地形形态特征的同震滑坡隐患自动提取，从图形形态学角度出发，以地形形态指标为约束，实现同震滑坡隐患区自动提取，如图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°之间，选择适当的坡度阈值对坡度图进行二值化处理。

针对高山地区背景坡度大、地势多起伏，特征提取得到的结果图中存在噪点现象，本文基于灰度形态学中的腐蚀操作对其进行去噪。灰度形态学中的腐蚀是定义在图像空间上的类似卷积的一种操作，首先选定结构元素矩阵模板，然后将结构元素矩阵模板在图像矩阵上平移，用图像矩阵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核大小及占空比阈值。除此之外，还可进一步挖掘隐患区裂隙特征的其他特征表现形式，扩充滑坡隐患识别地形表征，更好地排查同震滑坡隐患。