一种在强不均匀磁场下磁共振图像扭曲的矫正方法

摘要

一种在强不均匀磁场下磁共振图像扭曲的矫正方法，涉及磁共振图像。模型上，通过二维傅里叶正交变化将不均匀场完备地描述为四维矩阵的形式，再将四维矩阵降维成二维的大型矩阵。方法上主要是迭代进行两个模块，一是给定不均匀场，建立并求解了基于正交基变换和压缩感知的l1范数最优化模型；二是给定磁矩密度分布图像，使用解缠绕算法和布谷鸟最优化搜索算法，实现了拟合不均匀场。对强不均匀场和局部不均匀场下磁共振快速序列成像扭曲的矫正有显著的作用。克服磁共振成像快速序列对强不均匀场的敏感和矫正不均匀场造成的图像扭曲。

说明书

技术领域

本发明涉及磁共振图像，尤其是涉及一种在强不均匀磁场下磁共振图像扭曲的矫正方法。

背景技术

磁共振成像具有良好临床安全性，广泛应用在医学诊断上。磁共振采集到的数据的空间称为k空间，磁共振图像的空间称为图像空间。自旋回波EPI序列是一种快速序列，只需要激发一个RF脉冲就可以采集整个k空间数据，能够有效地克服磁共振成像速度慢的缺点。传统的磁共振图像重建方法是对采集到的二维数据进行傅里叶变化。然而，EPI成像对不均匀场很敏感，傅里叶变化重建的图像在不均匀场区域通常会存在扭曲。

Zibetti(ZibettiMVW,DePierroAR.ANewDistortionModelforStrongInhomogeneityProblemsinEcho-PlanarMRI[J].IEEETransactionsonMedicalImaging,2009,28(11):1736-1753.)提出了一种强不均匀场问题的模型，通过点扩散函数和离散分辨率的思想，可以在已知不均匀场的情况下重建图像。但是他通过调制离散分辨率水平来简化和加速算法的同时，也损失了许多细节信号，导致重建图像残留比较严重的伪影。

Donoho(DonohoDL.Compressedsensing[J].IEEETransactionsonInformationTheory,2006,52(4):1289-1306.)提出了压缩感知模型，若图像在变换域内是稀疏的，并且采样方式和稀疏变化的相干性足够小，就能够通过低于奈奎斯特采样率的数据来重建图像。Lustig(LustigM,DonohoD,PaulyJM.SparseMRI:TheapplicationofcompressedsensingforrapidMRimaging[J].MagneticResonanceinMedicine,2007,58(6):1182-1195.)将压缩感知应用磁共振成像中。目前压缩感知在磁共振图像重建的主要应用包括超分辨重建，欠采样重建和去噪。

在求取复数相位时候，求得的值总是限定在[-π,π]之间。事实上，超过这个区间的相位值会被取余从而被限定在区间内，由此产生了相位缠绕现象。Maier(MaierF,FuentesD,WeinbergJS,etal.RobustphaseunwrappingforMRtemperatureimagingusingamagnitude-sortedlist,multi-clusteringalgorithm[J].MagneticResonanceinMedicine,2015,73(4):1662-1668.)提出了一种相位解缠绕方法，能够有效地从限定在[-π,π]的相位恢复到解缠绕的平滑的相位值。

2009年，Yang(YangX-S,DebS,EngineeringOptimisationbyCuckooSearch[J].MathematicalModelling&NumericalOptimisation,2010,1(4):330-343)提出了一种新的启发式智能算法算法——布谷鸟搜索。其主要思想是通过模拟布谷鸟寄生育雏的行为以及Lévy飞行机制来有效地解决最优化问题。

共轭梯度法是一种技术成熟的求解线性方程组的方法。ShewchukJR.(ShewchukJR,Anintroductiontotheconjugategradientmethodwithouttheagonizingpain[M].Carnegie-MellonUniversity.DepartmentofComputerScience.1994.)详细介绍了共轭梯度法的原理和实现方法。

发明内容

本发明的目的在于提供一种在强不均匀磁场下磁共振图像扭曲的矫正方法。

本发明包括以下步骤：

1)计算不均匀场演化时间：t(k_x,k_y)表示自从施加RF脉冲后到k空间采样位置为(k_x,k_y)的不均匀场演化时间时间，表示为：

其中t_c为施加RF脉冲后到数据采样前的时间差值，EPI序列中t_c数值很小，可以忽略不计。k_x,k_y的为k空间笛卡尔系正交的离散的坐标变量。k空间的正交两个坐标方向分别称为行坐标轴方向(k_x方向)和列坐标轴方向(k_y方向)。对应的,x,y为图像空间中正交的离散的两个坐标变量，图像空间正交的两个方向也称为相位编码方向(y坐标轴方向)和频率编码方向(x坐标轴方向)。n(k_y)表示列坐标轴方向采集点的坐标值为k_y对应着的列数编号，T_ro为采集一列k空间数据的时间，τ_pe为相位编码梯度单个周期内的施加时长。G_ro为频率编码梯度，γ为旋磁比。

2)正交基变化系数矩阵C的构造函数：磁矩的演化相位为：

其中B_inh(x,y)为不均匀场图，系数矩阵C(k_x,k_y,h_x,h_y)通过对部分二维反傅里叶变化获得，即

其中h_x,h_y为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N。那么k_x,k_y,h_x,h_y,x,y离散点的个数也都分别等于N。令将四维的矩阵C(k_x,k_y,h_x,h_y)的转变为二维矩阵其中L_x为频率编码方向上的视野。n(k_y)和n(h_y)分别表示k空间列坐标轴方向上采集位置为k_y和h_y对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(B_inh)，其中为的矩阵形式，B_inh为不均匀场B_inh(x,y)的矩阵形式。函数的输入为不均匀场B_inh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。

3)矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CFρ。其中为采集到的k空间数据，为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的向量。此时C＝C(B_aut)表示是真实存在的而未知的不均匀场B_aut对应的系数矩阵。

4)初始化拟合的不均匀场为B_inh＝0，初始化重建图像ρ＝F^Hs。此时的ρ是受不均匀场影响而扭曲的图像。其中F^H的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。

5)更新拟合不均匀场：使用解缠绕算法计算图像ρ的解缠绕相位b_inh(参见文献：MaierF,FuentesD,WeinbergJS,etal.RobustphaseunwrappingforMRtemperatureimagingusingamagnitude-sortedlist,multi-clusteringalgorithm[J].MagneticResonanceinMedicine,2015,73(4):1662-1668.)，并使用布谷鸟算法求解系数p(参见文献：YangX-S,DebS,EngineeringOptimisationbyCuckooSearch[J].MathematicalModelling&NumericalOptimisation,2010,1(4):330-343)。布谷鸟的最优化函数设置为其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(B_inh+pb_inh)表示用矩阵C的构造函数和不均匀场B_inh+pb_inh来构造C矩阵。||·||₂表示2范数。设置中间变量B_T＝B_inh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位b_inh和系数p的乘积，表示为：

B_inh＝B_T+pb_inh(4)

若则进行步骤10)，否则转至步骤6)；

6)用已拟合的不均匀场B_inh来构造系数矩阵C＝C(B_inh)。建立基于正交基变换和压缩感知的最优化模型其中α＝Ψ_ρ表示为ρ经过小波变换后的系数。

7)设置模型参数β,λ，β取值范围在[5,1.5]，λ取值范围与β相关，为[20β,100β]。ρ⁽ⁱ⁾,ν⁽ⁱ⁾,α⁽ⁱ⁾表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν⁽⁰⁾＝α⁽⁰⁾＝0，ρ⁽⁰⁾＝ρ，构造矩阵

8)计算中间变量其中Ψ^H,C^H分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x(参见文献：ShewchukJR,Anintroductiontotheconjugategradientmethodwithouttheagonizingpain[M].Carnegie-MellonUniversity.DepartmentofComputerScience.1994.)。通过以下三个式子更新ρ⁽ⁱ⁺¹⁾，α⁽ⁱ⁺¹⁾，ν⁽ⁱ⁺¹⁾；

ρ⁽ⁱ⁺¹⁾＝F^Hx(5)

$>{\alpha }^{(i+1)}=max\left\{\right|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}|-\frac{1}{\beta },0\}\frac{{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}}{\left|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}\right|}---\left(6\right)$>

ν⁽ⁱ⁺¹⁾＝ν⁽ⁱ⁾+(Ψρ⁽ⁱ⁺¹⁾-α⁽ⁱ⁺¹⁾)(7)

若迭代次数i≥100或者则ρ＝ρ^(k+1)，并进行步骤9)，否则i值增加1，重复进行步骤8)；

9)若则进行步骤10)，否则转至步骤5)。

10)最终输出矫正扭曲后的图像ρ。

本发明的有益效果是：不均匀场场图采用拟合的方式获得的，无需采集参考场，避免额外的采集时间。由于傅里叶基正交完备，给定的不均匀场信息可以被完整的表达在系数矩阵C中，从而获得良好的矫正结果。同时，不均匀场的拟合方式对局部信息有良好的拟合效果，对产生的较强的局部的扭曲有明显矫正效果。

附图说明

图1是EPI序列图。PE表示相位编码方向，FE表示频率编码方向。SS表示选层方向。在传统的EPI序列基础上增加不均匀场示意图。不均匀场B_inh在整个采样过程可视为是不变的。

图2是采集k空间数据的轨迹图。

图3是未矫正的磁共振图像幅值图和相位图。(a)(b)是柠檬图像的幅值图和相位图；(c)(d)是鼠脑图像的幅值图和相位图。

图4是矫正后磁共振图像的幅值图和相位图。(a)(b)是柠檬图像的幅值图和相位图；(c)(d)是鼠脑图像的幅值图和相位图。

具体实施方式

本发明包括以下步骤：

1)计算不均匀场演化时间：t(k_x,k_y)表示自从施加RF脉冲后到k空间采样位置为(k_x,k_y)的不均匀场演化时间时间，表示为

其中t_c为施加RF脉冲后到数据采样前的时间差值，EPI序列中t_c数值很小，可以忽略不计。k_x,k_y的为k空间笛卡尔系正交的离散的坐标变量。k空间的正交两个坐标方向分别称为行坐标轴方向(k_x方向)和列坐标轴方向(k_y方向)。对应的,x,y为图像空间中正交的离散的两个坐标变量，图像空间正交的两个方向也称为相位编码方向(y坐标轴方向)和频率编码方向(x坐标轴方向)。n(k_y)表示列坐标轴方向采集点的坐标值为k_y对应着的列数编号，T_ro为采集一列k空间数据的时间，τ_pe为相位编码梯度单个周期内的施加时长。G_ro为频率编码梯度，γ为旋磁比。

采集k空间数据的轨迹图参见图2。

2)正交基变化系数矩阵C的构造函数：磁矩的演化相位为

其中B_inh(x,y)为不均匀场图，系数矩阵C(k_x,k_y,h_x,h_y)通过对部分二维反傅里叶变化获得，即

其中h_x,h_y为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N。那么k_x,k_y,h_x,h_y,x,y离散点的个数也都分别等于N。令将四维的矩阵C(k_x,k_y,h_x,h_y)的转变为二维矩阵其中L_x为频率编码方向上的视野。n(k_y)和n(h_y)分别表示k空间列坐标轴方向上采集位置为k_y和h_y对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(B_inh)，其中为的矩阵形式，B_inh为不均匀场B_inh(x,y)的矩阵形式。函数的输入为不均匀场B_inh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。

3)矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CF_ρ。其中为采集到的k空间数据，为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的向量。此时C＝C(B_aut)表示是真实存在的而未知的不均匀场B_aut对应的系数矩阵。

4)初始化拟合的不均匀场为B_inh＝0，初始化重建图像ρ＝F^Hs。此时的ρ是受不均匀场影响而扭曲的图像。其中F^H的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。

5)更新拟合不均匀场：使用解缠绕算法计算图像ρ的解缠绕相位b_inh(MaierF,FuentesD,WeinbergJS,etal.RobustphaseunwrappingforMRtemperatureimagingusingamagnitude-sortedlist,multi-clusteringalgorithm[J].MagneticResonanceinMedicine,2015,73(4):1662-1668.)，并使用布谷鸟算法求解系数p(YangX-S,DebS,EngineeringOptimisationbyCuckooSearch[J].MathematicalModelling&NumericalOptimisation,2010,1(4):330-343)。布谷鸟的最优化函数设置为其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(B_inh+pb_inh)表示用矩阵C的构造函数和不均匀场B_inh+pb_inh来构造C矩阵。||·||₂表示2范数。设置中间变量B_T＝B_inh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位b_inh和系数p的乘积，表示为

B_inh＝B_T+pb_inh(4)

若则进行步骤10，否则转至步骤6

6)用已拟合的不均匀场B_inh来构造系数矩阵C＝C(B_inh)。建立基于正交基变换和压缩感知的最优化模型其中α＝Ψ_ρ表示为ρ经过小波变换后的系数。

7)设置模型参数β,λ，β取值范围在[5,1.5]，λ取值范围与β相关，为[20β,100β]。ρ⁽ⁱ⁾,ν⁽ⁱ⁾,α⁽ⁱ⁾表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν⁽⁰⁾＝α⁽⁰⁾＝0，ρ⁽⁰⁾＝ρ，构造矩阵

8)计算中间变量其中Ψ^H,C^H分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x(ShewchukJR,Anintroductiontotheconjugategradientmethodwithouttheagonizingpain[M].Carnegie-MellonUniversity.DepartmentofComputerScience.1994.)。通过以下三个式子更新ρ⁽ⁱ⁺¹⁾，α⁽ⁱ⁺¹⁾，ν⁽ⁱ⁺¹⁾

ρ⁽ⁱ⁺¹⁾＝F^Hx(5)

$>{\alpha }^{(i+1)}=max\left\{\right|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}|-\frac{1}{\beta },0\}\frac{{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}}{\left|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}\right|}---\left(6\right)$>

ν⁽ⁱ⁺¹⁾＝ν⁽ⁱ⁾+(Ψρ⁽ⁱ⁺¹⁾-α⁽ⁱ⁺¹⁾)(7)

若迭代次数i≥100或者则ρ＝ρ^(k+1)，并进行步骤9，否则i值增加1，重复进行步骤8。

9)若进行步骤10，否则转至步骤5。

10)最终输出矫正扭曲后的图像ρ。

以下实施例使用7T(299.8MHz)/160-mm孔径MRI系统对柠檬和鼠脑进行成像。使用EPI序列实验来采集数据，序列图如图1所示，实验参数记录见表1。

表1

实验参数柠檬实验鼠脑实验频率编码方向视野L_x(mm)8045频率编码梯度G_ro(T·m^-1)0.07340.1305相位编码方向视野L_y(mm)8045相位编码梯度G_pe(T·m^-1)0.02450.0326采样点数N_x,N_y12864

第一步，计算不均匀场演化时间：t_c为施加RF脉冲后到数据采样前的时间差值，EPI序列中t_c数值很小，可以忽略不计,所以t_c设置为0。两个方向数值分辨率N＝128，旋磁比设置为水的旋磁比γ＝2π×42.576×10⁶。计算其他参数的值其中T_ro为采集一列k空间数据的时间，τ_pe为相位编码梯度单个周期内的施加时长。k_x,k_y的为k空间笛卡尔系正交离散的坐标变量。k_x,k_y取值范围是$>k_x\in [\frac{-64}{L_x},\frac{-63}{L_x}...0,...\frac{62}{L_x},\frac{63}{L_x}],k_y\in [\frac{-64}{L_y},\frac{-63}{L_y}...0,...\frac{62}{L_y},\frac{63}{L_y}].$>n(k_y)表示列坐标轴方向采集坐标值为k_y对应着的列数编号，代入t_c,T_ro,τ_pen(k_y),k_x,G_ro,γ的数值得到从施加RF脉冲后到k空间采样位置为(k_x,k_y)的时间，即不均匀场的演化时间t(k_x,k_y)，表示为：

第二步，正交基变化系数矩阵C的构造函数：磁矩的演化相位为

B_inh(x,y)为不均匀场图，系数矩阵C(k_x,k_y,h_x,h_y)通过对的部分二维反傅里叶变化获得，即

其中h_x,h_y为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N＝128。那么k_x,k_y,h_x,h_y,x,y离散点的个数也都分别等于128。令将C(k_x,k_y,h_x,h_y)的转变为二维矩阵其中n(k_y)和n(h_y)分别表示相位编码方向上采集位置为k_y和h_y对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵C的过程的写成函数的形式，写作C＝C(B_inh)，其中的矩阵形式，B_inh为不均匀场B_inh(x,y)的矩阵形式。函数的输入为不均匀场B_inh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。

第三步，矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CFρ。其中为采集到的k空间数据，为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的列向量。此时C＝C(B_aut)表示是真实存在的而未知的不均匀场B_aut对应的系数矩阵。

第四步，初始化拟合的不均匀场为B_inh＝0，重建图像ρ＝F^Hs。其中F^H的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。此时的ρ是受不均匀场影响扭曲的图像，如图3所示。

第五步，更新拟合不均匀场：使用解缠绕算法获得图像ρ的解缠绕相位b_inh，并使用布谷鸟算法求解系数p。布谷鸟的最优化函数设置为其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(B_inh+pb_inh)表示用矩阵C的构造函数和不均匀场B_inh+pb_inh来构造系数矩阵，||·||₂表示2范数。设置中间变量B_T＝B_inh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位b_inh和系数p的乘积，表示为

B_inh＝B_T+pb_inh(4)

若则进行第十步，否则转至第六步

第六步，用已拟合的不均匀场B_inh来构造系数矩阵C＝C(B_inh)。建立基于正交基变换和压缩感知的最优化模型其中α＝Ψ_ρ表示为ρ经过小波变换后的系数。

第七步，设置模型参数β＝8,λ＝400。ρ⁽ⁱ⁾,ν⁽ⁱ⁾,α⁽ⁱ⁾表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν⁽⁰⁾＝α⁽⁰⁾＝0，ρ⁽⁰⁾＝ρ，构造矩阵$>A=(I+\frac{\lambda }{\beta }{C}^{H}C).$>

第八步，计算中间变量其中Ψ^H,C^H分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x。通过以下三个式子更新ρ⁽ⁱ⁺¹⁾，α⁽ⁱ⁺¹⁾，ν⁽ⁱ⁺¹⁾

ρ⁽ⁱ⁺¹⁾＝F^Hx(5)

$>{\alpha }^{(i+1)}=max\left\{\right|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}|-\frac{1}{\beta },0\}\frac{{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}}{\left|{\mathrm{\Psi \rho }}^{(i+1)}+v^{\left(i\right)}\right|}---\left(6\right)$>

ν⁽ⁱ⁺¹⁾＝ν⁽ⁱ⁾+(Ψρ⁽ⁱ⁺¹⁾-α⁽ⁱ⁺¹⁾)(7)

若迭代次数i≥100或者则ρ＝ρ^(k+1)，并进行第九步，否则i值增加1，重复进行第八步。

第九步，若进行第十步，否则转至第五步。

第十步，最终输出矫正扭曲后的图像ρ，如图4所示。

本发明克服磁共振成像快速序列对强不均匀场的敏感和矫正不均匀场造成的图像扭曲。模型上，通过二维傅里叶正交变化将不均匀场完备地描述为四维矩阵的形式，再将四维矩阵降维成二维的大型矩阵。方法上主要是迭代进行两个模块，一是给定不均匀场，建立并求解了基于正交基变换和压缩感知的范数最优化模型；二是给定磁矩密度分布图像，使用解缠绕算法和布谷鸟最优化搜索算法，实现了拟合不均匀场。对强不均匀场和局部不均匀场下磁共振快速序列成像扭曲的矫正有显著的作用。