(57)【特許請求の範囲】
【請求項１】
入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析部と、
前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析部により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定部と、
前記入力時系列信号と、前記周期変動分析部により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定部により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成部と、をそなえた、
ことを特徴とする、信号処理装置。
【請求項２】
前記予測精度推定部が、
前記モデル作成部により出力される、前記入力時系列信号の前記所定時間後の値の算出過程に関する情報を用いて、前記モデル設計情報を推定する、
ことを特徴とする、請求項１記載の信号処理装置。
【請求項３】
前記周期変動分析部が、
前記入力時系列信号の一部と前記入力時系列信号との相関解析を行なう相関解析部と、
前記相関解析部での解析結果を用いて、該周期変動及び前記周期変動の信頼度を推定する周期推定部と、
前記周期推定部により推定した前記周期変動及び前記周期変動の信頼度を用いて、前記入力時系列信号の前記所定時間後までの周期を予測する周期予測部と、
前記相関解析部での前記相関解析の結果と、前記周期予測部により予測した前記周期とを用いて、前記同位相時刻情報と前記同位相時刻情報の信頼度とを推定する同位相時刻推定部と、をそなえた、
ことを特徴とする、請求項１または２に記載の信号処理装置。
【請求項４】
前記モデル作成部が、構造や特性の異なる複数のサブモデルを用いて、前記予測モデルを作成する、
ことを特徴とする、請求項１〜３のいずれか１項に記載の信号処理装置。
【請求項５】
前記サブモデルが、少なくともＳＡＲＩＭＡ（Ｓｅａｓｏｎａｌ  Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル，季節調整指数平滑法モデル，非線型ＡＲＩＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル及びニューラルネットワーク，ファジィのソフトコンピューティングモデルのいずれかである、
ことを特徴とする、請求項４記載の信号処理装置。
【請求項６】
前記モデル作成部が、
前記予測精度推定部により推定した前記入力時系列信号の前記所定時間後の値の予測の精度を用いて、前記サブモデルを変更する、
ことを特徴とする、請求項４または５に記載の信号処理装置。
【請求項７】
入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出し、
前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定し、
前記入力時系列信号と、前記算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記推定した該モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力する、
ことを特徴とする、信号処理方法。
【請求項８】
患者の患部に放射線を照射する放射線治療装置であって、
前記放射線を発生する放射線発生部と、
前記放射線発生部から発生した前記放射線の照射範囲を所望の形状に形成するコリメータ部と、
前記患部の位置を測定する測定部と、
前記測定部により測定された前記患部の位置に関する入力時系列信号を用いて、前記患部の所定時間後の位置の予測と前記予測位置の精度情報とを推定する請求項１〜６のいずれか１項に記載の信号処理装置と、
前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線の照射位置と照射範囲とを算出し、前記算出結果を基に前記コリメータ部を駆動制御する駆動制御部と、をそなえた、
ことを特徴とする、放射線治療装置。
【請求項９】
前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線発生部での前記放射線の発生タイミングを制御するタイミング制御部をそなえた、
ことを特徴とする、請求項８記載の放射線治療装置。
【請求項１０】
患者の患部に放射線を照射する放射線治療装置であって、
前記放射線を発生する放射線発生部と、
前記患部の位置を測定する測定部と、
前記測定部により測定された前記患部の位置に関する入力時系列信号を用いて、前記患部の所定時間後の位置の予測と前記予測位置の精度情報とを推定する請求項１〜６のいずれか１項に記載の信号処理装置と、
前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線発生部での前記放射線の照射タイミングを制御する制御部と、をそなえた、
ことを特徴とする、放射線治療装置。
【請求項１１】
入力時系列信号の所定時間後の値を予測する信号処理装置において、前記予測機能をコンピュータに実現させるための信号処理プログラムであって、
前記信号処理プログラムが、
入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析機能と、
前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定機能と、
前記入力時系列信号と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定機能により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成機能と、を前記コンピュータに実現させる、
ことを特徴とする、信号処理プログラム。
【請求項１２】
入力時系列信号の所定時間後の値を予測する信号処理装置において、前記予測機能をコンピュータに実現させるための信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記信号処理プログラムが、
入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析機能と、
前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定機能と、
前記入力時系列信号と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定機能により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成機能と、を前記コンピュータに実現させる、
ことを特徴とする、信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理装置、信号処理方法、信号処理プログラム及び信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体並びに上記信号処理装置をそなえた放射線治療装置に関する。放射線治療装置には、例えば、患者の患部の位置に正確に放射線を照射する放射線治療装置が含まれる。
【背景技術】
【０００２】
患者の患部（例えば、腫瘍など）に放射線（例えば、陽子線や重粒子線等のイオンビームまたはＸ線など）を照射して治療を施す放射線治療装置が知られている。
上記の放射線治療装置により患者の患部に放射線を照射する場合、例えば、患者の呼吸などにより患部の位置や形状が変動することがある。
これに鑑み、患者の正常組織への影響を抑制しつつ患部組織へ放射線を照射すべく、患部組織が放射線の照射位置に移動した場合にだけ、断続的に放射線を照射する方法がある。【０００３】
しかし、このような方法では、治療時間が長引くため、患者に強いる負担が大きい。また、一般的に、照射する放射線のエネルギー量の総和が等しい場合、照射時間が短時間であるほど治療効果が高くなるため、上記方法では、生物学的効果（治療効果）が低くなるおそれがある。
そこで、患部位置の動きを予測することにより、放射線の照射位置を患部位置に追随させて照射する方法が知られている。【０００４】
なお、放射線治療装置に関する既存の技術として、透視画像取得用の照射タイミングおよび治療用放射線の起動タイミングを制御することにより、微妙な濃淡差が確保された透視画像を取得する方法や、時系列データに基づいて治療放射線照射時点での最も確からしい患部位置を予測することにより、患部に精度良く安全に放射線を照射する方法、また、取得された照射対象の透視画像を特定の評価因子について基準画像と照合することにより、照射対象に対して信頼性の高い照射制御を行なう方法などが知られている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００６−５１１９９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上述した従来方法では、例えば、患部位置の観測信号を入力として得られる入力時系列信号に自己回帰モデルを単純適用して予測を行なう。このような予測は、例えば、極めて短時間（例えば、０．１５秒程度）後の患部位置を予測する際に用いられる。
また、上記方法では、放射線の照射位置を患部位置に追随させるべく、例えば、放射線の照射範囲を形成するコリメータの通過窓の形状及び大きさなどを制御する。【０００７】
しかしながら、コリメータの制御に要する時間が大きくなると、自己回帰モデルで想定している予測時間内にコリメータを所望の形状及び大きさなどに制御することができず、結果、患部位置の動きに追随した照射を実現できない場合がある。
本発明は、以上のような課題に鑑み創案されたもので、患部位置の入力時系列信号を正確且つより中長期に予測することを目的とする。【０００８】
なお、前記目的に限らず、後述する実施形態に示す各構成により導かれる作用効果であって、従来の技術によっては得られない作用効果を奏することも他の目的の一つとして位置付けることができる。
【課題を解決するための手段】
【０００９】
例えば、以下の手段を用いる。
（１）入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析部と、前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析部により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定部と、前記入力時系列信号と、前記周期変動分析部により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定部により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成部と、をそなえた信号処理装置を用いることができる。【００１０】
（２）ここで、前記予測精度推定部が、前記モデル作成部により出力される、前記入力時系列信号の前記所定時間後の値の算出過程に関する情報を用いて、前記モデル設計情報を推定するようにしてもよい。
（３）また、前記周期変動分析部が、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行なう相関解析部と、前記相関解析部での解析結果を用いて、該周期変動及び前記周期変動の信頼度を推定する周期推定部と、前記周期推定部により推定した前記周期変動及び前記周期変動の信頼度を用いて、前記入力時系列信号の前記所定時間後までの周期を予測する周期予測部と、前記相関解析部での前記相関解析の結果と、前記周期予測部により予測した前記周期とを用いて、前記同位相時刻情報と前記同位相時刻情報の信頼度とを推定する同位相時刻推定部とをそなえていてもよい。【００１１】
（４）さらに、前記モデル作成部が、構造や特性の異なる複数のサブモデルを用いて、前記予測モデルを作成してもよい。
（５）また、前記サブモデルが、少なくともＳＡＲＩＭＡ（Ｓｅａｓｏｎａｌ  Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル，季節調整指数平滑法モデル，非線型ＡＲＩＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル及びニューラルネットワーク，ファジィなどのソフトコンピューティングモデルのいずれかであってもよい。【００１２】
（６）さらに、前記モデル作成部が、前記予測精度推定部により推定した前記入力時系列信号の前記所定時間後の値の予測の精度を用いて、前記サブモデルを変更するようにしてもよい。
（７）また、入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出し、前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定し、前記入力時系列信号と、前記算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度情報と、前記推定した該モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力する信号処理方法を用いることができる。【００１３】
（８）さらに、患者の患部に放射線を照射する放射線治療装置であって、前記放射線を発生する放射線発生部と、前記放射線発生部から発生した前記放射線の照射範囲を所望の形状に形成するコリメータ部と、前記患部の位置を測定する測定部と、前記測定部により測定された前記患部の位置に関する入力時系列信号を用いて、前記患部の所定時間後の位置の予測と前記予測位置の精度情報とを推定する上記の信号処理装置と、前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線の照射位置と照射範囲とを算出し、前記算出結果を基に前記コリメータ部を駆動制御する駆動制御部と、をそなえた放射線治療装置を用いることができる。【００１４】
（９）また、上記の放射線治療装置は、前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線発生部での前記放射線の発生タイミングを制御するタイミング制御部をそなえていてもよい。
（１０）さらに、患者の患部に放射線を照射する放射線治療装置であって、前記放射線を発生する放射線発生部と、前記患部の位置を測定する測定部と、前記測定部により測定された前記患部の位置に関する入力時系列信号を用いて、前記患部の所定時間後の位置の予測と前記予測位置の精度情報とを推定する上記の信号処理装置と、前記信号処理装置により予測した前記患部の所定時間後の位置と前記精度情報とを用いて、前記放射線発生部での前記放射線の照射タイミングを制御する制御部と、をそなえた放射線治療装置を用いることができる。【００１５】
（１１）また、入力時系列信号の所定時間後の値を予測する信号処理装置において、前記予測機能をコンピュータに実現させるための信号処理プログラムであって、前記信号処理プログラムが、入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析機能と、前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定機能と、前記入力時系列信号と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定機能により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成機能と、を前記コンピュータに実現させる、信号処理プログラムを用いることができる。【００１６】
（１２）さらに、入力時系列信号の所定時間後の値を予測する信号処理装置において、前記予測機能をコンピュータに実現させるための信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体であって、前記信号処理プログラムが、入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析機能と、前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定機能と、前記入力時系列信号と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定機能により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成機能と、を前記コンピュータに実現させる、信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体を用いることができる。
【発明の効果】
【００１７】
上記本発明によれば、入力時系列信号の将来の情報を正確かつ中長期に予測することが可能となる。また、放射線治療においては、患部への正確かつ連続的な照射が可能になる。
【図面の簡単な説明】
【００１８】
【図１】放射線治療装置の概要を示す図である。
【図２】図１に示す放射線治療装置の機能ブロック図である。
【図３】（Ａ）はＭＬＣ２１０の概観を示す図であり、（Ｂ）はＭＬＣ２１０の構成の一例を示す図である。
【図４】ＯＢＩ２２０の構成の一例を示す図である。
【図５】制御部１の構成の一例を示す図である。
【図６】（Ａ）〜（Ｃ）は周期変動分析部１２の動作例を示す図である。
【図７】（Ａ）〜（Ｄ）は相関解析部１２１の動作例を示す図である。
【図８】（Ａ）〜（Ｃ）は周期推定部１２３の動作例を示す図である。
【図９】周期予測部１２４の動作例を示す図である。
【図１０】（Ａ）〜（Ｄ）は同位相時刻推定部１２５の動作例を示す図である。
【図１１】（Ａ）及び（Ｂ）はモデル作成部１１の動作例を示す図である。
【図１２】（Ａ）及び（Ｂ）は予測精度推定部１６の動作例を示す図である。
【図１３】患部２６０の所定の観測位置の３次元プロットを示す図である。
【図１４】図１３に示す３次元プロットの各軸成分を示す図である。
【図１５】図１４に示す各軸成分についての各サブモデルでの対象時系列の実測値とその予測値を示す図である。
【図１６】（Ａ）は図１４に示す各軸成分の実測値及び予測値を示す図であり、（Ｂ）及び（Ｃ）を予測精度推定部１６による予測精度推定情報を示す図である。
【図１７】図１４に示す各軸成分についての各信号成分の実測値及び従来方法による予測値を示す図である。
【図１８】本例の予測方法と従来の予測方法との比較結果を示す図である。
【図１９】本例の予測方法の一例を説明するフローチャートである。
【図２０】本例の予測方法の一例を説明するフローチャートである。
【図２１】本例の予測方法の一例を説明するフローチャートである。
【図２２】本例の予測方法の一例を説明するフローチャートである。
【図２３】本例の予測方法の一例を説明するフローチャートである。
【図２４】本例の予測方法の一例を説明するフローチャートである。
【図２５】本例の予測方法の一例を説明するフローチャートである。
【図２６】本例の予測方法の一例を説明するフローチャートである。
【図２７】本例の予測方法の一例を説明するフローチャートである。
【図２８】本例の予測方法の一例を説明するフローチャートである。
【図２９】周期固定の季節成分を説明する図である。
【図３０】周期が変動する季節成分を説明する図である。
【図３１】本例の予測方法と従来の予測方法との比較結果の一例を示す図である。
【発明を実施するための形態】
【００１９】
以下、図面を参照して実施の形態を説明する。ただし、以下に示す実施形態は、あくまでも例示に過ぎず、以下に示す実施形態で明示しない種々の変形や技術の適用を排除する意図はない。即ち、本実施形態は、その趣旨を逸脱しない範囲で種々変形（各実施形態を組み合わせる等）して実施することができる。
【００２０】
〔１〕−実施形態
（１．１）装置概要
図１に放射線治療装置の概要を示す。この図１に示す放射線治療装置は、例示的に、制御部１と、クライストロン１００と、冷却装置１１０と、導波管１２０と、真空ポンプ１３０と、電子銃１４０と、加速管１５０と、偏向電磁石１６０と、ターゲット１７０と、平坦化フィルタ（散乱フォイル）１８０と、をそなえる。また、図１に示す放射線治療装置は、例示的に、モニタ線量計１９０と、照射野形成部〔ジョー、（ＪＡＷ）〕２００と、ＭＬＣ（Ｍｕｌｔｉ  Ｌｅａｆ  Ｃｏｌｌｉｍａｔｏｒ）２１０と、ＯＢＩ（Ｏｎ−Ｂｏａｒｄ  Ｉｍａｇｅｒ）２２０と、ベッド２３０と、ＥＰＩＤ（Ｅｌｅｃｔｒｏｎｉｃ  Ｐｏｒｔａｌ  Ｉｍａｇｉｎｇ  Ｄｅｖｉｃｅ）２４０と、パルスモジュレータ２５０と、をそなえる。
【００２１】
本例では、例えば、制御部１が、ＯＢＩ２２０により測定した患部２６０の位置情報（位置，形状及び大きさなど）から、所定時間後の患部２６０の位置情報を予測する。そして、制御部１は、前記予測結果に基づき、例えば、ＭＬＣ２１０を駆動制御することにより、所定時間後の患部２６０の位置，形状及び大きさなどに応じた放射線照射を実現する。
【００２２】
このため、図１に例示した放射線治療装置は、図２に例示する機能ブロックを有する。
ここで、パルスモジュレータ２５０は、パルス電圧を発生する電圧発生源である。本例のパルスモジュレータ２５０は、例えば、トリガーパルス発生器（図示省略）からのトリガーパルスの入力を契機として、所定のパルス電圧を発生する。
クライストロン１００は、電子ビームと高周波との相互作用を利用して高周波を増幅して大電力のマイクロ波を出力する電子管である。本例のクライストロン１００は、例えば、パルスモジュレータ２５０から受けたパルス電圧により電子ビームを発生させ、この電子ビームを速度変調して間隙から出力することにより、加速管１５０で電子を加速するためのエネルギーとなるマイクロ波パルス電力を出力する。
【００２３】
冷却装置１１０は、クライストロン１００を冷却する機能を具備する。冷却の方法としては、例えば、クライストロン１００の本体部を覆う外装部と前記本体部との間に空洞を設け、この空洞に冷却水を流入して冷却する方法などがある。
導波管１２０は、クライストロン１００により出力されたマイクロ波パルス電力を加速管１５０へ伝達する。導波管１２０には、例えば、中空導波管を用いることができるが、そのほか、テフロン（登録商標）などの誘電体を用いた同軸ケーブルを用いることもできる。
【００２４】
真空ポンプ１３０は、容器内から気体を排出して真空状態を形成するためのポンプである。本例では、例えば、電子銃１４０の電子線放出部の近傍部分，加速管１５０の内部空間及び偏向電磁石１６０の近傍部分に真空状態を形成する。これらの各部分を真空状態に保つことで、クライストロン１００により出力されたマイクロ波パルス電力の損失や、電子銃１４０から発射された電子線のエネルギー損失を低減することができる。なお、真空ポンプ１３０には、例えば、加速管１５０の内部などへの侵食度合いが少ないドライポンプを用いることができるが、そのほか、油を用いるウェットポンプを用いることもできる。
【００２５】
電子銃１４０は、自由電子を所定の金属材料などから取り出して電子線を放出する。例えば、固体中の電子を高熱や高電界により空間に放出させ、放出した電子に電界を加えて加速するとともに、電子レンズにより電子線をビーム状に集束させて照射する。電子銃１４０には、例えば、量子力学的トンネル効果により空間に電子を取り出す電界放射型電子銃を用いることができるが、そのほか、熱エネルギーにより電子放出を行なう熱電子放出型電子銃などを用いることもできる。
【００２６】
加速管１５０は、クライストロン１００にて出力され導波管１２０により伝達されるマイクロ波パルス電力を用いて、電子銃１４０から放出された電子線のエネルギーを増幅する。電子線は、マイクロ波パルス電力により加速されて加速管１５０から放出される。
偏向電磁石１６０は、所定の方向に磁界を発生させることにより、加速管１５０から放出された電子線の軌道を修正する。本例の偏向電磁石１６０は、例えば、ターゲット１７０へ電子線を照射させるべく、例示的に、３個の偏向電磁石１６０を用いて電子線の軌道を修正する。なお、偏向電磁石１６０の数は、図１に示す例に限定されない。
【００２７】
ターゲット１７０は、電子線を照射されることにより、放射線を照射する材料である。
ターゲット１７０には、例えば、タングステン（Ｔｇ）を用いることができる。なお、加速された電子線の照射により放射線を発生することのできる材料であれば、その他の材料を用いてもよい。
即ち、パルスモジュレータ２５０，クライストロン１００，冷却装置１１０，導波管１２０，真空ポンプ１３０，電子銃１４０，加速管１５０，偏向電磁石１６０及びターゲット１７０は、放射線を発生する放射線発生部の一例として機能する。
【００２８】
平坦化フィルタ１８０は、ターゲット１７０から照射された放射線の線量分布を減衰して、平坦化させるフィルタである。本例の平坦化フィルタ１８０は、例えば、ターゲット１７０から照射された放射線の線量レベルを波長に応じて減衰し、放射線の線量レベルを所定のレベル（例えば、人体への放射が可能な線量レベル）に制御する。
モニタ線量計１９０は、放射線の線量レベルを計測する。本例のモニタ線量計１９０により、例えば、平坦化フィルタ１８０を通過後の放射線の線量レベルが所定のレベルにまで減衰されているかどうかを計測、確認することができる。もし、放射線の線量レベルが所定のレベルまで減衰されていない場合は、例えば、クライストロン１００でのマイクロ波パルス電力などを制御することなどにより、線量レベルを調節することが可能である。
【００２９】
ＪＡＷ２００は、放射線を遮断する材料からなる、照射野形成部材である。ＪＡＷ２００は、例えば、複数の板状部材で構成され、板状部材の各位置を決定することにより、ターゲット１７０から照射された放射線の一部を通過させるとともに、その他を遮断することにより、放射線の照射野をおおまかに形成する。例えば、放射線が、ＪＡＷ２００の後段に位置するＭＬＣ２１０への入力が可能な範囲に収まるように、放射線の照射野を絞ることができる。
【００３０】
ＭＬＣ（コリメータ部）２１０は、ターゲット１７０から照射される放射線の照射範囲を所望の形状に形成する。本例のＭＬＣ２１０は、例えば、制御部１にて予測された所定時間後の患部２６０の位置と予測精度情報とに応じて、ＪＡＷ２００を通過後の放射線の照射範囲を変化させることができる。
図３（Ａ）及び図３（Ｂ）にＭＬＣ２１０の概観及び構成例を示す。
【００３１】
図３（Ａ）に例示するように、ＭＬＣ２１０は、例えば、放射線を遮断する板状部材（リーフ）の対を複数有して構成される。各リーフが、紙面左右方向にそれぞれ独立して変位することで、協動して、所望の形状に放射線照射野〔図３（Ａ）の中央の空間部分〕３０４を形成する。ＪＡＷ２００を通過した放射線が、ＭＬＣ２１０を紙面垂直方向に通過することにより、放射線照射野３０４の形状に応じた範囲で放射線を患部２６０に照射することができる。
【００３２】
ここで、ＭＬＣ２１０は、図３（Ｂ）に例示するように、複数のリーフ３０３と、駆動部３０２−１，３０２−２と、駆動制御部３０１とを有する。
リーフ３０３は、放射線を遮断可能な板状部材である。本例では、例えば、対向する２個のリーフ３０３が対となり、複数のリーフ対を形成する。各リーフ３０３は、駆動部３０２−１，３０２−２によりそれぞれ独立的に駆動される。なお、リーフ３０３の数は図示する例に限定されない。
【００３３】
駆動部３０２−１，３０２−２は、各リーフ３０３を独立的に変位させる。例えば、駆動部３０２−１，３０２−２は、各リーフ３０３に個別に設けたモータなどにより、各リーフ３０３を駆動して、図中のリーフ移動方向に沿って変位させる。図３（Ｂ）に示す例では、駆動部３０２−１により、紙面右側のリーフ３０３がそれぞれ駆動され、また、駆動部３０２−２により、紙面左側のリーフ３０３がそれぞれ駆動される。また、例えば、駆動部３０２−１，３０２−２により、各リーフ３０３を連動して変位させることもできる。例えば、あるリーフ３０３の変位をトリガとして、その周辺のリーフ３０３を連動して変位させるようにしてもよい。
【００３４】
駆動制御部３０１は、駆動部３０２−１，３０２−２を制御して、リーフ３０３のそれぞれの位置を決定する。各リーフ３０３が、それぞれの所望の位置に独立して位置決め制御されることにより、所望の形状を有する放射線照射野３０４が形成される。
本例の駆動制御部３０１は、例えば、制御部１による予測結果（所定時間後の患部２６０の予測位置及び予測精度情報）に基づいて、駆動部３０２−１，３０２−２を制御し、各リーフ３０３の位置を決定する。これにより、制御部１により予測された所定時間後の患部２６０の位置，形状及び大きさに応じて、所望の形状の放射線照射野３０４を形成することが可能となるので、患部２６０の予測位置，形状及び大きさなどに適した放射線照射を実現することができる。
【００３５】
ＯＢＩ（測定部）２２０は、患部２６０の位置を測定する。例えば、測定対象にＸ線を照射して、測定対象の位置，形状及び大きさなどに関する位置情報を出力することができる。
本例のＯＢＩ２２０は、例えば、患者の患部２６０にＸ線を複数の角度から照射することにより得られる複数の位置情報に基づき、患部２６０の位置，形状及び大きさなどの観測値（実測値）を算出する。ＯＢＩ２２０により実測された患部２６０の位置，形状及び大きさなどに関する位置情報は、制御部１に入力され、後述する予測処理に用いられる。
【００３６】
ここで、図４にＯＢＩ２２０の構成例を示す。この図４に示すように、ＯＢＩ２２０は、例示的に、Ｘ線照射部３０５と、センサーパネル（Ｓｅｎｓｏｒ  Ｐａｎｅｌ）３０６とを有する。
Ｘ線照射部３０５は、ベッド２３０上の患者の患部２６０に対してＸ線を照射する。なお、Ｘ線は、放射線の一種であり、例えば、その波長が１ｐｍ〜１０ｎｍ程度の電磁波である。
【００３７】
センサーパネル３０６は、表面にＸ線を検出するセンサを有し、照射されたＸ線の位置座標を出力する。本例のセンサーパネル３０６は、例えば、患部２６０を通過後の上記Ｘ線をイメージング領域（Ｉｍａｇｉｎｇ  Ａｒｅａ）にて受信し、その受信位置の二次元座標を算出、出力する。
Ｘ線照射部３０５及びセンサーパネル３０６は、図４中の回転軸を中心として、一体的に回転しうるように構成されており、本例のセンサーパネル３０６は、例えば、図４中の符号Ａで示す位置（Ｘ線照射軸と回転軸とが＋４５°の角度をなす位置）で算出した患部２６０の２次元座標Ｐ１（ｘ１，ｙ１）と、図４中の符号Ｂで示す位置（Ｘ線照射軸と回転軸とが−４５°の角度をなす位置）で算出した患部２６０の２次元座標Ｐ２（ｘ２，ｙ２）とに基づき、患部２６０の３次元座標Ｐ（Ｘ，Ｙ，Ｚ）を算出する。
【００３８】
例えば、Ｘ線照射部３０５の照射点から患部２６０までの距離をＳＡＤ（ｓｏｕｒｃｅ−ａｘｉｓ  ｄｉｓｔａｎｃｅ）とし、Ｘ線照射部３０５のＸ線照射点からセンサーパネル３０６の表面までの距離をＳＤＤ（ｓｏｕｒｃｅ−ｄｅｔｅｃｔｏｒ  ｄｉｓｔａｎｃｅ）とすると、患部２６０の位置座標Ｐ（Ｘ，Ｙ，Ｚ）は、次式（１）で表される。
【００３９】
【数１】
0005604306.tif 000002
【００４０】
上記の式（１）により算出された患部２６０の３次元座標Ｐ（Ｘ，Ｙ，Ｚ）は、例えば、患部２６０の現在位置を示す実測値（入力信号）として制御部１へ入力され、制御部１は、この実測値の時系列情報を用いて、所定時間後の患部２６０の３次元座標を予測することができる。
ベッド２３０は、患部２６０の位置が放射線照射に適した位置及び高さになるよう、患者を支える器具である。このため、本例のベッド２３０は、例えば、位置及び高さ調節可能な支持部とマット部とで構成される。
【００４１】
ＥＰＩＤ２４０は、患部２６０の位置、形状及び大きさを可視化する。例えば、可視化された患部２６０の形状などに基づき、施術者により、ベッド２３０などの位置及び高さが決定される。
制御部１は、ＯＢＩ２２０により得られる患部２６０の位置座標Ｐ（Ｘ，Ｙ，Ｚ）の時系列情報を用いて、所定時間後の患部２６０の位置座標の予測と予測精度情報とを推測する。また、予測結果と予測精度情報とを用いて、各種の制御を行なう。例えば、本例の制御部１は、前記予測結果と予測精度情報とを用いて、ＭＬＣ２１０を制御して、所定時間後の患部２６０の位置，予測精度情報，大きさ及び形状などに応じた形に放射線照射野３０４を形成する。これにより、ＭＬＣ２１０の物理的な駆動（位置決め）制御に要する時間に先だって、ＭＬＣ２１０を予測制御できるので、患部２６０の位置に追従して放射線照射を行なうことが可能となる。
【００４２】
即ち、本例の制御部１は、ＯＢＩ２２０により測定した患部２６０の位置の時系列情報を用いて、患部２６０の所定時間後の位置と予測精度情報とを予測する信号処理装置の一例として機能する。また、制御部１は、予測した患部２６０の所定時間後の位置と予測位置の精度情報とを用いて、放射線の照射位置と照射範囲とを算出し、前記算出結果を基にＭＬＣ２１０を駆動制御する駆動制御部の一例として機能する。
【００４３】
また、制御部１は、ＭＬＣ２１０の上記駆動制御に加え、予測した患部２６０の所定時間後の位置と予測精度情報とを用いて、放射線の照射タイミングを制御するようにしてもよい。例えば、制御部１は、電子銃１４０での電子線の発射タイミングを制御することにより、放射線の発生タイミングを制御するようにしてもよいし、クライストロン１００により出力するマイクロ波電力量を制御することで、ターゲット１７０における放射線の発生タイミングを制御するようにしてもよい。
【００４４】
即ち、本例の制御部１は、上記方法により予測した患部２６０の所定時間後の位置と予測位置の精度情報とを用いて、放射線の発生タイミングを制御するタイミング制御部の一例として機能する。
なお、制御部１は、上記駆動制御及びタイミング制御のいずれか一方のみを行なうようにしてもよい。例えば、制御部１がタイミング制御のみを行なう場合は、放射線の照射位置を固定し、上記の予測結果に基づいて、患部２６０がその照射位置に移動するときを狙って、放射線を照射するので、上述した例と同様の効果を得られる上、ＭＬＣ２１０を備えなくてもよいので、装置を単純化することができ、さらに装置コストを削減することが可能となる。
【００４５】
以下、制御部１の構成例及び予測方法について説明する。
（１．２）制御部１
制御部１は、主要成分が周期時変の季節ダイナミクス（周期的に繰り返される時系列信号）である任意の時系列信号を入力として、その所定時刻経過後の予測値と、その予測精度に関する情報（予測精度情報）とを出力する。
【００４６】
図５に例示する制御部１は、モデル作成部１１と、周期変動分析部１２と、変量推定部１３と、キャッシュメモリ１４と、キャッシュメモリ１５と、予測精度推定部１６とを有する。
ここで、制御部１への入力は、観測開始時刻から現在時刻ｔｃまでの患部２６０の位置の変化に対応する時系列信号
【００４７】
0005604306.tif 000003
である。ここで、
【００４８】
【数２】
0005604306.tif 000004
であり、また、
0005604306.tif 000005
は時刻τにおけるｊ（ｊ＝１〜Ｄｉｍ）次元方向の患部の位置であり、Ｄｉｍは位置情報の次元数である。
一方、制御部１の出力は、任意の予測期間をＬとした場合の予測対象時刻（ｔｃ＋Ｌ）における患部２６０の位置
0005604306.tif 000006
の予測値
【００４９】
【数３】
0005604306.tif 000007
と、その予測値の精度情報である誤差量の推定値
【００５０】
【数４】
0005604306.tif 000008
と、単位時間当たりの位置変動量の推定値
【００５１】
【数５】
0005604306.tif 000009
とである。
【００５２】
（１．３）モデル作成部１１
モデル作成部１１は、制御部１への入力である上記の時系列信号と、周期変動分析部１２の出力の一部である、予測対象時刻の位相と同位相の同位相時刻情報と、その同位相時刻情報の信頼度と、予測精度推定部１６の出力の一部であるモデル設計情報とを用いて、予測モデルを作成し、前記予測モデルを用いて、入力時系列信号の所定時刻経過後の予測値と、その予測値を算出する過程の詳細情報とを出力する。なお、モデル作成部１１は、フィルタ及びモデルの種類やパラメータ値が異なる、周期時変の影響を考慮した複数のサブモデルを作成し、各サブモデルを合成することにより前記予測モデルを作成してもよい。また、モデル作成部１１は、予め作成された複数のサブモデルを用いて、前記予測モデルを作成してもよい。
【００５３】
図５及び図１１（Ａ）に例示するモデル作成部１１は、サブモデル１１１−１〜１１１−Ｓ（Ｓはサブモデルの総数）と合成処理部１１２とを有する。
本例のモデル作成部１１は、制御部１への入力時系列信号
0005604306.tif 000010
と、周期変動分析部１２の出力の一部である、予測対象時刻（ｔｃ＋Ｌ）の信号
0005604306.tif 000011
の位相と同位相の入力時系列信号における時刻
【００５４】
【数６】
0005604306.tif 000012
【００５５】
【数７】
0005604306.tif 000013
と、その信頼度
【００５６】
【数８】
0005604306.tif 000014
【００５７】
【数９】
0005604306.tif 000015
と、予測精度推定部１６の出力の一部であるモデル設計情報
【００５８】
【数１０】
0005604306.tif 000016
とを入力として、予測対象時刻（ｔｃ＋Ｌ）の予測値
【００５９】
【数１１】
0005604306.tif 000017
【００６０】
【数１２】
0005604306.tif 000018
と、その予測値算出の過程の詳細情報
【００６１】
【数１３】
0005604306.tif 000019
【００６２】
【数１４】
0005604306.tif 000020
【００６３】
【数１５】
0005604306.tif 000021
とを出力する。ここで、
0005604306.tif 000022
は、サブモデルｓ（ｓ＝１，２，…，Ｓ）のｊ次元方向の出力値である。さらに、関数
0005604306.tif 000023
は、合成処理部１１２による加算・乗算処理に対応し、モデル設計情報
0005604306.tif 000024
を用いて決定され、各サブモデル１１１−１〜１１１−Ｓの出力の合成として、予測対象時刻（ｔｃ＋Ｌ）の予測値を算出する。ここで、
0005604306.tif 000025
は、サブモデルｓのモデル設計情報であり、
0005604306.tif 000026
は、サブモデルｓ出力
0005604306.tif 000027
の出力算出過程の情報（算出過程情報）と関数
0005604306.tif 000028
の設計情報の一部である。さらに、その要素である
0005604306.tif 000029
は、モデルパラメータの情報や出力算出過程や関数設計の具体的な情報である。
即ち、本例のモデル作成部１１は、入力時系列信号と、周期変動分析部１２により算出される同位相時刻情報及び前記同位相時刻情報の信頼度と、予測精度推定部１６により推定されるモデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成部の一例として機能する。
【００６４】
また、モデル作成部１１は、制御部１への入力である主要成分が周期時変の季節ダイナミクスである任意の時系列信号と、周期変動分析部１２の出力の一部である予測対象時刻の位相と同位相の同位相時刻情報とその信頼度情報と、予測精度推定部１６の出力の一部であるモデル設計情報とを用いて、フィルタ及びモデルの種類やパラメータ値が異なる複数のサブモデルと、それらを合成するモデルとを作成し、入力時系列信号の所定時刻経過後の予測値と、その予測値算出の過程の詳細情報を出力することもできる。
【００６５】
さらに、本例のモデル作成部１１は、予測精度推定部１６により推定される予測精度情報を用いて、サブモデルと、それらを合成するモデルの設計を変更することができる。
【００６６】
（１．３．１）サブモデル１１１−１〜１１１−Ｓ
サブモデル１１１−１〜１１１−Ｓは、モデル設計情報（フィルタ及びモデルの種類、また、そのパラメータの情報）を用いて作成され、入力時系列信号（またはその特定成分）の所定時刻経過後の値を予測するモデルである。サブモデルの数は、フィルタ及びモデルの種類や、そのパラメータの種類の組み合わせ総数により決定される。
【００６７】
サブモデル１１１−１〜１１１−Ｓの設計において、フィルタ処理後の時系列信号の季節ダイナミクスの周期が時間変化（時変）する場合には、周期変動分析部１２からの入力である予測対象時刻の位相と同位相の入力時系列信号における同位相時刻情報とその信頼度情報とを用いて、その周期時変が原因で予測誤差が大きくなる問題を回避する。
各サブモデルで用いるフィルタの種類には、例えば、移動平均フィルタ，差分フィルタ，バンドパスフィルタ及びカルマンフィルタなどのうち少なくともいずれかを用いることができる。
【００６８】
サブモデル１１１−１〜１１１−Ｓには、例えば、ＳＡＲＩＭＡ（Ｓｅａｓｏｎａｌ  Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル，季節調整指数平滑法モデル，非線型ＡＲＩＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ）モデル及びニューラルネットワーク，ファジィなどのソフトコンピューティングモデルなどの少なくともいずれかを用いることができる。
【００６９】
また、モデルの設計には、例えば、ユール・ウォーカー法，バーグ法，最小二乗法及び勾配法の設計アルゴリズムと、任意の設計用パラメータと評価関数と、ＡＩＣ（赤池情報量規準）や最尤法の評価アルゴリズムとを用いることができる。
図１１（Ｂ）に例示するサブモデルｓ（１１１−１〜１１１−Ｓ）は、モデルの種類にＳＡＲＩＭＡモデルを用いて、さらに、予測対象時刻（ｔｃ＋Ｌ）の信号と位相が同位相の入力時系列信号における時刻情報
0005604306.tif 000030
を用いて周期時変の季節ダイナミクスを効果的にモデル化する工夫と、予測対象時刻（ｔｃ＋Ｌ）を直接的に予測する工夫とを実現する例である。ここで、入力時系列信号
0005604306.tif 000031
を入力とした場合の出力は次式のように求まる。
【００７０】
【数１６】
0005604306.tif 000032
ここで、上記の式（２）の導出方法について説明する。
上記の式（２）は、季節成分の周期が変動する時系列データを記述するための時系列モデルの式である。
まず、上記の式（２）の導出方法を理解しやすくするために、いくつかの言葉の定義と、ＡＲ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｏｎ：自己回帰），ＭＡ（Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ：移動平均），ＡＲＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｏｎ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ：自己回帰移動平均），ＡＲＩＭＡ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｏｎ  Ｉｎｔｅｇｒａｔｅｄ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ：自己回帰和分移動平均），ＳＡＲＩＭＡの各モデルの説明とを行なう。
まず、時系列モデルとは、時系列データの振る舞いの規則を記述する関係式のことをいう。また、定常過程とは、確率分布が変化しない確率過程のことをいう。さらに、定常時系列とは、定常過程により生成される時系列データのことをいい、時系列のどの区間においても、平均値と分散が変化しないという特徴を有する。なお、定常時系列に対する、ＡＲ，ＡＲＭＡモデルの設計方法は確立されている。
また、トレンド成分とは、時系列の平均値を、長時間継続して変動させる成分のことをいう。さらに、季節成分とは、時系列の平均値を、周期的に変動させる成分のことをいい、非定常時系列とは、定常時系列以外の時系列のことをいう。この非定常時系列には、例えば、定常時系列に、トレンド成分や季節成分などが混ざった時系列がある。
まず、ＡＲモデルについて説明する。
ＡＲモデルは、定常時系列ｚ（ｔ）（ｔ＝０，１，２，・・・）の各値を、それより過去の時系列データの線形方程式として次式（３）のように表す。
ｚ（ｔ）＝φ１ｚ（ｔ−１）＋φ２ｚ（ｔ−２）＋・・・＋φｐｚ（ｔ−ｐ）＋ｅ（ｔ）  ・・・（３）
ここで、ｐはモデルの次数、φｉ（ｉ＝１，２，・・・，ｐ）はモデル係数、ｅ（ｔ）は時刻ｔにおけるモデル誤差である。
上記の式（３）の左辺から、ｅ（ｔ）以外を右辺へ移項すると、以下の式（４）が得られる。
ｚ（ｔ）−φ１ｚ（ｔ−１）−φ２ｚ（ｔ−２）−・・・−φｐｚ（ｔ−ｐ）＝ｅ（ｔ）  ・・・（４）
さらに、遅延演算子Ｂとして、例えば、次式（５）
Ｂｄｚ（ｔ）＝ｚ（ｔ−ｄ）  ・・・（５）
を満足するＢを用いて、上記の式（４）を表記すると、以下の式（６），（７）が得られる。
φｐｚ（ｔ）＝ｅ（ｔ）  ・・・（６）
φｐ≡（１−φ１Ｂ−φ２Ｂ２−・・・−φｐＢｐ）  ・・・（７）
次に、ＭＡモデルについて説明する。
ＭＡモデルは、定常時系列ｚ（ｔ）の各値を、それまでに発生したノイズの時系列データｅ（ｔ）の線形方程式として次式（８）のように表す。
ｚ（ｔ）＝ｅ（ｔ）−θ１ｅ（ｔ−１）−θ２ｅ（ｔ−２）−・・・−θｑｅ（ｔ−ｑ）  ・・・（８）
ここで、ｑはモデルの次数、θｊ（ｊ＝１，２，・・・，ｑ）はモデル係数である。
上記の式（８）を、遅延演算子Ｂを用いて表記すると、以下の式（９），（１０）が得られる。
ｚ（ｔ）＝θｑｅ（ｔ）  ・・・（９）
θｑ≡（１−θ１Ｂ−θ２Ｂ２−・・・−θｑＢｑ）  ・・・（１０）
さらに、ＡＲＭＡモデルについて説明する。
ＡＲＭＡモデルは、定常時系列ｚ（ｔ）の各値を、それより過去の時系列データと、それまでに発生したノイズの時系列データｅ（ｔ）との線形方程式として次式（１１）のように表す。
ｚ（ｔ）＝φ１ｚ（ｔ−１）＋φ２ｚ（ｔ−２）＋・・・＋φｐｚ（ｔ−ｐ）＋ｅ（ｔ）−θ１ｅ（ｔ−１）−θ２ｅ（ｔ−２）−・・・−θｑｅ（ｔ−ｑ）  ・・・（１１）
ここで、上記の式（１１）の左辺から、ノイズに関する項以外を右辺へ移項すると、次式（１２）が得られる。
ｚ（ｔ）−φ１ｚ（ｔ−１）−φ２ｚ（ｔ−２）−・・・−φｐｚ（ｔ−ｐ）＝ｅ（ｔ）−θ１ｅ（ｔ−１）−θ２ｅ（ｔ−２）−・・・−θｑｅ（ｔ−ｑ）  ・・・（１２）
さらに、上記の式（１２）を、遅延演算子Ｂ，φｐ，θｑを用いて表記すると、次式（１３）が得られる。
φｐｚ（ｔ）＝θｑｅ（ｔ）  ・・・（１３）
次に、ＡＲＩＭＡモデルについて説明する。
ＡＲＩＭＡモデルは、トレンド成分を含む非定常時系列ｙ（ｔ）（ｔ＝１，２，・・・）の各値を記述するためのモデルである。
非定常時系列ｙ（ｔ）に含まれるトレンド成分は、所定の回数ｄ（≧１）回、遅れ１の時系列データとの差分をとることで除去することができる。即ち、このときの残差ｚ（ｔ）は次式（１４）で表される。
ｚ（ｔ）＝（１−Ｂ）ｄｙ（ｔ）  ・・・（１４）
ここで、適切に差分をとった残差ｚ（ｔ）の時系列は定常時系列となるため、ＡＲＭＡモデルによりモデル化することができる。即ち、次式（１５）が成立する。
φｐｚ（ｔ）＝θｑｅ（ｔ）  ・・・（１５）
したがって、ＡＲＩＭＡモデルは、上記の式（１５）に式（１４）を代入することにより得られる次式（１６）で表される。
φｐ（１−Ｂ）ｄｙ（ｔ）＝θｑｅ（ｔ）  ・・・（１６）
次に、ＳＡＲＩＭＡモデルについて説明する。
ＳＡＲＩＭＡモデルは、トレンド成分と季節成分とを含む非定常時系列ｘ（ｔ）（ｔ＝１，２，・・・）の各値を記述するためのモデルである。なお、以下の説明では、季節成分の周期ｃがｃ＞＞１である場合を考える。
非定常時系列ｘ（ｔ）に含まれるトレンド成分は、ＡＲＩＭＡモデルの説明にて示したように、ｄ（≧１）回、遅れ１の時系列データとの差分をとることで除去することができる。このときの残差ｗ（ｔ）は次式（１７）で表される。
ｗ（ｔ）＝（１−Ｂ）ｄｘ（ｔ）  ・・・（１７）
また、トレンド成分が除去された後の非定常時系列ｗ（ｔ）に含まれる季節成分は、次の［１］〜［３］のいずれかの方法により除去することができる。
方法［１］として、例えば、Ｄ（≧０）回、遅れｃの時系列データとの差分をとることにより、このときの残差ｖ（ｔ）は次式（１８）で表される。
ｖ（ｔ）＝（１−Ｂｃ）Ｄｗ（ｔ）  ・・・（１８）
また、方法［２］として、例えば、同位相のデータのみを抽出した時系列データ（即ち、時間間隔ｃでサンプリングした時系列データ）よりＡＲモデルを作成し、その差分をとる。即ち、このときの残差ｖ（ｔ）は次式（１９），（２０）で表される。
ｖ（ｔ）＝ΦＰｗ（ｔ）  ・・・（１９）
ΦＰ≡（１−Φ１Ｂｃ−Φ２Ｂ２ｃ−・・・−ΦＰＢＰｃ）  ・・・（２０）
さらに、方法［３］として、例えば、上記の方法［１］及び［２］を併用する方法がある。即ち、このときの残差ｖ（ｔ）は次式（２１）で表される。
ｖ（ｔ）＝ΦＰ（１−Ｂｃ）Ｄｗ（ｔ）  ・・・（２１）
ここで、季節成分が除去された残差時系列ｖ（ｔ）は定常時系列となるため、ＡＲＭＡモデルでモデル化することができる。即ち、ｖ（ｔ）を用いて次式（２２）が成立する。
φｐｖ（ｔ）＝θｑｅ（ｔ）  ・・・（２２）
したがって、ＳＡＲＩＭＡモデルは、上記の式（２２）に、式（１７）及び式（２１）を代入することにより、次式（２３）のように表される。
φｐΦＰ（１−Ｂｃ）Ｄ（１−Ｂ）ｄｘ（ｔ）＝θｑｅ（ｔ）  ・・・（２３）
次に、以上を踏まえた上で、上記の式（２）について説明する。
先述のように、通常のＳＡＲＩＭＡモデルは、季節成分の周期ｃが常に一定であると仮定している。
これはつまり、図２９に例示するように、季節成分の時系列において、時間間隔ｃ毎に同位相の箇所が出現することを意味する。上述したＳＡＲＩＭＡモデルでは、式（１８）〜式（２１）により、前記季節成分を除去している。
一方、本例の予測処理が予測対象とする、腫瘍位置などの入力時系列信号の季節成分の周期は、図３０に例示するように、変動しているので、図２９に示した例とは異なり、時間間隔ｃ毎に同位相の箇所が出現するとは限らない。
そのため、腫瘍位置などの入力時系列信号に対しては、上記の式（１８）〜式（２１）を用いて、当該入力時系列信号に含まれる季節成分を除去することができない。
そこで、本例では、上記の式（１８）〜式（２１）を改良し、季節成分の周期が変動する非定常時系列に対しても適切に機能するようにした。これが、上述した式（２）である。
上記の改良点としては、例えば、式（１８）のＤの値を０とすることにより、式（１８）の処理を行なわないようにした。
また、式（１９）及び式（２０）の処理に用いる、同位相のデータのみを抽出した時系列データを、通常のＳＡＲＩＭＡモデルのように時間間隔ｃでサンプリングすることにより作成するのではなく、周期変動分析部１２により算出される同位相時刻情報を用いて作成するようにした。
なお、上記の式（２）を用いて有効な時系列予測を行なうためには、モデル係数φ，θ，Φと、モデル次数ｐ，ｑ，ｄの値とを、適切に設計することが求められる。
本例では、既述のように、モデル係数φ，θ，Φの設計方法として、例えば、ユール・ウォーカー法，バーグ法，最小二乗法及び勾配法などの方法を用いることができ、また、モデル次数ｐ，ｑ，ｄの設計方法として、例えば、ＡＩＣや最尤法などの方法を用いることができる。
【００７１】
なお、上述のように、Ｂは時間パラメータを遅延させる演算子である。例えば、時間パラメータｔを持つ変数ｘに前から作用させることで、
【００７２】
【数１７】
0005604306.tif 000033
の関係が成り立つ。また、
0005604306.tif 000034
は、制御部１への入力時系列信号
0005604306.tif 000035
と、モデル設計情報
0005604306.tif 000036
とから、サブモデルｓで対象とするフィルタ処理後の時系列信号
0005604306.tif 000037
を算出する関数である。ここで、
0005604306.tif 000038
は任意定数である。また、モデル次数ｐ、ｑ、ｄと、モデルパラメータ
【００７３】
【数１８】
0005604306.tif 000039
は、
0005604306.tif 000040
の信頼度情報
0005604306.tif 000041
と、モデル設計情報
0005604306.tif 000042
と、ＡＩＣや最尤法の評価アルゴリズムとを用いて決定する。
例えば、始めに、ある任意定数
0005604306.tif 000043
に対して、
【００７４】
【数１９】
0005604306.tif 000044
の条件を満たす場合は、入力時系列信号の季節ダイナミクスが予測に適していないと判断して、
【００７５】
【数２０】
0005604306.tif 000045
を出力してモデル作成処理を終了する。ここで、
0005604306.tif 000046
は、出力算出過程の情報
0005604306.tif 000047
の要素情報であり、値を０（出力値の更新の無しに対応）とする。
一方、もし上記の条件が満たされない場合は、モデルパラメータの設計アルゴリズムを
【００７６】
【数２１】
0005604306.tif 000048
とする。また、モデル次数ｐ、ｑ、ｄをＡＩＣや最尤法で決定した出力
0005604306.tif 000049
を算出し、出力算出過程の要素情報を、出力値の更新の有りに対応する
【００７７】
【数２２】
0005604306.tif 000050
とする。本例において、ｐ＝０、ｄ＝０の場合、サブモデルｓの出力は
【００７８】
【数２３】
0005604306.tif 000051
となる。なお、ここで、
0005604306.tif 000052
が実数の場合は、入力を四捨五入して出力する関数
0005604306.tif 000053
を用いて、
0005604306.tif 000054
と
0005604306.tif 000055
との線形近似により、
0005604306.tif 000056
を求める。
即ち、サブモデル１１１−１〜１１１−Ｓは、モデル設計情報（フィルタやモデルの種類、また、そのパラメータの情報）を用いることで、入力時系列信号（またはその特定成分）の所定時刻経過後の値を予測する予測モデルの一例として機能する。
また、サブモデル１１１−１〜１１１−Ｓは、異なる構造や特性を有していてもよい。
【００７９】
（１．３．２）合成処理部１１２
合成処理部１１２は、サブモデル１１１−１〜１１１−Ｓの出力を加算・乗算処理することにより、制御部１への入力時系列信号の予測対象時刻（ｔｃ＋Ｌ）の予測値を算出する。図５および図１１（Ａ）に例示する合成処理部１１２は、全てのサブモデル１１１−１〜１１１−Ｓの出力
【００８０】
【数２４】
0005604306.tif 000057
を用いて、制御部１への入力時系列信号の予測対象時刻（ｔｃ＋Ｌ）の予測値を次式のように算出する。
【００８１】
【数２５】
0005604306.tif 000058
ここで、関数
0005604306.tif 000059
の具体的な仕様は、必要に応じて、モデル設計情報
0005604306.tif 000060
を考慮しながら、任意に決定することができる。
（１．４）周期変動分析部１２
周期変動分析部１２は、入力時系列信号の周期変動が原因で、制御部１の予測精度が低下してしまう問題を回避するために有効な情報を分析・算出する。
図５及び図６（Ａ）に例示する周期変動分析部１２は、相関解析部１２１と、キャッシュメモリ１２２と、周期推定部１２３と、周期予測部１２４と、同位相時刻推定部１２５と、時間遅れ作用部１２６と、キャッシュメモリ１２７と、差分処理部１２８とを有する。
【００８２】
本例の周期変動分析部１２は、時系列信号
0005604306.tif 000061
を入力として、予測対象時刻（ｔｃ＋Ｌ）の信号
0005604306.tif 000062
の位相と同位相となる入力時系列信号の位置の時刻
【００８３】
【数２６】
0005604306.tif 000063
【００８４】
【数２７】
0005604306.tif 000064
と、その信頼度
【００８５】
【数２８】
0005604306.tif 000065
【００８６】
【数２９】
0005604306.tif 000066
と、周期の予測精度情報
【数３０】
0005604306.tif 000067
【数３１】
0005604306.tif 000068
とを算出する。ここで、
0005604306.tif 000069
は、予測対象時刻（ｔｃ＋Ｌ）の季節ダイナミクスの位相と同位相となる入力時系列の位置の時刻とのうち、ｍ番目に大きい値であり、図６（Ｂ）に例示するような関係がある。また、
0005604306.tif 000070
は、
0005604306.tif 000071
の信頼度であり、
0005604306.tif 000072
は任意の正定数である。また、
0005604306.tif 000073
は、時刻τの周期（推定値）
0005604306.tif 000074
を時刻（τ−Ｌ）までの時系列信号
0005604306.tif 000075
の情報を用いて予測した場合の予測値
0005604306.tif 000076
の精度情報であり、図６（Ｃ）に例示するような関係がある。
即ち、本例の周期変動分析部１２は、入力時系列信号の或る成分（例えば、主要成分）を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析部の一例として機能する。
【００８７】
つまり、本例の周期変動分析部１２は、主要成分が周期時変の季節ダイナミクスである制御部１への入力時系列信号に対して、その入力時系列信号の周期変動が原因で、その所定時間経過後の予測精度が低下してしまう問題を回避するために有効な情報を算出することができる。
【００８８】
（１．４．１）相関解析部１２１
相関解析部１２１は、入力時系列信号の相関解析により、入力時系列信号の周期変動を推定または予測するために必要な相関解析結果を算出する。
【００８９】
図５および図７（Ａ）に例示する相関解析部１２１は、時系列信号
0005604306.tif 000077
を入力として、その相関解析の結果
【数３２】
0005604306.tif 000078
【数３３】
0005604306.tif 000079
を出力する。ここで、
0005604306.tif 000080
は、入力時系列信号の通常の季節ダイナミクスで起こり得ると想定する周期変動の最大値の推定値である。また、
0005604306.tif 000081
の値は、
【００９０】
【数３４】
0005604306.tif 000082
により算出され、この値は図７（Ｂ）、（Ｃ）の太波線部分の時系列信号
0005604306.tif 000083
と、図７（Ｃ）の実線部分の時系列信号
0005604306.tif 000084
との相関係数であり、入力時系列信号に季節ダイナミクスが含まれる場合は、ｋの値を変化させることで図７（Ｄ）に例示するように変化する。
即ち、本例の相関解析部１２１は、入力時系列信号の一部と前記入力時系列信号とのする相関解析を行なう相関解析部の一例として機能する。例えば、本例の相関解析部１２１は、入力時系列信号に対する相関解析により、入力時系列信号の周期変動を推定または予測するために用いる相関解析結果を算出して、キャッシュメモリ１２２へ当該相関解析結果を出力することができる。
【００９１】
（１．４．２）キャッシュメモリ１２２
キャッシュメモリ１２２は、相関解析部１２１の出力である相関解析結果
0005604306.tif 000085
を格納するメモリであり、過去一定期間に格納された情報のうち任意のものを出力する。
図５に例示するキャッシュメモリ１２２は、時刻τ＝１から現在時刻ｔｃまでに格納された全情報の時系列
0005604306.tif 000086
を出力する。
【００９２】
（１．４．３）周期推定部１２３
周期推定部１２３は、入力時系列信号の相関解析結果の時系列情報を用いて、入力時系列信号の周期の変動の推定値の時系列情報と、また、その推定値の信頼度の時系列情報とを算出する。
【００９３】
図５及び図８（Ａ）に例示する周期推定部１２３は、キャッシュメモリ１２２の出力である相関解析結果の時系列情報
0005604306.tif 000087
を入力として、制御部１への入力時系列信号の周期の推定値の時系列情報
0005604306.tif 000088
【数３５】
0005604306.tif 000089
と、その信頼度の時系列情報
0005604306.tif 000090
【数３６】
0005604306.tif 000091
とを出力する。ここで、入力時系列信号のｊ次元方向の時刻τの位相（図８（Ｂ）の星印に対応）の周期
0005604306.tif 000092
とは、時系列信号
0005604306.tif 000093
の季節ダイナミクス成分の時刻τから直近１周期分の長さに対応し、また、
0005604306.tif 000094
はその値の信頼度であり、図８（Ｂ）に例示するように、ｔｓ＝τからｔｓ＝τ＋ｍａｘ（Ｐｊ）までの範囲で求まる値の重み平均として、それぞれｊ次元方向の時系列信号に対する相関解析結果
0005604306.tif 000095
より、次式で算出する。
【００９４】
【数３７】
0005604306.tif 000096
【数３８】
0005604306.tif 000097
【数３９】
0005604306.tif 000098
【数４０】
0005604306.tif 000099
【数４１】
0005604306.tif 000100
ここで、
0005604306.tif 000101
はｋを０から
0005604306.tif 000102
まで徐々に変化させた場合に
0005604306.tif 000103
の値が最初に凸ピークとなるｋの値を出力する関数である（図８（Ｃ）を参照）。
即ち、本例の周期推定部１２３は、前記相関解析部１２１での解析結果を用いて、入力時系列信号の周期変動と、前記周期変動の信頼度とを推定する周期推定部の一例として機能する。例えば、本例の周期推定部１２３は、キャッシュメモリ１２２の出力である入力時系列信号の相関解析結果の時系列情報により、入力時系列信号の周期の変動の推定値の時系列情報と、またその推定値の信頼度の時系列情報とを算出することができる。
【００９５】
（１．４．４）周期予測部１２４
周期予測部１２４は、入力時系列信号の現時刻ｔｃまでの周期の変動の推定値の時系列情報と、その推定値の信頼度の時系列情報を用いて、時刻（ｔｃ＋１）から予測対象時刻（ｔｃ＋Ｌ）までの周期を予測する。
図５及び図９に例示する周期予測部１２４は、周期推定部１２３の出力である入力時系列信号の現時刻ｔｃまでの周期の変動の推定値の時系列情報
0005604306.tif 000104
とその信頼度の時系列情報
0005604306.tif 000105
とを入力として、現時刻ｔｃにおける
0005604306.tif 000106
の予測時系列情報
0005604306.tif 000107
【００９６】
【数４２】
0005604306.tif 000108
を出力する。ここで、
0005604306.tif 000109
は時刻ｔｃにおいて予想した、時刻（ｔｃ＋１）のｊ次元方向の信号の周期
0005604306.tif 000110
の予測値である。この値は、例えば、周期予測部１２４の入力である周期の推定値
0005604306.tif 000111
のデータうち、信頼度が任意定数
0005604306.tif 000112
以上である、すなわち適正なデータ
【００９７】
【数４３】
0005604306.tif 000113
を用いて、次式のような自己回帰モデルにより求める。
【００９８】
【数４４】
0005604306.tif 000114
ここで、
0005604306.tif 000115
は、時系列信号
0005604306.tif 000116
の通常の季節ダイナミクスで起こり得る周期変動の最小値の推定値である。また、
0005604306.tif 000117
と
【００９９】
0005604306.tif 000118
とはモデルの次数と自己回帰係数であり、例えば、ユール・ウォーカー法，バーグ法，最小二乗法及び勾配法などのモデル設計アルゴリズムと、ＡＩＣや最尤法などのモデル評価アルゴリズムとにより決定することができる。
即ち、本例の周期予測部１２４は、周期推定部１２３により推定した前記周期変動及び前記周期変動の信頼度を用いて、前記入力時系列信号の前記所定時間後までの周期を予測する周期予測部の一例として機能する。例えば、本例の周期予測部１２４は、周期推定部１２３の出力である、入力時系列信号の現時刻ｔｃまでの周期の変動の推定値の時系列情報と、その推定値の信頼度の時系列情報を用いた、ユール・ウォーカー法やバーグ法や最小二乗法や勾配法などのモデル作成アルゴリズムと、ＡＩＣや最尤法などのモデル評価アルゴリズムとにより予測モデルを作成し、その予測モデルを用いて予測を行なうことで、時刻（ｔｃ＋１）から予測対象時刻（ｔｃ＋Ｌ）までの周期を予測することができる。
【０１００】
（１．４．５）同位相時刻推定部１２５
同位相時刻推定部１２５は、制御部１への入力時系列信号と、入力時系列信号の相関解析結果の時系列情報と、入力時系列信号の周期の予測対象時刻までの予測結果の時系列情報を用いて、入力時系列信号の周期変動が原因でモデル作成部１１の予測精度が低下してしまう問題を回避するために有効な情報である、予測対象時刻の位相と同位相の入力時系列信号の位相の同位相時刻情報とその信頼度情報を出力する。
【０１０１】
図５及び図１０（Ａ）に例示する同位相時刻推定部１２５は、入力時系列信号
0005604306.tif 000119
と、キャッシュメモリ１２２の出力である相関解析結果の時系列情報
0005604306.tif 000120
と、入力時系列信号の周期の予測対象時刻までの予測結果の時系列情報
0005604306.tif 000121
とを入力として、予測対象時刻の位相と同位相の入力時系列信号の位相の時刻
【０１０２】
【数４５】
0005604306.tif 000122
【０１０３】
【数４６】
0005604306.tif 000123
と、その信頼度
【０１０４】
【数４７】
0005604306.tif 000124
【０１０５】
【数４８】
0005604306.tif 000125
とを出力する。ここで、
0005604306.tif 000126
は、予測対象時刻（ｔｃ＋Ｌ）の季節ダイナミクスの位相の予測値と同位相の時系列信号
0005604306.tif 000127
の時刻のうちｍ番目に大きい値であり、図６（Ｂ）に例示するような関係がある。また、
0005604306.tif 000128
は、
0005604306.tif 000129
の信頼度であり、
0005604306.tif 000130
は任意の正定数である。
0005604306.tif 000131
及び
0005604306.tif 000132
の算出方法を以下に説明する。
【０１０６】
まず、任意の定数
0005604306.tif 000133
に対して、
【０１０７】
【数４９】
0005604306.tif 000134
の条件を満たすかどうかを検証し、条件を満たす場合は、直近の季節ダイナミクスが予測に適さないと判断し、
【０１０８】
【数５０】
0005604306.tif 000135
【０１０９】
として、同位相時刻推定部１２５の処理を終了する。
一方、上記の条件を満たさない場合は、予測対象時刻（ｔｃ＋Ｌ）の位相（図１０（Ｂ）の一番右の星印）と同位相の入力時系列信号（図１０（Ｂ）の星印）の位相の時刻のうち、現時刻ｔｃ以下で、かつ、ｔｃに最近傍（図１０（Ｂ）の四角で囲まれた一番左の星印）の時刻が、図１０（Ｂ）に例示するような関係があることから、次式のように算出する。
【０１１０】
【数５１】
0005604306.tif 000136
【数５２】
0005604306.tif 000137
【数５３】
0005604306.tif 000138
【０１１１】
ここで、
0005604306.tif 000139
は、入力を四捨五入して出力する関数であり、
0005604306.tif 000140
は、
【０１１２】
【数５４】
0005604306.tif 000141
を満たすｉの最小値である。また、
0005604306.tif 000142
の信頼度を
【数５５】
0005604306.tif 000143
とする。
さらに、ｍ＝２，３，…，ＭＴにおける
0005604306.tif 000144
（図１０（Ｄ）の左から１、２番目の星印の時刻）と、その信頼度
0005604306.tif 000145
とは、適切なｔｓの範囲で求まる値の重み平均として、次式で算出する。
【０１１３】
【数５６】
0005604306.tif 000146
0005604306.tif 000147
ここで、
0005604306.tif 000148
は、図１０（Ｃ）に例示するように、ｋを０から
0005604306.tif 000149
まで徐々に変化させた場合に、
0005604306.tif 000150
の値がｎ回目に形成する凸ピークのｋの値を出力する関数である。
即ち、本例の同位相時刻推定部１２５は、相関解析部１２１での相関解析の結果と、周期予測部１２４により予測した前記周期とを用いて、前記同位相時刻情報と前記同位相時刻情報の信頼度とを推定する同位相時刻推定部の一例として機能する。例えば、本例の同位相時刻推定部１２５は、制御部１への入力時系列信号と、キャッシュメモリ１２２の出力である入力時系列信号の相関解析結果の時系列情報と、周期予測部１２４の出力である入力時系列信号の周期の予測対象時刻までの予測結果の時系列情報とを用いて、制御部１の予測精度が低下してしまう問題を回避するために有効な情報である、予測対象時刻の位相と同位相の入力時系列信号における同位相時刻情報と、その信頼度情報とを推定することができる。
【０１１４】
（１．４．６）時間遅れ作用部１２６
時間遅れ作用部１２６は、時間パラメータを持つ任意の入力情報に対して、その時間を任意の期間だけ遅らせた情報を出力する。
図５に例示する時間遅れ作用部１２６は、周期予測部１２４の出力である時刻（ｔｃ＋１）から予測対象時刻（ｔｃ＋Ｌ）までの周期の予測情報
0005604306.tif 000151
を入力として、その時間パラメータをＬだけ遅らせた
0005604306.tif 000152
を出力する。
【０１１５】
（１．４．７）キャッシュメモリ１２７
キャッシュメモリ１２７は、時間遅れ作用部１２６の出力である周期の予測値の時系列情報
0005604306.tif 000153
を格納するメモリであり、過去一定期間に格納された情報のうち任意のものを出力する。
図５に例示するキャッシュメモリ１２７は、過去にキャッシュされた情報のうち、周期の予測期間がＬである予測結果の時系列情報
0005604306.tif 000154
を出力する。
【０１１６】
（１．４．８）差分処理部１２８
図５に例示する差分処理部１２８は、周期推定部１２３の出力の一部である
0005604306.tif 000155
と、キャッシュメモリ１２７の出力である
0005604306.tif 000156
とについて差分処理を施し、その予測誤差（精度）の時系列情報
【０１１７】
【数５７】
0005604306.tif 000157
を出力する。
即ち、本例の差分処理部１２８は、前記周期推定部１２３により推定した前記周期変動と、前記周期予測部１２４により予測した前記周期とを用いて、前記予測した周期の精度情報を推定することができる。
【０１１８】
（１．５）変量推定部１３
変量推定部１３は、時間パラメータを持つ任意の入力情報に対して、単位時間当たり変化量を推定し出力する。
【０１１９】
図５に例示する変量推定部１３は、制御部１への入力時系列信号
0005604306.tif 000158
を入力として、次式で推定されるその単位時間当たりの変化量の時系列情報を出力する。
0005604306.tif 000159
【０１２０】
【数５８】
0005604306.tif 000160
【０１２１】
【数５９】
0005604306.tif 000161
【０１２２】
（１．６）キャッシュメモリ１４
キャッシュメモリ１４は、モデル作成部１１の出力である出力算出過程の情報
0005604306.tif 000162
を格納するメモリであり、過去一定期間に格納された情報のうち任意のものを出力する。
図５に例示するキャッシュメモリ１４は、時刻τ＝１から現在時刻ｔｃまでに格納された全情報の時系列
0005604306.tif 000163
を出力する。
【０１２３】
（１．７）キャッシュメモリ１５
キャッシュメモリ１５は、モデル作成部１１の出力である制御部１への入力時系列信号の予測対象時刻（ｔｃ＋Ｌ）の予測値
0005604306.tif 000164
を格納するメモリであり、過去一定期間に格納された情報のうち任意のものを出力する。
図５に例示するキャッシュメモリ１５は、時刻τ＝１から現在時刻ｔｃまでに格納された全情報の時系列
0005604306.tif 000165
を出力する。
【０１２４】
（１．８）予測精度推定部１６
予測精度推定部１６は、制御部１への入力時系列信号と、変量推定部１３の出力である入力時系列信号の単位時間当たりの変化量の時系列情報と、モデル作成部１１の出力である入力時系列信号の予測対象時刻の予測値とその算出過程の時系列情報と、周期変動分析部１２の出力である入力時系列信号の予測対象時刻の位相と同位相の同位相時刻情報とその信頼度情報と、予測対象時刻までの周期変動の予測精度の時系列情報を用いて、過去の予測における誤差量推移や各サブモデル１１１−１〜１１１−Ｓの有効性（予測における貢献）を評価し、また、入力時系列信号の予測対象時刻の単位時間当たりの変化量を推定することで、入力時系列信号の予測対象時刻の予測値の精度情報と、対象予測時刻の最適な予測のための各サブモデル１１１−１〜１１１−Ｓと合成処理部１１２のモデル設計情報を出力する。
【０１２５】
予測精度推定部１６が、入力時系列信号の予測対象時刻の予測値の精度情報と、対象予測時刻の最適な予測のためのモデル設計情報を推定するアイデアを以下に述べる。
一般的に、時系列信号の単位時間当たりの変化量が大きい場合は、予測精度が低く、遅い場合は予測精度が高くなる。また、その変化量の違いにより、予測に適するモデルの種類やモデル設計用パラメータも異なる。すなわち、予測対象時刻の時系列信号の単位時間当たりの変化量を推定することは、予測対象時刻の予測精度を推定し、また、モデルの種類やモデル設計用パラメータを選択するための重要な情報になりうる。そこで、本例では、変量推定部１３の出力である、入力時系列信号の単位時間当たりの変化量の時系列情報を用いて、予測対象時刻の時系列信号の単位時間当たりの変化量を推定する。
【０１２６】
また、一般的に、入力時系列信号の主要成分が季節ダイナミクスの場合、予測対象時刻の予測精度と、予測対象時刻の位相と同位相の過去の位相の時刻における予測精度に一定の相関がある。そこで、本例の予測精度推定部１６は、例えば、入力時系列信号の予測対象時刻の位相と同位相の時刻に対応する同位相時刻の信頼度を重みとした場合の、入力時系列信号とその予測対象時刻の予測値とその算出過程の時系列情報と先に推定した予測対象時刻の時系列信号の単位時間当たりの変化量の情報により各サブモデル１１１−１〜１１１−Ｓと合成処理部１１２を様々な設計パラメータを用いた場合の重み付き平均誤差量を算出し、その評価と、先に推定した予測対象時刻の時系列信号の単位時間当たりの変化量の情報を用いて、最適なモデル設計情報を推測する。また、その最適なモデル設計情報に対応する重み付き平均誤差量に、必要ならば任意の補正を行なった結果を、予測値の精度情報の推測値とする。
【０１２７】
図５および図１２（Ａ）に例示する予測精度推定部１６は、入力時系列信号
0005604306.tif 000166
と、変量推定部１３の出力である入力時系列信号の単位時間当たりの変化量の時系列情報
0005604306.tif 000167
と、モデル作成部１１の出力である、これまでの入力時系列信号に対する予測値の時系列情報
0005604306.tif 000168
とその算出過程の時系列情報
0005604306.tif 000169
と、周期変動分析部１２の出力である入力時系列信号の予測対象時刻の位相と同位相の同位相時刻情報
0005604306.tif 000170
と、その信頼度情報
0005604306.tif 000171
と、予測対象時刻までの周期変動の予測精度の時系列情報
0005604306.tif 000172
とを入力として、入力時系列信号の予測対象時刻の予測値の精度情報である予測の誤差量の推定値
【０１２８】
【数６０】
0005604306.tif 000173
と時系列信号の単位時間当たりの変化量の推定値
【０１２９】
【数６１】
0005604306.tif 000174
と、予測対象時刻の最適な予測のための各サブモデル１１１−１〜１１１−Ｓと合成処理部１１２のモデル設計情報
0005604306.tif 000175
とを出力する。
例えば、まず、任意に用意したＵ種類のモデル設計情報候補のうち、ｕ番目の候補
【０１３０】
0005604306.tif 000176
に対する、入力時系列信号の予測対象時刻の誤差量と入力時系列信号の単位時間当たりの変化量の推定値を以下のように定義する。
【０１３１】
【数６２】
0005604306.tif 000177
【数６３】
0005604306.tif 000178
ここで、
0005604306.tif 000179
と
0005604306.tif 000180
とはそれぞれ、サブモデル１１１−１〜１１１−Ｓと合成処理部１１２とに依存する関数であり、また、その引数である
【０１３２】
【数６４】
0005604306.tif 000181
は予測精度推定部１６への全入力情報である。
次に、全候補のなかで予測誤差量が最小になる候補の番号
【０１３３】
【数６５】
0005604306.tif 000182
を選出する。ただし、
【０１３４】
【数６６】
0005604306.tif 000183
である。この選出候補の番号より、モデル設計情報が
【０１３５】
【数６７】
0005604306.tif 000184
のように決まり、また、予測対象時刻の誤差量と時系列信号の単位時間当たりの変化量の推定値も次式のように求まる。
【０１３６】
【数６８】
0005604306.tif 000185
【数６９】
0005604306.tif 000186
ここで、
0005604306.tif 000187
と
0005604306.tif 000188
とはそれぞれ任意の定数である。
さらに、より具体的な例を示すために、サブモデル１１１−１〜１１１−Ｓの数Ｓを１と仮定する。この場合、例えば、入力時系列信号の予測対象時刻の位相と同位相の入力時系列の位相の時刻付近における実測値の重み付き平均量を用いて、予測対象時刻の誤差量と時系列信号の単位時間当たりの変化量の推定値を次式のように算出してもよい。
【０１３７】
【数７０】
0005604306.tif 000189
【数７１】
0005604306.tif 000190
【数７２】
0005604306.tif 000191
【数７３】
0005604306.tif 000192
【数７４】
0005604306.tif 000193
ここで、
0005604306.tif 000194
と
0005604306.tif 000195
とは、それぞれ通常の入力時系列信号に対する予測で起こりえる誤差量と時系列信号の単位時間当たりの変化量の最大値との推定値である。また、図１２（Ｂ）に例示する関数
0005604306.tif 000196
は予測対象時刻の周期予測の誤差の影響を考慮した任意の重み関数であり、その引数ｄは次式のように決まる。
【０１３８】
【数７５】
0005604306.tif 000197
ここで、
0005604306.tif 000198
は整数の集合であり、
0005604306.tif 000199
は入力を四捨五入して出力する関数である。
また、別の具体例としては、入力時系列信号の軌跡のうち、予測対象時刻の位置の近傍点である過去の時刻付近における実測値の重み付き平均量を用いて、予測対象時刻の誤差量と時系列信号の単位時間当たりの変化量の推定値とを次式のように算出してもよい。
【０１３９】
【数７６】
0005604306.tif 000200
【数７７】
0005604306.tif 000201
【数７８】
0005604306.tif 000202
【数７９】
0005604306.tif 000203
ここで、
0005604306.tif 000204
は任意の正定数である。
なお、これらの具体例のようにサブモデル１１１−１〜１１１−Ｓの数Ｓが１であり、モデルの詳細が一意に決まる場合は、
【数８０】
0005604306.tif 000205
としてもよい。
即ち、本例の予測精度推定部１６は、入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析部１２により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度を用いて、前記入力時系列信号の所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定部の一例として機能する。
【０１４０】
また、本例の予測精度推定部１６は、制御部１への入力時系列信号と、変量推定部１３の出力である入力時系列信号の単位時間当たりの変化量の時系列情報と、モデル作成部１１の出力である入力時系列信号の予測対象時刻の予測値とその算出過程の時系列情報と、周期変動分析部１２の出力である入力時系列信号の予測対象時刻の位相と同位相の同位相時刻情報とその信頼度情報と予測対象時刻までの周期変動の予測精度の時系列情報とを用いて、過去の予測における誤差量推移や各サブモデル１１１−１〜１１１−Ｓの有効性を評価し、また、入力時系列信号の予測対象時刻の単位時間当たりの変量を推定することで、入力時系列信号の予測対象時刻の予測値の精度情報と、対象予測時刻の最適な予測のための各サブモデル１１１−１〜１１１−Ｓと合成処理部１１２のモデル設計情報とを推定することができる。
【０１４１】
さらに、本例の予測精度推定部１６は、モデル作成部１１により出力される、前記入力時系列信号の前記所定時間後の値の算出過程に関する情報を用いて、前記モデル設計情報を推定することができる。
また、本例の予測精度推定部１６は、上記予測モデルを用いて入力時系列信号の所定時間後の値を予測した場合の予測精度情報を推定することができる。
（１．９）本制御の流れの一例
ここで、上記制御の流れの一例を、図１９〜図２８に例示するフローチャートを用いて説明する。
まず、図１９に例示するように、上記制御が開始されると（ステップＳ１）、制御部１は、入力時系列信号の所定時間後の予測に用いる初期パラメータを決定する（ステップＳ２）。この初期パラメータには、例えば、現時刻としてｔｃ＝１，モデル設計情報として
0005604306.tif 000206
及び各種の任意定数としてＬ，ＴＥ，Ｍｔが含まれる。ここで、ＴＥは、予測終了時刻を表す。
そして、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000207
が制御部１に入力されると（ステップＳ３）、制御部１は、周期変動分析部１２へ処理を移行する（ステップＳ４）。
ステップＳ４において、処理が周期変動分析部１２へ移行されると、図２０に例示するように、周期変動分析部１２が処理を開始する（ステップＳ１２）。
そして、周期変動分析部１２は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000208
を入力されると（ステップＳ１３）、相関解析部１２１，周期推定部１２３，周囲予測部１２４，同位相時刻推定部１２５及び時間遅れ作用部１２６へ処理を順次に移行させ（ステップＳ１４〜Ｓ１８）、周期変動分析部１２での処理を終了する（ステップＳ１９）。
ステップＳ１４において、処理が相関解析部１２１へ移行されると、図２１に例示するように、相関解析部１２１が処理を開始する（ステップＳ２０）。
そして、相関解析部１２１は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000209
を入力されると（ステップＳ２１）、前記入力時系列信号の一部と前記入力時系列信号との相関解析の結果
0005604306.tif 000210
を算出して出力する（ステップＳ２２）。
さらに、相関解析部１２１は、キャッシュメモリ１２２に保存された相関解析の結果の時系列データ
0005604306.tif 000211
を出力して（ステップＳ２３）、相関解析部１２１での処理を終了する（ステップＳ２４）。
次に、ステップＳ１５において、処理が周期推定部１２３へ移行されると、図２２に例示するように、周期推定部１２３が処理を開始する（ステップＳ２５）。
そして、周期推定部１２３は、相関解析部１２１から、例えば、上記相関解析の結果の時系列データ
0005604306.tif 000212
を入力されると（ステップＳ２６）、前記入力時系列信号に含まれる季節成分の周期の推定値の時系列データ
0005604306.tif 000213
と、その信頼値の時系列データ
0005604306.tif 000214
とを算出して出力し（ステップＳ２７）、周期推定部１２３での処理を終了する（ステップＳ２８）。
次に、ステップＳ１６において、処理が周期予測部１２４へ移行されると、図２３に例示するように、周期予測部１２４が処理を開始する（ステップＳ２９）。
そして、周期予測部１２４は、周期推定部１２３から、前記入力時系列信号に含まれる季節成分の周期の推定値の時系列データ
0005604306.tif 000215
と、その信頼値の時系列データ
0005604306.tif 000216
とを入力されると（ステップＳ３０）、前記入力時系列信号に含まれる季節成分の周期を、予測対象時刻（ｔｃ＋Ｌ）まで予測した時系列データ
0005604306.tif 000217
を算出して出力し（ステップＳ３１）、周期予測部１２４での処理を終了する（ステップＳ３２）。
次に、ステップＳ１７において、処理が同位相時刻推定部１２５へ移行されると、図２４に例示するように、同位相時刻推定部１２５が処理を開始する（ステップＳ３３）。
そして、同位相時刻推定部１２５は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000218
を入力され（ステップＳ３４）、また、相関解析部１２１から上記相関解析の結果の時系列データ
0005604306.tif 000219
を入力され（ステップＳ３５）、さらに、周期予測部１２４から上記周期予測結果の時系列データ
0005604306.tif 000220
を入力されると（ステップＳ３６）、予測対象時刻の前記入力時系列信号（例えば、予測対象時刻の腫瘍位置など）
0005604306.tif 000221
に含まれる季節成分の位相と同位相である時刻
0005604306.tif 000222
と、その信頼度
0005604306.tif 000223
とを算出して出力し（ステップＳ３７）、同位相時刻推定部１２５での処理を終了する（ステップＳ３８）。なお、予測対象時刻の前記入力時系列信号
0005604306.tif 000224
に含まれる季節成分の位相と同位相である時刻は、例えば、新しいものから順にＭｔ個算出されるようにしてもよい。
次に、ステップＳ１８において、処理が時間遅れ作用部１２６へ移行されると、図２５に例示するように、時間遅れ作用部１２６が処理を開始する（ステップＳ３９）。
そして、時間遅れ作用部１２６は、周期予測部１２４から上記周期予測結果の時系列データ
0005604306.tif 000225
を入力されると（ステップＳ４０）、当該入力時系列
0005604306.tif 000226
の時間パラメータをＬだけ遅らせた
0005604306.tif 000227
を算出して出力し（ステップＳ４１）、時間遅れ作用部１２６での処理を終了する（ステップＳ４２）。
周期変動分析部１２での処理が終了すると、図１９に戻り、制御部１は、変量推定部１３へ処理を移行させる（ステップＳ５）。
ステップＳ５において、処理が変量推定部１３へ移行されると、図２６に例示するように、変量推定部１３が処理を開始する（ステップＳ４３）。
そして、変量推定部１３は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000228
を入力されると（ステップＳ４４）、当該入力時系列信号の単位時間当たりの変化量の時系列データ
0005604306.tif 000229
を算出して出力し（ステップＳ４５）、変量推定部１３での処理を終了する（ステップＳ４６）。
変量推定部１３での処理が終了すると、図１９に戻り、制御部１は、モデル作成部１１へ処理を移行させる（ステップＳ６）。
ステップＳ６において、処理がモデル作成部１１へ移行されると、図２７に例示するように、モデル作成部１１が処理を開始する（ステップＳ４７）。
そして、モデル作成部１１は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000230
を入力され（ステップＳ４８）、また、同位相時刻推定部１２５から、予測対象時刻の前記入力時系列信号（例えば、予測対象時刻の腫瘍位置など）
0005604306.tif 000231
に含まれる季節成分の位相と同位相である時刻
0005604306.tif 000232
と、その信頼度
0005604306.tif 000233
とを入力され（ステップＳ４９）、さらに、制御部１からモデル設計情報
0005604306.tif 000234
とを入力されると（ステップＳ５０）、当該モデル設計情報
0005604306.tif 000235
をもとに、サブモデル１１１−１〜１１１−Ｓを設計し、他の入力情報
0005604306.tif 000236
及び
0005604306.tif 000237
を用いて、前記入力時系列信号の所定時間後の予測時系列信号（例えば、予測腫瘍位置など）
0005604306.tif 000238
を算出し出力するとともに、その算出過程に関する詳細情報として
0005604306.tif 000239
を出力し（ステップＳ５１）、モデル作成部１１での処理を終了する（ステップＳ５２）。
モデル作成部１１での処理が終了すると、図１９に戻り、制御部１は、予測精度推定部１６へ処理を移行させる（ステップＳ７）。
ステップＳ７において、処理が予測精度推定部１６へ移行されると、図２８に例示するように、予測精度推定部１６が処理を開始する（ステップＳ５３）。
そして、予測精度推定部１６は、制御部１から、例えば、腫瘍位置の時系列データとして、入力時系列信号
0005604306.tif 000240
を入力され（ステップＳ５４）、また、変量推定部１３から当該入力時系列信号の単位時間当たりの変化量の時系列データ
0005604306.tif 000241
を入力され（ステップＳ５５）、さらに、同位相時刻推定部１２５から、予測対象時刻の前記入力時系列信号（例えば、予測対象時刻の腫瘍位置など）
0005604306.tif 000242
に含まれる季節成分の位相と同位相である時刻
0005604306.tif 000243
と、その信頼度
0005604306.tif 000244
とを入力され（ステップＳ５６）、また、モデル作成部１１から、キャッシュメモリ１４および１５を経由して、過去に行なわれた予測結果の時系列データ
0005604306.tif 000245
と、当該予測の算出過程に関する詳細情報
0005604306.tif 000246
とを入力されると（ステップＳ５７，Ｓ５８）、上記の各入力情報を用いて、予測に最適なモデル設計情報
0005604306.tif 000247
を決定し出力するとともに、予測腫瘍位置
0005604306.tif 000248
の誤差量の予測値
0005604306.tif 000249
と、予測腫瘍位置
0005604306.tif 000250
における単位時間当たりの変化量の予測値
0005604306.tif 000251
とを算出して出力し（ステップＳ５９）、予測精度推定部１６での処理を終了する（ステップＳ６０）。
予測精度推定部１６での処理が終了すると、図１９に戻り、制御部１は、腫瘍位置の予測値
0005604306.tif 000252
と、予測腫瘍位置
0005604306.tif 000253
における単位時間当たりの変化量の予測値（腫瘍速度の予測値）
0005604306.tif 000254
と、予測腫瘍位置
0005604306.tif 000255
の誤差量の予測値
0005604306.tif 000256
とを出力する（ステップＳ８）。
そして、制御部１は、現時刻ｔｃをインクリメントすることにより、時刻を更新し（ステップＳ９）、前記更新した時刻が予測終了時刻ＴＥよりも小さいかどうかを判断する（ステップＳ１０）。
ここで、制御部１が、前記更新した時刻が予測終了時刻ＴＥよりも小さいと判断した場合（ステップＳ１０のＹｅｓルート）、図１９中のステップＳ３〜ステップＳ９までの処理を繰り返す。
一方、制御部１が、前記更新した時刻が予測終了時刻ＴＥ以上であると判断した場合（ステップＳ１０のＮｏルート）、上記予測処理を終了する（ステップＳ１１）。
【０１４２】
（１．１０）実験例
上記の制御部１による予測機能の有効性を、実際の患部２６０の観測位置の時系列信号を用いた実験例により説明する。
【０１４３】
図１３は、ある患部２６０の観測位置の時系列信号を３次元空間上にプロットした図である。また、この時系列信号の次元方向ごとの各成分をグラフ化したものが図１４に例示する時系列信号となる。ここで、縦軸のスケールはミリメートル（ｍｍ）であり、横軸のスケールは０．０３３秒／ｓｔｅｐである。これより、入力時系列信号中の主要成分が、振る舞いが単調かつ安定的な季節ダイナミクスであることがわかる。なお、相関解析によりこの入力時系列信号の周期が時変であることがわかっている。
【０１４４】
本例では、入力時系列信号に対して周期の長さで移動平均を取ることでトレンド成分を分離し、また、残差時系列信号に対して、周期の長さの時間遅れの値で差分をとることで、季節成分とその他の定常成分に分離した。各成分の時系列信号を、サブモデル１１１−１〜１１１−Ｓ（Ｓ＝３）で対象とするフィルタ処理後の時系列信号として、その３０ステップ〔Ｌ＝３０（約１秒）〕先の値を予測した結果を図１５に、また、それら出力の合成として入力時系列信号を予測した結果を図１６（Ａ）に示す。ここで、青線が実測値、赤線がその予測結果に対応し、図１５中のｘｊＴ，ｘｊＰ，ｘｊＡ，ｊ＝１，２，３は、ｊ次元方向データのトレンド成分、季節成分、定常成分に対応する各サブ時系列である。比較のため、従来の一般的なＳＡＲＩＭＡモデルを用いて、入力時系列信号を予測した結果も図１７に示す。ここで、緑線がＡＲＭＡモデルによる予測結果に対応する。なお、各図では、見やすいように横軸を２０００ステップ〜４０００ステップの区間に限定している。
【０１４５】
図１５の結果より、季節成分に対する予測の精度が、トレンド成分と定常成分の予測精度に比べて良好であることがわかる。本例では、予測精度推定部１６において、過去の予測精度実績に対する評価によりこの状況を検知し、トレンド成分の予測モデルをゼロ次ホールドモデルに、また、全体としての予測結果に季節成分の予測結果が占める割合を大きくし、定常成分の予測結果が占める割合を小さくするモデル設計情報が作成される。図１６（Ａ）に示される結果より、全ての成分を合成することで得られる入力時系列信号の予測結果の精度が、図１７に示す従来の一般的なＳＡＲＩＭＡモデルによるものと比べて、良好であることがわかる。
【０１４６】
また、本例の予測精度推定部１６において推定した予測対象時刻の誤差量と時系列信号の単位時間当たりの変化量の時系列情報を、それぞれ図１６（Ｂ）と（Ｃ）に示す。ここで、青線が実測値、赤線がその推定値に対応する。これらの結果より、本例の予測精度推定部による推定結果が、実測値と高い相関があるのがわかる。
また、予測期間の長さＬを変化させた場合の平均予測誤差量の変化を図１８に示す。図１８の横軸は予測間隔（Ｌ）を示し、縦軸は予測誤差の平均値を示す。なお、図中「Ｚｅｒｏ−ｏｒｄｅｒ  ｈｏｌｄ」と示した結果は、入力時系列信号の最新のデータをそのまま予測値として用いた場合の誤差を表している。この図１８に示すように、５ステップ（０．１６７秒）以上の中長期的な予測では、本例の予測方法が効果的であることが分かる。
【０１４７】
このように、本例の制御部１は、患部２６０の動きの長期（１秒程度）予測とその予測精度情報の推定が可能であるので、例えば、ＭＬＣ２１０の照射野形成制御に要する時間に先だって、ＭＬＣ２１０を駆動制御することが可能となる。結果、患部２６０の動きに追従した連続的な放射線照射を実現することができる。
ところで、入力時系列信号の３７００ステップ付近に、特異的なダイナミクスの変化が生じていることが分かる。本例では、予測精度推定部１６において、最新の予測と対応する過去の予測の結果情報（誤差量と時系列信号の単位時間当たりの変化量）データの照合・解析を行なうため、このような特的なダイナミクスの変化を、例えば、誤差量と時系列信号の単位時間当たりの変化量の推定値が想定する最大値となる（図１６（Ｂ）、（Ｃ）の該当部分の出力が飽和する）ことで検知可能である。すなわち、特異的なダイナミクスの変化が長期に継続して患部２６０の動きの正確な予測と追随による放射線照射が困難な場合には、その検知情報により放射線の照射を停止し、患部２６０以外への正常な細胞への照射を回避することができる。
【０１４８】
これらの結果より、「主要成分が周期時変の季節ダイナミクスである任意の時系列信号を、周期時変の影響を考慮して、さらに、過去の予測結果の情報をもとに、複数のサブモデルを合成としてモデル設計を行ない予測とその精度を推定する」という本例の方法が有効であることがわかる。
また、本例の周期変動分析部１２は、相関解析部１２１の処理
【数８１】
0005604306.tif 000257
における
0005604306.tif 000258
と
0005604306.tif 000259
の値とを任意に複数設定した場合の処理を、同時並行して行なうこともできる。
これは、相関解析部１２１において、相関解析の計算式やデータ範囲を任意に複数設定し、各設定の
0005604306.tif 000260
と、
0005604306.tif 000261
と、
【数５７】
0005604306.tif 000262
とを、同時並行で算出することを意味する。
また、予測誤差推定部は１６、周期変動分析部１２が同時平行で算出した、
0005604306.tif 000263
と、
0005604306.tif 000264
と、
【数５７】
0005604306.tif 000265
の情報とを用いて、より多様なモデル設計情報
0005604306.tif 000266
を推定可能である。
この場合、例えば、モデル作成部１１は、フィルタ及びモデルの種類やパラメータ値だけでなく、
0005604306.tif 000267
と、
0005604306.tif 000268
の値との組み合わせが異なる、多様なサブモデル１１１−１〜１１１−Ｓを作成可能である。
ここで、
0005604306.tif 000269
と、
0005604306.tif 000270
とを複数設定してサブモデルを多様化する上記改良を施す前の本例の予測方法を「本例の予測方法（１）」とし、上記改良を施した後の本例の予測方法を「本例の予測方法（２）」とすると、予測期間長と予測誤差との関係は図３１に示すようになる。この図３１に示すように、従来の予測方法の一例である季節調整指数平滑化法やＳＡＲＩＭＡモデルによる予測方法に比して、本例の予測方法（１）が効果的であり、さらに、本例の予測方法（１）よりも本例の予測方法（２）の方が効果的であることが分かる。
【０１４９】
〔２〕その他
なお、上述した制御部１としての機能は、コンピュータ（ＣＰＵ，情報処理装置，各種端末を含む）が所定のアプリケーションプログラム（信号処理プログラム）を実行することによって実現されてもよい。
【０１５０】
即ち、上記プログラムは、入力時系列信号の所定時間後の値を予測する制御部１において、前記予測機能をコンピュータに実現させるための信号処理プログラムであって、前記信号処理プログラムが、入力時系列信号の或る成分を周期が時間変化する時系列信号と仮定することにより、前記入力時系列信号の一部と前記入力時系列信号との相関解析を行ない、前記相関解析の結果を用いて、前記入力時系列信号の周期変動及び前記周期変動の信頼度を推定し、前記周期変動及び前記周期変動の信頼度を用いて前記入力時系列信号の所定時間後までの周期を予測し、前記相関解析の結果と前記予測した周期とを用いて、前記入力時系列信号の前記所定時間後の位相と同位相となる前記入力時系列信号の時刻である同位相時刻情報及び前記同位相時刻情報の信頼度とを算出する周期変動分析機能と、前記入力時系列信号と、前記入力時系列信号の単位時間当たりの変化量の時系列情報と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度とを用いて、前記入力時系列信号の前記所定時間後の値を予測するためのモデル設計情報と、前記モデル設計情報による前記入力時系列信号の前記所定時間後の値の予測の精度とを推定する予測精度推定機能と、前記入力時系列信号と、前記周期変動分析機能により算出した前記同位相時刻情報及び前記同位相時刻情報の信頼度と、前記予測精度推定機能により推定した前記モデル設計情報とを用いて、前記入力時系列信号の前記所定時間後の値を予測する予測モデルを作成し、前記予測モデルを用いて前記入力時系列信号の前記所定時間後の値を予測して出力するモデル作成機能と、を前記コンピュータに実現させる、信号処理プログラムの一例である。
【０１５１】
また、上記プログラムは、例えば、フレキシブルディスク，ＣＤ（ＣＤ−ＲＯＭ，ＣＤ−Ｒ，ＣＤ−ＲＷなど），ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−ＲＡＭ，ＤＶＤ−Ｒ，ＤＶＤ−ＲＷ，ＤＶＤ＋Ｒ，ＤＶＤ＋ＲＷなど）等のコンピュータ読取可能な記録媒体に記録された形態で提供されうる。この場合、コンピュータはその記録媒体から信号処理プログラムを読み取って内部記憶装置または外部記憶装置に転送し格納して用いることができる。また、そのプログラムを、例えば、磁気ディスク，光ディスク，光磁気ディスク等の記憶装置（記録媒体）に記録しておき、その記憶装置から通信回線を介してコンピュータに提供するようにしてもよい。
【０１５２】
ここで、コンピュータとは、ハードウェアとＯＳ（オペレーティングシステム）とを含む概念であり、ＯＳの制御の下で動作するハードウェアを意味している。また、ＯＳが不要でアプリケーションプログラム単独でハードウェアを動作させるような場合には、そのハードウェア自体がコンピュータに相当する。ハードウェアは、少なくとも、ＣＰＵ等のマイクロプロセッサと、記録媒体に記録されたコンピュータプログラムを読み取るための手段とをそなえている。
【０１５３】
上記信号処理プログラムとしてのアプリケーションプログラムは、上述のようなコンピュータに、制御部１としての機能を実現させるプログラムコードを含んでいる。また、その機能の一部は、アプリケーションプログラムではなくＯＳによって実現されてもよい。
なお、本実施形態としての記録媒体としては、上述したフレキシブルディスク，ＣＤ，ＤＶＤ，磁気ディスク，光ディスク，光磁気ディスクのほか、ＩＣカード，ＲＯＭカートリッジ，磁気テープ，パンチカード，コンピュータの内部記憶装置（ＲＡＭやＲＯＭなどのメモリ），外部記憶装置等や、バーコードなどの符号が印刷された印刷物等の、コンピュータ読み取り可能な種々の媒体を利用することもできる。
【０１５４】
また、上述した放射線治療装置及び制御部１の各構成及び各処理は、必要に応じて取捨選択してもよいし、適宜組み合わせてもよい。
例えば、上述した例では、予測精度推定部１６が、モデル作成部１１により出力される、入力時系列信号の所定時間後の値の算出過程に関する情報を用いて、モデル設計情報を推定する例について説明したが、この推定処理を省略することもできる。
また、上述した例では、入力時系列信号をトレンド成分，季節成分，定常成分に分離し、それぞれについて予測モデルを適用することで、入力時系列信号の予測を行なったが、入力時系列信号から上記のうちいずれかの成分を選択的に分離して予測処理に用いるようにしてもよい。
【０１５５】
例えば、入力時系列信号に含まれる他の信号成分よりも周期性の強い季節成分のみを分離して所定の予測モデルを適用するようにしても、上記と同様の効果を得ることができる。
また、本例の信号処理方法は、例えば、波の高さや、心臓の拍動に起因する位置変動を予測するシステムに適用することもできる。例えば、心臓の拍動に起因する位置変動を上記信号処理方法により予測することにより、通常の状態からどの程度ずれているか（不整脈など）を警告するシステムなどを実現することが可能である。
【０１５６】
さらに、この他、電力需要予測，インターネットの通信量及びサーバ負荷などの予測に適用することもできる。
【産業上の利用可能性】
【０１５７】
以上詳述したように、本発明によれば、従来のように患部位置の観測信号を入力として得られる入力時系列信号に自己回帰モデルを単純適用して予測を行なう場合に比して、患部位置の入力時系列信号を正確且つより中長期に予測することができる。したがって、患部位置の動きに追随した放射線照射を実現することが可能であり、放射線照射技術分野において極めて有用と考えられる。
【符号の説明】
【０１５８】
１  制御部
１１  モデル作成部
１２  周期変動分析部
１３  変量推定部
１４  キャッシュメモリ
１５  キャッシュメモリ
１６  予測精度推定部
１００  クライストロン
１１０  冷却装置
１１１−１  サブモデル１
１１１−２  サブモデル２
１１１−Ｓ  サブモデルＳ
１２０  導波管
１２１  相関解析部
１２２  キャッシュメモリ
１２３  周期推定部
１２４  周期予測部
１２５  同位相時刻推定部
１２６  時間遅れ作用部
１２７  キャッシュメモリ
１２８  差分処理部
１３０  真空ポンプ
１４０  電子銃
１５０  加速管
１６０  偏向電磁石
１７０  ターゲット
１８０  平坦化フィルタ（散乱フォイル）
１９０  モニタ線量計
２００  ＪＡＷ（照射野形成部）
２１０  ＭＬＣ（マルチリーフコリメータ）
２２０  ＯＢＩ
２３０  ベッド
２４０  ＥＰＩＤ
２５０  パルスモジュレータ
３０１  駆動制御部
３０２−１，３０２−２  駆動部
３０３  リーフ
３０４  放射線照射野
３０５  Ｘ線照射部
３０６  センサーパネル
【図２】
0005604306.tif 000271
【図４】
0005604306.tif 000272
【図５】
0005604306.tif 000273
【図６】
0005604306.tif 000274
【図７】
0005604306.tif 000275
【図８】
0005604306.tif 000276
【図９】
0005604306.tif 000277
【図１０】
0005604306.tif 000278
【図１１】
0005604306.tif 000279
【図１２】
0005604306.tif 000280
【図１３】
0005604306.tif 000281
【図１４】
0005604306.tif 000282
【図１５】
0005604306.tif 000283
【図１６】
0005604306.tif 000284
【図１７】
0005604306.tif 000285
【図１８】
0005604306.tif 000286
【図１９】
0005604306.tif 000287
【図２０】
0005604306.tif 000288
【図２１】
0005604306.tif 000289
【図２２】
0005604306.tif 000290
【図２３】
0005604306.tif 000291
【図２４】
0005604306.tif 000292
【図２５】
0005604306.tif 000293
【図２６】
0005604306.tif 000294
【図２７】
0005604306.tif 000295
【図２８】
0005604306.tif 000296
【図２９】
0005604306.tif 000297
【図３０】
0005604306.tif 000298
【図３１】
0005604306.tif 000299
【図１】
0005604306.tif 000300
【図３】
0005604306.tif 000301