CN101718848A - 一种优化的磁共振灌注成像图像中血容量量化方法 - Google Patents

Abstract

本发明涉及一种优化的磁共振灌注成像图像中血容量量化方法，技术特征在于：在磁共振灌注成像图像的颅内主中央动脉选取具有典型双峰形状的动脉输入函数。将造影剂浓度时间序列信号与动脉输入函数做去卷积运算，得到具有近似单峰形状的残余函数。利用gamma函数拟合所选取的动脉输入函数，得到动脉输入函数的首次通过。将拟合得到的动脉输入函数的首次通过与残余函数做卷积运算，从而重建出初始时间-显影剂浓度曲线中的首次通过分量。利用重建的首次分量来量化血容量，抑制量化误差。本发明将上述首次通过重建算法与造影剂渗漏效应矫正算法相结合，进一步抑制血容量量化误差。

Description

一种优化的磁共振灌注成像图像中血容量量化方法

技术领域

本发明涉及一种优化的磁共振灌注成像图像中血容量量化方法，是一种减小在利用灌注核磁共振成像量化血容量中因造影剂渗漏效应及血液循环所带来的量化误差的方法。属于磁共振灌注成像中血容量的医学图像处理。

背景技术

现代医学影像学已不满足于对病变形态变化的研究，而向反映组织和器官生理与病理变化甚至功能成像方面发展。脑核磁共振灌注成像是脑功能成像方法之一，它通过研究脑组织的血流灌注状态以及组织血管化程度来获取微循环信息从而揭示肿瘤的生物学行为及预后情况。因此灌注成像在脑部肿瘤的诊断、鉴别诊断、疗效评价及判断预后中发挥主要作用。局部血容量rCBV是灌注成像中一个简便而有效的分析参数，它反映了脑肿瘤微血管指标，与组织学上血管紧密度及肿瘤分级具有显著相关性。然而，rCBV的量化正比于血液首次通过所致的时间-造影剂浓度曲线的积分。二次及后续的多次通过将会带来显著的过估计误差。另一方面，造影剂经血管渗漏而产生的附加T2*顺磁效应使得上述误差进一步放大。

近些年来研究学者已经提出一些旨在去除血液循环引起的rCBV量化误差的方法，但仍存在很多问题。方法1：人为的在时间-造影剂浓度曲线中选取固定的截断点，仅对时间-造影剂浓度曲线进行部分积分。其特点在于实现简单，无需附加计算。缺点在于截断的时间-造影剂浓度曲线无法反映其真实形状。且由于延迟及通过时间差异的存在，使得选取的截断点无法适用于图像中所有的像素点；方法2：基于独立成分分析ICA的方法。该方法利用ICA将时间-造影剂浓度分解为首次通过，二次及多次通过的组合，以及噪声分量。其缺点在于需要人工的干预来确定哪一部分是首次通过分量，无法适用于自动的计算整幅图像像素的血容量；方法3：利用gamma函数拟合时间-造影剂浓度曲线，得到其首次通过分量。该方法对图像质量要求较高，适用于血容量较大、信噪比较高的图像区域。缺点在于在血容量较低的区域，信噪比降低，拟合误差较大。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种优化的磁共振灌注成像图像中血容量量化方法，提出一种用于磁共振灌注成像图像中重建时间-造影剂浓度曲线中首次通过分量的方法，并与造影剂渗漏效应矫正算法相结合，以达到有效抑制rCBV量化中的误差的目的。

本发明的思想在于：血液循环将会在时间-造影剂浓度曲线中形成多峰，传统的动脉输入函数为单峰的gamma函数。时间-造影剂浓度曲线与动脉输入函数做去卷积运算后，得到的残余函数中存在多峰。不同与传统的动脉输入函数选取方法，本发明在磁共振灌注成像图像中颅内主中央动脉选取具有典型双峰形状的动脉输入函数。将造影剂浓度时间序列信号与动脉输入函数做去卷积运算，得到具有近似单峰形状的残余函数。利用gamma函数拟合所选取的动脉输入函数，得到动脉输入函数的首次通过。将拟合得到的动脉输入函数的首次通过与残余函数做卷积运算，从而重建出初始时间-显影剂浓度曲线中的首次通过分量。利用重建的首次分量来量化血容量，抑制量化误差。同时，造影剂泄漏效应在脑肿瘤磁共振灌注成像中一个不可避免的因素，而造影剂渗漏将在颅内组织中带来附加的顺次效应，带来血容量的过估计。本发明将上述首次通过重建算法与造影剂渗漏效应矫正算法相结合，进一步抑制磁共振灌注成像图像中的处理中血容量量化误差。

技术方案

一种优化的磁共振灌注成像图像中血容量量化方法，其特征在于步骤如下：

步骤1：根据 $C(x,y,z,t)=-\frac{1}{\mathrm{TE}}\ln \left(\frac{S(x,y,z,t)}{S_0(x,y,z)}\right),$ 将图像中三维像素的时间-灰度信号曲线转换为时间-造影剂浓度曲线，其中：(x，y，z)为三维像素在图像矩阵中的索引值，C(x，y，z，t)为索引值(x，y，z)处三维像素的时间-造影剂浓度曲线，TE为磁共振灌注成像图像的回波时间，S(x，y，z，t)为索引值(x，y，z)处三维像素的时间-灰度信号曲线，S₀(x，y，z)为造影剂到达之前索引值(x，y，z)处三维像素的基准灰度值，ln(·)为对数运算符；所述的基准灰度值为该像素时间-灰度信号曲线中造影剂达到之前的时间点的均值；

步骤2：根据造影剂渗漏效应矫正算法，矫正三维像素的时间-造影剂浓度曲线中因造影剂渗漏效应带来的附加顺磁效应，得到矫正后的时间-造影剂浓度曲线C_corr(t)；

步骤3：在图像中脑主中央动脉邻域选取一个像素，将该像素的时间-造影剂浓度曲线作为动脉输入函数f(t)；所述选取像素的时间-造影剂浓度曲线具有典型双峰形状；

步骤4：采用标准奇异值分解算法，以步骤3得到的f(t)为输入函数，对步骤2得到的C_corr(t)做去卷积运算，得到残余函数R(t)；

步骤5：采用gamma函数以非线性最小二乘算法拟合步骤3中得到的动脉输入函数f(t)，得到f(t)的首次通过f₁(t)；

步骤6：将步骤5中得到的f₁(t)与步骤4中得到的R(t)做卷积运算以重建步骤2中的C_corr(t)的首次通过分量C₁(t)；

步骤7：根据 $\mathrm{rCBV}\∝{\∫}_{t_0}^t{C}_1\left(\mathrm{\τ}\right)\mathrm{d\τ},$ 利用步骤6中得到的C₁(t)来量化该像素点处的局部血容量rCBV，其中rCBV为局部血容量，t₀为造影剂到达时刻，τ为积分常数。

步骤2中的造影剂泄漏效应校正算法采用文献“Boxerman，J.L.，K.M.Schmainda，and R.M.Weisskoff，Relative cerebral blood volume maps corrected for contrastagent extravasation significantly correlate with glioma tumor grade，whereasuncorrected maps do not.AJNR Am J Neuroradiol，2006.27(4)：p.859-67.”中公布的算法。

有益效果

本发明提出的一种优化的磁共振灌注成像图像中血容量量化方法，首先，随着磁共振成像设备的不断更新，获取的灌注成像数据的信噪比越来越高。信噪比的提高有效的提高去卷积运算的精度，从而提高重建时间-造影剂浓度曲线的首次通过分量的精度。其次，动脉输入函数与普通像素点的时间-造影剂浓度曲线相比，具有更高的信噪比，利用gamma函数拟合动脉输入函数，在准确的反映动脉输入函数的真实形状的同时，有效的提取其首次通过分量。最后，与造影剂渗漏效应矫正算法相结合，更近一步抑制血容量量化误差。

本发明相对于其它方法具有以下优点：1、重建的时间-造影剂浓度的首次通过更接近其真实形状，使得rCBV的量化更加接近真实值；2、与造影剂渗漏效应矫正算法相结合，有效的抑制rCBV的过估计。

附图说明

图1：本发明方法的基本流程图。

图2：模拟的时间-造影剂浓度曲线、其首次通过以及重建的首次通过示例。

图3：血容量的量化值与真实值之间的变化关系。

图4：两张同一位脑胶质细胞瘤患者手术后的两次对比度增强结构扫描的磁共振灌注成像图。

图5：有无误差抑制的血容量图对比

(a)无误差抑制的血容量图；

(b)利用本发明方法所计算的血容量图。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本实施例中：首先将灌注成像中三维像素的时间-灰度曲线转换为时间-造影剂浓度曲线。然后，利用造影剂渗漏矫正算法去除时间-造影剂浓度曲线中的造影剂渗漏效应。其次，在颅内主中央动脉邻域选取具有典型双峰形状的动脉输入函数，将所得的动脉输入函数输入函数，对前一步所得的时间-造影剂浓度曲线做去卷积运算，得到具有近似单峰形状的残余函数。再次，利用gamma函数拟合所选取的动脉输入函数，得到其首次通过分量。最后，将动脉输入函数的首次通过与残余函数做卷积运算，重建出时间-造影剂浓度曲线的首次通过。利用重建的首次通过，量化rCBV。

所述的造影剂渗漏矫正算法如下。首先，在图像中选取正常组织区域并计算区域中所有三维像素的时间-造影剂浓度曲线的均值，将其看作被试个体正常的时间-造影剂浓度曲线。其次，受造影剂渗漏效应影响的时间-造影剂浓度曲线被建模为正常时间-造影剂浓度曲线与其在时间轴上积分的线性加权组合。最后，通过线性最小二乘估计算法计算出最优权重，从而达到矫正时间-造影剂浓度中的造影剂渗漏效应的目的。

所述的去卷积运算采用标准奇异值分解sSVD算法。动脉输入函数的gamma函数拟合采用非线性最小二乘估计算法。

本实施例我们用C++语言实现了一个原型系统。图像数据的来源为接受放射疗法的脑胶质细胞瘤病人的磁共振灌注成像。成像方式为梯度回波(gradient-echo)平面回波(echo planar)扫描。重复时间TR＝2000ms，回波时间TE＝90ms。

本发明整个流程可以参考附图1，具体的实施步骤如下：

1.时间-灰度信号曲线——时间-造影剂浓度曲线转换

时间-灰度信号曲线到时间-造影剂浓度曲线的转换根据以下公式进行：

$C(x,y,z,t)=-\frac{1}{\mathrm{TE}}\ln \left(\frac{S(x,y,z,t)}{S_0(x,y,z)}\right)---\left(1\right)$

式中(x，y，z)为某三维像素在图像矩阵中的索引值；C(x，y，z，t)为(x，y，z)处三维像素的时间-造影剂浓度曲线；TE为回波时间；S(x，y，z，t)为(x，y，z)处三维像素的时间-灰度信号曲线；S₀(x，y，z)为造影剂到达之前(x，y，z)处像素的基准灰度值，计算中取该像素时间-灰度信号曲线中造影剂达到之前的时间点的均值；ln(·)为对数运算符。

2.造影剂渗漏效应矫正

受造影剂渗漏效应影响的时间-造影剂浓度曲线被建模为正常时间-造影剂浓度曲线与其在时间轴上积分的加权线性组合，如式(2)所示：

$C_{\mathrm{eff}}(x,y,z,t)=K_1(x,y,z)\·\stackrel{\‾}{C_{\mathrm{nor}}\left(t\right)}-K_2(x,y,z){\∫}_0^t\stackrel{\‾}{C_{\mathrm{nor}}\left(\mathrm{\τ}\right)}\mathrm{d\τ}---\left(2\right)$

式中(x，y，z)为三维像素在图像矩阵中的索引；C_eff(x，y，z，t)为该像素处受造影剂渗漏效应影响的时间-造影剂浓度曲线，根据式(1)计算得到；K₁(x，y，z)和K₂(x，y，z)为加权系数；τ为积分常数；C_nor(t)反映了该被试个体磁共振灌注成像中不受造影剂渗漏效应影响像素的时间-造影剂浓度曲线的形状，计算中，C_nor(t)为该被试个体颅内不受肿瘤病变影响的正常区域像素的时间-造影剂浓度曲线的均值。式(2)中，C_eff(x，y，z，t)及C_nor(t)为已知变量，K₁(x，y，z)和K₂(x，y，z)为未知变量。利用线性最小二乘估计方法计算出K₁(x，y，z)和K₂(x，y，z)的最优解，进而该像素处受造影剂渗漏效应影响的时间-造影剂浓度曲线可被矫正为：

$C_{\mathrm{corr}}(x,y,z,t)=C_{\mathrm{eff}}(x,y,z,t){+K}_2(x,y,z){\∫}_0^t\stackrel{\‾}{C_{\mathrm{nor}}\left(\mathrm{\τ}\right)}\mathrm{d\τ}---\left(3\right)$

式中C_corr(x，y，z，t)为矫正后的时间-造影剂浓度曲线。C_eff(x，y，z，t)为该像素处受造影剂渗漏效应影响的时间-造影剂浓度曲线；K₂(x，y，z)为前述步骤中加权系数的最优解；τ为积分常数；C_nor(t)为该被试个体颅内不受肿瘤病变影响的正常区域像素的时间-造影剂浓度曲线的均值。

3.求取残余函数的去卷积方法

颅内组织的血液动力学模型为：

$C_{\mathrm{corr}}\left(t\right)=\mathrm{CBF}\·C_a\left(t\right)\⊗R\left(t\right)---\left(4\right)$

式中C_corr(t)为矫正后的时间-造影剂浓度曲线；CBF为血流量参数；R(t)为残余函数；

为卷积运算符；C_a(t)为动脉输入函数。计算中，在颅内主中央动脉邻域选取具有典型双峰形状的时间-造影剂浓度曲线作为动脉输入函数。

将式(4)改写为矩阵形式：

$\left[\begin{array}{c}{C}_{\mathrm{corr}}\left(t_0\right)\\ C_{\mathrm{corr}}\left(t_1\right)\\ M\\ C_{\mathrm{corr}}\left(t_{N-1}\right)\end{array}\right]=\mathrm{CBF}\·\left[\begin{array}{cccc}{C}_a\left(t_0\right)& 0& \mathrm{\Λ}& 0\\ C_a\left(t_1\right)& C_a\left(t_0\right)& \mathrm{\Λ}& 0\\ M& M& O& M\\ C_a\left(t_{N-1}\right)& C_a\left(t_{N-2}\right)& \mathrm{\Λ}& C_a\left(t_0\right)\end{array}\right]\×\left[\begin{array}{c}R\left(t_0\right)\\ R\left(t_1\right)\\ M\\ R\left(t_{N-1}\right)\end{array}\right]---\left(5\right)$

式(5)中t₀为造影剂到达时刻；N为时间-造影剂浓度曲线的总长度；

式(5)可简写为矩阵形式：

c＝CBF·A·r      (6)

式(6)中c为时间-造影剂浓度序列；A为动脉输入函数矩阵；r为残余函数序列。r可表达为：

CBF·r＝A^-1·c    (7)

我们采用标准奇异值分解(sSVD)算法求取A^-1。首先利用奇异值分解算法将A分解为：

A＝U·S·V^T    (8)

则

$A^{-1}=V\·\frac{1}{S}\·U^T---\left(9\right)$

4.动脉输入函数的gamma函数拟合

动脉输入函数的首次通过的形状近似为gamma函数，可表达为：

$C_a^{1\mathrm{st}}\left(t\right)=\left\{\begin{array}{cc}0& t\≤t_0\\ C_0{(t-t_0)}^a{e}^{-(t-t_0)/b}& t>t_0\end{array}\right.---\left(10\right)$

式(10)中t₀为造影剂到达时刻；C_a ^1st(t)为动脉输入函数序列的首次通过；C₀，a和b为gamma函数的参数，决定了gamma函数的形状。利用最小二乘算法求得C₀，a和b的最优解，即可得到动脉输入函数的首次通过C_a ^1st(t)。

5.重建时间-造影剂浓度序列的首次通过

在求得残余函数R(t)及动脉输入函数的首次通过C_a ^1st(t)后，时间-造影剂浓度的首次通过可通过下式重建：

$C^{1\mathrm{st}}\left(t\right)=C_a^{1\mathrm{st}}\left(t\right)\⊗R\left(t\right)---\left(11\right)$

式(11)中C^1st(t)为时间-造影剂浓度的首次通过；R(t)为残余函数。为卷积运算符；

最后，利用重建的时间-造影剂浓度的首次通过来量化局部血容量(rCBV)：

$\mathrm{rCBV}\∝{\∫}_{t_0}^t{C}^{1\mathrm{st}}\left(\mathrm{\τ}\right)\mathrm{d\τ}---\left(12\right)$

式(12)中rCBV为局部血容量，t₀为造影剂到达时刻；C^1st(t)为时间-造影剂浓度的首次通过；τ为积分常数；

我们利用模拟实验及真是数据来测试本发明所提出的抑制磁共振灌注成像血容量量化误差的方法。

1.模拟数据

在模拟实验中，首先利用gamma函数来模拟动脉输入函数的首次通过，如式(10)所示。参数为C₀＝4.5，a＝3.0，b＝1.5，t₀＝15s(秒)。利用负指数函数来模拟残余函数 $R\left(t\right)=e^{-\frac{t}{\mathrm{MTT}}},$ 其中MTT为平均通过时间，取MTT＝4s。利用式(4)模拟时间-造影剂浓度序列的首次通过C^1st(t)。为了模拟时间-造影剂浓度序列中血液循环效应，引入血液循环传递函数 $H\left(t\right)=e^{-(t-T_0)/T_d}.$ 其中T₀为血液循环延迟时间，T_d为指数函数衰减系数。取T₀＝30s，T_d＝12s。带血液循环的时间-造影剂浓度序列通过下式来模拟：

$C\left(t\right)=C^{1\mathrm{st}}\left(t\right)+C^{1\mathrm{st}}\left(t\right)\⊗H\left(t\right)---\left(13\right)$

由式(13)模拟的C(t)经式(1)模拟出时间-灰度序列(取S₀＝1000)。利用均值为0，方差为S₀/SNR的正态分布序列来模拟真实图像数据中的噪声，其中SNR为信噪比。图2给出了模拟的时间-造影剂浓度曲线、其首次通过以及我们重建的首次通过的一个示例(SNR＝100)。从图中可以看出，用本发明方法所重建出的首次通过能够非常准确的贴近其真实形状。图3给出了当模拟的CBV在[0.05 0.40]之间变化时，用本发明的方法量化的CBV与其真实值之间的变化关系。从图中可以看出，当CBV在0.05～0.40之间变化时，从重建的首次通过量化的CBV非常接近其真实值。

2.脑胶质细胞瘤患者灌注成像数据

图4给出了一个脑胶质细胞瘤患者在接受放射疗法过程中的两次对比度增强结构扫描。其中，上一行属于第一次扫描，下一行属于第二次扫描，时间间隔为2个月。从图中可以看出，高亮区域的面积呈显著减小趋势，表明病灶区域没有呈现出复发或恶化。然而，由于造影剂渗漏及血液循化的影响，采用现有技术计算出的血容量图中病灶区域(图5a白色方框区域)的血容量要远高于其对侧区域(图5a黑色方框区域)的血容量值，暗示着病灶区域的复发或恶化，与前述结论相悖。

利用本发明中所述的抑制误差算法计算出的血容量图如图5b所示。图5b中，病灶区域(图5b白色方框区域)的血容量略小于其对侧区域(图5b黑色方框区域)的血容量图，暗示着病灶区域并不存在复发或者恶化，与根据对比度增强结构扫描得出的结论相符。脑胶质细胞瘤患者数据的试验结果表明，本发明所述的抑制磁共振灌注成像血流量量化误差的方法，能有效的减少治疗效果评价及预后过程中的误检。

Claims (2)

1.一种优化的磁共振灌注成像图像中血容量量化方法，其特征在于步骤如下：

步骤1：根据 $C(x,y,z,t)=-\frac{1}{\mathrm{TE}}\ln \left(\frac{S(x,y,z,t)}{S_0(x,y,z)}\right),$ 将图像中三维像素的时间-灰度信号曲线转换为时间-造影剂浓度曲线，其中：(x，y，z)为三维像素在图像矩阵中的索引值，C(x，y，z，t)为索引值(x，y，z)处三维像素的时间-造影剂浓度曲线，TE为磁共振灌注成像图像的回波时间，S(x，y，z，t)为索引值(x，y，z)处三维像素的时间-灰度信号曲线，S₀(x，y，z)为造影剂到达之前索引值(x，y，z)处三维像素的基准灰度值，ln(^·)为对数运算符；所述的基准灰度值为该像素时间-灰度信号曲线中造影剂达到之前的时间点的均值；

步骤2：根据造影剂渗漏效应矫正算法，矫正三维像素的时间-造影剂浓度曲线中因造影剂渗漏效应带来的附加顺磁效应，得到矫正后的时间-造影剂浓度曲线C_corr(t)；

步骤3：在图像中脑主中央动脉邻域选取一个像素，将该像素的时间-造影剂浓度曲线作为动脉输入函数f(t)；所述选取像素的时间-造影剂浓度曲线具有典型双峰形状；

步骤4：采用标准奇异值分解算法，以步骤3得到的f(t)为输入函数，对步骤2得到的C_corr(t)做去卷积运算，得到残余函数R(t)；

步骤5：采用gamma函数以非线性最小二乘算法拟合步骤3中得到的动脉输入函数f(t)，得到f(t)的首次通过f₁(t)；

步骤6：将步骤5中得到的f₁(t)与步骤4中得到的R(t)做卷积运算以重建步骤2中的C_corr(t)的首次通过分量C₁(t)；

步骤7：根据 $\mathrm{rCBV}\∝{\∫}_{t_0}^t{C}_1\left(\mathrm{\τ}\right)\mathrm{d\τ}$ ，利用步骤6中得到的C₁(t)来量化该像素点处的局部血容量rCBV，其中rCBV为局部血容量，t₀为造影剂到达时刻，τ为积分常数。

2.根据权利要求1所述的优化的磁共振灌注成像图像中血容量量化方法，其特征在于：步骤2中的造影剂泄漏效应校正算法采用文献“Boxerman，J.L.，K.M.Schmainda，and R.M.Weisskoff，Relative cerebral blood volume mapscorrected for contrast agent extravasation significantly correlate withglioma tumor grade，whereas uncorrected maps do not.AJNR Am J Neuroradiol，2006.27(4)：p.859-67.”中公布的算法。

Images