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

Abstract

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

Description

本発明は、包括的には放射線治療に関し、より具体的には、粒子線群により組織に送達される放射線量の最適化に関する。

放射線量最適化の最も単純なものでは、線量行列

は、放射線標的線量がｎ個のボクセルでｍ個のビームからなる群によって送達されることを示す。ここで、ボクセルは、組織体積の小さな立方体（ｍｍ^３）を表す。非負のビーム荷重のベクトル

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

を得る。通常、線量行列Ａは、比較的高く（ｎ＞＞ｍ）、疎である。問題は、標的線量パターン

を所与として、線量不足がないことおよび最小の線量過剰を保証するビーム荷重を決めることである。これは、線形問題（ＬＰ）である。

小さな問題（ｍが〜１０^３個のビームおよびｎが〜２ｍ個のボクセル）は、現在、２ＧＨｚのＰｅｎｔｉｕｍ(登録商標)４プロセッサで、おおよそ１秒で解くことができる。しかし、ＬＰ解法時間は、一般に、問題のサイズの三乗として拡大するため、その方法は、現在、はるかに大きなサイズを有する実際の臨床問題に対しては、実用的ではない。

ＬＰおよび二次錐（ＳＯＣＰ）法は、線量送達誤差を最小化する処方に役立つため、魅力的である。しかし、問題は、僅か数千の変数および制約に制限され、解くために数時間を要する。そのため、問題は、目的関数のＬ_２ノルムの方を選ぶ、すなわち、線量不足を許容せずに、二乗和誤差を最小化する。これは、不等式問題（ＬＳＩ）を有する最小二乗である。

ＬＳＩは、ＬＰに対する上限の最適化として是認することができる。二乗和線量過剰を最小化することにより、僅かにより平滑である解が与えられるが、線量過剰がより大きいため、それほど安全ではない。このＬＳＩは、行列Ａが全行ランクを有し、ＱＰが厳密に凸であると仮定して、二次計画法（ＱＰ）の特に上手く機能する事例である。ここで、ＬＳＩおよびＱＰは、同義で用いられる。

高速近似法

従来技術では、線量不足制約を回避し、単純に二乗和誤差

または加重二乗和誤差

を低減しようとする方法を用いることが一般的である。

問題のサイズに起因して、解は、勾配、共役勾配、もしくは準ニュートン法、または好ましくは、メモリ制限付きブロイデン−フレッチャー−ゴルトファルプ−シャノ（Ｌ−ＢＦＧＳ）法を介して反復的に計算される。物理的に不可能な解を回避するために、これらの方法は、通常、負であるあらゆるビーム荷重ｘ_ｉがゼロにリセットされるように変更される。それは、最適を見つけること、または線量不足を回避することを保証しないが、臨床での使用に十分に上手く機能するように思われる。

勾配射影等のより原理的な方法は、非負性および最適性を保証することができるが、処理時間を増大させるため、より低速である。このため、小さな問題を数秒で解くグラフィックプロセッサユニット（ＧＰＵ）での並列解が用いられることがある。

本発明の実施の形態は、粒子線群が送達される放射線量を最適化する方法を提供する。この方法は、非負最小二乗（ＮＮＬＳ）問題を解くための非常に高速な乗法更新族および二次計画法（ＱＰ）の特殊なクラスを用いる。問題の変換を通じて、反復更新によって、最小距離問題（ＬＤＰ）、最小二乗不等式問題（ＬＳＩ）、および一般二次計画法（ＱＰ）を解くことができる。

動機を与える用途は、制約行列を実際のメモリに記憶できないほど大きな線形問題（ＬＰ）である、放射線量最適化である。

本方法は、ＬＰのＬ１目的に、Ｌ２ノルムを用いて上限を設け、ＱＰを取得し、乗法反復を介して様々な厳密解および近似解を展開する。

本発明の実施形態による方法の流れ図である。

本発明の実施の形態は、粒子線治療ビームの二乗誤差最適性および非負性荷重を保証するかなり高速の方法を提供する。

「ペンシルビーム」治療の場合、約ｎ〜１０^６個のボクセルがあり、ｍ〜１０^４個のビームがある。その結果、計算時間は、大方、ボクセル数によって決まり、ボクセルに基づく計算を最小化するように、問題を再公式化することが望ましい。方法のステップは、当該技術分野において既知のメモリおよび入出力インタフェースに接続されたプロセッサで実行することができる。また、ステップに対応する手段もプロセッサ内で実行することができる。したがって、本発明による装置は、手段で構成される。

図１に示されるように、方法への入力は、ビームモデル１０１、標的線量１０２、および危険に曝されている臓器（ＯＡＲ）１０３の位置を含む。

第１のステップ１１０は、正規化モデルおよび非正規化標的線量を決定する。

ビームモデルのグラム行列

が事前に決定される（１１５）。ただし、^Ｔは、転置演算子である。次に、同じ変換を標的線量に適用して、

を取得する。内積空間内のベクトル群のグラム行列は、内積のエルミート行列である。

標的線量がサブサンプリングされて（１２０）、ビーム群に対する初期強度値が決定される。

ここでは、不等式制約は、無視される。目的関数は、制約なしの最小二乗（ＬＳ）問題

として書き換えられ、これは、元のＬＳＩよりも２桁小さい。通常の形態では、数値精度がいくらか失われているが、これは、線量行列Ａの疎性により軽減されることに留意されたい。

制約のない最小二乗問題は、ハード不等式制約Ａｘ≒ｂをソフト二乗誤差制約Ａｘ≧ｂで置換する。この問題は、非負であり（Ａ，ｂ≧０）、過制約である（ｎ＞＞ｍ）ため、ＬＳ最適解は、２つの欠点を有する：いくつかのＬＳ最適ビーム荷重は、かなり負になり（Ｘ_ｉ＜０）、いくつかのボクセルは、線量不足になる（Ａ_ｉｘ＜ｂ_ｉ）。

本発明による方法は、負のビーム荷重という第１の問題を完全になくす。非負最小二乗（ＮＮＬＳ）

問題が与えられる場合、

とする。

乗法更新１３０

は、任意の正のｈおよび正の確定したＱについて、任意の正の荷重ベクトルｘ＞０から最適に収束する縮小である。この式は、ｘの各要素を、勾配に関連する比率で乗算する。実際には、収束は、すべての乗数が、数値許容誤差内で１に等しい場合に達成される（１５０）。この場合、真であれば、ビーム群の強度値が出力される（１６０）。そうではなく、偽の場合、強度値は、次の反復前に結果を乗算される（１７０）。

線量最適化の状況では、

であり、したがって、反復を単一の行列−ベクトル積に簡単化することができ、その場合、ビーム荷重毎に乗算１３０および除算１４０を行い、粒度の細かい並列化を受け入れるようにする。反復のコストは、グラム行列Ｑ内の非ゼロエントリの数により支配され、非ゼロの数は、消失しない重複を有するビーム対の数である。

線量不足の第２の問題は、最小二乗解が誤差を平滑に分散させ、その結果、治療すべき組織体積の境界上のボクセルが、健康な組織に対応する隣接ボクセルへの線量過剰を最小に抑える一手段として、線量不足になるために生じる。式（２．２）のＮＮＬＳ問題の場合、線量不足の回避は、単なるソフト制約である。問題を解くことを根本的に難しくすることなく、その制約を強化するいくつかのヒューリスティックがある。
（１）Ｑに

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

がある。
（３）下式であり、

（４）ビームの減衰プロファイルに減衰がおおよそ一致する、正の値の「スカート」により囲まれた標的組織体積を変更する。

厳密解

線量不足の回避が、ハード制約として設定される場合、問題は、ボクセル線量の非常に広い空間内での計算を設定するとともに、ランクが不十分な行列とともに動作することが必要になり、これらの条件は、両方とも、反復ソルバの収束を遅くさせるため、はるかに厄介になる。他方、ヒューリスティックは必要ない。

このセクションでは、線量不足がないことを保証し、乗法更新を介して解を提供する元のＬＳＩ／ＱＰの再公式化を考慮する。

合計線量の最小化

線量不足禁止制約（no-underdose constraint）を受ける合計線量‖Ａｘ‖_１を最小化することが望ましい場合がある。一般性を失うことなく、非負のＡは、単位列和を有するため、‖Ａｘ‖_１＝‖ｘ‖_１である。‖ｘ‖_２≧‖ｘ‖_１であるため、合計線量の上限を以下の最小距離問題（ＬＤＰ）

において最小化することができ、これは、等価のＮＮＬＳ形態

を有することがわかっており、等価のＮＮＬＳ形態は、ＮＮＬＳ反復を介して解くことができる。ＬＤＰおよびＮＮＬＳの最適は、元の制約に矛盾がない場合、

によって関係付けられ、そうでない場合、右辺はゼロである。このＬＤＰでは、目的‖ｘ‖＋２がより平滑な解を好むため、式（１．２）の元のＱＰよりも大きな線量過剰が可能なことに留意されたい。また、ＬＤＰは、非常に大きく、ｄｉｍ（ｕ）＋ｄｉｍ（ｖ）＝＃ボクセル＋＃未知のビームである。しかし、これでは、速度のために行列Ａの疎性を利用することができない。

ＬＳＩを回転させてＬＤＳ形態にする

元のＬＳＩは、

となるように、変数をｚに変更することにより、僅かにより複雑なＬＤＰに変換することによって、厳密に解くことができる。

次に、ＬＳＩは、等価のＬＤＰとして書き換えることができる。Ａの「薄い」ＱＲ分解は、ＵＲ＝Ａであり、ただし、Ｕはユニタリであり、Ｒは対角上で上三角かつ非負である。ＱＰは、Ｌ_２ノルムを保持した状態で変数

を変更して書き換えられて、最小

を取得し、

であり、†は、疑似逆を示す。

対応するＮＮＬＳは、

であり、最適は、

によって関係付けられ、このとき

である。

対応する乗法更新は、線量不足なし、負の荷重なし、かつ最小二乗和線量過剰で線量問題を解く。しかし、前のＬＤＰ／ＮＮＬＳが疎であったのに対して、このＮＮＬＳは、密であり、反復毎にかなり多くの更なる計算が必要とされる。また、問題が２回変換されることによる精度の損失もある。

ＱＰ双対を解く

厳密な凸二次計画が主形態

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

の形態をとり、

である。

主最適および双対最適は、

によって関係付けられ、双対は、ＮＮＬＳ反復により解くことができる。ｘ≧０が、制約Ａｘ≧ｂに明示的に符号化されない限り保証されないことに留意されたい。これを式（１．２）の線量最適化ＬＳＩ／ＱＰに対して特殊化するために、

および

を設定し、

を解く。
ここで、制約ｘ≧０を組み込むために

および

である。逆数（Ａ^ＴＡ）^−１は、上述した畳み込み分解を用いて効率的に計算することができる。（Ａ^ＴＡ）の条件付けが不十分である可能性が高いため、注意しなければならない。

反復を介して縮小した主ＱＰを解く

放射線量の誤差は、矩形の部分領域内の線量値が欠けることに関して、比率の分母に、追加のペナルティを表す補助行列−ベクトル積を追加することにより、標的線量の部分領域に縮小される。

特に、元のＱＰの制約群を、ｋ≦２ｍボクセルでのみ、例えば、境界上にあり、ビーム中心から最も離れたボクセルでのみ、線量不足を禁止するように制限することを考える。非負のＱ、ｈ、Ａ、ｂを有する制限付きＱＰおよびすべてのｉ、Ｑ_ｉｊ＞０、ｈ_ｉ＞０について、補助行列−ベクトル積

が提供され、ただし、ｂ’’は、ｋ個の制約付きボクセルに対応するｂの部分集合であり、Ａ’’は、線量行列Ａの行の対応する部分集合である。

補助行列−ベクトル積

は、潜在価格（shadow price）としても知られる双対変数を制約し、この変数は、制約が有効な場合には解において正であり、その他の場合はゼロである。本明細書でのような制約付き最適化では、潜在価格は、１単位だけ制約を緩和することにより得られる最適問題の最適解の目的値の変更を表す。

この方法は、ハード制約を任意のボクセル群に対する最大累積線量または最小累積線量に対して課すためにも用いることができる。λの更新は、λに関するラグランジアンの勾配に沿って小さなステップに変更することができる。

反復の計算詳細

乗法反復手順は、線形複雑性を有し、正円錐内の任意の箇所から線形率で収束する。一定数の精度ビットが、各反復で得られる。乗法更新は、すべての非負制約が暗黙的に強制されるという利点を有する。しかし、数値誤差が、非常に小さな強度値をゼロに丸め、その後、強度がゼロ境界に「固定」されたままである場合もある。

したがって、ステップ１３０において、非常に小さな正の値０＜ε＜＜１が各強度値に加算されて、乗算前に正の強度値を生成する。これにより、反復は、ゼロ境界を免れる。変数毎の更新は、メッセージパッシングプロシージャとしてベクトルに並行して実行することができるか、または逐次実行することができ、この場合、ビーム荷重Ｘ_ｉの更新は、Ｑ_ｉｊ≠０である場合には、常に、荷重Ｘ_ｊの更新から利益を受ける。この結果、収束がいくらか高速になる。

この用途では、±２．５％の誤差の線量を施すことが珍しくないため、早期終了が推奨され、したがって、７ビットを越えるビーム荷重ｘのいかなる精度も無駄になる。

線量最適化の計算詳細

可分畳み込み（separable convolution）

行列Ｑ＝Ａ^ＴＡの形成は、通常、これらの方式の中で最も時間のかかる演算である。多くの場合、その全体を回避することができる。鍵となる洞察は、Ｑ_ｉｊが２つのビームの内積であり、したがって、畳み込みにより計算することができることである。これを可能な限り単純なビームモデルにおいて利用することができる。ビームは、標準偏差σを有するｘｙ平面ガウス関数と、ｚでのブラッグ曲線との外積である。結果として、ｐ、ｑ、ｒだけずれた２つのビームの３Ｄ畳み込みは、１Ｄ畳み込みの外積に分解される。

上記式において、

は単純に

の逆数である。したがって、行列Ｑのエントリは、自動畳み込みブラッグ曲線と２つの１Ｄガウス関数との外積である。

この構造を利用して、恒等式の再形成を介してｘに対するＱの動作を高速で求めることができる。

この恒等式を外積毎に１回用いることにより、行の係数（Ｂ_Ｇ）^２だけ、すなわち、平面内のビームスポットの数だけ、スカラー乗算の総数が低減される。利点として、ｘ＋ｙ＋ｚボクセルの体積の場合、ベクトルＱ_Ｂ、Ｑ_Ｇ、Ｑ_Ｇ内のｘ＋ｙ＋ｚ個のみのエントリが、グラム行列Ｑ内の（ｘｙｚ）^２個のエントリに代えて用いられる。

同じ論理により、１Ｄガウス関数を用いてｘ軸およびｙ軸に沿って３Ｄ標的線量ｂを逐次畳み込み、次に、ブラッグ曲線を用いてｚ軸に沿って逐次畳み込むことにより、ベクトルｈ＝Ａ^ｔｂを効率的に形成することができる。

有界領域内の加重誤差

ガウス畳み込みが、有限ｘｙ平面全体に広がり、これが、１があらゆる箇所での最小二乗誤差であることを暗示することに留意する。標的線量は、略あらゆる箇所で暗黙的にゼロであるため、有限の広がりは、標的の境界での減衰に寄与する。この寄与は、無視できるほどに小さく見えるが、場合によっては、その積分は、ある小さな体積を有界とする。

ただし、Φ（）は、誤差関数である。ここでも、これは非負の外積である。

加重誤差

放射線量の誤差は、矩形の部分領域内の線量値が欠けることに関して、比率の分母に、追加のペナルティを表す仕様にない補助的な加重誤差セクション行列−ベクトル積の第２の式を追加することにより、部分領域に縮小される。

グラム行列Ｑの有界形態は、加重二乗誤差

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

である。

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

例えば、ｗ_ｉ＝ｃ＞１により、ある境界ボックス、例えば、危険に曝されている臓器の周囲の境界ボックスを均一に過度に加重することにより、第２の被加数は、単純に、上述したように、Ｃ^２−１を乗算される有界積分Ｑ行列である。小数の個々のピクセルを過度に加重するために、

とする。

そして、Ｗ’がｄｉａｇ（ｗ’）の非ゼロ列を含むものとすると、第２の被加数は、

として書くことができる。

の列は、ｈと同じ方法で事前に計算することができ、ｘに対するその動作は、

として最も効率的に計算される。

外積乗算、ｘに対する２つの被加数の動作、および和を用いるために、結果が別個に求められることに留意されたい。

［発明の効果］
本発明の実施の形態は、二次計画法（ＱＰ）を用いて、放射線治療の粒子線を最適化する極めて大きく複雑な問題を解く。

ＱＰは、ネットワーク最適化、機械学習、物流、制御、および他の多くの用途に広く普及している。ＱＰは、負になり得ない、多くの場合にジュール、時間、または金銭で示される数量を推定する場合に特に有用である。

ＱＰを解く方法は、第二次世界大戦から集中的に開発され、少なくとも１つのノーベル賞がその分野の研究に結びついている（Ｈａｒｒｙ Ｍａｒｋｏｗｉｔｚ−１９９０）。近年、強力な方法の大きなコレクションがあるが、問題のサイズは、手に負えなくなっている。

本発明は、線形率収束を用いて、線形時間、線形空間内でＱＰを解く新しく非常に簡易な不動点を提供する。これは、ＱＰが＞１０^４個の変数および＞１０^６個の制約を有することができる粒子線の最適化に対して特に有用であることがわかる。

不動点は、その簡易性に起因して、非常に高速であり、並列化可能である。従来技術によるＱＰでは、問題を解くために数時間を要した。本明細書に記載の方法は、ビデオフレームレートで問題を解くことができる。

本発明について、好ましい実施形態の例として説明したが、他の様々な適合および変更を、本発明の趣旨および範囲内で行うことができることを理解されたい。したがって、添付の特許請求の目的は、本発明の真の趣旨および範囲内にあるすべての変形および変更を包含することである。

Claims (12)

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

Images