頭頸癌放療計劃劑量分布的預(yù)測方法：基于深度學(xué)習(xí)的算法

滕 琳，王 斌，馮前進

1南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，2廣東省醫(yī)學(xué)圖像處理重點實驗室，廣東 廣州 510515；3上?？萍即髮W(xué)生物醫(yī)學(xué)工程學(xué)院，上海 201220

頭頸癌是一類癌癥的統(tǒng)稱，包括口腔癌、咽癌和鼻癌等眾多類型［1］，其發(fā)病率在國內(nèi)位居第九，對人類健康產(chǎn)生極大的威脅［2］。由于頭頸部危機器官眾多，解剖關(guān)系復(fù)雜，為了殺死腫瘤細(xì)胞，同時保護周圍正常組織，調(diào)強放射治療（IMRT）成為首選治療方式［3］。為了使計劃靶區(qū)（PTV）達到高劑量輻射，危機器官（OARs）接收低劑量輻射甚至零劑量輻射，放療物理師需要根據(jù)自身經(jīng)驗反復(fù)調(diào)整治療參數(shù)，從而獲得最優(yōu)的放療計劃。此過程需要一周左右的時間，從而延誤了病情的最佳治療時間［4-7］。

近幾年，隨著機器學(xué)習(xí)的迅速發(fā)展［8］，特別是深度學(xué)習(xí)，涌現(xiàn)出很多基于深度學(xué)習(xí)的方法，實現(xiàn)放療計劃自動設(shè)計的目標(biāo)［9-12］。即利用CT等圖像作為輸入信息，通過卷積神經(jīng)網(wǎng)絡(luò)CNN，根據(jù)目標(biāo)任務(wù)自動提取有用特征，通過損失函數(shù)監(jiān)督網(wǎng)絡(luò)學(xué)習(xí)，網(wǎng)絡(luò)輸出信息為劑量分布圖。舉例來說，Liu等［13］設(shè)計了一個級聯(lián)網(wǎng)絡(luò)架構(gòu)，將第一網(wǎng)絡(luò)預(yù)測的不同尺度特征圖作為第二網(wǎng)絡(luò)的輸入，可以充分提取并整合全局和局部特征。然而網(wǎng)絡(luò)關(guān)注的是全局空間學(xué)習(xí)，導(dǎo)致ROIs等局部區(qū)域的劑量預(yù)測不準(zhǔn)確。為了解決此問題，Ngyuen等［14］提出可微分的DVH損失函數(shù)，將DVH劑量學(xué)參數(shù)作為網(wǎng)絡(luò)優(yōu)化目標(biāo)，監(jiān)督網(wǎng)絡(luò)在ROIs等局部區(qū)域的特征提取。然而DVH劑量學(xué)參數(shù)的計算要求重復(fù)遍歷整幅圖像，導(dǎo)致劑量預(yù)測的效率低下。為了進一步提升劑量預(yù)測精確性，Zhang 等［15］提出距離場的概念，即計算所有體素到PTV邊界的距離，將PTV的相對距離與劑量分布的關(guān)系等先驗信息作為網(wǎng)絡(luò)輸入。然而，上述方法的網(wǎng)絡(luò)輸入信息（CT圖像和ROIs分割掩膜）缺失高能輻射是通過射線束形式傳送等先驗信息，導(dǎo)致射線束等局部區(qū)域預(yù)測不準(zhǔn)確。

為了解決上述方法存在的問題，本研究提出新穎的射線束劑量肢解學(xué)習(xí)方法，引入射線束分割掩膜和ROIs內(nèi)劑量分布特征作為先驗信息訓(xùn)練CNN，進一步提升局部區(qū)域的劑量預(yù)測精確性。具體來說，本文方法分為兩個階段，第一階段，利用全局劑量網(wǎng)絡(luò)（GDN）得到全局空間的粗糙劑量分布圖。同時，利用預(yù)定義的射線束角度和PTV分割掩膜，擬合射線束分割掩膜；第二階段，進一步精細(xì)化第一階段預(yù)測的粗糙劑量分布圖，通過以下四種新穎方法：（1）通過射線束分割掩膜，將粗糙劑量分布圖肢解為多個子劑量分布圖，每個子劑量分布圖代表某個射線束方向的劑量分布；（2）通過射線束投票機制，將多個子劑量分布圖聚合為一個劑量分布圖；（3）通過DVH校正，即基于值的和基于評價指標(biāo)的DVH損失函數(shù)，將ROIs內(nèi)劑量分布特征作為損失函數(shù)，約束網(wǎng)絡(luò)在ROIs內(nèi)的特征提??；（4）引入邊界增強約束，計算預(yù)測和真實的三維梯度圖之間的MAE，監(jiān)督網(wǎng)絡(luò)在劑量邊界區(qū)域的特征提取。

1 材料和方法

1.1 實驗數(shù)據(jù)

利用來自O(shè)penKBP-2020 挑戰(zhàn)公開的頭頸癌放療計劃數(shù)據(jù)集［17］，包括340例患者。其中，200例患者作為訓(xùn)練集，40例作為驗證集，100例作為測試集。每個患者包括CT圖像，PTV分割掩膜，OAR分割掩膜，有效劑量值區(qū)域分割掩膜和真實劑量分布圖。劑量分布圖是通過9例等間距共面射線束（0°，40°，…，320°）傳送高能輻射形成的，每個射線束的能量是6 MV，使用步進式靜態(tài)調(diào)強放療技術(shù)IMRT。所有數(shù)據(jù)大小為128×128×128 mm3，體素近似為3.906 mm×3.906 mm×2.5 mm。

1.2 方法

本研究提出的方法分為兩個階段（圖1），第一個階段為全局空間劑量學(xué)習(xí)，預(yù)測粗糙劑量分布圖；第二階段為劑量精細(xì)化階段，精細(xì)化局部區(qū)域的劑量分布，即射線束方向，ROIs以及劑量邊界區(qū)域。方法的實施細(xì)節(jié)見1.3節(jié)。

圖1 方法流程Fig.1 Framework of the proposed method.A: Global-wise dose learning.B: Beam-wise dose refinement.C: Global-wise dose refinement.

1.2.1 全局劑量學(xué)習(xí) 為了得到全局空間粗糙劑量分布圖，將CT圖像，OAR分割掩膜，PTV分割掩膜拼接為多個通道信息，輸入到GDN（網(wǎng)絡(luò)細(xì)節(jié)見1.3節(jié)），預(yù)測全局空間粗糙劑量分布圖，訓(xùn)練過程使用平均絕對誤差（MAE）損失函數(shù)進行監(jiān)督。然而，由于未考慮射線束與ROIs之間的空間幾何關(guān)系，導(dǎo)致劑量分布圖在局部區(qū)域預(yù)測不準(zhǔn)確，例如射線束區(qū)域和ROIs。為了解決此問題，本文提出劑量精細(xì)化方法，包括基于射線束的劑量精細(xì)化和全局劑量精細(xì)化。

1.2.2 基于射線束的劑量精細(xì)化 利用擬合的射線束分割掩膜，精細(xì)化射線束方向上的劑量分布。具體步驟如下：

1.2.2.1 射線束分割掩膜的形成 IMRT是體外放射治療的常用技術(shù)之一，頭頸癌放療一般要求7～9個不同角度共面射線束對PTV進行照射。因此，根據(jù)預(yù)定義的射線束角度和PTV的邊界信息，擬合放療中的多個角度共面射線束，稱為“射線束分割掩膜”。由于數(shù)據(jù)集是三維的，不同切片上的腫瘤位置和大小不同。為了提升擬合射線束的準(zhǔn)確性，首先我們擬合二維射線束分割掩膜，然后將其堆疊為三維射線束分割掩膜。具體來說，以二維切片的形式從PTV中提取腫瘤邊界，得到二維邊界切片（Ei），根據(jù)多個預(yù)定義角度（θi），在二維邊界切片Ei上按照各預(yù)設(shè)角度θi勾勒兩條邊界切線，兩條切線之間的區(qū)域為該預(yù)設(shè)角度的射線束分割掩膜；將擬合的全部二維射線束分割掩膜進行堆疊，得到三維射線束分割掩膜（圖2）。

圖2 不同方法預(yù)測的劑量分布圖的可視化結(jié)果Fig.2 Visualization of the dose distribution map predicted using different state-of-the-art methods.The red boxes indicate the obviously improved regions.The color closer to blue indicates a lower dose value;a color closer to red indicates a higher dose value.

1.2.2.2 基于射線束的劑量學(xué)習(xí) 在實際放療中，為了使PTV內(nèi)劑量分布均勻，不同角度的射線束傳送不同的劑量分布，導(dǎo)致CNN同時學(xué)習(xí)不同角度射線束內(nèi)的劑量分布是困難的。為了解決此問題，本文利用擬合的射線束分割掩膜，將粗糙劑量分布圖肢解為多個沿著射線束方向的子劑量分布圖。因此，全局劑量預(yù)測任務(wù)轉(zhuǎn)變?yōu)橐幌盗械木植孔觿┝款A(yù)測任務(wù)，使網(wǎng)絡(luò)更容易提取每個不同射線束路徑下的劑量分布特征。將上述GDN預(yù)測的粗糙劑量分布圖，CT圖像，OARs分割掩膜，PTV分割掩膜和擬合的射線束分割掩膜拼接成多通道信息，輸入射線束劑量網(wǎng)絡(luò)（BDN），預(yù)測多個和輸入數(shù)據(jù)相同大小的子劑量分布圖，每個子劑量分布圖代表某射線束方向的劑量分布，訓(xùn)練過程使用MAE進行監(jiān)督，其定義為公式（1）：

其中，N代表子劑量分布圖的個數(shù)，即N=9，P代表預(yù)測劑量分布圖，G為真實劑量分布圖。

1.2.2.3 聚合多個子劑量分布 為了將多個局部空間的子劑量分布圖聚合為一個全局空間的劑量分布圖，本文提出多射線束投票機制，即聚合后的劑量分布圖中的劑量值由多個子劑量分布進行投票決定。具體來說，若某一體素由多個子劑量分布經(jīng)過，則此體素的劑量值由經(jīng)過此體素的多個子劑量分布通過平均投票決定。若某一體素由單個子劑量分布經(jīng)過，則此體素的劑量值由該單個子劑量分布中相應(yīng)體素位置的劑量值決定。最后，通過MAE損失函數(shù)，監(jiān)督聚合后的全局劑量分布圖學(xué)習(xí)，其定義為公式（2）：

其中，N為體素的個數(shù)，為預(yù)測劑量分布圖，Y為真實劑量分布圖。

1.2.3 全局劑量精細(xì)化 為了進一步精細(xì)化ROIs內(nèi)劑量分布，本文提出新穎的DVH校正和邊界增強約束方法。

1.2.3.1 DVH 校正 在臨床放療中，放療物理師通過DVH評價放療計劃質(zhì)量，其橫坐標(biāo)代表劑量值，縱坐標(biāo)代表接收相應(yīng)劑量值照射的相對體積［18］。為了利用該信息，Nguyen 等［14］提出了可微分的DVH損失函數(shù)，計算預(yù)測和真實DVH劑量學(xué)參數(shù)的體積之差（即DVH縱坐標(biāo)）。然而，通過選取劑量閾值，遍歷整幅圖像計算DVH劑量學(xué)參數(shù)，導(dǎo)致計算精度和效率受所取的劑量閾值的影響。舉例來說，選取的劑量閾值間隔小，計算精度提升，計算效率下降；選取的劑量閾值間隔大，計算效率提升，計算精度下降。為了解決此問題，本研究提出了基于值的DVH損失函數(shù)，其定義為公式（3）：

其中，N為ROIs的個數(shù)，R為排序操作，G為真實劑量分布圖，P為預(yù)測劑量分布圖，Ws為第s 個ROI分割掩膜。具體來說，首先利用ROIs分割掩膜提取ROIs內(nèi)的劑量值，然后通過排序操作對其劑量值進行排序，最后計算排序后預(yù)測和真實結(jié)果的劑量值之差。值得注意的是，DVH體積信息（縱坐標(biāo)）的計算原理與劑量值排序的操作原理十分相似。本研究提出的方法在利用DVH先驗信息時，只需遍歷一次圖像的計算，沒有涉及超參數(shù)的選擇，因此更有效和更精確。

除此之外，放療物理師更關(guān)注DVH上的關(guān)鍵點評價放療計劃的可行性，例如分別指1%體積的PTV所接受的劑量，95%體積的PTV所接受的劑量，99%體積的PTV所接受的劑量，0.1cc體積的OAR所接受的最大劑量，OAR區(qū)域內(nèi)的平均劑量。為了有效利用這些關(guān)鍵點，本研究提出基于評價指標(biāo)的DVH損失函數(shù)，進一步提升ROIs內(nèi)的預(yù)測精確性，其定義為公式（4）：

1.2.3.2 邊界增強約束 劑量分布圖存在銳利的邊界，主要原因包括：（1）高能輻射通過射線束形式進行傳送，射線束內(nèi)劑量值高，射線束外劑量值低，劑量值變化大區(qū)域形成銳利邊界；（2）為了殺死腫瘤細(xì)胞，同時保護周圍正常組織，PTV區(qū)域內(nèi)劑量值高，其周圍區(qū)域劑量值低，劑量值呈現(xiàn)PTV到周圍區(qū)域迅速下降的趨勢。然而，由于提取特征的感受野較大，使卷積神經(jīng)網(wǎng)絡(luò)提取銳利的邊界特征是困難的，導(dǎo)致預(yù)測結(jié)果在邊界區(qū)域過于平滑。受啟發(fā)于Tan等［19］提出的利用多任務(wù)學(xué)習(xí)加強銳利邊界的學(xué)習(xí)，預(yù)測劑量分布作為主要任務(wù)，預(yù)測邊界圖作為輔助任務(wù)。由于提取特征的感受野較大，銳利的邊界特征較難提取，導(dǎo)致CNN不能精確預(yù)測邊界圖。另外，引入一個單獨的解碼器進行邊界圖的預(yù)測，增加計算參數(shù)量。因此，本研究引入三維梯度損失函數(shù)，約束CNN 更加關(guān)注劑量分布圖中邊界區(qū)域的學(xué)習(xí)。具體來說，應(yīng)用Sobel算子，提取預(yù)測劑量分布圖和真實劑量分布圖的梯度圖，其代表著臨近體素的劑量值之差。計算預(yù)測和真實的梯度圖之間的MAE，其定義為公式（5）：

其中，Y為真實劑量分布圖，為預(yù)測劑量分布圖，Sx,Sy,Sz分別代表x,y,z維度上的Sobel操作。

總之，在BDN的訓(xùn)練期間，總損失函數(shù)為：

其中，α1，α2，α3為超參數(shù)，用于平衡不同損失函數(shù)的貢獻。

1.3 實驗

所有實驗以Pytorch框架搭建，NVIDIA TITANV GPU 12GB條件下運行。批量大小為2，使用Adam優(yōu)化器［20］，初始學(xué)習(xí)率為10-4。為了避免過擬合現(xiàn)象發(fā)生，我們使用在線數(shù)據(jù)增廣策略，包括隨機旋轉(zhuǎn)（-20°～20°）和沿著Z 軸隨機翻轉(zhuǎn)。CT 圖像強度值歸一化到0～1。GDN和BDN都是基于3D U-Net［21］網(wǎng)絡(luò)結(jié)構(gòu)。

為了公平比較，利用OpenKBP-2020挑戰(zhàn)［17］提供的量化指標(biāo)Dose score和DVH score。Dose score 計算預(yù)測劑量分布圖和真實劑量分布圖之間的體素級MAE，DVH score計算預(yù)測結(jié)果和真實結(jié)果之間的DVH評價指標(biāo)MAE。Doses score和DVH score 越低，說明模型的預(yù)測性能越好。另外，我們也應(yīng)用DVH 曲線驗證ROIs內(nèi)劑量預(yù)測準(zhǔn)確性。

為了驗證本文方法優(yōu)越性，與目前該領(lǐng)域4 種SOTA方法（V-Net［22］，DCNN［15］，HD U-net［23］，C3D［13］）進行對比。C3D在OpenKBP-2020挑戰(zhàn)中排名第一。利用公開代碼復(fù)現(xiàn)后3種方法，復(fù)現(xiàn)V-Net方法使用與本文方法相同的訓(xùn)練參數(shù)。

為了驗證本文方法有效性，對各組成部分進行消融實驗分析。在級聯(lián)網(wǎng)絡(luò)的基礎(chǔ)上，引入本文提出的新穎方法，包括3部分：（1）基于射線束的劑量預(yù)測（BDP）（2）DVH 校正：基于值的DVH 損失函數(shù)（LvDVH）和基于評價指標(biāo)的DVH損失函數(shù)（LcDVH），（3）邊界增強約束（Le）。

2 結(jié)果

2.1 對比實驗

針對本文方法的對比實驗分析，量化結(jié)果顯示，本文方法在Dose score和DVH score上取得最小值，顯著優(yōu)于其他SOTA方法（表1）。可視化結(jié)果顯示，本文方法預(yù)測劑量分布圖更精確，尤其在射線束區(qū)域和ROIs（圖2），另外，本文方法預(yù)測DVH曲線更逼近真實DVH曲線（圖3）。

表1 目前該領(lǐng)域SOTA方法在Dose score 和DVH score上的量化結(jié)果分析Tab.1 Comparison of Dose scores and DVH scores of the SOTA methods with the proposed

圖3 不同方法預(yù)測的DVH曲線圖的可視化結(jié)果Fig.3 Visualization of dose volume histogram predicted using different state-of-the-art methods.Dashed line represents the prediction,and solid line is the ground truth.

2.2 消融實驗

針對本文方法的消融實驗分析，量化結(jié)果顯示（表2），在級聯(lián)網(wǎng)絡(luò)的基礎(chǔ)上，引入BDP后可大幅提升劑量預(yù)測的精確性（Dose score提升了46.1%）。在此基礎(chǔ)上，增加LvDVH后DVH score提升了10.8%，使ROIs內(nèi)的劑量預(yù)測更精確，并顯著提升劑量預(yù)測的計算效率。增加LcDVH后Dose score提升了12.2%，并且DVH score提升了20.2%。另外，引入Le后，劑量預(yù)測精確性得到進一步的提升。

表2 對于本研究方法在Dose score和DVH score上的消融實驗結(jié)果分析Tab.2 Dose score and DVH score of the proposed method in ablation experiments

3 討論

與4種SOTA方法（V-Net［22］，DCNN［15］，HD U-net［23］，C3D［13］）的對比結(jié)果顯示，本文方法在量化指標(biāo)Dose score和DVH score上取得最優(yōu)結(jié)果，通過配對t檢驗，驗證本文方法帶來的提升具有統(tǒng)計學(xué)意義?？梢暬Y(jié)果顯示本文方法預(yù)測的劑量分布圖和DVH曲線圖，與真實結(jié)果更加逼近，更加滿足臨床放療要求，從而引導(dǎo)放療物理師設(shè)計治療計劃等作用，提升放療計劃質(zhì)量和治療效率。

本研究針對基于射線束的劑量預(yù)測（BDP）的消融實驗結(jié)果可知，相比于簡單的級聯(lián)網(wǎng)絡(luò)，BDP可以顯著提升劑量預(yù)測的精確性，說明射線束分割掩膜的引入，便于網(wǎng)絡(luò)提取到更多局部區(qū)域的劑量分布特征。針對DVH校正的消融實驗結(jié)果可知，引入ROIs內(nèi)劑量分布等先驗特征，使模型在ROIs 內(nèi)的預(yù)測精確性得到提升。另外，相比于Nguyen等［14］提出的可微分DVH損失函數(shù)，計算時間從“6 min/輪次”降為“2 min/輪次”，說明可以有效提升劑量預(yù)測效率。針對邊界增強約束的消融實驗結(jié)果可知，引入三維梯度損失函數(shù)，可以進一步提升模型的預(yù)測性能。綜上消融實驗證明，本文方法中的各個組成部分是有效的，可以實現(xiàn)放療計劃自動預(yù)測的目標(biāo)。

為了進一步提升劑量預(yù)測模型性能，實現(xiàn)“All in One”的一站式全自動放療，本文方法存在一定的局限性，主要包括以下幾個方面：（1）數(shù)據(jù)集單一：僅使用了頭頸癌放療計劃數(shù)據(jù)集，模型的泛化性有限；（2）擬合的射線束分割掩膜具有一定的誤差：實際臨床放療中，射線束是通過點源形式形成，然而我們以線源的形式擬合射線束分割掩膜，導(dǎo)致擬合出的射線束分割掩膜形狀和位置與真實射線束存在一定的差異。因此，在未來的工作中，首先，我們會擴大數(shù)據(jù)集改進所提出方法的泛化性。然后，如何更加精確地擬合射線束分割掩膜，將成為我們未來工作的重點。