説明

ビーム群により送達される放射線量を最適化する方法および装置

放射線量は、正規化された形態のビーム群のモデルおよび標的線量を提供することにより、最適化される。グラム行列がモデルから決定される。標的線量がサブサンプリングされて、ビーム群の初期強度値が決定される。次に、以下のステップが、収束するまで繰り返される。非常に小さな正の値0<ε<<1を各強度値に加算して、強度値がゼロよりも大きいことを保証する。各強度値にグラム行列を乗算して、積を求め、積を要素毎に除算して正規化された標的線量にし、対応する比率を求める。比率が、数値許容誤差内ですべて1に近い場合、ビーム群の初期値を出力する。そうでない場合、次の反復前に、強度値に比率を乗算される。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、包括的には放射線治療に関し、より具体的には、粒子線群により組織に送達される放射線量の最適化に関する。
【背景技術】
【0002】
放射線量最適化の最も単純なものでは、線量行列
【数1】

は、放射線標的線量がn個のボクセルでm個のビームからなる群によって送達されることを示す。ここで、ボクセルは、組織体積の小さな立方体(mm)を表す。非負のビーム荷重のベクトル
【数2】

は、ビームを線形結合して、累積ボクセル放射線量のベクトル
【数3】

を得る。通常、線量行列Aは、比較的高く(n>>m)、疎である。問題は、標的線量パターン
【数4】

を所与として、線量不足がないことおよび最小の線量過剰を保証するビーム荷重を決めることである。これは、線形問題(LP)である。
【数5】

【0003】
小さな問題(mが〜10個のビームおよびnが〜2m個のボクセル)は、現在、2GHzのPentium(登録商標)4プロセッサで、おおよそ1秒で解くことができる。しかし、LP解法時間は、一般に、問題のサイズの三乗として拡大するため、その方法は、現在、はるかに大きなサイズを有する実際の臨床問題に対しては、実用的ではない。
【0004】
LPおよび二次錐(SOCP)法は、線量送達誤差を最小化する処方に役立つため、魅力的である。しかし、問題は、僅か数千の変数および制約に制限され、解くために数時間を要する。そのため、問題は、目的関数のLノルムの方を選ぶ、すなわち、線量不足を許容せずに、二乗和誤差を最小化する。これは、不等式問題(LSI)を有する最小二乗である。
【数6】

【0005】
LSIは、LPに対する上限の最適化として是認することができる。二乗和線量過剰を最小化することにより、僅かにより平滑である解が与えられるが、線量過剰がより大きいため、それほど安全ではない。このLSIは、行列Aが全行ランクを有し、QPが厳密に凸であると仮定して、二次計画法(QP)の特に上手く機能する事例である。ここで、LSIおよびQPは、同義で用いられる。
【0006】
高速近似法
【0007】
従来技術では、線量不足制約を回避し、単純に二乗和誤差
【数7】

または加重二乗和誤差
【数8】

を低減しようとする方法を用いることが一般的である。
【0008】
問題のサイズに起因して、解は、勾配、共役勾配、もしくは準ニュートン法、または好ましくは、メモリ制限付きブロイデン−フレッチャー−ゴルトファルプ−シャノ(L−BFGS)法を介して反復的に計算される。物理的に不可能な解を回避するために、これらの方法は、通常、負であるあらゆるビーム荷重xがゼロにリセットされるように変更される。それは、最適を見つけること、または線量不足を回避することを保証しないが、臨床での使用に十分に上手く機能するように思われる。
【0009】
勾配射影等のより原理的な方法は、非負性および最適性を保証することができるが、処理時間を増大させるため、より低速である。このため、小さな問題を数秒で解くグラフィックプロセッサユニット(GPU)での並列解が用いられることがある。
【発明の概要】
【0010】
本発明の実施の形態は、粒子線群が送達される放射線量を最適化する方法を提供する。この方法は、非負最小二乗(NNLS)問題を解くための非常に高速な乗法更新族および二次計画法(QP)の特殊なクラスを用いる。問題の変換を通じて、反復更新によって、最小距離問題(LDP)、最小二乗不等式問題(LSI)、および一般二次計画法(QP)を解くことができる。
【0011】
動機を与える用途は、制約行列を実際のメモリに記憶できないほど大きな線形問題(LP)である、放射線量最適化である。
【0012】
本方法は、LPのL1目的に、L2ノルムを用いて上限を設け、QPを取得し、乗法反復を介して様々な厳密解および近似解を展開する。
【図面の簡単な説明】
【0013】
【図1】本発明の実施形態による方法の流れ図である。
【発明を実施するための形態】
【0014】
本発明の実施の形態は、粒子線治療ビームの二乗誤差最適性および非負性荷重を保証するかなり高速の方法を提供する。
【0015】
「ペンシルビーム」治療の場合、約n〜10個のボクセルがあり、m〜10個のビームがある。その結果、計算時間は、大方、ボクセル数によって決まり、ボクセルに基づく計算を最小化するように、問題を再公式化することが望ましい。方法のステップは、当該技術分野において既知のメモリおよび入出力インタフェースに接続されたプロセッサで実行することができる。また、ステップに対応する手段もプロセッサ内で実行することができる。したがって、本発明による装置は、手段で構成される。
【0016】
図1に示されるように、方法への入力は、ビームモデル101、標的線量102、および危険に曝されている臓器(OAR)103の位置を含む。
【0017】
第1のステップ110は、正規化モデルおよび非正規化標的線量を決定する。
【0018】
ビームモデルのグラム行列
【数9】

が事前に決定される(115)。ただし、は、転置演算子である。次に、同じ変換を標的線量に適用して、
【数10】

を取得する。内積空間内のベクトル群のグラム行列は、内積のエルミート行列である。
【0019】
標的線量がサブサンプリングされて(120)、ビーム群に対する初期強度値が決定される。
【0020】
ここでは、不等式制約は、無視される。目的関数は、制約なしの最小二乗(LS)問題
【数11】

として書き換えられ、これは、元のLSIよりも2桁小さい。通常の形態では、数値精度がいくらか失われているが、これは、線量行列Aの疎性により軽減されることに留意されたい。
【0021】
制約のない最小二乗問題は、ハード不等式制約Ax≒bをソフト二乗誤差制約Ax≧bで置換する。この問題は、非負であり(A,b≧0)、過制約である(n>>m)ため、LS最適解は、2つの欠点を有する:いくつかのLS最適ビーム荷重は、かなり負になり(X<0)、いくつかのボクセルは、線量不足になる(Ax<b)。
【0022】
本発明による方法は、負のビーム荷重という第1の問題を完全になくす。非負最小二乗(NNLS)
【数12】

問題が与えられる場合、
【数13】

とする。
【0023】
乗法更新130
【数14】

は、任意の正のhおよび正の確定したQについて、任意の正の荷重ベクトルx>0から最適に収束する縮小である。この式は、xの各要素を、勾配に関連する比率で乗算する。実際には、収束は、すべての乗数が、数値許容誤差内で1に等しい場合に達成される(150)。この場合、真であれば、ビーム群の強度値が出力される(160)。そうではなく、偽の場合、強度値は、次の反復前に結果を乗算される(170)。
【0024】
線量最適化の状況では、
【数15】

であり、したがって、反復を単一の行列−ベクトル積に簡単化することができ、その場合、ビーム荷重毎に乗算130および除算140を行い、粒度の細かい並列化を受け入れるようにする。反復のコストは、グラム行列Q内の非ゼロエントリの数により支配され、非ゼロの数は、消失しない重複を有するビーム対の数である。
【0025】
線量不足の第2の問題は、最小二乗解が誤差を平滑に分散させ、その結果、治療すべき組織体積の境界上のボクセルが、健康な組織に対応する隣接ボクセルへの線量過剰を最小に抑える一手段として、線量不足になるために生じる。式(2.2)のNNLS問題の場合、線量不足の回避は、単なるソフト制約である。問題を解くことを根本的に難しくすることなく、その制約を強化するいくつかのヒューリスティックがある。
(1)Qに
【数16】

を追加することにより、治療すべきボクセルへの線量の分散に明示的にペナルティを課す。ただし、Sは、これらのボクセルを選択する二値行列であり、Iは、恒等行列であり、1は、すべて1の行列であり、tは、そのようなボクセルの数であり、τは、分散に対して線量誤差を平衡させる加重項である。
(2)境界ボクセルでの線量誤差が、その他の場所よりもはるかに重くペナルティが課されるように、二乗誤差の和を再加重し、同様に、治療体積外部のボクセルの重みを軽くする。すなわち、あるボクセル誤差加重ベクトルwおよび集合
【数17】

がある。
(3)下式であり、
【数18】

(4)ビームの減衰プロファイルに減衰がおおよそ一致する、正の値の「スカート」により囲まれた標的組織体積を変更する。
【0026】
厳密解
【0027】
線量不足の回避が、ハード制約として設定される場合、問題は、ボクセル線量の非常に広い空間内での計算を設定するとともに、ランクが不十分な行列とともに動作することが必要になり、これらの条件は、両方とも、反復ソルバの収束を遅くさせるため、はるかに厄介になる。他方、ヒューリスティックは必要ない。
【0028】
このセクションでは、線量不足がないことを保証し、乗法更新を介して解を提供する元のLSI/QPの再公式化を考慮する。
【0029】
合計線量の最小化
【0030】
線量不足禁止制約(no-underdose constraint)を受ける合計線量‖Ax‖を最小化することが望ましい場合がある。一般性を失うことなく、非負のAは、単位列和を有するため、‖Ax‖=‖x‖である。‖x‖≧‖x‖であるため、合計線量の上限を以下の最小距離問題(LDP)
【数19】

において最小化することができ、これは、等価のNNLS形態
【数20】

を有することがわかっており、等価のNNLS形態は、NNLS反復を介して解くことができる。LDPおよびNNLSの最適は、元の制約に矛盾がない場合、
【数21】

によって関係付けられ、そうでない場合、右辺はゼロである。このLDPでは、目的‖x‖+2がより平滑な解を好むため、式(1.2)の元のQPよりも大きな線量過剰が可能なことに留意されたい。また、LDPは、非常に大きく、dim(u)+dim(v)=#ボクセル+#未知のビームである。しかし、これでは、速度のために行列Aの疎性を利用することができない。
【0031】
LSIを回転させてLDS形態にする
【0032】
元のLSIは、
【数22】

となるように、変数をzに変更することにより、僅かにより複雑なLDPに変換することによって、厳密に解くことができる。
【0033】
次に、LSIは、等価のLDPとして書き換えることができる。Aの「薄い」QR分解は、UR=Aであり、ただし、Uはユニタリであり、Rは対角上で上三角かつ非負である。QPは、Lノルムを保持した状態で変数
【数23】

を変更して書き換えられて、最小
【数24】

を取得し、
【数25】

であり、†は、疑似逆を示す。
【0034】
対応するNNLSは、
【数26】

であり、最適は、
【数27】

によって関係付けられ、このとき
【数28】

である。
【0035】
対応する乗法更新は、線量不足なし、負の荷重なし、かつ最小二乗和線量過剰で線量問題を解く。しかし、前のLDP/NNLSが疎であったのに対して、このNNLSは、密であり、反復毎にかなり多くの更なる計算が必要とされる。また、問題が2回変換されることによる精度の損失もある。
【0036】
QP双対を解く
【0037】
厳密な凸二次計画が主形態
【数29】

で与えられる場合、二次計画のラグランジュ双対は、
【数30】

の形態をとり、
【数31】

である。
【0038】
主最適および双対最適は、
【数32】

によって関係付けられ、双対は、NNLS反復により解くことができる。x≧0が、制約Ax≧bに明示的に符号化されない限り保証されないことに留意されたい。これを式(1.2)の線量最適化LSI/QPに対して特殊化するために、
【数33】

および
【数34】

を設定し、
【数35】

を解く。
ここで、制約x≧0を組み込むために
【数36】

および
【数37】

である。逆数(AA)−1は、上述した畳み込み分解を用いて効率的に計算することができる。(AA)の条件付けが不十分である可能性が高いため、注意しなければならない。
【0039】
反復を介して縮小した主QPを解く
【0040】
放射線量の誤差は、矩形の部分領域内の線量値が欠けることに関して、比率の分母に、追加のペナルティを表す補助行列−ベクトル積を追加することにより、標的線量の部分領域に縮小される。
【0041】
特に、元のQPの制約群を、k≦2mボクセルでのみ、例えば、境界上にあり、ビーム中心から最も離れたボクセルでのみ、線量不足を禁止するように制限することを考える。非負のQ、h、A、bを有する制限付きQPおよびすべてのi、Qij>0、h>0について、補助行列−ベクトル積
【数38】

が提供され、ただし、b’’は、k個の制約付きボクセルに対応するbの部分集合であり、A’’は、線量行列Aの行の対応する部分集合である。
【0042】
補助行列−ベクトル積
【数39】

は、潜在価格(shadow price)としても知られる双対変数を制約し、この変数は、制約が有効な場合には解において正であり、その他の場合はゼロである。本明細書でのような制約付き最適化では、潜在価格は、1単位だけ制約を緩和することにより得られる最適問題の最適解の目的値の変更を表す。
【0043】
この方法は、ハード制約を任意のボクセル群に対する最大累積線量または最小累積線量に対して課すためにも用いることができる。λの更新は、λに関するラグランジアンの勾配に沿って小さなステップに変更することができる。
【0044】
反復の計算詳細
【0045】
乗法反復手順は、線形複雑性を有し、正円錐内の任意の箇所から線形率で収束する。一定数の精度ビットが、各反復で得られる。乗法更新は、すべての非負制約が暗黙的に強制されるという利点を有する。しかし、数値誤差が、非常に小さな強度値をゼロに丸め、その後、強度がゼロ境界に「固定」されたままである場合もある。
【0046】
したがって、ステップ130において、非常に小さな正の値0<ε<<1が各強度値に加算されて、乗算前に正の強度値を生成する。これにより、反復は、ゼロ境界を免れる。変数毎の更新は、メッセージパッシングプロシージャとしてベクトルに並行して実行することができるか、または逐次実行することができ、この場合、ビーム荷重Xの更新は、Qij≠0である場合には、常に、荷重Xの更新から利益を受ける。この結果、収束がいくらか高速になる。
【0047】
この用途では、±2.5%の誤差の線量を施すことが珍しくないため、早期終了が推奨され、したがって、7ビットを越えるビーム荷重xのいかなる精度も無駄になる。
【0048】
線量最適化の計算詳細
【0049】
可分畳み込み(separable convolution)
【0050】
行列Q=AAの形成は、通常、これらの方式の中で最も時間のかかる演算である。多くの場合、その全体を回避することができる。鍵となる洞察は、Qijが2つのビームの内積であり、したがって、畳み込みにより計算することができることである。これを可能な限り単純なビームモデルにおいて利用することができる。ビームは、標準偏差σを有するxy平面ガウス関数と、zでのブラッグ曲線との外積である。結果として、p、q、rだけずれた2つのビームの3D畳み込みは、1D畳み込みの外積に分解される。
【数40】

【0051】
上記式において、
【数41】

は単純に
【数42】

の逆数である。したがって、行列Qのエントリは、自動畳み込みブラッグ曲線と2つの1Dガウス関数との外積である。
【数43】

【0052】
この構造を利用して、恒等式の再形成を介してxに対するQの動作を高速で求めることができる。
【数44】

【0053】
この恒等式を外積毎に1回用いることにより、行の係数(Bだけ、すなわち、平面内のビームスポットの数だけ、スカラー乗算の総数が低減される。利点として、x+y+zボクセルの体積の場合、ベクトルQ、Q、Q内のx+y+z個のみのエントリが、グラム行列Q内の(xyz)個のエントリに代えて用いられる。
【0054】
同じ論理により、1Dガウス関数を用いてx軸およびy軸に沿って3D標的線量bを逐次畳み込み、次に、ブラッグ曲線を用いてz軸に沿って逐次畳み込むことにより、ベクトルh=Abを効率的に形成することができる。
【0055】
有界領域内の加重誤差
【0056】
ガウス畳み込みが、有限xy平面全体に広がり、これが、1があらゆる箇所での最小二乗誤差であることを暗示することに留意する。標的線量は、略あらゆる箇所で暗黙的にゼロであるため、有限の広がりは、標的の境界での減衰に寄与する。この寄与は、無視できるほどに小さく見えるが、場合によっては、その積分は、ある小さな体積を有界とする。
【数45】

ただし、Φ()は、誤差関数である。ここでも、これは非負の外積である。
【0057】
加重誤差
【0058】
放射線量の誤差は、矩形の部分領域内の線量値が欠けることに関して、比率の分母に、追加のペナルティを表す仕様にない補助的な加重誤差セクション行列−ベクトル積の第2の式を追加することにより、部分領域に縮小される。
【0059】
グラム行列Qの有界形態は、加重二乗誤差
【数46】

を最小化する場合に有用になる。通常の形態と共に動作するには、加重積は、
【数47】

である。
【0060】
不都合なことに、wの大半の選択が、構造体Qへの外積分解の使用を回避する。しかし、いくつかの有用な場合、加重されたQを2つの行列に分けることにより、外積を得ることができる。
【数48】

【0061】
例えば、w=c>1により、ある境界ボックス、例えば、危険に曝されている臓器の周囲の境界ボックスを均一に過度に加重することにより、第2の被加数は、単純に、上述したように、C−1を乗算される有界積分Q行列である。小数の個々のピクセルを過度に加重するために、
【数49】

とする。
【0062】
そして、W’がdiag(w’)の非ゼロ列を含むものとすると、第2の被加数は、
【数50】

として書くことができる。
【0063】
【数51】

の列は、hと同じ方法で事前に計算することができ、xに対するその動作は、
【数52】

として最も効率的に計算される。
【0064】
外積乗算、xに対する2つの被加数の動作、および和を用いるために、結果が別個に求められることに留意されたい。
【0065】
[発明の効果]
本発明の実施の形態は、二次計画法(QP)を用いて、放射線治療の粒子線を最適化する極めて大きく複雑な問題を解く。
【0066】
QPは、ネットワーク最適化、機械学習、物流、制御、および他の多くの用途に広く普及している。QPは、負になり得ない、多くの場合にジュール、時間、または金銭で示される数量を推定する場合に特に有用である。
【0067】
QPを解く方法は、第二次世界大戦から集中的に開発され、少なくとも1つのノーベル賞がその分野の研究に結びついている(Harry Markowitz−1990)。近年、強力な方法の大きなコレクションがあるが、問題のサイズは、手に負えなくなっている。
【0068】
本発明は、線形率収束を用いて、線形時間、線形空間内でQPを解く新しく非常に簡易な不動点を提供する。これは、QPが>10個の変数および>10個の制約を有することができる粒子線の最適化に対して特に有用であることがわかる。
【0069】
不動点は、その簡易性に起因して、非常に高速であり、並列化可能である。従来技術によるQPでは、問題を解くために数時間を要した。本明細書に記載の方法は、ビデオフレームレートで問題を解くことができる。
【0070】
本発明について、好ましい実施形態の例として説明したが、他の様々な適合および変更を、本発明の趣旨および範囲内で行うことができることを理解されたい。したがって、添付の特許請求の目的は、本発明の真の趣旨および範囲内にあるすべての変形および変更を包含することである。

【特許請求の範囲】
【請求項1】
ビーム群により送達される放射線量を最適化する方法であって、各ビームは、粒子線治療ビームであり、該方法は、
前記ビーム群のモデルおよび標的線量を提供するステップと、
前記モデルから、グラム行列を含む、通常形態のモデルおよび標的線量を決定するステップと、
前記標的線量をサブサンプリングして、前記ビーム群の初期強度値を決定し、収束するまで、以下のステップ:
非常に小さな正の値0<ε<<1を各強度値に加算して前記強度値がゼロよりも大きいことを保証する、加算するステップ、
強度値のベクトルに前記グラム行列を乗算して積を求める、乗算するステップ、
要素毎に、前記積を除算して標的線量形態にし、対応する比率を求める、除算するステップ、および
前記比率が、数値許容誤差内ですべて1に近いか否かを判断し、真の場合、前記ビーム群の前記強度値を出力し、そうではなく、偽の場合、次の繰り返しの前に、前記強度値に前記比率を乗算する、判断するステップ、
を繰り返す、サブサンプリングするステップと、
を含む、ビーム群により送達される放射線量を最適化する方法。
【請求項2】
前記グラム行列の各要素は、ビーム対毎に、該ビーム対により送達される前記放射線量の内積を表す、請求項1に記載の方法。
【請求項3】
可能な限り、前記グラム行列を分解して、外積にし、該外積の係数を前記ビーム群の前記強度値に直接適用する、分解するステップをさらに含む、請求項1に記載の方法。
【請求項4】
矩形の部分領域内の線量値が欠けることに関して、前記比率の分母に、追加のペナルティを表す補助行列−ベクトル積を追加することにより、前記放射線量の誤差を前記標的線量の前記部分領域に縮小するステップをさらに含む、請求項1に記載の方法。
【請求項5】
前記比率の分子に、ハード放射線制約を確保する潜在価格を表す補助行列−ベクトル積を追加することにより、ハード放射線制約を前記放射線量の任意の部分領域に対して課すステップをさらに含む、請求項1に記載の方法。
【請求項6】
前記標的線量に対して規定される制約毎の潜在価格を計算するために、前記モデルを双対形態にし、該双対モデルにサブサンプリングステップを適用するステップをさらに含む、請求項1に記載の方法。
【請求項7】
ビーム群により送達される放射線量を最適化する装置であって、各ビームは粒子線治療ビームであり、該装置は、
前記ビーム群のモデルおよび標的線量を提供する手段と、
前記モデルから、グラム行列を含む、通常形態のモデルおよび標的線量を決定する手段と、
前記標的線量をサブサンプリングする手段であって、前記ビーム群の初期強度値を決定し、収束するまで、以下のステップ:
非常に小さな正の値0<ε<<1を各強度値に加算し、前記強度値がゼロよりも大きいことを保証する、加算するステップ、
強度値のベクトルに前記グラム行列を乗算し、積を求める、乗算するステップ、
要素毎に、前記積を除算し、標的線量形態にし、対応する比率を求める、除算するステップ、および
前記比率が、数値許容誤差内ですべて1に近いか否かを判断し、真の場合、前記ビーム群の前記強度値を出力し、そうではなく、偽の場合、次の繰り返しの前に、前記強度値に前記比率を乗算する、判断するステップ、
を繰り返す、サブサンプリングする手段と、
を備える、ビーム群により送達される放射線量を最適化する装置。
【請求項8】
前記グラム行列の各要素は、ビーム対毎に、該ビーム対により送達される前記放射線量の内積を表す、請求項7に記載の装置。
【請求項9】
可能な限り、前記グラム行列を分解して、外積にし、該外積の係数を前記ビーム群の前記強度値に直接適用する、分解する手段をさらに備える、請求項7に記載の装置。
【請求項10】
矩形の部分領域内の線量値が欠けることに関して、前記比率の分母に、追加のペナルティを表す補助行列−ベクトル積を追加することにより、前記放射線量の誤差を前記標的線量の前記部分領域に縮小する手段をさらに備える、請求項7に記載の装置。
【請求項11】
前記比率の分子に、ハード放射線制約を確保する潜在価格を表す補助行列−ベクトル積を追加することにより、ハード放射線制約を前記放射線量の任意の部分領域に対して課す手段をさらに備える、請求項7に記載の装置。
【請求項12】
前記モデルを双対形態にし、該双対モデルにサブサンプリングステップを適用して、前記標的線量に対して規定される制約毎の潜在価格を計算する、双対形態にし、適用する手段をさらに備える、請求項7に記載の装置。

【図1】
image rotate


【公表番号】特表2013−500048(P2013−500048A)
【公表日】平成25年1月7日(2013.1.7)
【国際特許分類】
【出願番号】特願2012−504586(P2012−504586)
【出願日】平成23年5月9日(2011.5.9)
【特許番号】特許第5063828号(P5063828)
【特許公報発行日】平成24年10月31日(2012.10.31)
【国際出願番号】PCT/JP2011/061113
【国際公開番号】WO2011/148800
【国際公開日】平成23年12月1日(2011.12.1)
【出願人】(597067574)ミツビシ・エレクトリック・リサーチ・ラボラトリーズ・インコーポレイテッド (484)
【住所又は居所原語表記】201 BROADWAY, CAMBRIDGE, MASSACHUSETTS 02139, U.S.A.
【Fターム(参考)】