基于磁共振偶极场空间分布的局部场估计方法

摘要

基于磁共振偶极场空间分布的局部场估计方法，涉及核磁共振定量磁化率成像。基于相位图的缠绕数；基于由相位图的幅值、梯度和缠绕数调制的高斯卷积核标准差σ(r)；基于空间信息将感兴趣区域划分为高磁化率变化区域和均匀磁化率区域；利用相位图的空间信息，即引入了相位图的幅值、梯度和归缠绕数逐个体素地调制高斯卷积核的标准差，进而其球心点及其他点的权重，该方法能够有效地去除强磁化率区域在不同成像方向下的背景场并准确地估算了局部场，同时保留了脑组织的完整性，其效果明显优于现有的方法。

说明书

技术领域

本发明涉及核磁共振定量磁化率成像，尤其是涉及基于磁共振偶极场空间分布的局部场估计方法。

背景技术

磁化率是物质的固有属性，能够反映生物组织在外加磁场作用下的磁化程度，是磁共振成像(Magnetic Resonance Imaging,MRI)的主要信号来源之一。由于生物组织的磁化率分布受到组织的铁含量、钙化、血氧饱和度等因素的影响，因此人体不同组织之间、正常组织与病变组织之间的磁化率不同。人脑中灰质有铁的沉积，表现为顺磁性，磁化率值为正；血管中含氧血红蛋白是抗磁性的，脱氧血红蛋白是顺磁性的；脑白质的主要成分是髓磷脂，表现为抗磁性，磁化率值为负。

MRI中常用磁敏感加权成像(Susceptibility Weighted Imaging，SWI)和R2*图像检测组织的磁化率，但这些方法的敏感度和特异性都较低，且SWI只是定性测定磁化率变化，R2*图像则易受强梯度边缘处的体素内自旋相位离散影响。定量磁化率成像(Quantitative Susceptibility Mapping,QSM)利用以往常被忽略的相位信息，经过相位解缠绕和去除背景场的预处理，得到局部磁场的变化特性，通过复杂的场到源反演计算获得组织的定量磁化率图^[1,2]。QSM技术可以有效克服SWI和R2*的上述缺陷，已经广泛应用于脑出血、多发性硬化症及多种退行性神经疾病的研究中，对这些疾病的临床诊断具有重要意义^[3-6]，与传统MRI成像技术形成互补，表现出很高的临床价值及应用前景。

QSM技术是现今医学影像领域的一个研究热点，其面临的主要问题之一是如何去除感兴趣区域(Volume of Interest，VOI)外的背景场对VOI内部场的影响，从而得到准确的局部场估计。局部场估计的准确与否直接影响着QSM的精度。根据背景场的物理原理，目前主要有两种重要的局部场估计模型：其一，为偶极场投影法(Projection onto DipoleField,PDF)^[7]，该方法假设由局部场与背景场在VOI内部是相互正交的，二者的内积为零，通过方程最优化解法即可得到局部场；其二，为相位数据复杂谐波伪影抑制法(Sophisticated>[8]，即假设背景场在VOI内满足调和函数性质，可用归一化的球均值卷积核求解得到局部场。由于采用球均值卷积核，SHARP方法在VOI边缘区域的体素因与VOI外的无效信号做卷积运算而导致与卷积核直径相等的边缘区域模糊，因此这些模糊在QSM重建中被去除。由于SHARP存在腐蚀脑组织边缘的缺陷，研究人员提出了解决拉普拉斯方程边界问题的LBV方法^[9]、变卷积核半径的V-SHARP法^[10]等相应的改进措施，均取得了一定的效果。

由于脑组织与鼻窦或颅骨的交界处的强磁化率值非常高，而且偶极场的空间分布与组织的几何结构、成像方向均有关，随着人脑成像角度的变化，该区域产生的强偶极场对周边的脑组织产生不同方向的强干扰，对局部场的正确估计产生严重影响^[11]。若采用上述的方法，均会在这些区域的局部场图上留下较为严重的相位残余，在后续的QSM重建中留下严重的磁化率伪影，掩盖正常组织结构细节。为消除此类的影响，上述的各种方法均需用掩模板将该区域的脑组织去除以减少相位残余，从而降低QSM伪影。然而，采用不同方向掩膜版会导致脑组织缺失严重，将使我们无法获取完整的大脑诊断信息，这不利于QSM技术的临床应用。因此，在去除强磁化率区域在多方向成像下的背景场，准确估算局部场时，上述方法均存在不足，或脑组织缺失严重，或相位残余严重，皆有待进一步提高。

参考文献：

[1]Shmueli,K.,de Zwart,J.A.,van Gelderen,P.,Li,T.Q.,Dodd,S.J.,Duyn,J.H.,2009.Magnetic susceptibility mapping of brain tissue in vivo using MRIphase data.Magn.Reson.Med.62,1510-1522.

[2]Haacke,E.M.,Liu,S.F.,Buch,S.,Zheng,W.L.,Wu,D.M.,Ye,Y.Q.,2015.Quantitative susceptibility mapping:current status and futuredirections.Magn.Reson.Imag.33,1-25.

[3]He,N.Ling,H.Ding,B.Huang,J.Zhang,Y.Zhang,Z.Liu,C.Chen,K.Yan,F,2015.Region-Specific Disturbed Iron Distribution in Early IdiopathicParkinson’s Disease Measured by Quantitative Susceptibility Mapping,HumanBrain Mapping,36,4407-4420.

[4]Wang,Y.,Liu,T.,2015.Quantitative susceptibility mapping(QSM):decoding MRI data for a tissue magnetic biomarker.Magn.Reson.Med.73,82-101.

[5]Bergen,J.,Hua,J.,Lim,I.A.L,Jones,C.K.,Margolis,R.L.,Ross,C.A.,P.C.M.van Zijl,Li,X,2016.Quantitative susceptibility mapping suggests alteredbrain iron in premanifest Huntington disease.Am.J.Neuroradiol.37,789-796.

[6]Lotfipour A K,Wharton S,Schwarz S T,et al,2012.High resolutionmagnetic susceptibility mapping of the substantia nigra in Parkinson'sdisease.[J].Journal of Magnetic Resonance Imaging Jmri,35,48–55.

[7]Liu,T.,Khalidov,I.,de Rochefort,L.,Spincemaille,P.,Liu,J.,Tsiouris,A.J.,2011.A novel background field removal method for MRI usingprojection onto diplole fields(PDF).NMR Biomed.24,1129-1136.

[8]Schweser,F.,Deistung,A.,Lehr,B.W.,Reichenbach,J.R.,2011.Quantitative imaging of intrinsic magnetic tissue properties using MRIsignal phase:an approach to in vivo brain iron metabolism？NeuroImage.54,2789-2807.

[9]Zhou,D.,Liu,T.,Spincemaille,P.,Wang,Y.,2014.Background fieldremoval by solving the Laplacian boundary value problem.NMR Biomed.27,312-319.

[10]Wu,B.,Li,W.,Guidon,A.,Liu,C.L.,2012.Whole brain susceptibilitymapping using compressed sensing.Magn.Reson.Med.67,137-147.

[11]Bao,L.Li,X.Cai,C.Chen,Z.,P.C.M.V.Zijl,2016.QuantitativeSusceptibility Mapping Using Structural Feature Based CollaborativeReconstruction(SFCR)in the Human Brain,IEEE Transactions on MedicalImaging.35,2040-2050.

发明内容

本发明的目的在于提供基于磁共振偶极场空间分布的局部场估计方法。

本发明包括如下步骤：

1)基于相位图的缠绕数，其模型为：

模型中，K(r)为相位图缠绕数，γ为旋磁比，TE为回波时间，B^Δ为偶极场的总场值，为解缠绕后的相位图，为缠绕相位图，W[·]为相位缠绕算子；

2)基于由相位图的幅值、梯度和缠绕数调制的高斯卷积核标准差σ(r)，其模型为：

模型中，C为正的调节常数，B^Δ为总场图，为总场图的梯度，ε为大于零的实数，K(r)∈[0,1]为归一化的相位缠绕数，由此建立的与逐个体素空间信息相关的卷积核表示为：

式中(x,y,z)为高斯核内点的空间坐标；

3)基于空间信息将感兴趣区域划分为高磁化率变化区域和均匀磁化率区域，其模型为：

模型中，为σ(r)倒数，记为为区域分割阈值，定义为其中，与σ_max分别表示的均值和最大值；δ为单位冲击函数，B_int为局部场，为临时变量，表示卷积运算；通过将感兴趣区域分为高磁化率变化区域D₁和均匀磁化率区域D₂，且分别与高斯卷积核和球均值卷积作卷积运算，将背景场去除，得到局部场。

在步骤1)中，所述相位图的缠绕数可包含不同场图变化在不同时间下的演化结果，可用于任何回波时间的局部场估计。

在步骤2)中，所述相位图的幅值、梯度和缠绕数组成的空间信息可逐体素调制高斯卷积核的标准差，当空间信息值越大表示了该区域存在强背景场，则得到一个越小的高斯标准差，即相应的高斯核球心的权重高而周边点值快速衰减，可更好地抑制背景场从而得到为准确的局部场估计。

在步骤3)中，所述基于空间信息将感兴趣区域划分为高磁化率变化区域和均匀磁化率变化区域，可以加快算法的运算速度，空间信息考虑了不同成像方向下的场图分量，可用于多方向的局部场估算。

本发明利用相位图的空间信息，即引入了相位图的幅值、梯度和归缠绕数逐个体素地调制高斯卷积核的标准差，进而其球心点及其他点的权重，该方法能够有效地去除强磁化率区域在不同成像方向下的背景场并准确地估算了局部场，同时保留了脑组织的完整性，其效果明显优于现有的方法。

附图说明

图1为不同成像方向的人脑数据图。(a)解缠绕相位图，(b)空间信息值图，(c)局部场图，(d)磁化率图。

图2为仿真数据实验图。(a)方向1的真实磁化率图、全场图和局部场图，(b)方向2的真实磁化率图、全场图和局部场图，(c)方向3的真实磁化率图、全场图和局部场图，(d)方向4的真实磁化率图、全场图和局部场图。

图3为仿真数据局部场估计结果。(a)方向1由PDF、R-SHARP和iRSHARP三种方法计算得到的局部场图，(b)方向2由PDF、R-SHARP和iRSHARP三种方法计算得到的局部场图，(c)方向3由PDF、R-SHARP和iRSHARP三种方法计算得到的局部场图，(d)方向4由PDF、R-SHARP和iRSHARP三种方法计算得到的局部场图。

图4为不同成像方向下的人脑局部场估计结果。(a)4个方向解缠绕后的相位图，(b)由PDF方法的局部场估计结果，(c)由R-SHARP方法局部场估计结果，(d)由iRSHARP方法局部场估计结果。

图5为不同成像方向下的人脑磁化率重建结果。(a)根据PDF估计的局部场重建的磁化率结果，(b)根据R-SHARP估计的局部场重建的磁化率结果，(C)根据iRSHARP估计的局部场重建的磁化率结果。

具体实施方式

以下实施例将结合附图对本发明作进一步的说明。

图1给出不同成像方向的人脑数据图。本发明能够有效地去除背景场，获得准确的局部场，同时保护了大脑组织的完整性。方法具体实施过程如下：

1)首先，进行数值仿真实验。创建一个128ⅹ128ⅹ64的矩阵，在矩阵中嵌入一个大椭球体用于仿真人脑，在大椭球体中再生成5个小的椭球体分别仿真鼻窦、血管、黑质、苍白球和尾状核，相应的磁化率为9.4ppm，0.3ppm，0.16ppm，0.1ppm和0.05ppm，脑区剩余的部位的磁化率设置为0ppm，在VOI外增加另外1个磁化率值为9.4ppm的椭球体用于模拟其它的背景场源，并在仿真数据中加入SNR＝40的高斯标准噪声。将鼻窦分别向左旋转10°、向右旋转20°及绕向后旋转10°，用于模仿不同的成像方向，如图2所示。

2)其次，选取Liu,T.于2011提出的PDF(Projection onto Dipole Field)方法和Fang,J.S.于2017年提出的R-SHARP(Sophisticated Harmonic Artifact Reduction forPhase Data Using Region Adaptive Kernel)方法作为对比，以说明本发明提出的基于磁共振偶极场空间分布的局部场估计算法能够有效地去除磁共振强磁化率区域在多方向成像下的背景场，可获得准确的局部场。在局部场图(图3)上，如箭头所示处，PDF方法在仿真鼻窦周围存在很强的相位残余，这些相位残余将产生严重的磁化率伪影；R-SHARP在第一个方向上，能够很好地将背景场去除，但是在其他的三个方向存在着较小的相位残余，如图3放大图所示；iRSHARP方法则在四个方向上都无相位残余，且运行时间远少于R-SHARP方法。

为验证本方法在局部场估算，采用均方根误差(RMSE)来定量评价三种方法的处理结果。其中，均方根误差的定义为：

其中，B_int表示计算得到的局部场图，B₀表示真实的局部场，n为计算点数。

通过表1可知，iRSHARP去局部场估计好于其他两种方法，且运行时间较R-SHARP短。

表1：仿真图实验的RMSE和运行时间对比

方法方向1方向2方向3方向4时间(s)PDF0.2250.2010.2020.22337.8R-SHARP0.1450.1450.1450.146122.1iRSHARP0.1440.1440.1440.14458.9

3)最后，进行人脑实验数据，进一步验证本发明方法在实际应用中的可行性。实验数据在飞利浦7T，32通道的人体成像仪上采集得到，采用三维梯度回波序列，成像参数为TR/TE1/ΔTE＝45ms，2ms，2ms，8回波，成像视野为220mm×220mm×110mm，层厚为1mm，数据矩阵为224×224×110。图4列出了不同成像方向下，PDF,R-SHARP和iRSHARP三种方法的局部场估计结果，PDF在鼻窦上方的脑组织留下明显的方向性相位残余，且在脑组织边缘区域也存在很强的相位残余；R-SHARP方法在方向1中，能够很好地去除背景场，但在其他三个方向中，均存在较弱的相位残余；iRSHARP方法在四个方向中，均能很好地去除背景场。将得到的局部场进行磁化率反演，如图5所示，在鼻窦周围的脑组织(如箭头所示)，PDF方法留下的残余相位产生了严重的磁化率伪影并会覆盖其他区域；R-SHARP方法在方向2、3、4余留的弱相位也产生肉眼可观测得到的磁化率伪影；而iRSHARP方法则在各个方向上很好地去除了相位残余信息，在重建的磁化率图上能够分辨脑组织结构细节。

磁共振定量磁化率成像中，大脑组织的偶极场分布与组织的几何结构、组织与主磁场的夹角有关，因而鼻窦、组织与颅骨交界等强磁化率区域会随着头部与主磁场的成像夹角的变化而对周边脑组织产生很强的不同方向的背景场干扰，从而影响局部场的准确估计。基于磁共振偶极场空间分布的局部场估计算法，利用相位图的梯度值、相位幅值和归一化的相位缠绕数组合成的空间信息检测得到大脑组织中的强磁化率区域。空间信息包含着偶极场在不同时间和不同方向下的演化结果，因此通过这些空间信息调制高斯卷积核的中心权重，可以有效地消除强磁化率区域在不同方向下的背景场并准确保留由组织产生的局部场信息，从而有效地抑制磁化率伪影。因此，在多方向的人脑数据实验中，本方法无需去除鼻窦周边的强磁化率组织并能够保证人脑组织的完整性，效果明显优于其他现有方法，具有潜在的临床应用价值。