在研究太阳图像的观测数据时，经常需要将不同设备的数据进行对比，从而发挥它们各自在不同波段、不同视场、不同分辨率上的优势。太阳光学成像观测包括全日面和高分辨观测两类。对于全日面图像，由于太阳几何尺寸和运动方式已知，如太阳的半径、太阳极轴角、自转等，因此各种观测设备均把自己的数据在预处理后归算到标准的太阳坐标系。然而对于局部高分辨像，由于观测设备的指向误差、视场旋转以及焦距的变化，使得直接依赖仪器参数确定准确的太阳视场坐标不精确，因此只能借助太阳图像上相似的太阳结构(如光球层对应光球层，色球层对应色球层)与全日面像进行转换，即视场匹配。然而全日面观测有大于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太阳望远镜像视场相同的太阳动力天文台的日震及磁场成像仪子图I_S；

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

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

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

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

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

(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.2 分块确定特征控制点

设I_R为待匹配图像，I_S为参考图像。将I_R划分N个特征区域，记为集合{f_Rⁱ, i=1, 2, ..., N}；I_S的尺寸为U × V，每个区域f_Rⁱ的尺寸为u × v，可知u＜U, v＜V，则参考图像I_S和图像f_Rⁱ的互相关函数Corr(x, y)定义为

$ Corr\left( {x, y} \right) = \frac{{\sum\limits_{\mathit{s}{\rm{ = 0}}}^{\mathit{u}{\rm{ - 1}}} {\sum\limits_{\mathit{t}{\rm{ = 0}}}^{\mathit{v}{\rm{ - 1}}} {\left[{{\mathit{I}_{\rm{S}}}\left( {\mathit{x}{\rm{ + }}\mathit{s}{\rm{, }}\;\mathit{y}{\rm{ + }}\mathit{t}} \right)-\overline {\mathit{f}_{\rm{S}}^\mathit{i}} } \right]\left[{\mathit{f}_{\rm{R}}^\mathit{i}\left( {\mathit{s}{\rm{, }}\;\mathit{t}} \right)-\overline {\mathit{f}_{\rm{R}}^\mathit{i}} } \right]} } }}{{\sqrt {\sum\limits_{\mathit{s}{\rm{ = 0}}}^{\mathit{u}{\rm{ - 1}}} {\sum\limits_{\mathit{t}{\rm{ = 0}}}^{\mathit{v}{\rm{ - 1}}} {{{\left[{{\mathit{I}_{\rm{S}}}\left( {\mathit{x}{\rm{ + }}\mathit{s}{\rm{, }}\;\mathit{y}{\rm{ + }}\mathit{t}} \right)-\overline {\mathit{f}_{\rm{S}}^\mathit{i}} } \right]}^2}\sum\limits_{\mathit{s}{\rm{ = 0}}}^{\mathit{u}{\rm{ - 1}}} {\sum\limits_{\mathit{t}{\rm{ = 0}}}^{\mathit{v}{\rm{ - 1}}} {{{\left[{\mathit{f}_{\rm{R}}^\mathit{i}\left( {\mathit{s}{\rm{, }}\;\mathit{t}} \right)-\overline {\mathit{f}_{\rm{R}}^\mathit{i}} } \right]}^2}} } } } } }}, $

(1)

其中，(s, t)为f_Rⁱ上的坐标；$ \overline {\mathit{f}_{\rm{S}}^\mathit{i}} {\rm{ = }}\frac{1}{{\mathit{uv}}}\sum\limits_{\mathit{s}{\rm{ = 0}}}^{\mathit{u}{\rm{ - 1}}} {\sum\limits_{\mathit{t}{\rm{ = 0}}}^{\mathit{v}{\rm{ - 1}}} {\left[{{\mathit{I}_{\rm{S}}}\left( {\mathit{x}{\rm{ + }}\mathit{s}{\rm{, }}\;\mathit{y}{\rm{ + }}\mathit{t}} \right)} \right]} } $；$ \overline {\mathit{f}_{\rm{R}}^\mathit{i}} {\rm{ = }}\frac{1}{{\mathit{uv}}}\sum\limits_{\mathit{s}{\rm{ = 0}}}^{\mathit{u}{\rm{ - 1}}} {\sum\limits_{\mathit{t}{\rm{ = 0}}}^{\mathit{v}{\rm{ - 1}}} {\left[{\mathit{f}_{\rm{R}}^\mathit{i}\left( {\mathit{s}{\rm{, }}\;\mathit{t}} \right)} \right]} } $。$ \overline {\mathit{f}_{\rm{S}}^\mathit{i}} $为I_S和f_Rⁱ覆盖的那个区域图像的平均值，求和是对I_S和f_Rⁱ覆盖的(s, t)进行。

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

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

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

$ {F_i} = F(\mathit{f}_{\rm{R}}^\mathit{i}{\rm{, }}\mathit{f}_{\rm{S}}^\mathit{i}), \;1 \le i \le N, $

(2)

F_i表示任意两个对应的局部特征区域(f_Rⁱ, f_Sⁱ)。

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

$ m = \frac{{\sum\nolimits_\mathit{D} {\left( {I' X} \right)} }}{{\sum\nolimits_\mathit{D} {I' } }}, \;\;n = \frac{{{\sum _D}\left( {I' Y} \right)}}{{{\sum _D}I' }}, $

(3)

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

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

$ {P_i} = P(\mathit{p}_{\rm{R}}^\mathit{i}{\rm{, }}\;\mathit{p}_{\rm{S}}^\mathit{i}), \;\;1 \le i \le N{\rm{, }} $

(4)

2.3 最小二乘求解转换参数

考虑平移、旋转和缩放的情况，假设控制点对(p_Rⁱ, p_Sⁱ)，对应的坐标分别为(x_A, y_A)和(x_B, y_B)，变换关系为

$ \begin{array}{l} {x_B} = \Delta x + \left( {1 + m} \right){x_A}{\rm{cos}}\theta - \left( {1 + m} \right){y_A}{\rm{sin}}\theta {\rm{, }}\\ {y_B} = \Delta y + \left( {1 + m} \right){x_A}{\rm{sin}}\theta + \left( {1 + m} \right){y_A}{\rm{cos}}\theta . \end{array} $

(5)

转换为矩阵形式，

$ {\left[\begin{array}{l} x\\ y \end{array} \right]_B} = \left[{\begin{array}{*{20}{c}} 1&0&{{x_A}}&{-{y_A}}\\ 0&1&{{y_A}}&{{x_A}} \end{array}} \right]\left[\begin{array}{l} \Delta x\\ \Delta y\\ \;a\\ \;b \end{array} \right], \;\begin{array}{*{20}{c}} {a = \left( {1 + m} \right){\rm{cos}}\theta }\\ {b = \left( {1 + m} \right){\rm{sin}}\theta } \end{array}. $

(6)

误差方程为

$ \left[\begin{array}{l} {V_{{x_A}}}\\ {V_{{y_B}}} \end{array} \right] = \left[{\begin{array}{*{20}{c}} 1&0&{{x_A}}&{-{y_A}}\\ 0&1&{{y_A}}&{{x_A}} \end{array}} \right]\left[\begin{array}{l} \Delta x\\ \Delta y\\ \;a\\ \;b \end{array} \right] - \left[\begin{array}{l} {x_B}-{x_A}\\ {y_B}-{y_A} \end{array} \right]{\rm{, }} $

(7)

其中，$ {\left[\begin{array}{l} \mathit{x}\\ \mathit{y} \end{array} \right]_\mathit{A}}$为转换前坐标；${\left[\begin{array}{l} \mathit{x}\\ \mathit{y} \end{array} \right]_\mathit{B}}$为转换后坐标；$ \left[\begin{array}{l} \Delta \mathit{x}\\ \Delta \mathit{y} \end{array} \right]$为平移参数；(1 + m)为缩放参数；θ为旋转参数。理论上，只要存在两组转换控制点，利用上述方程，可联立解出一组转换四参数；存在3个或3个以上的控制点时，可利用最小二乘法求解四参数X=(Δx, Δy, m, θ)^T；然后将待匹配图像转换到参考图像。

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

3 实验结果与分析

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

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

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

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

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

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

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

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

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

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

4 结论

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

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

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

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