CN108693491A

CN108693491A - 稳健的定量磁化率成像系统和方法

Info

Publication number: CN108693491A
Application number: CN201810291820.9A
Authority: CN
Inventors: 王乙; 刘喆; 奇永旭; 阿列克谢·迪莫夫; 文焱; 张景威; 帕斯卡·斯宾瑟梅里
Original assignee: Cornell University
Current assignee: Cornell University
Priority date: 2017-04-07
Filing date: 2018-04-03
Publication date: 2018-10-23
Anticipated expiration: 2038-04-03
Also published as: US10890641B2; US20230160987A1; US11635480B2; CN108693491B; US20180321347A1; US20210132170A1

Abstract

提供示例性定量磁化率图的方法，系统和计算机可访媒介从复数磁共振成像数据来生成组织磁性特性的图像，通过使用贝叶斯推断方法最小化一个由数据保真项和两个正则化项组成的成本函数。数据保真项直接由复数磁共振成像数据构成。第一个先验是从已知形态的匹配结构或信息内容构建的。第二先验是从具有近似均匀且已知的磁化率值和解剖图像上的特征的区域构建的。定量磁化率图可以通过最小化成本函数来确定。因此，根据该示例性实施例，可以提供用于确定与至少一个结构相关联的磁化率信息的系统，方法和计算机可访媒介。

Description

稳健的定量磁化率成像系统和方法

相关申请的交叉引用

本申请要求了2017年4月7日提交的美国临时申请62/482,864号的优先权，该美国临时申请在此被全文引用。

技术领域

本发明涉及磁共振成像，并且更具体地涉及使用磁共振成像中收集的信号来定量地对物体固有物理特性成像。

背景技术

磁共振成像(MRI)中的定量磁化率图(QSM)已经获得越来越多的临床和科学兴趣。QSM在表征和量化化学成分如铁，钙和造影剂(包括钆和超顺磁性氧化铁(SPIO) 纳米粒子)方面表现出前景。这些化合物的组织成分可能在各种神经系统疾病中发生改变，如帕金森病，阿尔茨海默氏病，中风，多发性硬化症，血色素沉着症和肿瘤以及整个身体其他疾病。QSM能够解开与磁性敏感化率(或简称为磁化率)相关的新信息，即基础组织的物理性质。由于生物体中铁和钙的普遍存在，它们在细胞功能中的积极参与，它们在肌肉骨骼系统中的重要作用，QSM通常对通过跟踪循环系统中和代谢活动的的铁来研究铁/钙的分子生物学非常有用，使用铁和钙作为替代标记。因此，铁，钙和造影剂引起的磁化率的精确成像将为临床研究人员探索人体的结构和功能以及临床医生更好地诊断各种疾病和提供相关治疗提供巨大的益处。

发明内容

以下描述用于收集和处理受试者的MRI信号的系统和方法的实现，以及重建受试者固有的物理特性(例如，磁化率)的影像。在一些实现中，对应于受试者的MR信号数据可以被转换成定量描绘受试者的结构和/或组成和/或功能的基于磁化率的图像。使用该定量磁化率图，可以生成受试者的一个或多个基于磁化率的图像并将其显示给用户。然后，使用者可以将这些图像用于诊断或实验目的，例如，至少部分地在图像上研究患者的结构和/或组成和/或功能调查和/或诊断和/或治疗各种病症或疾病。由于与其他磁化率成像技术相比，所描述的实施方式中的一个或多个可以导致基于磁化率的图像具有更高的质量和/或准确度，所以这些实施方式中的至少一些可以用于改善用户对受试者结构的理解和/或组成和/或功能，并可用于提高任何由此产生的医学诊断或实验分析的准确性。

下面描述的一些实施方式可用于评估磁偶极子反演，同时允许并入预处理，各向异性加权，原始对偶解法，基于化学组成的脂肪补偿，和强化已知区域的近似恒定磁化率，以改善磁敏度图质量，其包括抑制同质区域中的人造信号变化和提供绝对磁化率值，这些是确保定量准确性，实验重复性和再现性所迫切关注的。

通常，这里公开的本发明的一个方面通过参考具有已知的同质化易感性的组织区域，例如脑中脑脊液(CSF)的零敏感性和体内完全充氧的主动脉，来实现组织磁化率的绝对定量。本发明的这一方面解决了纵向和交叉中心研究以及需要绝对量化的研究中零参考的实际问题。这里公开的发明还能够减少由与偶极场模型不兼容的场数据所引起的阴影伪影，此外还减少了由数据中的粒状噪声引起的拖尾伪影。

通常，这里公开的本发明的另一方面是能够进行准确的组织结构匹配，在对最佳匹配的数值搜索中使用取向依赖性L1范数和其导数的原始-对偶计算。所公开的预处理全场反演方法绕过背景场去除步骤，其稳健性可能受到相位解缠，颅骨剥离和边界条件的影响。

通常，这里公开的本发明的另一方面是实现精确的成像方法，通过校正传统T2 和T2^*方法的脂肪，纤维化和水肿的混淆影响来量化肝脏铁，通过结合脑血流量来量化脑氧代谢率，量化身体中的骨矿化和其他钙化，和量化心室血液氧合水平和心肌内出血。

这些方面的实施可以包括以下一个或多个特征。

在一些实施方式中，通过应用R2^*阈值和连通性，基于低R2^*值将脑室内的脑脊液分割。

在一些实施方式中，关于磁化率分布的先验信息的确定是基于L1范数和从所获取的物体图像而估计的物体结构信息。

在一些实施方式中，所述物体包括大脑的皮质，壳核，苍白球，红核，黑质或丘脑底核，或腹部中的肝，并且产生所述物体的一个或多个图像包括产生描绘皮层，壳核，苍白球，红核，黑质或丘脑底核或腹部肝脏的一幅或多幅图像。

在一些实施方式中，所述物体包括多发性硬化病灶，癌性病灶，钙化病灶或出血中的至少一个病灶，并且生成所述物体的一个或多个图像包括生成一个或多个图像描绘所述多发性硬化病灶，癌性病灶，钙化病灶或出血中的至少一个病灶。

在一些实施方式中，操作还包括基于一个或多个图像量化在对比增强的磁共振成像过程中引入到物体的造影剂的分布。

在一些实施方式中，组织化学成分是脂肪和水，并且使用适当的初始化从复数MRI数据迭代地计算体素中的磁场，水分和脂肪分数。

根据本发明的一个方面，提供了一种用于给一个物体生成一个或多个图像的方法，法包括：

接收由磁共振扫描仪获得的磁共振成像数据，其中磁共振成像数据是复数的并且包括：

与物体组织有关的幅度和相位信息或实部和虚部信息，以及高斯噪声；

根据组织成分对复数磁共振成像数据建模来计算物体组织感受的磁场；

根据一个图像特征识别组织中的具有近似均匀磁化率的感兴趣区域；

基于计算的磁场来估计物体的磁化率分布，其中物体磁化率分布的估计包括：

确定成本函数，至少涉及与磁共振成像数据中的高斯噪声的似然函数相关联的数据保真项，与关于磁化率分布的一个结构先验信息相关联的第一正则化项，和与强制组织感兴趣区域内均匀磁化率相关联的第二正则化项；

基于最小化成本函数来确定一个磁化率分布；和

在显示装置上呈现基于所确定敏化率分布而生成的物体的一个或多个图像。

根据本发明的另一个方面，提供了一种用于对物体产生一个或多个图像的系统，包括：

处理器；

通信地耦合到处理器的图形输出模块；

通信地耦合到处理器的输入模块；

编码有计算机程序的非暂时性计算机存储介质，所述程序包括指令，所述指令在由处理器执行时使所述处理器执行包括以下步骤的操作：

基于最小化成本函数来确定一个磁化率分布；和

在附图和下面的描述中阐述了一个或多个实施例的细节。从说明书和附图以及权利要求中，其他特征和优点将显而易见。

附图说明

图1展示了不同CG迭代步数(Iter)后使用传统PDF+LFI，未预处理的TFI以及预处理的TFI重建的人体脑部QSM。为了得到同PDF+LFI的50步迭代结果相似的QSM 图像，未预处理的TFI需要300步迭代。预处理过程将TFI所需迭代数目降至50。

图2显示了仿真实验中利用PDF+LFI,LBV+LFI以及TFI测量0.1ppm点源磁化率的误差分布图。当点源靠近ROI边界时，使用TFI的误差要显著低于PDF+LFI的误差。

图3中的(a)是仿真实验中使用的真实磁化率图。右侧图片显示使用不同预处理权重P_B时的拟合误差(b)及重建磁化率与真实磁化率之间的误差(c)。

图4展示了在体模实验(a)和人体实验(b)中，通过最小化重建磁化率与标准值之间的误差来确定最优的预处理权重P_B。

图5展示了使用COSMOS，PDF+LFI，SSQSM，微分QSM以及本文提出的TFI重建的健康人体QSM。同SSQSM和微分QSM相比，PDF+LFI与TFI能够更完整的保留上矢状窦，在脑后部边界处尤为明显(由空心箭头标明)。在大脑底部区域(由实心箭头标明)，TFI的重建结果较PDF+LFI更加平滑。

图6展示了图5中使用COSMOS，PDF+LFI，TFI，微分QSM以及SSQSM重建结果中针对苍白球，壳核，丘脑及尾状核测量的磁化率平均值与标准差。这五种方法对壳核，丘脑及尾状核的测量结果较为相近；同种方法相比，SSQSM明显低估了苍白球的磁化率。

图7展示了对一例脑出血患者的重建结果。从上到下分别为：T2^*加权幅值图，使用PDF+LFI，SSQSM，微分QSM以及本文提出的TFI重建的QSM。

图8从上到下依次展示了幅值图，整体场，使用TFI重建的脑部QSM，以及分别使用微分QSM及TFI重建的全头部QSM。从左到右依次为的轴向，矢状位和冠状位图像。全头部TFI结果清晰地描绘出头骨，窦腔空气以及皮下脂肪的分布。利用示例 ROI测量的磁化率均值为：鼻腔空气7.38ppm，头骨-1.36ppm，脂肪0.64ppm。

图9是在理论部分中定义的各向异性权重矩阵P_⊥ξ对相对于对应幅度图像(绿色)的梯度的χ梯度的垂直分量进行惩罚。因此，它促进了两个向量之间的平行定向。

图10(a)是二值各向同性加权(左起第一列)；图10(b)是各向异性权重(左起第二至第四列)。对于各向同性加权，非对角线分量全部为零。值显示在[-1， 1]的范围内。因此，图10(b)中的对角线和非对角线分量分别在[0,1]和[-1,1] 的范围内。

图11是计算的磁化率图。图11(a)康奈尔是MRI研究实验室数据(λ/2＝1000)；图11(b)是QSM 2016重建挑战数据(λ/2＝235)。左侧的第一列，第二列和最后一列分别是轴向，矢状和冠状视图。在各向异性加权中，在轴向视图中表现为斑点的条纹(拖尾)伪影以及矢状和冠状视图中的条纹基本上被抑制。这些可以在整个图像中看到(全局)。箭头也放置在伪影显眼的地方。请注意，一旦放大图像，就会更清楚地看到伪影。

图12显示了通过全体指标评估QSM准确性。将RMSE(第一列,图12a)，HFEN (第二列,图12a)，SSIM(第三列,图12b)和DoA(第四列,图12b)与正则化参数(各向异性权和各向同性权重)。(顶部第一行)关于COSMOS的康奈尔MRI研究实验室数据。(从上至下第二行)QSM 2016重建挑战与COSMOS相关的数据。(排名第三)QSM 2016重建挑战数据相对于χ₃₃。

图13显示了通过基于本方法和参考标准之间的ROI的线性回归评估QSM准确性。斜率和R²值与正则化参数(各向异性加权和各向同性加权)作图。(顶部第一行) 关于COSMOS的康奈尔MRI研究实验室数据。(从上至下第二行)QSM 2016重建挑战与COSMOS相关的数据。(排名第三)QSM 2016重建挑战数据相对于χ₃₃。

图14显示了具有中心差分格式(GNCG-C)的高斯-牛顿共轭梯度(GNCG)，具有正向差分的GNCG(GNCG-F)和关于基础真值的原始-对偶算法(PD)的收敛行为。尽管PD与GNCG-F一样快且∈＝10^-6，但GNCG不能达到PD所能达到的精度(实和虚线之间的差距)。

图15显示了具有不同窗口级别和输入相对差异场(RDF)的真值(GT)位于第一行。GT由分段平滑区域(白质)和分段恒定区域(其余)组成。第二行显示来自 PD，GNCG-F和GNCG-C的计算解。最后一行显示GT与不同算法计算结果之间的差值图像(GT图像中对应黄色边界框的特写图像)。

图16显示了与COSMOS相比，来自GNCG-F和PD的离体数据的计算的磁化率图。

图17显示了通过线性回归和布兰德-奥特曼图对不同算法与COSMOS之间进行精度评估。这些图是从离体脑块中产生的。

图18显示了与COSMOS相比，来自GNCG-F和PD的体内脑数据的计算的磁化率图。

图19显示了通过线性回归和布兰德-奥特曼图对不同算法与COSMOS之间进行精度评估。这些图是从体内大脑数据集中产生的。

图20显示了(a)如果适当选择减小因子γ，则qNewton连续方法有助于加速算法。注意，当γ≥0.3时，从黑色曲线上方的曲线可以看出，算法收敛较慢。(b) 我们观察到所有的虚曲线最终都比黑曲线的相对误差更低。原因在于每种连续方法 (不同的γ)的∈最终值略小于10^-6(1.0×10^-7，6.25×10^-7，1.0×10^-7, 5.90×10^-7,7.63×10^-7,9.13×10^-7在图例底部)。在第3500次迭代中，我们计算了相同的顺序的相对误差：8.93％，8.78％，8.83％，8.77％，8.84％，8.85％和8.88％。请注意，这些数字都大于PD算法产生的数字：8.64％和8.66％。虽然 PD显示出快速收敛，但在500次迭代后会出现超调。这种行为可以通过调整步长来控制，如图14所示。

图21显示了仿真结果：第一行展示了标准图 χ₃₃(a)，以f_aniso为输入分别使用MEDI(b)和QSM0(c)重建的QSM，和以 f_iso为输入分别使用MEDI(d)和QSM0(e)重建的QSM。第二行展示了图(b-e)中各QSM与χ₃₃的差值图(f-i)。使用各向异性模型产生的场作为输入时，QSM0减少了黑质与丘脑底核附近的低信号阴影状伪影(箭头所示)。

图22的(a)是脑室分割展示：脑室脑脊髓液自动分割过程的示例：先对(第一列)进行阈值化生成二值掩码图(第二列)，然后通过分析连通关系提取脑室脑脊髓液区域。另外该图还显示了使用不同R值时，RMSE的曲线(b)及幅值与QSM图像 (c)。分割后的CSF ROI被覆盖至幅值图上。RMSE在R＝5时达到最小。

图23显示了使用MEDI及QSM0重建的两例MS患者的脑部QSM。低信号伪影在 QSM0结果中得到抑制(箭头所示)。

图24的(a)显示了针对全部23处病灶相对于NAWM的QSM测量结果的散点图(左图)与布兰德-奥特曼(右图)分析结果。图24(b)显示了增强和非增强性新MS 病灶的MR图像。A)T1w+Gd，B)T2w和C)患有复发缓解型MS的44岁女性的QSM。在T1w+Gd中发现两个增强病灶(A，B，箭头)。一个是壳增强(A，白色箭头)，另一个是结节增强(A，黑色箭头)。壳增强病灶出现轻度QSM高信号(C，白色方块)，结节状出现QSM等信号(C，黑色方块)。D)T1w+Gd，E)T2w和F)QSM在一名35岁复发缓解型MS患者中的表现。与前6个月前的MRI相比，在T1w+Gd 和T2w图像中发现两个新的非增强性病灶(D，E，箭头)。两个病灶均出现QSM高信号，伴有明亮的边缘(F，箭头)。(c)相对于NAWM磁化率预测病灶增强状态的接受操作者特征曲线的。引导模型的AUC＝0.9594和来自折刀交叉验证的ROC1的 0.9530。

图25显示了在脂肪水钆模型上测量的磁化率和R2^*图。a-d：测量的磁化率图； e-h：测量的R2^*图；i：测得的溶液总体磁化率值与变化的Gd浓度的关系；j：脂肪贡献(在零Gd的小瓶中测量)与溶液总磁化率对变化的脂肪浓度的关系；k：脂肪校正的磁化率对Gd浓度；l：测量的R2^*值与变化的Gd浓度的关系。PDFF，质子密度脂肪分数。

图26显示了来自患有肝脏血管瘤的患者的T2加权图像(a)，水图(b)，脂肪图 (c)，无脂肪校正的磁化率图(d)，脂肪校正的磁化率图(e)和R2^*图(f)。病灶 (箭头)显示脂肪含量降低(c)，这大大降低了R2^*(f)，但在QSM(d对e)上效果得到了很好的补偿。

图27显示了来自疑似肝癌患者的T2加权图像(a)，水图(b)，脂肪图(c)，无脂肪校正的磁化率图(d)，脂肪校正磁化率图(e)和R2^*图(f)。正常肝脏中的大量脂肪(c)增加R2^*(f)，但不影响QSM(d对e)。肿瘤(红色圆圈)脂肪少(c)，因此其R2^*(f)很少。

图28显示了具有肝动脉化疗栓塞介入治疗后的多发性肝转移瘤患者的T2加权图像(a)，水图(b)，脂肪图(c)，无脂肪校正的磁化率图(d)，具有脂肪校正的磁化率图(e)和R2^*图(f)。转移瘤残留的碘油(红色圆圈)被视为脂肪增加(c) 并强烈增加R2^*(f)，但对QSM影响不大(d对e)。

图29显示了在肝脏局灶性病变患者中，(a)肝脏脂肪百分比与R2^*测量之间存在线性相关性，但(b)肝脏磁化率值与脂肪百分比之间没有相关性。

图30显示了水，脂肪，局部场，磁化率和R2^*图来自三个不同程度的铁超载的受试者(a-o)。对于轻度铁负荷过重的患者，左叶也存在纤维化，但不影响QSM，但增加R2^*(第二列箭头)。在这组患者中，无肝纤维化的肝脏区域(实心正方形)肝脏磁化率与R2^*之间存在线性相关性(p)，但肝纤维化区域(实心园形)则不存在线性相关性(p)。

图31(a)显示了在没有病变的ROI中肝脏磁化率与R2^*之间存在线性相关性；图 31(b)显示了病灶内的ROI没有线性相关性。

图32的上图显示了获取一个超短回波和一个常规回波图像的双回波梯度回波采集。两个回波在连续的TR之间移动，以获得四个独特的回波时间值(在相同TR中采集的回波使用相同的线型显示)。图32的下面的图像比较两个获取的回波的幅度和相位(在线圈组合期间第一个回波的相位被设置为零以消除不同通道的相位偏移)。

图33显示了QSM与CT对猪蹄体模ROI值的线性回归。在估计的磁化率和亨氏单位值之间观察到良好的相关性。

图34显示了比较CT图像(a,c)与重建QSM(b,d)的对应平面和超短回波幅度 (e,f)。注意皮质骨的均匀抗磁性显示以及高HU和低磁化率值区域之间的总体上好的对应关系。

图35显示了来自健康志愿者股骨的磁场图重建结果(a)，磁化率(b)和超短回波幅度图(c)。皮质骨的整体均匀抗磁性被成功显示。注意场图和QSM(*)中的可见小梁，以及股四头肌肌腱的强抗磁性(**)。有关所选区域的更多详细信息 (虚线)，请参阅图36。

图36显示了比较脂肪和水分数图，场图，QSM和不同信号模型的估计衰减率，包括单峰和多峰脂肪谱以及常见和纯水R2^*建模。

图37显示了使用所提出的技术重建膝关节QSM(右)的局部场(左)和薄板MIP。对股骨，胫骨和腓骨皮质区域实现了良好的描绘。还要注意股骨中小梁和骨骺线(*) 的描述。重建后的QSM患有膝关节周围的模糊。

图38显示了使用MLV方法的三个随机选择的受试者的平均L曲线。λ(方程 42)被选择为500。

图39显示了使用咖啡因激发和MLV方法从健康受试者的轴向，矢状和冠状切片中3D CMRO2图的实例。相应的反转准备SPGR T1w图像显示在右侧。这两种方法的CMRO2图显示出良好的一致性，并具有良好的灰白质对比度，与反转准备的SPGR T1w图像一致。

图40显示了使用基于刺激(左)和基于MLV的(右)方法重构的第二个受试者的OEF和CMRO2图的示例。

图41显示了布兰德-奥特曼图比较了使用基于刺激和基于MLV的方法重建的CMRO2和OEF图。在所有比较中检测到小偏差(两次测量平均值的<4％)，没有统计学意义(p>0.05)。

图42显示了来自2DBH方法的QSM与来自3DNAV方法的QSM具有类似的RV对LV 对比度。当由于屏气不一致(如箭头所示)导致注册错误时，由此产生的QSM有重大伪影，难被可靠的理解。由于2DBH中的SNR比3DNAV中的低，所以来自2DBH的 QSM看起来比来自3DNAV的QSM有更多噪音。

图43显示了用于组织磁化率成像过程的示例。

图44显示了计算机系统框图的示例。

各附图中相似的附图标记表示相同的元件。

具体实施方式

以下描述用于收集和处理受试者的MRI信号的系统和方法的实现，以及重建受试者固有的物理特性(例如，磁化率)的图像。一些公开的实施方式克服了定量敏感性映射(QSM)中的限制，其中计算速度和精度以及贝叶斯先验信息影响QSM图像质量和实际使用。在一些实施方式中，为了重建受试者的磁化率图像找到最佳地适合可用场数据和组织结构先验信息的磁化率分布，使用预处理和原始对偶优化方法来加速迭代计算过程，并且使用原始对偶优化方法来精确迭代计算过程。在一些实施方式中，使用各向异性边缘权重和使用侧脑室中脑脊液(CSF)的接近零磁化率的附加信息来改进用于估计磁化率的贝叶斯方法中的组织结构先验信息，所述附加信息可以自动地从大脑图像分割出。这些改进的实施使得QSM技术稳健且准确，将QSM 的使用从大脑延伸至包括颅骨，骨骼，心脏和肝脏，以及量化与有氧能量代谢相关的脱氧血红素。

QSM基于最小化适合于测量的组织磁场b(r)(单位B0)与偶极场的成本函数，将b中的噪声近似为高斯，数据保真项是

该积分形式b＝d*χ允许数字估计磁化率解，经常在与先验信息相关的正则化下。从麦克斯韦方程的物理学中可知，积分形式可以代表其统治偏微分方程(PDE)的解。对于磁化率解的数学分析，可直接从麦克斯韦方程和洛伦兹球校正得到的等效PDE 能提供了有用的见解：

方程2是关于磁化率χ(z轴是时间)的波动方程和关于场b的拉普拉斯方程。在没有任何误差的连续空间中，b(r)＝d*χ(r)，B(k)＝D(k)X(k)，可以获得没有拖尾的精确解：如果D(k)＜ε/2，η_∈＝1；否则为0。任何不能与偶极子场拟合的场数据都称为不相容场数据v(r)，它可能来自各种原因，包括噪声，数字化误差和各向异性源。这个不相容场会导致伪影,根据波传播解：

这里是波传播器，在不相容源上下的魔角锥上具有较大的值。包括噪声，离散误差和各向异性源的偶极不相容部分产生由波传播器以B0场方向为时间轴定义的伪影。颗粒噪声导致拖尾和连续错误导致拖尾和阴影伪影。

这里描述的发明利用上述数学事实来改善QSM速度和准确度。拖尾伪影的特征是边缘沿着在k空间中魔角(54.70)或图像空间中的互补魔角，这些几乎与组织边缘完全不同。因此，它们可以通过基于组织结构信息的正则化项在数值优化期间被最小化，以识别没有拖尾条纹的解，其是形态使得偶极子反演(MEDI)方法。使用边缘梯度和L1范数对组织结构信息进行数值评估。该计算先前通过将L1范数的导数近似为连续函数并忽略边缘梯度方向来执行。使用原始对偶优化方法并根据边缘梯度方向使用各向异性加权，可以在准确性和速度两方面改善该近似阴影伪影具有较低空间频率内容，因此可以通过引入低空间频率的成本函数项来最小化。此外，可以使用特定的组织结构和阴影结构信息来制定惩罚项。一种特殊的组织结构是脑室内的CSF几乎是纯水，具有很少细胞内含物。因此，脑室磁化率图应该接近均匀，并且任何与均匀性的偏差都应该被视为阴影伪影，在数值优化过程中通过结合正则化项而受到惩罚。一个具体的阴影结构信息是在大范围内平滑分布的误差源所产生的阴影，例如在区域上方和下方强烈的白质各向异性分量，遍布整个区域，并且由低空间频率的分量组成。因此，可以使用k²加权(拉普拉斯算子)作为预处理，在数值优化期间避免克雷洛夫空间中的低空间频率，对阴影伪影进行惩罚，或者使用正则化项惩罚低空间频率，其相应空间波长是物体大小和半物体大小，其成本函数与估计的白质各向异性贡献相关。在另一实施方式中，根据脑图谱模拟的作为偶极不相容部分的白质各向异性效应被调整为使得分割的脑室区域中的CSF磁化率均匀，同时使脑QSM中的阴影伪影最小化。

根据可用的磁场数据和组织先验信息寻找最佳的磁化率分布的迭代过程可以使用预处理来极大地加速。这对于搜索大范围的可能磁化率值特别有用。结果，本发明中描述的QSM系统和方法的改进实现了快速，准确和稳健的QSM，改进了大脑中的应用，并且使得能够应用于大脑外的其他器官。

稳健的QSM改善了许多应用，包括脑QSM具有阴影和条纹伪影的减少以及均匀的脑室CSF作为自动和一致的零基准，全脑QSM无颅骨剥离和骨矿化，其中骨抗磁化率与CT衰减系数成正比，可用于正电子发射断层影像的衰减校正，实现不含R2/R2^*混淆误差的肝铁含量成像，大脑氧耗量代谢率成像，包括氧合水平和心肌内出血在内的心脏磁化率成像，以及血管壁斑块磁化率成像来分解斑块钙化和斑块内出血。

1.用于QSM的预处理整体场反演方法

在这一节中，我们描述了一种针对QSM的改进。该改进算法使用预处理加快QSM 的收敛速度，并对背景场去除和组织场反演两部分进行了整合。

1.1.介绍

定量磁敏感图(QSM)可以从梯度回波信号相位提取组织磁化率的空间分布。目前的QSM方法包括两个步骤：i)背景场去除以确定由组织产生的局部场，以及ii) 从局部场反演到组织磁化率。其对中央大脑区域，特别是对于中脑核中的铁沉积的磁化率重建较为稳健。但是，依然存在着若干技术挑战。

一个主要的挑战是由于背景场去除方法中固有假设造成的背景和组织场之间不精确的分离。由于组织-空气交界存在显著的磁化率差异，这个问题在大脑边界附近尤为明显。为了避免背景场和局部场的各自拟合，有学者基于前向信号方程的偏微分形式提出了基于拉普拉斯算子的QSM方法。拉普拉斯操作隐式地消除了背景场。然而，拉普拉斯算子的实际实现需要在针对误差增大的稳健性以及可视皮质脑组织的完整性之间进行权衡。对大脑区域的必要的侵蚀可能会影响大脑边界处的某些结构的可视化，例如上矢状窦。

QSM问题中另一个挑战是在感兴趣区域(ROI)内磁化率的动态变化范围可以很大，这通常导致条纹伪影。例如，脑出血(ICH)与周围组织之间的磁化率差异可能超过1.6ppm。使用ICH信号的非线性QSM模型可以减少这些条纹伪影，但不会消除它们。针对这一挑战，近期研究通过分离强(如ICH)和弱(正常脑组织)磁化率来源的拟合过程，从而防止与ICH有关的伪影渗入正常大脑。然而，这些方法需要仔细选择正则化参数或用于检测ICH的阈值以最小化后续弱磁化率反演中的伪影。

在这项工作中，我们使用适合具有大动态范围的磁化率的广义逆问题，来解决这两个挑战。我们提出了一个预处理QSM来执行整体场反演(TFI)，大大减少了与不精确背景场去除相关的误差传播，并且抑制了脑出血QSM中的条纹伪影。

1.2.理论

对于QSM，整体磁场通常被分解为两个分量：局部场f_L被定义为由给定感兴趣区域M内的磁化率χ_L(局部磁化率)产生的磁场，而背景场f_B被定义为由M以外的磁化率χ_B(背景磁化率)产生的磁场:

f＝f_B+f_L＝d*(χ_B+χ_L) [4]

这里*是卷积算子，d是由洛伦兹球校正的单位偶极子产生的场。分量f_B和f_L可以独立估计：χ_B或f_B的估计被称为背景场去除。以往人们提出了各种背景场去除方法，诸如高通滤波(HPF)，偶极子场投影(PDF)或基于拉普拉斯算子的方法。之后通常使用正则化，从局部场f-f_B获得χ_L。这一步被称为局部场反演(LFI)。背景场中的错误会传播到随后的局部场中，导致最终磁化率中的误差。尽管近期研究允许在局部场反演期间更新背景场，但为了首先估计背景磁化率，它需要先验的磁化率图谱以及共同配准。

这里，我们提出了整体场反演(TFI)算法来同时估计χ_B与χ_L:

这里使用与传统QSM反演问题相同的形式。这里χ＝χ_L+χ_B表示整个图像空间中的磁化率分布。数据加权项w可以通过计算从MRI数据到整体场中的误差传播，从幅度图像导出。这里使用了加权总变差项来抑制由偶极子核d中零值的存在而引起的条纹伪影。

迭代优化算法，如共轭梯度(CG)在求解方程5时的收敛速度会十分缓慢，如图1中健康脑部结果所示。方程5在χ＝0处线性化并使用CG求解。在早期迭代过程中，部分背景场被拟合为局部场，产生不合理的局部磁化率图(图1d)。为了解决这个问题，我们使用先验辅助的预处理来改善收敛性。如果解是具有均值和协方差矩阵的高斯随机向量，则使用近似协方差矩阵Γ的右预处理矩阵P(见下文)，即P^HP≈Γ，会提高收敛速度。对于TFI，二元对角预处理矩阵P构造如下:

其中i表示体素，M表示组织ROI。它被设计成使得体外和体内M体素之间的矩阵P^HP的差异近似等于局部和背景区域(包括骨骼和空气)之间的磁化率差异。我们进一步改进P，通过对R2^*图进行阈值化来提取M内的强磁化率(如脑内出血)，如果假设具有高R2^*的体素具有强磁化率。因此，P被定义为:

预处理的整体场反演问题变为：

最终磁化率为χ^*＝Py^*。方程8由高斯-牛顿算法求解，具体过程如下:

这里F为傅里叶变换,为k空间偶极子核而G为有限差值算子。在第 n步牛顿迭代，Φ在当前y_n处被线性化:

使用弱微分，可以通过求解如下方程获得dy_n:

或者:

而这可以使用共轭梯度(CG)方法求解。注意到我们使用且∈＝10^-6来避免除0。终止条件为：迭代数到达100 或相对残差低于0.01。当相对变化||dy_n||₂/||y_n||₂小于0.01时终止高斯-牛顿方法。

最近提出的一些方法类似地避开背景场去除，直接从整体场中估计磁化率。这些方法基于信号方程的偏微分形式(方程4):

这里，Δ为拉普拉斯算子。注意到拉普拉斯操作去除了来自背景源χ_B的贡献。其中两种方法为：

i)单步QSM(SSQSM):

这里M_G为从幅值图中获得的二值掩码。该问题可以使用预处理的共轭梯度算法来高效地求解。

ii)微分QSM:

这里和

基于拉普拉斯的方法实现了在一步内重建QSM。然而，拉普拉斯算子的实际实现需要对ROIM进行侵蚀。侵蚀的程度取决于用于计算拉普拉斯的核函数的宽度。

1.3.材料与方法

为了开发和评估所提出的预处理整体场反演方法(TFI)，我们进行了数值仿真，体模实验，使用COSMOS作为参照的健康受试实验以及脑出血的患者的脑部实验。我们探索了全脑QSM。

方程7中的P_B是利用经验性的最优情况确定，以确保在有限的CG迭代内，高斯- 牛顿求解器能够达到相对于标准磁化率的最小误差。然后这个P_B值被用于其他类似的数据集。在这项工作中，我们使用了我们将所提出的TFI方法在体模和人体实验中的性能与SSQSM和微分QSM进行比较。在这项工作中，我们使用具有可变半径的球形平均值(SMV)来实现SSQSM和微分QSM中的拉普拉斯算子，其中核直径从d变化到3个体素。两种方法中的ROI掩模也相应地被侵蚀。如之后所述，ROI最中心处的核直径d会被优化。

数值仿真

PDF+LFI和TFI的准确性。我们构建了一个大小为80×80×64的数字脑模型。除0.1ppm的单点磁化率源之外，脑组织磁化率设定为0ppm。这里将背景磁化率设为0，并使用前向模型f＝Md*χ计算整体磁场。点源χ_S的磁化率值通过LFI及未预处理(P_B＝1)的TFI分别进行估计，并与真值χ_ST(0.1ppm)进行比较。对于顺序反演方法，PDF和LBV被用于进行背景场去除，而MEDI被用于局部场反演。对于每种方法，通过在整个ROI上移动点源来重复该过程，从而生成每种方法误差图来显示估计误差|χ_S-χ_ST|/χ_ST随点源位置的空间变化。对于PDF+LFI，LBV+LFI 和TFI，正则化参数λ设为10^-3。

预处理矩阵P对TFI的影响.我们构建了一个大小为256×256×98的数字脑模型，并模拟了不同脑组织的磁化率：白质(WM)-0.046ppm,灰质(GM)0.053ppm, 苍白球(GP)0.19ppm,及脑脊髓液(CSF)0。为了模拟ROI外部的空气，背景磁化率设为9ppm。然后使用前向模型，从真实磁化率分布χ_T(图3a)计算整体磁场，并在此基础上添加高斯白噪声(SNR＝200)。我们在TFI方法中使用了4种不同的预处理矩阵(Eq.6)：P_B＝1(未预处理)，5，30及80。在重建过程中，我们通过拟合残差||M(f-d*χ)||₂以及重建磁化率M_χ与真实值MχT之间方均根误差(RMSE)来衡量收敛性能。这里n是M内的体素数目。正则化参数λ设为10^-3。

体模实验

为了检测TFI方法的准确度，我们将5个装有钆溶液的橡胶球置于1％的琼脂糖体模中。每个球体内标准磁化率χ_ref依次为0.049,0.098,0.197,0.394及0.788 ppm。测量得到的磁化率值以背景琼脂糖作为基准。我们在3T(美国GE)上使用多回波梯度回波序列扫描此体模，并使用未经流动补偿的单方向(前/后方向)笛卡尔采集方式。采集参数为:FA＝15°,FOV＝13cm,TE1＝3.6ms,TR＝71ms,#TE ＝8,ΔTE＝3.5ms,采集矩阵＝130×130×116,体速大小＝1×1×1mm³,带宽＝±31.25kHz，总扫描时间13min。

接下来，我们在同一组体模数据中使用PDF+LFI,TFI,SSQSM及微分QSM，然后使用基于图分割的相位解缠绕方法及非线性场估计方法来生成整体磁场。我们优化了方程7中的P_B。通过针对一定范围的P_B值重复TFI方法，并在第300次CG迭代结束时最小化每个气球的平均估计磁化率χ与其对应标准磁敏度χ_ref之间的误差正则化参数λ通过最小化PDF+LFI的 Err_phantom来确定，其随后也被用于TFI。对于SSQSM和微分QSM，通过最小化 Err_phantom选取λ和可变半径SMV的最大核心直径d。

健康受试的脑部实验

我们在3T(美国GE)上使用多回波梯度回波序列采集了5位健康受试的脑部数据，并使用未经流动补偿的单方向采集方式。本文中所有工作经过了本机构审核会的批准。采集参数为:FA＝15°,FOV＝24cm,TE₁＝3.5ms,TR＝40ms,#TE＝ 6,ΔTE＝3.6ms,采集矩阵＝240×240×46,体素大小＝1×1×3mm³,带宽＝ ±62.5kHz，总扫描时间9min。对每个受试扫描3个不同方向，使用COSMOS重建的QSMχ_Lref作为标准。对于PDF+LFI,TFI,SSQSM及微分QSM,仅一个方向上的数据被采用。ROI由BET算法自动得到。我们使用其中一位受试的数据，通过最小化第 300次CG迭代结束时，估计的磁化率χ_L与标准值χ_Lref之间的RMSE来确定最优的P_B，并用于其他4名受试的数据。CSF 被取作本文中人体试验的基准组织。我们测量了苍白球(GP),壳核(PT),丘脑(TH) 及尾状核(CN)的ROI中的磁化率值用于量化分析。我们计算了每个受试中的相对误差其中χ_GP为该受试GP的平均磁化率值。正则化参数λ通过最小化PDF+LFI中Err_invivo来确定，并用于TFI。对于SSQSM和微分QSM，通过最小化Err_invivo选取λ和可变半径SMV的最大核心直径d。

脑出血病患的脑部实验

我们采集了18例患者的人体数据，每个患者都有脑出血(ICH)，包括癌性病变的脑瘤，以及钙化病变。我们使用多回波梯度回波序列在3T(美国GE)利用未经流动补偿的单方向采集方式进行采集。采集参数为：FA＝15°,FOV＝24cm, TE₁＝5ms,TR＝45ms,#TE＝8,ΔTE＝5ms,采集矩阵＝256×256×28～52, 体素大小＝1×1×2.8mm³,带宽＝±31.25kHz，并行成像R＝2.扫描时间正比于扫描层数(约10片/min)。对每一例患者，我们使用了PDF+LFI,TFI,SSQSM及微分QSM来估计含有脑出血的QSM。对于TFI,我们使用方程4并使用之前COSMOS 人体实验中的P_B。R2^*阈值设为用来在此初步研究中区分出血部位与周围的脑组织。正则化参数λ通过L-曲线确定并用于PDF+LFI,SSQSM及微分QSM。其中TFI与PDF+LFI使用相同的λ。而SSQSM和微分QSM分别使用同先前COSMOS 实验相同的最大核心直径d。为了量化ICH附近的伪影，我们计算了PDF+LFI与TFI 结果中，每处ICH附近5mm之内区域的平均磁化率值。非ICH区域的标准差(SD)降低值也被用于考量ICH相关的伪影减弱的程度。这里，标准差降低值定义为 R＝(SD(χ_PDF|M_non-ICH)-SD(χ_TFI|M_non-ICH))/SD(χ_PDF|M_non-ICH),其中 χ|M_non-ICH代表非ICH区域的磁化率。

人体试验：全头部实验

由于本文提出的方法不再分离背景场去除和局部场反演，因此可以通过将大脑外部的组织(例如头皮，肌肉和口腔软组织)包含到方程8中的ROI掩码M中来生成整个头部的磁化率图。我们在3T(美国GE)上使用多回波梯度回波序列，利用单方向三维辐射采集方式获得了一个健康的受试大脑数据，以取得较大的覆盖空间。这里未使用流量补偿。采集参数为：FA＝15°，FOV＝26cm，TE1＝1ms，TR＝34ms， #TE＝9，ΔTE＝3ms，采集矩阵＝256×256×256，体素尺寸＝1×1×1mm³，带宽＝±62.5kHz，投影数量＝30000，采集时间为15分钟，重建方法为重网格化。TFI 中使用先前COSMOS实验确定的P_B。方程8中的ROI掩码M通过对幅度图像I进行阈值化来确定：M：＝I＞0.1max(I)。我们使用基于图分割的相位解缠和化学位移补偿方法来估计整体磁场，因为相位缠绕是以2π为单位的整位运算，并且图分割技术最适合处理这种离散整位操作。水脂分离问题可以分两步进行估算：1)通过将体素成分简化为水或脂肪，生成相位解缠，磁化率和化学位移的初始估计，2)将初始估计值输入至完整的水脂信号模型通过数值优化来微调磁化率和水脂比例。由于水脂分离不是一个凸问题，基于图分割的这种首步处理为水脂数值优化过程提供了一个稳健的初始化，进而收敛到一个合理的解。

此ROI也应用在微分QSM中。我们同样使用TFI重建了仅包含脑部的QSM作为比较对象，其中掩码使用BET算法获得。我们在脑部TFI，微分QSM及全脑TFI重建中，使用L-曲线分析选择正规化参数λ。

本文中的所有计算都在一台具有6核CPU(Intel Core i7)和32GB内存的台式机的MATLAB环境中完成。

1.4.实验结果

数值仿真

PDF+LFI和TFI的准确性。在图2中，除了ROI的边界外，PDF+LFI方法显示的估计误差小于10％：当点源距边界距离小于4个体素时，PDF+LFI误差至少为 40％。相比之下，LBV+LFI和TFI的最大误差在整个ROI(包括边界)中均为4.8％。

预处理矩阵P对TFI的影响。图3b显示，对于P_B＞1，与未预处理(P_B＝1) TFI相比，测量和估计的整体场之间的差异下降得更快。图3c显示，对于足够大的 CG迭代数目(>1000)，对于所有P_B值，预处理的求解器均可以收敛于相对于标准磁化率误差Err_simu＜0.04。另一方面，在100步CG迭代时，同P_B＝5或80相比，P_B＝ 30实现了较小的Err_simu(图3c)。

体模实验

正则化参数λ在PDF+LFI及TFI中为2×10^-2,在SSQSM中为5×10^-1，在微分QSM中为5×10^-3。最大内核直径d在SSQSM中为13，在微分QSM中为3。最优P_B为10(图4a).各球的估计磁化率(y)与标准磁化率(x)间的线性回归结果：PDF+LFIy＝0.978x-0.004(R＝1.000)；TFI y＝0.973x+0.004(R＝1.000)；SSQSMy＝0.833x+0.019(R＝0.999)；微分 QSMy＝0.971x+0.004(R＝1.000)。Err_phantom为:PDF+LFI:0.0158,TFI: 0.0050,SSQSM:0.0350,微分QSM:0.0010。

健康受试的脑部实验

图1显示了在不同的CG迭代步数中使用PDF+LFI，未预处理的TFI(P_B＝1) 和预处理的TFI(P_B＝30)重建的脑QSM。使用预处理(P_B＝30)的TFI，在50次 CG迭代时生成的结果与使用PDF+LFI获得的局部组织磁化率相近，而未预处理的 TFI需要300次迭代才能达到类似的结果。通过最小化Err_invivo确定P_B＝30(图4b)。平均重建时间：PDF+LFI为232秒，TFI为270秒，SSQSM为50秒，微分QSM为 284秒。对于PDF+LFI和TFI，正规化参数λ：PDF+LFI和TFI为10^-3，SSQSM为 5×10^-2，微分QSM为1×10^-3。SSQSM的最大内核直径d为17，微分QSM为3。

图5显示，所有方法的QSM都与靠近大脑中心的COSMOS结果一致。同时，与PDF +LFI相比，TFI改善了大脑底部QSM的平滑性(图5中的实线箭头)。使用PDF+LFI 和TFI比使用SSQSM和微分QSM能够更好地显示上矢状窦，特别是在大脑后部皮层 (图5中的空心箭头)。图6显示了深部灰质合团的磁化率测量结果。使用COSMOS 作为参考，与微分QSM，PDF+LFI或TFI相比，SSQSM明显低估了GP的磁化率。这种低估现象也发生在其他受试数据中。与COSMOS测量结果相比，SSQSM对GP磁化率的平均相对低估量为27.4％，而PDF+LFI，TFI和微分QSM分别为5.4％，9.6％和8.9％。

脑出血病患的脑部实验

所有出血患者脑部图像的TFI重建中，P_B＝30。PDF+LFI的平均重建时间为 328秒，TFI为325秒，SSQSM为50秒，微分QSM为344秒。对于PDF+LFI和TFI，正规化参数λ为10^-3，SSQSM为2×10^-3，微分QSM为5×10^-4。SSQSM的最大内核直径d为17，微分QSM为3。图7显示的例子中，与使用PDF+LFI，SSQSM或微分 QSM相比，本文提出的预处理TFI方法在ICH部位周围减少了阴影伪影。特别的，我们观察到GP结构在去除阴影伪影后更加显著(如图7中箭头所示)。考虑到ICH周围的阴影伪影表现为负磁化率，ICH附近的平均磁化率的增加预示着伪影的减少。标准偏差的降低同样印证了阴影伪影的减少。

人体试验：全头部实验

将整个头部QSM与幅值/相位图像和脑QSM进行比较(图8)。对于脑部的TFI，重建时间为13.4分钟，微分QSM为19分钟，全头部TFI为14.2分钟。对于脑部TFI 和微分QSM，正则化参数λ均为10^-3，而对于全头部TFI为2.5×10^-3。结果显示全脑TFI QSM的脑内部分与仅脑部的TFI一致，尽管后者在脑的顶部具有稍好的平滑性，如在矢状和冠状位视图中所见。同时，全头部QSM提供了更多的脑外结构的磁化率分布，例如颅骨，鼻窦和皮下脂肪。由于ROI是通过对阈值图像进行确定的，全头部QSM中脑边界处同脑部TFI的QSM相比更显著，尤其是在脑干和小脑部位。另外，全头部QSM清晰地区分了头骨和鼻窦空气，而由于这两个区域中的MR信号丢失，幅度图像难以对他们进行区分。如图8所示，利用鼻窦空气，颅骨和脂肪的示例ROI，平均磁化率的测量值为：蝶状窦7.38ppm，颅骨-1.36ppm，脂肪0.64ppm。这些磁化率值与文献一致。当比较微分QSM和全头部TFI时，脑内部分有着类似的磁化率图，而微分QSM无法显示颅骨和鼻窦空气的磁化率。

1.5.结论

我们的数据表明，针对QSM的整体场反演(TFI)消除了对背景场消除和局部场反演(LFI)的分别处理的需要，并且预处理可以加速TFI收敛。与传统的PDF+LFI 方法相比，TFI在相似的重建时间内，在高磁化率区域附近能够域提供更稳健的QSM，例如空气或含铁血黄素的出血部位。我们还证明了TFI能够在不需要脑分割的情况下生成整个头部的QSM。

背景场去除和LFI通过顺序连接的方式，使用相同的麦克斯韦方程处理相同的数据，其中需要假设或正则化来区分背景场和局部场。例如用PDF进行背景场去除时，需要假设由所有可能的局部单位偶极子场f_dL构建的希尔伯特空间L与由所有可能的背景单位偶极子场f_dB构建的希尔伯特空间B之间具有正交关系。但是，当局部单位偶极子接近ROI边界时，这种正交性假设会失效，进而导致PDF中的误差。由任何其他背景场去除方法引入的类似错误会传播到随后进行的LFI，并产生不准确的局部磁化率。然而，这种利用麦克斯韦方程将MRI数据分离为背景场去除和LFI两个单独问题的方式是不必要的。我们提出的TFI消除了这种分离并避免了相关的错误传播。仿真(图2)和体内(图5)结果显示了TFI相对于PDF+LFI的改进，特别是当局部源接近ROI边界时。注意到LBV+LFI在分离局部和背景场方面也优于PDF+ LFI，但为了避免将ROI边界处的相位噪声引入拉普拉斯算子，LBV需要精确的ROI 掩模，这对于体内大脑QSM可能是具有挑战性的。

预处理对于通过TFI实际性能是必要的。图1表明，对于大脑QSM，未预处理的 TFI需要300次CG迭代才能收敛到与PDF+LFI相似的解，但后者只需要50次迭代。在这里我们介绍了一种针对TFI问题的先验增强的右预处理矩阵(方程7)。类似的先验增强的预处理矩阵被用于MR动态成像中，其中右预处理矩阵被构造为平滑滤波器，反映了线圈灵敏度通常具有空间平滑性。在QSM问题中，我们通过在预处理矩阵(方程7)中将较大的权重(P_B>1)分配给强大的磁化率源(例如空气，头骨或出血)，达到区分强磁化率和弱磁化率源(例如正常脑组织)的目的。这个预处理矩阵P旨在通过P^HP≈Γ逼近协方差矩阵Γ，这里假设方程5中的解χ是高斯随机向量 χ～N(0，Γ)。克雷洛夫子空间迭代求解器(如CG)的收敛性可以使用这个预处理矩阵来改进。这在图1中显示的结果中得到了证实，其中预处理将TFI所需的CG迭代次数减少了6倍。所提出的预处理矩阵不同于另一预处理矩阵，即将局部场反演问题中的系统矩阵通过对角矩阵来近似，并在求逆之后用作预处理矩阵。因此，可能在 TFI中将两种预处理矩阵结合起来实现进一步加速。

由于我们的预处理矩阵是针对大动态范围磁化率设计的，它直接适用于脑出血(ICH)患者。如图7所示，预处理的TFI有效地抑制了出血部位附近的伪影，并提高了ICH的QSM质量。以前在处理大范围磁化率方面的工作基于分段常数模型–显式或者隐式地使用强边缘正则化。在这些方法中，需要为强磁化率和弱磁化率仔细选择不同的正则化参数。否则，欠正则化会导致出血部位附近的条纹伪影，或者过度的正则化会在弱磁性区域中牺牲精细细节。我们提出的TFI方法利用预处理而不是多层正则化来改善收敛速度，减少解中的伪影。这消除了组合来自多个局部场反演实例的QSM的需要。另一方面，所描述的预处理可应用于LFI以抑制ICH的QSM 中的出血位点附近的伪像，类似于预处理的TFI。对于全头部QSM(图8)，我们提出的预处理TFI能够同时产生脑内(弱磁化率)和脑外(强磁化率)的QSM，其中由于过度正则化，脑内部分欠佳。

基于拉普拉斯算子的方法，如SSQSM和微分QSM使用背景场f_B的谐波性质(ROI 内的Δf_B＝0)来消除相位解缠和消除背景场，从而通过单个步骤估计局部磁化率成为可能。通过省略信噪比权重及使用L2范数的正则化，SSQSM能够实现进一步加速。然而，基于拉普拉斯算法的方法会导致脑部区域的侵蚀：拉普拉斯算子使用有限差分算子或球形核算子来实现，二者都需要侵蚀ROI掩模。因此，可以在SSQSM和微分QSM中观察到上矢状窦的侵蚀(图5)。PDF+LFI和TFI无需涉及拉普拉斯算子，避免了这种侵蚀问题。此外，与微分QSM，PDF+LFI和TFI相比，SSQSM在体模和体内实验中对磁化率值有着严重的低估。对于文中使用的正则化参数λ附近范围内的一系列值也观察到了这种低估(结果未显示)。原因可能是L2范数正则化。与L1 范数正则化相比，L2正则正则化被证明低估了磁化率。此外，微分QSM无法估计颅骨和鼻窦空气的磁化率(图8)。这与以前的文献是一致的。原因在于，ROI是通过对幅度图像进行阈值化确定的，因此不包括被认为是“背景”的颅骨和鼻窦空气。因此，拉普拉斯操作去除了由颅骨和鼻窦空气产生的磁场。相反，TFI保留它们产生的场并描绘它们的磁化率分布(图8)。

本文通过全头部QSM实验证明了该方法定量描绘抗磁性和顺磁性材料的能力。类似地，这可应用于动脉粥样硬化成像以解析和定义在颈动脉中发现的斑块中的钙化 (具有负磁化率的抗磁性)和斑块内出血(具有正磁化率的顺磁性)，冠状动脉，主动脉和其他血管壁。本文中的多回波梯度回波3D序列可被修改，以包括脑电门控，从而减少心动周期中脉动变化的影响，另外可以包括流量补偿以最小化流动产生的相位。斑块钙化和斑块内出血成分的QSM将解决钙化和出血混杂的问题。这两者在幅度图像中均呈现低信号。这种对斑块内出血的测量会是决定斑块易损性的主要因素。

全头部QSM提供了颅骨钙化空间分布的量化，而QSM与CT上测量的衰减系数被显示成线性比例。在QSM上测量的这种颅骨和其他颅内钙化可用于校正磁共振引导聚焦超声(MRgFUS)中的超声衰减和散射，并校正MR/PET图片的伽马衰减，同时形成MRI和正电子发射断层扫描(PET)的图片。MRgFUS和MR/PET两者都需要基于骨骼空间定量分布的校正，这仅用标准MRI是不可行的，而额外的CT成像会导致成本和辐射的增加。通过QSM定义的骨骼信息可以消除MRgFUS和MR/PET中对CT 的需求。

本文对预处理TFI的初步实现方式可以从以下方面进行改进。首先，预处理矩阵 P的当前选择是经验性的，而其性能可以通过对患者提供有针对性的模型或者参数来改善。进一步的工作可能集中于利用磁化率分布的物理或统计模型，更好地近似协方差矩阵Γ来确定P。QSM可以通过诸如SSQSM等方法进行初步的快速计算，用来建模该分布和构造预处理矩阵。其次，我们构建预处理矩阵时需要R2^*图来提取具有强磁化率的脑内区域，如ICH。如果R2^*信息不可用，特别是仅获取单一回波图像时，预处理的TFI可能无法有效抑制基于ICH的阴影伪影。在未来的工作中，我们将着重于将其他对比度(幅值，T1加权或T2加权图像)结合到预处理矩阵的构建过程中。第三，对于全头部QSM，颅骨具有非常短的T2^*，使得多回波GRE序列未能捕捉足够的颅骨MR信号用于场估计。在未来的工作中，使用超短回波时间(UTE)脉冲序列来提供颅骨中的相位信息可能会实现更精确的磁化率图估计。来自UTE序列的幅度信号也可以用于区分颅骨和充满空气的鼻窦，进而通过为这些区域分配不同的权重 P_B来实现更有效的预处理。最后，所提出的TFI方法使用线性信号模型，但可以类似地扩展到非线性方程式，以绕过相位解缠步骤并改进噪声建模。但是，非线性和非凸的目标函数会具有多个局部最小值。考虑到ROI之外的期望磁化率相当大，如空气，本工作中使用的全零初始猜测可能会导致求解器收敛到不正确的局部最小值。未来的工作将集中在构建对这个问题有效的初始猜测。目前的线性解可能会提供这样的初步猜测。

总之，我们描述了用于QSM重建的预处理整体场反演(TFI)算法。它消除了背景场去除和局部场反演的连续进行的需要，并且在磁化率图边界处实现更高的精度。我们还引入了一个预处理矩阵来提高收敛性。实验表明，预处理TFI能够抑制脑出血QSM中的条纹伪影。TFI还能够描绘整个头部的磁化率分布，包括脑组织，颅骨，鼻窦和皮下脂肪。

2.各向异性权重结构先验：

在本节中，我们提出通过考虑结构边缘取向来构造结构先验的改进。

2.1简介：

通过对磁场-磁化率反演问题的数学分析得知，QSM中的条纹伪影起源于场的偶极不相容部分。这种不兼容的部分包括噪声，离散化引起的误差，各向异性源或化学位移。此外，无条纹的解可以只从偶极兼容的部分获得。然而，在实践中噪声是无所不在的，离散化对于数值计算是不可或缺的，并且在单角度采集中不能可靠地测量各向异性源。因此，图像空间中的正则化在计算最小条纹解中起着至关重要的作用。

其中，形态学相似的偶极反演(MEDI)是一种通过在候选磁敏度图和与磁敏度图相关联的幅度图之间施加空间一致性(边缘信息)来增强组织结构的方法。具体而言，结构先验作为各向同性矩阵并入全变差(TV)正则化中。这旨在：1)保存幅度图像与候选磁化率图之间的共同边缘-组织结构；2)惩罚其他地方的强梯度– 条纹伪影。因此，其性能在很大程度上取决于如何设计这样的结构先验(矩阵)。本文重点讨论这种结构先验的设计，它考虑了以前在各向同性加权中缺失的方向信息。

本文的其余部分组织如下：在理论中，我们考虑结构先验的各向异性权重(AW)，其中考虑到方向一致性。在方法和结果中，我们证明了AW在两个体内脑部MRI数据集中的有效性，其中COSMOS或磁敏度张量χ₃₃的分量可用作参考。在讨论和结论中，我们对我们的结果进行了一些讨论。MATLAB源代码将在 http://weill.cornell.edu/mri/上提供。

2.2理论

在空间连续设置中，在感兴趣区域Ω内的正则化偶极子反转问题被写为以下变分问题：

其中从测量的局部场b：来估计磁化率分布χ：(L¹-可积)。有界变差(BV)的容许解空间BV(Ω；R)允许存在跳跃(尖锐边缘)。数据保真项中的映射w：L¹→L¹补偿了非均匀相位噪声(SNR加权)。这里，d和*分别是单位偶极子核和卷积运算。正则项Reg(χ)是这样设计的，它可以有效地惩罚源自场b的偶极不相容部分的条纹伪影。例如，在MEDI中，正则项表述如下：

其中δ_x(r)，δ_y(r)和δ_z(r)分别是在r＝(x，y，z)处的x，y和z方向的边缘指示函数，即如果在r₀＝(x₀，y₀，z₀)处存在沿着x轴的边缘，则δ_x＝(x₀，y₀，z₀)＝1。否则 δ_x(r₀)＝0。同样，定义了δ_y(r)和δ_z(r)。由于这些边缘指示函数是在假定两幅图像中的边缘共存的情况下根据与未知磁敏度图相关联的幅度图像计算出的，所以正则化矩阵Reg(χ)仅惩罚不期望存在组织结构(边缘)的区域。请注意，λ>0 是一个平衡数据保真度和正则项的调整参数。

为了考虑边缘方位信息，各向异性权重可以定义如下：

其中ξ(r)是边缘指标δ：乘以幅值图像Mag(r)在r处的梯度向量，如下所示：

请注意，∨表示布尔或运算。矩阵P_⊥ξ是秩为2的正交投影算子，它消除了ξ方向上的分量。换句话说，我们修改正则化项，使得它惩罚与幅度图像梯度正交的磁化率分布梯度的分量，如图9所示。因此，我们促进幅度和磁化率边缘之间的平行取向-而不仅仅是并存，这增强了两幅图像之间的形态学一致性。

数据保真度项不仅是凸的，而且正则化项Reg(χ)和Reg^⊥(χ)在χ中都是凸的。因此，可以使用一些成熟的算法用于这种凸优化问题，如高斯-牛顿共轭梯度法，分裂布雷格曼方法和原始对偶算法，以有效地将能量最小化[12]。在本文中，我们使用高斯-牛顿共轭梯度法与索伯列夫空间假设，它解决了与能量函数第一变分相关的欧拉-拉格朗日方程[12]。Reg(χ)的第一变分如下所示：

其中梯度算子应该在BV空间的分布意义上理解。请注意，上述第一变分[16]在BV中没有明确定义，但是在索伯列夫空间中，这是高斯-牛顿共轭梯度法背后的假设。 Reg^⊥(χ)的第一变分推导如下：

同样，这个第一变分是在索伯列夫空间中定义的。

2.3素材和方法

这里我们描述实验细节来评估TV正则化中各向同性权重和各向异性权重。

两名健康受试者的两项体内脑部MRI数据被用于实验验证使用COSMOS或χ₃₃作为参考的各向异性加权的影响。这些数据集来自康奈尔MRI研究组 (http://weill.cornell.edu/mri/pages/qsm.html)以及QSM 2016重建挑战 (http://www.neuroimaging.at/qsm2016/pages/QSM-challenge.php)。对于成像参数和相位预处理，我们请读者浏览网站和参考文献中的内容。

数值调节设为∈＝10^-6；内部迭代的最大数(嵌套CG循环)被设置为100；外迭代的最大数设置为100；嵌套CG循环的容限设置为0.01；并且相对更新范数 ||χ^k+1-χ^k||/||χ^k||被设置为0.01，其中χ^k是第k次外迭代处的解。请注意，停止标准比默认设置更加严格，以确保收敛-通常外部迭代的最大数为10，相对更新范数为0.1。该算法在Matlab R2016a上运行，CPU为64GB内存的英特尔酷睿i7 3.30GHz。

在各向同性加权中，边缘指示函数δ_x,δ_y,和δ_z由以下过程确定：1)计算和对于r∈Ω；2) 将边缘指示函数定义为参数t的特征函数如下：和3)发现这样的t:

然后，所得到的二值指示函数δ_x(r),δ_y(r),和δ_z(r)用于[13]和[16]。

在各向异性加权中，我们需要确定ξ(r)×(ξ₁(r)，ξ₂(r)，ξ₃(r))^T。这归结为如[15]中定义的计算δ：一旦δ_x,δ_y,和δ_z通过上述各向同性加权的程序计算，我们将δ确定为δ_x(r)∨δ_y(r)∨δ_z(r)对于r∈Ω.。使用康奈尔MRI研究小组的数据，这两种加权方法可视化为图10。

至于λ/2的值，1000和235分别根据偏差原理从康奈尔MRI研究实验室数据和QSM2016重建挑战数据经验性地选择。然后，对于康奈尔MRI研究组数据，在范围为从700到1350(700,750，...，1350)，对于QSM 2016重建挑战数据，在范围从164.5到305.5(164.5,176.25，...，305.5)实施各向异性和各向同性加权的MEDI算法。

以下评估指标用于QSM准确性评估：RMSE(均方根误差)，HFEN(高频误差范数) 和SSIM(结构相似性指数)。我们使用在 http://www.neuroimaging.at/qsm2016/pages/qsm-challenge.php获得的RMSE和 HFEN MATLAB代码，而对于SSIM，我们在图像处理工具箱中使用了MATLAB函数(SSIM) 和默认设置。另外，我们使用了我们所说的对齐度来衡量两幅图像的对齐程度，它定义如下：令u和v是标量值图像(例如中的函数)；并让n_u，n_v分别表示 u和v表面法线的矢量场。假设是r∈Ω处两个矢量场之间的夹角。然后，两个矢量场的对准度(DoA)被定义为

其中|·|表示中的Lebesgue测度，即体积。请注意，数量DoA(u，v)的下界为0(当两个矢量场几乎处处垂直时)，并且其上界为1(当两个矢量场平行时- 无论是0弧度还是π弧度)。因此，数量DoA(u，v)测量两个图像u和v在[0,1] 范围内的对齐程度(较高表示较好的对齐)。

为了根据不同大脑区域的平均体素值来评估QSM的准确性，19位感兴趣的圆形区域(每个区域的测量值为14.5mm3)由具有超过20年经验的神经放射学家在选定的脑切片上绘制在纹状体(尾状核，壳核，苍白球)，丘脑，皮质灰质，胼胝体压部，额叶白质，枕叶白质和侧脑室脑脊液(CSF)。在MATLAB中计算每个ROI的平均体素值，并相对于CSF平均值记录。对于正规化参数λ的各种值，记录用于MEDI 导出的ROI值与COSMOS或χ33导出的ROI值之间的线性回归的斜率和确定系数R2。这些拟合参数在MEDI的各向同性或各向异性加权之间进行比较。

2.4结果

对于康奈尔MRI研究小组的数据和λ/2(700至1350)的值，各向异性和各向同性加权的CPU运行时间分别为大约38分钟和44分钟。对于QSM 2016重建挑战数据和λ/2(164.5至305.5)的值，各向异性和各向同性加权的CPU运行时间大约为10分钟。

图11显示了康奈尔MRI研究实验室数据(λ/2＝1000)和QSM 2016重建挑战数据(λ/2＝235)在轴向，矢状和冠状视图中各向异性和各向同性加权的计算磁化率图。在磁化率图上，我们看到各向同性加权中的条纹伪影比各向异性加权更多。这些伪影在轴向视图中表现为块状斑点(见图11(a)和(b)左侧的第一列)，矢状和冠状视图中的表现为条纹(参见图11左侧的第二和第三列(a)和(b))。箭头置于这些伪影显眼的地方；请注意，在矢状面和冠状位图像上显示条纹。这些伪影在相对于参考图像的各向异性加权(COSMOS或χ_33)中被显著抑制。

图12显示了RMSE(越低越好)，HFEN(越低越好)，SSIM(越高越好)和DoA(越高越好)，其中所有测量均相对于正则化参数λ/2绘制。从上到下，每行对应于关于COSMOS的康奈尔MRI研究组数据，关于COSMOS和χ₃₃的QSM 2016重建挑战数据。从图12中可以清楚地看出，通过各向异性加权实现的QSM精度的改善在正则化参数范围内的两个数据集上是一致的。在整个正则化参数范围内，康奈尔MRI研究组数据的绝对改善RMSE平均为9.4％，HFEN平均为8.4％，SSIM平均为0.4％，DoA平均为0.8％。在整个正则化参数范围内，QSM 2016重建挑战数据的绝对改善关于 COSMOS的RMSE平均为8.7％，HFEN平均为5.2％，SSIM平均为0.4％，DoA平均为 0.5％；关于χ₃₃的RMSE平均为8.7％，HFEN平均为5.5％，SSIM为0.4％，DoA为 0.3％。

通过基于各向同性和各向异性方法和参考标准(COSMOS或χ₃₃)之间的ROI内的体素值的线性回归，以评估不同方法的QSM精确度。图13显示了康奈尔MRI研究组数据(关于COSMOS-从顶部开始的第一行)和QSM 2016重建挑战数据(关于COSMOS 和χ₃₃，分别从顶部起第二和第三行)中的斜率和R²关于λ的变化曲线。在正则化参数范围内，各向异性权重始终比各向同性权重有更高的斜率:平均而言，对于康奈尔MRI研究组数据，各向异性和各向同性权重的斜率分别为0.949和0.947。对于 QSM 2016重建挑战数据，各向异性和各向同性权重下χ₃₃的斜率分别为0.927和 0.919。就COSMOS而言，各向异性和各向同性加权的斜率分别为0.875和0.866。

2.5结论

我们已经将各向异性权重视为QSM逆问题中的新结构先验。与各向同性权重相反，各向异性加权强化候选磁化率图与与磁化率图相关联的幅度图之间的方位一致性。实验显示QSM精度在多种正则化参数下关于量化测量和图像质量均有改善。

在本文中考虑的各向异性权重的概念在关于图像处理中的各向异性扩散的开创性专著中并且在矢量值图像处理中被广泛研究，其矢量值图像处理的主要焦点是检索不同通道的共同边缘(例如，RGB)和抑制随机噪声。最近，以“定向总变分” 为命名的研究开展了一项类似的图像处理工作，另外一项类似的工作已用于快速多磁共振成像重建，其中T₁加权图像中的结构信息被用于重建T₂加权图像(稀疏k空间样本)。本文所提出的方法与这些工作类似；但是这一方法在定量敏感性分析中是全新的。

在图11中，我们观察到各向同性加权比各向异性加权更多的视觉伪影，这明显与参考图像(COSMOS和χ₃₃)不一致。一个直观的解释如下：考虑两个二维矢和r＝r₀，其中θ(r₀)是r＝r₀. 与x轴的角度。以下各项在各向异性权重中受到惩罚：

因此，当θ(r₀)＝π/4和5π/4时，没有惩罚；当θ(r₀)＝3π/4和7π/4时，惩罚最大化。另一方面，在各向同性加权中，以下是正则化：

其中δ_x(r₀)和δ_y(r₀)分别是x和y方向的边缘指示函数。因此，我们有|sinθ(r₀)|,|cosθ(r₀)|,和0分别惩罚x，y和x，y两个方向边缘。显然，这不考虑方向一致性-例如如果检测到x和y方向边缘，则允许的每个可能的方向，即θ(r₀)没有惩罚。这种在各向异性加权中的“自适应惩罚”似乎抑制了源于定向不一致性的这种视觉伪影。

康奈尔大学MRI研究组数据和QSM 2016重建挑战数据的各向异性和各向同性权重之间的CPU时间差分别约为6分钟和几秒。各向异性权重在康奈尔MRI研究小组数据上运算更快，在QSM 2016重建挑战数据上和各向同性权重一样快。各向异性加权涉及3×3全矩阵乘法，而不是3×3对角线矩阵乘法以进行各向同性加权，导致在高斯-牛顿共轭梯度法的每次外迭代时额外增加几秒钟。然而，通常情况下，各向异性权重以较少的外迭代次数收敛得更快，这解释了各向异性权重与各向同性权重一样快的事实。

设计一个考虑到QSM背后的生物物理模型的结构先验是基于正则化的反演方法取得成功的基石。正如本文所见，各向异性加权方法的性能优于各向同性加权，其实现非常简单。我们相信所提出的方法将是朝着最佳场分解方法迈出的重要一步；尽管目前未知。

总之，将各向异性权重纳入MEDI的正则化项可以通过考虑各向同性权重中缺少的方向一致性来提高QSM的准确性。

3.原始对偶算法

在本节中，我们将介绍QSM问题的原始对偶算法，它可以提高计算L1范数项的导数的准确性和速度。

3.1.简介

贝叶斯最大概率估计已被越来越多地用于磁共振成像中，以重建来自不完整且含噪声数据的图像，它使用似然估计对噪声建模，同时使用先验概率弥补欠采样，其在图像重建中的应用常常需要解决最小化问题。高斯-牛顿共轭梯度(GNCG)算法可以解决此类最小化问题。定量磁化率成像(QSM)就是这种情况，其中，一个结构先验概率被用于确定含最小条纹噪声(形态学使能偶极子反演，MEDI)的解。结构先验强加了已知组织图像和目标组织磁化率图之间边缘(梯度)的空间一致性。

本文的目的是研究L1范数项的具体实现对QSM质量的影响。首先，典型的GNCG 实现采用小幅修改先验项的方法平滑L1范数，从而求解近似的能量函数。虽然这种调节的影响已经在图像去噪方面进行了研究，但它对QSM的影响尚不明确。我们将比较这种算法和原始对偶算法。原始对偶算法已用于计算贝叶斯磁共振成像和QSM 中的全变分。其次，GNCG方法中，中心差分算法已经被用来计算梯度，这可能会导致棋盘伪影。我们将与另一种方法，即前向差分方案进行比较。因此，本文将使用具有正向梯度实现的Chambolle-Pock算法来评估QSM的原始对偶算法，并将其与具有中心差分算法实现的GNCG方法进行比较。

3.2.理论

在这一部分，我们将概述QSM能量最小化问题的空间连续公式，以描述各向异性全变分的对偶问题，同时不使用GNCG中的近似导数。原始对偶算法的细节可以在 MEDI工具箱中找到(http://weill.cornell.edu/mri/pages/qsm.html)。

A MEDI的空间连续公式

目前，QSM已被定义为用于估计磁化率分布的离散能量最小化问题。它的连续形式可以用下面的变分模型来表示：

在这里，磁化率函数u：(L¹-可积的)通过测量到的磁场强度来估计。 b：(Ω是李普希兹域).含边界变分的解空间允许突变的存在(边界突变).在离散形式下，目标函数含两个能量项,数据约束项和正则项. 数据约束项包含一个线性算子T，由Tu(r)＝(d*u)(r)，r∈Ω定义,其中d和* 是单位偶极核和卷积算子。我们用L¹代表绝对可积函数(L¹-可积的)空间。映射W：L¹→L¹弥补了不均匀的相位噪声(SNR加权).正则项是各向异性加权的全变分空间正则，同时使用从对应磁场图提取的边缘图 M：梯度算子D可以从分布的角度理解.参数λ＞0调节数据约束项和正则项,可以通过L曲线方法的差异性原理得到.基于MEDI的空间连续性方程,通过BV去模糊理论，我们可以确定它存在一个最小值解。在这个公式中，磁化率u是r的函数，我们需要将其最小化。随后，在QSM文献中常用的磁化率χ通过离散u给出。我们的目标是计算这个最小值解，我们从GNCG算法开始。

高斯-牛顿对偶梯度(GNCG)

GNCG方法假定u是可微分的(理论上我们假定u存在于索伯列夫空间 W^1，1),这允许我们推导出含能量项的欧拉-拉格朗日方程:

其中(WT)^*是WT的共轭.然而当时欧拉-拉格朗日方程会退化，所以出于计算目的，我们需要增加数值条件:在以及∈＞0但很小时被替换。这会导致这个条件与L¹- 可积函数中对TV的定义不同，其中图像边缘由离散点集描述。因此，这个算法解决的是包含条件参数∈的稍有不同的问题。

原始-对偶公式

与含离散点函数的分布导数类似,我们可以在L¹空间对可积函数得出一个统一的对各向异性加权TV的定义：

其中是有紧支撑的平方可积函数的集合。 ||ξ(r)||_∞＝max{|ξ^x(r)|，|ξ^y(r)|，|ξ^z(r)|}并且M^*是M的共轭.ξ是包含逐点最大范数约束的对偶变量。向量的散度记为div,集合的上确界记为sup.注意在上述对偶形式下u可以是不连续的.因为是自共轭的，而且对于 H：L²→(-∞，+∞]且α＞0的条件,(αH)^*＝αH^*(·/α)，数据项可以用如下方式对偶:

其中是在平方可积函数空间下定义的内积，即,因此，最小化问题[18]可以被[20]和[21]重新定义，产生了所谓的原始-对偶公式:

其中，集中的示性函数ι_P(ξ)有如下定义：

其中凸集

这是一个有凸约束的鞍点问题.在存在鞍点满足解的假设下[22]，可以应用可证明的收敛算法来计算这样一个最优点。其中，我们利用Chambolle-Pock提出的快速原始-对偶(PD)算法。在下文中，我们提供了[22]的离散设置和PD算法的实现细节来计算最佳磁化率图。请注意，GNCG现在可以被看作是原始的方法。

离散设置

要计算鞍点问题的最优解[22],首先我们要进行离散化设置。我们考虑一个三维笛卡尔坐标系中的栅格，大小是N_x×N_y×N_z:

{(ih_x，jh_y，kh_z)：1≤i≤N_x，1≤i≤N_y，1≤i≤N_z}，

其中h_x,h_y,h_z表示体素大小,(i，j，k)表示位置.接下来我们可以定义标量积

对于u∈H,离散梯度有如下定义

其中

离散向量是在有限维实数空间K＝H×H×H中的向量，其中标量积如下定义：

注意，使用纽曼边界条件的前向差分方法被用于计算离散梯度。现在我们定义一个线性映射M：K→K使得

对所有i，j，和k成立。随后我们可以得到[20]的离散对偶问题：

相应的，和M^*：K→K是和M的共轭。即, L^∞-范数||·||_∞由如下定义:

从和的离散设定,如下的共轭算子的离散形式容易得到：

注意离散梯度的共轭算子是负梯度，即它由含狄利克雷边界条件的后向梯度算法进行离散化。

原始-对偶算法

相应地，函数空间和L²被K和H替代,同时，凸集[23]采用如下公式进行离散化:

随后，[22]的离散形式可以容易地得到。这种离散原始对偶公问题可以由一个可证明收敛鞍点的算法有效地求解。在这里，我们采用一阶原始-对偶算法，它的超松驰迭代收敛性是由Chambolle-Pock证明的。该算法基本上被认为是相对于原始变量 u的交替梯度下降方法，和耦合<W(Tu-b)，v>_H的逼近算子，以及对于对偶变量ξ 和ν的梯度下降方法，以及耦合ι_P(ξ)和的逼近算子。

我们使用<W(Tu-b)，v>_H定义G，随后放缩函数的逼近算子τG,τ＞0由下式给出：

其中和表示快速傅里叶变换(FFT)和反傅里叶变换IFFT。是的复共轭，算子⊙表示逐点乘积。注意对这种算子在每次循环中，只需要进行一次FFT和一次IFFT。被提前计算。另一方面，放缩函数σι_P(ξ)和的逼近算子如下给出：

对所有i，j，和k成立。最后，原始-对偶算法表示如下:

对n＝0，...，N-1。在每个循环中，要更新v，需要一次FFT和一次IFFT。在初始化过程中，我们选择(u⁰，ξ⁰)∈H×K和set原始和对偶参数τ和σ 满足来保证算法收敛性。由于和我们有:

3.3.实验和方法

这里我们描述实现和评估使用前向梯度的原始-对偶QSM方法的实验细节。

数据集介绍

我们使用大脑仿真，钆(Gd)模型，4个健康受试者的体内脑部MRI数据和离体大脑组织用于实验验证，在仿真和钆仿体实验中，使用真实值作为对照，在在体成像和脑组织实验中使用COSMOS作为对照。脑部仿真和一组在体脑部MRI数据来自 http:// weill.cornell.edu/mri/pages/qsm.html。剩下的在体数据是从中抽取。对于相位预处理，在每个体素中进行相位的时间解缠，然后在TE上的每个体素上对时间解缠的相位进行加权最小二乘拟合。然后，我们应用幅度图导航空间解缠算法。然后通过将投影到偶极子场(PDF)算法提取原磁场。请注意，模拟大脑的白质是分段光滑的，大脑的其他部分(灰质等)是分段恒定的。

将来自三个死亡脑组织(对照组)和三个死后帕金森病大脑组织(美国耶鲁-纽黑文医院)的切下的脑组织固定并包埋在1％琼脂糖凝胶中。使用标准8通道头部线圈使用3T临床MRI系统(美国GE)扫描脑块。使用以下参数获取多回波3D梯度回波(GRE)数据：0.4mm各向同性分辨率，TE＝4.3ms，ΔTE＝4.8ms，#TE＝11 ，TR＝74.2，BW＝62.5kHz，FA＝20°。我们采集了三个共面的方向(Δφ＝60°)。对于相位处理方面，我们使用了时间拟合方法和空间相位解缠算法。

前向差分和中心差分

用h_x，h_y，和h_z表示离散宽度，或体素大小,(i，j，k)表示离散体素 (ih_x，jh_y，kh_z)在图象域中的位置；随后，基于中心差分的用于梯度算子的关于u的各向异性TV算法有如下定义：

另一个基于前向差分的用于梯度算子的关于u的各向异性TV算法有如下定义：

随后，各向异性TV的离散形式可以轻易得到。GNCG使用中心差分(GNCG-C)和前向差分(GNCG-F)实现。PD使用前向差分实现。

GNCG和PD的参数设置

对GNCG,数值条件被设置为∈＝10^-6,内循环最大次数(嵌套CG循环)是 100,最大外循环次数是20,嵌套CG循环和相关的范数||u^k+1-u^k||/||u^k||,其中u^kt是第k次外循环的解的误差允许范围，都被设置为1％.对PD算法,相关范数被设置为0.0001％，同时我们设置来保证收敛,其中对离体大脑切片对其余数据集这些算法在Matlab R2015b上编写，运行于Intel Core i7 3.30GHz CPU及64GB内存。此外,PD使用支持GPU 计算的MATLAB(未优化)编写，在NVIDIA GPUs(NVIDIAGeForce GTX 750Ti和 TITAN X)运行。

收敛情况

在仿真大脑实验中GNCG-C,GNCG-F,和PD使用了多种不同的参数(对GNCG是∈，对PD是τ和σ)。这些算法的停止标准被设置为比前面描述的更严格，以方便调查每个算法的收敛性。最后，在优化过程中第k次循环，我们会计算相关范数 ||u^k-u^*||₂/||u^*||₂,其中u^*是真实值。

定量分析

在体和离体COSMOS图像被读入ITK-SNAP图像处理系统，随后感兴趣区域 (ROIs)被标记出来。在体ROIs皮层，深部灰质结构(红核，黑质，尾状核，壳核，苍白球)和白质区域(胼胝体压部，额叶白质，枕叶白质)。离体数据ROIs包括皮层，壳核，苍白球，红核，黑质和丘脑底核。

通过拟合回归线，我们进行了线性回归,X＝k₁X_COSMOS+b,其中X是ROI 区域对GNCG-F或PD的平均值，而X_COSMOS是对应的COSMOS值集合。斜率k₁, 截矩b,和决定系数R²在每次回归中都被记录。为了进一步评估COSMOS和两种重建方法之间的一致性，我们绘制了Bl和-Altman图，以确定平均差异(偏差)和 95％的一致性上限。

3.4.结果

如下是我们的结果，依次是收敛分析，钆造影实验，体外大脑成像和体内组织成像。

收敛性分析

图14展示了算法的相对误差,数值条件取10^-8(高精度)的PD比GNCG-F和 GNCG-C收敛更快，和数值条件取10^-6的GNCG-F收敛一样快，但后者精度低.。3000 次循环后,对GNCG-F(∈＝10^-6)和GNCG-F(∈＝10^-8)相对误差分别为8.93％和 8.84％。然而GNCG-F(∈＝10^-8)没有达到PD的精度:3000次循环后,对PD(L³) 和PD((L³+L4⁾/2)相对误差分别为8.62％和8.62％。正如上文所说,PD算法在 1000次循环时同样有最高的精度,对GNCG-C(ε＝10^-6)，GNCG-C(ε＝10^-8), GNCG-F(ε＝10^-6),GNCG-F(ε＝10^-8)，PD(L³)和PD((L³+L⁴)/2),相对误差分别为10.92％,10.96％,8.96％,9.45％,8.67％,和8.51％，其中

对GNCG和PD算法每次循环需要的浮点运算数几乎相同。这里，迭代是指用于 GNCG的嵌套CG循环的单次循环(#操作＝2个FFT+2个iFFT+3个点乘积+2个 SAXPY+2个GAXPY)和用于PD的算法的一个循环(#操作＝2个FFT+2iFFTs+ 4SAXPYs+2GAXPYs)。请注意，对于GNCG，一旦嵌套CG循环终止，将从外部循环添加附加操作(1SAXPY和1GAXPY)。对于GNCG，每个外环的最大内部迭代次数为100次，实际内部迭代次数从外循环到外循环和从数据集到数据集(50～100)不等。

对GNCG-C(ε＝10^-6),GNCG-C(ε＝10^-8),GNCG-F(ε＝10^-6),GNCG-F (ε＝10^-8),PD(L³),和PD((L³+L⁴)/2)算法,单次循环CPU运行时间大约为0.79 秒,0.80秒,0.75秒,0.75秒,0.75秒,和0.75秒。

图15展示了在GNCG-C中出现的棋盘伪影在GNCG-F中被去除,同时GNCG-F 的误差相比于PD进一步降低。

钆仿体

GNCG-F的准确性如下:对λ/2＝150,线性回归斜率＝1.06,R²＝0.79；对λ/2＝300,斜率＝1.06,R²＝0.92.PD的准确性如下:对λ/2＝150,斜率＝1.04,R²＝0.87；对λ/2＝300,斜率＝1.04,R²＝0.81。对GNCG-F和PD算法，单次循环CPU运行时间(钆仿体大小＝130×130×116)大约0.36秒,对PD算法， GPU(GTX 750Ti)运行时间大约0.03秒。对GNCG-F总循环次数为958(10次外循环).对PD总循环次数为800。

离体脑组织成像

图16和17展示了GNCG-F和PD在离体脑组织成像上的结果和定量分析。对 λ/2＝600,GNCG-F和PD斜率是0.79(R²＝0.91)和0.81(R²＝0.91)。对 λ/2＝1000,GNCG-F和PD斜率是0.81(R²＝0.91)和0.82(R²＝0.91)。对大小为512×512×136的脑组织,GNCG-F和PD算法的单循环CPU运行时间是大约 4.15秒和4.05秒。PD算法的GPU(TITAN-X)单循环运行时间是大约0.22秒。对 GNCG-F总循环次数为2000(10次外循环)。对PD总循环次数为2100。

在体脑部MRI

图18和19展示了GNCG-F和PD在在体脑MRI成像上的结果和定量分析。对 λ/2＝600,GNCG-F和PD斜率是0.99(R²＝0.87)和0.99(R²＝0.87)。对 λ/2＝1000,GNCG-F和PD斜率是0.99(R²＝0.87)和0.98(R²＝0.88)。对图 18中大小为230×230×46的物体,GNCG-F和PD算法的单循环CPU运行时间均是大约0.29秒。PD算法的GPU(GTX 750Ti)单循环运行时间是大约0.05秒。对GNCG-F 总循环次数为1304(14次外循环),对PD总循环次数为1000。

3.5.结论

在本文中，我们展示了QSM中正则化项的各种实现对收敛速度和准确性的影响。与使用原始对偶(PD)算法的QSM鞍点问题的连续形式相比，使用GNCG求解方法所必需的数值调节过程导致更慢的收敛。正则化项中梯度的中心差分近似的使用可能大大减少在使用正向差近似时的棋盘伪影。

Chambolle-Pock(原始-对偶)算法已应用于QSM背景下的总体广义变分(TGV)，而分裂布雷格曼方法可替代地用于快速QSM。我们的贡献包括以下几个方面：首先，与使用简化数据保真度项(最小二乘法)不同，我们考虑涉及SNR权重的MEDI问题 (加权最小二乘法)，这对于减少QSM中的噪声误差是必要的，但会显著增加计算成本，因为每个CG迭代和PD迭代增加了两个3D傅立叶变换。其次，我们使用加权各向异性TV而不是TGV来强化已知组织结构和靶向磁化率图(形态学使能偶极子反演， MEDI)之间的结构一致性。这种物理MEDI已经与其他现有的临床数据方法进行了定量比较，与未加权的TV或TGV定义的通用稀疏性相比具有优势。第三，我们提供了 QSM变分方法的观点。

注意到[18]中的QSM问题和图像去模糊问题有如下的不同:首先，与图像去模糊中的传统模糊算子相反，QSM中的前向算子(即具有SNR权重的偶极子内核)在远离k空间原点的地方不会衰减到零，而是衰减到由魔术角度(≈54.7°)定义的主磁场展开的圆锥面上。前向算子的这种根本性差异导致QSM中某些奇特的谱特性,这种特性被称为离散皮卡德条件。其次，TV正则项的作用是不同的。通过仅提取场的偶极子兼容部分，我们可以避免条纹伪影。虽然这种提取方法尚不清楚，但文献[20] 中的各向异性加权TV正则能够在保留底层组织结构的同时有效地处理条纹伪影。另一方面，图像去模糊中的TV正则化项在减少附加性高斯白噪声的同时在边缘保持平滑中起作用。

我们在这里的工作说明了标准GNCG中与条件约束有关的收敛性和准确性比PD差(图14)。对该序列，使用GNCG的准确性可以通过将∈从10^-6减小到10^-8进行提高,但代价是显着减慢算法的收敛速度。这是因为由嵌套CG循环解决的线性系统变得病态。更确切地说，CG迭代的收敛速度由线性系统的特征值(频谱)的位置确定,同时条件参数∈的值影响频谱的位置。因此，从计算角度，将∈设置为非常小的值，如～10^-10来提高精度是不切实际的。这使得在非光滑凸优化问题上GNCG 要差于PD。一种解决方法是使用准牛顿(qNewton)连续性方法。即，在每次外循环中，我们通过规则∈＝γ·∈设置∈其中γ是一个递减参数。我们初始化∈＝0.1然后设定γ在0到1之间变动。图20显示了γ取0.01,0.05,0.1,0.3,0.5,和 0.8时的相对误差。

相比于GNCG更快的收敛PD(图14)可以用在每次循环中达到指定精度(例如相对误差和残差)需要的浮点运算次数(FLOPS)来解释。我们对在每次循环中PD和GNCG需要的FLOPS进行计数，其中N＝N_xN_yN_z是体素数目。具体的讲，GAXPY操作(矩阵向量乘法)本身通常是一个计算瓶颈，耗时是O(N²)，其中PD耗时只有GNCG 的3/4。然而，在PD和GNCG中使用的GAXPY操作实际所需的FLOPS数量只有O(N)，因为涉及的矩阵非常稀疏。这些矩阵是从前向/中心差分算法导出的对角线或梯度矩阵(除了第一行和最后一行，每行有两个非零分量)。随后FFT和IFFT操作决定了GNCG和PD算法的运行时间，这些算法是可以并行化的。因此，两个算法的耗时均为O(NlogN)，使用GPU可以进一步加速。在实践中，GNCG和PD的成本远远高于诸如存储器层次和数据之间的通信等操作。然而，由于决定GNCG和PD(数值条件和停止规则)的准确性的其他因素，在计算方面，#FLOPS为我们提供了以下经验法则：在低停止标准下，GNCG算法似乎比PD更快，因为1)GNCG的嵌套CG迭代运行至多100次迭代(仅在前几个外循环中)和2)算法通常在10次外迭代内终止。然而，当我们需要高精度的时候，与PD算法相比，GNCG并不适合，因为嵌套的CG循环需要更多次迭代来达到更小的误差和更多次外部循环。最后，我们发现PD算法的收敛性受步长σ的影响，相应的，虽然当前的参数设置保证了PD算法的全局收敛，但通过自适应调度步长还是有一定的提升空间。请注意，收敛速度的差异直接转化为总计算时间，因为每次迭代的CPU运行时间在GNCG和PD 之间几乎相同。

通过前向差分方法对梯度运算中的棋盘伪影的抑制可以从图15中GNCG-C和GNCG-F之间的区别看出。这可以直观地理解如下.考虑5×5大小的2D棋盘图案： u_1，1＝1，u_1，2＝0，u_1，3＝1，u_1，4＝0，u_1，5＝1，u_2，1＝0，u_2，2＝1，u_2，3＝0，u_2，4＝1，u_2，5＝ 0，u_3，1＝1，u_3，2＝0，u_3，3＝1，u_3，4＝0，u_3，5＝1，u_4，1＝0，u_4,2＝1，u_4,3＝0，u_4，4＝1，u_4，5＝ 0，u_5,1＝1，u_5，2＝0，u_5,3＝1，u_5，4＝0，u_5，5＝1

.我们设置h_x＝h_y＝1，我们由如下定义的中心差分方法计算各向异性的TV项：

前向差分估计有如下的定义：

对十矩阵中嵌套的3×3万格，找们有J_c＝0和J_f＝12.这里没有中心差分方法中的惩罚项，这是因为u在x和y方向上的不断波动确实存在，这种波动被前向差分方法捕捉到了。

GNCG算法中的误差可以用胡伯正则项分析解释.在索伯列夫假设下，在 MEDI中的欧拉-拉格朗日方程的数值条件概念与胡伯正则项有着紧密的联系：

其中α＞0是决定对较小的α的二次项正则和对较大的α的各向异性TV项正则权重的系数。用胡伯-正则项替换掉[18]中的被积函数，我们可以得到当 |MDu|₁≤α时的欧拉-拉格朗日方程：

时间推进算法可被用于解决这个问题；我们的总任务包含求解一个扩散方程，其中扩散系数是1/2α.事实上，这个扩散过程可以在|MDu|₁≤α的区域中出现.现在，阔率一个欧拉-拉格朗日方程，其中分母上的条件系数是∈，例如，在理论部分提到的使用时间推进算法来解决这个问题，我们注意到如果那么对一些∈和α，因此，在图.15中展示的光滑区域GT-GNCG-F(真实值和GNCG-F的差异图)很可能是胡伯-正则项导致的扩散过程的结果。由于胡伯正则项的对偶形式允许我们得到一个相邻算子的手收敛表达,PD算法可以被直接应用。

最近在非光滑凸优化方面的技术进步可以进一步改善QSM精度。基于稀疏诱导准则的能量最小化方法，例如TV-半范数，已经成为病态问题的可行解决方案。然而，使用忠实地反映成像数据的基本物理性质的能量函数依然是最好的选择。否则，计算修改后的能量函数的最优解可能导致我们不需要的解，例如在对数据和先验项的噪声建模中所看到的。在体外和体内数据分析中观察到的MEDI和COSMOS之间的差异也可以由忽略白质磁化率各向异性和错误配准的误差来解释。

总之，我们展示了QSM中正则化项具体实现方法对成像效果的影响，因为它涉及收敛速度和准确性。在标准GNCG算法中数值调节相关的较慢的收敛速度和准确度通过使用不再需要修改QSM中的能量函数的原始-对偶算法得到显着改善。此外，我们已经表明，使用梯度算子的前向差分方法抑制了由中心差分方法引起的棋盘伪影。在从仿体到体内人类大脑的各种类型的数据中，使用前向梯度的原始-对偶算法与当前GNCG相比显示出更好的收敛性和精确性。

4.阴影伪影的抑制：QSM的零参考与多发性硬化症(MS)病灶磁化率的测量

在本节中，我们描绘了一种通过提供均匀的脑脊髓液作为自动零参考的方法，来对QSM进行改进，同时抑制一般的阴影伪影。

4.1.介绍

当前的定量磁化率重建(QSM)方法中磁化率的量化基于某一水平，称为零磁化率参考。对于脑实质QSM，CSF由于其与纯水的化学相似性而被广泛用作零参考。然而，在当前的QSM中，CSF通常存在非均匀性，这增加了参考磁化率的不确定性。文献表明脑脊液流动和白质各向异性可能是导致这种脑脊液磁化率变化的原因。

考虑到在当前贝叶斯QSM(形态学相似的偶极子反演方法，MEDI)中使用组织结构先验知识来处理来自于场到源逆问题中零空间的条纹伪影，我们认为CSF磁化率的不均匀性也属于某一零空间，并建议使用先验知识约束CSF内的均匀性。在本文中，我们提供的初步数据表明，针对QSM的一致和自动化的零参考(称为QSM0)，CSF中磁化率的变化可以被最小化。此外，QSM0能够带来图像质量的改善。

4.2.理论基础

在先前文献(MEDI)中的QSM优化问题中的成本函数包括对应于复数磁共振数据中高斯噪声的数据保真度项和反映关于磁化率分布的结构先验信息的正则化项，后者反映了组织边缘和磁化率图中边缘之间的差异的稀疏性：

这里χ为磁化率图,*为同偶极子d的卷积操作,w为噪声权重,b为组织中经过Lorentz矫正后的磁场，或者简称为组织场或者磁场，为梯度算子，M_G为幅值图像中获得的二值边缘掩码。这里提出的QSM0将额外的正则化项引入到成本函数中，以在感兴趣组织的区域中施加已知的均匀性约束，并求解以下优化问题：

这里相对于MEDI(方程24)多了一项额外的L2-正则化项,而CSF作为具有已知均匀磁化率的感兴趣组织。M_CSF代表的感兴趣的侧脑室区域(ROI)中，CSF平均磁化率为按照均匀CSF磁化率的假设，这一项对脑室中CSF的磁化率变化进行抑制。因此，方程25中的第二个正则化项在感兴趣的区域中施加已知均匀磁化率的约束。为了在足够SNR情况下成像，方程24与25右边第一项可以近似为

我们提出了一项利用脑部ROI掩码M与图来获得侧脑室掩码M_CSF的自动化流程。CSF因其化学成分与水大致相同，具有低R1,R2及R2^*的特性。因而利用CSF 在图像中的生理特性，例如其均匀的R2^*值，通过阈值化的方法从梯度回波的MRI幅值图像或者对应生成的R2^*图像中分割得到。多处区域可以通过施加连通约束进行整合。具体实现过程如下:

(a)阈值化这里假设CSF的低于某一阈值R.

(b)确定脑部质心:,这里N为M重体素数目.

(c)确定脑中心区域:

(d)分析连通性:将分割为各自连通部分M_ci(6-邻域)

然后将最大的3个部分合并:

(e)分析连通性:将分割为各自连通部分M_i(6-邻域)然后将所有与M_cCSF有交集的部分合并:

该过程引入了两个额外参数:λ₂及R,他们会在之后通过数值仿真进行确定。

我们使用高斯牛顿共轭梯度求解优化问题(方程25)。在最后一步中，我们将从磁化率图中减去以实现零参考。

4.3.方法与材料

数值仿真

我们构建了一个数值脑部体模来模拟由各向异性磁化率源产生的CSF的不均匀现象。矩阵大小为256×256×126,体素大小1×1×1.5mm³。原始磁化率张量从一个12-方向人体脑部扫描的数据使用磁化率张量重建(STI)方法进行重建获得。手动确定侧脑室区域(Y.Y.,具有4年经验)，并作为CSF的区域。依据 CSF具有均匀性因而不会导致各向异性的磁化率的假设，CSF的磁化率在各磁化率张量分量中设为0。之后，局部组织场由在B₀＝[0,0,1]情况下通过以下两种方式生成：首先,利用各向异性的前向模型从χ₁₃,χ₂₃,χ₃₃生成组织场f_aniso:

其中F和X为f和χ各自的傅里叶变换；接下来,利用各向同性的前向模型仅从χ₃₃生成组织场f_iso,

最后,在两个场中加入高斯噪声(SNR＝100)。MEDI及QSM0被用于两个场的重建，并针对每种方法计算重建的QSM与χ₃₃之间的方均根误差(RMSE), E_rr＝||χ-χ₃₃||₂。MEDI(方程24)和QSM0(方程25)中的λ₁通过最小化Err来确定。类似地，λ₂及R在QSM0中也通过最小化Err确定。

人体试验:多发性硬化症

8名患有多发性硬化症(MS)的患者的脑部在3T(美国GE)使用多回波GRE进行成像。这项工作的所有研究都经过我们的机构审查委员会的批准。成像序列是：T2w 快速自旋回波序列(FA＝90°,FOV＝24cm,TE＝86ms,TR＝5250ms,层厚度＝3 mm,采集矩阵＝416×256×68)和3D T2^*w抑制梯度回波GRE(FA＝20°,FOV＝24 cm,TE1＝4.3ms,TR＝57ms,#TE＝11,ΔTE＝4.8ms,采集矩阵＝512×512×68, 体素大小＝0.47×0.47×2mm³,带宽＝±62.5kHz,使用带流量补偿的单方向采集方式，并行成像系数R＝2，整体扫描时间～10min)。我们使用非线性场估计及基于图分割的相位解缠绕方法估计整体磁场，然后使用偶极场投影(PDF)进行背景场去除以获得局部场。接下来，MEDI和QSM0在相同的局部场上执行。在两种方法中，我们测量了脑室CSF磁化率的标准偏差。对于病灶磁化率的测量，一名神经放射科医师(Y.Z.，具有4年经验)绘制了T2w图像上每个病灶的ROI。这些图像已经与T2^*w扫描的幅度图像进行了配准。另外我们在所识别病灶的对侧上的具有正常表现白质 (NAWM)中绘制ROI，作为病灶磁化率的参考区域。我们使用线性回归和布兰德-奥特曼图来测量和分析病灶与NAWM之间的磁化率均值差异，从而量化MEDI与QSM0之间的一致性。我们还对包括皮层，深灰质结构(红核，黑质，尾状核，壳核，苍白球) 和白质区域(胼胝体压部，额叶白质，枕叶白质)进行ROI分析。

为了自动测量MS病灶敏感性值，T2加权(T2w)图像上的明亮病灶信号可用于分割病灶。白质可以在T1加权图像上实现自动分割。种子点可以相应地自动放置在T2w 的病灶质量中心。我们可以使用迭代区域生长和阈值(IRGT)算法，在IRGT体积停止增长并追踪整个病灶的那一步迭代来计算病灶体积，随后开始跳跃某个大值并泄漏到周围的结构。

临床研究

我们检查了2011年8月至2015年1月MS患者的磁共振图像，其中每例至少包含两次连续的磁共振扫描，包括T2加权(T2w)，Gd增强T1加权(T1w+Gd)和GRE成像。QSM 是利用GRE数据通过对相位进行解卷积的自动化方式进行重建，其中偶极子核负责连接组织磁化率与从磁共振相位信号中估计的磁场。我们在两次连续的磁共振上比较了病灶，并确定了至少有一个新的T2w病变的患者，即在之前的脑部磁共振中不存在的病灶(距基线磁共振<1年的随访磁共振)。然后将所有新的病灶依据T1w+Gd图像上的增强和非增强性进行了分组。

所有检查均采用3T MR扫描仪(Signa HDxt；美国GE)及八通道头部线圈。每个患者的序列为：(a)T2w快速自旋回波，使用钆之前(b)与之后(c)进行3D反转恢复准备的T1w快速抑制梯度回波，(d)三维T2^*w抑制GRE序列。多回波GRE序列的成像参数如下：TR 57ms；回波数目11；首回波时间4.3ms；回波间隔4.8ms；翻转角 20°；带宽244kHz；视野范围24cm；矩阵大小416×320；层厚度2mm。GRE序列在Gd注射前进行。总成像时间为16分30秒。

QSM是使用形态学相似的偶极子反演(MEDI)从GRE数据重建。通过使用FSL软件(英国牛津的FLIRT软件，)将通过其他模式获得的图像配准至QSM。

在将所有新T2w病灶与先前的磁共振进行比较后，三位神经放射学医生(JC，AG 和GC，分别具有18，9和8年的经验)使用T1w+Gd图像将这些病灶分类为增强或非增强的病灶。他们还将QSM上的所有病灶分类为相对于邻近白质的高信号和等信号病灶。如对病灶分类出现差异则通过多数结果进行解决。

一名神经放射学医生(Y.Z.，具有4年经验)在对Gd增强分类不知情的情况下，绘制了T2w图像上每个局部病灶的区域。在T1w和T2w图像上没有异常信号的白质区域被鉴定为具有正常表现白质(NAWM)。我们在T2w图像上在已识别病灶的对侧镜面位置处的NAWM区域中，选择具有相似形状和尺寸的感兴趣区域(ROI)作为零参考。然后，使用半自动软件将病灶和NAWM参考的ROI覆盖在QSM图像上以评估病灶磁化率值。我们仔细检查保证ROI内不含静脉或伪影。

使用相对磁化率作为区分增强性与非增强性病灶的手段，接受者操作特征(ROC)曲线被用于确定敏感性，特异性和最佳截断敏感性值(单位十亿分之一(ppb))，引导模型估计被用于AUC，而95％置信区间用来评估方差。刀切交叉验证技术用于评估模型的预测性能。使用广义估计方程(GEE)来预测三种病变类型的QSM值：结节，壳和非增强。该模型使用高斯分布和可交换的相关结构假设来解释每个患者的多个病灶。利用每位患者重复测量结果，GEE分析还用于预测增强和非增强病灶的QSM值。所有统计分析均使Windows版SPSS(版本16.0；SPSS美国芝加哥)进行。P<0.05被认为有统计学意义。我们还计算了增强病变患者识别的准确性。

4.4.实验结果

数值仿真

图21显示了数值仿真中QSM及相对χ₃₃的对应误差图。使用f_iso作为输入, MEDI和QSM0获得相对于χ₃₃较小的RMSEs(分别为25.4％和25.8％)。当输入含有各向异性的磁化率贡献时(f_aniso),QSM0的RMSE(73.7％)比MEDI(77.3％)低；脑室CSF附近10mm区域的RMSE分别为：QSM0 65.1％，MEDI 74.3％。这反映了对应于在丘脑下核附近的低信号伪影的减少(图21f和21g)。我们使用λ₁＝0.001, λ₂＝0.1及R＝5s^-1。分割的脑室区域内磁化率的标准偏差分别为：MEDI 36.5ppb， QSM0 11.8ppb。

图22a展示了使用R＝5s^-1时的CSF自动分割过程。图22b与22c显示了选择不同R时的RMSE,ROI掩码与QSM。最优R(R＝5)通过最小化的RMSE值 (73.7％)确定。当R过大时(R＝12),脑室外部区域被包含至分割结果中，而当R 过小时(R＝1)，算法无法完整捕捉脑室结构，特别是小脑前方的第4脑室。当R从5 增大至8的过程中，RMSE始终保持在75％以下，且QSM质量没有显著变化。

人体试验：多发性硬化症

图23显示了对于两名MS患者使用MEDI和QSM0重建的脑QSM。使用λ₁＝0.001, λ₂＝0.1及R＝5s^-1。在图23中用箭头表示的脑中心部分，低信号阴影伪影被抑制。在QSM0中，脑室内CSF的标准偏差比MEDI低5倍。在磁化率测量中，对于来自8名患者的所有23个病灶，QSM0显示出与MEDI(k＝1.11,R＝0.93)的强相关性(图24a)。布兰德-奥特曼分析表明，QSM0相对于MEDI具有1.3ppb偏差与[-8.7,11] ppb的一致区间。

临床研究

从符合条件的482名MS患者中，我们确定了54名至少有一处新T2w病灶的患者；总共有133个新的T2w病灶。(一名患者因GRE图像上的运动伪影而被排除)。其余54名患者(12名男性和40名女性)的平均年龄为34.7±8.1(范围20-52岁)。这些患者的病程为0至18岁(5.71±4.51)，扩展的残疾状态量表评分(EDSS)为0 至6。

在T1w+Gd中，来自33位患者的133个病灶中的86个(64.7％)被鉴定为增强性，而来自25位患者的47个(35.3％)未增强性。三位读者的意见完全一致。对增强性病灶，与邻近的白质相比，QSM中86例中有69例(80.2％)等信号，17 例(19.8％)高信号。根据其对T1w+Gd的增强作用，将增强病灶类型分为69个结节形和17个壳形。17个高信号增强病灶中的13个是壳增强性的。所有47例非增强性病灶在QSM上均为高信号，但其中4例(8.5％)信号仅稍高。图24b中显示了示例图像。

相对于NAWM的病灶的平均磁化率为：非增强性病灶20.26±7.55ppb(平均值 ±SD)，增强性病灶(结节和壳)2.49±6.39ppb。在三种病灶类型之间的病灶磁化率GEE分析中，结节增强(β＝-19.6,95％CI＝-23.5至-15.8，p＝<0.0001)和壳增强病灶(β＝-13.5,95％CI＝-19.0至-8.0，p＝<.0001)与非增强性病灶相比具有显著更低的磁化率值。在增强与非增强病灶之间磁化率GEE分析中，与非增强病灶相比，增强病灶的磁化率明显较低(β＝-17.2,95％CI＝-20.2至-11.2， p＝0.0001)。病灶磁化率模型的交换相关系数为0.12。

图24c显示了由病灶的平均相对磁化率值构建的ROC曲线。根据QSM测量的磁化率值，曲线下交叉验证的面积(AUC)为0.9530(95％CI：0.9201至0.9859)，而经过引导后AUC为0.9594(95％CI：0.9305至0.9884)。11.2ppb作为相对磁化率截断值，用以区分增强性和非增强性病灶，其敏感性和特异性分别为88.4％和 91.5％。

4.5.结论

本文提出QSM0方法，基于目前单方向MEDI QSM，使用自动化的脑室分割及加入针对脑室中脑脊髓液磁化率变化的L2正则化项。我们的数据表明，QSM0方法能够减少伪影，包括CSF磁化率的不均匀性及低信号伪影。这些伪影可能是由MEDI中的各向异性的磁化率引起的。同时，QSM0不影响MEDI的组织对比度，如MS病变相对于 NAWM的对比度。因此，QSM0能够实现一致的CSF零参考。

由于测量的组织场对于磁化率中任何均匀的变化是不变的，在QSM中需要零参考来提供磁化率的绝对量化。因为CSF在化学上与纯水相似，所以它是零参考的常用选择，并且常在侧脑室部位通过手动绘制的方法来选取对应ROI。然而，脑室的大小在不同受试间存在很大差异，并且脑脊液在QSM中常常显得不均匀。因此，测量的 CSF值会对ROI选择敏感，从而为零参考带来了不确定性。本文描述的QSM0方法从两个方面解决了这个挑战：首先，它消除了手动选择CSF ROI的需要，因为它能够自动分割脑室并用作零参考区；其次，它在重建期间保证了CSF均匀性，从而降低对ROI的敏感性，正如CSF中较小的标准偏差所示。

本文描述的针对CSF的正则化项来源于对当前QSM中CSF磁化率不均匀性起源的研究，即其原因可能是侧脑室中脑脊髓液的流动或者周围白质束的各向异性特征。这项工作通过模拟具有各向同性与各向异性效应的组织场，在数值模型中研究了白质各向异性在CSF不均匀性中的效果。如图1和图2所示，来自各向异性分量(χ₁₃及χ₂₃)的场在QSM中产生平滑的阴影伪影，如图21f和21g所示。由于QSM0保证了CSF的均一性，在模拟结果的深部脑核中的阴影伪影被显著抑制(图21g)。该模拟结果与体内实验结果能够很好地吻合(图23)。以这种方式，QSM0中的L2正则化不仅减少CSF磁化率变化，还减少了CSF外部的阴影伪影。

针对脑脊液磁化率变化的L2正则化在概念上类似于使用加权全变化作为组织结构先验知识，从而减少条纹伪影。测量的场数据中的粒状噪声误差，通过偶极子内核的零空间进行波传播，表现为条纹伪影。一般的贝叶斯策略是通过使用惩罚与解剖(通常为幅值)图像中的组织结构不一致的结构(边缘)的先验知识项，来抑制伪影。类似地，CSF磁化率不均匀性和阴影伪影可以被认为是由标量偶极场模型中的空间平滑误差引起的伪影。相应的，本QSM0工作中使用了先验项通过惩罚CSF磁化率变化来减少这些伪影。

QSM0可能直接改善对深部灰质核团的描绘，包括丘脑底核，这是深部脑激励手术的重要目标。由于QSM已经用于许多临床应用，因此检查QSM0对CSF外部相对对比度的改变程度变得尤为重要。在测量MS病灶和NAWM之间的相对磁化率差异时(图 24)，MEDI和QSM0之间的强烈一致性表明QSM0未引入组织对比度的显著变化。这是我们希望看到的，因为零参考值不应改变组织对比度。因此，文献中报道的使用MEDI 测量的组织对比度分析可以直接用于QSM0。

我们的临床数据表明，QSM和T2w一起能够准确识别在连续磁共振的新病灶中，未经Gd注射的MS患者中的增强性病灶。这可能作为一个潜在的临床应用，利用我们观察到的MS病灶从Gd增强性变为非增强性时磁化率的值迅速增加。QSM在MS中的一个应用是，在定期监测MS患者期间的连续磁共振中，QSM可以替代Gd增强来评估炎症活性。T1w+Gd的增强是目前评估持续CNS炎症，监测并优化炎症抑制治疗的标准方法。在最初的炎症反应之后，BBB打开，免疫细胞渗入脑中约3周时间；因此，T1w+Gd只能为病灶病理分析提供短的时间窗口。在此期间，小胶质细胞/巨噬细胞(m/M)吸收并降解髓磷脂碎片，这在QSM上反映为活性病变初期磁化率的变化不显著。然而，在BBB封闭后，免疫细胞在脑组织中保持活性。例如，m/M去除抗磁髓磷脂片段，与此同时或之后，具有顺磁性铁的m/M细胞聚集在周边和病灶内以进一步促进炎症。由于这些原因，髓磷脂碎片的清除和铁积累可能有助于QSM 观察到病变磁化率的增加。MS病变可以在几年内呈高信号，通常在QSM上有明亮的边缘。这些明亮的边缘可以被解释为铁。因此，MS中的另一项QSM应用是，通过在 Gd增强基础上包括QSM，在MS患者的磁共振方案中提供对MS中早期病灶动态的更详细的了解。

在本文提出的自动分割过程中，我们根据当前研究中的案例选择参数R。对于一些具有不同脑室解剖结构的病例，例如胎儿和小儿成像，可以进一步调整参数以对 CSF进行更稳健的分割。另外可以扩展当前算法，使得其中对CSF脑室中的“阴影伪影”的有效抑制可以被扩展至整个大脑区域。

这里M_LF为从全脑区域中提取低空间频率分量的算子。正如方程25中包含脑室 CSF的一个均匀的掩码一样，作为方程27的一种变种，M_LF可以集中于脑部上/下，左/右或是前/后半球的空间变化的低通分量。

总之，本文提出的QSM0利用CSF的自动分割与脑室CSF磁化率的平滑性约束，减弱了QSM中由各向异性磁化率产生的伪影，进而为QSM的测量值提供了一致的零参考。

5.QSM测量肝脏铁含量

在本节中,我们描述了一个用QSM方法来测量肝脏铁含量的方法，它可以补偿R2^*的混杂因素包括纤维化、脂肪和水肿。

5.1.引言

肝脏是唯一一个总是随所有类型的系统性铁沉积导致铁含量增加的器官，而系统性铁沉积可能来自输血、膳食铁吸收的增加、或者两者兼有。肝脏铁含量的测量可以作为一种评价身体铁含量多少的参考方法，可应用于对遗传和后天铁沉积疾病的诊断和管理。铁沉积导致在功能和运输系统中铁的量变化很小，几乎所有多余的铁都被螯合成顺磁性铁蛋白和含铁血黄素。磁共振对顺磁性铁非常敏感，因此已经成为铁沉积定量和去铁治疗监控的主要的非侵入性手段。横向弛豫率R2(＝1/T2), 或更快并且更敏感的R2^*(＝1/T2^*)很容易采用商用磁共振脉冲序列得到，可用以量化肝脏铁含量(LIC)。在使用R2和R2^*测量LIC方法的主要挑战是,他们易受到脂肪、纤维化和其他细胞结构变化等因素的影响,因为这些因素也会对R2和R2^*产生贡献。这些对R2和R2^*的混杂效应对水的微环境(混杂因素相互作用)和体素内混杂因素的空间分布均非常敏感，因此很难进行模型量化。因此,R2和R2^*是铁和混杂因素的复杂函数,从而很难从R2和R2^*值提取出铁的含量。应该注意的是,脂肪抑制并不影响脂肪对水的R2和R2^*的贡献。一般来说,纤维化和其他变化对R2和R2^*的贡献是很难量化和补偿。

采用多回波梯度回波磁共振成像的模图数据可得到R2^*,而用相位数据重建生成的定量磁化率图像(QSM)可以直接测量磁化率而且没有R2^*中的开花伪影。最近，通过考虑脂肪对测量信号相位的贡献，QSM应用于腹部器官的研究。磁化率是一种内在的组织特性，与肝脏铁含量线性相关。脂肪对组织磁化率的贡献可以根据化学成分线性补偿。其他组织结构细胞的变化，如纤维化和水肿，不影响磁化率。因此,在本研究中QSM可以作为一种避免R2^*的混杂因素的测量肝脏铁含量的方法。

5.2方法和材料

原理

迄今为止，QSM算法已经通过使用一种正则化惩罚边界方法成功地最小化了来自微小噪声源的条带伪影。根据贝叶斯法则,我们提出使用生物学的先验知识减少阴影伪影，具体为：在解剖结构图像中(如T1和T2加权像)信号强度均匀的感兴趣区域，铁含量也应该是均匀的(M_U，可以采用均匀强度的特征自动分辨出并通过区域变量计算探测得到)。例如,在重型地中海贫血患者,通常没有肝硬化而且铁沉积均匀分布在整个肝脏。就这种疾病而言，纤维化对磁化率的影响可以忽略不计,但是对T2 的贡献却非常大(图C2b)。因此,脂肪校正后,我们希望磁化率与M_U非常均匀，这样可使用L2范数拟合算法。此外,主动脉(感兴趣区M_A，可根据解剖图像的特征自动分辨出，尤其飞时法血管造影法的高亮信号，这方面我们具备丰富的经验)可作为零基线,因为它相对水的磁化率可以忽略不计。因此,我们提出以主动脉中磁化率是均匀的为先验知识，将平均磁化率<χ>_A设置为0。我们在大脑QSM获得的初步数据表明,侧脑室处脑脊液的磁化率变化最小，可以作为零基准以进一步提高图像质量。那么在以零基准和没有运动、流动、噪声伪影和其他错误的条件下，肝脏QSM可以表示为：

M_e为来源于模图的两值边缘掩膜(105)。脂肪磁化率可以通过脂肪图像使用 χ_fat＝X_OA·ρ_fat/ρ_OA得到，其中，ρ_fat/ρ_OA是皮下脂肪区域的ρ_OA和平均ρ_fat的比例图像。在信噪比较高的情况下,e^-ib-e^-i(d*χ)≈-i(b-d*χ),线性化第一项为偶极子模型拟合场数据。人类脂肪组织主要由油酸组成，体积磁化率X_OA＝0.75ppm。在脂肪图像中，纯脂肪可以使用最大信号区域自动识别。脂肪对磁化率的影响可以根据脂肪比例和纯脂肪体积磁化率进行补偿，这样可生成定量铁图像。肝铁生物材料的绝对的肝脏铁含量可以通过脂肪校正的磁化率χ_c并采用铁蛋白和含铁血黄素特定质量磁化率X＝1.6·10^-6m³/Kg_Fe得到:LIC＝χ_c/X。

一旦如上所述从脂肪校正的QSM估计肝脏铁浓度，则可以根据各种成像参数和铁浓度的校准曲线C(C(χ_c))估计肝脏铁浓度对R2^*的贡献。相对周围健康肝脏组织的R2^*，增加来自铁和纤维化。因此，肝纤维化生物材料可以通过从铁校正的R2^*来定量绘制：

正如我们在基于组织块的常数模型中发现的一样，方程28的L2生物学限制条件可以大大提高QSM重建的信噪比。因此,基于肝脏形态学先验知识的肝脏QSM在测量肝脏铁含量方面是非常稳健和准确的。

水模

为了对脂肪中钆(Gd)含量进行验证,我们制作了四组脂肪、水、和钆的混合物的水模，主要由24个28毫米直径的不同脂肪和水浓度的气球组成，其中，质子密度脂肪比例(PDFF)分别为：0,14.3,28.9和43.7％；Gd浓度分别为：1.25,2.5,5.0，7.5, 和10.0mmol/L。

每个脂肪-水-钆的溶液由蛋黄酱和钆混合配制。假定在293K下Gd-DTPA摩尔磁化率为326ppm L/mol，这样在室温情况下，钆溶液浓度为1.25,2.5,5.0,7.5,和10.0 mmol/L分别对应的磁化率为0.41,0.81,1.63,2.45和3.26ppm.气球被置一个注满水的塑料长方体容器中。

被试者

9个没有已知肝脏疾病的健康对照组、9个已知或疑似肝铁过载患者和9个肝脏局灶性病变患者参与这个研究，这个研究符合HIPAA的研究规范并得到伦理委员会的同意。所有参与者均提供书面知情同意书。铁沉积病人有遗传性血色素沉着症、输血的含铁血黄素沉着症、白血病或骨髓增生异常综合征。

磁共振成像

脂肪-水-钆的水模在临床3.0T磁共振成像系统(MAGNETOM Trio Tim；德国西门子)上采用头颈组合线圈扫描。水模的QSM图像和R2^*图像使用相同的3D单极型读梯度多回波梯度回波序列扫描的数据重建得到，具体参数如下：TR＝13.1ms,TE₁＝ 1.24ms,ΔTE＝2.08ms,回波数＝6,接收带宽＝1028Hz/像素,FA＝3°,采集矩阵＝256×208×48,空间分辨率＝1.5mm×1.5mm×2.0mm。

已知没有肝脏疾病和已知或怀疑肝铁沉积的病人的数据在临床1.5T磁共振成像系统(MAGNETOM Aera；德国西门子)上采用脊柱和体部组合线圈扫描得到。肝脏的QSM 图像和R2^*图像使用相同的3D双极型读出多回波梯度回波序列扫描的数据重建得到，具体参数如下：TR＝10.0ms,TE₁＝1.44ms,ΔTE＝1.36ms,回波数＝6,接收带宽＝1065Hz/像素,FA＝6°,采集矩阵＝224×168×44,空间分辨率＝1.7mm×1.7 mm×5.0mm.此外,为了减少采集时间，在前后方向采用加速因子为2的并行采集和椭圆采样。最终扫描时间为14秒，可在一次屏气期间完成。

肝脏局灶性病变患者的数据在临床3.0T磁共振成像系统(MAGNETOM Verio；德国西门子)上采集得到。肝脏的QSM图像和R2^*图像使用相同的3D单极型读出多回波梯度回波序列扫描的数据重建得到，集体参数如下：TR＝14.0ms,TE₁＝1.49ms,ΔTE ＝2.06ms,回波数＝6,接收带宽＝1502Hz/像素,FA＝10°,采集矩阵＝ 256×176×26,空间分辨率＝1.5mm×1.5mm×5.0mm。最终扫描时间为20秒，可在一次屏气期间完成。

QSM和R2^*重建

相位图像从复数磁共振数据中提取得到。首先，对场图中由于双极型读梯度梯度不理想引起的奇数和偶数回波之间的相位差异进行测量，并以多项式建模并进行校正；然后，校正后的复数磁共振信号根据组织水和脂肪的化学成分表示为：水信号ρ_w为在一个体素中的水的比例；脂肪信号，ρ_f为在一个体素中的脂肪比例，另外，相位因子对应于磁化率产生的磁场，因此，对相位图像使用基于条件跳跃的图像切割法同步进行相位解缠绕和化学位移测量,从而得到水脂比例和场图，并以此水脂信号模型中迭代计算的初始值来精确调整场图。得到的场然后使用偶极场投影法(PDF)去除背景场，得到的局部场通过形态学偶极反演算法(MEDI) 计算出磁化率图像。得到的场也可采用没有背景场去除的总场进行反演处理。

对于肝脏的有效R2^*图像,则采用文伯格-马夸特算法对多脂肪峰的多回波数据的模图进行非线性拟合得到。

数据分析

肝脏磁化率通过在QSM图像上使用感兴趣的区域(ROI)方式进行测量。背部肌肉作为零磁化率的参考。测量时ROI的区域分别有：在远离主要血管的大小为12cm²的区域，在背部肌肉附近大小为4cm²的区域。在病灶区域ROI大小为较病灶大小稍小。

脂肪对肝脏磁化率的贡献被评估并且从测量的肝脏磁化率中去除，从而得到脂肪矫正后的QSM图像。由于接收线圈引起的强度不均匀首先通过一个水平集方法进行评估，然后被用作对脂肪图像的矫正。

所有的统计学分析，使用Windows 17.0版本的SPSS进行，另外，对肝脏磁化率和R2^*也进行了线性回归分析。

5.3结果

图25显示了从脂肪-水-钆样品中测得的磁化率和R2^*分布图。钆溶液的磁化率正比于钆浓度，因此根据线性回归斜率得到Gd的摩尔磁化率(PDFF＝0,14.3,28.9,和 43.7％分别为321,318,294,和227ppm L/mol)(图25i&25k)，而脂肪对于溶液磁化率的贡献是Gd浓度为零的水模的磁化率测量值(图25i)。对于Gd浓度小于2.5mmol/L的含脂溶液，脂肪贡献的磁化率正比于脂肪的浓度(图25j)。R2^*的值与 Gd浓度有如下复杂关系：在PDFF＝0时，为线性关系；在PDFF为14.3％时，为二阶相关；在更高的浓度的PDFF时，为严重的非线性相关(图25l)。

在所有的9个对照组和18个患者中的17个人上成功进行水脂分离。失败的一列患者是因为过高的肝脏铁沉积，这导致在梯度回波成像中只有第一个回波有可探测的信号。患者的局灶性病变包有如下病变：3例血管瘤、1例肝细胞瘤、3例疑似胆管瘤或转移性瘤、1例多发性原位恶性瘤、1例为多发性肝脏转移性瘤介入治疗后。局灶性病变对R2^*有显著性影响而对脂肪校正后的QSM只有很少的影响(详情见图 26-28)。

图26显示了一例肝血管瘤病人的T2加权图像、水分布图、脂肪分布图、没有脂肪校正的QSM磁化率图、脂肪校正后的磁化率图和R2^*图。病灶在T2加权图中为高信号(图26a)，在R2^*图中信号较低(图26f)。病灶区域R2^*的平均值为13s^-1，而周围正常组织中为42s^-1，这种R2^*中病灶和邻近正常组织之间的明显对比在空间上与脂肪分布图是相关的(图26c对26f)。脂肪校正前的平均磁化率值在病灶中为 0.162ppm，在邻近肝组织中为0.203ppm(图26d)。但是在脂肪校正后，磁化率值未发现有可分辨的对比度(图26e)。

图27显示了疑似肝细胞瘤患者的T2加权图像、水分布图、脂肪分布图、没有脂肪校正的磁化率图、脂肪校正后的磁化率图和R2^*图。在T2加权图(图27a)和R2^*分布图(图27f)中，病灶均为明显的高信号。在R2^*分布图中，病灶与邻近正常的肝组织有明显对比，分别为32s^-1和121s^-1，该对比与脂肪分布图空间上是相关(图27c 对27f)。在病灶中，脂肪矫正后的平均磁化率值为0.081ppm，而邻近的肝脏正常组织中，平均磁化率值逐渐增加到0.455ppm，在空间中没有明显对比(图27e)。

图28显示了一位进行碘油介入治疗后的多发性肝脏转移瘤患者的T2加权图像、水分布图、脂肪分布图、没有脂肪校正的磁化率图、脂肪校正后的磁化率图和R2^*图。该转移瘤在R2^*分布图上显示了复杂的异质性和信号的衰减(图28f)，但是在QSM图上仅存在边缘的变化(图28d&28e)。在治疗部位脂肪矫正后的QSM中有残留的亮点出现(图28e)。

在局灶性病变的患者中，肝脏的R2^*值和脂肪比例具有相关性(图29a,R²＝0.787,P＝0.001)，但是肝脏磁化率值和脂肪比例之间却没有显著相关性(图29b,R²＝0.193, P＝0.236)。

图30显示了从没有铁沉积到中度铁沉积的3例被试的水分布图、脂肪分布图、局部场图、磁化率图和相应的R2^*图。随着铁含量的增加,磁化率和R2^*值也相应的增加。对于有轻度铁沉积的病人,在肝左叶中疑似纤维化,但这并不影响QSM的值但导致R2^*的值的增加(第二列箭头)。这组患者中,在没有纤维化的肝脏区域中，磁化率和R2^*值线性相关(R²＝0.996,P＝0.039,图30p)，而纤维化的肝区域不存在该相关性。

对于没有病灶组织和血管结构的ROI的测量,显示肝脏磁化率和R2^*测量值线性相关(R²＝0.797,P<0.0001,图31a)，而对于病灶组织，上述相关并不存在(P＝0.51,图 31b)。

5.4结论

27位被试的肝脏铁沉积分布数据说明了肝脏疾病会影响R2^*的测量，但不会影响磁化率的测量，这些疾病包括囊肿、纤维化、血管瘤、肝细胞癌、胆管癌、以及脂肪浸润。对于没有脂肪和出血的病灶，其R2^*可能减小，但是在QSM图像上是混在背景中。QSM图像上脂肪对磁化率的贡献可根据线性的组织成分模型得以校正。我们的结果说明了QSM能消除在定量肝脏铁含量时R2^*测量的混杂因素。

目前对肝脏铁含量的测量的主要挑战就是对R2和R2^*的混杂因素的校正，比如脂肪、纤维化和其他与铁过载引起的组织损伤相关的细胞变化。通过利用脂肪独特的化学位移，铁存在的脂肪分数仍能被很好定量，但是脂肪对R2和R2^*的影响则取决于脂肪和水互相影响的微观环境和脂肪在一个体素内的空间分布，这很难建模并补偿(图25l)。同样,纤维化对R2和R2^*的影响取决于水和纤维化的相互作用的微环境。从根本上讲,从MRI信号幅度可获得R2和R2^*，其依赖于细胞结构,但方式复杂很难定量。

幸运的是，通过对磁共振信号的相位严格的生物物理建模可测定组织磁化率，这就是定量磁化率成像。因为分子电子云属性线性正比于分子浓度,标准线性化学成分分解可以应用于组织磁化率,使用直接的生物物理定义就能使铁定量得以实现,而 R2和R2^*建模则不能对铁进行定量。这种情况类似于在对大脑耗氧代谢率(CMRO2)成像时对脱氧血红素的量化，QSM很大程度的简化了不理想的传统经验非线性R2^*模型。同样,QSM可以通过对其他细胞源引起的磁化率贡献进行校正，从而精确量化铁的含量。基于纤维化和其他混杂因素的磁化率相对于水几乎为0的属性，这种校正得以简化。

与其他文献一致,我们也观察到肝转移瘤、血管瘤(图26)或恶性肿瘤(图27)的R2^*值降低，而这些病变区域的磁化率值却和邻近的正常肝脏组织类似。这主要是由于在病灶区域水含量的增加导致R2^*的衰减，但是对磁化率的影响却很小。肝脏的脂肪导致了R2^*的增加(图27&28)。在治疗区域脂肪矫正后的QSM图像中，残余亮点的出现是由于在治疗后残存的碘油所致(图28)。在疑似肝细胞瘤被试中，发现特别是在邻近病灶的肝脏组织中，平均R2^*为121s^-1而且平均磁化率为0.455ppm，这些较高的值表明了铁沉积的存在(图27)。脂肪矫正后的QSM图显示，从病灶外围的正常组织到病灶中心，磁化率逐渐降低，这可能和在肿瘤边界的血管和炎症、以及肿瘤中心坏死的细胞较多有关。

在本研究中，肝脏组织被建模为顺磁性的铁(水模中的钆络合物)和脂肪。通过信号相位分析,脂肪的化学位移可通过对水脂信号模型的迭代以量化脂肪含量，继而得到QSM图像。不同磁化率的脂肪和水的混合导致了在像素内额外的场不均匀性，从而增大了表观R2^*值，而R2^*依赖于成像参数如像素大小、场强和回波时间。脂肪对R2^*具有强大的贡献可以解释病灶中R2^*值和脂肪百分比之间的相关性(图29a)。如预期的一样，在脂肪校正后的QSM图像中，测得的磁化率值和脂肪百分比(图29b) 及R2^*值(图31b)均没有相关性。在感兴趣区避开局部病变和血管结构，QSM和R2^*具有很强的相关性(图30p&31a)。

初步的实现还有一些如下不足之处可以改进。首先,研究中的患者数量相对较小。本研究中患有局灶病变的患者可以验证QSM能够克服R2^*在定量铁含量方面的混杂因素。未来应该继续招募一些由于弥散性铁沉积引起的弥散性肝纤维化病人。第二,QSM 没有考虑多峰脂肪模型。在水脂分离最终微调的步骤中，可以采用多峰IDEAL(迭代分解与回波改变和最小二乘拟合)方法。在一个实施例中，脂肪和水分数以及磁场的迭代估计可以用基于图形切割的同时相位解缠和化学位移确定来初始化。另一种初始化方法是数据采集包含同相回波(用于单峰模型的迪克森方法)，其允许估计用于初始化的R2^*和磁场。第三,对非常高的肝铁沉积患者不能进行水脂分离，而且信号衰减很快。基于放射采集超短TE技术或者偏移回波GRE序列可以用于采集有较快 T2^*信号衰减的病人的数据。第四,与免疫组织化学没有相关性,这将帮助我们解释QSM 和R2^*在对多发性硬化症白质病变研究中类似的结果。第五,在水模中Gd和脂肪的浓度较高时，似乎有较复杂的影响导致非线性磁化率(图25i)。Gd浓度较低时的数据和线性化学组分是一致。第六,屏气采集对某些病人是困难的。自由呼吸方法可以使用沿着层编码方向重复采集K空间中心数据来得到呼吸和心脏导航信号，这样可以在不增加时间的情况下实现自由呼吸导航采集。使用3D黄金比例编码顺序可用以确定如何合理采集分布于kx和ky平面、层编码平面和所有的回波的交错螺旋数据。螺旋轨迹自身具有流动补偿。回波维中冗余的信息可用于基于字典的重建以对每个回波分别进行重建。

总之,肝脏铁含量的R2^*图像会受到脂肪、纤维化、水肿和其他细胞性疾病的混杂影响，然而，QSM在测量肝脏铁含量可以避免这些混杂因素。

6.定量磁化率成像(QSM)测量骨矿化

在本节中，我们描述一种基于脂肪和皮质骨信号的详细的生物物理模型的QSM方法来测量骨矿化。

6.1.引言

磁化率是一种可以在MRI中观察到的基本组织特性。密集钙化的组织，如骨，具有很强的抗磁性。定量磁化率成像(QSM)可以提供定量的、可重复的磁化率源图像，用于评估大多数组织的健康状况，但在骨中QSM受到限制。考虑到骨密度对评估绝经后妇女和老年人骨折风险的重要性，QSM可能会成为一种非侵入性的、无电离辐射的、有用的骨健康状况成像的诊断工具。

骨中QSM面临着很大的挑战，因为它需要完整地测量整个感兴趣区域(ROI)的相位，而在梯度回波(GRE)常规的回波时间成像中，皮质骨区域信号通常非常低。虽然皮质骨中含有丰富的水分(约占体积15％)，但它大多以“结合”的形式存在，与晶体矿物结构或胶原基质相连。因此，骨中结合形式的水的表观横向弛豫时间非常短(T2^*约为300μs)，导致在常规的MRI序列上不能获得准确的相位用于QSM重建。由于这些限制，先前的工作中QSM在肌肉骨骼的应用或者侧重于软骨，或者利用分段约束估计骨的磁化率。另外一个难题是骨髓中既有脂肪质子又有水质子。

在本项工作中，我们使用TE<<1ms的超短回波时间(UTE)脉冲序列测量骨中MRI 信号。我们通过研究化学位移和R2^*成分来正确地模拟UTE和常规GRE中骨的MRI信号。最后，我们提出了一种用于从UTE-GRE数据中定量计算磁敏感的MRI方法。

6.2.理论

磁场估计是QSM中最初始的步骤。对于主要由磁化率引起相位累积的解剖区域(如在脑中)中QSM，可以从MRI相位中估计磁场。对于脂肪含量高的区域(如骨髓) 中QSM，需要根据水和脂肪的组织化学成分对复杂的MRI信号进行建模。磁化率产生的不均匀场需通过稳健的水脂信号分离来估计。

在多回波GRE序列中，包含多种成分的体素其信号-时间规律可以表示为如下的常规形式：

这里，ρ_k是体素内第k种成分在时间t＝0处的复数信号，是第k种成分第n个谱峰的相对振幅，和分别是第k种成分中第n个谱峰对应的化学位移和R2^*衰减率，是由质子经历的磁化率源引起的空间变化场,或磁化率磁场,或磁场。从MRI信号中提取多个参量对噪声非常敏感，减少等式[29]中的参数数量是稳健估计的关键步骤。我们提出下面的较少参量的模型：

在许多临床应用中，水和脂肪是主要的成分，因此等式[29]变成为

这里，是脂肪第n个谱峰的相对幅度，ρ_f和ρ_w是时间t＝0时脂肪和水的复数幅度，和分别是水和脂肪第n个谱成分的纵向衰减率。

假定所有脂肪谱峰R2^*都是接近的，则等式[30]可简化为

当前，通常忽略掉体素内水和脂肪R2^*的差异：

不幸地，在骨组织中，水是通过多孔矿物质基质分布，而脂肪在骨髓中倾向于聚集在一起，由此我们知道因此，使用等式[32]来拟合骨中MRI信号可能会在估计脂肪分数、磁化率场和QSM时出现非常大的误差。

因而，我们建议不要使用等式[32]中简化的R2^*。在皮质骨体素中，水信号是最主要的成分。在我们的实验中(见下文)，最大回波时间比包括骨髓在内的周围软组织的值短得多。所以，骨中的水是唯一在采集过程中经历显著衰减的成分。于是，我们进一步提出了一种骨特异的R2^*模型，其模型参数减少如下：

6.3.方法和材料

脉冲序列。在临床3T扫描仪(Excite HD，美国GE)上采用TE≥40μs的放射状采样的3D GRE UTE序列进行图像采集。该序列利用非选择性硬脉冲(100μs持续时间的矩形脉冲)来实现体积激发，且在每个TR中进行两次读出来加速信号采集(图 32)。在连续的TR中，对回波时间进行偏移以取得四个不同回波时间的回波，其中有两个超短回波时间的回波。

猪体膜实验。制作一个将猪蹄(长度16cm，厚度6.5cm)嵌入到1％琼脂糖凝胶 (容器高度20cm，平均直径11cm)的体膜，并且为了进行相关比较分别用一台MR 扫描仪(3T，Excite HD，美国GE)和一台临床CT系统(LightSpeed Xtra，美国 GE)来对体膜进行成像。MR成像参数如下：TE＝0.04,0.24,3.0,4.0ms，TR＝12ms， FA＝15°，FOV＝18cm，每个回波采集32000条放射状投影线，体素大小 0.7×0.7×1.4mm³，BW＝±62.5kHz，总采集时间＝13分钟。CT数据使用以下参数进行采集：120kVp，200mA，0.625mm层厚，512×512矩阵，达到接近各向同性的分辨率。

志愿者MR成像。在获得书面知情同意后，扫描7名健康志愿者(6名男性，1名女性年龄范围25-32岁，平均年龄＝30)。所有检查均根据机构审查委员会批准的协议进行，并严格服从HIPAA标准。在3T MRI系统上使用8通道发射-接收膝盖线圈和8通道发射-接收手腕线圈来采集图像(11个股骨远端，1个桡骨远端和1个胫骨近端)。使用的3D UTE GRE序列参数如下：TE＝0.04,0.24,3.0,4.0ms，TR＝12ms (在该实验中，两个TR分别获得第一/第三和第二/第四回波)，FA＝15°，FOV＝18cm， 32000放射状投影线，体素大小＝0.7×0.7×1.4mm³，所有回波总采集时间＝13 分钟。

图像处理和分析。利用网格化来重建每个放射状MR数据，网格化是将放射状采集的信号插值到笛卡尔网格。在本文中，使用凯泽-贝塞尔内核和最小-最大插值法的NUFFT来进行网格化。由于我们的放射状轨迹采样密度不均匀，因此在进行网格化之前需要对测量的信号进行密度补偿。每个数据集中每个线圈和回波分别单独重建。

对于每个线圈，每个后续回波的相位都减去第一个回波的相位，以消除线圈敏感度产生的相位。然后对于每个回波，相位相减后的复合图像在所有线圈上进行求总和。这个总和信号的相位用作对应回波的相位。在这个过程之后，第一个回波的相位就变为零。这个过程用数学的表达可以写成：

其中是对应于第m个回波的复数数据，是这个回波的幅度，是它的相位。n用于表示第n个线圈，I^c是I的共轭。

随后使用基于IDEAL的方法来获得水、脂肪和不均匀性场。对于这个迭代方法，为了避免收敛到局部最小值，用一个相当接近于真实解的不均匀性场作为初始化猜测是很重要的。在本文中，假定只有单一的化学位移谱峰f＝-3.5ppm·γB⁰Hz，使用基于图形切割的相位解缠和化学位移去除来进行初始场图估计。通过单指数拟合获得初始的R2^*估计。假设使用不同脂肪化学谱峰模型-单峰模型(Δf＝-3.46ppm) 和多峰模型)，复数数据可以根据方程 [33]和[32]来进行拟合。然后将得到的不均匀性场b用基于形态学的偶极子解卷积方法(MEDI)进行求逆磁化率成像，其中偶极场投影(PDF)方法用于背景场去除。

将猪蹄的CT图像进行重采样并使用FSL工具箱中的FLIRT算法将其配准到重建的磁化率图上。对配准后的图像分析ROI中CT(以亨氏单位(HU))与MR(计算的敏感性值)的相关性。为了分析它们的相关性，手动勾画层面内的58个感兴趣区域包括骨小梁、皮质区域的掌骨、皮质区域的指骨和屈肌腱。ROI中QSM(肌肉组织作为参考)和CT体积平均值用于做相关分析。

6.4.结果

猪体膜结果。成功重建出了样品的磁化率图。图34显示了CT和QSM图像的对照以及ROI分析的结果(图33)。我们可以看到，QSM中的抗磁区域与CT图像中的高 HU值区域之间有良好的对应关系，这一点与得到的强的线性相关结果(相关系数R²＝ 0.77)是一致的。

志愿者结果。使用本文提出的水特异的和脂肪多峰模型的方法在所有志愿者中都能够成功重建QSM(参见图35中的示例)。图36显示了用不同技术估计的场图及其对应的计算出的磁化率图。从使用传统的估计方法计算获得的参数图中可以观察到骨和肌腱区域的脂肪分数系统性高估(高达80％，与健康骨中可忽略的脂质含量相比)。这些误差出现在场图及随后的磁化率图中，导致骨和肌腱的磁化率值出现很大的误差–尤其是股骨皮质区域，当假设脂肪峰为单一峰值模型时或者脂肪和水之间用共同的R2^*时，骨的磁化率错误地表现为顺磁性。还应该指出的是，我们提出的方法使小梁骨具有最佳的可视效果。

图37显示一名健康受试者膝关节重建的QSM的最小强度投影图。可以注意到的是股骨和胫骨干厚皮质区域都显示出均匀抗磁性。在整个FOV中，QSM很好地描绘了小梁，显示出骺板清晰的外观。

志愿者结果显示皮质骨磁化率值在个体内和个体间都有好的一致性(个体间平均值-1.4±0.2ppm)。

6.5.结论

我们的数据表明了骨磁化率成像的可行性，包括皮质骨，骨小梁和骨髓。通过利用放射状采集的UTE梯度回波(GRE)序列获取快速衰减的骨中水信号并建立新颖的骨特异的R2^*模型以实现准确的水脂分离来完善骨中QSM。体外猪蹄的结果表明QSM 值与CT有极好的相关性。通过与当前脂肪和水同一R2^*的标准信号模型重建的QSM 进行比较，骨特异的R2^*模型降低了皮质骨内的脂肪分数估计的误差并且在7名健康志愿者的皮质骨中得出了对所有骨组织成分都有很好的描述的合理磁化率值。

使用传统的IDEAL方法从UTE和常规回波GRE数据中估计得场都是不容易的，倾向于受到噪声传播和R2^*偏差的干扰。从GRE数据中稳健获得生物物理参量要求使用最小数量参量的适当的信号模型。因为将骨中快速信号衰减错误解释为由于高脂肪分数发生的散相，所以常规的GRE信号模型在水和脂肪中使用同样的R2^*往往不能正确地估计场图，因而就会产生具有严重错误值的骨QSM(例如顺磁性骨)。我们意识到，在本研究所获得的回波时间中，脂肪信号的R2^*衰减可忽略不计，但在钙化基质内分散的水的快速R2^*衰减不能忽略。这就需要我们利用骨特异的R2^*模型，该模型准确地反映了潜在的组织物理学特性，将信号衰减归因于水成分。猪蹄和志愿者数据都表明，这种骨特异的R2^*模型与迭代最小二乘拟合方法相结合，可以提供化学成分分布的可靠估计以及生物学上可信的场图。此外，应用脂肪多谱峰模型进一步提高了水脂分离的质量，因此也进一步提高了磁化率图的质量。

以前的研究表明，骨的磁化率可以通过局部磁场的高度正则化分段求逆来估计，这种方式当组织磁化率存在局部变化是不可取的。掩膜的磁化率迭代计算也可用于估计组织中强磁化率源的分布，但迭代趋向于迫使高磁化率结构均匀和低估。尤其是，当这些尝试并没有恰当地考虑MRI信号产生潜在的脂肪和水生物物理学，也没有考虑使用如本文所述UTE采集骨信号。对于密度远低于皮质骨的松质骨，传统的回波时间可以直接用于QSM。

骨密度(BMD)评估对于骨质疏松症的诊断是极其重要的。目前，BMD使用双能X 线骨密度仪(DXA)进行测量，这种测量方式可能会受到包括退行性变化、骨尺寸和重叠的解剖结构等因素的干扰。定量CT(QCT)被提出用于克服DXA中的限制，但它需要更高的X射线剂量。还有许多可用于表征骨组织的非侵入性非电离辐射的方法。定量骨超声特征是骨质疏松患者临床评估的重要方式。一些研究报道使用沿着长骨轴的超声波的轴向传播速度作为监测骨质疏松患者的骨折愈合和骨评估的生物学标志物。基于MRI的技术(例如ZTE和UTE)依赖于骨中有机基质上成像的结合水。获得的这些信号可以提供关于骨矿化密度和骨矿化孔隙度的信息；这需要分离有着显著差异的驰豫常数的结合水和自由水信号成分。使用双成分分析的UTE成像已经成功应用于先前的体外和体内皮质骨成像。

常规MRI可能可用于评估骨质量但不能用于评估骨矿物化密度。这里提出的基于MRI的骨QSM有可能提供没有X射线辐射的BMD评估。因此，我们的结果值得在今后进行研究，在骨质疏松症患者中比较骨QSM与DXA，以建立准确的没有X线辐射的骨密度评估标准用于预测骨折风险和指导骨治疗。由于QSM作为生物标志物与其他非电离辐射测量方法的综合比较超出了本文的范围，因此需要另外的研究以进一步评估骨QSM的临床的实用性。

UTE在QSM中的初步实现可以在考虑以下的情况下进行改进。使用贝叶斯MRI的数据采集加速策略可以大大减少扫描时间。初步实现的QSM还假定，在采集最后一个回波时，水是唯一显著衰减的物质。虽然在存在骨基质结合水的许多实际情况下这种假设是成立的，但是这种假设可能在强的场不均匀性的条件下不成立(例如，在腹腔附近的成像)。这些场变化可导致强的自旋散相，因此导致软组织区域中快速 R2^*衰减。在这种情况下，可能需要高阶的匀场。在本文中没有考虑到的额外的限制因素包括由于脂肪的偏共振频率造成的图像畸变。由于带宽的限制，在本研究中使用的放射状采集方式，偏共振频率表现出幅度和相位图像中水和脂组织边界的模糊。预计该模糊会影响最终的QSM结果，因而可能需要应用适用于非笛卡尔采集的偏共振校正方法。最后，正如之前报道的那样，QSM易于低估强的磁化率源，并且可能需要进一步调整成像方案(矩阵大小，分辨率)来最小化估计误差。

本研究中获得的结果与先前在文献中报道的结果的对比表明本文提出的方法可能低估了骨的磁化率。例如，在作者报告中，体外实验中获得的骨的磁化率 χ_bone≈-2.4ppm。此外，在体内磁化率估计中，不同重建方法得到的磁化率值有相当宽范围χ_bone∈{-1.8，-2.3}ppm。最值得注意的是，这种低估可以在关节周围观察到(图37)。虽然这个问题需要进一步研究，但是这种低估的确切性质可能与部分容积效应和上述技术的局限性有关。

与CT衰减系数成正比例的骨的磁化率可直接用于校正正电子发射断层成像(PET)重建中的衰减。这对于组合的MRI和PET成像设备(MR/PET)特别有用。

总而言之，使用UTE和常规TE的梯度回波采集方式与骨特异的R2^*信号模型的组合对于整个骨横截面的定量磁化率图是可行的。

7.QSM脑氧代谢率成像

上述技术可以实现许多不同的应用。例如，以下描述了使用一种或多种这些技术在数值仿真，模型实验和人脑成像中的应用。

7.1简介

有氧呼吸是大脑能量供应的主要方式：大脑占成年人体重的2％，却消耗了全身氧气供应的20％。因此，氧缺乏症很容易损伤脑组织，如阿尔茨海默病(AD)、多发性硬化及缺血脑卒中的缺氧。定量CMRO2成像(qmCMRO2)对评估这些脑部疾病特别有价值。例如，在缺血性脑卒中治疗中，识别可以挽救的缺血半暗带是关键。采用脑血流量(CBF)减少来定义缺血梗死核心是有问题的。扩散/灌注不匹配模型经证实并不能改善临床结果。梗塞周围区域中的氧提取分数(OEF)升高预示组织的最终结局也存在不足。CMRO2(＝CBF*OEF*动脉氧含量)直接预示着神经元的死亡，可以克服这些不足。另一个例子是AD有氧糖酵解的改变。有氧糖酵解的测量部分取决于氧代谢的定量测量。AD在人口老龄化中的日益增长引发了对定量CMRO2研究的兴趣的增加。

除了检查疾病，定量CMRO2还是研究基础大脑功能的有用工具。例如，发现脑默认网络可以说是我们在过去二十年中对脑功能知识中最重要的进步，并且它依赖于使用定量CMRO2成像。

华盛顿大学的研究人员率先开发了基于¹⁵O示踪剂的PET成像技术，该技术仍然是定量脑循环和脑代谢成像的金标准。基于¹⁵O-氧气，¹⁵O-一氧化碳，¹⁵O-水和¹⁸F- 脱氧葡萄糖(FDG)的PET提供了测量健康志愿者和病人脑循环和脑代谢成像的唯一真正定量成像方法。尽管血氧水平依赖性(BOLD)功能磁共振成像被广泛用于脑功能的定性评估，但定量成像还存在困难。不幸的是，¹⁵O的半衰期只有短短2分钟，需要现场生产，限制了其在临床和研究中的广泛应用。

为了解决¹⁵O PET的局限性，已经通过MRI研究进行CMRO2定量成像。CBF可以通过空间选择性射频脉冲RF来产生以产生动脉自旋标记(ASL)或使用造影剂标记，但qmCMRO2需要从MRI计算额外的脱氧血红素浓度[dH]。到目前为止，一直集中在对[dH]有复杂依赖性的幅度信号上进行研究。一个有见识的想法是通过使用速度选择性标记来估计小静脉中的T2来计算[dH]，例如氧摄取和组织消耗定量成像 (QUIXOTIC)及其延伸的方法。但是，这种方法存在流速捕获，灵敏度差以及动脉血和脑脊液污染等问题。[dH]依赖的MRI幅度信号已被严格地模拟为近似血管系统，即定量BOLD(qBOLD)方法，但是在模拟中存在误差，因而用于估计[dH]的结果同样也差。最广泛研究的方法是采用经校准的fMRI(cfMRI)，其经验地模拟梯度回波(GRE) R2^*信号的衰减。来自R2^*的[dH]估计需要在两个血管干预期间采集数据，典型的方法是缺氧和高碳酸血症，这使得cfMRI在研究中难以执行并且难以在临床中使用。总之，这些基于MRI的定量CRMO2技术敏感度差，使用不确定的生物学假设，或需要可能复杂的血管干预。此外，基于MRI的定量CMRO2成像尚未使用PET进行验证。

7.2原理

QSM消除了当前基于的校准fMRI(cfMRI)中的经验参数

cfMRI是CMRO2定量成像目前应用最广泛的MRI方法。在cfMRI中，通过将GRE 幅度信号建模为指数衰减，即来评估静脉脱氧血红素浓度[dH]：

这里v是脑静脉脑血容量；r₀是[dH]＝0时的信号衰减。公式34中的a和β是经验参数，没有明确的物理意义，因为它们取决于场强和血管几何形状，并且需要使用血管干预例如缺氧和高碳酸血症的额外数据采集。因此，采用cfMRI方法的定量CMRO2成像在临床实践中很难执行。

我们通过利用通常被丢弃的GRE相位信号来解决这个问题，该信号可以用于定量磁化率成像(QSM)。用于估计[dH]的磁化率χ(QSM)模型去掉了公式34中的经验参数a和β：

χ＝v·X_dH·[dH]+η₀ [35]

这里X_dH＝151.054ppb ml/μmol是dH分子的磁化率；η₀是[dH]＝0时的磁化率。因此，QSM可以明显简化基于cfMRI实验的定量CMRO2成像。

从GRE相位数据的QSM分析中定量CMRO2成像(公式35)需要标准化处理，可以通过血管干预采用qBOLD来描述GRE幅度数据和最小局部变化(MLV)来实现CMRO2 定量成像。通过贝叶斯推断来优化解决CMRO2成像中qBOLD中血管干预和MLV的不确定性。因此，定量CMRO2成像可以从GRE幅度、相位数据以及ASL数据中获得，这些数据几乎可以在每台MRI扫描仪上轻松获取。

使用QSM cfMRI+MLV+Leenders实现定量CMRO2成像(qmCMRO2)

利用贝叶斯计算CMRO2时，我们将数据噪声和先前的不确定性都近似为高斯分布。qmCMRO2需要测量脑血流量(CBF≡f)，可通过ASL获得；以及需要计算静脉脱氧血红蛋白的浓度[dH]和动脉脱氧血红蛋白[dH]_a之间的差异，即氧摄取分数 (OEF≡σ＝([dH]-[dH]_a)/H)，其中H＝7.53μmol/ml，组织血液中血红素的浓度(血细胞Hct＝0.3567(直窦))。对于98％动脉氧饱和度，[dH]_a＝0.02H ＝0.15μmol/ml。此外，当K＝n_HbM_HbC_HbHctX_dH＝.0084L/mol并且 χ₀＝Kv[dH]_a/H+η₀时，CMRO2方程和公式35变为：

CMRO₂≡ε＝Hfσ，χ＝Kvσ+χ0. [36]

为了从ASL数据f和GRE相位得出的χ来优化估计CMRO2ε，我们考虑三个生物学先验假设：1)Leenders假设相关的脑静脉血流量v与流速数据f:v＝L(f)。2) 在比布罗德曼区域或具有类似细胞构筑的区域小得多的模块中，灰质和白质可能有均匀的但不同的CMRO2值。CMRO2的这种最小局部变异(MLV)与¹⁵O PET CMRO2文献和fMRI激活图一致。3)灰/白质CMRO2的总和等于在大静脉中测量的全局CMRO2。这三个先验假设可以被纳入能量最小化问题来估计ε(和v)，即：

这里，w_p是QSM中噪声的倒数乘以f，Γ＝K/H,是在一个模块中ε的组织特异性平均值，是静脉中测量的全局CMRO2平均值。公式37右边的第一项是公式 36中QSM数据保真度，第二项是血液流速和血流量之间的Leenders关系，第三项是 MLV，第四项是大静脉测量的氧饱和度的一致性。我们通过v＝L(f), (模块常数ε)。来进行初步测试。然后将公式37简化为通过平移模块来减少块形伪影，我们获得了非常令人鼓舞的不需要血管干预的MLV结果，与咖啡因血管干预后的结果相当(图39)。除了灰/白质分割之外，MLV还可以包括分割每个模块内的病灶和正常出现的组织。我们计划通过根据差异原理和L曲线分析仔细选择标准化参数来求解公式37。预处理技术可以用于加速求解。 Broyden-Fletcher-Goldfarb-Shanno-Bound(L-BFGS-B)型方法将利用OEF的生理特性：σ∈[0，1]。

先前的Leenders先验假设可能适用于研究健康受试者的神经生理学，并可能使用格拉布关系扩展到血管干预等情况。病理状态可能会使血流量和血容量减少，因此，需要替代方法来估计静脉CBVv。动态对比增强MRI仅提供总CBV。血管空间成分(VASO)方法可以提供v但在患者中应用可能存在问题。我们将探讨qBOLD方法，以克服Leenders方法的限制。

采用QSM cfMRI+qBOLD+MLV方法进行qmCMRO2成像优化

qBOLD方法使用先前的血管内随机方向的血流，可以估计大多数健康和病理组织体素内的静脉。qBOLD已应用于使用多个非对称自旋回波或GRE来生成CMRO2。qBOLD 仅利用信号幅度来计算多个参数：与CMRO2相关的σ和v，以及与CMRO2无关的 S₀，R₂，χ₀(在下面的公式38中定义)。由于qBOLD信噪比差，血管系统的假设存在错误，我们提出通过QSM利用相位数据和通过qBOLD的幅度数据。qBOLD假设以及 MLV和其他标准化中的不确定性可以通过贝叶斯统计最优化地生成公式37，方法如下：

这里，为qBOLD模型，G＝(4π/3)γB₀Δχ₀Hct＝51.6Hz，S₀为未知的自旋密度信号，R_c通过场不均匀性(包括体素敏感性函数效应)来求解，R₂为未知的组织T2弛豫率，R_b(σ)为依赖于脱氧血红蛋白[dH]和的T2弛豫率，为标准化，包括公式37里的MLV、S₀以及R₂，χ₀解剖结构的平滑。MLV的想法可以扩展到信号直方图空间，例如，对所有体素具有类似信号行为的比如类似的时间衰减趋势，施加最小的方差的或相似的σ和R_2。公式38将用QSM中的非线性高斯-牛顿最优化法来求解。我们还将研究使用原始对偶算法来加速计算。

7.3材料和方法

该研究经当地伦理审查委员会批准。收集健康志愿者(n＝11，其中女性1人，平均年龄34±12岁)，使用8通道脑接收线圈在3T磁共振扫描仪(HDxt，美国GE) 上行脑MRI扫描。指示所有受试者在MRI前24小时避免摄入咖啡因或酒精。

在口服200mg咖啡因之前和之后25分钟分别行MRI扫描。扫描方案包括3D快速自旋回波(FSE)动脉自旋标记(ASL)序列，3D多回波梯度回波(SPGR)序列和T1w SPGR序列(BRAVO)。

3DFSE ASL序列参数为：20cm FOV，3mm的各向同性分辨率，标记时间1500ms，标记后延迟时间1525ms，读出带宽62.5kHz，8个交错的螺旋采样，每次512个读出点，扫描层数35，TE 10.5ms，TR 4796ms，激励次数3和扫描时间6min。

3D SPGR序列参数包括：平面内分辨率0.52mm，层厚2mm，与3D FSE ASL序列相同的扫描范围，7个等距回波，TE1＝4.3ms，回波间隔7.9ms，TR 56.6ms，读出带宽62.5kHz，翻转角15°，扫描时间7分钟。脉冲序列在所有三个方向采用流动补偿技术。

反转T1w SPGR序列参数为：1.2mm各向同性分辨率，与3D FSE ASL序列相同的扫描范围，TE 2.92ms，TR 7.68ms，准备时间450ms，读出带宽±19.5kHz，翻转角 15°，扫描时间2分钟。

CBF和QSM重建

使用FuncTool软件包(美国GE)从ASL数据生成CBF图(ml/100g/min)。对于QSM重建，使用GRE相位的自适应二次拟合，然后在偶极场上进行投影，以获得局部场强图。使用MEDI算法从局部场强图计算得到磁化率图。磁化率值以从脑脊液(CSF)中的水作为参考(见下文)。使用FSL FLIRT算法将所有图像与QSM图进行配准。

CMRO₂和OEF重建

根据氧分子的质量守恒，脑氧代谢率CMRO2(μmol/100g/min)可以表示为：

CMRO₂＝CBF·([dH]_v-[dH]_a) [39]

其中CBF是脑血流量(ml/100g/min)，[dH]_v、[dH]_a分别是静脉和动脉血液中的脱氧血红素浓度(μmol/ml)。

氧摄取分数(OEF)是从动脉血中摄取的氧的百分比，可以表示为：

组织血液中血红素的比容[H]＝7.53μmol/ml，假设组织血液中红细胞的比容Hct_tissue＝0.357，这是通过假设直窦中的红细胞比容Hct_ss＝0.47，大静脉与脑组织的红细胞比例约0.759计算而来。假设动脉氧饱和度为98％，[dH]_a＝0.15μmol/ml。

计算组织血液中血红素浓度([H])

其中4代表每个血红蛋白Hb含有4个血红素，Hct_tissue＝0.357，ρ_RBC，Hb＝0.34g/ml为红细胞中血红蛋白的浓度，M_Hb＝64450×10^-6g/μmol是脱氧血红蛋白dHb的摩尔质量。

χ_o计算

χ_o＝CBV·X_ba+CBV_a·X_dH，mol·[dH]_a

这里CBV_a＝CBV-CBV_v＝0.23CBV为动脉血容积分数，[dH]_a＝0.15μmol/ml是假设氧饱和度为98％时的动脉血的脱氧血红蛋白含量，X_ba是100％氧饱和是的磁化率。

X_ba＝Hct_tissueδ_Hb·X_oHb+(1-Hct_tissueδH_b)X_p＝-108.3ppb

其中δ_Hb是红细胞内血红蛋白的容积分数。

这里ρ_RBC，Hb＝0.34g/ml是一个红细胞内血红蛋白的质量分数，ρ_Hb＝1.335g/ml是纯聚合物中的血红蛋白密度，X_oHb是脱氧血红蛋白的容积磁化率， X_oHb＝-4πρ_Hb(0.587×10^-6ml/g)-(-9035ppb)＝-813ppb。 0.587×10^-6ml/g是相对于珠蛋白采用CGS单位的球蛋白的磁化率。4π是指从cgs 单位转换到si单位。-9035ppb是相对于珠蛋白的水容积磁化率，X_p＝-37.7ppb是氧合血红蛋白和血浆的容积磁化率。公式39、40中的[dH]_v和相应的CMRO2可以通过QSM来测量。

这里χ是使用QSM测量得到的磁化率，χ₀反应纯氧血和动脉脱氧血红蛋白的磁化率。 CBV_v是静脉血的容积分数，大约为总的脑血容量CBV的77％。 X_dH，mol＝151.054ppb ml/μmol是脱氧血红蛋白的摩尔磁化率。

公式41中，CMRO₂和χ_nb是两个未知数。在之前的研究中，测量了基础状态下及脑干预状态下的磁化率及脑血流量。因此，CMRO₂,χ_nb,以及随后的OEF可以求解得到。然而，像咖啡因的摄入之类的CBF干预非常耗时，在病人身上应用存在许多困难。

为解决这一局限性，将QSM和CBF图分割成小的模块，并使用T1加权图像将每个模块进一步分割成两种组织类型(灰质和白质)。在每个模块和每种组织类型中， χ_nb和CMRO2被假定为恒定的。这假定在一个小的局部区域内，相同的组织类型具有相似的生物和化学组成。因此，每个模块和组织类型的CMRO2和χ_nb可以以线性系统格式组织成未知数x的向量：

这里ψ_i＝CBV_v，i·X_dH，mol/CBF_i，φ_i＝χ_i-CBV_i·(X_ba+X_dH，mol·[dH]_a)。指数 i＝1，...，N代表模块内某一确定组织内第i个体素，N为模块内的体素总数。 X_ba＝-108.3ppb为总氧合血液的磁化率。

随后，该算法将被定义为最小局部方差(MLV)。为了减少计算错误并提高拟合速度，因而应用了正确的预处理技术。新的方法定义为A_Px_P＝b，其中A_P＝AP， x_P＝P^-1x。P为预处理，定义为P＝diag(CMRO_2，ub，avg(χ))，其被设计为允许x_P的体素具有相似的幅度顺序。CMRO_2，ub是使用公式39和40在OEF＝1时在每个模块内计算的CMRO2的上限。avg(χ)表示模块中具有相同组织类型的体素的磁化率的平均值。

对于确定的组织类型(灰质或白质)的解决方案是通过使用具有生理界限的有限存储Broyden-Fletcher-Goldfarb-Shanno-Bound(L-BFGS-B)算法使以下函数最小化来获得CMRO2图：

A是diag(A_P)，对于确定的组织类型和沿对角线的每个模块的对角矩阵A_P。x, b及CMRO₂分别是模块内确定组织类型的x_P,b及CMRO₂的串联的列向量。在第二项中，是给定组织类型的所有模块的平均CMRO2。该项基于质量守恒强加了全脑生理性限制：由CMRO2图计算的全局CMRO2应与从直窦计算的全局CMRO2相似。这里，是根据公式41计算的。全脑血流量通过对全脑体素的CBF平均值来估计。直窦静脉中的磁化率通过一个有10年经验(AG)的放射科医生绘制的直窦ROI计算而来，因为，这些直窦ROI中只有静脉血，并依据CBV＝CBV_v＝100％， χ_nb＝0，Hct_ss＝0.47。λ是全局约束条件下的标准化加权参数，根据三个随机选择的研究对象中通过L曲线分析获得的平均值选择的。

使用FSL FAST算法在反转T1w SPGR图像上获取全脑灰质(GM)，白质(WM)和 CSF模板。这些模板被进一步分成6x 6x 6mm³各向同性的模块。这种基于模块的计算是分别对GM和WM计算，在每个模块内仅选择那些属于相应组织类型的体素。当求解公式38时，这允许适当的整合。为了减小由于特定模块造成的潜在误差，将模块沿(3D)对角线方向平移4次，每次长度为1.5mm。数据被重新处理并且平均以后从而获得最终的图像。

CMRO2，OEF和χ_nb也按如下方法求解。简单来说，假设χ_nb在干预前后保持不变，对基础状态和干预状态下每个体素的磁化率和CBF值进行配对，以使用提出的具有相同OEF/CMRO2范围和全局CMRO2约束的LBFGS-B算法为每个体素求解OEF和χ_nb。 CMRO2在两种状态之间的变化率在直窦中测量并最小化。

对于该可行性研究，CBV是由CBF和CBV的线性相关关系计算而来，基于¹⁵O-PET 在健康受试者(n＝34)灰质和白质测量的结果，CBV＝(0.0723CBF+1.144)/100。使用0.38的格拉布’s指数从基础状态CBV估算咖啡因干预后的CBV，是因为来自 PET数据的经验性的CBF和CBV关系数据是从没有血管干预的健康受试者中获取的，可能不适合计算咖啡因干预后CBV的计算。如文献所述，干预后CBV的变化假定在均发生动脉内。

统计学分析

与前面相同的神经放射医生(A.G.)将皮质区域的灰质(CGM)分割成双侧前(ACA)、大脑中(MCA)和大脑后动脉(PCA)血管区域(VT)。对于每个ROI，测量CMRO2和 OEF值的平均值和标准差。使用先前提出的方法使用布兰德-奥特曼和配对样本t检验比较不同血管区域VT中的CMRO2和OEF值。P<0.05被认为统计学差异显著。所有数据均以均数±标准差表示。

数值仿真

进行数学模拟以检测是否有以下因素造成CMRO2计算的错误：1)CBF错误，2) 部分容积效应，3)QSM错误。重建一个与本研究相同的矩阵大小为8×8×5的模块。上半部分和下半部分分别设置为白质(WM)和灰质(GM)。输入为CMRO2，χ_nb和CBF (参考)：GM的CMRO2:181,WM的CMRO2:138(μmol/100g/min),GM的χ_nb:0.02, WM的χ_nb:0.45(ppb)。CBF值被赋予基于体素的正常随机值，WM为49.5±4.3， GM为64.8±3.6(ml/100g/min)。这些值是受试者的平均值，计算CBF的标准差，全脑GM和WM标准差并除以5，以获得小的脑区中相同组织类型的合理的拟合结果。随后用给定的CMRO2，χ_nb和CBF值以基于体素的方式计算QSM。1)为了研究 CMRO2图中的CBF误差，我们考虑了两种情况：1-1)低信噪比CBF和1-2)CBF偏差。对于低信噪比情况，噪声被添加到CBF(SNR：4.5)，随后执行MLV。这被重复 100次以获得平均值和标准差。CBF的SNR是由在尾状核测量的平均值与标准差的比值计算而来，4.5±1.7(n＝5)。2)为了研究单个体素内的部分容积效应，将GM 边缘的QSM和CBF值分别与直接相邻的WM的QSM和CBF值进行平均。这意味着GM 边缘的体素实际上由一半GM和另一半WM组成，但被视为GM。3)为了估计CMRO2 计算中QSM的误差，将噪声添加到QSM(SNR 3.1，与CBF SNR相同的方式计算)。重复100次以获得平均值和标准差。对于所有受试者，本研究中使用的全局CMRO2约束权重λ设置为500。

7.4结果

图38显示了使用MLV方法的三个随机选择的受试者的平均L曲线。λ(公式 42)选择为500。在CGM中，在咖啡因干预后CBF降低16.3±7.0ml/100g/min (p<0.01)。使用基于咖啡因干预的方法和基于MLV的方法测得CGM中的CMRO2值分别为142±16.5和139±14.8μmol/100g/min。相应的OEFbase分别为 33.0±4.0％和31.8±3.2％，OEFchal分别为39.0±4.9％和39.3±4.1％。

图39显示了一位健康志愿者使用基于咖啡因干预和基于MLV的方法重建得到的全脑CMRO2图，紧接着为反转T1w SPGR图像作为解剖参考。CMRO2图显示良好的一致性，灰白质对比度与解剖图像基本相符。

图40显示了第二个受试者用咖啡因干预和MLV方法重建的OEFbase，OEFchal和基础状态CMRO2图。OEF图也显示出不同重建方法之间的良好一致性。很容易理解， OEF的全局增加与咖啡因引起的全局CBF减少一致。

图41显示了使用咖啡因干预和基于MLV的方法得到的CMRO2和OEF在不同血管区域(VT)中的布兰德-奥特曼图。这些比较显示两次测量的平均值的差异<4％，没有统计学意义(p>0.05)。

数值仿真表明CMRO2的计算误差受CBF的低信噪比和CBF的偏移影响。在低信噪比CBF(SNR 4.5)的情况下，GM的CMRO2误差为6.0％，WM的误差为4.4％。在 CBF偏移(低估10％)的情况下，GM和WM的相对CMRO2误差为10％。就GM边缘的部分容积效应而言，GM的相对CMRO2差异为3.0％，WM为2.2×10-5％。在QSM 噪声转移到CMRO2(SNR3.1)的情况下，由此产生的CMRO2误差对于WM为0.1％， GM为1.5％。

7.5结论

我们的初步结果说明基于QSM和CBF使用MLV方法进行无创CMRO2成像而无需血管干预的可行性。CMRO2测量结果在健康受试者中使用血管干预和MLV方法之间一致性好。去除血管干预的程序明显减少了扫描方案的复杂性，整体检查时间以及使用改变血流量的药物的潜在风险。

之前报道的基于R2^*的CMRO2成像方法受R2^*对[dH]建模的限制。需要对这些模型参数进行校准，其可能随扫描参数(如场强)而变化。因此，定量CMRO2成像需要在三种脑状态(常规氧状态，高碳酸血症和高压氧状态)下测量。然而，将兴奋剂用于大脑产生三种不同的状态是非常麻烦的。磁化率与[dH]具有完全确定的线性关系，消除了R2^*模型的参数校准，并将CMRO2成像所需的大脑状态数量减少为两个。需要两个大脑状态来区分血液和非血液组织的磁化率值，并且可以使用血管干预来实现，例如通过服用咖啡因或使用过度通气的方法。但是，即使是单一的干预措施，在实践中也可能会很麻烦。使用MLV定量CMRO2成像可以在没有这样的干预的情况下实现，使得CMRO2在临床实践中更加可行。

在本研究中使用的MLV使用6×6×6mm³的各向同性模块，这样使得各向同性3x3x3mm³的CBF体素均在每个模块内。1909年，布罗德曼根据其细胞结构将人类大脑皮层分为44个区域，后来证明这些区域与各种神经功能有关。在这项研究中，大脑皮层分为4000多个区域。因此，假设相同的模块和相同组织类型内的CMRO2 和χ_nb是恒定的被认为是合理的。然后MLV方法实现多个数据点唯一地计算每个模块内的两个未知数CMRO2和χ_nb。从CBF和QSM输入到CMRO2输出存在大量的噪声转移，可以使用全局CMRO2约束来降低，即公式42中的标准化。

本项研究也有几个局限。在MLV中，假设χ_nb和CMRO2在每个模块和每种组织类型内是恒定的。因此，MLV计算出的CMRO2值是模块恒定的并且取决于准确的组织分割。通过考虑与生理状态相关的大脑分割可以进一步改进MLV。选择6×6×6mm³各向同性模块的理由是，对于这组数据，每个模块都包含2×2×2CBF和12×12×3 QSM，两者都被认为足够大以避免不能够定性，但又不是太大以避免不能保持组织均匀性的假设。为了在这些需求之间保持平衡，选择了一个合适的模块的大小。为了研究模块大小的影响，我们对一个受试者使用了不同的模块大小(2×2×2mm³和 9×9×9mm³)。对于2mm，6mm和9mm模块大小，皮层GM的CMRO2值分别为 119.0±107.5,123.5±83.2和124.8±64.2μmol/100g/min。随着模块大小的增加，CMRO2的标准差减小。如本文所做的那样，对几个子模块的平移进行平均，可以抑制所产生的块状表现。尽管如此，由于对每种组织区域分别进行计算，所以不同组织区域之间仍然存在急剧的转变(图39)。组织分割对于疾病状态下可能存在困难，这可能会导致错误。关于不同病人潜在的复杂的病理状态下的进一步研究，例如获取病灶内不同区域之间的变异存在一定困难，这也需要进一步研究。MLV另一个局限是使用基于模块的方法可能会平滑具有高度脱氧特性的静脉。如果相同组织类型的某些体素具有较多的高度脱氧静脉，则可能被静脉脱氧含量较低的其他体素平滑。避免这个问题的一种可能的方法是在实施MLV时从灰质和白质模板中除掉高度脱氧的大静脉。缺乏对Hct测量是另一个局限。我们假设Hctss＝0.47。我们通过重复计算人Hct范围的预期下限(0.38)和上限(0.52)来量化由于个体间变异引起的最大CMRO2计算误差。对于皮层GM，Hct＝0.47和两个界限病例之间的最大CMRO2 差异为2.45±0.89(n＝11)。另一个局限是基于QSM的CMRO2需要测量CBV，如同基于R2^*的方法。我们MLV的可行性研究通过CBF来估算CBV，假设来自标准PET 数据的经验公式，并假设静脉与总CBV之间的比率为0.77，如以前的R2^*和QSM研究一样。尽管这两个假设都与健康人研究的文献一致，但在某些患者中可能无效。此外，由ASL测量CBF的准确性可能受到WM或大静脉中的低信噪比和长动脉传输时间的影响。我们的数值仿真表明，CMRO2计算会受到低信噪比CBF(WM为4.4％误差， SNR 4.5)和CBF值误差(WM和GM均为10％误差，存在10％低估)的影响。这可以通过使用动态对比增强(DCE)或诸如血管空间占用(VASO)的非对比增强MRI技术，基于多期ASL的CBVa或定量BOLD方法来直接测量CBV来解决。另一个局限性是使用直窦的OEF来计算全局CMRO2约束。颈静脉对估计全局CMRO2更为理想，但对于同一视野和/或线圈内的图像更为复杂。此外，脑外QSM成像仍然具有挑战性，因为需要进行额外的相位预处理，例如水脂分离。因此，假定直窦和颈静脉具有相似的氧合状态，使用直窦进行测量，因为使用¹⁵O PET时健康受试者的全脑OEF显示均匀。由于在采集过程中采用流速补偿以及场强的自适应二次拟合(Xu B et.al., MRM 2014)，QSM的流速效应被最小化。最后，本研究中的受试者年轻且健康，该技术在一般患者群体中的验证和适用性仍有待进一步研究。

总之，我们的研究证明了基于QSM和CBF使用MLV进行无需血管干预的无创性CMRO2成像的可行性。在健康受试者中的体内成像表明，用所提出的MLV方法获得的CMRO2和OEF图与使用基于血管干预的方法获得的CMRO2和OEF图一致性较好。

8.心脏的QSM

在本节中，我们使用导航门控数据采集来呈现心脏磁化率源的QSM。

8.1.介绍

混合静脉血氧饱和度(SvO2)测量通过血管系统循环后返回至右心室腔的氧合血红蛋白的百分比。它反映了全身氧运送和全身氧消耗之间的平衡，使其成为广泛用于管理危重病人的心肺功能的重要测量方法。目前，SvO2需要靠侵入性心脏导管插入术来测量。因此，一种使用非侵入性MRI来测量SvO2的方法可能对改善患者护理有一定的价值。

定量磁化率成像(QSM)是一种能够量化抗磁性和顺磁性材料之间差异的MRI成像技术。由于铁中不成对电子的数量不同，氧合血红蛋白呈反磁性，脱氧血红蛋白呈顺磁性。QSM能利用这点来量化身体内血液的氧饱和度水平。然而，活体QSM由于三个主要挑战尚未应用于心脏。这三个挑战是：心脏和呼吸运动，胸部和心脏中的脂肪化学位移，以及背景场。在这里，我们将描述克服这些挑战以成像活体心脏QSM 的方法，并用它来量化左心室腔血液和右心室腔血液中氧饱和度水平之间的相对差异。

8.2。材料和方法

心脏和呼吸运动：

心脏和呼吸运动一直是心脏MRI的挑战，对于心脏QSM更是如此，因为QSM数据采集需要长的射频重复时间(TR)。为了抑制心脏运动的影响，我们使用了ECG触发序列，并且将舒张期内的采集时间限制在200ms以内。为了抑制呼吸运动，我们要求志愿者在扫描过程中进行屏气。志愿者在最大呼气时保持屏气，以确保每个切片在同一呼吸位置成像。

我们还开发了一种通过导航门控能够在自由呼吸的情况下采集数据的方法。在1.5T扫描仪(美国GE)上开发了用ECG触发的导航门控多回波3D FGRE序列来采集活体心脏QSM的数据。使用铅笔光束导航回波来检测横膈膜位置，并且使用高效的2 相位相序自动窗口选择(PAWS²)门控算法来实时控制基于横膈膜位置的数据采集。在每次心跳都会播放两次导航回声，一次是在采集之前，另一次是在采集之后。如果由这两个导航回波检测到的两个横膈膜位置之间的差异大于2mm，那么在该心跳中获取的数据会被拒绝。在健康志愿者心脏上，我们对比了2DBH和3DNAV方法获得的 QSM。

2DBH方法的扫描参数如下：1^stTE/TE/#TE/TR/BW＝1.5ms/2.2ms/5 /15ms/±62.5kHz，采集分辨率＝192x192，75％相位FOV，20片，切片厚度5mm，每心跳10次TR。3DNAV的扫描参数和2DBH的相同，只是用0.75NEX和0.75部分切片编码获得了24个切片。

我们招募了五名健康志愿者，每位志愿者都同意IRB批准的研究协议。我们在五名健康志愿者上用2DBH和3DNAV序列进行扫描。在2DBH序列扫描中他们被要求在最大呼气时执行屏气；在3DNAV序列扫描时他们能自由呼吸。

我们使用了化学位移更新方法³来处理复数数据来获得整体场，该方法通过基于图形切割的相位解缠和脂肪水分离方法⁴的输出进行初始化。然后使用整体场反演的方法重建QSM。

一个成功的偶极子反演需要一个连续的三维场作为输入，所以为了确保这个条件得到满足，采集每个切片时志愿者必须保持一致的屏气位置。为了模拟切片错位的后果，我们回顾性地从一组数据中的中间切片上进行平面移动(x和y方向均为 10mm)，然后计算QSM。

脂肪化学位移：

体内的脂肪会产生人造相移或化学位移。如果在重建QSM之前未消除这些化学位移，则产生的QSM将具有严重的条纹伪影。为了准确地去除脂肪化学位移，我们首先采用了幅度图像引导图切割方法，用于相位解缠和化学位移的初始估计。联合场图估计和水脂分离问题可以表示为以水和脂肪作为组织化学成分的复杂MRI数据的以下非线性最小二乘拟合问题。

其中s(TE_n)是在某个TE_n回波时间测得的MR信号,是衰减率，f_F是推测的脂肪化学位移(在1.5T扫描仪f_F≈-220HZ，在3T扫描仪f_F≈-440HZ)， b是场图(或场不均匀性)，ρ_W是水分量，ρ_F是脂肪分量。上述目标函数是一个具有多个局部最小值的非凸函数。这种最小化算法对b的初始估计特别敏感。水和脂肪产生的复合信号可以表示为其中是s(TE_n)的数量级，而是的相位。为简单起见，我们假设ρ_W和ρ_F是实数。一个有效的可以用一个单一的物种拟合来估计，这里我们研究和真正的f_s之间的关系。我们旨在最小化以下残差函数，它使用单个物种模型来拟合两个物种的信号

对于这个非线性最小二乘拟合问题，计算解析解很困难，所以我们在这里提供一个近似的解析解。简单地说，我们使用泰勒展开到近似b的一阶，并将残差函数的导数设为零来得到近似解

其中)，k是整数，并且k/ΔTE卷绕由有限的动态范围引起。真正的场图可以写成

假设实际场图f_s在空间上是连续的，并且模型误差比相位卷绕或化学位移小，那么一组合适的k和m可以解决f中的不连续性

其中表示梯度算子。这种关于k和m的联合优化可以使用带有条件跳跃移动的图形切割来有效解决。场图中的模型误差受ψ(TE_n，ρ_W，ρ_F)限制。从 ψ(TE_n，ρ_W，ρ_F)的近似解中，我们发现更大的最后回声TE_N导致较小的拟合误差。另外，如果ρ_W＝0或ρ_F＝0，则ψ(TE_n，ρ_W，ρ_F)＝0。也就是说，如果体素由纯水或纯脂肪组成，那则拟合的模型误差为零。在单个物种拟合中引入的真实模型误差将在方法部分的数值模拟详细描述。真实模型误差随着TE_N的增大而减小；当 TE_N大于17.6ms时，模型误差ψ(TE_n，ρ_W，ρ_F)小于15Hz，平均误差大小为2.67Hz。因此，上述同时相位解缠和化学位移结果b可以用作IDEAL的良好初始猜测。

以上结果用于上述公式43的信号拟合，以获得水图和脂肪图。这种基于图形切割的同步相位解缠和化学位移估定来初始化迭代化学位移更新方法以消除来自脂肪的化学位移。

背景场和偶极子反演：

由于肺T2很短，在GRE图像中来自肺的信号通常是不可靠的。因此，心脏的心肌紧邻没有可靠信号的区域，即具有产生背景场的磁化率源的区域。用于去除或抑制背景场的现有技术方法通常通过估计过高或侵蚀而在感兴趣区边缘产生误差。为了确保在最终的QSM中精确抑制背景场的影响，我们使用了整体场反演方法：

其中，χ是整个视野中的磁化率图，b是由上述脂肪水分解估计的磁化率产生的整体场，包括来自感兴趣组织和背景的磁化率源，d是偶极子核，是梯度算子 λ是正则化参数，M_G是边缘掩码，并且w是SNR加权。P是提高算法收敛速度的预处理器。预处理器是一个二值掩码：在感兴趣区内部是1，而在外部是大于1。我们从经验中选择P的外部值为30。

8.3.结果

五名志愿者中的三名能够成功地在同一个呼吸位置上屏气，从而成功地用2DBH方法获得QSM。相比之下，3DNAV方法成功地在所有5名志愿者上获得QSM。而且在这五位志愿者里从3DNAV方法的得到的QSM的图像质量始终优于从2DBH方法上得到的QSM图像质量。

如图42所示，当所有屏气都保持在同一个呼吸位置，来自2DBH方法的QSM与来自3DNAV方法的QSM具有类似的右心室腔对左心室腔对比度。当由于屏气位置不一致时(如箭头所示)，由此产生的QSM受到重大伪影的影响，从而变得不可靠。因为2DBH方法的信噪比比3DNAV方法的低，所以来自2DBH的QSM看起来比来自3DNAV 的QSM拥有更多的噪声。

从能够在同一个呼吸位置上屏气的三名志愿者的2DBH QSM上测量的右心室腔血液和左心室腔血液之间的平均磁化率差异为270±18ppb，这相当于79.6±1.1％的氧合水平。从五名志愿者的3DNAV QSM上测量的右心室腔血液和左心室腔血液之间的平均磁化率差异为262±13ppb，其为80±1％SvO2。

8.4。结论

我们的初步数据演示了2DBH方法在活体心脏QSM成像的实际运用性有限，因为其很容易受到不一致屏气位置的影响。为了理解这个漏洞，我们在一组成功地用 2DBH方法获取的数据上模拟单个切片错位，然后在最后的QSM中发现了显着增加的伪像。因此，哪怕只有少量单片错位可能在QSM中产生严重伪像。3DNAV方法消除了屏气的需要，并且由此产生的心脏QSM比2DBH方法来的心脏QSM一致性更好，所以3DNAV方法在临床环境中是一个更可行的选择。

综上所述，采用自由呼吸导航门控进行3D采集是一种有前途的活体心脏QSM数据采集的方法。

9.设备系统的实现

某些情况下，上述的一种或多种定量磁化率成像技术可以使用图43所示的程序800实现。例如，程序800可用于被试组织的磁化率成像，如病人(或病人的某个部位)或实验样品(如水模、一种材料或组织样本等)。

在一些实现过程中，程序800可以将被试的MR信号转换为组织的磁化率图进而定量描绘被试的结构和组成。例如，某些情况下程序800可以获得被试的MR数据，进一步数据处理可得到被试的定量磁化率图像。采用这种定量磁化率图，被试的一幅或多幅基于磁化率的图像得以产生并呈现给用户，这些图像可以用于临床诊断或其他实验目的，例如研究被试的结构、组成和功能，或者用于诊断各种疾病，又或者基于部分图像进行各种疾病的治疗。相比其他磁化率成像技术，程序800能够生成更高品质和准确性的磁化率图像，因此程序800的执行可以提高用户对被试的结构、组织和功能的认识，同样可以用来提高医疗诊断或实验分析的准确性。

程序800首先采集被试的磁共振(MR)数据(步骤810)。MR数据采集对象可以是病人(例如，病人的全部组织或某部分组织，如病人身体的特定部位)。另外，也可以采集一个或多个样品的MR数据(例如，水模，一种或多种材料，一种或多种组织的样本或其他对象)。

在具体实施中，MR数据可以通过使用MRI扫描仪扫描并选择合适的脉冲序列获得。例如：使用梯度回波序列采集单回波或多回波(例如：二、三、四、五甚至更多回波)MR数据。扫描时可以使用各种不同的扫描参数，如：使用3T磁共振扫描仪 (例如，美国GE EXCITEHD MR扫描仪)采用三维多回波扰相梯度回波序列进行扫描，脉冲序列可采用不同的回波时间(TE)(例如：4回波采集，回波时间分别为3.8ms、 4.3ms、6.3ms和16.3ms)进行间隔采集，其他成像参数如下：重复时间(TR)＝22ms；体素大小＝0.5*0.5*0.5mm³；矩阵大小＝256*64*64；带宽(BW)＝±31.25kHz；翻转角＝15°。又如，在3T磁共振扫描仪上使用三维多回波扰相梯度回波序列采用如下参数进行扫描：回波时间(TE)＝5ms；回波间隔(ΔTE)＝5ms；回波个数(#TE)＝8；重复时间(TR)＝50ms；翻转角(FA)＝20°；带宽(BW)＝±62.50kHz；视野大小(FOV) ＝21.8×24×12.8cm³；矩阵大小＝232×256×64。上述给出的扫描序列和序列参数的例子仅用以说明使用。在具体实施过程中，还可以使用很多其他序列和参数，这取决于各种不同因素(例如，扫描区域的大小、已知或假设的被试的磁化率值范围、扫描时间考虑、设备自身的限制等等)。

获取被试MR数据后，程序800可基于MR数据确定磁场大小(步骤820)。如上所述，MRI的相位值与被试体内磁化率分布导致的磁场和组织成分的化学位移息息相关。

磁场可通过对组织各成分共同决定的复数MR数据进行拟合得到，其中每个信号组分包括表示信号衰减的负实数部分和代表依赖于磁场和化学位移的相位的虚数部分(步骤830)。这一复数信号模型的拟合采用迭代数值优化算法进行。数值优化的稳健初始值评估可以使用一个简化的单一模型和图像分割分区平滑模型解决化学位移和相位缠绕问题。

利用MR数据确定磁场之后，程序800进一步建立磁场和磁化率之间的关系(步骤840)。如上所述，磁场和磁化率之间的关系可以表示为一个给定的位置的磁场与该位置的磁化率之间的关系。在某些情况下，这种关系可以用积分形式或微分形式表达。微分形式中，给定位置的磁场与该位置的磁化率之间的关系可以概括为一个方程，磁场经拉普拉斯运算等于三分之一的磁化率的拉普拉斯运算减去磁化率的二阶导数。

在确定了磁场和磁化率之间的关系后，程序800然后通过先验知识确定组织的磁化率分布(步骤850)。一种评估组织磁化率分布的方法是利用组织的R2^*值。一个区域的R2^*值越高，表示该区域的磁化率值也越高(例如出血)，此可作为预处理值以加速优化算法的收敛。预处理的一个例子是将全部图像区域分为高磁化率区(包括出血、空气区域和背景)和正常组织磁化率区域。这种预处理减小了磁化率搜索范围，从而加速收敛。另一种估计组织磁化率分布的方法是使用低(接近零)R2^*区域确定纯水区域，如大脑脑室内脑脊液、主动脉的含氧动脉血、腹壁纯脂肪区域。由于水和脂肪的磁化率已知，因此可被用作参考零值(水)，然后使用正则化最小方差算法得到组织的绝对磁化率。

获得磁化率初始信息和信号的噪声特性后，程序800继续通过之前的信息和噪声特性评估被试或其部分区域的磁化率分布(步骤860)。如上所述，评估被试的磁化率分布包括确定对应的磁化率分布的价值函数、磁场、感兴趣区的掩模。价值函数包括一个基于数据的噪声特性的数据保真项，用积分或微分形式表达磁场和组织磁化特性的关系；一个正则项表示先验信息。被试的磁化率分布评估可以通过确定一个特定的磁化率分布得到，这一特定分布通过最小化一个或多个价值函数来确定。如上所述，在某些情况下，数值上是可以确定的。

估算出被试磁化率分布后，程序800然后生生一幅或多幅磁化率分布图像(步骤870)。这些图像可以在合适的电子显示设备上显示(如液晶显示器，LED显示器，CRT 显示器，或其他可以显示图像的设备)或通过其他合适的物理介质进行显示(例如，印刷、蚀刻、绘画、或呈现在一张纸、塑料或其他材料上)。

在某些情况下，通过将一种颜色映射到一个特定的磁化率值或特定范围的磁化率值，可以得到直观的彩色磁化率分布图像。因此，亦可生成二维或三维的图像，其中图像的每个像素或体素对应于被试特定的空间位置，该像素或体素的颜色表示该位置的磁化率值。

一般情况下，运用色标表示一个范围内的磁化率值，色标包括彩色梯度或灰色梯度(如灰度)。在某些情况下，色标包括渐变的色彩或渐变的灰色色调，色标的一端对应特定磁化率值范围的最低值(例如，任意选择一个窗口的值)，另一端对应该磁化率范围的最高值。例如，灰度色阶，包括纯白色和纯黑色之间灰度，纯白色可以用来表示在一个特定的窗口磁化率最高值，纯黑色表示该窗口的磁化率最低值。彩色/灰度色阶和磁化率之间的其他对应关系也可以实现。例如，在某些情况下，介于纯白色和纯黑色之间的灰度梯度，纯白色可以表示在一个特定的任意窗口中较低磁化率值，纯黑色可以表示该窗口的磁化率最高值。虽然上述例子是以灰度色调为背景，但实际中彩色色阶和磁化率值之间的也可以建立类似的对应关系。

磁化率和颜色之间可以线性映射(例如，磁化率值的绝对变化对应了颜色的成比例变化)、对数映射(例如，磁化率指数变化对应了颜色的线性变化)、或任何其他映射均可。以上仅仅给出了一些说明性的实例用于表示彩色色阶与磁化率值之间的映射。在具体实现中，可以根据实际情况建立其他的映射。

在某些情况下，根据先验知识可以完成磁化率分布估计的正则化。例如，使用以上所述的一种或多种技术来正则化磁化率分布的估算。在某些情况下，系统评估磁化率分布所依据的先验知识包括与被试相关的确定边界条件。

上述技术使用计算机系统进行实现。图43是使用的计算机系统900的结构图，如，执行程序800。在一些实施中，计算机系统900可通过通信连接到另一个计算机系统(例如，另一个计算机系统900)，以接收数据(例如，MRI数据集)并使用上述一种或多种上述技术进行数据分析。

系统900包括处理器910，内存920，存储设备930，和输入/输出设备940，每个组件910，920，930和940之间可以互联。例如，使用系统总线950。处理器910 能够执行系统900内指令。有些情况下，处理器910是单线程的处理器；另一些情况下，处理器910是多线程处理器；还有些时候，处理器910是量子计算机。处理器910能够处理存储器920或存储设备930中的指令，其能执行上述一种或多种技术。

内存920将信息存储在系统900内。一些情况下，存储器920是可读计算机媒介；另外一些情况中，存储器920是易失性内存；还有些情况中，存储器920是非易失性存储器单元。

存储设备930能够为系统900提供大容量存储。有些情况中，存储装置930是非暂时性的可读计算机介质。在实际情况中，存储设备930可能包含硬盘装置、光盘装置、固态硬盘、闪存驱动器、磁带或其他大容量存储设备。有时，存储装置930 可能是一个云存储设备，例如，一个逻辑存储设备包括分布在一个网络并可使用网络访问的多个物理存储设备。在有些应用中，存储设备可以存储长期数据。输入/输出设备940为系统900提供输入/输出操作。在一些实例中，输入/输出设备940包括一个或多个网络接口，例如，一个以太网卡、串行通信设备；又如，一个RS-232 端口、无线接口设备；再如，一个802.11卡、一个3G无线上网卡、4G无线调制解调器等。网络接口设备允许系统900进行传输操作，例如，发送和接收数据。在一些实施方案中，输入/输出装置包括驱动装置用来接收输入数据并发送输出数据到其他输入/输出设备，如键盘、鼠标、打印机、传感器(例如，测量元件或系统相关性能的传感器，测量环境相关属性的传感器，或者其他类型的传感器)和显示装置960。另外，还可使用移动计算设备、移动通信设备和其他设备。

一个计算机系统可通过一个或多个处理装置执行指令完成上述所述的处理和功能，例如，存储、维护和展示结果。这些指令包括，解释说明如脚本指令、执行代码或其他存储在计算机可读介质中的指令。一个计算机系统可以实现分布式网络，如：服务器群，或一组广泛分布的服务器；或包括多个分布式设备的单一的虚拟设备，其彼此之间可协调运行。例如，一个设备可以控制其他设备，或配套设备在一系列协同的规则或协议下运行，或设备与另一个设备协调运行。多个分布式设备的协调运行就好像一个单一设备的运行。

图43中已经描述了一个处理系统实例，主体的执行和上述功能的运行还可以在其他类型的数字电子电路或计算机软件、固件或硬件上实现，包括本规范描述的结构及其结构等价物或它们的组合。本规范中所描述的主体的实现，例如执行一个或多个上述过程，可以用一个或多个计算机程序产品来实现，如在一个有形的程序载体上进行计算机程序指令模块的编码，一个可读计算机介质用于执行或控制处理系统。可读计算机介质可以是一个可读的机械存储设备、可读的存储介质或者存储设备、影响可读传播信号的组合设备、或它们中的一个或多个的组合。

“处理模块”涵盖了所有的仪器、设备、以及用于处理数据的机器，包括一个可编程处理器、一台电脑、或多个处理器或计算机。处理模块除了包括硬件外，还包括就某一问题计算程序而创建执行环境的代码，例如构成处理器的固件协议栈的代码、数据库管理系统、操作系统、或这些中一个或多个的结合。

一个计算机程序(也叫程序、软件、应用软件、脚本、执行逻辑或代码)可以用任何形式的编程语言编写，包括编译或解释性语言，声明或程序语言。它可以任何形式实现，包括作为一个独立的程序或模块、组件、子程序或其他适用于计算环境的单元。一个计算机程序并不一定对应于某一个特定文件系统，它可以存储在一个保存了其他程序或数据的文件中(例如，存储在标记语言文档的一个或多个脚本中)，也可存储在一个单一程序文件或多个文件中(例如，文件存储了一个或多个模块，子程序或部分代码)。计算机程序可以在一台计算机上执行，也可以在多台计算机上执行，这些计算机位于同一个网址或者分布在通过通信网络相互连接的多个网址上。

适合存储程序指令和数据的计算机可读介质包括各种形式的非易失性或易失性存储器、媒体和存储设备，包括半导体存储设备(如EPROM、EEPROM)、闪存设备、磁盘(例如，内部硬盘、可移动磁盘或磁带)、磁光盘、CD-ROM和DVD-ROM光盘。处理器和内存可以附属或并入专用逻辑电路。有时服务器是一台通用计算机，有时它是一个定制的专用电子设备，有时它是这些东西的一个组合。它的实现包括后端组件，例如，一个数据服务器和中间组件；一个应用程序服务器或前端组件；一台具有图形用户界面或Web浏览器的客户端计算机，通过它用户可以与主体进行互动；后端或前端组件、中间件的任意组合。该系统的组件可以通过任何形式或数字数据媒介通信，如通信网络。通信网络包括局域网(“LAN”)和广域网(“WAN”)，例如，互联网。

上述分开执行的某些功能可以组合成单一执行。反言之，一个包含很多功能的单次执行也可分别执行或在任何子组合中实现。

如上所述的操作顺序可以改变。在某些情况下，多任务并行处理可能更有利，但上述系统组件的分离实现并不是必须。

以上叙述了本发明大量的实例。然而，在不脱离主旨的前提下，可以进行各种修改。因此，其它实现方法包含在下述申明范围内。

Claims

1.一种用于给一个物体生成一个或多个图像的方法，所述方法包括：

基于最小化成本函数来确定一个磁化率分布；和

2.权利要求1的方法，其中组织选自：脑脊髓液，血液，软组织，肌肉或脂肪。

3.权利要求1的方法，其中特征是低R2^*值,感兴趣组织区域是通过应用R2^*阈值和连通性来确定的。

4.权利要求1的方法，成本函数最小化涉及预处理。

5.权利要求1的方法，其中物体包括大脑的皮层，壳核，苍白球，红核，黑质和脑丘脑底核中的一个或多个，或腹部肝脏，或胸部心脏，或身体骨骼中，或病变组织的病灶。

6.权利要求1的方法，其中组织成分由脂肪和水组成，并且使用适当的初始化从复数磁共振成像数据迭代地计算体素中的磁场，水分和脂肪分数。

7.权利要求6的方法，其中迭代计算初始化使用基于图形切割的同时相位解缠和化学位移确定或使用同相回波。

8.权利要求1的方法，其中定量生物材料图得以生成，通过校正脂肪对磁化率分布的贡献来产生铁图，或通过校正铁对R2^*的贡献来产生纤维化图，或通过利用短TE的磁共振成像数据和利用骨特异性R2^*对磁共振成像数据建模来生产骨矿物图。

9.权利要求1的方法，其中通过包括磁共振成像幅度数据对脱氧血红素浓度的依赖性模型来生成组织氧提取分数图。

10.权利要求9的方法，其中通过在类似细胞构筑的体素块之间或者在具有类似信号行为的体素之间施加最小方差来对计算进行正则化。

11.一种用于对物体产生一个或多个图像的系统，所述系统包括：

处理器；

通信地耦合到处理器的图形输出模块；

通信地耦合到处理器的输入模块；

基于最小化成本函数来确定一个磁化率分布；和

12.权利要求11的系统，其中组织选自：脑脊髓液，血液，软组织，肌肉或脂肪。

13.权利要求11的系统，其中特征是低R2^*值,感兴趣组织区域是通过应用R2^*阈值和连通性来确定的。

14.权利要求11的系统，成本函数最小化涉及预处理。

15.权利要求11的系统，其中物体包括大脑的皮层，壳核，苍白球，红核，黑质和脑丘脑底核中的一个或多个，或腹部肝脏，或胸部心脏，或身体骨骼中，或病变组织的病灶。

16.权利要求11的系统，其中组织成分由脂肪和水组成，并且使用适当的初始化从复数磁共振成像数据迭代地计算体素中的磁场，水分和脂肪分数。

17.权利要求16的系统，其中迭代计算初始化使用基于图形切割的同时相位解缠和化学位移确定或使用同相回波。

18.权利要求11的系统，其中定量生物材料图得以生成，通过校正脂肪对磁化率分布的贡献来产生铁图，或通过校正铁对R2^*的贡献来产生纤维化图，或通过利用短TE的磁共振成像数据和利用骨特异性R2^*对磁共振成像数据建模来生产骨矿物图。

19.权利要求11的系统，其中通过包括磁共振成像幅度数据对脱氧血红素浓度的依赖性模型来生成组织氧提取分数图。

20.权利要求19的系统，其中通过在类似细胞构筑的体素块之间或者在具有类似信号行为的体素之间施加最小方差来对计算进行正则化。