JPWO2016031174A1 - シミュレーション装置、シミュレーション方法、および、記憶媒体 - Google Patents

Abstract

本発明は、非理想的、非連続的または特異性がある観測データを考慮し、広域にわたって高分解能かつ高精度なシミュレーションを行う。シミュレーション装置１００は、システムモデル２１、データ選択処理部３０、複数の観測モデル３１、事後分布生成部４０、事後分布統合部５０および判定部５１を備える。システムモデル２１は、状態ベクトルの時間発展を計算する。データ選択処理部３０は、複数の観測データを選択する。観測モデル３１は、システムモデル２１からの状態ベクトルを観測データとの関係性に基づき変換する。事後分布生成部４０は、観測モデル３１からの状態ベクトルおよび選択された観測データを基に、全観測データに基づく第１の事後分布または欠如した観測データに基づく第２の事後分布を生成する。事後分布統合部５０は、第１および第２の事後分布を統合する。判定部５１は、第２の事後分布と統合後の事後分布との何れを用いるか判定する。

Description

本発明は、実世界で起きている現象や仮説的な状況を数理的にモデル化してコンピュータ上で数値的に計算するシミュレーション技術に関する。

シミュレーションは、実世界で起きている現象や仮説的な状況を数理的にモデル化してコンピュータ上で数値的に計算するものである。数理的にモデル化することにより、時間や空間を自在に設定して計算することが出来る。このようなシミュレーションにより、現実の結果を得ることが困難な状況（例えば、観測が困難な場所の状況）や、未来に起きうることを予測することが出来る。また、計算の条件を意図的に変化させることにより、現実では見ることが難しい状況での特性や振る舞いを調べることができる。そのようなシミュレーション結果は、因果関係の理論的解明や設計、計画などの指標として役立てることができる。

特に、シミュレーションは、現実に得られる観測データの数が少なく、空間的にも時間的にも偏った分布をしている状況について、広域にわたって連続的に状態を把握および理解したい場合に有効である。しかしながら、シミュレーションは、あくまでも現実を数理的に模擬したものに過ぎないため、その精度は、どれだけ現実を理解して忠実に模擬出来ているかに依存する。したがって、上記のように現実の観測データが少なく、また理解が不完全な現象を対象とする場合は、モデルに不完全性が含まれることとなる。さらに、計算は離散的に行われるため、広域にわたって連続的に把握するためには、対象領域を細分化して大量に計算する必要がある。しかし実用上は、許される計算時間や計算機資源に応じて、不完全性のある計算条件を設定せざるを得ない。これらの不完全性は、シミュレーションの精度を悪化させる。

そこで、このような不完全性のある条件下でシミュレーションの精度を向上させる仕組みとして、データ同化が知られている。データ同化は、現実で得られた観測データを数値シミュレーションに取り入れる手法である。同じ数理的なモデルによるシミュレーションを用いても、前述した内在する不完全性、与えられる初期条件、または、境界条件等によって、様々な結果が現れてしまう。データ同化は、その様々なシミュレーション結果の中から、現実で得られた観測データをもっとも上手く説明するものを探しだし、同時にモデルや条件を更新する。

データ同化は、地球科学や海洋学で用いられることが多いが、特に気象学での発展は著しく、すでに日々の気象予報の精度向上に貢献しているだけでなく、新たな手法も提案され続けている。これは、気象に関係する観測技術の向上により、得られる観測データが多種多様になっていること、そして、空間的にも時間的にもデータ量が増加していることにも起因する。

気象予報の中でも、急発達する雷雲や雨等を正確に予測するためには、シミュレーションの空間的、時間的な高精度化が課題である。このような課題に対応する関連技術として、特許文献１に記載された気象予測装置が知られている。この気象予測装置は、地上に多数設置されていて頻繁に観測可能なＧＰＳ受信機による可降水量データ、ドップラーレーダによる風向風速データ、レーダアメダスによる雨量強度データを利用する。そして、この気象予測装置は、リアルタイムまたは準リアルタイムで計測されるこれらのデータを取り込んで、３次元変分法によりデータ同化する。また、このような課題に対応する他の関連技術として、特許文献２に記載された同時化装置およびメッシュ化装置が知られている。この同時化装置は、複数の観測器からのデータが非同期である場合に、これらの観測データを、内挿処理によって時間軸上で再構成することにより、同じ時刻の観測データを示すよう同時化する。また、このメッシュ化装置は、同時化された複数地点の観測データを、対象地域内において、水平空間上で一定距離間隔のメッシュ点（格子点）上に位置するように再配置する。

また、気象や海洋環境の急変により起こる自然災害や人為的な事故などの後に、その土地の状態を広域的に正確に把握することが求められる。このためには、広域的な観測手段を用いて、なおかつ高精度に、その土地の状態を観測および推定することが課題となる。このような課題に対応する関連技術として、特許文献３に記載された手法が知られている。この手法は、シミュレーションを活用した手法ではないが、３時期以上の衛星画像を用いて、対応する土地の状態を表す特徴量を推定する。詳細には、この手法は、現地で直接観測や調査を行うのに十分な時間やコストがない場合や、現地に近づくことが危険な場合に、人工衛星に搭載された合成開口レーダ（Synthetic Aperture Radar：ＳＡＲ）などによる画像データを用いる。そして、この手法は、衛星ＳＡＲ画像を用いて、迅速、広域、そして安全に土地の状態を把握する。その具体例として、特許文献３には、津波発生後の沿岸地域における土壌塩分濃度や排水性を、津波発生前、直後、数か月後の３時期の衛星画像ごとの土壌水分量指標値の変化に基づき推定する例が記載されている。

また、衛星ＳＡＲ画像と、作物の収量予測モデルとを用いて、広範囲の地域における水稲圃場の収量予測を、生育期前半に少ない労力で高精度に行う方法の例が、特許文献４に記載されている。一般に、可視や近赤外領域での太陽光の反射光強度を観測する衛星搭載の光学センサは、雲があると観測が出来ないといったように、天候の影響を大きく受ける。そこで、この手法は、雲の影響を受けないマイクロ波（Ｘバンド：波長３．１ｃｍ）でのＳＡＲ画像データを用いる。そして、この手法は、得られたＳＡＲ画像データと、草丈や茎数といった作物の生育状態を表す量との相関関係から、収量予測式を回帰分析で求めることにより、収量予測を行う。

また、実世界のシミュレーションとは異なるが、観測データの分析に数理モデルが使われる場合がある。例えば、ミリ波レーダなどにより所定領域内に物体が存在するか否かを判定する用途では、観測データのみでの正確な判定は、困難であった。これは、特に物体が歩行者であると、レーダの反射強度が非常に小さい（ＳＮ比:signal to noise ratioが小さい）ためと、歩行者の様々な姿勢変化に伴って反射強度が時々刻々と変化するためである。このような問題に対応する関連技術として、特許文献５に記載された手法が知られている。この手法は、予め対象物の検出位置に対する反射強度の分布から、歩行者信号やノイズ信号の特徴量モデルを作成しておく。そして、この手法は、実際の観測データをこれらモデルと照らし合わせることで、非理想的な観測データから存在・不存在・不明の各種状態を確率的に推定する。さらに、特許文献５には、存在・不存在・不明の各種状態について、複数センサから得られた複数状態の確率を統合する手法についても記載されている。

特開２０１０−６０４４３号公報 特開２０１０−６０４４４号公報 特開２０１３−８４０７２号公報 国際公開第２０１１／１０２５２０号 特開２００７−９３４６１号公報

上述した各関連技術では、不完全性のあるシミュレーションの高精度化や、非理想的な観測データからの状態推定には、取得する観測データの種類や数を増やすことが、解決手段の起点となっている。しかしながら、増やす観測データがすべて理想的であるとは限らないため、非理想的な観測データを考慮した処理方法が必要不可欠である。例えば、特許文献１には、観測データの取得後に異常値を除去する異常値除去装置の記載がある。この異常値除去装置は、同じ時刻の気象モデルによる算出値と比較して一定以上の差がある観測データを、異常値と判定する。しかし、この局所的でモデル計算値という単独の情報に基づく異常値判定では、急激に変化する状態や、特異性のある環境・現象を考慮することができない。また、特許文献５に記載された関連技術は、観測データをモデルと照らし合わせた後に確率的に統合するため、より多くの情報を元に、特異性などの可能性も含めて各種状態を推定可能である。しかしながら、この関連技術は、確率の統合方法として一般的な重み付け平均を主として用いており、対象となるのが同種類・同次元の観測データに限られるという問題がある。また、この関連技術は、統合結果の判定や修正を可能とする仕組みについては考慮していないという問題もある。

また、不完全性のあるシミュレーションの高精度化や、非理想的な観測データからの状態推定には、取得する観測データの種類や数を増やすことに加え、得られた観測データを予め加工することが考えられる。特許文献２に記載の同時化装置やメッシュ化装置は、得られた観測データのみに基づいて統計的に時間および空間補間を行うが、観測データが連続的で相関が高いことが前提である。また特許文献３には、異なる３時期の衛星画像の統合により土地の状態を推定する手法の記載がある。しかしながら、特許文献３は、衛星データが最適な時期に取得できていなかった場合や、土壌表面に溜まった水の影響などによりデータが特異的な値を示した場合の対処について言及していない。このように、特許文献２および３に記載された関連技術は、何らかの非連続性や特異性がある観測データについては考慮していない。また特許文献４に記載の関連技術は、理想的な観測データの取得時間間隔が気候の影響により不安定な光学センサデータを用いずに、データ取得の確実性を優先してマイクロ波センサデータを用いている。そして、この関連技術は、そのようにして得られたデータから回帰的（帰納的）に収量予測モデルを作成するため、データに非連続性や特異性がある場合に不適切なモデルを生成する、という懸念がある。

次に、特許文献１および２に記載の、典型的なデータ同化を使ったシミュレーションにおける課題について、図面を参照して説明する。図１５に、典型的なデータ同化を使ったシミュレーションを説明する模式図を示す。図１５において、横軸は時間変化を表し、縦軸は観測値を表す。なお、観測値は、実際にセンサなどによって測定される値である。また、状態変数は、シミュレーションモデル上で時間発展を計算するための変数である。ここで、観測値は、観測頻度やセンサ精度の影響などを受けるため、対象とする量の真の値であるとは限らない。そこで、図１５では、この未知である真の値を仮想的に破線で表している。また、観測値は、状態変数と観測モデルとによって関係付けられるため、状態変数の時間変化を計算することで、観測値のシミュレーションが可能である。しかし前述したように、モデルに内在する不完全性、与えられる初期条件、また、境界条件に不確実性がある場合は、シミュレーションが真の値を再現することは難しい。そこでデータ同化は、意図的に擾乱を加えて確率的にシミュレーションをして、得られた様々なシミュレーション結果の中から、現実で得られた観測データをもっとも上手く説明するものを探しだす。図１５では、複数の計算値（シミュレーション）を表す実線で、このプロセスを模式的に表している。各シミュレーションは、観測値によって修正された値を次ステップの初期値としてシミュレーションを継続する。図１５から明らかなとおり、次の観測値が得られるまでの時間が長いと、その間のシミュレーション値を修正することは出来ないので、真の値からの誤差が積算されて大きくなる。また、得られた観測値が不適切、もしくは誤差を多く含む場合は、そのデータによって修正された値が次ステップの初期値となるので、次の確かな値が得られるまでの間のシミュレーション値は誤差が大きくなる。

本発明は、上記の問題に鑑みてなされたものである。すなわち、本発明は、非理想的な観測データや、非連続性または特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行う技術を提供することを目的としている。

上記目的を達成するために、本発明のシミュレーション装置は、シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとを入力として取得する入力手段と、前記初期状態およびパラメータを基に、前記状態ベクトルの時間発展をシミュレーションするシステムモデルと、前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択するデータ選択処理手段と、選択された複数の観測データにそれぞれ対応する複数の観測モデルであって、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換して出力する複数の観測モデルと、前記複数の観測モデルから出力される状態ベクトルおよび前記データ選択処理手段で選択された観測データに基づいて、前記状態ベクトルの事後分布を生成し、前記データ選択処理手段で選択された全ての観測データに基づく事後分布を第１の事後分布として出力し、１つ以上欠如した観測データに基づく事後分布を第２の事後分布として出力する事後分布生成手段と、前記第１の事後分布と前記第２の事後分布とを統合する事後分布統合手段と、前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定する判定手段と、前記判定手段によって決定された事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力するとともに、該状態ベクトルの時系列を出力する出力手段と、を備える。

また、本発明のシミュレーション方法は、シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとが入力されると、前記初期状態およびパラメータを基に、システムモデルを用いて前記状態ベクトルの時間発展をシミュレーションし、前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択し、選択された複数の観測データにそれぞれ対応する複数の観測モデルを用いて、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換し、前記複数の観測モデルから出力される状態ベクトルおよび選択された前記観測データに基づいて、前記状態ベクトルの事後分布を生成し、選択された全ての観測データに基づく第１の事後分布と、１つ以上欠如した観測データに基づく第２の事後分布とを統合し、前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定し、決定した事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力し、決定した事後分布および前記第１の事後分布からなる状態ベクトルの時系列を出力する。

また、本発明の記憶媒体は、シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとを入力として取得する入力ステップと、システムモデルを用いて、前記初期状態およびパラメータを基に、前記状態ベクトルの時間発展をシミュレーションするシステムモデル計算ステップと、前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択するデータ選択処理ステップと、選択された複数の観測データにそれぞれ対応する複数の観測モデルを用いて、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換して出力する観測モデル計算ステップと、前記複数の観測モデルから出力される状態ベクトルおよび前記データ選択処理ステップで選択された観測データに基づいて、前記状態ベクトルの事後分布を生成し、前記データ選択処理ステップで選択された全ての観測データに基づく事後分布を第１の事後分布として出力し、１つ以上欠如した観測データに基づく事後分布を第２の事後分布として出力する事後分布生成ステップと、前記第１の事後分布と前記第２の事後分布とを統合する事後分布統合ステップと、前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定する判定ステップと、前記判定ステップで決定された事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力するとともに、該状態ベクトルの時系列を出力する出力ステップと、をコンピュータ装置に実行させるコンピュータ・プログラムを記憶している。

本発明は、非理想的な観測データや、非連続性または特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行う技術を提供することができる。

本発明の第１の実施の形態としてのシミュレーション装置の構成を示すブロック図である。 本発明の第１の実施の形態としてのシミュレーション装置のハードウェア構成の一例を示すブロック図である。 本発明の第１の実施の形態としてのシミュレーション装置における状態ベクトルと観測データとの関係を模式的に表す図である。 本発明の第１の実施の形態としてのシミュレーション装置のシミュレーション開始時の動作を説明するフローチャートである。 本発明の第１の実施の形態としてのシミュレーション装置のシミュレーション動作を説明するフローチャートである。 本発明の第１の実施の形態の効果を模式的に説明する図である。 本発明の第２の実施の形態としてのシミュレーション装置の構成を示すブロック図である。 本発明の第２の実施の形態における観測データの時系列変化および計算格子空間の関係を模式的に表す図である。 本発明の第２の実施の形態としてのシミュレーション装置の事後分布統合において、バリオグラムを推定した結果の一例を示す図である。 本発明の第２の実施の形態としてのシミュレーション装置の事後分布統合動作を説明するフローチャートである。 本発明の第３の実施の形態としてのシミュレーション装置の構成を示すブロック図である。 本発明の第３の実施の形態における観測データの時系列変化および計算格子空間の関係を模式的に表す図である。 本発明の第４の実施の形態としてのシミュレーション装置の構成を示すブロック図である。 本発明の第４の実施の形態における観測データおよび計算格子空間の関係を模式的に表す図である。 関連技術のデータ同化を用いたシミュレーションを模式的に表す図である。

以下、図面を参照して本発明の実施の形態について説明する。

（第１の実施の形態）
本発明の第１の実施の形態としてのシミュレーション装置１００について説明する。シミュレーション装置１００は、物理法則に基づいた連続時間・空間の偏微分方程式を解いて時間発展を追うシミュレーションに適用可能である。そのような偏微分方程式には、例えば、運動を記述する運動方程式、流体を記述するナビエーストークス方程式、熱変化を記述する熱力学方程式、津波を記述する浅水波方程式などがある。また、シミュレーション装置１００は、有限要素法を用いたシミュレーションにも適用可能である。なお、本実施の形態では、シミュレーション対象となる系は、時間変化を追う状態ベクトルが、実際の観測データと何らかの関係式で結ばれる系、すなわち、シミュレーション結果と観測データとが比較可能な系であるものとする。

まず、シミュレーション装置１００の構成を図１に示す。図１において、シミュレーション装置１００は、システムモデル２１と、データ選択処理部３０と、ｍ個の観測モデル３１（第１の観測モデル３１−１〜第ｍの観測モデル３１−ｍ）と、事後分布生成部４０と、事後分布統合部５０と、判定部５１とを備える。また、シミュレーション装置１００は、入力部１０と、出力部６０とを備える。なお、ｍは、２以上かつＭ以下の整数であり、Ｍは、２以上の整数である。また、図１において、シミュレーション装置１００は、各機能ブロック間で入出力されるデータの格納領域として、事前分布記憶部２２と、第１の事後分布記憶部４１ａと、第２の事後分布記憶部４１ｂと、統合後事後分布記憶部５２とを含む。なお、事前分布記憶部２２は、本発明のシステムモデルの一部の一実施形態を構成する。また、第１の事後分布記憶部４１ａおよび第２の事後分布記憶部４１ｂは、本発明の事後分布生成部の一部の一実施形態を構成する。また、統合後事後分布記憶部５２は、本発明の事後分布統合部の一部の一実施形態を構成する。

ここで、シミュレーション装置１００のハードウェア構成の一例を図２に示す。図２において、シミュレーション装置１００は、ＣＰＵ（Central Processing Unit）１００１と、ＲＡＭ（Random Access Memory）１００２と、ＲＯＭ（Read Only Memory）１００３と、ハードディスク等の記憶装置１００４と、入力装置１００５と、出力装置１００６とによって構成される。この場合、入力部１０は、入力装置１００５と、ＲＯＭ１００３または記憶装置１００４に記憶されたコンピュータ・プログラムをＲＡＭ１００２に読み込んで実行するＣＰＵ１００１とによって構成される。また、システムモデル２１、データ選択処理部３０、観測モデル３１、事後分布生成部４０、事後分布統合部５０および判定部５１は、次のように構成される。すなわち、これらの機能ブロックは、ＲＯＭ１００３または記憶装置１００４に記憶されたコンピュータ・プログラムをＲＡＭ１００２に読み込んで実行するＣＰＵ１００１によって構成される。また、出力部６０は、出力装置１００６と、ＲＯＭ１００３または記憶装置１００４に記憶されたコンピュータ・プログラムをＲＡＭ１００２に読み込んで実行するＣＰＵ１００１とによって構成される。なお、シミュレーション装置１００およびその各機能ブロックのハードウェア構成は、上述の構成に限定されない。

次に、シミュレーション装置１００の各機能ブロックの詳細について説明する。

まず、入力部１０について説明する。入力部１０は、シミュレーションで扱う観測領域の状態ベクトルの初期状態およびパラメータと、Ｍ種類の観測データ（第１〜第Ｍの観測データ）とを取得する。ここで、Ｍ種類の各観測データは、センサ等による観測値である。また、Ｍ種類の各観測データは、他の一部または全部の観測データに対して次元が異なっていてもよいし、同一であってもよい。例えば、入力部１０は、記憶装置１００４に記憶されたこれらの情報を取得してもよい。また、入力部１０は、これらの情報の記憶装置１００４における格納位置情報を、入力装置１００５を介して取得することにより、これらの情報を取得してもよい。

次に、システムモデル２１について説明する。システムモデル２１は、入力部１０で取得された初期状態およびパラメータを基に、状態ベクトルの時間発展をシミュレーションする。ここで、対象とする現実の現象の時間発展は、連続時間・空間の偏微分方程式で表されるが、シミュレーション実施のためには、シミュレーションの対象とする領域を時間・空間上で離散化する必要がある。シミュレーション装置１００では、観測領域における現実の現象の時間発展を追うため、状態変数の組からなる状態ベクトルを用いる。なお、状態変数の数は、シミュレーションの目的に応じて定めればよく、いくつであっても構わない。本実施の形態では、状態変数が２つである例を中心に説明する。この場合、２つの状態変数Ｕ，Ｖは、ベクトルξを使うと、ξ＝（Ｕ，Ｖ）とも表される。

時間上の離散化は、ある時刻ｔでの状態変数ξｔから１ステップ進めて時刻ｔ＋１での状態変数ξｔ＋１を求めることにより実現される。なお、以降の説明において、時刻とは、シミュレーションにおける１ステップを示し、例えば時刻ｔ−１とは、時刻ｔより１ステップ前を意味する。また、以降、シミュレーションにおけるステップを、時間ステップとも記載する。

また、空間上の離散化は、２次元空間を格子状に区切ることを仮定し、ある基準点から数えてｋ番目の格子点で定義される時刻ｔでの状態変数を（ξ_ｋ）ｔと表すことにより実現される。ただしξ_ｋ＝（Ｕ_ｋ,Ｖ_ｋ）で表される。つまり、状態ベクトルは、シミュレーションの対象とする領域内で離散化された格子点ごとの状態変数からなる。対象とする領域を表す最後の格子点番号をＬとすると、時刻ｔでの状態変数の組は、次式（１）に示す状態ベクトルＸ_ｔで表される。

・・・（１）
なお、式中の記号Ｔは、転置を意味する。この状態ベクトルの次元は、１格子点あたりの状態変数の数と格子点数Ｌとの積となる。（１）式の場合は、２Ｌ次元となる。

また、システムモデル２１は、離散化された時間・空間での、時刻ｔ−１での状態ベクトルＸ_ｔ−１から、時刻ｔでの状態ベクトルＸ_ｔへの更新操作を行う。更新操作の写像をｆとすると、システムモデル２１は、

・・・（２）
という関係式で記述される。ここで、θは、モデルの計算に必要な各種パラメータベクトルを表す。また、Ｖ_ｔは、時刻ｔでのシステムノイズを表す。システムノイズＶ_ｔは、モデル内の不完全性の効果を数値的に表すために、状態ベクトルに作用を及ぼす確率的な駆動項として導入される。また、写像ｆは、対象とする現象によって線形、または非線形であっても良い。また、（２）式から明らかなように、時刻ｔでの状態ベクトルＸ_ｔが時刻ｔ−１の状態ベクトルＸ_ｔ−１で陽に定義されている必要はない。すなわち、本実施の形態のシステムモデル２１は、時刻ｔ−１での状態ベクトルＸ_ｔ−１を入力として、時刻ｔでの状態ベクトルＸ_ｔを計算できれば良い。

以下に、システムモデル２１による状態ベクトルの時刻ｔ−１から時刻ｔへの更新操作について詳細に説明する。ここでは、まず、不完全性のあるシミュレーションを取り扱うための、アンサンブル近似について説明しておく。以降、システムモデル２１の状態ベクトルＸ_ｔおよびシステムノイズＶ_ｔを、写像ｆの不完全性や入力されるパラメータθの不完全性を反映して、確定値Ｘ_ｔ、Ｖ_ｔではなく、確率分布ｐ（Ｘ_ｔ）、ｐ（Ｖ_ｔ）として扱う。これら確率分布を、それぞれＮ個のアンサンブルｉの集合、

・・・（３）

によって近似することを、アンサンブル近似という。したがって、実際のシステムモデル２１は、この各アンサンブルｉの時間発展

・・・（４）

をすべてのアンサンブルについて求める。これによって、時刻ｔでの状態ベクトルＸ_ｔの確率分布ｐ（Ｘ_ｔ）は、Ｎ個のアンサンブル

・・・（５）
によって近似される。この（４）式で表されるアンサンブルの計算は、アンサンブルごとに独立している点が特徴である。したがって、システムモデル２１は、Ｎ回繰り返して計算をしても良いし、Ｎ個の並列計算機によって１度に計算を行ってもよく、計算リソースに応じて柔軟に計算方法を変えることが出来る。なお、以降では、状態ベクトルＸ_ｔの確率分布をｐ（Ｘ_ｔ）と示し、事前分布とも呼ぶ。

このように、システムモデル２１をアンサンブルごとに独立に用いることで、時刻ｔ−１での状態ベクトルの確率分布ｐ（Ｘ_ｔ−１）から、時刻ｔでの状態ベクトルの確率分布ｐ（Ｘ_ｔ）を算出することができる。また、システムモデル２１は、算出した確率分布ｐ（Ｘ_ｔ）を事前分布として、後述のｍ個の観測モデル３１に出力する。例えば、システムモデル２１は、算出した事前分布ｐ（Ｘ_ｔ）を、ｍ個の観測モデル３１によって読み出し可能な事前分布記憶部２２等に格納してもよい。

次に、データ選択処理部３０および観測モデル３１について説明する。データ選択処理部３０は、状態ベクトルに関連する情報に基づいて、第１〜第Ｍの観測データのうち、利用するｍ種の複数の観測データを選択する。そして、データ選択処理部３０は、選択したｍ種の観測データを、後述の事後分布生成部４０に出力する。

ここで、本実施の形態では、状態ベクトルに関連する情報は、システムモデル２１からデータ選択処理部３０に制御信号ＣＴＬ０として入力されるものとする。制御信号ＣＴＬ０は、例えば、状態ベクトルＸ_ｔの次元数や、状態変数に関するその他の情報などを含んでいてもよい。この制御信号ＣＴＬ０に含まれる情報を元に、データ選択処理部３０は、利用するｍ種類、時刻ｔでの観測データＯＢＳ_１〜ＯＢＳ_ｍを選択する。そして、データ選択処理部３０は、選択した観測データＯＢＳ_１〜ＯＢＳ_ｍを、事後分布生成部４０へ出力する。

また、データ選択処理部３０は、システムモデル２１で設定された状態ベクトルおよび各観測データに関する情報（例えば、物理量や次元）を比較することにより、選択したｍ種の観測データにそれぞれ対応するｍ個の観測モデル３１を生成してもよい。例えば、データ選択処理部３０は、状態ベクトルＸ_ｔと観測データＯＢＳ_１〜ＯＢＳ_ｍとの間を関係付けるためのｍ個の制御信号ＣＴＬ１〜ＣＴＬｍにより、ｍ個の観測モデル３１を生成してもよい。

ここで、観測モデル３１の生成とは、状態ベクトルＸ_ｔと観測データＯＢＳ_１〜ＯＢＳ_ｍとの間を関係付ける観測モデル式を生成することである。これらの状態ベクトルと観測データとの関係を、模式的に図３に示す。このような関係を表す観測モデル式は、例えば、次式（６）で表される。

・・・（６）
この場合、データ選択処理部３０は、（６）式における写像ｈ_１〜ｈ_ｍおよび各観測データで考慮すべきノイズ量ｗ_１〜ｗ_ｍの情報を、制御信号ＣＴＬ１〜ＣＴＬｍとして、観測モデル３１−１〜３１−ｍにそれぞれ出力すればよい。これにより、ｍ個の観測モデル３１が生成される。

ここで、（６）式におけるノイズ量ｗ_１〜ｗ_ｍおよび観測データＯＢＳ_１〜ＯＢＳ_ｍは、もし理想的に全ての格子点１〜Ｌで観測データが得られたとすると、それぞれ、Ｌ次元の列ベクトルとなる。なお、データ選択処理部３０は、各観測データの分散値やノイズ量を元に、その観測データに対応する観測モデル３１のノイズ量を設定してもよい。

また、（６）式においてＥ_１〜Ｅ_ｍは、システムモデル２１の格子点１〜Ｌと、実際に観測データが得られた観測点の解像度とを対応付ける行列である。例えば、状態ベクトルＸ_ｔの各格子点ｋにおける状態変数ξ_ｋ＝（Ｕ_ｋ，Ｖ_ｋ）が２つであるとすると、Ｅ_１〜Ｅ_ｍは、最大２Ｌ×２Ｌ次元の行列となる。

なお、一般に格子点ｋにおける状態変数ξ_ｋ＝（Ｕ_ｋ，Ｖ_ｋ）において、Ｕ_ｋおよびＶ_ｋは、それぞれ異なる物理量である。そのため、Ｕ_ｋおよびＶ_ｋと観測データとの関係を、同一の観測モデル式の写像ｈ_１〜ｈ_ｍで関係付けることは出来ない。以降では、状態変数のうち、Ｕ_ｋ（ｋ＝１〜Ｌ）に対する構成を例として説明する。この場合、上記の行列Ｅ_１〜Ｅ_ｍは、２Ｌ×Ｌ次元となる。例えば、システムモデル２１の格子点１〜Ｌについて、観測データＯＢＳ_１がすべて得られたと仮定する。この場合、Ｅ_１は２Ｌ×Ｌ次元となり、次式（７）の行列で表される。

・・・（７）
この行列は、ｊ（ｊは１以上Ｌ以下の整数）行目の｛１＋２（ｊ−１）｝列目の要素のみ１を持つ。

また、図３に示す観測データＯＢＳ_２のように、格子点１〜Ｌのうち一部の点でデータが観測されていない場合を考える。＊は、観測された格子点（観測点）を表す。この場合、行列Ｅ_２は、（７）式において、観測されなかった格子点に対応する行の要素を０に変更した行列で表される。このとき、（６）式では、左辺においてｈ_２が作用する次元がＬより小さくなり、観測データＯＢＳ_２と同一次元となる。また、（７）式は、例えば格子点ｋにおける状態変数が（Ｕ_ｋ，Ｖ_ｋ，…Ｚ_ｋ）のようにＪ個ある場合には、ｊ行目の｛１＋Ｊ（ｊ−１）｝列目の要素のみ１を持つこととなる。このため、行列Ｅ_１〜Ｅ_ｍは、状態変数の数がいくつであっても表現可能である。また、（６）式における写像ｈ_１〜ｈ_ｍは、状態変数および観測データの関係に依って、線形または非線形であっても良い。したがって（６）式で表される演算は、例えば第ｊの観測モデル３１−ｊの場合、（４）式で計算される時刻ｔの状態ベクトルＸ_ｔが入力されると、

・・・（８）
が出力される、と表すことができる。したがって、全ｍ種の観測モデル３１は、状態ベクトルＸ_ｔに対してそれぞれ（８）式の演算を行い、全ｍ種の変換された状態ベクトルＸ_ｔを、事後分布生成部４０へ出力する。なお（２）式または（４）式、および、（８）式の組は、状態空間モデルと呼ばれる。

なお、上述のＥ_１〜Ｅ_ｍの例は、システムモデル２１の格子点（Ｌ次元）と観測モデルの格子点（Ｌ次元）とが一致する場合を想定しているが、実用上では、一致しない場合も想定される。この場合、行列Ｅ_１〜Ｅ_ｍにおいて、実際に観測データが得られた観測点が、例えば近傍の格子点の値の加重平均や重み付け平均となるように、各要素の値を変更すればよい。このように、上述のＥ_１〜Ｅ_ｍは、状態変数の格子点と、複数の観測データの観測点の解像度とを、１対１、加重平均、または、重み付き和などにより、観測データ毎に関係付けることを表している。

このように、ｍ個の観測モデル３１は、データ選択処理部３０によって選択されたｍ個の観測データにそれぞれ対応している。そして、各観測モデル３１は、システムモデル２１から出力される状態ベクトルを、観測データおよび状態ベクトルの関係性を表す（８）式に基づき、所定の状態ベクトルに変換する。また、各観測モデル３１は、変換した状態ベクトルを、事後分布生成部４０に対して出力する。変換した状態ベクトルは、時刻ｔにおけるｍ種の変換された状態ベクトルＸ_ｔの事前分布を有する。

事後分布生成部４０は、ｍ個の観測モデル３１から出力される状態ベクトルおよびデータ選択処理部３０で選択された観測データに基づいて、状態ベクトルの事後分布を生成する。そして、事後分布生成部４０は、生成した事後分布のうち、データ選択処理部３０で選択されたｍ種全ての観測データに基づく事後分布を第１の事後分布とする。また、事後分布生成部４０は、データ選択処理部３０で選択されたｍ種の観測データのうち１つ以上欠如した観測データに基づく事後分布を第２の事後分布とする。また、事後分布生成部４０は、生成した第１の事後分布および第２の事後分布を、事後分布統合部５０等に出力する。例えば、事後分布生成部４０は、生成した第１の事後分布を、事後分布統合部５０等によって読み出し可能な第１の事後分布記憶部４１ａに格納してもよい。また、事後分布生成部４０は、生成した第２の事後分布を、事後分布統合部５０等によって読み出し可能な第２の事後分布記憶部４１ｂに格納してもよい。また、事後分布生成部４０は、生成した第１の事後分布を、システムモデル２１および出力部６０に対しても出力する。この場合、例えば、事後分布生成部４０は、生成した第１の事後分布を、システムモデル２１および出力部６０によって読み出し可能な統合後事後分布記憶部５２に格納してもよい。

ここで、事後分布生成部４０による事後分布の生成処理について詳細に説明する。事後分布生成部４０には、時刻ｔにおけるｍ種の変換されたＸ_ｔの事前分布と、観測データＯＢＳ_１〜ＯＢＳ_ｍとが入力される。ここで、一般に、事前分布ｐ（ｘ）と観測データの分布ｐ（ｙ）が入力されたときの事後分布ｐ（ｘ｜ｙ）は、ベイズの定理より、

・・・（９）
と表される。ここで、分子のｐ（ｙ｜ｘ）は、状態変数ｘの観測値ｙへの当てはまり度合の指標である尤度と呼ばれる。尤度ｐ（ｙ｜ｘ）は、観測モデル３１が（８）式で表されるように写像ｈとノイズ量ｗに分離できる場合は、

・・・（１０）
によって計算される量を用いることができる。ただし、ｒは、ノイズ量ｗの密度関数である。なお、（１０）式において、右辺を、ｙとｈ（ｘ）を用いた関数ＬＨとして再定義した。さらに、ｍ種類の観測値ｙ＝｛ｙ_１，ｙ_２，・・・ｙ_ｍ｝が得られた場合の尤度ｐ（ｙ_１：ｍ｜ｘ）は、乗法定理を再帰的に使うことで、

・・・（１１）
と積の形で表される。ここで、（１１）式における１項目のｐ（ｙ_１｜ｙ_１：０，ｘ）は、観測データがないときのｙ_１の確率であり、すなわち、ｙ_１が得られたときのｘの尤度ｐ（ｙ_１｜ｘ）となる。また、第２項のｐ（ｙ_２｜ｙ_１：１，ｘ）は、ｙ_１が得られたときのｙ_２の確率である。ただし、それぞれの観測データは、異なるセンサなどで取得されており、ｙ_１とｙ_２の同時分布が存在しない。このため、第２項は、結果的に、ｙ_２が得られたときのｘの尤度ｐ（ｙ_２｜ｘ）となる。したがって、この場合の（９）式で表される事後分布は、

・・・（１２）
と表される。ここで分母のＺは規格化定数とする。この関係を用いると、格子点ｋにおける状態変数Ｕ_ｋの事後分布は、事前分布をｐ（Ｕ_ｋ）と表すと、観測データｙがＯＢＳ_１〜ＯＢＳ_ｍとｍ種得られていることから、

・・・（１３）
と表される。分子は、（１２）式で表されるとおり、各観測データでの尤度の積と、事前分布ｐ（Ｕ_ｋ）との積である。さらに、それぞれの尤度は（１０）式で表されるので、（１３）式の事後分布は、

・・・（１４）
と表される。このように、事後分布生成部４０は、格子点ｋにおける状態変数Ｕｋの事後分布を、ｍ個の観測データＯＢＳ_１〜ＯＢＳ_ｍと写像ｈ_１〜ｈ_ｍとにより計算されるｍ種の尤度ＬＨと、事前分布ｐ（Ｕｋ）とから算出する。同様にして、事後分布生成部４０は、格子点１〜Ｌの全てについて、（１３）式すなわち（１４）式により事後分布を計算する。

ただし、事後分布生成部４０は、ｍ種の観測データのうち１点以上に欠損がある格子点については、（１３）式の代わりに次式（１５）を用いる。例えば、図３に示したＯＢＳ_２のように、一部の格子点でしか観測データが得られていない場合がある。この場合、事後分布生成部４０は、（１３）式の分子に表されるような全ての観測データに対する尤度の積を算出できない。例えば、格子点ｋ’において、ｍ−１種の観測データしか得られないことを想定すると、事後分布は、

・・・（１５）
すなわち、

・・・（１６）
と表され、分子を構成する尤度の数がｍ−１に減少する。なお、（１５）式および（１６）式において、“ｍ−１”は、ｍ種のうち少なくとも１種の観測データが得られていないことを表しており、得られていない（欠損している）観測データの数を１つに限定するものではない。

このように、事後分布生成部４０は、各状態変数について格子点毎に、事後分布を生成する。以降、各状態変数についての格子点毎の事後分布を、状態変数・格子点毎の事後分布、とも記載する。そして、事後分布生成部４０は、全ての観測データに基づき（１３）式を用いて算出した事後分布を第１の事後分布として出力する。また、事後分布生成部４０は、ｍ種のうち少なくとも１種が欠如した観測データから（１５）式を用いて算出した事後分布を、第２の事後分布として出力する。

ここで、ある事前分布ｐ（ｘ）が平均μ０、分散Ｖ_ｐｒｉｏの正規分布に従うと仮定し、ｎ個の観測値ｙ_１，ｙ_２，．．．ｙ_ｎも平均μ、分散Ｖの正規分布に従うと仮定する。この場合、ベイズの定理（９）式で計算される事後分布ｐ（ｘ｜ｙ）も正規分布となり、その分散Ｖ_ｐｏｓｔは、

・・・（１７）
と表される。これより、事後分布の計算に用いられる観測値の数が多いほど、分散は減少する、すなわち事後分布の確度が向上することがわかる。

ここで、第１および第２の事後分布は、それぞれ正規分布とは限らないが、第２の事後分布は、第１の事後分布に比べて取り込んだ観測データの数が少ない。そこで、（１３）式の第１の事後分布ｐ（Ｕ_ｋ｜ＯＢＳ_１：ｍ）の分散をＶａｒ（ｐ（Ｕ_ｋ｜ＯＢＳ_１：ｍ））と表し、（１５）式の第２の事後分布ｐ（Ｕ_ｋ’｜ＯＢＳ_{１：ｍ−１}）の分散をＶａｒ（ｐ（Ｕ_ｋ’｜ＯＢＳ_{１：ｍ−１}））と表す。すると、次式（１８）の不等式が成り立つ。

・・・（１８）
次に、事後分布統合部５０について説明する。事後分布統合部５０は、第１の事後分布と第２の事後分布とを統合する。詳細には、事後分布統合部５０は、第２の事後分布が算出された状態変数・格子点のそれぞれについて、第１の事後分布および第２の事後分布を統合して新たな事後分布を算出する。また、事後分布統合部５０は、統合後の新たな事後分布を判定部５１に出力する。本実施の形態では、事後分布は、アンサンブルの集合によって近似されているため、事後分布統合部５０は、第１の事後分布を近似するアンサンブルおよび第２の事後分布を近似するアンサンブルを、所定の割合で重ねあわせることにより統合すればよい。

具体的には、事後分布統合部５０は、第１の事後分布記憶部４１ａおよび第２の事後分布記憶部４１ｂより、前述した第１の事後分布および第２の事後分布を入力として取得する。ここで、（１８）式の関係より、第２の事後分布の１つｐ（Ｕ_ｋ’｜ＯＢＳ_{１：ｍ−１}）は、少なくとも第１の事後分布の１つｐ（Ｕ_ｋ｜ＯＢＳ_１：ｍ）よりも分散が大きい（すなわち確度が低い）。そこで、事後分布統合部５０は、第２の事後分布が算出された各状態変数・格子点について、第１の事後分布および他の第２の事後分布を統合して新たな事後分布を算出する。例えば、格子点ｊについて第２の事後分布ｐ（Ｕ_ｊ｜ＯＢＳ_{１：ｍ−１}）が算出されていたとする。この場合、事後分布統合部５０は、格子点ｊについて、新たな事後分布ｐ’（Ｕ_ｊ｜ＯＢＳ_１：ｍ）を、ｇをある関数として、次式（１９）により算出する。

・・・（１９）
ここで、πは、関数ｇを決定付けるパラメータセットを表す。また、ｋは、第１の事後分布が生成された格子点を表す。また、ｉは、第２の事後分布が生成された他の格子点を表す。なお、ｉ≠ｊである。以降、（１９）式における確率分布ｐ’のダッシュ（’）は、事後分布統合部５０での統合後の確率分布であることを示す。事後分布統合部５０は、このようにして新たに算出された事後分布ｐ’（Ｕ_ｊ｜ＯＢＳ_１：ｍ）と、元の第２の事後分布ｐ（Ｕ_ｊ｜ＯＢＳ_{１：ｍ−１}）とを、判定部５１に出力する。

次に、判定部５１について説明する。判定部５１は、第２の事後分布と、統合後の事後分布とのうち、いずれを用いるかを決定する。詳細には、判定部５１は、第２の事後分布が生成された状態変数・格子点について、元の第２の事後分布と、統合後の事後分布とのうち、いずれを事後分布として用いるかを決定する。具体的には、判定部５１は、決定した事後分布を、統合後事後分布記憶部５２に保存すればよい。統合後事後分布記憶部５２には、前述のように、第１の事後分布が保存されている。これにより、統合後事後分布記憶部５２には、各状態変数・格子点について、第１の事後分布、または、決定された事後分布が保存されることになる。

例えば、判定部５１は、第２の事後分布および統合後の事後分布の各分散値に基づいて、いずれを用いるかを決定してもよい。具体的には、判定部５１には、事後分布統合部５０によって新たに算出された統合後の事後分布ｐ’（Ｕ_ｊ｜ＯＢＳ_１：ｍ）と、元の第２の事後分布ｐ（Ｕ_ｊ｜ＯＢＳ_{１：ｍ−１}）とが入力される。これらの事後分布は、どちらも格子点ｊにおける事後分布である。例えば、判定部５１は、式（１８）と同様に、これらの事後分布の分散を計算して比較してもよい。この場合、統合後の事後分布ｐ’（Ｕ_ｊ｜ＯＢＳ_１：ｍ）の分散の方が小さければ、判定部５１は、統合後の事後分布を選択して出力する。また、判定部５１は、選択した事後分布を、統合後事後分布記憶部５２に保存する。

一方、統合後の事後分布ｐ’（Ｕ_ｊ｜ＯＢＳ_１：ｍ）の分散の方が大きい場合、判定部５１は、統合後の事後分布の分散が元の第２の事後分布の分散より小さくなるまで、（１９）式における関数ｇのパラメータπを変化させて繰り返し計算を行ってもよい。例えば、関数ｇが重み付け平均の場合、判定部５１は、その重み付けの係数を変化させてもよい。また、判定部５１は、パラメータπに、ある事前分布ｐ（π_ｐｒｉｏ）を仮定し、（４）式の分散値を観測値としてベイズの定理（９）式を用いることにより、分散が最小となるパラメータπの事後分布ｐ（π_ｐｏｓｔ）を求めても良い。パラメータπを変化させることにより、統合後の事後分布の分散が元の第２の事後分布の分散より小さくなった場合、判定部５１は、その統合後の事後分布を選択して統合後事後分布記憶部５２に保存する。パラメータπを変化させても分散が小さくならない場合、判定部５１は、元の第２の事後分布を選択し、統合後事後分布記憶部５２に保存する。

このようにして、統合後事後分布記憶部５２では、第１の事後分布と、判定部５１によって選択された統合後の事後分布または第２の事後分布とによって、時刻ｔにおけるすべての格子点ｋ（ｋ＝１〜Ｌ）における状態変数Ｕ_ｋの事後分布が揃うこととなる。

次に、出力部６０について説明する。シミュレーションを継続する場合、出力部６０は、判定部５１によって選択された事後分布および第１の事後分布からなる時刻ｔの状態ベクトルを、システムモデル２１に入力する。そして、システムモデル２１は、この時刻ｔの事後分布を用いて、次の時間ステップである時刻ｔ＋１での事前分布を算出する。また、出力部６０は、シミュレーションの結果として、判定部５１によって選択された事後分布および第１の事後分布からなる状態ベクトルの時系列を、出力装置１００６等に出力する。

前述のように、統合後事後分布記憶部５２では、第１の事後分布と、判定部５１によって選択された統合後の事後分布または第２の事後分布とによって、時刻ｔにおけるすべての格子点ｋ（ｋ＝１〜Ｌ）における状態変数Ｕ_ｋの事後分布が揃っている。出力部６０は、統合後事後分布記憶部５２に記憶された各状態変数・格子点での事後分布をシステムモデル２１に入力するとともに、その時系列を出力すればよい。

なお、上述してきたように、各機能ブロックの構成を、状態変数Ｕ_ｋ（ｋ＝１〜Ｌ）を例として説明したが、各機能ブロックは、他の状態変数（例えば、Ｖ_ｋ（ｋ＝１〜Ｌ））についても同様に構成される。

以上のように構成されたシミュレーション装置１００の動作について、図面を参照して説明する。

まず、シミュレーション装置１００が、シミュレーション開始時に実行する動作について、図４を用いて説明する。なお、図４において、シミュレーション開始をステップの基準、すなわち時刻ｔ＝１とする。

図４において、まず、システムモデル２１は、時間・空間上で離散化されたシミュレーションを行うために、時間ステップ、格子点、および、時間発展を求める状態変数を決める（ステップＳ１０１）。例えば、時間ステップや格子点は、所要の精度に基づくか、または、計算が収束するように、適切な値を選べばよい。

次に、入力部１０は、第１〜第Ｍの観測データを取得する（ステップＳ１０２）
次に、データ選択処理部３０は、システムモデル２１で設定された状態変数に関する情報を参照し、第１〜第Ｍの観測データのうち、利用するｍ種の観測データを選択する（ステップＳ１０３）。

次に、データ選択処理部３０は、状態変数とｍ種の観測データとをそれぞれ関係づける関係式、および含まれるノイズ量を設定し、観測モデル３１−１〜３１−ｍを生成する（ステップＳ１０４）。例えば、データ選択処理部３０は、観測データの種類、性質、物理量、観測データおよび状態変数の次元数等に基づいて、関係式およびノイズ量を設定すればよい。これにより、ｍ個の観測モデル３１が生成される。

以上で、シミュレーション装置１００は、シミュレーション開始時の動作を終了する。

次に、シミュレーション装置１００が、シミュレーションを実行する動作について、図５を用いて説明する。

図５において、まず、入力部１０は、状態ベクトルの初期状態を示すアンサンブルおよびパラメータを取得し、システムモデル２１に対して出力する（ステップＳ２０１）。

次に、システムモデル２１は、次の時間ステップのアンサンブル、すなわち事前分布を計算し、事前分布記憶部２２に保存する（ステップＳ２０２）。

ここで、入力部１０は、この時間ステップの時刻において、第１〜第ｍの観測データの少なくともいずれかが得られているか否かを判断する（ステップＳ２０３）。

ここで、いずれの観測データも得られていない場合（ステップＳ２０３でＮｏ）、システムモデル２１は、事前分布記憶部２２に保存した次の時間ステップの事前分布を用いて、再度ステップＳ２０２を実行し、時間ステップをもう１つ進める計算を行う。

なお、いずれかの観測データが得られていても、この時間ステップではデータを修正しないことが指定されている場合も、ステップＳ２０３でＮｏと判断されるものとする。

一方、第１〜第ｍの観測データの少なくともいずれかが得られておりデータを修正する場合（ステップＳ２０３でＹｅｓ）、観測モデル３１−１〜３１−ｍは、事前分布記憶部２２に記憶された事前分布をそれぞれ変換する（ステップＳ２０４）。

なお、このとき、同一のシミュレーション中では、観測モデル３１−１〜３１−ｍとして、シミュレーション開始時にステップＳ１０４で生成されたものが基本的に用いられる。ただし、同一のシミュレーション中であっても、観測データの振る舞いが大きく変化する場合や、シミュレーションの計算がうまくいかないような例外的な場合には、新たに前述のステップＳ１０４が実行されてもよい。この場合、このステップＳ２０４では、新たに生成されたｍ個の観測モデル３１により変換を行ってもよい。

次に、事後分布生成部４０は、生成されたｍ種の変換された事前分布と、この時間ステップの時刻におけるｍ種の観測データとに基づいて、状態変数・格子点ごとに事後分布を生成する（ステップＳ２０５）。これにより、元の事前分布が修正される。

次に、事後分布生成部４０は、ステップＳ２０５で生成した状態変数・格子点ごとの事後分布が、ステップＳ１０３で選択されたｍ種すべての観測データによる事後分布であるか、一部が欠損した観測データによる事後分布であるかを判定する（ステップＳ２０６）。

ここで、ｍ種全ての観測データに基づく事後分布であると判定した場合、事後分布生成部４０は、この事後分布を、第１の事後分布として、第１の事後分布記憶部４１ａに保存する（ステップＳ２０７）。

また、この場合、事後分布生成部４０は、この第１の事後分布を、その状態変数・格子点での事後分布として、統合後事後分布記憶部５２にも保存する（ステップＳ２０８）。

一方、ステップＳ２０６において、この事後分布が、ｍ種の観測データのうち一部に欠損がある観測データによる事後分布であると判定した場合、事後分布生成部４０は、この事後分布を、第２の事後分布とし、その分散値Ｖ０を算出する（ステップＳ２０９）。そして、事後分布生成部４０は、この第２の事後分布およびその分散値Ｖ０を、第２の事後分布記憶部４１ｂに記憶する。

次に、事後分布統合部５０は、第２の事後分布が生成された各状態変数・格子点について、第１の事後分布および第２の事後分布を統合することにより、分散値Ｖが最小となる新たな事後分布（統合後の事後分布）を算出する（ステップＳ３００）。

具体的には、事後分布統合部５０は、対象の状態変数・格子点について、パラメータセットπを変化させながら、（１９）式を用いて統合後の事後分布を繰り返し算出し、分散値Ｖが最小となるπを探索すればよい。例えば、事後分布統合部５０は、最小二乗法やベイズの定理を用いて探索を行っても良い。ここで、探索された分散の最小値をＶｍｉｎとする。

次に、判定部５１は、分散の最小値Ｖｍｉｎと、統合前の分散値Ｖ０とを比較する（ステップＳ３０１）。

ここで、統合後の分散の最小値Ｖｍｉｎの方が小さければ、判定部５１は、統合後の事後分布を、その状態変数・格子点での新たな事後分布とし（ステップＳ３０２）、統合後事後分布記憶部５２に保存する（ステップＳ２０８）。

一方、統合後の分散の最小値Ｖｍｉｎの方が小さくならない場合、判定部５１は、統合を中断し、第２の事後分布を、その状態変数・格子点での事後分布とし（ステップＳ３０３）、統合後事後分布記憶部５２に保存する（ステップＳ２０８）。

次に、シミュレーションが規定の時間または規定のステップに達していなければ（ステップＳ３０４でＮｏ）、シミュレーション装置１００は、ステップＳ２０２からの動作を繰り返す。すなわち、システムモデル２１は、統合後事後分布記憶部５２に保存された各状態変数・格子点の事後分布を入力としてステップＳ２０２を実行し、次ステップの計算を開始する。

一方、シミュレーションが規定の時間または規定のステップに達していれば（ステップＳ３０４でＹｅｓ）、出力部６０は、統合後事後分布記憶部５２に保存された各状態変数・格子点の事後分布の時系列を出力し、シミュレーション動作を終了する。

次に、本発明の第１の実施の形態の効果について述べる。

本発明の第１の実施の形態としてのシミュレーション装置は、非理想的な観測データや、非連続性や特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行うことができる。

その理由について説明する。本実施の形態では、システムモデルが、状態ベクトルの時間発展をシミュレーションする。また、データ選択処理部が、Ｍ種の観測データのうちｍ種を選択する。そして、ｍ種の観測データにそれぞれ対応するｍ個の観測モデルが、システムモデルが算出した次ステップの状態ベクトルの事前分布を、観測データと状態ベクトルとの関係性に基づいて変換する。そして、事後分布生成部が、変換されたｍ種の事前分布と、選択されたｍ種の観測データとに基づいて、状態変数・格子点毎の事後分布を生成する。そして、事後分布生成部が、生成された事後分布のうち、ｍ種すべての観測データに基づくものを第１の事後分布として、ｍ種のうち一部に欠損のある観測データに基づくものを第２の事後分布とする。そして、事後分布統合部が、第２の事後分布が生成された状態変数・格子点について、第１の事後分布および第２の事後分布を統合して新たな事後分布を生成する。そして、判定部が、第２の事後分布および新たな事後分布のいずれを選択するかを決定し、決定した事後分布を、その状態変数・格子点の統合後の事後分布とする。そして、システムモデルは、第１の事後分布および統合後の事後分布からなる状態ベクトルの事後分布を入力として、次ステップの状態ベクトルを算出するからである。

これにより、本実施の形態は、一部の観測データが不適切または誤差を多く含むような場合であっても、事後分布の統合によって、他の観測データも考慮して修正することができる。もしくは、本実施の形態は、そのような観測データを修正に使わなくなるため、誤差の増大を防ぐことができる。また、本実施の形態は、測定頻度が低い観測データについても、測定頻度が高い観測データも考慮して修正することができるため、より高精度なシミュレーションが可能となる。

このような効果について、模式的に示した図６を用いて説明する。図６において、横軸は時間変化を表す。また、シミュレーションモデル上で時間発展を計算する変数を状態変数とし、実際にセンサなどによって測定される値を観測値としている。これは、関連技術のデータ同化によるシミュレーションを模式的に表した図１５と同様である。ただし、図６では、データ処理部によって選択された観測値が２種（観測値１および観測値２）である場合について示している。本実施の形態は、ｍ種の観測データにそれぞれ対応する異なるｍ種の観測モデルを用いることにより、異なる種類の観測データを同種の状態変数と関係付けることができる。図６では、観測値１に対応する観測モデルを観測モデル１として表し、観測値２に対応する観測モデルを観測モデル２と表している。なお、観測値１については、関連技術による図１５のシミュレーションにおける観測値と同種であり同一の値が得られていることを想定した。

この図６および図１５を比較すると、図１５に示した関連技術では、観測値１のシミュレーションにおいて、観測値１の測定頻度の影響を受けるため誤差が積算されてしまっていた。これに対して、図６では、観測値１よりも観測値２がより短い周期で観測されている。このため、観測値による修正の効果が観測値１のシミュレーション値にも表れている。これは、観測値１および観測値２のシミュレーション値とも、観測モデルの差はあっても、同種の状態変数から生成されているためである。したがって、本実施の形態は、観測頻度が高い観測値２だけでなく観測頻度が低い観測値１に関しても、関連技術よりも高精度なシミュレーションを可能としている。また、図１５の関連技術に比べて、不適切または誤差を多く含む観測値（図６における△で示した観測値１）が得られた場合について説明する。このような場合であっても、本実施の形態は、事後分布の統合によって、そのような観測値１のシミュレーション値を他方の観測値２も考慮して修正することができ、また、そのような観測値１を修正に使わなくなる。このため、本実施の形態は、誤差の増大を防ぐことができる。

このように、本実施の形態は、単独では、データが少ない、または空間的・時間的に偏った分布をしている観測データであっても、それらの観測データを複数種類用いることで、シミュレーションをより広域にわたって高分解能かつ高精度とすることができる。今後は、観測技術の進歩やＭ２Ｍ（Ｍａｃｈｉｎｅ−ｔｏ−Ｍａｃｈｉｎｅ）のように大量のセンサからの情報収集により、より多種多様で大量の観測データが集まることが予想される。本実施の形態は、そのようなより多種多様で大量データが集まる状況において、観測データの特性で精度が律速されていた関連技術と比べて、複数の観測データからの情報を統合的に用いることで、より有効なシミュレーションを行うことができる。

（第２の実施の形態）
次に、本発明の第２の実施の形態について図面を参照して説明する。本実施の形態は、空間上は離散的だが値の精度が高い観測データと、空間上は連続的だが精度が十分でない観測データを利用したシミュレーションに適用できる。以下では、本発明のシミュレーション装置を用いて、土壌水分量のシミュレーションを行う具体例について説明する。なお、本発明の第２の実施の形態において参照する各図面において、本発明の第１の実施の形態と同様の構成およびステップには同一の符号を付して、本実施の形態における詳細な説明を省略する。

まず、本発明の第２の実施の形態としてのシミュレーション装置２００の構成を、図７に示す。図７において、シミュレーション装置２００は、本発明の第１の実施の形態としてのシミュレーション装置１００におけるシステムモデル２１として、土壌モデル２２１を適用した構成である。また、シミュレーション装置２００は、本発明の第１の実施の形態としてのシミュレーション装置１００におけるｍ個の観測モデル３１として、２種類の観測データに対応する２個の観測モデル２３１−１〜２３１−２を適用した構成である。また、シミュレーション装置２００は、本発明の第１の実施の形態としてのシミュレーション装置１００に対して、事後分布統合部５０に替えて事後分布統合部２５０を備える点が異なる。

また、本実施の形態では、本発明における初期状態として土壌初期状態と、本発明におけるパラメータとして地形・気象パラメータとが適用される。また、Ｍ＝２個の観測データとして、土壌水分データＯＢＳ_１、および衛星データＯＢＳ_２の２種類が適用される。

ここで、本実施の形態で用いる２種類の観測データ、土壌水分データＯＢＳ_１、および衛星データＯＢＳ_２について説明する。図８に、土壌水分データＯＢＳ_１、および衛星データＯＢＳ_２について、その対象とする計算格子空間（１〜９の９格子）の時系列変化（ｔ−３〜ｔの４ステップ）を模式的に示す。図８において、塗りつぶし部は、観測データが取得された格子点を示す。なお、本実施の形態では、これら２種類の観測データの空間的な取得範囲や間隔は、計算格子空間と同様とみなす。また、データ取得点が格子内で局所的であっても、各格子内における値は一様であるとみなす。

土壌水分データＯＢＳ_１は、例えば、地中に埋めて誘電率から算出する誘電率土壌水分センサから得られた観測値であってもよい。その他、土壌水分データＯＢＳ_１は、他のセンサから得られたものであってもよい。この土壌水分データＯＢＳ_１の特徴は、センサが設置された点の値しか計測出来ないため、空間上は離散的だが、土壌水分に相当する物理量を直接計測しているため高精度である、という点である。図８では、格子点番号１，３，８の３点のみにセンサが設置されていることを想定している。

また、衛星データＯＢＳ_２は、例えば、Ｔｅｒｒａ衛星搭載のセンサＡＳＴＥＲ（Ｔｅｒｒａ／ＡＳＴＥＲ）から得られるリモートセンシングデータであってもよい。詳細には、Ｔｅｒｒａ／ＡＳＴＥＲによる近赤外線（バンド３，０．７８−０．８６μｍ）および短波長赤外（バンド４，１．６００−１．７００μｍ）波長での太陽光に対する反射光強度を表すデータが、衛星データＯＢＳ_２として適用可能である。その他、衛星データＯＢＳ_２には、他の方式や波長のデータを適用してもよい。この衛星データＯＢＳ_２の特徴は、近赤外および短波長赤外波長領域における太陽光に対する地表での反射光強度を２次元的な画像データとして取得することが出来るため、空間上は連続的である、という点である。ただし、この衛星データＯＢＳ_２は、得られたデータに基づいて、これらの波長での反射光強度や反射率などと、地面表層土壌の含水率との間の有意な相関を利用して推定されるため、間接的な値となり、精度が不十分となる可能性がある。

次に、土壌モデル２２１について説明する。土壌モデル２２１は、本発明の第１の実施の形態におけるシステムモデル２１の一例である。土壌モデル２２１は、対象とする土壌の物理的な性質、例えば傾斜や水はけの度合などと、降水量などの気象状態をパラメータとして、土壌水分量などの時空間変化を算出するものである。土壌モデル２２１には、例えば、ＬＳＭ（ＬＡＮＤ−ＳＵＲＦＡＣＥ ＭＯＤＥＬ）を適用してもよい。また、土壌モデル２２１には、農業向けの意思決定支援システムＤＳＳＡＴ（Ｄｅｃｉｓｉｏｎ Ｓｕｐｐｏｒｔ Ｓｙｓｔｅｍ ｆｏｒ Ａｇｒｏｔｅｃｈｎｏｌｏｇｙ Ｔｒａｎｓｆｅｒ）の土壌モジュールなどを適用してもよい。

事後分布統合部２５０は、本発明の第１の実施の形態と同様に、第２の事後分布が生成された状態変数・格子点について、第１の事後分布および第２の事後分布を統合して新たな事後分布を算出する。加えて、事後分布統合部２５０は、統合処理の際に、算出済みの各事後分布の空間的な相関性に基づいて生成されたモデルを用いてもよい。モデルとしては、例えば、共分散関数や、バリオグラム関数が適用可能である。ただし、本発明における事後分布統合部は、各事後分布の空間的な相関性に基づく他のモデルを用いてもよい。また、この場合、事後分布統合部２５０は、統合に用いるモデルの演算を特徴付けるパラメータを、算出済みの各事後分布の空間的な相関性に基づいてベイズ更新してもよい。空間的な相関性に基づくモデルを用いた処理およびそのパラメータをベイズ更新する処理については、以下の動作の説明において具体例とともに説明する。

以上のように構成されたシミュレーション装置２００の動作の具体例について説明する。

まず、土壌モデル２２１は、土壌初期状態および地形・気象パラメータを、入力部１０を介して取得し、格子点ｋにおける土壌水分量ＳＭ_ｋを状態変数に設定する（図４のステップＳ１０１、図５のＳ２０１）。

ここで、土壌水分量シミュレーションにおける時刻ｔでの状態ベクトルは、図８に記載の格子点ｋ（ｋ＝１〜９）における状態変数が土壌水分量ＳＭ_ｋのみとすれば、

・・・（２０）
と表せる。ここでは、１格子点における状態変数を土壌水分量のみと設定した例を中心に説明する。ただし、状態変数としては、時間変化する動的な変数や、値を推定したい量に加えて、静的な変数を適用可能である。また、状態変数は、シミュレーション対象の現象やシステムモデル、目的などに応じて選べばよい。状態変数は、（２）式に示したように、ある時刻の状態ベクトルが１ステップ前の状態ベクトルおよび土壌モデル２２１によって生成可能であればよい。また、状態変数が増えると計算量が増加するため、状態変数は、許された計算資源に応じて適切に設定すべきである。

次に、データ選択処理部３０は、２種の観測データを取得し（図４のステップＳ１０２）、利用するｍ種の観測データとして、土壌水分データＯＢＳ_１および衛星データＯＢＳ_２を選択する（ステップＳ１０３）。

次に、データ選択処理部３０は、土壌水分データＯＢＳ_１に対応した第１の観測モデル２３１−１と、衛星データＯＢＳ_２に対応した第２の観測モデル２３１−２の２つを生成する（ステップＳ１０４）。

ここで、土壌水分データＯＢＳ_１が状態変数ＳＭと同じ次元であり、観測値のノイズがガウス（正規）分布に従う場合を想定する。また、図８に示したように、計算格子点と観測格子点とが一致していることを想定する。この場合、（７）式で表される行列は、単位行列となる。したがって、観測データＯＢＳ_１と状態変数とは、（８）式の観測モデル式に基づいて、

・・・（２１）
という線形な関係で表される。ここで観測ノイズｗ_１は、例えば平均０、分散σ１のガウス分布とすることが出来る。このように、データ選択処理部３０は、（２１）式で表される第１の観測モデル２３１−１を生成する。

また、衛星データＯＢＳ_２について、観測される近赤外および短波長赤外波長での反射光強度または反射率と、土壌水分量とが非線形な関数ｈ_２で関係付けられるとする。ただし観測格子点は同様に計算格子と一致しているとする。この場合、観測データＯＢＳ_２と状態変数とは、（８）式の観測モデル式に基づいて、

・・・（２２）
という非線形な関係で表される。ここでも観測ノイズｗ_２は、例えば平均０、分散σ２のガウス分布とすることが出来る。このように、データ選択処理部３０は、（２２）式で表される第２の観測モデル２３１−２を生成する。

次に、土壌モデル２２１は、シミュレーションの開始点で、土壌初期状態に基づくアンサンブル（（２０）式においてｔ＝０）、地形・気象パラメータ、およびシステムノイズのアンサンブルを取得する。そして、土壌モデル２２１は、（４）式で示されたアンサンブルごとの時間発展式によって、ｔ＝１の状態ベクトルの事前分布を計算し、事前分布記憶部２２に保存する（図５のステップＳ２０２）。

次に、時刻ｔ＝１で観測データＯＢＳ_１およびＯＢＳ_２が得られているとする（ステップＳ２０３でＹｅｓ）。そこで、観測モデル２３１−１および２３１−２は、事前分布記憶部２２に保存された、時刻ｔ＝１での状態ベクトルのアンサンブルを、（２１）式および（２２）式を用いて変換する（ステップＳ２０４）。

次に、事後分布生成部４０は、格子点ごとに、（９）式のベイズの定理により事後分布を計算する（ステップＳ２０５）。ただし、図８に示したように、格子点１、３、８については、２つの観測値ＯＢＳ_１、ＯＢＳ_２が得られているのに対し、格子点２、４、５、６、７、９については、１つの観測値ＯＢＳ_２のみが得られている。したがって、事後分布生成部４０は、前者の格子点１、３、８については、データ選択処理部３０によって選択されたすべての観測データを用いて、（１３）式に基づく次式（２３）を用いて第１の事後分布を計算する。

・・・（２３）
ここで、（２３）式では、ｉ＝１、３、８である。また、ＯＢＳ１ｉ、ＯＢＳ２ｉは、格子点ｉで得られた観測データＯＢＳ_１およびＯＢＳ_２をそれぞれ表すものとする。（２３）式により計算された格子点１、３、８の各第１の事後分布は、第１の事後分布記憶部４１ａに保存される（ステップＳ２０６でＹ、Ｓ２０７）。また、これらの格子点１、３、８の各第１の事後分布は、統合後事後分布記憶部５２にも保存される（ステップＳ２０８）。

また、事後分布生成部４０は、後者の格子点２、４、５、６、７、９については、データ選択処理部３０によって選択された観測データのうち１つが欠如しているので、（１５）式に基づく次式（２４）を用いて第２の事後分布を計算する。

・・・（２４）
ここで、（２４）式では、ｊ＝２、４、５、６、７、９である。（２４）式により計算された格子点２、４、５、６、７、９の各第２の事後分布は、第２の事後分布記憶部４１ｂに保存される（ステップＳ２０６でＮ、Ｓ２０９）。

次に、事後分布統合部２５０は、（２３）式および（２４）式により計算された第１および第２の事後分布を統合する（ステップＳ３００）。

ここでは、（１９）式に示した事後分布統合の関数ｇとして、例えば、周辺格子点の事後分布の線形結合を考える。例えば、図８に示す格子点２については、観測データが衛星データＯＢＳ_２のみであるので、（２４）式で表される第２の事後分布が第２の事後分布記憶部４１ｂに保存されている。この事後分布を、格子点２以外の格子点の事後分布の線形結合で表すと、

・・・（２５）
と表される。ここで、α１〜α９は（１９）式におけるパラメータセットπに相当する重み係数である。以降、（２５）式における確率分布ｐ’のダッシュ（’）は、事後分布統合部２５０での統合後の確率分布であることを示す。このとき、（２５）式は、未知の格子点２における値を、周辺格子点の値とのある確率的相互関係、すなわち空間的相関性に基づいて決定する、いわゆるクリギング（Ｋｒｉｇｉｎｇ）法と同様とみなすことが出来る。ただし、各格子点における値は確定値ではなく、（２３）式および（２４）式から求まる事後分布である。すなわち、格子点ｋの位置ｒ_ｋと距離γ離れた点ｒ_ｋ＋γとの間の、土壌水分量の事後分布ｐ（ＳＭ｜ＯＢＳ）の空間的な相関性を表す共分散関数

・・・（２６）
が求まれば、（２５）式の重み係数α１〜α９、すなわちパラメータセットπも求まる。なお、（２６）式において、ＳＭ（ｘ）は、位置ｘにおける格子点の状態変数ＳＭを表す。また、ＯＢＳは、ｍ種の観測データを表す。これは、例えば次式（２７）に示すような、単純型クリギング方程式を解くことで求めることができる。なお、本発明において、事後分布統合部が事後分布を統合する際に用いる関数ｇのパラメータπを求める手法は、この手法に限定されず、他の手法であってもよい。

・・・（２７）
次に（２６）式の共分散関数を求める動作を説明する。共分散関数Ｃ（γ）とバリオグラム関数Ｖ（γ）の間には簡単な関係

・・・（２８）
が成り立つため、いずれか一方を求めれば良い。以下では、事後分布統合部２５０において、バリオグラム関数Ｖ（γ）を先に求める場合について説明する。バリオグラムも、共分散関数と同様に、格子ｋの位置ｒ_ｋと距離γ離れた点ｒ_ｋ＋γとの間の確率的相互作用、すなわち空間的な相関性を表す。図９の左側に、バリオグラム推定結果の一例を示す。この例は、統合によって算出する格子点以外の格子点のバリオグラムを算出した結果に対し、指数型のバリオグラムモデル

・・・（２９）
にあてはめて、そのパラメータであるζを推定した。ここで、ζは、バリオグラムを特徴付ける３種のパラメータセットを表しており、一般にナゲットτ^２、レンジφ、シルσ^２と呼ばれる。ここでは、これらのパラメータのうち、レンジφとナゲットτ^２について、ベイズの定理（９）式を用いて推定を行った結果を示している。具体的には、レンジφに一様な事前分布、ナゲットτ^２は０に近い値が予想されたため、指数型の事前分布をそれぞれ仮定し、実際に算出されたバリオグラムの結果から、ベイズの定理によって事後分布をそれぞれ求めた。その結果の例を図９の右側に示した。この図から明らかなように、事後分布には最大値が見られ、この最大値は、算出されたバリオグラムを最も良く再現するパラメータの値、最尤値とみなすことができる。図９左側の曲線（推定値）は、このパラメータの下で（２９）式に基づいて描いたものである。以上より、バリオグラムＶ（γ）の関数が算出できたので、（２８）式より共分散関係も算出できる。なお、ここではベイズの定理を使ったパラメータ推定方法を示したが、あくまでも例示であって他の方法を用いても良い。また、図９は推定結果の一例であって、図８に記載の格子空間（１〜９）に対応するものではない。

したがって、（２６）式で表す共分散関数Ｃ（γ）が求まったので、例えば（２７）式に示す単純型クリギング方程式を解くことができる。これより、事後分布の統合を記述する（２５）式の係数、すなわちパラメータセットπが算出されるので、事後分布統合部２５０は、格子点２の統合後の事後分布ｐ’（ＳＭ_２｜ＯＢＳ_１，ＯＢＳ_２）を求めることができる。また、事後分布統合部２５０は、第２の事後分布が生成された他の格子点ｋについても、統合後の事後分布ｐ’（ＳＭ_ｋ｜ＯＢＳ_１，ＯＢＳ_２）を同様にして求める。

ここで、ステップＳ３００における事後分布統合部２５０の統合動作の詳細を、図１０に示す。なお、図１０には、第２の事後分布が算出されたある格子点を対象とした統合動作を示している。

図１０では、まず、事後分布統合部２５０は、対象の格子点以外の格子点について、バリオグラムまたは共分散を計算する（ステップＳ４０１）。

次に、事後分布統合部２５０は、ステップＳ４０１で計算したバリオグラムまたは共分散にあてはまる関数を定義する（ステップＳ４０２）。

次に、事後分布統合部２５０は、ステップＳ４０２で定義した関数のパラメータに事前分布を仮定する（ステップＳ４０３）。

次に、事後分布統合部２５０は、ステップＳ４０３で仮定したパラメータの事前分布を、バリオグラムまたは共分散の計算結果に基づきベイズの定理を用いて更新することにより、パラメータの事後分布を求める（ステップＳ４０４）。

次に、事後分布統合部２５０は、ステップＳ４０４で得られたパラメータ事後分布を用いて、共分散関数を導出する（ステップＳ４０５）。

次に、事後分布統合部２５０は、クリギング方程式により、対象の格子点以外の格子点の事後分布を統合する際の重み付け係数（パラメータセットπ）を求める（ステップＳ４０６）。

次に、事後分布統合部２５０は、ステップＳ４０６で用いたパラメータセットπを用いて、対象の格子点以外の格子点の事後分布を統合する（ステップＳ４０７）。

このようにして、図５のステップＳ３００において、統合後の事後分布が算出された。

以降、シミュレーション装置２００は、ステップＳ３０１〜Ｓ３０４、Ｓ２０８を本発明の第１の実施の形態と同様に実行する。これにより、第２の事後分布が生成された各格子点について、統合後の事後分布または第２の事後分布が、統合後事後分布記憶部５２に保存される。そして、土壌モデル２２１は、統合後事後分布記憶部５２に保存されたこの時間ステップの事後分布からなる状態ベクトルを用いて、次の時間ステップの計算を継続する。

次に、本発明の第２の実施の形態の効果について述べる。

本発明の第２の実施の形態としてのシミュレーション装置は、非理想的な観測データや、非連続性や特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行うことができる。

その理由について説明する。本実施の形態では、本発明の第１の実施の形態と同様の構成に加えて、次の構成を備えるからである。すなわち、事後分布統合部が、第２の事後分布が生成された格子点について、第１の事後分布および第２の事後分布を統合する際に、算出済みの各事後分布の空間的な相関性に基づいて生成されたモデルを用いるからである。また、事後分布統合部が、統合に用いるモデルの演算を特徴付けるパラメータを、算出済みの各事後分布の空間的な相関性に基づいてベイズ更新するからである。

これにより、本実施の形態は、単独では、データが少ない観測データや、空間的に偏った分布をしている観測データであっても、そのような観測データを複数種類用いることで、格子点毎の事後分布をより精度よく統合することができる。その結果、本実施の形態は、シミュレーションをより広域にわたって高分解能かつ高精度とすることができる。

なお本発明の第２の実施の形態において、システムモデルに土壌モデルを適用し、複数の観測データに土壌センサデータと衛星データとを適用し、土壌水分量のシミュレーションを行う例について説明した。この他、本実施の形態は、他の対象について他のシステムモデルや観測データを用いて実施することも可能である。例えば、本実施の形態は、システムモデルに気象モデルを適用し、複数の観測データに気象センサデータおよび衛星データを適用してもよい。

（第３の実施の形態）
次に、本発明の第３の実施の形態について図面を参照して説明する。本実施の形態は、複数の観測データ間で観測格子間隔が異なる場合や、取得時間間隔が異なる場合のシミュレーションに適用できる。以下では、本発明のシミュレーション装置を用いて、作物生育（成長）についてシミュレーションを行う具体例について説明する。なお、本発明の第３の実施の形態において参照する各図面において、本発明の第１の実施の形態と同様の構成およびステップには同一の符号を付して、本実施の形態における詳細な説明を省略する。

まず、本発明の第３の実施の形態としてのシミュレーション装置３００の構成を、図１１に示す。図１１において、シミュレーション装置３００は、本発明の第１の実施の形態としてのシミュレーション装置１００におけるシステムモデル２１として、作物モデル３２１を適用した構成である。また、シミュレーション装置３００は、本発明の第１の実施の形態としてのシミュレーション装置１００におけるｍ個の観測モデル３１として、２種類の観測データに対応する２個の観測モデル３３１−１〜３３１−２を適用した構成である。また、シミュレーション装置３００は、本発明の第１の実施の形態としてのシミュレーション装置１００に対して、事後分布統合部５０に替えて事後分布統合部３５０を備える点が異なる。

また、本実施の形態では、本発明における初期状態として土壌初期状態と、本発明におけるパラメータとして地形・気象・作物パラメータとが適用される。また、Ｍ＝２個の観測データとして、２種類の衛星データ（リモートセンシングデータ）が適用される。

ここで、本実施の形態で用いる２種類の観測データ、衛星データＯＢＳ_１、および衛星データＯＢＳ_２について説明する。

図１２に示すように、本実施の形態では、第１の衛星データは、高頻度かつ低空間分解能であり、第２の衛星データは、低頻度かつ高空間分解能である。なお、図１２には、第１の衛星データＯＢＳ_１および第２の衛星データＯＢＳ_２について、対象とする計算格子空間（１〜１６の１６格子）での時系列変化（ｔ−３〜ｔの４ステップ）を模式的に示している。また、図１２において、塗りつぶし部は、データが取得された格子点を示している。

高頻度かつ低空間分解能な第１の衛星データは、例えば、Ｔｅｒｒａ衛星またはＡＱＵＡ衛星搭載のセンサＭＯＤＩＳ（Ｔｅｒｒａ・ＡＱＵＡ／ＭＯＤＩＳ）から得られるデータであってもよい。詳細には、Ｔｅｒｒａ・ＡＱＵＡ／ＭＯＤＩＳによる可視赤バンド（波長０．５８−０．８６μｍ）および近赤外バンド（波長０．７２５−１．１００μｍ）での太陽光に対する反射光強度のデータが、第１の衛星データとして適用可能である。このような第１の衛星データは、取得する地域の緯度に依存するが、基本的に毎日取得可能である。しかしながら、このような第１の衛星データは、地上での空間分解能は、約２５０ｍと低い。

また、低頻度かつ高空間分解能な第２の衛星データとしては、例えば、ＬＡＮＤＳＡＴ衛星、ＰＬＥＩＡＤＥＳ衛星、または、ＡＳＮＡＲＯ衛星などから得られる観測データがある。これらにより取得される波長域は、第１の衛星データとして取得される波長とほぼ同一である。このような第２の衛星データの取得頻度と地上分解能は、ＬＡＮＤＳＡＴ衛星の場合、８〜１６日間隔で約３０ｍ、ＰＬＥＩＡＤＥＳ衛星、ＡＳＮＡＲＯ衛星の場合、２〜３日間隔で約２ｍである。

なお、作物の生育状態を示す植生指標として一般的に用いられているＮＤＶＩ（Ｎｏｒｍａｌｉｚｅｄ Ｄｉｆｆｅｒｅｎｃｅ Ｖｅｇｅｔａｔｉｏｎ Ｉｎｄｅｘ：正規化植生指標）は、前述の２つのバンド（可視赤および近赤外）の反射率から算出できる。ただし、観測データとして取得される波長域は、必ずしもこれらのバンドに限定されるものではない。本実施の形態では、作物モデル３２１は、作物の生育状態を表す量として、葉面積指数（Ｌｅａｆ Ａｒｅａ Ｉｎｄｅｘ：ＬＡＩ）を算出する。ＬＡＩは、植生指数ＮＤＶＩと相関があることが知られている。このようなＬＡＩは、土壌初期状態や地形・気象・作物パラメータが作物モデル３２１へデータが入力されることで算出される。

また、本実施の形態が、前述した本発明の他の実施の形態と異なる点の１つは、２種の観測データＯＢＳ_１およびＯＢＳ_２の格子間隔が異なる点である。そこで、ＯＢＳ_１およびＯＢＳ_２の少なくともいずれか一方の観測データと一致するように、作物モデル３２１の計算格子が設定される。なお、他方の観測データについては、（７）式に示した観測モデル式におけるベクトルを、例えば近傍の格子点の値の加重平均や重み付け平均となるように変更すれば良い。また、本実施の形態が、前述した本発明の他の実施の形態と異なる点のもう１つは、２種の観測データＯＢＳ_１およびＯＢＳ_２の取得時間の間隔が異なる点である。そこで、事後分布統合部３５０は、取得頻度の少ない観測データＯＢＳ_２から得られた事後分布を、時間的な相関関係に基づいて推定することにより、取得頻度の多い観測データＯＢＳ_１から得られた事後分布の取得時間に合わせて統合する。

以下、図１２に示した２種の観測データＯＢＳ_１およびＯＢＳ_２の例を用いて、シミュレーション装置３００の動作について具体的に説明する。

まず、作物モデル３２１は、図１２に記載の格子点ｋ（ｋ＝１〜１６）における状態変数として、葉面積指数ＬＡＩ_ｋを設定する（図４のステップＳ１０１）。なお、１格子点における状態変数の設定は、変数の時間に対する依存性や、推定したい未知量の数に応じて選べばよい。

次に、データ選択処理部３０は、２種の観測データを取得し（ステップＳ１０２）、利用するｍ種の観測データとして、第１の衛星データＯＢＳ_１および第２の衛星データＯＢＳ_２を選択する（ステップＳ１０３）。

次に、データ選択処理部３０は、第１の衛星データＯＢＳ_１に対応した第１の観測モデル３３１−１と、第２の衛星データＯＢＳ_２に対応した第２の観測モデル３３１−２の２つを生成する（ステップＳ１０４）。

ここで、図１２を参照すると、第１の観測データＯＢＳ_１については、１つの観測データ取得格子点（図中の塗りつぶし部）の値を、重なりがある４つの計算格子点の平均値と対応づけることができる。したがって、観測モデル３３１−１は、４点の観測格子点における第１の観測データＯＢＳ_１と、状態変数ＬＡＩ_ｋ（ｋ＝１〜１６）との関係により、

・・・（３０）
と表される。

また、観測モデル３３１−２は、第２の観測データＯＢＳ_２が取得される格子点と計算格子点とが１対１に対応することから、単位行列を用いて、

・・・（３１）
と表される。ここで、Ｈ_１およびＨ_２は、それぞれの観測データと状態変数ＬＡＩとを対応づける写像ｈおよび格子点を対応づける行列を含む写像である。また、ｗ_１およびｗ_２は、観測ノイズであり、例えば平均０、分散σのガウス分布などに設定可能である。なお、（３０）式および（３１）式の観測モデル３３１−１および３３１−２は、（８）式で表される観測モデル式の具体例である。

次に、作物モデル３２１は、土壌初期状態および地形・気象・作物パラメータを取得し、シミュレーションにおける次ステップでの状態ベクトルの事前分布を計算する（図５のステップＳ２０１、Ｓ２０２）。そして、観測モデル３３１−１および３３１−２は、（３０）式および（３１）式によって事前分布を変換する（ステップＳ２０３）。ここでは、図５に示す動作が適宜繰り返された後のステップＳ２０３において、図１２の時刻ｔ−１における変換された事前分布ｐ（ＬＡＩ_ｋ）が計算されたものとする。

次に、事後分布生成部４０は、図１２の時刻ｔ−１における格子点ｋ（ｋ＝１〜１６）の状態変数ＬＡＩｋの事後分布を生成する。時刻ｔ−１では、観測データがＯＢＳ_１およびＯＢＳ_２の全て得られていることから、事後分布生成部４０は、第１の事後分布として、

・・・（３２）
を計算し、第１の事後分布記憶部４１ａに保存する（ステップＳ２０５、Ｓ２０６でＹｅｓ、Ｓ２０７）。ここで、ＬＨは（１３）式で表される尤度を計算する関数で、分母のＺは規格化定数である。

また、ステップＳ２０３において、図１２の時刻ｔにおける変換された事前分布ｐ（ＬＡＩ_ｋ）が計算された場合についても説明する。この場合、時刻ｔでは、第２の観測データＯＢＳ_２が欠如しているため、事後分布生成部４０は、第２の事後分布として、

・・・（３３）
を計算し、第２の事後分布記憶部４１ｂに保存する（ステップＳ２０５、Ｓ２０６でＮｏ、Ｓ２０９）。

次に、事後分布統合部３５０は、（３２）式で表された第１、および（３３）式で表された第２の事後分布を統合する。具体的には、（１９）式に示した事後分布統合の関数ｇとして、ここでは、現時刻の第２の事後分布と、異なる時刻の第１の事後分布から時間的な相関関係に基づいて推定された事後分布との線形結合が適用される。例えば、図１２の時刻ｔでの事後分布は、上述のように、観測データが第１の衛星データＯＢＳ_１のみであるので、（３３）式で表される第２の事後分布が、第２の事後分布記憶部４１ｂに保存されている。そこで、事後分布統合部３５０は、この時刻ｔの第２の事後分布と統合するための時刻ｔの第１の事後分布を、時刻ｔより前の過去の時刻において生成した事後分布から時間的な相関関係により推定し生成する。時刻ｔ＝１からｔ−１までの事後分布から推定された時刻ｔの第１の事後分布を、ｐ（ＬＡＩ_ｋ｜ＯＢＳ_１，ＯＢＳ_２）_{１：ｔ−１}と表す。そして、事後分布統合部３５０は、時刻ｔの第２の事後分布ｐ（ＬＡＩ_ｋ｜ＯＢＳ_１）_ｔと、時間的な相関関係に基づいて推定された時刻ｔの第１の事後分布ｐ（ＬＡＩ_ｋ｜ＯＢＳ_１，ＯＢＳ_２）_{１：ｔ−１}とを次式（３４）により統合する。

・・・（３４）
ここで、α_０、β_０は、（１９）式におけるパラメータセットπに相当する重み係数である。また、（３４）式において、確率分布ｐ’のダッシュ（’）は、事後分布統合部３５０での統合後の確率分布であることを示す。

ここで、時刻ｔより前（ｔ−１，ｔ−２，ｔ−３，．．．）の事後分布から時間的な相関関係に基づいて、時刻ｔの事後分布ｐ（ＬＡＩ_ｋ｜ＯＢＳ_１，ＯＢＳ_２）_{１：ｔ−１}を推定する処理の具体例について説明する。一般に、時刻ｔより前（ｔ−１，ｔ−２，ｔ−３，．．．）の値から時刻ｔの値を推定する方法として、いわゆる自己回帰（ＡＲ：autoregressive）モデル

・・・（３５）
を適用することができる。ここで、例として、ＡＲモデルｆ_ＡＲが線形に記述される場合を考える。また、第１の衛星データＯＢＳ_１および第２の衛星データＯＢＳ_２のいずれもが観測されて第１の事後分布が生成されている時刻で、かつ、時刻ｔより前の時刻をｔ−ｉとする（図１２ではｉ＝１，３）。このような時刻ｔ−ｉにおける第１の事後分布からのみ推定する場合の（３５）式は、具体的に、

・・・（３６）
と表される。

また、第２の観測データＯＢＳ_２が得られず、第２の事後分布が生成されている時刻（図１２では、ｔ−２，ｔ−４，．．．）の事後分布も考慮する場合を考える。この場合、時刻ｔより前の時刻について統合後事後分布記憶部５２に保存されている統合後の事後分布を用いて、（３５）式は、具体的に、

・・・（３７）
と表される。ただし、（３７）式において、時刻ｔ−ｉは、第１の事後分布が算出されている時刻、時刻ｔ−ｊは、第２の事後分布が算出されている時刻をそれぞれ表す。図１２の場合、ｉ＝１，３、ｊ＝２であり、ｉ≠ｊである。本実施の形態では、（３７）式において、統合後事後分布記憶部５２に保存されている統合後の事後分布を用いる場合を想定している。このため、図１１において、統合後事後分布記憶部５２から事後分布統合部３５０に対して、統合後の事後分布に関する情報を伝達するデータ経路を矢印で示している。

このようにして推定した時刻ｔの事後分布ｐ（ＬＡＩ_ｋ｜ＯＢＳ_１，ＯＢＳ_２）_{１：ｔ−１}を用いて、事後分布統合部３５０は、（３４）式により統合を行う（ステップＳ３００）。

以降、シミュレーション装置３００は、ステップＳ３０１〜Ｓ３０４、Ｓ２０８を本発明の第１の実施の形態と同様に実行する。これにより、各格子点において第２の事後分布が生成された時刻ｔにおいて、各格子点について統合後の事後分布または第２の事後分布が、統合後事後分布記憶部５２に保存される。そして、作物モデル３２１は、統合後事後分布記憶部５２に保存された時刻ｔの事後分布からなる状態ベクトルを用いて、次の時間ステップの計算を継続する。

そして、規定時間になると、シミュレーション装置３００は、動作を終了する。

次に、本発明の第３の実施の形態の効果について述べる。

本発明の第３の実施の形態としてのシミュレーション装置は、非理想的な観測データや、非連続性や特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行うことができる。

その理由について説明する。本実施の形態では、本発明の第１の実施の形態と同様の構成に加えて、次の構成を備えるからである。すなわち、事後分布統合部が、第２の事後分布が生成された格子点について、事後分布を統合する際に、過去に算出済みの事後分布の時間的な相関性に基づいて生成されたモデルを用いるからである。

このように、本実施の形態は、第２の事後分布が生成された時刻ｔより過去に算出済みの事後分布から時間的な相関関係に基づいて時刻ｔの事後分布を推定し、推定した事後分布を用いて時刻ｔでの統合後の事後分布を算出する。これによって、本実施の形態は、取得頻度が異なる複数の観測データを、観測頻度がより高い観測データの取得タイミングに合わせて統合することが可能となる。すなわち、本実施の形態では、時刻ｔにおける事前分布（シミュレーション結果）が、より確からしい統合後の事後分布により短い間隔で修正されることになる。その結果、本実施の形態は、次の時間ステップ以降の値を予測する際の誤差を小さくすることが出来る。

なお、本実施の形態において、システムモデルに作物モデルを適用し、複数の観測データとしていずれも衛星データを適用した例について説明したが、システムモデルおよび観測データの種類や内容を限定するものではない。例えば、本実施の形態では、システムモデルとして水動態・流体モデルを適用し、観測データとして河川の水位センサデータおよび衛星データを適用することも可能である。このように、本実施の形態は、高頻度だが局所離散的な観測データと、低頻度だが高分解能で面的な観測データとの組み合わせについて、適宜対応するシステムモデルを用いて適用することができる。

（第４の実施の形態）
次に、本発明の第４の実施の形態について図面を参照して説明する。本実施の形態では、本発明のシミュレーション装置を用いて、降水量についてシミュレーションを行う具体例について説明する。なお、本発明の第４の実施の形態は、本発明の第２の実施の形態における計算格子空間を３次元に拡張したものである。また、本発明の第４の実施の形態において参照する各図面において、本発明の第２の実施の形態と同様の構成およびステップには同一の符号を付して、本実施の形態における詳細な説明を省略する。

まず、本発明の第４の実施の形態としてのシミュレーション装置４００の構成を図１３に示す。図１３において、シミュレーション装置４００は、本発明の第２の実施の形態としてのシミュレーション装置２００において、土壌モデル２２１に替えて気象モデル４２１を適用した構成である。また、シミュレーション装置４００は、本発明の第２の実施の形態としてのシミュレーション装置２００において、２つの観測モデル２３１−１および２３１−２に替えて、２つの観測モデル４３１−１および４３１−２を適用した構成である。また、本実施の形態では、本発明における初期状態として気象値初期状態と、本発明におけるパラメータとして地形パラメータとが適用される。また、Ｍ＝２個の観測データとして、ＧＰＳ可降水量データＯＢＳ_１と、音波レーダデータＯＢＳ_２とが適用される。

ここで、本実施の形態で想定する２種類の観測データ、ＧＰＳ可降水量データＯＢＳ_１、および音波レーダデータＯＢＳ_２について説明する。ＧＰＳ可降水量とは、ＧＰＳ（Global Positioning System）衛星から放射された電波がＧＰＳ受信機に届くまでに、大気中の水蒸気が多いほど到達時間が遅れる性質を利用し、大気中の水蒸気量の鉛直積算量を見積もったデータである。ＧＰＳ可降水量は、局地的豪雨の発生するタイミングの推定や、１度の降雨での総雨量の推定の精度向上に役立っている。ＧＰＳ可降水量は、地上側には、ＧＰＳ受信機を配備すれば良いため、陸上の面内においては比較的高密度化がし易いという特徴がある。一方で、ＧＰＳ可降水量は、鉛直方向については、あくまでも鉛直方向の積算値であるため、空間分布を適切に表現することが困難である。これに対して、音波レーダを用いると、水蒸気量の高度依存性を測定することができる。例えば、鉛直上向きに音波を出し、大気の乱流散乱エコーを受けると、そのエコーは、大気屈折率の高度勾配に依存する。また、さらに、大気屈折率の高度勾配は、水蒸気量の高度勾配に強く依存する。したがって、このエコーを観測することで、水蒸気量の高度依存性が測定可能となる。

このような２種の観測データと状態変数との関係を表す観測モデル４３１について、図１４を用いて説明する。図１４において、計算格子点１〜８は、３次元空間上に配置されている。また、ここでは、状態ベクトルは、図１４に記載の格子点ｋ（ｋ＝１〜８）における状態変数として、降水量ＲＡＩＮ_ｋのみが設定される。

ここでは、ＧＰＳ可降水量データＯＢＳ_１については、１つの観測格子点での値は、ｘｙ平面の座標が同じで、ｚ（鉛直）座標が異なる２つの計算格子点の積算値と対応づけることができる。図１４では、ＯＢＳ_１は、ｘｙ座標が同一の計算格子点１および５の積算値と、２および６の積算値と、４および８の積算値とそれぞれ対応付け可能な各観測格子点で取得されたデータである。したがって、第１の観測データＯＢＳ_１と、状態変数ＲＡＩＮ_ｋ（ｋ＝１〜８）との関係を表す観測モデル４３１−１は、次式（３８）で表される。なお、（３８）式は、（８）式で表される観測モデルの具体例である。

・・・（３８）
次に、音波レーダデータＯＢＳ_２については、１つの観測データ取得格子点の値は、ｚ（鉛直）座標が同じ、すなわち同一平面内の４つの計算格子点の平均値と対応づけることができる。図１４では、ＯＢＳ_２は、ｚ座標が同一の計算格子点１〜４の平均値と、５〜８の平均値とそれぞれ対応付け可能な観測点で取得されたデータである。したがって、第２の観測データＯＢＳ_２と、状態変数ＲＡＩＮ_ｋ（ｋ＝１〜８）との関係を表す観測モデル４３１−２は、次式（３９）で表される。なお、（３９）式は、（８）式で表される観測モデルの具体例である。

・・・（３９）
以上のように構成されたシミュレーション装置４００は、本発明の第２の実施の形態としてのシミュレーション装置２００と略同様に動作する。

つまり、気象モデル２４１は、気象値初期状態および地形パラメータを元に計算された次の時間ステップの状態ベクトルの事前分布を算出する（図５のステップＳ２０１〜Ｓ２０２）。そして、上述の２つの観測モデル４３１−１〜４３１−２は、事前分布をそれぞれ変換する（ステップＳ２０３〜Ｓ２０４）。そして、事後分布生成部４０は、各格子点での第１の事後分布または第２の事後分布を生成する（ステップＳ２０５〜Ｓ２０７、Ｓ２０９）。このとき、格子点３、７では、第１の観測データＯＢＳ_１が観測されていないので、第２の事後分布が生成される。また、それ以外の格子点では、第１の事後分布が生成される。

そして、事後分布統合部２５０および判定部５１は、第２の事後分布が生成された格子点３、７について、統合後の事後分布を生成し、生成した統合後の事後分布または元の第２の事後分布のいずれかを決定する（ステップＳ３００〜Ｓ３０３）。そして、気象モデル４２１は、各格子点についての第１の事後分布または決定された事後分布からなる状態ベクトルを用いてシミュレーションを続ける。また、出力部６０は、規定時間になると（ステップＳ３０４でＹｅｓ）、状態ベクトルの時系列を出力し、動作を終了する。

このように、本発明の第４の実施の形態としてのシミュレーション装置は、観測データが単純に格子点と１対１に対応をつけられない場合や、３次元空間のシミュレーションであっても適用可能である。このような場合であっても、本実施の形態は、適切な観測モデルを用いることにより、非理想的な観測データや、非連続性や特異性がある観測データを考慮しながら、広域にわたって高分解能かつ高精度なシミュレーションを行うことができる。

なお、上述した本発明の各実施の形態において、シミュレーション装置の各機能ブロックが、記憶装置またはＲＯＭに記憶されたコンピュータ・プログラムを実行するＣＰＵによって実現される例を中心に説明した。これに限らず、各機能ブロックの一部、全部、または、それらの組み合わせが専用のハードウェアにより実現されていてもよい。

また、上述した本発明の各実施の形態において、シミュレーション装置の機能ブロックは、複数の装置に分散されて実現されてもよい。

また、上述した本発明の各実施の形態において、各フローチャートを参照して説明したシミュレーション装置の動作を、本発明のコンピュータ・プログラムとしてコンピュータの記憶装置（記憶媒体）に格納しておいてもよい。そして、係るコンピュータ・プログラムを当該ＣＰＵが読み出して実行するようにしてもよい。そして、このような場合において、本発明は、係るコンピュータ・プログラムのコードあるいは記憶媒体によって構成される。

また、上述した各実施の形態は、適宜組み合わせて実施されることが可能である。

以上、上述した実施形態を模範的な例として本発明を説明した。しかしながら、本発明は、上述した実施形態には限定されない。即ち、本発明は、本発明のスコープ内において、当業者が理解し得る様々な態様を適用することができる。

この出願は、２０１４年８月２７日に出願された日本出願特願２０１４−１７２３７１を基礎とする優先権を主張し、その開示の全てをここに取り込む。

１００、２００、３００、４００ シミュレーション装置
１０ 入力部
２１ システムモデル
２２１ 土壌モデル
３２１ 作物モデル
４２１ 気象モデル
２２ 事前分布記憶部
３０ データ選択処理部
３１、２３１、３３１、４３１ 観測モデル
４０ 事後分布生成部
４１ａ 第１の事後分布記憶部
４１ｂ 第２の事後分布記憶部
５０、２５０、３５０ 事後分布統合部
５１ 判定部
５２ 統合後事後分布記憶部
６０ 出力部
１００１ ＣＰＵ
１００２ ＲＡＭ
１００３ ＲＯＭ
１００４ 記憶装置
１００５ 入力装置
１００６ 出力装置

Claims (10)

1. シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとを入力として取得する入力手段と、
   前記初期状態およびパラメータを基に、前記状態ベクトルの時間発展をシミュレーションするシステムモデルと、
   前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択するデータ選択処理手段と、
   選択された複数の観測データにそれぞれ対応する複数の観測モデルであって、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換して出力する複数の観測モデルと、
   前記複数の観測モデルから出力される状態ベクトルおよび前記データ選択処理手段で選択された観測データに基づいて、前記状態ベクトルの事後分布を生成し、前記データ選択処理手段で選択された全ての観測データに基づく事後分布を第１の事後分布として出力し、１つ以上欠如した観測データに基づく事後分布を第２の事後分布として出力する事後分布生成手段と、
   前記第１の事後分布と前記第２の事後分布とを統合する事後分布統合手段と、
   前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定する判定手段と、
   前記判定手段によって決定された事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力するとともに、該状態ベクトルの時系列を出力する出力手段と、
   を備えたシミュレーション装置。
2. 前記データ選択処理手段は、前記システムモデルで設定された状態ベクトルおよび前記各観測データにそれぞれ関する情報を比較することにより、前記各観測データに対応する前記観測モデルを生成することを特徴とする請求項１に記載のシミュレーション装置。
3. 前記データ選択処理手段は、前記各観測データに対応する前記観測モデルのノイズ量を設定することを特徴とする請求項１または請求項２に記載のシミュレーション装置。
4. 前記事後分布統合手段は、前記第１の事後分布および前記第２の事後分布を統合する処理において、算出済みの各事後分布の相関性に基づいて生成されたモデルを用いることを特徴とする請求項１から請求項３のいずれか１項に記載のシミュレーション装置。
5. 前記事後分布統合手段は、前記統合に用いる前記モデルの演算を特徴付けるパラメータを、算出済みの各事後分布の相関性に基づいてベイズ更新することを特徴とする請求項４に記載のシミュレーション装置。
6. 前記判定手段は、前記第２の事後分布および前記統合後の事後分布の各分散値に基づいて、いずれを用いるかを決定することを特徴とする請求項１から請求項５のいずれか１項に記載のシミュレーション装置。
7. 前記状態ベクトルは、シミュレーションの対象領域内で離散化された格子点ごとの状態変数からなり、
   前記観測モデルは、前記状態変数の格子点と、前記複数の観測データの観測点の解像度とを、前記観測データ毎に関係付けることを特徴とする請求項１から請求項６のいずれか１項に記載のシミュレーション装置。
8. 前記各状態変数の確率分布が、離散化されて独立に計算されるアンサンブルの集合によってそれぞれ近似され、
   前記事後分布統合手段は、前記アンサンブルの集合によって近似される前記各状態変数の確率分布を、所定の割合で重ねあわせることにより統合することを特徴とする請求項１から請求項７のいずれか１項に記載のシミュレーション装置。
9. シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとが入力されると、
   前記初期状態およびパラメータを基に、システムモデルを用いて前記状態ベクトルの時間発展をシミュレーションし、
   前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択し、
   選択された複数の観測データにそれぞれ対応する複数の観測モデルを用いて、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換し、
   前記複数の観測モデルから出力される状態ベクトルおよび選択された前記観測データに基づいて、前記状態ベクトルの事後分布を生成し、
   選択された全ての観測データに基づく第１の事後分布と、１つ以上欠如した観測データに基づく第２の事後分布とを統合し、
   前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定し、
   決定した事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力し、
   決定した事後分布および前記第１の事後分布からなる状態ベクトルの時系列を出力するシミュレーション方法。
10. シミュレーションにおける状態ベクトルの初期状態およびパラメータと、複数の観測データとを入力として取得する入力ステップと、
    システムモデルを用いて、前記初期状態およびパラメータを基に、前記状態ベクトルの時間発展をシミュレーションするシステムモデル計算ステップと、
    前記システムモデルにおける前記状態ベクトルに関連する情報に基づいて、前記複数の観測データのうち利用する複数の観測データを選択するデータ選択処理ステップと、
    選択された複数の観測データにそれぞれ対応する複数の観測モデルを用いて、それぞれ、前記システムモデルから出力される状態ベクトルを、前記観測データおよび前記状態ベクトルの関係性に基づき変換して出力する観測モデル計算ステップと、
    前記複数の観測モデルから出力される状態ベクトルおよび前記データ選択処理ステップで選択された観測データに基づいて、前記状態ベクトルの事後分布を生成し、前記データ選択処理ステップで選択された全ての観測データに基づく事後分布を第１の事後分布として出力し、１つ以上欠如した観測データに基づく事後分布を第２の事後分布として出力する事後分布生成ステップと、
    前記第１の事後分布と前記第２の事後分布とを統合する事後分布統合ステップと、
    前記第２の事後分布および前記統合後の事後分布のいずれを用いるかを決定する判定ステップと、
    前記判定ステップで決定された事後分布および前記第１の事後分布からなる状態ベクトルを前記システムモデルに入力するとともに、該状態ベクトルの時系列を出力する出力ステップと、
    をコンピュータ装置に実行させるコンピュータ・プログラムを記憶した記憶媒体。

Images