"複数の連続した時系列画像を取り込むデータ入力手段と、前記時系列画像を蓄積するデータ蓄積手段と、前記時系列画像の持つ初期エネルギーを計算する初期エネルギー計算手段と、前記時系列画像から、移流拡散方程式とナビエストークス方程式により前記時系列画像よりも未来の予測画像を生成する予測手段と、前記時系列画像の初期エネルギーと前記予測手段が生成した予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように前記予測画像を補正する最適化計算手段とを備えることを特徴とする画像生成装置。請求項1記載の画像生成装置において、前記初期エネルギー計算手段と前記最適化計算手段とは、画像の動きベクトルから波数エネルギーと運動エネルギーと発散渦エネルギーを計算し、画像の輝度分布からエッジエネルギーを計算し、前記波数エネルギーと前記運動エネルギーと前記発散渦エネルギーと前記エッジエネルギーの平均エネルギーを当該画像の持つエネルギーとすることを特徴とする画像生成装置。請求項2記載の画像生成装置において、前記初期エネルギー計算手段と前記最適化計算手段とは、前記波数エネルギーを計算する際に、連続的な波数とエネルギーとの関係をガンマ分布により近似することにより、前記波数エネルギーを計算することを特徴とする画像生成装置。請求項2または3記載の画像生成装置において、前記最適化計算手段は、前記誤差が最小になるまで前記移流拡散方程式の拡散係数と前記ナビエストークス方程式の粘性係数および密度を変化させることにより、前記予測画像を補正することを特徴とする画像生成装置。複数の連続した時系列画像を取り込むデータ入力ステップと、前記時系列画像を蓄積するデータ蓄積ステップと、前記時系列画像の持つ初期エネルギーを計算する初期エネルギー計算ステップと、前記時系列画像から、移流拡散方程式とナビエストークス方程式により前記時系列画像よりも未来の予測画像を生成する予測ステップと、前記時系列画像の初期エネルギーと前記予測ステップで生成した予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように前記予測画像を補正する最適化計算ステップとを備えることを特徴とする画像生成方法。請求項5記載の画像生成方法において、前記初期エネルギー計算ステップと前記最適化計算ステップとは、画像の動きベクトルから波数エネルギーと運動エネルギーと発散渦エネルギーを計算し、画像の輝度分布からエッジエネルギーを計算し、前記波数エネルギーと前記運動エネルギーと前記発散渦エネルギーと前記エッジエネルギーの平均エネルギーを当該画像の持つエネルギーとすることを特徴とする画像生成方法。請求項6記載の画像生成方法において、前記初期エネルギー計算ステップと前記最適化計算ステップとは、前記波数エネルギーを計算する際に、連続的な波数とエネルギーとの関係をガンマ分布により近似することにより、前記波数エネルギーを計算することを特徴とする画像生成方法。請求項6または7記載の画像生成方法において、前記最適化計算ステップは、前記誤差が最小になるまで前記移流拡散方程式の拡散係数と前記ナビエストークス方程式の粘性係数および密度を変化させることにより、前記予測画像を補正することを特徴とする画像生成方法。請求項5乃至8のいずれか1項に記載の各ステップをコンピュータに実行させることを特徴とする画像生成プログラム。以上のように、従来の画像生成方法の第1の問題点は、画像がぼやけてしまうことである。第2の問題点は、予測過程において、パラメータの最適化基準がなかったことである。第3の問題点は、テクスチャと物理的な性質を予測過程で両立させることができなかったことである。第4の問題点は、数多くの過去の時系列画像を用いる必要があることである。本発明は、上記課題を解決するためになされたもので、予測画像のぼやけを抑制することができ、時系列画像のテクスチャ性と物理特性とに関するエネルギーを保存することができ、モデルパラメータを自動的に最適化することができる画像生成装置、画像生成方法およびプログラムを提供することを目的とする。本発明の画像生成装置は、複数の連続した時系列画像を取り込むデータ入力手段と、前記時系列画像を蓄積するデータ蓄積手段と、前記時系列画像の持つ初期エネルギーを計算する初期エネルギー計算手段と、前記時系列画像から、移流拡散方程式とナビエストークス方程式により前記時系列画像よりも未来の予測画像を生成する予測手段と、前記時系列画像の初期エネルギーと前記予測手段が生成した予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように前記予測画像を補正する最適化計算手段とを備えることを特徴とするものである。また、本発明の画像生成装置の1構成例において、前記初期エネルギー計算手段と前記最適化計算手段とは、画像の動きベクトルから波数エネルギーと運動エネルギーと発散渦エネルギーを計算し、画像の輝度分布からエッジエネルギーを計算し、前記波数エネルギーと前記運動エネルギーと前記発散渦エネルギーと前記エッジエネルギーの平均エネルギーを当該画像の持つエネルギーとすることを特徴とするものである。また、本発明の画像生成装置の1構成例において、前記初期エネルギー計算手段と前記最適化計算手段とは、前記波数エネルギーを計算する際に、連続的な波数とエネルギーとの関係をガンマ分布により近似することにより、前記波数エネルギーを計算することを特徴とするものである。また、本発明の画像生成装置の1構成例において、前記最適化計算手段は、前記誤差が最小になるまで前記移流拡散方程式の拡散係数と前記ナビエストークス方程式の粘性係数および密度を変化させることにより、前記予測画像を補正することを特徴とするものである。また、本発明の画像生成方法は、複数の連続した時系列画像を取り込むデータ入力ステップと、前記時系列画像を蓄積するデータ蓄積ステップと、前記時系列画像の持つ初期エネルギーを計算する初期エネルギー計算ステップと、前記時系列画像から、移流拡散方程式とナビエストークス方程式により前記時系列画像よりも未来の予測画像を生成する予測ステップと、前記時系列画像の初期エネルギーと前記予測ステップで生成した予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように前記予測画像を補正する最適化計算ステップとを備えることを特徴とするものである。また、本発明の画像生成プログラムは、画像生成方法の各ステップをコンピュータに実行させることを特徴とするものである。本発明によれば、複数の連続した時系列画像を取り込み、時系列画像の持つ初期エネルギーを計算し、時系列画像から、移流拡散方程式とナビエストークス方程式により時系列画像よりも未来の予測画像を生成し、時系列画像の初期エネルギーと予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように予測画像を補正することにより、偏微分方程式における物理パラメータを最適化計算により決定することで、不自然な画像のぼやけを抑制することができる。また、本発明では、統計的な乱流エネルギーモデルを導入することで、物理的な性質を保持することができる。すなわち、本発明では、時系列画像のテクスチャ性と物理特性とに関するエネルギーを保存することができる。また、本発明では、2枚の時系列画像からでも連続的に複数の新しい予測画像を生成することができる。また、本発明では、連続的な波数とエネルギーとの関係をガンマ分布により近似することにより、画像から散在的な波数とエネルギーとの関係しか得られない場合でも、連続的な値を得られるようになり、エネルギーの誤差計算をし易くすることができる。本発明は、与えられたテクスチャ性と物理特性とに関するエネルギーが保存されるように、過去から現在までの時系列画像と近未来の予測画像との間で最適化計算を行う方法を提供する。画像の持つ総エネルギーは、運動エネルギーと波数エネルギーと発散渦エネルギーとエッジエネルギーの和である。本発明では、客観的に画像の持つエネルギーを定義し、エネルギー保存則を導入している。本発明では、予測画像においてもエネルギーが保存されるように、移流拡散方程式における拡散係数と、ナビエストークス方程式における粘性係数および密度の3つのパラメータについて自動的に最適化する方法により画質の崩れを抑制することができる。以下、本発明の実施の形態について図面を参照して説明する。図1は本発明の実施の形態に係る画像生成装置の構成を示すブロック図、図2は画像生成装置の処理の流れを説明するフローチャートである。本実施の形態の画像生成装置は、データ入力部100と、データ蓄積部110と、初期エネルギー計算部120と、予測部130と、最適化計算部140と、表示部150とから構成される。画像生成装置の処理を図2を用いて簡単に説明すると、データ入力部100は、カメラから複数の連続した時系列画像を取り込みステップS100、データ蓄積部110は、データ入力部100が取り込んだ時系列画像のデータを保存するステップS101。初期エネルギー計算部120は、過去から現在までの時系列画像の持つ初期エネルギーを計算するステップS102。予測部130は、時系列画像から、移流拡散方程式とナビエストークス方程式により時系列画像よりも未来の予測画像を生成するステップS103。最適化計算部140は、時系列画像の初期エネルギーと予測部130が生成した予測画像のエネルギーとの誤差を計算し、この誤差が最小になるように予測画像を補正するステップS104。表示部150は、予測部130が生成し最適化計算部140が補正した予測画像を表示するステップS105。本実施の形態では、2枚の連続した時系列画像の入力を前提とし、この2枚の時系列画像から未来の1枚以上の時系列画像を予測することを目的とする。予測計算開始前の画像のもつ初期エネルギーを次のように定義する。式1に示すEng1ijnは波数エネルギー、式2に示すEng2ijnは運動エネルギー、式3に示すEng3ijnは発散渦エネルギー、式4に示すEng4ijnはエッジエネルギーである。これらの波数エネルギーEng1ijn、運動エネルギーEng2ijn、発散渦エネルギーEng3ijn、およびエッジエネルギーEng4ijnは、画素i,j毎に、ある時刻nで定義される。kijnは時刻nにおいてx方向の座標がiでy方向の座標がjの画素i,jにおける波数、uijは画素i,jにおける速度動きベクトル、Iijは画素i,jにおける画像輝度分布である。また、wijは画素i,jにおける動きベクトルのx成分、vijは画素i,jにおける動きベクトルのy成分、Ixijは画素i,jにおける画像輝度分布のx成分、Iyijは画素i,jにおける画像輝度分布のy成分である。波数エネルギーEng1ijn、運動エネルギーEng2ijn、および発散渦エネルギーEng3ijnは、画像のもつ動きに関わるエネルギーであり、エッジエネルギーEng4ijnは、画像の見掛けのテクスチャ性に関するエネルギーである。速度uは、オプティカルフロー法やパターンマッチング法などにより計算することができる。オプティカルフロー法については、例えば文献Optical_flow,http://en.wikipedia.org/wiki/Optical_flowに開示されている。本実施の形態では、2枚の連続した時系列画像から速度uを計算するが、1枚の画像しかない場合は、人為的に設定された速度uを活用する。エッジとは、画像輝度分布Iの1次微分の大きさである。画素i,jでの平均エネルギーについては、式10に示すように、4つのエネルギーEng1ijn,Eng2ijn,Eng3ijn,Eng4ijnの和をとり、この和を1枚の画像の全画素数ΩΩはx方向M画素とy方向N画素の積M×Nで割ることにより求めることができる。このとき、式11に示すように、重み付け係数αccは1~4の和を1とする。なお、波数kの推定方法については、後述する図4の説明の中で述べる。図3はステップS102の初期エネルギー計算処理とステップS103の予測処理とステップS104の最適化計算処理のアルゴリズムの概要を示す図である。ステップS200図2のステップS102において、初期エネルギー計算部120は、データ入力部100が取り込んだ2枚の時系列画像から動きベクトルを画素毎に計算し、式1~式4に示した4つの初期エネルギーを画素毎に計算する。図3におけるIijn1は時刻n-1において取り込まれた時系列画像の画素i,jにおける画像輝度分布、Iijnは時刻nにおいて取り込まれた時系列画像の画素i,jにおける画像輝度分布である。ステップS201図2のステップS103,S104において、予測部130は、式12に示す移流拡散方程式と式13に示すナビエストークスNS方程式とを交互に時間発展させていくことにより、時刻nより未来の新しい画像を予測する。ここでは、流体変化を伴う自然現象を対象とした時系列画像を想定し、流体モデルを仮定して、式14に示す非圧縮性条件を課した。式12、式13におけるυは粘性係数、ρは密度、ηは拡散係数である。これら粘性係数υ、密度ρ、拡散係数ηは未知数なので、最適化計算処理で推定することになる。図3に示すハットυ、ハットρ、ハットηは、それぞれ粘性係数υ、密度ρ、拡散係数ηの最適値である。予測部130は、式6に示した移流拡散方程式を時間発展させる。また、NS方程式は、MACMarkerandCell法により非線形に解かれる。MAC法については、例えば文献大宮司久明,数値流体力学大全第13章非圧縮性流れの解法-MAC型解法,http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/Chapter13.pdfに開示されている。こうして、予測部130は、移流拡散方程式により未来の画像輝度分布Iをあらかじめ定められた時間ステップ毎および画素毎に計算すると共に、NS方程式により未来の速度uを時間ステップ毎および画素毎に計算することにより、未来の画像を時間ステップ毎に予測する。図3の例では、時刻n1,n2,,npにおける画像の画素i,jの輝度分布Iと速度uとが計算されている。次に、最適化計算部140は、時刻nの入力時系列画像の初期エネルギーEnと時刻npの予測画像のエネルギーEnpとの誤差|ε|を式15により画素毎に計算する。この計算では、|En-En1|,|En-En2|,,|En-Enp|のように、常に時刻nにおける初期エネルギーと時刻nよりも後の各時刻n1,n2,,npにおけるエネルギーとの誤差|ε|を計算する。粘性係数υ、密度ρ、拡散係数ηの全ての組み合わせについて誤差|ε|を計算すると膨大な計算となることから、予め粘性係数υ、密度ρ、拡散係数ηの範囲を、10.0≦υ≦50.0、0.1≦ρ≦10.0、0.0≦η≦500.0に設定した。NS式については、無次元化変換により、モデルパラメータをレイノルズ数に集約させることができ、この場合は、より少ないモデルパラメータでの最適化が可能である。こうして、最適化計算部140は、予測部130が予測した予測画像の平均エネルギーを初期エネルギーと同様に計算し、初期エネルギー計算部120が計算した入力時系列画像の初期エネルギー平均エネルギーと予測画像の平均エネルギーとの誤差を計算し、誤差が最小になるまで粘性係数υ、密度ρ、拡散係数ηを変化させていくことを、予測の時間ステップ毎および画素毎に行う。図3の例では、時刻n1,n2,,npの予測画像毎および画素毎に、このような処理を行うことになる。最適化計算部140は、時刻n1の各画素について誤差が最小になると、この時刻n1についての最適化計算処理が完了したものと見なして、誤差が最小になったときの粘性係数υ、密度ρ、拡散係数ηで時刻n1の予測画像を確定し、次の時刻n2の各画素について最適化計算処理を行う。図4A~図4Cは予測過程における乱流エネルギーモデルの適用方法を説明する図である。流体力学分野でよく知られているコルモゴロフエネルギーがあり、本実施の形態ではその特性を応用して波数エネルギーを求めた。コルモゴロフエネルギーについては、例えば文献ThayerWatkins,KolmogorovsTheoryofTurbulenceandtheEnergySpectrum,SanJoseStateUniversity,http://www.applet-magic.com/kolmo.htmに開示されている。図4Aに示す分布は、横軸に波数k、縦軸にコルモゴロフエネルギーEkをとったものである。図4Aの400は入力時系列画像の初期エネルギー、401は予測画像のエネルギー、400は初期エネルギーと予測画像のエネルギーとの誤差であるエネルギーロスである。図4Aに示した分布では、低周波数側にピークが現れる。波数kは2πと周波数の積で与えられる。初期エネルギー計算部120は、入力時系列画像の速度uijnを時刻nで得る。そして、初期エネルギー計算部120は、図4Bに示すように周波数変換を高速フーリエ変換FFTにより行い、速度uijnを波数エネルギーEknに変換する。最適化計算部140は、予測過程において、図4Bに示すように時刻npの速度uijnpをFFTにより波数エネルギーEknpに変換する。速度uijnpは、予測過程ではNS方程式より得られたものである。このとき、画像から計算される波数分布は、図4Cに示すように散在する。以後の誤差計算を式15を用いて行うときは、波数kの値が連続でないと、ゼロの値を用いた誤差計算となってしまう。即ち、何らかの補間された波数kの値が必要となる。このため、図4Cに示すようなエネルギー分布について、連続的な波数kとエネルギーとの関係を式16に示すようなガンマ分布により近似している。近似された分布を用いることで、予測過程におけるエネルギーの誤差を計算しやすくすることができる。式16におけるshは形状パラメータ、scはスケールパラメータである。形状パラメータshは整数である。ここで、x>0,sh>0,sc>0となる。Γshはガンマ関数である。次に、式15において、初期と予測過程での誤差エネルギーεに基づいて、補正されるべき速度が推定できる。即ち、式17のように、誤差エネルギーεを逆変換IFFTすることで速度に換算される。また、パラメータの組み合わせにより、初期エネルギーよりも予測過程でエネルギーが増加もしくは減少する。このことを考慮すると、速度の増減を考慮する必要がある。符号関数sgnεは、εの符号に応じて1,0,-1のいずれかの値を返す。即ち、ε>0のときはsgnε=1、ε=0のときはsgnε=0、ε<0のときはsgnε=-1である。図5A~図5Dを用いて、従来技術による予測画像と本実施の形態による予測画像との違いについて説明する。ここでは、画像の例として水面の画像を用いている。図5Aは実画像正解画像を示す図である。図5Bは図5Aの画像よりも過去の画像を入力時系列画像として本実施の形態により予測した画像を示す図である。本実施の形態では、2枚の連続した入力時系列画像を用いて、1枚先の画像1ステップ後の時刻の画像を予測している。一方、図5Cは図5Aの画像よりも過去の画像を入力時系列画像として従来技術により予測した画像を示す図である。ここでは、従来技術として非特許文献1に開示された3次元自己回帰モデルを適用した方法を用いている。この従来技術では、10枚の入力画像を用いて、1枚先の画像を予測している。3次元自己回帰モデルの係数は10枚の画像を入力し、最小二乗法計算により係数を決定した。本実施の形態では、図5Aに示した正解画像と図5Bに示した予測画像との誤差は13となり、従来技術では、正解画像と図5Cに示した予測画像との誤差は30となった。本実施の形態で適用されたフローの推定結果を図5Dに示す。本実施の形態によれば、実際の正解画像に近い光沢感、透明感を持つ予測画像が得られた。一方、従来技術では、予測画像がぼやけてしまうという欠点が現れている。本実施の形態と従来技術との画質の差は明らかである。図6A、図6B、図7A、図7B、図8A、図8Bは本実施の形態による気象レーダ降水パターン画像の予測例を示す図である。ここでは、粘性係数υのみを変化させて予測画像を求めている。速度場は渦度により評価した。図6Aは入力時系列画像を示しており、図6Bは図6Aの入力時系列画像に速度場を描いた画像を示している。図6B中の矢印が速度uを表している。図7Aは図6Aの画像から1時間後の正解画像を示しており、図7Bは正解画像と同時刻の予測画像を示している。このときの粘性係数υは10000である。図8A、図8Bは図6Aの画像から5時間後の予測画像を示している。図8Aの場合の粘性係数υは0.01、図8Bの場合の粘性係数υは1.5である。図9に示す粘性係数υと渦度との関係から分かるとおり、粘性係数υが小さいほど、速度uへの感度が高くなる。粘性係数υが最適化されない場合、台風の渦の速度場が予測の間も残ってしまい、降水パターンは不自然に収縮してしまい、並進運動がなくなる。以上のように、本実施の形態では、偏微分方程式における物理パラメータを最適化計算により決定することで、不自然な画像のぼやけを抑制することができる。また、本実施の形態では、統計的な乱流エネルギーモデルを導入することで、物理的な性質を保持することができる。また、本実施の形態では、2枚の時系列画像からでも連続的に複数の新しい予測画像を生成することができる。本実施の形態で説明した画像生成装置は、例えばCPU、記憶装置およびインタフェースを備えたコンピュータとこれらのハードウェア資源を制御するプログラムによって実現することができる。このようなコンピュータを動作させるための画像生成プログラムは、フレキシブルディスク、CD-ROM、DVD-ROM、メモリカードなどの記録媒体に記録された状態で提供される。CPUは、読み込んだ画像生成プログラムを記憶装置に書き込み、このプログラムに従って本実施の形態で説明した処理を実行する。"