毛凱娜，張鵬程，桂志國(guó)

1.中北大學(xué)信息與通信工程學(xué)院，山西太原030051；2.生物醫(yī)學(xué)成像與影像大數(shù)據(jù)山西省重點(diǎn)實(shí)驗(yàn)室（中北大學(xué)），山西太原030051

前言

惡性腫瘤己成為嚴(yán)重危害人類健康的一類疾病，目前治療腫瘤的三大主要手段為放射治療、手術(shù)治療和化學(xué)藥物治療。與其它兩種治療方式相比，放射治療能最大限度將照射劑量集中在靶區(qū)，在不同程度上保留被治療器官的功能，改善患者的愈后生存質(zhì)量，從而成為眾多學(xué)者的研究重點(diǎn)。近年來(lái)隨著放射物理學(xué)、計(jì)算機(jī)科學(xué)技術(shù)、醫(yī)學(xué)影像技術(shù)等學(xué)科的更新與發(fā)展，出現(xiàn)了許多新的放療技術(shù)，目前臨床研究較為廣泛的一項(xiàng)技術(shù)是調(diào)強(qiáng)放射治療（IMRT）。在IMRT 中，逆向計(jì)劃系統(tǒng)的方案優(yōu)化過(guò)程中需要頻繁計(jì)算射線在病人體內(nèi)的劑量沉積，劑量計(jì)算的精度決定方案優(yōu)化的質(zhì)量，劑量計(jì)算的速度則影響方案優(yōu)化的時(shí)間。因此，精確快速的劑量計(jì)算是放射治療系統(tǒng)的核心。

近年來(lái)，劑量計(jì)算的研究重點(diǎn)是半解析法和解析法。解析法包括蒙特卡羅法［1］（Monte Carlo,MC）和求解粒子輸運(yùn)方程的方法［2］，雖然它們能夠精確地模擬射線在病人體內(nèi)的能量沉積過(guò)程，但是由于計(jì)算量復(fù)雜，還是無(wú)法滿足方案優(yōu)化對(duì)劑量計(jì)算速度的要求，因此很多腫瘤治療系統(tǒng)仍然采用半解析法進(jìn)行逆向計(jì)劃系統(tǒng)的方案優(yōu)化。常用的半解析法是基于核模型的劑量計(jì)算方法，包括筆形束核劑量計(jì)算方法［3］和點(diǎn)核劑量計(jì)算方法［4］，其中基于點(diǎn)核劑量計(jì)算方法的計(jì)算精度較高，但計(jì)算速度較慢；基于筆形束核的劑量計(jì)算方法雖然計(jì)算精度不高，但計(jì)算速度快，因此目前在臨床上仍然被廣泛使用和研究。1992年，Ahnesj? 等［3］提出一種指數(shù)函數(shù)表示的筆形束核模型，利用三角場(chǎng)積分將主要?jiǎng)┝颗c散射劑量分開(kāi)計(jì)算。2006年和2007年，Whitton 等［5］和Zhou等［6］分別提出可編程門陣列（Field Programmable Gate Array,FPGA）的加速劑量計(jì)算方法，利用硬件實(shí)現(xiàn)劑量計(jì)算的加速。2008年，Tillikainen 等［7］提出在球殼坐標(biāo)系下計(jì)算劑量分布的筆形束核模型，同時(shí)用6個(gè)指數(shù)函數(shù)表示筆形束核的側(cè)向分量，增加劑量計(jì)算的精度。2010年，Lu［8］提出一種混合劑量計(jì)算方法，能夠兼容筆形束核劑量計(jì)算速度和點(diǎn)核劑量計(jì)算精度，不過(guò)計(jì)算過(guò)程還是比較復(fù)雜。2011年，Chen 等［9］利用GPU 對(duì)點(diǎn)核劑量計(jì)算方法進(jìn)行加速，但還是沒(méi)有達(dá)到降低算法本身的計(jì)算復(fù)雜度的目的。

由于筆形束核劑量計(jì)算中，計(jì)算量最大的部分在于重復(fù)利用射線追蹤法計(jì)算筆形束核側(cè)向分量中軸線與體素的相交情況，因而減少射線追蹤法的使用，降低軸線與體素相交的計(jì)算次數(shù)就可以降低算法復(fù)雜度。本文結(jié)合文獻(xiàn)［7］方法，提出一種降低筆形束核劑量計(jì)算復(fù)雜度的算法。在球殼坐標(biāo)系下，以碰撞點(diǎn)為中心計(jì)算劑量分布，每一層碰撞點(diǎn)處軸線與體素的相交情況都是一致的。因此每條射束只需計(jì)算初始層軸線與體素的相交情況，再根據(jù)體素的相對(duì)位置和校正因子得到其它層軸線與體素的相交情況。在這種情況下，軸線與體素相交次數(shù)的計(jì)算量減少，繼而算法復(fù)雜度降低，計(jì)算速度也會(huì)大大提高。

1 獲取核參數(shù)

1.1 球殼坐標(biāo)系

本文使用的模擬照射源為點(diǎn)源。由于在直角坐標(biāo)系中，點(diǎn)源發(fā)出的射線不是平行入射，因而傳統(tǒng)算法根據(jù)發(fā)散射線的偏轉(zhuǎn)角度，適當(dāng)旋轉(zhuǎn)中心軸處的筆形束核從而得到對(duì)應(yīng)角度的筆形束核［10-11］。然而旋轉(zhuǎn)筆形束核是非常耗時(shí)的方法。但是在球殼坐標(biāo)系下［7］，來(lái)自點(diǎn)源的射線垂直于每一個(gè)球殼，仍然可以按照平行射束的方式計(jì)算劑量分布，避免對(duì)筆形束核的旋轉(zhuǎn)，從而減少計(jì)算時(shí)間。

筆形束核的坐標(biāo)系以正交基的形式給出。筆形束核按照固定頂角進(jìn)行采樣［7］。在這種采樣方式下，由直角坐標(biāo)系到球殼坐標(biāo)系的映射關(guān)系表示為：

式中，x=(xx,xy,xz)表示直角坐標(biāo)系中的點(diǎn)，p=(px,py,pz)表示在球殼坐標(biāo)系中對(duì)應(yīng)的點(diǎn)，:代表映射關(guān)系，從直角坐標(biāo)系映射到球殼坐標(biāo)系。

1.2 筆形束核參數(shù)的獲取

通常將筆形束核分為兩部分：深度方向分量和側(cè)向分量。深度方向分量表示射束β在模體深度為pz的球殼層上釋放的所有能量，均勻介質(zhì)中的深度分量表示為：

其中，Φ 表示模體表層的能量注量分布，hβ表示基于光束β的多能核［7］。

在非均勻介質(zhì)中根據(jù)有效深度和相對(duì)電子密度［12］對(duì)深度分量進(jìn)行校正：

其中，ρ(p')表示射束與模體在實(shí)際深度層碰撞點(diǎn)p'處的相對(duì)電子密度：

其中，pz'表示射束入射到模體內(nèi)的等效深度，用d(X)表示：

其中，X為射束在模體內(nèi)入射的實(shí)際路徑。

側(cè)向分量fβ(r,pz)表示在距離射束中心軸長(zhǎng)度為r的體素中接收到的由該射束在深度為pz的球殼上釋放的部分能量分量：

在本文中我們使用了由Tillikainen等［7］提出的側(cè)向散射函數(shù)來(lái)模擬側(cè)向分量：

其中，μi是衰減因子，將在1 到200 mm 內(nèi)按等對(duì)數(shù)間隔取值得到。ci(pz)是側(cè)向散射函數(shù)的權(quán)重參數(shù)，通過(guò)對(duì)fβ(r,pz)和kβ(r,pz)進(jìn)行最小二乘法擬合得到。r是碰撞點(diǎn)處軸線與體素相交的有效長(zhǎng)度，通過(guò)射線追蹤法計(jì)算得到。

在非均勻介質(zhì)中，需要根據(jù)有效深度對(duì)應(yīng)球殼層碰撞點(diǎn)處軸線與體素的相交長(zhǎng)度以及沉積點(diǎn)的相對(duì)電子密度對(duì)側(cè)向分量進(jìn)行校正：

其中，ρ(p)是能量沉積點(diǎn)p處的相對(duì)電子密度。本文方法r值無(wú)需在有效深度下進(jìn)行校正，可以直接通過(guò)初始球殼層軸線與體素的相交情況得到。下文會(huì)詳細(xì)介紹。

最后將筆形束核中的深度分量Iβ(pz)、衰減因子μi、權(quán)重參數(shù)ci(pz)保存，為后續(xù)劑量計(jì)算的使用提供方便。

1.3 劑量計(jì)算

本文使用Ahnesj?等［4］提出的筒串卷積算法對(duì)筆形束核進(jìn)行疊加來(lái)計(jì)算模體內(nèi)各點(diǎn)處所有沉積的能量。同時(shí)使用等效路徑［13］的方法來(lái)校正非均勻介質(zhì)中的劑量分布。為了方便在計(jì)算機(jī)上實(shí)現(xiàn)，使用8個(gè)離散的筒串方向［7］進(jìn)行能量疊加。射束β在模體內(nèi)p點(diǎn)處沉積的能量表示為：

射野內(nèi)所有射束在p點(diǎn)處沉積的總能量為：

其中，S為射野內(nèi)所有射束的集合。

將單位質(zhì)量中接收到的能量定義為劑量，因此在直角坐標(biāo)系下x點(diǎn)處的劑量為：

其中，J(M(x))是兩坐標(biāo)系之間映射關(guān)系M的雅克比行列式，雅克比行列式的值為：

2 本文算法

2.1 算法改進(jìn)部分

在筆形束核劑量計(jì)算方法中，射束進(jìn)入模體之后的所有能量按筆形束核的概率分布模型進(jìn)行沉積。在筆形束核側(cè)向分量的校正中，傳統(tǒng)的筆形束核劑量計(jì)算是以能量沉積點(diǎn)為中心計(jì)算劑量分布，該方法需要重復(fù)利用射線追蹤法［14］計(jì)算射野內(nèi)所有射束在不同方向上軸線與體素的相交長(zhǎng)度，計(jì)算過(guò)程復(fù)雜，計(jì)算量較大，比較耗時(shí)。因而提出以碰撞點(diǎn)為中心計(jì)算劑量分布的改進(jìn)方法，根據(jù)射束在不同層碰撞點(diǎn)處能量散射方式的相似性，通過(guò)只計(jì)算初始球殼層碰撞點(diǎn)處的不同方向軸線與體素的相交情況的方式來(lái)減少射線追蹤的計(jì)算次數(shù)，進(jìn)而減少計(jì)算量。

在球殼坐標(biāo)系下，隨著模體深度的增加，不同球殼層體素的實(shí)際體積在增大，軸線與體素的實(shí)際相交長(zhǎng)度也相應(yīng)增加。但是其它層碰撞點(diǎn)處軸線與體素的相交長(zhǎng)度與初始層碰撞點(diǎn)處軸線與體素的相交長(zhǎng)度的比值是一個(gè)與體素深度位置相關(guān)的值，將這個(gè)比值作為校正因子提前計(jì)算。本算法計(jì)算射束在初始球殼層碰撞點(diǎn)處的不同方向軸線與體素的相交情況并保存。然后在劑量計(jì)算中根據(jù)射束與球殼層碰撞點(diǎn)的位置計(jì)算出相對(duì)偏移位置，從而獲得當(dāng)前球殼層碰撞點(diǎn)處與軸線相交體素的密度值；同時(shí)利用校正因子，得到軸線與體素相交的有效長(zhǎng)度。對(duì)于一條射束，該方法只計(jì)算了一層軸線與體素的相交情況，大大減少了射線追蹤的計(jì)算次數(shù)，并且在最大射野面積下存儲(chǔ)的一條射束8 個(gè)不同方向軸線與體素的相交比例的數(shù)據(jù)內(nèi)存大小僅為0.2 M，減少計(jì)算過(guò)程中的大量?jī)?nèi)存方面的開(kāi)銷。

對(duì)于筆形束核劑量計(jì)算模型，在直角坐標(biāo)系中將模體看做N=n×n×n個(gè)大小相等的立方體組成，每個(gè)立方體代表一個(gè)體素，將體素的中心點(diǎn)坐標(biāo)(xx,xy,xz)作為體素的索引坐標(biāo)。同時(shí)將模體中每個(gè)體素的密度定義在一個(gè)三維數(shù)組中。在計(jì)算時(shí)通過(guò)式（1）的變換將體素坐標(biāo)對(duì)應(yīng)到球殼坐標(biāo)系中。利用射線追蹤法計(jì)算初始層某一方向軸線與體素相交的有效長(zhǎng)度為：

其中，ρs(px,py,pz0)表示與軸線相交的體素(px,py,pz0)的密度值，lengths(px,py,pz0)表示在初始層軸線與體素(px,py,pz0)的相交長(zhǎng)度，n表示與軸線相交的體素個(gè)數(shù)。

在計(jì)算其它層軸線與體素的相交長(zhǎng)度時(shí)，需要首先根據(jù)碰撞點(diǎn)的位置計(jì)算出相對(duì)偏移位置，獲得待計(jì)算球殼層碰撞點(diǎn)處與軸線相交的各體素的密度值，然后利用校正因子，得到待計(jì)算球殼層對(duì)應(yīng)方向軸線與體素相交的有效長(zhǎng)度：

此處的ρs(px,py,pz)表示計(jì)算層與軸線相交的各體素的密度值。C(k)是與球殼深度有關(guān)的校正因子：

其中，pz表示待計(jì)算球殼層的深度，pz0表示初始球殼層的深度。

2.2 算法步驟

本文的算法設(shè)計(jì)流程如圖1所示，算法步驟如下：（1）輸入劑量計(jì)算所需的模體三維密度信息和射野信息。根據(jù)輸入的數(shù)據(jù)信息，利用MC 方法精確模擬筆形束核［15］，對(duì)于同一個(gè)治療頭，只需生成一次核模型即可。（2）核參數(shù)計(jì)算與保存。擬合式（6）和式（7）獲得衰減因子μi和權(quán)重參數(shù)ci( )pz，將筆形束核中的深度分量Iβ(pz)、衰減因子μi和權(quán)重參數(shù)ci(pz)存入查找表。在球殼坐標(biāo)系下，計(jì)算筆形束核模型中每條射束在初始球殼層碰撞點(diǎn)處軸線與體素的相交情況以及校正因子并保存在數(shù)組中。對(duì)于一條射束，只需利用射線追蹤法計(jì)算一次初始層碰撞點(diǎn)處不同方向軸線與體素的相交情況即可。（3）沿射束前進(jìn)方向確定射束與體素發(fā)生碰撞的位置，根據(jù)該碰撞點(diǎn)的位置讀取存儲(chǔ)的筆形束核模型參數(shù)、體素密度信息和校正因子，計(jì)算在該碰撞點(diǎn)處各軸線與體素的相交長(zhǎng)度，進(jìn)而得到該碰撞點(diǎn)處的能量沉積。（4）將步驟（3）中得到的球殼坐標(biāo)系下的能量分布轉(zhuǎn)換為直角坐標(biāo)系下的劑量分布，并與傳統(tǒng)算法和MC算法進(jìn)行比較得出結(jié)果。

圖1 算法流程圖Fig.1 Algorithm flowchart

3 實(shí)驗(yàn)結(jié)果分析與討論

3.1 實(shí)驗(yàn)材料

本文針對(duì)不同的射野大小使用的源皮距為100 cm，使用MC 方法的劑量計(jì)算是利用開(kāi)源軟件DOSXYZnrc 實(shí)現(xiàn)的［16］，傳統(tǒng)筆形束核劑量計(jì)算是根據(jù)文獻(xiàn)［7］方法自己編程實(shí)現(xiàn)。

實(shí)驗(yàn)?zāi)ｓw采用Tillikainen 等［7］提出的假設(shè)模體：（1）水模體；（2）肺阻塊模體，在水模體中距離表層5 cm 的位置放置一個(gè)厚10 cm、密度為ρw= 0.3 g/cm3的擋塊，該擋塊距離射野中心軸2 cm；（3）骨阻塊模體，在水模體中距離表層5 cm 的位置放置一個(gè)厚5 cm、密度為ρw= 1.85 g/cm3的擋塊，該擋塊距離射野中心軸2 cm。

3.2 水模體結(jié)果

在水模體中，我們使用改進(jìn)算法、MC 算法和傳統(tǒng)算法分別計(jì)算不同射野大小下［（3×3）、（5×5）、（10×10）cm2］的深度劑量曲線以及射野大小為（10×10 cm2深度為模體表層下5、10、20 cm 位置處的側(cè)向劑量曲線。

在圖2中，F(xiàn)AAA 表示改進(jìn)算法得到的結(jié)果，MC表示MC 方法計(jì)算得到的結(jié)果，AAA 表示傳統(tǒng)算法計(jì)算得到的結(jié)果，F(xiàn)s 表示射野大小，D 表示側(cè)向劑量曲線所在深度位置。不同算法得到的深度劑量曲線和側(cè)向劑量曲線分別用不同顏色和不同線型表示。

圖2a、圖2b 和圖2c 表示3 種算法分別在不同射野大?。郏?×3）、（5×5）、（10×10）cm2］下的深度劑量曲線，圖2d 為兩種算法在射野為（10×10）cm2的情況下，模體內(nèi)不同深度5、10、20 cm 處的側(cè)向劑量分布曲線。兩種筆形束核算法與MC 方法的對(duì)比結(jié)果如表1所示，在不同射野面積大小下兩種算法的深度劑量的平均絕對(duì)誤差相差不大，算法的計(jì)算精度基本一致。

表2是傳統(tǒng)算法與改進(jìn)算法分別在水模體中不同射野大小下的劑量計(jì)算時(shí)間。兩種算法均在同一臺(tái)計(jì)算機(jī)上進(jìn)行，系統(tǒng)配置為：64 位，Windows 7，CPU：Intel（R）Core（TM）i5-3450。由表2可知，射野面積越大，計(jì)算時(shí)間越長(zhǎng)。與傳統(tǒng)算法相比，改進(jìn)算法的劑量計(jì)算所用時(shí)間大大減少。結(jié)合表1，兩種算法的計(jì)算精度基本一致，但改進(jìn)算法的劑量計(jì)算速度提升約2.7倍，大大減少了劑量計(jì)算的時(shí)間。

3.3 肺阻塊模體和骨阻塊模體結(jié)果

在肺阻塊模體中，我們使用傳統(tǒng)算法、改進(jìn)算法和MC方法分別計(jì)算射野大小為（10×10）cm2深度為10、16 cm處的側(cè)向劑量分布。在骨阻塊模體中，我們使用相同的算法分別計(jì)算射野大小為（10×10）cm2、深度為7.5和11.0 cm處的側(cè)向劑量分布。

圖2 3種算法在水模體中劑量計(jì)算結(jié)果Fig.2 Results of dose calculation in water phantom using 3 different algorithms

表1 兩種算法水在模體中不同射野大小下深度劑量與MC法的對(duì)比誤差（%）Tab.1 Errors of depth-dose curves of two methods compared with Monte Carlo(MC)algorithm in different field sizes of water phantom(%)

如圖3所示，圖3a 和圖3b 分別為肺阻塊模體中距離模體表層10和16 cm位置處的側(cè)向劑量曲線，圖3c 和圖3d 分別為骨阻塊模體中距離模體表層7.5 和11.0 cm位置處的側(cè)向劑量曲線。結(jié)合表3，在射野范圍內(nèi)，兩種算法在相同模體、相同深度處的側(cè)向劑量分布的誤差基本一致，并且都在允許范圍內(nèi)［17］。但是，由于電離不平衡現(xiàn)象，在兩種介質(zhì)的相交界面處和射野邊緣均出現(xiàn)了較大的誤差。

表4是傳統(tǒng)算法與改進(jìn)算法分別在兩種非均勻模體中不同深度處的側(cè)向劑量分布計(jì)算時(shí)間。兩種算法的計(jì)算環(huán)境與表2相同。與傳統(tǒng)算法相比，改進(jìn)算法大大減少劑量計(jì)算時(shí)間，同時(shí)改進(jìn)算法在這兩種模體中的劑量計(jì)算速度提升約2.6倍左右。

表2 水模體中兩種算法在不同射野大小下的劑量計(jì)算時(shí)間Tab.2 Time for dose calculation of two methods in different field sizes of water phantom

4 結(jié)論

圖3 3種算法在肺阻塊模體和骨阻塊模體中射野為（10×10）cm2的不同深度處的側(cè)向劑量分布Fig.3 Lateral dose distribution of 3 algorithms at different depths of field size of(10×10)cm2 of lung block phantom and bone block phantom

表3 兩種算法在不同模體不同深度處側(cè)向劑量分布與MC法對(duì)比誤差（%）Tab.3 Errors of lateral dose profiles of two algorithms compared with MC algorithm at different depths of different phantoms(%)

在本文中，采用以射束與模體的碰撞點(diǎn)為中心的方法計(jì)算劑量分布，在球殼坐標(biāo)系下利用粒子與模體碰撞后在不同球殼層散射方式的一致性，通過(guò)計(jì)算射束在初始球殼層碰撞點(diǎn)處軸線與體素的相交情況以及與深度有關(guān)的校正因子，得到其它層碰撞點(diǎn)處軸線與體素的相交情況，減少射線追蹤法的利用次數(shù)，降低劑量計(jì)算的復(fù)雜度，節(jié)省計(jì)算時(shí)間。

在不同模體內(nèi)計(jì)算劑量分布的實(shí)驗(yàn)結(jié)果表明，改進(jìn)算法與傳統(tǒng)算法相比，計(jì)算精度基本一致，在水模體和非均勻模體中的計(jì)算速度提升將近3倍，為劑量計(jì)算節(jié)省很多時(shí)間。同時(shí)，在球殼坐標(biāo)系下以碰撞點(diǎn)為中心計(jì)算劑量分布也為后續(xù)利用FPGA 和GPU對(duì)算法進(jìn)行加速的研究提供思路。

表4 兩種算法不同模體不同深度處的側(cè)向劑量曲線計(jì)算時(shí)間（s）Tab.4 Time for lateral dose calculation of two algorithms compared with MC algorithm at different depths of different phantoms(s)