結合噪聲去除的極大似然圖像復原

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

（1.解放軍信息工程大學 地理空間信息學院，鄭州450052；2.中國科學院 國家天文臺，北京100012；3.中國人民解放軍61175部隊，南京210049）

0 引 言

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

作者針對Benvenuto算法的不足，首先，對圖像進行去噪預處理；其次，加入新的PSF 迭代公式，并隨迭代過程動態更新PSF，同時對目標估計圖像去噪，在降低對PSF初始估計依賴的同時增強算法的抗噪性；最后，將估計的PSF 代入維納濾波以進一步提高復原圖像的質量。

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

通常情況下，傳感器成像平面圖像信號的數學模型可表示為：

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

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

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

為了計算理想目標圖像的極大似然估計，對式（4）取負對數可得：

將式（5）兩邊對待估像素值求導數，可得：

在式（1）確定的成像模型中，并不能確保各像點的像素值都為正值，這是因為當背景的灰度值接近或等于零時，由于加性高斯噪聲的污染導致噪聲點的像素值可能為負值。然而，圖像復原的目標就是得到觀測目標圖像的最優非負估計，由此可得極大似然算法計算步驟：

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

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

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

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

2 圖像去噪算法簡介與對比

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

2.1 TV 去 噪

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

2.2 BLS－GSM 去噪

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

2.3 NLM 去 噪

2005年，Buades等［7］提出NLM 去噪算法，該方法結合圖像中的全局信息綜合考慮圖像噪聲模型（用全部像素的加權平均來估計去噪后的像素值），既消除了傳統鄰域濾波器中出現的偽影，又可保持邊緣細節，同時利用圖像自身的自相似性能較好地恢復圖像。

2.4 BM3D去噪

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

2.5 對比分析

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

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

圖1 噪聲標準差為10時的去噪情況Fig.1 Standard deviation when noise is 10

表1 噪聲標準差為10時的去噪情況Table 1 Denoising results when noise variance is 10

表2 噪聲標準差為30時的去噪情況Table 2 Denoising results when noise variance is 30

圖2 噪聲標準差為30時的去噪情況Fig.2 Standard deviation when noise is 30

3 改進的極大似然算法

由第2 節討論可知，Benvenuto等［4］提出的基于混合噪聲模型的極大似然算法存在以下兩方面問題：迭代過程中沒有考慮噪聲的影響，隨著迭代次數的增加，穩定性變差，噪聲放大現象加劇；假定PSF已知，且不隨迭代過程動態更新，導致復原結果嚴重依賴于PSF 初始估計。在上一節中提出用相應的去噪算法對含噪圖像作預處理，并在迭代過程中對估計的目標圖像進行去噪處理以解決該算法在迭代過程中對噪聲敏感的問題。在本節中，將引入新的PSF 迭代公式，以避免該算法對PSF初始估計的依賴。

對于未知參數b和σ2，采用矩估計法計算其估計值［9］，分別記為和

根據矩估計法得到的極大似然算法參數值，結合式（10）得到一個對稱的PSF迭代公式為：

綜合主成分得分為正的有長海縣、嵊泗縣、岱山縣、普陀區、玉環縣、洞頭縣、平潭縣和東山縣，其中長海縣、岱山縣和普陀區的經濟發展綜合實力最高；綜合主成分得分為負的有長島縣、崇明縣和定海區，這三海島縣經濟發展水平欠發達，其中定海區和長島縣的經濟發展水平最低。

將式（10）和式（13）進行交替迭代，PSF 的值會隨著迭代過程動態更新，有效克服了原算法對PSF初始估計的依賴，實現真正的盲解卷積。另外，在交替迭代過程中，除了加入PSF 正性約束和歸一化約束外，同時對目標估計進行去噪處理，進一步保證了算法的穩定性及抗噪性。最后，由式（10）和式（13）進行迭代得到的復原圖像強度分布較集中，這是因為極大似然算法是基于概率統計的方法，只能最大似然地復原圖像，而不能很好地保持圖像細節。為了克服該缺陷，考慮在算法迭代完成后，將估計的最終PSF代入維納濾波算法，以進一步提高復原圖像質量。改進后的算法流程圖如圖3所示。

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

4 實驗結果與分析

為驗證本文算法的穩定性和抗噪性，分別選取空間點源目標和擴展目標作為模擬數據進行圖像復原實驗，圖像大小為256×256，測試平臺為Matlab7.14。采用帶有泊松高斯混合噪聲模型的長曝光大氣湍流模型對原始圖像進行降質處理，其中，光學系統長曝光OTF為：

式中：λ表示波長；f 表示望遠鏡焦距；r0表示大氣相干長度；表示頻率。

為客觀評價算法的復原性能，引入有參考圖像評價指標相關系數和峰值信噪比以及無參考圖像評價指標ImageQ 對復原結果進行評價。其中，ImageQ 綜合考慮復原前后圖像棱邊區局部方差的平均增加量和平坦區局部方差的減少量來評價棱邊恢復程度以及圖像偽像的消除程度（原理見文獻［10］）。

峰值信噪比主要評價復原算法的抗噪性，其值越大，說明復原圖像受噪聲的影響越小，其定義為：

相關系數通過計算復原圖像與理想圖像對應像點的差值，反映兩者之間的相似度，相關系數趨近于1，則表明復原結果與理想圖像越接近，復原效果越好；反之，兩者之間的偏差越大，復原效果越差。其定義為：

4.1 模擬點源目標復原實驗

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

表3 雙星圖像評價指標對比Table 3 Evaluation scales comparison of double－star image

4.2 模擬擴展目標復原實驗

圖4 模擬點源目標復原結果對比Fig.4 Result comparison of simulated point source object

表4 擴展目標圖像評價指標對比Table 4 Evaluation scales comparison of satellite image

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

圖5 模擬擴展目標復原結果對比Fig.5 Result comparison of simulated satellite image

5 結束語

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

［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，Fatemi 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，Foi A，Katkovnik V，et al.Image denoising by sparse 3Dtransform－domain collaborative filtering［J］.IEEE Trans Image Process，2007，16（8）：2080－2095.

［9］周宏潮，朱炬波，王正明.混合泊松－高斯分布模型的參數估計［J］.中國空間科學技術，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］陳波.自適應光學圖像復原理論與算法研究［D］.鄭州：解放軍信息工程大學，2008：63－68.Chen Bo.The theory and algorithms of adaptive optics image restoration［D］.Zhengzhou：PLA Information Engineering University，2008：63－68.