涉及磁共振（MR）成像模拟的方法与流程

本发明涉及一种涉及磁共振(mr)成像模拟的方法。

背景技术：

0、发明背景

1、本发明总体上涉及不需要使用真实mri扫描仪的高级mri模拟。更具体地，本发明涉及一种用于磁共振(mr)成像模拟的方法，其中使用优选的计算机模型，使得可以有利的方式分析对象的特定理论体积。

技术实现思路

1、上述后一目的是通过涉及磁共振(mr)成像模拟的方法实现的，所述方法包括：

2、-将脉冲序列划分为具有rf激励的部分和不具有rf激励的部分；-对于被称为梯度部分的不具有有源rf的部分中的每个部分，对该梯度部分执行压缩，并且然后将驱动该梯度的信号表示为读出时间点之前的面积累积；以及

3、-然后执行模拟。

4、本发明涉及一种方法，其中聚焦于不具有rf激励的部分。这进一步意味着该方法涉及分离出这些非rf部分，这在关闭rf发射时实现。

5、此外，从上文应当理解，根据本发明的方法包括对梯度部分进行压缩，并且然后将驱动梯度的信号表示为读出时间点之前的面积累积，这意味着对序列的时间压缩被定义为在读出时间点之前的面积累积(参见以下(a，h)对)而不是时间信号。与现今已知和使用的方法相比，这是一个新颖的特征和概念。

6、对本发明和相关技术背景的描述

7、用于快速mri模拟的基于面积的序列压缩

8、mri模拟的基础

9、mri技术基于在存在强恒定磁场(主磁场)的情况下以及当被较小的变化磁场激发时，对软组织中某些原子质子(自旋)的磁化的时间演化的视觉解释。

10、因此，mri基于自旋动力学的隐式使用，并且在存在外部变化磁场的情况下的磁化的演变通过布洛赫(bloch)方程[1]在宏观尺度上建模。对mri现象的模拟依赖于求解布洛赫方程(或包含其它物理特性的升级的衍生版本，例如，用于扩散的布洛赫-托里(bloch-torrey)方程[2])。

11、用于mri目的的布洛赫方程的基本形式是：

12、dmx(r，t)/dt＝γ(bz(r，t)my(r，t)-by(r，t)mz(r，t))-mx(r，t)/t2(r)

13、dmy(r，t)/dt＝γ(bx(r，t)mz(r，t)-bz(r，t)mx(r，t))-my(r，t)/t2(r)

14、dmz(r，t)/dt＝γ(by(r，t)mx(r，t)-bx(r，t)my(r，t))-(mz(r，t)-m0z(r))/t1(r)

15、其中：

16、r是″自旋″或质子的坐标；

17、t是时间；

18、γ是旋磁比；

19、mx(r，t)、my(r，t)、mz(r，t)是磁化m(r，t)的笛卡尔分量；

20、dm_(r，t)/dt是磁化的给定分量的时间导数；

21、bx(r，t)、by(r，t)、bz(r，t)是外部磁场b(r，t)的笛卡尔分量；

22、t2(r)是组织相关性自旋-自旋弛豫时间；

23、t1(r)是组织相关性自旋-晶格驰豫时间；并且

24、m0z(r)是初始(或平衡)磁化。

25、mri的基础

26、mri基于一系列事件。

27、平衡：通常在z方向上施加恒定的强主磁场b0，其将自旋对准到初始状态m0z(平衡)，并且根据旋磁比产生进动。关于z，在上文和下文中，z方向被认为是主磁场的方向。

28、rf激励：施加具有进动频率下的横向(x和y)磁分量(bx和by)的短射频(rf)场，以便扰乱自旋的磁化，并且将其带到特定的mx、my和mz组合。

29、编码和采集：在关闭rf之后，施加z分量中的梯度场(随位置变化的bz)以生成所谓的″空间编码″，使得每个体素具有与其位置相关的进动，并且因此存在mx和my磁化(横向磁化)的差异，该差异随后被用于创建图像。

30、这些横向磁化的组合效应辐射出具有x和y分量的rf场，该rf场使用法拉第定律被线圈接收，并且被变换成电信号s(t)＝s(mx(t)，my(t))，其中s是接收线圈灵敏度的函数。

31、弛豫：(可以是全部或部分)没有施加外场，使得横向(xy)磁化将由于t2而自然衰减，而mz磁化将在t1恒定的情况下朝向m0z恢复。

32、mri使用这些基本事件的精确定时和重复序列来获取时间信号s(t)。使用编码的知识将这些信号解析成图像，并且所应用的信号的类型将使信号或多或少地取决于t1和t2，从而导致基于组织特性的自然对比度。

33、mri模拟的具体细节

34、假定mri在z方向上的恒定强静态场b0下操作，bz(t)＝b0+δbz(t)(其中b0是恒定场值，并且δbz(t)是由于外部激励或不均匀性和非线性引起的变化)，则布洛赫方程的模拟的典型表示基于旋转参考系[4]，其中从bz去除b0的进动效应以简化公式。

35、bx和by场是由于rf激励而产生的，而bz场是由于梯度场的效应和其它效应而产生的，诸如由于物理效应(易感性、化学位移)或硬件缺陷而产生的不均匀性，使得一旦移除b0场的效应，磁场的z分量即为：bz(r，t)＝δbz(r)+gx(t)gx(r)+gy(t)gy(r)+gz(t)gz(r)

36、其中

37、δbz(r)是非理想因素；

38、gx(r)gy(r)gz(r)是由梯度线圈产生的场图；并且

39、gx(t)gy(t)gz(t)是驱动梯度线圈的mri序列的时间信号。

40、在这一点上，还应当注意，在更高级的情况下，梯度场图也可以是时变的：gx(r，t)gy(r，t)gz(r，t)。

41、布洛赫方程是一种常微分方程(ode)[3]，并且可被表示为dm/dt＝f(m，b，t1，t2，m0z)，其中f是函数，并且时间和位置相关性是隐含的。

42、在给定初始解的情况下，可通过时间离散化对ode进行数值求解，其中时间被分成小段，在这些小段期间假设b场是恒定的，并且通过对先前解进行时间积分，可迭代地计算每个时间步长处的磁化的解。

43、存在可应用的不同的众所周知的方法[3]。例如，使用显式前向差分方案，其中dt＝t(n)-t(n-1)，并且已知t(n-1)处的m(t)，可求解差分方程(m(t+dt)-m(t))/dt＝f(m(t)，b，t1，t2，m0z)来计算m(t+dt)。

44、另选地，通过应用矩阵指数，可获得恒定场时间步长的分析解。给定布洛赫方程的结构，该操作可被分解为一系列旋转(取决于b场)和指数衰减[4，5]，以压缩形式表示为：

45、m(t+dt)＝e(dt)r(dt)m(t)+eeq(dt)

46、其中对于每个时间步长：

47、m(t)是磁化矢量，m(t)＝[mx(t)，my(t)，mz(t)]′，其中[]′表示转置，e(dt)是具有取决于t1和t2的指数衰减的3x3矩阵，e(dt)＝diag(exp(-dt/t2)，exp(-dt/t2)，exp(-dt/t1))

48、eeq(dt)是纵向恢复矢量，eeq(dt)＝[0，0，m0z*(1-exp(-dt/t1))]’，其中m0z是平衡磁化。

49、并且r(dt)是3x3旋转矩阵，取决于b场。

50、在笛卡尔分量中，在间隔dt期间指数衰减操作e(dt)和eeq(dt)的效应是mx(t+dt)＝exp(-dt/t2)*mx(t)

51、mx(t+dt)＝exp(-dt/t2)*my(t)

52、mz(t+dt)＝m0z+(mz-m0z)*exp(-dt/t1)

53、其中t1和t2是组织相关性弛豫常数。

54、r(dt)中的操作包括一系列旋转，该旋转的角度取决于间隔dt期间的b场。典型的表示[4]将r(dt)分解为一系列旋转

55、

56、其中每个旋转可以矩阵形式表示

57、

58、

59、相同操作的其它表示和近似是可能的[6]，这在计算方面对于某些情形可能是有利的。

60、在任何情况下，数值布洛赫模拟依赖于时间离散化为具有恒定b场的小间隔，以及基于间隔开始处的磁化的数值或分析时间积分。

61、需注意，这是一个顺序过程，其中新的解需要计算先前的解。

62、对于变化的场，近似的精度取决于时间步长，其中误差随着间隔的持续时间变小而变小。但是由于需要将每个步长应用于可能数百万个体素，所以精细离散化会导致昂贵的模拟。在这些场长时间保持恒定的情况下，可用单个(较大)时间步长获得解而没有任何精度损失[4]。

63、需注意，对于大多数mri模拟，需要磁化的唯一相关时间点是在采集步骤(被称为读出)期间，或者在序列的开始和结束处。这些是在实践中为了精确模拟而需要计算的时间点的一小部分。

64、出于其中想要模拟含有生物组织的3d视场(fov)的实际目的，还应用空间离散化，其中3d fov被划分为小体素。体素是包含具有相同组织特性的多个核的小体积，并且由于它们的接近性，假设经历相同的外场。因此，假设b、m、t1、t2和m0z在被分配了坐标r和体积的每个体素处是恒定的。该体素的磁化可利用以上过程来模拟。关于上文，应当注意，本发明还涵盖在单个体素内实际上有多于一种组织类型的情况。

65、用于mri的典型布洛赫模拟基于以上提供的数值方法，其中时间在对时域的小离散化中被离散为小步长。

66、本发明的具体实施方案和相关细节

67、本发明涉及一种包括基于面积的序列压缩的方法。根据本发明，本发明人提出应用该序列的替代表示。作为预处理步骤，将序列划分成多个部分，分离成具有rf激励并因此具有时变bx、by和bz场的部分和不具有rf激励的部分(序列的编码、采集和弛豫部分)，其中仅存在时变bz场，即bx(t)＝by(t)＝0。这可以通过对具有驱动发射线圈的非零信号的时间步长进行识别和分组来完成。在尚未使用时间离散化来表示序列的情况下，可预先应用时间离散化。

68、对于不具有有源rf的部分(被适当地称为梯度部分)中的每个部分，驱动梯度的信号可被表示为读出时间点之前的面积a的累积：a＝f(幅度，h，t)，其中f是函数，并且h＝t(读出)-t0，并且到该部分的初始时间且仅存储感兴趣时间(读出)的这些值对(a，h)。面积a可以通过时间信号的积分来计算(使用黎曼积分以分析方式或以数值方式进行)。

69、基于上文，根据本发明的一个具体实施方案(以下被称为第一实施方案方面)，对该梯度部分的压缩通过预先计算(a，h)对来执行，或者通过累积相位和时间来即时执行，其中a是驱动梯度的作为读出时间点t(读出)之前的累积面积的信号，并且其中h是h＝t(读出)-t0，其中t是时间。

70、与上文一致，根据本发明的另一个具体实施方案，驱动梯度的作为读出时间点之前的累积面积的信号a被呈现为信号幅度随时间的积分，优选地诸如通过a＝sum_i(dt*amplitude(t_i))在数值上近似，其中t_i＝t(i)，dt＝t(i)-t(i-1)，i处于离散区间[1，读出]中。

71、一旦已经将每个梯度部分处的序列表示为一组面积(a)和周期(h)，则可以执行模拟。需注意，γ*a*梯度的乘积提供了由于对应信号而累积的相位，并且因此由于bz场而产生的旋转对应于：

72、

73、其中r是″自旋″或质子的坐标，δbz(r)是位置相关的非线性场效应，gx(r)、gy(r)和gz(r)是由于梯度线圈而产生的磁场(bz)的z分量的空间分布，并且是时间周期h期间的位置相关的旋转角，其中ax、ay和az是通过梯度信号gx(t)、gy(t)和gz(t)中的每个梯度信号的时间积分计算的面积。

74、而且，在这种情况下，能够计算每个读出点处的磁化的解，给定梯度部分的开始(t0)处的初始磁化(m0)，其中由旋转和指数衰减组成的单个步长应用于梯度序列部分的开始处的初始磁化m0。

75、基于上文，根据本发明的一个具体实施方案，该方法通过执行以下步骤来执行：

76、对于在时间t处的每次读出，其中h＝t-t0：

77、计算衰减和旋转：

78、e(h)、eeq(h)

79、将它们应用于m(t0)以获得在读出时间t处的解：

80、

81、其中e(h)和eeq(h)是在时间周期h期间的弛豫常数t1和t2的效应，并且是由于相位的累积而产生的旋转，与bz场直接相关。

82、在这种情况下，由于仅具有场的z分量，所以旋转是围绕z轴的旋转：

83、

84、因此，可以在笛卡尔分量中表示为

85、

86、mz(t)＝m0z+(mz(t0)-m0z)*exp(-h/t1)

87、其中mx(t0)、my(t0)和mz(t0)是梯度部分的开始(t0)处的初始条件。

88、需注意，在没有压缩的情况下，需要替代地以迭代方式计算解：

89、m(t0+dt)＝e(dt)r(dt)m(t0)

90、m(t0+2dt)＝e(dt)r(dt)m(t0+dt)

91、...

92、m(t0+h)＝e(dt)r(dt)m(t0+h-dt)

93、而且，需注意，迭代过程等效于将信号划分为持续时间dt的分段恒定步长，并且应用与该小步长相关联的区域的旋转。

94、r(dt)＝r(rho)

95、其中

96、rho(r)＝γ*dt*(δbz(r)+gx(t＝ti)*gx(r)+gy(t＝ti)*gy(r)+gz(t＝ti)*gz(r))

97、在某种意义上，传统(迭代)的方法在概念上类似于使用每个步长的单独旋转和衰减来应用黎曼积分。

98、此外，下文提供了根据上文公开的本发明的第一方面的实施方案的概念的一组含义和优势：

99、1)需要预先计算以针对每个感兴趣点(该部分的读出和结束点)，将序列gx(t)、gy(t)和gz(t)重新表示为具有面积和持续时间的(a，h)对的集合。这非常高效且便宜，因为它只需要处理序列，并且每个序列只需要进行一次。

100、2)能够任意地控制面积的计算精度。对于分析积分是精确的；-对于数值积分，可细化对时间信号的离散化，而不会影响实际模拟的性能。

101、3)代替需要对具有dt离散化的大量时间步长应用数值模拟，只需要对感兴趣点(读出)中的每个感兴趣点应用一次旋转和衰减。这对于实际波形(梯形)或更高级的编码(例如采用正弦波形的螺旋)尤其相关，并且为了高精度将需要精细的离散化。这导致性能大幅提升，特别是因为实际模拟可能必须将这些操作应用于数百万个体素上。

102、4)每次读出时的解仅取决于梯度部分的开始处的初始解m(t0)及其(a，h)。因此，可以采用多个计算机机器或计算机内核来并行地计算每次读出的解。这进一步提升了潜在的性能改进。基于此，根据本发明的一个具体实施方案，一旦已经预先计算了(a，h)对，就并行执行针对每个读出时间点的模拟。

103、5)仍然需要以常规方式模拟rf部分。并且这些部分的模拟需要按顺序：每个梯度部分的初始m(t0)仅在先前部分已被模拟时可用。

104、6)磁场不均匀性或其它效应的影响可被解释为恒定(或时变，如果适用)信号，并且可如上文所示针对每次读出被包括为γ*h*δbz(r)

105、根据本发明的第二实施方案方面，存在磁化的替代表示，并且因此存在所应用的布洛赫方程的替代表示。代替以笛卡尔坐标表示系统，提出使用极坐标，将横向磁化表示为时间相关的相量。因此，根据本发明的一个实施方案，所述方法还包括将横向磁化表示为时间相关的相量，并且因此以极坐标表示系统。

106、对于幅度(mm)和相位可应用下式：

107、

108、其中mxy是横向磁化，并且i是复合单位。

109、在笛卡尔分量和极坐标分量之间存在明确的关系：

110、

111、

112、mz＝mz。

113、使用这种表示具有一些优势。与上文公开的第一实施方案方面的情况类似，当正在模拟不具有有源rf的时间点(bx(t)＝by(t)＝0)时，用于步长dt的模拟的计算简化为

114、

115、mm(t+dt)＝mm(t)*exp(-dt/t2(r))

116、mz(t+dt)＝m0z+(mz(t)-m0z)*exp(-dt/t1(r))。

117、基于上文，根据本发明的另一个实施方案，该方法包括：如果时间点不是感兴趣时间点，则避免通过累积时间来计算指数，优选地通过使用以下要求：

118、如果读出：

119、

120、mm＝mm(t0)*exp(-h/t2(r))

121、mz＝m0z+(mz(t0)-m0z)*exp(-h/t1(r))

122、否则：

123、

124、h+＝dt。

125、因此，根据另一个实施方案，该方法包括跟踪该相位相对于空间方向的空间导数，优选地作为根据下式的额外变量：

126、

127、

128、

129、这对于一些相关情况是有用的，即捕获体素上的相位的变化以去除假回波(这需要在体素上进行相位积分)，或者计算需要这些导数的扩散效应[6]。

130、下文提供了根据上文公开的本发明的该第二实施方案方面的概念的含义和优势。

131、1)在某些情况下，不能应用如上所述的根据第一实施方案方面的序列压缩，并且该方法可以是有用的。例如，在具有运动或流动的模型中[5]，体素的位置随时间变化。替代应用更复杂的方法，即使体素位置r随时间变化，也可很容易地累积所产生的相位。利用小的时间步长，可非常准确地跟踪该移动，从而生成不是感兴趣点的许多时间点。每当需要计算磁化(感兴趣点)时，可使用累积的相位和时间步长来生成信号。

132、2)性能改进：在不能应用根据第一方面的方法的此类情况下，由于不需要计算指数或旋转(其需要评估sin和cos)，所以计算量大大减少。这种操作比简单的乘积和求和昂贵得多。

133、3)通过跟踪相位随时间的变化，可以很容易地解释由于体素上的相位变化而产生的效应。利用笛卡尔坐标，实际相位丢失，并且只能在[0，2π]或[-π，π]区间中获得。这可能产生难以处理的相位缠绕，而这可能会限制扩散和其它相位效应的正确计算。

134、应当理解，根据本发明，如上所解释的概念和方法可独立应用，或者可以组合应用。

135、当组合时，相位累积直接来自于每个信号的预先计算的面积a。

136、

137、其中ax、ay和az是在周期h期间梯度信号的面积积分。

138、并且在这种情况下，将只需要计算读出时间点：

139、对于每个读出：

140、

141、mm＝mm(t0)*exp(-h/t2(r))

142、mz＝m0z+(mz(t0)-m0z)*exp(-h/t1(r))

143、其中mm(t0)、和mz(t0)是梯度部分的开始(t0)处的初始条件。