CN101661611B - 基于贝叶斯非局部均值滤波器的实现方法 - Google Patents

Abstract

本发明公开了一种基于贝叶斯非局部均值滤波器的实现方法，它属于图像处理技术领域，主要克服现有SAR图像去斑结果中细节保持与平滑程度的矛盾问题。其实现过程是：(1)对输入SAR图像的一个像素点，根据距离公式计算该像素点与其搜寻区域内所有像素点的距离；(2)根据计算的所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值；(3)根据根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到该像素点的恢复值；(4)按照上述步骤，计算出每个像素点的恢复值，得到去斑图像。本发明相对于其它的一些经典的去斑方法能够更好平滑斑点噪声同时保持SAR图像的边缘和纹理细节，可用于对SAR图像的去斑处理。

Description

基于贝叶斯非局部均值滤波器的实现方法

技术领域

本发明属于图像处理技术领域，具体地说是一种滤波器，可用于对SAR图像的去斑处理。

背景技术

合成孔径雷达SAR所成的图像具有全天候、全天时、高分辨率和强大的穿透能力等特点，因此，这种图像被广泛的应用到目标识别，变换检测和水面监视。然而，SAR图像被乘性噪声所腐蚀，这种噪声来自后向散射雷达反射的连续干扰。这种斑点噪声毁坏了SAR图像辐射测量的分辨率，同时影响到背景分析的性能和理解任务。

降斑方法的目标就是在去除噪声的同时保留图像的特征信息，像纹理，边缘和点状目标。但是由于斑点噪声的乘性背景，这一目标的实现非常困难。目前，大量斑点滤波的方法已经被提出。空域统计类去斑方法一般先假定噪声的乘性模型，然后基于像素邻域窗口的局部统计特性进行滤波。基于空域滤波的典型算法有Lee滤波，Kuan滤波，Frost滤波及他们的增强版。这些方法在同质区域取均值，对变化较快的点采取保留的策略，其缺点是或者过平滑图像的纹理或者不能有效的去除边缘周围的噪声。图像结构信息如边缘、线性体、点等目标会在一定程度上被模糊或滤除。在最近几年，针对加性高斯噪声的许多去噪方法被用来对SAR图像去噪，例如小波软阈值的方法。这些方法首先通过对数运算将斑点噪声的乘性模型转为加性模型，然后直接利用自然图像去噪算法。基于小波、Contourlet等变换的SAR图像去斑算法由于其多尺度、多分辨等优点获得了广泛的应用。但是变换域去斑算法本质上仍是基于固定窗口的滤波，在图像的边缘、线等区域会产生吉布斯现象。

发明内容

本发明的目的在于克服已有技术的不足，提出了一种基于贝叶斯非局部均值滤波的实现方法，以实现对SAR图像降斑中边缘和平滑区域的兼顾，提高图像去斑效果。

为实现上述目的，本发明适于强度SAR图像滤波器的实现方法，包括如下步骤：

(1)对输入强度SAR图像的一个像素点，根据距离公式计算该像素点与其搜索邻域内所有像素点的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)),$

其中，L为强度SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

(2)根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，h是一个平滑参数，h²＝1/L，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，N×N表示以像素点x_i为中心的搜寻区域的大小；

(3)根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),$

其中，权值w(x_i，x_j)满足：0≤w(x_i，x_j)≤1， $\underset{j}{\mathrm{\Σ}}w(x_i,x_j)=1;$ v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值；

(4)按照上述步骤，计算出每个像素点的恢复值，用求得的恢复值取代原强度SAR图像中所有点的灰度值，得到强度SAR图像的去斑图像。

为实现上述目的，本发明适于幅度SAR图像滤波器的实现方法，包括如下步骤：

1)对输入幅度SAR图像的一个像素点，根据距离公式计算该像素点与其搜索邻域内所有像素点的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)),$

其中，L为幅度SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

2)根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，h是一个平滑参数，h²＝1/L，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，N×N表示以像素点x_i为中心的搜寻区域的大小；

3)根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),$

其中，权值w(x_i，x_j)满足：0≤w(x_i，x_j)≤1， $\underset{j}{\mathrm{\Σ}}w(x_i,x_j)=1;$ v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值；

4)按照上述步骤，计算出每个像素点的恢复值，用求得的恢复值取代原幅度SAR图像中所有点的灰度值，得到幅度SAR图像的去斑图像。

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

1.本发明是在空域中进行的，实现过程简单，且可以并行实现。

2.本发明由于利用贝叶斯框架和SAR图像斑点模型推导得到的一种新的距离公式，能够准确地计算SAR图像中以待估计像素点x_i为中心的块和以搜寻区域内像素点x_j为中心的块之间的相似性，进而能准确的计算出待估计像素点x_i的恢复值。

3.本发明由于使用了新提出的距离公式，能够准确计算待估计像素点x_i的恢复值，因此，能够较好平滑斑点噪声，同时保持和恢复SAR图像的边缘和纹理细节。

附图说明

图1是本发明的流程图；

图2是本发明实验输入的两幅SAR图像；

图3是用现有增强Lee滤波算法对输入的两幅SAR图像进行去斑的结果图；

图4是用现有增强Gamma MAP滤波算法对输入的两幅SAR图像进行去斑的结果图；

图5是用现有小波软阈值滤波算法对输入的两幅SAR图像进行去斑的结果图；

图6是用现有原始非局部均值滤波算法对输入的两幅SAR图像进行去斑的结果图；

图7是用本发明对输入的两幅SAR图像进行去斑的结果图；

图8是图6中两幅去斑结果图的局部放大图；

图9是图7中两幅去斑结果图的局部放大图。

具体实施方式

参照附图1，本发明给出以下两种实施例。

实施例1：适于强度SAR图像滤波器的实现方法，包括如下步骤：

步骤1，构造强度SAR图像的距离公式。

1.1)根据强度SAR图像的噪声模型I＝RS_I，及斑点噪声S_I概率密度函数： $p_{S_I}\left(s\right)=\frac{L^L{s}^{L-1}{e}^{-\mathrm{Ls}}}{\mathrm{\Γ}\left(L\right)},(s\≥0),$ 得到条件概率密度函数为：

$p_{I|R}\left(x\right|y)=\frac{1}{\mathrm{\Γ}\left(L\right)}{\left(\frac{L}{y}\right)}^L{x}^{L-1}{e}^{-\mathrm{Lx}/y},$

其中，I为图像强度的观测值，R是图像强度的真实值，S_I表示斑点噪声，L表示强度图的视数，Γ(L)表示sigma函数；

1.2)利用得到的条件概率密度函数，假设所有条件概率之间相互独立，得到v(x_i)和v(x_j)之间的条件概率：

$p\left(v\left(x_i\right)\right|v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Π}}}p\left(v^{\left(m\right)}\left(x_i\right)\right|v^{\left(m\right)}\left(x_j\right))$

$={\left(\frac{L^L}{\mathrm{\Γ}\left(L\right)}\right)}^{M\×M}\exp (-\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{{\mathrm{Lv}}^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)})\exp \left(\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}\ln \frac{{\left(v^{\left(m\right)}\left(x_i\right)\right)}^{L-1}}{{\left(v^{\left(m\right)}\left(x_j\right)\right)}^L}\right)$

$\∝\exp (-\frac{1}{1/L}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))),$

$\∝\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)))$

其中，h为平滑参数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

1.3)把得到的v(x_i)和v(x_j)两个向量之间的条件概率带入到现有的非局部均值滤波的贝叶斯公式：

$\mathrm{BNL}\left(v\left(x_i\right)\right)=\frac{\frac{1}{N\×N}\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)p\left(v\left(x_i\right)\right|v\left(x_j\right))}{\frac{1}{N\×N}\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}p\left(v\left(x_i\right)\right|v\left(x_j\right))},$

$=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j)$

得到基于贝叶斯非局部均值滤波的权值公式：

$w(x_i,x_j)=\frac{p\left(v\left(x_i\right)\right|v\left(x_j\right))}{\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}p\left(v\left(x_i\right)\right|v\left(x_j\right))},$

$\∝\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)))$

该权值公式中，归一化常数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))),$

权值公式中幂指数因子： $\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))$ 起到衡量待估计像素点x_i和搜寻区域内像素点x_j相似性的作用，因此，它就是我们所构造的适应强度图像的距离公式，即：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right));$

步骤2，输入强度SAR图像，计算其中一个像素点与其搜索区域内所有像素点的距离。

1.1)输入强度图像如附图2(b)所示，以该SAR图像其中一个待估计像素点x_i为中心取M×M大小的块，实验中，M的值取为7，即块的大小是7×7；

1.2)以像素点x_i为中心取出一个大小为N×N的搜寻区域，实验中，N的值取为21，即搜寻区域的大小为21×21；

1.3)以搜寻区域中每一个像素点x_j为中心取M×M大小的块；

1.4)利用本发明构造的距离公式计算强度SAR图像中像素点x_i和x_j之间的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)).---\left(1\right)$

步骤3，根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值。

搜寻区域内的像素点x_j的权值计算公式如下：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),---\left(3\right)$

其中，w(x_i，x_j)为像素点x_j相对于x_i的权值，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right)))---\left(4\right)$

h为平滑参数，由距离公式的构造过程可以发现，理论上h²＝1/L，但是，这时并不能取得满意的去噪结果。由于强度图像噪声的方差σ²＝1/L，所以假定参数h和噪声标准差σ之间是一种线性正比关系：h＝kσ。实验中发现，对于输入的四视强度SAR图像取 $k=\sqrt{4.5}$ 能获得满意的去噪结果。

步骤4，根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),---\left(5\right)$

其中，w(x_i，x_j)为像素点x_j相对于x_i的权值，且满足：0≤w(x_i，x_j)≤1， $\underset{j}{\mathrm{\Σ}}w(x_i,x_j)=1;$ v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值。

步骤5，重复步骤2到步骤4，计算出每个像素点的恢复值，用求得的恢复值取代原强度SAR图像中所有点的灰度值，得到强度SAR图像的去斑图像，如图7(b)所示。

实施例2：适于幅度SAR图像滤波器的实现方法，包括如下步骤：

步骤A，构造幅度SAR图像的距离公式。

A.1)根据幅度SAR图像的噪声模型A＝RS_A，及斑点噪声S_A概率密度函数： $p_{S_A}\left(s\right)=\frac{{2L}^L{s}^{2L-1}{e}^{-{\mathrm{Ls}}^2}}{\mathrm{\Γ}\left(L\right)},(s\≥0),$ 得到条件概率密度函数为：

$p_{A|R}\left(x\right|y)=\frac{2L^L}{\mathrm{\Γ}\left(L\right)}\frac{x^{2L-1}}{y^{2L}}{e}^{-L{(x/y)}^2},$

其中，A为图像幅度的观测值，R是图像幅度的真实值，S_A表示斑点噪声，L表示强度图的视数，Γ(L)表示sigma函数；

A.2)利用得到的条件概率密度函数，假设所有条件概率之间相互独立，得到v(x_i)和v(x_j)之间的条件概率：

$p\left(v\left(x_i\right)\right|v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Π}}}p\left(v^{\left(m\right)}\left(x_i\right)\right|v^{\left(m\right)}\left(x_j\right))$

$={\left(\frac{{2L}^L}{\mathrm{\Γ}\left(L\right)}\right)}^{M\×M}\exp (-L\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}{\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2)\exp \left(\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}\ln \frac{{\left(v^{\left(m\right)}\left(x_i\right)\right)}^{2L-1}}{{\left(v^{\left(m\right)}\left(x_j\right)\right)}^{2L}}\right)$

$\∝\exp (-\frac{1}{1/L}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)))$

$\∝\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))),$

其中，h为平滑参数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

A.3)把得到的v(x_i)和v(x_j)两个向量之间的条件概率带入到现有的非局部均值滤波的贝叶斯公式：

$\mathrm{BNL}\left(v\left(x_i\right)\right)=\frac{\frac{1}{N\×N}\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)p\left(v\left(x_i\right)\right|v\left(x_j\right))}{\frac{1}{N\×N}\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}p\left(v\left(x_i\right)\right|v\left(x_j\right))},$

$=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j)$

得到基于贝叶斯非局部均值滤波的权值公式：

$w(x_i,x_j)=\frac{p\left(v\left(x_i\right)\right|v\left(x_j\right))}{\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}p\left(v\left(x_i\right)\right|v\left(x_j\right))},$

$\∝\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))),$

该权值公式中，归一化常数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))),$

权值公式中幂指数因子： $\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right))$ 起到衡量待估计像素点x_i和搜寻区域内像素点x_j相似性的作用，因此，它就是我们所构造的适应幅度图像的距离公式，即：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)).$

步骤B，输入幅度SAR图像，计算其中一个像素点与其搜索区域内所有像素点的距离。

B.1)输入幅度图像如附图2(a)所示，以该SAR图像其中一个待估计像素点x_i为中心取M×M大小的块，实验中，M的值取为7，即块的大小是7×7；

B.2)以像素点x_i为中心取出一个大小为N×N的搜寻区域，实验中，N的值取为21，即搜寻区域的大小为21×21；

B.3)以搜寻区域中每一个像素点x_j为中心取M×M大小的块；

B.4)利用本发明构造的距离公式计算幅度SAR图像中像素点x_i和x_j之间的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)).---\left(1\right)$

步骤C，根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值。

搜寻区域内的像素点x_j的权值计算公式如下：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))3,---\left(3\right)$

其中，w(x_i，x_j)为像素点x_j相对于x_i的权值，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),---\left(4\right)$

h为平滑参数，由距离公式的构造过程可以发现，理论上h²＝1/L，但是，这时并不能取得满意的去噪结果。由于幅度图像噪声的方差σ²＝(4/π-1)/L，所以假定参数h和噪声标准差σ之间是一种线性正比关系：h＝kσ。实验中发现，对于输入的二视幅度SAR图像取k＝6能获得满意的去噪结果。

步骤D，根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),---\left(5\right)$

其中，w(x_i，x_j)为像素点x_j相对于x_i的权值，且满足：0≤w(x_i，x_j)≤1， $\underset{J}{\mathrm{\Σ}}w(x_i,x_j)=1;$ v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值。

步骤E，重复步骤B到步骤D，计算出每个像素点的恢复值，用求得的恢复值取代原幅度SAR图像中所有点的灰度值，得到幅度SAR图像的去斑图像，如图7(a)所示。本发明效果可以通过以下实验进一步证实：

一.实验条件和内容

实验条件：实验所使用的输入图像是图2所示，其中，图2(a)是两视幅度SAR图像，图2(b)是四视强度SAR图像。图2中，标号1、2、3、4和5标示出了待评价去斑效果的同质区域。实验中，各种滤波方法都是使用matlab语言编程实现。

实验内容：在上述实验条件下，分别使用窗口大小为5×5的增强Lee滤波，窗口大小为7×7的增强Gamma MAP滤波，原始非局部均值滤波方法和本发明进行实验。其中，小波软阈值方法使用的是DB8小波基，原始非局部均值滤波方法是先将图进行对数变换，这时噪声就变成标准差为σ′的加性噪声，再使用7×7大小的块，21×21大小的搜寻区域进行滤波，其中平滑参数h＝0.3σ′。

二.实验结果

使用窗口大小为5×5的增强Lee滤波方法的去斑结果如图3所示，其中图3(a)是两视幅度SAR图像去斑结果，图3(b)是四视强度SAR去斑结果。从结果图可以看出，该滤波方法使点目标亮度被压缩，使边缘模糊且不能保持边缘的连续性，同时同质区域不够平滑；

使用窗口大小为7×7的增强Gamma MAP滤波方法去斑结果如图4所示，其中图4(a)是两视幅度SAR图像去斑结果，图4(b)是四视强度SAR去斑结果。从去斑结果图可以看出，该滤波方法在同质区域比增强Lee滤波更平滑，但是，它在点目标保护和边缘连续性保持上仍和增强Lee滤波一样都存在缺陷；

使用小波软阈值滤波方法去斑结果如图5所示，其中图5(a)是两视幅度SAR图像去斑结果，图5(b)是四视强度SAR去斑结果。从结果图可以看出，该滤波方法在同质区域具有较好的平滑性能，但是，它却不能保持边缘的连续性，因为在边缘处总是出现吉布斯现象；

使用原始非局部均值滤波方法去斑结果如图6所示，其中图6(a)是两视幅度SAR图像去斑结果，图6(b)是四视强度SAR去斑结果。该两幅去斑结果图的细节放大图分别如图8(a)和8(b)所示。从图6的去斑结果图及其细节放大图可以看出，该滤波方法对同质区域具有一定的平滑能力，且能一定程度上保护图像中的纹理，相对于前三种方法，该滤波方法能相对较好的保护点目标和保持边缘的连续性，但是在边缘处仍然有不连续的地方；

使用本发明的去斑结果如图7所示，其中图7(a)是两视幅度SAR图像去斑结果，图7(b)是四视强度SAR去斑结果。该两幅去斑结果图的细节放大图分别如图9(a)和9(b)所示。从图7的去斑结果图及其细节放大图可以看出，对比于上述的各种滤波方法，本发明能够很好的保护点目标，且在平滑同质区域的同时还能保护同质区域的纹理信息，同时，还能保持边缘的连续，另外，对植被等细节信息保持方面也有很好的性能。

用图2(a)中标记为1、2、3和图2(b)中标记为4、5的五个同质区域的均值、标准差和等效视数ENL作为去斑结果的定量评价指标。上述四种现有的滤波方法去斑结果和本发明的去斑结果在五个标记的同质区域上均值、标准差和ENL列在表1中。

表1  各种去斑结果对比

从表1中可以发现，现有的增强Lee和增强Gamma MAP两种滤波方法具有较好的均值保持能力，但是，它们的标准差和ENL两项评价指标远落后于其它几种滤波方法；小波软阈值滤波方法的对比前两种方法有较低的标准差，但是其均值保持能力较差；原始非局部均值滤波方法和小波软阈值方法在各种评价指标上各有优劣，可谓旗鼓相当；本发明相对前面所述的各种方法，其标准差和ENL两项评价指标都占有绝对的优势，本发明在均值保持上比增强Lee和增强Gamma MAP两种滤波方法要差一些，但是要比小波软阈值和原始非局部均值两种滤波方法要稍好一点。

以上实验结果表明，本发明相对于其它的一些经典的去斑方法具有更好的性能，能够更好平滑斑点噪声同时保持SAR图像的边缘和纹理细节。

Claims (2)

1.一种基于贝叶斯非局部均值滤波器的实现方法，包括如下步骤：

(1)对输入强度SAR图像的一个像素点，根据距离公式计算该像素点与其搜索邻域内所有像素点的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}+\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{L-1}{L}\ln (v^{\left(m\right)}\left(x_i\right)\left)\right),$

其中，L为强度SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，实验中，M的值取为7，N值取为21，即M小于N，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

(2)根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，h是一个平滑参数，h²＝1/L，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，N×N表示以像素点x_i为中心的搜寻区域的大小；

(3)根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),$

其中，权值w(x_i，x_j)满足：0≤w(x_i，x_j)≤1，

v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值；

(4)按照上述步骤，计算出每个像素点的恢复值，用求得的恢复值取代原强度SAR图像中所有点的灰度值，得到强度SAR图像的去斑图像。

2.一种基于贝叶斯非局部均值滤波器的实现方法，包括如下步骤：

1)对输入幅度SAR图像的一个像素点，根据距离公式计算该像素点与其搜索邻域内所有像素点的距离：

$d(v\left(x_i\right),v\left(x_j\right))=\underset{m=1}{\overset{M\×M}{\mathrm{\Σ}}}({\left(\frac{v^{\left(m\right)}\left(x_i\right)}{v^{\left(m\right)}\left(x_j\right)}\right)}^2+2\ln \left(v^{\left(m\right)}\left(x_j\right)\right)-\frac{2L-1}{L}\ln \left(v^{\left(m\right)}\left(x_i\right)\right)),$

其中，L为幅度SAR图像视数，x_i是待估计的像素点，x_j是以x_i为中心的大小为N×N的搜寻区域中的像素点，v(x_i)是以x_i为中心的M×M大小的块的灰度值所形成的向量，实验中，M的值取为7，N值取为21，即M小于N，v(x_j)是以x_j为中心的M×M大小的块的灰度值所形成的向量，v^(m)(x_i)表示向量v(x_i)的第m个元素，v^(m)(x_j)表示向量v(x_j)的第m个元素；

2)根据计算出的搜寻区域内所有像素点的距离，使用权值计算公式，计算搜寻域内所有像素点的权值：

$w(x_i,x_j)=\frac{1}{Z\left(x_i\right)}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，h是一个平滑参数，h²＝1/L，Z(x_i)是归一化系数：

$Z\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\frac{1}{h^2}d(v\left(x_i\right),v\left(x_j\right))),$

其中，N×N表示以像素点x_i为中心的搜寻区域的大小；

3)根据计算出的搜寻区域内所有像素点的权值，对搜寻区域内所有像素点进行加权平均，得到待估计像素点x_i的恢复值：

$\mathrm{BNLv}\left(x_i\right)=\underset{j=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(x_j\right)w(x_i,x_j),$

其中，权值w(x_i，x_j)满足：0≤w(x_i，x_j)≤1，

v(x_j)为原图中x_j处的灰度值，BNLv(x_i)为加权平均后x_i点的恢复值；

4)按照上述步骤，计算出每个像素点的恢复值，用求得的恢复值取代原幅度SAR图像中所有点的灰度值，得到幅度SAR图像的去斑图像。

Images