CN102136134A - 基于mrf先验的sar图像去斑方法 - Google Patents

Abstract

一种图像处理技术领域基于MRF先验的SAR图像去斑方法，包括如下步骤：1.输入一副待去斑的SAR图像；2.选取中心像素块；3.确定搜索窗；4.选取相似块；5.计算相似性权值；6.遍历搜索窗；7.计算待估计像素点的恢复值；8.遍历整幅SAR图像；9.恢复整幅SAR图像；10.迭代去斑；11.获得最终去斑结果。本发明充分利用了图像的先验信息，相对于目前空域滤波的去斑方法，能够更好的平滑斑点噪声，同时保持SAR图像的边缘和纹理细节，可用于对SAR图像的去斑处理。

Description

基于MRF先验的SAR图像去斑方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种雷达图像处理技术的滤波方法，该方法使用马尔科夫(MRF)随机场提取先验信息，以贝叶斯非局部均值(BNL)滤波为去斑方法，可用于对合成孔径雷达(SAR)图像进行去斑处理。

背景技术

目前，基于空域滤波的去斑效果较好的算法有原始非局部均值方法，改进sigma滤波方法和PPB滤波方法。

Buades等人在“Buades，B.Coll，and J.M.Morel.A non-local algorithm for image denoising.Proc.IEEE CVPR，2005，2：60-65.”中提出了一种非局部均值去斑滤波方法。该方法不是用图像中单个像素的灰度值进行比较，而是对该像素周围的窗口中整个灰度的分布状况进行比较，根据灰度分布的相似性贡献权值。实践证明，使用该方法获得的灰度分布相似性权值不准确，导致图像结构信息如边缘、线性体、点等目标会在一定程度上被模糊或滤除。

J.S.Lee等人在“J.S.Lee，J.H.Wen，T.L.Ainsworth，K.S.Chen，and A.J.Chen，“Improved sigma filter for speckle filtering of SAR imagery，”IEEE Trans.Geosci.Remote Sens.，vol.47，no.1，Jan 2009，202-213.”中提出了改进sigma去斑滤波方法。该方法使用新的sigma范围选取窗口内的像素点，然后对选取的像素点采用最小均方误差(MMSE)估计中心点像素的值。实际使用中证明该方法在一定程度上压缩了强反射点目标的亮度，且不能很好平滑同质区域的噪声，导致去斑效果不理想。

C.A.Deledalle等人在“C.A.Deledalle，L.Denis，and F.Tupin.Iterative weighted maximum likelihood denoising with probabilistic patch-based weights.IEEE Trans.Image Processing，vol.18，no.12，Dec.2009，2661-2672.”中提出了一种在最大似然框架下迭代去斑的滤波方法。该方法在最大似然估计框架下得到加权平均公式，并定义块权值的计算公式，通过迭代方法逐步修改先验信息，最终收敛至最佳去斑结果。此方法在使用中需迭代25次，因此其复杂度很高，降低了该方法的实用性。

发明内容

本发明的目的在于克服现有技术的不足，先使用MRF模型提取图像的先验信息，再根据贝叶斯非局部均值滤波框架对图像去斑的方法，提出了一种基于MRF先验的BNL均值滤波器对SAR图像空域去斑方法，以实现对SAR图像降斑中边缘和平滑区域的兼顾，提高图像去斑效果。

本发明的具体步骤包括如下：

(1)输入一副待去斑的SAR图像；

(2)选取中心像素块：在待去斑的SAR图像中，逐行扫描选取一个像素点为待估计像素点，以待估计像素点为中心，以固定长度为块半径，选取一个正方形的中心像素块；

(3)确定搜索窗：以选取的待估计像素点为中心，以固定长度为搜索窗半径，选取一个正方形的搜索窗；

(4)选取相似块：在搜索窗中逐行扫描选取一个像素点作为相似点，以此相似点为中心，选取一个与中心像素块等大小的像素块作为相似块；

(5)计算相似性权值：

5a)计算相似性条件概率：依据条件概率密度函数，计算在步骤(4)所选取的相似块的条件下，中心像素块出现的条件概率；

5b)根据基于MRF先验的相似性度量公式，计算先验相似性：

$g(\·)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right)$

其中，g(·)为先验相似性；β为平滑参数；η_s，k和η_t，k分别为u_s，k和u_t，k的四邻域u_i为η_s，k中的像素点；u_j为η_t，k中的像素点

5c)将相似性条件概率与先验相似性相乘获得相似性权值；

(6)遍历搜索窗：判断搜索窗内的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤(4)、(5)，计算出整个搜索窗中所有的相似块与中心块的相似性权值；

(7)计算待估计像素点的恢复值：根据步骤(6)中求出的相似性权值，将获得的每个相似块的相似性权值和其对应的相似点的灰度值，加权求和得出待估计像素点的恢复值；

(8)遍历整幅SAR图像：判断整幅SAR图像中的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤(2)、(3)、(4)、(5)、(6)、(7)，计算出整幅SAR图像中所有点的恢复值；

(9)恢复整幅SAR图像：根据步骤(8)中获得的恢复值，恢复出整幅SAR图像；

(10)迭代去斑：使用最新恢复出的SAR图像计算相似性权值，将获得的相似性权值加权求和恢复出整幅图像；

(11)获得最终去斑结果：重复执行步骤(10)，直到去斑图达到预想要求，中止迭代。

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

1.本发明是在空域中直接对像素点进行处理的，实现过程简单，且可以并行实现。

2.本发明由于利用BNL均值滤波方法和SAR图像斑点模型以及MRF先验模型推导得到的一种新的相似性权值，能够准确地计算出SAR图像中以待估计像素点为中心的中心块和以搜寻窗内的相似点为中心的相似块之间的相似性，进而能准确的恢复出待估计像素点，较好的保持了图像中的线目标和点目标，并能很好的平滑同质区域的噪声，从而提高了图像的去斑效果。

3.本发明由于充分利用了图像中的先验信息，因此能够准确的计算相似性权值，从而能够准确计算待估计像素点的恢复值，因此，能够较好平滑斑点噪声，同时保持和恢复SAR图像的边缘和纹理细节，图像的去斑效果明显增强。

附图说明

图1为本发明的流程图；

图2为基于空域滤波的去斑方法的效果图；

图3为本发明的效果图。

具体实施方式

下面结合附图对本发明做进一步的描述。

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

步骤1，输入一副待去斑的SAR图像。

步骤2，选取中心像素块。

选取中心像素块主要是为了确定待估计像素点的领域信息。

1a)选取待估计像素点。在待去斑的SAR图像v中，逐行扫描选取一个像素点t作为待估计像素点；

1b)确定中心像素块。以待估计像素点为中心，以固定长度2～5像素点为块半径，选取一个大小为M×M正方形的中心像素块。本发明实施例中以像素点t为中心，选取一个大小为7×7的中心像素块v_t。

步骤3，确定搜索窗。

由于搜索窗的大小直接决定了相似点的选取范围，它必须包含较多的相似点。确定搜索窗时，以选取的待估计像素点t为中心，以固定长度7～15像素点为搜索窗半径，选取大小为N×N的搜索窗。本发明实施例中选取大小为21×21的搜索窗。

步骤4，选取相似块。

在步骤2获得的搜索窗中逐行扫描选取一个像素点s作为相似点，以此相似点为中心，选取一个大小为M×M相似块u_s。相似块的大小一定要与中心相似块的大小相同，它描述了相似点s的领域信息。

步骤5，计算相似性权值。

5a)计算相似性条件概率。根据条件概率密度函数，计算在已知相似块u_s的条件下，中心像素v_t块出现的条件概率。相似性条件概率部分决定了搜索窗中的相似点在加权平均中对恢复待估计相似点的贡献，根据强度和幅度的SAR图像，它有不同的计算公式。

适用于强度SAR图像的条件概率密度函数为：

$p\left(v_t\right|u_s)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}\ln v_{t,k}))$

适用于幅度SAR图像的条件概率密度函数为：

$p\left(v_t\right|u_s)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}\ln v_{t,k}))$

其中，v_t为中心像素块；u_s为相似块；

为平滑参数

v_t，k为v_t中逐行扫描的第k个像素点

u_s，k为u_s中逐行扫描的第k个像素点

L为SAR图像视数

5b)计算先验相似性。先验相似性主要是充分利用了相似块与中心像素块之间的先验信息，根据它们之间的先验信息来计算它们之间的相似度。

计算先验相似性需使用先验信息，其相似性度量公式如下：

$g(\·)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right)$

其中，g(·)为先验相似性；

β为平滑参数

η_s，k和η_t，k分别为u_s，k和u_t，k的四邻域

u_i为η_s，k中的像素点

u_j为η_t，k中的像素点

5c)将相似性条件概率与先验相似性相乘获得相似性权值w(t，s)。相似性权值描述了相似点与待估计像素点之间的相似度。根据强度和幅度的SAR图像，它有不同的计算公式。

适用于强度SAR图像的相似性权值公式为：

$w(t,s)=p\left(v_t\right|u_s)\×g(\·)$

$=\frac{1}{Z\left(t\right)}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

其中，Z(t)为归一化因子：

$Z\left(t\right)=\underset{s=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

在实施例中

β＝1×10⁵能获得满意的去斑效果。

适用于幅度SAR图像的相似性权值公式为：

$w(t,s)=p\left(v_t\right|u_s)\×g(\·)$

$=\frac{1}{Z\left(t\right)}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

其中，Z(t)为归一化因子：

$Z\left(t\right)=\underset{s=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

在实施例中

β＝3×10⁵能获得满意的去斑效果。

步骤6，遍历搜索窗。判断搜索窗内的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤4、5，计算出整个搜索窗中所有的相似块与中心块的相似性权值。

步骤7，计算待估计像素点的恢复值。根据步骤6中求出的相似性权值，将获得的每个相似块的相似性权值和其对应的相似点的灰度值，进行加权求和得出待估计像素点的恢复值。

$\mathrm{MBNLv}\left(t\right)=\underset{s=1}{\overset{N\×N}{\mathrm{\Σ}}}v\left(s\right)w(t,s)$

其中，MBNLv(t)为加权平均后t点的恢复值

v(s)为原图中s处的灰度值

w(t，s)为像素点s相对于t的权值，且满足：0≤w(t，s)≤1，

步骤8，遍历整幅SAR图像。判断整幅SAR图像中的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤2、3、4、5、6、7，计算出整幅SAR图像中所有点的恢复值。

步骤9，恢复整幅SAR图像：根据步骤8中获得的恢复值，获得整幅SAR图像的第一次恢复值u。

步骤10，迭代去斑。使用最新恢复出的SAR图像u，重复步骤2、3、4、5、6、7、8、9，恢复出整幅图像的更新值。

步骤11，获得最终去斑结果：重复执行步骤10，直到去斑图达到预想要求，中止迭代。

下面结合图2、3对本发明的效果做进一步的描述

图2(a)和(b)表示原始非局部均值滤波方法对两幅SAR图像的去斑效果图，图2(c)和(d)表示改进sigma滤波方法对两幅SAR图像的去斑效果图，图2(e)和(f)表示PPB滤波方法对两幅SAR图像的去斑效果图。

图3(a)和(b)表示输入待去斑的SAR图像，其中图3(a)表示待去斑的两视幅度SAR图像，图3(b)表示待去斑的四视强度SAR图像。图3(a)和(b)中，标号1、2、3、4和5标示出了待评价去斑效果的同质区域。图3(c)和(d)表示本发明两次迭代对两幅真实SAR图像去斑的效果图，为了证明我们滤波方法有较好的去斑效果，分别使用窗口大小为5×5的原始非局部均值滤波方法、改进sigma滤波方法和PPB滤波方法与本章滤波方法实验结果进行对比。将去斑后区域1、区域2、区域3、区域4和区域5对应的均值、标准差、等效视数ENL作为评价指标，比较各种空域滤波方法的优劣。其比较结果如表1所示：

表1各种去斑结果对比

从表1可以发现，本发明在均值、标准差和等效视数各方面都有较佳的表现，本方法的去斑效果显示出了同质区域斑点抑制能力和一定的保护点线等强反射目标的能力。同时，本发明与基于边缘、同质区域、点目标分类的空域滤波方法相比，在对细节信息，像线和纹理进行恢复时比其它方法表现出了更好的效果。本发明最明显的优点在于仅仅利用了图像中的先验信息，也能把点目标的亮度保持的相当好。从图3(c)和(d)中可以看出，本发明对同质区域的平滑性能很好，因此可能会因为过平滑而丢失细小纹理信息，但是还是能把大片区域的地物轮廓恢复出来。另外，如果对细小纹理感兴趣可以减少迭代次数，降低对同质区域的平滑力，就能保护好细小纹理结构信息。相对于其它的一些空域的去斑方法具有更好的性能，能够更好平滑斑点噪声同时保持SAR图像的边缘和纹理细节。

Claims (7)

1.一种基于MRF先验的SAR图像去斑方法，包括如下步骤：

(1)输入一副待去斑的SAR图像；

(2)选取中心像素块：在待去斑的SAR图像中，逐行扫描选取一个像素点为待估计像素点，以待估计像素点为中心，以固定长度为块半径，选取一个正方形的中心像素块；

(3)确定搜索窗：以选取的待估计像素点为中心，以固定长度为搜索窗半径，选取一个正方形的搜索窗；

(4)选取相似块：在搜索窗中逐行扫描选取一个像素点作为相似点，以此相似点为中心，选取一个与中心像素块等大小的像素块作为相似块；

(5)计算相似性权值：

5a)计算相似性条件概率：依据条件概率密度函数，计算在步骤(4)所选取的相似块的条件下，中心像素块出现的条件概率；

5b)根据基于MRF先验的相似性度量公式，计算先验相似性：

$g(\·)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right)$

其中，g(·)为先验相似性；β为平滑参数；η_s，k和η_t，k分别为u_s，k和u_t，k的四邻域u_i为η_s，k中的像素点；u_j为η_t，k中的像素点

5c)将相似性条件概率与先验相似性相乘获得相似性权值；

(6)遍历搜索窗：判断搜索窗内的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤(4)、(5)，计算出整个搜索窗中所有的相似块与中心块的相似性权值；

(7)计算待估计像素点的恢复值：根据步骤(6)中求出的相似性权值，将获得的每个相似块的相似性权值和其对应的相似点的灰度值，加权求和得出待估计像素点的恢复值；

(8)遍历整幅SAR图像：判断整幅SAR图像中的点是否搜索完，若搜索完则进行下一步骤，若没有搜索完则重复步骤(2)、(3)、(4)、(5)、(6)、(7)，计算出整幅SAR图像中所有点的恢复值；

(9)恢复整幅SAR图像：根据步骤(8)中获得的恢复值，恢复出整幅SAR图像；

(10)迭代去斑：使用最新恢复出的SAR图像计算相似性权值，重复步骤2、3、4、5、6、7、8、9，恢复出整幅图像的更新值；

(11)获得最终去斑结果：重复执行步骤(10)，直到去斑图达到预想要求，中止迭代。

2.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤(2)中的固定长度为2～5个像素点。

3.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤(3)中的固定长度为7～15个像素点。

4.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤5a)中适用于强度SAR图像的条件概率密度函数为：

$p\left(v_t\right|u_s)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}\ln v_{t,k}))$

其中，v_t为中心像素块，u_s为相似块，

为平滑参数，v_t，k为v_t中逐行扫描的第k个像素点，u_s，k为u_s中逐行扫描的第k个像素点，L为SAR图像视数。

5.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤5a)中适用于幅度SAR图像的条件概率密度函数为：

$p\left(v_t\right|u_s)=\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}\ln v_{t,k}))$

其中，v_t为中心像素块，u_s为相似块，

为平滑参数，v_t，k为v_t中逐行扫描的第k个像素点，u_s，k为u_s中逐行扫描的第k个像素点，L为SAR图像视数。

6.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤5c)中适用于强度SAR图像的相似性权值公式为：

$w(t,s)=p\left(v_t\right|u_s)\×g(\·)$

$=\frac{1}{Z\left(t\right)}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

其中，Z(t)为归一化因子：

$Z\left(t\right)=\underset{s=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}}{u_{s,k}}+\ln u_{s,k}-\frac{L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

7.根据权利要求1所述的一种基于MRF先验的SAR图像去斑方法，其特征在于：所述步骤5c)中适用于幅度SAR图像的相似性权值公式为：

$w(t,s)=p\left(v_t\right|u_s)\×g(\·)$

$=\frac{1}{Z\left(t\right)}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

其中，Z(t)为归一化因子：

$Z\left(t\right)=\underset{s=1}{\overset{N\×N}{\mathrm{\Σ}}}\exp (-\underset{k=1}{\overset{M\×M}{\mathrm{\Σ}}}(\frac{1}{\tilde{h}}(\frac{v_{t,k}^2}{u_{s,k}}+\ln u_{s,k}-\frac{2L-1}{L}{\ln v}_{t,k})+\frac{1}{\mathrm{\β}}|\underset{u_i\∈{\mathrm{\η}}_{s,k}}{\mathrm{\Σ}}{(u_{s,k}-u_i)}^2-\underset{u_j\∈{\mathrm{\η}}_{t,k}}{\mathrm{\Σ}}{(u_{t,k}-u_j)}^2\left|\right))$

Images