基于先驗知識的頭頸部腫瘤調強放療計劃多目標優化方法

劉勁光，劉國才

1.湖南大學電氣與信息工程學院，湖南長沙410082；2.湖南工程學院計算科學與電子學院，湖南湘潭411104

前言

惡性腫瘤嚴重威脅人們的生命健康，手術、放療和化療是3 種臨床上主要的治療方法。放射治療因其適應證廣泛、療效較好，在腫瘤的治療中有著無可替代的重要地位。據報道大約有70%的惡性腫瘤病人需要接受放療，在美國、日本等國家，接受放療者約占當年新發病例的50%～60%，目前仍有上升趨勢。放射治療的主要目標是使腫瘤靶區（包括GTV、CTV、PTV）接收到足夠的高能輻照劑量以徹底消滅腫瘤細胞，同時又要使腫瘤周邊的正常組織受到盡可能小的照射劑量，以最大限度地避免正常組織受到損傷，減少放療副作用。目前臨床上廣泛采用的逆向調強放療計劃（IMRT）是使照射的劑量在整個腫瘤靶區內按照處方劑量形成梯度分布，同時又要使周邊的危及器官（OAR）滿足相應的臨床劑量約束條件。因此，腫瘤IMRT 計劃本質上是一個多目標優化問題，但因該多目標優化問題的目標函數很多（頭頸部腫瘤通常有20多個），且待優化的參數個數非常多（通常大于1 000 個），目前沒有精確求其Pareto 最優解集（Pareto optimal solution set）的有效方法。臨床上通常將多目標的IMRT 計劃優化問題轉化為單目標優化問題求解，尋找臨床上可接受的次優解。目前臨床上IMRT計劃通常采用兩步優化方法［1］：

步驟1：臨床放射物理師根據各靶區及其危及器官的臨床劑量約束條件構造一組適當的目標函數fi，i=1,2,…,n，并根據經驗對每個目標函數賦予適當的權重wi，i=1,2,…,n，將這組目標函數加權求和形成一個單目標優化問題。通過優化該單目標可以獲取各照射野的強度分布或通量圖，這一步也稱為通量圖優化（Fluence Map Optimization,FMO）。

步驟2：將各方向照射野的通量圖分解為一組臨床上醫用直線加速器可執行的照射子野或口徑的加權疊加，各子野的形狀由多葉準直器（MLC）形成。各子野的權重通常稱為子野機器跳數（MU），由直線加速器的出束率與出束時間實現。

兩步優化的主要優點是第一步的優化過程，即通量圖優化，其目標函數是凸函數或可轉化為凸函數，可利用很多方法求解。通量圖優化模型如下：

其中，X是待優化參數向量，表示各射野單元野的權重，即照射強度，在區域Ω 內取值。D是劑量貢獻矩陣，其元素Dij表示單位強度照射下的第i個單元野對第j個體素劑量沉積的貢獻大小。d是空間感興趣區域的劑量分布向量。由于權重因子wi，i=1,2,…,n，無直接臨床意義，僅用于調整各目標函數的相對重要程度，故要選擇一組合適的權重因子非常困難，需要設計者在計劃優化過程中根據劑量分布的情況和劑量體積直方圖等對優化結果進行可視化評估，多次修正權重因子并反復優化，直至找到臨床可接受的放療計劃為止。

為了避免多次修正模型（1）中的權重因子，可直接求解多目標優化模型（2），即多個目標函數同時最優化并選取最優解。

多目標優化的目的是使各目標函數同時最優，達到最小值。然而由于各目標函數間的相互沖突，一般情況下不可能所有目標函數同時達到最小。因此，工程上通常求其Pareto 最優解集［2］。任何一個Pareto 最優解均滿足條件：如果其中某個目標函數值減小，則必定導致至少一個其它的目標函數值增大。全體Pareto 最優解所對應的目標函數值集合構成了目標函數值空間的一個Pareto前沿。

目前逼近Pareto前沿的方法較多，主要有夾心算法和進化算法兩大類［3-4］。若所有目標函數都是凸函數，則可以采用夾心算法［5-6］逼近多目標優化問題（2）的Pareto 前沿。由凸分析理論［7］可知，在Pareto 前沿上可逐步選取一組適當的點，這組點稱為錨點，連接這組錨點所構成的超凸多邊形成為Pareto 前沿的一個“上界面”，而由通過這組錨點的切平面所構成的超凸多邊形成為Pareto 邊界的“下界面”，Pareto 前沿被“夾在”以上兩個“上、下界面”之間。因此，這兩個“上、下界面”之間的最大距離可用來度量它們逼近Pareto前沿的程度。

當找到Pareto 前沿后，IMRT 計劃設計者只要通過可視化的導航便可以選出臨床上最優的Pareto解［8-9］，即最優的臨床放療計劃。因此，逼近Pareto 前沿是多目標優化的核心。但是直接逼近模型（2）中的Pareto 前沿一般比較困難，其原因主要有幾點：首先，目標函數的個數較多，特別對頭頸部腫瘤，模型（2）中的目標函數往往超過20 個。由于逼近Pareto前沿所需要的錨點個數隨目標函數的增加呈指數增長并且收斂速度也會隨目標函數的增加而變慢，因此，對頭頸部腫瘤的多目標優化問題，直接逼近Pareto 前沿幾乎無法實現；其次，待優化參數X 的個數也非常大，頭頸部腫瘤一般超過1 000。為減少計算量，提高優化速度，一方面可以重新定義多目標調強放療計劃優化的目標函數，通過減少目標函數的個數來降低Pareto前沿的維數；另一方面也可在優化過程中添加適當的約束條件［10］，限定Pareto前沿的范圍，即僅逼近有臨床意義的部分Pareto前沿。但如何選取約束條件目前尚無具體標準，這是因為選取的約束條件如過于苛刻，有可能導致無滿足約束條件的放療計劃，因而無法生成Pareto前沿；反之，如選取的約束條件過于寬松，可能導致生成的Pareto前沿過大，在逼近過程中需要確定更多的放療計劃［11］，即需要更多Pareto前沿上的錨點。

已有文獻表明腫瘤IMRT計劃對同類腫瘤的新患者放療計劃設計具有一定的指導意義［12-13］，如Chanyavanich等［12］將新的前列腺腫瘤病人與已有的100個前列腺病人樣本影像進行匹配，找出影像互信息最大的樣本，將該樣本的放療計劃作為新病人的放療計劃。實驗結果表明該方法能使靶區PTV受到的照射劑量與傳統計劃相當，同時對危及器官的保護程度也類似，但無需重新進行放療計劃的優化，因而基于先驗信息的放療計劃優化比傳統IMRT方法計算速度更快。此外，Tol等［13］也得到類似的結果。

本文首先根據同類腫瘤放療劑量的先驗知識，對各危及器官的線性EUD 按相關系數進行R 型聚類分析［14］，將危及器官分為3類；其次，將所有屬于同類危及器官的線性EUD 均值看作目標函數，聯合腫瘤靶區的約束條件，構造了一個維度較低的多目標優化模型；最后，通過基于角解的增強夾心算法逼近該低維、局部的Pareto 前沿，并通過可視化導航方法確定一個實時的調強放療計劃。對比傳統的多目標調強放療計劃優化，本文方法的創新點主要有兩個：其一，通過對頭頸部腫瘤的各危及器官進行R型聚類分析，并按類構造目標函數，顯著減少了多目標IMRT計劃優化的目標函數個數，降低了Pareto前沿的復雜度；其二，在逼近Pareto 前沿時提出了基于角解的增強夾心算法，能比原增強夾心算法生成更完整的近似Pareto前沿。實驗結果表明，本文方法可以將頭頸部腫瘤IMRT計劃多目標優化模型中的20～21個目標函數減少到6～7 個，并能有效確定其近似的Pareto 前沿，找到與臨床實際放療計劃質量相當、部分指標甚至更優的放療計劃。

1 基于先驗知識的多目標IMRT計劃優化

1.1 頭頸部腫瘤IMRT計劃優化靶區目標函數的構造

對腫瘤靶區，主要考慮均勻劑量、最大劑量和最小劑量3類約束條件。

參考Rivera 等［15-16］的目標函數構造方式，靶區目標函數定義為式（3）：

其中Target 為計劃優化時所考慮的靶區，如GTV、CTV、PTV。NTarget表示該靶區的體素總數，di為靶區第i個體素的劑量值。

1.2 頭頸部腫瘤IMRT計劃優化危及器官目標函數的構造

對頭頸部腫瘤放療計劃優化，我們在計劃設計過程中主要考慮Brain Stem、Left eye、Left lens、Right lens、Right eye、L_temporal lobe、R_temporal lobe、Left TM_joint、Right TM_joint、Spinal cord、Larynx、Mandible、Left optic nerve、Right optic nerve、L_parotid、R_parotid、Optic chiasm、Pituitary 共18 個危及器官，IMRT 計劃目的是使腫瘤受到更高輻射劑量照射的同時，保證所有危及器官的并發癥（Normal Tissue Complication Probability,NTCP）最小。而NTCP 可以看成等效均勻劑量（Equivalent Uniform Dose,EUD）的函數［17］，EUD 又可進一步由危及器官劑量的平均值和最大值線性組合來擬合［8,18］，即線性EUD。因此NTCP 最小本質上要求危及器官的線性EUD最小。我們采用線性EUD作為危及器官的目標函數，即

其中，αOAR為系數因子，其具體取值與危及器官的性質有關（表1），如對串行器官腦干、脊椎等，αOAR的取值較小，而對并行器官如雙顳、腮腺等，αOAR的取值較大。

表1 系數因子αOAR的取值Tab.1 The value of coefficient factor αOAR

如直接將頭頸部腫瘤中每個危及器官對應的目標函數fOAR看成目標函數，聯合腫瘤靶區約束條件進行放療計劃優化，則優化目標有20多個，要生成多目標優化模型（2）中的Pareto 前沿幾乎無法實現。因此，我們基于放療計劃的先驗知識，回顧性地計算了71 位同類患者IMRT 計劃中所有危及器官的線性EUD 函數值，再對所有危及器官按線性EUD 進行R型聚類分析［19］，其基本步驟為：

首先，各危及器官單獨為一類，對線性EUD劑量，計算危及器官兩兩之間的Pearson相關系數；其次，將相關系數最大的兩類器官合并為一新類，同時剔除該器官所在的原始類；再次，計算新類與其他類之間的相關系數，重新生成新類，直至所有危及器官聚為一類為止。其中，類與類之間的相關系數利用類平均法進行計算，即對兩類E1和E2，則兩類之間的相關系數為：

其中，r(i,j)為危及器官i、j之間的相關系數，i、j分別屬于類E1、E2。如考慮危及器官A1～A6 的R 型聚類分析，假設編號依次為1～6，危及器官之間的相關系數矩陣見表2。

表2 相關系數矩陣Tab.2 Correlation coefficient matrix

進行R型聚類分析時首先每個器官各成一類，即聚為6類的結果為{{1},{2},{3},{4},{5},{6}}，根據器官之間的相關系數矩陣及類平均法的計算公式（5）可知，類{1},{2}之間的相關系數最大，為0.763，因而{1},{2}兩類可合并為一個類{1,2}。所以R 型聚類分為5 類的結果為{{1,2},{3},{4},{5},{6}}，計算其余各類與類{1,2}之間的相關系數。

對比5類之間的相關系數，類{3},{4}之間的相關系數最大，為0.708。因此，進一步聚為4 類的結果為{{1,2},{3,4},{5},{6}}，計算其余各類與類{3,4}之間的相關系數。對比4類之間的相關系數，類{1,2},{3,4}之間的相關系數最大，為0.407，因此，進一步聚為3 類的結果為{{1,2,3,4},{5},{6}}。依此類推，直至所有危及器官全部被合成一類為止，其聚類的樹狀圖見圖1。

圖1 R型聚類分析的樹狀圖Fig.1 Tree diagram of R-mode cluster analysis

對回顧性的71 例頭頸部腫瘤患者，放療計劃設計過程中共考慮了18 個主要危及器官。按R 型聚類分析方法，其聚類的樹狀圖見圖2。

圖2 危及器官的R型聚類分析樹狀圖Fig.2 Dendrogram of R-mode clustering for organs-at-risk

從圖2可知，如將所有危及器官分為3 類，則具體的分類為：

一類：Left eye,Left lens,Right lens,Right eye；二類：Spinal cord,Brain stem；三類：Pituitary,Left optic nerve,Right optic nerve,L-parotid,R-parotid,Mandible,Larynx,Left TM-joint,Right TM-joint Ltemporal lobe,R-temporal lobe,Optic chiasm。

對于每一類，我們將所有危及器官的線性EUD平均值看作多目標IMRT計劃優化的目標函數，即

其中，ni表示各類中危及器官的個數。因此，在多目標IMRT 計劃優化過程中，如對所有危及器官按類構造目標函數，則可將原多目標IMRT計劃優化中的18個危及器官目標函數縮減為3個，極大程度減少了待優化的目標函數個數。

1.3 多目標優化的Pareto前沿近似方法

對放療計劃的多目標優化模型（2），如直接生成其對應的Pareto 前沿，需要更多的Pareto 最優解集，且大多解對應的放療計劃無臨床意義，因此我們考慮對腫瘤靶區添加適當的先驗約束，使生成的Pareto前沿僅包含臨床有意義的部分。參考Craft等［20］前期的工作，對腫瘤靶區PTV，定義斜坡函gPTV如下：

其中，NPTV表示靶區PTV 的體素總數，di為PTV 的第i個體素劑量值，并定義先驗約束條件如下：

將上述先驗約束條件（8）加入到模型（2）中，得多目標IMRT計劃優化模型（9）：

1.3.1 目標函數值范圍界定為生成模型（9）中的Pareto 前沿，現有的夾心算法需要對各目標函數進行歸一化，即將目標函數的取值范圍歸一化到［0,1］之內，因而需要事先確定模型中各目標函數的取值范圍。目前，確定目標函數取值范圍的方法是令各目標函數值分別最小，通過優化尋找到一組初始錨點，再根據各錨點的取值確定目標函數的取值范圍［6］。該方法所確定的范圍可能只是其實際取值范圍的一部分，因而生成的Pareto 前沿可能不完整。如，考慮下面3個目標函數的優化問題。

其中，函數fi=xi，i=1,2,3，該多目標優化問題對應的Pareto前沿見圖3。

圖3 多目標優化問題（10）的Pareto前沿Fig.3 Pareto frontier of multi-criteria optimization(10)

如通過令各目標函數f1、f2、f3分別最小來確定各目標函數的取值范圍，可以得到3 個錨點，分別為（0,0.5,0.5）、（0.5,0,0.5）、（0.5,0.5,0）。由這3 個錨點所確定的目標函數取值范圍僅為［0,0.5］，并非目標函數的實際取值范圍［0,1］。為盡可能生成完整的Pareto 前沿，本文利用Singh 等［21］在多目標函數值空間降維過程中提出的角解方法，確定各目標函數的取值范圍，即在多目標優化模型中，依次令每個目標函數最小，同時，也讓剩余n-1 個目標函數之和最小，得2n個Pareto 最優解，這些解稱為角解。目標函數的取值范圍由這2n個角解確定，如對上述3 個目標函數的多目標優化問題，分別令f1、f2、f3、f2+f3、f1+f3、f1+f2最小，求得6個角解，見表3。

根據表3可知，目標函數f1、f2、f3的取值范圍都是［0,1］，與實際取值范圍一致。

表3 多目標優化的6個角解Tab.3 Six corner solutions of multi-criteria optimization

1.3.2 Pareto 前沿近似方法為求解模型（9），我們聯合Rennen 等［6］提出的增強夾心算法與Singh 等［21］生成角解的方法，提出了逼近Pareto 前沿的新方法，即基于角解的增強夾心算法，該算法基本思路如下：

step1 生成初始錨點；

在模型（9）中，依次令每個目標函數最小，同時，也讓剩余n-1 個目標函數之和最小，得模型（9）Pareto前沿上的2n個角解，即初始錨點：

zA1,…,zAn和zM1,…,zMn

step2 添加輔助點，生成凸殼，形成“上界面”；

設所有錨點組成的集合為Z，令ziub為Z 中第i個目標函數的最大值。為使錨點生成凸殼面時各面的法向量非負，對集合Z中的每個錨點z=(z1,...,zn)T，添加n個輔助點dj(z)，j=1,2,…,n，如下：

其中，θ＞0是任意給定非常小的正數。

所有輔助點的集合記為E，錨點Z與輔助點E共同生成的凸殼為IPS，則該凸殼面形成了Pareto 前沿的一個“上界面”（圖4）。

step3 生成“下界面”；

對集合Z 中每個錨點z，Pareto 前沿上存在過該點的切平面λTp=λTz，其中λ是該切平面的法向量，令

則OPS 的邊界構成了Pareto 前沿的一個“下界面”（圖4）。

圖4 基于角解的增強夾心算法示意圖Fig.4 Schematic diagram of enhanced sandwich algorithm based on corner solutions

step4 確定“上、下界面”誤差；

對凸殼面IPS，我們按照Solanki 等［22］的 算 法xnise2提取相應的超平面及法向量。

步1：在IPS 中找出有n個點屬于Z的超平面；所有超平面的集合記作IPSnld，并在IPS中剔除相應的超平面；

步2：在IPS 中找出有n-1 個點屬于Z的超平面，如這n-1 個點屬于IPSnld中的某個超平面，則剔除這n-1個點所在的超平面；

步3：更新n，讓n=n-1，如n≥1，返回步2，否則進入步4；

步4：修正IPSnld。在IPSnld中，如存在多個超平面，其屬于Z中的錨點相同，則將這多個超平面合并成一個過這些相同錨點的超平面，其法向量取原多個超平面法向量的平均值；

步5：確定“上、下界面”誤差。對IPSnld中的每個超平面，設其對應的方程為wTp=b，求解模型（11），得最優解:

計算該超平面到“下界面”之間的距離e：

其中，ε為參數向量，當無特別要求時常取ε的每個分量為對應目標函數取值區間的長度，特別當目標函數的取值范圍為［0,1］時，此時ε的取值為ε=(1,1,…,1)T，如同時w的各分量之和等于1，則誤差e=b-wT變為Solanki 等報道的夾心算法中的逼近誤差。

“上、下界面”間的誤差e*定義為IPSnld中所有超平面到“下界面”之間距離的最大值。

step5算法終止條件；

如“上、下界面”誤差e*小于事先給定的閾值，則算法終止，輸出所有錨點的集合Z。否則，從IPSnld中選取一個具有最大誤差e*的超平面F*，設該面F*的方程為wT0p=b0，進入step6。

step6 生成新的錨點；

在模型（9）中，對各目標函數賦權，權重為超平面F*的法向量w0，將原多目標優化問題可轉化為單一目標函數的最優化模型（13）。

求解模型（13），得F(d)的最優解為如z*在超平面F*上，令超平面F*到“下界面”的距離為0，更新“上、下界面”之間的誤差e*，返回step5，否則進入step7。

step7 更新錨點集合

更新錨點的集合Z=Z∪{z*}，返回step2。

根據step1～step7，基于角解的增強夾心算法的具體流程見圖5。

圖5 Pareto前沿近似流程圖Fig.5 Flow chart of Pareto frontier approximation

1.4 Pareto前沿導航

由上述基于角解的增強夾心算法，可以在Pareto前沿上尋找到一組錨點，這些點的凸殼構成了Pareto前沿的近似。設計者可以通過實時導航方法［9,23-24］，并根據臨床放療劑量評估指標（劑量分布情況、劑量體積直方圖等），可視化地從眾多Pareto 最優IMRT計劃中快速選出一個導航放療計劃作為臨床最終的執行計劃。本文利用Craft［25］的Pareto 前沿導航軟件，聯合放療計劃軟件CERR［26］進行實時導航選取合適的放療計劃。

2 實驗結果

2.1 鼻咽癌IMRT計劃

本文回顧性選取71 例鼻咽癌患者，并從中選取兩例患者進行對比實驗驗證分析。患者臨床IMRT計劃的PET/CT 影像和IMRT 計劃數據均來源于湖南省腫瘤醫院放療中心。兩例患者的放療靶區與危及器官（腦干等）由臨床放療專家通過商用放療計劃系統Nucletron Oncentra MasterPlan v3.2 手動勾畫完成。臨床IMRT計劃的射束角度為0o、51o、103o、155o、206o、257o、309o。

患者1 的放療靶區包括GTVPET、CTV、PTV，所有靶區都考慮最大劑量、最小劑量約束，同時對靶區GTVPET考慮均勻劑量約束，具體約束條件見表4。放療計劃設計優化過程考慮的危及器官有17 個，按R型聚類分析結果，共分為3類：

表4 患者1的靶區劑量約束Tab.4 Target dose constraints for patient 1

1類（OAR1）：Left eye,Left lens,Right lens,Right eye；

2類（OAR2）：Spinal cord,Brain stem；

3類（OAR3）：Optic chiasm,Left optic nerve,Right optic nerve,Mandible,Larynx,Left TM-joint,Right TM-joint,L-temporal lobe,R-temporal lobe,Lparotid,R-parotid

對這3類危及器官中的每一類，其對應的目標函數按式（6）構造，共3 個目標函數。對該患者，IMRT計劃設計優化過程中考慮的靶區共有3個，對應的目標函數按式（3）構造，共3 個目標函數。因此，患者1的放療計劃優化問題，其實是6個目標的多目標優化問題。

患者2 的放療靶區包括PET-GTVnx、PETGTVnx+1.0、UP of CTV、Down of CTV，所有靶區都考慮最大劑量、最小劑量約束，同時對靶區PETGTVnx考慮均勻劑量約束，具體約束條件見表5。放療計劃設計優化過程考慮的危及器官有17個。對比患者1，在第3 類危及器官中增加了Pituitary，但減少了危及器官larynx。患者2 的放療計劃優化問題，其實是7個目標的多目標優化問題。

表5 患者2的靶區劑量約束Tab.5 Target dose constraints for patient 2

2.2 Pareto前沿逼近

對上述兩例患者，采用基于角解的增強夾心算法逼近模型（9）中的Pareto 前沿。在逼近過程中，預設允許誤差為5%。對患者1，優化過程中共生成了57 個錨點；對患者2，優化過程中共生成了116 個錨點。這些錨點的凸殼面構成了Pareto 前沿的一個逼近。由于這些都是6 或7 維的錨點，無法在笛卡爾直角坐標系中直觀地顯示，我們繪制出了其平行坐標圖［27］，即橫坐標代表錨點的各個分量，縱坐標代表各分量的取值，平行坐標圖中的每根折線代表一個錨點，見圖6和圖7。可以看出對每個目標函數的取值，都存在一定的取值區間，因而可以進行放療計劃的導航并進行可視化的放療計劃評估。

圖6 患者1所有錨點的平行坐標Fig.6 Parallel coordinates of all anchor points for patient 1

圖7 患者2所有錨點的平行坐標Fig.7 Parallel coordinates of all anchor points for patient 2

2.3 放療計劃導航與評估

對已有的Pareto 前沿近似曲面，使用Craft 等的Pareto 前沿導航軟件，聯合放療計劃軟件CERR 進行驗證評估，可在Pareto前沿的近似曲面上快速選擇一個導航計劃。圖8顯示患者1 的6 個目標導航界面。中間的可滑動選擇控件對應上述6 個目標函數的取值，中間的可滑動小方塊可以上下移動，向上移動表明該目標函數值變大，反之表明縮小目標函數的取值。圖中所有小方塊所在位置對應一個導航放療計劃，其劑量分布如圖9所示。

圖8 導航界面（從左至右,依次為GTVPET、CTV、PTV、OAR1、OAR2、OAR3對應的目標函數）Fig.8 Navigation interface(from left to right,the criteria corresponding to GTVPET,CTV,PTV,OAR1,OAR2,OAR3,respectively)

為驗證該方法的有效性，本文將選取的導航計劃與臨床計劃分別從劑量分布、劑量體積直方圖、靶區劑量分布的定量指標、危及器官劑量分布的定量指標4個方面進行對比分析，具體方法如下。

2.3.1 劑量分布對比對比患者1 臨床計劃的劑量分布圖10，導航計劃（圖9）與臨床計劃有著近似的劑量分布，劑量都呈現由靶區向外圍逐漸遞減，導航計劃的最大劑量為82.5 Gy，比臨床計劃的76.3 Gy 稍高。靶區GTVPET的劑量分布均勻在處方劑量68 Gy附近，圖中"+"顯示了GTVPET（中心紫色線）內部同一個位置導航計劃與臨床計劃的劑量值，分別為68.75 Gy和71.67 Gy，在該處導航計劃的劑量更接近GTVPET的處方劑量68 Gy。

圖9 導航計劃對應的劑量分布Fig.9 Dose distribution of navigation plan

圖10 臨床計劃對應的劑量分布Fig.10 Dose distribution of the clinical plan

2.3.2 劑量體積直方圖對比為進一步評估靶區與各危及器官的劑量分布情況，我們繪制出兩例患者導航計劃與臨床計劃的靶區和危及器官DVH 圖。圖11是患者1各個靶區劑量的DVH 圖。從該圖可以看出導航計劃與臨床計劃有著相似的劑量分布，且對靶區GTVPET，導航計劃有更均勻的劑量分布。圖12～14是患者1 的3 類危及器官劑量分布DVH 圖。對第1類危及器官，除左晶體（Left lens）外，其它危及器官的劑量分布導航計劃都要優于臨床計劃（圖12）；對第2 類危及器官，脊髓（Spinal cord）的劑量分布兩種計劃劑量分布相當，但對腦干（Brain stem）的劑量分布，臨床計劃的結果更好（圖13）；對第3類危及器官，絕大多數危及器官的劑量分布在導航計劃下要比臨床計劃更優（圖14）。

圖11 靶區的DVH圖Fig.11 Dose-volume histogram of target areas

圖12 第1類危及器官劑量分布的DVH圖Fig.12 Dose-volume histogram of organs-at-risk of the first category

圖13 第2類危及器官劑量分布的DVH圖Fig.13 Dose-volume histogram of organs-at-risk of the second category

圖14 第3類危及器官劑量分布的DVH圖Fig.14 Dose-volume histogram of organs-at-risk of the third category

患者2 進一步驗證了導航計劃能生成與臨床相當且對部分危及器官保護更好的放療計劃。圖15是患者2 的一個導航計劃中各個靶區劑量分布的DVH圖。對靶區UP of CTV 和Down of CTV，導航計劃將劑量的平均值由臨床計劃的69.10和65.27 Gy提升到了71.02 和66.36 Gy；對靶區PET-GTVnx，雖然導航計劃的劑量平均值只有74.04 Gy，相對臨床計劃的劑量平均值（75.67 Gy）較低，但導航計劃的劑量分布更均勻，且更接近處方劑量74 Gy。圖16是同一導航計劃主要危及器官的DVH 圖。對大多數的危及器官，導航計劃的劑量分布要優于臨床計劃，但少數危及器官例外，如對危及器官L-parotid、R-parotid，導航計劃的最大劑量要比臨床計劃大。

圖15 患者2靶區的DVH圖Fig.15 Dose-volume histogram of target areas in patient 2

圖16 導航計劃部分危及器官劑量分布的DVH圖Fig.16 Dose-volume histogram of some organs-at-risk in navigation plan

2.3.3 靶區均勻性指數與適形度指數對比為進一步定量評估靶區劑量的分布情況，我們按下面公式分別計算出兩例患者的各個靶區在導航計劃和臨床計劃下的均勻性指數（HI）［28］和適形度指數（CI）［29］。

其中，D5和D95表示靶區5%與95%體積所得到的絕對照射劑量，Dp為靶區的處方劑量，TVRI表示被等劑量面包繞的靶區體積，TV 為靶區的體積。HI 反應靶區的劑量分布均勻程度，值越小越好，而CI指數的值越大，說明被等劑量面包繞的靶區體積越大。

實驗中兩例患者HI和CI的計算結果見表6和表7。對患者1 的所有靶區，導航計劃的HI 指數值更小，因而靶區劑量分布具有更好的均勻性，但對靶區GTVPET，導航計劃下的CI 指數稍低（表6）。對患者2的靶區UP of CTV、Down of CTV，導航計劃具有更好的CI；對靶區PET-GTVnx，導航計劃的HI 指數值更小，因而對該靶區導航計劃下的劑量分布具有更好的均勻性，但對其它評價指標，導航計劃要比臨床計劃差（表7）。也許通過適當更換導航計劃，可以進一步消除這種現象。

表6 患者1靶區定量評估指標HI與CITab.6 Quantitative evaluation indexes HI and CI of target areas in patient 1

表7 患者2靶區定量評估指標HI與CITab.7 Quantitative evaluation indexes HI and CI of target areas in patient 2

2.3.4 危及器官的并發癥對比為進一步定量評估危及器官的劑量分布情況，我們根據Lyman-Kutcher-Burman（LKB）模型［30-31］計算出部分主要危及器官的并發癥（NTCP）：

其中，式中vj表示危及器官的第j個體素占總器官的體積比，K表示危及器官的體素總個數，dj為第j個體素受照劑量，n、m、TD50為參數，具體取值見表8。對實驗中考慮的所有危及器官，表9列出了其中13 個主要危及器官導航計劃與臨床計劃的NTCP 結果。從表9的NTCP 結果可以看出，導航計劃對大多數危及器官的保護要好于臨床計劃，如實驗中的右腮腺（R-parotid），對比臨床計劃，兩例患者的NTCP 都有不同程度的減少，患者1 的NTCP 值由原來的0.12 減少到0.11，減少了8.3%；患者2 的NTCP 值由原來的0.084 減少到0.032，減少了61.9%。但也有個別危及器官的NTCP增大，如患者2中的左腮腺（L-parotid），NTCP的值從臨床計劃的0.07增加到0.13。

表8 NTCP中參數n,m,TD50的取值Tab.8 Values of parameters n,m and TD50 in NTCP

表9 主要危及器官的NTCPTab.9 NTCP of primary organs-at-risk

3 討論

本文回顧性對71 例頭頸部腫瘤患者按R 型聚類分析方法，將放療計劃設計過程中的18 個主要危及器官聚為3類，因而能極大程度地降低優化目標函數的個數，但具體聚為多少類合適仍需進一步探索。雖然，有眾多的評價指標幫助我們從聚類的樹狀圖中選擇最佳的聚類個數，但最佳聚類個數與評價指標的選取有極大關系，甚至對類別個數已知的聚類問題，不同的評價指標會得到完全不一樣的聚類個數［32］。為能有效逼近Pareto前沿，最終的聚類個數選取不能太大，否則進行多目標IMRT 計劃優化時對應的Pareto 前沿維度會非常高，無法生成其近似曲面。因此，我們選擇Dunny 指標作為評價指數，對圖2中聚類結果的樹狀圖分別計算了聚為2、3、4、5 類的Dunny 指標結果［33-34］，按指標值越大越好的原則，選取將所有危及器官聚為3類作為最佳的聚類個數。

為便于計算各危及器官的線性EUD 劑量，我們選取了適當的系數因子α。雖然α的確定與危及器官是串行或并行相關，但缺少大量的實驗數據及相關資料佐證，如α 取值為1，則危及器官的線性EUD 變為該器官劑量的平均值；如α的取值為0，則危及器官的線性EUD 變為最大劑量值。因此，當α 取值介于0和1 之間時，危及器官的線性EUD 介于其最大劑量和平均劑量之間。為進一步探索系數因子α 對優化結果的影響，對實驗患者中的第1 例，我們將表1中系數因子α 的取值做適當修改，即系數因子小于0.5的取為0.05，大于0.5 的取為0.95，在模型（13）中讓各目標函數的權重相等并求出不同系數因子下的兩個放療計劃。對比系數因子按表1取值的計劃（plan1），修改系數因子α 后的計劃（plan2）對優化結果的影響并不大。修改系數因子的主要差別，體現在是約束危及器官劑量的最大值還是危及器官劑量的平均值上。例如，對危及器官Right TM-joint，將系數因子由表1中的0.2 修正為0.05 后，該危及器官劑量的最大值由計劃plan1 中的61.9 Gy 降低到計劃plan2 中的60.7 Gy，但對該危及器官，系數因子取0.2 比0.05 可能更合理，因為系數因子的這兩個取值，目標函數同樣都是側重考慮最大劑量約束，但前者更大程度地考慮了該危及器官平均劑量的約束。因此，計劃plan1 中的劑量平均值明顯要低于計劃plan2，分別為35.76 和42.37 Gy。類似的情形在危及器官Mandible中再次出現，系數因子取0.7 比0.95 在適當增大危及器官Mandible的平均劑量情況下（22.65和21.27 Gy），可使該危及器官的最大劑量值由plan1 中的71.3 Gy降低到計劃plan2中的67.1 Gy。

4 總結

將頭頸部腫瘤臨床IMRT 計劃先驗知識集成到放療計劃多目標優化模型中，通過對各危及器官按線性EUD 進行R 型聚類分析，可將頭頸部腫瘤IMRT計劃的多目標優化模型中的20～21 個目標函數減少到6～7 個目標函數，因而極大程度簡化了其Pareto 前沿。同時，為進一步有效逼近其Pareto 前沿，在增強夾心算法的基礎上，我們提出了基于角解的增強夾心算法，使生成的Pareto 前沿近似曲面更完整。在Pareto 前沿的近似曲面上，借助可視化導航，放療計劃設計者可以快速地確定一個導航計劃，并能實時觀察該導航計劃的劑量分布、DVH 圖等劑量信息。對比臨床計劃，實驗結果表明該方法能生成與臨床IMRT 計劃質量相當、部分指標更優的放療計劃，且能根據需要實時調節靶區與各危及器官的劑量分布。但本文在優化過程中對腫瘤靶區僅考慮將均勻劑量、最大及最小劑量作為約束條件，并未考慮靶區的體積約束情形，因此劑量體積約束的多目標優化問題是我們進一步需要研究的內容之一。