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

(1. 昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院，云南 昆明 650500; 2. 昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650500)

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

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

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

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

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

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

1　Demons方法

1.1　Demons 配準(zhǔn)的原理和流程

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

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

(4)

(4)式可以簡化為

v·r=f-r,

(5)

(6)

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

(7)

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

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

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

(8)

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

基于圖像灰度互信息的配準(zhǔn)方法，最大的不足是忽略了圖像的空間信息。就太陽圖像配準(zhǔn)過程而言，對相同目標(biāo)所成圖像可能具有不同的灰度、不同的分辨率，甚至圖像本身大小也不同，但是對于同一目標(biāo)，其邊界是確定的，不會(huì)發(fā)生很明顯的變化；當(dāng)兩幅圖像配準(zhǔn)時(shí)，對應(yīng)位置像素點(diǎn)的梯度矢量將具有相同或者相反的指向?；谝陨戏治觯荻然バ畔⒉捎肞arzen窗方法進(jìn)行計(jì)算。

2　實(shí)驗(yàn)數(shù)據(jù)

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

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

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

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

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

圖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個(gè)數(shù)據(jù)集采用Demons方法計(jì)算第2幀相對于第1幀的逐像素位移量，即光流場，然后根據(jù)相應(yīng)的時(shí)間間隔歸算到逐點(diǎn)的橫向速度場。圖2顯示了數(shù)據(jù)集1中一個(gè)小區(qū)域的兩幀圖以及用Demons方法計(jì)算出的位移矢量。米粒的細(xì)微運(yùn)動(dòng)可以清晰地反映在圖上。

圖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)跟蹤法和微分仿射速度估計(jì)方法進(jìn)行比對，在分別采用不同方法得到光流場后，采用雙線性插值將第2幀與第1幀進(jìn)行非剛性配準(zhǔn)，并利用圖像結(jié)構(gòu)相似度指數(shù)(Structural Similarity Index Measurement, SSIM)[16]測試配準(zhǔn)的精度。較高的配準(zhǔn)精度即意味著較高的光流場(速度場)測量精度。

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

(9)

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

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

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

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

Table1TheSSIMindexafternon-rigidregistrationwithdifferentmethod

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

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

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

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

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

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

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

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

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

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

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

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

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

5　結(jié)　論

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

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

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

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

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

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

參考文獻(xiàn)：

[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ū)動(dòng)源—太陽風(fēng)暴研究 [J]. 中國科學(xué): 地球科學(xué), 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.