一种基于高性能DSP的MRI梯度计算模块设计

潘文宇，张富，罗海，周荷琴

中国科学技术大学自动化系，合肥，230027

磁共振成像（Magnetic Resonance Imaging，MRI）具有分辨率高、成像参数多、组织敏感性丰富、可任意层面断层成像和对人体无电离辐射等特点[1]，作为一种先进的医学成像方式被广泛用于临床诊断中。在MRI系统中，对成像物体的编码定位是由梯度子系统在空间中施加X、Y、Z三维线性梯度磁场来实现的。梯度子系统由梯度计算模块、数模转换模块、功率放大器和梯度线圈等组成。梯度计算是其中的关键模块，作用是根据用户设定的成像序列和参数，计算得到实际输出的梯度波形，它需要具备高速和高精度运算的能力，并能实现角度旋转变换、涡流补偿和匀场补偿等功能。因为梯度磁场的线性度与磁共振成像定位的精确度直接相关，定位不准会在图像中产生伪影，所以梯度计算模块设计的优劣对图像质量有着重要的影响。

早期的梯度计算模块由微处理器、数字乘法器和用于波形预加重的模拟电路组成[2]，由于在设计中使用了模拟器件，系统稳定性较差。近年来国内研究机构在梯度模块的研发上也取得了一定的进展，我们实验室曾提出一种用计算机进行梯度计算的方法[3]，在成像扫描前将序列中的全部梯度波形一次性计算完成，下载到梯度发生板卡的板载内存中，再经DAC输出到梯度放大器。这种方案实现了梯度计算的全数字化，简化了硬件设计，稳定性好，但对于板载内存的容量要求很高，且不利于参数的实时调整。文献[4]提出了一种基于DSP的梯度计算方案，具有设计简单、算法易于调整的优点，但受制于DSP的顺序执行架构的影响，计算速度不高。文献[5]提出的基于FPGA的设计，在一定程度上提升了运算能力，但在运算过程中使用的是16bit和24bit的数据，精度不够高，且使用FPGA实现指数等浮点运算不如DSP编程方便。

针对上述问题，本文提出了一种基于高性能DSP的MRI梯度计算模块的设计。该设计具有以下特点：充分运用单指令多数据（Single Instruction Multiple Data，SIMD）功能，在算法上对计算程序进行并行优化，使得梯度计算的速度明显提升；在计算过程中使用32bit/40bit的数据精度，提高了运算精度；梯度波形的预加重可根据用户需要灵活配置时间常数和幅度参数的个数，通常使用4～6组预加重参数，如有特殊需求，可以支持10组或更多的参数设置。此外，还支持多种梯度波形，包括常见的数学函数波形以及用户自定义的波形。

1 设计方法

1.1 硬件结构

我们设计的MRI谱仪梯度子系统的硬件结构如图1所示，以DSP芯片为核心构成的梯度计算模块是其中的关键部件。DSP从脉冲序列控制模块获取梯度时序和相关参数，经过计算后得到X、Y、Z三路梯度数据和B0通道的补偿数据。DSP输出的数据以并行方式传输给FPGA，进行各通道数据的分流和并-串转换后再输出给DAC，转换得到的模拟信号经梯度放大器后驱动梯度线圈，产生所需的梯度磁场和B0匀场补偿。采用这样的硬件结构，梯度计算模块使用单片DSP实现，充分利用了DSP运算速度快、精度高和编程方便的优点，且算法通用性好，具有良好的可移植性。

图1 MRI梯度子系统的硬件结构图Fig.1 Hardware structure diagram of the MRI gradient subsystem

梯度计算模块中选用了Analog Device公司的ADSP-21369[6]。这是一款SHARC架构的高性能DSP，具有400MHz主频，基于SIMD内核，支持32bit定点和32bit/40bit浮点运算，峰值性能可达2.4GFLOPS。ADSP-21369片上存储器包括2Mb的SRAM和6Mb的ROM，还提供了32bit的外部存储器接口，可以方便地与SDRAM、SRAM和FLASH等存储器直接相连。我们在DSP上外接了2Mb的FLASH和512Mb的SDRAM，其中FLASH用来存储DSP的启动程序，SDRAM用来缓存梯度参数和梯度计算的结果，梯度计算程序则存储在DSP内部的ROM中。

1.2 算法设计

在磁共振成像中，用户编写成像序列时使用的层选梯度(Slice Gradient,Gs)、相位梯度(Phase Gradient,Gp)和读出梯度(Read Gradient,Gr)称为逻辑梯度，而梯度子系统实际输出的三维空间X、Y、Z通道的梯度Gx、Gy、Gz称为物理梯度。梯度计算模块的算法流程如图2所示。

图2 梯度计算模块的算法流程图Fig.2 Algorithm flow chart of the gradient calculation module

梯度计算模块的输入为逻辑梯度波形描述、梯度时序、幅度调节、旋转变换、预加重、一阶匀场等参数信息。梯度计算首先计算出逻辑梯度Gs、Gp和Gr的具体波形，再经过逻辑比例缩放、旋转变换、物理梯度增益调节而得到物理梯度基础波形；然后再对波形进行预加重和一阶匀场处理；为保障系统的稳定性和安全性，还需要对波形的幅度进行溢出检查，如果超出最大量程则进行幅度限制，并产生中断请求；最后得到实际输出的物理梯度Gx、Gy、Gz和B0通道的匀场补偿。

1.2.1 物理梯度的基础波形计算

物理梯度的基础波形计算包括逻辑梯度波形计算、逻辑梯度比例缩放、旋转变换和物理梯度增益调节几个部分。梯度计算模块首先根据输入的波形描述和时序要求计算得到MRI序列所需的逻辑梯度波形Gs、Gp和Gr，计算公式如下：

其中n为当前时刻，psi、ppi和pri为逻辑梯度描述参数i=1,2, ...k。

然后根据用户需要，使用设定的比例对逻辑梯度进行缩放，比例因子记为Ls、Lp和Lr，分别对应层选、相位和读出梯度。逻辑梯度缩放的计算公式如下：

式中L为逻辑比例矩阵，通过改变Ls可以调节成像的选层厚度，改变Lp和Lr则可以改变相位编码方向或者读出方向(也称频率编码方向)的视野(Field of View，FOV)。使用这种比例缩放的方式，用户可以在序列编写完成后方便地调节层厚和FOV，而不用在序列中重新计算梯度。

解剖学中定义了三个标准断面：横断面、矢状面和冠状面。MRI的特点之一是可以任意方向断层成像。当成像层面不是标准断面时，物理梯度与逻辑梯度的方向不一致，则需要将逻辑梯度进行旋转变换，投影到物理梯度方向上。用户给出的旋转变换参数为逻辑梯度方向与物理梯度方向之间的夹角，设它们的夹角如表1所示。

表1 逻辑梯度与物理梯度方向之间的夹角Tab.1 Angles between logical and physical gradient directions

根据所给的夹角参数，可以计算出梯度的旋转投影矩阵，进而计算出相应的物理梯度，如式(3)所示，其中R为旋转投影矩阵。

以上的梯度计算均基于X、Y、Z三个通道梯度强度相等的假设，即给予相同的输入数据能产生同样大小的梯度磁场。但在实际系统中，由于工艺和环境等因素的差别，三个通道的强度存在差异性，因此我们需要对物理梯度的强度进行增益调节，以平衡不同的梯度通道，计算公式为：

式中P为增益调节矩阵，Px、Py、Pz分别为X、Y、Z通道的增益调节因子。

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

其中矩阵M称为总转换矩阵，其计算公式为：

在实际设计中，我们先根据输入参数将M计算出来，再与逻辑梯度[Gs,Gp,Gr]T相乘得到物理梯度的基础波形。这与分别将三个矩阵P，R，L和逻辑梯度相乘比较，可以减少乘法运算的次数。另外，在矩阵乘法运算中，为了充分利用DSP在硬件结构上对循环结构、循环寻址、多操作数寻址（在一个时钟周期内同时存取两个操作数）、多运算指令（在一个时钟周期内同时执行一次加法、一次减法和一次乘法）的支持，我们使用汇编语言编程来直接实现矩阵乘法。与用C语言编程相比，采用汇编的实现方式优化了程序的执行效率，极大地提高了计算速度。

1.2.2 预加重和匀场补偿

由于外部磁场环境的影响，将物理梯度的基础波形直接输出并不能够在实际系统中得到理想的梯度磁场。根据法拉第电磁感应定律，随时间变化的磁场在其周围的金属体内会激发变化的涡旋电场，金属中的载流子在该电场的作用下运动而形成电流，称为涡流[7]。在MRI系统中，在梯度线圈上施加梯度电流产生时变的梯度磁场时，必然在周围的导体中感应出涡流，涡流的存在会严重影响梯度磁场的变化，使其波形产生畸变，如果不进行修正，则会在成像中形成伪影，影响图像质量。

根据涡流的L-R电路模型[8]，磁场中感应出的涡流场可以表示为：

上式中的e(t)为涡流的冲激响应，可以表示为多个指数函数的累加：

式中H(t)为阶跃函数，如式(9)所示，Ai为幅度常数，Ti为时间常数，m是所包含的涡流环路个数。

为了消除涡流的影响，在梯度计算中需要对梯度电流进行预加重处理，在梯度上升沿和下降沿施加与涡流作用相反的过电流。将连续系统离散化后由式(7)和式(8)可以推导出预加重电流的形式为：

其中yi(n)为补偿单个涡流环路的预加重电流，x(n)为输入的物理梯度基础波形，将p(n)叠加到x(n)上即可得到预加重后的波形。一般来说，m的值越大，即幅度和时间常数的个数越多，涡流补偿的效果越好，但运算量也越大，通常情况下m的取值为4～6。当用户有特殊需求时，本设计的运算能力足以支持m取10或更大的值。

在实际的MRI系统中，X、Y、Z任何一路梯度信号的输出均会对主磁场B0产生影响，使主磁场产生一定的偏移，偏移的程度与梯度脉冲的强度、涡流的影响以及线圈的设计等因素有关。为了对梯度磁场造成的主磁场偏移进行补偿，我们需要计算补偿信号Bp，用以驱动B0匀场线圈。本文使用了一个补偿函数pB(n)来定义Bp，在默认的情况下pB(n)定义为X、Y、Z三路梯度涡流补偿的和，如下式所示：

其中px(n)、py(n)、pz(n)的计算方法如式(10)所示。这样的设计主要是考虑了涡流的影响。近年来，出现了一些对B0补偿的新研究[9]。本设计支持用户修改式(11)中的幅度和时间常数，也提供了用户自定义补偿函数pB(n)的接口，以便用户根据实际情况考虑更多影响因素和优化补偿效果。

在理想的磁共振成像中，主磁场的静磁场强度应该在成像范围内具有严格的均匀性，但是由于制造工艺和周围环境的影响，静磁场强度的均匀性往往存在一定的偏差，因此需要在梯度线圈中叠加直流偏置，起到一阶匀场的作用。记在X、Y、Z和B0通道施加的直流匀场分量为Sx、Sy、Sz和SB，则上述对物理梯度波形的预加重和一阶匀场补偿的计算流程可以综合用图3表示，相应的计算公式为：

在上述预加重和匀场补偿中，需要进行大量的指数乘法和加法运算，如果使用SISD（单指令单数据）方式来编写程序，则计算时间开销较大。在实际设计中，我们充分利用了ADSP-21369的SIMD内核，在SIMD模式下进行预加重数值的计算。这样两个涡流环路的预加重可并行计算，明显缩短了计算时间，改善了算法效率。

图3 梯度波形的预加重和一阶匀场补偿流程图Fig.3 Flow chart of pre-emphasis and first-order shimming for gradient waveforms

梯度计算模块的最后一个环节是对每个通道将要输出的梯度和匀场波形的幅值进行检查，如果超过了梯度和匀场功放的输入量程，则将该通道的输出限定为满量程的最大值，同时向上层的序列控制器提出中断请求。设梯度功放的量程为[-Gm,Gm]，Gm＞0，以X通道为例，则梯度计算模块最后的输出如式(16)所示。

2 实验结果及分析

为了验证本文设计的梯度计算模块算法的有效性，我们使用ADI公司的VisualDSP++ 5.0来完成梯度计算程序的开发和调试，并在ADSP-21369的开发板上进行了实验。

梯度计算的实验结果以图形化的方式输出。因为实际系统中梯度波形会由驱动接口电路将梯度电流转换为电压方式输出，所以纵坐标以电压为单位。特殊梯度波形的输出如图4所示，(a)为正弦波形，(b) 的梯度下降沿包含一段双曲线波形。图5显示的是单个物理梯度波形补偿前后的对比，可见预加重和一阶匀场补偿算法对单通道数据的有效性。图6显示的是一个自旋回波(Spin Echo，SE)序列的物理梯度各通道输出波形，图中的Gx、Gy、Gz分别对应Gr、Gp、Gs，Bp为B0补偿通道的输出，涡流补偿使用4组预加重参数，幅度常数为A1=4.5，A2=2.5，A3=1.4，A4=1.2，时间常数为T1=1.0ms，T2=0.1ms，T3=0.01ms，T4=0.001ms，一阶匀场的参数设定为Sx= -1mV，Sy= 0，Sz=2mV，SB=1mV。图6的结果表明，本算法可以计算MRI成像序列中所有通道的梯度数据，并施加相应的预加重和补偿处理，在功能上满足实际MRI系统对梯度计算模块的要求。

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

图5 预加重和匀场补偿前后的梯度波形Fig.5 Gradient waveforms before/after pre-emphasis and shimming

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

为了提高计算速度，我们在编程实现时对梯度计算程序的各个环节进行了优化。实验中将DSP的主频设为400 MHz，计算了所有通道的2000个波形点，再取平均得到计算所有通道一个波形点的时间，优化前后的计算时间对比如表2所示。实验结果表明，优化后的整体计算时间减少了约90%，计算速度大幅提升，实际计算一个梯度波形数据点的时间约为0.332μs，而目前所用的PCM1704(DAC)转换一个数据点的周期为1.302μs(对应768 KHz的转换率)，这表示在DSP计算完一个数据点后，还有0.970μs的空闲时间，这给用户提供了足够的升级空间：可以更新梯度计算算法，计算更多的波形预加重参数(可支持10组以上)或者自定义新的计算功能，也可以随着电子技术的发展升级使用转化率更高的DAC(最高支持3 MHz)。另一方面，梯度计算模块在输入和内部计算的过程中使用32 bit/40 bit的数据精度，最后输出时才转换成PCM1704所需的24 bit数据，这使计算结果的精度更高。因此，本文的设计在速度和精度上都能满足实际MRI系统的需求。

表2 程序优化前后的梯度计算时间对比Tab.2 Comparison of gradient calculation time before/after program optimization

3 结论

本文提出了一种基于SHARC DSP的MRI梯度计算模块的设计方法。该模块使用一片高性能DSP完成梯度计算有关的旋转、预加重和匀场补偿等全部功能，能够方便地与序列控制模块和梯度输出模块相连接。在DSP开发板上进行的实验结果表明，所设计的模块在功能、速度和精度上都能满足MRI梯度计算的要求，且还有足够的余量可供用户实现更多的自定义功能和后续的系统升级，为研制数字化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技术的梯度波形发生器 [J].波谱学杂志, 2009, 26(1): 44-50

[5]肖亮, 汤伟男, 王为民.基于单片FPGA的磁共振成像梯度计算模块 [J].波谱学杂志, 2010, 27(2): 163-171.

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

[7]刘正敏, 周荷琴, 武海澄.磁共振成像系统的一种快速涡流补偿方法 [J], 中国医疗器械杂志, 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.