滑坡因发生频率高、隐蔽性强、分布广、造成损失严重而受到高度重视^[1]。四川省茂县核桃坪地区受汶川-茂县断裂影响，褶皱及断裂等不良构造发育，同时由于新构造运动与强降雨作用，斜坡变形迹象突出，威胁人民生命财产安全^[2]，迫切需要对该地区的关键斜坡进行地质灾害的长期监测预警。

合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术可以获取地表cm级甚至mm级形变，并且可实现全天候、不受光照和气候等因素影响持续观测，已经成为监测滑坡变形不可或缺的技术手段^[3]。SBAS-InSAR技术通过将数据分割为若干小基线子集，获得空间基线差分干涉图，从而降低空间去相干影响，同时结合其他SAR数据集计算滑坡表面形变速率，得到空间高密度的形变量变化结果^[4-5]。

利用离散元颗粒流方法(particle flow code, PFC)对边坡进行数值模拟，可分析滑坡体破坏形式和评价预测其稳定程度。周健等^[6]利用PFC方法，以边坡土体性质为对照进行数值模拟，结果表明边坡土体随着颗粒粘性的减小，其破坏形式从脆性逐渐向塑性演化。对边坡进行离散元颗粒流数值模拟，分析滑坡从萌芽、发展到破坏的运动过程和破坏机理，已成为边坡防治及空间预测的重要手段。

已有大量研究利用单一的遥感手段或数值模拟方法对滑坡灾害进行分析，而将两者进行有效结合的研究较少。本文以构造运动复杂的四川省茂县核桃坪滑坡为研究区，使用2019~2021年覆盖研究区的62景Sentinel-1A影像，基于SBAS技术获得滑坡区变形发展特征，并结合离散元数值模拟方法预测边坡的变形破坏趋势和破坏后运动过程。同时，对滑坡开展针对性的现场复核工作，以检验SBAS-InSAR技术和PFC数值模拟结果的准确性。

1 研究区概况

茂县核桃坪滑坡由汶川地震诱发，且一直处于活跃状态，其位于岷江支沟永和沟左岸，属渭门乡核桃村下组，坐标103°50′48″E、31°45′28″N(图 1(a))，属于典型的山区峡谷地貌和高原型季风气候区，年平均气温变化幅度约11.1 ℃，年平均降水量484.1 mm，主要集中在4~9月，最大日降水量75.2 mm。滑坡区相对高差约80 m，坡度28°，主滑方向334°，具阶梯状的纵向剖面，滑坡长165 m，宽195 m，滑坡面积为3.3×10⁴ m²，平均厚度6.1 m，总体积约为18.6×10⁴ m³，属于中型土质滑坡(图 1(b))。

研究区出露的地层岩性主要包括震旦系(Z)、志留系茂县群(Smx)、泥盆系(D)、三叠系侏倭组(T₃zh)和杂谷脑组(T₂z)白云岩、石英砂岩、千枚岩和板岩等(图 1(c))。滑坡区岩层产状为348°∠55°，岩层倾向与山体坡向趋于一致。滑坡坡脚堆积物主要为后缘的老崩坡积物。

核桃坪滑坡所在区域属于龙门山华夏系构造带，发育多条断层，其中茂汶断裂对滑坡起控制作用。2008年汶川地震在研究区地震烈度为Ⅶ度，严重破坏坡体上房屋等基础设施，滑坡后缘产生多条拉裂缝，坡体东部出现小型滑塌现象，地震已导致滑坡进入局部失稳状态，并在此后反复活动。为避免滑坡再遇强降雨、地震等灾害失稳垮塌，故对核桃坪滑坡稳定性进行监测与分析。

2 InSAR变形特征分析 2.1 数据源及处理过程

本文使用欧空局(European space agency, ESA) Sentinel-1A卫星影像数据，数据集的空间覆盖范围如图 2所示。所用影像属于128轨道，涵盖2019-02-26~2021-02-27期间62景数据图像，轨道方向为升轨，地面分辨率为5 m×20 m，幅宽为250 km×250 km，成像模式为IW宽幅模式，极化方式为VV，并采用30 m分辨率的ALOS World 3D DEM以减少地形相位产生的影响。

使用ENVI软件SARscape模块的SBAS算法完成Sentinel-1A 62景升轨数据的干涉处理。具体流程为：1)连接图生成，选择超级主影像时间为2020-03-28，以单视复数(single look complex, SLC)图像为基准进行配准后，根据设定的时空基线阈值进行配对；2)将配对的影像进行干涉处理，完成去平、滤波、相位解缠等一系列工作，去除相干性较低的数据^[7-8]；3)在无相位跃变和干涉条纹的位置选择地面控制点(ground control point, GCP)并结合DEM，对干涉结果进行轨道精化与重去平，优化相位；4)经过残余相位估算和形变速率估算可去除大气相位得到完整的形变速率；5)对上述结果进行地理编码，将SAR坐标下的时序形变信息转换到WGS-84大地坐标系，获取视线向地表形变点的形变速率^[9]。具体处理流程如图 3所示。

2.2 InSAR变形特征

茂县核桃坪滑坡整体植被覆盖密度较小，仅在坡体中后部植被相对茂盛，SBAS处理结果中相干点较连续，影像失相干现象较少，干涉效果明显，可观测到明显的形变速率。采用上述SAR数据集获得核桃坪滑坡2019-02-26~2021-02-27期间雷达视线向(LOS)形变速率，结果如图 4所示。当坡体表面形变点朝向卫星传感器方向运动时形变速率为正值，当形变点背离卫星传感器方向运动时形变速率为负值^[10]。InSAR处理结果显示，核桃坪滑坡后缘视线向形变量为负值，表明滑坡减损体沿滑动面下滑，并且为整个滑坡变形最强烈的区域(图 4红色虚线框)，最大形变速率为-78.4 mm/a。滑坡中部(图 4黄色区域)形变速率为-34.7~-5.5 mm/a。滑坡前缘位置(图 4蓝色虚线框)形变速率出现正值，为滑坡的滑覆堆积体(加积体)，形变速率为9.1~30.1 mm/a。结果表明，滑坡整体向坡脚永和沟方向滑动，滑坡减损体表现为向下沉降，滑覆堆积体因崩滑体堆积、推挤而相对抬升，符合滑坡变形破坏的运动规律。

2.3 InSAR结果验证

为验证InSAR数据的准确性，开展多次实地工程地质调查进行复核，现场坡体的局部变形范围及平均变形速率与InSAR监测结果高度吻合，表明监测数据准确性较高。实地踏勘结果表明，核桃坪滑坡体呈扇面型，滑坡顶端滑壁明显(图 5)，主滑方向为334°。坡体表面耕地分布较广，发育阶梯状陡坎，高度数米不等，延伸可达几十米到数百米。坡体还分布少量冲沟及大量剪切裂缝和拉裂缝，根据延伸方向不同可将裂缝分为3组，其中两组剪切裂缝分别位于滑坡两侧(图 5(b)、(c))，走向约为345°和331°，另一组为坡顶的不规则拉裂缝。整个坡体出现3个次级滑动变形体H1~H3(图 5(a))，每个次级滑动体的分界较为明显。坡体前缘由于受雨水冲蚀和临空面等作用，发育少量冲沟及崩滑变形体，并且可见明显的树木倾倒现象，树龄大多为10 a左右，说明核桃坪滑坡处于持续变形状态(图 5(g)、(h))。

3 离散元模拟 3.1 数值模型建立

采用PFC2D程序对滑坡变形和位移进行计算分析，对滑坡变形破坏的危险区域、滑动过程进行预测和模拟。为更好地分析滑坡不同部位变形运动特征，通过观察滑坡的外部形态，可将滑坡分为滑坡前缘、滑坡中部、滑坡后缘3个部分，根据斜坡剖面图(图 1(d))建立概化地质模型(图 6)，模型底边宽180 m，左侧边界高112 m，右侧边界高44 m。颗粒最小半径为0.03 m，最大半径为0.06 m，模型共生成颗粒19 008个。

采用PFC程序模拟边坡时，颗粒细观特性参数的选取是能否成功模拟的关键^[11]。由于细观参数与岩石宏观力学参数之间无直接的对应转换关系，需要通过数值实验，如单轴压缩实验、双轴压缩实验等模拟得出应力应变曲线，不断调整参数使数值实验结果与边坡岩土体的变形模量、泊松比、强度参数等结果一致，从而标定出合适的细观参数。本文通过资料收集和综合类比确定岩土体的宏观力学参数(表 1)，采用单轴压缩实验标定岩土体的细观特性参数(表 2)。

3.2 模拟结果

本文旨在通过将岩体强度参数折减直到边坡失稳，观察边坡的变形破坏特征，模拟采用的岩体参数无法完全真实地反映边坡实际岩体力学性质，模型的坡体边界线只作为坡体变形的参考。通过调研确定颗粒流参数折减方法，即同时折减切向粘结强度、法向粘结强度和摩擦系数^[12]。

滑坡在发展演化过程中坡体要经历孕育阶段-蠕动阶段-变形阶段-滑动阶段，各发展阶段对应不同的变形速度，根据InSAR监测的坡体变形速率可确定PFC对应的模拟阶段，因此可将二者结合对比分析。通过分析SBAS-InSAR和PFC结果可知，两者的坡体变形范围基本一致，图 7(b)为运行至80 000步时边坡X方向位移图，从图中可以看出，坡体已发生明显变形，最大位移量达到9 m。同时可以观察到边坡各部位的变形特征各不相同，变形体上部较缓区域以地面开裂为主；坡面陡缓转换处岩层产生强错动变形，同时伴随有地面开裂现象；下部较陡区域表现为浅表层塌滑和错动，与现场实地调查相吻合。

当运行至120 000步时，滑坡剪切面已经被贯通，坡体发生整体滑动，坡脚土体滑出数十米(图 8)，变形体H1对应位置堆积物明显减少，呈下降趋势，而变形体H3对应区域有堆积现象，呈隆起趋势。为研究斜坡滑动过程并预测坡体不同部位的发展演化情况，在PFC建模过程中对坡体颗粒进行标记，布置控制监测点，位置如图 6所示。通过各测量点的位移-时程曲线(图 9)可知，坡体在90 000步之前曲线斜率较小且基本保持一致，表明该阶段坡体因地震及降雨等因素变形启动后保持匀速变形，由InSAR监测结果可知目前坡体处于该阶段，且位于坡体中部的C和D监测点的变形速率明显大于位于前缘和后缘的A和D点，与InSAR解译结果一致。若遇到强降雨或地震，坡体极有可能进入加速变形阶段，该阶段位移加速度增大，由于坡体运动过程存在不协调变形，导致各部位位移量差距增大，促使后缘产生更多拉张裂隙，而中、前缘移动后形成临空面，为岩土体进一步失稳提供变形空间，进而坡体将出现从下到上的牵引式破坏。

4 变形破坏趋势分析 4.1 滑坡发展演化机制

核桃坪滑坡的强烈变形区分布在坡体减损带和滑覆堆积区，滑坡两翼则以缓慢匀速变形为主，滑坡现场实地调查形变特征与SBAS-InSAR显示的形变速率结果(图 4)及PFC数值模拟结果(图 7)具有高度一致性。为分析降雨对滑坡变形速率的影响，选取滑坡减损体形变速率较大的形变点A，得到其时序形变曲线，与2019-01~2021-02研究区降雨量进行对比分析(图 10)。结果表明，在2019-06~08、2020-06~09集中降雨期，A点变形速率增大，累积形变量突增；降雨相对较少时，A点变形速率变缓，因此A点形变轨迹与该区降水强度存在良好响应，说明降雨是诱发核桃坪滑坡变形乃至最终破坏的主要原因。

核桃坪滑坡属于中型顺层土质滑坡，坡顶发育明显拉裂缝，其LOS向形变速率特征及PFC模拟结果表明，核桃坪滑坡属于蠕滑-拉裂变形破坏机制，边坡的似均质体特性以及顺层结构为蠕滑变形提供了有利条件。坡体前、中、后三部分变形速率差异明显，说明导致滑坡发育台坎的原因为坡体的显著差异性形变。结合现场地质调查和InSAR形变速率结果进行分析，滑坡减损带沉降明显，加积带存在抬升现象，因滑坡下伏千枚岩与上覆土体的物质特性存在较大差异，因此滑坡潜在滑动面为第四系松散崩坡积物与千枚岩岩层的土岩结合面，该面以上为剪切蠕变带。滑坡变形破坏进程主要受滑动面控制，当上覆堆积体的下滑力超过滑动面的实际抗剪阻力时，坡体将失稳迅速滑落。当前滑坡正以匀速沿滑动面向冲沟及坡前临空面方向滑动，破坏形式表现为坡体局部滑塌和前缘散落堆积。

4.2 滑坡局部变形分析

对3处次级滑动变形体(H1~H3)进行细致调查发现，H1因坡脚土体崩落临空面变形使抗滑力降低导致滑动变形；H2周围可见冲沟，表明该滑动变形体为集中降雨引起，范围相对较小；H3因变形体前部临空面变形导致出现牵引式滑动变形，且降雨和人类工程活动因素促使变形体下滑，该变形体体量最大，可反映滑坡整体变形趋势，因此对滑动变形体H3进行具体分析。

H3滑动变形体位于滑坡右翼，滑动方向与坡向一致，变形体上可见次滑壁和新鲜崩滑物，表明该变形体正在发生快速滑动。结合InSAR形变速率分析H3变形体的滑动特征，取滑动变形体主滑方向上前、中、后1^#~5^#监测点，绘制其2019~2021年滑坡监测点变形-时间曲线(图 11)。由图可知，位于变形体顶部的1^#监测点累积形变量为负值，表明该处发生沉降，且形变速率和累积形变量均最大；2^#监测点位于滑坡中部，仅发生小幅沉降，形变速率在0附近波动；3^#~5^#监测点均位于滑坡前缘滑覆堆积区，累积形变量均为正值，其中3^#、4^#监测点形变趋势基本一致，但4^#监测点形变量更大，表明该位置堆积作用更强烈。5^#监测点在2020-04由前期的快速抬升转为沉降，说明该位置的前缘堆积体在快速加积后，向下沉降速率大于堆积速率，表现为向下发生缓慢的滑动变形，这也可以解释滑坡前缘不断延伸的现象，与野外实地调查结果基本吻合。

4.3 趋势分析

将PFC模拟结果(图 7(a))和SBAS年均变形速率(图 7(b))进行对比可知，PFC模拟得到的边坡变形特征与InSAR监测结果具有一致性，变形最大值均位于滑坡中部，并向两侧递减。基于以上工作，可进一步预测核桃坪滑坡的变形破坏趋势，在一定条件下，坡体变形将经历3个阶段：1)坡体最初表现为表面蠕滑，上覆土层向坡下弯曲变形，坡体后缘产生拉应力；2)后缘拉裂，可形成面向坡下的台阶，也是坡体当前所处阶段；3)随着剪切变形进一步发展，潜在剪切面发生剪切扰动，中部剪应力集中部位可被扰动扩容，从而使斜坡下半部分逐渐隆起。随着变形体发生滑动，后缘明显下沉，拉裂面从初期的张开变为渐趋闭合，变形进入破坏阶段，一旦潜在剪切面被间断贯通，则发生滑坡。

尽管本文结合InSAR和离散元数值模拟来预测边坡的变形破坏趋势，由于深部岩土体的变形破坏程度、降雨量、区域构造作用无法精确掌握，理想化建模固然与现实存在差异，因此预测结果的可靠性仍然需要时间验证以及其他案例佐证。

5 结语

本文结合SBAS-InSAR测量技术与PFC离散元颗粒流方法，对四川省茂县核桃坪滑坡形变特征进行分析，并模拟滑坡未来变形破坏趋势，经野外实地调查分析，得到以下结论：

1) InSAR处理结果表明，核桃坪滑坡2019-02~2021-02坡体后缘视线向形变速率为负值，表明滑坡减损体沿滑动面下滑，最大形变速率为-78.4 mm/a，滑坡中部最大形变速率为-34.7 mm/a，滑坡前缘滑覆堆积区形变速率为正值，最大形变速率为30.1 mm/a。综合地形因素以及现场地质调查分析认为，核桃坪滑坡现阶段主滑方向为334°，破坏形式为坡体局部滑塌和前缘散落堆积，由于滑坡中、后缘崩坡积物的推挤作用，使得坡体表现为沿滑动面向下的匀速滑移运动。

2) PFC模拟结果表明，核桃坪滑坡中部台坎位置变形程度最大，向前后两侧递减，与InSAR监测结果一致。通过对模拟过程中坡体不同部位位移进行监测，发现若坡体进入加速变形阶段，其后缘将产生更多拉张裂隙，而中、前缘移动后形成的临空面可为岩土体进一步失稳提供变形空间，进而坡体将出现从下到上的牵引式破坏。

3) 预测核桃坪滑坡的失稳趋势可分为3个阶段，分别为初期的表层蠕滑，坡体后缘产生拉应力；中期的后缘拉裂，形成台坎；并预测后期滑坡受潜在剪切面扰动作用，当剪切面被间断贯通，则发生滑坡。