姜 超，耿則勛，劉立勇，潘映峰

（1.解放軍信息工程大學(xué) 地理空間信息學(xué)院，鄭州450052；2.中國(guó)科學(xué)院 國(guó)家天文臺(tái)，北京100012；3.中國(guó)人民解放軍61175部隊(duì)，南京210049）

0 引 言

光學(xué)探測(cè)裝置（如航空航天光學(xué)成像傳感器和地基光學(xué)望遠(yuǎn)鏡）透過(guò)大氣對(duì)感興趣目標(biāo)進(jìn)行觀測(cè)時(shí)通常會(huì)受到大氣湍流和噪聲的干擾。這種干擾會(huì)造成觀測(cè)目標(biāo)圖像質(zhì)量嚴(yán)重退化，影響后續(xù)的目標(biāo)提取和識(shí)別等操作。為了解決大氣湍流引起的圖像模糊問(wèn)題，早期的研究人員提出了逆濾波、卡爾曼濾波以及維納濾波等方法，但上述濾波算法要求退化圖像的點(diǎn)擴(kuò)散函數(shù)（Point spread function，PSF）極大似然估計(jì)算法（Maxium likelihood，ML），其中PSF是精確已知的。然而，建立一個(gè)精確的數(shù)學(xué)模型來(lái)描述高度隨機(jī)的大氣湍流運(yùn)動(dòng)是非常困難的。因此，在PSF未知或部分已知的情況下，通常采用盲解卷積技術(shù)對(duì)目標(biāo)圖像和點(diǎn)擴(kuò)散函數(shù)進(jìn)行估計(jì)。常用的盲解卷積算法有：迭代盲解卷積（Iterative blind deconvolution，IBD）、非負(fù)支持域遞歸逆濾波（Non－negativity and support constraints recursive inverse filtering，NASRIF）、模擬退火算法（Simulated annealing，SA）以及極大似然估計(jì)算法。其中，IBD 算法需要對(duì)PSF的支持域進(jìn)行較緊的約束且算法的收斂性不夠好，通常只能得到一個(gè)近似的解［1］。SA 算法雖然具有全局收斂性，但因該算法計(jì)算量過(guò)大，而嚴(yán)重限制了其適用范圍［2］。NAS－RIF 算法需要先確定目標(biāo)的支持域并假定理想圖像非負(fù)，其對(duì)背景均勻的退化圖像能取得較好的復(fù)原效果，但該算法對(duì)噪聲敏感［3］?；贛L的盲解卷積算法是建立在真實(shí)的成像模型基礎(chǔ)上的一種簡(jiǎn)單實(shí)用的最優(yōu)化方法；通過(guò)對(duì)解施加約束可以對(duì)該算法進(jìn)行擴(kuò)展，從而推導(dǎo)出其他算法。2008 年，Benvenuto等［4］在傳統(tǒng)ML算法的基礎(chǔ)上提出了基于混合噪聲模型的ML 算法并證明了解的存在性。但該算法沒(méi)有充分考慮迭代過(guò)程中的噪聲影響，同時(shí)假定PSF已知且在迭代過(guò)程中保持不變，從而導(dǎo)致復(fù)原結(jié)果對(duì)噪聲敏感且依賴于PSF的初始估計(jì)。

作者針對(duì)Benvenuto算法的不足，首先，對(duì)圖像進(jìn)行去噪預(yù)處理；其次，加入新的PSF 迭代公式，并隨迭代過(guò)程動(dòng)態(tài)更新PSF，同時(shí)對(duì)目標(biāo)估計(jì)圖像去噪，在降低對(duì)PSF初始估計(jì)依賴的同時(shí)增強(qiáng)算法的抗噪性；最后，將估計(jì)的PSF 代入維納濾波以進(jìn)一步提高復(fù)原圖像的質(zhì)量。

1 混合噪聲模型的極大似然算法

通常情況下，傳感器成像平面圖像信號(hào)的數(shù)學(xué)模型可表示為：

式中：yi是成像平面的第i 個(gè)像素值；和分別表示目標(biāo)輻射和背景輻射產(chǎn)生的光電子數(shù)量，均服從泊松分布，其數(shù)學(xué)期望分別為（Hx）i和bi，令總光電子數(shù)量，服從泊松分布；為讀出噪聲，這里假定其為均值r＝0、方差σ2的高斯白噪聲。根據(jù)泊松分布原理有：

傳感器輸出圖像時(shí)，ni受到讀出噪聲的污染，最終的像素值yi服從如下概率分布：

由式（2）（3）可得泊松－高斯混合噪聲模型下圖像像素值的概率密度函數(shù)為：

為了計(jì)算理想目標(biāo)圖像的極大似然估計(jì)，對(duì)式（4）取負(fù)對(duì)數(shù)可得：

將式（5）兩邊對(duì)待估像素值求導(dǎo)數(shù)，可得：

在式（1）確定的成像模型中，并不能確保各像點(diǎn)的像素值都為正值，這是因?yàn)楫?dāng)背景的灰度值接近或等于零時(shí)，由于加性高斯噪聲的污染導(dǎo)致噪聲點(diǎn)的像素值可能為負(fù)值。然而，圖像復(fù)原的目標(biāo)就是得到觀測(cè)目標(biāo)圖像的最優(yōu)非負(fù)估計(jì)，由此可得極大似然算法計(jì)算步驟：

Step1 給定初始值，滿足條件x0≥0。

Step2 由上式給定的初值按下式迭代計(jì)算。

通常情況下，式（9）要經(jīng)過(guò)多次迭代方能得到較好的復(fù)原結(jié)果，但其計(jì)算量非常大，使得傳統(tǒng)極大似然算法效率低下。文獻(xiàn)［4］提出了式（9）的等效近似模型，將參數(shù)b 和σ2引入該式，以表達(dá)泊松－高斯混合噪聲的影響，最終的迭代公式為：

大量實(shí)驗(yàn)表明，用式（10）形成的極大似然迭代算法對(duì)符合泊松高斯混合噪聲模型的退化圖像有較好的重建效果。但是隨著迭代次數(shù)的增加，算法收斂的穩(wěn)定性快速下降，同時(shí)，由于迭代過(guò)程中沒(méi)有考慮噪聲的影響，明顯出現(xiàn)噪聲放大的現(xiàn)象。此外，該算法假設(shè)PSF 已知，且在整個(gè)迭代過(guò)程中保持不變，所以PSF初值的選取嚴(yán)重影響算法收斂的穩(wěn)定性。在PSF未知的情況下，很難得到目標(biāo)圖像的最優(yōu)估計(jì)。本文將針對(duì)Benvenuto算法的不足，加入去噪算法對(duì)含噪圖像進(jìn)行預(yù)處理，并且在算法迭代過(guò)程中對(duì)估計(jì)目標(biāo)圖像進(jìn)行去噪，降低噪聲對(duì)算法穩(wěn)定性的影響；其次，在迭代過(guò)程中加入PSF 的迭代公式，同時(shí)對(duì)其進(jìn)行正性約束和歸一化約束，使得復(fù)原結(jié)果不嚴(yán)重依賴于PSF初值的選取，保證復(fù)原結(jié)果的可靠性；最后，將最終估計(jì)的PSF代入維納濾波，以避免估計(jì)圖像強(qiáng)度集中。

2 圖像去噪算法簡(jiǎn)介與對(duì)比

圖像復(fù)原的實(shí)質(zhì)就是恢復(fù)退化圖像中的高頻信息，但是噪聲也屬圖像中的高頻信息，對(duì)含噪圖像進(jìn)行復(fù)原操作將不可避免地產(chǎn)生噪聲放大的現(xiàn)象。為了解決該問(wèn)題，在復(fù)原操作前對(duì)含噪圖像進(jìn)行去噪預(yù)處理是一種可行的方案。本文簡(jiǎn)要介紹當(dāng)前使用最為廣泛的四種去噪算法（總變分去噪（Total variation，TV）、貝葉斯最小二乘－高斯尺度混合去噪（Bayes least squares－gaussian scale mixtures，BLS－GSM）、非 局 部 均 值 去 噪（Nonlocal means，NLM）以及塊匹配3D 濾波去噪（Block－matching and 3Dfiltering，BM3D）），并對(duì)算法的去噪性能進(jìn)行對(duì)比分析。詳細(xì)的算法原理見(jiàn)文獻(xiàn)［5－8］。

2.1 TV 去 噪

Rudin等［5］于1992年提出了TV 去噪模型，該模型采用能量函數(shù)的最小化問(wèn)題來(lái)對(duì)圖像去噪進(jìn)行建模，并在去噪處理中引入偏微分方程的各向異性擴(kuò)散方程，在抑制圖像噪聲的同時(shí)可以保持邊緣，較好地解決了抑制噪聲和圖像邊緣保持之間的矛盾。

2.2 BLS－GSM 去噪

Portilla等［6］于2003年提出BLS－GSM 去噪算法，該去噪算法主要包括3個(gè)步驟［3－4］：①將圖像在不同尺度和方向進(jìn)行金字塔分解；②對(duì)分解后的每一層進(jìn)行去噪；③用處理后的分層重構(gòu)圖像并進(jìn)行金字塔反變換得到去噪圖像。

2.3 NLM 去 噪

2005年，Buades等［7］提出NLM 去噪算法，該方法結(jié)合圖像中的全局信息綜合考慮圖像噪聲模型（用全部像素的加權(quán)平均來(lái)估計(jì)去噪后的像素值），既消除了傳統(tǒng)鄰域?yàn)V波器中出現(xiàn)的偽影，又可保持邊緣細(xì)節(jié)，同時(shí)利用圖像自身的自相似性能較好地恢復(fù)圖像。

2.4 BM3D去噪

2007年，Dabov等［8］提出了一種理想的圖像去噪方法——BM3D。該方法不僅利用圖像的自相似性和冗余性等信息，而且結(jié)合了變換域的閾值方法，是一種多尺度、非局部的去噪技術(shù)，廣泛應(yīng)用于圖像和視頻的去噪。

2.5 對(duì)比分析

本文主要考察算法的去噪性能和處理效率，因此分別記錄4 種算法的峰值信噪比（Peak signal－to－noise ratio，PSNR）和處理 時(shí)間。為了比較在不同噪聲強(qiáng)度下的去噪性能，對(duì)原始圖像分別施加標(biāo)準(zhǔn)差為10和30的隨機(jī)噪聲，最終的結(jié)果如圖1、圖2及表1、表2所示。

通過(guò)上述圖表可知，以上4種算法在去噪質(zhì)量上都能達(dá)到較好的效果，除TV 算法的PSNR值略低外，其余3 種算法的PSNR 值相近；在處理效率上，由于NLM 算法不僅需要計(jì)算鄰域相似性的權(quán)值，而且圖像中的每一像素點(diǎn)要與所有像素點(diǎn)的鄰域一一比較，所以計(jì)算量非常大，導(dǎo)致效率低下。另外3種算法效率相當(dāng)，其中BM3D和TV 略有優(yōu)勢(shì)。綜上所述，考慮復(fù)原算法的效率和重建質(zhì)量，本文采用BM3D 和TV 算法作為去噪預(yù)處理算法。

圖1 噪聲標(biāo)準(zhǔn)差為10時(shí)的去噪情況Fig.1 Standard deviation when noise is 10

表1 噪聲標(biāo)準(zhǔn)差為10時(shí)的去噪情況Table 1 Denoising results when noise variance is 10

表2 噪聲標(biāo)準(zhǔn)差為30時(shí)的去噪情況Table 2 Denoising results when noise variance is 30

圖2 噪聲標(biāo)準(zhǔn)差為30時(shí)的去噪情況Fig.2 Standard deviation when noise is 30

3 改進(jìn)的極大似然算法

由第2 節(jié)討論可知，Benvenuto等［4］提出的基于混合噪聲模型的極大似然算法存在以下兩方面問(wèn)題：迭代過(guò)程中沒(méi)有考慮噪聲的影響，隨著迭代次數(shù)的增加，穩(wěn)定性變差，噪聲放大現(xiàn)象加?。患俣≒SF已知，且不隨迭代過(guò)程動(dòng)態(tài)更新，導(dǎo)致復(fù)原結(jié)果嚴(yán)重依賴于PSF 初始估計(jì)。在上一節(jié)中提出用相應(yīng)的去噪算法對(duì)含噪圖像作預(yù)處理，并在迭代過(guò)程中對(duì)估計(jì)的目標(biāo)圖像進(jìn)行去噪處理以解決該算法在迭代過(guò)程中對(duì)噪聲敏感的問(wèn)題。在本節(jié)中，將引入新的PSF 迭代公式，以避免該算法對(duì)PSF初始估計(jì)的依賴。

對(duì)于未知參數(shù)b和σ2，采用矩估計(jì)法計(jì)算其估計(jì)值［9］，分別記為和

根據(jù)矩估計(jì)法得到的極大似然算法參數(shù)值，結(jié)合式（10）得到一個(gè)對(duì)稱的PSF迭代公式為：

綜合主成分得分為正的有長(zhǎng)?？h、嵊泗縣、岱山縣、普陀區(qū)、玉環(huán)縣、洞頭縣、平潭縣和東山縣，其中長(zhǎng)海縣、岱山縣和普陀區(qū)的經(jīng)濟(jì)發(fā)展綜合實(shí)力最高；綜合主成分得分為負(fù)的有長(zhǎng)島縣、崇明縣和定海區(qū)，這三海島縣經(jīng)濟(jì)發(fā)展水平欠發(fā)達(dá)，其中定海區(qū)和長(zhǎng)島縣的經(jīng)濟(jì)發(fā)展水平最低。

將式（10）和式（13）進(jìn)行交替迭代，PSF 的值會(huì)隨著迭代過(guò)程動(dòng)態(tài)更新，有效克服了原算法對(duì)PSF初始估計(jì)的依賴，實(shí)現(xiàn)真正的盲解卷積。另外，在交替迭代過(guò)程中，除了加入PSF 正性約束和歸一化約束外，同時(shí)對(duì)目標(biāo)估計(jì)進(jìn)行去噪處理，進(jìn)一步保證了算法的穩(wěn)定性及抗噪性。最后，由式（10）和式（13）進(jìn)行迭代得到的復(fù)原圖像強(qiáng)度分布較集中，這是因?yàn)闃O大似然算法是基于概率統(tǒng)計(jì)的方法，只能最大似然地復(fù)原圖像，而不能很好地保持圖像細(xì)節(jié)。為了克服該缺陷，考慮在算法迭代完成后，將估計(jì)的最終PSF代入維納濾波算法，以進(jìn)一步提高復(fù)原圖像質(zhì)量。改進(jìn)后的算法流程圖如圖3所示。

圖3 改進(jìn)的極大似然算法流程圖Fig.3 Flow chart of improved ML algorithm

4 實(shí)驗(yàn)結(jié)果與分析

為驗(yàn)證本文算法的穩(wěn)定性和抗噪性，分別選取空間點(diǎn)源目標(biāo)和擴(kuò)展目標(biāo)作為模擬數(shù)據(jù)進(jìn)行圖像復(fù)原實(shí)驗(yàn)，圖像大小為256×256，測(cè)試平臺(tái)為Matlab7.14。采用帶有泊松高斯混合噪聲模型的長(zhǎng)曝光大氣湍流模型對(duì)原始圖像進(jìn)行降質(zhì)處理，其中，光學(xué)系統(tǒng)長(zhǎng)曝光OTF為：

式中：λ表示波長(zhǎng)；f 表示望遠(yuǎn)鏡焦距；r0表示大氣相干長(zhǎng)度；表示頻率。

為客觀評(píng)價(jià)算法的復(fù)原性能，引入有參考圖像評(píng)價(jià)指標(biāo)相關(guān)系數(shù)和峰值信噪比以及無(wú)參考圖像評(píng)價(jià)指標(biāo)ImageQ 對(duì)復(fù)原結(jié)果進(jìn)行評(píng)價(jià)。其中，ImageQ 綜合考慮復(fù)原前后圖像棱邊區(qū)局部方差的平均增加量和平坦區(qū)局部方差的減少量來(lái)評(píng)價(jià)棱邊恢復(fù)程度以及圖像偽像的消除程度（原理見(jiàn)文獻(xiàn)［10］）。

峰值信噪比主要評(píng)價(jià)復(fù)原算法的抗噪性，其值越大，說(shuō)明復(fù)原圖像受噪聲的影響越小，其定義為：

相關(guān)系數(shù)通過(guò)計(jì)算復(fù)原圖像與理想圖像對(duì)應(yīng)像點(diǎn)的差值，反映兩者之間的相似度，相關(guān)系數(shù)趨近于1，則表明復(fù)原結(jié)果與理想圖像越接近，復(fù)原效果越好；反之，兩者之間的偏差越大，復(fù)原效果越差。其定義為：

4.1 模擬點(diǎn)源目標(biāo)復(fù)原實(shí)驗(yàn)

模擬雙星圖像的大氣相干長(zhǎng)度r0＝5cm，添加方差σ2＝0.05的高斯白噪聲，理想點(diǎn)源目標(biāo)圖像和模擬退化圖像分別如圖4（a）（b）所示。Benvenuto算法的結(jié)果如圖4（c）所示，本文算法的復(fù)原結(jié)果如圖4（d）（e）所示；評(píng)價(jià)指標(biāo)如表3所示。

表3 雙星圖像評(píng)價(jià)指標(biāo)對(duì)比Table 3 Evaluation scales comparison of double－star image

4.2 模擬擴(kuò)展目標(biāo)復(fù)原實(shí)驗(yàn)

圖4 模擬點(diǎn)源目標(biāo)復(fù)原結(jié)果對(duì)比Fig.4 Result comparison of simulated point source object

表4 擴(kuò)展目標(biāo)圖像評(píng)價(jià)指標(biāo)對(duì)比Table 4 Evaluation scales comparison of satellite image

從圖4、圖5及表3、表4可以發(fā)現(xiàn)，本文算法的各項(xiàng)評(píng)價(jià)指標(biāo)較Benvenuto算法均有不同程度的提高，證明了本文算法的可靠性和有效性。其中，結(jié)合TV 去噪的復(fù)原算法在PSNR 和相關(guān)系數(shù)指標(biāo)上均優(yōu)于結(jié)合BM3D 去噪的復(fù)原算法，建議對(duì)復(fù)原圖像細(xì)節(jié)信息恢復(fù)要求較高的情況下使用；結(jié)合BM3D 去噪的復(fù)原算法在復(fù)原圖像視覺(jué)效果上有明顯優(yōu)勢(shì)，背景噪聲去除效果最好，但同時(shí)也損失了部分目標(biāo)細(xì)節(jié)信息，建議在對(duì)復(fù)原圖像視覺(jué)效果要求較高的情況下使用。

圖5 模擬擴(kuò)展目標(biāo)復(fù)原結(jié)果對(duì)比Fig.5 Result comparison of simulated satellite image

5 結(jié)束語(yǔ)

由于Benvenuto提出的基于泊松高斯混合噪聲模型的極大似然算法存在兩個(gè)方面的缺陷：沒(méi)有充分考慮迭代過(guò)程中的噪聲影像；要求PSF 已知，且在迭代過(guò)程中保持不變。這使得復(fù)原結(jié)果對(duì)噪聲敏感且嚴(yán)重依賴于PSF 的初始估計(jì)。作者針對(duì)以上問(wèn)題，提出結(jié)合噪聲去除的改進(jìn)極大似然算法。首先，采用去噪算法對(duì)含噪圖像進(jìn)行預(yù)處理；其次，加入新的PSF 迭代公式并隨迭代過(guò)程動(dòng)態(tài)更新，同時(shí)對(duì)估計(jì)目標(biāo)圖像進(jìn)行去噪，以減輕噪聲對(duì)迭代過(guò)程的影響；最后，將估計(jì)的PSF代入維納濾波，進(jìn)一步改善復(fù)原質(zhì)量。實(shí)驗(yàn)證明，本文算法有效克服了復(fù)原過(guò)程中的噪聲放大以及對(duì)PSF初始估計(jì)的依賴性，復(fù)原質(zhì)量有明顯改善，可以作為一種有效的復(fù)原方法用于觀測(cè)圖像的高清晰重建。

［1］Ayers C R，Dainty J C.Iterative blind deconvolution method and its application［J］.Optics Letters，1988，13（7）：547－549.

［2］Maccallum B C.Blind deconvolution by simulated annealing［J］.Optics Communication，1990，75（2）：101－105.

［3］Kundur D，Hatzinakos D.Blind image deconvolution［J］.IEEE Single Processing Magazine，1996，13（3）：43－64.

［4］Benvenuto F，Camera A L.The study of an iterative method for the reconstructionof images corrupted by Poisson and Gaussian noise［J］.Inverse Problems，2008，24（3）：35016－35035.

［5］Rudin L I，Osher S，F(xiàn)atemi E.Nonlinear total variation based noise removal algorithms［J］.Physica D：Nonlinear Phenomena，1992，60：259－268.

［6］Portilla J，Strela V，Wainwright M J，et al.Image denoising using Gaussian scale mixtures in the wavelet domain［J］.IEEE Trans Image Processing，2003，12（11）：1338－1351.

［7］Buades A，Coll B，Morel J M.A non local algorithm for image denoising［J］.IEEE Computer Vision and Pattern Recognition，2005，2：60－65.

［8］Dabov K，F(xiàn)oi A，Katkovnik V，et al.Image denoising by sparse 3Dtransform－domain collaborative filtering［J］.IEEE Trans Image Process，2007，16（8）：2080－2095.

［9］周宏潮，朱炬波，王正明.混合泊松－高斯分布模型的參數(shù)估計(jì)［J］.中國(guó)空間科學(xué)技術(shù)，2005，25（2）：1－5.Zhou Hong－chao，Zhu Ju－bo，Wang Zheng－ming.Parametric estimation of mixed Poisson－Gaussian distribution model［J］.Chinese Space Science and Technology，2005，25（2）：1－5.

［10］陳波.自適應(yīng)光學(xué)圖像復(fù)原理論與算法研究［D］.鄭州：解放軍信息工程大學(xué)，2008：63－68.Chen Bo.The theory and algorithms of adaptive optics image restoration［D］.Zhengzhou：PLA Information Engineering University，2008：63－68.