磁共振成像具有无辐射、对比度高等优势，所以在临床^[1]、神经科学^[2]等诸多领域得到广泛应用。梯度系统在成像序列执行期间需要快速切换，实现对任意组织进行成像。根据电磁感应定律，梯度系统的切换会在梯度线圈附近的导体内产生涡流，造成梯度磁场的畸变，结果在图像上会产生伪影等^[3,4]。

       为了减小梯度涡流的影响，可以将梯度电流波形预先进行补偿以平衡梯度磁场波形的失真，从而得到标准的梯度磁场波形。涡流波形预补偿方法，根据补偿波形的极性大致可分为过补偿与欠补偿两种。过补偿通常又称为预加重，由Jehenson等^[5]提出并使用LR模型分析了预加重的原理，通过对标准的梯度波形和补偿波形进行叠加得到补偿后波形。对于已经具有自屏蔽梯度线圈的磁共振系统，屏蔽线圈中的电流可能过大，造成实际的梯度磁场波形产生过冲，需要叠加与标准波形极性相反的欠补偿波形，对梯度场进行欠补偿^[3]。

       预补偿方法，通常采用模拟或数字两种方式来产生补偿波形。早期补偿波形通过RC模拟电路^[5]方法实现，通过调节机械电位器改变补偿波形的形状和幅值。虽然具有较高幅度分辨率，但调节较为麻烦。使用数字电位器替代机械电位器，调节更加方便，但仍属于模拟电路实现方法，由于电阻等模拟元件受环境因素影响较大，所以实际输出波形与期望的预设波形一致性较差。另外预补偿模拟电路一旦确定，就只能用于过补偿或者欠补偿，不能对补偿波形进行切换^[6]。

       随着微电子技术的发展，现多使用数字涡流预补偿方法，通过数字信号处理器（digital signal processor，DSP)^[7]、现场可编程门阵列（field programmable gate array，FPGA)^[8]等数字器件计算预补偿波形数值，然后将补偿波形叠加到标准梯度波形上，再通过数模转换器（digital-to-analog converter，DAC）进行转换得到补偿后的模拟梯度波形。基于DSP的方案计算精度较高，预补偿算法实现方便，但同时仍需使用FPGA实现DAC控制、信号通讯等配套逻辑功能，增加了开发难度和成本。黄朝晖等^[8]实现了在FPGA内实时计算波形值的数字预补偿方法，但采用单片16 bit的DAC用于输出补偿后的总波形，导致补偿波形和标准波形幅度的精度均不会高于DAC的最高精度，由于幅度精度不够，不能较好地对涡流进行补偿。同时由于实现数字预补偿花费的时间较长，使得从数据输入到输出波形之间有4 μs的传输延迟。除了提高幅度的输出精度外，由于成像序列中的射频脉冲和各个梯度波形之间都有非常严格的时序关系，因此要尽可能缩短因实现数字预补偿算法所造成的梯度波形输出延时，并确保该延时时间固定而不发生抖动，这样才能保证在实时输出梯度波形的同时，避免因每次梯度波形输出时间的不一致而造成图像伪影^[9]。王鹤等^[10]为解决梯度输出在时间上的抖动，将梯度波形发生器的输出时钟反馈送回脉冲序列发生器，以保证触发信号与波形输出时钟的相干性，但这样不仅需要增加额外的控制电路，同时导致脉冲序列控制逻辑更加复杂。

       针对以上问题，笔者基于课题组设计的高精度梯度卡，设计了具有五组时间常数的梯度波形数字预补偿算法，使用FPGA对涡流预补偿波形进行实时的高速计算，输出了高精度、低延时的数字预补偿波形，并且可以任意产生过补偿和欠补偿两种不同的波形，提高了现有磁共振系统的成像效果。

1　材料与方法

1.1　涡流预补偿方案

       涡流将造成标准梯度波形（即梯度功放的输入电压波形，或者梯度功放的输出电流波形）所产生梯度场的畸变（图1)^[11]。图1A中用实线表示输入的标准梯度波形，用点线表示实际产生的梯度磁场波形。在涡流的影响下，实际梯度磁场波形的变化时间被拉长，在梯度下降时因不能及时关闭而产生拖尾，如果不进行补偿的话将使图像产生伪影。通过对标准梯度波形进行预补偿，补偿后波形变化幅度将大于原波形，在涡流抑制下，能够得到与标准梯度磁场波形相近的效果。如图1B所示，其中实线表示具有过补偿的梯度波形，虚线表示实际产生的梯度磁场波形。

       如图2所示，预补偿后的梯度波形（图2C)可以看成是由标准波形（图2A）和补偿波形（图2B）叠加而成，其中补偿波形可由标准波形数值实时计算得到，再分别由两片DAC独立输出标准波形和补偿波形，然后在模拟域进行信号叠加，得到预补偿后的梯度电压波形。这样处理的好处是：首先，在脉冲序列设计时，完全按照标准的梯度波形进行设计，不需要考虑涡流的影响；其次，大大提高了预补偿波形的分辨率，能够精确产生具有预补偿的模拟电压波形。

       预补偿波形采用多组e指数求和的形式来表示^[12,13]，公式如下：

       其中ai是第i组的幅值常数，τi是第i组的时间常数，k是e指数组合的组数，本研究采用五组参数进行拟合。

       式（1)为模拟表达形式，为实现FPGA运算，将补偿波形表示成数字形式如下：

       式(2)中的时间参数，Δt为数字信号输出的时间间隔，通过Δt将模拟波形的连续时间变量t转换为数字波形离散时间变量n，即n=t/Δt。式（2)中的幅值参数Ai与式(1)的幅值常数ai保持一致。

       考虑单e指数情况，X [n]表示标准的输入梯度波形，Yi[n］表示第i组输出波形信号，可求得差分方程为：

       其中，Ai和Ei分别表示幅值和时间参数。通过式(3）可以看出，任何时刻的输出信号，都是对前一时刻的输出和输入信号计算得到，这个过程可以通过FPGA设计相应的逻辑来完成，从而达到实时计算的目的。将五组单e指数的预补偿波形值Yi[n］如式(2）相加，则可得到最终的预补偿波形值Y[n]。这就是图2中补偿波形由标准波形实时求解的计算过程，下一节介绍本方案如何通过FPGA实现此过程。

点击查看大图

图1  预补偿前后梯度波形示意图。A：涡流影响波形示意图；B：预补偿波形示意图
Fig. 1  Diagram of gradient waveform before and after eddy-current precompensation. Figure A shows the effect of eddy-current, where solid line represents standard trapezoidal gradient waveform and dot line represents the actual gradient field. Figure B shows precompensated waveform, where solid line represents precompensated waveform and dashed line represents the net gradient field.

点击查看大图

图2  波形叠加示意图。A：标准波形；B：涡流预补偿波形；C：涡流补偿后波形
Fig. 1  Superposition of waveforms. Add standard waveform (A) and eddy-current precompensation waveform (B) to generate compensated waveform (C).

点击查看大图

点击查看大图

点击查看大图

1.2　涡流预补偿方案的FPGA实现

       本设计使用FPGA实现涡流预补偿计算过程以及与计算模块相连的逻辑控制，如DAC控制、标准波形数值读取等。通过Simulink和DSP Builder搭建预补偿数字运算模型，然后直接生成硬件描述语言(hardware description language，HDL)，从而简化了FPGA开发过程^[8]。

       图3为式(3）的单e指数预补偿算法实现框图。幅值参数Ai、时间参数Ei、标准波形X [n]是单e指数计算模块的输入端口，补偿波形Yi [n]是此模块的输出端口，这些端口的位宽是根据实际希望达到的预补偿精度而设定。在DSP Builder中需设定各个运算信号的位宽，图3中使用双线箭头表示数据格式的转换，S[xx].[xx]代表有符号小数，US [xx]. [xx］代表无符号小数，小数点前后分别表示整数部分与小数部分。因为本研究使用的是20 bit DAC器件，所以标准波形X [n]数据格式是20位表示的有符号定点纯小数，X [n］范围为(-1，1)。幅值参数Ai极性的正和负分别表示预补偿是属于过补偿还是欠补偿，Ai的大小与补偿波形幅值的大小成正比，本文用百分比的形式表示Ai，其数据格式为24位表示的有符号定点纯小数，范围是(-1，1)。预补偿时间常数τi的大小反映补偿波形的衰减快慢，τi越大，补偿波形衰减就越慢，用于补偿持续时间更长的涡流。τi一般在1～1000 ms量级，将此代入时间参数公式对Ei进行计算，并采用32位无符号定点纯小数进行表示。

       在定点计算中，乘法与加法会增加位宽，需要对数据进行截位以保证计算的持续运行，如图3中的两端双线箭头S[4].[64]∶S[1].[32]所示。在磁共振成像序列设计中，出于对被试安全等因素考虑，梯度场切换率存在上限^[14]，标准梯度波形X [n]的切换率不能剧烈变化。实际上相邻标准梯度波形值的变化不会超过满幅度波形的50%，即|X[n]-X [n-1]|<0.5，又因为Ai与Ei的取值范围限制，所以Yi [n]的取值范围为(-1，1)，因此预补偿计算中数据的高位截断并不会对计算结果产生影响。此外，计算过程的数据位宽远大于最终波形数据的输出位宽，所以数据的低位截断对计算结果精度影响较小。在DSP Builder中通过延迟模块Z^-1实现X [n]的标准波形延迟，设置求和模块的流水线级数实现补偿波形Yi [n-1]的延迟，最终实现式(3）的实时计算。

       图4是五组e指数涡流预补偿计算框图，其中虚线框部分代表图3虚线框中的单e指数计算。除了标准波形和五组补偿参数外，输入端口还有匀场偏移量和两个DAC的零点漂移校正量。DAC_X和DAC_Y为预补偿计算模块的输出端口，直接送给后续的DAC器件产生模拟波形。将五组预补偿波形在数字域求和，范围是(-5，5)，通过乘上因子0.2归一化至(-1，1)。这样补偿波形DAC_Y也使用20位有符号定点纯小数表示，与原始波形数据格式统一，充分发挥DAC性能。由于完成全部预补偿计算需要两个时钟周期，所以将标准梯度波形也延迟两个周期以确保两路波形时序正确，图中用X'[n]表示延时后的标准梯度波形。此外，在送入DAC之前，两路波形均需加上独立的漂移校正量，用于微调DAC器件的零点漂移。对于DAC_X还需加上匀场偏移量，用于产生输入梯度线圈的匀场电流以实现匀场^[15]。因为漂移校正量和匀场偏移量在一次成像实验中均为定值，所以仅需在梯度波形输出前求和即可实现。在本设计中，使用的DAC型号是Analog Devices公司的AD5791，其位宽是20 bit，输出建立时间为1 μs。为了提高梯度波形的时间分辨率，设计中DAC的更新频率始终保持1 MHz^[16]。

       补偿波形的幅值相较于原始波形较小^[3]，按照实际测试，本文设定其不超过梯度满幅度电压波形的10%，即通过调节DAC后的分压电路，将DAC_Y的最大输出调至满幅度电压波形的10%。DAC_X的输出电压范围是±4.5 V，DAC_Y的输出电压范围是±0.5 V，均使用独立的20 bit DAC。这样，对于预补偿部分，其补偿波形的幅度分辨率为1.0 V/2²⁰≈0.954 μV。按照传统数字涡流预补偿设计，若使用一片DAC产生整个预补偿波形，由于10.0 V/2^23.3≈0.954 μV，因此为了达到同样的幅度分辨率，需要的DAC位数约为23.3 bit。也就是说，本设计采用两片相同型号的20 bit DAC器件，通过模拟信号叠加的办法，获得了优于23 bit的预补偿波形幅度分辨力性能。

       本方案梯度补偿数值计算模块使用DSP Builder搭建多级流水线(pipeline）结构，可以合理分配FPGA资源。标准梯度波形X [n]输入更新后，需要两个时钟周期完成计算，然后输出补偿波形Y[n]，之后每个时钟周期均更新并输出Y[n]。

点击查看大图

图3  单e指数涡流预补偿算法框图
Fig. 3  Scheme of eddy precompensation described by single exponential decay term.

点击查看大图

图4  五组e指数涡流预补偿计算框图
Fig. 4  Scheme of eddy precompensation described by 5 exponential decay terms.

1.3　实时计算与同步触发

       本方案基于实验室自行研制的磁共振谱仪系统，预先将梯度波形数据存储在梯度卡上，然后由序列控制卡产生触发信号，来实时读取梯度波形数据，并进行实时预补偿计算。由于DAC器件的更新速率为1 μs，如果使用1 MHz固定时钟来更新DAC输出，由于在脉冲序列执行时所产生的梯度触发信号可能会在一个时钟周期(1 μs)中的任意时刻到来，这样将导致不确定的补偿波形输出延时（图5)，两次触发距下一个更新时钟上升沿的时间t1和t2不等，差距最大为1 μs，这样就造成了每次实际输出梯度波形信号的时间不一致，这会在图像中产生伪影。

       所以在FPGA设计中，每当触发信号上升沿来临时，都应重置DAC输出时钟(1 MHz)，这样DAC输出时钟与触发信号之间的相干性可以得到改善。实际上，当触发信号来临时，并不会立刻更新DAC输出时钟的上升沿，而应为预补偿计算预留一段时间。

       由于FPGA对预补偿波形计算采用的是二级流水线结构设计，在标准梯度波形数据X [n］更新两个时钟周期后，补偿波形Y[n]开始更新，之后每个时钟周期Y[n]均更新一次。为了缩短预补偿计算耗时，预补偿计算使用40 MHz时钟，这样预补偿计算耗时仅为50 ns。同时采用此40 MHz时钟实现对DAC器件的数据写入和相关逻辑控制，数据写入后同步复位DAC的1 MHz更新时钟。由于DAC更新时钟是对40 MHz时钟分频得到，同时从预补偿数据计算开始到第一个DAC更新时钟的产生，这之间的时间间隔每次都是相等的，因此其与触发信号间延时的不确定度降低到了25 ns以内，这样每次模拟波形的输出时间都非常稳定。

       需要注意的是，由于DAC更新频率与预补偿计算数值更新频率的倍数关系，对于实时计算得到的Y[n]，每40个数据只有一个会由DAC转化为模拟波形。另外，由于预补偿实时计算时钟为40 MHz，因此需要将时间参数计算公式中的时间间隔Δt设置为25 ns。

点击查看大图

图5  延时误差示意图
Fig. 5  Diagram of uncertain delay time.

2　结果

2.1　不同精度涡流预补偿方案对比

       为了说明高精度波形输出对于涡流预补偿的重要性，笔者在MATLAB中搭建预补偿计算模型对同一标准波形进行计算，通过设定不同的输出位宽来显示不同精度预补偿方案。如图6所示，图中实线表示补偿计算的实际波形值；虚线表示本方案中送入20 bit DAC的补偿波形数据；点线表示传统数字预补偿方案，即采用一片16 bit DAC输出补偿波形与标准波形之和。可以看出，本方案与实际波形吻合度较高，可以更逼真地输出补偿波形，而传统方案的输出波形幅度分辨率较低，与实际波形吻合度较差。需注意的是，两种方案的输出更新速率相同，图中相邻的纵向参考虚线之间，两种方案均输出了32次波形数据。FPGA对预补偿波形进行了精确的计算，虽然预补偿波形数据一直在改变，但由于传统方案的幅度分辨率较低，在较长的一段时间内，DAC输出的模拟波形可能仍未改变，从而造成实际输出波形与计算得到的预补偿波形之间的误差，因此达不到期望的梯度预补偿效果。从图6同样可以看出，如果DAC器件的更新时间较慢的话，即使有非常高的幅度分辨率，也不能精确输出预补偿模拟波形。本研究采用了目前国内外比较高端的DAC器件，具有较高的幅度分辨率和较快的更新时间。

点击查看大图

图6  不同精度涡流预补偿对比图
Fig. 6  Output value of eddy-current precompensation in different resolution.

2.2　涡流预补偿前后性能对比

       为了验证预补偿效果，笔者在上海康达卡勒幅医疗科技有限公司生产的OPM35I永磁型磁共振成像系统上进行了实验。使用常规圆柱形水模进行测试，用本研究中的新梯度卡替换原系统谱仪中的Z梯度卡，采用相位演化速度的测量方法^[17,18]，获得了预补偿前后的涡流时间曲线（图7)。图中实线表示预补偿前涡流场强度随时间的变化，虚线表示预补偿后涡流场强度随时间的变化，可以看到，通过对涡流进行预补偿，涡流的影响已经得到了很好的抑制。计算得到的五组涡流预补偿参数如表1所示。

点击查看大图

图7  预补偿前后涡流场对比图
Fig. 7  Gradient of eddy-current before and after precompensation.

点击查看表格

表1  涡流预补偿参数
Tab. 1  Exponential decay terms of eddy-current precompensation

2.3　成像结果对比

       为进一步验证涡流预补偿效果，本文在OPM35I系统上进行了圆柱形水模的成像实验，使用新梯度卡替代了原系统谱仪的梯度卡(Z方向）。本实验采用SE序列进行横断面成像，层厚8 mm，重复时间380 ms，回波时间21 ms，采样矩阵256×256，采样带宽25 kHz。图8A为预补偿前的水模图像，图8B为经过预补偿后得到的水模图像。可以看出，由于涡流影响，预补偿前的图像存在明显伪影，成像质量较低；预补偿后涡流引起的伪影已基本消失，成像质量得到较大提升。通过成像实验说明，本设计可以有效的降低涡流对图像的影响。

点击查看大图

图8  涡流预补偿前后水模图像对比。A：涡流预补偿前水模成像；B：涡流预补偿后水模成像
Fig. 8  Image of phantom before and after eddy-current precompensation. A: Image of phantom before eddy-current precompensation. B: Image of phantom after eddy current pre compensation.

3　讨论

       基于本实验室新设计的双DAC高精度梯度卡，笔者提出了一种新的数字涡流预补偿方法，通过FPGA实时计算后，将补偿波形与标准梯度波形分别使用两片独立的DAC进行数模转换，然后通过模拟电路实现波形叠加。本方法较大地提高了预补偿波形的幅度分辨率，采用20 bit的DAC器件获得了优于23 bit的预补偿幅度分辨率，解决了单片DAC器件分辨率不足的问题。另一方面，本研究在不引入额外接口设计的条件下，通过预补偿波形的同步触发设计和高速计算，大大缩短了补偿波形输出相对于触发信号之间的延时，同时考虑DAC器件的输入数据传输时间后，该传输延时仍只有0.83 μs，并且将该延时的不确定误差控制在25 ns以内。在梯度卡与脉冲序列发生器使用同一个相干时钟的情况下，使用本研究同步设计方法可以完全消除该传输延迟的不确定误差。

       本方法与商业磁共振成像系统OPM35I实现了兼容，在该系统上进行了水模成像实验。实验结果表明将涡流预补偿波形的精度显著提高，可以有效降低涡流，提升成像效果。因此，本技术可以提升临床实验的图像质量，减小因涡流引起的图像伪影。下一步，本技术将应用于人体成像实验，进行临床验证，并计划作为磁共振成像谱仪的一部分，应用于成像系统产品的研制中。

       本研究的不足之处在于，在进行磁共振成像实验的验证环节，对涡流的测量与补偿计算，主要考虑了涡流的时间相关性，对于涡流的空间分布仅处理了与梯度场分布相同的一阶涡流，其余涡流项未作补偿，后续将对高阶涡流项的补偿展开进一步研究，以便进一步对本研究的数字涡流预补偿方法进行验证。