时域磁共振的参数图确定的制作方法

时域磁共振的参数图确定
1.本专利公开内容涉及用于基于根据施加的脉冲序列激励样本之后从样本发射的时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布的方法和设备，以及用于执行该方法的计算机程序产品。
2.磁共振成像(mri)是用于许多应用的成像方式，并且具有许多可以调整的序列参数和许多可以观察到的成像参数以提取例如不同种类的生物信息。常规的mri图像重建涉及获取k空间信号并且对获取的数据执行快速傅立叶逆变换(eft)。常规的mri成像慢，这是因为对于要测量的每个参数(例如t1或t2)，要使用具有不同设置的mri设备获取单独的mri测量。扫描可能需要多达30分钟至45分钟。
3.时域中的磁共振自旋断层扫描(mr-stat)是用于直接从时域数据获得mr图像的定量方法。特别地，mr-stat是用于使用来自单个短扫描的数据获得多参数定量mr图的框架。
4.如wo 2016/184779 a1中所描述的，解决了大规模优化问题，其中通过使用不精确的高斯牛顿(gauss-newton)方法将基于布洛赫(bloch)的体积信号模型直接拟合至时域数据来同时执行信号的空间定位和组织参数的估计。高度并行化、无矩阵的不精确高斯牛顿重建算法可以用于解决针对高分辨率扫描的大规模优化问题。
5.g.w
ü
bbeler等人，在“a large-scale optimization method using a sparse approximation of the hessian for magnetic resonance fingerprinting”,siam journal on imaging sciences,vol.10,nr.3,18july 2017,pp.979-1004中，描述了用于磁共振指纹(mrf)的最小二乘法，其中使用了黑塞矩阵(hessian)的稀疏近似。
6.虽然测量时间大大减少至几分钟的量级，但是用于将信号模型拟合至时域数据的计算时间对于单个2d切片来说是约1小时。
7.本专利公开内容的目的之一是用于改进时域mr信号至定量mr图的转换。
8.根据第一方面，提供了一种方法，优选地为计算机实现的方法，该方法用于基于时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布，该tdmr信号是在根据施加的脉冲序列激励样本之后从样本发射的，该方法包括：
9.i)确定tdmr信号模型以近似发射的时域磁共振信号，其中，tdmr信号模型被配置成取决于tdmr信号模型参数，所述tdmr信号模型参数包括：
10.样本内的至少一个组织参数的空间分布；以及
11.施加的脉冲序列；
12.ii)估计tdmr信号模型参数的初始集合；
13.iii)基于作为输入的tdmr信号模型参数的初始集合，确定表示发射的tdmr信号与tdmr信号模型之间的差异的残差向量；
14.iv)通过下述方式确定tdmr信号模型参数的更新集合
15.v)计算残差向量的梯度；
16.vi)计算tdmr信号模型的稀疏黑塞矩阵，所计算的稀疏黑塞矩阵具有使用施加的脉冲序列计算的稀疏模式；以及
17.vii)基于所计算的梯度和所计算的黑塞矩阵更新tdmr信号模型参数的初始集合；
18.viii)针对作为输入的tdmr信号模型参数的更新集合，重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)直到发射的tdmr信号与tdmr信号模型之间的差异小于预定阈值为止或者直到完成预先确定的重复的次数为止，从而获得tdmr信号模型参数的最终更新集合；以及
19.ix)从tdmr信号模型参数的最终更新集合中提取至少一个组织参数的空间分布。
20.当重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)时，使用tdmr信号模型参数的更新集合代替tdmr信号模型参数的初始集合。以这种方式，tdmr信号模型参数的更新集合迭代地收敛至tdmr信号模型参数的最终更新集合。通常，可以使用其中上述差异小于预定阈值和预先确定的重复的次数表示示例的任何停止标准。因此，重复150步骤可以完成直到满足停止标准为止。
21.在本公开内容中使用“稀疏”来指示作为矩阵的黑塞矩阵包含为零或可忽略的元素(例如元素的总数量的一半以上)。例如，稀疏黑塞矩阵的稀疏性是黑塞矩阵的零值和可忽略元素的数量除以黑塞矩阵的元素的总数量。在步骤vi)中的计算的黑塞矩阵的稀疏性基于施加的脉冲序列。更具体地，黑塞矩阵具有稀疏模式，例如包括非零的和不可忽略的元素的对角带的稀疏模式。选择施加的脉冲序列，使得黑塞矩阵是具有特定稀疏模式的稀疏黑塞矩阵。
22.在本方法中，与不精确的高斯牛顿方法相比，计算时间减少了至少10倍。在不精确的高斯牛顿方法中，为了近似黑塞矩阵，必须计算信号相对于每个参数的偏导数，这需要大量的计算时间。相反，在本方法中，对于tdmr信号模型参数的每个集合，作为稀疏黑塞矩阵的黑塞矩阵被预先计算并存储，使得不需要针对每个参数计算这些偏导数。此外，本方法的内存需求可以低至不精确的高斯牛顿方法的内存需求的0.04％。
23.在g.w
ü
bbeler等人中，稀疏黑塞矩阵的近似是独立于施加的脉冲序列的近似，这是因为p
lt
p
l
≈(m
l
/n)in的近似致使稀疏黑塞矩阵总是仅具有稀疏黑塞矩阵的每个对角带的单一参数。另一方面，在本技术中，稀疏黑塞矩阵是基于施加的脉冲序列来确定的，并且允许带宽大于1。这致使所需的计算时间至少减少一个或者甚至两个数量级。
24.注意，虽然确定了残差向量，但是该残差向量也可以被命名为残差函数。
25.优选地，在将tdmr信号模型的黑塞矩阵计算为稀疏黑塞矩阵时，黑塞矩阵的稀疏模式基于施加的脉冲序列。
26.在实施方式中，施加的脉冲序列包括梯度编码模式和/或射频激励模式，并且黑塞矩阵的稀疏模式是基于梯度编码模式和/或射频激励模式确定的。
27.在实施方式中，黑塞矩阵被计算为残差向量的雅可比行列式(jacobian)和残差向量的雅可比行列式的厄尔米特(hermitian)转置的乘积。代替残差向量的雅可比行列式，它将等同于采用tdmr信号模型的雅可比行列式和对应的厄尔米特转置。梯度编码模式确定k空间采样轨迹，例如笛卡尔、径向、螺旋。
28.在实施方式中，施加的脉冲序列被配置成产生笛卡尔获取、径向获取或螺旋获取中的任一个。对于这些类型的脉冲序列，稀疏模式适当地包括黑塞矩阵中不可忽略的元素的多个对角带。其他元素优选为零或可忽略的。
29.根据实施方式，施加的脉冲序列的梯度编码模式被配置成产生笛卡尔获取，使得对应的点扩散函数仅在相位编码方向上传播。还可能存在两个点扩散函数，例如对于3d获
取，每个点扩散函数与两个相位编码方向中的对应一个相关联。这两个相位编码方向中的每一个是不同的。
30.在实施方式中，施加的脉冲序列被配置成产生变化的翻转角。
31.优选地，施加的脉冲序列的射频激励模式被配置成产生平滑变化的翻转角，使得对应的点扩散函数在宽度方向上在空间上受到限制。平滑变化可以指示其中rf激励的幅值随时间变化有限量的序列。在k空间的(或每个k空间的)采样期间，两个连续rf激励之间的变化量小于预先确定量，优选地小于5度。
32.更优选地或替选地，施加的脉冲序列的射频激励模式被配置成在k空间或多个k空间的每个k空间的采样内产生平滑变化的翻转角，使得对应的点扩展函数在宽度方向上在空间上受到限制。
33.在实施方式中，至少一个组织参数包括t1弛豫时间、t2弛豫时间、t2*弛豫时间和质子密度，或者它们的组合中的任一个。
34.在实施方式中，计算残差向量的梯度包括将残差向量与残差向量的雅可比行列式的厄尔米特转置相乘。代替残差向量的雅可比行列式，它将等同于采用tdmr信号模型的雅可比行列式和对应的厄尔米特转置。
35.在实施方式中，tdmr信号模型是基于布洛赫的体积信号模型。根据阅读本公开内容将明显的是，该方法可以用替选的体积信号模型来执行。
36.在实施方式中，黑塞矩阵中的稀疏模式由施加的脉冲序列的梯度编码模式来确定。可以选择施加的脉冲序列以获得合适稀疏的黑塞矩阵，特别是具有合适的稀疏模式的黑塞矩阵。
37.在实施方式中，黑塞矩阵被计算为包括多个对角带的稀疏黑塞矩阵。稀疏模式可以包括对角带的数量。换言之，黑塞矩阵的不可忽略的元素的对角带表示一种稀疏模式。
38.在实施方式中，黑塞矩阵中的对角带的数量由施加的脉冲序列的梯度编码模式确定。作为第一示例，确定笛卡尔获取为每个模型参数提供一个对角带的合适稀疏模式。作为第二示例，确定径向获取和螺旋获取为每个模型参数提供4个对角带的合适稀疏模式。
39.在实施方式中，多个对角带中的任一个的宽度由施加的脉冲序列的射频激励模式来确定。
40.在实施方式中，多个对角带中的任一个的宽度大于1，其中，多个对角带中的任一个的宽度优选地在2至55的范围内，多个对角带中的任一个的宽度更优选地在3至8的范围内，多个对角带中的任一个的宽度最优选地为3。宽度可以指示稀疏黑塞矩阵的每行的数字元素/系数。
41.在实施方式中，稀疏黑塞矩阵对于tdmr信号模型参数的每个集合是恒定的。优选地，通过基于稀疏黑塞矩阵和残差函数的梯度求解线性方程组来获得tdmr信号模型参数变化的集合。
42.在实施方式中，更新的步骤vii)包括：
43.通过求解稀疏黑塞矩阵与tdmr信号模型参数变化相乘等于求负的残差函数的梯度的线性方程组，获得tdmr信号模型参数变化的集合，以及
44.将所获得的tdmr信号模型参数变化与tdmr信号模型参数被更新的集合相加以获得tdmr信号模型参数的更新集合。
45.根据本专利公开内容的第二方面，提供了一种用于基于时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布的设备，该tdmr信号是在根据施加的脉冲序列激励样本之后从样本发射的，该设备包括处理器，该处理器被配置成：
46.i)确定tdmr信号模型以近似发射的时域磁共振信号，其中，tdmr信号模型被配置成取决于tdmr信号模型参数，该tdmr信号模型参数包括：
47.样本内的至少一个组织参数的空间分布；以及
48.施加的脉冲序列；
49.ii)估计tdmr信号模型参数的初始集合；
50.iii)基于作为输入的tdmr信号模型参数的初始集合，确定表示发射的tdmr信号与tdmr信号模型之间的差异的残差向量；
51.iv)通过下述方式确定tdmr信号模型参数的更新集合
52.v)计算残差向量的梯度；
53.vi)计算tdmr信号模型的稀疏黑塞矩阵，所计算的稀疏黑塞矩阵具有使用施加的脉冲序列计算的稀疏模式；以及
54.vii)基于所计算的梯度和所计算的黑塞矩阵更新tdmr信号模型参数的初始集合；
55.viii)针对作为输入的tdmr信号模型参数的更新集合，重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)直到发射的tdmr信号与tdmr信号模型之间的差异小于预定阈值为止或者直到完成预先确定的重复的次数为止，由此获得tdmr信号模型参数的最终更新集合；以及
56.ix)从tdmr信号模型参数的最终更新集合中提取至少一个组织参数的空间分布。
57.如将明显的，根据第二方面的设备特别地被配置成用于应用上面和/或下面所描述的方法步骤中的任一个或多个。另外，将明显的是，针对本文中所描述的方法和方法步骤所提到的任何优点都适用于该设备，并且针对该设备所提到的优点适用于对应的方法和方法步骤。
58.根据另一方面，提供了包括计算机可执行指令的计算机程序产品，当该程序在计算机上运行时，该计算机可执行指令用于根据上面和/或下面所描述的任一个实施方式的任一个步骤来执行任一个方法的方法。
59.根据另一方面，提供了包括计算机可执行指令的计算机程序，当该程序在计算机上执行时，该计算机可执行指令用于根据上面和/或下面所描述的任一个实施方式的任一个步骤来执行该方法。
60.根据另一方面，提供了被编程成执行上面和/或下面所描述的方法中的任一个实施方式的一个或更多个步骤的计算机设备或其他硬件设备。
61.根据又一方面，提供了以机器可读和机器可执行形式编码程序的数据存储设备，以用于执行上面和/或下面所描述的方法的任一个实施方式的一个或更多个步骤。
62.附图用于说明本公开内容的设备的当前优选的非限制性示例性实施方式。在结合附图阅读的情况下，根据以下详细描述，本公开内容的特征和目的的以上内容和其他优点将变得更加明显，并且将更好地理解各方面和实施方式，在附图中：
63.图1是示出确定空间分布的示例性方法的实施方式的流程图；
64.图2是示出确定图1的方法的更新的tdmr信号模型参数集合的实施方式的流程图；
65.图3是针对其中使用笛卡尔采样策略的施加的脉冲序列的稀疏黑塞矩阵的示意图；
66.图4是针对其中使用非笛卡尔采样策略的施加的脉冲序列的稀疏黑塞矩阵的示意图；
67.图5是示出示例施加的脉冲序列的(顶部)翻转角和(底部)对应的相位编码步骤与激励数量的图；
68.图6是稀疏黑塞矩阵的彩色示意图，其中颜色指示稀疏黑塞矩阵的每个元素的幅值的log
10
；
69.图7是另一稀疏黑塞矩阵的灰度示意图，其中灰色的阴影指示稀疏黑塞矩阵的每个元素的幅值的log
10
；
70.图8是相关技术方法与本专利公开内容的示例性方法的目标函数的对数与计算时间的对数的对数图；
71.图9是(左)如通过本公开内容的示例性方法获得的两个重建的体内t1图和t2图的集合，(中)相关技术方法的两个重建的体内t1图和t2图的集合，以及(右)比较方法的获得的平均t1值和平均t2值的比较；
72.图10是示出更新图1的方法的tdmr信号模型参数集合的实施方式的流程图；
73.图11是示出更新图1的方法的tdmr信号模型参数集合的另一实施方式的流程图；
74.图12是根据本专利公开内容的用于执行确定空间分布的示例性设备的框图；
75.图13是执行样本的时域磁共振测量和分析的方法的流程图；
76.图14是用于执行样本的时域磁共振测量的示例性设备的框图；以及
77.图15是本专利公开内容的方法的目标函数的对数与计算时间的对数的对数图，其中，稀疏黑塞矩阵的对角带的宽度是变化的。
78.参照图1，示出了用于基于时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布的方法100，该时域磁共振tdmr信号是在根据施加的脉冲序列激励样本之后从样本发射的。该方法100包括下述步骤：
79.i)确定110tdmr信号模型以近似发射的时域磁共振信号，其中，tdmr信号模型被配置成取决于tdmr信号模型参数，该tdmr信号模型参数包括：
80.样本内的至少一个组织参数的空间分布；以及
81.施加的脉冲序列；
82.ii)估计120tdmr信号模型参数的初始集合；
83.iii)基于作为输入的tdmr信号模型参数的初始集合，确定130表示发射的tdmr信号与tdmr信号模型之间的差异的残差向量；
84.iv)基于tdmr信号模型参数的输入集合，确定140tdmr信号模型参数的更新集合。tdmr信号模型参数的输入集合在第一迭代时是tdmr信号模型参数的初始集合，并且之后是tdmr信号模型参数的更新集合。信号模型参数因此被迭代地更新以获得拟合结果。确定140包括：
85.v)计算142残差向量的梯度；
86.vi)计算144tdmr信号模型的黑塞矩阵作为稀疏黑塞矩阵；以及
87.vii)基于计算的梯度和计算的黑塞矩阵更新146tdmr信号模型参数的输入集合，
从而获得tdmr信号模型参数的更新集合；
88.viii)针对tdmr信号模型参数的更新集合，重复150步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)，直到确定152发射的tdmr信号与tdmr信号模型之间的差异小于预定阈值为止或者直到完成预先确定的重复的次数为止，从而获得tdmr信号模型参数的最终更新集合。通常，可以使用任何停止标准。因此，重复150可以完成直到满足停止标准为止。
89.此后，执行从tdmr信号模型参数的最终更新集合中提取160至少一个组织参数的空间分布的步骤ix)。
90.提取的空间分布允许获得mr分布图，或者换言之，mr图像例如样本的切片等。
91.当重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)时，使用tdmr信号模型参数的更新集合代替tdmr信号模型参数的初始集合。以这种方式，tdmr信号模型参数的更新集合迭代地收敛至tdmr信号模型参数的最终更新集合。
92.在将tdmr信号模型的黑塞矩阵计算为稀疏黑塞矩阵时，稀疏黑塞矩阵可能具有基于施加的脉冲序列的稀疏模式。
93.在实施方式中，基于梯度编码模式和/或射频激励模式来确定黑塞矩阵的稀疏模式。
94.在实施方式中，稀疏黑塞矩阵的计算144是通过将黑塞矩阵计算为残差向量的雅可比行列式和残差向量的雅可比行列式的厄尔米特转置的乘积来完成的。代替残差向量的雅可比行列式，它将等同于采用tdmr信号模型的雅可比行列式和对应的厄尔米特转置。
95.在实施方式中，残差向量的梯度的计算142包括将残差向量与残差向量的雅可比行列式的厄尔米特转置相乘。代替残差向量的雅可比行列式的厄尔米特转置，它将等同于采用tdmr信号模型的厄尔米特转置。
96.在实施方式中，黑塞矩阵的计算144是通过将黑塞矩阵计算为包括多个对角带的稀疏黑塞矩阵来完成的。稀疏模式可以包括对角带的数量。换言之，黑塞矩阵的不可忽略的元素的对角带表示一种类型的稀疏模式。其他元素为零或以其他方式可忽略。例如，如果元素的(绝对)值小于阈值，则该元素可以被认为是可忽略的。该阈值可以是预先确定的。
97.在实施方式中，黑塞矩阵中的对角带的数量由施加的脉冲序列的梯度编码模式确定。作为第一示例，确定笛卡尔获取为每个模型参数提供一个对角带的合适稀疏模式。作为第二示例，确定径向和螺旋获取为每个模型参数提供4个对角带的合适稀疏模式。
98.在实施方式中，多个对角带中的任一个的宽度由施加的脉冲序列的射频激励模式确定。
99.作为示例，如图2所示，在包括稀疏黑塞矩阵的计算144的tdmr信号模型参数的更新集合的确定140中，可以包括下述：确定210不可忽略的黑塞矩阵的元素/系数。不可忽略可以指示例如残差向量的雅可比行列式与残差向量的雅可比行列式的厄尔米特转置的乘积的每个元素的绝对值大于预先确定的值，例如0.1、0.2、0.3、0.4、0.5、
……
、1。可以选择预先确定的值以允许足够准确的参数图的重建。确定210的步骤是可选的，因为这可以针对施加的脉冲序列预先确定。然后，仅对不可忽略的元素进行评估220和存储230。所得矩阵是对真实黑塞矩阵的稀疏近似黑塞矩阵由于对于本方法和施加的脉冲序列，仅一小部分系数是不可忽略的，因此可以快速计算和存储该近似。由于稀疏模式(例如，不可
忽略的sah系数的索引)与获取策略之间存在明确的关系，因此人们预先知道需要计算的系数。所有其他(可忽略的)系数都不被计算。这致使节省大量的计算和内存。
100.梯度编码模式确定了sah矩阵的稀疏模式。对于笛卡尔获取，点扩散函数仅在相位编码方向上传播。这导致sah矩阵中的带状对角块。图3示出了用于笛卡尔获取的针对sah矩阵300的确定的稀疏模式的示例。白色元素指示相应元素是可忽略的。然后，这些元素可能不会在其中使用sah的任何计算步骤中进行评估。黑色元素指示相应元素是不可忽略的。轴302和轴304指示sah矩阵300的相应的行和列。轴302和轴304可以被认为是[每个体素的参数]
×
[体素]。在该示例中，存在由虚线分隔的3
×
3块310。sah矩阵300具有不可忽略的元素的对角带320。块的数量通常为n
×
n，其中n为要重建的参数类型的数量，例如，如果重建t1、t2和质子密度，则n＝3。
[0101]
图4以关于图3类似的方式示出了用于非笛卡尔获取例如螺旋和/或径向获取的sah矩阵400的确定的稀疏模式的示例。在这种情况下，点扩散函数在两个方向上传播，因此它导致针对每个块的若干个带420、422和424。类似于轴302和轴304，轴402和轴404指示sah矩阵400的相应的行和列。存在由虚线分隔的3
×
3块410。sah矩阵400具有不可忽略的元素的主对角带420以及对称子带422和424。由单独地标记的元素形成的线的厚度或宽度示出相应的sah矩阵元素的绝对数量的值。
[0102]
从数据获取的角度来看，非笛卡尔序列更有效，而笛卡尔序列具有对于硬件缺陷(例如涡流)更加稳健的益处。对于适合于在mr-stat中使用的梯度平衡序列尤其是这种情况。
[0103]
射频激励(翻转角列)与稀疏模式和平滑度之间的示例关系
[0104]
翻转角(rf)激励的平滑度确定了稀疏模式中的带的宽度，例如如以上针对图4所描述的。对于平滑的翻转角列，仅考虑薄带将足够。这意味着sah矩阵的每行仅需要评估和存储很少的元素/系数。因此，当翻转角列平滑时，以上附图中的带是薄的。平滑的示例是关于每个激励的翻转角的变化等于或小于5度。对于快速变化的翻转角列，情况相反，也就是说需要更宽的带。这意味着需要评估和存储sah矩阵的每行的更多系数。例如，图3和图4中的带320、420、422和424将因此更宽。
[0105]
翻转角列的所需的平滑度取决于样本、用于进行测量的设备、其设置等，但是优选在翻转角中具有一些程度的变化，以能够准确地量化参数例如t1、t2和质子密度。翻转角列的一个示例在图5中示出。该图5指示了用于多参数定量序列的翻转角列的平滑度的合适度的示例。特别地，图5示出了具有不同翻转角和线性笛卡尔采样策略的瞬态2d平衡梯度回波序列。
[0106]
现在返回至图1的方法100，并且现在参照图10，在优选实施方式中，更新146的步骤vii)包括获得520tdmr信号模型参数变化的集合，该tdmr信号模型参数变化的集合是通过求解510基于稀疏黑塞矩阵和残差函数的梯度的线性方程组而获得的。
[0107]
现在参照图11，在更优选的实施方式中，更新146的步骤vii)包括：通过求解630稀疏黑塞矩阵与tdmr信号模型参数变化相乘等于求负的残差函数的梯度的线性方程组，来获得620tdmr信号模型参数变化的集合，并且将获得的tdmr信号模型参数变化与tdmr信号模型参数被更新的集合相加640，以获得tdmr信号模型参数的更新集合。
[0108]
参照图12，设备700是用于基于根据施加的脉冲序列激励样本之后从样本发射的
时域磁共振tdmr信号确定样本内的至少一个组织参数的空间分布的设备的实施方式，该设备700包括被配置成执行以上所描述的方法中的任何一种或更多种的处理器710。设备700可以包括用于存储执行方法步骤所需的任何模型、参数和/或其他数据的存储介质720。存储介质720还可以存储可执行代码，当由处理器710执行时，该可执行代码执行如以上所描述的方法步骤中的一个或更多个。设备700还可以包括用于接收用户输入的网络接口730和/或输入/输出接口740。尽管示出为单个设备，但是设备700还可以被实现为设备例如用于并行计算的系统或超级计算机等的网络。
[0109]
参照图13，示出了包括激励1310样本、接收1320信号以及在时域中分析1330信号的方法1300。分析1330可以包括方法100。
[0110]
如图14所示，方法步骤1310和方法步骤1320可以由mri设备1400执行。设备1400包括梯度线圈1420、用于瞬时激励样本1401从而使样本1401发射发射信号的激励设备1400、以及用于接收发射信号的接收线圈1403。mri设备1400可以从样本1401获得信号。
[0111]
激励设备1400和接收线圈1403在mri领域中是公知的，并且本文中不再对其进行更详细地描述。设备1400还包括控制包括激励设备1400和接收线圈1403的系统的部件的处理器1405。处理器1405还可以控制用于输出图像和/或状态信息的显示器1404，以及用于接收来自用户的命令和辅助信息的输入设备1406例如触摸屏、鼠标和键盘。处理器1405可以包括被配置成协作以执行任务例如计算的多个处理设备。替选地，处理器1405由单个处理设备组成。这样的处理设备例如中央处理单元(cpu)、控制器或fpga在本领域中是已知的。为了不使描述模糊，在附图和该描述中省略了mri设备的一些公知的元件。
[0112]
设备1400还包括存储装置，例如存储器1407。存储器1407可以被配置成在处理器1405的控制下存储从接收线圈1403接收的信号1412。重建的图像数据1413可以由设备1400本身生成，并且然后包括用于执行如本文中所描述的方法的所有必要特征，例如设备700的所有特征。图像数据也可能由外部设备例如设备700重建。重建是通过由本文中所描述的方法例如方法100处理接收到的信号1412来完成的。存储器1407还可以存储计算机代码以使处理器1405执行其任务。例如，计算机代码可以包括用于基于诸如图1所示出的获取方案操作数据获取的图像获取模块1408。该图像获取模块1408可以使激励设备1400瞬时激励样本1401，从而使样本发射发射信号。此外，图像获取模块1408可以使处理器1405接收并且存储来自接收线圈403的发射信号，该接收线圈403接收由样本401发射的信号。
[0113]
示例实现方式
[0114]
时域中的mr自旋断层扫描(mr-stat)是用于使用来自单个短扫描的数据获得多参数定量mr图的框架。解决了大规模优化问题，其中通过将基于布洛赫的体积信号模型直接拟合至时域数据来同时执行信号的空间定位和组织参数的估计。在相关技术方法中，使用高度并行化、无矩阵的高斯牛顿重建算法可以用于解决针对高分辨率扫描的大规模优化问题。
[0115]
理论
[0116]
在本节中，示出了关于相关技术方法的mr-stat框架，并且总结了先前如何解决重建过程中的一些计算挑战。此后，对于施加的脉冲序列的示例，从理论上推导出黑塞矩阵具有在本方法中采用的稀疏结构以加速mr-stat重建。
[0117]
标量(实数和复数两者)以小写字母表示，向量以粗体小写字母表示，并且矩阵以
粗体大写字母表示。
[0118]
mr-stat框架
[0119]
其中空间坐标为r＝(x，y，z)的单个自旋等色线m＝(m
x
，my，mz)的时间演化可以通过布洛赫方程来描述。该时间演变取决于所使用的脉冲序列(例如，rf激励脉冲和作用在自旋等色线上的空间梯度)并且也取决于其与mr相关的生物物理组织特性θ＝(t1，t2，...)。设与脉冲序列相关联的梯度轨迹由k(t)表示。
[0120]
设m＝m
x
+imy为旋转坐标系中自旋等色线的横向分量。将通过mr系统获得的解调时域信号s建模为视场v内所有自旋等色线的横向磁化强度的体积积分。也就是说，
[0121]
s(t)＝∫vm(θ(r)，t)e-2πik(t)
·rdr.
ꢀꢀ
(1)
[0122]
在完全采样的稳态序列中，横向磁化强度失去其时间依赖性，并且可以使用fft来恢复定性图像。在瞬态序列的更一般情况下，不能再直接使用fft以在测量的时域(或k空间)数据与图像空间之间进行变换。因此，在瞬态情况下，执行以下操作。首先，视场v被离散化为nv个体素，并且方程(1)变为
[0123][0124]
此处，mj为体素j中的磁化强度。对于梯度平衡序列，可以通过对单个自旋等色线的布洛赫方程的数值积分来计算mr。
[0125]
设n
t
为通过mr系统的接收器(例如接收线圈1403)获取的样本的总数量，并且设表示采样时间。如果将体素j中的磁化强度向量mj限定为：
[0126][0127]
并且对于体素j的梯度编码向量为：
[0128][0129]
那么可以将离散化的信号向量限定为：
[0130][0131]
此处，
⊙
表示逐分量向量乘法。信号向量s是上面所描述的tdmr信号模型的示例。现在设n
p
表示每个体素的不同参数的数量(例如，包括t1、t2以及质子密度的实部和虚部)。则s取决于n：＝nv×np
不同的参数。所有参数以参数{j+knv|k＝0...，n
p-1}是与体素j相关联的参数这样的方式连接成单个向量
[0132]
给出测量的时域样本的向量这是如以上所描述的从样本发射的tdmr信号的示例，可以将残差向量限定为
[0133]
r(α)＝d-s(α)
ꢀꢀ
(6)
[0134]
残差向量r是针对方法100上面所描述的残差向量的示例。可以将目标函数限定为：
[0135][0136]
如以上所描述的，目标函数f可以用于以绝对术语描述发射的tdmr信号与tdmr信号模型之间的差异。参数图α
*
是如以上所描述的针对组织参数的空间分布的示例，该参数图α
*
通过数值求解获得
[0137]
α
*
＝argmin
α
f(α)，
ꢀꢀ
(8)
[0138]
受到由布洛赫方程和针对参数可达到的区间表示的物理约束。
[0139]
不精确的高斯牛顿方法
[0140]
注意，方程(8)是非线性优化问题。这样的问题可以使用基于牛顿的方法来解决。基于牛顿的方法从初始猜测α开始，并且然后通过求解线性系统获得参数空间中的更新步骤p
[0141]
h(α)p＝-g(α).
ꢀꢀ
(9)
[0142]
此处，是目标函数的梯度，并且是被限定为下述的黑塞矩阵：
[0143][0144]
其中在mr-stat中直接应用牛顿方法的困难是问题的固有大规模。即使对于2d问题，参数n的数量也可以是以106的量级。因此，在当今的计算机体系结构上，显式地形成黑塞矩阵或其逆矩阵是不可行的。已经开发了对牛顿方法的多种修改，其使用针对(逆)黑塞矩阵的近似而不是(逆)黑塞矩阵本身。一种适用于大规模优化问题的这样的方法是有限内存的bfgs方法(“l-bfgs”)。该方法直接近似逆黑塞矩阵，但是仅需要存储有限数量的来自先前迭代的梯度向量。替选地，最小二乘问题中通常使用的技术是用于将黑塞矩阵近似为其中是被限定为下述的雅可比行列式矩阵：
[0145][0146]
jh是j的厄尔米特转置，并且是实部算子。即使矩阵j通常太大而无法存储在计算机内存中，也可以计算jv和jhv形式的矩阵向量乘积而无需显式存储j。给出用于计算这些矩阵向量乘积的能力，然后可以使用共轭梯度技术以迭代方式求解其中h被替换的方程9中的线性系统。而不是求解方程9至任意精度，在该内部循环中执行的迭代的次数——即，针对α的每个设置求解方程9——是有限的，致使所谓的不精确高斯牛顿方法。该不精确的高斯牛顿方法优于l-bfgs方法，并且可以用于重建高分辨率的参数图。针对该无矩阵高斯牛顿方法的伪算法1如下：
[0147]
伪算法1的开始
[0148]
要求：初始猜测α
[0149]
(“外层循环”)
[0150]
while！收敛do
[0151]
1.计算残差：r(α)＝d-s(α)
[0152]
2.计算梯度：
[0153]
3.(“内层循环”)以无矩阵方式求解线性系统：
[0154]
4.更新参数：α＝α+p
[0155]
end while
[0156]
伪算法1的结束
[0157]
在实践中，观察到无矩阵高斯牛顿方法在参数空间p中产生良好的更新步骤，但是该方法确实需要针对每个矩阵向量乘法重新计算j的列。也就是说，在内部循环的每次迭代中，计算信号相对于n参数中的每一个的偏导数。这些计算形成了重建算法的计算瓶颈。
[0158]
稀疏黑塞矩阵近似
[0159]
针对与所选择的施加的脉冲序列对应的tdmr信号，具有稀疏结构，例如当使用笛卡尔获取时为带状。然后可以预先计算和存储(每次外部迭代一次)。然后，随后与稀疏相乘的计算量是可忽略的，并且方程9中的线性系统可以通过迭代算法快速求解以获得p，致使更快的mr-stat重建。
[0160]
的稀疏性可以表述如下。矩阵包含有关不同参数之间的依赖性的信息。在不精确的高斯牛顿方法中，该矩阵被反转以便“重新聚焦”参数依赖性。对于完全采样的稳态序列，已经通过梯度编码消除了与不同体素相关联的参数之间的依赖性。也就是说，如果参数j和k属于同一个体素(mod(j-k，nv)＝0)，则矩阵项仅为非零。因此，将由对角矩阵的个块组成，对角矩阵的一个块针对每对不同的参数类型。
[0161]
在瞬态序列的情况下，获取的数据可以被解释为来自顶部上具有扰动的稳态序列的数据。如果对这种数据应用fft，扰动将致使与图像空间中的点扩散函数(扰动的fft)的卷积。也就是说，源自一个体素的信号可能泄漏至其他体素中。然而，在笛卡尔采样策略的情况下，与获取相关联的点扩散函数在空间上(良好近似)受限于相位编码方向。如果除此之外采用平滑变化的翻转角(致使平滑的信号响应)，点扩散函数的宽度也受到限制。矩阵则更加稀疏。
[0162]
稀疏模式的形式推导的示例如下。首先，体素j中的磁化强度可以被表示为mj(其取决于rf激励脉冲和体素j中的参数，例如t1，t2，......)与包含来自空间编码梯度的相位项的向量grj的逐分量乘积。
[0163]
1.为了简化记号，假定每个体素仅具有一个相关联的参数，使得可以使用相同的索引来表示体素j以及其对应的参数αj。由于grj与组织参数无关，由此得出
[0164][0165]
在使用规则的、完全采样的稳态序列的情况下，mj实际上与时间无关，并且每个元素具有要由mj表示的相同的值。
[0166][0167]
由此得出，jhj是对角矩阵。
[0168]
2.对于每个体素具有多个参数的另一示例(n
p
＞1)，矩阵jhj将由对角矩阵的个块组成。
[0169]
3.现在，对于另一示例，其中再次为了简单起见，n
p
＝1，替代地使用瞬态序列。内积可以被写为测量数据集(tdmr信号)上的内积的和，如下：
[0170][0171]
针对方程14中术语如果读数是笛卡尔的，由此得出
[0172][0173]
也就是说，消除了与具有不同x坐标的体素相关联的参数之间的依赖性。
[0174]
4.假定r
j，x
＝r
k，x
，方程14的方括号中的另一项可以被写为内积并且使用帕赛瓦尔定理，由此得出：
[0175]
[0176]
如果施加的脉冲序列的翻转角列是平滑的，其示例在图5中示出，则信号响应(从样本发射的tdmr信号)以及信号响应w.r.t.参数的偏导数也是平滑的。而且傅立叶变换具有有限的带宽。对于沿y方向具有足够大间距的体素，内积将然后为零。因此，仅针对与沿相位编码方向(相同的x坐标)在相同线上并且彼此接近的体素相关联的参数存在相关性。这致使矩阵jhj的显著稀疏性。
[0177]
5.针对情况n
p
＝1，并且假定坐标被排序为[(x1，y1)，(x1，y2)，...，(x2，y1)，(x2，y2)，...])，矩阵jhj将为带状矩阵。针对n
p
＞1，矩阵jhj将由带状矩阵(由于对称性，其中仅n
p
(n
p
+1)/2)是唯一的)的个块组成。
[0178]
在图6中，针对小的数值模型(针对其可以计算完整的jhj)，示出了jhj的幅值的对数图。在图6中可以看到，所示出的jhj具有块结构(6
×
6块，一个块针对每对不同的参数)并且每个块近似地是稀疏带状矩阵。在图7中，以与图6相同的比例示出了针对更多参数(并且因此存在更多块)的jhj的幅值的类似对数图，但是以灰度示出。
[0179]
在伪算法2中，示出了作为上面所描述的方法100的示例的稀疏高斯牛顿方法的示例。
[0180]
伪算法2的开始
[0181]
要求：初始猜测α
[0182]
(“外层循环”)
[0183]
while！收敛do
[0184]
1.计算残差：r(α)＝d-s(α)
[0185]
2.计算梯度：g＝jhr
[0186]
3.计算稀疏黑塞矩阵：
[0187]
4.(“内层循环”)以无矩阵方式求解线性系统：
[0188]
5.更新参数：α＝α+p
[0189]
end while
[0190]
伪算法2的结束
[0191]
在该示例中，可以看到，在步骤4中，需要少得多的计算时间，这是因为该方法不需要针对内层循环中的每个矩阵向量乘法重新计算j的列。也就是说，内层循环的每次迭代都涉及求解线性系统
[0192]
在本文中，r是残差向量的示例，d是发射的tdmr信号的示例，s是tdmr信号模型的示例，g是残差函数的所计算梯度的示例，是稀疏黑塞矩阵的示例，jh是残差向量(或tdmr信号模型)的雅可比行列式的厄尔米特转置的示例，j是残差向量(或tdmr信号模型)的雅可比行列式的示例，p是tdmr信号模型参数变化的集合的示例，以及α是tdmr信号模型参数(更新的或初始的)集合的示例。
[0193]
在伪算法2中，步骤1是方法100中确定130残差向量的步骤iii)的示例。伪算法2的步骤2是方法100的计算142残差向量的梯度的步骤v)的示例。伪算法2的步骤3是方法100的计算144作为稀疏黑塞矩阵的tdmr信号模型的黑塞矩阵的步骤vi)的示例。伪算法2的步骤4和步骤5是方法100的更新146tdmr信号模型参数的输入集合的示例。特别地，伪算法2的步
骤4是通过求解630稀疏黑塞矩阵与tdmr信号模型参数变化相乘等于求负的图11的残差函数的梯度的线性方程组来获得620tdmr信号模型参数变化的集合(此处为p)的示例，以及伪算法2的步骤5是将获得的tdmr信号模型参数变化与tdmr信号模型参数被更新的集合相加640以获得图11中的tdmr信号模型参数的更新集合的示例。
[0194]
示例
[0195]
示例脉冲序列
[0196]
相关技术的无矩阵高斯牛顿方法和本公开内容的稀疏高斯牛顿(gn)方法对合成生成的数据、对来自凝胶体模的数据以及对体内大脑数据进行了比较。
[0197]
在所有情况下，使用具有平滑地改变每个tr(重复时间)和线性笛卡尔k空间填充的翻转角的2d梯度平衡的瞬态脉冲序列，如图5所示出的。对于合成实验，该序列用于模拟具有7.8秒的(模拟)扫描时间的数据。对于凝胶体模和体内大脑实验，瞬态序列在1.5飞利浦-英格尼亚(philips-ingenia)mr系统上实现，并且用于扫描具有13通道接收头线圈的健康志愿者。扫描时间为6.8秒，具有1
×1×
5mm3的扫描分辨率。
[0198]
根据模拟/获取的数据，重建t1,t2,和质子密度图。对于无矩阵gn方法，由于时间约束，内部迭代的次数被限制为10。对于稀疏gn方法，执行重建，其中，针对每个参数对计算15个次对角线的数量，致使与针对体内大脑数据集的完整矩阵幅值相比需要0.04％的存储。重建算法是用开源julia编程语言编写的，并且在使用64核的计算集群上运行。
[0199]
结果
[0200]
在图8中，示出了针对无矩阵和稀疏gn方法的收敛曲线。如可以看到的，获得与无矩阵gn方法相比的关于重建时间减少的数量级。
[0201]
图9示出了针对稀疏gn方法的重建的体内t1和t2图，以及与人类大脑的无矩阵gn方法的比较。在右侧示出了针对公知的参数gm左(在大脑的左侧的灰质)、gm右(在大脑的右侧的灰质)、wm左(在大脑的左侧的白质)和wm右(在大脑的右侧的白质)的平均t1值(顶部)和平均t2值(底部)。使用稀疏gn方法重建时间为50分钟，而对于无矩阵gn方法重建时间为320分钟。
[0202]
即使在稀疏gn方法中使用了jhj的近似(其本身近似黑塞矩阵)，优化算法也收敛至相同的参数图。事实上，该方法甚至可以导致更好的更新步骤(如图8中看到的)，因为对于内层循环的迭代的次数不受限制，允许本方法更好地考虑目标函数的曲率。
[0203]
根据以上将明显的是，对角带的宽度例如图3和图4的确定的稀疏黑塞矩阵中的对角带的宽度是基于施加的脉冲序列或换言之使用施加的脉冲序列来确定。其附加的一个优点是当宽度(或带宽)大于1时，计算时间减少。该效果在图15中示出，其中，可以看到宽度为1致使的计算时间比当施加的宽度为5时的计算时间长至少50倍或者甚至大于100倍。与图15中的(带)宽度(bw)5对应的图15的数据类似于图8中所示出的针对稀疏mr-stat的数据。例如，图15中所示出的计算可以使用线性笛卡尔采样和平滑rf列(例如图5中的一个)来完成，但是将明显的是，也可以使用其他类型的采样和rf列。
[0204]
本领域技术人员将容易认识到，可以通过编程的计算机来执行各种以上所描述的方法的步骤。在本文中，一些实施方式还旨在涵盖程序存储设备例如数字数据存储介质，其是机器或计算机可读的并且编码机器可执行的或计算机可执行的指令程序，其中所述指令执行所述以上所描述的方法的一些或全部步骤。程序存储设备可以是例如数字存储器、磁
存储介质例如磁盘和磁带、硬盘驱动器或光学可读数字数据存储介质。实施方式还旨在涵盖被编程为执行所述以上所描述的方法的步骤的计算机。
[0205]
图中所示出的各种元件的功能——包括被标记为“单元”、“处理器”或“模块”的任何功能块——可以通过使用专用硬件以及能够执行与适当软件相关联的软件的硬件来提供。当由处理器提供时，功能可以由单个专用处理器、由单个共享处理器或由多个单独的处理器提供，其中多个单独的处理器中的一些可以是共享的。此外，术语“单元”、“处理器”或“控制器”的明确使用不应被解释为专指能够执行软件的硬件，并且可以隐含地包括但不限于数字信号处理器(dsp)硬件、网络处理器、专用集成电路(asic)、现场可编程门阵列(fpga)、用于存储软件的只读存储器(rom)、随机存取存储器(ram)和非易失性存储装置。还可以包括常规和/或定制的其他硬件。类似地，图中所示出的任何开关仅是概念性的。它们的功能可以通过程序逻辑的操作、通过专用逻辑、通过程序控制和专用逻辑的交互来执行、或者甚至手动地执行，特定技术可由实现者如根据上下文更具体地理解来选择。
[0206]
本领域技术人员应当理解，本文中的任何框图表示体现发明的原理的说明性电路的构思图。类似地，将理解，任何流程图、流向图、状态转换图、伪代码等表示可以基本上在计算机可读介质中表示并且因此由计算机或处理器执行的各种过程，无论是否明确地示出这样的计算机或处理器。
[0207]
尽管上面已经结合具体实施方式阐述了所描述的方法和设备的原理，但是应当理解，该描述仅是通过示例的方式进行的，而不是作为对由所附权利要求书确定的保护的范围的限制。
[0208]
本公开内容还包括以下条款：
[0209]
1.一种用于基于时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布的方法，所述时域磁共振tdmr信号是在根据施加的脉冲序列激励所述样本之后从所述样本发射的，所述方法包括：
[0210]
i)确定tdmr信号模型以近似发射的时域磁共振信号，其中，所述tdmr信号模型被配置成取决于tdmr信号模型参数，所述tdmr信号模型参数包括：
[0211]
所述样本内所述至少一个组织参数的空间分布；以及
[0212]
所述施加的脉冲序列；
[0213]
ii)估计tdmr信号模型参数的初始集合；
[0214]
iii)基于作为输入的所述tdmr信号模型参数的初始集合，确定表示所述发射的tdmr信号与所述tdmr信号模型之间的差异的残差向量；
[0215]
iv)通过下述，基于所述参数的输入集合确定tdmr信号模型参数的更新集合
[0216]
v)计算所述残差向量的梯度；
[0217]
vi)计算所述tdmr信号模型的黑塞矩阵作为稀疏黑塞矩阵，其中，所述黑塞矩阵的稀疏模式基于所述施加的脉冲序列；以及
[0218]
vii)基于所计算的梯度和所计算的黑塞矩阵更新所述tdmr信号模型参数的初始集合；
[0219]
viii)使用作为所述输入的所述tdmr信号模型参数的更新集合，重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)直到所述发射的tdmr信号与所述tdmr信号模型之间的差异小于预定阈值为止或者直到完成预先确定的重复的次数为止，从而获得tdmr信号模型
参数的最终更新集合；以及
[0220]
ix)从所述tdmr信号模型参数的最终更新集合中提取所述至少一个组织参数的空间分布。
[0221]
2.根据条款1的方法，其中，所述施加的脉冲序列包括梯度编码模式和/或射频激励模式，以及其中，所述黑塞矩阵的所述稀疏模式基于所述梯度编码模式和/或所述射频激励模式来确定。
[0222]
3.根据条款1或2的方法，其中，将所述黑塞矩阵计算为所述残差向量的雅可比行列式与所述残差向量的雅可比行列式的厄尔米特转置的乘积。
[0223]
4.根据前述条款中任一项的方法，其中，所述施加的脉冲序列被配置成产生笛卡尔获取、径向获取或螺旋获取中的任一个。
[0224]
5.根据条款4的方法，其中，所述施加的脉冲序列的所述梯度编码模式被配置成产生笛卡尔获取，使得对应的点扩散函数仅在相位编码方向上传播。
[0225]
6.根据前述条款中任一项的方法，其中，所述施加的脉冲序列被配置成产生变化的翻转角。
[0226]
7.根据条款6的方法，其中，所述施加的脉冲序列的所述射频激励模式被配置成产生平滑变化的翻转角，使得对应的点扩散函数在宽度方向上在空间上受到限制。
[0227]
8.根据前述条款中任一项的方法，其中，所述至少一个组织参数包括t1弛豫时间、t2弛豫时间、t2*弛豫时间和质子密度，或者它们的组合中的任一个。
[0228]
9.根据前述条款中任一项的方法，其中，计算所述残差向量的梯度包括将所述残差向量与所述残差向量的雅可比行列式的厄尔米特转置相乘。
[0229]
10.根据前述条款中任一项的方法，其中，所述tdmr信号模型是基于布洛赫的体积信号模型。
[0230]
11.根据前述条款中任一项的方法，其中，所述黑塞矩阵中的所述稀疏模式由所述施加的脉冲序列的所述梯度编码模式来确定。
[0231]
12.根据前述条款中任一项的方法，其中，将所述黑塞矩阵计算为稀疏黑塞矩阵，所述稀疏黑塞矩阵包括多个对角带。
[0232]
13.根据条款12的方法，其中，所述黑塞矩阵中的所述对角带的数量由所述施加的脉冲序列的所述梯度编码模式来确定。
[0233]
14.根据条款12或13的方法，其中，所述多个对角带中的任一个的宽度由所述施加的脉冲序列的所述射频激励模式来确定。
[0234]
15.根据前述条款中任一项的方法，其中，所述稀疏黑塞矩阵对于tdmr信号模型参数的每个集合是恒定的。
[0235]
16.根据前述条款中任一项的方法，其中，更新的步骤vii)包括：
[0236]
通过求解所述稀疏黑塞矩阵与所述tdmr信号模型参数变化相乘等于求负的所述残差函数的梯度的线性方程组，获得tdmr信号模型参数变化的集合，以及
[0237]
将所获得的tdmr信号模型参数变化与所述tdmr信号模型参数被更新的集合相加以获得所述tdmr信号模型参数的更新集合。
[0238]
17.一种用于基于时域磁共振tdmr信号来确定样本内的至少一个组织参数的空间分布的设备，所述时域磁共振tdmr信号是在根据施加的脉冲序列激励所述样本之后从所述
样本发射的，所述设备包括处理器，所述处理器被配置成：
[0239]
i)确定tdmr信号模型以近似发射的时域磁共振信号，其中，所述tdmr信号模型被配置成取决于tdmr信号模型参数，所述tdmr信号模型参数包括：
[0240]
所述样本内所述至少一个组织参数的空间分布；以及
[0241]
所述施加的脉冲序列；
[0242]
ii)估计tdmr信号模型参数的初始集合；
[0243]
iii)基于所述tdmr信号模型参数的初始集合，确定表示所述发射的tdmr信号与所述tdmr信号模型之间的差异的残差向量；
[0244]
iv)通过下述方式确定tdmr信号模型参数的更新集合
[0245]
v)计算所述残差向量的梯度；
[0246]
vi)计算所述tdmr信号模型的黑塞矩阵作为稀疏黑塞矩阵，其中，所述黑塞矩阵的稀疏模式基于所述施加的脉冲序列；以及
[0247]
vii)基于所计算的梯度和所计算的黑塞矩阵更新所述tdmr信号模型参数的初始集合；
[0248]
viii)针对所述tdmr信号模型参数的更新集合，重复步骤iii)、步骤iv)、步骤v)、步骤vi)和步骤vii)直到所述发射的tdmr信号与所述tdmr信号模型之间的差异小于预定阈值为止或者直到完成预先确定的重复的次数为止，由此获得tdmr信号模型参数的最终更新集合；以及
[0249]
ix)从所述tdmr信号模型参数的最终更新集合中提取所述至少一个组织参数的空间分布。
[0250]
18.一种包括计算机可执行指令的计算机程序产品，当所述程序在计算机上运行时，所述计算机可执行指令用于执行所述条款1至16中的任一项的方法。