基于L-BFGS-B局部極小化的自適應(yīng)尺度CLEAN算法

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

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

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

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

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

1 理論基礎(chǔ)

1.1 干涉成像原理

天空亮度分布可以分解為奇對(duì)稱(chēng)成分和偶對(duì)稱(chēng)成分，在射電天文望遠(yuǎn)鏡進(jìn)行動(dòng)態(tài)范圍成像時(shí)，使用一對(duì)cosine 和 sine 相關(guān)器的組合，即復(fù)相關(guān)器，來(lái)輸出天空?qǐng)D像的可見(jiàn)度：

(1)

(2)

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

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

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

(3)

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

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

(4)

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

1.2 H?gbom CLEAN算法

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

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

(5)

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

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

(6)

其中，*表示卷積運(yùn)算；

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

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

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

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

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

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

xi+1=xi-αiHi▽fi。

(7)

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

(8)

Vi=I-ρiyisiT；

(9)

si=xi+1-xi；

(10)

yi=▽fk+1-▽fk。

(11)

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

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

(12)

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

步驟2 通過(guò)直接法、共軛梯度法或?qū)ε挤ㄓ?jì)算搜索方向dk；

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

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

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

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

(13)

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

(14)

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

(15)

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

(16)

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

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

(17)

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

Irestored=Bclean*Imodel+Iresidual。

(18)

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

本文采用L-BFGS-B算法對(duì)分量進(jìn)行更準(zhǔn)確的擬合。對(duì)比Hg-CLEAN算法，該算法通過(guò)對(duì)最小化目標(biāo)函數(shù)部分的優(yōu)化，能夠?qū)μ炜樟炼确植甲龀龈珳?zhǔn)的建模，使重建后的圖像更加接近真實(shí)天空?qǐng)D像。

3 試驗(yàn)結(jié)果分析與算法評(píng)估

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

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