一种利用时空信息的磁共振动脉自旋标记序列部分容积校正方法与流程

本发明属于医学影像的部分容积效应校正方法技术领域，特别涉及一种利用时空信息的磁共振动脉自旋标记序列的部分容积校正方法。

背景技术：
磁共振(MagnaticResonanceImaging,MRI)灌注成像是用来反映组织的微血管分布及血流灌注情况的MRI检查技术。动脉自旋标记技术(ArterialSpinLabeling,ASL)，通过磁化血液中水分子作为内源性对比剂，来评价脑血流灌注情况，可完全无创检测，可重复性高，且可获得绝对定量的脑血流(CerebralBloodFlow,CBF)，完全不受血脑屏障的影响。ASL技术以磁标记的动脉血作为内源性对比剂，在成像平面上游，利用反脉冲标记动脉血的质子；延迟一段时间后，待标记的血液进入组织，血液与组织进行物质交换后成像(即标记图像，label)。label图像包括静态组织和标记动脉血的信号。为了消除静态组织的信号，进行另外一次未标记血成像(即控制图像，control)，control图像只包括静态组织信号。将label图像与control图像减影，所得的差值图像只与流入成像平面的标记动脉血有关。ASL图像的缺点在于它对噪声敏感。一般说来，label图像与control图像的差值大约是control像灰度值的1～2％。因此，一对label/control图像的差值不能够有效的反映灌注情况，需要大约60对数据来获得合理的信噪比(SignaltoNoiseRatio,SNR)。同时，作为一组配对图像，label和control图像需要快速、连续获取。所以，快速成像方法(如EPI方法)普遍应用于ASL成像，但这类方法的使用却降低了图像的空间分辨率，使图像受部分容积(PartialVolume,PV)效应影响严重。目前，针对ASL图像的PV校正算法很少，Asllani等(AsllaniI,etal.MagnResonMed,2008；60(6):1362-1371.)首先使用局部区域的线性回归(LinearRegression,LR)方法对图像进行校正，但该方法却使局部区域平滑，易引起CBF计算的误差。之后，Liang等(LiangX,etal.MagnResonMed,2013；69(2):531-537.)采用最小截平方和方法对平滑效应进行后处理。同时，Chappell等(ChappellMA,etal.MagnResonMed,2011；65(4):1173-1183.)对多个反转时间的ASL序列的PV进行校正。虽然上述方法已可以在空间域对ASL数据进行一定程度的PV校正，但ASL序列中的时间信息并未被充分使用。

技术实现要素：
针对目前ASL序列的PV校正方法中，时间信息利用率低的问题，本发明提供了一种利用时空信息的ASL序列PV校正方法，该方法充分利用ASL序列中的时间和空间信息，对序列进行有效PV校正，从而大幅提高ASL图像质量。本发明是通过以下技术方案来实现：一种利用时空信息的磁共振动脉自旋标记序列的部分容积校正方法，包括以下步骤：(1)采集被试者的MRI数据，包括结构像和动脉自旋标记序列；(2)将结构像和动脉自旋标记序列进行配准；(3)利用SPM软件对结构像进行分割，获取灰质、白质和脑脊液的概率分布图像；(4)利用线性回归方法对动脉自旋标记序列进行空间部分容积校正，并利用校正结果作为最大期望算法的初始值，在时间方向上进行部分容积校正，充分利用动脉自旋标记序列的空间和时间信息，获得准确的校正结果。步骤(1)中，采集的结构像为T1或T2序列，以及动脉自旋标记序列。步骤(2)所述将结构像和动脉自旋标记序列进行配准，是利用SPM提供的MNI坐标系作为中间值来进行配准。步骤(4)利用动脉自旋标记的四维数据，沿着时间轴，将具有相同三维空间坐标i的体素构成一个时间向量；假定该向量中所有元素相互独立，且其中的灰质和白质符合高斯分布，利用最大期望算法估计混合组织中灰质和白质灌注信号，具体求解方法如下：对于某一空间体素i，磁矩表示为：ΔMi＝PiGMΔMiGM+PiWMΔMiWM(1)其中，PiGM和PiWM是混合体素中灰质和白质的概率；ΔMiGM和ΔMiWM代表灰质和白质的灌注信号，灌注信号是通过label/control图像的差值ΔM来描述，且某一空间体素的脑血流与ΔM和M0图像的比值ΔM/M0等价；则根据动脉自旋标记成像的两腔室模型得到：ftissue＝(ΔM/M0)·Ftissue；其中，Ftissue是与血脑相关的血流灌注系数；M0图像是与灌注图像序列分别扫描的，采用T1序列的形式，并通过与灌注图像配准，转换成与灌注图像相同大小；故假设M0图像未受到部分容积效应的影响，得到如下关系式：其中，PGM和PWM通过动脉自旋标记序列与同一被试的结构图像进行配准后获得，FGM和FWM与图像的成像参数有关；ΔMiGM和ΔMiWM为灰质和白质的灌注信号；沿时间轴，将具有相同三维空间坐标i的体素构成一个时间向量{Yit,t＝1,…,T}，T代表了时间向量的维数；在观察值Yit中，灰质和白质的分量表示为：Yit＝XitGM+XitWM(3)其中，XitGM和XitWM分别是以均数为和方差为和的随机变量，假设所有的T个体素是相互独立的，则有：将灌注模型和式(4)结合，得到：其中，SiGM和SiWM分别代表了灰质和白质灌注信号的方差；假设XitGM和XitWM服从高斯分布，则公式(4)变换为:在最大期望算法中，第t个体素的观察值Yit是一个不完整的随机变量；XitGM和XitWM代表的是第t个体素中完整的混合组织信息，是一个完整变量；在以公式(1)为条件进行积分，建立不完整变量{Yit}与完整变量{XitGM}和{XitWM}之间的概率分布的关系，如下式：其中，i代表某一体素三维空间位置，{Yit}代表具有相同空间位置i的体素所构成一个时间向量；{XitGM}和{XitWM}分别为观察值{Yit}中灰质和白质的分量。式(6)的最优解获取步骤如下：采用最大期望算法求取完全变量条件概率分布的最大期望，E-step是用于对变量的对数似然估计，M-step用于求取期望最大值：E-step：计算条件概率期望值p(X|Θ)，其中条件期望表示为：根据条件期望的推导，得出：M-step：通过n+1次迭代使条件概率期望值最大化，混合组织模型中的均值通过最大化条件概率的n+1次迭代求出：步骤(4)具体操作如下：1)对于每幅三维灌注图像，使用线性回归方法进行空间部分容积校正；2)基于线性回归方法的空间部分容积校正结果，将具有相同空间位置i的体素构成时间向量；3)对于时间向量，利用最大期望方法进行时间轴上的混合组织估计；4)对空间位置i+1，重复第2)和第3)步，直至整幅灌注图像都被校正。与现有技术相比，本发明具有以下有益的技术效果：本发明的基于时空信息的ASL序列部分容积校正方法，首先，采集被试者的MRI数据，将采集的结构像和ASL序列进行配准；然后，利用SPM软件对结构像进行分割，分别生成灰质(graymatter,GM)、白质(whitematter,WM)和脑脊液(cerebrospinalfluid,CSF)的概率分布图像，能够为后续的PV校正提供GM和WM的概率分布信息，有利于估计混合组织的参数；再次，利用LR方法对ASL序列进行空间PV校正；最后，利用LR方法的空间校正结果作为最大期望(expectationmaximization,EM)算法的初始值，并在时间方向上进行校正，从而准确估计脑血流值，有助于后续的数据分析。该方法(EM-LR)充分利用了动脉自旋标记序列的空间和时间信息，弥补了EM算法对初始值敏感的问题，加速了EM迭代算法的收敛速度。本发明提供的利用时空信息的部分容积校正方法，不仅适用于ASL序列，也适用于其他具有时间信息影像数据的部分容积校正。附图说明图1为本发明利用时空信息的磁共振动脉自旋标记序列部分容积校正方法的步骤流程图；图2为使用未校正、EM方法、LR方法和EM-LR方法对模拟数据进行校正得到的GM脑血流分布图的中间层结果图；图3为在不同噪声条件下，GM的CBF感兴趣区域(ROI)分析的结果；其中，(a)噪声强度为2.5时的ROI分析，(b)噪声强度为6.5时的ROI分析；图4为EM、LR和EM-LR方法，分别在在高灌注和低灌注区域中心线的校正结果；其中，(a)和(b)分别为在噪声强度为2.5时，高灌注和低灌注区域的中心线上GM的CBF变化，(c)和(d)分别为在噪声强度为6.5时，高灌注和低灌注区域的中心线上GM的CBF变化；图5对于模拟仿真2，LR、EM和EM-LR的校正结果。具体实施方式下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。参见图1，本发明公开的一种利用时空信息的磁共振动脉自旋标记序列部分容积校正方法，包括下述步骤：(1)进行被试的准备工作后，采集被试的MRI数据，包括结构像和ASL序列；(2)将结构像和ASL序列进行配准；(3)利用SPM(StatisticalParametricMapping)软件对结构像进行分割，获取GM、WM和CSF的概率分布图像；(4)利用LR方法对ASL序列进行空间PV校正，并将其作为EM算法的初始值，在时间方向上进行PV校正(EM-LR算法)，从而充分利用ASL序列的空间和时间信息，获得准确的校正结果。上述方法中，步骤(1)所述采集的结构像可以是T1或T2序列。根据动脉血反转标记的方法不同，ASL技术可大体分为两种类型:连续式动脉自旋标记(ContinuousArterySpinLabeling,CASL)和脉冲式动脉自旋标记(PulsedArterySpinLabeling,PASL)。本发明可针对任一类型的ASL序列进行PV校正，一般情况下，ASL序列的重复扫描数据大约为60幅，以便获取合理的信噪比。步骤(2)所述的结构像和ASL序列进行配准，利用SPM提供的MNI(MontrealNeurologicalInstitute)坐标系作为中间值来进行配准。使用SPM读入结构像和ASL数据时，会产生图像坐标系与MNI坐标系之间的转换矩阵，MNI坐标对应的是标准模板，可以用于结构像和ASL数据的配准。步骤(3)所述的结构像分割，利用SPM软件提供的标准流程进行分割，分别生成GM、WM和CSF的概率分布图像。分割获取的概率分布图像，可为后续的PV校正提供GM和WM的概率分布信息，有利于估计混合组织的参数。步骤(4)所述的基于EM算法的时间方向PV校正。利用ASL的四维数据，沿着时间轴，将具有相同三维空间坐标i的体素可以构成一个时间向量。假定该向量中所有元素相互独立，且其中的GM和WM符合高斯分布，利用EM算法估计混合组织中GM和WM灌注信号，具体过程如下：由于ASL图像的空间分别率低，在一个混合体素中可能包含了GM、WM和CSF的共同作用。一般认为，CSF对ASL序列中的CBF没有贡献，因此，对于某一空间体素i，磁矩可以表示为ΔMi＝PiGMΔMiGM+PiWMΔMiWM(1)其中，PiGM和PiWM是混合体素中GM和WM的概率；ΔMiGM和ΔMiWM代表GM和WM的信号。由ASL图像成像原理可知，灌注信号是通过label/control图像的差值ΔM来描述的，而且某一体素的CBF与ΔM和M0图像的比值(ΔM/M0)等价。因此，根据ASL成像的两腔室模型可以得到，ftissue＝(ΔM/M0)·Ftissue，其中Ftissue是与血脑相关的血流灌注系数。M0图像是与灌注图像序列分别扫描的，采用的是T1序列的形式，并通过与灌注图像配准，转换成与灌注图像相同大小，因此M0图像受PV效应的影响很小，故本研究假设M0图像未受到PV效应的影响。由此，可以得到如下关系式：在上述脑血流测量公式中，PGM和PWM可通过ASL序列与同一被试的结构图像进行配准后获得(步骤3)，FGM和FWM与图像的成像参数有关，接下来就是如何计算GM和WM的灌注信号ΔMiGM和ΔMiWM。为了获取ΔMiGM和ΔMiWM，利用ASL的四维数据，沿着时间轴，将具有相同三维空间坐标i的体素构成一个时间向量{Yit,t＝1,…,T}，T代表了时间向量的维数。在观察值Yit中，GM和WM的分量可表示为XitGM和XitWM，因此，Yit＝XitGM+XitWM(3)其中，XitGM和XitWM分别是以均数为和方差为和的随机变量。假设所有的T个体素是相互独立的，则将灌注模型和上述公式结合，可以得到其中SiGM和SiWM分别代表了GM和WM信号的方差。假设XitGM和XitWM服从高斯分布，则公式(4)可变换为：在EM算法中，第t个体素的观察值Yit是一个不完整的随机变量。而XitGM和XitWM代表的是第t个体素中完整的混合组织信息，因此是一个完整变量。在以公式(3)为条件进行积分，就可以表示不完整变量{Yit}与完整变量{XitGM}和{XitWM}之间的概率分布的关系，为了获取方程(6)的最优解，EM算法是用来求取完全变量条件概率分布的最大期望，其中E-step是用于对变量的对数似然估计，M-step用于求取期望最大值。E-step：计算条件概率期望值p(X|Θ)，其中条件期望可以表示为：根据条件期望的推导，可以得出：M-step：通过n+1次迭代使条件概率期望值最大化。混合组织模型中的均值可以通过最大化条件概率的n+1次迭代求出，即通过上述方法，可以利用ASL序列的时间信息，获取位置i的GM和WM的灌注信号ΔMiGM和ΔMiWM，但由于EM算法对初始值比较敏感，使用相对准确的初始值，可加速算法的收敛，提高EM算法的计算准确性。为此，本发明提出两种初始值设置方案：1)EM方法：利用文献中查到的空间先验信息作为算法的初始值，即E-step时，采用对EM方法进行初始化，并利用EM方法(公式8-15)进行时间轴上的混合组织估计。2)EM-LR方法：将基于LR方法的空间校正和基于EM算法的时间校正相结合，利用LR方法的校正结果为EM算法提供初始值。具体过程如下：①对于每幅三维灌注图像，使用LR方法(AsllaniI,etal.MagnResonMed,2008；60(6):1362-1371.)进行空间PV校正；②基于LR方法的空间校正结果，具有相同空间位置i的体素构成时间向量；③对于时间向量，利用EM方法(公式8-15)进行时间轴上的混合组织估计；④对空间位置i+1，重复第2)和第3)步，直至整幅灌注图像都被校正。如上所述，就完成了利用时空信息的ASL序列PV校正，该方法充分利用了ASL序列的时间和空间信息，弥补了EM算法对初始值敏感的问题，加速了EM迭代算法的收敛速度。模拟仿真1：利用模拟仿真生成高低灌注区域，使用本专利公开的方法和传统方法对模型进行PV校正，从而评价本方法的有效性。仿真模型生成步骤：在使用SPM软件对磁共振T1序列进行标准化和分割的基础上，获取GM、WM和CSF的概率分布图像。ASL序列采用如下方式进行模拟：1)对GM和WM图像进行标准化，图像大小为60×72×60，体素大小为3×3×3mm3；2)WM的CBF值设置为ΔM＝20ml/100g/min；3)GM的CBF值设置为ΔM＝60ml/100g/min，并选取2个半径为5个体素的球状区域，分别设置为低灌注和高灌注区，并定为30和90ml/100g/min；4)通过公式(1)，获取整幅三维灌注图像；5)采用两种高斯噪声强度2.5(文献中常用的值)和6.5(文献中报道的ASL序列最强的噪声)叠加到三维灌注成像中；6)根据常规ASL序列的重复扫描次数(60次)，重复第1-5步，生成两组四维ASL序列。仿真结果：参照图2，使用未校正、EM方法、LR方法和EM-LR方法对模拟数据进行校正。从图中可以看出，EM、LR和EM-LR方法优于未校正的结果，尤其在GM和CSF的交界处，但LR方法得到的结果在高/低灌注区域的边缘具有明显的平滑作用。参照图3，感兴趣区域(RegionofInterest,ROI)分析，其中，(a)噪声强度为2.5时的ROI分析，(b)噪声强度为6.5时的ROI分析。为了定量比较不同方法的CBF测量准确性，采用GM平均CBF在GM概率分布中的一致性来进行评价。基于GM的概率分布图像，选取9个区域，即PGM在[10％～20％],[20％～30％],……[90％～100％]范围，并分别计算每个范围内GM的平均CBF值。图3展示了不同噪声水平下，模拟数据的ROI分析结果。与图2的观察结果类似，EM、LR和EM-LR方法得到的结果与真实值接近，但在GM概率分布较低的区域，EM-LR方法和LR方法获得的结果较好。参照图4，为了研究EM、LR、EM-LR方法在边缘保持方面的作用，分别在高灌注和低灌注区域选取两个穿过中心的剖面。其中，(a)和(b)分别为在噪声强度为2.5时，高灌注和低灌注区域的中心线上GM的CBF变化，(c)和(d)分别为在噪声强度为6.5时，高灌注和低灌注区域的中心线上GM的CBF变化。图4展示了在剖面与图2所示的中间层相交剖线上，真实值与EM、LR和EM-LR方法的估计值之间的比较。可以看出，EM和EM-LR方法具有较好的边缘保持特性，可以在有效保持细节的同时准确估计脑血流值。模拟仿真2：为了研究本发明方法对小区域和/或轻微灌注改变区域的校正准确性，将“模拟仿真1”的第3)步修改为：GM的CBF值设置为ΔM＝60ml/100g/min，并在其中，①选取一个半径为5个体素的球状区域，设置ΔM＝65ml/100g/min(轻微高灌注)；②选择一个3×3×3的立方体区域，设置ΔM＝55ml/100g/min(小区域，轻微低灌注)；③选择一个2×2×2的立方体区域，设置ΔM＝65ml/100g/min(小区域，且轻微高灌注)。仿真结果：参照图5，虽然所研究的区域比较小，且灌注改变不大，但使用EM方法和EM-LR方法都可以获得很好的校正结果，而LR方法却很难对小的灌注异常区域进行校正。通过上述两个模型仿真可知，本发明公开的EM(使用经验值作为EM算法的初始值)和EM-LR(使用LR方法的空间校正结果作为EM算法的初始值)可以对ASL数据进行PV校正，且对于区域较小的灌注异常具有特异性。但由于人的个体差异，经验值可能会与真实值差异较大，所以利用LR方法作为EM算法的初始值(EM-LR)更加合理。通过仿真实验证明，本发明提供的利用时空信息的部分容积校正方法，可以对动脉自旋标记序列进行有效的部分容积校正，同时对于范围较小的灌注异常区域具有特异性。同时，该方法也适用于其他具有时间信息的影像数据的部分容积校正。