磁共振定量磁化率成像（quantitative susceptibility mapping，QSM）通过梯度回波序列（gradient recalled echo ，GRE)来获取复数数据，然后利用相位信息重建得到局部磁化率分布图。QSM发展至今，其已经在脑出血、多发性硬化症和帕金森综合征等多种脑神经疾病研究中得到了应用^[1,2,3,4]。

       QSM的重建过程主要包括场图拟合、背景场去除和磁化率反演三个步骤^[5,6,7]。准确的场图可以提高磁化率分布图的图像质量，使得磁化率测量更加准确，核团显示更清晰，进而提升QSM应用的可靠性。场图可以使用一个回波的相位图直接计算出来^[8]，也可以使用多个回波的相位图进行拟合得到^[9,10]。有研究表明，使用5~ 10个回波拟合出的场图更准确，重建得到的磁化率分布图信噪比更高，并且可以获得较好的颅内核团磁化率对比度^[10,11,12]。

       对于多回波场图拟合，参与拟合回波数的增多会增强数据的抗噪能力，从而提升场图计算的准确性。然而，并不是回波数目越多场图拟合就越准确，这主要受到两个方面的限制：一是由于回波信号强度受T2弛豫的影响，回波数目越多势必定会导致后面的回波信号强度降低，有时甚至降至噪声水平；二是由于噪声、部分感兴趣区（如颅底）磁化率不均匀、血液流动、主磁场本身的波动等因素，每获得一幅相位图都会伴随着误差，相位误差过大会影响场图拟合的准确性。多回波场图拟合可以采用线性拟合或非线性拟合方法，而非线性拟合方法相较线性拟合法可以更好地抑制噪声、减少解缠绕误差的传递，是一种更加鲁棒的拟合方法^[13,14]。

       为了减少低信噪比的数据点对场图拟合结果的影响，我们先前的研究提出了一种截断加权最小二乘法(weighted linear least-square，WLS)，有效改善了颅底部位QSM的图像质量^[15]。但该算法采用模图来间接判断数据点相位的可靠性，有时判断不太准确。因此笔者提出了一种多回波场图拟合的改进算法，其通过场图拟合的残差来判断数据点相位的可靠性，以提高场图拟合的准确性，从而提高磁化率分布图的图像质量。

1　材料与方法

1.1　算法原理

1.1.1　非线性场图拟合

       通过多回波梯度回波序列获得的磁共振信号可以表示为下述公式：

       其中，r为空间位置，TEj为回波时间，a(r, TEj)为每个回波信号的幅值，γ为旋磁比，γΔB(r）为所求场图。采用非线性最小二乘法拟合可以表述为下式^[5]：

       其中N为回波个数。这个公式一般可采用高斯牛顿法或共轭梯度法进行迭代求解，迭代初始值可对前3个回波进行线性拟合来获得。

点击查看大图

点击查看大图

1.1.2　改进的场图拟合算法

       在多回波的情况下，数据的可靠性会对拟合结果的准确性有很大影响，因此很有必要在拟合前对相位数据的可靠性进行检查。笔者提出的场图拟合改进算法对回波相位可靠性进行了检查，主要包括以下几个步骤。

       第一步，对所有回波的相位进行时间域解缠绕^[13,14,16]。

       第二步，对前三个回波相位进行线性拟合得到初始相位φ0′和角频率ω′，具体可表述为如下线性方程^[17]：

       其中，φ0为初始相位，ω为角频率，TEj为第j个回波的回波时间，φj为第j个回波的相位。

       第三步，奇异点判断。采用拟合得到的初始相位和角频率计算前三个回波的拟合相位与测量相位之间的平均差值，如下式：

       设置一个误差阈值th，如果>th，则认为该点相位误差严重，在场图上为奇异点，对奇异点之后再做进一步处理。

       第四步，非奇异点回波可靠性判断。利用下式计算第j(N≥j> 3)个回波的拟合相位与测量相位的差值δj：

       如果δj>th，则认为该回波相位可靠性低，舍弃该回波及其后面的所有回波信号。

       第五步，非奇异点场图拟合。利用公式（2)对非奇异点剩余的回波信号进行拟合得出场图^[5]。

       第六步，奇异点场图处理。一般认为磁化率随空间位置是缓慢变化的，则空间各位置的局部磁场也应该是缓慢变化的，所以可利用奇异点周围的可靠点的局部磁场的平均值来替代前面场图上奇异点的局部磁场。

点击查看大图

点击查看大图

点击查看大图

1.2　数据采集

       此次研究共招募了15名(8名女性，7名男性）健康被试，平均年龄为（24.5±1.5)岁，实验前均要求被试仔细阅读磁共振检查知情同意书，在获得被试同意后进行磁共振扫描。其中5名被试[3名女性，2名男性，平均年龄为（25.8±1.3)岁］在西门子3.0 T磁共振成像系统Trio Tim上完成扫描，其信号接收采用12通道头线圈；另外10名被试[5名女性，5名男性，平均年龄为（23.8± 1.1）岁］则在西门子3.0 T磁共振成像系统Prisma上完成扫描，其信号接收采用20通道头线圈。两台不同仪器上设置的扫描序列和扫描参数相同。QSM扫描采用三维多回波梯度回波序列，读出梯度为单极梯度，具体参数为：重复时间（TR)60 ms，第一个回波时间（TE1) 6.8 ms，回波间隔时间(ΔTE)6.8 ms，回波数目(N)8，翻转角（FA) 15° ，视野（FOV)240 mm×180 mm，矩阵(matrix)384×288，体素大小（voxel size)0.625 mm×0.625 mm×2.000 mm，层数（slices)96，采样带宽(readout bandwidth)240 Hz/pixel，并行成像加速因子为2。为了评价本研究算法的降噪效果，对所有被试连续进行了两次梯度回波数据采集。

1.3　数据处理

       QSM重建包括如下步骤：首先，由不同回波的原始模图求均方根得到组合模图；其次，对多回波的原始相位图分别采用常规非线性拟合方法和本文提出的改进算法拟合得到场图；然后，采用区域生长法对场图进行空间解缠绕^[18]；再者，采用偶极场投影(projection onto dipole fields，PDF）方法去除背景场^[6]；最后，采用基于形态学的偶极子反演（morphology enabled dipole inversion，MEDI)方法反演得到磁化率分布图^[7]。

1.4　数据分析

       为了检验改进的多回波场图拟合算法是否能提高QSM重建结果的图像质量，对采用常规非线性拟合方法和本文提出的改进算法所得到的磁化率分布图的噪声水平进行对比评价。在磁化率分布图上利用ITK-SNAP(www.itksnap.org）软件在核团显示最清晰的层面上分别勾画出左右双侧的黑质、红核、尾状核、苍白球、壳核等感兴趣区（region of interest，ROI)，如图1所示。将连续两次扫描得到的磁化率分布图相减得到磁化率差图，然后将前面勾画出的ROI加载到磁化率差图上，并计算出每个ROI内的噪声。一个ROI内的噪声定义为在磁化率差图上该ROI内所有体素的标准偏差^[11,19]。

       采用非参数配对样本Wilcoxon秩和检验方法检验实验组和对照组的磁化率分布图的噪声水平差异(P<0.05表示差异具有统计学意义），所用分析软件为SPSS 20.0。

点击查看大图

图1  颅底部位（A）及基底节部位（B）深部核团ROI示意图。SN：黑质；RN：红核；CN：尾状核；GP：苍白球；PU：壳核；R：右侧；L：左侧
图2  一名在Trio Tim系统上完成扫描的健康被试颅底部位采用常规算法和改进算法得到的场图以及对应的磁化率分布图。A：采用常规算法得到的场图；B：采用改进算法得到的场图；C：采用常规算法得到的磁化率分布图；D：采用改进算法得到的磁化率分布图
图3  1名在Trio Tim系统上完成扫描的健康被试基底节部位采用常规算法和改进算法得到的场图以及对应的磁化率分布图。A：采用常规算法得到的场图；B：采用改进算法得到的场图；C：采用常规算法得到的磁化率分布图；D：采用改进算法得到的磁化率分布图
Fig. 1  The ROIs on the susceptibility map at the skull base level (A) and basal ganglia (B). The ROIs were drawn manually in the right (R) and left (L) red nucleus (RN), substantia nigra (SN), caudate nucleus (CN), globus pallidus (GP) and putamen (PU).
Fig. 2  Field maps (A&B) and susceptibility maps (C&D) at the skull base level from a healthy subject with the conventional field fitting algorithm (A&C) and improved field fitting algorithm (B&D). The subject was scanned on the 3T Trio Tim system.
Fig. 3  Field maps (A&B) and susceptibility maps (C&D) in the basal ganglia from a healthy subject with the conventional field fitting algorithm (A&C) and improved field fitting algorithm (B&D). The subject was scanned on the 3T Trio Tim system.

2　结果

       图2显示了1名在Trio Tim系统上完成扫描的健康被试采用常规场图拟合算法和改进算法后得到的其颅底部位的场图（图2A、B）以及对应的磁化率分布图（图2C、D)。采用常规场图拟合算法得到的场图在核团区域（红色箭头所指）和一些非核团区域（黄色箭头所指）有明显的伪影，磁化率分布图的对应区域也非常模糊，无法很好地区分黑质和红核；而采用改进算法后得到的场图更平滑，所对应的磁化率分布图的组织边界更清晰，更容易区分出黑质和红核。图3显示了一名在Trio Tim系统上完成扫描的健康被试基底节部位在场图拟合时采用常规算法和改进算法后得到的场图（图3A、B）以及对应的磁化率分布图（图3C、D)。从图中可以看出，两种场图拟合算法得到的场图在基底节部位区别较小，对应的磁化率分布图均能很好地显示尾状核、苍白球和壳核等核团，而采用改进算法得到的磁化率分布图在一些非核团区域（黄色箭头所指）显示更平滑，伪影更少。

       在Trio Tim系统上完成扫描的5名被试各个灰质核团在采用两种拟合算法得到的磁化率分布图上的噪声平均值与标准差如表1所示。在场图拟合时采用改进算法后，颅底的核团黑质（双侧）、红核（右侧）与基底节的核团苍白球（左侧）的噪声水平显著下降(P<0.05，差异具有统计学意义），下降幅度为27%～67%；颅底与基底节部位的其他核团的噪声水平均有所下降，颅底的红核（左侧）与基底节的苍白球（右侧）的噪声水平下降幅度较大，但差异不具有统计学意义。

       在Prisma系统上完成扫描的10名被试各个灰质核团在采用两种拟合算法得到的磁化率分布图上的噪声平均值与标准差如表2所示。在场图拟合时采用改进算法后，颅底的核团黑质（双侧）、红核（右侧）的噪声水平显著下降(P<0.05，差异具有统计学意义），下降幅度为15%～35%；颅底与基底节部位的其他核团如红核（左侧）、尾状核（右侧）、苍白球（左侧）的噪声水平略有下降，尾状核（左侧）、苍白球（右侧）、壳核（双侧）的噪声水平略有上升，但差异均不具有统计学意义。

点击查看表格

表1  Trio Tim系统上完成扫描的被试采用常规算法和改进算法得到的磁化率分布图上各个核团的噪声水平对比
Tab. 1  Susceptibility noise (mean±standard deviation) of ROIs in deep gray matter nuclei from the subjects scanned on the Trio Tim system with conventional field fitting algorithm and improved field fitting algorithm

点击查看表格

表2  Prisma系统上完成扫描的被试采用常规算法和改进算法得到的磁化率分布图上各个核团的噪声水平对比
Tab. 2  Susceptibility noise (mean±standard deviation) of ROIs in deep gray matter nuclei from the subjects scanned on the Prisma system with conventional field fitting algorithm and improved field fitting algorithm

3　讨论

       针对在QSM重建过程中多回波拟合场图时某些部位（如颅底）原始相位数据线性度较差的情况，笔者提出了一种改进的多回波场图拟合算法。对比采用常规算法和改进算法的磁化率分布图噪声水平后发现：对于在Trio Tim系统上采集的数据，改进算法可明显提高重建后的磁化率分布图的图像质量（图2，图3)，颅底的黑质（双侧）、红核（右侧）和基底节的苍白球（左侧）的噪声水平显著下降，其他如红核（左侧）、尾状核（双侧）、苍白球（右侧）和壳核（双侧）的噪声水平虽有所下降，但差异不具有统计学意义（表1)；对于在Prisma系统上采集的数据，磁化率分布图的图像质量没有明显提升，但从噪声水平定量对比结果可以看出，颅底的黑质（双侧）、红核（右侧）的噪声水平显著下降，其他如红核（左侧）、尾状核（右侧）、苍白球（左侧）的噪声水平略有下降，尾状核（左侧）、苍白球（右侧）、壳核（双侧）的噪声水平略有上升，但差异均不具有统计学意义（表2)。

       通过对两组不同磁共振成像系统（Trio Tim和Prisma)上的数据的结果分析发现，两组数据结果都显示在颅底的黑质（双侧）、红核（右侧）的噪声水平显著下降，并且在其他核团区域噪声水平无显著上升，这说明本文提出的改进算法对Trio Tim和Prisma系统上的数据都有良好的降噪效果。Trio Tim系统上的数据结果显示磁化率分布图的图像质量有明显提升，并且苍白球（左侧）的噪声水平显著下降，在黑质（双侧）、红核（右侧）和苍白球（左侧）的噪声水平下降幅度为27%～67%，而Prisma系统上的数据结果显示相应区域的噪声水平下降幅度为15%～35%，这说明仪器性能（如主磁场的稳定性、均匀性）的提升增强了数据采集时的抗噪能力。

       改进的场图拟合算法通过对原始相位数据线性度的检查，舍去部分或者全部不可靠的回波信号，提高了场图拟合的准确性。场图拟合的准确性降低了磁化率分布图的噪声水平，有效提升了图像质量。而我们先前提出的截断WLS算法虽然有效提高了颅底部位的磁化率分布图的图像质量，但磁化率分布图上基底节部位的深部核团（尾状核、苍白球、壳核）的噪声略有上升^[15]，因此笔者改进的算法在提高QSM图像质量方面更加全面。

       在评价算法的可靠性方面，本文仅对比研究了采用常规算法和改进算法所得到的磁化率分布图的噪声水平差异，结果显示改进算法能显著地降低噪声。但本文并未对磁化率值的准确性进行相应比较，其主要原因是：年龄和性别会影响人体大脑深部核团的磁化率^[20]，而且QSM测量得到的磁化率还受到扫描参数不同^[21]和重建方法不同^[22]等因素的影响，因此磁化率测量值目前还没有一个黄金标准。

       本文提出的改进算法还有以下两点不足：一是判断原始相位数据线性度的阈值是通过多次试验得出的经验值，当采用不同的扫描方案时，该阈值是否适用仍有待研究；二是奇异点的相位计算方法在组织边界处的适用性较差。理论上来说，阈值越小则对原始相位数据的线性度要求越高，场图拟合应该越准确，但如果原始相位数据本身线性度较差，则会产生大量的奇异点，而奇异点的相位计算方法又要求其周围的可靠点数尽可能多，所以阈值不宜太小；然而，阈值也不能太大，阈值过大会降低对相位线性度的约束，不能有效地抑制相位误差过大的点，期望在以后的研究中能给出一个客观、具有普遍适用性的阈值。对奇异点的处理是建立在磁化率缓慢变化的假设上的，但在组织边界处磁化率可能出现突变，此时奇异点的相位计算方法不再适用，期望在以后的研究中给出更可靠的解决方法。

       综上所述，本文提出的场图拟合改进算法通过对场图拟合前原始相位数据线性度的检查和奇异点的判断，提高了场图拟合的准确性，降低了磁化率分布图的噪声水平，提高了图像质量，使得核团的勾画更加容易，更有利于神经退行性疾病铁沉积方面的研究。