基于啟發式信息的非凸放療規劃模型的求解方法

1.淮北師范大學 計算機科學與技術學院，安徽 淮北 235000

2.南陽師范學院 計算機與信息技術學院，河南 南陽 473061

3.連云港市第二人民醫院 醫學影像科，江蘇 連云港 222005

1.淮北師范大學 計算機科學與技術學院，安徽 淮北 235000

2.南陽師范學院 計算機與信息技術學院，河南 南陽 473061

3.連云港市第二人民醫院 醫學影像科，江蘇 連云港 222005

1 引言

調強放射治療是目前對惡性腫瘤進行治療的有效手段之一，這種治療手段通過可調的照射角度、可調的各照射野射線強度，在待治療區域上形成可調的放射線劑量分布，保證腫瘤殺滅區域的劑量達到治療水平，而非腫瘤區域的放射劑量控制在發生并發癥的水平以下，危險區域的放射劑量則控制在最低水平。這樣強度可調的三維形狀適應性很強的放療手段改變了過去傳統放療并發癥較高的弊端，已經在歐美等發達國家對頭頸部癌、前列腺癌、肺癌、腦癌等癌癥的治療中廣泛使用，目前國內的一些大中城市的醫院也在推廣普及這種治療技術。

調強放療醫療器械的生產廠商主要還是國外的知名公司，如通用、瓦里安等。該類器械的一個關鍵技術就是如何從醫生設定計劃逆向求解得到放療規劃結果，即前期通過醫生對待治療區域器官、照射靶區的勾勒后，形成劑量要求和劑量約束條件，送入放療逆向求解系統中進行優化求解，最終得到硬件設備可以施行的照射方案。其中逆向優化過程的一個關鍵環節就是得到各照射野的照射分布圖[1-3]。一個好的照射分布圖應該具備以下兩個方面的要求。（1）待照射靶區的劑量接收值應與醫生設定的計劃劑量值盡可能相等，且保證劑量分布的一致性。（2）待照射靶區周邊組織的劑量接收值應越低越好，從而在靶區三維輪廓上形成劑量快速下降的分布趨勢，特別是一些危險的串行組織需要嚴格限制劑量上限，而一些并行組織則可以對劑量約束適當放松。

目前已經有大量的文獻正關注于物理劑量模型的建立和相應的求解算法[4-5]，考慮到求解算法的高效性要求，目前的逆向計劃大多轉化為線性或二次規劃模型[6-8]。在這些模型中，劑量逼近程度一般采用實際劑量與設定劑量的絕對差值或絕對差值的平方進行描述；并廣泛采用劑量體積曲線的多點限制模式，即通過劑量體積圖曲線的調節指導放療優化，讓正常組織滿足一定的劑量體積比要求的情況下，盡可能地讓待治療區域達到放射劑量要求。

非常可惜的是，帶有劑量體積限制的數學規劃模型為非凸規劃問題[9-11]，為了求解這樣一類較難的數學問題，近幾年來國內外科學家提出了許多求解方法，取得了一定的成效，但離獲得最優解仍有較大距離。

Cotrutz等[12]使用基于體素權重系數的方法來幫助求解，在他們的方法中，首先對線性規劃問題或二次規劃問題進行初始求解，并獲得初始劑量分布，然后按該劑量分布，根據低劑量體素更加容易抑制為更低劑量、高劑量體素更加容易給予更高劑量這樣的啟發式信息為指導，重新對體素的權重進行分配，再次進行優化；重復上述過程直到達到優化限定條件。這樣的啟發式信息有一定的直觀依據，但是形而上的想法缺乏科學依據。本文提出一種新的更加科學的啟發式方法對上述非凸規劃問題進行求解，并在四個病例上進行了相關實驗。實驗結果證實，新的啟發式方法能夠獲得更優的逆向計劃。

2 方法

2.1 強度照射分布圖優化問題

整個待治療區域的人體組織被分為兩個部分，且用體素進行離散化，一部分為待照射靶區T={T1，T2，…，Tu}，相應的體素數目為…，STu，另一部分為正常組織和危險器官組織N={N1，N2，…，Nv}，相應的體素數目為…，。照射分布圖用像素進行離散化，按照照射野的順序以及分布圖逐行的順序，用列向量X進行描述，FT=為靶區劑量作用矩陣，為正常組織和危險器官組織的劑量作用矩陣。任意一個元素 fi，j表示第j個像素采用單位強度照射時對第i個體素所產生的劑量沉積值。 DT=FT·X和 DN=FN·X代表整個待照射區域在強度分布圖X下的三維劑量分布。設定醫生指定的劑量為，…，，正常組織以及危險區域組織的劑量上限為，…，

給出物理劑量模型所廣泛采用的線性約束二次規劃模型如下：

2.2 劑量體積曲線限制條件

放射科臨床醫生常用的一個概念為劑量體積曲線，這個曲線能夠從宏觀上描述某個器官的劑量分布生物學特征。劑量體積曲線的橫軸為劑量值，縱軸為體積百分比，曲線上的每一個點表示某個器官中低于某個劑量的組織體積占整個器官組織體積的百分比。因此該曲線為一條遞減曲線。

而相應的劑量體積曲線限制即要求曲線在某個點的上方或某個點的下方。以正常組織為例，劑量體積曲線限制條件一般限定曲線在一個點或幾個點的下方。一般描述為某個器官組織接受? Gy劑量或更高劑量的體素數目不超過該器官組織總體素數目的η%，其中“Gy”為放射劑量單位。直觀上，即要求曲線在點( ) ?，η的下方。

圖1 劑量體積曲線及約束點

設正常組織和危險器官組織 N1，N2，…，Nv有劑量體積限定條件 DVC1，DVC2，…，DVCv，其中為劑量體積限定條件 DVCi中的約束點個數，為相應的約束點。由此形成以下的帶有劑量體積約束的數學規劃問題：

對于臨床醫生常用的一個概念被引入到逆向計劃中卻帶來了求解過程的極大困難，因為上述模型為組合優化模型。

2.3 正規化空間內的空間距離啟發信息

上述組合優化模型的求解為非凸規劃問題的求解，獲得該問題的最優解是NP難問題。目前國外一些參考文獻的求解方法是采用基于體素權重系數的問題轉換求解思路，即通過對各體素設定一定的權重，從而將上述問題轉換為凸規劃問題進行求解，當然，還有一些思路則回避劑量體積曲線約束，考慮采用其他約束取代這類非凸約束，從而轉換為凸規劃問題。對于后者，認為約束條件的轉換已經基本違背了臨床醫生的劑量設定初衷，因此本文不予考慮。對于問題轉換思路，本文予以采納，即將上述問題轉換為凸規劃從而獲得轉換問題的最優解，但是對于體素排序的思路，本文提出新的啟發式想法，即從原來的劑量高低排序方式轉換為二次規劃正規化空間內的空間距離遠近排序方式。具體算法流程如圖2所示。

圖2 優化算法流程圖

根據Cholesky分解，Q=RTR，帶入 x=R-1y，上述模型轉換為：

對于線性約束二次規劃問題：

設(R-1)Tc=e，BR-1=G，上述線性約束二次規劃問題有以下的正規化形式：

下面通過一個簡單的示例來對劑量高低排序方式以及二次規劃正規化空間內的空間距離遠近排序方式進行比較。

該示例只有一個規劃靶區，并為一個凹的形狀，另一個危險組織區域為圓形，并一定程度地被規劃靶區包圍。照射方向為五個均勻分布的方向，靶區設定的劑量為單位劑量，危險組織區域的最大劑量不能超過0.5，正常組織劑量不超過0.8，危險組織區域的劑量體積曲線要求接受高于0.3劑量的體積大小不超過50%（如圖3所示）。

圖3 示例初始劑量分布示意圖

求解了不帶有劑量體積曲線限制條件的規劃問題，并按照劑量高低排序方式以及二次規劃正規化空間內的空間距離遠近排序方式對危險組織區域進行了排序，排序結果如圖4所示，顯然，基于二次規劃正規化空間內的空間距離遠近排序方式的排序結果更加自然，遠離靶區的體素具有更高的優先級被納入線性約束中去，從而達到所要求的劑量體積曲線約束。

圖4 危險組織區域的正規化空間內空間距離和劑量排序比較示意圖

2.4 實驗設定

為了驗證所提方法的有效性，使用四個臨床實驗數據進行了測試，劑量設定值見表1～4。本文中實驗強度矩陣的計算基于QIB劑量計算引擎。所有實驗均在Matlab計算平臺下進行。

表1 頭頸部測試例的劑量設定值

表2 前列腺測試例的劑量設定值

表3 肺部測試例的劑量設定值

3 實驗結果和分析

表4 口咽部測試例的劑量設定值

表5 優化結果比較

圖5 頭頸部測試例的劑量體積曲線

在四個測試病例上，對基于劑量高低排序方式的啟發信息以及基于正規化空間內的空間距離遠近排序方式的啟發信息進行了實驗對比，表5中給出了兩種方法所得到的目標函數值。圖5～8中給出了兩種方法得到的靶區各照射目標以及非靶區各器官的劑量體積曲線比較結果（標有實心圖例的實線是本文方法的結果，標有對應的虛心圖例的虛線是基于劑量高低排序方式方法的結果）。圖9～12中還給出了兩種方法得到的斷層劑量分布圖的比較結果。實驗結果顯示，在達到同樣的劑量體積曲線限制條件的約束情況下，采用提出的啟發信息，病例一的目標函數值降低了（（26.9－44.9）/44.9=－40%，病例二降低了（54.5－ 84.0）/84.0=－35%，病例四降低了（17.2－23.5）/23.5=－27%），除病例三以外降低幅度均較大。通過比較靶區各目標的劑量體積曲線，可以看到實線比虛線下降得更為陡峭，也就是說采用本文提出的啟發信息，靶區的劑量分布更加均勻，且更接近于所設定的處方劑量值，通過斷層劑量分布圖可以明顯地得到這個結論。

圖6 前列腺測試例的劑量體積曲線

圖7 肺部測試例的劑量體積曲線

圖8 口咽部測試例的劑量體積曲線

圖9 頭頸部測試例的某斷層劑量分布圖

圖10 前列腺測試例的某斷層劑量分布圖

圖11 肺部測試例的某斷層劑量分布圖

圖12 口咽部測試例的某斷層劑量分布圖

4 結束語

本文對調強放療數學規劃問題中的非凸規劃提出一種新的啟發式思想，利用這種思想得到的非凸規劃問題求解方法能夠獲得更優的優化結果。在四個病例上的實驗結果證實了所提新方法的有效性。此外，這種啟發思想能夠在其他非凸數學規劃問題的求解中同樣得到廣泛應用。

[1]Bortfeld T.IMRT：a review and preview[J].Physics in Medicine and Biology，2006，51：363-379.

[2]Hatano K，Araki H，Sakai M，et al.Current status of Intensity-ModulatedRadiation Therapy（IMRT）[J].International Journal of Clinical Oncology，2007，12：408-415.

[3]Yu C X，Amies C J，Svatos M.Planning and delivery of intensity-modulated radiation therapy[J].Medical Physics，2008，25：5233-5241.

[4]Shepard D M，Ferris M C，Olivera G H，et al.Optimizing the delivery of radiation therapy to cancer patients[J].Siam Review，1999，41：721-744.

[5]Lim J，Ferris M C，Wright S J，et al.An optimization framework for conformal radiation treatment planning[J].INFORMS Journal on Computing，2007，19：366-380.

[6]Morrill S M，Lane R G，Wong J A，et al.Dose-volume considerations with linear programming optimization[J].Medical Physics，1991，18：1201-1210.

[7]Breedveld S，Storchi P R M，Keijzer M，et al.Fast，multiple optimizations of quadratic dose objective functions in IMRT[J]. Physics in Medicine and Biology，2006，51：3569-3579.

[8]Censor Y，Ben-Israel A，Xiao Y，et al.On linear infeasibility arising in intensity-modulated radiation therapy inverse planning[J].LinearAlgebra and ItsApplications，2008，428：1406-1420.

[9]Deasy J O.Multiple local minima in radiotherapy optimization problems with dose-volume constraints[J].Medical Physics，1997，24：1157-1161.

[10]Wu Q，Mohan R.Multiple local minima in IMRT optimization based on dose-volume criteria[J].MedicalPhysics，2002，29：1514-1527.

[11]Wu C，Jeraj R，Mackie T R.The method of intercepts in parameter space for the analysis of local minima caused by dose-volume constraints[J].Physics in Medicine and Biology，2003，48：149-157.

[12]Cotrutz C，Xing L.Using voxel-dependent importance factors for interactive DVH-based dose optimization[J].Physics in Medicine and Biology，2002，47：1659-1669.

基于啟發式信息的非凸放療規劃模型的求解方法

張棟冰1，蘭義華2，萬金鑫3

ZHANG Dongbing1,LAN Yihua2,WAN Jinxin3

1.School of Computer Science and Technology,Huaibei Normal University,Huaibei,Anhui 235000,China
2.School of Computer and Information Technology,Nanyang Normal University,Nanyang,Henan 473061,China
3.Medical Imaging Department,Lianyungang Second People’s Hospital,Lianyungang,Jiangsu 222005,China

For the inverse planning process of intensity-modulated radiotherapy—the non-convex mathematical programming for the fluence map optimization with dose volume constraints,this paper presents a novel and scientific heuristic information—distance priority value in normalized space.Compared to the traditional dose ordering heuristic information,this new method leads to better solutions.A phantom example and four test cases show the effectiveness of the proposed method.

non-convex fractional programming;heuristic algorithm;intensity-modulated radiation therapy technology;linear constrained quadratic programming

針對調強放療逆向優化過程中的關鍵環節——各照射野的強度照射分布圖在帶有劑量體積曲線限制條件下的非凸數學規劃問題，提出了一種新穎的更加科學的啟發式信息——正規化空間內的空間距離排序值。與傳統的劑量排序啟發式信息相比較，新方法可以得到更好的解。一個簡單示例和四個測試病例表明了該方法的有效性。

非凸數學規劃；啟發式求解；調強放療；線性約束二次規劃

A

TP391

10.3778/j.issn.1002-8331.1302-0021

ZHANG Dongbing,LAN Yihua,WAN Jinxin.Method based on heuristic information for solving non-convex model of radiotherapy planning.Computer Engineering and Applications,2013,49（11）：265-270.

國家自然科學基金（No.61102117）。

張棟冰（1974—），男，講師，研究領域為圖像處理、算法優化；蘭義華，男，博士，講師，研究領域：圖像處理、算法優化；萬金鑫，男，研究員級高級工程師，研究方向：醫學圖像診斷方法。E-mail：291508366@qq.com

2013-02-04

2013-04-15

1002-8331（2013）11-0265-06