基于栈式自编码器的动态PET图像示踪剂动力学宏参数估计方法与流程

本发明属于pet成像技术领域，具体涉及一种基于栈式自编码器的动态pet图像示踪剂动力学宏参数估计方法。

背景技术

正电子发射断层扫描(positronemissiontomography，pet)是一种核医学功能成像技术，在生物医学研究和临床诊断中发挥着至关重要的作用。pet成像一般采用短寿命的放射性同位素(如¹⁸f、¹¹c等)来标记生物分子，这些生物分子在生物体内的代谢活动中衰变产生正电子，产生的正电子与体内的负电子发生湮灭反应产生一对方向相反能量为511kev的伽马光子对，伽马光子被环型的探测器所捕捉，从而得到发射数据，然后通过计算机分析构建体内示踪剂浓度的图像。

相比于静态pet成像，动态pet成像能够提供在不同的时间帧下的放射性示踪剂的活度分布图，反应出了生物组织在一段时间内的代谢过程，通过动力学模型建模，可以得到组织器官的局部血流量，代谢速率和结合能等功能参数，这些生理参数能够定量直观地表达生物组织的代谢情况，对生物体内组织器官是否异常具有重要的意义。在pet动态成像中，为了能够定量地分析生物体的生命活动状态，通常结合动力学模型来获得表征示踪剂在生命体内的代谢活动的量化参数，这些宏参数信息对于分析生物体内生理状态具有重要的意义，因此，如何快速准确地估计出这些参数就变得尤为重要。

目前，估计生理参数的方法主要分为两类：基于房室模型的参数估计和基于图模型的参数估计。前者依赖于房室模型的数量，且估计的结果稳定性不高；后者不依赖于房室模型的数量，但是需要事先确定示踪剂属于可逆的房室模型还是不可逆的房室模型，且该方法还依赖于拟合方法，拟合方法的好坏直接决定了生理参数的准确度。因此，在不依赖于房室模型的条件下，能够快速且准确地估计动力学宏参数是研究pet成像的一个非常关键的问题。

技术实现要素：

鉴于上述，本发明提供了一种基于栈式自编码器的动态pet图像示踪剂动力学宏参数估计方法，该方法不依赖于房室模型，且又能快速准确地获得动态pet图像。

一种基于栈式自编码器的动态pet图像示踪剂动力学宏参数估计方法，包括如下步骤：

(1)利用探测器对注入有放射性示踪剂的生物组织进行探测，采集得到对应探测器各晶体块的符合计数向量，进而构建动态pet的符合计数矩阵；

(2)根据pet成像方程以及符合计数矩阵，通过ml-em(极大似然期望最大化)算法重建出动态的pet浓度分布图像序列，进而从pet浓度分布图像序列中提取出每一像素点对应的tac(time-activitycurve，时间放射性活度曲线)；

(3)根据房室模型建立pet浓度分布图像的动力学模型，进而利用动力学模型计算获得每一tac对应的动力学宏参数；

(4)根据步骤(1)～(3)对不同的生物组织进行多次探测及重建计算，从而得到大量样本，每个样本即包含一条tac及其对应的动力学宏参数；

(5)将所有样本分为训练集、验证集和测试集，利用训练集样本对栈式自编码器进行训练及微调，得到动态pet图像宏参数估计模型，所述栈式自编码器由多个自编码器级联组成；

(6)将测试集样本中的tac逐一输入至动态pet图像宏参数估计模型中，从而估计出动态pet图像的动力学宏参数。

进一步地，所述步骤(2)中的pet成像方程表达式如下：

y＝dx+nosie

其中：d为系统矩阵，y为符合计数矩阵，x为pet浓度分布图像，nosie为包含反射和散射符合事件的测量噪声矩阵。

进一步地，所述步骤(2)中通过以下表达式从pet浓度分布图像序列中提取每一像素点对应的tac：

xi＝[x1i,x2i,x3i,...,xki]^t

其中：xi为第i个像素点对应的tac，xji为第j帧pet浓度分布图像中第i个像素点的浓度值，i和j均为自然数且1≤i≤m²，1≤j≤k，k为pet浓度分布图像序列的总帧数，m×m为pet浓度分布图像的尺寸大小，t表示转置。

进一步地，所述步骤(3)中计算tac对应动力学宏参数的具体过程为：首先利用动力学模型通过非线性最小二乘法计算出tac对应的动力学微参数[k1,k2,k3,k4]，进而根据动力学微参数[k1,k2,k3,k4]通过以下公式计算出tac对应的动力学宏参数[bp,dvt,ki]：

bp＝k3/k4

dvt＝k1(1+bp)/k2

ki＝k1k3/(k2+k3)

其中：bp为束缚势能，dvt为分布容积，ki为净示踪剂摄取常数。

进一步地，所述步骤(5)中的自编码器由输入层、隐藏层、输出层组成，栈式自编码器中前一个自编码器的隐藏层即作为后一个自编码器的输入层。

进一步地，所述自编码器的函数模型如下：

h＝g(wt+b)

z＝f(w'h+b')

其中：t、h和z分别为自编码器的输入层、隐藏层和输出层，w和b为输入层与隐藏层之间的权重参数和偏置参数，w'和b'为隐藏层与输出层之间的权重参数和偏置参数，g()为输入层与隐藏层之间的激活函数，f()为隐藏层与输出层之间的激活函数。

进一步地，所述步骤(6)中对栈式自编码器进行训练的具体过程如下：

首先，对于栈式自编码器中的第一个自编码器，将训练集样本中的tac作为该自编码器的输入层，以该自编码器输出层与输入层之差的二范数作为目标函数，通过反向传播和梯度下降法不断计算更新该自编码器中的权重参数和偏置参数，直至目标函数值小于设定的阈值，则该自编码器训练完成；

对于栈式自编码器中的其他自编码器，则使前一个自编码器的隐藏层作为当前自编码器的输入层，以当前自编码器输出层与输入层之差的二范数作为目标函数，通过反向传播和梯度下降法不断计算更新当前自编码器中的权重参数和偏置参数，直至目标函数值小于设定的阈值，则当前自编码器训练完成；

当所有自编码器训练完成后，将各自编码器的隐藏层堆叠一起形成栈式自编码器，进而利用训练集样本中的tac作为栈式自编码器的输入，使样本中的动力学宏参数作为栈式自编码器的输出真值，再次对栈式自编码器进行训练；

最后，利用验证集样本对训练完成的栈式自编码器进行验证及微调，得到最终的动态pet图像宏参数估计模型。

本发明首次将深度学习的思想引入动态pet示踪剂动力学宏参数估计中，其过程主要可以分为训练和估计两个阶段。在训练阶段，将动态pet示踪剂的浓度分布图像作为输入，动态pet示踪剂的宏参数作为标签训练自编码器，进而构建栈式自编码器；在估计阶段，将动态pet示踪剂的浓度分布图输入到训练好的栈式自编码器中即可估计出示踪剂的动力学宏参数。本发明从数据驱动的角度实现了在无房室模型的依赖下，快速且准确地估计动态pet图像的动力学宏参数。

附图说明

图1为本发明动力学宏参数估计方法的流程示意图。

图2为大脑phantom模板图。

图3为自编码器的构建过程示意图。

图4为栈式自编码器结构示意图。

图5(a)为示踪剂fdg的宏参数真值图像。

图5(b)为采用本发明方法估计得到的示踪剂fdg宏参数值图像。

图6(a)为示踪剂ace的宏参数真值图像。

图6(b)为采用本发明方法估计得到的示踪剂ace宏参数值图像。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

本发明基于栈式自编码器的动态pet图像示踪剂动力学宏参数估计方法，整体框架如图1所示，具体包括如下步骤：

s1.将示踪剂注入到生物组织中，利用探测器对注入有放射性药剂的生物组织进行探测，采集得到对应探测器各晶体块的符合计数向量，进而构建动态pet的符合计数矩阵y。

s2.根据pet成像原理，通过ml-em算法求解以下方程得到示踪剂的pet浓度分布图像x：

y＝gx+noise

其中：g为系统矩阵，noise为测量噪声矩阵(包括反射符合事件和散射符合事件)。

然后将不同帧的pet图像x按照像素点顺序排列成如下形式，作为训练集的输入数据。

x＝[x1,x2,···,xm×m]

其中：xi＝[x1i,x2i,···,xki]^t，k为动态pet图像总共的帧数，xki为第k帧pet图像的第i个像素值，m×m为图像的分辨率即图像的像素总数。

s3.根据房室模型建立动力学模型，使用非线性二乘法获得微参数[k1,k2,k3,k4]，进一步根据以下公式获得宏参数[bp,dvt,ki]，并将其作为训练集的标签数据。

其中：bp为束缚势能，dvt为分布容积，ki为净示踪剂摄取常数；示踪剂浓度分布图像中的每一个像素点都有对应的宏参数，因此宏参数可表示为o＝[o1,o2,···,om×m]，其中表示第i个像素点处的三个宏参数值。

s4.利用上述训练集训练栈式自编码器。栈式自编码器则是由多个自编码器级联组成，自编码器由输入层、隐藏层和输出层所组成，其中从输入层到隐藏层部分称为编码器，隐藏层到输出层部分称为译码器；栈式自编码器中前一个自编码器的隐藏层的输出作为后一个自编码器的输入。

首先设定学习速率α，每个自编码器的隐藏层节点数n,l,t，栈式自编码器的层数s，初始化每个自编码器的权重参数w,w'和偏置参数b，设定训练阈值，自编码器的函数模型如下：

h＝f(wx+b)

o＝f(w'h+b')

其中：h为隐藏层的输出，o为输出层的输出，f为激活函数，w和b分别为编码器的权重和偏置参数，w'和b'分别为译码器的权重和偏置参数。

s5.构建多个自编码，将x＝[x1,x2,···,xm×m]作为第一个自编码器的输入数据，o＝[o1,o2,···,om×m]作为输出层的真值，构建过程如图3所示，利用反向传播算法和梯度下降算法不断更新权重参数w,w'和偏置参数b。

对于栈式自编码器中的第一个自编码器，将训练集的输入数据作为该自编码器的输入，以自编码器的输出层的输出与输入量之差的二范数作为目标函数，通过反向传播算法和梯度下降法不断地计算和更新该自编码器中的权重和偏置参数，直到目标函数的值小于设定的阈值时，自编码器训练完成。

对于栈式自编码器中其他自编码器，将前一个自编码器的隐藏层作为下一个自编码器的输入层，以该自编码器的输入量和输出量之差的二范数作为目标函数，通过反向传播算法和梯度下降法不断地计算和更新该自编码器中的权重和偏置参数；其中目标函数的表达式如下：

其中：表示第i个像素点的宏参数真值，表示第i个像素点的预测值，θ表示栈式自编码器中的权重。

s6.将训练好的自编码器级联在一起就得到了栈式自编码器，如图4所示，然后将x＝[x1,x2,···,xm×m]和o＝[o1,o2,···,om×m]作为栈式自编码器的输入和真值，再次进行训练，微调权重参数w,w'和偏置参数b。

s7.给定新的示踪剂的发射数据根据ml-em算法求解出pet浓度分布将其排列成

s8.将的每一列都输入到训练好的栈式自编码器中，输出得到估计值即获得了该示踪剂的pet图像动力学宏参数。

以下我们采用大脑phantom模拟数据验证本发明的准确性，采用的大脑phantom如图2所示，设定的放射性核素及药物为¹⁸f-fdg和¹¹c-acetate。我们使用蒙特卡洛模拟了3种不同的采集时间的动态pet图像，每种采集时间又分为了3种不同的采集间隔，对于每种采集间隔，我们模拟了3种不同的计数率，分别为1×10⁵、1×10⁶和5×10⁶，详细数据如表1所示：

表1

我们模拟的pet图像分辨率为64×64，总帧数为28帧，主要评价指标包括偏差bias，方差variance：

其中：ui，分别表示第i个像素点真值、第i个点像素的估计值和平均的像素估计值，表2和表3分别展示了fdg和ace两种示踪剂的动力学宏参数估计结果。

表2

fdg

表3

ace

从表2和表3可以看出，本发明方法得到的结果与真值之间的差异很小，无论是哪种示踪剂，其偏差基本都小于0.5，方差基本都小于0.05，验证了本发明方法的有效性。图5(a)～图5(b)以及图6(a)～图6(b)直观地展示了两种示踪剂的宏参数图，图5(a)～图5(b)展示了示踪剂为fdg，采样时间为50分钟，采样间隔为5×24s的估计结果，其中图5(a)表示宏参数真值，图5(b)表示本发明方法获得的宏参数值。图6(a)～图6(b)展示了示踪剂为ace，采样时间为50分钟，采样间隔为5×24s估计结果，其中图6(a)表示宏参数真值，图6(b)表示本发明方法获得的宏参数值，从两幅图中可以看出本发明能够准确地估计动态pet图像的动力学宏参数。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。