结合噪声去除的极大似然图像复原

姜 超，耿则勋，刘立勇，潘映峰

（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.