CN105931195A - 一种合成孔径雷达图像噪声抑制方法 - Google Patents

Abstract

本发明公开了一种合成孔径雷达图像相干噪声的抑制方法，包括以下步骤：(1)将合成孔径雷达原始图像，取参考图像块；对于每一个参考图象块，进行相似块匹配，找出备选相似块并按列堆积成噪声污染的观测矩阵；(2)根据相干噪声统计特性建立能量泛函模型；(3)求解得到低秩矩阵L：(4)对于每一个参考图像块，将其备选相似图像块对应的低秩矩阵L的各列堆成替换图像块；对于同一图像块，将多个替换图像块进行加权平均处理后得到噪声抑制结果图像。本发明提供的方法，考虑到高斯型噪声和泊松型噪声，兼顾相干斑噪声去除和图像细节保存。

Description

一种合成孔径雷达图像噪声抑制方法

技术领域

本发明属于图像处理领域，更具体地，涉及一种合成孔径雷达图像噪声抑制方法。

背景技术

合成孔径雷达(即SAR)成像是一种回波相干叠加成像，不可避免会引入相干斑噪声。严重影响SAR图像的视觉质量，并极大制约了后续的特征提取、目标识别等解译处理。相干斑噪声是一种乘性噪声，其抑制困难在于有效抑制相干斑噪声的同时较完整地保存图像原有的结构信息。

SAR图像相干斑抑制算法按照所采用的图像数据格式可分为针对复数数据的方法和针对实数数据的方法两大类。

目前针对复数SAR数据的相干斑噪声抑制算法相对较少，主要代表算法有如下两类：第一类是建立在吉布斯先验模型基础上的贝叶斯估计方法；第二类是以最大熵原理为基础的规则化方法。

针对实数域的相干斑抑制算法则比较多，大致可分为如下三类：(1)基于估计理论的相干斑噪声抑制算法。(2)基于变换域的相干斑噪声抑制算法。(3)基于变分正则化理论的相干斑噪声抑制算法，将相干斑噪声抑制问题看作反问题，通过引入正则化约束来确保恢复结果的存在、唯一和受噪声干扰较小，因而正则化方法对于相干斑噪声抑制来说是一种行之有效的方法。传统矩阵低秩近似问题能量泛函建模为关于L的核范数最小化问题：

${\min }_L\frac{1}{2}\left|\right|P-L|{|}_F^2+\mathrm{\ρ}|\left|L\right|{|}_{*}$

该模型的数据项仅有这一项，它假设噪声的概率分布是高斯型的，然而实际情况并非尽然，因此该方法效果有待改进。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种合成孔径雷达图像相干噪声的抑制方法，其目的在于通过同时考虑高斯型噪声和泊松型噪声，将相干斑噪声抑制问题表示为从受噪声污染的相似块观测矩阵估计其低秩近似矩阵的逆问题，利用了图像具有非局部结构自相似性这一特性，在尽可能滤除噪声的同时保留细节，由此解决现有的噪声抑制方法噪声滤除效果不理想或者细节损失较严重的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种合成孔径雷达图像相干噪声的抑制方法，包括以下步骤：

(1)将合成孔径雷达原始图像，按照预设的边长和步长，取参考图像块；对于每一个参考图象块，在其搜索框内计算其与各相似图像块的相似度，从而进行相似块匹配；根据预设的相似度阈值，找出备选相似块并按列堆积成噪声污染的观测矩阵；

(2)根据相干噪声统计特性建立能量泛函模型如下：

$E\left(L\right)={\mathrm{\λ}}_1\left|\right|P-L|{|}_1+\frac{{\mathrm{\λ}}_2}{2}||P-L|{|}_F^2+\mathrm{\ρ}\left|\right|L|{|}_{*}$

其中，P(x,y)为噪声污

染的观测矩阵，L(x,y)为原始的低秩矩阵，E(L)表示关于待恢复低秩矩阵L的能量泛函；对泊松噪声建模的数据项λ₁||P-L||₁，对高斯噪声建模的数据项及对矩阵L的低秩正则化项ρ||L||_*；

(3)求解步骤(2)获得的关于待恢复低秩矩阵L的能量泛函模型，得到低秩矩阵L：

(4)对于每一个参考图像块，将其备选相似图像块对应的低秩矩阵L的各列堆成替换图像块；对于同一图像块，作为多个参考图像块的备选相似块可能得到多个替换图像块，将所述多个替换图像块进行加权平均处理后得到噪声抑制结果图像。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其步骤(1)所述相似度根据欧氏距离、曼哈顿距离或马氏距离计算。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其步骤(3)所述关于待恢复低秩矩阵L的能量泛函模型，可根据交换方向乘子算法或分裂Bregman算法求解。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其采用交换方向乘子算法具体步骤为：

引入辅助变量Z，得到能量泛函模型如下：

$E(L,Z)={\mathrm{\&lambda;}}_1\left|\right|Z|{|}_1+\frac{{\mathrm{\&lambda;}}_2}{2}||P-L|{|}_F^2+\mathrm{\&rho;}\left|\right|L|{|}_{*}+<Y,Z-P+L>+\frac{{\mathrm{\&lambda;}}_3}{2}||Z-P+L|{|}_F^2$

其中，辅助变量Z作为替代项替换第一数据项的P-L，是Z替换P-L引入的二次惩罚项，Y是拉格朗日乘子，λ₃是惩罚参数；

将引入辅助变量Z后的能量泛函模型分解为三个子问题交替迭代求解，得到低秩近似后的矩阵L，所述三个子问题分别为：

第一子问题：固定Z和Y，求解L；

第二子问题：固定L和Y，求解Z；

第三子问题：固定L和Z，求解Y。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其所述三个子问题的对应迭代目标函数为：

$\left\{\begin{array}{c}\hat{L}=D_{\mathrm{\&rho;}/{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}(P-\frac{{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}(Z+\frac{1}{{\mathrm{\&lambda;}}_3}Y))={\mathrm{US}}_{\mathrm{\&rho;}/{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}\left(\mathrm{\&Sigma;}\right)V^T.\\ \hat{Z}=S_{{\mathrm{\&lambda;}}_1/{\mathrm{\&lambda;}}_3}(P-L-\frac{1}{{\mathrm{\&lambda;}}_3}Y).\\ \hat{Y}=Y+\mathrm{\&mu;}(Z-P+L);\widehat{\mathrm{\&mu;}}=\mathrm{\&alpha;}\mathrm{\&mu;},(\mathrm{\&alpha;}>1).\end{array}\right.$

其中，S()为软阈值收缩算子，其公式为D()为奇异值阈值收缩算子，其公式为D_b(a)＝US_b(Σ)V^T,a＝UΣV^T。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其步骤(1)所述边长为8像素。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其步骤(1)所述步长为6像素。

优选地，所述合成孔径雷达图像相干噪声的抑制方法，其步骤(1)所述搜索框半径为21像素。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明与现有技术相比，具有如下优点：

兼顾相干斑噪声去除和图像细节保存。本发明考虑到高斯型噪声和泊松型噪声，同时将相干斑噪声抑制问题表示为从相干斑噪声图像中估计真实图像的逆问题，利用了图像自身具有非局部结构自相似性这一特性知识。既能有效地去除相干斑噪声，又能较好的保存图像细节。

优选方案，算法收敛速度快。本发明利用交替方向乘子(ADMM)方法进行数值优化求解，有效地加快算法的收敛速度。

优选方案，噪声适应性强。本发明采用混合数据项来捕捉不同类型噪声的统计特性，通过调节算法接口参数可适应多类型混合噪声。

附图说明

图1是本发明提供的合成孔径雷达想干噪声抑制方法的流程示意图；

图2是本发明实施例提供的合成孔径雷达相干噪声抑制方法的流程示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明提供的合成孔径雷达图像相干噪声抑制方法，如图1所示，包括以下步骤：

(1)将合成孔径雷达原始图像，按照预设的边长和步长，取参考图像块；对于每一个参考图象块，在其搜索框内计算其与各相似图像块的相似度，从而进行相似块匹配；根据预设的相似度阈值，找出备选相似块并按列堆积成噪声污染的观测矩阵；

所述相似度，可根据欧氏距离、曼哈顿距离或马氏距离计算。

(2)根据相干噪声统计特性建立待恢复低秩矩阵L的能量泛函模型如下：

首先，对观测矩阵退化建模并建立低秩约束的噪声抑制能量泛函模型：

噪声污染的观测矩阵P退化建模为：

P＝L+N

其中，L为待估计低秩矩阵，N为噪声矩阵。

进而推导出关于待恢复低秩矩阵L的能量泛函模型：

$\begin{array}{c}E\left(L\right)={\mathrm{\&lambda;}}_1\left|\right|P-L|{|}_1+\frac{{\mathrm{\&lambda;}}_2}{2}||P-L|{|}_F^2+\mathrm{\&rho;}\left|\right|L|{|}_{*}\\ ={\mathrm{\&lambda;}}_1\underset{\mathrm{\&Omega;}}{\&Integral;}|P-L|dxdy+\frac{{\mathrm{\&lambda;}}_2}{2}\underset{\mathrm{\&Omega;}}{\&Integral;}{(P-L)}^{2}dxdy+\mathrm{\&rho;}\underset{i}{\&Integral;}{\mathrm{\&sigma;}}_i\left(L\right)d_i\end{array}$

其中，P(x,y)为噪声污染的观测矩阵，L(x,y)为原始的低秩矩阵，E(L)表示关于待恢复低秩矩阵L的能量泛函，λ₁，λ₂，ρ是正则化参数，Ω表示图像像素坐标集合。能量泛函第一项表示对泊松噪声建模的数据项，第二项表示对高斯噪声建模的数据项，第三项表示对矩阵L的低秩正则化项，σ_i(L)表示矩阵L的第i个奇异值。第二行为具体展开表达式。

本发明在传统低秩近似模型基础上加入数据项λ₁||P-L||₁，将泊松型噪声一并考虑进来，可有效处理合成孔径雷达成像这类受混合型噪声污染的图像。

所建立的低秩约束的噪声抑制能量泛函模型包括如下三项：对泊松噪声建模的数据项λ₁||P-L||₁，对高斯噪声建模的数据项及对矩阵L的低秩正则化项ρ||L||_*。混合数据项的使用提升了算法对不同类型噪声的稳定性。

(3)求解步骤(2)获得的关于待恢复低秩矩阵L的能量泛函模型，得到低秩矩阵L。

所述关于待恢复低秩矩阵L的能量泛函模型，可根据交换方向乘子算法(ADMM算法)或分裂Bregman算法求解。

优选采用ADMM算法，具体如下：

为便于低秩矩阵L的后续求解，引入辅助变量Z分离出数据项λ₁||P-L||₁中的L，得到能量泛函模型如下：

$E(L,Z)={\mathrm{\&lambda;}}_1\left|\right|Z|{|}_1+\frac{{\mathrm{\&lambda;}}_2}{2}||P-L|{|}_F^2+\mathrm{\&rho;}\left|\right|L|{|}_{*}+<Y,Z-P+L>+\frac{{\mathrm{\&lambda;}}_3}{2}||Z-P+L|{|}_F^2$

其中，辅助变量Z作为替代项替换第一数据项的P-L，是Z替换P-L引入的二次惩罚项，Y是拉格朗日乘子，λ₃是惩罚参数，用于约束引入的辅助变量。

将步骤引入辅助变量Z后的能量泛函模型分解为三个子问题交替迭代求解，得到低秩近似后的矩阵L，所述三个子问题分别为：

第一子问题：固定Z和Y，求解L；得到：

利用奇异值软阈值公式可得：

$\hat{L}=D_{\mathrm{\&rho;}/{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}(P-\frac{{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}(Z+\frac{1}{{\mathrm{\&lambda;}}_3}Y\left)\right)={\mathrm{US}}_{\mathrm{\&rho;}/{\mathrm{\&lambda;}}_2+{\mathrm{\&lambda;}}_3}\left(\mathrm{\&Sigma;}\right)V^T,$

其中，

第二子问题：固定L和Y，求解Z；得到：

利用软阈值公式可得：

第三子问题：固定L和Z，求解Y；得到：

$\hat{Y}=Y+\mathrm{\&mu;}(Z-P+L)$

$\widehat{\mathrm{\&mu;}}=\mathrm{\&alpha;}\mathrm{\&mu;},(\mathrm{\&alpha;}>1)$

其中，μ和α为两个常数用于更新拉格朗日乘子Y，分别表示L，Z，Y，u的更新结果，S()为软阈值收缩算子，其公式为D()为奇异值阈值收缩算子，其公式为D_b(a)＝US_b(Σ)V^T,a＝UΣV^T。

三个子问题迭代更新具体步骤如下：

A、初始化辅助变量Z和拉格朗日乘子Y为零矩阵，迭代次数k＝1；

B、将辅助变量Z、拉格朗日乘子Y代入子第一子问题，求解得到低秩近似后的矩阵L^k；

C、将低秩近似矩阵L^k代入第二子问题求解Z和第三子问题求解Y对应的迭代目标函数，得到Z^k和Y^k；

D、若或迭代次数k达到最大值，ε为预设阈值，则L^k为最终低秩近似矩阵，结束；否则，进入步骤E；

E、更新k＝k+1，进入步骤B。

(4)对于每一个参考图像块，将其备选相似图像块对应的低秩矩阵L的各列堆成替换图像块；对于同一图像块，作为多个参考图像块的备选相似块可能得到多个替换图像块，将所述多个替换图像块进行加权平均处理后得到噪声抑制结果图像。

所述加权平均处理，其权重根据替换图像块相应的低秩矩阵L的秩确定，其具体计算公式如下：

wei＝(maxRank-r+1)/maxRank,

maxRank＝min(row,col)；

其中，wei为权重，row、col为L的行、列数，r为L的秩。

以下为实施例：

本发明根据相干斑噪声的特点，提出包括L₁范数数据项λ₁||P-L||₁、L₂范数数据项及核范数正则化项ρ||L||_*的能量泛函，通过交替方向乘子法(ADMM)求解实现。本发明实施例提供的合成孔径雷达图像相干噪声抑制方法，如图2所示，包括以下步骤：

(1)将合成孔径雷达原始图像按横、纵间隔均为6个像素的步长，取出8×8的参考图像块，利用图像的非局部结构自相似性质，在半径为21的搜素窗内对各参考图像块计算相似度进行相似块匹配，取出最相似的80个相似图像块并按列堆积成64×80的噪声污染观测矩阵P；

(2)根据相干斑噪声的统计特建立能量泛函模型：

$\begin{array}{c}E\left(L\right)={\mathrm{\&lambda;}}_1\left|\right|P-L|{|}_1+\frac{{\mathrm{\&lambda;}}_2}{2}||P-L|{|}_F^2+\mathrm{\&rho;}\left|\right|L|{|}_{*}\\ ={\mathrm{\&lambda;}}_1\underset{\mathrm{\&Omega;}}{\&Integral;}|P-L|dxdy+\frac{{\mathrm{\&lambda;}}_2}{2}\underset{\mathrm{\&Omega;}}{\&Integral;}{(P-L)}^{2}dxdy+\mathrm{\&rho;}\underset{i}{\&Integral;}{\mathrm{\&sigma;}}_i\left(L\right)d_i\end{array}$

其中，P(x,y)为噪声污染的观测矩阵，L(x,y)为原始的低秩矩阵，E(L)表示关于待恢复低秩矩阵L的能量泛函，λ₁，λ₂，ρ是正则化参数，Ω表示图像像素坐标集合。能量泛函第一项表示对泊松噪声建模的数据项，第二项表示对高斯噪声建模的数据项，第三项表示对矩阵L的低秩正则化项，σ_i(L)表示矩阵L的第i个奇异值。第二行为具体展开表达式。

(3)交替方向乘子法(ADMM)求解引入辅助变量Z,建立新的能量泛函：

$E(L,Z)={\mathrm{\&lambda;}}_1\left|\right|Z|{|}_1+\frac{{\mathrm{\&lambda;}}_2}{2}||P-L|{|}_F^2+\mathrm{\&rho;}\left|\right|L|{|}_{*}+<Y,Z-P+L>+\frac{{\mathrm{\&lambda;}}_3}{2}||Z-P+L|{|}_F^2$

其中，辅助变量Z作为替代项替换第一数据项的P-L，是Z替换P-L引入的二次惩罚项，Y是拉格朗日乘子，λ₃是惩罚参数，用于约束引入的辅助变量。

(4)将步骤(3)获得的能量泛函模型分解为三个子问题交替迭代求解，得到低秩近似后的矩阵L，其中，

(a)第一子问题为：固定Z、Y，求解L，得到：

利用奇异值软阈值公式可得：其中，

(b)第二子问题为：固定L、Y，求解Z，得到：

利用软阈值公式可得：

(c)第三子问题为：固定L、Z，求解Y，得到：

$\hat{Y}=Y+\mathrm{\&mu;}(Z-P+L)$

$\widehat{\mathrm{\&mu;}}=\mathrm{\&alpha;}\mathrm{\&mu;},(\mathrm{\&alpha;}>1)$

其中，μ和α为两个常数用于更新拉格朗日乘子Y，分别表示L，Z，Y，u的更新结果，S()为软阈值收缩算子，其公式为D()为奇异值阈值收缩算子，其公式为D_b(a)＝US_b(Σ)V^T,a＝UΣV^T。

三个子问题迭代更新具体步骤如下：

A、初始化辅助变量Z和拉格朗日乘子Y为零矩阵，迭代次数k＝1；

B、将辅助变量Z、拉格朗日乘子Y代入子第一子问题，求解得到低秩近似后的矩阵L^k；

C、将低秩近似矩阵L^k代入第二子问题求解Z和第三子问题求解Y对应的迭代目标函数，得到Z^k和Y^k；

D、若或迭代次数k达到最大值，ε为预设阈值，则L^k为最终低秩近似矩阵，结束；否则，进入步骤E；

E、更新k＝k+1，进入步骤B。

(5)将各参考图像块对应的低秩化矩阵L的各列重新堆成8×8的图像块，置换所述参考图像块，并进行加权平均处理，输出噪声抑制结果图像。所述加权平均处理，其权重根据替换图像块相应的低秩矩阵L的秩确定，其具体计算公式如下：

wei＝(maxRank-r+1)/maxRank,

maxRank＝min(row,col)；

其中，wei为权重，row、col为L的行、列数，r为L的秩。

在本方法建模过程中引入了数据项权重参数λ₁和λ₂，正则化项参数ρ，以及在ADMM优化求解过程中引入的惩罚参数λ₃。本发明针对不同强度的相干斑噪声图像示例列出多组与之对应的参数，其只作为示例不作为对本发明的限制。算法外循环最大迭代次数设为20次，内循环最大迭代次数设为5次，相似块尺寸w＝8，相似块个数nblk＝80，参考块采样间隔step＝6，搜索窗口半径s＝21，拉格朗日乘子更新系数μ和α分别设为1和1.05。在轻噪声下的λ₁＝0.01，λ₂＝0.8，λ₃＝3，ρ＝0.1；在中度噪声下的λ₁＝0.01，λ₂＝0.8，λ₃＝3，ρ＝0.1；在重度噪声下的λ₁＝0.01，λ₂＝0.8，λ₃＝3，ρ＝0.2。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种合成孔径雷达图像相干噪声的抑制方法，其特征在于，包括以下步骤：

(1)将合成孔径雷达原始图像，按照预设的边长和步长，取参考图像块；对于每一个参考图象块，在其搜索框内计算其与各相似图像块的相似度，从而进行相似块匹配；根据预设的相似度阈值，找出备选相似块并按列堆积成噪声污染的观测矩阵；

(2)根据相干噪声统计特性建立能量泛函模型如下：

其中，P(x,y)为噪声污

染的观测矩阵，L(x,y)为原始的低秩矩阵，E(L)表示关于待恢复低秩矩阵L的能量泛函；对泊松噪声建模的数据项λ₁||P-L||₁，对高斯噪声建模的数据项 及对矩阵L的低秩正则化项ρ||L||_*；

(3)求解步骤(2)获得的关于待恢复低秩矩阵L的能量泛函模型，得到低秩矩阵L：

(4)对于每一个参考图像块，将其备选相似图像块对应的低秩矩阵L的各列堆成替换图像块；对于同一图像块，作为多个参考图像块的备选相似块可能得到多个替换图像块，将所述多个替换图像块进行加权平均处理后得到噪声抑制结果图像。

2.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，步骤(1)所述相似度根据欧氏距离、曼哈顿距离或马氏距离计算。

3.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，步骤(3)所述关于待恢复低秩矩阵L的能量泛函模型，可根据交换方向乘子算法或分裂Bregman算法求解。

4.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，采用交换方向乘子算法具体步骤为：

引入辅助变量Z，得到能量泛函模型如下：

其中，辅助变量Z作为替代项替换第一数据项的是Z替换P-L引入的二次惩罚项，Y是拉格朗日乘子，λ₃是惩罚参数；

将步骤引入辅助变量Z后的能量泛函模型分解为三个子问题交替迭代求解，得到低秩近似后的矩阵L，所述三个子问题分别为：

第一子问题：固定Z和Y，求解L；

第二子问题：固定L和Y，求解Z；

第三子问题：固定L和Z，求解Y。

5.如权利要求4所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，所述三个子问题的对应迭代目标函数为：

其中，S()为软阈值收缩算子，其公式为D()为奇异值阈值收缩算子，其公式为D_b(a)＝US_b(Σ)V^T,a＝UΣV^T。

6.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，步骤(1)所述边长为8像素。

7.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特 征在于，步骤(1)所述步长为6像素。

8.如权利要求1所述的合成孔径雷达图像相干噪声的抑制方法，其特征在于，步骤(1)所述搜索框半径为21像素。