Removal of image motion blur caused by distortion in TDI remote-sensing image

近年来，海洋遥感仪器的发展趋向于大视场和高分辨率。文中研究的新型海洋探测器采用离轴三反射式光学系统和面阵TDI推扫成像。离轴三反射式光学系统的光路无中间遮拦，提高了可用光学视场，因此被广泛应用于大视场相机。而引入TDI技术的面阵推扫成像与单线阵推扫成像相比，前者可以在积分时间内对同一景物进行多次曝光，通过多级信号的累加提高图像的信噪比，从而满足高分辨率成像的需求。然而大视场相机的离轴三反系统存在较大的光学畸变，在TDI成像过程中,空间变化的畸变导致图像累加之后出现像移模糊。

目前，关于抑制畸变对TDI成像影响的研究，大多是对成像系统进行改进, 李富强等^[1]通过弧形拼接探测器焦平面，补偿畸变引起的像移。任秉文等^[2]推导了畸变和像移的数学关系，从光学系统的设计上，给出了相关参数的设计约束。然而对于不同的探测器而言，成像系统的结构各不相同，这类方法的通用性较差，因此，研究相应的图像复原算法是十分必要的。

图像复原是根据成像原理，由模糊图像恢复出清晰图像的过程。根据退化函数是否已知，图像复原可以分为非盲复原和盲复原。在实际应用中，退化函数往往无法直接精确获取，因此近年来许多学者开展了对图像盲复原算法的研究。Fregus等^[3]基于自然图像的梯度分布，利用变分贝叶斯方法估计模糊核，采用Richardson-Lucy算法反卷积清晰图像，但是变分贝叶斯运算复杂效率低，另外Richardson-Lucy复原容易产生振铃效应；Krishnan等^[4]在图像复原模型中引入一种尺度不变的正则化项，较完整的保留了图像的高频分量，但是模糊核约束项较为单一，复原图像细节处存在残留模糊；李正周等^[5]基于空间目标图像梯度的稀疏特性，采用L₀范数提取显著边缘估计模糊核，较好的保留了复原图像的纹理信息，但是仅用L₀范数对图像梯度进行约束，不利于提取细小的图像结构，影响模糊核的估计精度；Shan等^[6]利用空间随机分布的噪声模型和局部区域的平滑约束改善了图像的振铃效应，但是复原图像整体不够清晰。洪汉玉等^[7]通过提取退化图像的多尺度形态学梯度，利用加权L₁范数测度筛选出适合模糊核估计的区域，使梯度信息的获取更准确，但是该方法不适用于图像局部区域梯度信息匮乏且模糊是空间变化的情况。Zhang等^[8]利用相对全变差模型提取潜像的显著边缘估计模糊核，虽然模糊核精度较高，但是耗时较长。以上这些方法利用图像梯度信息作为先验信息，对于空间变化的模糊，不适用于所有的局部区域，尤其是遥感图像往往存在梯度信息匮乏的区域，如海洋和荒漠，此时结构特征较难提取，单一的稀疏约束将影响模糊核的估计精度。

文中针对现有去模糊方法先验约束单一和遥感图像梯度信息匮乏影响模糊核估计精度的问题，提出根据图像退化原理和成像先验信息预估初始模糊核，将其作为先验信息构建图像半盲复原模型细化模糊核。另外提出一种多向各异性全变差模型提取图像结构，保留了图像的细小结构有助于提取显著边缘。对比其他方法，文中方法提高了模糊核的估计精度，具有更好的去模糊效果。

光学系统的畸变是指视场位置不同的物点在焦平面上成像的放大率不一致，导致图像变形，对于一般的成像光学系统，畸变只影响成像的形状，不影响图像的清晰度^[2]。但是当采用TDI成像方式时，由于畸变和视场位置有关，目标在各级线阵不同位置的成像畸变都不一致，导致各级线阵的成像存在相对像移，累加之后出现模糊，如图1所示，由此可见畸变导致的TDI像移模糊具有空间变化性。

Figure 1.  Schematic diagram of TDI image motion blur caused by optical distortion

由文中分析可知，不同积分时刻的图像产生了相对像移，相当于图像发生了相对运动，根据运动模糊理论^[9]，此时的PSF近似于像点位置的概率密度函数(Probability Density Function, PDF)，当位移是二维的且平行于焦平面时，PDF可以表示为速度的倒数，因此PSF的表达式为：

式中：$T$为曝光时间；$N$为积分级数；${{\boldsymbol{V}}_i}(x,y)$为像点$(x,y)$在第$i$次积分时刻的像移速度矢量，通过对像移路径求导可得。初始PSF的求解步骤如下：

(1)计算N个积分时刻下像点的实际位置；首先采样理想位置，以第一级线阵和最后一级线阵的积分时刻为成像的始末，采样像点在每级线阵成像时的位置，固定初始像点位置$({x_1},{y_1})$，则第$n$个积分时刻的理想像点位置为：

式中：$w$为像元的宽和高。

(2)计算像点的实际位置；通过CODE V仿真的离轴三反射光学系统获取畸变数据，采用三次多项式拟合畸变模型^[10]。将公式(3)采样的理想位置代入模型，得到实际坐标序列，拟合为关于t的曲线，如图2(a)所示。

Figure 2.  Solution process of initial PSF

(3)对像移路径曲线求导，得到不同积分时刻的像移速度矢量，代入公式(1)得到PSF；以像素为单位离散化PSF函数，如图2(b)所示。

(4)将离散的PSF写成矩阵形式；以初始像点为中心位置，根据像移方向和位置，将PSF值填入矩阵中，对矩阵进行归一化和中心化，结果如图2(c)所示。

由于遥感图像数据量较大，直接用空间变化的复原方法需要很大的存储量，因此选择分块复原的方式。推扫图像的一行即一帧，不考虑其他引起成像模糊的因素时，帧间模糊相同，PSF的变化性体现在帧内各像点。通过计算帧内各像点PSF的相似度划分区域^[11]。由于畸变在垂直推扫方向上呈现连续渐变性，因此推及所有帧，图像被划分为连续的若干矩形子块。对子块内像点的PSF求和取平均，作为该子块的初始模糊核。

上述采样过程的离散性和拟合模型的固有误差导致初始模糊核含有较多噪声，需要对其进一步优化。文中提出一种图像半盲复原算法，将初始模糊核作为输入，求取粗略的中间潜像，利用潜像的显著边缘进一步细化模糊核。同时将初始模糊核用于构建正则项约束，增强结果的可靠性。采用基于图像金字塔分解的多尺度迭代框架，实现模糊核估计的逐步细化^[12]，整体流程如图3所示。

Figure 3.  Framework of semi-blind algorithm for estimating blur kernel

第一次迭代时的中间潜像由初始模糊核恢复，后续迭代过程中的潜像由中间模糊核恢复。中间潜像的估计函数采用全变分(TV)作为正则项，可以保留中间潜像的丰富纹理，表达式如下：

式中：${\;\beta _3}$为TV正则项的权重，取值为0.1；$\nabla {\rm{ = }}\left( {{\partial _x},{\partial _y}} \right)$为$x$和$y$方向一阶导数。采用迭代重加权最小二乘法(Iterative reweighed least square, IRLS)^[13]求解，外层IRLS迭代2次，新的权重根据更新的${\boldsymbol{I}}$计算，内层共轭梯度算法(Conjugate Gradient, CG)迭代5次。

图像的显著边缘信息有助于模糊核的估计，而显著边缘的提取来自可靠的图像结构，常用相对全变差模型分离图像结构^[12]，定义如下：

式中：${{\boldsymbol{I}}_s}$为分离的图像结构；${{\boldsymbol{I}}^*}$为中间潜像；$\lambda $为平滑系数，经验值取0.02；该式中一阶导数算子只考虑了$x$和$y$方向，忽略了其他方向的细小结构。文中提出一种多方向权重异性的全变差模型，将两条对角线方向的梯度变化也考虑在内，另外在每个方向开辟方形窗口，以窗口区域内的像素共同确定权值，改进后的模型定义如下：

式中：${\partial _{{{xy}}}}$和${\partial _{{{yx}}}}$分别为两个对角线方向的一阶导数算子；${w^ * }{\rm{ = }}\left\{ {{w_1},{w_2},{w_3},{w_4}} \right\}$表示四个方向的权重，定义如下：

式中：${\boldsymbol{r}}$为权值算子，定义如下：

式中：$p$为目标像素索引；$q$为以$p$为中心的正方形区域$R(p)$内的像素索引；$\Theta {\rm{ = }}\left\{ {{\partial _{{x}}},{\partial _{{y}}},{\partial _{{{xy}}}},{\partial _{{{yx}}}}} \right\}$；$\varepsilon $用于避免分母为0的情况，同时保证结构的平滑性，经验值取0.4。公式(5)是凸优化问题，采用快速傅里叶变换(FFT)直接求解：

与传统方法相比，改进后的图像结构提取模型能够保留更多的细小结构，如图4所示。

Figure 4.  Comparison of image structure extraction

得到图像结构后，用shock滤波器进行锐化^[13]：

式中：${\widetilde {\boldsymbol{I}}_{{\rm{s}}t}}$为第 $t$次迭代时的锐化图像；$\Delta $为拉普拉斯算子；$ \mathrm{sign}(\cdot )$为符号函数；最后筛选${\widetilde {\boldsymbol{I}}_{\rm{s}}}$中的显著边缘用于模糊核估计，表达式如下：

式中：${\widetilde {\boldsymbol{I}}_{\rm{s}}}$为锐化后的图像结构；$M(D,\tau )$为单位二值掩膜函数^[12-13]，定义如下：

式中：$D{\rm{ = }}{\left\| {\nabla {{\widetilde {\boldsymbol{I}}}_{\rm{s}}}} \right\|_2}$；$\tau $的取值要保证至少有$\sqrt {{k^2}WH} $个像素点被选中；${k^2}$和$W \times H$分别为模糊核和输入图像的大小。

为了保证模糊核的稀疏性和平滑性，选取L₁范数正则项$||{\boldsymbol{H}}|{|_1}$和梯度约束的正则项$||\nabla {\boldsymbol{H}}||_2^2$，利用显著边缘${\boldsymbol{\tilde I}}$的梯度信息$\nabla {\boldsymbol{I}}$构建保真项，模糊核估计模型如下：

式中：$\nabla {\boldsymbol{G}}$为模糊图像的梯度；${\;\beta _1}$和${\;\beta _2}$为正则项系数；该模型的正则项仅从模糊核的特性方面进行了约束，迭代过程中依赖于图像的信息，应该利用更多关于探测器的先验信息进行约束^[14-15]。由第2节分析可知，初始模糊核是根据探测器的畸变数据和TDI成像参数求得，先验条件可以通过光学仿真和硬件测试获取。因此将初始模糊核作为先验核，加入模糊核估计模型作为先验约束，使估计结果贴合成像退化过程，避免过度依赖图像梯度信息。初始模糊核的约束项设计为$||{{\boldsymbol{H}}_{{0}}} - {\boldsymbol{H}}||_2^2$，最终的模糊核优化估计模型为：

利用Split-Bregman算法求解上式，引入辅助变量${{\boldsymbol{b}}_{{0}}}$和惩罚变量${{\boldsymbol{t}}_{{0}}}$，公式(13)表示为：

式中：${{\boldsymbol{b}}_{{0}}} = {\boldsymbol{H}}$;$\eta $为惩罚项的正参数。将公式(14)分解成以下几个表达式：

利用变分原理和快速傅里叶变换解公式(17)：

式中：$ F(\cdot )$和$ {F}^{-1}(\cdot )$分别表示傅里叶变换和逆变换；${{\boldsymbol{I}}^{\rm{T}}}$为${\boldsymbol{I}}$的共轭转置。采用软阈值法求解公式(15)^[16]。考虑到模糊核的物理意义，迭代过程中需对其阈值进行动态约束：

模糊核已知后，采用非盲反卷积模型去除图像模糊，非盲反卷积的病态性比盲复原小很多，采用L_P正则项可以取得较好的效果^[4,17]。图像复原模型如下：

式中：$p$取值0.8；该函数非凸，故采用IRLS算法求解。

为了更好的抑制边界振铃效应，迭代复原之前，对模糊图像进行边界延拓和渐变处理^[18]。复原后的子块去除拓展区域再拼接得到最终的复原图像。

文中研究的探测器光学系统视场14.4°×0.84°，前置望远镜采用离轴三反射式结构，边缘视场存在较大的畸变。焦平面像元大小30 μm×30 μm，线阵感光像元数为2 048。由于TDI级数越多畸变像移模糊越明显，为了更好的观察算法效果，TDI级数设为96。根据文中提出的分块方式，TDI图像被分为10个子区域，为了简化验证过程，选取边缘视场的两个相邻子区域验证算法，推扫帧数设为300，实验的图像大小为300×352，模糊核大小为15×15，用于对比的其他方法模糊核设置为同样大小，分块方式保持一致。

为了验证文中算法的有效性，使用三组仿真图像进行测试，并与Zhang的算法^[8]、Krishnan的算法^[4]作比较，结果如图5所示，从左往右依次是模糊图像、Zhang算法、Krishnan算法和文中方法的去模糊图像。红框中的是振铃波纹，绿框中的是下方放大的细节。

Figure 5.  Deblurring results of three images by different methods

在主观视觉上，三种复原方法都能取得较好的去模糊效果。文中复原后的清晰程度和Zhang的相当，但是Zhang的复原图像边缘出现了振铃；三组实验中，Krishnan的方法均消除了振铃效应，但是通过放大的细节图可以发现，Krishnan的复原图像存在残留模糊，对比之下，文中方法的复原结果不仅清晰自然，且无振铃效应。

在客观评价上，使用无参考评价指标拉普拉斯梯度值(Laplacian Gradient, LG)和离散余弦变换值(Discrete Cosine Transform, DCT)作为评价指标。LG是图像空域参数，利用拉普拉斯算子计算图像各像素点的梯度，算法简洁快速。DCT是图像频域参数，对图像清晰度变化敏感，通过计算图像高频分量的多少评价模糊程度，两种评价指标都是值越大代表图像越清晰，计算公式如下：

式中：${\boldsymbol{I}}$为待评价图像，其尺寸为$M \times N$；${\boldsymbol{C}}(i,j)$为离散余弦变换的系数矩阵。客观评价结果如图6所示，可以直观的看到文中方法复原后的三幅图像评价指标都得到提高，其中LG值平均提升了67.5%，DCT值平均提升了14.2%。三幅图像的指标都优于Zhang和Krishnan，客观评价与主观分析相符，证明文中方法有效。

Figure 6.  Comparison of objective evaluation

文中方法在估计模糊核时引入先验模糊核的约束项，在图像结构信息不足时，能够取得比其他方法更稳定的模糊核估计。选取如图7(a)所示的港口沿岸水面目标红外图像进行验证，由于水面多为平坦区域，大结构信息较少，此时Zhang估计的模糊核图7(f)与先验模糊核图7(e)相比，能量产生了较大的扩散；Krishnan的估计结果图7(g)，其能量分布呈现出不同方向的压缩和扩散；相比之下文中估计的模糊核图7(h)更加贴近先验模糊核的能量分布形式。由于模糊核的估计与先验模糊核差异较大，Zhang和Krishnan去模糊后的图像存在明显的振铃和伪影，如图7(b)和图7(c)所示。对比之下，文中结果图7(d)仍然能清晰地识别水面船只。

Figure 7.  Blur kernel estimation and deblurring of surface target image

为了定量测试文中方法估计模糊核的鲁棒性，统计了8个不同图像的估计模糊核与先验模糊核的相似度，相似度计算公式如下：

式中：$Cov({{\boldsymbol{H}}_1},{{\boldsymbol{H}}_2})$为模糊核${{\boldsymbol{H}}_{{1}}}$、${{\boldsymbol{H}}_2}$的协方差；$D({{\boldsymbol{H}}_{{1}}})$、$D({{\boldsymbol{H}}_2})$分别为${{\boldsymbol{H}}_{{1}}}$、${{\boldsymbol{H}}_2}$的方差；${\;\rho _{{{\boldsymbol{H}}_1}{{\boldsymbol{H}}_2}}}$取值0~1，值越大，代表相似度越高。统计结果如表1所示，文中对每个样本估计的模糊核与先验模糊核的相似度都高于对比方法；文中模糊核相似度均值最高且标准差Std最小，说明对于不同的图像，文中估计的模糊核能够稳定在先验模糊核周围。另外根据公式(21)~(22)计算了每个复原图像的清晰度指标，文中结果均优于对比方法。证明了文中方法估计的模糊核是稳定且有效的，具有较好的鲁棒性。

在算法的运行速度上，统计了4个不同大小的图像块计算文中半盲算法估计模糊核的时间，并且和Zhang、Krishnan的方法对比，皆在i7-7700 CPU @3.60 GHz/8 GB的平台上，利用MATLAB 2018a编程实现。结果如表2所示。

由结果可见，文中方法与Krishnan相比耗时较长，因为后者直接对图像做全局处理，没有提取图像结构信息；与Zhang方法相比，图像尺寸越大，文中速度优势越明显，因为在先验模糊核的约束下只需要迭代较少的次数，另外文中在计算过程中多次利用FFT算法加快了求解速度。

针对畸变引起大视场空间相机TDI成像模糊的问题，提出根据畸变数据和成像原理预估初始模糊核，将其作为先验信息，构建半盲复原模型估计最终模糊核的方法。将TDI成像过程中的畸变像移看作是多帧相对运动，通过拟合像移路径求解初始模糊核。提出一种多向权重异性的全变差模型提取图像边缘指导模糊核的细化。构建含有初始模糊核的正则项约束，使得模糊核的估计不容易受图像内容的影响。采用多尺度迭代方法细化模糊核，提高了模糊核估计精度，实验表明文中方法兼具有效性和鲁棒性，处理速度较快。对于提高大视场空间相机的TDI成像质量具有参考意义。