CN102810202B - 基于分数阶差分加权的图像多步残差反馈迭代滤波方法 - Google Patents

Abstract

本发明公开了一种基于分数阶差分加权的图像多步残差反馈迭代滤波方法，该滤波方法包括分数阶奇异性指标计算单元、分数阶权重矩阵计算单元以及多步残差反馈滤波单元。本方法首先估计每个像素点对应的分数阶奇异性指标；然后依据分数阶差分格式的系数计算方式生成分数阶加权系数矩阵；最后，执行多步残差反馈滤波迭代，利用多个去噪残差图像的分数阶加权组合更新生成中间待去噪图像，并对其进行全变差滤波生成迭代去噪图像序列，迭代收敛图像为最终去噪图像。本方法得到的迭代序列可以快速收敛于具有较高峰值信噪比的去噪图像，迭代结果对于迭代终止条件敏感性及依赖程度小，可以在有效抑制图像噪声的同时，较好保持图像的纹理等细节信息。

Description

基于分数阶差分加权的图像多步残差反馈迭代滤波方法

技术领域

本发明涉及图像处理领域中图像噪声抑制的滤波技术，特别是一种基于分数阶差分加权的图像多步残差反馈迭代滤波方法。

背景技术

图像在采集、传输和显示过程中，不可避免地会产生噪声污染，噪声抑制是图像分割、识别与目标检测等许多图像后续过程的预处理过程。而在一些应用领域，例如在医学和遥感图像处理中，图像强、弱边缘的检测以及纹理分析都是重要的图像分析手段，因此在抑制噪声的同时有效保持图像的边缘、纹理等结构非常重要。

目前,国际上关于抑制噪声过程中的图像细节保持问题已有许多研究。2004年,Chamolle提出的全变差滤波的快速投影迭代算法(ChambolleA.Analgorithmfortotalvariationminimizationandapplications.JournalofMathematicalImagingandVision,2004，20(1)：89-97)计算量小，收敛速度快，能够较好地保持图像的边缘,但是容易出现边缘阶梯效应,纹理等细节信息丢失严重。虽然非局部TV正则化方法(Bresson,X.Ashortnotefornonlocal-TVminimization.June2009,technicalreport)可以较好地保持图像纹理细节，但这种方法计算量较大。Osher等人提出的迭代正则化方法(S.Osher,M.Burgeretal.AnIterativeRegularizationMethodforTotalVariationBasedImageRestoration.MultiscaleModelingandSimulation,2005,4(2):460-489)可以较好地保持图像的细节，并且计算量较小，因此受到了很大的关注。该方法实际上是一种单步残差反馈迭代滤波方法，在迭代过程中，残差图像(即在滤波过程中被去掉的部分)中的图像纹理等细节信息包括噪声都被逐步加回到去噪图像中去，最终所得到的迭代序列收敛到原含噪声图像，因此可以较好地保持图像的纹理细节信息。然而，为了防止过多的噪声被加回到去噪图像中，该迭代方法必须算法收敛之前在适当的时候终止。然而，由于该方法收敛速度非常快，因此迭代结果对于迭代终止步的选择非常敏感，这大大影响了该方法在实际应用中的稳定性。

发明内容

本发明目的是提供一种基于分数阶差分加权的图像多步残差反馈迭代滤波方法。该滤波方法基于分数阶差分格式，对迭代过程中所产生的多幅残差图像进行加权组合，并反馈到原始待去噪图像中生成新的待去噪图像进行去噪，通过选择合适的参数，所得到的迭代序列的峰值信噪比可以很稳定地保持在一个较高的水平，可以通过经典的迭代序列相对误差就可以很方便地控制迭代终止条件，有利于实际应用。该方法算法计算量小，收敛速度快，可以在去除图像噪声的同时，有效保持图像的纹理等细节成分，得到具有高峰值信噪比和良好视觉效果的去噪图像。

实现本发明目的的技术解决方案为：本发明基于分数阶差分加权的图像多步残差反馈迭代滤波方法包括分数阶奇异性指标计算单元、分数阶权重矩阵计算单元以及多步残差反馈滤波单元。本方法首先在分数阶奇异性指标计算单元中估计每个像素点对应的分数阶奇异性指标；然后在分数阶权重计算单元中，依据分数阶差分格式的系数计算方式生成分数阶加权系数矩阵；最后，执行多步残差反馈滤波单元，利用多个去噪残差图像的分数阶加权组合更新生成中间待去噪图像，并对其进行全变差滤波生成迭代去噪图像序列，迭代收敛图像为最终去噪图像。

1.1所述的分数阶奇异性指标计算单元,包括如下步骤:

①全变差滤波预处理。输入一幅大小为M×N的待去噪图像f，应用全变差滤波方法得到初始滤波图像u₁；

②计算残差图像的局部方差。计算初始残差图像v₁＝f-u₁,逐个计算残差图像中像素(i,j)所在局部大小为K×K(5≤K≤13,K为奇数)窗口内的局部方差LV(i,j)，其计算关系为:

$LV(i,j)=\frac{1}{K^2}\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1^2(p,q)-\frac{1}{K^4}{\left(\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1\right(p,q\left)\right)}^2,i=1,2,\mathrm{...},M;j=1,2,\mathrm{...},N$

其中W_i,j＝[i-(K-1)/2；i+(K-1)/2]×[j-(K-1)/2；j+(K-1)/2]是残差图像中以(i,j)为中心,大小为K×K的窗口。

③确定分数阶奇异性指标分布区间。计算局部方差LV(i,j)(i＝1,2,…,M；j＝1,2,…,N)的平均值将分数阶奇异性指标分布区间记为[c₁,c₂]，其中c₁,c₂的取值规则为:

(i)当取c₁＝0.4,c₂＝0.6；

(ii)当取c₁＝0.6,c₂＝0.8；

(iii)当取c₁＝0.8,c₂＝1.0；

其中:为噪声标准差，由下面的经验估计公式进行估计

$\widetilde{\mathrm{\σ}}=\frac{Midium\left(v_{1,HH}\right)}{0.6745}$

式中Midium(v_1,HH)为残差图像v₁在Db4小波小波分解的最高频的HH子带小波系数的幅度中值。

④计算分数阶奇异性指标。利用Sigmoid函数，将局部方差LV(i,j)的值限定到在分数阶奇异性指标分布区间[c₁,c₂]上，并作为对图像内各点分数阶奇异性指标α(i,j)的估计，其具体计算关系为

$\mathrm{\α}(i,j)=c_1+\frac{c_2-c_1}{1+e^{-\left(LV\right(i,j)-\stackrel{\‾}{LV})/\mathrm{\σ}\left(LV\right)}},$

其中σ(LV)分别为局部方差LV(i,j)(i＝1,2,…,M；j＝1,2,…,N)的标准差。

1.2所述的分数阶权重矩阵计算单元

本发明中的残差反馈是利用离散分数阶微分来实现的,在本发明中,图像中每一点(i,j)处的分数阶差分阶数,即1.1中所述的分数阶奇异性指标α(i,j)可以不同，图像中的每一点对应的分数阶微分组合权重系数也不一样，因此本发明图像中各点的分数阶权重系数构成了若干个空间变化的分数阶权重矩阵，本计算单元主要用于计算和存储这些分数阶权重矩阵。

设立L(2≤L≤5)个大小为M×N的矩阵存储装置,存储L个分数阶权重矩阵w₁,w₂,…,w_L。权重矩阵w_l中每个元素w_l(i,j)的计算公式为：

$w_l(i,j)={(-1)}^l{C}_l^{\mathrm{\α}(i,j)},$

其中Γ(x)为Gamma函数，且当a(i,j)≤l-1时，

1.3所述的多步残差反馈迭代滤波单元,包括如下步骤：

步骤1：设立L(2≤L≤5)个大小为M×N的矩阵存储装置,存储当前第k+1步前面的L个残差图像其中k为迭代次数；k＝1时 $v_1^{\left(1\right)}=v_1,v_2^{\left(1\right)}=0,v_2^{\left(1\right)}=\mathrm{0...},v_L^{\left(1\right)}=0;$

步骤2:执行迭代,对L个残差图像进行线性组合，并从原含噪声图像f中减去该残差线性组合，得到第k步(k≥1)中间待去噪图像g_k,其计算关系为:

$g_k=f-{\mathrm{\Σ}}_{j=1}^L{w}_j\·v_j^{\left(k\right)}$

其中每一个的运算都是矩阵中对应元素的点对点的乘法。

步骤3：对第k步(k≥1)中间图像g_k，利用全变差滤波方法进行快速去噪，得到新的中间去噪图像u_k+1及其残差图像ERR＝g_k-u_k+1；

步骤4：判断迭代终止条件||u_k+1-u_k||₂/||u_k|||₂≤ε₂，如果满足迭代终止条件，则输出最终的去噪结果u_k+1；否则，令k:＝k+1，更新残差图像存储装置中的数据,即转到迭代计算过程中的步骤2继续计算；

本发明与现有技术相比，其显著优点：本发明方法能够实现高视觉质量的快速去噪。本发明方法得到的迭代序列可以快速收敛并稳定地保持在较高峰值信噪比水平上，从而大大降低了迭代结果对于迭代终止条件敏感性及依赖程度。本方明方法计算量小，收敛速度快，可以在有效抑制图像噪声的同时，较好保持图像的纹理等细节信息，得到具有良好视觉效果的去噪图像。本发明在遥感图像处理、医学影像处理等都有广泛的应用前景。

附图说明

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

图2是分数阶奇异性指标计算单元算法及数据流程图。

图3是多步残差反馈迭代滤波单元算法及数据流程图。

图4是本发明的实验测试图像，(a)纯卡通含噪声图像，(b)纯纹理含噪声图像，(c)含噪声的自然图像。

图5是本发明方法在不同分数阶阶数下的去噪图像峰值信噪比比较，(a)纯卡通图像去噪峰值信噪比，(b)纯纹理图像去噪峰值信噪比。

图6是本发明方法不同分数阶阶数下与与传统方法的纯卡通图像的去噪比较，(a)Bregman迭代，(b)α＝0.8，(c)α＝0.6，(d)α＝0.4，(e)ROF。

图7是本发明方法不同分数阶阶数下与与传统方法的纯纹理图像的去噪比较，(a)Bregman迭代，(b)α＝0.8，(c)α＝0.6，(d)α＝0.4，(e)ROF。

图8是本发明方法与传统方法的自然图像去噪比较，(a)传统迭代正则化方法(这里的去噪图像是迭代中具有最高峰值信噪比的去噪图像u₃)，(b)非局部TV正则化方法，(c)本文自适应方法。

具体实施方式

下面结合附图对本发明做进一步说明。

首先，结合附图1至图3给的本发明方法的整体算法结构以及具体单元算法流程,详细介绍本发明的使用方法。图1中给出了本发明的基本流程及功能模块，依次为分数阶奇异性指标计算单元、分数阶权重矩阵计算单元以及多步残差反馈滤波单元。

1.1所述的分数阶奇异性指标计算单元,如图2所示,按照下面的具体方法进行计算:

①全变差滤波预处理。输入一幅大小为M×N的待去噪图像f，应用全变差滤波方法得到初始滤波图像u₁。在这一步，所采用的，在该快速投影迭代算法中，设定前后相邻两次迭代相对误差小于给定的精度(一般设定相对误差精度为10^-3)时，迭代终止，并将此时的迭代结果作为全变差滤波得到的去噪图像。在本发明中所有述的全变差滤波方法都采用Chamolle的快速投影迭代算法。

②计算残差图像的局部方差。残差图像每个像素(i,j)局部方差LV(i,j)的计算关系为:

$LV(i,j)=\frac{1}{K^2}\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1^2(p,q)-\frac{1}{K^4}{\left(\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1\right(p,q\left)\right)}^2,i=1,2,\mathrm{...},M;j=1,2,\mathrm{...},N$

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

③确定分数阶奇异性指标分布区间。其具体实现过程是：

步骤1，利用Db4小波对残差图像v₁进行一次小波分解，得到最高频的HH子带小波系数，对该子带小波系数幅度进行排序后取中值记为Midium(v_1,HH)，然后估计图像的噪声标准差其估计公式为

$\widetilde{\mathrm{\σ}}=\frac{Midium\left(V_{1,HH}\right)}{0.6745}$

步骤2，计算局部方差的平均值

步骤3，确定分数阶奇异性指标分布区间[c₁,c₂]，其中c₁,c₂的取值规则为:

(i)当取c₁＝0.4,c₂＝0.6；

(ii)当取c₁＝0.6,c₂＝0.8；

(iii)当取c₁＝0.8,c₂＝1.0；

④计算分数阶奇异性指标。利用所得到的残差图像的局部方差及分数阶奇异性指标分布区间[c₁,c₂]，计算图像内各点分数阶奇异性指标α(i,j)，其具体计算关系为

$\mathrm{\α}(i,j)=c_1+\frac{c_2-c_1}{1+e^{-5\×\left(LV\right(i,j)-\stackrel{\‾}{LV})/\mathrm{\σ}\left(LV\right)}},$

其中σ(LV)分别为局部方差LV(i,j)(i＝1,2,…,M；j＝1,2,…,N)的标准差。

1.2所述的分数阶权重矩阵计算单元，按照下面的具体方法进行计算：

设立L(2≤L≤5)个大小为M×N的矩阵存储装置,存储L个分数阶权重矩阵w₁,w₂,…,w_L。权重矩阵w_l中每个元素w_l(i,j)的计算公式为：

$w_l(i,j)={(-1)}^l{C}_l^{\mathrm{\α}(l,j)},$

其中Γ(x)为Gamma函数，且当a(i,j)≤l-1时，

1.3所述的多步残差反馈迭代滤波单元,其算法即数据流程如图3所示，具体计算按照下面的方法进行：

步骤1：设立L(2≤L≤5)个大小为M×N的矩阵存储装置,存储当前第k+1步前面的L个残差图像其中k为迭代次数；k＝1时 $v_1^{\left(1\right)}=v_1,v_2^{\left(1\right)}=0,v_2^{\left(1\right)}=\mathrm{0...},v_L^{\left(1\right)}=0;$

步骤2:计算中间待去噪图像得到中间待去噪图像g_k,其计算关系为:

$g_k=f-{\mathrm{\Σ}}_{j=1}^L{w}_j\·v_j^{\left(k\right)}$

其具体计算过程是：首先，对存贮装置中的对应数据进行点乘运算，具体为这里的每一个的运算都是矩阵中对应元素的点对点的乘法；然后从原含噪声图像矩阵f中减去将点乘运算之后得到的L个矩阵，得到中间待去噪图像g_k。

步骤3：对中间图像g_k，利用Chambolle快速投影迭代算法进行全变差滤波，得到新的中间去噪图像u_k+1及其残差图像ERR＝g_k-u_k+1；

步骤4：判断迭代终止条件||u_k+1-u_k||₂/||u_k|||₂≤ε(本发明中预设ε＝10^-4)，如果满足迭代终止条件，则输出最终的去噪结果u_k+1；否则，令k:＝k+1，更新残差图像存储装置中的数据,即转到迭代计算过程中的步骤2继续计算；

在本次实例中，我们考虑3种图像(图4)：纯卡通图像(256×256)、纯纹理图像(256×256)以及既有卡通部分也有纹理部分的一般自然图像(256×256)。本发明方法中的正则化参数λ＝0.01，本发明方法的整体迭代的迭代误差限ε＝10^-4。作为衡量去噪效果，本实验中采用峰值信噪比(PSNR)作为客观标准，这里PSNR定义为：

$PSNR\left(u\right)=10\lg \frac{\underset{1\≤i,j\≤N}{\max }|{\left(u_0\right)}_{i,j}{|}^2}{\frac{1}{N^2}\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\[{\left(u_0\right)}_{i,j}-u_{i,j}\]}^2}$

图5给出了在不同固定分数阶奇异性指标下(即对图像中所有的点都取相同的分数阶奇异性指标)，本发明方法对于纯卡通图像和纯纹理图像的峰值信噪比的改善。从峰值信噪比的变化图可以看到，传统的单步残差反馈迭代正则化方法的去噪图像的峰值信噪比开始上升非常快，在达到最大值之后，迅速衰减到很低的程度，而且这个过程只需要很少的迭代步骤，这就说明想要得到具有较高峰值信噪比的图像，迭代终止条件选择非常重要。而本发明方法可以明显改善这一状况。对于纯卡通图像，当α位于0.5附近时，峰值信噪比上升之后，随着迭代进行经过非常缓慢的下降之后，很快趋于稳定保持在较高的水平，在后面的迭代过程中变化非常小。而对于纯纹理图像，当α比1.0略小的时候，峰值信噪比上升时候，也很快稳定在较高的水平。这说明选择合适的参数时，我们可以很方便地得到具有较高峰值信噪比的去噪图像。图6-图7给出了不同参数下的去噪图像的比较。对于纯卡通图像，当α位于1.0附近时，过多的噪声被加回到去噪图像中；当α更加靠近0.0时，噪声被很好地抑制，但从残差图像中可以看到图像边缘损失也是比较严重的；而当α位于0.5附近时，既可以较好抑制噪声，也可以较好保持图像的边缘结构。对于纯纹理图像，相对于α＝1.0，即传统的迭代正则化方法，在当α位于1.0附近且比1.0时，已经可以去除大部分噪声，虽然在图像中仍然保留了一些噪声，但由于纹理对于噪声有视觉上的掩盖效应，因此这部分噪声对于视觉效果影响并不是太大；但当α较小时，会损失过多的纹理细节。对于试验一般自然图像(图4(c)),我们根据1.1所述方法计算分数阶奇异性指标，这样各象素点对应的分数阶奇异性指标是自适应的。实验结果表明，本发明能在抑制图像非纹理区域中的噪声的同时，可以较好地保持图像纹理区域中的纹理信息。

Claims (1)

1.一种基于分数阶差分加权的图像多步残差反馈迭代滤波方法，包括分数阶奇异性指标计算、分数阶权重系数矩阵计算以及多步残差反馈迭代滤波；

1.1所述的分数阶奇异性指标计算，包括如下步骤：

步骤1.1.1：全变差滤波预处理：输入一幅大小为M×N的待去噪图像f，应用全变差滤波方法得到初始滤波图像u₁；

步骤1.1.2：计算初始滤波残差图像v₁＝f-u₁的每个像素点的局部方差LV(i,j)；残差图像的每个像素点的局部方差LV(i,j)满足如下关系：

$LV(i,j)=\frac{1}{K^2}\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1^2(p,q)-\frac{1}{K^4}{\left(\underset{(p,q)\∈W_{i,j}}{\Σ}{v}_1(p,q)\right)}^2,i=1,2,\mathrm{...},M;j=1,2,\mathrm{...},N$

其中：v₁＝f-u₁，W_i,j＝[i-(K-1)/2；i+(K-1)/2]×[j-(K-1)/2；j+(K-1)/2]是一个以(i,j)为中心,大小为K×K的窗口，5≤K≤13,K为奇数；

步骤1.1.3：确定分数阶奇异性指标分布区间[c₁,c₂]；按照如下方法进行确定：计算局部方差LV(i,j)的平均值i＝1,2,…,M；j＝1,2,…,N，c₁,c₂的取值规则为:

(i)当取c₁＝0.4,c₂＝0.6；

(ii)当取c₁＝0.6,c₂＝0.8；

(iii)当取c₁＝0.8,c₂＝1.0；

其中:为噪声标准差，由下面的经验估计公式估计

$\widetilde{\mathrm{\σ}}=\frac{Midium\left(v_{1,HH}\right)}{0.6745}$

式中Midium(v_1,HH)为残差图像v₁在Db4小波小波分解的最高频的HH子带小波系数的幅度中间值；

步骤1.1.4：计算图像每个像素点对应的分数阶奇异性指标α(i,j)；分数阶奇异性指标α(i,j)的具体计算关系为

$\mathrm{\α}(i,j)=c_1+\frac{c_2-c_1}{1+e^{-5\×\left(LV\right(i,j)-\stackrel{\‾}{LV})/\mathrm{\σ}\left(LV\right)}},$

其中σ(LV)为局部方差LV(i,j)的标准差，i＝1,2,…,M；j＝1,2,…,N；

1.2所述的分数阶权重系数矩阵计算，通过设立L个大小为M×N的矩阵存储装置,2≤L≤5，存储计算得到的L个分数阶权重系数矩阵w₁,w₂,…,w_L，分数阶权重系数矩阵计算中权重系数矩阵w_l中每个元素w_l(i,j)都是按照分数阶差分格式中的系数计算方式计算,1≤l≤L，其计算公式为：

$w_l(i,j)={(-1)}^l{C}_l^{\mathrm{\α}(i,j)},$

其中Γ(x)为Gamma函数，且当a(i,j)≤l-1时，

1.3所述的多步残差反馈迭代滤波,包括如下步骤：

步骤1.3.1：设立L个大小为M×N的矩阵存储装置,2≤L≤5，在第k+1步计算开始时，存储前一步即第k步到第k+1-L步迭代中产生的L个相邻的残差图像其中k为迭代次数；k＝1时 $v_1^{\left(1\right)}=v_1,v_2^{\left(1\right)}=0,v_2^{\left(1\right)}=0...,v_L^{\left(1\right)}=0;$

步骤1.3.2:按照分数阶差分格式，对得到的L个残差图像利用分数阶权重系数矩阵w₁,w₂,…,w_L进行线性组合，并从原含噪声图像f中减去该残差线性组合，得到中间待去噪图像g_k,其计算关系为:

$g_k=f-{\Σ}_{j^{\′}=1}^L{w}_{j^{\′}}\·v_{j^{\′}}^{\left(k\right)}$

其中每一个的运算都是矩阵中对应元素的点对点的乘法；

步骤1.3.3：对中间待去噪图像g_k，利用全变差滤波方法进行快速去噪，得到新的中间滤波图像u_k+1及其残差图像ERR＝g_k-u_k+1；

步骤1.3.4：判断迭代终止条件||u_k+1-u_k||₂/||u_k|||₂≤ε，其中ε为人工设定的相对误差精度要求，如果满足迭代终止条件，则输出最终的去噪图像u_k+1；否则，令k:＝k+1，更新残差图像存储装置中的数据,即转到步骤1.3.2继续计算。