本文基于點核疊加的劑量計算模型,提出了一種由點核疊加構建的筆形束核進行劑量計算的方法,這種方法可大大提高調強放射治療(IMRT)優化迭代過程的劑量計算速度,使得基于點核疊加技術的計劃系統得以集成直接孔隙的逆向調強技術。通過對模體及臨床實際病例進行試驗,結果表明此方法在提高劑量計算速度的同時,也保證了迭代過程劑量計算的精度,從而保證了與使用精確劑量計算模型得到的優化結果相一致,完全可用于調強優化過程中的劑量計算。在基于點核疊加劑量計算模型的計劃系統中使用該方法,可以避免重新研發筆形束劑量計算模型以及產品維護造成的成本大量增加。
引用本文： 鄒豪, 王玉. 基于點核構建的筆形束模型在逆向調強中的應用. 生物醫學工程學雜志, 2014, 31(1): 166-171. doi: 10.7507/1001-5515.20140033 復制
引言
調強放射治療(intensity modified radiation therapy,IMRT)[1]是通過調整多個照射野內射束的強度分布，在靶區得到高度適形劑量的同時使周圍正常組織盡可能受到較少的照射，從而達到提高治療增益比的目的。IMRT對提高腫瘤的局部控制率以及降低正常組織并發癥發生率起到了非常大的作用，因此臨床應用也越來越多。
調強優化的算法根據其優化求解方式，基本可以分為確定算法和隨機算法兩類。確定算法主要包括梯度法[2-3] 、二次規劃、線性規劃、迭代法等。隨機算法包括模擬退火[4]和遺傳算法[5]等。每種方法都有各自的優點和缺點：梯度法優化速度快，但容易陷入局部極小點；模擬退火和遺傳算法是全局優化算法，可以避免陷入局部極小點，但是迭代次數多，優化速度慢。
IMRT從實現方式上可分為兩大類：一類是連續或螺旋式的斷層放療，此類又可分為將射束準直為窄條形或用機架旋轉進行射束調強兩種，此類調強也叫扇束調強；第二大類是基于多葉準直器(multileaf collimator,MLC)的射束調強，此類調強也叫錐束調強。錐束調強可分為靜態和動態兩種，其區別是出束時葉片是否在伴隨運動。在靜態錐束調強中，調強束是用很多MLC 的均勻束分段迭加實現的，所以也叫“停-射”(Step and Shoot)模式；
傳統的調強計劃使用的是兩步法，即首先對每個射野進行強度分布優化，得到射野強度分布圖，然后將強度分布圖轉換為子野序列。兩步法優點是可以支持動態錐束調強，缺點是靶區適形度不好。而且對于形狀不規則的靶區，如鼻咽癌等，兩步法優化得到的子野數及機器跳數(monitor units,MU)非常多，所需治療時間很長。而直接機器參數優化法(direct machine parameter optimization,DMPO)通過直接計算子野的形狀及權重來得到一個可實施的調強治療計劃。它的特點是物理師可以預先設置計劃的子野數，從而有效地減少子野數量，減少治療時間，提高加速器的工作效率。同時由于MLC葉片位置可直接優化得到，避免了子野分割造成的適形度誤差。因此，近年來基于隨機算法的直接子野優化方法越來越受到重視，其中尤以基于模擬退火優化算法的直接孔隙優化(direct aperture optimization,DAO)方法為代表。但由于隨機算法迭代次數很多，而每次迭代都需要更新劑量分布，因此提高迭代過程的劑量計算速度對優化時間影響最大。在迭代過程中，為了提高劑量計算速度，一般會采用筆形束劑量計算模型，但由于基于點核疊加的劑量計算模型具有更高的精度，因此研究如何基于點核疊加構建筆形束進行優化迭代過程的劑量計算有非常大的價值，在保證系統兼容性的前提下，有望提高迭代過程中劑量計算的精度和計算速度。
本文基于已有的點核疊加模型的計劃系統，在調強優化迭代過程，提出了一種由點核疊加構建的筆形束進行劑量計算的方法，在提高了劑量計算速度的同時，也保證了迭代劑量計算的算法精度及其與精確計算結果的一致性，具體見2.3節，使得DMPO在基于點核疊加劑量計算模型系統中得以實現，避免了重新研發筆形束劑量計算模型的大量工作。
1 直接孔徑優化
直接孔徑調強是一步調強法，利用模擬退火算法直接優化子野形狀和對應照射強度，從而得到計劃所需的子野信息，這種方法比傳統的兩步法減少了子野數和MU，提高了實施照射的效率。
直接孔徑調強主要有三方面工作：① 模擬退火優化方法；② 劑量計算；③ 目標函數與約束。
1.1 模擬退火優化方法
模擬退火算法是20世紀80年代發展起來的一種用于求解大規模優化問題的隨機搜索算法，此問題的自變量組合狀態空間是MLC葉片的位置和此子野的強度。優化的目標函數是處方劑量與模體劑量之間的差值函數。設定初始溫度、子野葉片初始狀態及權重，自變量狀態改變時計算模體的劑量分布，然后得到目標函數值，與上個狀態的目標函數值進行比較，如果目標函數值減小，接受這個狀態，否則以概率P接受。進行一定次數的狀態變換后，降低溫度，重復上述步驟，直至溫度或目標函數值達到停止條件。
影響模擬退火算法速度和質量有三個關鍵因素：① 狀態跳躍方式；② 接受方式；③ 退火機制。
根據以上三個因素的不同可以構造不同的模擬退火算法。
① 跳躍方式
葉片位置和大小的改變服從高斯分布，高斯寬度的改變如下：
| $\sigma =1+\left( A-1 \right){{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{Step}}},$ |
式中A為初始寬度，nsucc為成功迭代次數，T0step為冷卻比率[4, 6]。
② 接受方式
如果葉片位置的改變與MLC的限制條件相抵觸，則重新調整葉片位置，直到可以接受為止。如果葉片位置可以接受，則計算劑量分布，確定目標函數的值，與上次迭代的目標函數值進行比價，如果減小，則接受這次迭代，如果增加，則接受概率P以標準玻爾茲曼模擬退火公式給出：
| $P=2B\frac{1}{1+{{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{prob}}}},$ |
式中B為初始概率，nsucc為成功迭代數，T0prob為降溫比率。
③ 退火機制
| ${{T}_{n}}={{\alpha }^{n}}{{T}_{0}},$ |
式中T0為初始溫度。
1.2 劑量計算
由于模擬退火算法需要進行至少幾千次的迭代，需要快速的劑量計算算法保證整個優化的速度，因此對于迭代過程中的劑量計算一般都采用筆形束模型。
筆形束算法的基本思想是利用單一方向的筆形束核做二維卷積。基于這種思想，發展了一系列的基于筆形束核的二維卷積算法[7-9]。由于筆形束核是二維卷積，與點核疊加相比，速度更快。關于迭代過程中基于筆形束模型的劑量計算將在第二節詳細敘述。
1.3 目標函數與約束
目標函數是將治療計劃輸入的約束簡化為某一函數數學表達式，優化過程即是尋找這一目標函數的極小值，其基本形式由公式(4) 給出：
| $\underset{m\ge 0}{\mathop{\min }}\,{{\sum\limits_{m}{{{\theta }_{m}}\left( d_{l}^{c}-d_{l}^{p} \right)}}^{2}},$ |
其中θm為感興趣區(region of interest,ROI)的權重，dlc為ROI體元計算得到的劑量，dlp為ROI體元處方劑量。
對于MLC物理限制，如葉片最大行程、相鄰葉片位置差等約束，如在新的狀態下超出物理限制即不接受。同樣地，對一些特殊器官，如脊髓，即使是一小部分的缺失也會造成整個器官的損傷，我們就要設置一個最大劑量確保劑量不能超出。圖 1是直接孔徑調強算法流程。
圖1
直接孔徑調強算法流程圖
Figure1.
Process flowchart of DAO IMRT algorithm
2 基于點核疊加構建筆形束模型
2.1 點核疊加模型
迭加的原理可以被用于外部射束劑量分布的計算[10-11]。利用這個方法，單能光子入射到一個均勻吸收介質中得到的劑量，可以由射束的密度分布與一個空間不變核卷積得到：
| $D\left( r \right)=\int\limits_{v}{T\left( {{r}'} \right)A\left( r-{r}' \right)}{{d}^{3}}{r}'$ |
式中D(r)表示的是吸收介質中點r處的劑量，T(r′)表示Terma，是主光子在點r′處的體元d3r′中作用釋放的所有能量，A(r-r′)為劑量伸展核。幾何位置關系如圖 2所示。
圖2
點核卷積/疊加原理
Figure2.
Point kernel convolution/superposition theory
單能光子在點r′處的Terma值T(r′)可以按照如下公式解析地計算：
| $T\left( {{r}'} \right)={{\left( \frac{{{r}_{0}}}{{{r}'}} \right)}^{2}}\left( \frac{\mu }{\rho } \right)E{{\Phi }_{0}}\exp \left( -\mu {{\left| {{r}'} \right|}_{true}} \right),$ |
式中μ為介質的線性衰減因子，E為光子能量，Φ0為沿r′方向上模體表面的光子能量注量，ρ為介質的密度，｜r′｜true為在r′方向上在模體中的實際長度。
2.2 筆形束核疊加模型
在不同的算法中，筆形束核的取法也不同，可以是解析核，或者直接由蒙特卡洛程序生成的線核。下面簡要介紹一下有限筆形束算法的思想[12]：
模型假設：
(1) 射束平行照射；
(2) 射束能被從幾何意義上分割成有限個有限尺寸大小的射束，叫做有限大小筆形束(finite size pecil beam,FSPB)；
(3 )FSPB由蒙特卡洛方法產生；
(4 )FSPB劑量分布的迭加產生整個射束的劑量分布；
(5) 每個FSPB 都有一個與它在整個射束中的位置有關的權重注量；
(6) 每個放射束都有一個能譜，可以由有限個離散的能量段加權來模擬。
平行射束被分成多個固定大小的窄射束(筆形束)，每個窄射束可以作為一個獨立的射束，根據它的能譜，就可以利用蒙特卡洛方法計算出筆形束在均勻模體中單位注量下的劑量分布，這個劑量分布就是窄射束對應的筆形束核。如圖 3所示，分割后的窄射束(I,J)對應不同的能量通量強度，但具有相同的筆形束核。因此在知道各個筆形束位置的注量后，通過對各個核的迭加就可以得到整體射束下的劑量分布。
圖3
筆形束核卷積/疊加原理
Figure3.
Pencil beam convolution/superposition theory
模體內每個體元的劑量Dp等于每個筆形束核對這個體元的劑量貢獻的和，這個求和的過程就是劑量迭加過程。假設射束被分成n×m個筆形束，則某點P的劑量為：
| ${{D}^{p}}=\sum\limits_{IJ=0}^{nm}{D_{IJ}^{p}{{W}_{IJ}}}$ |
這個權重WIJ可以理解為筆形束(I,J)的能量注量強度。
2.3 基于點核構建筆形束模型
利用蒙特卡洛方法產生的單能筆形束核或解析核進行劑量疊加運算模型，需要重新處理多能核、核傾斜、密度修正等問題。對于那些已經建立的基于點核疊加劑量計算模型的計劃系統，重新建立一套筆形束模型的研發投入是非常大的。筆形束模型假設條件與點核疊加模型有區別，尤其是對于非均勻修正區別很大。兩種模型的自身差別會造成基于筆形束模型優化得到的劑量分布結果與系統原有的基于點核模型得到的結果有一定差別，這種差別會因為迭代過程的劑量計算模型簡化被放大，也就造成優化結果不精確，從而無法穩定達到優化目標。
圖 4(a)是射野針對整個靶區適形后得到的MLC成野，圖 4(b)是在射野MLC成野范圍內基于MLC葉片寬度對能量通量進行離散處理。單位能量通量下，每一個離散單元野Bij在受照射體內都對應一個劑量分布Kij。二級準直器成野到射野MLC適形野的跟隨野，如圖 4(a)二級準直器開野位置，MLC成封閉野狀態，在單位能量通量下，在受照射體內對應劑量分布S。
圖4
基于MLC適形野離散能量通量
(a)MLC適形野;(b)能量通量離散
Figure4. Discrete energy fluence based on MLC conform field(a) MLC conform field; (b) discrete energy fluence
考慮到本射野內的調強子野是本射野適形野的子集，對于不同子野下受照射體內點的劑量分布計算如下：
| ${{D}^{p}}={{S}^{p}}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{M}_{ij}},$ |
式中Dp為模體內P點處體元的劑量值，Sp為MLC封閉狀態下模體內P點處體元的劑量值。Kij為離散單元野Bij對模體內P點處體元的劑量貢獻。Sij為MLC封閉狀態下離散單元野Bij對模體內P點處體元的劑量貢獻。Mij對應離散單元野Bij狀態，且
| ${{M}_{ij}}=\left\{ \begin{matrix} 0\left( {{B}_{ij}}被MLC葉片遮擋 \right) \\ 1\left( {{B}_{ij}}未被MLC葉片遮擋 \right) \\ \end{matrix} \right.$ |
考慮到調強優化過程迭代次數較多，為了減少迭代過程的劑量計算時間，一般Kij大小限制在Bij沿射束投影范圍內的體元或小范圍擴展。此外考慮到迭代過程基于簡化的筆形束模型劑量計算是有誤差的，為了減小與精確模型劑量計算的誤差，在迭代過程中基于當前的子野形狀會進行精確模型的劑量計算，得到在受照射體內對應的劑量分布Cij。基于當前子野下精確模型的劑量分布，下一次迭代的劑量計算為
| ${{D}^{n+1}}=C_{ij}^{n}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{L}_{ij}},$ |
式中Cijn為第n次精確劑量計算得到的模體內P點處體元的劑量值，Dn+1為第n次精確劑量計算后的下一次迭代模體內P點處體元的劑量值，Lij對應第n+1次迭代離散單元野Bij狀態，且
| ${{L}_{ij}}\left\{ \begin{matrix} -1,\left( {{B}_{ij}}由打開變為遮擋 \right) \\ 0,\left( {{B}_{ij}}狀態未改變 \right) \\ 1,\left( {{B}_{ij}}由遮擋變為打開 \right) \\ \end{matrix} \right.$ |
3 試驗數據與結果
基于實際的患者計劃數據，我們分別用點核卷積疊加和本文提出的方法做了相應的劑量計算，并對兩種劑量計算得到的結果做了對比，對比結果如圖 5所示。
圖5
劑量比較
(a)中心軸深度劑量比較;(b)X軸離軸劑量比較;(c)Y軸離軸劑量比較
Figure5. Results of dose comparison(a) results of central axis depth dose comparison; (b) results of X off axis dose comparison; (c) results of Y off axis dose comparison
圖 5(a)表示經過等中心中心軸深度劑量。圖 5(b)表示經過等中心X軸離軸劑量分布。圖 5(c)表示經過等中心Y軸離軸劑量分布。其中X、Y軸為與中心軸形成笛卡兒坐標系的垂直軸，縱坐標表示本文2.3節提供的方法計算得到的結果與精確點核疊加模型計算結果的比值乘以100%。綠色曲線為點核卷積疊加得到的結果，紅色曲線為本文提出的方法得到的結果。
為了進一步驗證本文提到的劑量計算方法在逆向優化設計中的應用，我們對模體數據進行了相應的實際優化計算。模體CT圖像數據如圖 6所示，其中，黃色區域為靶區計劃靶體積(plan target volume,PTV)，綠色區域為危及器官Cord。
圖6
模體CT圖像數據
Figure6.
Phantom CT image
在優化過程中我們選擇了8個野進行照射，每個野設定為6個子野。優化結果的劑量體積直方圖(dose volume histogram,DVH)如圖 7所示，其中紅色曲線為優化前DVH曲線，藍色曲線為優化后DVH曲線。從圖中可以看出，危及器官和靶區的DVH曲線都有了相應的改進，尤其是危及器官的改進非常明顯。
圖7
優化前后DVH結果
Figure7.
DVH results before and after optimization
為了進行定量評估，對本文提及的劑量計算方法和精確點核疊加劑量計算的時間進行了對比試驗。患者CT圖像數據圖像層數為76，層間距0.5 cm，劑量網格0.4 cm，計劃按Gantry角度分布設計了8個射野，每個射野為10×10方野照射，分別用精確點核疊加及本文2.3方法進行劑量計算。軟件運行環境為臺式機Windows XP系統，硬件配置為Intel? CoreTM2 Duo CPU E8400@3.00GHz，內存2 GB，表 1是相關對比數據。
表1
劑量計算時間比較
Table1.
Comparison of dose calculation time
從表 1中數據可以看出，本文提及的方法基本上可以把迭代過程中單野的劑量計算控制在1 s以內。
4 結論
逆向調強的優化過程是基于目標函數確定子野形狀及權重的最優解，在迭代過程中對于劑量計算并不要求完全精確，只要基于劑量計算結果的目標函數值與要求的收斂預期一致，就可接近或達到優化目標的設定。
考慮到上述原因及在優化結束后可基于優化結果進行精確劑量計算后評估優化結果是否滿足要求，我們認為在優化迭代過程中的劑量計算更重要的是保證一定精度的前提下讓劑量計算速度更快，這樣才能提高逆向調強的優化速度。從試驗數據及結果看，本文提出的方法可以保證優化迭代過程中的劑量計算精度滿足上述要求。使得直接機器參數優化法(DMPO)在基于點核疊加劑量計算模型系統中得以實現，避免了重新研發筆形束劑量計算模型的大量工作。
引言
調強放射治療(intensity modified radiation therapy,IMRT)[1]是通過調整多個照射野內射束的強度分布，在靶區得到高度適形劑量的同時使周圍正常組織盡可能受到較少的照射，從而達到提高治療增益比的目的。IMRT對提高腫瘤的局部控制率以及降低正常組織并發癥發生率起到了非常大的作用，因此臨床應用也越來越多。
調強優化的算法根據其優化求解方式，基本可以分為確定算法和隨機算法兩類。確定算法主要包括梯度法[2-3] 、二次規劃、線性規劃、迭代法等。隨機算法包括模擬退火[4]和遺傳算法[5]等。每種方法都有各自的優點和缺點：梯度法優化速度快，但容易陷入局部極小點；模擬退火和遺傳算法是全局優化算法，可以避免陷入局部極小點，但是迭代次數多，優化速度慢。
IMRT從實現方式上可分為兩大類：一類是連續或螺旋式的斷層放療，此類又可分為將射束準直為窄條形或用機架旋轉進行射束調強兩種，此類調強也叫扇束調強；第二大類是基于多葉準直器(multileaf collimator,MLC)的射束調強，此類調強也叫錐束調強。錐束調強可分為靜態和動態兩種，其區別是出束時葉片是否在伴隨運動。在靜態錐束調強中，調強束是用很多MLC 的均勻束分段迭加實現的，所以也叫“停-射”(Step and Shoot)模式；
傳統的調強計劃使用的是兩步法，即首先對每個射野進行強度分布優化，得到射野強度分布圖，然后將強度分布圖轉換為子野序列。兩步法優點是可以支持動態錐束調強，缺點是靶區適形度不好。而且對于形狀不規則的靶區，如鼻咽癌等，兩步法優化得到的子野數及機器跳數(monitor units,MU)非常多，所需治療時間很長。而直接機器參數優化法(direct machine parameter optimization,DMPO)通過直接計算子野的形狀及權重來得到一個可實施的調強治療計劃。它的特點是物理師可以預先設置計劃的子野數，從而有效地減少子野數量，減少治療時間，提高加速器的工作效率。同時由于MLC葉片位置可直接優化得到，避免了子野分割造成的適形度誤差。因此，近年來基于隨機算法的直接子野優化方法越來越受到重視，其中尤以基于模擬退火優化算法的直接孔隙優化(direct aperture optimization,DAO)方法為代表。但由于隨機算法迭代次數很多，而每次迭代都需要更新劑量分布，因此提高迭代過程的劑量計算速度對優化時間影響最大。在迭代過程中，為了提高劑量計算速度，一般會采用筆形束劑量計算模型，但由于基于點核疊加的劑量計算模型具有更高的精度，因此研究如何基于點核疊加構建筆形束進行優化迭代過程的劑量計算有非常大的價值，在保證系統兼容性的前提下，有望提高迭代過程中劑量計算的精度和計算速度。
本文基于已有的點核疊加模型的計劃系統，在調強優化迭代過程，提出了一種由點核疊加構建的筆形束進行劑量計算的方法，在提高了劑量計算速度的同時，也保證了迭代劑量計算的算法精度及其與精確計算結果的一致性，具體見2.3節，使得DMPO在基于點核疊加劑量計算模型系統中得以實現，避免了重新研發筆形束劑量計算模型的大量工作。
1 直接孔徑優化
直接孔徑調強是一步調強法，利用模擬退火算法直接優化子野形狀和對應照射強度，從而得到計劃所需的子野信息，這種方法比傳統的兩步法減少了子野數和MU，提高了實施照射的效率。
直接孔徑調強主要有三方面工作：① 模擬退火優化方法；② 劑量計算；③ 目標函數與約束。
1.1 模擬退火優化方法
模擬退火算法是20世紀80年代發展起來的一種用于求解大規模優化問題的隨機搜索算法，此問題的自變量組合狀態空間是MLC葉片的位置和此子野的強度。優化的目標函數是處方劑量與模體劑量之間的差值函數。設定初始溫度、子野葉片初始狀態及權重，自變量狀態改變時計算模體的劑量分布，然后得到目標函數值，與上個狀態的目標函數值進行比較，如果目標函數值減小，接受這個狀態，否則以概率P接受。進行一定次數的狀態變換后，降低溫度，重復上述步驟，直至溫度或目標函數值達到停止條件。
影響模擬退火算法速度和質量有三個關鍵因素：① 狀態跳躍方式；② 接受方式；③ 退火機制。
根據以上三個因素的不同可以構造不同的模擬退火算法。
① 跳躍方式
葉片位置和大小的改變服從高斯分布，高斯寬度的改變如下：
| $\sigma =1+\left( A-1 \right){{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{Step}}},$ |
式中A為初始寬度，nsucc為成功迭代次數，T0step為冷卻比率[4, 6]。
② 接受方式
如果葉片位置的改變與MLC的限制條件相抵觸，則重新調整葉片位置，直到可以接受為止。如果葉片位置可以接受，則計算劑量分布，確定目標函數的值，與上次迭代的目標函數值進行比價，如果減小，則接受這次迭代，如果增加，則接受概率P以標準玻爾茲曼模擬退火公式給出：
| $P=2B\frac{1}{1+{{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{prob}}}},$ |
式中B為初始概率，nsucc為成功迭代數，T0prob為降溫比率。
③ 退火機制
| ${{T}_{n}}={{\alpha }^{n}}{{T}_{0}},$ |
式中T0為初始溫度。
1.2 劑量計算
由于模擬退火算法需要進行至少幾千次的迭代，需要快速的劑量計算算法保證整個優化的速度，因此對于迭代過程中的劑量計算一般都采用筆形束模型。
筆形束算法的基本思想是利用單一方向的筆形束核做二維卷積。基于這種思想，發展了一系列的基于筆形束核的二維卷積算法[7-9]。由于筆形束核是二維卷積，與點核疊加相比，速度更快。關于迭代過程中基于筆形束模型的劑量計算將在第二節詳細敘述。
1.3 目標函數與約束
目標函數是將治療計劃輸入的約束簡化為某一函數數學表達式，優化過程即是尋找這一目標函數的極小值，其基本形式由公式(4) 給出：
| $\underset{m\ge 0}{\mathop{\min }}\,{{\sum\limits_{m}{{{\theta }_{m}}\left( d_{l}^{c}-d_{l}^{p} \right)}}^{2}},$ |
其中θm為感興趣區(region of interest,ROI)的權重，dlc為ROI體元計算得到的劑量，dlp為ROI體元處方劑量。
對于MLC物理限制，如葉片最大行程、相鄰葉片位置差等約束，如在新的狀態下超出物理限制即不接受。同樣地，對一些特殊器官，如脊髓，即使是一小部分的缺失也會造成整個器官的損傷，我們就要設置一個最大劑量確保劑量不能超出。圖 1是直接孔徑調強算法流程。
圖1
直接孔徑調強算法流程圖
Figure1.
Process flowchart of DAO IMRT algorithm
2 基于點核疊加構建筆形束模型
2.1 點核疊加模型
迭加的原理可以被用于外部射束劑量分布的計算[10-11]。利用這個方法，單能光子入射到一個均勻吸收介質中得到的劑量，可以由射束的密度分布與一個空間不變核卷積得到：
| $D\left( r \right)=\int\limits_{v}{T\left( {{r}'} \right)A\left( r-{r}' \right)}{{d}^{3}}{r}'$ |
式中D(r)表示的是吸收介質中點r處的劑量，T(r′)表示Terma，是主光子在點r′處的體元d3r′中作用釋放的所有能量，A(r-r′)為劑量伸展核。幾何位置關系如圖 2所示。
圖2
點核卷積/疊加原理
Figure2.
Point kernel convolution/superposition theory
單能光子在點r′處的Terma值T(r′)可以按照如下公式解析地計算：
| $T\left( {{r}'} \right)={{\left( \frac{{{r}_{0}}}{{{r}'}} \right)}^{2}}\left( \frac{\mu }{\rho } \right)E{{\Phi }_{0}}\exp \left( -\mu {{\left| {{r}'} \right|}_{true}} \right),$ |
式中μ為介質的線性衰減因子，E為光子能量，Φ0為沿r′方向上模體表面的光子能量注量，ρ為介質的密度，｜r′｜true為在r′方向上在模體中的實際長度。
2.2 筆形束核疊加模型
在不同的算法中，筆形束核的取法也不同，可以是解析核，或者直接由蒙特卡洛程序生成的線核。下面簡要介紹一下有限筆形束算法的思想[12]：
模型假設：
(1) 射束平行照射；
(2) 射束能被從幾何意義上分割成有限個有限尺寸大小的射束，叫做有限大小筆形束(finite size pecil beam,FSPB)；
(3 )FSPB由蒙特卡洛方法產生；
(4 )FSPB劑量分布的迭加產生整個射束的劑量分布；
(5) 每個FSPB 都有一個與它在整個射束中的位置有關的權重注量；
(6) 每個放射束都有一個能譜，可以由有限個離散的能量段加權來模擬。
平行射束被分成多個固定大小的窄射束(筆形束)，每個窄射束可以作為一個獨立的射束，根據它的能譜，就可以利用蒙特卡洛方法計算出筆形束在均勻模體中單位注量下的劑量分布，這個劑量分布就是窄射束對應的筆形束核。如圖 3所示，分割后的窄射束(I,J)對應不同的能量通量強度，但具有相同的筆形束核。因此在知道各個筆形束位置的注量后，通過對各個核的迭加就可以得到整體射束下的劑量分布。
圖3
筆形束核卷積/疊加原理
Figure3.
Pencil beam convolution/superposition theory
模體內每個體元的劑量Dp等于每個筆形束核對這個體元的劑量貢獻的和，這個求和的過程就是劑量迭加過程。假設射束被分成n×m個筆形束，則某點P的劑量為：
| ${{D}^{p}}=\sum\limits_{IJ=0}^{nm}{D_{IJ}^{p}{{W}_{IJ}}}$ |
這個權重WIJ可以理解為筆形束(I,J)的能量注量強度。
2.3 基于點核構建筆形束模型
利用蒙特卡洛方法產生的單能筆形束核或解析核進行劑量疊加運算模型，需要重新處理多能核、核傾斜、密度修正等問題。對于那些已經建立的基于點核疊加劑量計算模型的計劃系統，重新建立一套筆形束模型的研發投入是非常大的。筆形束模型假設條件與點核疊加模型有區別，尤其是對于非均勻修正區別很大。兩種模型的自身差別會造成基于筆形束模型優化得到的劑量分布結果與系統原有的基于點核模型得到的結果有一定差別，這種差別會因為迭代過程的劑量計算模型簡化被放大，也就造成優化結果不精確，從而無法穩定達到優化目標。
圖 4(a)是射野針對整個靶區適形后得到的MLC成野，圖 4(b)是在射野MLC成野范圍內基于MLC葉片寬度對能量通量進行離散處理。單位能量通量下，每一個離散單元野Bij在受照射體內都對應一個劑量分布Kij。二級準直器成野到射野MLC適形野的跟隨野，如圖 4(a)二級準直器開野位置，MLC成封閉野狀態，在單位能量通量下，在受照射體內對應劑量分布S。
圖4
基于MLC適形野離散能量通量
(a)MLC適形野;(b)能量通量離散
Figure4. Discrete energy fluence based on MLC conform field(a) MLC conform field; (b) discrete energy fluence
考慮到本射野內的調強子野是本射野適形野的子集，對于不同子野下受照射體內點的劑量分布計算如下：
| ${{D}^{p}}={{S}^{p}}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{M}_{ij}},$ |
式中Dp為模體內P點處體元的劑量值，Sp為MLC封閉狀態下模體內P點處體元的劑量值。Kij為離散單元野Bij對模體內P點處體元的劑量貢獻。Sij為MLC封閉狀態下離散單元野Bij對模體內P點處體元的劑量貢獻。Mij對應離散單元野Bij狀態，且
| ${{M}_{ij}}=\left\{ \begin{matrix} 0\left( {{B}_{ij}}被MLC葉片遮擋 \right) \\ 1\left( {{B}_{ij}}未被MLC葉片遮擋 \right) \\ \end{matrix} \right.$ |
考慮到調強優化過程迭代次數較多，為了減少迭代過程的劑量計算時間，一般Kij大小限制在Bij沿射束投影范圍內的體元或小范圍擴展。此外考慮到迭代過程基于簡化的筆形束模型劑量計算是有誤差的，為了減小與精確模型劑量計算的誤差，在迭代過程中基于當前的子野形狀會進行精確模型的劑量計算，得到在受照射體內對應的劑量分布Cij。基于當前子野下精確模型的劑量分布，下一次迭代的劑量計算為
| ${{D}^{n+1}}=C_{ij}^{n}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{L}_{ij}},$ |
式中Cijn為第n次精確劑量計算得到的模體內P點處體元的劑量值，Dn+1為第n次精確劑量計算后的下一次迭代模體內P點處體元的劑量值，Lij對應第n+1次迭代離散單元野Bij狀態，且
| ${{L}_{ij}}\left\{ \begin{matrix} -1,\left( {{B}_{ij}}由打開變為遮擋 \right) \\ 0,\left( {{B}_{ij}}狀態未改變 \right) \\ 1,\left( {{B}_{ij}}由遮擋變為打開 \right) \\ \end{matrix} \right.$ |
3 試驗數據與結果
基于實際的患者計劃數據，我們分別用點核卷積疊加和本文提出的方法做了相應的劑量計算，并對兩種劑量計算得到的結果做了對比，對比結果如圖 5所示。
圖5
劑量比較
(a)中心軸深度劑量比較;(b)X軸離軸劑量比較;(c)Y軸離軸劑量比較
Figure5. Results of dose comparison(a) results of central axis depth dose comparison; (b) results of X off axis dose comparison; (c) results of Y off axis dose comparison
圖 5(a)表示經過等中心中心軸深度劑量。圖 5(b)表示經過等中心X軸離軸劑量分布。圖 5(c)表示經過等中心Y軸離軸劑量分布。其中X、Y軸為與中心軸形成笛卡兒坐標系的垂直軸，縱坐標表示本文2.3節提供的方法計算得到的結果與精確點核疊加模型計算結果的比值乘以100%。綠色曲線為點核卷積疊加得到的結果，紅色曲線為本文提出的方法得到的結果。
為了進一步驗證本文提到的劑量計算方法在逆向優化設計中的應用，我們對模體數據進行了相應的實際優化計算。模體CT圖像數據如圖 6所示，其中，黃色區域為靶區計劃靶體積(plan target volume,PTV)，綠色區域為危及器官Cord。
圖6
模體CT圖像數據
Figure6.
Phantom CT image
在優化過程中我們選擇了8個野進行照射，每個野設定為6個子野。優化結果的劑量體積直方圖(dose volume histogram,DVH)如圖 7所示，其中紅色曲線為優化前DVH曲線，藍色曲線為優化后DVH曲線。從圖中可以看出，危及器官和靶區的DVH曲線都有了相應的改進，尤其是危及器官的改進非常明顯。
圖7
優化前后DVH結果
Figure7.
DVH results before and after optimization
為了進行定量評估，對本文提及的劑量計算方法和精確點核疊加劑量計算的時間進行了對比試驗。患者CT圖像數據圖像層數為76，層間距0.5 cm，劑量網格0.4 cm，計劃按Gantry角度分布設計了8個射野，每個射野為10×10方野照射，分別用精確點核疊加及本文2.3方法進行劑量計算。軟件運行環境為臺式機Windows XP系統，硬件配置為Intel? CoreTM2 Duo CPU E8400@3.00GHz，內存2 GB，表 1是相關對比數據。
表1
劑量計算時間比較
Table1.
Comparison of dose calculation time
從表 1中數據可以看出，本文提及的方法基本上可以把迭代過程中單野的劑量計算控制在1 s以內。
4 結論
逆向調強的優化過程是基于目標函數確定子野形狀及權重的最優解，在迭代過程中對于劑量計算并不要求完全精確，只要基于劑量計算結果的目標函數值與要求的收斂預期一致，就可接近或達到優化目標的設定。
考慮到上述原因及在優化結束后可基于優化結果進行精確劑量計算后評估優化結果是否滿足要求，我們認為在優化迭代過程中的劑量計算更重要的是保證一定精度的前提下讓劑量計算速度更快，這樣才能提高逆向調強的優化速度。從試驗數據及結果看，本文提出的方法可以保證優化迭代過程中的劑量計算精度滿足上述要求。使得直接機器參數優化法(DMPO)在基于點核疊加劑量計算模型系統中得以實現，避免了重新研發筆形束劑量計算模型的大量工作。
