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

摘要

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

主权项

一种基于优化解空间搜索法的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{\Delta \Phi }}_k^{PS}(x_{r_p},x_{s_p})={\Phi }_k^{PS}\left(x_{r_p}\right)-{\Phi }_k^{PS}\left(x_{s_p}\right)---\left(1\right)$其中，p为小于或等于Delaunay三角网络总边数的正整数，r_p、s_p均是小于或等于N的正整数；步骤四：建立目标函数：$\begin{array}{c}{\gamma }_p[▿\Delta Q(x_{r_p},x_{s_p}),\Delta v(x_{r_p},x_{s_p})]\\ =\left|\frac{1}{K}\underset{k=1}{\overset{K}{\Sigma }}\exp \right\{j[{\mathrm{\Delta \Phi }}_k^{PS}(x_{r_p},x_{s_p})-B_k^{\perp }·▿\Delta Q(x_{r_p},x_{s_p})-C_v·T_k·v(x_{r_p},x_{s_p})]\left\}\right|\end{array}---\left(2\right)$式中，$▿\Delta Q(x_{r_p},x_{s_p})=\Delta Q\left(x_{r_p}\right)-\Delta Q\left(x_{s_p}\right)---\left(3\right)$$\Delta Q\left(x_{r_p}\right)=C_q\left(x_{r_p}\right)·\Delta q\left(x_{r_p}\right),\Delta Q\left(x_{s_p}\right)=C_q\left(x_{s_p}\right)·\Delta q\left(x_{s_p}\right)---\left(4\right)$$C_q\left(x_{r_p}\right)=\frac{4\pi }{{\mathrm{\lambda R}}_{r_p}{\mathrm{sin\alpha }}_{r_p}},C_q\left(x_{s_p}\right)=\frac{4\pi }{{\mathrm{\lambda R}}_{s_p}{\mathrm{sin\alpha }}_{s_p}}---\left(5\right)$$\Delta 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\pi }{\lambda }---\left(7\right)$其中，j是虚数单位；为第k幅差分干涉相位图像的垂直空间基线；和分别为PS点和的高程误差值；T_k为第k幅辅图像与主图像的成像时间间隔；和分别为PS点和沿雷达视线方向的线性形变速率值；和分别为PS点和与获取主图像的雷达对应的斜距；和分别为PS点和与获取主图像的雷达对应的当地入射角；λ为雷达波长；γ_p为相位相干系数；相邻PS点的二次差分相位由高程误差项、线性形变项和非线性因素构成，即：${\mathrm{\Delta \Phi }}_k^{PS}(x_{r_p},x_{s_p})=B_k^{\perp }·▿\Delta Q(x_{r_p},x_{s_p})+C_v·T_k·\Delta v(x_{r_p},x_{s_p})+{\mathrm{\Delta \Phi }}_k^{PS,res}(x_{r_p},x_{s_p})---\left(8\right)$式中，和为高程误差和线性形变对二次差分相位的贡献项；为非线性形变、大气延迟和失相关噪声的非线性因素对二次差分相位的贡献项；非线性形变、大气延迟和失相关噪声的不确定性使得无法直接测量，相位缠绕使得式(8)无法使用线性最小二乘法直接求解；因为建立目标函数：${\gamma }_p[▿\Delta Q(x_{r_p},x_{s_p}),\Delta v(x_{r_p},x_{s_p})]=\left|\frac{1}{K}\underset{k=1}{\overset{K}{\Sigma }}\exp \right\{{\mathrm{j\Delta \Phi }}_k^{PS,res}(x_{r_p},x_{s_p})\left\}\right|---\left(9\right)$步骤五：对目标函数γ_p在二维空间中以搜索步长进行搜索，使得γ_p＞T_γ且最小，记作和其中T_γ是给定的阈值；步骤六：以和作为迭代初值，使用Levenberg‑Marquardt算法对目标函数γ_p进行局部优化，获得使γ_p取得极大值的和步骤七：根据和解算地表形变量和高程误差。