基于DS-InSAR技术的金沙江流域贡觉地区滑坡与地裂缝形变特征

盛 磊，张 露，杜玉玲，邓文杰，闫世勇*

(1. 中国矿业大学 自然资源部国土环境与灾害监测重点实验室,江苏 徐州 221116; 2. 中国矿业大学 环境与测绘学院,江苏 徐州 221116; 3. 中国科学院空天信息创新研究院, 北京 100094)

0 引 言

地处中国西部山区的金沙江流域，地形复杂，岩体破碎，再加上金沙江断裂带等构造活动影响，极易导致滑坡和地裂缝等地质灾害，严重威胁周边地表建筑安全，对沿线基础工程施工与运营安全形成了严重威胁。水准测量、GPS测量和三维激光扫描技术等传统的高精度地表形变监测手段，难以在植被覆盖茂密、人迹罕至的高山峡谷区展开对滑坡等地质灾害的大范围监测。光学遥感技术虽然凭借其范围广、低成本等优势被广泛应用于地质灾害隐患排查，但是易受气候条件影响，且对缓慢变形的地质灾害隐患难以识别。

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)技术，作为一种广域厘米级形变测量技术，为地震、火山、滑坡、地裂缝等地质灾害监测提供了一种有效手段。但由于其受到时间、空间失相干和大气延迟等因素影响，难以提取毫米级地表变形，所以后期发展出了永久式散射体差分干涉测量(Persistent Scatterers Interferometry,PSI)技术和小基线集差分干涉测量(Small Baseline Subset Interferometry,SBAS)技术，有效弥补了传统InSAR技术的不足。但相比城市地区，金沙江流域属于高山峡谷区，该区域由于缺少人工建筑物等永久式散射体点，所以PSI技术仅能获取少量测量点；SBAS技术易受浓密植被与大气延迟影响导致失相干现象，无法识别到足够的高相干点目标，导致难以全面准确地反映区域内滑坡与地裂缝形变信息。而且，金沙江流域高差巨大，地形复杂，会导致强烈的大气垂直分层效应，使得传统的时空滤波方法在该区域难以有效去除大气延迟误差。而融合分布式散射体差分干涉测量(Distributed Scatterers InSAR,DS-InSAR)技术通过识别稀疏植被、裸土等分布式散射体，将分布式散射体与永久式散射体融合后采用传统时序方法进行解算并获取形变信息。DS-InSAR技术在高山峡谷区地质灾害形变监测中，能够有效提高低相干区域的相干性，减弱大气延迟误差，减小相位梯度对相位解缠的影响，有助于提升形变监测结果的可靠性。

本文基于欧洲航天局Sentinel-1升降轨数据，采用DS-InSAR技术，利用InSAR通用型大气改正服务(Generic Atmospheric Correction Online Service for InSAR,GACOS)补偿大气延迟相位，获取金沙江流域贡觉地区2019年2月至2021年7月高精度地表形变分布图，并在此基础上结合地形和降雨等数据分析了重点区域滑坡与地裂缝形变特征，为该区域地质灾害防治以及沿线交通工程安全建设与运营提供参考。

1 研究区概况与数据来源

研究区位于西藏自治区昌都市贡觉县沙东乡(30.30°N，98.91°E)至雄松乡(30.50°N，98.91°E)内的金沙江流域(图1)。该区域内海拔最高为5 054 m，最低为2 587 m，最大高差为2 476 m，受金沙江水流侵蚀以及SN向金沙江断裂带的影响，地形复杂，孕育多处滑坡和地裂缝等地质灾害。

图1 金沙江流域贡觉地区遥感影像Fig.1 Remote Sensing Images of Gongjue Area in Jinsha River Basin

考虑到研究区内地形起伏，阴影、叠掩等几何畸变情况严重[图1(c)、(d)]，本文联合Sentinel-1A升轨与降轨数据进行时序InSAR处理。笔者选用覆盖研究区的64景升轨和69景降轨IW模式VV极化数据(表1)，分别以2020年3月14日和2020年3月21日影像为参考影像，采取多主影像策略构建差分干涉图，相应时空基线如图2所示；采用分辨率为90 m的航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission,SRTM)DEM作为参考DEM，用于在差分干涉时去除地形相位。InSAR通用型大气改正服务用于补偿大气延迟相位。

表1 SAR数据集参数和时空基线组合情况

图2 升轨和降轨时空基线图Fig.2 Spatio-temporal Baselines of Ascending and Descending Orbits

2 分析方法

2.1 DS-InSAR形变提取

DS-InSAR技术核心思想是通过识别研究区各分辨单元的统计同质像元，并采取一定方法进行相位优化，将优化质量高的像元作为分布式散射体融入传统时序InSAR技术，实现低相干区的高精度地表形变监测。

本文采取多主影像策略构建差分干涉图，将振幅离差指数(DA)低于阈值(0.6)的像元确定为高相干点候选点，用于与后续得到的分布式散射体候选点融合。基于配准后SAR影像振幅信息，采取FaSHPS方法识别各像元周围一定大小矩形窗口(15×15)内的同质像元(SHPs)数量，将数量高于20的像素作为分布式目标初步候选点；之后在各初步候选点上建立相干矩阵，采用特征值分解方法选取最大特征值对应的特征向量作为优化相位，并计算分布式散射体候选点时间相干性，将时间相干性高于阈值(本文Sentinel-1升轨和降轨阈值分别为0.65和0.75)的像元确定为最终的分布式散射体候选点。

分布式散射体候选点时间相干性计算公式为

(1)

将分布式散射体候选点与高相干候选点融合处理，并进行相位稳定性分析，将时间相干系数高于阈值(Sentinel-1升轨和降轨阈值分别为0.7和0.8)的像元作为最终的高相干点。

时间相干系数计算公式为

(2)

对缠绕相位去除残余DEM误差后，进行3D相位解缠得到高相干点解缠相位值，此时相位中往往包含较为严重的大气延迟误差。可基于InSAR通用型大气改正服务计算出各干涉对大气延迟相位并从解缠相位中去除，补偿大气延迟误差。接着，在高相干点建立关于相位变化速率与解缠相位的观测方程。由于采用多主影像策略导致方程秩亏，所以通过奇异值分解(SVD)求解方程广义最小二乘解，经过相位与形变转换后得到监测时段内研究区形变速率和时序形变量。图3为本文DS-InSAR技术数据处理流程。

图3 DS-InSAR数据处理流程Fig.3 Processing Flow of DS-InSAR Data

2.2 坡度向形变转换

由于视线向形变速率难以直观反映滑坡真实运动状态，所以将研究区视线向形变投影至各点最大坡度方向，获取滑坡坡度向形变速率。图4为Sentinel-1升轨雷达视线向(LOS)和坡度向几何关系图，其中P点为滑坡上一点，为卫星飞行方向角(方位向与北方向夹角)，为卫星雷达入射角，为滑坡方位角，为P点最大坡度角，为视坡夹角(视线向与坡度向夹角)。假设滑坡运动沿坡度向向下为正，根据式(3)将雷达视线向形变速率投影得到坡度向形变速率，其中和分别代表雷达视线向和坡度向形变速率。雷达视线向形变速率为正值，表示目标沿着靠近卫星方向运动；雷达视线向形变速率为负值，表示目标沿着远离卫星方向运动。

图4 Sentinel-1雷达视线向和坡度向几何关系Fig.4 Geometric Relationship Between LOS Direction of Sentinel-1 and Slope Direction

3 结果分析与讨论

3.1 整体形变特征

基于DS-InSAR技术和相应的Sentinel-1升降轨数据，分别获取了金沙江流域贡觉地区地表年平均形变速率监测结果(图5)。由图5可见，升轨和降轨数据共监测到22处明显地质灾害隐患，D2、D5、D6、D8、D9、D16区域为非重叠区域。其中，升轨数据共监测到17处明显变形区域，雷达视线向最大形变速率为-254 mm·年；降轨数据共监测到19处明显变形区域，雷达视线向最大形变速率为-215 mm·年。监测结果显示，沉降最严重的区域位于沙东乡(D1区域)，其沉降速率和沉降区域范围最大，发生地质灾害的可能性较大。对比图5(a)和(b)可以看出，升轨和降轨监测结果均能获取到足够数量和密度的监测点，较好地刻画了滑坡沉降中心和沉降范围。从图5(a)可以看出，金沙江东岸的升轨数据监测结果受到阴影、叠掩等几何畸变影响，导致监测点稀少，丢失有效形变信息，如D2、D8、D9、D16等区域，而这些区域在图5(b)中形变明显。由于升轨与降轨数据对不同坡度和坡向滑坡的形变敏感度不同，所以D1、D5和D11等区域在升降轨结果中也表现为不同的形变特征。上述分析表明，结合升轨与降轨数据不仅可以有效避免单一轨道因几何畸变影响导致监测盲区，而且可以弥补单一轨道方向对部分滑坡形变不敏感的不足，从而降低地质灾害隐患监测漏检率，提高监测结果准确性。为进一步验证形变结果的可靠性和分析地质灾害特点，有必要针对典型滑坡与地裂缝形变信息开展详细分析。

(3)

3.2 典型区域形变特征

3.2.1 沙东乡沙东滑坡

图5 DS-InSAR年平均形变速率图Fig.5 Images of DS-InSAR Annual Average Deformation Rate

沙东滑坡位于沙东乡格果村，地处金沙江西岸，滑坡体总面积约为10.3 km，整体坡度为20°～40°，坡向为南向与北向(图5中D1区域)。沙东滑坡光学影像如图6(a)所示，滑坡后缘高约3 790 m，存在明显陡坎；滑坡中部为地势平缓的平台区，为格果村行政驻地，局部区域存在裂缝；滑坡前缘高约2 680 m，受金沙江水流冲刷侵蚀影响，变形迹象明显。由图6(b)可知，滑坡变形体整体雷达视线向平均形变速率超过20 mm·年，最大形变速率为-254 mm·年。结合光学影像目视解译和形变分布特征，可将沙东滑坡分为北部S1和南部S2两个子变形体[图6(b)]。变形体S1最大形变速率为-254 mm·年，强变形区位于变形体中部至前缘部分，后缘存在两个较强变形区；变形体S2强变形区域位于坡体中部，最大形变速率为-226 mm·年，坡体后缘形变速率整体上略大于前缘。图6(c)为该滑坡沿坡度向平均形变速率图，坡度向形变速率最大为279.6 mm·年。

H为高程；图(c)坡度向形变速率为正值，表示滑坡沿最大坡度方向向下滑动；黑色箭头表示目标坡度向，其长短与形变速率成正比图6 沙东滑坡光学影像与形变速率图Fig.6 Optical and Deformation Rate Images of Shadong Landslide

沿变形体S1和S2主要滑动方向分别绘制沙东滑坡后缘至前缘的剖面线a—a′和b—b′[图6(a)]，得到剖面线平均形变速率[图7(a)和(b)]。从图7(a)可以看出：变形体S1在坡体后缘存在一处较强变形区域，位于剖面线0.38 km处，形变速率最大为-75.2 mm·年；平台区形变速率较小，整体速率处于-50～-30 mm·年之间；坡体前缘存在两处强变形区域，分别位于剖面线2.0 km和2.4 km处，相应最大形变速率分别为-217和-151 mm·年。从图7(b)可以看出：变形体S2形变速率沿剖面线向坡脚方向先逐渐增大，在坡体中部形变速率达到最大值，约为-162 mm·年，随后形变速率减少，在坡体前缘形变速率出现两处局部最大值，分别是-144 mm·年(剖面线1.61 km处)和-92.1 mm·年(剖面线1.86 km处)，在剖面线2.24 km处形变速率由负变正，之后呈现抬升趋势，原因可能是坡体内部蠕滑导致坡脚部分发生反翘，有待实地进一步验证。

图7 沙东滑坡形变特征Fig.7 Deformation Characteristics of Shadong Landslide

图7(c)给出了变形体S1和S2强变形区域中心A点和B点累积形变量曲线及相应时段的月降雨量数据。由图7(c)可知，A点和B点雷达视线向累积形变量分别为-217.9 mm和-405.5 mm，形变速率与降雨量有一定相关性。月降雨量较小时，如2019年11月至2020年4月以及2020年11月至2021年4月，A点和B点形变速率较小；月降雨量较大时，如2019年5月至10月、2020年5月至10月以及2021年2月至6月，A点和B点形变速率较大。随着2021年5月降雨量相比往年同期降雨量显著增加，A点和B点形变速率也呈现出加速趋势。因此，沙东滑坡变形受降雨因素影响较为显著，滑坡整体处于强变形状态，存在失稳可能，应对该滑坡加以重点关注。

3.2.2 敏都乡果巴滑坡

果巴滑坡位于敏都乡果巴村，地处金沙江西岸，滑坡体总面积约1.54 km，整体坡度处于10°至40°之间，坡向为东(图5中D7区域)。果巴滑坡光学影像如图8(a)所示，滑坡后缘高约3 470 m，主要为村间道路，坡体中部为果巴村居民地和耕地；坡体前缘高约2 650 m，发育一条长约340 m的地裂缝，坡脚受金沙江水流冲刷侵蚀，碎屑垮落，变形迹象明显。

图8(b)为果巴滑坡降轨雷达视线向平均形变速率图，最大形变速率为61.9 mm·年。监测结果显示，果巴滑坡雷达视线向形变速率由坡体后缘至坡体前缘逐渐增大，坡体后缘以及坡体中部南侧居民区无明显形变，坡体中部北侧存在较强形变，整体强变形区位于坡体前缘南部。果巴滑坡坡度向平均形变速率如图8(c)所示，最大形变速率为187.4 mm·年，强变形区域位于坡体前缘南部，同时坡体前缘北部(即地裂缝北侧)也存在较强形变迹象，对附近道路安全存在一定威胁。

图8 果巴滑坡光学影像与形变速率图Fig.8 Optical and Deformation Rate Images of Guoba Landslide

由图9(a)主滑方向雷达视线向形变速率剖面线可知：果巴滑坡滑动速率沿剖面线随距离增加而逐渐增大；在剖面线0.6 km处，坡体形变速率减小；而后形变速率继续增加，在剖面线1.8 km处达到最大值(56.4 mm·年)，随后又有所减小。由果巴滑坡强变形区域中心点A处的雷达视线向时序累积形变量[图9(b)]可以看出，2019年2月至2021年7月A点累积形变量整体呈逐渐增加趋势，最大累积形变量为112.8 mm，而相应时段降雨数据则与果巴滑坡形变不存在显著的相关性。果巴滑坡形变可能主要受金沙江水流常年冲刷侵蚀影响，且滑坡中部旱季农田灌溉渗水等人为活动因素对滑坡形变也存在一定影响。

图9 果巴滑坡形变特征Fig.9 Deformation Characteristics of Guoba Landslide

3.2.3 雄松乡上缺所地裂缝群

雄松乡上缺所地裂缝群位于雄松乡上缺所村北侧(图5中D12区域)，受金沙江断裂带构造活动影响，发育超过10条长短不一的地裂缝，对区域内房屋、道路以及农田造成了不同程度的破坏。雄松乡上缺所地裂缝群光学影像如图10(a)所示，走向近EW向，西起上缺所村西北侧，东至金沙江，该区域内主要发育有3条较大地裂缝G1、G2、G3，地裂缝最大长度可达1.4 km，最大宽度可达10 m。

图10 上缺所地裂缝群光学影像与形变速率图Fig.10 Optical and Deformation Rate Images of Shangquesuo Ground Fissure Group

图10(b)显示了雄松乡上缺所地裂缝群雷达视线向平均形变速率分布图，最大值为-178.8 mm·年，形变速率沿地裂缝走向逐渐增加。由图10(b)可知：上缺所地裂缝群西部区域较为稳定，无明显形变；中部和东部区域形变严重，强变形区域大致呈椭圆状，长轴与地裂缝群主走向平行。

图11 上缺所地裂缝群剖面线累积形变量Fig.11 Cumulative Deformations of Profile Lines of Shangquesuo Ground Fissure Group

沿平行于地裂缝走向绘制剖面线a—a′，沿垂直于地裂缝走向等间隔绘制剖面线b—b′、c—c′和d—d′，得到上缺所地裂缝群沿各剖面线累积形变量(图11)。从图11可以看出，雄松乡上缺所地裂缝群形变量变化方向与地裂缝走向较为一致。上缺所地裂缝群累积形变量沿剖面线a—a′先增加后减少，形变量最大位置位于剖面线0.72 km处，累积形变量达-307.64 mm；剖面线b—b′累积形变量最大值为-248.31 mm，位于地裂缝G1处，说明上缺所地裂缝群在该剖面线方向的形变主要受地裂缝G1影响，地裂缝G2、G3影响较弱；剖面线c—c′累积形变量最大值位于地裂缝G1与G2之间，为-314.74 mm，两地裂缝之间累积沉降量均超过-300 mm，且两侧形变曲线呈现较好的对称性，说明该剖面线方向形变主要受地裂缝G1、G2共同影响，且二者影响程度大致相当；剖面线d—d′累积形变量最大值为-292.94 mm，位于地裂缝G3处，同时在地裂缝G1处存在累积形变量较大值(-241.45 mm)，说明该剖面线方向形变受地裂缝G3和G1影响，并且前者影响最大。综上所述，雄松乡上缺所地裂缝群形变变化方向与地裂缝走向一致，不同地裂缝的不同部分对两侧形变影响程度有所不同，但形变量最大位置均位于地裂缝两侧。

上缺所地裂缝群形变中心A点时序累计形变量与该区域月降雨量如图12所示。当月降雨量较大时，A点形变速率较小，如2019年5月至9月以及2020年7月至2020年10月；月降雨量较小时，A点形变速率较大，如2019年10月至2020年6月以及2020年11月至2021年5月。由此可见，雄松乡上缺所地裂缝群形变对降雨量的变化响应复杂，且在时间上存在一定的滞后性。

图12 上缺所地裂缝群沉降中心时间序列累积形变量Fig.12 Time Series Cumulative Deformation of Subsidence Center of Shangquesuo Ground Fissure Group

4 结 语

本文采用融合分布式目标的DS-InSAR技术，利用Sentinel-1升降轨数据提取了金沙江流域贡觉地区滑坡与地裂缝形变速率和时间序列累积形变量信息，共探测到22处明显地质灾害隐患，其中沙东乡附近滑坡雷达视线向最大形变速率可达到-254 mm·年；在此基础上，进一步结合降雨数据分析了重点滑坡与地裂缝形变特征。

(1)融合分布式目标的DS-InSAR技术应用于同期升降轨数据，不仅能够克服研究区高相干点稀少导致的低相干影响，而且减弱了高山峡谷区叠掩和阴影等几何畸变影响，能够全面准确地获取研究区滑坡和地裂缝等地质灾害隐患的形变信息，为进一步开展滑坡和地裂缝等灾害分析提供了良好的技术方法和数据支撑。

(2)金沙江流域内滑坡形变受降雨影响程度不同，强降雨会加速部分滑坡形变，易导致滑坡失稳造成堵江，如沙东滑坡等，因此，十分有必要对该区域内滑坡开展持续监测。雄松乡上缺所地裂缝群形变速率主要沿地裂缝走向变化，地表不同位置的形变受地裂缝影响程度有所不同，但最大形变位置均靠近地裂缝两侧，且地裂缝群形变速率变化相对于月降雨量变化存在滞后性。

(3)本文仅获取和分析了一维雷达视线向以及坡度向形变信息，后续将考虑联合多个轨道和区域地形特征获取多维形变信息，并将综合利用高分辨率光学影像解译与现场调查进一步分析，为该地区滑坡和地裂缝等地质灾害普查与监测提供依据。