基于动态对比增强MRI的生理参数定量统计优化方法与流程

本发明涉及到一种基于动态对比增强mri的生理参数定量统计优化方法。

背景技术

动态磁敏感对比增强磁共振成像(简称dsc-mri)是评估血管内指标动力学的重要工具，这个过程通常包括将γ变量函数拟合到观察的浓度-时间曲线。目前，传统上有两种广泛的处理方法。一种是非线性优化方法，表现出高精度，但通常在计算上是很繁重的并且需要可靠的初始值以保证成功，而另一种是对数线性最小二乘法(ll-ls)方法表现得更加稳定和高效，并且不受初始值问题的影响，但其使数据统计复杂化，并且意味着噪声不再符合高斯分布，估计时会产生偏差，特别是当存在少量数据或异常值时。

随着研究的推进，需要在整个大脑中以更精细的体素比例来确定局部脑血容量(简称cbv)，以便更准确评估血液动力学。但因此，采样间隔就会相对较长，在曲线拟合的第一阶段只收集5或6个有效数据点，而由于现有的两种方法表现不佳，因此有必要开发一种结合两种方法优点的替代方法，即使在数据点很少的情况下仍然表现出令人满意的性能。

技术实现要素：

本发明要克服现有技术的上述缺点，提出一种基于动态对比增强mri的生理参数定量统计优化方法。

在实际情况中，所采集到的信号都是有很大的噪声，且有效采集点也会较少。因此，我们提出一种新的统计优化算法，通过将起初的曲线拟合问题转化为γ概率分布拟合问题，就是将浓度-时间曲线看成一个以时间为随机变量的γ分布的随机样本。然后通过最大似然估计(mle)解出最优估计量，从而解决上述问题。

本发明是一种基于动态对比增强mri的生物参数定量统计优化方法，具体步骤如下：

步骤1，建立信号强度和时间的关系模型：

其中ci(t)是造影剂浓度，k是未知的比例系数，te是成像序列的回波时间，s(t)是时间点t处的mri信号强度，s0是施用造影剂之前的初始mri信号强度。

步骤2，将信号强度和时间关系转化为浓度和时间关系式。根据稀释理论，相对脑血容量(rcbv)与浓度-时间曲线(δr2^*(t))下面积成正比，因此rcbv公式如下：

其中s(t)/s0表示参考体素中的相对信号衰减，脑组织的密度ρ＝1.04g/ml，校正因子kh是指感兴趣的体素和参考体素之间的血细胞比容差异，一般取1。

步骤3，考虑实际情况中信号传递延迟和残余造影剂的影响，用γ-函数对曲线建模：

其中t0是指定区域应用造影剂的时间，a，α和β是确定函数形状的参数。

步骤4，计算rcbv的值。由于rcbv正比于浓度-时间曲线所围成的面积，故将浓度-时间积分可得：

步骤5，将浓度-时间关系重新构建，表示成γ概率分布模型：

其中，arcbv表示rcbv的值，是γ分布，t>0，α>0，β>0。

步骤6，结合γ分布，运用mle对模型进行估计，建立似然函数。如下：

其中y＝(y(t1),y(t2),…,y(tk),…)，k＝1,...,n，对应在独立同分布离散随机样本集x＝(x1,x2,...,xn)中n个不同取值t1,t2,...,tn发生的概率，f(y(ti)|α,β)表示α，β未知情况下的γ概率密度函数。

步骤7，对公式(6)进行mle估计，得到公式(7)、(8)，化简并联立方程得公式(9)。可得：

其中

步骤8，根据thom提出的关于γ函数的渐近线展开的方法，对公式(7)中伽马函数展开：

其中，bk是伯努利数，b1＝1/6,b2＝1/30,…，rm是m项后的余项，而当α≥1时，rm就可以被忽略，α越大，误差就越小。

步骤9，取m＝1，将ψ(α)展开式带入mle估计方程并解出估计值：

其中，加下缀以区别上述的符号。

步骤10，结合和公式(5)、(8)以及(12)，就可推导得出相关待定系数：

通过获得的离散时间数据点，就可确定相关函数形状参数，从而得到估计模型。

步骤11，将待定系数带入公式(4)，得出模型并作图。可得：

步骤12，收集数据点，代入本发明所提出的优化方法，建立估计模型。

步骤13，定义参数变量和评价因子。

步骤14，根据评价因子和参数变量对三种方法进行评价。

本发明的有益效果是：无论信噪比的高低和时间分辨率的优劣，本发明方法始终表现出很优异的估计性能，并且估计结果变现的很稳定，灵敏度也很高；本发明通过将浓度-时间曲线转变为随时间变化的γ分布的随机样本，本发明将原始曲线拟合问题转化为γ分布统计估计问题，然后采用mle来解决这个问题，因此，本发明不受初始值的影响，并且考虑了整个样本所有观测值，产生的误差很小，可以提供更稳定和准确的曲线分析，十分适合低信噪比，时间分辨率差的情况下进行曲线估计分析。

附图说明

图1为信号强度和时间的关系图。

图2为浓度和时间关系图。

图3为浓度-时间的γ概率分布模型。

图4为ll-ls方法、非线性方法的和本发明方法的对比图。

图5为本发明方法μ值的统计估计图。

图6为本发明方法σ值的统计估计图。

图7为ll-ls方法μ值的拟合估计图。

图8为ll-ls方法σ值的拟合估计图。

图9为非线性方法μ值的拟合估计图。

图10为非线性方法σ值的拟合估计图。

图11是本发明方法的流程图。

具体实现方式

下面结合附图进一步说明本发明方法。

本发明是一种基于动态对比增强mri的生物参数定量统计优化方法，为验证所提出方法的有效性和优越性，基于浓度-时间曲线，并在不同信噪比和时间分辨率下比较三种拟合算法的性能。具体实现方式如下:

步骤1，建立信号强度和时间的关系模型，(如图1)。关系模型如下：

其中ci(t)是造影剂浓度，k是未知的比例系数，te是成像序列的回波时间，s(t)是时间点t处的mri信号强度，s0是施用造影剂之前的初始mri信号强度。

步骤2，将信号强度和时间关系转化为浓度和时间关系式，(如图2)。根据稀释理论，相对脑血容量(rcbv)与浓度-时间曲线(δr2^*(t))下面积成正比，因此rcbv公式如下：

其中s(t)/s0表示参考体素中的相对信号衰减，脑组织的密度ρ＝1.04g/ml，校正因子kh是指感兴趣的体素和参考体素之间的血细胞比容差异，一般取1。

步骤3，考虑实际情况中信号传递延迟和残余造影剂的影响，用γ-函数对曲线建模：

其中t0是指定区域应用造影剂的时间，a，α和β是确定函数形状的参数。

步骤4，计算rcbv的值。由于rcbv正比于浓度-时间曲线所围成的面积，故将浓度-时间积分可得：

步骤5，将浓度-时间关系重新构建，表示成γ概率分布模型，(如图3)。模型如下：

其中，arcbv表示rcbv的值，是γ分布，t>0，α>0，β>0。

步骤6，结合γ分布，运用mle对模型进行估计，建立似然函数。如下：

其中y＝(y(t1),y(t2),…,y(tk),…)，k＝1,...,n，对应在独立同分布离散随机样本集x＝(x1,x2,...,xn)中n个不同取值t1,t2,...,tn发生的概率，f(y(ti)|α,β)表示α，β未知情况下的γ概率密度函数。

步骤7，对公式(6)进行mle估计，得到公式(7)、(8)，化简并联立方程得公式(9)。可得：

其中

步骤8，根据thom提出的关于γ函数的渐近线展开的方法，对公式(7)中伽马函数展开：

其中，bk是伯努利数，b1＝1/6,b2＝1/30,...，rm是m项后的余项，而当α≥1时，rm就可以被忽略，α越大，误差就越小。

步骤9，取m＝1，将ψ(α)展开式带入mle估计方程并解出估计值：

其中，加下缀以区别上述的符号。

步骤10，结合和公式(5)、(8)以及(12)，就可推导得出相关待定系数：

通过获得的离散时间数据点，就可确定相关函数形状参数，从而得到估计模型。

步骤11，将待定系数代入公式(4)，得出模型。可得：

步骤12，收集数据点，代入公式(14)，建立估计模型，同时也将ll-ss和非线性模型作图(如图4)。用1.5-tesla扫描仪获取图像，cbv成像序列由30个数据节点构成。

步骤13，定义参数变量和评价因子。所选的参数变量和评价因子如下：

其中ymax表示曲线的最大值，σ是指噪声的标准差。时间分辨率δt取值范围为[0.2，3.2]，步长为0.2s。

其中是指250个生成曲线的rcbv的均值，fi是指第i个估计曲线的rcbv。

步骤14，根据评价因子和参数变量对三种方法进行评价(如图5-10)。通过对比图5-10，ll-ls方法和非线性方法在snr(即有效信号)很低的时候，随着时间分辨率的增加，μ和σ的波动范围都会超过50％(当估计参数不确定性超过50％，一般认为拟合失败)，拟合失败。这两种方法在缺乏可靠数据，拟合性能很差。而无论信噪比的高低和时间分辨率的优劣，本发明方法始终表现出很优异的估计性能，并且估计结果变现的很稳定，灵敏度也很高。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。