太陽(yáng)磁場(chǎng)觀測(cè)中相關(guān)位移疊加算法的比較*

陳 潔，馮 松，鄧 輝，季凱帆

(昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650500)

對(duì)太陽(yáng)而言，微弱磁場(chǎng)只有幾高斯[1]，要得到如此高的靈敏度，需要通過(guò)CCD長(zhǎng)時(shí)間積分或者多幀疊加獲取足夠信噪比的圖像。但如果使用長(zhǎng)積分或者直接短曝光圖像疊加，由于大氣抖動(dòng)和風(fēng)的影響使得圖像存在線性位移，無(wú)法得到幾高斯的分辨率，甚至?xí)绊懭醮艌?chǎng)的探測(cè)[1]。例如懷柔太陽(yáng)觀測(cè)基地35 cm太陽(yáng)磁場(chǎng)望遠(yuǎn)鏡在單幀觀測(cè)時(shí)間為20 ms時(shí)，只能觀測(cè)到百高斯以上的強(qiáng)磁場(chǎng)，要得到幾高斯的弱磁場(chǎng)，至少需要在一個(gè)極性方向上疊加256幀[2]。在太陽(yáng)觀測(cè)中，通常使用相關(guān)跟蹤技術(shù)消除圖像的線性位移，雖然該技術(shù)非常有效，但它需要對(duì)望遠(yuǎn)鏡硬件做較大的調(diào)整，而懷柔太陽(yáng)磁場(chǎng)望遠(yuǎn)鏡結(jié)構(gòu)緊湊，在硬件上沒(méi)有調(diào)整的空間，因此，需要從軟件上考慮。最佳的方法是將觀測(cè)得到的圖像通過(guò)圖像相關(guān)計(jì)算偏移量，然后進(jìn)行位移疊加[2]的方法實(shí)時(shí)消除圖像的線性位移。

本文通過(guò)比較和分析相位相關(guān)和互相關(guān)算法的圖像配準(zhǔn)算法，對(duì)模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的處理，認(rèn)為相位相關(guān)算法對(duì)噪聲更加敏感，而互相關(guān)算法更適應(yīng)于磁場(chǎng)的位移疊加觀測(cè)。為了比較兩種算法分別配準(zhǔn)后的圖像效果，選取觀測(cè)差分圖像的平均值和均方差、配準(zhǔn)圖像的最大相關(guān)值、疊加圖像的熵和磁場(chǎng)圖像的能量五個(gè)統(tǒng)計(jì)量分別評(píng)價(jià)。

1 相關(guān)算法

圖像配準(zhǔn)技術(shù)廣泛應(yīng)用于工業(yè)、天文圖像處理、醫(yī)學(xué)影像學(xué)、航空航天等領(lǐng)域，圖像匹配算法是這些技術(shù)的核心。常用的圖像匹配算法可以分為兩大類(lèi)：基于灰度相關(guān)的匹配算法和基于特征相關(guān)的匹配算法[3-6]。

基于灰度相關(guān)的圖像匹配算法是把一幅圖像的其中一個(gè)像素的灰度鄰域作為模板在另外一幅圖像中搜索具有相同灰度值分布的對(duì)應(yīng)鄰域。這種算法具有簡(jiǎn)單易實(shí)現(xiàn)和能夠獲得較高定位精度等特點(diǎn)。常用的算法有相位相關(guān)算法、互相關(guān)算法、序貫相似性檢測(cè)匹配算法(SSDA)、交互信息算法等；基于特征相關(guān)的匹配算法從待匹配的圖像中提取特征，通常利用空間位置相對(duì)不變的特性(例如邊緣、輪廓、紋理等)進(jìn)行匹配。其優(yōu)點(diǎn)是大大壓縮了圖像的匹配信息，計(jì)算量小、速度快。但由于它的匹配性能很大程度取決于特征的提取，因此匹配的精度不高。這種算法常用于幾何特征單一的圖像。常用的有基于興趣點(diǎn)的匹配、基于矩陣特征的匹配等[3-6]。

雖然圖像匹配算法眾多，但太陽(yáng)觀測(cè)圖像具有灰度變化較明顯和幾何特征并不單一的特點(diǎn)，而且在實(shí)際觀測(cè)中，由于光子噪聲和CCD讀出噪聲的影響，特征點(diǎn)的提取非常困難。因此，基于特征相關(guān)的匹配算法并不適合太陽(yáng)磁場(chǎng)觀測(cè)圖像的處理。而交互信息算法目前主要用于醫(yī)學(xué)圖像處理，而序貫相似性檢測(cè)匹配算法在匹配時(shí)需要選擇閥值，匹配的精度很大程度上依賴(lài)閥值的選擇[7-8]。因此，暫時(shí)不考慮交互信息算法和序貫相似性檢測(cè)匹配算法。互相關(guān)算法和相位相關(guān)算法均可以通過(guò)快速傅里葉變換將時(shí)域問(wèn)題轉(zhuǎn)化為頻域問(wèn)題，可以滿足實(shí)時(shí)處理的需求。另外通過(guò)對(duì)實(shí)測(cè)圖像進(jìn)行匹配，發(fā)現(xiàn)懷柔太陽(yáng)磁場(chǎng)觀測(cè)基地短時(shí)間采集的圖像并不存在明顯的圖像畸變、像場(chǎng)旋轉(zhuǎn)和比例尺變化，只存在微小的線性位移，這些微小的線性位移主要來(lái)自風(fēng)、跟蹤誤差和大氣抖動(dòng)的一階成分[1,4]。因此選取相位相關(guān)算法和互相關(guān)算法作為懷柔太陽(yáng)磁場(chǎng)觀測(cè)基地的候選匹配算法。

1.1 相位相關(guān)算法(Phase Correlation, PC)

假設(shè)f(x，y)、g(x，y)為兩張要進(jìn)行相關(guān)的圖像，其中g(shù)(x，y)是f(x，y)平移(x0，y0)后的圖像。

g(x,y)=f(x-x0,y-y0)

(1)

由傅里葉變換的平移性

G(u,v)=F(u,v)e-j(ux0+vy0)

(2)

(2)式中的G(u，v)和F(u，v)是g(x，y)、f(x，y)分別經(jīng)過(guò)傅里葉變換得到的結(jié)果。

f(x，y)和g(x，y)的互功率譜[8]為：

(3)

式中，G*(u,v)是G(u，v)取共軛得到的。對(duì)(3)式求逆傅里葉變換，然后確定峰值的位置(x0，y0)，峰值代表兩個(gè)圖像的相關(guān)度，f(x，y)與g(x，y)的線性位移為(x0，y0)。

1.2 互相關(guān)算法(Cross Correlation, CC)

互相關(guān)算法是一種基于統(tǒng)計(jì)的配準(zhǔn)方法，由Rosenfeld在1982提出[9]。對(duì)于給定大小為(m，n)的圖像g(x，y)、h(x，y)，兩者的相關(guān)函數(shù)的定義如下：

(4)

從上面的相關(guān)函數(shù)的定義式中可以證明：

F{Rgh(s,t)}=F{h(x,y)}×F*{g(x,y)}=H(u,v)×G*(u,v)=R(u,v)

(5)

(5)式中，H(u,v)是h(x,y)的傅里葉變換形式；G(u,v)是g(x,y)的傅里葉變換形式。對(duì)G(u,v)取共軛，得到G*(u,v)，然后與H(u,v)進(jìn)行點(diǎn)乘，即可獲得相關(guān)函數(shù)的傅里葉變換R(u,v)。將R(u,v)做傅里葉逆變換，并得到逆變換矩陣中的最大值所在的位置，即為兩幅圖相關(guān)性最強(qiáng)的位置。

一般來(lái)講，這兩種算法對(duì)噪聲都有一定的容忍度，文[10]作者從理論上對(duì)這兩種算法進(jìn)行過(guò)仔細(xì)的推導(dǎo)，認(rèn)為無(wú)論對(duì)于高頻噪聲還是低頻噪聲，互相關(guān)算法相比較而言都優(yōu)于相位相關(guān)算法。懷柔基地也一直采用互相關(guān)算法進(jìn)行圖像配準(zhǔn)[2]。但針對(duì)太陽(yáng)觀測(cè)圖像，尚未有仔細(xì)的研究，因此本文從模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)兩個(gè)方面對(duì)這兩種算法進(jìn)行了系統(tǒng)的比對(duì)。

2 模擬圖像處理實(shí)驗(yàn)

實(shí)驗(yàn)采用的計(jì)算機(jī)配置為4 GB內(nèi)存，英特爾奔騰雙核的中央處理器(Central Processing Unit, CPU)，主頻為2.7 GHz，顯卡為英偉達(dá)公司的GT430。模擬圖像為中心在不同位置的多組64×64的負(fù)高斯圖像，類(lèi)似于黑子圖像，并添加均值為0、標(biāo)準(zhǔn)差不同的高斯分布噪聲。采用相位相關(guān)和互相關(guān)算法分別計(jì)算其偏移量，并與理論值進(jìn)行比對(duì)。在沒(méi)有噪聲的情況下，相位相關(guān)和互相關(guān)算法測(cè)試得出的結(jié)果和實(shí)際的偏移量完全一致。但隨著噪聲的增加，兩種算法的結(jié)果出現(xiàn)明顯差異，對(duì)每組噪聲水平進(jìn)行了100組測(cè)量，測(cè)試結(jié)果如圖1。其中X軸表示加入的峰值信噪比(PSNR)，定義為：

(6)

式中，max表示圖像顏色的最大值；MSE(Mean Square Error)表示均方誤差。

式中，M、N分別表示模擬圖像的尺寸；I(i,j)和K(i,j)表示含噪聲圖像和原始圖像。

圖1 不同噪聲水平下相位相關(guān)和互相關(guān)算法處理模擬圖像的對(duì)比結(jié)果

Fig.1 Performance evaluations of simulated image processings with the CC and PC for various noise levels

圖1中的Y軸表示100組相同噪聲水平的測(cè)試中，相位相關(guān)和互相關(guān)中與實(shí)際的偏移量相同的次數(shù)。

從(6)式分析得出，峰值信噪比越小時(shí)，均方誤差越大，表明加入的噪聲越多，信噪比越低。從圖1可以明顯看到，當(dāng)PSNR<10 db時(shí)，相位相關(guān)和互相關(guān)算法均無(wú)法測(cè)出實(shí)際的偏移量。當(dāng)PSNR=15 db時(shí)，兩個(gè)算法可以測(cè)準(zhǔn)少量的值，隨著峰值信噪比的增大，結(jié)果越來(lái)越理想。當(dāng)峰值信噪比從20 db到25 db時(shí)，互相關(guān)算法的測(cè)試結(jié)果迅速提高，然后維持在一個(gè)較為理想的水平。而相位相關(guān)算法提高的較為緩慢，到PSNR=50 db以后才能沒(méi)有誤差。

從模擬數(shù)據(jù)可以看到，互相關(guān)算法比相位相關(guān)算法有較好的抗噪能力，這和前人的理論結(jié)果[10]是一致的。

3 太陽(yáng)磁場(chǎng)實(shí)測(cè)數(shù)據(jù)處理實(shí)驗(yàn)

實(shí)驗(yàn)采用的樣本來(lái)自于懷柔35 cm太陽(yáng)磁場(chǎng)望遠(yuǎn)鏡的實(shí)測(cè)圖像，CCD和望遠(yuǎn)鏡的參數(shù)見(jiàn)表1。獲得了9個(gè)觀測(cè)時(shí)段圖像，每組圖像分左右旋光，合計(jì)18組圖像，每組包括100幅短曝光圖像，用于疊加，這些圖像的峰值信噪比約為33 db。分別使用相位相關(guān)和互相關(guān)算法對(duì)各組圖像進(jìn)行測(cè)試。由于實(shí)測(cè)圖像并沒(méi)有理論偏移量，因此需要對(duì)處理結(jié)果進(jìn)行多方位的比較和評(píng)價(jià)。

表1 CCD和望遠(yuǎn)鏡的參數(shù)

評(píng)價(jià)圖像處理結(jié)果是一個(gè)較為復(fù)雜的問(wèn)題，直接用肉眼從匹配后的圖像中很難分辨出哪個(gè)算法更有優(yōu)勢(shì)。為了進(jìn)行統(tǒng)計(jì)分析，選取了5個(gè)統(tǒng)計(jì)量，分別是差分圖像平均值μ、差分圖像均方差σ、配準(zhǔn)圖像的最大相關(guān)值、位移疊加后圖像的熵、磁場(chǎng)圖像的能量分別對(duì)兩種算法的結(jié)果進(jìn)行評(píng)價(jià)。

3.1 差分圖像平均值和均方差

以第1幅圖像為基準(zhǔn)圖像，然后計(jì)算后續(xù)圖像位移后的圖像與此基準(zhǔn)圖像的差分圖像。理論上，如果沒(méi)有噪聲，將配準(zhǔn)圖按偏移量移動(dòng)后，應(yīng)與基準(zhǔn)圖像對(duì)齊，每個(gè)對(duì)應(yīng)點(diǎn)像素的值相同，兩張圖相減的殘差圖的所有像素點(diǎn)值為0。同時(shí)由于配準(zhǔn)圖和基準(zhǔn)圖受不同的光照條件或是曝光精度的影響，因此采用兩幅圖像相除的辦法得到差分圖像。差分圖像平均值μ表示經(jīng)過(guò)平移后的配準(zhǔn)圖和基準(zhǔn)圖的所有像素比例的均值，這個(gè)值越接近1，說(shuō)明匹配過(guò)程中圖像越穩(wěn)定。而均方差σ表示經(jīng)過(guò)平移后的配準(zhǔn)圖和基準(zhǔn)圖的所有像素比例的均方差，這個(gè)值越小，說(shuō)明匹配的效果越好。

首先計(jì)算未經(jīng)過(guò)位移的各個(gè)圖像的差分均值和均方差，然后計(jì)算經(jīng)過(guò)相位相關(guān)和互相關(guān)算法位移以后差分圖像的均值和均方差，并和未位移差分圖像測(cè)量的結(jié)果進(jìn)行歸一化處理(圖2)。

從圖2(a)可以得出，相位相關(guān)和互相關(guān)算法的殘差均值基本一致，互相關(guān)算法要稍微小一些。從圖2(b)可以得出，互相算法的均方差要明顯小于相位相關(guān)算法的結(jié)果。這說(shuō)明相位相關(guān)算法得到的結(jié)果的差分圖總體起伏大于互相關(guān)算法。

圖2 相位相關(guān)和互相關(guān)算法均值(a)、均方差(b)的對(duì)比圖

Fig.2 The averages (a) and standard deviations (b) for the errors in the results of processings with the CC and PC

3.2 配準(zhǔn)圖像的極大相關(guān)性

從上述的分析可以得出，不管是相位相關(guān)算法還是互相關(guān)算法得出的偏移量坐標(biāo)，均為傅里葉逆變換后的最大值坐標(biāo)，即為相關(guān)系數(shù)的最大值，所以也可以通過(guò)比較相位相關(guān)和互相關(guān)算法配準(zhǔn)圖像所求得的相關(guān)系數(shù)最大值的大小來(lái)評(píng)價(jià)圖像處理結(jié)果[11]。

圖3(a)是18組數(shù)據(jù)的相位相關(guān)和互相關(guān)算法進(jìn)行配準(zhǔn)疊加后相關(guān)系數(shù)值與未位移而疊加的圖像直接計(jì)算的相關(guān)系數(shù)的比值，每組數(shù)據(jù)是由100張圖配準(zhǔn)后的均值。從圖3(a)中得出互相關(guān)算法的相關(guān)性值均要略大于相位相關(guān)算法的相關(guān)性值，也就是通過(guò)互相關(guān)算法位移疊加后得到相關(guān)系數(shù)大于相位相關(guān)算法的結(jié)果。

3.3 疊加圖像的熵

熵(Entropy)是由Rudolf Clausius提出的[12]，用來(lái)描述體系的混亂程度。后來(lái)Shannon將熵的概念應(yīng)用到信息學(xué)理論中，對(duì)于離散的隨機(jī)信號(hào)I，有N種隨機(jī)值，且發(fā)生的概率為pi，則信號(hào)的信息熵為：

(8)

圖像的熵值表示圖像信源的平均信息量，由于熵值代表系統(tǒng)的無(wú)序程度，當(dāng)系統(tǒng)越不穩(wěn)定，熵值越大，則圖像中灰度概率分布越均勻，圖像的細(xì)節(jié)越明顯，圖像的紋理越清晰。反之，圖像的熵值越小，圖像的紋理越模糊。也就是如果位移疊加的效果越好，結(jié)果圖像的細(xì)節(jié)越清晰，圖像的熵越高。

圖3(b)統(tǒng)計(jì)的是18組數(shù)據(jù)的經(jīng)過(guò)兩種算法進(jìn)行位移疊加后的圖像的熵和未經(jīng)過(guò)位移而直接疊加的圖像的熵的比值。

圖3 相位相關(guān)和互相關(guān)算法的相關(guān)值(a)、熵值(b)對(duì)比圖

Fig.3 The correlation coefficients (a) and entropy values (b) (relative to the respective true values) for the results of processings with the CC and PC

從圖3(b)可以看出，互相關(guān)算法得到的疊加圖像的熵要比相位相關(guān)算法得到的熵的比值要大，說(shuō)明互相關(guān)算法配準(zhǔn)后的圖像紋理更為清晰，細(xì)節(jié)更多。

3.4 磁場(chǎng)圖像的能量

首先根據(jù)左右旋光的觀測(cè)資料，采用不同的位移疊加算法，得到不同的磁場(chǎng)圖。

圖4統(tǒng)計(jì)的是9組疊加的磁場(chǎng)圖的能量比較，從圖4中可以得出，除第3組和第6組相位相關(guān)后的磁場(chǎng)圖的能量略大于互相關(guān)算法外，其他組互相關(guān)算法得到的磁場(chǎng)圖比相位相關(guān)算法得到的磁場(chǎng)圖的能量要大。由于沒(méi)有對(duì)原始圖像進(jìn)行暗場(chǎng)和平場(chǎng)改正，所以CCD上面的系統(tǒng)誤差會(huì)影響最后的能量計(jì)算結(jié)果。但總體上看，互相關(guān)算法還是優(yōu)于相位相關(guān)算法。

4 討論與下一步工作

圖4 相位相關(guān)和互相關(guān)算法的能量對(duì)比圖

Fig.4 The image energy values (relative to the respective true values) for the results of processings with the CC and PC

針對(duì)兩種常用的相關(guān)位移算法，從模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的對(duì)比結(jié)果上看，互相關(guān)算法在處理太陽(yáng)圖像時(shí)均優(yōu)于相位相關(guān)算法。對(duì)于實(shí)測(cè)數(shù)據(jù)，通過(guò)多種評(píng)價(jià)指標(biāo)對(duì)差分圖像、位移疊加圖像和最終的磁場(chǎng)圖像進(jìn)行了比較。從算法原理上看，互相關(guān)算法的計(jì)算量也小于相位相關(guān)算法。目前圖形處理器(Graphic Processing Unit, GPU)技術(shù)已經(jīng)非常成熟，利用GPU技術(shù)進(jìn)行快速傅里葉變換并配以相應(yīng)的目標(biāo)處理窗口，是完全可以達(dá)到實(shí)時(shí)處理的要求。

通過(guò)實(shí)驗(yàn)可以看到，兩種算法都依賴(lài)于圖像的信噪比，但同時(shí)也發(fā)現(xiàn)他們還依賴(lài)于圖像灰度的分布情況，因此僅僅通過(guò)峰值信噪比并不能精確確定匹配精度。如何通過(guò)原始圖像的灰度分布估計(jì)匹配精度是需要進(jìn)一步研究的問(wèn)題。

但是無(wú)論是互相關(guān)算法和相位相關(guān)算法，都是通過(guò)獲取相關(guān)系數(shù)的最大值得到位移量的，因此是一個(gè)整數(shù)象元的圖像匹配算法。這對(duì)于目前的地面觀測(cè)，例如懷柔基地的磁場(chǎng)觀測(cè)圖像，這個(gè)精度是可以滿足要求的。但隨著觀測(cè)質(zhì)量的進(jìn)一步提高和空間天文學(xué)的發(fā)展，圖像匹配向亞象元方向發(fā)展是一個(gè)必然的趨勢(shì)，也是值得深入研究的問(wèn)題。事實(shí)上，互相關(guān)算法由于得到是一個(gè)相關(guān)系數(shù)的輪廓，而相位相關(guān)算法一般得到的是一個(gè)近似脈沖的信號(hào)，因此通過(guò)互相關(guān)算法更有利于亞象元的測(cè)量，這也是進(jìn)一步準(zhǔn)備開(kāi)展的工作。

致謝：衷心感謝中國(guó)科學(xué)院國(guó)家天文臺(tái)懷柔觀測(cè)基地鄧元勇研究員、林佳本博士等的大力支持和協(xié)助。

[1]林佳本, 鄧元勇, 胡柯良. 局部相關(guān)跟蹤算法在高分辨太陽(yáng)磁場(chǎng)觀測(cè)中的應(yīng)用與實(shí)現(xiàn)[J]. 計(jì)算機(jī)工程與應(yīng)用, 2006(27): 203-205.

Lin Jiaben, Deng Yuanyong, Hu Keliang. Application of LCT in solar magnetic field observation[J]. Computer Engineering and Applications, 2006(27): 203-205.

[2]Deng Yuanyong, Zhang Bin. Application of a large capacity real-time image acquisition system. I. “Deep integration” magnetogram of quiet region on the Sun[J]. Astrophysics Reports, 1999(34): 15-19.

[3]Heipke C. Overview of image matching technique[J]. OEEPE Official Publications, 1993(10): 173-189.

[4]林佳本, 鄧元勇, 胡柯良, 等. 實(shí)時(shí)相關(guān)跟蹤圖像處理系統(tǒng)[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2006, 3(4): 345-354.

Lin Jiaben, Deng Yuanyong, Hu Keliang, et al. Real-time lmage-processing system with correlation tracking[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2006, 3(4): 345-354.

[5]Barbara Zitova, Jan Flusser. Image registration methods: a survey, Image and Vision Computing[J]. Image and Vision Computing, 2003(21): 977-1000.

[6]Brown L G. A survey of image registration techniques[J]. ACM computing Surveys, 1992, 24(4): 325-376.

[7]Barnea D I, Silverman H F. A class of algorithm for fast digital image registration[J]. IEEE Transactions on Computers, 1972, C-21(2): 179-186.

[8]Srinivasa B Reddy, Chatterji 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.

[9]Lewis J P. Fast normalized cross-correlation[J]. Circuits, System and Signal Processing, 2009(28): 819-843.

[10]Manduchi R, Mian G A. Accuracy analysis for correlation-based image registration algorithms[J]. Circuits and Systems, 1993(1): 834-837.

[11]Huynh-Thu, Ghanbari M. Scope of validity of PSNR in image/video quality assessment[J]. Electronics Letters, 2008, 44(13): 800-801.

[12]A Rényi. On measures of entropy and information[J]. Mathematics of Statistics and Probability, 1961(1): 547-561.