在现代化的图像应用中，例如航天图像，生物医学成像，计算摄影和显微成像等，单幅图像的盲反卷积已经引起了广泛的关注，对于整图唯一的模糊核$\mathit{\boldsymbol{k}}$，图像运动模糊可表示为

$\mathit{\boldsymbol{B}} = \mathit{\boldsymbol{L}} \otimes \mathit{\boldsymbol{k + \xi }}$

(1)

式中，⊗表示卷积操作，$\mathit{\boldsymbol{\xi }}$表示加性高斯白噪声，$\mathit{\boldsymbol{L}}$和$\mathit{\boldsymbol{B}}$分别表示清晰图像和观测得到的模糊图像。单幅图像盲反卷积就是由$\mathit{\boldsymbol{B}}$中同时估计出$\mathit{\boldsymbol{L}}$和$\mathit{\boldsymbol{k}}$，这是一个极具挑战性的不适定(ill-posed)反问题。

为了减轻盲反卷积的病态性，往往需要将图像先验与模糊核先验作为正则项加入图像盲去模糊模型中。Fergus^[1]指出自然图像梯度服从重尾分布，并且假设模糊核服从混合指数分布，在其盲去模糊模型中以混合高斯模型刻画重尾分布，然后变分推演估计出模糊核。Shan^[2]以分段函数刻画图像的重尾分布先验，同时为了抑制振铃效应引入局部先验，提高了模糊核的估计精度。Krishnan^[3]提出一种非自然图像先验，在其盲去模糊模型中以l₁/l₂范数约束图像梯度稀疏性，以l₁范数约束模糊核稀疏性取得了较好的效果。之后Kotera^[4]和Xu^[5]分别提出基于l_0.3范数和基于逼近l₀范数的MAP型盲去模糊算法。Shao^[6]引入l₀-l₂正则项同时约束待估图像与模糊核，并利用交替迭代、算子分裂和增广拉格朗日法求解其盲去模糊模型，导出一种快速图像盲去模糊算法。Daniele和Paolo^[7]重新审视了基于l₁范数的盲去模糊方法并且提升了模糊核的估计精度，详细地分析了所提方法成功的相关细节。Zuo^[8]采用l_p范数约束图像梯度并采用数据驱动方式为盲去模糊的各个阶段学习出相应先验。此外，在某种特定条件下，模糊核可以被假设为分段且足够光滑的曲线，因此有学者^{[5, 9-10]}认为模糊核服从高斯分布并在其盲去模糊模型中引入l₂范数约束模糊核解空间。

在运动图像盲复原中，图像强边缘有助于提高模糊核估计精度^{[6, 10-11]}。为此，本文先结合图像梯度稀疏性，利用其l₀范数作为正则项估计出中间图像的强边缘。运动图像的模糊核具有稀疏性，以l₀范数约束模糊核，同时，为确保模糊核中元素不发生突变，即模糊核的连续性，采用l₂范数约束模糊核的梯度，此外，模糊核元素往往具有归一化特性，在此基础上以强边缘指导模糊核估计。结合以上先验，本文提出了基于强边缘的盲去模糊方法。为进一步减轻非盲去模糊过程中的振铃效应且更多地保留图像细节，分别采用TV范数和l_p范数作为正则项求解清晰图像，然后再将求得的结果取平均，作为最终的清晰图像。实验结果表明，提出的盲去模糊算法能够准确地估计模糊核和复原图像。

目前，众多盲去模糊方法采用先估计中间清晰图像，然后利用中间清晰图像的梯度信息估计模糊核。研究表明^{[5, 12]}直接利用图像像素信息估计获得的模糊核精度不高，使用图像梯度信息则可提高模糊核估计精度。为此，采用先估计清晰图像的梯度，然后由梯度图像估计模糊核以提高其估计精度。考虑到图像梯度的稀疏性先验，引入l₀范数约束图像梯度；考虑到模糊核的稀疏性，引入l₀范数约束模糊核；为确保模糊核元素不发生突变，即模糊核具有连续性，引入l₂范数约束模糊核梯度；结合以上先验以及模糊核元素的归一化特性，提出了基于l₀-l₂范数的多先验协同盲去模糊模型，即

$\begin{array}{l} \mathop {\min }\limits_{x,k} \frac{1}{2}\left\| {\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \gamma {\left\| \mathit{\boldsymbol{x}} \right\|_0} + \\ \quad \quad \quad {\eta _1}{\left\| \mathit{\boldsymbol{k}} \right\|_0} + {\eta _2}\left\| {\nabla \mathit{\boldsymbol{k}}} \right\|_2^2\\ \quad \quad \quad \quad {\rm{s}}{\rm{.t}}.\;\sum\limits_i {{k_i} = 1} \end{array}$

(2)

式中，$\mathit{\boldsymbol{x}}$表示$\nabla \mathit{\boldsymbol{L = }}{\left( {{\nabla _h}\mathit{\boldsymbol{L}},{\nabla _v}\mathit{\boldsymbol{L}}} \right)^{\rm{T}}}$，$\mathit{\boldsymbol{y}}$表示$\nabla \mathit{\boldsymbol{B = }}{\left( {{\nabla _h}\mathit{\boldsymbol{B}},{\nabla _v}\mathit{\boldsymbol{B}}} \right)^{\rm{T}}}$，${{\nabla _h}}$，${{\nabla _v}}$分别表示水平方向的梯度算子和垂直方向上的梯度算子；${k_i}$表示模糊核$\mathit{\boldsymbol{k}}$中第$i$个元素；$\gamma $，${\eta _1}$，${\eta _2}$表示正则化系数。式(2) 包含4项：第1项是数据保真项，其作用是使清晰图像和模糊核之间的卷积尽可能逼近观测得到的模糊图像；第2项是约束梯度图像的l₀正则项，其作用是使估计出来的梯度图像保留强边缘，丢弃一些弱分量，这有利于模糊核的估计；第3和第4项是模糊核的l₀、l₂范数约束项，这两项确保估计得到的模糊核具有一定的稀疏性的同时保持了模糊核的连续性。

Xu^[5]和Pan^[10]指出并非所有的边缘都有利于模糊核的估计，因此采用Pan^[10]提出的自适应加权方式对待估清晰图像的不同区域采用不同的正则化参数，以抑制微小细节，从而恢复出有利于模糊核估计的强边缘。于是盲去模糊模型式(2) 改写为

$\begin{array}{l} \mathop {\min }\limits_{x,k} \frac{1}{2}\left\| {\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \gamma \mathit{\boldsymbol{\kappa }} \circ {\left\| \mathit{\boldsymbol{x}} \right\|_0} + \\ \quad \quad \quad {\eta _1}{\left\| \mathit{\boldsymbol{k}} \right\|_0} + {\eta _2}\left\| {\nabla \mathit{\boldsymbol{k}}} \right\|_2^2\\ \quad \quad \quad \quad {\rm{s}}{\rm{.t}}.\;\sum\limits_i {{k_i} = 1} \end{array}$

(3)

式中，$\mathit{\boldsymbol{\kappa }} \circ {\left\| \mathit{\boldsymbol{x}} \right\|_0} = \sum\limits_p {\mathit{\boldsymbol{\kappa }}\left( p \right){{\left\| {\mathit{\boldsymbol{x}}\left( p \right)} \right\|}_0}} $，$\mathit{\boldsymbol{\kappa }}$为自适应权值，$\mathit{\boldsymbol{\kappa }}$=exp(－ |$\mathit{\boldsymbol{r}}$|^0.8)，$\mathit{\boldsymbol{r}}$的表达式为

$\mathit{\boldsymbol{r}}\left( p \right) = \frac{{\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( q \right)} } \right\|_2^2}}{{\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( q \right)} } \right\|_2^2 + 0.5}}$

(4)

式中，$\mathit{\boldsymbol{r}}\left( p \right)$是$\mathit{\boldsymbol{r}}$在$p$处的一个元素值；${N_h}\left( p \right)$是以像素点$p$为中心边长为$h$的正方形窗口。0.5的作用是为了防止在图像的平坦区域出现一个过大的$\mathit{\boldsymbol{r}}\left( p \right)$值。表达式$\mathit{\boldsymbol{r}}\left( p \right)$中分子$\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( p \right)} } \right\|_2^2$的作用是计算${N_h}\left( p \right)$窗口中图像结构的强度，当像素点$p$处于纹理丰富区域则这个窗口中像素梯度值正负差异较大，则在计算求和项时${\nabla \mathit{\boldsymbol{B}}\left( p \right)}$大多数会出现相互抵消情况，这时分子$\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( p \right)} } \right\|_2^2$会比较小，求得的$\mathit{\boldsymbol{r}}\left( p \right)$较小；当像素点$p$处于平坦区域，则${N_h}\left( p \right)$窗口内像素值变化较慢其梯度方向大体一致且梯度值较小，分子$\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( p \right)} } \right\|_2^2$较小，此时$\mathit{\boldsymbol{r}}\left( p \right)$较小；如果$p$处于强边缘区域，则其梯度方向大体一致且梯度值较大，这时分子$\left\| {\sum\limits_{q \in {N_h}\left( p \right)} {\nabla \mathit{\boldsymbol{B}}\left( p \right)} } \right\|_2^2$较大，$\mathit{\boldsymbol{r}}\left( p \right)$较大。根据以上分析可知，$\mathit{\boldsymbol{r}}$值可判断像素点是否是纹理丰富区域或者平坦区域，较大的$\mathit{\boldsymbol{r}}\left( p \right)$值表示$p$所在窗口内的梯度方向大体一致，该像素点属于强边缘区域的可能性较大，相应的$\mathit{\boldsymbol{\kappa }}\left( p \right)$=exp(－|$\mathit{\boldsymbol{r}}\left( p \right)$| ^0.8)就比较小，此时对属于该强边缘区域的像素点惩罚较小；相反，$\mathit{\boldsymbol{r}}\left( p \right)$值较小则表示该像素所在窗口是一个纹理丰富区域或者平坦区域，这时$\mathit{\boldsymbol{\kappa }}\left( p \right)$的值会比较大。综上所述，上述优化模型对处于纹理丰富区域或平坦区域的像素点的正则参化数较大，对该区域惩罚力度较大，优化过程中会平滑掉该区域，而保留强边缘区域的像素点。

直接对模型式(3) 优化$\mathit{\boldsymbol{x}}$、$\mathit{\boldsymbol{k}}$较为困难，通常采用交替迭代的策略进行优化。因此优化问题式(3) 等价于两个交替优化的子问题，即

$\mathop {\min }\limits_x \frac{1}{2}\left\| {\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \gamma \mathit{\boldsymbol{\kappa }} \circ {\left\| \mathit{\boldsymbol{x}} \right\|_0}$

(5)

$\begin{array}{l} \mathop {\min }\limits_k \frac{1}{2}\left\| {\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\eta _1}{\left\| \mathit{\boldsymbol{k}} \right\|_0} + {\eta _2}\left\| {\nabla \mathit{\boldsymbol{k}}} \right\|_2^2\\ \quad \quad \quad \quad \quad {\rm{s}}{\rm{.t}}.\;\sum\limits_i {{k_i} = 1} \end{array}$

(6)

在模型式(5) 中固定$\mathit{\boldsymbol{k}}$优化$\mathit{\boldsymbol{x}}$，求解得到$\mathit{\boldsymbol{x}}$后固定$\mathit{\boldsymbol{x}}$再优化模型式(6)，交替优化模型式(5) 和模型式(6) 直到收敛，即求得了模型式(3) 的解。

式(5) 含l₀范数，是一个高度非凸最优化问题，求解l₀范数是一个NP难问题。为此采用增广拉格朗日方法，引入辅助变量$\mathit{\boldsymbol{w}}$将原问题转换为多个易求的子问题进行求解。先将问题式(5) 转换为等式约束的最小化问题，即

$\begin{array}{l} \mathop {\min }\limits_{x,w} \frac{1}{2}\left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \gamma \mathit{\boldsymbol{\kappa }} \circ {\left\| \mathit{\boldsymbol{w}} \right\|_0}\\ \quad \quad \quad \;{\rm{s}}{\rm{.t}}.\mathit{\boldsymbol{w}} = \mathit{\boldsymbol{x}} \end{array}$

(7)

式中，$\mathit{\boldsymbol{x}} \otimes \mathit{\boldsymbol{k}} \equiv \mathit{\boldsymbol{Kx}},\mathit{\boldsymbol{K}} \in {{\bf{R}}^{M \times M}}$是由模糊核$\mathit{\boldsymbol{k}}$组成的卷积矩阵；$\mathit{\boldsymbol{M}}$是图像像素点个数；使用增广拉格朗日法将等式约束最小化问题转换为无约束最小化问题，即

$\begin{array}{l} \left( {{\mathit{\boldsymbol{w}}^{l + 1}},{\mathit{\boldsymbol{x}}^{l + 1}}} \right) = \mathop {\min }\limits_{x,w} \frac{1}{2}\left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + \gamma \mathit{\boldsymbol{\kappa }} \circ {\left\| \mathit{\boldsymbol{w}} \right\|_0} + \\ \quad \quad \quad \quad {\left( {\mathit{\boldsymbol{\mu }}_x^l} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right) + \frac{{\gamma _x^l}}{2}\left\| {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right\|_2^2 \end{array}$

(8)

式中，$\gamma _x^l$是第$l$次迭代中拉格朗日惩罚参数，${\mathit{\boldsymbol{w}}^{l + 1}}$、${\mathit{\boldsymbol{x}}^{l + 1}}$分别是第$l$+1次迭代的$\mathit{\boldsymbol{w}}$和$\mathit{\boldsymbol{x}}$。$\mathit{\boldsymbol{\mu }}_x^l$是第$l$次迭代中等式约束条件$\mathit{\boldsymbol{w}}$= $\mathit{\boldsymbol{x}}$的拉格朗日乘子，它的更新公式为

$\mathit{\boldsymbol{\mu }}_x^{l + 1} = \mathit{\boldsymbol{\mu }}_x^l + \gamma _x^l\left( {{\mathit{\boldsymbol{x}}^{l + 1}} - {\mathit{\boldsymbol{w}}^{l + 1}}} \right)$

(9)

原则上拉格朗日惩罚参数$\gamma _x^0$初始值是很小的正数，然后按照$\gamma _x^{l + 1} = {\rho _x} \times \gamma _x^l$方式递增，式中更新系数${\rho _x}$>1促使随着迭代的进行$\mathit{\boldsymbol{w}}$将越来越逼近$\mathit{\boldsymbol{x}}$。

式(8) 可转换为两个优化问题间的交替迭代，即

${\mathit{\boldsymbol{w}}^{l + 1}} = \mathop {\min }\limits_w \left[ {\begin{array}{*{20}{l}} {\frac{{\gamma _x^l}}{2}\left\| {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right\|_2^2 + }\\ {{{\left( {\mathit{\boldsymbol{\mu }}_x^l} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right) + }\\ {\gamma \mathit{\boldsymbol{\kappa }} \circ {{\left\| \mathit{\boldsymbol{w}} \right\|}_0}} \end{array}} \right]$

(10)

${\mathit{\boldsymbol{x}}^{l + 1}} = \mathop {\min }\limits_x \left[ {\begin{array}{*{20}{l}} {\frac{1}{2}\left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + }\\ {{{\left( {\mathit{\boldsymbol{\mu }}_x^l} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right) + }\\ {\frac{{\gamma _x^l}}{2}\left\| {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{w}}} \right\|_2^2} \end{array}} \right]$

(11)

式(10) 经整理后得

${\mathit{\boldsymbol{w}}^{l + 1}} = \mathop {\min }\limits_w \left[ {\begin{array}{*{20}{l}} {\frac{{\gamma _x^l}}{2}\left\| {\mathit{\boldsymbol{w}} - \left( {{\mathit{\boldsymbol{x}}^l} + \mathit{\boldsymbol{\mu }}_x^l/\gamma _x^l} \right)} \right\|_2^2 + }\\ {\gamma \mathit{\boldsymbol{\kappa }} \circ {{\left\| \mathit{\boldsymbol{w}} \right\|}_0}} \end{array}} \right]$

(12)

式(12) 可以通过硬阈值方法求解，即

${\mathit{\boldsymbol{w}}^{l + 1}} = {\Theta _{{\rm{Hard}}}}\left( {{\mathit{\boldsymbol{x}}^l} + \mathit{\boldsymbol{\mu }}_x^l/\gamma _x^l,{{\left( {2\gamma \mathit{\boldsymbol{\kappa }}/\gamma _x^l} \right)}^{\frac{1}{2}}}} \right)$

(13)

式中，Θ_Hard定义为Θ_Hard($a$，$b$)=$a$×($a$ >$b$)。式(11) 的解可以对其求导并令导数等于${\bf{0}}$得到，即

$\begin{array}{l} \quad \left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \gamma _x^l\mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{x}}_x^l = \\ {\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{y + }}\gamma _x^l\left( {\mathit{\boldsymbol{w}} - \mathit{\boldsymbol{\mu }}_x^l/\gamma _x^l} \right) \end{array}$

(14)

式中，$\mathit{\boldsymbol{I}}$表示单位矩阵，通过快速傅里叶变换在频域求解，即

${\mathit{\boldsymbol{x}}^{l + 1}} = {F^{ - 1}}\left( {\frac{{\overline {F\left( \mathit{\boldsymbol{k}} \right)} F\left( \mathit{\boldsymbol{y}} \right) + \gamma _x^lF\left( {\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{\mu }}_x}/\gamma _x^l} \right)}}{{\overline {F\left( \mathit{\boldsymbol{k}} \right)} F\left( \mathit{\boldsymbol{k}} \right) + \gamma _x^l}}} \right)$

(15)

式中，$F\left( \circ \right)$和${F^{ - 1}}\left( \circ \right)$分别表示傅里叶正变换和傅里叶反变换，$\overline {F\left( \circ \right)} $表示取共轭操作。交替迭代求解式(10) 和式(11) 直到$\mathit{\boldsymbol{x}}$收敛到局部最佳值或者迭代达到最大的迭代次数。然后，固定$\mathit{\boldsymbol{x}}$转而优化式(6) 估计模糊核$\mathit{\boldsymbol{k}}$。中间清晰图像的求解过程可总结为：

输入：模糊图像$\mathit{\boldsymbol{B}}$和模糊核$\mathit{\boldsymbol{k}}$并设置$\mathit{\boldsymbol{\mu }}_x^0 = {\bf{0}},l = 0,\gamma = 1{{\rm{E}}^{ - 2}},\gamma _x^0 = 1{{\rm{E}}^{ - 2}},{\max _{{\gamma _x}}} = {2^8},{\mathit{\boldsymbol{x}}^0} = \mathit{\boldsymbol{y}},{\rho _x} = 2.25$迭代：

1) 通过计算式(13) 更新辅助变量${\mathit{\boldsymbol{w}}^{l + 1}}$;

2) 通过计算式(15) 更新中间清晰图像${\mathit{\boldsymbol{x}}^{l + 1}}$;

3) 更新拉格朗日乘子

$\mathit{\boldsymbol{\mu }}_x^{l + 1} = \mathit{\boldsymbol{\mu }}_x^l + \gamma _x^l\left( {{\mathit{\boldsymbol{x}}^{l + 1}} - {\mathit{\boldsymbol{w}}^{l + 1}}} \right)$

4) 更新惩罚参数

$\gamma _x^{l + 1} = {\rho _x} \times \gamma _x^l$

5) 检查算法收敛条件：

条件1：${\left\| {{\mathit{\boldsymbol{x}}^{l + 1}} \otimes \mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_\infty } ＜ \varepsilon $与

${\left\| {{\mathit{\boldsymbol{w}}^{l + 1}} - {\mathit{\boldsymbol{x}}^{l + 1}}} \right\|_\infty }＜\varepsilon $

条件2：$\gamma _x^{l + 1} \ge {\max _{{\gamma _x}}}$

6) 迭代次数加1

$l = l + 1$

直到：满足收敛条件。

输出：中间梯度图像$\mathit{\boldsymbol{x}}$。

上述过程中的2个收敛条件，只要满足条件1或条件2中的任何一个，算法都将停止迭代，收敛条件1中l_∞范数表示向量中各元素分量绝对值的最大值。

由式(5) 估计获得清晰边缘图像$\mathit{\boldsymbol{x}}$后，可用于指导模糊核$\mathit{\boldsymbol{k}}$的估计。式(2) 中的模糊核$\mathit{\boldsymbol{k}}$可以通过求解式(6) 最小化问题来获得。同样式(6) 也是个NP难问题。采用半二次函数法求解，引入辅助变量$\mathit{\boldsymbol{g}}$将约束最小化问题转换为无约束最小化问题，即

$\begin{array}{l} \quad \mathop {\min }\limits_{k,g} \frac{1}{2}\sum\limits_{d \in \Lambda } {\left\| {{\mathit{\boldsymbol{X}}_d}\mathit{\boldsymbol{k}} - {\mathit{\boldsymbol{y}}_d}} \right\|_2^2} + {\eta _1}{\left\| \mathit{\boldsymbol{g}} \right\|_0} + \\ {\eta _2}\left\| {\nabla \mathit{\boldsymbol{k}}} \right\|_2^2 + \frac{{{\delta _1}}}{2}\left\| {\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \frac{{{\delta _2}}}{2}\left\| {{{\bf{1}}^{\rm{T}}}\mathit{\boldsymbol{k}} - 1} \right\|_2^2 \end{array}$

(16)

式中，$d \in {\bf{\Lambda }} \buildrel \Delta \over = \left\{ {h,v} \right\}$，${\mathit{\boldsymbol{X}}_h}$，${\mathit{\boldsymbol{X}}_v}$分别表示估计得到的水平梯度图像和垂直梯度图像形成的卷积矩阵；${\mathit{\boldsymbol{y}}_h}$，${\mathit{\boldsymbol{y}}_v}$分别表示${\nabla _h}\mathit{\boldsymbol{B}}$和${\nabla _v}\mathit{\boldsymbol{B}}$的向量形式；${\bf{1}}$表示元素全为1的列向量；${\eta _1}$, ${\eta _2}$, ${\delta _1}$, ${\delta _2}$是正则参数。

采用交替迭代的方式求解式(16)。首先将其转化为两个易求的子问题，即

${\mathit{\boldsymbol{g}}^{l + 1}} = \mathop {\min }\limits_g \frac{{{\delta _1}}}{2}\left\| {\mathit{\boldsymbol{g}} - {\mathit{\boldsymbol{k}}^l}} \right\|_2^2 + {\eta _1}{\left\| \mathit{\boldsymbol{g}} \right\|_0}$

(17)

$\begin{array}{l} {\mathit{\boldsymbol{k}}^{l + 1}} = \mathop {\min }\limits_k \left[ {\frac{1}{2}} \right.\sum\limits_{d \in \Lambda } {\left\| {{\mathit{\boldsymbol{X}}_d}\mathit{\boldsymbol{k}} - {\mathit{\boldsymbol{y}}_d}} \right\|_2^2} + {\eta _2}\left\| {\nabla \mathit{\boldsymbol{k}}} \right\|_2^2 + \\ \quad \quad \quad \left. {\frac{{{\delta _1}}}{2}\left\| {\mathit{\boldsymbol{k}} - {\mathit{\boldsymbol{g}}^{l + 1}}} \right\|_2^2 + \frac{{{\delta _2}}}{2}\left\| {{{\bf{1}}^{\rm{T}}}\mathit{\boldsymbol{k}} - 1} \right\|_2^2} \right] \end{array}$

(18)

式中，${\mathit{\boldsymbol{g}}^{l + 1}}$，${\mathit{\boldsymbol{k}}^{l + 1}}$分别是第$l$+1次迭代的$\mathit{\boldsymbol{g}}$和$\mathit{\boldsymbol{k}}$。

式(17) 可以使用硬阈值方法求得，即

${\mathit{\boldsymbol{g}}^{l + 1}} = {\Theta _{{\rm{Hard}}}}\left( {{\mathit{\boldsymbol{k}}^l},{{\left( {2{\eta _1}/{\delta _1}} \right)}^{\frac{1}{2}}}} \right)$

(19)

式(18) 是一个二次函数优化问题，它具有闭合形式的解，即

$\begin{array}{l} \left( {\sum\limits_{d \in \Lambda } {\mathit{\boldsymbol{X}}_d^{\rm{T}}{\mathit{\boldsymbol{X}}_d} + 2{\eta _2}{\nabla ^{\rm{T}}}\nabla + {\delta _2}\mathit{\pmb{11}}^ {\rm{T}} + {\delta _1}\mathit{\boldsymbol{I}}} } \right){\mathit{\boldsymbol{k}}^{l + 1}} = \\ \quad \quad \quad \quad \quad \sum\limits_{d \in \Lambda } {\mathit{\boldsymbol{X}}_d^{\rm{T}}{\mathit{\boldsymbol{y}}_d}} + {\delta _2}{\bf{1}} + {\delta _1}{\mathit{\boldsymbol{g}}^{l + 1}} \end{array}$

(20)

$\begin{array}{l} {\mathit{\boldsymbol{k}}^{l + 1}} = {\left( {\sum\limits_{d \in \Lambda } {\mathit{\boldsymbol{X}}_d^{\rm{T}}{\mathit{\boldsymbol{X}}_d} + 2{\eta _2}{\nabla ^{\rm{T}}}\nabla + {\delta _2}\mathit{\pmb{11}}^ {\rm{T}} + {\delta _1}\mathit{\boldsymbol{I}}} } \right)^{ - 1}} \times \\ \quad \quad \quad \quad \quad \left( {\sum\limits_{d \in \Lambda } {\mathit{\boldsymbol{X}}_d^{\rm{T}}{\mathit{\boldsymbol{y}}_d}} + {\delta _2}{\bf{1}} + {\delta _1}{\mathit{\boldsymbol{g}}^{l + 1}}} \right) \end{array}$

(21)

式中，${\bf{1}}$是$M$×1的向量，$\mathit{\pmb{11}}^ {\rm{T}}$是$M$×$M$的分块循环(BCCB)矩阵，它与向量相乘可以利用分块循环矩阵可对角化性质，类似于式(14) 可使用快速傅里叶变换在频域快速求得。交替求解式(17) 和式(18) 直到$\mathit{\boldsymbol{k}}$收敛到局部最优值。理论上当${\delta _1}$、${\delta _2}$趋向于无穷时式(13) 的解就是原问题的解，然而本文在实现算法时采取${\delta _1}$逐渐增大^[13]，达到一个较大的数值${{\max }_{{{\delta }_{1}}}}$后迭代停止，${\delta _2}$则简单的取定值。模糊核求解过程可总结为：

输入：模糊图像$\mathit{\boldsymbol{B}}$和中间清晰图像并设置${\eta _1}$=5E^－4，${\eta _2}$=1E^－2，${\delta _1}$=1E^－3，${\delta _2}$=1E^－5，${{\max }_{{{\delta }_{1}}}}$=1E⁺⁵，模糊核${{\mathit{\boldsymbol{k}}}^{0}}$。

迭代：

1) 通过计算式(19) 更新辅助变量${\mathit{\boldsymbol{g}}^{l + 1}}$；

2) 通过计算式(21) 更新中间清晰图像${\mathit{\boldsymbol{k}}^{l + 1}}$；

3) 更新惩罚参数${\delta _1}$：${\delta _1}$=2×${\delta _1}$；

4) 检查算法收敛条件：${{\delta }_{1}}\ge {{\max }_{{{\delta }_{1}}}}$；

5) 迭代次数加1：$l\text{=}l+1$；

直到：满足收敛条件。

输出：模糊核$\mathit{\boldsymbol{k}}$。

上述过程中第4) 步收敛条件成立则停止迭代并输出模糊核。

获得模糊核$\mathit{\boldsymbol{k}}$后采用非盲去模糊算法恢复最终的清晰图像。Krishnan等人^[14]采用超拉普拉斯先验约束清晰图像，获得的清晰图像中保留图像的主要细节并且有效地抑制振铃效应。使用超拉普拉斯先验对模糊图像进行非盲去模糊会平滑掉一些细小纹理，使得估计图像的局部区域太过平滑；而使用全变差^[15]作为正则约束项，估计获得的图像能保留图像的细小纹理但边界往往出现噪声和振铃。本文结合上述两种算法，采用优势互补策略^[16]，分别使用上述两种先验估计两幅清晰图像${{\mathit{\boldsymbol{L}}}_{1}}$、${{\mathit{\boldsymbol{L}}}_{2}}$，然后取两者平均得到最终的清晰图像。

${{\mathit{\boldsymbol{L}}}_{1}}$可由最小化超拉普拉斯先验约束的能量泛函式获得，即

$\mathop {\min }\limits_{{L_1}} \left\| {{\mathit{\boldsymbol{L}}_1} \otimes \mathit{\boldsymbol{k}} - {\mathit{\boldsymbol{I}}_s}} \right\|_2^2 + \rho {\left\| {\nabla {\mathit{\boldsymbol{L}}_1}} \right\|^\alpha }$

(22)

式中，${{\mathit{\boldsymbol{I}}}_{s}}$表示模糊图像使用冲击滤波器增强图像后的图像。0.5≤$\alpha $≤0.8，$\rho $是权重。$\alpha $越小则估计获得的图像越平滑，可以通过调整$\alpha $来得到令人满意复原效果。本文算法实现时$\alpha $取0.8。

${{\mathit{\boldsymbol{L}}}_{2}}$可以的计算公式为

$\mathop {\min }\limits_{{L_2}} \left\| {{\mathit{\boldsymbol{L}}_2} \otimes \mathit{\boldsymbol{k}} - {\mathit{\boldsymbol{I}}_s}} \right\|_2^2 + \mu {\left\| {\nabla {\mathit{\boldsymbol{L}}_2}} \right\|_{TV}}$

(23)

${\left\| {\nabla {\mathit{\boldsymbol{L}}_2}} \right\|_{TV}} = \sqrt {{{\left( {{\partial _x}{\mathit{\boldsymbol{L}}_2}} \right)}^2} + {{\left( {{\partial _y}{\mathit{\boldsymbol{L}}_2}} \right)}^2}} $

(24)

式中，$\mu $是正则参数，${{\left\| \nabla {{\mathit{\boldsymbol{L}}}_{2}} \right\|}_{TV}}$是各向同性全变分(isotropic total variation)范数。由于模糊图像尤其是模糊程度较大的图像丢失了大量的细节，这不利于图像的重建。因此，有必要使用冲击滤波器增强有意义的边缘，然后求解式(23) 得到清晰图像。

最终的清晰图像$\mathit{\boldsymbol{L}}$计算公式为

$\mathit{\boldsymbol{L}} = \left( {{\mathit{\boldsymbol{L}}_1} + {\mathit{\boldsymbol{L}}_2}} \right)/2$

(25)

为了检验本文算法的有效性，对Levin^[17]测试集和实际拍摄的模糊图像分别进行仿真测试，并同现有算法进行比较。为确保实验的公平性，实验比较的各种方法的参数按照原作者论文描述设置，在Levin数据集上计算各算法的盲反卷积误差累计直方图时，采用与Zuo等人^[8]相同的非盲去模糊算法，其余实验所有方法采用的非盲去模糊算法均采用第3部分中描述的非盲去模糊方法。

Levin^[17]的测试集包含32幅模糊图像，由4幅清晰图像分别与8个尺寸范围从13×13像素到27×27像素大小不同的模糊核卷积产生的。图 1给出了这4幅清晰图像和8个真实的模糊核。

图 2给出了Levin^[17]测试图像的盲去模糊误差比值的累积直方图。其中，盲去模糊误差比值定义^[17]为

$RDE = \frac{{\sum\limits_{a = 1}^A {\sum\limits_{b = 1}^B {{{\left\| {u\left( {a,b} \right) - \bar u\left( {a,b} \right)} \right\|}^2}} } }}{{\sum\limits_{a = 1}^A {\sum\limits_{b = 1}^B {{{\left\| {u\left( {a,b} \right) - {u_r}\left( {a,b} \right)} \right\|}^2}} } }}$

(26)

式中，$A\times B$表示图像大小，$u\left( a,b \right)$表示原始清晰图像在$\left( a,b \right)$处的像素值，$\bar{u}\left( a,b \right)$是借助于估计得到的模糊核采用非盲去模糊算法复原得到的像素值，${{u}_{r}}\left( a,b \right)$表示采用真实模糊核采用反卷积算法恢复的像素值。$RDE$的值越小则盲反卷积的效果越好。对于盲去模糊是否成功，文献[17]指出当盲去模糊误差比值$RDE$小于3时则认为盲去模糊是成功的，即最终非盲反卷积得到的清晰图像是可以接受的。由图 2可知，本文方法在对32幅模糊图像进行盲反卷积其误差比值$RDE$不大于2.6，均小于3，即本文算法在Levin测试集上的盲反卷积成功率为100 %。

此外，图 3给出了本文算法与其他盲去模糊算法在Levin^[17]的测试集图像im02上的8幅模糊图像上估计得到的模糊核。图 3中从上到下依次对应真实的模糊核，Fergus算法^[1]，Cho和Lee算法^[18]，Levin算法^[12]，Xu和Jia算法^[5]，Zuo算法^[8]，Pan算法^[10]以及本文算法估计得到的模糊核。对比图 3中的真实模糊核发现，本文算法估计的模糊核与其他算法估计出来的模糊核相比拥有更准确的支撑和较少的孤立噪点。

图 4针对测试图像-04和模糊核-04生成的模糊图像，给出了本文算法与其他几种盲去模糊算法的最终反卷积得到的清晰图像及估计出的模糊核。为了从客观上验证算法的有效性，同时计算了反卷积得到的清晰图像的峰值信噪(PSNR)。图 4中从左到右从上到下依次为模糊图像，已知真实模糊核下的使用本文图像非盲去模糊算法得到的图像，Fergus算法^[1]，Pan算法^[10]，Levin算法^[12]，Krishnan算法^[3]，Zuo算法^[8]，Xu算法^[5]以及本文算法。对比真实图像及对应的模糊核可知，Fergus算法^[1]不能很好地从大尺度的模糊图像估计出模糊核，针对im04-ker04的模糊图像的盲反卷积复原图像完全失败。对比真实模糊核-04，Krishnan算法^[3]，Levin算法^[12]，估计的模糊核与真实的模糊核之间存在比较大偏差；Pan算法^[10]获得的模糊核存在明显的拖尾，而本文方法估计的模糊核更接近于真实模糊核。从视觉效果上，本文算法以及文献[5, 8, 12]的最终复原清晰图像优于文献[1, 3, 10]。从客观评价来看，本文算法估计出清晰图像的PSNR为30.59 dB，在对比算法中具有最高的PSNR。

为了进一步验证本文算法的有效性，采用实际彩色模糊图像作为实验数据，并与文献[1, 3-5, 8, 10, 18-19]所提算法进行比较。获得模糊核后，将使用式(25) 得到最终的清晰图像。图 5—图 7列出了不同模糊程度的模糊图像及各种算法估计获得的模糊核和最终反卷积得到的清晰图像。由图 5—图 7可知，本文算法可以较好地去除实际彩色图像的运动模糊，具有较高的视觉质量，这得益于模糊核的准确估计。结合清晰图像和估计的模糊核发现，本文算法估计的模糊核具有更准确的支撑和更少的噪声，这也使获得的清晰图像具有更清晰的细节以及更为锐化的边缘。除Cai^[19]算法外，3幅图像模糊核依次设置为31×31像素，55×55像素，51×51像素。由图 5图像的盲去模糊结果可知，Kotera算法^[4]估计的模糊核含有较多的噪点，导致最终的清晰图像含有大量的振铃；Fergus算法^[1]，Pan算法^[10]以及Krishnan算法^[3]估计的清晰图像在视觉效果上优于Kotera算法^[4]估计的清晰图像，但Fergus算法^[1]估计的模糊核仍有一定程度的拖尾；而Krishnan算法^[3]和Pan算法^[10]的模糊核含有明显的孤立噪点。Cai算法^[19]获得的模糊核类似于脉冲函数，其获得的“清晰图像”较为模糊；Xu和Jia^[5]获得的模糊核有少量孤立噪声存在，其清晰图像左上方出现了振铃。相比之下只有Zuo算法^[8]获得清晰图像在视觉上与本文最为接近，但Zuo算法获得的图像右上角的“小白块”处有轻微的振铃现象。由图 6和图 7模糊图像的复原实验可知，Kotera算法^[4]和Krishnan算法^[3]复原的清晰图像视觉效果较差；Cho和Lee^[18]算法也不能从wall模糊图像中较好地估计出模糊核，导致最终的清晰图像出现大量的伪迹。Zuo算法^[8]复原wall图像时也出现了伪迹，复原roma图像出现严重的重影，极大的影响了视觉效果。由比较可知，本文算法均能从3幅模糊图像中获得较好的清晰图像和模糊核。实验结果表明，本文算法可以从实际模糊图像中较准确地估计出模糊核，使得最终的清晰图像具有较高的视觉质量。

为了更好地评价本文算法，表 1给出了本文算法与其他各种盲去模糊算法在不同尺寸模糊图像上的计算时间。由表 1可知，本文算法采用l₀范数作为正则约束项进行优化，过程比较耗时，因此在计算时间上并不是最优的。后续研究将通过C+ +优化算法代码并结合GPU并行计算提升算法运行效率。

在基于正则化的图像盲复原框架下，提出了一种基于强边缘的盲去模糊算法。采用自适应加权l₀约束图像，隐式地估计出图像强边缘，用于指导模糊核估计；利用l₀约束模糊核像素和范数l₂约束模糊核梯度，确保模糊核的稀疏性和连续性；最后把模糊核归一化先验加入盲去模糊模型，并给出了最终的优化模型及相应的求解算法。实验结果表明，相比于其他算法，本文算法获得的模糊核具有更准确的支撑和较少的噪点，获得的清晰图像具有较优的视觉效果。

本文盲去模糊模型所含参数较多，实际调参比较困难，后期工作将通过数据驱动的方式为模型自动学习正则化参数，进一步提高模糊核估计精度；此外，本文算法的模糊核需要人工设定模糊大小，并不能由模糊图像估计出模糊核的尺度，因此有必要研究模糊核尺寸自适应调整的盲去模糊算法。同时，算法代码优化和并行实现也是后续的研究内容之一。