甲居滑坡位于四川省丹巴县聂呷乡甲居村，滑坡长1 200 m，宽1 000 m，滑体平均厚度22 m，体积约为2 600万m³，属特大型牵引式中深层土质滑坡^[1]。为监测该滑坡的活动性，文献[2-4]采用GPS和InSAR技术进行了周期性观测和融合处理。但由于GPS时间和空间采样率均较低，难以监测该滑坡的时间发育特征，而且研究区域地形起伏较大，特别是植被茂密导致传统差分InSAR技术在该地区的时空失相干非常严重，形变监测结果精度有待进一步提高，限制了差分InSAR技术的应用。

人工角反射器干涉技术 (corner reflector-InSAR) 通过在监测体布设一定数量的人工散射体来克服时空失相干对差分InSAR技术的影响，从而用来高精度探测地形复杂区域的形变。目前，常用的CR-InSAR解算方法有Lambda算法^[5-8]和最小费用流算法等，这些算法要求CR点能够成功解缠，但是，很多时候两个CR点间的形变相位超过雷达半波长，相位解缠变得异常困难。为此，Xia等提出CR-InSAR周期图算法^[9-10]，针对每个CR点建立一个周期图，这样对于每个CR点，其相位时间序列可以看成是一个关于时间间隔和垂直基线的二维谐振信号，因此，形变估计问题就转化成频谱估计问题^[11-19]。但是该算法进行参数估计时，对干涉对数量、各干涉对时间和空间基线分布的均匀性都有较高要求，而且在获取非线性形变时经常采用滤波的方法，导致形变结果较为平滑，非线性形变不明显。

本文借鉴InSAR小基线集算法和CR-InSAR周期图算法, 采用一种非线性CR-InSAR算法，对甲居滑坡进行了非线性CR-InSAR试验，通过与GPS监测形变结果一致性的比较，表明该算法比原CR-InSAR算法获取的滑坡形变结果更为合理。

1 CR-InSAR解算原理及算法

假设有M+1个角反射器，获取N+1幅SAR影像，按照传统小基线集InSAR处理方法，选取其中一幅影像作为主图像，可以组成N个SAR干涉对，每个干涉对包含了M+1个CR点对应的相对于主影像的差分相位。通过SAR轨道信息和外部DEM去除地形相位后，第i个CR点在第m个差分干涉图中的差分相位表示为：

式中，R_im是对应的形变；B_m是垂直基线；h_i是第i个CR点的高程改正数；θ_m是入射角；φ_{im, atom}是大气延迟相位，φ_{im, noise}是噪声相位。

如果将CR点形变看成近似线性形变，则式 (1) 可以写成：

式中，V_i是第i个CR点的形变速率；T_m是第m个干涉图的时间间隔。当两个CR点间的距离为1 km时，大气延迟相位可以认为是空间强相关的，如果选取其中一个比较稳定的CR点作为参考点，则其他的CR点相对于参考CR点的差分相位可以表示为:

1.1   基于解缠相位的CR-InSAR算法

对式 (3) 采用最小费用流或Lambda算法解缠得到：

式中，φ_{i, ref}为解缠相位，对每个CR点的解缠相位对N个干涉图采用最小二乘准则 (或奇异值分解) 解算，获取每个CR点的累计形变时间序列结果，该方法的应用前提是要求CR点相位能够准确解缠。但当两CR点间的形变相位差大于|π|时，无法准确解缠，从而导致形变时间序列错误。

1.2   基于二维周期图的线性形变算法

文献[9]等提出基于二维周期图的线性形变解法，如果将CR点形变看成线性形变，选择位于稳定区域的CR点为参考点，假设其形变速率为0，建立如下模型：

式中，Δφ_i为第i个CR点相对于参考点的差分干涉相位；V_i为第i个CR点的线性形变速率；h_i为第i个CR点相对于参考点的高程改正数，若两个CR点间的距离不是很远，大气延迟相位忽略不计，式 (5) 可以写成：

式中，f_iv=2V_i/λ；t=T_m；f_ih=2h_i/λR_msinθ_m；b=B_m；φ_{i, noise}是独立的高斯噪音，均值为0。

为了求解形变速率相关的频率f_iv以及和高程改正值相关的频率f_ih，对于每个CR点构造一个周期图：

当E_{(fv, fh)}达到最大值时，根据对应的频率f_iv和f_ih很容易得到所要求的形变速率值以及高程改正值。然而干涉对数量、干涉对时间基线和空间基线分布的均匀性都会对参数估计的精度产生影响，而且采用滤波的手段获取非线性形变会导致形变结果线性较为明显。

1.3   基于相位相干性的非线性CR-InSAR算法

本文采用一种改进的基于相位相干性的非线性CR-InSAR算法，定义相位相干性为：

该约束条件相当于一个最优函数，依据一定的先验信息确定解空间的大小、位置和搜索步长，即可较快地获得参数的最佳估计值。当估计参数越接近真值时，γ就越大，也就是说γ越大，参数求解的精度就越高。

因此, 给定一个二维解空间，在解空间内搜索γ的最大值，与最大值对应的形变速率值和高程改正数值就是所求解的参数，CR点相位相干性的三维图如图 1所示。

图  1  CR点相位相干性的三维图

Figure  1.  3D Phase Coherence of CR Points

下载: 全尺寸图片 幻灯片

将与形变速率和高程改正数有关的相位从差分干涉相位中减掉，可得到残余相位：

依据最小二乘准则将非线性形变相位当成未知参数进行求解，由于大气相位在时间域上表现为高频信号，而在空间域上表现为低频信号。因此，当CR点距离小于1 km时，将其和噪声相位一起当成随机误差处理，可以得到：

式中，为非线性形变相位，S_nl为非线性形变。

因此，CR点的形变可以表示为：

基于相位相干性的非线性CR-InSAR数据处理流程如图 2所示。

图  2  基于相位相干性的非线性CR-InSAR数据处理流程图

Figure  2.  Flowchart of Nonliner CR-InSAR Data Processing Based on Phase Coherence

下载: 全尺寸图片 幻灯片

2 模拟实验

为了验证本文采用算法的可靠性，首先对模拟数据进行了实验分析，模拟线性形变速率为-0.4 mm/d，高程改正数为50 m，非线性形变采用正弦函数进行模拟，分别采用Lambda算法、二维周期图算法以及本文采用的算法进行求解，并与模拟真值进行对比，对比结果如图 3所示。

图  3  模拟实验结果对比图

Figure  3.  Comparison of Simulation Results

下载: 全尺寸图片 幻灯片

从图 3可以看出，由于模拟数据的形变尺度较大，导致Lambda算法存在明显解缠误差，而本文采用基于相位相干性的算法求解得到线性形变速率为-0.387 mm/d，高程改正数为65.9 m，而且与二维周期图算法求解结果相比，非线性形变也恢复较好。从对比结果来看，本文算法不需要相位解缠即可获取CR点的形变信息，求解结果与模拟真值较为一致，对大尺度形变梯度具有较强探测能力，现将该算法应用于甲居滑坡CR点形变监测分析。

3 数据处理方案与结果分析

图 4为四川甲居滑坡地形图，自2007年布设了10台角反射器，30个GPS监测点，每年进行4期GPS观测。结果显示，该滑坡形变具有分区特征，可分为强变形区、弱变形区和潜在形变区。为此获取了覆盖研究区域自2007-08至2008-06的9景Envisat ASAR数据，开展InSAR监测与试验。

在InSAR数据处理中采用90 m分辨率的SRTM DEM，该DEM在山区的平面精度为20 m，高程精度为20 m^[11]。

图  4  甲居滑坡地质和形变监测图

Figure  4.  Geological and Deformation Monitoring Chart of Jiaju Landslide

下载: 全尺寸图片 幻灯片

首先通过小基线集InSAR来检验研究区域的相干性以及结果中可能存在的误差。在此基础上对比两种CR-InSAR解算方法，并对解算结果与同期GPS监测结果进行比较，验证本文方法的可靠性。

3.1   小基线集InSAR结果

本文以Gamma软件作为差分InSAR干涉处理平台，共生成36个干涉对，基线分布如图 5所示。图 6(a)为生成的36个干涉对的相干图；图 6(b)为对应的解缠图。

图  5  甲居地区获取的SAR数据基线分布图 (20080104)

Figure  5.  Baseline of Envisat ASAR Data (20080104)

下载: 全尺寸图片 幻灯片

图  6  36个干涉图相干图和解缠图

Figure  6.  Coherence Chart and Unwrapping Chart of 36 Interferograms

下载: 全尺寸图片 幻灯片

从相干图和解缠图中可以看出，研究区域由于植被等因素的影响，相干性很差，导致解缠误差难以消除，因此，本文仅对相干性最好的两个干涉对 (2008-01-04~2008-02-08和2008-02-08~2008-03-14) 进行常规差分处理，获得的形变结果如图 7所示，两个干涉对均存在明显形变信息，最大沉降量3 cm左右 (35 d)。

图  7  两个差分形变图

Figure  7.  Two Differential Deformation Chart

下载: 全尺寸图片 幻灯片

3.2   基于解缠相位的CR-InSAR解算

由于角反射器具有很强的后向散射特性，它在SAR图像中表现为一个光斑状的点目标，很容易在强度图中将其与周围的相干地物区分开来，图 8为安装在甲居滑坡上的角反射器在强度图中的近似位置。

图  8  甲居滑坡角反射器在强度图中的显示

Figure  8.  Display of CRs in Intensity Map

下载: 全尺寸图片 幻灯片

在获得CR点的近似位置周围开一个窗口，通过多倍内插计算其幅度的峰值的精确位置，从而获取CR点亚像素级的精确坐标。图 9为某一CR点子像元级坐标幅度图。

图  9  某一CR点的子像元级幅度图

Figure  9.  Sub-pixel Level Amplitude Diagram of CR

下载: 全尺寸图片 幻灯片

直接提取35个干涉对解缠图中的解缠相位，采用最小二乘解算获取CR点的形变时间序列，结果如图 10所示。

图  10  基于SBAS的CR形变时间序列图

Figure  10.  Time Series Deformation of CR Based on SBAS

下载: 全尺寸图片 幻灯片

从图 10可以看出，形变结果在±1 cm之间变动，造成此现象的原因是由于SAR干涉对存在解缠误差 (图 6(b))，导致CR解算结果存在解缠误差，形变时间序列结果不正确。

3.3   基于二维周期图的CR-InSAR算法

根据式 (5)~式 (7)，以位于稳定区的S₉号CR点为参考点，其他9个CR点的形变时间序列如图 11所示。

图  11  基于二维周期图法的CR形变时间序列图

Figure  11.  Time Series Deformation of CR Based on 2D Periodogram

下载: 全尺寸图片 幻灯片

从图 11中可以看出，基于二维周期图法的CR-InSAR算法能够有效避免相位解缠误差，但该方法由于使用滤波的方法分离非线性形变，导致非线性形变有一定损失。

3.4   基于相位相干性的非线性CR-InSAR解算方法

与§3.2不同，本文算法仅提取CR点对应的差分缠绕相位，基于式 (8)~式 (11)，解算出CR点的形变时间序列。

以位于稳定区的S9号CR点为参考点，其他9个CR点的形变时间序列如图 12所示，累计形变结果如表 1所示。

图  12  基于相位相干性的非线性CR形变时间序列图

Figure  12.  Nonlinear Time Series Deformation of CR Based on Phase Coherence

下载: 全尺寸图片 幻灯片

表  1  甲居滑坡角反射器累计形变结果表

Table  1.  Accumulative Deformation of CR Points

形变区域	稳定区		北侧滑坡	南侧滑坡	弱形变区	强形变区
CR点号	S5	S4	S1	S8	S7	S6	S2	S3	S0
累计形变量/cm	-1.8	0.4	-12.7	-3.9	-4.8	-10.3	-9.7	-12.4	-3.9

下载: 导出CSV 

| 显示表格

从上述3种算法解算结果来看，非线性形变都不是特别明显，分析原因主要是由于CR点间距离较近，非线性形变在空间域上表现出一定的相关性，两个CR点间干涉相位做差分的同时，非线性形变也会有一定的损失，再采用滤波方法分离非线性形变，会导致非线性形变进一步损失，因此我们根据大气相位在空间域和时间域表现出的不同特性，将其和噪声相位一起当成随机误差处理，采用最小二乘的方法，将非线性形变当成参数求解，尽可能减小非线性形变的损失。

从上面CR点的累计形变结果可以看出，位于参考点附近的S₄号CR点比较稳定，累计形变量0.4 cm，其他CR点均出现明显沉降形变，其中位于北侧滑坡体上沿南面的S₅号CR点累计形变量-1.8 cm，位于北侧滑坡体上的S₁号CR点累计形变量-12.7 cm，位于南侧滑坡体上的S₈号CR点累计形变量-3.9 cm，位于弱形变区的S₇号CR点累计形变量为-4.8 cm，而位于强形变区的S₆、S₂、S₃和S₀号CR点累计形变量分别达到-10.3 cm、-9.7 cm、-12.4 cm和-3.9 cm。

在甲居滑坡上布设的10个角反射器中有6个的周围布设了GPS监测站，为了更好地验证本算法的准确性，将基于本文CR-InSAR解算结果和基于解缠相位解算结果与GPS解算结果进行比较，如表 2所示。

表  2  基于相位相干性CR解算结果、基于解缠相位CR结果与GPS解算结果对比表

Table  2.  Comparison Between CR Results Based on Phase Coherence and Phase Unwrapping with GPS Method

GPS点号	G10	G8	G12	G13	G16	G15
GPS解算结果/cm	-9.5	-4.6	-4.4	-3.3	-11.0	-3.6
CR点号	S1	S7	S2	S6	S3	S0
基于相位相干性结果	-12.7	-4.8	-9.7	-10.3	-12.4	-3.9
与GPS差值	-3.2	-0.2	-5.3	-7.0	-1.4	-0.3
基于解缠相位结果	0.9	-0.2	0.2	-0.3	-0.4	0.4
与GPS差值	10.4	4.4	4.6	3.0	10.6	4.0

下载: 导出CSV 

| 显示表格

从表 2第5行和第7行可以看出，基于解缠相位的解算结果与GPS解算结果相差较大，而基于相位相干性的CR-InSAR解算结果与GPS解算结果在反映甲居滑坡形变特征上有较好的一致性，但S₁、S₂和S₆号CR点的解算结果与GPS解算结果相差较大，对于与GPS结果存在差异，原因有以下几个方面：①在监测点的位置上，GPS点和CR点之间存在一定的距离，我们将CR点附近的GPS点作为参照，提供一定的先验信息；②在结果的获取时间上，由于没有对应的GPS监测结果，只获取2006-08~2009-11共39个月的累计形变结果，然后取平均求得对应InSAR数据覆盖的10个月的结果，也可能导致存在差异；③在变形特征上，甲居滑坡存在一定的水平形变，InSAR对水平方向形变不敏感，而且GPS在高程方向上的监测精度较低 (1 cm)，将GPS上三维形变结果投影到视线方向时也存在一定的误差。但从CR点整体解算结果来看，较为真实地反映了该滑坡的形变特征。