用于磁共振成像系统的梯度涡流校正方法和装置与流程

本发明主要涉及磁共振成像技术领域，尤其涉及一种用于磁共振成像系统的梯度涡流校正方法和装置。

背景技术：

在磁共振成像过程(magnetic resonance imaging,MRI)中，梯度线圈中的电流随时间快速切换，会在周围导体结构中产生涡流，涡流会产生一个空间和时间上都不断变化的磁场，使成像区域内的梯度磁场发生畸变，最终影响MRI图像质量。

为了消除涡流对梯度磁场的影响，一种方法是以产生反向梯度磁场的附加线圈作为自屏蔽线圈，通过减少线圈与磁体之间的相互作用，以抑制涡流的产生。但是这种方法不能完全消除在梯度线圈周围导体中产生的涡电流，因而需要对涡电流进行进一步的校正。

因此期望提供一种快速、准确地拟合涡流场模型，并利用涡流场模型对涡流进行校正的解决方案。

技术实现要素：

针对上述现有技术中存在的问题，本发明提出一种用于磁共振成像系统的梯度涡流校正方法和装置。在涡流场系统参数估计过程中，该方法和装置具有准确、高效、稳健等优点。能够有效解决传统计算方法中目标函数繁琐，初始猜测不准确性以及涡流成分数量难确定等问题。

本发明提供了一种用于磁共振成像系统的梯度涡流校正方法，主要包括以下步骤：建立涡流场模型，所述涡流场模型使用拉普拉斯变换来表示涡流场；采集涡流测量数据并使用所述涡流场模型对所述涡流测量数据进行处理以获取所述涡流场模型的初始参数值；根据所述涡流测量数据和所述初始参数值拟合所述涡流 场模型的参数值，以获得标定涡流场模型；以及根据所述标定涡流场模型对梯度涡流进行校正。

优选地，所述使用所述涡流场模型对涡流测量数据进行处理以获取所述涡流场模型的初始参数值包括：对所述涡流测量数据进行拉普拉斯逆变换运算以获取所述涡流场模型的参数的拉普拉斯谱；提取所述拉普拉斯谱的谱峰信息；以及根据所述拉普拉斯谱的谱峰信息确定所述初始参数值。

优选地，所述拉普拉斯谱为稀疏拉普拉斯谱。

优选地，所述拉普拉斯谱的谱峰信息包括谱峰的位置、边界和/或强度。

优选地，所述根据所述拉普拉斯谱的谱峰信息确定所述初始参数值包括：预设一阈值；确定所述拉普拉斯谱中谱峰强度的绝对值大于和/或等于所述阈值的可信谱峰信息；根据所述可信谱峰信息确定所述初始参数值。

优选地，所述根据所述可信谱峰信息确定所述初始参数值包括：选择与所述可信谱峰的峰值相对应的参数值作为所述初始参数值。

优选地，所述根据所述可信谱峰信息确定所述初始参数值包括：确定所述可信谱峰的边界值；计算位于所述可信谱峰的边界内的峰值的总强度；选择与所述可信谱峰的峰值的总强度相对应的参数值作为初始参数值。

优选地，所述使用所述涡流场模型对所述涡流测量数据进行处理包括：使用正则化对所述涡流测量数据进行拉普拉斯逆变换运算。

优选地，所述根据所述涡流测量数据和所述初始参数值拟合所述涡流场模型的参数值包括：根据所述涡流测量数据和所述初始参数值，采用非线性最小二乘拟合所述涡流场模型的参数。

优选地，所述采用非线性最小二乘拟合所述涡流场模型的参数是通过对所述涡流场模型的参数施加双边线性不等式约束来实现的。

本发明还提供了一种用于磁共振成像系统的梯度涡流校正装置，包括：用于建立涡流场模型的模块，所述涡流场模型使用拉普拉斯变换来表示涡流场；用于采集涡流测量数据并使用所述涡流场模型对所述涡流测量数据进行处理以获取所述涡流场模型的初始参数值的模块；用于根据所述涡流测量数据和所述初始参数值拟合所述涡流场模型的参数值，以获得标定涡流场模型的模块；以及用于根据所述标 定涡流场模型对梯度涡流进行校正的模块。

在涡流场系统参数估计过程中，采用本发明的梯度涡流校正方法和装置具有准确、高效、稳健等优点。能够有效解决传统计算方法中目标函数繁琐，初始参数值估计不准确以及涡流成分数量难确定等问题。

注意到，虽然本说明书描述了以上方面，但本领域技术人员能够理解，在执行本发明的方法和装置时可根据实际需要省略其中一个或多个方面，也落在本发明的范围内。

附图说明

图1是根据本发明提供的梯度涡流校正方法的流程图；

图2示出了不同惩罚强度下所得的拉普拉斯谱图的结果；

图3是根据本发明的梯度涡流拟合曲线图。

具体实施方式

为让本发明的上述目的、特征和优点能更明显易懂，以下结合附图对本发明的具体实施方式作详细说明。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其它不同于在此描述的其它方式来实施，因此本发明不受下面公开的具体实施例的限制。

图1是本发明提供的梯度涡流校正方法的流程图。请参考图1，本发明的梯度涡流校正方法包括以下步骤：

执行步骤S1，建立涡流场模型，所述涡流场模型使用拉普拉斯变换来表示涡流场。

为了提高运算效率和准确性，本实施例中使用拉普拉斯变换来表示涡流场，如下所示：

其中，k＝1/τ，d_m为第m个涡流采样数据，γ为氢核的旋磁比，G_x为测试梯度大小，，T_A为测试梯度脉冲下降沿持续时间，T_m为第m个涡流采样数据的延迟时间，m表示测量时间点，α为幅度参数，τ为时间参数。

线性部分x(k)为：

非线性部分k为：k＝1/τ。

式(1)的形式表明，d_m是x(k)的拉普拉斯变换，故求解x(k)只需对d_m进行拉普拉斯逆变换(Inverse Laplace Transformation，ILT)即可。

虽然，上述以式(1)为例对本发明中涡流场模型进行了说明，但本领域技术人员将知晓其他使用拉普拉斯变换来表示涡流场的涡流场模型也是在本发明的保护范围内，只要涡流场模型中线性部分与非线性部分相互独立且符合拉普拉斯变换形式即可达到优化计算的技术效果。

执行步骤S2，采集涡流测量数据并使用所述涡流场模型对所述涡流测量数据进行处理以获取所述涡流场模型的初始参数值。

在本实施例中，采用以下方法获取涡流场模型的初始参数值，包括：

(1)对所述涡流测量数据进行拉普拉斯逆变换运算以获取所述涡流场模型的参数的拉普拉斯谱；

(2)提取所述拉普拉斯谱的谱峰信息；

(3)根据所述拉普拉斯谱的谱峰信息确定所述初始参数值。

下面以式(1)为例进行说明，与离散傅里叶变换(Discrete Fourier Transform，DFT)相似，拉普拉斯逆变换运算在数值上可通过求解一个线性反演问题来实现。首先需要将反演基离散化，即取离散的k值覆盖其可能的取值范围。由于k可能跨越多个数量级，可采用对数等间隔分布的k值。为避免反演得到的拉普拉斯谱分辨率过低，一般取几十到上百个离散的k值。假设n(n＝1,2,…,N)表示离散k值的序号，则式(1)可改写为：

上式若以线性代数的形式表达则为：

d＝Cx (3)

其中，d＝[d1,d2,…,dm]^T，x＝[x1,x2,…,xn]^T为列向量，上标T表示转置， 而C是一个m行n列的系数矩阵，此处亦称拉普拉斯变换算子，其元素定义为：

由于式(3)是欠定的(m<n)，因此直接采用x＝C\d或x＝C⁺d(上标+表示Moore-Penrose广义逆)求解式(2)的结果是病态的。这类问题可采用正则化来求解，即将以上所述的k的离散取值代入下式来求取x，从而得到x关于k的分布函数：

其中，R(x)是惩罚项，用于附加特定的约束条件，如平滑性、稀疏性等。为了得到实践上能够用来修正梯度场的涡流场系统参数，有必要要求涡流成分数目尽可能有限。此时，应当考虑采用稀疏性的约束条件，即l₁-范数惩罚：

R(x)＝λ‖x‖₁ (6)

上式中，λ是拉格朗日乘数，或称正则化参数，用于调节惩罚力度，在回归效果和稀疏度之间寻求平衡。其适当的取值可通过某种后验策略实现，如L-曲线法或对谱峰稀疏度的限定等。此处的l₁-范数惩罚具有强制稀疏特性，可用于重建x关于k的稀疏拉普拉斯谱。举例而言，可使用“Sparse MRI:The Application of Compressed Sensing for Rapid MR Imaging,M.Lustig,D.Donoho,和J.M.Pauly,Magn.Reson.Med.,2007,58,1182”中描述的方法来通过l₁-范数惩罚获得稀疏拉普拉斯谱。

图2示出了不同惩罚强度下所得的拉普拉斯谱图的结果。

通过以上步骤，可以得到x关于k的分布函数，进而可得到α关于τ的分布函数。这样做的好处是，与DFT相似，可排除人为干预，在无需猜测的情况下得到成分及其相应含量的谱图。

根据以上步骤得到x关于k的稀疏拉普拉斯谱后，提取所述拉普拉斯谱中的谱峰信息，包括谱峰的位置、强度等信息，然后根据提取的谱峰信息确定涡流场模型中的初始参数值。理论上可以用所有非零的x结合其相应的k得到涡流场系统参数α和τ，然后用于梯度涡流的补偿。但实际上由于原始数据包含噪声，计算得到的稀疏拉普拉斯谱可能带有成分含量较小的伪成分。另一方面，由于谱图的分辨率有限，离散取值的反演基不可能恰好是涡流的特征成分。因此，需要对稀疏拉普拉斯谱进行处理，从中获取关键的信息。

首先，提取稀疏拉普拉斯谱中的可信谱峰。具体为，预设一阈值，确定所述拉普拉斯谱中谱峰强度的绝对值大于和/或等于预设阈值的可信谱峰信息，这样可以移除部分由于噪声影响而计算出的伪成分，得到L个可信谱峰，最后根据获取的可信谱峰信息确定涡流场模型中的初始参数值。

在本实施例中，根据可信谱峰信息确定涡流场模型中参数x和k的初始值，具体为，选择与所述可信谱峰的峰值相对应的x值，以及x值所对应的k值作为所述初始参数值。在本发明的其他实施例中，采用以下方式确定涡流场模型中参数x和k的初始值。

首先确定可信谱峰的边界值，读取每个可信谱峰(如第l个峰)的最大x_l值所对应的k_l值，记为k0_l。从k0_l开始，分别向零和正无穷两个方向进行检索。此时谱峰强度的绝对值逐渐减小，当强度绝对值小于阈值，或强度绝对值突然增大(表明谱峰有重叠)时，停止检索，将此时的k值记为klb_l(零方向)或kub_l(正无穷方向)。然后计算位于可信谱峰的边界内的峰值强度和，例如将klb_l与kub_l之间的x值进行积分，记为x0_l，即为该谱峰的总强度。通过上述操作，可将每个可信谱峰的位置(k0)、边界(klb与kub)以及相应的总强度等信息提取出来。

当谱峰信息提取完毕，可利用各个可信谱峰的总强度[x0_l]和[k0_l]作为初始猜测(即，初始参数值)，对式(1)进行非线性最小二乘拟合。

通过以上步骤，涡流成分数难以确定以及未知量的合理初始猜测难以给出的问题可以得到较好的解决。但由于未知量之间的跨度较大，尤其对于不同的涡流成分来说，其k值可能相差多个数量级。这种情况下实施无约束优化，原本较小的未知量在迭代过程中比较容易偏离其可能的正确取值范围。

通过申请人的大量数值试验，观察到在x关于k的稀疏拉普拉斯谱中，每个谱峰的两个边界klb_l和kub_l很好地界定了k_l的取值范围，而相应总强度x0_l与x_l的真实值比较接近，可取xlb_l＝μ·x0_l以及xub_l＝ν·x0_l作为x_l的边界。关于系数μ和ν，通过申请人的大量实验得到，如果x0_l>0，则0<μ<1，1<ν<+∞；如果x0_l<0，则1<μ<+∞，0<ν<1，可得出比较理想的拟合结果。一般可取μ＝0.5和ν＝2(x0_l>0)或μ＝2和ν＝0.5(x0_l<0)来生成x_l的上下界。这样，所有待求解的未知量都有明确的取值范围，使得初始猜测既是局部极优，同时也是全局极优，使优化过程更为稳健。

基于式(2)，涡流系统参数的估计可通过求解如下问题实现：

其中，k_l＝1/τ_l，d_m为第m个涡流采样数据，γ为氢核的旋磁比，G_x为测试梯度大小，，T_A为测试梯度脉冲下降沿持续时间，T_m为第m个涡流采样数据的延迟时间，m(m＝1,2,3,…,M)表示测量时间点，l(l＝1,2,3,…,L)表示涡流成分的序号，α_l为幅度参数，τ_l为时间参数。

关于式(7)所描述的非线性最小二乘问题，可利用通过以上步骤得到的谱峰信息对涡流场系统参数施加双边线性不等式约束来进一步提高算法的准确性。双边线性不等式约束通过内点法实现，即在原目标函数的基础上加入对数障碍函数(logarithm barrier function)作为惩罚项。此时的目标函数为：

其中，惩罚项P(x,k)为：

β为拉格朗日乘数，其取值是通过最小残差来确定的。举例而言，该过程可通过“Numerical Recipes,W.H.Press,B.P.Flannery,S.A.Teukolsky,W.T.Vetterling,Cambridge University Press,New York,1988”中描述的数值方法实现。

经过以上所述的约束优化，可以求解出所有x和k值。计算k的倒数可以得到梯度涡流场系统参数τ。如上所述，使用x和τ可计算出涡流场系统参数α。由此，梯度涡流场系统参数α和τ被估计出来，以用于梯度涡流补偿。

图3示出了根据本发明的涡流拟合曲线图，其中参考曲线是实际测得的涡流场，拟合曲线是通过本发明的方法拟合所述涡流场模型得到的涡流场。从图3可看出根据本发明的方法能够准确地拟合实际涡流场。

以上描述了本发明的多个方面，但本领域技术人员能够理解，在本发明的执行中可根据实际需要省略其中一个或多个方面，也落在在本发明的范围内。举例而言，在获取拉普拉斯谱之后也可使用其他方法提取初始参数值，例如按照 现有技术中的人工初始猜测，从而在设计复杂度与技术效果之间达到折衷，这也是本发明所构想到的。

各个方面的其他组合方式也是本领域技术人员在本说明书的基础上能够构想到的。

通过本发明的方法，解决了现有计算方法中目标函数繁琐，初始猜测准确性以及涡流成分数量的确定等问题。

应理解，所公开的方法中各步骤的具体次序或阶层是示例性过程的解说。基于设计偏好，应该理解，这些方法中各步骤可以重新编排次序或组合。

本领域技术人员将领会，结合本文中所公开的实施例描述的各种解说性逻辑块、模块、电路、和算法步骤可被实现为电子硬件、计算机软件、或两者的组合。为清楚地解说硬件与软件的这一可互换性，各种解说性组件、块、模块、电路、和步骤在上面是以其功能性的形式作一般化描述的。此类功能性是被实现为硬件还是软件取决于具体应用和施加于整体系统的设计约束。技术人员对于每种特定应用可用不同的方式来实现所描述的功能性，但这样的实现决策不应被解读成导致脱离了本发明的范围。

虽然本发明已参照当前的具体实施例来描述，但是本技术领域中的普通技术人员应当认识到，以上的实施例仅是用来说明本发明，在没有脱离本发明精神的情况下还可作出各种等效的变化或替换，因此，只要在本发明的实质精神范围内对上述实施例的变化、变型都将落在本申请的权利要求书的范围内。