CN102360498A - 图像超分辨率重建方法 - Google Patents

Abstract

本发明公开了一种图像超分辨率重建方法，包括如下步骤：步骤(1)输入用向量表示的m帧针对同一场景的低分辨率图像序列；步骤(2)采用双线性插值法对所有低分辨率图像进行处理，获取初始高分辨率图像估计，即初值估计

步骤(3)将初值估计

代入迭代公式进行迭代计算；步骤(4)判断是否满足迭代收敛条件，若满足则停止迭代，输出所得高分辨率图像。本发明提供的图像超分辨率重建方法在视频和图像处理、医学成像、遥感图像等领域中有着广泛的应用前景。

Description

图像超分辨率重建方法

技术领域

本发明属于图像信号处理技术领域，具体涉及一种洛伦兹函数结合自适应正则化的图像超分辨率重建方法。

背景技术

超分辨率图像重建技术是运用信号与信息处理方法来提高分辨率的一种有效手段。如果有多帧关于同一场景的欠采样图像，这些图像之间可能有空间平移、偏转、缩放等几何差异，通过这些图像信息之间的互补、融合及图像先验信息，同时去除噪声以及模糊，则可以得到高分辨率图像。此方法克服了物理条件的限制，克服了图像插值等方法的局限性，利用多帧图像信息弥补空间分辨率的不足；改善由于图像降晰、离散化引起的空间分辨率下降；发掘现有图像数据(如多角度、多时相、多平台遥感图像，序列图像等)的潜力；突破了图像获取手段的空间分辨率极限；效果显著，应用方便，代价低廉。

目前，超分辨率图像重建技术研究中，基于正则化的重建方法和基于统计的重建方法是较为常用、效果也较好的主流方法。

基于正则化的重建方法是结合图像和模糊的先验知识对问题进行正则化处理，根据退化图像序列的前向模型，构造正则化最小化泛函，求解最小化泛函以得到高分辨率图像。该方法的优点是可直接加入先验约束、确保解的存在和唯一、收敛稳定性高、降噪能力强。但运算量大、收敛较慢，高分辨率图像的细节容易被平滑掉。

基于统计的重建方法是把超分辨率重建问题看成统计估计问题，常见有最大后验概率估计(MAP，maximun a posteriori)和最大似然估计(ML，maximum likelihoodestimation)两种算法。MAP算法是在已知低分辨率图像序列的条件下，使出现高分辨率图像的后验概率达到最大。实际上，MAP算法就是将目标函数最小化，并在目标函数中引入不同的先验模型决定其收敛性。ML算法则可以认为是MAP算法在等概率先验模型下的特例。

基于正则化的重建方法和基于MAP框架的超分辨率重建本质上都是通过引入若干先验知识，将不适定性重建问题转化为适定性重建问题，同时实现数据保真和边缘保持的高分辨率图像重建。由于数字图像处理领域中关于图像的先验知识具体对应的为图像模型，因此对基于正则化理论的超分辨率重建算法和基于最大后验概率框架的超分辨率重建算法的先验知识研究具有共通之处。

发明内容

发明目的：本发明的目的在于提供一种新的自适应正则化图像超分辨率重建算法，即洛伦兹函数结合自适应正则化的图像超分辨率重建算法。

为了实现本发明方法，本发明公开了一种图像超分辨率重建方法，其特征在于，包括如下步骤：

步骤(1)输入用向量表示的m帧针对同一场景的低分辨率图像序列；

步骤(2)采用双线性插值法对所有低分辨率图像进行处理，获取初始高分辨率图像估计，即初值估计

步骤(3)将初值估计

代入迭代公式进行迭代计算，其中，β为设定控制步长，设定值范围为0～1，n表示迭代次数，

表示经过第n次迭代的高分辨率图像估计，G_n表示经过第n次迭代计算得到的梯度；当n＝0时，公式变为

则可将带入。

根据下式在每次迭代中计算梯度G_n，

$G_n=[\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{H}_k^T\·\mathrm{\η}\·{\mathrm{\Ψ}}_{\mathrm{lor}}^{\′}(H_k{\hat{Z}}_n-Y_k)]\·\frac{1}{\mathrm{\α}}+{\mathrm{\λ}}_k\·C^T\·{\mathrm{\Ψ}}_{\mathrm{lor}}^{\′}\left(C{\hat{Z}}_n\right),$

其中，

参数T₁是洛伦兹函数的半值半宽，参数T₁为控制阈值，设定值范围为1～12之间的实数；C表示拉普拉斯算子，C＝[1/8 1/8 1/8；1/8 -1 1/8；1/81/8 1/8]；α和η为常数系数，α设定值范围为0～100之间的实数，η设定值范围为100～500之间的实数；H_k为等效降质因子；λ_k为自适应正则化参数。

根据下式在每次迭代中计算自适应正则化参数λ_k：

${\mathrm{\λ}}_k={\left|\right|{\mathrm{\Ψ}}_{\mathrm{lor}}(Y_k-H_k{\hat{Z}}_n)\left|\right|}^2/(2{\left|\right|Y_k\left|\right|}^2-{\left|\right|{\mathrm{\Ψ}}_{\mathrm{lor}}\left(C{\hat{Z}}_n\right)\left|\right|}^2);$

其中，

参数T2是洛伦兹函数的半值半宽，为控制阈值设定值范围为1～12之间的实数；

步骤(4)判断是否满足迭代收敛条件，若满足则停止迭代，输出所得高分辨率图像；收敛条件设定为 $\mathrm{\&rho;}={\left|\right|{\hat{Z}}_{n+1}-{\hat{Z}}_n\left|\right|}^2/{\left|\right|{\hat{Z}}_n\left|\right|}^2<5\&times;{10}^{-7}.$

本发明步骤(1)中输入的低分辨率图像序列是满足如下条件建立高分辨率图像的降质模型：针对带限连续的景物图像，成像系统以奈奎斯特频率采样，得到理想高分辨率图像Z；在成像过程中，只存在运动、降质因素、噪声影响；高分辨率图像Z依次经过变形、模糊、下采样、添加噪声后得到低分辨率观测图像序列Y_k；

图像降质模型用下式表示：

对单帧高分辨率图像进行降质处理，依次进行运动平移、模糊、降采样、添加噪声；重复该步骤，获取m帧关于同一景物的不同的低分辨率图像，设H_k＝D_kB_kW_k，称为等效降质因子；

所述Z是理想高分辨率图像，大小为qN₁×qN₂像素，Z由高分辨率图像像素以向量形式按字典顺序(是已知的一种向量的排列方法)排列而成，即

N＝qN₁×qN₂，为向量Z中的第N个分量，N表示向量Y_k中的分量个数，q为采样比率；

Y_k是低分辨率像序列，每一帧图像的大小都是N₁×N₂像素，Y_k由低分辨率图像像素以向量形式按字典顺序排列而成，即Y_k＝[Y_k，l，Y_k，2，……Y_k，M]^T，k＝1，2，…，m，M＝N₁×N₂，M表示向量Y_k中的分量个数；

D_k表示下采样矩阵，矩阵大小为(N₁N₂)²×q²N₁N₂；

B_k表示模糊矩阵，矩阵大小为q²N₁N₂×q²N₁N₂；

W_k表示获取图像时的全局或局部几何变形，旋转，矩阵大小是q²N₁N₂×q²N₁N₂；

V_k表示成像系统的噪声影响，设为均值为零，方差为σ²的高斯白噪声。

有益效果：本发明在正则化方法框架中，通过引入洛伦兹函数构建数据保真项和正则化算子项，可以较好的解决图像的出界点问题，且图像的边缘特性保持较好；其次，采用了自适应正则化参数算法，对正则化参数进行自适应选择，避免了试错法选择参数的随机性。通过实验证明，本发明针对不同图像与不同噪声情况均能达到较好的重建质量，能较好的去除不同种类的噪声，保持边缘能力强，具有良好的适应性和鲁棒性。

附图说明

图1a和图1b是原始高分辨率图像，分别是国际标准测试图像“Cameraman”与一帧卫星遥感测试图像。

图2a～图2e是“Cameraman”图像在无噪声情况下的实验结果对比图。

图3a～图3e是“Cameraman”图像在高斯噪声情况下的实验结果对比图。

图4a～图4e是“Cameraman”图像在泊松噪声情况下的实验结果对比图。

图5a～图5e是“Cameraman”图像在椒盐噪声情况下的实验结果对比图。

图6a～图6e是卫星遥感测试图像在无噪声情况下的实验结果对比图。

图7a～图7e是卫星遥感测试图像在高斯噪声情况下的实验结果对比图。

图8a～图8e是卫星遥感测试图像在泊松噪声情况下的实验结果对比图。

图9a～图9e是卫星遥感测试图像在椒盐噪声情况下的实验结果对比图。

图10为本发明方法流程图。

具体实施方式

本发明方法公开了以下步骤：(1)输入数据为用向量表示的已知运动参数的m帧针对同一场景的低分辨率图像序列。

(2)根据运动参数对序列图像进行配准，并采用常见的双线性插值法对输入数据进行处理，获取初始高分辨率图像估计，即初值估计

(3)得到初值估计后，将其代入迭代公式进行迭代计算：

${\hat{Z}}_{n+1}={\hat{Z}}_n-\mathrm{\&beta;}\&CenterDot;G_n$

式中，β为自定义控制步长，n表示迭代次数，G_n表示梯度。

(4)根据式

在每次迭代中计算步骤(3)中的梯度：

$G_n=[\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{H}_k^T\&CenterDot;\mathrm{\&eta;}\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}(H_k{\hat{Z}}_n-Y_k)]\&CenterDot;\frac{1}{\mathrm{\&alpha;}}+{\mathrm{\&lambda;}}_k\&CenterDot;C^T\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}\left(C{\hat{Z}}_n\right)$

其中，

参数T₁是洛伦兹函数的半值半宽，用作自定义控制阈值；C表示拉普拉斯算子，C＝[1/8 1/8 1/8；1/8 -1 1/8；1/8 1/8 1/8]；α、η为自定义常数系数；H_k为等效降质因子；λ_k为自适应正则化参数；

(5)根据下式在每次迭代中计算步骤(4)中的自适应正则化参数λ_k：

${\mathrm{\&lambda;}}_k={\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}(Y_k-H_k{\hat{Z}}_n)\left|\right|}^2/(2{\left|\right|Y_k\left|\right|}^2-{\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}\left(C{\hat{Z}}_n\right)\left|\right|}^2);$

其中，

参数T₂是洛伦兹函数的半值半宽，用作自定义控制阈值，log表示对数函数，x表示输入变量。

(6)判断是否满足迭代收敛条件，若满足则停止迭代，输出所得高分辨率图像。

本发明首先要建立高分辨率图像的降质模型，从数学建模的角度刻画实际成像系统获取数字图像的过程，即高分辨率图像退化为低分辨率图像的过程。设景物图像是带限连续的，成像系统首先以奈奎斯特频率采样，得到理想高分辨率图像Z。设在成像过程中，只存在运动、降质因素、噪声影响，则高分辨率图像Z依次经过变形、模糊、下采样、添加噪声后得到低分辨率观测图像序列Y_k。图像观测模型可用下式表示：

Y_k＝D_kB_kW_kZ+V_k $\&ForAll;k=\mathrm{1,2},...,m$ ①

式中，Z是理想高分辨率(HR，High Resolution)图像，大小为qN₁×qN₂像素，

N＝qN₁×qN₂，是高分辨率图像像素以向量形式按字典顺序排列而成，表示单个像素，N表示向量Y_k中的分量个数，

为向量Z中的第N个分量，q为采样比率，。Y_k是低分辨率(LR，Low Resolution)图像序列，每一帧图像的大小都是N₁×N₂像素，Y_k＝[Y_k，1，Y_k，2，……Y_k，M]^T，k＝1，2，…，m，M＝N₁×N₂，是低分辨率图像像素以向量形式按字典顺序堆积而成。D_k表示下采样矩阵，矩阵大小为(N₁N₂)²×q²N₁N₂。B_k表示模糊矩阵，包括成像系统的光学模糊、运动模糊、传感器点扩散函数等因素，矩阵大小为q²N₁N₂×q²N₁N₂。W_k表示获取图像时发生的全局或局部几何变形，旋转等，矩阵大小是q²N₁N₂×q²N₁N₂。V_k表示成像系统的噪声影响，设为均值为零，方差为σ²的高斯白噪声。为简便起见，可用下式表示：

Y_k＝H_kZ+V_k $\&ForAll;k=\mathrm{1,2},...,m$ ②

其中H_k＝D_kB_kW_k，称为等效降质因子，概括了成像系统的降质因素。

图像超分辨重建过程是从低分辨率序列图像估计高分辨率图像的过程。由式①和式②可知，重建方程不满足Hadamard良态问题条件。即超分辨率重建时对输入的观测数据误差非常敏感，微小误差和少量噪声可能会导致重建图像与真实图像偏离甚远。本发明中采取正则化方法，添加解的限制条件来得到稳定解，即通过建立代价函数：

L_λ(Z)＝Ψ(Y-HZ)+λ·Φ(Z)       ③

使得：

$\hat{Z}=\underset{\hat{Z}}{\arg \min }{L}_{\mathrm{\&lambda;}}\left(\hat{Z}\right)$ ④

其中，Ψ(·)为某种误差范数，Φ(·)为某种正则化算子，λ是正则化参数，

是高分辨率重建图像估计。

式③中等号右边第一项为数据逼近项，用误差范数表示实际观测值与模型观测值之间的相似度的度量，引入高频能量保持细节，确保解逼近于真实解；第二项是正则化算子，引入先验知识约束解的范围，增加解的稳定性。λ是正则化参数，控制前述两项对解的贡献程度。正则化方法的难点在选择合适的误差范数以保证解的准确性；构造恰当的正则化算子(先验知识)以获得稳定性；确定合理的正则化参数以平衡重建图像的平滑程度与保持边缘能力。

本发明同时采用洛伦兹函数构造数据逼近项与正则化算子项。洛伦兹函数及其影响函数定义如下：

${\mathrm{\&Psi;}}_{\mathrm{lor}}\left(x\right)=\log [1+\frac{1}{2}{\left(\frac{x}{T}\right)}^2],$ ${{\mathrm{\&Psi;}}^{\&prime;}}_{\mathrm{lor}}\left(x\right)=\frac{2x}{2T^2+x^2}$ ⑤

其中，参数T是洛伦兹函数的半值半宽，用作控制阈值。

采用洛伦兹函数构建数据逼近项，其在图像的平坦区域，平滑能力强；在图像的边缘部分，对于出界点不敏感，具有良好的边缘保持能力。参数T值的改变将影响洛伦兹函数作为误差范数的作用。误差范数的作用是保真，保护高频能量以使近似解逼近于真实解。当T值增大时，函数峰值减小，函数形状变低变宽，意味着更加严厉的惩罚出界点，并抑制噪声。反之亦然。

本发明在构建正则化算子项时，考虑了真实图像的局部性质的差异。噪声对图像的不同特性区域的视觉影响是不同的：在图像的边缘附近方差较大，噪声带来的影响较小；在图像的平坦区域方差较小，噪声带来的不良视觉影响较大。

传统吉洪诺夫正则化算子(高斯先验)的构建是以图像的平滑性模型作为先验知识，即假设图像是平滑连续的。但是实际图像通常是不连续的，存在着诸多边缘突变，这些不连续信息包含了图像的重要内容。这种方法往往使得重建的图像丢失细节、模糊边缘，显得过于平滑。人的视觉恰恰对于边缘变化更加敏感，故此以平滑性为前提先验的正则化方法存在局限性。由于限制重建结果中的高频分量所带来的边缘振铃效应也是需要克服的。

洛伦兹先验是基于边缘图像的分布特性，能获得更好的边缘保持效果，且对于噪声的惩罚力度较大。Fuderer证明了边缘图像的直方图近似服从洛伦兹分布函数，该函数已成功运用于解决图像处理中的消噪及去振铃等问题。

式⑤中，T值代表洛伦兹函数的半值半宽。通过计算可知，T值越大，函数将越宽，峰值将越小。这表明通过调节T值大小，可以适应不同噪声分布特性的图像。若高分辨率图像具有丰富的细节，则应设置较小的T值，以保留更多的边缘细节。反之亦然。

本发明引入自适应方法确定正则化参数，针对图像特性采取不同的正则化策略，可以更好的平衡近似解的精确程度与不适定性。选取原则为根据图像特点，在图像平坦区域选择较大的正则化参数，以抑制噪声，对人的视觉影响不大的区域施加平滑性约束；在图像边缘细节丰富区域选择较小的正则化参数，以保持边缘并增加细节，保护重要信息不丢失。避免了依赖人为经验确定或需要事先通过大量计算得到的缺点。

由于引入先验知识时，洛伦兹函数是对边缘图像先验分布特性的建模，故本发明首先引入拉普拉斯算子提取图像的边缘，然后引入洛伦兹函数构建误差范数项及正则化算子。构建正则化代价函数：

$L_{\mathrm{\&lambda;}}\left(Z\right)=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&Psi;}}_{\mathrm{lor}}(Y_k-H_{k}Z)+\mathrm{\&lambda;}\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}\left(\mathrm{CZ}\right)$ ⑥

使得：

$\hat{Z}=\underset{\hat{Z}}{\arg \min }{L}_{\mathrm{\&lambda;}}\left(\hat{Z}\right)$ ⑦

则当代价函数最小时，得到最优估计。其中，Ψ_lor(·)为洛伦兹函数，λ是正则化参数，Z是高分辨率重建图像估计，C是拉普拉斯算子。

可得，本发明采用的自适应正则化参数为：

${\mathrm{\&lambda;}}_k={\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}(Y_k-H_k{\hat{Z}}_n)\left|\right|}^2/(2{\left|\right|Y_k\left|\right|}^2-{\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}\left(C{\hat{Z}}_n\right)\left|\right|}^2)$ ⑧

本发明采用最速下降法进行求解。最速下降法是一种基于梯度的迭代算法。基于梯度的迭代方法必须从解的一个初始估计开始。最速下降法是从初始估计

开始，沿性能曲面最陡的方向向下搜索曲面的最低点。曲面的最陡下降方向是曲面的负梯度方向，该迭代搜索过程首先从曲面上某个初始点出发，沿该点负梯度方向搜索至第1点，然后再从第一点开始，重复前面的过程，直至寻找到曲面的最低点。该方法的迭代公式为：

${\hat{Z}}_{n+1}={\hat{Z}}_n-\mathrm{\&beta;}\&CenterDot;G_n$ ⑨

其中，β为步长，G_n为梯度。

对代价函数求梯度，得

$G_n=\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{H}_k^T\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}(H_k{\hat{Z}}_n-Y_k)+{\mathrm{\&lambda;}}_k\&CenterDot;C^T\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}\left(C{\hat{Z}}_n\right)$ ⑩

之后将求得的梯度带入迭代公式求解。

为验证本发明性能，以下给出了针对两幅不同特性的图像分别在四种不同的噪声情况下的仿真实验结果，并与双线性插值算法(后文简称为双线性插值，是一种常规图像插值算法)、基于高斯先验的自适应正则化方法(后文简称高斯方法，是一种基于Schultz等人提出的MAP算法的常用改进方法；为方便实验比较，本发明对其进行了改进，加入了自适应正则化参数算法)进行了对比，来阐述本发明的技术效果(本发明方法简称为洛伦兹方法)。

实验中，图1a是原始高分辨率图像，是国际标准测试图像“Cameraman”，图2a～图5e为国际标准测试图像，且具有大量的平坦区域。图1b是原始高分辨率图像，是一帧卫星遥感测试图像，图6a～图9e为卫星遥感图像，且具有丰富的细节与边缘突变。每幅图像模拟四种噪声情况，分别为无噪声情况、高斯噪声情况、泊松噪声情况、椒盐噪声情况。

图2a～图2e，分别是“Cameraman”图像在无噪声情况下的实验结果对比图。其中图2a为低分辨率图像，图2b为采用双线性插值方法得出的高分辨率图像，图2c为采用高斯方法得出的高分辨率图像，图2d为采用洛伦兹方法得出的高分辨率图像，图2e为三种重建方法的评价指标与迭代次数的关系对比图。从图2b～图2d的重建图像上看，洛伦兹方法重建图像2d最为清晰，边缘保持较好；高斯方法重建图像2c由于平滑过度，稍显模糊；双线性插值方法重建图像2b则模糊不清。从图2e的曲线可见，洛伦兹方法曲线在达到最优值后仍然平稳，高斯方法曲线迅速达到最优值，但随着迭代次数的增加，评价指标值劣化明显。表明洛伦兹方法具有较强的稳定性。

图3a～图3e，分别是“Cameraman”图像在高斯噪声情况下的实验结果对比图。其中图3a为低分辨率图像，图3b为采用双线性插值方法得出的高分辨率图像，图3c为采用高斯方法得出的高分辨率图像，图3d为采用洛伦兹方法得出的高分辨率图像，图3e为三种重建方法的评价指标与迭代次数的关系对比图。从图3b～图3d的重建图像上看，洛伦兹方法重建图像3d稍为清晰，边缘保持较好；高斯方法重建图像3c则略为模糊；双线性插值方法重建图像3b则模糊不清，人物背景中布满噪声。从图3e的曲线可见，洛伦兹方法曲线稳定在高斯方法曲线上方，表明该算法优于高斯方法，且输出值稳定。两条曲线在达到最优值以后，随着迭代次数的增加，均较为稳定，没有出现明显的劣化现象。

图4a～图4e，分别是“Cameraman”图像在泊松噪声情况下的实验结果对比图。其中图4a为低分辨率图像，图4b为采用双线性插值方法得出的高分辨率图像，图4c为采用高斯方法得出的高分辨率图像，图4d为采用洛伦兹方法得出的高分辨率图像，图4e为三种重建方法的评价指标与迭代次数的关系对比图。从图4b～图4d的重建图像上看，洛伦兹方法重建图像4d清晰度较高，消除了图像中人物背景部分的大部分噪声，且人物细节边缘保持较好；高斯方法重建图像4c则比较模糊，且残留噪声较明显。双线性插值图像4b中人物更加模糊不清，噪声也较明显。从图4e的曲线可见，两种方法达到最优值的迭代次数相近，且继续迭代计算并没有造成指标值的下降，较为稳定。

图5a～图5e，分别是“Cameraman”图像在椒盐噪声情况下的实验结果对比图。其中图5a为低分辨率图像，图5b为采用双线性插值方法得出的高分辨率图像，图5c为采用高斯方法得出的高分辨率图像，图5d为采用洛伦兹方法得出的高分辨率图像，图5e为三种重建方法的评价指标与迭代次数的关系对比图。从图5b～图5d的重建图像上看，洛伦兹方法重建图像5d清晰度较高，基本完全消除了图像中的椒盐噪声，且人物细节边缘保持较好；高斯方法重建图像5c则由于过于强调平滑导致重建图像比较模糊，在去除噪声的同时，也损失了许多图像的细节。双线性插值图像5b质量基本无提高。从图5e的曲线可见，洛伦兹方法曲线与高斯方法曲线在垂直方向上相隔的距离表明采用洛伦兹方法较高斯方法对指标值的改善更加明显。洛伦兹方法曲线在达到最优值后仍然平稳，高斯方法曲线迅速达到最优值，但随着迭代次数的增加，评价指标值劣化明显。表明洛伦兹方法具有较强的稳定性。

图6a～图6e，分别是卫星遥感测试图像在无噪声情况下的实验结果对比图。其中图6a为低分辨率图像，图6b为采用双线性插值方法得出的高分辨率图像，图6c为采用高斯方法得出的高分辨率图像，图6d为采用洛伦兹方法得出的高分辨率图像，图6e为三种重建方法的评价指标与迭代次数的关系对比图。从图6b～图6d的重建图像上看，洛伦兹方法重建图像6d最为清晰，边缘保持较好；高斯方法重建图像6c由于平滑过度，稍显模糊；双线性插值图像6b则模糊不清。从图6e的曲线可见，洛伦兹方法曲线在达到最优值后仍然平稳，高斯算法曲线虽然很快达到最优值，但随着迭代次数的增加，评价指标值出现明显的劣化。表明洛伦兹方法具有较强的稳定性。本次实验与“Cameranman”图像实验中相同噪声情况相比，高斯方法收敛速度提高5.6倍，洛伦兹方法提高22倍。再比较两次实验的PSNR值，高斯方法与双线性插值相比两次分别提高2.61倍、2.58倍；洛伦兹方法与双线性插值相比两次分别提高42倍、3.7倍。可知，卫星遥感测试图像实验提高的倍数略小，但相差并不大。两次实验参数设置几乎相同，表明收敛速度与图像特性有关。两图对比可见，“Hohai”图像具有更多的边缘与细节，“Cameraman”图像则具有大量的平坦区域。表明在无噪声情况下，针对具有较多细节的图像，洛伦兹方法具有更佳的重建效果。

图7a～图7e，分别是卫星遥感测试图像在高斯噪声情况下的实验结果对比图。其中图7a为低分辨率图像，图7b为采用双线性插值方法得出的高分辨率图像，图7c为采用高斯方法得出的高分辨率图像，图7d为采用洛伦兹方法得出的高分辨率图像，图7e为三种重建方法的评价指标与迭代次数的关系对比图。从图7b～图7d的重建图像上看，洛伦兹方法重建图像7d最清晰，边缘保持较好；高斯方法重建图像7c则略为模糊；，双线性插值图像7b则模糊不清。从图7e的曲线可见，高斯方法与洛伦兹方法达到最优值后，但随着迭代次数的增加，指标值出现劣化。而在“Cameraman”图像实验中相同噪声情况下，迭代曲线并无明显劣化。表明针对具有较多细节的图像重建中，细节成分可能被误认为噪声成分，随着迭代次数的增加被逐渐清除，从而造成重建图像质量的下降。这是由于高斯方法在重建时过度强调平滑约束的结果。

图8a～图8e，分别是卫星遥感测试图像在泊松噪声情况下的实验结果对比图。其中图8a为低分辨率图像，图8b为采用双线性插值方法得出的高分辨率图像，图8c为采用高斯方法得出的高分辨率图像，图8d为采用洛伦兹方法得出的高分辨率图像，图8e为三种重建方法的评价指标与迭代次数的关系对比图。从图8b～图8d的重建图像上看，洛伦兹方法重建图像8d清晰度较高，各房屋的边缘较为清晰；高斯方法重建图像8c则比较模糊，且有少许重影；双线性插值图像8b则模糊不清。从图8e的曲线可见，两条曲线出现了交叉，交叉点恰好在高斯方法取得最优值附近。表明高斯方法首先取得最优值。但随着迭代次数的增加，洛伦兹方法能够继续取得更优的指标值。

图9a～图9e，分别是卫星遥感测试图像在椒盐噪声情况下的实验结果对比图。其中图9a为低分辨率图像，图9b为采用双线性插值方法得出的高分辨率图像，图9c为采用高斯方法得出的高分辨率图像，图9d为采用洛伦兹方法得出的高分辨率图像，图9e为三种重建方法的评价指标与迭代次数的关系对比图。从图9b～图9d的重建图像上看，洛伦兹方法重建图像9d清晰度较高，基本完全消除了图像中的椒盐噪声，可较好的分辨房屋与道路的细节边缘；高斯方法重建图像9c则由于过于强调平滑导致重建图像比较模糊，在去除噪声的同时，也损失了许多图像的细节；双线性插值图像9b基本没有图像质量上的提高。从图9e的曲线可见，与前次实验相比，曲线对比情况类似，表明采用洛伦兹方法较高斯方法对指标值的改善更加明显，且重建性能更为稳定。

以下结合实验和附图对本发明的具体实施方式进行说明。

(1)首先将原始高分辨率图像模拟全局平移，假设运动参数已知；然后进行高斯模糊，并进行2倍的下采样，最后添加噪声，生成10帧低分辨率图像。为充分论证本发明算法的鲁棒性及对不同情况的适应性，分别对低分辨率图像添加高斯噪声、泊松噪声与椒盐噪声，并分次进行实验。

(2)根据运动参数对序列图像进行平移配准，采用双线性插值法获取低分辨率图像的对应的初始高分辨率图像估计，即初值估计

(3)根据迭代公式

进行迭代计算；

(4)根据式

在每次迭代中计算梯度；

(5)根据自适应正则化参数公式

在每次迭代中计算λ_k；

(6)判断是否满足迭代收敛条件，若满足则停止迭代，输出所得高分辨率图像。

关于实验中所涉及的参数：误差范数中洛伦兹函数T值设为T₁，正则化算子中洛伦兹函数T值设为T₂。α、η为常数系数，根据经验取值。步长β控制收敛速度的快慢，步长的选择通过试错法来决定。为充分验证其收敛性能，收敛条件设定为某迭代次数。

通过实验，证明本发明算法的重建质量无论从主观感受上还是客观指标上都好于高斯方法，重建图像的清晰度与边缘保持能力较好。其次，本发明算法具有较强的适应性与鲁棒性，在多种情况下均能保持良好的性能，得到较好的结果。尤其在椒盐噪声情况下，几乎能够完全消除噪声，重建图像质量接近原始高分辨率图像。本发明优越性在于对正则化参数进行自适应选择，无需按传统方法处理候选值，避免了试错法选择参数的随机性。

本发明采用国际常用的图像评价标准：峰值信噪比PSNR与均方根误差RMSE。PSNR值单位为dB(分贝)，数值越大表示图像质量越好，RMSE值单位为像素，数值越小表示图像质量越好。

表1和表2中具体描述了针对同一幅原始高分辨图像，在4种不同的噪声情况下，分别采用3种不同实验方法得出的超分辨率重建图像之间的质量比较；同时给出了相关实验参数。

由表可见，在同种噪声情况下，洛伦兹方法的PSNR值与RMSE值均优于另2种方法。特别的，在椒盐噪声的情况下，洛伦兹方法的PSNR值与RMSE值大幅优于另2种方法。对比表明，洛伦兹方法对于不同类型的图像，不同类型的噪声均具有较强的适应性与鲁棒性，且图像重建质量较好。

表1给出了采用本发明方法和其他常用方法对国际标准测试图像“Cameraman”的实验关键参数指标值与实验结果客观评价指标的比较结果。

表1

表2给出了采用本发明方法和其他常用方法对卫星遥感测试图像的实验关键参数指标值与实验结果客观评价指标的比较结果。

表2

本发明提供了一种图像超分辨率重建方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (6)

1.一种图像超分辨率重建方法，其特征在于，包括如下步骤：

步骤(1)输入用向量表示的m帧针对同一场景的低分辨率图像序列；

步骤(2)采用双线性插值法对所有低分辨率图像进行处理，获取初始高分辨率图像估计，即初值估计

步骤(3)将初值估计

代入迭代公式进行迭代计算，其中，β为设定控制步长，设定值范围为0～1，n表示迭代次数，

表示经过第n次迭代的高分辨率图像估计，G_n表示经过第n次迭代计算得到的梯度；

步骤(4)判断是否满足迭代收敛条件，若满足则停止迭代，输出所得高分辨率图像。

2.根据权利1要求所述的图像超分辨率重建方法，其特征在于，根据下式在每次迭代中计算梯度G_n，

$G_n=[\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{H}_k^T\&CenterDot;\mathrm{\&eta;}\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}(H_k{\hat{Z}}_n-Y_k)]\&CenterDot;\frac{1}{\mathrm{\&alpha;}}+{\mathrm{\&lambda;}}_k\&CenterDot;C^T\&CenterDot;{\mathrm{\&Psi;}}_{\mathrm{lor}}^{\&prime;}\left(C{\hat{Z}}_n\right),$

其中，

参数T₁是洛伦兹函数的半值半宽，为控制阈值，设定值范围为1～12之间的实数；C表示拉普拉斯算子，C＝[1/8 1/8 1/8；1/8 -1 1/8；1/8 1/8 1/8]；α和η为常数系数，α设定值范围为0～100之间的实数，η设定值范围为100～500之间的实数；H_k为等效降质因子；λ_k为自适应正则化参数。

3.根据权利2要求所述的图像超分辨率重建方法，其特征在于，根据下式在每次迭代中计算自适应正则化参数λ_k：

${\mathrm{\&lambda;}}_k={\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}(Y_k-H_k{\hat{Z}}_n)\left|\right|}^2/(2{\left|\right|Y_k\left|\right|}^2-{\left|\right|{\mathrm{\&Psi;}}_{\mathrm{lor}}\left(C{\hat{Z}}_n\right)\left|\right|}^2);$

其中，

参数T₂是洛伦兹函数的半值半宽，参数T₂设定值范围为1～12之间的实数。

4.根据权利1或3要求所述的图像超分辨率重建方法，其特征在于，步骤(1)中输入的低分辨率图像序列是满足如下条件建立高分辨率图像的降质模型：针对带限连续的景物图像，成像系统以奈奎斯特频率采样，得到理想高分辨率图像Z；在成像过程中，只存在运动、降质因素、噪声影响；高分辨率图像Z依次经过变形、模糊、下采样、添加噪声后得到低分辨率观测图像序列Y_k。

5.根据权利4要求所述的图像超分辨率重建方法，其特征在于，图像降质模型用下式表示：

对单帧高分辨率图像进行降质处理，依次进行运动平移、模糊、降采样、添加噪声；重复该步骤，获取m帧关于同一景物的不同的低分辨率图像，设H_k＝D_kB_kW_k，称为等效降质因子；

所述Z是理想高分辨率图像，大小为qN₁×qN₂像素，Z由高分辨率图像像素以向量形式按字典顺序排列而成，即

N＝qN₁×qN₂，

为向量Z中的第N个分量，N表示向量Y_k中的分量个数，q为采样比率；

Y_k是低分辨率像序列，每一帧图像的大小都是N₁×N₂像素，Y_k由低分辨率图像像素以向量形式按字典顺序排列而成，即Y_k＝[Y_k，1，Y_k，2，……Y_k，M]^T，k＝1，2，…，m，M＝N₁×N₂，M表示向量Y_k中的分量个数；

D_k表示下采样矩阵，矩阵大小为(N₁N₂)²×q²N₁N₂；

B_k表示模糊矩阵，矩阵大小为q²N₁N₂×q²N₁N₂；

W_k表示获取图像时的全局或局部几何变形，旋转，矩阵大小是q²N₁N₂×q²N₁N₂；

V_k表示成像系统的噪声影响，设为均值为零，方差为σ²的高斯白噪声。

6.根据权利1要求所述的图像超分辨率重建方法，其特征在于，步骤(4)中所述收敛条件ρ设定为 $\mathrm{\&rho;}={\left|\right|{\hat{Z}}_{n+1}-{\hat{Z}}_n\left|\right|}^2/{\left|\right|{\hat{Z}}_n\left|\right|}^2<5\&times;{10}^{-7}.$

Images