基于重复轨道的差分合成孔径雷达干涉测量(differential interferometric SAR,DInSAR)技术始于1989年^[1]，由于在测量1992年Landers地震形变场^{[2, 3]}的成功应用而受到广泛关注。但是，该技术易受空间失相干、时间失相干^[4]和大气干扰^[5]等因素的影响，难以进行常态化的实际应用。为了克服D-InSAR技术的这些限制，自20世纪90年代末，一些时间序列InSAR技术被提出。它们总体上可以概括为两类：以永久散射体干涉(permanent scatterer^{[6, 7, 8]}、persistent scatterer interferometry^[9]或PS-InSAR)为代表的单一主影像时间序列InSAR技术和以小基线集技术(small baseline subset interferometry或SBAS InSAR)^{[10, 11, 12]}为代表的多主影像时间序列InSAR技术。时间序列InSAR技术对上述3种制约因素均有良好的免疫力，目前已经取代传统的D-InSAR技术在火山、地震、滑坡等领域得到大量应用^{[13, 14, 15, 16]}。由于城市地区拥有密集的天然点目标，时间序列InSAR技术在城市地面沉降测量方面已经逐步实现工程化应用。国外最典型的例子是欧空局的Terrafirma项目^[17]。在国内，时间序列InSAR技术也被应用于多个城市的地面沉降监测^{[18, 19, 20, 21, 22, 23]}。

但是，时间序列InSAR技术中仍有一些问题有待完善，其中之一就是形变模型。现有的时间序列InSAR处理中，都将形变模式假设为以线性形变为主，然后在后续处理中再恢复出残余的非线性形变分量。事实上，当点目标的形变呈现高度的非线性且在空间不相关的时候，形变反演结果会有很大误差^[7]。目前国内外关于形变模型的研究较少，文献[24—26]曾采用三角函数模型对一类特殊的形变——周期性形变进行了研究^{[24, 25, 26]}。本文将对现有时间序列InSAR处理中的形变模型问题展开研究。首先从原理上分析线性形变模型的不足；然后根据魏尔斯特拉斯逼近定理，提出基于多项式的形变模型，并给出相应的时间序列InSAR形变解算方法；最后介绍多项式形变模型的试验结果及分析。

1 线性形变模型解算的原理分析 1.1 时间序列InSAR相位模型

时间序列InSAR处理中，第i个点目标上的干涉相位(假定地形相位已经通过已有的DEM除去了)可以表示为^{[7, 8]}

式中，wrap(·)表示相位缠绕算子；φ_i^k表示点目标上的干涉相位观测值；φ_i,atmo^k表示大气不均匀引起的干涉相位；φ_i,orb^k为基线误差引起的干涉相位；φ_i,topo-res^k为DEM残差引起的干涉相位；φ_i,noise^k表示失相干噪音。φ_i,topo-res^k与DEM残差Δh_i之间具有线性关系^[8] 式中，B_i,⊥^k是第k个干涉图的垂直基线距，系数a₁由该点的斜距、入射角以及雷达波长决定。地表位移引起的干涉相位φ_i,defo^k可以分解为线性形变相位和非线性形变相位，如下式所示 式中，v_i为线性形变速率；t_i^k为第k个干涉图的影像时间差。这样，干涉相位模型可进一步表示为 式中，非线性形变相位Δφ_{i,,non-linear}^k、基线误差相位φ_i,orb^k和大气相位φ_i,atmo^k都具有较强的空间相关特性，通过两个相邻点目标i、j间的相位差分运算可以基本消除这3项的贡献，此干涉相位差为 式中，Δh_i,j表示两点目标间的高程残差之差；Δv_i,j表示两点目标间的线性形变速率之差异；n^k为相位缠绕因子，残差相位 是相邻点目标间的大气贡献相位、基线误差相位、非线性沉降相位和噪声贡献相位的差分之和。联立全部干涉图(假设共有K个)上i、j点目标间的如式(5)所示的干涉相位差方程可以组成一个方程组，形成时间序列InSAR的干涉相位模型。

1.2 相位模型的求解原理

干涉相位模型是一个不定方程组，没有确定解。求解方法有二维周期图法^[7]和空间搜索法^{[11, 23]}。二者在本质上等价的^[11]，都是基于以下的时态相干因子取得最大值来确定得到了Δh_i,j、Δv_i,j的最佳解

式中，J为虚数单位。这个求解方法的原理是：时态相干因子对应着复平面内K个单位矢量之和，当这些单位矢量有相似的幅角主值时(如图 1(a)所示)，其矢量和有较大的模，极限情况是这些矢量的辐角主值Δφ_i,j,res^k 完全相等，此时γ_i,j取得最大值1；当这些辐角主值完全背离时(如图 1(b)所示)，这些矢量相加可能相互抵消，当完全抵消时γ_i,j取得最小值0。辐角主值相似意味着相位方程解的取值Δh_i,j,Δv_i,j使相位模型与观测值的契合度最佳。

上述求解原理隐含着一个重要的假设，即

只有这个条件成立，Δφ_i,j,res^k才能成为单位矢量exp(JΔφ_i,j,res^k)的幅角主值，而不仅是幅角。若式(8)不满足，如Δφ_i,j,res^k=1.5π，则按照式(6)使干涉相位残差等于1.5π的一组解(Δh_i,j,Δv_i,j)和使干涉相位残差等于－0.5π的另一组解将取得同样的γ_i,j值，此时利用时态相位相干因子最大化得到的模型解将不唯一。因此，式(8)是干涉相位模型能够得到唯一最佳解的一个关键必要条件。同时，从图 1可以看到，相位模型得到唯一最佳解的另一个必要条件是：在所有干涉图上的相邻点目标残差相位Δφ_i,j,res^k的值相近。这个条件在时间序列InSAR处理中通常被当做检验点目标质量的判据，即将相位模型解算出的Δh_i,j、Δv_i,j代回式(6)可以计算出每个干涉图的残差相位Δφ_i,j,res^k，所有干涉图上的残差相位的标准差必须小于某个预设的阈值。如果某一对相邻点目标的残差相位不满足此条件，则可判断其中的一个点目标质量不高，须将此点去掉^[18] 。

一般情况下，式(8)是容易满足的。按式(6)定义的残差相位的4个分量中，由于点目标的特性，噪声相位很小。 由于大气相位和基线误差相位的高度空间相关性，邻近点的Δφ_i,j,atmo^k及Δφ_i,j,orb^k也会很小。 但是，如果真实形变的非线性非常强烈，则相邻点目标之间的非线性形变相位之差可能会较大，从而导致条件式(8)不满足。这就是形变模型的重要性之所在，也是线性形变模型的主要弱点。

2 时间序列InSAR的多项式形变模型 2.1 多项式形变模型的表达

针对上述的以线性为主的形变模型的不足，提出如下多项式形变模型

式中，u_i,p分别对应着N阶多项式形变模型中的线性速率、加速度、三次速度项等；φ_i,res-def^k表示N阶多项式模拟后的残余形变引起的干涉相位。 这个改进的形变模型基于著名的魏尔斯特拉斯逼近定理^[27]：在有界区间上定义的任何连续函数，在此区间内都可以用多项式一致逼近。而地面形变完全可以看成是一个定义于研究时段上的连续函数。因此，从理论上它可以被多项式一致逼近。

相应的，两个邻近点目标i、j之间的干涉相位之差可表示为

式中，残差相位Δφ_i,j,res^k为

由于多项式模型能比线性模型更精确地逼近真实形变，式(9)中的残余形变相位φ_i,res-def^k将会是一个很小的量，因而此时的残差相位将更容易满足式(8)。 图 2显示某水准点观测到的1997—2006年间的地面沉降量(每年1个记录)，以及分别用一次多项式和三次多项式进行拟合的状况。对于这种非线性比较明显的形变，相邻点目标之间的非线性形变相位之差可能较大，从而导致式(8)不成立。而用高次多项式对形变过程进行模拟，由于残余形变相位会很小，则会确保式(8)成立从而得到正确的模型解。

2.2 多项式形变模型的求解

采用多项式形变模型后，干涉相位方程中的形变参数从线性形变速率差扩展为N阶多项式的系数Δu_i,j,1,…,Δu_i,j,N。解算方法也可以采用类似的空间搜索法。首先通过关于形变和DEM精度的先验知识给定待求解参数的变化范围，然后让待求解参数(Δu_i,j,1,…,Δu_i,j,N,Δh_i,j)在各自的变化范围内按照一定的步长进行遍历，对每一组取值计算由式(7)定义的时态相干因子γ_i,j，对应于γ_i,j最大值的取值即为待求参数的最佳解。

计算完相邻点目标之间形变模型参数和高程残差之间的增量Δu_i,j,1,…,Δu_i,j,N,Δh_i,j后，利用一个形变状态和高程残差已知的点作为参考点，通过基于不规则三角网建立的最小二乘增量集成方法^[22]或以时态相干因子作为权重的加权平均法^[11]可以计算得到每一个点目标的多项式形变模型和高程残差的绝对量。然后，需要对残余形变相位和大气相位进行估计。PS-InSAR和SBAS InSAR对此步骤采用的是非常相似的处理，详细处理可参考文献[7, 11, 19]。

3 应用实例分析

以山西省太原市地面沉降为例，进行多项式形变模型的应用研究。所用的SAR数据为2003-08-17—2009-02-01的23景ENVISAT ASAR影像，具体时间和空间基线分布见表 1。 研究区范围为太原市，覆盖面积约为40 km×30 km＝1200 km²(图 3(a))。

太原市是我国最缺水的省会城市之一，大规模的地面沉降开始于20世纪70年代，20世纪80年代为地面沉降快速发展阶段，形成了4个沉降中心，90年代地面沉降持续急剧扩张，进入21世纪，地面沉降呈现出不均衡发展的态势^[28]。

采用小基线集InSAR技术反演太原市地表形变。按照时间间隔小于3年、垂直基线小于400 m的原则，形成了93个小基线干涉图。利用SRTM DEM去除了地形相位。按照93个干涉图的相干系数平均值高于设定阈值(本例中为0.3)的条件进行点目标初选。在相邻点相位模型解算中，通过对时间序列残差相位的均方差及时态相干因子分别设定阈值(本例中设定残差相位均方差小于1.0 rad，时态相干因子大于0.7)，进一步去掉质量不高的点目标。为了对比分析，分别采用三次多项式形变模型和线性模型进行了计算，两种方法进行高相干点相位模型解算的时间分别为296.3 min及92.7 min。图 3(b)、(c)、(d)分别显示了不同阶多项式系数值(分别对应于线性速率、加速度、三次速度项)。将三阶多项式形变和残余形变累加在一起可以得到点目标在整个时段内的平均形变速率或累计形变量。图 4(a)、(b)分别显示了采用两种形变模型得到的总形变量。比较发现，采用线性形变模型，最终获得8111个点目标，而三次模型得到了8755个点目标。如图 4(a)下部黑色虚线方框中包含的点目标(形变量均在400 mm以上)，线性模型方法都没有检测出来。为了进一步说明，将图 4(a)黑色虚线方框区域及其中两个点目标的沉降历史一并显示在图 4(c)中。可以看出，这两个点的形变演化曲线均呈现明显的非线性。在采用线性形变模型的解算中，这些点因其残差相位不满足条件式(8)无法得到正确的模型解而被去掉了。三次形变模型能更精确地模拟这些点的形变，因此能正确解算出这些点的形变量。另一方面，两种模型中绝大部分点重合的事实也表明，研究区内绝大部分点目标的形变状况符合线性模型。

为了进行定量的精度评价，利用从太原市有关部门得到的30个红色水准点(如图 3(a)所示)2004—2006年的实测沉降数据为参考，将两种模型分别得到的同一时段内的形变量(已投影到竖直方向)进行了比较，结果列于表 2。可以看到，仅在5个点(第5、6、13、15、30号点)，三次多项式模型得到的形变量稍不如线性模型得到的形变量准确，而在其他的点目标上，三次多项式模型都更准确，在某些点目标上(如第16、21、26号点)，三次多项式模型更是大幅提高了形变测量的精度。 从这30个水准点上统计，线性模型的形变中误差为6.11 mm，而三次模型的中误差缩小到了3.17 mm。

4 结 论

理论上，高阶多项式总能比低阶的多项式更精确地逼近时间序列观测值。但这不意味着在时间序列InSAR形变监测中，多项式形变模型的阶数越高越好。模型阶数越高，相位模型中需要求解的未知数越多，求解中容易出现解的发散，而且计算速度也会大大降低。在对研究区的形变状况有一定先验知识的情况下，可以根据该区典型的形变曲线确定合适的多项式阶数,阶数等于该形变曲线的拐点数加1。 一般情况下，认为三阶多项式是一个稳健的选择。但是，对于如文献[25]所研究的周期性形变，必须采用正弦或余弦函数模型。在采用更高阶的形变模型时，如何进行模型参数的稳定、快速求解是未来需要进一步研究的问题。