潘潤寧

摘? 要： 針對射電天文圖像重建中由于越來越大的數(shù)據(jù)集產(chǎn)生的保證與提高圖像重建質(zhì)量的問題，提出一種基于壓縮感知技術(shù)的射電干涉圖像凸優(yōu)化重構(gòu)算法。該算法對原始對偶算法引入迭代收縮算子和閾值模型，迭代收縮算子在每一步迭代過程中利用前兩點(diǎn)進(jìn)行組合更新，迭代過程更為精確，通過閾值模型克服不連續(xù)性產(chǎn)生的偽吉布斯效應(yīng)，使得重建射電圖像光滑并保留了有用高頻信息。實(shí)驗(yàn)結(jié)果表明，改進(jìn)算法在信噪比與重建精度方面優(yōu)于PD算法，得到了更好的效果。

關(guān)鍵詞： 射電圖像重建; 壓縮感知; 原始對偶算法; 射電天文觀測; 凸優(yōu)化算法; 閾值去噪

中圖分類號： TN911.73?34? ? ? ? ? ? ? ? ? ? ? ?文獻(xiàn)標(biāo)識碼： A? ? ? ? ? ? ? ? ? ? ? ? ?文章編號： 1004?373X（2020）09?0046?05

An improved radio image reconstruction algorithm based on compressed sensing

PAN Running

（MOE Key Laboratory of Cognitive Radio and Information Processing， Guilin University of Electronic Technology， Guilin 541004， China）

Abstract： For the poor image quality due to growth of data sets increasingly in the radio astronomical image reconstruction， a radio interference image convex optimization reconstruction algorithm based on the compressed sensing technology is proposed to ensure and improve the image quality. For the proposed algorithm， an iterative shrinkage operator and a threshold model are introduced into the primal dual algorithm. On the basis of the first two points， the iterative shrinkage operator implements combination update in each step of iteration to obtain more accurate iteration process， and the pseudo?Gibbs effect caused by discontinuity is overcome by the threshold model， which makes the reconstructed radio image smooth and useful high?frequency information retained. The experimental results show that the improved algorithm is superior to the PD algorithm in terms of signal?to?noise ratio （SNR） and reconstruction accuracy， and achieves a better effect.

Keywords： radio image reconstruction; compressed sensing; primal dual algorithm; radio astronomical observation;? convex optimization algorithm; threshold denoising

0? 引? 言

宇宙一直是人類進(jìn)行探索的重要目標(biāo)，從探知某一顆星體地表的結(jié)構(gòu)，以至于到外太空之中是否存在著其他的生命文明，是否擁有可以被人類利用的物資資源等，這些信息可以由天文圖像獲得。到目前為止，射電天文圖像重建中比較成熟的技術(shù)路線采用的算法大部分是以潔化算法[1]和最大熵算法[2]作為核心思想而進(jìn)行改進(jìn)的算法，而這兩種算法并沒有在理論和實(shí)際上解決欠采樣問題。文獻(xiàn)[3]發(fā)表的壓縮感知理論為解決該問題提供了可能性。壓縮感知引入了一種新的優(yōu)于傳統(tǒng)奈奎斯特采樣的信號重建框架，相比于傳統(tǒng)成像方法，在計算速度、靈活性和數(shù)據(jù)保真度方面有較大優(yōu)勢[4]。文獻(xiàn)[5]利用壓縮感知框架研究宇宙射線成像算法，比傳統(tǒng)算法獲得了更好的性能。文獻(xiàn)[6]提出基于壓縮感知的SparseRI，多功能圖像重構(gòu)技術(shù)。文獻(xiàn)[7]提出一種基于壓縮感知理論新的解卷積的方法，用于恢復(fù)點(diǎn)源或展源圖像，獲得比傳統(tǒng)算法更好的性能。文獻(xiàn)[8]提出一種用于射電圖像重建的原始對偶算法（PD），通過原始問題和對偶問題間的迭代得到重建結(jié)果。本文在基于PD算法的基礎(chǔ)上，提出一種改進(jìn)的算法，對射電圖像得到了更好的重建質(zhì)量，提高了重建圖像的精度。

1? 射電干涉成像

射電干涉測量天體是由一系列的陣列望遠(yuǎn)鏡來虛擬綜合等效口徑及高角分辨率進(jìn)行實(shí)現(xiàn)的。陣列望遠(yuǎn)鏡測量天空亮度的強(qiáng)度數(shù)據(jù)，即射電干涉的能見度數(shù)據(jù)是在空間頻率域中。Van?Cittert?Zernike理論指出了可見度[V（u，v）]和源亮度函數(shù)[I（l，m）]之間的關(guān)系：

[V（u，v）=I（l，m）e-i2π（ul+vm）dldm] （1）

式中：[（u，v）]是空間頻率域的坐標(biāo);[（l，m）]是圖像域的坐標(biāo)。但是在現(xiàn)實(shí)測量中，采樣是不完全的。對源亮度函數(shù)利用采樣函數(shù)[S（u，v）]進(jìn)行[u?v]覆蓋后，得到可見度函數(shù)[V]與測量亮度[I（l，m）]的關(guān)系：

[I（l，m）=S（u，v）V（u，v）e-i2π（ul+vm）dudv]? （2）

通過離散重寫式（2），得到一個線性測量模型：

[y=Φx+n]? （3）

式中：[x]為天空亮度;[y]為測量值;測量算子[Φ∈C?M×N]，映射圖像域與[u?v]能見度空間;[n]為噪聲。

由于望遠(yuǎn)鏡的物理位置以及采樣方案的局限性，使得式（3）成為一個不適定逆問題。測量算子[Φ=GMFZ]，其中，[G]為卷積差值算子矩陣，引入分塊矩陣[M]，[F]為過采樣傅里葉算子，[Z]為預(yù)補(bǔ)償差值縮放算子。壓縮感知理論指出，絕大多數(shù)信號通過一個字典[Ψ]能夠有稀疏表示，[Ψ?x]只包含有少量的非零元素。從而可以通過解決以下不適定逆問題恢復(fù)信號[x]：

[minxΨ?x1s.t.? ? y-Φx2≤ε]? ?（4）

式中：[Ψ?]為[Ψ]的伴隨算子;[ε]為噪聲[l2]范數(shù)的上界。在此框架下，采用SDMM算法能夠?qū)?shù)據(jù)進(jìn)行分布式處理，但每次迭代都需要求解線性方程組，這樣使得恢復(fù)大圖像時需要的計算量很大。相比之下，PD算法的梯度算子、接近算子以及線性算子都可以單獨(dú)使用，這與其他方案相比，在計算上具有很大的優(yōu)勢。

2? 原始?對偶算法基本思想

現(xiàn)在引入求解約束優(yōu)化問題的方法，即原始?對偶算法，其思想是通過對偶問題來對原始問題進(jìn)行求解，就必須使得原始問題與對偶問題的最優(yōu)值相等?？紤]如下約束最優(yōu)化問題：

[min f0（x）s.t.? ? hi（x）=0，? ? i=1，2，…，m? ? ? ? ? gj（x）≤0，? ? j=1，2，…，n] （5）

式（5）稱為約束最優(yōu)化問題的原始問題，將其構(gòu)建為拉格朗日函數(shù)如下：

[L（x，λ，v）=f0+i=1mλifi（x）+i=1nvihi（x）]? ?（6）

那么，在滿足約束條件下有：

[minxP（x）=minxmaxλ，vL（x，λ，v）=minxf（x）] （7）

即原始的優(yōu)化問題與[minxP（x）]等價，則原始問題的最優(yōu)值為：

[p*=minxP（x）]? ?（8）

然后考慮原始問題的對偶問題如下：

[maxλ，vD（λ，v）=maxλ，vminxL（x，λ，v）]? ? ?（9）

定義對偶問題的最優(yōu)解為：

[d*=maxλ，vD（λ，v）] （10）

當(dāng)原始問題和對偶問題都有最優(yōu)解時，得到：

[d*=maxλ，vminxL（x，λ，v）≤minxmaxλ，vL（x，λ，v）=p*]? （11）

當(dāng)原始問題的最優(yōu)值與對偶問題的最優(yōu)值相等，即[p*=d*]時，就可以用求解對偶問題來對原始問題進(jìn)行求解。

3? 改進(jìn)的射電圖像重建算法

3.1? 射電干涉成像的原始?對偶問題

文獻(xiàn)[8]提出了一種新的凸優(yōu)化算法，將凸優(yōu)化算法引入到射電干涉成像之中，在壓縮感知的框架下，可以將成像任務(wù)寫成如下形式：

[minxf（x）+l（Ψ?x）+h（Φx）]? （12）

式中[f=ιC（z）]，[ιC（z）]是凸集[C?RN]的指示函數(shù)：

[ιC（z）=0，? ? z∈C+∞，? ? z?C]? （13）

[h（z）=ιB（z），B=z∈C?M：z-y2≤ε]表示數(shù)據(jù)保真度項，該項對殘差約束為低于由噪聲水平[ε]定義的[l2]范數(shù)邊界，[l=?1]表示給定的完備字典[Ψ]的稀疏先驗(yàn)性。

將原始對偶算法的思想應(yīng)用到式（12）定義的最小化任務(wù)中，得到原始問題和對偶問題如下：

[minxf（x）+γl（Ψ?x）+h（Φx）]? ?（14）

[minu，vf*（（-Ψu-Φ?v）+1γl*（u）+h*（v））] （15）

式中：[γ]是一個附加的調(diào)整參數(shù);[u]，[v]是對偶變量，算法收斂到一個庫恩塔克點(diǎn)[（x，u，v）]，[x]即為原始問題的解，[（u，v）]則為對偶問題的解。[γ]的取值不會影響優(yōu)化問題。解決成像問題即在原始問題式（14）與對偶問題式（15）之間迭代交替求出最優(yōu)解。

3.2? 改進(jìn)的原始?對偶算法

3.2.1? 迭代收縮算子

具有向前?向后迭代結(jié)構(gòu)[9]的原始對偶算法框架，對于解決式（14）和式（15）具有較好的并行與可擴(kuò)展性，在存儲器需求和每次迭代的計算負(fù)載方面有較大的靈活性，但成像精確度不足，并且在低迭代次數(shù)時算法的重建結(jié)果信噪比較差。為了提高大視場下高分辨率天文圖像的成像質(zhì)量，針對存在的缺陷問題，在PD算法[8]基礎(chǔ)上提出改進(jìn)算法1。在提升圖像重建結(jié)果信噪比的同時，考慮低迭代次數(shù)時信噪比過差的問題。

算法1

1.? 輸入[x（0），x（0），v（0），u（0），v（0），α，τ，κ，?，σ]

2.? 重復(fù)[t=1，2，…]

3.? [p（t）=MFZx（t-1）]

4.? [v（t）=（I-ΡB）（v（t-1）+Gp（t））]

5.? [v（t）=v（t-1）+α（v（t）-vt-1）]

6.? [v（t）=G?v（t）]

7.? [u（t）=（I-ΓκΨs）（u（t-1）+Ψ?x（t-1））]

8.? [u（t）=u（t-1）+α（u-u（t-1））]

9.? [u（t）=Ψu（t）]

10. [x（t）=PC（x（t-1）-τ（Z?F??M?v（t）+σu（t）））]

11. [x（t）=x（t）+tk-1tk（x（t）-x（t-1））]

12. [x（t）=x（t-1）+α（x（t）-x（t-1））]

13. [x=2x（t）-x（t-1）]

14. 直至收斂

根據(jù)FISTA[10]技術(shù)，引入迭代收縮算子：

[tk+1=λ1λk+1+λ21λ2k+1+4t2k2，? ? ? ? ? ? ?t1=λ1， （λk）k≥1?ε，1， ε∈0，1] （16）

式中[λ]為可調(diào)節(jié)參數(shù)[11]，可以使得迭代收縮算子擁有更大的靈活性。

步驟10，步驟11中運(yùn)用迭代收縮算子，不單是通過前一次迭代求得函數(shù)最小點(diǎn)[xk-1]，而且在迭代過程中利用前兩點(diǎn)[xk，xk-1]進(jìn)行特定的線性組合更新得出[x]，由此每次得到的迭代結(jié)果更為準(zhǔn)確。

3.2.2? 改進(jìn)的閾值模型

引入文獻(xiàn)[12]中的閾值模型對原算法進(jìn)行改進(jìn)。由于軟閾值模型構(gòu)造的稀疏先驗(yàn)函數(shù)存在原始信號估計不準(zhǔn)確的問題，[Ψ?x]的系數(shù)由固定的閾值進(jìn)行收縮，具有恒定誤差[13]，會丟失一些有用的高頻信息，可能使得重構(gòu)圖像存在邊緣模糊等現(xiàn)象，損失圖像的細(xì)節(jié)，也可能導(dǎo)致重建系數(shù)和原始系數(shù)間產(chǎn)生偏差。鑒于軟閾值自身所具有的缺陷，使用改進(jìn)的閾值模型為：

[ΓT（z）=sgn（z）z-Texpz-Ta，? ? z≥T0，? ? z

式中：[a]是一個調(diào)節(jié)因子，使得該閾值模型具有靈活可調(diào)性，針對含噪信號進(jìn)行改變，可以獲得去噪效果更優(yōu)的閾值函數(shù)曲線;[sgn（）]是符號函數(shù)。新的閾值模型能夠克服不連續(xù)性產(chǎn)生的偽吉布斯效應(yīng)，使得重構(gòu)圖像更加光滑，有更好的去噪效果。

改進(jìn)的原始?對偶算法如算法1所示。所有的對偶變量在步驟4～步驟9中并行更新。步驟4通過應(yīng)用共軛函數(shù)[h*]的接近算子，將殘差約束到[l2-ball]上以實(shí)現(xiàn)數(shù)據(jù)保真度，這就接受了一個封閉形式的解投影到[l2-ball]上：

[PB（z）?εz-yz-y2+y，? ? z-y2>εz，? ? z-y2≤ε]? ? （18）

步驟10中通過對偶變量[v（t），u（t）]對原始變量進(jìn)行類子梯度更新，投影到凸集[C]上：

[PC（z）??（z），? ? ?（z）>00，? ? ?（z）≤0]? ?（19）

然后送到步驟11中用迭代收縮算子更新，最后在步驟12中得到[x（t）]。其中，[κ]為無標(biāo)度參數(shù)，[?，σ，τ]為更新步長，[?=1Φ2S]，[σ=1Ψ2S]，[τ=0.49]，[α]為松弛因子。

4? 仿真實(shí)驗(yàn)

為了測試改進(jìn)算法對射電天文圖像的重構(gòu)性能，下面通過實(shí)驗(yàn)仿真進(jìn)行評估，并與其他三種算法（交替方向乘子法（ADMM）、同時同向乘子法（SDMM）以及原始對偶算法（PD））進(jìn)行比較。實(shí)驗(yàn)中，使用高斯采樣隨機(jī)生成[u?v]覆蓋，為了模擬RI中通常存在的不完全采樣，使用逆高斯分布在高頻數(shù)據(jù)中引入空穴。稀疏先驗(yàn)使用SARA小波集合[14]，基集合分為[nd=9]個單獨(dú)基，第一個是狄拉克基，其他8個基為Daubechies小波的前8個Db1～Db8。小波等級設(shè)置[n=4]。插值內(nèi)核矩陣[G]使用[8×8]Kaiser?Bessel進(jìn)行插值，平均近似均勻分布的頻率值，以估計與每個可見度相關(guān)的頻率處的值。過采樣因子設(shè)置為4。由于宇宙噪聲服從高斯分布，實(shí)驗(yàn)中通過引入高斯白噪聲對原始數(shù)據(jù)進(jìn)行干擾，輸入信噪比為[PISNR=20logy02n2]，[y0]為潔凈的測量矢量，計算可見度[M]下進(jìn)行圖像重構(gòu)的平均信噪比。算法的重構(gòu)性能由信噪比，即真實(shí)天空的標(biāo)準(zhǔn)偏差[σx]與估計模型的標(biāo)準(zhǔn)偏差[σx-x]的比值來定義：

[PSNR=20lgσxσx-x]? ?（20）

對于實(shí)驗(yàn)中的算法，其終止標(biāo)準(zhǔn)為全局殘差范數(shù)[j=1ndyi-Φjx（t）22]低于閾值[ε2=（2M+34M）σ2] ，或迭代滿足[l2]約束條件，并且[δ≤δ]解范數(shù)變化很小時終止，其中，[δ=x（t）-x（t-1）2x（t）2]。

實(shí)驗(yàn)中測試用圖如圖1所示。其中，以[lg]為標(biāo)度。

首先研究在不同可見度[M]下本文算法以及其他三種算法的重構(gòu)性能。以M31 galaxy圖像為例，如表1所示，輸入數(shù)據(jù)分為16塊，ISNR為20 dB，設(shè)置最大迭代次數(shù)為100次，得到可見度[M]分別在10[N]，5[N]，2[N]，1[N]時重建結(jié)果的信噪比。結(jié)果表示在各個可見度設(shè)置下，本文算法均優(yōu)于PD，ADMM，SDMM三種算法的重構(gòu)質(zhì)量?？梢姸葹?[N]時，重構(gòu)得到的SNR優(yōu)于其他算法1 dB左右;可見度為10[N]時，重構(gòu)后的SNR比其余三種算法約有2 dB左右的增益，隨著可見度數(shù)據(jù)的增加，優(yōu)勢更加明顯。

接下來研究本文算法隨著重構(gòu)迭代次數(shù)變化得到的SNR結(jié)果的變化情況。以M31 galaxy和Cygnus A圖像為例，輸入數(shù)據(jù)分為16塊，ISNR為20 dB，可見度[M=][10N]。如圖2所示，本文算法在迭代穩(wěn)定后，信噪比均高于PD，ADMM，SDMM三種算法，能夠得到更好的射電圖像質(zhì)量。在較低迭代次數(shù)時，較之改進(jìn)前的PD算法，信噪比有更顯著的提升，并且約為SMDD算法的2倍。

重建結(jié)果如圖3所示。以W28 supernova remnant為例，圖3a）為本文算法重建結(jié)果（29.36 dB）;圖3b）為PD算法重建結(jié)果（27.12 dB）。參數(shù)設(shè)置在可見度[M=2N]，輸入數(shù)據(jù)分為16塊，ISNR為20 dB。兩個重建效果圖顯示本文算法較PD算法獲得的質(zhì)量更高，圖像更清晰，邊緣表現(xiàn)更為細(xì)膩，能保留更多的高頻信息，在信噪比上提高了2 dB，有較少偽影。而PD算法則丟失了較多圖像細(xì)節(jié)，光源外沿模糊。本文算法在使邊緣輪廓較PD算法更為清晰的基礎(chǔ)上，能保留更多的射電圖像細(xì)節(jié)，在結(jié)構(gòu)內(nèi)部區(qū)域中也擁有更好的表現(xiàn)。

5? 結(jié)? 語

為了提高射電圖像重建質(zhì)量，本文提出了一種改進(jìn)的基于壓縮感知的射電圖像重建算法。該算法在原始對偶算法的基礎(chǔ)上，利用迭代收縮算子和閾值模型精確重構(gòu)原始信號并提高算法魯棒性。實(shí)驗(yàn)仿真結(jié)果表明，本文算法有效地提高了重構(gòu)圖像的精度和質(zhì)量。

參考文獻(xiàn)

[1] H?GBOM J A. Aperture synthesis with a non?regular distribution of interferometer baselines [J]. Astronomy & astrophysics supplement， 2011， 15（15）： 417.

[2] FOCSA A， TOMA S A， DATCU M. Maximum entropy image reconstruction applied to C?band ground based synthetic aperture radar [C]// 2017 IEEE International Geoscience and Remote Sensing Symposium （IGARSS）. Fort Worth， TX， US： IEEE， 2017： 3437?3440.

[3] CANDES E J， ROMBERG J， TAO T. Robust uncertainty principles： exact signal reconstruction from highly incomplete frequency information [J]. IEEE transactions on information theory， 2006， 52（2）： 489?509.

[4] WIAUX Y， JACQUES L， PUY G， et al. Compressed sensing imaging techniques for radio interferometry [J]. Monthly notices of the Royal Astronomical Society， 2017， 395（3）： 1733?1742.

[5] WIAUX Y， PUY G， VANDERGHEYNST P. Compressed sen?sing reconstruction of a string signal from interferometric observations of the cosmic microwave background [J]. Monthly notices of the Royal Astronomical Society， 2010， 402（4）： 2626?2636.

[6] WENGER S， MAGNOR M， PIHLSTR?M Y， et al. SparseRI： a compressed sensing framework for aperture synthesis imaging in radio astronomy [J]. Publications of the Astronomical Society of the Pacific， 2010， 122（897）： 1367?1374.

[7] LI F， BROWN S， CORNWELL T J， et al. The application of compressive sampling to radio astronomy： II. Faraday rotation measure synthesis [J]. Astronomy & astrophysics， 2012， 531（2）： 1261?1268.

[8] ONOSE A， CARRILLO R E， REPETTT A， et al. Scalable splitting algorithms for big?data interferometric imaging in the SKA era [J]. Monthly notices of the Royal Astronomical Society， 2016， 462（4）： 4314?4335.

[9] KOMODAKIS N， PESQUET J C. Playing with duality： an overview of recent primal?dual approaches for solving large?scale optimization problems [J]. IEEE signal processing magazine， 2015， 32（6）： 31?54.

[10] CHEN Shaoli， YANG Min. An adaptive fast iterative shrin?kage threshold algorithm [C]// 2017 29th Chinese Control and Decision Conference （CCDC）. Chongqing， China： IEEE， 2017： 2190?2194.

[11] ZIBETTI M V W， HELOU E S， PIPA D R. Accelerating over?relaxed and monotone fast iterative shrinkage?thresholding algorithms with line search for sparse reconstructions [J]. IEEE transactions on image processing， 2017， 26（7）： 3569?3578.

[12] ZHANG Songjie， MA Yuliang， ZHANG Qizhong， et al. Motor imagery electroencephalogram de?noising method based on EEMD and improved wavelet threshold [C]// 2018 Chinese Control And Decision Conference （CCDC）. Shenyang， China： IEEE， 2018： 5589?5594.

[13] 易清明，陳明敏，石敏.一種改進(jìn)的小波去噪方法在紅外圖像中應(yīng)用[J].計算機(jī)工程與應(yīng)用，2016，52（1）：173?177.

[14] CARRILLO R E， MCEWEN J D， WIAUX Y. Sparsity avera?ging reweighted analysis （SARA）： a novel algorithm for radio?interferometric imaging [J]. Monthly notices of the Royal Astronomical Society， 2012， 426（2）： 1223?1234.