基于多源遥感技术的红层滑坡识别与监测研究

王 庆 芳，郑 志 军，董 继 红，余 天 彬，刘 文，黄 细 超

(1.成都市地质环境监测站，四川 成都 610041； 2.自然资源部成都地质灾害野外科学观测研究站，四川 成都 610041； 3.四川省地质调查院 稀有稀土战略资源评价与利用四川省重点实验室，四川 成都 610081；4.四川省智慧地质大数据有限公司，四川 成都 610081)

0 引 言

滑坡作为一种频繁发生，破坏最为严重的自然灾害之一，具有隐蔽性强、破坏性大、突发性高等特点，严重危害着人民群众的生命财产安全，同时也制约着区域经济的发展[1-3]。红层是一种外观以红色为主色调的陆相碎屑岩沉积地层，龙泉山是典型的低山丘陵红层区。地层主要是侏罗系、白垩系，岩性以红色、紫红色、砖红色砂岩、粉砂岩、泥岩、砾岩为主[4-5]，该区域的植被覆盖率较高，人类工程活动强烈。地层具有强度低、亲水性强、易崩解、软化、膨胀等特性，易发生规模较大的顺层岩性滑坡[6-8]。据调查显示，自2008年以来，该区域发生多起大规模的滑坡灾害，如雷打石滑坡、五家坟滑坡、油榨房滑坡等。区内滑坡具有分布广、隐蔽性强、突发性高、识别难度大等特点。因此很有必要利用遥感技术对该区域的滑坡隐患开展大范围的早期识别研究，对典型的滑坡开展时序监测，分析变形规律及特征。

光学遥感技术因结果准确率高、解译精度高、获取信息丰富、可以回溯历史形变等特点广泛用于滑坡识别、解译分析等工作[9-10]，但是该技术自动化程度低、工作量大。而合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)技术因其非接触、不受云雾天气影响、空间覆盖范围广、监测精度高等优势，广泛用于地表形变监测[11-13]。张路等[14]以四川省丹巴县为例，利用InSAR技术获取的形变速率图及差分干涉图成功识别出17处滑坡隐患，同时总结了影响时序InSAR技术监测效果的主要因素。戴可人等[15]通过SBAS技术获取的形变速率图进行隐患的识别圈定，并讨论了研究区域内Sentinel-1数据可视性分布情况。朱赛楠等[16]利用光学遥感技术、InSAR技术结合物探及地表位移监测等多个手段对金沙江流域色拉滑坡的特征、变形过程、形变机理进行研究，并对其变形趋势进行了预判。李媛茜等[17]基于InSAR技术结合野外资料成功圈定了白龙江流域114处滑坡隐患，分析了滑坡发育特征，并通过对牙豁口滑坡结合无人机数据进行分区研究，证实了InSAR技术在滑坡变形监测的有效性和可靠性。同时一些学者也利用Stacking InSAR技术及SBAS InSAR技术进行滑坡隐患的大范围识别及监测，取得了不错的结果[18-19]。

基于此，本文以四川省成都市龙泉山区域为实验区域，收集了覆盖该区域的ALOS PALSAR-2数据、Sentinel-1数据及高分光学数据，开展多源遥感技术在低山丘陵区红层滑坡隐患的早期识别，并对比不同波段SAR数据的隐患识别情况；基于SBAS技术结合无人机数据及现场调查结果分析了典型滑坡的变形规律及特征，并给出了相应建议。

1 研究区概况及数据来源

1.1 研究区概况

龙泉山位于四川省成都市东侧，西接成都平原区，东部接川东丘陵区，地势中部较高沿龙泉山山脊展布，东西两侧低，本次研究区主要分布在龙泉山中段(见图1)，研究区地貌类型主要为构造侵蚀低山地貌，属于川东低山丘陵区。区内属中亚热带湿润气候区，四季分明、湿润温暖，降雨丰沛，雨季多集中于6～9月，降雨量约占全年的75%。调查分析发现降雨条件对区域的滑坡灾害具有较强的控制作用，滑坡发生的时间主要集中在强降雨期。

图1 研究区位置Fig.1 Geographical location of the study area

1.2 数据收集与处理

本次研究采用多源遥感技术综合识别，利用多源SAR数据和高分光学卫星影像开展区域滑坡隐患早期识别。收集了覆盖该区域2017～2020年的ALOS PALSAR-2降轨数据。同时也收集了可以免费获取的Sentinel-1数据进行研究区内典型滑坡隐患的形变规律监测分析(见表1)。利用AW3D30 DSM数据来消除InSAR干涉处理中的地形相位及辅助SAR影像进行地理编码。Sentinel-1卫星的POD精密轨道星历(POD Precise Orbit Ephemerides)数据被用来辅助Sentinel-1数据的预处理和基线误差改正。

表1 选用的遥感数据信息

同时在数据处理中对与高程不相干的随机大气扰动误差，利用李振洪教授团队研制的通用型卫星雷达在线大气改正系统(Generic Atmospheric Correction Online Service for InSAR，GACOS)来去除[20]。

光学卫星数据主要包括国产高分二号和北京二号数据，用来进行滑坡隐患的光学遥感解译和识别验证。

2 低山丘陵区滑坡识别监测方法

结合工作区滑坡发生的特殊性质，如分布广泛、隐蔽性强、空间尺度较小、受降雨影响大以及突发性强等特点，在基于InSAR技术用于滑坡隐患识别的相关研究和前期数据处理结果分析基础上，本文结合收集到的SAR数据及其他相关资料，针对低山红层区突发性滑坡隐患识别、监测研究，给出了本文的滑坡隐患识别与监测处理流程(见图2)。该流程的关键是对基于地表形变速率解译的滑坡隐患从光学影像中寻找依据，通过SBAS技术判定其变形量及变形规律，通过野外调查、核查识别结果的准确性，并对灾害库进行更新，同时积累该区域的解译标识库。

图2 滑坡隐患识别与监测流程Fig.2 Flow chart of landslide hidden danger identification and monitoring

2.1 Stacking InSAR技术用于滑坡隐患编目

目前开展大范围滑坡隐患早期识别工作，普遍采用Stacking InSAR技术，该技术具有使用影像数量少、结果获取效率高等优势。因此本文在对比不同波段数据获取结果解译情况与已有滑坡隐患库基础上，建议在进行Stacking计算中采用ALOS-2数据与Sentinel-1数据相结合的形式进行。

其中Stacking技术[21]是由Sandwell等在1998年提出，该技术原理是对解缠相位进行加权平均解算，以达到削弱空间上不相关的噪声的影响，从而获取研究区域的平均形变速率[22-23]。主要处理步骤如图3所示。

图3 STACKING-INSAR技术流程Fig.3 Stacking-InSAR technology process

2.2 基于光学遥感技术的滑坡隐患识别

结合工作区地质灾害特征，对滑坡隐患点的识别，可采用对比分析方法和直接识别方法进行判定。对比分析方法首先利用遥感影像对斜坡失稳前地质环境条件和变形迹象进行解译，然后对工作区内地质环境条件相似或者变形迹象明显的地质灾害进行遥感解译。直接识别方法的关键是建立遥感影像识别标志，遥感解译标志包括了地质灾害发育的地质环境条件和斜坡变形特征解译两个方面，如表2所列。

表2 滑坡光学遥感解译识别标志

2.3 SBAS InSAR技术用于滑坡监测分析

对于典型滑坡隐患采用SBAS InSAR技术进行长时间序列监测，获取其变形趋势及变形量等信息，因为Sentinel-1数据时间分辨率高(时间间隔为6/12 d)、免费获取，所以本文利用Sentinel-1数据基于SBAS技术开展典型滑坡隐患的时间序列监测分析。

SBAS技术是由Berardino等于2002年提出的，基本原理是假设在监测时间t0～ts内获取了S+1幅SAR影像，根据干涉组合条件可形成M幅干涉组合，符合式(1)：

(1)

一般根据实际情况设置一定的垂直基线阈值和时间基线阈值对M的值域进行限制，通过这种形式可以降低因垂直基线和时间基线过大引起的失相干现象[23-24]。利用差分干涉计算获得差分干涉相位，用于第j幅干涉相位是根据tA和tB(tA

3 结果分析

3.1 基于Stacking InSAR技术的滑坡隐患早期识别与效果分析

分别利用ALOS-2数据和Sentinel-1数据进行Stacking InSAR处理。对ALOS-2数据进行距离向和方位向采用2×4的多视处理来提高信噪比，由于ALOS-2数据量较少，因此采用全组合的形式进行差分干涉处理，并进行Stacking InSAR解算获得沿雷达卫星视线向(Line of Sight，LOS)年平均形变速率结果(见图4)。对于Sentinel-1数据，在处理中采用8×2进行距离向和方位向多视处理，设置空间基线不大于200 m，时间基线不大于48 d进行干涉对组合，进而提高相干性和运算效率，并通过Stacking InSAR技术获取了形变速率结果。

图4 利用ALOS-2数据获取沿LOS向地表形变速率Fig.4 Surface deformation rate along the LOS direction obtained by ALOS-2 data

根据形变速率图中的条纹形态判识疑似的地灾隐患[29]，结合地形条件及光学影像，对基于不同数据源获取的地表形变监测结果进行滑坡隐患的识别解译，其中ALOS-2数据共识别出9处隐患点(见图4)，野外核查发现8处有明显变形，1处无明显变形。Sentinel-1升轨数据识别出8处，经野外验证5处为滑坡隐患，Sentinel-1降轨数据识别出6处，经野外验证2处为滑坡隐患，Sentinel-1数据经验证正确的隐患点与ALOS-2数据识别的隐患点位置一致，具体滑坡隐患信息如表3所列。

表3 滑坡隐患信息

选择了3处经野外验证为滑坡隐患进行多源遥感数据识别对比分析(见图5)。其中狮子山滑坡和宋家湾滑坡为新增隐患点。通过基于ALOS-2数据获取了这3处滑坡隐患的解译信息，可以看到明显的变形信息。相应地，在Sentinel-1升轨数据中可以看到较为明显的形变信息。其中狮子山滑坡和鄢家沟滑坡的变形信息与ALOS-2数据获取的变形信息具有较好的一致性，只是Sentinel-1是升轨数据，ALOS-2是降轨数据，所以在速率图中图斑颜色正好相反。这3处滑坡隐患在Sentinel-1降轨数据中没有变形信息，故没有展示。

宋家湾滑坡隐患是基于ALOS-2数据获取的新增隐患点，该滑坡在Sentinel-1升轨数据中变形信息并不明显，为此依据图2流程进行验证分析。图6为该滑坡隐患对应的光学影像，从图6中可以看到该滑坡隐患具备地质灾害发育的地质环境条件，但是斜坡变形特征并不明显。为此，提取基于InSAR监测变形强烈的部分，再基于SBAS InSAR技术提取其时间序列变形信息(见图7)。从图7可以看到：该滑坡一直处于变形状态，变形具有明显的周期性，每年7月前后处于加速变形状态，变形时间为工作区雨季，同时2018年1月至2020年12月累积变形量达到了132 mm，佐证了该点为滑坡隐患点，后期通过野外调查验证，确定该点为一新增滑坡隐患点。

图6 宋家湾滑坡光学卫星影像解译信息(北京2号)Fig.6 Optical satellite image interpretation information of Songjiawan landslide(satellite of Beijing 2)

图7 宋家湾滑坡时序点变形信息Fig.7 Deformation information of Songjiawan landslide sequence point

3.2 典型滑坡隐患监测

选择利用时序InSAR技术识别的一处新增隐患点——郑家房子滑坡为例进行隐患监测分析。该滑坡地理位置为104°13′21.93301″E，30°23′0.09891″N；滑坡形态整体呈簸箕形，滑坡前缘高程685 m，后缘高程770 m，高差85 m，纵长约260 m，横宽约150 m，坡向约245°，坡度15°～20°，滑坡平面面积为6.25万 m2，按10 m厚度估算，估算方量为6.25×105m3，为一大型缓倾滑坡。地层岩性为侏罗系中统沙溪庙组(J2s)浅紫红、黄灰色厚-块状中细粒岩屑砂岩、岩屑长石砂岩与紫红色、灰紫色粉砂质泥岩、粉砂岩互层，滑坡结构为一顺向坡，属于典型的低山红层区顺层滑坡。

同时获取了该区域的无人机正射影像图(见图8)，从光学影像图中可以看到滑坡后缘呈圈椅状，可见早期滑动形成的洼地缓坡平台；滑坡两侧以负地形凹地为界，可见沟道微地貌；滑坡前缘以缓坡平台为界，堆积体外凸挤向临空面，陡缓交界线明显。

图8 郑家房子滑坡光学遥感解译(据2021年无人航测数据)Fig.8 Interpretation of Zhengjiafangzi landslide by optical remote sensing(according to unmanned aerial survey data in 2021)

通过对此滑坡进行现场核查，发现滑坡变形区域与监测识别吻合性较好。滑坡变形区主要集中在左侧中部居民区处，房屋局部发生拉裂变形，坡体局部发生溜滑。

据访问坡体上部右侧边界处的一户居民，可知该区域斜坡在2018年至2021年底，每年降雨期坡表均有不同程度的溜滑，该居民房屋无变形，屋后堡坎局部发生溜滑。现场调查坡体中部发育裂缝，裂缝走向NE35°，裂缝延伸长度10～25 m，宽10～20 cm，错落高度5～15 cm。该滑坡为一土质滑坡，斜坡坡向与滑动方向基本一致，属于龙泉山典型的易滑斜坡区。

为了对郑家房子滑坡的变形特征进行分析，采用SBAS InSAR技术对Sentinel-1升轨数据进行处理，获得了该滑坡沿LOS向地表形变速率结果(见图9)，正值表示滑坡位移靠近卫星运动方向变形，负值表示远离卫星飞行方向。

图9 郑家房子滑坡形变速率Fig.9 Deformation rate of Zhengjiafangzi landslide

将图9与ALOS-2数据结果对比后发现，两者在强变形区域具有较好的一致性，图9中形变速率最大达到了-40 mm/a。在滑坡体上选择特征点提取其在时间域的累积形变量结果(见图10)。可以看到特征点的变形信息具有明显的规律性，每年的5～7月之间发生变形，变形速率最大，其他时间段基本趋于稳定，其变形时间与当地的雨季吻合。其中B、C点形变速率最大，累积变形量也最大，其在2018年1月至2020年12月累积变形量达到134 mm。结合现场调查情况可知该滑坡变形主要集中在每年的汛期，在强降雨作用下，容易发生大规模滑坡，威胁坡体居民区、村道，因此需要在雨季加强持续性监测。

图10 郑家房子滑坡体监测点累计形变曲线Fig.10 Cumulative deformation curve of monitoring points of Zhengjiafangzi landslide mass

4 讨 论

4.1 低山丘陵区雷达数据适用性及SAR数据选取分析

3.1节中提到Sentinel-1数据和ALOS-2数据识别滑坡隐患的数目及准确率不一致，经分析主要是受差分干涉相干性的影响。ALOS-2数据搭载L波段的传感器，波长较长，对地表植被具有较好的穿透性。Sentinel-1数据搭载C波段，使得在植被茂密区域相干性差，尤其是夏季更甚。分别在5月与8月，每个月选取2景覆盖研究区域时间间隔12 d的Sentinel-1数据进行差分干涉处理，获取了相干系数对比图(见图11)。从图11中可以看到时间间隔12 d之内，相关系数图在山区相干性较差，大部分区域相关系数低于0.2，因此在进行相位解缠的时候，绝大部分会被掩膜掉，同时滑坡易发区域位于山区，因此会造成监测区域缺少有效信息。结合图7可知该区域的滑坡隐患变形主要发生在5～7月，但是这段时间因为植被茂密使得基于Sentinel-1数据进行差分干涉处理，获取到的监测信息效果较差。

图11 不同时间段相关系数(间隔12 d)Fig.11 Correlation coefficient diagram for different time periods(interval of 12 days)

为了提高该区域的隐患识别能力和监测效果，经过上面的对比和结合文献[30]，认为在该区域进行滑坡隐患的早期识别，应选择分辨率更高或者波长更长的SAR数据，比如日本的ALOS-2卫星数据，阿根廷的SAOCOM-1A卫星数据，以及即将发射的国产SAR卫星。而Sentinel-1数据因高时间分辨率、免费获取可以用来对典型滑坡隐患开展时序监测分析。

4.2 Stacking InSAR技术适用性分析

在本文中用Stacking 技术进行隐患识别的前提是假设识别对象变形为线性形变，进而获取了地表形变结果。但是对该区域的滑坡通过SBAS InSAR技术进行时序监测，发现其变形具有一定的周期性，这与假设相矛盾，仍然可以识别出滑坡隐患的原因有：① 收集到的ALOS-2数据集，在夏季有影像数据，使得可以很好获取降雨期间的变形信息。② Sentinel-1数据在5～8月虽然相干性较差，但是时间分辨率高，差分干涉组合多，因此在加权平均过程中可以获取到有效变形信息。

需要指出来的是Stacking技术在该区域进行滑坡隐患识别时，可能因为平均计算使得获取的地表形变形变速率结果值偏小，部分变形速率较小的滑坡隐患存在漏判现象，如文中的宋家湾滑坡，在基于Sentinel-1数据的Stacking结果中该滑坡隐患变形效果并不明显，但是在累积形变量中均有较强烈的表现。

5 结 论

本文以龙泉山区域为例，采用时序InSAR技术利用多源SAR数据，结合光学遥感技术利用高分光学数据对该区域滑坡隐患开展识别研究，基于SBAS技术对典型滑坡隐患开展时序监测分析，并对识别的隐患进行野外验证。主要有以下结论：

(1) 利用Stacking技术、结合光学遥感技术成功识别出8处滑坡隐患，并通过野外验证验证了其可靠性，发现ALOS-2数据的识别效果优于Sentinel-1数据，也证明了时序InSAR技术用于该区域的滑坡隐患识别具有可行性。

(2) 对基于ALOS-2数据识别的郑家房子滑坡，通过时序监测发现其最大形变速率为-40 mm/a，最大累积形变量为134 mm，并且其变形规律与汛期相吻合。通过野外调查发现，该滑坡在强降雨作用下，容易发生大规模失移，威胁坡体居民区、村道，因此需要在雨季加强持续性监测，同时结果表明多种时序InSAR技术相结合用于滑坡隐患识别监测具有较强可行性和可靠性。

(3) 通过对比不同波段SAR数据的识别结果及影响因素分析，发现影响该区域识别效果主要是干涉相干性，因此建议该区域采用分辨率更高或者波长较长的SAR数据开展滑坡隐患识别工作。

(4) 通过讨论Stacking InSAR技术的适用性，分析了部分滑坡隐患漏判的原因，建议采用多种时序InSAR技术相结合形式开展滑坡隐患识别与监测工作。