郭金森GUO Jin-sen；李傳龍LI Chuan-long；王薇WANG Wei

（中國原子能科學(xué)研究院，北京102413）

0 引言

在利用加速器打靶的韌致輻射光子束進(jìn)行治療時，劑量計(jì)算軟件為加速器TPS 治療計(jì)劃系統(tǒng)的基礎(chǔ)。任何劑量計(jì)算程序都需要用到光子能譜數(shù)據(jù)。然而在實(shí)際的治療系統(tǒng)中，有些情況下很難直接測量能譜數(shù)據(jù)。因此針對醫(yī)用加速器，不同的研究者一直在探求獲取準(zhǔn)確能譜的方法[1]。主流方法有：一是利用蒙特拉羅方法模擬加速器治療機(jī)頭得到能譜，如Deng J 等人[2]在2000 年利用EGS4 程序模擬了來自瓦里安Clinac2100c 和2300c/d 加速器的4、6和15MV 光子束能譜。2015 年Fátima Padilla-Cabal 等人[3]利用MCNPX 和EGS 程序建立了醫(yī)科達(dá)加速器的精確模型并做了劑量計(jì)算。這種方法能夠保證所建能譜的精度，但前提是需要精確地知道準(zhǔn)直器、均整器、電離室等部件的尺寸、材料等數(shù)據(jù)。二是利用射束穿過不同介質(zhì)的透射數(shù)據(jù)得到X 射線能譜，典型如[4，5]等人的工作，該方法并不適用于高能光子射束。三是在射束中軸測量PDD 數(shù)據(jù)得到精確測量的前提下[6]，利用數(shù)學(xué)方法進(jìn)行能譜重建。如2001 年，Deng J[7]等人在研究了從電子束中心軸百分比深度劑量（PDD）曲線，采取“隨機(jī)蠕變”（random creep）算法推導(dǎo)了臨床電子和光子能譜，但該方法有可能收斂至局部最優(yōu)解。2010 年北京大學(xué)姚杏紅[8]等人采用Cimmino 方法重建瓦里安600CD 直線加速器X 射線能譜，該方法基于雙源模型，結(jié)果比較精確，但大射野情況下重建能譜出現(xiàn)了病態(tài)解。2013 年陳元華[9]等人針對瓦里安15MV 光子束利用遺傳算法優(yōu)化重建了光子能譜，得到的能譜與蒙特卡洛模擬得到的能譜具有較好的一致性。2013 年劉娟[10]等人根據(jù)利用模擬退火算法重建了西門子醫(yī)用加速器6MV X射線能譜。

模擬退火算法為化工、冶金工業(yè)術(shù)語，現(xiàn)在已經(jīng)作為最優(yōu)化算法用于各個行業(yè)。模擬退火算法的核心是在爬山算法的基礎(chǔ)上增加Metropolis 準(zhǔn)則，即從當(dāng)前狀態(tài)i→新狀態(tài)j 時，即使新狀態(tài)j 在某個規(guī)則下優(yōu)于狀態(tài)i，也以一定的概率接受狀態(tài)j，否則仍保留狀態(tài)i。模擬退火算法有一定的概率跳出局部最優(yōu)解從而找到全局最優(yōu)解。本文擬利用模擬退火算法，結(jié)合實(shí)驗(yàn)PDD 數(shù)據(jù)，對國內(nèi)某型醫(yī)用加速器6MV、18MV 能量下X 射線能譜進(jìn)行重建，重建結(jié)果和蒙特卡羅模擬得到的能譜進(jìn)行對比以驗(yàn)證可靠性。并為后續(xù)劑量計(jì)算程序開發(fā)奠定基礎(chǔ)。

1 材料與方法

1.1 問題的數(shù)學(xué)描述

由測量PDD 和單能光子PDD 求解光子能譜的問題，可以描述為解線性方程組問題。假設(shè)加速器出射能譜可以離散化為E1，E2…En，其中E1~En為每個離散化能量區(qū)間的平均能量。n 組能量的按照注量加權(quán)的權(quán)重為ω1，ω2，…ωn，則應(yīng)有如下如下線性方程組成立：

其中dij代表深度為i 處，Ei能量下的單能PDD 的劑量數(shù)據(jù)。Dj為深度為Hj處不同權(quán)重下的單能PDD 加權(quán)相加后的合成PDD 數(shù)據(jù)，在ω1，ω2，…ωn為實(shí)際能譜情況下，D1…Dm即為測量PDD 數(shù)據(jù)。該方程組一般為超定方程組，其解不止一個。因此需要找到最符合物理意義的一組解，即為能譜數(shù)據(jù)?？梢姡顑?yōu)化的能譜可以使合成PDD 數(shù)據(jù)最接近測量PDD，也即相似度最高。為描述合成PDD（用量F 表示）和測量PDD（用M 表示）x 相似度，引入相關(guān)系數(shù)：

進(jìn)行能譜重建。其中f（x）為目標(biāo)函數(shù)，F(xiàn)i為深度為Hi處的合成PDD，Mi為深度為Hi處的測量PDD。

1.2 單能光子PDD 以及能譜模擬

為獲?。?）中的單能PDD，采用蒙特卡羅模擬程序BEAMnrc 的劑量計(jì)算程序xyznrc 進(jìn)行模擬。模體為30cm×30cm×30cm 厚的水模，密度為1.0g/cm3。源皮距為100cm，射野為5cm×5cm。沿射野中心軸取1.5cm×1.5cm×0.1cm 體素共299 個。針對6MV 計(jì)算了17 個能量點(diǎn)的單能PDD 數(shù)據(jù)，針對18MV 則計(jì)算了23 個能量點(diǎn)，即公式（1）中n=23。所有模擬結(jié)果誤差最大均不超過0.4%。

為驗(yàn)證重建能譜的準(zhǔn)確性，利用蒙特卡羅模擬程序EGSnrc 建立了加速器治療頭模型，獲取了治療頭6MV 和18MV 的出射能譜。加速器治療頭結(jié)構(gòu)包括靶、初級準(zhǔn)直器、均整器、電離室，次級準(zhǔn)直器（上、下兩個）等結(jié)構(gòu)。通過控制上下次級準(zhǔn)直器控制射野大小為5cm×5cm，分別模擬了6MeV 與18MeV 窄電子束入射情況下的韌致輻射譜。見圖2，圖3。

1.3 利用模擬退火算法重建光子能譜

在熱力學(xué)上，一塊被加熱至高溫的物體的降溫過程被稱之為“退火”。退火過程滿足Metropolis 準(zhǔn)則，即溫度為T時，出現(xiàn)能量差為dE 的降溫概率P（E）為：

即溫度越高，出現(xiàn)一次能量差為dE 的降溫的概率就越大；溫度越低，相應(yīng)概率就越小。

模擬退火算法的基本流程圖見圖1。

圖1 模擬退火算法流程圖

模擬退火算法基本流程為：

①對光子各能量箱賦予一定的權(quán)重，經(jīng)驗(yàn)表明該權(quán)重不可過分偏離實(shí)際權(quán)重，否則可能導(dǎo)致收斂過慢甚至收斂至病態(tài)解。

②開始迭代過程，每一次迭代在每個能量箱上加一個小的隨機(jī)量，即

其中，η 為攝動系數(shù)，取0.00001，εi取[-1，1]之間的隨機(jī)數(shù)。

③如果新的目標(biāo)函數(shù)f（x）（見公式（3））小與當(dāng)前目標(biāo)函數(shù)，則接受新解和新的目標(biāo)函數(shù)。如果新目標(biāo)函數(shù)大于當(dāng)前目標(biāo)函數(shù)，則以概率P（C）=e-dC/T接受新解和新的目標(biāo)函數(shù)。容易看出隨著T 逐漸降低，越難接受“壞”的新解。迭代終止條件為目標(biāo)函數(shù)低于某一個截止值，或模擬退火溫度低于某個截止值。

模擬退火降溫過程由初始溫度T 及溫度控制參數(shù)a表示，如下

這里初始溫度取1×10-9，a 取0.999，每個溫度出迭代100 次，取目標(biāo)函數(shù)截止值Scut為4.6×10-5，即合成PDD 和測量PDD 之間的相關(guān)系數(shù)C=1-Scut即0.99954。

2 結(jié)果

圖2、圖3 分別為優(yōu)化前后6MV、18MV 能譜和MC 模擬能譜對比圖。二者均作了歸一化處理?？梢钥闯鰞?yōu)化能譜和MC 模擬能譜峰位相同，譜形基本一致。計(jì)算得到相關(guān)系數(shù)均為0.99 以上。另外可以看出，圖2 優(yōu)化后能譜尾端有微小抬升，懷疑此處算法陷入局部最優(yōu)解。但低能光子由于權(quán)重小對總劑量貢獻(xiàn)較低。而高能部分兩條能譜基本一致。圖3 優(yōu)化后能譜和MC 模擬能譜除個別點(diǎn)外也基本一致。

圖2 優(yōu)化后6MV 能譜和MC 模擬能譜對比圖

圖3 優(yōu)化后18MV 能譜和MC 模擬能譜對比圖

圖4-圖5 為優(yōu)化前后6MV、18MV PDD 數(shù)據(jù)相對于測量數(shù)據(jù)的偏差。可以看出優(yōu)化后絕大部分點(diǎn)相對誤差在0.5%以下。

圖4 優(yōu)化前后后6MV PDD 相對測量數(shù)據(jù)的偏差

圖5 優(yōu)化前后18MV PDD 相對測量數(shù)據(jù)的偏差

3 討論

本研究利用模擬退火算法，基于蒙特卡羅模擬的單能光子PDD 數(shù)據(jù)和測量得到的PDD，對醫(yī)用電子加速器6MV、18MV X 射線能譜進(jìn)行重建。計(jì)算得到的能譜與蒙卡程序直接模擬治療頭得到的能譜形狀基本一致。

模擬退火算法具有全局搜索性，相比簡單的爬山算法，有更大的概率得到全局最優(yōu)解。需要指出的是，選擇合適的初值仍然是有必要的，首先可以減小搜索時間，其次進(jìn)一步降低收斂至局部最優(yōu)解的風(fēng)險。今后的工作將分為兩步進(jìn)行：①調(diào)整算法，嘗試在迭代計(jì)算初期以更大的步長進(jìn)行搜索，在目標(biāo)函數(shù)f（x）滿足一定條件時進(jìn)行小步長精細(xì)化搜索，提高算法魯棒性；②利用建成能譜開發(fā)卷積劑量計(jì)算程序。