0.   引言

2008年5月12日汶川Ms 8.0级特大地震后效应显著，震后滑坡与崩塌频发，产生了大量松散固体物质，在极端气候触发下震后泥石流灾害活跃，震后泥石流物源主要分为同震崩滑堆积体及强震作用形成的震裂山体物源^[1～3]。2010年文家沟“8•13”泥石流中，强暴雨冲刷同震滑坡碎屑流堆积体形成高容重泥石流，造成6人死亡^[4]；2013年汶川县桃关沟震后崩滑物源在连续降雨中受水流侧蚀影响持续坍滑补给，最终暴发“7•10”泥石流^[5]；2019年汶川县登溪沟同震崩滑坡堆积物源在降雨汇流影响下最终起动，暴发“8•20”泥石流^[6]。蟹子沟同震形成的震裂山体分别于2008、2010年雨季发生滑塌，最终暴发泥石流并形成堰塞湖^[7]。

地震会导致坡体变形并长期处于蠕变状态，随着震后风化和降雨侵蚀影响，一旦被暴雨激发极易作为物源参与泥石流起动。震后坡体的滑动或蠕变可借助干涉合成孔径雷达技术（InSAR）进行毫米级精度全天候地表形变监测^[8]。目前InSAR技术凭借其高时空分辨率、广域覆盖能力及毫米级监测精度，已成为地质灾害研究的核心手段^[9]。在灾害预警研究方面，首先通过长周期时序监测捕捉地表形变的动态演化特征^[10]，进而基于空间形变场匹配，对比地表形变区与InSAR监测结果的空间对应关系，实现老滑坡危险区快速圈定与新发育潜在不稳定区的识别^{[11 − 12]}。在震后灾害研究方面，结合质量守恒原理与InSAR技术，可推演滑坡体形变深度并揭示地震前后形变模式的突变特征^[13]，最终通过形变深度结合形变速率正演滑体滑面运动状态，并建立灾害演化模型，可为预警阈值设定提供定量依据^[14]。针对同震崩滑体作为震区泥石流主要物源的特征，InSAR技术可通过物源形变监测，快速锁定泥石流暴发的物源起动诱因^[15]，进而基于反演深度确定此类物源储量。因此准确识别震裂物源位置及变形量对震后泥石流动力学成灾机制具有重要意义。

本文以唐家山滑坡区大水沟泥石流为例，在现场调查查明泥石流发育特征基础上，基于SBAS-InSAR技术，采用Sentinel-1卫星遥感数据对其流域内的物源分布进行综合遥感识别，对识别的物源区进行时序形变监测，然后以物源识别结果为依据估算可起动物源区储量，综合判定大水沟可起动物源存在的可能性，为大水沟上游拟建唐家山水库风险规避提供参考。

1.   研究区概况

大水沟泥石流位于四川省绵阳市北川羌族自治县曲山镇，通口河右岸，原唐家山滑坡左侧边界（图1）。泥石流流域面积0.68 km²，主沟长1.37 km，最大高程1613 m，最小高程709 m，相对高差904 m，平均纵坡降581‰。

图  1  研究区概况图

Figure  1.  Location map of the study area

下载: 全尺寸图片 幻灯片

1.1   地质环境条件

北川县位于青藏高原东缘向四川盆地过渡地带，中部为侵蚀构造中山地貌。研究区地处侵蚀构造中山的东南边缘，地形陡峻，相对高差超900 m，受区域地质环境因素影响，区内以暴雨型泥石流为主。大水沟流域地势南高北低，沟道呈“V”型，谷坡25°～45°，沟道整体顺直，发育有两条支沟，流域分区界线不明显，根据流域特征将其划分为形成区、流通区和堆积区三部分（图1）。

研究区处于松潘—甘孜造山带前缘冲断带，呈现大规模、多层级的叠瓦状冲断推覆构造特征，地质构造复杂。地表主要出露基岩为寒武系下统清平组（Є1c）硅质岩、砂岩、泥灰岩、泥岩，第四系堆积物以冲积、坡积物为主，主要分布于河床、两岸坡顶和坡脚、小型冲沟及局部地形较缓部位。

1.2   气象条件

研究区属于山地亚热带湿润季风气候，四季分明，雨量充沛。区域内多年平均气温15.6 °C，最高气温36.1 °C，最低气温−4.5 °C，年温差较大。研究区多年平均降雨量1355.4 mm，降雨主要集中于5—10月，占全年降雨的91.1%，历年一日最大降雨量323.4 mm。

大水沟属中山地貌，地形坡度较大，支沟发育，利于汇水，加之震后沟内松散物源丰富，大量降雨使松散物源稳定性降低，造成该地泥石流频发，每年雨季大水沟均会暴发山洪或泥石流。

1.3   “5•12”地震震后地质灾害发育特征

2008年“5•12”汶川Ms 8.0级地震导致大水沟流域内形成同震崩塌、滑坡等大量松散物源，同时唐家山滑坡后壁产生裂缝，形成震裂山体。震后大水沟于同年6月10日遭遇5年一遇暴雨并暴发泥石流，造成通口河堵塞近4 h（图3a），同年9月24日再次遭遇100年一遇暴雨，暴发更大规模泥石流，造成通口河堵塞17 h，水位抬升近10 m^[16]。

图  3  大水沟历史暴发泥石流图片

Figure  3.  Historical Debris Flow events in the Dashuigou debris flow

下载: 全尺寸图片 幻灯片

地震导致原唐家山滑坡后壁发育多条裂缝（图2），包括LF1走向355°，延伸长度约70 m；LF1-1走向30°，延伸长度约80 m；LF2走向100°，延伸长度约100 m；LF3走向120°，延伸长度约200 m；LF4走向50°～80°，延伸长度约250 m。顶部震裂山体多年持续下错，裂缝不断延展，且震裂山体表层滑塌现象发育，流域内松散物质每逢雨季即会随水流冲出，形成山洪或泥石流，大水沟2013年、2018年、2020年和2021年均冲出不等规模泥石流，导致道路损坏，造成财产损失（图3b）。

图  2  大水沟沟域上游右岸地震震裂裂缝分布位置

Figure  2.  Distribution of earthquake-induced fractures on the right bank of the upper Dashuigou catchment

下载: 全尺寸图片 幻灯片

2.   基于SBAS-InSAR的可起动物源识别

2.1   SAR数据来源

SBAS-InSAR技术是将影像进行配准后，通过设置适当的时空基线阈值，生成差分干涉图后相位解缠，并利用奇异值分解法（SVD）获取地表形变时间序列最小二乘解，可以减少时空失相干和大气相位误差，从而提高地表变形测量的精度。

本文InSAR数据主要选取欧空局（ESA）所提供重访周期12天的Sentinel-1雷达卫星2018-2024年C波段单视复数影像（SLC）作为影像数据，其主要采用垂直极化（VV）方式，影像数据参数见表1。由于大水沟流域较小且流域内大部分坡体朝向N～NW，升轨对该区域敏感度较低，因此本文采用234景Sentinel-1降轨SLC影像数据进行分析。为避免时空失相干问题导致最终生成干涉对的效果较差，需要设置较小的时空基线。因Sentinel-1回访周期为12天，研究区地形起伏较大，变化较快，故设置时间基线阈值为36天、空间基线阈值为200 m，时空基线图见图4。

表  1  研究区SAR影像数据参数

Table  1.  Parameters of SAR image in the study area

数据源	Sentinel-1
影像获取时间	20180109-20241028
波段/波长/（cm）	C/5.6
空间分辨率/（m）	2.7×22
重访周期/（d）	12
极化方式	VV
影像数量	234
数据类型	IW
轨道号	62

下载: 导出CSV 

| 显示表格

图  4  SBAS-InSAR时空基线分布图

Figure  4.  Space-time baseline distribution of SBAS-InSAR

下载: 全尺寸图片 幻灯片

2.2   数据处理

通过HyP3在线平台^[17]和ESIS^[18]对覆盖研究区的Sentinel-1A 降轨数据进行SBAS-InSAR处理，处理具体流程如下，InSAR技术流程图见图5：

图  5  SBAS-InSAR技术流程图

Figure  5.  Processing flow chart of SBAS-InSAR data

下载: 全尺寸图片 幻灯片

（1）通过HyP3平台对Sentinel-1影像数据进行处理。HyP3平台依托云端服务，集成多种InSAR算法，可快速处理大量影像数据并生成差分干涉图，具有使用方便、占用内存小、运算时间短的优点^[19]。该平台提供20×4、10×2和5×1三种多视，本文选择分辨率最高的5×1多视进行处理，设置完毕上传处理任务。该任务主要基于ISCE2进行处理，处理流程如下：首先以第一幅影像为主影像完成影像配准后，进行5×1多视处理以优化影像信噪比，生成初始相位差分干涉图；随后采用Goldstein-Werner功率谱滤波器进行相位噪声抑制，为兼顾噪声抑制与细节平衡，滤波系数设置为0.5，同时结合水掩膜消除水体区域相位突变；进一步利用内置SNAPHU的MCF算法解缠干涉相位，并基于哥白尼GLO-30 DEM及ESA精密轨道数据同步校正椭球相位偏差与地形误差；最终对解缠相位进行UTM/WGS84地理编码，输出可供下载的小基线差分干涉图。

（2）通过ESIS对HyP3在线平台处理完成的相位干涉图进行时序处理。ESIS基于MinPy算法进行时序处理，该算法提出创新方法改正相位解缠误差，因具有计算效率快、测量精度高、可靠性强的优点被广泛应用于地表形变监测^[20]。由于该步骤在本地处理，为提高处理效率和结果精度，故需要对每幅差分干涉图的研究区域进行裁剪对齐。主要处理过程如下：首先在保证平均相干性更高的前提下保留尽可能多的干涉图，选取相干系数阈值为0.6，优选高相干性干涉图；其次由于研究区为深切河谷地貌，对流层大气延迟主要受地形影响，故使用高度相关性方法校正对流层延迟；然后以线性方式去除大尺度趋势误差，提升数据质量；最终在去除上述相位误差的基础上，选取相干性较高的位置作为参考点，开展时间序列分析。

（3）识别验证泥石流物源位置及形变量：根据卫星参数与坡度坡向空间几何关系转换将雷达视线LOS向的形变速度$V_{LOS}$转换成沿坡向和法向的形变速度^[21]：

式中：$V_{LOS}$——雷达视线方向形变速度/（mm·a⁻¹）；

$V_{slope}$——沿坡向形变速度/（mm·a⁻¹）；

$V_{norm}$——沿法向形变速度/（mm·a⁻¹）；

$\phi$——斜坡坡度/（°）；

$\theta$——入射角/（°）；

$\alpha$——斜坡坡向/（°）；

${\alpha }_i$——轨道方向与正北方向的夹角/（°）。

为避免$LS$和$LN$值趋近于0导致速度趋于无穷大，当$LS$和$LN$处于−0.3～0.3区间时，根据正负分别设置为−0.3或0.3^[22]。为减少外部因素对结果的影响并直观展现变形发生区域，根据Sentinel-2卫星提供的植被分布数据，结合研究区的正射影像，确定二维形变转换后的沿坡向和法向形变速度$V_{slope}$和$V_{norm}$中低可视区和植被茂密区的边界，在ArcGIS中采用掩膜去除。考虑到物源变形通常沿坡向下的原因，本文去除其中速度为正的形变点，将可视化分析结果与现场调查发现的松散物源位置对比，确定物源位置。坡体变形主要受重力影响，假设监测期间坡体密度不变，沿坡体走向的变形忽略不计，则可根据质量守恒法推算通过坡向速度与法向速度计算坡体变形厚度^[23]：

式中：$V_{norm}$——沿法向的形变速度/（mm·a⁻¹）；

$f$——流变系数；

$\mathrm{\Delta }x$——沿坡向数据采样间距/m；

$V_{slope}$——沿坡向的形变速度/（mm·a⁻¹）；

$h$——变形深度/m。

对物源区进行时间序列形变监测，综合判定物源变形失稳起动可能性，最终实现对大水沟泥石流沟可起动物源的识别分析。

3.   结果与分析

根据上述SBAS-InSAR流程处理得到大水沟2018—2024年LOS向、沿坡向和法向年平均地表形变速度结果（图6）。LOS向形变速度场表示该区域2018年1月至2024年10月期间地表形变速度在−158.58～62.87 mm·a⁻¹之间。其中大水沟流域内形变速度在−120.28～48.41 mm·a⁻¹之间，二维形变转换后沿坡向和法向形变区域基本与LOS向形变区域一致，变形主要集中于原唐家山滑坡后壁侧、大水沟流域中上游及下游两侧斜坡零星滑塌处（图6中1～3区域），沟内最大变形速度出现在原唐家山滑坡后壁裂缝LF1处，沿坡向和法向形变速度分别达到−226.28 mm·a⁻¹和−159.97 mm·a⁻¹，其余裂缝位置变形也均较大，与现场调查的结果一致。

图  6  大水沟流域可起动物源SBAS-InSAR监测结果

Figure  6.  SBAS-InSAR monitoring results of potential debris flow sources in the Dashuigou catchment

下载: 全尺寸图片 幻灯片

大水沟流域可起动物源识别结果及变形厚度分布见图7。本次共识别到12处发生变形的物源区，其中除3号和9号物源区外，其余物源区坡表均已出现不同程度的滑塌，以6号物源区滑塌最为严重，11号、12号物源区坡表也存在明显滑塌现象。10号、11号物源区内出露基岩破碎现象显著，并且都分布有大面积坡面物源（图8）。

图  7  大水沟流域可起动物源识别与厚度估算结果

Figure  7.  Identification and thickness estimation of potential debris flow sources in the Dashuigou catchment:

下载: 全尺寸图片 幻灯片

图  8  物源区坡表滑塌现象（正射影像）

Figure  8.  Slope surface slumping in source areas （orthophoto imagery）

下载: 全尺寸图片 幻灯片

物源区厚度箱线图见图9，再以95%置信区间估算其极端条件下可能释放的最大潜在可起动物源量，即完全失稳时的理论最大值，结果见表2。其中最大厚度18 m分布于8号物源区，估算储量38.73×10⁴ m³，为原唐家山滑坡后壁残留震裂山体，整体临空面朝向通口河，若沿裂缝发生大面积失稳则直接滑向通口河，仅沿大水沟侧的表层滑塌作为松散物源参与泥石流活动。位于分水岭附近的1、2、5号物源区，物源厚度相对较薄，仅发生浅层变形，估算储量分布在0.15×10⁴ m³～0.39×10⁴ m³不等。其余沟道附近的物源区变形厚度变化较大，大部分集中在1～7 m左右，估算储量分布在0.12×10⁴ m³～6.06×10⁴ m³不等，滑塌现象显著，是泥石流物源的主要补给区。估算误差来自多个方面：（1）InSAR技术的时空分辨率限制导致微地形突变识别能力不足，造成局部厚度估算值存在偏差，且暴雨期滑塌速率超出监测时效易引发滞后误差；（2）MintPy虽已消除大尺度、大气及解缠误差，但残留误差仍会影响形变反演精度；（3）仅采用降轨单轨道数据进行计算亦会引入额外误差。

图  9  估算变形厚度箱线图

Figure  9.  Box plot of estimated deformation thickness

下载: 全尺寸图片 幻灯片

表  2  各可起动物源区厚度及总储量计算结果

Table  2.  Estimated deformation thickness and total volume of each potential source area

物源区编号	厚度范围/m	物源区面积/m2	物源总储量/×104m3
1	1～2	2,192.2	0.39（0.26，0.52）
2	1～3	968.3	0.15（0.05，0.26）
3	1～2	1,128.6	0.12（0.06，0.18）
4	2～7	8,910.2	3.56（2.54，4.58）
5	1～2	2,115.3	0.32（0.13，0.51）
6	2～7	9,083.2	4.09（2.68，5.48）
7	2～7	14,084.9	6.06（4.90，7.22）
8	5～18	43,031.8	38.73（35.83，41.63）
9	1～4	6668.3	1.78（1.41，2.14）
10	2～6	11,317.2	3.40（2.55，4.24）
11	1～7	13,703.9	5.48（4.42，6.54）
12	1～5	4,154.4	1.25（0.83，1.68）
合计		117,358.3	65.33（55.66，74.98）

下载: 导出CSV 

| 显示表格

为进一步验证估算结果的准确性，将InSAR反演结果与图7c所示物探测线数据进行对比分析（图10）。结果显示，估算的物源变形厚度与物探电性分界在整体趋势上呈现出良好一致性，表明估算结果总体合理。图中存在两处局部差异：（1）在测线320−370 m段，InSAR估算厚度未反映该处电性异常，通过正射影像发现该区域存在浅层冲沟，推测差异原因在于暴雨冲刷导致冲沟下切形成破碎带充水，其电性异常可被物探识别，但受限于InSAR的时空分辨率，难以捕捉此类局部突变特征；（2）测线410 m以后区域为6号物源滑塌堆积区，表层快速滑动导致InSAR出现失相干现象，造成厚度估算值异常波动，而现场调查与物探则通过电阻异常可清晰揭示该区域的堆积特征，对比可知InSAR失相干会导致已有崩滑堆积体的厚度估算误差增大，估算值相比实际值偏小。

图  10  InSAR估算厚度与物探结果对比

Figure  10.  Comparison of InSAR thickness estimates with geophysical results

下载: 全尺寸图片 幻灯片

考虑到泥石流物源起动通常与降雨密切相关，而研究区年降雨量高达1355.4 mm，雨季暴雨频发，是导致大水沟流域内坡表滑塌的重要原因之一。为探究降雨事件对物源稳定性的影响，对6、8、9、10四处物源区分别选取1个特征点（PS1～PS4）进行时序形变监测，并将时序监测点LOS向累积位移-时间曲线与唐家山地区降雨量进行对比分析，结果见图11。结果表明，物源形变位移量与降雨量呈明显正相关，2018年极端暴雨事件中，各PS点变形量显著增加，最大达180 mm左右；2020年受极端暴雨影响再次发生强烈变形，并在之后两年间持续发生形变，至2022年雨季，原唐家山滑坡后壁裂缝处监测点PS1累积位移量已达到500 mm左右，其他各PS点位移量也达300 mm左右。经现场证实，大水沟泥石流在2018、2020和2021年都有固体物质冲出，造成道路中断，但未对河道产生显著影响。

图  11  典型PS点LOS向累积形变量

Figure  11.  Cumulative LOS deformation of typical PS points

下载: 全尺寸图片 幻灯片

除已识别的物源区外，大水沟沟道中广泛分布有沟道物源（见图12a、b），下游沟道内堆积有工程弃渣（见图12c），堆渣处因平硐坍塌产生的崩滑物源（见图12d），在极端暴雨条件下，这些未成功识别到的物源亦存在起动的可能，经野外调查确定该部分物源总储量约21.98×10⁴ m³，与前述InSAR技术识别结果结合可得大水沟流域内可起动物源总储量为87.31×10⁴ m³。若极端暴雨状况下起动，泥石流不仅会沿沟道冲入通口河道，还会由于下游沟道狭窄而向两侧漫流，增加对下游侧坡道的侵蚀，最终增加水库的泄洪压力，对水库的施工运营产生不利影响，因此需要在工程前对其进行治理。

图  12  未识别到的物源及弃渣

Figure  12.  Unidentified debris flow sources and waste deposits

下载: 全尺寸图片 幻灯片

上述InSAR技术对部分物源的识别并未取得良好的效果，即使下游沟道依然呈现显著堆积特征，但却难以将沟道内物源与周围特征点作出明显区分。这种物源识别的显著局限性主要来源于几个方面：（1）地形效应干扰：下游沟道岸坡陡峻，局部近直立，导致雷达侧视成像出现几何畸变，可能出现叠掩与阴影现象，显著增加相位解缠误差；（2）人类活动干扰：人工弃渣与工程引起的崩滑物源缺乏前期形变积累，人类高频扰动、松散异质性及短时突发特性引发时空失相干，破坏了干涉相位连续性；（3）植被散射差异：上游因人类活动植被以灌丛居多，雷达波可轻易监测地表形变，但中下游植被以大型灌木为主，植被冠层对微波信号的衰减作用会显著降低地表真实形变的监测精度^[24]。（4）暴雨事件时效性限制：极端暴雨诱发上游物源快速起动形成山洪/泥石流时，沟道剧烈侵蚀过程超出InSAR监测的时间分辨率阈值。这些局限性易使某些区域的物源受相位梯度减弱、形变信号微弱等因素影响，造成物源空间边界漏检现象，尤其在沟谷深切及地形陡峻区域泥石流物质的剧烈活动会进一步降低InSAR识别敏感度。

因此，为提高物源识别体系的可靠性，本文融合了实地调查及无人机航摄构建三维地形约束模型，并联合光学遥感综合开展物源识别。此外在以后的研究中，针对InSAR识别局限性问题，可考虑将高分光学影像与GNSS地表移动数据结合，通过多源数据空间配准与特征级融合，增强物源识别能力，同时结合无人机LiDAR点云数据重建沟道三维结构，可以修正因InSAR几何畸变导致的物源空间定位偏差，从而更精准预测灾害发展趋势，在此基础上进行地物分类提升物源的空间辨识度，重点探索多平台数据的协同解译方法，结合降雨过程中物源稳定性分析不同降雨条件引发的物源起动量差异，实现InSAR监测体系从“形变监测”向“灾害预警”的功能跃升，最终构建“空-天-地”一体化泥石流监测预警系统。

4.   结论

（1）SBAS-InSAR共识别12处具有典型变形特征的震裂山体崩滑和坡面物源，分布于原唐家山滑坡后壁侧、流域中上游冲沟发育处及下游零星滑塌处。

（2）大水沟流域内不同物源区的厚度及储量差异显著，8号物源区厚度最大达到18 m左右，物源储量38.73×10⁴ m³；其余物源区变形厚度多集中于1～8 m左右，储量分布在0.12×10⁴ m³～6.06×10⁴ m³不等。针对可识别的物源，已有崩滑堆积松散物源估算的误差更大，相比实际值会偏小。

（3）极端暴雨事件对大水沟物源形变的增加具有显著影响，受两次极端暴雨作用，原唐家山滑坡后壁位移量达500mm，丰富的物源起动补给，表示大水沟是一条典型暴雨型泥石流沟。

（4）结合InSAR监测与野外调查结果，最终确定大水沟流域内可起动物源总储量为87.31×10⁴m³。大水沟流域内存在丰富的物源补给，一旦起动可能导致河道堵塞，对拟建水库供水工程构成威胁。