黃敏，林穎，陳軍波

(中南民族大學(xué) 生物醫(yī)學(xué)工程學(xué)院，武漢 430074)

1 引 言

磁共振成像技術(shù)是一種無電離輻射，多層面多參數(shù)的重要醫(yī)學(xué)成像方式，已廣泛應(yīng)用于臨床和科研[1]。但MRI成像數(shù)據(jù)采集掃描過長，使該項檢查的適應(yīng)癥大為減少。因此，提高磁共振成像速度一直是磁共振成像領(lǐng)域要解決的核心目標(biāo)之一。

傳統(tǒng)的快速成像方法是通過提高梯度場性能來加快 k 空間數(shù)據(jù)采集速度。然而, 受制于硬件制約和生理因素,這類方法加速效果已經(jīng)接近極限。近年誕生的壓縮感知(compressed sensing,CS)新理論[2-3],提供了一種不需要遵循采樣定理的欠采樣新方式，可明顯縮短數(shù)據(jù)采集時間，并可從極少的采樣數(shù)據(jù)中恢復(fù)出原始信號。

MRI圖像在適當(dāng)?shù)淖儞Q域中(如小波變換域)是可壓縮的，且 MRI 采集的原始數(shù)據(jù)不是直接的圖像信息，而是k空間頻率域信息，因此，MRI可以采用壓縮感知的采集方式[4]。

要在MR掃描儀上實現(xiàn)壓縮感知MRI新技術(shù)，需要解決兩個關(guān)鍵問題：一是在掃描儀中實現(xiàn)壓縮感知序列，二是對欠采樣數(shù)據(jù)進(jìn)行壓縮感知重建。針對后一個問題，國內(nèi)外已經(jīng)進(jìn)行了大量的圖像重建方法的仿真。但在掃描儀上最終實現(xiàn)該技術(shù)的關(guān)鍵問題——壓縮感知序列的設(shè)計和實現(xiàn)，還很少涉及。序列的設(shè)計與實現(xiàn)也是目前國內(nèi)MR掃描儀自主研發(fā)的重點。

本研究主要就壓縮感知磁共振脈沖序列的設(shè)計、編程和實現(xiàn)方法進(jìn)行研究，設(shè)計不同的加速因子和序列，對掃描得到的數(shù)據(jù)進(jìn)行離線重建，用重建結(jié)果來驗證序列的正確性。

2 壓縮感知MRI序列原理

壓縮感知MRI與傳統(tǒng)MRI的信號采樣不同， 它只需測量少部分隨機(jī)線性組合信號，通過非線性壓縮感知方法從測量值中高準(zhǔn)確度地重建出圖像。將壓縮感知理論應(yīng)用到MRI序列中，可通過隨機(jī)欠采樣序列的實現(xiàn)，采集少量的K空間數(shù)據(jù)，提高成像速度，減少患者在掃描過程中的不適，更好地服務(wù)臨床。

壓縮感知MRI序列主要是要確定欠采樣的方式?，F(xiàn)有的全采樣方式主要有二維笛卡爾軌跡、三維笛卡爾軌跡、二維螺旋軌跡等。在這些基礎(chǔ)上進(jìn)行欠采可分為一維隨機(jī)欠采和二維隨機(jī)欠采。隨機(jī)欠采的相關(guān)研究有很多，如Lustig[4-6], Wang Haifeng[7]等人提出的1-D隨機(jī)采樣[4]，2-D隨機(jī)采樣[4,7]，螺旋隨機(jī)采樣[6]等。一維隨機(jī)欠采主要是在相位編碼方向隨機(jī)欠采，二維隨機(jī)欠采是在相位編碼方向和頻率編碼方向同時進(jìn)行隨機(jī)欠采。

設(shè)計和實現(xiàn)壓縮感知MRI序列可以有多種實現(xiàn)”基”?？梢栽谧孕夭?SE)序列基礎(chǔ)上實現(xiàn)，也可以在螺旋序列及其它序列基礎(chǔ)上實現(xiàn)。國內(nèi)MR掃描儀由于梯度系統(tǒng)及其它硬件配備很難實現(xiàn)螺旋等快速采集方式，從這種現(xiàn)狀考慮，我們基于自旋回波序列的基礎(chǔ)上進(jìn)行設(shè)計。

在實際儀器上，對于笛卡兒采集的二維成像：由于譜儀采集數(shù)據(jù)都是等時間間隔采集數(shù)據(jù)，要在頻率編碼上實現(xiàn)欠采樣必須重新調(diào)整采樣數(shù)據(jù)的例程；另外，由于TR的限制，頻率編碼方向上即便實行了欠采樣，對掃描時間也沒有實質(zhì)上的減少，只有相位編碼方向的欠采樣才能明顯減少采集時間。因此，本文進(jìn)行1D隨機(jī)欠采的設(shè)計及實現(xiàn)。

對于1D加速的壓縮感知MRI序列,與常規(guī)自旋回波的序列不同的是：相位編碼梯度幅值不再均勻變化，而是施加‘偽隨機(jī)’的變間隔梯度。壓縮感知MRI脈沖序列及回波信號采集見圖1。

圖1 壓縮感知MRI序列圖

若1D欠采樣加速因子為F，實際采集數(shù)據(jù)時間為全采集的1/F。一維相位編碼方向上隨機(jī)采樣使K空間的中心附近密集采集，周圍稀疏采集[5]。變密度采樣的概率密度函數(shù)為：

(1)

其中r的取值范圍為(0，1)，p為多項式的系數(shù)，r為采樣點到k空間中心的距離；

通過以上概率密度函數(shù)設(shè)計出欠采樣軌跡，根據(jù)設(shè)計的欠采樣模式即可得到相應(yīng)的相位編碼次序，從而在序列中實現(xiàn)。而CS-SE序列即將SE序列的相位編碼次序改變成設(shè)計的隨機(jī)相位編碼次序即可。

3 壓縮感知MRI序列設(shè)計

本設(shè)計是在深圳貝斯達(dá)醫(yī)療器械有限公司磁共振設(shè)備上完成的，由于各方面的局限性，只能實現(xiàn)笛卡爾軌跡采集，并且只能在相位編碼方向上欠采。

由CS 理論可知：對 K 空間數(shù)據(jù)在相位編碼方向上進(jìn)行隨機(jī)欠采，再采用非線性重建即能重構(gòu)出圖像。我們先通過MATLAB設(shè)計相位編碼隨機(jī)欠采模式，得到不同分辨率和不同加速因子條件下，K空間欠采樣的相位編碼序號，并將其做成表格存儲，然后將對應(yīng)采樣行寫入序列中。

CS-SE序列是利用MR Solutions公司的磁共振譜儀對應(yīng)的PPL (Pulse Program Language) 語言，來設(shè)計和編寫的，主要對射頻脈沖 RF，選層梯度 Gs，相位編碼梯度 Gp 和頻率編碼梯度 Gr 進(jìn)行設(shè)計，包括各脈沖的形狀，脈沖寬度，時序等的控制[8]。

3.1 射頻脈沖選擇

磁共振2D成像要選擇性激發(fā)具有一定厚度的層面，需要采用軟脈沖來實現(xiàn)選層。選層需要在頻域上為一理想窄帶的矩形形狀最佳，時域無限長sinc信號可滿足此要求。但實際使用時通常選擇一定數(shù)量的主旁瓣，以激發(fā)不同區(qū)域，常用3瓣或5瓣sinc信號。

MR Solutions公司的磁共振譜儀底層軟件中，已經(jīng)存在著多種固定類形的RF脈沖，如sinc型、gauss型及正切型的波形，并且各種不同波形的RF脈沖也根據(jù)帶寬大小的不同而進(jìn)行分類編號。設(shè)計RF脈沖需要調(diào)用這些脈沖類型。如果序列開發(fā)者覺得已有RF脈沖不滿足應(yīng)用，也可根據(jù)自己的需求在波形編輯器中編輯生成想要的波形并保存到RF波形文件RFstd.seq中再進(jìn)行調(diào)用即可。在程序開始處需要進(jìn)行RF脈沖波形調(diào)用的聲明，再用對應(yīng)的語句調(diào)出可能使用到的波形，如下所示：

#use RF1 "c:smisseqlibRFstd.seq" pf1

NEWSHAPE_MAC(1,pf1,"3lobe_sinc_3kHz",1332, 3000);NEWSHAPE_MAC(7, pf7,"gauss_240Hz" , 7500, 240);NEWSHAPE_MAC(8,pf8,"5sinc_7_5cos_6ms",6000, 1500)

3.2 梯度脈沖設(shè)計

MR儀器上使用的梯度脈沖是梯形狀的正脈沖或者負(fù)脈沖，正負(fù)脈沖的上升和下降速度相同。為了輸出不同的梯度，通常是改變梯度的強(qiáng)度及持續(xù)時間，各個方向上的梯度脈沖的強(qiáng)度大小需要定義，并存儲到對應(yīng)的梯度矩陣中，而持續(xù)時間也是在程序中定義的。梯度強(qiáng)度大小使用CREATE_MATRIX(mat, slice, phase, read)命令來定義，mat是定義的矩陣名稱，slice, phase, read是固定的三個方向上的梯度值，以相位編碼梯度矩陣pe_mat為例，定義如下：

CREATE_MATRIX(pe_mat, gs_on*gs_comp_s, gp_on*gp_var, gr_on*gr_comp_s)

delay( caldelay, us );

gs_on、gp_on、gr_on只能取0或者1，取0則對應(yīng)方向上這段時間內(nèi)的梯度值為0，取1則輸出對應(yīng)的梯度值，caldelay是程序中創(chuàng)建矩陣及計算矩陣所需要的等待時間。梯度方向上的梯度輸出同樣需要定義，用到MR3040_InitList()。譜儀底層軟件中已有定義好的正負(fù)梯度波形，可以直接調(diào)用POSPULSE_HOLD、NEGPULSE_HOLD來組成各個方向上的梯度，也可用 MR3040_OUTPUT來分別輸出各個梯度，具體代碼如下：

MR3040_SetListAddress(0);

slice_list_1 = MR3040_InitList();

POSPULSE_HOLD()

if (slice_rephase = 0 & flow_comp_on == 0) {NEGPULSE(tref, clock)}

if (slice_rephase == 0 & flow_comp_on == 1) {NEGPULSE((2*tref + tramp), clock)

POSPULSE_SEC((tsel90/2-tramp/2), clock)}

slice_list_2 = MR3040_InitList();

POSPULSE_HOLD()

if ( flow_comp_on) {NEGPULSE_HOLD()}

ref_list = MR3040_InitList();

if ( flow_comp_on == 0) {POSPULSE(tref, clock)}

else {NEGPULSE_SEC_HOLD()

POSPULSE_SEC(tsel90/2, clock)}

read_list = MR3040_InitList();

POSPULSE_HOLD()

phase_list = MR3040_InitList();

POSPULSE(tref, clock)

slice_list_1和slice_list_2分別定義了90°和180°射頻脈沖輸出時對應(yīng)的選層梯度的脈沖波形；read_list和ref_list 定義了頻率編碼通道輸出的梯度波形；phase_list定義了相位編碼梯度的輸出波形。

3.3 梯度脈沖設(shè)計

為了找到適合不同場強(qiáng)掃描儀的壓縮感知欠采樣加速因子F，我們實現(xiàn)了F分別為1.5、2、2.5、3、4、5、8等的CS-MRI序列，實際采樣矩陣大小為256×256/F，并取整數(shù)。實現(xiàn)過程分為以下兩步：(1)用MATLAB產(chǎn)生1D變密度隨機(jī)采樣模板的概率密度函數(shù)及對應(yīng)的采樣模式；(2)用PPL語言實現(xiàn)(1)中用MATLAB設(shè)計的采樣模式。

以欠采加速因子F=3為例，實際采樣矩陣大小則為256×256/3(即256×84)，取r=0.15，p=6及r=0.1，p=11得到兩種采樣模式分別見圖2(a)、(b)，圖中水平方向是頻率編碼，垂直方向是相位編碼，也就是相應(yīng)的k空間填充順序，白色表示采集，黑色不采集。

圖2 壓縮感知欠采樣的K空間填充順序

為了在盡可能節(jié)省時間的情況下得到較好的圖像結(jié)果，就要對不同部位進(jìn)行多種不同欠采樣加速因子選取的實驗。我們先在MATLAB中得到欠采樣加速因子對應(yīng)的采樣模式，再將采樣模式的相位編碼行號寫入到SMIS譜儀可讀取的pph文件中，記為gp.pph。在程序開始時對該文件進(jìn)行聲明，并在程序中進(jìn)行相應(yīng)的調(diào)用。在PPL程序中具體的實現(xiàn)如下：

#include "C:smisincludegp.pph"

phase_encode_loop:

current_view = current_view + 1;

if (no_views=256){gp_var=gp_var-gp_inc;}

if(no_views=84){gp_var=(-gp_inc*(gp_mul2[current_view]-128));}

if (no_views=103){gp_var=(-gp_inc*(gp_mul3[current_view]-128));}

if(no_views=128){gp_var=(-gp_inc*(gp_mul5[current_view]-128));}

if(no_views=171){gp_var=(-gp_inc*(gp_mul4[current_view]-128));}

if (current_view < no_views) goto phase_encode_loop;

3.4 時序設(shè)計

MRI序列中施加的各個脈沖的順序及各脈沖間的延遲時間的長短都會對最終的圖像產(chǎn)生一定的影響，為了保證數(shù)據(jù)的準(zhǔn)確性，需要對脈沖時序進(jìn)行嚴(yán)格要求，下面以讀出梯度的時序為例進(jìn)行說明。

starttimer();

ret = tramp * 10 ;

MR3040_SelectMatrix( aq_mat );

MR3040_Start(CHANNEL_R);

frequency_buffer(1);

reset_frequency();

waittimer(ret);

acquire( sample_period, no_samples );

MR3040_Continue(CHANNEL_R);

delay(tramp,us);

用starttimer來提示計時開始，ret為定義的計時時間的大小，tramp是梯度脈沖上升沿及下降沿的持續(xù)時間，ret、tramp都是以us為單位的，MR3040_SelectMatrix用來選擇即將輸出的梯度強(qiáng)度矩陣，aq_mat是當(dāng)前選擇的讀梯度強(qiáng)度的大小矩陣，CHANNEL_R 指的是讀梯度通道，MR3040_Start(CHANNEL_R)則是執(zhí)行開啟讀梯度通道命令，frequency_buffer頻率緩存器，可存儲多達(dá)到10個頻率，buffer 0指的是發(fā)射頻率，buffer 1指的是接收機(jī)頻率，當(dāng)前選擇存儲了接收機(jī)頻率，并需要用reset_frequency()來重置一下，才能將頻率切換過來，starttimer()與waittimer()是配套使用的，兩者之間的持續(xù)時間即是定義的ret，包括了代碼執(zhí)行時間，用acquire來提示進(jìn)行數(shù)據(jù)的采集，sample_period是兩個數(shù)據(jù)點之間的采樣間隔，no_samples 是采樣點數(shù)，在參數(shù)設(shè)置時設(shè)定的，MR3040_Continue用來結(jié)束梯度的輸出。

時序的精確性體現(xiàn)在程序開始時用MR3040_Clock(clock) 開啟定時器，計時開始時，配套運用starttimer和waittimer進(jìn)行時間的精確控制，及運用delay命令進(jìn)行延時，及程序代碼執(zhí)行時間也被考慮計算進(jìn)去。

程序編寫完成后進(jìn)行編譯，生成的參數(shù)界面反映了在脈沖序列程序編寫過程中設(shè)置的參數(shù)范圍，一旦設(shè)置超出范圍，就會進(jìn)行錯誤提示。同時界面上的各個參數(shù)的數(shù)值與實際數(shù)據(jù)之間有各自不同的對應(yīng)關(guān)系，有的是程序編寫者人為定義的，有的是譜儀軟件里固定的，若需要使用這些值，則需要進(jìn)行一個轉(zhuǎn)換。序列設(shè)計者可進(jìn)一步根據(jù)自己的需求修改參數(shù)，如頻率編碼數(shù)、相位編碼數(shù)、層數(shù)、層厚、層間距、平均次數(shù)等等，來獲得自己想要的序列。

4 結(jié)果與討論

我們在深圳貝斯達(dá)醫(yī)療器械有限公司的MR掃描儀(包括0.5T永磁型及1.5T超導(dǎo)型)上進(jìn)行了序列的編寫和調(diào)試，并進(jìn)行了CS-SE序列采集試驗。

對于0.5T MR掃描儀試驗數(shù)據(jù)，設(shè)計了欠采樣率為1/5、1/3、1/2.5、1/2、2/3、1共六種情況下的2種不同采樣模式，采集得到膝蓋和頭部的原始數(shù)據(jù)。成像參數(shù)為：(1)頭部數(shù)據(jù)：TR=520 ms，TE=20 ms，Number of Samples=256，SLICE_THICKNESS=5 mm，BW=25 KHz，SLICE_SEPARATION=1 mm，F(xiàn)OV=224 mm，NEX=2，TA=4 min26s； (2)膝蓋數(shù)據(jù)：TR=400 ms，TE=20 ms，Number of Samples=256，SLICE_THICKNESS=5 mm，BW=25 KHz， SLICE_SEPARATION=1 mm，F(xiàn)OV=224 mm，NEX=3，TA=5 min7s。CS-SE序列中的相位編碼步Number of Views，分別取52、84、103、128、171、256。

對全采集數(shù)據(jù)(256行)，直接采用快速傅立葉變換進(jìn)行重建。對各組壓縮感知欠采數(shù)據(jù)(小于256行)，按照MR掃描儀MRD文件的存儲格式進(jìn)行數(shù)據(jù)讀取，采用本實驗室在VC6.0環(huán)境下編寫的壓縮感知重建程序進(jìn)行離線重建[9]，該程序是基于M. Lustig的非線性共軛梯度下降法，采用迭代方式求取最優(yōu)值進(jìn)行圖像重建的。

1D隨機(jī)欠采樣采用了k-空間中心附近密集，周圍稀疏的變密度采集方式，而k空間中心半徑r及多項式次數(shù)p的選取對成像效果影響很大。對不同的欠采模式1和2下的頭部數(shù)據(jù)，采用相同的Db小波結(jié)合全變差分的壓縮感知方法進(jìn)行重建，當(dāng)?shù)螖?shù)取72次，重建出一幅256×256圖像需要時間為203 s，得到重建結(jié)果見圖3。從左至右欠采樣加速因子分別為 5、3、2.5、2、1.5、1。由于采樣模式2采用的中心半徑r=0.1比采樣模式1的r=0.15要小的多，低頻區(qū)域損失信息相對較多，加上儀器本身的噪聲，在同欠采樣率下，重建效果要差的多。由此可見：為了能達(dá)到最優(yōu)的結(jié)果，首先要考慮中心區(qū)域應(yīng)取k空間最大值的多少比例，以及各區(qū)域的采樣密度的分配，才能設(shè)計出合理的欠采樣模式。

圖3頭部CS-SE序列成像結(jié)果

(a)采樣模式1；(b)采樣模式2

Fig3ImagingresultofCS-SEsequenceonbrain

(a)Model 1; (b)Model 2

圖像稀疏變換基有很多，如有限差分、小波基、curvelet基、Harr基、Contourlets基等等。不同部位的成像數(shù)據(jù)的稀疏度是不同的，選用不同的稀疏變換基，對圖像重建速度和效果有很大的影響。我們對同一采樣模式下加速因子為3的頭部兩層的數(shù)據(jù)以及膝蓋數(shù)據(jù)，采用多種不同的稀疏變換后重建圖像，結(jié)果見圖4。從左至右稀疏方式分別為全變差分、5層Db小波、5層Symmlet小波、5層Battle小波、8層Db小波。迭代次數(shù)選用45，重建出256×256圖像所需要時間分別為：頭部13、94、92、118 s，膝蓋13、92、100、118 s。注：最右邊的是全采樣，采用快速傅立葉變換重建。由此可見：對于稀疏程度高的膝蓋數(shù)據(jù)，只需選擇簡單的差分運算便可得到很快且較好的結(jié)果；對于稀疏程度相對較低的頭部數(shù)據(jù)，采用5層Db4的稀疏變換方式較優(yōu)。

當(dāng)然，迭代次數(shù)同樣影響著圖像重建效果。迭代次數(shù)越多，重建結(jié)果也就越接近全采樣值。不同的稀疏變換下的迭代時間也各不相同，如：全變差分變換每增加一次迭代時間增加0.29 s，而5層的Db4變換一次迭代要2 s。所以，在不同稀疏變換的選擇下，若是重建效果差不多，則盡可能選擇重建時間短的變換；而且對于不同部位的數(shù)據(jù)，最佳的迭代次數(shù)也會不同。我們還在1.5T 超導(dǎo)MR掃描儀上進(jìn)行了序列調(diào)試和采集數(shù)據(jù)的試驗，成像參數(shù)同0.5T掃描儀。當(dāng)加速因子分別為5、4、3、2.5、2，1時(即采樣行數(shù)為52、64、84、103、128，256)，頭部兩個部位以及膝蓋的重建結(jié)果見圖5。根據(jù)上述掃描數(shù)據(jù)及重建結(jié)果，我們認(rèn)為：在貝斯達(dá)公司的掃描儀上運用CS-SE序列后，綜合考慮掃描時間及圖像效果兩個方面，頭部合適的欠采樣加速因子為2，膝蓋則為4。此時將采集并重建圖像的時間分別縮短到SE序列的1/2及1/4，掃描CS-SE序列獲得膝蓋橫斷位數(shù)據(jù)需要1 min 17 s，頭部橫斷位數(shù)據(jù)需要2 min 23 s，全采膝蓋需要5 min 7 s，全采頭部需要4 min 26 s，CS-SE大大節(jié)省了數(shù)據(jù)采集時間。不同部位，由于圖像本身的稀疏程度不同，膝蓋部位可采用比頭部高的欠采樣加速因子。

圖4不同稀疏變換的重建結(jié)果

(a) 頭部層面1 ；(b) 頭部層面2；(c) 膝蓋
Fig4Reconstructionresultfromdifferentsparsetransforms

(a)Slice 1 on brain; (b)Slice 2 on brain;(c)On knee

圖5 1.5T超導(dǎo)MR掃描儀CS-SE序列成像結(jié)果

由此可見：實際儀器采集CS數(shù)據(jù)可以實施的加速因子，與很多文獻(xiàn)中采用‘仿真’的理想數(shù)據(jù)得到的可欠采樣十倍甚至幾十倍的結(jié)論相差甚遠(yuǎn)。因此，我們必須考慮實際采集數(shù)據(jù)中不可預(yù)見的各種噪聲，不能以仿真數(shù)據(jù)作為推斷。壓縮感知成像技術(shù)要真正用于商業(yè)掃描儀上，還得進(jìn)行大量的試驗。

5 結(jié)論

目前，國內(nèi)MR掃描儀的核心部分(包括譜儀，脈沖序列，線圈等)的自主研發(fā)還處于薄弱環(huán)節(jié)，壓縮感知序列更沒有應(yīng)用于商用掃描儀上。本研究介紹的基于壓縮感知的磁共振序列，對序列的設(shè)計和編程，圖像的重建進(jìn)行了介紹，可為國內(nèi)MR開發(fā)領(lǐng)域在脈沖序列設(shè)計及實現(xiàn)方面提供參考。

CS-SE序列最終要嵌入到商用掃描儀，必須對于不同部位分析出最佳的欠采樣因子，并直接固定重建程序中不同部位對應(yīng)的稀疏方法以及最佳迭代次數(shù)，這樣才能既實現(xiàn)加速掃描，獲得重建質(zhì)量好的MR圖像。我們下一步將研究更優(yōu)的壓縮感知重建算法來適應(yīng)更大的欠采樣加速因子，并將CS重建算法集成到MR掃描儀上，進(jìn)行在線重建，爭取將壓縮感知成像應(yīng)用到臨床成像序列中。

感謝深圳貝斯達(dá)醫(yī)療器械公司、湖北省醫(yī)學(xué)信息分析及腫瘤診療重點實驗室，以及國家民委腦認(rèn)知實驗室的支持！