CN104459696B - 一种基于平地相位的sar干涉基线精确估计方法 - Google Patents

Abstract

本发明公开了一种基于平地相位的SAR干涉基线精确估计方法，该方法根据SAR干涉基线与平地相位的几何关系，建立严密的SAR干涉基线估计的平差模型，将基线估计的几何结构从起伏变化的地球表面转移到规则的参考椭球面上，并提出一种结合岭估计和截断奇异值分解的迭代未知参数求解算法。该SAR干涉基线精确估计方法是测量平差理论和InSAR理论的有机结合，其包括以下步骤：步骤1：数据预处理；步骤2：系统性信号去除；步骤3：模型建立；步骤4：模型解算。该SAR干涉基线精确估计方法原理直观，易于编程实现及应用，是一种稳健可靠的SAR干涉基线估计算法，可以用于提高InSAR技术测量的精度。

Description

一种基于平地相位的SAR干涉基线精确估计方法

技术领域

本发明涉及涉及合成孔径雷达干涉测量(InSAR，Synthetic Aperture RadarInterferometry)技术，具体涉及一种基于平地相位的SAR干涉基线精确估计方法。

背景技术

经过二十多年的发展，合成孔径雷达干涉测量(InSAR)技术的理论逐渐完善，并在地形测绘和地表形变监测得了广泛的应用。近年来，其应用领域已从地震、火山、城市地表沉降等大范围宏观监测逐步拓展到桥梁、大坝、建筑物等基础设施的精细监测以及城市高精度三维重建，前景十分广阔。

然而，InSAR测量精度受多种误差的影响，如大气延迟误差、DEM误差、随机相位误差以及基线估计误差等。在所有的这些误差中，基线估计误差对InSAR测量精度的影响是系统性的，尤其严重。基线误差不仅影响平地相位的去除，同时会影响地形相位与地表高程的转换系数，对差分InSAR会影响地形相位的去除，而对InSAR技术会直接影响生成DEM产品的精度。目前，基线初始估计方法，受卫星轨道状态矢量精度等因素的影响，难以获取高精度的基线估计结果。现有的基线估计方法包括：

1)基于卫星轨道状态矢量的基线估计方法，利用卫星瞬时坐标的差分可以计算出干涉基线，但是该方法直接受到卫星轨道状态矢量精度的影响。根据误差传播定律可知，该方法估计的基线精度低于卫星轨道状态矢量的精度。如TerraSAR、COSMO-SkyMed等SAR卫星装备了GNSS定位装置，其精密轨道状态矢量的精度约为5cm，而快速轨道的精度仅为10cm；ERS-1/2、ENVISAT等的轨道状态矢量精度约为10～15cm；而ALOS/PALSAR等卫星的轨道状态矢量精度更差。

2)基于影像配准偏移量的基线估计方法，根据主从影像在斜距向的偏移量与斜距差的关系，将距离向的偏移量估计结果转换成斜距差，并结合SAR成像几何条件计算SAR干涉基线。该方法受影像配准算法精度的制约，其基线估计结果精度相对较低，通常在米级或分米级。例如，假设SAR影像配准精度达到1/10像元，则对于3米分辨率的TerraSAR，其斜距差计算精度约为30cm，由此计算得到基线精度较低。

3)基于干涉图条纹率的基线估计方法，其利用InSAR干涉图条纹率与SAR空间儿何的关系确定SAR干涉基线。在实际应用中，由于缺少精确的地形坡度信息，通常选取地形平坦的区域进行条纹频率和基线估计。然而该方法也存在缺陷，一方面，在很多情况下难以找到地形平坦的区域；另一方面，该方法的基线估计精度受FFT计算条纹频率的窗口大小和窗口内干涉相位质量的影响。由于SAR干涉基线在整幅干涉图中不是一个固定值，即基线随时间的变化而变化，基于干涉图条纹率的基线估计方法通过对一定窗口大小(如512×512像素)的干涉图计算平均条纹频率，进而计算出选取区域的平均基线。当条纹频率的估计窗口越大，条纹频率估计越准确，但由此计算得到的基线由于包含了更多的干涉图样本，其与局部真实基线的偏差越大；当条纹频率的估计窗口越小，可以获取与局部真实基线偏差较小的估计结果，但是其估计精度受干涉图质量以及样本容量过小的影响。

4)最小二乘基线精确估计方法，其利用地面控制点(GCP)精化SAR干涉基线的方法，根据已知的解缠相位和地面控制点高程建立观测方程，采用最小二乘方法进行基线参数的求解，通常称之为最小二乘基线精确估计方法。但是很多情况下，我们没法获得研究区域的地面控制点。有研究学者根据控制点为斜距椭球面、多谱勒锥面和地球椭球面的交点这一关系，采用迭代方法计算出地面控制点高程，进而采用上述最小二乘基线估计方法求解基线，然而迭代计算得到的控制点高程的精度有限。

根据以上对现有的InSAR干涉基线估计方法的分析，可以看出这些方法均存在固有的缺陷，无法获取高精度的基线估计结果，导致InSAR干涉测量技术无法达到其应有的理论精度，从而影响了该技术的推广与应用。

因此，有必要设计一种能够获取高精度的、且不依赖于地面控制点的SAR干涉基线精确估计方法。

发明内容

本发明所要解决的技术问题是突破现有SAR干涉基线估计方法所存在的缺陷的约束，提供一种基于平地相位、且无需地面控制点的SAR干涉基线精确估计方法，该基于平地相位且无需地面控制点的SAR干涉基线精确估计方法易于实施，应用广泛，有助于提高InSAR干涉测量技术的精度。

发明的技术解决方案如下：

一种基于平地相位的SAR干涉基线精确估计方法，其特征在于，包括以下步骤：

步骤1：数据预处理：对获取的N景SAR单视复数影像进行配准和重采样，并利用主、从影像的轨道状态矢量信息进行SAR干涉基线初始化估计；采用外部DEM进行差分干涉处理获得差分干涉图，并对差分干涉图进行滤波、解缠，获得解缠的差分干涉相位图φ_diff，unw；N为整数，且N≥2；

步骤2：系统性信号去除：采用先验信息或时序InSAR形变监测方法去除解缠的差分干涉相位图中的形变信号φ_def，同时去除与高程相关的误差信号φ_err，topo，获取残余解缠差分干涉相位图φ_res，unw；

此处的残余解缠差分干涉相位图φ_res，unw包括平地相位残余φ_flat，res、大气延迟相位φ_atm、DEM误差相位φ_topo及相位噪声φ_noise，有

φ_res，unw＝φ_flat，res+φ_atm+φ_topo+φ_noise； (1)

步骤3：模型建立：利用初始化基线估计参数计算初始平地相位φ_flat，init，将φ_flat，init与φ_res，unw相加获得平差模型的观测值；根据φ_res，unw计算观测值的方差，并建立起用于参数平差的函数模型和用于后续参数平差定权的随机模型；【随机模型和函数模型是一起用于后续的参数平差。】

步骤4：模型解算：基于测量平差计算，求得SAR干涉基线的精确估计结果。

所述的步骤3中，参数平差的函数模型由下式表征：

其中L为观测方程的观测向量；X为待求SAR干涉基线的参数向量；

X＝[B_c，0，B_n，0，α_c，α_n，φ₀]^T，参数向量的初始近似值为参数向量的改正数x＝X-X⁰；B_c，0和B_n，0分别为影像中心基线的C和N分量，α_c和α_n分别为基线在C和N方向的变化率，φ₀为参考点的解缠平地相位【参考点通常选取为靠近影像中心的某个相干性较高的像素点】；

为参数向量X的非线性函数；为f(X)的一阶偏导数；其中为主影像的名义斜距(即相对于椭球面的斜距，详见附图1)，为从影像的名义斜距，且 为基线在TCN坐标系下矢量分解：

根据基线的线性变化模型，又可以写成下面的公式：

B_c(t)＝B_c，0+α_c·t

B_n(t)＝B_n，0+α_n·t； (4)

Δ为L的随机误差向量；

以上各个公式中，仅B_c，0，B_n，0，α_c，α_n和φ₀为待求基线参数和参考点平地相位，其余的为已知量或可以由已知量计算得到的中间过程量；

该参数平差的函数模型的误差方程矩阵表达方式为：

式中，V表示观测值改正数向量，x＝X-X⁰为参数改正数向量，l＝L-f(X⁰)，误差方程的设计矩阵；

函数模型对应的随机模型采用方阵D来表征【此处随机模型可以用方阵来表征】，D为对角阵，仅有N个待估计的非零元素，对应于每个观测值L_i(i＝1，2，...，N)的方差有：

公式如下：

式中，为第i个观测值L_i(i＝1，2，...，N)的方差；非对角性元素均为零；

μ为窗口内所有像元值的平均值，即

此处φ_res，unw(h，w)为计算的窗口内的像元值，其从φ_res，unw提取得到。式中2H+1和2W+1为计算局部方差的窗口高度和宽度。

在后续的模型解算中，函数模型确定了基线参数的数学求解表达式，而随机模型D用于设置初始权重矩阵，即权重矩阵为P＝D^-1，D^-1为随机模型D的逆矩阵，通过设置不同质量观测量的权重，以降低误差对函数模型解算的影响。

对所述的参数平差的函数模型采用迭代方法求解，其步骤为：

步骤1：模型初始化：令权阵P＝D^-1为随机模型D的逆矩阵，岭参数k设置为远小于1的非零值，且基线参数向量的初始值其中前四个值为通过轨道状态矢量计算得到的初始基线参数，并设置参考点平地相位初始值为φ₀＝0；

步骤2：对式L＝f(X)+Δ在X⁰处进行泰勒级数展开【L＝f(X)+Δ是非线 性函数，需要进行泰勒级数展开，得到计算观测方程系数矩阵A和向量l，并组成法方程，并计算加权残差平方和SWST₁：

其中，ΔL_i＝L_i-f_i(X⁰)，以及令SWST_old＝SWST₁；

步骤3：在法方程的系数阵对角线所有元素上分别加岭参数k，组成岭估计法方程式，并计算相应法方程系数阵的广义逆；

步骤4：根据岭估计参数平差方法计算参数改正数x的估计值对式处进行泰勒级数展开，计算新的观测方程系数矩阵A′和向量l′，并根据式计算加权残差平方和SWST_new；

步骤5：比较SWST_nww和SWST_old的大小，按以下方法判断是否接受本次迭代所得参数向量改正数，并输出加权残差平方和

如果SWST_new＜SWST_old，则由式计算更新以及更新设计矩阵A＝A′和向量l＝l′，且调整岭估计参数，同时根据残差更新权矩阵P，以及设置SWST_old＝SWST_new；否则，保持X⁰，A，l和P不变，将调整岭参数；

步骤6：判断迭代终止条件：如果迭代次数为1，则定义令并直接跳转到步骤3进行第二次迭代；如果迭代次数大于1，则判断迭代终止条件是否满足，具体终止条件设置为

且 (9)

连续满足两次【即第n次和第n+1次迭代过程中均满足公式(9)的要求】，或者迭代次数达到设定的最大次数(如20次)；如果终止条件满足，则输出参数向量估计结果否则，令并跳转到步骤3继续迭代；

通过上述的迭代参数求解过程，基线参数向量的估计结果已精确获得，其包含了基线模型的参数估计值和即SAR干涉基线的精确估计结果，这四个基线模型参数能用于后续的高精度InSAR数据处理。

定义：残余解缠差分干涉相位图(φ_res，unw)包括平地相位残余φ_flat，res、大气延迟φ_atm、DEM误差φ_topo及相位噪声φ_noise，即

所述步骤1包括如下儿个子步骤：

步骤a：最优主影像选取及从影像配准到最优主影像

首先，根据N景(N≥2)SAR单视复数影像的时间和空间基线分布情况，选择用于配准的最优参考主影像，具体的最优准则为

函数f(*，*)的定义为

其中，m为选取的最优主影像编号，s_m为m被选为主影像时所有从影像的编号， 和分别表示主影像为m、从影像为s_m对应的空间垂直基线和时间基线，B_c和T_c分别为空间和时间临界基线，α和β分别为空间基线和时间基线对最优主影像选择的指数因子。

然后，基于灰度匹配算法将所有从影像与最优主影像配准，并根据配准转换关系将所有从影像重采样到主影像成像空间。

步骤b：InSAR数据差分处理及差分干涉图滤波与解缠

首先，将上述配准到同一成像空间的N景(N≥2)SAR影像进行短基线组合， 得到相干性较好的M个InSAR干涉对(N-1≤M≤(N-1)(N-2)/2)。

然后，采用外部DEM对上述M个InSAR干涉对进行差分干涉处理，并对得到的M个差分干涉图进行滤波以及相位解缠，从而得到M个解缠差分干涉图。

所述的步骤2中：对于差分干涉图中可能存在的形变信号，可以采用以下两种可用策略对其去除：

1)利用研究区域的先验信息，如地面实测数据等将干涉图中的形变信号去除；

2)采用时序InSAR形变监测技术(如PS-InSAR、SBAS-InSAR技术)结合数据本身解算出干涉图中的形变信号，并将其去除。

对于差分干涉图中的与高程相关的误差信号，包括与高程相关的大气延迟以及由于基线存在误差导致的地形残余，采用线性回归分析方法，从差分干涉图中分离出与高程相关的误差信号。

所述步骤3包括如下儿个子步骤：

步骤a：建立SAR干涉基线与平地相位之间的严密数学关系式

如附图1示，SAR干涉基线定义为卫星两次对同一目标成像时，其位置构成的空间矢量，即有

式中，为基线矢量在空间直角坐标系的表达，和 为卫星的三维坐标，T为向量转置操作。而在实际的软件编程实现，通常采用基线在TCN坐标系下的投影矢量为了便于推导，我们先给出TCN坐标系的定义。如附图2示，卫星平台在空间直角坐标系中的三维坐标矢量为，则卫星平台的定义：T轴为卫星飞行速度的切线方向，N轴为卫星平台位置与坐标原点O的连线方向(指向坐标原点)，C轴为垂直于T轴和N轴且构成右手坐标系的方向。根据上述定义，可以得到空间直角坐标系XYZ到TCN坐标系的转换矩阵为

式中，向量和分别为T轴、C轴和N轴在空间直角坐标系中的归一化矢量。则与的转换关系为

由于SAR干涉基线在整幅干涉图中不是一个固定值，因此，我们采用常用的线性模型对基线分量进行建模，即

B_c(t)＝B_c，0+α_c·t (16)

B_n(t)＝B_n，0+α_n·t

式中，B_c，0和B_n，0为别为参考时刻SAR干涉基线模型在C和N方向的分量，即基线的常量，通常设置影像中心时间为参考时间点(则此时B_c，0和B_n，0分别为影像中心基线的C和N分量)，α_c和α_n分别为基线在C和N方向的变化率，t为影像任意行相对参考时刻的时间累积量。待求模型参数为B_c，0，B_n，0，α_c和α_n。此处，考虑到经过影像配准，基线在T方向的分量已经被消除，因此令B_t(t)＝0，且在后续的基线精化中不再考虑。即有

根据附图1所示的InSAR成像儿何关系，准确的解缠平地相位可以表示为

式中，λ为雷达波长，为斜距在TCN坐标系中的矢量分解，φ₀为参考点的解缠平地相位。根据余弦定理有

且

上述(18)-(20)即为SAR干涉基线与平地相位之间的严密数学关系式，通过求解这一关系式的未知参数，即可以得到SAR干涉基线的精确估计结果。式中，B_c，n，B_n，0，α_c，α_n和φ₀为待求基线参数和参考点平地相位，其余的为已知量 或可以由已知量计算得到的中间过程量。

步骤b：确定解缠的平地相位

从步骤b可以知道，只要解缠的平地相位φ_{flat，unw，true}可以获得，则SAR干涉基线可严密解算出来。根据步骤1的分析，我们已经获得了平地相位的残余部分φ_res，unw，并利用初始化基线估计初始平地相位φ_flat，init，将其与φ_res，unw相加获得真实平地相位的观测值，即

φ_flat，unw＝φ_{flat，unw，true}+φ_reeor (21)

式中，φ_error＝φ_atm+φ_topo+φ_noise为观测值φ_flat，unw中含有的误差相位。

步骤c：建立SAR干涉基线精确估计的平差函数模型和随机模型

假定有N个观测值，令观测方程的观测向量L＝φ_flat，unw，其维数为N×1，待求SAR干涉基线参数向量X＝[B_c，0，B_n，0，α_c，α_n，φ₀]^T，以及 为参数向量X的非线性函数，则式(18)可改写为

L＝f(X)+Δ (22)

其中，Δ为观测量L的随机误差向量。式(22)为非线性方程，在本文中我们将其进行线性化转为线性方程求解。在近似值处对式(22)进行泰勒级数展开，保留至一次项，即得SAR干涉基线参数估计的函数模型

将式(23)改写为误差方程并采用矩阵形式表达，有

式中，V观测值改正数向量，参数改正数向量x＝X-X⁰，l＝L-f(X⁰)，A为非线性方程f(X)的一阶偏导数构成的矩阵.

关于与函数模型(式23)对应的随机模型的建立，即观测向量L的先验方差阵D(L)的估计(为了便于公式推导，将D(L)简写为D)，在此假设观测向量L之间是相互独立的，则D为对角阵，仅有N个待估计的非零元素，对应于每个观测值L_i(i＝1，2，...，N)的方差此时有

对于的估计，利用去除系统性信号之后得到的残余解缠差分干涉相位图φ_res，unw进行估计，具体的实现过程：以L_i所处的行列位置为中心，取一个宽度和高度分别为2W+1和2H+1的窗口，W、H为大于1的正整数，用选取的窗口范围内的所有像元作为样本计算其方差作为的估计，即

式中，μ为窗口内所有像元值的平均值，即

由于φ_res，unw中含有长波长的轨道误差相位，这一趋势相位的存在导致φ_res，unw不满足平稳条件假设，在此，我们采用局部平面多项式拟合方法获取这一趋势相位，并将其从φ_res，unw中去除，然后采用去除趋势相位之后的相位值代入式(26)和(27)计算

所述步骤4基于测量平差计算，求解SAR干涉基线的具体实现如下：

根据上述步骤3，我们已得到了平差的函数模型和随机模型，根据最小二乘原理，在

V^TPV＝min (28)下，可以得到法方程为

A^TPAx＝A^TPl (29)

式中，为观测向量L的权矩阵，当取单位权方差时，P＝D^-1，即P为D的逆矩阵D^-1。

从式(29)可以看出，由于待求参数个数为5，所以需要至少5个观测量，才能求解法方程。此时，x的最优线性无偏估计结果为

从而可求得基线参数向量X最优线性无偏估计

通过分析设计矩阵A，发现基线参数向量的初始值X⁰非常接近X时，则A的前四列中每个元素的分子将会比较小，加上除以分母中的斜距，前四列会高度相关，A趋向于复共线，从而导致法方程病态，此时最小二乘的解将不再是最优解。因此，引入岭估计解决法方程病态问题，法方程改写为

(A^TPA+KI)x＝A^TPl (32)

其中，k为岭参数，为大于0的常数，I为单位阵。应当注意到，引入岭估计可以改善法方程(式30)的病态条件，但是如果k值过小，矩阵A^TPA+KI仍可能会存在奇异的情形。在此，我们引入截断奇异值分解(TSVD)求解上述岭估计的法方程。令矩阵A^TPA+KI的SVD分解为

A^TPA+kI＝USV^T (33)

其中，U和V为大小5×5的酉矩阵，对角阵S＝Diag(s₁，s₂，s₃，s₄，s₅)的对角元素为相应的奇异值，且s₁≥s₂≥s₃≥s₄≥s₅≥0。TSVD的阂值设置为

thres＝T·max{s_i}i＝1，2，3，4，5 (34)

式中，T为阂值比例系数(如10^-6)，max{s_i}为求取最大值操作。在条件s_r＞thres的限制下，矩阵S的广义逆为

矩阵A^TpA+KI的广义逆为

因此，x的估计结果为

将式(37)代入式(31)，即可以求得参数向量X的估计值

根据前面的分析，在准确给定权阵P和岭参数k的情况下，可以很容易求解上述SAR干涉基线估计平差模型，并得到基线参数向量的估计值而在本文中，权阵P和岭参数k无法一开始就准确给出，同时，方程式(24)的线性化也会存在误差，因此，我们采用迭代方法求解上述基线估计模型，其具体步骤为：

Step 1：模型初始化：令权阵P＝D^-1为随机模型D的逆矩阵，岭参数k设置为远小于1的非零值，且基线参数向量的初始值其中前四个值为通过轨道状态矢量计算得到的初始基线参数，并设置参考点平地相位初始值为φ₀＝0；

Step 2：对式L＝f(X)+Δ在X⁰处进行泰勒级数展开【L＝f(X)+Δ是非线性函数，需要进行泰勒级数展开，得到】，计算观测方程系数矩阵A和向量l，并组成法方程，并计算加权残差平方和SWST₁：

其中，ΔL_i＝L_i-f_i(X⁰)，以及令SWST_old＝SWST₁；

Step 3：在法方程的系数阵对角线所有元素上分别加岭参数k，组成岭估计法方程式，并计算相应法方程系数阵的广义逆；

Step 4：根据岭估计参数平差方法计算参数改正数x的估计值对式处进行泰勒级数展开，计算新的观测方程系数矩阵A′和向 量l′，并根据式计算加权残差平方和SWST_new；

Step 5：比较SWST_new和SWST_old的大小，按以下方法判断是否接受本次迭代所得参数向量改正数并输出加权残差平方和

如果SWST_new＜SWST_old，则由式计算更新以及更新设计矩阵A＝A′和向量l＝l′，且调整岭估计参数，同时根据残差更新权矩阵P，以及设置SWST_old＝SWST_new；否则，保持X⁰，A，l和P不变，将调整岭参数；

Step 6：判断迭代终止条件：如果迭代次数为1，则定义令并直接跳转到步骤3进行第二次迭代；如果迭代次数大于1，则判断迭代终止条件是否满足，具体终止条件设置为

且 (38)

连续满足两次【即第n次和第n+1次迭代过程中均满足公式(38)的要求】，或者迭代次数达到设定的最大次数(如20次)；如果终止条件满足，则输出参数向量估计结果否则，令并跳转到步骤3继续迭代；

通过上述的迭代参数求解过程，基线参数向量的估计结果【说明：X包含了所需求估计的四个基线参数B_c，0，B_n，0，α_c和α_n，同时包括一个常量φ₀，即X＝[B_c，n，B_n，0，α_c，α_n，φ₀]^T。在测量平差中，符号X表示真值，表示真值的估计结果，同时也在实施例中的估计结果列于表2中】已精确获得，其包含了基线模型的参数估计值和即SAR干涉基线的精确估计结果，这四个基线模型参数将可以用于后续的高精度InSAR数据处理。

有益效果：

本发明的基于平地相位的SAR干涉基线精确估计方法，该方法根据SAR 干涉基线与平地相位的儿何关系，建立了严密的SAR干涉基线精确估计的数学方程式，将基线估计的儿何结构从起伏变化的地球表面转移到规则的参考椭球面上，解决了SAR干涉基线精确估计依赖于地面控制点的问题，有效地避免了地形起伏对基线估计的影响。

在所述的步骤1中，对N景(N≥2)SAR单视复数影像的配准，根据影像集的时间和空间基线分布情况，采用最优准则选择影像配准的最优参考主影像，有利于SAR影像的精确配准。

在所述的步骤2中，采用先验信息或时序InSAR形变监测方法去除解缠的差分干涉相位图中的形变信号φ_def，有利于避免形变信号对SAR干涉基线估计的影响。

在所述的步骤2中，在去除形变信号的基础上，采用线性回归分析方法去除与高程相关的误差信号φ_err，topo，获取残余解缠差分干涉相位图，有利于避免与高程相关的系统性信号对SAR干涉基线估计的影响。

在所述的步骤3中，采用泰勒级数对SAR干涉基线精确估计的非线性数学方程式在近似值处展开，建立参数平差的函数模型，将非线性方程组转为线性方程组求解，有利于降低参数求解的复杂性。

在所述的步骤3中，采用局部统计方法计算观测量的方差，进而建立参数平差的随机模型；顾及观测量非平稳特性，采用局部平面多项式拟合方法，去除观测量的非平稳分量的趋势项，使随机模型的建立更严密。

在所述的步骤4中，参数的测量平差计算同时采用了岭估计方法和截断SVD分解求逆方法，有效解决了平差模型法方程系数矩阵病态甚至秩亏问题，有利于提高SAR干涉基线参数估计的可靠性。

在所述的步骤4中，采用迭代思想并结合岭估计和和截断SVD分解求逆方法求解SAR干涉基线模型参数，在迭代过程中岭估计的作用在于以有偏估计代替无偏估计，以加速迭代收敛，并随着迭代次数的增加，基线估计结果逐渐趋近于真值，岭参数逐渐比例缩小，逐渐弱化有偏估计的作用；而截断SVD分解方法用于求法方程系数阵的逆矩阵。当法方程系数阵秩亏时，截断SVD分解方法计算得到的是法方程系数阵的广义逆；当方法方程系数病态时，通过 截断SVD分解方法求逆，有利于提高参数解算的可靠性。

所述的步骤4的迭代求解参数的过程中，以观测值残差的绝对值更新观测值权的权重。与较大的残差绝对值对应的观测量，其在下轮迭代的的权重将会减小；与较小的残差绝对值对应的观测量，其在下轮迭代的权重将会增大，该处理思想有利于减弱和消除观测值中含有粗差对基线参数解算的影响。

在所述的步骤4中，采用加权残差平方和作为迭代收敛指标，若本轮迭代得到的加权残差平方和比上一轮的小，则本次迭代收敛，若本轮迭代发散，则通过调整岭估计参数重新求解基线参数向量改正数，并进入下轮迭代；同时采用前后两次迭代对应的加权残差平方和变化及相对变化量作为迭代终止条件，本轮迭代的加权残差平方和比上一轮迭代的加权残差平方和小，且相对变化量小于设定的阂值这两个条件连续满足两次，则终止迭代。采用这些指标比采用参数变化量作为迭代收敛和终止条件，在所求解的参数具有不同量纲更合理，且易于编程实现。

本发明的基于平地相位的SAR干涉基线精确估计方法是一种基于平地相位且无需地面控制点的SAR干涉基线精确估计方法，该方法根据SAR干涉基线与平地相位的儿何关系，建立了严密的SAR干涉基线精确估计的平差模型，将基线估计的儿何结构从起伏变化的地球表面转移到规则的参考椭球面上，解决了SAR干涉基线精确估计依赖于高面控制点的问题，且有效地避免了地形起伏到基线估计的影响；该发明提出了一种结合岭估计和截断奇异值分解(TSVD，Truncated Singular Value Decomposition)的迭代参数估计算法，解决了平差模型中法方程系数阵奇异和病态引起参数估计不可靠地、不稳定的问题。该SAR干涉基线精确估计方法原理直观，易于编程实现及应用扩展，是一种稳健可靠的SAR基线估计算法，可以提高InSAR技术测量的精度。

附图说明

图1为InSAR成像儿何条件示意图

图2为卫星平台TCN坐标系示意图

图3为基于平地相位且无需地面控制点的SAR干涉基线精确估计算法流程图；

图4为实施案例中步骤3数据处理的中间结果图。(a)编码到SAR坐标系下的DEM；(b)滤波后的差分干涉图；(c)滤波后的相干图；(d)解缠的差分干涉图

图5为实施案例中步骤4数据处理的中间结果图。(a)与高程相关的误差信号；(b)去除系统性信号的残余解缠差分干涉相位图；(c)去除系统性信号的残余(缠绕)差分干涉相位图

图6为实施案例中步骤5得到的解缠平地相位图，其作为平差模型的观测值用于后续的平差计算。

图7为实施案例中步骤8采用新的基线估计结果进行重新差分干涉图处理得到的差分干涉图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

如图1-7，一种基于平地相位的SAR干涉基线精确估计方法，其特征在于，包括以下步骤：

步骤1：数据预处理：对获取的N景SAR单视复数影像进行配准和重采样，并利用主、从影像的轨道状态矢量信息进行SAR干涉基线初始化估计；采用外部DEM进行差分干涉处理获得差分干涉图，并对差分干涉图进行滤波、解缠，获得解缠的差分干涉相位图φ_diff，unw；N为整数，且N≥2；

步骤2：系统性信号去除：采用先验信息或时序InSAR形变监测方法去除解缠的差分干涉相位图中的形变信号φ_def，同时去除与高程相关的误差信号φ_err，topo，获取残余解缠差分干涉相位图φ_res，unw；

此处的残余解缠差分干涉相位图φ_res，unw包括平地相位残余φ_flat，res、大气延迟相位φ_atm、DEM误差相位φ_topo及相位噪声φ_noise，有

φ_res，unw＝φ_flat，res+φ_atm+φ_topo+φ_noise； (1)

步骤3：模型建立：利用初始化基线估计参数计算初始平地相位φ_flat，init，将φ_flat，init与φ_res，unw相加获得平差模型的观测值；根据φ_res，unw计算观测值的方差，并建立起用于参数平差的函数模型和用于后续参数平差定权的随机模型；【随机模型和函数模型是一起用于后续的参数平差。】

步骤4：模型解算：基于测量平差计算，求得SAR干涉基线的精确估计结果。

所述的步骤3中，参数平差的函数模型由下式表征：

其中L为观测方程的观测向量；X为待求SAR干涉基线的参数向量；

X＝[B_c，0，B_n，0，α_c，α_n，φ₀]^T，参数向量的初始近似值为参数向量的改正数x＝X-X⁰；B_c，0和B_n，0分别为影像中心基线的C和N分量，α_c和α_n分别为基线在C和N方向的变化率，φ₀为参考点的解缠平地相位【参考点通常选取为靠近影像中心的某个相干性较高的像素点】；

为参数向量X的非线性函数；为f(X)的一阶偏导数；其中为主影像的名义斜距(即相对于椭球面的斜距，详见附图1)，为从影像的名义斜距，且 为基线在TCN坐标系下矢量分解：

根据基线的线性变化模型，又可以写成下面的公式：

B_c(t)＝B_c，0+α_c·t

B_n(t)＝B_n，0+α_n·t； (4)

Δ为L的随机误差向量；

以上各个公式中，仅B_c，0，B_n，0，α_c，α_n和φ₀为待求基线参数和参考点平地相位，其余的为已知量或可以由已知量计算得到的中间过程量；

该参数平差的函数模型的误差方程矩阵表达方式为：

式中，V表示观测值改正数向量，x＝X-X0为参数改正数向量，l＝L-f(X0)，误差方程的设计矩阵；

函数模型对应的随机模型采用方阵D来表征【此处随机模型可以用方阵来表征】，D为对角阵，仅有N个待估计的非零元素，对应于每个观测值L_i(i＝1，2，...，N)的方差有：

公式如下：

式中，为第i个观测值Li(i＝1，2，...，N)的方差；非对角性元素均为零；

μ为窗口内所有像元值的平均值，即

此处φ_res，unw(h，w)为计算的窗口内的像元值，其从φ_res，unw提取得到。

式中2H+1和2W+1为计算局部方差的窗口高度和宽度。

在后续的模型解算中，函数模型确定了基线参数的数学求解表达式，而随机模型D用于设置初始权重矩阵，即权重矩阵为P＝D^-1，D^-1为随机模型D的 逆矩阵，通过设置不同质量观测量的权重，以降低误差对函数模型解算的影响。

对所述的参数平差的函数模型采用迭代方法求解，其步骤为：

步骤1：模型初始化：令权阵P＝D^-1为随机模型D的逆矩阵，岭参数k设置为远小于1的非零值，且基线参数向量的初始值其中前四个值为通过轨道状态矢量计算得到的初始基线参数，并设置参考点平地相位初始值为φ₀＝0；

步骤2：对式L＝f(X)+Δ在X⁰处进行泰勒级数展开【L＝f(X)+Δ是非线性函数，需要进行泰勒级数展开，得到】，计算观测方程系数矩阵A和向量l，并组成法方程，并计算加权残差平方和SWST₁：

其中，ΔL_i＝L_i-f_i(X⁰)，以及令SWST_old＝SWST₁；

步骤3：在法方程的系数阵对角线所有元素上分别加岭参数k，组成岭估计法方程式，并计算相应法方程系数阵的广义逆；

步骤4：根据岭估计参数平差方法计算参数改正数x的估计值对式处进行泰勒级数展开，计算新的观测方程系数矩阵A′和向量l′，并根据式计算加权残差平方和SWST_new；

步骤5：比较SWST_new和SWST_old的大小，按以下方法判断是否接受本次迭代所得参数向量改正数并输出加权残差平方和

如果SWST_new＜SWST_old，则由式计算更新以及更新设计矩阵A＝A′和向量l＝l′，且调整岭估计参数，同时根据残差更新权矩阵P，以及设置SWST_old＝SWST_new；否则，保持X⁰，A，l和P不变，将调整岭参数；

步骤6：判断迭代终止条件：如果迭代次数为1，则定义令并直接跳转到步骤3进行第二次迭代；如果迭代次数大于1，则判断迭代终止条件是否满足，具体终止条件设置为

且 (9)

连续满足两次【即第n次和第n+1次迭代过程中均满足公式(9)的要求】，或者迭代次数达到设定的最大次数(如20次)；如果终止条件满足，则输出参数向量估计结果否则，令并跳转到步骤3继续迭代；

通过上述的迭代参数求解过程，基线参数向量的估计结果已精确获得，其包含了基线模型的参数估计值和即SAR干涉基线的精确估计结果，这四个基线模型参数能用于后续的高精度InSAR数据处理。

为了使本发明的目的、技术方案及其优点更加清楚明白，利用覆盖2010年舟曲泥石流的两景PALSAR数据，具体参数见附表1，对本发明进行深入详细地说明。应当理解，此处所描述的具体实施案例仅用以解释本发明，并不用于限定本发明。

本发明具体实施案例以单个干涉对来阐述，对于多个干涉对的实施方式与单个干涉对的区别仅在于利用时序InSAR形变监测方法去除所有干涉对中的形变信号，这已有很多相关资料可以参考，该部分不是本发明所涉及的难点及重点，为了简洁起见，实施案例以单个干涉图开展。

基于上述所采用的舟曲的PALSAR数据，本发明的具体实施步骤如下：

步骤1：利用SAR影像灰度配准技术，将从影像与主影像进行配准，配准精度优于1/10像素，并将从影像重采样到主影像成像空间，为了降低重采样过程对信号的损失，具体的重采样算法的核函数使用二维sinc函数，用于采样计算的点数设置为11×11个；

步骤2：利用主影像和从影像的轨道状态矢量信息进行初始化基线估计，获得基线参数的初始化估计结果，即和列于附表2中；

步骤3：利用步骤2得到的初始化基线参数，以及结合外部DEM对主从影像进行差分干涉数据处理。

处理过程中，对InSAR原始干涉图进行多视，使其空间分辨率变为30米；外部DEM采用的是美国NASA发布的v4.1版本SRTM(空间分辨率为90米)，并将其过采样到30米空间分辨率，且编码到SAR坐标系使其与原始干涉图的分辨率和维数均一致，附图4(a)显示了编码到SAR坐标系下的DEM。

采用Goldstein滤波算法对差分干涉图进行二维频率域滤波，滤波窗口设置为32×32像素，滤波因子α＝0.65，从而得到滤波的差分干涉图，见附图4(b)；

采用最小网络费用流方法对滤波后的差分干涉图进行二维相位解缠，解缠所采用的加权因子为滤波后的相干图，见附图4(c)。为了避免低相干区域对解缠的影响，在解缠的过程中，去除相干性低于0.4的像元点，最后获得解缠的差分干涉图φ_diff，unw，见附图4(d)；

步骤4：根据已获得的研究区域的形变先验信息，发生滑坡、泥石流的区域较小(见附图4(b)的黑色椭圆圈出的区域)，且在主、从影像获取的时间间隔内没有明显的大范围地表形变活动，可以认为图中波长较长的条纹是由基线误差引起的平地相位未去除干净而产生的，即这些条纹为平地相位残余。因此，在本实施案例中，人工的扣除附图4(b)所示的黑色椭圆范围即可以排除形变信号对SAR干涉基线估计的影响；对于与高程相关的误差信号φ_err，topo的去除是在扣除形变信号的基础上进行，将解缠的相位与地面高程进行回归分析，求得回归方程系数，然后根据回归方程及所求得的回归方程系数计算出与高程相关的误差信号，见附图5(a)，并将这一误差信号从解缠的相位图中去除，从而得到残余解缠差分干涉相位图φ_res，unw，见图5(b)；

步骤5：利用初始化基线估计参数计算初始平地相位φ_flat，init，将其与φ_res，unw相加获得平差模型的观测值，即有φ_unw，unw＝φ_flat，init+φ_res，unw，见附图6；

步骤6：根据公式(22)到(27)建立基线估计平差模型的函数模型和随机模型；

步骤7：采用前面所述的迭代参数求解策略，并基于测量平差计算，求解出基线参数向量X的估计结果即

从而得到基线模型参数的估计结果，即和列于附表2中。

步骤8：采用步骤7得到的基线精确估计结果进行重新差分干涉图，得到新的差分干涉图显示于附图7中。

从附图7中可以看出，图4(b)滤波的差分干涉图中的长波长轨道误差条纹已经完全被去除，同时潜在的滑坡区域的形变也突显出来(见附图7中黑色椭圆范围)，这说明了基线估计结果已经足够准确，也验证了本发明的SAR干涉基线精确估计方法的有效性。

表1 舟曲实验区PALSAR数据具体参数

表2 SAR干涉基线初始估计和精确估计结果

Claims (3)

1.一种基于平地相位的SAR干涉基线精确估计方法，其特征在于，包括以下步骤：

步骤1：数据预处理：对获取的N景SAR单视复数影像进行配准和重采样，并利用主、从影像的轨道状态矢量信息进行SAR干涉基线初始化估计；采用外部DEM进行差分干涉处理获得差分干涉图，并对差分干涉图进行滤波、解缠，获得解缠的差分干涉相位图φ_diff,unw；N为整数，且N≥2；

步骤2：系统性信号去除：采用先验信息或时序InSAR形变监测方法去除解缠的差分干涉相位图中的形变信号φ_def，同时去除与高程相关的误差信号φ_err,topo，获取残余解缠差分干涉相位图φ_res,unw；

此处的残余解缠差分干涉相位图φ_res,unw包括平地相位残余φ_flat,res、大气延迟相位φ_atm、DEM误差相位φ_topo及相位噪声φ_noise，有

φ_res,unw＝φ_flat,res+φ_atm+φ_topo+φ_noise； (1)

步骤3：模型建立：利用初始化基线估计参数计算初始平地相位φ_flat,init，将φ_flat,init与φ_res,unw相加获得平差模型的观测值；根据φ_res,unw计算观测值的方差，并建立起用于参数平差的函数模型和用于后续参数平差定权的随机模型；

步骤4：模型解算：基于测量平差计算，求得SAR干涉基线的精确估计结果。

2.根据权利要求1所述的基于平地相位的SAR干涉基线精确估计方法，其特征在于，所述的步骤3中，参数平差的函数模型由下式表征：

$L=f\left(X^0\right)+\frac{\∂f\left(X\right)}{\∂X}{|}_{X=X^0}\·x+\mathrm{\Δ};---\left(2\right)$

其中L为观测方程的观测向量；X为待求SAR干涉基线的参数向量；

X＝[B_c,0,B_n,0,α_c,α_n,φ₀]^T，参数向量的初始近似值为参数向量的改正数x＝X-X⁰；B_c,0和B_n,0分别为影像中心基线的C和N分量，α_c和α_n分别为基线在C和N方向的变化率，φ₀为参考点的解缠平地相位；

为参数向量X的非线性函数；为f(X)的一阶偏导数；其中为主影像的名义斜距，为从影像的名义斜距，且 为基线在TCN坐标系下矢量分解：

$\begin{array}{c}\stackrel{\→}{B}{\left(t\right)}_{t,c,n}={\[0,B_c\left(t\right),B_n\left(t\right)\]}^T\\ =0\·\stackrel{\→}{t}+B_c\left(t\right)\·\stackrel{\→}{c}+B_n\left(t\right)\·\stackrel{\→}{n}\end{array};---\left(3\right)$

根据基线的线性变化模型，又可以写成下面的公式：

B_c(t)＝B_c,0+α_c·t

B_n(t)＝B_n,0+α_n·t； (4)

Δ为L的随机误差向量；

以上各个公式中，仅B_c,0,B_n,0,α_c,α_n和φ₀为待求基线参数和参考点平地相位，其余的为已知量或可以由已知量计算得到的中间过程量；

该参数平差的函数模型的误差方程矩阵表达方式为：

$\begin{array}{c}V=\frac{\∂f\left(X\right)}{\∂X}{|}_{X=X^0}\·x-l\\ =Ax-l;\end{array}---\left(5\right)$

式中，V表示观测值改正数向量，x＝X-X⁰为参数改正数向量，l＝L-f(X⁰)，误差方程的设计矩阵；

函数模型对应的随机模型采用方阵D来表征，D为对角阵，仅有N个待估计的非零元素，对应于每个观测值L_i(i＝1,2,…,N)的方差有：

公式如下：

${\mathrm{\σ}}_{L_i}^2=\frac{1}{(2H+1)\×(2W+1)-1}\underset{h=-H}{\overset{H}{\Σ}}\underset{w=-W}{\overset{W}{\Σ}}{\left({\mathrm{\φ}}_{res,unw}\right(h,w)-\mathrm{\μ})}^2;---\left(7\right)$

式中，为第i个观测值L_i(i＝1,2,…,N)的方差；非对角性元素均为零；

μ为窗口内所有像元值的平均值，即

式中2H+1和2W+1为计算局部方差的窗口高度和宽度。

3.根据权利要求2所述的基于平地相位的SAR干涉基线精确估计方法，其特征在于，对所述的参数平差的函数模型采用迭代方法求解，其步骤为：

步骤1:模型初始化：令权阵P＝D^-1为随机模型D的逆矩阵，岭参数k设置为远小于1的非零值，且基线参数向量的初始值其中前四个值为通过轨道状态矢量计算得到的初始基线参数，并设置参考点平地相位初始值为φ₀＝0；

步骤2:对式L＝f(X)+Δ在X⁰处进行泰勒级数展开，计算观测方程系数矩阵A和向量l，并组成法方程，并计算加权残差平方和SWST₁：

${\mathrm{SWST}}_1=\underset{i=1}{\overset{N}{\Σ}}{({\mathrm{\ΔL}}_i\·\frac{P_{ii}}{P_{sum}})}^2---\left(8\right)$

其中，ΔL_i＝L_i-f_i(X⁰)，以及令SWST_old＝SWST₁；

步骤3:在法方程的系数阵对角线所有元素上分别加岭参数k，组成岭估计法方程式，并计算相应法方程系数阵的广义逆；

步骤4:根据岭估计参数平差方法计算参数改正数x的估计值对式L＝f(X)+Δ在处进行泰勒级数展开，计算新的观测方程系数矩阵A′和向量l′，并根据式计算加权残差平方和SWST_new；

步骤5:比较SWST_new和SWST_old的大小，按以下方法判断是否接受本次迭代所得参数向量改正数并输出加权残差平方和

如果SWST_new<SWST_old，则由式计算更新以及更新设计矩阵A＝A′和向量l＝l′，且调整岭估计参数，同时根据残差更新权矩阵P，以及设置SWST_old＝SWST_new；否则，保持X⁰，A，l和P不变，将调整岭参数；

步骤6:判断迭代终止条件：如果迭代次数为1，则定义令并直接跳转到步骤3进行第二次迭代；如果迭代次数大于1，则判断迭代终止条件是否满足，具体终止条件设置为

连续满足两次，或者迭代次数达到设定的最大次数；如果终止条件满足，则输出参数向量估计结果否则，令并跳转到步骤3继续迭代；

通过上述的迭代参数求解过程，基线参数向量的估计结果已精确获得，其包含了基线模型的参数估计值和即SAR干涉基线的精确估计结果，这四个基线模型参数能用于后续的高精度InSAR数据处理。