近年来，由于人类活动与自然环境变化的影响，中国西部地区滑坡等地质灾害频发，并有逐年加重的趋势，给人民生命财产安全带来了巨大威胁^[1]。由于地质环境脆弱且地形复杂，黄土高原地区饱受水土流失与滑坡灾害的困扰。通常情况下，非常陡峭的黄土坡体也可以长时间保持稳定。由于黄土孔隙度大，一旦遇水，即使平缓的黄土坡体也有可能出现变形甚至失稳^[2]。1983-03-07，甘肃东乡洒勒山突发巨型灾难性高位远程滑坡，滑坡向前推进了1 000余米，最终造成237人死亡^[3]。甘肃永靖黑方台地区作为刘家峡库区、盐锅峡库区和八盘峡库区移民安置地，受到常年累月的灌溉影响，年均发生3~5次滑坡灾害，给当地居民带来严重的困扰^[4]。陕西泾阳地区由于人类工程活动密集频繁，自1976年以来发生滑坡40余次^[5-6]。因此，急需对黄土滑坡分布进行识别与长期监测，及时发现险情，为防灾减灾提供必要的依据。

2017年茂县新磨村滑坡和2018年金沙江白格滑坡等破坏力巨大的突发滑坡事件发生后，国家对滑坡泥石流等重大自然灾害的识别、监测、预警提出了迫切需求。地质灾害防治工作的重心也从被动避灾和治灾逐渐向主动防灾转移，事先发现和识别滑坡等灾害隐患，是主动有效防灾的重要前提。由于中国西部地区地形起伏大，传统地面测量手段受成本和覆盖范围的限制，难以实现大规模灾害普查。合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术具有不受云雾影响、覆盖范围大和监测微小地表形变的能力，可以为地质灾害普查提供重要的支持。

中国在2000年将InSAR技术引入滑坡监测中，并在三峡库区链子崖滑坡和新滩滑坡等典型滑坡体上布设人工角反射器，开展形变监测可行性分析与精度验证^[7-8]。随着合成孔径雷达(synthetic aperture radar, SAR)影像的不断积累与处理算法的成熟，永久散射体干涉测量技术^[9-10]、小基线数据集算法^[11-13]与SqueeSAR^[14]等算法利用时间序列SAR数据集中的稳定散射点目标测量地表变形，精度可达毫米级。上述技术已在三峡库区^[15]、金沙江乌东德水电站库区^[16]、白龙江上游^[17]、大渡河上游丹巴县^[18-19]等地区的滑坡识别和监测中得到了成功应用。

本文获取了2016-01至2018-08时间段内覆盖甘肃永靖黑方台及其周边区域的升降轨Sentinel-1数据集，开展了时序InSAR分析。通过与GPS观测资料对比，验证了InSAR结果的有效性，并进一步深入分析了典型区域的变形特点与影响因素，此外还定量分析探讨了温度波动与桥梁形变之间的关联关系。

1 时序InSAR分析技术

InSAR技术可以测量地表微小变形，但由于时间去相干、空间去相干的影响，且干涉相位中包含轨道误差、数字高程模型(digital elevation model, DEM)误差和大气传播延迟等干扰成分，其应用和测量精度往往受到限制。要想高精度地获取地表形变观测，应尽量避免去相干因素的影响，有效分离出干扰项。时序InSAR分析技术就是在这样的背景下发展起来的。SAR影像中的像素可以划分为永久散射体(persistent scatterer，PS)主导和分布式散射体(distributed scatterer，DS)主导两大类。典型的PS像素主要包含人工地物和裸露岩石等，相位稳定性受时间和空间去相干的影响较小，是时序InSAR分析的首选对象。DS像素中包含低矮植被和裸地等，相位稳定性相对较低，只能在短时间内保持相干性。因此，在对DS目标进行分析时，需充分利用短时间基线的干涉图，并通过滤波来进一步提高干涉相位的稳定性。

在对获取的SAR数据集进行主影像选取、配准、重采样之后，采用短时间基线和空间基线的组合策略生成序列差分干涉图。对得到的差分干涉图进行滤波，进一步增强干涉图中DS目标的相位稳定性。

在时序分析中，通常希望保持PS像素的相位不变，对DS像素滤波来提高相位稳定性^{[14, 20]}。传统的窗口滤波不区分像素的类型，直接进行等权或加权平均来获取中心像素的相位，不能实现上述目的。局部范围内相似的地物目标在SAR影像中对应的回波也是相似的，在数学上也服从相似的统计分布，这样具有统计相似性的点就可以称为同质点。因此，在自然场景中，PS像素分布稀疏，相应的同质点少，而DS像素则相反。因此，可以通过检验同质点数量的多少来实现PS和DS像素的初步区分。在SAR影像数量足够的情况下，可通过假设检验来选取每一个像素的同质点集。本文采用快速同质点选取方法进行同质点集的选取，该方法将高斯分布条件下的假设检验问题转化为置信区间估计，有效提高了计算效率^[21-22]。对同质点数量大于20的像素，对其干涉相位应用加权滤波，低于此阈值的像素保持干涉相位不变。滤波后，干涉图中每个像素的相位值φ可以表示为：

$ \varphi=W\left\{\varphi_{\mathrm{disp}}+\varphi_{\mathrm{orb}}+\varphi_{\mathrm{atm}}+\varphi_{\mathrm{dem}}+\varphi_{n}\right\} $

(1)

式中，等式右边各分量分别为地表形变φ_disp、轨道误差φ_orb、大气φ_atm、高程误差φ_dem和噪声φ_n；W{·}表示缠绕运算。前3项相位分量在空间上相关。对干涉图滤波之后，可以通过振幅离差阈值的方法实现PS和DS的初选^[23]，并通过相位稳定性测度确定最终的PS和DS集合。相位稳定性可以通过初始的缠绕相位来估计。首先利用带通滤波估计缠绕相位中的空间相关分量，利用高程误差与基线之间的线性关系估计高程误差相位。将上述相位从原始的缠绕相位中移除后，按照式(2)，基于相位残差φ_rⁱ估计每个像素的相位稳定性γ：

$ \gamma=\frac{1}{N}\left|\sum\limits_{i=1}^{N} \mathrm{e}^{j \varphi_{r}^{i}}\right| $

(2)

式中，N表示干涉图数量。根据预设的相位稳定性阈值确定最终点集，并进行离散时空三维相位解缠。获取解缠相位后，就可以根据各相位分量的不同特性逐一加以分离。

解缠相位中的轨道趋势可以通过一个双线性多项式进行拟合。数字高程模型(dlgital elevation model, DEM)误差则是采用最小二乘法，通过构建基线和高程误差之间的线性关系模型来估算。地形起伏较大区域的大气扰动分量可以分为垂直分层和湍流混合。其中，垂直分层信号可以通过高程与大气之间的线性关系获取，湍流混合可以通过组合滤波的方法求得^[24]。将上述3个分量从解缠相位中去除后，对剩余相位进行奇异值分解，即可获得时间序列形变反演结果。

时序InSAR技术形变反演结果的精度主要由SAR影像数量、波长和时间基线等因素决定。其理论精度可用式(3)^[25]表示为：

$ \delta_{\Delta v}^{2} \approx\left(\frac{\lambda}{4 \pi}\right)^{2} \frac{\delta_{\varphi}^{2}}{M \delta_{B_{t}}^{2}} $

(3)

式中，δ_Δv²表示时序InSAR技术获取的形变速率精度；δ_φ²为相位离差；λ为波长；M为SAR影像数量；δ_Bt²为时间基线离差。

2 实验区与实验数据

本文实验区位于黄土高原西南缘甘肃兰州西固区与永靖县交界处，距离兰州市区约60 km，这里也是湟水河汇入黄河的地方。从上世纪50年代起，此地区陆续修建了刘家峡水电站、盐锅峡水电站和八盘峡水电站^[26]。伴随着水电工程的建设运行，实验区中的黑台、方台以及赵家台等适宜耕种的平坦台塬被开发成为移民安置点。黄土高原降水少，蒸发量大，因此当地村民大量抽取黄河水对台塬地区进行灌溉，以保证作物生长。一般一年灌溉4~6次，7月下旬的灌水量最大^[27]。与此同时，台塬地区地下水位快速上涨，土壤中水分的增加破坏了原有的应力结构，导致了大量滑坡的发生，引起了各方的重视。滑坡的存在给当地交通运输安全和居民人身财产安全造成了巨大威胁。其中，黑台和方台组成的黑方台滑坡群最受瞩目，截至1989年，黑方台的累计沉降已经达到1.6~1.8 m^[26]。自1968年以来，有据可查的滑坡记录约200次^[4]。图 1给出了本文实验区的位置和范围。

本文实验数据为2016-01至2018-08获取的覆盖实验区的两组Sentinel-1序列数据集。其中，升轨数据集总共有63景，中心入射角为46.1°；降轨数据集有54景，入射角为36.9°。两组数据集的干涉图组合如图 2所示。采用覆盖实验区30 m分辨率的SRTM(shuttle radar topography mission)DEM作为参考高程数据，用于辅助去除地形相位。同时，还获取了当地一个GPS站点2017-11-27-2018-04-20的变形观测资料，用于和InSAR处理结果进行比对验证。

3 滑坡监测实验结果与分析 3.1 平均形变速率分布图

图 3给出了利用Sentinel-1升降轨数据集获取的实验区雷达视线向(line of sight, LOS)平均形变速率分布图。由于植被覆盖较少，因此两个数据集均获取到了密集的观测点。其中，升轨数据集获取的点密度为3 642个/km²，降轨数据集的点密度为2 841个/km²。两个数据集探测到的变形区域分布高度一致，主要集中分布在黑方台、赵家台和三马台等灌溉台塬的边缘区域。其中，黑方台区域的形变最为显著，最大形变速率约60 mm/a；盐锅峡大坝与八盘峡大坝未探测到形变信号。除此之外，由于黄河穿越本实验区，实验区内修建了许多大大小小的桥梁，例如河口特大桥、盐锅峡黄河大桥和兰新铁路黄河大桥等。本文实验在河口特大桥上探测到温度相关形变信号，具体情况将通过时序变形分析加以探讨。

3.2 党川滑坡群变形特点

党川滑坡群位于黑台西北侧。图 4(a)和图 4(b)给出了2002-04-20以及2018-10-13获取的该区域光学影像。从两个时相的影像对比可以看出，党川滑坡地表发生了巨大变化，黑台的边缘区域更加支离破碎，出现了4个明显的陡坡，如图 4(b)的白色虚线所示。图 4(c)和图 4(d)给出了Sentinel-1升降轨数据集获取的视线向平均形变速率图。其中，正值表示朝向卫星移动，负值表示远离卫星运动。由于观测几何的差异，升轨数据集可以直接照射到党川滑坡群，而降轨数据集上党川滑坡群的部分区域位于阴影区域，因此可以看到降轨数据集在滑坡顶部探测到的点非常稀疏。本文获取了一个安装于黑台边缘GPS站点GP203的时序形变监测值，位置如图 4(c)和图 4(d)中的白色正方形所示。实验中将GPS时序观测转换到雷达视线向^[28]，图 5给出了GP203与邻近PS观测值的对比结果。由于GPS和PS点之间不完全重合，并且两种观测都有一定的误差，因此两种测量结果并不完全一致。将观测日期相同的GPS值与升降轨InSAR观测值作差，差值的均方根分别为1.4 mm和1.0 mm，升降轨InSAR数据与GPS观测之间差值的标准差分别为1.7 mm和4.9 mm。总体来讲，GP203观测和时序InSAR技术的观测结果具有较好的一致性。除了GP203之外，本文还选取了图 3(a)中所示的P₁和P₂两个点进行时序分析，并尝试进一步分析滑坡变形的影响因素。

3.3 黄土滑坡变形的影响因素分析

黄土台塬坡体稳定的影响因素主要是持续灌溉^{[4, 29]}，降雨、泉水与其他人类活动带来的影响范围非常有限。图 6给出了P₁和P₂点的时序变形。由于每年7月下旬的秋灌水量最大^[27]，从图 6(a)和图 6(b)中可以看到，7月之后坡体的变形趋势都会发生变化。虽然整个时序上坡体一直都在变形，但秋灌之后有一定程度的加速，在图 5(a)和图 6(a)中表现最为明显。升降轨观测到的P₁点形变趋势和量级之间的差异主要是由两个轨道的观测几何不同带来的地表三维变形在不同数据集上的投影差异所造成的。由于秋季集中灌溉，地下水位急速上涨，黄土层中土壤水分含量相应增加，降低了坡体强度，导致坡体变形甚至失稳滑动。

同时，可以发现黑台北侧磨石沟滑坡群周围的PS点形变呈现线性趋势。如图 6(c)和图 6(d)中的P₂点形变序列所示，该点在2016-2018年间的变形趋势整体稳定。由于黑台地势北高南低，现场调查未见地下水露出^[26]。因此，此区域变形的影响因素主要为降雨和灌溉。但在本文研究时间段内，未发现这两种因素对滑坡变形趋势的显著影响。

3.4 河口特大桥变形特征分析

虽然本文的主要研究对象是黄土滑坡，InSAR时序分析也可以识别到研究区域内黄河大桥的温度相关变形特征。大型人工建筑例如高楼、桥梁等一般都会受到热胀冷缩的影响，因此，环境温度的变化也会引起建筑表面相应的变形，变形的程度主要跟建筑材料的物理特性有关^[30]。本文以河口黄河特大桥为例进行说明。河口黄河特大桥主跨360 m，为一座斜拉索大桥。图 7(a)和图 7(b)中给出了从Sentinel-1升降轨数据集获取的桥上B点的时序变形与温度差之间的关系，可以看到，温差变大，桥梁的变形也随之线性增加。图 7(c)和图 7(d)将形变和温度差进行线性拟合，Sentinel-1升降轨观测的形变与温度差之间的关系可以表示为：

$ \left\{\begin{array}{l}{\Delta d_{\mathrm{asc}}=-2.63 \Delta T+11.69} \\ {\Delta d_{\mathrm{desc}}=-0.46 \Delta T+2.16}\end{array}\right. $

(4)

式中，Δd_asc、Δd_desc分别表示升、降轨形变量；ΔT表示温度差。在升降轨反演结果中，二者之间的线性相关度绝对值(R²)分别达到了0.92和0.73，一定程度上也验证了InSAR观测的可靠性。

4 结语

本文利用时间序列InSAR方法对升降轨的Sentinel-1星载SAR序列数据进行分析，探测识别出了2016—2018年期间甘肃永靖黑方台及其周边灌溉台塬的不稳定区域，分析了典型滑坡体的形变演化特征。通过InSAR观测结果与GPS观测资料的对比，验证了InSAR处理方法的有效性。通过时序分析发现，持续的灌溉行为是本区域坡体失稳的主要诱发因素。因此，灌溉时期需密切注意台塬边缘坡体的稳定性。与此同时，本文发现研究区域内的跨黄河大桥受到热胀冷缩的影响，大桥变形和季节性温度变化呈现高度的线性相关。随着Sentinel-1卫星星座的业务化运行，雷达数据不断积累，SAR数据可提供丰富的研究资料。由于卫星观测几何的限制，未来有必要将多轨道数据观测进行融合，获取地表三维时序形变场，更好地表征地面目标的形变模式和特点。