基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法及装置和应用

技术领域

本发明涉及磁共振成像的技术领域，具体涉及基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法及装置和应用。

背景技术

动态对比增强磁共振成像(dynamic contrast-enhanced magnetic resonanceimaging,DCE-MRI)技术，通过将时间序列MRI信号拟合到示踪剂动力学模型，能提供大量的与临床诊断息息相关的病理信息。这种技术在临床研究和诊断中得到了广泛的应用。然而，从DCE-MRI数据中估算出的药代动力学参数受到多种因素的影响，如信噪比(SNR)、原生T1、起始时间、动脉输入函数和拟合算法等。例如，广泛使用的线性示踪剂动力学模型之一的Patlak模型，其拟合结果很容易受到信噪比的影响。此外，某些非线性示踪剂动力学模型，如Tofts模型和扩展Tofts模型(eTofts)，具有非凸性质，这进一步增加了问题的复杂性。具体来说，可能存在多个的药代动力学参数组合，这些参数组合可能在拟合测量到的时间序列MRI信号上具有相似的拟合残差，这进一步增加了估计结果的不确定性。

为了解决这些问题，研究者提出了多种药代动力学参数估计策略。其中，非概率方法采用优化算法来寻找最优的药代动力学参数，使观测到的DCE-MRI数据与示踪剂动力学模型模型相拟合，例如最小二乘估计方法。然而，这些药代动力学参数估计结果可能表现出大的方差，并且不提供关于不确定性的信息。贝叶斯估计方法是一种经典的概率方法，已在DCE-MRI数据的后处理中被使用，并在一致性和准确性方面显示出出色的结果。药代动力学参数的后验分布可以指示真实参数在参数空间任意位置出现的概率，其中标准差是参数估计不确定性的描述符。根据贝叶斯定理，后验分布的直接计算需要一个观测模型和一个先验分布，这些条件大部分是未知的，不适当的观测模型或先验假设可能会降低估计的准确性。此外，贝叶斯规则中积分的数值计算也是比较棘手的问题。

估计药代动力学参数的后验分布的另一种方法是通过一类预定义的分布来近似后验分布，该类分布具有固定数量的未知参数(例如，均值和偏差)，并使用期望最大化算法来求解分布参数的解析解(KELM B M,MENZE B H,NIX O,等.Estimating KineticParameter Maps From Dynamic Contrast-Enhanced MRI Using Spatial PriorKnowledge[J/OL].IEEE Transactions on Medical Imaging,2009,28(10):1534-1547.https://doi.org/10.1109/TMI.2009.2019957)。此外，神经网络模型也可以被用来缓解求解过程中复杂的计算问题，并具有同时估计药代动力学参数后验分布的均值和标准差的能力(BLIESENER Y,ACHARYA J,NAYAK K S.Efficient DCE-MRIParameter andUncertainty Estimation Using a Neural Network[J/OL].IEEE Transactions onMedical Imaging,2020,39(5):1712-1723.)。然而，如已有的文献所述(ORTON M R,COLLINSD J,WALKER-SAMUEL S,等.Bayesian estimation of pharmacokineticparameters for DCE-MRI with a robust treatment of enhancement onset time[J/OL].Physics in Medicine and Biology,2007,52(9):2393-2408.和BLIESENER Y,ACHARYA J,NAYAK K S.Efficient DCE-MRIParameter and Uncertainty EstimationUsing a Neural Network[J/OL].IEEE Transactions on Medical Imaging,2020,39(5):1712-1723.)，真实的后验分布与假定的后验分布之间的差距可能会影响药代动力学参数的估计准确性，从而使得估计出的标准差(不确定性)也可能是不准确的。当假设药代动力学参数遵循更一般的未知分布时，为了准确估计药代动力学参数，更精确地近似后验分布是非常必要的。

因此，开发一种高效且灵活的方法，用于更好得近似和估计药代动力学参数的真实后验分布，被认为是至关重要的。

发明内容

本发明的目的在于提供基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法及装置，该方法及装置可以提高药代动力学参数分布的准确估计。

本发明采取以下技术方案：

基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法，所述方法包括以下步骤：

(1)通过采集定量T1图像数据和DCE-MRI数据，计算血管内造影剂浓度变化曲线C

(2)设定所要估计的药代动力学模型参数的范围，并在所设定的范围内进行均匀随机采样，生成药代动力学参数组；根据C

(3)构建基于流模型的参数分布估计模型FPDEN：包括参数回归网络和流模型网络，将

(4)采集定量T1图像数据和DCE-MRI数据，并计算血管内造影剂浓度曲线和对应的组织内造影剂浓度曲线，将其作为FPDEN中参数回归模型的输入，输出对应的药代动力学参数的均值和标准差。

在本发明中，所述药代动力学模型可以为eTofts和Patlak等多种药代动力学模型。本发明提供的方法在多种药代动力学模型上具有通用性，可根据需要任选。

进一步地，在步骤(2)中，设定药代动力学模型参数的范围为K

进一步地，在步骤(2)中，从C

使用eTtofts模型作为药代动力学模型，生成仿真DCE-MRI信号，具体为：

eTofts模型的信号由公式

进一步地，在步骤(2)中，对仿真DCE-MRI信号加入莱斯噪声。

在步骤(2)中，所述参数回归网络估计参数分布均值和标准差，流模型网络用于学习PK参数(药代动力学参数)的分布形状，并结合回归网络估计的均值和标准差计算给定采样值在所预测分布上的概率值。

进一步地，在步骤(3)中，所述参数回归网络包括一维卷积神经网络CNN层、双向长短时记忆网络Bi-LSTM层、平均层和全连接层；将两个通道的数据

进一步地，在步骤(3)中，所述流模型网络使用了改编后的RealNVP网络结构，包括6个仿射耦合层。

进一步地，在步骤(3)中，所述参数回归网络和流模型网络的训练使用最大似然估计loss函数，训练目标是使得给定的PK参数在预测分布上的log似然值最大。

所述步骤(3)具体为：将步骤(2)生成的组织内造影剂浓度曲线和对应的血管内造影剂浓度曲线作为参数回归网络的输入，输出均值和标准差，该曲线对应的PK参数作为一个从参数真实分布上的采样值，使用输出的均值和标准差对该采样值进行归一化，然后输入到流模型网络中，计算其log似然值。

在步骤(3)中，模型的训练使用最大似然估计loss函数，模型的训练目标是使得给定的PK参数在预测分布上的log似然值最大。网络训练到收敛后，得到优化的模型参数。

所述步骤(4)具体为：使用采集得到的T1定量图像和采集到的对应的DCE-MRI信号，计算血管内造影剂浓度曲线。然后，逐个像素计算组织的造影剂浓度曲线。对于每个像素，将血管内造影剂浓度曲线和组织造影剂浓度曲线作为步骤(4)获得的回归模型的输入，输出对应的药代动力学参数的均值和标准差。

本发明还提供了一种上述基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法在胶质瘤分级上的应用。

进一步地，所述应用为非疾病诊断为目的的胶质瘤分级上的应用。

胶质瘤分级可用于科学研究，如通过揭示不同级别胶质瘤在基因表达和信号通路活性方面的差异，帮助研究人员进行更精细的患者分层和数据标准化，进而在药物开发、治疗策略优化、资源分配和多学科交叉研究等方面提供有力的支持。

进一步地，所估计药代动力学参数分布的均值作为待估计参数实际估计值，分布的标准差作为待估计参数的不确定性指标，用于过滤不可靠的实际估计值，过滤后的实际估计值应用在胶质瘤分级上。

在本发明中，所估计的药代动力学参数的分布包含两个指标，均值和标准差；其中，使用均值作为待估计参数的实际值，进行胶质瘤分级，而标准差就是不确定性，用于甄别所估计的参数是否可靠，不确定性较大的估计值会被排除从而不参与胶质瘤分级的计算，从而提高使用所估计药代动力学参数作为胶质瘤分级依据时的准确性。

进一步地，所述胶质瘤分级为区分低级与高级胶质瘤、区分III级与IV级胶质瘤。

本发明还提供了一种基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计装置，包括：

图像采集及处理模块，通过采集定量T1图像数据和DCE-MRI数据，计算血管内造影剂浓度变化曲线C

数据集构建模块，设定所要估计的药代动力学模型参数的范围，并在所设定的范围内进行均匀随机采样，生成药代动力学参数组；根据C

模型构建及训练模块，建基于流模型的参数分布估计模型FPDEN：包括参数回归网络和流模型网络，将

检测模块，采集定量T1图像数据和DCE-MRI数据，并计算血管内造影剂浓度曲线和对应的组织内造影剂浓度曲线，将其作为FPDEN中参数回归模型的输入，输出对应的药代动力学参数的均值和标准差。

本发明还提供了计算设备，包括存储器和一个或多个处理器，所述存储器中存储有可执行代码，所述一个或多个处理器执行所述可执行代码时，用于实现上述基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法。

本发明还提供了一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时，用于实现上述基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法。

现有技术中，PK参数的估计方法存在多种不确定性和可变性，这影响了其在临床研究和诊断中的应用效果。为了解决这一问题，本发明提供了基于流模型的DCR-MRI药代动力学参数分布估计方法及装置，所采用的基于流模型的参数分布估计模型FPDEN，用于自适应地学习和近似药代动力学(PK)参数的后验分布：本发明提供的方法中采用了流模型，该模型能够从数据集中学习一系列可逆转换，将复杂的参数分布转化为具有解析解的简单分布，从而实现PK参数分布的准确估计。此外，本发明还引入了重新参数化策略，以提高参数分布估计模型的可行性。

因此，本发明具有以下优点：

(1)高度准确的参数估计：本发明采用流模型以数据驱动的方式学习PK模型参数的后验分布，明显优于依赖预定义后验分布的传统方法，从而确保了PK参数的更为准确估计。

(2)增强的不确定性评估：本发明所估计的药代动力学参数分布中的标准差可作为一个可靠的不确定性指标，这有助于有效识别和排除不可靠的参数(均值)结果；进一步可以为医疗决策提供坚实的依据。

(3)本发明提供的方法及装置还可以为下游任务(例如，体内数据的胶质瘤分级)提供有价值的不确定性指标，从而提高其性能。如卓越的分类性能：在胶质瘤世界卫生组织(WHO)分级任务中，本方法能够有效地区分不同的分级，用于科学研究，也可以为临床医生提供了更为精确的诊断参考。

附图说明

图1为本发明所提出的框架和网络结构示意图。

图2为不同噪声水平下的平均不确定性与估计的平均绝对误差的关系。

图3为实施例所提供方法估计的参数图、相应的不确定性图和参数SNR图示例。

图4为实施例所提供方法估计的不确定性可用来提高胶质瘤分类性能。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

本发明提出的模型基于Pytorch框架实现。

本发明提供的基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法主要包括以下步骤：

S1：通过采集定量T1图像数据和DCE-MRI数据，计算得到每个被试者的血管内造影剂浓度变化曲线C

S1-1：将被试者置于3T磁共振成像系统中，以头部中心为扫描中心点，进行头部的图像采集。本实施例共采集了53个被试者的磁共振数据。

S1-2：在注射造影剂之前，设置flash3D序列以获取T1图像。使用了双翻转角double flip angle fast low angle shot(FLASH)技术，轴向切片数量为80。体素大小为0.9×0.9×1.5mm

S1-3：造影剂注射阶段，首先在48秒内按0.1mmol/kg体重给予碘酸钆二甲胺(gadopentetate dimeglumine；Gd-DTPA；Bayer，Berlin，Germany)。随后，用15mL盐水以2.0mL/s的速度进行冲洗。同时，使用了三维(3D)CAIPIRINHA-TWIST技术，其体素大小为0.9×0.9×1.5mm

S1-4：根据实验DCE-MRI数据计算出血管内造影剂浓度变化曲线C

进一步地，为了保证计算结果的准确性，所采集到的3D的T1定量图像需使用配准工具将其配准到DCE-MRI的第四帧上。同样地，对于DCE-MRI除第四帧以外的其他时间帧数据，也需要被配准到到第四帧上，以减少扫描过程被试头部微动对扫描结果的影响。

S2：生成仿真数据集：设定所要估计的药代动力学模型参数的范围，并在所设定的范围内进行均匀随机采样，生成药代动力学参数组；根据C

S2-1：设定药代动力学模型参数的范围为K

S2-2：从C

S2-3：使用eTtofts模型生成DCE-MRI数据，其中，eTofts模型的信号由公式

S2-4：向合成的信号中添加不同强度的莱斯噪声，使造影剂注射前DCE-MRI信号的SNR在21.8dB和36.98dB之间变化(平均SNR为27.80dB)。这参考了数据集中53位患者的全脑平均SNR。这可以使本发明提供的模型对一系列的SNR保持鲁棒性。

S2-5：合成500,000个DCE-MRI信号。从所合成的信号中，选取400,000个信号作为训练数据用于训练和验证，100,000个信号作为测试数据评估模型性能。从训练数据中抽取20％，即80,000个信号，设为验证集，以监控训练过程并防止过拟合。具体来说，使用46个被试的C

S3：构建基于流模型的参数分布估计模型(FPDEN)：FPDEN包括一个参数回归网络和流模型网络。

其中，参数回归网络估计参数分布均值和标准差，流模型网络用于学习参数的分布形状，并计算给定采样值在所预测分布上的概率值。

具体地，使用Python语言的Pytorch工具包，构造如图1所示的基于流模型的参数分布估计网络模型：图1中的(a)展示整体框架和网络运行过程的最大似然估计的计算方式；图1中的(b)描述了基于LSTM的回归网络结构；图1中的(c)展示了组成流模型的仿射耦合层的结构。

在本实施例中，FPDEN模型结构主要包括一个参数回归网络和一个流模型网络：如图1中的(b)所示的回归网络首先将两个通道的数据C

流模型网络结构使用了改编后的RealNVP网络结构，它由6个仿射耦合层构成，每个耦合层的结构如图1中的(c)所示。其中s

其中，图1中的Regression network为回归网络，training pathway为训练路径，Inference pathway为推理路径，Learned distribution为学习到的分布，Basedistribution为基本分布，Affi

具体地，流模型网络使用real-valued non-volume preserving(RealNVP)模型构建。这种模型由一系列简单映射组成，可以表示为：

其中，g

z

其中，⊙表示元素逐个乘积，s和t是尺度和平移函数。这种转换的雅可比行列式可以表示为：

由于其三角形状，其行列式可以有效地计算为：

S4：网络训练：，将

具体且：将步骤(2)生成训练集中的组织内造影剂浓度曲线和对应的血管内造影剂浓度曲线作为参数回归网络的输入，输出均值和标准差，该曲线对应的PK参数作为一个从参数真实分布上的采样值，使用输出的均值和标准差对该采样值进行归一化，然后输入到流模型中，计算其log似然值。模型的训练使用最大似然估计loss函数，模型的训练目标是使得给定的PK参数在预测分布上的log似然值最大。网络训练到收敛后，得到优化的模型参数。

回归网络和流模型是以端到端的方式同时训练的。在训练过程中，回归模型与标准化流模型均进行了80,000次迭代训练，采用Adam优化器，其学习率设为0.0005，动量值β1和β2分别为0.9和0.999。在所有网络训练中，批次大小均设定为128。

网络使用最大似然估计loss函数进行训练，并使用了重参数化技术，其表达式为：

其中，

由于参数回归网络的不确定性分支可能为负数。为了避免这种情况，定义回归网络的输出为

S5：PK参数估计：使用采集得到的T1定量图像和采集到的对应的DCE-MRI信号，，计算血管内造影剂浓度曲线。然后，逐个像素计算组织的造影剂浓度曲线。对于每个像素，将血管内造影剂浓度曲线和组织造影剂浓度曲线作为步骤(4)获得的回归模型的输入，输出对应的药代动力学参数的均值和标准差。

对基于流模型的参数分布估计模型FPDEN的评估：为了证明所提出的基于流模型的分布估计网络的有效性，将基于流模型的最大似然估计(MLE)损失方法

另外，本发明还对比了参数的分布假设对回归性能的影响，以证明本发明的准确性。这是通过将从流模型中学习到的分布替换为多元高斯分布来实现的。这种替换导致在模型训练中使用

表1展示了在测试数据集上，使用不同损失函数训练的回归模型结果的定量比较结果。从表1中可观察到，使用MLE损失

表1使用不同loss训练的参数回归模型的均方误差(MSE)

对基于流模型的参数分布估计模型FPDEN的不确定性验证：为了评估回归网络输出的不确定性，考虑每个参数的σ作为结果的不确定性程度指标。采用Monte-Carlo(MC)模拟法在模拟的世界卫生组织(WHO)IV级胶质瘤组织(K

本实施例模拟了WHO IV级胶质瘤组织的DCE-MRI信号，并向其中加入指定强度的莱斯噪声来模拟不同的SNR，并对每个SNR情况进行了100次噪声实现以计算MAE和平均不确定性。图2中的(a)-(c)的左列展示了不确定性与SNR之间的关系。右列中，点代表每种不确定性情况下的预测MAE(每个平均不确定性来自对应的SNR)。虚线表示使用线性回归拟合估计平均误差到输出的不确定性的结果，其相关系数(R^2)标记在图2中。

图2中的(a)-(c)揭示，随着SNR的降低，三个参数K

其中，图2中的Mean Uncertainty为平均不确定性，MAE(Mean Absolute Error)为平均绝对误差。

图3展示了含有活体数据的胶质瘤区域的参数图。参数SNR图通过将估计的参数视为信号，将估计的不确定性视为噪声水平来计算。V

其中，图3中的Estimated为估计出的参数值，Uncertainty为参数不确定性，SNR为参数的信噪比。

对基于流模型的参数分布估计模型FPDEN在活体数据上验证。使用步骤(1)所说明的磁共振扫描参数采集已被临床确诊为胶质瘤的53个被试者DCE-MRI数据，其中WHO II级有11例被定义为低级别胶质瘤，WHO III级11例和WHO IV级31例共42例被定义为高级别胶质瘤。然后由专业人员画出肿瘤区域作为ROI。计算ROI内PK参数的平均值作为胶质瘤分级的指标。

模型估计的不确定性可以作为使用DCE-MRI数据进行胶质瘤分级人物的参数筛选的依据，如图4中的(a)所示，展示了利用不确定性过滤(uncertainty filtering，UF)的数据处理流程。当使用UF时，ROI内只有部分被选取的像素PK参数被用来计算PK平均值。在此UF过程中，每个体素内的每个参数的不确定性都被独立评估：置信度较高的参数被保留，而置信度较低的参数则被过滤掉。为了表示相对不确定性，一个变异系数(COV，coefficientof variation)被定义并通过公式

其中，图4中的Parameter Map为参数图，Uncertainty Mask为不确定性掩码，Ucertainty Map为不确定性图，Masked Parameter Map为筛选后参数图，Sensitivity为敏感性，Specificity为特异性。

使用不确定性过滤(UF)可提升胶质瘤分级的分类性能。如图4(b)和(c)所示，分别展示了使用V

基于同样的发明构思，本发明实施例还提供了一种基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计装置，包括：

图像采集及处理模块，通过采集定量T1图像数据和DCE-MRI数据，计算血管内造影剂浓度变化曲线C

数据集构建模块，设定所要估计的药代动力学模型参数的范围，并在所设定的范围内进行均匀随机采样，生成药代动力学参数组；根据C

模型构建及训练模块，建基于流模型的参数分布估计模型FPDEN：包括参数回归网络和流模型网络，将

检测模块，采集定量T1图像数据和DCE-MRI数据，并计算血管内造影剂浓度曲线和对应的组织内造影剂浓度曲线，将其作为FPDEN中参数回归模型的输入，输出对应的药代动力学参数的均值和标准差。

基于同样的发明构思，本发明实施例还提供了一种计算设备，包括一个或多个处理器，存储器中存储有可执行代码，处理器执行可执行代码时，用于实现上述实施例中的基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法。以软件实现为例，作为一个逻辑意义上的装置，是通过其所在任意具备数据处理能力的设备的处理器将非易失性存储器中对应的计算机程序指令读取到内存中运行形成的从硬件层面而言，除处理器、内存、网络接口、以及非易失性存储器之外，实施例中装置所在的任意具备数据处理能力的设备通常根据该任意具备数据处理能力的设备的实际功能，还可以包括其它硬件，对此不再赘述。

基于同样的发明构思，本发明实施例还提供一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时，实现上述实施例中的基于流模型的动态增强对比磁共振成像中药代动力学参数分布估计方法：计算机可读存储介质可以是前述任一实施例所述的任意具备数据处理能力的设备的内部存储单元，例如硬盘或内存。所述计算机可读存储介质也可以是任意具备数据处理能力的设备，例如所述设备上配备的插接式硬盘、智能存储卡(Smart Media Card，SMC)、SD卡、闪存卡(Flash8 Card)等。进一步的，所述计算机可读存储介质还可以既包括任意具备数据处理能力的设备的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述任意具备数据处理能力的设备所需的其它程序和数据，还可以用于暂时地存储已经输出或者将要输出的数据。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。