基于高斯混合模型Varian NovalisTX直線加速器盆腔腫瘤放療擺位誤差分布預測模型的構建

丘敏敏，鐘嘉健，歐陽斌，肖振華，鄧永錦

（中山大學附屬第一醫(yī)院放射治療科，廣東廣州510080）

放療以其無創(chuàng)、精確、療效確切的優(yōu)點，成為惡性腫瘤治療中的主要手段之一［1］，而調(diào)強放射治療技術（intensity modulated radiation therapy，IMRT）在保證靶區(qū)劑量覆蓋滿足要求的同時，有效降低了周圍危及器官劑量受量，在放療中得到了廣泛應用［2］。IMRT的實施要求患者在每次治療中體位嚴格精確重復。放療過程中腫瘤大小形狀變化、患者體質(zhì)改變、呼吸運動以及治療師擺位習慣，都會引入治療擺位誤差，影響IMRT放療的準確實施［3-5］。圖像引導放射治療（image guided radiotherapy，IGRT）技術的應用可較好地控制擺位誤差。通過參考IGRT影像驗證和調(diào)整患者放療體位，在一定程度上提高了IMRT精確性［6］。然而IGRT所產(chǎn)生X線的劑量不容忽視。據(jù)研究，按不同部位采用不同的錐形束CT（cone beam computed tomography，CBCT）掃描協(xié)議，產(chǎn)生的劑量達0.1～2.0 cGy；而千伏級二維影像（2 dimensional kilo-volt image，2D-KV）按不同部位產(chǎn)生劑量只有0.1～1.2 cGy［7］。另一方面，IGRT技術必然會占用治療機的一定時間。尤其是目前國內(nèi)普遍存在設備不足、患者數(shù)量較多的情況下，較多數(shù)放療單位只會每周1次甚至只在療程第1分次做IGRT，很難保證每次IMRT治療擺位精度。本文選取了2016年1月至2017年1月在中山大學附屬第一醫(yī)院放療科Varian NovalisTX直線加速器治療病例中，療程每分次均做IGRT的30例盆腔腫瘤患者的臨床資料及其IGRT擺位誤差數(shù)據(jù)，采用高斯混合模型（Gaussian mixed model，GMM）建模方法，在此先驗IGRT誤差數(shù)據(jù)基礎上建立盆腔腫瘤放療擺位誤差分布預測模型，對擺位誤差進行了定量描述和預測分析，為缺少IGRT情況下盆腔腫瘤放療擺位誤差控制及腫瘤計劃靶區(qū)外擴大小提供參考。

1 材料與方法

1.1 病例資料

選取中山大學附屬第一醫(yī)院放療科2016年1月至2017年1月期間，采用Varian NovalisTX直線加速器放療，且療程內(nèi)每分次均做IGRT的30例盆腔腫瘤患者病例。所有患者充分理解治療及數(shù)據(jù)采集過程，同意參與數(shù)據(jù)采集，并均簽署知情同意書。獲取其臨床資料和每分次治療的IGRT擺位誤差數(shù)據(jù)，其中20例直腸癌和10例宮頸癌，分期均為T2-4N0-2M0。20例直腸癌患者中15例為男性，5例為女性?；颊吣挲g分布為32～89歲，中位年齡為54歲。所有患者體位均采用真空袋，頭先進俯臥位固定，使用Philips 16排大孔徑螺旋CT掃描，掃描層厚5 mm。由Varian Eclipse計劃系統(tǒng)進行靶區(qū)勾畫和計劃設計，并在Varian NovalisTX直線加速器上實施放療。患者療程內(nèi)IGRT包括每周1次CBCT和4次2D-KV，采用機載影像系統(tǒng)（On Board Image，OBI）獲取每次治療的擺位誤差數(shù)據(jù)，按床升降方向（Vertical，Vrt）、進出方向（Longitudinal，Lng）和左右方向（Lateral，Lat）進行記錄，共收集了770次IGRT誤差數(shù)據(jù)進行模型構建。

1.2 盆腔腫瘤放療擺位誤差分布預測模型的構建

1.2.1 誤差分布預測模型構建使用到的概念（1）單高斯模型高斯分布 有時也被稱為正態(tài)分布（Normal distribution），是一種在自然界大量的存在的、最為常見的分布形式。高斯分布定義為：若隨機變量X服從一個數(shù)學期望為μ、方差為σ2的高斯分布，則記為N（μ，σ2），其概率密度函數(shù)為：

式中參數(shù)μ表示均值，均值對應正態(tài)分布的中間位置，參數(shù)σ表示標準差。

單高斯模型在二維平面上的分布如圖1。

圖1 二維單高斯分布Fig.1 Two-dimensional single Gaussian distribution

（2）GMM模型 若事物分布由多個高斯分布線性組合而成，則需要引入GMM模型。其定義為，假設數(shù)據(jù)服從高斯混合分布，則其概率分布模型具有如下形式：

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

稱為第k個模型。

理論上GMM模型可以擬合出任意類型的分布，通常用于解決同一集合下的數(shù)據(jù)包含多個不同的分布的情況。GMM模型在二維平面上的分布如圖2。

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

（3）EM（Expectation Maximization）算法 EM算法即期望最大化算法，在統(tǒng)計學中被用于尋找依賴于不可觀察隱性變量的概率模型的參數(shù)最大似然估計，是一種解決存在隱含變量優(yōu)化問題的有效方法。因GMM模型函數(shù)很難通過展開求偏導方式處理，優(yōu)化問題麻煩，故通常采用EM算法求解其參數(shù)。EM算法為迭代型算法，在每一次的迭代過程分為求期望（expectationZ，E）步驟和最大化（maximization，M）步驟。E步驟假設知道各模型參數(shù)（可初始化或基于上一步迭代結果）估計每個高斯模型的隱變量；M步驟基于估計的隱變量返回確定模型參數(shù)，重復EM步驟直至收斂。具體算法實現(xiàn)如下：①將GMM模型改寫成為有3個參數(shù)π，μ，ε待求解的形式：

②定義聚類數(shù)為K，對每個分類k設置πk，μk和εk的初始值，然后計算（4）式的對數(shù)似然函數(shù)；③E step：根據(jù)當前的πk，μk和εk計算后驗概率γ（znk）：

④M step：根據(jù)E step計算的γ（znk）在計算新的πk，μk和εk：

⑤計算（4）式中的對數(shù)似然函數(shù)lnp（x│π，μ，ε）〗；⑥檢查參數(shù)是否收斂或?qū)?shù)似然函數(shù)是否收斂，若不收斂則返回步驟3。

（4）K-means算法 means算法是一種常用的聚類算法，它認為由一組數(shù)據(jù)點構成的一個聚類中，聚類內(nèi)部點之間的距離應該小于數(shù)據(jù)點與聚類外部的點之間的距離。在特定條件下K-means和GMM可互相用對方思想表達，K-means算法可視作GMM的一種特殊形式，而GMM模型提供了更強的描述能力。相比于K-means而言GMM每一次迭代計算量較大，故實際運用中通常先采用K-means算法獲取初始聚類結果，然后將其聚類數(shù)和聚類中心作為初值傳給GMM模型進行更細致迭代。K-means具體實現(xiàn)算法如下：①隨機選取k個中心點；②遍歷所有數(shù)據(jù)，將每個數(shù)據(jù)劃分到最近的中心類內(nèi)；③計算每個聚類的平均值，并作為新的中心點；④重復②-③，直到收斂或執(zhí)行最大次數(shù)完成。其中Kmeans的聚類數(shù)k通?？捎肊lbow法進行確定，即采用枚舉法將數(shù)據(jù)按從小到大的k值進行分類，獲取每次分類的方差百分比（Percentage of variance explained）并作圖，取曲線中拐點（Elbow）對應k值作為最佳聚類數(shù)。K-means初始聚類代碼如下（R語言實現(xiàn)）：

Kmeans聚類，k=2，3，4...，取 elbow 的 k值為最佳聚類數(shù)，并獲取其初始聚類中心my_km＜-kmeans（data，k）

1.2.2 誤差分布預測模型構建過程 ①將患者IGRT擺位誤差數(shù)據(jù)按Vrt、Lng、Lat 3個方向記錄，并將其轉化為三維矩陣數(shù)據(jù)保存和處理；②采用K-means算法對原始擺位誤差數(shù)據(jù)進行聚類，獲取其聚類數(shù)和初始聚類中心；③將K-means算法獲取的聚類數(shù)和聚類中心作為初始值傳遞入GMM模型，采用EM算法進行迭代，求解確定GMM模型的參數(shù)，并分析該模型的臨床意義。GMM模型的部分編程代碼如下（MATLAB實現(xiàn)）：將K-means算法獲取的聚類數(shù)和聚類中心作為初始值傳遞入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 數(shù)據(jù)處理軟件

本研究誤差數(shù)據(jù)使用Microsoft Office Excel 2007進行數(shù)據(jù)收集統(tǒng)計，R語言編程進行K-means初始聚類，MATLAB R2015a編程進行GMM模型的構建和參數(shù)求解。

2 結果

2.1 原始誤差矩陣

30例盆腔腫瘤患者，總共采集770組擺位誤差數(shù)據(jù)，統(tǒng)計結果如表1所示。

表1 原始擺位誤差數(shù)據(jù)統(tǒng)計Table 1 Raw set-up errors statistics（mm）

將患者IGRT擺位誤差數(shù)據(jù)按Vrt、Lng、Lat三個方向進行記錄并將其轉化為三維矩陣，得原始誤差數(shù)據(jù)矩陣分布如圖3所示。

2.2 最佳聚類數(shù)k

采用K-means法對原始三維矩陣誤差數(shù)據(jù)進行聚類，通過Elbow法確定的最佳聚類數(shù)k值為4（圖4）。

將誤差數(shù)據(jù)按k=4進行聚類，可得初始聚類結果如圖5所示。

圖3 原始誤差數(shù)據(jù)三維矩陣分布Fig.3 Three-dimensional matrix distribution of raw set-up errors

圖4 Elbow法確定最佳聚類數(shù)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算法獲取的聚類數(shù)和聚類中心作為初始值傳遞入GMM模型，采用EM法進行迭代，獲取擺位誤差分布預測GMM模型，其聚類效果及聚類中心分布如圖6所示。

求解得到的GMM誤差分布預測模型參數(shù)如下：

圖6 GMM聚類效果Fig.6 Clustering results by GMM

各誤差中心坐標（即GMM模型的均值μ）如表2所示；誤差模型的協(xié)方差矩陣（即GMM模型σ）如表3所示；各誤差中心概率（即GMM模型的系數(shù)α）如表4所示。

表3 GMM模型協(xié)方差Table 3 Covariance of GMM model

表4 擺位誤差中心點概率Table 4 Center point probability of set-up errors

2.4 GMM模型的參數(shù)的分析

擺位誤差主要往4個中心點（μ1～μ4）的方向集中。各個中心的空間坐標值可反映該中心內(nèi)的點平均偏移方向及偏移量，如μ1的Vrt為4.28 mm，表明擺位在Vrt方向平均偏移量較大；同理μ2的Lat為-0.25 mm，表明該中心內(nèi)點在Lat方向偏移量較小。從總體中心分布數(shù)據(jù)上可以看出，所有擺位在Vrt方向偏移（-3.88～4.28 mm）和Lng方向偏移（-2.41～1.54 mm）都較大，而Lat方向上偏移較?。?1.85～0.72 mm）。擺位誤差中心的概率（系數(shù)α）反映了誤差分布落于該中心方向及附近的可能性。由個中心概率可知，擺位誤差往μ 2和μ4方向偏移可能性（0.301、0.310）較μ1和μ3的（0.190、0.196）更大。協(xié)方差矩陣反映了標準差的大小。據(jù)該GMM模型協(xié)方差參數(shù)可知，擺位誤差的統(tǒng)計標準差可達5.2 mm。

3 討論

GMM建模法通過將一個事物分解為若干基于高斯概率密度函數(shù)（正態(tài)分布曲線）的方式構建模型，以達到精確量化地描述事物特性的目的。理論上無論觀測數(shù)據(jù)集以何種規(guī)律分布，都可以通過由單一高斯模型線性組合的GMM模型進行擬合［8-9］。本文將盆腔腫瘤放療先驗IGRT擺位誤差數(shù)據(jù)集視為高斯混合形式分布的前提下，通過統(tǒng)計建模方法求解其GMM模型分布函數(shù)，實現(xiàn)對誤差數(shù)據(jù)分布規(guī)律的定量描述和預測分析，為臨床擺位和治療誤差控制提供了參考。由本文GMM模型參數(shù)可知：擺位誤差主要往μ1～μ4 4個中心點方向集中；中心坐標值表明在Vrt和Lng方向誤差偏移較大，Lat方向偏移較?。粡母怕噬戏治稣`差往μ2和μ4方向偏移的可能性較μ1和μ3大。

據(jù) Stroom［10］及 Van Herk［11］的研究定義，治療擺位誤差包括每個患者放療分次間及分次內(nèi)Vrt、Lng及Lat各軸向誤差，同時又分為系統(tǒng)誤差和隨機誤差。其中系統(tǒng)誤差為所有分次擺位誤差的平均值，隨機誤差為所有分次擺位誤差的標準差。在此理論基礎上，本文通過統(tǒng)計學及建模方法，對先驗IGRT誤差數(shù)據(jù)分布進行更深入研究，認為擺位偏差不只是簡單的在三個軸向方向上誤差，而是更趨向于向空間內(nèi)幾個確定中心方向集中偏移分布，且各中心方向上的偏移分布概率不同。

治療過程體位的固定方式、腫瘤大小和形狀變化、患者體質(zhì)變化、呼吸運動以及治療師擺位習慣不同，都會引入治療擺位誤差［3-5］。尤其在盆腔腫瘤的擺位治療中，因盆腔生理結構以非剛性肌肉及內(nèi)臟器為主，且其骨性結構活動度較大，無論采用何種體位固定方式都難以保證體位的高度重復；同時因呼吸運動引起的胸腹部運動，以及盆腔內(nèi)腸道器官蠕動等，均致使盆腔部位治療擺位精度較其他部位更低。宋偉男等［12］報道直腸癌調(diào)強放療時擺位誤在左右、頭腳、前后方向分別為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誤差數(shù)據(jù)也表現(xiàn)了接近的結果。

IMRT技術的廣泛應用保證了靶區(qū)劑量覆蓋滿足要求的同時，降低了周圍危及器官劑量受量，其順利實施要求患者在每次治療中體位嚴格精確重復。IGRT技術可較好地控制擺位誤差，提高分次IMRT精確性。然而IGRT所產(chǎn)生X線的劑量不容忽視，同時也導致放療工作效率下降。本文通過對完整IGRT治療的患者先驗誤差數(shù)據(jù)進行統(tǒng)計建模分析，呈現(xiàn)誤差分布規(guī)律并做誤差概率預測，給出了盆腔部腫瘤放療擺位中最可能偏移的方向及誤差大小，可為缺少IGRT的日常擺位誤差控制提供參考。

另一方面，PTV外擴邊界的確定是放射治療中的關鍵問題，合理的PTV邊界既要保證包含靶區(qū)的可能運動區(qū)域，又要盡可能的降低靶區(qū)附近正常組織的器官受量。因此擺位誤差是確定CTV到PTV外放距離的重要因素［14-15］。本文的研究結果認為，PTV的外放不能單純的從Vrt、Lng和Lat三個軸向上考慮，更應該在其4個偏移中心的方向及其方差綜合進行外擴，具體而言，需在各中心方向上采用不均勻外擴，且將方差偏移量包含進去。

本文模型尚存在待完善之處。首先模型的原始數(shù)據(jù)取自Varian NovalisTX直線加速器的IGRT誤差數(shù)據(jù)，其最小精度為1mm，且沒有考慮角度旋轉誤差。其次所有病例的體位固定方式均為真空袋固定俯臥位，其它體位固定的盆腔擺位誤差不能確定由該GMM模型統(tǒng)一描述預測，但該模型建模方法可為任何部位誤差分析提供參考。此外因臨床數(shù)據(jù)有限，本文僅采集了30例病例進行分析，日后可收集更多數(shù)據(jù)進行驗證完善。

IMRT技術的廣泛應用保證了靶區(qū)劑量覆蓋滿足要求的同時，降低了周圍危及器官劑量受量，其順利實施要求患者在每次治療中體位嚴格精確重復。IGRT技術可較好地控制擺位誤差并提高分次IMRT精確性，然而IGRT所產(chǎn)生X線的劑量不容忽視，同時也會占用機器較多時間。本文選取了Varian NovalisTX直線加速器治療的30例盆腔腫瘤患者臨床資料及其IGRT擺位誤差數(shù)據(jù)，采用GMM建模法構建盆腔腫瘤放療擺位誤差分布預測模型，對擺位誤差分布進行了定量描述和預測分析，為缺少IGRT情況下盆腔腫瘤放療擺位誤差控制及腫瘤計劃靶區(qū)外擴大小提供參考。