SBAS-InSAR技术在特大型滑坡形变监测中的应用——以西藏自治区庞村滑坡为例

许东丽， 王 涛， 达瓦桑杰， 杨爱平， 李金锁， 朱建东

(1.四川地质矿产勘查开发局207地质队， 乐山 614000； 2.四川嘉源蓉创地质科技有限公司， 成都 610000；3.西藏自治区地质矿产勘查开发局第二地质大队, 拉萨 850003)

中国是滑坡灾害频发国家，尤其是中国西部地区的滑坡，往往规模大、机制复杂，危害大[1]。滑坡在形成发展过程中，往往伴随着地表形变，通过对这些形变的监测是识别及预测滑坡的重要手段。目前，滑坡监测的主要方法有精密水准测量、全球定位系统(global positioning system，GPS)、位移监测等，但这些方法监测覆盖范围小、工作量大，不适用于大型滑坡监测。近年来，随着合成孔径雷达(interferometric synthetic aperture radar,InSAR)技术的发展，研究者们拓展了其应用范围，将原本主要用于地面沉降监测的技术——星载InSAR用于滑坡监测，取得了一系列良好的结果。

在InSAR变形监测中，早期应用的是D-InSAR[2-5](differential InSAR)技术，在D-InSAR基础上发展了长时间序列InSAR技术并运用在滑坡监测中，主要包括永久散射技术PS-InSAR(permanent scatterer InSAR)和小基线集技术(SBAS-InSAR)。其中，PS-InSAR适用于城市等干涉条件和辐射比较稳定的区域[6]，牛全福等[7]使用PS-InSAR技术对兰州市地面形变进行监测；而SBAS-InSAR是一种基于多主影像的InSAR时序方法，只利用时空基线较短的干涉对来提取地表形变信息，对数据的数量要求较低，适用于山区等干涉条件不稳定的地区。因此，目前主要将其运用于非城区的线状工程以及地质灾害(滑坡为主)的形变监测。文献[8-11]使用SBAS-InSAR技术对川东地区、金沙江流域滑坡和鲜水河地质灾害进行识别和监测;文献[12-13]分别使用SBAS-InSAR技术对地铁及高速公路进行沉降监测，以上研究均证明了SBAS-InSAR在低相干区也具有准确估计形变的能力。

综上所述，中外众多学者采用InSAR技术对不同类型工程和灾害的形变特征进行研究，取得了较满意的研究成果。但将InSAR技术应用于青藏高原特大型滑坡形变特征的研究相对较少。因此以规模大、植被覆盖率低、处于变形加速阶段的西藏自治区山南市庞村特大型滑坡为研究对象，采用SBAS-InSAR技术，选取2018年1月—2020年7月共75景Sentinel-1A数据，对滑坡区域的形变速率、累计形变量、时序形变特征及空间形变特征进行分析，并与地面监测站实测的位移数据进行精度对比，验证了SBAS-InSAR技术应用于青藏高原特大型滑坡形变监测的可靠性，以期为青藏高原特大型滑坡的监测工作提供有益参考。

1 滑坡特征与数据源

1.1 滑坡基本特征

庞村特大型滑坡位于西藏自治区山南市隆子县的加玉乡东侧，加玉河右岸(图1)。该滑坡为一古滑坡，滑坡整体地形呈南高北低，上陡下缓，坡度20°～40°，局部地段可达50°以上，平均坡度约23°，整体面积约0.3 km2，整体体积约1.1×107m2，属于特大型滑坡。滑坡平面形态呈圈椅状，剖面形态呈“阶梯”形，平面形态呈舌形，主滑方向约14°。滑坡体后缘拉张裂缝已基本贯通并形成下错陡坎，基岩出露，最大下错达5.2 m，滑坡后壁可见明显擦痕；滑坡前缘剪出口因G219改道开挖，出现不同程度垮塌(图2)。当地村民于2019年4月发现滑坡裂缝，随着雨季和地震的影响，裂缝持续增大，滑坡整体处于不稳定状态。

图1 研究区域卫星影像图Fig.1 Satellite imagery of the study area

1.2 数据源

1.2.1 雷达数据

本文中数据源选用Sentinel-1A卫星数据，Sentinel-1A卫星是由欧空局发射的免费对地观测雷达卫星，轨道高度693 km，多极化，影像采用TOPS(terrain observation with progressive scans in azimuth)模式成像，重访周期为12 d，雷达波段为C波段，波长5.6 cm。数据空间分辨率30 m，本文中所采用的数据时间跨度为2018年1月3日—2020年7月3日，总计75景影像，数据来源(https://qc.sentinel1.eo.esa.int/)。数据的具体时序信息如表1所示，覆盖范围如图1所示。

表1 数据时间序列列表Table 1 Time series list of data

同时，由于研究区域所处区域地形起伏较大，为消除地形相位对形变相位的影响，本文中采用30 m分辨率的SRTM(shuttle radar topography mission)数据作为辅助数据，计算地形相位分量。

1.2.2 地面监测数据

相关管理部门于2019年8月初在滑坡体上安装4台一体化GNSS监测站以及3台拉线式裂缝计，对滑坡的位移状况进行监测。本次搜集该地面监测仪器从2019年8月1日—2020年7月3日的地面监测数据，覆盖InSAR数据的后半段，对SBAS-InSAR解译结果进行精度验证。地面监测仪器安装位置如图3所示。

1.2.3 降雨量数据

本次实验还搜集了滑坡周边日降雨量数据，该数据的采集地点位于滑坡威胁区外，由雨量计测得。数据时间为2019年8月1日—2020年7月3日，同样覆盖InSAR数据的后半段。

2 数据处理

2.1 时序InSAR形变监测原理

InSAR技术能够获取大面积、高精度的数字高程模型，可用来监测地震、火山、地面沉降等造成的微小地表形变信息。用于形变分析时，InSAR技术可获得监测区域厘米级甚至毫米级的形变量，因此许多研究者将其用于滑坡的缓慢形变监测，分析滑坡的滑动趋势，为滑坡有效治理提供可靠的技术依据。

在滑坡监测领域，最早使用的是D-InSAR技术，它获取同一研究区域不同时间获取的三幅或以上的影像进行干涉处理，获得干涉像对，再对干涉像对进行差分处理，消除地形起伏和几何畸变影响，最终获得该时段的滑坡形变信息。其原理如下：设形变前卫星在A1、A2处获取影像，斜距为r1、r2,θ为卫星在A1位置成像的侧视角；形变发生后在A3位置获取影像，得到斜距r3、空间基线B2,δ1、δ2为基线B1、B2的基线倾角，B‖、B⊥和B′‖、B′⊥分别为空间基线B1、B2沿雷达视线和垂直于雷达视线方向上的分量，ΔRd为地面点的形变量(图4)。则D-InSAR形变监测原理可表达为

图4 D-InSAR原理Fig.4 Principle of D-InSAR

r1r2≅B1sin(θ-δ1)≅B‖

(1)

由式(1)可得，干涉测量相位差与视线方向基线分量成正比。

将在位置A1和A2处获得的地表未形变前的影像进行干涉处理，相位仅包含地形分量，相位差为

(2)

式(2)中：λ为波长；P为重复轨道，P=2。

地表发生形变后在位置A3获取影像，将其与未发生形变时在A1处获取的影像进行干涉处理，得到相位包含地形和形变分量，相位差为

(3)

设视线方向的形变量为ΔRLOS,则由式(1)～式(3)，可得形变量引起的干涉图相位差为

(4)

式(4)中：ΔφLOS可根据相位解缠获得，则根据式(1)～式(4)可计算沿雷达视线方向的形变量ΔRLOS。

由式(4)可得，当ΔRLOS=λ/2时，ΔφLOS=2π，因此当形变量为半波长时，即可引起一个周期的相位变化，说明差分相位对地表形变较为敏感，可用于探测微小形变。

由D-InSAR形变监测原理可知，该技术仅仅使用三期数据进行差分处理，适用于分析缓慢形变滑坡较短时间周期内的形变，当滑坡形变量较大和监测周期较长时，使用D-InSAR技术会导致空间失相干和时间失相干，降低形变监测精度。为了解决以上问题，研究者们提出了长时间序列InSAR技术，主要包括PS-InSAR技术和SBAS-InSAR技术，这两种方法对时空失相干和大气扰动两方面进行限制，显著提高了形变监测能力。其中PS-InSAR利用时间相干原理，对影像数量要求高，适用于城区等获取时间规律、连续或干涉条件和辐射比较稳定的区域；SBAS-InSAR技术利用空间相干原理，对数据数量要求低。

由于本文中研究区域面积较大、地形起伏大、滑坡持续时间长，干涉条件不稳定，因此选用SBAS-InSAR进行形变监测和分析。

2.2 SBAS处理流程

SBAS技术是短基线差分干涉方法，将所有的SAR影像依据空间和时间基线分成不同的短基线子集，各子集内进行差分干涉，再对缠绕的干涉相位进行解缠，根据各相干像元的相位与观测时间的关系，利用奇异值分解方法链接各差分干涉图，抑制DEM(digital elevation model)和大气相位对形变相位的影响，获得最小二乘解，其处理流程见图5。

图5 SBAS-InSAR处理流程Fig.5 Data processing flow of SBAS-InSAR

3 研究结果与分析

3.1 地表形变速率与累计形变量

3.1.1 庞村滑坡及其周边区域形变情况

搜集Sentinel-1A降轨数据75景，基于GAMMA软件处理平台，使用SBAS-InSAR数据处理方法，获得了庞村滑坡及其周边区域的雷达视线(line of sight, LOS)方向形变速率图(图6)，其中负值表示形变方向为沿雷达视线并远离卫星方向。在 2018年1月3日—2020年7月3日期间，庞村及其周边区域最大形变速率约为-128.09 mm/年。

图6 庞村滑坡及其周边地区地表形变速率图Fig.6 The surface deformation rate of Pangcun landslide and the surrounding area

3.1.2 地表形变速率与累计形变量

选取庞村滑坡处的解译结果进行详细分析，结果如图7所示，在2018年1月3日—2020年7月3日期间，庞村滑坡区域的最大形变速率为-38.182 mm/年。从解译结果中可以看出滑坡形变主要发生在右侧；在形变区域存在明显的形变分区，形变速率较大的地区主要位于滑坡右侧的中下部，形变速率大于23.55 mm/年，其中形变速率最大的区域主要有三处，形变速率均大于33.340 mm/年(图7)。

依据形变速率，对滑坡区域2018年1月3日—2020年7月3日的累计形变量进行统计分析，可得滑坡区域内最大累计形变量为110.393 mm(图8)。

图8 庞材滑坡地表累计形变量图Fig.8 The cumulative displacement of Pangcun landslide

3.2 滑坡形变空间分布特征

根据形变速率以及累计形变量图，可将滑坡体划分为3个不同形变程度的分区，如图9所示。

图9 庞材滑坡地表滑坡形变分区Fig.9 The deformation zone of Pangcun landslide

为进一步研究滑坡体形变的空间分布特点，根据初步划定的形变分区，沿滑坡主滑方向取两条剖面线A—A′和B—B′(图10),长度分别为895 m和876 m，从起点开始，每隔25 m取剖面点进行形变速率和累计形变量分析，结果如图11所示。

图10 剖面线分布图Fig.10 Section line distribution

图11 剖面各点形变速率与累计形变量Fig.11 Deformation rate and cumulative displacement at each point of the section line

在剖面线A—A′上，形变速率和累计形变量最大的点距离起点600 m，形变速率达26.036 mm/年，累计形变量73.54 mm；剖面线B—B′上，形变速率和累计形变量最大的点位于距起点500 m处，形变速率和累计形变量分别为28.786 mm/年，累计形变量80.086 mm。

由剖面点形变速率和累计形变量分析结果可知，剖面上形变最大的点位于滑坡前缘区域，与实地调查结果一致，该区域主要为G219国道改造，人工修路开挖造成边坡失稳，形成陡坎并下滑。

3.3 滑坡形变时间分布特征

根据滑坡形变分区结果，在重点形变区(Ⅳ、Ⅴ、Ⅵ表示)内选取两个特征点P1、P2(图12)进行时序分析，研究滑坡形变的时间变化特点，结果如图13所示。

图12 重点形变区及特征点分布Fig.12 Distribution of the main deformation areas and feature points

由图13可以看出，滑坡在两个时段存在明显的形变加速，分别为2018年11月29日—2019年1月22日(Ⅰ期)和2019年10月1日—2019年11月18日(Ⅱ期)，其余时期为均匀沉降。

根据实地调查结果以及降雨量数据：Ⅰ期的形变加速主要是受2018年10月对G219国道的改造工作以及隆子县2018年12月1日地震(28.12°N， 92.62°E，震源深度5 km，据庞村滑坡23 km)影响，国道G219改造工作对滑坡前缘边坡开挖，造成边坡失稳，并在地震的作用下加速下滑。Ⅱ期形变加速主要受2019年7月19日错那5.6级地震(27.67°N， 92.89°E，震源深度10 km，据庞村滑坡72 km)以及研究区域连续降雨(图13)，降低了岩土体的稳定性，滑坡体加速形变。

图13 特征点累计形变量时序图Fig.13 Time series diagram of cumulative displacement of feature points

3.4 精度验证

精度验证部分，将本文中SBAS-InSAR解译结果与地面监测结果进行对比，得到如下结果。

(1)SBAS与GNSS监测结果。如图14所示，将地面GNSS监测站监测结果与SBAS在相应区域的解译结果进行对比发现：四个区域均处于持续形变阶段，并且仪器GNSS-4处的累计形变量比其他三处小，SBAS解译结果可得到相同结论。

图14 SBAS-InSAR与GNSS结果对比Fig.14 Comparison of SBAS-InSAR and GNSS results

(2)SBAS与拉线式位移计监测结果。如图15所示，将地面拉线式位移计监测结果与SBAS在相应区域的解译结果进行对比发现：LF-1与LF-2两个区域处于持续形变阶段，LF-3在前期形变后处于稳定阶段，SBAS解译结果可得到相同结论；同时结合监测点位分布图(图3：LF-1LF-2处于滑坡后缘右侧，LF-3处于滑坡后缘左侧)，验证了SBAS解译结果以及根据解译结果划分的形变分区与地面拉线式位移监测结果的一致性，即滑坡形变主要集中在滑坡右侧。

图15 SBAS-InSAR与拉线式位移计结果对比Fig.15 Comparison of SBAS-InSAR and linear displacement results

由两种对比结果可以看出，地面监测结果在数值上比SBAS解译结果大，这是由于地面监测的主要为滑坡滑动方向的形变量，而SBAS解译结果为沿着雷达视线方向的形变量；并且由于Sentinel-1A数据的地面分辨率较低，区域内形变结果代表了该对应像素分辨率内地物的平均形变量值。但二者在总体形变趋势、形变的空间分布以及时间分布特征上基本一致。

因此运用SBAS-InSAR对免费的Sentinel-1数据进行解译处理，所得的滑坡形变结果可为滑坡形变监测与治理提供可靠的技术依据；并且相对于地面监测站来说，成本更低，结果较为准确，可在部分无地面监测站点的大型滑坡监测中，使用SBAS解译结果代替地面监测结果。

4 结论

以西藏自治区庞村滑坡为研究对象，采用SBAS-InSAR技术进行特大型滑坡形变监测，获得了研究区域2018年1月3日—2020年7月3日时间段内的地表形变速率与累计形变量图，并对结果进行分析，得出了以下结论。

(1)在2018年1月3日—2020年7月3日期间，滑坡区域的最大形变速率为-38.182 mm/年。滑坡形变主要发生在滑坡右侧中下部。

(2)在2018年1月3日—2020年7月3日期间，滑坡区域内最大累计形变量为110.393 mm，滑坡区域存在明显的分界线，可分为三个重点形变区。

(3)空间分布上，剖面线A—A′和B—B′的起点和终点均不存在形变，形变最大的区域分别位于距离剖面A—A′和B—B′起点600 m和500 m处。

(4)时间分布上，在2018年1月3日—2020年7月3日期间，滑坡表面存在两个明显的形变加速时段，分别为2018年12月29日—2019年1月22日和2019年10月1日—11月18日，其余时期为均匀沉降。

(5)精度验证上，本文中SBAS-InSAR解译结果与地面监测所获得的形变结果基本一致，验证了使用免费的Sentinel-1A数据，基于SBAS-InSAR技术对大型滑坡进行形变监测的有效性。

因此，使用免费的Sentinel-1A数据，基于SBAS-InSAR技术对西藏自治区庞村特大型滑坡进行形变监测与分析，充分展示了该滑坡形变的时间和空间分布特征，为探究滑坡的滑动机理以及后续滑坡治理工程提供了可靠的技术依据；将形变监测结果与地面调查以及地面监测结果进行对比，验证了解译结果的可靠性，可将监测结果进行由单点到面域的转换；并且在部分无地面监测站点的滑坡灾害中，可将SBAS-InSAR监测结果在一定程度上代替地面监测结果，为后续类似的特大型滑坡监测提供新的思路和参考。