联合InSAR和SSA的膨胀土边坡形变特征分析——以南水北调工程为例

朱 武，窦 昊，殷那政，程意清，张双成，张 勤

1. 长安大学地质工程与测绘学院，陕西 西安 710054； 2. 西部矿产资源与地质工程教育部重点实验室 陕西 西安 710054； 3. 生态地质与灾害防控自然资源部重点实验室，陕西 西安 710054

膨胀土是一种高塑性黏性土，土体内含有较多的蒙脱石等黏土矿物，具有很强的亲水性。同时，土体的含水量变化会引起体积变化，具体表现为吸水剧烈膨胀、失水剧烈收缩及反复胀缩形变。除胀缩性外，膨胀土还具有裂隙性和超固结性两个重要特性，次生裂隙的产生与超固结作用都会破坏土体的稳定性[1]。南水北调中线工程于2014年12月正式通水，从丹江口水库调水，为沿线的郑州、石家庄、北京、天津等14座城市提供生产和生活用水，是一项缓解我国华北平原水资源短缺、优化水资源配置的重大战略工程[2]。南水北调中线总干渠约有1/3段穿越膨胀土(岩)区域，干渠通水后其周边地区土壤内的含水量必然升高，而潮湿态的膨胀土在经历降水和蒸发时其体积和强度都会发生变化[3]。特别地，如果在裂隙发育密集的范围内有强降水渗入，裂隙周边的岩土体也会吸水膨胀而逐渐向周边扩散，使渠道表层岩土体达到饱和状态，抗剪强度降低，在表层引发膨胀土浅层滑坡。因此，膨胀土的存在增加了南水北调工程沿线边坡失稳的可能性，为了避免膨胀土渠段滑坡对南水北调工程通水产生影响，对膨胀土区域进行形变监测和特征分析是非常有必要的。

GNSS和应变计等地表观测技术是目前膨胀土形变监测的主要方法，但是其只能获取离散的单点观测结果，对于大面积未知区域的形变探测和监测则存在局限性。合成孔径雷达干涉测量技术(synthetic aperture radar interferometry，InSAR)作为雷达遥感的一个重要分支，具有低成本、大范围、高精度等优势，能够全天时全天候对地表进行观测，已经成为滑坡识别[4]、地面沉降[5]、DEM重建[6]以及地震变形监测[7]的有效手段。目前，已有一些学者利用InSAR技术开展了膨胀土形变监测方面的研究。文献[8]基于TerraSAR-X数据获取了荷兰某农田短期内的膨胀土形变，并定性分析了黏土储水量与地表位移的线性关系。文献[9]采用PS-InSAR方法获取了荷兰某膨胀土堤坝长达6 a的时间序列形变，并结合外部数据与PS点的形变序列详细分析了降水、气温及土壤类型对堤坝形变的影响。文献[10]利用多时相InSAR技术获取了法国某试验基地的膨胀土垂向形变，并结合SMOS卫星的土壤表面湿度数据估计了黏土层深度和厚度的变化。国内应用InSAR技术监测膨胀土形变的研究相对较少。文献[11]利用PS-InSAR技术获取了南宁市区的地表形变信息，重点关注下伏基岩以第三系膨胀土为主的地铁修建区，研究表明地铁沿线沉降均在安全范围内。文献[3]结合南水北调中线工程辉县段膨胀土在运营期间存在的问题，分析了PS-InSAR技术在该渠段膨胀土形变监测中存在的难点和重点。文献[12]应用多时相InSAR技术对南水北调中线地表形变进行研究，并结合降水和土壤湿度变化分析了辉县市段膨胀岩土的形变机理。

尽管已有学者利用InSAR技术开展了膨胀土的形变监测研究，多数针对施工过程中回填的平坦地表，较少开展膨胀土边坡的研究，且未对胀缩形变进行有效分离。基于此，本文收集了2017年3月12日—2021年4月8日的112景Sentinel-1升轨SAR数据，采用多时域InSAR技术获取了南水北调中线某一膨胀土渠道边坡的时间序列形变，并通过奇异谱分析(singular spectrum analysis,SSA)分解得到形变序列的周期项、趋势项和噪声项，开展膨胀土形变与区域内气温和降水之间的关系研究，解析膨胀土的形变特征，为南水北调中线干渠的边坡稳定性评价提供科学依据。

1 研究区域和数据

1.1 研究区域

研究区域位于南水北调中线工程某一膨胀土边坡，该区域地处暖温带，属大陆性季风气候，四季分明，降水充沛，光照充足，夏季降水最多，冬季降水最少，每年最高气温出现在7月，最低气温出现在1月。该渠段地处黄河以南，相对南水北调中线其他渠段而言，区域内降水较多、地下水位较高，且地质条件复杂，沿渠存在大量膨胀土[13]。图1展示了研究区域的地理位置、地形变化及膨胀土分布。图1中标识的P1—P2段为高填方和高挖方膨胀土区域。图1(b)为美国SRTM数字高程模型与德国TanDEM-X数字高程模型之间的差值，反映了该区域2001—2015年的地形变化。可以看到，南水北调渠道大部分为挖方区域，最大挖深约40 m。同时，部分穿越河流的渠道采取填充处理，最大填充高度约为20 m，如图1(b)中的P2点附近。图1(c)展示了渠道沿线P1—P2的地质剖面，可以看到整个沿线分布着厚度不等的膨胀土，其中大部分为强膨胀土，这样的条件增加了南水北调工程沿线边坡失稳的可能性[14]。

图1 研究区域DEM变化及P1—P2点地质剖面Fig.1 Data coverage, DEM changes and geological profile along P1—P2 over the study area

1.2 数据收集

为了分析研究区的时空形变信息，获取了112景哨兵1号(Sentinel-1)升轨SAR数据。其成像模式为IW单视复数(single look complex，SLC)影像，成像方式为TOPS(terrain observation by progressive scans)，时间为2017年3月12日—2021年4月8日，具体参数见表1。同时，为了分析降水、气温与InSAR形变之间的关系，从河南省气象数据探测中心收集了研究区附近的月平均降水量和月平均气温数据。数据采集时间为每日2、8、14、20时4次气象观测，然后计算月平均降水量与月平均气温。

表1 SAR影像参数

2 研究方法

对于获取的112景SAR影像，主要开展差分干涉、时间序列和奇异谱分析3个步骤的处理，详细流程如图2所示。在差分干涉环节，首先进行多视处理，本文试验按10(距离向)×2(方位向)的比例进行多视处理，从而抑制斑点噪声，并使得距离向和方位向具有近似的空间分辨率[15]。其次，由于Sentinel-1 TOPS模式影像在采集时天线会发生旋转，为了避免相邻burst之间发生相位跳变，基于SAR影像相位和强度信息估计距离向和方位向的偏移量，随后拟合偏移量多项式模型，并进行迭代处理，直至配准精度达到千分之一个像元[16]。然后，在选取适当的时空基线生成干涉图后，依次去除干涉图中的地形相位、大气相位、噪声相位及趋势项，将形变相位分离出来并完成相位解缠[17]。由于解缠误差会影响时间序列形变的估计精度，在时序分析前需要剔除含有解缠误差的干涉基线，为此试验中选择利用相位闭合环方法去除含有解缠误差的干涉图，选取最优的干涉对网络。时间序列形变采用短基线(SBAS)干涉测量方法进行求解。最后，为了揭示研究区域膨胀土的形变特征，开展了基于奇异谱分析(SSA)的时间序列形变分解研究，并分析形变与降水、气温之间的关系。

图2 方法流程Fig.2 Flowchart of used method

2.1 相位闭合环选取干涉图

在进行InSAR相位解缠处理时，由于去相干、地形复杂等原因引起的干涉图相位不一致，增加了相位解缠的困难，不可避免地产生相位解缠误差，从而降低了形变参数的解算精度。因此在进行时序分析之前，需要生成最优的干涉网络以提高结果的可靠性。基于干涉网络的冗余性，本文采用了相位闭合环的方法识别并校正解缠误差[18-20]。假设3景SAR影像φ1、φ2和φ3，以及与之相关的3幅解缠相位φ12、φ23和φ13，则相位环闭合差被定义为[16]

φ123=φ12+φ23-φ13

(1)

式中，φ12为SAR影像φ1和φ2之间的解缠相位；φ23为φ2和φ3之间的解缠相位；φ13为φ1和φ3之间的解缠相位；φ123为相位环闭合差。理论上，若3个干涉对φ12、φ23和φ13均不存在解缠误差，则φ123等于0[21]；如果其中干涉图含有解缠误差时，则φ123为2π的整数倍。

本文基于相位闭合环对解缠后的干涉图进行筛选。由于自动校正解缠误差可能会产生更加难以识别的误差，导致时间序列形变精度估计误差，因此保守起见，逐干涉对而非逐像素评估环相位闭合度[22]，计算每幅干涉图相位环闭合差的均方根误差(root mean square,RMS)，将RMS大于给定阈值(1.5 rad)的干涉对从干涉图网络中剔除。图3所示为2017年3月—2021年4月Sentinel-1数据干涉生成的时空基线图，该时间段的数据冗余最多，更适合按照给定阈值自动进行筛选。

图3 时空基线Fig.3 Spatio-temporal baseline network

2.2 时间序列形变估计

小基线集(SBAS)技术[23]是一种时序InSAR分析方法。该方法基于多个主影像，选取较短的时空基线进行干涉以减弱InSAR处理中的去相关问题，最后合并若干个干涉对子集求解形变时间序列及形变速率。目前，SBAS-InSAR方法已经得到多方面的扩展应用[24-26]，可以满足大部分研究区域与数据条件的需要。对于覆盖研究区的N+1景SAR影像，基于短基线原则获取M幅干涉图。假设两景SAR影像的采集时刻分别为ta和tb，并进行差分干涉生成第i幅干涉图，则其任一点的差分干涉相位可以表示为[27]

(2)

式中，dta和dtb分别为ta和tb时刻相对于某一参考时刻的累计位移。本文试验中采用三次模型[28]。如果以t0作为初始时刻，以除t0时刻外任意时刻ti(i=1,2,…,N)的干涉相位φ(ti)作为未知数，以经过相位解缠解算的差分干涉相位φ(tj)(j=1,2,…,M)作为观测值，则有

(3)

φ(ti)=[φ(t1)φ(t2)…φ(tN)]

(4)

φ(tj)=[φ(t1)φ(t2)…φ(tM)]

(5)

将式(2)—式(5)联立可获得M个方程组

Aφ=φ

(6)

式中，A为M×N阶系数矩阵。由于SBAS方法在生成干涉图时按一定的时空基线进行组合，有可能出现不连续的干涉图子集，因此系数矩阵A可能秩亏，导致方程组(6)产生无穷多解。这时需要采用奇异值分解(SVD)方法，对多个短基线干涉图子集联合求解，获得系数矩阵A的广义逆，求出速度矢量在最小范数意义上的最小二乘解，进而求得形变时间序列。

2.3 基于SSA的形变分解

由于环境、系统等的影响，由时序InSAR技术所获取的原始形变时间序列中不可避免地含有噪声误差，限制了形变信息的精度与可靠性。因此，将可靠的形变信号从原始形变时间序列中提取出来并确定其变化趋势，对地表形变的时序分析具有重要意义。SSA方法建立在Karhumen-Loeve分解理论的基础之上，是近年来飞速发展起来的一种研究非线性时间序列的方法。该方法无须先验信息，不受正弦波假定约束，通过对一维时间序列数据进行主成分分析，构造出原始时间序列的轨迹矩阵并在此基础上分解、动力重构，再结合经验正交函数，最终提取出趋势项和周期项等有效信息[29-30]。一维时间序列x=(x1,x2,…,xS)的SSA过程可分为分解和重构两个步骤。首先，选择合适的分解窗口L(1

(7)

然后，采用奇异值分解算法(SVD)对轨迹矩阵XL×N进行分解得到

(8)

(9)

最后，根据式(10)对角平均化重构时间序列(y1,y2,…,yS)

(10)

基于时间序列分解原理，本文首先采用SSA算法得到了InSAR时间序列各分量特征值，随后计算了时间序列各分量特征值的贡献率。基于趋势形变平滑、周期形变规律明显、残差符合白噪声特点等原则，依据特征值累积贡献率，将原始监测数据划分为趋势项、周期项和噪声项。

3 结果和分析

3.1 形变空间分布特征

利用收集的112景Sentinel-1A数据获取了观测时段的平均形变速率，如图4(a)所示，其中绿色代表稳定区域，黄色及红色代表下沉区域，蓝色代表上升区域。由图4中可以看出，InSAR观测点位主要位于附近村庄建筑物密集区域，以及南水北调两侧人工修筑的边坡(图4(b))，大部分农田则表现为失相干。空间分布上，除了东北角之外，南水北调干渠沿线大部分区域表现为抬升形变，最大年抬升速率约为18 mm/a。而干渠沿线东北角表现为下沉形变，最大年下沉速率约为15 mm/a。除此之外，干渠东侧区域也表现为较为明显的下沉，最大年下沉速率达到22 mm/a。进一步分析图4(a)和图1(b)，发现研究区的形变空间分布与地形变化相关，即挖方区域表现为抬升形变，填充区域表现为下沉形变。同时发现，形变量大小与挖填方深度相关，挖方越深抬升形变越大，填方越深下沉形变越大。

图4 研究区地表形变速率和挖方渠道现场Fig.4 Surface deformation rate of study area and photo of excavation channel

3.2 形变时间分布特征

为了进一步调查研究区膨胀土的时间序列分布特征，绘制了图4(a)中A、B、C和D4个点2017年3月—2021年4月近4年的时间序列形变，如图5所示。A点位于干渠堤岸边坡上，挖方深度为35 m，可以看到整体表现为明显的抬升形变，观测时段内累计抬升达到48 mm，时间上经历先缓慢下沉后加速抬升再缓慢抬升的变化特征。B点同样为挖方区域，但是挖方深度只有25 m，相应的累计抬升量小于A点，观测时段内约为20 mm，时间变化特征与A点相似，经历缓慢下沉后加速抬升再缓慢抬升。C点位于膨胀土填方渠道上，填方厚度为15 m，总体表现为持续的下沉形变，累计下沉量为63 mm，时间上具有一定的波动性，每年夏季表现为加速下沉特征。D点位于干渠东侧填方区域，填方厚度为25 m，累计下沉量达到90 mm，时间波动性与C点相似。通过A、B、C和D点的时间序列分析，进一步证实了挖方抬升、填充下沉的形变特征，且形变量大小与填挖深度正相关，时间序列上表现为一定的波动性。同时发现在相同的挖填方深度情况下，填方下沉速度大于挖方抬升速度，如图5中的B点。

图5 研究区地表时间序列形变Fig.5 Time series deformation of the study area

3.3 时间序列形变分解和重构

膨胀土具有吸水剧烈膨胀、失水剧烈收缩及反复胀缩形变的特性，对其时间序列形变进行分解和重构，提取趋势项、周期性和噪声项，对分析膨胀土形变规律具有重要意义。为此，本文利用SSA方法，基于趋势形变平滑、周期形变规律明显、残差符合白噪声特点的原则，对图5中的A点深挖方抬升形变和D点高填方下沉形变进行了分解和重构，如图6所示。考虑到本文时间序列观测个数为97及填挖方形变特征，设置A、D点奇异谱分析窗口长度l为12。图6(a)为A点趋势形变，其主要反映了上部膨胀土卸载情况下，下部土壤自然膨胀的过程，表现为抬升形变，占据图5中A点总体形变量的75%。图6(c)为A点周期形变，为降雨、气温、蒸发、地下水、荷载等外部因素引起的膨胀土形变，占据图5中A点总体形变量的12%。可以发现，除了2017年之外，A点在2018—2021年表现为每年4至9月抬升形变，10月至来年3月下沉形变的周期性形变规律。这一变化与大气降雨时间强相关，可以推测挖方膨胀土在每年4至9月集中降雨阶段遇水膨胀，每年10月至来年3月少雨阶段失水收缩。图6(e)为A点噪声形变，其主要由SAR成像几何、数据处理等引起的误差项，具有随机性的特点，占据总体形变量的13%，具体分析中应去除此部分贡献。

图6 基于SSA的时间序列分解和重构Fig.6 Decomposition and reconstruction of time series deformation based on SSA method

与A点相比较，D点呈现出相反的趋势形变(图6(b))，主要为未压实的填充膨胀土自然固结而产生的下沉形变，占据图5中D点总体形变量的81%。图6(d)中D点的周期项形变也主要受降雨、温度等外在因素的影响，发现2017年和2019年具有较大波动，2018年和2020年波动相对较小，占据图5中D点总体形变量的6%。图6(f)为A点噪声形变，占据总体形变量的13%。

4 形变成因分析

膨胀土边坡渠坡形变是内外因综合作用的结果。结合图1(c)以及相关资料分析，研究区膨胀土主要成分为伊利石，其次是蒙脱石，此外还含有少量的高岭石、绿泥石和水云母等[31-32]。该类膨胀土的渗透系数和自由膨胀率较大，具有较强的涨缩潜势。因此，土体物理力学性质等内因对形变起主导作用。对于挖方膨胀土边坡，由于膨胀土的超固结性和胀缩性，渠道开挖过程中会产生卸荷回弹，在垂直方向产生抬升形变，如图5中的A、B点。对于填方膨胀土边坡，由于膨胀土的超固结性，填方体的自然沉降固结导致沉降产生，如图5中的C、D点。由于干渠两侧的膨胀土本身含水量较大，土体的膨胀性不是很强，因此抬升形变量较小。

膨胀土边坡形变的外部原因主要是降雨、蒸发、温度、地下水等引起的膨胀土土体含水率变化，进而导致膨胀土渠道胀缩。已有试验结果表明，每次干湿循环(降雨与蒸发交替)，膨胀土边坡均累积了向坡下的沉降和水平位移。其中，降雨决定土体湿化程度及入渗深度，是渠坡形变的最直接气候因素。降雨入渗有一定滞后性，且由于地表径流、水分蒸发等原因，入渗量要小于降雨量，因此由于降雨引起的形变具有一定的滞后性。图7(a)为深挖方A点周期形变与月平均降雨量之间关系，可以看到随着降雨量的逐渐增大，A点抬升形变逐渐增大，且形变会滞后2～3个月，之后随着降雨量的减小形变开始下沉。气温变化会引起土体温度变化，进一步引起含水率变化，从而导致膨胀土渠坡形变。图7(b)为深挖方A点周期形变与月平均气温之间关系，可以看到温度逐渐升高形变逐渐增大，温度逐渐降低形变也逐渐降低，且也表现出一定的时滞性，表明温度与形变具有较强的相关性。地下水是影响渠坡形变的另一重要因素，地下水不仅可加快结构面软化，使得滑面抗滑力降低，还能在底滑面上提供扬压力、在后缘拉裂面提供静水压力，导致渠坡形变和滑坡启动。由于本文中的地下水位略高于渠道水位，因此渠道形变外部因素主要是降雨和温度。

图7 膨胀土边坡形变与降雨和温度之间的关系Fig.7 The relationship between deformation and precipitation and temperature

5 结 论

本文以南水北调中线工程某一膨胀土边坡为例，采用多时域InSAR和SSA技术提取了膨胀土边坡趋势项和周期项形变，并结合膨胀土物理特性以及降雨、温度数据，分析了膨胀土形变特征。得到如下初步结论：

(1) 挖方膨胀土边坡表现为抬升形变，而填方边坡表现为下沉形变。利用112景Sentinel-1A数据获取的年平均形变速率显示南水北调干渠沿线大部分挖方区域表现为抬升形变，最大年抬升速率约为18 mm/a，而干渠部分填充区域表现为下沉形变，最大年下沉速率约为15 mm/a。

(2) 膨胀土边坡时间序列形变呈现一定的周期性变化，具有明显的胀缩特性。通过SSA对InSAR获取的时间序列形变进行分解，发现挖方和填方膨胀土边坡周期项形变具有明显的胀缩性，具体表现为挖方膨胀土每年4至9月出现抬升形变，10月至来年3月出现下沉形变。同时发现，膨胀土形变量大小与填挖深度正相关，挖方越深抬升形变越大，填方越深下沉形变越大。在相同的挖填方深度情况下，填方下沉速度大于挖方抬升速度。

(3) 卸载回弹和土体固结是膨胀土边坡形变的内在因素，降雨和温度则是形变的外部因素。考虑到膨胀土的超固结性和胀缩性，渠道开挖过程中的卸荷回弹是挖方膨胀土形变的内因，而填方体的自然固结是填方膨胀土下沉形变的内因。通过分析膨胀土形变与降雨、温度之间的关系，发现膨胀土周期性形变与降雨、温度强相关，其造成的形变会时滞2～3个月。