一种基于非线性优化策略的SBAS-DInSAR方法

摘要

本发明公开了一种基于非线性优化策略的SBAS-DInSAR方法，涉及雷达技术领域。该方法无需相位解缠即可实现对方程的解算，通过建立缠绕的差分干涉相位的非线性模型，构造与地表形变速率相关的优化目标函数，然后通过拟牛顿方法解算方程并提取地表形变信息，实现了一种新的SBAS-DInSAR地表形变测量方法。该方法避免了相位解缠算法引入的误差，同时提高了形变演算效率；在仅利用缠绕相位的情况下仍能获得高精度的地表形变结果，为地表形变测量提供了一种新的途径。

说明书

技术领域

本发明涉及雷达技术领域，具体地说，是指一种基于非线性优化策略的小基线集合成孔 径雷达差分干涉测量技术。

背景技术

小基线集合成孔径雷达差分干涉测量(SBAS-DInSAR)技术是一种将SAR图像按照小 基线原则组合成若干个子集，即：集合内基线距较小而集合间基线距较大，再利用在长时间 间隔内地表散射特性保持稳定的高相干点相位信息，建立高相干点的差分相位模型，通过模 型的解算来获取地表形变的测量技术。如何有效解算高相干点形变速率是保证高精度地表形 变测量的关键因素之一。

传统的SBAS方法利用解缠后的差分干涉相位，结合线性最小二乘法，并利用奇异值分 解(SVD)方法联合小基线集求得最小范数最小二乘解。然而，该方法需要相位解缠，虽已 有多种相位解缠算法，但各有其弊端，且在时间跨度较大的时序图像中，由于地表形变情况 复杂，往往导致解缠算法在某些区域解缠效果不理想，相位解缠的可靠性将直接影响此方法 的有效性。同时，解缠算法时间复杂度极高，使得该方法效率较低。

因此，针对SBAS-DInSAR地表形变测量技术，如何实现高精度、高效率的解算结果是 十分必要的。

发明内容

本发明的目的是针对已有技术的不足之处，提出一种基于非线性优化策略的 SBAS-DInSAR方法。该方法无需相位解缠即可实现对方程的解算，通过建立缠绕的差分干涉 相位的非线性模型，构造与地表形变速率相关的优化目标函数，然后通过拟牛顿方法解算方 程并提取地表形变信息，实现了一种新的SBAS-DInSAR地表形变测量方法。

本发明的技术方案如下：

基于非线性优化策略的SBAS-DInSAR方法，包括以下几个步骤：

步骤一：输入包含地表形变信息的长时间SAR图像序列，设共有K+1幅SAR图像，获 取时刻分别为t₀，t₁，…，t_K，其中K是正整数；对SAR图像进行预处理，包括配准、干涉、去 除参考面相位、去除地形相位，获得差分干涉相位图像序列；设所述差分干涉相位图像序列 包括M幅差分干涉相位图像，其中M是正整数；第j幅差分干涉相位图像中主图像和辅图像 的获取时刻分别为和j是小于或等于M的正整数，IE_j和IS_j是小于或等于K的非负 整数，且IE_j＞IS_j；

步骤二：对所述差分干涉相位图像序列，提取高相干点，设高相干点的数量为H，其中 H是正整数；第h个高相干点记为x_h，第j幅差分干涉相位图像中第h个高相干点的非缠绕 差分干涉相位记为其中h是小于或等于H的正整数；

步骤三：对第j幅差分干涉相位图像的第h个高相干点x_h有：

${\mathrm{\Delta \Phi }}_j^{\mathrm{uw}}\left(x_h\right)={\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IE}}_j})-{\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IS}}_j})$

其中，和分别为高相干点x_h在时刻和时刻由于地表形变引起的 相位；

步骤四：建立目标函数：

$f\left({\zeta }_h\right)=\frac{1}{M}{l}^T·\exp \left[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\right]$

式中，

其中，i是虚数单位；(D)_pq为矩阵D第p行第q列的元素，p为大于零且小于或等于M 的正整数，q为大于零且小于或等于K的正整数,t_q和t_q-1分别是第q幅和第q-1幅SAR图像 的获取时刻，IS_p和IE_p分别是第p幅差分干涉相位图像中主图像和辅图像的序号，第p幅差 分干涉相位图像由获取时刻为和的SAR图像干涉生成；λ为雷达波长；ζ_h为目标函数 f(ζ_h)的自变量；为由1构成的列向量；的第j个分量表示第j幅差分 干涉相位图像中第h个高相干点的缠绕差分干涉相位；

步骤五：通过拟牛顿法对所述目标函数进行第一次优化，提取高相干点的线性形变速率；

步骤六：以线性形变速率作为迭代初值，对所述目标函数进行第二次优化，计算总形变 速率。

本发明的优点在于：

本发明提出了一种新的SBAS-DInSAR地表形变速率解算方法，该方法无需相位解缠， 避免了相位解缠算法引入的误差，同时提高了形变演算效率；在仅利用缠绕相位的情况下仍 能获得高精度的地表形变结果，为地表形变测量提供了一种新的途径。

附图说明

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

图2是本发明实施例中仿真生成的SAR图像；

图3是本发明实施例中仿真生成的差分干涉相位图像；

图4是本发明实施例中高相干点提取位置示意图；

图5是本发明实施例中传统SBAS方法和改进后的SBAS方法所获得的形变速率估计值 与理论值的比较；

图6是本发明实施例中传统SBAS方法和改进后的SBAS方法计算误差比较。

具体实施方式

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

本发明是一种基于非线性优化策略的SBAS-DInSAR方法，方法流程如图1所示，包括 以下几个步骤：

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

输入K+1幅SAR图像，获取时刻分别为t₀，t₁，…，t_K，其中K是正整数；按照小基线原则 组合成若干个子集，也即根据相干性要求限制适当的空间基线距和时间基线距，将基线距满 足条件的图像配准后进行干涉处理，生成M幅干涉图像，其中M是正整数且 去除参考面相位和地形相位后，获得M幅差分干涉相位图像；第j幅 差分干涉相位图像中主图像和辅图像时刻分别为和j是小于或等于M的正整数，IE_j和IS_j是小于或等于K的非负整数，且IE_j＞IS_j；

步骤二：对差分干涉相位图像序列，提取高相干点；

根据K+1幅SAR图像，基于振幅信息提取高相干点；在M幅差分干涉相位图像中提取 每个高相干点的差分干涉相位设高相干点的数量为H，第h个高相干点记为x_h， 第j幅差分干涉相位图像中第h个高相干点的非缠绕差分干涉相位记为其中h是 小于或等于H的正整数；

步骤三：对第j幅差分干涉相位图像的第h个高相干点x_h有：

${\mathrm{\Delta \Phi }}_j^{\mathrm{uw}}\left(x_h\right)={\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IE}}_j})-{\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IS}}_j})---\left(1\right)$

其中，和分别为点x_h在时刻和时刻由于地表形变引起的相位；

步骤四：建立目标函数：

$f\left({\zeta }_h\right)=\frac{1}{M}\mathrm{\Sigma exp}\left[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\right]=\frac{1}{M}{l}^T·\exp \left[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\right]---\left(2\right)$

式中，

其中，i是虚数单位；Σ·表示对矩阵·的所有元素进行求和；(D)_pq为矩阵D第p行第q 列的元素，p为大于零且小于或等于M的正整数，q为大于零且小于或等于K的正整数,t_q和 t_q-1分别是第q幅和第q-1幅SAR图像的获取时刻，IS_p和IE_p分别是第p幅差分干涉相位图 像中主图像和辅图像的序号，第p幅差分干涉相位图像由获取时刻为和的SAR图像干 涉生成；λ为雷达波长；ζ_h为目标函数f(ζ_h)的自变量；为由1构成的列向量；的第j个分量表示第j幅差分干涉相位图像中第h个高相干点的缠绕差分干涉相位；

根据式(1)可以建立差分相位与地表形变速率的数学模型，定义平均相位形变速率 V(x_h,t_k)为：

$V(x_h,t_k)=\frac{{\Phi }^{\mathrm{uw}}(x_h,t_k)-{\Phi }^{\mathrm{uw}}(x_h,t_{k-1})}{t_k-t_{k-1}}---\left(6\right)$

结合式(1)和式(6)，可获得：

$\underset{k={\mathrm{IS}}_j+1}{\overset{{\mathrm{IE}}_j}{\Sigma }}(t_k-t_{k-1})V(x_h,t_k)={\mathrm{\Delta \Phi }}_j^{\mathrm{uw}}\left(x_h\right)---\left(7\right)$

根据地表形变速率v(x_h,t_k)与平均相位形变速率的关系：

$V(x_h,t_k)=\frac{4\pi }{\lambda }v(x_h,t_k)---\left(8\right)$

以及式(7)，可获得地表形变速率满足的关系式：

$\frac{4\pi }{\lambda }D·v_h={\mathrm{\Delta \Phi }}_h^{\mathrm{uw}}---\left(9\right)$

式中，D的定义如式(3)所示，

实际获得的差分干涉相位往往是缠绕的，即式(9)中的非缠绕相位无法直接获取， 因此式(9)无法使用线性最小二乘法直接求解；为了避免复杂的相位解缠，将式(9)改写 为：

$\frac{4\pi }{\lambda }D·v_h={\mathrm{\Delta \Phi }}_h^w+2\mathrm{\pi INT}---\left(12\right)$

其中，为由整数构成的列向量，表示缠绕差分干涉相位值；将式(12)变换 到复数域，进行如下变换：

$\exp (i\frac{4\pi }{\lambda }D·v_k)=\exp \left[i({\mathrm{\Delta \Phi }}_k^w+2\mathrm{\pi INT})\right]=\exp \left(i{\mathrm{\Delta \Phi }}_k^w\right)---\left(13\right)$

即：

$\exp \left[i(\frac{4\pi }{\lambda }D·v_k-{\mathrm{\Delta \Phi }}_k^w)\right]=l---\left(14\right)$

其中，为由1构成的列向量；

构建目标函数如下：

$f\left({\zeta }_h\right)=\frac{1}{M}\mathrm{\Sigma exp}\left[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\right]=\frac{1}{M}{l}^T·\exp \left[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\right]---\left(15\right)$

当目标函数|f(ζ_k)|取得最大值时，即可获得地表形变速率v_h，即：

$v_h=\arg \underset{{\zeta }_h}{\max }\left|f\left({\zeta }_h\right)\right|=\arg \underset{{\zeta }_h}{\max }\left|\frac{1}{M}\mathrm{\Sigma exp}\right[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\left]\right|---\left(16\right)$

估计形变速率过程是一个反复迭代的最优化求解过程，在不断迭代和优化的过程中逼近 所求真值，实现模型的解算；

步骤五：通过拟牛顿法对目标函数进行第一次优化，提取高相干点的线性形变速率；

令：

g(ζ_h)＝-|f(ζ_h)|²＝-f(ζ_h)f^*(ζ_h)  (17)

则式(16)等价于：

$v_h=\arg \underset{{\zeta }_h}{\min }g\left({\zeta }_h\right)=\arg \underset{{\zeta }_h}{\min }\{-{\left|\frac{1}{M}\mathrm{\Sigma exp}\right[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\left]\right|}^2\}---\left(18\right)$

拟牛顿法是一种求解非线性优化问题的方法，式(18)共有K个待优化参数；待优化参 数数量多、目标函数形式复杂等因素使迭代初值的选取对最终计算结果影响甚大，因此首先 应确定如何选取迭代初值；本发明采取的方法是首先认为v_h的各个分量相等，且均等于线性 形变速率，利用式(18)求解线性形变速率

在计算线性形变速率时，假定v_h的各个分量均等于线性形变速率则式(18)的最 优化问题等价于：

$v_h^{\mathrm{lin}}=\arg \underset{{\zeta }_h^{\mathrm{lin}}}{\min }g({\zeta }_h^{\mathrm{lin}}·l)=\arg \underset{{\zeta }_h^{\mathrm{lin}}}{\min }\{-{\left|\frac{1}{M}\mathrm{\Sigma exp}\right[i(\frac{4\pi }{\lambda }D·{\zeta }_h^{\mathrm{lin}}·l-{\mathrm{\Delta \Phi }}_h^w)\left]\right|}^2\}---\left(19\right)$

式中，为目标函数的自变量；

拟牛顿法迭代求解过程中，需使用目标函数的一阶导数信息；对目标函数求导：

$\frac{\partial g({\zeta }_h^{\mathrm{lin}}·l)}{\partial {\zeta }_h^{\mathrm{lin}}}=-\frac{\partial f({\zeta }_h^{\mathrm{lin}}·l)}{\partial {\zeta }_h^{\mathrm{lin}}}{f}^{*}({\zeta }_h^{\mathrm{lin}}·l)-f({\zeta }_h^{\mathrm{lin}}·l)\frac{{\partial f}^{*}({\zeta }_h^{\mathrm{lin}}·l)}{\partial {\zeta }_h^{\mathrm{lin}}}=-2\mathrm{Re}\left[f^{*}({\zeta }_h^{\mathrm{lin}}·l)\frac{\partial f({\zeta }_h^{\mathrm{lin}}·l)}{\partial {\zeta }_h^{\mathrm{lin}}}\right]---\left(20\right)$

其中，

其中，ο表示矩阵对应元素相乘；使用式(19)、式(20)和式(21)，结合拟牛顿法即可求解线 性形变速率

步骤六：以线性形变速率作为迭代初值，对目标函数进行第二次优化，计算总形变速率；

经过第一次优化获得线性形变速率后，以作为第二次优化的迭代初值，对v_h的 各个分量进行优化，计算总形变速率；

使用拟牛顿法求解式(18)，在拟牛顿法求解过程中需要使用目标函数g(ζ_h)的一阶导数：

$\frac{\partial g\left({\zeta }_h\right)}{\partial {\zeta }_h}=-2\mathrm{Re}\left[f^{*}\left({\zeta }_h\right)\frac{\partial f\left({\zeta }_h\right)}{\partial {\zeta }_h}\right]---\left(22\right)$

$\frac{\partial f\left({\zeta }_h\right)}{\partial {\zeta }_h}=\frac{i}{M}{\left\{\exp \right[i(\frac{4\pi }{\lambda }D·{\zeta }_h-{\mathrm{\Delta \Phi }}_h^w)\left]\right\}}^T\left(\frac{4\pi }{\lambda }D\right)---\left(23\right)$

使用式(18)、式(22)和式(23)，结合拟牛顿法即可求解地表形变速率v_h。

实施例：

本发明是一种基于非线性优化策略的SBAS-DInSAR技术，具体实施例为：

步骤一：以雷达天线重复轨道飞行模式对设定的场景仿真生成K+1＝17幅SAR图像。仿 真场景中心有一直径为700m的圆锥状山峰，初始高度为119.96m，其峰高以每年0.01m的平 均速率增长，除山峰外场景中其余点的高程值始终为零。将山峰高度达到120m的时刻记为 第0年，时间跨度为8年，以0.5年的采样间隔生成SAR图像序列。在场景中布设110个后 向散射系数较大的散射点作为高相干点，如图2所示。

按照小基线干涉图像对组合原则确定主辅图像，主图像不唯一。本发明是在时间基线不 大于3年、空间垂直基线不超过200米、相关系数不小于0.2的条件下，组成小基线集，将 每对主辅图像配准后做干涉处理，并去除参考面相位和地形相位，获得M＝26幅差分干涉图 像，如图3所示。图像干涉对组合如表1所示。

表1小基线干涉图像对组合参数

步骤二：在SAR图像中基于幅值法提取高相干点的位置，高相干点数量H＝110，对应于 场景中布置的110个后向散射系数较大的点。高相干点提取位置如图4所示。

步骤三：对第j幅差分干涉相位图上第h个高相干点x_h，计算其差分干涉相位：

${\mathrm{\Delta \Phi }}_j^{\mathrm{uw}}\left(x_h\right)={\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IE}}_j})-{\Phi }^{\mathrm{uw}}(x_h,t_{{\mathrm{IS}}_j}).$

步骤四：根据式(2)，对第h个高相干点，构建待优化的目标函数f(ζ_h)。

步骤五：将目标函数f(ζ_h)的最大化问题转化为目标函数g(ζ_h)的最小化问题，通过拟牛 顿法对目标函数进行第一次优化，获得高相干点处的线性形变速率

步骤六：步骤五求得线性形变速率后，以作为迭代初值，通过拟牛顿法对目标 函数g(ζ_h)进行第二次优化，计算总形变速率。

以第0至0.5年的形变速率为例，需要相位解缠的传统SBAS方法和本发明提出的基于 非线性优化策略的改进SBAS方法分别与理论值比较，如图5所示。可以看出，传统SBAS 方法和改进后的SBAS方法解算值与理论值都十分吻合。比较传统SBAS方法与改进的SBAS 方法的计算误差，如图6所示，可以看出，两种方法的计算精度相当，解算误差均在9×10^-4m/a 以下。

表2是传统SBAS方法和改进后的SBAS方法的算法耗时比较结果，可以看出，改进的 SBAS方法无需进行相位解缠，计算效率高于传统的SBAS方法。例如，使用枝切法对512×512 的DInSAR图像进行相位解缠，需花费约35s的时间(运行环境：CPU主频2.3GHz，内存4G， Matlab R2013b)，SVD分解求解线性方程组需约25ms的时间，而且当图像噪声水平很大或失 相关现象较严重时，相位解缠结果可能非常不可靠；而改进的SBAS方法在同样的计算条件 下求解最优化问题仅需约2.3s的时间，高于传统SBAS方法SVD分解求解线性方程组耗费的 时间，但远低于传统SBAS方法相位解缠付出的时间代价，算法效率提升14倍以上。

表2传统SBAS方法和改进SBAS方法算法耗时比较

本发明提出一种基于非线性优化策略的SBAS-DInSAR方法，在无需相位解缠的情况下， 直接建立差分干涉相位与形变速率的非线性模型，构造与地表形变信息相关的优化目标函数， 通过拟牛顿法进行一次优化，提取高相干点的线性形变速率，并以该线性形变速率作为二次 优化的迭代初值，再次使用拟牛顿法提取高相干点的总形变速率。在与传统方法计算精度相 当的情况下，大大提高了算法效率，实现了一种新的SBAS-DInSAR技术。通过实施例分析， 进一步详述了本发明方法的实施过程，验证了本发明方法的正确性以及高效性。