改進(jìn)的射電干涉成像凸優(yōu)化算法

溫 珊， 王俊義， 勞保強(qiáng)

(1.桂林電子科技大學(xué) 信息與通信學(xué)院,廣西 桂林 541004；2.中國(guó)科學(xué)院 上海天文臺(tái)，上海 200030)

射電干涉測(cè)量因?yàn)楦呓嵌确直媛屎挽`敏度在射電宇宙成像中起著關(guān)鍵作用[1]。然而，由于實(shí)際的測(cè)量條件的限制，在觀測(cè)期間獲得的能見度測(cè)量的數(shù)量受到陣列中天線對(duì)數(shù)量的限制。為了從不完整的數(shù)據(jù)中重建真實(shí)的天空亮度分布，解決這類不適定線性反演問題的圖像重建方法一直在發(fā)展。為了獲得更好的圖像質(zhì)量，作為魯棒框架的壓縮感知(CS)理論被提出來，處理從射電干涉測(cè)量中不完全線性測(cè)量的稀疏信號(hào)恢復(fù)。壓縮感知引入了一種優(yōu)于傳統(tǒng)奈奎斯特采樣范式的信號(hào)采集和重構(gòu)框架，該框架提出了采集信號(hào)技術(shù)的優(yōu)化和信號(hào)重建的非線性迭代算法[2]。Wiaux等[3]首次提出壓縮感知框架應(yīng)用于射電干涉測(cè)量中的圖像解卷積，在該框架中，壓縮感知方法在圖像保真度、靈活性和計(jì)算速度上優(yōu)于傳統(tǒng)成像方法。最近，Onose等[4]開發(fā)了求解RI成像問題的新算法—交替方向乘子法(ADMM)和原始對(duì)偶算法(PD)。它們依賴于近端分裂和前后迭代，在多數(shù)據(jù)，先驗(yàn)和圖像空間中并行運(yùn)行復(fù)雜的類似于CLEAN的迭代。鑒于此，基于交替方向乘子的凸優(yōu)化算法的數(shù)學(xué)框架，提出了一種改進(jìn)的凸優(yōu)化算法，它保留了與ADMM相同的計(jì)算負(fù)擔(dān)，但重構(gòu)質(zhì)量在理論上和實(shí)驗(yàn)中證明優(yōu)于前者。仿真和實(shí)驗(yàn)結(jié)果表明，該方法提高了重建圖像的精度和質(zhì)量。

1 基于壓縮感知的射電干涉成像

Donoho等[2]提出了針對(duì)稀疏特征信號(hào)的壓縮感知理論。理論表明，它以低于奈奎斯特采樣定理的速率對(duì)信號(hào)進(jìn)行采樣而不丟失信息，并且可以完整恢復(fù)信號(hào)。大多數(shù)自然信號(hào)在某些基也是稀疏或可壓縮的，例如天文圖像。壓縮感知的關(guān)鍵在于由向量x表示的復(fù)雜估值信號(hào)被認(rèn)為具有稀疏表示，x=Ψα，其中α∈￡D包含只有幾個(gè)非零元素，在字典Ψ∈￡M×N，例如小波基的集合，或更一般地說，一個(gè)過度完整的框架[5]。

干涉成像被定義為求解測(cè)量模型y=Φx+n的不適定逆問題，其中y∈￡M表示測(cè)量矢量，Φ∈￡M×N，其中M

(1)

其中:ε為噪聲的l2范數(shù)的上界；Ψ?為Ψ的伴隨算子。由于x是一個(gè)強(qiáng)度像，所以也假設(shè)了正先驗(yàn)性。數(shù)據(jù)保真度通過殘差低于噪聲ε限定。在這里，用最接近的凸松弛l1范數(shù)替代稀疏性的模型非凸的l0范數(shù)。由于l1最小化問題是一個(gè)凸問題，因此可以用有效的凸優(yōu)化算法來解決。

2 交替方向乘子法

為了解決式(1)中的不適定逆問題，引入Boyd等[8]提出的用于求解約束優(yōu)化問題的方法——交替方向乘子法(ADMM)。ADMM是旨在將對(duì)偶上升的可分解性與乘子方法的優(yōu)越收斂性相結(jié)合的一種算法。該算法通過式(2)解決問題：

minf(x)+g(x) subject toAx+Bz=c,

(2)

其中：變量x∈Rn且z∈Rm，A∈Rp×n，B∈Rp×m，c∈Rp。x和z是優(yōu)化變量。f(x)+g(x)是被由包含變量x的f(x)和涉及變量z的g(x)組成最小化的目標(biāo)函數(shù)。Ax+Bz=c是一個(gè)最小化約束。增廣的拉格朗日表達(dá)式為

Lp(x,z,y)=f(x)+g(z)+yT(Ax+Bz-c)+

(3)

(4)

uk+1=uk+Axk+1+Bzk+1-c。

(5)

為解決問題(4)，簡(jiǎn)單而有效的方法是引入Gauss-Seidel提出的交替最小化思想，即固定其他變量并交替變量x和z之間的最小化。式(4)可表示為：

(6)

(7)

式(6)和式(7)闡述了ADMM的基本思想。ADMM的關(guān)鍵思想是在x和z上使用單個(gè)Gauss-Seidel，而不是通常的聯(lián)合最小化。在每次迭代中，x和z以交替方式更新，然后在x更新前z更新后完成雙變量u更新。當(dāng)收斂條件滿足時(shí)迭代停止，然后x和z獲得最優(yōu)解。ADMM結(jié)合了乘子和交替最小化算法的優(yōu)點(diǎn)成為解決約束優(yōu)化問題的一個(gè)有效的框架，例如圖像重建。在接下來的章節(jié)中，將回顧基于交替方向乘子的凸優(yōu)化算法結(jié)構(gòu)[4]，以解決射電干涉成像中出現(xiàn)的凸優(yōu)化任務(wù)，并提出一種用于射電干涉成像改進(jìn)的基于交替方向乘子法的對(duì)偶前后向算法。

3 改進(jìn)的凸優(yōu)化射電干涉成像算法

3.1 基于交替方向乘子法的對(duì)偶前后向算法

交替方向乘子法是一種求解優(yōu)化問題的計(jì)算框架，但不具有任何內(nèi)在的并行結(jié)構(gòu)。為實(shí)現(xiàn)并行化處理，Onose等人提出了一種新的凸優(yōu)化算法—基于交替方向乘子法的對(duì)偶前后向算法。式(1)中的重構(gòu)任務(wù)可以使用指示函數(shù)作為凸最小化問題來重新定義：

=z,

(8)

其中，g(z)=iB(z)，B={z∈RM:‖Φx-y‖2≤ε}。f(x)=‖Ψ?x‖1+iC(x)，iC(x)是凸集C?RN的指示函數(shù)：

(9)

函數(shù)iC(x)為復(fù)原結(jié)果引入了正定性要求?！?x‖1表示先驗(yàn)稀疏性，g(z)=iB(z)是數(shù)據(jù)保真度項(xiàng)，該項(xiàng)約束殘差低于由噪聲水平ε定義的l2邊界。ADMM利用下面的增廣拉格朗日函數(shù)[8]：

(10)

式(10)對(duì)函數(shù)f和g的平滑性沒有明確的假設(shè)。在算法上，它們通過在每個(gè)原始變量x和z上的最小化之間交替來分割最小化步驟，接著通過梯度上升實(shí)現(xiàn)乘子λ的最大化。式(10)中x的最小化不接受閉合形式解[4]。因此，通過執(zhí)行梯度步驟和鄰近步驟來實(shí)現(xiàn)前向-后向迭代結(jié)構(gòu)。前向-后向迭代結(jié)構(gòu)通過線性化當(dāng)前點(diǎn)x(t)處的增廣拉格朗日的二次懲罰項(xiàng)并添加近鄰項(xiàng)來解決x上的最小化問題。近鄰ADMM的迭代更新方程由下式給出:

z(t+1)=proxγh(y-Φx(t)-λ(t)),

(11)

x(t+1)=proxμγf(x(t)-μΦH(λ(t)+Φx(t)-

y+z(t+1))),

(12)

λ(t+1)=λ(t)+β(Φx(t+1)-y+z(t+1)),

(13)

其中proxf為函數(shù)f的近鄰算子[9]，定義如下：

(14)

其中β和μ為合適的步長(zhǎng)。正向梯度步驟通過在與殘差的l2范數(shù)的梯度相反的方向上執(zhí)行遞減梯度步驟來實(shí)現(xiàn)，類似于CLEAN算法的主循環(huán)。隨后是向后步驟，將近似算子逼近非光滑函數(shù)f，類似于CLEAN的次循環(huán)[4]。它的解由固定點(diǎn)方程表征：

x(t+1)=

(15)

在給定的迭代中，它使用FB迭代執(zhí)行梯度步驟和在給定基Ψ?中執(zhí)行軟閾值操作。前向(梯度)步驟在于與殘差的l2范數(shù)的梯度相反的方向更新。它類似于CLEAN的一個(gè)主要循環(huán)。后向(軟閾值)步驟包括將閾值以上的所有Ψ?x系數(shù)的絕對(duì)值減小到閾值以上，并將閾值以下的那些系數(shù)設(shè)置為零。該步驟與CLEAN的次循環(huán)相似，其中軟閾值與環(huán)路增益因子類似。

3.2 改進(jìn)的基于交替方向乘子法的對(duì)偶前后向算法

3.2.1 基于交替方向乘子法的兩步法

為了提高重建高分辨率天文圖像的質(zhì)量，提出一種改進(jìn)的基于交替方向乘子法的對(duì)偶前后向法，它保留了ADMM的計(jì)算簡(jiǎn)便性，但重建質(zhì)量明顯提高。改進(jìn)算法稱之為IADMM。該算法的設(shè)計(jì)方法是仿照文獻(xiàn)[10]的思想，通過對(duì)式(15)中解x進(jìn)行二次加權(quán)更新，其定義如下：

(16)

在算法1中詳細(xì)描述了所得到的改進(jìn)算法。該算法由前向步驟和后向步驟組成。前向步驟對(duì)應(yīng)于梯度步驟，后向步驟通過近鄰算子執(zhí)行隱式子梯度步驟。整個(gè)算法流程為：

算法1IADMM

輸入：y,F,M,Z,G,Ψ,κ,ρ,d

輸出：重建圖像x

1：初始值x(0),h(0),g(0),s(0),d0

2：循環(huán)t=1,2,

3：p(t)=MFZx(t)

4：h(t)=Gρ(t)+g(t-1)

5：g(t)=g(t-1)+ρ(Gp(t)-h(t))

6：s(t)=G?(Gp(t)+h(t)-g(t))

10：直至收斂

11：Γ(Z,κ)

12：循環(huán)t=1,2,

13：初始值：q(0),z(0),η

14：Re=ηq(k-1)+Ψ?z(k-1)

16：z(k)=z(k-1)-Ψq(k)

17：直至收斂

18：返回z(k)

3.2.2 復(fù)雜度

4 仿真結(jié)果

(17)

4.1 參數(shù)選擇

4.2 結(jié)果與分析

實(shí)驗(yàn)通過從模擬不完整可見度復(fù)原測(cè)試圖像來評(píng)估IADMM的重建性能。在所有測(cè)試中，使用的測(cè)試圖像是M31galaxy的256×256圖像，W28 supernova remnant的1024×1024圖像和Cygnus A的477×1025圖像。亮度值在區(qū)間[0.01,1]，在圖1中以log10標(biāo)度顯示。

圖1 測(cè)試圖像

IADMM與所有其他基準(zhǔn)算法進(jìn)行比較的SNR結(jié)果如圖2所示，其顯示了IADMM、PD、ADMM和SDMM的SNR隨M31測(cè)試圖像的迭代次數(shù)而變化的情況。在這個(gè)測(cè)試中，輸入數(shù)據(jù)被分成4個(gè)塊，小波級(jí)別被設(shè)置為4，可見度設(shè)置為M=10N。可見度被高斯噪聲干擾產(chǎn)生20 dB的ISNR。參數(shù)κ為歸一化閾值并且控制收斂速度。當(dāng)κ=10-3到κ=10-5通常產(chǎn)生良好且一致的性能。實(shí)驗(yàn)結(jié)果表明，對(duì)于所有不同的迭代次數(shù)，IADMM都優(yōu)于其他方法。在達(dá)到收斂之前，隨著迭代次數(shù)的增加，SNR增加。如圖2所示，改進(jìn)算法IADMM產(chǎn)生的SNR至少在前期迭代中大約是SDMM的2倍。此外，對(duì)于所有次數(shù)的迭代，IADMM都比ADMM增益1～2 dB。值得注意的是，迭代次數(shù)對(duì)IADMM和ADMM的SNR變化率影響較小，這表明它們具有更好的魯棒性。實(shí)驗(yàn)表明，IADMM比現(xiàn)有技術(shù)的方法具有更好的重建質(zhì)量。

圖2 以M31為測(cè)試圖像的4個(gè)算法SNR隨迭代次數(shù)變化的函數(shù)

接下來繼續(xù)從不同的可見度M對(duì)IADMM、ADMM、PD和SDMM算法的性能進(jìn)行研究。M31和Cygnus A的重建結(jié)果分別如表1、2所示，它們分別顯示了當(dāng)所使用的可見度M的數(shù)量是10N，5N，2N，1N和0.5N所重建得到的SNR。實(shí)驗(yàn)中，輸入數(shù)據(jù)被分成16個(gè)塊，參數(shù)κ=10-3。結(jié)果表明，IADMM在所有不同的可見度方面均優(yōu)于其他方法。如表1、2所示，PD和SDMM在M31和Cygnus A中顯示出相似的重建結(jié)果。值得注意的是，表1中M31的結(jié)果顯示IADMM的SNR有明顯提高，對(duì)于較大的可見度覆蓋范圍，比ADMM提高了2 dB的增益，同時(shí)又比PD和SDMM提高1 dB左右。隨著可見性數(shù)量的增加，改進(jìn)算法的優(yōu)點(diǎn)更加明顯。對(duì)于更復(fù)雜的圖像Cygnus A，改進(jìn)算法的優(yōu)勢(shì)更明顯，在表2中結(jié)果顯示IADMM的SNR ADMM提高約3 dB。此外，與PD和SDMM相比，改進(jìn)后的算法在所有可見度下均可實(shí)現(xiàn)約2 dB的增益。值得注意的是，增加更多的數(shù)據(jù)能提高重建的SNR，因?yàn)樵肼暺骄酶谩５?，?dāng)增加更可見度數(shù)據(jù)時(shí)，SNR增益會(huì)稍微停滯，主要是因?yàn)轭l率平面中的某些點(diǎn)仍未覆蓋。

表1 以M31為測(cè)試圖像的4個(gè)算法SNR在不同可見度M的結(jié)果

最后，從視覺評(píng)估上比較IADMM與ADMM的重建性能。在圖3中，(a)和(b)顯示的分別是IADMM算法下重建圖像和誤差圖像，(c)和(d)顯示的分別是ADMM算法下重建圖像和誤差圖像。圖3給出對(duì)于M=N，κ=10-3和輸入信噪比為20 dB的W28的重建結(jié)果：IADMM(28.82 dB)和ADMM(26.74 dB)。圖3的結(jié)果表明，IADMM實(shí)現(xiàn)了更好的SNR，與其他方法相比，IADMM仍然提高達(dá)到了2 dB左右增益。很顯然，2種方法都能為W28帶來好復(fù)原效果，都能夠恢復(fù)光源周圍的模糊區(qū)域。從視覺上看，IADMM產(chǎn)生的重建圖像在背景區(qū)域中具有較少的偽影，并且在結(jié)構(gòu)化的內(nèi)部區(qū)域中比在其他方法中具有較小的誤差。

表2 以Cygnus A為測(cè)試圖像的4個(gè)算法SNR在不同可見度M的結(jié)果

圖3 W28的重建結(jié)果

5 結(jié)束語(yǔ)

針對(duì)射電干涉成像中產(chǎn)生的凸優(yōu)化任務(wù)，提出了一種改進(jìn)的射電干涉成像凸優(yōu)化算法?；贏DMM算法高度可并行化的結(jié)構(gòu)，在ADMM的正向(梯度)步驟中通過對(duì)復(fù)原結(jié)果引入加權(quán)二次更新。其設(shè)計(jì)思想為更新方程取決于前兩個(gè)結(jié)果估計(jì)，而不僅僅是前一個(gè)。整個(gè)算法并行地執(zhí)行近鄰和梯度更新，其可以看作由在多數(shù)據(jù)、先驗(yàn)和圖像空間中并行運(yùn)行近鄰分裂和FB更新組成。通過u-v覆蓋模擬和使用不同星系的可見度來研究改進(jìn)算法。實(shí)驗(yàn)結(jié)果表明，改進(jìn)后的算法不僅可保持與ADMM相同的并行實(shí)現(xiàn)結(jié)構(gòu)和相同的計(jì)算負(fù)擔(dān)，而且可實(shí)現(xiàn)更好的噪聲圖像重構(gòu)。