CN103675733B - 一种基于不均匀磁场拟合线形的自动搜索匀场方法 - Google Patents

Abstract

本发明公开了一种基于不均匀磁场拟合线形的自动搜索匀场方法，采用梯度回波脉冲序列测量空间磁场的不均匀像，并以此拟合表征当前磁场均匀性的虚拟谱图线形及其相关性能指标；根据匀场线圈所构建的多维空间向量，将虚拟谱图线形的相关性能指标作为得到的磁场均匀性好坏的评价标准，依据此评价标准进行多维下降单纯形法迭代搜索，并最终返回最优匀场电流。本发明在匀场过程中能充分考虑频谱的半高宽和对称性等性能指标；同时，和常规的采用梯度回波成像技术的梯度匀场相比，本匀场技术不必考虑匀场线圈场图的测量和拟合，降低了通过梯度成像来匀场的难度和限定要求；自动拟合搜索步径和收敛条件摆脱了繁琐的手动设置、提高了匀场的速度和效率。

Description

一种基于不均匀磁场拟合线形的自动搜索匀场方法

技术领域

本发明涉及核磁共振谱仪的匀场技术，更具体涉及一种基于不均匀磁场拟合线形的自动搜索匀场方法。适用于高分辨的核磁共振实验匀场。

背景技术

对于现代液体核磁共振(NMR)谱仪来说，分辨率是一个重要的指标，而影响谱图分辨率的主要因素之一就是磁场空间均匀性。由于静磁场的不均匀性引起的核磁共振信号展宽必定导致实验结果的重大误差，因此从核磁共振现象被发现以来，如何提高静磁场的均匀性（匀场）一直是一项非常重要的研究工作。

随着谱仪技术近几十年的发展，现阶段的自动匀场技术主要分为两类：自动搜索匀场技术和梯度匀场技术。自动搜索匀场方法最早出现于1968年(R.R.Ernst.Measurementandcontrolofmagneticfieldhomogeneity,Thereviewofscientificinstruments.39,1968,998-1012)，其基本思路是利用单纯形（Simplex）搜索等算法来优化氘（D,或²H）锁电平值或FID面积大小等评价标准从而补偿磁场。这就引入了两个问题：1）锁电平或是FID面积作为方便、直观的人为判断磁场好坏的依据，并不一定能达到良好的匀场效果（频谱线形标准）。事实证明，较高的氘锁电平和较大的FID面积都不一定对应于较窄的谱图峰宽；2）匀场搜索迭代的范围和收敛条件必须依靠人为设置。由于不同样品对补偿电流的敏感性不同，不合适的迭代步径和收敛条件都会使得搜索时间过长或是过短，均不利于获得好的匀场结果，因此繁琐的手动设置既增加了操作的难度也浪费了时间。

1994年Zijl基于核磁共振成像(MRI)技术的原理提出了梯度匀场方法(P.C.M.VanZijl,S.Sukumar,M.O’NeilJohnson,atal.Optimizedshimmingforhigh-resolutionNMRusingthree-dimensionalimage-basedfieldmapping,JournalofmagneticresonanceA.111,1994,203-207)，梯度匀场的最大优点就是可以通过梯度回波脉冲序列测量出磁场不均匀性的空间分布，然后基于这种空间分布进行匀场。2006年，Weiger在传统梯度匀场的基础上又更进了一步，提出采用梯度匀场后的剩磁场（ResidualField）的像拟合虚拟频谱图，通过模拟的线形指标衡量并优化当前磁场(MarkusWeiger,ThomasSpeck,MichaelFey.Gradientshimmingwithspectrumoptimisation,Journalofmagneticresonance.182,2006,38-48)。通常情况，梯度匀场必须同时获得磁场不均匀性的像和所用匀场线圈的单位电流作用的磁场影响（匀场线圈场图）。匀场线圈场图的制作过程稍显复杂且耗费较多时间。同时，对匀场线圈的空间影响进行成像必须在一个较均匀的磁场环境，均匀性较差的磁场容易使得匀场线圈场图测量数据的畸变从而影响拟合。

本专利采用一种基于不均匀磁场拟合线形的自动搜索匀场方法，该技术能自动拟合搜索迭代的步径和收敛条件，无需调试和校验，节省了手动设置的时间；另一方面由于通过测量空间磁场的像来拟合谱图线形，使得优化过程充分考虑频谱的半高宽和对称性等性能指标。再则，和常规的采用梯度回波成像技术的梯度匀场相比，此匀场方式可不必考虑匀场线圈场图的测量和拟合，加快了匀场速度。

发明内容

本发明的目的是提供了一种基于不均匀磁场拟合线形的自动搜索匀场方法，为了实现上述目的，本发明采用以下技术方案：

一种基于不均匀磁场拟合线形的自动搜索匀场方法，包括以下步骤：

步骤1、选择匀场线圈组：选择需要进行室温匀场的匀场线圈组；

步骤2、拟合匀场线圈的搜索步径和收敛条件：利用谱图吸收线形的二次矩大小I和匀场线圈的电流增量x的二次函数关系I＝Ax²+Bx+C，其中A、B、C为二次函数系数，拟合各匀场线圈的搜索步径d_j和收敛条件r_j；

步骤3、拟合表征当前磁场不均匀性的频谱图及线形指标：利用梯度回波技术获得不均匀磁场的图像数据，根据该图像数据拟合出表征当前磁场不均匀性的吸收线形谱图F，即虚拟线形，并最终计算出该虚拟线形的相关性能指标Crit；

步骤4、采用多维下降单纯形法迭代搜索：利用N个匀场线圈构建N维空间[X₁,X₂,...,X_N]，将谱图吸收线形的相关性能指标Crit作为得到的磁场均匀性好坏的评价标准Q，并利用Crit越小则磁场越均匀的性质进行空间点的多维下降单纯形法迭代搜索；

步骤5、迭代结束返回匀场电流：多维下降单纯形法迭代结束，最终得到最大评价值Q并返回最优匀场电流。

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

步骤2.1、选定某一匀场线圈a，a=j，j∈{1～N}，设定其匀场电流增量x_a=0并直接采用单脉冲序列进行采样，从而得到原始数据；

步骤2.2、利用傅里叶变换和相位校正处理采集到的原始数据，得到吸收线形的谱图数据；

步骤2.3、计算当前谱图的吸收线形二次矩的大小I₀，依据公式：

$I=\frac{\∫\mathrm{}S\left(\mathrm{\Ω}\right){\mathrm{\Ω}}^2\mathrm{d\Ω}}{\∫S\left(\mathrm{\Ω}\right)\mathrm{d\Ω}}-{\left(\frac{\∫S\left(\mathrm{\Ω}\right)\mathrm{\Ωd\Ω}}{\∫S\left(\mathrm{\Ω}\right)\mathrm{d\Ω}}\right)}^2$

其中，I为吸收线形的二次矩大小，S(Ω)是谱峰信号，Ω为谱图频率坐标；

步骤2.4、重复步骤2.1～2.3分别计算预设匀场线圈电流增量为x_a=h和x_a=-h时的吸收线形二次矩I₊和I_-；

步骤2.5、根据（0,I₀）、（h,I₊）、（_-h,I_-）三点拟合二次矩大小I和匀场线圈电流增量x的二次函数关系I＝Ax²+Bx+C，求得系数A和B，A＝(I₊+I_--2I₀)(2h²)，B＝(I₊-I_-)(2h)；

步骤2.6、计算当前匀场线圈a的搜索步径d_a和收敛条件r_a，依据公式：

$d_a=\frac{-B}{\left|B\right|}\·\sqrt{\frac{I_0}{A}}$

$r_a=\left|\frac{B}{2A}\right|$

其中，系数A和B的大小由步骤2.5求得，I₀为匀场线圈电流增量为x_a=0时谱图的吸收线形二次矩的大小；

步骤2.7、重复步骤2.1～2.6拟合其他匀场线圈的搜索步径d_j和收敛条件r_j。

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

步骤3.1、根据步骤1所选的室温匀场线圈组设定梯度回波的脉冲序列；

步骤3.2、根据样品差异调整射频脉冲产生的翻转角α；

步骤3.3、执行脉冲序列进行梯度回波采样并获得采样数据；

步骤3.4、对所有采样数据进行傅里叶变换获得表征磁场像的相位数据(r)，其中r表示为采样数据点的空间位置坐标，(r)和脉冲序列中两次的成像回波时间有关，对应于( ₁(r),TE1),( ₂(r),TE2)；同时设置相位数据的有效范围；

步骤3.5、计算步骤3.4中设置的有效范围内各点的两次成像相位差Δ _o(r)，并对相位差Δ _o(r)进行相位解缠，从而获得表征磁场均匀性的像数据

B₀(r)＝unwrap(Δ _o(r))，其中unwrap表示对相位数据进行解缠操作；

步骤3.6、根据步骤3.4的第一次成像回波时间TE1求得的相位数据 ₁(r)拟合表征当前射频场区域均匀性的像数据B₁(r)，

步骤3.7、拟合表征当前磁场均匀性的像数据B₀(r)的频谱加权统计直方图His(Ω)；

步骤3.8、采用洛伦兹线形L(Ω)和步骤3.7拟合的频谱统计直方图His(Ω)直接作卷积计算，获得表征当前磁场均匀性的虚拟线形F(Ω)；

步骤3.9、计算虚拟线形F(Ω)的相关性能指标Crit。

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

步骤3.6.1、计算相位数据的绝对值谱数据Abs（ ₁(r)），并获得绝对值谱最大的强度Max(Abs( ₁(r)))；

步骤3.6.2、计算射频场均匀性的像数据B₁(r)＝Abs( ₁(r))/Max(Abs( ₁(r)))。

如上所述的步骤3.7的频谱加权统计直方图His(Ω)是基于以下公式：

$\mathrm{His}\left(\mathrm{\Ω}\right)={\∫}_r\mathrm{\δ}(\mathrm{\Ω}-\mathrm{\γ}\·B_0\left(r\right))\·w\left(r\right)\mathrm{dr}$

式中，Ω为谱图频率坐标，δ为狄拉克函数，γ为采样核的旋磁比，离散分布各点的权重影响w(r)由以下公式获得：

w(r)＝sin(B₁(r)·α)·B₁(r)

式中，α为射频脉冲产生的翻转角，B₁(r)为射频场均匀性的像数据。

如上所述的步骤3.8中洛伦兹线形L(Ω)由下式得到：

$L\left(\mathrm{\Ω}\right)=\frac{\mathrm{\λ}}{\mathrm{\π}}\·\frac{1}{{\mathrm{\λ}}^2+(\mathrm{\Ω}-{\mathrm{\Ω}}_0{)}^2}$

式中，Ω为谱图频率坐标，理想情况下洛伦兹线形可以表示为关于Ω₀＝0轴对称，且线形的半高宽λ＝1/(π·T₂)，T₂为样品的自旋-自旋弛豫时间。

如上所述的步骤3.9中的Crit为Crit_symm或Crit_half或Crit_env，

其中Crit_symm由以下步骤获得：

比较虚拟线形F(Ω)最大峰强度点左右两边频率对称位置的强度差求得虚拟线形的对称性Crit_symm：

${\mathrm{Crit}}_{\mathrm{symm}}=\frac{\mathrm{\Σ}|F({\mathrm{\Ω}}_{\max }+\mathrm{\Ω})-F({\mathrm{\Ω}}_{\max }-\mathrm{\Ω})|}{\mathrm{\ΣF}\left(\mathrm{\Ω}\right)}$

式中，Ω_max表示为最大峰强度对应的谱图频率坐标；

其中Crit_half和Crit_env由以下步骤获得：

选择典型洛伦兹线形作为包络线E(Ω)，将虚拟线形F(Ω)包络在其内，关于包络线的求解可通过最小化惩罚函数P_func获得：

$P_{\mathrm{func}}=\mathrm{\&Sigma;}{P}_{\mathrm{func}}\left(\mathrm{\&Omega;}\right)=\left\{\begin{array}{cc}{\left(\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\right)}^2,& \mathrm{Diff}\left(\mathrm{\&Omega;}\right)\&GreaterEqual;0\\ {(200\&CenterDot;|\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\left|\right)}^2,& \mathrm{Diff}\left(\mathrm{\&Omega;}\right)<0\end{array}\right.$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线和虚拟线形的差；

采用包络法评价磁场均匀性直接选择计算包络线E(Ω)的半高宽Crit_half；

计算包络线和虚拟线形的包络误差Crit_env：

${\mathrm{Crit}}_{\mathrm{env}}=\frac{\mathrm{\&Sigma;}\left|\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\right|}{\mathrm{\&Sigma;F}\left(\mathrm{\&Omega;}\right)}$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线和虚拟线形的差。

如上所述的步骤4中的多维下降单纯形法的迭代搜索包括以下步骤：

步骤4.1、利用当前磁场状态下N个匀场线圈的电流值和步骤2求得的各匀场线圈的搜索步径d_j建立(N+1)×N的空间矩阵M(x_ij)，具体如下所示：

其中，M(x_ij)中单一矩阵元素x_ij代表各匀场线圈具有的电流大小；矩阵的列编号依次对应于所选的N个匀场线圈，i∈{0～N}，j∈{1～N}；矩阵中每个行向量(x_i1,x_i2,...,x_ij,...,x_iN)，均构成一个空间点，则初始化建立的N+1个空间点表示通有不同匀场电流后的N+1个磁场状态：行编号0表示无匀场线圈电流变化的初始磁场状态，行编号1～N分别对应于第1～N个线圈增加电流增量大小为d_j后的磁场状态；

步骤4.2、利用虚拟线形的相关性能指标Crit越小则磁场越均匀的性质，将Crit作为搜索得到磁场均匀性好坏的评价标准Q，Q=1/Crit，则对于具有不同匀场电流的空间各点(x_i1,x_i2,...,x_ij,...,x_iN),i∈{0～N}，j∈{1～N}，均有相应的评价标准Q(i)；

步骤4.3、根据步骤4.1初始化N+1个空间点并拟合出各点的评价标准大小Q(i);

步骤4.4、依据Q的大小排序N+1个点，其中数值最小的点为W，数值最大的点为B，中间点设为H；当N+1为奇数时有1个H点，当N+1为偶数时有2个H点，分别表示为H₁和H₂两点，两者大小比较为Q(H₁)≥Q(H₂)；

步骤4.5、拟合除W外的所有N个空间点的中心，设为M，并产生由W指向M的向量；

步骤4.6、沿W至M向量方向拟合W的反射点R，R视为W关于M呈中心对称的点，并计算其评价值Q(R)；

步骤4.7、比较Q(R)和Q(B)的大小，若Q(R)>Q(B)，则继续沿方向扩展至反射扩展点RE，RE视为方向上M关于R的反射点，并评价Q(RE)，进入步骤4.8；若Q(R)≤Q(B)，则跳至步骤4.9；

步骤4.8、比较Q(RE)和Q(R)的大小，两者较大的点取代W点，并跳至步骤4.13；

步骤4.9、比较Q(R)和Q(H)的大小，根据H点的个数分为以下两种情况：

当N+1为奇数时，若Q(R)>Q(H)，则将R取代W，并跳至步骤3.13，若Q(R)≤Q(H)，则跳至步骤3.10;

当N+1为偶数时，若Q(R)>Q(H₂)时，则将R取代W，并跳至步骤3.13，若Q(R)≤Q(H₂)，则跳至步骤3.10;

步骤4.10、比较Q(R)和Q(W)的大小，若Q(R)>Q(W)，则将R取代W，并跳至步骤4.13；若Q(R)≤Q(W)，取M和W的中点C，即收缩点C，并求解其评价值Q(C)，进入步骤4.11;

步骤4.11、比较Q(C)和Q(W)的大小，若Q(C)>Q(W)，则C取代W，并跳至步骤4.13；若Q(C)≤Q(W)，则迭代无更新点，跳至步骤4.12；

步骤4.12、由于迭代无更新点，选择步骤4.4中数据空间的评价值第二差的点W+1作为W，并返回步骤4.5重新开始；若仍无更新点，则继续选择点W+2、W+3，...直至B作为W返回操作步骤4.5；若采用最大的点B开始迭代仍无更新点，则截取除点B外的所有N个点和B的中点，设为全收缩点FC_j，j∈{1～N}，并使用各全收缩点取代原来的N个点，跳至步骤4.13；

步骤4.13、判断迭代是否收敛，其标准为，更新空间点后的矩阵中，任一列向量的任意两个矩阵元素的差均应小于等于步骤2.6中获得的收敛终止条件rj，表示如下：

若收敛判断为是，则终止匀场迭代，并返回当前最优Q(B)对应的匀场线圈电流大小；若收敛判断为否，则返回步骤4.4，继续迭代搜索。

如上所述的步骤2中所有频率Ω的积分范围为不低于最大谱峰高度的15%的频率区域；步骤2中匀场线圈预设电流增量h范围为5～50mA。

如上所述的步骤3.4中相位数据的有效范围为不低于采样数据的绝对值谱最大谱峰高度的25%的采样数据区域。

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

1、可以自动拟合搜索算法的初始迭代步径和收敛条件。现有的自动搜索匀场技术必须依靠实验操作者对搜索的迭代步径进行调试，线圈电流的迭代步径通常为针对不同样品的氘锁电平值或是FID面积大小的5%～15%敏感变化，而收敛条件则为0.5%～1%。本发明克服了繁琐的人为干预，直接通过拟合谱图吸收线形的二次矩和匀场线圈电流的函数关系进行自动化设置。

2、可以有效考虑谱图的线形来进行匀场。频谱的半高宽和对称性等是判断匀场效果最本质的指标，因此根据磁场不均匀性的像拟合的虚拟线形可以作为判断磁场好坏的直接依据，避免了采用锁电平值或FID面积大小匀场时出现的判断误差。

3、避免了制作匀场线圈场图的繁琐步骤。由于梯度匀场必须依靠制作匀场线圈场图来配套进行匀场，因此对基础磁场提出较高的均匀度要求，同时需要花费大量的制作场图的时间，降低了匀场的效率。本匀场技术方案采用直接获取梯度回波的像来拟合评价磁场均匀性的指标，一方面借助和梯度匀场相同的脉冲序列实现方式来有效获取当前磁场状况，另一方面又避免了复杂的匀场线圈的场图拟合，节省了匀场时间。

附图说明

图1是一种基于不均匀磁场拟合线形的自动搜索匀场方法流程图。

其中选择匀场线圈组1、拟合匀场线圈的搜索步径和收敛条件2、拟合表征当前磁场不均匀性的频谱图及线形指标3、采用多维下降单纯形法（Simplex）迭代搜索4、迭代结束返回匀场电流5。

图2是拟合表征当前磁场不均匀性的频谱图及线形指标的流程图。

其中通过梯度回波测得当前磁场不均匀性的像数据3.1～3.5、计算射频场均匀性的像数据3.6、拟合当前磁场不均匀性的频谱加权统计直方图3.7、根据理想状况洛伦兹线卷积拟合虚拟频谱线形3.8、计算线形的相关性能指标3.9。

图3是一种基于不均匀磁场拟合线形的自动搜索匀场方法中多维下降单纯形法（simplex）。

图4是多维下降单纯形法中提到的相关操作。

以匀场线圈数N为6、空间点个数N+1=7为例，分别为W(Worst)、W+1、W+2、H(Half)、B-2、B-1和B(Best)；M是除W点外所有点的中心点；向量的方向表示单纯形法操作方向，图示多维下降单纯形法包括：反射(Reflect,R)、反射扩展(ReflectExpand,RE)、收缩(Contract,C)和全收缩(FullContract,FC)操作。

图5（a）、图5（b）、图5（c）和图5（d）是本匀场方法针对1H谱的匀场过程中的效果对比图。匀场分为两步，分别采用加弛豫试剂的水溶液（0.1mg/mLGdCl₃，1%H₂O，99%D₂O）和1%氘代氯仿的丙酮溶液（1%chloroform，Acetone-d₆）进行粗调优化和细调优化。

图5（a）为采用加弛豫试剂的水溶液（0.1mg/mLGdCl₃，1%H₂O，99%D₂O）粗调优化前；

图5（b）为采用加弛豫试剂的水溶液（0.1mg/mLGdCl₃，1%H₂O，99%D₂O）粗调优化后；

图5（c）为1%氘代氯仿的丙酮溶液（1%chloroform，Acetone-d₆）细调优化前；

图5（d）为1%氘代氯仿的丙酮溶液（1%chloroform，Acetone-d₆）细调优化后。

具体实施方式

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

实施例1：

如图1流程图所示，一种基于不均匀磁场拟合线形的自动搜索匀场方法，包括以下步骤：

步骤1、选择匀场线圈组：选择需要进行室温匀场的匀场线圈组（j=1～N，N为所选匀场线圈的个数，例如:Z1/Z2/Z3/Z4/Z5/Z6，N=6）；

步骤2、拟合匀场线圈的搜索步径和收敛条件：利用谱图吸收线形的二次矩大小I和匀场线圈的电流增量x的二次函数关系I＝Ax²+Bx+C（其中A、B、C为二次函数系数）拟合各匀场线圈的搜索步径d_j和收敛条件r_j包括以下步骤：

步骤2.1、选定某一匀场线圈a（a=j，j∈{1～N}），设定其匀场电流增量x_a=0并直接采用单脉冲序列进行采样，从而得到原始数据；

步骤2.2、利用傅里叶变换和相位校正处理采集到的原始数据，得到吸收线形的谱图数据；

步骤2.3、计算当前谱图的吸收线形二次矩的大小I₀，依据公式：

$I=\frac{\&Integral;\mathrm{}S\left(\mathrm{\&Omega;}\right){\mathrm{\&Omega;}}^2\mathrm{d\&Omega;}}{\&Integral;S\left(\mathrm{\&Omega;}\right)\mathrm{d\&Omega;}}-{\left(\frac{\&Integral;S\left(\mathrm{\&Omega;}\right)\mathrm{\&Omega;d\&Omega;}}{\&Integral;S\left(\mathrm{\&Omega;}\right)\mathrm{d\&Omega;}}\right)}^2$

其中，I为吸收线形的二次矩大小，S(Ω)是谱峰信号，Ω为谱图频率坐标；

步骤2.4、重复步骤2.1～2.3分别计算预设匀场线圈电流增量为x_a=h和x_a=_-h时的吸收线形二次矩I₊和I_-；

步骤2.5、根据（0,I₀）、（h,I₊）、（_-h,I_-）三点拟合二次矩大小I和匀场线圈电流增量x的二次函数关系I＝Ax²+Bx+C，求得系数A和B，A＝(I₊+I_--2I₀)(2h²)，B＝(I₊-I_-)/(2h)；

步骤2.6、计算当前匀场线圈a的搜索步径d_a和收敛条件r_a，依据公式：

$d_a=\frac{-B}{\left|B\right|}\&CenterDot;\sqrt{\frac{I_0}{A}}$

$r_a=\left|\frac{B}{2A}\right|$

其中，系数A和B的大小由步骤2.5求得，I₀为匀场线圈电流增量为x_a=0时谱图的吸收线形二次矩的大小；

步骤2.7、重复步骤2.1～2.6拟合其他匀场线圈的搜索步径d_j和收敛条件r_j。

步骤3、拟合表征当前磁场不均匀性的频谱图及线形指标：利用梯度回波技术获得不均匀磁场的图像数据，根据该图像数据拟合出表征当前磁场不均匀性的吸收线形谱图F（虚拟线形），并最终计算出该虚拟线形的相关性能指标Crit，其中如附图2所示包括以下步骤：

步骤3.1、根据步骤1所选的室温匀场线圈组设定梯度回波的脉冲序列，分为以下两种情况：1）当匀场组只含有Z轴方向（轴向或On-Axis）的匀场线圈（如Z1～Z7），则采用一维梯度回波（PFGSTE）脉冲序列并需要两次成像采样；2）当匀场组含有非Z轴（径向或Off-Axis）的匀场线圈（如X，ZX，X2-Y2，XY，X3等）采用倾斜三维梯度回波（Ramped3DGRE）脉冲序列，此时X和Y方向具有相位编码，相位编码数均设定为N_X=N_Y=4；

读入脉冲序列默认参数：一维梯度回波（PFGSTE）脉冲序列和倾斜三维梯度回波（Ramped3DGRE）脉冲序列的参数设置分别如下表1和表2所示；

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

	氘（2H）核	氢（1H）核
			采样谱宽sw（kHz）	20	40
采样时间acq（ms）	12.8	12.8
			第一次成像回波时间TE1（ms）	50	5
第二次成像回波时间TE2（ms）	200	20

表2倾斜三维梯度回波（Ramped3DGRE）脉冲序列参数设置

步骤3.2、根据样品差异调整参数：根据样品中氘代含量的差别设置累加次数ns=1～4；调节射频脉冲产生的翻转角α、调节接收机增益（确保信号强度大且不溢出）；调节脉冲前等待时间d1，确保每次梯度回波信号形状无变形且强度无减弱；调节脉冲序列中两次的成像回波时间（包括第一次成像回波时间TE1和第二次成像回波时间TE2），使得前后两次回波信号强度比满足大于2/3等。

步骤3.3、执行脉冲序列进行梯度回波采样并获得采样数据：当采用一维梯度回波脉冲序列时，需要进行两次回波采集；当执行倾斜三维梯度回波脉冲序列时，还需要进行X与Y方向的相位编码N_X×N_Y，则总回波数据采集为2×N_X×N_Y；

步骤3.4、对所有采样数据进行傅里叶变换获得表征磁场像的相位数据φ(r)（r表示为采样数据点的空间位置坐标，φ(r)和脉冲序列中两次的成像回波时间有关，对应于( ₁(r),TE1),( ₂(r),TE2)）；同时设置相位数据的有效范围；

步骤3.5、计算步骤3.4中设置的有效范围内各点的两次成像相位差Δ _o(r)，并对相位差Δ _o(r)进行相位解缠，从而获得表征磁场均匀性的像数据B₀(r)＝unwrap(Δ _o(r))（其中unwrap表示对相位数据进行解缠操作）；

步骤3.6、根据步骤3.4的第一次成像回波时间TE1求得的相位数据 ₁(r)拟合表征当前射频场区域均匀性的像数据B₁(r)，具体步骤如下：

1）计算相位数据的绝对值谱数据Abs( ₁(r))，并获得绝对值谱最大的强度Max(Abs( ₁(r)));

2）计算射频场均匀性的像数据B₁(r)＝Abs( ₁(r))/Max(Abs( ₁(r)))；

步骤3.7、拟合表征当前磁场均匀性的像数据B₀(r)的频谱加权统计直方图His(Ω)：

His(Ω)＝∫_rδ(Ω-γ·B₀(r))·w(r)dr

式中，Ω为谱图频率坐标，δ为狄拉克函数，γ为采样核的旋磁比，离散分布各点的权重影响w(r)可由以下公式获得：

w(r)＝sin(B₁(r)·α)·B₁(r)

式中，α为射频脉冲产生的翻转角，B₁(r)为射频场均匀性的像数据；

步骤3.8、采用理想情况的洛伦兹线形L(Ω)和步骤3.7拟合的频谱统计直方图His(Ω)直接作卷积计算，获得表征当前磁场均匀性的虚拟线形F(Ω)；

步骤3.9、计算虚拟线形F(Ω)的相关性能指标Crit，主要包括如下：

1）虚拟线形的对称性

比较虚拟线形F(Ω)最大峰强度点左右两边频率对称位置的强度差求得虚拟线形的对称性Crit_symm：

${\mathrm{Crit}}_{\mathrm{symm}}=\frac{\mathrm{\&Sigma;}|F({\mathrm{\&Omega;}}_{\max }+\mathrm{\&Omega;})-F({\mathrm{\&Omega;}}_{\mathrm{nax}}-\mathrm{\&Omega;})|}{\mathrm{\&Sigma;F}\left(\mathrm{\&Omega;}\right)}$

式中，Ω_max表示为最大峰强度对应的谱图频率坐标；

2）频谱包络（envelope）法

包络法一般选择典型洛伦兹线形作为包络线E(Ω)，将目标波形（虚拟线形F(Ω)）包络在其内，关于包络线的求解可通过最小化惩罚函数P_func获得：

$P_{\mathrm{func}}=\mathrm{\&Sigma;}{P}_{\mathrm{func}}\left(\mathrm{\&Omega;}\right)=\left\{\begin{array}{c}{\left(\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\right)}^2,\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\&GreaterEqual;0\\ {(200\&CenterDot;|\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\left|\right)}^2,\mathrm{Diff}\left(\mathrm{\&Omega;}\right)<0\end{array}\right.$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线和虚拟线形的差；

采用包络法评价磁场均匀性可以直接选择计算包络线E(Ω)的半高宽Crit_half；也可以计算包络线和虚拟线形的包络误差Crit_env：

${\mathrm{Crit}}_{\mathrm{env}}=\frac{\mathrm{\&Sigma;}\left|\mathrm{Diff}\left(\mathrm{\&Omega;}\right)\right|}{\mathrm{\&Sigma;F}\left(\mathrm{\&Omega;}\right)}$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线和虚拟线形的差；

利用虚拟线形F(Ω)的相关性能指标Crit作为搜索评价指标时，可以在Crit_symm、Crit_env和Crit_half中根据谱图线形要求选择其中一个；

步骤4、采用多维下降单纯形法（Simplex）迭代搜索：利用N个匀场线圈构建N维空间[X₁,X₂,...,X_N]，将谱图吸收线形的相关性能指标Crit作为得到的磁场均匀性好坏的评价标准Q，并利用Crit越小则磁场越均匀的性质进行空间点的多维下降单纯形法（Simplex）迭代搜索，其中多维下降单纯形法的迭代搜索如附图3所示包括以下步骤：

步骤4.1、利用当前磁场状态下N个匀场线圈的电流值(c₀₁,c₀₂,...,c_0j,...,c_0N)和步骤2求得的各匀场线圈的搜索步径d_j建立(N+1)×N的空间矩阵M(x_ij)，具体如下所示：

其中，M(x_ij)中单一矩阵元素x_ij代表各匀场线圈具有的电流大小；矩阵的列编号依次对应于所选的N个匀场线圈；矩阵中每个行向量(x_i1,x_i2,...,x_ij,...,x_iN)，i∈{0～N}均构成一个空间点，则初始化建立的N+1个空间点表示通有不同匀场电流后的N+1个磁场状态：行编号0表示无匀场线圈电流变化的初始磁场状态，行编号1～N分别对应于第1～N个线圈增加电流增量大小为d_j，后的磁场状态，j∈{1～N}；

步骤4.2、利用虚拟线形的相关性能指标Crit越小则磁场越均匀的性质，选择步骤3.9中的一种指标（如包络线E(Ω)的半高宽Crit_half）作为搜索得到磁场均匀性好坏的评价标准Q（这里令Q=1/Crit），则对于具有不同匀场电流的空间各点(x_i1,x_i2,...,x_ij,...,x_iN),i∈{0～N}，j∈{1～N}，均有相应的评价标准Q(i)；

步骤4.3、根据步骤4.1初始化N+1个空间点并拟合出各点的评价标准大小Q(i);

步骤4.4、依据Q的大小排序N+1个点，其中数值最小的点为W，数值最大的点为B，中间点设为H；当N+1为奇数时有1个H点，当N+1为偶数时有2个H点，分别表示为H₁和H₂两点，两者大小比较为Q(H₁)≥Q(H₂)；

步骤4.5、拟合除W外的所有N个空间点的中心，设为M，并产生由W指向M的向量；

步骤4.6、沿W至M向量方向拟合W的反射点R（R视为W关于M呈中心对称的点），并计算其评价值Q(R)；

步骤4.7、比较Q(R)和Q(B)的大小，若Q(R)>Q(B)，则继续沿方向扩展至反射扩展点RE（RE视为方向上M关于R的反射点），并评价Q(RE)，进入步骤4.8；若Q(R)≤Q(B)，则跳至步骤4.9；

步骤4.8、比较Q(RE)和Q(R)的大小，两者较大的点取代W点，并跳至步骤4.13；

步骤4.9、比较Q(R)和Q(H)的大小，根据H点的个数分为以下两种情况：

当N+1为奇数时，若Q(R)>Q(H)，则将R取代W，并跳至步骤3.13，若Q(R)≤Q(H)，则跳至步骤3.10;

当N+1为偶数时，若Q(R)>Q(H₂)时，则将R取代W，并跳至步骤3.13，若Q(R)≤Q(H₂)，则跳至步骤3.10;

步骤4.10、比较Q(R)和Q(W)的大小，若Q(R)>Q(W)，则将R取代W，并跳至步骤4.13；若Q(R)≤Q(W)，取M和W的中点C（收缩点），并求解其评价值Q(C)，进入步骤4.11;

步骤4.11、比较Q(C)和Q(W)的大小，若Q(C)>Q(W)，则C取代W，并跳至步骤4.13；若Q(C)≤Q(W)，则迭代无更新点，跳至步骤4.12；

步骤4.12、由于迭代无更新点，选择步骤4.4中数据空间的评价值第二差的点W+1作为W，并返回步骤4.5重新开始；若仍无更新点，则继续选择点W+2、W+3，...直至B作为W返回操作步骤4.5；若采用最大的点B开始迭代仍无更新点，则截取除点B外的所有N个点和B的中点，设为全收缩点FC_j，j∈{1～N}，并使用各全收缩点取代原来的N个点，跳至步骤4.13；

步骤4.13、判断迭代是否收敛，其标准为，更新空间点后的矩阵中，任一列向量的任意两个矩阵元素的差均应小于等于步骤2.6中获得的收敛终止条件r_j，表示如下：

若收敛判断为是，则终止匀场迭代，并返回当前最优Q(B)对应的匀场线圈电流大小(x_B1,x_B2,...,x_Bj,...,x_BN)，(B或Best表示最优评价标准的空间点编号)；若收敛判断为否，则返回步骤4.4，继续迭代搜索。

步骤5、迭代结束返回匀场电流：Simplex匀场迭代结束，最终得到最大评价值Q并返回最优匀场电流。

步骤2中所有频率Ω的积分范围为不低于最大谱峰高度的15%的频率区域。

步骤2中匀场线圈预设电流增量h范围为5～50mA。

如上所述的步骤3.4中相位数据的有效范围为不低于采样数据的绝对值谱最大谱峰高度的25%的采样数据区域。

如上所述的步骤3.8中洛伦兹线形L(Ω)由下式得到：

$L\left(\mathrm{\&Omega;}\right)=\frac{\mathrm{\&lambda;}}{\mathrm{\&pi;}}\&CenterDot;\frac{1}{{\mathrm{\&lambda;}}^2+{(\mathrm{\&Omega;}-{\mathrm{\&pi;}}_0)}^2}$

式中，Ω为谱图频率坐标，理想情况下洛伦兹线形可以表示为关于Ω₀＝0轴对称，且线形的半高宽λ＝1/(π·T₂)，T₂为样品的自旋-自旋弛豫时间；

如图5所示，是采用本匀场方法针对1H谱的匀场过程中的效果对比图，匀场分为以下两阶段：

第一阶段采用加弛豫试剂的水溶液（0.1mg/mLGdCl₃，1%H₂O，99%D₂O）进行粗调匀场，初始匀场线圈电流全设为0A（冷场）：图5a表示匀场前的测试结果，4.7ppm峰为重水中的¹H信号峰，其50%/0.55%/0.11%处的线宽为252.93/457.89/462.18（Hz），线形较宽，谱峰严重下凹且严重不对称；图5b表示采用本方法匀场后的测试结果：4.7ppm峰为重水中的¹H信号峰，其50%/0.55%/0.11%处的线宽为9.17/104.23/195.38（Hz），线宽较窄且对称性良好。

本阶段匀场速度较快，采用以下步骤对27组匀场线圈进行优化：

1）优化轴向匀场线圈组Z1\Z2\Z3\Z4\Z5\Z6\Z7的电流；

2）优化径向低阶匀场线圈组Z1\Z2\X\Y\ZX\ZY\XY\X2-Y2的电流；

3）优化径向匀场线圈组X\Y\ZX\ZY\Z2X\Z2Y\Z3X\Z3Y\Z4X\Z4Y的电流；

4）优化径向匀场线圈组XY/X2-Y2/ZXY/Z(X2-Y2)/Z2XY/Z2(X2-Y2)的电流；

5)优化径向匀场线圈组X/Y/X3/Y3/ZX3/ZY3的电流。

第二阶段替换样品采用1%氘代氯仿的丙酮溶液（1%chloroform，Acetone-d₆）进行细调匀场：图5c表示匀场前的测试结果，8.05ppm峰为¹H氘代氯仿信号峰，其50%/0.55%/0.11%处的线宽为4.43/31.72/34.53（Hz），谱峰底部隆起有宽包且偏向左侧不对称；图5d表示采用本方法匀场后的测试结果：8.05ppm峰为¹H氘代氯仿信号峰，其50%/0.55%/0.11%处的线宽为0.39/6.16/9.92（Hz），线宽很窄且对称性好，达到探头的性能指标（≤0.80/8.00/15.00Hz）。

本阶段由于细调优化匀场，只针对轴向及低阶线圈进行匀场：

1）优化轴向匀场线圈组Z1\Z2\Z3\Z4\Z5\Z6的电流；

2）优化径向低阶匀场线圈组Z1\Z2\X\Y\ZX\ZY\XY\X2-Y2的电流。

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

Claims (10)

1.一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，包括以下步骤：

步骤1、选择室温匀场线圈组：选择需要进行室温匀场的室温匀场线圈组；

步骤2、拟合匀场线圈的搜索步径和收敛条件：利用谱图吸收线形的二次矩大小I和匀场线圈的电流增量x的二次函数关系I＝Ax²+Bx+C，其中A、B、C为二次函数系数，拟合各匀场线圈的搜索步径d_j和收敛条件r_j；

步骤3、拟合表征当前磁场不均匀性的频谱图及线形指标：利用梯度回波技术获得不均匀磁场的图像数据，根据该图像数据拟合出表征当前磁场不均匀性的虚拟线形F(Ω)，并最终计算出该虚拟线形F(Ω)的相关性能指标Crit；

步骤4、采用多维下降单纯形法迭代搜索：利用N个匀场线圈构建N维空间[X₁,X₂,...,X_N]，将虚拟线形F(Ω)的相关性能指标Crit作为得到的磁场均匀性好坏的评价标准Q，并利用相关性能指标Crit越小则磁场越均匀的性质进行空间点的多维下降单纯形法迭代搜索；

步骤5、迭代结束返回匀场电流：多维下降单纯形法迭代结束，最终得到最大评价值并返回最优匀场电流。

2.根据权利要求1所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤2包括以下步骤：

步骤2.1、选定某一匀场线圈a，a＝j，j∈{1～N}，设定匀场线圈a的电流增量x_a＝0并直接采用单脉冲序列进行采样，从而得到原始数据；

步骤2.2、利用傅里叶变换和相位校正处理采集到的原始数据，得到谱图吸收线形；

步骤2.3、计算当前谱图吸收线形的二次矩大小I₀，依据公式：

$I_0=\frac{\&Integral;S\left(\mathrm{\&Omega;}\right){\mathrm{\&Omega;}}^{2}d\mathrm{\&Omega;}}{\&Integral;S\left(\mathrm{\&Omega;}\right)d\mathrm{\&Omega;}}-{\left(\frac{\&Integral;S\left(\mathrm{\&Omega;}\right)\mathrm{\&Omega;}d\mathrm{\&Omega;}}{\&Integral;S\left(\mathrm{\&Omega;}\right)d\mathrm{\&Omega;}}\right)}^2$

其中，I₀为谱图吸收线形的二次矩大小，S(Ω)是谱峰信号，Ω为谱图频率坐标；

步骤2.4、重复步骤2.1～2.3分别计算预设匀场线圈a的电流增量为x_a＝h和x_a＝-h时的谱图吸收线形的二次矩大小I₊和I_-；

步骤2.5、根据(0,I₀)、(h,I₊)、(-h,I_-)三点拟合谱图吸收线形的二次矩大小I和匀场线圈a的电流增量x_a的二次函数关系I＝Ax_a ²+Bx_a+C，求得系数A和B，A＝(I₊+I_--2I₀)/(2h²)，B＝(I₊-I_-)/(2h)；

步骤2.6、计算当前匀场线圈a的搜索步径d_a和收敛条件r_a，依据公式：

$d_a=\frac{-B}{\left|B\right|}\&CenterDot;\sqrt{\frac{I_0}{A}}$

$r_a=\left|\frac{B}{2A}\right|$

其中，系数A和B的大小由步骤2.5求得，I₀为匀场线圈a的电流增量为x_a＝0时谱图吸收线形的二次矩大小；

步骤2.7、重复步骤2.1～2.6拟合其他匀场线圈的搜索步径d_j和收敛条件r_j。

3.根据权利要求1所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3包括以下步骤：

步骤3.1、根据步骤1所选的室温匀场线圈组设定梯度回波的脉冲序列；

步骤3.2、根据样品差异调整射频脉冲产生的翻转角α；

步骤3.3、执行脉冲序列进行梯度回波采样并获得采样数据；

步骤3.4、对所有采样数据进行傅里叶变换获得表征磁场像的相位数据其中r表示为采样数据点的空间位置坐标，相位数据和脉冲序列中两次的成像回波时间有关，对应于同时设置相位数据的有效范围；

步骤3.5、计算步骤3.4中设置的有效范围内各点的两次成像的相位差并对相位差进行相位解缠，从而获得其中unwrap表示对相位差进行解缠操作；

步骤3.6、根据步骤3.4的第一次成像回波时间TE1求得的相位数据拟合表征当前射频场区域均匀性的像数据B₁(r)，

步骤3.7、拟合B₀(r)的频谱加权统计直方图His(Ω)；

步骤3.8、采用洛伦兹线形L(Ω)和步骤3.7拟合的频谱加权统计直方图His(Ω)直接作卷积计算，获得虚拟线形F(Ω)；

步骤3.9、计算虚拟线形F(Ω)的相关性能指标Crit。

4.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3.6包括以下步骤：

步骤3.6.1、计算相位数据的绝对值谱数据并获得绝对值谱最大的强度

步骤3.6.2、计算表征当前射频场区域均匀性的像数据

5.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3.7的频谱加权统计直方图His(Ω)是基于以下公式：

His(Ω)＝∫_rδ(Ω-γ·B₀(r))·w(r)dr

式中，Ω为谱图频率坐标，δ为狄拉克函数，γ为采样核的旋磁比，离散分布各点的权重影响w(r)由以下公式获得：

w(r)＝sin(B₁(r)·α)·B₁(r)

式中，α为射频脉冲产生的翻转角，B₁(r)为表征当前射频场区域均匀性的像数据。

6.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3.8中洛伦兹线形L(Ω)由下式得到：

$L\left(\mathrm{\&Omega;}\right)=\frac{\mathrm{\&lambda;}}{\mathrm{\&pi;}}\&CenterDot;\frac{1}{{\mathrm{\&lambda;}}^2+{(\mathrm{\&Omega;}-{\mathrm{\&Omega;}}_0)}^2}$

式中，Ω为谱图频率坐标，理想情况下洛伦兹线形L(Ω)表示为关于Ω₀＝0轴对称，且洛伦兹线形L(Ω)的半高宽λ＝1/(π·T₂)，T₂为样品的自旋-自旋弛豫时间。

7.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3.9中的相关性能指标Crit为虚拟线形F(Ω)的对称性Crit_symm或包络线E(Ω)的半高宽Crit_half或包络误差Crit_env，

其中虚拟线形F(Ω)的对称性Crit_symm由以下步骤获得：

比较虚拟线形F(Ω)最大峰强度点左右两边频率对称位置的强度差求得虚拟线形F(Ω)的对称性Crit_symm：

${\mathrm{Crit}}_{symm}=\frac{\mathrm{\&Sigma;}|F({\mathrm{\&Omega;}}_{\max }+\mathrm{\&Omega;})-F({\mathrm{\&Omega;}}_{max}-\mathrm{\&Omega;})|}{\mathrm{\&Sigma;}F\left(\mathrm{\&Omega;}\right)}$

式中，Ω_max表示为最大峰强度对应的谱图频率坐标；

其中包络线E(Ω)的半高宽Crit_half和包络误差Crit_env由以下步骤获得：

选择典型洛伦兹线形L(Ω)作为包络线E(Ω)，将虚拟线形F(Ω)包络在其内，关于包络线E(Ω)的求解可通过最小化惩罚函数P_func获得：

$P_{func}={\mathrm{\&Sigma;P}}_{func}\left(\mathrm{\&Omega;}\right)=\left\{\begin{array}{cc}{\left(Diff\right(\mathrm{\&Omega;}\left)\right)}^2,& Diff\left(\mathrm{\&Omega;}\right)\&GreaterEqual;0\\ {(200\&CenterDot;|Diff\left(\mathrm{\&Omega;}\right)\left|\right)}^2,& Diff\left(\mathrm{\&Omega;}\right)<0\end{array}\right.$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线E(Ω)和虚拟线形F(Ω)的差；

采用包络法评价磁场均匀性直接选择计算包络线E(Ω)的半高宽Crit_half；

计算包络线E(Ω)和虚拟线形F(Ω)的包络误差Crit_env：

${\mathrm{Crit}}_{env}=\frac{\mathrm{\&Sigma;}\left|Diff\left(\mathrm{\&Omega;}\right)\right|}{\mathrm{\&Sigma;}F\left(\mathrm{\&Omega;}\right)}$

式中Diff(Ω)＝E(Ω)-F(Ω)表征包络线E(Ω)和虚拟线形F(Ω)的差。

8.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤4中的多维下降单纯形法的迭代搜索包括以下步骤：

步骤4.1、利用当前磁场状态下N个匀场线圈的电流值和步骤2求得的各匀场线圈的搜索步径d_j建立(N+1)×N的空间矩阵M(x_ij)，具体如下所示：

其中，M(x_ij)中单一矩阵元素x_ij代表各匀场线圈具有的电流大小；矩阵的列编号依次对应于所选的N个匀场线圈，i∈{0～N}，j∈{1～N}；矩阵中每个行向量(x_i1,x_i2,...,x_ij,...,x_iN)，均构成一个空间点，则初始化建立的N+1个空间点表示通有不同匀场电流后的N+1个磁场状态：行编号0表示无匀场线圈电流变化的初始磁场状态，行编号1～N分别对应于第1～N个线圈增加电流增量大小为d_j后的磁场状态；

步骤4.2、利用虚拟线形F(Ω)的相关性能指标Crit越小则磁场越均匀的性质，将相关性能指标Crit作为搜索得到磁场均匀性好坏的评价标准Q，Q＝1/Crit，则对于具有不同匀场电流的空间各点(x_i1,x_i2,...,x_ij,...,x_iN),i∈{0～N}，j∈{1～N}，均有相应的评价值Q(i)；

步骤4.3、根据步骤4.1初始化N+1个空间点并拟合出各点的评价值Q(i)；

步骤4.4、依据评价标准Q的大小排序N+1个点，其中数值最小的点为W，数值最大的点为B，中间点设为H；当N+1为奇数时有1个H点，当N+1为偶数时有2个H点，分别表示为H₁和H₂两点，两者大小比较为评价值Q(H₁)≥评价值Q(H₂)；

步骤4.5、拟合除W外的所有N个空间点的中心，设为M，并产生由W指向M的向量；

步骤4.6、沿W至M向量方向拟合W的反射点R，R视为W关于M呈中心对称的点，并计算其评价值Q(R)；

步骤4.7、比较评价值Q(R)和评价值Q(B)的大小，若评价值Q(R)>评价值Q(B)，则继续沿方向扩展至反射扩展点RE，RE视为方向上M关于R的反射点，并计算评价值Q(RE)，进入步骤4.8；若评价值Q(R)≤评价值Q(B)，则跳至步骤4.9；

步骤4.8、比较评价值Q(RE)和评价值Q(R)的大小，两者较大的点取代W点，并跳至步骤4.13；

步骤4.9、比较评价值Q(R)和评价值Q(H)的大小，根据H点的个数分为以下两种情况：

当N+1为奇数时，若评价值Q(R)>评价值Q(H)，则将R取代W，并跳至步骤3.13，若评价值Q(R)≤评价值Q(H)，则跳至步骤3.10；

当N+1为偶数时，若评价值Q(R)>评价值Q(H₂)时，则将R取代W，并跳至步骤3.13，若评价值Q(R)≤评价值Q(H₂)，则跳至步骤3.10；

步骤4.10、比较评价值Q(R)和评价值Q(W)的大小，若评价值Q(R)>评价值Q(W)，则将R取代W，并跳至步骤4.13；若评价值Q(R)≤评价值Q(W)，取M和W的中点C，即收缩点C，并求解其评价值Q(C)，进入步骤4.11；

步骤4.11、比较评价值Q(C)和评价值Q(W)的大小，若评价值Q(C)>评价值Q(W)，则C取代W，并跳至步骤4.13；若评价值Q(C)≤评价值Q(W)，则迭代无更新点，跳至步骤4.12；

步骤4.12、由于迭代无更新点，选择步骤4.4中数据空间的评价值第二差的点W+1作为W，并返回步骤4.5重新开始；若仍无更新点，则继续选择点W+2、W+3，...直至B作为W返回操作步骤4.5；若采用最大的点B开始迭代仍无更新点，则截取除点B外的所有N个点和B的中点，设为全收缩点FC_j，j∈{1～N}，并使用各全收缩点取代原来的N个点，跳至步骤4.13；

步骤4.13、判断迭代是否收敛，其标准为，更新空间点后的矩阵中，任一列向量的任意两个矩阵元素的差均应小于等于步骤2.6中获得的收敛终止条件r_j，表示如下：

|x_aj-x_bj|≤r_j,a,b∈[0,N],j∈{1～N}

若收敛判断为是，则终止匀场迭代，并返回当前最优评价值Q(B)对应的匀场线圈电流大小；若收敛判断为否，则返回步骤4.4，继续迭代搜索。

9.根据权利要求2所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤2中所有谱图频率坐标Ω的积分范围为不低于最大谱峰高度的15％的频率区域；步骤2中预设匀场线圈的电流增量h范围为5～50mA。

10.根据权利要求3所述的一种基于不均匀磁场拟合线形的自动搜索匀场方法，其特征在于，所述的步骤3.4中相位数据的有效范围为不低于采样数据的绝对值谱最大谱峰高度的25％的采样数据区域。