楊 婕，桂志國

（1.中北大學(xué) 電子測試技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，山西 太原030051；2.山西中醫(yī)學(xué)院 醫(yī)藥管理學(xué)院，山西 太原030619）

0 引 言

調(diào)強(qiáng)放射治療（Intensity Modulation Radia-tion Therapy，IMRT）被譽(yù)為21世紀(jì)放射治療發(fā)展的方向，是目前治療惡性腫瘤的主要手段之一.它通過調(diào)節(jié)照射腫瘤靶區(qū)（Planning Target Vol-ume，PTV）射線束的強(qiáng)度，使得靶區(qū)接受的劑量最大，靶區(qū)周圍危及器官（Organ At Risk，OAR）和正常組織（Normal Tissues，NT）接受的劑量最小，從而提高治療增益比，改善病人生存質(zhì)量.逆向計(jì)劃是IMRT的關(guān)鍵環(huán)節(jié)，其主要過程是依據(jù)醫(yī)生給出PTV處方劑量和OAR耐受劑量，綜合考慮各種因素確定優(yōu)化目標(biāo)函數(shù)，再由計(jì)算機(jī)通過各種優(yōu)化算法計(jì)算每個(gè)射野的最佳射束強(qiáng)度分布，最終得到較優(yōu)的治療計(jì)劃.逆向計(jì)劃的質(zhì)量直接影響到治療的精度和療效.

IMRT逆向計(jì)劃中優(yōu)化目標(biāo)函數(shù)分為物理目標(biāo)函數(shù)和生物目標(biāo)函數(shù).物理目標(biāo)函數(shù)描述的是處方劑量與所計(jì)算劑量分布之間的差值，如平均劑量函數(shù)、最小劑量函數(shù)、最大劑量函數(shù)和劑量-體積（Dose-Volume，DV）準(zhǔn)則等［1］.生物目標(biāo)函數(shù)描述的是方案治療效果，是腫瘤或正常組織內(nèi)的非線性放射生物效應(yīng)，如等效均一劑量（Equivalent Uniform Dose，EUD）、腫瘤控制率（Tumor Control Probability，TCP）和正常組織并發(fā)癥發(fā)生概率（Normal Tissue Complication Probability，NTCP）.物理目標(biāo)函數(shù)的具有定義簡單，便于應(yīng)用的優(yōu)點(diǎn)，但具有一定的局限性，不能準(zhǔn)確地預(yù)測腫瘤或正常組織內(nèi)的非線性放射生物效應(yīng)，生物目標(biāo)函數(shù)已成為一種替代物理目標(biāo)函數(shù)的選擇［2］.有文獻(xiàn)［3-5］表明，生物目標(biāo)函數(shù)可以得到更好的劑量分布.

IMRT逆向計(jì)劃優(yōu)化方法分為單目標(biāo)法和多目標(biāo)法.單目標(biāo)法是將多個(gè)優(yōu)化目標(biāo)通過加權(quán)求和化為單一函數(shù)再使用成熟算法優(yōu)化，如L-BFGS（Limited Memory BFGS）算法.該方法的優(yōu)點(diǎn)是收斂速度較快，優(yōu)化時(shí)間短；缺點(diǎn)是對于非凸函數(shù)得到的可能是局部極值，各目標(biāo)函數(shù)的權(quán)重賦值是個(gè)反復(fù)試誤（Trial and Error）的過程［6］，每次計(jì)算只能得到一個(gè)優(yōu)化解.多目標(biāo)法是采用進(jìn)化算法直接進(jìn)行優(yōu)化，如NSGA-Ⅱ（Non-Dominated Sorting Genetic Algorithm-Ⅱ）算法.其優(yōu)點(diǎn)是一次計(jì)算就可以得到Pareto最優(yōu)解集合，具有更好的靈活性，缺點(diǎn)是響應(yīng)時(shí)間長，收斂速度慢等.

針對上述情況，本文將物理-生物混合目標(biāo)函數(shù)引入到IMRT逆向計(jì)劃中，利用L-BFGS算法對NSGA-Ⅱ算法進(jìn)行加速，并對NSGA-Ⅱ算法的快速非支配排序策略進(jìn)行了修改.將本文算法與單一物理目標(biāo)和其他優(yōu)化方法在10名前列腺腫瘤患者中進(jìn)行了對比實(shí)驗(yàn).實(shí)驗(yàn)結(jié)果表明，利用新算法所得放射治療計(jì)劃質(zhì)量更高，所需計(jì)算時(shí)間更短.

1 本文方法

逆向計(jì)劃優(yōu)化的目的是得到最佳的放療方案，在對靶區(qū)產(chǎn)生不可恢復(fù)性摧毀的同時(shí)保護(hù)周圍關(guān)鍵組織和正常組織.為了得到更好的放療方案，本文對常見優(yōu)化算法［7-11］進(jìn)行了深入研究，在此基礎(chǔ)提出一種新算法.首先確定物理-生物的混合目標(biāo)函數(shù)，以物理函數(shù)最小-最大劑量函數(shù)控制靶區(qū)的劑量分布，以物理函數(shù)DV準(zhǔn)則和生物函數(shù)廣義EUD（generalized EUD，gEUD）［12］控制危及器官的劑量分布.然后僅選擇靶區(qū)和危及器官的平均劑量函數(shù)作為優(yōu)化目標(biāo)簡化混合優(yōu)化目標(biāo)，利用L-BFGS對其優(yōu)化.將L-BFGS優(yōu)化所得到的解與隨機(jī)生成的值合并作為NSGA-Ⅱ的第一代種群.最后利用NSGA-Ⅱ優(yōu)化混合目標(biāo)函數(shù)，得到多個(gè)劑量分布方案供臨床使用.

1.1 混合優(yōu)化目標(biāo)函數(shù)

通常逆向計(jì)劃優(yōu)化采用的是物理函數(shù)目標(biāo)，計(jì)算每個(gè)體素劑量值與處方劑量差的平方，但其不能準(zhǔn)確反映器官組織受到不均勻劑量照射時(shí)的生物效應(yīng).因此為了保證放射治療計(jì)劃的質(zhì)量，本文在逆向計(jì)劃中提出一種帶約束的物理-生物混合目標(biāo)函數(shù)：使用最小-最大劑量函數(shù)fmin-max（D）作為約束條件控制PTV的劑量，確保靶區(qū)的高劑量分布；利用DV準(zhǔn)則函數(shù)fDV（V50）作為約束條件控制危及器官內(nèi)劑量達(dá)到50 Gy的體積；利用gEUD函數(shù)fgEUD（D）和DV準(zhǔn)則函數(shù)fDV（V75）作為目標(biāo)函數(shù)，盡量降低危及器官內(nèi)平均劑量和高劑量的分布.其數(shù)學(xué)模型如下

式中：Di為第i個(gè)體素的劑量值；N為器官組織內(nèi)的所有體素點(diǎn)；α（-∞≤α≤+∞）為劑量體積效應(yīng)因子，靶區(qū)組織的α≤1，正常組織的α＞1；TD5／5是危及器官最小耐受劑量；V75指器官組織內(nèi)接收到劑量75 Gy的體積；Dmean為靶區(qū)處方劑量或正常組織的平均耐受劑量；Cd是適當(dāng)放寬靶區(qū)劑量的指數(shù)；V50指器官組織內(nèi)接收到劑量50 Gy的體積；Vmax%指DV約束中對應(yīng)的體積上限.

由于DV準(zhǔn)則函數(shù)和gEUD函數(shù)不是凸函數(shù)，求導(dǎo)過程較復(fù)雜，為了便于后期使用L-BFGS算法優(yōu)化，本文將上述混合模型簡化，將其改為只含有平均劑量函數(shù)的無約束優(yōu)化目標(biāo)函數(shù).數(shù)學(xué)模型如下

1.2 混合優(yōu)化算法

逆向計(jì)劃優(yōu)化常見算法分為確定型和概率型，其中具有代表性的是L-BFGS算法和NSGA-Ⅱ算法.L-BFGS算法是BFGS算法的改進(jìn)，是擬牛頓法優(yōu)化算法的一種，屬于確定型算法.其采用目標(biāo)函數(shù)的一階梯度信息構(gòu)造一系列的正定矩陣逼近大規(guī)模的Hessian矩陣，從而在計(jì)算復(fù)雜度和存儲(chǔ)空間都有明顯降低，適合求解大規(guī)模的優(yōu)化問題.NSGA-Ⅱ算法是NSGA的改進(jìn)，屬于進(jìn)化算法的一種，屬于概率型算法.其通過增加精英保留策略、密度估計(jì)策略和快速非支配排序策略，改善了NSGA算法的缺點(diǎn).它以在解決多目標(biāo)優(yōu)化問題時(shí)的優(yōu)良性得到人們的廣泛關(guān)注.

優(yōu)化目標(biāo)函數(shù)時(shí)，如果單一使用L-BFGS算法首要需要解決的問題就是對目標(biāo)函數(shù)求導(dǎo)，而逆向計(jì)劃優(yōu)化中大部分函數(shù)都是非凸的，其求導(dǎo)過程較復(fù)雜；如果單一使用NSGA-Ⅱ算法直接優(yōu)化，又存在由于變量維數(shù)過高得不到滿意解和響應(yīng)時(shí)間過長的問題.因此，本文將兩種算法進(jìn)行混合，先以L-BFGS算法對簡化的凸函數(shù)優(yōu)化目標(biāo)函數(shù)進(jìn)行優(yōu)化，再以NSGA-Ⅱ算法對全部目標(biāo)函數(shù)進(jìn)行優(yōu)化.主要流程如下

Stepl：輸入CT數(shù)據(jù)；

Step2：判斷是否運(yùn)行L-BFGS算法5次，滿足則跳轉(zhuǎn)到Step5，否則執(zhí)行Step3；

Step3：確定優(yōu)化函數(shù)；

Step4：運(yùn)行L-BFGS算法優(yōu)化目標(biāo)函數(shù)f，得到優(yōu)化解，并保存至數(shù)組X；

Step5：確 定 多 目 標(biāo) 優(yōu) 化 函 數(shù)，fgEUD（DOAR）和fDV（VOAR75），約束條件為fmin-max（DPTV）和fDV（VOAR50）；

Step6：隨機(jī)產(chǎn)生NSGA-Ⅱ算法295個(gè)初始解，并與X合并；

Step7：運(yùn)行NSGA-Ⅱ算法優(yōu)化目標(biāo)；

Step8：得到Pareto解集，結(jié)束.

本文需采用L-BFGS算法優(yōu)化的目標(biāo)函數(shù)只有平均劑量函數(shù)，給出其一階導(dǎo)數(shù).

本文采用NSGA-Ⅱ算法優(yōu)化的混合目標(biāo)，將fmin-max（DPTV和fDV（VOAR50）作為約束條件.

逆向計(jì)劃優(yōu)化過程屬于帶約束的多目標(biāo)優(yōu)化問題.NSGA-Ⅱ算法在處理約束條件時(shí)是將約束條件改成目標(biāo)函數(shù)進(jìn)行優(yōu)化，得到優(yōu)化解后再通過決策挑選出適用解，這樣無疑加大了優(yōu)化的工作量.因此本文對NSGA-Ⅱ算法的快速非支配排序策略進(jìn)行了改動(dòng)，假設(shè)a和b是兩個(gè)解，在判斷其是否為非劣解時(shí)首先要求滿足約束條件，再判斷目標(biāo)函數(shù)最小值.具體改動(dòng)為

2 實(shí)驗(yàn)及結(jié)果分析

為了驗(yàn)證算法有效性，本文對10例前列腺腫瘤病人的實(shí)際數(shù)據(jù)，分別使用混合目標(biāo)混合優(yōu)化法（M-M）、物理目標(biāo)混合優(yōu)化法（S-M）、混合目標(biāo)單一NSGA-Ⅱ法（M-S）和混合目標(biāo)無約束混合法（M-NM）進(jìn)行了對比實(shí)驗(yàn)，并利用DVH圖、TCP［13］、NTCP［13］和HI［14］指數(shù)對不同方案優(yōu)化算法的結(jié)果進(jìn)行了評價(jià).

器官組織CT（見圖1）由患者仰臥位采取.PTV為前列腺，OAR1是直腸（Rectum），OAR2是膀胱（Bladder）.

圖1 盆腔器官組織CT Fig.1 CT of pelvic organs

CERR開源軟件采用5個(gè)6 MeV共弧面照射靶目標(biāo)，機(jī)架角度依次為：36°，100°，180°，260°，324°，勾畫PTV和OAR得到劑量效應(yīng)矩陣.

L-BFGS優(yōu)化目標(biāo)函數(shù)為f（D）=ζ1·fmean（DOAR2）+ζ2fmean（DOAR1+ζ3fmean（DPTV），靶區(qū)的均勻處方劑量為78 Gy，危及器官膀胱和直腸的耐受劑量設(shè)定為最小耐受劑量TD5／5：60 Gy，60 Gy.ζi是服從［0，1］的隨機(jī)數(shù)，并做歸一化處理.

NSGA-Ⅱ優(yōu) 化 目 標(biāo) 函 數(shù) 為fgEUD（DOAR2），fgEUD（DOAR1），fDV（VOAR275）和fDV（VOAR175）.約束條件為fmin-max（DPTV），fDV（VOAR250）和fDV（VOAR150）.其中靶區(qū)處方劑量為78 Gy，危及器官受到的g EUD為60 Gy，gEUD模型中直腸和膀胱α取值分別為8和8［15-16］.約束條件是VOAR150%=50%，VOAR150%=50%，Cd=2.

本文對比實(shí)驗(yàn)所使用的方案均在MATLAB平臺上編程實(shí)現(xiàn).NSGA-Ⅱ算法種群規(guī)模為300，進(jìn)化代數(shù)為600，變異率0.1，交叉率為［0，1］上的隨機(jī)數(shù).第一代父代種群中，包括5組L-BFGS算法得到的優(yōu)化解，及295組隨機(jī)數(shù).

2.1 與單一物理目標(biāo)的對比

在一例前列腺腫瘤患者上測試了本文目標(biāo)與物理目標(biāo)對優(yōu)化結(jié)果的影響.物理目標(biāo)選擇為fmax（DOAR2），fmax（DOAR1），fDV（VOAR275）和fDV（VOAR175）.約束條件為fmin-max（DPTV），fDV（VOAR250）和fDV（VOAR150）.其中fmax（D）函數(shù)定義如下

式中：DmaxO危及器官的耐受劑量.優(yōu)化設(shè)置靶區(qū)均勻劑量為78 Gy，危及器官接受劑量為60 Gy.其他設(shè)置同混合目標(biāo)相同.兩種目標(biāo)函數(shù)均采用本文混合算法優(yōu)化.

優(yōu)化結(jié)果如圖2所示，圖2（a）是采用M-M法優(yōu)化得到的所有DVH曲線圖，圖2（b）是采用S-M法優(yōu)化得到的所有DVH曲線圖，不同線型分別對應(yīng)靶區(qū)、直腸和膀胱.從圖中可以看出，采用混合目標(biāo)函數(shù)得到的靶區(qū)劑量分布更加均勻，即有90%的靶區(qū)達(dá)到最小劑量74 Gy，且最大劑量不超過處方劑量的110%（85.8 Gy）.但對于單一物理目標(biāo)函數(shù)，劑量上限超過了85.8 Gy，存在熱點(diǎn)，影響了治療效果.圖2（c）為分別選擇兩個(gè)優(yōu)化目標(biāo)得到結(jié)果中TCP指數(shù)最高的一組對比，可以看出在確保PTV高劑量的同時(shí)，混合目標(biāo)函數(shù)明顯降低了ORA內(nèi)高劑量的分布.同時(shí)與表1危及器官的DV約束準(zhǔn)則進(jìn)行比較，混合目標(biāo)所得到的DVH曲線在高劑量區(qū)域更加符合要求.

分別計(jì)算兩種目標(biāo)函數(shù)優(yōu)化后得到數(shù)據(jù)的TCP、NTCP及HI，選擇其中TCP較高的前5項(xiàng).數(shù)據(jù)見表2.從表中可以看出由于兩種目標(biāo)函數(shù)均使用fmin-max（DPTV）作為約束條件，因此所得到的TCP和HI指數(shù)幾乎相同，但用單一的物理目標(biāo)函數(shù)與混合目標(biāo)函數(shù)相比，其對危及器官高劑量分布的控制較差，直腸和膀胱的NTCP均有增大，尤其是直腸更加明顯.可見本文提出的混合優(yōu)化函數(shù)，加入了g EUD生物目標(biāo)函數(shù)，在TCP相同的情況下，NTCP有了明顯降低，更好地保護(hù)了危及器官，治療效果更好.

本文在10例前列腺腫瘤患者中測試了混合目標(biāo)與單一物理目標(biāo)對優(yōu)化結(jié)果的影響.分別計(jì)算兩種目標(biāo)函數(shù)優(yōu)化后得到數(shù)據(jù)的TCP、NTCP及HI，每個(gè)病例均選擇TCP最大的一組數(shù)據(jù)，結(jié)果見表3.表中結(jié)果表明混合目標(biāo)對于不同的病人在TCP沒有減少的情況下，NTCP均有略微下降.

圖2 不同的目標(biāo)函數(shù)在一例前列腺患者上所得結(jié)果 Fig.2 Results of patients with prostate using different objective function

表1 危及器官的DV約束條件 Tab.1 DV constraint conditions for endanger organs

表2 混合目標(biāo)函數(shù)和單一物理目標(biāo)函數(shù)在一例前列腺患者上所得TCP、NTCP和HI值 Tab.2 TCP，NTCP and HI values obtained from the mixed objective function and a single physical target function in the same patient with prostate

表3 混合目標(biāo)函數(shù)和單一物理目標(biāo)函數(shù)在10例前列腺患者上所得TCP、NTCP和HI值 Tab.3 TCP，NTCP and HI values obtained from the mixed objective function and a single physical target function in 10 patients with prostate

2.2 與單一NSGA-Ⅱ算法和無約束混合算法的對比

在同一例前列腺腫瘤患者上，使用相同的混合目標(biāo)分別測試了M-M、M-S和M-NM.M-M參數(shù)設(shè)置與上文相同.M-S參數(shù)設(shè)置為種群300，遺傳代數(shù)600，變異率0.1，交叉率為［0，1］上的隨機(jī)數(shù).M-NM將混合目標(biāo)中的約束條件轉(zhuǎn)換為優(yōu)化目標(biāo)函數(shù)，參數(shù)設(shè)置為種群300，遺傳代數(shù)600，變異率0.1，交叉率為［0，1］上的隨機(jī)數(shù)，L-BFGS算法參數(shù)設(shè)置同M-M算法設(shè)置.優(yōu)化結(jié)果見圖3.圖3（a）是M-M優(yōu)化得到的所有DVH曲線圖，圖3（b）是M-S優(yōu)化得到的所有DVH曲線圖，圖3（c）是M-NM優(yōu)化得到的所有DVH曲線圖.圖3（d）為分別選擇三組優(yōu)化結(jié)果中TCP指數(shù)最高的一組對比.

圖3 不同的優(yōu)化法在一例前列腺患者上所得結(jié)果 Fig.3 Results of different optimization methods used a patient with prostate

從圖3中可以看出，相比單一NSGA-Ⅱ算法，在相同的種群和迭代次數(shù)設(shè)置下，混合算法能快速得到理想的DVH曲線，保證靶區(qū)的劑量分布，而單一NSGA-Ⅱ算法，靶區(qū)高劑量分布達(dá)不到90%的覆蓋率.而與無約束混合算法相比，相同的靶區(qū)分布下，危及器官內(nèi)高劑量的分布更少.因此，本文提出的混合算法優(yōu)化效率更高.在新算法中，雖然需要事先使用梯度算法優(yōu)化目標(biāo)函數(shù)，但優(yōu)化目標(biāo)已簡化，因此計(jì)算時(shí)間并無明顯增加.

3 結(jié) 語

本文提出了一個(gè)新的逆向計(jì)劃優(yōu)化算法.先構(gòu)造了一個(gè)基于物理-生物準(zhǔn)則并帶有約束條件的混合優(yōu)化目標(biāo).再使用混合算法進(jìn)行優(yōu)化.為了驗(yàn)證新算法的有效性，本文算法在10例前列腺腫瘤患者上與基于物理準(zhǔn)則的優(yōu)化目標(biāo)、單一優(yōu)化算法和無約束優(yōu)化算法進(jìn)行了對比實(shí)驗(yàn).從大量的實(shí)驗(yàn)數(shù)據(jù)中，可以看出，混合目標(biāo)混合優(yōu)化算法，可以在確保靶區(qū)劑量均勻分布的前提下，降低危及器官內(nèi)的劑量分布，提高優(yōu)化質(zhì)量.同時(shí)，新算法避免了權(quán)值的選擇，并一次得到多個(gè)DVH曲線，醫(yī)生可以具體要求進(jìn)行個(gè)性化選擇，提高了臨床治療的靈活度，具有較好的可行性.但由于優(yōu)化目標(biāo)函數(shù)和變量較多，在加速的情況下優(yōu)化時(shí)間仍較長.下一步將繼續(xù)嘗試為優(yōu)化算法加速.

［1］Wu Q，Mohan R，Niemierko A，et al.Optimization of intensity-modulated radiotherapy plans based on the equivalent uniform dose［J］.Int J Radiat Oncol Biol Phys，2002，52（1）：224-235.

［2］Stavrev P，Hristov D，Warkentin B，et al.Inverse treatment planning by physically constrained minimization of a biological objective function［J］.Med Phys，2003，30（11）：2948-2958.

［3］Li A X，Alber M，Deasy J O，et al.The use and QA of biologically related models for treatment planning：short report of the TG-166 of the therapy physics committee of the AAPM［J］.Med Phys，2012，39（3）：1386-1409.

［4］Dirscherl T，Alvarez-Moret J，Bogner L.Advantage of biological over physical optimization in prostate cancer［J］.Z Med Phys，2011，21（3）：228-235.

［5］Diot Q，Kavanagh B，Timmerman R，et al.Biological-based optimization and volumetric modulated arc therapy delivery for stereotactic body radiation therapy［J］.Med Phys，2012，39（1）：237-245.

［6］Langer M，et al.Comparison of mixed integer programming and fast simulated annealing for optimizing beam weights in radiation therapy［J］.Med.Phys.，1996，23：957-964.

［7］林琳.IMRT逆向計(jì)劃中多目標(biāo)優(yōu)化算法進(jìn)化策略的研究［D］.浙江：浙江工業(yè)大學(xué)，2009.

［8］盛大寧.IMRT逆向計(jì)劃中的混合多目標(biāo)梯度算法研究［D］.安徽：合肥工業(yè)大學(xué)，2010.

［9］閔志方，宋恩民，金人超，等.用于逆向計(jì)劃設(shè)計(jì)的整合優(yōu)化方法［J］.華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2010，38（1）：5-8.Min Zhifang，Song Enmin，Jin Renchao，et al.Integrated optimization method to design inverse planning［J］.Journal of Huazhong University of Science and Technology（Nature Science Edition），2010，38（1）：5-8.（in Chinese）

［10］Daly-Schveitzer N，Juliéron M，Gan Tao Y，et al.Intensity-modulated radiation therapy（IMRT）：Toward a new standard for radiation therapy of head and neck cancer［J］.European Annals of Otorhinolaryngology，Head and Neck Diseases，2011，128（5）：241-247.

［11］Kondoh T，Kashima H，Yang J F，et al.Dynamic optical modulation of an electron beam on a photocathode RF gun：Toward intensity-modulated radiation therapy（IMRT）［J］.Radiation Physics and Chemistry，2008，77（10）：1142-1147.

［12］Niemierko A.A generalized concept of equivalent uniform dose（EUD）（abstract）［J］.Med Phys，1999，26（6）：1110.

［13］Withers H R，McBride W H.Biologic basis of radiation therapy［M］.Philadelphia：Lippincott-Raven，1998.

［14］Murshed H.Dose and volume reduction for normal lung using intensity-modulated radiotherapy for advanced-stage non-small-cell lung cancer［J］.Int J Radiat Oncol Biol Phys，2004，58（4）：1258-1267.

［15］Burman C，Kutcher G J，Emami B，et al.Fitting of normal tissue tolerance data to an analytic function［J］.Int J Radiat Oncol Biol Phys，1991，21：123-135.

［16］Emami B，Lyman J，Brown A，et al.Tolerance of normal tissue to therapeutic irradiation［J］.Int J Radiat Oncol Biol Phys，1991，21：109-122.