迭代去卷積算法

張麗娟，楊進華，蘇偉，姜成昊，王曉坤，譚芳

(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 系統設備參數引入估計模型?；诓ㄇ跋辔恍畔⒌腜SF 重構算法是將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)