丘敏敏，鐘嘉健，歐陽(yáng)斌，肖振華，鄧永錦

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

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

1 材料與方法

1.1 病例資料

選取中山大學(xué)附屬第一醫(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歲。所有患者體位均采用真空袋，頭先進(jìn)俯臥位固定，使用Philips 16排大孔徑螺旋CT掃描，掃描層厚5 mm。由Varian Eclipse計(jì)劃系統(tǒng)進(jìn)行靶區(qū)勾畫和計(jì)劃設(shè)計(jì)，并在Varian NovalisTX直線加速器上實(shí)施放療?；颊忒煶虄?nèi)IGRT包括每周1次CBCT和4次2D-KV，采用機(jī)載影像系統(tǒng)（On Board Image，OBI）獲取每次治療的擺位誤差數(shù)據(jù)，按床升降方向（Vertical，Vrt）、進(jìn)出方向（Longitudinal，Lng）和左右方向（Lateral，Lat）進(jìn)行記錄，共收集了770次IGRT誤差數(shù)據(jù)進(jìn)行模型構(gòu)建。

1.2 盆腔腫瘤放療擺位誤差分布預(yù)測(cè)模型的構(gòu)建

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

式中參數(shù)μ表示均值，均值對(duì)應(yīng)正態(tài)分布的中間位置，參數(shù)σ表示標(biāo)準(zhǔn)差。

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

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

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

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

稱為第k個(gè)模型。

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

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

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

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

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

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

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

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

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

2 結(jié)果

2.1 原始誤差矩陣

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

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

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

2.2 最佳聚類數(shù)k

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

將誤差數(shù)據(jù)按k=4進(jìn)行聚類，可得初始聚類結(jié)果如圖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 擺位誤差中心點(diǎn)坐標(biāo)Table 2 Center point coordinates ofset-up errors

2.3 擺位誤差分布預(yù)測(cè)GMM模型

將K-means算法獲取的聚類數(shù)和聚類中心作為初始值傳遞入GMM模型，采用EM法進(jìn)行迭代，獲取擺位誤差分布預(yù)測(cè)GMM模型，其聚類效果及聚類中心分布如圖6所示。

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

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

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

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

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

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

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

3 討論

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

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

治療過程體位的固定方式、腫瘤大小和形狀變化、患者體質(zhì)變化、呼吸運(yùn)動(dòng)以及治療師擺位習(xí)慣不同，都會(huì)引入治療擺位誤差［3-5］。尤其在盆腔腫瘤的擺位治療中，因盆腔生理結(jié)構(gòu)以非剛性肌肉及內(nèi)臟器為主，且其骨性結(jié)構(gòu)活動(dòng)度較大，無論采用何種體位固定方式都難以保證體位的高度重復(fù)；同時(shí)因呼吸運(yùn)動(dòng)引起的胸腹部運(yùn)動(dòng)，以及盆腔內(nèi)腸道器官蠕動(dòng)等，均致使盆腔部位治療擺位精度較其他部位更低。宋偉男等［12］報(bào)道直腸癌調(diào)強(qiáng)放療時(shí)擺位誤在左右、頭腳、前后方向分別為1.1±2.3、2.1±5.0、-1.1±2.2，許峰等［13］通過CBCT圖像與計(jì)劃CT圖像配準(zhǔn)的方法得到盆腹部靶中心左右、頭腳、前后方向的誤差為-0.8±2.1、-0.3±5.9、0.1±2.6。本文選取分析的先驗(yàn)IGRT誤差數(shù)據(jù)也表現(xiàn)了接近的結(jié)果。

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

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

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

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