CN108416750B - 一种图像恢复方法 - Google Patents

Abstract

本发明涉及一种图像恢复方法，具体的为一种基于分数阶交叠组合稀疏全变分的图像恢复方法。结合传统模型和分数阶全变分得出分数阶交叠组稀疏全变分正则化去模糊问题的目标函数，进而得出增广拉格朗日目标函数，将求解增广拉格朗日目标函数转换为三个子问题进行求解。为了提高图像还原的运算速度，我们将图像的横向、纵向差分矩阵运算建模为卷积运算，结合周期性边界条件，将二维快速傅里叶变换巧妙应用到图像复原问题中，利用频域上的点乘操作代替空域上大型矩阵运算，从而提高运算效率。本发明在充分考虑图像所有局部和非局部特征的基础上，提高了抗噪声的鲁棒性、平滑区域与边缘区域之间的差异性和对图像边缘的保护。

Description

一种图像恢复方法

技术领域

本发明涉及图像处理领域，尤其涉及一种图像恢复方法。

背景技术

数字图像处理在生物医学、航空航天、机器视觉等众多领域有着广泛的应用。由于成像设备的限制和图像传输过程中的不利因素的影响，我们得到的图像往往出现不同程度的降质，如图像的模糊、图像噪声、图像部分信息丢失等。这对于后期的图像处理，如图像检测、图像分割、图像匹配等，产生不利的影响。如何减弱或者去除降质因素对图像质量的影响，以满足实际应用的需求，成为人们极为关注的问题。

图像恢复是利用获得的降质图像及某些先验信息，恢复和重建清晰图像，从而改善图像质量的技术。在医学成像、空间探索、遥感成像、物质检测、军事目标识别等诸多领域有着广泛的应用。合理的图像恢复模型并配合高性能的求解算法具有非常重要的意义。

Rudin等人于1992年在Physica D刊物上发表了题为“Nonlinear totalvariation based noise removal algorithms”的论文中首次开创性地提出将全变分(Total Variation,TV)变换应用到图像去噪，图像去模糊等图像复原问题中，较好地重构图像。然而，在一阶TV模型中存在两个突出问题，第一个问题，TV模型假设图像是分片光滑的，TV去噪模型在保持图像边缘方面具有非常明显的优势，但是也容易带来阶梯效应。第二个问题，传统TV模型中仅仅考虑了横纵两个方向的梯度信息，并没有充分将每个点的邻域梯度信息利用起来。

发明内容

为了解决上述问题，本发明旨在提供一种图像恢复方法，将一种新的去模糊模型应用在自然图像去模糊中，交叠组稀疏正则项的提出是为了解决阶梯效应问题，分数阶全变分为解决邻域信息挖掘不够充分的问题。然后利用交替方向乘子法(AlternatingDirection Method of Multipliers，ADMM)进行求解。同时，为提高计算效率，在交替方向乘子法中引入快速傅里叶变换(Fast Fourier Transform，FFT)，将图像域图像差分操作变换到频域，这样可以有效避免大矩阵相乘运算。

具体方案如下：

一种图像恢复方法，包括以下步骤：

S1：输入观察图像；

S2：对图像模型进行初始化设定，设定参数：停止阈值tol、二次惩罚项的惩罚系数β、正则参数μ、二维图像元素的个数M、迭代次数阈值Nit、迭代次数n＝0、循环次数k＝0，设定

S3：计算分数阶交叠组稀疏全变分正则化去模糊问题的目标函数的增广拉格朗日目标函数：

S4：将求解增广拉格朗日目标函数转换为F子问题、X₁子问题和X₂子问题三个子问题进行求解，所述三个子问题的迭代表达式为：

其中k为迭代次数；

S5：求解F子问题的计算公式：

S6：根据优化最小化算法(Majorization Minimization,MM)得出X_i(i＝1,2)子问题的迭代表达式为：

根据迭代次数阈值Nit将上述X_i(i＝1,2)子问题的迭代表达式进行Nit次迭代求解X₁和X₂；

S7：计算对偶变量和

所述对偶变量的计算公式为：

S8：判断是否满足：如果满足，进入S9，如果不满足，设定k＝k+1，返回S3；

S10：输出恢复图像。

进一步的，步骤S2中设定参数tol＝10^-4、β＝1、μ＝1、M＝20、Nit＝10。

进一步的，步骤S3中所述增广拉格朗日目标函数的求解过程为：

S31：根据传统模型和分数阶全变分得出分数阶交叠组稀疏全变分正则化去模糊问题的目标函数为：

S32：设定和X_i∈R^N×N(i＝1,2)，将目标函数修改为：

S33：利用增广拉格朗日方法，求得目标函数的增广拉格朗日目标函数为：

其中Λ_i(i＝1,2)为拉格朗日乘子，β>0为惩罚参数；

S34：利用交替方向乘子法，引入分裂变量的对偶变量与二次惩罚项β，

令设定

则增广拉格朗日目标函数修改为：

进一步的，步骤S5中所述F子问题的计算公式通过以下步骤得出：

S51：利用快速二维傅里叶变换将时域图像差分操作变换到频域，计算F子问题的频域表达式：

其中表示x的频谱，符号表示不同矩阵之间相同位置元素点乘；

S52：令频域表达式的一阶差分结果为0，利用二维傅里叶反变换得到F子问题的计算公式为：

其中表示二维傅里叶反变换，除法表示不同矩阵之间相同位置元素点除。

进一步的，步骤S6中所述X_i(i＝1,2)子问题的迭代表达式计算方法为：

根据优化最小化算法得出X_i(i＝1,2)子问题的迭代表达式为：

其中为一个对角线矩阵，对角线上的元素定义为：

其中Λ的对角元通过Matlab的内置函数conv2计算获得；

迭代循环计算公式中的表示单位矩阵，表示在第k+1次外循环中MM算法的第n次代，的计算公式为：

本发明采用如上技术方案，提出一种基于分数阶交叠组合稀疏全变分的图像恢复方法。分数阶全变分将传统整数阶TV模型延伸到分数阶模型，不仅考虑到图像的局部特性，还将图像的非局部特性考虑进去，提高了抗噪声的鲁棒性。而交叠组合稀疏全变分则将每个像素的全变分梯度推广为组合梯度，从而提高平滑区域与边缘区域之间的差异性，利用这种方法提高对图像边缘的保护。为了提高图像还原的运算速度，我们将图像的横向、纵向差分矩阵运算建模为卷积运算，结合周期性边界条件，将二维快速傅里叶变换巧妙应用到图像复原问题中，利用频域上的点乘操作代替空域上大型矩阵运算，从而提高运算效率。

附图说明

图1所示为本发明实施例一的流程示意图。

图2所示为本发明实施例一中的恢复前的图像示意图。

图3所示为本发明实施例一中的恢复后的图像示意图。

具体实施方式

为进一步说明各实施例，本发明提供有附图。这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理。配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点。图中的组件并未按比例绘制，而类似的组件符号通常用来表示类似的组件。

现结合附图和具体实施方式对本发明进一步说明。

本发明实施例提供了一种图像恢复方法，如图1所示，其为本发明实施例一所述的图像恢复方法的流程示意图，所述方法可包括以下步骤：

S1：输入观察图像。

如图2所示，图2(a)为原图像，图2(b)为模糊核，为图2(c)为被模糊的图像，图2(d)为被模糊并且加高斯噪声的图像，该实施例中所述观察图像为图2(c)的被模糊的图像，对其进行图像复原，尽量使输出的恢复图像逼近原图。

S2：对图像模型进行初始化设定，设定参数tol、β、μ、M、Nit、n＝0、k＝0，设定

其中，tol是停止阈值，该实施例中设定为tol＝10^-4。β表示二次惩罚项的惩罚系数，该实施例中设定为β＝1。μ表示正则参数，该实施例中设定为μ＝1。M表示二维图像元素的个数，该实施例中设定为M＝20。Nit表示迭代次数阈值，该实施例中设定Nit＝10，n表示迭代次数，k表示循环次数。

S3：计算交叠组稀疏全变分正则化去模糊问题的目标函数的增广拉格朗日目标函数：

所述增广拉格朗日目标函数的求解过程为：

S31：根据传统各向异性TV去模糊模型和分数阶全变分得出分数阶交叠组稀疏全变分正则化去模糊问题的目标函数为：

所述目标函数的求解过程如下：

S311：传统各向异性TV去模糊模型如下：

其中，G∈R^N×N表示被模糊并高斯噪声污染的观测图像；F∈R^N×N表示由去噪模型恢复出的图像；||·||₂表示欧式L₂范数；R_ATV(F)表示基于各向异性TV的稀疏约束项；式(3)第一项称为保真项，第二项μR_ATV(F)对真实图像的先验假设称为正则项，μ为正则参数，用来衡量保真项与正则项之间的权重；H∈R^N×N表示模糊核，*为卷积算符。

传统各向异性TV(Anisotropic TV,ATV)定义如下：

R_ATV(F)＝||K_h*F||₁+||K_v*F||₁ (4)

其中K_h＝[-1,1]表示横向卷积核，K_v＝[-1,1]^T表示纵向卷积核。

S312：交叠组稀疏全变分正则化去模糊问题的目标函数为：

其中和分别表示横向和竖向组合梯度，考虑交叠组合稀疏的正则项表示为在图像恢复过程中，若考虑到图像邻域M×M个点的信息(M为组合的大小)，定义：

二维图像的M×M个元素构成的组块定义为：

其中， 向下取整，表示小于或等于x的最大整数；元素V_i,j位于矩阵的中心位置，选取的组将中心围起来。当M＝1时，交叠组合稀疏模型就退化为ATV模型。采用这种组合稀疏的方式作为正则项，充分考虑到图像邻域信息，从而有效抑制了“阶梯效应”。

S313：分数阶全变分是全变分的另一种延伸，它将传统的差分算子扩展到分数阶差分算子，更充分挖掘邻域梯度信息。τ阶全变分去模糊去噪模型可以描述为：

其中和分别代表了水平方向和垂直方向的τ阶梯度算子。有限的分数阶差分定义为：

其中广义二项式系数为：

伽马函数为：

为简单起见，可以从以下式子得到递归：

其中，当：τ＝1，时分数阶模型退化为ATV模型。

S314：根据传统模型和分数阶全变分得出交叠组稀疏全变分正则化去模糊问题的目标函数为：

S32：设定和X_i∈R^N×N(i＝1,2)，将式(13)的交叠组稀疏全变分正则化去模糊问题的目标函数修改为：

S33：利用增广拉格朗日方法(Augmented Lagrangian Multiplier method,ALM)解决形如式(14)的约束优化问题，将其转换为无约束问题，其增广拉格朗日目标函数表示为：

其中Λ_i(i＝1,2)为拉格朗日乘子，β>0为惩罚参数。

S34：对变量X₁,X₂联合求解最值问题相对困难，为了解决上述难点，我们采用了交替方向乘子法。ADMM的思想是将原问题中X₁,X₂分开独立求解，从而收敛到原模型的最优解。为此，引入分裂变量的对偶变量(即拉格朗日系数)与二次惩罚项，令并对目标函数补配方项则有：

S4：由于变量之间相互独立，将求解增广拉格朗日目标函数转化为三个子问题进行单独求解，于是得到如下的迭代格式，k为迭代次数：

下面，逐一给出求解每一个子问题的方法:

S5：求解F子问题的计算公式：

对于F子问题，直接从空域求解，要将目标函数中的卷积算子写成Toeplitz矩阵形式，并且将G和F进行列化操作，导致涉及的矩阵规模为计算复杂度极高。为了有效避免大矩阵相乘运算带来的计算复杂度，我们引入快速二维傅里叶变换将时域图像差分操作变换到频域。F子问题的频域表达式为：

式(19)中表示x的频谱，符号表示不同矩阵之间相同位置元素点乘。

由于式(19)等式右边的函数是光滑的凸函数，可令其一阶差分结果为0，从而求解最优值。利用二维傅里叶反变换得到F子问题的最优解，整理得：

式(20)中，表示二维傅里叶反变换；除法表示不同矩阵之间相同位置元素点除。

S6：根据优化最小化算法(Majorization Minimization,MM)得出X₁和X₂子问题的迭代表达式为：

根据迭代次数阈值Nit将上述X₁和X₂子问题的迭代表达式进行Nit次迭代求解X₁和X₂；

对于X_i(i＝1,2)子问题，为了求解包含交叠组合稀疏的最优解，我们采用了优化最小化(Majorization Minimization,MM)算法。式(17)中的X₁,X₂子问题可由下列迭代循环求解，

是一个对角线矩阵，对角线上的元素定义如下，

Λ的对角元可由Matlab的内置函数conv2计算获得。式(21)中，表示单位矩阵；表示在第k+1次外循环中MM算法的第n次迭代，

S7：计算对偶变量和

根据交替方向乘子法，对偶变量的更新方式如下：

S8：判断是否满足：如果满足，进入S9，如果不满足，设定k＝k+1，返回S3；

S10：输出恢复图像。

如图3所示，图3(a)为使用传统各向异性ATV方法的恢复图像，图3(b)使用分数阶TV方法的恢复图像，为图3(c)为使用交叠组稀疏TV方法的恢复图像，图3(d)为使用该实施例中方法的恢复图像，对比各种正则项计算得到的恢复图像，可以看到，ATV方法阶梯效应最为严重，分数阶TV方法中的阶梯效应依然存在，交叠组合稀疏方法的阶梯效应得到一定的缓解，而我们提出的方法恢复的图像则相对于其他方法的恢复图像更好，阶梯效应去除得更为明显。

本发明实施例提出一种基于分数阶交叠组合稀疏全变分的图像恢复方法。分数阶全变分将传统整数阶TV模型延伸到分数阶模型，不仅考虑到图像的局部特性，还将图像的非局部特性考虑进去，提高了抗噪声的鲁棒性。而交叠组合稀疏全变分则将每个像素的全变分梯度推广为组合梯度，从而提高平滑区域与边缘区域之间的差异性，利用这种方法提高对图像边缘的保护。为了提高图像还原的运算速度，我们将图像的横向、纵向差分矩阵运算建模为卷积运算，结合周期性边界条件，将二维快速傅里叶变换巧妙应用到图像复原问题中，利用频域上的点乘操作代替空域上大型矩阵运算，从而提高运算效率。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。

Claims (4)

1.一种图像恢复方法，其特征在于：包括以下步骤：

S1：输入观察图像；

S2：初始化设定：设定参数：停止阈值tol、二次惩罚项的惩罚系数β、正则参数μ、二维图像元素的个数M、迭代次数阈值Nit、迭代次数n＝0、循环次数k＝0，设定

S3：计算分数阶交叠组稀疏全变分正则化去模糊问题的目标函数的增广拉格朗日目标函数：

其中G∈R^N×N表示观测图像，F∈R^N×N表示恢复出的图像，H∈R^N×N表示模糊核，*为卷积算符，||·||₂表示欧式L₂范数，K_h＝[-1,1]表示横向卷积核，K_v＝[-1,1]^T表示纵向卷积核，和分别代表了水平方向和垂直方向的τ阶梯度算子，为X₁,X₂的对偶变量；

其中，系数可以从下式递归得到：

函数定义为：

其中为二维图像的M×M个元素构成的组块，即：

其中， 向下取整，表示小于或等于x的最大整数，元素V_i,j位于矩阵的中心位置；

S4：将求解增广拉格朗日目标函数转换为F子问题、X₁子问题和X₂子问题三个子问题进行求解，所述三个子问题的迭代表达式为：

其中k为迭代次数；

S5：求解F子问题的计算公式：

其中，表示x的频谱，符号ο表示不同矩阵之间相同位置元素点乘，表示二维傅里叶反变换；除法表示不同矩阵之间相同位置元素点除；

S6：根据优化最小化算法得出X_i(i＝1,2)子问题的迭代表达式为：

其中，是一个对角线矩阵，对角线上的元素定义如下，

表示单位矩阵；表示在第k+1次外循环中优化最小化算法的第n次迭代，

根据迭代次数阈值Nit将上述X_i(i＝1,2)子问题的迭代表达式进行Nit次迭代求解X₁和X₂；

S7：计算对偶变量和

所述对偶变量的计算公式为：

S8：判断是否满足：如果满足，进入S9，如果不满足，设定k＝k+1，返回S3；

S10：输出恢复图像。

2.根据权利要求1所述的图像恢复方法，其特征在于：步骤S2中设定参数tol＝10^-4、β＝1、μ＝1、M＝20、Nit＝10。

3.根据权利要求1所述的图像恢复方法，其特征在于：步骤S3中所述增广拉格朗日目标函数的求解过程为：

S31：根据传统模型和分数阶全变分得出分数阶交叠组稀疏全变分正则化去模糊问题的目标函数为：

S32：设定和X_i∈R^N×N(i＝1,2)，将目标函数修改为：

S33：利用增广拉格朗日方法，求得目标函数的增广拉格朗日目标函数为：

其中Λ_i(i＝1,2)为拉格朗日乘子，β＞0为惩罚参数；

S34：利用交替方向乘子法，引入分裂变量的对偶变量与二次惩罚项β，

令设定

则增广拉格朗日目标函数修改为：

4.根据权利要求1所述的图像恢复方法，其特征在于：步骤S5中所述F子问题的计算公式通过以下步骤得出：

S51：利用快速二维傅里叶变换将时域图像差分操作变换到频域，计算F子问题的频域表达式：

其中表示x的频谱，符号ο表示不同矩阵之间相同位置元素点乘；

S52：令频域表达式的一阶差分结果为0，利用二维傅里叶反变换得到F子问题的计算公式为：

其中表示二维傅里叶反变换，除法表示不同矩阵之间相同位置元素点除。