基于臨床質子加速器的DeepPlan筆形束模型構建

張增鵬，陶寅，劉紅冬，陳志，徐榭，4，裴曦，4

1.中國科學技術大學核科學技術學院，安徽合肥，230025；2.中國科學技術大學核醫學物理研究所，安徽合肥230025；3.中山大學附屬腫瘤醫院放療科，廣東廣州510060；4.安徽慧軟科技有限公司，安徽合肥230088

前言

相比光子，質子放射治療具有獨特的Bragg 峰，在靶區給予更高劑量的同時，降低對腫瘤后端器官的輻射損傷［1-2］，因此在全球范圍高速發展。根據國際粒子治療協作組（Particle Therapy Co-Operative Group,PTCOG）的統計，截止2020年3月份，全球正在運營的質子治療中心已達91 個，我國正在運營的質子放療中心僅有2家，仍有7家正在建設中［3-4］。可以看到，隨著越來越多質子中心的投入使用，質子治療在放射治療中的比重會逐漸增大［5］。與此同時，與質子治療相配套的治療計劃系統也會越來越受到重視。目前，質子治療計劃系統幾乎被國外各大廠商壟斷，如瑞典RaySearch 的RayStation，美國Varian 的Eclipse，瑞典Elekta 的Xio 和荷蘭Philips 的Pinnacle等［6-7］。我國質子放療計劃研究起步較晚，國內擁有自主知識產權的放療系統主要有中國科學技術大學放射醫學物理中心研發的DeepPlan 和中國科學院核能安全技術研究所研發的KylinRay［8-9］。兩者質子模塊均采用筆形束算法，通過質子束在體模中的縱向積分深度劑量和橫線束斑展寬計算體模中各體素的劑量沉積。目前由于缺乏質子放療臨床數據，國產質子放療計劃系統開展的臨床數據對比實驗較少。

本研究旨在研究DeepPlan 質子模塊中建立佛羅里達大學質子放療中心（University of Florida Health Proton Therapy Institute，UFHPTI）質子加速器筆形束模型，通過對比UFHPTI質子加速器在筆形束掃描模式下的臨床實驗數據，驗證模型構建的準確性并與UFHPTI配套的商用放療計劃軟件RayStation就前列腺癌臨床數據進行劑量對比。

1 材料與方法

1.1 DeepPlan質子模塊

1.1.1 模塊簡介DeepPlan 是由中國科學技術大學放射醫學物理中心自主研發的放療計劃系統［9-12］。其質子模塊采用筆形束算法，公式如下：

其中，D（x,y,z,x0,y0,E）表示能量為E，中心坐標為（x0,y0）的質子束在模體位置（x,y,z）的劑量沉積，IDD是質子束在縱向的積分深度劑量，利用蒙卡軟件TOPAS模擬得到，w表示核散射權重，Gmcs表示多庫倫散射的高斯分布，Gnuc表示核散射的高斯分布，σmcs是能量為E的質子束由于多庫倫散射在x軸或y軸方向上高斯分布的標準差。σnuc是能量為E的質子束由于核散射在x軸或y軸方向上高斯分布的標準差，這些參數根據Soukup等［13］、Gottschalk等［14-15］的工作計算得到。

DeepPlan 質子模塊的工作流程為：首先通過數據讀取得到模體信息和質子束信息，包括模體位置、大小、體素材料、等中心點位置、質子能量、入射方向、初始束斑大小和質子束權重等，并導入計算所需要的數據庫，如上述計算所需要的IDD 數據庫、多庫倫散射數據庫、核散射數據庫和核散射權重數據庫等。接著光線追蹤通過設定兩個位置點P1、P2計算質子束穿過體素的幾何距離和等效水深度。最后劑量計算根據等效水深度得到計算所需的數據參數，對每個質子束利用公式（1）～（3）得到體模的總劑量。這3個部分實現筆形束掃描下質子劑量的快速計算。

1.1.2 模型構建

1.1.2.1 數據讀取數據庫構建的準確性是劑量計算的前提條件，比如IDD 數據庫構建時不同初始能散和角度發散質子束的IDD 是不一樣的。圖1是用蒙卡軟件TOPAS模擬的200 MeV質子在不同能散下的IDD數據曲線。由圖可知，相同能量不同能散質子束的IDD 曲線在Bragg 峰前端有一定差距，能散越小，IDD 曲線越低，這樣根據式（1）計算的劑量會比實際情況偏低。本研究結合UFHPTI 質子加速器信息并結合劉紅冬等［16］的研究，插值計算得到加速器輸出各種能量的參數信息，構建包括IDD在內的數據庫。

圖1 200 MeV質子束在不同能散的積分深度劑量Fig.1 Integrated depth dose of 200 MeV proton beam in different energy spreads

1.1.2.2 光線追蹤光線追蹤是劑量計算的重要一環，算法程序根據Siddon［17］、Han 等［18］的工作實現。該算法認為構成體模的每一個體素不是獨立的，而是一系列等距平行平面的正交集焦點體積，加上質子束穿過體模時的兩個位置點信息，可實現位置點連線上各體素幾何距離和等效水深度的快速計算。需要注意的是在實際筆形束掃描模型中，質子加速器通過X 和Y 方向掃描磁鐵實現對質子束位置的精確控制，因而有X 和Y 兩個焦點，如圖2所示。以UFHPTI質子加速器為例，其X、Y 兩個焦點與等中心平面的距離分別為2 329 和1 934 mm。因此，在設置兩個位置點進行光線追蹤時必須考慮X、Y 方向之間的差異。

1.1.2.3 劑量計算劑量計算過程需要考慮的因素較多，包括空氣處理和散射修正等，公式如下：

圖2 筆形束掃描模型虛擬源示意圖Fig.2 Diagram of virtual source in pencil beam scanning model

其中，式（4）用于在空氣中進行劑量計算，主要是因為在空氣中等效水深度為0，且在空氣中只有多庫倫散射沒有核散射。式（5）用于在非空氣體模中進行劑量計算，與式（1）相比，該公式增加了多庫倫散射修正項δmcs和核散射修正項δnuc，表達式如上所示。在（8）～（9）式中多庫倫散射修正因子cmcs和核散射修正因子cnuc分別是與多庫倫散射標準差σmcs和核散射標準差σnuc有關的常數項。

圖3是劑量計算示意圖，在參數導入之后首先應該判斷各體素是否為空氣，判別方法為計算各體素的質量組織本領之和T，若T＜10-6，則認為在空氣中，采用式（4）進行劑量計算；否則，采用式（5）進行劑量計算。

圖3 劑量計算示意圖Fig.3 Diagram of dose calculation

1.1.3 測試環境本文完成DeepPlan質子模塊臨床測試的計算機配置如下：操作系統為Windows10 專業版（64 位），CPU 為Intel(R)CoreTMi7-87100K、六核十二線程、16 G內存，GPU 為NVIDIA GeForce GTX10606GB。DeepPlan 質子模塊有c++和c#兩個版本，本研究基于c#版本在軟件Visual Studio 2017上運行。

1.2 UFHPTI質子加速器臨床實驗數據測量

1.2.1 IDD 數據測量UFHPTI 使用IBA 公司的BPC探測器（Bragg Peak Chamber,model 34070,PTW,Freiburg,Germany）在IBA 水箱中進行積分深度的測量，共有30 組質子能量，分別為：226.35、223.91、219.06、214.12、209.16、204.54、199.63、195.1、190.18、185.74、180.9、176.21、171.89、167.79、163.02、158.09、153.04、149.07、144.81、140.63、136.07、131.51、127.38、122.46、118.69、115.09、110.04、106.57、102.36、99.59 MeV 等。針對測試數據，在DeepPlan質子模塊中計算上述30組能量的質子束在水模中的劑量沉積，體模大小為100 mm×100 mm×400 mm，每個體素的大小為1 mm×1 mm×1 mm，對垂直于質子束方向的每個平面進行劑量求和，利用最大值進行歸一即可得到在縱向的積分深度劑量曲線。

1.2.2 空氣中束斑尺寸測量UFHPTI 使用IBA 公司的Lynx PT 閃爍探測器進行束斑尺寸的測量。30 組能量在上面已經提過，5 個測量位置分別為距離等中心平面20、10、0、-10、-20 cm。為了將二維數據轉成一維，需要將測試數據導入IBA的MyQA軟件提取X或Y方向上的橫向劑量分布，并做歸一處理。

針對測試數據，在DeepPlan 質子模塊中計算30組能量的質子束在體模為空氣中的劑量沉積，平面大小為40 mm×40 mm，平面精度和Lynx PT 的精度保持一致，為0.25 mm×0.25 mm，將計算結果在X 或Y方向求和得到一維數據，并做歸一處理。

1.2.3 SOBP 和橫向劑量的測量為保證劑量分布在縱向和橫向上與靶區保持適形，需要使用不同能量、不同權重、不同位置的質子束。在SOBP 測量時UFHPTI 使用IBA 公司的多層電離室Zebra 測量中心軸上的絕對劑量；在橫向劑量測量時UFHPTI 使用Lynx PT 閃爍探測器進行測量。針對測量數據，在DeepPlan 質子劑量模塊根據提供的質子束能量信息、權重信息和位置信息等，在水模中進行劑量計算，水模大小為200 mm×200 mm×400 mm，體素大小為2 mm×2 mm×2 mm。

1.2.4 前列腺癌的劑量計算RayStation 是由RaySearch 公司于2008年開始研發，2011年市場發布的治療計劃系統［19］。RayStation 的技術主要有輪廓的勾畫、形變配準及定量評估、自適應治療、劑量追蹤、多標準優化和后備計劃等。作為與UFHPTI質子加速器相配套的治療計劃系統，RayStation 在筆形束模型下的質子劑量計算有兩種算法，一種是筆形束解析算法，另一種是蒙卡算法。筆形束算法比蒙卡算法計算速度快，但準確性低于蒙卡算法：對于大量非均勻介質體模的計算和所有質子束使用射程移位器的情況等，主要以蒙卡算法為主。本文以RayStation 蒙卡算法的計算結果進行對比。在UFHPTI 的兩個前列腺癌臨床放療計劃中，相關信息如下：病例1的體模大小為512 mm×512 mm×361 mm，體素大小為1.171 9 mm×1.171 9 mm×1.000 0 mm，270°和90°方向的質子束總數分別為1 039 和977；病例2的體模大小為512 mm×512 mm×293 mm，體素大小為1.171 9 mm×1.171 9 mm×1.000 0 mm，270°和90°方向的質子束總數分別為1 307 和1 303。由于RayStation 中給出的三維劑量信息體素大小為2 mm×2 mm×2 mm，為了進行gamma 通過率對比，需要對病例1 和2 的體模按相同體素大小進行重采樣。最后使用PTW 公司的VeriSoft 軟件對二者的計算結果進行gamma分析。

2 結果與討論

2.1 IDD結果比較

為展示方便，圖4所示將30 組IDD 數據從高到低分為6 組（a～f 組），每組5 個能量進行比較，每組能量的IDD 曲線都利用該能量下最大值進行歸一，深度方向進行適當截取，DeepPlan 質子模塊計算結果用DeepPlan 表示，UFHPTI質子加速器臨床實驗數據用Measure 表示。從圖中可以看到，DeepPlan 質子模塊計算結果與UFHPTI 質子加速器臨床實驗數據在Bragg 峰前后曲線相吻合。兩者平均相對誤差為0.01%，少數能量如圖4f 組中的99.59、102.36 MeV等，DeepPlan 計算結果與實驗數據在最大值位置附近出現一個0.23%的相對誤差。出現這個情況的原因可能是測試數據最小深度間隔為3 mm，精度較低，導致在測量峰值最大位置時出現一定誤差。

圖4 30組不同能量下IDD結果對比Fig.4 Comparison of integrated depth dose under different energies

2.2 空氣中束斑尺寸結果對比

如圖5所示，以能量為199.63 MeV 的質子束為例，藍色曲線表示UFHPTI質子加速器的臨床測試結果，紅色曲線表示DeepPlan 質子模塊的計算結果，圖中包含5個不同位置的對比結果，束斑寬度由小到大分別為距離等中心平面20、10、0、-10、-20 cm 處的結果。對于某一位置的計算結果和測量數據，導入MATLAB R2018b 版本中，使用軟件自帶的高斯擬合函數求出束斑大小。圖6是不同能量質子束在等中心平面束斑大小的對比結果，大部分相對誤差在0.5%以內，只有163.02 MeV 比實際測試略低一些，相對誤差為1.13%。除此之外，表1展示距離等中心平面20、10、-10 和-20 cm 處的計算結果。從表中可知，DeepPlan 質子模塊的計算結果與臨床實驗數據大部分相對誤差在0.3%以內，平均相對誤差為0.15%。可以看到，由于質子束在空氣和體模中的輸運不同，在劑量計算時將二者分開可使結果更加準確。

圖5 199.63 MeV質子束在不同位置結果對比Fig.5 Comparison of 199.63 MeV proton beams at different positions

圖6 不同能量質子束在等中心平面的束斑尺寸Fig.6 Spot size of proton beams with different energies in the isocenter plane

表1 30組不同能量質子束在空氣不同位置的束斑比較Tab.1 Comparison of the spot size of 30 groups of proton beams with different energies at different positions in the air

2.3 SOBP和橫向劑量結果對比

圖7是縱向和橫向的結果對比，黑色曲線表示UFHPTI 質子加速器的臨床測量結果，紅色曲線表示DeepPlan 質子模塊的計算結果。由圖7a 在縱向擴展布拉格峰的對比中，可以看到實現射程為200 mm，展寬100 mm 目的。DeepPlan 劑量計算結果與臨床實驗數據平均相對誤差1.07%，最大相對誤差3.91%。圖7b 展示的是等中心平面橫向劑量分布的對比，等中心平面在距離水模表面15 cm 的位置，可以看到平坦區域達到100 mm 目的。DeepPlan 劑量計算結果與臨床實驗數據平均相對誤差1.92%，最大相對誤差4.09%。

圖7 縱向和橫向結果對比Fig.7 Comparison of longitudinal and transverse results

2.4 前列腺癌結果對比

DeepPlan與RayStation的結果在該平面二維gamma通過率達到了97.1%，但在末端右下角未通過，出現一些冷點和熱點，見圖8。可能原因是相比DeepPlan解析算法，RayStation采用蒙卡算法，對非均勻介質部分計算結果更好，這是后續工作需要改進的地方。表2具體展示了兩個前列腺癌病例在兩個射野方向的三維gamma通過率，測試標準為3 mm/3%。由表2可知每個射野方向在大于10%處方劑量區域里的gamma通過率都達到95%以上，其中病例1在兩個方向的gamma通過率分別為96.4%和97.5%，病例2在兩個方向的gamma通過率分別為99.3%和98.9%。

圖8 DeepPlan與RayStation關于前列腺癌的二維劑量結果比較Fig.8 Comparison of two-dimensional dose calculation results of DeepPlan and RayStation on prostate cancer

表2 DeepPlan與RayStation關于前列腺癌的gamma分析結果Tab.2 Gamma analysis results of DeepPlan and RayStation on prostate cancer

3 結語

本文對DeepPlan 質子模塊在數據讀取、光線追蹤和劑量計算等方面建模，通過與UFHPTI質子加速器在30組IDD對比、30組空氣中質子束斑發散對比、SOBP 和橫向劑量分布對比，驗證DeepPlan 質子模塊模型構建的準確性。同時以兩個前列腺癌的放療計劃為實例，與商用放療計劃系統RayStation 進行gamma 分析，表明DeepPlan 質子模塊在筆形束掃描模式下可以初步應用于臨床前列腺癌的劑量計算。后續還將繼續對治療計劃系統DeepPlan 質子模塊進行改進，包括對非均勻介質的修正和加入射程移位器的劑量計算等，對更多類型的病例進行測試分析，進一步增強DeepPlan 質子模塊在臨床應用的可行性。

感謝UFHPTI提供的臨床實驗數據。