2020, Vol. 41
2. 东北大学 医学影像智能计算教育部重点实验室, 辽宁 沈阳 110169;
3. 东北大学 计算机科学与工程学院, 辽宁 沈阳 110169;
4. 深圳技术大学 健康与环境工程学院, 广东 深圳 518118
2. Key Laboratory of Intelligent Computing in Medical Image, Ministry of Education, Northeastern University, Shenyang 110169, China;
3. School of Computer Science & Engineering, Northeastern University, Shenyang 110169, China;
4. College of Health Science and Environmental Engineering, Shenzhen Technology University, Shenzhen 518118, China
自旋平面回波成像(echo-planar imaging, EPI)是一种快速磁共振成像技术, 由Mansfield在1977年提出[1].EPI可在一次(单次激发EPI)或多次射频脉冲激发(多次激发EPI)后通过读出梯度场的连续正反向切换, 在回波链中产生梯度回波.由于长回波链会导致相位误差的不断累积, 因此EPI对梯度缺陷和其他因素, 如涡流、磁场不均匀等引起的相位误差非常敏感, 造成图像伪影及几何畸变.单次激发EPI成像常用于灌注加权成像、脑功能成像及扩散加权(diffusion-weighted, DW)成像.DW-EPI是一种基于EPI技术的磁共振功能成像方法, 广泛用于肿瘤、脑梗塞及中风等疾病的诊断与研究中.在DW-EPI成像过程中, 由于梯度的快速切换会在周围的磁性材料中产生涡流, 造成梯度编码错误, 产生图像伪影[2-5].另外由于涡流梯度场引起的相位误差累积到扩散梯度上, 会对信号产生调制[6-8].因此有效的涡流补偿对DW-EPI十分重要.
梯度预加重是一种常用的涡流补偿方法[9-11].为了去除涡流的影响, 得到理想的梯度场, 梯度预加重方法要根据测量到的涡流建立涡流的数学、物理模型, 计算出表征涡流特性的多个时间常数和涡流幅度值, 在梯度线圈原有梯度电流控制信号的基础上, 谱仪在输出时根据时间常数和涡流幅度值, 增加一部分预加重电流控制信号, 这部分信号产生的磁场与涡流场相反, 来抵消涡流场的影响, 从而得到理想的梯度场波形.常用的涡流测量方法有:通过分析不同空间位置球状水模的自由感应衰减信号来测量梯度涡流[12-14], 该方法对水模的摆位精度要求比较高, 做一次校正的时间比较长; 通过涡流测试序列, 采集在不同梯度延迟下的信号, 通过分析信号相位斜率测量出涡流[15].该方法速度快, 并在3T磁共振成像系统进行了验证.但在永磁磁共振成像系统上, 由于信噪比低, 信号易受噪声的影响, 使得该方法难以实现.
本文针对永磁磁共振成像系统信噪比低, 而自旋回波(spin echo, SE)序列信噪比高、抗干扰能力强的特点, 提出一种全自动梯度涡流预加重算法, 该算法校正速度快, 不需要水模精确摆位, 且不易受到噪声的影响.本方法通过测量涡流对回波中心位置(回波信号的最大模值)的影响来定量分析涡流, 为了使回波中心偏移对涡流梯度场更加敏感, 本算法使用了正、负两个涡流测试梯度, 使回波中心向两个相反的方向偏移.针对DW-EPI序列, 本文还提出一种相位补偿方法, 用于校正梯度预加重补偿后剩余涡流对扩散梯度的影响.
1 模型和算法 1.1 涡流测量设y(n)为需要产生的预加重电流控制信号, 它由理想梯度电流控制信号x(n)产生, Y(z)和X(z)分别为y(n)和x(n)的z变换, 涡流预加重系统的传输函数为H(z)=Y(z)/X(z).
由于y(n)可以用一系列具有不同时间常数、不同幅度的e指数函数的线性和来表示, x(n)为迅速变化的梯度电流, 可以看作是单位阶跃响应.这样, 涡流预加重系统的差分方程为
|
(1) |
其中:Ai为涡流幅度值; Di=e-1/Ti, Ti为涡流时间常数.要实现一个这样的涡流预加重系统, 最重要的是得到反映涡流特性的幅度值A和时间常数T.
涡流产生的梯度场累加到读梯度上, 会使SE序列的回波中心发生偏移, 基于这一原理, 涡流测试时序如图 1所示.该时序基于SE序列开发, 关闭相位编码梯度和层选梯度, 在90°射频脉冲开始之前在读梯度方向加了一个涡流测试梯度G.涡流测试梯度G到90°射频脉冲之间的时间延迟可以变化.图 1中的虚线为涡流所产生的梯度场的衰减曲线, E1为90°射频脉冲到180°射频脉冲时间段涡流梯度场对时间的积分, E2为180°射频脉冲到回波中心时间段涡流梯度场对时间的积分.
|
图 1 涡流测试时序图 Fig.1 The pulse diagram of eddy currents quantification |
从磁共振成像的基本原理可知, 当读梯度B的面积加上E2等于F的面积加上E1时, 恰好是回波中心的位置.由于E1大于E2, 要使该等式成立, 就需要增大梯度B的面积, 即回波中心会从理想位置向右移.基于相同的原理, 如果涡流测试梯度G为负梯度, 回波中心会向左移.如果没有涡流的影响, 施加正负涡流测试梯度所对应的回波中心会重合.
由于涡流可以表示为多个时间常数T和幅度A的e指数函数的和:
|
(2) |
为了推导回波中心位置偏移和涡流的关系, 以一个时间常数为例:
|
(3) |
忽略回波中心偏移, 式(3)可以近似为
|
(4) |
其中：μ为回波时间的一半; C1为常数项.
同理,
|
(5) |
正负涡流测试梯度所引起的回波中心位置差为
|
(6) |
式(6)表明, 回波中心偏差和涡流之间有近似的线性关系, 也就是说可以通过正、负涡流测试梯度引起的回波偏移来表示涡流的大小.
得到由回波中心偏移所代表的涡流衰减曲线后, 可以根据式(2)计算涡流预加重系统参数T和A.
1.2 剩余涡流补偿为了消除经梯度预加重补偿后剩余涡流对DW-EPI扩散梯度的影响, 本文对180°射频脉冲后的扩散梯度进行了相位误差补偿, 如图 2所示, D表示扩散梯度.
|
图 2 剩余涡流引起的相位误差补偿 Fig.2 Phase error compensation due to residual eddy currents |
图 2中, t0为90°射频脉冲开始时间, t1为第一个扩散梯度的开始时间, t2为第一个扩散梯度的结束时间, t3为180°射频脉冲开始时间, t4为第二个扩散梯度的开始时间, t5为第二个扩散梯度的结束时间, t6为产生的回波中心时间.
在产生回波之前, 由剩余涡流累积的相位应该为0, 因此对第二个扩散梯度补偿的相位为
|
(7) |
其中,
|
(8) |
|
(9) |
|
(10) |
|
(11) |
采用本文所提出的方法, 在0.35T开放式永磁磁共振成像系统中进行了实验.该系统支持5个梯度预加重单元, 实验采用4通道膝关节接收线圈及细棍状硫酸铜模体.模体摆放在磁体中心, 长轴方向与所校正的涡流方向一致.采集到的数据使用交互式数据语言IDL进行数据分析及指数拟合.
涡流测试序列的参数:重复时间为1000ms, 回波时间为40ms, 视野为1000mm, 采样点数为256.涡流测试梯度的强度为±15mT/m.在涡流测试梯度结束到90°射频脉冲开始之前, 采用20个可变的延迟时间(范围为0~500ms), 分别加正、负涡流测试梯度进行回波采集.为了提高计算精度, 消除涡流对回波形状的影响, 本文采用回波的质心来代替回波中心.最后采用最小二乘法对采集到的0~500ms内的涡流根据式(2)进行拟合.由于最小二乘法对初值敏感, 因此在拟合时采用了分段拟合的办法来确定拟合初始值.
DW-EPI序列采用单激发自旋回波EPI序列, 时序如图 3所示.在180°射频脉冲两侧分别加上扩散梯度(虚线), 扩散梯度加在频率编码Gx、相位编码Gy与选层Gz方向上.
|
图 3 DW-EPI脉冲时序图 Fig.3 Pulse diagram of DW-EPI |
为了定量分析涡流校正对图像的几何畸变及伪影的影响, 本文以国家食品药品监督管理局发布的“YY/T 0482—2010医用成像磁共振设备主要图像质量参数的测定”为依据, 采用如下方法进行定量分析.
对于伪影, 采用3个测量值:1)在伪影处画感兴趣区域(region of interest, ROI)，找到最大伪影值Gmax; 2)在测量模体内画ROI, 计算信号平均值Savg; 3)在背景处画ROI, 计算背景噪声平均值Navg.伪影与信号的比值g为
|
(12) |
设模体的实际直径为d, 对于几何畸变, 计算模体在水平、垂直以及左右倾斜45°角的直径的平均值
|
(13) |
图 4为Gx, Gy和Gz三个方向的涡流补偿情况.经过涡流补偿后, Gx, Gy和Gz三个方向的涡流均降低至1个像素以内.表 1为涡流补偿前后最大涡流值的变化情况.从表 1可以看出, 经过涡流补偿后各方向的涡流值均有不同程度的下降, Gx方向的最大涡流值约降低1/15, Gy方向的最大涡流值约降低1/18, 涡流最大的Gz方向约降低了1/107.
|
图 4 涡流补偿结果 Fig.4 The eddy currents compensation results (a)—Gx方向; (b)—Gy方向; (c)—Gz方向. |
|
表 1 涡流补偿前后的最大涡流值 Table 1 Maximum value of eddy currents with or without compensation |
为了提高涡流测量的准确性、灵敏性, 消除其他梯度涡流对测量的影响，应注意以下几方面:1)图 1中的涡流测试梯度G的持续时间要足够长, 让梯度上升沿产生的涡流有充分的时间衰减, 不会对梯度下降沿的涡流产生影响.在实验时, 梯度G的持续时间设置为2s.2)对图 1中的散相梯度F, 应提高其梯度的强度, 缩短梯度的持续时间, 这样梯度上升沿产生的涡流会和下降沿产生的涡流相互抵消.由于本算法采用了正负涡流测试梯度相减的方法, 对未抵消的涡流也会被减掉, 因而散相梯度F不会对涡流测量产生影响.3)对图 1中的读梯度B, 应尽可能降低梯度的强度, 增加梯度的持续时间, 这样很小的梯度变化就会产生很大的回波中心偏移.
由于1ms以内的涡流变化复杂, 在1ms以内设置过多的时间延迟(采样点)会降低自动拟合算法的稳定性, 而且本系统支持的时间常数有限, 因此虽然经过涡流补偿, 系统三个方向的涡流均有下降, 但短时间尤其是1ms以内仍然有剩余涡流存在.如Gx方向的第一个采样点, 经过涡流补偿后仍有约6%的涡流存在.因此DW-EPI序列对剩余涡流引起的相位误差进行补偿是必要的.
图 5a, 图 5c为关闭涡流补偿后的扩散梯度加权图像, 图 5b, 图 5d为经过涡流补偿后的扩散梯度加权图像.依据前文所述的定量分析方法, 在图 5a, 图 5b上分别统计出背景噪声强度的平均值、最大伪影强度值, 以及信号强度的平均值.在图 5c，图 5d上, 圆心以标准SE图像为参照, 按照采集到的模体外轮廓画4条直径, 这些直径同心, 直线间夹角为45°, 根据直径的平均值来计算几何畸变.
|
图 5 涡流引起的伪影及畸变 Fig.5 Artifacts and geometrical distortion caused by eddy currents (a)—无涡流补偿时伪影测量信息; (b)—有涡流补偿时伪影测量信息; (c)—无涡流补偿时几何畸变测量信息; (d)—有涡流补偿时几何畸变测量信息. |
涡流引起的几何畸变及伪影定量值如表 2所示.从表中可以看出, 打开涡流补偿后伪影会降低1/3, 几何畸变会降低5%.
|
|
表 2 涡流补偿前后的伪影和几何畸变定量对比 Table 2 Quantitative comparison of artifacts and geometric distortion with or without eddy currents compensation |
没有经过涡流补偿的模体扩散加权图像, 由于涡流的影响, 使不同扩散梯度方向的图像产生了不同方向的位移、错切和拉伸, 如图 6所示.图 6a, 图 6b不但产生了错切变形, 图像中心也发生了偏移, 图 6c有拉伸变形.由于扩散梯度在各方向产生的图像形变不一样, 这样在3个扩散梯度图像合成时, 由于图像无法对齐, 因而会产生图像畸变和伪影.
|
图 6 DW-EPI不同方向扩散加权图 Fig.6 DW-EPI diffusion-weighted images of three directions (a)—扩散梯度加在Gx方向的模体图像; (b)—扩散梯度加在Gy方向的模体图像; (c)—扩散梯度加在Gz方向的模体图像; (d)—SE序列的模体图像. |
1) 本文针对永磁磁共振成像系统信噪比低、涡流校正困难这一问题, 提出了一种涡流预加重补偿算法.由于所提出的算法基于SE序列开发, 因此该算法具有SE序列的优点:信噪比高，算法的稳定性更好.
2) 对于剩余涡流场对DW-EPI扩散梯度累积的相位误差, 本文提出了一种方法进行相位补偿.
3) 从DW-EPI的模体图像可以看出, 本文提出的方法能有效去除图像几何畸变, 降低图像伪影.
| [1] |
Mansfield P. Multi-planar image formation using NMR spin echoes[J]. Journal of Physics C:Solid State Physics, 1977, 10(3): L55-L58. DOI:10.1088/0022-3719/10/3/004 |
| [2] |
Bennett L H, Wang P S, Donahue M J. Artifacts in magnetic resonance imaging from metals[J]. Journal of Applied Physics, 1996, 79(8): 4712-4714. DOI:10.1063/1.361649 |
| [3] |
Horsfield M A. Mapping eddy current induced fields for the correction of diffusion-weighted echo planar images[J]. Magnetic Resonance Imaging, 1999, 17(9): 1335-1345. DOI:10.1016/S0730-725X(99)00077-6 |
| [4] |
Aliotta E, Moulin K, Ennis D B. Eddy current-nulled convex optimized diffusion encoding(EN-CODE) for distortion-free diffusion tensor imaging with short echo times[J]. Magnetic Resonance in Medicine, 2018, 79(2): 663-672. DOI:10.1002/mrm.26709 |
| [5] |
Shen Y, Larkman D J, Counsell S, et al. Correction of high-order eddy current induced geometric distortion in diffusion-weighted echo-planar images[J]. Magnetic Resonance in Medicine, 2004, 52(5): 1184-1189. DOI:10.1002/mrm.20267 |
| [6] |
Jeong E K, Kim S E, Parker D L. High-resolution diffusion-weighted 3D MRI, using diffusion-weighted driven-equilibrium(DW-DE) and multishot segmented 3D-SSFP without navigator echoes[J]. Magnetic Resonance in Medicine, 2003, 50(4): 821-829. DOI:10.1002/mrm.10593 |
| [7] |
Mueller L, Wetscherek A, Kuder T A, et al. Eddy current compensated double diffusion encoded(DDE) MRI[J]. Magnetic Resonance in Medicine, 2015, 77(1): 328-335. |
| [8] |
Papadakis N G, Martin K M, Pickard J D, et al. Gradient preemphasis calibration in diffusion-weighted echo-planar imaging[J]. Magnetic Resonance in Medicine, 2000, 44(4): 616-624. DOI:10.1002/1522-2594(200010)44:4<616::AID-MRM16>3.0.CO;2-T |
| [9] |
Ahn C B, Cho Z H. Analysis of the eddy-current induced artifacts and the temporal compensation in nuclear magnetic resonance imaging[J]. IEEE Transactions on Medical Imaging, 1991, 10(1): 47-52. |
| [10] |
Tao S, Weavers P T, Trzasko J D, et al. Gradient pre-emphasis to counteract first-order concomitant fields on asymmetric MRI gradient systems[J]. Magnetic Resonance in Medicine, 2016, 77(6): 2250-2262. |
| [11] |
Stich M, Wech T, Slawig A, et al. Gradient waveform pre-emphasis based on the gradient system transfer function[J]. Magnetic Resonance in Medicine, 2018, 80(4): 1521-1532. DOI:10.1002/mrm.27147 |
| [12] |
Liu Q, Hughes D G, Allen P S. Quantitative characterization of the eddy current fields in a 40 cm bore superconducting magnet[J]. Magnetic Resonance in Medicine, 1994, 31(1): 73-76. DOI:10.1002/mrm.1910310112 |
| [13] |
Chinchali A P.Eddy current correction in PC-MRI: an analysis of local and global static tissue fitting techniques[D]. Los Angeles: University of California, 2018.
|
| [14] |
Robertson S, Hughes D G, Liu Q, et al. Analysis of the temporal and spatial dependence of the eddy current fields in a 40cm bore magnet[J]. Magnetic Resonance in Medicine, 1992, 25(1): 158-166. DOI:10.1002/mrm.1910250116 |
| [15] |
Schmithorst V J, Dardzinski B J. Automatic gradient preemphasis adjustment:a 15-minute journey to improved diffusion-weighted echo-planar imaging[J]. Magnetic Resonance in Medicine, 2002, 47(1): 208-212. DOI:10.1002/mrm.10022 |