"多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を行って評価指標値を導出し、当該評価指標値と予め設定された閾値とを比較して前記照射計画を作成する照射計画作成方法であって、前記逐次近似繰り返し演算は、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値を導出することを特徴とする照射計画作成方法。請求の範囲第1項に記載の照射計画作成方法であって、前記第3条件は、前記ペンシルビームの進行方向と平行なz方向と、前記z方向に対して垂直な面に平行な一の方向であるx方向と、前記垂直な面に平行であり、かつ前記x方向に対して垂直な方向であるy方向と、の各方向について設定される予測誤差に比例したペナルティー係数を用いることを特徴とする照射計画作成方法。請求の範囲第1項に記載の照射計画作成方法であって、前記逐次近似繰り返し演算は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を行うことを特徴とする照射計画作成方法。請求の範囲第3項に記載の照射計画作成方法であって、前記第4条件は、照射するペンシルビームが前記スポットに高い線量を照射する確率とともに増大するように設定したリスク指標を用いることを特徴とする照射計画作成方法。請求の範囲第1項から請求の範囲第4項のうちのいずれか一項に記載の照射計画作成方法であって、前記照射部位および前記非照射部位と、照射門数および前記照射方向と、に基づいて、前記照射対象に対するペンシルビームの照射位置を決定するペンシルビーム照射位置決定ステップと、前記ペンシルビーム照射位置決定ステップで決定した前記照射位置で重ね合わせるペンシルビームの線量分布を作成するペンシルビーム線量分布作成ステップと、前記ペンシルビームの重みが入力される重み入力ステップと、前記逐次近似繰り返し演算により前記評価指標値を導出する評価指標値導出ステップと、前記評価指標値導出ステップで導出した前記評価指標値が前記閾値よりも低いか、または、前記評価指標値導出ステップに設定された管理値が予め設定された制限値を超えた場合は、前記評価指標値の導出を終了するとともに、前記入力されたペンシルビームの重みを照射パラメータとして出力する出力ステップと、前記評価指標値導出ステップで導出した前記評価指標値が前記閾値以上、かつ、前記評価指標値導出ステップに設定された前記管理値が予め設定された前記制限値以下の場合は、重みを異なる数値に更新し、再度、前記重み入力ステップに入力する重み更新ステップと、を含むことを特徴とする照射計画作成方法。請求の範囲第5項に記載の照射計画作成方法であって、前記ペンシルビーム線量分布作成ステップが、前記ペンシルビームの線量分布の作成とともに、前記ペンシルビームの進行方向の線量分布に関する偏微分フィルターおよび当該進行方向に対して垂直な方向の線量分布に関する偏微分フィルターを作成し、前記評価指標値導出ステップが、作成したこれらの偏微分フィルターを用いて前記評価指標値を導出することを特徴とする照射計画作成方法。多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を行って評価指標値を導出し、当該評価指標値に基づいて前記照射計画を作成する照射計画作成装置であって、前記逐次近似繰り返し演算を行う逐次近似繰り返し演算手段が、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値の導出することを特徴とする照射計画作成装置。請求の範囲第7項に記載の照射計画作成装置であって、前記逐次近似繰り返し演算手段は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を行うことを特徴とする照射計画作成装置。多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、コンピュータを、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を実行させて評価指標値を導出し、当該評価指標値に基づいて前記照射計画を作成する逐次近似繰り返し演算手段として機能させる照射計画作成プログラムであって、前記逐次近似繰り返し演算手段は、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値の導出を実行させることを特徴とする照射計画作成プログラム。請求の範囲第9項に記載の照射計画作成プログラムであって、前記逐次近似繰り返し演算手段は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を実行させることを特徴とする照射計画作成プログラム。しかしながら、前記したように10000を超えるペンシルビームスポットを有する照射計画において、非特許文献1に提案された手法でロバストな計画を得るための計算を行うと、数日例えば2~3日というような、実用に耐えない長い計算時間を要するという問題がある。そのため、研究レベルのものは別として、市販の照射計画作成装置にロバストな計画を得るための手法が適用された例はない。本発明は前記問題に鑑みてなされたものであり、従来よりも短い時間で比較的ロバストな計画を作成決定することのできる照射計画作成方法、照射計画作成装置、および照射計画作成プログラムを提供することを課題とする。前記課題を解決した本発明は、次の1~10のとおりである。1本発明は、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を行って評価指標値を導出し、当該評価指標値と予め設定された閾値とを比較して前記照射計画を作成する照射計画作成方法であって、前記逐次近似繰り返し演算は、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値を導出することを特徴とする。2本発明は、前記1に記載の照射計画作成方法であって、前記第3条件は、前記ペンシルビームの進行方向と平行なz方向と、前記z方向に対して垂直な面に平行な一の方向であるx方向と、前記垂直な面に平行であり、かつ前記x方向に対して垂直な方向であるy方向と、の各方向について設定される予測誤差に比例したペナルティー係数を用いるのが好ましい。3本発明は、前記1に記載の照射計画作成方法であって、前記逐次近似繰り返し演算は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を行うのが好ましい。4本発明は、前記3に記載の照射計画作成方法であって、前記第4条件は、照射するペンシルビームが前記スポットに高い線量を照射する確率とともに増大するように設定したリスク指標を用いるのが好ましい。5本発明は、前記1から前記4のうちのいずれか一つに記載の照射計画作成方法であって、前記照射部位および前記非照射部位と、照射門数および前記照射方向と、に基づいて、前記照射対象に対するペンシルビームの照射位置を決定するペンシルビーム照射位置決定ステップと、前記ペンシルビーム照射位置決定ステップで決定した前記照射位置で重ね合わせるペンシルビームの線量分布を作成するペンシルビーム線量分布作成ステップと、前記ペンシルビームの重みが入力される重み入力ステップと、前記逐次近似繰り返し演算により前記評価指標値を導出する評価指標値導出ステップと、前記評価指標値導出ステップで導出した前記評価指標値が前記閾値よりも低いか、または、前記評価指標値導出ステップに設定された管理値が予め設定された制限値を超えた場合は、前記評価指標値の導出を終了するとともに、前記入力されたペンシルビームの重みを照射パラメータとして出力する出力ステップと、前記評価指標値導出ステップで導出した前記評価指標値が前記閾値以上、かつ、前記評価指標値導出ステップに設定された前記管理値が予め設定された前記制限値以下の場合は、重みを異なる数値に更新し、再度、前記重み入力ステップに入力する重み更新ステップと、を含むのが好ましい。6本発明は、前記5に記載の照射計画作成方法であって、前記ペンシルビーム線量分布作成ステップが、前記ペンシルビームの線量分布の作成とともに、前記ペンシルビームの進行方向の線量分布に関する偏微分フィルターおよび当該進行方向に対して垂直な方向の線量分布に関する偏微分フィルターを作成し、前記評価指標値導出ステップが、作成したこれらの偏微分フィルターを用いて前記評価指標値を導出するのが好ましい。7本発明は、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を行って評価指標値を導出し、当該評価指標値に基づいて前記照射計画を作成する照射計画作成装置であって、前記逐次近似繰り返し演算を行う逐次近似繰り返し演算手段が、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値の導出することを特徴とする。8本発明は、前記7に記載の照射計画作成装置であって、前記逐次近似繰り返し演算手段は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を行うのが好ましい。9本発明は、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するにあたり、コンピュータを、入力されたペンシルビームの重みについて、予め入力された前記照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられるという第1条件と、予め入力された前記非照射部位に対しては、予め入力された線量制限以下に抑えるという第2条件と、を含む逐次近似繰り返し演算を実行させて評価指標値を導出し、当該評価指標値に基づいて前記照射計画を作成する逐次近似繰り返し演算手段として機能させる照射計画作成プログラムであって、前記逐次近似繰り返し演算手段は、前記第1条件および前記第2条件に加えてさらに、予め入力された照射方向から前記照射対象に与えられる前記線量分布の勾配を抑制するという第3条件を用いて前記評価指標値の導出を実行させることを特徴とする。10本発明は、前記9に記載の照射計画作成プログラムであって、前記逐次近似繰り返し演算手段は、前記第1条件、前記第2条件および前記第3条件に加えてさらに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという第4条件を用いて前記評価指標値の導出を実行させるのが好ましい。本発明に係る照射計画作成方法、照射計画作成装置、および照射計画作成プログラムのいずれによっても、従来よりも短い時間で比較的ロバストな計画を作成決定することができる。以下に、適宜図面を参照して本発明に係る照射計画作成方法、照射計画作成装置、および照射計画作成プログラムの一実施形態について詳細に説明する。[照射計画作成方法]まず、本発明の一実施形態に係る照射計画作成方法について説明する。本発明の一実施形態に係る照射計画作成方法は、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う腫瘍などの照射部位および放射線の照射を行うべきでないOARなどの非照射部位を含んでなる照射対象に対して予め入力された線量処方となるように線量分布を形成させるための照射計画を作成決定する。照射計画の作成決定は概ね、図1に示すように、情報入力ステップS1、ペンシルビーム照射位置決定ステップS2、ペンシルビーム線量分布作成ステップS3、重み入力ステップS4、評価指標値導出ステップS5と、出力ステップS6という手順で行われる。前記した逐次近似繰り返し演算による評価指標値の導出は、評価指標値導出ステップS5で行われ、導出した評価指標値と、予め設定された閾値とを比較して照射計画を作成するものである。評価指標値導出ステップS5の詳細は後に説明する。情報入力ステップS1情報入力ステップS1は、照射計画を得るために必要な情報を入力するステップである。照射計画を作成するために必要な情報としては、例えば、照射部位および非照射部位と、照射門数および照射方向と、前記した照射部位に対する目標線量および前記した非照射部位に対する線量制限すなわち、線量処方などがある。このような照射計画を作成するために必要な情報は、照射計画者によって入力される。照射計画者は、予め標的となる照射対象が撮影されたCT画像等を基に医学的所見に基づいてペンシルビームを照射すべき部位照射部位、ペンシルビームを照射すべきではない部位非照射部位を特定し、特定した照射部位に対して効果的と思える照射門数およびその照射方向を決定し、さらに、照射部位に対する目標線量と非照射部位に対する線量制限とを決定し、情報として照射計画作成装置電子計算機に入力する。照射部位は、一般にPTVとして設定されるものであり、例えば、腫瘍等のできた部位および腫瘍が浸潤している可能性のある部位等を含めた領域が設定される。また、非照射部位は、OARとして設定されるものであり、例えば、脳や視神経等の重要な臓器等が設定される。照射部位に対する目標線量および非照射部位に対する線量制限は、照射計画者により適宜設定される。ペンシルビーム照射位置決定ステップS2ペンシルビーム照射位置決定ステップS2では、予め入力された照射部位および非照射部位と、照射門数および照射方向と、に基づいて、照射対象に対するペンシルビームの照射位置を決定する。ペンシルビーム線量分布作成ステップS3次に行うペンシルビーム線量分布作成ステップS3では、ペンシルビーム照射位置決定ステップS2で決定したペンシルビームの照射位置で重ね合わせるペンシルビームの線量分布を作成する。つまり、このステップではペンシルビーム線量分布カーネルを作成する。なお、ペンシルビーム照射位置決定ステップS2とペンシルビーム線量分布作成ステップS3における前記した内容は、公知の重粒子線スキャニング照射用治療計画計算エンジンを用いることにより行うことができる。なお、このペンシルビーム線量分布作成ステップS3では、ペンシルビーム線量分布カーネルの作成とともに、ペンシルビームの進行方向であるz方向の線量分布に関する偏微分フィルターと、このz方向に対して垂直な面上におけるx方向およびy方向の線量分布に関する偏微分フィルターを作成する。つまり、ペンシルビームの軸方向および当該軸方向に垂直な横方向についての、線量勾配に関するカーネル線量勾配カーネルを作成する。ペンシルビーム線量分布作成ステップS3で作成した線量勾配カーネルを後記する評価指標値導出ステップS5で用いるようにすれば、後記するように、評価指標値導出ステップS5で繰り返し行われる線量重畳積分の演算ループを用いて各門からターゲット照射対象に与えられる線量分布の勾配を導出することができ、評価指標値導出ステップS5における評価指標値を効率良く導出することができる。評価指標値導出ステップS5で繰り返し行われる演算ループを用いるということは、新たな演算ループを追加するわけではないので、各門から照射対象に与えられる線量分布の勾配を導出するにあたって計算時間の著しい延長が生じないことを意味する。なお、線量勾配カーネルの作成については後に詳述する。重み入力ステップS4次に行う重み入力ステップS4では、ペンシルビームの重みが入力される。重みの初期値は、情報入力ステップS1で予め入力された照射部位に対する目標線量および前記非照射部位に対する線量制限に基づいて任意に設定される。後に詳述するように、より良い解を得るためペンシルビームの重みは、後記する評価指標値導出ステップS5で導出される評価指標値が特定の条件を満たすまで新たな値に更新されて演算が繰り返される。評価指標値導出ステップS5次に行う評価指標値導出ステップS5については、出力ステップS6の説明を行った後に説明する。出力ステップS6出力ステップS6では、評価指標値導出ステップS5で決定または更新したペンシルビームの重みを照射パラメータとして出力する。評価指標値導出ステップS5の一態様について先に言及した評価指標値導出ステップS5の一態様としては、重み入力ステップS4で入力したペンシルビームの重みについて、ペンシルビーム線量分布作成ステップS3で作成したペンシルビーム線量分布カーネルと、好ましくは同ステップで作成した線量勾配カーネルと、を用いて線量重畳積分および偏微分重畳積分を行い、次の第1条件および第2条件に加えて、さらに、次の第3条件を用いて逐次近似繰り返し演算を行い、評価指標値を導出し、当該評価指標値と予め設定された閾値とを比較して照射計画を作成することが挙げられる。第1条件:入力されたペンシルビームの重みについて、予め入力された照射部位に対しては、予め入力された目標線量に対して必要十分な線量が与えられる。第2条件:予め入力された非照射部位に対しては、予め入力された線量制限以下に抑えることができる。第3条件:予め入力された照射方向から照射対象に与えられる線量分布の勾配を抑制して評価指標値を導出することができる。つまり、評価指標値導出ステップS5の一態様では、線量勾配を抑制するように評価指標値を導出し、これにより、よりロバストな照射計画を作成する。前記した評価指標値導出ステップS5の一態様を具現する数式については、後ほど説明する。前記した一態様による評価指標値導出ステップS5は、図2に示すように、導出した評価指標値Eが、予め設定された閾値Cよりも低いか、または、評価指標値導出ステップS5に設定された管理値nが予め設定された制限値Nを超えた場合判定ステップS51においてYes、前記した評価指標値Eの導出を終了するとともに、入力されたペンシルビームの重みを照射パラメータとして出力するため、前記した出力ステップS6に移行する。ここで、管理値nとしては、評価指標値導出ステップS5における計算の繰り返し回数や、計算に要した時間などが挙げられる。制限値Nとしては、計算の繰り返し回数の上限値として予め設定された回数や、計算時間の上限として予め設定された時間などが挙げられる。閾値Cおよび制限値Nは、照射計画者によって任意に設定することができる。他方、導出した評価指標値Eが、予め設定された閾値C以上、かつ、評価指標値導出ステップS5に設定された管理値nが予め設定された制限値N以下の場合判定ステップS51においてNo、重み更新ステップS52で管理値nを更新するとともにn=n1、ペンシルビームの重みを更新して重み入力ステップS4に戻り、再度、評価指標値導出ステップS5を行う。ペンシルビームの重みの更新は、後記する式7に基づいて行うとよい。評価指標値導出ステップS5の一態様を具現する数式について評価指標値導出ステップS5で導出する評価指標値Eは、式1で求めることができる。なお、式1の目的関数fが評価指標値Eに相当する。式1において、wlは、各照射野lのビームの重みであり、wは、各照射野lのビームの重みwlの合計である。Tは、ターゲットである照射部位を示し、Oは、OARなどの非照射部位を示す。また、H’は、H’r=rHrと表され、との間に示される数式rの計算値がゼロより大きいときは1、ゼロ以下のときは0と定義されるヘヴィサイドステップ関数である。QToは、ターゲットである照射部位に指定したオーバードーズ過大線量に対するペナルティー係数影響度係数であり、Diは、位置iにおける全ての照射野lからの合計線量であり、DTmaxは、ターゲットである照射部位に照射される最大線量である。QTuは、ターゲットである照射部位に指定したアンダードーズ過小線量に対するペナルティー係数であり、DTminは、ターゲットである照射部位に照射される最小線量である。QOoは、OARなどの非照射部位におけるオーバードーズに対するペナルティー係数であり、DOmaxは、OARなどの非照射部位に許容される最大線量である。Nfieldは、入射される照射野の数であり、Dilは、ビームの重みwlで照射野lから位置iに照射される線量である。zは、z方向であり、xは、x方向であり、yは、y方向である。そして、Qx、Qy、Qzは、各方向の線量勾配へのペナルティー係数である。かかるペナルティー係数は、その方向への予測誤差に比例して設定される。これらの記号は、以下に示す各数式においても同じ内容を意味する。当然ながら、DiとDilは、式2の関係を満たす。式1の第一項は前記第1条件に相当するものであり、PTVへの線量制限目標線量を記述したものである。第二項は前記第2条件に相当するものであり、OARへの線量制限を記述したものである。そして、第三項は前記第3条件に相当するものであり、照射野単位におけるPTV内での線量勾配制限を記述したものである。線量分布Diは、ペンシルビームの重ね合わせで表され、式3で表現できる。式3において、Nspotは、照射される全スポット数であり、各照射野lでのスポットとは式4の関係を満たす。なお、jは、j番目に照射されるビームであることを示す。そして、式2および式4を用いて式3を展開すれば、式5となる。そして、式2~5を用いて式1を展開すれば、式6となる。式6によれば、各照射野lでの線量勾配抑制項第三項が、ペンシルビームの線量応答dijの偏微分の重畳積分という形で表現されている。すなわち、評価指標値導出ステップS5で用いられる線量分布計算の重畳積分と同じ演算ループで計算できることがわかる。ここで、スキャニング照射法における線量最適化の問題とは、式6で与えられる目的関数fを最小にするような3次元でのスポット強度分布線量分布を具現するビームの重みwを非負の条件のもとで求めることに帰着する。目的関数fの最小値を導出する方法として、逐次近似法を用いるのが一般的である。逐次近似法で最適な線量分布を求めるにはk回目の近似値wkから次の近似値wk1を式7で求める。なお、式7において、hとμは、それぞれ修正ベクトルの方向ベクトルとその長さ大きさを表す。効率よく目的関数fに最小値を与えるビームの重みwを求めるためには、hとμの求め方が鍵となる。なお、方向ベクトルhの計算アルゴリズムには、最急降下法moststeepestdescentmethod、共役勾配法conjugategradientmethod、ニュートン法Newtonmethod、改良ニュートン法modifiedNewtonmethod、擬似ニュートン法quasiNewtonmethodなどがある。線量最適化における逐次近似解法で方向ベクトルhを求めるためのいずれのアルゴリズムにおいても、少なくとも目的関数fの傾きを求める必要がある。本発明の場合、目的関数fの傾きは、式8で求めることができる。式8によれば、式6と同様に、線量勾配抑制項第三項には、ペンシルビームの線量応答dijの偏微分の重畳積分が現れるのみである。従って、ペンシルビーム線量分布作成ステップS3における線量分布計算の重畳積分と同じ演算ループで計算できることがわかる。各ペンシルビームの線量応答dijは、積分線量分布を与えるペンシルビームのz方向と、z方向に対して垂直方向であるx方向およびy方向に、式9のように因数分解することができる。ここで、式9のdxとdyは、規格化したガウス分布としてそれぞれ式10、式11のように表現することができる。式10、式11において、x方向およびy方向への拡がりを表すσxとσyは、レンジシフター厚と深さの関数としてパラメータ化し、線量計画の線源データ図3参照として組み込まれる。計算機内での線量重畳積分では、拡がりσの関数として予めこれをメモリテーブル化し、ガウスフィルターとして記憶手段に記憶させておき、これをCT画像内でのスポット位置に応じて導出されるσx、σyのインデックスとして参照できるようにすることで適宜に処理することができるようになる。一方、積分線量分布を表すz方向成分は、例えば、平行平板電離箱での測定値を補完して線源データ図4参照として記憶手段に記憶させておくとよい。これらを踏まえ、式6に現れるペンシルビームの線量応答dijのx方向およびy方向に対する偏微分を求めると、式12、式13のようになる。ここで、式12、式13に現れるは、ガウス分布の偏微分である。従って、例えば、図5Aに示すガウス分布は、図5Bに示すような関数形状のガウス偏微分分布となる。つまり、式6、式8に現れるペンシルビームのx方向およびy方向への偏微分の重畳積分は、線量重畳積分での着目方向についてのガウスフィルターを図5Bに示したようなガウス分布の偏微分フィルターに置き換えるのみでよいことになる。従って、照射計画の計算コードでは、拡がりσの関数として、ガウスフィルターと同時にガウス偏微分フィルターをメモリテーブル化しておけばよいことになる。なお、このガウス偏微分フィルターが線量勾配フィルター線量勾配カーネルとなる。他方、ペンシルビームのz方向への偏微分は、式14により求めることができる。式14の第一項および第二項については、図3に示したとおり、x方向およびy方向への拡がりσx、σyは、深さとレンジシフターの関数としてパラメータ化され、逐次近似法を用いた重粒子線スキャニング照射用治療計画計算エンジン例えば、IDOSEに組み込まれている。従って、式14で初めて現れる関数は、x方向およびy方向への拡がりの深さ方向への偏微分である。この関数形を概観すれば、図6Aのような拡がりに対して、図6Bのようになる。重粒子線スキャニング照射用治療計画計算エンジンには、線源データとしてx方向およびy方向への拡がりσx、σyが、深さとレンジシフターの関数としてパラメータ化されている。従って、これらの偏微分は、式15、式16として容易に計算することができる。これらは、各ペンシルビームに対してペンシルビーム線量分布カーネルを作成する際に、σx、σyと同じように、偏微分を割り当てることによって偏微分フィルター線量勾配カーネルを作成することができる。次に、式14の第三項をみれば、積分線量dzのz方向の偏微分が現れる。これは、図7Aに示すような積分線量dzに対して、図7Bのような関数形になる。x方向およびy方向への拡がりと同様、通常の線量計算と線量最適化を行うことができる。線量最適化は、重粒子線スキャニング照射用治療計画計算エンジンに積分線量dzが線源データとして登録されているので、各ペンシルビームに対してペンシルビーム線量分布カーネルを作成する際に式17に従って偏微分を割り振れば偏微分フィルター線量勾配カーネルを作成することができる。以上に説明したように、評価指標値導出ステップS5の一態様によれば、ペンシルビームの重みの更新を行う逐次近似法による演算ループと式1を備えることにより、より良い評価指標値E最小の目的関数fを具現するペンシルビームの重みを導出することができる。その結果として、重ね合わせるペンシルビームの線量勾配を抑制した、比較的ロバストな計画を計算時間の延長を招くことなく短時間で得ることができる。また、ペンシルビーム線量分布作成ステップS3で作成したペンシルビーム線量分布カーネルとその偏微分カーネルを用い、線量の重畳積分と偏微分の重畳積分を同一の演算ループで行うことにより、効率良く比較的ロバストな計画を導出することができる。図8aおよびbに、評価指標値導出ステップS5の一態様を適用した照射計画作成方法によって作成した結果の一例を等線量線で示すとともに、図9aおよびbに、ロバストでない従来の照射計画作成方法によって作成した結果の一例を等線量線で示す。図8のaおよびbでは、線量勾配抑制項第三項が用いられているので、二つの照射野が重なり合う部分の線量勾配が抑制されている。そのため、飛程の誤差やセットアップの誤差が生じても影響を受け難く、これらの誤差に強い。これに対し、図9のaおよびbでは、線量勾配抑制項第三項を用いていないので、二つの照射野が重なり合う部分に大きな線量勾配が現れている。そのため、飛程の誤差やセットアップの誤差が生じると影響を受け易く、これらの誤差に弱い。評価指標値導出ステップS5の他の態様そして、評価指標値導出ステップS5の他の態様は、前記した第1条件、第2条件および第3条件に加えて、さらに、次の第4条件を用いて逐次近似繰り返し演算を行い、評価指標値Eを導出し、当該評価指標値Eと予め設定された閾値Cとを比較して照射計画を作成する。第4条件:予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いない。なお、かかる誤差としては、後記するように、飛程の誤差およびセットアップの誤差のうちの少なくとも一方が挙げられる。つまり、評価指標値導出ステップS5の他の態様では、線量勾配を抑制するとともに、非照射部位に与える線量が高くなるスポットを用いないように評価指標値Eを導出し、これにより、さらにロバストな照射計画を作成する。評価指標値導出ステップS5の他の態様を具現する数式について次に、評価指標値導出ステップS5の他の態様について説明する。誤差の大きさは非等方性であり、該当する照射野に固有である。これらの誤差は前記したようにz方向に平行な成分とz方向と垂直な成分に分解される。飛程の誤差を招く不確実性以下飛程の不確実性という。は、z方向に平行な成分であることが多い。飛程の不確実性には複数の因子があり、例えば、CTのアーチファクト、患者の体重増加または体重減少、HounsfieldユニットHUから阻止能への変換などがある。飛程の不確実性は水等価深さに比例し、極端な例では、飛程の不確実性δRは±5.0になる場合もある。セットアップの誤差による不確実性以下セットアップの不確実性という。については、照射対象が各方向に等しく移動し得ることから、等方性の不確実性を考慮する必要がある。等方性のセットアップの誤差は、各照射野方向について独立的なセットアップの誤差を前提とすれば、z方向に平行な成分と、z方向に垂直な成分つまり、x方向およびy方向の両方の成分と、に分割することができる。z方向に平行なシフトは、患者の前面の空隙の増減を引き起こすが、その結果は線量分布の最小限の変動にとどまる。従って、セットアップの不確実性に関しては、各治療ビームに垂直なシフトを考慮すれば充分である。ここで、x方向のセットアップの不確実性ΔSxおよびy方向のセットアップの不確実性ΔSyについては、RutsHPandLomaxAJDonutshapedhighdoseconfigurationforprotonbeamradiationtherapyStrahlentherOncol2005181p4953に基づけば、それぞれ±3.0mmの誤差が想定される。前記したように、強度変調型の重粒子線治療の主な目的は、iCTVClinicalTargetVolume;臨床標的容積内に必要十分な線量を与え、iiOARへの望ましくない意図しない照射を防止してダメージを線量制限以下に抑えることである。ロバストな照射計画作成方法は、各々の目的のアルゴリズムを組み合わせることによって両方の目的を同時に達成する。前記iの目的には、InaniwaTKanematsuNFurukawaTandNodaKRobustdualfieldoptimizationofscannedionbeamsagainstrangeandsetupuncertaintiesPhysMedBiolsubmittedに記載された方法を用いることができる。前記iiの目的のアルゴリズムを以下に説明する。通常は、照射計画において隣接するOARに高い線量を照射するリスクのあるペンシルビームはごく一部である。IMITやIMPTでは、そのようなリスクのあるペンシルビームの線量寄与を少なくし、別の入射角度で照射されるペンシルビームによって補填することができる。そこで、線量分布の最適化においてOARに高い線量を照射するリスクのあるペンシルビームにペナルティーを科すことが考えられる。これを具現したのが、式18である。この式18は、評価指標値導出ステップS5の他の態様を具現する数式であり、前記式1に、予め設定される範囲内で誤差が生じた場合に、非照射部位に与える線量が高くなるスポットを用いないという第4条件を第四項として追加したものである。式18では、j番目の1つのペンシルビームのリスクを特徴づける量的尺度であるリスク指標Pjを線量最適化に導入している。式18において、Nlは、照射野lにおけるペンシルビームの数であり、wjは、j番目のペンシルビームの重みであり、wkは、k番目のペンシルビームの重みであり、QpbRとQpbSは、前記追加された項の重要度を調整するペナルティー係数であり、PjRは、j番目のペンシルビームの飛程の不確実性のリスク指標であり、PjSは、j番目のペンシルビームのセットアップの不確実性のリスク指標である。その他の記号は式1の例による。式18は、QpbR=QpbS=0とすることによって、飛程の誤差のみを考慮する照射計画を得ることができ、Qx=Qy=Qz=QpbR=QpbS=0とすることによって、飛程の誤差およびセットアップの誤差を考慮しない従来の照射計画を得ることができる。あるペンシルビームを用いることのリスクは、そのペンシルビームがOARに高い線量を照射する確率とともに増大する。そこで、j番目に照射される1つのペンシルビームのリスクを定量化するため、式18に示したように、PjRとPjSの2つのリスク指標を導入している。ここで、図10を参照してリスク指標PjRとリスク指標PjSの算出について説明する。図10は、馬蹄形のCTVが円柱状のOARを囲んでいるRTOGベンチマークファントムを用いて、リスク指標PjRとリスク指標PjSを算出する手順を概略的に示している。なお、図10は、aの太い矢印で示されているように、0°の角度でペンシルビームが照射されることを前提としている。つまり、太い矢印がz方向になり、図10の紙面上で当該z方向と垂直な方向がx方向になる。飛程の誤差とセットアップの誤差に起因するCTVの周辺領域への低線量の照射に対処するため、図10bに示すように、2つのリスク指標の算出における第1のステップとして、標的容積をCTVとその周辺領域を取り囲むボリュームに拡大し、FTVFieldspecificTargetVolume;照射野固有標的容積を生成する。ここで、図10bに示されるΔRは、飛程の誤差であり、CTV境界の水等価深さRと、想定される飛程の不確実性δRから、ΔR=R×δRとして算出できる。飛程の不確実性のリスク指標PjRの導出における第2のステップとして、‘飛程の不確実性に対しリスクのある該当照射野に固有の重要臓器ROAR’を生成する。ここでは、図10cに示すように、ROARはOARそのものである。そのため、ROAR全体にリスク係数F=1.0を適用する。第3のステップとして、図10dに示すように、リスク指標PjRの算出において、スキャンするペンシルビームの横方向幅を考慮して半影部penumbraをROARに追加する。半影部のリスク係数Fは、ROARの境界からの距離lの関数として次の式19によって決定する。式19において、σは、ブラッグピークにおけるペンシルビームの横方向幅である。σは、例えば、5mmなどとすることができるが、ペンシルビームの横方向幅に応じて任意に設定し得る。なお、式19において、距離lが負の値である場合は、ROARの内側を意味し、距離lが正の値である場合は、ROARの外側を意味する。次に、FTV全体にわたって各々のブラッグピークの位置を決定する。まず、リスク指標PjRを算出するため、飛程の誤差ΔRjをブラッグピーク位置毎に水等価深さRjに基づいてΔRj=Rj×δRにて算出する。また、飛程の誤差ΔRが存在するため、j番目のペンシルビームについて、OARに照射される正規化された線量としてリスク指標PjRを式20のように定義する。式20において、djは、z=0から1.0までの飛程を正規化したj番目のペンシルビームの線量の平面積分であり、Fは、ROARとその半影部に適用するリスク係数であり、zmaxは、z方向に沿った計算領域の最大の深さである。また、正規化定数2ΔRjは、ゼロでないFがz方向で取り得る長さであり、Djは、式21によって算出される拡大ブラッグピークである。なお、式21において、dは、飛程rのペンシルビームの平面積分線量であり、djは、dに等しい。また、式21においては、飛程の誤差は等しく、Rj-ΔRj≦r≦RjΔRjの範囲であると想定している。セットアップの不確実性のリスク指標PjSの算出における第2のステップとして、‘セットアップの不確実性に対しリスクのある該当照射野に固有の重要臓器SOAR’を生成する。セットアップの誤差ΔSが存在するため、リスク指標PjRの定義と同様に、j番目のペンシルビームについて、OARに照射される正規化された線量としてリスク指標PjSを定義する。従って、例えば、図10c’に示すように、OARがセットアップの誤差ΔSによって横方向x方向に拡大された領域としてSOARを定義する。さらに、リスク指標PjSを算出するにあたって、スキャンするペンシルビームの横方向幅を考慮し、図10d’に示すように、半影部をSOARに追加する。SOARとその半影部にリスク係数F’を適用する手順は、ROARとその半影部における手順と同じでよい。そのため、セットアップの不確実性のリスク指標PjSは、式22のように算出することができる。式22において、F’は、セットアップの不確実性のためにSOARとその半影部に適用されるリスク係数であり、VOARは、OARの容積である。また、式22の分母は、SOARを、F’がとり得るゼロではないz方向の長さに対応する容積VOARの球であると想定した場合における、当該SOARの平均コード長さを表す規格化定数である。以上に説明したように、評価指標値導出ステップS5の他の態様によれば、ペンシルビームの重みの更新を行う逐次近似法による演算ループとリスク指標を用いた式18を備えることにより、より良い評価指標値最小の目的関数fを具現するビームの重みを導出することができる。その結果として、重ね合わせるペンシルビームの線量勾配を抑制するとともに、非照射部位に与える線量が高くなるスポットを用いないようにした、さらにロバストな計画を計算時間の延長を招くことなく短時間で得ることができる。図11aおよびbに、評価指標値導出ステップS5の他の態様を適用した照射計画によって作成した結果の一例を示すとともに、図12aおよびbに、ロバストでない従来の照射計画によって作成した結果の一例を示す。図11のaおよびbでは、線量勾配抑制項第三項とともに、予め設定される範囲内で誤差が生じた場合に、前記非照射部位に与える線量が高くなるスポットを用いないという項第四項を用いている。そのため、二つの照射野が重なり合う部分の線量勾配が抑制され、また、OARの直前の線量分布が低線量に抑制されている。その結果、飛程の誤差やセットアップの誤差が生じても影響を受け難く、これらの誤差に強い。これに対し、図12のaおよびbでは、線量勾配抑制項第三項を用いていないので、二つの照射野が重なり合う部分の線量勾配が大きい。また、第四項を用いていないので、ビームの入射方向から見てOARの直前の線量分布が高線量となっている。そのため、飛程の誤差やセットアップの誤差が生じると影響を受け易く、これらの誤差に弱い。[照射計画作成装置]次に、本発明に係る照射計画作成装置について説明する。本発明に係る照射計画作成装置は、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して予め入力された線量処方となるように線量分布を形成させるための照射計画を作成する。本発明に係る照射計画作成装置は、コンピュータを前記したように機能させるため、評価指標値を算出するにあたって、既に説明した第1条件から第3条件または第1条件から第4条件に基づいた逐次近似繰り返し演算を行う。本発明の照射計画作成装置1の一実施形態としては、図13に示すように、ペンシルビーム照射位置決定手段2と、ペンシルビーム線量分布作成手段3と、重み入力手段4と、評価指標値導出手段5と、出力手段6と、を含んでなる。これらの各手段は、コンピュータにおいて、図示しないハードディスクドライブなどの記憶装置に格納された、後記する本発明に係る照射計画作成プログラムをCPUCentralProcessingUnitが実行し、本発明の照射計画作成方法で説明した各ステップの内容を実現することにより達成される。つまり、照射計画作成装置1の各手段は、前記した照射計画作成方法の各ステップに対応するものである。すなわち、ペンシルビーム照射位置決定手段2とペンシルビーム照射位置決定ステップS2とが対応し、ペンシルビーム線量分布作成手段3とペンシルビーム線量分布作成ステップS3とが対応し、重み入力手段4と重み入力ステップS4とが対応し、評価指標値導出手段5と評価指標値導出ステップS5とが対応し、判定手段51と判定ステップS51とが対応し、重み更新手段52と重み更新ステップS52とが対応し、出力手段6と出力ステップS6とが対応する。従って、各手段の内容についてのここでの詳細な説明は省略する。なお、照射計画作成装置1は、RAMやROM、ハードディスクドライブいずれも不図示などの、通常のコンピュータが有する各種装置を備えており、また、様々な表示を行う表示手段や、照射計画を得るために必要な情報を入力するためのキーボード、マウス、CT画像等を入力するための入力ポートといった情報入力手段不図示を備えている。[照射計画作成プログラム]本発明に係る照射計画作成プログラムは、コンピュータを、多門照射により強度を変調させつつペンシルビームをスキャンさせてスポット毎に重ね合わせ、放射線の照射を行う照射部位および放射線の照射を行うべきでない非照射部位を含んでなる照射対象に対して、予め入力された線量処方となるように線量分布を形成させるための照射計画を作成するように機能させる。本発明に係る照射計画作成プログラムは、コンピュータを前記したように機能させるため、評価指標値を算出するにあたって、既に説明した第1条件から第3条件または第1条件から第4条件に基づいた逐次近似繰り返し演算を行う。本発明の照射計画作成プログラムの一実施形態としては、コンピュータに、ペンシルビーム照射位置決定ステップ、ペンシルビーム線量分布作成ステップ、重み入力ステップ、評価指標値導出ステップ、判定ステップ、重み更新ステップ、出力ステップを実行させる。これらのステップは順に、本発明に係る照射計画作成方法で説明したペンシルビーム照射位置決定ステップS2、ペンシルビーム線量分布作成ステップS3、重み入力ステップS4、評価指標値導出ステップS5、判定ステップS51、重み更新ステップS52、出力ステップS6に相当するものである。従って、ここでの詳細な説明は省略する。本発明に係る照射計画作成プログラムは、CD-ROM、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体不図示に記録し、記録媒体駆動装置不図示によって、当該記録媒体から照射計画作成プログラムを読み出して図示しない記憶手段にインストールして実行するようにしてもよい。また、照射計画作成装置として機能するコンピュータクライアントが通信ネットワークなどの通信手段を備えている場合、本発明に係る照射計画作成プログラムが通信ネットワークを介して接続された他のコンピュータサーバに記憶され、当該コンピュータから通信ネットワークを介して照射計画作成プログラムをダウンロードして実行したり、サーバに記憶された本発明に係る照射計画作成プログラムを実行させたりして比較的ロバストな計画を得るようにしてもよい。また、この場合、数値解析した結果をサーバに備えられた記憶手段不図示に記憶することとしてもよい。本発明に係る照射計画作成方法、照射計画作成装置、および照射計画作成プログラムの有用性を示すべく、次のような検討を行った。[第一実施例]第一実施例では、IMRTのベンチマーク試験で用いられるRTOGRadiationTherapyOncologyGroupベンチマークファントム図14参照に対し、135°、225°コプラナーIMIT照射を従来法による照射計画すなわち、ロバストでない照射計画比較例1と、ペンシルビーム照射位置決定ステップS2と、ペンシルビームの線量分布の作成とともに、ペンシルビームの進行方向および当該進行方向に対して垂直な方向つまり、ペンシルビームの進行方向に対して垂直な面上の一の方向と、当該一の方向に対して垂直な方向の線量分布に関する偏微分フィルターを作成するペンシルビーム線量分布作成ステップS3と、重み入力ステップS4と、評価指標値導出ステップS5と、出力ステップS6とを含む、本発明の照射計画作成方法を適用した照射計画実施例1とを比較した。なお、計画を立案する前提として、馬蹄形のPTVの周りに、x,y,zの全ての方向にセットアップの誤差マージンを4mm設けた図15参照。図16に、第1条件と第2条件のみを用いる従来の照射計画作成方法比較例1で計画した結果を示す。図16Aは、従来の照射計画による合成臨床線量分布を示す図であり、図16Bは、図16Aに対応するPTVおよびOARの線量容積ヒストグラムDoseVolumeHistogramDVHを表す。また、図16C、図16Dは、225°、135°方向からそれぞれ与えられる線量分布を等線量線で表示した図である。図16C、図16Dに示すように、90、80および70ラインが非常に近接しており、二つの照射野が重なり合う部分に大きな線量勾配が現れている。そのため、飛程の誤差やセットアップの誤差に強くないことがわかる。一方、式1によって第1条件と第2条件に加えて、さらに第3条件を用いた本発明の照射計画作成方法実施例1で計画した結果を図17に示す。図17Aは、本発明の照射計画作成方法ロバストな照射計画作成方法実施例1で計画した合成臨床線量分布を示す図であり、図17Bは、図17Aに対応するPTVおよびOARのDHVである。図17Bに示すように、本発明の照射計画作成方法で計画すると、線量分布の悪化、すなわち、PTVでの線量不均一やOARへの線量付与増大が無いことがわかる。図17C、図17Dには、225°、135°方向から与えられる線量分布を示す。図17C、図17Dに示すとおり、本発明の照射計画作成方法で計画すると、二つの照射野が重なり合う部分でも線量勾配が抑制されており、飛程の誤差やセットアップの誤差に強いことがわかる。3.0GBのRAMを装着したDellPrecision690workstationを用いて照射計画の作成を試みたところ、第1条件から第3条件を用いる実施例1は2分30秒を要した。これに対し、第1条件と第2条件のみを用いる比較例1では同じスペックのコンピュータで2分を要した。このように、本発明の照射計画作成方法によれば、現在一般的に行われている従来法による照射計画ロバストでない照射計画に比して、大幅な計算時間の延長を招くことなく、比較的ロバストな計画を得ることができた。また、既に述べているとおり、非特許文献1に記載されたアプローチによってロバストな計画を得るためには数日例えば2~3日を要するところ、本発明の照射計画作成方法によれば、僅か2~3分で比較的ロバストな計画を作成決定することができた。[第二実施例]第二実施例では、極大照射野におけるパッチ照射の例として、照射対象である子宮がんに対し、垂直2門の前照射すなわち、照射対象である子宮がんに対して正面および側面からの前照射を行う場合の照射計画を行った。図18A、図18Bには、第1条件と第2条件のみを用いる線量勾配抑制項第3条件を用いない従来の照射計画作成方法によって作成したパッチ照射の照射計画比較例2における合成線量分布を示す。また、合成線量分布の基となった各ペンシルビームBeam1,Beam2の線量分布を図18C、図18Dに示す。両分布ともに等線量線が密になっており、急峻な線量勾配が生じていた。次に、セットアップの誤差により、Beam1とBeam2間の距離が4mm縮まった場合に生じる線量分布を図18E、図18Fに示す。このように各照射方向からのビームの間隔が縮まった場合には、2ヶ所に顕著なオーバードーズが認められ、その領域に115以上の高線量が付与されていた。この領域は、図18C、図18Dにおいて急峻な線量勾配が現れた位置に対応する。また、セットアップの誤差により、Beam1とBeam2間の距離が4mm離れた場合に生じる線量分布を図18G、図18Hに示す。このように各照射方向からのビームの間隔が離れた場合には、同じ領域で90以上のアンダードーズが生じていた。なお、各線量分布に対応するDVHを図18Iに示した。これに対し、ペンシルビーム照射位置決定ステップS2と、ペンシルビームの線量分布の作成とともに、ペンシルビームの進行方向および当該進行方向に対して垂直な方向の線量分布に関する偏微分フィルターを作成するペンシルビーム線量分布作成ステップS3と、重み入力ステップS4と、式1を用いた評価指標値導出ステップS5と、出力ステップS6とを含むつまり、第1条件から第3条件を含む本発明の照射計画作成方法を適用したパッチ照射の照射計画実施例2における合成線量分布を図19A、図19Bに示す。また、合成線量分布の基となった各照射方向からのビームBeam1,Beam2の線量分布を図19C、図19Dに示す。両分布ともに等線量線がほぼ均等に離れ、緩やかな線量勾配が生じていた。次に、セットアップの誤差により、Beam1とBeam2間の距離が4mm縮まった場合に生じる線量分布を図19E、図19Fに示す。このように各照射方向からのビームの間隔が縮まった場合には、105程度のオーバードーズが認められるが、従来法に比べ、オーバードーズが小さかった。また、セットアップの誤差により、Beam1とBeam2間の距離が4mm離れた場合に生じる線量分布を図19G、図19Hに示す。このように各照射方向からのビームの間隔が離れた場合には、同じ領域で95程度のアンダードーズが生じていたが、従来法に比較してアンダードーズの程度は小さかった。なお、各線量分布に対応するDVHを図19Iに示した。3.0GBのRAMを装着したDellPrecision690workstationを用いて照射計画の作成を試みたところ、第1条件から第3条件を用いる実施例2は2分30秒を要した。これに対し、比較例2では同じスペックのコンピュータで2分を要した。このように、本発明の照射計画作成方法によれば、現在一般的に行われている従来法による照射計画ロバストでない照射計画に比して、計算時間の延長を招くことなく、比較的ロバストな計画を得ることができた。また、既に述べているとおり、非特許文献1に記載されたアプローチによってロバストな計画を得るためには数日例えば2~3日を要するところ、本発明の照射計画作成方法によれば、僅か2~3分で比較的ロバストな計画を作成決定することができた。[第三実施例]第三実施例では、図20に示すRTOGベンチマークファントム様形状を対象にして、さらに、予め設定される範囲内で誤差が生じた場合に、非照射部位に与える線量が高くなるスポットを用いないという項第四項を用いたロバストな最適化方法の有効性について検討した。図20に示すRTOGベンチマークファントム様形状は、患者に見立てた直径232mm、高さ80mmの円柱体の中心に、直径30mm、高さ40mmの円柱状のOARが設けられている。かかるOARは、内径36mm、外径80mm、高さ40mmの馬蹄形のCTVによって囲まれている。かかるRTOGベンチマークファントム様形状は、均一で水に等価であるものと想定した。なお、ボクセルの大きさはΔx=Δy=Δz=2.0mmである。そして、図20中の破線で示すように315°、0°、45°で照射される3つのビームを検討した。すなわち、各照射野は、これらの角度から入射されるビームによって形成される。InaniwaTFurukawaTSatoSTomitaniTKobayashiMMinoharaSNodaKandKanaiTDevelopmentoftreatmentplanningforscanningirradiationatHIMACNuclInstrumMethodsPhysResB266p21942198に記載のペンシルビームのアルゴリズムを線量の計算に使用した。第三実施例では、146MeVの陽子ビームを用いるIMPTを想定した。146MeVの陽子ビームのブラッグピークの最大深さは水中で150mmである。第三実施例では、アイソセンターから600mm上流に設置したレンジシフターを使ってビームの飛程をシフトさせた。また、第三実施例におけるペンシルビームの幅は、レンジシフターがない状態で患者に見立てた円柱体の入射位置におけるσを4.0mmと想定し、レンジシフターの挿入によるビームの広がりを線量の計算時に考慮した。そして、各照射野について、ビーム進行方向z方向とビーム進行方向に対して垂直な面上の横方向x方向およびy方向に3mmの間隔で規則正しく並べられた方形グリッド上にブラッグピークを配置するようにした。前記した条件のもと、7つの異なる照射計画計画No.1~7を作成し、最適化のパラメータDTmax、DTmin、DOmax、QTo、QTuおよびQOoをそれぞれ、3.0Gy、3.0Gy、1.0Gy、1.0、2.0および2.0に固定した。計画No.1~7におけるその他の最適化のパラメータの概略を表1に示す。今回の検討では、ペナルティー係数であるQx、Qy、Qz、QpbRおよびQpbSを暫定的に、表1に示す値に設定した。表1に示されるQx、Qy、Qz、QpbRおよびQpbSは、前記式1および前記式18で説明した事項に関するパラメータである。δRは、飛程の不確実性として設定されるパラメータであり、ΔSxは、x方向のセットアップの不確実性として設定されるパラメータであり、ΔSyは、y方向のセットアップの不確実性として設定されるパラメータである。また、Tは、各計画の作成に要した計算時間であり、Npbは、各計画の作成に使用したスポットの総数である。第三実施例の結果および検討は、次のとおりである。1飛程の不確実性について前記したように、強度変調型の重粒子線治療の主な目的は、iCTVClinicalTargetVolume;臨床標的容積内に線量の範囲を維持し、iiOARへの望ましくない意図しない照射を防止することである。前記第四項を用いたロバストな照射計画作成方法における飛程の不確実性に対する感応度つまり、影響の受け易さを検討するため、表1中の計画No.1~3を比較した。計画No.1は、不確実性を全く考慮しない従来の照射計画である。計画No.2は、iiのみを考慮して飛程の不確実性に対処する照射計画である。そして、計画No.3は、iおよびiiを考慮して飛程の不確実性に対処する照射計画である。図21a~cは、計画No.1で作成決定された線量分布を示している。図21aおよびbは、0°および45°の各ビームの線量寄与を示す。なお、315°のビームからの線量寄与は、ファントムが対称であるため、45°のビームの線量の鏡像となるが、図21では図示を省略している。図21cは、計画No.1で計画された3つすべての照射野からの線量寄与を重ね合わせた合成線量分布を示している。図21cに示されているように、計画No.1は、飛程の誤差が全く生じない場合、CTVに対する極めて均一な線量分布とOARの良好な保全が実現できている。これは、最適化において、OARの直前に配置されたブラッグピークに高い重みを付与し、個々のビームによって、著しく不均一な線量分布が与えられている。つまり、計画No.1では、OARの直前に高線量が付与され、OARとCTVの間に急峻な線量勾配を形成し、これによってOARを良好に保全している。しかしながら、実際には飛程の不確実性誤差が存在するため、計画No.1で期待される特徴は低下する可能性がある。そこで、飛程の不確実性を分析するため、計画No.1を幾何学的な摂動とともに再計算した。なお、摂動とは、主要な運動にわずかな撹乱を加えること、およびそのための付加項をいう。再計算は、計画通りの密度を含めてファントムの実効密度を-5から5に、1ずつ飛程を意図的にずらして行った。これら11の再計算された線量分布に相当する線量容積ヒストグラムDVHを図22aに示す。飛程の不確実性に対する従来のIMPT計画計画No.1の感応度が、図22a中に明瞭に観察された。具体的には、CTVに対する著しいアンダードーズとともにOARに対する著しいオーバードーズが観察された。計画No.1の場合、CTV中の95線量D95は、54.81.64Gyであり、OARに照射された最大線量は、最も不良な場合には103.23.10Gyに達した。次いで、図21d~fは、計画No.2で作成決定された線量分布を示している。図21dおよび3eは、0°および45°の各ビームからの線量寄与を示す。なお、計画No.1の場合と同様、315°のビームからの線量寄与については図示を省略している。図21fは、計画No.2で計画された3つすべての照射野からの線量寄与を重ね合わせた合成線量分布を示している。計画No.2では、飛程の不確実性に対して高いリスク指標PjRを有するペンシルビームにペナルティーを科した。そのため、OARの直前に配置されるブラッグピークを有するペンシルビームは、線量最適化においては好適に回避されている。この計画No.2では、OARを保全するため、OARの直前で止まるように配置されるペンシルビームの末端部の線量減衰の代わりに、OARの横に配置されるペンシルビームの側面の線量減衰を使用している。また、計画No.2では、FTVを生成することによって、高い線量の領域をz方向にわずかに拡大し、飛程の不確実性に起因するCTVの周辺領域におけるアンダードーズを補っている。計画No.2におけるこれらの特徴は、飛程を拡大した場合のOARのオーバードーズと、飛程を縮小した場合のCTVのアンダードーズの両方を回避し得る。しかしながら、各ビームの線量分布は、z方向には不均一のままである。そのため、計画No.2についても、11の異なる実効密度となるよう、計画No.1と同様の摂動で線量分布を再計算した。図22bにこれらの線量分布に関するDVHを示す。図22bに示すように、飛程の不確実性に対する照射計画の感応度は、計画No.2ではかなり低減されているが、各ビームの不均一な線量分布のため、CTVのDVHは、計画No.2の計画通りに線量分布が形成される場合と比べるとやや劣化した。計画No.2の場合、CTVのD95は、90.52.72Gyであり、OARに照射された最大線量は、最も不良な場合で74.32.23Gyであった。次いで、図21g~iは、計画No.3で作成決定された線量分布を示している。図21gおよびhは0°および45°の各ビームからの線量寄与を示す。なお、計画No.1の場合と同様、315°のビームからの線量寄与については図示を省略している。図21iは、計画No.3で計画された3つすべての照射野からの線量寄与を重ね合わせた合成線量分布を示している。計画No.3の照射計画は、図21に示すように、定性的には計画No.2と同様である。つまり、両計画とも、OARを保全するため、OARの直前で止まるように配置されるペンシルビームの末端部の線量減衰の代わりに、OARの横に配置されるペンシルビームの側面の線量減衰を使用している。なお、計画No.3では、線量勾配抑制項が効いているため、高線量領域は各ビームのz方向にやや拡大され、z方向に広く均一な線量分布が形成されている。計画No.3についても、11の異なる実効密度となるよう、計画No.1,2と同様の摂動で線量分布を再計算した。図22cにこれらの線量分布に関するDVHを示す。図22cに示すように、飛程の不確実性に対する照射計画の感応度は、計画No.3でかなり低減された。また、計画No.2とは対照的に、飛程の不確実性が存在してもCTV内に均一な線量分布が形成されていた。計画No.3の場合、CTV内のD95は、93.82.81Gyであり、OARに照射された最大線量は最も不良な場合で67.22.02Gyであった。2セットアップの不確実性について前記第四項を用いたロバストな照射計画作成方法におけるセットアップの不確実性に対する感応度を検討するため、表1中の計画No.1、4、5を比較した。計画No.1は、前記したように不確実性を全く考慮しない従来の照射計画である。計画No.4は、前記したiiのみを考慮してセットアップの不確実性に対処する照射計画である。計画No.5は、前記したiおよびiiを考慮してセットアップの不確実性に対処する照射計画である。図23a~cは、図21a~cと同様に計画No.1で作成決定された線量分布を示している。セットアップの不確実性の分析のため、計画No.1を幾何学的な摂動とともに再計算した。なお、セットアップの不確実性は、3つの照射野の各々の線量分布に生じ得る。そのため、x方向およびy方向の両方向に±3mm意図的に移動させて計算した。つまり、計画通りのセットアップ位置をこれに含めると、各照射野について9つの線量分布となる。従って、各照射野について合成線量分布の取り得る組み合わせは93=729通りとなるので、これらについて再計算を行った。図24aに729通りの線量分布のDVHを示す。ペンシルビームの側面の線量減衰は、末端部の線量減衰ほど急峻ではないため、セットアップの不確実性に対する従来の照射計画の感応度は、飛程の不確実性に対する感応度ほど重大ではない。しかし、セットアップの誤差はまだ、計画No.1における線量分布の劣化につながる可能性がある。図24aに示されているように、計画No.1の場合、CTV内のD95は、85.22.56Gyであり、OARに照射された最大線量は、最も不良な場合で69.02.07Gyに達した。次いで、図23d~fは、計画No.4で作成決定された線量分布を示している。図23dおよびeは、0°および45°の各ビームからの線量寄与を示す。なお、従前と同様の理由で、315°のビームからの線量寄与については図示を省略している。図23fは、計画No.4で計画された3つすべての照射野からの線量寄与を重ね合わせた合成線量分布を示している。FTVでセットアップの誤差によるCTVの周辺領域内のアンダードーズに対処するため、高線量領域が横断的に拡大した。計画No.4では、セットアップの不確実性に対して高いリスク指数PjSを有するペンシルビームにペナルティーを科した。そのため、OARの横に近接して配置されるペンシルビームは、線量最適化によってペナルティーが科せられている。このため、計画No.4では、ブラッグピークの末端部の線量減衰が、CTVとOARの間の線量勾配を生成するために好適に用いられている。しかし、各ビームによって与えられる線量分布は、図23dおよびeに示されるように横断的に不均一である。そこで、セットアップの誤差の取り得る729通りの組み合わせで計画No.4の線量分布を再計算した。これらの線量分布のDVHを図24bに示す。セットアップの誤差が存在するため、OARに高い線量を照射するリスクのあるペンシルビームは線量最適化においてペナルティーが科せられている。そのため、OARに照射される最大の線量は、最も不良な場合でも53.61.61Gyに抑制することができた。しかしながら、CTV内のD95は、各ビームにより与えられる不均一な線量分布が原因で83.92.52Gyに減少した。次いで、図23g~iは、計画No.5で作成決定された線量分布を示している。図23gおよびhは、0°と45°の各ビームからの線量寄与を示す。なお、前述と同様の理由で、315°のビームからの線量寄与については図示を省略している。図23iは、計画No.5で計画された3つすべての照射野を重ね合わせた合成線量分布を示している。図23iに示す照射計画では、ブラッグピークの末端部の線量減衰が、CTVとOAR間の線量勾配を生成するために好適に用いられており、最適化において、OARの横に近接して配置されるペンシルビームに低い重みが適用されている。計画No.5のこれらの特徴は、計画No.4の特徴と同様である。しかしながら、計画No.4と対照的に、各ビームにより与えられた線量分布は、図23gおよびhに示すように、横断的に均一であった。前述と同様に、計画No.5をセットアップの誤差の取り得る729通りの組み合わせについて再計算した。これらの線量分布のDVHを図24cに示す。図24cに示すように、セットアップの不確実性に対する線量計画の感応度は、計画No.5でも低減していた。図24cに示されているように、計画No.5の場合、CTV内のD95は、90.12.70Gyであり、OARに照射された最大線量は、最も不良な場合で55.51.67Gyであった。3飛程およびセットアップの不確実性について前記1に記載したように、照射計画において飛程の不確実性を考慮する場合は、ペンシルビームの側面の線量減衰がCTVとOARの間の線量勾配の生成に好適に用いられた。また、各照射野の線量分布はz方向に平坦化された。しかしながら、前記2に記載したように、線量計画においてセットアップの不確実性を考慮する場合は、ブラッグピークの末端部の線量減衰がCTVとOARの間の線量勾配の生成に好適に用いられ、各照射野の線量分布は横断的に平坦化された。これらの特徴はまさに相反する。従って、飛程とセットアップの両者の不確実性を線量計画において同時に考慮する場合は、これら2つの相反する目的の妥協点を見出すことになる。この点を検討するため、計画No.3、計画No.6、計画No.7の3つの照射計画を比較した。これらの計画において、飛程の不確実性の最適化パラメータQz、QpbRおよびδRをそれぞれ1×102、1×104、5に設定した。また、セットアップの不確実性のペナルティー係数Qx、QyおよびQpbSは、その概要を表1に示すように設定した。すなわち、各々、計画No.3では0、0、0とし、計画No.6では8.0、8.0、5×103、計画7で32、32、2×104とした。これらの計画における0°および45°の各ビームからの線量寄与を図25に示す。計画No.3においては、OARの直前に配置されるブラッグピークに高い重みを付与することを回避している。そのため、図25aおよびdに示すように、OARの横に配置されるペンシルビームの線量分布が高くなり、谷が形成されたようになっている。計画No.6および計画No.7においては、セットアップの不確実性のペナルティー係数を増加させている。そのため、図25b、eおよび図25c、fに示されるように、谷は横の線量の形状が平坦化され塗りつぶされている。これらのシミュレーションの結果は、最も不良な場合の最適化アルゴリズムを適用して、脊髄の周囲を模擬した馬蹄形の標的体積についての線量分布を最適化したPflugfelderDWilkensJJandOelfkeUWorstcaseoptimizationamethodtoaccountforuncertaintiesintheoptimizationofintensitymodulatedprotontherapyPhysMedBiol53p16891700に記載された結果と同様である。4計算時間3.0GBのRAMを装着したDellPrecision690workstationを用いた計画No.1~7の7つの照射計画の計算時間は表1のTに示したとおりである。表1に示したように、従来の照射計画である計画No.1の計算には4分を要した。一方、飛程の不確実性のみを考慮した計画No.3の計算には9分を要した。また、セットアップの不確実性のみを考慮した計画No.5の計算には8分を要した。そして、飛程とセットアップの両方の不確実性を同時に考慮した計画No.6および計画No.7の計算には11分を要した。表1には、各照射計画で照射されたペンシルビームの総数Npbも示している。計画No.2で要した計算時間と計画No.1で要した計算時間の比は、これらの計画におけるNpbの比と同じ程度である。この関係は、計画No.4と計画No.1の間にもみられる。これらの事実は、高いリスク指標を有するペンシルビームにペナルティーを科す項、すなわち式18の最後の第四項に要する計算時間はほとんど無視できることを意味する。他方、計算時間が、計画No.2に比べて計画No.3で長くなり、計画No.4に比べて計画No.5で長くなるのは、目的関数に線量勾配抑制項第三項を導入したためである。以上に説明したように、第三実施例では、本発明に係る照射計画作成方法を、水に等価でかつ均一であると想定したRTOGベンチマークファントム様形状に適用した。その結果、飛程の誤差およびセットアップの誤差に対して照射計画の感応度を低減したロバストな照射計画を得ることができた。[結論]第一実施例から第三実施例で説明したように、従来の照射計画は飛程の誤差やセットアップの誤差に対して非常に敏感で影響を受け易いものであった。しかし、本発明に係る照射計画作成方法では、予め入力された照射方向から照射対象に与えられる線量分布の勾配を抑制するという第3条件第三項、または当該第3条件第三項に加えて、さらに、予め設定される範囲内で誤差が生じた場合に、非照射部位に与える線量が高くなるスポットを用いないという第4条件第四項を用いた逐次近似繰り返し演算で照射計画を作成決定することによって、飛程の誤差やセットアップの誤差が生じてもその影響を著しく受け難くすることができた。例えば、飛程の誤差が生じた場合、従来の照射計画におけるCTV内の95線量は54.8であったが、本発明の方法では93.8に増加した。一方で、従来の照射計画では、OARに照射された最大線量は最も不良な場合103.2になったが、本発明の方法では67.2に低減した。他方、セットアップの誤差が生じた場合、従来の照射計画におけるCTV内の95線量は85.2であったが、本発明の方法では90.1に増加した。一方で、従来の照射計画では、OARに照射された最大線量は最も不良な場合69.0になったが、本発明の方法では55.5に低減した。本発明に係る照射計画作成方法、装置およびプログラムは、飛程の誤差およびセットアップの誤差の影響を受け難いロバストな照射計画を短時間で得ることができた。"