張茂林，周建鋒

(1. 清華大學(xué)工程物理系天體物理中心，北京　100084；2. 粒子技術(shù)與輻射成像教育部重點(diǎn)實(shí)驗(yàn)室(清華大學(xué))，北京　100084)

?

最大熵與MP-CLEAN方法對(duì)擴(kuò)展源圖像重建的比較

張茂林1,2，周建鋒1,2

(1. 清華大學(xué)工程物理系天體物理中心，北京100084；2. 粒子技術(shù)與輻射成像教育部重點(diǎn)實(shí)驗(yàn)室(清華大學(xué))，北京100084)

摘要：在射電綜合成像領(lǐng)域，通常需要用退卷積的方法補(bǔ)全頻域稀疏的采樣。由于擴(kuò)展源的頻域信息更為豐富，要補(bǔ)全這些信息相對(duì)于點(diǎn)源來(lái)說(shuō)更難，因此擴(kuò)展源的圖像重建是射電綜合成像領(lǐng)域的一大難點(diǎn)。為了探索射電干涉擴(kuò)展源圖像重建方法的特點(diǎn)，將最大熵與一種加速的CLEAN方法(文中稱(chēng)之為Multi-Point CLEAN，MP-CLEAN)對(duì)擴(kuò)展源的干涉陣模擬數(shù)據(jù)的圖像重建進(jìn)行了比較。通過(guò)比較，發(fā)現(xiàn)對(duì)于同樣的觀測(cè)數(shù)據(jù)，兩種方法都能較好地重建圖像，但MP-CLEAN方法的旁瓣祛除效果與重建效果優(yōu)于最大熵方法，而且在模擬數(shù)據(jù)重建中MP-CLEAN方法的總體速度比最大熵快3倍以上。最后，在討論部分通過(guò)研究?jī)煞N方法中參數(shù)的選擇對(duì)重建結(jié)果的影響，發(fā)現(xiàn)最大熵方法比MP-CLEAN方法對(duì)參數(shù)選擇的依賴(lài)性弱，這表明最大熵方法有更好的魯棒性。

關(guān)鍵詞：圖像重建；射電干涉陣；Multi-Point CLEAN；最大熵方法；擴(kuò)展源

射電天文是天文學(xué)的一個(gè)分支，射電天文顧名思義是通過(guò)探測(cè)射電波觀測(cè)天體，從而研究天體的一門(mén)科學(xué)。射電波的波長(zhǎng)相對(duì)于可見(jiàn)光長(zhǎng)得多，這使得射電觀測(cè)的角分辨率受到很大的限制，以6 cm的射電波為例，其波長(zhǎng)是可見(jiàn)光(波長(zhǎng)400至760 nm)的十萬(wàn)倍左右，也就是說(shuō)，一個(gè)1 m直徑的光學(xué)望遠(yuǎn)鏡的分辨率相當(dāng)于直徑約100 km的射電望遠(yuǎn)鏡，顯然分辨率是射電天文發(fā)展的一個(gè)很大的阻礙。從19世紀(jì)40年代起，干涉成像技術(shù)的發(fā)展克服了這種局限。干涉成像技術(shù)的創(chuàng)始人Martin Ryle也因此獲得了1974年的諾貝爾獎(jiǎng)。接下來(lái)的數(shù)十年，一些大型的射電干涉陣列(如One-Mile Telescope、 Green Bank Interferometer、 Westerbork Synthesis Radio Telescope、 Very Large Array)陸續(xù)建成，大大推動(dòng)了射電天文學(xué)的發(fā)展。如今的射電干涉望遠(yuǎn)鏡已經(jīng)跨入世界上分辨能力最高的望遠(yuǎn)鏡行列，例如阿塔卡瑪毫米/亞毫米波陣列望遠(yuǎn)鏡(Atacama Large Millimeter Array, ALMA)的空間分辨能力已經(jīng)是哈勃空間望遠(yuǎn)鏡(Hubble Space Telescope)的5倍左右。

干涉陣是利用干涉原理得到可見(jiàn)度數(shù)據(jù)(Visibility Data)，這里的可見(jiàn)度數(shù)據(jù)可等效理解為圖像的空間頻譜數(shù)據(jù)。由于干涉是來(lái)自?xún)蓚€(gè)信號(hào)的相干，所以得到的可見(jiàn)度數(shù)據(jù)與觀測(cè)天線的數(shù)量及相對(duì)位置分布有關(guān)。受材料、建造技術(shù)及資金等限制，現(xiàn)實(shí)中只能建造天線數(shù)目和尺寸都有限的干涉陣。但是，隨著干涉技術(shù)的不斷發(fā)展和完善，人們已經(jīng)開(kāi)始籌建越來(lái)越大的陣列，例如我國(guó)的 “天籟計(jì)劃” 等，由此也帶來(lái)了數(shù)據(jù)處理的挑戰(zhàn)[1]。干涉陣的排布對(duì)成像有較顯著的影響，這方面國(guó)內(nèi)已有所研究[2]。由于射電干涉陣的觀測(cè)數(shù)據(jù)是對(duì)圖像頻域空間的稀疏采樣，因此如果直接對(duì)獲取的可見(jiàn)度數(shù)據(jù)進(jìn)行傅里葉逆變換，得到的圖像往往很模糊，通常存在很強(qiáng)的旁瓣。為了重建清晰的圖像，通常需要進(jìn)行反卷積(Deconvolution)。在射電天文圖像重建領(lǐng)域，常用的幾種反卷積方法有CLEAN方法[3-4]及其衍生方法(如Multi-scale CLEAN[5]、Multi-resolution CLEAN[6]等)和最大熵方法[7]及其衍生方法(如Multi-scale MEM[8]等)。射電干涉成像中，擴(kuò)展源的圖像重建是其中的難點(diǎn)，普通的CLEAN方法雖然對(duì)簡(jiǎn)單的點(diǎn)源重建效果較好，但對(duì)于復(fù)雜的大面積彌散源來(lái)說(shuō)重建時(shí)計(jì)算量大，收斂速度慢，重建穩(wěn)定性差，得到的結(jié)果通常不理想[9]。具有代表性的適用于擴(kuò)展源干涉陣數(shù)據(jù)的重建方法主要有兩種：一種是在CLEAN的基礎(chǔ)上發(fā)展起來(lái)的多尺度方法，以Cornwell的Multi-Scale CLEAN[5]為代表，一種是以最大熵方法[7](MEM)為代表。本文使用一種簡(jiǎn)單的加速CLEAN的方法，與CLEAN方法每次選擇殘余圖中的最亮點(diǎn)不同，每次選擇大于最亮點(diǎn)γ(0<γ<1)倍的多個(gè)點(diǎn)，而不僅僅是一個(gè)[10]，這種方法本文稱(chēng)為Multi-Point CLEAN方法(簡(jiǎn)稱(chēng)MP-CLEAN方法)。這種方法既加速了CLEAN，又在一定程度上利用了相鄰像素值之間的物理聯(lián)系，適合進(jìn)行復(fù)雜擴(kuò)展源的重建。本文之所以選擇用最大熵方法與MP-CLEAN方法進(jìn)行對(duì)比，是由于MP-CLEAN方法的簡(jiǎn)單性，它并不屬于多尺度方法，同樣地，最大熵方法也不屬于多尺度的方法，二者同屬于擴(kuò)展源重建的初級(jí)方法。因此可以將這兩種比較初級(jí)的方法進(jìn)行比較，通過(guò)對(duì)比兩種擴(kuò)展源重建方法，可以了解兩種方法的特點(diǎn)，并為進(jìn)一步完善和改進(jìn)射電干涉擴(kuò)展源圖像重建方法提供啟發(fā)。

1模型的建立

為了簡(jiǎn)化模型，本文中認(rèn)為天圖的可見(jiàn)度數(shù)據(jù)與天圖的空間頻譜等效，即直接使用傅里葉頻譜數(shù)據(jù)作為可見(jiàn)度數(shù)據(jù)。任意給定一個(gè)亮度分布I(x,y)，其可見(jiàn)度可以用(1)式表示：

(1)

式中，V(u,v)為可見(jiàn)度數(shù)據(jù)；I(x,y)為亮度分布圖像，即可見(jiàn)度是圖像的傅里葉變換。射電望遠(yuǎn)鏡的觀測(cè)等效于在天圖的頻域進(jìn)行采樣。如果準(zhǔn)確獲得可見(jiàn)度的全部采樣，可以直接利用傅里葉逆變換求得I(x,y)，如(2)式：

(2)

然而實(shí)際數(shù)據(jù)往往是對(duì)可見(jiàn)度數(shù)據(jù)的部分采樣，采樣函數(shù)如(3)式：

(3)

(4)

在重建時(shí)，已知數(shù)據(jù)是采樣得到的可見(jiàn)度數(shù)據(jù)和采樣函數(shù)(在圖像域表達(dá)就是已知臟圖和點(diǎn)擴(kuò)展函數(shù))，而需要重建真實(shí)圖像I(x,y)。

如果考慮觀測(cè)中存在的噪聲，可見(jiàn)度觀測(cè)數(shù)據(jù)可以用下式表示：

(5)

式中，V(uk)表示可見(jiàn)度觀測(cè)數(shù)據(jù)；I(xi)表示真實(shí)亮度分布；uk是UV平面(頻域)上第k個(gè)采樣點(diǎn)的位置向量；xi是圖像平面第i個(gè)像素的位置向量；nk表示UV平面上第k個(gè)采樣點(diǎn)的噪聲。有了(5)式描述的關(guān)系，要求出真實(shí)亮度分布，需要用反卷積方法。

2重建方法

2.1最大熵方法

在引入最大熵方法之前先給出一種熵的定義，文[11]于1991年提出了一種熵的形式，如(6)式：

(6)

式中，Hg表示圖像的(相對(duì))熵；b表示待重建圖像；m表示圖像模型；i表示第i個(gè)像素。

利用這種熵的定義，結(jié)合觀測(cè)數(shù)據(jù)，就能求解真實(shí)亮度分布。在觀測(cè)數(shù)據(jù)的約束下求滿足熵最大的解，可以利用拉格朗日乘子法構(gòu)造目標(biāo)函數(shù)：

(7)

式中，χ2為觀測(cè)可見(jiàn)度數(shù)據(jù)與重建圖像的可見(jiàn)度數(shù)據(jù)各點(diǎn)偏差大小的歸一化平方和；E(χ2)為χ2的期望值，χ2具體表達(dá)式為

(8)

使目標(biāo)函數(shù)取得最大值的圖像b就是待求的解，函數(shù)的極值問(wèn)題有多種求解方法，這里不妨轉(zhuǎn)化為求目標(biāo)函數(shù)梯度的0點(diǎn)，即可利用牛頓法求解，迭代步長(zhǎng)的具體表達(dá)式由下式得到[10]：

(9)

式中，n從0開(kāi)始取值，b(0)即為已知給定的先驗(yàn)圖像m，本文所做模擬中的先驗(yàn)圖像m是一張所有點(diǎn)取值為1的圖像；J為目標(biāo)函數(shù)；α取常數(shù)；q取點(diǎn)擴(kuò)展函數(shù)的主峰峰值。

實(shí)現(xiàn)步驟[12]：

(1)給出初始重建圖像b0，圖像模型m，這里令b0及m所有像素的值等于1；

(2)對(duì)重建圖像的傅里葉域用采樣函數(shù)進(jìn)行采樣，得到重采樣數(shù)據(jù)Vk′；

(3)利用重采樣數(shù)據(jù)、觀測(cè)采樣數(shù)據(jù)及采樣函數(shù)的有關(guān)參數(shù)，代入(9)式，求出迭代步長(zhǎng)Δb；

(4)更新重建圖像：b=b+Δb；

(5)判斷迭代終止條件是否達(dá)到，即歸一化卡方值開(kāi)始小于等于1。如果達(dá)到終止條件，則輸出圖像b即為重建圖像；如果未達(dá)到，則返回第(2)步；

(6)將重建圖像卷積與最大基線長(zhǎng)對(duì)應(yīng)分辨率的高斯函數(shù)即得到最終的最大熵重建結(jié)果。

2.2MP-CLEAN方法

在射電天文圖像重建中，CLEAN方法是最常用的反卷積方法之一，文[10]于1980年提出了一種可以同時(shí)對(duì)多個(gè)點(diǎn)進(jìn)行CLEAN來(lái)加速CLEAN的方法，為了描述和比較方便，本文稱(chēng)之為Multi-Point CLEAN(MP-CLEAN)。在應(yīng)用MP-CLEAN時(shí)，需要一條經(jīng)驗(yàn)性假設(shè)，即有采樣較好的數(shù)據(jù)對(duì)應(yīng)點(diǎn)擴(kuò)散函數(shù)的主瓣比旁瓣高很多，即使有旁瓣的疊加，不同亮源同時(shí)作用在同一位置的概率不大，形成假源的峰值也不會(huì)太大。這就為同時(shí)選擇多個(gè)真實(shí)的源提供了可能，因此可以在每一步迭代中選擇大于最亮值γ(0<γ≤1)倍的多個(gè)點(diǎn)，而不僅僅是一個(gè)點(diǎn)[8]來(lái)加速CLEAN。每一次迭代γ允許的取值范圍可能不一樣，但是，文中為了簡(jiǎn)便，γ在同一次重建中都取一個(gè)較為保守的常數(shù)(注γ≡1時(shí)為CLEAN方法)，參數(shù)的選取在文中第4部分詳細(xì)討論。下面是MP-CLEAN方法的實(shí)現(xiàn)步驟[10]：

(1)給出一個(gè)初始模型M0，并令M0所有元素為0；

(2)利用采樣函數(shù)求出由模型得到的可見(jiàn)度數(shù)據(jù)；

(3)將模型可見(jiàn)度數(shù)據(jù)與觀測(cè)數(shù)據(jù)作差，再做傅里葉逆變換，取實(shí)部得到殘差圖R；

(4)找出殘差圖R中的最大值MAX；

(5)找出殘差圖中大于MAXγ(0<γ≤1)倍的點(diǎn)，作為MP圖(與模型尺寸一樣的圖，圖中的點(diǎn)除了找出的這些點(diǎn)與原圖位置及像素值一樣，其余像素點(diǎn)全為0)；

(6)更新模型：在原模型上加MP圖的β(0<β<1)倍；

(7)判斷是否達(dá)到迭代終止條件: 殘差圖對(duì)應(yīng)的可見(jiàn)度數(shù)據(jù)的均方根小于等于估計(jì)噪聲水平(即卡方小于等于1時(shí)終止迭代)。如果未達(dá)到終止條件，則返回第(2)步。

(8)將結(jié)果卷積與最大基線長(zhǎng)對(duì)應(yīng)分辨率的高斯函數(shù)即得到MP-CLEAN重建結(jié)果。

3數(shù)值模擬及分析

3.1模擬數(shù)據(jù)的生成

模擬數(shù)據(jù)的生成從干涉陣的排列及組合開(kāi)始，采樣函數(shù)由天線陣列及觀測(cè)等因素確定，本文模擬所用的采樣函數(shù)由排列在同一條直線上的天線組成，其位置參考韋斯特博克綜合孔徑射電望遠(yuǎn)鏡(Westerbork Synthesis Radio Telescope, WSRT)的天線位置設(shè)置。觀測(cè)時(shí)利用地球自身的旋轉(zhuǎn)，一條基線對(duì)應(yīng)的點(diǎn)在采樣空間掃過(guò)一條弧線。為了簡(jiǎn)便，用一條弧線表示這種觀測(cè)。觀測(cè)時(shí)間越長(zhǎng)，掃過(guò)弧線長(zhǎng)度越長(zhǎng)，考慮到地球自轉(zhuǎn)的周期性，12小時(shí)等效180°的掃描。由于實(shí)信號(hào)對(duì)應(yīng)的可見(jiàn)度數(shù)據(jù)具有對(duì)稱(chēng)性，因此掃描過(guò)程等效于繞著中心直徑上的點(diǎn)旋轉(zhuǎn)采樣，這里掃描的弧線是橢圓弧，橢圓的長(zhǎng)軸與短軸之比為1∶sinθ(這里的θ為被觀測(cè)源所在的赤緯，文中θ取60°為例)。

選一張擴(kuò)展源作為原圖(圖1左上圖：圖像尺寸為1 200 × 1 200，這里每個(gè)像素大小設(shè)為1 arcsec × 1 arcsec)做傅里葉變換，即得到可見(jiàn)度(圖1右上圖)，再利用采樣函數(shù)(圖1左下圖)對(duì)可見(jiàn)度采樣并加上噪聲就可以得到模擬的可見(jiàn)度觀測(cè)數(shù)據(jù)(圖1右下圖)。干涉陣的最大基線長(zhǎng)度決定了其分辨率(在波長(zhǎng)一定的情況下)，所選的觀測(cè)波長(zhǎng)λ約為6 cm，得到的分辨率約為4.4 arcsec × 5.1 arcsec。由于觀測(cè)數(shù)據(jù)分辨率的限制，為了統(tǒng)一分辨率，在后面的重建圖像與真圖(真圖與原圖不同，真圖是原圖卷積高斯函數(shù)得到的分辨率為4.4 arcsec × 5.1 arcsec)的比較中，統(tǒng)一卷積一個(gè)歸一化的二維高斯函數(shù)。

由于模擬可見(jiàn)度數(shù)據(jù)是復(fù)數(shù)，因此選擇的噪聲也必須取復(fù)數(shù)形式，其幅度取服從均方差為5(單位：Jansky, Jy)、均值為0的高斯分布的絕對(duì)值，其相位取[0, 2π]上的均勻分布。在后面的模擬中假設(shè)知道噪聲的大小(即預(yù)估的噪聲水平是準(zhǔn)確的)。

由于采樣函數(shù)沒(méi)有對(duì)可見(jiàn)度的原點(diǎn)(直流分量)進(jìn)行采樣，即總流量需要估計(jì)，而這里模擬中流量估計(jì)認(rèn)為準(zhǔn)確，即在頻域采樣時(shí)對(duì)可見(jiàn)度的原點(diǎn)多采了一點(diǎn)。數(shù)值模擬中輸入信息是圖1右下圖對(duì)應(yīng)的采樣得到的可見(jiàn)度數(shù)據(jù)(即模擬觀測(cè)數(shù)據(jù))、估計(jì)的噪聲、估計(jì)的流量及采樣函數(shù)。

圖1左上：模擬選用的擴(kuò)展源，單位：Jy；左下：12 h觀測(cè)對(duì)應(yīng)的采樣函數(shù)，又稱(chēng)為UV覆蓋(uv-coverage)，兩坐標(biāo)軸單位都是kλ(1kλ≈60 m)；右上：原圖可見(jiàn)度數(shù)據(jù)的幅度的對(duì)數(shù)，右下：模擬觀測(cè)(采樣)得到的可見(jiàn)度數(shù)據(jù)(已加噪聲)的幅度的對(duì)數(shù)

Fig.1Top-left panel: extended source chosen for this paper (unit: Jy); Bottom-left panel: the UV coverage of a 12h observation; Top-right panel: the logarithm of the visibility amplitude to the original image; Bottom-left panel: the logarithm of the sampled visibility

利用模擬數(shù)據(jù)，可以通過(guò)最大熵或者M(jìn)P-CLEAN等方法對(duì)真圖進(jìn)行重建。然而由于采樣數(shù)據(jù)本身的缺失，對(duì)于模擬中用到的復(fù)雜擴(kuò)展源而言，任何方法重建的結(jié)果都不可能恢復(fù)到與復(fù)雜的原始圖像完全一致，因?yàn)槿魏畏淳矸e方法都可看成在缺失的頻域采樣處插值，重建結(jié)果只能是基于一定的模型對(duì)真圖的近似估計(jì)。

3.2模擬重建比較

本文所用模擬中α、β、γ3個(gè)參數(shù)的選取：α=5、β=0.1、γ=0.6(參數(shù)的選取可參看第4部分的討論，選擇這些參數(shù)的基本原則是使得收斂的卡方盡可能地接近1且保證重建的穩(wěn)定性)。

固定上面的幾個(gè)參數(shù)，就可以分別用最大熵及MP-CLEAN對(duì)模擬的可見(jiàn)度數(shù)據(jù)進(jìn)行重建，得到的重建結(jié)果如圖2。圖2展示了最大熵方法及MP-CLEAN方法利用8 h和12 h的模擬觀測(cè)數(shù)據(jù)重建的結(jié)果?？梢园l(fā)現(xiàn)直接利用采樣數(shù)據(jù)作傅里葉逆變換得到的臟圖(Dirty map)與真圖(Ture image)差別很大，主要體現(xiàn)在旁瓣很多，而且有條紋狀結(jié)構(gòu)。MP-CLEAN方法和最大熵方法的重建結(jié)果在不同程度上抑制了旁瓣，正中間的面源及它的延展結(jié)構(gòu)都恢復(fù)得較好，尤其是MP-CLEAN方法，重建的結(jié)果很清楚，同時(shí)淹沒(méi)在旁瓣中的點(diǎn)源也得到了較好的恢復(fù)，但是，最大熵方法依舊有一部分源沒(méi)能恢復(fù)，而且存在微弱可見(jiàn)的旁瓣。從殘余圖來(lái)看(殘余圖是重建圖像頻譜重采樣與觀測(cè)得到的可見(jiàn)度數(shù)據(jù)做差之后傅里葉逆變換得到)，最大熵的殘余圖中還能依稀看見(jiàn)真圖中間面源的影子，也可以看到許多黑色點(diǎn)源的殘跡，而MP-CLEAN方法得到的殘余圖中幾乎看不到這些蹤跡，這說(shuō)明MP-CLEAN重建更勝一籌。對(duì)比真圖、臟圖以及重建結(jié)果，總體來(lái)說(shuō)兩種方法的重建效果較好，重建結(jié)構(gòu)比較清楚，比臟圖的效果要好很多。

為了更直觀地分析它們之間的差別，在RA=601 arcsec處取一剖面(下文中所有的剖面都是指此位置處的剖面)如圖3。表1是從剖面的殘差圖(此處殘差圖是指真圖減去重建圖像得到的圖，注意與前面的殘余圖有區(qū)別)圖3中得到的幾個(gè)統(tǒng)計(jì)量(由于隨機(jī)噪聲的存在，對(duì)于不同次的重建，表1至表3中的所有統(tǒng)計(jì)量具有一定的變化幅度，但均值變化幅度小于0.5 Jy，標(biāo)準(zhǔn)差變化幅度一般不超過(guò)10%，由于這種變化很小，并不影響兩種方法結(jié)果的對(duì)比分析，因此本文只取某一次重建結(jié)果進(jìn)行分析)，均值越接近0且標(biāo)準(zhǔn)差越小，表明越接近真圖。最大熵與MP-CLEAN方法重建的圖像均值與臟圖相比小得多，說(shuō)明兩種方法的重建都保證了總流量恢復(fù)得很好。

從標(biāo)準(zhǔn)差來(lái)看，MP-CLEAN的結(jié)果無(wú)論是8 h還是12 h觀測(cè)中，其標(biāo)準(zhǔn)差都要比最大熵小。從剖面圖可以看出，即使在源較弱的圖像邊緣，MP-CLEAN的結(jié)果與真實(shí)圖像吻合得很好，而最大熵結(jié)果還存在微弱的旁瓣難以去除，說(shuō)明MP-CLEAN的旁瓣抑制效果比最大熵要好。

圖2MP-CLEAN及最大熵方法對(duì)8 h及12 h模擬觀測(cè)的重建結(jié)果。左上：真圖；右上：8 h觀測(cè)對(duì)應(yīng)的臟圖；中間和下面各4幅子圖，分別對(duì)應(yīng)8 h和12 h觀測(cè)數(shù)據(jù)的重建結(jié)果。重建結(jié)果中的4幅子圖分別為：右上，最大熵重建結(jié)果；左上，最大熵結(jié)果的殘余圖；左下，MP-CLEAN重建結(jié)果；右下：MP-CLEAN重建的殘余圖

Fig.2The reconstruction results of the 8h and 12h observations from the MP-CLEAN and the MEM. Top-left panel: the True image; Top-right panel: the dirty image made from the 8h observation. The two four-subplots at the middle and bottom of the image correspond to the 8h and 12h observation respectively. For each four subplots, Middle-left panel: the MEM result; Middle-right panel: the residual of the MEM result; Bottom-left panel: the MP-CLEAN result; Bottom-right panel: the residual of the MP-CLEAN result

表1重建結(jié)果剖面殘差圖的均值及標(biāo)準(zhǔn)差

Table 1The mean and standard deviations of profile

differences from the reconstruction

圖3MP-CLEAN及最大熵方法對(duì)8 h(上圖)及12 h(下圖)觀測(cè)的重建結(jié)果的剖面殘差圖

(在剖面RA=601 arcsec處真圖減去臟圖、最大熵及MP-CLEAN的結(jié)果)

Fig.3The profile difference map of reconstruction results of the 8h (Top subplot) and 12h

(bottom subplot) observation (true image minus other images atRA=601 arcsec)

除了重建的效果，再看兩種方法的重建速度，圖4給出了最大熵和MP-CLEAN對(duì)12 h觀測(cè)數(shù)據(jù)進(jìn)行重建的χ2及殘差(絕對(duì)值)最大值隨迭代步數(shù)及時(shí)間變化的關(guān)系。從圖4中χ2與迭代步數(shù)之間的關(guān)系可以發(fā)現(xiàn)，MP-CLEAN只需要不到70步，χ2即達(dá)到1，而最大熵則需要140多步。從χ2隨時(shí)間的變化曲線可以發(fā)現(xiàn)，MP-CLEAN只需約15 s，而最大熵則需要55 s左右。結(jié)合迭代步數(shù)及時(shí)間，可以算出MP-CLEAN平均每步運(yùn)算所需時(shí)間為0.2 s左右，而最大熵大約要0.4 s。由此可見(jiàn)MP-CLEAN的每一次迭代的速度約為最大熵方法的兩倍，而MP-CLEAN需要的重建步數(shù)是最大熵方法的一半，最終導(dǎo)致最大熵方法所需的重建時(shí)間是MP-CLEAN的3倍多(接近4倍)。

圖4χ2及殘差絕對(duì)值最大值隨迭代步數(shù)及時(shí)間變化的關(guān)系曲線

Fig.4Theχ2and the maximum absolute residual versus the iteration or time

為了進(jìn)一步檢驗(yàn)這兩種方法的重建特點(diǎn)，選擇更多的原圖進(jìn)行模擬觀測(cè)，然后重建，采樣時(shí)間為12 h，采樣陣列保持不變。為了體現(xiàn)與上面的擴(kuò)展源不一樣，選用的源一個(gè)是規(guī)則的源(月球)，另一個(gè)是兩個(gè)不規(guī)則星系的原圖。重建得到的結(jié)果如圖5，從圖5的重建結(jié)果不難發(fā)現(xiàn)，與前面分析的重建特點(diǎn)一致：流量恢復(fù)均較好，都能較好地重建，而且MP-CLEAN對(duì)旁瓣的祛除更好，其剖面殘差的方差更小。限于篇幅原因，關(guān)于這些源重建速度的曲線并未給出，這兩個(gè)源的實(shí)驗(yàn)結(jié)果是：第1個(gè)源(Source 1)的重建中，MP-CLEAN的迭代時(shí)間不到15 s，而最大熵則花了近50 s(是MP-CLEAN的3倍多)；第2個(gè)源(Source 2)的重建中，MP-CLEAN的迭代時(shí)間不到30 s，而最大熵的迭代時(shí)間卻超過(guò)了150 s(是MP-CLEAN的5倍多)。模擬結(jié)果表明MP-CLEAN比最大熵快3倍以上，與前面的結(jié)果保持一致。

4討論

在3.2節(jié)已經(jīng)給出了最大熵及MP-CLEAN的重建結(jié)果，并且進(jìn)行了詳細(xì)的對(duì)比分析。然而，最大熵與MP-CLEAN方法的重建效果與重建參數(shù)的選擇有關(guān)，下面是利用前面8 h的觀測(cè)數(shù)據(jù)在不同參數(shù)選擇情況下的重建，通過(guò)對(duì)比和分析重建結(jié)果，進(jìn)一步了解兩種重建方法的特點(diǎn)。

首先改變最大熵方法中的參數(shù)α，得到圖6中不同α情況下χ2的變化曲線，可以看到如果不斷地增大α，χ2會(huì)迅速地降到1，重建結(jié)果的標(biāo)準(zhǔn)差較大，而如果不斷減小α到0附近，那么重建結(jié)果與真圖就有較大的整體偏差。這可以從表2中看出明顯的規(guī)律：隨著α的增大，均值的變化趨向0，然后穩(wěn)定在0附近；從標(biāo)準(zhǔn)差來(lái)看，隨著α的增大，標(biāo)準(zhǔn)差逐漸增大。

圖5最大熵及MP-CLEAN方法對(duì)于不同擴(kuò)展源的重建結(jié)果。中心線左邊為第1個(gè)源的重建結(jié)果及在RA=377 arcsec

處的剖面殘差曲線；中心線右邊為第2個(gè)源的重建結(jié)果及在Dec=601 arcsec處的剖面殘差曲線

Fig.5The reconstruction results from the MEM and MP-CLEAN method for different sources. The left panel is the reconstruction

results of Source 1 and the plot of the profile difference map atRA=377 arcsec; the right panel is the reconstruction

results of Source 2 and the plot of the profile difference map atDec=601 arcsec

這與α作用在觀測(cè)數(shù)據(jù)上有關(guān)，隨著α增大，觀測(cè)數(shù)據(jù)這個(gè)約束的權(quán)重越大，所以均值會(huì)收斂到0，而熵的權(quán)重就相對(duì)變小，熵本身是一種平滑性約束，熵權(quán)重減小，意味著重建結(jié)果的振蕩劇烈，所以標(biāo)準(zhǔn)差逐漸變大。在α=0.2時(shí)，不難發(fā)現(xiàn)旁瓣被抑制了，但是無(wú)論從圖6中χ2的收斂曲線，還是從表2中的均值來(lái)看，重建結(jié)果與真圖有整體上的偏離。從圖6中χ2的變化曲線可以看到，α=0.2時(shí)χ2收斂的值遠(yuǎn)大于1，從表2中的標(biāo)準(zhǔn)差來(lái)看，α=0.2時(shí)的結(jié)果振蕩最小，恰恰是熵作為一種平滑性約束的體現(xiàn)。重建結(jié)果的剖面殘差圖如圖7，在圖7中，除了α=0.2時(shí)重建結(jié)果與真圖有整體的偏移(從表2中對(duì)應(yīng)的均值可以看出)，難以看出重建結(jié)果之間較明顯的差別，這也說(shuō)明最大熵重建結(jié)果對(duì)參數(shù)的選擇有較好的魯棒性。

下面再看MP-CLEAN方法在γ(0<γ<1)取不同值時(shí)的重建情況，從圖8中的χ2隨迭代的變化曲線可以看到，γ越大則MP-CLEAN方法對(duì)擴(kuò)展源的重建速度越慢，圖像重建結(jié)果中的振蕩也變明顯；γ越小則重建結(jié)果越平滑，但是重建圖像中存在假源的可能性越大，γ越小假源的數(shù)目和強(qiáng)度的期望值都會(huì)變多和變大，γ太小(如圖8中γ=0.1)χ2難以降下來(lái)，只能收斂到大于1的值。

重建結(jié)果的剖面殘差圖如圖9，可以發(fā)現(xiàn)γ取不同的值對(duì)重建結(jié)果有較顯著的影響，對(duì)比圖9與圖7不難發(fā)現(xiàn)，最大熵方法在α取不同值時(shí)，重建結(jié)果之間的差別并不太大，而從表2與表3中的標(biāo)準(zhǔn)差對(duì)比中可以定量地說(shuō)明這一點(diǎn)，因此MP-CLEAN方法對(duì)影響其重建結(jié)果的參數(shù)γ的敏感程度較高，而最大熵則相對(duì)魯棒性更好。從表3中可以發(fā)現(xiàn)MP-CLEAN在不同參數(shù)取值下均值都接近0，但是方差則呈現(xiàn)中間(γ取0.5～0.7)低，兩側(cè)高的總體趨勢(shì)，這說(shuō)明MP-CLEAN的參數(shù)選擇受兩方面的約束，一個(gè)約束是重建結(jié)果中的假源應(yīng)當(dāng)數(shù)目少而且強(qiáng)度弱，另一個(gè)約束是結(jié)果的穩(wěn)定性。前一個(gè)約束要求γ盡可能地接近1，因?yàn)檫@樣重建結(jié)果含有假源的可能性才會(huì)變小，后一個(gè)約束則要求γ盡可能接近0，在這兩個(gè)約束的相互作用下，最后就呈現(xiàn)了表3中的方差中間低、兩側(cè)高的現(xiàn)象。采樣越多越全，γ可以選取的范圍越大，越有利于MP-CLEAN的重建。如果采樣太少，γ選擇的區(qū)間就會(huì)越小，甚至可能找不到合適的γ，此時(shí)再使用MP-CLEAN方法進(jìn)行復(fù)雜擴(kuò)展源的重建可能難以得到好的結(jié)果。

圖6不同α取值下χ2的自然對(duì)數(shù)

隨迭代變化的關(guān)系曲線

Fig.6The natural logarithm ofχ2versus the

iteration whenαvaries

圖7不同α取值下最大熵對(duì)8 h觀測(cè)的

重建結(jié)果剖面的殘差圖

Fig.7The profile differences of the MEM reconstruction

results within the 8h whenαvaries

圖8不同γ取值下χ2的自然對(duì)數(shù)

隨迭代變化的關(guān)系曲線

Fig.8The natural logarithm ofχ2versus the

iteration whenγvaries

圖9不同γ取值下MP-CLEAN對(duì)8 h觀測(cè)的

重建結(jié)果剖面的殘差圖

Fig.9The profile differences of the MP-CLEAN recon-

struction results within the 8h whenγvaries

表2最大熵方法在不同α取值下重建圖像的

剖面殘差的均值及標(biāo)準(zhǔn)差

Table 2The mean and standard deviation of profile

differences using MEM whenαvaries

表3MP-CLEAN在不同γ取值下重建圖像

的剖面殘差的均值及標(biāo)準(zhǔn)差

Table 3The mean and standard deviation of profile diff-

erences using the MP-CLEAN whenγvaries

由于β的取值對(duì)MP-CLEAN的影響主要體現(xiàn)在速度上，對(duì)重建圖像效果的影響不大(當(dāng)然β很大時(shí)也可能出現(xiàn)一定的不穩(wěn)定性)，因此β對(duì)MP-CLEAN重建的影響在此不再贅述。

5總結(jié)

本文利用最大熵方法與MP-CLEAN方法對(duì)模擬擴(kuò)展源的可見(jiàn)度數(shù)據(jù)進(jìn)行重建，可以發(fā)現(xiàn)：

(1)最大熵方法有著比MP-CLEAN更嚴(yán)格的數(shù)學(xué)迭代形式，MP-CLEAN方法則帶有一定的經(jīng)驗(yàn)性；

(2)對(duì)于同樣的采樣數(shù)據(jù)，MP-CLEAN對(duì)擴(kuò)展源的重建比最大熵方法重建效果好(使得χ2收斂到1附近時(shí)各參數(shù)的取值)，主要體現(xiàn)在旁瓣降得更低，重建圖像與真圖之間差別更小；

(3)MP-CLEAN方法的重建速度比最大熵快，從模擬結(jié)果來(lái)看，總體上所花時(shí)間不到最大熵的三分之一；

(4)最大熵方法及MP-CLEAN方法的重建效果與參數(shù)的選擇有關(guān)，最大熵方法對(duì)參數(shù)的依賴(lài)性要小于MP-CLEAN方法，即最大熵方法的魯棒性?xún)?yōu)于MP-CLEAN方法。

總而言之，就本文所用的模擬數(shù)據(jù)及重建結(jié)果而言，綜合模擬重建圖像的效果及重建的速度，MP-CLEAN比最大熵方法有更佳的重建效果及更快的重建速度。然而，由于MP-CLEAN對(duì)可見(jiàn)度數(shù)據(jù)采樣的要求較高，而且需要滿足源的旁瓣作用在同一位置的概率較小這種經(jīng)驗(yàn)性的假設(shè)，在采樣很差(點(diǎn)擴(kuò)散函數(shù)旁瓣峰值接近主峰峰值時(shí))或者在經(jīng)驗(yàn)性假設(shè)不成立時(shí)(例如亮源的某些對(duì)稱(chēng)性分布導(dǎo)致旁瓣疊加很?chē)?yán)重)，MP-CLEAN方法從原理上已經(jīng)不能適用，而最大熵方法理論上依舊能夠用于重建，這也是最大熵方法的魯棒性比MP-CLEAN好的體現(xiàn)之一。但MP-CLEAN方法在一定程度上挖掘了CLEAN方法在同一尺度上的潛能，可以為完善和改進(jìn)可見(jiàn)度采樣較好的擴(kuò)展源圖像重建方法提供更好、更快捷的途徑。

參考文獻(xiàn)：

[1]田?？? 徐洋, 陳學(xué)雷, 等. 射電干涉陣 GPU 相關(guān)器的研究初探[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2014, 11(3): 209-217.

Tian Haijun, Xu Yang, Chen Xuelei, et al. A preliminary study on GPU-based correlators for a radio interferometer array[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2014, 11(3): 209-217.

[2]徐永華, 汪敏, 郝龍飛, 等. 太陽(yáng)低頻射電干涉陣的構(gòu)建仿真[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2013, 10(3): 242-248.

Xu Yonghua, Wang Min, Hao Longfei, et al. Simulations of a low-frequency solar radio interferometry array[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2013, 10(3): 242-248.

[3]H?gbom J A. Aperture synthesis with a non-regular distribution of interferometer baselines[J]. Astronomy and Astrophysics Supplement, 1974, 15: 417-426.

[4]Schwarz U J. Mathematical-statistical description of the iterative beam removing technique (Method CLEAN) [J]. Astronomy and Astrophysics, 1978, 65: 345-356.

[5]Cornwell T J. Multiscale CLEAN deconvolution of radio synthesis images[J]. IEEE Journal of Selected Topics in Signal Processing, 2008, 2(5): 793-801.

[6]Wakker B P, Schwarz U J. The multi-resolution CLEAN and its application to the short-spacing problem in interferometry[J]. Astronomy and Astrophysics, 1988, 200: 312-322.

[7]Gull S F, Skilling J. Maximum entropy method in image processing[J]. IEE Proceedings F, 1984, 131(6): 646-659.

[8]Pantin E, Starck J L. Deconvolution of astronomical images using the multiscale maximum entropy method[J]. Astronomy and Astrophysics Supplement, 1996, 118(3): 575-586.

[9]Cornwell T J. A method of stabilizing the clean algorithm[J]. Astronomy and Astrophysics, 1983, 121(2): 281-285.

[10]Clark B G. An efficient implementation of the algorithm ‘CLEAN’[J]. Astronomy and Astrophysics, 1980, 89(3): 377-378.

[11]Skilling J, Gull S F. Bayesian maximum entropy image reconstruction[J]. Institute of Mathematical Statistics, 1991, 20: 341-367.

[12]Cornwell T J, Evans K F. A simple maximum entropy deconvolution algorithm[J]. Astronomy and Astrophysics, 1985, 143(1): 77-83.

CN 53-1189/PISSN 1672-7673

A Comparison of the MEM and the MP-CLEAN Methods for the Image Reconstruction of Extended Sources with Simulated Visibility Data

Zhang Maolin1,2, Zhou Jianfeng1,2

(1. Department of Engineering Physics and Center for Astrophysics, Tsinghua University, Beijing 100084, China,Email: zhangml13@mails.tsinghua.edu.cn; 2. Key Laboratory of Particle & Radiation Imaging(Tsinghua University) Ministry of Education, Beijing 100084, China)

Key words:Image reconstruction; Radio interferometer array; Multi-Point CLEAN; Maximum Entropy Method; Extended source

Abstract:In the radio synthesis imaging area, the deconvolution is used to fill up the un-sampled points at the Fourier domain. The algorithms like CLEAN and MEM often play an important role in the deconvolution. As extended sources have more complex visibility, the reconstruction for extended sourcesis difficult. In this paper we make a comparison of the MEM (Maximum Entropy Method) with the speed up CLEAN method (named Multi-Point CLEAN in this paper, or MP-CLEAN for short) for extended source reconstruction based on the simulated visibility data. The simulated data are produced by sampling the Flourier Frequency domain for a true image with the sample function generated from a simulated interferometer array. By analyzing and comparing, it is easy to find that both of these two methods work well on the image reconstruction while the MP-CLEAN is better than the MEM for getting rid of the side lobes. The MP-CLEAN is at least 3 times faster than the MEM in our simulations.A deeper analyzing of the selection of parameters shows that the MEM is more independent on the parameters than the MP-CLEAN, which indicates that the MEM is more robust.

基金項(xiàng)目：國(guó)家自然科學(xué)基金 (11173038, 11373025)；清華大學(xué)自主科研計(jì)劃 (20111081102) 資助.

收稿日期：2015-04-16；

修訂日期：2015-06-09

作者簡(jiǎn)介：張茂林，男，碩士. 研究方向：天體物理. Email: zhangml13@mails.tsinghua.edu.cn

中圖分類(lèi)號(hào)：P164

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào):1672-7673(2016)01-0100-11