CN103236046A - 基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法 - Google Patents

Abstract

本发明公开了一种基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法。该方法通过分数阶全变差正则化加性噪声去噪和残差图像加权反馈两个步骤的交替迭代来实现相干斑滤波，利用对噪声标准差、图像卡通形态成分及其相应的残差图像的局部方差的估计，计算每个像素点归属于图像边缘、纹理和平滑三种形态的模糊隶属度，在此基础上提出模型参数的自适应计算方法，并简化分数阶差分的计算，提出一种分数阶自适应相干斑滤波方法。本方法能有效抑制噪声和“阶梯效应”，较好地保持图像边缘和纹理细节，滤波图像具有良好的视觉效果。本方法计算速度快，算法参数自适应计算，具有良好的实用性，在遥感、合成孔径雷达以及医学成像等领域具有广泛应用前景。

Description

基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法

技术领域

本发明属于图像处理领域，具体地说属于图像乘性噪声（相干斑）抑制的滤波技术，特别是一种基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法。

背景技术

在遥感、合成孔径雷达、核磁共振等相干成像过程中，图像不可避免地受到乘性噪声（相干斑）的污染。相干斑噪声严重降低了图像质量，影响了图像的判读、分类以及进一步处理，因此相干斑噪声抑制是许多图像后处理前的必要步骤。在抑制图像相干斑噪声的过程中，保持图像的边缘、纹理等重要几何结构是非常重要的。同时，为适应实际应用中对不同的图像的处理，需要对去噪模型和算法的参数进行自适应计算，并对算法的稳定性及计算效率提出较高的要求。因此研究能够良好保持图像边缘、纹理等结构形态的自适应相干斑滤波方法，具有重要实用价值和广泛的应用前景。

目前，国际上关于相干斑抑制过程中的图像细节保持问题已有许多研究，其中基于全变差正则化的方法在边缘保持方面具有较好的效果，Aubert G.等人针对服从Gamma分布的乘性噪声所提出的AA模型（Aubert G.and Aujol F.A variational approach to removing multiplicative noise.SIAM Journal on Applied Mathematics,2008,68(4):925-946）是这类方法的代表，也是近年来使用较多的乘性噪声抑制模型。但该模型对于图像的纹理结构保持并不理想，而且全变差正则化容易导致去噪图像出现“阶梯效应”。2010年，张军等人在AA模型基础上提出了基于分数阶全变差正则化的分数阶AA模型（张军，韦志辉.SAR图像去噪的分数阶多尺度变分PDE模型及自适应算法.电子与信息学报,2010,32(7):1654-1659），这种方法由于采用了分数阶全变差正则化，在纹理保持和“阶梯效应”抑制方面有比较好的效果。

针对模型计算问题，虽然Huang Y.等人针对AA模型提出的HNW算法（Huang Y.,Ng M.and Wen Y.A new total variation method for multiplicative noise removal.SIAM Journal onImaging Sciences,2009,2(1):22–40）从一定程度上降低了模型的计算复杂度，但是它在每一个迭代步中都需要求解两个特定的非线性方程，算法的收敛速度慢，计算量比较大，而且模型的参数需要人工调节，这大大影响了该方法在实际问题中的应用。分数阶AA模型虽然在纹理保持和“阶梯效应”抑制方面有较好的效果，但由于分数阶微分计算的复杂性以及模型本身具有的强非线性不可微性，模型的计算也非常困难，目前所采用的梯度下降算法收敛速度非常慢，计算量很大。

针对加性噪声去噪问题，张军等人提出了分数阶全变差正则化加性噪声去噪方法（ZhangJ.Wei Z and Xiao L.Adaptive Fractional-order Multi-scale Method for Image Denoising.Journal ofMathematical Imaging and Vision，2012，43(1):39–49）,对目前常用的求解全变差正则化图像去噪问题的Chambolle投影算法（Chambolle A.An algorithm for total variation minimization and applications.Journal of Mathematical Imaging and Vision,2004,20(1),89-97）进行了推广，提出了一种推广Chambolle算法来求解分数阶全变差正则化加性噪声去噪问题，大大提高了计算效率。但这种方法只能解决加性噪声去噪问题，目前还没有用于相干斑乘性噪声抑制问题。

另外，虽然在张军的这两篇基于分数阶全变差图像去噪方法的文章中，对模型参数进行都进行了自适应计算，在计算参数时考虑了图像纹理、卡通（含平滑及边缘）形态成分的差异，但文章中所采用的参数估计方法受到纹理影响，容易产生较大误差，对于图像纹理形态和卡通形态成分的分类，以及分数阶差分阶数的计算都是采用简单的硬阈值方法，容易导致分类不准确，而且处理后的图像的灰度在不同图像形态成分的交界处有明显阶跃，视觉效果不够理想。

发明内容

本发明目的是提供一种基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法。

实现本发明目的的技术解决方案为：一种基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，通过以下两个步骤的交替迭代计算来实现图像相干斑噪声抑制：

步骤一：对于当前待处理的含噪声图像，利用分数阶全变差正则化加性噪声去噪方法进行去噪，得到当前的去噪图像和加性残差图像；

步骤二：将步骤一中得到的加性残差图像乘以一个权重矩阵之后，反馈加回到原始含噪声图像中，得到新的待处理的含噪声图像；

该交替迭代过程用下面的公式表示：

$\left\{\begin{array}{c}{u}^{\left(n\right)}=\arg \underset{u(i,j)>0}{\min }\{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{{\left({\mathrm{\Δ}}_x^{\mathrm{\α}(i,j)}u\right)}_{i,j}^2+{\left({\mathrm{\Δ}}_y^{\mathrm{\α}(i,j)}u\right)}_{i,j}^2}+\frac{\mathrm{\γ}}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(g_{i,j}^{\left(n\right)}-u_{i.j})}^2\},\\ g^{(n+1)}=f+w^{\left(n\right)}\·(f-u^{\left(n\right)}),\end{array}\right.---\left(1\right)$

其中，n＝1,2,...，公式（1）中各符号含义为：

f：大小为M×N的原始含噪声图像；

g⁽ⁿ⁾：第n步迭代中的待处理的大小为M×N的含噪声图像，其初始值为g⁽¹⁾＝f；

u⁽ⁿ⁾：第n步迭代计算得到的大小为M×N的去噪图像，f-u⁽ⁿ⁾为相应的加性残差图像；

w⁽ⁿ⁾：第n步迭代中所采用的大小为M×N的权重矩阵，为保证迭代法的收敛性，要求w⁽ⁿ⁾中的每一个元素w⁽ⁿ⁾(i,j)∈(-1,1)；

α(i,j)：图像中(i,j)处像素点的分数阶奇异性指标，该指标表示模型中在(i,j)处计算分数阶差分时的差分阶数，α＝[α(i,j)]_M×N是一个大小为M×N的矩阵，称为分数阶奇异性指标矩阵；

：图像u在像素点(i,j)处的沿垂直方向的α(i,j)阶分数阶差分；

：图像u在像素点(i,j)处的沿水平方向的α(i,j)阶分数阶差分；

γ：为一个正数，是分数阶全变差正则化加性噪声去噪模型的正则化参数。

公式（1）中的分数阶全变差正则化加性噪声去噪模型的正则化参数γ＝σ，其中σ为原始含噪声图像f中所含相干斑噪声的标准差。

图像中每个像素点(i,j)处的分数阶奇异性指标α(i,j)按照如下方式进行计算：

$\mathrm{\α}(i,j)=T_{\mathrm{area}}(i,j)\·[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\max }-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}\·1.6+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\min }-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}\·1.4]$      (2)

$+C_{\mathrm{area}}(i,j)\·[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\max }-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}\·1.4+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\min }-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}\·1.2]+E_{\mathrm{area}}(i,j),$

公式（2）中各符号的含义为：

T_area(i,j)：图像像素点(i,j)的纹理形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像纹理形态成分的可能性的度量；

C_area(i,j)：图像像素点(i,j)的平滑形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像平滑形态成分的可能性的度量；

E_area(i,j)：图像像素点(i,j)的边缘形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像边缘形态成分的可能性的度量；

Var_loc：对乘性残差图像v^C＝f/u^C的局部方差的估计，其中f为原始含噪声图像，u^C是对图像f中的卡通形态成分，这里的卡通形态成分包含图像的平滑形态成分与边缘形态成分；

${\mathrm{Var}}_{\mathrm{loc}}^{T\min }:\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\≠0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{T\max }:\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\≠0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{C\min }:\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\≠0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{C\max }:\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\≠0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\}.$

所述的相干斑噪声的标准差σ、图像卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差Var_loc，按照下面的方法进行计算：

首先，按照如下公式对图像所含相干斑噪声的标准差进行初步估计：

$\widetilde{\mathrm{\σ}}=\frac{{\mathrm{Mid}}_{\mathrm{HH}}}{0.6745}---\left(3\right)$

其中，Mid_HH为图像v_soft＝f/exp(W_soft(logf))经过小波分解后得到的最高频HH子带小波系数的幅度中值，这里W_soft(·)表示小波软阈值运算；

其次，固定参数

α(i,j)＝1，1≤i≤M,1≤j≤N，按照公式（1）所述迭代格式进行迭代计算，每当得到新的去噪图像u⁽ⁿ⁾，权重矩阵w⁽ⁿ⁾的每个元素按照下面公式进行计算：

$w_{i,j}^{\left(n\right)}=\frac{{\mathrm{\λ}}_{i,j}^{\left(n\right)}-{\mathrm{\λ}}_{\min }^{\left(n\right)}}{{\mathrm{\λ}}_{\max }^{\left(n\right)}-{\mathrm{\λ}}_{\min }^n}\·0.1+\frac{{\mathrm{\λ}}_{i,j}^{\left(n\right)}-{\mathrm{\λ}}_{\max }^{\left(n\right)}}{{\mathrm{\λ}}_{\min }^{\left(n\right)}-{\mathrm{\λ}}_{\max }^{\left(n\right)}}\·0.01-1---\left(4\right)$

其中， ${\mathrm{\λ}}_{i,j}^{\left(n\right)}=\frac{1}{\mathrm{MN}{\widetilde{\mathrm{\σ}}}^2}{\left[u_{i,j}^{\left(n\right)}\right]}^{-2}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}\left|\right[g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)}\left]\right[f_{k,l}-u_{k,l}^{\left(n\right)}\left]\right|$ ，1≤i≤M,1≤j≤N， ${\mathrm{\λ}}_{\min }^{\left(n\right)}=\min \{{\mathrm{\λ}}_{i,j}^{\left(n\right)},1\≤i\≤M,1\≤j\≤N\},$ ${\mathrm{\λ}}_{\max }^{\left(n\right)}=\max \{{\mathrm{\λ}}_{i,j}^{\left(n\right)},1\≤i\≤M,1\≤j\≤N\};$

在这个迭代过程中，每当得到新的去噪图像u⁽ⁿ⁾时，计算相应的乘性残差图像v⁽ⁿ⁾＝f/u⁽ⁿ⁾的局部方差矩阵

，该矩阵每个元素采用如下公式进行计算：

${\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}(i,j)=\frac{1}{K^2}\underset{(p,q)\∈W_{i,j}}{\mathrm{\Σ}}{[v^{\left(n\right)}(p,q)-\frac{1}{K^2}\underset{(p,q)\∈W_{i,j}}{\mathrm{\Σ}}{v}^{\left(n\right)}(p,q)]}^2---\left(5\right)$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]是一个以(i,j)为中心的大小为K×K的窗口，K为奇数；

在得到局部方差矩阵

后，计算该矩阵所有元素平均值，并记为

；然后计算该矩阵中满足

的所有元素的平均值，并记为

；当满足

和

两个条件之一时，迭代终止，并得到最终估计的相干斑噪声的标准差σ、图像f中的卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差矩阵Var_loc如下

$\left\{\begin{array}{c}\mathrm{\σ}=\min \{\widetilde{\mathrm{\σ}},\frac{1}{2}(M_{\mathrm{lv}-\mathrm{low}}^{\left(n\right)}+M_{\mathrm{lv}-\mathrm{low}}^{(n-1)})\},\\ u^C=u^{\left(n\right)},\\ {\mathrm{Var}}_{\mathrm{loc}}={\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}.\end{array}\right.---\left(6\right)$

图像中(i,j)处像素点的边缘形态模糊隶属度E_area(i,j)按照下面步骤进行计算：

首先，针对依据权利要求4计算得到的卡通形态成分u^C，使用Canny边缘检测算子进行边缘检测，得到二值图像E(u^C)，在边缘处值为1，其余部分值为0；

然后，将二值图像E(u^C)与一个标准化Gauss模板G进行卷积，得到图像 ${\mathrm{GE}}_{u^C}=E\left(u^C\right)*G;$

最后，计算图像中(i,j)处像素点的边缘形态模糊隶属度如下：

$E_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}1,& {\mathrm{GE}}_{u^C}(i,j)>0,\\ 0,& {\mathrm{GE}}_{u^C}(i,j)=0.\end{array}\right.---\left(7\right)$

图像中(i,j)处像素点的纹理形态模糊隶属度T_area(i,j)和平滑形态模糊隶属度C_area(i,j)按照下面步骤进行计算：

首先，利用依据权利要求4计算得到的局部方差矩阵Var_loc，对图像中(i,j)处像素点的纹理形态隶属度进行初步估计

${\tilde{T}}_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}1,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)\&GreaterEqual;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\\ 0,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)<\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\end{array}\right.---\left(8\right)$

其中mean(Var_loc)为局部方差矩阵Var_loc中所有元素的平均值；

然后，计算图像中的每个像素点(i,j)的纹理形态模糊隶属度T_area(i,j)和平滑形态模糊隶属度C_area(i,j)，其计算公式为：

$\left\{\begin{array}{c}{T}_{\mathrm{area}}(i,j)=\frac{1}{K^2}(1-E_{\mathrm{area}}(i,j))\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{\tilde{T}}_{\mathrm{area}}(p,q),\\ C_{\mathrm{area}}(i,j)=1-E_{\mathrm{area}}(i,j)-T_{\mathrm{area}}(i,j).\end{array}\right.---\left(9\right)$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]为以(i,j)为中心的大小为K×K的窗口，K为奇数，E_area(i,j)为依据权利要求5计算得到的像素点(i,j)处的边缘形态模糊隶属度。

对于任意一幅图像u，在图像像素点(i,j)处的α(i,j)阶分数阶差分和按照下面的公式计算：

${\left({\mathrm{\&Delta;}}_x^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i-k,j),{\left({\mathrm{\&Delta;}}_y^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i,j-k)---\left(10\right)$

其中，k＝0,1,...,L为L+1个分数阶差分系数，这里

k＝0,1,...,L-1，

，L和L_max为满足2≤L＜L_max≤min{M,N}的正数，

为广义二项式系数,Γ(x)为Gamma函数，且当α(i,j)≤k-1时， $C_k^{\mathrm{\&alpha;}(i,j)}=0.$

在按照迭代公式（1）进行计算过程中，当得到新的去噪图像u⁽ⁿ⁾时，按照下面方式对加权矩阵w⁽ⁿ⁾进行更新：

$w_{i,j}^{\left(n\right)}=(1-C_{\mathrm{area}}(i,j))\&CenterDot;[\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)})C_3}{{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}}+\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)})C_2}{{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}}]+C_{\mathrm{area}}(i,j)\&CenterDot;C_1-1---\left(11\right)$

其中 $\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}=\underset{(k,l)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}\left|(g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)})(f_{k,l}-u_{k,l}^{\left(n\right)})\right|/\left[K^2{\mathrm{\&sigma;}}^2{\left(u_{i,j}^{\left(n\right)}\right)}^2\right]\\ {\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}=\min \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ {\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}=\max \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right.$ ，这里σ是噪声整体标准差，W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]为以(i,j)为中心的大小为K×K的窗口，K为奇数，C₁,C₂,C₃为3个参数。

参数C₁,C₂,C₃满足0＜C₁＜C₂≤C₃＜2，其计算公式为

$\left\{\begin{array}{c}{C}_1=1-e^{-0.01/\left({2\mathrm{\&sigma;}}^2\right)},\\ C_2=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)-C_1\},\\ C_3=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)+C_1\},\end{array}\right.---\left(12\right)$

其中σ和Var_loc是按照权利4计算得到的相干斑噪声标准差和局部方差矩阵， $\left\{\begin{array}{c}\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{1,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ \mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right..$

本发明与现有技术相比，其显著优点：（1）本发明方法能够自适应地实现服从Gamma分布的相干斑噪声抑制。（2）本发明方法提出的迭代法所得到迭代序列可以快速并稳定地收敛，运算时间较少。（3）可以在有效抑制图像噪声的同时，较好保持图像的边缘、纹理等细节信息，并有效地抑制“阶梯效应”，得到具有良好视觉效果的去噪图像。本发明方法在光学遥感、合成孔径雷达以及医学成像等都有广泛的应用前景。

附图说明

图1是本发明方法整体算法结构及数据流程图。

图2是相干斑滤波预处理单元算法及数据流程图。

图3是图像形态模糊隶属度计算单元算法及数据流程图。

图4是分数阶差分运算预处理单元算法及数据流程图。

图5是自适应相干斑滤波单元算法及数据流程图。

图6是本发明的实验测试图像。

图7是本发明方法与HNW方法对含有方差为σ²＝0.03的Gamma乘性噪声的Barbara图像去噪得到的去噪图像、去噪图像平滑及纹理区域比较。

图8是本发明方法与HNW方法对含有方差为σ²＝0.03的Gamma乘性噪声的Lena图像去噪得到的去噪图像、去噪图像部分平滑及纹理区域比较。

图1至图5给出了本发明方法的整体算法结构以及具体计算单元的算法及数据流程图。图6为实施例的试验图像，图7至图8为本发明方法与其他方法的相干斑滤波图像及局部放大图的比较。表1给出了本发明方法与其他方法在噪声抑制、结构保持以及运算效率方面的客观指标比较。

具体实施方式

本发明方法考虑服从Gamma分布的相干斑乘性噪声滤波问题，通过分数阶全变差正则化加性噪声去噪和残差图像加权反馈两个步骤的交替迭代来实现相干斑滤波，利用对噪声标准差、图像卡通形态成分及其相应的残差图像的局部方差的估计，计算每个像素点归属于图像边缘、纹理和平滑三种形态的模糊隶属度，在此基础上提出模型参数的自适应计算方法，并简化分数阶差分的计算，提出一种分数阶自适应相干斑滤波方法。

本发明方法所指的相干斑是指服从Gamma分布的乘性噪声。利用算子分裂技术，将相干斑噪声抑制分解为现有的分数阶全变差加性噪声抑制和残差图像加权反馈两个步骤，然后对这两个步骤进行交替迭代来实现相干斑抑制，其交替迭代格式为：

$\left\{\begin{array}{c}{u}^{\left(n\right)}=\arg \underset{u(i,j)>0}{\min }\{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\sqrt{{\left({\mathrm{\&Delta;}}_x^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}^2+{\left({\mathrm{\&Delta;}}_y^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}^2}+\frac{\mathrm{\&gamma;}}{2}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(g_{i,j}^{\left(n\right)}-u_{i.j})}^2\},\\ g^{(n+1)}=f+w^{\left(n\right)}\&CenterDot;(f-u^{\left(n\right)}),\end{array}\right.---\left(1\right)$

其中f是相干斑噪声污染图像，迭代序列{u⁽ⁿ⁾}收敛的结果就是最终的滤波图像。

该迭代格式中的第一个等式是以g⁽ⁿ⁾为待处理的含噪声图像的分数阶全变差正则化加性噪声抑制变分模型，

分别表示图像u在像素点(i,j)处的沿垂直方向和水平方向的α(i,j)阶分数阶差分，这里α(i,j)称为分数阶奇异性指标，矩阵α＝[α(i,j)]_M×N称为分数阶奇异性指标矩阵。特别地，当α(i,j)≡1时，第一个等式就是经典的全变差正则化加性噪声抑制变分模型。参数γ＞0为正则化参数，在本发明中，取γ＝σ，其中σ为相干斑噪声的标准差。

在该迭代格式中的第二个等式中，将f-u⁽ⁿ⁾看作是每次迭代得到的加性残差图像，通过乘以一个权重矩阵w⁽ⁿ⁾，将加权后的残差图像反馈到原含噪声图像中去，得到新的待处理的含噪声图像g⁽ⁿ⁺¹⁾，这一步称为残差图像加权反馈。理论分析表明，该迭代格式所得到的迭代序列{u⁽ⁿ⁾}收敛的一个充分条件是权重矩阵w⁽ⁿ⁾中的每个元素都满足条件-1＜w⁽ⁿ⁾(i,j)＜1，在本发明中我们将对w⁽ⁿ⁾的每一个元素加上这一限制性条件。

迭代格式中的第一个等式所表示的分数阶全变差正则化加性噪声去噪变分模型，可以用在本发明技术背景中所提的推广Chambolle投影算法进行求解，而第二个等式是可以直接计算的。因此，本发明的核心内容是针对迭代格式中的参数，包括乘性噪声整体标准差σ、分数阶奇异性指标矩阵α以及权重矩阵w⁽ⁿ⁾提出自适应计算方法。一方面，参数自适应使得算法可以适用于不同的图像，提高算法适用性；另一方面，通过参数的自适应选择，对于图像的边缘、纹理和平滑等不同形态成分进行区分和处理，在有效地去除相干斑噪声的同时，保持图像纹理细节并有效抑制“阶梯效应”，得到视觉效果良好的滤波图像。另外，本发明还对模型求解中所涉及的分数阶差分的近似计算提出改进方法，进一步减少算法的计算量。

本发明所提供的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法由四个处理单元构成，分别是：相干斑滤波预处理单元、图像形态模糊隶属度计算单元、分数阶差分运算预处理单元和自适应相干斑滤波计算单元。

1.1相干斑滤波预处理单元

在实际处理中，我们得到的是一幅已经受到相干斑噪声污染的图像，但是噪声的标准差等参数却是未知的，需要通过计算来进行估计。另外，为更好地在图像的平滑部分抑制噪声的同时，保持图像的边缘、纹理等形态成分，也需要对这些不同形态的图像成分进行区分。

本单元的主要功能就是估计相干斑噪声的整体标准差σ、图像的卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差矩阵Var_loc。本单元所计算得到的这些参数是为本发明中的相干斑滤波中的参数自适应计算所做的预备工作，因此本单元称为相干斑滤波预处理单元。

所述相干斑滤波预处理单元，由下列计算步骤构成：

步骤1.输入一幅大小为M×N的待处理的原始含噪声图像f；

步骤2.按照如下公式对噪声的整体标准差进行初步估计：

$\widetilde{\mathrm{\&sigma;}}=\frac{{\mathrm{Mid}}_{\mathrm{HH}}}{0.6745}$

其中，Mid_HH为图像v_soft＝f/exp(W_soft(logf))通过小波分解后得到的最高频HH子带小波系数的幅度中值，这里W_soft(·)表示小波软阈值运算。

步骤3.设定初值u⁽⁰⁾＝f，g⁽¹⁾＝f，对n＝1,2,...,进行如下迭代：

步骤4.对待处理的含噪声图像g⁽ⁿ⁾，取定

，α(i,j)＝1，1≤i≤M,1≤j≤N，用现有的分数阶全变差加性噪声去噪方法进行去噪，得到新的去噪图像u⁽ⁿ⁾；

步骤5.计算乘性残差图像v⁽ⁿ⁾＝f/u⁽ⁿ⁾,然后对于残差图像v⁽ⁿ⁾中的每个像素点(i,j),计算(i,j)处的局部方差，其计算公式为:

${\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}(i,j)=\frac{1}{K^2}\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{[v^{\left(n\right)}(p,q)-\frac{1}{K^2}\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{v}^{\left(n\right)}(p,q)]}^2$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]为以(i,j)为中心的大小为K×K(7≤K≤15,K为奇数)的窗口。

步骤6.计算

的平均值，并记为

计算满足

的所有像素点的局部方差的平均值，并记为

步骤7.当满足两个迭代终止条件之一时，迭代终止，输出最终估计的相干斑噪声的标准差σ、图像f中的卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差矩阵Var_loc如下

$\left\{\begin{array}{c}\mathrm{\&sigma;}=\min \{\widetilde{\mathrm{\&sigma;}},\frac{1}{2}(M_{\mathrm{lv}-\mathrm{low}}^{\left(n\right)}+M_{\mathrm{lv}-\mathrm{low}}^{(n-1)})\},\\ u^C=u^{\left(n\right)},\\ {\mathrm{Var}}_{\mathrm{loc}}={\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}.\end{array}\right.$

若不满足上述迭代终止条件的任何一个，则转入步骤8；

步骤8.计算矩阵λ⁽ⁿ⁾，它在每个像素点(i,j)处的计算公式为：

${\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}=\frac{1}{\mathrm{MN}{\widetilde{\mathrm{\&sigma;}}}^2}{\left[u_{i,j}^{\left(n\right)}\right]}^{-2}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right[g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)}\left]\right[f_{k,l}-u_{k,l}^{\left(n\right)}\left]\right|,$

步骤9.计算权重矩阵w⁽ⁿ⁾，其每个元素按如下公式计算：

$w_{i,j}^{\left(n\right)}=\frac{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\min }^{\left(n\right)}}{{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\min }^n}\&CenterDot;0.1+\frac{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}}{{\mathrm{\&lambda;}}_{\min }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}}\&CenterDot;0.01-1$

其中 ${\mathrm{\&lambda;}}_{\min }^{\left(n\right)}=\min \{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)},1\&le;i\&le;M,1\&le;j\&le;N\},$ ${\mathrm{\&lambda;}}_{\max }^{\left(n\right)}=\max \{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)},1\&le;i\&le;M,1\&le;j\&le;N\}.$

步骤10.计算新的待处理的含噪声图像

g⁽ⁿ⁺¹⁾＝f+w⁽ⁿ⁾·(f-u⁽ⁿ⁾),

令n:＝n+1，转步骤4进行下一轮迭代。

1.2图像形态模糊隶属度计算单元

本发明所指的图像形态主要是指图像中的边缘、纹理以及平滑三种形态成分，这里分别用E_area(i,j)、T_area(i,j)和C_area(i,j)表示图像中(i,j)处像素点属于边缘、纹理和平滑形态的可能性（即隶属度），为了使各种形态隶属度之间有一个相对平缓的过渡，本发明计算三种形态成分的模糊隶属度。

所述的图像形态模糊隶属度计算单元，由下列计算步骤构成：

步骤1.对于1.1所述的相干斑滤波预处理单元得到的卡通形态成分u^C，首先使用Canny边缘检测算子进行边缘检测，得到关于边缘的二值图像E(u^C)（在边缘处取值为1，其余部分取值为0）；

步骤2.将二值图像E(u^C)与一个窗口大小为3×3进行标准化Gauss模板G_3×3进行卷积，得到图像 ${\mathrm{GE}}_{u^C}=E\left(u^C\right)*G_{3\&times;3}$

步骤3.计算图像中(i,j)处像素点的边缘形态模糊隶属度如下：

$E_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}{\mathrm{GE}}_{u^C}(i,j),& {\mathrm{GE}}_{u^C}(i,j)>0\\ 0,& {\mathrm{GE}}_{u^C}=0\end{array}\right.$

步骤4.利用1.1所述的相干斑滤波预处理单元得到的局部方差矩阵Var_loc，对图像中(i,j)处像素点的纹理形态隶属度进行初步估计

${\tilde{T}}_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}1,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)\&GreaterEqual;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\\ 0,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)<\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\end{array}\right.$

其中mean(Var_loc)为局部方差矩阵Var_loc中元素的平均值。

步骤5.计算图像中的每个像素点(i,j)的纹理形态模糊隶属度T_area(i,j)和平滑形态模糊隶属度C_area(i,j)，其计算公式为:

$\left\{\begin{array}{c}{T}_{\mathrm{area}}(i,j)=\frac{1}{K^2}(1-E_{\mathrm{area}}(i,j))\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{\tilde{T}}_{\mathrm{area}}(p,q),\\ C_{\mathrm{area}}(i,j)=1-E_{\mathrm{area}}(i,j)-T_{\mathrm{area}}(i,j).\end{array}\right.$

其中W_i,j为以(i,j)为中心的大小为K×K的窗口，与相干斑滤波预处理单元的步骤5所采用的窗口相同。

1.3分数阶差分运算预处理单元

本发明利用推广的Chambolle算法求解分数阶全变差正则化加性噪声去噪模型时，在不同像素点处采用不同的分数阶差分阶数，即取不同的分数阶奇异性指标α(i,j)。另外，在每次迭代中都需要计算分数阶差分

和

，其计算公式为：

${\left({\mathrm{\&Delta;}}_x^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i-k,j),{\left({\mathrm{\&Delta;}}_y^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i,j-k)$

其中，，k＝0,1,2,...,L为分数阶差分系数，为提高计算精度，一般取L≥20。但在本发明中，我们对该计算格式进行改进，只需要取L＝4即可，大大节省了计算量。

本单元的主要功能就是计算各个像素点的分数阶奇异性指标和分数阶差分系数，这些参数计算之后就存储起来，以方便在后面的分数阶差分运算中直接调用。

所述的分数阶差分运算预处理单元，由下列计算步骤构成：

步骤1.计算分数阶奇异性指标矩阵α，在每一点(i,j)处的分数阶奇异性指标α(i,j)按如下公式计算：

$\mathrm{\&alpha;}(i,j)=T_{\mathrm{area}}(i,j)\&CenterDot;[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\max }-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}\&CenterDot;1.6+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\min }-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}\&CenterDot;1.4]$

$+C_{\mathrm{area}}(i,j)\&CenterDot;[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\max }-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}\&CenterDot;1.4+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\min }-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}\&CenterDot;1.2]+E_{\mathrm{area}}(i,j),$

其中， ${\mathrm{Var}}_{\mathrm{loc}}^{T\min }=\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\},$

${\mathrm{Var}}_{\mathrm{loc}}^{T\max }=\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\},$

${\mathrm{Var}}_{\mathrm{loc}}^{C\min }=\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{C\max }=\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\}.$

步骤2.计算并存储五个分数阶差分系数矩阵

，k＝0,1,2,3,4，其每个元素按如下公式计算：

$\left\{\begin{array}{c}{w}_k^{\mathrm{\&alpha;}(i,j)}={(-1)}^k{C}_k^{\mathrm{\&alpha;}(i,j)},k=\mathrm{0,1,2,3}\\ w_4^{\mathrm{\&alpha;}(i,j)}={\mathrm{\&Sigma;}}_{k=4}^{\min (M,N,50)}{(-1)}^k{C}_k^{\mathrm{\&alpha;}(i,j)}\end{array}\right.,$

其中

为广义二项式系数,Γ(x)为Gamma函数，且当α(i,j)≤k-1时， $C_k^{\mathrm{\&alpha;}(i,j)}=0.$

1.4自适应相干斑滤波单元

本单元是具体执行相干斑滤波的单元，前面1.1—1.3所述计算单元都是为本单元的运算计算相关参数。

所述的自适应相干斑滤波单元，按照下述步骤来实现：

步骤1.输入一幅大小为M×N的待去噪图像f，输入利用1.1—1.3单元计算的导的噪声整体标准差σ和局部方差矩阵Var_loc，图像形态隶属度矩阵E_area，T_area和C_area，分数阶奇异性指标矩阵α和分数阶差分系数矩阵矩阵

步骤2.按照如下公式计算参数C₁,C₂和C₃：

$\left\{\begin{array}{c}{C}_1=1-e^{-0.01/\left({2\mathrm{\&sigma;}}^2\right)},\\ C_2=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)-C_1\},\\ C_3=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)+C_1\},\end{array}\right.$

其中 $\left\{\begin{array}{c}\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{1,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ \mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right..$

步骤3.设定初值u⁽⁰⁾＝f，g⁽¹⁾＝f，进行如下迭代：

步骤4.对待处理的去噪图像g⁽ⁿ⁾，将噪声整体标准差σ作为正则化参数，将分数阶奇异性指标矩阵α和分数阶差分系数矩阵矩阵

作为输入，利用推广的Chambolle算法对g⁽ⁿ⁾进行分数阶全变差正则化加性噪声去噪，得到新的去噪图像u⁽ⁿ⁾，对于给定的误差限ε＞0，若满足迭代终止条件

$\frac{{\left|\right|u^{\left(n\right)}-u^{(n-1)}\left|\right|}_2}{{\left|\right|u^{(n-1)}\left|\right|}_2}<\mathrm{\&epsiv;},$

则迭代终止，并输出最终的相干斑滤波图像u⁽ⁿ⁾；否则，进入步骤5；

步骤5.计算矩阵λ⁽ⁿ⁾，该矩阵第(i,j)元素的计算公式

${\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}=\frac{1}{K^2{\mathrm{\&sigma;}}^2{\left(u_{i,j}^{\left(n\right)}\right)}^2}\underset{(k,l)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}\left|(g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)})(f_{k,l}-u_{k,l}^{\left(n\right)})\right|,$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]是以(i,j)为中心的大小为K×K(3≤K≤9,K为奇数)的窗口；

步骤6.对加权矩阵w⁽ⁿ⁾进行更新，其每个元素

按照如下公式计算：

$w_{i,j}^{\left(n\right)}=(1-C_{\mathrm{area}}(i,j))\&CenterDot;[\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)})C_3}{{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}}+\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)})C_2}{{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}}]+C_{\mathrm{area}}(i,j)\&CenterDot;C_1-1,$

其中 $\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}=\min \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ {\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}=\max \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right..$

步骤7.计算新的待处理的含噪声图像

g⁽ⁿ⁺¹⁾＝f+w⁽ⁿ⁾·(f-u⁽ⁿ⁾),

令n:＝n+1，转步骤4进行下一轮迭代。

下面结合实施例对本发明做进一步说明。

1)实验条件：

实验所用的计算环境为Intel i7-2620M2.7GHz双核CPU，内存为8GB1333MHz，编程平台为Matlab R2011a。实验所用的测试图像为大小为512×512的国际标准测试图像Barbara图像和Lena图像。

2)实验内容：

我们首先在原始Barbara图像和Lena图像上分别加入均值为1，方差分别为0.03和0.1的服从Gamma分布相干斑噪声，然后利用HNW算法和本发明提出的方法来进行滤波。在使用本发明方法时，分为两种方式：

（a）不使用1.3所述的分数阶差分运算预处理单元中计算的分数阶奇异性指标矩阵，而是直接另α＝1（其中1为大小为512×512的全1矩阵），此时本发明方法和HNW算法求解的都是完全相同的模型，即AA模型；

（b）按照1.3所述的分数阶差分运算预处理单元计算分数阶奇异性指标矩阵α，此时求解的是分数阶全变差正则化模型。

本实验中从三个指标对算法进行比较：

（1）SSIM，即结构相似度，该指标主要衡量去噪方法对于图像边缘、纹理等几何结构的保持性能；

（2）PSNR，即峰值信噪比，该指标主要用于衡量算法的去噪性能；

（3）CPU运算时间，该指标主要用于衡量算法的计算速度和计算效率。

本实验中所有算法的迭代终止条件都是前后两次迭代结果的相对误差小于10^-4。本发明方法是自适应方法，算法中的参数都是自适应计算的，不需要人工干预，而HNW算法并不是自适应算法，在本实验中，我们调节HNW算法参数，使得当满足迭代终止条件时，去噪图像能够达到最高的结构相似度（SSIM）。

表1中给出了利用几种方法的性能指标比较。为了更好地比较运算时间，所给出了CPU时间是20次计算的平均结果（20次计算中，含噪声图像和参数都没有变化，因此SSIM和PSNR是一样的，但CPU时间会略有波动）。

表1本发明方法与HNW算法在SSIM、PSNR和CPU运算时间的比较

实验结果表明：在对AA模型的求解中，本发明方法在SSIM和PSNR改善方面，比HNW算法有所提高，而在CPU运算时间方面，本发明方法有非常明显的优势，运算时间只有HNM算法的1/5～1/6，运算效率大大提高，而且由于本发明方法的参数都是自适应计算，因此本发明方法的比HNM算法具有更好的实用性。当算法中的分数阶奇异性指数也自适应计算时，本发明方法在运算时间上也要少于HNW算法，而且在SSIM和PSNR方面有比较明显的改善。

图7-图8给出了三种方法对于含有σ²＝0.03的Gamma乘性噪声的实验图像的去噪图像及其局部图像的比较。图7和图8的第一列是去噪图像，第二列给出了平滑区域的局部图像，第三列给出了纹理区域的局部图像。图7和图8中的第一行是HNM方法得到的结果，第二行是本发明方法取α≡1时得到的结果，第三行本发明方法当分数阶奇异性指标矩阵自适应计算时的得到的结果。

实验结果表明，HNM方法在取得最高SSIM时，在平滑区域的去噪效果并不十分理想，而且容易出现“阶梯效应”，参见图7（a2）以及图8（a2）。本发明的方法在可以在较好抑制图像平滑区域噪声的同时，有效保持图像纹理成分，参见图7(b3)、图7(c3)、图8(b3)以及图8(c3)。当分数阶奇异性指标矩阵α≡1时，本发明方法在图像平滑区域还是有一些“阶梯效应”，参见图7(b2)和图8(b2)；但当分数阶奇异性指标也自适应计算时，本发明方法对于“阶梯效应”的抑制是比较好的，参见图7(c2)和图8(c2)。

本实验的结果表明，本发明方法能够自适应地实现服从Gamma分布的相干斑噪声抑制，所提出的迭代法所得到迭代序列可以快速并稳定地收敛，运算时间较少，所得到的去噪图像可以在有效抑制图像噪声的同时，较好保持图像的边缘、纹理等细节信息，并有效地抑制“阶梯效应”，得到具有良好视觉效果的去噪图像。

Claims (9)

1.一种基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于，通过以下两个步骤的交替迭代计算来实现图像相干斑噪声抑制：

步骤一：对于当前待处理的含噪声图像，利用分数阶全变差正则化加性噪声去噪方法进行去噪，得到当前的去噪图像和加性残差图像；

步骤二：将步骤一中得到的加性残差图像乘以一个权重矩阵之后，反馈加回到原始含噪声图像中，得到新的待处理的含噪声图像；

该交替迭代过程用下面的公式表示：

$\left\{\begin{array}{c}{u}^{\left(n\right)}=\arg \underset{u(i,j)>0}{\min }\{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\sqrt{{\left({\mathrm{\&Delta;}}_x^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}^2+{\left({\mathrm{\&Delta;}}_y^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}^2}+\frac{\mathrm{\&gamma;}}{2}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(g_{i,j}^{\left(n\right)}-u_{i.j})}^2\},\\ g^{(n+1)}=f+w^{\left(n\right)}\&CenterDot;(f-u^{\left(n\right)}),\end{array}\right.---\left(1\right)$

其中，n＝1,2,...，公式（1）中各符号含义为：

f：大小为M×N的原始含噪声图像；

g⁽ⁿ⁾：第n步迭代中的待处理的大小为M×N的含噪声图像，其初始值为g⁽¹⁾＝f；

u⁽ⁿ⁾：第n步迭代计算得到的大小为M×N的去噪图像，f-u⁽ⁿ⁾为相应的加性残差图像；

w⁽ⁿ⁾：第n步迭代中所采用的大小为M×N的权重矩阵，为保证迭代法的收敛性，要求w⁽ⁿ⁾中的每一个元素w⁽ⁿ⁾(i,j)∈(-1,1)；

α(i,j)：图像中(i,j)处像素点的分数阶奇异性指标，该指标表示模型中在(i,j)处计算分数阶差分时的差分阶数，α＝[α(i,j)]_M×N是一个大小为M×N的矩阵，称为分数阶奇异性指标矩阵；

：图像u在像素点(i,j)处的沿垂直方向的α(i,j)阶分数阶差分；

：图像u在像素点(i,j)处的沿水平方向的α(i,j)阶分数阶差分；

γ：为一个正数，是分数阶全变差正则化加性噪声去噪模型的正则化参数。

2.根据权利要求1所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：公式（1）中的分数阶全变差正则化加性噪声去噪模型的正则化参数γ＝σ，其中σ为原始含噪声图像f中所含相干斑噪声的标准差。

3.依据权利要求1所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于，图像中每个像素点(i,j)处的分数阶奇异性指标α(i,j)按照如下方式进行计算：

$\mathrm{\&alpha;}(i,j)=T_{\mathrm{area}}(i,j)\&CenterDot;[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\max }-{\mathrm{Var}}_{\mathrm{loc}}^{T\min }}\&CenterDot;1.6+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{T\min }-{\mathrm{Var}}_{\mathrm{loc}}^{T\max }}\&CenterDot;1.4]$      (2)

$+C_{\mathrm{area}}(i,j)\&CenterDot;[\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\max }-{\mathrm{Var}}_{\mathrm{loc}}^{C\min }}\&CenterDot;1.4+\frac{{\mathrm{Var}}_{\mathrm{loc}}(i,j)-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}{{\mathrm{Var}}_{\mathrm{loc}}^{C\min }-{\mathrm{Var}}_{\mathrm{loc}}^{C\max }}\&CenterDot;1.2]+E_{\mathrm{area}}(i,j),$

公式（2）中各符号的含义为：

T_area(i,j)：图像像素点(i,j)的纹理形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像纹理形态成分的可能性的度量；

C_area(i,j)：图像像素点(i,j)的平滑形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像平滑形态成分的可能性的度量；

E_area(i,j)：图像像素点(i,j)的边缘形态模糊隶属度，该模糊隶属度是对像素点(i,j)属于图像边缘形态成分的可能性的度量；

Var_loc：对乘性残差图像v^C＝f/u^C的局部方差的估计，其中f为原始含噪声图像，u^C是对图像f中的卡通形态成分，这里的卡通形态成分包含图像的平滑形态成分与边缘形态成分；

${\mathrm{Var}}_{\mathrm{loc}}^{T\min }:\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{T\max }:\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|T_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{C\min }:\min \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\};$

${\mathrm{Var}}_{\mathrm{loc}}^{C\max }:\max \left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;0,i=\mathrm{1,2},...,M,j=\mathrm{1,2},...,N\}.$

4.依据权利要求1、2或3所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：所述的相干斑噪声的标准差σ、图像卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差Var_loc，按照下面的方法进行计算：

首先，按照如下公式对图像所含相干斑噪声的标准差进行初步估计：

$\widetilde{\mathrm{\&sigma;}}=\frac{{\mathrm{Mid}}_{\mathrm{HH}}}{0.6745}---\left(3\right)$

其中，Mid_HH为图像v_soft＝f/exp(W_soft(logf))经过小波分解后得到的最高频HH子带小波系数的幅度中值，这里W_soft(·)表示小波软阈值运算；

其次，固定参数

，α(i,j)＝1，1≤i≤M,1≤j≤N，按照公式（1）所述迭代格式进行迭代计算，每当得到新的去噪图像u⁽ⁿ⁾，权重矩阵w⁽ⁿ⁾的每个元素按照下面公式进行计算：

$w_{i,j}^{\left(n\right)}=\frac{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\min }^{\left(n\right)}}{{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\min }^n}\&CenterDot;0.1+\frac{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}}{{\mathrm{\&lambda;}}_{\min }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\max }^{\left(n\right)}}\&CenterDot;0.01-1---\left(4\right)$

其中， ${\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}=\frac{1}{\mathrm{MN}{\widetilde{\mathrm{\&sigma;}}}^2}{\left[u_{i,j}^{\left(n\right)}\right]}^{-2}\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right[g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)}\left]\right[f_{k,l}-u_{k,l}^{\left(n\right)}\left]\right|$ ，1≤i≤M,1≤j≤N， ${\mathrm{\&lambda;}}_{\min }^{\left(n\right)}=\min \{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)},1\&le;i\&le;M,1\&le;j\&le;N\},$ ${\mathrm{\&lambda;}}_{\max }^{\left(n\right)}=\max \{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)},1\&le;i\&le;M,1\&le;j\&le;N\};$

在这个迭代过程中，每当得到新的去噪图像u⁽ⁿ⁾时，计算相应的乘性残差图像v⁽ⁿ⁾＝f/u⁽ⁿ⁾的局部方差矩阵，该矩阵每个元素采用如下公式进行计算：

${\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}(i,j)=\frac{1}{K^2}\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{[v^{\left(n\right)}(p,q)-\frac{1}{K^2}\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{v}^{\left(n\right)}(p,q)]}^2---\left(5\right)$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]是一个以(i,j)为中心的大小为K×K的窗口，K为奇数；

在得到局部方差矩阵

后，计算该矩阵所有元素平均值，并记为

然后计算该矩阵中满足

的所有元素的平均值，并记为

当满足

和两个条件之一时，迭代终止，并得到最终估计的相干斑噪声的标准差σ、图像f中的卡通形态成分u^C以及相应的乘性残差图像v^C＝f/u^C的局部方差矩阵Var_loc如下

$\left\{\begin{array}{c}\mathrm{\&sigma;}=\min \{\widetilde{\mathrm{\&sigma;}},\frac{1}{2}(M_{\mathrm{lv}-\mathrm{low}}^{\left(n\right)}+M_{\mathrm{lv}-\mathrm{low}}^{(n-1)})\},\\ u^C=u^{\left(n\right)},\\ {\mathrm{Var}}_{\mathrm{loc}}={\mathrm{Var}}_{\mathrm{loc}}^{\left(n\right)}.\end{array}\right.---\left(6\right).$

5.依据权利要求1或3所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：图像中(i,j)处像素点的边缘形态模糊隶属度E_area(i,j)按照下面步骤进行计算：

首先，针对依据权利要求4计算得到的卡通形态成分u^C，使用Canny边缘检测算子进行边缘检测，得到二值图像E(u^C)，在边缘处值为1，其余部分值为0；

然后，将二值图像E(u^C)与一个标准化Gauss模板G进行卷积，得到图像 ${\mathrm{GE}}_{u^C}=E\left(u^C\right)*G;$

最后，计算图像中(i,j)处像素点的边缘形态模糊隶属度如下：

$E_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}1,& {\mathrm{GE}}_{u^C}(i,j)>0,\\ 0,& {\mathrm{GE}}_{u^C}(i,j)=0.\end{array}\right.---\left(7\right).$

6.依据权利要求1或3所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：图像中(i,j)处像素点的纹理形态模糊隶属度T_area(i,j)和平滑形态模糊隶属度C_area(i,j)按照下面步骤进行计算：

首先，利用依据权利要求4计算得到的局部方差矩阵Var_loc，对图像中(i,j)处像素点的纹理形态隶属度进行初步估计

${\tilde{T}}_{\mathrm{area}}(i,j)=\left\{\begin{array}{cc}1,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)\&GreaterEqual;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\\ 0,& {\mathrm{Var}}_{\mathrm{loc}}(i,j)<\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}\right),\end{array}\right.---\left(8\right)$

其中mean(Var_loc)为局部方差矩阵Var_loc中所有元素的平均值；

然后，计算图像中的每个像素点(i,j)的纹理形态模糊隶属度T_area(i,j)和平滑形态模糊隶属度C_area(i,j)，其计算公式为：

$\left\{\begin{array}{c}{T}_{\mathrm{area}}(i,j)=\frac{1}{K^2}(1-E_{\mathrm{area}}(i,j))\underset{(p,q)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}{\tilde{T}}_{\mathrm{area}}(p,q),\\ C_{\mathrm{area}}(i,j)=1-E_{\mathrm{area}}(i,j)-T_{\mathrm{area}}(i,j).\end{array}\right.---\left(9\right)$

其中W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]为以(i,j)为中心的大小为K×K的窗口，K为奇数，E_area(i,j)为依据权利要求5计算得到的像素点(i,j)处的边缘形态模糊隶属度。

7.依据权利要求1所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：对于任意一幅图像u，在图像像素点(i,j)处的α(i,j)阶分数阶差分

和

按照下面的公式计算：

${\left({\mathrm{\&Delta;}}_x^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i-k,j),{\left({\mathrm{\&Delta;}}_y^{\mathrm{\&alpha;}(i,j)}u\right)}_{i,j}=\underset{k=0}{\overset{L}{\mathrm{\&Sigma;}}}{w}_k^{\mathrm{\&alpha;}(i,j)}u(i,j-k)---\left(10\right)$

其中，

,k＝0,1,...,L为L+1个分数阶差分系数，这里

k＝0,1,...,L-1，

，L和L_max为满足2≤L＜L_max≤min{M,N}的正数，

为广义二项式系数,Γ(x)为Gamma函数，且当α(i,j)≤k-1时，

8.依据权利要求1所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：在按照迭代公式（1）进行计算过程中，当得到新的去噪图像u⁽ⁿ⁾时，按照下面方式对加权矩阵w⁽ⁿ⁾进行更新：

$w_{i,j}^{\left(n\right)}=(1-C_{\mathrm{area}}(i,j))\&CenterDot;[\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)})C_3}{{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}}+\frac{({\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)})C_2}{{\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}-{\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}}]+C_{\mathrm{area}}(i,j)\&CenterDot;C_1-1---\left(11\right)$

其中 $\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}=\underset{(k,l)\&Element;W_{i,j}}{\mathrm{\&Sigma;}}\left|(g_{k,l}^{\left(n\right)}-u_{k,l}^{\left(n\right)})(f_{k,l}-u_{k,l}^{\left(n\right)})\right|/\left[K^2{\mathrm{\&sigma;}}^2{\left(u_{i,j}^{\left(n\right)}\right)}^2\right]\\ {\mathrm{\&lambda;}}_{\mathrm{ET}\min }^{\left(n\right)}=\min \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ {\mathrm{\&lambda;}}_{\mathrm{ET}\max }^{\left(n\right)}=\max \left\{{\mathrm{\&lambda;}}_{i,j}^{\left(n\right)}\right|C_{\mathrm{area}}(i,j)=\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right.$ ，这里σ是噪声整体标准差，W_i,j＝[i-(K-1)/2;i+(K-1)/2]×[j-(K-1)/2;j+(K-1)/2]为以(i,j)为中心的大小为K×K的窗口，K为奇数，C₁,C₂,C₃为3个参数。

9.依据权利要求1或8所述的基于图像形态模糊隶属度的分数阶自适应相干斑滤波方法，其特征在于：参数C₁,C₂,C₃满足0＜C₁＜C₂≤C₃＜2，其计算公式为

$\left\{\begin{array}{c}{C}_1=1-e^{-0.01/\left({2\mathrm{\&sigma;}}^2\right)},\\ C_2=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)-C_1\},\\ C_3=\min \{1.9,{2C}_1\&CenterDot;\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)/\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)+C_1\},\end{array}\right.---\left(12\right)$

其中σ和Var_loc是按照权利4计算得到的相干斑噪声标准差和局部方差矩阵，

$\left\{\begin{array}{c}\mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^{\mathrm{ET}}\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{1,1}\&le;i\&le;M,1\&le;j\&le;N\}\\ \mathrm{mean}\left({\mathrm{Var}}_{\mathrm{loc}}^C\right)=\mathrm{mean}\left\{{\mathrm{Var}}_{\mathrm{loc}}(i,j)\right|C_{\mathrm{area}}(i,j)\&NotEqual;\mathrm{0,1}\&le;i\&le;M,1\&le;j\&le;N\}\end{array}\right..$

Images