磁化率是一种描述物质内在磁介质属性的物理量,反映的是物质磁化强度与外界所施加的磁场之间的关系。磁量图(QSM)是一种可在体内对组织宏观磁化率定量的磁共振图像(MRI)技术。QSM可清楚显示具有顺磁性或者反磁性组织的解剖结构。QSM算法是解决从测量到的组织磁场中通过非常规的单位偶极子反卷积获得磁化率分布这一病态问题的方法。研究者们提出了多种基于多角度扫描和面向临床单角度扫描的QSM算法来解决这个病态问题。本文通过介绍QSM的基本概念、系统分析当前多种QSM算法、简要介绍QSM的研究进展,旨在促进QSM技术发展和临床应用。
引用本文： 王帅, 段昶, 张萍, 王春梅, 刘想, 李鸿升, 程建. 磁共振图像中磁量图技术研究进展. 生物医学工程学杂志, 2015, 32(5): 1131-1134. doi: 10.7507/1001-5515.20150200 复制
引言
磁共振成像(magnetic resonance imaging,MRI)是一种无辐射、无侵入性、具有良好软组织对比度的检测方法，二十世纪中后期以来，其临床应用价值日益凸显。常用的MRI成像序列有T1、T2、T2*等加权图像。当磁性物体被置于主磁场中，可引起主磁场幅度和相位的改变，从而导致MRI信号强度和相位的变化。常用的T1、T2、T2*加权图像均来源于幅度信息，而使得相位信息长期以来被研究人员忽视。近年来研究表明相位数据中含有丰富的信息，通过使用数据后处理算法，处理结果可用于临床磁性物质的检测与定量，具有广阔的应用前景。
磁化率(magnetic susceptibility)作为一种生物组织固有的物理属性，可为疾病诊断提供新的组织对比机制。不同组织的磁性不同，如铁、钆是顺磁性物质，钙是反磁性物质[1-2]，引起的MRI相位改变也不同。磁量图(quantitative susceptibility mapping,QSM)就是从磁性物质引起的相位扰动中反解出的磁性物质分布图[3-5]。因此QSM 可用于顺磁性物质和反磁性物质的鉴别。要求得QSM，必须从测量得到的组织磁场中反卷积一个偶极子核(单位磁偶极子场)；而这个偶极子核在k空间(k-space)有一个约54.7°的三维零面，所以求取QSM是一个具有挑战性的病态逆问题[3]。学者们提出了很多解决病态反卷积的方法，这些方法主要分为两类：一类是基于多角度多次扫描数据进行重建；另一类是采用单角度扫描数据进行重建。本文主要通过分析、讨论当前QSM需解决的科学问题、关键步骤和预处理算法等内容，从反卷积算法的角度探讨了各类QSM方法和目前的研究热点，并对今后的发展方向予以展望。
1 磁量图的模型
组织局部磁场δb(r) (简写为δb) 相对于主磁场B0的分布在图像域(简称为r空间，r-space)可以用偶极子核(Dipole Kernel)d(r) =(1/4π)(3cos2θ-1)/|r|3(在r-space简写为d)和组织磁化率分布χ(r)(简写为χ)的卷积来建模[3]:
| $ {\boldsymbol{\delta} _b} = \boldsymbol{d} \otimes \chi $ |
其中δb从MRI扫描的相位图中获得，式(1)中χ即所求的QSM。在k-space中组织的磁化率模型可以用矩阵表示为:
| $ {\Delta _b} = \boldsymbol{D} \cdot \boldsymbol{X} $ |
这里,F为矩阵算子，Fδb代表δb 的离散傅里叶变换，·代表矩阵点乘，X=Fχ；D=Fd为k-space的偶极子核，对应位置的元素为1/3-kz2/k2，k=|k|，k为代表k-space坐标的向量，kz代表k向量在B0方向的分量。
QSM的求解是基于式(1)或者式(2)。式(1)可以看出要从测量值δb得到χ是反卷积问题[3, 6]；对应式(2)是除法：
| $ X = {\boldsymbol{D}^{ - 1}} \cdot {\Delta _b} $ |
D有一个等于零的面，这个面称为圆锥面(cone surface)，如图 1 所示，所以不能通过除法直接求得D-1。这个圆锥面也使得无穷多解可以满足式(1)和(2)，所以求解QSM的分布是一个病态过程。QSM算法的核心正是解决这个病态问题。
图1
40×40×40体素的圆锥面示意图
Figure1.
Cone surface diagram shown in 40×40×40 voxels volume
2 磁量图的关键步骤
从MRI数据中得到QSM的主要步骤包括：解卷褶、背景场移除和反卷积。一个典型的QSM算法-形态学约束偶极子反卷积(morphology enabled dipole inversion，MEDI)[4]的步骤如下：①采集得到MRI的幅度图和相位图；②从幅度图当中得到解剖图P(通常是所有幅度图的平方和开方)和先验加权项M(这个算法中M是不可忽略的梯度项)；③从采集数据(包括相位图和幅度图)当中得到数据加权项W；④从采集数据中得到总场f，并且进行背景场移除得到局部场fL；⑤反卷积QSM算法得到QSM解[2, 4]。
QSM结果的质量不仅依靠反卷积算法，也依靠局部磁场信息的质量。QSM对通过梯度回波(gradient echo，GRE)扫描得到的相位数据的第一步处理是解卷摺2π相位混淆。其他研究领域中的二维或者三维相位解卷摺方法都可以用在QSM的相位解卷摺中。但从QSM原理可知QSM是基于三维的模型，所以三维相位解卷摺方法更适用于QSM的相位预处理。
为了移除感兴趣区域(volume of interest，VOI)磁源之外的背景场干扰，研究人员提出了多类背景场移除方法：第一类是将总磁场分解为一系列基本函数，这些函数中的一部分可认为是背景场，其代表方法是偶极子场投影法(projection onto dipole field，PDF)[7]；第二类是使用谐波函数的均值特性，代表方法有相位数据谐波伪影消除法(sophisticated harmonic artifact reduction for phase data，SHARP)[5]、迭代球均值法(iterative sphere mean Value，iSMV)[8]、规则化约束的相位数据谐波伪影消除法(regularization enabled SHARP，RESHARP)[9]、扩展相位数据谐波伪影消除法(externed-SHARP)[10]；第三类是直接求解拉普拉斯方程，代表方法是拉普拉斯边界法(Laplacian boundary value，LBV)[11]。
经过预处理和一系列反卷积必要的数据准备，通过反卷积就能得到QSM解，即组织磁化率分布。
3 磁量图反卷积算法
3.1 基于多角度多次扫描数据的QSM算法
此类算法采用多角度多次扫描的数据来重建QSM，主要思想是采用过样技术来弥补单次采集偶极子的零区域数据缺失问题。其代表性算法是多角度采样磁化率计算法(calculation of susceptibility through multiple orientation sampling，COSMOS)[3]。
3.2 单角度QSM算法
虽然多角度多次扫描数据可以很好地约束QSM的逆问题，但是实际临床使用中存在难以操作、有配准误差等问题。临床使用中可以接受的扫描方案是采用单次自然体位的扫描。为了开发能临床使用的单角度QSM方法，学者们提出了多种单角度QSM算法，本文也着重对其进行归类和讨论：
3.2.1 基于偶极子核修正的算法
这类算法对偶极子核绝对值大于某个预设值的区域采用直接除法来得到磁化率分布，对绝对值小于这个值(圆锥面附近，cone region，CR)的区域直接赋于不等于零的值，再在k-space做除法。此类方法中具有代表性的是：基于k-space经验值的直接截断法(truncated K-space division,TKD)[12]。TKD改变了偶极子核本身，在理想情况下仍会带来QSM解的偏离；同时此类算法没有考虑噪声的影响，在k-space除法当中会带来噪声放大作用。
3.2.2 基于偶极子核近似的算法
针对TKD在CR区域采用直接替换原来偶极子核的值的武断做法，研究人员尝试用其他办法来对CR区域做处理。加权k-space导数法(weighted k-space derivative，WKD)[13]在做k-space除法之前，在CR区域对磁化率分布和偶极子核同时做求导处理，使得算法在偶极子核的圆锥面上的求解不再包含零域。但由于在圆锥面每一边，D的符号是一致的，所以D的导数将小于D;这将导致该方法更倾向于在CR区域放大噪声,无法得到最佳估计。
3.2.3 基于通用数学先验条件的贝叶斯算法
通用数学先验条件可以用来约束最小化模型，所以多位学者提出使用通用数学条件约束QSM，从而得到唯一解。这些使用在QSM中的先验条件有：①基于梯度L2范数的空间平滑代价函数[6, 14]；②基于梯度的L1范数的稀疏代价函数[6, 15]；③基于全变差范数的稀疏代价函数[16]；④基于小波域的稀疏代价函数[16]；⑤基于两个稀疏函数混合的代价函数[15-16]。该类方法是将数学上的先验条件直接使用在QSM求解当中，并没有使用磁化率分布自身的特性，而且部分算法并没有考虑噪声不是高斯白噪声的条件。
3.2.4 基于物理结构先验条件的贝叶斯算法
通用数学先验条件并没有很好地体现和限定真实磁化率分布自身的特性。如果可以得到某种磁化率分布的特性，那么就可以用这种分布特性来指导图像重建。此类QSM方法的典型例子就是在QSM重建方法中使用特定物理解剖信息的先验知识[4, 14-15, 17]，代表方法是MEDI[4, 15, 17] 方法。Liu等[18]已经证明在高斯白噪声和给定精确限定下，MEDI方法能得到全局最优解。
4 磁量图研究总结与展望
综上所述，QSM作为一种有别于传统的自旋质子密度、T1、T2加权成像的新组织对比方法，扩展了MRI对磁性物质的检测能力，具有广阔的临床应用前景。QSM已经应用在测量脱氧血红蛋白或含铁血黄素沉积；跟踪氧化铁纳米粒子沉积、量化钆等外源性造影剂；检测钙沉积、测量神经退行性疾病中的铁沉积、脑血管病及脑外伤中的微出血灶、脑内静脉血管畸形[2, 19-22] 等科研和临床试验中。研究人员对QSM各个步骤的研究，尤其是伴随着基于单角度扫描的QSM反卷积方法中的从基于偶极子核的探索式方法到基于物理结构贝叶斯先进方法的发展，使得QSM解的精度越来越高，从而使得QSM的应用范围越来越广。
QSM的进一步发展将依靠先进数学和物理方法的使用，同时也需要研究和讨论有关信噪比、场强、分辨率、回波时间、反演的正则化项。未来的QSM研究将包括以下几个方面：①更精确的三维解卷褶算法研究。虽然当前已经存在针对人类脑组织的稳健二维或者三维解卷褶算法，但对于存在陡峭相位梯度的地方解卷摺仍然会失败。②先进的背景场移除方法研究。当前各种背景场移除方法虽然取得了较好的效果，但是并不能完全移除VOI之外的背景场。所以如何更精确地移除VOI之外的背景场是下一步研究的一个方向。③更精确的反卷积方法研究。当前的研究对反卷积算法的规范化、算法与正则化项的关系，以及反卷积结果与信噪比、场强等扫描参数的关系仍然不清楚。④磁化率对比机制研究。当前研究对磁化率对比在体内存在的机制仍然不清楚；当前研究表明非磁化率敏感效应对观测到的相位和局部场也有贡献，但对这个相关度并不清楚。⑤更精确的磁化率MRI成像序列研究。当前的GRE相位信号的扰动不止单纯由磁化率产生，如何分离其他机制产生的影响需要更精确的相位信号分离方法，这需要设计更先进的MRI序列。
引言
磁共振成像(magnetic resonance imaging,MRI)是一种无辐射、无侵入性、具有良好软组织对比度的检测方法，二十世纪中后期以来，其临床应用价值日益凸显。常用的MRI成像序列有T1、T2、T2*等加权图像。当磁性物体被置于主磁场中，可引起主磁场幅度和相位的改变，从而导致MRI信号强度和相位的变化。常用的T1、T2、T2*加权图像均来源于幅度信息，而使得相位信息长期以来被研究人员忽视。近年来研究表明相位数据中含有丰富的信息，通过使用数据后处理算法，处理结果可用于临床磁性物质的检测与定量，具有广阔的应用前景。
磁化率(magnetic susceptibility)作为一种生物组织固有的物理属性，可为疾病诊断提供新的组织对比机制。不同组织的磁性不同，如铁、钆是顺磁性物质，钙是反磁性物质[1-2]，引起的MRI相位改变也不同。磁量图(quantitative susceptibility mapping,QSM)就是从磁性物质引起的相位扰动中反解出的磁性物质分布图[3-5]。因此QSM 可用于顺磁性物质和反磁性物质的鉴别。要求得QSM，必须从测量得到的组织磁场中反卷积一个偶极子核(单位磁偶极子场)；而这个偶极子核在k空间(k-space)有一个约54.7°的三维零面，所以求取QSM是一个具有挑战性的病态逆问题[3]。学者们提出了很多解决病态反卷积的方法，这些方法主要分为两类：一类是基于多角度多次扫描数据进行重建；另一类是采用单角度扫描数据进行重建。本文主要通过分析、讨论当前QSM需解决的科学问题、关键步骤和预处理算法等内容，从反卷积算法的角度探讨了各类QSM方法和目前的研究热点，并对今后的发展方向予以展望。
1 磁量图的模型
组织局部磁场δb(r) (简写为δb) 相对于主磁场B0的分布在图像域(简称为r空间，r-space)可以用偶极子核(Dipole Kernel)d(r) =(1/4π)(3cos2θ-1)/|r|3(在r-space简写为d)和组织磁化率分布χ(r)(简写为χ)的卷积来建模[3]:
| $ {\boldsymbol{\delta} _b} = \boldsymbol{d} \otimes \chi $ |
其中δb从MRI扫描的相位图中获得，式(1)中χ即所求的QSM。在k-space中组织的磁化率模型可以用矩阵表示为:
| $ {\Delta _b} = \boldsymbol{D} \cdot \boldsymbol{X} $ |
这里,F为矩阵算子，Fδb代表δb 的离散傅里叶变换，·代表矩阵点乘，X=Fχ；D=Fd为k-space的偶极子核，对应位置的元素为1/3-kz2/k2，k=|k|，k为代表k-space坐标的向量，kz代表k向量在B0方向的分量。
QSM的求解是基于式(1)或者式(2)。式(1)可以看出要从测量值δb得到χ是反卷积问题[3, 6]；对应式(2)是除法：
| $ X = {\boldsymbol{D}^{ - 1}} \cdot {\Delta _b} $ |
D有一个等于零的面，这个面称为圆锥面(cone surface)，如图 1 所示，所以不能通过除法直接求得D-1。这个圆锥面也使得无穷多解可以满足式(1)和(2)，所以求解QSM的分布是一个病态过程。QSM算法的核心正是解决这个病态问题。
图1
40×40×40体素的圆锥面示意图
Figure1.
Cone surface diagram shown in 40×40×40 voxels volume
2 磁量图的关键步骤
从MRI数据中得到QSM的主要步骤包括：解卷褶、背景场移除和反卷积。一个典型的QSM算法-形态学约束偶极子反卷积(morphology enabled dipole inversion，MEDI)[4]的步骤如下：①采集得到MRI的幅度图和相位图；②从幅度图当中得到解剖图P(通常是所有幅度图的平方和开方)和先验加权项M(这个算法中M是不可忽略的梯度项)；③从采集数据(包括相位图和幅度图)当中得到数据加权项W；④从采集数据中得到总场f，并且进行背景场移除得到局部场fL；⑤反卷积QSM算法得到QSM解[2, 4]。
QSM结果的质量不仅依靠反卷积算法，也依靠局部磁场信息的质量。QSM对通过梯度回波(gradient echo，GRE)扫描得到的相位数据的第一步处理是解卷摺2π相位混淆。其他研究领域中的二维或者三维相位解卷摺方法都可以用在QSM的相位解卷摺中。但从QSM原理可知QSM是基于三维的模型，所以三维相位解卷摺方法更适用于QSM的相位预处理。
为了移除感兴趣区域(volume of interest，VOI)磁源之外的背景场干扰，研究人员提出了多类背景场移除方法：第一类是将总磁场分解为一系列基本函数，这些函数中的一部分可认为是背景场，其代表方法是偶极子场投影法(projection onto dipole field，PDF)[7]；第二类是使用谐波函数的均值特性，代表方法有相位数据谐波伪影消除法(sophisticated harmonic artifact reduction for phase data，SHARP)[5]、迭代球均值法(iterative sphere mean Value，iSMV)[8]、规则化约束的相位数据谐波伪影消除法(regularization enabled SHARP，RESHARP)[9]、扩展相位数据谐波伪影消除法(externed-SHARP)[10]；第三类是直接求解拉普拉斯方程，代表方法是拉普拉斯边界法(Laplacian boundary value，LBV)[11]。
经过预处理和一系列反卷积必要的数据准备，通过反卷积就能得到QSM解，即组织磁化率分布。
3 磁量图反卷积算法
3.1 基于多角度多次扫描数据的QSM算法
此类算法采用多角度多次扫描的数据来重建QSM，主要思想是采用过样技术来弥补单次采集偶极子的零区域数据缺失问题。其代表性算法是多角度采样磁化率计算法(calculation of susceptibility through multiple orientation sampling，COSMOS)[3]。
3.2 单角度QSM算法
虽然多角度多次扫描数据可以很好地约束QSM的逆问题，但是实际临床使用中存在难以操作、有配准误差等问题。临床使用中可以接受的扫描方案是采用单次自然体位的扫描。为了开发能临床使用的单角度QSM方法，学者们提出了多种单角度QSM算法，本文也着重对其进行归类和讨论：
3.2.1 基于偶极子核修正的算法
这类算法对偶极子核绝对值大于某个预设值的区域采用直接除法来得到磁化率分布，对绝对值小于这个值(圆锥面附近，cone region，CR)的区域直接赋于不等于零的值，再在k-space做除法。此类方法中具有代表性的是：基于k-space经验值的直接截断法(truncated K-space division,TKD)[12]。TKD改变了偶极子核本身，在理想情况下仍会带来QSM解的偏离；同时此类算法没有考虑噪声的影响，在k-space除法当中会带来噪声放大作用。
3.2.2 基于偶极子核近似的算法
针对TKD在CR区域采用直接替换原来偶极子核的值的武断做法，研究人员尝试用其他办法来对CR区域做处理。加权k-space导数法(weighted k-space derivative，WKD)[13]在做k-space除法之前，在CR区域对磁化率分布和偶极子核同时做求导处理，使得算法在偶极子核的圆锥面上的求解不再包含零域。但由于在圆锥面每一边，D的符号是一致的，所以D的导数将小于D;这将导致该方法更倾向于在CR区域放大噪声,无法得到最佳估计。
3.2.3 基于通用数学先验条件的贝叶斯算法
通用数学先验条件可以用来约束最小化模型，所以多位学者提出使用通用数学条件约束QSM，从而得到唯一解。这些使用在QSM中的先验条件有：①基于梯度L2范数的空间平滑代价函数[6, 14]；②基于梯度的L1范数的稀疏代价函数[6, 15]；③基于全变差范数的稀疏代价函数[16]；④基于小波域的稀疏代价函数[16]；⑤基于两个稀疏函数混合的代价函数[15-16]。该类方法是将数学上的先验条件直接使用在QSM求解当中，并没有使用磁化率分布自身的特性，而且部分算法并没有考虑噪声不是高斯白噪声的条件。
3.2.4 基于物理结构先验条件的贝叶斯算法
通用数学先验条件并没有很好地体现和限定真实磁化率分布自身的特性。如果可以得到某种磁化率分布的特性，那么就可以用这种分布特性来指导图像重建。此类QSM方法的典型例子就是在QSM重建方法中使用特定物理解剖信息的先验知识[4, 14-15, 17]，代表方法是MEDI[4, 15, 17] 方法。Liu等[18]已经证明在高斯白噪声和给定精确限定下，MEDI方法能得到全局最优解。
4 磁量图研究总结与展望
综上所述，QSM作为一种有别于传统的自旋质子密度、T1、T2加权成像的新组织对比方法，扩展了MRI对磁性物质的检测能力，具有广阔的临床应用前景。QSM已经应用在测量脱氧血红蛋白或含铁血黄素沉积；跟踪氧化铁纳米粒子沉积、量化钆等外源性造影剂；检测钙沉积、测量神经退行性疾病中的铁沉积、脑血管病及脑外伤中的微出血灶、脑内静脉血管畸形[2, 19-22] 等科研和临床试验中。研究人员对QSM各个步骤的研究，尤其是伴随着基于单角度扫描的QSM反卷积方法中的从基于偶极子核的探索式方法到基于物理结构贝叶斯先进方法的发展，使得QSM解的精度越来越高，从而使得QSM的应用范围越来越广。
QSM的进一步发展将依靠先进数学和物理方法的使用，同时也需要研究和讨论有关信噪比、场强、分辨率、回波时间、反演的正则化项。未来的QSM研究将包括以下几个方面：①更精确的三维解卷褶算法研究。虽然当前已经存在针对人类脑组织的稳健二维或者三维解卷褶算法，但对于存在陡峭相位梯度的地方解卷摺仍然会失败。②先进的背景场移除方法研究。当前各种背景场移除方法虽然取得了较好的效果，但是并不能完全移除VOI之外的背景场。所以如何更精确地移除VOI之外的背景场是下一步研究的一个方向。③更精确的反卷积方法研究。当前的研究对反卷积算法的规范化、算法与正则化项的关系，以及反卷积结果与信噪比、场强等扫描参数的关系仍然不清楚。④磁化率对比机制研究。当前研究对磁化率对比在体内存在的机制仍然不清楚；当前研究表明非磁化率敏感效应对观测到的相位和局部场也有贡献，但对这个相关度并不清楚。⑤更精确的磁化率MRI成像序列研究。当前的GRE相位信号的扰动不止单纯由磁化率产生，如何分离其他机制产生的影响需要更精确的相位信号分离方法，这需要设计更先进的MRI序列。
