CN107942272B - 主磁场不均匀下的分数域磁共振成像方法 - Google Patents

Abstract

本发明公开了一种主磁场不均匀下的分数域磁共振成像方法，属于磁共振成像(MRI)领域。该方法通过测量待成像组织的某一层面上若干个点的磁场强度大小，建立成像区域磁场强度分布的模型，然后建立磁场的多项式模型，按照测量的磁场中是否存在明显的二阶分量将该多项式模型分为二阶多项式模型和高阶多项式模型；再将这两个模型分别代入磁共振的自由感应衰减(FID)信号中，对于二阶模型可以用分数阶傅里叶变换工具进行求解成像物体某一层上的自旋密度函数，对于高阶模型需要通过求解代数方程的方法得到成像物体某一层面上的自旋密度函数，取模得到被测物的灰度图像。该方法能仅通过数学计算而有效的避免主磁场不均匀所造成的伪影、图像模糊等不良影响。

Description

主磁场不均匀下的分数域磁共振成像方法

技术领域

本发明涉及一种主磁场不均匀下的磁共振成像方法，属于核磁共振成像(MRI)技术领域。

背景技术

在磁共振领域，MRI用灰度值把组织自旋密度ρ作为空间坐标的函数表达为ρ(x,y,z)。设一静磁场强度为B₀，在这磁场上加一线性梯度场

因此在样品中沿梯度方向，不同位置就有不同的共振频率，通过拉莫尔公式得到

其中

在实际中，线性梯度场是由梯度线圈g_x,g_y,g_z来产生的，而这些梯度场的是又开关电流控制的。按照上面三维的梯度场成像则成出来的像是一个三维的图像。为了让得到一个二维的灰度图像，假设主磁场B₀方向为z轴方向，而梯度场

即在z轴方向没有梯度变化，在实际中，梯度场为G＝(g_x,g_y,0)。

在样品上加一个90^°射频脉冲时，产生的共振自由感应衰减(FID)信号，该信号在时域上可以表达为

S(t)＝∫ρ(r)ei(Δω+γG·r)tdr,Δω为频偏(1)为了能用傅里叶变换得到自旋密度函数，必须将检测得到的FID信号变成虚拟二维信号的形式，为了得到二维时间的信号，需要用计算机控制高速继电器以快速开关梯度线圈中的电流，在t₁时间段内加梯度场g_x，然后关闭g_x磁场，打开g_y，在t₂时间段内加g_y梯度线圈，这样即可得到FID信号，然后依次重复该实验，在重复中保持t₂时间不变，不断变化t₁时间，这样就可以得到一个虚构的二维信号S(t_x,t_y)。在上面实验的过程中主磁场一直保持加载，因此检测到的二维信号为

展开即为

在主磁场一个常数即为均匀的情况下，若要用检测到的FID信号S(t_x,t_y)得到自旋密度函数ρ(x,y)，只要令

则FID信号的二维形式除了相移之外信号与自旋密度函数为一对傅里叶变换对，因而对二维的FID信号进行二维傅里叶变换即可得到ρ(x,y)，取模值即可得到成像的灰度值。目前MRI技术主要是在主磁场均匀情况下，再外加一个线性梯度场，对检测到的FID信号通过傅里叶变换技术实现对自旋密度或弛豫时间T1，T2表达为灰度值。但是在主磁场不均匀的情况下会使液体吸收线加宽，最终成的图像会产生伪影，分辨率低，信噪比低等不良效果。在通常情况下补偿磁场不均匀多带来的成像的负面效应是通过人工调试，加磁片等方式补偿的，另一方面为了提高磁场的均匀性一般还采用高强度的磁场强度，但是这两种方式都极大的增加了硬件成本，同时也浪费了人力资源。

分数域信号处理是新型的信号处理方式，分数域信号处理包括分数阶傅里叶变换(FRFT)，多项式相位的信号处理以及多种分数阶的信号处理工具。本发明将主磁场不均匀下的磁共振成像问题转化为分数域信号处理的问题，最终解决主磁场不均匀下的磁共振成像问题。

发明内容

本发明的目的为解决磁共振成像中磁场不均匀的问题，提出了一种主磁场不均匀下的磁共振成像方法，该方法能补偿主磁场不均匀所造成的图像伪影，还能增加成像后图像的清晰度。

本发明基于主磁场模型，使得主磁场不均匀的全部信息都参与了成像方法。具体为：首先测量虚拟二维的FID信号，确定带有磁场强度参数的磁共振虚拟二维FID信号的表达式，然后测量并计算图像每个像素点的磁场强度大小，最后通过矩阵运算等数学方法将被成像物体某一层面的自旋密度函数计算出来，该层面的位置和大小由选层编码梯度和射频脉冲的带宽决定。取自选密度函数的模值即为图像的灰度值。

本发明的技术方案分为以下步骤：

步骤一：确定带有磁场强度参数的磁共振虚拟二维FID信号的表达式。在磁场强度不均匀的情况下，对被成像物体的其中一个层面进行二维磁共振成像，得到带有磁场强度参数的虚拟二维FID磁共振信号S(t_x,t_y)，将其表达为如下形式：

其中g_x和g_y分别表示频率编码梯度和相位编码梯度，t_x和t_y分别表示频率编码梯度和相位编码梯度加载的时间，x表示成像区域在磁共振仪中被成像物体层面上的横坐标，y表示成像区域在磁共振仪器中被成像物体层面上的纵坐标，γ表示磁旋比，ρ(x,y)表示被成像物体在坐标(x,y)处的自旋密度大小，B(x,y)表示坐标(x,y)处主磁场强度的大小。

步骤二：采集磁场强度数据，建立多项式主磁场模型。

采集主磁场强度数据，并拟合曲面，选择若干个被成像物体成像层面上的点进行磁场强度测量。根据采集到的磁场数据，建立多项式主磁场模型：

其中，参数a_k，b_k,B₀均为多项式的系数，用数学上的插值方法得到。其中n为多项式的阶数，取法为：若被采集到的主磁场不均匀的数据中有明显的二阶分量，则取2，否则，可以取大于2的任意正整数。

步骤三：求自旋密度函数。

(1)当主磁场不均匀模型的多项式的阶数n为2时；将主磁场模型(6)代入二维FID信号的表达式(5)中，得到：

整理后为

整理虚拟二维FID信号，使其变成标准的二维FRFT。因而在(3)式中定义如下参数：

所述ρ为表示二维分数阶变量的二维参数，α为表示二维分数阶的两个角度的二维参数。

将上述参数代入(7)式，得到

令

则(7)式变为二维FRFT的标准表达形式：

S_α(ρ)＝C_α(ρ)∫∫ρ(x,y)exp[i(xTAx-xTBρ)]dx(11)

其中，(·)T表示矩阵或向量的转置运算，

对得到的二维FRFT(11)进行反变换，得到连续的自旋密度函数的表达式。二维FRFT的反变换表达式为

ρ(x,y)＝∫∫S_α(ρ)C_α(-ρ)exp[i(xTAx+xTBρ)]dρ(12)再由(10)式得到

对(12)式右边进行二重积分的变量代换，再将(8)代入，得到连续的自旋密度函数的计算公式为

其中

J|是变量代换产生的雅可比行列式。可以知道A与B只与t_x,t_y有关。

(2)当主磁场不均匀模型的多项式为高阶时；得到带有磁场强度参数的虚拟二维FID磁共振信号S(t_x,t_y)，将其表达为如下形式：

该表达式没有二阶时的反变换表达式，因此自旋密度函数ρ(x,y)没有像(14)一样的显示表达。

步骤四：求离散的自旋密度函数，取模得被测物的灰度图像。该步骤需要将步骤三的式(14)或式(15)中连续的自选密度函数ρ(x,y)表达式进行离散化，最终的目的是将其求出来。

首先，对于连续的自旋密度的表达式(14)，分别从时间t_x,t_y和空间x,y进行离散化。故而按照t_x,t_y的采样间隔Δt_x,Δt_y采样，并且采样点数分别取N和M，x,y的采样间隔取Δx和Δy，采样点数同样分别取N和M，这样将(14)离散化后的自旋密度函数的表达式为

其中Ω＝Δt_xΔt_yγ²(a₁g_y+g_xb₁+g_xg_y)，为一常数；Δx,Δy是相邻像素点的坐标间隔，x_kl＝(kΔx,lΔy)'是第k行第l个像素点的坐标；t_xn,t_ym分别表示在时间t_x和t_y上的第n个和第m个采样时刻；A_nm，B_nm，ρ_nm，α_xnm，α_ynm分别表示在相应的t_xn,t_ym时刻A，B，ρ，α_x，α_y的值。

上面Δx,Δy的取值原则为：通过磁共振过程中的射频脉冲带宽、频率编码梯度及相位编码梯度的大小确定成像区域的视野(FOV)，即成像方形区域的大小，Δx为将该方形区域x轴方向的边长N等分得到的间隔，Δy为将该方形区域y轴方向的边长M等分得到的间隔。

为了符合傅里叶变换时间间隔与频率间隔之间的关系，(16)式中时间间隔Δt_x和Δt_y分别取值为

令

则得到的二维离散的

即为通过二维分数阶傅里叶变换补偿磁场不均匀性而得到的磁共振自旋密度函数，FID信号S(t_xn,t_ym)为测量的已知值，因此可以得到成像物体一个层面上的自旋密度函数，取模则为待测图像的灰度值。至此，完成了主磁场二阶模型下的磁共振成像。

其次，对于步骤二中自旋密度表达式(15)，也需要同样按照上面的方法进行离散化，取Δx,Δy分别表示x轴和y轴方向离散化的间隔，N，M分别表示x轴和y轴方向离散点的个数(N，M为自然数，任意取值)，则空间离散化得到

x_n,y_m分别表示第n和第m个空间采样点处的坐标值，Δx,Δy和Δt_x,Δt_y的取法与上面提到的原则一样，取时间t_x和t_y的离散化点数也同样分别为N，M个，因此时间离散化后表达式为

式中，1≤k≤N,1≤l≤M k,l分别表示第k个和第l个时间点。离散的虚拟二维FID磁共振信号

总共有N×M个，为已知值；所需求解的离散自旋密度函数ρ(x_n,y_m)总共有N×M个；磁场强度B(x_n,y_m)总共有N×M个。

将主磁场强度的矩阵数据代入步骤二中的时间离散化后表达式中，得到

式中，

B(x_n,y_m)，N，M，Δx，Δy，Δt_x，Δt_y均为已知。等式(20)右边的每一个数的值是待成像物体自旋密度函数的线性组合的形式，总共有NM个方程，本发明通过解方程组的方法求出N×M个自旋密度函数ρ(x_n,y_m)的值，为了将该方程组写成矩阵方程的形式，作如下处理：

令

则(20)式变为

将这个线性方程组写成矩阵Ax＝b的形式。

令

为N×M的矩阵，其中的元素

将该矩阵按列顺序放入大小为NM×1的列向量b中；令ρ为N×M的矩阵，其中的元素ρ_nm＝ρ(x_n,y_m)，x为将ρ矩阵按列顺序放入的大小为NM×1的列向量；令A为NM×NM的矩阵，其中的元素为(21)中的φ(n,m,k,l)，排列形式满足矩阵的乘法规则，使其保证Ax＝b符合(11)式的乘法关系。

求解方程Ax＝b，得到x＝A-¹b，将x按列写成大小为N×M的矩阵，取模便得到被成像物体自旋密度成的灰度图，至此完成了主磁场的任意高阶模型下的磁共振成像方法。

有益效果

本发明通过数据处理的方式进行核磁共振主磁场不均匀补偿，避开了现在核磁共振设备面对主磁场的不均匀而采取贴磁片补偿磁场不均匀等硬件的方法，极大的降低了硬件的成本，由于提出主磁场离散模型的MRI技术完全是自适应的成像方法，只要知道磁体在任意放置情况下在活体组织某一层上的主磁场强度的分布数据，和测量得到的虚拟二维的FID信号，就可以按照该成像方法进行成像，因此这种方法避开了通常情况下为了使磁场变得均匀或者接近均匀常常雇佣专业人员进行磁场均匀调试的麻烦，这样也降低了专业人员配置费用。相对于原始的贴磁片或者人工调试等原始主磁场不均补偿方法，该方法具有自动化程度高，自适应能力强，对磁场不均匀程度要求低等特点。

附图说明

图1示出本发明的实施流程图。

图2示出在主磁场不均匀情况下、二阶模型时，基于传统傅里叶变换的头部切片成像结果图2(a)，以及基于本发明的实施本发明方法的磁共振成像结果图2(b)。

图3示出在主磁场不均匀情况下、高阶模型时，基于传统傅里叶变换的头部切片成像结果图3(a)，以及基于本发明的实施本发明方法的磁共振成像结果图3(b)。

具体实施方式

下面对本发明方法的实现做详细说明：

由于本发明中模拟主磁场的多项式的阶数需要按照主磁场数据中是否有明显的二阶分量而分两种情况，因此在该实施方式中，也分主磁场二阶模型和高阶模型两种情况进行分别说明。

本发明是一种主磁场不均匀下的分数域磁共振成像方法，其流程如图1所示，以像素点为N×M＝128×128的脑部切片的磁共振成像为例，其实现过程如下：

若磁场为二阶模型，则按照如下方法进行成像。

步骤一：在磁场中加入线性梯度场，梯度场强度为g_x＝1Hz/cm,g_y＝12Hz/cm，获取头部切片的K空间数据。接下来将主磁场模型代入公式(1)得到磁共振成像的模型公式

其中磁旋比γ＝42.6Hz/T。

通过变量代换将成像模型转化为标准的二维分数阶傅里叶变换的标准形式：

S_α(ρ)＝C_α(ρ)∫∫ρ(x,y)exp[i(x^TAx-x^TBρ)]dx (23)

步骤二：对于静磁场强度为0.36T的磁共振仪，根据空间分布位置采集磁场数据，对于成像层面的磁场建立二阶多项式模型：

B(x,y)＝p_2xx²+p_2yy²+p_1xx+p_1yy+p₀ (24)

其中B₀＝0.36T，a₂＝-2.432×10^-6,b₂＝-5.946×10^-6,a₁＝2.919×10^-3,b₁＝2.335×10^-2.

步骤三：通过分数阶傅里叶变换理论来求用于成像的自旋密度函数ρ(x,y)，即对(24)式进行反变换。分数阶傅里叶变换的反变换理论为

并且没有交叉项的二维分数阶傅里叶变换具有可分离的核，可以将其分解为两次一维分数阶傅里叶变换，由此很容易得到(23)式的反变换为

ρ(x,y)＝∫∫S_α(ρ)C_α(-ρ)exp[i(x^TAx+x^TBρ)]dρ (25)再次对上式进行数学变换，使之变形为对时间t_x,t_y进行积分：

从而得到了连续自旋密度函数ρ(x,y)的值。

步骤四：对(26)式进行离散化，由于(26)式是在二维时间上对FID信号与一个函数乘积做积分，故而对时间进行采样。设对t_x,t_y分别的采样间隔为Δt_x＝5.761×10^-3,Δt_y＝4.801×10^-4，像素坐标的间隔为Δx＝0.2,Δy＝0.2，则离散化后的自旋密度函数的表达式为

由此所得到的离散的自旋密度，就是通过二维分数阶傅里叶变换补偿磁场不均匀性而得到的磁共振自旋密度函数，取模则得到成像灰度值。

利用上述成像算法，对于静磁场为0.36T的磁共振仪，在主磁场不均匀情况下，基于传统傅里叶变换的头部切片成像结果如图2(a)所示，基于磁场二阶模型下的磁共振成像结果如图2(b)所示。由两幅图的比较，可以明显的看到在主磁场不均匀的情况下，与传统傅里叶变换的成像结果相比，二维FRFT成像结果明显减少了图像畸变和伪影。

若磁场为高阶模型，则按照如下方法进行成像。

确定带有磁场强度参数的磁共振虚拟二维信号的表达式，磁场强度测量和计算图像灰度三步骤组成。本成像方法采用计算机模拟的磁场不均匀性，然后用Matlab中自带的脑切片128×128的矩阵数据进行说明。

首先，写出基于128×128的虚拟二维FID信号的表达式。按照发明内容步骤二的结论，该离散的虚拟二维FID信号的表达式为：

式中，1≤k≤128,1≤l≤128。在该表达式中取相位编码梯度g_y为23T/cm，频率编码梯度g_x为23T/cm，γ为42.6Hz/T像素点的间隔Δx和Δy都为0.36cm，按照(7)式的要求，需取时间采样间隔为Δt_x＝1.39×10^-4s，Δt_y＝1.39×10^-4s，两者的采样点数都取128个点，即N＝M＝128。

再次，得到不均匀的磁场数据，然后进行高阶多项式模拟。为了模拟不均匀的磁场数据，用计算机程序生成128×128的随机数据B，使得成像层面上第nm个像素处的磁场强度B_nm＝0.36T+δ，其中δ为在[-0.10.1]上均匀分布的随机变量。

在本实施方案中，认为这些点B_nm便是高阶多项式模拟然后按照时间间隔和空间间隔采样得到的点，这样便得到了N×M＝128×128的磁场不均匀数据。将主磁场强度的矩阵数据代入离散的虚拟二维的FID信号表达式中，如

其中1≤k≤128,1≤l≤128。

最后，计算图像的灰度值。可以知道，磁共振成像算法就转化成了已知所有

和和主磁场强度数据B(x_n,y_m)的情况下求出ρ(x_n,y_m)的问题，显然等式(10)左边的每一个数的值是待成像组织自旋密度线性组合的形式，可知这样总共有128^{^2}个方程，而未知数也有128^{^2}个，因此可以通过解方程组的形式将其解出来，容易将其写成矩阵Ax＝b的形式。

令

为128×128的矩阵，其中的元素

将该矩阵按列顺序放入大小为NM×1的列向量b中；令ρ为128×128的矩阵，其中的元素ρ_nm＝ρ(x_n,y_m)，x为将ρ矩阵按列顺序放入的大小为NM×1的列向量；令A为128^{^2}×128^{^2}的矩阵，其中的元素为(11)中的φ(n,m,k,l)，排列形式需满足矩阵的乘法规则，使其保证Ax＝b符合(29)式的乘法关系。

求解方程Ax＝b，得到x＝A^-1b，将x写成按列写成大小为128×128的矩阵，取模便得到被成像物体自旋密度成的灰度图像。

用本发明成像的结果如图3(b)所示。为了对比该方法的有效性，本发明将该实验中得到128×128的

按照传统的傅里叶变换方法直接进行FFT得到的结果与本发明中的方法进行对比，可以发现由本方法得到的磁共振图像如图3(b)明显比由FFT得到的图像如图3(a)清晰，而且伪影也很少。

Claims (9)

1.主磁场不均匀下的分数域磁共振成像方法，其特征在于，包括以下步骤：

S1：确定带有磁场强度参数的磁共振虚拟二维FID信号的表达式；

S2：采集磁场强度数据，通过测量成像物体上的磁场数据来拟合磁场的多项式模型，建立多项式主磁场模型：

式(1)中，n为多项式的阶数，若测量的磁场中有明显的二阶分量，则取n＝2；否则，取n>2；参数a_k，b_k,B₀均为多项式的系数，用数学上的插值方法得到，x表示成像区域在磁共振仪中被成像物体层面上的横坐标，y表示成像区域在磁共振仪器中被成像物体层面上的纵坐标，B(x,y)表示坐标(x,y)处主磁场强度的大小；

S3：求自旋密度函数；

S4：求离散的自旋密度函数，取模得被测物的灰度图像。

2.如权利要求1所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，主磁场强度不均匀的情况下，所述步骤S1中的磁共振虚拟二维FID信号的表达式为：

式(2)中g_x和g_y分别表示频率编码梯度和相位编码梯度，t_x和t_y分别表示频率编码梯度和相位编码梯度加载的时间，S(t_x,t_y)表示时间t_x、t_y时的磁共振虚拟二维FID信号，γ表示磁旋比，ρ(x,y)表示被成像的物体在坐标(x,y)处的自旋密度大小。

3.如权利要求1所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，在所述步骤S3中，当n＝2时，连续的自旋密度函数的计算式为：

式(3)中，t_x和t_y分别表示频率编码梯度和相位编码梯度加载的时间，S(t_x,t_y)表示时间t_x、t_y时的磁共振虚拟二维FID信号，γ表示磁旋比，ρ(x,y)表示被成像的物体在坐标(x,y)处的自旋密度大小，且：

其中，g_x和g_y分别表示频率编码梯度和相位编码梯度，ρ表示二维分数阶变量的二维参数，α表示二维分数阶的两个角度α_x和α_y的二维参数，a₁、b₁分别表示主磁场模型B(x,y)中多项式中系数a_k、b_k在k＝1时的值，定义：

其中，a₂、b₂分别表示主磁场模型B(x,y)中多项式中系数a_k、b_k在k＝2时的值，当n>2时，含自旋密度函数的虚拟二维FID磁共振信号S(t_x,t_y)表达为：

4.如权利要求3所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，在所述步骤S4中，求离散的自旋密度函数的方法包括将S3中式(3)或式(4)中连续的自旋密度函数ρ(x,y)表达式进行离散化的步骤，具体为：

当n＝2时，对式(3)进行离散化；

分别从时间t_x,t_y和空间x,y进行离散化；按照t_x,t_y的采样间隔Δt_x,Δt_y采样，并且采样点数分别取N和M，x,y的采样间隔取Δx和Δy，采样点数同样分别取N和M，离散化后的自旋密度函数的表达式为：

式(5)中Ω＝Δt_xΔt_yγ²(a₁g_y+g_xb₁+g_xg_y)，为一常数；Δx,Δy是相邻像素点的坐标间隔，x_kl＝(kΔx,lΔy)是第k行第l个像素点的坐标；t_xn,t_ym分别表示在时间t_x和t_y上的第n个和第m个采样时刻；A_nm，B_nm，ρ_nm，α_xnm，α_ynm分别表示在相应的t_xn,t_ym时刻A，B，ρ，α_x，α_y的值；

当n>2时，对式(4)进行离散化；

取Δx,Δy分别表示x轴和y轴方向离散化的间隔，N，M分别表示x轴和y轴方向离散点的个数，取时间t_x和t_y的离散化点数分别为N，M个，因此时间离散化后表达式为

1≤k≤N,1≤l≤M (6)

x_n,y_m分别表示第n和第m个空间采样点处的坐标值。

5.如权利要求4所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，所述式(6)右边的每一个数的值是待成像物体自旋密度函数的线性组合的形式，总共有N×M个方程，通过解方程组的方法求出N×M个自旋密度函数ρ(x_n,y_m)的值；式(6)变为：

其中，

6.如权利要求5所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，所述式(7)中离散的虚拟二维FID磁共振信号

总共有N×M个，为已知值，所需求解的离散自旋密度函数ρ(x_n,y_m)总共有N×M个，需要参与计算的磁场强度点总共有N×M个。

7.如权利要求5所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，所述式(7)写成Ax＝b的形式，求解得x＝A^-1b，其中，自旋密度函数为x中的元素，将x中的元素按照行排列到N×M大小的矩阵中，然后每个元素取模便是图像的灰度值。

8.如权利要求4～7中任一权利要求所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，所述的的N，M是任意取的，取Δx,Δy的原则为：通过射频脉冲的带宽，频率编码梯度，相位编码梯度的大小确定成像区域的视野，Δx为该方形区域x轴方向的边长N等分的间隔，Δy为该方形区域y轴方向的边长M等分的间隔。

9.如权利要求4～7中任一权利要求所述的主磁场不均匀下的分数域磁共振成像方法，其特征在于，所述的时间间隔Δt_x和Δt_y分别为：

Images