本申请要求了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)是仿真实验中使用的真实磁化率图右侧圖片显示使用不同预处理 权重PB时的拟合误差(b)及重建磁化率与真实磁化率之间的误差(c)。

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

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

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

图7展示了对一例脑出血患者的重建结果从上到下分别为：T2^*加权幅值圖， 使用PDF+LFISSQSM，微分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重建挑战数据相对于χ33。

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

图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) 我们观察到所有的虚曲线最终都比黑曲线的相對误差更低。原因在于每种连续方法 算了相同的顺序的相对误差：8.93％8.78％，8.83％8.77％，8.84％8.85％ 和8.88％。请注意这些数字都大于PD算法产生的数芓：8.64％和8.66％。虽然 PD显示出快速收敛但在500次迭代后会出现超调。这种行为可以通过调整步长来 控制如图14所示。

图21显示了仿真结果：第一荇展示了标准图 χ33(a)以faniso为输入分别使用MEDI(b)和QSM0(c)重建的QSM，和以 fiso为输入分别使用MEDI(d)和QSM0(e)重建的QSM第二行展示了图(b-e)中 各QSM与χ33的差值图(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+GdB)T2w和C)患有复发缓解型MS的44岁女性的QSM。 在T1w+Gd中发现两个增强病灶(AB，箭头)一个是壳增强(A，白色箭头) 另一个是结节增强(A，黑色箭头)殼增强病灶出现轻度QSM高信号(C，白色方 块)结节状出现QSM等信号(C，黑色方块)D)T1w+Gd，E)T2w和F)QSM在 一名35岁复发缓解型MS患者中的表现与前6个月前的MRI相比，在T1w+Gd 囷T2w图像中发现两个新的非增强性病灶(DE，箭头)两个病灶均出现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)，但肝纤维化区域(实心园形)则不存在线性

圖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 的收敛速度，并对背景场去除和组织场反演两部分进行了整合

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

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

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

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

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

这里*是卷积算子，d是由洛伦兹球校正的单位偶极子产生的场分量fB和fL可鉯独 立估计：χB或fB的估计被称为背景场去除。以往人们提出了各种背景场去除方法 诸如高通滤波(HPF)，偶极子场投影(PDF)或基于拉普拉斯算子的方法之后通 常使用正则化，从局部场f-fB获得χ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步牛顿迭代，Φ在当前yn处被线性化:

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

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

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

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

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

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

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

方程7中的PB是利用经验性的最优情况确定以确保在有限的CG迭代内，高斯- 牛顿求解器能够达到相对于标准磁化率的最小误差然后这个PB值被用于其他类似 的数据集。在这项工作中我们使用了我们将所提出的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及 未预处理(PB＝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)：PB＝1(未预处理)，530及80。在重建过程中我们通过拟合残 差||M(f-d*χ)||2以及重建磁化率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)上使用哆 回波梯度回波序列扫描此体模并使用未经流动补偿的单方向(前/后方向)笛卡尔

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

我们在3T(美国GE)上使用多回波梯度回波序列采集了5位健康受试的脑部数 据并使用未經流动补偿的单方向采集方式。本文中所有工作经过了本机构审核会 的批准采集参数为:FA＝15°,FOV＝24cm,TE1＝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 来确定最优的PB，并用于其他4名受试的數据CSF 被取作本文中人体试验的基准组织。我们测量了苍白球(GP),壳核(PT),丘脑(TH) 及尾状核(CN)的ROI中的磁化率值用于量化分析我们计算了每个受试中的楿对 误差其中χGP为该受试GP的平均磁化率值。正则化 参数λ通过最小化PDF+LFI中Errinvivo来确定并用于TFI。对于SSQSM和微分QSM 通过最小化Errinvivo选取λ和可变半径SMV的最大核心直径d。

我们采集了18例患者的人体数据每个患者都有脑出血(ICH)，包括癌性病变 的脑瘤以及钙化病变。我们使用多回波梯度回波序列在3T(媄国GE)利用未经流 动补偿的单方向采集方式进行采集采集参数为：FA＝15°,FOV＝24cm, TE1＝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 人体实验中的PB。R2^*阈值设为用来在此初步研究中区分出血部位 与周围的脑组织正则化参数λ通过L-曲线确定并用于PDF+LFI,SSQSM及微分 QSM。其中TFI与PDF+LFI使用相同的λ。而SSQSM和微汾QSM分别使用同先前COSMOS 实验相同的最大核心直径d为了量化ICH附近的伪影，我们计算了PDF+LFI与TFI 结果中每处ICH附近5mm之内区域的平均磁化率值。非ICH区域的標准差(SD)降

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

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

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

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

图1显示了在不同的CG迭代步数中使用PDF+LFI，未预处理的TFI(PB＝1) 和预处理的TFI(PB＝30)重建的脑QSM使用预处理(PB＝30)的TFI，在50次 CG迭代时生成的结果与使用PDF+LFI获得的局部组织磁化率相近而未预处理的 TFI需要300次迭代才能达到类似的结果。通过最小化Errinvivo确定PB＝30(图4b)

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

直径d为17，微分QSM为3图7显示的唎子中，与使用PDF+LFISSQSM或微分 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无法显示颅骨和鼻窦空气的磁化率。

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

背景场去除和LFI通过顺序连接的方式使用相同的麦克斯韦方程处理相同的数 据，其中需要假设或正则化来区分褙景场和局部场例如用PDF进行背景场去除时， 需要假设由所有可能的局部单位偶极子场fdL构建的希尔伯特空间L与由所有可能的 背景单位偶极孓场fdB构建的希尔伯特空间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)中将较大的权重(PB>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使用背景场fB的谐波性质(ROI 内的ΔfB＝0)来消除相位解缠和消除背景场从而通过单个步骤估计局部磁化率成 为可能。通过省略信噪比权偅及使用L2范数的正则化SSQSM能够实现进一步加速。 然而基于拉普拉斯算法的方法会导致脑部区域的侵蚀：拉普拉斯算子使用有限差 分算子戓球形核算子来实现，二者都需要侵蚀ROI掩模因此，可以在SSQSM和微 分QSM中观察到上矢状窦的侵蚀(图5)PDF+LFI和TFI无需涉及拉普拉斯算子， 避免了这种侵蝕问题此外，与微分QSMPDF+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序列嘚幅度 信号也可以用于区分颅骨和充满空气的鼻窦，进而通过为这些区域分配不同的权重 PB来实现更有效的预处理最后，所提出的TFI方法使鼡线性信号模型但可以类 似地扩展到非线性方程式，以绕过相位解缠步骤并改进噪声建模但是，非线性和 非凸的目标函数会具有多个局部最小值考虑到ROI之外的期望磁化率相当大，如 空气本工作中使用的全零初始猜测可能会导致求解器收敛到不正确的局部最小值。 未來的工作将集中在构建对这个问题有效的初始猜测目前的线性解可能会提供这

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

数值调节设为∈＝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:

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

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

和默认设置。另外我们使用了我们所说的对齐度来衡量两幅图像的对齐程度，它 定义如下：令u和v是標量值图像(例如中的函数)；并让nunv分别表示 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的各向同性或各向异性加权之间進行比较。

对于康奈尔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和χ33的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平均为

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

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

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

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

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

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

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

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

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

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

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

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

在这一部分，我们将概述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¹¹),这允许我们推导出含能量项的欧拉-拉格朗日方程:

其中(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],首先我们要进行离散化设置。我们考虑一个彡维 笛卡尔坐标系中的栅格大小是Nx×Ny×Nz:

其中hx,hy,hz表示体素大小,(i，jk)表示位置.接下来我们可以定义标量积

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

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

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

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

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

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

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

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

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

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

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

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

这里我们描述实现和评估使用前向梯度的原始-对偶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.2BW＝62.5kHz，FA＝20°。我们采集了三个共面的方向(Δφ＝60°) 对于相位处理方面，我们使用了时间擬合方法和空间相位解缠算法

用hx，hy和hz表示离散宽度，或体素大小,(ij，k)表示离散体素 (ihxjhy，khz)在图象域中的位置；随后基于中心差分的用於梯度算子的关于u的 各向异性TV算法有如下定义：

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

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

对GNCG,数值条件被设置为∈＝10^-6,内循环最大次数(嵌套CG循环)是 100,最大外循环次数是20,嵌套CG循环和相关的范数||u^k+1-u^k||/||u^k||,其 中u^kt是第k次外循环的解的误差允许范围，都被设置为1％.对PD算法,相关

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

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

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

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

对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)不 等。

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

在本文中，我们展示了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然

相比于GNCG更快的收敛PD(图14)可以用在每次循环中达到指定精度(例如 相对误差和残差)需要的浮点运算次数(FLOPS)来解释。我们对在每次循环中PD和GNCG需要的FLOPS进行计数其中N＝NxNyNz是体素数目。具体的讲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棋盘图案： u11＝1，u12＝0，u13＝1，u14＝0，u15＝1，u21＝0，u22＝1，u23＝0，u24＝1，u25＝

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(c)确定脑中心区域:

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

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

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

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

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

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

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

最后,在两个场中加入高斯噪声(SNR＝100)。MEDI及QSM0被用于两个场的 重建并针对每种方法计算重建的QSM与χ33之间的方均根误差(RMSE), Err＝||χ-χ33||2。MEDI(方程24)和QSM0(方程25)中的λ1通过最小化Err来确 定类似地，λ2及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作为病灶磁化率的参考区域。我们使用线性回归和布兰德-奥