朱永貴，劉 平，叢 佳

(中國傳媒大學(xué)理學(xué)院，北京100024)

Candes等［1］和Donoho［2］于2006年提出了壓縮感知(Compressed Sensing)理論。根據(jù)壓縮感知理論，如果稀疏信號(hào)或圖像的測量數(shù)據(jù)滿足不相關(guān)性［3］，那么稀疏信號(hào)或圖像可以從非常有限的抽樣數(shù)據(jù)中重構(gòu)出來，而這些測量數(shù)據(jù)被稱作原始稀疏信號(hào)或圖像的壓縮。對(duì)于二維圖像可以按一列一列地排列形成一個(gè)向量信號(hào)。用u表示由圖像形成的長度為N的向量信號(hào)。如果A表示一個(gè)M×N(M≤N)的觀測矩陣，那么在壓縮感知中b=Au是信號(hào)u的觀測數(shù)據(jù)即壓縮。原始信號(hào)u的重構(gòu)可以通過求解最優(yōu)化問題(1)實(shí)現(xiàn)。式中:代表‖u‖0中不等于零的分量個(gè)數(shù)。求解問題(1)是一個(gè)NP難問題［4］，因此該問題的求解依賴于計(jì)算機(jī)進(jìn)行數(shù)值計(jì)算難以實(shí)現(xiàn)。Trzasko等［5］于2009年提出了L0擬范數(shù)的同倫逼近壓縮感知重構(gòu)方法，這一方法能夠從抽樣率非常低的頻域數(shù)據(jù)中精確重構(gòu)出MR原始圖像。由于該方法涉及連續(xù)性求解過程，其計(jì)算時(shí)間很長，并且該方法目前僅能用于特殊類型的MR圖像重構(gòu)上，而不能用于其它類型如自然圖像、視頻圖像上，因此該方法不具有普遍性。盡管L0問題是NP難問題，但Candes等［6］證明了在一定條件下，問題(1)可以轉(zhuǎn)化為如下L1范數(shù)最小化問題

L1最小化問題是一個(gè)凸最優(yōu)化問題，利用最優(yōu)化方法可以求解。目前求解最優(yōu)化問題(2)的主要流行方法有內(nèi)點(diǎn)法［7］、不動(dòng)點(diǎn)連續(xù)性方法FPC［8－9］、Bregman迭代方法［10］、分裂 Bregman迭代法［11］、偏微分方程方法［12］和共軛梯度下降法［13］。最近我們也提出了求解問題(2)的凸分裂算法［14］。

當(dāng)0＜p＜1時(shí)，用于壓縮感知求解的最優(yōu)化問題變?yōu)槿缦路峭棺顑?yōu)化問題

同L0最小化問題一樣，問題(3)是非凸最優(yōu)化問題，但從計(jì)算角度看Lp(0＜p＜1)問題比L0問題容易處理。與問題(2)相比，問題(3)無論從理論分析還是數(shù)值計(jì)算，都非常困難。目前求解Lp(0＜p＜1)非凸問題已經(jīng)成為壓縮感知研究領(lǐng)域的一個(gè)熱點(diǎn)，引起了眾多科學(xué)家的重視。2007年Jung等［15］和Ye等［16］使用FOCUSS方法求解了最優(yōu)化問題(3)，取得了比較好的壓縮感知重構(gòu)結(jié)果。同年美國洛斯阿拉莫斯國家實(shí)驗(yàn)室的Chartrand教授［17］建立了求解Lp(0＜p＜1)問題的Fourier基最小化方法，這一方法已成為壓縮感知研究的主要方法。這一重要成果促進(jìn)了許多專家研究壓縮感知中求解Lp問題新的最優(yōu)化方法。例如Sidky等［18］使用Lp(0＜p＜1)方法利用非常少的 Views頻域抽樣精確重構(gòu)出原始圖像。Chartrand等［19］給出了壓縮感知的重加權(quán)迭代法。Chartrand［20］于2009年將其非凸最優(yōu)化的Fourier基算法推廣到核磁共振成像(Magnetic Resonance Imaging，MRI)領(lǐng)域中，得到了MR圖像重建非凸壓縮感知的快速算法。

本文是利用半二次罰函數(shù)法的算子分裂思想給出了求解壓縮感知Lp(0＜p＜1)非凸最優(yōu)化問題的分裂算法。其基本思想是將原問題分裂成兩個(gè)子問題:X子問題和Y子問題。然后對(duì)形成的X子問題通過光滑函數(shù)求導(dǎo)的方法給出X問題的閉形式解，對(duì)于Y子問題通過不動(dòng)點(diǎn)方法迭代求解。最后通過這兩個(gè)子問題的交替求解形成原問題的交替求解數(shù)值方法。我們從理論上分析和證明了Lp(0＜p＜1)非凸問題的分裂算法的收斂性。通過對(duì)X子問題和Y子問題的交替求解，能夠在誤差容許的范圍內(nèi)求得圖像壓縮感知Lp問題的精確重構(gòu)解。由于該方法的快速性使它不僅能夠進(jìn)行一維稀疏信號(hào)和二維稀疏圖像的精確重構(gòu)，而且還可應(yīng)用于核磁共振(MR)圖像的快速精確重構(gòu)。

1 非凸Lp最優(yōu)化問題的求解方法

在不考慮噪聲的情況下非凸Lp最優(yōu)化問題即為最優(yōu)化問題(3)。根據(jù)最優(yōu)化理論和凸分析知識(shí)，可以將問題(3)化為如下無約束形式

式中:‖·‖是指Euclid范數(shù)，x∈RN。

在含有噪聲污染的情況下，非凸Lp最優(yōu)化問題為如下形式

式中:n為高斯白噪聲。問題(5)的無約束形式如下

在MR圖像壓縮感知中，A=PFΦ－1，P∈RM×N，F(xiàn)是離散Fourier變換矩陣，Φ表示正交基如小波變換基、余弦變換或稀疏字典等。在數(shù)值模擬實(shí)驗(yàn)部分，將利用核磁共振圖像進(jìn)行壓縮感知重構(gòu)實(shí)驗(yàn)，并采用正交小波基作為稀疏基。

我們已在文獻(xiàn)［14］中利用半二次罰函數(shù)法，給出了壓縮感知凸最優(yōu)化問題的分裂算法。本文仍采用半二次罰函數(shù)法的思想將問題(6)變?yōu)?/p>

當(dāng)β→∞時(shí)，模型(7)等價(jià)于模型(6)，所以對(duì)于足夠大的β可以用式(7)的解去逼近式(6)的解。

還可以從設(shè)計(jì)階段來合理控制建筑結(jié)構(gòu)設(shè)計(jì)對(duì)造價(jià)成本的影響。設(shè)計(jì)階段需要充分關(guān)注建筑物的空間高度、使用功能以及內(nèi)部平面設(shè)計(jì)等方面，通過對(duì)這些方面的科學(xué)合理設(shè)計(jì)可以促進(jìn)結(jié)構(gòu)設(shè)計(jì)工作的有效開展。例如，在高層建筑中，隨著建筑層數(shù)的增多，水平荷載也成為建筑結(jié)構(gòu)設(shè)計(jì)的控制因素。

下面給出x－子問題和y－子問題的求解方法。對(duì)于x－子問題，由于目標(biāo)函數(shù)是光滑函數(shù)，所以通過求導(dǎo)的方法可以得到其最優(yōu)解

而y－子問題是一個(gè)非凸問題，可以利用不動(dòng)點(diǎn)方法求解，現(xiàn)在介紹具體的求解過程。

根據(jù)Lp范數(shù)的定義可以將y－子問題化為如下形式

因?yàn)閥的各分量y1，y2，…，yN與x的各分量x1，x2…，xN是相互獨(dú)立的，因此子問題(9)又可以變?yōu)榉至啃问?/p>

從而可以求解yj的不動(dòng)點(diǎn)迭代公式

下面證明和分析迭代算子是收縮的。

定理1 對(duì)于0＜p＜1，實(shí)數(shù)域R中的任何兩個(gè)點(diǎn)y(n+1)和y(n)，對(duì)于足夠大的β有

證明 不適一般性，不妨設(shè)y(n+1)和y(n)都大于0，則有

由于0＜p＜1，β足夠大，所以

證畢。

2 數(shù)值模擬

在數(shù)值模擬實(shí)驗(yàn)中，采用信噪比(SNR)和相對(duì)誤差(Relative Error)評(píng)價(jià)壓縮感知圖像重構(gòu)的質(zhì)量。信噪比和相對(duì)誤差的定義分別如下

式中:x和x0各自代表計(jì)算重構(gòu)的圖像和原始圖像。

在下面的數(shù)值模擬實(shí)驗(yàn)中，所采用的噪聲是均值為0，方差為0.01的Gauss白噪聲，稀疏基采用正交Harr小波。并且用非凸分裂算法進(jìn)行數(shù)值模擬實(shí)驗(yàn)過程中，假設(shè)圖像的初始值x(0)=0，令μ=1000，p=1/3，取迭代終止準(zhǔn)則為:

圖1(a)、(b)分別給出220×220的心臟MR原始圖像和對(duì)其進(jìn)行77 Views頻域抽樣圖。對(duì)圖1采用非凸分裂算法進(jìn)行重構(gòu)實(shí)驗(yàn)，其結(jié)果如圖2所示。

圖1 原始圖像和頻域空間中的77 ViewsFig.1 Original images and 77 Views in frequency space

圖2 非凸分裂算法重構(gòu)的圖像Fig.2 Reconstructions by non-convex splittingmethod

計(jì)算得到非凸分裂算法重構(gòu)的大腦MR圖像的信噪比和相對(duì)誤差分別為 17.802 1 dB和0.134 2。

現(xiàn)在對(duì)大腦MR圖像進(jìn)行121Views頻域抽樣(見圖3(a))，并采用凸分裂算法進(jìn)行重構(gòu)實(shí)驗(yàn)，重構(gòu)結(jié)果如圖3(b)所示。凸分裂算法重構(gòu)的大腦MR圖像的信噪比和相對(duì)誤差分別為15.361 2 dB和0.187 1。從重構(gòu)圖像的結(jié)果可以看出非凸分裂算法重構(gòu)的圖像質(zhì)量明顯好于凸分裂算法重構(gòu)的圖像，但非凸分裂算法的抽樣率僅是凸分裂算法的67.77%。

圖3 頻域空間中的121Views和凸分裂算法重構(gòu)的圖像Fig.3 121Views and Rec.by convex splitting method

下面給出的身體MR圖像(見圖4(a))，對(duì)其進(jìn)行抽樣率為23.56%的笛卡爾線形抽樣(見圖4(b))。

圖4 原始圖像和頻域空間中的笛卡爾線型抽樣Fig.4 O riginal images and Cartesian frequency sam p ling

現(xiàn)采用非凸分裂方法進(jìn)行重構(gòu)實(shí)驗(yàn)，重構(gòu)結(jié)果如圖5所示。其信噪比為15.532 9 dB，相對(duì)誤差為0.124 1。

圖5 非凸分裂算法重構(gòu)的圖像Fig.5 Reconstructions by non-convex splittingmethod

對(duì)身體MR圖像進(jìn)行抽樣率為40.87%的笛卡爾線形頻率抽樣(見圖6(a))，并采用凸分裂算法進(jìn)行重構(gòu)實(shí)驗(yàn)，重構(gòu)結(jié)果如圖6(b)。重構(gòu)圖像的信噪比為13.421 3 dB，相對(duì)誤差為0.213 5。

圖6 40.87%的笛卡爾抽樣和凸分裂算法重構(gòu)的圖像Fig.6 40.87%sam pling and Rec.by convex method

從重構(gòu)圖像的結(jié)果可以看出非凸分裂算法重構(gòu)的圖像質(zhì)量明顯好于凸分裂算法重構(gòu)的圖像，而且抽樣率僅是凸分裂算法的57.65%。

3 結(jié)束語

本文給出了求解壓縮感知中稀疏圖像重構(gòu)Lp問題的非凸分裂方法。通過數(shù)值模擬實(shí)驗(yàn)對(duì)非凸分裂方法與凸分裂方法進(jìn)行了比較，實(shí)驗(yàn)結(jié)果表明:

(1)非凸分裂方法重構(gòu)的圖像比凸分裂方法重構(gòu)的圖像具有更高的精度和更低的頻域抽樣率。

(2)通過主觀觀察，在圖像壓縮感知中非凸分裂方法較凸分裂方法具有更好的圖像重構(gòu)質(zhì)量。

本文給出的壓縮感知中稀疏圖像重構(gòu)的非凸分裂方法可以降低稀疏重構(gòu)的抽樣率，能夠用于核磁共振快速成像，在醫(yī)學(xué)成像領(lǐng)域中具有一定的應(yīng)用價(jià)值。

［1］Candes E J，Romberg J，Tao T.Robust uncertainty principles:Exact signal reconstruction from highly incompletee frequency information［J］.IEEE Trans Inform Theory，2006，52(2):489－509.

［2］Donoho D L.Compressed sensing［J］.IEEE Trans Inform Theory，2006，52(4):1289－1306.

［3］Candes E J，Romberg J.Sparsity and incoherence in compressive sampling［J］.Inverse Problems，2007，23 (3):969－985.

［4］Natarajan B K.Sparse approximation solutions to linear systems［J］.SIAM JComput，1995，24(2):227－234.

［5］Trzasko J，Manduca A，Borisch E.Highly undersampled magnetic resonance image reconstruction via homotopic L0－minimization［J］IEEE Transactions on Medical Imaging，2009，28(1):106－121.

［6］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.

［7］Kim S J，Koh K，Lustig M.An interior－pointmethod for large－scale l1－regularized least squares［J］.IEEE Trans.on Selected Topics in Signal Processing，2007，1 (4):606－617.

［8］Hale E，Yin W，Zhang Y.Fixed－point continuation for l1 minimization:Methodology and convergence［J］.SIAM Journal on Optimization，2008，19(3):1107－1130.

［9］Hale E，Yin W，Zhang Y.Fixed－point continuation applied to compressed sensing:Implementation and numerical experiments［J］.Journal of Computational Mathematics，2010，28(2):170－194.

［10］Yin W，Osher S，Goldfarb D.Bregman iterative algori thm for l1－minimization with applications to compressed sensing［J］.SIAM J.Image Sciences，2008，1(1):143－168.

［11］Goldstein T，Osher S.The split bregman method for L1－regularized problems［J］.SIAM J.Image Sciences，2009，2(2):323－343.

［12］He L，Chang TC，Osher S，et al.MR image reconstruction by using the iterative refinementmethod and nonlinear inverse scale space methods［R］.UCLA CAM Report，2006:6－35.

［13］Lustig M，Donoho D，Pauly J.Sparse MRI:The application of compressed sensing for rapid MR imaging［J］. Magnetic Resonance in Medicine，2007，58(6):1182－1195.

［14］朱永貴，楊曉蘭.稀疏MR圖像重構(gòu)的快速算法［J］.中國圖象圖形學(xué)報(bào)，2011，16(9):1736－1744.

Zhu Yong－gui，Yang Xiao－lan.Fast reconstructionmethod for sparse MR image［J］.Journal of Image and Graphics，2011，16(9):1736－1744.

［15］Jung H，Ye J，Kim E.Improved k－t blask and k－t sense using focuss［J］.Phys Med Biol，2007，52(11):3201－3226.

［16］Ye J C，Tak S，Han Y et al.Projection reconstruction MR imaging using FOCUSS［J］.Magnetic Resonance in Medicine，2007，57(4):764－775.

［17］Chartrand R.Exact reconstruction of sparse signals via nonconvexminimization［J］.IEEE signal Letters，2007，14(10):707－710.

［18］Sidky E Y，Chartrand R，Pan X.Image reconstruction from few views by non－convex optimization［R］. IEEE transaction Medical Imaging Conference Record，2007.

［19］Chartrand R，Yin W.Iteratively reweighted algorithms for compressive sensing［C］//33rd International Conference on Acoustics Speech，and Signal Processing，2008.

［20］Chartrand R.Fast algorithms nonconvex compressive sensing:MRI reconstruction from very few data［R］. IEEE International Symposium on Biomedical Imaging，2009.