图像超分辨率(super-resolution)重建技术是针对采集到的同一场景中多幅较低分辨率(LR)图像，利用这些图像中所附带的信息结合重建算法估计出一幅或多幅高分辨率(HR)图像的过程. 该技术是对于图像退化或降质的一个逆过程，高分辨率图像和退化或降质后的低分辨率图像之间存在某种转化关系，从数学的角度考虑，这种转化关系可以由一个变换矩阵的线性系统和一个加性噪声来表示. 而变换矩阵的单一性体现了超分辨率重建的病态性，从而产生不适定问题. 对于这种存在不唯一解的病态问题，通常有效的办法就是引入正则化. 从实现域角度考虑，正则化方法主要分为频域正则化和空域正则化. 在频域正则化使用的算法中，Tsai等^[1]在1984年提出了一个重要的算法，该算法易于实现，但缺乏空域先验知识导致后续发展受阻. 而空域算法能够较为灵活地适应复杂的降质模型，如Henry Stark^[2]提出的凸集投影法(projection on to convex sets, POCS)，利用非线性退化模型，能较好地利用低分辨率图像的高频信息复原出一帧HR图像. 而最大后验概率估计(maximum a posteriori, MAP)^[3]和迭代反向投影法(iterative back projection, IBP)^[4]，受到配准精度、降质模型误差和信号噪声的影响，放大倍数受到限制，一般仅能达到1.6倍. 但是这类算法具有较好的信号还原能力，在图像配准精度较高的遥感图像领域取得了广泛应用和良好的重建效果. 此外，基于学习的超分辨率算法，通过样本库中的先验知识来约束重建过程，结合机器学习与统计学习的理论和方法，为重建提供更多的信息来源，从而提升了超分辨率的主客观质量. 近年来由于其良好的重建效果受到了广泛的关注，成为图像超分辨率研究的热点方向.

与此同时，正则项的模型对重建的效果起着至关重要的作用. Tikhonov正则化^[5]是一个经典的超分辨率正则化方法，该方法改善了噪声影响，但不能很好地保留重建图像的细节信息. 为了抑制去噪过程中丢失边缘信息的弊端，Rudin等^[6]提出了总变分正则化(total variation, TV)模型. 在该模型的基础上，Farsiu等^[7]进一步提出了Bilateral Total Variation (BTV)模型，此模型较TV模型而言有更好的鲁棒性，并且能够增强边缘信息.

Yang等^[8]将稀疏表示理论引入超分辨图像重建，使得高分辨率图像块之间的线性关系可以直接由低分辨率图像块之间的线性关系得到. 研究表明，若图像重建具备重建约束以及稀疏先验信息，即可使用稀疏表示重建HR图像. 然而，大量稀疏表示的方法中使用的字典是通过训练不相关的图像对，并且形成一个体量较大的字典对，因而字典训练时间较长. 为此，文献[9]对该方法进行了改进，提出训练过完备字典对可以极大地缩短训练时间. 在训练过程中，不兼容的高频细节训练集对重建结果容易造成影响，文献[10]中基于LR图像本身的自相似性超分辨率重建的方法有效地解决了这一问题.

事实上，大量的研究集中于通过凸优化获取较为稀疏的解. L₀正则化算法理论上可以得到精确的稀疏解，但是求解L₀是一个NP-hard问题，故难以求得结果. 而L₁正则化算法求得的解稀疏性并不理想. 文献[11]采用L_1/2范式代替L₁范式，并证明了L_1/2正则化算法较之后者具有更加稀疏的解集以及更优的鲁棒性. 在此基础上，本文提出了一种新的稀疏表示正则化超分辨重建方法，结合L_1/2非凸优化，并对迭代加权最小二乘法(Iterative Reweighted Least Squares，IRLS)进行转化，使其适用于非凸优化问题，建立一种噪声自适应正则化参数求解模型. 实验结果表明，本文提出的L_1/2稀疏正则化超分辨率重建方法比双三次插值和L₁稀疏表示方法能得到更小的均方根误差，和更好的去噪效果.

1 基于L_1/2正则化自适应图像超分辨率重建算法

研究超分辨率重建过程的第一步，需要建立一个合理的图像退化模型，以此观测LR图像序列与HR图像之间的变化. 广泛采用以下形式表示

${\mathit{\boldsymbol{y}}_k} = \mathit{\boldsymbol{DB}}{\mathit{\boldsymbol{F}}_k}\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{N}}_k},\quad k = 1,2, \cdots ,K.$

(1)

其中，y_k为第k幅LR图像，x为估计的HR图像，D、B、F_k分别对应第k幅图像的下采样算子，模糊化算子以及位移算子，N_k表示随机加性噪声.

针对退化模型对x的估计，通常使用正则化方法解决该逆过程转化所产生的不适定问题. 由此，模型可变化成式(2).

${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over {\mathit{\boldsymbol{x}}}} } \!=\! \arg \min \left\{ \sum\nolimits_{k \!=\! 1}^K {\rho (\mathit{\boldsymbol{DB}}{\mathit{\boldsymbol{F}}_k}\mathit{\boldsymbol{x}}\! -\! {\mathit{\boldsymbol{y}}_k}) \!\!+\!\! \lambda \psi (\mathit{\boldsymbol{x}})} \right\} ,k \!=\! 1,2,3, \cdots K.$

(2)

${\rho (\mathit{\boldsymbol{DB}}{\mathit{\boldsymbol{F}}_k}\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{y}}_k})}$ 为重建模型中描述待估计的HR图像与LR图像之间相似程度的数据保真项. ${\psi (\mathit{\boldsymbol{x}})}$ 则为超分辨率重建的正则项，该项中包含重建所需的先验信息并以此约束重建结果. λ项作为正则项权重系数，可用来调节正则项和保真项之间的权重关系，且过高的λ将会使图像边缘变得模糊，而过低将会影响去噪的表现.

1.1 图像的稀疏表示

大量研究表明，若图像重建具备重建约束以及稀疏先验信息，即可使用稀疏表示重建HR图像. 其中，图像超分辨率重建是一个病态问题，对于给定的LR输入图像可以有较多的HR图像满足其重建约束. 而稀疏先验信息可以进一步解决该问题的不适定性，只需假设每个HR图像块 ${\mathit{\boldsymbol{x}}_i}(i = 1,2, \cdots ,k)$ 可以由一个合适的过完备字典D_x线性表示，因此满足该条件. 相应地，可以产生以下模型

${\mathit{\boldsymbol{x}}_i} = \mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{\alpha }}_i},{\rm{s.t.}}\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\| \leqslant M,{\mathit{\boldsymbol{\alpha }}_i} \in {R^M}(i = 1,2,\cdots ,k).$

(3)

其中，α_i为HR图像块x_i对应的稀疏系数，D_x为过完备字典. 该式表明，根据过完备字典D_x和相应的稀疏系数α_i可以复原HR图像块x_i. 由文献[12]可知，稀疏系数α_i也可以通过输入的LR图像所获取的LR图像块y_k和LR字典D_y得到，因此可以用一个相同的稀疏系数α_i共同作用于LR图像块与其相应的HR图像块. 为了简化重建过程的求解，可以采用y_k最稀疏的表示如式(4)所示.

$\mathop {\min }\limits_{\{ {{\mathit{\boldsymbol{\alpha }}}_i}\} } {\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_0},{\rm{s.t.}}\left\| {\mathit{\boldsymbol{E}}\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{y}}}{\mathit{\boldsymbol{\alpha }}_i} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{y}}_k}} \right\|_2^2 \leqslant \varepsilon .$

(4)

式中E为特征提取操作矩阵， ${\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_0}$ 涉及L₀范式的求解，是一个NP-hard问题，理论上可以利用贪婪算法或门限法求得，但前者所需的时间代价过高，且无法保证收敛至全局最优；后者对噪声较为敏感，通常不能很好地收敛到全局最小. 从可行性角度考虑，可以使用L₁范式代替L₀范式，从而将问题演变为凸集优化问题. 然而由于重尾性，L₁范式也存在自身的弱点，并且稀疏性较差. 根据文献[13]提出的L_p (0<p<1)范式的框架，进一步延伸至非凸优化问题. 自然地，如何从0<p<1中选取一个最优的值使得稀疏性和可行性都达到较好水平成为焦点. 文献[11]通过研究提出了L_1/2范式，且验证了其稀疏性和鲁棒性较L₁范式更优，在可行性方面，对比L₀范式，它更易求解. 综合以上，本文选取L_1/2范式重建HR图像. 结合分析，问题模型(4)可转变成

$\mathop {\min }\limits_{\{ {\mathit{\boldsymbol{\alpha }}_i}\} } \left\| {\mathit{\boldsymbol{E}}\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{y}}}{\mathit{\boldsymbol{\alpha }}_i} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{y}}_k}} \right\|_2^2 + \lambda \left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_{1/2}^{1/2}.$

(5)

模型(5)建立在局部图像块的理论基础上，忽略了邻近块之间的重叠部分，故无法确保邻近块间的兼容性^[14]. 在此，添加一个控制邻近块相似性约束，致使由y_k重建所得的HR图像充分接近于前一步重建的邻近块. 因此，模型(5)可改良为以下优化问题

$\begin{array}{l}\mathop {\min }\limits_{\{ {\mathit{\boldsymbol{\alpha }}_i}\} } \left\| {\hat{D}{\mathit{\boldsymbol{\alpha }}_i} - {{\hat{y}}_k}} \right\|_2^2 + \lambda \left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_{1/2}^{1/2},\\[7pt]\hat{D} = \left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{E}}\mathit{\boldsymbol{D}}_{bm{y}}}\\[7pt]{{\mathit{\boldsymbol{T}}_k}\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{x}}}}\end{array}} \right],{{\hat{y}}_k} = \left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{y}}_k}}\\[7pt]{{\mathit{\boldsymbol{z}}_k}}\end{array}} \right].\end{array}$

(6)

其中，T_k为第k次作用的图像块与前一次重建所得的HR图像的重叠部分的矩阵，z_k为前一次重建的HR图像重叠部分的像素值.

1.2 字典训练模型

稀疏表示的超分辨率重建通常可以分为两个部分：字典学习和图像重建. 根据式(3)和(4)分析，在进行稀疏表示时，需要获取式中LR字典D_y和HR字典D_x两个矩阵. 而该步骤是字典学习阶段的关键. 一种直观的方法是采集保留高低分辨率对应关系的图像块对^[15]，但涉及的字典规模较大，易造成计算时效不佳的后果. 文献[16]中运用基于L_1/2正则化的字典训练模型，结合K-SVD算法得到一对高低分辨率字典对，并对样本中高低分辨率图像块进行稀疏表示. 实验证明，K-SVD算法随着字典的增大，图像整体重建效果有一定提升的趋势，但稀疏度的增强对字典整体误差没有明显的作用，因此，需要前期合理设置字典大小和稀疏度. 文献[9]中提出，训练一个更加完备的字典对可以加快计算速度. 通过采集的HR图像块 $X = {\rm{\{ }}{\mathit{\boldsymbol{x}}_1}{\rm{,}}{\mathit{\boldsymbol{x}}_2}{\rm{,}} \cdots {\rm{,}}{\mathit{\boldsymbol{x}}_n}{\rm{\} }}$ 和对应的LR图像块 $Y = {\rm{\{ }}{\mathit{\boldsymbol{y}}_1}{\rm{,}}{\mathit{\boldsymbol{y}}_2}{\rm{,}} \cdots {\rm{,}}{\mathit{\boldsymbol{y}}_n}{\rm{\} }}$ ，合并成一组高低分辨率图像块对 $P = {\rm{\{ }}X,Y{\rm{\} }}$ ，最终使得整体的HR图像块与LR图像块有相同的稀疏表示形式，即对应各自的字典有相同的稀疏表示系数. 针对两个字典D_x、D_y，可以利用文献[9]中拼接字典的方法获取. 结合文献[17]改进该方法，得出以下字典模型

$\mathop {\min }\limits_{{{\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{x}}}},{{\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{y}}}},\{ {{\mathit{\boldsymbol{\beta }}}_i}\} _{i = 1}^n} \sum\limits_{i = 1}^n {(\left\| {{{\mathit{\boldsymbol{x}}}_i} - {{\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{x}}}}{{\mathit{\boldsymbol{\beta }}}_i}} \right\|_2^2 + \left\| {{{\mathit{\boldsymbol{y}}}_i} - {{\mathit{\boldsymbol{D}}}_{{\mathit{\boldsymbol{y}}}}}{{\mathit{\boldsymbol{\beta }}}_i}} \right\|_2^2 + \gamma {{\left\| {{{\mathit{\boldsymbol{\beta }}}_i}} \right\|}_1}).} $

(7)

其中，γ为正则化系数，且γ>0. 通过该模型得到的字典表征了图像块的基本原型，如轮廓边缘等信息，相对于直接采集图像块训练字典，由于拼接字典模型具有完备的特性不仅时效显著，而且有更好的实现效果.

1.3 自适应正则化参数

通过指定正则化参数λ并结合已获取的高低分辨率字典，对采样的图像块对求解式(6)所表示的优化问题可以重建HR图像块. 在求得最优解前，需要确定合适的正则化参数λ的值. 正则化参数起着平衡保真项和稀疏项的作用. 一般地，正则化参数λ可由多次实验比较不同参数设定所得到的结果取最优值，以此固定λ. 然而，正则化参数的选取受输入数据噪声的影响，针对不同的输入低分辨率图像块，其各自的噪声水平也是不同的. 因此，自适应地调整正则化参数λ的数值是非常有必要的.

考虑到细节信息的贡献权重，文献[17]利用迭代加权最小二乘法(Iterative Reweighted Least Squares, IRLS)转化式(6)的最优化问题，通过迭代加权正则项，并结合约束条件表征自适应正则化参数. 由于IRLS方法通常基于p $\geqslant$ 1，对于本文中0<p<1的情况需要进行转换，因此，在选取自适应参数的过程中，本文用L₂范式替代L_1/2范式作用于正则项，使得IRLS方法有效作用于非凸优化问题，并继续对式(6)进行加权处理，转化成新的目标函数

$\mathop {\min }\limits_{\{ {\mathit{\boldsymbol{\alpha }}_i}\} } \left\| {\hat{D}{\mathit{\boldsymbol{\alpha }}_i} - {{\hat{y}}_k}} \right\|_2^2 + \lambda \sum\limits_{j = 1}^n {{w_i}\mathit{\boldsymbol{\alpha }}_{ij}^2} .$

(8)

其中 ${\mathit{\boldsymbol{\alpha }}_i} = ({\mathit{\boldsymbol{\alpha }}_{i1}},{\mathit{\boldsymbol{\alpha }}_{i2}}, \cdots ,{\mathit{\boldsymbol{\alpha }}_{in}})$ ，w_i为前一次的 $\mathit{\boldsymbol{\alpha }}_{ij}^{\left( {k - 1} \right)}$ 迭代所得的权值，根据IRLS方法, 推出如下具体迭代形式

$w_j^{(k)} = {\left[{(\mathit{\boldsymbol{\alpha }}_{ij}^{(k - 1)})^2} + \varepsilon \right]^{ - 3/4}}.$

(9)

同理，第k次迭代时，可计算

$\arg \mathop {\min }\limits_{\{ {\mathit{\boldsymbol{\alpha }}_i}\} } \left\| {\hat{D}{\mathit{\boldsymbol{\alpha }}_i} - {{\hat{y}}_k}} \right\|_2^2 + \lambda \sum\limits_{j = 1}^n {w_j^{(k)}\mathit{\boldsymbol{\alpha }}_{ij}^2} .$

(10)

在稀疏表示中，正则化参数对保证数据的稀疏性起着至关重要的作用. 因此，可用 $\lambda ({\mathit{\boldsymbol{\alpha }}_i})$ 代替固定的λ值进行自适应调整，同时建立一个新的加权模型 $S({\mathit{\boldsymbol{\alpha }}_i})$ .

$S({\mathit{\boldsymbol{\alpha }}_i}) = \arg \mathop {\min }\limits_{\{ {\mathit{\boldsymbol{\alpha }}_i}\} } \left\| {\hat{D}{\mathit{\boldsymbol{\alpha }}_i} - {{\hat{y}}_k}} \right\|_2^2 + \lambda ({\mathit{\boldsymbol{\alpha }}_i})\sum\limits_{j = 1}^n {w_j^{(k)}\mathit{\boldsymbol{\alpha }}_{ij}^2} .$

(11)

文献[17]提出，定义式中的正则化参数 $\lambda ({\mathit{\boldsymbol{\alpha }}_i})$ 与 $S({\mathit{\boldsymbol{\alpha }}_i})$ 呈线性相关，即 $\lambda ({\mathit{\boldsymbol{\alpha }}_i}) = \mu S\left( {{\mathit{\boldsymbol{\alpha }}_i}} \right)$ . 其中，μ作为约束系数控制凸度. 为了保证 $S({\mathit{\boldsymbol{\alpha }}_i})$ 对于任一α_i都是凸的，需要使μ满足

$\mu \leqslant (1/\sum\limits_{j = 1}^n {{w_i}\mathit{\boldsymbol{\alpha }}_{ij}^2} ).$

(12)

在得到 $\lambda ({\mathit{\boldsymbol{\alpha }}_i})$ 的关联项后，通过转化式(11)，可以生成 $\lambda ({\mathit{\boldsymbol{\alpha }}_i})$ 的表达式

$\lambda ({\mathit{\boldsymbol{\alpha }}_i}) = (\left\| {\hat{D}{\mathit{\boldsymbol{\alpha }}_i} - {{\hat{y}}_k}} \right\|_2^2)/(\frac{1}{\mu } - \sum\limits_{j = 1}^n {{w_i}\mathit{\boldsymbol{\alpha }}_{ij}^2} ).$

(13)

由此，可以求解得到自适应正则化参数 $\lambda ({\mathit{\boldsymbol{\alpha }}_i})$ .

1.4 算法实现过程

步骤一：输入LR测试图像Y，一组HR训练图像 ${{\rm{\{ }}{X_{{j}}}{\rm{\} }}_{{j}}}$ ，以及放大因子；

步骤二：从 ${{\rm{\{ }}{X_{{j}}}{\rm{\} }}_{{j}}}$ 提取训练块，并结合其对应的LR图像块建立高低分辨率图像块对 $P = {\rm{\{ }}X,Y{\rm{\} }}$ ；

步骤三：使用式(7)训练高低分辨率字典对D_x、D_y；

步骤四：用插值法对Y进行上采样，具体实现方法：对Y中每个LR图像块y_i用式(6)求解最优化问题，其中λ根据本文第四部分描述的方法，用式(13)替代，且通过式(3)重建第i幅HR图像块x_i；

步骤五：计算所有已重建的局部HR图像块x_i在重叠区域的平均值，合并成整体的HR图像；

步骤六：输出最终的HR图像.

2 实验结果与分析

本文实验通过使用本文方法、双三次插值法以及文献[9]中的方法进行图像超分辨率重建，获得一系列HR图像并对比实验效果. 一般而言，图像块的尺寸可以任意设定，但不同的尺寸会影响重建的效果以及优化时间：尺寸越小，重建结果越精细，且更耗时. 因此，重叠区域的尺寸需要设定一个合理的值. 根据本文提出的方法，实验首先需要设定字典训练的相关参数，其中字典大小设为1 000，图像块的尺寸为9，邻近块的重叠区域水平和垂直方向的大小合理设置成6个像素. 所有处理的图像裁去纹理区域并舍弃平滑部分. 由于人眼对视觉效果的区分主要靠自身对亮度的敏感性，故本文算法实现对于彩色图片只采用亮度通道^[18]. 而对于色度通道，使用双三次插值的方法进行上采样操作.

评价一个图像重建效果，不仅需要从肉眼观察上做出初步直观的判断，从重建图像本身的性能量化指标中更能客观反映结果图像的质量. 本文使用图像处理方法最为常用的评价指标：峰值信噪比(PSNR)和结构相似度(SSIM). 相应地，两个指标均在亮度通道计算得出. 实验选取的所有测试图片如图1所示. 为节约篇幅，本文仅使用Lena、Butterfly的实验结果作演示.

图2中a为重建图像的结果，图2中b中采集的局部细节信息，重建结果对应的PSNR和SSIM值由表1给出.

在与不同图像重建方法的横向比较后，需要对本文提出的采用自适应正则化参数方法效果进行纵向比较. 可通过设定不同的λ比较可得. 根据文献[9]的方法，对实验中的局部块设定一组λ=0.1作为对照，同时，设置另一组λ=0.25，实验结果如图3所示，评价指标值见表2.

以上实验结果表明，无论在直观的视觉角度还是计算数值上，本文提出的L_1/2稀疏正则化方法对于图像重建的效果明显优于其他对照的算法. 另外，不同的正则化参数较大程度影响重建的结果，基于自适应选取的正则化参数的方法较固定的参数选取具有更强的鲁棒性，同时实验也表明，一个合适的正则化参数可以重建出更加清晰的图像，其结果较为接近测试的原始图像.

3 结束语

本文提出了一种基于L_1/2稀疏表示的自适应正则化图像超分辨率重建算法. L_1/2非凸优化方法较L₁正则化可以获得更加稀疏的解，基于该方法的图像超分辨率重建能够提升生成的HR图像的质量. 除此之外，一个合适的正则化参数对于正则化方法至关重要，且不同的输入图像块对应的最佳正则化参数值也不同. 通过增加控制邻近块相似性的约束，建立具有邻近块相似性约束的L_1/2稀疏正则化模型，可以确保邻近块间的兼容性；采用加权L₂范式代替L_1/2范式，对迭代加权最小二乘法(IRLS)进行转化，使其可用于求解自适应参数. 建立的模型可以针对不同的低分辨率图像块自适应地调整正则化参数. 实验结果验证了本文方法的有效性. 改进字典训练方法，提高字典训练效率将是下一步研究的重点.