張 利，米立功，肖一凡，楊江河，李桂城,衛(wèi)星奇，寧娜文，李丹楊，賀春林

(1.貴州大學(xué)大數(shù)據(jù)與信息工程學(xué)院，貴州 貴陽(yáng) 550025；2. 黔南民族師范學(xué)院物理與電子科學(xué)學(xué)院，貴州 都勻 558000；3. 西華師范大學(xué)計(jì)算機(jī)學(xué)院，四川 南充 637002；4. 湖南文理學(xué)院 湖南 常德 415000)

射電干涉測(cè)量通過(guò)干涉技術(shù)能夠獲得等價(jià)于最長(zhǎng)基線長(zhǎng)度的陣列口徑，這種測(cè)量技術(shù)大大提高了望遠(yuǎn)鏡的分辨率.世界已經(jīng)有很多著名的干涉陣，如甚大天線陣(VLA)、阿塔卡馬大型毫米波/亞毫米波天線陣(ALMA)等.本世紀(jì)初以來(lái)，中國(guó)也加大對(duì)干涉陣的建設(shè)和資助力度，如已經(jīng)運(yùn)行的“宇宙第一縷曙光”探測(cè)項(xiàng)目(21CMA)和建設(shè)階段的國(guó)際合作大項(xiàng)目——平方公里射電陣(SKA)等.這些望遠(yuǎn)鏡陣都有很好的分辨率，能夠有效分辨遙遠(yuǎn)的天體.然而很多干涉陣的點(diǎn)擴(kuò)展函數(shù)(Point Spread Function, PSF)具有較高的旁瓣，因此移去點(diǎn)擴(kuò)展函數(shù)的模型影響是有必要的.即使像SKA這種望遠(yuǎn)鏡數(shù)量眾多且設(shè)計(jì)優(yōu)秀的望遠(yuǎn)鏡陣，如果想要得到一個(gè)高動(dòng)態(tài)范圍的圖像，反卷積重建仍然是必要的.

1 基本問(wèn)題

Van Citter-Zernike理論[1]指出，可見(jiàn)度函數(shù)Vt與天體圖像It存在傅里葉變換關(guān)系.由于測(cè)量是在空間頻率域進(jìn)行的，因此噪聲也在空間頻率域.測(cè)量的可見(jiàn)度函數(shù)Vm可以表示成

Vm=M(Vt+n0).

(1)

其中，n0為空間頻率域噪聲，M為采樣函數(shù).M屬于非笛卡爾采樣，在數(shù)據(jù)處理中，會(huì)重新采樣到笛卡爾坐標(biāo)系中.它的值為0或1，0代表未采樣的點(diǎn)，1代表采樣的點(diǎn).由于望遠(yuǎn)鏡的地理位置、觀測(cè)源的位置和觀測(cè)時(shí)間等一系列原因，一般情況下，會(huì)有大量的點(diǎn)測(cè)量不到.數(shù)據(jù)重建就是通過(guò)一些方法來(lái)估計(jì)出這些測(cè)量不到的空間頻率.在空間域，測(cè)量的可見(jiàn)度函數(shù)Vm和采樣函數(shù)M分別對(duì)應(yīng)臟圖Id和點(diǎn)擴(kuò)展函數(shù)(臟束)B的傅里葉反變換.由卷積理論知道，臟圖是點(diǎn)擴(kuò)展函數(shù)與天體圖像和噪聲的卷積，即

Id=It*B+n*B，

(2)

式中，*表示卷積運(yùn)算.未采樣的頻率在空間上對(duì)應(yīng)“不可見(jiàn)分布”[1].如果p對(duì)應(yīng)采樣頻率所對(duì)應(yīng)的分布，那么p*B=0.如果Ii是卷積等式(2)的解，那么Ii+kp也是這個(gè)方程的解.因?yàn)閗是任意的，所以有無(wú)數(shù)多個(gè)解滿足卷積方程(2).數(shù)據(jù)重建的目的是找到一個(gè)最合理的“不可見(jiàn)分布”.

2 常用算法及其問(wèn)題

目前有許多用于射電干涉圖像重建的算法，如最大熵算法[2-3]、CLEAN算法[4-8]、壓縮感知重建算法[9-11]等.其中一些基于函數(shù)模型分解的算法被廣泛使用，如CLEAN算法.基于函數(shù)模型分解的算法主要有兩大類(lèi)：基于無(wú)尺度(Delta)函數(shù)模型分解的算法[4-8]和基于尺度函數(shù)模型分解的算法[12-13].基于Delta函數(shù)模型分解的算法擅長(zhǎng)于處理致密源，而基于尺度函數(shù)模型分解的算法則更擅長(zhǎng)于處理延展源.基于尺度函數(shù)模型分解的算法也能處理致密源，所以這類(lèi)算法的性能常常要優(yōu)于基于Delta函數(shù)模型分解的算法.非常典型的基于尺度函數(shù)模型分解的算法有2種：多尺度CLEAN算法[12]和自適應(yīng)尺度像素分解算法[13].多尺度CLEAN算法基于球體波函數(shù)的分解，使用匹配濾波技術(shù)來(lái)尋找一個(gè)合適的分量.由于技術(shù)的限制，函數(shù)的尺度僅有幾個(gè)，且這幾個(gè)尺度是用戶指定的，在重建過(guò)程中無(wú)法修改，所有的分量將被分解到這幾個(gè)固定的分量；因此圖像分量的尺度是無(wú)法確定的，需要更多的尺度來(lái)逼近潛在的真實(shí)圖像.同時(shí)，多尺度CLEAN算法有預(yù)計(jì)算過(guò)程，預(yù)計(jì)算有較高的內(nèi)存要求，如果圖像過(guò)大或使用太多的尺度，可能會(huì)因內(nèi)存不足而造成計(jì)算失敗.自適應(yīng)尺度像素分解算法是基于高斯函數(shù)的分解，使用非線性最優(yōu)方法來(lái)實(shí)現(xiàn)尺度的自適應(yīng)，很好地解決了多尺度CLEAN算法因固定尺度所帶來(lái)的問(wèn)題，有非常好的圖像重建效果.同時(shí)，自適應(yīng)尺度像素分解算法追求更加稀疏的分解，所以重復(fù)使用非線性最優(yōu)方法來(lái)優(yōu)化前面所得到的參數(shù).然而盡管使用了加速的方法，這種算法的計(jì)算量還是非常大，以至于沒(méi)有寫(xiě)入軟件中.因此，設(shè)計(jì)一種新的算法以平衡計(jì)算的復(fù)雜度和性能，具有重要的意義.

3 自適應(yīng)尺度重建算法

為了解決計(jì)算復(fù)雜度和性能的平衡問(wèn)題，筆者提出了一種新的反卷積重建算法.算法基于自適應(yīng)尺度高斯函數(shù)的分解，使用Levenberg-Marquardt非線性最優(yōu)方法[14]來(lái)實(shí)現(xiàn)尺度的自適應(yīng)能力.這種二階的非線性最優(yōu)方法也保證了模型分量的最優(yōu).對(duì)于每一個(gè)模型分量，算法需要最優(yōu)化目標(biāo)函數(shù)

(3)

(4)

其中，f(x,y;p)是高斯函數(shù)，g是循環(huán)增益.該算法中截?cái)嗟母咚购瘮?shù)被使用，這有效地解決了高斯函數(shù)支撐集太大的問(wèn)題.算法的基本過(guò)程為：

Step 1 匹配濾波找到一個(gè)合適的尺度和位置；

Step 2 使用非線性優(yōu)化方法進(jìn)行高斯函數(shù)擬合，找到最優(yōu)解；

Step 3 更新最優(yōu)解到模型；

Step 4 計(jì)算殘差；

Step 5 滿足最大迭代次數(shù)或達(dá)到噪聲水平時(shí)停止.

自適應(yīng)尺度重建算法被設(shè)計(jì)用于柵格化的數(shù)據(jù)，然而測(cè)量的射電干涉數(shù)據(jù)并不在柵格點(diǎn)上，需要很多的運(yùn)算來(lái)進(jìn)行數(shù)據(jù)柵格化.因此結(jié)合CLEAN框架，通過(guò)使用CLEAN算法的主循環(huán)(Major Cycle)，有效利用成圖(Imaging)部分的算法，如卷積插值和重采樣等，使算法能處理非柵格化的直接測(cè)量數(shù)據(jù)，也能有效地糾正一些小循環(huán)中產(chǎn)生的誤差[1].該算法可以在CLEAN框架里作為次循環(huán)(Minor Cycle)，主要為了找到最優(yōu)的模型，具體流程如圖1所示.

圖1 基于自適應(yīng)尺度重建算法重建圖像的流程Fig. 1 Diagram of Reconstruction Algorithm Based on Adaptive Scale

與自適應(yīng)尺度像素分解算法一樣，自適應(yīng)尺度重建算法也使用了非線性最優(yōu)方法來(lái)優(yōu)化模型分量；所不同的是，自適應(yīng)尺度像素分解算法重復(fù)地優(yōu)化前面已經(jīng)優(yōu)化的參數(shù)，而自適應(yīng)尺度重建算法不再重復(fù)優(yōu)化前面已得到的參數(shù)，大大減少了算法的計(jì)算復(fù)雜度.相比于多尺度CLEAN算法，自適應(yīng)尺度重建算法的尺度是自適應(yīng)的，能夠有效解決固定尺度分解所帶來(lái)的問(wèn)題.同時(shí)，該算法沒(méi)有多尺度CLEAN算法的預(yù)計(jì)算步驟，可有效地解決多尺度CLEAN算法的內(nèi)存開(kāi)銷(xiāo)問(wèn)題.

4 數(shù)值模擬

自適應(yīng)尺度重建算法可應(yīng)用于非柵格化數(shù)據(jù).數(shù)據(jù)利用通用天文軟件程序包(Common Astronomy Software Applications, CASA)來(lái)模擬，使用VLA的B陣型在L波段模擬觀測(cè)到的射電源M31，模擬的觀測(cè)數(shù)據(jù)格式和真實(shí)的觀測(cè)數(shù)據(jù)格式是相同的，測(cè)量的量也是一致的.相對(duì)于真實(shí)數(shù)據(jù)，模擬得到的數(shù)據(jù)更加完善，模擬的數(shù)據(jù)相當(dāng)于校準(zhǔn)之后的真實(shí)觀測(cè)數(shù)據(jù).對(duì)于測(cè)試圖像重建算法而言，與使用真實(shí)數(shù)據(jù)的函數(shù)有同樣的效果.筆者使用殘差圖像的均方根誤差(Root Mean Square Error, RMSE)和動(dòng)態(tài)范圍(Dynamic Range, DR)[15]比較了自適應(yīng)尺度重建算法和H?gbom CLEAN算法及多尺度CLEAN算法.H?gbom CLEAN算法和多尺度CLEAN算法是在射電干涉天文圖像處理中常用的兩種算法，分別是基于Delta函數(shù)分解的反卷積算法代表和基于多尺度函數(shù)分解的反卷積算法代表.

殘差圖像的均方根誤差能夠有效地表示殘差的水平，均方根誤差越小，說(shuō)明殘差水平越低，圖像中的信息被恢復(fù)越多.殘差圖像的均方根誤差

(5)

其中，N是圖像的像素點(diǎn)數(shù)，‖·‖ln是ln范數(shù).

動(dòng)態(tài)范圍是天文圖像處理中常用的判斷圖像質(zhì)量的依據(jù)，越高的動(dòng)態(tài)范圍意味著圖像更好地被重建.動(dòng)態(tài)范圍被定義為重建圖像的最大值與殘差圖像的非源區(qū)域的均方根誤差之比，即

(6)

其中，Ire表示重建的圖像，RMSEp表示殘差圖像的非源區(qū)域的均方根誤差.

本數(shù)值模擬實(shí)驗(yàn)所處理的射電源M31圖像原圖及其點(diǎn)擴(kuò)展函數(shù)圖像和臟圖如圖2所示.

圖2 自適應(yīng)尺度重建算法所處理的M31圖像Fig. 2 M31 Image Processed by Improved Algorithm

從圖2(a)可以看出，射電源M31是一個(gè)非常復(fù)雜的源，包含了大量子結(jié)構(gòu)，能有效地測(cè)試各重建算法的性能.采集到的數(shù)據(jù)表明，射電源M31點(diǎn)擴(kuò)展函數(shù)(圖2(b))的最大負(fù)旁瓣是-0.068 7.從圖2(b)可以看出，旁瓣結(jié)果復(fù)雜且分布很廣.被點(diǎn)擴(kuò)展函數(shù)模糊化之后的臟圖為圖2(c).測(cè)試圖像原圖的大小為512×512像素，為了使顯示更加清晰便于比較，成圖僅顯示處理結(jié)果的中央256×256像素部分.

在成圖過(guò)程中所有算法均使用Robust加權(quán)函數(shù).Robust加權(quán)函數(shù)能夠有效抑制旁瓣水平，是一種常用的加權(quán)函數(shù).H?gbom CLEAN算法使用的參數(shù)是：迭代次數(shù)30 000，循環(huán)增益0.1.多尺度CLEAN算法使用的參數(shù)是：迭代次數(shù)10 000，循環(huán)增益0.3，尺度分別為0，4，12，36，108像素.自適應(yīng)尺度重建算法使用的參數(shù)是：迭代次數(shù)7 000，循環(huán)增益0.7.由于自適應(yīng)尺度重建算法是基于自適應(yīng)尺度的，所以不需要用戶指定尺度.各算法的處理結(jié)果如圖3所示.

圖3 典型CLEAN算法與自適應(yīng)尺度重建算法的結(jié)果對(duì)比Fig. 3 Comparison of Results from Typical CLEAN Algorithm and the Improved Algorithm

從圖3(a)中可以看出，H?gbom CLEAN算法重建的模型圖像由點(diǎn)組成，這是由于H?gbom CLEAN算法是基于Delta函數(shù)分解的，即圖像被分解為Delta函數(shù)集合.Delta函數(shù)集合是通過(guò)找到殘差圖像中的最大絕對(duì)值點(diǎn)，然后乘以循環(huán)增益得到Delta函數(shù)的位置和系數(shù).這種分解不能有效地表示延展結(jié)構(gòu)，所以使用這類(lèi)算法處理延展源時(shí)，殘差圖像中會(huì)殘留“平臺(tái)”式的相關(guān)結(jié)構(gòu)，而且需要大量的分量來(lái)表示圖像.從這個(gè)測(cè)試中可以看出，盡管H?gbom CLEAN算法使用了遠(yuǎn)大于其他兩種基于尺度函數(shù)分解的算法的迭代次數(shù)，然而結(jié)果卻是最差的.多尺度CLEAN算法使用了幾個(gè)尺度函數(shù)來(lái)表示圖像，較好地解決了基于Delta函數(shù)模型的算法不能有效表示圖像的問(wèn)題.從圖3(b)也能看出，相對(duì)于H?gbom CLEAN算法，多尺度CLEAN算法重建的模型圖像更加接近于真實(shí)圖像，殘差圖像中所留下的相關(guān)結(jié)構(gòu)也更小(圖3(e)).然而基于多尺度函數(shù)模型分解的算法將圖像強(qiáng)行分解到僅有的尺度當(dāng)中，這需要更多的分量來(lái)表示圖像.基于自適應(yīng)尺度函數(shù)模型的算法能有效地解決這一問(wèn)題，使用更少的分量來(lái)表示圖像.從圖3(c)中可以看出，相對(duì)于H?gbom CLEAN算法和多尺度CLEAN算法，筆者提出的自適應(yīng)尺度重建算法模型圖像更加接近于真實(shí)圖像.殘差圖像(圖3(f))也能有效證明這一結(jié)論.圖3還表明，在這三種算法中，自適應(yīng)尺度重建算法的殘差相關(guān)結(jié)構(gòu)是最少的，說(shuō)明自適應(yīng)尺度函數(shù)模型能夠更好地表示圖像.表1列出三種算法重建圖像的RMSE和DR數(shù)值結(jié)果.

表1 三種方法重建結(jié)果的數(shù)值比較

從表1可以看出，在這三種算法中，自適應(yīng)尺度重建算法有最小的殘差均方根誤差和最高的動(dòng)態(tài)范圍，進(jìn)一步說(shuō)明自適應(yīng)尺度重建算法的性能優(yōu)于基于Delta函數(shù)模型算法和基于多尺度函數(shù)模型算法.

5 結(jié)語(yǔ)

多尺度CLEAN算法是一種基于多尺度函數(shù)模型的重建算法，廣泛應(yīng)用于射電干涉天文圖像的重建，但是由于尺度是固定的，內(nèi)存消耗嚴(yán)重，因此不宜處理大圖像.自適應(yīng)尺度像素分解算法利用非線性最優(yōu)化使尺度具有自適應(yīng)能力，解決了固定尺度不能有效表示延展結(jié)構(gòu)問(wèn)題，圖像質(zhì)量較多尺度CLEAN算法有很大提升.然而由于該算法重復(fù)優(yōu)化前面已經(jīng)優(yōu)化了的參數(shù)，盡管每一分量更加準(zhǔn)確，但是計(jì)算量卻大大增加，以致無(wú)法寫(xiě)入軟件中.筆者提出了一種改進(jìn)的自適應(yīng)尺度重建算法，不再重復(fù)地優(yōu)化前面已經(jīng)優(yōu)化了的參數(shù)，后面的優(yōu)化會(huì)糾正前面優(yōu)化的不準(zhǔn)確性.自適應(yīng)尺度重建算法是基于CASA和Python編程實(shí)現(xiàn)的，為了能夠處理非柵格化數(shù)據(jù)，將自適應(yīng)尺度重建算法與CLEAN框架相結(jié)合，計(jì)算的復(fù)雜度較自適應(yīng)尺度像素分解算法大大減少，重建結(jié)果優(yōu)于基于多尺度函數(shù)模型算法.