楊 盼，曾曙光，劉 鎖，鄭 勝，林鋼華，何慧靈

(1. 三峽大學(xué)理學(xué)院，湖北 宜昌 443000；2. 中國科學(xué)院太陽活動重點實驗室，北京 100012)

太陽磁場在太陽活動中扮演了極為關(guān)鍵的角色。為了預(yù)測耀斑和日冕物質(zhì)拋射等太陽爆發(fā)活動，須監(jiān)測從太陽光球到日冕的磁場和流場結(jié)構(gòu)及其演化過程[1]。目前，國內(nèi)外已經(jīng)建有多個太陽磁場觀測設(shè)備。如國內(nèi)懷柔太陽觀測基地(Huairou Solar Observing Station, HSOS)的35 cm太陽磁場望遠(yuǎn)鏡(Solar Magnetic Field Telescope, SMFT)[2]、國外太陽動力天文臺(Solar Dynamics Observatory, SDO)的日震及磁場成像儀(Helioseismic and Magnetic Imager, HMI)[3]等。不同太陽觀測站間存在觀測時間、觀測范圍、觀測設(shè)備等差異，為了方便科學(xué)研究，提高太陽物理數(shù)據(jù)的使用效率，有必要找到不同觀測站間太陽圖像自動配準(zhǔn)的方法。另外，當(dāng)兩觀測設(shè)備的觀測范圍不同時，還需要實現(xiàn) “小圖像” 在 “大圖像” 中的自動定位。例如，懷柔太陽觀測基地的太陽磁場望遠(yuǎn)鏡得到的是局部太陽磁場圖像，然而，太陽動力學(xué)天文臺的日震及磁場成像儀得到的是全日面太陽磁場圖像。前者的觀測范圍遠(yuǎn)遠(yuǎn)小于后者，這就需要找到一種能夠?qū)⒕植刻枅D像在全日面圖像中快速定位的方法。

目前已有一些關(guān)于太陽圖像配準(zhǔn)方面的研究工作，如文[4]提出了基于修正矩的亞像素太陽圖像配準(zhǔn)算法，文[5]采用基于信息熵與尺度不變特征變換(Scale invariant feature Transform, SIFT)算法對太陽圖像進(jìn)行了配準(zhǔn)。這些都是針對同一觀測站的序列圖像進(jìn)行的，其成像設(shè)備位置固定，圖像內(nèi)容在短時間內(nèi)變化較小，位移、旋轉(zhuǎn)和仿射變換較少。不同太陽觀測站的太陽圖像由于在觀測時間、觀測范圍、觀測設(shè)備等方面存在較大差異，需要適應(yīng)力更強(qiáng)的圖像配準(zhǔn)方法進(jìn)行配準(zhǔn)。已有的圖像配準(zhǔn)方法很多，如基于灰度的配準(zhǔn)方法[6-7]、基于特征的配準(zhǔn)方法[8-9]等等?；诨叶鹊呐錅?zhǔn)方法具有精度高的優(yōu)點，但計算復(fù)雜度高，且對目標(biāo)的旋轉(zhuǎn)、形變、遮擋以及灰度變化等很敏感?；诔叨炔蛔兲卣鼽c匹配的圖像配準(zhǔn)方法[10-11]計算量小，對灰度的變化、圖像形變及遮掩等有較好的適應(yīng)能力。太陽磁場圖像含有許多內(nèi)部結(jié)構(gòu)，特征點較明顯，本文提出了基于尺度不變特征點匹配的太陽磁場圖像配準(zhǔn)與定位方法，開展了懷柔太陽觀測基地的局部太陽磁場圖像在太陽動力天文臺日震及磁場成像儀的全日面太陽磁場圖像中自動配準(zhǔn)與定位的研究，并做了檢測實驗證明方法的有效性。

1 算法原理

尺度不變特征點匹配算法對局部太陽圖像在全日面太陽中自動配準(zhǔn)與定位，主要包括尺度不變特征提取算法的特征點檢測、特征點初步匹配、特征點精確匹配、圖像變換矩陣求解4個步驟。

1.1 尺度不變特征變換特征點檢測

1.1.1 尺度空間極值檢測

尺度空間理論是檢測不變特征的基礎(chǔ)。文[12]證明了高斯卷積核是實現(xiàn)尺度變換的唯一變換核。

一幅二維圖像在不同尺度下的尺度空間可以表示為圖像與高斯核卷積：

L(x,y,σ)=G(x,y,σ)*I(x,y)，

(1)

其中，(x,y)為圖像點的像素坐標(biāo)；I(x,y)為圖像灰度值；G(x,y,σ)為高斯核函數(shù)。

(2)

其中，σ為尺度空間因子，它是高斯正態(tài)分布的方差，反映了圖像被平滑的程度；L(x,y,σ)為圖像的尺度空間。

為高效地在尺度空間檢測穩(wěn)定的極值，文[10]使用尺度空間中高斯差分算子(Difference of Gaussian, DoG)的極值作為判斷依據(jù)。高斯差分算子定義如下：

(3)

其中，k為相鄰兩個尺度空間的比例因子。

高斯差分算子局部極值的檢測是將每個像素的高斯差分算子值跟同一尺度的周圍鄰域8個像素和相鄰尺度對應(yīng)位置的周圍鄰域9 × 2個像素總共26個像素進(jìn)行比較。僅當(dāng)被檢測點的高斯差分算子值大于26個像素點或小于26個像素點時，才將其判斷為極值點并保存，以進(jìn)行后續(xù)計算。

1.1.2 確定關(guān)鍵點位置及尺度

通過擬合三維二次函數(shù)以精確確定關(guān)鍵點的位置和尺度，從而可得到原圖像尺度不變特征提取算法的候選點集合。其中，對比度低的點對噪聲比較敏感，而位于邊緣上的點難以準(zhǔn)確定位。為了確保尺度不變特征提取算法特征點的穩(wěn)定性，須從尺度不變特征提取算法候選點集合將這兩類點剔除。

1.1.3 關(guān)鍵點方向確定

通過確定關(guān)鍵點的方向，可以使特征描述符以與方向相關(guān)的方式構(gòu)造，從而使算子具有旋轉(zhuǎn)不變性。關(guān)鍵點的方向利用鄰域像素的梯度分布特性確定。對每個高斯圖像，每個點L(x,y)的梯度方向θ(x,y)可以由下式計算[10]：

(4)

1.1.4 特征向量生成

為了確保描述符的旋轉(zhuǎn)不變性，首先需旋轉(zhuǎn)坐標(biāo)，保持坐標(biāo)方向與極值點方向一致；然后將平面內(nèi)可能分布的360°梯度方向劃為8個方向范圍，每個方向范圍各包含45°。在子區(qū)域為4 × 4的情況下，共計有4 × 4 × 8=128個數(shù)據(jù)，即生成128維描述符。方向與信息聯(lián)合不僅增強(qiáng)了抗噪性，還消除了尺度變化、旋轉(zhuǎn)和變形的影響[5,10]。

1.2 特征點初步匹配

完成特征點檢測后，可利用尋找同名點的方法對兩幅圖像特征點進(jìn)行粗匹配[5,10]。以歐式距離d作為圖像間的相似性度量以確定同名點。歐式距離二維表示為

(5)

首先在待配準(zhǔn)圖像中找出兩個與參考圖像中對應(yīng)某極值點的歐氏距離最短的點；然后以兩個值中小的值作為分子，大的值作為分母，若該分?jǐn)?shù)值在給定的閾值范圍內(nèi)，那么建立對應(yīng)點的同名點關(guān)系，稱為匹配點。閾值的選取將直接影響同名點的數(shù)量。閾值越大，獲得的匹配點越多，可能的誤匹配也越多；閾值過小時，匹配點會很少。大量實驗表明，閾值在0.4～0.6之間最佳[5]。

1.3 特征點精確匹配及變換矩陣求解

在完成特征點粗匹配后，可利用隨機(jī)抽樣一致性算法[9]消除錯誤的匹配點，從而提高配準(zhǔn)精度。隨機(jī)抽樣一致性算法，可以在一組包含 “外點” 的數(shù)據(jù)集中，采用不斷迭代的方法，尋找最優(yōu)參數(shù)模型，不符合最優(yōu)模型的點，被定義為 “外點”。該算法最早由文[9]提出，在圖像配準(zhǔn)以及拼接上得到廣泛應(yīng)用。

一般地，兩幅圖像間的變換可以由變換矩陣表示[13]：

(6)

其中，(x′,y′)和(x,y)分別為兩幅圖像中對應(yīng)點的像素坐標(biāo)；hi(i=0, 1, 2, ..., 7)為變換矩陣的8個未知系數(shù)。

隨機(jī)抽樣一致性算法可對特征點進(jìn)行精確匹配并求變換矩陣。首先隨機(jī)選擇4組對應(yīng)點組成一個隨機(jī)樣本并計算變換矩陣H，對假設(shè)的每組對應(yīng)計算距離d，計算與H一致的內(nèi)點數(shù)，選擇有最多內(nèi)點數(shù)的H，在數(shù)目相等時，選擇內(nèi)點標(biāo)準(zhǔn)方差最小的那個解[13]。最后根據(jù)圖像間變換矩陣H，可以計算出局部太陽圖像的4個頂點及中心點在全日面圖像中的像素坐標(biāo)，從而完成局部太陽圖像在全日面太陽圖像中的配準(zhǔn)與定位。

2 算法設(shè)計及實驗分析

2.1 算法設(shè)計

基于尺度不變特征匹配的太陽磁場圖像配準(zhǔn)與定位方法流程如下：

(1)對圖像進(jìn)行預(yù)處理，包括降采樣、黑白反相、對比度增強(qiáng)等，對圖像降采樣到1 024 × 1 024左右；

(2)通過尺度不變特征提取算法提取圖像特征點；

(3)采用尋找同名點的方法對兩幅圖像特征點進(jìn)行粗匹配；

(4)利用隨機(jī)抽樣一致性算法去除錯誤的匹配點，并估計兩幅圖像之間的最佳變換矩陣H；

(5)根據(jù)仿射變換矩陣，實現(xiàn)局部太陽圖像在全日面太陽圖像中的初步定位；

(6)利用初步定位的圖像區(qū)域和懷柔局部圖像進(jìn)行精確配準(zhǔn)。

算法的具體流程如圖1。

圖1 太陽圖像配準(zhǔn)與定位流程圖
Fig.1 Flow chart of solar image registration and positioning

2.2 實驗結(jié)果與分析

2.2.1 不同觀測站太陽磁場圖像的配準(zhǔn)與定位

以懷柔太陽觀測基地的局部太陽磁場圖像在SDO/HMI的全日面太陽磁場圖像中的自動匹配與定位展開實驗。圖2給出了原始的全日面太陽圖像及局部太陽圖像。其中，圖2(a)給出日震及磁場成像儀[3]得到的全日面太陽磁場圖像，拍攝時間2012年11月12日05時48分00秒。圖2(b)給出35 cm太陽磁場望遠(yuǎn)鏡[2]得到的局部太陽磁場圖像，拍攝時間2012年11月12日05時55分34秒。兩者的拍攝時間非常接近。值得指出的是，兩者的觀測范圍和分辨率有很大差異。懷柔太陽觀測基地的局部太陽圖像僅僅是全日面太陽圖像中很小的一部分。全日面圖像的尺度是4 096 × 4 096，懷柔站局部圖像的尺度830 × 992。

圖2 原始的太陽磁場圖像
Fig.2 Original solar magnetic field image (Full-disk solar magnetic field image in (a) and local solar magnetic field image (in(b)))

圖2中兩幅圖像的對比度不夠。在進(jìn)行后續(xù)的圖像配準(zhǔn)前，有必要通過對比度增強(qiáng)等對兩幅圖像進(jìn)行預(yù)處理。另外，由于日震及磁場成像儀的全日面圖像分辨率太高，為了提高配準(zhǔn)速度，先做降采樣處理，將尺度從4 096 × 4 096降低到1 024 × 1 024。最后，部分日震及磁場成像儀的全日面太陽磁場圖像的正負(fù)極性與懷柔觀測站的相反，因此要做黑白反相處理。經(jīng)過降采樣、黑白反相及對比度增強(qiáng)處理后的全日面太陽磁場圖像如圖3(a)。經(jīng)過對比度增強(qiáng)后的局部太陽圖像如圖3(b)。對比圖2(b)和圖3(b)可以發(fā)現(xiàn)，預(yù)處理后的太陽圖像比原始圖像更加清晰，這有利于后續(xù)的特征提取與圖像配準(zhǔn)。

圖3 經(jīng)過預(yù)處理后的太陽磁場圖像
Fig.3 Full-disk solar magnetic field image (in(a)) and local solar magnetic field image (in(b)) after pre-treatment

隨后，采用尺度不變特征提取算法提取圖像特征點。圖4給出了經(jīng)過預(yù)處理后的局部太陽圖像的尺度不變特征檢測結(jié)果。一共檢測到了204個尺度不變特征點。由圖4可見，這些特征點較好地代表了局部太陽圖像的特征。

在對圖3中的兩幅圖像進(jìn)行尺度不變特征檢測后，需要將兩幅圖的特征點配準(zhǔn)。圖3中的全日面圖像一共檢測到了3 447個尺度不變特征點，局部太陽圖像一共檢測到了204個尺度不變特征點。接下來，采用尋找同名點的方法以及隨機(jī)抽樣一致性算法找出兩幅圖中相互配準(zhǔn)的特征點對。如圖5，一共找到了20個相互配準(zhǔn)的特征點對(圖中用紅色直線連接)。根據(jù)這些配準(zhǔn)特征點對，可初步估計兩幅圖像之間的變換矩陣H1[13]，具體如(7)式：

(7)

根據(jù)變換矩陣H1可以計算局部圖像的4個頂點及中心點在全日面圖像中的像素位置，從而完成了局部圖像在全日面圖像中的粗定位(圖5)。

圖4 經(jīng)過預(yù)處理后的局部太陽圖像的尺度不變特征檢測

Fig.4 Scale-invariant feature detection of localized solar image after pre-treatment

圖5 局部圖像在全日面圖像中的配準(zhǔn)與定位
Fig.5 Registration and positioning of local image in full-disk solar image

對日震及磁場成像儀全日面圖像進(jìn)行降采樣處理可以大大提高計算速度，但同時也會降低匹配與定位精度。為此，將日震及磁場成像儀全日面原始圖像中綠色四邊形對應(yīng)的最小外接矩形區(qū)域(見圖5)與懷柔局部太陽圖像做精確匹配，即可確定局部太陽圖像在全日面圖像中的精確范圍。精確匹配后的結(jié)果如圖6，一共找到了25個相互配準(zhǔn)的特征點對(圖中用紅色直線連接)。根據(jù)這些配準(zhǔn)特征點對，可估計兩幅圖像之間的最佳變換矩陣H2，具體如(8)式：

(8)

由變換矩陣H2可以看出，兩圖像間主要存在縮放、旋轉(zhuǎn)、平移等幾何變換，不存在透射變換，且可以計算局部圖像的4個頂點在全日面圖像中的精確位置分別為(2 909, 2 704)、(2 371, 2 484)、(2 560, 2 045)、(3 097, 2 265)(圖6中用綠色四邊形標(biāo)出)，從而完成了局部圖像在全日面圖像中的精確定位。另外，還可以計算局部太陽圖像相對于全日面圖像的旋轉(zhuǎn)角度約為157.9°(順時針)；水平方向放大比例為1.71，垂直方向放大比例為1.74。

圖6 局部圖像在全日面圖像中的精確配準(zhǔn)與定位
Fig.6 Accurate registration and positioning of local image in full-disk solar image

圖7給出了局部放大后的全日面太陽圖像(圖7(a))與經(jīng)過變換矩陣H變換后的局部太陽圖像(圖7(b))。由圖可見，兩圖的形狀、角度等基本一致，但細(xì)節(jié)方面略有差異。兩圖形狀、角度等基本一致證明了本文算法的準(zhǔn)確性。兩圖細(xì)節(jié)方面略有差異是因為兩者是由不同觀測站得到的，在觀測時間、望遠(yuǎn)鏡參數(shù)等方面存在差異所致。但這些細(xì)微的差異不會影響局部太陽圖像在全日面圖像中的配準(zhǔn)與定位。

圖7 全日面太陽磁場圖像與經(jīng)過變換后的局部太陽磁場圖像
Fig.7 Full-disk solar magnetic field image (in(a)) and transformed local solar magnetic field image (in(b))

2.2.2 不同時間段的匹配結(jié)果

為了驗證算法的有效性，對71組不同時段的太陽圖像進(jìn)行配準(zhǔn)，并給出相應(yīng)的配準(zhǔn)點對數(shù)與正確配準(zhǔn)點對數(shù)，結(jié)果如圖8。其中，紅色實線圓圈為實際的配準(zhǔn)點對數(shù)，綠色虛線點為正確配準(zhǔn)的點對數(shù)，藍(lán)色虛線五角星為配準(zhǔn)的正確率(即正確配準(zhǔn)的點對數(shù)與實際配準(zhǔn)的點對數(shù)之比)，全日面圖像分辨率均為1 024 × 1 024。由圖8可以看出，大部分情況下正確配準(zhǔn)的點對數(shù)與實際的配準(zhǔn)點對數(shù)相等，即配準(zhǔn)準(zhǔn)確率為100%；少數(shù)情況下正確配準(zhǔn)的點對數(shù)小于實際的配準(zhǔn)點對數(shù)，即存在一定的錯誤配準(zhǔn)。出現(xiàn)錯誤配準(zhǔn)的原因在于：懷柔局部太陽圖像和日震及磁場成像儀全日面太陽圖像，一個是地面拍攝的，一個是空間拍攝的。在大氣湍流等干擾下，懷柔局部太陽圖像會變得模糊，尺度不變特征變換算法能檢測到的特征點減少，進(jìn)而導(dǎo)致配準(zhǔn)點對數(shù)和配準(zhǔn)準(zhǔn)確率的下降。圖9給出了局部太陽圖像較模糊時的配準(zhǔn)情況，當(dāng)懷柔局部太陽圖像較為模糊時，配準(zhǔn)點對數(shù)與配準(zhǔn)準(zhǔn)確率較低。

圖8 不同時間太陽磁場圖像的配準(zhǔn)情況
Fig.8 The registration of the solar magnetic field image at different time

圖9 局部太陽磁場圖像較模糊時的配準(zhǔn)情況
Fig.9 Registration results of a blurry local solar magnetic field image

2.2.3 仿射變換對太陽磁場圖像配準(zhǔn)結(jié)果的影響

不同觀測站的太陽磁場圖像間除了常見的縮放、旋轉(zhuǎn)和平移變換，很多情況下還存在比較明顯的仿射變換。本節(jié)研究仿射變換對太陽磁場圖像配準(zhǔn)結(jié)果的影響，通過定量分析并計算配準(zhǔn)實驗的誤差。從原始日震及磁場成像儀全日面圖像中(圖2(a))摳出一塊矩形子圖像。然后，將矩形子圖像進(jìn)行一定的仿射變換，最后，將仿射變換后的子圖像與降采樣后的全日面圖像(圖3(a))做配準(zhǔn)，將變換矩陣計算得到的子圖像在全日面圖像中的位置與子圖像在全日面圖像中的已知位置相減，即可準(zhǔn)確地計算配準(zhǔn)誤差。配準(zhǔn)誤差η定義為

(9)

其中，(xi0,yi0)和(xi,yi)分別為子圖像的4個頂點在全日面圖像中的已知位置和經(jīng)過變換矩陣預(yù)測的位置。

仿射變換包含平移、旋轉(zhuǎn)、縮放、錯切等變換。圖10給出了縮放為0.6倍的情況下，配準(zhǔn)誤差(單位：像素)隨錯切系數(shù)和旋轉(zhuǎn)角度的變化情況。由圖10可以看出，當(dāng)錯切系數(shù)較小時，配準(zhǔn)誤差較小，而當(dāng)錯切系數(shù)增大到一定程度如超過0.3后，配準(zhǔn)誤差將迅速增大；旋轉(zhuǎn)角度對配準(zhǔn)結(jié)果影響較小。圖11給出了旋轉(zhuǎn)角度為45°的情況下，配準(zhǔn)誤差隨錯切系數(shù)和縮放系數(shù)的變化情況。由圖11可以看出，放大2倍以內(nèi)，對配準(zhǔn)誤差的影響較小，超出此范圍后配準(zhǔn)誤差將明顯增大。這是因為縮小變換需要丟棄一些圖像信息，而放大變換則需要插值來增大圖像。當(dāng)縮小導(dǎo)致丟棄的信息太多或放大插值的信息太多時對特征點的檢測與配準(zhǔn)造成影響。綜上所述，本文采用的尺度不變特征變換配準(zhǔn)算法對旋轉(zhuǎn)、縮放、錯切變換等均具有一定的魯棒性。為了更好地反映仿射變換對配準(zhǔn)精度的影響，本節(jié)的結(jié)果都是通過初步定位得到的，通過精確定位可以大大減小配準(zhǔn)誤差。

圖10 配準(zhǔn)誤差隨錯切系數(shù)和旋轉(zhuǎn)角度的變化
Fig.10 The change of registration error varies with the shear factor and the rotation angle

圖11 配準(zhǔn)誤差隨錯切系數(shù)和縮放系數(shù)的變化
Fig.11 The change of registration error varies with the shear factor and zoom factor

3 結(jié) 論

采用基于尺度不變特征點匹配的方法對懷柔觀測站的局部太陽磁場圖像在太陽動力天文臺日震及磁場成像儀的全日面圖像中的配準(zhǔn)與定位問題進(jìn)行了研究，結(jié)果表明，該算法對光照、旋轉(zhuǎn)、尺度縮放等有很強(qiáng)的魯棒性，對不同觀測站得到的具有不同觀測條件、不同觀測尺度等的太陽磁場圖像的配準(zhǔn)有效且提高了太陽物理數(shù)據(jù)的使用效率。

[1] 劉忠, 鄧元勇, 季海生, 等. 中國地基大太陽望遠(yuǎn)鏡[J]. 中國科學(xué): 物理學(xué) 力學(xué) 天文學(xué), 2012, 42(12): 1282-1291.

Liu Zhong, Deng Yuanyong, Ji Haisheng, et al. China Large Earth-based Solar Telescope[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2012, 42(12): 1282-1291.

[2] 艾國祥, 胡岳風(fēng). 太陽磁場望遠(yuǎn)鏡的工作原理[J]. 天文學(xué)報, 1986(2): 91-98.

Ai Guoxiang, Hu Yuefeng. On principle of solor magnetic field telescope[J]. Acta Astronomica Sinica, 1986(2): 91-98.

[3] Schou J, Scherrer P H, Bush R I, et al. Design and ground calibration of the Helioseismic and Magnetic Imager(HMI) Instrument on the Solar Dynamics Observatory (SDO) [J]. Solar Physics, 2012, 275(1-2): 229-259.

[4] 陳潔. 太陽圖像的亞像素配準(zhǔn)算法的研究與應(yīng)用[D]. 昆明: 昆明理工大學(xué), 2013.

[5] 岳昕, 尚振宏, 強(qiáng)振平, 等. 基于信息熵與SIFT算法的天文圖像配準(zhǔn)[J]. 計算機(jī)科學(xué), 2015, 42(6): 57-60.

Yue Xin, Shang Zhenhong, Qiang Zhenping, et al. Astronomical image registration combining information entropy and SIFT algorithm[J]. Computer Science, 2015, 42(6): 57-60.

[6] Goshtasby A, Stockman G C, Page C V. A region-based approach to digital image registration with subpixel accuracy[J]. IEEE Transactions on Geoscience & Remote Sensing, 1986, GE-24(3): 390-399.

[7] Pratt W K. Correlation techniques of image registration[J]. IEEE Transactions on Aerospace & Electronic Systems, 1974, AES-10(3): 353-358.

[8] Bouchiha R, Besbes K. Automatic remote-sensing image registration using SURF[C]// The 3rd International Conference on Machine Vision. 2013: 406-410.

[9] Fischler M A, Bolles R C. Random sample consensus: a paradigm for model fitting with a pplications to image analysis and automated cartography[J]. Readings in Computer Vision, 1987, 24(6): 726-740.

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

[11]Lowe D G. Object recognition from local scale-invariant features[C]// The Proceedings of the Seventh IEEE International Conference on Computer Vision. 1999.

[12]Koenderink J J. The structure of images[J]. Biological Cybernetics, 1984, 50(5): 363-370.

[13]趙小川. 現(xiàn)代數(shù)字圖像處理技術(shù)提高及應(yīng)用案例詳解[M]. 北京: 北京航空航天大學(xué)出版社, 2012.