CN103824273A - 基于复合运动和自适应非局部先验的超分辨率重建方法 - Google Patents

Abstract

本发明提供一种基于复合运动和自适应非局部先验的超分辨率重建方法，包括如下步骤：在p幅低分辨率图像中选取参考帧图像和非参考帧图像；采用全局运动参数和局部光流的方法进行图像配准，得到非参考帧图像相对于参考帧图像的运动场m_k(x)，利用m_k(x)构建出运动变换矩阵M_k；计算参考帧图像的插值图像、非局部先验参数h_i，j和欧式阈值；计算每个像素与其他像素的相似度权重w_NLM[i，j；s，t]，利用w_NLM构建关于高分辨率图像X的非局部权重矩阵S；利用运动变换矩阵M_k和非局部权重矩阵S求解目标泛函得到重建的高分辨率估计图像。相对于现有技术，本发明通过采用复合运动模型，有效地解决了目前运动估计的计算量大，可伸缩性不强、精度不高的缺点，采用自适应的非局部先验减少了重建图像的失真。

Description

基于复合运动和自适应非局部先验的超分辨率重建方法

技术领域

本发明涉及图像处理、计算机视觉的技术领域，特别是涉及图像清晰化领域中的一种提高图像的空间分辨率的基于复合运动和自适应非局部先验的超分辨图像重建方法。

背景技术

在成像领域中，高空间分辨率的图像一直是追求的目标之一。高空间分辨率的图像充分地记录了物体的细节信息，能够给人和计算机的推理、判断、决策提供更加丰富的信息。因此，在许多成像应用中，高分辨率的图像通常是非常必要的，比如：视频监控、医学诊断、军事侦察、遥感等应用。提高图像的空间分辨率可以通过两种途径，“硬件途径”和“软件途径”。常规情况下，人们主要通过改进高精度的CCD和CMOS传感器等硬件设备来获取高分辨的图像。但是，单纯地通过改进硬件设施来提高分辨率会受到诸多限制，比如传感器的电荷转移率、热噪声，光学镜头的瑞利熵，以及硬件费用等限制。因此，出于硬件设备和经济上的考虑，“软件途径”成为了一种更加可行的方案。1984年Tsai和Huang首先提出超分辨重建问题便是从“软件途径”，通过发展理论、算法来提高图像的空间分辨率，已成为图像处理领域最为活跃的研究方向之一。近年来，一些学者将超分辨率重建分为基于多帧的超分辨率重建和基于学习的超分辨率重建。本发明是一种基于多帧的超分辨重建方法。

对基于多帧的超分辨率重建又可以分为基于频域和基于空域的重建方法。最新的基于频域的重建方法是1999年Rhee和Kang提出的基于离散余弦变换重建方法，以及Chan在2003年提出的基于小波变换的重建方法。基于频域的重建方法优势是理论简洁、计算简单，但劣势在于多帧之间的相对运动只能是全局相似运动，且只适用于图像的模糊是线性时不变的情况。因此针对这些劣势，基于空域的重建方法产生了许多经典的超分辨重建方法，如非均匀内插法、迭代反投影法、凸集投影法、最大似然法、最大后验法、混合最大后验／凸集投影法。本发明是一种基于最大后验估计的重建方法。

在基于最大后验的估计方法中，运动估计和先验项设计是两项非常重要的任务。在大多数情况下，图像超分辨率所需的运动矢量为未知，为了求解这些运动矢量，可以有两种思路：一种思路是先求解这些运动矢量然后再进行超分辨率重建，虽然这种分开求解的方式计算非常简单，但是都各自存在着很大的局限性；另一种思路是运动矢量和超分辨图像进行联合求解，这类方法比单独分开求解方法的运动估计精确，重建效果好，但缺点是求解速度慢，很难用于实际应用。针对运动矢量的估计，又可分为基于参数模型的和基于光流的运动估计。基于参数模型的运动估计比较简单，但伸缩性低；而基于光流的运动估计伸缩性高，但估计精度低，重建效果不理想。为此，本发明基于分开求解运动矢量和超分辨图像的思路，提出了一种具有很高伸缩性和估计精度的、低运算量的复合运动模型。

针对先验项，也即正则项设计，国内外文献中已经提出了许多先验项，包括Tikhonov先验、Huber先验、TV先验、BTV先验等。这些先验都是基于图像邻域差分进行描述的，未能准确地描述出自然图像的先验信息，导致重建图像失真。最近，基于自然图像中存在着大量冗余重复的图像结构的事实，国际上提出了一种非局部先验，并成功应用于图像反卷积。但是，由于非局部先验中存在做许多需要人为调节的参数，不能做到参数的自适应。

发明内容

针对现有技术存在的问题，本发明的主要目的在于提供一种降低计算量、提高图像精度的基于复合运动和自适应非局部先验的超分辨率重建方法。

为实现上述目的，本发明提供了一种基于复合运动和自适应非局部先验的超分辨率重建方法的实施例，高分辨率图像X经过降质过程获取p幅低分辨率观测图像Y_k(k＝，…，p)，每幅观测图像的大小为m×n，该利用p幅低分辨率观测图像Y_k(k＝1，…，p)重建高分辨率估计图像的超分辨率重建方法包括如下步骤(1)至步骤(5)：

(1)在p幅低分辨率观测图像中选取参考帧图像Y_ref(1≤ref≤p)和非参考帧图像Y_k(k＝1，…，ref-1，ref+1，…，p)，针对参考帧图像和非参考帧图像之间的亚像素运动采用全局参数运动和局部光流的复合运动模型，参考帧图像与非参考帧图像之间的关系表示为： $Y_{\mathrm{ref}}\left(x\right)=Y_k\left(m_k\left(x\right)\right)=Y_k(m_k^g(x;{\mathrm{\θ}}_k)+m_k^l\left(x\right))=Y_k\left(m_k^g(x;{\mathrm{\θ}}_k)\right)+{\mathrm{\ϵ}}_k\left(x\right),$ 其中m_k表示二维运动场，

表示全局参数运动，

为局部光流运动，θ_k为全局运动参数，

表示用非参考帧图像预测的参考帧图像，ε_k(x)表示残差图像；

(2)求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)和局部光流

，采用全局运动参数θ_k＝(a₀，a₁，₂，a₃，a₄，a₅)和局部光流

的方法进行图像配准，得到非参考帧图像相对于参考帧图像的运动场m_k(x)，利用m_k(x)构建出运动变换矩阵M_k；

(3)计算参考帧图像Y_ref的r倍插值图像、非局部先验参数_i，j(0≤i<rm，0≤j<rn)和相似图像的欧式阈值；

(4)利用非局部先验参数、欧式阈值，且以插值图像

作为高分辨图像X的初始图像，计算出插值图像

中每个像素点(i，j)(其中0≤i<rm，0≤j<rn)与其他像素点(s，t)(其中0≤s<rm，0≤t<rn)的相似度权重w_NLM[i，j；s，t]，利用相似度权重w_NLM构建非局部权重矩阵S；

(5)利用运动变换矩阵M_k和非局部权重矩阵S求解目标泛函 $X=\arg \min [\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|(I-S)\left|\right|}_{\mathrm{\&rho;}}^{\mathrm{\&rho;}}],$ 其中B_k为观测图像Y_k对应的降晰函数，M_k为观测图像Y_k相对于参考帧图像的亚像素运动，非局部权重矩阵S是一个自适应高分辨率图像X的非局部均值滤波器，且ρ＞0，采用共轭梯度迭代法最小化目标泛函，得到重建的高分辨率估计图像。

进一步地，该步骤(1)针对参考帧图像和非参考帧图像之间的亚像素运动采用全局参数运动和局部光流的复合运动模型，其中运动向量为m_k(x)＝[m_k，u(x)m_k，v(x)]，x＝[x_ux_v]，二维运动场m_k表示为m_k(x)＝m_k ^g(x)+m_k ^l(x)＝m_k ^g(x；θ_k)+m_k ^l(x)，其中

表示全局参数运动，

为局部光流运动，θ_k为全局运动参数。

更进一步地，该步骤(2)求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)包括如下步骤(21A)至步骤(22A)：

(21A)采用a₀，a₁，a₂，a₃，a₄，a₅六个参数的仿射变换作为全局参数运动模型： $m_{k,u}^g(x;{\mathrm{\&theta;}}_k)=a_0+a_1{x}_u+a_2{x}_v,m_{k,v}^g(x;{\mathrm{\&theta;}}_k)=a_3+a_4{x}_u+a_5{x}_v;$

(22A)建立最小二乘标准

求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)。该步骤(22A)建立最小二乘标准

求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)包括如下步骤(221A)至步骤(224A)：(221A)将θ_k写成θ_k+Δθ；(222A)对最小二乘目标泛函进行泰勒展开得到关于Δθ的线性函数；(223A)对展开后的函数进行一系列运算操作得到 $\mathrm{\&Delta;\&theta;}=H^{-1}\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T[Y_{\mathrm{ref}}\left(x\right)-Y_k\left(m_k^g(x;{\mathrm{\&theta;}}_k)\right)],$ 其中 $H=\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{\&PartialD;m_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T\underset{x}{\mathrm{\&Sigma;}}[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}],$ $\&dtri;Y_k=\frac{{\&PartialD;Y}_k}{{\&PartialD;m}_k^g(x;{\mathrm{\&theta;}}_k)},$ $\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}=\left(\begin{array}{cccccc}{x}_u& 0& x_v& 0& 1& 0\\ 0& x_u& 0& x_v& 0& 1\end{array}\right);$ (224A)判断是否符合||Δθ||≤∈，阈值0≤∈≤0.01，若不符合，则θ_k←θ_k+Δθ，并返回步骤(223A)；若符合，则表明参数θ_k收敛，停止迭代，得到全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)，最后可以得到全局参数运动

更进一步地，该步骤(2)求解局部光流包括如下步骤(21B)至步骤(24B)：

(21B)利用图像灰度恒常性假设和图像梯度恒常性假设得到关于局部光流

的数据置信度能量函数 $E_{\mathrm{Data}}(u,v)=\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&Psi;}({|Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-Y_k\left(\stackrel{\&OverBar;}{x}\right)|}^2+|\&dtri;Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-\&dtri;Y_k{|}^2)d\stackrel{\&OverBar;}{x},$ 其中， $\stackrel{\&OverBar;}{x}={(x,t)}^T={(x_u,x_v,1)}^T,v={(m_k^g{\left(x\right)}^T,0)}^T,\mathrm{\&omega;}={(m_k^l{\left(x\right)}^T,1)}^T={(u,v,1)}^T,$ 函数 $\mathrm{\&Psi;}\left(s^2\right)=\sqrt{s^2+{\mathrm{\&tau;}}^2},$ 其中0<τ≤0.01；

(22B)根据图像分段平滑性假设得到平滑惩罚函数 $E_{\mathrm{Smooth}}(u,v)=\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&Psi;}\left({|\&dtri;u|}^2{+|\&dtri;v|}^2\right)d\stackrel{\&OverBar;}{x};$

(23B)得到整个能量函数为E(u，v)＝E_Data+αE_Smooth，其中α>0为正则化参数；

(24B)采用非线性数值求解方法对整个能量函数E(u，v)＝E_Data+αE_Smooth求最优值便得到局部光流

的解。

更进一步地，该步骤(2)利用m_k(x)构建出运动变换矩阵M_k包含如下步骤(21C)至步骤(23C)：

(21C)计算非参考帧图像相对于参考帧图像的运动场m_k(x)＝m_k ^g(x)+m_k ^l(x)，并将运动场m_k(x)进行r倍线性插值，得到插值后的运动场m_k(x)，该r为最终重建的高清图像的放大倍数；

(22C)计算相对位移Δx_k＝m_k(x)-x＝(Δx_k，u-Δx_k，v)，以及d_k＝Δx_k，u-floor(Δx_k，u)，e_k＝Δx_k，v-floor(Δx_k，v)，其中操作符floor(.)表示取小于或者等于指定值的最大整数；

(23C)计算运动变换矩阵M_k中的每个元素的值：M_k(j*m+i，floor(Jx_k，u)+x_u+ceil(Jx_k，v+x_v+1)*m)＝d_k*(1-e_k)，M_k(j*m+i，ceil(Δx_k，u)+xu+ceil(Δx_k，v+x_v+1)*m)＝d_k*e_k，M_k(j*m+ⁱ，floor(Δx_k，u)+x_u+floor(Δx_k，v+x_v+1)*m)＝(1-d_k)*(1-e_k)，M_k(j*m+i，ceil(Jx_k，u)+x_u+floor(Jx_k，v+x_v+1)*m+x)＝(1-d_k)*e_k，其中操作符ceil(.)表示取大于或者等于指定值的最小整数，对于矩阵M_k的第j*m+i行中的除这四列以外的其他元素值为零。

更进一步地，该步骤(3)计算非局部先验参数h_i，j，

其中std(N_i，j)是以像素点(i，j)为中心，搜索区域为N_i，j的标准差，β为大于零的常数(1<β<5)，σ²为参考帧图像Y_ref的噪声方差，r为图像放大倍数，该参考帧图像Y_ref的噪声方差σ²估计为： ${\widehat{\mathrm{\&sigma;}}}^2=\frac{1}{n-2}\underset{i=2}{\overset{n-1}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&epsiv;}}}_i^2,$ 其中n为图像的像素点总数， ${\widehat{\mathrm{\&epsiv;}}}_i=\sqrt{\frac{4}{5}}(Y_{\mathrm{ref}}\left(x_i\right)-\frac{1}{4}\underset{x_j\&Element;N\left(x_i\right)}{\mathrm{\&Sigma;}}{Y}_{\mathrm{ref}}\left(x_j\right)),$ N(x_i)为像素点x_i的四邻域。

更进一步地，该步骤(3)中相似图像的欧式阈值采用的最大加权欧式距离为4σ²，也即： ${\left|\right|R_{k,l}\stackrel{.}{X}-R_{i,j}\stackrel{.}{X}\left|\right|}_{2,a}^2\&le;4{\mathrm{\&sigma;}}^2\&ap;4{\widehat{\mathrm{\&sigma;}}}^2.$

更进一步地，该步骤(4)利用非局部先验参数、欧式阈值，且以插值图像

作为高分辨图像X的初始图像，计算出插值图像

中每个像素点(i，j)与其他像素点(s，t)的相似度权重w_NLM[i，j；s，t]，利用相似度权重w_NLM构建非局部权重矩阵S包含如下步骤(41)至步骤(43)：

(41)计算每个以坐标(i，j)为中心的像素块与其邻域N(i，j)内以坐标(s，t)为中心的像素块的相似度

其中R_i，j是以(i，j)为中心像素点的图像块提取算子，表示两个图像块的加权欧式距离，其中a>0是高斯核函数的标准差，h_i，j为滤波器平滑参数，依赖于图像的噪声大小以及图像本身，f为依赖于两个中心像素点的几何距离的常函数；

(42)归一化像素(i，j)与邻域N(i，j)内的所有像素相似度值， $\mathrm{\&omega;}[i.j;s,t]=\frac{{\mathrm{\&omega;}}_{\mathrm{NLM}}[i,j;s,t]}{{\mathrm{\&Sigma;}}_{(s,t)\&Element;N(i,j)}{\mathrm{\&omega;}}_{\mathrm{NLM}}[i,j;s,t]};$

(43)构建非局部权重矩阵 $S_{j*m+i,s*m+t}=\left\{\begin{array}{cc}\mathrm{\&omega;}[i,j;s,t]& (s,t)\&Element;N(i,j)\\ 0& \mathrm{otherwise}\end{array}\right..$

更进一步地，该步骤(5)中当

时，迭代终止，其中n为迭代次数。

相对于现有技术，首先，本发明首次将非局部先验应用于图像的多帧超分辨率重建，提出了一种非局部先验的参数自适应求解方法；其次，通过采用复合运动模型，有效地解决了目前运动估计的计算量大，可伸缩性不强、精度不高的缺点；再次，采用自适应的非局部先验更能准确地自动地描述图像的先验信息，减少了重建图像的失真。为此，该方法做到了可伸缩的高精度的运动估计和准确的自适应的先验描述，能够更有效地应用于实际的图像特征清晰化工程。

附图说明

图1为本发明的基于复合运动和自适应非局部先验的超分辨率重建方法的实施例的流程图

图2为本发明的求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)的实施例的流程图

图3为本发明的求解局部光流

的实施例的流程图

具体实施方式

下面结合附图，详细说明本发明的具体实施方式。

图像的超分辨率重建，即将观测图像恢复成理想图像。观测图像即一系列的低分辨率图像，理想图像即所求的高分辨率图像。给定一定场景的高分辨率图像X，经过一系列几何运动、光学模糊、亚采样以及附加噪声的降质过程，产生P幅低分辨率观测图像Y_k，用一个常用的图像观测模型描述理想图像与观测图像之间的关系，此观测模型为：Y_k＝DB_kM_kX+n_k，k＝1，…，p，其中M_k为运动变化矩阵，B_k为模糊矩阵，D为下采样矩阵，n_k为附加噪声。

基于上述观测模型，本发明采用最大后验方法对理想图像X进行估计。p幅低分辨率图像表示为Y＝[Y₁ ^T，Y₂ ^T，…，Y_p ^T]^T，基于最大后验估计理论的超分辨重建问题表示为：X＝argmaxP(X|Y)，经过简单的运算操作，可以得到： $X=\arg \max P\left(Y\right|X)P\left(X\right)=\arg \max \underset{k=1}{\overset{P}{\mathrm{\&Sigma;}}}P\left(Y_k\right|X)P\left(X\right).$ 此处P(Y_k|X)＝P(n_k)表示了观测模型噪声的类型，通常的假设噪声是均值为0，方差为σ_k ²的高斯噪声，即 $P\left(Y_k\right|X)=\frac{1}{C_1}\exp \{-\frac{{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2}{{{2\mathrm{\&sigma;}}_k}^2}\},$ 其中C₁为常数。图像先验概率密度一般形式为：

其中C₂为常数，η为控制参数，U(X)为关于图像X的先验能量函数。经过简单的化简，最终可以得到如下的代价函数： $X=\arg \min [\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2+\mathrm{\&lambda;U}\left(X\right)],$ 其中

在此，本发明假设图像的模糊核已知，也即模糊矩阵B_k已知。

根据 $X=\arg \min [\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2+\mathrm{\&lambda;U}\left(X\right)]$ 可知，要获得理想图像X，需要求解运动变换矩阵M_k和先验能量U(X)，即图像配准和自适应非局部先验设计。最终采用常用的共轭梯度法对目标泛函 $X=\arg \min [\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2+\mathrm{\&lambda;U}\left(X\right)]$ 进行求解获得理想图像X。

因此，本发明主要包括：(1)在多帧低分辨率图像配准时采用全局运动参数加局部光流的方法进行图像配准，并构建运动变换矩阵M_k；(2)设计自适应非局部先验，构建非局部权重矩阵S；(3)建立目标泛函，并采用共轭梯度法进行求解。

如图1所示，为本发明的基于复合运动和自适应非局部先验的超分辨率重建方法的实施例的流程图。超分辨率重建方法的实施例包括如下步骤S1至步骤S8：

S1、获取p幅由高分辨率图X经过降质过程形成的低分辨率观测图像Y_k，每幅观测图像的大小为m×n。得到图像观察模型Y_k＝DB_kM_kX+n_k，k＝1，…，p，其中M_k为运动变化矩阵，B_k为模糊矩阵，D为下采样矩阵，n_k为附加噪声，X为高分辨率图像，k为获取的p幅低分辨率图像的编号。

S2、在p幅低分辨率图像中选取参考帧图像Y_ref(1≤ref≤p)和非参考帧图像Y_k(k＝1，…，ref-1，ref+1，…，p)，针对参考帧图像和非参考帧图像之间的亚像素运动采用全局参数运动加局部光流的复合运动模型，参考帧图像与非参考帧图像之间的关系表示为： $Y_{\mathrm{ref}}\left(x\right)=Y_k\left(m_k\left(x\right)\right)=Y_k(m_k^g(x;{\mathrm{\&theta;}}_k)+m_k^l\left(x\right))=Y_k\left(m_k^g(x;{\mathrm{\&theta;}}_k)\right)+{\mathrm{\&epsiv;}}_k\left(x\right),$ 其中m_k表示二维运动场，

表示全局参数运动，

为局部光流运动，也即非参数局部补偿运动，θ_k为全局运动参数，

表示用非参考帧图像预测的参考帧图像，ε_k(x)表示残差图像。在全局参数运动和局部光流的复合运动模型中，运动向量为m_k(x)＝[m_k，u(x)m_k，v(x)]，估计此运动向量的目的是用来构建运动变换矩阵M_k，其中x=[x_ux_v]，二维运动场m_k表示为 $m_k\left(x\right)=m_k^g\left(x\right)+m_k^l\left(x\right)=m_k^g(x;{\mathrm{\&theta;}}_k)+m_k^l\left(x\right),$ 其中

表示全局参数运动，

为局部光流运动，θ_k为全局运动参数。

s3、求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)和局部光流

采用全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)和局部光流

的方法进行图像配准，得到非参考帧图像相对于参考帧图像的运动场m_k(x)，利用m_k(x)构建出运动变换矩阵M_k：

求解全局运动参数θ_k。(a₀，a₁，a₂，a₃，a₄，a₅)时，首先建立最小二乘标准 $\underset{{\mathrm{\&theta;}}_k}{\min }\underset{x}{\mathrm{\&Sigma;}}{[Y_k\left(m_k^g(x;{\mathrm{\&theta;}}_k)\right)-Y_{\mathrm{ref}}\left(x\right)]}^2,$ 其中 $m_k^g(x;{\mathrm{\&theta;}}_k)=\left[m_{k,u}^g(x;{\mathrm{\&theta;}}_k)m_{k,v}^g(x;{\mathrm{\&theta;}}_k)\right],$ $m_{k,u}^g(x;{\mathrm{\&theta;}}_k)=a_0+a_1{x}_u+a_2{x}_v,m_{k,v}^g(x;{\mathrm{\&theta;}}_k)=a_3+a_4{x}_u+a_5{x}_v;$ 然后将θ_k写成θ_k+Δθ，以增量的方式求解θ_k，也即θ_k←θ_k+Δθ。对最小二乘目标泛函进行泰勒展开得到关于Δθ的线性函数，对展开后的函数进行一系列运算操作得到 $\mathrm{\&Delta;\&theta;}=H^{-1}\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T[Y_{\mathrm{ref}}\left(x\right)-Y_k\left(m_k^g(x;{\mathrm{\&theta;}}_k)\right)],$ 其中 $H=\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{\&PartialD;m_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T\underset{x}{\mathrm{\&Sigma;}}[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}],$ $\&dtri;Y_k=\frac{{\&PartialD;Y}_k}{{\&PartialD;m}_k^g(x;{\mathrm{\&theta;}}_k)},$ $\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}=\left(\begin{array}{cccccc}{x}_u& 0& x_v& 0& 1& 0\\ 0& x_u& 0& x_v& 0& 1\end{array}\right).$ 当运动参数增量Δθ小于某个阈值∈的时候(0≤E≤0.01)，也即||Δθ||≥∈，参数θ_k收敛，停止迭代，得到全局运动参数θ_k。(a₀，a₁，a₂，a₃，a₄，a₅)，最后可以得到全局参数运动否则更新θ_k，也即θ_k←θ_k+Δθ，并重新计算Δθ。

求解局部光流

时，首先利用图像灰度恒常性假设和图像梯度恒常性假设得到关于局部光流

的数据置信度能量函数 $E_{\mathrm{Data}}(u,v)=\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&Psi;}({|Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-Y_k\left(\stackrel{\&OverBar;}{x}\right)|}^2+|\&dtri;Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-\&dtri;Y_k{|}^2)d\stackrel{\&OverBar;}{x},$ 其中， $\stackrel{\&OverBar;}{x}={(x,t)}^T={(x_u,x_v,1)}^T,v={(m_k^g{\left(x\right)}^T,0)}^T,\mathrm{\&omega;}={(m_k^l{\left(x\right)}^T,1)}^T={(u,v,1)}^T,$ 函数

0<τ≤0.01；然后根据图像分段平滑性假设得到平滑惩罚函数

最终得到整个能量函数为E(u，v)＝E_Data+αE_Smooth，其中α>0为正则化参数。采用非线性数值求解方法对整个能量函数求最优值便得到局部光流

的解。

最终得到非参考帧图像相对于参考帧图像的运动场m_k(x)＝m_k ^g(x)+m_k ^l(x)，并以此构建运动变换矩阵M_k。首先将运动场m_k(x)进行r倍(最终重建的高清图像的放大倍数)线性插值，得到插值_后的运动场m_k(x)。然后计算相对位移Δx_k＝m_k(x)-x＝(Δx_k，u-Δx_k，v)，以及d_k＝Δx_k，u-floor(Δx_k，u)，e_k＝Δx_k，v-floor(Δx_k，v)，其中操作符foor(.)表示取小于或者等于指定值的最大整数。最后计算运动变换矩阵M_k中的每个元素的值：M_k(j*m+i，floor(Δx_k，u)+x_u+ceil(Δx_k，v+x_v+1)*m)＝d_k*1-e_k)，M_K(j*m+i，ceil(Δx_k,u)+z_u+ceil(Δx_k，v+x_v+1)*m)＝d_k*e_k，M_k(j*m+i，floor(Δx_k，u)+x_u+floor(Δx_k，v+x_v+1)*m)＝(1-d_k)*(1-e_k)，M_k(j*m+i，ceil(Δx_k，u)+x_u+floor(Δx_k，v+x_v+1)*m+x)＝(1-d_k)*e_k，其中操作符ceil(.)表示取大于或者等于指定值的最小整数，对于矩阵M_k的第j*m⁺i行中的除这四列以外的其他元素值为零。

步骤S2、S3完成了图像配准和运动矩阵M_k的构建。

S4、计算参考帧图像Y_ref的r倍双三次插值图像

、非局部先验平滑参数h_i，j和相似图像的欧式阈值。由于理想图像X事先未知，非局部权重矩阵便根据低分辨率的参考帧图像Y_ref的插值图像

估计得到。因此，平滑参数h_i，j是一个关于参考帧图像Y_ref的噪声、图像数据本身、下采样倍数的函数。

其中std(N_i，j)是以像素点(i，j)为中心，搜索区域为N_i，j的标准差，为大于零的常数(1<β<5)，σ²为参考帧图像的噪声方差，r为图像放大倍数。该参考帧图像Y_ref的噪声方差σ²估计为：

其中n为图像的像素点总数，

N(x_i)为像素点x_i的四邻域。最后采用欧式阈值进行不相似像素去除，欧式阈值(即相似度阈值)采用的最大加权欧式距离为4σ²，也即： ${\left|\right|R_{k,l}\stackrel{.}{X}-R_{i,j}\stackrel{.}{X}\left|\right|}_{2,a}^2\&le;4{\mathrm{\&sigma;}}^2\&ap;4{\widehat{\mathrm{\&sigma;}}}^2.$

S5、计算每个像素点(i，j)(其中0≤i<rm，0≤j<rn)的与其他像素点(s，t)(其中0≤s<rm，0≤t<rn)的相似度权重w_NLM[i，j；s，t]，构建关于以插值图像

作为初始高分辨图像X的非局部权重矩阵S包含如下步骤：计算每个像素点(i，j)与其邻域N(i，j)内的像素(s，t)的相似度

其中R_i，j是以(i，j)为中心像素点的图像块提取算子，一般为q×q的方块(q＝5，7，9…)，

表示两个图像块的加权欧式距离，其中a>0是高斯核函数的标准差，h_i，j为滤波器平滑参数，依赖于图像的噪声大小以及图像本身，f为依赖于两个中心像素点的几何距离的常函数(单调不增函数)；归一化像素(i，j)与邻域N(i，j)内的所有像素相似度值， $\mathrm{\&omega;}[i.j;s,t]=\frac{{\mathrm{\&omega;}}_{\mathrm{NLM}}[i,j;s,t]}{{\mathrm{\&Sigma;}}_{(s,t)\&Element;N(i,j)}{\mathrm{\&omega;}}_{\mathrm{NLM}}[i,j;s,t]};$ 最终构建非局部权重矩阵 $S_{j*m+i,s*m+t}=\left\{\begin{array}{cc}\mathrm{\&omega;}[i,j;s,t]& (s,t)\&Element;N(i,j)\\ 0& \mathrm{otherwise}\end{array}\right..$

步骤S4、S5完成了非局部先验参数的自适应求解和非局部权重矩阵的S构建。

S6、利用运动变换矩阵M_k和非局部权重矩阵S求解目标泛函 $X=\arg \min [\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}{\left|\right|Y_k-{\mathrm{DB}}_k{M}_{k}X\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|(I-S)\left|\right|}_{\mathrm{\&rho;}}^{\mathrm{\&rho;}}],$ 其中ρ＞0。采用共轭梯度迭代法最小化目标泛函；

S7、判断是否符合迭代终止条件：

其中n为迭代次数，若不符合，则返回步骤S6，若符合，则表明解收敛，迭代终止，则进入步骤S8；

S8、得到重建的高分辨率估计图像。

步骤S6至步骤S8完成了求解目标泛函。

如图2所示，为本发明的求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)的实施例的流程图，包括如下步骤S31A至步骤S32A：

S31A、采用a₀，a₁，a₂，a₃，a₄，a₅六个参数的仿射变换作为全局参数运动模型： $m_{k,u}^g(x;{\mathrm{\&theta;}}_k)=a_0+a_1{x}_u+a_2{x}_v,m_{k,v}^g(x;{\mathrm{\&theta;}}_k)=a_3+a_4{x}_u+a_5{x}_v;$

S32A、建立最小二乘标准

求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)，包括如下步骤S321A至步骤S324A：

S321A、由于最小二乘目标泛函关于θ_k是非线性的需将θ_k写成θ_k+Δθ的形式，采用增量迭代方式求解θ_k，也即θ_k←θ_k+Δθ。并对最小二乘目标泛函进行泰勒展开得到关于Δθ的线性函数；S322A、对展开后的函数进行运算操作得到 $\mathrm{\&Delta;\&theta;}=H^{-1}\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T[Y_{\mathrm{ref}}\left(x\right)-Y_k\left(m_k^g(x;{\mathrm{\&theta;}}_k)\right)],$ 其中 $H=\underset{x}{\mathrm{\&Sigma;}}{[\&dtri;Y_k\frac{\&PartialD;m_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}]}^T\underset{x}{\mathrm{\&Sigma;}}[\&dtri;Y_k\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}],$ $\&dtri;Y_k=\frac{{\&PartialD;Y}_k}{{\&PartialD;m}_k^g(x;{\mathrm{\&theta;}}_k)},$ $\frac{{\&PartialD;m}_k^g}{{\&PartialD;\mathrm{\&theta;}}_k}=\left(\begin{array}{cccccc}{x}_u& 0& x_v& 0& 1& 0\\ 0& x_u& 0& x_v& 0& 1\end{array}\right);$

S323A、判断是否符合||Δθ|≤∈，阈值0≤∈≤0.01，若不符合，则θ_k←θ_k+Δθ，并返回步骤S322A，若符合，则表明参数θ收敛，停止迭代，进入步骤S324A；S324A、得到全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)。

如图3所示，为本发明的求解局部光流

的实施例的流程图。在得到全局运动参数θ＝(a₀，a₁，a₂，a₃，a₄，a₅)，也即全局运动场后，如果仿射运动已经可以很好地近似物体的运动模型，此时便可以省略对运动场进行补偿的局部光流场的求解。但如果全局仿射模型未能精确地近似整个运动场，便需要求解局部光流场。然而，已经求解的全局参数运动场已经在很大程度上接近整体运动场，有助于接下来的局部光流场的求解，使得求解的局部光流场更加精确。这就体现了本发明提出的复合运动模型的良好的可伸缩性、高精度的运动估计、以及低的运算量。为了求得局部光流

的解，本发明利用Thomas Brox提出的光流求解法。ThomasBrox光流求解法主要基于三个假设：图像灰度恒常性假设、图像梯度恒常性假设、图像分段平滑性假设。

求解局部光流

包括如下步骤S31B至S34B：

S31B、利用图像灰度恒常性假设和图像梯度恒常性假设得到关于局部光流

的数据置信度能量函数 $E_{\mathrm{Data}}(u,v)=\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&Psi;}({|Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-Y_k\left(\stackrel{\&OverBar;}{x}\right)|}^2+|\&dtri;Y_k(\stackrel{\&OverBar;}{x}+v+\mathrm{\&omega;})-\&dtri;Y_k{|}^2)d\stackrel{\&OverBar;}{x},$ 其中， $\stackrel{\&OverBar;}{x}={(x,t)}^T={(x_u,x_v,1)}^T,v={(m_k^g{\left(x\right)}^T,0)}^T,\mathrm{\&omega;}={(m_k^l{\left(x\right)}^T,1)}^T={(u,v,1)}^T,$ 函数 $\mathrm{\&Psi;}\left(s^2\right)=\sqrt{s^2+{\mathrm{\&tau;}}^2},$ 0<τ≤0.01：

S32B、根据光流的图像分段平滑性假设得到平滑惩罚函数 $E_{\mathrm{Smooth}}(u,v)=\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&Psi;}\left({|\&dtri;u|}^2{+|\&dtri;v|}^2\right)d\stackrel{\&OverBar;}{x};$

S33B、得到整个能量函数为E(u，v)＝E_Data+αE_Smooth，其中α>0为正则化参数；

S34B、采用非线性数值求解方法对整个能量函数求最优值便得到局部光流

的解。

以上介绍了一种基于复合运动和自适应非局部先验的超分辨率重建方法。本发明并不限定于以上实施例，任何未脱离本发明技术方案，即仅仅对其进行本领域普通技术人员所知悉的改进或变更，均属于本发明的保护范围之内。

Claims (10)

1.一种基于复合运动和自适应非局部先验的超分辨率重建方法，高分辨率图像X经过降质过程获取p幅低分辨率观测图像Y_k(k＝1，…，p)，每幅观测图像的大小为m×n，其特征在于，所述利用p幅低分辨率观测图像Y_k(k＝1，…，p)重建高分辨率估计图像的超分辨率重建方法包括如下步骤： 

(1)在p幅低分辨率观测图像中选取参考帧图像Y_ref(1≤ref≤p)和非参考帧图像Y_k(k＝1，…，ref-1，ref+1，…，p)，针对参考帧图像和非参考帧图像之间的亚像素运动采用全局参数运动和局部光流的复合运动模型，参考帧图像与非参考帧图像之间的关系表示为： 

其中m_k表示二维运动场，

表示全局参数运动，

为局部光流运动，θ_k为全局运动参数，

表示用非参考帧图像预测的参考帧图像，ε_k(x)表示残差图像； 

(2)求解全局运动参数θ_k＝(a₀，₁，a₂，a₃，a₄，a₅)和局部光流采用全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)和局部光流

的方法进行图像配准，得到非参考帧图像相对于参考帧图像的运动场m_k(x)，利用m_k(x)构建出运动变换矩阵M_k； 

(3)计算参考帧图像Y_ref的r倍插值图像X、非局部先验参数h_i，j(0≤i<rm，0≤j<rn)和相似图像的欧式阈值； 

(4)利用非局部先验参数、欧式阈值，且以插值图像

作为高分辨图像X的初始图像，计算出插值图像中每个像素点(i，j)与其他像素点(s，t)的相似度权重w_NLM[i，j；s，t]，其中0≤i<rm，0≤j<rn，0≤s<rm，0≤t<rn，利用相似度权重w_NLM构建非局部权重矩阵S； 

(5)利用运动变换矩阵M_k和非局部权重矩阵S求解目标泛函 

其中B_k为观测图像Y_k对应的降晰函数，M_k为观测图像Y_k相对于参考帧图像的亚像素运动，非局部权重矩阵S是一个自适应高分辨率图像X的非局部均值滤波器，且ρ>0，采用共轭梯度迭代法最小化目标泛函，得到重建的高分辨率估计图像。 

2.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(1)针对参考帧图像和非参考帧图像之间的亚像素运动采用全局参数运动和局部光流的复合运动模型，其中运动向量为m_k(x)＝[m_k，u(x)m_k，v(x)]，x＝[x_ux_v]，二维运动场m_k表示为

其中

表示全局参数运动，

为局部光流运动，θ_k为全局运动参数。 

3.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(2)求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)包括如下步骤： 

(21A)采用a₀，a₁，a₂，a₃，a₄，a₅六个参数的仿射变换作为全局参数运动模型：

(22A)建立最小二乘标准

求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)。 

4.根据权利要求3所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(22A)建立最小二乘标准 

求解全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)包括 

如下步骤： 

(221A)将θ_k写成θ_k+Δθ； 

(222A)对最小二乘目标泛函进行泰勒展开得到关于Δθ的线性函数； 

(223A)对展开后的函数进行一系列运算操作得到 其中 

(224A)判断是否符合||Δθ||≤e，阈值0≤∈≤0.01，若不符合，则θ_k←θ_k+Δθ，并返回步骤(223A)；若符合，则表明参数θ_k收敛，停止迭代，得到全局运动参数θ_k＝(a₀，a₁，a₂，a₃，a₄，a₅)，最后可以得到全局参数运动。

5.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(2)求解局部光流

包括如下步骤： 

(21B)利用图像灰度恒常性假设和图像梯度恒常性假设得到关于局部光流

的数据置信度能量函数 

其中， 

函数 

其中0<τ≤0.01； 

(22B)根据图像分段平滑性假设得到平滑惩罚函数 

(23B)得到整个能量函数为E(u，v)＝E_Data+αE_Smooth，其中α>0为正则化参数； 

(24B)采用非线性数值求解方法对整个能量函数E(u，v)＝E_Data+αE_Smooth求最优值便得到局部光流

的解。 

6.根据权利要求5所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(2)利用m_k(x)构建出运动变换矩阵M_k包含如下步骤： 

(21C)计算非参考帧图像相对于参考帧图像的运动场m_k(x)＝m_k ^g(x)+m_k ^l(x)，并将运动场m_k(x)进行r倍线性插值，得到插值后的运动场m_k(x)，所述r为最终重建的高清图像的放大倍数； 

(22C)计算相对位移Δx_k＝m_k(x)-x＝(Δx_k，u-Δx_k，v)，以及d_k＝Δx_k，u-floor(Δx_k，u)，e_k＝Δx_k，v-floor(Δx_k，v)，其中操作符floor(.)表示取小于或者等于指定值的最大整数； 

(23C)计算运动变换矩阵M_k中的每个元素的值：M_k(j*m+i，floor(Δx_k，u)+x_u+ceil(Δx_k，v+x_v+1)*m)＝d_k*(1-e_k)，M_k(j*m+i，ceil(Δx_k，u)+x_u+ceil(Δx_k，v+x_v+1)*m)＝d_k*e_k，M_k(j*m+i，floor(Δx_k，u)+x_u+floor(Δx_k，v+x_v+1)*m)＝(1-d_k)*(1-e_k)，M_k(j*m+i，ceil(Δx_k，u)+x_u+floor(Δx_k，v+x_v+1)*m+x)＝(1-d_k)*e_k，其中操作符ceil(.)表示取大于或者等于指定值的最小整数，对于矩阵M_k的第j*m+i行中的除这四列以外的其他元素值为零。 

7.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨 率重建方法，其特征在于：所述步骤(3)计算非局部先验参数h_i，j， 

其中std(N_i，j)是以像素点(i，j)为中心，搜索区域为N_i，j的标准差，β为大于零的常数(1<β<5)，σ²为参考帧图像Y_ref的噪声方差，r为图像放大倍数，所述参考帧图像Y_ref的噪声方差σ²估计为：

其中n为图像的像素点总数， 

N(x_i)为像素点x_i的四邻域。 

8.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(3)中相似图像的欧式阈值采用的最大加权欧式距离为4σ²，也即：

9.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(4)利用非局部先验参数、欧式阈值，且以插值图像

作为高分辨图像X的初始图像，计算出插值图像

中每个像素点(i，j)与其他像素点(i，j)的相似度权重w_NLM[i，j；s，t]，利用相似度权重w_NlM构建非局部权重矩阵S包含如下步骤： 

(41)计算每个以坐标(i，j)为中心的像素块与其邻域N(i，j)内以坐标(i，j)为中心的像素块的相似度 

，其中R_i，j是以(i，j)为中心像素点的图像块提取算子，表示两个图像块的加权欧式距离，其中a>0是高斯核函数的标准差，h_i，j为滤波器平滑参数，依赖于图像的噪声大小以及图像本身，f为依赖 于两个中心像素点的几何距离的常函数； 

(42)归一化像素(i，j)与邻域N(i，j)内的所有像素相似度值， 

(43)构建非局部权重矩阵

10.根据权利要求1所述的基于复合运动和自适应非局部先验的超分辨率重建方法，其特征在于：所述步骤(5)中当

时，迭代终止，其中n为迭代次数。 

Images