技术领域
[0001] 本
发明涉及一种确定磁化率梯度图的方法、一种
计算机程序产品以及一种通过
磁共振成像确定对象的磁化率梯度图的设备。
[0002] 当前，尤其是在医学诊断领域中广泛使用了成像MR（磁共振）方法，该方法利用
磁场和核自旋之间的相互作用，以便形成2维或3维图像，因为对于软组织成像，它们在很多方面优于其他成像方法，不需要
电离辐射并且通常是无创的。
背景技术
[0003] 根据一般的MR方法，患者身体或一般对象，必须要布置在强的均匀磁场中，磁场的方向同时定义测量所依据的
坐标系的轴（通常为z轴）。磁场独立于磁场强度针对个体核自旋产生不同的能级，可以通过施加规定
频率（所谓的拉莫尔频率或MR频率）的交变
电磁场（RF场）激励核自旋（自旋共振）。从微观
角度讲，个体核自旋的分布产生整体磁化强度，在磁场垂直于也称为纵轴的z轴延伸的同时，通过施加适当频率的电磁脉冲（RF脉冲）可以使磁化强度偏离平衡状态，使得磁化强度关于z轴进行进动。进动运动描绘出锥形表面，锥形的孔径角称为翻转角。翻转角的大小取决于所施加电磁脉冲的强度和持续时间。对于所谓的90度脉冲的情况，自旋从z轴偏斜到横向平面（翻转角90度）。
[0004] 在终止RF脉冲之后，磁化强度弛豫回初始平衡状态，其中以第一时间常数T1（自旋-晶格或纵向弛豫时间）再次建立z方向的磁化强度，垂直于z方向的方向上的磁化强度以第二时间常数T2（自旋或横向弛豫时间）弛豫。可以利用接收RF线圈探测到磁化强度的变化，在MR设备的检查体积之内布置并定向所述接收RF线圈，从而在垂直于z轴的方向上测量磁化强度的变化。例如，在施加90度脉冲之后，横向磁化强度的衰减伴随着由局部磁场不均匀性诱发的核自旋从具有相同
相位的有序状态到所有相角均匀分布的状态（去相）的过渡。可以利用重新聚焦脉冲，例如180度脉冲补偿去相。这样在接
收线圈中产生回波
信号（自旋回波）。
[0005] 为了在身体中实现空间分辨，在均匀磁场上
叠加沿三个
主轴延伸的线性
磁场梯度，导致自旋共振频率的线性空间相关性。那么接收线圈拾取的信号包含可能与身体中不同
位置相关联的不同频率分量。经由接收线圈获得的信号数据对应于空间频率域，被称为k空间数据。k空间数据通常包括利用不同相位编码采集的多条线。通过收集若干样本对每条线进行数字化。利用傅里叶变换将一组k空间数据变换成MR图像。
发明内容
[0006] 具有与周围偏离的磁化率的对象引起主磁场的局部不均匀性。这适用于金属对象，诸如手术器械、
植入物或其他设备、含
铁物质，如去
氧的血液和自然发生的组织之内的铁沉积、或基于氧化铁的
造影剂或标记细胞。这还适用于被成像身体之内的空白区域，因为在空白区域（填充气体）和周围组织之间有显著的磁化率差异。这种效应的利用是用于以下不同MR成像应用的重要工具：从检测造影剂（例如SPIO）和自然发生的组织边界到
定位设备，例如
导管或植入物。
[0007] 通常通过T2或T2*加权序列进行磁化率
对比度增强的MR成像。利用这些序列，在局部磁场干扰一侧由信号损失生成对比度。在这些已知技术生成的图像中，不能将由于场不均匀性造成的暗图像特征与导致如低自旋
密度或极短的T2时间的信号损失的其他效应导致的特征区分开。
[0008] 已经提出了几种MR序列以将这种暗信号转换成正对比度。不同的方法展示出明亮的正对比度图像，但需要关于场干扰强度的先验知识以便优化正
图像对比度。这一点是利用几种为从梯度回波图像计算正对比度而开发的后期处理方法而避免的。在采集基于梯度回波的图像期间，磁化率梯度局部改变所施加的成像梯度。除了已知的几何失真之外，这还导致k空间中受影响回波的移动。确定针对每个成像
体素（3D体素）的这种移动导致磁化率诱发的梯度图。
[0009] 为了确定磁化率梯度图，当前使用两种
算法。第一种算法是“初始磁化率梯度绘图”（SGM）算法，其依赖于小项（small term）傅里叶变换的计算。例如，在WO 2007/122527A2中详细论述了这种方法。以更高
分辨率产生磁化率梯度图的第二种方法依赖于截短傅里叶变换的使用。这种算法在下文中称为“真实分辨率SGM”，该算法在Dahnke H，Liu W，BowtellR，Frank JA，“High Resolution Positive Contrast via Post-Processing from Conventional 3D Imaging”，Int Soc Magn Reson Med 2008；16 1513中详细描述。它基于Chen N，Oshio K，Panych LP，“Application of k-space energy spectrum analysis to susceptibility field mapping and distortion correction in gradient-echo EPI”，NeuroImage 2006；31：609-622中介绍的k空间能谱分析方法。
[0010] 用于真实分辨率SGM的算法与非常高的计算时间相关联，这阻碍了将它们应用于临床程序。利用上文引用的Chen等人工作中提出的算法，对于64×64×1图像的2D图的计算时间通常为11秒左右，这意味着，例如，对于典型的大脑临床数据集（256×256×150），计算时间大约为40小时。上文引用的Dahnke等人的工作中提出的真实分辨率SGM改善可以将计算时间减少大约20小时，不过对于临床应用而言仍然太长。
[0011] 因此，容易理解，需要一种改进的方法，用于以高
质量和高速确定真实分辨率的磁化率梯度图。因此本发明的目的还有提供一种MR设备，其能够以高质量和高速度确定真实分辨率磁化率梯度图。
[0012] 根据本发明，公开了一种确定磁化率梯度图的方法。该方法包括采集磁共振k空间数据集，所述数据集包括磁共振回
波数据，以及将所采集的k空间数据集转换成图像数据集I（x,y，z），其中对于给定k空间方向，例如x，该方法包括如下步骤：
[0013] a)沿所述给定k空间方向对所述图像数据集执行1D离散傅里叶变换，所述1D傅里叶变换得到一组傅里叶项，所述组的每个傅里叶项都与沿给定k空间方向（x）的k空间指标（index）相关联，
[0014] b)选择在-N/2和N/2之间沿所述k空间方向平稳变化的截短值（kx），其中N是那个方向上所述k空间数据集的尺寸，
[0015] c)通过如下方式递归地生成更新图像数据集Ikx1（x,y，z）：确定所述组傅里叶项中相关联的k空间指标数对应于所述截短值的傅里叶项，并从紧接的先前图像数据集生成1
步骤中生成的图像数据集Ikx-1（x,y，z）减去所确定的傅里叶项，所述相减得到新的更新图
1
像数据集Ikx（x,y,z），
[0016] d)计算新的更新图像数据集的大小Mkx1（x,y，z），
[0017] e)利用增大或减小的截短值（kx）重复步骤b)到d)，
[0018] f)针对给定k空间方向，根据新的更新图像数据集之间的大小变化计算回波移动（m），
[0019] g)根据所述回波移动计算沿所述给定k空间方向的磁化率梯度（Gsu）。
[0020] 优选地，针对所有k空间方向执行上述步骤，其中，之后可以将不同的磁化率梯度图组合成一个大小梯度图。
[0021] 与用于确定磁化率梯度图的真实分辨率SGM方法的
现有技术相比，根据本发明的方法能够显著降低计算时间而不会影响质量，这允许将这种技术应用于临床实践中。例如，在将上述临床数据集用作相应k空间数据集的情况下，可以实现大约10分钟的计算时间。
[0022] 于是，根据本发明，用顺序算法替代现有技术的真实分辨率SGM技术所知的重复完整尺寸3D傅里叶变换的应用，在所述顺序算法中使用递归关系在每次
迭代时更新截短傅里叶变换的计算。所得的算法在大致Nx×Ny×Nz×（Nx+Ny+Nz）次操作中计算3D SGM，从而能够在合理时间内处理大的3D数据集，例如，如在神经学应用中遇到的。
[0023] 通常，SGM技术的输入优选是利用回波序列，例如梯度回波序列，获得的3D，也可能是2D复图像。已知梯度-回波序列通常对B0不均匀性，尤其是对磁化率的局部变化敏感。由于磁化率诱发的B0不均匀性导致的几种效应或伪影是已知的，例如几何失真、急剧信号衰减和回波移动。在SGM技术中，使用从B0不均匀性梯度得到的局部回波移动导出这些梯度的图。
[0024] 由以下方程给出由磁化率诱发的B0不均匀性梯度Gsu导致的局部回波移动m：
[0025]
[0026] 其中G是在三个空间方向（r：读出，e：相位编码，s：切片）之一上施加的成像梯度，τ是施加梯度G的时间间隔（针对读出方向的
停留时间），TE是序列的预期回波时间。
[0027] 从这个方程，显然看出，对于给定的梯度Gsu，回波移动随着回波时间TE线性增大。
[0028] 可以将上文所示方程中的关系反转以获得简化的公式，其将k空间中的回波移动和磁化率梯度Gsu的强度相关：
[0029]
[0030] 其中Gsu是磁化率诱发的B0不均匀性梯度，G是在三个空间方向（r：读出，e：相位编码，s：切片）之一上施加的成像梯度（或梯度增量），m是回波移动，τ是施加梯度G的时间间隔（针对读出方向的停留时间），TE是序列的期望回波时间。
[0031] 真实分辨率SGM方法通过局部计算（即针对每个体素）每个空间方向上的回波移动量而进行。这个步骤是利用k空间谱分析实现的。这种分析的原理在之后的沿X方向计算回波移动的情况下例示。用S（kx,ky，kz）表示所采集的k空间信号，Nx表示沿X方向的
采样矩阵的尺度，真实分辨率SGM算法执行以下步骤：
[0032] -对于kx=-Nx/2...0的每个值，将所有样本S（k,ky，kz），k这种截短的量表示为Skx。然后计算Skx 的3D傅里叶变换的大小，表示为Mkx。
[0033] -针对每个体素计算如下和：M1（x,y，z）=ΣkxMkx1（x,y，z）。
[0034] -针对kx=+Nx/2...0重复同样的程序，其产生M2（x,y,z）。
[0035] -将回波移动计算为归一化差异：（M2（x,y，z）-M1（x,y，z））/M（x,y，z），其中M是S的3D傅里叶变换的大小（无任何截短）。
[0036] 在图1中示出了没有B0不均匀性梯度的情况下，图2示出了在有B0不均匀性梯度的情况下，根据截短程度kx的量Mkx的演变。观察到在第一种情况下针对kx=0的信号下降100，以及在第二种情况下针对kx≠0的信号下降200。那么，差异M2（x,y，z）-M1（x,y，z）与X方向上的回波移动mx成比例（参见图2中的
附图标记202），其允许利用上文所示方程（2）计算磁化率梯度Gsu,x。
[0037] 从关于现有技术真实分辨率SGM的以上描述可以理解，在真实分辨率SGM方法中重复执行所有截短的所采集k空间信号的完整尺寸3D傅里叶变换，因此耗费大量时间。
[0038] 根据本发明，用顺序算法代替所采集k空间信号的截短和后继3D傅里叶变换，在顺序算法中使用递归关系。可以利用以下方程理解这种情况：
[0039]
[0040]
[0041]
[0042] 在这里，I（x,y，z）表示3D复图像，Ikx1（x,y，z）表示截短信号Skx1（k,ky，kz）的3D傅里叶变换，S（k,y，z）表示图像I（x,y,z）的1D傅里叶变换（沿X方向）的第k傅里叶系数，其中k=0对应于“空值”频率项。必须要指出，可以通过信号S（kx,ky，kz）的3D傅里叶变换获得3D复图像I（x,y，z）。
[0043] 上文所示的方程提供了一种有效率地计算M1的计算中涉及的项Mkx1（x,y，z）。为了计算M2，对应的关系写为：
[0044] I0(x,y,z)=I(x,y,z)
[0045]
[0046]
[0047] 必须要指出，优选针对所有体素（x,y，z）执行上述方法。优选地，逐个体素进行Ikn1和傅里叶项之间的相减。还必须要指出，针对每个截短值，仅执行一次方程（3）和（4）中指数的计算和存储，例如，通过明确地评估和存储对应的余弦和正弦项。然后在对所有体素迭代时，可以重新使用这些指数。这些指数仅取决于截短值，例如，在X方向上的情况下，其取决于“x”，而非其他两个指标“y”和“z”。
[0048] 这种算法在大致（Nx+Ny+Nz）·Nx·Ny·Nz次运算中计算3D SGM。作为对比，初始算法（使用重复的3D傅里叶变换）对于同样的计算需要（Nx+Ny+Nz）·（Nx·Ny·Nz）·ln（Nx·Ny·Nz）次运算，不采用递归关系，仅使用1D傅里叶变换，需要Nx·Ny·Nz·（Nx·ln（Nx）+Ny·ln（Ny）+Ny·ln（Ny））次运算。
[0049] 在以上描述中，计算M1和M2，稍后使用其以利用M2和M1之间的归一化差异来计算回波移动。这对应于本发明的
实施例，其中在第一次循环中，在步骤b中，所述截短值从0平稳地增大到沿所述k空间方向的k空间数据集尺寸的一半，在步骤c中得到第一新的更新图像数据集并在步骤d中得到第一大小。此外，在本实施例中，在第二次循环中，在步骤b中，所述截短值从0平稳地减小到沿所述k空间方向的k空间数据集尺寸的负一半，在步骤c中得到第二新的更新图像数据集并在步骤d中得到第二大小。此外，在步骤f中，根据所有第二新的更新图像数据集的所有第二大小之和与所有第一新的更新图像数据集的所有第一大小之和之间的差异计算回波移动。
[0050] 不过，必须要指出这仅仅是实施例。通常，在步骤f中比较针对不同k空间截短值的大小变化并检测对于每个体素发生强度变化的k空间指标就足够了。通过检测强度发生变化的k空间指标，可以确定每个体素的k空间中的移动。
[0051] 利用0和沿给定k空间方向的k空间数据集尺寸一半之间的截短值结合根据第一和第二更新图像数据集所有大小之和之间的差异（即M2-M1）计算回波移动的特征仅仅是鲁棒地确定回波移动。因此，为此目的，分别针对k空间的两个一半，从+k最大到k0（k0=0）并从另一方向-k最大到k0，执行k空间截短。不过，可以结合步骤c中论述的顺序和递归算法使用检测针对每个体素发生强度变化的k空间指标的任何其他技术。
[0052] 根据本发明的另一实施例，在计算第一大小或第二大小的每个单独步骤之后，将所述实际计算的第一大小或第二大小加到迄今为止计算的所有第一大小和第二大小之和，得到新的迄今为止（即到目前为止）计算的所有第一大小或第二大小之和。即在每次更新1 2 1 2
Mkx 和Mkx 之后，计算M 和M 总和。这限制了所提出算法的存储要求。无需在
存储器中保存每个大小计算步骤之后的单个大小，而是简单地将新计算的大小加到先前确定的大小，以数学方式表达为M新=M新+M实际计算的。
[0053] 根据本发明的另一实施例，该方法还包括针对所有剩余的k空间方向重复步骤a到g。例如，这些k空间方向包括读出梯度方向、相位编码方向或切片编码方向之一。因此，通过重复所有步骤，可以提供总的梯度图。
[0054] 根据本发明的另一实施例，3D和1D离散傅里叶变换是快速傅里叶变换（FFT）。因此，进一步加快了算法。
[0055] 根据本发明的另一实施例，在第一次执行步骤c（即递归算法步骤）的情况下，所述方法包括通过对所采集的k空间数据集进行3D傅里叶变换来生成初始图像数据集，其中在这种情况下，在步骤c中，先前的图像数据集生成步骤通过生成所述初始图像数据集来给出。换言之，一开始对所采集的k空间数据集进行3D傅里叶变换，得到初始图像数据集，将其用作在步骤c中递归生成更新图像数据集的起点。在第一次（“更新”）图像生成时，接着从初始图像数据集减去所确定的傅里叶项，可以将更新的图像数据集用于进一步递归步骤。
[0056] 根据本发明的另一实施例，所采集的数据集包括具有不同回波时间的多重回波的磁共振回波数据，其中针对每个回波时间分别确定磁化率梯度图而得到临时（temporary）磁化率梯度图，其中所述方法还包括通过组合所有确定的临时磁化率梯度图计算总磁化率梯度图。例如，多重回波为双重回波。这样允许进一步改善所得磁化率梯度图的质量，不会在采集和
数据处理时间中导致任何可
感知的延迟。
[0057] 为了更深入地理解这种情况，参考方程（2），应用本方程得到了不取决于回波时间TE的特定选择的磁化率梯度的测量。不过，TE的选择对所计算的磁化率梯度图的SNR（
信噪比）有影响。假设用于计算Gsu的梯度回波图像特征在于SNR等于S/σ（S和σ分别是局部信号幅度和噪声标准偏差），可以显示出在计算回波移动m中的噪声标准偏差与σ/S成比例。从方程（2）最终能够得出结论，对于给定的采集参数G和τ，磁化率梯度图的噪声标准偏差与σ/（S·TE）成比例。
[0058] 因此，增大回波时间对SGM的SNR具有直接影响：回波时间加倍使得磁化率梯度图的SNR加倍。不过，通过横向弛豫时间T2*，信号振幅S也受到TE影响。尽管通过增加TE获得的噪声降低通常超过由于指数T2*衰减导致的信号损失，但如果B0不均匀性梯度非常大，可以观察到非常迅速衰减的信号。在这种情况下，在所施加的TE处可能局部不再有可用信号来计算磁化率梯度。
[0059] 为了受益于长回波时间提供的高SNR，同时避免由于非常大B0不均匀性导致的完全信号损失，根据实施例，使用双重回波方法（或通常多重回波方法）。除了在“通常”回波时间TE1处采集的梯度回波图像，额外采集在短回波时间TE0处的梯度回波图像。通常，因为TE1足够长，所以可以在双重回波扫描中这样做而无任何扫描时间损失。示出了TE0和TE1选择的一个范例是，TE1至少比TE0大三倍。
[0060] 从两个梯度回波图像S0和S1，可以通过应用上述SGM技术计算两个磁化率梯度图Gsu,0和Gsu,1。不过，通常可以结合这种多重回波（例如双重回波）方法使用任何SGM技术。通过组合两个图Gsu,0和Gsu,1，获得最佳的磁化率梯度图Gsu。
[0061] 尽管有很多不同方式组合这两幅图，但根据本发明的实施例，可以使用线性组合，选择权重w0和w1，以使最终图中的噪声最小化。这样产生了以下方程：
[0062] Gsu(x,y,z)=w0(x,y,z)·Gsu,0(x,y,z)+w1(x,y,z)·Gsu,1(x,y,z)
[0063]
[0064]
[0065] 权重w0和w1取决于局部信号幅度S0（x,y，z）和S1（x,y，z）。在方程（5）中，已经假设S0和S1的噪声标准偏差相同，如在双重回波方法中那样，如果利用除了回波时间之外相同的序列参数采集两幅图像就是这种情况。否则，应当相应地
修改w0和w1的定义。
[0066] 利用针对权重w0和w1的方程（5），针对每个
像素在信号强度和回波时间之间实现最佳自动平衡。与两幅图Gsu,0和Gsu,1的每幅相比，所得的磁化率梯度图Gsu具有改善的SNR，受到信号缺失的影响比Gsu,1更小。由于非常快衰减信号导致的伪影可以被完全避免的程度，取决于最小可实现的回波时间TE0（通常大约为2-3ms）、梯度回波图像的空间分辨率和B0不均匀性的强度。
[0067] 为了克服这些限制，根据本发明的实施例，可以利用超短回波时间序列，像UTE（例如，在Rahmer J，Boernert P，Groen J，Bos C.3D Radial ultrashort echo-time imaging with T2 adapted sampling.Magn Reson Med 2006;55:1075-1082中描述的）进行第一回波的采集。
[0068] 有可能扩展所提出的方法以采集多重回波图像。如用于T2*绘图那样，可以为此应用多重回波梯度回波序列。或者，可以在增加的回波时间下分别采集图像，其允许在选择回波时间中有更多的灵活性。然后针对每个回波时间TEi计算磁化率梯度图Gsu,i。组合单个的图Gsu,i以产生最佳图Gsu是通过根据下式使方程（5）一般化实现的：
[0069]
[0070]
[0071] 多重回波方法具有以下好处，即回波时间的特定选择几乎对所产生的SGM没有影响。理想地，用户只能决定要花费多少扫描时间，并由系统自动选择回波时间的数量和间距。这种设置确保了SGM结果的最大用户内和用户间可再现性。
[0072] 在另一方面中，本发明涉及一种确定磁化率梯度图的方法，该方法包括采集磁共振k空间数据集，所述数据集包括磁共振回波数据，其中所采集的数据集包括具有不同回波时间的多重回波的磁共振回波数据，其中针对每个回波时间分别确定磁化率梯度图，得到临时磁化率梯度图，其中所述方法还包括通过组合所有确定的临时磁化率梯度图计算总磁化率梯度图，即，可以将这种方法与任何类型的SGM技术一起使用。
[0073] 在另一方面中，本发明涉及一种计算机程序产品，包括计算机可执行指令，以执行上述方法。
[0074] 在另一方面中，本发明涉及一种用于通过磁共振成像确定对象的磁化率梯度图的设备，其中所述设备包括
数据采集系统，用于采集磁共振k空间数据集，所述数据集包括磁共振数据，其中所述设备被布置为将所采集的k空间数据集转换成图像数据集I（x,y，z），其中对于给定k空间方向，例如x，所述设备还被布置为：
[0075] a)沿给定k空间方向对所述图像数据集执行1D离散傅里叶变换，所述1D傅里叶变换得到一组傅里叶项，所述组的每个傅里叶项都与沿给定k空间方向的k空间指标相关联，
[0076] b)选择在-N/2和N/2之间沿所述k空间方向平稳变化的截短值（kx），其中N是那个方向上所述k空间数据集的尺寸。
[0077] c)通过如下方式递归地生成（公式3和4）图像数据集Ikx1（x,y，z）：确定所述组傅里叶项中相关联的k空间指标对应于所述截短值的傅里叶项，并从紧接的先前图像数据1
集生成步骤中生成的图像数据集Ikx-1（x,y，z）减去所确定的傅里叶项，所述相减得到新的
1
更新图像数据集Ikx（x,y，z），
[0078] d)计算所述新的更新图像数据集的大小Mkx，
[0079] e)利用增大或减小的截短值（kx）重复步骤b)到d)，
[0080] f)针对给定k空间方向根据新的更新图像数据集之间的大小变化计算回波移动（m），
[0081] g)根据所述回波移动计算沿所述给定k空间方向的磁化率梯度（Gsu）。
[0082] 优选地，针对所有k空间方向执行上述步骤，其中，之后可以将不同的磁化率梯度图组合成一个大小梯度图。
[0083] 在另一方面中，本发明涉及一种用于通过磁共振成像确定对象的磁化率梯度图的设备，其中所述设备包括数据采集系统，用于采集磁共振k空间数据集，所述数据集包括具有不同回波时间的多重回波的磁共振回波数据，其中所述设备被布置为：
[0084] a)针对每个回波时间分别确定磁化率梯度图，得到临时磁化率梯度图，[0085] b)通过组合所有确定的临时磁化率梯度图计算总磁化率梯度图。
附图说明
[0086] 附图披露了本发明的优选实施例。不过要理解，附图仅仅是为了例示的目的，并不作为定义本发明的限度。在附图中：
[0087] 图1图示了在没有b0不均匀性梯度的情况下根据截短程度的截短信号幅度的演变；
[0088] 图2图示了在有b0不均匀性梯度的情况下根据截短程度的截短信号幅度的演变；
[0089] 图3是图示了根据本发明的方法各步骤的
流程图；
[0090] 图4图示了3T情况下双重回波磁化率梯度绘图采集的各种图像；
[0091] 图5示出了根据本发明的MR扫描器。
具体实施方式
[0092] 图1和2图示了在没有（图1）和有（图2）B0不均匀性梯度时，根据截短程度kx的2
截短信号幅度Mkx 的演变。在图1中，对于kx=0观察到信号下降100，在图2中，对于kx值不等于0，观察到信号下降200。这导致回波移动202，如上文详细所述，允许利用方程2计算磁化率梯度Gsu，x。通常，减去图1和2中所示强度曲线下方的面积并归一化到体素强度导致针对每个体素的回波移动。在所有空间维度中，优选针对所有体素执行这个程序能够产生磁化率梯度图。
[0093] 图3是图示了用于确定磁化率梯度图的个体方法步骤的流程图。该方法开始于步骤300，采集磁共振k空间数据集，其中这个数据集包括磁共振回波数据。在步骤302中，将数据集变换成3D图像（I），继之以步骤304，其中设置期望的k空间方向。
[0094] 之后，在步骤306中，沿着步骤304中设置的k空间方向对图像数据集执行1D离散傅里叶变换，例如，利用快速傅里叶变换。这得到一组傅里叶项，其中所述组的每个傅里叶项都与沿给定k空间方向的相应k-空间指标相关联。在公式3和4中，这些项对应于傅里叶系数S（k,y，z）乘以相应的指数。
[0095] 在步骤308中该方法继续进行，其中kn和h被设置为1。在这里，n指示步骤304中设置的空间方向，即在本实施例中n例如是x、y或z方向。需要指出的是，可以使用任何适当方向作为空间方向，甚至独立于执行MR扫描的任何空间方向。kn表示截短值。
[0096] 由于在结合图3所述的实施例中，k空间的两半都是分别分析的，所以相应指示符h具有1和2二者之间的值。因此，在h=1时，kn的值（对于n=x）介于-Nx/2...0之间，对于h=2，kx的值介于+Nx/2...0之间，Nx为x方向（或通常在期望的k空间方向n上）采样矩阵的尺度。
[0097] 此外，在步骤308中，Ikn1和Ikn2被设置成与I相等。在公式3和4中，这对应于在这些公式中执行的第一步。
[0098] 在执行步骤308之后，该方法继续步骤310，这是递归过程，其中计算截短信号的3D傅里叶变换。详细地，在步骤310中，从对应于kn=0的截短信号的先前计算的3D傅里叶变换中减去对应于kn=1的项。在步骤308中将kn=0时截短信号的这种3D傅里叶变换定义为3D傅里叶变换I，即3D复图像。最终这得到了新的更新3D傅里叶变换。必须要提到的是，步骤310之内可能有额外的内部循环，其存在于所有体素（x,y，z）上的迭代。如上所
1
述，可以逐个像素地进行Ikn 和傅里叶项之间的相减。
[0099] 在步骤312中，确定在步骤310中计算的所述3D傅里叶变换的大小并加到先前步h骤中计算的大小上。由于此时不存在先前步骤，所以假设“先前大小”为0。因此，大小Mknh
（h=1或h=2，kn=0）对应于Ikn 的大小（h=1或h=2且kn=0）。
[0100] 在步骤314中，将kn加1，在步骤316中，检查kn的绝对值是否大于Nx/2（广义地为Nn/2）。假设情况不是这样的，该方法返回到步骤310，在步骤310中对于新的kn，确定更新图像数据集。现在，使用在先前步骤中采集的图像数据集并从该图像数据集减去与新的kn对应的项。这样得到了新的更新图像数据集，在步骤312中将该图像数据集用于更新大小的计算。步骤312中的这种计算包括将在步骤310中确定的更新图像数据集的大小加到在先前循环中在步骤312中计算的大小。
[0101] 在步骤314中，再次增加kn，而在步骤316中，再次检查kn的绝对值是否大于Nx/2（广义地为Nn/2）。
[0102] 在kn的绝对值大于Nn/2时，这表示对于k空间的第一半，己经分析了所有的k空间值。于是，必须要针对第二半执行用于大小计算的同样分析，使得在步骤318中，首先检查是否已经这样做了。由于此刻不是这种情况，该方法继续步骤320，其中将h设置为2，将kn设置为-1。接下来，如上所述，迭代继续在步骤310和316之间进行循环。
[0103] 必须要指出，在替代实施例中，可以在算法的最内层循环中同时进行步骤“h=1”和“h=2”。
[0104] 最后，在成功分析期望k空间方向中k空间中所有截短值时，将该方法引导到步骤322，其中计算针对期望k空间方向n的回波移动mn。上文已经详细论述了这种计算。基本上，计算k空间第二半的所有更新图像数据集的所有大小之和与k空间第一半的所有更新图像数据集的所有大小之和之间的差异并由在步骤300中采集的数据集的3D傅里叶变换大小进行归一化。这样允许利用方程1计算回波移动并进一步计算期望的磁化率梯度。
[0105] 该方法在步骤324中继续，其中检查k空间中的下一个方向是否可用于分析。如果是这种情况，该方法跳回步骤304，其中设置另一k空间方向并继续利用步骤306-324进行循环。在没有任何k空间方向留下用于分析的情况下，该方法在步骤326中结束，在该步骤中最后计算磁化率梯度图。
[0106] 图4示出了3T情况下的双回波磁化率梯度绘图。在步骤a中，在短回波时间TE0处采集第一梯度回波图像S0。图4b中所示的图像示出了对应的磁化率梯度图Gsu,0，图4c中示出了对应的加权图w0，这是根据上文论述的方程5计算的。
[0107] 图4d、e和f中示出了在更长回波时间TE1采集的第二梯度回波图像S1、对应的磁化率梯度图Gsu,1和对应的权重图w1。与用于采集图4a的第一梯度回波图像S0的回波时间TE0相比，第二梯度回波图像S1的回波时间TE1大约大三倍。
[0108] 最终图像4g示出了作为具有权重w0和w1的Gsu,0和Gsu,1线性组合获得的组合磁化率梯度图Gsu,。在0（黑）和1（白）之间缩放权重图。在大磁化率梯度的绘示（箭头）中实现了显著改善。它们对应于第二梯度回波图像中的信号几乎消失但仍然存在于第一回波中的位置。
[0109] 在图5中，将根据本发明的MR成像设备1示为方
框图。设备1包括用于生成静止且基本均匀的主磁场的一组主
磁性线圈2，以及用于叠加额外磁场的三组
梯度线圈3、4和5，额外磁场具有可控强度并具有沿
选定方向的梯度。常规上，将主磁场的方向标记为z方向，将垂直的两个方向标记为x和y方向。由电源11为梯度线圈3、4和5加电。
[0110] 成像设备1还包括RF发射天线6，用于向身体7发射射频（RF）脉冲。天线6耦合到用于生成和调制RF脉冲的
调制器9。还提供了用于接收MR信号的RF接收器，该接收器可以与发射天线6相同或分开。如果发射天线6和接收器在物理上是如图1所示的同一天线，布置发送-接收
开关8以将接收的信号与要发射的脉冲分开。将接收到的MR信号输入到解调器10中。用于梯度线圈3、4和5的发送-接收开关8、调制器9和电源11受到控制系统12控制。控制系统12控制馈送给天线6的RF信号的相位和幅度。
[0111] 控制系统12通常是具有存储器和程序控制的微计算机。解调器10耦合到重建装置14，例如计算机，用于将接收的信号转换成能够在例如视觉显示单元15上看到的图像。
[0112] 此外，MR成像设备1包括一组三个匀场线圈（未示出）。从匀场电源经由分别的匀场通道流经匀场线圈的匀场
电流产生辅助磁场。由控制系统12控制匀场电流的强度以优化主磁场的均匀性。
[0113] 为了实际实施本发明，MR设备1包括用于执行上述方法的程序。例如可以由重建装置14或附着于设备1的另一计算机或
硬件部件执行所述程序。
