综合InSAR 技术和多源SAR 数据在滑坡变形监测中的应用—以吉林治新村滑坡为例

于海明，张熠斌，方向辉，徐思瑜，徐誉维，张旭晴

（1.吉林省地质环境监测总站（吉林省地质灾害应急技术指导中心），吉林 长春 130021；2.吉林省航测遥感院，吉林 长春 130061；3.吉林大学地球探测科学与技术学院，吉林 长春 130026）

0 引言

新世纪以来，人类活动对环境和气候的影响逐渐增大，生态环境及各种地质条件也随之发生了较大的改变。中国构造运动频繁，地形结构复杂，地质灾害频发，滑坡、崩塌等地质灾害呈现出分布广、高易发等特点。地质灾害严重威胁群众的生命财产安全，并限制区域经济发展[1-2]。滑坡是斜坡岩石体以及大面积冰雪覆盖体沿着贯通的剪切面所发生的体位滑动现象，由于气候、天气、地形、交通、通讯等因素的影响，很难做到早期预警和提前防范[3]。由于环境限制，传统测量手段如水准测量、GPS 测量难以大范围开展地质灾害监测[4]；光学遥感技术易受气候条件影响，并且难以实现对缓慢变形地质灾害隐患的识别[5-6]。

合成孔径雷达干涉测量技术（interferometricsynthetic aperture radar，InSAR）作为一种新型地表形变监测技术，具有空间分辨高、时间分辨率高、精度高以及大范围空间连续覆盖等优势，已成为地质灾害监测的重要手段之一，被众多学者用于地震[7]、滑坡灾害应急排查[8]、滑坡发育监测[9]、地面沉降[10]等的形变监测研究。法国学者Fruneau 等[11]于1996 年首次将InSAR 技术应用于滑坡监测领域，采用差分雷达干涉技术（differential interferometric synthetic aperture radar，D-InSAR）技术处理ERS-1 数据获取了6 景差分干涉图，监测结果清晰的呈现了滑坡的形变特征，且与实测结果一致。由于不同时期获取的SAR 影像间的相干性较差或失相干，以及大气效应的影响制约了D-InSAR 技术在滑坡形变监测领域的应用，而差分干涉测量短基线集时序分析技术（small baseline subset InSAR,SBAS-InSAR ）[12-13]可实现对目标区域的连续时空监测，意大利学者Ferretti 等[14]提出的永久散射体差分干涉测量技术（persistent scatterer interferometric synthetic aperture radar，PS-InSAR）以及Berardino 等[15]SBAS-InSAR，能够获取地面较高精度的地形数据以及地表时序形变数据，较大程度提升了SAR数据干涉测量的精度。姚佳明等[16]利用SBAS-InSAR技术对煤层开采诱发的地表形变模式进行了研究，证实了SBAS-InSAR 技术监测地表形变的准确性；赵富萌等[17]采用SBAS-InSAR 技术识别了Karakorum 公路沿线滑坡点；Zhang 等[18]利用了传统经验模型与SBASInSAR 监测相结合的方式成功实现了甘肃省中部黄河黑泰河潜在滑坡的分析评价。InSAR 技术已被广泛应用于滑坡地质灾害监测，但单一SAR 数据监测的局限性限制了对滑坡形变的监测效果。C 波段的Sentinel-1A 影像时间基线小，影像数据多，适合应用时序InSAR技术监测滑坡时序形变特征，但在植被覆盖区因相干性过低而获得的监测点稀少；L 波段的ALOS-2 数据时间基线长，影像数量少，适合利用D-InSAR 技术进行滑坡形变特征监测。采用少量长波段ALOS-2 影像和大量短波段Sentinel-1A 影像结合进行滑坡形变监测，能够在具有一定植被覆盖度的山区探测到较为明显的滑坡地表形变[19-20]。

治新村滑坡是吉林市内非常典型的土质滑坡，在2017 年7 月强降雨过程中发生了滑动，滑坡体滑到村民住宅墙根处，造成住宅后墙变形，雨水进入民房内，但未造成人员伤亡。以往对该滑坡的监测，只是传统的人工地面调查，基于InSAR 技术对治新村滑坡进行形变监测的研究却鲜有报道。本文选取2017 年1 月6 日—2017 年12 月8 日的27 景Sentinel-1A 数据，使用SBASInSAR 技术对治新村滑坡及周边区域进行了地表形变监测，并进行形变特征及趋势分析。考虑到滑坡所在区域植被覆盖情况及东北地区季节性积雪对SAR 影像相干性的影响，选取2 景穿透性更强的L 波段ALOS-2 数据监测滑坡体变形特征，综合运用SBAS-InSAR 和DInSAR 技术，以及C 波段的Sentinel-1A 数据和L 波段ALOS-2 数据，以期验证监测结果的可靠性。研究结果对治新村滑坡灾害防治具有重要的指导意义，并可为后期类似滑坡地质灾害监测和研究提供参考依据。

1 研究区概况及数据

1.1 治新村滑坡

治新村滑坡位于吉林省吉林市船营区大绥河镇治新村12 社（图1）。滑坡前缘高程306.8 m，后缘高程322.8 m，长约327 m，宽约48 m，滑坡面积达15 696 m2。平面形态为不规则扇形，滑动面剖面形态为凹形，产状：151°∠33°。土质滑体，平均厚1.2 m。地层岩性为残坡积物砂砾质土，滑床岩性为范家屯组板岩、粉砂岩。滑坡陡坎、后壁发育状况一般，尚可辩认。侧边界不发育，前缘、剪出口发育一般，尚可辩认。拉张裂缝、剪切裂缝较发育，裂缝宽度一般在0.4 m，长1～5 m 不等，规模较小，贯通性较差，见有树木歪斜现象。滑坡体在雨水的浸润作用下，沿基覆界面滑动形成滑坡，为自然滑坡。目前，该滑坡现处于蠕变阶段，不稳定。

图1 研究区位置图Fig.1 Location map of study area

1.2 SAR 影像和DEM 数据

文中用于SBAS-InSAR 的C 波段数据为欧空局于2014 年发射的Sentinel-1A 卫星影像，该卫星轨道高度693 km，重访周期12 d，覆盖范围达到42 500 km2，方位向分辨率为13.98 m，斜距向分辨率为 2.33 m，雷达波长为 5.6 cm，具有条带模式（strip map，SM）、宽幅干涉模式（interferometric wide，IW）、极宽模式（extra-wide swath，EW）和波模式（wave mode）4 种成像模式；用于D-InSAR 处理的L 波段数据采用日本陆地观测卫星ALOS PALSAR-2 影像，L 波段的ALOS PALSAR-2 数据具有较强的穿透力，可以很好地监测地壳运动，获取较准确的数据。

本文选用了C 波段（2017-01-06—2017-12-08）27 景降轨影像（表1），成像模式为IW 模式，极化方式为VV；L 波段（2016-07-26—2017-08-22）2 景影像（表2），成像模式为ScanSAR 扫描模式，极化方式为HH。DEM 数据可用于去除地形相位以及高程误差相位计算，是SAR 数据处理的辅助数据，选用NASA 和NIMA 联合发布的分辨率为90 m 的SRTM3 数据。

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

表2 ALOS-2 影像信息Table 2 ALOS-2 image information

2 研究方法

数据处理使用ENVI SARscape 平台。采用SBASInSAR、D-InSAR 方法分别对C 波段和L 波段SAR 数据进行处理，SBAS-InSAR 对配准好的多幅SAR 影像通过设定的时间基线和空间基线阈值选取合适的干涉对组合，并基于相干目标进行相位解缠和参数解算获取长时间序列的形变信息，相较于D-InSAR，SBAS-InSAR可以分析更大时间序列数据，更好地克服了时空失相干的影响，并减小了大气延迟误差、地形误差、高程误差和其余噪声误差。SBAS-InSAR、D-InSAR 数据处理流程如图2 所示。

图2 SAR 数据处理流程图Fig.2 Workflow of SAR data processing

SBAS-InSAR 是一种基于分布式目标的时间序列分析技术，将时空基线较短的影像两两配对成干涉像对，应用奇异值分解的方法求解单点的形变相位方程，解算高程误差及形变速率，通过残余相位对大气相位和非线性形变进行反演后获得该时间段内的形变时间序列。本文通过设置时间基线、空间基阈值控制Sentinel-1A 数据集生成干涉像对的数量，时间基线长度为 60 d，空间基线阈值为理论空间基线值的45%，C 波段27 景SAR 影像共获取干涉像对98 对，其中主影像日期为 2017年3 月31 日，干涉像对以及基线连接情况如图3 所示。然后基于干涉像对进行SLC 影像配准，生成干涉图，经过Goldstein 滤波处理提高干涉条纹的清晰度。选取近30 个控制点进行轨道精炼和重去平处理，选择automatic refinement 方法消除可能的斜坡相位。历经2 次SBAS 反演，精确估计且去除地形残余相位、大气效应相位，最后结合DEM 数据进行地理编码后获取各期累计形变图和平均形变速率图。

图3 SBAS-InSAR 时空基线图与像对连接图Fig.3 SBAS-InSAR space-time baseline diagram and image pair connection diagram

D-InSAR 技术是对2 景SAR 影像进行差分干涉，同时需要外部的DEM 数据模拟地形相位，进而去除掉地形相位的影响，仅保留形变相位。首先估算2 景ALOS-2 影像基线情况，包括像对的时间基线、空间基线、多普勒频移、相位代表的高程变化值等，作为InSAR 形变监测的前提，然后基于DEM 数据对主辅影像配准和干涉处理，获取干涉相位图，本文选取2016年7 月26 日影像为主影像，2017 年8 月22 日影像为辅影像。采用Goldstein 滤波提高干涉条纹的清晰度，减少由于时空基线引起的失相干的噪声，使用最小费用流（minimum cost flow，MCF）进行相位解缠，解决相位模糊度的问题。在相对稳定区域选择GCP 进行轨道精炼与重去平处理，采用automatic refinement 方法消除可能的斜坡相位。最后将经过绝对校准和解缠的相位结合，并将合成相位通过地理编码转换到制图坐标系统。

3 结果与分析

3.1 SBAS-InSAR 监测结果

由于治新村滑坡低矮植被覆盖和冬季积雪影响，对于C 波段的Sentinel-1A 数据，干涉测量时失相干较严重，获取的监测点主要位于滑坡后缘和滑坡前缘居民区（图4）。SBAS-InSAR 监测结果显示，滑坡后缘斜坡在2017 年发生沉降，而山谷村落地表发生抬升。滑坡后缘斜坡LOS 向平均沉降速率为2.88 mm/a，山谷村落LOS 向平均抬升速率为19.99 mm/a。

图4 SBAS-InSAR 监测结果Fig.4 SBAS-InSAR monitoring results

为进一步获取滑坡后缘斜坡及山谷居民区地表形变特征，对滑坡后缘监测点P1、P2、P3（图4）进行累计形变量时序分析，同时综合分析居民区最大抬升点和平均抬升累计形变特征，分析滑坡对居民区的影响。图5 说明了滑坡后缘3 个监测点（P1、P2、P3）沿雷达视线方向（LOS 向）的累计形变特征。滑坡后缘在2017年1 月6 日—4 月12 日期间处于缓慢抬升阶段，至4 月12 日累计抬升达10.95 mm；4 月12 日—7 月5 日发生剧烈沉降后处于稳定状态；7 月5 日—7 月29 日期间滑坡后缘地表沉降达12.47 mm，此后地表抬升8.21 mm 后处于相对稳定的起伏变化，最终累计地表累计抬升4.81 mm。同时，图4 展示了受滑坡威胁的山谷村落监测点平均累计形变量变化图，山谷村落地表在1 月6 日—4 月24 日期间同滑坡后缘形变趋势相似，但在4 月24 日之后便持续抬升，至12 月8 日平均累计抬升达19.59 mm，最大抬升点累计抬升达36.62 mm。

图5 SBAS-InSAR 监测点时序特征Fig.5 Temporal characteristics of SBAS-InSAR monitoring points

2017 年7 月13 日8 时—14 日8 时，吉林市遭受罕见暴雨天气影响；7 月19 日8 时—21 日8 时，吉林地区再次出现大范围暴雨、大暴雨天气，结合图4 可以看出，治新村滑坡后缘在1 月6 日—7 月5 日期间地表处于相对稳定状态，而在7 月5 日—8 月10 日出现剧烈的沉降抬升变化，7 月5 日—7 月29 日期间的高强度降雨导致了此次滑坡灾害的发生，具体表现为滑坡后缘先沉降后抬升。因此，雨季是治新村滑坡监测与防治的重点时期。

3.2 D-InSAR 监测结果

D-InSAR 监测结果表明，监测期间在治新村滑坡斜坡上存在5 处主要变形体（图6、表3）。1 号变形体位于治新村滑坡汇水区西南侧斜坡顶端，变形体面积达5 318 m2，监测期间最大沉降量为58.6 mm，最小沉降量为29.4 mm，平均沉降量达43.2 mm；2 号变形体位于滑坡汇水区西侧斜坡中上部，变形体面积达17 973 m2，监测期间最大沉降量为50.6 mm，最小沉降量为24.2 mm，平均沉降量达36.4 mm；3 号变形体位于滑坡汇水区东侧斜坡顶端，变形体面积为4 636 m2，监测期间最大沉降量为68.7 mm，最小沉降量为29.8 mm，平均沉降量达43.5 mm；4 号变形体位于滑坡汇水区东侧斜坡顶端，与3 号变形体北侧相接，变形体面积3 043 m2，监测期间最大沉降量为71.8 mm，最小沉降量为31.3 mm，平均沉降量达49.9 mm；5 号变形体位于滑坡汇水区西侧斜坡北部顶端，较接近居民区，变形体面积达11 281 m2，监测期间最大沉降量为47.6 mm，最小沉降量为20.1 mm，平均沉降量达38.3 mm。

表3 斜坡变形体Table 3 Deformation of slope

图6 D-InSAR 监测结果Fig.6 D-InSAR monitoring results

汇水区两侧斜坡变形体分布及变形体沉降量显示，西侧斜坡大面积处于不稳定状态，并向斜坡下端逐步延伸；东侧斜坡虽然顶端变形体形变量较大，但是变形体面积较小，整个东侧斜坡大部分区域由于植被覆盖程度较好，处于相对稳定状态。因此，治新村滑坡威胁主要来源于汇水区植被覆盖程度较差的西侧斜坡。

根据斜坡两侧变形体特征及斜坡地形特征，针对剖线AA′、BB′、CC′、DD′作形变量及高程变化剖面分析，对比二者变化特征（图7）。从图7 可以看出，区域形变累计量与斜坡高程存在明显负相关，即斜坡高处相对沉降量较大，而随着斜坡高程降低，沉降量亦逐渐减小，局部居民区由于滑坡物质累积，地表呈现抬升现象。

图7 治新村滑坡剖面特征Fig.7 Profile characteristics of Zhixincun landslide

4 结论

（1）SBAS-InSAR 技术可以监测滑坡形变演化特征，D-InSAR 技术可以监测滑坡形变体特征，二者结合进行监测可以更全面地反映滑坡时空演化态势。受研究区低矮植被和冬季积雪的影响，通过C 波段Sentinel-1A 影像获取的滑坡体监测点较少，使用L 波段ALOS-2 数据则可以很好的解决影像失相干的问题。

（2）滑坡后缘在监测期间发生了明显形变，该区域在7 月份强降雨过程中发生了先沉降后抬升的形变，期间沉降量达12.47 mm，监测期间平均沉降速率为2.88 mm/a。同时，山谷居民区地表也在降雨后处于持续抬升阶段，至12 月8 日平均累计抬升达19.59 mm，最大抬升点累计抬升达36.62 mm，监测期间平均抬升速率19.99 mm/a。因此，雨季是治新村滑坡灾害防治的重点时期。

（3）基西侧斜坡大面积处于不稳定状态，有向斜坡下端逐步延伸的趋势，最大变形体面积17 973 m2，平均沉降量36.4 mm；东侧斜坡变形体面积小于西侧，但变形体形变量较大，平均沉降量49.9 mm。结合实际调查，东侧斜坡大部分区域由于植被覆盖程度较好，处于相对稳定状态。因此，治新村滑坡灾害的主要威胁来源于汇水区植被覆盖程度较差的西侧斜坡，且形变同地形相关明显，斜坡顶端沉降量几乎是最大的，随着斜坡高度降低，沉降量亦逐渐降低，但山谷居民区由于滑坡物质累积则呈现地表抬升。