CN103236044A - 基于泰勒展开积分退化模型的微扫描图像重构方法 - Google Patents

Abstract

本发明提出了一种基于泰勒展开积分退化模型的微扫描图像重构方法，首先基于Taylor展开建立了焦平面探测器的积分退化模型，进而到积分退化函数，利用积分退化函数构建的Wiener滤波器，对过采样图像进行滤波；本发明在进行泰勒（Taylor）展开得到积分退化函数的过程中，考虑了焦平面探测器的占空比，克服了现有亚像元重构算法中假设焦平面探测器的占空比为100%的限制；利用积分退化函数构建的Wiener滤波器对过采样图像进行滤波，抑制了焦平面探测器的积分退化效应，在帧间过采样图像重构的基础上进一步提高图像质量，从而进一步提高微扫描成像系统的成像性能。

Description

基于泰勒展开积分退化模型的微扫描图像重构方法

技术领域

本发明涉及一种基于泰勒展开积分退化模型的微扫描图像重构方法，属于数字图像处理中的图像重构领域。 

背景技术

随着科学技术的发展，CCD(电荷耦合器件)、CMOS(互补金属-氧化物-半导体)以及红外焦平面阵列等面阵探测器越来越多地应用在工业生产、军事侦察、卫星遥感等诸多领域。随着加工工艺的不断发展，CCD等面阵探测器的像元尺寸逐渐变小，面阵规模逐渐变大。像元尺寸的越小，越接近脉冲采样，然而其集光能力越差，散粒噪声越严重，这已经成为成像系统分辨力的瓶颈。此外，相对于CCD和CMOS，红外焦平面阵列由于材料的限制，像元尺寸更大，所以其积分效应更为明显。由于焦平面探测器的像元具有有限大小，因此对应像素值与辐射到整个像元上的能量和成正比，这使像元无法区分更精细的局部图像，即像元对局部能量分布的积分过程使得采集图像发生模糊退化，我们将这种退化现象称之为积分退化效应。 

为了克服焦平面探测器离散采样的限制，国内外先后出现了多种多帧图像超分辨率重建技术，例如频域法、凸集投影法、最大后验概率法、非均匀插值法以及学习法。然而，这些算法复杂性高，难以实现实时处理。微扫描成像系统是一种利用微扫描机构获得多帧具有帧间微位移的成像装置，它在不改变原有成像系统基本组成的情况下，只在光路中添加精密的微扫描装置即可用最低的成本提高系统分辨率。 

针对微扫描成像系统获得的图像，已经提出了多种亚像元重构算法。大多 数亚像元重构算法都假设焦平面探测器的占空比为100％，或者假设局部图像灰度相等。然而，这些假设都与实际情况不符。由于焦平面探测器占空比的存在，使得上述假设出现大量噪声，甚至使得重构算法失效。为此，考虑占空比条件下的微扫描图像重构是亟待解决的问题。 

发明内容

本发明提出了一种基于泰勒展开积分退化模型的微扫描图像重构方法，首先基于Taylor展开建立了焦平面探测器的积分退化模型，进而到积分退化函数，利用积分退化函数构建的Wiener滤波器进行滤波，抑制了焦平面探测器的积分退化效应。 

该方法包括如下步骤： 

第一步：利用基于离散傅里叶变换DFT的亚像元图像配准算法，获得多帧微扫描图像的帧间位移量； 

第二步：在步骤一获得的帧间位移量的基础上，利用基于帧间差分过采样技术，将多帧微扫描图像重建出一幅过采样图像g_os[m，n]； 

第三步：建立焦平面探测器的积分退化模型，进行泰勒(Taylor)展开，得到积分退化函数H[u，v]，用其构建Wiener滤波器M_W[u，v]，经直流分量归一化之后得到归化后Wiener滤波器M′_W[u，v]，用其对第二步得到的过采样图像g_os[m，n]进行Wiener滤波，复原出重构图像

完成微扫描图像的重构；其中[u，v]为频域坐标，u为频域上的水平坐标，v为频域上的竖直坐标； 

其中得到积分退化函数H[u，v]的过程为： 

假设f(x，y)为经过光学系统成像在焦平面探测器上的模拟图像，g[m，n]为经焦平面探测器积分和采样后得到的退化图像，则 

$\begin{array}{c}g[m,n]=\frac{1}{\mathrm{pq}}{\∫}_{\mathrm{nb}-q/2}^{\mathrm{nb}+q/2}{\∫}_{\mathrm{ma}-p/2}^{\mathrm{ma}+p/2}f(x,y)\mathrm{dxdy}\\ =\frac{1}{\mathrm{pq}}{\∫}_{-q/2}^{q/2}{\∫}_{-p/2}^{p/2}f(x+\mathrm{ma},y+\mathrm{nb})\mathrm{dxdy}\end{array}m=\mathrm{1,2},...,M,n=\mathrm{1,2},...,N---\left(1\right)$

式中，a和b为焦平面探测器的水平和竖直像元间隔，p为像元在水平方向的尺寸，q为像元在竖直方向的尺寸，M×N为焦平面探测器像面规格，M为焦平面探测器像面在水平方向的尺寸，N为焦平面探测器像面在竖直方向的尺寸，(x，y)是模拟图像坐标系中的坐标，[m，n]是离散图像坐标系中的坐标，m为离散图像坐标系中的水平坐标，最大取值范围为焦平面探测器像面在水平方向的尺寸，n为离散图像坐标系中的竖直坐标，最大取值范围为焦平面探测器像面在竖直方向的尺寸； 

对式(1)中f(x+ma，y+nb)关于f[m，n]进行2阶泰勒展开得 

$f(x+\mathrm{ma},y+\mathrm{nb})\≅f[m,n]+f_m[m,n]x+f_n[m,n]y$

$+\frac{f_{m^2}[m,n]x^2+f_{n^2}[m,n]y^2+2f_{\mathrm{mn}}[m,n]\mathrm{xy}}{2}---\left(2\right)$

将式(2)代入式(1)得 

$g[m,n]\≅f[m,n]+\frac{p^2}{24}{f}_{m^2}[m,n]+\frac{q^2}{24}{f}_{n^2}[m,n]---\left(3\right)$

式中，

和分别表示图像f[m，n]在m和n方向上的2阶偏导数，f_m[m，n]和f_n[m，n]分别表示图像f[m，n]在m和n方向上的1阶偏导数，f_mn[m，n]表示图橡f[mm，m]先在m方向上取1阶偏导数再在n方向上取1阶偏导数；为了抑制噪声，在数值上用式(4)近似

和

$\left\{\begin{array}{c}{f}_{m^2}[m,n]=\left\{f\right[m-1,n-1]-2f[m,n-1]+f[m+1,n-1]\\ +2\left(f\right[m-1,n]-2f[m,n]+f[m+1,n\left]\right)\\ +f[m-1,n+1]-2f[m,n+1]+f[m+1,n+1]\}/4\\ f_{n^2}[m,n]=\{\left(f\right[m-1,n-1]-2f[m-1,n]+f[m-1,n+1\left]\right)\\ +2\left(f\right[m,n-1]-2f[m,n]+f[m,n+1\left]\right)\\ +\left(f\right[m+1,n-1]-2f[m+1,n]+f[m+1,n+1\left]\right)\}/4\end{array}\right.---\left(4\right)$

将式(4)代入式(3)，再进行DFT得 

$G[u,v]=F[u,v][1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)]---\left(5\right)$

式中，F[u，v]和G[u，v]分别表示f[m，n]和g[m，n]的DFT，α＝p/a和β＝q/b分别表示焦平面探测器在水平和竖直方向的占空比； 

因此，焦平面探测器的积分退化函数为H[u，v] 

$H[u,v]=1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)---\left(6\right)$

Wiener滤波过程为： 

基于积分退化函数H[u，v]，构建Wiener滤波器M_W[uv] 

$M_W[u,v]=\frac{H^{*}[u,v]}{H^2[u,v]+\mathrm{\Γ}}---\left(7\right)$

式中H^*[u，v]为H[u，v]的复共轭；Γ为设定的常数；将M_W[u，v]用直流分量归一化后得到归一化后Wiener滤波器M′_W[u，v] 

$M_W^{\′}[u,v]=\frac{M_W[u,v]}{M_W\left[\mathrm{0,0}\right]}---\left(8\right)$

式中，M_W[0，0]是M_W[u，v]在u＝0，v＝0情况下的值，即M_W[u，v]的直流分量； 

步骤303、利用式(8)描述的Wiener滤波器M′_W[u，v]，对第二步得到的过采样图像g_os[m，n]进行滤波，获得具有重建图像

式中，表示DFT运算，

表示DFT逆运算。 

本发明的有益效果： 

本发明在考虑了焦平在探测器的占空比条件下建立焦平面探测器的积分退化模型，并对其进行泰勒(Taylor)展开得到积分退化函数，克服了现有亚像元重构算法中假设焦平面探测器的占空比为100％的限制；利用积分退化函数构建的Wiener滤波器进行滤波，抑制了焦平面探测器的积分退化效应，在帧间过采样图像重构算法的基础上进一步提高图像质量，从而进一步提高微扫描成像系统的成像性能。 

附图说明

图1是本发明的流程图； 

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图及具体实施例对本发明做进一步详细描述。 

如图1，下面以2×2微扫描图像重构为例，对本发明做进一步说明： 

步骤一、利用基于DFT的亚像素图像配准算法，计算2、3、4帧微扫描图像相对于第1帧图像的帧间微位移量； 

步骤二、在步骤一获得的帧间位移量的基础上，利用基于帧间差分的过采 样重构算法，将4帧微扫描重建出一帧过采样图像g_os[m，n]； 

步骤三、基于Taylor展开建立焦平面探测器的积分退化模型，用其构建Wiener滤波器M_W[u，v]，经直流分量归一化后，对过采样图像g_os[m，n]进行Wiener滤波，抑制图像中因焦平面探测器有限像元尺寸引起的积分退化效应，复原出更高图像质量的高分辨率图像。 

步骤301、建立焦平面探测器的积分退化模型，获得积分退化函数H[u，v]。 

假设f(x，y)为经过光学系统成像在焦平面探测器上的模拟图像，g[m，n]为经焦平面探测器积分和采样后得到的退化图像，则 

$\begin{array}{c}g[m,n]=\frac{1}{\mathrm{pq}}{\∫}_{\mathrm{nb}-q/2}^{\mathrm{nb}+q/2}{\∫}_{\mathrm{ma}-p/2}^{\mathrm{ma}+p/2}f(x,y)\mathrm{dxdy}\\ =\frac{1}{\mathrm{pq}}{\∫}_{-q/2}^{q/2}{\∫}_{-p/2}^{p/2}f(x+\mathrm{ma},y+\mathrm{nb})\mathrm{dxdy}\end{array}m=\mathrm{1,2},...,M,n=\mathrm{1,2},...,N---\left(1\right)$

式中，a和b为焦平面探测器的水平和竖直像元间隔，p×q为像元尺寸，M×N为焦平面探测器像面规格，(x，y)是模拟图像坐标系中的坐标，[m，n]是离散图像坐标系中的坐标。 

对式(1)中f(x+ma，y+nb)关于f[m，n]进行2阶泰勒展开得 

$f(x+\mathrm{ma},y+\mathrm{nb})\≅f[m,n]+f_m[m,n]x+f_n[m,n]y$

$+\frac{f_{m^2}[m,n]x^2+f_{n^2}[m,n]y^2+2f_{\mathrm{mn}}[m,n]\mathrm{xy}}{2}---\left(2\right)$

将式(2)代入式(1)得 

$g[m,n]\≅f[m,n]+\frac{p^2}{24}{f}_{m^2}[m,n]+\frac{q^2}{24}{f}_{n^2}[m,n]---\left(3\right)$

式中，

和

分别表示图像f[m，n]在m和n方向上的2阶偏导数。为了抑制噪声，可以用邻域加权平均计算偏导数，在数值上可用式(4)近似  $f_{m^2}[m,n]$ 和 $f_{n^2}[m,n],$

$\left\{\begin{array}{c}{f}_{m^2}[m,n]=\left\{f\right[m-1,n-1]-2f[m,n-1]+f[m+1,n-1]\\ +2\left(f\right[m-1,n]-2f[m,n]+f[m+1,n\left]\right)\\ +f[m-1,n+1]-2f[m,n+1]+f[m+1,n+1]\}/4\\ f_{n^2}[m,n]=\{\left(f\right[m-1,n-1]-2f[m-1,n]+f[m-1,n+1\left]\right)\\ +2\left(f\right[m,n-1]-2f[m,n]+f[m,n+1\left]\right)\\ +\left(f\right[m+1,n-1]-2f[m+1,n]+f[m+1,n+1\left]\right)\}/4\end{array}\right.---\left(4\right)$

将式(4)代入式(3)，再进行DFT得 

$G[u,v]=F[u,v][1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)]---\left(5\right)$

式中，F[uv]和G[u，v]分别表示

和g[m，n]的DFT，[u，v]是频域坐标；α＝p/a和β＝q/b分别表示焦平面探测器在水平和竖直方向的占空比。 

因此，焦平面探测器的积分退化函数H[u，v]为 

$H[u,v]=1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)---\left(6\right)$

步骤302、基于式(6)描述的积分退化模型，构建Wiener滤波器M_W[u，v] 

$M_W[u,v]=\frac{H^{*}[u,v]}{H^2[u,v]+\mathrm{\Γ}}---\left(7\right)$

式中，H[u，v]为退化函数；H^*[u，v]为H[u，v]的复共轭；Γ为常数，它并不随频率变化，而随着占空比、噪声等级以及图像内容变化；[u，v]是频域坐标。将M_W[u，v]用直流分量归一化后得到归一化后Wiener滤波器M′_W[u，v] 

$M_W^{\′}[u,v]=\frac{M_W[u,v]}{M_W\left[\mathrm{0,0}\right]}---\left(8\right)$

式中，M_W[0，0]是M_W[u，v]在u＝0，v＝0情况下的值，即M_W[u，v]的直流分量。 

步骤303、利用直流分量归一化后的Wiener滤波器M′_W[u，v]，对过采样图 像g_os[m，n]进行滤波，获得具有更好图像质量的重建图像

式中，

表示DFT运算，

表示DFT逆运算。由于过采样过程将采样间隔缩小1倍，Wiener滤波器M′_W[u，v]中的占空比应提高1倍，即α_OS＝2α，β_OS＝2β。 

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。 

Claims (1)

1.一种基于泰勒展开积分退化模型的微扫描图像重构方法，其特征在于，

第一步：利用基于离散傅里叶变换DFT的亚像元图像配准算法，获得多帧微扫描图像的帧间位移量；

第二步：在步骤一获得的帧间位移量的基础上，利用基于帧间差分过采样技术，将多帧微扫描图像重建出一幅过采样图像g_os[m，n]；

第三步：建立焦平面探测器的积分退化模型，进行泰勒(Taylor)展开，得到积分退化函数H[u，v]，用其构建Wiener滤波器M_W[u，v]，经直流分量归一化之后得到归一化后Wiener滤波器M′_W[u，v]，用其对第二步得到的过采样图像g_os[m，n]进行Wiener滤波，复原出重构图像

完成微扫描图像的重构；其中[u，v]为频域坐标，u为频域上的水平坐标，v为频域上的竖直坐标；

其中得到积分退化函数H[u，v]的过程为：

假设f(x，y)为经过光学系统成像在焦平面探测器上的模拟图像，g[m，n]为经焦平面探测器积分和采样后得到的退化图像，则

$\begin{array}{c}g[m,n]=\frac{1}{\mathrm{pq}}{\∫}_{\mathrm{nb}-q/2}^{\mathrm{nb}+q/2}{\∫}_{\mathrm{ma}-p/2}^{\mathrm{ma}+p/2}f(x,y)\mathrm{dxdy}\\ =\frac{1}{\mathrm{pq}}{\∫}_{-q/2}^{q/2}{\∫}_{-p/2}^{p/2}f(x+\mathrm{ma},y+\mathrm{nb})\mathrm{dxdy}\end{array}m=\mathrm{1,2},...,M,n=\mathrm{1,2},...,N---\left(1\right)$

式中，a和b为焦平面探测器的水平和竖直像元间隔，p为像元在水平方向的尺寸，q为像元在竖直方向的尺寸，M×N为焦平面探测器像面规格，M为焦平面探测器像面在水平方向的尺寸，N为焦平面探测器像面在竖直方向的尺寸，(x，y)是模拟图像坐标系中的坐标，[m，n]是离散图像坐标系中的坐标，m为离散图像坐标系中的水平坐标，最大取值范围为焦平面探测器像面在水平方向的尺寸，n为离散图像坐标系中的竖直坐标，最大取值范围为焦平面探测器像面在竖直方向的尺寸；

对式(1)中f(x+ma，y+nb)关于f[m，n]进2阶泰勒展开得

$f(x+\mathrm{ma},y+\mathrm{nb})\≅f[m,n]+f_m[m,n]x+f_n[m,n]y$

$+\frac{f_{m^2}[m,n]x^2+f_{n^2}[m,n]y^2+2f_{\mathrm{mn}}[m,n]\mathrm{xy}}{2}---\left(2\right)$ 将式(2)代入式(1)得

$g[m,n]\≅f[m,n]+\frac{p^2}{24}{f}_{m^2}[m,n]+\frac{q^2}{24}{f}_{n^2}[m,n]---\left(3\right)$

式中，

和

分别表示图像f[m，n]在m和n方向上的2阶偏导数，f_m[m，n]和f_n[m，n]分别表示图像f[m，n]在m和n方向上的1阶偏导数，f_mn[m，n]表示图像f[m，n]先在m方向上取1阶偏导数再在n方向上取1阶偏导数；为了抑制噪声，在数值上用式(4)近似

和

$\left\{\begin{array}{c}{f}_{m^2}[m,n]=\left\{f\right[m-1,n-1]-2f[m,n-1]+f[m+1,n-1]\\ +2\left(f\right[m-1,n]-2f[m,n]+f[m+1,n\left]\right)\\ +f[m-1,n+1]-2f[m,n+1]+f[m+1,n+1]\}/4\\ f_{n^2}[m,n]=\{\left(f\right[m-1,n-1]-2f[m-1,n]+f[m-1,n+1\left]\right)\\ +2\left(f\right[m,n-1]-2f[m,n]+f[m,n+1\left]\right)\\ +\left(f\right[m+1,n-1]-2f[m+1,n]+f[m+1,n+1\left]\right)\}/4\end{array}\right.---\left(4\right)$

将式(4)代入式(3)，再进行DFT得

$G[u,v]=F[u,v][1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)]---\left(5\right)$

式中，F[u，v]和G[u，v]分别表示f[m，n]和g[m，n]的DFT，α＝p/a和β＝q/b分别表示焦平面探测器在水平和竖直方向的占空比；

因此，焦平面探测器的积分退化函数为H[u，v]

$H[u,v]=1-\frac{{\mathrm{\α}}^2}{6}{\sin }^2\left(\frac{u}{M}\mathrm{\π}\right){\cos }^2\left(\frac{v}{N}\mathrm{\π}\right)-\frac{{\mathrm{\β}}^2}{6}{\sin }^2\left(\frac{v}{N}\mathrm{\π}\right){\cos }^2\left(\frac{u}{M}\mathrm{\π}\right)---\left(6\right)$

Wiener滤波过程为：

基于积分退分函数H[u，v]，构建Wiener滤波器M_W[u，v]

$M_W[u,v]=\frac{H^{*}[u,v]}{H^2[u,v]+\mathrm{\Γ}}---\left(7\right)$

式中H^*[u，v]为H[u，v]的复共轭；Γ为设定的常数；将M_W[u，v]用直流分量归一化后得到归一化后Wiener滤波器M′_W[u，v]

$M_W^{\′}[u,v]=\frac{M_W[u,v]}{M_W\left[\mathrm{0,0}\right]}---\left(8\right)$

式中，M_W[0，0]是M_W[u，v]在u＝0，v＝0情况下的值，即M_W[u，v]的直流分量；

步骤303、利用式(8)描述的Wiener滤波器M′_W[u，v]，对第二步得到的过采样图像g_os[m，n]进行滤波，获得具有重建图像

式中，

表示DFT运算，

表示DFT逆运算。