CN110780246B - 一种基于射频场空间分布加权的梯度匀场方法 - Google Patents
一种基于射频场空间分布加权的梯度匀场方法 Download PDFInfo
- Publication number
- CN110780246B CN110780246B CN201911059277.0A CN201911059277A CN110780246B CN 110780246 B CN110780246 B CN 110780246B CN 201911059277 A CN201911059277 A CN 201911059277A CN 110780246 B CN110780246 B CN 110780246B
- Authority
- CN
- China
- Prior art keywords
- magnetic field
- shimming
- radio frequency
- coil
- field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000009826 distribution Methods 0.000 title claims abstract description 96
- 238000000034 method Methods 0.000 title claims abstract description 37
- 230000008859 change Effects 0.000 claims abstract description 47
- 230000003068 static effect Effects 0.000 claims abstract description 40
- 238000001208 nuclear magnetic resonance pulse sequence Methods 0.000 claims abstract description 22
- 239000000523 sample Substances 0.000 claims abstract description 17
- 238000001228 spectrum Methods 0.000 claims abstract description 11
- 238000005070 sampling Methods 0.000 claims description 42
- 238000003384 imaging method Methods 0.000 claims description 28
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000010586 diagram Methods 0.000 claims description 7
- 230000000694 effects Effects 0.000 claims description 7
- 230000009471 action Effects 0.000 claims description 6
- 238000011156 evaluation Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 4
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000012512 characterization method Methods 0.000 claims description 2
- 230000003595 spectral effect Effects 0.000 abstract description 5
- 238000005481 NMR spectroscopy Methods 0.000 description 5
- 238000005259 measurement Methods 0.000 description 5
- 230000005284 excitation Effects 0.000 description 4
- 230000035945 sensitivity Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000002904 solvent Substances 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/28—Details of apparatus provided for in groups G01R33/44 - G01R33/64
- G01R33/38—Systems for generation, homogenisation or stabilisation of the main or gradient magnetic field
- G01R33/387—Compensation of inhomogeneities
- G01R33/3875—Compensation of inhomogeneities using correction coil assemblies, e.g. active shimming
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/28—Details of apparatus provided for in groups G01R33/44 - G01R33/64
- G01R33/32—Excitation or detection systems, e.g. using radio frequency signals
- G01R33/36—Electrical details, e.g. matching or coupling of the coil to the receiver
- G01R33/3607—RF waveform generators, e.g. frequency generators, amplitude-, frequency- or phase modulators or shifters, pulse programmers, digital to analog converters for the RF signal, means for filtering or attenuating of the RF signal
Landscapes
- Physics & Mathematics (AREA)
- Condensed Matter Physics & Semiconductors (AREA)
- General Physics & Mathematics (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明提出了一种基于射频场空间分布加权的梯度匀场方法，包括以下步骤：S1，测量探头射频线圈中心的射频场空间分布；S2，根据选择的匀场线圈组合设定一维梯度匀场的脉冲序列并调制序列参数；S3，测量并拟合匀场线圈的场图；S4，将射频场空间分布针对匀场线圈的场图进行加权；S5，测量并拟合待匀场的静磁场场图；S6，将射频场空间分布针对待匀场的静磁场场图进行加权；S7，计算匀场线圈的电流改变量；S8，通过剩磁场空间分布谱学特性判断是否需要迭代。本发明通过直接测量探头的实际射频场空间分布，加权于梯度匀场的静磁场场图和匀场线圈场图拟合计算匀场电流，通过剩磁场空间分布的谱学特性建立虚拟谱图并评价线形指标，最终优化迭代直至收敛。
Description
技术领域
本发明涉及核磁共振波谱仪的梯度匀场技术，具体地说涉及一种基于射频场空间分布加权的梯度匀场方法。
背景技术
梯度匀场是核磁共振波谱技术中最直接、有效的自动匀场方法。借助梯度场空间编码的方式实现相位差成像，测量并拟合待匀场的静磁场场图
和各匀场线圈产生的磁场分量的影响(匀场线圈的场图ωShim(r,j))，从而使得匀场的过程变成从数学上处理一个线性方程：
式中，N表示匀场线圈的个数，NP表示磁场空间有效像素点r的个数。这样，所需的匀场线圈加载电流增量x(j)(简称匀场电流)可以通过最小二乘法拟合剩磁场空间的场图ωResidual(r)得到全局最优解。因此，如何获取有效的空间磁场分布来准确表征静磁场场图和匀场线圈场图是计算出符合真实磁场均匀性的匀场电流的关键。
核磁共振仪器对场图的测量依赖于探头射频线圈的设计结构，其决定了射频场空间分布(B1(r))，进而影响激发和被检测信号的质量。现代核磁共振波谱仪通常以增加射频线圈的长度来满足对于高检测灵敏度的性能诉求，由于射频场激发在两端衰减较快，容易形成自旋密度较低的边缘区域，使得较弱的有用信号埋没在噪声中，最终导致相位数据畸变失真和拟合计算的匀场电流失效。
为了考虑射频场空间分布对场图测量的影响，现有仪器假定梯度场对激发区域的(归一化)投影轮廓Bprofile(r)为射频场空间分布，通过加权Bprofile(r)的最小二乘法计算匀场电流。该方法考虑了射频场两端信号衰减过快对相位数据的影响，弱化边缘数据的拟合权重，一定程度上改善磁场测量造成的匀场电流计算误差。但是，这种方式仍然存在以下问题：(1)梯度场对信号有展宽效应导致投影轮廓Bprofile(r)是无法真实表征射频场空间分布B1(r)，Weiger【Markus Weiger et al..Gradient shimming with spectrumoptimisation.《Journal of Magnetic Resonance》.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表示获取插值后的一维射频场空间分布B1 int(r)的最大强度值；
步骤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)2]中，通过最小化惩罚函数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为通过磁场空间测量仪(RF B1 field measurement)采集出探头射频线圈中心沿轴线方向的射频场强度RF(k),有效范围为[-20,20]mm，图中所示的测量的三条射频场强度RF(k)曲线分别取自长度为(RF_coil length)为15mm、18mm和21mm的射频线圈。
图3表示随着射频线圈长度(RF_coil length)的增加，梯度线圈所产生的非线性梯度场(Grad shape)导致轮廓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-ⅠPlus 400MHz核磁共振谱仪，其中配备的匀场线圈为中科牛津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)。可以采用奇异值分解(Singular ValueDecomposition，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)2]中，通过最小化惩罚函数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继续迭代。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种基于射频场空间分布加权的梯度匀场方法，其特征在于，包括以下步骤：
S1，测量探头射频线圈中心的射频场空间分布；
S2，根据选择的匀场线圈组合设定一维梯度匀场的脉冲序列并调制序列参数；
S3，测量并拟合匀场线圈的场图；
S4，将射频场空间分布针对匀场线圈的场图进行加权；
S5，测量并拟合待匀场的静磁场场图；
S6，将射频场空间分布针对待匀场的静磁场场图进行加权；
S7，计算匀场线圈的电流改变量；
S8，通过剩磁场空间分布谱学特性判断是否需要迭代。
2.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤1具体包括：
步骤1.1、通过磁场空间测量仪采集出探头射频线圈中心沿轴线方向的射频场强度RF(k)，其中采集总点数为k；
步骤1.2、计算一维射频场空间分布
3.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤2具体包括：
步骤2.1、根据所选的N个匀场线圈的组合设定梯度回波的脉冲序列，假定匀场组合含有Z轴方向的匀场线圈时，采用一维梯度回波脉冲序列并需要两次成像采样；
步骤2.2、读入脉冲序列默认参数：一维梯度回波脉冲序列的参数。
4.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤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。
5.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤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表示有效像素点。
6.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤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表示对相位数据进行解缠操作。
7.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤6具体包括：
对待匀场的静磁场场图
进行加权：
其中r表示有效像素点，Weig(r)表示射频场空间分布的权重函数。
8.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤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。
9.如权利要求1所述的一种基于射频场空间分布加权的梯度匀场方法，其特征在于，步骤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 )2]中，通过最小化惩罚函数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继续迭代。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911059277.0A CN110780246B (zh) | 2019-11-01 | 2019-11-01 | 一种基于射频场空间分布加权的梯度匀场方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911059277.0A CN110780246B (zh) | 2019-11-01 | 2019-11-01 | 一种基于射频场空间分布加权的梯度匀场方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110780246A CN110780246A (zh) | 2020-02-11 |
| CN110780246B true CN110780246B (zh) | 2022-05-24 |
Family
ID=69388270
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911059277.0A Active CN110780246B (zh) | 2019-11-01 | 2019-11-01 | 一种基于射频场空间分布加权的梯度匀场方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110780246B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112244814B (zh) * | 2020-10-22 | 2022-09-16 | 无锡鸣石峻致医疗科技有限公司 | 一种单边磁体磁共振的脂肪定量方法及其系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6049206A (en) * | 1996-08-19 | 2000-04-11 | National Research Council Of Canada | Compensation for inhomogeneity of the field generated by the RF coil in a nuclear magnetic resonance system |
| CN103675733A (zh) * | 2013-11-26 | 2014-03-26 | 中国科学院武汉物理与数学研究所 | 一种基于不均匀磁场拟合线形的自动搜索匀场方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE102004002009B4 (de) * | 2004-01-14 | 2006-07-06 | Siemens Ag | Verfahren zum Betrieb eines Magnetresonanzsystems, Magnetresonanzsystem und Computerprogrammprodukt |
| DE602004012290T2 (de) * | 2004-11-27 | 2009-03-05 | Bruker Biospin Ag | Verfahren zum automatischen Shimmen für die Kernspinresonanzspektroskopie |
| US7605589B2 (en) * | 2006-04-10 | 2009-10-20 | Bruker Biospin Ag | Method for automatic shimming for nuclear magnetic resonance spectroscopy |
| US9568573B2 (en) * | 2014-03-17 | 2017-02-14 | Siemens Healthcare Gmbh | Methods and systems for automated magnetic field shimming |
-
2019
- 2019-11-01 CN CN201911059277.0A patent/CN110780246B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6049206A (en) * | 1996-08-19 | 2000-04-11 | National Research Council Of Canada | Compensation for inhomogeneity of the field generated by the RF coil in a nuclear magnetic resonance system |
| CN103675733A (zh) * | 2013-11-26 | 2014-03-26 | 中国科学院武汉物理与数学研究所 | 一种基于不均匀磁场拟合线形的自动搜索匀场方法 |
Non-Patent Citations (1)
| Title |
|---|
| Gradient shimming with spectrum optimisation;Markus Weiger;《ScienceDirect》;20061231;第38-48页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110780246A (zh) | 2020-02-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dowell et al. | Fast, accurate, and precise mapping of the RF field in vivo using the 180 signal null | |
| Shcherbakova et al. | PLANET: an ellipse fitting approach for simultaneous T1 and T2 mapping using phase‐cycled balanced steady‐state free precession | |
| US7348775B2 (en) | Method for automatic shimming for nuclear magnetic resonance spectroscopy | |
| Klassen et al. | Robust automated shimming technique using arbitrary mapping acquisition parameters (RASTAMAP) | |
| US5218299A (en) | Method for correcting spectral and imaging data and for using such corrected data in magnet shimming | |
| US7605589B2 (en) | Method for automatic shimming for nuclear magnetic resonance spectroscopy | |
| Heidemann et al. | VD‐AUTO‐SMASH imaging | |
| Uecker et al. | Image reconstruction by regularized nonlinear inversion—joint estimation of coil sensitivities and image content | |
| Malik et al. | Tailored excitation in 3D with spiral nonselective (SPINS) RF pulses | |
| US5539316A (en) | Shimming method for NMR magnet having large magnetic field inhomogeneities | |
| Brodsky et al. | Characterizing and correcting gradient errors in non‐cartesian imaging: are gradient errors linear time‐invariant (LTI)? | |
| CN104297709B (zh) | 一种基于正则化磁场分布图像重建的梯度匀场方法 | |
| US6268728B1 (en) | Phase-sensitive method of radio-frequency field mapping for magnetic resonance imaging | |
| CN112881959B (zh) | 一种用于磁共振成像的梯度涡流补偿方法及系统 | |
| CN102612657A (zh) | 针对磁共振成像的图像强度校正 | |
| US20150362576A1 (en) | Metal resistant mr imaging | |
| CN110780246B (zh) | 一种基于射频场空间分布加权的梯度匀场方法 | |
| Weiger et al. | Gradient shimming with spectrum optimisation | |
| Splitthoff et al. | SENSE shimming (SSH): a fast approach for determining B0 field inhomogeneities using sensitivity coding | |
| Renou et al. | Radio-frequency pulse calibration using the MISSTEC sequence | |
| CN111830450B (zh) | 确定和消除磁共振设备中射频脉冲与选层梯度之间的时延的方法 | |
| US20170242085A1 (en) | Method and apparatus for determining a b1 field map in a magnetic resonance scanner | |
| Zelaya et al. | Direct visualisation of B1 inhomogeneity by flip angle dependency | |
| Nöth et al. | B1 mapping using an EPI‐based double angle approach: A practical guide for correcting slice profile and B0 distortion effects | |
| Belsley et al. | Optimal flip angles for in vivo liver 3D T1 mapping and B1+ mapping at 3T |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20230315 Address after: No.128, Guanggu 7th Road, Donghu New Technology Development Zone, Wuhan City, Hubei Province, 430000 Patentee after: WUHAN ZHONGKE NIUJIN WAVE SPECTRUM TECHNOLOGY CO.,LTD. Address before: No.128, Guanggu 7th Road, Donghu Development Zone, Wuhan City, Hubei Province, 430000 Patentee before: WUHAN ZHONGKE KAIWU TECHNOLOGY Co.,Ltd. |
|
| TR01 | Transfer of patent right | ||
| PE01 | Entry into force of the registration of the contract for pledge of patent right |
Denomination of invention: A gradient homogenization method based on weighted RF field spatial distribution Granted publication date: 20220524 Pledgee: Qingshan Branch of Wuhan Rural Commercial Bank Co.,Ltd. Pledgor: WUHAN ZHONGKE NIUJIN WAVE SPECTRUM TECHNOLOGY CO.,LTD. Registration number: Y2024980037072 |
|
| PE01 | Entry into force of the registration of the contract for pledge of patent right |