国内低场磁共振成像系统大多采用永磁磁体，相较于高场超导系统，永磁系统场强低（场强一般为0.5 T及以下），并且磁场均匀性和稳定性都比较差，磁场强度受温度的影响也比较大。主磁场非均匀性对成像的影响程度一直是研究人员、技术人员与使用者关心的问题，以往主要靠定性分析，缺乏直观的图像表达。

       自1984年Bittoun等^[1]建立第一个基于Bloch方程的MRI模拟器以来，迄今国内外对其进行了大量广泛深入的研究，主要集中在以下几个方面：模拟器从1D到2D、3D推广及应用^[2,3]；组织解剖模型的建立^[4,5,6]，根据需求可分为离散组织模型和模糊组织模型；为了加快计算速度，模拟器算法的优化及并行化^[7,8,9]。这些研究取得了良好的效果，使其可以应用于低成本大规模磁共振成像教学、磁共振实验设计、MRI序列编写、参数优化、系统伪影等问题的解决^{[10,11,12,13]}。不过，对于低场中主磁场不均匀性的模拟一直没有系统的研究与模拟。

       本文实现一种基于Bloch方程的MRI模拟器，通过对主磁场不均匀性分布进行建模，提出了对主磁场不均性的模拟仿真方案。以自旋回波（spin echo ，SE）和梯度回波（gradient echo，GRE）序列成像方法为例对主磁场的不均匀性进行了模拟，并得到了预期的实验结果。

1　材料和方法

1.1　模拟器

       模拟器构成框图如图1所示，包括虚拟物体、静磁场定义、MRI序列、磁化强度矢量计算核心，计算出射频（radio frequency，RF）信号（k空间信号），加入噪声，最后用快速傅立叶变换（fast Fourier transform，FFT）重建图像。

       虚拟物体是对真实自旋系统的离散化描述。虚拟物体的每一个体素包含一组用于通过Bloch方程计算其自旋磁化矢量的物理量：质子密度(ρ)、两个弛豫常量T1和T2。

       MRI序列主要包含3个事件：自由进动、RF脉冲激发、信号采集。在MRI模拟器中通过排列组合以上3种事件、设置事件参数得到不同MRI序列。MRI序列实施过程中的磁化强度矢量计算是基于简化的3D Bloch方程的解，Bloch方程描述了自旋磁化矢量随时间的变化，如式（1）所示：

       其中M0是自旋磁化矢量的平衡常数，由质子密度决定；(T1, T2)是弛豫时间；γ是物体成分的旋磁比。磁场强度由式（2）表示：

       式中B⁰是主磁场强度，是主磁场不均匀性，表示施加的梯度场，表示RF脉冲，表示空间坐标。

       模拟计算核心通过旋转矩阵和指数扩展执行Bloch方程的离散时间解，得到演变的磁化强度矢量。对于RF信号采集，一次激发可以得到一维离散复数信号，可以填充k空间的一条相位编码线。RF信号的一点S [t]是对应时间点虚拟物体所有体素的磁化强度矢量之和，如式(3)所示：

点击查看大图

图1  MRI模拟器构成框图
Fig. 1  The overview frame of the MRI simulator

点击查看大图

点击查看大图

点击查看大图

1.2　系统参数设置

       创建一个2D虚拟水模，对模拟器的2D模拟成像效果进行实验验证。2D虚拟水模模型：FOV 250 mm × 250 mm，采样矩阵和重建矩阵同为256 × 256；水模大小150 mm×150 mm，弛豫时间T1=100 ms，T2=100 ms（参考实验用水模：1L水中加入3.6 g NaCl和1.995 g CuSO4·5H2O）；主磁场强度B0为0.4 T；在k空间加入信噪比为40 db的高斯白噪声。分别用SE序列和GRE序列进行实验，SE序列参数为：TR 300 ms，TE 20 ms，采样时间间隔为48 us；GRE序列参数为：TR 80 ms，TE 11.6 ms，射频脉冲翻转角为70°，采样时间间隔为36 us。虚拟成像中不考虑涡流的影响，梯度的上升下降时间设置为零。在实际成像实验中除将梯度的上升下降时间均设为532 us外，其他参数设置与上述参数相同。

1.3　主磁场不均匀性建模

       定义为主磁场B0的不均匀性，是空间位置坐标。主磁场不均匀会影响质子的自旋频率，造成一定的频率偏移（），并加速散相，加快信号的衰减。在本模拟器中，主磁场本身的不均匀性可以通过一个空间位置矩阵给出，通过建模给出主磁场的变化与空间位置的关系，给矩阵赋值。

       用∆B0表示系统的匀场精度，(x,y)表示体素点的坐标位置，图像大小为N×N。永磁系统中主磁场的非均匀分布可以用基于球坐标的解析法分解为一系列线性无关的空间谐波的线性组合^[14]，在直角坐标系中可以展开为一个关于空间坐标的多项式。据此假定主磁场不均匀性分两种模式：线性分布、抛物线形式分布。

       模式一：线性分布函数如式（4）所示，当场不均匀性沿x轴线性变化时，i代表体素的x坐标；当场不均匀性沿y轴线性变化时，i代表体素的y坐标。

       模式二：抛物线形式的分布函数如式(5)所示：

点击查看大图

点击查看大图

2　结果

2.1　主磁场均匀条件下的模拟成像

       设置∆B0=0 ppm，图2给出了实际实验和虚拟MRI的SE和GRE序列成像结果。在实际成像中，因为圆柱形水模中存在气泡，实验结果中存在缺口。

       实际实验中因为射频场的不均匀、梯度的非线性、接收线圈的敏感度变化等因素得到的磁共振图像信号强度并不是均匀的。相比而言，因为虚拟成像中设定的都是理想化的实验条件，所以图像各像素点的信号强度基本相同。

点击查看大图

图2  SE和GRE序列的MRI图像。A：虚拟成像结果；B：实验结果
Fig. 2  MRI images obtained with SE and GRE sequences. A: the simulated imaging results; B: the experimental imaging results.

2.2　主磁场不均匀条件下的模拟成像

       图3分别给出了模式一下SE和GRE序列的虚拟图像，其中A、C表示场不均匀性沿x轴（频率编码方向）线性变化时的图像，B、D表示场不均匀性沿y轴（相位编码方向）线性变化时的图像。从图3可以看出，随∆B0增大图像几何失真逐渐增强。当∆B0增加到50 ppm，图像几何失真会非常明显。现在常用的永磁系统匀场精度能够达到30 ppm，对图像的影响一般视觉上是不可见的。脉冲序列成像对频率编码方向的场不均匀性较敏感，对相位编码方向的场不均匀性不敏感。

       实际实验中，分别将x方向的梯度值设定为正常梯度的99%、y方向梯度不变，y方向的梯度值设定为正常梯度的101%、x方向梯度不变，将SE序列和GRE序列的成像与相应的虚拟成像比较，结果如图4所示。它们的几何失真幅度和模式几乎相同，可见模拟器对线性主磁场不均匀具有良好的模拟效果。

       在场不均匀度相同的情况下，图像的畸变会随采样带宽增大而减小。图5给出了在模式一x方向上的匀场精度为125 ppm情况下，采样带宽分别为20.83 kHz和41.67 kHz时SE序列的虚拟成像结果。在图像中，频率编码方向水模大小从20.83 kHz采样时的129像素变为41.67 kHz采样时的121像素，图像畸变减小。

       图6给出了模式二下的SE和GRE序列虚拟成像。图像的信号强度发生变化，并且GRE序列中信号变化比SE序列明显。由于180°重聚相位脉冲存在，SE脉冲序列成像对场不均匀性引起的信号衰减不敏感。

点击查看大图

图3  模式一下的虚拟成像。A、B为SE序列图像，C、D为GRE序列图像；A、C为不均匀性沿x轴线性变化，B、D为不均匀性沿y轴线性变化；deltaB0表示匀场精度
Fig. 3  The virtual MRI with various shimming accuracy (deltaB0) of model one. A, B were obtained by SE sequence; C, D were obtained by GRE sequence; The main field inhomogeneity linearly related to the corresponding axis, A and C were x axis, B and D were y axis.

点击查看大图

图4  虚拟实验（A）与实际成像实验（B）对比。(1)、(2)为SE序列成像，(3)、(4)为GRE序列成像；(1)、(3)为改变x方向的梯度，(2)、(4)为改变y方向的梯度
Fig. 4  Comparison of (A) experimental and (B) simulated images. (1) and (2) were SE sequence images, (3)and (4) were GRE sequence images; (1) and (3) were images with shimming of x_axis, (2) and (4) were images with shimming of y_axis.

点击查看大图

图5  采样带宽对MRI图像的影响。采样带宽：A：20.83 kHz，B：40.67 kHz
图6  模式二下的虚拟成像。A：SE序列图像，B：GRE序列图像，deltaB0表示匀场精度
Fig. 5  Impact of the bandwidth to MRI. the bandwidth of A is 20.83 kHz, the bandwidth of B is 40.67 kHz.
Fig. 6  The virtual MRI with various shimming accuracy (deltaB0) of model two. A was obtained by SE sequence, B was obtained by GRE sequence.

3　讨论

       通过与实际实验成像的对比可得本模拟器对线性主磁场不均匀性具有较好的模拟效果。不过在模拟中，只考虑了在空间上的场不均匀性，而没有提到随时间变化的场不均匀性。考虑磁场随时间的变化需要为磁场分布增加一个时间维度，增加了建模的复杂度，操作比较困难，但这可以是以后的一个研究方向，以模拟实际磁共振成像中主磁场随温度等周围环境变化的情况。

       实际MRI过程中，图像还会受到RF的不均匀性、梯度的非线性及采集线圈的空间敏感度等因素的影响，因为对它们的直接模拟不容易实现，本文没有涉及。为了更真实准确地模拟MRI成像过程，可以考虑将这些因素通过一些应用的选择实现允许用户加入到模拟器中。

       基于简化的3D Bloch方程的离散时间解来计算k空间磁化矢量的方法相当消耗时间。由式（3）可见，对于包含体素数为N²、k空间栅格大小同样为N²的情况，算法的复杂度为0 (N⁴)。由此可见随着虚拟物体包含的体素增加，算法的计算量增加会越来越快。对于上面的2D SE序列成像实验，256 × 256虚拟物体k空间磁化矢量的计算时间大概为50 min，所用的处理器为Intel Core i3，主频为2.4 GHz。同时由式（3）还可以看出，磁化矢量计算核心具有很好的可并行性，利用GPU、Open-CL等方法实现程序并行化可以大幅度的减少程序运行时间。这也是模拟器的主流研究方向之一。

       本文提出的MRI模拟器是一种基于MRI事件的模拟器，它可以通过事件的排列仿真多种MRI序列；可并行化的磁化矢量计算核心算法，为程序的并行化发展提供了广阔的可能性；通过对主磁场不均匀性进行空间建模，并以此为带入条件进行模拟成像，提高了模拟器的对场不均匀性模拟的准确性；可以通过设置与低场磁共振成像设备相同的参数进行对实际情况的模拟仿真，有助于工程与科研中分析主磁场不均匀对成像的影响，帮助研究人员更好地进行主磁场矫正。