基于分步子空间映射的无标记膈肌运动预测算法

余 航 李晨阳 余绍德 冯冬竹 许录平

呼吸运动(Respiratory motion)是放射治疗(Radiation therapy,RT)中,导致误差和不确定性的重要因素[1−2].一方面,现代放射治疗技术可以达到毫米级的传输精度,如调强放疗(Intensity-modulated radiation therapy,IMRT)[3];另一方面,呼吸运动可以导致胸部和腹部的肿瘤移动多达35 毫米[4].在此情况下[5],如果不对呼吸运动进行有效控制和补偿,要么肿瘤得不到足够的放射计量,要么肿瘤周围的健康组织将会受到放射的伤害[6−7].因此,有效减少呼吸对器官和肿瘤运动的影响,对整个放射治疗的进程和效果显得尤为重要[8−10].

针对这一问题,研究者提出了多种不同的解决方法,包括:1)屏气法[11−13];2)浅呼吸法[11−19];3)呼吸门控法[15−16];4)动态目标跟踪法[17−21].呼吸门控法和动态目标跟踪法都需要实时了解肿瘤的位置.呼吸门控法通过在呼吸周期的预定阶段打开放射光束,从而减少正常组织的照射量.动态目标跟踪法可以让光束和肿瘤进行同步的运动,实现连续照射,是目前最为有效的呼吸管理技术.肿瘤的位置可以通过直接方式和间接方式来获得.直接方式通过放射线成像系统对肿瘤或植入体内的标志物进行成像,该方法可以获得人体内的结构信息,但会给患者带来额外的辐射负担[10,22−27],并且需要进行昂贵的、创伤性手术,因此目前还没有广泛应用.

间接方式是使用呼吸信息对内部解剖结构的运动进行建模和预测.Hoisak 等[27]通过肺活量测定和实时位置跟踪系统,研究肺部肿瘤运动与呼吸运动的相关性.该方法的局限是,只使用了线性相关系数来度量肿瘤和参考物在一个方向上的运动相关情况.Paganelli 等[8]从核磁共振图像(Magnetic resonance imaging,MRI)中提取SIFT (Scale-invariant feature transform)图像特征,对基于图像特征的跟踪方法与基于标志物的跟踪方法进行对比.该方法采用二次多项式进行建模,需要在人体表面放置标志物,并计算45个图像特征的均值,较为繁琐.文献[9]通过实验表明,相比于其它解剖结构,膈肌与呼吸运动具有更强的相关性.根据这一结论,很多研究者采用膈肌作为参考点,对肿瘤运动情况进行研究.文献[28]针对不同患者,提出了一种基于膈肌的患者特异性呼吸模型,取得了良好的实验结果.

另一方面,人体内部结构的运动情况,可以通过机体外部的观察来进行预测.最近的研究[9]表明,人体整个外表面的运动与内部结构的运动具有很强的相关性.据此,文献[29−30]使用K-means 算法对人体胸腹部表面的运动轨迹进行聚类,以产生呼吸运动的参考点.但是由于K-means 算法易受到初始值和数据结构的影响,该方法容易受到噪声的干扰.Malinowski 等[1]将三个标记物固定在患者贴身背心上,据此来研究人体表面与肿瘤位置的相互运动关系.Ozhasoglu 等[4]通过使用光学监测装置观察附着在胸部和腹部的红外发光二极管(LED),来研究呼吸运动对肺和胰腺肿瘤运动的影响.这两种方法[1,4]均需在人体外表面放置标记物,因此较为繁琐,易受到操作者的影响.

对体内膈肌和体外胸腹表面之间关系的准确、系统的研究,不仅可以提高呼吸运动管理的准确性,而且可以显著降低放射治疗的复杂性,但相关工作仍然较少,且要么需要患者佩戴标记物,较为繁琐,要么模型建立比较简单.考虑到这两个器官具有不同的运动模式,其数据具有不同的分布结构,形成不同的子空间,无法直接进行准确的映射.针对这一问题,本文提出了一种新的分步子空间映射(Twostep subspace mapping,TSSM)算法,通过体外无标记胸腹部表面的测量来预测体内膈肌的运动情况.该方法无需患者佩戴标记物,且采用多种线性、非线性模型,系统地对体内膈肌和体外胸腹表面运动情况进行对比研究.本文首先针对4D CT 图像,利用三维图像分割方法,对体内膈肌和体外胸腹部表面的位移进行准确的测量.为了解决域自适应(Domain adaptation)问题,TSSM 首先利用主成分分析法(Principal component analysis,PCA)为每组数据构建特征子空间.然后通过线性岭回归(Ridge regression)优化过程,将膈肌数据与胸腹部表面数据连接起来,得到一个子空间映射模型.基于得到的映射模型,该方法可以通过体外胸腹部表面数据,对体内膈肌的运动情况进行预测.为了研究数据的非线性相关性,本文进一步将TSSM 扩展到基于核的分步子空间映射算法(Kernel TSSM,kTSSM),包括多项式核和高斯核.本文给出了两种算法的解析解,其不需要递归迭代的操作,具有运算速度快的特点.对比实验表明,这种分步映射的策略可以有效解决跨空间数据的预测问题,具有很好的准确性和鲁棒性,优于经典的线性模型和ANN模型,本文也进一步给出了内部膈肌与外部胸腹部表面之间的运动关系模型.

本文的其余部分安排如下:第1 节具体介绍了所提方法和数据,实验在第2 节中介绍,第3 节给出结论和讨论.

1 方法

1.1 内部膈肌位移提取

当获得患者的4D CT 图像后,本文设计了一种三维图像分割方法对其进行分割.在4D CT 图像中,由于具有不同的组织结构,身体区域的体素值(Voxel)均大于800,而背景区域和肺部区域主要为空气,其体素值在50 以下.考虑到图像的背景、身体区域和肺部区域具有显著的灰度差别,本文采用阈值法对其进行初步分割.具体地,本文设置阈值(本文设为500),首先将身体区域从背景中分割出来;然后根据区域面积,将图像中的最大目标区域作为身体区域;最后进一步通过形态学方法,去除背景中的孤立区域并填充身体区域中的孔,从而得到身体区域的分割结果.在获得身体区域之后,由于肺部区域在身体区域内部,并与身体区域的灰度值明显不同,因此可以进一步执行阈值法分割出肺部区域,此时分割算法会将肺部区域作为目标,将身体区域作为背景.在最终的三维分割结果中,本文将背景体素设置为0,身体区域体素设置为1,肺部区域体素设置为3.图1 显示了第55、65 和75横断平面(Axial plane)的分割结果.图2(a)显示4D CT 第1 相位的三维分割结果,我们可以进一步地将肺部区域和身体区域分开,如图2(b)、图2(c)所示.

图1 三维图像分割结果横断平面展示图,背景体素值为0,身体区域的体素值为1,肺部体素值为3Fig.1 3D segmentation images shown on the axial plane,where the background voxels is set to 0,the voxels of body area to 1,and the voxels of lungs to 3

图2 三维图像分割结果,及其对应的肺部区域和身体区域Fig.2 A 3D segmentation result and the corresponding separated 3D lungs and body

膈肌位于肺部下面,并且与肺部的下边界同步运动.考虑到难以直接通过4D CT 观测到膈肌结构,因此本文使用肺部的下边界来代表膈肌.本文工作面向的是肝癌立体定向体部放疗(Stereotactic body radiation therapy,SBRT),其主要受到右肺的影响.因此,本文设计了右肺自动确定方法,来提取右肺区域,具体操作如下:首先,在肺部区域的三维分割结果上(如图2(b)),从左到右分别计算各个矢状平面(Sagittal axis)上肺部区域的面积(体素标签为3 的数目),这样以左右方向为横轴,以肺部面积为纵轴,便可得到沿左右方向的肺部面积曲线,如图3 所示.可以注意到,由于肺部的对称性,在两肺之间肺部面积值比较小,因此肺部面积曲线像字母 “M”一样具有两个峰,每个峰对应一个肺.因此,选择肺部面积曲线两个峰之间的最低点,便可以确定右肺所包含的切片.

图3 肺部面积沿左右方向的曲线Fig.3 The curve of lung area along the left/right direction

接下来,在右肺的三维分割图像中,可以通过从下到上(解剖坐标系下)确定肺部区域的最低体素(在分割图像中,身体体素值为1,肺部区域体素值等于3),从而获得右肺下表面的三维图像.在肺部下表面,有一些高于膈肌的体素,这些体素可以通过形态学运算进行去除.最后,可以获得右肺下表面膈肌的三维空间坐标,表示为:

1.2 外部胸腹部表面位移提取

本节从三维分割结果的身体区域(如图2(c)所示)中获取胸腹部表面,进而计算胸腹部表面的位移.首先,沿前/后方向找出身体的边界体素(体素值为1),获得胸腹部表面.图4(a)显示了胸腹部表面在前/后方向的透视图,其中像素灰度代表距最前端冠状平面的距离.我们进一步对初始胸腹部表面进行处理,首先在得到的初始胸腹部表面图上(如图4(a)),沿左/右方向找到最宽的位置,该位置对应人体的肩膀,我们去除肩膀之上的区域.考虑到人体是一个类似的圆柱体,其在冠状面的投影,将会包含正面的胸腹表面和人体两侧的区域,因此,为了更加准确地获得胸腹部表面,本文运用腐蚀操作,来进一步地去掉人体两侧的区域.当人体进行呼吸时,人体两侧区域会沿左/右方向伸缩,其在冠状面的投影有显著的变化,且位于初始胸腹部表面的边缘区域.因此,在对不同相位的初始胸腹部表面进行腐蚀后,得到的人体胸腹表面大小不一致,考虑到是同一个病人,为了统一大小,本文对所有相位的胸腹表面取交集,我们可以得到如图4(b)所示的患者胸腹部掩膜.将该掩膜与初始胸腹部表面取交集,我们可以进一步获得最终胸腹部表面,如图4(c)所示,其表示为:

图4 胸腹部表面Fig.4 The thoracoabdominal surface

相比于胸腹部的面积,胸腹部位移是非常小的,因此本文采用ICP 算法(Iterative closest point)[31]来计算胸腹部表面的位移.ICP 算法的输出包括一个三维旋转矩阵和一个平移向量.通过测试发现,人体胸腹部运动的旋转角度小于0.5 度,其位移主要来源于平移,因此可以根据ICP 算法输出的平移向量,来计算胸腹部表面的位移,表示为:

1.3 分步子空间映射

在提取胸腹表面的位移数据x ∈X和膈肌的位移数据y ∈Y之后,本节将要在它们之间建立映射F:X →Y,这样对于任意新输入测试数据,可以通过观测胸腹表面的位移数据xTest,对膈肌的位移数据进行预测yTest=F(xTest).然而,胸腹表面和膈肌是两种不同的解剖器官,其运动模式也不相同,所产生的两种数据具有不同的分布结构.例如,在实验中发现,有的病人膈肌位移很大,而其胸腹表面的运动却很小.因此难以在这两种数据集之间直接构造映射.为了解决这一问题,本节提出一种新的分布子空间映射算法(Two-step subspace mapping,TSSM).

本文将数据分为训练数据和测试数据,x=xTrain∪xTest,y=yTrain∪yTest.假设器官结构相似的人,其胸腹部表面与膈肌有相似的运动情况.基于此假设,本文定义两种距离:1)子空间域内距离:dintra(xi,xj)或dintra(yi,yj),i,j=1,2,···,N,计算在同一个子空间的两个数据之间的距离;2)子空间域间距离dinter(xi,yj),i,j=1,2,···,N,计算来自不同子空间的两个数据之间的距离.TSSM 算法的核心是找到一种映射,使同一患者数据的子空间域间距离dinter(xi,yi),i=1,2,···,N小于不同患者的子空间域间距离dinter(xi,yj),ij.子空间域内距离dintra能够准确度量两个器官运动的相异性.这样对于新输入的数据,我们通过对观测数据xTest作映射来近似表示yTest.下面给出子空间的构造方法和映射方法.

膈肌位移数据和胸腹部表面位移数据的维数均为K维,因此可以生成K维子空间,选择一组数据作为子空间的基,其他数据由这些基来表示,但这样会导致数据冗余或不完整.为了克服这个问题,本文采用PCA 来构建子空间[32−34].PCA 通过最大均值差异最小(Maximum mean discrepancy,MMD)的准则,找到包含数据本质特征的公共子空间.

首先,对膈肌的位移数据和胸腹部表面的位移数据进行归一化,表示为K维向量的形式,如下所示:

其中,λ ≥0 是控制收缩量的规则化参数.对式(6)关于β求导,并使等式等于0,则可以获得最优表达式,如下所示:

图5 TSSM 算法的流程图Fig.5 The flowchart of the propose TSSM

与支持向量机(Support vector machine,SVM)相似,仅使用线性模型进行子空间映射,无法充分地描述数据之间的关系.为了解决这个问题,本文进一步将TSSM 推广到基于核的分步子空间映射算法(kTSSM).令ϕ为非线性空间 Rx到高维空间 Zx的映射,核函数为K(xi,xj)=〈ϕ(xi),ϕ(xj)〉.可以注意到,预测数据可通过变换后数据ϕ() 的线性组合来表示.因此,基于核函数的非线性岭回归优化变换可以表示为:

现有关于体内组织器官与体外基准点相对运动的研究[1−2,4,7−10],主要集中在线性模型或二次指数模型,为了更加详细、全面地研究体内膈肌和体外胸腹表面的相对运动关系,以进行准确的膈肌位移预测,本文除了构造分步线性映射模型,并进一步将其扩展到了基于非线性核的映射模型.本文采用了两种最常见的核函数:多项式核和高斯核,其中多项式核可以实现将低维的输入空间映射到高维的特征空间,其并不只限于二次函数,通过升维操作可以提高建模的准确性,但多项式核的缺点是参数较多,容易出现过拟合现象.高斯核函数是最常用的核函数之一,可以将数据映射到无限维空间.考虑到本文医用数据样本较少,维数较高,且没有数据分布的先验知识,本文进一步选择高斯核对体内膈肌和体外胸腹表面的相对运动关系进行度量.在实验中,本文将对线性模型和两种非线性模型进行对比研究,通过设置不同的参数,给出膈肌与胸腹部表面之间运动关系的准确度量.

2 实验

2.1 患者情况与数据采集

本研究中使用的临床数据由Philips Brilliance CT Big Bore Scanner 设备获得,在这项研究中采用了来自20 位患者的4D CT 数据集,每个患者包括10个相位(具体为平均呼吸周期的0%,10%,···,90%).每个胸部CT 的大小为512×512×185个体素(Voxel),尺寸分别为1.1719 mm×1.1719 mm×3 mm (分别对应于A/P、I/S和L/R三个方向).在采集数据时,病人按正常状态、有规律的呼吸.

在实验中,我们选择相关研究中的经典线性模型和非线性模型,与本文所提三种模型进行对比,具体如下:OLS (Ordinary least squares)多变量线性回归模型[1]、人工神经网络模型(Artificial neural network,ANN)、TSSM、多项式核的kTSSM 和高斯核的kTSSM.本文采用三个指标,对5 种算法的性能进行定量评估:均方误差(Mean-square error,MSE)、R2误差和平均绝对百分比误差(Mean absolute percentage error,MAPE),计算方法如下:

其中,yi是模型的实际值,是模型实际值的均值,是模型的预测值.MSE 衡量模型的预测均方误差,MSE 的值越低,预测结果越好.R2比较模型的优劣性,消除数据分布域对实验结果的影响.R2的范围为(−∞,1],值越大预测效果越好.MAPE 衡量模型的预测值与实际值的平均偏差,消除数据取值范围的影响,值越低预测结果越好.

在实验中,本文采用K-fold 策略将整个数据集分为训练数据和测试数据,其中80%的随机样本作为训练数据,其余样本进行模型测试.为了评估该算法的鲁棒性,在训练数据和测试数据中分别添加高斯噪声σ={0.1,0.2,0.3,0.4,0.5},从而分析不同级别的高斯噪声对算法的影响.本文算法参数选择如下:对于岭回归优化,λ={10−10,10−9,···,1010};在高斯核K(x,y)=中,σkernel的值等于{0.1,0.2,0.4,0.8,1.6,3.2,6.4,12.8};在多项式核K(x,y)=(xTy+c)d中,d={0.5,1,···,11},c={−10,−9,···,10};对于ANN,本文采用两层前馈神经网络,其隐层节点数依次选择{2,3,4,···,29,30}.对于每个算法的每组参数,本文单独运行算法100 次,选择每个算法的最优参数组合,并获得统计结果.所有算法的运行环境相同,均为Inter Xeon 3.6 GHz,32 GB 内存,Windows10 和MATLAB R2019b 平台.

2.2 膈肌和胸腹部表面的位移

图6 膈肌和胸腹部表面的位移,每一种颜色的曲线对应一个患者的数据Fig.6 The displacement of diaphragm and thoracoabdominal surface,where each color corresponds to a specific patient＇s data

由于性别、年龄、体型和情绪等都会对身体器官的运动产生影响,不同患者器官运动的幅度不尽相同.因此,根据式(5)对膈肌和胸腹部表面的位移进行归一化,结果如图7 所示.从图7 可以看出,不同患者胸腹部表面在前/后方向,膈肌在下/上运动方向的数据分布更加一致.

图7 归一化后的膈肌和胸腹部表面的位移,每一种颜色的曲线对应一个患者的数据Fig.7 The displacement of diaphragm and thoracoabdominal surface after normalization,where each color corresponds to a specific patient＇s data

2.3 三个方向的预测

在本节中,将根据胸腹部表面在三个方向的运动,来预测膈肌在三个方向上的运动.本文将同一患者三个方向10个相位的数据连接在一起,形成一个包含30个元素的向量,胸腹部表面位移向量为输入,膈肌位移向量为输出,预测结果如表1 所示.其中每个单元有两个数据,前一个值代表100 次独立运行的平均值,括号中的值代表100 次独立运行的标准偏差.本实验记录MSE 和MAPE 的最小值,记录R2的最大值,加粗标记代表5 种算法中的最佳结果.在表1 中,从上到下,噪声等级依次增加.从表1 可以看出,同为线性模型,TSSM 算法优于OLS 模型;非线性模型ANN 的表现较OLS 模型要更好.所提的三个算法的预测结果相近,MSE值和R2值均令人满意,但多项式核的kTSSM 算法获得多数的最佳预测值.值得一提的是,由于向量取值的不平衡性,两种组织器官在有的相位和方向上移动范围大,在其它相位和方向上移动范围小,取值范围小的预测值会增加MAPE 值,因此使得MAPE 值均超过100%.

表1 在三个方向上的预测结果Table 1 Prediction results in three directions

图8 显示5 种算法预测结果,其对应100 次独立运行中,MSE 的中间值对应的结果.图8(a)为 OLS模型,图8(b)为ANN,图8(c)为TSSM,图8(d)为多项式内核的kTSSM,图8(e) 为高斯内核的kTSSM.在每个子图中,有4个预测结果(占总数的20%),其中横轴表示10个相位,纵轴表示实际位移值,红色曲线表示实际值,蓝色曲线表示预测值.在横轴上,L/R表示左/右方向,I/S表示下/上方向,而A/P表示前/后方向.从图8(a)可以看出,OLS 模型在下/上方向上可以预测膈肌运动,但是在某些相位,OLS 模型在左/右和前/后两个方向上的预测效果不佳,仍然存在误差.同为线性模型,TSSM 算法具有更好的预测效果.与OLS 模型相比,非线性模型ANN 的表现更好,而本文设计的kTSSM (多项式)在大多数情况下,取得了最佳的预测结果.

图8 三个方向的预测结果.该结果对应于100 次运行MSE 的中位数值Fig.8 Prediction results corresponding the median MSE value in three directions after 100 independent runs

图9 显示了100 次独立运行的统计箱图,其中横轴表示三个方向10个相位,纵轴表示预测值和实际值之间的绝对误差.图9 子图中的每个方框都包含100 次独立运行的结果,在每个方框上,中心的标记表示中值,而框的底部和顶部边缘分别表示第25个和第75个百分点,“+”符号代表异常值.从图9我们可以清楚地看到,图9(a)中纵轴的范围比其他子图中的大得多,图9(b)中纵轴取值范围居中,而图9(c)～图9(e)的误差绝对值比较小,这也说明本文所提的方法具有更优、更稳定的预测性能.从图9(c)～图9(e)可以看到,下/上方向上的统计值低于其他方向上的统计值;对于下/上方向,位于中间相位(5、6、7 相位)的预测误差较小,这表明当患者处于满吸气状态时,预测效果会更好.

图9 三个方向各相位预测误差的统计箱图Fig.9 Boxplot of the prediction performance on every phase in three directions

2.4 主方向的预测

在本节中,仅使用胸腹部表面的前/后方向,对膈肌的下/上方向的运动进行预测.表2 显示了100 次独立运行后的预测结果.图10 显示5 种方法对应MSE 中间值的预测结果.图11 为每个相位的预测误差绝对值的统计箱图.从表2 可以看出,TSSM 算法的性能优于OLS 模型.ANN 算法结果优于OLS 模型.在没有左/右和前/后方向的数据影响的情况下,TSSM 算法的性能得到了很大的提高,而且MAPE 值已降至100%以下,证明了TSSM算法的良好性能.从图10 可以看出,本文所提方法均取得了令人满意的预测结果.图11 可以看出,在第5、6、7 相位上的预测结果,要好于在其它相位上的预测结果,这也进一步证明,当患者处于满吸气状态时,从外部体表能更加准确地预测体内结构的运动.

图10 主方向的预测结果,该结果对应于100 次运行MSE 的中位数值Fig.10 Prediction results corresponding the median MSE values in the principal direction after 100 independent runs

图11 主方向各相位预测误差的统计箱图Fig.11 Boxplot of the prediction performance on every phase in the principal direction

表2 在主方向上的预测结果Table 2 Prediction results in principal direction

2.5 噪声的影响

本节验证算法的鲁棒性.噪声对三个方向预测性能的影响如图12 所示.在图12 的子图中,横轴表示高斯噪声的等级,纵轴表示算法100 次独立运行后指标的平均值.图13 是噪声对主方向预测的影响,可以看出,随着噪声等级的增加5 种算法的性能都会下降.但是与OLS 模型相比,ANN 算法和TSSM 算法的下降速度较慢且稳定,这表明ANN,TSSM 和kTSSM 算法比OLS 模型对噪声的鲁棒性更高;对于准确性和稳定性,尤其是在MAPE度量方面,具有多项式核的kTSSM 算法在5 种算法中表现最好.

图12 噪声对三个方向预测性能的影响Fig.12 Influence of noise on the prediction in three directions

图13 噪声对主方向预测性能的影响Fig.13 Influence of noise on the prediction in principal direction

2.6 关于算法参数的讨论

TSSM 和kTSSM 取得最优值时,参数设置如表3 所示.对于线性模型TSSM,其参数只有岭回归一个参数λ,其取值范围为λ={10−10,10−9,···,1010}.在绝大多数情况下,λ的最优设置为1,只有两次为0.1,这两次对应λ=1 的预测结果与最优结果接近.因此,TSSM 对膈肌位移的预测是比较稳定的,受参数影响较小,可直接设置为1 即可.对于kTSSM(高斯)模型,其参数为λ和σkernel,其中λ的取值范围与TSSM 相同,σkernel={0.1,0.2,0.4,0.8,1.6,3.2,6.4,12.8}. 可以看到,在绝大多数情况下λ=0.1,σkernel=0.1 或0.2,只相差一个取值间隔.在实验结果中我们发现,当σkernel=0.1 或σkernel=0.2 时,kTSSM (高斯)模型的预测结果非常相近,具有很好的稳定性.kTSSM (高斯)模型的参数设置,可按如下策略:当获得的数据没有噪声或者比较干净时,可直接设置为λ=0.1,σkernel=0.2;当数据噪声比较 大 时,参 数 可 直 接 设 置 为λ=0.1,σkernel=0.1.对于kTSSM (多项式)模型,其具有三个参数λ、d和c,其中λ的取值范围与TSSM 相同,d={0.5,1,···,11},c={−10,−9,···,10}.可以看到,相比于TSSM 和kTSSM (高斯)模型,kTSSM (多项式)模型的参数变化浮动较大,特别是同时进行膈肌三个位移方向的预测时,岭回归参数λ、指数d随着噪声的增加逐渐降低,最后维持在λ=10 和d=1.这主要是因为在没有噪声或噪声较小时,kTSSM(多项式)模型发生了过拟合现象;当噪声增加时,降低了kTSSM (多项式)模型对数据的敏感性,因而最优参数可以达到相对稳定.kTSSM (多项式)模型对于主方向的预测,受参数影响较小,比较稳定,λ在1 上下浮动,d在0.5 和1 之间选择,参数c设置为4 即可.从kTSSM (多项式)模型的参数d可以看出,胸腹部表面与膈肌的相对位移,近似满足线性模型(d=1),这一点也可以通过TSSM 的参数稳定性得到验证.尽管kTSSM (多项式)模型在实验中取得了最好的预测结果,但是考虑到kTSSM(多项式)模型具有更多的参数需要设置,并容易出现过拟合现象,本文认为在实际中,只需要采用线性模型TSSM 即可达到准确的预测结果.

表3 TSSM 和kTSSM 取得最优预测结果时的参数设置Table 3 Parameter setting of TSSM and kTSSM for optimal results

本文进一步使用TSSM 模型研究PCA 算法中不同的阈值ε对算法性能的影响.我们选择不同的阈值ε={1,10−1,10−2,10−3,10−4,10−5,10−6},确保由PCA 构造子空间的维数,从1 维逐渐变到原始数据的维数,然后在不同的噪声水平(σ={0.1,0.2,0.3,0.4,0.5}) 下,单独运行TSSM 模型(λ=1) 100 次.计算结果如图14 所示,其中图14(a)为三个方向预测的统计结果,图14(b)为主方向预测的统计结果.每个子图中,横坐标表示不同的阈值ε,纵坐标依次为子空间维数、MSE 值、R2值和MAPE值.可以看到,当ε=1 时,此时PCA 构造的子空间只有1 维;随着ε的增大,子空间维数逐渐增大至原始数据的维数,即PCA 没有对原始数据进行降维.随着ε的增大,MSE 值逐渐减小,这主要是由于子空间映射时,维数越多,信息量损失越小,因此可以得到更加准确的预测结果.同时注意到,MSE值、R2值和MAPE值的变化范围非常小,特别是R2值和MAPE 值仅出现小范围的变化,可近似看作稳定.即使当PCA构造的子空间只有1 维时,TSSM模型也能取得令人满意的预测结果,这也说明PCA构造的子空间能够有效提取原始数据的主要成分,具有很强的稳定性.为了对比验证所提算法的有效性,本文在实验中选择ε=10−4,此时PCA 构造的子空间保留了原始数据绝大部分信息,信息损失量较少,获得预测 结果较为准确.

图14 PCA 中的阈值 ε 对TSSM 预测模型的影响Fig.14 The influence of ε in PCA on the prediction performance of TSSM

3 讨论和结论

呼吸运动是放射治疗(RT)中导致误差和不确定性的重要因素[1−2].本文首先采用三维图像分割技术,对4D CT 图像进行分割,在不使用标记物的情况下,通过提取膈肌的质心以准确度量膈肌的位移,使用ICP 算法计算胸腹部表面的位移.为了解决跨空间预测问题,本文提出了一种新的分步子空间映射算法(TSSM),通过构造特征子空间,并在高维子空间中进行映射,从而可以有效地提高跨空间数据预测的准确性.

对20 位患者真实数据的实验表明,本文提出的方法可以有效地描述膈肌和胸腹部表面的运动情况.相比于通过原始空间的直接预测,TSSM 算法在准确性和鲁棒性方面均优于经典线性的OLS 模型和非线性的ANN 模型.考虑到具有多项式核的kTSSM 算法比TSSM 算法具有更多的参数,容易出现过拟合现象,因此在实际应用中,我们建议直接应用TSSM 算法,线性模型已经可以获得令人满意的预测效果.值得一提的是,我们发现在病人吸气处于最高点的位置时,预测效果会更好.所提算法具有解析解,总计算时间(训练加上测试)在0.6 ms以下,可以达到实时处理速度,运行效率较高,这将有助于提高放射治疗中门控技术和跟踪技术的效率和精度.