1 引言

目前,太阳多波段成像观测是实测太阳物理的重要手段,观测图像可能来自于不同仪器或者同一仪器的不同观测系统。为了综合运用多波段成像数据对太阳大气活动现象进行分析与研究,图像之间的视场匹配是一个必须解决的问题。2009年,Berkebile-Stoiser等^[1]在对太阳微耀斑的精细结构进行多波段分析时,以太阳黑子为特征用互相关算法对不同仪器的图像进行视场匹配。2010年,Yurchyshyn等^[2]在研究太阳小尺度磁浮现的色球特征时,认为在Hα波段和Ca II H波段观测的亮结构是具有空间对应性的磁亮点,从而以磁亮点为参考将美国大熊湖1.6 m太阳望远镜的Hα色球像与日本空间望远镜Hinode的Ca II H色球像进行视场匹配。2007年,Almeida等^[3]在研究太阳冕环的光球足点时也需要解决不同仪器设备所观测到的光球像、紫外像以及磁图等之间的视场匹配问题,他们同样通过寻找相似结构并利用互相关算法或试错法对不同的图像进行匹配。上述方法通常是在一定的物理模型假设下,以不同波段图像中的相似结构为特征进行视场匹配,这样得到的匹配精度较低(1″~2″),而且完全依赖图像的特征。此外,常用的互相关算法也只能解决图像之间的平移差异,并不能精确地处理由于空间分辨率(即尺度因子)差异或旋转差异造成的图像视场不匹配。因此,对于具有亚角秒空间分辨率的太阳图像而言,需要更精确的方法来对图像进行视场匹配,从而使各波段观测的太阳像达到与其空间分辨率相当的视场匹配精度。

1 m新真空太阳望远镜(NVST)作为中国新一代地面高分辨率太阳观测设备之一,于2011年正式投入观测。NVST中一个重要的终端设备是多波段高分辨率成像系统,它提供了对太阳大气从光球层到色球层的高空间分辨率(0.1″~0.3″)和高时间分辨率(约为10 s)的成像观测。该成像系统包括5个独立的通道,分别是两个光球宽带通道[TiO(705.8 nm)通道和G-band(430 nm)通道],以及三个色球窄带通道[Hα(656.28 nm)通道、He I(1083 nm)通道和Ca II H(393.3 nm)通道]^[4]。为了解决NVST多波段高分辨率成像系统各通道太阳图像视场之间的匹配问题,实现各波段图像高精度的视场匹配,提出了一套基于针孔阵列光阑的方法,即针孔光阑视场定标。本文详细叙述了该方法的原理和如何运用该方法对TiO通道和Hα通道进行视场定标,从而实现了两通道图像之间0.1″精度的视场匹配。

2 NVST多波段视场定标方案

2.1 NVST多波段成像系统光路简介

将NVST多波段成像系统放置于一个直径为6 m的可旋转光学平台上,图1(a)是NVST建筑剖面示意图^[5],红色圆圈处是放置各种终端仪器设备的光学平台所处位置;图1(b)是光学平台示意图。图2是多波段成像系统的光路示意图,其中F3是望远镜的焦平面。来自望远镜的光在F3之后经由场镜、准直镜以及分束镜等分出进入5个独立的成像通道,从成像通道出射的光分别成像在各通道中的探测器上。值得注意的是,在各个波段观测到的太阳图像存在很大差异,很难利用相似的结构特征实现视场的匹配。同时,由于焦比的差异,不同通道观测数据的空间分辨率会略有不同,而且根据观测目标不同的动力学演化时标,所选探测器的型号以及所采用的数据采集帧频也有所不同。以上因素都会导致各通道观测到的太阳图像存在一定的视场差异,很难达到与其空间分辨率相当的视场匹配精度。下文将对NVST多波段成像系统各通道太阳图像之间的视场匹配方法作详细介绍。

图 1. (a) NVST建筑剖面图;(b)光学平台示意图

Fig. 1. (a) NVST architectural section; (b) diagram of optical platform

下载图片 查看所有图片

图 2. NVST多波段高分辨率成像观测系统光路示意图

Fig. 2. Schematic of optical path of NVST multiband high resolution imaging observation system

下载图片 查看所有图片

2.2 视场匹配变换模型与参数计算方法

首先,NVST多波段成像系统在经过高精度装调之后,探测器靶面的法线与光路严格重合,各通道视场之间仅存在二维平面上的几何变换关系;其次,多波段成像系统的观测视场较小(<3'),根据NVST多通道实际光路的ZEMAX仿真分析得到整个视场的畸变小于0.157%^[6],于是可以不考虑大视场的视场高阶畸变^[7-9]。在上述条件下,可以用图像之间的旋转、放缩以及平移的仿射变换来描述各通道视场的不匹配性,称之为RST变换。用旋转因子θ、尺度因子m,以及平移因子t_x和t_y 4个参数来表示这种变换关系,即四参数模型。在该模型的假设下,如果有足够多的高精度特征点,就可利用图像点匹配的原理并通过最小二乘估计求解变换四参数。

设r(x_r,y_r)为参考通道中的特征点,p(x_p,y_p)为变换通道与参考通道对应的特征点,两点之间的匹配变换可表示为

r(xr,yr)=p[tx+m(xrcosθ-yrsinθ),ty+m(xrsinθ+yrcosθ)],(1)

(1)式对应线性方程的矩阵表达式为

xryr1=mcosθmsinθtx-msinθmcosθty001xpyp1。(2)

一般地,(2)式可表示为

xryr1=abc-bad001xpyp1+δxδy1,(3)

式中δ_x和δ_y分别表示仿射变换残余误差,a、b、c和d为矩阵四参数。矩阵四参数和RST变换四参数的关系为

θ=arctanbam=a2+b2tx=cty=d。(4)

若有N个对应特征点,对(3)式进行扩充可得

xr1yr1xr2yr2︙xrNyrN=xp1yp110-yp1xp101xp2yp210-yp2xp201︙︙︙︙xpNypN10-ypNxpN01abcd+δx1δy1δx2δy2︙δxNδyN,(5)

式中(x_r1,y_r1),(x_r2,y_r2),…,(x_rN,y_rN)分别为N个对应特征点。根据(5)式,运用最小二乘法求解矩阵四参数a, b, c和d,再由(4)式得到仿射变换四参数θ,m,t_x和t_y,并对参数误差作出估计。

2.3 针孔阵列光阑

基于点特征的图像匹配方法中一个重要的问题是精确探测和提取图像中的匹配特征点。在这方面前人已经研究了很多算法,例如Moravec^[10]于1977年提出基于图像灰度自相关函数的特征点检测算子;Lowe^[11]于1999年提出了尺度不变特征变换的特征点检测方法;You等^[12]在2000年提出基于小波变换的图像匹配算法等。以上方法不仅需要图像的内容,而且往往要求图像之间有较高的相似度。对相似度较低的图像,提取对应特征点的有效性会大大降低。对于NVST多波段成像观测而言,各通道采集到的图像包括太阳光球像和色球像,不同波段观测到的太阳大气活动差异很大,甚至可能完全不同,故不能直接通过上述方法检测和提取特征点。为此,在望远镜的F3处放置一个二维针孔阵列光阑,采集针孔阵列像,以针孔像的质心作为多通道视场定标的匹配特征点来求解变换参数。

在针孔光阑的设计过程中,以TiO通道为标准对孔径大小、孔间间距以及加工精度提出相应的要求。为了保证针孔像质心求解的准确度,要求针孔像半宽需要在10 pixel左右,根据TiO通道相对于F3的垂轴放大率R(R≈0.58)以及探测器像元大小S(S=6.5 μm)计算可得所需设计的针孔孔径大小ϕ=0.1 mm。为使针孔光阑加工误差不影响后续的计算,孔径大小以及孔间间距的误差均需控制在1 pixel以下,于是针孔光阑实际的加工精度需在0.01 mm左右。为在保证计算精度的前提下获得尽可能多的针孔,孔间距设计为3 mm。最终针孔阵列光阑设计为11×11点阵的矩形阵列,孔间距为(3±0.01) mm,针孔大小为(0.1±0.01) mm,可覆盖整个视场。图3给出了直径为50 mm的圆形针孔阵列光阑示意图。图3中红点表示针孔偏离标准阵列的位置,据此可标记出各通道针孔像的对应关系,并克服各通道图像的镜像或90°旋转的对称性。为了能明显定位标记针孔像的位置而又不对其他针孔像产生干扰,所设计的标记针孔点偏离阵列0.8 mm。

图 3. 针孔阵列光阑的示意图

Fig. 3. Schematic of pinhole array aperture

下载图片 查看所有图片

3 双通道视场定标与数据分析

3.1 双通道视场定标的数据采集

双通道视场定标过程中,TiO通道和Hα通道均使用英国Andor公司的Zyla CMOS相机,靶面包含2560×2160个像元,每个像元大小为6.5 μm×6.5 μm。由于NVST是一架地平式望远镜,为了消除像场旋转的影响,NVST采取了使光学平台与像场同步旋转的消旋方案来保持图像稳定^[13-14]。因此,在双通道视场定标中,需考虑光学平台的转动对视场匹配四参数的影响。针孔像的采集包括:1)光学平台处于静止状态下观测一组针孔阵列像(约5 min),用来计算特定位置处的视场匹配四参数;2)旋转光学平台模拟消旋过程,在光学平台处于旋转状态下观测一组针孔阵列像(约7 min),研究消旋过程中匹配参数是否发生变化以及变化的大小,以确定一天之内视场定标的次数,保证通道间太阳图像的匹配精度。

为了将匹配参数应用到实测太阳图像上并验证视场匹配的精度,在采集完针孔阵列像之后采集了一组太阳图像。需要强调的是,在太阳图像的采集过程中,由于地球大气特性的不断变化以及望远镜镜体晃动等因素的影响,各通道数据的视场也会不断变化^[15-16]。为了克服这些因素的影响并有效验证视场定标对双通道视场匹配的效果,采用了同步观测模式采集太阳图像,即通过探测器外触发工作模式实现两通道图像的同步采集,以保证两通道具有相同的曝光时间^[17]。选取了日面中心附近的一个小活动区作为观测目标,分别进行Hα线心色球像(656.28 nm)与TiO像的同步观测(2016-10-06 03:46:52—03:52:03 UT);Hα远线翼光球像(656.28~656.30 nm)与TiO像的同步观测(2016-10-06 03:53:00—03:58:45 UT)。

图4是双通道观测到的针孔阵列像的灰度反像图,图4(a)是Hα通道针孔阵列像,图4(b)是对应的TiO通道针孔阵列像。从图4可看出Hα通道的视场内包含了9×9阵列的针孔像,而TiO通道的视场只包含8×8阵列的针孔像,直观地说明了两个通道的视场放大率存在明显的差异。

图 4. (a) Hα通道获得的针孔阵列像的灰度反像图;(b) TiO通道获得的针孔阵列像的灰度反像图

Fig. 4. (a) Grayscale anti-image of pinhole array image obtained by Hα channel; (b) grayscale anti-image pinhole array obtained by TiO channel

下载图片 查看所有图片

3.2 双通道RST变换四参数计算及误差分析

采用质心法求解双通道各针孔的中心坐标值来得到匹配控制点。为了进一步降低随机噪声对计算结果的影响,对得到的450组针孔阵列像进行叠加平均,得到一组平均针孔阵列像,以此计算匹配控制点。从图4可看出,实验中的实际有效针孔数目为64,因此可得到64个有效控制点。由64个控制点计算得到Hα通道相对于TiO通道的RST变换四参数值以及相应的标准差分别为θ=(0.00299±8.807×10^-5) rad,m=1.155±1.017×10^-4,t_x=(-198.214±0.171) pixel,t_y=(-167.255±0.171) pixel。

经计算可得双通道视场匹配的综合误差为0.594 pixel,对应太阳表面0.031″,完全满足0.1″的匹配精度要求。

为了进一步考察误差的分布,对Hα通道针孔阵列像作RST变换。比较变换后各针孔质心位置与相应的TiO通道针孔质心位置,计算它们的残余误差(即各对应质心位置之间的坐标差)。图5是对针孔阵列像进行RST变换之后各针孔质心位置残余误差的矢量分布。从图5可以看出,质心位置残差在整个视场中的分布不均匀,残差畸变分布,残差的标准差为0.59 pixel,对应太阳表面0.031″;最大残差位于右上角,可达1.46 pixel,对应太阳表面0.076″,依然满足匹配误差精度要求。残差畸变有可能是由光路上的一些光学元件造成的光学像差所导致的。

图 5. Hα通道针孔阵列像经过四参数RST变换视场定标后的残余误差分布

Fig. 5. Residual error distribution of Hα channel pinhole array image after RST transformation field calibration of four parameters

下载图片 查看所有图片

由于光学平台存在载荷形变,在消旋过程中会产生图像视场的漂移^[18]。为了考察光学平台处于不同位置时视场定标参数的变化情况,在光学平台转动的同时连续采集针孔阵列像,计算得到一系列RST变换四参数,由此反映出匹配四参数在整个消旋过程中的变化情况。图6是光学平台消旋过程中(从东限位0°到西限位约190°)四参数的变化情况,图中有几段数据点缺失,这是由于观测过程中有云遮挡,图像信噪比过低,无法正确计算这段时间内的针孔质心位置。图6(a)是旋转因子的变化情况,图6(b)是尺度因子的变化情况,图6(c)是x轴平移因子的变化情况,图6(d)是y轴平移因子的变化情况。表1给出了四参数变化的均值和标准差。

图6中红线是对四参数变化趋势的拟合结果,通过拟合结果得到在消旋过程中,四参数变化造成的视场匹配差异不超过0.05″。因此,在对实测太阳图像进行视场匹配时,可忽略光学平台的转动对四参数的影响,选取某一固定位置计算匹配参数,并以此对全天的观测数据进行视场定标。

图 6. 消旋过程中RST变换四参数的变化情况。(a)旋转因子;(b)尺度因子;(c) x轴平移因子;(d) y轴平移因子

Fig. 6. Changes of four parameters after RST transformation in de-rotation process. (a) Rotation actor; (b) scaling factor; (c) x-axis translation factor; (d) x-axis translation factor

下载图片 查看所有图片

图 7. (a)经过RST变换之后的Hα线心色球像;(b)与图7(a)对应的TiO同步帧图像;(c)经过RST变换之后的Hα远线翼光球像(656.28~656.30 nm);(d)与图7(c)对应的TiO同步帧图像

Fig. 7. (a) Hα line center chromosphere image after RST transformation; (b) TiO synchronous frame image corresponding to figure 7(a); (c) Hα far line-wing photosphere image after RST transformation (656.28~656.30 nm); (d) TiO synchronous frame image corresponding to figure 7(c)

下载图片 查看所有图片

图 8. 经RST变换之后Hα远线翼光球像与对应TiO 同步帧图像之间的局部相关偏移量矢量分布图

Fig. 8. Local correlation offset vector distribution of Hα far line-wing photosphere image after RST transformation compared with the corresponding TiO synchronous frame image

下载图片 查看所有图片

4 双通道太阳数据的匹配结果

用所得四参数对Hα线心色球图像和Hα远线翼光球图像进行RST变换,并分析它们与各自对应的同步TiO像的匹配效果。图7(a)是Hα通道经过RST变换之后的线心色球像,图7(b)是与图7(a)对应的同步TiO像,图7(c)是Hα通道经过RST变换之后的远线翼光球像(656.28~656.30 nm),图7(d)是与图7(c)对应的同步TiO像。因为视场中存在黑子,所以可直观地看到两通道视场匹配的效果,但在实际的科学观测中Hα通道的观测像为线心色球像,且通常与TiO光球像没有明显的相似性。由于Hα远线翼数据与TiO数据都为太阳光球像,两者具有直接的可比性。为了定量验证针孔光阑视场定标对双通道实测太阳图像的视场匹配精度,利用该组数据进行计算分析。将图7(c)、(d)两帧同步数据进行局部相关位移计算,用相关偏移量反映视场定标匹配对实测太阳图像的匹配精度。图8是局部相关偏移量的矢量图,所有相关偏移量在x方向和y方向上的平均值分别为0.346 pixel和0.470 pixel。从图8中可以看出局部相关偏移量的分布同样是不均匀的。经过对240组同步数据进行局部相关位移量计算,得出每组局部相关偏移量的平均值分布在0.3~0.5 pixel之间,精度小于0.03″,局部最大偏移量在1~2 pixel之间,精度小于0.1″。造成局部偏移量偏大的因素包括视场定标后残余的视场畸变以及双通道太阳图像采集时地球大气湍流在不同波段对太阳图像影响不同导致图像的局部视场偏移存在差异。

5 结论

提出了NVST多波段成像系统各通道之间的视场匹配方法,即针孔光阑视场定标方法,并在Hα通道和TiO通道中进行了一系列实验与分析,得出以下结论:

1) 在有64个有效控制点的情况下,视场匹配精度可达0.031″。

2) 双通道经过视场定标之后仍存在残余的视场匹配畸变,且匹配残差的标准差为0.031″,在视场边缘最大可达0.076″,但仍然在匹配精度要求的范围内。

3) 光学平台在消旋过程中处于不同位置时,计算得到的视场匹配参数会有所不同,但最大差值不超过0.05″。因此,可以不考虑不同位置处匹配参数的差异,即可用光学平台处于某一固定位置处的匹配参数对全天的观测数据进行视场定标匹配。

4) 以对太阳Hα远线翼光球像与同步观测TiO像的匹配为准,用局部相关偏移量验证视场定标在实测数据上的效果,同样得到偏移量的平均值在0.03″以下,局部区域会有0.1″的偏移。

因此,可认为在四参数模型的假设下,用所提方法对双通道进行视场定标之后,即便不考虑光路自身高阶像差带来的视场畸变和地球大气湍流造成的图像局部偏移在不同观测波段的差异,也能实现两通道之间小于0.1″的视场匹配精度。所提方法不依赖于图像的特征结构,可在实际观测中高精度地匹配不同波段结构特征完全不同的图像。该方法不仅适用于NVST多波段成像系统的观测,也可为其他类似望远镜的视场匹配提供借鉴。