射電天文圖像的反卷積算法研究?

張 利 徐 龍 米立功 馬家君

(1 貴州大學(xué)大數(shù)據(jù)與信息工程學(xué)院貴陽(yáng)550025)

(2 中國(guó)科學(xué)院太陽(yáng)活動(dòng)重點(diǎn)實(shí)驗(yàn)室北京100012)

(3 黔南民族師范學(xué)院物理與電子科學(xué)學(xué)院都勻558000)

1 引言

射電干涉測(cè)量通過(guò)若干小口徑望遠(yuǎn)鏡實(shí)現(xiàn)大口徑分辨率,并且在天空亮度分布的空間頻率域進(jìn)行采樣[1].由于干涉陣的天線數(shù)量有限,所以僅部分空間頻率被采樣,這使得望遠(yuǎn)鏡的點(diǎn)擴(kuò)展函數(shù)(也稱為臟束,dirty beam)具有不可忽略的旁瓣.這些旁瓣使得測(cè)量圖像模糊,從而限制成圖的動(dòng)態(tài)范圍.

在連續(xù)且完全采樣的情況下,由van Cittert-Zernike理論[1]可知,真實(shí)的天空亮度分布Itrue(x,y)與可見(jiàn)度函數(shù)(visibility)Vtrue(u,v)是一個(gè)傅里葉變換對(duì),即

其中j為虛數(shù)單位.然而實(shí)際測(cè)量是離散且不完全的,即對(duì)天空亮度分布的傅里葉空間進(jìn)行了不完全采樣,且存在噪聲Vnoise(u,v).測(cè)量的可見(jiàn)度函數(shù)Vmeasure(u,v)可表示為

其中S(u,v)為傅里葉空間采樣函數(shù),采樣點(diǎn)上的值為1,非采樣點(diǎn)上的值為0.通常定義臟圖(dirty image)Idirty(x,y)、點(diǎn)擴(kuò)展函數(shù)B(x,y)分別為

其中F?1表示傅里葉反變換.由卷積理論可知,傅里葉空間的乘積等價(jià)于像空間的卷積.則有

其中?表示卷積運(yùn)算,Inoise(x,y)=F?1(S(u,v)Vnoise(u,v)).因此,臟圖是設(shè)備點(diǎn)擴(kuò)展函數(shù)與天空亮度分布的卷積,同時(shí)包含噪聲的影響.點(diǎn)擴(kuò)展函數(shù)B(x,y)的旁瓣使得圖像模糊.反卷積的目的是去除旁瓣的影響,使圖像細(xì)節(jié)變得更加清晰.

在射電綜合成圖領(lǐng)域里,大致有3類反卷積:CLEAN算法[2],最大熵MEM算法[3?4]和壓縮感知算法[5?6].然而,CLEAN算法是最常用的,在射電干涉圖像處理軟件中是一個(gè)標(biāo)準(zhǔn)的組件.在1974年,由H?gbom提出的CLEAN算法[2]的主要目的是為了去除點(diǎn)擴(kuò)展函數(shù)的旁瓣影響.后來(lái)很多學(xué)者在不同的應(yīng)用場(chǎng)景下對(duì)算法進(jìn)行了改進(jìn).比較有代表性的是Clark CLEAN[7],S-C CLEAN[8],多尺度CLEAN[9?10]和自適應(yīng)尺度像素分解算法(Asp CLEAN)[11?13].Clark CLEAN算法使用了快速傅里葉變換和截?cái)嗟狞c(diǎn)擴(kuò)展函數(shù),加速了反卷積過(guò)程.在S-C CLEAN算法之前,CLEAN算法僅在像空間進(jìn)行反卷積,Schwab和Cotton將CLEAN算法擴(kuò)展到傅里葉空間,這種改進(jìn)有利于減少反卷積過(guò)程中的誤差累積[8].大多數(shù)的CLEAN算法將天空亮度分布分解為一系列delta函數(shù).對(duì)于分開(kāi)較好的致密源,這些方法的性能非常好,但是如果處理延展源,重建的圖像通常會(huì)帶有偽影(artefacts).為解決這個(gè)問(wèn)題,Cornwell[9]在CLEAN算法中引入多尺度基函數(shù).將天空亮度分布分解成一系列尺度基函數(shù),這些基函數(shù)的尺度也包括0,即包括delta函數(shù).所以多尺度CLEAN算法對(duì)于致密源和延展源均有效.然而,尺度的大小需要用戶指定且尺度數(shù)目受限于計(jì)算機(jī)內(nèi)存.為解決這個(gè)問(wèn)題,Bhatnagar和Cornwell提出了自適應(yīng)尺度像素分解算法[11].這種算法通過(guò)擬合來(lái)實(shí)現(xiàn)基函數(shù)尺度的自適應(yīng),從而能夠重建出更高質(zhì)量的圖像.但這種算法使用了有效集的優(yōu)化方式[11],使得反卷積的時(shí)間很長(zhǎng).文獻(xiàn)[12]通過(guò)引入解析高斯反卷積對(duì)這個(gè)算法進(jìn)行了優(yōu)化和加速.這種改進(jìn)方法在反卷積過(guò)程中需要將帶有旁瓣的點(diǎn)擴(kuò)展函數(shù)逼近成一個(gè)高斯函數(shù),所以僅適用于點(diǎn)擴(kuò)展函數(shù)旁瓣較小的情形.針對(duì)這些問(wèn)題,本文提出一種基于尺度基函數(shù)的CLEAN算法.

2 基于尺度基函數(shù)的反卷積算法

在CLEAN類反卷積算法中,每次迭代通過(guò)找到最可能的真實(shí)分量來(lái)不斷地逼近真實(shí)的天空亮度分布.在本文算法中,通過(guò)一系列高斯基函數(shù)來(lái)逼近天空亮度分布,

其中K為組成天空亮度分布的分量數(shù)目,?為逼近的天空亮度分布和原始天空亮度分布之間的誤差,ai為第i個(gè)分量的幅度,(xi,yi)為第i個(gè)分量的中心位置,(rix,riy)為第i個(gè)分量的x和y方向的寬度.本算法采用如下方法找到各個(gè)最優(yōu)分量.

(1)找到初始的參數(shù)(αi0,xi0,yi0,rix0,riy0).首先使用幾個(gè)高斯函數(shù)平滑殘差圖像(第1次迭代時(shí)是臟圖),然后從前一步得到的平滑殘差圖像中找到全局極大值,將其對(duì)應(yīng)的高斯函數(shù)的寬度作為高斯分量的寬度(半高全寬).

其中,Gk為1個(gè)Toeplitz矩陣,每行元素組成1維形式的高斯函數(shù),下一行是通過(guò)循環(huán)位移上一行的一個(gè)元素得到的,k是平滑的殘差圖像的數(shù)量,是第i次迭代的殘差圖像.

平滑時(shí)一般選取4個(gè)高斯函數(shù)即可.如果數(shù)目太少,則無(wú)法找到好的初始值.如果數(shù)目太多,會(huì)顯著增加計(jì)算量.另外,為了使描述更加簡(jiǎn)潔,在表示圖像的符號(hào)中,帶坐標(biāo)的表示圖像矩陣,如(x,y);不帶坐標(biāo)的表示圖像向量,如.兩者表示相同的圖像,只是表示不同.

(2)通過(guò)Levenberg-Marquardt最小化算法[14]優(yōu)化初始分量參數(shù)來(lái)獲得最優(yōu)分量參數(shù)(αi,xi,yi,rix,riy).在這里,需要最小化目標(biāo)函數(shù)

其中為第i個(gè)高斯分量.

(3)更新模型.在計(jì)算最優(yōu)分量時(shí),對(duì)高斯寬度乘以循環(huán)增益g.然后將最優(yōu)分量加入到模型圖像中.

其中為高斯寬度乘以循環(huán)增益g后的分量.循環(huán)增益可以優(yōu)化模型分量,防止過(guò)分估計(jì),實(shí)現(xiàn)更好的重建.正如其他CLEAN算法一樣,需要用戶指定其大小.對(duì)于尺度敏感CLEAN算法,一般g<0.5.

(4)計(jì)算殘差.

重復(fù)上面步驟,直到最大迭代次數(shù)或噪聲水平時(shí)停止.正如其他CLEAN算法一樣,最后的重建圖像ICLEAN等于潔束(CLEAN beam,通常是高斯函數(shù))與模型圖像Imodel的卷積與殘差I(lǐng)residual之和,即

其中Imodel,Iresidual分別是最后的模型圖像和殘差圖像.BCLEAN是一個(gè)Toeplitz矩陣,其行元素由潔束的離散值組成.

在初始參數(shù)計(jì)算階段,需要通過(guò)平滑當(dāng)前殘差圖像來(lái)找到一個(gè)合適的初始參數(shù).在自適應(yīng)尺度像素分解算法[11]中,這種平滑通過(guò)卷積實(shí)現(xiàn),是一個(gè)比較耗時(shí)的部分.在本文的算法中,當(dāng)分量尺度較大時(shí)使用卷積平滑的方式來(lái)找到初始參數(shù),當(dāng)分量尺度較小時(shí)直接使用殘差圖像中絕對(duì)值最大的點(diǎn)作為初始位置參數(shù)和幅度參數(shù),并使用點(diǎn)擴(kuò)展函數(shù)主瓣寬度作為初始寬度參數(shù).本算法中,當(dāng)最近連續(xù)10次的分量寬度平均值在x和y方向同時(shí)小于點(diǎn)擴(kuò)展函數(shù)的主瓣寬度時(shí),即

下一次計(jì)算初始參數(shù)時(shí),適用分量尺度較小的方法.實(shí)驗(yàn)證明它能夠有效地估計(jì)較小的初始分量,從而避免更多的計(jì)算量.在測(cè)試中發(fā)現(xiàn),在反卷積過(guò)程中,小尺度的分量占據(jù)較大比例,這種方法加速了反卷積過(guò)程.在本算法中,首先找到初始的分量參數(shù),然后通過(guò)擬合的方式對(duì)初始參數(shù)進(jìn)行優(yōu)化.在自適應(yīng)尺度像素分解算法[11]中,使用有效集的方式來(lái)實(shí)現(xiàn)分量的更好擬合.盡管只有少部分的分量在有效集中,然而參數(shù)的成倍增加使得擬合時(shí)間大大增加,同時(shí)也增加擬合失敗的風(fēng)險(xiǎn).在本文算法中,并不是使用有效集的方式,即前面的分量不被再次擬合優(yōu)化.每次擬合參數(shù)并不增加,這顯著加速了反卷積過(guò)程.本質(zhì)上,文獻(xiàn)[11]追求搜索空間的正交化,從而實(shí)現(xiàn)對(duì)源更加稀疏的表示.然而這需要付出更大的時(shí)間代價(jià).本文算法并未追求搜索空間的正交化,同時(shí)使用新的擬合參數(shù)計(jì)算方法.兩者同時(shí)使得本文算法的速度明顯提升.在文獻(xiàn)[12]中,反卷積過(guò)程中點(diǎn)擴(kuò)展函數(shù)是被近似成一個(gè)高斯函數(shù),這會(huì)導(dǎo)致分量估計(jì)不準(zhǔn)確.當(dāng)點(diǎn)擴(kuò)展函數(shù)的旁瓣水平越高時(shí),分量估計(jì)誤差越明顯.然而本文的方法使用完整的點(diǎn)擴(kuò)展函數(shù)來(lái)進(jìn)行反卷積,所以分量估計(jì)更加準(zhǔn)確.Levenberg-Marquardt算法結(jié)合了高斯-牛頓算法和梯度下降算法,具有魯棒的性能,有利于分量的準(zhǔn)確擬合.

CLEAN類算法通過(guò)迭代方法求解線性等式系統(tǒng)[15?16],文獻(xiàn)[1,15]已經(jīng)很好地證明了CLEAN類算法的收斂性,其收斂條件為:(1)點(diǎn)擴(kuò)展函數(shù)B(矩陣形式下)必須是對(duì)稱的;(2)B必須是正定或半正定的;(3)臟圖Idirty必須在B的值域(range)內(nèi).本文算法使用了不同的分量提取方法,但算法的結(jié)構(gòu)與CLEAN算法相同,同時(shí)保持了CLEAN算法的最速下降思想,因此文獻(xiàn)[15]的收斂證明適用于本文算法.具體細(xì)節(jié)請(qǐng)參考文獻(xiàn)[15].

3 算法實(shí)驗(yàn)

本文利用天文通用軟件包(Common Astronomy Software Applications,CASA)和Python語(yǔ)言實(shí)現(xiàn)了上述算法.為了測(cè)試算法的反卷積性能和避免其他因素影響反卷積的測(cè)試效果,我們使用CASA模擬The Karl G.Jansky Very Large Array(JVLA)B陣型對(duì)圖像M51(如圖1(a))進(jìn)行觀測(cè),獲得模擬觀測(cè)數(shù)據(jù).觀測(cè)波段為L(zhǎng)波段,1 GHz帶寬和32個(gè)通道.我們將模擬的可見(jiàn)度函數(shù)加入高斯噪聲,使得臟圖(圖1(b))中噪聲的均方根誤差為5×10?5Jy.在測(cè)試中,使用魯棒加權(quán)對(duì)測(cè)量數(shù)據(jù)進(jìn)行柵格化.圖1和圖2中的單位均為央斯基每像素(Jy/pixel).從圖中可以看出,臟圖中天體的細(xì)節(jié)變得模糊不清.圖1(c)為本文算法從臟圖中重建的圖像,可以看出已經(jīng)很好地重建出原始圖像中細(xì)節(jié),包括致密源和延展源.圖2(a0)和(a1)分別顯示了H?gbom CLEAN(Hg-Clean)算法[2]的模型圖像和殘差圖像,可以看出該算法對(duì)延展結(jié)構(gòu)的表示不佳和殘差圖像中有明顯延展信號(hào).這源于delta函數(shù)無(wú)法很好地逼近延展結(jié)構(gòu).從圖2可以看出,多尺度CLEAN(Ms-Clean)算法[9]的結(jié)果優(yōu)于Hg-Clean算法的結(jié)果.這是因?yàn)槎喑叨菴LEAN算法使用尺度基函數(shù)來(lái)逼近潛在的真實(shí)圖像.然而它的尺度數(shù)量是有限的,一些延展結(jié)構(gòu)不可避免地被分解成預(yù)先設(shè)定的尺度.自適應(yīng)尺度像素分解算法[11]很好地解決了尺度的自適應(yīng)的問(wèn)題,相對(duì)于多尺度CLEAN算法具有更好的性能.圖2(d1)顯示的是殘差圖像,可以看出殘差圖像中沒(méi)有明顯的信號(hào),說(shuō)明本文算法能夠很好地分離源和噪聲.這得益于Levenberg-Marquardt算法很好的擬合性能.同時(shí)從圖2可以看出,自適應(yīng)尺度像素分解算法[11]的重建結(jié)果與本文算法類似,但是從表1和表2可以發(fā)現(xiàn)本文算法的速度更快.在測(cè)試中,我們發(fā)現(xiàn)循環(huán)增益使用0.2–0.7,可以獲得不錯(cuò)的性能.本文算法中的循環(huán)增益取值明顯高于基于delta函數(shù)分解的CLEAN算法中循環(huán)增益的取值,這是因?yàn)楸疚乃惴ㄊ褂昧藬M合技術(shù)對(duì)初始參數(shù)進(jìn)行了再一次優(yōu)化,而基于delta函數(shù)分解的CLEAN算法中并沒(méi)有進(jìn)行再優(yōu)化的過(guò)程.

本文算法的優(yōu)點(diǎn)是,能夠重建出高質(zhì)量圖像的同時(shí)提高了反卷積速度.為了比較不同加權(quán)情況下的反卷積速度,分別使用自然加權(quán)、均勻加權(quán)和魯棒加權(quán)計(jì)算臟圖.圖像尺寸為256×256像素.在空間頻率域使用不同的加權(quán)機(jī)制會(huì)得到不同的點(diǎn)擴(kuò)展函數(shù),同時(shí)得到不一樣的臟圖.一般來(lái)說(shuō),均勻加權(quán)得到的點(diǎn)擴(kuò)展函數(shù)旁瓣較少,自然加權(quán)得到的點(diǎn)擴(kuò)展函數(shù)旁瓣較高,而魯棒加權(quán)得到的點(diǎn)擴(kuò)展函數(shù)的旁瓣水平介于兩者之間.在當(dāng)前比較典型的普通計(jì)算機(jī)(Intel(R)Core(TM)i7-3770 CPU@3.4 GHz,4.00 GB RAM)上,本文算法(Python實(shí)現(xiàn))處理256×256大小的圖像大概在幾分鐘完成.正如其他算法一樣,參數(shù)不一樣時(shí),迭代次數(shù)不一樣,因而反卷積所用的時(shí)間也會(huì)不一樣.

圖1 本文算法處理M51的結(jié)果.(a)真實(shí)圖像;(b)使用魯棒加權(quán)形成的臟圖;(c)重建的模型圖像;(d)對(duì)應(yīng)的重建圖像.Fig.1The results of M51 processed by our algorithm.(a)The true image;(b)the dirty image with robust weighting function;(c)the reconstructed model image;(d)the corresponding restored image.

圖2 典型CLEAN算法的結(jié)果比較.(a0)Hg-Clean算法[2]的模型圖像;(a1)Hg-Clean算法的殘差圖像;(b0)多尺度CLEAN算法[9]的模型圖像;(b1)多尺度CLEAN算法的殘差圖像;(c0)Asp-Clean算法[11]的模型圖像;(c1)Asp-Clean算法的殘差圖像;(d0)本文算法的模型圖像;(d1)本文算法的殘差圖像.Fig.2 The results compared with that of typical CLEAN algorithm.(a0)The modeled image from the Hg-Clean algorithm[2];(a1)the residual image from the Hg-Clean algorithm;(b0)the modeled image from the Ms-Clean algorithm[9];(b1)the residual image from the Ms-Clean algorithm;(c0)the modeled image from the Asp-Clean algorithm[11];(c1)the residual image from the Asp-Clean algorithm;(d0)the modeled image from our algorithm;(d1)the residual image from our algorithm.

表1中比較了不同加權(quán)機(jī)制下的性能.從表中可以發(fā)現(xiàn):(1)本文算法的分量數(shù)目多一些.這是由于自適應(yīng)尺度像素分解算法采用了有效集,部分分量被多次優(yōu)化,有利于減少分量數(shù)目.(2)單次迭代所花時(shí)間縮短了3倍以上,同時(shí)單次迭代所花的時(shí)間也有較大的差異.這種差異部分來(lái)源于自適應(yīng)尺度像素分解算法中的有效集大小的不同.當(dāng)其他參數(shù)確定時(shí),有效集的大小取決于點(diǎn)擴(kuò)展函數(shù)的旁瓣水平.(3)整體上,本文算法速度提高了3倍左右.

表1 在不同加權(quán)機(jī)制下的性能比較Table 1 Performances with different weighting functions

在不同圖像大小情況下(其他參數(shù)一致,并使用均勻加權(quán)),對(duì)比了本文算法與自適應(yīng)尺度算法之間的反卷積速度.從表2可以看出,對(duì)于256×256,512×512,1024×1024,2048×2048的圖像大小,本文算法的反卷積速度均有提升.另外,速度的提升并不是隨圖像的大小而線性增加的,這可能是由于CLEAN反卷積過(guò)程本身是一個(gè)非線性過(guò)程導(dǎo)致的.

表2 在不同圖像大小情況下的速度比較Table 2 Runtime with different image sizes

4 總結(jié)與討論

自適應(yīng)尺度像素分解算法在射電綜合成圖中可以獲得當(dāng)前最好的重建效果,然而其計(jì)算時(shí)間卻很長(zhǎng).本文針對(duì)這個(gè)問(wèn)題,在自適應(yīng)尺度像素分解算法基礎(chǔ)上,使用不同的初始值計(jì)算方法和擬合策略.實(shí)驗(yàn)表明,本算法在獲得高質(zhì)量重建圖像的同時(shí),反卷積速度有較為明顯的提升.

反卷積研究的是如何有效地找到分量來(lái)逼近真實(shí)的天空亮度分布,是射電綜合成圖中的關(guān)鍵技術(shù).大量研究給出了不同的反卷積方法,然而這些方法往往適用于某一場(chǎng)景.比如,Hg-Clean算法適用于很好分開(kāi)的致密源場(chǎng)景,而自適應(yīng)尺度像素分解算法對(duì)延展源更有效.盡管自適應(yīng)尺度像素分解算法解決了尺度自適應(yīng)問(wèn)題,提高了重建圖像的質(zhì)量,然而導(dǎo)致了時(shí)間開(kāi)銷顯著增加[11].本文提出了一種保持與文獻(xiàn)[11]具有同等重建質(zhì)量且有速度提升的算法,應(yīng)用場(chǎng)景與文獻(xiàn)[11]一致.然而,針對(duì)不同的應(yīng)用場(chǎng)景,當(dāng)前并沒(méi)有相對(duì)普適的算法,所以提出更有效的方法逼近各種場(chǎng)景的天空亮度分布仍然是下一步的研究方向.