射電天文圖像恢復的改進方法*

衛 明， 袁梅宇， 張秋明， 王 赟

(1.昆明理工大學 信息與自動化學院 ，云南 昆明 650000；2.陽泉市高新創業服務中心，山西 陽泉 045000)

0 引 言

目前,在射電天文圖像重構算法中成熟的技術和方案大多采用潔化算法[1]和最大熵方法[2]作為核心思想的改進方法。潔化算法作為處理點源圖像時恢復效果好于最大熵方法；而在處理展源時，最大熵方法優于潔化算法。很多天文圖像處理軟件包均采用這兩種算法，但對于欠采樣問題兩者并未從根本上解決。

壓縮感知理論是近年來由Candes，Romberg，Tao等人[3]新發表的關于信號處理的一種新的稀疏信號重構技術，在射電天文成像系統中對于采樣函數是稀疏的。Donoho D L等人[4]曾多次將壓縮感知技術應用于廣域射電望遠鏡干涉成像之中，并用該方法解決了宇宙微波背景中的宇宙射線成像問題。為保證所成圖像的精確度，本文使用日像儀天線接收太陽電磁波得到數據并對成像觀測數據仿真，通過對輸入圖像傅里葉變換與對應位置相乘處理得到該圖的頻域數據，再經過傅里葉變換得到了模糊圖像。 壓縮感知的采樣過程在稀疏域上進行，并且采樣的點數很少；在重構方面，利用壓縮感知理論中應用的最小范數，在大量實驗中展示了其良好的效果。利用最大熵方法恢復了圖像。本文方法取得了不錯的天文圖像恢復重構效果。

1 壓縮感知基礎理論

1.1 信號稀疏表示

(1)

x為只有K個基向量的一個線性組合，并且滿足K?N。K為式(1)中si的非零個數，剩下的(N-K)個分量均為零。由于信號的稀疏特性，則式(1)中存在標準正交基 ，使得其系數只有很少的大系數和很多的小系數?？蓧嚎s的信號可以很好地被K稀疏信號逼近[8]。

獲得全部的N長度的可壓縮采樣信號X；基于s=ΨTX，計算在完備集下的變換系數{si}；確定K個大系數分量，丟棄(K-N)小系數；對K個大系數分量進行編碼，并記錄下該K個大系數分量的位置信息以便于重構信號。

這種采樣然后壓縮的構架模式本身存在著一些缺陷：1)采集到的數據長度N很大，但最終所確定的K個大系數數目可能會非常?。?)編碼器必須計算出所有的N長度的變換系數{si}，保留K個最大的，丟棄其余的；3)編碼器不但要記錄K個大系數，還要記錄其位置信息。

y=ΦX=ΦΨs=Θs

(2)

式中Θ=ΦΨ為一個大小為M×N的矩陣。

觀測過程是非自適應的，即觀測矩陣Φ的選取與信號X不相關。

設計一個觀測矩陣Φ和一個重構算法，要求：觀測矩陣是一個穩定的矩陣，使得在信號降維過程中不會損失重要的信息；重構算法是從y=Θs中，在已知y和Θ的情況下找到滿足y=Θs的最稀疏的s。

1.2 穩定觀測矩陣設計

為解決已知K個非零值的問題。線性方程組M×N有唯一解的充分條件是：對于任意一個具有與s中K個非零元素完全相同的N維向量ν，滿足

(3)

從式(3)中可以看出，矩陣Θ要必須要保證K稀疏向量的度量尺度長度在一定范圍內變化。

在實際運用中并不清楚列向量s中K個非零元素的具體位置。然而，對于K稀疏或者可壓縮信號有解的充分且必要的條件正是式(3)，即約束等距條件(restricted isometry property,RIP)[9];另外，一種保證穩定性的方法是保證觀測矩陣Φ和稀疏標準正交基Ψ是非相干的，即向量{φj}不能稀疏表示向量{ψj}，反之亦然。

在壓縮感知理論中，通常會為了避免關于選擇具有約束等距條件的矩陣，而選擇一個隨機矩陣來作為觀測矩陣Φ。如矩陣中每個元素φj,i均互相獨立且遵循同一分布的隨機變量(independent and identically distributed,IID)，比如，觀測矩陣Φ的元素均服從期望為0，方差為1/N的高斯分布[10]，即可保證觀測信號y是原始信號X的M個不同的隨機加權的線性組合。另外，矩陣元素為±1的隨機Rademacher矩陣也滿足約束等距條件，具有普遍性[11]。

1.3 信號重構算法

為了在轉換零空間中，找出原始信號的稀疏解向量s：

最小l2范數重構方法

(4)

由于式(4)中的l2范數不能反映出解的稀疏特性，因此,直觀的思想在解向量中直接找到最稀疏的那組解，即最小l0范數

(5)

由式(5)中可以看出，只用M=K+1個高斯觀測值，任意的優化解都能夠高概率直接重構K稀疏信號[12]。但求解式(5)是一個不穩定，而且NP難問題。

2 最大熵原理

圖像熵的定義有很多形式[13～16],其中,Frieden提出的形式主要用于天文圖像上。熵最大時表示這些事件具有最大的不確定度，這在客觀上要求事件是等概率分布,使圖像的灰度趨于平滑,圖像最大熵的恢復問題的公式表達為

(6)

由圖像熵的定義可知，I(x,y)不能為負,否則H(x,y)沒有意義，因此，使圖像最大熵作為判據，圖像的正性將這個先驗知識自然地得到了保證。

3 壓縮感知與最大熵結合的圖像恢復方法

3.1 最大熵方法

(7)

(8)

式中f(I)為I(x,y)的熵函數。對式(9)I(x,y)求導有

(9)

滿足式(9)的解取代了式(7)的約束條件，獲得了關于λ(x,y)的非線性等式。稱函數J(x,y)為“帶限”函數，通過式(9)能找到滿足熵值最大的圖像I(x,y)，即

I(x,y)=f'-1[J(x,y)]≡g[J(x,y)]

(10)

熵函數的定義為f(I)=-IlnI，則g[J(x,y) ]=e(-1-J)。因此，圖像函數I(x,y)即為一個關于“帶限”函數J(x,y)的非線性函數。

3.2 最大熵方法與壓縮感知相結合的方法理論

討論一維熵函數

H=∑-xlnx

(11)

由函數-xlnx在x∈(0,1]區間的圖像知曲線光滑，且在(0,+∞)區間內可導，在x=e-1時-xlnx達到最大值。

為了使函數滿足凸函數，將熵函數變形為xlnx，求其導數為1+lnx，其在x∈(0,5]區間取值示意圖如圖1所示。圖像中1+lnx只含有一個零點并且xlnx在零點左側是減函數，零點右側是增函數，因此，函數xlnx具有全局最小值。

圖1 函數 1+lnx的曲線圖像

用泰勒公式將xlnx展開在x0=0.4展開，有

(12)

圖2 函數y與y1對比示意

取得區間為[0.01,1]步長為0.01的MATLAB仿真圖。兩曲線在該區間下的均方誤差為0.053 7。求解最大熵就轉換為了求解一個二次規劃問題，即

(13)

因此，最大熵條件是作為一個目標函數，與壓縮感知方法相結合后的目標函數變為

(14)

4 實 驗

4.1 最大熵與壓縮感知相結合方法實驗

實驗分為3幅不同的輸入圖片和3個不同的采樣函數，即其效果和評價如圖3和表1～表3,即采樣函數為22條線，大型螺旋線、中心螺旋線，仿真實驗結果分別為峰值信噪比(PSNR)和結構相似度(SSIM)。

表1 采樣函數為22條線仿真實驗結果

表2 采樣函數為大型螺旋線實驗結果

圖3 采樣函數為中心螺旋線仿真實驗結果

圖像臟圖PSNR/dB恢復圖PSNR/dB臟圖SSIM/%恢復圖SSIM/%Sun117.364217.596871.4172.39Sun223.039923.215998.8399.04Sun319.902219.356394.0994.49

4.2 最大熵方法和最大熵與壓縮感知結合方法對比

對于兩種方法的對比結果如表4～表6。

表4 采樣函數為22條線2種算法的恢復效果數據對比

表5 采樣函數為大型螺旋線2種算法的恢復效果數據對比

表6 采樣函數為中心螺旋線2種算法的恢復效果數據對比

4.3 實驗結果分析

總體上看,用壓縮感知和最大熵相結合的方法恢復圖像的效果無論是PSNR還是SSIM上均較最大熵方法略勝一籌。不足之處在于,對于最后一組實驗即183點中心螺旋線構成的采樣函數中,sun1和sun3圖像的PSNR和SSIM均較傳統的最大熵恢復效果差，主要是因為采樣點數太少，采樣率為0.279 %，采樣條件比較苛刻，另外，又由于采樣點主要集中在零頻附近幾乎所有的采樣點均為低頻采樣。

5 結 論

應用最大熵方法和壓縮感知方法，對熵函數進行變換，用泰勒公式展開,化為二次型的形式，然后將這個二次型與稀疏解最小化作為目標函數，對圖像進行恢復，對比實驗數據，圖像的效果有所提高，證明了方法的有效性。

參考文獻:

[1] Monnier J D.An introduction to closure phases[C]∥Principles of Long Baseline Stellar Interferometry,2000:203.

[2] H?gbom J A.Aperture synthesis with a non-regular distribution of interferometer baselines[J].Astronomy & Astrophysics Supplement,1974,15(15):417.

[3] Narayan R,Nityananda R.Maximum entropy image restoration in astronomy[J].Annual Review of Astronomy & Astrophysics,2003,24(1):127-170.

[4] Donoho D L.Compressed sensing[J]. IEEE Transactions on Information Theory,2006,52(4):1289-1306.

[5] Dimakis A G,Smarandache R,Vontobel P O.LDPC codes for compressed sensing[J].Mathematics,2011,58(5):3093-3114.

[6] Wenger S,Magnor M,Pihlstr?m Y,et al.Sparse RI: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] Bajwa W U,Sayeed A,Nowak R.Compressed sensing of wireless channels in time,frequency,and space[C]∥Proc of 42nd Asilomar Conf Signals,System and Computers,Pacific Grove,CA,2008:2048-2052.

[8] 陳建明,張 彬.低功耗無線傳感器能量供應裝置的探索[J].傳感器與微系統,2010,29(11):48-50.

[9] Wilson T L,Rohlfs K,Hüttemeister S.Tools of radio astronomy[M].Berlin:Springer,2009.

[10] 金 堅,谷源濤,梅順良.壓縮采樣技術及其應用[J].電子與信息學報,2010,32(2):470-475.

[11] 張春花,劉方愛,申志遠,等.一種新的異構無線傳感器網絡分簇算法[J].傳感器與微系統,2013,32(6):143-146.

[12] Wiaux Y,Puy G,Boursier Y,et al.Spread spectrum for imaging techniques in radio interferometry[J].Monthly Notices of the Royal Astronomical Society,2009,400(2):1029-1038.

[13] Burg J P.Maximum entropy speetral analysis[C]∥Proc of the 37th Meeting of the Society of Exploration Geophysists,1967.

[14] Burg J P.A new analysis technique for time series data[C]∥Proc of the NATO Advanced Study Institute on Signal Proeessing with Emphasis on Underwater Aeousties,1965:12-23.

[15] Frieden B R.Rcstoring with maximum like lihood and maximum entropy[J].Journal of Optical Society of America,1972,62(3):511-518.

[16] Gull S F,Skilling J.Maximum entropy methods in image proces-sing[C]∥IEE Proc of Special Issue Commun Randar Signal Process,1984:842-850.