周正東 陳元華 劉 娟
(南京航空航天大學(xué)核科學(xué)與工程系,南京 210016)
基于遺傳算法的醫(yī)用直線加速器光子能譜重建方法
周正東 陳元華 劉 娟
(南京航空航天大學(xué)核科學(xué)與工程系,南京 210016)
基于光子束中軸百分深度劑量(PDD),探討研究了基于遺傳算法的醫(yī)用直線加速器光能能譜精確重建方法.首先,利用蒙特卡洛模擬仿真醫(yī)用直線加速器治療頭,獲得6 MeV光子束的模擬能譜以及單能光子中軸PDD數(shù)據(jù);然后,根據(jù)測(cè)量得到的中軸PDD數(shù)據(jù)以及模擬得到的單能光子中軸PDD數(shù)據(jù),運(yùn)用遺傳算法優(yōu)化求解重建光子能譜.實(shí)驗(yàn)結(jié)果表明:重建能譜與蒙特卡洛模擬得到的能譜具有良好的一致性,平均相對(duì)誤差為3.03%;根據(jù)重建能譜計(jì)算得到的中軸PDD數(shù)據(jù)與測(cè)量得到的中軸PDD數(shù)據(jù)之間的平均相對(duì)誤差為1.0%,與蒙特卡洛模擬得到的中軸PDD數(shù)據(jù)之間的平均相對(duì)誤差為2.0%.由此可見,利用所提方法進(jìn)行光子束能譜重建可靠有效.
光子能譜;遺傳算法;百分深度劑量;蒙特卡羅模擬
腫瘤放射治療是利用電離輻射(放射線)治療腫瘤的一種方法.部分腫瘤病人只需要放療即可治愈;但對(duì)于有些腫瘤病人,無論是手術(shù)治療還是化療,為了降低風(fēng)險(xiǎn),在實(shí)施過程中往往還需要進(jìn)行輔助放療.隨著放療技術(shù)的發(fā)展,其臨床應(yīng)用日益廣泛,幾乎可用于全身各個(gè)部位的癌癥治療,成為癌癥治療的重要手段之一.
在放射治療中,提高腫瘤的控制率并降低正常組織的并發(fā)癥概率,一直是放療技術(shù)研究的熱點(diǎn).國際輻射單位度量委員會(huì)(ICRU)24號(hào)與42號(hào)報(bào)告[1-2]指出,臨床上要求腫瘤靶區(qū)的劑量不準(zhǔn)確度必須小于5%.而高精度劑量計(jì)算方法(如Col-lapsed cone劑量計(jì)算方法)需要根據(jù)醫(yī)用加速器治療頭所產(chǎn)生的射束能譜進(jìn)行計(jì)算,只有當(dāng)射束的能譜信息足夠精確時(shí),才能保證劑量計(jì)算的精確性.因此,獲得醫(yī)用加速器中X射線的能譜對(duì)于研究放療患者體內(nèi)的劑量分布非常重要.
直接測(cè)量加速器治療頭產(chǎn)生的光子束能譜是非常困難的.利用蒙特卡洛仿真軟件EGS進(jìn)行模擬計(jì)算,則需要準(zhǔn)確的加速器治療頭的結(jié)構(gòu)信息,該信息往往不易獲得,且模擬計(jì)算需要耗費(fèi)大量時(shí)間,因而難以應(yīng)用于實(shí)際放療過程中.在該領(lǐng)域中,研究者們更關(guān)注于利用容易獲得的劑量測(cè)量數(shù)據(jù),根據(jù)數(shù)學(xué)模型快速重建出精確的能譜.Deasy等[3]研究發(fā)現(xiàn)電子束百分深度劑量(percent depth dose,PDD)曲線與能譜之間存在很強(qiáng)的關(guān)聯(lián)性;Luo等[4]探索了由電子束PDD推導(dǎo)電子能譜的方法;Faddegon等[5]根據(jù)深度劑量分布來計(jì)算電子能譜;Deng等[6]根據(jù)PDD數(shù)據(jù),運(yùn)用隨機(jī)微變法對(duì)能譜進(jìn)行重建,但重建結(jié)果過度依賴于初始條件;張松柏等[7]采用衰減法對(duì)X射線能譜進(jìn)行重建,計(jì)算結(jié)果與EGS4模擬結(jié)果符合較好;姚杏紅等[8]根據(jù)PDD數(shù)據(jù),采用雙源模型對(duì)X射線能譜進(jìn)行計(jì)算,結(jié)果顯示,雙源模型更符合臨床實(shí)際情況.
根據(jù)PDD數(shù)據(jù)進(jìn)行能譜重建實(shí)際上是一個(gè)大規(guī)模優(yōu)化計(jì)算問題.相比于經(jīng)典的優(yōu)化方法,現(xiàn)代優(yōu)化方法(如遺傳算法)具有獨(dú)特的優(yōu)勢(shì).本文研究了基于遺傳算法的能譜重建方法,根據(jù)PDD數(shù)據(jù)求解重建醫(yī)用直線加速器光子束的能譜.結(jié)果表明,運(yùn)用遺傳算法可以得到全局最優(yōu)解,并可提高搜索效率,進(jìn)而提高能譜重建的精度與效率.
本文研究對(duì)象是西門子電子直線加速器中的6 MeV光子束.根據(jù)模擬所得的單能光子束PDD數(shù)據(jù)與實(shí)驗(yàn)測(cè)得的標(biāo)準(zhǔn)水體模中的PDD數(shù)據(jù),運(yùn)用遺傳算法進(jìn)行優(yōu)化計(jì)算,得到6 MeV光子束的能譜.
本文選用的射野大小為10 cm×10 cm.PDD的實(shí)際測(cè)量數(shù)據(jù)是在標(biāo)準(zhǔn)條件下利用電離室測(cè)量水箱中垂直于射野平面中心的一列體素得到的.至上而下共取n=380個(gè)體素,每個(gè)體素沿Z軸方向的厚度為1 mm,測(cè)量結(jié)果表示為行向量D.
將0~6 MeV平均分成24個(gè)區(qū)間,以各個(gè)區(qū)間的均值Ei(i=1,2,…,24)表示該區(qū)間的能量,即最終選取的單能為 0.125,0.375,…,5.875 MeV,共計(jì)24個(gè)單能光子束.利用蒙特卡羅仿真軟件EGS計(jì)算標(biāo)準(zhǔn)水體模中的單能光子束PDD數(shù)據(jù)[9].利用DOSXYZnrc子程序模擬光子輸運(yùn)過程[10],得到10 cm×10 cm方形射野中單能光子束在源皮距(SSD)為100 cm時(shí)產(chǎn)生的深度劑量,并以此作為能譜重建的基礎(chǔ)數(shù)據(jù).模擬計(jì)算時(shí),電子和光子的截止能量分別取0.7和0.01 MeV.模體大小設(shè)置為42 cm×42 cm×38 cm,介質(zhì)為水.垂直射束中心軸取長(zhǎng)2.0 cm、寬2.0 cm、中心點(diǎn)厚0.1 cm的一列體素作為劑量沉積體素.模擬歷程數(shù)設(shè)置為1×109,統(tǒng)計(jì)誤差小于0.1%.
將模擬得到的24組單能光子PDD數(shù)據(jù)用矩陣A來表示,以便于利用Matlab軟件強(qiáng)大的矩陣運(yùn)算能力進(jìn)行快速運(yùn)算.Aij(i=1,2,…,24;j=1,2,…,380)表示能量為Ei的單能光子在深度j處的PDD數(shù)據(jù).
加速器治療頭所產(chǎn)生的光子并非單能,在通過源皮距100 cm的過程中,會(huì)與不同介質(zhì)發(fā)生各種相互作用,最終射束中光子能量會(huì)呈現(xiàn)出一種特定的連續(xù)譜分布.該射束是由所有單能成分按各自的權(quán)重疊加而成的.因此,測(cè)量的PDD數(shù)據(jù)是所有單能光子PDD數(shù)據(jù)按上述權(quán)重疊加的結(jié)果,即
式中,f(E)dE為能量E~E+ΔE的X射線中光子所占射束全部能量范圍內(nèi)光子數(shù)的百分?jǐn)?shù);P(E)為該能量區(qū)間內(nèi)光子的PDD數(shù)據(jù);Pm為測(cè)量的PDD數(shù)據(jù).
實(shí)際中不可能通過模擬得到所有單能光子的PDD數(shù)據(jù).故將0~6 MeV平均分成24個(gè)能量區(qū)間,以各個(gè)區(qū)間的均值代替該區(qū)間的能量,簡(jiǎn)化為24個(gè)單能成分,以求和代替積分來計(jì)算光子能譜,則
式中,W={w1,w2,…,w24}為光子能譜,是一個(gè)含有24個(gè)元素的行向量,其中的每個(gè)元素代表對(duì)應(yīng)單能成分的權(quán)重為單能光子PDD 矩陣,共計(jì)24行380列;D={P1,P2,…,P380}為測(cè)量百分深度劑量,是一個(gè)含有380個(gè)元素的行向量.由此可知,W×A表示線性疊加后的PDD數(shù)據(jù),是一個(gè)含有380個(gè)元素的行向量,與D同維.
運(yùn)用某種數(shù)學(xué)方法求解式(2),便可得到光子束的能譜.
實(shí)際劑量測(cè)量數(shù)據(jù)中存在誤差及計(jì)算的離散化,故式(2)屬于病態(tài)線性問題,微小的誤差將會(huì)被放大,使得反演結(jié)果遠(yuǎn)遠(yuǎn)偏離問題的真實(shí)解.該問題通常采用迭代優(yōu)化的方法進(jìn)行求解.
本文運(yùn)用遺傳算法對(duì)光子束能譜進(jìn)行優(yōu)化重建.遺傳算法是通過模擬生物在自然環(huán)境中的遺傳和進(jìn)化過程而形成的一種自適應(yīng)全局優(yōu)化概率搜索方法,其基本特征是:在種群的不斷演化過程(即求解過程)中,通過借鑒生物界“物競(jìng)天擇,適者生存”的原則,找到滿意或最優(yōu)的解.該算法一般包括3個(gè)基本操作:選擇、雜交和變異.選擇起著向?qū)У淖饔?,使搜索朝著搜索空間的最優(yōu)區(qū)域內(nèi)進(jìn)行,充分利用群體內(nèi)當(dāng)前所具有的有效信息,將搜索重點(diǎn)放在具有較高適應(yīng)度值的個(gè)體上,向最優(yōu)的方向進(jìn)行演化;雜交對(duì)群體內(nèi)出現(xiàn)的現(xiàn)有信息進(jìn)行重組,發(fā)現(xiàn)與環(huán)境更為適應(yīng)的個(gè)體,尋找可能的最優(yōu)區(qū)域;變異給群體帶來新的遺傳基因,恢復(fù)由于選擇而失去的個(gè)體多樣性.遺傳算法利用了概率搜索技術(shù),進(jìn)行解空間的多點(diǎn)搜索,具有極強(qiáng)的容錯(cuò)能力,易獲得全局最優(yōu)解.
針對(duì)能譜重建問題,定義如下的適應(yīng)度函數(shù):
由于劑量測(cè)量存在誤差,計(jì)算過程以求和代替積分,W×A會(huì)落在D附近一個(gè)很小的范圍內(nèi),無法滿足W×A=D嚴(yán)格的等式關(guān)系.為防止問題的解變成病態(tài)解或者無法收斂,引入一個(gè)合適的容差因子ε,構(gòu)成如下不等式:
當(dāng)相鄰數(shù)代的解均滿足式(4)時(shí),迭代終止.
利用遺傳算法進(jìn)行能譜重建的詳細(xì)步驟如下:
①初始化.隨機(jī)選取初始化個(gè)體數(shù)目為40的群體,群體內(nèi)每個(gè)個(gè)體以實(shí)數(shù)編碼形式進(jìn)行表示,產(chǎn)生一系列的個(gè)體其中q=1,2,…,n,k=1,2,…,40,n為最優(yōu)問題的變量個(gè)數(shù),k為個(gè)體編號(hào),[aq,bq]為問題解向量的取值范圍.
②適應(yīng)度計(jì)算.利用式(3)計(jì)算種群中每個(gè)個(gè)體的適應(yīng)度值,以此表示進(jìn)化的好壞程度,為下一步選擇操作提供數(shù)據(jù)和方向.
③選擇.根據(jù)個(gè)體的適應(yīng)度函數(shù)值,從父代中選擇出2個(gè)個(gè)體;適應(yīng)度函數(shù)值越大,則被選中的概率越大.
④ 交叉.將選擇出的2個(gè)個(gè)體以0.8的雜交概率進(jìn)行雜交,產(chǎn)生2個(gè)新的個(gè)體.
⑤ 變異.對(duì)新產(chǎn)生的2個(gè)個(gè)體按0.1的變異概率進(jìn)行變異.
⑥再生.接受新個(gè)體,并判斷是否完成新群體的生成.如果沒有,則返回步驟③.
⑦演化終止條件檢驗(yàn).若滿足不等式(4),則停止演化;否則,轉(zhuǎn)步驟②,重新進(jìn)行演化.每一次進(jìn)化過程均產(chǎn)生新一代的群體.當(dāng)演化停止時(shí),群體內(nèi)個(gè)體表示的解即是所求的最優(yōu)解.
將基于遺傳算法計(jì)算得到的能譜與蒙特卡洛模擬得到的能譜進(jìn)行比較,結(jié)果如圖1所示.圖中,Ss表示運(yùn)用EGS模擬加速器治療頭經(jīng)Beamdp程序分析獲得的模擬能譜;Sr表示根據(jù)本文方法優(yōu)化計(jì)算獲得的重建能譜.由圖可知,2條曲線的峰位均是0.625 MeV,峰值相對(duì)誤差為2.9%,平均相對(duì)誤差為3.03%,半高寬基本一致,即這2條曲線的峰值、峰位和半高寬均吻合良好.
圖1 能譜比較
根據(jù)上述重建能譜Sr及單能光子PDD數(shù)據(jù),可計(jì)算得到射束的PDD數(shù)據(jù).將重建能譜作為EGS模擬子程序DOSXYZnrc的輸入源,進(jìn)行EGS模擬,即可得到模擬的PDD數(shù)據(jù).將上述2個(gè)結(jié)果及測(cè)量的PDD數(shù)據(jù)進(jìn)行比較,結(jié)果如圖2所示.圖中,Pm為實(shí)際測(cè)量的PDD數(shù)據(jù),Pc為根據(jù)重建能譜計(jì)算得到的PDD數(shù)據(jù),Ps為以重建能譜作為DOSXYZnrc輸入源進(jìn)行模擬得到的PDD數(shù)據(jù).由圖可知,Pm和Pc的峰位相差2 mm,平均相對(duì)誤差為1.0%;Pc和Ps的峰位相差1 mm,平均相對(duì)誤差為2.0%.根據(jù)重建能譜計(jì)算得到的PDD數(shù)據(jù)與測(cè)量值及模擬值相比,除表層和最底層少數(shù)體素的相對(duì)誤差略大外,其余各點(diǎn)相對(duì)誤差均很小,半高寬基本一致.
圖2 PDD數(shù)據(jù)的比較
綜上所述,采用遺傳算法重建能譜,對(duì)初始化條件依賴小,能夠獲得全局最優(yōu)解,計(jì)算結(jié)果可靠精確.
本文根據(jù)測(cè)量得到的PDD數(shù)據(jù)以及蒙特卡羅模擬獲得的PDD數(shù)據(jù),運(yùn)用遺傳算法重建醫(yī)用直線加速器能譜,從能譜及PDD兩個(gè)方面對(duì)重建的結(jié)果進(jìn)行了評(píng)價(jià).實(shí)驗(yàn)結(jié)果表明,該方法可靠有效,具有臨床應(yīng)用價(jià)值.遺傳算法適用于處理多變量問題,在搜索空間中同時(shí)處理種群中的多個(gè)個(gè)體,可提高搜索效率,并通過設(shè)置合適的演化參數(shù),可避免陷入局部最優(yōu)解.由于遺傳算法隱含并行性,因此可以采用GPU并行計(jì)算來進(jìn)一步提高能譜重建的效率.
[1] International Commission on Radiation Units and Measurements.Determination of absorbed dose in patient irradiated by beams of X or gamma rays in radiotherapy procedures,ICRU Report No.24[R].Washington DC:ICRU,1976.
[2] International Commission on Radiation Units and Measurements.Use of computers in external beam radiotherapy procedure with high-energy photons and electrons,ICRU Report No.42[R].Washington DC:ICRU,1987.
[3] Deasy J O,Almond P R.The spectral dependence of electron central-axis depth-dose curves[J].Med Phys,1994,21(9):1369-1376.
[4]Luo Z,Jette D.On the possibility of determining an effective energy spectrum of clinical electron beams from percentage depth dose(PDD)data of broad beams[J].Phys Med Biol,1999,44(8):177-182.
[5] Faddegon B A,Blevis I.Electron spectra derived from depth dose distributions[J].Med Phys,2000,27(3):514-526.
[6] Deng J,Jiang S B,Pawlicki T,et al.Derivation of electron and photon energy spectra from electron beam central axis depth dose curves[J].Phys Med Biol,2001,46(5):1429-1449.
[7]張松柏,黃斐增,韓樹奎,等.6MV醫(yī)用電子直線加速器軔致輻射譜研究[J].原子能科學(xué)技術(shù),2003,37(3):206-210.
Zhang Songbai,Huang Feizeng,Han Shukui,et al.Investigation of bremsstrahlung spectra of 6MV medical linear accelerator[J].Atomic Energy Science and Technology,2003,37(3):206-210.(in Chinese)
[8]姚杏紅,黃斐增,吳昊,等.利用深度劑量數(shù)據(jù)重建放射治療X射線能譜[J].原子能科學(xué)技術(shù),2010,44(4):464-468.
Yao Xinghong,Huang Feizeng,Wu Hao,et al.Reconstruction of X-ray spectra for radiation therapy[J].Atomic Energy Science and Technology,2010,44(4):464-468.(in Chinese)
[9]王妙琪.利用蒙地卡羅法模擬西門子醫(yī)用直線加速器輸出之光子射束特性[D].臺(tái)北:陽明大學(xué)生物醫(yī)學(xué)影像暨放射科學(xué)系,2009.
[10] Ionizing Radiation Standards National Research Council of Canada.DOSXYZnrc users manual,NRCC Report No.PIRS-794 rev B[R].Ottawa,Canada:NRCC,2009.
Reconstruction method of photon spectra of medical linear accelerator based on genetic algorithm
Zhou Zhengdong Chen Yuanhua Liu Juan
(Department of Nuclear Science and Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China)
A method for precise reconstruction of photon spectra of medical linear accelerator based on genetic algorithm was investigated according to percentage depth dose(PDD)data of photon beam central axis.First,the spectra of 6 MeV photon beams and PDD data of mono-energetic photon beams were calculated by Monte Carlo simulation of accelerator treatment head.Then,the genetic algorithm was employed to reconstruct photon spectra according to both measured and simulated PDD data.The experimental results show that the reconstructed spectra is in good agreement with that calculated by Monte Carlo simulation,the average relative error being 3.03%.The derived PDD data calculated from the reconstructed spectra agree well with both measured and Monte Carlo simulated PDD data,and the average relative errors are 1.0%and 2.0%,respectively.Therefore,the proposed method is reliable and effective for photon spectra reconstruction.
photon spectra;genetic algorithm;percent depth dose;Monte Carlo simulation
R318;TP391
A
1001-0505(2012)06-1085-04
10.3969/j.issn.1001 -0505.2012.06.012
2012-07-10.
周正東(1969—),男,博士,副教授,zzd_msc@nuaa.edu.cn.
國家自然科學(xué)基金資助項(xiàng)目(20975052).
周正東,陳元華,劉娟.基于遺傳算法的醫(yī)用直線加速器光子能譜重建方法[J].東南大學(xué)學(xué)報(bào):自然科學(xué)版,2012,42(6):1085-1088.[doi:10.3969/j.issn.1001 -0505.2012.06.012]