CN111481827B - 定量磁化率成像及用于dbs潜在刺激靶区定位的方法 - Google Patents

Abstract

本发明公开了一种定量磁化率成像及用于DBS潜在刺激靶区定位的方法，针对现有的QSM重建技术在磁化率反演时得到的磁化率结果准确度不足，有些位置存在尖锐的伪影及QSM图像缺少解剖学标记，无法应用到DBS手术导航系统上，不适用于DBS临床治疗的问题，一方面通过在QSM重建过程中采用性能更优的STAR‑QSM算法，来提高QSM图像质量；另一方面利用T1加权图像的解剖信息与QSM图像进行线性融合，克服T1加权图像在脑深部核团没有解剖结构的缺陷，使最终的融合图像能在提供精确的深核靶区定位效果的同时可以很好的融合到现有的DBS定位软件系统中，满足临床需要。

Description

定量磁化率成像及用于DBS潜在刺激靶区定位的方法

技术领域

本发明属于医学成像及手术导航的设计领域，尤其涉及一种定量磁化率图像的重建及脑深部电刺激时治疗靶点的可视化方法。

背景技术

脑深部电刺激(DBS)是通过电极控制电流电压大小，直接对脑部某些核团进行电刺激的一种治疗方法。DBS在帕金森病等运动障碍性疾病的治疗中有着广泛的应用。精确的DBS电极植入是影响DBS治疗结果的关键。DBS手术需要对脑部潜在的深核刺激靶区(GPi、GPe、STN等)精准定位。靶区定位偏差将影响电极在脑内的放置路径和刺激位点，可能造成一些重要神经环路的结构完整性被破坏等后果，影响患者的术后效果。因此，实现脑深层核团的精准定位对于脑深部刺激手术有着重要意义。

实现靶区的精准定位要求靶区具有高的空间分辨率和组织对比度。靶区定位的方法可以分为两大类。一类是利用已经实现脑区精细划分的脑模板，将脑模板配准到个体的磁共振影像上。这种方法操作上比较简单，但是核团的位置、边界、大小是否精确取决于配准结果，很难实现核团的精准定位。另一类利用个体的磁共振影像，直接观察个体的靶区情况。这要求靶区和周围组织有很高的对比度。针对第二类方法已经提出了一些增强靶区和周围组织对比度的方法。比如利用超强磁场(7T)扫描可以得到空间分辨率很高的影像，在T2加权影像上可以清晰分辨靶区的位置和边界。这样的代价是7T扫描花费的时间与低的磁场(1.5T、3T)相比要长，扫描时间的延长导致患者更容易因头部移动而出现伪影；另外7T超高磁场在临床上还未普及。较低的主磁场扫描获得的磁共振影像，如T1加权、T2加权、R2*，靶区和周围组织的对比度仍不理想。

置于磁场中的物质将产生与主磁场同向或反向的磁场，引起物质周围磁场分布发生变化。描述这一性质的物质属性为磁化率。磁化率是物质的固有属性。定量磁化率成像(QSM)可以从梯度回波(GRE)相位图像中定量地计算提取出磁化率数值，是一种新的磁共振成像模态。QSM已经被应用于测量脑部铁沉积、探究膝关节软骨结构、白质纤维的各向异性等生物医学研究，正在受到越来越多的关注。

从相位图像到磁化率图像，QSM还需要经过相位卷缠绕、去除背景场、偶极子反演三个步骤。

磁共振机器采集到的磁共振相位图像，其相位大小被限制在[0,2π)。真实的相位值被缠绕在这个范围内，真实值和测量值相比相差2π的整数倍。因此需要进行相位解缠绕，恢复出真实的相位信息。

根据回波时间(TE)和旋磁比可以计算出感兴趣区域内磁场的变化。感兴趣区域内的磁场分布是区域内局部场和区域外背景场的叠加。目前通用的做法是基于背景场的先验知识，从总磁场中去除背景场。常用的方法包括PDF、SHARP、VSHARP等。

局部场与磁化率之间满足

其中ψ(r)是在图像域中的局部场与主磁场的比值，d(r)是偶极子核，χ(r)表示磁化率的空间分布。局部场是偶极子核函数和磁化率空间分布的卷积。需要指出的是，这里忽略了噪声的影像。通过傅里叶变换将上式转换到k空间，

Ψ(k)＝D(k)·χ(k) (2)

其中kz是k空间在z方向上的分量。从局部场中推导出磁化率数值是一个病态的逆问题：在k空间的偶极子核函数在与主磁场方向成54.7°的魔角时，会出现零值的锥面，因此直接求导数得到的磁化率是不可靠的，需要引入额外的信息进行求解。磁化率的计算在数学上可以看成是求解病态的逆问题。添加正则项可以唯一确定磁化率数值。许多研究人员在正则项的构造和形式上进行了很多工作，目前应用于磁化率反演的正则化方法包括L1范数、L2范数，L1+L2范数等。

脑深层核团相比周围组织铁沉积含量更高，这一特性在QSM图像上有着直观的体现。相比于T2加权、T1加权等磁共振模态，在QSM图像上脑深层核团和周围组织有更高的对比度和更好的可视化结果。但是QSM在经过一系列的处理之后，缺少了解剖学结构(如脑壳)，现有广泛使用的DBS手术导航系统无法识别QSM图像。

一方面，QSM重建技术在磁化率反演时仍存在挑战，得到的磁化率结果准确度不足，有些位置存在尖锐的伪影，不利于临床应用。另一方面，QSM图像缺少解剖学标记，无法应用到DBS手术导航系统上，不适用于DBS临床治疗。

发明内容

本发明的目的是提供一种定量磁化率成像及用于DBS潜在刺激靶区定位的方法，可消除QSM图像中的伪影，提高图像质量，且使QSM图像适用于DBS手术系统，实现刺激靶区的精准定位。

为解决上述问题，本发明的技术方案为：

一种定量磁化率成像的方法，包括：

S1：采用3D GRE序列采集脑部图像，所述脑部图像包括相位图像；

S2：基于拉普拉斯算法，对所述相位图像进行解缠绕处理；

S3：采用VSHARP算法去除解缠绕后的相位图像中的外部磁源干扰，得到目标区域内的局部场；

S4：采用STAR-QSM算法，反演得到所述局部场的第一磁化率；

S5：根据局部场与磁化率的数学关系，得到所述第一磁化率对应的第一局部场；

S6：所述目标区域内的局部场减去所述第一局部场，得到第二局部场；采用STAR-QSM算法，反演得到所述第二局部场的第二磁化率；所述第二磁化率小于所述第一磁化率；

S7：将所述第一磁化率与所述第二磁化率相加，得到重建的定量磁化率图像。

根据本发明一实施例，所述步骤S4进一步包括：

所述STAR-QSM算法构造的求解所述第一磁化率的函数包括保真项及约束项，所述保真项采用L2范数，所述约束项采用L1范数；其中，求解所述第一磁化率的函数如下所示：

其中，FT^-1表示傅里叶逆变换，k_z是k空间在z方向上的分量，χ(k)表示k空间的磁化率分布，χ||₂表示保真项，||W·G·χ||₁表示约束项，λ表示权重系数，W_Gx、W_Gy、W_Gz表示x、y、z三个方向上的梯度算子，G_x、G_y、G_z表示x、y、z三个方向上加权因子。

根据本发明一实施例，所述第一磁化率的权重系数为0.1，所述第二磁化率的权重系数为1×10^-5。

根据本发明一实施例，所述步骤S5进一步包括：

局部场与磁化率之间满足如下数学关系：

其中ψ(r)是在图像域中的局部场与主磁场的比值，d(r)是偶极子核，χ(r)表示磁化率的空间分布；

通过傅里叶变换，得到如下公式：

Ψ(k)＝D(k)·χ(k)

其中，k_z是k空间在z方向上的分量。

一种用于DBS潜在刺激靶区定位的方法，包括：

A1：采用本发明一实施例中的定量磁化率成像的方法，得到QSM图像；

A2：采用3D T1序列采集脑部的T1加权图像；

A3：将所述QSM图像的幅度图像去脑壳处理，得到只包含脑组织的幅度图像；

A4：将所述T1加权图像去脑壳处理，得到只包含脑组织的T1加权图像；

A5：将去除脑壳后的幅度图像与T1加权图像进行线性配准；

A6：将线性配准的T1加权图像的像素归一化；

A7：将归一化后的T1加权图像与所述QSM图像融合并添加脑壳信息，使融合后的图像适用于DBS中深层核团的定位。

根据本发明一实施例，所述步骤A5进一步包括：

将去除脑壳的幅度图像与去除脑壳的T1加权图像处于同一空间；

将所述T1加权图像作为参考图像，将所述幅度图像作为移动图像，以预设的自由度进行对齐，生成一个4×4的二维仿射变换矩阵；

将所述二维仿射变换矩阵作用于所述QSM图像，使所述QSM图像与所述T1加权图像处于同一坐标系。

本发明由于采用以上技术方案，使其与现有技术相比具有以下的优点和积极效果：

1)本发明一实施例中的定量磁化率成像方法，针对现有的QSM重建技术在磁化率反演时，得到的磁化率结果准确度不足，重建后的QSM图像存在尖锐的伪影，不利于临床应用的问题，通过在偶极子反演时将加强磁化率(第一磁化率)及较弱磁化率(第二磁化率)分别求解，最终从局部场中推导出的磁化率为第一磁化率与第二磁化率的和，解决了因QSM算法的数学模型，当两个磁化率的数值相差较大时，在偶极子反演过程中造成伪影的问题。采用本发明的定量磁化率成像方法，可明显消除QSM图像中的伪影，提高QSM的图像质量。

2)本发明一实施例中的用于DBS潜在刺激靶区定位的方法，针对现有的QSM图像缺少解剖学标记，无法应用到DBS手术导航系统上的问题，通过将消除伪影的QSM图像的幅度图像和T1加权图像分别去脑壳，然后将QSM图像与T1加权图像进行线性配准，得到仿射变化矩阵，将仿射变化矩阵作用于QSM，将T1加权图像进行归一化，归一化的T1加权图像与QSM结合并添加脑壳信息，使得到的图像中深核靶区和周围有很高的组织对比度，实现刺激靶区的精准定位。

附图说明

图1为本发明一实施例中的T2加权图像和QSM在定位STN时的图像对比示意图；

图2为本发明一实施例中的T2加权图像和QSM在定位GPi时的图像对比示意图；

图3为本发明一实施例中的定量磁化率成像方法示意图；

图4为本发明一实施例中的偶极子反演算法STAR-QSM的流程图；

图5为本发明一实施例中的定量磁化率成像方法的效果图；

图6为本发明一实施例中的用于DBS潜在刺激靶区定位的方法示意图；

图7为本发明一实施例中的QSM与T1加权图像线性加权处理的效果图；

图8为本发明一实施例中的QSM与T1加权图像线性加权后，与其他模态图在DBS靶点定位效果上的对比图。

具体实施方式

以下结合附图和具体实施例对本发明提出的一种定量磁化率成像及用于DBS潜在刺激靶区定位的方法作进一步详细说明。根据下面说明和权利要求书，本发明的优点和特征将更清楚。

实施例一

在DBS治疗过程中，能否对一些潜在刺激靶区的深层核团(如GPi、GPe、STN等)进行精准定位，将影响最终手术结果。传统的磁共振成像技术很难提供令人满意的对比度，而QSM图像则具有很好的深层核团可视化效果。图1是T2加权图像和QSM图像显示深层核团STN能力的对比，可以看出，QSM图像(图1右)可以清晰显示STN，而T2加权图像(图1左)则不能。从图2也可以看出，相比于T2加权图像(图2右)，QSM图像(图2左)的GPi与周围组织的对比度更高。目前QSM图像有着很好的应用前景，但是重建QSM图像的质量仍存在着问题，如包含尖锐的伪影和噪声。这些问题限制了QSM的进一步应用。

针对上述问题，本实施例提出了一种性能优化的定量磁化率成像方法，可在QSM图像的重建过程中，消除明显的伪影，提高QSM图像的质量。如图3所示，该定量磁化率成像方法包括：

S1：采用3D GRE序列采集脑部图像，该脑部图像包括相位图像；

S2：基于拉普拉斯算法，对相位图像进行解缠绕处理；

S3：采用VSHARP算法去除解缠绕后的相位图像中的外部磁源干扰，得到目标区域内的局部场；

S4：采用STAR-QSM算法，反演得到该局部场的第一磁化率；

S5：根据局部场与磁化率的数学关系，得到第一磁化率对应的第一局部场；

S6：用局部场减去第一局部场，得到第二局部场；采用STAR-QSM算法，反演得到第二局部场的第二磁化率；该第二磁化率小于第一磁化率；

S7：将第一磁化率与第二磁化率相加，得到重建的定量磁化率图像。

具体的，QSM图像的重建基本可以归纳为相位解缠绕、去除背景场及偶极子反演三个步骤。在对QSM图像进行反演之前，需采集源图像。本实施例采用3D GRE序列采集脑部图像是现有的一项技术，在此不作详细的描述。通过3D GRE序列采集脑部图像包括幅度图像和相位图像，其中，对相位图像进行相位解缠绕处理。

由于相位解缠绕的算法众多且各有千秋，或基于解缠相位图与真实相位图越接近越好的标准；或基于解缠相位减去真实相位得解缠误差，越小越好的标准；或基于解缠相位重新缠绕到区间，与原缠绕相位比较，越接近越好；亦或是基于算法运行时间，越少越好的标准。通常，磁共振机器采集到的磁共振相位图像，其相位大小被限制在[0,2π)。真实的相位值被缠绕在这个范围内，真实值和测量值相比相差2π的整数倍n。如下所示：

φ(i,j)＝ψ(i,j)+2πn(i,j)

其中，φ(i,j)表示真实相位，ψ(i,j)表示缠绕相位，(i,j)表示实域像素坐标。为了从缠绕相位通过解缠绕操作得到真实相位，通常是求得每个像素点对应的整数n值，再根据上述公式得到每个像素点继而整幅雷达图像的真实相位。通常，解缠算法是通过计算相邻像素点的相位差值，当该差值超过一定阈值时称为相位跳跃。而求差分的运算就相当于求偏导的一个过程，因此我们可以通过求上述公式的的二维拉普拉斯来近似遵循上述过程，因此上述公式就可表示为：

其中，表示向后二维拉普拉斯算子，/>表示向前二维拉普拉斯算子。

而缠绕相位与解缠相位的关系如下：

可以通过缠绕相位，计算出解缠相位。得到的解缠相位为真实的相位信息，该相位信息为感兴趣区域(目标区域)内、外磁源共同作用的叠加，因此，需要进行去背景场处理，消除感兴趣区域外磁源引发的磁场变化。

去除背景场需要先验知识。VSHARP假定背景场在感兴趣区域内满足谐和函数，该谐和函数遵循球面均值定理(spherical mean value)。球面均值定理认为，满足拉普拉斯方程的函数，对于任意半径的球形核，球面对应的函数平均值等于在球体中心的函数值。利用这一性质可消除来自外部磁源的相位贡献。这一方法在SHARP和VSHARP等方法中已经得到应用。VSHARP与SHARP相比，VSHARP的球形核半径在不同的区域半径是变化的。

本实施例对解缠绕后的相位图像使用VSHARP算法，去除外部磁源的干扰，得到感兴趣区域的局部场。

由于从局部场中推导出磁化率数值是一个病态的逆问题：在k空间的偶极子核函数在与主磁场方向成54.7°的魔角时，会出现零值的锥面，因此直接求导数得到的磁化率是不可靠的，需要引入额外的信息进行求解。磁化率的计算在数学上可以看成是求解病态的逆问题。添加正则化项可以唯一确定磁化率数值。本实施例在求解磁化率时采用了正则化的方法贴加约束项。

对相位图像去除背景场后，需进行偶极子反演，重建QSM图像。由于传统的偶极子反演方法得到的磁化率结果准确度不足，重建后的QSM图像存在尖锐的伪影，不利于临床应用的问题，本实施例提出了一种新的偶极子反演算法，即STAR-QSM(STreaking ArtifactsReduction for Quantitative Susceptibility Mapping)算法。如图4所示，该算法包括：

计算感兴趣区域的全部局部场对应的较强的磁化率χ1(第一磁化率)；

根据公式得到磁化率χ1对应的部分局部场(第一局部场)；

将全部局部场减去部分局部场，得到剩余的局部场(第二局部场)；

计算剩余的局部场对应的较弱的磁化率χ2(第二磁化率)；

最终得到的磁化率为χ＝χ1+χ2。

具体的，在计算感兴趣区域的全部局部场对应的较强的磁化率χ1时，STAR-QSM算法构造的求解函数包括保真项及约束项，其中保真项采用L2范数，保证数据的连续性；约束项采用L1范数，用于计算加权梯度的总变化量；该求解函数如下所示：

χ＝min{||FT^-1(D(k)·χ(k))-x||₂+λ||W·G·χ||₁}

其中，正则化项的表达式为：

其中，FT^-1表示傅里叶逆变换，k_z是k空间在z方向上的分量，χ(k)表示k空间的磁化率分布，χ||₂表示保真项，||W·G·χ||₁表示约束项，λ表示权重系数，W_Gx、W_Gy、W_Gz表示x、y、z三个方向上的梯度算子，G_x、G_y、G_z表示x、y、z三个方向上加权因子。

根据上述公式得到的较强的磁化率χ1，其权重系数取为λ＝0.1。

由于磁化率与局部场之间存在严格的数学关系，即公式 其中ψ(r)是在图像域中的局部场与主磁场的比值，d(r)是偶极子核，x(r)表示磁化率的空间分布；通过傅里叶变换，得到如下公式：Ψ(k)＝D(k)·χ(k)，其中，/>k_z是k空间在z方向上的分量。

根据上述原理，可根据磁化率χ1反推出其对应的第一局部场的大小。

将之前计算磁化率χ1的局部场(全部局部场)减去第一局部场，得到剩余的局部场(第二局部场)，采用与计算磁化率χ1相同的方法，计算第二局部场对应的较弱的磁化率χ2(第二磁化率)。与磁化率χ1不同的是，计算磁化率χ2时，其权重系数取为λ＝1×10^-5。

如图5所示为本实施例中的定量磁化率成像方法的效果图。在进行STAR-QSM算法反演时，将通过3D GRE序列采集脑部图像中的幅度图像去除脑壳，做成掩膜，得到最终提取的磁化率χ。

通过分别求解较强磁化率χ1及较弱磁化率χ2，解决QSM重建过程中的相位变化范围较大的问题。本实施例的STAR-QSM算法可处理变化范围大的相位图像，如相位为0.3ppm的图像。由于最终输出的磁化率为较强磁化率χ1与较弱磁化率χ2的和，与其他算法相比，重建的QSM图像消除了明显的伪影，提高了图像质量。

实施例二

上述实施例一中描述了一种QSM图像的重建方法，可定量提取局部场中的磁化率，消除QSM图像中的伪影，提高QSM的图像质量。上述对QSM图像的处理，均是为了使QSM图像适用于DBS系统中，辅助DBS手术。

为了使DBS系统得到的图像中深核靶区和周围有很高的组织对比度，实现刺激靶区的精准定位，本实施例提出了一种用于DBS潜在刺激靶区定位的方法，将上述实施例一中的重建的QSM图像与T1加权图像结合，既利用QSM图像优秀的深核靶区可视化能力，又结合T1加权图像提供的皮层结构信息，使融合后的图像适用于DBS手术系统，真正实现QSM的临床实用价值。

具体的，如图6所示，该方法包括：

A1：采用如实施例一所述的定量磁化率成像的方法，得到QSM图像；

A2：采用3D T1序列采集脑部的T1加权图像；

A3：将所述QSM图像的幅度图像去脑壳处理，得到只包含脑组织的幅度图像；

A4：将T1加权图像去脑壳处理，得到只包含脑组织的T1加权图像；

A5：将去除脑壳后的幅度图像与T1加权图像进行线性配准；

A6：将线性配准的T1加权图像的像素归一化；

A7：将归一化后的T1加权图像与QSM图像融合并结合脑壳信息，使融合后的图像适用于DBS中深层核团的定位。

具体的，通过3D GRE序列采集QSM重建所需的原始幅度图像和相位图像，采用3DT1w序列采集T1加权图像。其中，3D GRE采用轴状面扫描，FOV＝256×256mm，采集矩阵为320×320，分辨率为0.75×0.75×1.5mm，重复时间(TR)为33ms，从11ms到29.7ms均匀的6个回波时间(TE)。3D T1w采用轴状面扫描，FOV＝256×256mm，采集矩阵为256×256，分辨率为1×1×1mm，TR＝6.9ms，TE＝3.4ms。

3D GRE的相位图像根据实施例一所述的步骤进行相位解缠绕、去除背景场、偶极子反演得到QSM图像。

然后，将3D GRE的幅度图像和T1加权图像利用FSL的bet函数去除脑壳。

由于QSM不包含脑壳，且QSM直接与T1加权图像配准的效果不好，本实施例采用间接的方式对QSM图像及T1加权图像进行线性配准，如图7所示。

去除脑壳的幅度图像首先与去除脑壳的T1加权图像进行线性配准，幅度图像将与T1加权图像处于同一空间，使幅度图像与T1加权图像在空间上对齐。在线性配准的过程中，T1加权图像为参考图像，幅度图像为移动图像，配准的自由度设为12。在线性配准的过程中，会生成一个4×4的二维仿射变换矩阵，将该二维仿射变换矩阵作用于根据上述实施例一的处理方法得到的QSM图像，使QSM图像与T1加权图像处于同一坐标系。其中，仿射变换是一种二维坐标到二维坐标之间的线性变换，可保持二维图像的平直性(即变换后直线还是直线不会弯曲)和平行性(即二维图像间的相对位置关系不变，平行线还是平行线，相交直线的夹角不变)。本实施例中的二维仿射变换矩阵可通过平移，尺度放大，旋转等操作作用于QSM图像。

虽然QSM图像与T1加权图像处于同一坐标系，但是，QSM图像与T1加权图像的数值范围不一致，为了便于两者的融合，本实施例将T1加权图像的像素进行归一化处理。

将归一化处理后的T1加权图像与QSM图像线性加权融合并填补脑壳信息，使融合后的图像适用于DBS手术系统，真正实现QSM的临床实用价值。

由于本实施例采用3D T1w序列采集T1加权图像，其得到的图像记为T1w图像，将T1w图像与QSM图像线性加权融合后，得到的混合图像记为Hybrid，该Hybrid的计算公式如下：

Hylbrid＝μT1w-QSM

其中，μ为权重系数，用于调整T1w图像的像素值范围，使融合后的图像过渡自然。本实施例中取μ＝0.4。

如图8所示，相比于CT图像(图8.A)、T2加权图像(图8.B)和T1加权图像(图8.C)，混合Hybrid图像(图8.D)可以清晰显示GPi和GPe的边界，取得了较好的定位效果。在确认DBS电极路径和刺激位点时，能提供较准确的空间位置信息(如图8.E-8.G所示)。

综上，本发明一方面在QSM重建过程中采用性能更优的STAR-QSM算法，以提高QSM图像质量；另一方面利用T1加权图像的解剖信息与QSM图像进行线性融合，克服QSM图像没有解剖结构的缺陷，使最终的融合图像能在提供精确的深核靶区定位效果的同时，满足临床需要。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式。即使对本发明作出各种变化，倘若这些变化属于本发明权利要求及其等同技术的范围之内，则仍落入在本发明的保护范围之中。

Claims (2)

1.一种用于DBS潜在刺激靶区定位的方法，其特征在于，包括：

A1：通过3D GRE序列采集重建QSM所需的原始幅度图像和相位图像，对相位图像进行相位解缠绕、去除背景场及基于STAR-QSM算法的偶极子反演处理后得到重建的QSM图像；

A2：采用3D T1序列采集脑部的T1加权图像；

A3：将所述重建的QSM图像去脑壳处理，得到只包含脑组织的幅度图像；

A4：将所述T1加权图像去脑壳处理，得到只包含脑组织的T1加权图像；

A5：将去除脑壳后的幅度图像与T1加权图像进行线性配准；

A6：将线性配准的T1加权图像的像素归一化；

A7：将归一化后的T1加权图像与所述重建的QSM图像融合，并添加脑壳信息，使融合后的图像适用于DBS中深层核团的定位软件的对影像的需求。

2.如权利要求1所述的用于DBS潜在刺激靶区定位的方法，其特征在于，所述步骤A5进一步包括：

将去除脑壳的幅度图像与去除脑壳的T1加权图像处于同一空间；

将所述T1加权图像作为参考图像，将所述幅度图像作为移动图像，以预设的自由度进行对齐，生成一个4×4的二维仿射变换矩阵；

将所述二维仿射变换矩阵作用于所述QSM图像，使所述QSM图像与所述T1加权图像处于同一坐标系。