劉曉曼，叢佳，朱永貴

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

1 引言

壓縮感知理論指出：如果信號(hào)在某個(gè)變換域是稀疏的或可壓縮的，就可以利用一個(gè)與變換基不相關(guān)的觀測(cè)矩陣將變換所得的高維信號(hào)投影到一個(gè)低維空間上，根據(jù)這些少量的觀測(cè)值，通過求解凸優(yōu)化問題就可以實(shí)現(xiàn)信號(hào)的精確重構(gòu)。此理論有效地減少了采集系統(tǒng)的復(fù)雜度，降低了系統(tǒng)的成本，加快了圖像采集速度。壓縮感知成像在醫(yī)療診斷方面起到了很重要的作用，特別是核磁共振成像的應(yīng)用。核磁共振成像技術(shù)(Nuclear Magnetic Resonance Imaging，簡(jiǎn)稱NMRI)，又稱磁共振成像(Magnetic Resonance Imaging，簡(jiǎn)稱MRI)，被廣泛應(yīng)用于醫(yī)學(xué)診斷中，極大地推動(dòng)了醫(yī)學(xué)、神經(jīng)生理學(xué)和認(rèn)知神經(jīng)科學(xué)的迅速發(fā)展。由大量資料調(diào)查表明，減少掃描時(shí)間即減少核磁共振成像所需要的時(shí)間是極其重要的，這就意味著在保證圖像質(zhì)量的情況下，要盡可能少地從頻域中收集觀測(cè)數(shù)據(jù)。然而，這看上去直接違背了長(zhǎng)期以來建立的奈奎斯特標(biāo)準(zhǔn)，即所獲得的觀測(cè)數(shù)據(jù)的數(shù)量必須至少與恢復(fù)圖像所需要的信息數(shù)量相匹配。這就意味著完美的圖像重構(gòu)將是不可能的，但壓縮傳感允許我們?cè)跊]有相關(guān)噪聲的情況下重構(gòu)圖像。

本文主要將不動(dòng)點(diǎn)連續(xù)性(Fixed Point Continuation，簡(jiǎn)稱FPC)[1][2]方法應(yīng)用到壓縮感知的核磁共振圖像中，通過不動(dòng)點(diǎn)定理高效的重構(gòu)出圖像，并通過數(shù)值實(shí)驗(yàn)證明，核磁共振圖像可以從全部數(shù)據(jù)的40%-50%抽樣中幾乎精確重構(gòu)。

2 預(yù)備知識(shí)

2.1 壓縮感知理論

設(shè)向量u代表一個(gè)圖像，一個(gè)壓縮算法，例如，通過找到某一變換φ(例如傅里葉或小波基)來壓縮這個(gè)圖像，使得φu=x是(近似)稀疏的。為了恢復(fù)u，我們使用相同的變換，并且設(shè)u=φ-1x。

壓縮感知的全部過程包括三步：編碼、傳感和解碼。在第一步中，通過一個(gè)線性變換R，u被編碼成一個(gè)維數(shù)較小的向量b=Ru，其中b的維數(shù)為m，x的維數(shù)為n，m

所以在壓縮感知的過程中，三大核心問題是運(yùn)用壓縮感知理論的關(guān)鍵問題。信號(hào)的可稀疏表示是壓縮感知的先驗(yàn)條件[3]。而在已知信號(hào)是可壓縮的前提下，壓縮感知過程可分為兩步：

(1)設(shè)計(jì)一個(gè)與變換基不相關(guān)的m×n(m?n)維測(cè)量矩陣對(duì)信號(hào)進(jìn)行觀測(cè)得到m維的測(cè)量向量。

(2)由m維的測(cè)量向量重構(gòu)信號(hào)。

2.2 信號(hào)的稀疏表示

稀疏的數(shù)學(xué)定義：信號(hào)x∈Rn在正交基ψ下的變換系數(shù)向量為?=ψTX，假如對(duì)于00，這些系數(shù)滿足：

(1)

則說明系數(shù)向量?在變換域ψ下是稀疏的。

廣義地，如果?中非零元素的個(gè)數(shù)為K，則信號(hào)X稱為K-稀疏的。

2.3 測(cè)量矩陣的構(gòu)造

測(cè)量矩陣，又稱傳感矩陣，從壓縮傳感整個(gè)過程可以看出，傳感矩陣起著相當(dāng)重要的作用。它的選擇合適與否，既關(guān)系到能否達(dá)到壓縮的目的，又直接關(guān)系到信號(hào)能否被精確重構(gòu)，此研究工作已有初步成果[4]。就矩陣選取而言，目前常用的傳感矩陣是隨機(jī)矩陣，如高斯矩陣、貝努力矩陣、傅里葉隨機(jī)測(cè)量矩陣、非相關(guān)測(cè)量矩陣等，它們已經(jīng)被驗(yàn)證滿足UUP特性[4]。

至于傳感矩陣的構(gòu)造，則首先需要研究傳感矩陣與重構(gòu)算法的關(guān)系。在重構(gòu)算法中，常用的匹配追蹤算法[5]、正交匹配追蹤算法等的重構(gòu)質(zhì)量與傳感矩陣有著密切的聯(lián)系。所以，如果能構(gòu)造出性質(zhì)好的矩陣，將大大簡(jiǎn)化重構(gòu)算法的步驟，并提高信號(hào)的重構(gòu)質(zhì)量。

2.4 重構(gòu)算法的設(shè)計(jì)

從本質(zhì)上講，壓縮感知理論中的信號(hào)重構(gòu)問題就是尋找欠定方程組的最稀疏解的問題。設(shè)u為傳統(tǒng)采樣得到的信號(hào)，長(zhǎng)度為n，而通過壓縮感知?jiǎng)t直接得到b，長(zhǎng)度為m(m

如果未知的其它信息，則上述逆問題很難求得唯一解。但若信號(hào)u是K稀疏的，設(shè)m≥K·log(n)，則上述方程組可以求解，通過下式求最稀疏的向量，從而獲得u：

(2)

其中‖·‖0表示u中非零元素的個(gè)數(shù)。然而，盡管從理論上來講，這種方法可以實(shí)現(xiàn)，但是在實(shí)際中不可行，因?yàn)檫@是一個(gè)組合問題(NP難問題)，計(jì)算量非常大。因此我們需要尋求合適的范數(shù)來近似求解上述問題。采用l1范數(shù)最小化求得的解，為稀疏向量，非常接近最小化l0范數(shù)所得的真實(shí)解。因此，我們可以選擇l1范數(shù)最小化來近似求解上述優(yōu)化問題，從而可將一個(gè)非凸優(yōu)化問題轉(zhuǎn)化為凸優(yōu)化問題。即

(3)

我們需要將信號(hào)變換到變換域考慮。設(shè)φ是壓縮基(如小波基)或緊框架，滿足規(guī)范正交，即φ*φ=I，作變換φu=x，顯然x是稀疏的。于是在這種情況下的求解方法為：

(4)

其中A=Rφ-1。

以上考慮的都是等式約束，然而實(shí)際中，測(cè)量過程可能會(huì)引入噪聲，這時(shí)約束條件(4)式中的Ax=b必須被放松，從而引出問題

(5)

或者問題

(6)

3 主要算法

3.1 最優(yōu)性和最優(yōu)解集

因?yàn)閒(x)為凸函數(shù)，定義X*為(6)的最優(yōu)解集。通過凸分析理論[6]可知(6)的最優(yōu)化條件為

x*∈X*?0∈SGN(x*)+μg(x*)

(7)

0是Rn中的零向量，g(x)=▽f(x)，SGN(x)為多重符號(hào)函數(shù)(集值函數(shù))

我們的算法基于算子分裂法。由凸分析理論[6]可以知道，求凸函數(shù)φ(x)的最小值問題相當(dāng)于找到次微分?φ(x)的零值。大多數(shù)情況下φ可以分裂成兩個(gè)凸函數(shù)φ=φ1+φ2。這意味著T可以分裂成兩個(gè)最大單調(diào)算子的和，即：T=T1+T2。對(duì)于τ>0，如果T2是單值的，I+τT1是可逆的，則可得到：

0∈T(x)?0∈(x+τT1(x))-(x-τT2(x))

?(I-τT2)x∈(I+τT1)x

?x=(I+τT1)-1(I-τT2)x.

(8)

式(8)得到了找到T的零點(diǎn)的向前向后分裂算法：

xk+1=(I+τT1)-1(I-τT2)xk

(9)

對(duì)于最小值問題(6)，對(duì)應(yīng)以上理論的具體函數(shù)表達(dá)式為

φ1(x)=‖x‖1，φ2(x)=μf(x)，

T1=?‖·‖1，T2=μ▽f

且對(duì)于(I+τT1)-1，通過[17]可知道，(I+τT1)-1是分段的收縮算子，即軟閾值。

3.2 不動(dòng)點(diǎn)算法

命題X*是(6)的最優(yōu)解集。對(duì)于任何標(biāo)量τ>0，x*∈X*當(dāng)且僅當(dāng)

(10)

解方程(10)，從而解決(6)的最小值問題，很自然考慮到不動(dòng)點(diǎn)的迭代

xk+1=sv·h(xk)v=τ/μ

(11)

雖然我們通過完全不同的逼近得到不動(dòng)點(diǎn)迭代上式，但后來我們發(fā)現(xiàn)這只是向前向后分裂法，當(dāng)T1(x)=?‖x‖1/μ，T2(x)=g(x)通過簡(jiǎn)單的計(jì)算可以得到：

sv=(I+τT1)-1，h=I-τT2。

(12)

其中κ(HEE)是HEE的條件數(shù)，滿足當(dāng)|E|很小時(shí)，可以使HEE小于H。這就意味著稀疏性可以加快收斂速度。

對(duì)于本文算法的觀測(cè)值b定義為

b=Axs+

(13)

3.3 增強(qiáng)的不動(dòng)點(diǎn)算法

為了收斂平穩(wěn)，在BB方法下直接自動(dòng)默認(rèn)為非單調(diào)線性搜索，即選擇一個(gè)α使迭代

xk+1=xk+a(sv°h(xk)-xk)

這種線性搜索結(jié)合Armijo型步長(zhǎng)條件，即得本文算法(FPC和BB steps方法和非單調(diào)線性搜索的結(jié)合算法)。

4 數(shù)值實(shí)驗(yàn)

在我們的實(shí)驗(yàn)中，由R估量得出的傅里葉系數(shù)并不是在隨機(jī)的頻率中選取的。我們是通過以下方式來選擇的：在k-space中，我們采取的采樣方式是沿著一定數(shù)量的從中心散開、呈輻射狀的直線來選取，即半徑抽樣。例如圖1、圖2分別顯示了在一個(gè)k-space中的22×2、22×5條輻射狀直線(是對(duì)210×210的大腦原始圖像進(jìn)行抽樣時(shí)的采樣情況)，即這些分別是44views、110views抽樣頻域圖。

圖1 采樣率20.95%

圖2 采樣率52.38%

所有實(shí)驗(yàn)都是在MATLAB 7.8.0版本上運(yùn)行的。執(zhí)行運(yùn)算的筆記本電腦是AMD turion X2 Dual-Core Mobile RM-72 處理器，3.00GB的內(nèi)存。在MATLAB中，基于FPC編碼[1]之上，并對(duì)FPC方法進(jìn)行加強(qiáng)，利用BB步和線性搜索方法，將其應(yīng)用到了真實(shí)的二維核磁共振圖像上。

我們?cè)谒膹埐煌姆叫味S核磁共振圖像上檢測(cè)我們的編碼。這四張圖像分別是：210×210的大腦原始圖像、220×220的腎臟動(dòng)脈MR原始圖像。在所有的檢測(cè)問題中，采用的噪聲是均值為0方差為0.01的高斯白噪聲。對(duì)這些圖像分別進(jìn)行44views、110views頻域抽樣，重構(gòu)的圖像如圖3、圖4所示。每組圖中圖(a)為原始圖像圖(b)、(c)為恢復(fù)后圖像。

在數(shù)值試驗(yàn)中，相對(duì)誤差(Relative Error)和信噪比(Signal to Noise Rations，簡(jiǎn)稱SNR)用來評(píng)估重構(gòu)圖像的質(zhì)量。相對(duì)誤差和信噪比的定義如下：

(14)

(15)

(a) (b) (c)圖3 (a)是原始大腦圖像；(b)、(c)是恢復(fù)后的圖像，它們的采樣率分別是20.95%、52.38%

(a) (b) (c)圖4 (a)是原始大腦圖像；(b)、(c)是恢復(fù)后的圖像，它們的采樣率分別是20.95%、52.38%

MRI圖像Views/采樣率Rerr/相對(duì)誤差SNR/信噪比(dB)迭代時(shí)間(秒)210×210的大腦MR原始圖像22×2/20.95%0.187314.71812.792422×3/31.43%0.165615.62002.979622×5/52.38%0.158715.98982.9484220×220的腎臟動(dòng)脈MR原始圖像22×2/20%0.157816.03982.839222×3/30%0.147216.64202.839222×5/50%0.105319.54893.2916220×220的胸部MR原始圖像22×2/20%0.114018.86072.823622×3/30%0.100020.00313.073222×5/50%0.096220.33782.9952256×256的心臟原始圖像22×2/17.19%0.257611.78144.258822×3/25.78%0.233912.61964.243222×5/42.97%0.204013.80624.2588

圖5 相對(duì)誤差與采樣率之間的關(guān)系

圖6 信噪比與采樣率之間的關(guān)系

圖7 迭代時(shí)間與采樣率之間的關(guān)系

對(duì)于同一個(gè)圖像，同樣的采樣率下，不同的迭代次數(shù)下，信噪比和相對(duì)誤差的變換。以210×210的大腦原始圖像在66views頻域抽樣為例，結(jié)果如表2：

表2

圖8 相對(duì)誤差與迭代次數(shù)之間的關(guān)系

圖9 信噪比與迭代次數(shù)之間的關(guān)系

圖10 迭代時(shí)間與迭代次數(shù)之間的關(guān)系

5 結(jié)論

在本篇論文中，基于壓縮感知理論，我們使用不動(dòng)點(diǎn)定理(FPC)的最小化模型，通過較少數(shù)量的觀測(cè)數(shù)據(jù)來恢復(fù)核磁共振圖像，使不動(dòng)點(diǎn)定理很好的運(yùn)用到了核磁共振圖像的重構(gòu)中，得到了一個(gè)有效的算法。在真實(shí)的幾幅正方形的核磁共振圖像上進(jìn)行的數(shù)值實(shí)驗(yàn)表明，這種算法可以在相對(duì)采樣率較小的情況下，用不到一分鐘的時(shí)間來恢復(fù)忠實(shí)于原圖的正方形圖像。通過對(duì)相對(duì)誤差、信噪比、迭代時(shí)間和迭代次數(shù)的對(duì)比，我們知道，在迭代15次左右基本達(dá)到效果。無需達(dá)到最大迭代次數(shù)。在采樣率在40%-50%時(shí)，信噪比和相對(duì)誤差成平穩(wěn)趨勢(shì)，達(dá)到精確重構(gòu)。而通過利用最優(yōu)化的技巧，例如線性搜索等，可使算法的速度加快。然而，本文并沒有對(duì)不動(dòng)點(diǎn)定理的收斂速度進(jìn)行系統(tǒng)的分析，對(duì)這種方法的理論分析也缺乏詳細(xì)的研究，這些問題都有待以后再解決。

[1] T H Elaine，Wotao Yin，Yin Zhang. A Fixed-Point Continuation Method for l1-Regularized Minimization with Applications to Compressed Sensing[R]. CAAM Technical Report TR07-07，2007.

[2] T H Elaine，Wotao Yin，Yin Zhang. Fixed-point Continuation Applied To Compressed Sensing：Implementation and Numerical Experiments[J]. Journal of Computational Mathematics，2010，28(2):170-194.

[3] 張雄偉，黃建軍，朱 濤.信息處理領(lǐng)域的創(chuàng)新理論——壓縮感知[J].軍事通信技術(shù)，2011，32(4)：83-97.

[4] 趙瑞珍.壓縮傳感與稀疏重構(gòu)的理論及應(yīng)用[DB/OL].http：//www.paper.edu.cn，中國科技論文在線.

[5] Neff R，ZAkhor A.Very Low bit rate video coding based on matching pursuits[J].IEEE Transactions on Circuits and Systems for Video Technology，1997，7(1)：158 -171.

[6] Rockafellar R. Convex Analysis[M].Princeton：Princeton University Press，1970.

[7] Chambolle A，DeVore R A，N-Y Lee，et al. Nonlinear wavelet image processing：variational problems，compression，and noise removal through wavelet shrinkage[J]. IEEE Transactions on Image Processing，1998，7(3)：319-335.

[8] Barzilai J， Borwein J. Two point step size gradient methods[J].IMA J Numer Anal，1988，8(1)：141-148.