張 利 徐 龍 米立功 馬家君

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

(2 中國科學院太陽活動重點實驗室北京100012)

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中為第i個高斯分量.

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

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

(4)計算殘差.

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

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

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

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

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

3 算法實驗

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

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

圖1 本文算法處理M51的結果.(a)真實圖像;(b)使用魯棒加權形成的臟圖;(c)重建的模型圖像;(d)對應的重建圖像.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算法的結果比較.(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中比較了不同加權機制下的性能.從表中可以發(fā)現(xiàn):(1)本文算法的分量數(shù)目多一些.這是由于自適應尺度像素分解算法采用了有效集,部分分量被多次優(yōu)化,有利于減少分量數(shù)目.(2)單次迭代所花時間縮短了3倍以上,同時單次迭代所花的時間也有較大的差異.這種差異部分來源于自適應尺度像素分解算法中的有效集大小的不同.當其他參數(shù)確定時,有效集的大小取決于點擴展函數(shù)的旁瓣水平.(3)整體上,本文算法速度提高了3倍左右.

表1 在不同加權機制下的性能比較Table 1 Performances with different weighting functions

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

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

4 總結與討論

自適應尺度像素分解算法在射電綜合成圖中可以獲得當前最好的重建效果,然而其計算時間卻很長.本文針對這個問題,在自適應尺度像素分解算法基礎上,使用不同的初始值計算方法和擬合策略.實驗表明,本算法在獲得高質量重建圖像的同時,反卷積速度有較為明顯的提升.

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