一种基于L的制作方法

一种基于l
1/2
范数正则化的sar干涉图生成方法
技术领域
1.本发明属于干涉合成孔径雷达成像和稀疏信号处理技术领域。


背景技术：

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


技术实现要素：

3.发明目的：为了解决上述现有技术存在的问题，本发明提供了一种基于l
1/2
范数正则化的sar干涉图生成方法。
4.技术方案：本发明提供了一种基于l
1/2
范数正则化的sar干涉图生成方法，具体包括如下步骤：
5.步骤1：采用sar获取同一目标场景的sar主图像和sar副图像，通过对副图像复矩阵取相位值获取sar副图像的相位；
6.步骤2：基于sar主、副图像的干涉相位，以及通过解决l
1/2
范数正则化问题建立 sar干涉图的生成模型；
7.步骤3：采用阈值迭代算法求解l
1/2
范数正则化问题，得到与恢复后的sar干涉图有关的小波系数矩阵
8.步骤4：对步骤3求得的进行离散小波逆变换，取变换后的角度值作为恢复后的干涉图像。
9.进一步的，所述步骤1中sar副图像的相位y为：
[0010][0011][0012]
θm＝exp{j(φ
m-φ
flat
)}
[0013]
u＝exp(-jφ
topo
)
[0014]
n＝θmu{exp(-jφ
noise
)-1}
[0015]
其中，φm为sar主图像ym的相位，φ
flat
为平地相位，φ
topo
为地形相位；φ
noise
＝φ
s2-φ
s1
，φ
s1
和φ
s2
分别表示ym和ys的散射相位，ys表示sar副图像，exp{.}为e的指数次方，
j 表示虚部。
[0016]
进一步的，所述步骤2中sar干涉图的生成模犁为：
[0017][0018]
y＝θmοr(x) n
[0019]
其中，表示将矢量形式的y重新构写成矩阵形式的y，y为幅度为1的二维sar 副图像矩阵；表示将矢量形式的θm重新构写成矩阵形式的θm，θm为含有sar 主图像相位和平地相位的矩阵；r(.)为离散小波逆变换；。表示哈达玛积；x＝w(u)，表示将矢量形式的u重新构写成矩阵形式的u，u表示含有重建后的sar主、副图像的干涉相位的矩阵；表示将矢量形式的n重新构写成矩阵形式的n，n表示噪声矩阵，β表示正则化参数，表示斐波那契范数，argmin{.}表示取最小值，w(.) 为离散小波变换。
[0020]
进一步的，所述步骤3具体为：
[0021]
步骤3.1：设置小波系数的初始值x
(0)
＝0；所述小波系数为sar主、副图像的干涉相位经过dwt变换后的小波系数；
[0022]
步骤3.2：计算第t次迭代计算时的残差估计值：
[0023][0024]
x
(t)
表示第t次迭代计算时的小波系数；
[0025]
步骤3.3：更新第t次迭代计算时的小波系数：
[0026]s(t)
＝x
(t)
δx
(t)
[0027]s(t)
表示更新后的第t次迭代计算时的小波系数；
[0028]
步骤3.4：更新第t次迭代计算时的正则化参数β
(t)
：
[0029]
β
(t)
＝|s
(t)
|
k 1
/μ
[0030]
其中，μ为预设的参数，|s
(t)
|表示s
(t)
的幅值分量，|s
(t)
|
k 1
表示将s
(t)
的幅值分量按照降序排列的第k 1个分量；
[0031]
步骤3.5：根据s
(t)
的值更新下一次迭代计算时的小波系数x
(t 1)
：
[0032]
f(s
(t)
，μβ)＝f(g，μβ)
[0033][0034][0035]
g＝s
(t)
[0036]
其中，μ为迭代参数，将f(g，μβ)的值作为x
(t 1)
的值；
[0037]
步骤3.6：计算第t 1次迭代计算的误差：
[0038]
resi＝||x
(t 1)-x
(t)
||f[0039]
步骤3.7：如果resi大于预设的阈值ε且迭代次数t小于等于预设的最大迭代次数，
则迭代次数t 1，并转步骤3.2；否则停止计算。
[0040]
有益效果：本发明可基于mf图像进行干涉图重建，采用本发明重建的干涉图具有更好的性能和更小的相位噪声；采用本发明获得的sar干涉图质量更高。
附图说明
[0041]
图1为本发明基于l
1/2
范数正则化的sar干涉图生成方法的实现流程图；
[0042]
图2为目标场景和天线之间的几何位置关系图；
[0043]
图3为本发明基于l
1/2
范数正则化的sar干涉图生成方法的迭代实现流程图；
[0044]
图4为基于mf图像、基于l1范数正则化和基于l
1/2
范数正则化的仿真结果图；其中 (a)为噪声标准差为1.0时基于mf图像的sar干涉图生成仿真结果图；(b)为基于l1范数正则化的sar干涉图生成仿真结果图；(c)为基于l
1/2
范数正则化的sar干涉图生成仿真结果图。
具体实施方式
[0045]
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
[0046]
如图1所示，本实施例提出了一种基于l
1/2
范数正则化的sar干涉图生成方法，具体实现步骤如下：
[0047]
步骤s1：sar干涉相位分析
[0048]
令sar平台获取目标场景的sar主副图像分别为ym和ys，它们的相位分别为φm和φs，根据insar成像几何关系，主副图像的相位可以分别表示为
[0049][0050][0051]
其中，r1和r2分别表示两个天线距离目标场景中心的斜距，φ
s1
和φ
s2
分别表示ym和 ys的散射相位，λc表示接收信号的波长。目标场景和天线之间的几何位置关系如图2所示，其中b为两个天线之后的距离，h为高度。
[0052]
根据insar理论，主副图像的相位差可以表示为
[0053][0054]
其中，δr＝r
2-r1，噪声相位φ
noise
是散射相位φ
s1
和φ
s2
的差，可以表示为
[0055]
φ
noise
＝φ
s2-φ
s1
ꢀꢀꢀ
(4)
[0056]
由于公式3中等号右侧的表达式可以由平地相位φ
flat
和地形相位φ
topo
两部分组成，那么主副图像的相位差可以表示为三种类型相位的和
[0057]
φ
m-φs＝φ
flat
φ
topo
φ
noise
ꢀꢀꢀ
(5)
[0058]
从公式5中可以看出，地形相位φ
topo
可以通过出去平地相位φ
flat
和噪声相位φ
noise
来获得。因此，要想获得高质量的干涉相位图，必须尽可能地抑制平地相位和噪声相位。
[0059]
步骤s2：基于l
1/2
范数正则化的sar干涉图生成模型构建
[0060]
将副图像ys的相位写成如下形式：
[0061][0062]
exp{.}为e的指数次方，j表示虚部。
[0063]
的幅值＝1，|ys|表示ys的幅值。
[0064]
将公式6重新表达为如下形式：
[0065]
y＝θmu n
ꢀꢀꢀ
(7)
[0066][0067]
θm＝exp{j(φ
m-φ
flat
)}
ꢀꢀꢀ
(9)
[0068]
u＝exp(-jφ
topo
)
ꢀꢀꢀ
(10)
[0069]
n＝θmu{exp(-jφ
noise
)-1}
ꢀꢀꢀ
(11)
[0070]
令表示将矢量形式的y重新构写成矩阵形式的y，y为幅度为1的二维sar副图像矩阵；表示将矢量形式的θm重新构写成矩阵形式的θm，θm为含有sar主图像相位和平地相位的矩阵；表示将矢量形式的u重新构写成矩阵形式的u，u表示含有重建后的干涉相位的矩阵，表示将矢量形式的n重新构写成矩阵形式的n， n表示噪声矩阵；
[0071]
因此，公式7可表示为
[0072]
y＝θmοu n
ꢀꢀꢀ
(12)
[0073]
其中ο表示哈达玛积运算符，r(.)为离散小波逆变换。
[0074]
针对式12，首先使用dwt将矩阵u稀疏表示，令w(
·
)和r(
·
)分别表示dwt(离散小波变换后)和idwt(离散小波逆变换，discrete wavelet inversion transform，简称idwt)，则x＝w(u)，矩阵x表示矩阵u离散小波变换后的小波系数。那么，式12 又可以表示为
[0075]
y＝θmοr(x) n
ꢀꢀꢀ
(13)
[0076]
通过解决l
1/2
范数正则化问题来重建干涉图像。
[0077][0078]
为与恢复后的干涉图像有关的小波系数矩阵，argmin{.}表示取最小值，β表示正则化，表示斐波那契范数。
[0079]
步骤s3：阈值迭代
[0080]
针对公式14中的l
1/2
范数正则化问题，本实施例采用阈值迭代算法来实现对干涉图的重建。阈值迭代算法输入副图像y以及包含主图像相位和平地相位的相位矩阵θm。设小波系数的初始值x
(0)
＝0；所述小波系数为sar主、副图像的干涉相位经过dwt变换后的小波系数，迭代参数为μ，通常设置为常数，误差参数为ε，最大迭代步数为t
max
。当满足条件t≤t
max
且迭代误差resi＞ε时，如图3所示，执行下列步骤。
[0081]
步骤s31：残差估算
[0082]
[0083]
符号(
·
)h表示共轭转置算子，x
(t)
表示第t次迭代计算时的小波系数。
[0084]
步骤s32：更新第t次迭代计算时的小波系数：
[0085]s(t)
＝x
(t)
δx
(t)
ꢀꢀꢀ
(16)
[0086]
其中，s
(t)
表示更新后的第t次迭代计算时的小波系数。
[0087]
步骤s33：史新第t次迭代计算时的正则化参数β
(t)
：
[0088]
β
(t)
＝|s
(t)
|
k 1
/μ
ꢀꢀꢀ
(17)
[0089]
其中，|s
(t)
|表示s
(t)
的幅值分量，|s
(t)
|
k 1
表示将s
(t)
的幅值分量按照降序排列的第 k 1个分量；
[0090]
步骤s34：half阈值收缩
[0091]
阈值算子f(s
(t)
，μβ)为
[0092]
f(s
(t)
，μβ)＝f(g，μβ)
ꢀꢀꢀ
(18)
[0093]
阈值函数f(g，μβ)表示为
[0094][0095]
g＝s
(t)
[0096][0097]
将f(g，μβ)的值作为更新后的x
(t 1)
的值，并将该x
(t 1)
的值作为下一次迭代计算时的小波系数；既f(g，μβ)＝x
(t 1)
。
[0098]
步骤s35：计算t 1步迭代误差
[0099]
resi＝||x
(t 1)-x
(t)
||fꢀꢀꢀ
(21)
[0100]
若满足条件t≤t
max
且resi＞ε，则继续迭代，即t＝t 1，重复上述步骤，ε为预设的阈值。若不满足条件，则迭代过程结束，输出x
(t)
，
[0101]
最后，通过idwt并取逆变换后的角度值作为恢复后的sar干涉图：
[0102][0103]
作为恢复后的sar干涉图，angle(.)表示取角度值。
[0104]
下面通过实验对本发明中提供的一种基于l
1/2
范数正则化的sar干涉图生成方法进行验证。图4为仿真干涉图像，对其加入标准差为1.0的高斯白噪声，得到如图4中的 (a)所示的干涉图，它代表基于mf图像处理得到的干涉图。图4中的(b)为基于l1范数正则化的sar干涉图；图4中的(c)为基于l
1/2
范数正则化的sar干涉图。对比基于 l1和l
1/2
范数正则化的sar干涉图像质量，不难看出基于l1和l
1/2
范数正则化的方法均能够有效抑制相位噪声，本发明所提出的l
1/2
范数正则化方法能够获得质量更好的干涉图，使其拥有比l1范数正则化方法更少的残差点和更小的均方误差值。
[0105]
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛
盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。