附加系统参数的多时相InSAR时空建模和形变估计

韦建超 李志伟 段梦 曹云梦 冯光财 胡俊

摘要轨道误差和长波大气延迟组成的系统误差是影响InSAR形变监测精度的重要因素之一.传统方法在空间域对干涉图的系统误差建模，容易导致长波形变和系统误差相混淆.本文在时空域利用附加系统参数对系统误差建模，同时根据观测值质量对差分相位观测值定权，采用附加系统参数的加权最小二乘法估计形变参数和系统误差，实现了长波形变和系统误差的分离.模拟实验结果表明，在形变与系统误差的空间变化特性完全一致的极端情况下，本文方法能实现两者的有效分离，估计的形变速率均方根误差比传统方法降低了98.8%.ASAR数据实验显示当形变尺度较小且分散分布时，本文方法和传统方法得到的结果相似;当形变在研究区内表现为长波变化时，本文方法比传统方法估计的形变结果更为稳健.关键词附加系统参数;多时相InSAR;时空建模;形变估计

中图分类号P237

文献标志码A

0引言

合成孔径雷达干涉（InSAR）为研究大范围地表形变提供了强有力的工具，近年来被广泛应用于板块运动[1-2]、地震[3-6] 、地下水抽取[7]等各种自然或人为因素导致的地表沉降研究中.当InSAR用于形变监测时，不同误差源会使测量结果产生偏差，本文主要讨论轨道误差和长波大气延迟情况下的形变估计.如果构成差分干涉图的两个观测时刻轨道参数含有误差，在干涉图中轨道误差会以系统性的残差条纹存在.完全去除干涉图轨道误差条纹要求轨道绝对精度优于1 mm [8]，而目前ERS和Envisat卫星精密轨道数据在法向、沿轨向和交轨向的轨道误差分别为7、24和18 cm [9]，ALOS和Radarsat卫星没有提供精密轨道数据，因此基于卫星定轨数据的轨道误差改正远不能满足测量精度要求.此外，构成干涉图的两个观测时刻大气延迟不一致，导致干涉图中存在附加的大气延迟相位，其主要包含三部分[10]：1）由气压、温度和相对湿度随高度变化引起的地形相关的大气延迟;2）由气压、温度和湿度横向变化而导致的长波（10千米以上）大气延迟;3）由大气边界层湍流导致的小尺度（几千米以下）大气延迟.其中高程相关的大气延迟可以通过高程相关的模型进行建模[11-13]，湍流延迟可通过建立随机模型进行处理[14-15].长波大气延迟在空间上的变化与轨道误差相似，均呈系统性的梯度变化，为InSAR观测的系统性误差.

目前针对InSAR观测的系统误差建模方法主要分为三类：1）频率域建模法，包括通过傅里叶变换[16]和小波分析[17]来改正干涉图的系统误差.该方法难以区分长波形变和系统误差.2）利用卫星轨道状态矢量信息来获取基线参数[18-19].该方法在高相干区域能估计轨道参数并在一定程度上改正轨道误差，但长波大气延迟依然存在.3）空间域建模法.该方法在空间域上建立干涉图相位和空间位置的一阶或高阶函数关系[8，20-21].該方法同样存在长波形变和系统误差不能区分的问题.而多时相InSAR除了能克服传统InSAR技术的失相干问题外，也为系统误差的建模和估计提供了可能[22].本文用附加系统参数对每景SAR影像的系统误差进行建模，对受湍流延迟和失相干噪声影响的观测值定权，用附加系统参数的加权最小二乘方法同时估计系统误差和形变参数.与现有方法相比，本文方法对研究区的形变类型不敏感，在长波形变时也能正确估计系统误差和形变，增加了形变参数估计的精度和可靠性.

1多时域InSAR的时空函数模型

根据合成孔径雷达干涉原理，点P在解缠绕干涉图中的差分相位可表示为

φp=φdef，p+φtopo_res，p+φorb，p+φatm，p+φdec_noi，p，  （1）

公式（1）中φdef，p为形变相位.假设形变的时空模型为

ldef=f（x，y，t|aτ）， （2）

其中x，y表示P的空间位置，t为时间，aτ为形变模型的τ个待求参数向量.因变量x，y可选，当形变在空间上的变化未知时，f（x，y，t|aτ）可简化为f（t|aτ）.在tm时刻的主影像和ts时刻的从影像构成的干涉图中，P的形变相位可表示为

φdef，p=4πλ（f（xp，yp，tm）-f（xp，yp，ts））aτ.  （3）

式（1）中地形残余相位φtopo_res，p可以写成：

φtopo_res，p=4π·Bperp，pλ·R·sinθ·ΔhP， （4）

其中ΔhP为P点的地形改正参数，Bperp，p为干涉图的垂直基线，θ为入射角，R为主影像雷达天线与地面目标之间的距离.式（1）中的φorb，p为轨道误差相位、φatm，p为大气延迟相位.φatm，p包含了地形相关的垂直分层延迟项φatm_ver，p

、长波延迟项φatm_long，p

和小尺度的湍流延迟项φatm_tur，p

[10]：

φatm，p=φatm_ver，p+φatm_long，p+φatm_tur，p，  （5）

其中垂直分层延迟用高程回归模型进行建模和改正[11];湍流延迟和失相干噪声φdec_noi，p用随机模型建模[14，15];而长波大气延迟项φatm_long，p

和轨道误差φorb，p在干涉图上具有系统性的长波变化，构成多时相InSAR的系统误差项.假设各成像时刻系统误差项互相独立，基于轨道误差的非线性特性[23]，各个独立观测时刻的系统误差用以下二维多项式建模：

trend=b1x+b2y+b3xPyP+…. （6）

为避免高阶多项式带来的龙格现象[21]，式（6）最多取至二次项，令m=[xyxyx2y2]，r=[b1b2b3b4b5]T，式（6）写成矩阵形式有：

trend=mr，  （7）

则p在干涉图中的系统误差相位为

φtrend，ms=trend，m-trend，s=m（rm-rs）. （8）

将式（3）、（4）和（8）代入式（1），得到包含形变模型参数、地形残差参数和附加系统参数的观测相位方程：

φ=4π·Bperpλ·R·sinθ·Δh+4πλ（f（x，y，tm）-

f（x，y，ts））aτ+m（rm-rs）+Δp，  （9）

其中Δp主要包含湍流延迟和失相干噪声的随机误差相位.

如果多时相InSAR包含由N景SAR影像组成的M幅干涉图，同时对干涉图中的S个目标点进行分析（如图1所示），则共有M×S个观测方程.写成矩阵形式有：

ΨMS×1=AMS×（S+τ）X（S+τ）×1+BMS×uN YuN×1+WMS×1，  （10）

其中Ψ为差分干涉相位观测值向量;W是包括大气湍流延迟和失相干噪声在内的观测噪声;X为包含地形改正向量Δh和形变模型参数向量aτ的待求参数向量;Y为由各个SAR影像的附加系统参数向量r构成的附加系统参数向量，下标u为式（7）所示的系统误差向量r中参数的个数;A和B分别为待求参数向量和附加系统参数向量的系数矩阵.如果随机模型为：

DΨ=σ20P-1Ψ， （11）

σ0为单位权中误差，PΨ为观测值的权阵.如果仅考虑相干性对观测质量的影响，可以根据相干性估计各差分相位观测值的方差[24]：

σ2i=（1-μ2i）/（2Lμ2i），  （12）

其中μi为i处的相干性，L为相干性估计的多视数，再根据方差定权.如果考虑湍流延迟以及失相干噪声对观测质量的影响，可以建立顾及湍流延迟和失相干噪声的方差协方差[14-15].最后按附加系统参数的加权最小二乘准则平差，则有：

ATPψA

ATPψB

BTPψA

BTPψB

=

ATPψΨ

BTPψΨ. （13）

令N11=ATPΨA，N12=ATPΨB，N21=BTPΨA，N22=BTPΨB，H=N22-N21N-111N12，可以求解形变参数和附加系统参数：

=

N-111+N-111N12H-1N21N-111

-N-111N12M-1

-H-1N21N-111H-1

ATPΨΨ

BTPΨΨ.  （14）

2實例分析

2.1模拟实验

形变在研究区内为长波形变的情形常见于地质构造引起的地表形变，或者研究区位于地震断裂带某一侧时等情形.本文模拟了极端情况下，即形变和系统误差在空间上的变化完全一致，都可以用式（6）表示的长波变化时的形变.在时序上，形变随时间线性变化时的情形.形变的时空模型为

f（x，y，t|aτ）=（a1x+a2y+a3xPyP）t.  （15）

对应的形变相位为

φdef=x（tm-ts）y（tm-ts）xy（tm-ts）·aτ，  （16）

而时序上互相独立的系统误差相位可表示为

φtrend=[…x…-x…y…-y…xy…-xy]·r.  （17）

可见即使形变参数向量aτ和附加系统参数向量r在空间上的变化完全一致，但两者在时空函数模型中的系数矩阵并不相同.当目标点S足够多时，即可同时估计出形变参数和系统误差系数.此外，本文模拟了湍流大气延迟、地形残差、轨道误差及失相干噪声，最后得到包含各种误差的两组完全一样的缠绕干涉图序列（图1）.接着按图2所示的实验流程，分别用空间域建模法（传统方法）和附加系统参数的时空函数模型（本文方法）估计系统误差系数和形变速率参数，并分别计算系统误差系数和形变速率参数的均方根误差（RMSE），作为参数估计精度的评价，在模拟实验中，对所有观测值取相同的权.由于传统方法在空间域建模时只能求出各干涉图的差分系统误差系数，为便于比较，先将本文方法得到的各个SAR影像的系统误差系数按干涉相对进行差分，得到与干涉图对应的差分系统误差系数后再与传统方法的结果比较.图3显示了干涉图在距离向x上的差分系统误差参数估值Δ1，干涉图在方位向y上的差分系统误差参数估值Δ2，以及二次项xy上的差分系统误差参数估值Δ3.本文方法得到的Δ1，Δ2和Δ3的RMSE值比传统方法分别降低了97.54%、92.45%和12.90%.图4显示了沉降参数的估计结果，传统方法得到的形变速率在空间变化趋势上与模拟值明显不符，RMSE为60.92 mm/a;而本文方法得到的形变速率与模拟值更接近，RMSE为0.67 mm/a.

2.2真实数据实验

本文选择中心位置位于39°50′N，116°38′E，面积约为4 000 km2的实验区（图5），该区域平均海拔为43.5 m，地势较为平坦.对覆盖实验区的33景时间跨度从2007年1月3日至2010年9月29日的ASAR影像，利用GAMMA软件进行如下处理：

1）以2009年1月7日获取的影像为公共影像，将其他所有影像进行配准并重采样至公共影像.

2）以382 m为垂直基线阈值，280 d为时间基线阈值对多时相InSAR进行构网，优化网形最终生成61对小基线差分干涉对（图6）.

3）用30 s空间分辨率的SRTM DEM数据消除地形相位;在方位向和距离向分别进行20和4的多视;采用改进的Goldstein滤波方法[25]对干涉图进行滤波，并采用最小费用流法对相干性高于0.7的区域进行解缠，最后得到多时相InSAR的解缠相位序列.

4）对解缠绕的差分相位观测值分别用传统方法和本文方法估计其形变速率参数，同时考虑了湍流延迟和失相干噪声对观测质量的影响，用文献[14]介绍的方法估计观测值的时空方差协方差阵用于参数估计的定权.由于相干性大于0.7的区域绝大部分位于平原地带（图7），且最大高差仅为438 m，因此不考虑随高度变化的大气垂直分层延迟影响.

最终得到的形变速率结果如图7所示，正值（红色）表示从影像相对于主影像在雷达视线方向的上升，负值（绿色）表示从影像相对于主影像下降.两种方法估计的形变速率都处在0～145 mm/a的区间，其中沉降区都主要分布于潮白河、温榆河和泃河流域的冲积和洪积扇平原上，空间分布不均匀.位于昌平、顺义、朝阳和通州的沉降区有逐渐连成一片同时向东扩张的趋势，初步解释为与城市扩张相关.中心城区、海淀区、房山区、丰台区与石景山区的地表沉降速率都在10 mm/a，而且较为稳定;昌平区的沙河镇与上庄镇出现了沉降漏斗，年沉降速率超过70 mm/a;朝阳区金盏乡一带沉降明显，最大平均沉降速率接近80 mm/a.图7c为两种方法估计的形变速率之差，差值在10 mm/a以内，两种方法得到形变速率在量级和空间格局上都非常接近，

本文接着选择位于沉降区一侧的一个局部区域作为局部研究区（图5和图7c中红色矩形区域），假设只有该区域能获得高相干观测值，其在图7a和7b中的形变速率向量分别用向量V1，V2表示.由于V1和V2的差值在6 mm/a以内，本文取=（V1+V2）/2作为该研究区的形变速率参考值（图8a），分别用传统方法和本文方法估计该区域的形变速率向量V′1 和V′2 （图8b、图8c）.可见当以该局部区域作为研究区时，传统方法估计的沉降速率V′1 与参考值存在一个长波变化的偏差.而本文方法得到的V′2 与参考值更为接近.如果用Vi和V′i 的相关性评估形变估计的稳健性，用V′i 的RMSE值来评估两种方法的精度，结果见图8.本文方法比空间域建模法在稳健性方面提高了100%，在RMSE值上降低了59.73%.

3讨论

模拟实验结果表明，当研究区域内的形变在空间上表现为长波缓慢形变时，本文方法得到的结果更加稳健，即使在形变和系统误差在空间上的变化完全一致的极端情况下，本文方法也能实现形变和系统误差的有效分离.而在真实数据实验中，当以整个干涉图覆盖区域作为研究对象时，由于形变相对于研究区域的尺度更小，且分散的分布在研究区的不同位置，

在整个研究区域内不是长波变化，使用传统方法和本文方法得到的形变速率较为接近.而当以标识的局部区域作为研究区时，形变相对于研究区域为大尺度的长波变化，传统方法在估计系统误差参数时明显受到形变的影响，一部分长波形变被误当成系统误差进行了补偿.而本文方法得到的形变速率在空间分布和数值上与参考值更为接近，说明本文方法受形变类型的影响更小，其形变参数的估计结果更稳健.造成这种结果的主要原因如下：

1）空间域建模法在估计系统参数时，只是在空间维度对每个干涉图利用二维多项式拟合，将形变、地形误差项、大气延迟和失相干等分量视为随机误差.而实际上这些量在空间上不一定满足零均值随机分布的特点，尤其是干涉图中的长波形变会对系统误差系数的估计造成干扰，易将长波形变误当成系统误差，造成系统误差的过估计.而附加系统参数的时空函数模型在时空维度上对系统误差和长波形变分别建模，即使两者在空间上的变化完全一致，也能根据两者在时间上的不同特性进行区分，从而实现将长波形变从系统误差中有效分离，避免形变的错误估计.

2）空間域建模法对每个干涉图单独估计系统参数时，假设各干涉图互相独立.而实际上干涉图之间包含有公共影像，不满足独立条件.附加系统参数法对每一个观测时刻的系统误差赋予一组独立的系统参数，差分组合的误差系数由观测时刻的系数计算得到，在数据的解算过程中增加了闭合差约束，因此得到的各干涉图的系统参数符合性更好.

3）空间域建模法没有考虑不同相干性区域差分相位的观测质量对结果的影响，而本文方法可以根据相干性、大气湍流延迟和失相干噪声对观测质量的影响来给相位观测值合理定权.因此，其得到的结果和实际更加相符.

总而言之，随着覆盖范围更小的高分辨率SAR影像（如TerraSAR）数据的应用越来越广，或者在多时相InSAR应用中可能存在可用于分析建模的高相干区域只覆盖局部的形变区域（如本文实例中选的红色矩形区域），此时形变在研究区域内可能呈现出长波的变化趋势.这种情况下传统方法会导致系统误差的过估计和形变的低估计，因此，有必要采用能区分长波形变和系统误差，对形变类型不敏感的本文方法来获得稳健的形变估计结果.

4结论

本文系统介绍了附加系统参数的多时相InSAR时空建模和形变估计方法.该方法给每个观测时刻的轨道误差和长波大气延迟赋予一组独立的附加系统参数建立模型，对形变建立包含形变参数的时空函数模型，在时空维建立多时相InSAR差分相位的观测模型，利用附加系统参数的加权最小二乘平差方法同时估计系统误差系数和形变参数.通过模拟和真实数据实验验证了附加系统参数的多时相InSAR时空函数模型在有效估计系统参数的同时，很好地保留了长波趋势的形变信号.与传统多时相InSAR去轨道误差方法相比，它增加了系统参数闭合差的约束，同时也减少了长波形变对系统参数估计的影响，还考虑了不同观测质量差分相位的权重.虽然本文方法还存在着只针对解缠绕相位、增加计算复杂性等不足，但随着覆盖率更小的高分辨率SAR影像的广泛应用，以及在多时相InSAR应用中可能会出现的较小范围的高相干区域，导致研究区域内形变在空间上呈长波变化的可能性增加，本文的研究增加了利用InSAR估计长波形变的稳健性.

參考文献

References

[1]Champenois J，Fruneau B，Pathier E，et al.Monitoring of active tectonic deformations in the Longitudinal Valley （Eastern Taiwan） using Persistent Scatterer InSAR method with ALOS PALSAR data[J].Earth and Planetary Science Letters，2012，337/338：144-155

[2]Bürgmann R，Hilley G，Ferretti A，et al.Resolving vertical tectonics in the San Francisco Bay Area from permanent scatterer InSAR and GPS analysis[J].Geology，2006，34（3）：221-224

[3]Massonnet D，Rossi M，Carmona C，et al.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature，1993，364（6433）：138-142

[4]Massonnet D，Feigl K，Rossi M，et al.Radar interferometric mapping of deformation in the year after the Landers earthquake[J].Nature，1994，369（6477）：227-230

[5]Zebker H A，Rosen P A，Goldstein R M，et al.On the derivation of coseismic displacement fields using differential radar interferometry：the Landers earthquake[J].Journal of Geophysical Research：Solid Earth，1994，99（B10）：19617-19634

[6]Feng G C，Hetland E A，Ding X L，et al.Coseismic fault slip of the 2008 Mw7.9 Wenchuan earthquake estimated from InSAR and GPS measurements[J].Geophysical Research Letters，2010，37（1）：73-78

[7]Lu Z，Danskin W R.InSAR analysis of natural recharge to define structure of a ground-water basin，San Bernardino，California[J].Geophysical Research Letters，2001，28（13）：2661-2664

[8]Hanssen R F.Radar interferometry：data interpretation and error analysis[M].Springer Science & Business Media，2001

[9]Scharroo R，Visser P.Precise orbit determination and gravity field improvement for the ERS satellites[J].Journal of Geophysical Research：Oceans，1998，103（C4）：8113-8127

[10]Bekaert D P S，Walters R J，Wright T J，et al.Statistical comparison of InSAR tropospheric correction techniques[J].Remote Sensing of Environment，2015，170：40-47

[11]Onn F，Zebker H A.Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a GPS network[J].Journal of Geophysical Research Atmospheres，2006，111（B9）：B09102

[12]Li Z W，Xu W B，Feng G C，et al.Correcting atmospheric effects on InSAR with MERIS water vapour data and elevation-dependent interpolation model[J].Geophysical Journal International，2012，189（2）：898-910

[13]Bekaert D P S，Hooper A，Wright T J.A spatially variable power law tropospheric correction technique for InSAR data[J].Journal of Geophysical Research：Solid Earth，2015，120（2）：1345-1356

[14]Cao Y M，Li Z W，Wei J C，et al.Stochastic modeling for time series InSAR：with emphasis on atmospheric effects[J].Journal of Geodesy，2018，92（2）：185-204

[15]Agram P S，Simons M.A noise model for InSAR time series[J].Journal of Geophysical Research：Solid Earth，2015，120（4）：2752-2771

[16]朱珺，丁曉利，许兵，等.基于条纹频率检测的干涉图去线性趋势精化方法[J].大地测量与地球动力学，2011，31（3）：138-141

ZHU Jun，DING Xiaoli，XU Bing，et al.Refine method for removing the linear trend in interferogram based on fringe frequency test[J].Journal of Geodesy and Geodynamics，2011，31（3）：138-141

[17]Shirzaei M，Walter T R.Estimating the effect of satellite orbital error using wavelet-based robust regression applied to InSAR deformation data[J].IEEE Transactions on Geoscience and Remote Sensing，2011，49（11）：4600-4605

[18]Kohlhase A O，Feigl K L，Massonnet D.Applying differential InSAR to orbital dynamics：a new approach for estimating ERS trajectories[J].Journal of Geodesy，2003，77（9）：493-502

[19]Bhr H，Hanssen R F.Reliable estimation of orbit errors in spaceborne SAR interferometry[J].Journal of Geodesy，2012，86（12）：1147-1164

[20]Ferretti A，Montiguarnieri A，Prati C，et al.InSAR principles：guidelines for SAR interferometry processing and interpretation[J].Journal of Financial Stability，2007，10（10）：156-162

[21]Xu B，Li Z W，Wang Q J，et al.A refined strategy for removing composite errors of SAR interferogram[J].IEEE Geoscience and Remote Sensing Letters，2014，11（1）：143-147

[22]Zhang L，Ding X L，Lu Z，et al.A novel multitemporal InSAR model for joint estimation of deformation rates and orbital errors[J].IEEE Transactions on Geoscience and Remote Sensing，2014，52（6）：3529-3540

[23]Liu G，Hanssen R F，Guo H D，et al.Nonlinear model for InSAR baseline error[J].IEEE Transactions on Geoscience and Remote Sensing，2016，54（9）：5341-5351

[24]Rodriguez E，Martin J M.Theory and design of interferometric synthetic aperture radars[J].IEE Proceedings F Radar and Signal Processing，1992，139（2）：147

[25]Li Z W，Ding X L，Huang C，et al.Improved filtering parameter determination for the Goldstein radar interferogram filter[J].ISPRS Journal of Photogrammetry and Remote Sensing，2008，63（6）：621-634

Spatio-temporal modeling and deformation estimation of multi-temporal

InSAR with additional systematical parameters

WEI Jianchao1，2LI Zhiwei2DUAN Meng2CAO Yunmeng2FENG Guangcai2HU Jun2

1School of Resource & Environment and Safety Engineering，Hunan University of Science and Technology，Xiangtan411201

2School of Geosciences and Info-Physics，Central South University，Changsha410083

AbstractThe systematic error of InSAR is composed of orbit error and long-wavelength atmospheric delay，which is one of the most important factors that affect the accuracy of InSAR deformation monitoring.Traditional method models the systematic error of each individual interferogram in spatial domain，so it is difficult to distinguish between the long-wavelength deformation and systematic error.In this study，we modeled the systematic errors in spatio-temporal domain，determined the weight matrix of the observations according to its quality，and estimated the deformation parameters and systematic error by the weighted least square method.The simulation results show that，even in the extreme cases where the spatial characteristics of deformation and systematic error are completely consistent，the proposed method can separate the deformation and systematic errors effectively.The deformation rate RMSE of the proposed method is 98.8% lower than that of traditional method.The experiment with ASAR shows that，the proposed method and traditional method obtain similar results when the deformation scale is small and scattered over the study area.However，when large scale long-wavelength deformation occurs in the study area，the results obtained by the proposed method are more robust than those obtained by traditional method.

Key wordsadditional systematical parameters;multi-temporal InSAR;spatio-temporal modeling;deformation estimation

收稿日期2019-10-15

資助项目国家自然科学基金（41474007，41404013，41222027）

作者简介韦建超，男，博士，研究方向为合成孔径雷达干涉测量的误差建模与参数反演.wjchao1608@163.com

1湖南科技大学资源环境与安全工程学院，湘潭，411201

2中南大学地球科学与信息物理学院，长沙，410083