CN101540043B - 单幅图像复原的解析迭代快速频谱外推方法 - Google Patents

Abstract

本发明公开了一种单幅图像复原的解析迭代快速频谱外推方法，包括迭代系统初始化过程和解析迭代过程。其中：初始化过程计算降质图像频谱、设置初始迭代解、系统初始参数；解析迭代过程通过不断精细的梯度信息阈值收缩估计图像的水平和垂直方向梯度信息，通过解析迭代频谱和梯度信息频谱解析预测得到预测频谱，通过图像频谱合成校正将降质图像频谱和预测频谱解析合成校正得到高分辨率图像。本发明方法基于快速傅立叶变换技术，通过截止频率以下频谱，外推截止频率以上频谱，以很小的复杂度达到图像快速复原。去模糊和抑制噪声能力明显优于常规方法，信噪比得到显著提高。该方法非常便于利用快速傅立叶变换DSP芯片组成的图像处理硬件实现。

Description

单幅图像复原的解析迭代快速频谱外推方法

技术领域

本发明属于针对空间不变成像系统单幅退化图像复原技术，特别是针对光学遥感单幅图像复原的解析迭代频谱外推方法。 

背景技术

全天候、全天时高分辨率遥感成像理论和技术在全世界范围内受到极大重视和飞跃发展，成为空间对地观测发展的热点。未来遥感探测成像器将进一步朝着高分辨率、高信噪比、轻量化方向发展。遥感成像应用效能(如目标识别)主要取决于获取的遥感图像的视觉质量和空间分辨率，而空间分辨率是衡量遥感图像的重要指标。另外，为了获得高分辨率图像，目前大多数高分辨率对地观测相机系统的有效载荷体积和重量庞大，因此在对地观测和深空探测领域对遥感探测成像器相机的小型化和轻量化提出了非常高的要求。 

理论上，光学遥感图像的高分辨率，依赖于CCD阵列单元的尺寸、CCD像元个数、光学口径等，但实际获取图像的分辨率还会受到系统光学模糊、运动模糊、系统噪声等多种误差因素的影响。提高光学遥感图像的清晰度和分辨率通常有两种途径：一是改进和更新硬件设备；二是通过建立数学和物理模型，利用数据处理技术提高图像的清晰度和分辨率。 

一般遥感成像系统的降质模型与系统点扩散函数h相关。记F(h)为h的傅立叶变换，并令 $(\mathrm{\ξ},\mathrm{\η})\∈[-\frac{1}{2},\frac{1}{2}]$ 为频域坐标，则F(h)称为调制传递函数(MTF)，它一般包含如下3部分。 

1)CCD的调制传递函数。在焦平面上CCD敏感像素是一个细小的[-c/2，c/2]²的窗口，其中C为相邻CCD像元的距离，可以用一个矩形函数来描述，这样可推导出CCD的MTF为： 

$H_s(\mathrm{\ξ},\mathrm{\η})=\frac{\sin \left(\mathrm{\π\ξc}\right)}{\mathrm{\π\ξc}}\·\frac{\sin \left(\mathrm{\π\ηc}\right)}{\mathrm{\π\ηc}}{e}^{-2\mathrm{\π}{\mathrm{\β}}_{1}c\left|\mathrm{\ξ}\right|}{e}^{-2\mathrm{\π}{\mathrm{\β}}_{2}c\left|\mathrm{\η}\right|}$

2)光学孔径系统的调制传递函数H_o。可看作是一个各向同性的低通滤波器，表示为 

$H_o(\mathrm{\ξ},\mathrm{\η})=e^{-2\mathrm{\π\αc}\sqrt{{\mathrm{\ξ}}^2+{\mathrm{\η}}^2}}$

3)运动引起的调制传递函数H_m。由于CCD敏感单元会在一定的采样时间内对到达的光子数目计数。而在采样时间假设系统运动距离为τ，则将产生运动模糊，其MTF为： 

$H_m(\mathrm{\&xi;},\mathrm{\&eta;})=\frac{\sin (\mathrm{\&pi;}<(\mathrm{\&xi;},\mathrm{\&eta;}),(d_1,d_2)>\mathrm{\&tau;})}{(\mathrm{\&pi;}<(\mathrm{\&xi;},\mathrm{\&eta;}),(d_1,d_2)>\mathrm{\&tau;})}$

一般而言，遥感成像光学系统可假设为圆对称系统，则总体MTF是各部分MTF的乘积。 

F(h)＝H_soH_ooH_m。 

假设u(x，y)为未降质的图像，从而降质观测图像u₀(x，y)的成像模型为： 

u₀(x，y)＝F^-1{F(u(x，y))oF(h)}_(x，y)+n    (1) 

或 $u_0(x,y)=u(x,y)\&CircleTimes;h(x,y)+n$

式中，(x，y)为空间坐标，(ξ，η)为频率域坐标，F和F^-1分别表示傅立叶正反变换，n表示噪声， 表示卷积算子。 

图像复原问题实质上从u₀中估计合适的u。Wiener滤波复原算法源于最小二乘估计，通过假定图像和噪声都是属于广义平稳随机过程，得到一种简单的频域复原算法，公式为： 

$u=F^{-1}\left(\frac{F{\left(h\right)}^{*}\mathrm{oF}\left(u_0\right)}{F{\left(h\right)}^{*}\mathrm{oF}\left(h\right)+S_n\left(n\right)/S_u\left(u\right)}\right)---\left(2\right)$

其中F(u)，F(h)为u和h的傅里叶变换，s_n(n)，s_u(u)分别为噪声和未失真图像的功率谱，“o”表示矩阵的元素相乘，而“*”表示复共轭。功率谱比值s_n(n)/s_u(u)起到了规整化作用。但两个功率谱常常难以估计，因此利用近似Wiener滤波 

$u=F^{-1}\left(\frac{F{\left(h\right)}^{*}\mathrm{oF}\left(u_0\right)}{F{\left(h\right)}^{*}\mathrm{oF}\left(h\right)+\mathrm{\&gamma;}}\right)---\left(3\right)$

式中γ是一个正数，在数值上最好为观测图像信噪比的倒数。 

受限制最小二乘复原方法(RLS)，通过Laplace算子生成的循环矩阵的规整化，得到了如下滤波公式 

$u=F^{-1}\left(\frac{F{\left(h\right)}^{*}\mathrm{oF}\left(u_0\right)}{F{\left(h\right)}^{*}\mathrm{oF}\left(h\right)+\mathrm{\&alpha;}{F}^{*}\left(C\right)\mathrm{oF}\left(C\right)}\right)---\left(4\right)$

其中C为Laplace算子生成的循环矩阵，F(C)为C的傅里叶变换。Wiener滤波算法简单、便于实现、在噪声强度不大(高信噪比)情形能够获得较好的复原效果。但随着噪声强度的增大，复原效果急剧下降并存在较多的“寄生波纹效应”；RLS具有Wiener滤波类似的频域公式，实现也简单，能够减少“寄生波纹效应”，但是容易出现过度平滑现象，不能很好的保持图像边缘。 

Lucy-Richardson算法是一种基于泊松分布的最大似然图像复原算法。LR算法可表示为 

$u^{k+1}(i,j)=u^k(i,j)(\frac{u_0(i,j)}{(h\&CircleTimes;u^k)(i,j)}\&CircleTimes;h(i,j))---\left(5\right)$

从以上分析可看出，该算法的超分辨能力来自于每步迭代中涉及的非线性处理，而且该算法的更新规则自动实现了有约束条件的频谱外推，可用于图像复原。Lucy-Richardson算法是一种乘积形式的迭代恢复方法，最初主要运用于医学和天文图像复原。由于噪声的存在，算法虽然能够外推部分高频分量，但同时也破坏了图像通带内的低频分量，出现了一些寄生波纹，图像复原能力有限，而且收敛速度较慢。 

近几年，随着图像处理中的变分偏微分方程的深入研究，出现了一系列新方法。图像复原属于Hardmard意义下的非适定数学反问题。基于Rudin-Osher连续全变差模型提出的图像正则化复原算法(1.T.Chan，S.Esedoglu，F.Park，A.Yip，Recent developmentsin total variation image restoration，in：Handbook of Mathematical Models in ComputerVision，Springer-Verlag，2005，pp.17-30.)是目前应用广泛的算法。该算法基于如下变分模型： 

$\underset{u}{\min }{\left|\right|u_0(x,y)-u(x,y)\&CircleTimes;h(x,y)\left|\right|}_2^2+\mathrm{\&alpha;}{\left|\right|u\left|\right|}_{\mathrm{TV}}---\left(6\right)$

其中连续全变差模型定义为 

${\left|\right|u\left|\right|}_{\mathrm{TV}}=\underset{\mathrm{\&Omega;}}{\&Integral;\&Integral;}|\&dtri;u|\mathrm{dxdy}=\underset{\mathrm{\&Omega;}}{\&Integral;\&Integral;}\sqrt{{\left(u_x\right)}^2+{\left(u_y\right)}^2}\mathrm{dxdy}$

利用变分法和梯度下降法(GD)转换为求解非线性方程 

$\frac{\&PartialD;u}{\&PartialD;t}=H^{*}(u_0-\mathrm{Hu})+\mathrm{\&alpha;}\mathrm{div}\left(\frac{\&dtri;u}{|\&dtri;u|}\right)---\left(7\right)$

在(7)中H为h对应的卷积算子，为避免 

等于零，需要引入参数ε＞0，将其替换为 

全变差图像复原模型能够减少“寄生波纹效应”也能较好保持图像边缘，但是会导致图像中出现“阶梯效应”。另外梯度下降算法存在如下缺点：1)收敛速度较慢；2)在求解过程中需要引入参数ε＞0，而且因为参数ε的引入迭代解不能严格收敛于(6)的精确解；3)每次迭代都要计算一次H^*H，由于线性空间不变的模糊算子的支集一般很大，当受到模糊影响的遥感图像尺寸很大时去模糊的计算代价更高。其他变分偏微分方程的图像复原方法往往存在类似的问题，大大限制了其在工程应用中的实用性。 

频谱外推方法也是图像复原的一条有效途径，该方法主要源于解析延拓理论：1)任何空域有界函数的傅立叶变换是解析函数；2)对任何解析函数，只要能准确知道它在有限区间内的部分信息，就可以唯一的确定整个函数。对一幅图像而言，由于其空域有界，因此其频谱函数必然解析。根据解析延拓理论，截止频率以上的频谱信息可以通过截止频率以下的频谱得以重建，从而能够实现图像的高分辨率复原。例如根据此思想提出的一种针对无源毫米波成像运用的最大似然频域校正超分辨算法。该算法以Lucy-Richardson算法作为主迭代，通过一种频域校正算法，用Wiener滤波器恢复的频谱代替通带内的频谱，在外推图像高频分量的同时，保证图像的低频分量不被破坏，对无源毫米波成像复原获得了应用(2.郑鑫，杨建宇，李良超.无源毫米波成像的最大似然频域校正超分辨算法.自动化学报，2009，Vol 35，N0.1，PP 28-33)。但是由于Wiener滤波和Lucy-Richardson算法本身存在“寄生波纹效应”明显和边缘清晰度不足等问题，直接将此两种算法耦合外推对与复杂的遥感图像效果处理并不理想。 

发明内容

本发明的目的在于提供一种单幅图像复原的解析迭代快速频谱外推方法，以具有解析计算的迭代频谱外推方式，保证图像通带内频谱和高频频谱有效合成，克服传统方法“寄生波纹效应”明显和边缘清晰度不足等问题，而且算法充分利用快速傅立叶变换技术，实现图像的快速复原。 

实现本发明目的的技术解决方案为：一种单幅图像复原的解析迭代快速频谱外推方法，包括迭代系统初始化过程和解析迭代过程，即 

1.1所述的迭代系统初始化过程包括：输入一幅待复原的M×N大小的降质图像u₀，计算降质图像频谱F(u⁰)；设置初始迭代解u⁽⁰⁾＝u₀，初始频谱校正图像w⁽⁰⁾＝u₀；设置迭代系统参数初始值：梯度信息阈值α⁽⁰⁾＞0，预测贡献度控制参数β⁽⁰⁾＞0，合成贡献度参数γ⁽⁰⁾＞1，图像频谱误差预设阈值ε≤10^-4；并设置空间不变成像系统的系统调制传递函数F(h)；

1.2 所述的解析迭代过程包括如下步骤： 

步骤1：首先利用第k-1次的频谱校正图像w^(k-1)，其中k≥1，分别计算其在水平方向和垂直方向的向前差分，并利用第k-1次的梯度信息阈值α^(k-1)作阈值收缩判断，得到第k次的水平方向差分估计 

和垂直方向差分估计 

步骤2：将第k次的水平方向差分估计 

作快速傅立叶变换(FFT)得到水平方向高频频谱 

同样第k次的垂直方向差分估计 

作快速傅立叶变换(FFT)得到垂直方向高频频谱 

然后，利用计算得到的水平方向高频频谱 

垂直方向高频频谱 以及第k-1次合成复原图像频谱F(u^(k-1))进行频谱预测得到第k次的预测频谱F(w^(k))； 

步骤3：利用图像频谱合成校正，将第k次的预测频谱F(w^k)和降质图像频谱F(u⁰)结合光学系统调制传递函数F(h)及其复共轭F(h)^*通过矩阵元素逐点相乘、相加和矩阵元素逐点相除的方式进行合成校正，得到合成校正频谱F(u^(k))； 

步骤4：采用2范数的平方 

计算第k次校正图像频谱F(u^(k))和第K-1次校正图像频谱F(u^(k-1))的误差，如果 

则转至步骤5，否则转至步骤6； 

步骤5：表明频谱外推没有收敛，则进行参数自适应校正，转至步骤1，k＝k+1； 

步骤6：判断是否达到收敛条件 如果到达则算法结束，输出复原图像，对合成校正图像F(u^(k))进行快速逆傅立叶变换得到复原图像。 

本发明与现有技术相比，其显著优点：(1)本发明能够实现高视觉质量的快速复原，能够很好的克服传统Wiener滤波算法的“寄生波纹效应”和约束最小二乘算法的“过度平滑现象”，边缘清晰度高。经大量统计实验结果表明，而本发明方法一般在1.56秒时间内，能够取得峰值信噪比(PSNR)比传统算法高1.0-2.0dB左右，改进信噪比高0.5-0.8dB左右，在抗噪性能、去模糊效果、实时性等三方面得到很好的平衡。(2)本发明方法在光学遥感成像系统、合成孔径雷达成像系统、电视制导、航空探测等都有广泛的应用前景，同时也为超分辨复原提供了新的途径。 

下面结合附图对本发明作进一步详细描述。 

附图说明

图1是本发明的整体迭代结构的流程图。 

图2是本发明解析迭代过程步骤1数据流程图。 

图3是本发明解析迭代过程步骤2数据流程图。 

图4是本发明解析迭代过程步骤3数据流程图。 

图5是本发明的实验测试图像和MTF曲线。 

图6是本发明与传统方法的一个复原实验对比效果图。 

图7是本发明与传统方法的一个复原实验中的边缘清晰度对比图。 

具体实施方式

本发明以成像系统光学调制传递函数、运动模糊函数为基础，通过不断精细的估计梯度信息、预测频谱和合成低高频频谱的方式进行频谱迭代外推，达到降质图像反演清晰图像的目的，基本思想源于解析延拓理论：1)任何空域有界函数的傅立叶变换是解析函数；2)对任何解析函数，只要能准确知道它在有限区间内的部分信息，就可以唯一的确定整个函数。对一幅图像而言，由于其空域有界，因此其频谱函数必然解析。根据解析延拓理论，截止频率以上的频谱信息可以通过截止频率以下的频谱得以重建，从而能够实现图像的高分辨率复原。 

结合图1，为实现空间不变成像系统降质图像高视觉质量的快速图像复原，本发明基于解析迭代频谱外推的方法包括迭代系统初始化过程和解析迭代过程。 

1.1所述的迭代系统初始化过程包括：①输入一幅待复原的M×N大小的降质图像u₀，计算降质图像频谱F(u⁰)；②初始迭代解u⁽⁰⁾＝u₀，初始频谱校正图像w⁽⁰⁾＝u₀；③设置迭代系统参数初始化值：梯度信息阈值α⁽⁰⁾＞0，预测贡献度控制参数β⁽⁰⁾＞0，合成贡献度参数γ⁽⁰⁾＞1，图像频谱误差预设阈值ε≤10^-4；④空间不变成像系统的 系统调制传递函数F(h)，其中F(h)为各部分MTF的乘积F(h)＝H_soH_ooH_m。 

1.2所述的解析迭代过程步骤如下 

步骤1：首先利用第k-1次的频谱预测图像w^(k-1)(k≥1)，分别计算其在水平方向和垂直方向的向前差分，并利用第k-1次的梯度信息阈值α^(k-1)作阈值收缩判断，得到第k次的水平方向和垂直方向信息估计。 

如图2所示该步骤的详细过程为，首先通过具有周期边界条件的差分算子计算第k-1次的频谱校正图像w^(k-1)在空间位置(i，j)的水平和垂直方向差分：当1＜i＜N，1＜j＜M时，计算第k-1次频谱预测图像w^(k-1)在空间位置(i，j)处的水平方向的差分  ${\left(D_1{w}^{(k-1)}\right)}_{i,j}=w_{i+1,j}^{(k-1)}-w_{i,j}^{(k-1)}$ 和垂直方向的差分 ${\left(D_2{w}^{(k-1)}\right)}_{i,j}=w_{i,j+1}^{(k-1)}-w_{i,j}^{(k-1)};$ 在i＝N，j＝M时的图像边界需采取周期边界条件计算水平方向差分 ${\left(D_1{w}^{(k-1)}\right)}_{N,j}=w_{1,j}^{(k-1)}-w_{N,j}^{(k-1)}$ 和垂直方向差分  ${\left(D_2{w}^{(k-1)}\right)}_{i,M}=w_{i,1}^{(k-1)}-w_{i,M}^{(k-1)}.$

然后计算第k-1次的频谱校正图像w^(k-1)在空间位置(i，j)的梯度模  $\left|\right|{\left(D{w}^{(k-1)}\right)}_{i,j}\left|\right|=\sqrt{{\left({\left(D_1{w}^{(k-1)}\right)}_{i,j}\right)}^2+{\left({\left(D_2{w}^{(k-1)}\right)}_{i,j}\right)}^2},$

最后利用第k-1次的梯度信息阈值α^(k-1)进行梯度阈值收缩得到水平方向和垂直方向的差分信息估计，得到第k次的水平方向差分估计(v₁)_i，j ^(k)和垂直方向差分估计(v₂)^i，j(k)。估计方法为： 

水平方向： ${\left(v_1\right)}_{i,j}^{\left(k\right)}=\max \left(\right|\left|{\left(D{w}^{(k-1)}\right)}_{i,j}\right||-{\mathrm{\&alpha;}}^{(k-1)},0)\&CenterDot;\frac{{\left(D_1{w}^{(k-1)}\right)}_{i,j}}{\left|\right|{\left(D{w}^{(k-1)}\right)}_{i,j}\left|\right|}$

垂直方向： ${\left(v_2\right)}_{i,j}^{\left(k\right)}=\max \left(\right|\left|{\left(D{w}^{(k-1)}\right)}_{i,j}\right||-{\mathrm{\&alpha;}}^{(k-1)},0)\&CenterDot;\frac{{\left(D_2{w}^{(k-1)}\right)}_{i,j}}{\left|\right|{\left(D{w}^{(k-1)}\right)}_{i,j}\left|\right|}$

步骤2：利用高频频谱预测得到第k次的预测频谱。如图3所示该步骤的详细过程为：将第k次的水平方向差分估计(v₁)_i，j ^(k)和垂直方向差分估计(v₂)^i，j(k)作快速傅立叶变换(FFT)得到水平方向高频频谱F(v₁ ^(k))和垂直方向高频频谱F(v₂ ^(k))；将水平方向差 分算子[-1，1]卷积核D₁和垂直方向差分算子 $\left[\begin{array}{c}-1,\\ 1\end{array}\right]$ 卷积核D₂作FFT；结合第k-1次复原图像频谱F(u^(k-1))，通过矩阵元素逐点相加、相乘和矩阵元素逐点相除的方式作预测更新得到第k次的预测频谱F(w^(k))。计算公式为： 

$F\left(w^{\left(k\right)}\right)=\{F\left(u^{(k-1)}\right)+{\mathrm{\&beta;}}^{(k-1)}(F{\left(D_1\right)}^{*}\mathrm{oF}\left(v_1^{\left(k\right)}\right)+F{\left(D_2\right)}^{*}\mathrm{oF}\left(v_2^{\left(k\right)}\right))\}/\{1+{\mathrm{\&beta;}}^{(k-1)}(F{\left(D_1\right)}^{*}\mathrm{oF}\left(D_1\right)+F{\left(D_2\right)}^{*}\mathrm{oF}\left(D_2\right))\}.$

其中：F(D₁)，F(D₂)分别为水平和垂直差分算子线性卷积核的快速傅立叶变换；β^k为正参数，称为预测贡献度控制参数，控制水平方向高频频谱和垂直方向高频频谱对预测频谱的贡献度。“o”表示矩阵元素逐点相乘，“/”表示矩阵元素逐点相除。 

步骤3：利用图像频谱合成校正。如图4所示，该步骤的详细过程为：将第k次的预测频谱F(w^k)和降质图像频谱F(u⁰)结合光学系统调制传递函数F(h)及其复共轭F(h)^*通过通过矩阵元素逐点相乘、相加和矩阵元素逐点相除的方式进行合成校正，得到合成校正频谱F(u^(k))。计算公式为： 

F(u^(k))＝{F(h)^*oF(u₀)+γ^(k-1)F(W^(k))}/{F(h)^*oF(h)+γ^(k-1)} 

其中：F(h)为光学系统调制传递函数，其可以包括CCD的调制传递函数、光学孔径系统的调制传递函数、运动引起的调制传递函数等多个传递函数的乘积形成的系统调制传递函数，F(h)^*为F(h)的复共轭。γ^(k)为取正值的参数，称为合成贡献度参数。同样“o”表示矩阵元素逐点相乘，“/”表示矩阵元素逐点相除。 

步骤4：采用2范数的平方||·||₂ ²计算第k次校正图像频谱F(u^(k))和第K-1次校正图像频谱F(u^(k-1))的误差，如果 ${\left|\right|F\left(u^{\left(k\right)}\right)-F\left(u^{(k-1)}\right)\left|\right|}_2^2/{\left|\right|F\left(u^{\left(k\right)}\right)\left|\right|}_2^2>\mathrm{\&epsiv;},$ 则转至步骤5，否则转至步骤6。 

步骤5：表明频谱外推没有收敛，则进行参数自适应校正，自适应调整参数，转至步骤1，k＝k+1。参数调整方法如下： 

1.梯度信息阈值α^(k-1)，在迭代过程中逐渐减少，本发明方法中经实验证实α⁰取特定的参数值10≤α⁽⁰⁾≤125，按照特定的递减顺序α^(k)/α^(k-1)＜1进行控制即满足α⁽⁰⁾＞α⁽¹⁾＞...α^(k-1)＞...＞0即可。 

2.本发明方法经研究实验测得贡献度参数β^(k)与需与α^(k)满足正比关系，β^(k)＝k·α^(k)，k∈(0.2，0.8)为实验参数。 

3.合成贡献度参数γ^(k)，在迭代过程中按照递增顺序进行控制γ^(k)/γ^(k-1)＞2且γ⁽⁰⁾∈[1，2 0]。 

步骤6：判断是否达到收敛条件 ${\left|\right|F\left(u^{\left(k\right)}\right)-F\left(u^{(k-1)}\right)\left|\right|}_2^2/{\left|\right|F\left(u^{\left(k\right)}\right)\left|\right|}_2^2\&le;\mathrm{\&epsiv;},$ 如果到达则算法结束，输出复原图像，对合成校正图像F(u^(k))进行快速逆傅立叶变换得到复原图像； 

u^(k)＝F^-1(F(u^(k))) 

步骤1、2、3是本发明的核心处理过程，均参与迭代更新过程。每个核心单元均具有解析形式的数学公式控制。步骤1在图像空域中进行，而步骤2和步骤3均在傅立叶域进行，整个频谱一次外推过程包含二次快速傅立叶正变换和一次快速傅立叶逆变换的算法复杂度。算法的迭代执行速度快、对一幅M×M大小的图像整体复杂度为O(M² log M)。 

下面结合图5、图6和图7所示，通过一个图像复原实施例及其效果评价来说明本发明方法的技术效果。 

在实施例中，如图5.a所示，实验图像为一幅清晰的370×370的光学遥感图像(参考图像)。如图5.b所示为实验中遥感成像系统的MTF模拟曲线。 

本发明方法中的迭代参数α⁽⁰⁾＝100，β⁽⁰⁾＝30，γ⁽⁰⁾＝10，迭代终止条件均取为  ${\left|\right|F\left(u^{\left(k\right)}\right)-F\left(u^{(k-1)}\right)\left|\right|}_2^2/{\left|\right|F\left(u^{\left(k\right)}\right)\left|\right|}_2^2\&le;{10}^4.$

实验中对比算法包括：MATLAB 7.1中Wiener滤波、约束最小二乘算法(RLS)，Lucy-Richardson算法；同时我们和全变差梯度最速下降法(TVGD)进行比较。TVGD采取中心差分的离散迭代格式 

$u^{n+1}=u^n+\mathrm{\&gamma;}\&CenterDot;[H^{*}(u_0-H{u}^{\left(n\right)})+\mathrm{\&alpha;}\mathrm{div}\left(\frac{\&dtri;u^{\left(n\right)}}{\sqrt{{|\&dtri;u^{\left(n\right)}|}^2+\mathrm{\&epsiv;}}}\right)]$

其中参数γ取为0.1，网格比τ取10^-3，Lagrange参数α取25.0。 

实施例在MATLAB 7.1平台模拟仿真实现，计算环境为Dell D620笔记本电脑、Intel1.83GHz CPU，1G内存。为了验证本发明方法的有效性。用于定量评估去模糊算法的性能指标采取改进信噪比(Δ_SNR)、峰值信噪比(PSNR)和相对误差(ReErr)。设u₀，u，u^*分别为M×N的降质图像、参考图像和去模糊图像，Δ_SNR、PSNR和ReErr定义如下： 

${\mathrm{\&Delta;}}_{\mathrm{SNR}}=10\log 10\{{\left|\right|u-u_0\left|\right|}_2^2/{\left|\right|u-u^{*}\left|\right|}_2^2\},$

PSNR＝10 log 10{M×N×255²/||u-u^*||²}， 

$\mathrm{ReErr}={\left|\right|u-u^{*}\left|\right|}_2^2/{\left|\right|u\left|\right|}_2^2$

一般而言改进信噪比Δ_SNR和峰值信噪PSNR比越大越好，而相对误差ReErr越小越好。 

图6和表1分别给出了不同算法和本发明方法的复原结果图像和各项性能指标。 

如图6.a所示为根据降质模型模拟生成的降质图像，其中加入的噪声为均值为0、方差为15高斯噪声。可看到降质图像模糊效应和噪声干扰非常严重，导致空间分辨率非常低，严重影响人们的视觉理解。如图6.b所示为Wiener滤波结果。Wiener滤波复原图像“寄生波纹效应”明显，复原效果较差。如图6.c所示为约束最小二乘(RLS)算法的滤波结果，如图6.d所示为Lucy-Richardson算法的滤波结果，可看到两种算法对于“寄生波纹效应”有所改善，但“阶跃”边缘和线状边缘保持效果不好(例如图6.c，图6.d中江面和城市的交界线)。如图6.e所示为TVGD算法结果，该算法能够极大的减少“寄生波纹效应”，视觉质量得到较好提高，但存在视觉上的“阶梯效应”(如图6.d)。如图6.f所示为本发明方法的复原图像，本发明方法复原图像很好的克服了“寄生波纹效应”，并且复原图像的细节在对比算法中最清晰，视觉质量最高。 

表一给出了本文发明方法和对比算法的性能指标情况。从复原图像质量客观评价指标来看，本发明方法复原图像质量的客观评价指数得到进一步改善。例如PSNR比TVAP提高1dB，而Δ_SNR比TVAP提高0.8dB。 

从算法所花费的CPU时间来看，Wiener滤波和RLS算法是非迭代的算法，时间最少，分别为0.14秒和0.31秒。Lucy-Richardson算法、TVGD算法和本发明方法都属于迭代算法，时间较长。其中TVGD算法的所用时间最长，达到84.281秒，说明算法收敛速度较慢；Lucy-Richardson算法花费3.08秒，而本发明方法仅用了1.56秒就达到优于TVGD算法的图像质量。 

表1：不同复原算法性能指标比较结果 

Method	PSNR	ReErr	ΔSNR	Time(s)
					降质图像	13.2716	0.009784	1.97	-
Wiener滤波	20.6431	0.008104	2.41	0.14
					RLS算法	21.1897	0.006456	3.78	0.31
Lucy-Richardson算法	20.2462	0.007314	3.23	3.08
					TVGD	22.0017	0.006124	3.81	84.281
本发明方法	23.4100	0.005125	4.84	1.56

通过对复原图像边缘细节点集比对情况能够很好的区分复原算法的有效性。如图7所示为通过Matlab中的Canny算子提取原始图像、降质图像和各种算法输出图像的边缘图进行分析。图7.a为参考图像的清晰边缘图，图7.b为降质图像的边缘图。图7.c为Wiener滤波图像边缘图，不难发现Wiener滤波图像边缘存在较多的“寄生波纹效应”引起的伪边缘；图7.d所示为RLS算法复原图像边缘，图7.e所示为Lucy-Richardson算法复原图像边缘，相比Wiener滤波算法，这两种算法伪边缘得到减少，但是也丧失了较多的边缘细节。图7.f为全变差梯度最速下降法(TVGD)复原图像边缘图，该算法边缘清晰度得到很大提高。图7.g为本发明方法的复原图像边缘图，进行综合对比，本发明方法对图像边缘的保持效果最好，细节得到很好的复原。因此本发明方法得到了很好的图像质量-计算时间性价比，实用性非常强。 

Claims (8)

1.一种单幅图像复原的解析迭代快速频谱外推方法，其特征在于包括迭代系统初始化过程和解析迭代过程，即

1.1 所述的迭代系统初始化过程包括：输入一幅待复原的M×N大小的降质图像u₀，计算降质图像频谱F(u⁰)；设置初始迭代解u⁽⁰⁾＝u₀，初始频谱校正图像w⁽⁰⁾＝u₀；设置迭代系统参数初始值：梯度信息阈值α⁽⁰⁾＞0，预测贡献度控制参数β⁽⁰⁾＞0，合成贡献度参数γ⁽⁰⁾＞1，图像频谱误差预设阈值ε≤10^-4；并设置空间不变成像系统的系统调制传递函数F(h)；

1.2 所述的解析迭代过程包括如下步骤：

步骤1：首先利用第k-1次的频谱校正图像w^(k-1)，其中k≥1，分别计算其在水平方向和垂直方向的向前差分，并利用第k-1次的梯度信息阈值α^(k-1)作阈值收缩判断，得到第k次的水平方向差分估计

和垂直方向差分估计

步骤2：将第k次的水平方向差分估计

作快速傅立叶变换(FFT)得到水平方向高频频谱

同样第k次的垂直方向差分估计

作快速傅立叶变换(FFT)得到垂直方向高频频谱

然后，利用计算得到的水平方向高频频谱

垂直方向高频频谱以及第k-1次合成复原图像频谱F(u^(k-1))进行频谱预测得到第k次的预测频谱F(w^(k))；

步骤3：利用图像频谱合成校正，将第k次的预测频谱F(w^k)和降质图像频谱F(u⁰)结合光学系统调制传递函数F(h)及其复共轭F(h)^*通过矩阵元素逐点相乘、相加和矩阵元素逐点相除的方式进行合成校正，得到合成校正频谱F(u^(k))；

步骤4：采用2范数的平方

计算第k次校正图像频谱F(u^(k))和第K-1次校正图像频谱F(u^(k-1))的误差，如果

则转至步骤5，否则转至步骤6；

步骤5：表明频谱外推没有收敛，则进行参数自适应校正，调整系统参数，转至步骤1，k＝k+1；

步骤6：判断是否达到收敛条件

如果到达则算法结束，输出复原图像，对合成校正图像F(u^(k))进行快速逆傅立叶变换得到复原图像。

2.根据权利要求1所述的单幅图像复原的解析迭代快速频谱外推方法，其特征在于解析迭代过程步骤1中梯度信息阈值收缩满足如下关系：

水平方向：

垂直方向：

其中梯度算子D采取向前或向后差分算子，如

${\left(\mathrm{Dw}\right)}_{i,j}:=\left(\begin{array}{c}{\left(D_{1}w\right)}_{i,j}\\ {\left(D_{2}w\right)}_{i,j}\end{array}\right):=\left(\begin{array}{c}{w}_{i+1,j}-w_{i,j}\\ w_{i,j+1}-w_{i,j}\end{array}\right),1<i<N,1<j<M,$

这里i＝N，j＝M时采取周期边界条件。

3.根据权利要求1所述的单幅图像复原的解析迭代快速频谱外推方法，其特征在于解析迭代过程的步骤2中频谱预测得到第k次的预测频谱，是将第k次的水平方向差分估计和垂直方向差分估计

作快速傅立叶变换(FFT)得到水平方向高频频谱

和垂直方向高频频谱

将水平方向差分算子[-1，1]卷积核D₁和垂直方向差分算子卷积核D₂作FFT；并将第k-1次复原图像频谱F(u^(k-1))通过矩阵元素逐点相加、相乘和矩阵元素逐点相除的方式作预测更新得到第k次的预测频谱F(w^(k))，计算关系为：

4.根据权利要求1所述的单幅图像复原的解析迭代快速频谱外推方法，其特征在于解析迭代过程的步骤3中利用图像频谱合成校正，是将第k次的预测频谱F(w^k)和降质图像频谱F(u⁰)结合光学系统调制传递函数F(h)及其复共轭F(h)^*通过矩阵元素逐点相乘、相加和矩阵元素逐点相除的方式进行合成校正，得到合成校正频谱F(u^(k))，计算关系为：

F(u^(k))＝{F(h)^*οF(u₀)+γ^(k-1)F(W^(k))}/{F(h)^*οF(h)+γ^(k-1)}。

5.根据权利要求1所述的单幅图像复原的解析迭代快速频谱外推方法，其特征在于迭代系统包含梯度信息阈值α^(k-1)、预测贡献度参数β^(k-1)、合成贡献度参数γ^(k-1)三个关键参数。

6.根据权利要求5所述的单幅图像复原的解析迭代快速频谱外推方法，其特征是梯度信息阈值参数α^k-1作如下调整：梯度信息阈值α^(k-1)，在迭代过程中逐渐减少，且α⁰取特定的参数值10≤α⁽⁰⁾≤125，按照特定的递减顺序α^(k)/α^(k-1)＜1进行控制即满足α⁽⁰⁾＞α⁽¹⁾＞...α^(k-1)＞...＞0即可。

7.根据权利要求5所述的单幅图像复原的解析迭代快速频谱外推方法，其特征是预测贡献度参数β^k作如下调整：预测贡献度参数β^(k)需与α^(k-1)满足β^(k)＝c·α^(k)的正比关系，其中k∈(0.2，0.8)为实验参数。

8.根据权利要求5所述的单幅图像复原的解析迭代快速频谱外推方法，其特征是合成贡献度参数γ^k作如下调整：在迭代过程中按照递增顺序进行控制γ^(k)/γ^(k-1)＞2且γ⁽⁰⁾∈[1，20]。

Images