基于InSAR技术的三峡库区巫山—奉节段潜在滑坡识别

徐 帅，王尚晓，牛瑞卿

(1.中国地质大学(武汉)地球物理与空间信息学院，湖北 武汉 430074；2.中国地质调查局武汉地质调查中心(中南地质科技创新中心)，湖北 武汉 430205)

滑坡是指斜坡上的岩土体在自然条件和人类工程活动等因素的影响下，受重力作用沿贯通破坏面整体或者分散地顺坡向下滑动的地质灾害现象，同时也包括处于不稳定状态，可能演化成滑坡的斜坡体[1]。常规的滑坡调查方法主要包括野外调查、目视解译航片等，采用这些方法对典型的大型滑坡进行识别与监测，往往会花费大量的人力和财力，且多数情况下不能及时、准确地确定滑坡范围，而且对滑动不明显、潜在的滑坡的监测和识别能力不足[2-3]。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)作为一种新型的空间地表测量技术，不受自然天气条件的影响，能够全天候、全天时地对地表进行观测[4]。目前该技术已日渐成熟，并运用于由地震、火山、地面沉降、滑坡等造成的地表形变的测量，其不仅能够获取地表发生的微小形变，精度可达到厘米级甚至是毫米级，而且具有分辨率高、时效性高、覆盖范围广等优势，已经成为国内外学者研究的热点[5-9]。

目前，将InSAR技术应用于滑坡灾害的研究主要集中在对特定研究区域已知的单体滑坡进行滑坡体的变形监测和形变特征的分析[10-14],也有不少学者尝试将该技术用于区域滑坡的早期识别以及潜在滑坡范围的确定[15-16]，但都是基于形变点速率的数学统计分析，未考虑到滑坡的形成是复杂的孕灾因子和外界致灾因子综合影响的过程，因此在滑坡的识别精度和潜在滑坡范围的确定方面仍存在一定的问题。另外，常规的永久散射体干涉测量(Permanent Scatterer InSAR,PS-InSAR)技术和小基线集干涉测量(Small Baseline InSAR，SBAS-InSAR)技术有一定的局限性。PS-InSAR技术需要的影像数据较多，可以不用过多考虑由于时间基线和空间垂直基线过大引起的去相关问题，能有效地提取非线性形变信息并消除大气相位的影响，但是处理方法较复杂，对数据的时间覆盖范围要求高；SBAS-InSAR技术可以极大地增加相干影像对的数目，在估计出平均形变速率之后，利用干涉解缠后的相位，通过空间域的低通滤波和时域的高通滤波去除大气效应的影响，最大限度地解决了时间、空间失相关和大气效应对InSAR测量的影响。在轨道精炼步骤中，引入稳定的地面控制点(Ground Control Point，GCP)可以有效地消除平地效应。然而，在选择GCP点时，很难找到质量好的点，而且也不可能在未知区域建立大量可识别的人工稳定性目标，导致在SBAS处理流程中过于主观性，没有定量的标准去评价所选的GCP点，最终干涉结果不理想，生成的形变点密度太低，导致没有足够的形变信息。

本文以三峡库区滑坡灾害多发的巫山—奉节段为研究区，充分利用PS-InSAR技术和SBAS-InSAR技术各自的优势，综合考虑该区域滑坡的形成机制，通过统计区域历史滑坡的发育规律，将形成滑坡的孕灾因子引入到InSAR技术中，实现了高密度、高相关性的形变点的生成，并成功筛选出概率高、置信度水平较高的潜在滑坡体的形变点，最后结合该形变点区域地形地貌特征和高分辨率光学影像识别出潜在滑坡的范围，实现了研究区潜在滑坡的早期识别，用以指导滑坡灾害的预测预警和野外地质灾害的实地调查。

1 研究区概况与SAR数据源

1.1 研究区概况

研究区位于三峡库区巫山—奉节段(见图1)，地理位置为109°33′E～110°11′E、30°45′N～23°28′N之间，长江横贯东西。该地区地处中国地形第二阶梯向第三阶梯的过渡地带，是川东褶皱与鄂西山地会合部位，以中、低山侵蚀峡谷地貌为主，地形起伏较大，地质条件复杂，地层岩性以二叠系下统、三叠系中统巴东组、三叠系下统大冶组和嘉陵江组、志留系中统罗惹坪群第三段以及第四系为主，岩层软硬相间，次级褶皱和断裂褶皱发育，极易孕育滑坡等地质灾害。

1.2 SAR数据源

本文采用的数据源是欧洲空间局新一代C波段的中高分辨率卫星哨兵一号A星(Sentinel-1A)SAR数据，选取2017年3月—2018年11月共48景渐进式扫描对地观测(Terrain Observation by Progressive Scans SAR，TOPSAR)成像模式下的单视复数干涉宽模式(IW)升轨数据，入射角为33.9°，重访问周期为12 d，空间分辨率为5 m(R)×20 m(A)，测绘带幅宽为250 km，极化方式为VV(下载网址为：https://vertex.daac.asf.alaska.edu/),并使用卫星成像21 d后的精密轨道数据进行轨道矫正，定位精度优于5 cm。数字高程模型(DEM)为30 m的SRTM DEM(下载网址为：http://gdex.cr.usgs.gov/gdex/)。

图1 研究区概况及历史滑坡位置图Fig.1 Overview of the study area and historical landslide location map

1.3 地形可视性分析

InSAR技术应用于区域的形变监测时，往往由于地形的影响而产生叠掩和阴影，对落入该区域的形变点造成较大的干扰，因此需要对其进行地形可视性分析。本文采用Cigna等提出的R指数计算方法[17]和Notti等提出的叠掩和阴影计算方法[18]，并结合SAR传感器采集数据的几何参数、视线(light of sight，LOS)方向的角度以及地形的坡度和坡向，对研究区SAR数据的地形可视性进行分析，提取出干扰严重的区域，用于掩膜最终落入该区域的不可信的点。具体计算公式如下：

R=sin[θ-β·sin(A)]

(1)

Visibility=R·Lv·Sh

(2)

上式中：R为地形效应综合指数；θ为视线入射角(°);β为地形坡度(°)；当升轨数据时，A=α-ε[其中,α为地形坡向(°)，ε为卫星飞行方位向与北方向的夹角(°)]；Lv为研究区的叠掩；Sh为研究区的阴影；Visibility为SAR数据的地形可视性。

先采用ArcGIS山体阴影模型和DEM，输入参数为视线(LOS)方位向(γ)和入射角的余角(90°-θ)，计算出研究区阴影(Sh),采用Notti等提出的叠掩计算方法[18]，设置参数为视线(LOS)方位向的补角(γ+180°)和入射角(θ)，计算出研究区的叠掩(Lv);然后利用重分类将计算得到的阴影(Sh)和叠掩(Lv)分别设为0，其他的值设为1。

根据SAR成像特点以及Cigna等在文献[17]中的大量统计试验可得：当R≥sin(θ)，且背向雷达卫星时，地形坡度大于视线入射角余角的区域为主动阴影；当0≤R≤sin(θ)时，出现雷达影像的透视收缩现象；当-1≤R<0，且面向雷达卫星时，地形坡度大于视线入射角的区域，出现雷达影像的主动叠掩现象；当R≥sin(θ)，且非主动阴影和叠掩区时，为地形可视性较好的区域，在这个区域内InSAR测量技术能够获取有效的地形监测结果。因此，本文定义：R≥sin(θ)(Sentinel-1A升轨视线入射角θ=39.6°，即R≥0.56)且非主动阴影和叠掩区为地形可视性较好的区域，从而提高了干涉形变点的可信度。

2 研究方法

2.1 结合PS-InSAR技术的GCP点选择

稳定的地面控制点(GCP)可以有效地估算残余相位，进而消除平地效应，提高结果的可靠性。本文首先利用PS-InSAR技术的优势，依次通过设置振幅离差阈值(不大于0.25)、相干性阈值(不小于0.70)以及空间域的低通滤波和时域的高通滤波，减小干涉相位中的大气延迟相位误差，获取更加可靠的稳定散射点的形变信息；然后，设置形变速率阈值(±1 mm/a)，筛选出相对稳定的PS点；最后，在保证了高相干性和较好稳定性的基础上，进行地理坐标系到SAR斜距坐标系的转换，作为SBAS处理流程中的GCP点，避免了选点的盲目性。

振幅离差(DΔA)是依据SAR影像上的同一位置像元在时间序列上强度值的离散特性进行点目标的识别，其计算公式如下：

(3)

式中：σΔA为干涉图幅度的标准差；μA为干涉图幅度的均值。

定义相干性系数γ的计算公式为

(4)

式中：M、S分别为构成干涉对的SAR影像，其中M为主影像，S为辅影像；“*”为复数共扼;m、n分别为选择窗口的大小；i、j则为像元的坐标值。

2.2 研究区滑坡发育的统计规律

滑坡的形成是受复杂的内外因素共同作用的结果，因此仅仅使用InSAR技术获得地表形变信息来作为滑坡灾害识别与监测的依据，是不具有较强说服力的。由于同一区域内滑坡形成的外部因素大致相同，而且受人为干扰严重，因此本文尝试通过统计与分析研究区的历史滑坡数据，归纳出该区域滑坡的孕灾因子变化规律[19]，并确定各孕灾因子不同等级的滑坡数及其对滑坡的概率贡献值，进而对干涉生成的形变点进行潜在滑坡的概率估算。从三峡库区地质灾害防治工作指挥部1∶10 000灾害地质图数据库中获取了研究区范围内103处历史滑坡的相关数据，该区域以土质滑坡为主，分别统计出滑坡孕灾因子各等级的滑坡数及其对滑坡的概率贡献值，见表1。

2.3 形变点空间域的异常值分析和聚类

通过SBAS-InSAR技术流程生成大量形变干涉点，采用 Anselin Local Moran’sI指数进行形变点空间域的异常值分析和聚类[20]，识别出具有统计

表1 滑坡孕灾因子各等级的滑坡数及其对滑坡的概率贡献值

显著性的高值(热点)和低值(冷点)的空间聚类以及形变数据集范围内的高异常值和低异常值，自动聚合形变点异常数据，确定出适当的空间分析范围，并纠正数据的空间依赖性，计算出形变聚类点置信水平p值和得分，其相关计算公式如下：

(5)

(6)

(7)

(8)

3 研究结果与分析

3.1 研究区潜在滑坡点筛选

根据研究区形变点空间域的异常值分析和聚类处理结果，研究区内共生成了85 439个高相干性的形变点，沿LOS方向的平均形变速率值的范围为-29.83～50.53 mm/a[见图2(a)]。形变速率的正负号表示形变点相对于卫星传感器视线的运动方向，正值表示形变点沿雷达视线方向向靠近传感器的方向移动，负值表示形变点沿雷达视线方向向远离传感器的方向移动，形变速率的绝对值则代表时间域内平均形变速率的大小。沿LOS方向的形变往往不能真实地反映地表变形情况，因此本文利用视坡夹角及余弦定理将沿雷达视线方向的形变速率转化为沿斜坡坡度方向的形变速率，进而客观地反映沿斜坡面的地表形变信息[21]。当雷达视线方向与斜坡面夹角接近90°时，其余弦值接近0，此时沿斜坡面的形变速率趋于无穷大。为了避免出现绝对值极大的异常值，采用Herrera等提出的以余弦值±0.3为合理的阈值，即沿斜坡面的形变速率不能大于沿雷达视线方向形变速率的3.33倍[22]。沿斜坡面的形变速率为正值或者负值分别代表了地表点移动方向沿坡面向上或者向下。依据滑坡体的运动规律，虽然在斜坡坡脚可能出现垂直方向上的正向移动，但水平方向上的位移方向仍应该与沿坡面向下方向保持一致，因此剔除沿斜坡面的形变速率为正值的形变点。

本文选取地形可视性较好范围内的形变点，结合研究区统计的7个滑坡孕灾因子的变化规律以及各孕灾因子对滑坡概率贡献的值，假设这7个孕灾因子对该区域滑坡的概率贡献相同，依据全概率公式(9)，计算出所选形变点的累计概率贡献值P，即各孕灾因子对形变点作为潜在滑坡点的累计概率贡献值P(B)，并筛选出累计概率贡献值大于50%的点作为潜在滑坡点，对这些形变点空间域进行异常值分析和聚类，通过设置p值的大小，选择出置信水平在95%以上的形变点，其空间范围分布见图2(b)。

(9)

图2 沿雷达视线方向上的平均形变速率图(a)、置信水平在95%以上的作为潜在滑坡点的形变点空间 分布图(b)、图2(b)的局部放大图(c)Fig.2 Average rate image of deformation along the line of sight of the SAR(a)，Spatial distribution map of deformation points as potential landslides with confidence levels above 95%(b)，Partial enlargement map (c)

式中:P(B)为各孕灾因子对形变点作为潜在滑坡点的累计概率贡献值;P(Ai)为每个独立孕灾因子对滑坡的概率贡献值(即1/n);P(B|Ai)表示每个孕灾因子类别在发生滑坡事件下的概率值。

3.2 研究区潜在的滑坡提取

结合潜在滑坡点区域的地形地貌特征以及空间域异常值分析和聚类处理结果，在沿雷达视线方向上的平均形变速率图[见图2(a)]上标定出15处概率高、置信水平较好的疑似潜在滑坡范围[见图2(b)]，借助多时相的高分辨率光学影像数据和Google Earth影像数据，对这些潜在的滑坡范围做进一步分析，发现有3处(编号为8、9、13)位于历史滑坡上或者附近，这些区域斜坡有明显的地表形变，存在不稳定状况，出现再次滑动的概率较高；有5处(编号为3、4、10、12、15)是人工修路或者建设房屋造成的地表形变；还有7处(编号为1、2、5、6、7、11)地表有不同程度的破坏，发育成滑坡的概率较大。人类工程活动和地表季节性运动等因素也会影响研究区域内潜在滑坡识别的概率，且由于雷达的可视性有限，并不能完全识别出所有的潜在高概率滑坡，本文利用高分一号影像识别出9处新的潜在滑坡范围，见图3。

图3 利用高分一号影像识别出的潜在滑坡范围标定图Fig.3 Calibration map of potential landslides range based on GF-1 image

4 结 论

本文在滑坡灾害多发的三峡库区巫山—奉节段，利用InSAR技术对潜在的滑坡灾害进行了早期识别研究，探索了合理的干涉测量技术，并结合滑坡孕灾因子和区域地形地貌条件，从滑坡的成因机制方面对识别出的潜在滑坡点进行概率估计，成功识别出概率高、置信水平较好的潜在滑坡范围，可为该区域滑坡灾害的预测预警和野外地质灾害的实地调查提供科学依据，得到如下主要结论：

(1) 联合PS-InSAR和SBAS-InSAR两种技术的优势，获取了高密度、相干性较好的研究区形变结果，共计85 439个形变点，沿LOS方向的平均形变速率值的范围为-29.83～50.53 mm/a。

(2) 充分考虑了SAR的地形可视性，计算了地形综合效应指数R，筛选出R≥0.56且非主动阴影和叠掩区为地形可视性较好的区域，从而提高了干涉形变点的可信度。

(3) 引入了滑坡灾害的7个孕灾因子(地形坡度、地形坡向、植被覆盖度、地形湿度指数、斜坡形态、距离断层距离和地层岩性)，统计分析各孕灾因子对研究区滑坡概率的贡献值，采用全概率公式统计分析每个形变点作为潜在滑坡点的概率贡献值，并通过形变点的Anselin Local Moran’sI指数值对作为潜在滑坡概率贡献值较高的形变点进行空间域异常值分析和聚类处理，获得15处概率高、置信水平较好的疑似潜在滑坡范围，再借助区域地形地貌特征和高分辨率光学影像识别出研究区高概率的9处新的潜在滑坡。

本研究受限于单一的升轨数据，对地形可视性差的区域无法进行有效识别，对于未知区域，由于缺乏地面的监测站点，因此无法完全提取出区域内所有的高概率潜在滑坡体。本文所做的研究并不是追求更高的形变数据精度，而是通过考虑形变数据的异常情况和孕灾因子对滑坡的概率贡献，获取可能发育成潜在滑坡的概率，进而有针对性地调查和监测高概率潜在的滑坡范围，以达到防治滑坡灾害发生的目的。