中国是地质灾害高发国家，各类地质灾害给人们的生命和财产造成了巨大的损失。在地质灾害监测和滑坡研究中，地表形变量是反映斜坡体当前稳定性及运动状态最直接的物理量^[1]。因此，监测地表形变可以为探测滑坡等地质灾害隐患提供重要信息。

合成孔径雷达干涉测量（interferometric synthetic aperture radar, InSAR）技术是近20年来发展起来的一门广泛应用于监测自然灾害和人类活动导致的地表环境变化的技术。InSAR技术因其全天候、大范围、高精度、强时效、具有穿透性的特点，而在地表形变观测和地形测绘中得到广泛关注和研究^[2]。近年来，国内外学者利用InSAR技术对地质灾害监测进行了大量研究。文献[3]利用时间序列的InSAR技术对四川丹巴地区的滑坡灾害隐患进行了早期识别，并成功识别出17处持续变形的不稳定斜坡^[3]。文献[4]利用Sentinel-1A数据对美国加州纳帕地震断层几何参数进行反演，并得到了其滑动分布。文献[5]综合运用卫星遥感、无人机航拍、InSAR技术对四川茂县叠溪镇新磨村滑坡特征与成因机制进行了研究。文献[6]通过永久散射体测量技术（permanent scatters InSAR，PS-InSAR）技术对意大利卡拉马尼科泰尔梅市的滑坡进行监测，获取了内符合精度达到毫米级的监测结果。文献[7]使用小基线集测量（small baseline subset InSAR，SBAS-InSAR）技术对挪威北部滑坡进行形变监测，并进行滑坡编目。

本文首先从InSAR技术的原理出发，分别简述了差分干涉测量技术（differential InSAR, D-InSAR)、SBAS-InSAR和PS-InSAR的基本原理，并联合使用多种InSAR技术对四川省龙门山-大渡河区域约60 000 km²的区域进行了地表形变探测; 然后通过对InSAR成果的目视解译，为滑坡、泥石流等地质灾害的筛查提供了靶目标；最后对解译出的形变区进行了不少于5%的野外实地验证。

1 利用多种InSAR技术探测地表形变

根据研究区域的位置、获取影像的数量以及不同数据处理方法的技术特点，本文联合使用了D-InSAR、SBAS-InSAR和PS-InSAR进行数据处理，并结合多种InSAR技术测量成果对形变区进行综合解译，以提高形变区识别的准确率。通过联合多种InSAR技术可对绝大多数变形尺度的地质灾害隐患进行探查和监测。

1.1 利用D-InSAR技术探测地表形变

D-InSAR技术是通过去除地形、平地、基线、大气等相位分量来获取形变信息的干涉测量方法。干涉处理后的相位信息${{\varphi }_{\text{int}}}$可表述为：

${{\varphi }_{\text{int}}}={{\varphi }_{\text{def}}}+{{\varphi }_{\text{topo}}}+{{\varphi }_{\text{flat}}}+{{\varphi }_{\text{atm}}}+{{\varphi }_{\text{orbit}}}+{{\varphi }_{\text{noise}}}$

(1)

式中，${{\varphi }_{\text{def}}}$、${{\varphi }_{\text{topo}}}$、${{\varphi }_{\text{flat}}}$、${{\varphi }_{\text{atm}}}$、${{\varphi }_{\text{orbit}}}$、${{\varphi }_{\text{noise}}}$分别为形变相位、地形相位、平地相位、大气相位、轨道相位和噪声相位。本文使用10 m分辨率的数字高程模型（digital elevation model，DEM）模拟${{\varphi }_{\text{topo}}}$后将其去除；${{\varphi }_{\text{flat}}}$、${{\varphi }_{\text{orbit}}}$和${{\varphi }_{\text{noise}}}$均属于系统误差，本文根据卫星飞行状态矢量和设计参数等去除；${{\varphi }_{\text{atm}}}$可通过建立与高程相关的大气模型去除。最终，得到的形变干涉相位${{\varphi }_{\text{def}}}$为：

${\varphi _{{\rm{def}}}} = \frac{{4{\rm{ \mathsf{ π} }}}}{\lambda }{\rm{\Delta }}R $

(2)

式中，λ为L波段的波长；ΔR为雷达视线向形变量，通过三角几何关系可将ΔR转换为水平或垂直方向的形变量。

利用D-InSAR技术仅需2~4幅SAR影像即可实现地表厘米级形变的探测。因此，D-InSAR技术被广泛应用于滑坡、泥石流、地面沉降和塌陷监测，尤其适用于地质灾害隐患的早期筛查和应急调查。但由于D-InSAR技术在实际应用中受时空失相干和大气相位延迟的影响，其在地质灾害监测方面存在较大的误差^[8]。

1.2 利用SBAS-InSAR技术探测地表形变

首先利用SBAS-InSAR技术将SAR影像按时空基线的大小组合成若干个集合，然后利用最小二乘法得到每个小集合的地表形变时间序列，最后利用奇异值分解法将多个小基线集联合起来求解，从而得到整个观测时间内的形变时间序列^[9-10]。

SBAS-InSAR技术可有效降低时空去相干和大气延迟误差，当数据量达到5期时即可获得较为可靠的厘米级监测成果，适用于滑坡、泥石流、地面沉降、地裂缝和地面塌陷的探查和监测^[11]。

1.3 利用PS-InSAR技术探测地表形变

PS-InSAR技术是在D-InSAR技术的基础上，利用相位稳定的点目标降低时空去相干和大气延迟，达到时序地表形变的高精度测量。在没有获取先验永久散射体（permanent scatters，PS）点的情况下，需要进行PS点的初选，一般采用基于振幅离差的倒数法进行^[12]：

$M = \frac{{{m_t}}}{{{\sigma _t}}}$

(3)

式中，M为振幅离差倒数；${m_t}$为时序振幅均值；${\sigma _t}$为时序振幅标准差。选取PS点后，进行DEM误差校正，使得邻近PS点相位差小于π，确保PS点相位值可以正确解缠。对解缠后的相位在时间域进行高通滤波、空间域进行低通滤波，评估大气、轨道、DEM残余相位的误差后，最终获得形变信息。

PS-InSAR技术通过PS点间组网差分有效降低了大气效应对数据处理成果的影响，监测精度可达毫米级，能实现对滑坡、崩塌、泥石流、地面沉降、地裂缝和地面塌陷的精确监测。但在实际应用中，该技术要求数据量至少达到16期才能获得较为可靠的结果，同时在植被覆盖较为茂盛的区域容易因失相干而导致无可用PS点，使有效监测范围大大减少^[13]。

2 研究区域与数据

龙门山断裂带、大渡河沿岸（简称龙门山-大渡河区域）为青藏高原东南缘和横断山脉的一部分，地面海拔479~7 439 m，分为川西北高原和川西山地两部分。川西北高原地势由西向东倾斜，分为丘状高原和高平原，区域内丘谷相间，谷宽丘圆，排列稀疏，广布沼泽。川西山地西北高、东南低，根据切割深浅可分为高山原和高山峡谷区。主要山脉为岷山、巴颜喀拉山、牟尼芒起山、大雪山、雀儿山、沙鲁里山。龙门山-大渡河区域地处青藏高原东缘龙门山逆冲推覆构造带和鲜水河断裂带，受印度板块挤压及四川盆地阻挡作用，形成主体上从向北至向东的右旋构造应力场，内部地块应力作用极为复杂。该区域地层岩性复杂，带内岩层自前震旦纪以来均有发育，大量出露岩浆岩及中-深成变质岩，变形变质强烈，褶皱、节理等相互交错，断层以高角度陡倾为主，劈理极为发育。历史上地震活跃度高，尤其是汶川、芦山等地震叠加影响后，大量山体岩层震裂松动，“内伤”严重。加上区内地形高陡，人工难以全面排查，大量城镇等人口聚居区地处河谷、沟谷沿岸，是地质灾害高风险区域，也是国土部门急需全面摸底、排查地质灾害隐患的重点区域。

本文编程获取了2017-11—2018-10分辨率为3 m的127景ALOS-2影像，完整覆盖了龙门山-大渡河区域地质灾害防治重点16县约60 000 km²的区域，研究区域及影像覆盖范围如图 1所示。ALOS-2卫星由日本空间局于2014年发射，获取L波段SAR影像，重访周期为14 d，主要用于测绘、灾害监测、环境和资源调查等。获取了2016-01—2018-02分辨率为5 m×20 m的96景Sentinel-1历史存档数据。Sentinel-1雷达卫星是欧洲航天局哥白尼计划（Global Monitoring for Environment and Security，GMES）中的地球观测卫星，由两颗搭载C波段合成孔径雷达的卫星组成，重访周期12 d。本文使用的DEM为1:10 000基础测绘成果数据，分辨率为10 m，采用2000国家大地坐标系，1985国家高程基准，3°分带，现势性为2013—2015年。

3 数据处理与误差分析 3.1 数据处理

本文采用GAMMA软件进行数据处理，对获取得到的ALOS-2数据使用D-InSAR和SBAS-InSAR技术进行处理，Sentinel-1数据使用PS-InSAR技术进行处理。21景ALOS-2数据可完整覆盖龙门山-大渡河区域，每个影像覆盖区获取了6~10景不等的时间序列影像。4景Sentinel-1数据可以完整覆盖龙门山-大渡河区域，每个影像覆盖区获取了24景时间序列影像。

利用D-InSAR技术进行数据处理时，选取较新的影像作为主影像，所有影像两两组合，其中，时间基线最长为224 d，空间基线最长为127.5 m。主辅影像配准精度控制在0.1像素以内，地理编码精度控制在0.2像素以内。本文通过建立与高程相关的大气模型去除大气误差，使用最小费用流的方法进行相位解缠。

利用SBAS-InSAR技术进行数据处理时，以综合相关系数为测度选取最佳主影像，所有影像与公共主影像分别组合。其中，时间基线最长为112 d，空间基线最长为111.5 m。

利用PS-InSAR技术进行数据处理时，主影像选取采用了与SBAS-InSAR处理过程相同的选取方法，使用GAMMA软件对主辅影像进行处理共生成160个干涉对，选用最小费用流方法进行相位解缠，并进行时间域的高通滤波、空间域的低通滤波、去除轨道及大气误差等操作。

如图 2所示，在龙门山-大渡河区域，利用D-InSAR技术共处理得到188个干涉对，利用SBAS-InSAR技术在21个影像覆盖区域内分别处理得到形变探测成果，利用PS-InSAR技术在4个影像覆盖区域内分别进行数据处理得到形变探测成果。联合多种InSAR测量成果对龙门山-大渡河区域进行形变区解译：（1）人工目视解译形变区；（2）利用多期D-InSAR成果对解译形变区进行组合验证；（3）若解译形变区所处位置存在SBAS-InSAR成果，则利用SBAS-InSAR成果进行验证；（4）若解译形变区所处位置存在PS-InSAR成果，则利用PS-InSAR成果进行验证；（5）利用三维地图的地形信息和植被信息进行综合判断；（6）确认解译结果。对所有处理成果重复步骤（1）~（6），共解译得到形变区840处，根据D-InSAR、SBAS-InSAR和PS-InSAR处理成果分别解译出620、54和166处形变。

3.2 误差分析

根据野外实地调查反馈结果，约有70%的解译形变区发现形变迹象，其中明显形变约占56%，轻微形变约占14%，验证了联合多种InSAR技术监测地表形变的可靠性。由于龙门山-大渡河区域地处四川省西部，高程悬殊、地貌复杂，受到大气误差、植被覆盖等因素的影响，解译的形变区存在误检的情况。InSAR监测地表形变的主要误差来源于大气误差、时空基线去相干、外部DEM误差等。本文中时空基线均控制在较小的范围内，时间基线最长不超过224 d，空间基线最长不超过127.5 m，故时空基线引起的去相干可以忽略不计。本文采用10 m分辨率的DEM数据作为去除地形相位和地理编码的参考数据，ALOS-2数据经过4×4多视后分辨率约为10 m，Sentinel-1数据分辨率为5 m×20 m，所用数据与DEM配准误差控制在0.2像素左右，故10 m分辨率的DEM数据在本文研究中是可以满足要求的。大气误差包括层流引起的大气延迟相位和由大气湍流引起的大气延迟相位。西南山区地形复杂，本文根据外部DEM采用拟合的方法将与地形相关的大气误差去除，然后再基于干涉图采用叠加和滤波的方式进一步减弱大气效应^[14]。

4 地表形变区时空分析 4.1 形变区空间分布特征分析

形变区空间分布特征将从形变区所处地层结构、形变区所属岩性和形变区所处活动断层3个方面进行研究。

从地层结构的角度看，龙门山-大渡河区域位于华南地层大区，包括巴颜喀拉地层区、扬子地层区和南秦岭-大别山地层区3个Ⅱ级区^[15]。如图 3所示，位于龙门山-大渡河区域西部的形变区属巴颜喀拉地层区之玛多-马尔康分区的金川小区，北部的形变区属南秦岭-大别山地层区之摩天岭分区的九寨沟小区，东部的形变区属扬子地层区之上扬子地层分区的九顶山小区、成都小区，南部的石棉、汉源地区属扬子地层区之康定地层分区。区内地层从前寒武系至第四系均有出露，岩性复杂。

从岩性的角度看，形变区的岩性以坚硬岩组和较软岩组居多，且岩性地质时代跨越较宽，从震旦系至白垩系基本均有分布。如图 3所示，坚硬岩组的形变区主要分布在北川-映秀断裂与茂县-汶川断裂之间及西侧地区，以及平武-青川断裂北侧地区，岩性以元古界、志留系、二叠系-侏罗系的花岗岩、闪长岩及侵入岩脉体等为主。松散岩组的形变区主要分布在龙门山山前冲洪积平原区，岩性以第四系河湖相松散砾石、粘土沉积为主。软弱岩组的形变区分布最少，主要沿龙门山山前主边界断裂呈狭长线性分布，岩性以侏罗系-白垩系泥岩、粉砂岩，以及上新统砾岩与粘土岩互层为主。

从主要活动断裂的构造角度看，龙门山-大渡河区域内断裂构造主要为龙门山构造带。如图 4所示，Ⅰ为西秦岭造山带，其中Ⅰ₁为中秦岭前陆盆地，Ⅰ₂为南秦岭陆缘裂陷盆地，Ⅰ₃为玛沁-略阳结合带，Ⅰ₄为北秦岭后陆逆冲-褶皱带；Ⅱ为松潘-甘孜造山带主体，其中Ⅱ₁为摩天岭逆冲-推覆带，Ⅱ₂为巴颜喀拉-马尔康弧形逆冲-滑脱叠置带，Ⅱ₃为丹巴-汶川弧形逆冲-滑脱叠置带，Ⅱ₄为石渠-雅江弧形逆冲-滑脱叠置带；Ⅲ为上扬子陆块，其中Ⅲ₁为米仓山前陆逆冲-褶皱带，Ⅲ₂为龙门山前陆逆冲-推覆带，Ⅲ₃为四川前陆盆地；F₂为舟曲-成县断裂，F₁₀为荷叶断裂，F₄为玛曲-略阳断裂，F₅为岷山-雪山断裂，F₉为平武-青川断裂，F₆为茂县-汶川断裂，F₈为安县-灌县断裂，F₁₁为北川-映秀断裂，F₁₂为鲜水河断裂，F₁₃为武山-天水断裂，F₁₄为米仓山前缘断裂。解译的形变区在空间分布上有成群、成带的特征，与龙门山构造带的走向相一致。形变区主要集中于该断裂构造内的龙门山后山断裂、龙门山中央断裂、龙门山前山断裂及山前隐伏断裂内。

4.2 典型区域形变分析

为验证联合多种InSAR技术监测地表形变的准确性和一致性，在形变发育最为密集的茂县范围内选取InSAR解译形变区进行对比分析。

4.2.1 形变区概况

选取的InSAR解译形变区位于茂县三龙乡，如图 5所示，解译形变区微地貌表现为缓坡，顶部高程为2 400 m，底部高程为1 810 m。该形变区所处构造为单斜结构，地震烈度Ⅷ度。该形变区位于表层土体为碎石土的斜坡，下覆基岩为千枚岩，坡度为20°~30°，剖面呈梯形，坡向正北向，周边未见有明显基岩出露。如图 5所示，形变斜坡的坡面左侧有茂县长龙黄草坪电站，斜坡下部有河流及道路，坡面中部有居民点及道路，形变区整体坡度约20°~32°。将坡面分为3个区域，a区域位于坡面最左侧，已形成类似于“滑梯”的凹槽，最大坡度为32°；中部b区最大坡度为28°，与左侧a区域衔接处最大坡度为29°；最右侧为c区，最大坡度为27°。

4.2.2 多种InSAR监测结果对比验证

1）D-InSAR处理结果。利用D-InSAR技术对ALOS-2数据进行处理，将时间基线限制在30 d内，共有10个干涉对在此区域有形变反应，如图 6(a)所示。在此期间形变区整体存在持续性变形，局部存在形变差异，形变范围左侧的变形共有7个干涉对有所反应，形变范围的中部与右侧也有不连续形变的反应，但形变反应较微弱。将时间基线限制在60 d内，共有8个干涉对在此区域有形变反应，如图 6(b)所示。在此期间形变区整体变形明显，与周边无变形区有明显差异。其中，有5个干涉对在形变范围的左侧存在明显形变反应，另外3个干涉对在相同区域也有不同程度的形变反应。

2）SBAS-InSAR处理结果。利用SBAS-InSAR技术对ALOS-2数据进行处理，升降轨数据的时间跨度为2017-11—2018-08。将SBAS-InSAR处理得到的沿视线向形变成果数据叠加到强度影像上，如图 7所示。其中，地表形变在升轨方向表现为抬升，在降轨方向表现为下降。分别在升降轨成果上选取a、b和c这3个区域内的点P₁、P₂和P₃，将升降轨数据获取的形变速率和累积形变量进行平均，并生成它们的日形变速率图和累积沉降量曲线，如图 8所示。形变区域中部存在明显形变，在形变区的左侧凹槽处，由于失相干导致点较为稀疏。其中，a区域为形变最大的区域，在P₁点处，累积形变量为-100.5 mm，由日形变速率求得的年形变速率达到126.5 mm/a；区域b形变量比区域a小，P₂点处累积形变量为-47.2 mm，年形变速率为59.4 mm/a；区域c内P₃点处累积形变量为-16.1 mm，年形变速率为20.3 mm/a，是形变斜坡内年形变速率最小处。

3）PS-InSAR处理结果。利用PS-InSAR技术对Sentinel-1数据进行处理，使用数据的时间跨度为2016-01—2018-02。将沿视线向形变PS成果点叠加至10 m分辨率的DEM上，如图 9所示，在此形变范围内共选出PS点482个，形变范围面积为3.7 km²。在a、b和c这3个区域内分别选取点P₁、P₂和P₃（与SBAS-InSAR成果上选取的点相同），并生成它们的日形变速率图和累积沉降量曲线，如图 10所示。其中，a区域为形变最大的区域，在P₁点处，累积形变量为-270.4 mm，由日形变速率求得的年形变速率可达138.9 mm/a；b区域形变量比a区域形变量小，P₂点处累积形变量为-130.7 mm，年形变速率为61.8 mm/a；c区域内P₃点的累积形变量为-51.9 mm，年形变速率为22.6 mm/a，是形变斜坡内年形变速率最小处。分析发现，越靠近左侧凹槽，年形变速率越大，反之则越小，由此推断形变斜坡发生的失衡可能与左侧凹槽相关。

综上所述，SBAS-InSAR和PS-InSAR在P₁、P₂和P₃处形变速率具有一致性，证明了D-InSAR、SBAS-InSAR和PS-InSAR技术获取结果的准确性、可靠性和一致性。

4.2.3 野外核查与验证

通过野外核查发现该形变区呈圈椅状，为一古滑坡堆积体，平均宽约750 m，长约1 200 m，厚约2~5 m，形变区规模约200万m³，形变区内存在局部形变，坡上公路内侧土体垮塌现象明显，中后部发育拉张裂缝。该地质灾害隐患点威胁1户6人，威胁财产约20万元。

如图 11（a）所示，区域a范围内的P₁点位于农户果园内，附近发育一条长度大于10 m的裂缝，呈弧形展布，坡体局部下挫，最大下挫深度约50 cm，坡体表部修建的水渠挫断。根据现场访问，该水渠修建于2017年11月，水渠最大下挫深度约18 cm，平均每月下挫深度约2~3 cm。实地信息勘察结果与D-InSAR、SBAS-InSAR、PS-InSAR成果反映的趋势基本一致。

如图 11（b）所示，区域b范围内的P₂点位于公路外侧，附近发现蓄水池池壁拉裂，裂缝宽约2~5 cm。蓄水池左侧边界可见滑壁，局部整体下滑30~50 cm，实地勘察结果与D-InSAR、SBAS-InSAR、PS-InSAR成果反映的趋势基本一致。

区域c范围内的P₃点附近未发现明显形变。

5 结语

龙门山-大渡河地区是中国地质灾害高发区域，本文联合多种InSAR测量技术以工程实践的方式探测了该地区约60 000 km²范围内的地表形变，并通过野外实地勘察对解译形变区进行了验证。本文介绍了D-InSAR、SBAS-InSAR和PS-InSAR 3种干涉测量技术在形变监测中的应用原理，并利用不同的技术方法对ALOS-2和Sentinel-1数据分别进行处理，共解译形变区840处。其中，根据D-InSAR处理成果解译形变620处，SBAS-InSAR处理成果解译形变54处，PS-InSAR解译成果166处，根据野外验证结果，发现约有70%的解译形变区有形变迹象。本文统计分析了形变区分布与地层结构、岩性和活动断层构造之间的关系，并对典型形变区的不同形变成果进行了对比分析，并结合野外实地调查结果验证了不同InSAR技术的一致性和准确性。本文成功探测了龙门山-大渡河区域的地表形变，为地质灾害隐患的探测和调查提供了靶区，验证了InSAR技术应用于大区域范围内地质灾害监测和预防的发展前景。

致谢: 感谢日本宇宙航空研究开发机构提供ALOS-2数据，欧洲航天局提供Sentinel-1数据。