最大熵與MP-CLEAN方法對擴展源圖像重建的比較

張茂林，周建鋒

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

?

最大熵與MP-CLEAN方法對擴展源圖像重建的比較

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

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

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

關鍵詞：圖像重建；射電干涉陣；Multi-Point CLEAN；最大熵方法；擴展源

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

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

1模型的建立

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

(1)

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

(2)

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

(3)

(4)

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

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

(5)

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

2重建方法

2.1最大熵方法

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

(6)

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

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

(7)

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

(8)

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

(9)

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

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

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

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

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

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

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

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

2.2MP-CLEAN方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由于采樣函數(shù)沒有對可見度的原點(直流分量)進行采樣，即總流量需要估計，而這里模擬中流量估計認為準確，即在頻域采樣時對可見度的原點多采了一點。數(shù)值模擬中輸入信息是圖1右下圖對應的采樣得到的可見度數(shù)據(jù)(即模擬觀測數(shù)據(jù))、估計的噪聲、估計的流量及采樣函數(shù)。

圖1左上：模擬選用的擴展源，單位：Jy；左下：12 h觀測對應的采樣函數(shù)，又稱為UV覆蓋(uv-coverage)，兩坐標軸單位都是kλ(1kλ≈60 m)；右上：原圖可見度數(shù)據(jù)的幅度的對數(shù)，右下：模擬觀測(采樣)得到的可見度數(shù)據(jù)(已加噪聲)的幅度的對數(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ù)，可以通過最大熵或者MP-CLEAN等方法對真圖進行重建。然而由于采樣數(shù)據(jù)本身的缺失，對于模擬中用到的復雜擴展源而言，任何方法重建的結(jié)果都不可能恢復到與復雜的原始圖像完全一致，因為任何反卷積方法都可看成在缺失的頻域采樣處插值，重建結(jié)果只能是基于一定的模型對真圖的近似估計。

3.2模擬重建比較

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

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

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

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

圖2MP-CLEAN及最大熵方法對8 h及12 h模擬觀測的重建結(jié)果。左上：真圖；右上：8 h觀測對應的臟圖；中間和下面各4幅子圖，分別對應8 h和12 h觀測數(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é)果剖面殘差圖的均值及標準差

Table 1The mean and standard deviations of profile

differences from the reconstruction

圖3MP-CLEAN及最大熵方法對8 h(上圖)及12 h(下圖)觀測的重建結(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對12 h觀測數(shù)據(jù)進行重建的χ2及殘差(絕對值)最大值隨迭代步數(shù)及時間變化的關系。從圖4中χ2與迭代步數(shù)之間的關系可以發(fā)現(xiàn)，MP-CLEAN只需要不到70步，χ2即達到1，而最大熵則需要140多步。從χ2隨時間的變化曲線可以發(fā)現(xiàn)，MP-CLEAN只需約15 s，而最大熵則需要55 s左右。結(jié)合迭代步數(shù)及時間，可以算出MP-CLEAN平均每步運算所需時間為0.2 s左右，而最大熵大約要0.4 s。由此可見MP-CLEAN的每一次迭代的速度約為最大熵方法的兩倍，而MP-CLEAN需要的重建步數(shù)是最大熵方法的一半，最終導致最大熵方法所需的重建時間是MP-CLEAN的3倍多(接近4倍)。

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

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

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

4討論

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

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

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

處的剖面殘差曲線；中心線右邊為第2個源的重建結(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

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

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

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

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

隨迭代變化的關系曲線

Fig.6The natural logarithm ofχ2versus the

iteration whenαvaries

圖7不同α取值下最大熵對8 h觀測的

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

Fig.7The profile differences of the MEM reconstruction

results within the 8h whenαvaries

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

隨迭代變化的關系曲線

Fig.8The natural logarithm ofχ2versus the

iteration whenγvaries

圖9不同γ取值下MP-CLEAN對8 h觀測的

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

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

struction results within the 8h whenγvaries

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

剖面殘差的均值及標準差

Table 2The mean and standard deviation of profile

differences using MEM whenαvaries

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

的剖面殘差的均值及標準差

Table 3The mean and standard deviation of profile diff-

erences using the MP-CLEAN whenγvaries

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

5總結(jié)

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

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

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

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

(4)最大熵方法及MP-CLEAN方法的重建效果與參數(shù)的選擇有關，最大熵方法對參數(shù)的依賴性要小于MP-CLEAN方法，即最大熵方法的魯棒性優(yōu)于MP-CLEAN方法。

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

參考文獻：

[1]田海俊, 徐洋, 陳學雷, 等. 射電干涉陣 GPU 相關器的研究初探[J]. 天文研究與技術——國家天文臺臺刊, 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]徐永華, 汪敏, 郝龍飛, 等. 太陽低頻射電干涉陣的構建仿真[J]. 天文研究與技術——國家天文臺臺刊, 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.

基金項目：國家自然科學基金 (11173038, 11373025)；清華大學自主科研計劃 (20111081102) 資助.

收稿日期：2015-04-16；

修訂日期：2015-06-09

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

中圖分類號：P164

文獻標識碼：A

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