基于自适应梯度先验的旋转模糊图像复原算法

王新春， 王隆隆， 莫波， 刘福祥， 亓贺

(1.北京理工大学 宇航学院, 北京 100081; 2.北京航天自动控制研究所 宇航智能控制技术国家级重点实验室, 北京 100081)

0 引言

弹载红外成像系统随弹体作高速旋转运动时，由于在曝光时间内目标场景与成像传感器之间产生较大的旋转角度，导致所获取的图像存在严重旋转模糊，对后续的目标识别与图像跟踪工作造成极大困扰。因此，有效地快速复原旋转运动模糊图像具有很重要的现实意义。

旋转运动模糊是一种空间非均匀的运动模糊形式，其模糊程度随旋转半径的增大而加剧。针对这一问题，Sawchuk[1]通过极坐标变换，将直角坐标系下的旋转运动模糊变为极坐标系下旋转角方向的平移运动模糊来进行恢复处理，但由于过程中反复进行坐标变换和灰度插值，算法的计算量相对较大且容易丢失原始图像信息。Joshi等[2]、Gupta等[3]对相机动作路径进行分解并构造相应的高维稀疏矩阵，将非均匀运动模糊过程模拟为高维稀疏矩阵和图像向量的乘积。通过迭代求解稀疏矩阵和清晰图像的最大后验概率分布函数，该方法可以取得非常好的复原效果。但是，由于在每一次迭代优化过程中需要存储和处理超高维度的稀疏矩阵，导致该算法计算耗时过长且内存占用过量，很难在工程实际中得到应用。目前，工程上常采用沿模糊路径像素提取的快速复原方法。如Hong等[4]基于Bresenham 路径提取空间不变的模糊像素，避免了坐标变换的复杂运算，同时使用约束最小二乘复原算法提高了图像的恢复效果，缺点是由于假设相邻像素差分极小，算法对图像边缘的复原存在着过正则化的问题。Zheng等[5]在相同算法提取像素的基础上，通过剔除频率域下模糊矩阵的“病态成分”，改进维纳滤波方法，避免了复原过程的过正则化问题，缺点是由于改进的维纳滤波仍需要噪声的先验信息，导致图像的复原结果对噪声变化比较敏感。罗院红等[6]利用最小偏差圆弧插补算法提取模糊路径像素，并通过方向约束的梯度加载复原算法提高算法的抗噪能力，但与约束最小二乘复原算法类似，该算法在图像边缘复原中同样存在着过正则化问题。以上方法存在的缺陷都可认为是由于复原方法对图像先验信息利用不足而导致的。

鉴于正则化方法在利用图像先验信息作为约束项，克服图像复原病态特性的良好应用[7-9]，本文摒弃相邻像素差分极小的假设，引入图像微分先验，提出一种新的基于自适应梯度先验的旋转模糊图像复原算法。该算法在沿运动模糊路径初步复原图像的基础上，通过设计边缘保持滤波器、梯度先验估计方法来估计图像的轮廓信息，并将其作为梯度先验知识进行正则化约束，有效地抑制了复原的过正则化问题，提高了复原结果的抗噪能力。同时，经过推导该算法可以实现频域下求解，大幅提升了运算速度，为其在工程应用提供了可能。

另外，由于当以不同半径圆循环提取图像中像素时，一部分像素不能被任何圆形路径覆盖，导致复原图像中必然出现空穴点。对于此问题，本文改进高密度椒盐噪声的消除算法[10-11]，设计提出了一种新的专门针对Bresenham算法空穴点的自适应中值滤波算法。与文献[4-6]中采用的中值滤波或邻域求均值的方法相比，该算法只是针对像素中的空穴点进行自适应的梯度滤波，避免了图像被过度平滑的问题，同时由于被处理像素数量的减少，插值效率也得以提升。

1 旋转模糊图像退化模型

电荷耦合器件(CCD)图像传感器的成像过程其实就是场景能量在曝光时间内，在CCD像平面的积累过程，而旋转运动造成的图像模糊则相应地可以认为是原始图像的像素灰度沿一系列同心圆弧积分的结果[12]。设原始图像为f(x,y)，旋转模糊图像为g(x,y)，则在曝光时间T内，实际采集的模糊图像g(x,y)与原图像f(x,y)的关系可表述为

(1)

采用极坐标的形式表示为

(2)

将(2)式中r用下标表示，并令s=rωt，ar=rωT，l=rθ，则(2)式可用一维形式表示为

(3)

为了在离散情况下考虑旋转运动模糊图像的复原问题，将(3)式进行离散化，用离散坐标i来表示l，有

(4)

式中：gr(i)和fr(i)分别为半径为r的圆周上像素的模糊灰度值和原始灰度值序列，其周期为Nr，Nr为圆周上像素的总数；ar为模糊长度，用像素个数表示，ar≥1.

定义对应半径为r时的一维点扩散函数，

(5)

则(4)式可用离散卷积的形式统一表示为

(6)

(6)式即为旋转运动模糊图像在半径为r圆周上的退化模型。从模型可以看出，原始图像像素灰度值沿以旋转中心为圆心的一系列同心圆路径在像平面内进行模糊，且同一圆周上各像素灰度的模糊路径是空间不变的，即具有相同的模糊长度和相同的点扩散函数，如图1所示。因此，旋转运动模糊图像的复原问题可以分解成一系列以同心圆路径(模糊路径)为局部域的模糊一致的线性复原问题。

图1 旋转运动模糊路径示意图Fig.1 Blurring path of rotary motion

2 基于自适应梯度先验的复原算法

将(6)式的旋转模糊退化模型写成矢量方程的形式：

g=Hf+η,

(7)

式中：g=[gr(0),gr(1),…,gr(Nr-1)]T为模糊图像灰度值序；f=[fr(0),fr(1),…,fr(Nr-1)]T为原始图像灰度值序列；H为由点扩散函数hr(i)构成的Nr×Nr循环矩阵，其元素值以及阶数随模糊路径的不同而发生变化，且H大多数情况下是不可逆的。这导致图像复原问题通常具有病态特性，不能直接对其逆过程进行求解。

从正则化方法的角度出发，这类逆问题的求解可以通过加入合适的正则项来约束解空间，将上述不适定问题改变为适定问题：

(8)

本文引入图像梯度的稀疏先验信息，定义正则化项和正则化参数为

(9)

式中：w(1)和w(2)分别为图像1阶梯度和2阶梯度的先验序列；λ(1)和λ(2)为其相应的加权比例系数。将(9)式代入到(8)式，有

(10)

对目标函数J(f)关于f求偏导，并令其等于0，

(11)

可得，

(12)

(13)

(14)

(13) 式、(14)式给出了一种基于图像微分先验的模糊图像快速复原算法。然而，由于在实际复原过程中，原始图像微分信息并不是已知的，导致该算法不能被直接使用。这里，本文提出一种新的基于自适应梯度先验的复原算法，它包括图像初步复原—边缘保持滤波—梯度先验估计—图像最终复原4个步骤，具体如下。

2.1 图像初步复原

设w(1)=0,w(2)=0，此时是将相邻像素的差分极小作为空间相关的约束条件，相应的复原方程退化为

(15)

2.2 边缘保持滤波

(16)

式中：k为反馈系数；u为两个相邻像素在转换域之间的距离[14]。

2.3 梯度先验估计

(17)

2.4 图像最终复原

对估计得到的先验项w(1)、w(2)进行傅里叶变换，代入(13)式、(14)式，得到最终图像复原结果。

3 基于查表决策的自适应中值滤波算法

为实时复原旋转运动模糊图像，在使用上述复原算法前，必须首先快速提取出沿模糊路径上各像素点的灰度信息。本文采用Bresenham路径提取方法，通过简单的误差判断对原像平面进行画圆扫描，提取得到位于特定圆弧路径上的像素点及其像素值。该方法可保证提取出的离散圆弧，在位置上是最佳逼近理想圆弧，各对应像素点的最大位置误差不超过0.5个像素，并且由于过程中避开了三角函数、开方等复杂运算，像素提取的在线计算时间大幅减少。

但是，当以不同半径的Bresenham路径提取图像中像素时，一部分像素不能被任何圆形路径覆盖，导致最终复原图像中必然会出现空穴点，如图2(a)所示。这些空穴点处的像素灰度值为0，可以被视为椒盐噪声，采用传统中值滤波或邻域求均值的方法来进行消除。但这一过程同时会使其他像素点的图像信息过度平滑，导致细节丢失。这里，本文提出了一种基于查表决策的梯度自适应中值滤波算法。该算法通过对空穴点像素进行查找、分类，并进行自适应的梯度滤波处理，避免了非空穴点像素的过度平滑问题，同时，由于被处理像素数量的减少，算法的计算量小，运算效率高。

图2 Bresenham路径提取导致的像素空穴点分布图Fig.2 Distribution diagram of missing pixels (holes) caused by Bresenham circular fetching

该算法计算过程主要由裁剪图像空穴点和自适应中值滤波两部分组成。具体如下：

1)遍历查找空穴点像素位置，并根据其邻域内的空穴点分布状况对其进行分类。由于像素空穴点的分布完全取决于Bresenham画圆算法，与图像自身的内容无关，如图2(b)所示。因此，这一过程可以离线计算，计算结果预存于内存单元。实际实验中，为保证邻域内像素之间具有较高的相关性，选用3×3的窗口大小对空穴点像素进行分类，分类结果如图3所示。

图3 图像空穴点分布情况归类统计Fig.3 Distribution statistics of missing pixels (holes)

图3中，黑色方格表示中心空穴点像素，红色方格和橙色方格表示其邻域内必然存在和可能存在的其他空穴点像素。

2)使用查表法，在相应的空穴点像素位置上进行自适应中值滤波操作。算法步骤如表1所示。

表1 自适应中值滤波算法步骤

图4 以空穴点为中心的3×3窗口示意图Fig.4 3×3 window of missing pixels

4 实验结果及分析

4.1 复原效果评价指标

除了直观的视觉图像复原效果，本文采用峰值信噪比(PSNR)、信噪比改善因子(ISNR)及均方根误差(RMSE)作为客观评价指标来评价图像复原前后的效果。PSNR，ISNR及RMSE定义分别为

(18)

(19)

(20)

式中：M、N分别为图像的长和宽。

一般说来，PSNR和ISNR数值越大，RMSE数值越小，则说明图像的复原效果越好，反之则说明复原效果越差。

4.2 仿真图复原效果比较

本文提出的自适应梯度先验复原算法是一种基于模糊路径提取的快速复原方法。首先，为验证其快速性，本文测试对比了该方法与Sawchuk[1]的极坐标转换复原(CTR)方法、Joshi等稀疏矩阵复原方法的图像复原效果和算法耗用时长。测试所用图像由旋转角度为15°的模糊退化图像叠加信噪比为40 dB的高斯白噪声而得，大小为145×145像素，计算机配置为Intel(R) Core(TM) i5-4200H 2.80 GHz，RAM 8.00 GB. 仿真对比结果如表2所示。

从表2可以看出：本文所提出的复原算法取得了不错的复原效果，PSNR 和ISNR 值都比较高，耗用时长也相对较小，基本能达到工程应用的实时性要求；极坐标复原方法进行了两次坐标变换和两次灰度插值，耗用时间相对较长，且由于离散过程极坐标的统一量化，导致图像外缘信息部分丢失，再次插值后复原图像出现一定程度的模糊问题；Joshi等[2]的复原方法通过对目标函数的迭代优化可以取得极好的复原效果，复原图像PSNR超过了30 dB，但由于迭代过程反复处理超高维的稀疏矩阵导致计算量过大，算法处理时间将近1 h ，基本无工程应用价值。

其次，考虑到环境噪声变化对复原算法的影响，本文对模糊退化图像叠加不同信噪比高斯噪声进行仿真实验。实验所用图像大小为435×435像素，模糊角度15°，信噪比分别为40 dB、30 dB、20 dB. 为突出本文算法的鲁棒性和抗噪能力，本文算法与文献[4-5，7]中的复原算法进行了实验对比，对比的图像复原效果及评价指标如表3～表5所示。

表2 不同类别图像复原算法的复原过程及效果比较

Tab.2 Restoration process and effect comparison of different types of image restoration algorithms

表3 不同噪声条件下的各算法图像复原效果

从表3～表5可以看到：采用文献[7]中梯度加载算法得到的复原图像存在明显的振铃效应，且这一现象随信噪比的降低而加剧，在信噪比为20 dB时，复原图像已基本不能识别，ISNR<0 dB；与之相比，由文献[4]中约束最小二乘算法得到的复原图像有了明显改善，PSNR和ISNR值有所提升，但其在低信噪比情况下仍存在一定程度的振铃效应；文献[5]中的维纳滤波算法剔除了频率域下模糊矩阵的“病态成分”，削弱了复原图像的振铃效应，但由于仍需要噪声先验，在给定滤波系数条件下对变化噪声的抑制能力较弱，且复原图像存在过度平滑，细节丢失的问题；本文算法在20 dB、30 dB、40 dB信噪比条件下都相对较好地复原了图像中的细节信息，并且由于合理估计并使用了图像梯度的先验信息，有效地减弱了复原图像的振铃效应，PSNR和ISNR值也明显超出其他3种算法，这说明本文复原算法具有较强的抗噪能力。

表4 不同噪声条件下的复原图像PSNR值

表5 不同噪声条件下的复原图像ISNR值

最后，为检验本文算法对不同旋转角度模糊图像的复原能力，本文设计了一系列不同旋转角度的模糊图像进行复原，并计算不同噪声条件下各复原图像的RMSE. 图5为信噪比分别为20 dB、30 dB、40 dB时复原图像RMSE随旋转角度的变化曲线。

图5 不同噪声条件RMSE随旋转角度的变化曲线Fig.5 RMSE versus rotational angle under different SNRs

一般说来，当RMSE>24时，复原图像的视觉效果非常差[4]。以此为临界条件，可以仿真计算不同噪声条件下，复原算法所能承受的最大旋转角度。从图5可以看出，本文复原算法在20 dB信噪比条件下可以复原最大旋转角度为35°的模糊图像，而当信噪比为40 dB和60 dB时，可复原模糊图像的最大旋转角度将超过60°，基本可覆盖工程应用中绝大部分的高速旋转运动模糊问题。图6为不同复原算法在信噪比为20 dB条件下，复原图像RMSE的对比曲线。从图6可以看出，相比于文献[4-5，7]中的复原算法，本文提出的算法可解决更大旋转角度模糊的图像复原问题，这表明本文算法具有更强的鲁棒性和适应能力。

图6 不同复原算法的RMSE随旋转角度的变化曲线Fig.6 RMSE versus rotational angle

4.3 实际复原效果

在实际模糊图像中，像素的模糊运动模型和干扰噪声与理论分析存在较大差异，为验证本文复原算法的实际可行性，本文将工业CCD相机固定到转速可调的旋转运动模拟平台，以采集的实际模糊图像对算法进行复原测试。算法所需的图像旋转中心可使用激光辅助确定，旋转角度信息则由CCD相机曝光时间和平台转速相乘而得。测试的图像复原效果如图7所示。图7(a)为实际采集的模糊图像，图7(d)为其最终的复原图像，通过这两幅图的对比可以看到，尽管在复原图像的边缘出现一定程度的振铃效应和细节模糊问题(由旋转圆心不稳定、旋转角度存在误差、图像噪声非高斯、模糊模型不匹配等因素导致)，但图形主体部分，如变形的支架、边缘处模糊的旋翼均得到了一定程度的复原，这表明本文提出的算法可用于实际旋转模糊图像的复原，具有一定的工程应用价值。

图7 实际旋转模糊图像复原效果Fig.7 Restored results of actual rotary blurred image

5 结论

针对弹载红外成像系统随弹体作高速旋转运动时图像存在的旋转模糊问题，本文提出了一种基于自适应梯度先验的旋转模糊图像复原算法。该算法首先对初步复原图像进行边缘保持滤波及降噪处理，估计出原始图像的轮廓信息；然后将其作为梯度先验知识进行正则化约束，对模糊图像展开反卷积复原；最后，采用设计提出的自适应中值滤波算法对复原图像的空穴点进行插值，避免复原图像的过度平滑问题。通过仿真实验对比了本文算法与其他复原算法的结果，对比分析结果表明：本文提出的复原算法实时性好，鲁棒性强，可适应更强的噪声条件和较大的旋转角度；同时，图像复原结果在噪声抑制、细节复原及削弱振铃效应方面也具有突出的优越性；对现实中采集的旋转模糊图像，本文算法实现了较大程度的复原，具有重要的工程应用价值。