CN104091064A - 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法 - Google Patents

Abstract

本发明公开了基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，包括步骤一：获得差分干涉相位图像序列；步骤二：对差分干涉相位图像序列，提取永久散射体点，构建Delaunay三角网络；步骤三：对第k幅差分干涉相位图像，计算每一对相邻PS点的二次差分相位；步骤四：建立待优化的目标函数；步骤五：对目标函数在二维解空间进行搜索；步骤六：使用Levenberg-Marquardt算法对目标函数进行局部优化；步骤七：解算地表形变量和高程误差。本发明不涉及信号采样问题，因此不受各干涉图像对时间和空间基线不均匀性影响，在缺乏先验知识的情况下仍能获得高精度的结果，为地表形变测量提供了新的思路和途径。

Description

基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法

技术领域

本发明涉及雷达技术领域，具体地说，是指一种基于优化解空间搜索法的永久散射体合成孔径雷达差分干涉测量技术的地表形变测量方法。

背景技术

永久散射体合成孔径雷达差分干涉测量(PS-DInSAR)是一种利用地表散射特性稳定的永久散射体(PS)相位信息，建立合理的相邻PS点差分相位模型，通过模型求解获得地表形变、高程误差等的测量技术。然而，大气延迟、非线性形变、相位缠绕等因素使模型求解极其困难，往往需要合理估计模型中相邻PS点形变速率之差Δv、高程误差之差等参数。参数估算过程的准确性，是保证高精度地表形变测量的关键因素之一。

二维周期图法是一种常用的PS-DInSAR参数估计方法，将差分相位看作是T-B空间(时间基线-空间基线空间)到空间的二维傅里叶变换值，利用频谱估计的方法估计Δv和在时间基线和空间基线分布不均匀、SAR图像数量较少的情况下，二维周期图法往往不能获得合理的估计结果，估计误差较大。

解空间搜索法是另一种估计Δv和的方法，首先在一定先验知识的指导下确定Δv和可能的取值范围作为解空间，然后以一定搜索步长计算解空间内所有点的相位相干系数，选取相位相干系数最大的点作为Δv和的估计值。在缺乏先验知识的情况下，为了提高参数估计精度，往往需要使用较大范围的解空间和较小的搜索步长，限制了算法效率的提高；此外，解空间内相位相干系数最大的点往往未必是Δv和的真实值，存在较大的估计误差。

综上，尚未有一种针对PS-DInSAR地表形变测量的高精度、高效率的参数估计方法。

发明内容

本发明的目的是针对已有技术的不足之处，提出一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法。在缺乏先验知识的情况下，通过以较大的搜索步长对较大搜索空间进行搜索，初步确定地表形变速率之差和高程误差之差粗值，在此基础上通过对相位相干系数进一步优化，在兼顾算法效率的同时，获得高精度的地表形变测量值，实现了一种新的PS-DInSAR地表形变测量参数估计方法。

本发明的技术方案如下：

基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，包括以下步骤：

步骤一：输入包含地表形变信息的长时间SAR图像序列，并对SAR图像进行预处理，包括配准、干涉、去除参考面相位、去除地形相位，获得差分干涉相位图像序列，设差分干涉相位图像序列包括K张图像；

步骤二：对差分干涉相位图像序列，提取永久散射体点，即PS点，设PS点的数量为N，第n个PS点记为x_n，第k幅差分干涉相位图像中第n个PS点的差分干涉相位记为其中n是小于或等于N的正整数，k是小于或等于K的正整数，据提取PS点的坐标构建Delaunay三角网络；

步骤三：从空间尺度上，对第k幅差分干涉相位图像，计算每一对相邻PS点的二次差分相位：

${\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})={\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{r_p}\right)-{\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{s_p}\right)$

其中，p为小于或等于Delaunay三角网络总边数的正整数，r_p、s_p均是小于或等于N的正整数；

步骤四：建立目标函数：

$\begin{array}{c}{\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]\\ =\left|\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\exp \right\{j[{\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})-B_k^{\⊥}\·\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})-C_v\·T_k\·\mathrm{\Δv}(x_{r_p},x_{s_p})]\left\}\right|\end{array}$

式中，

$\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})=\mathrm{\ΔQ}\left(x_{r_p}\right)-\mathrm{\ΔQ}\left(x_{s_p}\right)$

$\mathrm{\ΔQ}\left(x_{r_p}\right)=C_q\left(x_{r_p}\right)\·\mathrm{\Δq}\left(x_{r_p}\right),\mathrm{\ΔQ}\left(x_{s_p}\right)=C_q\left(x_{s_p}\right)\·\mathrm{\Δq}\left(x_{s_p}\right)$

$C_q\left(x_{r_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{r_p}\sin {\mathrm{\α}}_{r_p}},C_q\left(x_{s_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{s_p}\sin {\mathrm{\α}}_{s_p}}$

$\mathrm{\Δv}(x_{r_p},x_{s_p})=v\left(x_{r_p}\right)-v\left(x_{s_p}\right)$

$C_v=\frac{4\mathrm{\π}}{\mathrm{\λ}}$

其中，j是虚数单位；为第k幅差分干涉相位图像的垂直空间基线；和分别为PS点和的高程误差值；T_k为第k幅辅图像与主图像的成像时间间隔；和分别为PS点和沿雷达视线方向的线性形变速率值；和分别为PS点和与获取主图像的雷达对应的斜距；和分别为PS点和与获取主图像的雷达对应的当地入射角；λ为雷达波长；γ_p为相位相干系数；

步骤五：对目标函数γ_p在二维空间中以一定的步长搜索，使得γ_p＞T_γ且最小，记作和其中T_γ是给定的阈值；

步骤六：以和作为迭代初值，使用Levenberg-Marquardt算法对目标函数γ_p进行局部优化，获得使γ_p取得极大值的和

步骤七：根据和解算地表形变量和高程误差。

本发明的优点在于：

(1)本发明提出了一种新的PS-DInSAR地表形变测量参数估计方法，该方法不涉及信号采样问题，因此不受各干涉图像对时间和空间基线不均匀性影响，在缺乏先验知识的情况下仍能获得高精度的结果，为地表形变测量提供了新的思路和途径。

(2)本发明采取两步优化方式，满足参数估计精度要求的同时提高了运算效率。

附图说明

图1是本发明的方法流程图；

图2是本发明实施例中仿真生成的差分干涉图像序列；

图3是本发明实施例中构建的三角网络示意图；

图4是本发明实施例中搜索空间内相位相干系数的取值；

图5是本发明实施例中未使用二次优化策略获得的形变速率估计值和真实形变速率的比较；

图6是本发明实施例中使用二次优化策略获得的形变速率估计值和真实形变速率的比较；

图7是本发明实施例中第5年形变量的计算结果。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，流程如图1所示，包括以下几个步骤：

步骤一：输入包含地表形变信息的长时间SAR图像序列，并对SAR图像进行预处理，包括配准、干涉、去除参考面相位、去除地形相位，获得差分干涉相位图像序列。

输入K+1幅SAR图像，选取一幅SAR图像作为主图像，其余K幅为辅图像；将辅图像与主图像配准，获得配准后的辅图像，并将每一幅辅图像与主图像进行干涉运算，组成K幅干涉相位图像，去除参考面相位和地形相位后获得K幅差分干涉相位图像。

步骤二：对差分干涉相位图像序列，提取永久散射体点，即PS点。

1、根据K幅差分干涉相位图像，基于幅值法等选取PS点，提取每个PS点的差分干涉相位设PS点的数量为N，第n个PS点记为x_n，第k幅差分干涉相位图像中第n个PS点的差分干涉相位记为其中n是小于或等于N的正整数，k是小于或等于K的正整数；

2、据提取PS点的坐标构建Delaunay三角网络，Delaunay三角网络所有边的集合记为S，设集合S中共有L个元素；S中的第p元素记为(r_p,s_p)，表示Delaunay三角网络第p条边的端点为PS点和其中r_p和s_p均为小于或等于N的正整数，p为小于或等于L的正整数。

步骤三：从空间尺度上，对第k幅差分干涉相位图像，计算Delaunay三角网络第p条边PS点(即相邻PS点)的二次差分相位：

${\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})={\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{r_p}\right)-{\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{s_p}\right)---\left(1\right)$

其中，p为小于或等于Delaunay三角网络总边数的正整数，r_p、s_p均是小于或等于N的正整数；

步骤四：建立目标函数：

$\begin{array}{c}{\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]\\ =\left|\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\exp \right\{j[{\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})-B_k^{\⊥}\·\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})-C_v\·T_k\·\mathrm{\Δv}(x_{r_p},x_{s_p})]\left\}\right|\end{array}---\left(2\right)$

式中，

$\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})=\mathrm{\ΔQ}\left(x_{r_p}\right)-\mathrm{\ΔQ}\left(x_{s_p}\right)---\left(3\right)$

$\mathrm{\ΔQ}\left(x_{r_p}\right)=C_q\left(x_{r_p}\right)\·\mathrm{\Δq}\left(x_{r_p}\right),\mathrm{\ΔQ}\left(x_{s_p}\right)=C_q\left(x_{s_p}\right)\·\mathrm{\Δq}\left(x_{s_p}\right)---\left(4\right)$

$C_q\left(x_{r_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{r_p}\sin {\mathrm{\α}}_{r_p}},C_q\left(x_{s_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{s_p}\sin {\mathrm{\α}}_{s_p}}---\left(5\right)$

$\mathrm{\Δv}(x_{r_p},x_{s_p})=v\left(x_{r_p}\right)-v\left(x_{s_p}\right)---\left(6\right)$

$C_v=\frac{4\mathrm{\π}}{\mathrm{\λ}}---\left(7\right)$

其中，j是虚数单位；为第k幅差分干涉相位图像的垂直空间基线；和分别为PS点和的高程误差值；T_k为第k幅辅图像与主图像的成像时间间隔；和分别为PS点和沿雷达视线方向的线性形变速率值；和分别为PS点和与获取主图像的雷达对应的斜距；和分别为PS点和与获取主图像的雷达对应的当地入射角；λ为雷达波长；γ_p为相位相干系数。

相邻PS点的二次差分相位由高程误差项、线性形变项和非线性因素等构成，即：

${\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})=B_k^{\⊥}\·\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})+C_v\·T_k\·\mathrm{\Δv}(x_{r_p},x_{s_p})+{\mathrm{\Δ\Φ}}_k^{\mathrm{PS},\mathrm{res}}(x_{r_p},x_{s_p})---\left(8\right)$

式中，和为高程误差和线性形变对二次差分相位的贡献项；为非线性形变、大气延迟、失相关噪声等非线性因素对二次差分相位的贡献项。非线性形变、大气延迟和失相关噪声的不确定性使得无法直接测量，相位缠绕使得式(8)无法使用线性最小二乘法直接求解。

研究表明，地表形变在空间上具有高度的自相关性，场景中相邻PS点的非线性形变之差是一个小量；同样，大气延迟也具有高度的空间自相关性，相邻PS点的大气延迟相近，相邻PS点的大气相位差很小；而PS点由于其相位和幅度变化稳定的特点，几乎不受时间和空间失相关的影响，因此相邻PS点的失相关噪声之差也是小量；因此一般是一个小量，且基于此，可建立目标函数：

${\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]=\left|\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\exp \right\{{\mathrm{j\Δ}}_k^{\mathrm{PS},\mathrm{res}}(x_{r_p},x_{s_p})\left\}\right|---\left(9\right)$

通过求解最优化问题：

$\{\▿\mathrm{\Δ}\tilde{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\tilde{v}(x_{r_p},x_{s_p})\}=\arg \max {\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]---\left(10\right)$

获得和根据式(8)可知，式(9)与式(2)是等价的。步骤五和步骤六将详述式(10)最优化问题的求解方法。

步骤五：对目标函数γ_p在二维空间中以一定的步长搜索，使得γ_p＞T_γ且最小，记作和其中T_γ是给定的阈值；

1、确定的搜索范围[Q_p,min,Q_p,max]，搜索步长Q_p,step；以及的搜索范围[v_p,min,v_p,max]，搜索步长v_p,step；

2、在所述的搜索空间内，选取使得γ_p＞T_γ且最小的和作为和的粗略估计值。由于PS点的散射特性稳定，相位相干系数γ_p一般较大，因此可选取一定阈值T_γ，要求相位相干系数γ_p＞T_γ，初步确定和的可能取值；若所监测的场景尺寸较小且形变缓慢，可认为其线性形变速率之差是一小量，因此在以上获得的可能取值内，选取使得最小的和可获得和的粗略估计值，记作记作和

步骤六：以和作为迭代初值，使用Levenberg-Marquardt算法对目标函数γ_p进行局部优化，获得使γ_p取得极大值的和

步骤五中获得的粗略估计值和的精度取决于搜索步长Q_p,step和v_p,step，且估计值和不一定是γ_p的真实极大值点。以和作为迭代初值，使用Levenberg-Marquardt算法求解目标函数γ_p的极大值点，获得精确的估计值和具体方法如下：

1、建立目标函数：

$f_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]=-{\left|{\mathrm{\γ}}_p\right[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})\left]\right|}^2---\left(11\right)$

将目标函数γ_p的最大化问题转换为目标函数f_p的最小化问题。

2、基于Levenberg-Marquardt算法，调用Matlab的最优化函数fminunc，以和作为迭代初始值，求解最优化问题获得精确的估计值和具体的调用方式如下：

$[\▿\mathrm{\Δ}\tilde{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\tilde{v}(x_{r_p},x_{s_p})]=f\min \mathrm{unc}(f_p,[\▿\mathrm{\Δ}\hat{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\hat{v}(x_{r_p},x_{s_p})],\mathrm{options})---\left(12\right)$

式中，options为函数fminunc的迭代参数。

步骤七：根据和解算地表形变量和高程误差。

1、剔除的相邻PS对，其中T_γ′为选定的阈值。若相位相干系数γ_p较小，则认为对应的相邻PS对的二次差分相位信息是不可靠的，故将其剔除；

2、选取某个PS点作为参考点，给出其形变速率和高程误差信息，使用最小二乘平差法计算每个PS点的线性形变速率和高程误差；

3、根据式(8)计算，计算相邻PS点的使用最小二乘平差法和滤波法从中提取每个PS点的非线性形变量；

4、结合每个PS点的线性形变速率、高程误差和非线性形变量，使用Kriging插值法计算SAR图像中各个点的形变量和高程误差。

实施例：

本发明是一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，具体实施例为：

步骤一：仿真生成K+1＝17幅SAR图像。在SAR图像的场景中心，有一直径为700m的圆锥形山峰，其高度以每年0.01m的速率增长，初始高度为119.96m，除山峰外场景中其余点的高程始终为零。以0.5年的采样间隔生成SAR图像序列，选取第4.5年的SAR图像作为主图像，此时山峰高度为120m，其余图像作为辅图像。将每张辅图像与主图像配准，并做干涉运算，去除参考面相位和地形相位，获得K幅差分干涉相位图像，如图2所示。在SAR图像对应的场景中布置110个后向散射系数较大的散射点作为永久散射体，即PS点。

步骤二：在SAR图像中基于幅值法提取PS点的位置，PS点数量N＝110，对应于场景中布置的110个后向散射系数较大的点。调用Matlab中的delaunay函数构建Delaunay三角网络，如图3所示，该网格共含有L＝289条边。

步骤三：根据式(1)计算相邻PS点(即Delaunay三角网络每条边的两个端点)的二次差分相位。

步骤四：根据式(2)，对第p对相邻PS点，构建待优化的目标函数

${\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})].$

步骤五和步骤六：选取的搜索范围为 $[Q_{p,\min },Q_{p,\max }]=[{-10}^{-1},{10}^{-1}]$ 米^-1，的搜索范围为 $[v_{p,\min },v_{p,\max }]=[-{10}^{-1},{10}^{-1}]$ 米/年。

图4为在上述搜索空间中，搜索步长Q_p,step＝10^-3米^-1，v_p,step＝10^-1米/年时，第46对相邻PS点的取值。若使用解空间搜索法对和进行估计，则应选取使相位相干系数γ_p的最大值点，此时和的估计值分别为0米^-1和-0.093米/年。实际上，第46对相邻PS点的和真实值分别为0米^-1和0米/年，解空间搜索法的估计值偏差较大。

本发明提出的优化解空间搜索法首先选定一定的阈值T_γ(此处选取)，初选出图4中γ_p＞T_γ的极大值点，并在这些极大值点中选取最小的极大值点作为初始估计值和此时计算获得的粗略估计值分别为0米^-1和0.001米/年，估计精度高于解空间搜索法；然后以此粗略估计值作为迭代初值，使用Levenberg-Marquardt算法二次寻找γ_p极大值的精确位置米^-1，米/年。图5和图6分别是实施例中使用Levenberg-Marquardt算法优化前后线性形变速率计算值和理论值的比较，平均误差分别为2.11×10^-4米/年和2.59×10^-6米/年。可以看出，使用Levenberg-Marquardt算法进行二次优化后，形变速率的计算值和理论值的吻合程度远高于优化前。

表1是未使用Levenberg-Marquardt算法二次优化时的算法耗时和形变速率估计误差(运行环境：2.30GHz4核8线程CPU、4G RAM、Matlab2013b)，可以发现，随着搜索步长的减小，估计误差越小，但耗时越长。而本发明提出的算法在搜索步长Q_p,step＝10^-3米^-1，v_p,step＝10^-1米/年时平均误差已可达2.59×10^-6米/年，算法耗时约11.8秒。可见，使用Levenberg-Marquardt算法二次优化方法可有效提高估计精度，且算法效率高。

表1 未使用二次优化算法时的算法耗时和形变速率估计误差结果

步骤七：根据和解算地表形变量和高程误差。

剔除 ${\mathrm{\γ}}_p[\▿\mathrm{\Δ}\tilde{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\tilde{v}(x_{r_p},x_{s_p})]\≤T_{\mathrm{\γ}}^{\′}=0.70$ 的相邻PS对，然后使用最小二乘平差、滤波等方法计算线性形变速率、高程误差、非线性形变量等，最后使用Kriging插值方法获得SAR图像中每一点的形变量。图7所示的是第5年形变量的计算结果。

本发明提出一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法。在缺乏先验知识的情况下，通过以较大的搜索步长对较大搜索空间进行搜索，初步确定地表形变速率之差和高程误差之差粗值，在此基础上通过对相位相干系数进一步优化，在兼顾算法效率的同时，获得高精度的地表形变测量值，实现了一种新的PS-DInSAR地表形变测量参数估计方法。通过实施例分析，详述了本发明方法的实施过程，验证了本发明方法的正确性以及高精度、高效率等优点。

Claims (5)

1.一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，包括以下几个步骤：

步骤一：输入包含地表形变信息的长时间SAR图像序列，并对SAR图像进行预处理，包括配准、干涉、去除参考面相位、去除地形相位，获得差分干涉相位图像序列，得到K幅差分干涉相位图像；

步骤二：对差分干涉相位图像序列，提取永久散射体点，即PS点，设PS点的数量为N，第n个PS点记为x_n，第k幅差分干涉相位图像中第n个PS点的差分干涉相位记为其中n是小于或等于N的正整数，k是小于或等于K的正整数，据提取PS点的坐标构建Delaunay三角网络；

步骤三：从空间尺度上，对第k幅差分干涉相位图像，计算Delaunay三角网络第p条边PS点的二次差分相位：

${\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})={\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{r_p}\right)-{\mathrm{\Φ}}_k^{\mathrm{PS}}\left(x_{s_p}\right)---\left(1\right)$

其中，p为小于或等于Delaunay三角网络总边数的正整数，r_p、s_p均是小于或等于N的正整数；

步骤四：建立目标函数：

$\begin{array}{c}{\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]\\ =\left|\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\exp \right\{j[{\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})-B_k^{\⊥}\·\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})-C_v\·T_k\·\mathrm{\Δv}(x_{r_p},x_{s_p})]\left\}\right|\end{array}---\left(2\right)$

式中，

$\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})=\mathrm{\ΔQ}\left(x_{r_p}\right)-\mathrm{\ΔQ}\left(x_{s_p}\right)---\left(3\right)$

$\mathrm{\ΔQ}\left(x_{r_p}\right)=C_q\left(x_{r_p}\right)\·\mathrm{\Δq}\left(x_{r_p}\right),\mathrm{\ΔQ}\left(x_{s_p}\right)=C_q\left(x_{s_p}\right)\·\mathrm{\Δq}\left(x_{s_p}\right)---\left(4\right)$

$C_q\left(x_{r_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{r_p}\sin {\mathrm{\α}}_{r_p}},C_q\left(x_{s_p}\right)=\frac{4\mathrm{\π}}{{\mathrm{\λR}}_{s_p}\sin {\mathrm{\α}}_{s_p}}---\left(5\right)$

$\mathrm{\Δv}(x_{r_p},x_{s_p})=v\left(x_{r_p}\right)-v\left(x_{s_p}\right)---\left(6\right)$

$C_v=\frac{4\mathrm{\π}}{\mathrm{\λ}}---\left(7\right)$

其中，j是虚数单位；为第k幅差分干涉相位图像的垂直空间基线；和分别为PS点和的高程误差值；T_k为第k幅辅图像与主图像的成像时间间隔；和分别为PS点和沿雷达视线方向的线性形变速率值；和分别为PS点和与获取主图像的雷达对应的斜距；和分别为PS点和与获取主图像的雷达对应的当地入射角；λ为雷达波长；γ_p为相位相干系数；

相邻PS点的二次差分相位由高程误差项、线性形变项和非线性因素等构成，即：

${\mathrm{\Δ\Φ}}_k^{\mathrm{PS}}(x_{r_p},x_{s_p})=B_k^{\⊥}\·\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p})+C_v\·T_k\·\mathrm{\Δv}(x_{r_p},x_{s_p})+{\mathrm{\Δ\Φ}}_k^{\mathrm{PS},\mathrm{res}}(x_{r_p},x_{s_p})---\left(8\right)$

式中，和为高程误差和线性形变对二次差分相位的贡献项；为非线性形变、大气延迟、失相关噪声等非线性因素对二次差分相位的贡献项；非线性形变、大气延迟和失相关噪声的不确定性使得无法直接测量，相位缠绕使得式(8)无法使用线性最小二乘法直接求解；

因为建立目标函数：

${\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]=\left|\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\exp \right\{{\mathrm{j\Δ}}_k^{\mathrm{PS},\mathrm{res}}(x_{r_p},x_{s_p})\left\}\right|---\left(9\right)$

通过求解最优化问题：

$\{\▿\mathrm{\Δ}\tilde{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\tilde{v}(x_{r_p},x_{s_p})\}=\arg \max {\mathrm{\γ}}_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]---\left(10\right)$

获得和根据式(8)可知，式(9)与式(2)是等价的；

步骤五：对目标函数γ_p在二维空间中以搜索步长进行搜索，使得γ_p＞T_γ且最小，记作和其中T_γ是给定的阈值；

步骤六：以和作为迭代初值，使用Levenberg-Marquardt算法对目标函数γ_p进行局部优化，获得使γ_p取得极大值的和

步骤七：根据和解算地表形变量和高程误差。

2.根据权利要求1所述的一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，所述的步骤二具体包括：

(1)、根据K幅差分干涉相位图像，基于幅值法等选取PS点，提取每个PS点的差分干涉相位设PS点的数量为N，第n个PS点记为x_n，第k幅差分干涉相位图像中第n个PS点的差分干涉相位记为其中n是小于或等于N的正整数，k是小于或等于K的正整数；

(2)、据提取PS点的坐标构建Delaunay三角网络，Delaunay三角网络所有边的集合记为S，设集合S中共有L个元素；S中的第p元素记为(r_p,s_p)，表示Delaunay三角网络第p条边的端点为PS点和其中r_p和s_p均为小于或等于N的正整数，p为小于或等于L的正整数。

3.根据权利要求1所述的一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，所述的步骤五具体包括：

(1)、确定的搜索范围[Q_p,min,Q_p,max]，搜索步长Q_p,step；以及的搜索范围[v_p,min,v_p,max]，搜索步长v_p,step；

(2)、在所述的搜索空间内，选取使得γ_p＞T_γ且最小的和作为和的粗略估计值，具体包括：选取阈值T_γ，要求相位相干系数γ_p＞T_γ，初步确定和取值；选取使得最小的和获得和的粗略估计值，记作和

4.根据权利要求1所述的一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，所述的步骤六具体包括：

(1)、建立目标函数：

$f_p[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})]=-{\left|{\mathrm{\γ}}_p\right[\▿\mathrm{\ΔQ}(x_{r_p},x_{s_p}),\mathrm{\Δv}(x_{r_p},x_{s_p})\left]\right|}^2---\left(11\right)$

将目标函数γ_p的最大化问题转换为目标函数f_p的最小化问题；

(2)、基于Levenberg-Marquardt算法，调用Matlab的最优化函数fminunc，以和作为迭代初始值，求解最优化问题获得精确的估计值和

具体的调用方式如下：

$[\▿\mathrm{\Δ}\tilde{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\tilde{v}(x_{r_p},x_{s_p})]=f\min \mathrm{unc}(f_p,[\▿\mathrm{\Δ}\hat{Q}(x_{r_p},x_{s_p}),\mathrm{\Δ}\hat{v}(x_{r_p},x_{s_p})],\mathrm{options})---\left(12\right)$

式中，options为函数fminunc的迭代参数。

5.根据权利要求1所述的一种基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法，所述的步骤七具体包括：

(1)、剔除的相邻PS对，其中T_γ′为选定的阈值；若相位相干系数γ_p较小，则认为对应的相邻PS对的二次差分相位信息是不可靠的，故将其剔除；

(2)、选取某个PS点作为参考点，给出其形变速率和高程误差信息，使用最小二乘平差法计算每个PS点的线性形变速率和高程误差；

(3)、根据式(8)计算，计算相邻PS点的使用最小二乘平差法和滤波法从中提取每个PS点的非线性形变量；

(4)、结合每个PS点的线性形变速率、高程误差和非线性形变量，使用Kriging插值法计算SAR图像中各个点的形变量和高程误差。