不動點方法及其在壓縮感知核磁共振成像中的應用

劉曉曼，叢佳，朱永貴

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

1 引言

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

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

2 預備知識

2.1 壓縮感知理論

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

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

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

(1)設計一個與變換基不相關的m×n(m?n)維測量矩陣對信號進行觀測得到m維的測量向量。

(2)由m維的測量向量重構信號。

2.2 信號的稀疏表示

稀疏的數學定義：信號x∈Rn在正交基ψ下的變換系數向量為?=ψTX，假如對于00，這些系數滿足：

(1)

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

廣義地，如果?中非零元素的個數為K，則信號X稱為K-稀疏的。

2.3 測量矩陣的構造

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

至于傳感矩陣的構造，則首先需要研究傳感矩陣與重構算法的關系。在重構算法中，常用的匹配追蹤算法[5]、正交匹配追蹤算法等的重構質量與傳感矩陣有著密切的聯系。所以，如果能構造出性質好的矩陣，將大大簡化重構算法的步驟，并提高信號的重構質量。

2.4 重構算法的設計

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

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

(2)

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

(3)

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

(4)

其中A=Rφ-1。

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

(5)

或者問題

(6)

3 主要算法

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

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

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

(7)

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

我們的算法基于算子分裂法。由凸分析理論[6]可以知道，求凸函數φ(x)的最小值問題相當于找到次微分?φ(x)的零值。大多數情況下φ可以分裂成兩個凸函數φ=φ1+φ2。這意味著T可以分裂成兩個最大單調算子的和，即：T=T1+T2。對于τ>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的零點的向前向后分裂算法：

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

(9)

對于最小值問題(6)，對應以上理論的具體函數表達式為

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

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

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

3.2 不動點算法

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

(10)

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

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

(11)

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

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

(12)

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

對于本文算法的觀測值b定義為

b=Axs+

(13)

3.3 增強的不動點算法

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

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

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

4 數值實驗

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

圖1 采樣率20.95%

圖2 采樣率52.38%

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

我們在四張不同的方形二維核磁共振圖像上檢測我們的編碼。這四張圖像分別是：210×210的大腦原始圖像、220×220的腎臟動脈MR原始圖像。在所有的檢測問題中，采用的噪聲是均值為0方差為0.01的高斯白噪聲。對這些圖像分別進行44views、110views頻域抽樣，重構的圖像如圖3、圖4所示。每組圖中圖(a)為原始圖像圖(b)、(c)為恢復后圖像。

在數值試驗中，相對誤差(Relative Error)和信噪比(Signal to Noise Rations，簡稱SNR)用來評估重構圖像的質量。相對誤差和信噪比的定義如下：

(14)

(15)

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

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

MRI圖像Views/采樣率Rerr/相對誤差SNR/信噪比(dB)迭代時間(秒)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的腎臟動脈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 相對誤差與采樣率之間的關系

圖6 信噪比與采樣率之間的關系

圖7 迭代時間與采樣率之間的關系

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

表2

圖8 相對誤差與迭代次數之間的關系

圖9 信噪比與迭代次數之間的關系

圖10 迭代時間與迭代次數之間的關系

5 結論

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

[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] 張雄偉，黃建軍，朱 濤.信息處理領域的創(chuàng)新理論——壓縮感知[J].軍事通信技術，2011，32(4)：83-97.

[4] 趙瑞珍.壓縮傳感與稀疏重構的理論及應用[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.