张丽娟，杨进华，苏伟，姜成昊，王晓坤，谭芳

(1.长春理工大学 光电工程学院，吉林 长春130022;2.长春工业大学 计算机科学与工程学院，吉林 长春130012;3.长春理工大学 信息化中心，吉林 长春130022)

0 引言

在地对空观测成像时，地基光学望远镜由于大气湍流的动态干扰、观测系统误差及成像路径的光子噪声等因素的影响，尤其是极弱光条件下，例如星体目标观测，自适应光学(AO)系统只能实现波前误差的部分校正，目标的高频信息仍受到严重的抑制和衰减［1-2］。所以，经AO 系统校正后的图像还必须进行图像复原的后处理工作，才能获得较清晰的图像。目前，有很多克服大气湍流影响的图像复原算法，如迭代盲去卷积(IBD)算法，基于Richardson-Lucy 迭代的IBD(RL-IBD)算法，基于Wiener 的IBD(Wiener-IBD)算法，基于期望值最大化(EM)算法等［2-5］。

本文针对大气湍流造成的AO 图像模糊，提出了一种与正则化相结合的改进期望值最大化(RTIEM)算法的多帧AO 图像复原。具体内容为，首先研究基于波前相位信息的点扩散函数(PSF)重建方法，然后研究基于图像功率谱密度(PSD)及约束图像支持域的去噪方法，接着建立基于RT-IEM 算法的AO 图像代价函数模型及参数估计方法，最后给出基于RT-IEM 算法的多帧AO 图像高清复原算法的实现过程。在AO 图像复原实验中，分别对模拟AO 图像、实际观测的双星图像和恒星图像进行去卷积，验证本文提出算法的有效性。

1 基于RT-IEM 算法的AO 图像复原

AO 系统对光学波前误差进行实时测量-控制-校正，图1给出了成像补偿的AO 系统工作原理框图。该系统包含3 个基本组成部分:波前探测器、波前校正器和波前控制器。其工作原理:首先用波前探测器实时探测光学系统中静态和动态波前畸变;然后由波前控制器实时处理后计算波前校正器的电压控制信号;最后在波前校正器上加载这些控制信号，使其产生与输入波前畸变大小相等、符合相反的共轭波面，校正后的光束接近于平面波，使成像质量接近衍射极限。

由于像差的存在和AO 系统孔径的限制，AO 系统观测到的目标图像严重降质，为了提高目标图像的分辨率，本文将AO 系统波前补偿或校正和图像处理技术相结合。

图1 成像补偿的AO 系统工作原理Fig.1 Schematic diagram of AO system for imaging compensation

1.1 基于波前相位信息的PSF 重构方法

本小节提出的PSF 重构方法属于模型估计法，该方法将图像的退化因素考虑在内，在应用时基于大气湍流的模型，并将AO 系统设备参数引入估计模型。基于波前相位信息的PSF 重构算法是将AO系统的望远镜入瞳函数引入PSF 估计，以傅里叶光学作为数学基础。

理论上，AO 系统闭环校正后的PSF 数学模型为

由于受光学系统本身存在的焦距误差或光学偏离的影响，修改入瞳函数为

式中:θ(·)是波前相位误差，是由光学偏离或焦距误差引起的。

由(3)式可知，点扩散函数依赖于θ(·). 受大气湍流的影响，相位误差随时间变化，所以随时间变化的PSF 一般形式为

式中:θt(·)是随时间变化的相位误差，由大气湍流引发的。

就AO 系统而言，观测图像的曝光时间与湍流波动时间相比较短。假定曝光时间内目标强度不变，则序列图像模型为

式中:x 表示象素点灰度值;函数f 表示理想的目标原图像;g 表示观测到的退化图像;θtm为第m 帧短曝光AO 图像的湍流波前相位误差。

1.2 基于图像PSD 及约束图像支持域的去噪方法

为了减少AO 图像的噪声干扰，国内学者提出基于总变分的方法抑制噪声［6-7］，文献［8］分析了AO 图像在对称支持域内减少噪声的方法，该方法的结果分析基于真实图像方差的比率。与上述去噪算法不同，本节提出的是基于图像PSD 及约束图像支持域的去噪方法，并用图像的噪声变化函数评价去噪的效果。

地基设备观测的AO 图像去噪过程实际上是一种傅里叶频谱图像的噪声校正过程，即在图像的对称支持域内，以CCD 读入噪声和AO 噪声为主，因此要研究噪声的PSD. 下面对观测图像g(x)进行支持域约束来抑制噪声，令gc(x)为支持域约束图像，

式中:w(x)为在图像支持域内的权重函数。定义图像的功率谱为

式中:G(u)= S(u)× N(u)，S(u)为信号函数，N(u)为噪声函数;带宽函数W(u)=G(u)×δ，在有限支持域内W(u)很窄，接近于G(u)宽度×δ，δ 为0 ～1 之间的常数。

根据J.M.Conan 提出的定义，一副图像的PSD是它协方差的傅里叶变换，定义噪声的功率谱PSD(u)nn和未失真图像的功率谱PSD(u)ff［9］为

本文以圆形作为目标图像支持域，该支持域的直径取值范围是8 ～185 像素. 定义图像的噪声变化函数

1.3 改进的EM 算法

文献［10］基于极大似然准则推导出了多帧湍流退化图像的PSF 和目标图像的迭代求解公式，该方法复原图像的解空间较大，本文结合正则化技术对极大似然估计算法进行必要的修改，引入符合AO 图像物理事实意义上的约束，采用RT-IEM 算法对多帧AO 图像迭代求解较为理想的目标图像估计和PSF 的估计{}.

1.3.1 RT-IEM 算法

EM 算法是由Dempster 等［11］提出的求参数极大似然估计的迭代算法，该算法交替地进行两个基本计算步骤:E 步(计算期望)和M 步(极大化计算)。本文将正则化方法与先验知识相结合并应用到EM 算法中，从而实现多帧AO 图像的高清复原。

基于目标图像的先验知识和PSF 的联合概率分布函数，定义极大化函数将得到目标图像估计和PSF 估计，即

式中:z 为归一化因子;J(f，h)为代价函数。因此该问题转换为最小化代价函数J(f，h)，即

下面建立AO 图像的代价函数及其参数估计的优化模型。

(13)式中，第1 项表示噪声代价函数，第2 项是PSF 的约束代价函数，第3 项是目标边缘保持约束代价函数;a 为先验信息;α 是一常系数;σ2n是混合噪声的方差;θt是波前相位误差;i、j 表示象素点坐标;γ(·)是正则化函数，与各点梯度有关的可变正则化系数。PSF 的功率谱密度PSDh(u)=，对(x)进行傅里叶变换得到是PSF 均值的傅里叶变换，PSDh(u)的估计可以通过AO 系统闭环控制数据估计。采用Mugnier［12］同向性目标边缘约束理论，得到

式中:γ(x)是正则化调节参数。根据文献［13］定义γ(Δx)如下:

式中:参数ξ 为可供选择的平滑控制系数，ξ 取值大时平滑作用大，取值小时平滑作用小;Δx 为点x(i，j)(其中i，j 表示象素点坐标，x 表示象素点灰度值)4 个方向的梯度幅值。显然，有0 ＜γ(Δx)≤1.

本文将所有的先验信息都同时应用到目标估计和PSF 估计中，为了得到目标和PSF 的最优估计值，首先构造多帧矢量化的一维成像方程:

结合(13)式和(18)式建立多帧AO 图像联合去卷积的代价函数:

式中:m 为帧数。

1.3.2 复原算法实现

本文结合EM 算法的基本思想［9］，通过有限次循环运算，找到最佳PSF 的估计，达到极大程度地恢复目标图像的目的。

设观测图像的第m 帧观测结果gm(y)的完全数据集为{(y|x)}，gm(y)服从Poisson 分布的，它的期望为

根据数理统计理论，完全数据与不完全数据对应在统计意义上是对应的［13］，因此，与完全数据对应的代价函数的期望如下:

在迭代过程处理中，假设已知第n-1 次迭代数据，用En(Jmulti(f，h|g，a))表示当前代价函数的期望值。为了极小化代价函数的期望值，可将(21)式分别对hnm(x)和fn(x)求导数，并令导数为0，即

具体推导过程省略，最后解方程并整理得

本文复原算法RT-IEM 具体实现的步骤如下:

步骤1:初始化操作。

1)根据1.2 节的方法，对观测图像进行支持域约束、去噪等，并用一次维纳滤波的结果作为目标估计的初值

3)设定外循环最大次数Max_ iteration.

4)设定目标估计和PSF 估计的内循环最大次数Max_ count.

5)设定PSF 估计的内循环计数变量h_ count，目标的内循环计数变量f_ count.

步骤2:参数值的计算。

1)根据1.2 节的公式，估计每帧降质图像的混合噪声方差σ2n.

2)根据1.2 节的方法和公式，计算PSF 的功率谱密度PSDh(u)和估计波前相位误差θt.

3)根据1.3.1 节的方案和公式，估计目标边缘正则化参数α 和γ(Δx).

步骤3:迭代计算:i =0，1，…，Max_ iteration.

1)PSF 的内循环计数变量归0，即h_ count=0.

2)PSF 估计循环过程，j=0，1，…，Max_count .

b)h_ count + +;j+ +.

c)判断循环变量j 的值，若j ＜Max_count，则继续;否则，转向本步骤3.

3)目标估计的内循环计数变量归0，即:f_count =0.

4)目标估计的循环过程，k = 0，1，…，Max_count.

b)f_ count + +;k+ +.

c)判断循环变量k 值，若k ＜Max_count，则继续;否则，转向本步骤5.

5)判断外层循环是否结束，若i ＞max_iteration，则转至步骤4.

6)i=i+1，返回本步骤1 继续循环。

2 AO 图像复原实验结果及分析

2.1 仿真图像复原实验

利用本文提出的算法，进行了降质斑点图像模拟实验，图像复原的质量采用均方根误差RMSE 和改进的ΔSNR 进行定量评价。

式中:N1和N2分别表示图像在x 轴和y 轴的像素总数;(x，y)表示二维空域坐标;f(x，y)表示理想图像在点(x，y)处的灰度值。

评价估计的精度，定义为

式中:h(x，y)是由(1)式计算得到的PSF 真正值;(x，y)是由RT-IEM 算法计算得到的PSF 估计值;Sh表示PSF 的支持域，其取值范围是以5 ～16 像素为直径的圆形区域。

在本文的仿真实验中，原图像为多目标远距离的斑点状的红外模拟成像，采用中国科学院光学重点实验室的图像退化软件模拟生成10 帧序列退化图像，然后添加高斯白噪声，使图像的信噪比SNR为20 dB.设定参数与云南天文台1.2 m 自适应光学望远镜相同，望远镜成像系统的主要参数:大气相干长度r0=13 cm，望远镜直径D =1.08 m，光学系统的综合焦距l =22.42 m，中心波长λ =700 nm，CCD像素 大 小7.4 μm. 实 验 将 对RL-IBD 算 法［14］、Wiener-IBD算法［15］以及本文提出的RT-IEM 算法进行比较。

图2是原图像、模拟生成多帧退化图像及3 种算法的复原结果对比，图像大小为217×218 像素，其中图2(a)为原目标图像。采用图像退化仿真软件生成10 帧序列湍流退化图像，如图2(b)、图2(c)、图2(d)(为了节省篇幅，只显示3 帧，其余省略). 图2(e)是Wiener-IBD 复原图像，迭代次数为800 次，RMSE=16.35. 图2(f)是RL-IBD 复原图像，迭代次数为800 次，RMSE=10.31. 图2(g)是本文提出的RT-IEM 算法复原图像，波前相位误差θt(·)的方差为6.14×10-3，所以=8.26×10-3，σ2n=1.1×10-5，本次实验中正则化参数α =1.25，正则化函数的参数ξ =1.34，外循环次数为100 次，PSF 内循环次数为4，图像内循环次数为2，即总循环次数为600 次，RMSE = 10.05. 比较图2(e)、图2(f)、图2(g)的ΔSNR 和精度εe，见表1所示。比较可知，本文提出的RT-IEM 算法复原效果优于Wiener-IBD 算法和RL-IBD 算法，而且需要的迭代次数更少。

表1 3 种算法的ΔSNR 和εe 比较Tab.1 ΔSNRs and εe values calculated by three algorithms

图2 实验模拟降质图像及3 种算法的复原效果对比Fig.2 Degraded images and the comparison of the restored images of three algorithms

2.2 双星图像复原实验

在没有理想图像参照的情况下，采用无参照图像质量客观评价，本文采用半高全宽(FWHM)［16］作为图像复原效果的客观评价指标，FWHM 是指成像峰值与半峰值之间像素距离的2 倍，包括x 方向和y 方向的FWHM，是天文观测领域较好的评价指标。

利用本文算法对无参照图像的双星图像进行了复原实验。实验中的双星图像由中国科学院云南天文台1.2 m AO 望远镜于2006年12月3日采集的。成像CCD 大小为320×240 像素阵列，成像系统的主要参数如下:大气相干长度r0=13 cm，望远镜全孔径D=1.03 m，成像观测波段为0.7 ～0.9 μm，中心波长λ=0.72 μm，CCD 像素大小为6.7 μm，光学成像焦距l =20 m. 复原实验将对RL-IBD 算法、Wiener-IBD 算法以及本文提出的RT-IEM 算法进行比较，本文采用帧选择技术［4］，从100 帧降质图像中遴选出40 帧作为盲解卷积的待处理图像。

图3是多帧退化图像及3 种算法的复原结果对比，其中图3(a)、图3(b)、图3(c)和图3(d)为观测的双星图像(为了节省篇幅，只显示4 帧，其余省略)，图3(e)是本文算法估计的PSF，图3(f)是本文算法复原图像，波前相位误差θt(·)的方差为6.04× 10-3，所以= 9.03× 10-3，σ2n= 1.32×10-5.本次实验中参数选择为:正则化参数α =1.15，正则化函数的参数ξ =1.34，外循环次数为50 次，PSF 内循环次数为4 次，图像内循环次数为2 次，即总循环次数为300 次，FWHM =6.17 像素，复原结果已经非常接近1.2 m AO 望远镜的衍射极限。图3(g)是RL-IBD 算法的复原效果，FWHM =5.62 像素. 图3(h)是Wiener-IBD 算法的复原效果，FWHM=4.76 像素. 图4为观测图像和3 种复原算法复原图像的能量示意图。比较可知，本文提出的算法复原效果优于Wiener-IBD 算法和RL-IBD算法，而且需要的迭代次数更少。

2.3 恒星观测实验

采用本文算法还进行了恒星图像复原实验，实验中的恒星图像由中国科学院云南天文台1.2 m AO 望远镜于2006年12月1日采集的。与双星图像的光学成像系统参数相同，本文算法实验过程中的参数选择类似。图5为多帧恒星观测图像及3 种算法的复原结果对比。图5(a)、图5(b)、图5(c)和图5(d)是观测的某恒星图像(为了节省篇幅只显示4 帧)。图5(e)为 Wiener-IBD 复原图像。图5(f)为RL-IBD 复原图像。图5(g)为本文算法估计的PSF，迭代初值= 9.56× 10-3，PSF 迭代6 次。图5(h)为本文算法复原图像。图5(i)为恒星观测图像的能量图。图5(j)为本文算法复原图像的能量图。图5(k)表示3 种算法复原图像对比度和分辨率的对比曲线。

图3 多帧观测图像及3 种算法的复原结果对比Fig.3 Multi-frame original images and restoration comparison of three algorithms

图4 观测图像和3 种复原算法复原图像的能量示意图Fig.4 Energy spectra of observed images and the comparison of restored images of three algorithms

由实验结果可以看出，本文算法去卷积后降质AO 图像的恢复质量得到明显提升，能量更为集中且分辨率更高。当然，在有限帧数及迭代次数情况下，是不可能将目标图像完全地恢复出来的，当AO图像帧数增加时，复原效果将越来越好，但计算的时间和计算量也随之增加。

图5 多帧恒星观测图像及3 种算法的复原结果对比Fig.5 The multi-frames observed stellar images and restoration comparison of three algorithms

3 结论

针对AO 图像的特点，在EM 算法的基础上，提出了改进的期望值最大化RT-IEM 算法用于AO 图像高清晰复原。仿真实验结果证明，RT-IEM 算法能准确辨识出AO 图像的PSF，PSF 的估算精度优于Wiener-IBD、RL-IBD 算法，且本文算法的迭代次数减少14.3%. 因此，本文的研究结果对实际AO 图像复原有一定的应用价值。

References)

［1］ 李大禹，胡立发，穆全全，等.高准确度LCOS 自适应光学成像系统的研究［J］. 光子学报，2008，37(3):506 -508.LI Da-yu，HU Li-fa，MU Quan-quan，et al. A high-resolution liquid crystal adaptive optics system［J］. Acta Photonica Sinica，2008，37(3):506 -508.(in Chinese)

［2］ 饶长辉，沈锋，姜文汉.自适应光学系统波前校正残余误差的功率谱分析方法［J］.光学学报，2000，20(1):68 -73.RAO Chang-hui，SHEN Feng，JIANG Wen-han. Analysis of closed-loop wavefront residual error of adaptive optical system using the method of power spectrum［J］. Acta Optica Sinica，2000，20(1):68 -73.(in Chinese)

［3］ 陈波.自适应光学图像复原理论与算法研究［D］. 郑州:解放军信息工程大学，2008:44 -46 CHEN Bo. The theory and algorithms of adaptive optics image restoration［D］. Zhengzhou:PLA Information Engineering University，2008:44 -46.(in Chinese)

［4］ Tian Y，Rao C H，Wei K. Postprocessing of adaptive optics images based on frame selection and multiframe blind deconvolution［C］∥Adaptive Optics Systems. Marseille，France:SPIE，2008.

［5］ Schulz T J. Nonlinear models and methods for space-object imaging through atmospheric turbulence［C］// The 1996 IEEE International Conference on Image Processing. Lausanne，Switzerland:IEEE，1996.

［6］ 路雅宁，郭雷，李晖晖. 带限剪切波变换与全变差结合的图像去噪［J］.光子学报，2013，42(12):1430 -1435.LU Ya-ning，GUO Lei，LI Hui-hui. Total variation based bandlimited sheralets transform for image denoising［J］. Acta Photonica Sinica，2013，42(12):1430 -1435.(in Chinese)

［7］ 赵杰，杨建雷.基于Cycle Spinning Contourlet 变换和总变分最小化的遥感图像去噪算法［J］.光子学报，2010，39(9):1659 -1660.ZHAO Jie，YANG Jian-lei. Remote sensing image denoising algorithm based on fusion theory using cycle spinning contourlet transform and total variation minimization［J］. Acta Photonica Sinica，2010，39(9):1659 -1660.(in Chinese)

［8］ Matson C L，Roggemann M C . Noise reduction in adaptive optics imagery with the use of support constraints［J］. Applied Optics，1995，34(5):767 -780.

［9］ Tyler D W，Matson C L. Reduction of nonstationary noise in telescope imagery using a support constraint ［J］. OSA Optics Express，1997，1(11):347 -350.

［10］ Zhang L J，Yang J H，Su W. Research on blind deconvolution algorithm of multiframe turbulence-degraded images［J］. Journal of Information and Computational Science，2013，10 (12):3625 -3633.

［11］ Dempster A P，Laird N M，Rubin D B. Maximum likelihood from incomplete data via the EM algorithm［J］. Journal of the Royal Statistical Society Series B，1977，39(1):2 -20.

［12］ Mugnier L M，Conan J M，Fusco T，et al. Joint maximum S posteriori estimation of object and PSF for turbulence degraded images［C］// Bayesian Inference for Inverse Problems. San Diego，California:SPIE，1998.

［13］ 洪汉玉.成像探测系统图像复原算法研究［D］. 武汉:华中科技大学，2004:51 -107.HONG Han-yu. Research on image restoration algorithms in imaging detection system［D］. Wuhan:Huazhong University of Science and Technology，2004:51 -107.(in Chinese)

［14］ Tsumuraya F. Iterative blind deconvolution method using Lucy’s algorithm［J］. Astronamy and Astrophysics，1994，282(2):699-708.

［15］ 张士杰，李俊山，杨亚威，等. 湍流退化红外图像降晰函数辨识［J］. 光学精密工程，2013，21(2):514 -520.ZHANG Shi-jie，LI Jun-shan，YANG Ya-wei，et al. Blur identification of turbulence-degraded IR images［J］. Optics and Precision Engineering，2013，21(2):514 -520.(in Chinese)

［16］ 陈波，程承旗，郭仕德，等. 自适应光学图像非对称图像迭代盲复原算法［J］. 强激光与粒子束，2011，23(2):313 -318.CHEN Bo，CHENG Cheng-qi，GUO Shi-de，et al. Unsymmetrical multi-limit iterative blind deconvolution algorithm for adaptive optics image restoration［J］. High Power Laser and Particle Beams，2011，23(2):313 -318.(in Chinese)