0 引言

地基光电探测系统是对深空目标进行探测与识别的主要技术手段之一，目标的成像特征是进行识别的主要特性.受成像系统、大气环境、动态噪声等诸多因素的影响，获取的空间目标图像产生严重的退化效应，影响系统对目标的识别与检测效率，因此需要利用图像反演技术对退化后的空间目标图像进行波后复原.

近年来，针对图像盲反演问题，研究学者提出许多有效的算法.FERGUS R等^[1]提出一种针对相机运动退化图像的盲图像反演方法，该方法利用混合高斯先验对图像进行约束，采用变分贝叶斯方法估计出模糊核，最后使用Richardson-Lucy算法反演出原始的图像.SHAN Q等^[2]提出了基于最大后验概率(Maximum A Posteriori，MAP)的图像盲反演方法采用连续的分段函数对图像梯度分布的重尾特性进行拟合，该方法可以恢复出图像中的复杂结构和边缘细节.KRISHNAN D等^[3]提出一种对图像梯度采用归一化稀疏先验的总变分盲反演方法，该方法可以快速反演出原始图像和点扩散函数.针对总变分盲反演模型进行求解时容易陷入局部最优解的问题，PERRONE D等^[4]提出一种新的迭代更新求解总变分反演模型的方法，可提高点扩散函数和图像的反演准确度.PAN J等^[5]提出了利用图像强度和梯度进行稀疏约束的盲反演方法，该方法对文本图像和自然图像都有较好的反演效果，但是运算速度较慢，耗时较长.LIN Z等^[6]针对图像中存在异常值时影响图像和点扩散函数估计效果的问题，提出了一种基于MAP框架的图像盲反演方法.该方法首先利用图像的显著边缘对模糊核进行估计，然后采用L_p范数对图像梯度进行约束从而实现图像反演，实验证明该方法对于非均匀退化图像和存在异常值的图像均具有较好的反演效果.郭从洲等^[7]根据模糊核的稀疏特性和自然光学图像的先验信息，提出了一种非凸高阶全变差的图像盲反演方法，对图像边缘保持效果较好.高伟哲等^[8]提出了一种基于指数律分布的快速图像盲反演方法，利用自然图像的指数律分布重建图像的频谱，采用多方向综合估计减少噪声，对G类点扩散函数造成的图像退化有较好的反演效果.仇翔等^[9]通过分析图像梯度分布和暗通道的稀疏特性，提出了一种基于L₀稀疏先验的改进正则化图像盲反演方法，采用近似线性映射矩阵来解决非凸优化最小化问题，该方法可有效抑制反演图像中的振铃效应.上述基于正则化的方法通过利用不同的先验信息构造盲反演模型中的正则约束项，并采用不同的求解算法对模型进行求解，从而估计出模糊核和原始图像，但还存在着对噪声敏感等问题.

随着稀疏表示理论的发展，进年来基于稀疏表示的图像盲反演方法取得了一系列进展.HU Z等^[10]利用图像自身的稀疏特性，提出一种自适应字典学习的单帧图像盲反演算法，但该方法仅适用于模糊核尺寸较小的退化图像.ZHANG H等^[11]提出一种基于字典稀疏表示的图像反演算法，该方法利用图像块的稀疏特性和字典的冗余特性可以有效反演出空间目标信号的结构特征，但字典稀疏表示技术需要对图像进行分块处理，利用反演后的图像块对图像进行重构时容易产生“边界伪像”.YU J等^[12]提出利用图像的稀疏表示和结构自相似特性实现模糊核估计和图像反演的盲反演方法.上述方法均能有效地反演出图像的结构特征，但对图像细节信息的重构能力有限.

本文针对目前基于稀疏表示的图像盲反演算法未考虑空间目标图像自身的特性，存在对空间目标图像细节信息恢复不理想、重构图像中易产生边界伪像等不足之处，利用空间目标图像自身的稀疏特性及其梯度图像的稀疏特性在反演模型中增加正则约束项，提出了一种基于稀疏表示的联合稀疏先验约束空间目标图像盲反演算法.模型中采用图像梯度的L₀范数和稀疏表示构造代价函数中的图像正则项，在提取有利于模糊核估计的图像显著边缘信息的同时，可有效提升算法的鲁棒性.模糊核具有稀疏性，采用L_p范数对模糊核进行约束，以提高模糊核的估计准确度.采用L₀范数和超拉普拉斯先验分别对图像在空域和梯度上进行稀疏约束，从而进一步恢复出图像中的纹理和边缘等细节信息.

1 联合稀疏先验约束盲反演方法

空间目标图像退化过程可以建模为

1 $ \mathit{\boldsymbol{g}} = \mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{n}} $

式中，u为原始的空间目标图像，g为退化后的空间目标图像，n为噪声，k为模糊核，$ \otimes $为卷积运算符.空间目标图像盲反演过程是利用观测到的退化图像g来推测原始的清晰图像u，由于空间目标图像退化过程中模糊核是未知的，因此需要利用模糊核和原始图像的先验信息构造盲反演模型中的正则约束项，从而对二者进行估计.

基于稀疏表示的空间目标图像反演方法利用图像自身的稀疏特性和字典的冗余特性可有效反演出目标信号的结构特征，同时对噪声具有较好的鲁棒性.其模型为

2 $ \left( {\mathit{\boldsymbol{\hat u}},\mathit{\boldsymbol{\hat k}}} \right) = \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u}},\mathit{\boldsymbol{k}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \omega \sum\limits_{i = 1}^L {{{\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|}_1}} + \lambda J\left( \mathit{\boldsymbol{u}} \right) + \gamma R\left( \mathit{\boldsymbol{k}} \right) $

式中，$\mathit{\boldsymbol{\hat u}}$为反演图像，$\mathit{\boldsymbol{\hat k}}$为模糊核估计值，i表示图像块序号，L为图像中总的图像块数量；R_i是提取图像块的矩阵；α_i为图像块的稀疏表示系数；过完备字典D∈R^n×k，通过K-SVD字典学习方法训练得到；函数J(u)和R(k)是利用原始空间目标图像u和模糊核k的先验信息构造的正则约束项；η、ω、ξ、μ、λ、γ为正则项系数.等式右边的第一项为数据保真项，用来保证反演图像尽可能地接近于原始图像；第二项和第三项用来保证反演图像块u_i在过完备字典D下的稀疏表示系数具有稀疏性；第四项为图像正则约束项；第五项为模糊核正则约束项.图像正则约束项和模糊核正则项需要利用图像和模糊核的先验信息进行构造，不同先验信息决定着反演图像质量的好坏和模糊核估计准确度.

处于深空背景下的空间目标图像，其大部分区域为平稳的黑色背景，而空间目标明亮且应具有明显的边缘和轮廓.图 1为深空背景下卫星、空间站图像及其强度分布直方图、梯度分布直方图.从强分布直方图 1(b)和(f)可看出，大部分像素强度为零，部分像素聚集在较大强度处，表明强度分布具有稀疏性.梯度分布直方图 1(d)和(h)的大部分区域为零，目标边缘等少部分区域取值较大，表明图像梯度具有稀疏特性.图 2为短、长两种曝光时间下的空间目标图像，及其强度和梯度分布直方图.对比强度分布直方图 2(b)和(f)可以看出，相比于短曝光时间，长曝光时间的图像强度整体上往高亮度偏移，但仍具有稀疏性，即图像强度分别主要集中在背景和目标处.相比于短曝光时间梯度分布直方图 2(d)，长曝光时间梯度分布直方图 2(h)没有显著变化，都具有稀疏特性.图 3是以地球为背景的空间目标图像，以及其强度和梯度分布直方图.强度分布直方图 3(b)为非零值的连续分布，分散分布在各强度级，表明该以地球为背景的空间目标图像强度不具有稀疏性；图 3(c)和图(d)表明其梯度也不具有很好的稀疏性.因此本文主要指深空背景下的空间目标.根据文献[13]可知，图像的细节信息不利于模糊核估计，而图像显著边缘具有极强的稀疏特性，可有效提高模糊核估计的准确度.图像梯度的L₀范数^[14]在平滑图像细节信息的同时，可以有效地保护图像边缘的稀疏性，因此采用图像梯度的L₀范数对图像梯度进行约束.

图 1. Intensity histogram and gradient histogram images of space object in deep space background

Fig. 1. Intensity histogram and gradient histogram images of space object in deep space background

下载图片 查看所有图片

图 2. Intensity histogram and gradient histogram images of space object at different exposure levels

Fig. 2. Intensity histogram and gradient histogram images of space object at different exposure levels

下载图片 查看所有图片

图 3. Intensity histogram and gradient histogram images of space object in ground background

Fig. 3. Intensity histogram and gradient histogram images of space object in ground background

下载图片 查看所有图片

在空间目标图像反演过程中，构建目标函数中正则项时选择的图像先验越能体现出原始图像梯度分布，反演的空间目标图像越接近于原始图像.图像反演中常采用高斯分布、Laplacian分布、Hyper-Laplacian分布等不同的分布先验对图像梯度分布的重尾特性进行拟合.图 4给出了三种先验对空间目标图像梯度分布的拟合效果，从图中可看出，原始空间目标图像的梯度分布具有稀疏特性，且Hyper-Laplacian分布^[15]对清晰的原始图像梯度的重尾分布特性拟合效果最好.因此采用Hyper-Laplacian分布先验作为空间目标图像的梯度分布先验，使用L_p(0 < p < 1)范数构造图像正则项对图像梯度分布的稀疏特性进行约束.

图 4. Fitting the gradient distribution of space object image with each prior

Fig. 4. Fitting the gradient distribution of space object image with each prior

下载图片 查看所有图片

因此，为了使反演图像具有较强的对比度且包含大量的细节信息，采用L_p范数和L₀范数对图像的梯度分布和空间域进行稀疏约束.为了获取有利于模糊核估计的图像显著边缘，从而提高模糊核估计准确度，采用图像梯度的L₀范数对图像梯度进行约束；为了保证模糊核的稀疏特性，采用L_p范数对模糊核进行约束，结合以上先验信息，以及模糊核元素的非负性和归一化特性，提出基于稀疏表示的联合稀疏先验约束图像盲反演模型，即

3 $ \begin{array}{*{20}{c}} {\left( {\mathit{\boldsymbol{\hat u}},\mathit{\boldsymbol{\hat k}}} \right) = \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u}},\mathit{\boldsymbol{k}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \omega \sum\limits_{i = 1}^L {{{\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|}_1}} + \xi {{\left\| \mathit{\boldsymbol{u}} \right\|}_0} + \mu \left\| {\nabla \mathit{\boldsymbol{u}}} \right\|_p^p + }\\ {\lambda {{\left\| {\nabla \mathit{\boldsymbol{u}}} \right\|}_0} + \gamma \left\| \mathit{\boldsymbol{k}} \right\|_p^p} \end{array} $

式中，等式右边的第四项对图像在空域上进行稀疏约束，以保证反演图像像素点具有显著对比度.第五项对图像梯度分布的稀疏特性进行约束，以保证反演图像包含丰富的细节信息.‖‖₀为L₀范数，‖‖_p^p为L_P范数.第六项为图像梯度的L₀正则项，可保留梯度图像的强边缘信息.$\nabla \mathit{\boldsymbol{u}}$为图像u的梯度，$\nabla $为梯度算子，水平方向的梯度算子为${\nabla _{\rm{h}}} = \left[ {1, - 1} \right]$，竖直方向的梯度算子为${\nabla _{\rm{v}}} = {\left[ {1, - 1} \right]^{\rm{T}}}$，实际中求取梯度通常用小区域模板卷积来近似计算，对水平梯度G_x和竖直梯度G_y各用一个模板，${\nabla _{\rm{h}}} = \left[ {1, - 1} \right]$为水平方向模板，${\nabla _{\rm{v}}} = {\left[ {1, - 1} \right]^{\rm{T}}}$为竖直方向模板.对梯度${\nabla _{\rm{u}}} = \left( {{G_x}, {G_y}} \right)$取模得到梯度大小.第七项为模糊核的正则项，用L_p范数进行约束，可保证模糊核的稀疏特性.η、ω、ξ、μ、λ、γ为正则项系数.

式(3)所示的空间目标图像盲反演模型中包含多个未知变量，可分解关于稀疏表示系数估计、模糊核估计和图像复原的子问题进行迭代求解，从而估计出模糊核和原始的空间目标图像.

2 超完备字典的稀疏表示

在稀疏表示模型中，一幅完整的图像u通常情况下可以被分解为许多图像块，通过图像块的稀疏表示可实现原图像的重构^[16].假设给定一幅大小为$\sqrt N \times \sqrt N $的图像u，过完备字典D∈ R^n×k，字典由大小为$\sqrt n \times \sqrt n $的图像块通过K-SVD字典学习方法训练得到.则图像u中的第i个图像块u_i可以用过完备字典D进行稀疏表示，即

4 $ {\mathit{\boldsymbol{u}}_i} = {\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i} $

式中，R_i是提取图像块的矩阵，其大小为n×N；α_i为图像块的稀疏表示系数.整幅图像u可表示为

5 $ \mathit{\boldsymbol{u}} = \sum\limits_{i = 1}^L {\mathit{\boldsymbol{R}}_i^{\rm{T}}\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} $

对于空间目标图像而言，由于其自身具有稀疏特性，因此将图像分解为各图像块进行稀疏表示时，稀疏表示系数α_i中的零元素越多，表明图像的稀疏性越强.过完备字典稀疏表示利用空间目标图像自身的稀疏特性可有效消除退化图像中的噪声信号，从而提高反演算法对噪声的鲁棒性.

通过式(3)对稀疏表示系数α_i进行求解时，固定原始图像u和模糊核k的值不变，可得到关于表示系数α_i的优化问题，表达式为

6 $ \left\{ {{{\mathit{\boldsymbol{\hat \alpha }}}_i}} \right\} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{\alpha }}_i}} \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \omega \sum\limits_{i = 1}^L {{{\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|}_1}} $

显然，式(6)是关于整幅原始图像u的稀疏表示矩阵的求解问题，因此可以将图像u分解为相对独立的图像块，于是将式(6)转化为对各个独立的图像块的稀疏系数进行求解，即

7 $ {{\mathit{\boldsymbol{\hat \alpha }}}_i} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{\alpha }}_i}} \eta \left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2 + \omega {\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|_1} $

式(7)可采用LARS算法^[17]进行快速求解.

图 5给出了空间目标图像在过完备字典D下的稀疏表示系数，从图中可看出，大小为8×8的空间目标图像块在过完备字典D下的稀疏表示系数中只有少部分不为0的值，图像可通过少量的原子进行线性组合来表示，因此证明了空间目标图像具有稀疏特性，可采用稀疏表示构造反演模型中的正则项.

图 5. Sparse representation of space object image

Fig. 5. Sparse representation of space object image

下载图片 查看所有图片

3 模糊核估计

为提高模糊核估计准确度，采用图像梯度的L₀范数对图像梯度进行约束，以获取有利于模糊核估计的图像显著边缘；采用L_p范数对模糊核进行约束，保证模糊核的稀疏特性，从而有效降低模糊核中的噪声.由式(3)可得关于模糊核估计的优化问题为

8 $ \left( {\mathit{\boldsymbol{\hat u}},\mathit{\boldsymbol{\hat k}}} \right) = \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u}},\mathit{\boldsymbol{k}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \omega \sum\limits_{i = 1}^L {{{\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|}_1}} + \lambda {\left\| {\nabla \mathit{\boldsymbol{u}}} \right\|_0} + \gamma \left\| \mathit{\boldsymbol{k}} \right\|_p^p $

式(8)可分解成关于u和k的子问题进行求解.对u子问题进行求解时，固定模糊核k和稀疏表示系数α_i的值不变，由式(8)可得关于u的优化问题为

9 $ \mathit{\boldsymbol{\hat u}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{u}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \lambda {\left\| {\nabla \mathit{\boldsymbol{u}}} \right\|_0} $

由于式(9)中有L₀正则化项，是一个非凸优化问题.因此采用半二次惩罚技术^[14]求解该模型.对$\nabla \mathit{\boldsymbol{u}}$引入辅助变量f=(u_x, u_x) ^T，因此式(9)可转化为

10 $ \left( {\mathit{\boldsymbol{\hat u}},\mathit{\boldsymbol{\hat f}}} \right) = \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u}},\mathit{\boldsymbol{f}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \beta \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{f}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{f}} \right\|_0} $

式中，β为正则约束项系数，当β→$\infty $时，式(10)的解将收敛于式(9)的解. u和f可以通过固定其余变量进行交替求解.

对f进行求解时，固定式(10)中的其余变量，可得f的代价函数为

11 $ \mathit{\boldsymbol{\hat f}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{f}} \beta \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{f}}} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{f}} \right\|_0} $

式(11)可通过分段的连续函数不断逼近进行求解^[14]，得

12 $ \mathit{\boldsymbol{f}} = \left\{ {\begin{array}{*{20}{c}} {\nabla \mathit{\boldsymbol{u}}}&{{{\left| {\nabla \mathit{\boldsymbol{u}}} \right|}^2} \ge \frac{\lambda }{\beta }}\\ 0&{{\rm{otherwise}}} \end{array}} \right. $

将得到的f固定不变，可得到每次迭代过程中u的代价函数为

13 $ \mathit{\boldsymbol{\hat u}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{u}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \beta \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{f}}} \right\|_2^2 $

式(13)是一个最小二乘问题，因此可通过快速傅里叶变换进行求解，即

14 $ \mathit{\boldsymbol{u}} = {{\rm{F}}^{ - 1}}\left( {\frac{{\overline {{\rm{F}}\left( \mathit{\boldsymbol{k}} \right)} {\rm{F}}\left( \mathit{\boldsymbol{g}} \right) + {\rm{F}}\left( {\eta \sum\limits_{i = 1}^L {\mathit{\boldsymbol{R}}_i^{\rm{T}}\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} } \right) + \beta {{\rm{F}}_f}}}{{\overline {{\rm{F}}\left( \mathit{\boldsymbol{k}} \right)} {\rm{F}}\left( \mathit{\boldsymbol{k}} \right) + {\rm{F}}\left( {\eta \sum\limits_{i = 1}^L {\mathit{\boldsymbol{R}}_i^{\rm{T}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{I}}} } \right) + \beta \overline {{\rm{F}}\left( \nabla \right)} {\rm{F}}\left( \nabla \right)}}} \right) $

式中，F(·)表示快速傅里叶变换；F^－1(·)表示快速傅里叶反变换；F(·)为F(·)的复共轭；${{\rm{F}}_f} = \overline {{\rm{F}}\left( {{\nabla _h}} \right)} {\rm{F}}\left( {{f_h}} \right) + \overline {{\rm{F}}\left( {{\nabla _v}} \right)} {\rm{F}}\left( {{f_v}} \right)$，${{\nabla _h}}$和${{\nabla _v}}$分别为水平方向和竖直梯度算子.

对模糊核进行估计时，利用式(9)得到的包含显著边缘的空间目标梯度图像后，由式(8)得到关于k的优化问题为

15 $ \begin{array}{l} \mathit{\boldsymbol{\hat k}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{k}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \gamma \left\| \mathit{\boldsymbol{k}} \right\|_p^p\\ {\rm{subject}}\;{\rm{to}}\;\;\mathit{\boldsymbol{k}} \ge 0,{\left\| \mathit{\boldsymbol{k}} \right\|_1} = 1 \end{array} $

p的取值范围为[0.5, 1]，本文中p取值为0.5；从文献[18]可知，直接用图像估计出的模糊核准确度不高，因此为了获得更加准确的模糊核估计值${\mathit{\boldsymbol{\hat k}}}$，采用梯度图像对模糊核进行估计，其数学模型为

16 $ \begin{array}{l} \mathit{\boldsymbol{\hat k}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{k}} \left\| {\mathit{\boldsymbol{k}} \otimes \nabla \mathit{\boldsymbol{u}} - \nabla \mathit{\boldsymbol{g}}} \right\|_2^2 + \gamma \left\| \mathit{\boldsymbol{k}} \right\|_p^p\\ {\rm{subject}}\;{\rm{to}}\;\;\mathit{\boldsymbol{k}} \ge 0,{\left\| \mathit{\boldsymbol{k}} \right\|_1} = 1 \end{array} $

式(16)是L_p范数优化问题，因此函数非凸的，可采用迭代重加权最小二乘(Iterative Re-weighted Least Squares，IRLS)算法^[19]进行求解.

通过交替迭代求解式(6)、(9)、(16)直至算法收敛，即可估计出空间目标退化图像的模糊核.

4 图像反演

为了保证反演图像的像素点间具有显著的对比度，同时保证图像中包含边缘和纹理等细节信息，分别采用L_p范数和L₀范数对图像的梯度分布和空间域进行稀疏约束，由式(3)可得图像反演阶段得模型为

17 $ \mathit{\boldsymbol{\hat u}} = \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u}},\mathit{\boldsymbol{k}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \omega \sum\limits_{i = 1}^L {{{\left\| {{\mathit{\boldsymbol{\alpha }}_i}} \right\|}_1}} + \xi {\left\| \mathit{\boldsymbol{u}} \right\|_0} + \mathit{\boldsymbol{\mu }}\left\| {\nabla \mathit{\boldsymbol{u}}} \right\|_p^p $

式(17)所示的空间目标图像的反演模型包含两个未知变量u和α_i，可采用交替迭代更新的策略进行求解.因此可将式(17)分解为关于u和α_i的两个子问题进行求解，通过固定其中一个变量，对另一个变量进行迭代更新.其中α_i子问题的求解方法式(7)已给出，下面介绍u子问题的求解方法.

固定α_i的值不变，对u进行求解，通过式(17)可得到关于u的子问题为

18 $ \mathit{\boldsymbol{\hat u}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{u}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \xi {\left\| \mathit{\boldsymbol{u}} \right\|_0} + \mu \left\| {\nabla \mathit{\boldsymbol{u}}} \right\|_p^p $

式(18)是一个非凸优化问题，分别对u和${\nabla _\mathit{\boldsymbol{u}}}$引入两个新的辅助变量ρ和z(z=(u_x, u_x) ^T)，可将代价函数式(18)转化为

19 $ \arg \mathop {\min }\limits_{\mathit{\boldsymbol{u,\rho ,z}}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + \xi {\left\| \rho \right\|_0} + v\left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\rho }}} \right\|_2^2 + \omega \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{z}}} \right\|_2^2 + \mathit{\boldsymbol{\mu }}\left\| \mathit{\boldsymbol{z}} \right\|_p^p $

式中，ν和ω为正则项惩罚系数，当ν→$\infty $和ω→$\infty $时，式(19)的解将收敛于式(17)的解.固定ν和ω的值不变，交替更新u、ρ和z可得到式(19)的最优解.

固定ρ和z，关于u的优化问题为

20 $ \arg \mathop {\min }\limits_\mathit{\boldsymbol{u}} \left\| {\mathit{\boldsymbol{k}} \otimes \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}} \right\|_2^2 + \eta \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} \right\|_2^2} + v\left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\rho }}} \right\|_2^2 + \omega \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{z}}} \right\|_2^2 $

式(20)可通过快速傅里叶变换进行求解，得

21 $ \mathit{\boldsymbol{u}} = {{\rm{F}}^{ - 1}}\left( {\frac{{\overline {{\rm{F}}\left( \mathit{\boldsymbol{k}} \right)} {\rm{F}}\left( \mathit{\boldsymbol{g}} \right) + {\rm{F}}\left( {\eta \sum\limits_{i = 1}^L {\mathit{\boldsymbol{R}}_i^{\rm{T}}\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\alpha }}_i}} } \right) + \omega {\rm{F}}\left( \nabla \right){\rm{F}}\left( \mathit{\boldsymbol{z}} \right) + v{\rm{F}}\left( \mathit{\boldsymbol{\rho }} \right)}}{{\overline {{\rm{F}}\left( \mathit{\boldsymbol{k}} \right)} {\rm{F}}\left( \mathit{\boldsymbol{k}} \right) + {\rm{F}}\left( {\eta \sum\limits_{i = 1}^L {\mathit{\boldsymbol{R}}_i^{\rm{T}}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{I}}} } \right) + \omega \overline {{\rm{F}}\left( \nabla \right)} {\rm{F}}\left( \nabla \right) + v}}} \right) $

固定ρ和u，关于z 的优化问题为

22 $ \mathit{\boldsymbol{\hat z}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{z}} \omega \left\| {\nabla \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{z}}} \right\|_2^2 + \mu \left\| \mathit{\boldsymbol{z}} \right\|_p^p $

对式(22)进行求导，可得到四次多项式为

23 $ {\mathit{\boldsymbol{z}}^4} - 3\left( {\nabla \mathit{\boldsymbol{u}}} \right){\mathit{\boldsymbol{z}}^3} + 3{\left( {\nabla \mathit{\boldsymbol{u}}} \right)^2}{\mathit{\boldsymbol{z}}^2} - {\left( {\nabla \mathit{\boldsymbol{u}}} \right)^3}\mathit{\boldsymbol{z}} + \frac{{{\mu ^5}}}{{27{\omega ^3}}} = 0 $

式(23)中四次多项式的解析解为式(22)的最优解，可采用文献[15]给出的方法进行求解.

固定u和z，关于ρ的优化问题为

24 $ \arg \mathop {\min }\limits_\mathit{\boldsymbol{\rho }} v\left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\rho }}} \right\|_2^2 + \xi {\left\| \mathit{\boldsymbol{\rho }} \right\|_0} $

式(24)可采用文献[14]的方法进行求解，得

25 $ \mathit{\boldsymbol{\rho }} = \left\{ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{u}}&{{{\left| \mathit{\boldsymbol{u}} \right|}^2} \ge \frac{\xi }{v}}\\ {\bf{0}}&{{\rm{otherwise}}} \end{array}} \right. $

通过交替迭代求解式(29)、(22)、(24)直至算法收敛，可得到式(17)所示的混合正则约束模型的解，即反演出清晰的空间目标图像.

5 实验结果与分析

为了验证所提算法的有效性和适用性，首先采用不同退化类型和不同退化的模拟退化图像进行仿真实验，然后采用真实的空间目标退化图像进行实验，并将本文算法反演效果与KRISHNAN D等^[3]、ZHANG H等^[11]、PERRONE D等^[4]、PAN J等^[5]和LIN Z等^[6]所提的算法进行对比.实验过程中，对空间目标图像的像素值进行归一化处理，将图像分成8×8大小的图像块进行稀疏表示，过完备字典D采用K-SVD字典学习方法对空间目标图像训练得到.通过主观评价和客观评价相结合的方式对各算法的图像盲反演效果进行比较.在模拟退化图像反演实验部分，采用的客观评价指标为灰度平均梯度(Gray Mean Gradient, GMG)、结构相似度(Structural Similarity, SSIM)^[20]和清晰度；在真实退化图像反演实验部分，采用的客观评价指标为图像熵和灰度平均梯度(GMG).

5.1 客观评价指标

假设原始图像为u (i, j)，反演图像为${\mathit{\boldsymbol{\hat u}}}$(i, j)，退化图像为g (i, j)，图像大小均为M×N.

SSIM的取值范围为[0, 1]，SSIM的值越接近于1，表明反演图像同原始图像的结构越相似，即反演图像的结构保持越好.SSIM的表达式为

26 $ {\rm{SSIM}} = \frac{{\left( {2{\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{u}}}{\mu _{\mathit{\boldsymbol{\hat u}}}} + {C_1}} \right)\left( {2{\mathit{\boldsymbol{\sigma }}_{\mathit{\boldsymbol{u\hat u}}}} + {C_2}} \right)}}{{\left( {\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{u}}^2 + \mathit{\boldsymbol{\mu }}_{\mathit{\boldsymbol{\hat u}}}^2 + {C_1}} \right)\left( {\mathit{\boldsymbol{\sigma }}_\mathit{\boldsymbol{u}}^2 + \mathit{\boldsymbol{\sigma }}_{\mathit{\boldsymbol{\hat u}}}^2 + {C_2}} \right)}} $

式中，μ_u和μ_{${\hat u}$}分别表示原始图像和反演图像的灰度均值，σ_u和σ${\hat u}$分别为这两幅图像的标准差，σ_{u${\hat u}$} 为两幅图像的协方差，C₁和C₂为两个大于零且较小的常数.

灰度平均梯度(GMG)能较好地反映图像的边缘和纹理等细节信息的变化，GMG值越大，表明反演图像中的纹理和细节信息越丰富，图像的反演效果越好.其表达式为

27 $ {\rm{GMG}} = \frac{1}{{\left( {M - 1} \right)\left( {N - 1} \right)}}\sum\limits_{i = 1}^{M - 1} {\sum\limits_{j = 1}^{N - 1} {\sqrt {\frac{{{{\left[ {\mathit{\boldsymbol{\hat u}}\left( {i + 1,j} \right) - \mathit{\boldsymbol{\hat u}}\left( {i,j} \right)} \right]}^2} + {{\left[ {\mathit{\boldsymbol{\hat u}}\left( {i + 1,j} \right) - \mathit{\boldsymbol{\hat u}}\left( {i,j} \right)} \right]}^2}}}{2}} } } $

图像熵反映图像中所含的信息量，熵值越大，表明反演图像中的信息越丰富.其表达式为

28 $ E = - \sum\limits_{i = 0}^{\max } {{P_i}{{\log }_2}{P_i}} $

式中，P_i为灰度值i的像素在图像中出现的概率，max为图像的灰度范围(一般情况下为[0, 255]).

5.2 模拟退化目标图像反演实验

首先采用尺寸为11×11，标准差为1至5的高斯模糊对海事卫星图像和空间站图像进行模糊退化，然后加入均值为0、标准差为0.001的高斯噪声，从而得到模糊且含噪的退化图像.图 6为高斯模糊核的标准差为2时，各算法对高斯模糊的海事卫星图像的反演效果.图 7为高斯模糊的标准差为1时，各算法对高斯模糊的空间站图像的反演效果.从图 6和7可看出，同对比算法相比，本文算法估计出的模糊核更加接近于真实的模糊核，且具有较好的稀疏特性.同时，本文算法得到的反演图像具有更清晰的轮廓，以及丰富的纹理和细节特征，图像视觉效果较好.

图 6. 海事卫星图像高斯退化反演结果比较Compare of inversion results of Gaussian degraded maritime satellite image

Fig. 6. 海事卫星图像高斯退化反演结果比较Compare of inversion results of Gaussian degraded maritime satellite image

下载图片 查看所有图片

图 7. 空间站图像高斯退化反演结果比较Compare of inversion results of Gaussian degraded space station image

Fig. 7. 空间站图像高斯退化反演结果比较Compare of inversion results of Gaussian degraded space station image

下载图片 查看所有图片

表 1和表 2是高斯模糊核尺寸为11×11，标准差为1至5时各算法反演图像的SSIM和GMG.从实验结果可以看出，本文算法反演图像的SSIM和GMG整体上优于对比算法，反演图像中保留了更多的细节信息，且算法对退化严重的图像也同样有效.

采用模糊尺度为9、13、17、21和25，方向为45°的运动模糊对海事卫星图像和空间站图像进行模糊退化，然后加入均值为0、标准差为0.001的高斯噪声，从而得到模糊且含噪的退化图像.图 8是模糊核的模糊尺度为21时，各算法对运动模糊的海事卫星图像的反演效果.图 9是模糊核的模糊尺度为9时，各算法对运动模糊的空间站图像反演效果.从图中可看出，本文算法估计出的模糊核更接近于真实的模糊核.从各算法的反演图像可看出，本文算法反演图像不含“伪像”，图像中的空间目标具有更清晰的边缘和纹理信息.图 8和9证明了对于运动退化的空间目标图像，本文算法得到的反演图像也具有较好的主观视觉效果.

图 8. 海事卫星图像运动退化反演结果比较Compare of inversion results of motion degraded maritime satellite image

Fig. 8. 海事卫星图像运动退化反演结果比较Compare of inversion results of motion degraded maritime satellite image

下载图片 查看所有图片

图 9. 空间站图像运动退化反演结果比较Compare of inversion results of motion degraded space station image

Fig. 9. 空间站图像运动退化反演结果比较Compare of inversion results of motion degraded space station image

下载图片 查看所有图片

表 3和表 4是运动模糊核的模糊尺度为9、13、17、21和25，方向为45°时各算法反演图像的SSIM和GMG，从实验结果可以看出，随着运动模糊核的模糊尺度的不断变大，图像的退化程度不断加剧，图像中的边缘细节信息丢失越多，因此增加了图像的反演难度，各算法反演图像的客观评价指标的值总体上呈下降趋势.对于不同退化程度的运动退化图像，本文算法反演图像的各项评价指标的值总体上优于对比算法.

5.3 真实退化目标图像反演实验

为了进一步验证所提算法的适用性，采用图 10(a)所示的美国Meade公司生产的LX90-ACF望远镜观测到的月球表面的图像和图 11(a)所示的土星退化图像进行实验，选取图像熵和GMG对反演图像的质量进行客观评价，这两个评价指标的值越高说明图像的细节越丰富，图像的质量也就越高.各算法对观测到的真实空间目标退化图像的反演图像如图 10和11所示.表 5给出了各算法得到的反演图像的客观评价结果.

图 10. 月球观测图像反演结果比较Compare of inversion results of lunar observation image

Fig. 10. 月球观测图像反演结果比较Compare of inversion results of lunar observation image

下载图片 查看所有图片

图 11. 土星退化图像反演结果比较Compare of inversion results of real Saturn degraded image

Fig. 11. 土星退化图像反演结果比较Compare of inversion results of real Saturn degraded image

下载图片 查看所有图片

从图 8、9和表 5可看出：本文算法反演得到的空间目标图像轮廓和纹理等细节信息清晰可辨，与对比算法相比，反演图像看起来更自然，图像的整体视觉效果更好.

6 结论

本文提出了一种基于稀疏先验的联合稀疏先验约束盲反演算法.在盲反演模型中采用稀疏表示和L₀范数对图像进行约束，可保证算法在图像退化严重时也能获得较好的图像边缘；分别对图像在梯度域和空间域进行稀疏约束，从而保证反演图像具有丰富的纹理等细节特征；采用L_p范数对模糊核进行约束，可保证模糊核的稀疏性.对本文算法与现有算法的反演效果进行实验对比，结果表明：对于不同退化程度的高斯退化和运动退化的空间目标图像，本文算法的SSIM和GMG总体上都高于对比算法，反演图像具有更多的边缘和纹理等细节信息，图像质量更好.对于真实的空间目标退化图像，本文算法反演图像的图像熵和GMG也高于对比算法.实验结果证明本文算法对于模拟和真实的空间目标图像都具有较好的反演效果.