基于高斯曲率和加权图总变分正则化的遥感图像盲去模糊算法

Blind remote sensing image deblurring algorithm based on Gaussian curvature and reweighted graph total variation

摘要

针对遥感图像在采集过程中出现的运动模糊现象，本文从图像曲面的几何性质和图像像素的代数性质出发，设计了一种基于高斯曲率和加权图总变分正则化的遥感图像盲去模糊算法。首先，将加权图总变分先验与高斯曲率先验相结合以获得骨架图像，骨架图像保留图像的梯度以及锐利边缘信息，并去除中间潜在图像中的有害结构；然后，利用骨架图像估计模糊核，进而利用非盲去模糊算法获得清晰图像；最后，在8张不同场景下的模糊遥感图像上进行仿真验证。结果表明，相比于其他先进的图像盲去模糊算法，本文提出的去模糊算法复原效果的峰值信噪比平均值分别高于对比算法2.76、1.84、3.11、2.79、3.35、2.76 dB，结构相似性平均值分别高于对比算法0.0792、0.0604、0.0873、0.0801、0.0997、0.0906。本文算法复原的遥感图像具有清晰的边缘轮廓和局部细节，提升了遥感图像的清晰度。

关键词

计算数学 / 遥感图像 / 高斯曲率 / 加权图总变分 / 图像盲去模糊

Key words

computational mathematics / remote sensing image / Gaussian curvature / reweighted graph total variation / blind image deblurring

0 引言

遥感成像^［1］过程中成像系统与目标在曝光时间内相对运动造成的光轴偏移、像面抖动、目标物体移动以及地球自转等因素都会导致采集图像出现模糊现象。由于成像系统中的扩散函数是无法获知的，因此对遥感图进行盲去模糊是非常必要的^［2］。与自然图像相比，遥感图像具有覆盖范围广、地物信息多、目标尺寸差异大、边缘模糊等特点。因此，如何构建高精度的单幅遥感图像盲去模糊模型具有重要意义。

单幅图像盲去模糊问题是在模糊核未知的情况下，进行清晰图像恢复。单幅图像盲去模糊算法主要包括基于边缘选择的方法、基于深度学习的方法和基于图像先验的方法。基于边缘选择的方法利用启发性的滤波器或阈值截断等方式，显式地选择主要边缘结构估计模糊核以及具有大尺度边缘、细小边缘被抑制的中间潜在图像^［3，4］。基于边缘选择的方法依赖于滤波器和阈值的选取，容易造成图像过度锐化和噪声放大。基于深度学习的图像盲去模糊算法可分为非端到端盲去模糊算法和端到端盲去模糊算法。非端到端盲去模糊算法利用神经网络估计模糊核，然后应用非盲去模糊算法获得相应的清晰图像；端到端盲去模糊图像算法避开计算模糊核的过程，直接通过网络模型以端到端的方式由输入的模糊图像得到清晰图像，如基于卷积神经网络的去模糊模型^［5］、基于条件生成对抗网络的去模糊模型^［6］、基于循环神经网络的去模糊模型^［7］等。基于深度学习的盲去模糊算法需要大量数据集进行训练并依赖于设备强大的计算力，成本较高。基于图像先验的盲去模糊算法通过对图像添加先验克服图像去模糊的病态问题，从而复原出清晰图像，如图像范数先验^［8］、梯度图像范数先验^［9，10］、暗通道先验^{［11，12］}、极通道先验^［13］、非线性通道先验^［14］、局部最大梯度先验^［15］、低秩先验^{［16，17］}、表面感知先验^［4，18］、块稀疏表示先验^{［19，20］}、图信号先验^{［21，22］}、局部像素稀疏先验^［23］等。

近年来，国内外很多学者针对遥感图像纹理信息丰富的特性，设计了各种遥感图像去模糊模型。Zhang等^［24］提出了基于局部二进制稀疏先验的遥感图像盲去模糊算法，该算法可以较完整地保留图像的纹理特征，复原后的遥感图像不会因为产生模糊而发生很大变化。Li等^［25］设计了基于参考的多级特征融合去模糊网络，从多级参考特征中传输高质量的纹理辅助图像去模糊。Lim等^［26］提出了基于图像范数和纹理感知的遥感图像盲去模糊方法，利用基于图像范数去模糊模型得到复原图像的大部分信息，利用纹理感知方法保留图像的锐利边缘和纹理细节。王琪瑶等^［27］提出了基于局部最大和最小强度稀疏先验的遥感图像盲去模糊算法。刘晨辉等^［28］利用拍摄过程中获得的不同积分时间遥感图像序列，从短积分时间图像中提取出轮廓边缘信息估计相邻长积分时间图像模糊核。

近年来，国内外很多学者将图像曲面的曲率滤波成功应用于图像去模糊领域。文献［4］提出了基于几何特性图像曲面感知策略的盲去模糊算法，该方法本质上是求解图像曲面的中值曲率最小化问题。文献［18］提出了基于图像曲面的高斯曲率范数和梯度图像范数正则化盲去模糊模型，验证了在经典基于梯度图像范数正则化去模糊模型中，引入高斯曲率正则化可以更好地保留清晰图像的边缘信息，并去除中间潜在图像中的有害结构信息。

基于图信号先验的图像去模糊方法将像素建模成具有反映像素间相似性的加权边节点，将图像理解为图上的信号，在图像去模糊领域取得了较好的成绩。文献［22］提出了基于加权图总变分先验的去模糊模型，即利用双模态边缘权重分布的加权图总变分先验，从模糊图像的骨架信息中恢复模糊核。

本文在文献［18］和文献［22］工作的基础上，设计了基于高斯曲率和加权图总变分正则化的遥感图像盲去模糊算法。高斯曲率先验可以更好地指导模糊核的估计，加权图总变分先验可以提高图像的锐度和细节。将这两种先验相结合，可以在保留边缘和细节信息的同时提高图像的锐度和清晰度，有效解决单一先验方法存在的问题。

1 理论部分

1.1 高斯曲率正则化

令 \mathit{z} = \mathit{z} ( x , y ) , ( x , y ) \in \Omega为图像 \mathit{z} ( x , y )所对应的图像曲面。令 \mathit{\psi } ( x , y , \mathit{z} ) = \mathit{I} ( x , y ) - \mathit{z}，则： \mathit{\psi } ( x , y , \mathit{z} ) = \mathrm{0}为曲面 \mathit{z} = \mathit{z} ( x , y )的隐函数表达式。图像曲面 \mathit{z} = \mathit{z} ( x , y )的曲率 K_{\mathrm{G}\mathrm{C}}为：

K_{\mathrm{G}\mathrm{C}} = \frac{\nabla \mathit{\psi }{H}^{\mathrm{*}}(\mathit{\psi })\nabla {\mathit{\psi }}^{\mathrm{T}}}{{‖\nabla \mathit{\psi }‖}_{\mathrm{2}}^{\mathrm{4}}} （1）

式中： {‖\nabla \mathit{\psi }‖}_{\mathrm{2}}为梯度向量 \nabla \mathit{\psi }的2-范数； H^{\mathrm{*}} ( \mathit{\psi } )为函数 \mathit{\psi } ( x , y , \mathit{z} )的Hessian矩阵的伴随矩阵。

文献［18］中，曲面 \mathit{z} = \mathit{z} ( x , y )的高斯曲率 K_{\mathrm{G}\mathrm{C}}^{\text{'}}为：

K_{\mathrm{G}\mathrm{C}}^{\text{'}} = \frac{{\mathit{z}}_{xx}{\mathit{z}}_{yy}-{\mathit{z}}_{xy}{\mathit{z}}_{yx}}{({\mathit{z}}_x^{\mathrm{2}}+{\mathit{z}}_y^{\mathrm{2}}{+\mathrm{1})}^{\mathrm{2}}} （2）

式中： {\mathit{z}}_{xx}、 {\mathit{z}}_{yy}、 {\mathit{z}}_{xy}、 {\mathit{z}}_{yx}、 {\mathit{z}}_x、 {\mathit{z}}_y分别为 \mathit{z}的二阶偏导数和一阶偏导数。

文献［18］基于式（2），引入结合梯度图像 L_{\mathrm{0}}范数和图像高斯曲率 L_{\mathrm{1}}范数的高斯曲率正则化因子 \phi，如下所示：

\phi ( \mathit{z} ) = \delta \left({‖{\mathit{z}}_x‖}_{\mathrm{0}}+{‖{\mathit{z}}_y‖}_{\mathrm{0}}\right) + {\iint }_{\Omega }^{} \left|K_{\mathrm{G}\mathrm{C}}\right|\mathrm{d}x\mathrm{d}y （3）

式中： \delta为正则化因子。对潜在清晰图像添加先验，这样有利于进行模糊核估计以及恢复潜在清晰图像的结构。

1.2 加权图总变分正则化

文献［22］引入加权图总变分正则化恢复保留目标强梯度且细节平滑的骨架图像。定义图像 \mathit{z}在邻近像素点 ( i , j )、 ( k , l )的局部邻接权重 w_{i,j}^{(k,l)}为：

w_{i,j}^{(k,l)} = \mathrm{e} \mathrm{x} \mathrm{p} \left(-{\left|z_{i,j}-z_{k,l}\right|}^{\mathrm{2}}/{\sigma }^{\mathrm{2}}\right) （4）

式中： z_{i,j}和 z_{k,l}分别为二维图像 \mathit{z}在像素点 ( i , j )和 ( k , l )处的像素值； \sigma为参数。

由式（4）可知， w_{i,j}^{(k,l)}值越大，像素点 ( i , j )与 ( k , l )处的像素值越接近。为简便起见，设图像 \mathit{z}包含 m \times n个像素点，设 V {}_{i,j}{}^{}为像素点 ( i , j )的局部邻域点集合，则图像 \mathit{z}的加权图总变分正则化因子 {‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}}为：

{‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}} = {\displaystyle \sum _{i=\mathrm{1}}^m} {\displaystyle \sum _{j=\mathrm{1}}^n}\sum _{(k,l)\in V_{i,j}}{w}_{i,j}^{(k,l)}\left|z_{i,j}-z_{k,l}\right| （5）

V_{i,j} = \{ ( i - \mathrm{1} , j ) , ( i + \mathrm{1} , j ) , ( i , j - \mathrm{1} ) , ( i , j + \mathrm{1} ) \} （6）

令 d = \left|z_{i,j}-z_{k,l}\right|，则由式（4）可知，正则化曲面 y为：

y = w_{i,j}^{(k,l)} \left|z_{i,j}-z_{k,l}\right|

当 d = \sigma / \sqrt[]{\mathrm{2}}时， y取得极大值；当 d = \mathrm{0}或 d = + \mathrm{\infty }时， y取得极小值。当 时，最小化式（5）意味着缩小 d的取值；当 d > \sigma / \sqrt[]{\mathrm{2}}时，最小化式（5）意味着放大 d的取值。因此，加权图总变分正则化因子可以有效促进清晰图像的双模边缘权值分布。

利用局部邻接权重 w_{i,j}^{(k,l)}定义图像 \mathit{z}在邻近像素点 ( i , j )、 ( k , l )点处局部Laplacian权重 {\gamma }_{i,j}^{(k,l)}为：

{\gamma }_{i,j}^{(k,l)} = \frac{w_{i,j}^{(k,l)}}{\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|z_{i,j}-z_{k,l}\right|,\epsilon \right\}} （7）

式中： \epsilon为非常小的正实数，目的是为了避免分母为零，本文 \epsilon取值为0.01。

文献［20］指出，正则化因子（5）可以用图Laplacian正则化因子

{{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}\sum _{(k,l)\in V_{i,j}}{\gamma }_{i,j}^{(k,l)}(z_{i,j}-z_{k,l})}^{\mathrm{2}} （8）

近似代替。

令 \mathrm{v} \mathrm{e} \mathrm{c} ( \mathit{z} )和 \mathrm{v} \mathrm{e} \mathrm{c} ( {\mathit{z}}^{\mathrm{T}} )分别为将二维矩阵 \mathit{z}和 {\mathit{z}}^{\mathrm{T}}按列展开得到的 n · m维列向量。对于 \mathrm{1} \le i \le m，定义 n \times n维三对角线矩阵 {\mathit{W}}_x^{(i)}，其主对角线元素为零，主对角元素之上和之下的对角线均由向量 ( {\gamma }_{i,\mathrm{1}}^{(i,\mathrm{2})} , {\gamma }_{i,\mathrm{2}}^{(i,\mathrm{3})} , \cdots , {\gamma }_{i,n-\mathrm{1}}^{(i,n)} {)}^{\mathrm{T}}生成；对于 \mathrm{1} \le j \le n，定义 m \times m维三对角线矩阵 {\mathit{W}}_y^{(j)}，其主对角线元素为零，主对角元素之上和之下的对角线均由向量 ( {\gamma }_{\mathrm{1},j}^{(\mathrm{2},j)} , {\gamma }_{\mathrm{2},j}^{(\mathrm{3},j)} , \cdots , {\gamma }_{m-\mathrm{1},j}^{(m,j)} {)}^{\mathrm{T}}生成。引入如下 ( m · n ) \times ( m · n )维矩阵 {\mathit{W}}_x和 {\mathit{W}}_y：

{\mathit{W}}_x = \left(\begin{array}{cccc}{W}_x^{(\mathrm{1})}& & & \\ & W_x^{(\mathrm{2})}& & \\ & & \ddots & \\ & & & W_x^{(m)}\end{array}\right)

{\mathit{W}}_y = \left(\begin{array}{cccc}{W}_y^{(\mathrm{1})}& & & \\ & W_y^{(\mathrm{2})}& & \\ & & \ddots & \\ & & & W_y^{(n)}\end{array}\right)

令 \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} ( · )表示生成对角线元素表示输入向量的方形对角线矩阵的算子，令 \mathrm{o} \mathrm{n} \mathrm{e} \mathrm{s} ( m · n , \mathrm{1} )表示分量均为1的 m · n维向量。利用 {\mathit{W}}_x和 {\mathit{W}}_y，定义 ( m · n ) \times ( m · n )维矩阵 {\mathit{L}}_x和 {\mathit{L}}_y分别为：

{\mathit{L}}_x = \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} ( {\mathit{W}}_x \mathrm{o} \mathrm{n} \mathrm{e} \mathrm{s} ( m · n , \mathrm{1} ) ) - {\mathit{W}}_x

{\mathit{L}}_y = \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} ( {\mathit{W}}_y \mathrm{o} \mathrm{n} \mathrm{e} \mathrm{s} ( m · n , \mathrm{1} ) ) - {\mathit{W}}_y

则当局部邻域点集合为式（6）时，图Laplacian正则化因子（式（8））可改写为：

\mathrm{v} \mathrm{e} \mathrm{c} {(\mathit{z})}^{\mathrm{T}} {\mathit{L}}_y \mathrm{v} \mathrm{e} \mathrm{c} ( \mathit{z} ) + \mathrm{v} \mathrm{e} \mathrm{c} ( {\mathit{z}}^{\mathrm{T}} {)}^{\mathrm{T}} {\mathit{L}}_x \mathrm{v} \mathrm{e} \mathrm{c} ( {\mathit{z}}^{\mathrm{T}} ) （9）

定义 m \times n维矩阵 {\mathit{\Gamma }}_x、 {\mathit{\Gamma }}_y分别为：

{\mathit{\Gamma }}_x = \left(\begin{array}{ccccc}{\gamma }_{\mathrm{1,1}}^{(\mathrm{1,2})}& {\gamma }_{\mathrm{1,2}}^{(\mathrm{1,3})}& \cdots & {\gamma }_{\mathrm{1},n-\mathrm{1}}^{(\mathrm{1},n)}& \mathrm{0}\\ {\gamma }_{\mathrm{2,1}}^{(\mathrm{2,2})}& {\gamma }_{\mathrm{2,2}}^{(\mathrm{2,3})}& \cdots & {\gamma }_{\mathrm{2},n-\mathrm{1}}^{(\mathrm{2},n)}& \mathrm{0}\\ ⋮& ⋮& & ⋮& ⋮\\ {\gamma }_{m,\mathrm{1}}^{(m,\mathrm{2})}& {\gamma }_{m,\mathrm{2}}^{(m,\mathrm{3})}& \cdots & {\gamma }_{m,n-\mathrm{1}}^{(m,n)}& \mathrm{0}\end{array}\right)

{\mathit{\Gamma }}_y = \left(\begin{array}{cccc}{\gamma }_{\mathrm{1,1}}^{(\mathrm{2,1})}& {\gamma }_{\mathrm{1,2}}^{(\mathrm{2,2})}& \cdots & {\gamma }_{\mathrm{1},n}^{(\mathrm{2},n)}\\ {\gamma }_{\mathrm{2,1}}^{(\mathrm{3,1})}& {\gamma }_{\mathrm{2,2}}^{(\mathrm{3,2})}& \cdots & {\gamma }_{\mathrm{2},n}^{(\mathrm{3},n)}\\ ⋮& ⋮& & ⋮\\ {\gamma }_{m-\mathrm{1,1}}^{(m,\mathrm{1})}& {\gamma }_{m-\mathrm{1,2}}^{(m,\mathrm{2})}& \cdots & {\gamma }_{m-\mathrm{1,1}}^{(m,\mathrm{1})}\\ \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\end{array}\right)

设 {\mathit{P}}_{{\pi }_x}和 {\mathit{P}}_{{\pi }_y}分别为对应置换 {\pi }_x = ( \mathrm{2,3} , \cdots , n - \mathrm{1,1} )和 {\pi }_y = ( m , \mathrm{1,2} , \cdots , m - \mathrm{1} )的置换矩阵，令 \nabla = ( {\nabla }_x , {\nabla }_y {)}^{\mathrm{T}}为梯度算子，经计算可知，对于固定的 {\mathit{L}}_x和 {\mathit{L}}_y，式（9）关于矩阵变量矩阵 \mathit{z}的偏导数为：

\begin{array}{c}\mathrm{2}(({\mathit{\Gamma }}_x{\mathit{P}}_{{\pi }_x})\otimes {\nabla }_x\mathit{z}-{\mathit{\Gamma }}_x\otimes {\nabla }_x\mathit{z}+\\ ({\mathit{\Gamma }}_y{\mathit{P}}_{{\pi }_y})\otimes {\nabla }_y\mathit{z}-{\mathit{\Gamma }}_y\otimes {\nabla }_y\mathit{z})\end{array} （10）

式中： \otimes为卷积算子。

2 本文算法

本文设计了结合高斯曲率正则化和加权图总变分正则化的遥感图像盲去模糊模型，即：

\begin{array}{c}\underset{\mathit{z},\mathit{k}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{b}-\mathit{k}\otimes \mathit{z}‖}_{\mathrm{2}}^{\mathrm{2}}+\beta {‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}}+\\ \mu {‖\mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda \phi (\mathit{z})\end{array} （11）

式中： \mathit{b}为模糊图像； \mathit{k}为模糊核； \beta \mathrm{、} \mu \mathrm{、} \lambda均为参数。

本文将求解式（11）分为估计骨架图像恢复和模糊核恢复两部分，即求解以下两个子问题：

\underset{\mathit{z}}{\mathrm{m}\mathrm{i}\mathrm{n}} \frac{\mathrm{1}}{\mathrm{2}} {‖\mathit{b}-\mathit{k}\otimes \mathit{z}‖}_{\mathrm{2}}^{\mathrm{2}} + \beta {‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}} + \lambda \phi ( \mathit{z} ) （12）

\underset{\mathit{k}}{\mathrm{m}\mathrm{i}\mathrm{n}} \frac{\mathrm{1}}{\mathrm{2}} {‖\mathit{b}-\mathit{k}\otimes \mathit{z}‖}_{\mathrm{2}}^{\mathrm{2}} + \mu {‖\mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}} （13）

2.1 骨架图像恢复

求解式（12）是非凸问题，故利用半二次分裂算法，引入与 \mathit{b}维数相同的辅助变量 \mathit{u}和参数 \alpha，将式（10）改写为：

\begin{array}{l}\underset{\mathit{z},\mathit{u}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{b}-\mathit{z}\otimes \mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}}+\beta {‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}}+\\ \alpha {‖\mathit{z}-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda \phi (\mathit{u})\end{array} （14）

骨架图像恢复可分为以下两个子问题：

\underset{\mathit{z}}{\mathrm{m}\mathrm{i}\mathrm{n}} \frac{\mathrm{1}}{\mathrm{2}} {‖\mathit{b}-\mathit{z}\otimes \mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}} + \beta {‖\mathit{z}‖}_{\mathrm{R}\mathrm{G}\mathrm{T}\mathrm{V}} + \alpha {‖\mathit{z}-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}} （15）

\underset{\mathit{u}}{\mathrm{m}\mathrm{i}\mathrm{n}} \alpha {‖\mathit{z}-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}} + \lambda \phi ( \mathit{u} ) （16）

对于固定 \mathit{u}和 \mathit{k}，模型（15）是非可微问题，故本文利用图Laplacian正则化方法求解优化问题（15）。利用式（9）将模型（15）转化为如下模型求解：

\begin{array}{l}\underset{\mathit{z}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{b}-\mathit{z}\otimes \mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}}+\beta (\mathrm{v}\mathrm{e}\mathrm{c}{(\mathit{z})}^{\mathrm{T}}{\mathit{L}}_y\mathrm{v}\mathrm{e}\mathrm{c}(\mathit{z})+\\ \mathrm{v}\mathrm{e}\mathrm{c}({\mathit{z}}^{\mathrm{T}}{)}^{\mathrm{T}}{\mathit{L}}_x\mathrm{v}\mathrm{e}\mathrm{c}({\mathit{z}}^{\mathrm{T}}))+\alpha {‖\mathit{z}-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}}\end{array} （17）

由式（10）可知，模型（17）的闭合解为线性系统的解。利用共轭梯度法求解如下线性系统：

\begin{array}{l}\mathrm{2}\beta \left(({\mathit{\Gamma }}_x{\mathit{P}}_{{\pi }_x})\otimes {\nabla }_x\mathit{z}-{\mathit{\Gamma }}_x\otimes {\nabla }_x\mathit{z}+({\mathit{\Gamma }}_y{\mathit{P}}_{{\pi }_y})\otimes \right.\\ \left.{\nabla }_y\mathit{z}-{\mathit{\Gamma }}_y\otimes {\nabla }_y\mathit{z}\right)+{\mathit{k}}^{\mathrm{T}}\otimes (\mathit{k}\otimes \mathit{z})+\\ \alpha \mathit{z}={\mathit{k}}^{\mathrm{T}}\otimes \mathit{b}+\alpha \mathit{u}\end{array} （18）

可得到估计的骨架图像 \mathit{z}。固定 \mathit{z}，本文利用文献［18］提出基于高斯曲率的方法求解问题（14）。图像 \mathit{u}可以理解为骨架图像 \mathit{z}经过融合梯度图像 L_{\mathrm{1}}范数和图像高斯曲率 L_{\mathrm{0}}范数的高斯曲率正则化因子获取修正后骨架图像。首先，将二维图像的像素点分为4个集合：黑色矩形 B_{\mathrm{R}}、白色矩形 W_{\mathrm{R}}、黑色圆形 B_{\mathrm{C}}和白色圆形 W_{\mathrm{C}}，如图1所示。利用算法1和算法2求解模型（16），以获得修正后的骨架图像。

算法1 最小投影算子算法

输入：骨架图像 \mathit{z} ( x , y )

计算：

d_{\mathrm{1}} = ( z ( x - \mathrm{1} , y ) + z ( x + \mathrm{1} , y ) ) / \mathrm{2} - z ( x , y )

d_{\mathrm{2}} = ( z ( x , y - \mathrm{1} ) + z ( x , y + \mathrm{1} ) ) / \mathrm{2} - z ( x , y )

d_{\mathrm{3}} = ( z ( x - \mathrm{1} , y - \mathrm{1} ) + z ( x + \mathrm{1} , y + \mathrm{1} ) ) / \mathrm{2} - z ( x , y )

d_{\mathrm{4}} = ( z ( x - \mathrm{1} , y + \mathrm{1} ) + z ( x + \mathrm{1} , y - \mathrm{1} ) ) / \mathrm{2} - z ( x , y )

d_{\mathrm{5}} = z ( x - \mathrm{1} , y ) + z ( x , y - \mathrm{1} ) - z ( x - \mathrm{1} , y - \mathrm{1} ) - z ( x , y )

d_{\mathrm{6}} = z ( x - \mathrm{1} , y ) + z ( x , y + \mathrm{1} ) - z ( x - \mathrm{1} , y + \mathrm{1} ) - z ( x , y )

d_{\mathrm{7}} = z ( x , y - \mathrm{1} ) + z ( x + \mathrm{1} , y ) - z ( x + \mathrm{1} , y - \mathrm{1} ) - z ( x , y )

d_{\mathrm{8}} = z ( x , y + \mathrm{1} ) + z ( x + \mathrm{1} , y ) - z ( x + \mathrm{1} , y + \mathrm{1} ) - z ( x , y )

确定 d_m，满足 \left|d_m\right| = \mathrm{m} \mathrm{i} \mathrm{n} \{ \left|d_i\right| , i = \mathrm{1,2} , \dots , \mathrm{8} \}

令 z ( x , y ) \leftarrow z ( x , y ) + d_m

输出：更新后骨架图像 \mathit{z} ( x , y )。

算法2 模型（13）的求解算法

输入：骨架图像 \mathit{z} ( x , y )，最大迭代次数 t_{\mathrm{m}\mathrm{a}\mathrm{x}}，参数 \alpha和 \lambda。

步骤1 将图像 \mathit{z} ( x , y )的像素点分类： B_{\mathrm{C}} \mathrm{、} B_{\mathrm{R}} \mathrm{、} W_{\mathrm{C}} \mathrm{、} W_{\mathrm{R}}。

步骤2 令 t \leftarrow \mathrm{1} {\mathit{u}}_{\mathrm{1}} ( x , y ) \leftarrow \mathit{z} ( x , y )。

步骤3 令 P为最小投影算子

while t \le t_{\mathrm{m}\mathrm{a}\mathrm{x}}do

\forall ( x , y ) \in B_{\mathrm{C}} , u_{n+\mathrm{1}/\mathrm{4}} ( x , y ) = P ( u_n ( x , y ) )

\forall ( x , y ) \in B_{\mathrm{R}} , u_{n+\mathrm{2}/\mathrm{4}} ( x , y ) = P ( u_{n+\mathrm{1}/\mathrm{4}} ( x , y ) )

\forall ( x , y ) \in W_{\mathrm{C}} , u_{n+\mathrm{3}/\mathrm{4}} ( x , y ) = P ( u_{n+\mathrm{2}/\mathrm{4}} ( x , y ) )

\forall ( x , y ) \in W_{\mathrm{R}} , u_{n+\mathrm{1}} ( x , y ) = P ( u_{n+\mathrm{3}/\mathrm{4}} ( x , y ) )

若

令 t \leftarrow \mathit{u}，算法中止。

否则令 t \leftarrow t + \mathrm{1}

end

输出：修正后的骨架图像 \mathit{u} ( x , y )。

2.2 模糊核恢复

为了抑制伪影现象，利用文献［11］的方法将模型（9）转换到梯度图像求解，由于修正后的骨架图像 \mathit{u}能更好地保留图像的边缘信息，将模型（9）中的骨架图像 \mathit{z}替换为修正后的骨架图像 \mathit{u}。综上，模型（13）可转化为以下模型：

\underset{\mathit{k}}{\mathrm{m}\mathrm{i}\mathrm{n}} \frac{\mathrm{1}}{\mathrm{2}} {‖\nabla \mathit{b}-\nabla \mathit{u}\otimes \mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}} + \mu {‖\mathit{k}‖}_{\mathrm{2}}^{\mathrm{2}} （19）

优化问题（19）的目标函数是二次函数，可以利用反卷积算法对模型（19）求解，本文利用快速傅里叶变换加速反卷积求解过程，即：

\mathit{k} = F^{-\mathrm{1}} \left(\frac{\overline{F(\nabla \mathit{u})}F(\nabla \mathit{b})}{\overline{F(\nabla \mathit{u})}F(\nabla \mathit{u})+\mu }\right) （20）

式中： F ( · )为快速傅里叶变换； F^{-\mathrm{1}} ( · )为快速傅里叶逆变换； \overline{F(·)}为快速傅里叶变换的共轭算子。

利用式（20）得到模糊核 k后，将 k中非负元素置为零，并将 k正规化使得其所有分量求和为1。

2.3 遥感图像盲去模糊算法

本文利Fergus等^［29］提出的基于图像金字塔的多尺度“由粗到细”策略，改善模糊核的估计精度。基于“由粗到细”策略的盲去模糊算法流程图如图2所示。首先对模糊图像进行降采样操作，在分辨率最低的灰度图像上利用本文算法估计骨架图像和模糊核；将估计的模糊核进行上采样并用于下一尺度的模糊遥感图像上，重复进行交替迭代优化，直到将本文提出的去模糊算法应用于最精细尺度的模糊图像上以获得最终估计的模糊核；最后，利用估计的模糊核和用非盲反卷积方法获得清晰图像。图3为在“由粗到细”盲去模糊过程中所得到的模糊核和骨架图像。本文算法流程如算法3所示。

算法3 遥感图像盲去模糊算法

输入：模糊遥感图像 \mathit{b}，最大迭代次数 N和 t_{\mathrm{m}\mathrm{a}\mathrm{x}}，参数 \alpha， \beta， \lambda。

步骤1 建立图像 \mathit{b}由低到高的金字塔图像 b_{\mathrm{1}} , b_{\mathrm{2}} , \cdots , b_J。

步骤2 for j = \mathrm{1,2} , \cdots , J

初始化 \mathit{u} \leftarrow b_j， i \leftarrow \mathrm{1}，参数 \alpha、 \beta、 \lambda取输入值

while i \le N do

利用共轭梯度法求当 \mathit{b} = b_j， \mathit{k} = k_j时式（17）的解 \mathit{z}

求解模型（15）得到修正后骨架图像 \mathit{u} ( x , y )

利用式（18）估计模糊图像 b_j对应的模糊核 k_j

令 i \leftarrow i + \mathrm{1} , \alpha \leftarrow \alpha / \mathrm{1.2} , \beta \leftarrow \beta / \mathrm{1.1} , \lambda \leftarrow \mathrm{0.9} \lambda

end

对 k_j上采样得到模糊核 k_{j+\mathrm{1}}

end

输出：估计的模糊核。

3 实验结果及分析

本文在处理器为Inter（R） Core（TM） CPU i5-8250U@1.60 GHz运行内存8 GB的64位的Windows10操作系统上，利用Matlab2020进行实验。选取来自AID数据集中空间分辨率为600×600的海滩、沙漠、学校、度假村、机场、山地、广场和高架桥遥感图像，如图4所示。对于这8张清晰遥感图像，利用Levin数据集^［30］中提供的8个模糊核生成8张模糊遥感图像，进行仿真实验。

本文采用峰值信噪比（Peak signal to noise ratio，PSNR）和结构相似性（Structural similarity，SSIM）验证本文方法的有效性。

峰值信噪比的表达式为：

\mathrm{P} \mathrm{S} \mathrm{N} \mathrm{R} = \mathrm{10} · \mathrm{l} \mathrm{g} \frac{L^{\mathrm{2}}}{\mathrm{M}\mathrm{S}\mathrm{E}} （21）

\mathrm{M} \mathrm{S} \mathrm{E} = \frac{\mathrm{1}}{N\times M} {\displaystyle \sum _{i=\mathrm{1}}^N} {\displaystyle \sum _{j=\mathrm{1}}^M}{s}_{ij}-{\widehat{s}}_{ij} （22）

式中： M \mathrm{、} N分别为 x \mathrm{、} y方向像素点的个数； s_{ij}和 {\widehat{s}}_{ij}分别为原始图像和去模糊后的重建图像在点 ( i , j )上的取值； L为图像灰度值的范围；MSE为均方误差。

结构相似性的公式为：

\mathrm{S} \mathrm{S} \mathrm{I} \mathrm{M} ( x , y ) = l ( x , y ) · c ( x , y ) · s ( x , y ) （23）

\left\{\begin{array}{l}l(x,y)=\frac{\mathrm{2}{\mu }_x{\mu }_y+C_{\mathrm{1}}}{{\mu }_x^{\mathrm{2}}+{\mu }_y^{\mathrm{2}}+C_{\mathrm{1}}}\\ c(x,y)=\frac{\mathrm{2}{\sigma }_x{\sigma }_y+C_{\mathrm{2}}}{{\sigma }_x^{\mathrm{2}}+{\sigma }_y^{\mathrm{2}}+C_{\mathrm{2}}}\\ s(x,y)=\frac{{\sigma }_{xy}+C_{\mathrm{3}}}{{\sigma }_x{\sigma }_y+C_{\mathrm{3}}}\end{array}\right. （24）

式中： {\mu }_x和 {\mu }_y分别为图像 x和 y的均值； {\sigma }_x和 {\sigma }_y分别为图像 x和 y的方差； {\sigma }_{xy}为图像 x和 y的协方差； C_{\mathrm{1}} \mathrm{、} C_{\mathrm{2}} \mathrm{、} C_{\mathrm{3}}均为常数。

本文选取文献［11］、文献［22］、文献［15］、文献［14］、文献［23］和文献［18］中的盲去模糊算法进行对比实验。为了公平比较，本文统一采用基于Laplacian先验的快速反卷积算法^［31］进行非盲去模糊，并且本文所有算法的非盲去模糊部分使用相同参数。进一步，对比实验中的盲去模糊算法是在原始代码的默认参数下进行的。

在数值实验中，本文算法的参数 \sigma、 \beta、 \mu的选取与文献［22］相同，即 \sigma = \mathrm{0.1}、 \beta = \mathrm{0.01}、 \mu = \mathrm{0.05}。另有参数 \alpha \mathrm{和} \lambda，其中 \alpha用于控制骨架图像恢复中的辅助变量 \mathit{u}与骨架图像 \mathit{z}之间的差异，不是直接关注于去模糊效果，而 \lambda用于控制高斯曲率正则化因子对骨架图像恢复的影响，直接关注于去模糊效果。引入的参数 \lambda和 \alpha设置如下： \lambda = \mathrm{0.5}、 \alpha \in \{ \mathrm{0.01,0.015,0.018} \}。当模糊核尺寸大于25时， \alpha = \mathrm{0.018}，当模糊核尺寸小于15时， \alpha = \mathrm{0.01}，其他情况 \alpha = \mathrm{0.015}。

图5为参数 \lambda在不同取值下去模糊后图像的PSNR和SSIM值，说明了参数 \lambda选取的有效性。图6为遥感图像盲去模糊算法运行时间。由图6可知，本文算法的运行时间略高于文献［22］和文献［23］算法的运行时间，低于其余4种对比算法的运行时间。

表1为模糊遥感图像去模糊后的PSNR值和SSIM值。由表1可以看出：除了机场图像，本文算法复原图像的SSIM值略低于文献［14］的结果，其他情况下本文算法复原图像的PSNR和SSIM结果均高于其他对比实验的结果。对于纹理信息丰富的山地图像，本文算法的PSNR值比其他6种算法依次高4.82、4.77、5.19、4.78、4.74、4.28 dB，SSIM值依次高0.2160、0.2397、0.2222、0.2131、0.2145、0.2097。由表1可见，本文算法对纹理特征丰富的遥感图像，去模糊效果提升更明显。

图7为模糊遥感图像去模糊效果对比图。由图7可见，经本文算法处理后的图像质量有明显提升，计算出的模糊核接近真实模糊核，丢失的纹理细节也得到了一定程度的恢复。

4 结束语

针对遥感图像在成像过程中出现的运动模糊现象，本文提出了基于高斯曲率和加权图总变分正则化的遥感图像盲去模糊算法。本文将高斯曲率正则化因子嵌入基于加权图总变分正则化的遥感图像盲去模糊模型中，以获得保留图像梯度信息和边缘信息，并丢弃潜在清晰图像中有害结构信息的骨架图像，利用骨架图像计算模糊核，从而得到清晰图像。与其他算法比，本文所提出的算法无论对于纹理信息丰富的遥感图像，还是对于具有较高的非局部自相似性的遥感图像，都取得了优于对比算法的结果。本文算法的不足之处在于对于尺寸比较小的模糊核，算法的优越性不明显。如何提高本算法对小尺寸模糊核造成的遥感图像模糊的处理能力是今后工作的重点。

参考文献

[1] 王李祺, 候宇超, 高翔, 基于深度学习与狮群SVM算法的遥感场景分类[J]. 吉林大学学报: 理学版, 2023, 61(4): 863-874.

[2] Wang Li-qi, Hou Yu-chao, Gao Xiang, et al. Remote sensing scene classfication based on deep learning and lion swarm SVM algorithm[J]. Journal of Jilin University(Science Edition), 2023, 61(4): 863-874.

[3] 张玉波, 王建阳, 韩爽, 一种非对称的轻量级图像盲与模糊网络[J]. 吉林大学学报: 理学版, 2023, 61(2): 362-370.

[4] Zhang Yu-bo, Wang Jian-yang, Han Shuang, et al. An asymmetric lightweight image blind deblurring network[J]. Journal of Jilin University(Science Edition), 2023, 61(2): 362-370.

[5] Lai W S, Ding J J, Lin Y Y, et al. Blur kernel estimation using normalized color-line priors[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Boston, USA, 2015: 64-72.

[6] Liu J, Yan M, Zeng T Y. Surface aware blind image deblurring[J]. IEEE transactions on Pattern Analysis and Machine Intelligence, 2021, 43(3): 1041-1055.

[7] Nah S, Kim T H, Lee K M. Deep multi-scale convolutional neural network for dynamic scene deblurring[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Hawaii, USA, 2017: 257-265.

[8] Kupyn O, Budzan V, Mykhailych M, et al. DeblurGAN: blind motion deblurring using conditional adversarial neworks[C]//IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Salt Lake City, USA, 2018: 8183-8192.

[9] Zhang J W, Pan J S, Ren J, et al. Dynamic scene deblurring using spatially variant recurrent neural networks[C]//IEEE/CVF Conference on Computer Vision and Pattern Recognition(CVPR), Salt Lake City, USA, 2018: 2521-2529.

[10] Krishnan D, Tay T, Fergus R. Blind deconvolution using a normalized sparsity measure[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Colorado Springs, USA, 2011: 233-240.

[11] Xu L, Zheng S C, Jia J Y. Unnatural L0 sparse representation for natural image deblurring[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Portland, USA, 2013: 1107-1114.

[12] Pan J, Su Z. Fast regularized kernel estimation for robust motion deblurring[J]. IEEE Signal Processing Letters, 2013, 20(9): 841-844.

[13] Pan J S, Sun D Q, Pfister H, et al. Deblurring images via dark channel prior[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017, 40(10): 2315-2328.

[14] 岳有军, 云赛, 王红君, 基于暗通道与低秩先验的运动模糊图像盲复原[J]. 电光与控制, 2019, 26(11): 95-98, 110.

[15] Yue You-jun, Yun Sai, Wang Hong-jun, et al. Blind restoration of motion-blurred image based on dark channel and low-rank prior[J]. Electronics Optics ＆ Control, 2019, 26(11): 95-98, 110.

[16] Yan Y Y, Ren W Q, Gou Y F, et al. Image deblurring via extreme channels prior[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Honolulu, USA, 2017: 6978-6986.

[17] Ge X, Tan J, Zhang L. Blind image deblurring using a non-linear channel prior based on dark and bright channels[J]. IEEE Transactions on Image Processing, 2021, 30: 6970-6984.

[18] Chen L, Fang F, Wang T, et al. Blind image deblurring with local maximum gradient prior[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition(CVPR), Long Beach, USA, 2019: 1742-1750.

[19] Ren Wen-qi, Cao Xiao-chun, Pan Jin-shan, et al. Image deblurring via enhanced low rank prior[J]. IEEE Transactions on Image Processing, 2016, 25(7): 3426-3437.

[20] Dong J X, Pan J S, Su Z X. Blur kernel estimation via salient edges and low rank prior for blind image deblurring[J]. Signal Processing: Image Communication, 2017, 58: 134-145.

[21] Ge X Y, Tan J Q, Zhang L, et al. Blind image deblurring with Gaussian curvature of the image surface[J]. Signal Processing: Image Communication, 2022, 100: No.116531.

[22] Tang Y, Xue Y, Chen Y, et al. Blind deblurring with sparse representation via external patch priors[J]. Digital Signal Processing, 2018, 78: 322-331.

[23] Yu J, Chang Z, Xiao C. Blur kernel estimation using sparse representation and cross-scale self-similarity[J]. Multimedia Tools and Applications, 2019, 78(13): 18549-18570.

[24] Berger P, Hannak G, Matz G. Graph signal recovery via primaldual algorithms for total variation minimization[J]. IEEE Journal of Selected Topics in Signal Processing, 2017, 11(6): 842-855.

[25] Bai Y, Cheung G, Liu X, et al. Graph based blind image deblurring from a single photograph[J]. IEEE Transactions on Image Processing, 2019, 28(3): 1404-1418.

[26] Wen F, Ying R D, Liu Y P, et al. A simple local minimal intensity prior and an improved algorithm for blind image deblurring[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2021, 31(8): 2923-2937.

[27] Zhang Z Y, Zheng L L, Piao Y J, et al. Blind remote sensing image deblurring using local binary pattern prior[J]. Remote Sensing, 2022, 14(5): No.1276.

[28] Li Z Y, Guo J Y, Zhang Y T, et al. Reference based multi-level features fusion deblurring network for optical remote sensing images[J]. Remote Sensing, 2022, 14(11): No.2520.

[29] Lim H, Yu S, Park K, et al. Texture aware deblurring for remote sensing images using L0-based deblurring and L2-based fusion[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 3094-3108.

[30] 王琪瑶, 胡琸悦, 李潇雁, 基于局部最大和最小强度先验的遥感图像盲去模糊[J]. 激光与光电子学进展, 2023, 60(4): 429-436.

[31] Wang Qi-yao, Hu Zhuo-yue, Li Xiao-yan, et al. Blind deblurring of remote sensing images based on local maximum and minimum intensity prior[J]. Laser and Optoelectronics Progress, 2022, 60(4): 429-436.

[32] 刘晨辉, 尹增山, 高爽. 基于遥感图像序列的运动去模糊算法研究[J]. 激光与光电子学进展, 2022, 59(8): 496-504.

[33] Liu Chen-hui, Yin Zeng-shan, Gao Shuang. Motion deblurring algorithm based on remote sensing image sequence[J]. Laser and Optoelectronics Progress, 2022, 59(8): 496-504.

[34] Fergus R, Singh B, Hertzmann A, et al. Removing camera shake from a single photograph[J]. ACM Transactions on Graphics, 2006, 25(3): 787-794.

[35] Levin A, Weiss Y, Durand F, et al. Understanding and evaluating blind deconvolution algorithms[C]//IEEE Conference on Computer Vision and Pattern Recognition, Miami, USA, 2009: No.10836014.

[36] Krishnan D, Fergus R. Fast image deconvolution using hyper-Laplacian priors[J]. Advances in Neural Information Processing Systems, 2009, 22: 1033-1041.

基金资助

国家自然科学基金项目(12171054)

吉林省教育厅科学技术研究项目(JJKH20230791KJ)

吉林省教育厅科学技术研究项目(JJKH20230788KJ)