一种基于射频场空间分布加权的梯度匀场方法与流程

本发明涉及核磁共振波谱仪的梯度匀场技术，具体地说涉及一种基于射频场空间分布加权的梯度匀场方法。

背景技术：

梯度匀场是核磁共振波谱技术中最直接、有效的自动匀场方法。借助梯度场空间编码的方式实现相位差成像，测量并拟合待匀场的静磁场场图和各匀场线圈产生的磁场分量的影响(匀场线圈的场图ωshim(r,j))，从而使得匀场的过程变成从数学上处理一个线性方程：

式中，n表示匀场线圈的个数，np表示磁场空间有效像素点r的个数。这样，所需的匀场线圈加载电流增量x(j)(简称匀场电流)可以通过最小二乘法拟合剩磁场空间的场图ωresidual(r)得到全局最优解。因此，如何获取有效的空间磁场分布来准确表征静磁场场图和匀场线圈场图是计算出符合真实磁场均匀性的匀场电流的关键。

核磁共振仪器对场图的测量依赖于探头射频线圈的设计结构，其决定了射频场空间分布(b1(r))，进而影响激发和被检测信号的质量。现代核磁共振波谱仪通常以增加射频线圈的长度来满足对于高检测灵敏度的性能诉求，由于射频场激发在两端衰减较快，容易形成自旋密度较低的边缘区域，使得较弱的有用信号埋没在噪声中，最终导致相位数据畸变失真和拟合计算的匀场电流失效。

为了考虑射频场空间分布对场图测量的影响，现有仪器假定梯度场对激发区域的(归一化)投影轮廓bprofile(r)为射频场空间分布，通过加权bprofile(r)的最小二乘法计算匀场电流。该方法考虑了射频场两端信号衰减过快对相位数据的影响，弱化边缘数据的拟合权重，一定程度上改善磁场测量造成的匀场电流计算误差。但是，这种方式仍然存在以下问题：(1)梯度场对信号有展宽效应导致投影轮廓bprofile(r)是无法真实表征射频场空间分布b1(r)，weiger【markusweigeretal..gradientshimmingwithspectrumoptimisation.《journalofmagneticresonance》.2006，第38-48页】提出一种通过轮廓的指数函数[bprofile(r)]^k近似b1(r)的方法，但是指数项k依赖于用户样品溶剂的选择和特定先验函数的设定(受检测信号灵敏度的影响)；(2)梯度场的非线性有可能导致轮廓bprofile(r)的畸变和位移，产生梯度场非线性的主要因素有梯度线圈设计本身的一阶线性不纯度以及真实不均匀静磁场对梯度场空间编码的干扰，对于核磁仪器来说，两者的影响都是无法消除的。

技术实现要素：

本发明提出一种基于射频场空间分布加权的梯度匀场方法，通过直接测量探头的实际射频场空间分布，加权于梯度匀场的静磁场场图和匀场线圈场图拟合计算匀场电流，通过剩磁场空间分布的谱学特性建立虚拟谱图并评价线形指标，最终优化迭代直至收敛。

本发明的技术方案是这样实现的：

一种基于射频场空间分布加权的梯度匀场方法，包括以下步骤：

s1，测量探头射频线圈中心的射频场空间分布；

s2，根据选择的匀场线圈组合设定一维梯度匀场的脉冲序列并调制序列参数；

s3，测量并拟合匀场线圈的场图；

s4，将射频场空间分布针对匀场线圈的场图进行加权；

s5，测量并拟合待匀场的静磁场场图；

s6，将射频场空间分布针对待匀场的静磁场场图进行加权；

s7，计算匀场线圈的电流改变量；

s8，通过剩磁场空间分布谱学特性判断是否需要迭代。

优选地，步骤1具体包括：

步骤1.1、通过磁场空间测量仪采集出探头射频线圈中心沿轴线方向的射频场强度rf(k)，其中采集总点数为k；

步骤1.2、计算一维射频场空间分布

优选地，步骤2具体包括：

步骤2.1、根据所选的n个匀场线圈的组合设定梯度回波的脉冲序列，假定匀场组合含有z轴方向的匀场线圈时，采用一维梯度回波脉冲序列并需要两次成像采样；

步骤2.2、读入脉冲序列默认参数：一维梯度回波脉冲序列的参数。

优选地，步骤3具体包括：

步骤3.1、初始化n个匀场线圈电流值value(j),j＝1,2,3...n，将初始化通有匀场线圈电流的磁场状态设为基础磁场状态0：value(1),value(2),...,value(n)，并进行梯度回波采样，获得成像回波数据：当采用一维梯度回波脉冲序列时，需要进行两次回波采集；

步骤3.2、设置所选的各个匀场线圈在采样过程中的电流变化量δchange(j),j＝1,2,3...n：其中，n为匀场线圈的个数；

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

与步骤3.1相同方式进行采样，并获得成像回波数据；

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

相位：(φ1(r,i),te1),(φ2(r,i),te2),r＝1,2,...,np；i＝0,1,2,...,n

其中，相位数据φ1和φ2分别和回波采样的时间te1和te2相关联；

i表征不同的磁场状态，当i＝0时为基础磁场状态0采样得到的信号强度和相位数据，当i＝1,2,3…n时表示依次改变n个匀场线圈电流值后的磁场状态1～n采样得到的信号强度和相位数据；

r表示有效像素点，r的个数np表征信号相位在z方向上有效点的个数；

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

phl(r,j)＝(φl(r,j)-φl(r,0)),r＝1,2,...np；j＝1,2,3...n；l＝1,2

其中j＝1,2,3…n表示对应的n个匀场线圈；l＝1,2表示单个匀场线圈所对应的两次扫描的成像回波时间te1和te2；r表示有效像素点。

步骤3.6、将步骤3.5中获得的表征各个匀场线圈单独作用的相位中的第二次成像回波时间te2和第一次成像回波时间te1对应得到的相位数据ph2和ph1作差，得到成像的相位差并进行相位解缠：

δφ21(r,j)＝unwrap(ph2(r,j)-ph1(r,j)),r＝1,2,...np；j＝1,2,3...n

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

步骤3.7、初始化各个匀场线圈场图ωshim(r,j)：

ωshim(r,j)＝δφ21(r,j)/(te2-te1),r＝1,2,...np；j＝1,2,3...n。

优选地，步骤4具体包括：

步骤4.1、对一维射频场空间分布b1(k)进行线性插值，得到插值后的一维射频场空间分布

其中，有效像素点r的个数np表征信号相位在z方向上有效点的个数，interp表示对测量的一维射频场空间分布b1(k)进行线性插值操作，满足插值后的点数等于np；

步骤4.2、拟合射频场空间分布的权重函数weig(r)：

其中r表示有效像素点，max表示获取插值后的一维射频场空间分布的最大强度值；

步骤4.3、对各个匀场线圈场图ωshim(r,j)进行加权：

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

优选地，步骤5具体包括：

步骤5.1、记录待匀场的n个匀场线圈的电流值currvalue(j),j＝1,2,3...n，将此时通有匀场线圈电流的磁场状态设为待匀场的磁场状态，进行采样并获得采样数据；

步骤5.2、对步骤5.1的采样数据进行傅里叶变换，获得一系列表征待匀场的静磁场的相位：

(phb01(r),te1),(phb02(r),te2),r＝1,2,...,np

其中，相位数据phb01，phb02分别表示和回波采样的时间te1，te2相关联；r表示有效像素点，r的个数np表征信号相位在z方向上有效点的个数。

步骤5.3、初始化待匀场的静磁场场图将步骤5.2获得的第二次成像回波时间te2和第一次成像回波时间te1对应得到的相位数据phb02和phb01作差，得到成像的相位差并进行相位解缠，最终获得表征待匀场的静磁场场图

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

优选地，步骤6具体包括：

对待匀场的静磁场场图进行加权：

其中r表示有效像素点，weig(r)表示射频场空间分布的权重函数

优选地，步骤7具体包括：

将加权后的各个匀场线圈场图表示为矩阵a(r,j)，加权后的待匀场的静磁场场图表示为向量b(r)，x(j),j＝1,2,3...n代表匀场线圈的电流改变量，则计算匀场线圈电流改变量简化为求解线性方程组a(r,j)·x(j)＝b(r)；可以采用奇异值分解方法得到匀场线圈的电流改变量x(j),j＝1,2,3...n。

优选地，步骤8中具体包括：

步骤8.1、拟合计算剩磁场空间分布的场图ωresidual(r)：

其中，ωshim(r,j)表示各个匀场线圈场图，表示待匀场的静磁场场图，x(j)表示得到的匀场线圈的电流改变量，j＝1,2,3…n表示对应的n个匀场线圈，r表示有效像素点；

步骤8.2、拟合表征剩磁场空间分布的场图ωresidual(r)的频率分布统计直方图hisg(ω)：

hisg(ω)＝∫rδ(ω-ωresidual(r))·w(r)dr,r＝1,2,...np

其中，ω为谱图频率坐标，δ为狄拉克函数，γ为采样核的旋磁比；

离散分布各点的权重影响w(r)可由以下公式得到：

w(r)＝sin(weig(r)·α)·weig(r),r＝1,2,...np

式中，α为步骤2.2中射频脉冲产生的翻转角，weig(r)表示步骤4.2的射频场空间分布的权重函数；

步骤8.3、采用理想情况的洛伦兹线形l(ω)和步骤8.2拟合的频率分布统计直方图hisg(ω)直接作卷积计算，获得剩磁场空间分布的虚拟线形f(ω)：

其中，洛伦兹线形l(ω)由以下公式得到：

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

步骤8.4、采用频谱包络法将剩磁场空间分布的虚拟线形f(ω)包络在一个典型的包络线函数e(ω)＝h1/[(h2·h2+(ω-h3)²]中，通过最小化惩罚函数pfunc求解包络线函数的系数项h1，h2和h3，并得到包络线函数e(ω)：

式中diff(ω)＝e(ω)-f(ω)表征包络线函数e(ω)和剩磁场空间分布的虚拟线形f(ω)的差。

步骤8.5、将包络线函数e(ω)的半峰宽fwhm作为迭代的评价指标；

步骤8.6、判断迭代是否收敛，即包络线函数e(ω)的半峰宽fwhm满足小于等于收敛终止条件crithalf：若收敛判断为是，则将计算的匀场线圈电流改变量x(j),j＝1,2,3...n分别加至各线圈原电流值上currvalue(j),j＝1,2,3...n并终止匀场迭代；若收敛判断为否，则将计算的匀场线圈电流改变量x(j),j＝1,2,3...n分别加至各线圈原电流值上currvalue(j),j＝1,2,3...n，并返回步骤s5继续迭代。

本发明产生的有益效果为：

(1)本发明有效考虑真实射频场对梯度匀场的影响，通过实际射频线圈产生的场分布降低边缘弱激发带来的相位数据畸变的权重，提高匀场的效率和准确性。特别是针对静磁场均匀性较差的情况，高阶匀场线圈的磁场敏感性不强，边缘相位畸变数据的抑制有助于降低高阶匀场电流的越界可能。

(2)优化梯度匀场的评价指标，常规梯度匀场的评价指标采用迭代最小化剩磁场空间分布(归零)的方式，探头射频场衰减的不一致导致场数据点的计算误差、降低匀场准确性。本方法采用剩磁场空间分布的谱学特性建立虚拟谱图并评价线形指标，可以有效建立起静磁场均匀性和诸如谱的半峰宽等评价指标的对应关系，提高匀场迭代的准确性和效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图。

图2为通过磁场空间测量仪(rfb1fieldmeasurement)采集出探头射频线圈中心沿轴线方向的射频场强度rf(k),有效范围为[-20,20]mm，图中所示的测量的三条射频场强度rf(k)曲线分别取自长度为(rf_coillength)为15mm、18mm和21mm的射频线圈。

图3表示随着射频线圈长度(rf_coillength)的增加，梯度线圈所产生的非线性梯度场(gradshape)导致轮廓bprofile(r)的畸变和位移。图3(a)表示15mm射频线圈长度的轮廓bprofile(r)，由于射频场rf(k)激发在梯度场的线性区间，因此轮廓bprofile(r)无畸变；图3(b)表示18mm射频线圈长度的轮廓bprofile(r)，边缘处受到梯度非线性影响；图3(c)表示21mm射频线圈长度的轮廓bprofile(r)，畸变较为严重。

图4表示加权前匀场线圈场图ωshim(r,j)和加权后匀场线圈场图的比较，射频场空间分布的权重函数weig(r)弱化了场图的边缘强度。(a)-(f)分别表示匀场线圈场图z1-z6。

图5射频场空间分布加权对(a)梯度匀场计算剩磁场空间分布ωresidual(r)的影响及(b)根据剩磁场空间分布拟合的虚拟线形f(ω)。

具体实施方式

下面将结合本发明实施例对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例使用中科波谱quantum-ⅰplus400mhz核磁共振谱仪，其中配备的匀场线圈为中科牛津23-shims、探头为中科牛津h,f/x,f-5。

一种基于射频场空间分布加权的梯度匀场方法，包括以下步骤，

s1，测量探头射频线圈中心的射频场空间分布：具体包括以下步骤，

步骤1.1、通过磁场空间测量仪采集出探头射频线圈中心沿轴线方向的射频场强度rf(k)，有效范围设定为-12.0mm～12.0mm，数据点采集步径设置为0.2mm，采集总点数为121，则k＝1,2,...121。

步骤1.2、计算一维射频场空间分布

s2，选择需要进行梯度匀场的n(本实例中n＝6)个匀场线圈：根据所选的匀场线圈组合自动设定一维梯度回波的脉冲序列；读入并调整默认的脉冲序列参数，具体如下：

步骤2.1、根据所选的室温匀场线圈组合设定梯度回波的脉冲序列，假定匀场组合含有z轴方向(轴向或on-axis)的匀场线圈(如z1～z6)时，采用一维梯度回波(pfgste)脉冲序列并需要两次成像采样。

步骤2.2、读入脉冲序列默认参数：一维梯度回波(pfgste)脉冲序列的参数如下表1所示；

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

s3，测量并拟合匀场线圈的场图，具体包括以下步骤：

步骤3.1、初始化n个匀场线圈电流值value(j),j＝1,2,3...n，将初始化通有匀场线圈电流的磁场状态设为基础磁场状态(磁场状态0：value(1),value(2),...,value(n))，并进行梯度回波采样，获得成像回波数据：当采用一维梯度回波脉冲序列时，需要进行两次回波采集。

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

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

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

步骤3.4、对步骤3.1中磁场状态0及步骤3.3中磁场状态1～n的所有采样数据进行傅里叶变换，获得一系列表征不同磁场状态的信号相位：

相位：(φ1(r,i),te1),(φ2(r,i),te2),r＝1,2,...,np；i＝0,1,2,...,n

其中，相位数据φ1和φ2分别和回波采样的时间te1和te2相关联；i表征不同的磁场状态，当i＝0时为基础磁场(磁场状态0)采样得到的信号强度和相位数据，当i＝1,2,3…n时表示依次改变n个匀场线圈电流值后(磁场状态1～n)采样得到的信号强度和相位数据；r表示有效像素点，r的个数np表征信号相位在z方向上有效点的个数(即z方向不低于最大信号强度25％的区域范围的点个数)。

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

phl(r,j)＝(φl(r,j)-φl(r,0)),r＝1,2,...np；j＝1,2,3...n；l＝1,2

其中j＝1,2,3…n表示对应的n个匀场线圈；l＝1,2表示单个匀场线圈所对应的两次扫描的成像回波时间te1和te2；r表示有效像素点。

步骤3.6、将步骤3.5中获得的表征各个匀场线圈单独作用的相位中的第二次成像回波时间te2和第一次成像回波时间te1对应得到的相位数据ph2和ph1作差，得到成像的相位差并进行相位解缠：

δφ21(r,j)＝unwrap(ph2(r,j)-ph1(r,j)),r＝1,2,...np；j＝1,2,3...n

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

步骤3.7、初始化各个匀场线圈场图(频率数据)ωshim(r,j)：

ωshim(r,j)＝δφ21(r,j)/(te2-te1),r＝1,2,...np；j＝1,2,3...n

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

s4、将射频场空间分布针对匀场线圈的场图进行拟合加权：具体包括以下步骤，

步骤4.1、对步骤1.2中一维射频场空间分布b1(k)进行线性插值，得到插值后的一维射频场空间分布

其中，interp表示对测量的一维射频场空间分布b1(k)进行线性插值操作，满足插值后的点数等于步骤3.4所述有效点数np。

步骤4.2、拟合射频场空间分布的权重函数weig(r)：

其中r表示有效像素点，max表示获取插值后的一维射频场空间分布的最大强度值。

步骤4.3、对步骤3.7中各个匀场线圈场图ωshim(r,j)进行加权：

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

s5，测量并拟合待匀场的静磁场场图：具体包括以下步骤，

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

步骤5.2、对步骤5.1的采样数据进行傅里叶变换，获得一系列表征待匀场的静磁场的相位：

(phb01(r),te1),(phb02(r),te2),r＝1,2,...,np

其中，相位数据phb01，phb02分别表示和回波采样的时间te1，te2相关联；r表示有效像素点，r的个数np表征信号相位在z方向上有效点的个数，等同于步骤3.4。

步骤5.3、初始化待匀场的静磁场场图(频率数据)将步骤5.2获得的第二次成像回波时间te2和第一次成像回波时间te1对应得到的相位数据phb02和phb01作差，得到成像的相位差并进行相位解缠，最终获得表征待匀场的静磁场场图

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

s6，将射频场空间分布针对待匀场的静磁场场图进行拟合加权：对步骤5.3中待匀场的静磁场场图进行加权：

其中r表示有效像素点，weig(r)表示步骤4.2的射频场空间分布的权重函数。

s7，最小二乘法拟合计算匀场线圈的电流改变量：将步骤4.3中加权后的各个匀场线圈场图表示为矩阵a(r,j)，步骤6中加权后的待匀场的静磁场场图表示为向量b(r)，x(j),j＝1,2,3...n代表匀场线圈的电流改变量，则计算匀场线圈电流改变量简化为求解线性方程组a(r,j)·x(j)＝b(r)。可以采用奇异值分解(singularvaluedecomposition，svd)方法得到匀场线圈的电流改变量x(j),j＝1,2,3...n。

s8，通过剩磁场空间分布的谱学特性判断是否需要迭代：具体包括以下步骤，

步骤8.1、拟合计算剩磁场空间分布的场图ωresidual(r)：

其中，ωshim(r,j)表示步骤3.7中各个匀场线圈场图，表示步骤5.3的待匀场的静磁场场图，x(j)表示步骤7得到的匀场线圈的电流改变量，j＝1,2,3…n表示对应的n个匀场线圈，r表示有效像素点。

步骤8.2、拟合表征剩磁场空间分布的场图ωresidual(r)的频率分布统计直方图hisg(ω)：

hisg(ω)＝∫rδ(ω-ωresidual(r))·w(r)dr,r＝1,2,...np

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

w(r)＝sin(weig(r)·α)·weig(r),r＝1,2,...np

式中，α为步骤2.2中射频脉冲产生的翻转角，weig(r)表示步骤4.2的射频场空间分布的权重函数。

步骤8.3、采用理想情况的洛伦兹线形l(ω)和步骤8.2拟合的频率分布统计直方图hisg(ω)直接作卷积计算，获得剩磁场空间分布的虚拟线形f(ω)：

其中，洛伦兹线形l(ω)由以下公式得到：

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

步骤8.4、采用频谱包络法将剩磁场空间分布的虚拟线形f(ω)包络在一个典型的包络线函数e(ω)＝h1/[(h2·h2+(ω-h3)²]中，通过最小化惩罚函数pfunc求解包络线函数的系数项h1，h2和h3，并得到包络线函数e(ω)：

式中diff(ω)＝e(ω)-f(ω)表征包络线函数e(ω)和剩磁场空间分布的虚拟线形f(ω)的差。

步骤8.5、将包络线函数e(ω)的半峰宽fwhm作为迭代的评价指标。

步骤8.6、判断迭代是否收敛，即包络线函数e(ω)的半峰宽fwhm满足小于等于收敛终止条件crithalf(fwhm≤crithalf，本实例中crithalf＝0.1hz)：

若收敛判断为是，则将计算的匀场线圈电流改变量x(j),j＝1,2,3...n分别加至各线圈原电流值上currvalue(j),j＝1,2,3...n并终止匀场迭代；

若收敛判断为否，则将计算的匀场线圈电流改变量x(j),j＝1,2,3...n分别加至各线圈原电流值上currvalue(j),j＝1,2,3...n，并返回步骤5继续迭代。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。