SBAS-InSAR技术在西藏江达县金沙江流域典型巨型滑坡变形监测中的应用

杨成业，张 涛，高 贵，卜崇阳，吴 华

（1.西藏大学工学院,西藏 拉萨 850000；2.西南交通大学地球科学与环境工程学院,四川 成都 611756）

0 引言

进入20世纪后，随着人类活动对环境和气候的影响逐渐增大，生态环境以及各种地质条件发生了极大的变化，滑坡、泥石流、崩塌等地质灾害频频出现。其中，滑坡是斜坡岩石体以及大面积冰雪覆盖体沿着贯通的剪切面所发生的体位滑动现象[1]，是一种常见地质灾害。另高山峡谷地带发生的滑坡灾害，以往由于气候、天气、地形、交通、通讯等因素，很难做到早期预警和提前防范。

星载合成孔径雷达（Synthetic Aperture Radar，SAR）作为一种主动式微波遥感观测手段，以其全天侯全天时的对地观测能力，已成为一种不可或缺的大地测量技术。尤其雷达干涉测量（Interferometry Synthetic Aperture Radar，InSAR）技术的兴起，更是极大的推进了SAR 技术在地质灾害监测的应用[2−3]。早期阶段，主要利用SAR 影像的灰度纹理特征进行滑坡区域的识别与提取工作[4]，随着SAR 技术不断完善，InSAR技术逐渐在滑坡灾害监测中占据了主导地位，可以利用SAR 影像间的相干性获取滑坡区的形变特征。

然而不同时段的SAR 影像间易出现低相关或失相关的现象以及大气对雷达波的延迟效应使得D-InSAR技术的发展受到了一定程度的制约[5−8]。针对这种现象，时序InSAR技术逐步兴起[9−11]，它可以实现对目标区域的连续时空观测，能够获取地面较高精度的地形数据以及地表时序形变数据，在很大程度上提升了SAR数据干涉测量的精度，与此同时，意大利Ferretti 等[12]提出了永久散射体差分干涉测量（Permanent Scatterer Interferometry SAR ，PSI）技术，Berardino 等[13]提出了小基线集技术（Small Baseline Subset,SBAS）。唐青松等[14]利用PS-InSAR技术对新磨村滑坡灾害发生前后的形变进行了监测，并根据滑坡灾害发生前后的SAR 影像获取了该区域的形变速率图；姚佳明等[15]利用SBASInSAR 数据采用雷达干涉测量技术对缓倾煤层开采诱发顺层岩体地表形变的模式进行了研究，证实了SBASInSAR技术识别地表形变的准确性；赵富萌等[16]利用SBAS-InSAR技术成功获取了Karakorum公路沿线滑坡灾害；Zhang 等[17]实现了一种将经验模型与SBASInSAR 数据相结合的滑坡潜力体积与面积评估的方法，成功对甘肃省中部黄河黑泰河潜在滑坡进行了分析评价。这些案例充分说明SBAS-InSAR技术在区域形变监测及获取形变参数中的应用潜力。

西藏江达县金沙江流域是滑坡的易发多发区，另外该区域滑坡会导致金沙江堵塞，对当地和下游基础设施建设带来极大的危害。鉴此，文中选取了2017年7月—2018年12月的23 景Sentinel-1A 影像数据，使用SBASInSAR技术对西藏江达县金沙江流域典型滑坡灾害进行形变特征和趋势分析，为该区域下一步地灾监测和防治提供参考。

1 技术原理与数据处理流程

时序InSAR技术获取地表目标连续形变信息的关键在于获取时序影像中散射特征表现稳定的散射体[10]，通常越稳定的地物，散射特性越强，其回波信号受干扰能力也就越强。通过对稳定散射体干涉相位的统计分析与建模，实现干涉相位中各个分量的有效估算与分离，在此基础上对各个分量进行拟合建模，进而获取相应的解算模型，最后在最小化其他干扰相位的情况下获取研究区域干涉相位中的绝对形变相位[18−21]。在时序InSAR 处理方法中，PS-InSAR技术对地面稳定散射体的密度要求较高[22−23]，考虑到研究区域人工建筑等稳定散射体密度有限，文中研究使用SBAS-InSAR技术获取研究区域的形变信息，SBAS-InSAR技术的主要处理流程如图1所示。

图1 SBAS-InSAR技术处理流程Fig.1 Processing flow of SBAS-InSAR technology

在干涉对组合和高相干点的选取上，SBAS-InSAR技术与其他时序InSAR技术有所不同[13,24]。为了最小化时空去相干噪声的影响，SBAS-InSAR技术采用时空基线阈值的方式对所有影像进行最优组合来获取M景干涉像对，这种方式在控制时间和空间的情况下有效的降低了组合影像之间地物变化的显著性，在提升干涉相干性方面具有一定的优势。在不考虑高程误差、大气延迟、噪声等干扰相位的情况下，每个干涉像对生成的干涉图所对应的干涉相位可表示为[10−11]：

式中：Δφi(x,y)、ΔφA(x,y)、ΔφB(x,y)—差分干涉相位、主影像相位、副影像相位信息；

TA、TB——表示获取主、副影像的时间；

x、y——表示影像在方位向、视线向的坐标；

λ——雷达波长。

经过解缠后时序SAR 影像的相位矩阵φT以及所有干涉图的差分干涉相位矩阵ΔφT可以表示为：

式中：φ(ti)——表示为每一景SAR 影像中的相位值，；

φT——为时序SAR 影像相位值。

将主副影像的干涉相位进行差分，并将干涉相位和差分干涉相位按照时间序列排列，主影像为IE，副影像为IS，将公式二进行转换从而得到：

其中形变相位可以简化含有N个未知数的M个方程：

式中：Δφ——差分干涉相位。

矩阵G[M×N]的每一行对应每一个差分干涉图，每一列对应不同的SAR影像数据，所以有矩阵中其他元素为0[25]，即可以将[M×N]阶矩阵G表示为：

当生成一系列干涉对处于同一小基线子集中时，M≥N且N为G的秩，通过最小二乘法即可得到相位矩阵：

对于不可能所有影像的干涉对都在同一个基线子集中的现象，将已有的SAR 数据依据阈值，组合为具有一定数量的子集，继而GTG就是一个降秩后的矩阵，从而可以降低噪声和相干性对干涉对的影响。对于矩阵为不满秩矩阵的现象采用奇异值分解的方法，即将奇异矩阵分解为：

式中：U——正交矩阵；

S——对角矩阵；

VT——平均相位的速率。

继而可以将相位转换为平均相位速度：

在将形变相位构建成方程组的基础上，利用奇异值分解的算法求得研究区域线性形变速率和高程误差相位；将干涉相位中的残余相位进行分离，并进行二次解缠，通过空间域滤波的方式估算并去除大气相位和噪声相位。干涉相位中的各个成份求出后，将原始相位时间序列与求出的各个干扰相位进行作差，最终求得研究区域的线性形变相位和非线性形变相位。

2 研究区域概况

江达县位于西藏自治区昌都市东部。平均海拔3 800 m，最低海拔为2 800 m，最高海拔为5 436 m。横断高山峡谷地貌，山峰陡峭，地形复杂多样，流水及冻融风化作用是塑造山地形态的主要外营力。地势由西北向东南倾斜，西北高东南低。江达县内河流分布广泛，其中金沙江流域分布在江达县东北以及东南县界（图2）。

图2 研究区范围Fig.2 Scope of the study area

金沙江河床纵断面的起伏不平，加上人类活动对河流地貌的影响，导致江达县成为地质灾害易发区域，特别是江达县东部金沙江流域地质灾害发生数量最多，文中选取了江达县东部主要灾害易发区域作为探测目标。

3 试验数据

文中选取的数据为欧洲航天局Sentinel-1A 卫星影像，该卫星载有C 波段合成孔径雷达，在不受时间、空间与天气影响的情况，可以连续提供时间间隔为12 d的全球卫星影像数据。文中研究使用了升轨Sentinel-1A 雷达数据，成像模式为IW，幅宽为250 km，极化方式为VV，影像数据斜距向分辨率为2.33 m，方位向分辨率为13.97 m，入射角为33.8°。时间跨度为2017年7月4日—2018年12月26日共23 景影像数据（表1）。数字高程模型作为SBAS-InSAR 处理流程的辅助数据，可以用于去除地形相位以及后续高程误差相位的计算，并采用由美国太空总署（NASA）和国防部国家测绘局（NIMA）联合发布的SRTM3 数据，其分辨率精度为90 m。选择时间基线长度为365 d，空间基线的阈值设置为理论空间基线的45%，在以上基础上23 景影像共获取了62 个干涉组合，其中主影像日期为2017年10月8日，干涉对以及基线连接情况如图3所示。

表1 Sentinel-1A 数据集Table 1 Sentinel-1A data set

图3 干涉对以及基线连接图Fig.3 Interference pairs and baseline connection diagram

4 典型滑坡形变监测结果与形变趋势分析

4.1 SBAS-InSAR 形变监测结果

在SBAS-InSAR 中，区域内低、失相干区越少，后续反演出的形变相位信息越可靠。换言之，基于SBASInSAR技术获取研究区域的地表形变监测结果的准确性受高相干点的影响较大。因此，根据相干系数阈值法选取时间跨度内的干涉影像中的较为稳定的散射单元（图4）。干涉图中像元点相干系数值的高低是由地物散射特性的变化所引起，林地、植被覆盖和河流区域表现为低相干或失相干，靠近或处于形变区域的像元点同样表现为较低的相干性，对应图4 中的红色区域，裸地和村落区域表现为高相干，对应图4 中的蓝色区域。根据相干点分布图可以看出，区域内中部和东南方向存在较多的像元空洞位置，但图4 中低相干区域较少，因此在一定程度上说明了文中选取相干点的可靠性。

图4 沃达村和白格村滑坡区域的相干点分布Fig.4 Distribution of coherent points in the landslide area of Woda and Baige Village

根据研究区内现有的研究资料，选取了沃达和白格2 处滑坡区域的谷歌影像（图5），与获取的2 处相干点的分布进行对比分析。沃达滑坡区域存在多处含有少量植被覆盖或裸地区域，由于该处滑坡属于古滑坡，年代久远，滑源区出现较为稀疏的植被覆盖，中部区域基本无植被覆盖，存在零星人工建筑（图5）。区域相干点的分布结果显示，在时间跨度内除金沙江沿岸出现的几处低相干区域，滑坡的中部及斜坡出现了较为明显的低相干特征。结合光学遥感影像分析，斜坡区域为植被稀疏地带，所以表现为明显低相干特征的区域可以排除湖泊和植被的可能性，属于高相干地物。在InSAR技术中，地物的形变是引起SAR 影像之间相干性变化的主要因素，同时斜坡地处滑源区下方，地表地物多为古滑坡堆积体所构成，表现为较松散的特征。因此结合斜坡所处的地理位置以及其地物特征分析，斜坡区域存在形变，图4（b）中沃达滑坡中部红色区域。

图5 白格滑坡与沃达滑坡区域地表覆盖特征Fig.5 Surface cover characteristics of Baige landslide and Woda landslide areas

Google 影像中显示白格滑坡区域的植被覆盖量较大，并且在2018年之前靠近后缘的位置已经出现多处明显的滑动条带，因此白格滑坡区域出现多处的失相干区域。相干系数图中显示靠近白格滑坡最顶端位置表现为高相干特征，结合Google 影像分析，见图5（b），滑坡后缘以上区域存在较多人工建筑（如公路、房屋等）以及裸地，因此在相干系数影像中表现为蓝色，而低相干区域主要集中在滑坡中部以及底端位置，见图4（b）。此外，因为滑坡发生垮塌时形成的大量碎屑，在进入金沙江内继续向前运动，将对岸植被冲刷殆尽，暴露出岩石碎屑等高相干散射体，因此白格滑坡对岸出现高相干区域。高相干区域的分布情况以及区域内的相干系数变化进一步验证了相干系数阈值法所保留的高相干区域的可靠性，为后续反演目标区域形变特征的可靠性奠定了基础。

在确定依据相干系数阈值法获取的相干区域准确性的基础上，对干涉相位中的残余地形相位以及形变相位进行了估算，共包括2 次形变反演。基于SBASInSAR技术获取的研究区域内的整体形变分布结果（图6），可以看出该区域在所选时间跨度内的最大累计形变量为−134.54 mm，最大抬升量为52.34 mm，图6 中颜色越接近蓝色的位置表示沉降量越大，对应地表发生塌陷或错位，靠近红色条带的区域表示地表出现的抬升现象，对应该位置处发生物体堆积或由于地壳运动导致的地表抬升。根据图6 中的像元点对应条带的颜色可以看出整个范围大部分区域处于稳定状态，中部和北部区域存在轻微变形，形变量分布在−30～10 mm。发生明显形变的区域主要分布在东部金沙江沿岸附近，最大形变量达到−134.121 1 mm，如图6 中的标注的沃达滑坡和白格滑坡位置。

图6 研究区的整体形变分布结果Fig.6 Results of the overall deformation distribution in the study area

前文根据Google 影像对2 滑坡区域内的地物分布情况进行了灾害解译[25]，根据解译结果，白格滑坡位置存在明显的地表位移现象，沃达滑坡由于地表地物变化不明显，在Google 影像中尚未显现出明显的形变错位特征。结合现有资料，白格滑坡分别于2018年10月与11月发生2 次大型山体垮塌，并导致金沙江主河道发生堵塞，属于典型的巨型滑坡。图6 中的形变探测结果为2017年7月4日至2018年12月26日的累计形变量，因此截止2018年12月26日，白格滑坡范围仍然存在明显的形变位移现象。岩比乡位置处的形变漏斗较为集中，时间跨度的累计形变量超过100 mm，该位置属于较为典型的巨型滑坡隐患区域，自1985年出现形变，至今未发生大规模的垮塌。结合Google 影像可以看出沃达滑坡范围存在多处植被覆盖或裸地区域，滑源区存在较为稀疏的植被覆盖，中部区域与形变漏斗分布的位置较为接近，区域内基本无植被覆盖，存在零星的人工建筑。鉴于其特殊的地理位置，一旦发生垮塌，势必造成金沙江的堵塞，进而引发灾害链，影响下游工程建设与生产生活的进行。

4.2 形变趋势分析

根据前文对江达县整体形变特征的分析结果，进一步对沃达和白格两处典型的巨型滑坡在时间跨度内的累计形变和时序形变特征进行了分析。

从沃达滑坡的累计形变图6（a）中可以看出，目前整个滑坡范围内滑源区处于稳定状态，存在明显形变的位置主要集中在中部和中前缘，且形变范围较为广泛，最大形变量超过−100 mm（图7 中紫色区域），紫色周边区域的形变量较多分布在−60～−100 mm。为核实形变解译数据的准确性，进行了实地地质调查，根据Google 影像和现场勘测照片，形变位置较为陡峭且发育多处冲沟，其中，表现为明显形变的区域及其附近发育出多条横向拉张裂缝，见图7（b）；最大裂缝宽度达到60 cm，滑坡前缘区域出现“马刀树”，见图7（c），根据树木存在时间可以看出该变形区域存在时间至少在5年以上；其周围坡度为30°～40°的区域最大错落超过2 m，见图7（d）。对比滑坡范围的形变分布特征和现场地表特征可以看出沃达滑坡中部和前缘位置存在滑坡复活迹象，目前处于形变活跃状态。

图7 沃达滑坡形变速率图及现场照片Fig.7 Woda landslide deformation rate map and site photos

白格滑坡在2018年发生2 次大量级滑动，在此之后区域内存在大量的残留体，区域内整体又表现出明显的变形，一半以上的范围形变量超过−30 mm（图8浅蓝色区域），滑坡后缘和中部存在多处形变量超过−90 mm 的区域（图8 蓝色和紫色区域）。根据图5 Google 影像显示，滑坡体在发生第1 次脱落前，区域内已存在明显的坡体破碎现象，由于后缘位置的坡体不断出现垮塌，碎屑物在向下滑行的过程中堆积至中下部区域，进而表现为明显的抬升现象（图8 黄色和红色区域）。结合现有的形变位移解译资料，白格滑坡具有分块和分级的形变特征，图8 中同样呈现出类似的变形迹象（图8 红色线圈）。根据目前形变活跃的分布位置可见形变区域位于滑坡的边缘位置，其稳定性受周边地质的影响，加上其自身是被扰动后的残留体，稳定性较差，在图8 中显示为明显的位移特征。结合现场勘测照片可以看出滑坡后缘位置存在多处裂缝和沉降现象见图8（a），与GD1 位置的形变相切合。实地地质调查复核结果与SBAS-InSAR 形变探测结果的一致性验证了采用SBAS-InSAR技术获取研究区域形变结果的准确性，同时可以看出白格滑坡范围在经历2 次崩塌滑动后，区域内表现明显的地表形变特征，最大变形已超过110 mm（图8 紫色区域）。且根据现场勘测的地表错位现象可以看出，滑坡区域的形变处于活跃状态。

图8 白格滑坡形变速率图及现场照片Fig.8 Baige landslide deformation rate map and site photos

降雨量一直被认为是滑坡形成的一个重要因素。雨水不仅可以侵蚀岩体裂隙，同时也增加了滑坡体的负重，从而导致滑坡区域的裂缝形变速率加大。在形变量增加的过程中，滑坡体局部区域产生挤压，进而造成了滑坡体中部和底部的隆起，当前缘位置不足以支撑堆积的岩体时，从而导致滑坡整体垮塌。因此，为进一步预测2 处典型滑坡区域的形变趋势，文中结合隐患点的时序形变特征与其所在时间跨度内的时序降雨量进行了对比分析（图9）。

分别在2 处滑坡范围内选取了3 个特征点（图7 和图8 红色三角形位置），根据特征点在时间跨度内的形变趋势可以看出2 处滑坡区域形变加速时间多发生在月降雨量超过90 mm 的月份。特征点形变特征和江达县境内气象监测站获取了2017年1月—2018年12月的月降雨量（图9），沃达滑坡区域内的特征点发生第1 次形变加速的时间为2018年5月，对应月降雨量112 mm，再次出现形变加速时间为2018年7月和9月，对应的月降雨量分别为290 mm 和190 mm；白格滑坡区域内的特征点发生第一次形变加速的时间为2017年10月，再次出现形变加速的时间分布在2018年5—9月，这一时间段对应江达县2018年的降雨高峰期。结合研究区的降雨量数据及其时序形变特征，发现滑坡区域形变量的大小与降雨量的多少密切相关，降雨次数密集或降雨量大的时间段，灾害区域的形变速率远大于降雨频率少的时间段。进一步采用两处滑坡的现有研究资料[24,26]对文中使用SBAS-InSAR技术所得研究区域形变结果进行的验证，研究资料与现有研究结论相一致，沃达与白格滑坡区域存在复活迹象。加上降雨量与形变量之间的正相关关系，雨季时将进一步增加2 处滑坡区域的形变位移量，进而增加沃达滑坡隐患与白格滑坡发生再次垮塌的风险。

图9 白格与沃达滑坡形变特征和区域内月降水量Fig.9 Baige and Woda landslides deformation characteristics and monthly rainfall in the region

在以上分析2 处典型巨型滑坡特征的基础上，根据研究区内形变量的大小，发现在金沙江下游波罗乡阿当村附近存在2 处隐患区域（图6、图10）。此2 处区域在时间跨度内的最大累计形变量均超过45 mm。在文中选取影像数据的时间跨度内，A 区域的累计形变量达到48 mm，其中该区域最大形变速率−53 mm/a，发生在2017年9月，当月降雨量达到168 mm；B 区域位于A 区域的下游，在研究区域内的整体形变探测结果中，该区域的累计形变量达到53 mm，其时序形变速率在2018年7月达到最大，数值为−45 mm/a，当月降雨量达到290 mm。图10 分别为2 处区域时序形变量达到最大时的分布情况，相比于沃达和白格2 处已知滑坡区域的形变特征，A 和B 区域的累计形变量较小，其形变范围较为分散，没有产生局部剧烈形变，这是其未发生垮塌的原因之一。但随着在降雨影响下形变量的增加，不排除2 处隐患区域发生坡体坠落的可能性，且2 处隐患均处于金沙江干流沿岸，灾害一旦发生必将对金沙江的流通造成极大地影响。

图10 滑坡隐患区域A、B 时序形变Fig.10 Time-series deformation of landslide hazards (A、B)

5 结论

（1）文中采用SBAS-InSAR技术监测西藏江达县金沙江流域沃达和白格典型巨型滑坡形变，结合遥感影像和现场勘查，证明了SBAS-InSAR技术是一种非常有效的滑坡形变监测手段。同时，随着观测次数的增多可以有效降低时空去相关对InSAR技术的影响，能准确恢复滑坡形变历史和分析灾变趋势。

（2）在验证和确保该技术有效的基础上，对西藏江达县内2 处典型的巨型古滑坡形变趋势进行分析，结果表明沃达滑坡中部及前缘位置存在存在明显变形，且部分形变量已超过100 mm；白格滑坡残留体同样存在较大的形变，尤其是边界区域已存在明显的形变漏斗，最大形变量已超过110 mm。同时结合降雨量表明降雨量对滑坡形变速率有非常密切的正相关关系。

（3）在金沙江下游江达县波罗乡阿当村存在2 处隐患区域，形变量均超过45 mm，相比于沃达和白格2 处已知滑坡区域的形变特征，其形变范围较为分散，没有产生局部剧烈形变，这是其未发生垮塌的原因。但应持续进行监测，在有明显灾变迹象时提前开展防治。