太阳高分辨观测图像与全日面像的高精度配准方法*

冯 涛，戴 伟 ，王 瑞，季凯帆,

(1. 昆明理工大学云南省计算机技术应用重点实验室，云南 昆明 650500；2. 中国科学院云南天文台，云南 昆明 650011)

在研究太阳图像的观测数据时，经常需要将不同设备的数据进行对比，从而发挥它们各自在不同波段、不同视场、不同分辨率上的优势。太阳光学成像观测包括全日面和高分辨观测两类。对于全日面图像，由于太阳几何尺寸和运动方式已知，如太阳的半径、太阳极轴角、自转等，因此各种观测设备均把自己的数据在预处理后归算到标准的太阳坐标系。然而对于局部高分辨像，由于观测设备的指向误差、视场旋转以及焦距的变化，使得直接依赖仪器参数确定准确的太阳视场坐标不精确，因此只能借助太阳图像上相似的太阳结构(如光球层对应光球层，色球层对应色球层)与全日面像进行转换，即视场匹配。然而全日面观测有大于31′ 的视场，但最高的空间像元分辨率仅为0.5″(空间太阳动力天文台的日震及磁场成像仪，Solar Dynamics Observatory/Helioseismic and Magnetic Imager, SDO/HMI)。而地面高分辨观测，如中国的1 m新真空太阳望远镜(New Vaccum Sloar Telescope, NVST)和美国大熊湖天文台的新太阳望远镜(New Solar Telescope, NST)，视场仅仅为1′ 到2′，但像元分辨率达到0.05″ 以上，这就意味着需要将两个空间分辨率相差10倍的图像进行配准。即便不考虑其他畸变，局部像的视场旋转(R)、比例尺(S)和中心位置坐标(Tx,Ty)是必须得到的4个参数。

图像配准是指在不同时间、不同条件下得到的位于不同坐标系下的同一场景的两幅或者多幅图像进行匹配的过程[1]。由于存在各种不同的情况，同一场景的多幅图像会在分辨率、位置、尺度、成像模式等方面存在很多差异，图像配准是找到一个几何变换的最优解，使得图像上相对应的点，或者至少是有意义的点达到匹配状态。目前在图像配准领域，常见的配准方法有基于特征的点匹配方法，即根据图像中的特征点，或者显著特征构造描述子确定的匹配点，然后构造转换方程进行配准。理论上这个方法如果有足够多和足够精确的控制点是可以很好地解决两个图像的几何变换，如旋转、缩放、平移等问题。但关键问题是如何确定这些点，这方面在图像处理界有很多研究，如Harris角点及其衍生的Harris-Laplace方法[2]，尺度不变特征(Scale-invariant feature transform, SIFT)[3]及其衍生的PCA-SIFT算法(Principal Com-ponent Analysis SIFT)[4]、局部特征描述子[5]、高速尺度不变特征算法(Very Fast SIFT, VF-SIFT)[6]，加速鲁棒性特征(Speeded Up Robust Features, SURF)[7-8]等。另外一种是基于区域信息的统计配准方法，主要通过图像区域的统计信息度量图像的相似程度，如互信息[9]，互相关函数[10-11]，模板匹配[12]，不变矩，相关匹配[13-14]等。但一般来讲这种方法主要解决图像间的平移问题，对于具有旋转和比例缩放的图像，又引申出在频率域的相位相关技术[15]，Fouier-Mellin[16-17]等。

然而，由于太阳观测图像和普通自然图像有很大的区别，其中的特征边界不清晰，结构相对模糊且存在很多噪声，同时动态范围大，运动速度较快，并且存在很多不同的活动现象，大部分运动现象为非刚性运动；同时太阳图像的图幅一般较大，有效目标小且缺少有效的纹理，不存在明显的梯度变化等特征，并且由于拍摄太阳图像的设备以及拍摄过程中的不确定性，采集的太阳图像还存在旋转、缩放及平移等情况，不同图像的空间分辨率经常不在同一尺度，这就为确定特征对应点带来非常大的困难。长期以来，同一观测设备的太阳时间序列图像的视场对齐一直采用基于互相关的统计配准，而不同观测设备之间的图像配准，还停留在利用一些显著的太阳特征，如黑子、气孔、大米粒构造对应点的方式。然而这些特征点数量少，结构不规则，难以保证位置测量精度，同时也难以使用计算机技术进行自动化测量。尤其是将高分辨局部像对应到全日面像的时候，这种问题更为突出。但从另一方面，通过望远镜机械和终端光学系统，对局部太阳像指向、比例尺和位置有初步估计，只是这些数值不很精确且随时间有一些变化。因此实际面临的是如何提高配准精度，进行精细图像匹配的问题。

基于这些问题， 本文提出了一种将统计信息和相关点匹配相结合的太阳图像配准算法，解决将不同分辨率下的局部高分辨观测数据和全日面数据进行高精度坐标转换的问题，其中包括测量图像指向角、像元比例尺和图像中心在全日面像上的位置。这一方法的基本思想是在初步匹配的基础上，首先将视场划分成大量重叠的局部区域代替点匹配中的控制点，通过基于统计的相关匹配确定对应区域间的亚像元位移，从而确定特征点对在同一坐标系中的位置；然后用最小二乘迭代求解整个视场的转换参数，并逐步迭代提高拟合的精度，成功将1 m太阳望远镜和新太阳望远镜的TiO观测数据与几乎同时刻(最长间隔11 s)的太阳动力天文台的日震及磁场成像仪连续谱进行了匹配，拟合结果的偏差在0.25″ 之内。

1 实验数据来源

本文实验使用的全日面图像来自太阳动力天文台(Solar Dynamics Observatory, SDO)上日震及磁场成像仪(Helioseismic and Magnetic Imager, HMI)采集的连续谱全日面像，而高分辨图像为1 m太阳望远镜和新太阳望远镜 的TiO光球像。SDO/HMI设备的空间分辨率约为0.504 3 arcsec/pixel，图像尺寸为4 096像元 × 4 096像元；1 m太阳望远镜和新太阳望远镜图像的像元比例尺为0.04 arcsec/pixel和0.03 arcsec/pixel。

图1(a)为1 m太阳望远镜的TiO观测数据，观测时间为2013-07-15 07:29:09，有效区域尺寸为2 304像元 × 1 920像元；图1(b)为SDO/HMI连续谱观测数据，观测时间为2013-07-15 07:29:15， 红框标记区域为1 m太阳望远镜在SDO/HMI中的大体位置。本文主要以此图像的匹配说明算法流程。

2 配准的方法

2.1 配准方法的流程

配准方法的主要流程如图2。图2可以描述为下列11个步骤：

(1)对1 m太阳望远镜图像的指向、比例尺和位置进行初步估计，得到旋转(Rotation)、缩放(Scaling)及平移(Translation)参数(下文称RST参数)；

(2)根据RST参数，对太阳动力天文台的日震及磁场成像仪全日面像进行旋转，并截取与1 m太阳望远镜像视场相同的太阳动力天文台的日震及磁场成像仪子图IS；

图1 观测数据
Fig.1 Observed data

(3)根据RST参数和太阳动力天文台的日震及磁场成像仪的比例尺，将1 m太阳望远镜像缩小至和太阳动力天文台的日球层磁场观测仪同一比例尺下，即对1 m太阳望远镜图像进行降采样；

(4)对降采样后的1 m太阳望远镜像进行高斯平滑，减小采样过程中的频率混叠。平滑后的1 m太阳望远镜降采样像用IR表示；

(5)将图像IR等间隔划分为N个重叠的子块，并将每个子块记为fiR；

(6)对每一个子块fiR和太阳动力天文台的日球层磁场观测仪子图IS利用归一化互相关方法求得fiR对应的每一个IS子块fiS；

(7)对每一组fiR和fiS，测量其亚像元偏移，并将子块中心设置为一组对应的控制点；

(8)筛选每组控制点，将控制点间距离大于一定像元的点剔除；

(9)根据剩下的控制点，建立仿射变换的转换方程，采用最小二乘求解RST参数；

(10)计算残差σ和迭代次数L，若满足{σ<0.5∧Δσ<0.1}；{σ<0.3}；{L>30,σ=min(σ1～30)}任一条件，则停止迭代；否则，返回步骤(2)，根据步骤(9)解得的RST参数重新对SDO/HMI和1 m太阳望远镜像进行迭代处理；

(11)最终得到迭代完成后的RST参数并进行配准。

在上述步骤中， (6)至(9)步为算法的核心，下文对该部分进行更为具体的描述。需要注意，开始迭代后(2)到(3)步中的RST参数会根据上一次迭代的结果进行下一次迭代，意味着算法会重新对图像进行平滑。

图2 算法流程图
Fig.2 Algorithm flow chart

2.2 分块确定特征控制点

设IR为待匹配图像，IS为参考图像。将IR划分N个特征区域，记为集合{fiR,i=1, 2, ...,N}；IS的尺寸为U×V，每个区域fiR的尺寸为u×v，可知u

(1)

根据上述公式，每个fiR和IS做归一化互相关，可得到一个互相关曲面，而峰值的坐标代表了fiR图像中心在IS中的整数像元精度的匹配位置。

在IS中根据上述位置可截取和fiR同样大小的特征区域，得到每一个fiR对应的fiS。由它们组成的集合记为fiS,i=1, 2, ...,N。

对于两组对应的特征区域的fiR和fiS，使用集合F描述fiR和fiS任意两个特征区域之间的关系，即

Fi=F(fiR,fiS),1≤i≤N，

(2)

Fi表示任意两个对应的局部特征区域(fiR,fiS)。

对于每一对特征区域(fiR,fiS)，计算标准的互相关函数。假设D为特征区域(fiR,fiS)间互相关函数分布的峰值局部区域，按照修正矩方法计算重心坐标(m,n)：

(3)

其中，I为曲面强度值，则I′=I-b，b=min(I)，I∈D。每一对(fiR,fiS)图像间互相关函数的局部重心坐标代表每一对特征区域(fiR,fiS)间的亚像元偏移。

根据亚像元偏移的情况， 可以建立每一对(fiR,fiS)图像中心点对应的坐标，得到和IS之间的N对坐标转换控制点。用集合P描述控制点对之间的关系，即

Pi=P(piR,piS),1≤i≤N.

(4)

2.3 最小二乘求解转换参数

考虑平移、旋转和缩放的情况，假设控制点对(piR,piS)，对应的坐标分别为(xA,yA)和(xB,yB)，变换关系为

xB=Δx+(1+m)xAcosθ-(1+m)yAsinθ，

yB=Δy+(1+m)xAsinθ+(1+m)yAcosθ.

(5)

转换为矩阵形式，

(6)

误差方程为

(7)

需要注意，在视场平面坐标转换时，需要对误配点对进行筛选，通过计算控制点间的距离。在实际匹配中，如果控制点之间的距离大于5像元，即2.5″ 时，则认为该点为误配点对，将其剔除，根据筛选后的控制点求解的RST参数，重新对太阳动力天文台的日震及磁场成像仪和1 m太阳望远镜像进行RST迭代，计算RST参数，最后得到迭代完成后的RST转换参数。

3 实验结果与分析

本文实验中，分块尺寸设置为30像元 × 30像元，块间隔为5像元。图3为利用上述方法，将1 m太阳望远镜高分辨像和太阳动力天文台的日震及磁场成像仪全日面像配准后重叠区域的对比图。

图3 1 m太阳望远镜高分辨像和太阳动力天文台的日球层磁场观测仪全日面像配准后重叠区域的对比图
Fig.3 Comparison overlapping area of NVST high-resolution image and full-disk Solar image after registration

从图3可以看到，配准后图像的对应区域除清晰度明显不同外，其他特征已经完全匹配。为了更好地观察配准效果，将配准区域叠加后，对边界区域进行放大，如图4，可以看到边界细节的信息匹配度良好，气孔、米粒、米粒暗径均较好地衔接。

根据配准结果和太阳动力天文台的日震及磁场成像仪的标准参数，表1显示了推算的这张1 m太阳望远镜图像极轴角和像元比例尺，以及视场中心在太阳动力天文台的日震及磁场成像仪全日面像中的位置。

表1 配准结果Table 1 Registration results

从表1可以知道平均残差为0.31太阳动力天文台的日震及磁场成像仪像元，意味着拟合结果的偏差约为0.15″。图5显示了这些残差的空间分布。从图5可以看出，残差分布并不是完全随机，显示了1 m太阳望远镜像存在畸变。这是由于地面高分辨观测图像重建后，依然还有一些低频的大气抖动残余未能消除，而这限制了最终拟合结果的精度。同时由于两幅图像的观测时间并不完全同步，有6 s的间隔，因此有些太阳图像特征已经发生了改变。

图4 配准区域的边界细节
Fig.4 The edge details of the registration area

对多张1 m太阳望远镜和新太阳望远镜观测的高分辨像和太阳动力天文台的日震及磁场成像仪全日面像进行了匹配，结果见表2。

从上述结果可以看到，最大的匹配残差约0.5像元，相当于0.25″。

实验过程中发现，不同参数的选取会对算法产生一定的影响，具体见图6～图8。

图6显示了初步估计的参数对迭代次数的影响。图6中，三角型(绿色)、空心圆型(蓝色)、米字型(红色)、十字型(黑色)分别代表估计角度为30°、35°、40°、45°时的迭代过程。可以看出，随着初值误差增大，为了修正误差使结果收敛在合理的范围内，算法的迭代次数相应增加。当初始估计角度达到45°(此时的误差约为30°)时，迭代次数由6次增加到42次才能收敛。

图5 配准后的向量场
Fig.5 The vector field after registration

表2 新太阳望远镜和1 m太阳望远镜的TiO图像与SDO/HMI全日面像的配准结果Table 2 Registration result of TiO image (NST and NVST) and SDO/HMI full-disk Solar image

图6 不同的初始参数
Fig.6 The different initial parameters

图7 不同的控制点筛选间距
Fig.7 The different filter distance between control points

图7显示了筛选控制点间距离的设置对算法收敛速度的影响。图7中，十字型(红色)、空心圆型(黑色)分别代表筛选阈值为5 pixel和10 pixel的迭代过程。可以看出，初值误差较大时，高阈值的迭代次数更少，意味着算法收敛更快，反之低阈值效果更理想。当角度误差在30°附近时，不同阈值的迭代次数相差28次，高阈值更优。角度误差在10°附近时，迭代次数仅相差9次，低阈值更优。

图8显示了不同的平滑参数对迭代次数的影响。图8中，米字型(绿色)、空心圆型(蓝色)、十字型(红色)分别表示平滑参数为0.5、1、1.5时的迭代过程。可以看出，平滑参数增大使迭代次数增加，并且拟合结果有一定的降低。当平滑参数为0.5时，算法达到收敛所需的迭代次数和结果优于其它情况。

图8 不同的平滑参数
Fig.8 The different filtering parameters

4 结 论

本文提出了一种针对太阳高分辨像和全日面像进行视场匹配的算法，从而解决了不同空间分辨率下、具有旋转角的太阳观测图像匹配问题。算法的基本思想是结合区域相关匹配和控制点最小二乘求解的方法。1 m太阳望远镜和新太阳望远镜高分辨图像与SDO/HMI全日面像的视场匹配实验都证明了两者图像拟合结果的偏差优于0.25″。通过匹配可以精确测量高分辨像的极轴角、像元比例尺等。

本文方法要计算大量的小区域相关，同时进行迭代计算，因此，时间开销较大，应用时需要注意不同参数对迭代次数的影响。除此之外，小尺寸子块会导致匹配点数增加，计算量增大，拟合精度会存在一定程度的下降；而大尺寸虽然对精度有一定提升，但是会导致匹配区域变小，从而影响整张图像的匹配效果，子块间距更多地影响计算速度。结合实际情况，我们认为子块尺寸在30像元、块间距5像元附近，算法的稳定性较合理。在应用中，可根据实际情况对参数进行调整。

此外，将高分辨TiO像通过与SDO/HMI全日面连续谱图像的高精度匹配，也意味着可以转换为和任意标准的全日面像进行匹配，如SDO/HMI磁图、SDO/AIA各波段观测，甚至和色球全日面像进行匹配。这意味着可以在0.25″ 的偏差下，研究小尺度结构在各个波段，以及在磁场上的对应情况。

这套方案可以用于任何具有相似图像结构的太阳图像匹配，如TiO和G波段、Hα偏带、10830，也可以用于Hα线心图像和全球日震网全日面像的配准，因此具有非常广阔的应用前景。

[1] Zitova B, Flusser J. Image registration methods: a survey[J]. Image and Vision Computing, 2003, 21(11): 977-1000.

[2] Mikolajczyk K, Schmid C. Scale & affine invariant interest point detectors[J]. International Journal of Computer Vision, 2004, 60(1): 63-86.

[3] Lowe D G. Distinctive image features from scale-invariant keypoints[J]. International Journal of Computer Vision, 2004, 60(60): 91-110.

[4] Yan K, Sukthankar R. PCA-SIFT: a more distinctive representation for local image descriptors[C]// IEEE Computer Society Conference on Computer Vision & Pattern Recognition. 2004: 506-513.

[5] Hou J, Qi N M, Kang J X. Image matching based on representative local descriptors[C]// International Conference on Advances in Multimedia Modeling. 2010: 303-313.

[6] Alhwarin F, Ristic′-Durrant D , Gräser A . VF-SIFT: Very Fast SIFT feature matching[C]// Dagm Conference on Pattern Recognition. 2010: 222-231.

[7] Bay H, Tuytelaars T, Gool L V. SURF: speeded up robust features[J]. Computer Vision & Image Understanding, 2006, 110(3): 404-417.

[8] 张开玉, 梁风梅. 基于改进SURF的图像配准关键算法研究[J]. 科学技术与工程, 2013, 13(10): 2875-2879.

Zhang Kaiyu, Liang Fengmei. Research on key algorithms of image registration based on improved SURF[J]. Science Technology and Engineering, 2013, 13(10): 2875-2879.

[9] Pluim J P W, Maintz J B A, Viergever M A. Multual-information-based registration of medical images: a survey[J]. IEEE Transactions on Medical Imaging, 2003, 22(8): 986-1004.

[10]Lewis J P. Fast normalized crossed-correlation[J]. Circuits Systems & Signal Processing, 2001, 82(2): 144-156.

[11]Lee S U, Heo Y S, Lee K M. Robust stereo matching using adaptive normalized cross-correlation[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2011, 33(4): 807-822.

[12] 高军, 李学伟, 张建, 等. 基于模板匹配的图像配准算法[J]. 西安交通大学学报, 2007(3): 307-311.

Gao Jun, Li xuewei, Zhang Jian. et al. Image registration algorithm based on template matching[J]. Journal of Xi′an Jiaotong University, 2007(3): 307-311.

[13] Sarvaiya J N, Patnaik S, Bombaywala S. Image registration using log-polar transform and phase correlation[C]// Tencon IEEE Region 10 Conference. 2009: 1-5.

[14] 辛登松, 彭建雄. SAR图像的改进相位相关配准方法[J]. 计算机与数字工程, 2011, 39(8): 126-128.

Xin Dengsong, Peng Jianxiong. Improved phase correlation registration method for SAR images[J]. Computer and Digital Engineering, 2011, 39(8): 126-128.

[15] Reddy B S, Chantterji B N. An FFT-based technique for translation, rotation, and scale-invariant image registration[J]. IEEE Transactions on Image Processing, 1996, 5(8): 1266-1271.

[16] Guo X X, Lu Y N, Xu Z W, et al. An application of Fourier-Mellin transform in image registration[C]// International Conference on Computer & Information Technology. 2005: 619-623.

[17] 王亮, 刘蓉, 张丽, 等. 基于Fourier-Mellin变换的气象卫星光谱图像配准[J]. 光谱学与光谱分析, 2013, 33(3): 855-858.

Wang Liang, Liu Rong, Zhang Li, et al. Spectral image registration of meteorological satellite based on Fourier-Mellin transform[J]. Spectroscopy and Spectral Analysis, 2013, 33(3): 855-858.