劉　輝，楊云飛，尚振宏，李潤鑫

(1. 昆明理工大學信息工程與自動化學院，云南 昆明 650500; 2. 昆明理工大學云南省計算機技術(shù)應用重點實驗室，云南 昆明 650500)

太陽觀測已經(jīng)全面進入高空間和高時間分辨率的時代。以日本HINODE衛(wèi)星[1]、中國1 m新真空太陽望遠鏡(New Vacuum Solar Telescope，NVST)[2]、美國Goode Solar Telescope[3]和瑞典Swedish Solar Telescope[4]為代表的高分辨望遠鏡已經(jīng)產(chǎn)生了海量的多波段觀測序列資料，為太陽物理的研究帶來重要的機遇和挑戰(zhàn)[5]。這些資料直觀地從強度上展現(xiàn)了太陽大氣中時刻發(fā)生的各種尺度的活動現(xiàn)象(太陽風暴、大小尺度活動以及存在于各種尺度中的磁流體動力學波和震蕩)的細節(jié)及變化情況。高時空分辨率的圖像不但能從形態(tài)上展示太陽大氣的強度變化，同時也為量化分析和定量描述帶來了可能。其中，通過強度像隨著時間的變化計算這些活動的橫向速度場已被廣泛應用于太陽光球、色球表面特征的動力學分析。

太陽光球和色球的橫向運動速度的典型值在每秒幾千米到幾十千米。對于目前的高分辨太陽圖像序列的空間分辨率和時間分辨率，這意味著要精確測量圖像兩幀之間亞像素級或者幾個像素級的位移量。在太陽圖像的數(shù)據(jù)處理領(lǐng)域，這方面已經(jīng)有了很多的發(fā)展。最為經(jīng)典的有局部相關(guān)跟蹤法(Local Correlation Tracking, LCT)[6]以及后來的基于傅里葉變換的局部相關(guān)跟蹤法(Fourier Local Correlation Tracking, FLCT)[7]。微分仿射速度估計(Differential Affine velocity estimator, DAVE)[8-9]等技術(shù)在引入磁場感應方程后，首先被用于磁圖非勢性研究，然后被應用于強度像的測量[10-11]。此外文[12-13]提出了一種Balltracking方法，用于光球超米粒結(jié)構(gòu)的研究。雖然這些方法都得到了一些很好的結(jié)果，但在測量精度上還不能匹配1 m太陽望遠鏡這種像元分辨率在0.05 arcsec、時間分辨率高達十多秒的高分辨圖像處理要求，因而有必要對相關(guān)方法進行更為深入的研究。

橫向速度場的測量實際上是通過測量一定時間間隔的兩幀圖像之間特征的位移量。由于太陽序列圖像是某種 “光” 強度隨著時間的變化，因此速度場是一種典型的光流場。從另外一方面，兩幀圖像之間的位移測量也可以看作是一種圖像非剛性的配準。而隨著計算機圖像處理技術(shù)的發(fā)展，光流場的測量和圖像非剛性配準技術(shù)已經(jīng)有了非常大的發(fā)展。

圖像非剛性配準方法可以由3部分描述：(1)聯(lián)系源圖像和目標圖像的空間變換；(2)測量目標圖像和源圖像相似性的相似性測度；(3)決定最優(yōu)變換參數(shù)的優(yōu)化方法。

非剛性配準方法主要有基于空間變換的配準方法和基于物理模型的配準方法兩大類。其中，由文[14]提出的基于熱力學中流體擴散理論的Demons非剛性配準方法具有處理非光滑變形、抗噪聲性能好等優(yōu)點，被廣泛應用于醫(yī)學圖像配準[15]，也非常符合太陽光球和色球?qū)拥拇帕黧w非均勻局部變形的物理特性。

本文用Demons方法測量1 m太陽望遠鏡高分辨光球和色球觀測數(shù)據(jù)的橫向速度場。測量結(jié)果表明該方法可以較為精確地測量太陽高分辨圖像的橫向速度場，優(yōu)于傳統(tǒng)的基于傅里葉變換的局部相關(guān)跟蹤法和微分仿射速度估計方法。

1　Demons方法

1.1　Demons 配準的原理和流程

經(jīng)典Demons算法由Thirion提出，即 “Demons-base” 算法[14]。在概念上，它和19世紀麥克斯韋(Maxwell)的實驗原理很相似。該算法的基本思想是假設運動目標的灰度不隨時間改變，那么非剛性配準可以看作是浮動圖像中各個像素向參考圖像逐步擴散的過程，浮動圖像各個像素的擴散速度由參考圖像的灰度梯度信息決定。假設參考圖像r、浮動圖像f都是圖像函數(shù)I[x(t),y(t),t]在某一時刻的 “快照”，圖像函數(shù)的灰度值保持常數(shù)，則此假設可被模型化為

I[x(t),y(t),t]=c,

(1)

f=I[x(t0),y(t0),t0],

(2)

r=I[x(t1),y(t1),t1].

(3)

在初始時刻t0，圖像函數(shù)I等于f，經(jīng)過一定時間到達t1時刻后，圖像函數(shù)被完全變形為r。圖像配準的過程就是要得出一個能驅(qū)動f中的各個像素點向r中對應像素點移動的向量場。為了得出驅(qū)動力，將(1)式偏微分，得到：

(4)

(4)式可以簡化為

v·r=f-r,

(5)

(6)

若將浮動圖像的梯度信息作為一種正內(nèi)力，參考圖像的梯度信息作為負內(nèi)力，利用這兩種力同時驅(qū)動形變，得到的擴散速度為

(7)

Demons算法利用圖像局部信息驅(qū)動形變，在全局范圍內(nèi)使形變連續(xù)，從而保證圖像的拓撲結(jié)構(gòu)。

1.2　基于梯度互信息的Demons方法

在Demons模型中，浮動圖像各個像素都可以自由移動，使得在浮動圖像中具有某個特定灰度值的所有像素有可能映射到參考圖像上的同一像素點，導致錯誤的配準結(jié)果，為了使圖像像素能夠沿著正確的方向移動，在(7)式的基礎(chǔ)上，采用梯度互信息進一步改進Demons算法，在原有擴散速度的基礎(chǔ)上，增加了兩幅圖像之間梯度互信息的作用，使兩幅圖像配準結(jié)束的同時梯度互信息也達到最大。改進后的Demons算法擴散速度模型為

(8)

其中，max[IMI(vn)]為兩幅圖像的梯度互信息；β為正常數(shù)，表示此項的權(quán)重。

基于圖像灰度互信息的配準方法，最大的不足是忽略了圖像的空間信息。就太陽圖像配準過程而言，對相同目標所成圖像可能具有不同的灰度、不同的分辨率，甚至圖像本身大小也不同，但是對于同一目標，其邊界是確定的，不會發(fā)生很明顯的變化；當兩幅圖像配準時，對應位置像素點的梯度矢量將具有相同或者相反的指向?；谝陨戏治?，梯度互信息采用Parzen窗方法進行計算。

2　實驗數(shù)據(jù)

1 m太陽望遠鏡作為中國最大的太陽望遠鏡和世界三大高分辨成像太陽望遠鏡之一，其優(yōu)良的觀測性能和高質(zhì)量的圖像已經(jīng)在很多方面得以證實。為檢驗Demons方法的可用性，選取了具有代表性的3個1 m太陽望遠鏡觀測資料數(shù)據(jù)集，包括兩組光球圖像和一組色球圖像，每組圖像為前后兩幀。通過計算這兩幀圖像的光流場或者對它們進行非剛性配準，可得到相應的速度場。

數(shù)據(jù)集1包括的兩幀圖像均為1 m太陽望遠鏡在TiO波段觀測的光球Level 1 + 重建像，觀測時間為2013年11月3日04時58分31秒(UT)和05時01分28秒(UT)，時間間隔2分57秒。視場大小為685 × 650像元，像元分辨率為0.04 arcsec。光球上典型的1 km/s的橫向速度將帶來約6像元的位移。意味著整體光流場是超像素運動。

數(shù)據(jù)集2與數(shù)據(jù)集1是同目標、同波段、同期的觀測，只是第1幀圖像的觀測時間為05時1分10秒(UT)，與第2幀相隔僅18秒，大約是數(shù)據(jù)集1的1/10，其他參數(shù)與數(shù)據(jù)集1相同。這樣1 km/s的橫向速度將帶來僅僅0.6像元的位移，即光流場為亞像素運動。

數(shù)據(jù)集3的兩幀圖像為1 m太陽望遠鏡在Hα波段的Level 1 + 重建像。觀測時間為2017年8月6日09時59分53秒和10時00分25秒，時間間隔為32秒。視場大小為238 × 225像元，像元分辨率0.13 arcsec。色球上5 km/s的橫向速度帶來約1.7像元的位移，總體光流場中包括亞像素和超像素的位移。

從圖1可以看到，光球像數(shù)據(jù)集1和2是一個太陽活動區(qū)的局部，特征結(jié)構(gòu)覆蓋了米粒、黑子本影和半影。數(shù)據(jù)集3是典型的色球纖維，圖像結(jié)構(gòu)本身就有很強的自相似性。

圖1(a)數(shù)據(jù)集1和2中的第2幀；(b)數(shù)據(jù)集3中的第2幀
Fig.1The second frame of dataset 1 and 2 (a) and the second frame of dataset 3 (b)

3　數(shù)據(jù)處理結(jié)果

對上述3個數(shù)據(jù)集采用Demons方法計算第2幀相對于第1幀的逐像素位移量，即光流場，然后根據(jù)相應的時間間隔歸算到逐點的橫向速度場。圖2顯示了數(shù)據(jù)集1中一個小區(qū)域的兩幀圖以及用Demons方法計算出的位移矢量。米粒的細微運動可以清晰地反映在圖上。

圖2數(shù)據(jù)集1中同一區(qū)域的變化及用Demons方法測量到的矢量場
Fig.2The same region of two frames in dataset 1 and the vector fields measured by Demons

為了測試Demons的結(jié)果并與目前太陽觀測資料處理中常用的基于傅里葉變換的局部相關(guān)跟蹤法和微分仿射速度估計方法進行比對，在分別采用不同方法得到光流場后，采用雙線性插值將第2幀與第1幀進行非剛性配準，并利用圖像結(jié)構(gòu)相似度指數(shù)(Structural Similarity Index Measurement, SSIM)[16]測試配準的精度。較高的配準精度即意味著較高的光流場(速度場)測量精度。

對于x，y兩張圖像，其圖像結(jié)構(gòu)相似度指數(shù)由下式求得：

(9)

為了逐像素比對非剛性配準后的兩幀圖像的相似度，采用11 × 11像素鄰域的窗口平均作為該點的局部圖像結(jié)構(gòu)相似度指數(shù)值。對于配準的結(jié)果，如果有更多點的圖像結(jié)構(gòu)相似度指數(shù)值接近1，說明該方法總體位移測量更為準確，配準更好。而圖像中所有像素點圖像結(jié)構(gòu)相似度指數(shù)值的平均作為兩幀圖像的平均相似度值，也代表整體測量水平。

表1給出了3個數(shù)據(jù)集在未配準前和使用基于傅里葉變換的局部相關(guān)跟蹤法、微分仿射速度估計和Demons方法分別測量并非剛性配準后與第1幀圖像比較的平均圖像結(jié)構(gòu)相似度指數(shù)值。

從表1可以看出，原始數(shù)據(jù)彼此的相似度不但取決于時間間隔，也取決于特征結(jié)構(gòu)的運動速度。對于色球數(shù)據(jù)，快速的運動導致在半分鐘內(nèi)圖像就有較大的變化。另外，使用基于傅里葉變換的局部相關(guān)跟蹤法配準后，圖像結(jié)構(gòu)相似度指數(shù)提升非常有限，即只測量到很少的運動。再則，微分仿射速度估計方法的測量結(jié)果有顯著的提升，但Demons方法更優(yōu)。

表1不同方法配準后的平均圖像結(jié)構(gòu)相似度指數(shù)值

Table1TheSSIMindexafternon-rigidregistrationwithdifferentmethod

數(shù)據(jù)集1數(shù)據(jù)集2數(shù)據(jù)集3未配準074100988406731FLCT074400988506881DAVE083320992107532Demons090720996508298

數(shù)據(jù)集1間隔時間較長，而數(shù)據(jù)集3的色球運動較快。在這兩個數(shù)據(jù)集中，特征結(jié)構(gòu)除了局部的運動外，還有新結(jié)構(gòu)的浮現(xiàn)和舊結(jié)構(gòu)的消失，因此圖像結(jié)構(gòu)相似度指數(shù)相對偏小。而數(shù)據(jù)集2中兩幀圖像時間間隔僅為18 s，且光球運動較為緩慢，圖像差異微小，圖像結(jié)構(gòu)相似度指數(shù)也明顯高于其他兩個數(shù)據(jù)集。但無論原始圖像情況如何，Demons方法對于不同的數(shù)據(jù)集均表現(xiàn)出良好的處理能力。

圖3顯示了對數(shù)據(jù)集1應用3種方法配準后的逐點圖像結(jié)構(gòu)相似度指數(shù)的直方圖概率分布?？梢钥吹紻emons方法配準后圖像結(jié)構(gòu)相似度指數(shù)值的整體分布更靠近1，明顯好于微分仿射速度估計和基于傅里葉變換的局部相關(guān)跟蹤法。同時基于傅里葉變換的局部相關(guān)跟蹤法的結(jié)果和未配準圖像非常接近，即測量效果并不顯著。這與表1的結(jié)果是一致的。數(shù)據(jù)集2和3的圖像結(jié)構(gòu)相似度指數(shù)概率分布也反映出同樣的趨勢。

從最終速度場的測量情況看，3個數(shù)據(jù)集都呈現(xiàn)同樣的趨勢： 基于傅里葉變換的局部相關(guān)跟蹤法基本上只能在位移較大的情況下有測量值，而微分仿射速度估計的測量值系統(tǒng)性小于Demons的測量值，大多數(shù)情況只有Demons方法測量精度的一半。

圖3數(shù)據(jù)集1的不同方法配準后的逐點圖像 結(jié)構(gòu)相似度指數(shù)概率分布

Fig.3The probability distribution of SSIM index with different registration method applying dataset 1

4　模擬數(shù)據(jù)測試

Demons方法雖然結(jié)果好于微分仿射速度估計和基于傅里葉變換的局部相關(guān)跟蹤法，但第2幀經(jīng)過非剛性配準后也并沒有做到完全匹配，即圖像結(jié)構(gòu)相似度指數(shù)并不等于1。其中最為重要的原因無論是光球還是色球，經(jīng)過一段時間間隔后，不但原有的米粒、黑子或者色球纖維產(chǎn)生位移，同時也有新特征的產(chǎn)生和舊結(jié)構(gòu)的消失，這不是僅從光流場能夠計算的。

為了進一步比較Demons和微分仿射速度估計的光流場測量精度，測試他們在不同太陽圖像(光球和色球)上亞像素和超像素位移的可靠性，設計了兩個模擬實驗，即分別將數(shù)據(jù)集1和數(shù)據(jù)集3中的第1幀圖像在X和Y方向整體移動一個固定量，從而產(chǎn)生第2幀剛性運動的圖像，然后分別再用兩個方法測量其與第1幀圖像的光流場。由于模擬的是全局性位移，因此其逐點光流場的均值和標準差即可以很好地反映測量的精度。

表2給出了整體位移為(5.7, -0.2)時，這兩個方法在X和Y方向上測量的統(tǒng)計結(jié)果。

表2　模擬數(shù)據(jù)使用微分仿射速度估計和Demons方法測量的統(tǒng)計結(jié)果Table2　ThestatisticsresultsofDAVEandDemonsforsimulationdata

可以看到，與預設值相比，Demons的均值和標準差表現(xiàn)都明顯好于微分仿射速度估計方法，逐點測量精度在0.1像素。而微分仿射速度估計不但標準差比Demons方法高一個量級，而且均值也有較大的系統(tǒng)偏差，尤其對于自相關(guān)較強的色球圖像。一般認為微分仿射速度估計在做超像素測量時誤差較大，但測試結(jié)果表明微分仿射速度估計受到圖像本身的特征結(jié)構(gòu)的影響也很大。這與微分仿射速度估計本質(zhì)上也是使用相關(guān)方法測量有關(guān)。基于傅里葉變換的局部相關(guān)跟蹤法由于能計算的像素太少，缺乏統(tǒng)計意義，因此被排除在外。

同時，由于是模擬整圖的剛性位移，因此也可以用標準的整幅圖像互相關(guān)方法直接測量兩組圖像的剛性位移。光球像的測量結(jié)果為(5.757 8, -0.181 4)，而色球像的結(jié)果為(5.744 6， -0.065 4)，也比用Demons方法逐點測量后的均值結(jié)果明顯偏差要大。這意味著Demons方法不但適用于非剛性的配準，同時也可以利用均值測量圖像剛性的位移，用于圖像的剛性配準。尤其對于自相似性較強的目標，如色球纖維、太陽邊緣等，這一方法有較強的優(yōu)勢。

5　結(jié)　論

本文詳細描述了Demons方法，并利用其測量不同目標和不同位移量的太陽高分辨像的橫向速度場。對測量結(jié)果進行圖像非剛性配準后，圖像結(jié)構(gòu)相似度指數(shù)有了明顯的提高，從而說明位移的測量更為準確，表明這一方法的可行性。通過與微分仿射速度估計和基于傅里葉變換的局部相關(guān)跟蹤法進行對比，Demons方法表現(xiàn)出較為明顯的優(yōu)勢。而光球和色球數(shù)據(jù)的亞像素和超像素模擬位移實驗表明，Demons逐點測量精度在0.1像素量級。這意味著在1 m太陽望遠鏡18 s間隔的光球像上可以得到200 m/s精度的橫向速度，而在32 s的色球像上可以得到500 m/s的速度精度。另外，對于圖像結(jié)構(gòu)自相關(guān)性較大的圖像配準，這一方法無論在非剛性和剛性配準方面相對于傳統(tǒng)互相關(guān)方法有明顯的優(yōu)勢。

對比研究同時也表明，基于傅里葉變換的局部相關(guān)跟蹤法基本不適合測量高分辨像，而微分仿射速度估計直接用于強度像的測量也難以得到精確的結(jié)果。

Demons方法在亞像素和超像素測量上都能獲得較高的精度，但是由于其算法的復雜性，使得計算機開銷較大，計算時間較長，并不太適合大尺寸高分辨圖像的整體運算。同時處理的圖像有較為明顯的邊界效應，因此在實際處理時，需要對關(guān)注的區(qū)域做一定的邊緣擴展。

必須指出的是，這種利用光流場測量橫向速度場從本質(zhì)上是假設物質(zhì)的運動表現(xiàn)為圖像強度的橫向變化，而且這也僅僅是在觀測平面上的投影，并不是真正的物質(zhì)流動方向。同時如果有新圖像結(jié)構(gòu)的產(chǎn)生和舊結(jié)構(gòu)的消失，從本質(zhì)上這種測量方式是有缺陷的，因此需要高時間分辨率的觀測，以減少這種變化。但更短的時間間隔意味著更小的位移量。

理論上，這一方法可以用于處理其他望遠鏡的高分辨觀測圖像，也可以用于一些低分辨數(shù)據(jù)，或者其他類型的圖像上，這些還需要進一步嘗試。

致謝：感謝中國科學院云南天文臺澄江觀測站和1 m太陽望遠鏡全體人員對論文工作的大力支持和幫助。

參考文獻：

[1]Tsuneta S, Ichimoto K, Katsukawa Y, et al. The Solar Optical Telescope for the Hinode mission: an overview[J]. Solar Physics, 2008, 249(2): 167-196.

[2]Liu Z, Xu J, Gu B Z, et al. New vacuum solar telescope and observations with high resolution[J]. Research in Astronomy and Astrophysics, 2014, 14(6): 705-718.

[3]Goode P R, Denker C J, Didkovsky L I, et al. 1.6 M solar telescope in Big Bear-the NST[J]. Journal of the Korean Astronomical Society, 2003, 36(S1): S125-S133.

[4]Scharmer G B, Bjelksjo K, Korhonen T K, et al. The 1-meter Swedish solar telescope[C]// Proceedings of SPIE. 2002: 341-350.

[5]汪景琇, 季海生. 空間天氣驅(qū)動源—太陽風暴研究 [J]. 中國科學: 地球科學, 2013, 43(6): 883-911.

Wang Jingxiu, Ji Haisheng. Space weather driving Source-solar storm research[J]. Scientia Sinica Terrae, 2013, 43(6): 883-911.

[6]November L J, Simon G W. Precise proper-motion measurement of solar granulation[J]. The Astrophysical Journal, 1988, 333: 427-442.

[7]Fisher G H, Welsch B T. FLCT: a fast, efficient method for performing local correlation tracking[C]// Subsurface and Atmospheric Influences on Solar Activity ASP Conference Series. 2008: 383.

[8]Schuck P W. Local correlation tracking and the magnetic induction equation[J]. The Astrophysical Journal, 2005, 632(1): L53-L56.

[9]Schuck P W, Chen J. Tracking magnetic footpoints with the magnetic induction equation[J]. The Astrophysical Journal, 2006, 646(2): 1358-1391.

[10]Wang S, Liu C, Deng N, et al. Sudden photospheric motion and sunspot rotation associated with the X2.2 flare on 2011 February 15[J]. Astrophysical Journal Letters, 2014, 782(2): L31-L36.

[11]Liu C, Xu Y, Cao W, et al. Flare differentially rotates sunspot on Sun′s surface[J]. Nature Communications, 2016, 7: 13104.

[12]Potts H E, Barrett R K, Diver D A. Balltracking: an highly efficient method for tracking flow fields[J]. Astronomy & Astrophysics, 2004, 424(1): 253-262.

[13]Potts H E, Diver D A. Automatic recognition and characterisation of supergranular cells from photospheric velocity fields[J]. Solar Physics, 2008, 248(2): 263-275.

[14]Thirion J P. Image matching as a diffusion process: an analogy with Maxwell′s demons[J]. Medical Image Analysis, 1998, 2(3): 243-260.

[15]Kroon D J, Slump C H. MRI modality transformation in demon registration[C]// Proceedings of the IEEE International Symposium on Biomedical Imaging. 2009: 963-966.

[16]Wang Z, Bovik A C, Sheikh H R, et al. Image quality assessment: from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612.