磁化率是物质的一种物理性质，它反映物质在外磁场中的磁化程度。当磁化率源置于外磁场时，会引起局部的磁场变化。在传统的磁共振成像中，磁化率引起的磁场变化往往被视为图像伪影的来源。然而，磁化率同样也是组织的一种内在特性，如果能有效地加以利用，可以为研究组织的结构和功能提供重要的信息。许多研究表明，利用生物组织内在的磁化率信息，可以有效地对组织的铁含量、钙化、血氧饱和度等进行定量测量^[1,2,3,4,5]。

       磁敏感加权成像（susceptibility weighting imaging, SWI)^[6,7]是利用组织的磁化率信息来增强磁共振图像对比度的技术。传统的磁共振成像大多利用磁共振信号的幅值信息而忽略了相位信息，SWI通过对相位信息进行一定的预处理来获取与组织磁化率信息密切相关的局部场图信息，并将该场图信息与T2^*加权幅值图相结合，提高T2^*加权图的对比度。虽然SWI已得到广泛的发展和应用^[8,9,10]，但由于磁化率引起的场图变化是非局部的，每个点的场图变化是由所有的磁化率源所引起的场图变化的累加，因此，很难直接通过场图信息定位真正的组织磁化率分布，也无法对磁化率进行定量的分析。针对这个问题，近年来发展出定量磁化率成像（quantitative susceptibility mapping, QSM)^{[11,12,13,14,15]}。与SWI相似，QSM也是首先利用磁共振成像的相位信息来获取组织局部的场图变化信息。与SWI不同的是，QSM不是利用场图对幅值图进行加权，而是直接通过场图与磁化率之间的物理关系来反演出磁化率分布图像。QSM已经在脑出血、多发性硬化症和帕金森综合症等多种脑神经疾病中得到了应用^[16,17,18]。

       虽然QSM具有广泛的应用前景，但是定量磁化率图像的重建比较复杂。一般QSM是通过梯度回波序列(gradient echo, GRE）来获取磁共振数据^[19,20]，而重建过程主要可以分为相位图的初始处理、背景场的去除和磁化率反演3个部分。本文将对定量磁化率成像的基本原理和重建流程进行概述，并对重建过程中每个步骤的不同方法进行介绍。此外，本文也将对当前定量磁化率成像的一些主要的临床应用进行介绍。

1　定量磁化率成像基础

       在磁共振成像中，空间分布的磁化率源在外磁场作用下会引起局部的磁场变化，场图与磁化率分布的关系可近似表示为一个偶极核(dipole kernel）跟磁化率分布的卷积^{[13, 21]}：

       其中r表示空间域中的坐标向量，χ(r)表示磁化率的分布，θr表示向量r与外磁场B0方向之间的夹角，b(r)表示由磁化率引起的场图变化，它是相对于外磁场B0测量的，即b(r)=[B(r)-B0]/B0。是描述磁化率与磁场关系的偶极核。空间域中磁场与磁化率的这种卷积关系在傅里叶域中可以表示为：

       其中k为傅里叶域中的坐标向量，kz为沿z方向的坐标，是偶极核在傅里叶域中的表现形式。可以看出，这种表示方式更为简单。同样，利用傅里叶变换也可以将式(1）表示成更简单的形式，即b=F-1DFχ，其中F和F-1分别表示傅里叶变换矩阵和傅里叶反变换矩阵。

       在MRI中，磁化率引起的场图变化随回波时间演化转变成相位信息。在GRE序列中，当回波时间为TE时，相位与场图的关系可以表示为：

       其中φ(r)为空间位置r处的相位，φ0(r）为与感应线圈相关的初始相位，γ为旋磁比。根据式(3)可知，场图b(r)可以通过不同回波的相位值进行拟合得到。由于组织外部（例如空气）也存在磁化率，因此，拟合所获得的场图不仅包括组织内部磁化率引起的局部场图，还包括外部磁化率引起的背景场图，也就是b(r)=bt(r)=bb(r)+bl(r)，其中bt(r)、bb(r）和bl(r）分别表示总场图、背景场图和局部场图。因此，在由b(r)计算组织内部的磁化率分布前，需要先去除背景场信息。

       磁共振采样获得的信号是包含实部虚部的复数数据，例如在GRE序列中，第j个回波的复数信号可用指数形式表示为：

       其中s(r,TEj）和a(r,TEj）分别表示第j个回波的复数信号和幅值，TEj为第j个回波的回波时间。由于反正切函数具有周期性，由式（4)计算得到的相位信息会存在缠绕现象，即：

       其中φw(r)为缠绕的相位，J(r)为缠绕的整数倍数。因此，在拟合得到场图后，需要进行空间的相位解缠绕。

       由式(1)～（5)可知，在磁共振成像中，组织内部磁化率图像的重建是一个逆向的后处理过程。首先从获取的相位图出发，对不同回波的相位图进行拟合，再经过相位解缠绕和背景场去除，得到与组织磁化率相关的场图信息，然后通过磁化率与场图的物理关系反演出最终的组织磁化率。图1给出了定量磁化率重建的基本流程。

点击查看大图

点击查看大图

点击查看大图

点击查看大图

点击查看大图

点击查看大图

图1  定量磁化率重建流程示意图。步骤1：数据预处理，产生幅值图和初始总场图并解相位缠绕；步骤2：背景场去除；步骤3：磁化率反演
Fig. 1  Schematics for quantitative susceptibility reconstruction. Step 1: Data pre-processing to generate magnitude image and unwrapped field; Step 2: Background field removing; Step 3: Susceptibility inversion.

2　数据预处理

2.1　场图拟合

       由式（4）可知，在多回波的情况下，初始的总场图b(r)可以采用以下的非线性最小二乘式子拟合得到^[22]：

       其中Ne为回波数。这个式子一般使用迭代的算法进行求解，例如高斯牛顿法、共轭梯度法等。除了式（6）这种非线性的拟合方法，还有线性的拟合方法^[23]，即先分别求出各个回波的相位信息，再根据不同回波信号的相位和回波时间利用线性最小二乘法进行拟合。线性的拟合方法在求解上更加简单，但由于磁共振相位的噪声分布比较复杂^[24]，因此通过线性方法拟合场图对噪声的抑制较不理想。相比之下，因为磁共振信号在实部和虚部的噪声分布都可以看作是正态分布，因此采用式（6）这种非线性最小二乘拟合可以更好地抑制噪声。

       此外，若由于采样时间等原因的限制只能获取单个回波数据，则无法通过多个回波数据的拟合来得到初始的场图，这时一般假设由感应线圈产生的初始相位φ0(r）为0，然后直接由b(r)= φ (r)/(γ B0TE)获取初始的场图信息。

点击查看大图

2.2　相位解缠绕

       由于相位缠绕的问题在很多领域都会遇到，因此有大量的方法被提出用于解决这个问题^{[25,26,27,28]}。在定量磁化率成像中，比较常用的两种方法是拉普拉斯法^[25]和区域生长法^[27]。在拉普拉斯法中，待求的未解缠绕相位经过拉普拉斯算子的作用后可以表示为缠绕相位的数学形式，即：

       因此可以得到φ= (∇²)^-1 (cos φw∇² sin φw-sin φw∇² cos φw) ，这里的拉普拉斯算子∇²和反拉普拉斯算法(∇²)^-1都可以通过傅里叶变换快速实现，因此拉普拉斯法的求解速度很快。而在区域生长法中，相位是利用相邻点之间相位的连续性进行逐点的解缠绕，从质量最优的相位点出发，根据其周围的相位点与其之间的相位差来增加或减少相应的整数倍的2π，然后逐渐扩展，直到整个区域都处理完成。这里的质量可以有多种衡量方法，一般是根据每个相位点所对应的幅值大小来确定其质量，这是因为幅值越大，其受噪声的影响就相对越小。与拉普拉斯法相比，由于区域生长方法是逐点进行处理，对磁化率这种三维数据的处理速度会慢很多。而拉普拉斯法虽然较快，但其所得到的结果往往过于平滑，因此在一些相位变化较明显的区域其准确性不如区域生长法。

点击查看大图

3　背景场去除

       背景场去除是定量磁化率重建中一个重要的步骤，能否准确的滤除总场图中的背景成分对于后面的磁化率反演十分重要。由于背景场总体较为平滑，所以最简单的方法是通过高通滤波的方式滤除总场图中的低频部分。然而，高通滤波法在滤除背景信息的同时，往往会导致局部磁场中部分低频信息也被滤除。对于组织边界处的高频背景场，高通滤波的效果也不理想。随着定量磁化率成像的发展，很多新的去背景场方法不断涌现^{[29,30,31,32,33]}，其中最典型的是基于磁场物理特性的复杂调和伪影去除法（sophisticated harmonic artifact reduction for phase data, SHARP)以及它的改进方法^{[29,30, 32]}。

       由于背景场来源于组织外部的磁化率，根据麦克斯韦方程可以得到背景场bb(r)在组织内部满足拉普拉斯方程：

       因此在数学上，满足拉普拉斯方程的函数称为调和函数，调和函数具有良好的均值性质，即调和函数bb(r）与非负且径向对称、标准化（即积分为1)的函数ρ(r)做卷积后，其值等于调和函数bb(r)本身，即bb(r) = (ρ ⊗ bb)(r)。利用这个均值性质可以得到总场图在组织内部有：

       也即

       其中δ(r)表示单位冲击响应函数。需要注意的是这个等式只在组织内部成立。利用傅里叶变换，可以把式(11)转换成

       其中M表示组织内部区域的二值掩膜矩阵（内部为1，外部为0)，矩阵C为卷积核(δ -ρ)的傅里叶形式。

       在原始的SHARP方法中^[32]，求解式（12)时将其近似成F^-1CFbl=MF^-1CFbt，然后利用截断的奇异值分解进行求解，求解完成后再乘上掩膜矩阵M，虽然不够精确，但其操作简单。另一种称为正则化的复杂调和伪影去除法（regularization-enable SHARP, RESHARP)^[30]则引入了一个关于场图的先验信息：局部磁场通常远小于背景磁场。利用这个先验信息以及公式（12)，RESHARP方法在求解时使用如下的正则化公式：

       其中第一项为公式（12)的最小二乘拟合形式，第二项则是引入的先验信息，λ表示正则项系数。这个式子求导后可以通过共轭梯度算法进行求解。除了RESHARP，还有另一种利用场图先验信息的去背景场方法叫做拉普拉斯边界值法(laplacian boundary value, LBV)^[29]。与SHARP和RESHARP利用调和函数的均值性质进行求解不同，LBV直接求解偏微分方程（8)或者（9)。由于偏微分方程在求解时需要边界条件，即bl或bb在边界处的值，而这些又是未知的，因此需要利用合理的近似值来代替。根据前面的先验信息知道，总场图中局部磁场的成分远小于背景磁场，因此可以假设在边界处有bb≈bb+bl=bt和bl≈0，有了这些近似的边界条件之后，LBV通过偏微分求解算法中经典的完全多网格算法（full multigrid, FMG)来求解方程（8)或(9)，从而得到背景场图或者局部场图。

       相比高通滤波法，SHARP、RESHARP和LBV这3种方法都是基于相应的物理和数学性质来进行背景场的去除，因此在效果上一般会更好，在定量磁化率重建中的应用也更为广泛。

点击查看大图

点击查看大图

点击查看大图

点击查看大图

点击查看大图

4　磁化率反演

       磁化率反演是定量磁化率重建中最重要的部分。通过式(2）知道，局部磁场在傅里叶域下可以表示成偶极核与磁化率相乘的结果。但由于偶极核D(k)在魔角附近的锥面区域(1 /3 = kz² /|k|²）存在0值，因此无法通过直接相除[即X(k) = b(k)/D(k)]的方式得到磁化率，需要有效的正则化来解决这个不适定的反问题。近年来大量的方法被提出来解决这个问题^{[1, 11, 13, 23, 34,35,36,37]}，其中通过多方向采样来计算磁化率^[34](calculation of susceptibility through multiple orientation sampling, COSMOS）是一种非常有效的方法，可以有效地抑制重建出的磁化率图像中的伪影，但使用COSMOS需要采样多个不同方向的数据，这在临床上是不现实的。单方向的方法在实际应用中更为广泛，大致可以分为两类，一类是基于傅里叶空间的阈值处理方法，一般称为傅里叶域法。另一类则是引入先验知识并利用最优化算法进行求解，一般称为空间域法。这里主要介绍这两类方法中比较常用的两个：阈值截断法（truncated k-space division, TKD)^[11]和总变分法（total variation, TV)^[1]。

       作为傅里叶域法的一种，TKD使用阈值来代替D(k)中那些小于阈值的数，从而使得原来无法直接进行的除法操作得以实现，其具体的表达式如下：

       其中thr表示所选取的阈值，sign(D(k)）为D(k）的正负号。由式(14）所得的χTKD(k)经过傅里叶反变换即可得到空间域的磁化率分布。TKD的优点是非常简单，缺点是不容易找到一个合适的阈值，过大的阈值会使整体的磁化率值变得太小，而过小的阈值又不能起到足够的伪影抑制作用。另一方面，阈值的截断会导致傅里叶域数据的不连续，因此很容易引入其它的伪影。

       TV法是一种空间域的正则化方法，由于组织中的磁化率具有比较平滑的特性，因此TV法引入了梯度的稀疏性这样一个正则项，并利用最优化的解法来重建磁化率图像，其最优化表达式为：

       其中第一项为磁化率与场图的拟合式，第二项为引入的梯度稀疏项，χ*和χ分别表示所求得的磁化率和待求的磁化率值，G表示空间域的梯度算子，包括3个方向，即G=[Gx, Gy, Gz]。参数λ则是用于控制正则项的权重。式(15)可以通过非线性共轭梯度算法进行迭代求解。由于引入了平滑的先验信息，因此TV法所获得的磁化率图像在伪影的抑制上要大大优于TKD等一些傅里叶域方法。但由于需要利用迭代的方式进行求解，其重建速度远慢于TKD。另一方面，TV法由于是对整个图像的梯度进行抑制，因此未能考虑组织内不同结构之间磁化率的差异，所得的磁化率图像往往在结构边界处比较模糊。针对这个问题，基于形态学相似性的反演方法^[35] (morphology enabled dipole inversion, MEDI）可以很好地解决。MEDI的基本原理与TV法相似，可以视为一种改良的TV法。MEDI的最优化表达式为：

       与原始的TV法相比，MEDI在两个范数项中各引入一个权重矩阵，其中Wd用于衡量不同场图数据点的可靠性，防止场图的噪声或误差传递到后面的磁化率图像中。Wd一般与幅值图成正比，这是因为幅值越大，其所对应的场图点受噪声的影响相对越小，可靠性越强，这与前面相位解缠绕中质量图的选择是类似的。Wg为是梯度权重矩阵，只在非边界处取1，而在边界处取0，因此可以将TV中的平滑性约束限制在不同结构的内部。Wg一般也是通过幅值图计算得到，这是因为不管是幅值图还是磁化率图，其生理解剖结构是一致的，因此幅值图中的边界位置等价于磁化率图中的边界位置。通过这两个权重矩阵，MEDI不仅抑制了磁化率图像中的伪影，也提高了图像在组织边界处的清晰度，从而大幅提高磁化率图像的质量。

点击查看大图

点击查看大图

点击查看大图

5　定量磁化率成像的临床应用

       随着QSM重建技术的发展，目前QSM已经在一些疾病的检测和诊断上得到了应用，受篇幅所限，本文只介绍其中主要的3种应用。

5.1　在脑部微出血检测上的应用

       由于对铁血黄素沉积引起的不均匀场的敏感性，T2^*加权的GRE图像和SWI是目前广泛使用的两种用于检测脑部微出血的技术，这些图像在出血点的地方往往表现为低信号。然而，低信号区与出血点的这种相关关系会随所用的成像参数的变化而改变，例如，对于同样一个出血点，不同的回波时间所产生的信号衰减不同，因此低信号区的形状和大小也会有所不同。这种对参数的依赖影响了脑出血检测的准确性。相比之下，QSM作为一种测量脑部微出血的方法，其所得到的定量磁化率图像可以有效克服对成像参数的依赖性^[16]。

5.2　在组织钙化等抗磁性物质检测上的应用

       与幅值图、相位图等常用MRI对比图相比，QSM的一个优点是它可以轻易地将抗磁性的组织钙化与顺磁性的物质（例如血铁黄素）区分开。在GRE幅值图中，钙化和慢性出血的地方都表现为低信号，因此通常很难区分低信号区域是由哪一种病变引起。而GRE的相位图虽然可用于检测抗磁性的钙化，但其诊断的准确性尚未得到很好的证明。另一方面，CT是一种广泛用于检测组织钙化的成像技术，但与MRI相比，CT对人体具有一定的放射性损害。随着QSM技术的发展，QSM可望替代CT成为检测组织钙化的一项无损技术。目前已经有临床研究表明，QSM在检测组织钙化方面优于一般的相位成像，且具有非常高的敏感度(90%）和精确性（95%)^[5]。

5.3　在基于铁含量异常的神经疾病上的应用

       帕金森病与多巴胺能细胞消亡和黑质致密部的铁积累有关系。最近的研究表明，由于对铁含量变化的敏感性，定量磁化率成像可用于辅助帕金森综合征的诊断和治疗^{[18, 38]}。与正常人相比，帕金森病病人在黑质致密部的磁化率值升高。相比R2和R2^*图，QSM对于帕金森病病人和正常人的区分更加敏感。因此，QSM可以为帕金森病的研究提供一种非常有效且更加定量的方法。同样，在多发性硬化病人的基底神经节和病灶区域，铁分布通常会异常增高，并且随病灶时间和炎症状态的变化而发生改变，因此QSM可用于测量多发性硬化病人脑部病灶和非病灶区的磁化率变化。最近的QSM研究表明^[17]，在最初的几年（大约4年），多发性硬化的病灶相对于表现正常的脑白质保持较高的磁化率，随着病变时间延长（大于7年），逐渐回归到与表现正常的脑白质相同的磁化率。相比R2^*图像，QSM对组织病变的敏感度更高，因此也可以用于疾病的早期预测。

6　总结及展望

       由于定量磁化率成像具有广泛的医学应用前景，近年来关于定量磁化率图像重建方法的研究十分热门。本文主要介绍了定量磁化率成像的基本原理以及重建流程，并阐述了每一个重建步骤目前使用的一些主要方法。同时也介绍了定量磁化率成像当前一些主要的临床应用。由于定量磁化率重建是一个连续的后处理过程，每一个步骤的处理结果都会对下一个步骤产生影响，前一个步骤的错误很容易在后续的步骤中被放大。因此，了解重建的每个步骤并选择合适的方法，不仅有助于高质量定量磁化率图像的获得，也有助于为临床应用提供更加准确有效的检测和诊断信息。