溫 珊， 王俊義， 勞保強

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

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

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

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

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

(1)

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

2 交替方向乘子法

為了解決式(1)中的不適定逆問題，引入Boyd等[8]提出的用于求解約束優(yōu)化問題的方法——交替方向乘子法(ADMM)。ADMM是旨在將對偶上升的可分解性與乘子方法的優(yōu)越收斂性相結合的一種算法。該算法通過式(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)組成最小化的目標函數(shù)。Ax+Bz=c是一個最小化約束。增廣的拉格朗日表達式為

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

(3)

(4)

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

(5)

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

(6)

(7)

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

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

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

交替方向乘子法是一種求解優(yōu)化問題的計算框架，但不具有任何內在的并行結構。為實現(xiàn)并行化處理，Onose等人提出了一種新的凸優(yōu)化算法—基于交替方向乘子法的對偶前后向算法。式(1)中的重構任務可以使用指示函數(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)為復原結果引入了正定性要求?！?x‖1表示先驗稀疏性，g(z)=iB(z)是數(shù)據(jù)保真度項，該項約束殘差低于由噪聲水平ε定義的l2邊界。ADMM利用下面的增廣拉格朗日函數(shù)[8]：

(10)

式(10)對函數(shù)f和g的平滑性沒有明確的假設。在算法上，它們通過在每個原始變量x和z上的最小化之間交替來分割最小化步驟，接著通過梯度上升實現(xiàn)乘子λ的最大化。式(10)中x的最小化不接受閉合形式解[4]。因此，通過執(zhí)行梯度步驟和鄰近步驟來實現(xiàn)前向-后向迭代結構。前向-后向迭代結構通過線性化當前點x(t)處的增廣拉格朗日的二次懲罰項并添加近鄰項來解決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)

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

x(t+1)=

(15)

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

3.2 改進的基于交替方向乘子法的對偶前后向算法

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

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

(16)

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

算法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 復雜度

4 仿真結果

(17)

4.1 參數(shù)選擇

4.2 結果與分析

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

圖1 測試圖像

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

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

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

表1 以M31為測試圖像的4個算法SNR在不同可見度M的結果

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

表2 以Cygnus A為測試圖像的4個算法SNR在不同可見度M的結果

圖3 W28的重建結果

5 結束語

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