基于高斯混合模型Varian NovalisTX直线加速器盆腔肿瘤放疗摆位误差分布预测模型的构建

丘敏敏，钟嘉健，欧阳斌，肖振华，邓永锦

（中山大学附属第一医院放射治疗科，广东广州510080）

放疗以其无创、精确、疗效确切的优点，成为恶性肿瘤治疗中的主要手段之一［1］，而调强放射治疗技术（intensity modulated radiation therapy，IMRT）在保证靶区剂量覆盖满足要求的同时，有效降低了周围危及器官剂量受量，在放疗中得到了广泛应用［2］。IMRT的实施要求患者在每次治疗中体位严格精确重复。放疗过程中肿瘤大小形状变化、患者体质改变、呼吸运动以及治疗师摆位习惯，都会引入治疗摆位误差，影响IMRT放疗的准确实施［3-5］。图像引导放射治疗（image guided radiotherapy，IGRT）技术的应用可较好地控制摆位误差。通过参考IGRT影像验证和调整患者放疗体位，在一定程度上提高了IMRT精确性［6］。然而IGRT所产生X线的剂量不容忽视。据研究，按不同部位采用不同的锥形束CT（cone beam computed tomography，CBCT）扫描协议，产生的剂量达0.1～2.0 cGy；而千伏级二维影像（2 dimensional kilo-volt image，2D-KV）按不同部位产生剂量只有0.1～1.2 cGy［7］。另一方面，IGRT技术必然会占用治疗机的一定时间。尤其是目前国内普遍存在设备不足、患者数量较多的情况下，较多数放疗单位只会每周1次甚至只在疗程第1分次做IGRT，很难保证每次IMRT治疗摆位精度。本文选取了2016年1月至2017年1月在中山大学附属第一医院放疗科Varian NovalisTX直线加速器治疗病例中，疗程每分次均做IGRT的30例盆腔肿瘤患者的临床资料及其IGRT摆位误差数据，采用高斯混合模型（Gaussian mixed model，GMM）建模方法，在此先验IGRT误差数据基础上建立盆腔肿瘤放疗摆位误差分布预测模型，对摆位误差进行了定量描述和预测分析，为缺少IGRT情况下盆腔肿瘤放疗摆位误差控制及肿瘤计划靶区外扩大小提供参考。

1 材料与方法

1.1 病例资料

选取中山大学附属第一医院放疗科2016年1月至2017年1月期间，采用Varian NovalisTX直线加速器放疗，且疗程内每分次均做IGRT的30例盆腔肿瘤患者病例。所有患者充分理解治疗及数据采集过程，同意参与数据采集，并均签署知情同意书。获取其临床资料和每分次治疗的IGRT摆位误差数据，其中20例直肠癌和10例宫颈癌，分期均为T2-4N0-2M0。20例直肠癌患者中15例为男性，5例为女性。患者年龄分布为32～89岁，中位年龄为54岁。所有患者体位均采用真空袋，头先进俯卧位固定，使用Philips 16排大孔径螺旋CT扫描，扫描层厚5 mm。由Varian Eclipse计划系统进行靶区勾画和计划设计，并在Varian NovalisTX直线加速器上实施放疗。患者疗程内IGRT包括每周1次CBCT和4次2D-KV，采用机载影像系统（On Board Image，OBI）获取每次治疗的摆位误差数据，按床升降方向（Vertical，Vrt）、进出方向（Longitudinal，Lng）和左右方向（Lateral，Lat）进行记录，共收集了770次IGRT误差数据进行模型构建。

1.2 盆腔肿瘤放疗摆位误差分布预测模型的构建

1.2.1 误差分布预测模型构建使用到的概念（1）单高斯模型高斯分布 有时也被称为正态分布（Normal distribution），是一种在自然界大量的存在的、最为常见的分布形式。高斯分布定义为：若随机变量X服从一个数学期望为μ、方差为σ2的高斯分布，则记为N（μ，σ2），其概率密度函数为：

式中参数μ表示均值，均值对应正态分布的中间位置，参数σ表示标准差。

单高斯模型在二维平面上的分布如图1。

图1 二维单高斯分布Fig.1 Two-dimensional single Gaussian distribution

（2）GMM模型 若事物分布由多个高斯分布线性组合而成，则需要引入GMM模型。其定义为，假设数据服从高斯混合分布，则其概率分布模型具有如下形式：

kkkk斯分布密度，θk=（μk，σk2），而

称为第k个模型。

理论上GMM模型可以拟合出任意类型的分布，通常用于解决同一集合下的数据包含多个不同的分布的情况。GMM模型在二维平面上的分布如图2。

图2 二维高斯混合分布Fig.2 Two-dimensional Gaussian mixture distribution

（3）EM（Expectation Maximization）算法 EM算法即期望最大化算法，在统计学中被用于寻找依赖于不可观察隐性变量的概率模型的参数最大似然估计，是一种解决存在隐含变量优化问题的有效方法。因GMM模型函数很难通过展开求偏导方式处理，优化问题麻烦，故通常采用EM算法求解其参数。EM算法为迭代型算法，在每一次的迭代过程分为求期望（expectationZ，E）步骤和最大化（maximization，M）步骤。E步骤假设知道各模型参数（可初始化或基于上一步迭代结果）估计每个高斯模型的隐变量；M步骤基于估计的隐变量返回确定模型参数，重复EM步骤直至收敛。具体算法实现如下：①将GMM模型改写成为有3个参数π，μ，ε待求解的形式：

②定义聚类数为K，对每个分类k设置πk，μk和εk的初始值，然后计算（4）式的对数似然函数；③E step：根据当前的πk，μk和εk计算后验概率γ（znk）：

④M step：根据E step计算的γ（znk）在计算新的πk，μk和εk：

⑤计算（4）式中的对数似然函数lnp（x│π，μ，ε）〗；⑥检查参数是否收敛或对数似然函数是否收敛，若不收敛则返回步骤3。

（4）K-means算法 means算法是一种常用的聚类算法，它认为由一组数据点构成的一个聚类中，聚类内部点之间的距离应该小于数据点与聚类外部的点之间的距离。在特定条件下K-means和GMM可互相用对方思想表达，K-means算法可视作GMM的一种特殊形式，而GMM模型提供了更强的描述能力。相比于K-means而言GMM每一次迭代计算量较大，故实际运用中通常先采用K-means算法获取初始聚类结果，然后将其聚类数和聚类中心作为初值传给GMM模型进行更细致迭代。K-means具体实现算法如下：①随机选取k个中心点；②遍历所有数据，将每个数据划分到最近的中心类内；③计算每个聚类的平均值，并作为新的中心点；④重复②-③，直到收敛或执行最大次数完成。其中Kmeans的聚类数k通常可用Elbow法进行确定，即采用枚举法将数据按从小到大的k值进行分类，获取每次分类的方差百分比（Percentage of variance explained）并作图，取曲线中拐点（Elbow）对应k值作为最佳聚类数。K-means初始聚类代码如下（R语言实现）：

Kmeans聚类，k=2，3，4...，取 elbow 的 k值为最佳聚类数，并获取其初始聚类中心my_km＜-kmeans（data，k）

1.2.2 误差分布预测模型构建过程 ①将患者IGRT摆位误差数据按Vrt、Lng、Lat 3个方向记录，并将其转化为三维矩阵数据保存和处理；②采用K-means算法对原始摆位误差数据进行聚类，获取其聚类数和初始聚类中心；③将K-means算法获取的聚类数和聚类中心作为初始值传递入GMM模型，采用EM算法进行迭代，求解确定GMM模型的参数，并分析该模型的临床意义。GMM模型的部分编程代码如下（MATLAB实现）：将K-means算法获取的聚类数和聚类中心作为初始值传递入GMM模型，采用EM算法进行迭代，获取摆位误差分布预测GMM模型。

while 1

%E步骤

Qt=E_step（data，mu，msigma，mp）；

%M步骤

［mu，msigma，mp］=M_step（Qt，data）；

loglik_nxt=loglike（data，mu，msigma，mp）；

if abs（（loglik_nxt/loglik_pre）-1）＜ loglik_threshold

break；

end

step=step+1；

%求似然值

loglik_pre=loglik_nxt；

end

1.3 数据处理软件

本研究误差数据使用Microsoft Office Excel 2007进行数据收集统计，R语言编程进行K-means初始聚类，MATLAB R2015a编程进行GMM模型的构建和参数求解。

2 结果

2.1 原始误差矩阵

30例盆腔肿瘤患者，总共采集770组摆位误差数据，统计结果如表1所示。

表1 原始摆位误差数据统计Table 1 Raw set-up errors statistics（mm）

将患者IGRT摆位误差数据按Vrt、Lng、Lat三个方向进行记录并将其转化为三维矩阵，得原始误差数据矩阵分布如图3所示。

2.2 最佳聚类数k

采用K-means法对原始三维矩阵误差数据进行聚类，通过Elbow法确定的最佳聚类数k值为4（图4）。

将误差数据按k=4进行聚类，可得初始聚类结果如图5所示。

图3 原始误差数据三维矩阵分布Fig.3 Three-dimensional matrix distribution of raw set-up errors

图4 Elbow法确定最佳聚类数k值Fig.4 Select the optimal number of clusters（k value）by Elbow method

图5 K-means初始聚类效果图Fig.5 Preliminary clustering results by K-means

表2 摆位误差中心点坐标Table 2 Center point coordinates ofset-up errors

2.3 摆位误差分布预测GMM模型

将K-means算法获取的聚类数和聚类中心作为初始值传递入GMM模型，采用EM法进行迭代，获取摆位误差分布预测GMM模型，其聚类效果及聚类中心分布如图6所示。

求解得到的GMM误差分布预测模型参数如下：

图6 GMM聚类效果Fig.6 Clustering results by GMM

各误差中心坐标（即GMM模型的均值μ）如表2所示；误差模型的协方差矩阵（即GMM模型σ）如表3所示；各误差中心概率（即GMM模型的系数α）如表4所示。

表3 GMM模型协方差Table 3 Covariance of GMM model

表4 摆位误差中心点概率Table 4 Center point probability of set-up errors

2.4 GMM模型的参数的分析

摆位误差主要往4个中心点（μ1～μ4）的方向集中。各个中心的空间坐标值可反映该中心内的点平均偏移方向及偏移量，如μ1的Vrt为4.28 mm，表明摆位在Vrt方向平均偏移量较大；同理μ2的Lat为-0.25 mm，表明该中心内点在Lat方向偏移量较小。从总体中心分布数据上可以看出，所有摆位在Vrt方向偏移（-3.88～4.28 mm）和Lng方向偏移（-2.41～1.54 mm）都较大，而Lat方向上偏移较小（-1.85～0.72 mm）。摆位误差中心的概率（系数α）反映了误差分布落于该中心方向及附近的可能性。由个中心概率可知，摆位误差往μ 2和μ4方向偏移可能性（0.301、0.310）较μ1和μ3的（0.190、0.196）更大。协方差矩阵反映了标准差的大小。据该GMM模型协方差参数可知，摆位误差的统计标准差可达5.2 mm。

3 讨论

GMM建模法通过将一个事物分解为若干基于高斯概率密度函数（正态分布曲线）的方式构建模型，以达到精确量化地描述事物特性的目的。理论上无论观测数据集以何种规律分布，都可以通过由单一高斯模型线性组合的GMM模型进行拟合［8-9］。本文将盆腔肿瘤放疗先验IGRT摆位误差数据集视为高斯混合形式分布的前提下，通过统计建模方法求解其GMM模型分布函数，实现对误差数据分布规律的定量描述和预测分析，为临床摆位和治疗误差控制提供了参考。由本文GMM模型参数可知：摆位误差主要往μ1～μ4 4个中心点方向集中；中心坐标值表明在Vrt和Lng方向误差偏移较大，Lat方向偏移较小；从概率上分析误差往μ2和μ4方向偏移的可能性较μ1和μ3大。

据 Stroom［10］及 Van Herk［11］的研究定义，治疗摆位误差包括每个患者放疗分次间及分次内Vrt、Lng及Lat各轴向误差，同时又分为系统误差和随机误差。其中系统误差为所有分次摆位误差的平均值，随机误差为所有分次摆位误差的标准差。在此理论基础上，本文通过统计学及建模方法，对先验IGRT误差数据分布进行更深入研究，认为摆位偏差不只是简单的在三个轴向方向上误差，而是更趋向于向空间内几个确定中心方向集中偏移分布，且各中心方向上的偏移分布概率不同。

治疗过程体位的固定方式、肿瘤大小和形状变化、患者体质变化、呼吸运动以及治疗师摆位习惯不同，都会引入治疗摆位误差［3-5］。尤其在盆腔肿瘤的摆位治疗中，因盆腔生理结构以非刚性肌肉及内脏器为主，且其骨性结构活动度较大，无论采用何种体位固定方式都难以保证体位的高度重复；同时因呼吸运动引起的胸腹部运动，以及盆腔内肠道器官蠕动等，均致使盆腔部位治疗摆位精度较其他部位更低。宋伟男等［12］报道直肠癌调强放疗时摆位误在左右、头脚、前后方向分别为1.1±2.3、2.1±5.0、-1.1±2.2，许峰等［13］通过CBCT图像与计划CT图像配准的方法得到盆腹部靶中心左右、头脚、前后方向的误差为-0.8±2.1、-0.3±5.9、0.1±2.6。本文选取分析的先验IGRT误差数据也表现了接近的结果。

IMRT技术的广泛应用保证了靶区剂量覆盖满足要求的同时，降低了周围危及器官剂量受量，其顺利实施要求患者在每次治疗中体位严格精确重复。IGRT技术可较好地控制摆位误差，提高分次IMRT精确性。然而IGRT所产生X线的剂量不容忽视，同时也导致放疗工作效率下降。本文通过对完整IGRT治疗的患者先验误差数据进行统计建模分析，呈现误差分布规律并做误差概率预测，给出了盆腔部肿瘤放疗摆位中最可能偏移的方向及误差大小，可为缺少IGRT的日常摆位误差控制提供参考。

另一方面，PTV外扩边界的确定是放射治疗中的关键问题，合理的PTV边界既要保证包含靶区的可能运动区域，又要尽可能的降低靶区附近正常组织的器官受量。因此摆位误差是确定CTV到PTV外放距离的重要因素［14-15］。本文的研究结果认为，PTV的外放不能单纯的从Vrt、Lng和Lat三个轴向上考虑，更应该在其4个偏移中心的方向及其方差综合进行外扩，具体而言，需在各中心方向上采用不均匀外扩，且将方差偏移量包含进去。

本文模型尚存在待完善之处。首先模型的原始数据取自Varian NovalisTX直线加速器的IGRT误差数据，其最小精度为1mm，且没有考虑角度旋转误差。其次所有病例的体位固定方式均为真空袋固定俯卧位，其它体位固定的盆腔摆位误差不能确定由该GMM模型统一描述预测，但该模型建模方法可为任何部位误差分析提供参考。此外因临床数据有限，本文仅采集了30例病例进行分析，日后可收集更多数据进行验证完善。

IMRT技术的广泛应用保证了靶区剂量覆盖满足要求的同时，降低了周围危及器官剂量受量，其顺利实施要求患者在每次治疗中体位严格精确重复。IGRT技术可较好地控制摆位误差并提高分次IMRT精确性，然而IGRT所产生X线的剂量不容忽视，同时也会占用机器较多时间。本文选取了Varian NovalisTX直线加速器治疗的30例盆腔肿瘤患者临床资料及其IGRT摆位误差数据，采用GMM建模法构建盆腔肿瘤放疗摆位误差分布预测模型，对摆位误差分布进行了定量描述和预测分析，为缺少IGRT情况下盆腔肿瘤放疗摆位误差控制及肿瘤计划靶区外扩大小提供参考。