CN1555097A - 面阵电荷耦合器件超分辨率成象技术中的快速算法 - Google Patents

Abstract

一种面阵电荷耦合器件超分辨率成象技术中的快速算法，本发明涉及利用对同一目标的有相互位移的多帧图象进行处理得到一帧分辨率更高图象的面阵电荷耦合器件超分辨率成象技术领域。其特征在于它建立了像素级的邻域双线性插值图象降质模型，采用四维非稀疏矩阵存储模型系数，避免了对大型稀疏矩阵的存储和处理；对亚像素位移量进行量化处理，大大降低了计算过程对存储空间的需求；给出改进Cimmino迭代算法，加快了迭代收敛速度，减小计算量，适合PC机或DSP等嵌入式系统的实现，为高分辨率面阵CCD的超分辨率成象问题应用于工程实际，研制超分辨率成象系统解决了其中的两个关键技术问题。

Description

面阵电荷耦合器件超分辨率成象技术中的快速算法

技术领域：

本发明涉及图象成象技术领域，具体涉及利用对同一目标的有相互位移的多帧图象进行处理得到一帧分辨率更高图象的面阵电荷耦合器件超分辨率成象技术的算法。

背景技术：

面阵电荷耦合器件(简称：CCD)超分辨率成象技术研究具有重要的意义，由于工艺限制，目前以电荷耦合器件为接收器的成象系统分辨率与许多实际需求总有一定的差距，超分辨率成象系统可利用已有电荷耦合器件获取高于自身分辨率的图象，可达到获取高分辨率图象的目的。特别是对红外成象系统，由于灵敏度问题，红外电荷耦合器件的像素一般较少，获取的图象分辨率较低，而在夜视军事侦察等领域，提高成象系统分辨率是非常必要的。由于高分辨率电荷耦合器件尤其是红外高像素电荷耦合器件常用于航天、国防等领域，国外对我国进行出口限制，因此研究利用现有电荷耦合器件获取更高分辨率图象的面阵电荷耦合器件超分辨率成象技术意义重大。

目前，提高电荷耦合器件成象分辨率的方法主要有以下四类：

1.改进电荷耦合器件制作工艺，增加电荷耦合器件像素数。随着技术水平的提高，电荷耦合器件器件的像素数逐年增加，但总是与人们的最高需求有差距；

2.多块电荷耦合器件几何拼接方法。可有效增加系统分辨率，但拼接工艺要求很高，系统结构十分复杂，而且增加了成象镜头的设计难度，视场角扩大，成为一个超大视场角光学系统，其镜头成象质量难以保证。

3.对电荷耦合器件像元进行形状和排列方式的重新设计，采用蜂窝状排列和多边形形状，使得单位面积上能够放置更多的与传统电荷耦合器件感光面积相同的感光单元，该技术工艺上还有难度，而且分辨率提高幅度也有限。

4.基于有相互位移的图象序列的超分辨率图象重建。该方法既克服了直接增加感光单元提高电荷耦合器件分辨率所面对的工艺难题，又能够充分利用多幅图象的不同信息达到获取更高分辨率图象。

基于序列图象的超分辨成象技术的主要难点集中在两个方面：①有相对亚像素级微位移的序列图象的获取方法；②超分辨图象的重建算法，尤其是可应用于工程实际的快速重建算法。

其中有相对亚像素级微位移的序列图象的获取方法主要包括：利用棱镜进行多路分光法，利用多块双折射晶体和偏振片的双折射法、压电陶瓷微位移法、微变焦法等。超分辨图象的重建算法也有多种，如MAP算法、POCS算法、IBP算法、反演解析法等。上述算法各有优缺点，但一个共同的问题是算法的计算量非常大，计算过程对存储空间的要求也非常大，对普通PC机或DSP处理系统等而言，处理起来非常困难，因此要应用于工程实际还有一定的距离。

发明内容：

本发明所要解决的技术问题是超分辨率重建算法应用于工程实际中，研制超分辨率成象系统的两个关键问题：重建算法计算量太大和占用存储空间太大的问题。

本发明是通过下面的技术方案解决上述技术问题的。其特征在于它包括下述步骤：

1)建立像素级双线性插值图象降质模型：设p帧观察图象分辨率为m×n，分辨率提高因子为L，建立包含空间移不变邻域模糊、加性噪声、欠采样等因素的像素级双线性插值图象降质模型，观察图象中的每一个像素值为该像素在高分辨率图象中对应点的一个领域内的像素值的加权平均，根据图象间的配准关系，得到观察图象中任一像素C，其坐标为(k₁，l₁′)，经配准后在理想图象X中对应点C′，其坐标为(k₂′，l₂′)，其中k₂′＝k₂+δ_k，l₂′＝l₂+δ_l，0≤δ_k，δ_l＜1，k₂，l₂为整数，再由双线性插值公式

$Y_{k_1,l_1}=(1-{\mathrm{\δ}}_k-{\mathrm{\δ}}_l+{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2,l_2}^{\′}+({\mathrm{\δ}}_k-{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2+1,l_2}^{\′}$

$+({\mathrm{\δ}}_l-{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2,l_2+1}^{\′}+\left({\mathrm{\δ}}_k{\mathrm{\δ}}_l\right)X_{k_2+1,l_2+1}^{\′}$

其中

$X_{k,l}^{\′}=\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}{h}_{i,j}{X}_{k+i,l+j}$ k＝k₂，k₂+1；l＝l₂，l₂+1

合并上两式得

$Y_{k_1,l_1}=\underset{i=-1}{\overset{2}{\mathrm{\Σ}}}\underset{j=-1}{\overset{2}{\mathrm{\Σ}}}{h}_{i,j}^{\′}{X}_{k_2+i,l_2+j}=h^{\′\′}\·\×X_{{\mathrm{\Ω}}_{k_2,l_2}}$

2)对线性方程组：

                               Ax＝b

其中A＝(a₁，a₂，Λ，a_m)^T∈K^mn；a₁，a₂，Λ，a_m∈Kⁿ；b＝(b₁，b₂，Λ，b_m)^T∈K^m；x＝(x₁，x₂，Λ，x_n)^T∈Kⁿ，K为数域，包含实数域R和复数域C

采用迭代式：

${x_j}^{k+1}=x_j^k+\frac{{\mathrm{\λ}}_k}{q_j}\underset{i=1}{\overset{q_j}{\mathrm{\Σ}}}\frac{b_{s_i}-a_{s_i}^T{x}^k}{{\left|\right|a_{s_i}\left|\right|}^2}{a}_{s_i,j}$ j＝1∶n；0＜λ_k＜2

3)将亚像素位移量δ_k和δ_l量化到小数点后一位或两位，则δ_k和δ_l的取值分别有10种或100种取值，组合起来有100种或10000种取值，对每一种取值计算方程组的系数，得到一张系数表，将其提前计算好存储起来，迭代求方程组时计算出量化后的δ_k和δ_l，然后通过查表得到方程组系数；

4)建立与X维数相同的权值矩阵W和修正矩阵ΔX，其中W中每个元素为对应修正矩阵ΔX中像素点在每次迭代中已被修改的次数，设g为迭代次数，算法步骤为：

(1)初始值任意，设为0：g＝0，X^g＝0；

(2)权值矩阵和修正矩阵置0矩阵：W＝Φ，ΔX＝Φ；

(3)对观察图象的每个像素计算在理想高分辨率图象中的对应点坐标，从而得到k₂、δ_k、l₂、δ_l，再算出h″、Ω_k2，l2，并且

${\mathrm{\ΔX}}_{i,j}=\frac{W_{i,j}}{W_{i,j}+1}\mathrm{\Δ}{X}_{i,j}+\frac{1}{W_{i,j}+1}[Y_{k_1,l_1}-h^{\′\′}\·\×{X^g}_{{\mathrm{\Ω}}_{k_2,l_2}}]h_{i,j}^{\′}/{\left|\right|h^{\′\′}(:)\left|\right|}_2^2---i,j\∈{\mathrm{\Ω}}_{k_2,l_2}$

             W_i，j＝W_i，j+1 $i,j{\∈\mathrm{\Ω}}_{k_2,l_2}$

(4)X^g+1＝X^g+λ_gΔX，g＝g+1

上述算法在并行处理环境中把观察图象的像素分为P组，P为并行计算的路数，对每一组都需定义权值矩阵和修正矩阵，然后对各组分别处理，得到相应的Wⁱ和ΔX^l(i＝1∶P)

$X_{k,l}^{g+1}=X_{k,l}^g+{\mathrm{\λ}}_g\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\left((W_{k,l}^l/\underset{j=1}{\overset{P}{\mathrm{\Σ}}}{W}_{k,l}^j)\mathrm{\Δ}{X}_{k,l}^i\right)$ k＝1∶Lm；l＝1∶Ln

本发明具有下列技术效果：

1、建立了像素级的邻域双线性插值图象降质模型，采用四维非稀疏矩阵存储模型系数，避免了对大型稀疏矩阵的处理，可降低一半的存储空间，大大加快寻址时间，加快处理速度；

2、对亚像素位移量进行量化处理，并将系数矩阵提前计算好，计算过程中通过实时查表的方式得到系数矩阵，大大降低了计算过程对存储空间字节数的需求，由G级降低到K级。

3、给出了改进Cimmino迭代算法，加快了迭代收敛速度，减小计算量，而且每次只处理其中一行，计算过程本身对内存要求也比较低；

通过上述改进，使超分辨率重建算法适合PC机或DSP等嵌入式系统的实现，为高分辨率面阵电荷耦合器件的超分辨率成象问题应用于工程实际，研制超分辨率成象相机解决了其中的两个关键技术问题。

附图说明：

图1为观察图象与理想高分辨率图象关系示意图；

图2为本发明像素级成象模型示意图；

图3为本发明非相容系统Kaczmarz算法迭代过程示意图；

图4为本发明非相容系统Cimmino算法迭代过程示意图；

图5为本发明与现有技术Cimmino算法收敛性比较图；

图6(a)为本发明算法实验结果一之低分辨率图象序列图；

图6(b)为本发明算法实验结果一之双线性插值结果图；

图6(c)为本发明算法实验结果一之双三次方插值结果图；

图6(d)为本发明算法实验结果一之经2次迭代处理后效果图；

图7(a)为本发明算法实验结果二之原始分辨率测试图象；

图7(b)为本发明算法实验结果二之低分辨率观察图象序列；

图7(c)为本发明算法实验结果二之双线性插值结果图；

图7(d)为本发明算法实验结果二之双三次方插值结果图；

图7(e)为本发明算法处理后效果图。

具体实施方式：

本发明是在研究现有技术的基础上进行改进的，它提出了像素级的邻域双线性插值图象降质模型，提出四维系数矩阵存储结构，通过量化亚像素位移量，建立系数表，将原模型系数的占用空间由G级降低到K级。同时提出改进Cimmino行处理迭代算法，给出实现空间域超分辨图象重建的快速算法，不仅使迭代收敛速度快，而且算法本身对处理系统内存的要求也非常低，更加速了算法的计算速度。对多CPU组成的并行系统，还给出了并行算法。

下面对算法步骤做出说明。

一、建立像素级双线性插值图象降质模型：

为便于对比和说明的方便，先给出一般的空域图象降质数学模型：

设p帧观察图象分辨率为m×n，分辨率提高因子为L，设低分辨率观察图象为一帧高分辨率理想图象经模糊、形变、加性噪声、降采样后的结果，并建立图象离散成象模型为：

                            Y＝H X+ E               (1)

其中 X为[L²mn×1]矩阵，其值为原理想图象[Lm×Ln]个像素点的逐行排列； Y为多帧观察图象逐行排列后形成的列向量，维数为pmn×1， E _k为系统噪声，维数与 Y相同，H为包含图象的几何形变、模糊、降采样等因素的系数矩阵，维数为pmn×L²mn。

当噪声信息未知时，(1)式一般为大规模、超定、稀疏线性方程组，而且由于量化误差、噪声等因素，方程组往往是不相容的，可以求方程组的最小二乘解作为目标高分辨率图象。

分析(1)式中系数矩阵的维数，设观察图象分辨率为1300×1300，帧数p＝6，当L＝2时H的维数为10140000×6760000，占用的存储空间为5.4×10¹⁴字节。即使考虑到H的稀疏性，用稀疏矩阵的结构处理该矩阵，设退化算子为3×3线性空间移不变邻域模糊算子，稀疏矩阵占用的空间也超过2.6×10⁹＝2.6G字节(稀疏矩阵每元素占16个字节)，而且由于稀疏矩阵的存储结构特点，计算中寻址非常慢，处理起来比非稀疏矩阵要慢得多，实践证明对稀疏矩阵中元素的操作比对普通矩阵的操作要慢1个数量级。因此在一般PC机或嵌入式系统上存储或计算都比较困难。

为此建立包含空间移不变邻域模糊、加性噪声、欠采样等因素，像素级的双线性插值图象降质模型。如图1所示，观察图象中的每一个像素值为该像素在高分辨率图象中对应点的一个领域内的像素值的加权平均。

设退化模型为线性空间移不变邻域模糊算子，设模糊算子为3×3矩阵h。

图象之间的几何关系可建立为图象配准的仿射变换六参数模型：

$\left[\begin{array}{c}{x}_i^{\′}\\ y_i^{\′}\end{array}\right]=\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]+\left[\begin{array}{c}{a}_3\\ b_3\end{array}\right]----\left(2\right)$

此模型中的参数a₁，a₂，a₃，b₁，b₂，b₃为两图象的配准参数。采用基于光流计算的金字塔形分层迭代算法计算出每一帧图象与参考图象之间的仿射变换参数。

如图2所示，根据仿射变换关系，可以得到观察图象中任一像素C(设坐标为(k₁，l₁′))经配准后在理想图象X中对应点C′，其坐标为(k₂′，l₂′)，其中k′₂＝k₂+δ_k，l₂′＝l₂+δ_l，0≤δ_k，δ_l＜1，k₂，l₂为整数。图中将图1中高分辨率图象的每个像素抽象为一个点，点旁的序号为其坐标。由双线性插值公式

$Y_{k_1,l_1}=(1-{\mathrm{\δ}}_k-{\mathrm{\δ}}_l+{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2,l_2}^{\′}+({\mathrm{\δ}}_k-{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2+1,l_2}^{\′}----\left(3\right)$

$+({\mathrm{\δ}}_l-{\mathrm{\δ}}_k{\mathrm{\δ}}_l)X_{k_2,l_2+1}^{\′}+\left({\mathrm{\δ}}_k{\mathrm{\δ}}_l\right)X_{k_2+1,l_2+1}^{\′}$

其中

$X_{k,l}^{\′}=\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}{h}_{i,j}{X}_{k+i,l+j}$ k＝k₂，k₂+1；l＝l₂，l₂+1       (4)

合并上两式得

$Y_{k_1,l_1}=\underset{i=-1}{\overset{2}{\mathrm{\Σ}}}\underset{j=-1}{\overset{2}{\mathrm{\Σ}}}{h}_{i,j}^{\′}{X}_{k_2+i,l_2+j}=h^{\′\′}\·\×X_{{\mathrm{\Ω}}_{k_2,l_2}}----\left(5\right)$

其中h_i，j′在h确定的情况下仅与δ_k，δ_l有关；h″和

为4×4矩阵，Ω_k2，l2为式中X的计算窗口，“·×”符号表示两个同型矩阵对应元素之积的和。

对观察图象中的每一个像素，都对应这样的一个4×4矩阵，这里采用四维矩阵将所有系数矩阵存储起来，矩阵维数为：pm×n×4×4，矩阵的前两维分别对应观察图象中每个像素点的位置。其占用空间为原模型稀疏存储方式的一半。

该模型将原来的大型稀疏方程组分解为一个个低维数的线性方程，充分利用了观察图象和高分辨率理想图象之间的邻域结构关系，建立四维系数矩阵，避免了对大型稀疏矩阵的存储和处理，降低了存储空间，加快了计算速度。

下面给出求解大型稀疏非相容线性方程组的行处理迭代算法。

二、对Cimmino行处理迭代算法进行改进

设线性方程组：

                                  Ax＝b            (6)

其中 $A={(a_1,a_2,\mathrm{\Λ},a_m)}^T\∈K^{\mathrm{mn};}$

a₁，a₂，Λ，a_m∈Kⁿ；

b＝(b₁，b₂，Λ，b_m)^T∈K^m；

x＝(x₁，x₂，Λ，x_n)^T∈Kⁿ：

K为数域，包含实数域R和复数域C。

在实际物理系统中，A往往为系统的模型，b为观测数据，x为被测量。由于方程的病态性、观察噪声的存在、巨大的维数、A无规律性结构等因素，(6)式往往为大型非相容稀疏线性方程组，用直接法和常规迭代法求解比较困难，而行处理方法则适合大规模(维数大于10⁵)、稀疏(非0值小于1％)、系数矩阵不对称或其它规律性结构等情况下的方程组求解问题。

目前常用的行处理方法主要有Kaczmarz算法和Cimmino算法。Kaczmarz算法是由任意初始值开始，依次向各个超平面投影，逐次得到下一个迭代值。算法迭代式为：

$x^{k+1}=x^k+{\mathrm{\&lambda;}}_k\frac{b_i-a_i^H{x}_k}{{\left|\right|a_i\left|\right|}_2^2}{a}_i---i=k\mathrm{mod}m+1;0<{\mathrm{\&lambda;}}_k<2;----\left(7\right)$

其中H符号表示共轭转置，‖·‖₂为向量的2范数，λ_k为松弛因子，初始值x⁰可任取。

算法收敛速度快，适合于相容方程组的求解，对不相容方程组，一般情况下不收敛，在解附近振荡，如图3所示，图中1、2、3分别为3个超平面。只有当松弛因子为足够小，才能保证收敛，但这时由于每步迭代修正量太小，收敛速度大大降低。

Cimmino算法是由当前迭代值向各个超平面投影，然后计算所有投影矢量的平均值作为迭代修正量，适合于不相容方程组，当x⁰∈R(A^H)时，收敛到系统的极小范数最小二乘解。算法迭代式为：

$x_j^{k+1}=x_j^k+\frac{{\mathrm{\&lambda;}}_k}{m}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{b_i-a_i^H{x}^k}{{\left|\right|a_i\left|\right|}^2}{a}_{i,j}---j=1:n;0<{\mathrm{\&lambda;}}_k<2----\left(8\right)$

迭代过程示意图如图4所示，其中λ_k＝1，设当前迭代值为x^k，分别向各超平面作正交投影，设垂足分别为A₁，A₂，A₃，将三矢量

相加得到

则取 $x^{k+1}=x^k+\stackrel{\&RightArrow;}{x^k{A}^{\&prime;}}/3,$ 依此类推。

分析式(8)，设λ_k＝1，每次迭代对x_j的修正量为当前点向个各超平面的投影矢量之和再除以系统行向量的个数m。当a_i，j＝0时，即当前超平面a_i与x_j无关，这时投影为空，不应计入投影矢量之和的计算，当Cimmino算法未考虑这个因素。对稀疏性方程组，由于A的大部分元素为0，实际有效投影数量远小于取平均时的基数m，因此每次迭代修正量远小于实际有效投影矢量的平均值，相当于λ_k被除以一个很大的数，故算法收敛速度很慢。

考虑上述因素，提出改进Cimmino迭代式为：

$x_j^{k+1}=x_j^k+\frac{{\mathrm{\&lambda;}}_k}{q_j}\underset{i=1}{\overset{q_j}{\mathrm{\&Sigma;}}}\frac{b_{s_i}-a_{s_i}^T{x}^k}{{\left|\right|a_{s_i}\left|\right|}^2}{a}_{s_i,j}---j=1:n;0<{\mathrm{\&lambda;}}_k<2----\left(9\right)$

其中q_j为A中第j列不为0的元素的个数；1＜s_i＜m，为A的第j列中第i个不为0的元素所在行的行数。

当系统为非稀疏性时，q_j＝m，s_i＝i，(9)式与(8)式相同。对稀疏系统，则q_j＜m，算法收敛速度明显加快。算法的实质与Cimmino算法完全相同，因此收敛性也相同。

分析(9)式，不仅求每步迭代修正量时求平均的基数与实际有效投影个数相同，而且计算投影时只考虑不为0的量，因此对稀疏系统，计算量大大减小。另外，由于x_j(j＝1∶n)的迭代计算之间互相独立，因此算法可以并行实现。

当然算法也具有行处理算法的共性特点：不改变原系数矩阵；不对原始系数矩阵做任何整体操作；每次迭代只需系数矩阵中的一行；每次迭代只需上次迭代值，与历史迭代结果无关。

实验证明，算法对稀疏方程组的求解收敛速度很快，对这里涉及的超分辨率图象重建问题，一般只需3～4次迭代即可。

三、超分辨率图象重建中的快速算法

根据上述算法，结合前述模型的特点，给出新的实现图象超分辨率重建的迭代算法：

首先建立与目标高分辨率图象X维数相同的权值矩阵W和修正矩阵ΔX，其中W中每个元素为对应修正矩阵ΔX中像素点在每次迭代中已被修改的次数。设g为迭代次数，算法分四步：

第1步：初始值任意，设为0：g＝0，X^g＝0；

第2步：权值矩阵和修正矩阵置0矩阵：W＝Φ，ΔX＝Φ；

第3步：对观察图象的每个像素计算在理想高分辨率图象中的对应点坐标，从而得到k₂、δ_k、l₂、δ_l，再算出h″、Ω_k2，l2，并且

${\mathrm{\&Delta;X}}_{i,j}=\frac{W_{i,j}}{W_{i,j}+1}\mathrm{\&Delta;}{X}_{i,j}+\frac{1}{W_{i,j}+1}[Y_{k_1,l_1}-h^{\&prime;\&prime;}\&CenterDot;\&times;{X^g}_{{\mathrm{\&Omega;}}_{k_2,l_2}}]h_{i,j}^{\&prime;}/{\left|\right|h^{\&prime;\&prime;}(:)\left|\right|}_2^2---i,j\&Element;{\mathrm{\&Omega;}}_{k_2,l_2}----\left(10\right)$

$W_{i,j}=W_{i,j}+1---i,j{\&Element;\mathrm{\&Omega;}}_{k_2,l_2}------\left(11\right)$

第4步：X^g+1＝X^g+λ_gΔX，g＝g+1，若‖ΔX‖＜ε则跳出，否则回到第2步，ε为终止迭代条件数。

上述算法为串行算法，在并行处理环境中则要对算法作适当修改，其中在第3步要把观察图象的像素分为P组，P为并行计算的路数，对每一组都需定义权值矩阵和修正矩阵，然后对各组分别处理，得到相应的Wⁱ和ΔXⁱ(i＝1∶P)，再将第4步修正公式改为：

$X_{k,l}^{g+1}=X_{k,l}^g+{\mathrm{\&lambda;}}_g\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\left((W_{k,l}^l/\underset{j=1}{\overset{P}{\mathrm{\&Sigma;}}}{W}_{k,l}^j)\mathrm{\&Delta;}{X}_{k,l}^i\right)---k=1:\mathrm{Lm};l=1:\mathrm{Ln}-----\left(12\right)$

分析(5)式，h″仅于δ_k、δ_l和h有关，其中h是已知的3×3模糊算子矩阵，而亚像素位移量δ_k、δ_l是在0和1之间的小数。考虑到有相对亚像素级微位移的序列图象的获取方法中的微位移精度往往是亚像素级的，即δ_k和δ_l的精度只到小数点后一位，因此对其进行量化，四舍五入到小数点后一位，则δ_k和δ_l的取值分别有10种可能，组合起来有100种可能取值，对每一种可能取值计算方程组的16个系数，可以得到一张系数表，提前计算好该系数表并存储起来，迭代求方程组时只要计算出量化后的δ_k和δ_l，就可查表得到方程组系数。其占用的存储空间为10×10×16×8＝12800字节，对比前述分辨率为1300×1300的例子，不到其稀疏矩阵存储方式的0.0005％。目前面阵CCD的像素数已超过千万像素级，若处理6帧3000×3000的图象，则原有方法的稀疏矩阵存储方式系数矩阵需要13.8G的存储空间，而本方法仍仅需12800字节，不及原来的万分之0.001。即使当微位移精度提高，δ_k和δ_l量化到小数点后2位，也只需要1.28M字节存储空间即可。

节省存储空间的意义在于使原来只是在原理上和理论上研究高分辨率面阵CCD的超分辨率成象技术可以应用于工程实际，使研制超分辨率成象系统成为可能。

需要说明的是，在许多CCD超分辨率成象技术应用中，低分辨率序列图象间的相互位移关系和系统模糊因子是已知的，或者说是先验知识，因此算法的1，2步可以省略。

实验结果：

图6和图7为实际实验效果对比，图6(a)为10帧观察图象序列中的任意4帧图象；图6(b)、图6(c)为对其中两帧图象分别进行双线性插值和双三次方插值的结果；图6(d)为本发明经2次迭代处理后的结果其中λ＝1.5。图7(a)为图象空间分辨率测试图原图，图7(b)为低分辨率观察图象中的任意四帧，图7(c)、图7(d)为将其中两帧经双线性和双三次插值放大的结果，图7(e)为本发明取松弛因子λ_k＝2时第2和第4次迭代之后的结果。图5为实验1中的算法收敛曲线，可见算法的收敛速度比较快，一般经3-4次迭代即可，而Cimmino算法的收敛速度非常慢，这是因为系数矩阵的稀疏因子约为0.018％，远小于1。

本发明中的迭代算法部分也可应用于其它类似大型稀疏非相容线性方程组的求解问题，如医学图象重建、图象复原等。

Claims (1)

1、一种面阵电荷耦合器件超分辨率成象技术中的快速算法，其特征在于它包括下述步骤：

1)建立像素级双线性插值图象降质模型：

设p帧观察图象分辨率为m×n，分辨率提高因子为L，建立包含空间移不变邻域模糊、加性噪声、欠采样等因素的像素级双线性插值图象降质模型，观察图象中的每一个像素值为该像素在高分辨率图象中对应点的一个领域内的像素值的加权平均，根据图象间的配准关系，得到观察图象中任一像素C，其坐标为(k₁，l′₁)，经配准后在理想图象X中对应点C′，其坐标为(k′₂，l′₂)，其中k′₂＝k₂+δ_k，l′₂＝l₂+δ_l，0≤δ_k，δ_l＜1，k₂，l₂为整数，再由双线性插值公式

$Y_{k_1,l_1}=(1-{\mathrm{\&delta;}}_k-{\mathrm{\&delta;}}_l+{\mathrm{\&delta;}}_k{\mathrm{\&delta;}}_l)X_{k_2,l_2}^{\&prime;}+({\mathrm{\&delta;}}_k-{\mathrm{\&delta;}}_k{\mathrm{\&delta;}}_l)X_{k_2+1,l_2}^{\&prime;}$

$+({\mathrm{\&delta;}}_l-{\mathrm{\&delta;}}_k{\mathrm{\&delta;}}_l)X_{k_2,l_2+1}^{\&prime;}+\left({\mathrm{\&delta;}}_k{\mathrm{\&delta;}}_l\right)X_{k_2+1,l_2+1}^{\&prime;}$

其中

$X_{k,l}^{\&prime;}=\underset{i=-1}{\overset{l}{\mathrm{\&Sigma;}}}\underset{j=-1}{\overset{l}{\mathrm{\&Sigma;}}}{h}_{i,j}{X}_{k+i,l+j}$                  k＝k₂，k₂+1；l＝l₂，l₂+1

合并上两式得

$Y_{k_1,l_1}=\underset{i=-1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{j=-1}{\overset{2}{\mathrm{\&Sigma;}}}{h}_{i,j}^{\&prime;}{X}_{k_2+i,l_2+j}=h^{\&prime;\&prime;}\&CenterDot;\&times;X_{{\mathrm{\&Omega;}}_{k_2,l_2}}$

其中h′_i，j在h确定的情况下仅与δ_k，δ_l有关；h″和

为4×4矩阵，Ω_k2，l2为式中X的计算窗口，“·×”符号表示两个同型矩阵对应元素之积的和；

2)对线性方程组：

                                  Ax＝b

其中A＝(a₁，a₂，Λ，a_m)^T∈K^mn；a₁，a₂，Λ，_am∈Kⁿ；b＝(b₁，b₂，Λ，b_m)^T∈K^m；x＝(x₁，x₂，Λ，x_n)^T∈Kⁿ，K为数域，包含实数域R和复数域C

采用迭代式：

${x_j}^{k+1}=x_j^k+\frac{{\mathrm{\&lambda;}}_k}{q_j}\underset{i=1}{\overset{q_j}{\mathrm{\&Sigma;}}}\frac{b_{s_i}-a_{s_i}^T{x}^k}{{\left|\right|a_{s_i}\left|\right|}^2}{a}_{s_i},j---j=1:n;0<{\mathrm{\&lambda;}}_k<2$

其中q_j为A中第j列不为0的元素的个数；1＜s_i＜m，为A的第j列中第i个不为0的元素所在行的行数；

3)将亚像素位移量δ_k和δ_l量化到小数点后一位或两位，则δ_k和δ_l的取值分别有10种或100种取值，组合起来有100种或10000种取值，对每一种取值计算方程组的系数，得到一张系数表，将其提前计算好存储起来，迭代求方程组时计算出量化后的δ_k和δ_l，然后通过查表得到方程组系数；

4)建立与X维数相同的权值矩阵W和修正矩阵ΔX，其中W中每个元素为对应修正矩阵ΔX中像素点在每次迭代中已被修改的次数，设g为迭代次数，算法步骤为：

(1)初始值任意，设为0：g＝0，X^g＝0；

(2)权值矩阵和修正矩阵置0矩阵：W＝Φ，ΔX＝Φ；

(3)对观察图象的每个像素计算在理想高分辨率图象中的对应点坐标，从而得到k₂、δ_k、l₂、δ_l，再算出h″、Ω_k2，l2，并且

${\mathrm{\&Delta;X}}_{i,j}=\frac{W_{i,j}}{W_{i,j}+1}{\mathrm{\&Delta;X}}_{i,j}+\frac{1}{W_{i,j}+1}\&lsqb;Y_{k_1{l}_1}-h^{\&prime;\&prime;}\&CenterDot;\&times;{X^g}_{{\mathrm{\&Omega;}}_{k_2,l_2}}\&rsqb;h_{i,j}^{\&prime;}/\left|\right|h^{\&prime;\&prime;}(:)|{|}_2^{2}i,j\&Element;{\mathrm{\&Omega;}}_{k_2,l_2}$

                 W_i，j＝W_i，j+1              $i,j\&Element;{\mathrm{\&Omega;}}_{k_2,l_2}$

(4)X^g+1＝X^g+λ_gΔX，g＝g+1

上述算法在并行处理环境中把观察图象的像素分为P组，P为并行计算的路数，对每一组都需定义权值矩阵和修正矩阵，然后对各组分别处理，得到相应的Wⁱ和ΔXⁱ(i＝1∶P)

$X_{k,l}^{g+l}=X_{k,l}^g+{\mathrm{\&lambda;}}_g\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\left((W_{k,l}^i/\underset{j=1}{\overset{P}{\mathrm{\&Sigma;}}}{W}_{k,l}^j){\mathrm{\&Delta;X}}_{k,l}^i\right),k=1:\mathrm{Lm};l=1:\mathrm{Ln}$

Images