(57)【特許請求の範囲】
【請求項１】
媒質を微小な厚さΔｚｎ  （ｎ＝１，２，３，．．．）を有する多数の層に分割し、入射するガウスビームのエネルギーＥ０  と粒子数Ｎ（Ｅ０  ）、２乗平均散乱角と発散と２乗平均半径の各初期値、および媒質の物性値分布を入力し、番号ｎを１とする第１のステップと、エネルギー減衰率が規定値となる厚さの増分Δｚｎ  とエネルギー減衰後のエネルギーＥｎ  と表面からの深さｚｎ  を計算する第２のステップと、核反応による粒子数ウェイトＷ（ｚｎ  ）を計算する第３のステップと、平均エネルギー損失を計算する第４のステップと、２乗平均散乱角を計算する第５のステップと、発散を計算する第６のステップと、２乗平均半径を計算する第７のステップと、エネルギーＥｎ  が規定値以下になったかを判定する第８のステップとを有し、前記エネルギーＥｎ  が規定値に達するまでは前記ｎを１つ増加させて第２から第７のステップを繰り返し実行し、規定値に達した後は計算したすべてのｎの前記平均エネルギー損失に対して、レンジストラグリングと入射ビームのエネルギー分散を考慮して実効エネルギー損失Ｓ（Ｅ０  ，ｚｎ  ）を計算する第９のステップと、線量分布Ｄ（ｘ，ｚｎ  ）を前記すべてのｎに対して計算する第１０のステップとを実行することを特徴とする粒子線線量分布の計算方法。
【請求項２】
粒子線として陽子線、電子線、または重粒子線を選んだことを特徴とする請求項１記載の粒子線線量分布の計算方法。
【請求項３】
発散を以下の式で逐次計算したことを特徴とする請求項１記載の粒子線線量分布の計算方法。
【数１】
0003620227.tif 000014
【請求項４】
２乗平均半径を以下の式で逐次計算したことを特徴とする請求項１記載の粒子線線量分布の計算方法。
【数２】
0003620227.tif 000015
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
この発明は、陽子線・電子線・重粒子線（炭素線、アルゴン線など）で代表される粒子線の線量分布、特に体内における線量分布を計算する方法に関するものである。
【０００２】
【従来の技術】
従来の陽子線の体内線量分布の計算方法は、例えば１９９５年発行の雑誌「フィジックス・イン・メディスン・アンド・バイオロジー」（Ｐｈｙｓｉｃｓ  ｉｎ  Ｍｅｄｉｃｉｎｅ  ａｎｄ  Ｂｉｏｌｏｇｙ）、４０巻、１０３１頁〜１０４３頁に記載されており、体内の線量分布Ｄ（ｘ，ｚ）は次式で与えられることが知られている。
【０００３】
【数３】
0003620227.tif 000002
【０００４】
ここに、入射エネルギーＥ０  、粒子数Ｎ（Ｅ０  ）の陽子線ガウスビーム（ガウスビームとは断面内の強度分布がガウス分布であるビームを指す）において、体内深さｚにおける単位深さ当たりの陽子線の実効エネルギー損失をＳ（Ｅ０  ，ｚ）、体内深さｚにおける陽子線ビームの２乗平均半径ｒ２（ｚ）ＡＶＥを
【０００５】
【数４】
0003620227.tif 000003
【０００６】
とした。また、ビーム進行方向と直交する平面において、ビーム軸からの距離をｘとした。なお、説明を簡単にするために、以下では入射ビームの粒子数Ｎ（Ｅ０  ）を１とする。
【０００７】
式（１）で線量分布を計算する場合、媒質内（例えば、体内）の実効エネルギー損失Ｓ（Ｅ０  ，ｚ）と陽子線ビームの２乗平均半径ｒ２（ｚ）ＡＶＥをあらかじめ計算しておく必要がある。前記文献では、Ｓ（Ｅ０  ，ｚ）は他の文献に記載されているデータテーブルに基づき計算し、陽子線ビームの２乗平均半径ｒ２（ｚ）ＡＶＥは以下の３つの式で順次計算することを提案している。
【０００８】
【数５】
0003620227.tif 000004
【０００９】
ここに、ｚｎ  （ｎ＝１，２，３，．．．）は体内の深さｚを離散的に与えた座標である。つまり、体内に陽子線ビームを照射したときに深さｚｎ  （ｎ＝１，２，３，．．．）において、ビームの２乗平均半径ｒ２（ｚ）ＡＶＥを逐次計算していくことになる。なお、式（２）において、
【００１０】
【数６】
0003620227.tif 000005
【００１１】
である。
【００１２】
【発明が解決しようとする課題】
従来の粒子線線量分布の計算方法では、複雑な式（２）、（３）、（４）を逐次計算する必要があり、人体内部の線量分布を計算する場合、その計算時間が増大する欠点があった。
【００１３】
この発明は、かかる問題点を解決するためになされたものであり、陽子線・電子線・重粒子線などの粒子線の線量分布、特に体内における線量分布を従来よりも高速・高精度に計算する方法を提供するものである。
【００１４】
【課題を解決するための手段】
この発明の第１の方法に係る粒子線線量分布の計算方法は、媒質を微小な厚さΔｚｎ  （ｎ＝１，２，３，．．．）を有する多数の層に分割し、入射するガウスビームのエネルギーＥ０  と粒子数Ｎ（Ｅ０  ）、２乗平均散乱角と発散と２乗平均半径の各初期値、および媒質の物性値分布を入力し、番号ｎを１とする第１のステップと、エネルギー減衰率が規定値となる厚さの増分Δｚｎ  とエネルギー減衰後のエネルギーＥｎ  と表面からの深さｚｎ  を計算する第２のステップと、核反応による粒子数ウェイトＷ（ｚｎ  ）を計算する第３のステップと、平均エネルギー損失を計算する第４のステップと、２乗平均散乱角を計算する第５のステップと、発散を計算する第６のステップと、２乗平均半径を計算する第７のステップと、エネルギーＥｎ  が規定値以下になったかを判定する第８のステップとを有し、前記エネルギーＥｎ  が規定値に達するまでは前記ｎを１つ増加させて第２から第７のステップを繰り返し実行し、規定値に達した後は計算したすべてのｎの前記平均エネルギー損失に対して、レンジストラグリングと入射ビームのエネルギー分散を考慮して実効エネルギー損失Ｓ（Ｅ０  ，ｚｎ  ）を計算する第９のステップと、線量分布Ｄ（ｘ，ｚｎ  ）を前記すべてのｎに対して計算する第１０のステップとを実行するものである。
【００１５】
この発明の第２の方法に係る粒子線線量分布の計算方法は、粒子線として陽子線、電子線、または重粒子線を選んだものである。
【００１６】
この発明の第３の方法に係る粒子線線量分布の計算方法は、発散を以下の式で逐次計算したものである。
【００１７】
【数７】
0003620227.tif 000006
【００１８】
この発明の第４の方法に係る粒子線線量分布の計算方法は、２乗平均半径を以下の式で逐次計算したものである。
【００１９】
【数８】
0003620227.tif 000007
【００２０】
【発明の実施の形態】
実施の形態１．
図１はこの発明の実施の形態１による粒子線線量分布の計算方法を示すフローチャートである。まず、各ステップについて説明する。
・ステップ１：入射するガウスビームのエネルギーＥ０  と粒子数Ｎ（Ｅ０  ）、２乗平均散乱角と発散と２乗平均半径の各初期値、および媒質の物性値分布を入力する。また、ｎ＝１、ｚ０  ＝０、粒子数ウェイトＷ（ｚ０  ）＝１とする。
・ステップ２：エネルギー減衰率が０．０１となる厚さの増分Δｚｎ  と減衰後のエネルギーＥｎ  と表面からの深さｚｎ  を計算する。
・ステップ３：核反応による粒子数ウェイトＷ（ｚｎ  ）を計算する。
・ステップ４：平均エネルギー損失Ｗ（ｚｎ  ）ｄＥ／ｄｚ（ｚｎ  ）を計算する。
・ステップ５：下記の式で表される２乗平均散乱角を計算する。
【００２１】
【数９】
0003620227.tif 000008
【００２２】
・ステップ６：下記の式で表される発散を計算をする。
【００２３】
【数１０】
0003620227.tif 000009
【００２４】
・ステップ７：下記の式で表される２乗平均半径を計算する。
【００２５】
【数１１】
0003620227.tif 000010
【００２６】
・ステップ８：エネルギーＥｎ  ＜０．１ＭｅＶの判定を行う。
・ステップ９：平均エネルギー損失に対して、レンジストラグリングと入射ビームのエネルギー分散を考慮して実効エネルギー損失Ｓ（Ｅ０  ，ｚｎ  ）を計算する。
・ステップ１０：線量分布Ｄ（ｘ，ｚｎ  ）を求める（ｎ＝１，２，．．．）。
【００２７】
次に、計算の詳細について説明する。ステップ１は粒子線たとえば陽子線が、人体のような媒質に照射される時の初期条件を入力するステップである。媒質の物性値分布とは、例えば人体を１立方センチメートルの立方体要素に分割して、各要素内の平均密度、平均原子番号、平均質量数、平均電離エネルギーなどで代表される分布を指す。以下の計算は媒質を多層平行平板で離散的に近似して、平行平板の厚さ方向の座標ｚをｚｎ  （ｎ＝０，１，２，．．．）として、各ｎに対して逐次計算を進めていく（図２）。このようにすれば、人体のような不均質な媒質も扱うことができる。ステップ２はエネルギーの減衰が１％になるように次の層の厚さの増分Δｚｎ  を決めて、Δｚｎ  だけ粒子が進行した位置におけるエネルギーと表面からの深さｚｎ  を計算する。なお、ここで１％としたのは、解の収束状況から最適な値として決定したが、さらに高速計算を望む場合は多少の精度低下を許せば、１％ないし１０％程度にしてもよい。エネルギーの減衰率を大きくすれば、厚さの増分Δｚｎ  が大きくなり、計算は高速化する。したがって、減衰率計算精度と計算速度の優先度に応じて任意に選択すればよいが、本計算法によれば実用上有効な範囲として、０．１％ないし１０％とすることができる。ステップ３は粒子が輸送中に媒質内で核反応を起こして、消滅することを考慮した公知の補正係数の計算である。粒子数ウェイトは最初は１であるが、粒子の消滅が起こると１より小さくなる。例えば水中３０ｃｍの輸送後、粒子が停止する場合は、最終的に０．７程度の値になる。この値は、陽子が水の中の酸素原子核に衝突して核反応する場合の反応断面積データに基づいて計算する方法が公知である。つまり、１００００個の粒子を照射すれば、停止位置３０ｃｍ深さでは平均的には７０００個に減ってしまう。ステップ４の平均エネルギー損失を計算する場合は、この粒子数ウェイトを掛けておく必要がある。ステップ５は２乗平均散乱角の計算である。χｃ２（ｚｎ  ）とＢ（ｚｎ  ）は散乱に関する公知のパラメータであり、この式は公知の式である。ステップ６とステップ７は発散および２乗平均半径の計算であり、従来例の式（３）、（４）に比べて、はるかに簡単な式で計算できる。ステップ６とステップ７で計算に用いる式は、前記文献に記載されている以下の式（５）〜（７）をｚに関して微分して、式変形を繰り返し行うことにより得ることができた新規な式である。
【００２８】
【数１２】
0003620227.tif 000011
【００２９】
ここに、式（５）の右辺の第１項は２乗平均散乱角の初期値であり、第２項の積分は媒質内での散乱角の増加分を与える。また、式（６）はビームの発散具合を表す。右辺の第１項はビームの発散の初期値である。第２項は２乗平均散乱角の初期値に起因する項、第３項は媒質内の散乱による効果である。式（７）は２乗平均ビーム半径を表し、右辺の第１項は２乗平均ビーム半径の初期値、第２項は発散の初期値による項、第３項は２乗平均散乱角の初期値による項、第４項は媒質内の散乱による項である。
ステップ８は粒子のエネルギーが規定値、例えば０．１ＭｅＶまで低下したかを判定する。０．１ＭｅＶに到達していない場合は、ｎを１つ増加させてステップ２に戻る。到達したら平均エネルギー損失をｎ＝１，２，．．．に対して並べて、レンジストラグリングと入射ビームのエネルギー分散を考慮して実効エネルギー損失Ｓ（Ｅ０  ，ｚｎ  ）を計算する（ステップ９）。この計算は次式のコンボリューションで与えられる。通常、ｓは陽子線の飛程の１ないし２％程度の値になる。
【００３０】
【数１３】
0003620227.tif 000012
【００３１】
ここに、ｈは平均エネルギー損失をあらわし、ｓはレンジストラグリングと入射ビームのエネルギー分散をガウス分布で近似したときの標準偏差である。最後に、ステップ１０では各ｚｎ  に対して式（１）を適用して線量分布Ｄ（ｘ，ｚｎ  ）を求める。
【００３２】
図３は、陽子線入射エネルギー１７７．２５ＭｅＶにおいて図１の計算方法で得られた深さ方向の実効エネルギー損失Ｓ（Ｅ０  ，ｚ）、すなわち相対線量と実験データ（黒丸）を媒質が水の場合に比較した一例である。図４は、陽子線入射エネルギー９４ＭｅＶにおいて、図１の計算方法で得られた２乗平均ビーム半径の平方根と実験データ（黒丸）を比較した一例である。ただし、この図では実験的に得られたビーム半径から式（７）の第１ないし３項の項を引き算補正して、第４項のみを抽出して計算結果と比較している。図３および図４は本発明の計算方法の妥当性を示している。従来例の式に比べて、本発明による２乗平均ビーム半径を計算するステップ５、６、７の式は積算・除算の回数が著しく少ないので、計算を高速化することができる。従って、従来の計算法と比べて目標精度を同程度とした場合、数倍の計算速度が得られる。
【００３３】
上記の一連のステップはこの発明を実施する一例に過ぎず、数学的に等価な操作に置き換えても同様な効果が得られる。例えば、ステップ１でｎ＝０として計算を開始して、計算ループの先頭位置（ステップ２の直前）でｎを１つ増加させてもよい（図５）。また、ループ計算を終了する判定基準の０．１ＭｅＶは前記文献で実用上妥当な値として採用された値であるが、目的に応じて変更してもよい。さらに、ステップ５の右辺の第１項のかわりに係数が多少異なる前記文献の１０３３頁記載の式（５）の右辺の２乗
【００３４】
【数１４】
0003620227.tif 000013
【００３５】
を用いてもよい。これらは、いずれも実験データを用いてカーブフィッティングしており、両者は係数が多少異なるが、いずれも同等の近似計算精度を有する。
【００３６】
【発明の効果】
以上のように、この発明の第１から第４の粒子線線量分布の計算方法によれば、媒質を微小な厚さΔｚｎ  （ｎ＝１，２，３，．．．）を有する多数の層に分割し、入射するガウスビームのエネルギーＥ０  と粒子数Ｎ（Ｅ０  ）、２乗平均散乱角と発散と２乗平均半径の各初期値、および媒質の物性値分布を入力し、番号ｎを１とする第１のステップと、エネルギー減衰率が規定値となる厚さの増分Δｚｎ  とエネルギー減衰後のエネルギーＥｎ  と表面からの深さｚｎ  を計算する第２のステップと、核反応による粒子数ウェイトＷ（ｚｎ  ）を計算する第３のステップと、平均エネルギー損失を計算する第４のステップと、２乗平均散乱角を計算する第５のステップと、発散を計算する第６のステップと、２乗平均半径を計算する第７のステップと、エネルギーＥｎ  が規定値以下になったかを判定する第８のステップとを有し、前記エネルギーＥｎ  が規定値に達するまでは前記ｎを１つ増加させて第２から第７のステップを繰り返し実行し、規定値に達した後は計算したすべてのｎの前記平均エネルギー損失に対して、レンジストラグリングと入射ビームのエネルギー分散を考慮して実効エネルギー損失Ｓ（Ｅ０  ，ｚｎ  ）を計算する第９のステップと、線量分布Ｄ（ｘ，ｚｎ  ）を前記すべてのｎに対して計算する第１０のステップとを実行するようにしたので、計算を高精度で高速化することができる。
【図面の簡単な説明】
【図１】この発明の実施の形態１による粒子線線量分布の計算方法を示すフローチャートである。
【図２】この発明の実施の形態１による粒子線線量分布の計算方法を説明する説明図である。
【図３】この発明の実施の形態１による粒子線線量分布の計算方法の妥当性を示す図である。
【図４】この発明の実施の形態１による粒子線線量分布の計算方法の妥当性を示す図である。
【図５】この発明の実施の形態１による粒子線線量分布の計算方法を示す他のフローチャートである。
【図１】
0003620227.tif 000016
【図２】
0003620227.tif 000017
【図３】
0003620227.tif 000018
【図４】
0003620227.tif 000019
【図５】
0003620227.tif 000020