張 利，肖一凡，米立功，盧 梅，趙慶超，王 蓓，劉 祥，張 明，謝 泉

(1.貴州大學 大數(shù)據(jù)與信息工程學院，貴州 貴陽 550025；2.中國科學院新疆天文臺，新疆 烏魯木齊 830011)

通過觀測所接收到的電磁輻射來認識宇宙結(jié)構(gòu)，理解宇宙起源及演化，是人類不斷追求的目標。在射電波段對天體和宇宙開展研究的天文學分支稱為射電天文學。近年來，射電天文學取得了顯著突破，包括脈沖星和中子星[1]、星系中的暗物質(zhì)[2]、由超大質(zhì)量黑洞驅(qū)動的射電星系[3]和類星體[4]等，這對于理解宇宙、恢復宇宙圖景具有重要意義。

研究宇宙演化的重要手段之一是研究射電天文圖像及其結(jié)構(gòu)變化。由于銀河系外天體(例如類星體和星系)的射電輻射到達地球時信號已經(jīng)非常微弱，只有借助具有較高靈敏度的大型射電干涉陣列才能檢測到這些信號。目前已建成的大型干涉陣列包括巨米波射電望遠鏡(giant metrewave radio telescope， GMRT)、默奇森寬場陣列望遠鏡(murchison widefield array，MWA)、低頻陣列射電望遠鏡(low frequency array，LOFAR)，此外還有若干正在建設的新型干涉陣列，例如平方公里射電陣(square kilometre array，SKA)。

隨著越來越多的巡天計劃的實施和大型干涉陣列的建成，人們能夠得到海量天文數(shù)據(jù)，但是望遠鏡數(shù)量有限和地球自轉(zhuǎn)綜合造成的稀疏采樣使射電天文圖像變得模糊。為了對源的結(jié)構(gòu)特性做更精準的分析，需要對觀測到的射電天文圖像進行去模糊處理。H?gbom[5]開創(chuàng)性地提出了CLEAN反卷積算法，通過迭代識別點源并消除點擴展函數(shù)(PSF)的影響，有效解決由于稀疏采樣引起的圖像模糊問題。該算法的提出具有超越時代的意義，并延伸應用到其它領域。CLARK算法[6]和Cotton-Schwab算法[7]對其作出改進，在提高反卷積速度的同時，減少算法誤差，對展源的處理更加精準。然而點源分解機制和點擴展函數(shù)的旁瓣導致包含延展源和復雜特征的重建圖像中出現(xiàn)條紋(stripes)。研究人員提出了一些改進算法解決該問題，例如，尺度敏感的CLEAN算法[8-12]，利用尺度基函數(shù)將天文圖像參數(shù)化，以表達像素之間的相關(guān)性，從而進一步消除圖像殘差，但對于重建圖像質(zhì)量仍有提升空間。本文提出一種自適應尺度CLEAN算法，通過局部極小化目標函數(shù)，尋找最優(yōu)分量，以提高反卷積的性能，從而更加精準地重建射電天文圖像。

1 理論基礎

1.1 干涉成像原理

天空亮度分布可以分解為奇對稱成分和偶對稱成分，在射電天文望遠鏡進行動態(tài)范圍成像時，使用一對cosine 和 sine 相關(guān)器的組合，即復相關(guān)器，來輸出天空圖像的可見度：

(1)

(2)

在滿足小視場成像和所有基線矢量共面的情況下，上式可消去w額外項，變?yōu)槎S傅里葉變換，然后通過逆傅里葉變換，得到天空亮度分布。

圖1 干涉儀成像中常用的(u,v,w)直角坐標系[13]Fig.1 The (u,v,w) rectangular coordinate system commonly used in interferometer imaging

受限于干涉陣列的天線數(shù)目和基線的長度范圍，所觀測的uv平面覆蓋不完全，因此無法獲得天空亮度分布的全部信息，傅里葉逆變換僅能得到臟圖Idirty(l,m)，

(3)

其中，S(u,v)表示采樣函數(shù)。采樣函數(shù)是干涉陣列的點擴展函數(shù)B(l,m)的傅里葉變換。根據(jù)卷積定理，上式可表示為：

Idirty(l,m)=I(l,m)*B(l,m)。

(4)

臟圖存在明顯的旁瓣等干擾，無法反映真實的天空亮度分布，因此利用反卷積算法，例如CLEAN算法、最大熵算法[14]，從臟圖中盡可能地對真實射電天文圖像進行重建。

1.2 H?gbom CLEAN算法

H?gbom提出的CLEAN算法將天空亮度分布假設為一系列點源集合，在圖像域中通過迭代潔化的過程，消除臟圖中的旁瓣效應。對H?gbom CLEAN(Hg-CLEAN)算法的具體實現(xiàn)過程如下：

步驟1尋找臟圖中最大絕對亮源，取其位置(xi,yi)和幅度ai作為初始參數(shù)；

(5)

其中，g表示循環(huán)增益；

步驟3更新第i次殘差圖像，Iiresidual=Idirty-B(l,m)*Iimodel。

(6)

其中，*表示卷積運算；

步驟4重復步驟1—3，直到迭代次數(shù)超過預設閾值或殘差圖像中的最大值小于規(guī)定的噪聲水平時停止。將潔束卷積模型圖像，再加上殘差圖像，得到恢復出來的圖像。

這種算法對于點源集合的處理具有很好的效果，然而點源無法表示延展源中像素間的相關(guān)性，因此不能很好地重建展源。

2 基于L-BFGS-B局部極小化的自適應尺度CLEAN算法

受CLEAN算法[5,12]和文獻[15]的啟發(fā)，本文提出一種L-BFGS-B局部極小化的自適應尺度CLEAN算法，和其它CLEAN算法的變種一樣，使用相同的框架[10]。本文引入的L-BFGS-B算法由L-BFGS算法[16]改進得到，是一種基于梯度投影的非線性優(yōu)化方法。該方法結(jié)合原有算法的Hessian、線搜索算法和信任域方法應用于更新過程。

圖2 本文算法框架Fig.2 The framework for our algorithm

L-BFGS-B算法的迭代公式為：

xi+1=xi-αiHi▽fi。

(7)

其中，αi表示步長調(diào)節(jié)因子，Hi為Hessian矩陣，▽fi為多維一階向量梯度。Hi通過公式在每次迭代中更新參量，

(8)

Vi=I-ρiyisiT；

(9)

si=xi+1-xi；

(10)

yi=▽fk+1-▽fk。

(11)

在本文算法的每一次迭代中，L-BFGS-B算法按照以下步驟尋找最優(yōu)解：

步驟1結(jié)合信任域方法，根據(jù)下列公式在點xi附近用二次型模型對目標函數(shù)進行局部擬合，計算柯西點的近似值，

(12)

其中，()T表示轉(zhuǎn)置運算，Hi表示在第i次迭代的Hessian矩陣；

步驟2 通過直接法、共軛梯度法或?qū)ε挤ㄓ嬎闼阉鞣较騞k；

步驟3依據(jù)問題的約束條件，沿搜索方向dk進行線性搜尋，計算步長調(diào)節(jié)因子αk，更新參數(shù)xk+1=xk+αkdk，以尋找函數(shù)的最小值；

步驟4采用L-BFGS-B算法更新Hessian矩陣Hk，同時檢查是否收斂。

重復步驟1—4，直到滿足L-BFGS-B算法的3個停止準則之一：達到最大迭代次數(shù)、目標函數(shù)的減少量變小時，投影梯度的范數(shù)足夠小時，循環(huán)終止。從而得到函數(shù)最小值，閾值的選擇取決于臟圖和點擴展函數(shù)的實際情況。

然后，利用一系列高斯函數(shù)對殘差圖像進行平滑處理，取結(jié)果的最大值作為初始參數(shù)，

(13)

其中，ωi表示第i次迭代中重新計算得到的高斯模型的尺寸，ωic表示第i次迭代中由最小化算法擬合殘差得到的高斯模型尺寸，ωb是點擴展函數(shù)的高斯模型尺寸的近似值。優(yōu)化后的幅值ai可由下式計算得到：

(14)

其中，aic和ab分別表示擬合殘差得到的幅值和點擴展函數(shù)的模型幅值。接著，計算模型尺度的大小，并根據(jù)高斯尺度得到第i次擬合后的模型分量：

(15)

使用該最優(yōu)分量對模型圖像進行更新，同時利用循環(huán)增益g對模型分量進行優(yōu)化，以建立更加精確的天空模型，

(16)

其中，Iimodel表示第i次迭代的模型圖像。接著利用模型圖像計算殘差圖像，更新第i次殘差圖像，

Iiresidual=Idirty-B(l,m)*Iimodel。

(17)

其中，*表示卷積運算。對模型和殘差圖像進行循環(huán)更新，直到滿足最大迭代次數(shù)或者達到噪聲水平時停止，得到最后的模型圖像Imodel和殘差圖像Iresidual，重建后的圖像Irestored表示如下公式：

Irestored=Bclean*Imodel+Iresidual。

(18)

其中，Bclean為擬合主瓣得到的干凈光束。改進后的算法流程如圖2所示。

本文采用L-BFGS-B算法對分量進行更準確的擬合。對比Hg-CLEAN算法，該算法通過對最小化目標函數(shù)部分的優(yōu)化，能夠?qū)μ炜樟炼确植甲龀龈珳实慕＃怪亟ê蟮膱D像更加接近真實天空圖像。

3 試驗結(jié)果分析與算法評估

本文使用天文通用軟件包(common astronomy software applications,CASA)模擬JVLA B陣型觀測仙女星系M31圖像，并用于驗證新算法的性能，圖像尺寸為256×256像素。圖3展示了本算法重建圖像前后的對比效果。對比圖3(a)中參考的干凈天空圖像，圖3(b)顯示了在點擴展函數(shù)旁瓣影響下的臟圖，其中包含大量的旁瓣，造成圖像中天體的細節(jié)結(jié)構(gòu)模糊不清。圖3(c)是本文算法實現(xiàn)的圖像重建效果，從圖中能夠觀察到，臟束的旁瓣效應基本被消除，天空圖像中的點源和展源附近的模糊得到較好的處理，說明本文算法能夠從臟圖中有效重建天空圖像。

(a)原始圖像；(b)臟圖；(c)重建圖像圖3 本算法重建圖像對比效果Fig.3 Contrast effect of the algorithm

相對于Hg-CLEAN算法，本文算法呈現(xiàn)出更好的圖像重建性能。圖4是本文算法與傳統(tǒng)的Hg-CLEAN算法對于同一圖像的處理結(jié)果，展示了對應的模型圖像和殘差圖像。從圖中可以看出，本文算法得到的殘差圖像中的結(jié)構(gòu)遠少于Hg-CLEAN算法，說明采用該算法得到的重建圖像，能夠恢復原始圖像中的大部分信息，同時能夠體現(xiàn)L-BFGS-B算法的擬合參數(shù)效果良好，提高了算法對天空亮度分布建模的精度；觀測最終的重建結(jié)果，發(fā)現(xiàn)Hg-CLEAN算法也能夠?qū)υ紙D像的部分信息進行重建，然而在最終的結(jié)果中仍帶有大量的偽影，并且在其殘差圖像中包含明顯的圖像結(jié)構(gòu)。

(a1)Hg-CLEAN算法的模型圖像;(a2)Hg-CLEAN算法的殘差圖像;(b1)本文算法的模型圖像; (b2)本文算法的殘差圖像圖4 不同算法的重建結(jié)果Fig.4 Reconstructed results from different algorithms

另外，從圖4(a2)中能夠觀察到明顯的延展信號，圖5展示了模型圖像中展源部分的處理效果。與Hg-CLEAN算法相比，本文提出的算法在展源處理方面效果更好，原因在于Hg-CLEAN算法將天空圖像分解為一系列delta函數(shù)，不能很好地表達圖像中的展源結(jié)構(gòu)，而本文算法能夠依據(jù)射電天文圖像中的展源結(jié)構(gòu)自適應地選取尺度參數(shù)，因此具有處理圖像中延展源的性能。

(a)原始圖像;(b)Hg-CLEAN算法的處理結(jié)果;(c)本文算法的處理結(jié)果圖5 展源部分的局部特征圖像Fig.5 The local feature image of the source part

同時，采用客觀指標對天空圖像重建效果做量化分析，選取峰值信噪比(peak signal-to-noise ratio,PSNR)[17]、結(jié)構(gòu)相似性(structural similarity,SSIM)[18]和均方根誤差(root mean square error,RMSE)[19]作為重建圖像的質(zhì)量評價指標，對Hg-CLEAN算法和本文提出的算法進行評價。為了對比不同算法之間的差異，在同樣運行環(huán)境下對不同算法進行測試，從表1可以看出，本文提出算法在PSNR、SSIM和RMSE指標上均有所提升。

表1 重建圖像質(zhì)量指標比較Tab.1 Reconstructed image quality index comparison

為了驗證算法的穩(wěn)定性，采用三種加權(quán)方案對原始圖像做處理，即自然加權(quán), 均勻加權(quán)和Briggs 加權(quán)。圖6展示了不同加權(quán)方案下圖像的結(jié)構(gòu)變化，觀察可以發(fā)現(xiàn)經(jīng)自然加權(quán)后的圖像中的旁瓣效應最為明顯，均勻加權(quán)后的圖像的旁瓣效應最弱。

(a)自然加權(quán);(b)均勻加權(quán);(c)Briggs加權(quán)圖6 不同加權(quán)方案下的臟圖Fig.6 Dirty images of different weighting schemes

圖7顯示了對于不同加權(quán)方案的重建結(jié)果。對比殘差圖像，本算法對不同加權(quán)下的臟圖均有較好的重建效果，尤其對于均勻加權(quán)和Briggs加權(quán)具有明顯優(yōu)勢。

(a1)自然加權(quán)的重建圖像;(a2)自然加權(quán)的殘差圖像;(b1)均勻加權(quán)的重建圖像;(b2)均勻加權(quán)的殘差圖像; (c1)Briggs加權(quán)的重建圖像;(c2)Briggs加權(quán)的殘差圖像圖7 不同加權(quán)方案下的重建效果對比Fig.7 Comparison of reconstruction effects under different weighting schemes

4 結(jié)語

CLEAN算法對于射電天文圖像處理領域具有重大意義，并能夠延伸到天文學相關(guān)的其他領域[20]，這系列算法對傳統(tǒng)的圖像處理算法造成沖擊，對圖像處理的相關(guān)領域也產(chǎn)生了深遠影響。本文基于目前已有的CLEAN算法[12]，結(jié)合L-BFGS-B算法對目標函數(shù)的最小值求解進行優(yōu)化，并將算法的重建結(jié)果與傳統(tǒng)的Hg-CLEAN算法[5]做比較，提高了原有算法的性能，使重建得到的射電天文圖像更加接近真實天空圖像，為測試該算法對不同觀測圖像的普適性，在以后的工作中，計劃擴展到不同陣列和不同科學目標的觀測成像中去。