CN101639537B - 基于方向波域混合高斯模型的sar图像噪声抑制方法 - Google Patents

Abstract

本发明公开一种基于方向波域混合高斯模型的SAR图像噪声抑制方法，涉及图像处理领域，主要克服现有去噪方法对SAR图像造成的模糊失真和改进方向波的陪集分解添加零值过多的问题。其实现步骤为：(1)对SAR图像进行卡通纹理模型分解；(2)纹理部分进行改进的方向波变换；(3)对陪集的高频子带系数进行混合高斯模型的拟合；(4)估计噪声方差；(5)利用贝叶斯最小均方误差估计无噪系数；(6)无噪系数进行改进的方向波逆变换，得到去噪后的纹理部分；(7)将卡通部分和去噪后的纹理部分相加，得到去噪后的SAR图像。本发明具有改进了去噪后SAR图像的清晰度、细节信息保持较好、均匀区域较平滑的优点，可用于滤除SAR图像中的相干斑噪声。

Description

基于方向波域混合高斯模型的SAR图像噪声抑制方法

技术领域

本发明属于图像处理领域，具体地涉及SAR图像噪声抑制方法，该方法可用于滤除SAR图像中的相干斑噪声。

背景技术

合成孔径雷达(Synthetic Aperture Radar，简称SAR)是一种能在各种天气条件下获得类似光学图像的高分辨率图像的全天候、全天时雷达。它的相干成像机理不可避免地会产生相干斑噪声，使其图像不能有效反映地物目标的散射特性。良好的斑点抑制方法要能够在减少噪声的同时尽可能保持SAR图像中的细节信息，去噪后的SAR图像在视觉效果和信息量方面都更有利于后续的计算机自动分析和理解。

近年来，SAR图像相干斑噪声抑制技术飞速发展，可分为成像前多视平滑技术和成像后滤波技术两大类，成像后处理技术中空域滤波和变换域方法得到了广泛地研究。对于空域滤波方法，如Lee，Gamma-MAP等基于乘性斑点噪声模型的自适应局部统计滤波方法得到了应用，这些滤波器算法的问题在于滤波窗口大小的选择。均匀的区域需要较大的滤波器窗口来增强斑点抑制效果，但大的滤波器窗口同时会降低图像的实际分辨率，造成图像中边缘和线状目标的模糊，很在均匀区域的平滑和细节信息的保持之间兼顾。变换域方法如小波方法通过对其系数进行阈值分析处理或者是建立系数统计模型实现对SAR图像噪声的抑制，虽然能较好保持边缘细节，但平滑能力有限。为了更好地挖掘如何更有效地表示SAR图像，多尺度几何分析是很好的思路，并且在这方面开展了许多富有成效的工作。孙强等人在文章：统计先验指导的非下采样Contourlet变换域SAR图像降斑中提出的降班方法取得了一定得成果，在边缘保持和后向散射系数保持的性能上优于以上几种方法。FlorianLuisier等人在文章：SURE-LET Multichannel Image Denoising：Interscale OrthonormalWavelet Thresholding中提出的去噪方法，在细节的保持上取得了较好的效果。但是现有的这些方法在细节的保持和去噪程度上不能达到比较好的平衡，为了达到更好的去噪效果的同时尽量地保持图像的细节，需要提出新的方法。

2005年由Vladan Velisavljevic等人提出的Directionlet变换是一种新的多尺度几何分析方法，它采用整数格变换，可分离滤波和临界采样结构设计，计算复杂度小，滤波器构造相对简单，能够提供SAR图像更好的稀疏表示，而且能够有效地捕捉SAR图像中的各向异性特征。但是采用的陪集分解操作由于加入的零值太多，对分解后得到的系数进行模型的拟合带来了困难。

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种基于Directionlet域混合高斯模型的SAR图像噪声抑制方法，使得图像去噪和图像细节的保持取得较好的平衡，以及对Directionlet分解系数进行较好的模型拟合，提高对SAR图像的噪声抑制的效果。

实现本发明技术目的技术方案，包括如下步骤：

(1)对SAR图像进行卡通纹理模型分解，选择最优的迭代次数，得到SAR图像的卡通部分和纹理部分；

(2)保持SAR图像卡通部分不变，将纹理部分进行改进的Directionlet变换；

(2a)选定SAR图像纹理部分的变换方向和队列方向，其斜率构成采样矩阵M_Λ：

$M_{\mathrm{\Λ}}=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]$ $a_1,a_2,b_1,b_2\⋐Z$

其中，沿向量d₁的方向称为纹理部分的变换方向，其斜率表示为r₁＝b₁/a₁，沿向量d₂的方向称为纹理部分的队列方向，其斜率表示为r₂＝b₂/a₂；

(2b)对纹理部分进行M_Λ采样，得到|det(M_Λ)|个陪集，再对得到的陪集分别进行同样的二次陪集分解，得到2|det(M_Λ)|个子陪集，将零点去除；

(2c)对每个子陪集沿变换和队列方向分别进行两次和一次的一维正交小波变换，得到相应的高频和低频系数子带，在低频子带的变换和队列方向再分别进行两次和一次一维正交小波变换，获得14个高频子带和1个低频子带系数为W_Y(i，j)，( $i,j\⋐Z$ 表示每个系数的坐标)；

(3)对变换后得到的各陪集的高频子带系数进行混合高斯模型的拟合，并采用EM算法估计混合模型中的参数；

(4)估计各陪集的高频子带系数的噪声方差

(5)利用Bayes最小均方误差估计法，计算各陪集的高频子带的无噪系数

(6)对无噪系数进行改进的Directionlet逆变换，得到去噪后的纹理部分；

(7)将去噪后的纹理部分与保持不变的卡通部分相加，得到去噪后的SAR图像。

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

1、本发明采用了卡通纹理模型对SAR图像进行分解，因此能够在有效地去除噪声的同时也很好地保持了图像细节等特征。大量的试验证实，在SAR图像去噪处理中，若直接将SAR图像进行去噪的处理，很难兼顾均匀区域的平滑和细节信息的保持，因此需要行之有效的预处理工作来尽可能的达到两者之间较好的平衡。

在SAR图像去噪中使用卡通纹理模型的优点在于：将卡通和纹理两部分分开后，首先SAR图像的一些细节信息可以保留在卡通部分，而不会像其它一些去噪方法被过于损失；其次对剩余纹理部分的噪声也可以有针对性的去除；在仿真实验中也发现，本发明提出首先对SAR图像进行卡通纹理模型的分解，在此基础上进行去噪处理得到的SAR图像细节保留更多，图像也更加清晰，同时也有效地提高了一些衡量去噪效果的指标。

2、本发明改进的Direcctionlet变换方案对后续的模型拟合更有利。

未改进的Direcctionlet变换在进行陪集分解时添加的零值的数目太多，如果直接将第一次分解得到的陪集进行后续的小波变换以及模型拟合，则会导致图像曲线在零点处峰值过尖，无法进行拟合。为了排除这种干扰，必须将零点去除，本发明采用了二次陪集分解，即在得到的一次分解陪集的基础上再进行一次陪集分解，然后找出二次陪集分解以后对应于原图像中的像素，并提取出来，避免了过多的零值干扰，这样就可以较好地进行模型拟合。

仿真实验结果表明：本发明应用于几幅SAR图像的去噪，能够得到更高的衡量去噪效果的指标和更好的视觉效果。

附图说明

图1是本发明的主要操作过程示意图；

图2是用本发明与已有方法对测试图像x波段等效视数是2的幅度SAR图像256×256的去噪效果对比图；

图3是用本发明与已有方法对测试图像ku波段等效视数是4的强度SAR图像256×256的去噪效果对比图。

具体实施方式

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

步骤1：对SAR图像进行卡通纹理模型分解，选择最优的迭代次数，得到SAR图像的卡通部分和纹理部分，该卡通纹理模型的数学表达式为：

f(x，y)＝u(x，y)+v(x，y)

其中u∈BV(R²)，由分片光滑的区域及清晰的边缘构成，被视为SAR图像的卡通部分，v就是剩余的纹理和噪声，称之为纹理部分。当v∈G，G是一个定义的Banach空间时，得到模型：

$\underset{u\∈\mathrm{BV}}{\inf }\{E\left(u\right)=\∫|\▿u|+\mathrm{\λ}{\left|\right|v\left|\right|}_u,f=u+v\}$

其中E(u)表示u的能量，通过在保证u∈BV(R²)，其能量泛函最小化时达到对图像f的卡通纹理分解。经过多次试验证明，一般最优的迭代次数为50～60次。

步骤2：保持SAR图像卡通部分不变，将纹理部分进行改进的Directionlet变换；选定SAR图像纹理部分的变换方向和队列方向，其斜率构成采样矩阵M_Λ：

$M_{\mathrm{\Λ}}=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]$ $a_1,a_2,b_1,b_2\⋐Z$

其中，沿向量d₁的方向称为纹理部分的变换方向，其斜率表示为r₁＝b₁/a₁，沿向量d₂的方向称为纹理部分的队列方向，其斜率表示为r₂＝b₂/a₂；对纹理部分进行M_Λ采样，得到|det(M_Λ)|个陪集，再对得到的陪集再分别进行同样的二次陪集分解，得到2|det(M_Λ)|个子陪集，将零点去除；对每个子陪集沿变换和队列方向分别进行两次和一次的一维正交小波变换，得到相应的高频和低频系数子带，在低频子带的变换和队列方向再分别进行两次和一次一维正交小波变换，获得14个高频子带和1个低频子带系数为W_Y(i，j)， $i,j\⋐Z$ 表示每个系数的坐标；

步骤3：对变换后得到的各陪集的高频子带系数进行混合高斯模型的拟合，采用EM算法估计混合模型中的参数。

假设W_X表示含有噪声的Directionlet系数，则其混合概率密度函数可表示为：

$p_{W_X}\left(W_X\right)=\underset{k=\mathrm{0,1}}{\mathrm{\Σ}}p(S=k)p\left(W_X\right|S=k)$

$\underset{k=\mathrm{0,1}}{\mathrm{\Σ}}p(S=k)=1$

其中p(S＝k)为混合参数，满足上式所示的关系；

$p\left(W_X\right|S=k)=\frac{1}{\sqrt{2\mathrm{\π}{d}_k}}\exp (-\frac{x^2}{2d_k})$

其中p(W_X|S＝k)为第k个高斯分量的概率密度函数，均值为1，方差为d_k。

步骤4：估计各陪集的高频子带系数的噪声方差

这里采用一种启发式的先验估计形式来估计噪声标准差：

${\mathrm{\σ}}_{W_B}=\mathrm{median}\left(\right|W_X|-\mathrm{median}\left(W_X\right))/0.6745$

其中，W_X表示含噪信号多尺度分解的高频子带系数。

步骤5：利用Bayes最小均方误差估计计算各陪集的高频子带的无噪系数

具体求解公式如下：

对于分解后的高频系数W_X(i，j)，无噪Directionlet系数W_Y的贝叶斯估计表达式为：

${\hat{W}}_Y=\mathrm{\η}{W}_X$

在Directionlet系数的概率密度函数为混合高斯函数的前提下，基于最小均方误差估计的系数η表示为：

$\mathrm{\η}=\underset{k=\mathrm{0,1}}{\mathrm{\Σ}}p(S=k|W_X)\frac{{\mathrm{\σ}}_{W_X}^2-{\mathrm{\σ}}_{W_B}^2}{{\mathrm{\σ}}_{W_X}^2}$

则结合上述两式，基于Bayes最小均方误差估计的无噪Directionlet系数为：

${\hat{W}}_Y=\underset{k=\mathrm{0,1}}{\mathrm{\Σ}}p(S=k|W_X)\frac{{\mathrm{\σ}}_{W_X}^2-{\mathrm{\σ}}_{W_B}^2}{{\mathrm{\σ}}_{W_X}^2}{W}_X$

这里

是相对于隐状态S的Directionlet噪声系数W_X的方差，

是估计的噪声系数的方差。其中：

$p(S=k|W_X)=\frac{p\left(W_X\right|S=k)p(S=k)}{p\left(W_X\right)}.$

步骤6：对滤波后的系数进行改进的Directionlet逆变换，得到去噪后的纹理部分。

步骤7：将保持不变的卡通部分和去噪后的纹理部分相加从而得到去噪后的SAR图像。

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

1、仿真条件

采用了SAR图像去噪中常用的x波段等效视数是2的幅度田野SAR图像，此图像是英国Bedfordshire地区，分辨率为3m，以及美国新墨西哥州Albuquerque地区，分辨率为1m，图像大小为256×256、ku波段等效视数是4的强度SAR图像，并且选取了具有代表性的5块区域计算衡量去噪效果的指标。分别用Lee滤波，Gamma-MAP滤波，小波软阈值，平稳小波软阈值，非下采样轮廓波阈值法，本发明去噪方法进行去噪处理，其中Lee和Gamma-MAP滤波器均采用3×3大小的滤波窗口。本发明方法对256×256大小的SAR图像使用了2层改进的Directionlet分解，同时选取2个变换方向，分别为45°和-45°方向，对应的采样矩阵为：

$M_1=\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right],$ $M_2=\left[\begin{array}{cc}-1& 1\\ 1& 1\end{array}\right]$

2、仿真结果分析

对田野SAR图像的仿真的结果如图2，其中：图2(a)为田野SAR图像，图2(b)为使用增强Lee滤波方法去噪后的结果，图2(c)为使用Gamma-MAP滤波方法去噪后的结果，图2(d)为使用小波软阈值方法去噪后的结果，图2(e)为使用平稳小波软阈值方法去噪后的结果，图2(f)为使用非下采样轮廓波阈值法方法去噪后的结果，图2(g)为使用本发明方法去噪后的结果。

对Albuquerque地区SAR图像的仿真结果如图3，其中：图3(a)为Albuquerque地区SAR图像，图3(b)为使用增强Lee滤波方法去噪后的结果，图3(c)为使用Gamma-MAP滤波方法去噪后的结果，图3(d)为使用小波软阈值方法去噪后的结果，图3(e)为使用平稳小波软阈值方法去噪后的结果，图3(f)为使用非下采样轮廓波阈值法方法去噪后的结果，图3(g)为使用本发明方法去噪后的结果。

从图2中可以看出，本发明去噪后不仅有较好的视觉效果，而且在边缘和细节信息的保持方面也优于其它几种去噪方法，有效地滤除了斑点噪声，提高了整体上的清晰度，在图像质量和视觉效果方面相比于其他几种方法均有所改善。如在图2(g)中经本发明去噪后SAR图像边缘和左下角的三个点目标比较清晰，均匀区域也得到了较好的平滑，而在图2(b)中Lee滤波和图2(c)中Gamma-MAP滤波去噪后图像的点目标模糊，图像细节信息丢失比较严重。在图2(d)中小波软阈值法去噪后的SAR图像在均匀区域上吉布斯现象明显，同时存在边缘模糊的现象。在图2(e)中平稳小波软阈值去除噪声的方法对SAR图像的均匀区域进行了比较好的处理，但是图像整体存在模糊失真的现象。在图2(f)中基于非下采样轮廓波阈值方法处理后的SAR图像，其均匀区域滤波不完全，划痕现象比较明显，边缘部分存在一定程度的失真。

从图3中可以看出，本发明方法在视觉效果的提高及图像质量的改善两方面均优于对比的方法，再次验证了其有效性。

本发明与已有方法对田野SAR图像的去噪性能指标对比表如表1所示：

表1各方法对田野SAR图像去噪性能对比

从表1中可以看到，本发明方法在区域1～3的的等效视数ENL的计算上取得了高于表中列举方法的值，表明本发明方法相干斑抑制程度较深，有效地去除了噪声；本发明方法在均值的指标上和原有SAR图像是最接近的，表明本发明在SAR图像的平均亮度的保持上优于以上已有的方法，而且本发明的标准差相比于其他方法是最低的，也表明了本发明方法的相干斑抑制效果明显；在均值比的指标上，本发明方法最接近1，表明本方法的辐射度失真是最小的；综上所述，本发明方法取得了较好的去噪指标和去噪效果。

本发明与已有方法对Albuquerque地区SAR图像的去噪性能指标对比表如表2。

表2各方法对Albuquerque地区SAR图像去噪性能对比

从表2中可以看出，本发明方法取得了优于表中所列举的方法的去噪指标，表明本发明方法的去噪效果更好。

总之，本发明在去噪性能和视觉效果方面较之小波等方法均有所改善，且较完整地保持了图像边缘等细节信息。

Claims (2)

1.一种基于方向波域混合高斯模型的SAR图像噪声抑制方法，包括如下步骤：

(1)对SAR图像进行卡通纹理模型分解，选择最优的迭代次数，得到SAR图像的卡通部分和纹理部分；

(2)保持SAR图像卡通部分不变，将纹理部分进行如下改进的Directionlet变换；

(2a)选定SAR图像纹理部分的变换方向和队列方向，其斜率构成采样矩阵M_Λ：

$M_{\mathrm{\Λ}}=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]a_1,a_2,b_1,b_2\⋐Z$

其中，沿向量d₁的方向称为纹理部分的变换方向，其斜率表示为r₁＝b₁/a₁，沿向量d₂的方向称为纹理部分的队列方向，其斜率表示为r₂＝b₂/a₂；

(2b)对纹理部分进行M_Λ采样，得到|det(M_Λ)|个陪集，再对得到的陪集分别进行同样的二次陪集分解，得到2|det(M_Λ)|个子陪集，将零点去除；

(2c)对每个子陪集沿变换和队列方向分别进行两次和一次的一维正交小波变换，得到相应的高频和低频子带，在低频子带的变换和队列方向再分别进行两次和一次一维正交小波变换，获得14个高频子带和1个低频子带系数均为W_Y(i，j)，

表示每个系数的坐标；

(3)对变换后得到的各子陪集的高频子带系数进行混合高斯模型的拟合，并采用EM算法估计混合高斯模型中的参数；

(4)估计各子陪集的高频子带系数的噪声方差

(5)利用Bayes最小均方误差估计法，计算各子陪集的高频子带的无噪系数

(6)对无噪系数进行改进的Directionlet逆变换，得到去噪后的纹理部分；

(7)将去噪后的纹理部分与保持不变的卡通部分相加，得到去噪后的SAR图像。

2.根据权利要求1所述的SAR图像噪声抑制方法，其中步骤(1)所述的对SAR图像进行卡通纹理模型分解，按如下步骤进行：

(2A)将SAR图像f∈L²(R²)表示为：

f(x，y)＝u(x，y)+v(x，y)

其中u∈BV(R²)，由分片光滑的区域及清晰的边缘构成，被视为SAR图像的卡通部分，v就是剩余的纹理和噪声，称之为纹理部分；

(2B)当v属于一个定义的Banach空间，得到模型：

$\underset{u\∈\mathrm{BV}}{\inf }\{E\left(u\right)=\∫|\▿u|+\mathrm{\λ}{\left|\right|v\left|\right|}_u,f=u+v\}$

其中E(u)表示u的能量，通过保证u∈BV(R²)时，进行迭代计算使得u的能量泛函最小化，以达到对图像f的卡通纹理模型分解，通过试验，最优的迭代次数一般为50～60次。

Images