结合L-Like曲线法与正则化模型实现多帧太阳斑点图像重建

Combining L-Like Curve Method with Regularization Model to Reconstruct Multi Frame Solar Speckle Images

摘要: 太阳斑点图像重建是天文观测领域中一个重要的研究问题。通过地基光学望远镜获取的天文图像受大气湍流和大气扰动的影响，会发生严重的模糊或降质，需要借助图像复原方法进行重建。正则化方法是图像盲复原的经典算法，它通过构造正则项实现模糊核与清晰图像的迭代求解，现有大多数正则化方法都是针对单幅图像进行处理，图像特征越多，重建效果越好，而单幅太阳斑点图由于特征不足会导致重建图像质量较差。本文结合多帧斑点图之间的互补关系建立了适合太阳斑点图重建的多帧盲复原模型，并结合L-like曲线法，对正则化参数进行计算。实验结果表明，本文方法能较好地实现太阳斑点图的复原，满足天文观测的要求。

Abstract: Solar speckle image reconstruction is an important research issue in the field of astronomical observation. Astronomical images acquired through ground-based optical telescopes will be affected by atmospheric turbulence and atmospheric disturbances, and will be seriously blurred or degraded, and need to be reconstructed using image restoration methods. The regularization method is a classic algorithm for blind image restoration. It constructs regularization terms to achieve iterative solution of fuzzy kernels and clear images. Most of the existing regularization methods deal with a single image. The more image features, the better the reconstruction effect, and the single solar speckle image could cause poor quality of the reconstructed image due to insufficient features. In this paper, a multi-frame blind restoration model suitable for the reconstruction of solar speckle images is established based on the complementary relationship between the multi-frame speckle images, and combined with the L-like curve method, the regularization parameters are calculated. Experimental results show that the proposed method can better restore the solar speckle image and meet the requirements of astronomical observation.

1. 引言

太阳斑点图像重建是天文图像处理的重要研究内容之一。在太阳高分辨观测中，受大气湍流和大气扰动的影响，利用地基光学望远镜获取的天文图像会发生严重的模糊或降质，需要进行高分辨图像重建 [1]。太阳斑点图重建技术可以分为两类：一类为频率重建算法，如斑点干涉术 [2]、斑点掩模法、Knox-Thompson法 [3] 等；另一类是直接在空间域中完成重建，如简单位移叠加法 [3] 和加权位移叠加法等。这些重建算法是基于大气物理模型，重建过程中需要较多的先验知识，如大气视宁度、斑点干涉函数等，并且需要较多的图像帧数，导致重建过程计算量大，不能满足天文观测的实时性业务处理需求。

太阳斑点图像中包含较多噪声信息，需要借助图像处理中的盲复原方法进行盲复原处理。盲复原是同时求解复原图像和未知的模糊核，未知量个数多余方程个数，属于严重不适定问题。正则化方法是Tikhonov等人在1977年提出 [4]，是求解不适定问题的经典方法，它是通过引入图像平滑性的正则化约束，将不适定问题转化成适定问题。Tikhonov等 [4] 利用高斯先验去拟合自然图像的重尾分布，构造关于图像梯度的正则项，对图像进行去模糊，虽然能够抑制噪声，但是去模糊效果不好，细节丢失严重；Rudin等人 [2] 提出了总变分正则化方法，使用拉普拉斯先验去拟合自然图像的梯度重尾分布，能够得到较好的复原效果，但会产生一定的振铃效应；Krishnan等人 [5] 提出了图像的标准化稀疏先验，对简单的模糊核进行估计；唐梦等人 [6] 在此基础上，首先对图像进行滤波预处理，降低噪声和突出边缘，再使用总变分正则化方法对图像进行非盲去卷积处理，计算过程较为复杂。M. Van Noort [7] 等人将盲反卷积应用于天文湍流退化图像的复原中，提出了太阳斑点图的盲反卷积算法。2006年，Y. V. Zhulina [8] 将极大似然估计图像复原算法和Ayers-Dainty的IBD算法相结合，提出了一种多帧迭代盲解卷积算法(MIBD-Multi-Frame Iterative Blind Deconvolution) [5]，该算法原理简单，能够处理多种不同类型PSF引起的图像衰退，同时在迭代过程中除了假设PSF和图像灰度值非负之外，不要求其他先验知识；但是，该算法收敛速度缓慢，并且对噪声非常敏感。

太阳斑点图的特征较为单一 [8]，如果采用传统的单帧重建算法会导致重建图像分辨率较低，复原效果较差。多帧重建算法可以充分利用图像帧间包含相似但不完全相同的信息进行互补，使重建图像获得更多的细节和信息，更加接近理想图像。

本文在Tikhonov经典正则化模型基础上，利用多帧斑点图之间的连续性特征建立了适合太阳斑点图重建的多帧盲复原模型，结合L-like曲线法，对模型中的正则化参数进行计算，并将其作为正则化的初始参数 [9]，经过迭代计算得到模糊核和重建图像。实验结果表明，本文方法计算量不大，只需少量帧就能复原出满足天文观测业务需求的清晰图像。

本文结构如下：第2部分对多帧图像盲复模型的建立过程及计算方法进行详细描述；第3部分结合L-like曲线法对模型中的正则化参数进行计算；第4部分给出相关实验结果，说明本文方法的可行性。

2. 多帧斑点图盲复原模型的建立

2.1. 单帧盲复原模型

在单帧图像盲复原计算中，最具代表性的正则化方法是Tikhonov在1977年提出的L₂正则约束，其模型如下 [6]：

$\underset{f,h}{\arg \min }\left({\Vert f\ast h-g\Vert }^2+\alpha {\Vert h\Vert }^2+\beta B\left(f\right)\right)$ (1)

其中f和h分别是待复原图像和待估计模糊核，g为观测的模糊图像，正则项 ${\Vert f\ast h-g\Vert }^2$ 为数据保真项，其值越小表明复原图像与原始清晰图像越接近；同样，模糊核的正则项 ${\Vert h\Vert }^2$ 越小，表明估计的模糊核越接近真实的模糊核， $B\left(f\right)$ 是图像f的其它信息，根据不同的先验知识进行定义， $\alpha$ 和 $\beta$ 为正则化参数，起着平衡正则项和数据项的作用。

在基于图像灰度和图像梯度的正则化模型(1)中， $B\left(f\right)$ 定义为 [10]：

$B\left(f\right)=\sigma P_0\left(f\right)+P_0\left(\nabla f\right)$ (2)

其中 $P_0\left(f\right)={\Vert f\Vert }_0$ 和 $P_0\left(\nabla f\right)={\Vert \nabla f\Vert }_0$ 分别是图像f和梯度图像 $\nabla f$ 中非零元素的个数， $\sigma$ 是权重系数。结合式(1)和(2)，图像盲复原的正则化模型可表示为：

$\underset{f,h}{\arg \min }\left\{{\Vert f\ast h-g\Vert }^2+\alpha {\Vert h\Vert }^2+\beta \sigma P_0\left(f\right)+P_0\left(\nabla f\right)\right\}$ (3)

图像复原的过程实质上是矩阵的反卷积运算，由于单帧太阳斑点图所包含的细节特征不足，在重建运算过程中噪声被放大，导致重建结果不理想，如图1所示。

(a) 模糊图像 (b) 利用式(3)重建的图像

Figure 1. Reconstruction result of single speckle image

图1. 单帧图重建结果

2.2. 多帧盲复原模型的推导

在天文观测中，太阳斑点图像的形成过程在等晕区是线性的或空间不变的，噪声主要是加性噪声，且独立于观测目标。等晕区里多帧观测图像的成像模型可表示为：

$g_i=f\ast h_i+n_i,i=1,2,3,\cdots ,k$ (4)

其中f是原始图像， $h_i$ 表示第i帧的模糊核， $n_i$ 是第i帧图像中的噪声， $g_i$ 是第i帧模糊图像。

对于同一观测目标的多帧图像，如果不考虑噪声的存在，每两帧图像有如下关系：

$g_i\ast h_j=g_j\ast h_i,i\ne j$ (5)

但在实际天文观测中，噪声存在并且是多种噪声的混合，所以式(5)中的 $h_i$ 并不是真实的点扩散函数，对 $h_i$ 加上修正算子s进行修正 [7] (s待定)，式(5)应写为：

$g_i\ast h_j\ast s=g_j\ast h_i\ast s,i\ne j$ (6)

综合式(4)和(6)，多帧图像退化模型为：

$\{\begin{array}{l}{g}_i=f\ast h_i+n_i,i=1,2,\cdots ,k\\ g_i\ast h_j\ast s=g_j\ast h_i\ast s,i\ne j\end{array}$

令

$A\left(\left\{h_i\right\}\right)=\underset{m,n\in i,m\ne n}{{\displaystyle \sum }}{\Vert g_m\ast h_n-g_n\ast h_m\Vert }^2$ (7)

其中 $\left\{h_i\right\}$ 表示多帧图像的模糊核。式(7)表示多帧图像复原后的差异，将其作为新的正则化约束项，在式(3)基础上得到多帧图像重建的正则化模型：

$\underset{f,\left\{h_i\right\}}{\text{argmin}}\left(\underset{i}{{\displaystyle \sum }}{\Vert f\ast h_i-g_i\Vert }^2+\alpha A\left(\left\{h_i\right\}\right)+\beta \sigma P_0\left(f\right)+P_0\left(\nabla f\right)\right)$ (8)

其中参数 $\alpha$， $\beta$ 是相应正则化项的权重。

对式(8)的求解分为如下两个子问题：

子问题1：求解模糊核 $\left\{h_i\right\}$

$\underset{\left\{h_i\right\}}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert f\ast h_i-g_i\Vert }^2+\alpha A\left(\left\{h_i\right\}\right)\right)$ (9)

子问题2：求解原始图像f

$\underset{f}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert f\ast h_i-g_i\Vert }^2+\beta \left(\sigma P_0\left(f\right)+P_0\left(\nabla f\right)\right)\right)$ (10)

以下对子问题的求解过程进行推导。

1) 问题1求解。在式(9)中，正则化约束 $A\left(\left\{h_i\right\}\right)$ 中的 $h_i$ 会受到噪声的干扰并不是真实的模糊核，根据式(6)，将其改写为

$A\left(\left\{h_i\right\}\right)=\underset{m,n\in i,m\ne n}{{\displaystyle \sum }}{\Vert g_m\ast h_n\ast s-g_n\ast h_m\ast s\Vert }^2$

由于s未知，问题1转化为求解 $\left\{h\right\},s$ 使得

$\underset{\left\{h_i\right\},s}{\text{argmin}}\left(\underset{i}{{\displaystyle \sum }}{\Vert f\ast h_i-g_i\Vert }^2+\alpha A\left(\left\{h_i\right\}\right)\right)$

成立。

首先根据 $h=\underset{h}{\arg \min }A\left(\left\{h_i\right\}\right)$ 得到h的初始估计 $h^0$，模糊核的退化模型可以看作：

$h^0=h\ast s+ϵ$ (11)

对式(11)的求解也可以使用正则化方法。加入正则化约束项构建目标优化问题：

$\underset{s,\left\{h_i\right\}}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert h_i\ast s-h_i^0\Vert }^2+\delta \phi \left(s\right)+\delta \phi \left(\left\{h_i\right\}\right)+\gamma {\omega }_p\left(\left\{h_i\right\}\right)\right)$ (12)

其中 $\phi (\cdot )$ 代表正约 [4]：

$\phi \left(h\right)={\displaystyle {\sum }_{X}p}\left(h\left(X\right)\right),p\left(t\right)=\{\begin{array}{l}t,\text{if}t\ge 0\\ +\infty ,\text{otherwise}\end{array}$

${\omega }_p(\cdot )$ 代表 $L_p$ 范式。类似式(8)的求解，对式(12)交替求解 $\left\{h_i\right\}$ 和s：

模糊核 $\left\{h_i\right\}$ 的求解使用

$\underset{\left\{h_i\right\}}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert h_i\ast s-h_i^0\Vert }^2+\delta \phi \left(\left\{h_i\right\}\right)+\gamma {\omega }_p\left(\left\{h_i\right\}\right)\right)$

s的求解使用

$\underset{s}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert h_i\ast s-h_i^0\Vert }^2+\delta \phi (\; s\; )\right)$

2) 问题2求解。根据问题1得到的模糊核估计 $\left\{{\hat{h}}_i\right\}$ 和s，对原始图像f进行求解：

$\underset{f}{\arg \min }\left(\underset{i}{{\displaystyle \sum }}{\Vert f\ast {\hat{h}}_i-g_i\Vert }^2+\beta \left(\sigma P_0\left(f\right)+P_0\left(\nabla f\right)\right)\right)$ (13)

借助傅里叶变换，得到式(13)的解为

$f=F^{-1}(\frac{{\displaystyle \sum \overline{F\left({\stackrel{^}{h}}_i\right)}F\left(g_i\right)}}{{\displaystyle \sum \overline{F\left({\stackrel{^}{h}}_i\right)}F\left({\hat{h}}_i\right)}+\beta \overline{F\left(P\right)}F(\; P\; )})$

3. 利用L曲线法计算正则项参数

在本文多帧正则化模型中，正则化参数有 $\alpha ,\beta ,\gamma ,\delta ,\sigma$，这些参数的选取不仅影响复原结果，还会影响迭代计算过程的收敛性，从图2中可以看出，错误的正则化参数不仅达不到去模糊的效果，还会导致图像质量更差。

(a) 原图像 (b) 参数合适 (c) 参数太小

Figure 2. The influence of regularization parameters on restoration results

图2. 正则项参数取值对复原结果的影响

经典的正则化参数计算方法有广义交叉验证法和L曲线法。其中，L曲线法是以log-log尺度来描述残差范数与解的范数曲线对比，根据对数尺度图形中是否出现明显的L形状确定近似参数的方法 [7]，更适合于图像的计算。

L曲线法的计算步骤如下：

$\eta ={\Vert h_i\Vert }^2,\rho ={\Vert f\ast h_i-g_i\Vert }^2$

两边取对数，得到：

$\bar{\eta }=\lg \eta =2\lg \Vert h_i\Vert ,\bar{\rho }=\lg \rho =2\lg \Vert f\ast h_i-g_i\Vert$

L曲线就是由点 $\left(\frac{\bar{\eta }}{2},\frac{\bar{\rho }}{2}\right)$ 拟合而成的，其曲率

$k=\frac{{\bar{\rho }}^{\prime }{\bar{\eta }}^{″}-{\bar{\rho }}^{″}{\bar{\eta }}^{\prime }}{{\left({{\bar{\rho }}^{\prime }}^2+{{\bar{\eta }}^{\prime }}^2\right)}^{3/2}}$ (14)

的最大值所对应的点即为式(8)的正则项参数 $\alpha$。

虽然L-曲线很容易定义，但是在实际情况中，通过公式(14)计算曲率是非常复杂和困难的，并不适用于本文模型的参数求解。本文使用另外一种方法来求取正则化参数，同样是通过曲线描述的方式，由于其形状与L-曲线类似，因此称为L-like曲线法。

L-like曲线法将模糊图像作为原始图像，模糊核作为其纹理部分，对本文模型中的正则化参数进行计算。由于实际应用中的噪声影响，为了在去模糊的同时能够保留边缘和细节纹理，使用下式对正则化参数 $\alpha$ 进行求解：

$\alpha =\underset{i,j}{\max }\left\{{\left(\left(u-v\right)\nabla u\right)}_{i,j}\right\},⥂i=1,2,\cdots ,M;j=1,2,\cdots ,N$ (15)

其中u和v分别表示模糊图像的初始结构部分和初始纹理部分。

先对原图进行结构纹理分解，得到初始的 $\alpha$ 值再代入本文正则化模型中，得到复原图像，再进行迭代计算，得到实验结果如图3所示。

(a) 模糊图像 (b) 初始的L-like图 (c) 复原图像 (d) 稳定的L-like曲线图

Figure 3. Image restoration and its L-like relationship curve

图3. 图像去模糊及其L-like关系曲线图

4. 实验结果分析

为验证本文模型的有效性，我们从云南天文台提供的2019年3月3日13点30分获取的200帧1912 × 2360太阳斑点观测图像中，先对图像进行预处理，再利用斑点干涉选帧方法选出20帧，分别对第1帧斑点图进行单帧重建，对20帧斑点图进行多帧重建，得到重建结果如图4所示。

(a) 第1帧 (b) 第10帧

Figure 4. Solar speckle reconstruction image

图4. 太阳斑点重建结果图

第1帧太阳斑点图像重建结果仅局部清晰，米粒轮廓模糊；20帧多帧重建结果相较于单帧重建结果轮廓大部分清晰，但边缘区域清晰度不够。通过实验我们发现这是由于相邻的两帧斑点图特征相似度不够所导致的。我们对同一组200帧太阳斑点图进行分块多帧重建，流程示意图如图5所示，重建步骤如下：

1. 通过计算每帧图像的功率谱大小，选取功率谱最大的一帧斑点图作为参考帧；

2. 对一组斑点图进行分块处理，分块大小为500 × 600；

3. 对于每一子块，通过计算每帧与参考帧之间的图像相似度ssim，并选取图像相似度最大的前10帧作重建图像，进行多帧重建；

4. 对每块重建结果进行拼接处理得到最终结果。

Figure 5. Multi-frame reconstruction process diagram

图5. 多帧重建过程中分块、重建再拼接流程示意图

(a) MIBD (b) 本文方法 (c) MIBD局部细节 (d) 本文方法局部细节

Figure 6. Reconstruction results of solar speckle image by different methods

图6. 不同方法对同组太阳斑点图像的重建结果

多帧分块重建结果如图6所示，我们在此基础上与多帧迭代盲反卷积算法，简单位移叠加算法进行对比重建实验。对于图像复原的效果，除了目视的主观评价之外，本文采用图像均方误差(MSE)，峰值信噪比(PSNR)，图像相似度(SSIM)来对其进行客观定量地评价，并且对两种算法的重建时间进行统计。表1给出了本文重建算法和其他不同重建算法对该组太阳斑点图像重建结果的客观评价指标及其重建所需时间。

Table 1. Evaluation of the first set of experimental results

表1. 该组实验结果的评价

通过以上实验结果可以看出，太阳斑点图的多帧盲复原结果相较于单帧盲复原结果从视觉上有明显的提升，细节部分得到增强，局部细节更加直观，能够较好地复原太阳斑点图，并且在时间复杂度上效果提升明显。综合表中的数据发现，本文方法具有比传统的多帧迭代盲反卷积算法，有更好的重建性能，在图像均方误差，峰值信噪比，图像相似度明显优于其余算法。

5. 总结

本文针对太阳斑点图的盲复原问题，研究单帧正则化盲复原算法，并在此基础上利用图像之间的关系建立了多帧图像重建模型，引入L-like曲线法计算正则化参数。实验表明，本文方法使用多帧图像包含的互补信息，能在实现图像重建的同时，保持边缘并抑制噪声,具有良好的抗噪性,有效地提高太阳斑点图重建质量，是一种耗时少、重建性能好的图像盲复原算法。现阶段本文算法，重建一个子块需要600s左右，时间效率较高，同时本文方法的重建结果图像中仍有部分区域存在一定程度的模糊，影响了重建结果的整体效果，需要在细节方面进一步改进，由于正则化模型的求解都是一个优化计算的过程，在分块重建过程中采用的是串行的方法，因而本文方法也存在算法时间效率不够理想的问题，需要进一步优化，下一步可以采用并行的方式对其进行重建加速。