图像在形成、传输和记录的过程中，不可避免会受到各种因素的影响，导致图像质量降低，典型体现在模糊、噪声、分辨率降低等方面。针对噪声的图像恢复技术被称为图像去噪，针对模糊的图像恢复技术被称为图像去模糊，针对分辨率降低的图像恢复技术被称为图像超分辨率重建。当降质条件中同时存在两种以上的降质因素时，图像恢复面临着更大的挑战。本文关注于模糊和分辨率降低同时存在时单幅降质图像的恢复问题，即单幅模糊图像超分辨率盲重建。
模糊降质因素出现的原因有很多，如相机和物体之间的相对运动、场景的深度不一、相机的聚焦变化等。模糊退化函数即模糊核，表示清晰场景图像在单点成像像素上累积的区域和权重，模糊核可分为全局一致模糊核和全局非一致模糊核，本文关注全局一致模糊核的降质情况。
单幅图像的超分辨率重建技术发展至今，已取得了显著的研究成果。它主要包括三个研究方向：基于插值的方法[1-2]、基于重建的方法[3-4]和基于学习的方法[5-7]。现有的绝大多数单幅图像超分辨率重建方法都是基于模糊核已知的假设，将图像先验的研究和设计作为研究重点。文献[8]首次进行了模糊核未知时超分辨率重建的性能分析，指出重建性能受图像先验模型选取和模糊核准确性两个关键因素影响，并经过理论分析和大量的实验证明，模糊核准确性对重建性能的影响最大。当估计的模糊核与真实模糊核之间存在较大误差时，超分辨率重建的性能急剧降低，这意味着在图像超分辨率盲重建时，需要重点关注模糊核的准确估计。
Levin等[9]指出，模糊核未知时，同时进行模糊核的估计和自然图像的恢复具有性能局限性，优化迭代中恢复的自然图像不利于模糊核的准确估计。因此，单幅图像超分辨率盲重建可分为两个独立的子过程：模糊核估计子过程和图像超分辨率非盲重建子过程，模糊核估计子过程是超分辨率盲重建的核心内容。在仅有的少量关于图像超分辨率盲重建的研究中，部分算法[10-12]基于参数化模型完成模糊核的估计，常用的模型为高斯函数，但高斯模型只适用于散焦模糊，对其他模糊情况下的重建问题失效；另一部分算法[13-14]则进行无参的模糊核估计子过程。Michaeli等[13]利用图像块跨尺度的冗余性，采用基于模糊核最大后验概率方法估计未知的模糊核，其搜索的最近邻块与参考块间的差异决定了估计模糊核的形态。Shao等[14]使用双l0-l2范数作为模糊核估计子过程的先验约束。由于lp范数在算法优化可靠性、稀疏重构信号效果等方面优于l0范数及l1范数[15-16]，因此本文使用双lp-l2范数约束模糊核估计子过程，并对图像梯度进行阈值判断获取自适应组合矩阵，将锐化图像lp-l2范数约束进行自适应组合，有效实现图像lp范数与l2范数约束功能的区域性分布。实验结果表明，本文方法能有效提高模糊核估计的准确性，进而提高低分辨率模糊图像超分辨率盲重建的质量。
1 模糊低分辨率图像降质模型模糊低分辨率图像中同时存在明显的模糊和低分辨率问题需要解决。模糊低分辨率图像对应的降质模型为：
| $\mathit{\boldsymbol{y}}=\left( \mathit{\boldsymbol{k}}*\mathit{\boldsymbol{x}} \right)\downarrow +\mathit{\boldsymbol{n}}$ | (1) |
其中：k表示模糊核；↓表示下采样操作；n为采集过程中的加性噪声；x表示原始的高分辨率图像；y是观测到的低分辨率模糊图像。k未知时，式(1) 的逆过程对应了超分辨率盲重建。本文重点研究超分辨率盲重建的核心内容，即模糊核估计子过程，图像超分辨率非盲重建子过程则采用基于非局部中心化稀疏表示的超分辨率方法[17]。模糊核估计子过程，实际也是式(1) 的逆求解过程。图像弱边缘及细微细节易受噪声和误差等因素影响，从而降低模糊核估计的准确性，因此在模糊核估计子过程中，需要为图像增加额外的非自然属性约束，平滑掉图像中的弱边缘和细微细节，只保留显著的强边缘。
基于最大后验概率方法，式(1) 对应的模糊核估计子过程可表示为：
| $\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right)} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}{J_\mathit{\boldsymbol{x}}}\left( \mathit{\boldsymbol{x}} \right) + {\alpha _\mathit{\boldsymbol{k}}}{J_\mathit{\boldsymbol{k}}}\left( \mathit{\boldsymbol{k}} \right)$ | (2) |
其中：K是模糊核k的矩阵表示；D是降低图像分辨率的下采样矩阵；Jx(x)表示对图像的约束先验，其中应包括图像非自然属性约束先验；Jk(k)为对模糊核的约束先验。模糊核估计子过程完成高分辨率锐化图像和模糊核的迭代交替优化，需要说明的是高分辨率锐化图像并非最终要重建的高分辨率图像。
2 基于自适应双lp-l2范数的超分辨率重建 2.1 双lp-l2范数约束的超分辨率重建模糊核估计子过程中增加的非自然属性先验包括：基于比值稀疏正则化的图像先验[18]，l0.3范数图像先验[19]，近似l0范数图像先验[20-21]，l0-l2范数图像先验[14, 22]等。Chartrand等[15]和Mazumder等[16]指出：基于lp范数优化算法在重构信号效果及可靠性方面优于l0及l1范数优化算法。因此，本文采用lp范数结合l2范数来约束锐化图像，使产生的高分辨率锐化图像边缘区域锐化，非边缘区域平滑，能够提高后续模糊核估计的准确性；同样地，模糊核也采用lp-l2范数约束，能消除模糊核估计过程中的孤立点和弱成分。本文提出的双lp-l2范数约束的超分辨率重建方法，其模糊核估计子过程可公式化表示为：
| $\begin{array}{l} \left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_p^p + \\ \quad \quad \quad \quad \quad \quad {\beta _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 \end{array}$ | (3) |
采用常规超分辨率方法(如A+[6])对低分辨率图像y进行超分辨率重建，能获取分辨率提高的图像
| $\begin{array}{l} \left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_p^p + {\beta _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \\ \quad \quad \quad \quad \quad \quad \quad {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \eta \left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 \end{array}$ | (4) |
正则化参数αx、βx、αk和βk的选取对模糊核估计准确性的影响较大。在式(4) 的迭代初期，需要加强双lp-l2先验的约束，从而能尽快完成图像的锐化、噪声和细节的平滑；而在迭代后期，则需要加强数据项的约束，使估计的模糊核与真实的模糊核逼近。因此我们采用了正则化参数迭代衰减的策略，使正则化参数按固定的衰减因子随迭代次数逐渐减小。
2.2 高分辨率锐化图像lp-l2范数约束的自适应组合对图像梯度进行l2范数约束能平滑图像，消除噪声，但不可避免会模糊边缘；lp(0＜p＜1) 范数约束能保持并锐化模糊图像中的边缘，有助于后续模糊核的估计。为了有效实现l2范数和lp范数功能的区域性分布，将式(4) 中高分辨率锐化图像lp-l2范数进行自适应组合，在平滑区域侧重l2范数的平滑和噪声抑制功能，在边缘区域侧重lp范数边缘锐化的功能。
为了完成lp-l2范数的自适应组合，按式(5) 构建自适应组合矩阵W：
| $\mathit{W}\left( {i,i} \right) = \left\{ {\begin{array}{*{20}{l}} {1,} & {\nabla \mathit{x}\left( i \right)thr}\\ {0,} & {\nabla \mathit{x}\left( i \right) ＜ thr} \end{array}} \right.$ | (5) |
其中：thr为梯度阈值，i为像素位置索引。按照矩阵W的指示，在图像梯度大于阈值thr的边缘区域使用lp范数约束，在图像梯度小于阈值thr的平滑区域使用l2范数约束，则式(4) 改写为：
| $\begin{array}{l} \left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}}} \right\|_p^p + \\ \quad {\beta _\mathit{\boldsymbol{x}}}\left\| {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \eta \left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 \end{array}$ | (6) |
其中I是单位矩阵。
3 算法流程 3.1 算法的交替优化流程本文采用交替迭代的优化方法完成式(6) 的求解。在第i次迭代时，首先固定模糊核ki-1，求解关于高分辨率锐化图像x的最小化问题：
| $\begin{array}{l} {\mathit{\boldsymbol{x}}_i} = {\rm{arg}}\mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{x}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}}} \right\|_p^p + \\ \quad \quad \quad \quad \quad {\beta _\mathit{\boldsymbol{x}}}\left\| {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \eta \left\| {{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 \end{array}$ | (7) |
然后固定高分辨率锐化图像xi，求解关于模糊核k的最小化问题：
| ${\mathit{\boldsymbol{k}}_i} = {\rm{arg}}\mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{k}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \eta \left\| {{\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2$ | (8) |
由于模糊核具有非负性和归一化约束，完成式(8) 的优化求解后，需要将k中负值进行归零，并将k归一化处理。
式(7) 和式(8) 都是非线性优化问题，本文采用罚函数技术与增广拉格朗日方法相结合的方式来完成它们的求解。式(7) 和式(8) 的具体实现细节将在3.2和3.3节中介绍。
3.2 自适应lp-l2范数约束的锐化图像优化方法通过罚函技术引入辅助变量u，式(7) 转换为等价的最小化问题：
| $\begin{array}{l} \quad \left( {{\mathit{x}_i},{\mathit{\boldsymbol{u}}_i}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{u}}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\mathit{\boldsymbol{Wu}}} \right\|_p^p\; + \\ \;\quad \quad \quad \quad \quad \quad {\beta _\mathit{\boldsymbol{x}}}\left\| {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \eta \left\| {{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2\\ {\rm{s}}.{\rm{t}}.\;\quad \mathit{\boldsymbol{Wu}} = \mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}} \end{array}$ | (9) |
其中将xi-1代入式(5) 求得W。式(9) 通过算子分裂完成了lp项和l2项的去耦合。增广拉格朗日方法通过求解式(10) 所示的无约束问题来获得式(9) 所示的等式约束问题的解。
| $\begin{array}{l} {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{u}}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\mathit{\boldsymbol{Wu}}} \right\|_p^p + {\beta _\mathit{\boldsymbol{x}}}\left\| {\left( {1 - \mathit{W}} \right)\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \\ \;\;\quad \quad \eta \left\| {{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 + \mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{x}}^{\rm{T}}\left( {\mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{Wu}}} \right) + \frac{{{\gamma _\mathit{\boldsymbol{x}}}}}{2}\left\| {\mathit{\boldsymbol{W}}\nabla \mathit{x} - \mathit{\boldsymbol{Wu}}} \right\|_2^2 \end{array}$ | (10) |
其中：μx为拉格朗日乘子，γx为罚参数。式(10) 的每次迭代过程可分为以下三个子步骤：
| $\begin{array}{l} \mathit{\boldsymbol{x}}_i^l = {\rm{arg}}\mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{x}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\beta _\mathit{\boldsymbol{x}}}\left\| {\left( {1 - \mathit{\boldsymbol{W}}} \right)\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \eta \left\| {{\mathit{\boldsymbol{K}}_{i - 1}}\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 + \\ \;\quad \quad \quad \quad \quad {\left( {\mathit{\boldsymbol{\mu }}_x^{l - 1}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{Wu}}_i^{l - 1}} \right) + \frac{{{\gamma _\mathit{x}}}}{2}\left\| {\mathit{\boldsymbol{W}}\nabla \mathit{\boldsymbol{x}} - \mathit{\boldsymbol{Wu}}_i^{l - 1}} \right\|_2^2 \end{array}$ | (11) |
| $\mathit{\boldsymbol{u}}_i^l = {\rm{arg}}\mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{u}} {\alpha _\mathit{\boldsymbol{x}}}\left\| \mathit{\boldsymbol{u}} \right\|_p^p + {\left( {\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{x}}^{l - 1}} \right)^{\rm{T}}}\left( {\nabla \mathit{\boldsymbol{x}}_i^l - \mathit{\boldsymbol{u}}} \right) + \frac{{{\gamma _x}}}{2}\left\| {\nabla \mathit{\boldsymbol{x}}_i^l - \mathit{\boldsymbol{u}}} \right\|_2^2$ | (12) |
| $\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{x}}^l = \mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{x}}^{l - 1} + {\gamma _\mathit{\boldsymbol{x}}}\left( {\mathit{\boldsymbol{x}}_i^l - \mathit{\boldsymbol{u}}_i^l} \right)$ | (13) |
式(11) 是简单的凸二次规划问题，具有闭合解；式(12) 为lp范数非凸稀疏编码问题，解决方法包括迭代重加权最小二乘法(Iteratively Reweighted Least Squares, IRLS)[15]、迭代重加权l1最小化(Iteratively Reweighted l1-minimization, IRL1) [23]、迭代阈值方法(Iteratively Thresholding Method-lp, ITM-lp)[24]、查找表[25]等。由于ITM-lp方法易于实现且收敛性能优异，本文采用ITM-lp方法求解式(12)。
3.3 lp-l2范数约束的模糊核优化方法同样地，通过罚函数法引入辅助变量v，并结合增广拉格朗日方法，式(8) 转变为：
| $\begin{array}{l} \left( {{\mathit{\boldsymbol{k}}_i},{\mathit{\boldsymbol{v}}_i}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{k}},\mathit{\boldsymbol{v}}} \lambda \left\| {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{v}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \\ \quad \quad \quad \quad \;\quad \quad \eta \left\| {{\mathit{\boldsymbol{X}}_i}\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 + \mathit{\mu }_\mathit{\boldsymbol{k}}^{\rm{T}}\left( {\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{v}}} \right) + \frac{{{\gamma _\mathit{\boldsymbol{k}}}}}{2}\left\| {\mathit{\boldsymbol{k}} - \mathit{\boldsymbol{v}}} \right\|_2^2 \end{array}$ | (14) |
式(14) 的每次迭代过程也包括与式(11) ~(13) 相似的三个子步骤，不再赘述。需要说明的是，正如文献[20]所述，图像梯度域的使用能够避免亮度变化对模糊核估计的影响，所以式(14) 中Xi、y和
为了验证算法的有效性，进行了模拟模糊图像超分辨率重建实验和真实模糊图像超分辨率重建实验。模拟实验中，从Berkeley Segmentation Database BSDS500中随机选取了8幅高分辨率测试图像，采用文献[9]提供的8个运动模糊核(记为Kenel 1~8)，分别将高分辨率测试图像模糊再2倍下采样，得到64幅不同模糊程度的低分辨率图像。选用的运动模糊核如图 1所示。
|
图 1 Levin图像库[9]提供的运动模糊核 Figure 1 Motion blur-kernels from the benchmark image dataset proposed by Levin[9] |
实验中本文算法相关参数的经验值设置如下：式(5) 中梯度阈值为0.05；式(4) 中参数λ为0.01，参数η为100，αx为4，βx为100，αk为10，βk为1，αx和βx的衰减因子为2/3，αk和βk的衰减因子为4/5；式(7) 中罚参数γx为100；式(11) 中罚参数γk为106；式(4) 对应的外部循环次数为10；式(7) 和式(11) 对应的内部循环次数都设为10。客观评价指标选取峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)、结构相似度(Structural Similarity Index Measure, SSIM)和差方和(Sum of Squared Difference, SSD)比[9]。PSNR和SSIM是常用的图像质量评价指标，它们的值越大表明重建图像质量越高；差方和比是估计模糊核重建得到的图像跟原始高分辨率图像之间的差方和除以真实模糊核重建得到的图像跟原始高分辨率图像之间的差方和，公式表示为:
| ${\left\| {{\mathit{\boldsymbol{x}}_\mathit{\boldsymbol{k}}} - \mathit{\boldsymbol{x}}} \right\|^2}/{\left\| {{\mathit{\boldsymbol{x}}_{{\rm{kgt}}}} - \mathit{\boldsymbol{x}}} \right\|^2}$ | (15) |
其中：xk表示由估计模糊核重建得到的高分辨率图像；xkgt表示由真实模糊核重建得到的高分辨率图像；x表示原始高分辨率图像。差方和比根据重建图像质量来度量模糊核估计的准确性，并排除了超分辨率重建方法对模糊核准确性评估的影响，差方和比值越小表明模糊核估计越准确，重建图像质量越高，理想情况下值为1。
模糊核估计子过程中，本文采用了自适应lp-l2范数来约束锐化图像的重建。为了说明这种lp-l2范数组合形式的有效性，将本文方法与使用单一范数的重建方法进行比较，式(16) 表示使用l2范数的方法，式(17) 表示使用lp范数的方法：
| $\begin{array}{l} \left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\beta _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_2^2 + \\ \quad \quad \quad {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \eta \left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 \end{array}$ | (16) |
| $\begin{array}{l} \left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \right) = {\rm{arg}}\mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{x}},\mathit{\boldsymbol{k}}} \lambda \left\| {\mathit{\boldsymbol{DKx}} - \mathit{\boldsymbol{y}}} \right\|_2^2 + {\alpha _\mathit{\boldsymbol{x}}}\left\| {\nabla \mathit{\boldsymbol{x}}} \right\|_p^p + \\ \quad \quad \quad {\alpha _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_p^p + {\beta _\mathit{\boldsymbol{k}}}\left\| \mathit{\boldsymbol{k}} \right\|_2^2 + \eta \left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{\hat x}}} \right\|_2^2 \end{array}$ | (17) |
图 2给出了本文方法与这两个方法的性能比较结果，从中可以看出，本文算法的各项客观评价指标都优于使用单一范数的方法，从而验证了自适应lp-l2范数组合形式约束图像重建的有效性。
|
图 2 本文方法与使用单一范数的重建方法的实验结果比较 Figure 2 Comparison results of super-resolution reconstruction methods with single norm and the proposed method |
进一步地，将本文方法与级联方法1(先去模糊[18]再超分辨率重建[5])、级联方法2(先去模糊[18]再超分辨率重建[7])以及文献[14]方法进行比较。表 1分别给出了每一种运动模糊核情况下，8幅对应低分辨率模糊图像超分辨率盲重建的平均PSNR值、平均SSIM值和平均差方和比值。由表 1可以看出，除了Kernel 7下本文方法客观指标略低于文献[14]方法的客观指标，其他情况下本文方法都呈现较优的客观评价性能。
|
表 1 超分辨率盲重建的平均峰值信噪比、平均结构相似度和平均差方和比值比较 Table 1 Average PSNR, SSIM and SSD ratio comparison of blind super-resolution reconstruction results |
图 3给出了64幅低分辨率模糊图像用不同方法重建时的差方和比累积直方图。图 3的横坐标表示差方和比参考值，纵坐标表示64幅图中重建差方和比小于当前参考值的图像数量百分比。每个差方和比参考值对应的图像数百分比越大，表明按当前差方和比参考值为评价标准，成功估计的模糊核越多，因此算法的性能越好。图 3中，本文算法在每个差方和比参考值上的图像数百分比都大于比较算法，由此可见，本文方法的性能是优于比较算法的。文献[9]指出，当差方和比小于等于3时，盲重建的图像是直观可信的，盲重建方法是成功的，级联方法1差方和比小于等于3的成功率为23%，级联方法2差方和比小于等于3的成功率为33%，文献方法差方和比小于等于3的成功率为90%，而本文方法差方和比小于等于3的成功率达到了92%。
|
图 3 差方和比的累积直方图 Figure 3 Cumulative histogram of SSD ratios |
图 4是模糊核为Kernel 6时图像Woman的盲重建结果。为了更直观地进行视觉比较，图 4在每幅图的左下角给出了局部边缘图的放大显示。与文献[14]方法相比，本文方法估计的模糊核更稳定、噪点更少。同时，模糊核估计的不完全准确性会导致边缘区域的振铃现象，而从四种方法的重建结果可以看出，本文方法的重建图像最接近原始高分辨率图像，振铃现象也比另外三种对比方法弱，从而验证了本文方法对模糊图像超分辨率盲重建的有效性。
|
图 4 模糊核Kernel 6时测试图像Woman的重建结果比较 Figure 4 Comparison of super-resolution reconstruction results on image Woman blurred using Kernel 6 |
图 5给出了一幅真实模糊图像的超分辨率盲重建结果。通过比较四种方法的重建结果图和放大的局部边缘图，可以得到跟模拟实验相同的结论：本文方法重建结果图边缘更清晰，振铃现象有效减少，从而验证了本文方法对模糊核的估计更准确。
|
图 5 真实模糊图像的超分辨率盲重建结果比较 Figure 5 Comparison of super-resolution reconstruction results on a real blurred image |
本文重点研究了超分辨率盲重建的核心内容，即模糊核的估计，非盲重建子过程采用了基于非局部中心化稀疏表示的超分辨率方法[17]。为了避免所用的非盲超分辨率重建方法对时间复杂度分析的影响，本文只进行了算法模糊核估计的运行时间比较。实验中比较了四种算法对模拟模糊实验中64幅测试图的平均模糊核估计时间，所有算法使用了相同的实验平台(Intel 4.0 GHz CPU, 16 GB内存和Matlab R2012a)。级联方法1和级联方法2都使用了文献[18]的去模糊方法，因此两者具有相同的模糊核估计时间，都为5 s；文献[14]方法的平均模糊核估计时间为139 s，本文方法的平均模糊核估计时间为182 s。与文献[14]相比，本文方法的模糊核估计时间有少量增加；但是，从图 6所示PSNR性能与模糊核估计的综合比较可以看出，本文方法以少量的时间损耗为代价，换取了性能的提升。
|
图 6 重建PSNR与模糊核估计时间的综合比较 Figure 6 PSNR versus running time of blur-kernel estimate |
为了改善低分辨率模糊图像的质量，本文提出了一种基于自适应双lp-l2范数的超分辨率盲重建方法。超分辨率盲重建分为两个独立的子过程：模糊核估计和超分辨率非盲重建。模糊核估计是高分辨率锐化图像重建和模糊核估计的交替优化过程，在此过程中，使用双lp-l2范数同时约束锐化图像和模糊核。使用lp-l2范数约束锐化图像，能滤除易受噪声及误差影响的弱边缘和细节，只保留重要的强边缘参与模糊核的估计；使用lp-l2范数约束模糊核，能消除核的孤立点，保持核的连续性。同时，使用图像梯度的阈值分割，实现锐化图像lp-l2范数约束的自适应组合，加强l2范数在非边缘区域的平滑能力和lp范数在边缘区域的锐化能力。结合估计到的模糊核，使用非局部中心化稀疏表示的非盲超分辨率方法完成高分辨率图像的重建。实验结果表明，本文方法能准确地估计出模糊核，重建的图像具有较优的质量。本文研究了针对全局一致模糊核的模糊图像重建问题，而具有全局非一致模糊核的模糊图像在实际生活中也是普遍存在的，针对这类图像的重建问题有待进一步研究。
| [1] | ZHANG X, WU X. Image interpolation by adaptive 2-D autoregressive modeling and soft-decision estimation[J]. IEEE Transactions on Image Processing, 2008, 17(6): 887-896. DOI:10.1109/TIP.2008.924279 |
| [2] | 季成涛, 何小海, 符耀庆, 等. 一种基于正则化的边缘定向插值算法[J]. 电子与信息学报, 2014, 36(2): 293-297. (JI C T, HE X H, FU Y Q, et al. An edge directed interpolation algorithm based on regularization[J]. Journal of Electronics & Information Technology, 2014, 36(2): 293-297.) |
| [3] | YAN Q, XU Y, YANG X, et al. Single image superresolution based on gradient profile sharpness[J]. IEEE Transactions on Image Processing, 2015, 24(10): 3187-3202. DOI:10.1109/TIP.2015.2414877 |
| [4] | REN C, HE X, TENG Q, et al. Single image super-resolution using local geometric duality and non-local similarity[J]. IEEE Transactions on Image Processing, 2016, 25(5): 2168-2183. DOI:10.1109/TIP.2016.2542442 |
| [5] | ZEYDE R, ELAD M, PROTTER M. On single image scale-up using sparse-representations[C]//Proceedings of the 2010 International Conference on Curves and Surfaces, LNCS 6920. Berlin:Springer-Verlag, 2010:711-730. |
| [6] | TIMOFTE R, DE SMET V, VAN GOOL L. A+:Adjusted anchored neighborhood regression for fast super-resolution[C]//ACCV 2014:Proceedings of the 2014 Asian Conference on Computer Vision, LNCS 9006. Berlin:Springer-Verlag, 2014:111-126. |
| [7] | LI T, HE X, TENG Q, et al. Rotation expanded dictionary-based single image super-resolution[J]. Neurocomputing, 2016, 216: 1-17. DOI:10.1016/j.neucom.2016.06.066 |
| [8] | EFRAT N, GLASNER D, APARTSIN A, et al. Accurate blur models vs. image priors in single image super-resolution[C]//ICCV 2013:Proceedings of the 2013 IEEE International Conference on Computer Vision. Washington, DC:IEEE Computer Society, 2013:2832-2839. |
| [9] | LEVIN A, WEISS Y, DURAND F, et al. Understanding and evaluating blind deconvolution algorithms[C]//CVPR 2009:Proceedings of the 2009 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC:IEEE Computer Society, 2009:1964-1971. |
| [10] | BÉGIN I, FERRIE F P. PSF recovery from examples for blind super-resolution[C]//ICIP 2007:Proceedings of the 2007 IEEE International Conference on Image Processing. Piscataway, NJ:IEEE, 2007:421-424. |
| [11] | WANG Q, TANG X, SHUM H. Patch based blind image super resolution[C]//ICCV 2005:Proceedings of the 2005 IEEE International Conference on Computer Vision. Piscataway, NJ:IEEE, 2005:709-716. |
| [12] | HE Y, YAP K-H, CHEN L, et al. A soft MAP framework for blind super-resolution image reconstruction[J]. Image and Vision Computing, 2009, 27(4): 364-373. DOI:10.1016/j.imavis.2008.05.010 |
| [13] | MICHAELI T, IRANI M. Nonparametric blind super-resolution[C]//ICCV 2013:Proceedings of the 2013 IEEE International Conference on Computer Vision. Washington, DC:IEEE Computer Society, 2013:945-952. |
| [14] | SHAO W-Z, ELAD M. Simple, accurate, and robust nonparametric blind super-resolution[C]//ICIG 2015:Proceedings of the 8th International Conference on Image and Graphics, LNCS 921. Berlin:Springer, 2015:333-348. |
| [15] | CHARTRAND R, YIN W. Iteratively reweighted algorithms for compressive sensing[C]//ICASSP 2008:Proceedings of the 2008 IEEE International Conference on Acoustics, Speech and Signal Processing. Piscataway, NJ:IEEE, 2008:3869-3872. |
| [16] | MAZUMDER R, FRIEDMAN J H, HASTIE T. SparseNet:Coordinate descent with nonconvex penalties[J]. Journal of the American Statistical Association, 2011, 106(495): 1125-1138. DOI:10.1198/jasa.2011.tm09738 |
| [17] | DONG W, ZHANG L, SHI G, et al. Nonlocally centralized sparse representation for image restoration[J]. IEEE Transactions on Image Processing, 2013, 22(4): 1620-1630. DOI:10.1109/TIP.2012.2235847 |
| [18] | KRISHNAN D, TAY T, FERGUS R. Blind deconvolution using a normalized sparsity measure[C]//CVPR 2011:Proceedings of the 2011 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC:IEEE Computer Society, 2011:233-240. |
| [19] | KOTERA J, ŠROUBEK F, MILANFAR P. Blind deconvolution using alternating maximum a posteriori estimation with heavy-tailed priors[C]//CAIP 2013:Proceedings of the 2013 International Conference on Computer Analysis of Images and Patterns, LNCS 8048. Berlin:Springer-Verlag, 2013:59-66. |
| [20] | XU L, ZHENG S, JIA J. Unnatural L0 sparse representation for natural image deblurring[C]//CVPR 2013:Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC:IEEE Computer Society, 2013:1107-1114. |
| [21] | 闫敬文, 彭鸿, 刘蕾, 等. 基于L_0正则化模糊核估计的遥感图像复原[J]. 光学精密工程, 2014, 22(9): 2572-2579. (YAN J W, PENG H, LIU L, et al. Remote sensing image restoration based on zero-norm regularized kernel estimation[J]. Optics and Precision Engineering, 2014, 22(9): 2572-2579.) |
| [22] | SHAO W-Z, LI H-B, ELAD M. Bi-l0-l2-norm regularization for blind motion deblurring[J]. Journal of Visual Communication and Image Representation, 2015, 33: 42-59. DOI:10.1016/j.jvcir.2015.08.017 |
| [23] | CANDÈS E J, WAKIN M B, BOYD S P. Enhancing sparsity by reweighted l1 minimization[J]. Journal of Fourier Analysis and Applications, 2008, 14(5/6): 877-905. |
| [24] | SHE Y. An iterative algorithm for fitting nonconvex penalized generalized linear models with grouped predictors[J]. Computational Statistics & Data Analysis, 2012, 56(10): 2976-2990. |
| [25] | KRISHNAN D, FERGUS R. Fast image deconvolution using hyper-Laplacian priors[C]//NIPS 2009:Proceedings of the 200923rd Advances in Neural Information Processing Systems. Cambridge, MA:MIT Press, 2009:1033-1041. |