CN111557663B - 一种基于交叉模态的人脑磁化率张量成像方法 - Google Patents

Abstract

本发明公开了一种基于交叉模态的人脑磁化率张量成像方法，在进行GRE序列的磁共振成像时，通过旋转头部姿态可以获得多个方向的相位数据，由此可推导出组织的磁化率及其各向异性信息，即磁化率张量成像。临床应用中扫描时间和成像空间限制了人脑成像的方向数和方向角，导致磁化率张量成像的抗干扰能力差，是条件数严重病态的逆问题。本发明利用张量谱分解技术，将高分辨率的磁化率成像和低分辨率的扩散成像结合，提出一种基于交叉模态的磁化率张量成像方法，使用在6个不同头部方向采集的相位数据，实现磁化率张量及导出量的准确计算，可用于脑组织微结构特性评定，并且导出的纤维方向图有助于探索脑白质和灰质内部纤维特征的活体无损检测新方法，对脑科学研究和脑部疾病诊断具有重要意义。

Description

一种基于交叉模态的人脑磁化率张量成像方法

技术领域

本发明涉及医学影像处理和分析技术领域，特别是指一种基于交叉模态的人脑磁化率张量成像方法。

背景技术

磁化率是物质的固有属性，能够反映生物组织在外加磁场作用下的磁化程度，物质的磁敏感性越强，磁化率越大。作为一种磁共振成像方法，磁化率成像利用梯度回波数据的相位信息，经过预处理得到局部磁场的变化特性；可以对组织的铁含量、钙化、血氧饱和度等进行有效的量化分析，对脑出血、多发性硬化症及帕金森综合症等脑神经退行性疾病的研究和诊断具有重要意义。磁化率和磁场分布关系在本质上是非局部的，它依赖于磁化率的空间分布和相对于主磁场的方向，从而表现出各向异性。研究表明纤维的髓鞘中圆柱排列的脂质分子是引起磁化率各向异性的主要原因，而基于磁化率的纤维束成像机理是脑白质沿纤维走向具有最小反磁性。与测量水分子扩散速率的扩散张量成像(DiffusionTensor Imaging,DTI)相比，磁化率成像具有高信噪比、高空间分辨率、以及适合于超高场成像的特性。

磁化率张量成像(Susceptibility Tensor Imaging,STI)为活体组织神经纤维示踪提供了一种新方法。Liu首先测量并定量描述了鼠脑白质中枢神经系统中磁化率的各向异性，采用STI实现了鼠脑纤维束的三维重建，随后在活体鼠的大脑、心脏和肾脏展开一系列基于磁化率信息的纤维束示踪研究。同时，研究人员也在探索STI的活体人脑应用，诸如磁化率张量计算、张量导出指标评估，全脑纤维追迹等。磁化率张量用一个2阶对称矩阵X来描述人体组织中磁化率的各向异性。在成像目标定义的坐标系中，测得的局部场图与磁化率之间的关系表达式为：

进行第i个方向的成像时，δB⁽ⁱ⁾表示采集的局部场图，H⁽ⁱ⁾表示主磁场在成像目标坐标系中的向量。空间频率向量k＝[k_x,k_y,k_z]且

F与F^-1分别表示傅里叶变换及其反变换。磁化率张量X描述单个体素内的磁化率信息：

已知采集到N个方向的相位数据，则磁化率张量与局部场的关系可展开成：

定义磁化率张量成像的系统矩阵

即

公式(3)可写成矩阵形式δB＝AX，由推导可知，磁化率张量成像需要通过旋转头部姿态，获得至少6个方向的相位数据。实际上，MRI成像仪的主磁场无法旋转，而头部成像线圈空间狭窄，人脑在其内部可实现的旋转角度很小，通常情况下与主磁场的夹角范围在0～30°，导致系统矩阵A的条件数较大，产生误差放大效应，使得磁化率张量成像的抗干扰能力差，成为一个具有病态条件数的逆问题。

现有的STI研究采用SVD方法求广义逆或直接应用最小二乘LSQR方法求解磁化率张量，为了提高解的稳定性和可靠性，所用数据的成像方向大于10，使得成像时间显著延长。受制于成像时间和患者配合度等问题，临床应用难以实现多方向、大旋转角的磁化率成像。综上所述，由于受到成像时间和成像空间的制约，人脑STI无法像DTI那样实现高角度分辨率的纤维束示踪。另一方面，STI处理流程包含多个复杂的步骤：多方向数据配准，相位解缠绕，去除背景场，建立系统矩阵，求解病态逆问题；使得磁化率张量的求解受到成像噪声和伪影、预处理过程中的偏差，以及图像配准误差等多种因素的影响。因此，活体人脑的磁化率张量成像及以此为基础的纤维束示踪工作至今还很少见，现阶段相关研究还面临着很多的问题和很大的挑战。

为了利用有限的采样方向获得令人信服的纤维束示踪，一些假设和先验信息被作为约束项引入到STI计算。例如，假设纤维是对称圆柱型，将磁化率张量简化为平均磁化率(mean magnetic susceptibility,MMS)和磁化率各向异性(magnetic susceptibilityanisotropy,MSA)两个指标。进一步的研究提出将磁化率各向异性限定在白质区域，并加入平均磁化率的正则项，以加强各向同性与各向异性组织的区分，但是这样可能导致纤维通路被切断。由于方向依赖性也表现在R2*图中，研究人员提出基于互相关的各向异性和联合特征向量的STI方法，结合幅度导出的弛豫张量和相位衍生的磁化张量，改善基于磁化率成像的纤维束追迹效果，在小鼠大脑、肾脏和心脏等部位成像中获得了较为精确的纤维追迹图。

发明内容

本发明的主要目的在于克服磁化率张量成像的方向局限性，减少所需采样方向的数目，实现准确的磁化率各向异性估计，推动磁化率张量成像技术的应用。

本发明采用如下技术方案：

一种基于交叉模态的人脑磁化率张量成像方法，利用张量谱分解技术，将高分辨率的磁化率成像和低分辨率的扩散成像相结合，采用从扩散张量成像提取的脑白质体素内纤维方向构造磁化率张量求解模型的正则化约束项；其中，基于交叉模态的磁化率张量成像模型的数学表达式为：

式中，磁化量张量X是矩阵形式为

其中χ_ij表示沿着不同坐标方向的磁化率数值，下角标[1,2,3]分别对应三维成像空间的坐标轴方向[x,y,z]，第一项是数据保真项，第二项用于约束脑区外体素的磁化率值趋于零，第三项是基于张量谱分解的正则化约束，δB＝[δB⁽¹⁾(k),...,δB⁽ⁱ⁾(k),...,δB^(N)(k)]^T,i＝1,…N且N＝6，第i个方向的相位数据中坐标向量k＝[k_x,k_y,k_z]的空间点的局部场图表示为δB⁽ⁱ⁾(k)，F和F^-1分别是傅里叶变换和傅里叶反变换算子，矩阵A是系统矩阵，M₁为二进制模板，且有

M₂是脑白质模板，α、β是标量参数，用于调节两个约束项的权重。

根据张量谱分解技术，所述磁化率张量X可以分解为矩阵Q、Q^T和Λ的乘积，即X＝QΛQ^T，其中，矩阵Q＝[υ_1STI,υ_2STI,υ_3STI]，列向量υ_1STI、υ_2STI和υ_3STI是从所述磁化率张量导出的3个特征向量，分别对应特征值λ_1STI、λ_2STI和λ_3STI，根据扩散张量成像DTI导出的纤维方向信息约束磁化率张量的求解，得出基于主特征向量约束的磁化率张量谱分解正则项

与矩阵Q相比，矩阵Q_S的第一个列向量采用υ_1DTI替换υ_1STI，

的矩阵展开形式如下

其中υ_1DTI是由扩散张量导出的主特征向量；υ_1STI、υ_2STI和υ_3STI是从磁化率张量导出的3个特征向量，所述矩阵Λ是由所述磁化率张量X的3个特征值λ_1STI、λ_2STI和λ_3STI组成的对角矩阵。

所述矩阵

其矩阵形式为

其中，系数

表示在成像坐标系下，不同头部方向的偶极核；上角标(i)是方向标记；下角标m,n取值为[1,2,3]，对应坐标轴方向[x,y,z]。

磁化率张量X导出的3个特征值遵循λ_1STI>λ_2STI>λ_3STI的关系，λ_1STI所对应的特征向量υ_1STI是磁化率张量成像的主特征向量，所述υ_1DTI的方向用于反映纤维走向。υ_1DTI是扩散张量成像的主特征向量，所述υ_1DTI的方向用于描述纤维走向。

由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：

(1)本发明依据STI和DTI的主特征向量方向与纤维束走向一致的理论，采用张量谱分解技术创建磁化率张量求解模型的正则化项，即利用扩散成像估计出的纤维走向对脑白质磁化率张量的主特征向量进行约束，利用有限的6个方向相位数据，实现全脑磁化率张量及其各项异性指标的准确计算，有效克服磁化率张量成像的方向数和方向角局限性，缩短了成像采集时间，有助于临床应用的推广。

(2)、STI和DTI具有不同的物理机制，与常规DTI成像较低的信噪比和空间分辨率相比，本发明将高空间分辨率的磁化率张量成像和低空间分辨率的扩散张量成像结合，以磁化率数据为主体，引入扩散信息进行局部约束，有效抑制导出图像中的伪影，保护组织结构对比度。同时能够区分出低分辨率扩散成像中无法识别的一些结构和特征，对基于影像学的大脑微结构研究具有重要意义。

附图说明

图1.基于交叉模态的人脑磁化率张量成像方法的处理流程图；

图2.不同方向的局部场图(横断面)，由多方向场图数据重建出的定量磁化率图COSMOS(横断面和矢状面)，和脑白质不同ROI区域的局部场数值的方向差异性；

图3.磁化率张量的计算结果及张量分量的差值图；χ₁₁、χ₂₂、χ₃₃是二阶张量主对角线上的3个元素；

图4.磁化率张量的导出量图与扩散张量导出量的对比；

图5.磁化率张量的主特征向量的方向彩色编码图，即人脑内神经纤维束走向图。

具体实施方式

以下通过具体实施方式对本发明作进一步的描述。

本发明的技术方案：

一种基于交叉模态的人脑磁化率张量成像方法，利用张量谱分解技术，将高分辨率的磁化率成像和低分辨率的扩散成像相结合，采用从扩散张量成像(Diffusion TensorImaging,DTI)提取的脑白质体素内纤维方向构造磁化率张量求解模型的正则化约束项，克服磁化率张量成像(Susceptibility Tensor Imaging,STI)固有的方向数和方向角局限性，实现应用6个头部方向采集的磁共振相位数据，准确获得磁化率张量各分量的合理数值解，有效抑制图像伪影，保护组织结构对比度。基于交叉模态的磁化率张量成像模型的数学表达式为：

式中，二阶张量X∈R^3×3是一个包含6个分量的对称矩阵，用于描述成像坐标系中体素内组织在不同方向的磁化率信息，其矩阵形式为

其中χ_ij表示沿着不同坐标方向的磁化率数值，下角标[1,2,3]分别对应三维成像空间的坐标轴方向[x,y,z]。本发明的成像模型中，第一项是数据保真项，第二项用于约束脑区外体素的磁化率值趋于零，第三项是基于张量谱分解的正则化约束。δB＝[δB⁽¹⁾(k),...,δB⁽ⁱ⁾(k),...,δB^(N)(k)]^T,i＝1,…N且N＝6，第i个方向的相位数据中坐标向量k＝[k_x,k_y,k_z]的空间点的局部场图表示为δB⁽ⁱ⁾(k)，F和F^-1分别是傅里叶变换和傅里叶反变换算子，A是系统矩阵。M₁为二进制模板，大脑内部区域体素的取值为1，外部区域体素取值为0，且有

M₂是脑白质模板，白质组织对应体素的值设置为1，否则为0。α、β是标量参数，用于调节两个约束项的权重，取值为α＝4、β＝1。

据张量谱分解技术，磁化率张量X可以分解为矩阵Q、Q^T和Λ的乘积，即X＝QΛQ^T，矩阵Λ是由磁化率张量X的3个特征值λ_1STI、λ_2STI和λ_3STI组成的对角矩阵，λ_1STI>λ_2STI>λ_3STI。矩阵Q＝[υ_1STI,v_2STI,v_3STI]∈R^3×3，列向量v_1STI、v_2STI和υ_3STI是从磁化率张量导出的3个特征向量，分别对应特征值λ_1STI、λ_2STI和λ_3STI，其中λ_1STI所对应的特征向量υ_1STI是磁化率张量成像的主特征向量。前期研究表明，虽然STI和DTI的生物物理机理不同，然而两者都可用于研究脑白质的微结构，并且扩散张量和磁化率张量的主特征向量理论上是与纤维束走向一致的，因此对于同个体素有υ_1STI≈υ_1DTI，其中υ_1DTI是由扩散张量导出的主特征向量。在此基础上，本发明采用扩散张量成像DTI导出的纤维方向信息约束磁化率张量的求解，设计出基于主特征向量约束的磁化率张量谱分解正则项

与矩阵Q相比，构造矩阵Q_S中的第一个列向量采用υ_1DTI替换υ_1STI，

的矩阵展开形式如下

基于张量谱分解的正则化约束项中，扩散成像估计出的纤维走向只作用于约束白质区域的磁化率张量计算，灰质及其他区域的磁化率张量计算不受扩散成像信息的约束。由此可见，本发明提出的磁化率张量成像模型是以磁化率成像信息为主体，通过融合脑白质的扩散成像信息，提高人脑磁化率张量成像的准确性。采用基于交叉模态的人脑磁化率张量成像方法，获得的磁化率信息及由此导出的各向异性评价指标具有良好的空间分辨率，能够用于辨识低分辨率DTI中无法识别的组织边界和纹理特征，可实现脑灰质和白质的组织微结构特性分析。

系统矩阵

其矩阵形式为

其中，系数

表示在成像坐标系下，不同头部方向的偶极核；上角标(i)是方向标记；下角标m,n取值为[1,2,3]，对应坐标轴方向[x,y,z]。

的详细计算方法在本领域内相关文献中可查。

以脑脊液的磁化率为参照，脑白质是抗磁性的，磁化率值为负；脑灰质是顺磁性的，磁化率值为正。磁化率的各向异性主要是由包裹神经纤维的髓鞘引起的，沿着纤维方向最具顺磁性。从磁化率张量X导出的3个特征值遵循λ_1STI>λ_2STI>λ_3STI的关系，λ_1STI所对应的特征向量υ_1STI是磁化率张量成像的主特征向量，其方向可用于反映纤维走向。的各向异性是由于水分子运动受到轴突对称性约束，水分子的扩散速率沿着纤维走向最大，υ_1DTI是扩散张量成像的主特征向量，其方向可用于描述纤维走向。

下面对本发明提出的方法进行实验验证。

采用三维流动补偿梯度回波GRE序列，在7T Philips Healthcare MRI上采集6个方向的相位图和幅值图，主要成像参数有：FOV＝224×224×100mm³,matrix size＝384×384×100,flip angle＝15°,TR/TE1/ΔTE＝27/4/4ms，回波个数6，头部相对于主磁场的旋转角度0°,6.6°,17.1°,20.8°,22.1°,26.3°，每个方向的扫描时间约为6min。采用单次激发平面回波成像EPI序列，在3T Philips Medical Systems上采集32个扩散梯度编码方向的扩散加权成像，主要成像参数有：FOV＝212×212×143mm³,matrix size＝256×256×65,TR/TE＝6800/67ms,b-value＝700s/mm²,采集2组无扩散加权成像，2组扩散加权成像，数据平均后以提高图像信噪比。

如图1所示，STI的处理流程包含多个步骤：制作去脑壳模板，绘制脑部感兴趣区域，图像配准，计算扩散主特征向量，建立磁化率成像的系统矩阵等。使用FSL BET软件进行去脑壳处理，利用头部正位方向的幅值图制作脑模版。配准使用FSL FLIRT工具包，选用对比度较好的第三个回波数据，将其他方向的数据配准到头正位方向，选用12自由度的仿射变换。将获得的配准矩阵应用于多回波数据的实部和虚部，随后从配准后的复数信号中提取幅值和相位信息。DWI图像配准到头部正位方向的GRE幅值图，使得两种模态的数据具有相同的matrix size。使用DSI Studio软件计算扩散张量，并且生成扩散成像导出量PEV、各向异性分数FA、平均扩散率MD。

相位数据的预处理包括相位解缠绕、多回波拟合、去除背景场。选用基于拉普拉斯分解的相位解缠方法和区域自适应的背景场去除算法iRSHARP，获得的局部场图表现出显著的磁化率各向异性，如图2所示。头部姿态的方向角标注在图的上方，在成像目标坐标系中，AP表示前后，FH表示上下，RL表示左右，感兴趣区域局部场的方向各向异性如柱状图所示，ROI轮廓标注于COSMOS图，本发明提出的磁化率张量成像方法(附图中标记为DRSTI)，采用最小二乘法求解，迭代收敛条件1×10^-4，正则化参数经优化设置为α＝10和β＝0.1。对比方法有：采用奇异值分解计算伪逆矩阵的STI*，无正则化约束的最小二乘法STI**，进行白质各向异性和MMS形态学约束的MMSR-STI，以上算法使用提出者共享的源代码实现。磁化率张量的导出量定义为MMS＝(λ₁+λ₂+λ₃)/3和MSA＝λ₁-(λ₂+λ₃)/2。

首先进行磁化率张量计算的实验验证，如图3所示。其中，χ₁₁、χ₂₂、χ₃₃是二阶张量主对角线上的3个元素。本发明的方法标记为DRSTI，对比方法分别标记为STI^*，STI^**，MMSR-STI。前4列为不同STI算法求得的张量分量图，第5列是本发明方法的分量差值图。四种STI算法求解的张量中，主对角线分量χ₁₁、χ₂₂、χ₃₃都表现出不同的磁化率对比度。然而，对比算法STI^*和STI^**的结果存在明显的磁化率伪影，特别是在两个半球内囊附近的区域；MMSR-STI方法的结果灰质和白质组织间具有清晰的对比度，但是张量分量图像的磁化率数值偏低。本发明的方法能够有效的抑制图像伪影，各区域的磁化率值分布合理；由张量分量相互之间的差值图像可见，磁化率各向异性符合白质沿着纤维走向最具顺磁性的假设，例如视神经纤维束(深灰色箭头指示)，内囊后肢纤维束(浅灰色箭头指示)和连接左右脑的胼胝体(黑色箭头指示)。其中，胼胝体的纤维方向近似平行χ₁₁，垂直于χ₂₂和χ₃₃，其纤维束的χ₁₁分量比χ₂₂和χ₃₃更具顺磁性。

其次进行磁化率张量导出量平均磁化率MMS和磁化率各向异性MSA的实验验证，如图4所示，其中，包括平均磁化率MMS和磁化率各向异性MSA。最后一列是扩散张量成像DTI的各向异性分数FA图和平均扩散速率MD图。MMS图以脑脊液的MMS值为参照基准，即将脑脊液的磁化率均值取为0。MMSR-STI和本发明方法生成的MMS图中，可以清楚地识别出各解剖结构；而STI*和STI**的MMS图中组织对比度相对较低并且计算误差较大。与DTI中扩散张量的特征值非负性不同，由于顺磁性(灰质)和抗磁性(白质)脑组织的磁化率值存在正负差异，STI中磁化率张量的特征值有正数也有负数，导致MMS图存在一定的平滑性。尽管如此，MMS图的结构分辨率仍高于DTI的平均扩散率MD图。所有MSA图的各向异性表现出一致趋势，并且纤维束密度大的ROI各向异性值较高。受较严重的噪声和伪影影响，STI*和STI**方法的MSA图不准确。MMSR-STI方法将各向异性限定在脑白质中，使得其MSA图只包含部分脑部结构。由图可见，本发明方法的MSA图的结构对比度与DTI成像的FA图有些不同，特别是MSA图中可观察到纹状体中的各向异性，包括尾状核、壳核以及苍白球，通常基于扩散的方法难以获得此类信息。

然后进行磁化率张量的主特征向量PEV的实验验证，如图5所示，其中，本发明的方法标记为DRSTI，对比算法分别标记为STI^*，STI^**，MMSR-STI。最后一列是由扩散张量成像DTI估计的纤维束走向图，方向彩色编码图中绿色对应成像坐标系的anterior-posterior前后向、红色对应right-left左右向、蓝色对应inferior-superior上下向。第一行横断面视图，第二行矢状面视图，与前面的实验结果相一致，由于伪影的影响，STI*和STI**方法生成的PEV图中，组织结构的方向局部连贯性和区域一致性不清晰，无法进行有效的方向判读。MMSR-STI的结果优于STI*和STI**，虽然横断面视图与DTI差异较大，但是矢状面视图中可以观察到与神经解剖结论相符的胼胝体压部和膝部(黑色箭头所示)。与MMSR-STI相比，CMSTI在左脑和右脑之间表现出更好的对称性，白质区域纤维束走向与DTI估计出的方向基本一致。与DTI纤维走向图相比，本发明的纤维方向图中空间分辨率更高，如轴向视图中白色箭头所示，能够更清晰的观察到灰质和白质的边界以及大脑皮层的纹理，有利于更好地识别大脑微结构。

上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。

Claims (3)

1.一种基于交叉模态的人脑磁化率张量成像方法，其特征在于：利用张量谱分解技术，将高分辨率的磁化率成像和低分辨率的扩散成像相结合，采用从扩散张量成像提取的脑白质体素内纤维方向构造磁化率张量求解模型的正则化约束项；其中，基于交叉模态的磁化率张量成像模型的数学表达式为：

式中，磁化量张量X是矩阵形式为

其中χ_ij表示沿着不同坐标方向的磁化率数值，下角标[1,2,3]分别对应三维成像空间的坐标轴方向[x,y,z]，第一项是数据保真项，第二项用于约束脑区外体素的磁化率值趋于零，第三项是基于张量谱分解的正则化约束，δB＝[δB⁽¹⁾(k),...,δB⁽ⁱ⁾(k),...,δB^(N)(k)]^T,i＝1,…N且N＝6，第i个方向的相位数据中坐标向量k＝[k_x,k_y,k_z]的空间点的局部场图表示为δB⁽ⁱ⁾(k)，F和F^-1分别是傅里叶变换和傅里叶反变换算子，矩阵A是系统矩阵，M₁为二进制模板，且有

M₂是脑白质模板，α、β是标量参数，用于调节两个约束项的权重；

根据张量谱分解技术，所述磁化率张量X可以分解为矩阵Q、Q^T和Λ的乘积，即X＝QΛQ^T，其中，矩阵Q＝[υ_1STI,υ_2STI,υ_3STI]，列向量v_1STI、v_2STI和v_3STI是从所述磁化率张量导出的3个特征向量，分别对应特征值λ_1STI、λ_2STI和λ_3STI，根据扩散张量成像DTI导出的纤维方向信息约束磁化率张量的求解，得出基于主特征向量约束的磁化率张量谱分解正则项

与矩阵Q相比，矩阵Q_S的第一个列向量采用υ_1DTI替换υ_1STI，

的矩阵展开形式如下：

其中υ_1DTI是由扩散张量导出的主特征向量；v_1STI、v_2STI和v_3STI是从磁化率张量导出的3个特征向量，所述矩阵Λ是由所述磁化率张量X的3个特征值λ_1STI、λ_2STI和λ_3STI组成的对角矩阵。

2.如权利要求1所述的一种基于交叉模态的人脑磁化率张量成像方法，其特征在于：所述矩阵

其矩阵形式为

其中，系数

表示在成像坐标系下，不同头部方向的偶极核；上角标(i)是方向标记；下角标m,n取值为[1,2,3]，对应坐标轴方向[x,y,z]。

3.如权利要求1所述的一种基于交叉模态的人脑磁化率张量成像方法，其特征在于：所述磁化率张量X导出的3个特征值遵循λ_1STI>λ_2STI>λ_3STI的关系，λ_1STI所对应的特征向量υ_1STI是磁化率张量成像的主特征向量，所述υ_1DTI的方向用于反映纤维走向，υ_1DTI是扩散张量成像的主特征向量，所述υ_1DTI的方向用于描述纤维走向。

Images