稀疏先驗型的大氣湍流退化圖像盲復(fù)原

周海蓉，田 雨，饒長輝*

稀疏先驗型的大氣湍流退化圖像盲復(fù)原

周海蓉1,2,3，田 雨1,2，饒長輝1,2*

1中國科學(xué)院自適應(yīng)光學(xué)重點(diǎn)實驗室，四川 成都 610209；2中國科學(xué)院光電技術(shù)研究所，四川 成都 610209；3中國科學(xué)院大學(xué)，北京 100049

圖像盲復(fù)原是僅從降質(zhì)圖像就恢復(fù)出模糊核和真實銳利圖像的方法，由于其病態(tài)性，通常需要加入圖像先驗知識約束解的范圍。針對傳統(tǒng)的圖像梯度2和1范數(shù)先驗不能真實刻畫自然圖像梯度分布的特點(diǎn)，本文將圖像梯度稀疏先驗應(yīng)用于單幀大氣湍流退化圖像盲復(fù)原中。先估計模糊核再進(jìn)行非盲復(fù)原，利用分裂Bregman算法求解相應(yīng)的非凸代價函數(shù)。仿真實驗表明，與總變分先驗(1范數(shù))相比，稀疏先驗有利于模糊核的估計、產(chǎn)生銳利邊緣和去除振鈴等，降低了模糊核的估計誤差從而提高了復(fù)原質(zhì)量。最后對真實湍流退化圖像進(jìn)行了復(fù)原。

自適應(yīng)光學(xué)；稀疏先驗；盲解卷積；分裂Bregman

1 引 言

地基望遠(yuǎn)鏡對太空目標(biāo)(如衛(wèi)星、恒星或太陽)觀測時，受地球大氣湍流影響，其成像分辨力遠(yuǎn)遠(yuǎn)低于衍射極限，通過自適應(yīng)光學(xué)(adaptive optics, AO)系統(tǒng)能實時補(bǔ)償湍流造成的波前畸變[1]，但由于系統(tǒng)自身存在誤差如變形鏡擬合誤差、時間帶寬誤差和非等暈誤差等，閉環(huán)校正后的圖像仍然存在殘差，需要進(jìn)行事后處理以進(jìn)一步提升圖像質(zhì)量。常用的事后處理方法有相位差法、斑點(diǎn)圖像重建法和盲解卷積法。相位差法需要一套額外的成像裝置，而斑點(diǎn)圖像重建法則是利用上百幀短曝光的圖像恢復(fù)出目標(biāo)的相位和振幅信息[2]。與前兩種方法相比，盲解卷積法僅需一幀或幾幀降質(zhì)圖像即可實現(xiàn)重建，方法更顯簡單靈活。

盲復(fù)原算法在求解策略上可分為聯(lián)合估計法和兩階段估計法，前者同時恢復(fù)出目標(biāo)和模糊函數(shù)，后者先估計出模糊函數(shù)再進(jìn)行非盲復(fù)原。Levin[3]表明傳統(tǒng)的最大后驗概率聯(lián)合估計法MAP,k(maximum a posteriori, MAP)傾向于得到無模糊解導(dǎo)致復(fù)原失敗，而MAP先估計出模糊核的兩階段方法會成功。無模糊解即復(fù)原圖像為退化圖像，復(fù)原模糊核為脈沖響應(yīng)。但MAP,k方法在實際實驗中取得成功[4]。為解決這個矛盾，Perrone[4]首先證明了Levin的正確觀點(diǎn)，其次表明MAP,k成功原因在于歸一化步驟是順序進(jìn)行的，因此MAP,k和MAP兩種求解方法都是可行的，盲解卷積問題的關(guān)鍵還是在于先驗的選擇。

盲解卷積方法的困難之處在于方程存在無窮解，如真解和無模糊解(no-blur solution)都是問題的解?？朔そ饩矸e病態(tài)性的方法是通過加入先驗來約束解的范圍，從而避免不需要的解。傳統(tǒng)的迭代盲解卷積算法(iterative blind deconvolution, IBD)使用正性約束和能量守恒約束進(jìn)行交替求解，然而算法沒有穩(wěn)定收斂性保證[5]。Lane[6]將先驗化為代價函數(shù)的一個正則化項，使用共軛梯度法進(jìn)行聯(lián)合求解，算法雖然克服了病態(tài)性問題但容易陷入局域解。Jefferies等[7]在Lane的基礎(chǔ)上進(jìn)一步優(yōu)化，加入帶寬有限約束和傅里葉模量約束。其后You和Kaveh[8]提出將圖像梯度的1范數(shù)作為圖像平滑約束，但由于1范數(shù)是各向同性的，在復(fù)原同時無法保留圖像細(xì)節(jié)，Chan和Wong[9]在此基礎(chǔ)上引入圖像總變分(total variation, TV)約束，并設(shè)計了交替最小化的求解方案。與傳統(tǒng)圖像盲復(fù)原算法相比，當(dāng)前算法都引入自然圖像的統(tǒng)計信息作為先驗，并廣泛采用MAP兩階段估計法進(jìn)行求解。Fergus[10]使用混合高斯模型對圖像梯度進(jìn)行建模，Krishnan[11]采用另一種圖像先驗1/2使得銳利圖像的代價函數(shù)最小，并采用迭代再加權(quán)最小二乘(iterative reweighted least squares, IRLS)算法求解模糊核。Kotera等[12]采用圖像梯度的稀疏范數(shù)l(0<<1)去除圖像中的運(yùn)動模糊。

真實的自然圖像梯度[3]表現(xiàn)為稀疏的長拖尾分布。與傳統(tǒng)的總變分先驗[9]相比，稀疏范數(shù)l更能真實地刻畫這種稀疏分布，常用于圖像去運(yùn)動模糊的盲復(fù)原領(lǐng)域[12]，本文將這種稀疏先驗引入大氣湍流降質(zhì)圖像的盲復(fù)原中。首先，稀疏先驗將產(chǎn)生非凸的代價函數(shù)，本文采用分裂Bregman優(yōu)化算法[13]和查表法[14]進(jìn)行數(shù)值求解。其次，采用Krishnan[14]的稀疏先驗方法復(fù)原出銳利圖像，避免了Fergus等[10]出現(xiàn)的振鈴現(xiàn)象。最后，分別對仿真的湍流退化圖像和真實觀測圖像進(jìn)行圖像復(fù)原，驗證了算法的有效性。

2 基于稀疏的盲解卷積復(fù)原算法

當(dāng)前盲復(fù)原方法路線主要分為基于概率統(tǒng)計的貝葉斯法、正則化法以及變分貝葉斯估計，貝葉斯估計是假設(shè)圖像、模糊核或噪聲服從某種分布作為先驗知識加入最大似然函數(shù)，先驗即轉(zhuǎn)化為代價函數(shù)中的正則化項，因此貝葉斯估計與正則化方法實際是等價的。用()表示真實圖像，()為模糊核，()為噪聲，()表示觀測圖像，“*”表示空域卷積，圖像的線性退化過程可表示為

2.1 圖像的統(tǒng)計模型

表1羅列了常用的圖像先驗?zāi)Ｐ?，其中?”表示梯度算子。傳統(tǒng)的盲復(fù)原算法采取的先驗是基于圖像的強(qiáng)度信息，如維納濾波迭代盲解卷積假設(shè)圖像強(qiáng)度服從高斯分布。其后學(xué)者采用圖像的梯度作為先驗，如You和Kaveh[6]使用梯度1范數(shù)(即Gaussian分布)，Chan[7]采用梯度的1范數(shù)(即Laplace分布)有利于保持邊緣和細(xì)節(jié)。對自然圖像的梯度進(jìn)行統(tǒng)計時，其對數(shù)分布直方圖如圖1(a)所示，其統(tǒng)計分布為非高斯的長拖尾分布。若采用參數(shù)模型進(jìn)行擬合，如圖1(b)表明Gaussian先驗和Laplace先驗均不能近似這種分布，而稀疏先驗(<1)更能近似真實圖像梯度的長拖尾分布，這種先驗有利于產(chǎn)生銳利的邊緣、減少噪聲并去除振鈴現(xiàn)象[13]。

表1 常用的圖像先驗?zāi)Ｐ?/p>

2.2 算法框架

將稀疏先驗和模糊核先驗加入代價函數(shù)后的最小化目標(biāo)函數(shù)為

直接求解式(3)是困難的，Jefferies等[7]使用共軛梯度聯(lián)合求解和，但Levin[3]表明這種聯(lián)合求解策略容易陷入無模糊，采用先估計模糊核再非盲復(fù)原的兩階段方法會成功。本文采用交替最小化方案[9]可將式(3)轉(zhuǎn)換為兩個子問題交替求解。Chan[15]證明了交替最小化方案的收斂性，其中模糊核子問題和圖像核子問題分別如下：

與總變分盲解卷積[9]類似，稀疏先驗導(dǎo)致代價函數(shù)式(4)非凸，無法得到閉式解，需要進(jìn)行迭代求得數(shù)值近似解。Levin[16]提出IRLS算法求解這種l范數(shù)最優(yōu)化問題，但使用共軛梯度法需要迭代上百次。Goldstein[13]提出分裂Bregman算法求解總變分盲解卷積，Krishnan[14]提出查表法(look-up table, LUT)將分裂Bregman算法推廣到求解含有l范數(shù)的非凸解卷積問題，算法速度比IRLS提高了幾個數(shù)量級。本文采用分裂Bregman算法和查表法來求解上述兩個子問題式(4)和式(5)。

2.2.1 模糊核子問題

首先，由分裂Bregman算法思想，通過引入輔助變量v=h將式(4)化為一個有約束最小化問題，再引入Bregman變量b形成二次懲罰項進(jìn)行迭代求解，為二次懲罰項的權(quán)重，最終的最小化目標(biāo)函數(shù)為

加入輔助變量后，可以推導(dǎo)出和v的頻域閉式解，同時更新Bregman變量b，總的迭代求解過程如下：

圖1 自然圖像統(tǒng)計模型。(a) 真實自然圖像的梯度分布；(b) 參數(shù)模型

Fig. 1 Natural images statistical model. (a) Real natural image gradient distribution; (b) Parametric model

和歸一化：

約束。

2.2.2 圖像子問題

進(jìn)行變量分離后，在頻域推導(dǎo)出的閉式解，同時更新Bregman變量，求解過程如下：

其中：

2.2.3子問題

模糊核子問題中的變量v是存在閉式解的，可以用軟閾值求解。而子問題是一個l范數(shù)非凸問題，其最小化方程為

等價于求解最小化問題：

當(dāng)=1時，化為一個1正則化問題[13]，其解為

也寫為

當(dāng)0<<1時，采用Krishnan[11]查表法(lookup table, LUT)可得到近似解。

2.2.4 避免陷入局域解的策略

稀疏先驗將導(dǎo)致卡通效應(yīng)，在算法迭代時以1.5倍增大保真項的權(quán)重，即逐漸減小正則化的強(qiáng)度，有利于降低復(fù)原圖像的分段平滑效應(yīng)。此外，在每個子問題迭代時，迭代次數(shù)不宜過多，防止算法陷入局域解中[4]，設(shè)置算法的迭代停止條件：

3 實驗結(jié)果與分析

3.1 大氣湍流退化圖像仿真實驗

大氣湍流是限制大型天文望遠(yuǎn)鏡成像分辨率的主要原因，它造成的隨機(jī)波前畸變可以用Zernike多項式表示。它的每一項都表示傳統(tǒng)光學(xué)系統(tǒng)的基本像差，如離焦、像散和慧差等，這種方法得到的數(shù)值模擬效果比較符合Kolmogrov譜的湍流模型[17]。波前相位中的Zernike系數(shù)近似服從零均值的高斯分布(圖2)。測試圖像為OCNR5衛(wèi)星(圖3(a))，大小為256 pixels×256 pixels，降質(zhì)圖像為圖3(b)，利用Zernike前65項多項式生成點(diǎn)擴(kuò)散函數(shù)(point spread function, PSF)如圖3(b)右下角所示。定義核估計的相對誤差為

其中：hr為真實核，he為估計核。數(shù)值越小表示估計的核與真實核越接近。假設(shè)真實圖像為f(x,y)，降質(zhì)圖像為g(x,y)，使用圖像信噪比(signal to noise ratio, SNR)衡量復(fù)原圖像的質(zhì)量，數(shù)值越大說明復(fù)原質(zhì)量越好，定義為

使用稀疏先驗盲復(fù)原算法對降質(zhì)圖像進(jìn)行復(fù)原，算法中的參數(shù)初始化為：，，，，。將復(fù)原結(jié)果與維納濾波迭代盲復(fù)原和總變分盲復(fù)原進(jìn)行對比，圖3(c)~3(h)為三種算法的復(fù)原結(jié)果，圖中右下角為對應(yīng)的估計模糊核。從圖3(c)可以看出，維納濾波盲復(fù)原的PSF估計誤差大，信噪比低，但仍恢復(fù)出目標(biāo)的大尺度輪廓?？傋兎置?fù)原(圖3(d))能保持更多的細(xì)節(jié)，估計的PSF結(jié)構(gòu)和形狀更接近真實的PSF，但仍存在殘余模糊。圖3(e)~(h)分別為使用p值為0.9、0.6、0.3和0.1的稀疏先驗進(jìn)行盲復(fù)原的結(jié)果，可以看出p的取值過大(p=0.9)時，出現(xiàn)了顯著的振鈴現(xiàn)象；而p的取值過小(p=0.1)時，核估計誤差增大，細(xì)節(jié)和輪廓被噪聲淹沒；p的取值為0.6時，復(fù)原圖像具有了豐富的細(xì)節(jié)和輪廓；當(dāng)p的取值為0.3時，核估計誤差最小，對應(yīng)的復(fù)原圖像信噪比值最高，復(fù)原質(zhì)量最好。表2給出了不同算法的核估計誤差和信噪比值。圖4給出了三種算法的收斂情況，維納濾波盲復(fù)原算法的核估計相對誤差先是快速下降，但收斂緩慢，加入分裂Bregman算法優(yōu)化后，總變分盲復(fù)原和稀疏先驗盲復(fù)原算法都能穩(wěn)定地收斂。

表2 不同算法的核估計誤差和信噪比

表2還給出了不同算法的運(yùn)行時間，維納濾波盲復(fù)原算法直接在頻域求得解析解，單個子問題不需要迭代，因此時間消耗少僅需10.34 s；在使用分裂Bregman算法進(jìn)行數(shù)值計算后，總變分盲解卷積需要17.30 s；而稀疏先驗盲解卷積算法由于圖像子問題非凸需要查表，時間消耗為24.34 s(=0.9)、27.78 s(=0.6)、25.12 s(=0.3)和28.25 s(=0.1)。三種算法均是在Matlab2017Ra平臺上運(yùn)行，測試計算機(jī)系統(tǒng)為Window7旗艦版64位操作系統(tǒng)，處理器為AMD A10 PRO-7800B R7，12 Compute Cores 4C+8G @3.50 GHz，內(nèi)存為8 G。

圖4 算法收斂性比較

3.2 p值的選取原則

圖5 不同p值下的復(fù)原圖像SNRs和核估計誤差RE

3.3 真實湍流退化圖像復(fù)原

真實退化圖像(圖6(a))為口徑80 cm的Munin望遠(yuǎn)鏡采集到的國際空間站(International Space Station, ISS)圖像，采集時間為2007年8月1日(來自www.tracking-station.de/images/images.html)。受大氣湍流影響，圖像明顯降質(zhì)，同時還含有大量噪聲。分別使用總變分盲復(fù)原、稀疏先驗=0.5和稀疏先驗=0.1進(jìn)行復(fù)原重建?？梢钥闯銮罢?圖6(b))的圖像質(zhì)量并未顯著改善，而稀疏先驗=0.5時(圖6(c))，空間站的太陽能板的邊緣更加清晰，其細(xì)節(jié)部分得到增強(qiáng)。而當(dāng)?shù)娜≈颠^小(=0.1)時，復(fù)原圖像(圖6(d))伴隨著顯著的噪聲放大。圖6(e)中分別給出了退化圖像、總變分盲復(fù)原、稀疏先驗=0.5和稀疏先驗=0.1的復(fù)原圖像的細(xì)節(jié)對比圖。

4 結(jié) 論

圖像總變分先驗無法真實表征自然圖像的梯度分布，盲去運(yùn)動模糊常使用圖像梯度的稀疏先驗l(0<<1)來近似這種分布。本文將稀疏先驗引入到大氣湍流退化圖像盲復(fù)原中，運(yùn)用分裂Bregman算法和查表法求解稀疏導(dǎo)致的非凸代價函數(shù)，同時在非盲步驟中采用稀疏先驗復(fù)原出銳利圖像，避免了振鈴現(xiàn)象。實驗表明，與總變分先驗相比，稀疏先驗有利于模糊核的估計，降低了核估計誤差，提升了復(fù)原質(zhì)量。同時，值(0<<1)的選取不宜過大或過小，取值范圍在[0.3，0.7]之間復(fù)原質(zhì)量較高。最后，真實退化圖像的復(fù)原結(jié)果表明了算法的有效性和穩(wěn)定性。

[1] 姜文漢. 自適應(yīng)光學(xué)技術(shù)[J]. 自然雜志, 2006, 28(1): 7–13.

Jiang W H. Adaptive optical technology[J]., 2006, 28(1): 7–13.

[2] Bao H, Rao C H, Tian Y,. Research progress on adaptive optical image post reconstruction[J].,2018, 45(3): 58–67.

鮑華, 饒長輝, 田雨, 等. 自適應(yīng)光學(xué)圖像事后重建技術(shù)研究進(jìn)展[J]. 光電工程, 2018, 45(3): 58–67.

[3] Levin A, Weiss Y, Durand F,. Understanding and evaluating blind deconvolution algorithms[C]//, Miami, 2009: 1964–1971.

[4] Perrone D, Favaro P. Total variation blind deconvolution: the devil is in the details[C]//,Columbus,2014: 2909–2916.

[5] Ayers G R, Dainty J C. Iterative blind deconvolution method and its applications[J]., 1988, 13(7): 547–549.

[6] Lane R G. Blind deconvolution of speckle images[J]., 1992, 9(9): 1508–1514.

[7] Jefferies S M, Christou J C. Restoration of astronomical images by iterative blind deconvolution[J]., 1993, 415(2): 862–874.

[8] You Y L, Kaveh M. A regularization approach to joint blur identification and image restoration[J]., 1996, 5(3): 416–428.

[9] Chan T F, Wong C K. Total variation blind deconvolution[J]., 1998, 7(3): 370–375.

[10] Fergus R, Singh B, Hertzmann A,. Removing camera shake from a single photograph[J]., 2006, 25(3): 787–794.

[11] Krishnan D, Tay T, Fergus R. Blind deconvolution using a normalized sparsity measure[C]//, Providence, 2011: 233–240.

[12] Kotera J, ?roubek F, Milanfar P. Blind deconvolution using alternating maximum a posteriori estimation with heavy-tailed priors[C]//, York, 2013: 59–66.

[13] Goldstein T, Osher S. The split bregman method for L1-regularized problems[J]., 2009, 2(2): 323–343.

[14] Krishnan D, Fergus R. Fast image deconvolution using hyper-Laplacian priors[C]//, Vancouver, 2009: 1033–1041.

[15] Chan T F, Wong C K. Convergence of the alternating minimization algorithm for blind deconvolution[J]., 2000, 316(1–3): 259–285.

[16] Levin A, Fergus R, Durand F,. Image and depth from a conventional camera with a coded aperture[J].,2007, 26(3): 70.

[17] Wang Q T, Tong S F, Xu Y H. On simulation and verification of the atmospheric turbulent phase screen with Zernike polynomials[J]., 2013, 42(7): 1907–1911.王齊濤, 佟首峰, 徐友會. 采用Zernike多項式對大氣湍流相位屏的仿真和驗證[J]. 紅外與激光工程, 2013, 42(7): 1907–1911.

Blind restoration of atmospheric turbulence degraded images by sparse prior model

Zhou Hairong1,2,3, Tian Yu1,2, Rao Changhui1,2*

1Key Laboratory of Adaptive Optics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;2Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;3University of Chinese Academy of Sciences, Beijing 100049, China

(a) Simulation: degraded and restoration; (b) Real observation: degraded and restoration

Overview:Atmospheric turbulence is a major factor limiting the imaging resolution of ground-based telescopes. Adaptive optics (AO) is commonly used to compensate for the wavefront distortion caused by turbulence to obtain higher resolution. However, due to the limitations of the system itself, such as the fitting error of the deformed mirror and the residual error caused by the time bandwidth, the closed-loop image still has residual errors, thus AO postprocessing technique is needed to further improve the image quality.

Blind deconvolution (BD) could recover a sharp image only from several degraded images. However, BD problems have the difficulties of ill-conditioned and infinite solutions, it is necessary to add prior knowledge to avoid undesired solutions. Traditional Wiener filtering based iterative blind deconvolution method assumes that the intensity of the image obeys the Gaussian distribution, while Chan et al. employ total variation prior, which assumes the gradient of the image obeys the Laplacian distribution. Given the fact that the gradient distribution of natural images is a sparse one with heavy tails, both of Gaussian and Laplacian model cannot approximate this sparse model greatly. Therefore, this paper draws on the image gradient sparse priori derived from blind motion deblurring, and applied it to the blind restoration of turbulence-degraded images.

In order to cope with the non-convex cost function caused by sparse prior, this paper uses split Bregman and look-up table method to solve effectively. Secondly, a two-step estimation strategy is adopted including kernel estimation and non-blind restoration. This paper employs the deconvolution method of Krishnan to reconstruct a sharp image with the kernel from the former step, and this strategy is essential to avoid the ringing effect reported by Fergus.

Firstly, the simulation experiment is carried out. To model the atmospheric turbulence degradation, the Zernike polynomials are used to generate the point spread function, and the OCNR5 satellite image is used as an object for observing. This paper adopts the relative error to evaluate the kernel estimation error and the signal to noise ratio (SNR) to evaluate image quality of degraded and restored ones. Both simulations and experiments on the real degraded images show that: 1) compared with the traditional Wiener filtering method and total variation prior, the sparse prior is beneficial to kernel estimation, produces sharp edges and removes ringing, thus improving the restoration quality. 2) After employing split Bregman optimization, restoration with sparse prior can converge rapidly and steadily, hence the proposed algorithm in this paper is robust and stable. 3) It is worth noting that the value ofshould not be too large or too small, smallerwill amplify the noise.

Citation: Zhou H R, Tian Y, Rao C H. Blind restoration of atmospheric turbulence degraded images by sparse prior model[J]., 2020, 47(7): 190040

Blind restoration of atmospheric turbulence degraded images by sparse prior model

Zhou Hairong1,2,3, Tian Yu1,2, Rao Changhui1,2*

1Key Laboratory of Adaptive Optics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;2Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;3University of Chinese Academy of Sciences, Beijing 100049, China

Blind image deconvolution is one method of restoring both kernel and real sharp image only from degraded images, due to its illness, image priors are necessarily applied to constrain the solution. Given the fact that traditional image gradient2and1norm priors cannot describe the gradient distribution of natural images, in this paper, the image sparse prior is applied to the restoration of single-frame atmospheric turbulence degraded images. Kernel estimation is performed first, followed by non-blind restoration and the split Bregman algorithm is used to solve the non-convex cost function. Simulation results show that compared with total variation priori, sparse priori is better at kernel estimation, producing sharp edges and removal of ringing, etc., which reducing the kernel estimation error and improving restoration quality. Finally, the real turbulence-degraded images are restored.

adaptive optics; sparse prior; blind deconvolution; split Bregman

TP391

A

10.12086/oee.2020.190040

: Zhou H R, Tian Y, Rao C H. Blind restoration of atmospheric turbulence degraded images by sparse prior model[J]., 2020,47(7): 190040

周海蓉，田雨，饒長輝. 稀疏先驗型的大氣湍流退化圖像盲復(fù)原[J]. 光電工程，2020，47(7): 190040

Supported by National Natural Science Foundation of China (11727805, 11703029)

* E-mail: chrao@ioe.ac.cn

2019-01-24；

2019-03-15

國家自然科學(xué)基金資助項目(11727805，11703029)

周海蓉(1994-)，女，碩士，主要從事自適應(yīng)光學(xué)圖像盲復(fù)原的研究。E-mail：zhouhairongwhu@foxmail.com

饒長輝(1971-)，男，博士，研究員，主要從事大口徑高分辨力光學(xué)成像望遠(yuǎn)鏡技術(shù)和系統(tǒng)研制工作。E-mail：chrao@ioe.ac.cn