CN114019507A - 一种基于l1/2范数正则化的sar干涉图生成方法 - Google Patents

Abstract

本发明公开了一种基于L_1/2范数正则化的SAR干涉图生成方法，具体为：步骤1：采用SAR获取同一目标场景的SAR主、副图像，通过对副图像复矩阵取相位值获取SAR副图像的相位；步骤2：基于SAR主、副图像的干涉相位，以及通过解决L_1/2范数正则化问题建立SAR干涉图的生成模型；步骤3：采用阈值迭代算法求解L_1/2范数正则化问题，得到与恢复后的干涉图像有关的小波系数矩阵；步骤4：对步骤3求得的小波系数矩阵进行离散小波逆变换，取变换后的角度值作为恢复后的干涉图像。采用本发明重建的干涉图像具有更好的性能和更小的相位噪声；采用本发明获得的SAR干涉图质量更高。

Description

一种基于L1/2范数正则化的SAR干涉图生成方法

技术领域

本发明属于干涉合成孔径雷达成像和稀疏信号处理技术领域。

背景技术

SAR是一种高分辨率成像雷达，相比于光学成像方法，可以全天时全天候获取大场景、宽测绘带的图像，是目前微波成像技术的一种重要的研究方向。传统的SAR成像技术局限于二维(Two-dimensional，简称2D)成像，而干涉SAR(SAR interferometry，简称InSAR)成像技术突破了这一局限，它利用两幅SAR复图像的相位差来恢复目标场景的地形图，从而获得观测目标的高度信息和数字高程模型(Digital elevation model，简称DEM)，因此被广泛应用于地形测量、地形勘探等领域。高质量的DEM对干涉图的处理精度有着较高的要求。目前，典型的相位噪声滤波技术，例如Lee滤波器、Goldstein 滤波器等，都将导致图像空间分辨率的降低。

发明内容

发明目的：为了解决上述现有技术存在的问题，本发明提供了一种基于L_1/2范数正则化的SAR干涉图生成方法。

技术方案：本发明提供了一种基于L_1/2范数正则化的SAR干涉图生成方法，具体包括如下步骤：

步骤1：采用SAR获取同一目标场景的SAR主图像和SAR副图像，通过对副图像复矩阵取相位值获取SAR副图像的相位；

步骤2：基于SAR主、副图像的干涉相位，以及通过解决L_1/2范数正则化问题建立SAR干涉图的生成模型；

步骤3：采用阈值迭代算法求解L_1/2范数正则化问题，得到与恢复后的SAR干涉图有关的小波系数矩阵

步骤4：对步骤3求得的

进行离散小波逆变换，取变换后的角度值作为恢复后的干涉图像。

进一步的，所述步骤1中SAR副图像的相位y为：

θ_m＝exp{j(φ_m-φ_flat)}

u＝exp(-jφ_topo)

n＝θ_mu{exp(-jφ_noise)-1}

其中，φ_m为SAR主图像y_m的相位，φ_flat为平地相位，φ_topo为地形相位；φ_noise＝φ_s2-φ_s1，φ_s1和φ_s2分别表示y_m和y_s的散射相位，y_s表示SAR副图像，exp{.}为e的指数次方，j 表示虚部。

进一步的，所述步骤2中SAR干涉图的生成模犁为：

Y＝Θ_mοR(X)+N

其中，

表示将矢量形式的y重新构写成矩阵形式的Y，Y为幅度为1的二维SAR副图像矩阵；

表示将矢量形式的θ_m重新构写成矩阵形式的Θ_m，Θ_m为含有SAR 主图像相位和平地相位的矩阵；R(.)为离散小波逆变换；。表示哈达玛积；X＝W(U)，

表示将矢量形式的u重新构写成矩阵形式的U，U表示含有重建后的SAR主、副图像的干涉相位的矩阵；

表示将矢量形式的n重新构写成矩阵形式的N，N表示噪声矩阵，β表示正则化参数，

表示斐波那契范数，argmin{.}表示取最小值，W(.) 为离散小波变换。

进一步的，所述步骤3具体为：

步骤3.1：设置小波系数的初始值X⁽⁰⁾＝0；所述小波系数为SAR主、副图像的干涉相位经过DWT变换后的小波系数；

步骤3.2：计算第t次迭代计算时的残差估计值：

X^(t)表示第t次迭代计算时的小波系数；

步骤3.3：更新第t次迭代计算时的小波系数：

S^(t)＝X^(t)+ΔX^(t)

S^(t)表示更新后的第t次迭代计算时的小波系数；

步骤3.4：更新第t次迭代计算时的正则化参数β^(t)：

β^(t)＝|S^(t)|_K+1/μ

其中，μ为预设的参数，|S^(t)|表示S^(t)的幅值分量，|S^(t)|_K+1表示将S^(t)的幅值分量按照降序排列的第K+1个分量；

步骤3.5：根据S^(t)的值更新下一次迭代计算时的小波系数X^(t+1)：

F(S^(t)，μβ)＝f(g，μβ)

g＝S^(t)

其中，μ为迭代参数，将f(g，μβ)的值作为X^(t+1)的值；

步骤3.6：计算第t+1次迭代计算的误差：

Resi＝||X^(t+1)-X^(t)||_F

步骤3.7：如果Resi大于预设的阈值ε且迭代次数t小于等于预设的最大迭代次数，则迭代次数t+1，并转步骤3.2；否则停止计算。

有益效果：本发明可基于MF图像进行干涉图重建，采用本发明重建的干涉图具有更好的性能和更小的相位噪声；采用本发明获得的SAR干涉图质量更高。

附图说明

图1为本发明基于L_1/2范数正则化的SAR干涉图生成方法的实现流程图；

图2为目标场景和天线之间的几何位置关系图；

图3为本发明基于L_1/2范数正则化的SAR干涉图生成方法的迭代实现流程图；

图4为基于MF图像、基于L₁范数正则化和基于L_1/2范数正则化的仿真结果图；其中(a)为噪声标准差为1.0时基于MF图像的SAR干涉图生成仿真结果图；(b)为基于L₁范数正则化的SAR干涉图生成仿真结果图；(c)为基于L_1/2范数正则化的SAR干涉图生成仿真结果图。

具体实施方式

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

如图1所示，本实施例提出了一种基于L_1/2范数正则化的SAR干涉图生成方法，具体实现步骤如下：

步骤S1：SAR干涉相位分析

令SAR平台获取目标场景的SAR主副图像分别为y_m和y_s，它们的相位分别为φ_m和φ_s，根据InSAR成像几何关系，主副图像的相位可以分别表示为

其中，R₁和R₂分别表示两个天线距离目标场景中心的斜距，φ_s1和φ_s2分别表示y_m和 y_s的散射相位，λ_c表示接收信号的波长。目标场景和天线之间的几何位置关系如图2所示，其中B为两个天线之后的距离，H为高度。

根据InSAR理论，主副图像的相位差可以表示为

其中，ΔR＝R₂-R₁，噪声相位φ_noise是散射相位φ_s1和φ_s2的差，可以表示为

φ_noise＝φ_s2-φ_s1 (4)

由于公式3中等号右侧的表达式可以由平地相位φ_flat和地形相位φ_topo两部分组成，那么主副图像的相位差可以表示为三种类型相位的和

φ_m-φ_s＝φ_flat+φ_topo+φ_noise (5)

从公式5中可以看出，地形相位φ_topo可以通过出去平地相位φ_flat和噪声相位φ_noise来获得。因此，要想获得高质量的干涉相位图，必须尽可能地抑制平地相位和噪声相位。

步骤S2：基于L_1/2范数正则化的SAR干涉图生成模型构建

将副图像y_s的相位写成如下形式：

exp{.}为e的指数次方，j表示虚部。

的幅值＝1，|y_s|表示y_s的幅值。

将公式6重新表达为如下形式：

y＝θ_mu+n (7)

θ_m＝exp{j(φ_m-φ_flat)} (9)

u＝exp(-jφ_topo) (10)

n＝θ_mu{exp(-jφ_noise)-1} (11)

令

表示将矢量形式的y重新构写成矩阵形式的Y，Y为幅度为1的二维SAR副图像矩阵；

表示将矢量形式的θ_m重新构写成矩阵形式的Θ_m，Θ_m为含有SAR主图像相位和平地相位的矩阵；

表示将矢量形式的u重新构写成矩阵形式的U，U表示含有重建后的干涉相位的矩阵，

表示将矢量形式的n重新构写成矩阵形式的N， N表示噪声矩阵；

因此，公式7可表示为

Y＝Θ_mοU+N (12)

其中ο表示哈达玛积运算符，R(.)为离散小波逆变换。

针对式12，首先使用DWT将矩阵U稀疏表示，令W(·)和R(·)分别表示DWT(离散小波变换后)和IDWT(离散小波逆变换，Discrete wavelet inversion transform，简称IDWT)，则X＝W(U)，矩阵X表示矩阵U离散小波变换后的小波系数。那么，式12 又可以表示为

Y＝Θ_mοR(X)+N (13)

通过解决L_1/2范数正则化问题来重建干涉图像。

为与恢复后的干涉图像有关的小波系数矩阵，argmin{.}表示取最小值，β表示正则化，

表示斐波那契范数。

步骤S3：阈值迭代

针对公式14中的L_1/2范数正则化问题，本实施例采用阈值迭代算法来实现对干涉图的重建。阈值迭代算法输入副图像Y以及包含主图像相位和平地相位的相位矩阵Θ_m。设小波系数的初始值X⁽⁰⁾＝0；所述小波系数为SAR主、副图像的干涉相位经过DWT变换后的小波系数，迭代参数为μ，通常设置为常数，误差参数为ε，最大迭代步数为T_max。当满足条件t≤T_max且迭代误差Resi＞ε时，如图3所示，执行下列步骤。

步骤S31：残差估算

符号(·)^H表示共轭转置算子，X^(t)表示第t次迭代计算时的小波系数。

步骤S32：更新第t次迭代计算时的小波系数：

S^(t)＝X^(t)+ΔX^(t) (16)

其中，S^(t)表示更新后的第t次迭代计算时的小波系数。

步骤S33：史新第t次迭代计算时的正则化参数β^(t)：

β^(t)＝|S^(t)|_K+1/μ (17)

其中，|S^(t)|表示S^(t)的幅值分量，|S^(t)|_K+1表示将S^(t)的幅值分量按照降序排列的第 K+1个分量；

步骤S34：Half阈值收缩

阈值算子F(S^(t)，μβ)为

F(S^(t)，μβ)＝f(g，μβ) (18)

阈值函数f(g，μβ)表示为

g＝S^(t)

将f(g，μβ)的值作为更新后的X^(t+1)的值，并将该X^(t+1)的值作为下一次迭代计算时的小波系数；既f(g，μβ)＝X^(t+1)。

步骤S35：计算t+1步迭代误差

Resi＝||X^(t+1)-X^(t)||_F (21)

若满足条件t≤T_max且Resi＞ε，则继续迭代，即t＝t+1，重复上述步骤，ε为预设的阈值。若不满足条件，则迭代过程结束，输出X^(t)，

最后，通过IDWT并取

逆变换后的角度值作为恢复后的SAR干涉图：

作为恢复后的SAR干涉图，angle(.)表示取角度值。

下面通过实验对本发明中提供的一种基于L_1/2范数正则化的SAR干涉图生成方法进行验证。图4为仿真干涉图像，对其加入标准差为1.0的高斯白噪声，得到如图4中的 (a)所示的干涉图，它代表基于MF图像处理得到的干涉图。图4中的(b)为基于L₁范数正则化的SAR干涉图；图4中的(c)为基于L_1/2范数正则化的SAR干涉图。对比基于 L₁和L_1/2范数正则化的SAR干涉图像质量，不难看出基于L₁和L_1/2范数正则化的方法均能够有效抑制相位噪声，本发明所提出的L_1/2范数正则化方法能够获得质量更好的干涉图，使其拥有比L₁范数正则化方法更少的残差点和更小的均方误差值。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

Claims (4)

1.一种基于L_1/2范数正则化的SAR干涉图生成方法，其特征在于：具体包括如下步骤：

步骤1：采用SAR获取同一目标场景的SAR主图像和SAR副图像，通过对副图像复矩阵取相位值获取SAR副图像的相位；

步骤2：基于SAR主、副图像的干涉相位，以及通过解决L_1/2范数正则化问题建立SAR干涉图的生成模型；

步骤3：采用阈值迭代算法求解L_1/2范数正则化问题，得到与恢复后的SAR干涉图有关的小波系数矩阵

步骤4：对步骤3求得的

进行离散小波逆变换，取变换后的角度值作为恢复后的干涉图像。

2.根据权利要求1所述的一种基于L_1/2范数正则化的SAR干涉图生成方法，其特征在于：所述步骤1中SAR副图像的相位y为：

θ_m＝exp{j(φ_m-φ_flat)}

u＝exp(-jφ_topo)

n＝θ_mu{exp(-jφ_noise)-1}

其中，φ_m为SAR主图像y_m的相位，φ_flat为平地相位，φ_topo为地形相位；φ_noise＝φ_s2-φ_s1，φ_s1和φ_s2分别表示y_m和y_s的散射相位，y_s表示SAR副图像,exp{.}为e的指数次方，j表示虚部。

3.根据权利要求2所述的一种基于L_1/2范数正则化的SAR干涉图生成方法，其特征在于：所述步骤2中SAR干涉图的生成模型为：

其中，

表示将矢量形式的y重新构写成矩阵形式的Y,Y为幅度为1的二维SAR副图像矩阵；

表示将矢量形式的θ_m重新构写成矩阵形式的Θ_m，Θ_m为含有SAR主图像相位和平地相位的矩阵；R(.)为离散小波逆变换；

表示哈达玛积；X＝W(U)，

表示将矢量形式的u重新构写成矩阵形式的U,U表示含有重建后的SAR主、副图像的干涉相位的矩阵；

表示将矢量形式的n重新构写成矩阵形式的N,N表示噪声矩阵，β表示正则化参数，

表示斐波那契范数，argmin{.}表示取最小值,W(.)为离散小波变换。

4.根据权利要求3所述的一种基于L_1/2范数正则化的SAR干涉图生成方法，其特征在于，所述步骤3具体为：

步骤3.1：设置小波系数的初始值X⁽⁰⁾＝0；所述小波系数为SAR主、副图像的干涉相位经过DWT变换后的小波系数；

步骤3.2：计算第t次迭代计算时的残差估计值：

X^(t)表示第t次迭代计算时的小波系数；

步骤3.3：更新第t次迭代计算时的小波系数：

S^(t)＝X^(t)+ΔX^(t)

S^(t)表示更新后的第t次迭代计算时的小波系数；

步骤3.4：更新第t次迭代计算时的正则化参数β^(t)：

β^(t)＝|S^(t)|_K+1/μ

其中，μ为预设的参数，|S^(t)|表示S^(t)的幅值分量，|S^(t)|_K+1表示将S^(t)的幅值分量按照降序排列的第K+1个分量；

步骤3.5：根据S^(t)的值更新下一次迭代计算时的小波系数X^(t+1):

F(S^(t),μβ)＝f(g,μβ)

g＝S^(t)

其中，μ为迭代参数，将f(g,μβ)的值作为X^(t+1)的值；

步骤3.6：计算第t+1次迭代计算的误差：

Resi＝||X^(t+1)-X^(t)||_F

步骤3.7：如果Resi大于预设的阈值ε且迭代次数t小于等于预设的最大迭代次数，则迭代次数t+1，并转步骤3.2；否则停止计算。

Images