針對傳統直接子野優化算法（DAO）收斂速度慢、易停滯、全局搜索能力低的缺點，本文提出一種基于梯度信息的直接子野優化方法（GDAO）。在 GDAO 中分別采用不同的優化方法對子野形狀和子野權重進行迭代優化。首先為提高子野形狀優化時每次搜索的有效性，對傳統模擬退火算法（SA）進行了改進，將梯度信息融合在 SA 算法中。采用基于梯度的 SA 法確定子野形狀，并在優化同時充分考慮多葉準直器（MLC）葉片間的約束條件，保證優化后的子野形狀滿足臨床放射治療的要求。之后再利用計算量少、迭代代價低、收斂快且穩定的梯度類具有求解大規模約束優化問題能力的帶約束最小存儲擬牛頓算法（L-BFGS-B）優化子野權重。實驗結果表明，與傳統 SA 算法相比，新算法計算時間減少了 15.90%，同時得到的治療方案靶區最低劑量提高了 0.29%，最高劑量降低了 0.45%；危及器官膀胱最高劑量降低了 0.25%；危及器官直腸最高劑量降低了 0.09%，說明在調強放射治療（IMRT）中采用 GDAO 方法直接優化子野，可在短時間內得到滿足臨床要求并可直接實施照射的治療方案，具有較好的臨床實用價值。
引用本文： 楊婕, 張鵬程, 張麗媛, 桂志國. 一種基于梯度信息的直接子野優化算法. 生物醫學工程學雜志, 2018, 35(3): 358-367. doi: 10.7507/1001-5515.201609041 復制
引言
調強放射治療[1]（intensity-modulated radiotherapy，IMRT）是當前最常用的臨床放射治療技術。它采用多個方向強度不均勻的的射野照射靶區（planning target volume，PTV），產生與 PTV 高度適形的劑量分布，同時減少對危及器官（organ at risk，OAR）和正常組織的損傷，明顯提高了治療精度。世界一流的治療計劃系統（treatment planning system，TPS）大部分都提供了 IMRT 功能，如 Elekta Inc.公司的 PrecisePLAN 系統、Philips Medical Systems 公司的 P3IMRT 系統和 Varian Medical Systems 公司的 Eclipse Integrated Treatment Planning 系統。傳統的 IMRT 一般采用“兩步法”：先優化注量圖（fluence map optimization，FMO）得到最優強度分布，再將最優強度分布分割成多葉準直器（multileaf collimator，MLC）序列。“兩步法”的缺點在于優化注量圖時并未考慮 MLC 約束，射野分割后得到的強度分布與最優強度分布之間的差別可達 5.7%[2]，這種差別降低了劑量分布的適形度，影響了放射治療質量。為了克服“兩步法”照射強度分布與最優強度分布之間不一致的缺點，直接子野優化（direct aperture optimization，DAO）[3]通過逆向計劃直接對子野的形狀和權重進行優化，得到滿足 MLC 約束的子野，降低了治療的不確定性，從而成為近幾年較熱門的 IMRT 優化方法。
DAO 屬于 NP-hard 的問題，其中子野形狀和權重是兩組相互耦合的參數，每個特定的子野形狀都要優化其對應的權重。為了解決這個問題，Shepard 等[3]采用模擬退火算法（simulated annealing，SA）對子野形狀和權重同時進行優化，他們證明采用 DAO 技術可以得到適形度較高的劑量分布。隨后 Earl 等[4]對這種方法進行了改進，只使用矩形的子野形狀來實現各射野方向的強度調強。李永杰等[5]以二維編碼的方式來表達每個子野形狀，提出了遺傳算法（genetic algorithm，GA）優化子野形狀以及共軛梯度算法（conjugate gradient，CG）優化子野權重的算法。Cotrutz 等[6]也采用遺傳算法直接對子野進行了優化。Preciado-Walters 等[7]、Romeijn 等[8]、Men 等[9]、Salari 等[10]、Bednarz 等[11]，以及 Salari 等[12]都利用列生成的方法直接生成子野形狀。Cassioli 等[13]提出一種利用梯度算法優化子野形狀的方法，但該方法只有在葉序變動范圍小的時候才有效，因此他們還特為算法設置了信任區域。王捷等[14]則先采用 CG 優化子野權重，然后再利用枚舉法優化子野形狀。以上方法均取得了較好的方案優化結果，但是仍存在一些問題。在子野形狀優化算法中，文獻[3-6]采用的是隨機搜索算法，優化過程中有一定概率跳出局部最小值，但不保證每次的優化結果為全局最優；文獻[13]采用局部梯度的方法，由于葉序位置與函數梯度并沒有嚴格的線性關系，所以不能實現全局最優；文獻[7-12]使用圖論的方法不能完全描述目標函數與 MLC 葉片位置的復雜關系，影響了放療方案的質量。因此如何更好地優化子野形狀，是一個亟待解決的問題。同時，在子野權重優化方面，也需有所改進。
綜上所述，本文提出一種基于梯度信息的直接子野優化方法（gradient-based direct aperture optimization，GDAO）用于 IMRT 逆向計劃。在 GDAO 中分別采用兩種不同的優化方法對子野形狀和子野權重進行迭代優化。首先，為了增加子野葉片位置的精度，采用基于梯度信息的 SA 算法優化子野形狀，根據葉片位置與最終劑量的關系提出一種含有梯度信息的狀態產生函數與含有梯度信息的狀態接受函數，并在優化同時充分考慮 MLC 葉片間的約束條件。最后再利用計算量少、迭代代價低、收斂快且穩定的梯度類具有求解大規模約束優化問題能力的帶約束最小存儲擬牛頓算法[15]（limited-memory BFGS for bound-constrained，L-BFGS-B）優化子野權重。
1 本文方法
1.1 目標函數
在 IMRT 治療過程中采用多個方向強度不均勻的射野照射 PTV。本文用 B 表示這些方向的射野集合，對于每個屬于該集合的子射束 b∈B 都被離散成一個矩陣。MLC 包括 M 對葉片，寬為 L，每對左右葉片的位置記作
和
，kb =
表示 MLC 對 b 的遮擋形成子野，則所有子野的集合即可表示為
。二進制數值
表示子野 kb 內第 m 行 j 列的值。
表示每個 kb 的權重。W 為劑量效應矩陣，代表單位強度射束在待照射體模內的劑量分布數據。劑量 D 是基于子野形狀和權重的線性函數。
目標函數的一般形式則可表示為：
 |
其中 l 是子目標函數個數，δl 是子目標函數的權重，fl（D）是子目標函數，對于同一個器官可以使用一個或多個子目標函數來控制其內的劑量分布。
放射治療中子目標函數可分為物理準則和生物準則。物理準則描述的是處方劑量與所計算劑量分布之間的差值，如平均劑量函數 fmean（D）、最小劑量函數 fmin（D）、最大劑量函數 fmax（D）和劑量-體積（dose-volume，DV）函數 fDVH（D）等[16]。具體表達形式如下：
 |
其中 Dmean 是 PTV 處方劑量或是 OAR 平均耐受劑量，Dmin 是 PTV 處方最小劑量，Dmax 是 OAR 耐受劑量。V 是器官內所有的體素點，Di 是器官第 i 個體素的劑量值。H（x）是階段函數。D1 是 DV 約束條件，
，即器官內接受到劑量大于 D1 的體素體積小于 V1。D2 是在當前劑量-體積直方圖（dose volume histogram，DVH）中體積為 V1 的體素接收到的實際劑量。
生物準則描述的是方案治療效果，是腫瘤或正常組織內的非線性放射生物效應，如等效均一劑量（equivalent uniform dose，EUD）函數 fgEUD（D）、腫瘤控制率（tumor control probability，TCP）函數 fTCP（D）和正常組織并發癥發生概率（normal tissue complication probability，NTCP）fNTCP（D）函數。這些準則大部分都是凸函數，非凸函數也可以轉換為等效凸函數[17-19]。
 |
其中 a 為反映劑量體積效應的參數因子。V0 是 PTV 內所有腫瘤克隆源性細胞數。
是在第 i 個體素單元內的腫瘤克隆源性細胞在受到 Di 劑量的照射后的存活率。
（x）是標準正態函數。D50（1）是體積 1 受照射時，由放射損傷引起的 NTCP 值為 50% 時所需的劑量。γ 是控制 NTCP 劑量效應曲線斜率的參數。
在這里，為了確保最終的子野可以被 MLC 照射，本文限制如下：① 對于同一行內的體素必須是連續的；② 相鄰葉片不能交錯。即：
 |
1.2 優化方法
子野形狀和權重是兩組相互耦合的參數，與劑量分布具有較為復雜的關系。為了簡化問題，本文分別采用兩種不同的優化方法對子野形狀和子野權重進行迭代優化。首先利用患者腫瘤特異的準直器截取束流方向（beam’s eye view，BEV）上的與靶區投影形狀一致的照射野作為子野初始形狀。然后優化子野形狀時，將子野中每行左右葉序位置作為變量，利用基于梯度信息的 SA 算法進行優化，最后將所有將個子野的權重作為變量利用 L-BFGS-B 算法進行優化。形狀優化和權重優化迭代進行。優化流程圖見圖 1。
圖1
GDAO 算法流程圖
Figure1.
Flow chart of the GDAO method
1.2.1 基于 SA 優化形狀
2002 年 Shepard 等[3]首次將 SA 算法應用于直接子野優化當中。SA 算法是基于 Monte Carlo 迭代求解策略的一種通用的隨機尋優算法，它在常規迭代算法的基礎上引入 Metropolis 準則，以一定的概率接受非最優解，在解大規模組合優化問題時展現出良好的性能，有較強的局部搜索能力。SA 算法最早由 Metropolis 等于 1953 年提出，1983 年 Kirkpatric 等將其應用于組合優化問題[20-21]。但隨著迭代次數的增加，SA 算法的收斂速度變慢，容易出現停滯，不利于在線優化，同時其算法性能與初始值有關而且參數敏感。
針對 SA 算法的以上不足，結合直接子野優化中葉片位置與劑量分布的關系，本文提出基于梯度信息的 SA 算法用于優化子野形狀。在文獻[2]和[3]中狀態產生函數采用服從寬度為 sigma 的高斯分布，本文根據葉片位置與最終劑量的關系提出一種含有梯度信息的狀態產生函數與含有梯度信息的狀態接受函數，確保每次搜索都是有效的，可快速收斂到最優解。
圖2
子野的形狀
Figure2.
Illustration of aperture shape
如圖 2 所示的子野形狀。在 m∈M 行時，如梯度為負數，意味著這行需要更多的照射，左右葉片需要向兩邊移動，保證更大的開口。相反，若梯度為正數，則意味要減少該行的照射，左右葉序要向中間移動，最終關閉。因此，本文用梯度決定每次迭代的變化步長，下一次迭代時
和
的值由式（5）確定。
 |
其中 gradtmt 是 m 行在迭代 t 時的梯度，rand（）為在[0，1]之間均勻分布的隨機數。
形成新的位置后，若
和
不滿足 MLC 約束，則直接拒絕新位置，并利用
和
已產生的位置代入式（5）重新計算新位置。否則，利用 Metropolis 準則對結果進行判斷：① 新位置的目標函數小于原有的最優函數值，直接接受新位置；② 大于原有的最優函數值，以式（6）中的概率接受新位置，否則拒絕。接受函數由子野 m 行在迭代 t 時的梯度 gradtmt（m∈M）和該行已接受使目標函數上升的解的次數 nacc 決定。當這行梯度值較大或接受次數較少，接受目標函數上升的解的概率較大；當梯度變小或接受次數變大，則接受目標函數上升的解的概率變小。隨著梯度的減少和接受上升解次數的增大，算法逐漸趨近于只能接受目標函數下降的解。
 |
其中 B 是初始概率，本文設置為 0.95。
綜上所述，給出利用改進后的 SA 優化子野形狀的具體步驟：
1）由 BEV 確定每個子野的初始形狀，得到 kb 中
和
的值，令
，并計算目標函數值
；
2）設置初始溫度 T（0），迭代次數 t = 1；
3）Do while T（t）>Tmin
4）for k = 1∶K
5）for m = 1∶M
6）對當前最優解 xbest 按照公式（5）計算產生新的解
；
7）如果 xnew 不滿足 MLC 約束，對
按照公式（5）產生新的解
；
8）計算新的目標函數值
，并計算目標函數值的增量
。
7）如果 ΔF < 0，則 xbest = xnew；執行 10）；
8）如果 ΔF > 0，則以公式（6）計算 P；
9）如果 rand（）<P，則 xbest = xnew；否則 xbest = xbest；
10）m = m + 1，轉 5）；
11）k = k + 1，轉 4）；
12）t = t + 1，轉 3）。
1.2.2 基于 L-BFGS-B 優化權重
所有子野形狀優化完畢后，本文利用梯度算法優化子野的權重。常見的梯度算法是 CG 算法，但一般每個治療計劃采用子野數目是 70～110[11]，有實驗數據表明當變量個數不超過 100 時，BFGS[22]法比 CG 效果好。L-BFGS[23]是 BFGS 的延伸，有效克服了 BFGS 占據內存過大的缺點，但其適用于求解大規模無約束優化問題，而直接子野優化中權重是帶有簡單邊界約束的。相較于 L-BFGS 算法，L-BFGS-B 算法具有求解約束優化問題的能力，在之前的工作中，我們已成功將其應用到 FMO 優化問題中[24]。因此，本文采用 L-BFGS-B 算法優化子野權重。具體步驟如下：
1）設定初始點
，定義初始有限內存矩陣 m 并令 q: = 0；
2）while 投影梯度
10–5；
3）計算廣義柯西點；
4）通過直接法計算搜索方向 dq；
5）在可行域內，沿著 dq 方向執行限定最大步長的線搜索來計算步長
，令
，計算
；
6）更新
、
、Lq 和 Rq，并且設置 θ =
；
7）令 q: = q + 1，轉 2）。
其中，
滿足充分下降條件，
是子空間上帶約束的二次型最小化的近似解，f 是函數值，gq 是梯度值，
是算法的修正對，
，
，
，
，
是投影梯度，Lq 和 Rq 是 m × m 的矩陣。
2 實驗結果與分析
本文將所提出的改進算法用于 10 例前列腺腫瘤病例，并與傳統 SA 算法進行了對比實驗。所有病例均來源于山西中醫學院附屬醫院腫瘤科前列腺腫瘤患者。X 射線計算機斷層成像（computed tomography，CT）數據由患者在仰臥位進行采集，然后在 CT 數據上勾畫靶區和 OAR。PTV 主要是前列腺（不考慮盆腔淋巴結），在臨床靶區（clinical target volume，CTV）的基礎上向后外擴 5 mm，在其余方向上外擴 10 mm。OAR 選取直腸壁和膀胱壁，手動圈定直腸輪廓和膀胱輪廓后分別外擴 5 mm。正常組織（Normal）是指將 PTV 外擴 5 cm 之外的所有組織。利用開源程序 CERR[25]讀取已勾畫好的患者 CT 數據，設置使用 5 個 6 MeV 射束分別以 34、100、180、260 和 324 機架角度照射靶目標并采用筆形束算法[26]計算劑量效應矩陣 W。實現算法的編程語言為 VC++（版本為 VS2013），計算機操作系統為 Windows10 專業版 64 位，處理器為英特爾 Core i3-4150@3.50GHz 雙核。
2.1 不同的優化目標與不同子野數目對優化結果的影響
在同一例前列腺腫瘤患者數據上測試了不同的優化目標與不同的子野數目對本文算法優化結果的影響。
首先目標函數設定為生物-物理混合目標函數[18]：OAR 的劑量分布利用 NTCP 函數
和
控制，靶區的劑量分布利用
和
控制，正常組織的劑量分布利用
控制。權重依次選取為 60、120、30、30 和 5。膀胱和直腸所使用的 NTCP 子目標函數中參數[27]分別為：
為 0.13，γ 為 0.11，D50（1）為 62 Gy；
為 0.13，γ 為 0.14，D50（1）為 80.7 Gy。靶區的最小劑量子目標函數的參數
為 74 Gy，PTV 的平均劑量子目標函數的參數
為 78 Gy。正常組織的 DVH 子目標函數約束條件為 V35 < 8%。分別設置子野數目為 25、35 和 45，優化結果如 圖 3 所示。圖 3a 顯示的是不同子野數目對應的生物目標函數值。其中實線代表采用 FMO 方法優化相同的目標函數未經過子野分割的最優結果。從圖中可以看出當子野數目為 45 時，優化結果最接近 FMO 最優結果。同時隨著迭代次數的增大，目標函數值趨于穩定。圖 3b、c、d 分別表示 PTV 的 DVH 曲線，從圖 3b 中可以看出，不同的子野數目在靶區 PTV 的劑量分布是相似的，但圖 3c 和圖 3d 的細節可以看出子野數目的越大，劑量分布在 PTV 的均勻性就越好。圖 3e 和圖 3f 顯示的是 OAR 的 DVH 曲線，可以看出不同的子野數目均很好地滿足了 Marks等[28]給出的臨床實踐指導標準，但子野數目較大時，OAR 收到的劑量相對少些。
圖3
生物目標下不同子野數目的優化結果
a. 目標函數值；b. PTV DVH 曲線；c. PTV 低劑量區；d. PTV 高劑量區；e. 直腸壁 DVH 曲線；f. 膀胱壁 DVH 曲線
Figure3. Results for NTCP-based cost functiona. objective function values; b. DVHs of PTV; c. low dose area of PTV; d. high dose area of PTV; e. DVHs of rectum; f. DVHs of bladder
保持其他參數不變，將目標函數設定為單一物理目標函數，OAR 膀胱和直腸的劑量分布分別由 3 個 DVH 子目標函數
和
控制。膀胱的 DV 約束條件分別為 V60 < 35%、 V37 < 25% 和 V75 < 15%，權重依次為 5、5 和 5。直腸的 DV 約束條件分別為 V50 < 50%、 V65 < 25% 和 V75 < 15%，權重依次為 5、8 和 30。同樣設置子野數目為 25、35 和 45，其優化結果如 圖 4 所示。圖 4a 顯示的是不同子野數目對應的物理目標函數值。其中實線代表采用 FMO 方法優化相同的目標函數未經過子野分割的最優結果。圖 4b、c、d 分別表示 PTV 的 DVH 曲線。圖 4e 和圖 4f 顯示的是 OAR 的 DVH 曲線，從圖 4 可以看出，本文算法對單一物理目標函數采用不同的子野數目得到的優化結果是相似的，都可找到滿足臨床要求的治療方案，而當子野數目為 45 時得到的目標函數值最接近 FMO 最優結果。
圖4
物理目標下不同子野數目的優化結果
a. 目標函數值；b. PTV DVH 曲線；c. PTV 低劑量區；d. PTV 高劑量區；e. 直腸壁 DVH 曲線；f. 膀胱壁 DVH 曲線
Figure4. Results for DVH-based cost functiona. objective function values; b. DVHs of PTV; c. low dose area of PTV; d. high dose area of PTV; e. DVHs of rectum; f. DVHs of bladder
以上結果表明，本文算法針對不同的目標函數，設置不同的子野數目均可得到滿足臨床需求的治療方案。
2.2 本文算法與傳統 SA 算法進行比較
將本文算法中優化子野形狀部分換為文獻[3]中的 SA 算法，進行對比試驗。目標函數選擇物理目標，子野個數設置為 45，其他參數與之前參數相同，優化結果對比如圖 5 所示。從圖 5a 中可以看出隨著迭代次數的增大本文算法趨于穩定，且目標函數值小于原 SA 算法的目標函數值。如圖 5b 所示，在每次迭代中，本文算法用于優化子野形狀總比原 SA 算法能更快地找到最小值。分析圖 5c，兩種算法靶區處方劑量覆蓋率均在 95% 以上，在靶區都體現出良好的均勻性；但在 OAR 部分本文算法所得到的 DVH 曲線明顯低于原 SA 算法所得到的 DVH 曲線。由此顯示，本文算法優化后的結果在提高靶區的劑量質量的同時，更好地保護了 OAR，治療增益比明顯提高。
圖5
兩種算法優化結果的對比
a. 迭代 50 次的目標函數值；b. 某次迭代中優化子野形狀得到的目標函數值；c. DVH 曲線
Figure5. Comparison of two algorithmsa. objective function values after 50 iterations; b. objective function values in one iteration; c. all DVHs
然后對 10 例前列腺腫瘤患者數據分別采用兩種算法進行優化。目標函數選為生物目標函數，權重設置見表 1，處方劑量等參數設置同上，子野個數設置為 45。優化后，靶區和 OAR 的平均劑量、最大劑量、最小劑量信息和算法運算時間如表 2 所示。采用 t 檢驗（檢驗水準為 0.05）分析 10 例患者分別用兩種算法優化后的結果顯示，GDAO 算法由于使用了梯度信息，優化效率比文獻[3]中傳統 SA 算法效率更高（P = 0.04），時間減少了 15.90%。同時，雖然兩者得到的優化方案非常相似（P > 0.05），但靶區最低劑量提高了 0.29%，最高劑量降低了 0.45%，OAR 膀胱最高劑量降低了 0.25%，OAR 直腸最高劑量降低了 0.09%。說明采用本文算法可在短時間內得到不差于原 SA 算法的優化結果，并且在保證 PTV 的劑量分布的同時，更好地控制了 OAR 的高劑量分布。美中不足的是新算法增加了部分低劑量分布（< 50 Gy），但低劑量分布對由放射治療而產生的正常組織并發癥的影響是有限的，所以并不影響放療方案質量。
表1
10 例前列腺癌癥患者所使用的生物目標函數權重
Table1.
Weights of cost function for the ten patients
表2
兩種算法的優化結果
Table2.
The optimization results of the two algorithms
3 結語
本文提出了一種 GDAO 優化方法。該方法首先采用含有梯度信息的 SA 算法對子野形狀進行優化，然后采用 L-BFGS-B 算法優化子野權重。該方法克服了原 SA 算法收斂慢、易停滯、全局搜索能力低的缺點，可以快速得到質量更好的最優解。實驗結果表明，可在短時間內得到滿足臨床要求并可直接實施照射的治療方案，具有較好的臨床實用價值。
引言
調強放射治療[1]（intensity-modulated radiotherapy，IMRT）是當前最常用的臨床放射治療技術。它采用多個方向強度不均勻的的射野照射靶區（planning target volume，PTV），產生與 PTV 高度適形的劑量分布，同時減少對危及器官（organ at risk，OAR）和正常組織的損傷，明顯提高了治療精度。世界一流的治療計劃系統（treatment planning system，TPS）大部分都提供了 IMRT 功能，如 Elekta Inc.公司的 PrecisePLAN 系統、Philips Medical Systems 公司的 P3IMRT 系統和 Varian Medical Systems 公司的 Eclipse Integrated Treatment Planning 系統。傳統的 IMRT 一般采用“兩步法”：先優化注量圖（fluence map optimization，FMO）得到最優強度分布，再將最優強度分布分割成多葉準直器（multileaf collimator，MLC）序列。“兩步法”的缺點在于優化注量圖時并未考慮 MLC 約束，射野分割后得到的強度分布與最優強度分布之間的差別可達 5.7%[2]，這種差別降低了劑量分布的適形度，影響了放射治療質量。為了克服“兩步法”照射強度分布與最優強度分布之間不一致的缺點，直接子野優化（direct aperture optimization，DAO）[3]通過逆向計劃直接對子野的形狀和權重進行優化，得到滿足 MLC 約束的子野，降低了治療的不確定性，從而成為近幾年較熱門的 IMRT 優化方法。
DAO 屬于 NP-hard 的問題，其中子野形狀和權重是兩組相互耦合的參數，每個特定的子野形狀都要優化其對應的權重。為了解決這個問題，Shepard 等[3]采用模擬退火算法（simulated annealing，SA）對子野形狀和權重同時進行優化，他們證明采用 DAO 技術可以得到適形度較高的劑量分布。隨后 Earl 等[4]對這種方法進行了改進，只使用矩形的子野形狀來實現各射野方向的強度調強。李永杰等[5]以二維編碼的方式來表達每個子野形狀，提出了遺傳算法（genetic algorithm，GA）優化子野形狀以及共軛梯度算法（conjugate gradient，CG）優化子野權重的算法。Cotrutz 等[6]也采用遺傳算法直接對子野進行了優化。Preciado-Walters 等[7]、Romeijn 等[8]、Men 等[9]、Salari 等[10]、Bednarz 等[11]，以及 Salari 等[12]都利用列生成的方法直接生成子野形狀。Cassioli 等[13]提出一種利用梯度算法優化子野形狀的方法，但該方法只有在葉序變動范圍小的時候才有效，因此他們還特為算法設置了信任區域。王捷等[14]則先采用 CG 優化子野權重，然后再利用枚舉法優化子野形狀。以上方法均取得了較好的方案優化結果，但是仍存在一些問題。在子野形狀優化算法中，文獻[3-6]采用的是隨機搜索算法，優化過程中有一定概率跳出局部最小值，但不保證每次的優化結果為全局最優；文獻[13]采用局部梯度的方法，由于葉序位置與函數梯度并沒有嚴格的線性關系，所以不能實現全局最優；文獻[7-12]使用圖論的方法不能完全描述目標函數與 MLC 葉片位置的復雜關系，影響了放療方案的質量。因此如何更好地優化子野形狀，是一個亟待解決的問題。同時，在子野權重優化方面，也需有所改進。
綜上所述，本文提出一種基于梯度信息的直接子野優化方法（gradient-based direct aperture optimization，GDAO）用于 IMRT 逆向計劃。在 GDAO 中分別采用兩種不同的優化方法對子野形狀和子野權重進行迭代優化。首先，為了增加子野葉片位置的精度，采用基于梯度信息的 SA 算法優化子野形狀，根據葉片位置與最終劑量的關系提出一種含有梯度信息的狀態產生函數與含有梯度信息的狀態接受函數，并在優化同時充分考慮 MLC 葉片間的約束條件。最后再利用計算量少、迭代代價低、收斂快且穩定的梯度類具有求解大規模約束優化問題能力的帶約束最小存儲擬牛頓算法[15]（limited-memory BFGS for bound-constrained，L-BFGS-B）優化子野權重。
1 本文方法
1.1 目標函數
在 IMRT 治療過程中采用多個方向強度不均勻的射野照射 PTV。本文用 B 表示這些方向的射野集合，對于每個屬于該集合的子射束 b∈B 都被離散成一個矩陣。MLC 包括 M 對葉片，寬為 L，每對左右葉片的位置記作
和
，kb =
表示 MLC 對 b 的遮擋形成子野，則所有子野的集合即可表示為
。二進制數值
表示子野 kb 內第 m 行 j 列的值。
表示每個 kb 的權重。W 為劑量效應矩陣，代表單位強度射束在待照射體模內的劑量分布數據。劑量 D 是基于子野形狀和權重的線性函數。
目標函數的一般形式則可表示為：
|
其中 l 是子目標函數個數，δl 是子目標函數的權重，fl（D）是子目標函數，對于同一個器官可以使用一個或多個子目標函數來控制其內的劑量分布。
放射治療中子目標函數可分為物理準則和生物準則。物理準則描述的是處方劑量與所計算劑量分布之間的差值，如平均劑量函數 fmean（D）、最小劑量函數 fmin（D）、最大劑量函數 fmax（D）和劑量-體積（dose-volume，DV）函數 fDVH（D）等[16]。具體表達形式如下：
|
其中 Dmean 是 PTV 處方劑量或是 OAR 平均耐受劑量，Dmin 是 PTV 處方最小劑量，Dmax 是 OAR 耐受劑量。V 是器官內所有的體素點，Di 是器官第 i 個體素的劑量值。H（x）是階段函數。D1 是 DV 約束條件，
，即器官內接受到劑量大于 D1 的體素體積小于 V1。D2 是在當前劑量-體積直方圖（dose volume histogram，DVH）中體積為 V1 的體素接收到的實際劑量。
生物準則描述的是方案治療效果，是腫瘤或正常組織內的非線性放射生物效應，如等效均一劑量（equivalent uniform dose，EUD）函數 fgEUD（D）、腫瘤控制率（tumor control probability，TCP）函數 fTCP（D）和正常組織并發癥發生概率（normal tissue complication probability，NTCP）fNTCP（D）函數。這些準則大部分都是凸函數，非凸函數也可以轉換為等效凸函數[17-19]。
|
其中 a 為反映劑量體積效應的參數因子。V0 是 PTV 內所有腫瘤克隆源性細胞數。
是在第 i 個體素單元內的腫瘤克隆源性細胞在受到 Di 劑量的照射后的存活率。
（x）是標準正態函數。D50（1）是體積 1 受照射時，由放射損傷引起的 NTCP 值為 50% 時所需的劑量。γ 是控制 NTCP 劑量效應曲線斜率的參數。
在這里，為了確保最終的子野可以被 MLC 照射，本文限制如下：① 對于同一行內的體素必須是連續的；② 相鄰葉片不能交錯。即：
|
1.2 優化方法
子野形狀和權重是兩組相互耦合的參數，與劑量分布具有較為復雜的關系。為了簡化問題，本文分別采用兩種不同的優化方法對子野形狀和子野權重進行迭代優化。首先利用患者腫瘤特異的準直器截取束流方向（beam’s eye view，BEV）上的與靶區投影形狀一致的照射野作為子野初始形狀。然后優化子野形狀時，將子野中每行左右葉序位置作為變量，利用基于梯度信息的 SA 算法進行優化，最后將所有將個子野的權重作為變量利用 L-BFGS-B 算法進行優化。形狀優化和權重優化迭代進行。優化流程圖見圖 1。
圖1
GDAO 算法流程圖
Figure1.
Flow chart of the GDAO method
1.2.1 基于 SA 優化形狀
2002 年 Shepard 等[3]首次將 SA 算法應用于直接子野優化當中。SA 算法是基于 Monte Carlo 迭代求解策略的一種通用的隨機尋優算法，它在常規迭代算法的基礎上引入 Metropolis 準則，以一定的概率接受非最優解，在解大規模組合優化問題時展現出良好的性能，有較強的局部搜索能力。SA 算法最早由 Metropolis 等于 1953 年提出，1983 年 Kirkpatric 等將其應用于組合優化問題[20-21]。但隨著迭代次數的增加，SA 算法的收斂速度變慢，容易出現停滯，不利于在線優化，同時其算法性能與初始值有關而且參數敏感。
針對 SA 算法的以上不足，結合直接子野優化中葉片位置與劑量分布的關系，本文提出基于梯度信息的 SA 算法用于優化子野形狀。在文獻[2]和[3]中狀態產生函數采用服從寬度為 sigma 的高斯分布，本文根據葉片位置與最終劑量的關系提出一種含有梯度信息的狀態產生函數與含有梯度信息的狀態接受函數，確保每次搜索都是有效的，可快速收斂到最優解。
圖2
子野的形狀
Figure2.
Illustration of aperture shape
如圖 2 所示的子野形狀。在 m∈M 行時，如梯度為負數，意味著這行需要更多的照射，左右葉片需要向兩邊移動，保證更大的開口。相反，若梯度為正數，則意味要減少該行的照射，左右葉序要向中間移動，最終關閉。因此，本文用梯度決定每次迭代的變化步長，下一次迭代時
和
的值由式（5）確定。
|
其中 gradtmt 是 m 行在迭代 t 時的梯度，rand（）為在[0，1]之間均勻分布的隨機數。
形成新的位置后，若
和
不滿足 MLC 約束，則直接拒絕新位置，并利用
和
已產生的位置代入式（5）重新計算新位置。否則，利用 Metropolis 準則對結果進行判斷：① 新位置的目標函數小于原有的最優函數值，直接接受新位置；② 大于原有的最優函數值，以式（6）中的概率接受新位置，否則拒絕。接受函數由子野 m 行在迭代 t 時的梯度 gradtmt（m∈M）和該行已接受使目標函數上升的解的次數 nacc 決定。當這行梯度值較大或接受次數較少，接受目標函數上升的解的概率較大；當梯度變小或接受次數變大，則接受目標函數上升的解的概率變小。隨著梯度的減少和接受上升解次數的增大，算法逐漸趨近于只能接受目標函數下降的解。
|
其中 B 是初始概率，本文設置為 0.95。
綜上所述，給出利用改進后的 SA 優化子野形狀的具體步驟：
1）由 BEV 確定每個子野的初始形狀，得到 kb 中
和
的值，令
，并計算目標函數值
；
2）設置初始溫度 T（0），迭代次數 t = 1；
3）Do while T（t）>Tmin
4）for k = 1∶K
5）for m = 1∶M
6）對當前最優解 xbest 按照公式（5）計算產生新的解
；
7）如果 xnew 不滿足 MLC 約束，對
按照公式（5）產生新的解
；
8）計算新的目標函數值
，并計算目標函數值的增量
。
7）如果 ΔF < 0，則 xbest = xnew；執行 10）；
8）如果 ΔF > 0，則以公式（6）計算 P；
9）如果 rand（）<P，則 xbest = xnew；否則 xbest = xbest；
10）m = m + 1，轉 5）；
11）k = k + 1，轉 4）；
12）t = t + 1，轉 3）。
1.2.2 基于 L-BFGS-B 優化權重
所有子野形狀優化完畢后，本文利用梯度算法優化子野的權重。常見的梯度算法是 CG 算法，但一般每個治療計劃采用子野數目是 70～110[11]，有實驗數據表明當變量個數不超過 100 時，BFGS[22]法比 CG 效果好。L-BFGS[23]是 BFGS 的延伸，有效克服了 BFGS 占據內存過大的缺點，但其適用于求解大規模無約束優化問題，而直接子野優化中權重是帶有簡單邊界約束的。相較于 L-BFGS 算法，L-BFGS-B 算法具有求解約束優化問題的能力，在之前的工作中，我們已成功將其應用到 FMO 優化問題中[24]。因此，本文采用 L-BFGS-B 算法優化子野權重。具體步驟如下：
1）設定初始點
，定義初始有限內存矩陣 m 并令 q: = 0；
2）while 投影梯度
10–5；
3）計算廣義柯西點；
4）通過直接法計算搜索方向 dq；
5）在可行域內，沿著 dq 方向執行限定最大步長的線搜索來計算步長
，令
，計算
；
6）更新
、
、Lq 和 Rq，并且設置 θ =
；
7）令 q: = q + 1，轉 2）。
其中，
滿足充分下降條件，
是子空間上帶約束的二次型最小化的近似解，f 是函數值，gq 是梯度值，
是算法的修正對，
，
，
，
，
是投影梯度，Lq 和 Rq 是 m × m 的矩陣。
2 實驗結果與分析
本文將所提出的改進算法用于 10 例前列腺腫瘤病例，并與傳統 SA 算法進行了對比實驗。所有病例均來源于山西中醫學院附屬醫院腫瘤科前列腺腫瘤患者。X 射線計算機斷層成像（computed tomography，CT）數據由患者在仰臥位進行采集，然后在 CT 數據上勾畫靶區和 OAR。PTV 主要是前列腺（不考慮盆腔淋巴結），在臨床靶區（clinical target volume，CTV）的基礎上向后外擴 5 mm，在其余方向上外擴 10 mm。OAR 選取直腸壁和膀胱壁，手動圈定直腸輪廓和膀胱輪廓后分別外擴 5 mm。正常組織（Normal）是指將 PTV 外擴 5 cm 之外的所有組織。利用開源程序 CERR[25]讀取已勾畫好的患者 CT 數據，設置使用 5 個 6 MeV 射束分別以 34、100、180、260 和 324 機架角度照射靶目標并采用筆形束算法[26]計算劑量效應矩陣 W。實現算法的編程語言為 VC++（版本為 VS2013），計算機操作系統為 Windows10 專業版 64 位，處理器為英特爾 Core i3-4150@3.50GHz 雙核。
2.1 不同的優化目標與不同子野數目對優化結果的影響
在同一例前列腺腫瘤患者數據上測試了不同的優化目標與不同的子野數目對本文算法優化結果的影響。
首先目標函數設定為生物-物理混合目標函數[18]：OAR 的劑量分布利用 NTCP 函數
和
控制，靶區的劑量分布利用
和
控制，正常組織的劑量分布利用
控制。權重依次選取為 60、120、30、30 和 5。膀胱和直腸所使用的 NTCP 子目標函數中參數[27]分別為：
為 0.13，γ 為 0.11，D50（1）為 62 Gy；
為 0.13，γ 為 0.14，D50（1）為 80.7 Gy。靶區的最小劑量子目標函數的參數
為 74 Gy，PTV 的平均劑量子目標函數的參數
為 78 Gy。正常組織的 DVH 子目標函數約束條件為 V35 < 8%。分別設置子野數目為 25、35 和 45，優化結果如 圖 3 所示。圖 3a 顯示的是不同子野數目對應的生物目標函數值。其中實線代表采用 FMO 方法優化相同的目標函數未經過子野分割的最優結果。從圖中可以看出當子野數目為 45 時，優化結果最接近 FMO 最優結果。同時隨著迭代次數的增大，目標函數值趨于穩定。圖 3b、c、d 分別表示 PTV 的 DVH 曲線，從圖 3b 中可以看出，不同的子野數目在靶區 PTV 的劑量分布是相似的，但圖 3c 和圖 3d 的細節可以看出子野數目的越大，劑量分布在 PTV 的均勻性就越好。圖 3e 和圖 3f 顯示的是 OAR 的 DVH 曲線，可以看出不同的子野數目均很好地滿足了 Marks等[28]給出的臨床實踐指導標準，但子野數目較大時，OAR 收到的劑量相對少些。
圖3
生物目標下不同子野數目的優化結果
a. 目標函數值；b. PTV DVH 曲線；c. PTV 低劑量區；d. PTV 高劑量區；e. 直腸壁 DVH 曲線；f. 膀胱壁 DVH 曲線
Figure3. Results for NTCP-based cost functiona. objective function values; b. DVHs of PTV; c. low dose area of PTV; d. high dose area of PTV; e. DVHs of rectum; f. DVHs of bladder
保持其他參數不變，將目標函數設定為單一物理目標函數，OAR 膀胱和直腸的劑量分布分別由 3 個 DVH 子目標函數
和
控制。膀胱的 DV 約束條件分別為 V60 < 35%、 V37 < 25% 和 V75 < 15%，權重依次為 5、5 和 5。直腸的 DV 約束條件分別為 V50 < 50%、 V65 < 25% 和 V75 < 15%，權重依次為 5、8 和 30。同樣設置子野數目為 25、35 和 45，其優化結果如 圖 4 所示。圖 4a 顯示的是不同子野數目對應的物理目標函數值。其中實線代表采用 FMO 方法優化相同的目標函數未經過子野分割的最優結果。圖 4b、c、d 分別表示 PTV 的 DVH 曲線。圖 4e 和圖 4f 顯示的是 OAR 的 DVH 曲線，從圖 4 可以看出，本文算法對單一物理目標函數采用不同的子野數目得到的優化結果是相似的，都可找到滿足臨床要求的治療方案，而當子野數目為 45 時得到的目標函數值最接近 FMO 最優結果。
圖4
物理目標下不同子野數目的優化結果
a. 目標函數值；b. PTV DVH 曲線；c. PTV 低劑量區；d. PTV 高劑量區；e. 直腸壁 DVH 曲線；f. 膀胱壁 DVH 曲線
Figure4. Results for DVH-based cost functiona. objective function values; b. DVHs of PTV; c. low dose area of PTV; d. high dose area of PTV; e. DVHs of rectum; f. DVHs of bladder
以上結果表明，本文算法針對不同的目標函數，設置不同的子野數目均可得到滿足臨床需求的治療方案。
2.2 本文算法與傳統 SA 算法進行比較
將本文算法中優化子野形狀部分換為文獻[3]中的 SA 算法，進行對比試驗。目標函數選擇物理目標，子野個數設置為 45，其他參數與之前參數相同，優化結果對比如圖 5 所示。從圖 5a 中可以看出隨著迭代次數的增大本文算法趨于穩定，且目標函數值小于原 SA 算法的目標函數值。如圖 5b 所示，在每次迭代中，本文算法用于優化子野形狀總比原 SA 算法能更快地找到最小值。分析圖 5c，兩種算法靶區處方劑量覆蓋率均在 95% 以上，在靶區都體現出良好的均勻性；但在 OAR 部分本文算法所得到的 DVH 曲線明顯低于原 SA 算法所得到的 DVH 曲線。由此顯示，本文算法優化后的結果在提高靶區的劑量質量的同時，更好地保護了 OAR，治療增益比明顯提高。
圖5
兩種算法優化結果的對比
a. 迭代 50 次的目標函數值；b. 某次迭代中優化子野形狀得到的目標函數值；c. DVH 曲線
Figure5. Comparison of two algorithmsa. objective function values after 50 iterations; b. objective function values in one iteration; c. all DVHs
然后對 10 例前列腺腫瘤患者數據分別采用兩種算法進行優化。目標函數選為生物目標函數，權重設置見表 1，處方劑量等參數設置同上，子野個數設置為 45。優化后，靶區和 OAR 的平均劑量、最大劑量、最小劑量信息和算法運算時間如表 2 所示。采用 t 檢驗（檢驗水準為 0.05）分析 10 例患者分別用兩種算法優化后的結果顯示，GDAO 算法由于使用了梯度信息，優化效率比文獻[3]中傳統 SA 算法效率更高（P = 0.04），時間減少了 15.90%。同時，雖然兩者得到的優化方案非常相似（P > 0.05），但靶區最低劑量提高了 0.29%，最高劑量降低了 0.45%，OAR 膀胱最高劑量降低了 0.25%，OAR 直腸最高劑量降低了 0.09%。說明采用本文算法可在短時間內得到不差于原 SA 算法的優化結果，并且在保證 PTV 的劑量分布的同時，更好地控制了 OAR 的高劑量分布。美中不足的是新算法增加了部分低劑量分布（< 50 Gy），但低劑量分布對由放射治療而產生的正常組織并發癥的影響是有限的，所以并不影響放療方案質量。
表1
10 例前列腺癌癥患者所使用的生物目標函數權重
Table1.
Weights of cost function for the ten patients
表2
兩種算法的優化結果
Table2.
The optimization results of the two algorithms
3 結語
本文提出了一種 GDAO 優化方法。該方法首先采用含有梯度信息的 SA 算法對子野形狀進行優化，然后采用 L-BFGS-B 算法優化子野權重。該方法克服了原 SA 算法收斂慢、易停滯、全局搜索能力低的缺點，可以快速得到質量更好的最優解。實驗結果表明，可在短時間內得到滿足臨床要求并可直接實施照射的治療方案，具有較好的臨床實用價值。
