基于分步子空間映射的無(wú)標(biāo)記膈肌運(yùn)動(dòng)預(yù)測(cè)算法

余 航 李晨陽(yáng) 余紹德 馮冬竹 許錄平

呼吸運(yùn)動(dòng)(Respiratory motion)是放射治療(Radiation therapy,RT)中,導(dǎo)致誤差和不確定性的重要因素[1?2].一方面,現(xiàn)代放射治療技術(shù)可以達(dá)到毫米級(jí)的傳輸精度,如調(diào)強(qiáng)放療(Intensity-modulated radiation therapy,IMRT)[3];另一方面,呼吸運(yùn)動(dòng)可以導(dǎo)致胸部和腹部的腫瘤移動(dòng)多達(dá)35 毫米[4].在此情況下[5],如果不對(duì)呼吸運(yùn)動(dòng)進(jìn)行有效控制和補(bǔ)償,要么腫瘤得不到足夠的放射計(jì)量,要么腫瘤周圍的健康組織將會(huì)受到放射的傷害[6?7].因此,有效減少呼吸對(duì)器官和腫瘤運(yùn)動(dòng)的影響,對(duì)整個(gè)放射治療的進(jìn)程和效果顯得尤為重要[8?10].

針對(duì)這一問(wèn)題,研究者提出了多種不同的解決方法,包括:1)屏氣法[11?13];2)淺呼吸法[11?19];3)呼吸門控法[15?16];4)動(dòng)態(tài)目標(biāo)跟蹤法[17?21].呼吸門控法和動(dòng)態(tài)目標(biāo)跟蹤法都需要實(shí)時(shí)了解腫瘤的位置.呼吸門控法通過(guò)在呼吸周期的預(yù)定階段打開(kāi)放射光束,從而減少正常組織的照射量.動(dòng)態(tài)目標(biāo)跟蹤法可以讓光束和腫瘤進(jìn)行同步的運(yùn)動(dòng),實(shí)現(xiàn)連續(xù)照射,是目前最為有效的呼吸管理技術(shù).腫瘤的位置可以通過(guò)直接方式和間接方式來(lái)獲得.直接方式通過(guò)放射線成像系統(tǒng)對(duì)腫瘤或植入體內(nèi)的標(biāo)志物進(jìn)行成像,該方法可以獲得人體內(nèi)的結(jié)構(gòu)信息,但會(huì)給患者帶來(lái)額外的輻射負(fù)擔(dān)[10,22?27],并且需要進(jìn)行昂貴的、創(chuàng)傷性手術(shù),因此目前還沒(méi)有廣泛應(yīng)用.

間接方式是使用呼吸信息對(duì)內(nèi)部解剖結(jié)構(gòu)的運(yùn)動(dòng)進(jìn)行建模和預(yù)測(cè).Hoisak 等[27]通過(guò)肺活量測(cè)定和實(shí)時(shí)位置跟蹤系統(tǒng),研究肺部腫瘤運(yùn)動(dòng)與呼吸運(yùn)動(dòng)的相關(guān)性.該方法的局限是,只使用了線性相關(guān)系數(shù)來(lái)度量腫瘤和參考物在一個(gè)方向上的運(yùn)動(dòng)相關(guān)情況.Paganelli 等[8]從核磁共振圖像(Magnetic resonance imaging,MRI)中提取SIFT (Scale-invariant feature transform)圖像特征,對(duì)基于圖像特征的跟蹤方法與基于標(biāo)志物的跟蹤方法進(jìn)行對(duì)比.該方法采用二次多項(xiàng)式進(jìn)行建模,需要在人體表面放置標(biāo)志物,并計(jì)算45個(gè)圖像特征的均值,較為繁瑣.文獻(xiàn)[9]通過(guò)實(shí)驗(yàn)表明,相比于其它解剖結(jié)構(gòu),膈肌與呼吸運(yùn)動(dòng)具有更強(qiáng)的相關(guān)性.根據(jù)這一結(jié)論,很多研究者采用膈肌作為參考點(diǎn),對(duì)腫瘤運(yùn)動(dòng)情況進(jìn)行研究.文獻(xiàn)[28]針對(duì)不同患者,提出了一種基于膈肌的患者特異性呼吸模型,取得了良好的實(shí)驗(yàn)結(jié)果.

另一方面,人體內(nèi)部結(jié)構(gòu)的運(yùn)動(dòng)情況,可以通過(guò)機(jī)體外部的觀察來(lái)進(jìn)行預(yù)測(cè).最近的研究[9]表明,人體整個(gè)外表面的運(yùn)動(dòng)與內(nèi)部結(jié)構(gòu)的運(yùn)動(dòng)具有很強(qiáng)的相關(guān)性.據(jù)此,文獻(xiàn)[29?30]使用K-means 算法對(duì)人體胸腹部表面的運(yùn)動(dòng)軌跡進(jìn)行聚類,以產(chǎn)生呼吸運(yùn)動(dòng)的參考點(diǎn).但是由于K-means 算法易受到初始值和數(shù)據(jù)結(jié)構(gòu)的影響,該方法容易受到噪聲的干擾.Malinowski 等[1]將三個(gè)標(biāo)記物固定在患者貼身背心上,據(jù)此來(lái)研究人體表面與腫瘤位置的相互運(yùn)動(dòng)關(guān)系.Ozhasoglu 等[4]通過(guò)使用光學(xué)監(jiān)測(cè)裝置觀察附著在胸部和腹部的紅外發(fā)光二極管(LED),來(lái)研究呼吸運(yùn)動(dòng)對(duì)肺和胰腺腫瘤運(yùn)動(dòng)的影響.這兩種方法[1,4]均需在人體外表面放置標(biāo)記物,因此較為繁瑣,易受到操作者的影響.

對(duì)體內(nèi)膈肌和體外胸腹表面之間關(guān)系的準(zhǔn)確、系統(tǒng)的研究,不僅可以提高呼吸運(yùn)動(dòng)管理的準(zhǔn)確性,而且可以顯著降低放射治療的復(fù)雜性,但相關(guān)工作仍然較少,且要么需要患者佩戴標(biāo)記物,較為繁瑣,要么模型建立比較簡(jiǎn)單.考慮到這兩個(gè)器官具有不同的運(yùn)動(dòng)模式,其數(shù)據(jù)具有不同的分布結(jié)構(gòu),形成不同的子空間,無(wú)法直接進(jìn)行準(zhǔn)確的映射.針對(duì)這一問(wèn)題,本文提出了一種新的分步子空間映射(Twostep subspace mapping,TSSM)算法,通過(guò)體外無(wú)標(biāo)記胸腹部表面的測(cè)量來(lái)預(yù)測(cè)體內(nèi)膈肌的運(yùn)動(dòng)情況.該方法無(wú)需患者佩戴標(biāo)記物,且采用多種線性、非線性模型,系統(tǒng)地對(duì)體內(nèi)膈肌和體外胸腹表面運(yùn)動(dòng)情況進(jìn)行對(duì)比研究.本文首先針對(duì)4D CT 圖像,利用三維圖像分割方法,對(duì)體內(nèi)膈肌和體外胸腹部表面的位移進(jìn)行準(zhǔn)確的測(cè)量.為了解決域自適應(yīng)(Domain adaptation)問(wèn)題,TSSM 首先利用主成分分析法(Principal component analysis,PCA)為每組數(shù)據(jù)構(gòu)建特征子空間.然后通過(guò)線性嶺回歸(Ridge regression)優(yōu)化過(guò)程,將膈肌數(shù)據(jù)與胸腹部表面數(shù)據(jù)連接起來(lái),得到一個(gè)子空間映射模型.基于得到的映射模型,該方法可以通過(guò)體外胸腹部表面數(shù)據(jù),對(duì)體內(nèi)膈肌的運(yùn)動(dòng)情況進(jìn)行預(yù)測(cè).為了研究數(shù)據(jù)的非線性相關(guān)性,本文進(jìn)一步將TSSM 擴(kuò)展到基于核的分步子空間映射算法(Kernel TSSM,kTSSM),包括多項(xiàng)式核和高斯核.本文給出了兩種算法的解析解,其不需要遞歸迭代的操作,具有運(yùn)算速度快的特點(diǎn).對(duì)比實(shí)驗(yàn)表明,這種分步映射的策略可以有效解決跨空間數(shù)據(jù)的預(yù)測(cè)問(wèn)題,具有很好的準(zhǔn)確性和魯棒性,優(yōu)于經(jīng)典的線性模型和ANN模型,本文也進(jìn)一步給出了內(nèi)部膈肌與外部胸腹部表面之間的運(yùn)動(dòng)關(guān)系模型.

本文的其余部分安排如下:第1 節(jié)具體介紹了所提方法和數(shù)據(jù),實(shí)驗(yàn)在第2 節(jié)中介紹,第3 節(jié)給出結(jié)論和討論.

1 方法

1.1 內(nèi)部膈肌位移提取

當(dāng)獲得患者的4D CT 圖像后,本文設(shè)計(jì)了一種三維圖像分割方法對(duì)其進(jìn)行分割.在4D CT 圖像中,由于具有不同的組織結(jié)構(gòu),身體區(qū)域的體素值(Voxel)均大于800,而背景區(qū)域和肺部區(qū)域主要為空氣,其體素值在50 以下.考慮到圖像的背景、身體區(qū)域和肺部區(qū)域具有顯著的灰度差別,本文采用閾值法對(duì)其進(jìn)行初步分割.具體地,本文設(shè)置閾值(本文設(shè)為500),首先將身體區(qū)域從背景中分割出來(lái);然后根據(jù)區(qū)域面積,將圖像中的最大目標(biāo)區(qū)域作為身體區(qū)域;最后進(jìn)一步通過(guò)形態(tài)學(xué)方法,去除背景中的孤立區(qū)域并填充身體區(qū)域中的孔,從而得到身體區(qū)域的分割結(jié)果.在獲得身體區(qū)域之后,由于肺部區(qū)域在身體區(qū)域內(nèi)部,并與身體區(qū)域的灰度值明顯不同,因此可以進(jìn)一步執(zhí)行閾值法分割出肺部區(qū)域,此時(shí)分割算法會(huì)將肺部區(qū)域作為目標(biāo),將身體區(qū)域作為背景.在最終的三維分割結(jié)果中,本文將背景體素設(shè)置為0,身體區(qū)域體素設(shè)置為1,肺部區(qū)域體素設(shè)置為3.圖1 顯示了第55、65 和75橫斷平面(Axial plane)的分割結(jié)果.圖2(a)顯示4D CT 第1 相位的三維分割結(jié)果,我們可以進(jìn)一步地將肺部區(qū)域和身體區(qū)域分開(kāi),如圖2(b)、圖2(c)所示.

圖1 三維圖像分割結(jié)果橫斷平面展示圖,背景體素值為0,身體區(qū)域的體素值為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 三維圖像分割結(jié)果,及其對(duì)應(yīng)的肺部區(qū)域和身體區(qū)域Fig.2 A 3D segmentation result and the corresponding separated 3D lungs and body

膈肌位于肺部下面,并且與肺部的下邊界同步運(yùn)動(dòng).考慮到難以直接通過(guò)4D CT 觀測(cè)到膈肌結(jié)構(gòu),因此本文使用肺部的下邊界來(lái)代表膈肌.本文工作面向的是肝癌立體定向體部放療(Stereotactic body radiation therapy,SBRT),其主要受到右肺的影響.因此,本文設(shè)計(jì)了右肺自動(dòng)確定方法,來(lái)提取右肺區(qū)域,具體操作如下:首先,在肺部區(qū)域的三維分割結(jié)果上(如圖2(b)),從左到右分別計(jì)算各個(gè)矢狀平面(Sagittal axis)上肺部區(qū)域的面積(體素標(biāo)簽為3 的數(shù)目),這樣以左右方向?yàn)闄M軸,以肺部面積為縱軸,便可得到沿左右方向的肺部面積曲線,如圖3 所示.可以注意到,由于肺部的對(duì)稱性,在兩肺之間肺部面積值比較小,因此肺部面積曲線像字母 “M”一樣具有兩個(gè)峰,每個(gè)峰對(duì)應(yīng)一個(gè)肺.因此,選擇肺部面積曲線兩個(gè)峰之間的最低點(diǎn),便可以確定右肺所包含的切片.

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

接下來(lái),在右肺的三維分割圖像中,可以通過(guò)從下到上(解剖坐標(biāo)系下)確定肺部區(qū)域的最低體素(在分割圖像中,身體體素值為1,肺部區(qū)域體素值等于3),從而獲得右肺下表面的三維圖像.在肺部下表面,有一些高于膈肌的體素,這些體素可以通過(guò)形態(tài)學(xué)運(yùn)算進(jìn)行去除.最后,可以獲得右肺下表面膈肌的三維空間坐標(biāo),表示為:

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

本節(jié)從三維分割結(jié)果的身體區(qū)域(如圖2(c)所示)中獲取胸腹部表面,進(jìn)而計(jì)算胸腹部表面的位移.首先,沿前/后方向找出身體的邊界體素(體素值為1),獲得胸腹部表面.圖4(a)顯示了胸腹部表面在前/后方向的透視圖,其中像素灰度代表距最前端冠狀平面的距離.我們進(jìn)一步對(duì)初始胸腹部表面進(jìn)行處理,首先在得到的初始胸腹部表面圖上(如圖4(a)),沿左/右方向找到最寬的位置,該位置對(duì)應(yīng)人體的肩膀,我們?nèi)コ绨蛑系膮^(qū)域.考慮到人體是一個(gè)類似的圓柱體,其在冠狀面的投影,將會(huì)包含正面的胸腹表面和人體兩側(cè)的區(qū)域,因此,為了更加準(zhǔn)確地獲得胸腹部表面,本文運(yùn)用腐蝕操作,來(lái)進(jìn)一步地去掉人體兩側(cè)的區(qū)域.當(dāng)人體進(jìn)行呼吸時(shí),人體兩側(cè)區(qū)域會(huì)沿左/右方向伸縮,其在冠狀面的投影有顯著的變化,且位于初始胸腹部表面的邊緣區(qū)域.因此,在對(duì)不同相位的初始胸腹部表面進(jìn)行腐蝕后,得到的人體胸腹表面大小不一致,考慮到是同一個(gè)病人,為了統(tǒng)一大小,本文對(duì)所有相位的胸腹表面取交集,我們可以得到如圖4(b)所示的患者胸腹部掩膜.將該掩膜與初始胸腹部表面取交集,我們可以進(jìn)一步獲得最終胸腹部表面,如圖4(c)所示,其表示為:

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

相比于胸腹部的面積,胸腹部位移是非常小的,因此本文采用ICP 算法(Iterative closest point)[31]來(lái)計(jì)算胸腹部表面的位移.ICP 算法的輸出包括一個(gè)三維旋轉(zhuǎn)矩陣和一個(gè)平移向量.通過(guò)測(cè)試發(fā)現(xiàn),人體胸腹部運(yùn)動(dòng)的旋轉(zhuǎn)角度小于0.5 度,其位移主要來(lái)源于平移,因此可以根據(jù)ICP 算法輸出的平移向量,來(lái)計(jì)算胸腹部表面的位移,表示為:

1.3 分步子空間映射

在提取胸腹表面的位移數(shù)據(jù)x ∈X和膈肌的位移數(shù)據(jù)y ∈Y之后,本節(jié)將要在它們之間建立映射F:X →Y,這樣對(duì)于任意新輸入測(cè)試數(shù)據(jù),可以通過(guò)觀測(cè)胸腹表面的位移數(shù)據(jù)xTest,對(duì)膈肌的位移數(shù)據(jù)進(jìn)行預(yù)測(cè)yTest=F(xTest).然而,胸腹表面和膈肌是兩種不同的解剖器官,其運(yùn)動(dòng)模式也不相同,所產(chǎn)生的兩種數(shù)據(jù)具有不同的分布結(jié)構(gòu).例如,在實(shí)驗(yàn)中發(fā)現(xiàn),有的病人膈肌位移很大,而其胸腹表面的運(yùn)動(dòng)卻很小.因此難以在這兩種數(shù)據(jù)集之間直接構(gòu)造映射.為了解決這一問(wèn)題,本節(jié)提出一種新的分布子空間映射算法(Two-step subspace mapping,TSSM).

本文將數(shù)據(jù)分為訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù),x=xTrain∪xTest,y=yTrain∪yTest.假設(shè)器官結(jié)構(gòu)相似的人,其胸腹部表面與膈肌有相似的運(yùn)動(dòng)情況.基于此假設(shè),本文定義兩種距離:1)子空間域內(nèi)距離:dintra(xi,xj)或dintra(yi,yj),i,j=1,2,···,N,計(jì)算在同一個(gè)子空間的兩個(gè)數(shù)據(jù)之間的距離;2)子空間域間距離dinter(xi,yj),i,j=1,2,···,N,計(jì)算來(lái)自不同子空間的兩個(gè)數(shù)據(jù)之間的距離.TSSM 算法的核心是找到一種映射,使同一患者數(shù)據(jù)的子空間域間距離dinter(xi,yi),i=1,2,···,N小于不同患者的子空間域間距離dinter(xi,yj),ij.子空間域內(nèi)距離dintra能夠準(zhǔn)確度量?jī)蓚€(gè)器官運(yùn)動(dòng)的相異性.這樣對(duì)于新輸入的數(shù)據(jù),我們通過(guò)對(duì)觀測(cè)數(shù)據(jù)xTest作映射來(lái)近似表示yTest.下面給出子空間的構(gòu)造方法和映射方法.

膈肌位移數(shù)據(jù)和胸腹部表面位移數(shù)據(jù)的維數(shù)均為K維,因此可以生成K維子空間,選擇一組數(shù)據(jù)作為子空間的基,其他數(shù)據(jù)由這些基來(lái)表示,但這樣會(huì)導(dǎo)致數(shù)據(jù)冗余或不完整.為了克服這個(gè)問(wèn)題,本文采用PCA 來(lái)構(gòu)建子空間[32?34].PCA 通過(guò)最大均值差異最小(Maximum mean discrepancy,MMD)的準(zhǔn)則,找到包含數(shù)據(jù)本質(zhì)特征的公共子空間.

首先,對(duì)膈肌的位移數(shù)據(jù)和胸腹部表面的位移數(shù)據(jù)進(jìn)行歸一化,表示為K維向量的形式,如下所示:

其中,λ ≥0 是控制收縮量的規(guī)則化參數(shù).對(duì)式(6)關(guān)于β求導(dǎo),并使等式等于0,則可以獲得最優(yōu)表達(dá)式,如下所示:

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

與支持向量機(jī)(Support vector machine,SVM)相似,僅使用線性模型進(jìn)行子空間映射,無(wú)法充分地描述數(shù)據(jù)之間的關(guān)系.為了解決這個(gè)問(wèn)題,本文進(jìn)一步將TSSM 推廣到基于核的分步子空間映射算法(kTSSM).令?為非線性空間 Rx到高維空間 Zx的映射,核函數(shù)為K(xi,xj)=〈?(xi),?(xj)〉.可以注意到,預(yù)測(cè)數(shù)據(jù)可通過(guò)變換后數(shù)據(jù)?() 的線性組合來(lái)表示.因此,基于核函數(shù)的非線性嶺回歸優(yōu)化變換可以表示為:

現(xiàn)有關(guān)于體內(nèi)組織器官與體外基準(zhǔn)點(diǎn)相對(duì)運(yùn)動(dòng)的研究[1?2,4,7?10],主要集中在線性模型或二次指數(shù)模型,為了更加詳細(xì)、全面地研究體內(nèi)膈肌和體外胸腹表面的相對(duì)運(yùn)動(dòng)關(guān)系,以進(jìn)行準(zhǔn)確的膈肌位移預(yù)測(cè),本文除了構(gòu)造分步線性映射模型,并進(jìn)一步將其擴(kuò)展到了基于非線性核的映射模型.本文采用了兩種最常見(jiàn)的核函數(shù):多項(xiàng)式核和高斯核,其中多項(xiàng)式核可以實(shí)現(xiàn)將低維的輸入空間映射到高維的特征空間,其并不只限于二次函數(shù),通過(guò)升維操作可以提高建模的準(zhǔn)確性,但多項(xiàng)式核的缺點(diǎn)是參數(shù)較多,容易出現(xiàn)過(guò)擬合現(xiàn)象.高斯核函數(shù)是最常用的核函數(shù)之一,可以將數(shù)據(jù)映射到無(wú)限維空間.考慮到本文醫(yī)用數(shù)據(jù)樣本較少,維數(shù)較高,且沒(méi)有數(shù)據(jù)分布的先驗(yàn)知識(shí),本文進(jìn)一步選擇高斯核對(duì)體內(nèi)膈肌和體外胸腹表面的相對(duì)運(yùn)動(dòng)關(guān)系進(jìn)行度量.在實(shí)驗(yàn)中,本文將對(duì)線性模型和兩種非線性模型進(jìn)行對(duì)比研究,通過(guò)設(shè)置不同的參數(shù),給出膈肌與胸腹部表面之間運(yùn)動(dòng)關(guān)系的準(zhǔn)確度量.

2 實(shí)驗(yàn)

2.1 患者情況與數(shù)據(jù)采集

本研究中使用的臨床數(shù)據(jù)由Philips Brilliance CT Big Bore Scanner 設(shè)備獲得,在這項(xiàng)研究中采用了來(lái)自20 位患者的4D CT 數(shù)據(jù)集,每個(gè)患者包括10個(gè)相位(具體為平均呼吸周期的0%,10%,···,90%).每個(gè)胸部CT 的大小為512×512×185個(gè)體素(Voxel),尺寸分別為1.1719 mm×1.1719 mm×3 mm (分別對(duì)應(yīng)于A/P、I/S和L/R三個(gè)方向).在采集數(shù)據(jù)時(shí),病人按正常狀態(tài)、有規(guī)律的呼吸.

在實(shí)驗(yàn)中,我們選擇相關(guān)研究中的經(jīng)典線性模型和非線性模型,與本文所提三種模型進(jìn)行對(duì)比,具體如下:OLS (Ordinary least squares)多變量線性回歸模型[1]、人工神經(jīng)網(wǎng)絡(luò)模型(Artificial neural network,ANN)、TSSM、多項(xiàng)式核的kTSSM 和高斯核的kTSSM.本文采用三個(gè)指標(biāo),對(duì)5 種算法的性能進(jìn)行定量評(píng)估:均方誤差(Mean-square error,MSE)、R2誤差和平均絕對(duì)百分比誤差(Mean absolute percentage error,MAPE),計(jì)算方法如下:

其中,yi是模型的實(shí)際值,是模型實(shí)際值的均值,是模型的預(yù)測(cè)值.MSE 衡量模型的預(yù)測(cè)均方誤差,MSE 的值越低,預(yù)測(cè)結(jié)果越好.R2比較模型的優(yōu)劣性,消除數(shù)據(jù)分布域?qū)?shí)驗(yàn)結(jié)果的影響.R2的范圍為(?∞,1],值越大預(yù)測(cè)效果越好.MAPE 衡量模型的預(yù)測(cè)值與實(shí)際值的平均偏差,消除數(shù)據(jù)取值范圍的影響,值越低預(yù)測(cè)結(jié)果越好.

在實(shí)驗(yàn)中,本文采用K-fold 策略將整個(gè)數(shù)據(jù)集分為訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù),其中80%的隨機(jī)樣本作為訓(xùn)練數(shù)據(jù),其余樣本進(jìn)行模型測(cè)試.為了評(píng)估該算法的魯棒性,在訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù)中分別添加高斯噪聲σ={0.1,0.2,0.3,0.4,0.5},從而分析不同級(jí)別的高斯噪聲對(duì)算法的影響.本文算法參數(shù)選擇如下:對(duì)于嶺回歸優(yōu)化,λ={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};在多項(xiàng)式核K(x,y)=(xTy+c)d中,d={0.5,1,···,11},c={?10,?9,···,10};對(duì)于ANN,本文采用兩層前饋神經(jīng)網(wǎng)絡(luò),其隱層節(jié)點(diǎn)數(shù)依次選擇{2,3,4,···,29,30}.對(duì)于每個(gè)算法的每組參數(shù),本文單獨(dú)運(yùn)行算法100 次,選擇每個(gè)算法的最優(yōu)參數(shù)組合,并獲得統(tǒng)計(jì)結(jié)果.所有算法的運(yùn)行環(huán)境相同,均為Inter Xeon 3.6 GHz,32 GB 內(nèi)存,Windows10 和MATLAB R2019b 平臺(tái).

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

圖6 膈肌和胸腹部表面的位移,每一種顏色的曲線對(duì)應(yīng)一個(gè)患者的數(shù)據(jù)Fig.6 The displacement of diaphragm and thoracoabdominal surface,where each color corresponds to a specific patient＇s data

由于性別、年齡、體型和情緒等都會(huì)對(duì)身體器官的運(yùn)動(dòng)產(chǎn)生影響,不同患者器官運(yùn)動(dòng)的幅度不盡相同.因此,根據(jù)式(5)對(duì)膈肌和胸腹部表面的位移進(jìn)行歸一化,結(jié)果如圖7 所示.從圖7 可以看出,不同患者胸腹部表面在前/后方向,膈肌在下/上運(yùn)動(dòng)方向的數(shù)據(jù)分布更加一致.

圖7 歸一化后的膈肌和胸腹部表面的位移,每一種顏色的曲線對(duì)應(yīng)一個(gè)患者的數(shù)據(jù)Fig.7 The displacement of diaphragm and thoracoabdominal surface after normalization,where each color corresponds to a specific patient＇s data

2.3 三個(gè)方向的預(yù)測(cè)

在本節(jié)中,將根據(jù)胸腹部表面在三個(gè)方向的運(yùn)動(dòng),來(lái)預(yù)測(cè)膈肌在三個(gè)方向上的運(yùn)動(dòng).本文將同一患者三個(gè)方向10個(gè)相位的數(shù)據(jù)連接在一起,形成一個(gè)包含30個(gè)元素的向量,胸腹部表面位移向量為輸入,膈肌位移向量為輸出,預(yù)測(cè)結(jié)果如表1 所示.其中每個(gè)單元有兩個(gè)數(shù)據(jù),前一個(gè)值代表100 次獨(dú)立運(yùn)行的平均值,括號(hào)中的值代表100 次獨(dú)立運(yùn)行的標(biāo)準(zhǔn)偏差.本實(shí)驗(yàn)記錄MSE 和MAPE 的最小值,記錄R2的最大值,加粗標(biāo)記代表5 種算法中的最佳結(jié)果.在表1 中,從上到下,噪聲等級(jí)依次增加.從表1 可以看出,同為線性模型,TSSM 算法優(yōu)于OLS 模型;非線性模型ANN 的表現(xiàn)較OLS 模型要更好.所提的三個(gè)算法的預(yù)測(cè)結(jié)果相近,MSE值和R2值均令人滿意,但多項(xiàng)式核的kTSSM 算法獲得多數(shù)的最佳預(yù)測(cè)值.值得一提的是,由于向量取值的不平衡性,兩種組織器官在有的相位和方向上移動(dòng)范圍大,在其它相位和方向上移動(dòng)范圍小,取值范圍小的預(yù)測(cè)值會(huì)增加MAPE 值,因此使得MAPE 值均超過(guò)100%.

表1 在三個(gè)方向上的預(yù)測(cè)結(jié)果Table 1 Prediction results in three directions

圖8 顯示5 種算法預(yù)測(cè)結(jié)果,其對(duì)應(yīng)100 次獨(dú)立運(yùn)行中,MSE 的中間值對(duì)應(yīng)的結(jié)果.圖8(a)為 OLS模型,圖8(b)為ANN,圖8(c)為TSSM,圖8(d)為多項(xiàng)式內(nèi)核的kTSSM,圖8(e) 為高斯內(nèi)核的kTSSM.在每個(gè)子圖中,有4個(gè)預(yù)測(cè)結(jié)果(占總數(shù)的20%),其中橫軸表示10個(gè)相位,縱軸表示實(shí)際位移值,紅色曲線表示實(shí)際值,藍(lán)色曲線表示預(yù)測(cè)值.在橫軸上,L/R表示左/右方向,I/S表示下/上方向,而A/P表示前/后方向.從圖8(a)可以看出,OLS 模型在下/上方向上可以預(yù)測(cè)膈肌運(yùn)動(dòng),但是在某些相位,OLS 模型在左/右和前/后兩個(gè)方向上的預(yù)測(cè)效果不佳,仍然存在誤差.同為線性模型,TSSM 算法具有更好的預(yù)測(cè)效果.與OLS 模型相比,非線性模型ANN 的表現(xiàn)更好,而本文設(shè)計(jì)的kTSSM (多項(xiàng)式)在大多數(shù)情況下,取得了最佳的預(yù)測(cè)結(jié)果.

圖8 三個(gè)方向的預(yù)測(cè)結(jié)果.該結(jié)果對(duì)應(yīng)于100 次運(yùn)行MSE 的中位數(shù)值Fig.8 Prediction results corresponding the median MSE value in three directions after 100 independent runs

圖9 顯示了100 次獨(dú)立運(yùn)行的統(tǒng)計(jì)箱圖,其中橫軸表示三個(gè)方向10個(gè)相位,縱軸表示預(yù)測(cè)值和實(shí)際值之間的絕對(duì)誤差.圖9 子圖中的每個(gè)方框都包含100 次獨(dú)立運(yùn)行的結(jié)果,在每個(gè)方框上,中心的標(biāo)記表示中值,而框的底部和頂部邊緣分別表示第25個(gè)和第75個(gè)百分點(diǎn),“+”符號(hào)代表異常值.從圖9我們可以清楚地看到,圖9(a)中縱軸的范圍比其他子圖中的大得多,圖9(b)中縱軸取值范圍居中,而圖9(c)～圖9(e)的誤差絕對(duì)值比較小,這也說(shuō)明本文所提的方法具有更優(yōu)、更穩(wěn)定的預(yù)測(cè)性能.從圖9(c)～圖9(e)可以看到,下/上方向上的統(tǒng)計(jì)值低于其他方向上的統(tǒng)計(jì)值;對(duì)于下/上方向,位于中間相位(5、6、7 相位)的預(yù)測(cè)誤差較小,這表明當(dāng)患者處于滿吸氣狀態(tài)時(shí),預(yù)測(cè)效果會(huì)更好.

圖9 三個(gè)方向各相位預(yù)測(cè)誤差的統(tǒng)計(jì)箱圖Fig.9 Boxplot of the prediction performance on every phase in three directions

2.4 主方向的預(yù)測(cè)

在本節(jié)中,僅使用胸腹部表面的前/后方向,對(duì)膈肌的下/上方向的運(yùn)動(dòng)進(jìn)行預(yù)測(cè).表2 顯示了100 次獨(dú)立運(yùn)行后的預(yù)測(cè)結(jié)果.圖10 顯示5 種方法對(duì)應(yīng)MSE 中間值的預(yù)測(cè)結(jié)果.圖11 為每個(gè)相位的預(yù)測(cè)誤差絕對(duì)值的統(tǒng)計(jì)箱圖.從表2 可以看出,TSSM 算法的性能優(yōu)于OLS 模型.ANN 算法結(jié)果優(yōu)于OLS 模型.在沒(méi)有左/右和前/后方向的數(shù)據(jù)影響的情況下,TSSM 算法的性能得到了很大的提高,而且MAPE 值已降至100%以下,證明了TSSM算法的良好性能.從圖10 可以看出,本文所提方法均取得了令人滿意的預(yù)測(cè)結(jié)果.圖11 可以看出,在第5、6、7 相位上的預(yù)測(cè)結(jié)果,要好于在其它相位上的預(yù)測(cè)結(jié)果,這也進(jìn)一步證明,當(dāng)患者處于滿吸氣狀態(tài)時(shí),從外部體表能更加準(zhǔn)確地預(yù)測(cè)體內(nèi)結(jié)構(gòu)的運(yùn)動(dòng).

圖10 主方向的預(yù)測(cè)結(jié)果,該結(jié)果對(duì)應(yīng)于100 次運(yùn)行MSE 的中位數(shù)值Fig.10 Prediction results corresponding the median MSE values in the principal direction after 100 independent runs

圖11 主方向各相位預(yù)測(cè)誤差的統(tǒng)計(jì)箱圖Fig.11 Boxplot of the prediction performance on every phase in the principal direction

表2 在主方向上的預(yù)測(cè)結(jié)果Table 2 Prediction results in principal direction

2.5 噪聲的影響

本節(jié)驗(yàn)證算法的魯棒性.噪聲對(duì)三個(gè)方向預(yù)測(cè)性能的影響如圖12 所示.在圖12 的子圖中,橫軸表示高斯噪聲的等級(jí),縱軸表示算法100 次獨(dú)立運(yùn)行后指標(biāo)的平均值.圖13 是噪聲對(duì)主方向預(yù)測(cè)的影響,可以看出,隨著噪聲等級(jí)的增加5 種算法的性能都會(huì)下降.但是與OLS 模型相比,ANN 算法和TSSM 算法的下降速度較慢且穩(wěn)定,這表明ANN,TSSM 和kTSSM 算法比OLS 模型對(duì)噪聲的魯棒性更高;對(duì)于準(zhǔn)確性和穩(wěn)定性,尤其是在MAPE度量方面,具有多項(xiàng)式核的kTSSM 算法在5 種算法中表現(xiàn)最好.

圖12 噪聲對(duì)三個(gè)方向預(yù)測(cè)性能的影響Fig.12 Influence of noise on the prediction in three directions

圖13 噪聲對(duì)主方向預(yù)測(cè)性能的影響Fig.13 Influence of noise on the prediction in principal direction

2.6 關(guān)于算法參數(shù)的討論

TSSM 和kTSSM 取得最優(yōu)值時(shí),參數(shù)設(shè)置如表3 所示.對(duì)于線性模型TSSM,其參數(shù)只有嶺回歸一個(gè)參數(shù)λ,其取值范圍為λ={10?10,10?9,···,1010}.在絕大多數(shù)情況下,λ的最優(yōu)設(shè)置為1,只有兩次為0.1,這兩次對(duì)應(yīng)λ=1 的預(yù)測(cè)結(jié)果與最優(yōu)結(jié)果接近.因此,TSSM 對(duì)膈肌位移的預(yù)測(cè)是比較穩(wěn)定的,受參數(shù)影響較小,可直接設(shè)置為1 即可.對(duì)于kTSSM(高斯)模型,其參數(shù)為λ和σkernel,其中λ的取值范圍與TSSM 相同,σkernel={0.1,0.2,0.4,0.8,1.6,3.2,6.4,12.8}. 可以看到,在絕大多數(shù)情況下λ=0.1,σkernel=0.1 或0.2,只相差一個(gè)取值間隔.在實(shí)驗(yàn)結(jié)果中我們發(fā)現(xiàn),當(dāng)σkernel=0.1 或σkernel=0.2 時(shí),kTSSM (高斯)模型的預(yù)測(cè)結(jié)果非常相近,具有很好的穩(wěn)定性.kTSSM (高斯)模型的參數(shù)設(shè)置,可按如下策略:當(dāng)獲得的數(shù)據(jù)沒(méi)有噪聲或者比較干凈時(shí),可直接設(shè)置為λ=0.1,σkernel=0.2;當(dāng)數(shù)據(jù)噪聲比較 大 時(shí),參 數(shù) 可 直 接 設(shè) 置 為λ=0.1,σkernel=0.1.對(duì)于kTSSM (多項(xiàng)式)模型,其具有三個(gè)參數(shù)λ、d和c,其中λ的取值范圍與TSSM 相同,d={0.5,1,···,11},c={?10,?9,···,10}.可以看到,相比于TSSM 和kTSSM (高斯)模型,kTSSM (多項(xiàng)式)模型的參數(shù)變化浮動(dòng)較大,特別是同時(shí)進(jìn)行膈肌三個(gè)位移方向的預(yù)測(cè)時(shí),嶺回歸參數(shù)λ、指數(shù)d隨著噪聲的增加逐漸降低,最后維持在λ=10 和d=1.這主要是因?yàn)樵跊](méi)有噪聲或噪聲較小時(shí),kTSSM(多項(xiàng)式)模型發(fā)生了過(guò)擬合現(xiàn)象;當(dāng)噪聲增加時(shí),降低了kTSSM (多項(xiàng)式)模型對(duì)數(shù)據(jù)的敏感性,因而最優(yōu)參數(shù)可以達(dá)到相對(duì)穩(wěn)定.kTSSM (多項(xiàng)式)模型對(duì)于主方向的預(yù)測(cè),受參數(shù)影響較小,比較穩(wěn)定,λ在1 上下浮動(dòng),d在0.5 和1 之間選擇,參數(shù)c設(shè)置為4 即可.從kTSSM (多項(xiàng)式)模型的參數(shù)d可以看出,胸腹部表面與膈肌的相對(duì)位移,近似滿足線性模型(d=1),這一點(diǎn)也可以通過(guò)TSSM 的參數(shù)穩(wěn)定性得到驗(yàn)證.盡管kTSSM (多項(xiàng)式)模型在實(shí)驗(yàn)中取得了最好的預(yù)測(cè)結(jié)果,但是考慮到kTSSM(多項(xiàng)式)模型具有更多的參數(shù)需要設(shè)置,并容易出現(xiàn)過(guò)擬合現(xiàn)象,本文認(rèn)為在實(shí)際中,只需要采用線性模型TSSM 即可達(dá)到準(zhǔn)確的預(yù)測(cè)結(jié)果.

表3 TSSM 和kTSSM 取得最優(yōu)預(yù)測(cè)結(jié)果時(shí)的參數(shù)設(shè)置Table 3 Parameter setting of TSSM and kTSSM for optimal results

本文進(jìn)一步使用TSSM 模型研究PCA 算法中不同的閾值ε對(duì)算法性能的影響.我們選擇不同的閾值ε={1,10?1,10?2,10?3,10?4,10?5,10?6},確保由PCA 構(gòu)造子空間的維數(shù),從1 維逐漸變到原始數(shù)據(jù)的維數(shù),然后在不同的噪聲水平(σ={0.1,0.2,0.3,0.4,0.5}) 下,單獨(dú)運(yùn)行TSSM 模型(λ=1) 100 次.計(jì)算結(jié)果如圖14 所示,其中圖14(a)為三個(gè)方向預(yù)測(cè)的統(tǒng)計(jì)結(jié)果,圖14(b)為主方向預(yù)測(cè)的統(tǒng)計(jì)結(jié)果.每個(gè)子圖中,橫坐標(biāo)表示不同的閾值ε,縱坐標(biāo)依次為子空間維數(shù)、MSE 值、R2值和MAPE值.可以看到,當(dāng)ε=1 時(shí),此時(shí)PCA 構(gòu)造的子空間只有1 維;隨著ε的增大,子空間維數(shù)逐漸增大至原始數(shù)據(jù)的維數(shù),即PCA 沒(méi)有對(duì)原始數(shù)據(jù)進(jìn)行降維.隨著ε的增大,MSE 值逐漸減小,這主要是由于子空間映射時(shí),維數(shù)越多,信息量損失越小,因此可以得到更加準(zhǔn)確的預(yù)測(cè)結(jié)果.同時(shí)注意到,MSE值、R2值和MAPE值的變化范圍非常小,特別是R2值和MAPE 值僅出現(xiàn)小范圍的變化,可近似看作穩(wěn)定.即使當(dāng)PCA構(gòu)造的子空間只有1 維時(shí),TSSM模型也能取得令人滿意的預(yù)測(cè)結(jié)果,這也說(shuō)明PCA構(gòu)造的子空間能夠有效提取原始數(shù)據(jù)的主要成分,具有很強(qiáng)的穩(wěn)定性.為了對(duì)比驗(yàn)證所提算法的有效性,本文在實(shí)驗(yàn)中選擇ε=10?4,此時(shí)PCA 構(gòu)造的子空間保留了原始數(shù)據(jù)絕大部分信息,信息損失量較少,獲得預(yù)測(cè) 結(jié)果較為準(zhǔn)確.

圖14 PCA 中的閾值 ε 對(duì)TSSM 預(yù)測(cè)模型的影響Fig.14 The influence of ε in PCA on the prediction performance of TSSM

3 討論和結(jié)論

呼吸運(yùn)動(dòng)是放射治療(RT)中導(dǎo)致誤差和不確定性的重要因素[1?2].本文首先采用三維圖像分割技術(shù),對(duì)4D CT 圖像進(jìn)行分割,在不使用標(biāo)記物的情況下,通過(guò)提取膈肌的質(zhì)心以準(zhǔn)確度量膈肌的位移,使用ICP 算法計(jì)算胸腹部表面的位移.為了解決跨空間預(yù)測(cè)問(wèn)題,本文提出了一種新的分步子空間映射算法(TSSM),通過(guò)構(gòu)造特征子空間,并在高維子空間中進(jìn)行映射,從而可以有效地提高跨空間數(shù)據(jù)預(yù)測(cè)的準(zhǔn)確性.

對(duì)20 位患者真實(shí)數(shù)據(jù)的實(shí)驗(yàn)表明,本文提出的方法可以有效地描述膈肌和胸腹部表面的運(yùn)動(dòng)情況.相比于通過(guò)原始空間的直接預(yù)測(cè),TSSM 算法在準(zhǔn)確性和魯棒性方面均優(yōu)于經(jīng)典線性的OLS 模型和非線性的ANN 模型.考慮到具有多項(xiàng)式核的kTSSM 算法比TSSM 算法具有更多的參數(shù),容易出現(xiàn)過(guò)擬合現(xiàn)象,因此在實(shí)際應(yīng)用中,我們建議直接應(yīng)用TSSM 算法,線性模型已經(jīng)可以獲得令人滿意的預(yù)測(cè)效果.值得一提的是,我們發(fā)現(xiàn)在病人吸氣處于最高點(diǎn)的位置時(shí),預(yù)測(cè)效果會(huì)更好.所提算法具有解析解,總計(jì)算時(shí)間(訓(xùn)練加上測(cè)試)在0.6 ms以下,可以達(dá)到實(shí)時(shí)處理速度,運(yùn)行效率較高,這將有助于提高放射治療中門控技術(shù)和跟蹤技術(shù)的效率和精度.