一种在强不均匀磁场下磁共振图像扭曲的矫正方法转让专利
申请号 : CN201610085326.8
文献号 : CN105738847B
文献日 : 2018-04-13
发明人 : 陈忠 , 庄孝星 , 蔡聪波
申请人 : 厦门大学
摘要 :
一种在强不均匀磁场下磁共振图像扭曲的矫正方法，涉及磁共振图像。模型上，通过二维傅里叶正交变化将不均匀场完备地描述为四维矩阵的形式，再将四维矩阵降维成二维的大型矩阵。方法上主要是迭代进行两个模块，一是给定不均匀场，建立并求解了基于正交基变换和压缩感知的l1范数最优化模型；二是给定磁矩密度分布图像，使用解缠绕算法和布谷鸟最优化搜索算法，实现了拟合不均匀场。对强不均匀场和局部不均匀场下磁共振快速序列成像扭曲的矫正有显著的作用。克服磁共振成像快速序列对强不均匀场的敏感和矫正不均匀场造成的图像扭曲。
权利要求 :
1.一种在强不均匀磁场下磁共振图像扭曲的矫正方法，其特征在于包括以下步骤：
1)计算不均匀场演化时间：t(kx,ky)表示自从施加RF脉冲后到k空间采样位置为(kx,ky)的不均匀场演化时间，表示为：其中tc为施加RF脉冲后到数据采样前的时间差值，EPI序列中tc数值很小，可以忽略不计；kx,ky的为k空间笛卡尔系正交的离散的坐标变量；k空间的正交两个坐标方向分别称为行坐标轴方向和列坐标轴方向；行坐标轴方向记为kx方向，列坐标轴方向记为ky方向；对应的,x,y为图像空间中正交的离散的两个坐标变量，图像空间正交的两个方向也称为相位编码方向和频率编码方向；相位编码方向记为y坐标轴方向，频率编码方向记为x坐标轴方向；
n(ky)表示列坐标轴方向采集点的坐标值为ky对应着的列数编号，Tro为采集一列k空间数据的时间，τpe为相位编码梯度单个周期内的施加时长；Gro为频率编码梯度，γ为旋磁比；
2)正交基变化系数矩阵C的构造函数：磁矩的演化相位为：
其中Binh(x,y)为不均匀场图，系数矩阵C(kx,ky,hx,hy)通过对 部分二维反傅里叶变化获得，即
其中hx,hy为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N；那么kx,ky,hx,hy,x,y离散点的个数也都分别等于N；令将四维的矩阵C(kx,ky,hx,hy)的转变为二维矩阵 其中 Lx为频率编码方向上的视野；n(ky)和n(hy)分别表示k空间列坐标轴方向上采集位置为ky和hy对应着的列数编号；根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(Binh)，其中 为 的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式；函数的输入为不均匀场Binh，输出为矩阵C；此函数称为的系数矩阵C的构造函数；
3)矩阵形式的采样模型：测量到的k空间数据用矩阵形式表示s＝CFρ；其中 为采集到的k空间数据， 为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量；按坐标次序是指将二维图像逐列首尾拼接为一维的向量；此时C＝C(Baut)表示是真实存在的而未知的不均匀场Baut对应的系数矩阵；
4)初始化拟合的不均匀场为Binh＝0，初始化重建图像ρ＝FHs；此时的ρ是受不均匀场影响而扭曲的图像；其中FH的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化；
5)更新拟合不均匀场：使用解缠绕算法计算图像ρ的解缠绕相位binh，并使用布谷鸟算法求解系数p；布谷鸟的最优化函数设置为 其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(Binh+pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造C矩阵；||·||2表示2范数；设置中间变量BT＝Binh；更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为：Binh＝BT+pbinh (4)
若 则进行步骤10)，否则转至步骤6)；
6)用已拟合的不均匀场Binh来构造系数矩阵C＝C(Binh)；建立基于正交基变换和压缩感知的最优化模型 其中α＝Ψρ表示为ρ经过小波变换后的系数；
7)设置模型参数β,λ，β取值范围在[5,1.5]，λ取值范围与β相关，为[20β,100β]；ρ(i),ν(i),α(i)表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…；其中ν为与α维度相同的中间变量；初始化i＝0，向量ν(0)＝α(0)＝0，ρ(0)＝ρ，构造矩阵
8)计算中间变量 其中ΨH,CH分别表示Ψ,C的共轭转置；使用
共轭梯度法求解线性方程Ax＝y得到方程的解x；通过以下三个式子更新ρ(i+1)，α(i+1)，ν(i+1)；
ρ(i+1)＝FHx (5)
ν(i+1)＝v(i)+(Ψρ(i+1)-α(i+1)) (7)若迭代次数i≥100或者 则ρ＝ρ(k+1)，并进行步骤9)，否则i值增加
1，重复进行步骤8)；
9)若 则进行步骤10)，否则转至步骤5)；
10)最终输出矫正扭曲后的图像ρ。
说明书 :
一种在强不均匀磁场下磁共振图像扭曲的矫正方法
技术领域
[0001] 本发明涉及磁共振图像，尤其是涉及一种在强不均匀磁场下磁共振图像扭曲的矫正方法。
背景技术
[0002] 磁共振成像具有良好临床安全性，广泛应用在医学诊断上。磁共振采集到的数据的空间称为k空间，磁共振图像的空间称为图像空间。自旋回波EPI序列是一种快速序列，只需要激发一个RF脉冲就可以采集整个k空间数据，能够有效地克服磁共振成像速度慢的缺点。传统的磁共振图像重建方法是对采集到的二维数据进行傅里叶变化。然而，EPI成像对不均匀场很敏感，傅里叶变化重建的图像在不均匀场区域通常会存在扭曲。
[0003] Zibetti(Zibetti MVW,De Pierro AR.A New Distortion Model for Strong Inhomogeneity Problems in Echo-Planar MRI[J].IEEE Transactions on Medical Imaging,2009,28(11):1736-1753.)提出了一种强不均匀场问题的模型，通过点扩散函数和离散分辨率的思想，可以在已知不均匀场的情况下重建图像。但是他通过调制离散分辨率水平来简化和加速算法的同时，也损失了许多细节信号，导致重建图像残留比较严重的伪影。
[0004] Donoho(Donoho DL.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4):1289-1306.)提出了压缩感知模型，若图像在变换域内是稀疏的，并且采样方式和稀疏变化的相干性足够小，就能够通过低于奈奎斯特采样率的数据来重建图像。Lustig(Lustig M,Donoho D,Pauly JM.Sparse MRI:The application of compressed sensing for rapid MR imaging[J].Magnetic Resonance in Medicine,
2007,58(6):1182-1195.)将压缩感知应用磁共振成像中。目前压缩感知在磁共振图像重建的主要应用包括超分辨重建，欠采样重建和去噪。
2007,58(6):1182-1195.)将压缩感知应用磁共振成像中。目前压缩感知在磁共振图像重建的主要应用包括超分辨重建，欠采样重建和去噪。
[0005] 在求取复数相位时候，求得的值总是限定在[-π,π]之间。事实上，超过这个区间的相位值会被取余从而被限定在区间内，由此产生了相位缠绕现象。Maier(Maier F,Fuentes D,Weinberg JS,et al.Robust phase unwrapping for MR temperature imaging using a magnitude-sorted list,multi-clustering algorithm[J].Magnetic Resonance in Medicine,2015,73(4):1662-1668.)提出了一种相位解缠绕方法，能够有效地从限定在[-π,π]的相位恢复到解缠绕的平滑的相位值。
[0006] 2009年，Yang(Yang X-S,Deb S,Engineering Optimisation by Cuckoo Search[J].Mathematical Modelling&Numerical Optimisation,2010,1(4):330-343)提出了一种新的启发式智能算法算法——布谷鸟搜索。其主要思想是通过模拟布谷鸟寄生育雏的行为以及Lévy飞行机制来有效地解决最优化问题。
[0007] 共轭梯度法是一种技术成熟的求解线性方程组的方法。Shewchuk JR.(Shewchuk JR,An introduction to the conjugate gradient method without the agonizing pain[M].Carnegie-Mellon University.Department of Computer Science.1994.)详细介绍了共轭梯度法的原理和实现方法。
发明内容
[0008] 本发明的目的在于提供一种在强不均匀磁场下磁共振图像扭曲的矫正方法。
[0009] 本发明包括以下步骤：
[0010] 1)计算不均匀场演化时间：t(kx,ky)表示自从施加RF脉冲后到k空间采样位置为(kx,ky)的不均匀场演化时间时间，表示为：
[0011]
[0012] 其中tc为施加RF脉冲后到数据采样前的时间差值，EPI序列中tc数值很小，可以忽略不计。kx,ky的为k空间笛卡尔系正交的离散的坐标变量。k空间的正交两个坐标方向分别称为行坐标轴方向(kx方向)和列坐标轴方向(ky方向)。对应的,x,y为图像空间中正交的离散的两个坐标变量，图像空间正交的两个方向也称为相位编码方向(y坐标轴方向)和频率编码方向(x坐标轴方向)。n(ky)表示列坐标轴方向采集点的坐标值为ky对应着的列数编号，Tro为采集一列k空间数据的时间，τpe为相位编码梯度单个周期内的施加时长。Gro为频率编码梯度，γ为旋磁比。
[0013] 2)正交基变化系数矩阵C的构造函数：磁矩的演化相位为：
[0014]
[0015] 其中Binh(x,y)为不均匀场图，系数矩阵C(kx,ky,hx,hy)通过对 部分二维反傅里叶变化获得，即
[0016]
[0017] 其中hx,hy为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N。那么kx,ky,hx,hy,x,y离散点的个数也都分别等于N。令将四维的矩阵C(kx,ky,hx,hy)的转变为二维矩阵 其中 Lx为
频率编码方向上的视野。n(ky)和n(hy)分别表示k空间列坐标轴方向上采集位置为ky和hy对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(Binh)，其中 为 的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
频率编码方向上的视野。n(ky)和n(hy)分别表示k空间列坐标轴方向上采集位置为ky和hy对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(Binh)，其中 为 的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
[0018] 3)矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CFρ。其中为采集到的k空间数据， 为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的向量。此时C＝C(Baut)表示是真实存在的而未知的不均匀场Baut对应的系数矩阵。
[0019] 4)初始化拟合的不均匀场为Binh＝0，初始化重建图像ρ＝FHs。此时的ρ是受不均匀场影响而扭曲的图像。其中FH的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。
[0020] 5)更新拟合不均匀场：使用解缠绕算法计算图像ρ的解缠绕相位binh(参见文献：Maier F,Fuentes D,Weinberg JS,et al.Robust phase unwrapping for MR
temperature imaging using a magnitude-sorted list,multi-clustering algorithm[J].Magnetic Resonance in Medicine,2015,73(4):1662-1668.)，并使用布谷鸟算法求解系数p(参见文献：Yang X-S,Deb S,Engineering Optimisation by Cuckoo Search[J].Mathematical Modelling&Numerical Optimisation,2010,1(4):330-343)。布谷鸟的最优化函数设置为 其中p的搜索范围为[-5,5]，|ρ|表示图像ρ
的模值，C(Binh+pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造C矩阵。||·||2表示2范数。设置中间变量BT＝Binh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为：
temperature imaging using a magnitude-sorted list,multi-clustering algorithm[J].Magnetic Resonance in Medicine,2015,73(4):1662-1668.)，并使用布谷鸟算法求解系数p(参见文献：Yang X-S,Deb S,Engineering Optimisation by Cuckoo Search[J].Mathematical Modelling&Numerical Optimisation,2010,1(4):330-343)。布谷鸟的最优化函数设置为 其中p的搜索范围为[-5,5]，|ρ|表示图像ρ
的模值，C(Binh+pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造C矩阵。||·||2表示2范数。设置中间变量BT＝Binh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为：
[0021] Binh＝BT+pbinh (4)
[0022] 若 则进行步骤10)，否则转至步骤6)；
[0023] 6)用已拟合的不均匀场Binh来构造系数矩阵C＝C(Binh)。建立基于正交基变换和压缩感知的最优化模型 其中α＝Ψρ表示为ρ经过小波变换后的系数。
[0024] 7)设置模型参数β,λ，β取值范围在[5,1.5]，λ取值范围与β相关，为[20β,100β]。(i) (i) (i)
ρ ,ν ,α 表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν(0)＝α(0)＝0，ρ(0)＝ρ，构造矩阵
ρ ,ν ,α 表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν(0)＝α(0)＝0，ρ(0)＝ρ，构造矩阵
[0025] 8)计算中间变量 其中ΨH,CH分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x(参见文献：Shewchuk JR,An
introduction to the conjugate gradient method without the agonizing pain[M].Carnegie-Mellon University.Department of Computer Science.1994.)。通过以下三个式子更新ρ(i+1)，α(i+1)，ν(i+1)；
introduction to the conjugate gradient method without the agonizing pain[M].Carnegie-Mellon University.Department of Computer Science.1994.)。通过以下三个式子更新ρ(i+1)，α(i+1)，ν(i+1)；
[0026] ρ(i+1)＝FHx (5)
[0027]
[0028] ν(i+1)＝ν(i)+(Ψρ(i+1)-α(i+1)) (7)
[0029] 若迭代次数i≥100或者 则ρ＝ρ(k+1)，并进行步骤9)，否则i值增加1，重复进行步骤8)；
[0030] 9)若 则进行步骤10)，否则转至步骤5)。
[0031] 10)最终输出矫正扭曲后的图像ρ。
[0032] 本发明的有益效果是：不均匀场场图采用拟合的方式获得的，无需采集参考场，避免额外的采集时间。由于傅里叶基正交完备，给定的不均匀场信息可以被完整的表达在系数矩阵C中，从而获得良好的矫正结果。同时，不均匀场的拟合方式对局部信息有良好的拟合效果，对产生的较强的局部的扭曲有明显矫正效果。
附图说明
[0033] 图1是EPI序列图。PE表示相位编码方向，FE表示频率编码方向。SS表示选层方向。在传统的EPI序列基础上增加不均匀场示意图。不均匀场Binh在整个采样过程可视为是不变的。
[0034] 图2是采集k空间数据的轨迹图。
[0035] 图3是未矫正的磁共振图像幅值图和相位图。(a)(b)是柠檬图像的幅值图和相位图；(c)(d)是鼠脑图像的幅值图和相位图。
[0036] 图4是矫正后磁共振图像的幅值图和相位图。(a)(b)是柠檬图像的幅值图和相位图；(c)(d)是鼠脑图像的幅值图和相位图。
具体实施方式
[0037] 本发明包括以下步骤：
[0038] 1)计算不均匀场演化时间：t(kx,ky)表示自从施加RF脉冲后到k空间采样位置为(kx,ky)的不均匀场演化时间时间，表示为
[0039]
[0040] 其中tc为施加RF脉冲后到数据采样前的时间差值，EPI序列中tc数值很小，可以忽略不计。kx,ky的为k空间笛卡尔系正交的离散的坐标变量。k空间的正交两个坐标方向分别称为行坐标轴方向(kx方向)和列坐标轴方向(ky方向)。对应的,x,y为图像空间中正交的离散的两个坐标变量，图像空间正交的两个方向也称为相位编码方向(y坐标轴方向)和频率编码方向(x坐标轴方向)。n(ky)表示列坐标轴方向采集点的坐标值为ky对应着的列数编号，Tro为采集一列k空间数据的时间，τpe为相位编码梯度单个周期内的施加时长。Gro为频率编码梯度，γ为旋磁比。
[0041] 采集k空间数据的轨迹图参见图2。
[0042] 2)正交基变化系数矩阵C的构造函数：磁矩的演化相位为
[0043]
[0044] 其中Binh(x,y)为不均匀场图，系数矩阵C(kx,ky,hx,hy)通过对 部分二维反傅里叶变化获得，即
[0045]
[0046] 其中hx,hy为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N。那么kx,ky,hx,hy,x,y离散点的个数也都分别等于N。令将四维的矩阵C(kx,ky,hx,hy)的转变为二维矩阵 其中 Lx为频
率编码方向上的视野。n(ky)和n(hy)分别表示k空间列坐标轴方向上采集位置为ky和hy对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(Binh)，其中 为 的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
率编码方向上的视野。n(ky)和n(hy)分别表示k空间列坐标轴方向上采集位置为ky和hy对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵的过程的写成函数的形式，写作C＝C(Binh)，其中 为 的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
[0047] 3)矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CFρ。其中为采集到的k空间数据， 为磁共振图像，F的作用是先将图像ρ二维傅里叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的向量。此时C＝C(Baut)表示是真实存在的而未知的不均匀场Baut对应的系数矩阵。
[0048] 4)初始化拟合的不均匀场为Binh＝0，初始化重建图像ρ＝FHs。此时的ρ是受不均匀场影响而扭曲的图像。其中FH的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。
[0049] 5)更新拟合不均匀场：使用解缠绕算法计算图像ρ的解缠绕相位binh(Maier F,Fuentes D,Weinberg JS,et al.Robust phase unwrapping for MR temperature imaging using a magnitude-sorted list,multi-clustering algorithm[J].Magnetic Resonance in Medicine,2015,73(4):1662-1668.)，并使用布谷鸟算法求解系数p(Yang X-S,Deb S,Engineering Optimisation by Cuckoo Search[J].Mathematical Modelling&Numerical Optimisation,2010,1(4):330-343)。布谷鸟的最优化函数设置为其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(Binh+
pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造C矩阵。||·||2表示2范数。设置中间变量BT＝Binh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为
pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造C矩阵。||·||2表示2范数。设置中间变量BT＝Binh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为
[0050] Binh＝BT+pbinh (4)
[0051] 若 则进行步骤10，否则转至步骤6
[0052] 6)用已拟合的不均匀场Binh来构造系数矩阵C＝C(Binh)。建立基于正交基变换和压缩感知的最优化模型 其中α＝Ψρ表示为ρ经过小波变换后的系数。
[0053] 7)设置模型参数β,λ，β取值范围在[5,1.5]，λ取值范围与β相关，为[20β,100β]。ρ(i),ν(i),α(i)表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν(0)＝α(0)＝0，ρ(0)＝ρ，构造矩阵
[0054] 8)计算中间变量 其中ΨH,CH分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x(Shewchuk JR,An introduction to the conjugate gradient method without the agonizing pain[M].Carnegie-Mellon University.Department of Computer Science.1994.)。通过以下三个式子更新ρ(i+1)，α(i+1)，ν(i+1)
[0055] ρ(i+1)＝FHx (5)
[0056]
[0057] ν(i+1)＝ν(i)+(Ψρ(i+1)-α(i+1)) (7)
[0058] 若迭代次数i≥100或者 则ρ＝ρ(k+1)，并进行步骤9，否则i值增加1，重复进行步骤8。
[0059] 9)若 进行步骤10，否则转至步骤5。
[0060] 10)最终输出矫正扭曲后的图像ρ。
[0061] 以下实施例使用7T(299.8MHz)/160-mm孔径MRI系统对柠檬和鼠脑进行成像。使用EPI序列实验来采集数据，序列图如图1所示，实验参数记录见表1。
[0062] 表1
[0063]实验参数 柠檬实验 鼠脑实验
频率编码方向视野Lx(mm) 80 45
频率编码梯度Gro(T·m-1) 0.0734 0.1305
相位编码方向视野Ly(mm) 80 45
相位编码梯度Gpe(T·m-1) 0.0245 0.0326
采样点数Nx,Ny 128 64
频率编码方向视野Lx(mm) 80 45
频率编码梯度Gro(T·m-1) 0.0734 0.1305
相位编码方向视野Ly(mm) 80 45
相位编码梯度Gpe(T·m-1) 0.0245 0.0326
采样点数Nx,Ny 128 64
[0064] 第一步，计算不均匀场演化时间：tc为施加RF脉冲后到数据采样前的时间差值，EPI序列中tc数值很小，可以忽略不计,所以tc设置为0。两个方向数值分辨率N＝128，旋磁比设置为水的旋磁比γ＝2π×42.576×106。计算其他参数的值 其中Tro为采集一列k空间数据的时间，τpe为相位编码梯度单个周期内的施加时长。kx,ky的为k空间笛卡尔系正交离散的坐标变量。kx,ky取值范围是
n(ky)表示列坐标轴方向采集坐标值
为ky对应着的列数编号，代入tc,Tro,τpe n(ky),kx,Gro,γ的数值得到从施加RF脉冲后到k空间采样位置为(kx,ky)的时间，即不均匀场的演化时间t(kx,ky)，表示为：
n(ky)表示列坐标轴方向采集坐标值
为ky对应着的列数编号，代入tc,Tro,τpe n(ky),kx,Gro,γ的数值得到从施加RF脉冲后到k空间采样位置为(kx,ky)的时间，即不均匀场的演化时间t(kx,ky)，表示为：
[0065]
[0066] 第二步，正交基变化系数矩阵C的构造函数：磁矩的演化相位为
[0067]
[0068] Binh(x,y)为不均匀场图，系数矩阵C(kx,ky,hx,hy)通过对 的部分二维反傅里叶变化获得，即
[0069]
[0070] 其中hx,hy为离散的频率坐标变量，磁共振图像两个方向数值分辨率设置为相同的数值，该数值即为N＝128。那么kx,ky,hx,hy,x,y离散点的个数也都分别等于128。令将C(kx,ky,hx,hy)的转变为二维矩阵 其中n(ky)和n(hy)分别表示相位编码方向上采集位置为ky和hy对应着的列数编号。根据公式(1)(2)(3)将构造系数矩阵C的过程的写成函数的形式，写作C＝C(Binh)，其中
的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀
场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
的矩阵形式，Binh为不均匀场Binh(x,y)的矩阵形式。函数的输入为不均匀
场Binh，输出为矩阵C。此函数称为的系数矩阵C的构造函数。
[0071] 第三步，矩阵形式的采样模型：测量到的k空间数据可以用矩阵形式表示s＝CFρ。其中 为采集到的k空间数据， 为磁共振图像，F的作用是先将图像ρ二维傅里
叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的列向量。此时C＝C(Baut)表示是真实存在的而未知的不均匀场Baut对应的系数矩阵。
叶变化后再按照坐标次序转化为一维向量。按坐标次序是指将二维图像逐列首尾拼接为一维的列向量。此时C＝C(Baut)表示是真实存在的而未知的不均匀场Baut对应的系数矩阵。
[0072] 第四步，初始化拟合的不均匀场为Binh＝0，重建图像ρ＝FHs。其中FH的作用是先按坐标次序的逆过程将一维向量s转变为二维矩阵，再对二维矩阵进行二维反傅里叶变化。此时的ρ是受不均匀场影响扭曲的图像，如图3所示。
[0073] 第五步，更新拟合不均匀场：使用解缠绕算法获得图像ρ的解缠绕相位binh，并使用布谷鸟算法求解系数p。布谷鸟的最优化函数设置为 其中p的搜索范围为[-5,5]，|ρ|表示图像ρ的模值，C(Binh+pbinh)表示用矩阵C的构造函数和不均匀场Binh+pbinh来构造系数矩阵，||·||2表示2范数。设置中间变量BT＝Binh。更新的拟合的不均匀场是已拟合不均匀场叠加上解缠绕相位binh和系数p的乘积，表示为
[0074] Binh＝BT+pbinh (4)
[0075] 若 则进行第十步，否则转至第六步
[0076] 第六步，用已拟合的不均匀场Binh来构造系数矩阵C＝C(Binh)。建立基于正交基变换和压缩感知的最优化模型 其中α＝Ψρ表示为ρ经过小波变换后的系数。
[0077] 第七步，设置模型参数β＝8,λ＝400。ρ(i),ν(i),α(i)表示ρ,ν,α在第i次迭代的数值，i＝0,1,2…。其中ν为与α维度相同的中间变量。初始化i＝0，向量ν(0)＝α(0)＝0，ρ(0)＝ρ，构造矩阵
[0078] 第八步，计算中间变量 其中ΨH,CH分别表示Ψ,C的共轭转置。使用共轭梯度法求解线性方程Ax＝y得到方程的解x。通过以下三个式子更新ρ(i+1)，α(i+1)，ν(i+1)
[0079] ρ(i+1)＝FHx (5)
[0080]
[0081] ν(i+1)＝ν(i)+(Ψρ(i+1)-α(i+1)) (7)
[0082] 若迭代次数i≥100或者 则ρ＝ρ(k+1)，并进行第九步，否则i值增加1，重复进行第八步。
[0083] 第九步，若 进行第十步，否则转至第五步。
[0084] 第十步，最终输出矫正扭曲后的图像ρ，如图4所示。
[0085] 本发明克服磁共振成像快速序列对强不均匀场的敏感和矫正不均匀场造成的图像扭曲。模型上，通过二维傅里叶正交变化将不均匀场完备地描述为四维矩阵的形式，再将四维矩阵降维成二维的大型矩阵。方法上主要是迭代进行两个模块，一是给定不均匀场，建立并求解了基于正交基变换和压缩感知的 范数最优化模型；二是给定磁矩密度分布图像，使用解缠绕算法和布谷鸟最优化搜索算法，实现了拟合不均匀场。对强不均匀场和局部不均匀场下磁共振快速序列成像扭曲的矫正有显著的作用。