CN101912262B - 磁共振成像定量参数计算方法及系统 - Google Patents

Abstract

一种磁共振成像定量参数计算方法及系统，该方法包括以下步骤：建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型；根据磁共振成像模型获取动脉和病灶区的纵向弛豫时间、获取动脉输入函数、造影剂浓度时间函数、血管空间容积分数；根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数；根据所得的残余函数测得定量参数。上述磁共振成像定量参数计算方法及系统，动脉输入函数是直接针对病人自动提取的，比以往通过正常人提取或者手动提取动脉输入函数的准确性高，这样计算得到的定量参数的准确性也提高了。

Description

磁共振成像定量参数计算方法及系统

【技术领域】

本发明涉及磁共振成像数据处理领域，特别涉及一种磁共振成像定量参数计算方法及系统。

【背景技术】

肿瘤、斑块等新生血管生长密集且生长迅速，新生血管生成在肿瘤的发展转移及斑块的不稳定化过程中起了重要作用。血管内皮生长因子(vascularendothelial factor，简称VEGF)是新生血管生成的主要因子，它能促进内皮细胞增殖，提高血管通透性。微血管通透性和新生血管生成密切相关，是评价肿瘤恶性程度、鉴别其复发、以及评估斑块稳定性等的主要指标。因此，测得微血管通透性的定量指标，则可量化新生血管的生成，从而可用于评估肿瘤的恶性程度及斑块的稳定性。

微血管通透性是与新生血管生成速率密切相关的重要指标。血管(血浆)空间和血管外细胞外空间(Extravascular Extracellular Space，EES)之间的体积转运系数可用K^trans表示，它(微血管内皮通透率与表面积的乘积)可定量(量纲为1/minute)地表述微血管通透性。如果能准确测得K^trans值，定量地检测出微血管通透性信息，则可定量评价斑块、肿瘤等新生血管的生成情况，对于早期诊断和综合评估斑块的易损性以及活体对肿瘤进行组织病理学分级具有极大的帮助。然而，最常用的计算K^trans值的Tofts数据分析模型是通过对正常人群采样推演获得其用于建模的动脉输入函数(Arterial Input Function，AIF)，不能准确反映每一患者本身的AIF，影响了结果的准确性与可靠性，并且其基于药物动力学数据分析模型进行结果处理的非线性拟合方法将AIF假设为一个单指数阶跃响应函数，导致了残余函数估计的准确性降低。而目前采样患者本身的AIF的研究通过手动选择动脉的ROI，很容易受到操作者的偏差和部分容积效应的影响，影响了结果的准确性。因此，降低了定量指标K^trans用于具有新生血管生成的肿瘤、斑块等评估的鲁棒性和可重复性。

【发明内容】

基于此，有必要提供一种提高测得定量参数的准确性的磁共振成像定量参数计算方法。

一种磁共振成像定量参数计算方法，包括以下步骤：

建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型；

根据采集的造影剂注入前受检组织的多个翻转角图像，分别计算动脉和病灶区的纵向弛豫时间；

根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数；

根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数；

根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度时间函数计算血管空间容积分数；

根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数；

根据所得的残余函数测得定量参数。

优选地，所述计算动脉和病灶区的纵向弛豫时间具体步骤包括：

在动脉和最大病灶区选取感兴趣区；

分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。

优选地，所述获取动脉输入函数的具体步骤包括：

选取动脉层，获取动脉层的动态增强图像；

根据所述动态增强图像，提取所述动脉层内每个体素的信号强度增强与时间关系曲线；

滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线；

对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数。

优选地，所述获取病灶区感兴趣区内的造影剂浓度的具体步骤包括：

获取最大病灶区所在层的动态增强图像，并在所述最大病灶区层中的病灶增强部位选取感兴趣区；

分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线；

根据所述磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。

优选地，所述计算残余函数的方法为反卷积法和奇异值分解法。

此外，还有必要提供一种提高测得定量参数的准确性的磁共振成像定量参数计算系统。

一种磁共振成像定量参数计算系统，该系统包括：

模型构建模块，建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型；

信息采集模块，采集造影剂注入前受检组织的多个翻转角图像和造影剂注入期间受检组织的动态增强图像；

处理模块，根据采集的造影剂注入前受检组织的多个翻转角图像分别计算动脉和病灶区的纵向弛豫时间；根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数；根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数；根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度时间函数计算血管空间容积分数；根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数；根据所得的残余函数测得定量参数。

优选地，所述建立的包含信号强度与造影剂浓度之间关系的磁共振成像模型为：

$S\left(t\right)=M_o\frac{(1-e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))})\sin \left(\mathrm{\α}\right)}{1-\cos \left(\mathrm{\α}\right)e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))}}{e}^{-\mathrm{TE}(R_{20}*+r_2*C_t\left(t\right))}$

其中，S(t)是信号强度，M₀是反映磁共振成像扫描仪自旋密度和硬件灵敏度的常量，TR是重复时间，TE是回波时间，α是翻转角，R₁₀和R₂₀ ^*分别是造影剂注入前纵向弛豫率和有效横向弛豫率，r₁和r₂ ^*分别是造影剂的纵向弛豫性和有效横向弛豫性，C_t(t)是组织内随时间变化的造影剂浓度。

优选地，所述处理模块计算动脉和病灶区的纵向弛豫时间具体为：在动脉和最大病灶区选取感兴趣区，分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。

优选地，所述处理模块包括：

选取模块，选取动脉层，获取所述动脉层的动态增强图像；

提取模块，根据所述动态增强图像，提取所述动脉层内每个体素的信号强度增强与时间关系曲线；

滤除模块，滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线；

计算模块，对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数。

优选地，所述处理模块进一步获取最大病灶区所在层的动态增强图像，在所述最大病灶区层中的病灶增强部位选取感兴趣区，再分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线，以及根据所述磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。

优选地，所述处理模块进一步采用反卷积法和奇异值分解法计算残余函数。

优选地，所述处理模块进一步对所得的残余函数采用最小二乘拟合法进行处理测得定量参数。

上述磁共振成像定量参数计算方法及系统，通过采集对病人注入造影剂前的多个翻转角图像计算动脉和病灶区的纵向弛豫时间，以及采用注入造影剂后的动态增强图像及计算的动脉区的纵向弛豫时间及磁共振成像模型获取动脉输入函数，该动脉输入函数是直接针对病人自动提取的，比以往通过正常人提取或手动提取的动脉输入函数的准确性高，然后通过该动脉输入函数进行计算得到的血管空间容积分数及残余函数和定量参数的准确性也提高了。

另外，采用反卷积法提高了残余函数估计的准确性。

【附图说明】

图1为一个实施例中磁共振成像定量参数计算方法的流程图；

图2为动态对比增强药物动力学数据分析模型的两室模型示意图；

图3为一个重复时间内的梯度回波成像序列示意图；

图4为一个实施例中计算动脉和病灶区的纵向弛豫时间的具体方法流程图；

图5A为一个实施例中大脑中动脉ROI内的平均SFC及其拟合曲线；

图5B为一个实施例中病灶区ROI内的平均SFC及其拟合曲线；

图6为造影剂经动脉流入组织后，由静脉汇集流出的说明示意图；

图7为一个实施例中获取AIF的具体方法流程图；

图8A为一个实施例中动脉区内的平均信号强度增强与时间关系曲线；

图8B为一个实施例中AIF曲线；

图9为一个实施例中获取造影剂浓度时间函数的具体方法流程图；

图10A为病灶区ROI内平均信号强度与时间关系曲线；

图10B为病灶区ROI内平均造影剂浓度与时间关系曲线；

图11为AIF的单指数阶跃响应示意图；

图12为一个实施例中测得的残余函数

曲线及其拟合曲线；

图13为各定量参数的Map图；

图14为一个实施例中磁共振成像定量参数计算系统的结构示意图；

图15为一个实施例中处理模块的内部结构示意图。

【具体实施方式】

如图1所示，在一个实施例中，一种磁共振成像定量参数计算方法，包括以下步骤：

S10，建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型。其中，动态对比增强(DynamicContrast-enhanced，DCE)药物动力学数据分析模型为两室模型，如图2所示。该模型描述了造影剂在血管(血浆)空间(中央室)和血管外细胞外空间(Extravascular Extracellular Space，EES)之间转运的动力学过程，具体模型如下式(1)所示：

$V_e\frac{{\mathrm{dC}}_e\left(t\right)}{\mathrm{dt}}=K^{\mathrm{trans}}(C_p\left(t\right)-C_e\left(t\right))---\left(1\right)$

另外，血管组织内每一像素中观察到的随时间变化的造影剂浓度C_t(t)同时取决于血管(血浆)内的造影剂浓度和EES内的造影剂浓度，因此，求解一阶微分方程(1)，得到如下关系式(2)：

$C_t\left(t\right)=V_e\·C_e\left(t\right)+V_p\·C_p\left(t\right)$ $\left(2\right)$

$=K^{\mathrm{trans}}{\∫}_0^t{C}_p\left(\mathrm{\τ}\right)e^{-K_{\mathrm{ep}}(t-\mathrm{\τ})}\mathrm{d\τ}+V_p\·C_p\left(t\right)$

其中，EES容积分数V_e由下式(3)给出：

V_e＝K^trans/K_ep    (3)

对公式(2)进行处理可表示为如下卷积方程(4)：

式(4)中

是卷积操作符，其中，

上述(1)～(5)式中各参数分别为：V_p为血管(血浆)空间容积分数，V_e为EES容积分数，V_e和V_p均为每单位组织容积的无量纲分数；K^trans(min^-1)为造影剂从血管(血浆)空间渗漏到EES的体积转运系数；K_ep(min^-1)为造影剂从EES返回到血管(血浆)空间的速率常数；C_p(t)表示血管(血浆)空间的造影剂浓度随时间变化的关系，该参数可作为动脉输入函数(Arterial Input Function，AIF)测得；C_e(t)表示EES内的造影剂浓度随时间变化的关系；

为残余函数，表示在时间t时刻仍存在于组织中的造影剂浓度，是一个依赖组织特征的生理参数的关系式表达。

由于新生血管不完整，使得注入血管内的造影剂渗漏到周围组织，造影剂颗粒未成对电子与组织内原子核质子的偶极-偶极作用，导致组织的纵向弛豫时间T₁缩短，因此，造影剂的渗漏造成局部组织弛豫率的增加，而增加的弛豫率表现为T₁加权信号的增强。其中，T₁加权像的特点为：组织的T₁越短，恢复越快，信号就越强；组织的T₁越长，恢复越慢，信号就越弱。因此，磁共振成像(Magnetic Resonance Imaging，简称MRI)模型中给出了梯度回波成像序列的信号增强S(t)与造影剂浓度之间的关系如下：

$S\left(t\right)=M_o\frac{(1-e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))})\sin \left(\mathrm{\α}\right)}{1-\cos \left(\mathrm{\α}\right)e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))}}{e}^{-\mathrm{TE}(R_{20}*+r_2*C_t\left(t\right))}---\left(6\right)$

其中，M₀是反映MRI扫描仪自旋密度和硬件灵敏度的常量；图像采集的各个参数：TR是重复时间、TE是回波时间、α是翻转角，R₁₀和R₂₀ ^*分别是造影剂注入前纵向弛豫率和有效横向弛豫率(s^-1mM^-1)，r₁和r₂ ^*分别是造影剂的纵向弛豫性和有效横向弛豫性(s^-1mM^-1)，也就是每单位造影剂浓度弛豫率的增加，C_t(t)是组织中随时间变化的造影剂浓度。

该式(6)表明了由适当的弛豫性常数r₁和r₂ ^*线性量度的C_t(t)可以对弛豫率增强效应以线性比例进行描述，如式(7)给出了造影剂浓度随T₁变化的关系：

R₁(t)＝R₁₀+r₁C_t(t)    (7)

其中，R₁(t)＝1/T₁，造影剂注入以前的基线T₁值为T₁₀≡R₁₀。

步骤S20，根据采集的造影剂注入前受检组织的多个翻转角图像，分别计算动脉和病灶区的纵向弛豫时间。根据造影剂注入前采集的受检组织的多个翻转角图像，在Matlab平台下编写程序应用拟合法计算纵向弛豫时间T₁，即基线T₁值。本实施例中，受检组织为脑组织。

如图3所示为一个重复时间(Time Repetition，TR)内的梯度回波成像序列示意图，其中，α为翻转角，FID(Free Induction Decay)为自由感应衰减，TI为反转时间，TE为回波时间，TR为重复时间。根据其性质和MRI物理可知，在无造影剂注入的情况下，T₁加权信号和翻转角α的关系可表示为下式：

$S_0=M_0\·\sin \mathrm{\α}\·\left(\frac{1-e^{-\mathrm{TR}/T_{10}}}{1-\cos \mathrm{\α}\·e^{-\mathrm{TR}/T_{10}}}\right)---\left(8\right)$

其中，用T₁₀表示造影剂注入之前的纵向弛豫时间T₁，即，基线T₁；M₀是反映MRI扫描仪自旋密度和硬件灵敏度的常量，与平衡纵向磁化强度成比例；α(＝2°，5°，10°，15°，20°，......)为翻转角。

在一个实施例中，如图4所示，计算动脉和病灶区的纵向弛豫时间的具体步骤包括：

步骤S200，在动脉和最大病灶区选取感兴趣区。参照常规增强T₁加权图像，在动脉和最大病灶区所在层选取适当的感兴趣区(Region Of Interest，简称ROI)，且多个翻转角图像的ROI选取位置和大小均一致。本实施例中，采集5个不同翻转角(2°、5°、10°、15°、20°)的图像。

步骤S201，分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。计算出两个ROI内的平均信号强度，分别建立该信号强度与翻转角之间的关系曲线，根据式(8)，应用非线性最小二乘拟合方法可分别得到动脉和病灶区的纵向弛豫时间，即基线T₁值。如图5A和图5B所示为动脉区和病灶区ROI内的平均信号强度与翻转角的关系曲线(Signal-Flip angle Curve，SFC)以及用非线性最小二乘拟合方法得到的拟合曲线，图5A为大脑中动脉ROI内的平均SFC及其拟合曲线，图5B为病灶区ROI内的平均SFC及其拟合曲线。

步骤S30，根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数。如图6所示为造影剂经动脉流入组织后，由静脉汇集流出的说明示意图，其中，脉管下方的曲线为相应脉管内的造影剂浓度与时间关系曲线示意图。动脉血管内信号强度增强与时间关系曲线具有与AIF类似的特征：团注造影剂到达时间(bolus arrival time，BAT)短、峰值(peak height，PH)高、到达峰值的时间(time to peak，TTP)短、初始斜率大、平均高度高等，并且信号强度增强与造影剂浓度存在近似正比例的关系，因此，可以根据动脉血管内信号强度增强与时间关系曲线的特征自动提取AIF。在一个实施例中，如图7所示，获取AIF的具体步骤包括：

步骤S300，选取动脉层，获取该动脉层的动态增强图像。手动选取合适的动脉所在层，并保证该层没有受到流动效应的影响，获取该动脉层的动态增强图像。

步骤S301，根据该动态增强图像，提取该动脉层内每个体素的信号强度增强与时间关系曲线。对一系列动态增强图像进行处理，提取该动脉层内每个体素的信号强度增强-时间曲线，即信号强度增强与时间的关系曲线。

步骤S302，滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线。本实施例中，计算每个体素信号强度增强-时间曲线的PH、TTP和初始斜率，然后选出PH最高的前20％的所有曲线，接着再选出初始斜率最大的前50％的曲线，滤除非血管组织的体素的信号强度增强与时间关系曲线，其中，曲线选取的百分比根据需要而设定。然后，再对滤除非血管组织的体素的信号强度增强与时间关系曲线后剩下的关系曲线求平均TTP，从中选出TTP大于这一平均TTP的所有曲线，滤除静脉血管体素的信号强度增强与时间关系曲线。

步骤S303，对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数。对滤除了非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线剩下的关系曲线求平均，再根据式(6)和式(7)的MRI模型以及动脉区的纵向弛豫时间将其转换为造影剂浓度随时间变化的曲线，从而得到一个平均的AIF。如图8A和8B所示，图8A表示动脉区内的平均信号强度增强与时间关系曲线，图8B为AIF曲线。

步骤S40，根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数。根据团注造影剂期间采集到的动态增强图像，在Matlab平台下编写程序建立病灶区ROI内的造影剂浓度与时间的关系曲线(CTC)。在一个实施例中，如图9所示，步骤S40具体包括：

步骤S400，选取最大病灶区所在层，获取最大病灶区所在层的动态增强图像，并在该最大病灶区所在层中病灶增强部位选取感兴趣区。获取最大病灶区所在层的动态增强图像后，应用软件在该最大病灶区所在层选取合适的ROI，以保证选取到病灶增强的部位，避免由于ROI的选取部位及大小对动态曲线造成的影响，且在所有动态时间点选取的ROI的位置与大小均一致，同一患者的动态增强图像和多个翻转角图像的病灶区ROI选取位置及大小也相同，这样保证获取的数据来源相同，提高后续计算的准确性。

步骤S401，分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线。依次分别计算出所有动态增强图像中每个动态增强图像的ROI内所有像素的平均值，根据计算的多个动态点的平均值建立病灶区ROI内的平均信号强度与时间关系曲线(signal intensity-time curve，STC)。

步骤S402，根据磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。根据式(6)与式(7)的MRI模型及病灶区的纵向弛豫时间T₁，将病灶区ROI内的平均信号强度与时间关系曲线转换为平均造影剂浓度与时间关系曲线CTC，如图10A和10B所示分别为病灶区ROI内平均STC及相对应的CTC。

步骤S50，根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度计算血管空间容积分数。在团注造影剂后的起初20分钟内，血管(血浆)空间内造影剂浓度的变化可用单指数衰减来很好地表述。因此，可将AIF简化为如下式(9)的一个单指数阶跃响应函数，如图11所示，其中C_p(t)为AIF，C_t(t)为组织内造影剂浓度随时间变化的关系曲线，C_p(0)是团注造影剂即刻(t＝0)血浆内的浓度，此刻造影剂的浓度达到最大值，V_p为血管(血浆)空间的分数容积，T_1/2是造影剂的排泄半衰期。

$C_p\left(t\right)=C_p\left(0\right)\·e^{-t/T_{1/2}}---\left(9\right)$

其中，

V_pC_p(0)＝C_t(0)    (10)

根据式(9)计算的AIF，求得其峰值，再根据到达该峰值的时间以及CTC，计算出C_t(0)，则可计算出V_p。

步骤S60，根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数。对式(4)所示的卷积方程可通过数学上的矩阵运算变换为下式(11)：

基于动态对比增强药物动力学数据分析模型和以上测得的病灶区感兴趣区内的造影剂浓度时间函数C_t(t)、动脉输入函数C_p(t)和血管空间容积分数V_p，采用反卷积法和奇异值分解(Singular Value Decomposition，SVD)法求解式(11)，则可建立残余函数

的曲线。采用反卷积法提高了估算残余函数的准确性，比非线性拟合方法能得到更准确的计算结果。

步骤S70，根据所得的残余函数测得定量参数。根据测得的残余函数

的曲线，在Matlab平台下编写程序完成拟合，测得定量参数K^trans、K_ep、V_e等值。具体步骤为：对所得的残余函数进行最小二乘拟合法测得定量参数。以式(5)所示的残余函数为拟合模型，应用非线性最小二乘拟合方法的Levenberg-Marquardt算法对该残余函数

的曲线进行拟合，以拟合结果goodness中的决定系数rsquare大于预定值作为拟合终止标准。本实施例中，该预定值为0.7。如图12所示为残余函数

的曲线及其拟合曲线。

另外，上述动态增强磁共振成像定量参数计算方法，还包括计算每个像素的参数值，建立各参数值的Map图。根据计算各定量参数的方法计算出每个像素的参数值，然后建立各参数值的Map图，以直观的显示肿瘤、斑块等的形态大小。如图13所示为各定量参数的Map图。

上述磁共振成像定量参数计算方法，通过采集对病人注入造影剂前的多个翻转角图像计算动脉和病灶区的纵向弛豫时间，以及采用注入造影剂后的动态增强图像及计算的动脉区的纵向弛豫时间及磁共振成像模型获取的动脉输入函数，该动脉输入函数是直接针对病人自动提取的，比以往通过正常人提取或手动提取动脉输入函数的准确性高，因而通过该动脉输入函数进行计算得到的血管空间容积分数及残余函数和体积转运系数等定量参数的准确性也提高了。

如图14所示，在一个实施例中，一种磁共振成像定量参数计算系统，包括模型构建模块10、信息采集模块20和处理模块30。

模型构建模块10建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型。动态对比增强(DynamicContrast-enhanced，DCE)药物动力学数据分析模型与磁共振成像模型的具体模型同上述动态增强磁共振成像定量参数计算方法中描述的一样，此处不再详述。

信息采集模块20采集造影剂注入前受检组织的多个翻转角图像和造影剂注入期间受检组织的动态增强图像。本实施例中，受检组织为脑组织。信息采集模块20采集造影剂注入前的血管组织的多个翻转角图像，本实施例中，采集5个不同翻转角(2°、5°、10°、15°、20°)的图像。

处理模块30根据采集的造影剂注入前血管组织的多个翻转角图像分别计算动脉和病灶区的纵向弛豫时间。处理模块30计算动脉和病灶区的纵向弛豫时间T₁的具体过程为：

(a)在动脉和最大病灶区选取感兴趣区。参照常规增强T₁加权图像，在动脉和最大病灶区所在层选取适当的感兴趣区(Region Of Interest，简称ROI)，且多个翻转角图像的ROI选取位置和大小均一致。本实施例中，采集5个不同翻转角(2°、5°、10°、15°、20°)的图像。

(b)分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。计算出两个ROI内的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，根据式(8)，应用非线性最小二乘拟合方法可分别得到动脉和病灶区的纵向弛豫时间，即基线T₁值。

处理模块30进一步根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数。在一个实施例中，如图15，处理模块30包括选取模块300、提取模块301、滤除模块302和计算模块303。其中，选取模块300选取动脉层，获取该动脉层的动态增强图像。手动选取合适的动脉所在层，并保证该层没有受到流动效应的影响，获取该动脉层的动态增强图像。

提取模块301根据该动态增强图像，提取该动脉层内每个体素的信号强度增强与时间关系曲线。对动态增强图像进行处理，提取该动脉层内每个体素的信号增强-时间曲线，即信号强度增强与时间关系曲线。

滤除模块302滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线。本实施例中，计算每个体素信号强度增强-时间曲线的PH、TTP和初始斜率，然后选出PH最高的前20％的所有曲线，接着再选出初始斜率最大的前50％的曲线，滤除非血管组织的体素的信号强度增强与时间关系曲线，其中，曲线选取的百分比根据需要而设定。然后，再对滤除非血管组织的体素的信号强度增强与时间关系曲线后剩下的关系曲线求平均TTP，从中选出TTP大于这一平均TTP的所有曲线，滤除静脉血管体素的信号强度增强与时间关系曲线。

计算模块303对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数。对滤除了非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线剩下的关系曲线求平均，再根据式(6)和式(7)的MRI模型以及动脉区的纵向弛豫时间将其转换为造影剂浓度随时间变化的曲线，从而得到一个平均的AIF。

处理模块30还根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数。具体过程为：

(a₁)获取最大病灶区所在层的动态增强图像，并在该最大病灶区所在层中的病灶增强部位选取感兴趣区。获取最大病灶区所在层的动态增强图像后，应用软件在该最大病灶区所在层选取合适的ROI，以保证选取到病灶增强的部位，避免由于ROI的选取部位及大小对动态曲线造成的影响，且在所有动态时间点选取的ROI的位置与大小均一致，同一患者的动态增强图像和多个翻转角图像的病灶区ROI选取位置及大小也相同，这样保证获取的数据来源相同，提高后续计算的准确性。该过程可以处理模块30的选取模块300实现。

(a₂)分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线。依次分别计算出所有动态增强图像中每个动态增强图像的ROI内所有像素的平均值，根据计算的多个动态点的平均值建立病灶区ROI内的平均信号强度与时间关系曲线(signal intensity-time curve，STC)。

(a₃)根据磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。根据式(6)与式(7)的MRI模型及病灶区的纵向弛豫时间T₁，将病灶区ROI内的平均信号强度与时间关系曲线转换为造影剂浓度与时间关系曲线CTC，如图10A和10B所示分别为病灶区ROI内平均STC曲线及相对应的CTC曲线。其中，(2)和(3)均可由处理模块30的计算模块303进行处理。

处理模块30还根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度计算血管空间容积分数，根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度、动脉输入函数及血管空间容积分数计算残余函数，以及根据所得的残余函数测得定量参数。该定量参数包括体积转运系数K^trans、造影剂从EES返回到血管(血浆)空间的速率常数K_ep、EES容积分数V_e等，具体过程是：处理模块30对动脉输入函数AIF处理为一个单指数阶跃响应函数简化计算，得到

$C_p\left(t\right)=C_p\left(0\right)\·e^{-t/T_{1/2}}---\left(9\right)$

其中，

V_pC_p(0)＝C_t(0)    (10)

其中C_p(t)为AIF，C_t(t)为组织内造影剂浓度随时间变化的关系曲线，C_p(0)是团注造影剂即刻(t＝0)血浆内的浓度，此刻造影剂的浓度达到最大值，V_p为血管(血浆)空间的分数容积，T_1/2是造影剂的排泄半衰期。

根据式(9)计算的AIF，求得其峰值，再根据到达该峰值的时间以及CTC，计算出C_t(0)，则可计算出V_p。

处理模块30对式(4)所示的卷积方程可通过数学上的矩阵运算变换为式(11)：

基于动态对比增强药物动力学数据分析模型和以上测得的病灶区感兴趣区内的造影剂浓度时间函数C_t(t)、动脉输入函数C_p(t)和血管空间容积分数V_p，处理模块30采用反卷积法和奇异值分解(Singular Value Decomposition，SVD)法求解式(11)，则可建立残余函数

的曲线。采用反卷积法提高了估算残余函数的准确性，比非线性拟合方法能得到更准确的计算结果。处理模块30再进一步对所得的残余函数采用最小二乘拟合法进行处理测得定量参数。以式(5)所示的残余函数为拟合模型，应用非线性最小二乘拟合方法的Levenberg-Marquardt算法对该残余函数

的曲线进行拟合，以拟合结果goodness中的决定系数rsquare大于预定值作为拟合终止标准。本实施例中，该预定值为0.7。

另外，处理模块30还计算每个像素的参数值，建立各参数值的Map图。根据计算各定量参数的方法计算出每个像素的参数值，然后建立各参数值的Map图，以直观的显示肿瘤、斑块等的形态大小。

上述磁共振成像定量参数计算方法及系统，通过实验得到了较好的效果。我们利用该方法分析了28例脑胶质瘤患者和2位正常人的数据，脑胶质瘤患者中男18例，女10例，年龄19-74岁，平均47.11±14.18岁。其中低级别(I级8例、II级6例)和高级别(III级6例、IV级8例)均为14例。所有实验数据均利用1.5T超导型MRI系统(Syngo MR 2002B，SIEMENS)和常规头部线圈采集。

患者均采用仰卧位头先进。对病灶进行定位后，先是对患者进行常规非增强MRI扫描，采集横断面、冠状或矢状面T₁加权图像和横断面T₂加权图像。T₁加权图像采用自旋回波(SE)序列，T₂加权图像采用快速自旋回波(TSE)序列，成像参数分别为TR/TE：450ms/10ms(T1)、4200ms/98ms(T2)，层厚5mm，间隔0.5mm，FOV：211mm×240mm，图像矩阵：256×208。

然后，在造影剂注入前，应用2D Turbo FLASH(fast low angle shot)成像序列采集2°、5°、10°、15°、20°等5个不同翻转角的T₁加权图像各40幅，对于每个翻转角，以最大病灶区所在层为中心采集10层图像，层厚3mm，层间距1mm，共采集4次，成像参数分别为TR/TE：199ms/1.05ms，FOV：211mm×260mm，图像矩阵：256像素×208像素。

动态组图像的采集在20°翻转角的条件下，采用2D Turbo FLASH扫描序列重复前面的采集过程，即采集病灶处10层图像，重复采集90组，每组采集时间为4s，共采集900幅DCE T₁加权图像。在采集过程中，从第6组采集开始用高压注射器经肘前静脉团注顺磁性造影剂药团Gd-DTPA(0.1mmol/kg，广州康臣药业有限公司)，注射速率为4ml/s，造影剂注射完后即刻以同样速度再注入同等容积生理盐水冲洗导管，以减少团注后造影剂的残留。成像参数与5个翻转角的T₁加权图像相同。计算结果如表1所示，四个病理分级与正常人的各定量参数的均值±标准偏差，0*表示为正常人，1～4表示为病人。本实施例中，采用SPSS11.5统计软件对所有患者的结果数据进行统计分析，对于所有参数，脑胶质瘤四个病理分级之间的差异比较采用多个独立样本比较的Kruskal-Wallis H检验，对于存在统计学差异的参数，再采用两独立样本比较的Mann-Whitney U检验进行四个病理分级之间的两两比较，P＜0.05作为所有参数的统计学差异标准。比较结果分别见表2和表3。从初步实验结果与病理结果比较可知，通过本发明方法测得的与微血管通透性相关的定量参数K^trans、V_e可用于手术或治疗前准确地对脑胶质瘤进行无创性的病理分级，从而为医生们提供一些精确的病理生理学参考信息，以便其选择更好的治疗(包括手术)方案，以提高治疗的成功率，达到提高患者的生存率、改善其生存质量的目的。

表1

表2

表3

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种磁共振成像定量参数计算方法，包括以下步骤：

建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型；

根据采集的造影剂注入前受检组织的多个翻转角图像，分别计算动脉和病灶区的纵向弛豫时间；

根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数；

根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数；

根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度时间函数计算血管空间容积分数；

根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数；

根据所得的残余函数测得定量参数；

所述获取动脉输入函数的具体步骤包括：

选取动脉层，获取动脉层的动态增强图像；

根据所述动脉层的动态增强图像，提取所述动脉层内每个体素的信号强度增强与时间关系曲线；

滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线；

对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数；

其中，所述滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线的步骤具体为：计算每个体素信号强度增强-时间曲线的峰值、到达峰值的时间和初始斜率，然后选出峰值最高的前20％的所有曲线，接着再选出初始斜率最大的前50％的曲线，滤除非血管组织的体素的信号强度增强与时间关系曲线，再对滤除非血管组织的体素的信号强度增强与时间关系曲线后剩下的关系曲线求平均到达峰值的时间，从中选出到达峰值的时间大于这一平均到达峰值的时间的所有曲线，滤除静脉血管体素的信号强度增强与时间关系曲线。

2.根据权利要求1所述的磁共振成像定量参数计算方法，其特征在于，所述计算动脉和病灶区的纵向弛豫时间具体步骤包括：

在动脉和最大病灶区选取感兴趣区；

分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。

3.根据权利要求1至2中任一项所述的磁共振成像定量参数计算方法，其特征在于，所述获取病灶区感兴趣区内的造影剂浓度时间函数的具体步骤包括：

选取最大病灶区所在层，获取最大病灶区层的动态增强图像，并在所述最大病灶区层中的病灶增强部位选取感兴趣区；

分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线；

根据所述磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。

4.根据权利要求1所述的磁共振成像定量参数计算方法，其特征在于，所述计算残余函数的方法为反卷积法和奇异值分解法。

5.一种磁共振成像定量参数计算系统，其特征在于，该系统包括：

模型构建模块，建立动态对比增强药物动力学数据分析模型及包含信号强度与造影剂浓度之间关系的磁共振成像模型；

信息采集模块，采集造影剂注入前受检组织的多个翻转角图像和造影剂注入期间受检组织的动态增强图像；

处理模块，根据采集的造影剂注入前受检组织的多个翻转角图像分别计算动脉和病灶区的纵向弛豫时间；根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及动脉区的纵向弛豫时间获取动脉输入函数；根据采集的造影剂注入期间受检组织的动态增强图像、磁共振成像模型及病灶区的纵向弛豫时间获取病灶区感兴趣区内的造影剂浓度时间函数；根据获取的动脉输入函数及病灶区感兴趣区内的造影剂浓度时间函数计算血管空间容积分数；根据动态对比增强药物动力学数据分析模型及获取的病灶区感兴趣区内的造影剂浓度时间函数、动脉输入函数及血管空间容积分数计算残余函数；根据所得的残余函数测得定量参数；

所述处理模块包括：

选取模块，选取动脉层，获取所述动脉层的动态增强图像；

提取模块，根据所述动脉层的动态增强图像，提取所述动脉层内每个体素的信号强度增强与时间关系曲线；

滤除模块，滤除非血管组织的体素及静脉血管的体素的信号强度增强与时间关系曲线；所述滤除模块还用于计算每个体素信号强度增强-时间曲线的峰值、到达峰值的时间和初始斜率，然后选出峰值最高的前20％的所有曲线，接着再选出初始斜率最大的前50％的曲线，滤除非血管组织的体素的信号强度增强与时间关系曲线，再对滤除非血管组织的体素的信号强度增强与时间关系曲线后剩下的关系曲线求平均到达峰值的时间，从中选出到达峰值的时间大于这一平均到达峰值的时间的所有曲线，滤除静脉血管体素的信号强度增强与时间关系曲线；

计算模块，对滤除后的信号强度增强与时间关系曲线求平均，再根据磁共振成像模型及动脉区的纵向弛豫时间获取平均的动脉输入函数。

6.根据权利要求5所述的磁共振成像定量参数计算系统，其特征在于，所述建立的包含信号强度与造影剂浓度之间关系的磁共振成像模型为：

$S\left(t\right)=M_0\frac{(1-e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))})\sin \left(\mathrm{\α}\right)}{1-\cos \left(\mathrm{\α}\right)e^{-\mathrm{TR}(R_{10}+r_1{C}_t\left(t\right))}}{e}^{-\mathrm{TE}(R_{20}*+r_2*C_t\left(t\right))}$

其中，S(t)是信号强度，M₀是反映磁共振成像扫描仪自旋密度和硬件灵敏度的常量，TR是重复时间，TE是回波时间，α是翻转角，R₁₀和R₂₀ ^*分别是造影剂注入前纵向弛豫率和有效横向弛豫率，r₁和r₂ ^*分别是造影剂的纵向弛豫性和有效横向弛豫性，C_t(t)是组织内随时间变化的造影剂浓度。

7.根据权利要求5所述的磁共振成像定量参数计算系统，其特征在于，所述处理模块计算动脉和病灶区的纵向弛豫时间具体为：在动脉和最大病灶区选取感兴趣区，分别计算动脉感兴趣区和病灶区感兴趣区的平均信号强度，建立该平均信号强度与翻转角之间的关系曲线，再分别计算得到动脉区和病灶区的纵向弛豫时间。

8.根据权利要求5至7中任一项所述的磁共振成像定量参数计算系统，其特征在于，

所述处理模块进一步获取最大病灶区所在层的动态增强图像，在所述最大病灶区所在层中的病灶增强部位选取感兴趣区，再分别计算所有动态增强图像中每个动态增强图像的感兴趣区内所有像素的平均值，并根据计算得出的多个动态点的平均值建立病灶区感兴趣区的平均信号强度与时间关系曲线，以及根据所述磁共振成像模型、病灶区的纵向弛豫时间及平均信号强度与时间关系曲线建立造影剂浓度与时间关系曲线。

9.根据权利要求5所述的磁共振成像定量参数计算系统，其特征在于，所述处理模块进一步采用反卷积法和奇异值分解法计算残余函数。

10.根据权利要求5所述的磁共振成像定量参数计算系统，其特征在于，所述处理模块进一步对所得的残余函数采用最小二乘拟合法进行处理测得定量参数。

Images