A non-local prior of infrared motion blurred image restoration method

成像模糊处理方法以观测图像来复原真实图像，是探测成像的逆计算过程^[1]，是弱小红外目标检测识别的预处理过程。图像模糊复原研究是计算机视觉的研究方向之一^[2-12]，尤其在运动模糊图像复原方面^[7-11]。图像运动模糊包括图像局部运动模糊(成像场景中的目标运动引起的图像模糊)和图像全局运动模糊(成像设备运动引起的图像模糊)两方面。局部运动模糊是成像过程中，成像场景中的目标运动产生的图像局部模糊；全局运动模糊则是由成像设备运动产生的图像整体模糊。图像的局部运动模糊复原研究以分析目标的运动信息，或者计算相邻图像帧之间的光流信息为主，获取目标运行信息，从而复原局部运动模糊图像；图像的全局运动模糊复原则研究成像设备的运动轨迹，通过估计成像全局运动信息，建立运动模糊成像的盲反卷积复原算法模型，从而依次复原图像和运动模糊核^[13-16]。运动模糊的盲反卷积复原算法模型的核心是盲反卷积先验约束的设计，本文作者在此基础上，通过研究分析红外运动模糊成像特征，建立了一种针对红外运动模糊成像特征的盲反卷积复原算法模型。

红外探测成像具有全天候、抗干扰等优点。与可见光探测成像相比，红外探测成像器件的热响应时间远大于电荷耦合器件的曝光时间，因此，红外图像的运动模糊现象更严重^[17]。

基于最大后验模型的正则化方法是一类经典的图像复原方法，该方法在贝叶斯公式基础上，利用图像结构信息，分别构建数据拟合项和正则化先验约束项^[18]。盲反卷积图像复原算法在最大后验方法的基础上，对成像模糊核进行设计，交替迭代过程中，约束复原图像与图像模糊核^[16]。非局部均值(non-local, NL)设计方法利用图像的非局部性，进行图像正则化先验建模^[19]。拉普拉斯分布是一种典型的稀疏先验建模方法之一^[1]。

通过红外成像特性分析和成像载体的运动分析，针对红外运动模糊图像复原问题，首先分析运动模糊成像特征，利用运动信息，计算初始图像模糊核，其次设计了一种图像非局部稀疏先验约束项，并结合图像模糊核的拉普拉斯先验约束，最终提出了一种盲反卷积框架下的图像复原方法。与现有红外运动模糊图像复原方法相比，该方法充分利用了成像特征，包括图像运动信息和载体运动条件下的红外图像成像特征(目标能量沿运动方向拖尾)，有利于进一步剖析运动模糊图像的复原。该方法以红外图像的运动特征为研究重点，具有实际应用意义。

探测成像热响应时间内系统平台的运动是产生图像退化的主要因素，具体表现为红外图像的运动模糊。运动导致成像器件接收的目标/场景高频信息随平台的运动方向和运动幅值分布，图像呈现出运动拖尾现象，尤其是热响应时间较长的长波非制冷成像器件。长波非制冷凝视红外探测成像系统的成像环境温度波动较大，在常温下，典型目标/场景的温度分布相似，导致非制冷长波红外图像的目标/场景对比度较小，与制冷型红外图像相比，成像质量更差；非制冷成像器件通过线扫方式成像，探测器件的热响应时间较长，且非制冷型红外成像系统作用距离较近，系统多为音速或亚音速，此时的大气热效应对非制冷红外成像系统影响较小。综上分析，运动因素是影响非制冷红外成像系统图像退化的主要因素，且红外图像运动模糊退化严重，需要深入分析红外图像运动模糊特征，建立相应的运动模糊图像退化复原模型，复原红外图像。本文中主要研究非制冷长波红外图像，其热响应时间较长，图像的运动模糊幅值较大，运动像移大。

红外探测成像过程中，红外成像器件的热响应时间是固定的和周期性的，在该时间内，平台的运动方向和运动幅值会导致视场内目标/场景辐射信息在器件单元的移动积累，从而呈现出图像高频信息随运动拖尾模糊。红外成像链路图如图 1所示。

Figure 1.  Infrared imaging link diagram

与采用局部特性设计的全变分^[16](total variation, TV)正则化先验不同，NL先验强调图像的非局部性，先验设计的图像区域比TV先验范围更大，且NL先验非常稀疏。NL先验能表示更多的图像信息，并且在求解过程中并没有增加太多的计算复杂度。NL稀疏先验约束U(f)表示为：

式中：ρ>1表示约束范数因子，求解过程中，取ρ=2；N(f)表示非局部滤波算子；f为原始图像。

平台运动引起红外图像退化的特征是：热响应时间内，红外辐射信息较强的目标等物体表现为高频信息，辐射信息变化不大的场景等物体为低频信息，高频信息受平台运动影响较大，运动模糊程度明显；低频信息受平台运动影响较小，运动模糊不明显。并且，成像光学系统的非均匀性放大了探测成像阵列接收的辐射信息分布不均匀，成像信息表现出不均匀的性质，平台运动对探测成像中心影响较大，成像四周影响较小，如图 2所示。

Figure 2.  a—the 35th motion blur image b—the 36th motion blur image

图 2为选中视频相邻两帧运动模糊图像的高频信息。图中的物体运动模糊现象非常严重，且运动模糊拖尾与运动幅值和方向密切相关。运动模糊图像呈现出随运动方向不同程度的运动模糊拖尾现象；同时，图像右下角的汉字和数字间隙小，运动模糊拖尾相互干扰。

图 3是红外成像系统在外界运动影响因素下的航向与俯仰方向连续多帧运动成像相关的角度运动分析图(纵向为运动角度，航向为图像帧)。从图 3可以初步计算出，图 2a在航向方向相邻图像帧之间的运动为-0.15°，俯仰方向相邻图像帧之间的运动为+0.19°，大约24个像素；图 2b在航向方向相邻图像帧之间运动为-0.18°，俯仰方向相邻图像帧之间运动为+0.30°，大约35个像素。

Figure 3.  Motion angle diagram in pitch and course direction

从图 3可以看出，平台红外成像系统成像模糊程度与平台运动密切相关，同一成像环境下，运动模糊是影响成像的主要因素，平台运动数据直接反映图像的模糊程度。利用运动数据信息，估计初始运动模糊核，是一种改善捷联系统成像模糊的有效方法。

接下来，推导平台运动条件下图像的模糊核。根据红外成像链路理论，在热响应时间内，对接收的辐射能量信息进行积分表示，即：

式中：$ \left(i^{\prime}, j^{\prime}\right)$表示像素点；T表示热响应时间；$ \boldsymbol{f}\left(i^{\prime}, j^{\prime}\right)$为原始图像；$ \boldsymbol{h}\left(i^{\prime}, j^{\prime} ; x, y, \boldsymbol{f}\left(i^{\prime}, j^{\prime}\right)\right)$表示成像等效能量转换过程，在成像热响应时间内，为有限值，因此绝对可积；$ \boldsymbol{g}\left(i^{\prime}, j^{\prime}\right)$表示观测图像，故上式可写成：

在探测成像系统的热响应时间内，平台运动假设近似为直线运动，并将直线运动分解到直角坐标系内^[20]。对上式左右两端进行傅里叶变换，将运动模糊图像在频域上表示为：

式中：N (u, v) 为噪声函数在频域内的表示；$ \boldsymbol{H}(u, v)$是频域上的等效模糊核函数，$ \boldsymbol{H}(u, v)=\int_0^T \exp \{-\mathrm{j} 2 \pi \times \left.\left[u i_0(t)+v j_0(t)\right]\right\} \mathrm{d} t$；$ \boldsymbol{G}(u, v)$为$ \boldsymbol{g}\left(i^{\prime}, j^{\prime}\right)$的频域观测图像表示；$ \boldsymbol{F}(u, v)$为$ \boldsymbol{f}\left(i^{\prime}, j^{\prime}\right)$原始图像的频域表示；i₀(t)和j₀(t)为x和y方向的像移运动。

等效模糊核函数与探测成像系统的热响应时间和系统平台运动相关，假设热响应时间T内，图像的像移幅值为L，像移方向角为θ，任意热响应时间t∈[0, T]内，x和y方向的像移运动i₀(t)和j₀(t)表示为：

式中：L_x和L_y为x和y方向的运动距离。因此，运动模糊核在像移幅值和方向条件下表示为：

假设运动模糊核大小为N×M，在像移幅值和方向条件下，对上式进行化简，图像模糊核的离散化表示为：

非局部滤波算子由迭代更新过程中的中间图像$ \tilde{f}$中任意以(k, l)像素为中心的图像块进行定义，N(f)具体定义为：

式中：$ \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, j^{\prime}\right)$和$ \omega\left(k, l ; i^{\prime}, j^{\prime}\right)$分别表示像素点(i′, j′)与中心像素点(k, l)的相似度度量和归一化相似度，相似度度量为非负，即0≤ω(k, l; i′, j′)≤1。本文中在定义非局部滤波算子时，利用了平台的运动参数(s_x, s_y)，将上式中的相似度度量$ \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, j^{\prime}\right)$分解表示为图像行方向与列方向的相似度度量权重函数，其数学表达式为：$ \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, j^{\prime}\right)=\sigma_x \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, \left.j^{\prime}\right)+\sigma_y \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, j^{\prime}\right)\right.$，σ_x和σ_y分别为航向方向和俯仰方向的运动幅度归一化参数；f (i′, j′)表示原始图像中心像素点(k, l)的邻域图像；相似度度量权重函数$ \omega_{\mathrm{NL}}\left(k, l ; i^{\prime}, j^{\prime}\right)$表示图像灰度值与空间位置的函数，其表达式为：

式中：第1项是图像灰度值算子，表示以(k, l)为中心像素的q×q图像块与以(i′, j′)为中心的q×q图像块在算子R作用下的指数分布；R_{k, l}是以(k, l)为中心像素点的图像块提取算子；a>0是核标准差参数；σ′是以(k, l)为中心和以(i′, j′)为中心的q×q图像块的图像灰度等级差异参数；R算子、参数a联合控制图像块图像灰度值的指数分布权重；第2项中的D(·)是关于图像块之间的空间距离函数，D(·)具有单调不增性质，通常取高斯函数、常数等，本文中在设计过程中对D函数取常数值处理。

将式(8)的非局部权重设计代入式(1)的非局部先验约束中，U(f)表示为：

式中：I和S分别表示单位矩阵和非局部稀疏权重矩阵。权重矩阵S是针对中间图像的块结构相似性进行设计的，在结构上，S满足稀疏性。实现过程中，首先对图像进行分块区域处理，根据图像分块，设计图像非局部稀疏先验约束，并将输入图像g作为初始值f₀。

结合图像的数据拟合、非局部稀疏先验约束设计以及图像模糊核的拉普拉斯稀疏先验约束，在盲反卷积框架下，得到红外图像运动模糊退化复原算法模型，其目标泛函表示为：

式中：f表示原始图像(需求解的)；g表示观测图像；h表示模糊核函数；U(f)表示图像先验约束；Q(h)表示图像模糊核约束；α₄和α₅为超参数；Ψ(h_{(i′, j′)})为拉普拉斯先验分布。

综上分析，利用分离变量法，将目标泛函中不同范数约束项进行分离，采用共轭梯度迭代法、增广拉格朗日乘子法(augmented Lagrange method, ALM)方法以及收缩法，得到图像的非局部稀疏先验约束算法模型下的运动模糊复原图像f ^*。改进的非局部稀疏先验约束算法流程如下(输入是运动模糊图像；输出是退化复原图像)：(a)计算图像的初始模糊核 h₀；利用系统平台运动信息，包括像移幅值L和像移方向角θ，根据式(7)计算图像运动模糊核初始值 H₀；对H₀进行快速傅里叶逆变换(inverse fast Fourier transform, IFFT) h₀=IFFT(H₀)，得到h₀；(b)对输入图像进行分块处理，设计图像的非局部稀疏先验约束项；设计图像非局部稀疏先验约束项时，利用相似度度量权重函数(见式(9))构造稀疏先验；(c)采用分离变量法，求解复原图像 F，并进行快速傅里叶逆变换，得到f_k和输入图像的模糊核 h_k；(d)固定图像模糊核h_k，采用共轭梯度法迭代求解复原图像F_k；(e)固定图像F_k，采用ALM法迭代求解模糊核 H_k+1；反复执行第(d)和第(e)步，直至满足迭代次数要求时，终止迭代，将迭代复原的图像 F和图像模糊核H进行傅里叶逆变换，得到最终的f和 h。

式(11)中，f和h均为未知量，在数学上具有无数组解，属于病态数学模型。采用交替迭代方法来求解式(11)，依次迭代求解 f和h。交替迭代方法的核心思想是先固定其中的第1个变量 h₀，迭代求解第2个变量f_k；然后固定求解的变量f_k，迭代求解第1个变量h_k，如此交替进行。h₀表示模糊核的初始值，来源于式(7)，即h₀ =IFT(H)；f则是以输入图像为初始值。由于式(11)中包含了对图像的非局部稀疏先验的L₂范数约束和对图像模糊核的L₁范数约束，因此不能直接对L(f, h)求偏导数。采用变量分离的方法，将L(f, h)分解成下述形式：

固定变量h_k，直接对 f进行求导数处理，并采用共轭梯度法求解f。关于L₁(f)的共轭梯度法表示为：

通过固定h_k变量，利用共轭梯度法，迭代求解图像 F。最后进行快速傅里叶逆变换f =IFFT(F)，得到迭代求解的f。

式(12)的L₂(h)中，在迭代求解F变量条件下，采用ALM方法，首先将L₂(h)替换成：

将约束条件下的泛函按照拉格朗日乘子法写成无约束泛函，即：

式中：β₄为超参数；a_h和υ_h是引入的关于h的中间变量。对式(16)进行分离变量处理，表示为：

为求解式(17)中关于h的泛函，对其进行求导数，并将其置为0，采用迭代的方式表示为：

为求解式(17)中关于υ_h的泛函，采用收缩的方式迭代求解[υ_{h, k+1}]，表示为：

中间变量a_{h, k+1}的迭代方式为：

令$ \boldsymbol{v}_{\boldsymbol{h}, 0}=0, \boldsymbol{a}_{\boldsymbol{h}, 0}=0, k=0$, 联立式(18)~式(20)，最终迭代求解 h。

在仿真实验设计中，模糊核的大小均为25 pixel×25 pixel，图像的分块大小分别为64 pixel×64 pixel，迭代次数为3，α₄=1, α₅=1, β₄=10。

实验结果分析方面，分别从主客观对比了TV正则化先验方法^[16]、双边全变分(bilateral total variation, BTV)正则化先验方法^[21]、Kotera’s方法和本文作者提出的基于运动信息的非局部稀疏先验约束方法。其中的400 pixel×400 pixel大小图像详见图 4~图 7。图 4~图 7是针对连续不同幅值的运动模糊红外图像。本文中提出的方法复原效果明显，复原结果图像中的目标较清晰。

Figure 4.  TV restoration result of Fig. 2a and Fig. 2b

Figure 5.  BTV restoration result of Fig. 2a and Fig. 2b

Figure 6.  Kotera's restoration result of Fig. 2a and Fig. 2b

Figure 7.  Result of non-local sparse prior restoration of motion information of Fig. 2a and Fig. 2b

本文中采用对图像像素概率模糊检测(cumulative probability of blur detection, CPBD)^[22]、图像结构相似度(structural similarity, SSIM)、峰值信噪比(peak signal-to-noise ratio, PSNR)方法量化复原结果，从客观角度分析评价了TV先验方法、BTV先验方法、Kotera’s方法以及本文中提出的方法。CPBD方法逐像素建立图像模糊概率模型，是一种非参考图像评估方法，该概率越大越好；SSIM表示图像之间的相似程度；PSNR通过统计计算复原图像与原始图像之间的像素差。上述3种客观评价方法量化对比见表 1和表 2。

从表 1和表 2可以看出，本文中提出的运动模糊图像复原方法的复原图像的客观评价较高，尤其是针对运动幅值较大的模糊图像，其中PSNR提高接近8%，结合图 4~图 7，基于运动信息的非局部稀疏先验图像复原方法针对性较强，对运动模糊退化图像具有明显的复原效果。

图 8是图 2a和图 2b的运动模糊核。可以看出，图 2b的运动幅值比图 2a更大。

Figure 8.  Image motion blur kernels of Fig. 2a and Fig. 2b

本文中以红外成像的运动特性为出发点，通过分析平台的运动特性、非制冷红外运动成像特性，通过推导分析运动模糊核、设计运动信息下的非局部稀疏先验约束，提出了一种盲反卷积框架下的非局部稀疏先验约束算法模型，用于处理运动模糊红外图像的退化复原，为红外成像系统应用研究提供了支撑。通过主观和客观对比，所提出的红外图像复原方法针对性较强，并且对于运动幅值较大的图像，从视觉上可以看出，复原效果明显；采用CPBD、SSIM以及PSNR的客观图像评价方法，尤其是运动幅值较大的红外图像，PSNR提高了接近8%。本文中提出的方法针对性较强，能有效地处理运动幅值较大的红外图像，具有实用意义。