一种基于剪切波和特征信息检测的太阳斑点图融合算法

A Fusion Algorithm of Solar Image Based on Shear Wave and Feature Information Detection

摘要: 太阳斑点图像重建是天文观测领域中一个重要的研究问题。通过地基光学望远镜获取的天文图像受大气湍流和大气扰动的影响，会发生严重的模糊或降质。剪切波变换是一种多尺度几何分析方法，它比传统的小波变换更符合人类视觉系统的感知特性，能更有效地表示和捕获图像中的边缘、纹理等几何特征，并能充分利用图像自身的几何特性实现对其更为“稀疏”的表示。该方法用于太阳斑点图像融合时存在高频信息缺失、边缘模糊等问题。本文结合剪切波变换和特征信息检测，对多帧太阳斑点图像进行融合处理。首先通过一种简单的剪切波变换域图像融合方法得到一幅初始的融合图像，然后根据所有斑点图与初始融合图像像素间的局部相似程度来获得每一斑点图的有效特征信息区域检测图，并据此将该斑点图中的所有像素分成有效特征信息、有效特征信息与无效信息之间的过渡以及无效信息，最终据此来指导剪切波变换域各子带系数的融合。实验结果表明，本文方法能较好地恢复高频信息，实现太阳斑点图高分辨率重建。

Abstract: The reconstruction of solar image is an important research problem in the field of astronomical observation. Astronomical images acquired through ground-based optical telescopes are affected by atmospheric turbulence and atmospheric disturbances, and will be severely blurred or degraded. Shearlet transform is a multi-scale geometric analysis method. It is more in line with the perceptual characteristics of the human visual system than traditional wavelet transform. It can more effectively express and capture geometric features such as edges and textures in the image, and can make full use of the image. Its own geometric characteristics realize its more “sparse” representation. This method has problems such as lack of high-frequency information and blurred edges when it is used in the fusion of solar images. In this paper, combining shearlet transform and feature in-formation detection, multi-frame solar images are fused. First, an initial fusion image is obtained through a simple shearlet transform domain image fusion method, and then the effective feature information area detection map of each image is obtained according to the local similarity between all the images and the initial fusion image pixels, and according to this, all pixels in the image are divided into effective feature information, the transition between effective feature information and invalid information, and invalid information, which ultimately guides the fusion of the coefficients of each sub-band in the shearlet transform domain. The experimental results show that the method in this paper can recover high-frequency information well and achieve high-resolution reconstruction of the solar image.

1. 引言

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

多尺度几何分析方法能有效地表示和捕获图像中的边缘、纹理等几何特征，并能充分利用图像自身的几何特性实现对其“稀疏”的表示。与传统的小波变换相比，多尺度几何分析方法 [3] 在图像处理方面己经表现出了极大的优越性。目前，剪切波变换(shearlet transform) [4] [5] [6] 已经取代轮廓波变换(contourlet transform) [7]，成为基于多尺度几何分析的图像处理领域新的研究热点，在多传感器图像融合以及图像去噪中都得到了成功地应用。2012年，郑红等 [8] 将剪切波变换应用到红外与可见光圈像的融合中，提出了一种剪切波变换域红外图像与可见光图像的融合算法，该算法采用基于区域显著性与相似性制定的融合规则，最终所得的融合图像能完整、清晰地表达红外图像的热福射信息与可见光图像的细节信息；同年，Geng等 [9] 将剪切波变换与脉冲耦合神经网络结合起来，提出了一种结合剪切波变换与脉冲耦合神经网络(PCNN)的通用图像融合方法；2013年，Shi等 [10] 将剪切波变换与脉冲错合神经网络(PCNN)应用于遥感图像的融合中，提出了一种结合二者的变换域遥感图像融合方法；2014年，Liu等 [11] 将剪切波变换应用于多聚焦图像及医学图像的融合中，提出了一种基于剪切波变换域邻域特征的图像融合方法；这些方法都取得了不错的融合效果。然而，上述方法在对太阳斑点图像进行融合的过程中都存在不同程度的误差，这种误差的存在会导致太阳斑点图像中部分有用信息的丢失。

本文提出一种基于特征信息检测的剪切波变换域图像融合算法。该算法先对所有待融合的太阳斑点图像进行剪切波变换，得到各自的变换系数；然后通过一种简单的变换域图像融合方法获得一幅初始的融合图像，并根据每一幅太阳斑点图像与上述初始融合图像局部像素灰度值之间的相似程度来判别每一源图像像素的特征信息情况，进而得到所有太阳斑点图像的局部特征检测图，并对上述局部特征检测图进行适当处理来提升局部特征检测的准确性；再根据每一太阳斑点图像的局部特征信息检测图将该斑点图像中的所有像素分成有效特征信息区域内部、有效特征信息与无效信息区域之间的过渡以及无效信息区域内部，并据此来指导剪切波变换域各子带系数的融合；最后对处理所得融合图像的剪切波变换系数进行剪切波变换的逆变换，从而得到最终的融合图像。

本文内容安排如下：第2部分，详细描述用于局部特征检测算法，剪切波变换系数的融合规则，以及具体的融合算法步骤；第3部分，给出相关实验结果，并与其它基于多尺度变换的融合算法进行对比，验证本文算法的可行性；第4部分，对本文方法进行总结，说明优点与不足。

2. 算法描述

2.1. 算法步骤

本章提出的基于特征信息检测的剪切波变换域太阳斑点图像融合算法的具体过程可以简单地描述如下：

1) 采用相同的参数，对已配准的所有斑点图 $A_k\left(k=1,2,\cdots ,K\right)$ 分别进行J尺度剪切波分解，得到各自的低频子带系数和不同尺度分解层上各方向的带通方向子带系数。

2) 根据2.2节介绍的特征信息检测算法对各斑点图 $A\_k\left(k=1,2,\cdots ,K\right)$ 进行有效特征信息区域检测。

3) 根据有效特征信息区域检测结果，结合局部特征信息增强，采用2.3节介绍的融合规则，分别对低频子带系数及各带通方向子带系数进行融合处理，得到融合图像的剪切波变换系数。

4) 对步骤(3)所得到的融合图像的剪切波系数进行剪切波逆变换，得到最终的融合图像。

本文算法流程图1如下：

Figure 1. Algorithm flow chart

图1. 算法流程图

2.2. 特征信息检测

假设有K (K > 2)帧太阳斑点图像(本文以5帧为例，如图2)，对其进行融合的目的是要通过这K帧部分清晰的太阳斑点图像获得一帧所有区域都清晰的图像。为了有效检测出太阳斑点图像中的有效特征信息区域和无效信息区域，需要先获得一帧所有区域都比太阳斑点图像中无效信息清晰的初始融合图像。本文所用的获取初始融合图像的方法如下；

1) 对已配准的太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 进行J尺度剪切波变换分解，得到各自的低频子带系数和带通方向子带系数。

2) 对低频子带采用“取平均”、带通方向子带采用“模值取大”的融合规则分别进行融合处理，得到融合图像的剪切波系数。

3) 对步骤(2)所得的融合图像的剪切波系数，进行J尺度剪切波逆变换，得到初始融合图像F。

Figure 2. Solar image

图2. 太阳斑点图

为了检测出太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 中的有效特征信息区域和无效信息区域，在获取上述初始融合图像F以后，还需要逐像素考虑每一幅太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 与上述初始融合图像F之间的相似性。对任意像素位置 $\left(i,j\right)$，如果第k幅斑点图 $A_k$ 与初始融合图像F在像素位置 $\left(i,j\right)$ 处局部相似性在所有K幅斑点图与初始融合图像F在该位置处的局部相似性中最大，则认为斑点图 $A_k$ 中像素位置 $\left(i,j\right)$ 处于其有效特征信息区域，否则认为像素位置 $\left(i,j\right)$ 处于斑点图 $A_k$ 的无效信息区域。局部平方误差可以在很大程度上反映两幅图像之间的局部相似程度，局部平方误差越小，说明这两幅图像在某局部范围内的相似性越离，因而可以用来度量太阳斑点图像与上述初始融合图像局部像素灰度值之间的相似性，进而据此判断出所有太阳斑点图像各像素的有效特征信息。有效特征信息区域检测的具体过程如下：

1) 分别计算初始融合图像F与太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 在每一像素领域内的局部平方误差 $S{E}_k\left(i,j\right)\left(k=1,2,\cdots ,K\right)$，

$S{E}_k\left(i,j\right)={\displaystyle \underset{m=-M}{\overset{m=M}{\sum }}{\displaystyle \underset{n=-N}{\overset{n=N}{\sum }}{\left(I_F\left(i+m,j+n\right)-I_k\left(i+m,j+n\right)\right)}^2}}$ (1)

其中 $I_k\left(i,j\right)\left(k=1,2,\cdots ,K\right)$ 和 $I_F\left(i,j\right)$ 分别为太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 和初始融合图像F在像素位置 $\left(i,j\right)$ 处的灰度值，M、N分别为邻域窗口的水平与垂直方向的半径，可以取0，1，2，3等。

2) 对每一像素位置 $\left(i,j\right)$，比较 $S{E}_k\left(i,j\right)\left(k=1,2,\cdots ,K\right)$，以确定像素 $\left(i,j\right)$ 在哪副太阳斑点图像中处于有效特征信息区域。从而对 $k=1,2,\cdots ,K$，可得如下特征信息检测图：

$Z\left(i,j\right)=\{\begin{array}{l}1,ifk=\arg \min S{E}_q\left(i,j\right)\\ 1\le q\le k.\left(k=1,2,\cdots ,K\right)\\ 0,otherwise\end{array}$ (2)

当太阳斑点图像 $A_k$ 中像素 $\left(i,j\right)$ 与初始融合图像F的对应像素之间的局部平方误差最小时，说明斑点图 $A_k$ 中的像素 $\left(i,j\right)$ 与初始融合图像F的对应像素之间具有最大的相似性，从而有理由认为：斑点图 $A_k$ 在像素 $\left(i,j\right)$ 处是位于有效特征信息区域的，否则就认为斑点图 $A_k$ 在像素 $\left(i,j\right)$ 处是位于无效信息区域的。在决策矩阵 $Z_k$ 中， $Z_k\left(i,j\right)=1$ 表示太阳斑点图 $A_k$ 的像素 $\left(i,j\right)$ 是位于有效特征信息区域的， $Z_k\left(i,j\right)=0$ 表示太阳斑点图 $A_k$ 的像素 $\left(i,j\right)$ 是位于无效信息区域的，获得的特征信息检测图如图3。

Figure 3. Feature information detection diagram

图3. 特征信息检测图

3) 有效特征信息区域边界像素的检测

用大小为 $M_1\times N_1$ 、元素全为1的滑动窗口与特征信息检测图 $Z_k\left(k=1,2,\cdots ,K\right)$ 进行卷积，重新定义 $Z_k\left(i,j\right)$ 为

$Z_k\left(i,j\right)={\displaystyle \underset{x=-\left(M_1-1\right)/2}{\overset{\left(M_1-1\right)/2}{\sum }}{\displaystyle \underset{y=-\left(N_1-1\right)/2}{\overset{\left(N_1-1\right)/2}{\sum }}{Z}_k\left(i+xj+x\right)}}$ (3)

则 $Z_k\left(i,j\right)=M_1{N}_1$ 表示像素 $\left(i,j\right)$ 位于斑点图 $A_k$ 有效特征信息区域内部， $0<Z_{k|k=1,2,\cdots ,K}\left(i,j\right)<M_1{N}_1$ 表示像素 $\left(i,j\right)$ 位于斑点图 $A_k$ 有效特征信息区域边界上， $Z_k\left(i,j\right)=0$ 表示像素 $\left(i,j\right)$ 位于斑点图 $A_k$ 无效信息区域内。

2.3. 融合规则

1) 低频子带系数的融合策略

图像经剪切波分解所得到的低频子带代表了原图像的近似信息，是图像轮廓的反映，低频子带系数融合的准确与否直接影响着融合图像的质量。本文使用2.2节得到的决策矩阵作为融合依据。

记 $I_L^{A_k}\left(i,j\right)\left(k=1,2,\cdots ,K\right)$ 和 $I_L^F\left(i,j\right)$ 分别为太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 和融合图像F在位置 $\left(i,j\right)$ 处的低频子带系数。由于 $Z_k\left(i,j\right)=M_1{N}_1\left(k=\arg {\max }_{1\le q\le K}{Z}_q\left(i,j\right)\right)$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 有效特征信息区域内部； $0<Z_k\left(i,j\right)<M_1{N}_1$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 有效特征信息区域与无效信息区域的边界上； $Z_k\left(i,j\right)=0$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 无效信息区域内。因此可直接将图像 $A_k$ 位于 $\left(i,j\right)$ 处的低频子带系数作为融合图像在对应位置的低频子带系数，从而可以在很大程度上避免传统基于像素以及基于区域的融合规则对低频子带系数的误选。据此，本文提出如下基于特征信息检测的低频子带系数融合规则：

$I_L^F\left(i,j\right)=I_L^{A_k}\left(i,j\right),k=\arg \max Z_q\left(i,j\right)$ (4)

2) 带通方向子带系数的融合策略

在多尺度变换域内，图像中物体的边缘和轮廓等显著特征通常都表现为模值较大的高频子带系数。本文提出一种基于特征信息检测和局部增强的带通方向子带系数融合策略。

记 $I_P^{A_k}\left(i,j\right)\left(k=1,2,\cdots ,K\right)$ 和 $I_P^F\left(i,j\right)$ 分别为太阳斑点图像 $A_k\left(k=1,2,\cdots ,K\right)$ 和融合图像F在位置 $\left(i,j\right)$ 处的带通方向子带系数。由于 $Z_k\left(i,j\right)=M_1{N}_1\left(k=\arg {\max }_{1\le q\le K}{Z}_q\left(i,j\right)\right)$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 有效特征信息区域内部； $0<Z_k\left(i,j\right)<M_1{N}_1$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 有效特征信息区域与无效信息区域的边界上； $Z_k\left(i,j\right)=0$ 表明像素 $\left(i,j\right)$ 位于太阳斑点图像 $A_k$ 无效信息区域内。因此，设置一个阈值W，使其满足 $0<W<M_1{N}_1$，当 $Z_k\left(i,j\right)>W$ 时，将图像 $A_k$ 位于 $\left(i,j\right)$ 处的带通方向子带系数叠加到融合图像F上，从而对像素 $\left(i,j\right)$ 处的有效特征信息进行增强。据此，本文采用基于特征信息检测和局部增强的带通方向子带系数融合规则，即：

$I_P^F\left(i,j\right)=I_P^{A_k}\left(i,j\right)+I_P^F\left(i,j\right),Z_k\left(i,j\right)>W$ (5)

3. 实验与结果分析

为了验证本文算法的有效性，选取三组太阳斑点图像进行融合实验，三组图像分别含有5幅大小为500 × 500的斑点图；并与基于DWT的简单融合算法、基于NSCT的简单融合算法、基于剪切波变换的简单融合算法进行对比融合实验。在对太阳斑点图像经多尺度分解后的各子带系数进行融合时，基于DWT与基于NSCT的简单融合算法均采用“低频系数取平均，高频系数模值取大”的融合规则。基于DWT的简单图像融合算法采用“db4”小波对斑点图进行４尺度小波分解；基于NSCT的简单图像融合算法对源图像进行４尺度NSCT分解，各尺度的分解方向数依次为2，2³，2³，2³；基于剪切波变换的简单融合算法与本文算法均对斑点图进行4尺度剪切波变换分解时，各尺度的分解方向数均取2² + 2。图4、图6与图8为用于实验的三组太阳斑点图像，图5、图7与图9为各算法对三组太阳斑点图像进行融合的结果图像。

对于图像融合的效果，除了目视的主观评价之外，本文采用灰度标准差(SD)、空间频率(SF)以及图像信息熵(IE)来对其进行客观定量地评价。表1、表2、表3分别给出了不同融合算法对三组太阳斑点图像融合结果的客观评价指标。

Figure 4. First set of images

图4. 第一组图像

Figure 5. Fusion results of the first set of images by different methods

图5. 不同方法对第一组图像的融合结果

Table 1. Evaluation of the first set of experimental results

表1. 第一组实验结果的评价

Figure 6. Second set of images

图6. 第二组图像

Figure 7. Fusion results of the second set of images by different methods

图7. 不同方法对第二组图像的融合结果

Table 2. Evaluation of the second set of experimental results

表2. 第二组实验结果的评价

Figure 8. Third set of images

图8. 第三组图像

Figure 9. Fusion results of the third set of images by different methods

图9. 不同方法对第三组图像的融合结果

Table 3. Evaluation of the third set of experimental results

表3. 第三组实验结果的评价

通过上述实验结果可以看出，本文方法得到的结果图像相对其它方法得到的结果图像从视觉上有较为明显的提升，更好的融合了各斑点图的局部特征，使其更加直观清晰。综合表中数据可以发现，本文方法具有比基于DWT、NSCT以及剪切波融合算法更好的融合性能，在灰度标准差、空间频率以及图像信息熵上明显优于其余三种算法。

4. 总结

本文针对太阳斑点图像融合问题，提出了一种基于剪切波变换和特征信息检测的太阳斑点图融合算法。实验表明：该融合算法能在融合图像中融入更多各斑点图中的有用信息，有效避免“伪影”等虚假信息的引入，从而有效提高融合图像的质量，是一种耗时少、融合性能好的图像融合算法。同时本文方法的结果图像中仍有一定程度的噪声，影响了重建结果的整体效果，需要在抑制噪声方面进一步改进。

致谢

非常感谢蒋慕容老师，在研究课题以及完成论文的过程中给自己的悉心指导，同时感谢各位同学在此期间给我的帮助和关心，在此向他们表达衷心的感激。