CN103077507A - 基于Beta算法的多尺度SAR图像降噪方法 - Google Patents

Abstract

本发明公开了一种基于Beta算法的多尺度SAR图像降噪方法，主要解决传统字典学习方法不适用SAR图像降噪的问题。其降噪过程为:对待降噪SAR图像进行像素分类，得到边缘图像A、纹理图像B和同质图像C；对这些图像分别进行重叠块提取和中心化操作，得到它们各自的训练样本集；初始化字典为DCT字典；用字典进行稀疏编码得稀疏系数矩阵；用稀疏系数矩阵更新字典的第k列；重复稀疏编码步骤和字典更新步骤共K次，得最终字典和最终稀疏系数矩阵，并将其相乘得降噪后的边缘图像A’、纹理图像B’和同质图像C’，进而得降噪后的SAR图像为：I’=A’+B’+C’。本发明具有在有效去除SAR图像中噪声的同时能保留图像的纹理和边缘信息的优点，可用于SAR图像目标识别。

Description

基于Beta算法的多尺度SAR图像降噪方法

技术领域

本发明属于图像处理领域，特别是一种涉及SAR图像降噪的方法，可应用于目标识别。

背景技术

合成孔径雷达SAR是一种高分辨率雷达体制，可应用于军事，农业，导航，地理监视等诸多领域。它与其它遥感成像系统，光学成像系统相比有很多差异。在军事目标识别方面，SAR图像被广泛的应用在目标检测领域，而SAR图像降噪则是从图像处理到图像分析的重要步骤，是目标分类和识别的基础。本质上SAR图像反映的是目标的电磁散射特性和结构特性，但由于SAR发射的是相干电磁波，当它向地面发射电磁波时，其后向散射回来的总回波并不完全由地物目标的散射系数决定，而是随着这些散射系数有很大的随机起伏，这种起伏的存在使SAR图像中出现许多斑点，这些斑点与较小的地物目标掺杂在一起严重影响图像的质量，这些斑点称为相干斑噪声。SAR成像的特殊性使得针对该类图像的降噪方法有别于普通光学图像：非光学成像方式导致SAR图像灰度级变化缓慢，因此光学图像常用的基于高斯白噪声假设的降噪方法对此类图像收效甚微。

针对SAR图像的以上特点，经典降噪方法有：

基于多视处理技术的降噪方法。它在SAR图像成像阶段对同一个场景获得L个独立的测量，并将这些测量进行平均，该技术虽然简单，但却是以牺牲图像分辨率为代价。因此，出现了各种对成像后的SAR图像进行滤波的降噪算法。

基于空域的降噪方法。主要包括Lee滤波、Frost滤波和Gamma-MAP滤波等，它们都假设真实场景是平稳的并使用乘性的相干斑模型,Lee滤波直接应用线性最小均方误差MMSE准则，得到一个一般的降噪表达式，通过计算这个表达式得到降噪结果；Frost滤波也使用了MMSE准则，但其发现场景的自相关函数与像素之间的空间距离呈负指数关系，基于此构建了一个自相关函数，将图像与该自相关函数进行卷积获得最终滤波结果。由于这两种滤波器并未涉及相干斑或场景的具体统计特性，因此Lopes A.等人假设真实场景的统计特性符合Gamma分布，相干斑的统计特性符合高斯分布，结合最大后验概率MAP准则求得抑斑结果，这就是Gamma-MAP滤波器，参见Lopes A,Nezry E,Touzi R,Laur H.Maximum a posteriorispeckle filtering and first order texture models in SAR images[C].Geoscience andRemote Sensing Symposium,1990 IGARSS′90′Remote Sensing Science for the Nineties′,10th Annual International.1990.2409-12.。上述的三种滤波方法虽然抑斑效果较好，但是它们都难以保持图像的细节特征，使图像边缘和线性目标模糊。原因在于它们都假设真实的SAR图像场景是平稳的，而这个假设在边缘和细节区域不成立。因此为了更好地满足“平稳”这个假设，科学家们做了以下改进：一种是将SAR图像的真实场景分为均匀和非均匀的，分别采用不同抑斑方法。例如Lopes A.等人提出的增强系列滤波器。另一种是通过统计对图像像素进行分类。例如Lee J.S.的增强Lee滤波器，它通过对图像的像素值进行统计并划定一个概率区间，将区间之外的像素作为噪声，经过滤波处理得到干净的图像，参见Lee JS,Wen JH,Ainsworth TL,Chen KS,Chen AJ.Improved sigma filter for speckle filtering of SAR imagery[J].Geoscience and Remote Sensing,IEEE Transactions on.2009,Vol.47(1).202-13.。

基于变换域的降噪方法。主要有小波变换、Bandelet变换、Curvelet变换和非下采样Contourlet变换等。这些变换域滤波方法相比于经典的空域滤波方法来说，图像的线性目标能得到较好的保持，但是由于噪声和图像边缘都是高频信号，因此抑斑后的图像在均匀区域和边缘附近常会出现伪吉布斯效应。随着“稀疏表示”热潮的席卷，变换域方法可以看成是将图像投影在一种或者几种基上，实现在用较少的系数表示图像的同时达到降噪的目的。但是由于这些基都有特定的方向，而一幅真实图像存在各种方向，因此变换域的稀疏表示方法不能很好地保持图像中的边缘信息。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于Beta算法的SAR图像多尺度降噪方法，在有效去除图像中的斑点噪声的同时保留图像的纹理和边缘信息，且在图像的均匀区域和边缘附近避免出现伪吉布斯效应。

实现本发明目的的技术方案是：使用像素分类将一幅SAR图像分为三幅图像，分别为边缘图像，纹理图像和同质图像。对这三幅图像分别使用Beta算法得到它们各自的稀疏表示系数和字典，将它们各自的稀疏表示系数和字典相乘得三幅降噪后图像，将这三幅降噪后图像相加即为原SAR图像的降噪结果图。其具体实现步骤如下：

(1)对待降噪SAR图像I进行像素分类，得到边缘图像A、纹理图像B和同质图像C；

(2)对边缘图像A进行重叠块提取和中心化操作，得到训练样本集

其中是长度为n=b₁×b₁的列向量，b₁为重叠块提取窗的边长，M=(N-b₁+1)²，N为边缘图像A的行数；

(3)将字典D初始化为离散正弦变换DCT字典，字典D大小为n×K，K=R×n，R为冗余度；

(4)用字典D对训练样本集

进行稀疏编码得到大小为K×M的稀疏表示系数矩阵

(5)用稀疏表示系数矩阵

更新字典D的第k列d_k；

(6)重复步骤(4)和步骤(5)共K次，得到最终字典

和最终稀疏表示系数矩阵

将最终字典

和最终稀疏表示系数矩阵

相乘，得降噪后边缘图像A′；

(7)对纹理图像B使用大小为b₂×b₂的重叠块提取窗，b₂>b₁，重复步骤(2)至步骤(6)得到降噪后纹理图像B′；

(8)对同质图像C使用大小为b₃×b₃的重叠块提取窗，b₃>b₂>b₁，重复步骤(2)至步骤(6)得到降噪后同质图像C′；

(9)将降噪后边缘图像A′、降噪后纹理图像B′和降噪后同质图像C′相加，得到降噪后的SAR图像I′。

本发明与现有的技术相比具有以下优点：

1、本发明由于使用像素分类和多尺度方法，能够很好保留SAR图像边缘和纹理区域的细节信息，保证了图像信息的完整性；

2、本发明使用Beta字典学习方法，很好地去除SAR图像中的噪声；

3、本发明是在空域上对SAR图像进行处理，避免了在图像的均匀区域和边缘附近出现伪吉布斯效应；

仿真结果表明，本发明方法较增强Lee滤波和Gamma-MAP滤波能更有效的进行SAR图像降噪。

附图说明

图1是本发明的流程图；

图2是本发明在一幅4视强度SAR图像上的仿真结果图；

图3是本发明在一幅6视强度SAR图像上的仿真结果图；

图4是本发明在一幅2视幅度SAR图像上的仿真结果图。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤一、对待降噪SAR图像I进行像素分类，得到边缘图像A、纹理图像B和同质图像C。

(1a)输入待降噪SAR图像I，根据公式1)计算该待降噪SAR图像I的局部方差图variance，并绘制局部方差图variance直方图，

$\mathrm{variance}(i,j)=\mathrm{\Σ}{(\frac{y[i,j]}{\stackrel{\‾}{y}[m,n]}-1)}^2---1)$

其中，variance(i,j)为局部方差图variance中的元素，m,n为窗口的大小，设置m=n=5，

是此窗口内所有元素的平均值，y[i,j]表示待降噪SAR图像I中的元素；

(1b)设置待降噪SAR图像I的标准方差为小方差阈值γ，根据局部方差图variance直方图，得到待降噪SAR图像I的大方差阈值V，V>γ，根据公式2)得到方差标记label1(i,j)：

$\mathrm{label}1(i,j)=\left\{\begin{array}{cc}255& \mathrm{variance}(i,j)\&GreaterEqual;V\\ 155& V>\mathrm{variance}(i,j)\&GreaterEqual;\mathrm{\&gamma;}\\ 55& \mathrm{variance}(i,j)<\mathrm{\&gamma;}\end{array}\right.---2)$

(1c)用3×3的窗口求该待降噪SAR图像I的均值图Mean1，用9×9的窗口求该待降噪SAR图像I的均值图Mean2，根据公式3)计算差异图difference并绘制差异图difference直方图，其中difference(i,j)为差异图difference中的元素为：

difference(i,j)=abs(Mean1(i,j)-Mean2(i,j))    3)

式中abs表示取绝对值，Mean1(i,j)表示均值图Mean1中的元素，Mean2(i,j)表示均值图Mean2中的元素；

(1d)根据差异图difference直方图确定两个差异阈值M₂和M₁，M₂>M₁，根据公式4)得到差异标记label2(i,j)：

$\mathrm{label}2(i,j)=\left\{\begin{array}{cc}255& \mathrm{difference}(i,j)\&GreaterEqual;M_2\\ 155& M_2>\mathrm{difference}(i,j)\&GreaterEqual;M_1\\ 55& \mathrm{difference}(i,j)<M_1\end{array}\right.---4)$

(1e)按公式5)将步骤(1b)和步骤(1d)的结果合并，得到边缘图像A、纹理图像B和同质图像C，即：

其中Max表示取方差标记label1(i,j)和差异标记label2(i,j)的最大值。

步骤二、对边缘图像A进行重叠块提取和中心化操作，得到训练样本集

其中

是长度为n=b₁×b₁的列向量，b₁为重叠块提取窗的边长，M=(N-b₁+1)²，N为边缘图像A的行数。

取滑动因子s=1，用大小为b₁×b₁的窗口对边缘图像A进行重叠块提取，得到重叠块向量集

其中yi是长度为n=b₁×b₁的列向量，M=(N-b₁+1)²，N为边缘图像A的行数，将Y减去其本身的均值得到训练样本集

是长度为n=b₁×b₁的列向量。

步骤三、将字典D初始化为离散正弦变换DCT字典，字典D大小为n×K，K=R×n，R为冗余度。

步骤四、用字典D对训练样本集

进行稀疏编码，得到大小为K×M的稀疏表示系数矩阵

(4a)初始化稀疏促进矩阵Z为大小等于K×K的全1矩阵；

(4b)初始化稀疏促进矩阵Z的第i列向量z_i的每个元素等于1的概率π为K×1大小的全0.5矩阵，根据公式6)更新概率π的第k个成分π_k，得到新的概率

的第k个成分

${\widehat{\mathrm{\&pi;}}}_k=\mathrm{Beta}(\frac{1}{K}+\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}},\frac{K-1}{K}+M-\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}})---6)$

其中Beta表示Beta概率密度函数，z_ik表示Z的第i列z_i的第k个成分；

(4c)根据公式7)更新稀疏促进矩阵Z的第i列的第k个成分z_ik，得到新的稀疏促进矩阵

的第i列的第k个成分

${\hat{z}}_{\mathrm{ik}}=\left\{\begin{array}{cc}1& {\widehat{\mathrm{\&pi;}}}_k={\mathrm{\&pi;}}_k\\ 0& {\widehat{\mathrm{\&pi;}}}_k\&NotEqual;{\mathrm{\&pi;}}_k\end{array}\right.---7)$

(4d)初始化稀疏表示系数矩阵α为大小等于K×M的全0矩阵，根据公式8)更新稀疏表示系数矩阵α，得到新的稀疏表示系数矩阵

$\widehat{\mathrm{\&alpha;}}=\arg \underset{\mathrm{\&alpha;}}{\min }{\left|\right|\hat{Z}\&CenterDot;\mathrm{\&alpha;}\left|\right|}_0+\mathrm{\&lambda;}\&CenterDot;\frac{{\left|\right|{\hat{Y}}^{^}-D\&CenterDot;Z\&CenterDot;\mathrm{\&alpha;}\left|\right|}_2^2}{{\left|\right|D\&CenterDot;\hat{Z}\&CenterDot;\mathrm{\&alpha;}\left|\right|}_2^2}---8)$

其中

为稀疏表示系数矩阵α取最小值，‖·‖₀为0范数，‖·‖₂为2范数，λ为调整因子取值为2。

步骤五、用稀疏表示系数矩阵

更新字典D的第k列向量d_k。

步骤六、重复步骤四和步骤五共K次，得到最终字典

和最终稀疏表示系数矩阵

将最终字典

和最终稀疏表示系数矩阵

相乘，得到降噪后边缘图像A′。

步骤七、对纹理图像B使用大小为b₂×b₂的重叠块提取窗，b₂>b₁，重复步骤二至步骤六，得到降噪后纹理图像B′。

步骤八、对同质图像C使用大小为b₃×b₃的重叠块提取窗，b₃>b₂>b₁，重复步骤二至步骤六，得到降噪后同质图像C′。

步骤九、将降噪后边缘图像A′、降噪后纹理图像B′和降噪后同质图像C′相加，得到降噪后的SAR图像I′。

本发明的的效果可通过以下仿真进一步说明：

1 仿真内容：应用本发明方法和增强Lee滤波方法、Gamma-MAP滤波方法分别对三幅SAR图像进行降噪实验，并从同质区域平滑程度、比值图像的均值和方差、边缘保持程度评价这些方法的各自性能。

2 仿真实验结果

仿真1 用本发明以及增强Lee滤波方法、Gamma-MAP滤波方法对一幅4视强度SAR图像进行降噪仿真，其结果如图2所示。其中：

图2(a)为原始4视强度SAR图像，在该SAR图像中标记A、B两个同质区域；

图2(b)为用增强Lee滤波方法对图2(a)进行降噪得到的结果；

图2(c)为用Gamma-MAP滤波方法对图2(a)进行降噪得到的结果；

图2(d)为用本发明对图2(a)进行降噪得到的结果。

图2的客观评价指标如表1，表1中计算了同质区域A和B的等效视数，等效视数越大说明同质区域降噪效果越好。

由表1可见，本发明的等效视数最大，表示本发明对图像的同质区域降噪效果最好。本发明的比值图均值和方差及边缘保持指数都更接近理想值，表明本发明能很好保留图像的边缘和纹理信息。

表1 图2的客观评价指标

仿真2 用本发明以及增强Lee滤波方法、Gamma-MAP滤波方法对一幅6视强度SAR图像进行降噪仿真，仿真结果如图3所示。其中：

图3(a)为原始6视强度SAR图像，在该SAR图像中标记同质区域A；

图3(b)为用增强Lee滤波方法对图3(a)进行降噪得到的结果；

图3(c)为用Gamma-MAP滤波方法对图3(a)进行降噪得到的结果；

图3(d)为用本发明对图3(a)进行降噪得到的结果。

图3的客观评价指标如表2，表2计算了同质区域A的等效视数，等效视数越大说明同质区域降噪效果越好。

表2 图3的客观评价指标

由表2可见，本发明的等效视数最大，且比值图的均值和方差及边缘保持指数都更接近理想值，说明本发明不仅能有效去除噪声且能很好保留图像的边缘和纹理信息，从图3(d)看出使用本发明对图像降噪避免了伪吉布斯效应。通过仿真1和仿真2表明本发明对不同视数的SAR图像都有很好的去噪效果。

仿真3 用本发明以及增强Lee滤波方法、Gamma-MAP滤波方法对2视幅度SAR图像进行降噪，其效果如图4所示。其中：

图4(a)为原始2视幅度SAR图像，在该SAR图像中标记同质区域A、B、C；

图4(b)为用增强Lee滤波方法对图4(a)进行降噪得到的结果；

图4(c)为用Gamma-MAP滤波方法对图4(a)进行降噪得到的结果；

图4(d)为用本发明对图4(a)进行降噪得到的结果。

图4的客观评价指标如表3，表3计算了同质区域A、B和C的等效视数，等效视数越大说明同质区域降噪效果越好。

表3 图4的客观评价指标

由表3可见，本发明无论是降噪能力还是边缘和纹理信息保留能力都优于增强Lee滤波方法和Gamma-MAP滤波方法，通过仿真1、仿真2和仿真3表明本发明对强度格式的SAR图像和幅度格式的SAR图像都有很好的去噪效果。

Claims (4)

1.一种基于Beta算法的多尺度SAR图像降噪方法，包括如下步骤：

1)对待降噪SAR图像I进行像素分类，得到边缘图像A、纹理图像B和同质图像C；

2)对边缘图像A进行重叠块提取和中心化操作，得到训练样本集

其中

是长度为n=b₁×b₁的列向量，b₁为重叠块提取窗的边长，M=(N-b₁+1)²，N为边缘图像A的行数；

3)将字典D初始化为离散正弦变换DCT字典，字典D大小为n×K，K=R×n，R为冗余度；

4)用字典D对训练样本集

进行稀疏编码，得到大小为K×M的稀疏表示系数矩阵

5)用稀疏表示系数矩阵

更新字典D的第k列d_k；

6)重复步骤4)和步骤5)共K次，得到最终字典

和最终稀疏表示系数矩阵

将最终字典

和最终稀疏表示系数矩阵

相乘，得降噪后边缘图像A′；

7)对纹理图像B使用大小为b₂×b₂的重叠块提取窗，b₂>b₁，重复步骤2)至步骤6)得到降噪后纹理图像B′；

8)对同质图像C使用大小为b₃×b₃的重叠块提取窗，b₃>b₂>b₁，重复步骤2)至步骤6)得到降噪后同质图像C′；

9)将降噪后边缘图像A′、降噪后纹理图像B′和降噪后同质图像C′相加，得到降噪后的SAR图像I′。

2.根据权利要求1所述的基于Beta算法的多尺度SAR图像降噪方法，其中步骤1)按如下步骤进行：

1a)根据公式(1)计算待降噪SAR图像I的局部方差图variance，并绘制局部方差图variance直方图，

$\mathrm{variance}(i,j)=\mathrm{\&Sigma;}{(\frac{y[i,j]}{\stackrel{\&OverBar;}{y}[m,n]}-1)}^2---\left(1\right)$

其中，variance(i,j)为局部方差图variance中的元素，m,n为窗口的大小，设置m=n=5，

是此窗口内所有元素的平均值，y[i,j]表示待降噪SAR图像I中的元素；

1b)设置待降噪SAR图像I的标准方差为小方差阈值γ，根据局部方差图variance直方图，得到待降噪SAR图像I的大方差阈值V，V>γ，根据公式(2)得到方差标记label1(i,j)：

$\mathrm{label}1(i,j)=\left\{\begin{array}{cc}255& \mathrm{variance}(i,j)\&GreaterEqual;V\\ 155& V>\mathrm{variance}(i,j)\&GreaterEqual;\mathrm{\&gamma;}\\ 55& \mathrm{variance}(i,j)<\mathrm{\&gamma;}\end{array}\right.---\left(2\right)$

1c)用3×3的窗口求该待降噪SAR图像I的均值图Mean1，用9×9的窗口求该待降噪SAR图像I的均值图Mean2，根据公式(3)计算差异图difference并绘制差异图difference直方图，其中difference(i,j)为差异图difference中的元素为：

difference(i,j)=abs(Mean1(i,j)-Mean2(i,j))    (3)

式中abs表示取绝对值，Mean1(i,j)表示均值图Mean1中的元素，Mean2(i,j)表示均值图Mean2中的元素；

1d)根据差异图difference直方图确定两个差异阈值M₂>M₁，根据公式(4)得差异标记label2(i,j)：

$\mathrm{label}2(i,j)=\left\{\begin{array}{cc}255& \mathrm{difference}(i,j)\&GreaterEqual;M_2\\ 155& M_2>\mathrm{difference}(i,j)\&GreaterEqual;M_1\\ 55& \mathrm{difference}(i,j)<M_1\end{array}\right.---\left(4\right)$

1e)按公式(5)将步骤1b)和步骤1d)的结果合并，得到边缘图像A、纹理图像B和同质图像C，

其中Max表示取方差标记label1(i,j)和差异标记label2(i,j)的最大值。

3.根据权利要求1所述的基于Beta算法的多尺度SAR图像降噪方法，其中步骤2)所述的对边缘图像A进行重叠块提取和中心化操作，按如下步骤进行：

取滑动因子s=1，用大小为b₁×b₁的窗口对边缘图像A进行重叠块提取，得到重叠块向量集

其中yi是长度为n=b₁×b₁的列向量，M=(N-b₁+1)²，N为边缘图像A的行数，将Y减去其本身的均值得到训练样本集

是长度为n=b₁×b₁的列向量。

4.根据权利要求1所述的基于Beta算法的多尺度SAR图像降噪方法，其中步骤4)所述的用字典D对训练样本集

进行稀疏编码，按如下步骤进行：

4a)初始化稀疏促进矩阵Z为大小等于K×K的全1矩阵；

4b)初始化稀疏促进矩阵Z的第i列向量z_i的每个元素等于1的概率π为K×1大小的全0.5矩阵，根据公式(6)更新概率π的第k个成分π_k，得到新的概率

的第k个成分

${\widehat{\mathrm{\&pi;}}}_k=\mathrm{Beta}(\frac{1}{K}+\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}},\frac{K-1}{K}+M-\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}})---\left(6\right)$

其中Beta表示Beta概率密度函数，z_ik表示Z的第i列z_i的第k个成分；

4c)根据公式(7)更新稀疏促进矩阵Z的第i列的第k个成分z_ik，得到新的稀疏促进矩阵

的第i列的第k个成分

${\hat{z}}_{\mathrm{ik}}=\left\{\begin{array}{cc}1& {\widehat{\mathrm{\&pi;}}}_k={\mathrm{\&pi;}}_k\\ 0& {\widehat{\mathrm{\&pi;}}}_k\&NotEqual;{\mathrm{\&pi;}}_k\end{array}\right.---\left(7\right)$

4d)初始化稀疏表示系数矩阵α为大小等于K×M的全0矩阵，根据公式(8)更新稀疏表示系数矩阵α，得到新的稀疏表示系数矩阵

$\widehat{\mathrm{\&alpha;}}=\arg \underset{\mathrm{\&alpha;}}{\min }{\left|\right|\hat{Z}\&CenterDot;\mathrm{\&alpha;}\left|\right|}_0+\mathrm{\&lambda;}\&CenterDot;\frac{{\left|\right|{\hat{Y}}^{^}-D\&CenterDot;Z\&CenterDot;\mathrm{\&alpha;}\left|\right|}_2^2}{{\left|\right|D\&CenterDot;\hat{Z}\&CenterDot;\mathrm{\&alpha;}\left|\right|}_2^2}---\left(8\right)$

其中

为稀疏表示系数矩阵α取最小值，‖·‖₀为0范数，‖·‖₂为2范数，λ为调整因子取值为2。

Images