胡德斌,张新明,孙浩,韩彦江,齐宏亮,唐大洋,路利军,吴湖炳,陈宏文

临床上正电子发射断层成像(positron emission tomography,PET)定量方法通常采用标准摄取值(standard uptake value,SUV)[1]。但SUV是一种半定量方法,当示踪剂在炎症病灶聚集时会导致假阳性,而在某些恶性程度不高的肿瘤会呈现假阴性[2-3]。针对这些问题,可以采用动态PET成像模式,通过动力学模型的应用,得到各组织器官的局部血流量、物质转运速率、代谢速度和受体结合率等参数图像,从而有效降低诊断的假阳性和假阴性[4-5]。在动力学模型的计算中18F-FDG PET的参数成像方法通常采用Patlak方法,其优点在于计算过程简单,计算速度快[6]。然而,针对18F-FDG PET Patlak参数成像模式需要扫描至少60 min,这会导致患者不能坚持或在扫描过程中带来较大的运动伪影,极大的限制了这项技术的应用。

当前,诸多方法相继被提出以实现18F-FDG PET短时间的Patlak参数成像,但短时间参数成像会引入高噪声和定量偏差的问题。为解决这些问题,其中一类方法是基于传统去噪方法,比如非局部均值去噪(non-local mean,NLM)方法[7],Wu等[8]利用该方法通过对不同时间段生成的参数图像进行去噪,发现使用20 min动态扫描的数据生成的参数图像去噪后能满足临床诊断效果。另一类方法是基于目前的深度学习方法,Huang等[9]通过采集200个60 min临床患者数据,每个患者数据生成两组数据,一组10 min的PET SUV图像,另一组60 min动态数据生成的Patlak参数图像,随后将PET SUV图像作为输入以及Patlak参数图像作为标签,训练卷积神经网络用于参数图像的生成。由于该方法是有监督的深度学习方法,需要首先获取大量的60 min动态扫描数据,对于大部分医院数据获取难度过高,限制在临床上应用。

目前无监督的深度学习方法应用越来越广泛,无监督方法最大的优势在于不需要大量的训练数据,利用本身的数据特点就可以实现数据处理。其中,深度图像先验网络(deep image prior,DIP)方法提出并应用于自然图像的去噪,能够取得良好的去噪效果[10]。为解决短时间Patlak参数成像图像噪声大问题,本文首先将DIP方法应用于Patlak参数图像去噪,为提升DIP方法的效果,本研究将原始DIP方法的网络输入由随机噪声图像改成PET活度图像,同时将深度学习网络采用编码器-解码器结构,使得处理后的短时间参数图像保持优异的定性和定量的效果,实现短时间参数成像,并与高斯去噪方法[11]和NLM方法进行比较。与此同时,本研究进一步比较不同短时间动态扫描协议下DIP方法去噪后参数图像的定量变化。

材料与方法

1.图像采集

仿真实验图像:本研究仿真1例Hoffman脑图像,包括灰质和白质,同时在脑图像中人为添加一个直径为10 mm的肿瘤。表1展示了用于仿真的动力学参数,具体参数值均来自文献[12]。图2a是大小为256×256的仿真脑PET图像,在仿真的动态扫描中扫描时间为60 min,然后在动态重建中总共分成92帧图像,具体分帧为2s/帧有30帧、5s/帧有12帧、10s/帧有6帧、30s/帧有4帧、60s/帧有25帧、120s/帧有15帧。根据表1中已知的动力学参数,我们采用Feng模型[13]可计算出不同组织的时间活度曲线(time activity curve,TAC),如图1b所示。

图1 仿真数据。a) 仿真Hoffman脑图像; b) 计算得到的不同组织时间活度曲线。

图2 动态扫描不同协议

表1 仿真的动力学参数设置

临床实验图像:本研究18F-FDG示踪剂用于临床采集经过审批(粤-ZHDY 2021001),临床数据是利用上海联影医疗生产的全身PET/CT uEXPLORER扫描的1例动态18F-FDG数据。患者注射6.43 mCi药物,数据从注射药物开始采集了60 min。动态重建中采用有序子集期望最大化方法(ordered subset expectation maximization,OSEM)重建,动态分帧分成92帧,具体分帧跟仿真数据一致。重建得到的PET图像大小为192×192×673,体素大小为3.125 mm×3.125 mm×2.89 mm。

2.实验方法

Patlak参数成像方法:Patlak方法是一种线性动力学参数估计方法,主要针对不可逆转的核素代谢,其表达式为:

(1)

式中,CT表示的是目标区域的活度值,CP表示动脉输入函数,是动脉血中示踪剂浓度,VE表示的是分布容积,tn表示帧动态PET活度图像帧测量的时间点,t*表示药物在体内分布达到动态平衡的时间点,Ki是组织摄取率常数,流动达到平衡后组织与血浆中药物浓度的比值,可用来作为定量评估指标,本文就是针对Ki参数图像进行处理。

一般来说,进行Patlak参数成像时需要扫描60 min PET动态数据(如图2所示Protocol 1),为了实现短时间参数成像,需要缩短动态扫描时间,目前有两种协议,如图2所示Protocol 2和3:①注射示踪剂先扫描一段时间T2,患者休息一段时间后再继续扫描另一段时间T3;②注射示踪剂一定时间后扫描一段时间T4。

考虑到公式(1)中的动脉输入函数需要从开始注射示踪剂到扫描时间段完整的信息,采用Protocol 1扫描时可完整的从降主动脉提取输入函数,Protocol 2和Protocol 3则需要利用模板输入函数[14]使用拟合方法[15]补齐缺失信息。Protocol 2相比Protocol 3多扫描了T2时间段,可以提供更多输入函数前期示踪剂在动脉中的个体化峰值信息。当我们利用3种协议T1、T3和T4时间段动态数据进行参数计算时,参数图像的噪声会随着扫描时间变短而增大。为了解决噪声增大问题,需要对参数图像进行去噪。本研究先利用Protocol 1分析T1时间变短情况下去噪方法的效果,再根据Protocol 2和Protocol 3分析不同短时间扫描方案的参数成像定量结果。

深度图像先验网络去噪方法:随着深度学习逐渐成为机器学习领域的研究热点, 深度学习为解决图像噪声问题提供了新的思路,具有更强大的去噪能力,更好的保持图像细节及纹理[16]。相比于有监督的深度学习网络,无监督的深度学习网络能够不依赖大量训练数据情况下达到优异的任务表现。DIP方法用于图像去噪,定义为:

x=f(|znoise)

(2)

(3)

其中x0表示噪声图像。

在初始DIP方法中,网络的输入是随机噪声图像,本文为了提升网络的去噪效果,提出条件深度图像先验方法,将参数图像计算对应时间段内多帧PET活度图像的平均活度图像视为先验信息,代替随机噪声图像作为网络输入,平均活度图像定义如下:

(4)

其中P(tn)表示在tn时间点的PET活度图像,j表示多帧PET活度图像开始的帧序号,k表示多帧PET活度图像结束的帧序号。此时整个参数图像去噪过程可以由如下的优化方程表示:

(5)

本文研究中深度学习网络采用编码器-解码器体系结构,总共网络分为6层,网络输入是PET平均活度图像,网络标签是噪声参数图像,网络输出是去噪后的参数图像。编码阶段网络降采样操作首先是的卷积层处理,随后以步长为2的降采样操作替代最大池化方法(max pooling),接着使用批量归一化方法(batch normalization,BN)[17]和泄漏整流线性单元(leaky rectified linear unit,LReLU)激活函数[18],针对降采样后的特征图像层的网络操作为的卷积层加上BN层和LReLU激活层处理,此时特征图像的图像维度保持不变。解码阶段为了去除棋盘伪影,本文提出将网络上采样操作采用双线性插值方法取代去卷积上采样方法,包含双线性升采样结构。上采样后的特征图像层的网络操作为首先进行的卷积层,BN层和LReLU激活层处理,随后又进行一次的卷积层,BN层和LReLU激活层处理。此外,我们使用跳跃连接以串联的方式连接编码器路径和解码器路径,以减少训练参数的数量并包含不同特征尺度的结构。

图3 提出的无监督深度学习去噪框架示意图

算法开发环境介绍:本文的Patlak参数成像算法使用PyCharm软件(版本为PC-222.3345.131)编程实现,深度学习网络采用pytorch深度学习框架,版本为pytorch 1.11.0,运行环境为一台DELL 2U机架式服务器,机型为PowerEdge R740,处理器为英特尔XE-ON十核银牌4210(2.2GHz),内存64GB,GPU型号为NVIDIA Geforce RTX 3090,显存为24GB,运行操作系统为Red Hat 4.8.5-44。

3.性能评价

针对去噪后的参数图像,采用偏离率(Bias)来定量计算跟真值的误差,变异系数(coefficient of variation,COV)计算图像的噪声。

偏离率表示目标值跟真值的偏差程度,定义为:

(6)

变异系数表示图像的噪声水平,定义如下:

(7)

4.参数优化

参数成像计算过程中分别使用不同时间段数据进行Patlak估计,具体时间段记为20～60 min,30～60 min,40～60 min,44～60 min。根据不同时间段数据进行Patlak计算得到的参数图像分别用3种不同的方法进行去噪,这时的噪声参数图像记为noisy image。图像数据进行高斯去噪,记为Gaussian;图像数据进行非局部均值去噪,记为NLM;方法2.2是本文提出的无监督的深度图像先验网络去噪,记为DIP。临床数据分析过程中因为20～60 min数据量充足,计算的Patlak参数图像信噪比高,同时参考Wu等[8]方法在数据处理中选取20～60 min的Patlak结果作为参考的真值图像,用于跟不同方法处理后结果进行比较分析。

临床实验中基于完整采集的60 min动态数据,我们模拟两种短时间动态采集协议(图2),第一种获取前5 min数据,再获取4个不同时间段数据,第二种直接获取4个不同时间段数据。考虑没有部分数据信息,动脉输入函数采用模板函数拟合方法[15]补齐缺失信息。两种协议生成的参数图像结果与采集60 min数据生成的输入函数生成的参数图像结果进行比较。

实验中由于不同时间段的参数图像噪声水平不一致,所以会调整不同的滤波参数。Gaussian方法的滤波参数为半高宽(full width at half maxima,FWHM),NLM方法的的搜索窗设为5×5,邻域窗设为3×3,通过设置控制滤波程度的参数来调控去噪效果,DIP方法中初始学习率设为0.01,迭代次数为控制噪声的参数。具体参数设置如表2所示。

表2 不同去噪算法的具体参数设置

结 果

图4是利用仿真数据使用不同去噪方法针对不同时间段参数图像去噪的结果。其中图4a一列图像是仿真的真值图像,图4b一列从上到下是使用20～60 min,30～60 min,40～60 min,44～60 min 4组时间段数据进行参数成像的结果,从图中可看到随着使用数据的时间段减短,参数图像的噪声增大。红色箭头所指的区域表示仿真的肿瘤区域,从图中可以看出40～60 min、44～60 min两组数据Gaussian和NLM方法不能很好识别肿瘤区域,而本文提出的DIP方法能够清晰的展示肿瘤细节。

图5给出了不同去噪方法在肿瘤区域计算的Bias和灰质、白质区域计算的COV变化曲线。曲线中的每个点表示不同的时间段数据计算的结果,从左到右为20～60 min,30～60 min,40～60 min,44～60 min。从图5a和图5b结果来看,经过不同去噪方法去噪后参数图像的COV值都会降低,其中NLM去噪效果好于Gaussian方法,而本文提出的DIP方法可达到最好的去噪效果。DIP方法去噪后的白质区域的COV值从44～60 min到20～60 min增加了11.71%,灰质区域增加了15.75%。而噪声图像(noisy image)的白质区域的COV值从44～60 min到20～60 min增加了265.50%,灰质区域增加了195.46%。与此同时,在不同时间段上Gaussian方法在去噪的同时相比于原始噪声参数图像会导致肿瘤区域的Bias值增加,NLM方法和DIP方法都会减少Bias。其中,本文提出的DIP方法Bias值最小,44～60 min和20～60 min的结果相比于真值的Bias值的分别为-6.39%和-7.41%,两个时间段Bias值的差值为1.02%,而噪声图像相应两个时间段的Bias值为-9.09%和-12.45%,相差3.36%。图5的结果可看出DIP方法能将44～60 min时间段计算的参数图像保持跟20～60 min在噪声水平和目标区域值偏差基本一致的效果。

图5 不同去噪方法在肿瘤区域和灰质及白质区域的Bias和COV变化曲线。a) 肿瘤区域跟灰质背景区域计算的Bias和COV结果; b) 肿瘤区域跟白质背景区域计算的Bias和COV结果。曲线中的每个点表示不同的时间段数据计算的结果,从左到右为20～60min,30～60min,40～60min,44～60min。 图6 不同时间段数据训练网络损失函数曲线。

图6展示了DIP方法中训练网络针对不同时间段临床数据的损失函数曲线,可看到随着网络迭代次数增加训练网络会逐渐收敛。其中20～60 min和30～60 min时间段数据在迭代1200次达到收敛,30～60 min和40～60 min时间段数据在迭代1000次和700次收敛,如箭头所示。

图7展示了44～60 min时间段参数图像使用DIP方法在700次、1200次、1600次网络迭代的结果,可观察到当网络达到收敛状态时,随着迭代次数增加,图像噪声会随之增大,存在过拟合现象,所以在选择网络训练的终止条件为如图6中损失函数达到最初收敛时刻的迭代次数。

图7 44～60min时间段DIP方法随迭代次数变化结果。a) 700次;b) 1200次;c) 1600次。

图8为临床数据不同方法去噪的结果。从上到下依次为20～60 min、30～60 min、40～60 min、44～60 min生成的参数图像。ROIA表示肿瘤区域,ROIB和ROIC表示背景区域,用于定量计算。图8展示了通过ROIA和ROIB及ROIC计算的Bias和COV结果。其中因为我们将Patlak方法在20～60 min时间段估计的参数图像作为参考的真值图像,所以在图9a和图9b中“noisy image”的第一个点Bias值是0。从定量结果来看,DIP方法相比于Gaussian和NLM方法在不同时间段都具有最小的Bias值,其中在20～60 min和30～60 min的Bias值基本一致,40～60 min和44～60 min的Bias值相比于20～60 min会增大5.04%和6.3%,而噪声图像44～60 min的Bias相比于20～60 min时增大15.54%。与此同时,通过ROIB和ROIC计算的COV值可看到DIP方法具有最好的去噪效果。相比于噪声图像在44～60 min时的COV结果,DIP方法可以ROIB和ROIC的COV值分别从0.44降至0.27和0.58降至0.26,下降率分别为162.96%和223.08%,可以有效的保证参数图像在44～60 min时的图像质量。

图8 不同去噪方法去噪的的临床参数图像。ROI A是代表肿瘤组织, ROI B是背景区域。a) 噪声图像; b) Gaussian方法; c) NLM方法; d) DIP方法。

图10展示了不同扫描协议下DIP方法去噪后图8所示肿瘤区域定量偏差,其中“Protocol 2-Protocol 1”为Protocol 2与Protocol 1之间的定量偏差,“Protocol 3-Protocol 1”为Protocol 3与Protocol 1之间的定量偏差。结果看随着数据时间段减少,两种协议定量偏差都逐渐增大。此外,Protocol 2的定量偏差小于Protocol 3,Protocol2的定量偏差范围为0.6%～2.88%,Protocol 3的定量偏差范围为1.34%～3.92%。

讨 论

相比于常规SUV定量方法,PET的参数图像能够提供组织摄取率和受体结合率等信息,可有效减少临床上的假阳性和假阴性诊断率。但是,18F-FDG PET参数成像过程需要至少扫描60 min动态数据,这大大限制了该方法在临床上的使用。针对Patlak参数成像方法,本文分析了使用20～60 min、30～60 min、40～60 min、44～60 min不同时间段数据生成的参数图像的噪声及定量变化。为解决噪声增大的问题,提出了深度图像先验网络用于参数图像去噪,结果显示DIP方法能够有效去除噪声,同时44～60 min数据生成的参数图像定量能达到20～60 min数据基本一致的效果。根据这一结果,我们分析了DIP方法在临床上两种短时间扫描协议应用效果,定量偏差<4%。因此我们可以设计具体短时间扫描协议,第一种协议从注射药物开始扫描5 min,随后休息39 min,然后再扫描16 min,整体时间可以将60 min扫描缩减为21 min。第二种协议采用基于模板的输入函数方法可以不用采集前5 min数据,直接从44～60 min时间段开始采集,这样可以将整体时间从60 min缩减为16 min。这两种扫描协议可以减少Patlak参数成像所需的采集时间,有利于参数成像方法在临床上应用。

深度图像先验网络是一种无监督深度学习方法,优点在于不需要大量的训练数据,针对本身数据就可以进行处理。效果来说设计合理的网络输入和深度学习网络可以达到优异的结果。但是,无监督深度学习方法存在过拟合的问题[19],在本文结果中可以观察到随着网络迭代次数的增加图像噪声增大。为解决这个问题,我们需要对迭代次数进行参数选择,通过提前结束迭代来防止过拟合。我们提出根据训练网络的损失函数收敛情况,达到最初收敛的迭代次数作为中止条件,所以在本文中可以看到不同时间段的数据的迭代次数是不一样的,分别为1200、1200、1000、700。

本研究进一步的工作需要采集多组临床数据进行诊断分析,通过实际的临床诊断结果来验证本文方法在临床的应用效果。综上所述,本文提出了一种基于深度图像先验网络的PET短时间参数成像方法。针对短时间参数图像噪声增加问题,本文提出的方法在仿真实验、临床数据实验中都展示了抑制噪声增加,保持图像细节和定量准确性的良好效果。