CN104297709B - 一种基于正则化磁场分布图像重建的梯度匀场方法 - Google Patents

Abstract

本发明公开了一种基于正则化磁场分布图像重建的梯度匀场方法。与现有的相位差成像的梯度匀场技术相比，本发明提出的图像重建的优化方法通过三次相位差扫描成像获取初始场图，采用加入“平滑”和“先验”势函数的罚似然估计方法建立场图数据的罚似然估计函数，并通过优化变换进行函数的最小化，从而获得最终的场图。本方法一方面能有效考虑像素点的连续性，抑制噪声和图像毛刺；另一方面，先验势函数能预估数据的目标形状，避免优化变换过程走向错误的方向，针对较不均匀的磁场引起的图像畸变和缠绕有很好的校正效果。本方法保证了场图的完整性、提高了信噪比，使得即使在磁场均匀性极差的条件(如“冷场”等)也能获得一定的匀场效果。

Description

一种基于正则化磁场分布图像重建的梯度匀场方法

技术领域

本发明涉及核磁共振波谱仪的梯度匀场技术，更具体涉及一种核磁共振波谱仪上基于正则化磁场分布图像重建的梯度匀场的方法。

背景技术

对于现代核磁共振波谱仪来说，静磁场的均匀性直接影响采集信号的分辨率。为了提高静磁场的均匀性，大多数的谱仪都会采用一系列通有电流的匀场线圈从而产生特定的空间磁场来消除静磁场的不均匀分量(主动匀场)。由于匀场线圈正交设计和空间磁场分布所具有的复杂性，传统的手动匀场方式往往需要实验者过多的人为干预，繁琐的操作和大量宝贵机时的耗费都不利于实验的顺利进行。最早的自动匀场方法出现于1968年(R.R.Ernst.Measurement and control of magnetic field homogeneity,The review of scientific instruments.39,1968,998-1012)，其通过判断采样信号的特性来搜索迭代匀场。随着梯度在核磁共振成像中的应用和发展，1994年Zijl首次提出了应用于核磁共振波谱仪的梯度匀场方法(P.C.M.Van Zijl,S.Sukumar,M.O’NeilJohnson,et al.Optimized shimming for high-resolution NMR using three-dimensionalimage-based field mapping,Journal of magnetic resonance A.111,1994,203-207)，由于所采用的梯度回波脉冲序列能够测量并拟合静磁场不均匀性的空间分布，梯度匀场被认为是核磁共振技术中最直接、有效的自动匀场方法。

梯度匀场必须首要获得静磁场不均匀性的像(静磁场场图)和匀场线圈的单位电流作用的磁场影响(匀场线圈场图)，能否准确表征所需的场图(包括静磁场场图和特定匀场线圈场图)是梯度匀场是否有效的关键。通常的磁场场图测量方法是，执行两次梯度回波脉冲序列采样，通过计算两次扫描的相位差和成像回波时间差来拟合场图的频率数据。这种直接通过相位差拟合场图数据的方式存在两个问题：1)没有考虑噪声的影响。通常来说，在场图的拟合计算时，噪声会被近似忽略，但是实际测量过程中，其是真实存在甚至是很大的(例如样品氘含量很弱的情况)，这样往往会使像的质量变差，特别是针对自旋密度较低的成像区域，容易使得较弱的有用信号埋没在噪声中。2)不适用于极差的磁场条件。传统的相位差重建算法为了不考虑采样期间静磁场空间不均匀性自身对信号产生的干扰，因此会在进行梯度匀场前将磁场调整到一个较为均匀的状态。但是对于极差的条件(如“冷场”环境等)，磁场不均匀性会使得最终拟合的场图分布不连续、畸变，从而直接影响匀场效果。

随着磁共振成像技术在医疗中的广泛应用，对于磁共振信号和图像质量要求越来越高，各种针对图像畸变伪影校正的研究快速发展。2008年，Funai利用成像的像素点自身具有连续性这一性质校正了噪声很大的图像数据(Amanda K.Funai,Jeffrey A.Fessler,Desmond T.B.Yeo,et al.Regularized Field Map Estimation in MRI,IEEETRANSACTIONS ON MEDICAL IMAGING.27,2008,1484-1494)，具体做法是多次扫描的相位差和成像回波时间差来拟合初始图像，利用罚似然估计(Penalized-Likelihood Estimation)的方法(加入“平滑”的势函数)抑制噪声影响，通过优化变换方式(Optimization Transformation,OT)来迭代优化并提高图像的质量。这种方式能有效校正场图像素点邻近区域的不连续性使其平滑，但是当该区域的磁场不均匀过大时，畸变、缠绕的信号也无法恢复成理想的状态，有时甚至会使优化变换过程沿着“畸变”点的方向使之更严重。

本申请提出一种核磁共振波谱仪上基于正则化磁场分布图像重建的梯度匀场方法，通过三次扫描相位差成像获取初始场图，利用罚似然估计的方法，以平滑和先验函数作为双重评价标准来迭代优化场图数据：一方面“平滑”函数能考虑像素点的连续性，抑制噪声和图像毛刺；另一方面，“先验”函数能预估数据的目标形状，避免了优化变换过程走向错误的方向，针对较不均匀的磁场引起的图像畸变和缠绕有很好的校正效果。采用本正则化算法重建磁场分布图像并进行梯度匀场，和常规梯度匀场相比，保证了场图的完整性、提高了信噪比，并且获得了很好的匀场效果。

发明内容

本发明的目的是提供一种磁共振波谱仪上基于正则化磁场分布图像重建的梯度匀场方法，为了实现上述目的，本发明采用以下技术方案：

一种基于正则化磁场分布图像重建的梯度匀场方法，包括以下步骤：

步骤1、根据需要进行梯度匀场的N个匀场线圈选择相应的梯度回波脉冲序列；

步骤2、对初始化匀场线圈的磁场状态执行梯度回波脉冲序列采样，根据回波采样数据初始化匀场线圈场图，并通过正则化磁场分布图像重建获得最终的匀场线圈场图；

步骤3、对待匀场的磁场状态执行梯度回波脉冲序列采样，根据回波采样数据初始化待匀场的静磁场场图，并通过正则化磁场分布图像重建获得最终的静磁场场图；

步骤4、将步骤2获得的最终的匀场线圈场图表示为矩阵A，将步骤3获得的最终的静磁场场图数据标识为向量b，求解线性方程组Ax＝b得到待匀场线圈的电流改变量；

步骤5、将计算的匀场线圈电流改变量分别加至各待匀场的匀场线圈电流值上，若磁场均匀度无改善或是均匀度满足预定设置的要求，则完成所有匀场；否则将当前N个待匀场的匀场线圈的电流值作为新的待匀场的N个匀场线圈的电流值返回步骤3。

如上所述的步骤2中包括以下步骤：

步骤2.1、初始化N个匀场线圈的电流值Value(j),j＝1,2,3...N，将初始化通有匀场线圈电流的磁场状态设为基础磁场状态，即磁场状态0，进行梯度回波采样，获得成像回波数据；

步骤2.2、设置所选的各个匀场线圈在采样过程中的电流变化量△Change(j)；

步骤2.3、将每个匀场线圈对应电流变化量△Change(j)依次叠加在基础磁场状态下对应的匀场线圈上，得到磁场状态1～N，并进行采样，得到叠加电流变化量后的成像回波数据；

步骤2.4、对磁场状态0～磁场状态N的成像回波数据进行傅里叶变换，得到各个磁场状态下的信号强度和相位：

信号强度：(s₁(r,i),TE₁),(s₂(r,i),TE₂),(s₃(r,i),TE₃)

相位：(φ₁(r,i),TE₁),(φ₂(r,i),TE₂),(φ₃(r,i),TE₃)

其中，信号强度s₁～s₃、相位数据φ₁～φ₃均和回波采样的时间TE₁～TE₃相关联；当i＝0时为基础磁场采样得到的信号强度和相位数据，i＝1,2,3…N时表示依次改变N个匀场线圈电流值后采样得到的信号强度和相位数据；r表示计算的有效像素点，r的个数为N_P；

步骤2.5、将磁场状态1～N状态下采样得到的信号强度和相位数据与磁场状态0采样得到的信号强度和相位数据分别做差，得到表征每个匀场线圈单独作用效果的信号强度y_k(r,j)和相位ph_k(r,j)；k＝1,2,3表示单个匀场线圈所对应的三次扫描的成像回波时间TE₁、TE₂和TE₃；

步骤2.6、将步骤2.5中的获得的表征各个匀场线圈单独作用的相位中的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据ph₂和ph₁作差，得到成像的相位差并进行相位解缠△φ₂₁(r,j)＝unwrap(ph₂(r,j)-ph₁(r,j))；

步骤2.7、根据以下公式初始化各个匀场线圈的场图(频率数据)：ω_Shim(r,j)＝△φ₂₁(r,j)/(TE₂-TE₁)；

步骤2.8、建立匀场线圈场图的罚似然估计函数Ψ_PL(ω_Shim(r,j))，然后通过优化变换的方式最小化罚似然估计函数Ψ_PL(ω_Shim(r,j))并获得最终的匀场线圈场图数据

如上所述的步骤2.8中建立罚似然估计函数Ψ_PL(ω_Shim(r,j))包括以下步骤：

步骤2.8.1、拟合各匀场线圈场图的最大似然估计函数Ψ_ML(ω_Shim(r,j))；

步骤2.8.2、利用二次函数凹凸性质建立场图数据ω_Shim的势函数，其中势函数包括平滑势函数R_smooth(r,j)和先验势函数R_prior(r,j)；

步骤2.8.3、将匀场线圈场图ω_Shim(r,j)的最大似然估计函数与势函数相加，建立ω_Shim(r,j)的罚似然估计函数Ψ_PL(ω_Shim(r,j))：

Ψ_PL(ω_Shim(r,j))＝Ψ_ML(ω_Shim(r,j))+α·R_smooth(r,j)+β·R_prior(r,j),j＝1,2,3…N

其中，α、β为势函数的系数。

如上所述的步骤3包括以下步骤：

步骤3.1、记录待匀场的N个匀场线圈的电流值CurrValue(j)，将此时通有匀场线圈电流的磁场状态设为待匀场的磁场状态，并获得采样数据；

步骤3.2、对步骤3.1的采样数据进行傅里叶变换，获得表征待匀场的静磁场像的信号强度和相位：

信号强度～((yb0₁(r),TE₁),(yb0₂(r),TE₂),(yb0₃(r),TE₃))

相位～((phb0₁(r),TE₁),(phb0₂(r),TE₂),(phb0₃(r),TE₃))

其中，信号强度yb0₁～yb0₃、相位数据phb0₁～phb0₃均分别和回波采样的时间TE₁～TE₃相关联；r表示计算的有效像素点，r的个数为N_P；

步骤3.3、将步骤3.2获得的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据phb0₂和phb0₁作差，得到成像的相位差并进行相位解缠，通过相位差和时间差(TE₂-TE₁)相除获得表征待匀场的静磁场场图的频率数据

${\mathrm{\ω}}_{B_0}\left(r\right)=unwrap\left(phb{0}_2\right(r)-phb{0}_1(r\left)\right)/({\mathrm{TE}}_2-{\mathrm{TE}}_1);$

步骤3.4、建立待匀场的静磁场场图的罚似然估计函数通过优化变换的方式最小化罚似然估计函数获得最终的静磁场场图数据

如上所述的步骤3.4中建立罚似然估计函数包括以下步骤：

步骤3.4.1、拟合待匀场的静磁场场图的最大似然估计函数

步骤3.4.2、利用二次函数凹凸性质建立待匀场的静磁场场图的势函数，其中势函数为平滑势函数

步骤3.4.3、将静磁场场图的最大似然估计函数与势函数相加，建立的罚似然估计函数

${\mathrm{\Ψ}}_{PL}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)={\mathrm{\Ψ}}_{ML}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)+{\mathrm{\α}}_{B_0}\·R_{smooth\_B_0}\left(r\right)$

其中，为势函数的系数。

与现有技术相比，本发明的优点和有益效果是：

1、能有效抑制噪声及磁场不均匀性对场图数据的影响。现有的梯度匀场技术无法从真实采样的角度考虑噪声或磁场不均匀性对场图数据的污染，导致最终匀场结果不理想。本申请充分利用采样期间不同成像回波时间对场图数据的影响(成像回波时间短，磁场不均匀性的影响小，信噪比低，但场图数据更完整；成像回波时间长，磁场不均匀性的影响大，场图数据可能有畸变，但信噪比高。)，通过将两次不同回波成像时间的场图数据合二为一达到优化的效果：首先初始化成像时间较短的数据，保证场图形状的完整性；再利用成像时间较长的数据来修正成像时间较短的数据，从而使得噪声被抑制。

2、能用于磁场均匀性极差的条件进行匀场。一般在使用梯度匀场前，会依靠手动匀场或是其他自动匀场方式(例如自动搜索匀场等)将磁场环境调整到一个相对均匀的状态，这是因为当静磁场环境变得很差时，磁场不均匀性自身产生的干扰和污染远远大于真实信号，使得场图的畸变和缠绕很严重，根本无法匀场。本申请采用加入“平滑”和“先验”势函数的罚似然估计方法，通过优化变换保证重建的图像无畸变且信噪比较高，使得即使在磁场均匀性极差的条件(如“冷场”等)也能获得一定的匀场效果。

附图说明

图1是梯度匀场方法的一般流程图。

图2(1)是采用正则化图像重建方式拟合场图的流程图。

图2(2)是直接通过相位差拟合场图的流程图。

图3(1)a～g、图3(2)a～g、图3(3)a～g、图3(4)a～g、图3(5)a～g是表示采用正则化图像重建方式拟合的Z1-Z7匀场线圈场图的过程。

如图所示，横坐标EffectNum表示场图有效区域的像素点N_P，纵坐标ShimFunc表示匀场线圈单位电流作用对应的空间磁场分布形状，单位为特斯拉(T)；实验样品为采用加弛豫试剂的水溶液(Auto Test Sample，0.3mg/ml GdCl₃,0.1％¹³C-METHANOL,1％H₂O，98.9％D₂O)，实验时磁体的静磁场为“冷场”状态(所有匀场线圈电流的配置值为零)。

图3(1)a～g为两次成像回波时间TE分别为0s和0.005s，通过相位差拟合的匀场线圈场图；

图3(2)a～g为两次成像回波时间TE分别为0s和0.1s，通过相位差拟合的匀场线圈场图；

通过拟合三次回波相位数据(TE＝[0,0.005,0.1]s),采用本正则化图像重建方法迭代优化匀场线圈场图；

图3(3)a～g为第一次优化迭代的结果，

图3(4)a～g为第二次优化迭代的结果，

图3(5)a～g为第三次优化迭代的结果。

图4(1)～图4(3)表示采用本正则化磁场分布图像重建的梯度匀场方法针对加弛豫试剂的水溶液样品(Auto Test Sample，0.3mg/ml GdCl₃,0.1％¹³C-METHANOL,1％H₂O，98.9％D₂O)匀场前后的谱图比较。

图4(1)表示“冷场”环境，极差的静磁场不均匀性测得的¹H谱图信号；

图4(2)表示采用相位差成像方式(两次成像回波时间分别为0s和0.005s)进行Z1-Z7匀场线圈的一维梯度匀场，图示为匀场后¹H谱图信号无改善；

图4(3)表示采用基于正则化磁场分布图像重建的一维梯度匀场(三次成像回波时间分别为0s、0.005s和0.1s)针对Z1-Z7匀场线圈匀场后，¹H谱图信号有明显改善。

具体实施方式

本实施例使用中科波谱WNMR-Ⅰ500MHz核磁共振谱仪，其中配备的匀场线圈为Oxford28-shims。下面结合附图，对本发明的技术方案作进一步的详细描述。

如图1流程图所示，一种基于正则化磁场分布图像重建的梯度匀场方法，包括以下步骤：

选择需要进行梯度匀场的N个匀场线圈1：根据所选的匀场线圈组合自动设定一维或三维梯度回波的脉冲序列；读入并调整默认的脉冲序列参数，具体如下：

步骤1.1、根据步骤1所选的室温匀场线圈组合设定梯度回波的脉冲序列，分为以下两种情况：1)当匀场组合只含有Z轴方向(轴向或On-Axis)的匀场线圈(如Z1～Z7)时，采用一维梯度回波(PFGSTE)脉冲序列并需要三次成像采样；2)当匀场组合包含非Z轴(径向或Off-Axis)的匀场线圈(如X，ZX，X2-Y2，XY，X3等)时，采用倾斜三维梯度回波(Ramped 3D GRE)脉冲序列，此时X和Y方向还需要进行相位编码，相位编码数通常设定为N_X＝N_Y＝8或16或32。

步骤1.2、读入脉冲序列默认参数：一维梯度回波(PFGSTE)脉冲序列和倾斜三维梯度回波(Ramped 3D GRE)脉冲序列的参数设置分别如下表1和表2所示；

表1一维梯度回波(PFGSTE)脉冲序列参数设置

表2倾斜三维梯度回波(Ramped 3D GRE)脉冲序列参数设置

步骤1.3、调整参数：根据样品中氘代含量的差别设置实验累加次数ns＝1～4；调节射频脉冲产生的翻转角α、调节接收机增益(确保信号强度大且不溢出)；调节脉冲前等待时间d1，确保每次梯度回波信号形状无变形且强度无减弱。

正则化图像重建拟合匀场线圈的场图2：具体包括以下步骤：

步骤2.1、初始化N个匀场线圈电流值Value(j),j＝1,2,3...N，将初始化通有匀场线圈电流的磁场状态设为基础磁场状态(磁场状态0：Value(1),Value(2),...,Value(N))，并进行梯度回波采样，获得成像回波数据：当采用一维梯度回波脉冲序列时，需要进行三次回波采集；当执行倾斜三维梯度回波脉冲序列时，还需要进行X与Y方向的相位编码N_X×N_Y，则总回波数据采集为3×N_X×N_Y。

步骤2.2、设置所选的各个匀场线圈在采样过程中的电流变化量△Change(j),j＝1,2,3...N：其中，N为匀场线圈的个数，匀场线圈的电流大小可参考3～100mA的范围进行数值调节；随着匀场线圈阶数的增加，在不超过允许的调节范围内，电流变化量应增大。

步骤2.3、将步骤2.2预设的电流变化量△Change(j)依次添加到对应的匀场线圈的电流值Value(j)上，即，将每个匀场线圈对应电流变化量△Change(j)依次叠加在基础磁场状态下对应的匀场线圈上，由于匀场线圈电流的变化导致静磁场状态的改变，磁场状态1～N分别表示由于单一匀场线圈的电流变化作用引起的静磁场状态改变，具体如下：

仍然采用与步骤2.1相同的采样方式进行采样，获得成像回波数据。

步骤2.4、对步骤2.1中磁场状态0及步骤2.3中磁场状态1～N的所有采样数据进行傅里叶变换，获得一系列表征不同磁场状态的信号强度和相位：

信号强度：(s₁(r,i),TE₁),(s₂(r,i),TE₂),(s₃(r,i),TE₃),i＝0,1,2,...,N

相位：(φ₁(r,i),TE₁),(φ₂(r,i),TE₂),(φ₃(r,i),TE₃),i＝0,1,2,...,N

其中，信号强度s₁～s₃、相位数据φ₁～φ₃均分别和回波采样的时间TE₁～TE₃相关联；i表征不同的磁场状态，当i＝0时为基础磁场(磁场状态0)采样得到的信号强度和相位数据，当i＝1,2,3…N时表示依次改变N个匀场线圈电流值后(磁场状态1～N)采样得到的信号强度和相位数据；r表示计算的有效像素点，r的个数N_P取决于采用哪种梯度匀场方式(一维或是三维)：当采用一维梯度匀场时，N_P就等于表征磁场状态的信号强度和相位在Z方向上有效点的个数N_z(N_z为Z方向不低于最大信号强度30％的区域范围的点个数)；当采用三维梯度匀场时，不仅需要考虑Z方向上有效点的个数，还需要考虑X和Y方向的相位编码数，因此N_P＝N_Z×N_X×N_Y。

步骤2.5、将磁场状态1～N下采样得到的信号强度和相位数据与基础磁场(磁场状态0)采样得到的信号强度和相位数据分别做差，得到表征每个匀场线圈单独作用效果的信号强度和相位：

表征N个匀场线圈单独作用效果的信号强度：

y_k(r,j)＝(s_k(r,j)-s_k(r,0)),j＝1,2,3...N,k＝1,2,3；

表征N个匀场线圈单独作用效果的相位：

ph_k(r,j)＝(φ_k(r,j)-φ_k(r,0)),j＝1,2,3...N,k＝1,2,3

其中j＝1,2,3…N表示对应的N个匀场线圈；k＝1,2,3表示单个匀场线圈所对应的三次扫描的成像回波时间TE₁、TE₂和TE₃；r表示有效像素点。

步骤2.6、将步骤2.5中的获得的表征各个匀场线圈单独作用的相位中的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据ph₂和ph₁作差，得到成像的相位差并进行相位解缠：

△φ₂₁(r,j)＝unwrap(ph₂(r,j)-ph₁(r,j)),j＝1,2,3…N；

其中，unwrap表示对相位数据进行解缠操作；j＝1,2,3…N表示对应的N个匀场线圈；r表示有效像素点。

步骤2.7、初始化各个匀场线圈的场图(频率数据)：

ω_Shim(r,j)＝△φ₂₁(r,j)/(TE₂-TE₁),j＝1,2,3...N

其中j＝1,2,3…N表示对应的N个匀场线圈；r表示有效像素点。

步骤2.8、对各个初始化的匀场线圈场图数据ω_Shim(r,j)进行正则化图像重建：首先利用罚似然估计的方法建立匀场线圈场图的罚似然估计(Penalized-LikelihoodEstimation)函数Ψ_PL(ω_Shim(r,j))，即将最大似然估计(Maximum-Likelihood Estimation)函数Ψ_ML(ω_Shim(r,j))加入“平滑”R_smooth(r,j)和“先验”R_prior(r,j)这两个势函数；然后通过优化变换(Optimization Transformation,OT)的方式来最小化匀场线圈的罚似然函数Ψ_PL(ω_Shim(r,j))，从而抑制噪声和数据畸变的影响、提高图像的质量并获得最终的匀场线圈场图数据

步骤2.8.1、拟合各匀场线圈场图的最大似然估计函数Ψ_ML(ω_Shim(r,j))：

$\begin{array}{c}{\mathrm{\Ψ}}_{ML}\left({\mathrm{\ω}}_{Shim}\right(r,j\left)\right)=\underset{r=1}{\overset{N_P}{\Σ}}\underset{m=1}{\overset{3}{\Σ}}\underset{n=1}{\overset{3}{\Σ}}|y_m(r,j)\·y_n(r,j)|\·w_{mn}(r,j)\\ \·\[1-\cos \left({\mathrm{ph}}_n\right(r,j)-{\mathrm{ph}}_m(r,j)-{\mathrm{\ω}}_{Shim}(r,j)\·({\mathrm{TE}}_n-{\mathrm{TE}}_m\left)\right)\],j=1,2,\mathrm{3...}N\end{array}$

其中，ω_Shim(r,j)表示初始化的匀场线圈场图数据；r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈；TE表示扫描的成像回波时间；y表示与扫描的成像回波时间相对应的信号强度；ph表示与扫描的成像回波时间相对应的相位数据；相关权重因子w_mn(r,j)为：

步骤2.8.2、利用二次函数(y＝x²/2)的“凹凸”性质建立场图数据ω_Shim(r,j)的势函数，主要包括平滑势函数R_smooth(r,j)和先验势函数R_prior(r,j)这两部分，

1)场图数据的平滑势函数R_smooth：

若为一维梯度匀场数据，只需要考虑Z方向连续性：

$R_{smooth}(r,j)=\underset{r=2}{\overset{N_P-1}{\Σ}}{(2\·{\mathrm{\ω}}_{Shim}(r,j)-{\mathrm{\ω}}_{Shim}(r-1,j)-{\mathrm{\ω}}_{Shim}(r+1,j\left)\right)}^2/2,j=1,2,\mathrm{3...}N$

其中，r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈。

若为三维梯度匀场数据，则在XYZ三个方向的连续性都需要考虑，需要将一维有效像素点r还原成三个维度方向上的有效像点

r→(a,b,c),a＝1,2,3...N_X,b＝1,2,3...N_Y,c＝1,2,3...N_Z，则平滑势函数R_smooth表示如下：

$\begin{array}{c}{R}_{smooth}(r,j)=R_{smooth}\left(\right(a,b,c),j)=\\ \underset{c=1}{\overset{N_Z-1}{\Σ}}\underset{b=1}{\overset{N_Y-1}{\Σ}}\underset{a=2}{\overset{N_X-1}{\Σ}}{(2\·{\mathrm{\ω}}_{Shim}(\left(a,b,c\right),j)-{\mathrm{\ω}}_{Shim}(\left(a-1,b,c\right),j)-{\mathrm{\ω}}_{Shim}(\left(a+1,b,c\right),j\left)\right)}^2/2\\ +\underset{c=1}{\overset{N_Z-1}{\Σ}}\underset{b=2}{\overset{N_Y-1}{\Σ}}\underset{a=1}{\overset{N_X-1}{\Σ}}{(2\·{\mathrm{\ω}}_{Shim}(\left(a,b,c\right),j)-{\mathrm{\ω}}_{Shim}(\left(a,b-1,c\right),j)-{\mathrm{\ω}}_{Shim}(\left(a,b+1,c\right),j\left)\right)}^2/2\\ +\underset{c=2}{\overset{N_Z-1}{\Σ}}\underset{b=1}{\overset{N_Y-1}{\Σ}}\underset{a=1}{\overset{N_X-1}{\Σ}}{(2\·{\mathrm{\ω}}_{Shim}(\left(a,b,c\right),j)-{\mathrm{\ω}}_{Shim}(\left(a,b,c-1\right),j)-{\mathrm{\ω}}_{Shim}(\left(a,b,c+1\right),j\left)\right)}^2/2\\ ,j=1,2,\mathrm{3...}N\end{array}$

其中，a、b和c分别表示XYZ三个维度方向的有效像素点；j＝1,2,3…N表示对应的N个匀场线圈。

2)场图数据的先验势函数R_prior：

由于各匀场线圈有特定的磁场分布形状，因此可以通过先验势函数R_prior来判断场图是否有效：通过最小二乘法拟合匀场线圈场图的目标函数f(r,j)：

argmin[ω_Shim(r,j)-f(r,j)]²,j＝1,2,3…N

则此时先验势函数为其中，r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈。

步骤2.8.3、将匀场线圈场图ω_Shim(r,j)的最大似然估计函数与势函数相加，建立ω_Shim(r,j)的罚似然估计函数Ψ_PL(ω_Shim(r,j))：

Ψ_PL(ω_Shim(r,j))＝Ψ_ML(ω_Shim(r,j))+α·R_smooth(r,j)+β·R_prior(r,j),j＝1,2,3…N

其中，α、β为势函数的系数，本实施例中α＝β＝5×10⁸；r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈。

步骤2.8.4、针对初始化的匀场线圈场图ω_Shim(r,j)采用优化变换(OT)中的SPS(separable paraboloidal surrogates)算法(Fessler J A,Erdogan H.A paraboloidalsurrogates algorithm for convergent penalized-likelihood emission image reconstruction.IEEE Nuclear Science Symposium and Medical Imaging Conference,2,1998,1132-1135.)进行优化获得最终的匀场线圈场图通过构造最大似然估计函数Ψ_ML(ω_Shim(r,j))的二次型代理函数并重组罚似然估计函数Ψ_PL(ω_Shim(r,j))，最小化罚似然估计函数Ψ_PL(ω_Shim(r,j))，从而得到优化后的匀场线圈场图

步骤2.8.4.1通过构造易于最大化的二次型代理函数(SPS,separableparabolidal surrogate)取代最大似然估计函数Ψ_ML(ω_Shim(r,j))，并对代理函数进行二阶Taylor展开：

其中，r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈；TE表示扫描的成像回波时间；y表示与扫描的成像回波时间相对应的信号强度；ph表示与扫描的成像回波时间相对应的相位数据；相关权重因子w_mn(r,j)为：

$\left|y_m(r,j)\right|\·\left|y_n(r,j)\right|/\underset{k=1}{\overset{3}{\Σ}}|y_k(r,j){|}^2.$

对代理函数进行二阶Taylor展开：

其中，s＝ω_Shim(r,j)，{s}_2π表示相位数据|s-2nπ|≤π,n＝1,2,3...，

步骤2.8.4.2、将构造的代理函数带入罚似然估计函数Ψ_PL(ω_Shim(r,j))进行计算，通过迭代最小化的优化方式，得到最终匀场线圈的场图

${\mathrm{\ω}}_{Shim}^{opt}(r,j)~{\mathrm{\ω}}_{Shim}^{l+1}(r,j)={\mathrm{\ω}}_{Shim}^l(r,j)-diag\left\{\frac{1}{d(r,j)+\mathrm{\ϵ}\·\mathrm{\λ}}\right\}\·\▿{\mathrm{\Ψ}}_{PL}\left({\mathrm{\ω}}_{Shim}^l\right(r,j\left)\right),j=1,2,\mathrm{3...}N$

其中，ε和λ为正则化系数，ε＝500，λ＝2×10⁸；r表示有效像素点；j＝1,2,3…N表示对应的N个匀场线圈。

如图3(1)a～g、3(2)a～g、3(3)a～g、3(4)a～g、3(5)a～g所示，是采用本正则化磁场分布图像重建的梯度匀场方法针对磁体“冷场”状态进行匀场线圈Z1-Z7场图的重建优化。

现有的梯度匀场通过计算两次成像回波的时间差和相位差拟合场图，由于基础磁场较差，场图数据极易受到噪声和磁场自身的不均匀性干扰，会出现以下两种情况：(1)当成像时间较短(图3(1)a～g)时，TE₂-TE₁＝0.005s，场图数据受磁场不均匀性影响小，形状基本能保证完整性，但由于成像时间过短，噪声大，特别是信号强度本身很小的高阶线圈(如Z5，Z6，Z7等)，极易使场图像素点被噪声所淹没；(2)当成像时间较长(图3(2)a～g)时，TE₂-TE₁＝0.1s，噪声对于信号的影响相应小，但是由于静磁场自身的不均匀性对采集数据影响过大，造成数据畸变严重，特别是边缘处的像素点有明显的变形和缠绕。

采用本优化重建算法能有效考虑成像时间长短对信号的影响，通过将正则化的罚似然估计函数和优化变换法结合起来迭代修正场图数据：(1)第一次迭代(图3(3)a～g)结束，由于低阶线圈(Z1，Z2，Z3)信号强度大，噪声、磁场不均匀性与信号本身比较影响较小，容易校正，但Z4-Z7线圈的场图仍有不同程度的较大的数据干扰和变形；(2)经过第二次迭代(图3(4)a～g)，Z4-Z7线圈的场图得到了一定程度的“平滑”，抑制了噪声和磁场不均匀性的污染；(3)第三次迭代(图3(5)a～g)结束时，Z4-Z7的数据得到进一步优化，所有匀场线圈的场图重建完成。

正则化图像重建拟合待匀场的静磁场场图3：具体包括以下步骤：

步骤3.1、记录待匀场的N个匀场线圈的电流值CurrValue(j),j＝1,2,3...N，将此时通有匀场线圈电流的磁场状态设为待匀场的磁场状态，使用与制作匀场线圈场图相同的脉冲序列进行采样并获得采样数据，方式同步骤2.2；

步骤3.2、对步骤3.1的采样数据进行傅里叶变换，获得一系列表征待匀场的静磁场像的信号强度和相位：

信号强度～((yb0₁(r),TE₁),(yb0₂(r),TE₂),(yb0₃(r),TE₃))

相位～((phb0₁(r),TE₁),(phb0₂(r),TE₂),(phb0₃(r),TE₃))

其中，信号强度yb0₁～yb0₃、相位数据phb0₁～phb0₃均分别和回波采样的时间TE₁～TE₃相关联；；r表示计算的有效像素点，r的个数N_P取决于采用哪种梯度匀场方式(一维或是三维)：当采用一维梯度匀场时，N_P就等于表征磁场状态的信号强度和相位在Z方向上有效点的个数N_z(为不低于Z方向点的最大信号强度30％的区域)；当采用三维梯度匀场时，不仅需要考虑Z方向上有效点的个数，还需要考虑X和Y方向的相位编码数，因此N_P＝N_Z×N_X×N_Y。

步骤3.3、初始化待匀场静磁场的场图将步骤3.2获得的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据phb0₂和phb0₁作差，得到成像的相位差并进行相位解缠，最终获得表征静磁场的频率数据：

${\mathrm{\ω}}_{B_0}\left(r\right)=unwrap\left(phb{0}_2\right(r)-phb{0}_1(r\left)\right)/({\mathrm{TE}}_2-{\mathrm{TE}}_1)$

其中r表示有效像素点，unwrap表示对相位数据进行解缠操作。

步骤3.4、同步骤2.8，对初始化的待匀场静磁场场图进行正则化图像重建：利用罚似然估计的方法建立罚似然估计函数即将最大似然估计函数加入“平滑”R_smooth(r)的势函数；通过优化变换(OT)的方式最小化静磁场场图的罚似然函数从而优化提高图像的质量并获得最终的静磁场场图数据

步骤3.4包括以下步骤：

步骤3.4.1、拟合待匀场的静磁场场图的最大似然估计函数

步骤3.4.2、利用二次函数凹凸性质建立待匀场的静磁场场图的势函数，其中势函数为平滑势函数

步骤3.4.3、将静磁场场图的最大似然估计函数与势函数相加，建立的罚似然估计函数

${\mathrm{\Ψ}}_{PL}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)={\mathrm{\Ψ}}_{ML}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)+{\mathrm{\α}}_{B_0}\·R_{smooth\_B_0}\left(r\right)$

其中，为势函数的系数。

最小二乘法拟合计算匀场线圈的电流改变量4：将匀场线圈的场图表示为矩阵A(r,j)，待匀场的静磁场场图表示为向量b(r)，x代表匀场线圈的电流改变量，则计算匀场线圈电流改变量简化为求解线性方程组Ax＝b。可以采用奇异值分解(Singular Value Decomposition，SVD)方法得到匀场线圈的电流改变量x(j),j＝1,2,3...N。

有时基础磁场均匀性很差，信号畸变或缠绕等可能被放大，造成方程组的病态，使得最终解被污染，甚至超过一定的设定范围(越界)。Tikhonov正则化参数引入可以有效减少求解线性方程组Ax＝b的误差，虽然会丧失一定的计算精度，但是提高了数据的准确性。

设置匀场线圈电流改变量并判断是否需要迭代5：将计算的匀场线圈电流改变量x(j),j＝1,2,3...N分别加至各线圈原电流值上CurrValue(j),j＝1,2,3...N，若磁场均匀度无改善或是均匀度满足预定设置的要求，则完成所有匀场；否则将当前N个匀场线圈的电流值作为待匀场的N个匀场线圈的电流值返回步骤3，可通过采样获得实际的谱图线形¹H指标(50％/0.55％/0.11％处的线宽)来判断匀场效果。

本实施例选取极差的“冷场”条件(Z1-Z7匀场线圈的电流值均设置为0A，CurrValue(j)＝0,j＝1,2,3...7)，由于静磁场的不均匀性偏差较大，使得测量的¹H谱图在4.7ppm水峰信号底部很宽且严重偏向右侧(图4(1))，其50％/0.55％/0.11％处的线宽为45.0/1598.0/2322.0Hz。

这种条件下，本实施例选取第一次回波成像时间TE₁和第二次回波成像时间TE₂分别为0s和0.005s进行传统相位差成像方式拟合匀场线圈的场图和静磁场场图，并求解匀场线圈电流改变量。如图3(1)a～g所示，采用传统相位差成像的方式获得的匀场线圈场图有明显的噪声误差，很难获得匀场效果(图4(2))。下表3(左)表示采用Tikhonov正则化SVD方式求解的匀场线圈电流改变量x，高阶Z7线圈电流改变量仍然出现越界(超过匀场线圈的可设电流值上限或下限)的情况，且设置相应配置值进行匀场，测量的¹H谱图在4.7ppm水峰信号基本无改善，50％/0.55％/0.11％处的线宽为62.0/1551.0/2160.0Hz。

采用本正则化优化重建算法拟合三次回波成像数据(TE₁＝0s、TE₂＝0.005s和TE₃＝0.1s)不仅能获得清晰完整的匀场线圈场图(如图3(5)a～g)，而且对于较差的磁场不均匀性有很好的校正效果，采用同上的Tikhonov正则化SVD方式求解匀场线圈电流改变量。如下表3(右)所示，匀场值数据均无越界，且设置相应配置值进行匀场后，测量的¹H谱图(图4(3))在4.7ppm水峰信号底部明显变窄，对称性改善，50％/0.55％/0.11％处的线宽为58.0/694.0/1017.0Hz。

表3采用两种梯度匀场方式针Z1-Z7匀场线圈校正后的匀场电流改变量比较

^*该Z7匀场线圈的电流改变值越界，超过可设电流值下限。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (5)

1.一种基于正则化磁场分布图像重建的梯度匀场方法，其特征在于，包括以下步骤：

步骤1、根据需要进行梯度匀场的N个匀场线圈选择相应的梯度回波脉冲序列；

步骤2、对初始化匀场线圈的磁场状态执行梯度回波脉冲序列采样，根据回波采样数据初始化匀场线圈场图，并通过正则化磁场分布图像重建获得最终的匀场线圈场图；

步骤3、对待匀场的磁场状态执行梯度回波脉冲序列采样，根据回波采样数据初始化待匀场的静磁场场图，并通过正则化磁场分布图像重建获得最终的静磁场场图；

步骤4、将步骤2获得的最终的匀场线圈场图表示为矩阵A，将步骤3获得的最终的静磁场场图数据标识为向量b，求解线性方程组Ax＝b得到待匀场线圈的电流改变量；

步骤5、将计算的匀场线圈电流改变量分别加至各待匀场的匀场线圈电流值上，若磁场均匀度无改善或是均匀度满足预定设置的要求，则完成所有匀场；否则将当前N个待匀场的匀场线圈的电流值作为新的待匀场的N个匀场线圈的电流值返回步骤3。

2.根据权利要求1所述的一种基于正则化磁场分布图像重建的梯度匀场方法，其特征在于，所述的步骤2中包括以下步骤：

步骤2.1、初始化N个匀场线圈的电流值Value(j),j＝1,2,3...N，将初始化通有匀场线圈电流的磁场状态设为基础磁场状态，即磁场状态0，进行梯度回波采样，获得成像回波数据；

步骤2.2、设置所选的各个匀场线圈在采样过程中的电流变化量△Change(j)；

步骤2.3、将每个匀场线圈对应电流变化量△Change(j)依次叠加在基础磁场状态下对应的匀场线圈上，得到磁场状态1～N，并进行采样，得到叠加电流变化量后的成像回波数据；

步骤2.4、对磁场状态0～磁场状态N的成像回波数据进行傅里叶变换，得到各个磁场状态下的信号强度和相位：

信号强度：(s₁(r,i),TE₁),(s₂(r,i),TE₂),(s₃(r,i),TE₃)

相位：(φ₁(r,i),TE₁),(φ₂(r,i),TE₂),(φ₃(r,i),TE₃)

其中，信号强度s₁～s₃、相位数据φ₁～φ₃均和回波采样的时间TE₁～TE₃相关联；当i＝0时为基础磁场采样得到的信号强度和相位数据，i＝1,2,3…N时表示依次改变N个匀场线圈电流值后采样得到的信号强度和相位数据；r表示计算的有效像素点，r的个数为N_P；

步骤2.5、将磁场状态1～N状态下采样得到的信号强度和相位数据与磁场状态0采样得到的信号强度和相位数据分别做差，得到表征每个匀场线圈单独作用效果的信号强度y_k(r,j)和相位ph_k(r,j)；k＝1,2,3表示单个匀场线圈所对应的三次扫描的成像回波时间TE₁、TE₂和TE₃；

步骤2.6、将步骤2.5中获得的表征各个匀场线圈单独作用的相位中的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据ph₂和ph₁作差，得到成像的相位差并进行相位解缠△φ₂₁(r,j)＝unwrap(ph₂(r,j)-ph₁(r,j))；

步骤2.7、根据以下公式初始化各个匀场线圈的场图：ω_Shim(r,j)＝△φ₂₁(r,j)/(TE₂-TE₁)；

步骤2.8、建立匀场线圈场图的罚似然估计函数Ψ_PL(ω_Shim(r,j))，然后通过优化变换的方式最小化罚似然估计函数Ψ_PL(ω_Shim(r,j))并获得最终的匀场线圈场图数据

3.根据权利要求2所述的一种基于正则化磁场分布图像重建的梯度匀场方法，其特征在于，所述的步骤2.8中建立罚似然估计函数Ψ_PL(ω_Shim(r,j))包括以下步骤：

步骤2.8.1、拟合各匀场线圈场图的最大似然估计函数Ψ_ML(ω_Shim(r,j))；

步骤2.8.2、利用二次函数凹凸性质建立场图数据ω_Shim的势函数，其中势函数包括平滑势函数R_smooth(r,j)和先验势函数R_prior(r,j)；

步骤2.8.3、将匀场线圈场图ω_Shim(r,j)的最大似然估计函数与势函数相加，建立ω_Shim(r,j)的罚似然估计函数Ψ_PL(ω_Shim(r,j))：

Ψ_PL(ω_Shim(r,j))＝Ψ_ML(ω_Shim(r,j))+α·R_smooth(r,j)+β·R_prior(r,j),j＝1,2,3…N

其中，α、β为势函数的系数。

4.根据权利要求3所述的一种基于正则化磁场分布图像重建的梯度匀场方法，其特征在于，所述的步骤3包括以下步骤：

步骤3.1、记录待匀场的N个匀场线圈的电流值CurrValue(j)，将此时通有匀场线圈电流的磁场状态设为待匀场的磁场状态，并获得采样数据；

步骤3.2、对步骤3.1的采样数据进行傅里叶变换，获得表征待匀场的静磁场像的信号强度和相位：

信号强度～((yb0₁(r),TE₁),(yb0₂(r),TE₂),(yb0₃(r),TE₃))

相位～((phb0₁(r),TE₁),(phb0₂(r),TE₂),(phb0₃(r),TE₃))

其中，信号强度yb0₁～yb0₃、相位数据phb0₁～phb0₃均分别和回波采样的时间TE₁～TE₃相关联；r表示计算的有效像素点，r的个数为N_P；

步骤3.3、将步骤3.2获得的第二次成像回波时间TE₂和第一次成像回波时间TE₁对应得到的相位数据phb0₂和phb0₁作差，得到成像的相位差并进行相位解缠，通过相位差和时间差(TE₂-TE₁)相除获得表征待匀场的静磁场场图的频率数据

${\mathrm{\ω}}_{B_0}\left(r\right)=unwrap\left(phb{0}_2\right(r)-phb{0}_1(r\left)\right)/({\mathrm{TE}}_2-{\mathrm{TE}}_1);$

步骤3.4、建立待匀场的静磁场场图的罚似然估计函数通过优化变换的方式最小化罚似然估计函数获得最终的静磁场场图数据

5.根据权利要求4所述的一种基于正则化磁场分布图像重建的梯度匀场方法，其特征在于，所述的步骤3.4中建立罚似然估计函数包括以下步骤：

步骤3.4.1、拟合待匀场的静磁场场图的最大似然估计函数

步骤3.4.2、利用二次函数凹凸性质建立待匀场的静磁场场图的势函数，其中势函数为平滑势函数

步骤3.4.3、将静磁场场图的最大似然估计函数与势函数相加，建立的罚似然估计函数

${\mathrm{\Ψ}}_{PL}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)={\mathrm{\Ψ}}_{ML}\left({\mathrm{\ω}}_{B_0}\right(r\left)\right)+{\mathrm{\α}}_{B_0}\·R_{smooth\_B_0}\left(r\right)$

其中，为势函数的系数。