一種基于高性能DSP的MRI梯度計算模塊設(shè)計

潘文宇，張富，羅海，周荷琴

中國科學(xué)技術(shù)大學(xué)自動化系，合肥，230027

磁共振成像（Magnetic Resonance Imaging，MRI）具有分辨率高、成像參數(shù)多、組織敏感性豐富、可任意層面斷層成像和對人體無電離輻射等特點(diǎn)[1]，作為一種先進(jìn)的醫(yī)學(xué)成像方式被廣泛用于臨床診斷中。在MRI系統(tǒng)中，對成像物體的編碼定位是由梯度子系統(tǒng)在空間中施加X、Y、Z三維線性梯度磁場來實現(xiàn)的。梯度子系統(tǒng)由梯度計算模塊、數(shù)模轉(zhuǎn)換模塊、功率放大器和梯度線圈等組成。梯度計算是其中的關(guān)鍵模塊，作用是根據(jù)用戶設(shè)定的成像序列和參數(shù)，計算得到實際輸出的梯度波形，它需要具備高速和高精度運(yùn)算的能力，并能實現(xiàn)角度旋轉(zhuǎn)變換、渦流補(bǔ)償和勻場補(bǔ)償?shù)裙δ?。因為梯度磁場的線性度與磁共振成像定位的精確度直接相關(guān)，定位不準(zhǔn)會在圖像中產(chǎn)生偽影，所以梯度計算模塊設(shè)計的優(yōu)劣對圖像質(zhì)量有著重要的影響。

早期的梯度計算模塊由微處理器、數(shù)字乘法器和用于波形預(yù)加重的模擬電路組成[2]，由于在設(shè)計中使用了模擬器件，系統(tǒng)穩(wěn)定性較差。近年來國內(nèi)研究機(jī)構(gòu)在梯度模塊的研發(fā)上也取得了一定的進(jìn)展，我們實驗室曾提出一種用計算機(jī)進(jìn)行梯度計算的方法[3]，在成像掃描前將序列中的全部梯度波形一次性計算完成，下載到梯度發(fā)生板卡的板載內(nèi)存中，再經(jīng)DAC輸出到梯度放大器。這種方案實現(xiàn)了梯度計算的全數(shù)字化，簡化了硬件設(shè)計，穩(wěn)定性好，但對于板載內(nèi)存的容量要求很高，且不利于參數(shù)的實時調(diào)整。文獻(xiàn)[4]提出了一種基于DSP的梯度計算方案，具有設(shè)計簡單、算法易于調(diào)整的優(yōu)點(diǎn)，但受制于DSP的順序執(zhí)行架構(gòu)的影響，計算速度不高。文獻(xiàn)[5]提出的基于FPGA的設(shè)計，在一定程度上提升了運(yùn)算能力，但在運(yùn)算過程中使用的是16bit和24bit的數(shù)據(jù)，精度不夠高，且使用FPGA實現(xiàn)指數(shù)等浮點(diǎn)運(yùn)算不如DSP編程方便。

針對上述問題，本文提出了一種基于高性能DSP的MRI梯度計算模塊的設(shè)計。該設(shè)計具有以下特點(diǎn)：充分運(yùn)用單指令多數(shù)據(jù)（Single Instruction Multiple Data，SIMD）功能，在算法上對計算程序進(jìn)行并行優(yōu)化，使得梯度計算的速度明顯提升；在計算過程中使用32bit/40bit的數(shù)據(jù)精度，提高了運(yùn)算精度；梯度波形的預(yù)加重可根據(jù)用戶需要靈活配置時間常數(shù)和幅度參數(shù)的個數(shù)，通常使用4～6組預(yù)加重參數(shù)，如有特殊需求，可以支持10組或更多的參數(shù)設(shè)置。此外，還支持多種梯度波形，包括常見的數(shù)學(xué)函數(shù)波形以及用戶自定義的波形。

1 設(shè)計方法

1.1 硬件結(jié)構(gòu)

我們設(shè)計的MRI譜儀梯度子系統(tǒng)的硬件結(jié)構(gòu)如圖1所示，以DSP芯片為核心構(gòu)成的梯度計算模塊是其中的關(guān)鍵部件。DSP從脈沖序列控制模塊獲取梯度時序和相關(guān)參數(shù)，經(jīng)過計算后得到X、Y、Z三路梯度數(shù)據(jù)和B0通道的補(bǔ)償數(shù)據(jù)。DSP輸出的數(shù)據(jù)以并行方式傳輸給FPGA，進(jìn)行各通道數(shù)據(jù)的分流和并-串轉(zhuǎn)換后再輸出給DAC，轉(zhuǎn)換得到的模擬信號經(jīng)梯度放大器后驅(qū)動梯度線圈，產(chǎn)生所需的梯度磁場和B0勻場補(bǔ)償。采用這樣的硬件結(jié)構(gòu)，梯度計算模塊使用單片DSP實現(xiàn)，充分利用了DSP運(yùn)算速度快、精度高和編程方便的優(yōu)點(diǎn)，且算法通用性好，具有良好的可移植性。

圖1 MRI梯度子系統(tǒng)的硬件結(jié)構(gòu)圖Fig.1 Hardware structure diagram of the MRI gradient subsystem

梯度計算模塊中選用了Analog Device公司的ADSP-21369[6]。這是一款SHARC架構(gòu)的高性能DSP，具有400MHz主頻，基于SIMD內(nèi)核，支持32bit定點(diǎn)和32bit/40bit浮點(diǎn)運(yùn)算，峰值性能可達(dá)2.4GFLOPS。ADSP-21369片上存儲器包括2Mb的SRAM和6Mb的ROM，還提供了32bit的外部存儲器接口，可以方便地與SDRAM、SRAM和FLASH等存儲器直接相連。我們在DSP上外接了2Mb的FLASH和512Mb的SDRAM，其中FLASH用來存儲DSP的啟動程序，SDRAM用來緩存梯度參數(shù)和梯度計算的結(jié)果，梯度計算程序則存儲在DSP內(nèi)部的ROM中。

1.2 算法設(shè)計

在磁共振成像中，用戶編寫成像序列時使用的層選梯度(Slice Gradient,Gs)、相位梯度(Phase Gradient,Gp)和讀出梯度(Read Gradient,Gr)稱為邏輯梯度，而梯度子系統(tǒng)實際輸出的三維空間X、Y、Z通道的梯度Gx、Gy、Gz稱為物理梯度。梯度計算模塊的算法流程如圖2所示。

圖2 梯度計算模塊的算法流程圖Fig.2 Algorithm flow chart of the gradient calculation module

梯度計算模塊的輸入為邏輯梯度波形描述、梯度時序、幅度調(diào)節(jié)、旋轉(zhuǎn)變換、預(yù)加重、一階勻場等參數(shù)信息。梯度計算首先計算出邏輯梯度Gs、Gp和Gr的具體波形，再經(jīng)過邏輯比例縮放、旋轉(zhuǎn)變換、物理梯度增益調(diào)節(jié)而得到物理梯度基礎(chǔ)波形；然后再對波形進(jìn)行預(yù)加重和一階勻場處理；為保障系統(tǒng)的穩(wěn)定性和安全性，還需要對波形的幅度進(jìn)行溢出檢查，如果超出最大量程則進(jìn)行幅度限制，并產(chǎn)生中斷請求；最后得到實際輸出的物理梯度Gx、Gy、Gz和B0通道的勻場補(bǔ)償。

1.2.1 物理梯度的基礎(chǔ)波形計算

物理梯度的基礎(chǔ)波形計算包括邏輯梯度波形計算、邏輯梯度比例縮放、旋轉(zhuǎn)變換和物理梯度增益調(diào)節(jié)幾個部分。梯度計算模塊首先根據(jù)輸入的波形描述和時序要求計算得到MRI序列所需的邏輯梯度波形Gs、Gp和Gr，計算公式如下：

其中n為當(dāng)前時刻，psi、ppi和pri為邏輯梯度描述參數(shù)i=1,2, ...k。

然后根據(jù)用戶需要，使用設(shè)定的比例對邏輯梯度進(jìn)行縮放，比例因子記為Ls、Lp和Lr，分別對應(yīng)層選、相位和讀出梯度。邏輯梯度縮放的計算公式如下：

式中L為邏輯比例矩陣，通過改變Ls可以調(diào)節(jié)成像的選層厚度，改變Lp和Lr則可以改變相位編碼方向或者讀出方向(也稱頻率編碼方向)的視野(Field of View，F(xiàn)OV)。使用這種比例縮放的方式，用戶可以在序列編寫完成后方便地調(diào)節(jié)層厚和FOV，而不用在序列中重新計算梯度。

解剖學(xué)中定義了三個標(biāo)準(zhǔn)斷面：橫斷面、矢狀面和冠狀面。MRI的特點(diǎn)之一是可以任意方向斷層成像。當(dāng)成像層面不是標(biāo)準(zhǔn)斷面時，物理梯度與邏輯梯度的方向不一致，則需要將邏輯梯度進(jìn)行旋轉(zhuǎn)變換，投影到物理梯度方向上。用戶給出的旋轉(zhuǎn)變換參數(shù)為邏輯梯度方向與物理梯度方向之間的夾角，設(shè)它們的夾角如表1所示。

表1 邏輯梯度與物理梯度方向之間的夾角Tab.1 Angles between logical and physical gradient directions

根據(jù)所給的夾角參數(shù)，可以計算出梯度的旋轉(zhuǎn)投影矩陣，進(jìn)而計算出相應(yīng)的物理梯度，如式(3)所示，其中R為旋轉(zhuǎn)投影矩陣。

以上的梯度計算均基于X、Y、Z三個通道梯度強(qiáng)度相等的假設(shè)，即給予相同的輸入數(shù)據(jù)能產(chǎn)生同樣大小的梯度磁場。但在實際系統(tǒng)中，由于工藝和環(huán)境等因素的差別，三個通道的強(qiáng)度存在差異性，因此我們需要對物理梯度的強(qiáng)度進(jìn)行增益調(diào)節(jié)，以平衡不同的梯度通道，計算公式為：

式中P為增益調(diào)節(jié)矩陣，Px、Py、Pz分別為X、Y、Z通道的增益調(diào)節(jié)因子。

由式（2）、（3）、（4）可得：

其中矩陣M稱為總轉(zhuǎn)換矩陣，其計算公式為：

在實際設(shè)計中，我們先根據(jù)輸入?yún)?shù)將M計算出來，再與邏輯梯度[Gs,Gp,Gr]T相乘得到物理梯度的基礎(chǔ)波形。這與分別將三個矩陣P，R，L和邏輯梯度相乘比較，可以減少乘法運(yùn)算的次數(shù)。另外，在矩陣乘法運(yùn)算中，為了充分利用DSP在硬件結(jié)構(gòu)上對循環(huán)結(jié)構(gòu)、循環(huán)尋址、多操作數(shù)尋址（在一個時鐘周期內(nèi)同時存取兩個操作數(shù)）、多運(yùn)算指令（在一個時鐘周期內(nèi)同時執(zhí)行一次加法、一次減法和一次乘法）的支持，我們使用匯編語言編程來直接實現(xiàn)矩陣乘法。與用C語言編程相比，采用匯編的實現(xiàn)方式優(yōu)化了程序的執(zhí)行效率，極大地提高了計算速度。

1.2.2 預(yù)加重和勻場補(bǔ)償

由于外部磁場環(huán)境的影響，將物理梯度的基礎(chǔ)波形直接輸出并不能夠在實際系統(tǒng)中得到理想的梯度磁場。根據(jù)法拉第電磁感應(yīng)定律，隨時間變化的磁場在其周圍的金屬體內(nèi)會激發(fā)變化的渦旋電場，金屬中的載流子在該電場的作用下運(yùn)動而形成電流，稱為渦流[7]。在MRI系統(tǒng)中，在梯度線圈上施加梯度電流產(chǎn)生時變的梯度磁場時，必然在周圍的導(dǎo)體中感應(yīng)出渦流，渦流的存在會嚴(yán)重影響梯度磁場的變化，使其波形產(chǎn)生畸變，如果不進(jìn)行修正，則會在成像中形成偽影，影響圖像質(zhì)量。

根據(jù)渦流的L-R電路模型[8]，磁場中感應(yīng)出的渦流場可以表示為：

上式中的e(t)為渦流的沖激響應(yīng)，可以表示為多個指數(shù)函數(shù)的累加：

式中H(t)為階躍函數(shù)，如式(9)所示，Ai為幅度常數(shù)，Ti為時間常數(shù)，m是所包含的渦流環(huán)路個數(shù)。

為了消除渦流的影響，在梯度計算中需要對梯度電流進(jìn)行預(yù)加重處理，在梯度上升沿和下降沿施加與渦流作用相反的過電流。將連續(xù)系統(tǒng)離散化后由式(7)和式(8)可以推導(dǎo)出預(yù)加重電流的形式為：

其中yi(n)為補(bǔ)償單個渦流環(huán)路的預(yù)加重電流，x(n)為輸入的物理梯度基礎(chǔ)波形，將p(n)疊加到x(n)上即可得到預(yù)加重后的波形。一般來說，m的值越大，即幅度和時間常數(shù)的個數(shù)越多，渦流補(bǔ)償?shù)男Ч胶茫\(yùn)算量也越大，通常情況下m的取值為4～6。當(dāng)用戶有特殊需求時，本設(shè)計的運(yùn)算能力足以支持m取10或更大的值。

在實際的MRI系統(tǒng)中，X、Y、Z任何一路梯度信號的輸出均會對主磁場B0產(chǎn)生影響，使主磁場產(chǎn)生一定的偏移，偏移的程度與梯度脈沖的強(qiáng)度、渦流的影響以及線圈的設(shè)計等因素有關(guān)。為了對梯度磁場造成的主磁場偏移進(jìn)行補(bǔ)償，我們需要計算補(bǔ)償信號Bp，用以驅(qū)動B0勻場線圈。本文使用了一個補(bǔ)償函數(shù)pB(n)來定義Bp，在默認(rèn)的情況下pB(n)定義為X、Y、Z三路梯度渦流補(bǔ)償?shù)暮?，如下式所示?/p>

其中px(n)、py(n)、pz(n)的計算方法如式(10)所示。這樣的設(shè)計主要是考慮了渦流的影響。近年來，出現(xiàn)了一些對B0補(bǔ)償?shù)男卵芯縖9]。本設(shè)計支持用戶修改式(11)中的幅度和時間常數(shù)，也提供了用戶自定義補(bǔ)償函數(shù)pB(n)的接口，以便用戶根據(jù)實際情況考慮更多影響因素和優(yōu)化補(bǔ)償效果。

在理想的磁共振成像中，主磁場的靜磁場強(qiáng)度應(yīng)該在成像范圍內(nèi)具有嚴(yán)格的均勻性，但是由于制造工藝和周圍環(huán)境的影響，靜磁場強(qiáng)度的均勻性往往存在一定的偏差，因此需要在梯度線圈中疊加直流偏置，起到一階勻場的作用。記在X、Y、Z和B0通道施加的直流勻場分量為Sx、Sy、Sz和SB，則上述對物理梯度波形的預(yù)加重和一階勻場補(bǔ)償?shù)挠嬎懔鞒炭梢跃C合用圖3表示，相應(yīng)的計算公式為：

在上述預(yù)加重和勻場補(bǔ)償中，需要進(jìn)行大量的指數(shù)乘法和加法運(yùn)算，如果使用SISD（單指令單數(shù)據(jù)）方式來編寫程序，則計算時間開銷較大。在實際設(shè)計中，我們充分利用了ADSP-21369的SIMD內(nèi)核，在SIMD模式下進(jìn)行預(yù)加重數(shù)值的計算。這樣兩個渦流環(huán)路的預(yù)加重可并行計算，明顯縮短了計算時間，改善了算法效率。

圖3 梯度波形的預(yù)加重和一階勻場補(bǔ)償流程圖Fig.3 Flow chart of pre-emphasis and first-order shimming for gradient waveforms

梯度計算模塊的最后一個環(huán)節(jié)是對每個通道將要輸出的梯度和勻場波形的幅值進(jìn)行檢查，如果超過了梯度和勻場功放的輸入量程，則將該通道的輸出限定為滿量程的最大值，同時向上層的序列控制器提出中斷請求。設(shè)梯度功放的量程為[-Gm,Gm]，Gm＞0，以X通道為例，則梯度計算模塊最后的輸出如式(16)所示。

2 實驗結(jié)果及分析

為了驗證本文設(shè)計的梯度計算模塊算法的有效性，我們使用ADI公司的VisualDSP++ 5.0來完成梯度計算程序的開發(fā)和調(diào)試，并在ADSP-21369的開發(fā)板上進(jìn)行了實驗。

梯度計算的實驗結(jié)果以圖形化的方式輸出。因為實際系統(tǒng)中梯度波形會由驅(qū)動接口電路將梯度電流轉(zhuǎn)換為電壓方式輸出，所以縱坐標(biāo)以電壓為單位。特殊梯度波形的輸出如圖4所示，(a)為正弦波形，(b) 的梯度下降沿包含一段雙曲線波形。圖5顯示的是單個物理梯度波形補(bǔ)償前后的對比，可見預(yù)加重和一階勻場補(bǔ)償算法對單通道數(shù)據(jù)的有效性。圖6顯示的是一個自旋回波(Spin Echo，SE)序列的物理梯度各通道輸出波形，圖中的Gx、Gy、Gz分別對應(yīng)Gr、Gp、Gs，Bp為B0補(bǔ)償通道的輸出，渦流補(bǔ)償使用4組預(yù)加重參數(shù)，幅度常數(shù)為A1=4.5，A2=2.5，A3=1.4，A4=1.2，時間常數(shù)為T1=1.0ms，T2=0.1ms，T3=0.01ms，T4=0.001ms，一階勻場的參數(shù)設(shè)定為Sx= -1mV，Sy= 0，Sz=2mV，SB=1mV。圖6的結(jié)果表明，本算法可以計算MRI成像序列中所有通道的梯度數(shù)據(jù)，并施加相應(yīng)的預(yù)加重和補(bǔ)償處理，在功能上滿足實際MRI系統(tǒng)對梯度計算模塊的要求。

圖4 特殊梯度波形Fig.4 Special gradient waveforms

圖5 預(yù)加重和勻場補(bǔ)償前后的梯度波形Fig.5 Gradient waveforms before/after pre-emphasis and shimming

圖6 SE序列的梯度各通道輸出波形圖Fig.6 Waveform of each gradient channel in SE sequence

為了提高計算速度，我們在編程實現(xiàn)時對梯度計算程序的各個環(huán)節(jié)進(jìn)行了優(yōu)化。實驗中將DSP的主頻設(shè)為400 MHz，計算了所有通道的2000個波形點(diǎn)，再取平均得到計算所有通道一個波形點(diǎn)的時間，優(yōu)化前后的計算時間對比如表2所示。實驗結(jié)果表明，優(yōu)化后的整體計算時間減少了約90%，計算速度大幅提升，實際計算一個梯度波形數(shù)據(jù)點(diǎn)的時間約為0.332μs，而目前所用的PCM1704(DAC)轉(zhuǎn)換一個數(shù)據(jù)點(diǎn)的周期為1.302μs(對應(yīng)768 KHz的轉(zhuǎn)換率)，這表示在DSP計算完一個數(shù)據(jù)點(diǎn)后，還有0.970μs的空閑時間，這給用戶提供了足夠的升級空間：可以更新梯度計算算法，計算更多的波形預(yù)加重參數(shù)(可支持10組以上)或者自定義新的計算功能，也可以隨著電子技術(shù)的發(fā)展升級使用轉(zhuǎn)化率更高的DAC(最高支持3 MHz)。另一方面，梯度計算模塊在輸入和內(nèi)部計算的過程中使用32 bit/40 bit的數(shù)據(jù)精度，最后輸出時才轉(zhuǎn)換成PCM1704所需的24 bit數(shù)據(jù)，這使計算結(jié)果的精度更高。因此，本文的設(shè)計在速度和精度上都能滿足實際MRI系統(tǒng)的需求。

表2 程序優(yōu)化前后的梯度計算時間對比Tab.2 Comparison of gradient calculation time before/after program optimization

3 結(jié)論

本文提出了一種基于SHARC DSP的MRI梯度計算模塊的設(shè)計方法。該模塊使用一片高性能DSP完成梯度計算有關(guān)的旋轉(zhuǎn)、預(yù)加重和勻場補(bǔ)償?shù)热抗δ?，能夠方便地與序列控制模塊和梯度輸出模塊相連接。在DSP開發(fā)板上進(jìn)行的實驗結(jié)果表明，所設(shè)計的模塊在功能、速度和精度上都能滿足MRI梯度計算的要求，且還有足夠的余量可供用戶實現(xiàn)更多的自定義功能和后續(xù)的系統(tǒng)升級，為研制數(shù)字化MRI譜儀提供了一種實用的梯度計算解決方案。

[1]Haacke EM, Brown RW, Thompson MR et al.Magnetic Resonance Imaging: physical principles and sequence design [M].John Wiley& Sons Publisher, 1999.

[2]Gach HM, Lowe IJ, Madio DP et al.A programmable pre-emphasis system [J].Magnetic Resonance in Medicine, 1998, 40: 427-431.

[3]Zhengmin Liu, Cong Zhao, Heqin Zhou et al.A Novel Digital magnetic resonance imaging spectrometer[A].Annual International Conference of the IEEE Engineering in Medicine and Biology -Proceedings[C], 2006, 280-283.

[4]戴祎棟, 寧瑞鵬, 李鯁穎.基于DSP技術(shù)的梯度波形發(fā)生器 [J].波譜學(xué)雜志, 2009, 26(1): 44-50

[5]肖亮, 湯偉男, 王為民.基于單片F(xiàn)PGA的磁共振成像梯度計算模塊 [J].波譜學(xué)雜志, 2010, 27(2): 163-171.

[6]ADSP-21367/ADSP-21368/ADSP-21369 SHARC Processors Data Sheet.Rev.E.Analog Devices, Inc., 2009.

[7]劉正敏, 周荷琴, 武海澄.磁共振成像系統(tǒng)的一種快速渦流補(bǔ)償方法 [J], 中國醫(yī)療器械雜志, 2005, 29(6):410-413.

[8]Matt A.Bernstein, Kevin Franklin King, Xiaohong Joe Zhou,Handbook of MRI Pulse Sequences [M].Academic Press, 2004.

[9]Terence W N, Scott McIntyre, Douglas L, et al.Compensation of gradient-induced magnetic field perturbations [J].Journal of Magnetic Resonance, 2008, 192:209-217.