0 引言

由于随机噪声、光模糊、采样频率和运动变形等因素的影响，成像设备最后得到的常常是视觉效果较差的低分辨率图像。超分辨率重建(Super-Resolution Reconstruction, SRR)是指利用一张或者多张低分辨率图像重构出高像素密度、视觉效果好的高分辨率图像。SRR在无需改变硬件设备的前提下提高了数字图像的质量，成本低廉、效果好，在医疗、遥感、军事等方面有很大的应用价值。

自Tsai等^[1]首先提出频率域的SRR算法后，许多经典空域算法框架及其改进方法被提出，如：迭代反投影(Iterative Backward Projection, IBP)^[2]、凸集投影(Projection Onto Convex Sets, POCS)^[3]、最大后验(Maximum A Posterior, MAP)^[4]等。由于频率域算法(利用傅里叶变换转换到频率域处理)只能用于简单的平移运动，而空滤算法因能够处理非线性的运动而成为研究热点，其中基于MAP的算法框架使用较为广泛。在基于MAP的算法中，常利用正则项来构造MAP的求解路径。许多文献提出了有效的正则项形式。例如：基于马尔可夫随机场(Markov Random Field, MRF)正则^[5]、Tikhonov正则^[6]、总变分(Total Variation, TV)正则^[7]以及双边总变分(Bilateral Total Variation, BTV)^[8]等。

然而在上述这些具有梯度罚函数的正则化方法中，都存在着抑制图像噪声与保护图像细节信息的矛盾。为了提高重建图像质量，必须提供更加准确的先验模型来识别噪声和边缘以达到抑制噪声和保护边缘的平衡。许多文献中提出自适应正则项参数的方法^[9-11]。此类方法将正则项中一个参数修改为根据先验知识而自适应调整的参数。还有一些文献提出根据先验知识修改代价式中罚项的加权系数^[12-14]。这些方法都结合区域先验信息实现参数的自适应选取，其改善效果都取决于选择的先验模型的准确性。此类方法与一般正则化方法相比提高了去除噪声保边缘能力，改善了重建图像的质量；但同时增加了计算负担，且其先验模型的准确性仍可提高。

本文针对正则化超分辨率重建过程噪声消除的问题，在广义总变分(General TV, GTV)^[15]的基础上提出一种自适应阈值去噪的方法以改善正则化技术超分辨率重建在噪声环境下的重建效果。本文先叙述了MAP算法框架和GTV正则项，将其作为本文的基建算法；然后在此算法的迭代重建过程中加入本文提出的阈值分割模型，并分析了该阈值模型中两种不同先验知识受噪声影响的误差；最后通过实验比较验证了本文自适应阈值分割的方法，相比单一GTV正则化重建和改进的自适应参数方法，重建图像质量得到改善并加快了迭代收敛。

1 GTV正则技术重建 1.1 MAP求解框架

SRR的目标是从观测到的已知退化图像或者图像序列中重建估计出未知的清晰无噪的高分辨率(High Resolution, HR)图像。图像的降质过程表示如下：

$ {\mathit{\boldsymbol{y}}_k} = \mathit{\boldsymbol{Hx}} + {\mathit{\boldsymbol{n}}_k} $

(1)

其中：y_k表示第k帧r₁×r₂大小的带噪声的低分辨率(Low Resolution, LR)图像；x表示大小为r₁×N×r₂×N的原始HR图像，N为采样因子；因为异构的LR图像其退化模型相同，所以H可用来表示扭曲、模糊、下采样合成的降质矩阵；n_k表示第k帧LR图像的加性噪声。因为H为变态或不可逆矩阵，所以其反问题的求解变得十分困难，往往会因微小的扰动而产生大量噪声。为了得到一幅视觉上可被接受的HR图像, 常采用MAP算法框架增加正则项来迭代求解。其公式表示如下：

$ J(\mathit{\boldsymbol{x}}) = \frac{\beta }{2}\sum\limits_{k = 1}^K {\left\| {{\mathit{\boldsymbol{y}}_k} - \mathit{\boldsymbol{Hx}}} \right\|_2^2 + \lambda \ln (p(\mathit{\boldsymbol{x}}))} $

(2)

式中：K为LR图像序列帧数；β、λ为经验参数；ln(p(x))为正则项。最终求解HR过程便为最小化代价函数J(x)的过程。

$ \mathit{\boldsymbol{\hat x}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{x}} (J(\mathit{\boldsymbol{x}})) $

(3)

其中${\mathit{\boldsymbol{\hat x}}}$表示最终重建的结果。

1.2 GTV正则项

GTV正则是在TV的基础上进行推广的正则项，它较BTV更加准确地推广了TV正则项，其计算式如下：

$ \begin{array}{l} GTV({x_0}) = \sum\limits_{l = {\rm{1}}}^p {\sum\limits_{m = {\rm{1}}}^p {{\alpha ^{\sqrt {{l^2} + {m^2}} }}} } \sum\limits_i {\sqrt {{{(\Delta _i^{lh}{x_0})}^2} + {{(\Delta _i^{mv}{x_0})}^2}} } = \\ \sum\limits_{l = {\rm{1}}}^p {\sum\limits_{m = {\rm{1}}}^p {{\alpha ^{\sqrt {{l^2} + {m^2}} }}} } {\left\| {{\Delta _{lm}}{x_0}} \right\|_1} \end{array} $

(4)

式中：Δ_i^lh、Δ_i^mv分别表示第i个像素在水平和垂直方向上的l阶和m阶差分；x₀为像素点；p、α为经验参数。α值越小越能保护图像细节，但是不能抑制噪声；其值越大则越能抑制噪声却不能保证边缘细节。GTV正则项中用欧几里得范数描述了像素间空间距离的关系，用双边差分作为梯度罚函数。

综合式(2)、(4)，最后GTV正则重建的代价式为：

$ \begin{array}{l} J(\mathit{\boldsymbol{x}}) = \frac{\beta }{2}\sum\limits_{k = 1}^K {\left\| {{\mathit{\boldsymbol{y}}_k} - \mathit{\boldsymbol{Hx}}} \right\|_2^2 + } \\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{\lambda }}\sum\limits_{l = 1}^p {\sum\limits_{m = 1}^p {{\alpha ^{^{^{\sqrt {{l^2} + {m^2}} }}}}{{\left\| {{\Delta _{lm}}\mathit{\boldsymbol{x}}} \right\|}_1}} } \end{array} $

(5)

2 阈值选取

因为在选取参考帧进行初始化估计的时候，LR图像原本为带有噪声的图片，所以在GTV重建代价式J(x)中的保真项虽然可以抑制重建过程出现的偏差，但并没有感知噪声的能力。而在罚项中每个点在邻域内进行差分计算，若该点在邻域内存在噪声或者本身就是噪声点，则罚项会过大，所以罚项确有感知噪声的能力。同时若该点为边缘点，也会相对增加罚项的值，这便是抑制噪声和保护边缘细节的矛盾。由于这种不协调性，需要一种准确的先验知识甄别图像中的噪声和边缘，使得在重建过程中能够抑制噪声而保护边缘。本文提出一种区分噪声值和模型数据的阈值分割模型，直接截断大于阈值的噪声点，然后重新插值。这样只改变了部分噪声点，从而保持了其他细节。阈值的获取公式如下：

$ \begin{array}{l} T({x_0}) = \\ C \cdot \exp \left[{-\omega \frac{{{{\left| {({x_0}-{x_{{\rm{mid}}}})} \right|}^2}}}{{\frac{1}{{P \times Q}}\sum\limits_{i = i-P}^{i + P} {\sum\limits_{j = j - Q}^{j + Q} {\left| {({x_{i, j}} - {x_0})} \right| + 1} } }}} \right] \end{array} $

(6)

式中：P × Q表示局部窗口的大小；x_mid为该局部窗口内的中值；ω为经验参数；x₀为像素点。此先验模型通过某像素点与其感兴趣区域分为以下几种情况进行比较获取:

1) 若此点是边缘点或者此区域比较平滑，则其值都将与区域中值接近，从而得到的很小平方差作为式(6) 中的分子可以扩大阈值。

2) 若该点为噪声点，其值与中值的较大差值则可以缩小阈值。式(6) 中, 与中值点求差的因子体现了对边缘点、噪声点和平滑区域的识别能力。

3) 若此点不是边缘点，但处在非平滑区域也要稍微扩大阈值，以防止此点被当作噪声点处理。式(6) 中采用平均绝对离差作为识别非平滑区域的先验。因为分母确保为大于1的值，所以若平均绝对离差越大，则获取的阈值稍微增大，若其值越小，则分子为主要影响源。

在许多文献中标准差一般是区分平滑和非平滑区域的重要先验知识，但是其计算过程较复杂、误差较大且其受噪声影响而产生的误差不稳定，这将会导致阈值计算不稳定并且不准确。根据式(6) 的特点，只需知处理点与感兴趣点变化大小，而不必衡量其整体细节复杂程度。这使得平均绝对离差适用于此先验模型。本文采用平均绝对离差可以在稳定噪声干扰的同时保持此先验知识，可以对标准差和平均绝对离差受噪声的影响加以分析。

设局部窗口有k个点，在无噪声干扰的情况下表示为x′=[x₁, x₂, …, x_k]，然后将其中一点用噪声点n代替，则有噪声情况下为x″=[x₁, x₂, …, x_k-1, n]。于是得到无噪声和有噪声下窗口的均值分别为：

$ {\mu _o} = \sum\limits_{i = 1}^k {{x_k} = \frac{1}{k}\sum\limits_{i = 1}^{k - 1} {{x_i} + \frac{{{x_k}}}{k}} } = t + {\mathit{\boldsymbol{s}}_0} $

(7)

$ {\mu _n} = \frac{1}{k}\sum\limits_{i = 1}^{k - 1} {{x_i} + \frac{n}{k}} = t + {\mathit{\boldsymbol{s}}_n} $

(8)

其中：μ₀、μ_n分别为无噪声和有噪声下均值。在均值中产生的误差为μ_e=|μ₀-μ_n|=|s₀-s_n|。平均绝对离差在两种情况下可表示为：

$ {\upsilon _0} = \sum\limits_{i = 1}^k {({x_i} - {x_0})} $

(9)

$ {\upsilon _n} = \sum\limits_{i = 1}^{k - 1} {({x_i} - {x_0}) + (n - {x_0})} $

(10)

其中：υ₀、υ_n分别表示无噪声下和有噪声下的平均绝对离差。所以其误差为:

$ {w_e} = \left| {{\upsilon _0} - {\upsilon _n}} \right| = \left( {{x_k} - n} \right)/k = {\mu _e} $

(11)

同理，方差在无噪声和有噪声两种情况下可分别表示为：

$ \begin{array}{l} {\sigma _0}^2 = \frac{1}{k}\sum\limits_{i = 1}^k {{{({x_i} - {\mu _0})}^2}} = \frac{1}{k}\sum\limits_{i = 1}^k {{{({x_i} - t - {\mathit{\boldsymbol{s}}_0})}^2}} = \\ \frac{1}{k}\sum\limits_{i = 1}^{k - 1} {{{({x_i} - t)}^2}} + \frac{1}{k}\left[{{{({x_k}-t)}^2}-2{\mathit{\boldsymbol{s}}_0}\sum\limits_{i = 1}^k {({x_i}-t)} + k\mathit{\boldsymbol{s}}_0^2} \right] = \\ h + \frac{1}{k}\left[{{{({x_k}-t)}^2}-2{\mathit{\boldsymbol{s}}_0}\sum\limits_{i = 1}^k {({x_i}-t)} + k\mathit{\boldsymbol{s}}_0^2} \right] \end{array} $

(12)

$ \begin{array}{l} {\sigma _n}^2 = \frac{1}{k}(\sum\limits_{i = 1}^{k - 1} {{{({x_i} - {\mu _n})}^2} + {{(n - {\mu _n})}^2}} ) = \\ \frac{1}{k}(\sum\limits_{i = 1}^{k - 1} {{{({x_i} - t - {\mathit{\boldsymbol{s}}_n})}^2}} + {(n - t - {\mathit{\boldsymbol{s}}_n})^2}) = \\ \frac{1}{k}\sum\limits_{i = 1}^{k - 1} {{{({x_i} - t)}^2}} + \frac{1}{k}({({x_k} - {\mathit{\boldsymbol{s}}_n})^2} + {(n - t - {\mathit{\boldsymbol{s}}_n})^2}) = h + \\ \frac{1}{k}\left[{{{(n-t)}^2}-2{\mathit{\boldsymbol{s}}_n}\sum\limits_{i = 1}^{k-1} {({x_i} - t) - 2{\mathit{\boldsymbol{s}}_n}(n - t)} + k\mathit{\boldsymbol{s}}_n^2} \right] \end{array} $

(13)

其中：σ₀、σ_n分别表示在无噪声和有噪声下的标准差。最后方差受噪声的干扰产生的误差化简后为：

$ {\sigma _e}^2 = \left| {\sigma _0^2 - \sigma _n^2} \right| = {\mu _e}\left| {(1 - 1/k)({x_k} + n) - 2t/k} \right| $

(14)

由式(14) 可知，在局部窗口中误差的大小受到窗口大小和窗口均值影响，使得产生的误差不稳定，最终导致阈值不稳定。由于k一般为9、25或者更大的正整数，于是式(14) 中表达式可以近似为：

$ \begin{array}{l} {\sigma _e}^2 \approx {\mu _e}\left| {{x_k} + n - 2{\mu _0}/k} \right| = {\mu _e}\left| {k{\mu _e} + 2n - 2{\mu _0}/k} \right| = \\ {\mu _e}\left| {(k - 2/k){\mu _0} - k{\mu _n} + 2n} \right| \approx \\ {\mu _e}\left| {k({\mu _0} - {\mu _n}) + 2n} \right| = {\mu _e}\left| {k{\mu _e} + 2n} \right| > \mu _e^2 \end{array} $

(15)

此处假设μ₀ > μ_n以消除绝对值进行近似计算，取μ₀ < μ_n进行近似计算可以得到同样的结果。通过推导可知，当窗口k越大标准差最终的误差也越大。而反观平均绝对离差，其值在噪声干扰下产生的误差大小不随窗口大小和窗口均值变化而变化；其误差值小于标准差；其计算式简单。所以使用平均绝对离差将比使用标准差作为区分平滑区和非平滑区的先验知识更加准确且其阈值模型性能更加稳定。

最后按式(16) 迭代重建直到满足收敛条件。

$ {\mathit{\boldsymbol{X}}_{n + 1}} = \kappa \left\{ {\begin{array}{*{20}{c}} {{\rm{median}}(\mathop W\limits_{{S_{{\rm{local}}}}} ({x_0}))}\\ {{\mathit{\boldsymbol{X}}_n} - \lambda \Delta J({\mathit{\boldsymbol{X}}_n})} \end{array}\begin{array}{*{20}{c}} {, J({\mathit{\boldsymbol{X}}_n}) > T({\mathit{\boldsymbol{X}}_n})}\\ {, J({\mathit{\boldsymbol{X}}_n}) \le T({\mathit{\boldsymbol{X}}_n})} \end{array}} \right. $

(16)

式中：X表示迭代重建的图像矩阵；J表示对X中对应像素点计算后的代价矩阵；T表示对X中像素点计算后的阈值矩阵；W表示X中对应像素点x₀的感兴趣区域；median表示取区域中值；ΔJ(X_n)表示根据最速梯度下降法^[16]对式(5) 进行梯度计算。式(16) 使J中大于阈值的对应像素点被消除，然后取局部窗口中值对这些点进行重新插值，并且不再参与本轮迭代重建，同时使不大于阈值的对应像素点根据最速梯度下降方法迭代重建。算法执行流程如下：

步骤1  选取一帧低分辨率图像作为参考帧，对其余帧进行运动配准。

步骤2  对参考帧进行样条插值放大作为迭代的初始估计。

步骤3  计算下一次迭代的代价函数J(X)和自适应阈值T(X)。

步骤4  将J(X)与T(X)比较，对大于阈值的点进行重新插值，对小于阈值的点进行最速梯度下降迭代计算，最终得到新的一轮HR估计图像。

步骤5   判断此轮迭代是否收敛，如果不收敛则回到步骤3继续迭代，如果收敛则输出本轮HR估计图像。

3 实验结果及分析

为了验证算法的效果和性能，所有算法都在PC(Intel I3 3.69GHz, RAM 4.0GB), Matlab 2014a环境下进行。实验选取尺寸均为256×256的Lena、Planet、Books共3幅图像作为原始图像，对其进行平移、模糊和降采样后分别加入均值为0、方差为0.01的高斯噪声和密度为0.02的椒盐噪声。在迭代重建中，局部窗口大小均为5×5，迭代步长为0.1，在误差达到‖X_n+1-X_n‖²/‖X_n‖²≤10^-7时收敛。实验中所选取高斯噪声的密度大、噪声块小，而椒盐噪声密度小，噪声块大，使得实验在两种噪声环境中具有对比性。

本文方法将与样条插值、GTV重建和文献[13]中自适应参数的方法进行对比，并选取峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)^[17]和平均结构相似度(Mean Structural Similarity Index, MSSIM)^[18]两种指标对重建结果进行综合评价。PSNR是评价画质中使用最普遍的量测方法，它能衡量重建图像对于原图的失真程度，若PSNR越大则重建效果越好，失真程度小。MSSIM提取场景中的结构信息进行对比，其值的变化与感知图像质量变化非常接近，若MSSIM越大，则两幅图越相似。计算式如下：

$ PSNR = 10\lg \left( {{s^2}/MSE} \right) $

(17)

$ MSSIM = \frac{1}{M}\sum\limits_{i = 1}^N {SSIM(i)} $

(18)

$ SSI{M_{xy}} = \frac{{(2{\mu _x}{\mu _y} + {c_1})(2{\sigma _{xy}} + {c_2})}}{{(\mu _x^2 + \mu _y^2 + {c_1})(\sigma _x^2 + \sigma _y^2 + {c_2})}} $

(19)

式中：s为图片中最大像素点的值；MSE为对比图像的均方差；M为局部窗口数；μ_x、μ_y、σ_x、σ_y分别为两个比较区域的均值和方差；σ_xy为协方差。计算平均结构相似度时要将原始图和重建图划分为多个对应局部区域, 先计算局部区域的结构相似度，然后求其平均值。

图 1、图 2分别为在高斯白噪声下和椒盐噪声下的重建结果。

从图 1~2视觉效果看，相比单一GTV重建和样条插值算法，自适应参数GTV方法与本文方法在高斯噪声中的重建结果感官效果更好，而本文的方法相比自适应参数方法显得更加平滑。在椒盐噪声下的重建图像中，自适应参数的方法减少了噪声点且减小了噪声块，而本文方法不仅减少噪声点，而且去除了部分噪声块。相比高斯噪声下的结果，本文方法在噪声点密度较小的图像中有更好去噪效果。

表 1、表 2分别为高斯噪声和椒盐噪声下，各方法对三张图片进行重建后的评价数据。

从表 1~2中可以看出，在三张图片中各方法的PSNR和MSSIM从上到下依次递增，说明本文方法的重建图像与原始图像更为相似。本文方法相比自适应参数方法：在高斯噪声中PSNR提升1.0左右，MSSIM提升0.01左右；而在椒盐噪声中PSNR提升4.0左右，MSSIM提升0.03左右。这是由于在高斯噪声中，噪声点密度较大，使其在局部区域中相互干扰造成的，而在椒盐噪声中，其较小的噪声密度使得一微小局部区域中的噪声点被去除后不会再受到其他噪声点的干扰。从迭代次数可以看出，本文方法相比GTV正则和自适应参数的重建方法，较大幅度地减少了迭代次数，表明直接截断噪声点然后填补合理的值的阈值分割方式，确能使噪声点不会在迭代过程中反复干扰，从而加快收敛。实验结果表明：本文方法提高了正则化技术在噪声环境下的重建能力。

4 结语

本文针对正则化超分辨率重建在噪声环境中效果欠佳的问题，以GTV正则项为基础展开研究，提出一种自适应阈值的去噪方法。该方法在识别噪声时，用较准确的先验知识作为GTV代价式的阈值选取模型，使噪声点从有用模型数据中分割出来从而达到去噪目的；在处理噪声时，直接在迭代过程中截断识别出的噪声点并重新插入合适的值，使得噪声信息不会在GTV重建过程中反复干扰从而加快了收敛，同时保持了图像信息。对于非噪声点和重新插值后的点可以继续进行迭代重建，从而在保持GTV超分辨率重建效果的同时增强了其去噪性能。相比于单一的GTV正则化重建和自适应正则参数的方法，本文方法提高了图像重建效果，加快了迭代重建的收敛，使得正则化超分辨率技术在噪声环境下能够达到更好的效果。

虽然该方法在实验中取得了一定效果，但是阈值选取模型依旧会受到噪声点的干扰，抑制噪声与保持边缘细节的矛盾依旧存在，且对截断的噪声处插入中值会使重建结果平滑。在未来的研究中将寻找更准确的先验知识或者利用学习的方式来提升阈值模型的性能，使得此自适应阈值分割模型更加逼近噪声和模型数据的分割函数；同时利用更好的插值方式，取代中值插值法以达到去平滑效果。