0 引言

受大气湍流、聚焦不准确以及相机与被拍摄场景之间存在相对位移等因素的影响，导致采集的图像出现模糊退化，如果假设这种退化是线性空不变的，则可用退化模型G=H⊗F+N来描述该退化过程，其中:G是模糊图像，H是模糊核，也称为点扩散函数(Point Spread Function, PSF)，Y表示二维卷积，F是清晰图像，N是加性噪声。模糊图像复原是根据采集到的模糊图像G估计清晰图像F。根据模糊核H是否已知，可将模糊图像复原分为非盲复原和盲复原两种类型。如果复原时模糊核H已知，则称为模糊图像非盲复原(解模糊)技术；如果模糊核H未知，则称为模糊图像盲复原技术。即使模糊核已知，由于存在噪声影响，模糊图像非盲复原仍然属于病态逆问题。非盲复原的经典算法包括维纳滤波^[1]、Richardson-Lucy迭代算法^[2-3]以及能较好地保持图像边缘的总变分(Total Variation, TV)复原算法^[4]等。对于盲复原技术，因模糊核未知而更加具有挑战性。

一般的模糊图像盲复原算法可分为模糊核估计和解模糊两个步骤。早期的模糊核估计算法研究中，通常假设PSF符合某个简单的参数模型，例如圆盘模糊、直线运动模糊、大气湍流模糊等，这种利用参数模型和估计模型参数以达到模糊核估计目的的方法，称为参数模型法。例如，Cannon^[5]利用运动模糊图像的频谱特点，从运动模糊图像中估计出运动模糊长度和运动模糊方向。Yitzhaky等^[6]提出利用旋转差分法检测运动模糊方向。Oliveira等^[7]利用改进的Radon变换估计散焦模糊及直线运动模糊参数。在一定条件下利用参数模型法能得到比较准确的参数估计，但此类方法对噪声比较敏感，且在模糊比较严重时参数估计的准确度会相应降低。实际拍摄中由于相机抖动或者被拍摄目标相对于相机的运动所带来的模糊远比前面所假设的这些退化模型复杂，仅仅用几个简单的参数模型远不能描述实际的退化过程。

2006年，Fergus等^[8]利用自然图像梯度的重尾分布特性，并用混合高斯模型对其拟合，而PSF采用指数混合模型进行拟合，最后采用变分贝叶斯方法求解，估计出任意形状的非规则运动模糊核。2008年，Shan等^[9]利用稀疏先验，对图像梯度采用分段函数进行拟合，从而估计出模糊核和清晰图像。Cho等^[10]首先通过Shock滤波器来预测运动模糊图像的显著边缘，再利用增强后的边缘梯度信息估计出模糊核，此方法的在计算效率上较之前面的方法有明显提高。Xu等^[11]提出两步模糊核估计算法，最后利用TV l₁模型估计清晰图像，该方法能估计出较大尺寸的模糊核。Krishnan等^[12]利用归一化稀疏先验测度估计出模糊核，并利用超拉普拉斯先验模型复原出清晰图像。

以上这些算法都是利用图像中的显著边缘信息估计模糊核，通常均能达到理想的复原效果，但当图像尺寸太小或图像中显著边缘信息过少时，则无法有效地估计出模糊核，从而导致最终复原出的图像质量不理想。另外，现有的模糊图像盲复原算法的处理流程一般都是先将彩色模糊图像转换成灰度图像，再利用图像的灰度信息估计出模糊核，最后根据估计到的模糊核对模糊图像的三个通道分别解模糊，复原出清晰图像。在彩色(三通道R、G、B)图像转换成灰度(单通道)图像时，很显然存在信息的丢失。针对上述问题，本文在张量框架下，将彩色模糊图像作为一个三阶张量，直接利用该三阶张量来估计模糊核(而不是将其转换成灰度图像后，利用灰度信息来估计模糊核)，探讨了一种基于张量总变分的模糊图像盲复原算法。

1 张量总变分

关于张量的讨论，最早可以追溯到1927年^[13]。2009年，Kolda等^[14]给出了一个关于高阶张量分解及应用的一个综述。2008年，Kilmer等^[15]定义了一个新的张量积，并在这个新的张量积的基础上提出一个新张量框架。2010年，Braman^[16]在Kilmer的基础上给出了系统的描述和证明。Kilmer等^[17-18]在这个新的张量框架下通过构造一个张量预决算子将LSQR (Least Squares QR-factorization)和共轭梯度法扩展到张量框架中，并将其应用到图像非盲复原中。需要指出的是，文献[18]只讨论了在模糊核已知时的图像解模糊算法，属于模糊图像非盲复原问题。本文利用三阶张量来表示彩色图像，拟在Kilmer提出的新张量框架下，探讨一种基于张量总变分的模糊图像盲复原算法。

1.1 新张量框架

下面给出新张量框架中的一些基本定义，详细的介绍可以参见文献[18]。

定义1 一个三阶张量Α∈R^m×p×n的第i个侧切片为Α_i∈R^m×n，第j个前切片为Α^(j)∈R^m×p。

定义2 一个三阶张量Α∈R^m×p×n和Β∈R^p×q×n的t-product积Α*Β是一个m×q×n大小的张量：

$ \boldsymbol{A} * \boldsymbol{B} = {\rm{fold}}({\rm{bcric}}(\boldsymbol{A}) \times {\rm{unfold}}(\boldsymbol{B})) $

(1)

其中bcirc (Α)是一个mn×pn大小的块循环矩阵:

$ {\rm{bcirc(}}\boldsymbol{A}{\rm{)}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{A}^{{\rm{(1)}}}}}&{{\boldsymbol{A}^{{\rm{(}}n{\rm{)}}}}}&{{\boldsymbol{A}^{{\rm{(}}n - 1{\rm{)}}}}}& \cdots &{{\boldsymbol{A}^{{\rm{(2)}}}}}\\ {{\boldsymbol{A}^{{\rm{(2)}}}}}&{{\boldsymbol{A}^{{\rm{(}}1{\rm{)}}}}}&{{\boldsymbol{A}^{{\rm{(}}n{\rm{)}}}}}& \cdots &{{\boldsymbol{A}^{{\rm{(3)}}}}}\\ \vdots&\vdots&\vdots &{}& \vdots \\ {{\boldsymbol{A}^{{\rm{(}}n{\rm{)}}}}}&{{\boldsymbol{A}^{{\rm{(}}n - 1{\rm{)}}}}}&{{\boldsymbol{A}^{{\rm{(}}n{\rm{ - 2)}}}}}& \cdots &{{\boldsymbol{A}^{{\rm{(1)}}}}} \end{array}} \right] $

(2)

unfold (Α)是一个mn×p大小的矩阵，fold是unfold的逆运算:

$ \left\{ \begin{array}{l} {\rm{unfold(}}\boldsymbol{A}{\rm{)}} = {\left[ {{\boldsymbol{A}^{\left( 1 \right)}}\;\;{\boldsymbol{A}^{\left( 2 \right)}}\;\; \cdots {\boldsymbol{A}^{\left( n \right)}}} \right]^{\rm{T}}}\\ {\rm{fold}}\left( {{\rm{unfold(}}\boldsymbol{A}{\rm{)}}} \right) = \boldsymbol{A} \end{array} \right. $

(3)

定义3 张量Α∈R^m×p×n的转置Α^T是一个p×m×n大小的张量:

$ \begin{array}{l} {\boldsymbol{A}^{\rm{T}}} = {\rm{fold}}\left( {{{\left[ {{{\left( {{\boldsymbol{A}^{\left( 1 \right)}}} \right)}^{\rm{T}}}\;\;{{\left( {{\boldsymbol{A}^{\left( 1 \right)}}} \right)}^{\rm{T}}}\;\;\; \cdots {{\left( {{\boldsymbol{A}^{\left( 1 \right)}}} \right)}^{\rm{T}}}} \right]}^{\rm{T}}}} \right) \end{array} $

(4)

定义4 一个n×n×p大小的单位张量Ι的第一个前切片是n×n大小的单位矩阵，其余前切片是全零矩阵。

定义5 对于张量Α∈R^n×n×p，若可以找到一个张量Β∈R^n×n×p使得Α*Β=Ι且Β*Α=Ι，则称张量Β为张量Α的逆。

1.2 张量总变分模型

1994年，Rudin等^[4]将总变分正则项引入到模糊图像复原中，复原出的图像能较好地保持图像的边缘，因此基于总变分模型的解模糊算法得到广泛的应用和推广。为了表示的简洁性，假设原始图像是一个m×n的灰度图像。从线性空间映射的角度来看，可将模糊退化用如下方程来描述：

$ \boldsymbol{b} = \boldsymbol{Ax} + \boldsymbol{e} $

(5)

其中：A∈R^mn×mn是模糊算子，x、b∈R^mn分别是原始灰度图像X和观测图像B的向量形式，e∈R^mn是加性噪声。通常，基于向量框架下的总变分模型可写成：

$ \mathop {\min }\limits_x \frac{\mu }{2}||\boldsymbol{Ax} - \boldsymbol{b}||_2^2 + \sum\limits_i {||{\boldsymbol{D}_i}\boldsymbol{x}|{|_2}} $

(6)

其中:μ是正则化参数，D_ix∈R²是图像x在像素i处的离散梯度，$ $是图像x的离散总变分。利用新张量框架中定义的t-product积，可将向量框架下的总变分模型扩展成张量框架下的张量总变分模型，即t-TV (tensor-Total Variation)模型^[19]

$ \mathop {\min }\limits_X \frac{\mu }{2}\left\| {\boldsymbol{A} * \boldsymbol{X} - \boldsymbol{B}} \right\|_2^2 + \sum\limits_i {{{\left\| {{\boldsymbol{D}_i} * \boldsymbol{X}} \right\|}_2}} $

(7)

其中:Α∈R^m×m×n是模糊算子张量，X∈R^m×c×n是原始灰度图像张量(c=1)或彩色图像张量(c=3), Β∈R^m×c×n是模糊图像张量，D_i∈R^m×m×n是差分算子张量。如果将图像的三个通道分别记为I₁，I₂，I₃，则可利用twist操作将彩色图像的三个通道的矩阵转换成一个三阶张量X，其中X(:, j, :)=twist (I_j))。

利用向量框架下的总变分模型处理彩色图像时，一般是将彩色图像分成三个单独的通道(R、G、B)独立完成，而张量总变分模型将彩色图像作为三阶张量，同时利用彩色图像三个通道的总变分对复原图像施加约束，从而复原出清晰图像。基于张量总变分模型图像解模糊算法与文献[18]中基于张量共轭梯度算法解模糊相比，复原图像的边缘更为清晰，振铃较少^[19]。

2 基于张量总变分的盲复原算法

一般的模糊图像盲复原算法可分两步完成：模糊核估计和图像解模糊。为了保证算法的有效性和高效性，可在一个coarse-to-fine (图像分辨率从低到高)的多尺度框架下，通过交替迭代估计得到模糊核。在此多尺度框架下，通过改变图像分辨率的方式，在不同尺度下估计模糊核。在最粗糙的级别上，模糊核的大小一般是3×3大小，将模糊图像B下采样作为估计的初始值，交替迭代求解出清晰图像X，并利用双线性插值将其上采后作为下一级的清晰图像的初始估计。估计到原始图像分辨率上的模糊核后，再用基于t-TV的解模糊算法来复原出清晰图像。基于张量总变分的图像盲复原算法的流程如图 1所示。

2.1 显著边缘信息提取

一幅图像可看成是由结构信息和纹理信息来构成。一般来说，结构信息对应图像中较大的目标，而纹理部分则对应图像中的细节信息，Aujol利用总变分模型来抽取图像中的结构信息^[20]，

$ \mathop {\min }\limits_x \frac{\mu }{2}||\boldsymbol{x} - \boldsymbol{b}||_2^2 + \sum\limits_i {||{\boldsymbol{D}_i}\boldsymbol{x}|{|_p}} $

(8)

通过调节参数μ可以得到图像在不同尺度上的结构信息。受此思想的启发，本文在张量框架下利用张量总变分抽取彩色图像中不同尺度下的主要结构信息(包含了显著边缘)，利用此部分信息来估计模糊核。张量框架下的显著边缘提取的目标函数为：

$ \mathop {\min }\limits_{{X_s}} \frac{\mu }{2}||{\boldsymbol{X}_s} - \boldsymbol{X}||_2^2 + \sum\limits_i {||{\boldsymbol{D}_i}{\boldsymbol{X}_s}|{|_2}} $

(9)

其中X_s∈R^m×3×n，X∈R^m×3×n。图 2中给出了不同尺度上的结构信息图像。得到包含显著边缘的图像后，再利用Shock滤波器^[21]对边缘进行增强。

2.2 模糊核估计

从模糊退化模型可以看出，图像中平缓区域对模糊核估计是无用的，而且在模糊核估计的初期迭代中估计到的模糊核与真实模糊核的误差很大，从而导致复原图像中有较强的振铃，而振铃通常出现在边缘附近灰度变化平缓的区域，因此本文算法选择在图像的梯度域中估计模糊核，在张量框架下建立相应的目标函数：

$ \min \left\| {\nabla {\boldsymbol{X}_s} * \boldsymbol{K} - \nabla \boldsymbol{B}} \right\|_2^2 + \gamma \left\| \boldsymbol{K} \right\|_2^2 $

(10)

其中:∇X_s∈R^3m×m×n，∇B∈R^3m×1×n，K∈R^m×1×n，∇表示差分运算。直接对式(10)求导，并令其等于0，即可得:

$ \boldsymbol{K} = \frac{{{\nabla _h}\boldsymbol{X}_s^{\rm{T}} * {\nabla _h}\boldsymbol{B} + {\nabla _v}\boldsymbol{X}_s^{\rm{T}} * {\nabla _v}\boldsymbol{B}}}{{{\nabla _h}\boldsymbol{X}_s^{\rm{T}} * {\nabla _h}{\boldsymbol{X}_s} + {\nabla _v}\boldsymbol{X}_s^{\rm{T}} * {\nabla _v}{\boldsymbol{X}_s} + \gamma }} $

(11)

其中∇_h和∇_v分别是水平方向和垂直方向的差分算子。

2.3 图像解模糊

利用新张量框架中的t-product积建立的张量总变分模型如式(7)所示。利用Wang等^[22]交替迭代最小化的思想，在张量框架下进行求解^[19]。首先，引入一个辅助张量W_i作为D_iX的近似，并引入W_i和D_iX的罚函数，式(7)被改写为：

$ \begin{array}{l} \mathop {\min }\limits_{X,W} \frac{\mu }{2}\left\| {\boldsymbol{A} * \boldsymbol{X} - \boldsymbol{B}} \right\|_2^2 + \sum\limits_i {{{\left\| {{\boldsymbol{W}_i}} \right\|}_2}} + \\ \;\;\;\;\;\;\frac{\beta }{2}\sum\limits_i {\left\| {{\boldsymbol{W}_i} - {\boldsymbol{D}_i} * \boldsymbol{X}} \right\|_2^2} \end{array} $

(12)

随着β→∞，式(12)的解逼近于式(7)的解。固定X求解W_i，则解为：

$ {\boldsymbol{W}_i} = \max \left\{ {||{\boldsymbol{D}_i} * \boldsymbol{X}|| - \frac{1}{\beta },0} \right\}\frac{{{\boldsymbol{D}_i}*\boldsymbol{X}}}{{||{\boldsymbol{D}_i}*\boldsymbol{X}||}} $

(13)

固定W_i，得到关于X的目标函数：

$ \mathop {\min }\limits_X \frac{\mu }{2}\left\| {\boldsymbol{A} * \boldsymbol{X} - \boldsymbol{B}} \right\|_2^2 + \frac{\beta }{2}\sum\limits_i {\left\| {{\boldsymbol{W}_i} - {\boldsymbol{D}_i} * \boldsymbol{X}} \right\|_2^2} $

(14)

直接对式(14)求导可得：

$ X = \frac{{(\mu /\beta ){\boldsymbol{A}^{\rm{T}}} * \boldsymbol{B} + \boldsymbol{D}_{\rm{v}}^{\rm{T}} * {\boldsymbol{W}_v} + \boldsymbol{D}_{\rm{h}}^{\rm{T}} * {\boldsymbol{W}_{\rm{h}}}}}{{(\mu /\beta ){\boldsymbol{A}^{\rm{T}}} * \boldsymbol{A} + \boldsymbol{D}_{\rm{v}}^{\rm{T}} * {\boldsymbol{D}_v} + \boldsymbol{D}_{\rm{h}}^{\rm{T}} * {\boldsymbol{D}_{\rm{h}}}}} $

(15)

3 实验分析

实验分为仿真模糊图像实验和真实模糊图像实验两部分。为验证算法的有效性，实验中选择文献[8]和文献[11]的算法进行对比，这两种算法是模糊图像盲复原算法中具有代表性的和复原效果较好的两种算法，作者均给出了算法的源代码或可执行程序，从而保证了算法对比的公正性。实验中显著边缘提取时μ可从0.005递增，从而得到不同尺度的边缘。模糊核估计中γ设为2。在最后图像解模糊时，μ的值为2 500。

3.1 仿真实验

实验中选择6张彩色清晰原图^[23]和4个模糊核^[24]，如图 3(从左至右，从上至下，第2张图为255×255大小，其余为401×401大小)和图 4所示，生成24张模糊图像。

文献[8]、文献[11]及本文算法的复原图像在平均峰值信噪比(Peak Signal-to-Noise Ratio，PSNR)上的对比结果，如表 1所示。从平均PSNR可以看出本文算法总体上优于其他两种算法(仅对#4模糊核生成的模糊图像，文献[8]、文献[11]与本文算法的复原图像平均PSNR值相近)。

图 5给出文献[8]、文献[11]及本文算法复原图像在主观视觉上的一个对比，可以看出本文算法复原出的图像振铃较少。图中右上角分别是真实模糊核、对比算法及本文算法分别估计出的模糊核。

3.2 真实模糊实验

在真实模糊实验中，选择Köhler等^[25]提供的真实模糊图像基准数据库中的图像编号为I1~I4，模糊核编号为K1~K7的28张真实模糊图像。基准数据库中的原图如图 6所示。需要指出的是，如利用基准数据库中模糊图像的全图(800×800大小)来进行模糊图像盲复原，文献[25]中所测试的几种盲复原算法均能得到较好的复原效果。表 2中给出了文献[8]和^[11]与本文算法的PSNR对比结果，从中可以看出在800×800大小时，文献[8]复原图像的PSNR稍低，而文献[11]和本文算法复原图像的PSNR相近。

但前面提到，目前大部分的模糊图像盲复原算法均是利用图像中的显著边缘信息估计模糊核，如果图像中显著边缘过少或图像的尺寸太小，则模糊核的估计往往会失效。实验中我们裁剪出模糊图像中的部分区域(201×201大小)来进行实验，28幅裁剪后的模糊图像如图 7所示。

表 3中给出的是本文算法与对比算法复原图像在平均PSNR上的对比结果，可以看出本文算法在平均PSNR上与这两种算法相比有明显的改善。文献[8]和文献[11]在部分模糊核条件下其模糊图像复原结果的PSNR值低于模糊图像自身的PSNR值，其原因是对这些图像的模糊核估计失效。对于模糊核编号为K3所对应的模糊图像，因模糊比较轻微(模糊图像自身的PSNR值较高)，所以本文算法在PSNR值上改善并不太明显。

图 8(a)是图像编号为I3，模糊核编号为K5的模糊图像，图 8(b)~(d)分别是文献[8]、文献[11]及本文算法的复原结果。图中右上角所示是三种算法估计的模糊核，显然文献[8]和文献[11]估计到的模糊核是失效的，从而导致其复原的图像出现很强的振铃，复原效果较差。

4 结语

本文在张量框架下，提出一种基于张量总变分的模糊图像盲复原算法，该算法利用新张量框架中的t-product积，将彩色图像作为一个三阶张量，有效地利用彩色图像中三个通道的边缘信息来估计模糊核，最后再利用张量总变分复原出清晰图像。实验结果表明了本文算法的有效性，特别是在处理尺寸较小的模糊图像时，

本文算法复原出的图像优于其他的算法。但是，需要指出的是：因本文算法主要是解决传统模糊图像复原算法在处理显著边缘过少或图像尺寸较小图像时，模糊核估计失效的问题，因此本文并没有讨论模糊核尺寸较大的情况(如数据库中的模糊核编号为K8~K11的模糊图像)，而文献[11]的算法则能较好地估计出模糊核尺寸较大的模糊核，并能复原出较为满意的清晰图像。下一步的工作可以针对张量框架下的模糊核尺寸较大情形的模糊图像复原算法展开研究。