CN103310427B - 图像超分辨率与像质增强方法 - Google Patents

Abstract

图像超分辨率与像质增强方法，涉及空间遥感成像系统中的数字图像处理技术领域，解决了现有维纳滤波图像复原法存在的可被利用优化元素少、像质提升空间有限而无法有效提高图像质量的问题，该方法如下：利用光线分束器在两个通道分别采集在焦像和离焦像，并采集N幅相对于在焦像在水平和竖直方向上均错位1/N个像元的错位在焦像；对在焦像和离焦像进行PD处理，采用共轭梯度算法求出像差系数，求出经PD处理复原后目标图像并计算目标函数，对错位在焦像和复原后目标图像进行维纳滤波处理，对得到两幅复原图像再进行亚像元融合得到一幅具有超高分辨率的复原图像。本发明实现了图像的超分辨率和高清晰度，分辨率提高了1.4倍。

Description

图像超分辨率与像质增强方法

技术领域

本发明涉及空间遥感成像系统中的数字图像处理技术领域，具体涉及一种图像超分辨率与像质增强方法。

背景技术

随着空间遥感技术的发展，对遥感成像系统的空间分辨率提出了更高的要求，尽管近年来探测器像元尺寸越做越小，像阵面积越来越大，但仍不能满足实际应用的要求。

目前对于卫星光学遥感成像有效载荷系统像质增强方法研究领域，仍然是采用传统的方法和单一的技术模式，事实上，受探测器像元尺寸的限制，现有的像质增强方法，如维纳滤波图像复原法，通过寻找一个使统计误差函数最小的估计值，由此得到图像的退化函数，根据图像退化函数得到目标图像的频谱，对目标图像的频谱进行逆傅里叶变换可求得目标物，即复原后的目标清晰像，但该方法可被利用优化元素少，像质提升空间有限，此外，空间对地光学遥感成像过程是极为复杂的物理过程，因此，探测器像元尺寸的限制、卫星姿态变化、成像系统的残余像差及大气传输、环境扰动（湍流和随机振动）等因素成为空间遥感成像系统分辨率提高的一个瓶颈，造成了图像质量的下降。

发明内容

为了解决现有的维纳滤波图像复原法存在的可被利用优化元素少、像质提升空间有限而无法有效提高图像质量的问题，本发明提供一种图像超分辨率与像质增强方法。

本发明为解决技术问题所采用的技术方案如下：

图像超分辨率与像质增强方法，该方法的条件和步骤如下：

步骤一、在光学成像系统中，扩展目标的光线经过望远物镜入射至光线分束器，分成两束分别在两个CCD接收图像，第一CCD在焦面位置上接收一幅在焦像，第二CCD在一个已知离焦量的离焦位置上接收一幅离焦像；使光线分束器相对于第一CCD的焦平面阵列方向进行扫描，采集N幅在第一CCD上形成的相对于第一CCD上的在焦像在水平和竖直方向上均错位1/N个像元的错位在焦像，N的取值范围是：2≤N≤10；

步骤二、将步骤一采集的一对在焦像和离焦像传送至DSP图像处理器，分别进行傅立叶变换得到在焦像和离焦像的傅立叶变换谱G₁(u,ν)和G₂(u,ν)，设定光瞳函数P₁和P₂，初始化泽尼克像差系数α_j，设定精度为ε，将光瞳函数P₁和P₂分别进行傅立叶变换得到其所对应的振幅分布函数，取振幅分布函数模的平方，分别计算出在焦像和离焦像的点扩散函数h₁和h₂，将h₁和h₂分别进行傅立叶变换得到在焦像和离焦像的光学传递函数H₁（u,v）和H₂（u,v），推导得出目标函数E(a)的表达式如式（1）所示：

$E\left(a\right)=\underset{u,v\∈{\mathrm{\χ}}_1}{\mathrm{\Σ}}\frac{{|G_1(u,v){{\hat{H}}_2}^{*}(u,v)-G_2(u,v){{\hat{H}}_1}^{*}(u,v)|}^2}{{\left|{\hat{H}}_1(u,v)\right|}^2+{\left|{\hat{H}}_2(u,v)\right|}^2}---\left(1\right)$

式（1）中，H₁ ^*（u,v）为H₁（u,v）的复共轭，H₂ ^*（u,v）为H₂（u,v）的复共轭；采用共轭梯度算法求解出泽尼克像差系数α_j，利用此像差系数α_j求得光学成像系统的像差分布，由此像差分布即可以得到一幅复原后的目标图像，再利用此像差系数α_j求解出光学成像系统的光学传递函数，即为图像的退化函数H(u,v)，对上述得到的复原后的目标图像进行维纳滤波处理，得到一幅二次复原后的目标清晰图像，同时对步骤一中得到的错位在焦像也进行维纳滤波处理，得到一幅复原后的错位在焦像；

步骤三、将步骤二中得到的复原后的目标图像和错位在焦像按照实际位置迭加在一起，得到的初步合成的图像中有反映图像真实信息的实际数据点，相邻实际数据点之间有1个待插值点，采用三次B样条插值算法估计出待插值点的像素值，在初步合成的图像的空间域进行插值处理，得到超分辨率复原图像。

步骤一中，所述光学成像系统还包括：主控计算机；与所述主控计算机相连的精密压电致动直线导轨；所述光线分束器固定在所述精密压电致动直线导轨上，所述主控计算机控制精密压电致动直线导轨运动，同时带动光线分束器移动。

步骤一中，在焦像和离焦像的光学传递函数H₁（u,v）和H₂（u,v）的表达式如式（2）所示：

$H_k(u,v)=C_k(u,v)\⊗C_k(u,v),k=\mathrm{1,2}---\left(2\right)$

式（2）中，代表相关运算，C_k(u,v)为光学成像系统的光学传递函数，其表达式如式（3）所示：

C_k(u,ν)=A(u,ν)exp{i[φ(u,ν)+△φ_k(u,ν)]}k=1,2（3）

式（3）中，A(u,ν)表示光瞳的二值函数，φ(u,ν)表示广义光瞳的相位差函数，表示成泽尼克多项式的形式如式（4）所示：

${\mathrm{\φ}}_n(u,v)=\frac{2\mathrm{\π}}{\mathrm{\λ}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{\mathrm{\α}}_j{z}_j---\left(4\right)$

式（4）中，△φ_k(u,ν)表示引入的离焦产生的相位差函数，

${\mathrm{\Δ\φ}}_k(u,v)=\frac{\mathrm{\π\Δ}{Z}_k}{\mathrm{\λ}{f}^2}(u^2+v^2),k=\mathrm{1,2}---\left(5\right)$

式（5）中，λ为入射光波长，f为光学成像系统的焦距，△Z_k为离焦距离，当k=1时，表示未离焦的情况，即△φ₀(u,ν)=0，当k=2时，Δφ₂(u，v)表示产生某一离焦距离时的相位差；离焦量的表达式如式（6）所示：

式（6）中，△Z_k为离焦距离，λ为入射光波长，F#为系统的F数，即光学成像系统焦距f与孔径直径D的比值。

步骤二中，所述共轭梯度算法的具体过程如下：（1）根据上述的目标函数E(a)的表达式，设定初始值x⁽⁰⁾及精度为ε；（2）计算g⁽⁰⁾为在点x⁽⁰⁾处目标函数E(a)的梯度函数（3）判断，若g⁽⁰⁾的绝对值小于设定的精度ε，则最终解x^*=x⁽⁰⁾，停止算法；（4）若g⁽⁰⁾的绝对值大于设定的精度ε，则令p⁽⁰⁾=-g⁽⁰⁾；（5）采用一维搜索获得最优补偿因子t_k；（6）令x^(k+1)=x^(k)+t_kp^(k)，计算在点x^(k+1)处目标函数E(a)的梯度函数值g^(k+1)；（7）判断，若g^(k+1)的绝对值小于设定的精度ε，则停止迭代，最终解x^*=x^(k+1)；（8）若g^(k+1)的绝对值大于设定的精度ε，则判断k是否等于n；（9）若k等于n，则令x⁽⁰⁾=x^(k+1)，g⁽⁰⁾=g^(k+1)，返回第（4）步继续迭代计算；（10）若k不等于n，则令a=||g^(k+1)||²/||g^(k)||²，p^(k+1)=-g^(k+1)+ap^(k)，返回第（4）步继续迭代计算，直到求解出这组泽尼克像差系数α_j。

步骤三中，所述采用三次B样条插值算法估计出待插值点的像素值的具体过程如下：三次B样条曲线方程的矩阵形式如式（10）所示：

$p_i\left(u\right)=\frac{1}{6}\left[\begin{array}{cccc}1& u& u^2& u^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{d}_i\\ d_{i+1}\\ d_{i+2}\\ d_{i+3}\end{array}\right]---\left(10\right)$

式（10）中，d_i为控制顶点，P_i(u)为对应于参数u的三次B样条曲线上的点；式（10）的三次规范B样条基函数的表达式分别如式（11）所示：

$\left\{\begin{array}{c}{N}_{\mathrm{0,3}}\left(u\right)=\frac{1}{6}{(1-u)}^3\\ N_{\mathrm{1,3}}\left(u\right)=\frac{1}{6}({3u}^3-6u^2+4)\\ N_{\mathrm{2,3}}\left(u\right)=\frac{1}{6}(-3u^3+3u^2+3u+1)\\ N_{\mathrm{3,3}}\left(u\right)=\frac{1}{6}{u}^3\end{array}\right.---\left(11\right)$

给定型值点Q_i，令P_i(1)=P_i+1(0)=Q_i，代入式（10）求解d_i，求解线性方程组：

d_i+4d_i+1+d_i+2＝6Q_i(12)

取d₀=d₁，d_n+1=d_n，代入式（12）求出d_i，并代入式（10）求得通过型值点Q_i的三次均匀B样条插值曲线，由实际数据点通过三次B样条插值得到B样条插值中点的插值，再由B样条插值中点的插值求均值得到待插值点的像素值，式（10）、式（11）和式（12）中，0≤u≤1；i=0，1，2，…，n，n为正整数。

本发明的有益效果是：

本发明基于光学遥感成像链模型和波前分析，将位相相异（PD）方法与亚像元图像融合处理结合在一起，利用多帧图像和移相技术实现图像超分辨率和像质增强，对存在探测器尺寸限制与波前像差扰动的光学成像系统的成像质量进行了全面优化；

1、本发明实现了图像的超分辨率和高清晰度，分辨率提高了1.4倍；

2、本发明不但能够实现超分辨率成像，还能反演出光学成像系统的静态残余像差，补偿和校正动态环境像差扰动（如大气湍流扰动和随机振动等动态像差）对光学成像系统成像质量的影响，如图像模糊等，增强和提高光学成像系统的成像质量，能够获得具有超分辨率和高清晰度特征的图像；

3、本发明所采用的位相相异（PD）方法集光学（硬件）及数字（软件）方法于一身，包含多种元素，像质提升空间大，能够从光学成像基本原理出发分析像质恶化原因，用硬件方法从根本改善光学成像系统成像质量的同时，利用已有成熟的数字图像处理方法，获得更加完善的像质；

4、本发明使用低分辨率的光学成像系统和探测器件（如CCD等）以获得超分辨率的图像，大大降低了高分辨率和高清晰度遥感成像设备的成本；

5、本发明的效果显著且易于实现，具有很好的应用前景和价值。

附图说明

图1为本发明中的特殊设计的光学成像系统的示意图；

图2为采集的在第一CCD上形成的一幅退化的低分辨率在焦像；

图3为采集的在第二CCD上形成的一幅离焦一个波长的低分辨率离焦像；

图4为采集的一幅相对于图2中的在焦像在水平和竖直方向上均错位1/2个像元的错位在焦像；

图5为PD算法流程图；

图6为共轭梯度算法流程图；

图7为对目标图像进行维纳滤波复原后的低分辨率复原图像；

图8为对图4中的错位在焦像进行复原后的错位低分辨率复原图像；

图9为三次B样条插值方法中的由控制顶点求B样条曲线的示意图；

图10为三次B样条插值方法中的由型值点求B样条曲线的示意图；

图11为亚像元重构示意图；

图12为对图7中的低分辨率复原图像和图8中的错位低分辨率复原图像进行亚像元融合后形成的大小为256×256的图像；

图13为对图2中的退化的低分辨率在焦像通过像元复制放大为大小为256×256的图像；

图14为对图7中的低分辨率复原图像通过像元复制放大为大小为256×256的图像。

图中：1、扩展目标，2、望远物镜，3、光线分束器，4、第一CCD，5、第二CCD，6、精密压电致动直线导轨，7、主控计算机，8、扰动波面，9、第一CCD的第一行，10、第二CCD的第一行，11、第一CCD的第二行，12、第二CCD的第二行，13、实际数据点，14、待插值点，15、B样条插值中点，16、像素值。

具体实施方式

以下结合附图对本发明作进一步详细说明。

本发明的一种图像超分辨率与像质增强方法是通过以下步骤实现的：

一、利用特殊设计的光学成像系统采集一对在焦像和离焦像及N幅错位在焦像

1、搭建特殊设计的光学成像系统

如图1所示，该光学成像系统包括：扩展目标1、望远物镜2、光线分束器3、第一CCD4、第二CCD5、精密压电致动直线导轨6和主控计算机7，扩展目标1发出的光线经过扰动波面8入射至望远物镜2，再经过望远物镜2出射至光线分束器3，经光线分束器3后分成两束出射，分别由通道一的第一CCD4和通道二的第二CCD5接收图像，第一CCD4和第二CCD5分别位于两束光线形成的焦平面位置上，上述的光线分束器3固定在精密压电致动直线导轨6上，而精密压电致动直线导轨6与主控计算机7相连，主控计算机7可以控制精密压电致动直线导轨6运动，并可以同时带动光线分束器3移动；

2、采集一对在焦像和离焦像

首先，在光线分束器3移动之前，采集在第一CCD4上形成的一幅在焦像，采集的退化的低分辨率在焦像如图2所示，图像大小为128×128；然后，根据事先计算或设定的离焦量，主控计算机7控制精密压电致动直线导轨6运动，同时使光线分束器3产生垂直于第一CCD4的焦平面阵列方向上的离焦位移动作，本实施方式中，设定离焦量为一个波长，使第二CCD5位于离焦一个波长的离焦位置上，采集此时在第二CCD5上形成的一幅离焦一个波长的离焦像，采集的低分辨率离焦像如图3所示，图像大小为128×128；如果连续获取在焦像和离焦像的时间序列，不但可以提高位相信息的提取精度，还可用于动态像差扰动的反演和图像复原；

3、采集N幅错位在焦像

主控计算机7控制精密压电致动直线导轨6运动同时带动光线分束器3移动，使光线分束器3产生相对于第一CCD4的焦平面阵列方向上的精密扫描运动，采集N幅在第一CCD4上形成的相对于第一CCD4上的在焦像（即图2所示的图像）在水平和竖直方向上均错位1/N个像元的错位在焦像，N次位移对应生成N幅错位在焦像，N的取值范围是：2≤N≤10，本实施方式中，N取2，即采集两幅在第一CCD4上形成的相对于第一CCD4上的在焦像在水平和竖直方向上均错位1/2个像元的错位在焦像，其中的一幅低分辨率错位在焦像如图4所示；光线分束器3在相对于第一CCD4的焦平面阵列方向上的精密扫描运动可以实现小于第一CCD4焦平面阵元（亚像元、1/N个像元）位移变化。

二、基于位相相异（PD）方法的图像复原

1、推倒目标函数E(a)的表达式

如图5所示，将采集的一对在焦像（如图2所示）和离焦像（如图3所示）传送至高速DSP图像处理器，高速DSP图像处理器选择型号为TMS320DM642的DSP芯片，分别进行傅立叶变换得到在焦像的傅立叶变换谱G₁(u,ν)和离焦像的傅立叶变换谱G₂(u,ν)，设定光瞳函数P₁和P₂，初始化泽尼克（zernike）像差系数α_j，设定终止条件，即最优化方法的设定精度ε，将光瞳函数P₁和P₂分别进行傅立叶变换得到其所对应的振幅分布函数，取振幅分布函数模的平方，分别计算出在焦像的点扩散函数（PSF）h₁和离焦像的点扩散函数（PSF）h₂，将这两个点扩散函数h₁和h₂分别进行傅立叶变换得到其对应的光学传递函数（OTF）H₁（u,v）和H₂（u,v），然后推导得出目标函数E(a)的表达式如式（1）所示，目标函数E(a)是一组泽尼克（zernike）像差系数α_j的函数：

$E\left(a\right)=\underset{u,v\∈{\mathrm{\χ}}_1}{\mathrm{\Σ}}\frac{{|G_1(u,v){{\hat{H}}_2}^{*}(u,v)-G_2(u,v){{\hat{H}}_1}^{*}(u,v)|}^2}{{\left|{\hat{H}}_1(u,v)\right|}^2+{\left|{\hat{H}}_2(u,v)\right|}^2}---\left(1\right)$

式（1）中，G₁(u,ν)为在焦像的傅立叶变换谱，G₂(u,ν)为离焦像的傅立叶变换谱，H₁(u,ν)为未离焦的光学成像系统的光学传递函数（OTF），H₁ ^*（u,v）为H₁（u,v）的复共轭，H₂(u,ν)为产生某一离焦距离时光学成像系统的光学传递函数（OTF），H₂ ^*（u,v）为H₂（u,v）的复共轭；H_k(u,ν)（包括H₁(u,ν)和H₂(u,ν)）是利用光学成像系统的相干传递函数（CTF）进行自相关运算得到的，其表达式如式（2）所示：

$H_k(u,v)=C_k(u,v)\⊗C_k(u,v),k=\mathrm{1,2}---\left(2\right)$

式（2）中，代表相关运算，光学成像系统的光学传递函数（CTF）的表达式如式（3）所示：

C_k(u,ν)=A(u,ν)exp{i[φ(u,ν)+△φ_k(u,ν)]}k=1,2（3）

式（3）中，A(u,ν)表示光瞳的二值函数，φ(u,ν)表示广义光瞳的相位差函数，表示为泽尼克（zernike）多项式的形式如式（4）所示：

${\mathrm{\φ}}_n(u,v)=\frac{2\mathrm{\π}}{\mathrm{\λ}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{\mathrm{\α}}_j{z}_j---\left(4\right)$

式（4）中，△φ_k(u,ν)表示引入的离焦产生的相位差函数，

${\mathrm{\Δ\φ}}_k(u,v)=\frac{\mathrm{\π\Δ}{Z}_k}{\mathrm{\λ}{f}^2}(u^2+v^2),k=\mathrm{1,2}---\left(5\right)$

式（5）中，λ为入射光波长，f为光学成像系统的焦距，△Z_k为离焦距离，当k=1时，表示未离焦的情况，即△φ₀(u,ν)=0，当k=2时，Δφ₂(u，v)表示产生某一离焦距离时的相位差；离焦量的表达式如式（6）所示：

式（6）中，△Z_k为离焦距离，λ为入射光波长，F#为系统的F数，即光学成像系统焦距f与孔径直径D的比值；

2、寻求使目标函数E(a)最小化的最优解

本实施方式中所采用的最优化方法—共轭梯度算法是一种非线性最优化方法，它是以二次目标函数为模型，用梯度来构造共轭方向的算法，可用于求解一般的非线性目标函数的无约束问题，为了使算法更具有鲁棒性，引入了“重新开始”策略；共轭梯度算法的流程如图6所示，（1）根据上述的目标函数E(a)的表达式，设定初始值x⁽⁰⁾及精度为ε；（2）计算g⁽⁰⁾为在点x⁽⁰⁾处目标函数E(a)的梯度函数（3）判断，若g⁽⁰⁾的绝对值小于设定的精度ε，则最终解x^*=x⁽⁰⁾，停止算法；（4）若g⁽⁰⁾的绝对值大于设定的精度ε，则令p⁽⁰⁾=-g⁽⁰⁾；（5）采用一维搜索获得最优补偿因子t_k；（6）令x^(k+1)=x^(k)+t_kp^(k)，计算在点x^(k+1)处目标函数E(a)的梯度函数值g^(k+1)；（7）判断，若g^(k+1)的绝对值小于设定的精度ε，则停止迭代，最终解x^*=x^(k+1)；（8）若g^(k+1)的绝对值大于设定的精度ε，则判断k是否等于n；（9）若k等于n，则令x⁽⁰⁾=x^(k+1)，g⁽⁰⁾=g^(k+1)，返回第（4）步继续迭代计算；（10）若k不等于n，则令a=||g^(k+1)||²/||g^(k)||²，p^(k+1)=-g^(k+1)+ap^(k)，返回第（4）步继续迭代计算，直到求解出这组泽尼克（zernike）像差系数α_j，得到的像差系数α_j即为光学成像系统的像差系数；

3、计算目标图像并估计退化函数

利用上述得到的泽尼克（zernike）像差系数α_j，可以得到光学成像系统的像差分布，由此像差分布即可以得到一幅复原后的目标图像；

再利用上述得到的泽尼克（zernike）像差系数α_j求解出对应的光学成像系统的光学传递函数（OTF），此光学传递函数（OTF）即为退化图像的退化函数H(u,v)；

4、对图像的维纳滤波处理

对上述得到的复原后的目标图像进行维纳滤波处理，得到一幅二次复原后的目标清晰图像，如图7所示；

对步骤一中得到的错位在焦像（如图4所示）进行维纳滤波处理，得到一幅复原后的错位在焦像（如图8所示）；

维纳滤波处理的具体过程如下：采用维纳滤波器寻找一个使统计误差函数最小的估计值统计误差函数的表达式如式（7）所示：

$e^2=E\left\{{(f-\hat{f})}^2\right\}---\left(7\right)$

式（7）中，E是期望值操作符，f为未退化的图像，为未退化的图像f的估计值，即退化图像，式（7）在频域可表示为如式（8）所示的形式：

$\hat{F}(u,v)=\left[\frac{H^{*}(u,v)}{{\left|H(u,v)\right|}^2+K}\right]W(u,v)---\left(8\right)$

式（8）中，W(u,v)是退化图像的频谱，H(u,v)为退化图像的退化函数，K是噪信功率比，是未知的，用一个常数代入，由上述得到的退化函数H(u,v)求得目标图像的频谱W(u,v)，对目标图像的频谱W(u,v)进行逆傅立叶变换可求得二次复原后的目标图像，对错位在焦像的维纳滤波处理的过程与上述相同。

利用PD方法，光学成像系统的波前探测与图像复原问题演变为寻求使目标函数最小化的最优解问题，经过PD方法，在存在光学像差的光学成像系统中，既可以获得复原后的目标清晰像，又可以得到光学成像系统的像差分布，去除了光学成像系统及大气湍流扰动等动态像差造成的图像模糊等。

三、亚像元图像融合

如图11所示，将上述得到的经过二次复原后的目标清晰像（如图7所示）和复原后的错位在焦像（如图8所示）按照实际位置迭加在一起，得到一幅初步合成的图像，初步合成的图像中既有反映图像真实信息的实际数据点13，在第一CCD4的第一行9与第一CCD4的第二行11之间、第二CCD5的第一行10与第二CCD5的第二行12之间，即在每两个同行或者同列的相邻实际数据点13之间又存在1个空缺点，即待插值点14，利用已有的图像真实信息估计出待插值点14的像素值，然后在初步合成的图像的空间域进行插值处理，即可复原出一幅具有超分辨率特征的新图像，待插值点14的像素值估计的越精确则复原图像的分辨率越高；

本实施方式中采用三次B样条插值算法估计出待插值点14的像素值，B样条曲线具有很强的曲线表达能力，在计算机图形学上有着广泛的应用，是重要的图像插值工具，一般情况下，K次B样条曲线的方程如式（9）所示：

$P\left(u\right)=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{d}_i{N}_{i,k}\left(u\right)---\left(9\right)$

式（9）中，d_i为控制顶点，i=0，1，2，…，n，n为正整数，将d_i顺序连接而成的折线称为B样条控制多边形，N_i,k(u)为K次规范B样条基函数，其中每一个称为规范B样条，简称B样条，它是由一个称为节点矢量的非递减参数的序列U:u₀≤u₁≤L≤u_n+k+1所决定的k次分段多项式，也就是k次多项式样条，P(u)是对应于参数u的B样条曲线上的点，对于三次B样条曲线方程可写为下面的矩阵形式，

$p_i\left(u\right)=\frac{1}{6}\left[\begin{array}{cccc}1& u& u^2& u^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{d}_i\\ d_{i+1}\\ d_{i+2}\\ d_{i+3}\end{array}\right]---\left(10\right)$

其中，0≤u≤1；i=0，1，2，…，n，n为正整数；式（10）的三次规范B样条基函数的表达式分别如式（11）所示：

$\left\{\begin{array}{c}{N}_{\mathrm{0,3}}\left(u\right)=\frac{1}{6}{(1-u)}^3\\ N_{\mathrm{1,3}}\left(u\right)=\frac{1}{6}({3u}^3-6u^2+4)\\ N_{\mathrm{2,3}}\left(u\right)=\frac{1}{6}(-3u^3+3u^2+3u+1)\\ N_{\mathrm{3,3}}\left(u\right)=\frac{1}{6}{u}^3\end{array}\right.---\left(11\right)$

由式（10）可知，每一曲线段是由四个顺序控制顶点定义的三次曲线，如果已知控制顶点d和节点矢量U，求曲线上的点P(u)，称为“正算”，如图9所示；如果己知曲线上的一系列点(又称为型值点)和节点矢量U，求出控制顶点d，再由控制顶点d求出这条曲线，称为“反算”，如图10所示；

本实施方式中采用的B样条插值算法属于第二种“反算”的过程，即已知曲线上的点（第一CCD4和第二CCD5的实际数据点13），反算出控制顶点d，再由控制顶点d计算出曲线，进一步计算出待插值点14的像素值，若给定型值点Q_i，其反算B样条曲线的控制顶点d即是令P_i(1)=P_i+1(0)=Q_i，i=0，1，2，…，n，n为正整数，代入式（10）求解d_i，求解线性方程组：

d_i+4d_i+1+d_i+2＝6Q_i(12)

式（12）中，i=0，1，2，…，n；n为正整数，对于非封闭曲线，式（12）共有n+3个未知数，但只有n+1个方程，因此需要再补充两个端点条件，可取d₀=d₁，d_n+1=d_n，将其代入式（12）求出d_i，并代入式（10）即可求得通过型值点Q_i的三次均匀B样条插值曲线，如图11所示，由1，2，3，4点（第一CCD4和第二CCD5的实际数据点13）通过B样条插值得到a点的插值，由5，6，7，8点得到b点的插值，由9，2，6，10点得到c点的插值，由11，3，7，12点得到d点的插值，然后再由a，b，c，d求均值最终得到待插值点14的像素值16，由于通过B样条插值得到的a，b，c，d点比CCD器件的实际数据点2，3，6，7更接近于待插值点14，因此，此方法可以得到更好的插值图像，图12为对图7中的低分辨率复原图像和图8中的错位低分辨率复原图像进行亚像元融合后形成的大小为256×256的图像；图13为对图2中的退化的低分辨率在焦像通过像元复制放大为大小为256×256的图像；图14为对图7中的低分辨率复原图像通过像元复制放大为大小为256×256的图像，经过比较，未经过处理的退化的低分辨率在焦像（如图13所示）只能分辨出9号线，而经过处理后的超分辨率图像（如图12所示）可以分辨出12号线，分辨率提高了1.4倍。

对复原后的图像进行插值处理以获得更高分辨率的清晰像，因为在复原前两幅图像在水平和竖直方向上均已错位1/N个像元，所以复原后的两幅图像在水平和竖直方向上同样错位1/N个像元，对两幅复原后的图像进行融合即可获得高分辨率的复原图像。

Claims (5)

1.图像超分辨率与像质增强方法，其特征在于，该方法的条件和步骤如下：

步骤一、在光学成像系统中，扩展目标的光线经过望远物镜入射至光线分束器，分成两束分别在两个CCD接收图像，第一CCD在焦面位置上接收一幅在焦像，第二CCD在一个已知离焦量的离焦位置上接收一幅离焦像；使光线分束器相对于第一CCD的焦平面阵列方向进行扫描，采集N幅在第一CCD上形成的相对于第一CCD上的在焦像在水平和竖直方向上均错位1/N个像元的错位在焦像，N的取值范围是：2≤N≤10；

步骤二、将步骤一采集的一对在焦像和离焦像传送至DSP图像处理器，分别进行傅立叶变换得到在焦像和离焦像的傅立叶变换谱G₁(u,ν)和G₂(u,ν)，设定光瞳函数P₁和P₂，初始化泽尼克像差系数α_j，设定精度为ε，将光瞳函数P₁和P₂分别进行傅立叶变换得到其所对应的振幅分布函数，取振幅分布函数模的平方，分别计算出在焦像和离焦像的点扩散函数h₁和h₂，将h₁和h₂分别进行傅立叶变换得到在焦像和离焦像的光学传递函数H₁(u,v)和H₂(u,v)，推导得出目标函数E(a)的表达式如式(1)所示：

$E\left(a\right)=\underset{u,v\∈{\mathrm{\χ}}_1}{\Σ}\frac{|G_1(u,v){H_2}^{*}(u,v)-G_2(u,v){H_1}^{*}(u,v){|}^2}{|H_1(u,v){|}^2+|H_2(u,v){|}^2}---\left(1\right)$

式(1)中，H₁ ^*(u,v)为H₁(u,v)的复共轭，H₂ ^*(u,v)为H₂(u,v)的复共轭；采用共轭梯度算法求解出泽尼克像差系数α_j，利用此像差系数α_j求得光学成像系统的像差分布，由此像差分布即可以得到一幅复原后的目标图像，再利用此像差系数α_j求解出光学成像系统的光学传递函数，即为图像的退化函数H(u,v)，对上述得到的复原后的目标图像进行维纳滤波处理，得到一幅二次复原后的目标清晰图像，同时对步骤一中得到的错位在焦像也进行维纳滤波处理，得到一幅复原后的错位在焦像；

步骤三、将步骤二中得到的复原后的目标图像和错位在焦像按照实际位置迭加在一起，得到的初步合成的图像中有反映图像真实信息的实际数据点，相邻实际数据点之间有1个待插值点，采用三次B样条插值算法估计出待插值点的像素值，在初步合成的图像的空间域进行插值处理，得到超分辨率复原图像。

2.根据权利要求1所述的图像超分辨率与像质增强方法，其特征在于，步骤一中，所述光学成像系统还包括：主控计算机；与所述主控计算机相连的精密压电致动直线导轨；所述光线分束器固定在所述精密压电致动直线导轨上，所述主控计算机控制精密压电致动直线导轨运动，同时带动光线分束器移动。

3.根据权利要求1所述的图像超分辨率与像质增强方法，其特征在于，步骤一中，在焦像和离焦像的光学传递函数H₁(u,v)和H₂(u,v)的表达式如式(2)所示：

$H_k(u,v)=C_k(u,v)\⊗C_k(u,v),k=1,2---\left(2\right)$

式(2)中，代表相关运算，C_k(u,v)为光学成像系统的光学传递函数，其表达式如式(3)所示：

C_k(u,ν)＝A(u,ν)exp{i[φ_k(u,ν)+Δφ_k(u,ν)]}k＝1,2(3)

式(3)中，A(u,ν)表示光瞳的二值函数，φ(u,ν)表示广义光瞳的相位差函数，表示成泽尼克多项式的形式如式(4)所示：

${\mathrm{\φ}}_k(u,\mathrm{\ν})=\frac{2\mathrm{\π}}{\mathrm{\λ}}\underset{j=1}{\overset{J}{\Σ}}{\mathrm{\α}}_j{z}_j---\left(4\right)$

式(4)中，λ为入射光波长，Z_j泽尼克多项式函数，α_j为泽尼克像差系数，J为整数；

式(3)中，Δφ_k(u,ν)表示引入的离焦产生的相位差函数，

${\mathrm{\Δ\φ}}_k(u,v)=\frac{{\mathrm{\π\ΔZ}}_k}{{\mathrm{\λf}}^2}(u^2+v^2),k=1,2---\left(5\right)$

式(5)中，λ为入射光波长，f为光学成像系统的焦距，ΔZ_k为离焦距离，当k＝1时，表示未离焦的情况，即Δφ₁(u,ν)＝0，当k＝2时，Δφ₂(u，v)表示产生某一离焦距离时的相位差；离焦量的表达式如式(6)所示：

式(6)中，ΔZ_k为离焦距离，λ为入射光波长，F#为系统的F数，即光学成像系统焦距f与孔径直径D的比值。

4.根据权利要求1所述的图像超分辨率与像质增强方法，其特征在于，步骤二中，所述共轭梯度算法的具体过程如下：(1)根据上述的目标函数E(a)的表达式，设定初始值x⁽⁰⁾及精度为ε；(2)计算g⁽⁰⁾＝▽f(x⁰)，g⁽⁰⁾为在点x⁽⁰⁾处目标函数E(a)的梯度函数▽f(x⁰)；(3)判断，若g⁽⁰⁾的绝对值小于设定的精度ε，则最终解x^*＝x⁽⁰⁾，停止算法；(4)若g⁽⁰⁾的绝对值大于设定的精度ε，则令p⁽⁰⁾＝-g⁽⁰⁾；(5)采用一维搜索获得最优补偿因子t_k；(6)令x^(k+1)＝x^(k)+t_kp^(k)，计算在点x^(k+1)处目标函数E(a)的梯度函数值g^(k+1)；(7)判断，若g^(k+1)的绝对值小于设定的精度ε，则停止迭代，最终解x^*＝x^(k+1)；(8)若g^(k+1)的绝对值大于设定的精度ε，则判断k是否等于n；(9)若k等于n，则令x⁽⁰⁾＝x^(k+1)，g⁽⁰⁾＝g^(k+1)，返回第(4)步继续迭代计算；(10)若k不等于n，则令a＝||g^(k+1)||²/||g^(k)||²，p^(k+1)＝-g^(k+1)+ap^(k)，返回第(4)步继续迭代计算，直到求解出这组泽尼克像差系数α_j。

5.根据权利要求1所述的图像超分辨率与像质增强方法，其特征在于，步骤三中，所述采用三次B样条插值算法估计出待插值点的像素值的具体过程如下：三次B样条曲线方程的矩阵形式如式(10)所示：

$p_i\left(u\right)=\frac{1}{6}\left[\begin{array}{cccc}1& u& u^2& u^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{d}_i\\ d_{i+1}\\ d_{i+2}\\ d_{i+3}\end{array}\right]---\left(10\right)$

式(10)中，d_i为控制顶点，P_i(u)为对应于参数u的三次B样条曲线上的点；式(10)的三次规范B样条基函数的表达式分别如式(11)所示：

$\left\{\begin{array}{c}{N}_{0,3}\left(u\right)=\frac{1}{6}{(1-u)}^3\\ N_{1,3}\left(u\right)=\frac{1}{6}(3u^3-6u^2+4)\\ N_{2,3}\left(u\right)=\frac{1}{6}(-3u^3+3u^2+3u+1)\\ N_{3,3}\left(u\right)=\frac{1}{6}{u}^3\end{array}\right.---\left(11\right)$

给定型值点Q_i，令P_i(1)＝P_i+1(0)＝Q_i，代入式(10)求解d_i，求解线性方程组：

d_i+4d_i+1+d_i+2＝6Q_i(12)

取d₀＝d₁，d_n+1＝d_n，代入式(12)求出d_i，并代入式(10)求得通过型值点Q_i的三次均匀B样条插值曲线，由实际数据点通过三次B样条插值得到B样条插值中点的插值，再由B样条插值中点的插值求均值得到待插值点的像素值，式(10)、式(11)和式(12)中，0≤u≤1；i＝0，1，2，…，n，n为正整数。