(57)【特許請求の範囲】
【請求項１】
処理部が、予め測定された耕地属性の時空間分布データを入力部により入力し、時刻毎に耕地属性の時空間分布データを入力データ格納部に記憶する入力ステップと、
処理部が、前記入力データ格納部を参照して耕地属性の時空間分布データを読み取り、対象とする空間フレームを、複数の分割領域に分割し、さらに分割領域を各セルに分割し、分割領域番号、空間位置（ｉ，ｊ）、時刻ｔに対応して各セルの耕地属性データＸｔ（ｉ，ｊ）を属性データ格納部に記憶する領域分割ステップと、
処理部が、前記属性データ格納部を参照し、対象とする分割領域について時刻ｔにおける位置（ｉ，ｊ）のセルの耕地属性データＸｔ（ｉ，ｊ）を読み取り；  処理部が、読み取った耕地属性データＸｔ（ｉ，ｊ）自身から、ベクトル量ＶＸｔ（ｉ，ｊ）＝｛Ｘｔ−ｎ（ｉ−ｋ，ｊ−ｍ）｜０≦ｎ≦ｄｎ，−ｄｉ≦ｋ≦ｄｉ’，−ｄｊ≦ｍ≦ｄｊ’｝（ただし、ｎ：時間ラグ、ｋ，ｍ：空間ラグ、ｄｎ：時間埋め込み次元、ｄｉ，ｄｉ’，ｄｊ，ｄｊ’：空間埋め込み次元、である。）で示される時空間埋め込みによって、Ｘｔ＋ｎ（ｉ，ｊ）＝Ｇ（ＶＸｔ（ｉ，ｊ））  を満たす局所ダイナミクスを表現記述するＧを、前記対象とする分割領域の耕地属性データに対して最適化及び再構成するモデル非依存法により同定し；  処理部が、同定された局所ダイナミクスを表現記述するＧを特定する各パラメータを、局所ダイナミクス識別子及び分割領域番号に対応してパラメータ格納部に記憶する局所ダイナミクス演算ステップと
を含む時空間ダイナミクスの同定方法。
【請求項２】
前記耕地属性データは、作物、収量、雑草量、水分、肥沃度、有機物、ｐＨ、土壌硬度、土壌分光特性、植物体分光特性、軟弱度等の耕地の時間的空間的な変動が生じる耕地属性であることを特徴とする請求項１に記載の時空間ダイナミクスの同定方法。
【請求項３】
パラメータ格納部には、局所ダイナミクスに影響するシステムパラメータがさらに記憶され、
前記局所ダイナミクス演算ステップは、システムパラメータをパラメータ格納部から読み取り、システムパラメータを用いて局所ダイナミクスを同定・再構成することを特徴とする請求項1又は２に記載の時空間ダイナミクスの同定方法。
【請求項４】
局所ダイナミクスの誤差値が所定範囲に収束しない場合、前記領域分割ステップによる分割領域の再設定又は前記局所ダイナミクス演算ステップによるパラメータの再設定を実行し、前記局所ダイナミクス演算ステップにより、再度局所ダイナミクスの同定・再構成を実行することを特徴とする請求項１乃至３のいずれかに記載の時空間ダイナミクスの同定方法。
【請求項５】
請求項１に記載された時空間ダイナミクスの同定方法により同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、同定・再構成された局所ダイナミクスを用い、所定時刻経過後の耕地属性データの時空間変動予測結果を演算し、予測された耕地属性データを、各分割領域番号、セル番号又は空間位置、時刻に対応して予測格納部に記憶する予測演算ステップと、
各前記格納部のいずれか又は複数を参照し、耕地空間フレーム又は分割領域の耕地属性データに関するマップを生成し、生成されたマップに関するデータを、マップ識別番号、耕地属性データ種類、分割領域番号、セル番号又は空間位置、時刻に対応してマップ格納部に記憶するマップ生成ステップと、
耕地空間フレーム又は分割領域のマップを表示又は出力する出力ステップと
を含む変動予測方法。
【請求項６】
請求項１に記載された時空間ダイナミクスの同定方法により同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、同定された局所ダイナミクスに関する各種統計量又は特徴量を演算又は求め、演算された又は求められた統計量又は特徴量を、分割領域番号、セル番号又は空間位置、時刻毎に統計量格納部に記憶する統計量演算ステップと、
統計量格納部を参照し、耕地空間フレーム又は分割領域の統計量に関するマップを生成し、生成されたマップに関するデータをマップ格納部に記憶するマップ生成ステップと、
耕地空間フレーム又は分割領域のマップを表示又は出力する出力ステップと
を含む変動予測方法。
【請求項７】
前記出力ステップでは、
（１）ある時間ステップ後の耕地属性の予測値を各セルに与え、これによって生成される空間フレームのマップ、
（２）時空間埋め込みによって得られた局所ダイナミクスから計算される統計量を領域毎に与え、これによって生成される耕地空間フレームのマップ、
（３）時空間埋め込みによって得られた局所ダイナミクスに基づいて、決定論的非線形予測を行い、その際に計算される統計量を分割領域毎に与え、これによって生成される耕地空間フレームのマップ、
（４）モデル依存法において、あらかじめユーザーによって与えられたモデルの式のパラメータを各領域毎にあたえ、これによって生成される耕地空間フレームのマップ、及び、
（５）得られた局所ダイナミクスを用いて、各領域の統計量を当該領域の特徴量として計算し、これらによって耕地空間フレーム全体を分級したマップ、
のいずれか又は複数のマップを生成することを特徴とする請求項５又は６に記載の変動予測方法。
【請求項８】
請求項１に記載された時空間変動予測方法における前記局所ダイナミクス演算ステップにより同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、前記予測演算ステップにより演算された耕地属性データの予測結果を、予測格納部から参照し、
ある時刻及び位置における耕地属性データと、将来のある時刻及び位置における希望する耕地属性データとに基づき、同定された局所ダイナミクスを用いて、耕地属性データの変動を制御することを特徴とする変動制御方法。
【請求項９】
処理部が、予め測定された耕地属性の時空間分布データを入力部により入力し、時刻毎に耕地属性の時空間分布データを入力データ格納部に記憶する入力ステップと、
処理部が、前記入力データ格納部を参照して耕地属性の時空間分布データを読み取り、対象とする空間フレームを、複数の分割領域に分割し、さらに分割領域を各セルに分割し、分割領域番号、空間位置（ｉ，ｊ）、時刻ｔに対応して各セルの耕地属性データＸｔ（ｉ，ｊ）を、属性データ格納部に記憶する領域分割ステップと、
処理部が、前記属性データ格納部を参照し、対象とする分割領域について時刻ｔにおける位置（ｉ，ｊ）のセルの耕地属性データＸｔ（ｉ，ｊ）を読み取り；  処理部が、読み取った耕地属性データＸｔ（ｉ，ｊ）自身から、ベクトル量ＶＸｔ（ｉ，ｊ）＝｛Ｘｔ−ｎ（ｉ−ｋ，ｊ−ｍ）｜０≦ｎ≦ｄｎ，−ｄｉ≦ｋ≦ｄｉ’，−ｄｊ≦ｍ≦ｄｊ’｝（ただし、ｎ：時間ラグ、ｋ，ｍ：空間ラグ、ｄｎ：時間埋め込み次元、ｄｉ，ｄｉ’，ｄｊ，ｄｊ’：空間埋め込み次元、である。）で示される時空間埋め込みによって、Ｘｔ＋ｎ（ｉ，ｊ）＝Ｇ（ＶＸｔ（ｉ，ｊ））  を満たす局所ダイナミクスを表現記述するＧを、前記対象とする分割領域の耕地属性データに対して最適化及び再構成するモデル非依存法により同定し；  処理部が、同定された局所ダイナミクスを表現記述するＧを特定する各パラメータを、局所ダイナミクス識別子及び分割領域番号に対応してパラメータ格納部に記憶する局所ダイナミクス演算ステップと
をコンピュータに実行させるための時空間ダイナミクスの同定プログラムを記録したコンピュータ読取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、時空間ダイナミクスの同定方法、変動予測方法及び変動制御方法、時空間ダイナミクスの同定プログラムを記録したコンピュータ読取り可能な記録媒体に係る。本発明は、特に、農業、林業、園芸、酪農、牧畜業等、施設園芸、施設農業、植物工場、養殖業における森林属性、草地属性等を含む耕地属性（耕地パラメータ）の空間分布について、その時間変動（時空間変動）の発生メカニズム及びダイナミクスを計測された耕地属性時空間データセットから同定し、これに基づく耕地属性の時空間変動の将来予測、対象耕地空間の領域分級、時空間制御を実現する耕地属性の時空間ダイナミクスの同定方法、変動予測方法及び変動制御方法、時空間ダイナミクスの同定プログラムを記録したコンピュータ読取り可能な記録媒体に関する。
【０００２】
ここで、耕地属性とは、例えば、収量、雑草量、水分、肥沃度、有機物、ＥＣ，ｐＨ，ｐＦ等の時空間的に変動する耕地の状態をいう。耕地属性の他にも、本発明は様々な状態変数に適用することができる。状態変数としては、例えば、収量・雑草・水分・肥料成分・病害虫分布・ＰＨ・スペクトル画像などおよびそれらのオブザーバ等である。なお、発明の詳細な説明では、主に、耕地属性について説明するが、本発明は、耕地生態系時空間ダイナミクスだけではなく、例えば、土壌パラメータ分布、森林生態系、草地生態系、顕微鏡画像内での細胞増殖、海洋・淡水生態系、ハウス・植物工場・閉鎖生物生産空間の生態系と気象変動、各種流体、細菌増殖の時空間パターンなどにも適応可能である。
【０００３】
【従来の技術】
近年、精密農業やリモートセンシングの発展によって、耕地パラメータのマップ作成が極めて容易に行われるようになり、大量のマップデータが生産されている。一方、これらの空間／時空間データに埋め込まれた生態学的・物理学的・力学的意味の理解は重要な課題である。
【０００４】
従来においては、耕地属性の空間分布データを取得し、単数または複数の耕地属性の空間分布データの重ね合わせにより、耕地の分級、予測を行っていた。従来は、耕地のある面積（領域）においてブロック化された各ブロック毎の耕地属性の時間的変動の解析が行われていたのみで、空間的相互作用を組み込んだ解析方法は未だ提案されていなかった。
【０００５】
【発明が解決しようとする課題】
しかしながら、従来の技術による分級や予測は、必ずしも十分な成果を挙げていない。その原因としては、実在の現象は、バックグラウンドとノイズと外乱と局所ダイナミクス等から構成されていることによる。従来では、特に局所ダイナミクスに基づかないため、現象発生のメカニズムが不明なままであり、必ずしも適当な予測結果を得ることができなかった。
【０００６】
本発明は、以上の点に鑑み、局所ダイナミクスに基づいて、従来では実現できなかった耕地属性の時空間変動ダイナミクスの同定、分級、予測、耕地属性の制御等を高精度で実現することを目的とする。本発明は、耕地属性のように時間的・空間的変動を引き起こすダイナミクスを局所的に同定し、得られたダイナミクスを用いて、耕地の分級、耕地の状態の将来予測、対象耕地の時空間制御を行うことを目的とする。
【０００７】
本発明は、局所ダイナミクスの同定に各種の手法を適用し、また均一ではない局所ダイナミクスに対し、領域分割の最適化で対応し、また、局所ダイナミクスの空間的クラスタリングを行うようにした、特に耕地属性についての時空間ダイナミクスの同定方法その変動予測方法及びその制御方法を提供することを目的とする。
【０００８】
【課題を解決するための手段】
本発明の第１の解決手段によると、
属性の時空間分布を入力し、時刻毎に属性の空間フレームを入力データ格納部に記憶する入力ステップと、
入力データ格納部を参照し、対象とする空間フレームを、複数の分割領域に分割し、さらに分割領域を各セルに分割し、各セルの属性データを、分割領域番号、セル番号又は空間位置、時刻に対応して属性データ格納部に記憶する領域分割ステップと、
属性データ格納部を参照し、あるセルについて、そのセルの周囲に近接する局所近傍のセルによって局所近傍グリッドを構成し、ある時刻における各セルの属性を、近接過去又は過去のある時刻における当該セル及びその近傍のセルの属性による写像として求めることにより、分割領域の局所ダイナミクスを同定し、局所ダイナミクスを特定する各パラメータを、局所ダイナミクス識別子及び分割領域番号に対応してパラメータ格納部に記憶する局所ダイナミクス演算ステップとを含む時空間ダイナミクスの同定方法及び時空間ダイナミクスの同定プログラムを記録したコンピュータ読取り可能な記録媒体を提供する。
【０００９】
本発明の第２の解決手段によると、
上述の時空間ダイナミクスの同定方法により同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、同定・再構成された局所ダイナミクスを用い、所定時刻経過後の属性データの時空間変動予測結果を演算し、予測された属性データを、各分割領域番号、セル番号又は空間位置、時刻に対応して予測格納部に記憶する予測演算ステップと、
各前記格納部のいずれか又は複数を参照し、耕地空間フレーム又は分割領域の属性データに関するマップを生成し、生成されたマップに関するデータを、マップ識別番号、属性データ種類、分割領域番号、セル番号又は空間位置、時刻に対応してマップ格納部に記憶するマップ生成ステップと、
耕地空間フレーム又は分割領域のマップを表示又は出力する出力ステップと
を含む変動予測方法を提供する。
【００１０】
本発明の第３の解決手段によると、
上述の時空間ダイナミクスの同定方法により同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、同定された局所ダイナミクスに関する各種統計量又は特徴量を演算又は求め、演算された又は求められた統計量又は特徴量を、分割領域番号、セル番号又は空間位置、時刻毎に統計量格納部に記憶する統計量演算ステップと、
統計量格納部を参照し、耕地空間フレーム又は分割領域の統計量に関するマップを生成し、生成されたマップに関するデータをマップ格納部に記憶するマップ生成ステップと、
耕地空間フレーム又は分割領域のマップを表示又は出力する出力ステップと
を含む変動予測方法を提供する。
【００１１】
本発明の第４の解決手段によると、
請求項１に記載された時空間変動予測方法における前記局所ダイナミクス演算ステップにより同定された局所ダイナミクスを特定する各パラメータをパラメータ格納部から参照し、前記予測演算ステップにより演算された属性データの予測結果を、予測格納部から参照し、
ある時刻及び位置における属性データと、将来のある時刻及び位置における希望する属性データとに基づき、局所ダイナミクスを特定するシステムパラメータを求め、そのシステムパラメータにより属性データを制御することを特徴とする変動制御方法を提供する。
【００１２】
【発明の実施の形態】
まず、耕地の時空間分布データについて説明する。
図１に、空間データについての説明図を示す。ここでは、耕地空間クレーム（Ｓ）を、耕地属性の時空間ダイナミクスの同定の対象とする。耕地空間フレーム（Ｓ）は、Ｑ個の分割領域（Ｐ）に分割される。さらに、分割領域（Ｐ）は、それぞれセル（Ｃ）に分割される。あるセル（図中、黒四角）について、その周囲の近接する局所近傍（近接近傍）のセルによって局所近傍グリッド（Ｎ）を構成する。各セルのある時刻t、位置（ｉ，ｊ）における状態Ｘｔ（ｉ、ｊ）は、近接過去又は過去のある時刻における当該セルの状態Ｘｔ−１（ｉ、ｊ）と、そのセルの局所近傍セルの状態Ｘｔ−１（ｉ−ｍ，ｊ−ｎ），  （−ｄｉ≦ｍ≦ｄｉ’，  −ｄｊ≦ｎ≦ｄｊ’）の非線形写像である局所ダイナミクスＦとして記述する。すなわち、局所ダイナミクスとは、領域内の各セルのある時刻の状態を、それに先立つ時刻の局所近傍グリッド（Ｎ）のセルの状態により支配するダイナミクスである。セルの状態は、即ち耕地属性データＸｔ（ｉ、ｊ）を示し、局所ダイナミクスＦを用いて次式のように表わされる。
Ｘｔ（ｉ，ｊ）＝Ｆ（Ｘｔ−１（ｉ−ｍ，ｊ−ｎ），（−ｄｉ≦ｍ≦ｄｉ’，  −ｄｊ≦ｎ≦ｄｊ’）      （１）
【００１３】
図２に、簡略化した局所ダイナミクスについての説明図を示す。ここでは、近接過去の時刻ｔ−１における近接するセルの状態（耕地属性）のみから、ある時刻ｔのセルの状態（耕地属性）が支配される場合を示す。このとき、（１）式で、
ｄｉ＝ｄｉ’＝ｄｊ＝ｄｊ’＝１となり、（１）式は次式となる。
Ｘｔ（ｉ，ｊ）＝Ｆ（Ｘｔ−１（ｉ＋ｍ，ｊ＋ｎ），（−１≦ｍ，ｎ≦１））  （２）
【００１４】
ここで局所ダイナミクス（局所力学系）の同定の一例について述べる。まず、２次元ＣＯＤＭ（Coupled Ordinary Equation）モデルについて説明する。今、２種ｘ，ｙの植物の空間内で競争を考える。また、空間をｉ×ｊ個のセルに分割する。ここで、Ｘｎ（ｉ、ｊ）を時刻ｎにおけるセル（ｉ、ｊ）のＸ種の成長量（質量）と考えると、次式によって、植物群落の時空間成長過程がモデル化できる。
【００１５】
【数１】
0003776719.tif 000002
【００１６】
ここで、位置（ｉ，ｊ）のセルにおける成長ｄｘｉ,ｊ／ｄｔは、単位時間前の周囲のセルの関数として定義される。但し、一例として、α＝１．５，β＝１．１，Ｋｘ＝３００．０，Ｋｙ＝２５０．０；ｂｍ，ｋ＝｛１．０（ｍ＝ｋ＝０），０．２（ｍ≠ｋ），０．０（ｍ＝ｋ≠０）｝，ｒｘ＝０．０２１，ｒｙ＝０．０２２  である。
よって、（１）式又は（２）式の関係を満たすような写像Ｆをすべてのセルに関して満足するように決定してやればよい。Ｆの同定法は次の２者に大別される。即ち、写像の関数形がユーザから明示的に与えられる場合をモデル依存法、データ自身から時空間埋込みによって再構成される場合をモデル非依存法と呼ぶ。
【００１７】
図３は、本発明に関する耕地時空間ダイナミクス同定・予測システムの構成図を示す。
このシステムは、処理部１０、記憶部２０、入力部３０、出力部４０を備える。処理部１０は、領域分割部１１、局所ダイナミクス演算部１２、統計量演算部１３、予測演算部１４、マップ生成部１５を有する。また、記憶部２０は、入力データ格納部２１、パラメータ格納部２２、統計量格納部２３、予測格納部２４、マップ格納部２５、属性データ格納部２６を有する。
【００１８】
入力部３０は、測定された耕地属性の時空間分布データをシステムへ取得する。入力データ格納部２１は、入力部３０により入力された耕地属性の空間分布データを、時刻毎に耕地空間フレーム番号に対応して、記憶する。領域分割部１１は、入力データ格納部２１を参照し、対象とする耕地空間フレーム（Ｓ）を、複数の分割領域（Ｐ）に分割し、さらに分割領域（Ｐ）を各セル（Ｃ）に分割する。各セル（Ｃ）の状態、即ち耕地属性データは、時刻、分割領域番号、セル番号、空間位置に対応して属性データ格納部２６に記憶される。局所ダイナミクス演算部１２は、属性データ格納部２６を参照し、分割された分割領域（Ｐ）の局所ダイナミクス（Ｆ）を同定する。パラメータ格納部２２は、同定した局所ダイナミクスを示す各式のパラメータデータを、局所ダイナミクス識別子、分割領域番号、属性データ種別に対応して記憶する。また、パラメータ格納部２２は、局所ダイナミクスの同定に影響するシステムパラメータとして、例えば、気温・降水量・日照・土壌タイプ・バックグラウンド等を記憶する。
【００１９】
予測演算部１４は、パラメータ格納部２２を参照し、同定された局所ダイナミクス（Ｆ）を用い、耕地属性の時空間データの予測結果を演算する。予測格納部２４は、予測演算部１４により演算された耕地属性データの予測結果を、各分割領域番号、セル番号、空間位置、時刻及び空間位置に対応して記憶する。統計量演算部１３は、同定された局所ダイナミクス（Ｆ）に基づいて各種統計量を演算する。統計量格納部２３は、統計量演算部１３により演算された統計量を、各分割領域番号、セル番号、時刻、空間位置毎に記憶する。
【００２０】
マップ生成部１５は、記憶部１０内の各格納部２１〜２５を適宜参照し、パラメータデータ、統計量、予測結果等に基づき、分割領域（Ｐ）の耕地属性に関するマップを生成する。マップ格納部２５は、マップ生成部１５により生成されたマップに関するデータを、各マップ識別番号、属性データ種類、分割領域番号、セル番号、時刻及び空間位置に対応して記憶する。出力部４０は、マップ格納部２５を参照し、耕地空間フレーム（Ｓ）のマップを表示又は出力する。また、出力部４０は、記憶部１０の各格納部を参照し、パラメータ、予測結果、統計量、分級マップ（Ｉ）等の中間結果や最終結果を適宜表示及び出力する。
【００２１】
つぎに、図４に、本発明に関する耕地時空間ダイナミクスの同定方法、変動予測方法及び変動制御方法の動作についてのフローチャートを示す。本発明では、分割領域（Ｐ）毎に局所ダイナミクスＦを同定することができる。さらに、複数の分割領域（Ｐ）に対する広域ダイナミクスを、求めた複数の局所ダイナミクスＦに基づき求めるようにしてもよい。
【００２２】
まず、入力部３０は、耕地属性データＸｔ（ｉ、ｊ）の時間空間分布データセットを、フロッピーディスク、ＣＤ−ＲＯＭ，  ＭＤ，  Ｚｉｐファイル、インターネットなどの適宜の媒体からを経て入力し、そのデータを入力データ格納部２１に格納する（Ｓ１０１）。領域分割部１０は、各種の媒体によって供給される耕地属性の時空間データセットから、離散時刻ｔ毎に構成される耕地空間フレーム（Ｓ）を構成し、領域を分割して分割領域（Ｐ）を得て、さらに、分割領域（Ｐ）をセル（Ｃ）に分割する（Ｓ１０５）。領域分割部１０は、分割領域（Ｐ）に関する耕地属性データＸｔ（ｉ，ｊ）  （  ０≦ｔ≦ｎｔ，  ０≦ｉ≦ｎｉ，  ０≦ｊ≦ｎｊ）を、属性データ格納部２６に記憶する。局所ダイナミクス演算部１２は、属性データ格納部２６から耕地属性データを読み取り、そのデータに対して、上述の（１）式又は（２）式の局所ダイナミクスＦを同定する（Ｓ１０７）。そのための方法は、主に、モデル非依存的方法とモデル依存的方法に大別される。
【００２３】
まず、モデル非依存的方法について説明する。このモデルでは、対象とする分割領域（Ｐ）の耕地属性データＸｔ（ｉ，ｊ）のみから局所ダイナミクスＦを同定・再構成する。局所ダイナミクスＦの同定方法としては、近接非線形写像法、空間埋め込み法、時空間埋め込み法、ヤコビ行列の推定、カオス、カオスニュ−ロ、ファジー、ファジーニューロ、ＲＳＭなどがある。また、モデル非依存法では、時空間埋め込み法により、得られた局所ダイナミクスＦによって、直接リアプノフ指数やフラクタル次元等の計算や、決定論的予測を実施できる。
【００２４】
また、モデル依存法では、局所ダイナミクスを記述するモデルをあらかじめ与える。このモデルは、例えば、微分方程式、偏微分方程式、差分方程式、セルオートマトン、ＣＭＬ（Coupled Map Lattice）、ファジーロジック、記号力学系などで与えられる。各分割領域（Ｐ）に対して領域内の全局所近傍グリッドに対して最適化する。ここで得られた各式パラメータをパラメータ格納部２２に格納しておき、これらを各分割領域（Ｐ）の特徴量として考える。マップ生成部１５では、これら特徴量によってマップを形成できる。また、モデルの方程式と、パラメータをパラメータ格納部２２に格納しておけば、随時に局所ダイナミクスＦを構成できる。これによって、時間経過に対して順方向にシミュレーションを行うことによって予測を行うことができる。これらモデル依存的方法・非依存的方法及び、局所ダイナミクスＦの同定方法は、公知の適宣の方法を採用することができる。局所ダイナミクス演算部１２は、同定された局所ダイナミクスＦを特定するためのパラメータをパラメータ格納部２２に記憶する。
【００２５】
ここで、一例として、近接非線形写像法についてのニューラルネットワークによるものについて説明する。
図５は、ニューラルネットワークによる局所ダイナミクス同定法の説明図である。領域内の時刻ｔ及びｔ＋１の局所近傍グリッド５１及び５３について、（１）式又は（２）式の局所ダイナミクスＦをニューラルネットワーク５２によって構成する。ニューラルネットワーク５２は、一例として、マルチレイヤーニューラルネットとする。対象セルの位置を（ｉ，ｊ）  として、入力層に耕地属性データＸｔ（ｉ＋ｍ，ｊ＋ｎ）、（−１≦ｍ，ｎ≦１）という９個の状態変数を入力する。またノイマン近傍を選択する場合には、｛Ｘｔ（ｉ−１，ｊ），  Ｘｔ（ｉ，ｊ），  Ｘｔ（ｉ＋１，ｊ），  Ｘｔ（ｉ，ｊ−１），  Ｘｔ（ｉ，ｊ＋１）｝という５個の状態変数の入力となる。分割領域（Ｐ）のサイズをＫ×Ｌセル個とすると、Ｋ×Ｌ個以下の局所近傍グリッド（Ｎ）を設定できる。
【００２６】
このようなニューラルネットワーク５２において、Ｋ×Ｌ対（若しくはそれ以下）の教師データＸｔ＋１（ｉ，ｊ）と入力データＸｔ（ｉ，ｊ）＝｛Ｘｔ（ｉ＋ｍ，ｊ＋ｎ）｜−１≦ｍ，ｎ≦１｝  によって学習を行う。誤差が設定値ｅを下回った場合に学習を終了し、学習によって得られた、重み係数、オフセット等からなるニューラルネットワーク５２のパラメータをベクトルデータ（特徴量）として分割領域（Ｐ）ごとに抽出する。なお、誤差を小さくするために領域分割（ステップＳ１０５）又は、パラメータ演算（ステップＳ１０７内）のいずれか又は両方の処理フィードバックをかけ、再分割することができる。
【００２７】
つぎに、時空間埋め込み法について説明する。対象とする分割領域（Ｐ）に対して、時空間埋め込みを（３）式に示す。ここで、「ＶＸ」は、ベクトル量である（以下同様）。
ＶＸｔ（ｉ，ｊ）＝｛Ｘｔ−ｎ（ｉ−ｋ，ｊ−ｍ）｜０≦ｎ≦ｄｎ，−ｄｉ≦ｋ≦ｄｉ’，−ｄｊ≦ｍ≦ｄｊ’｝      （３）
時間ラグ：ｎ
空間ラグ：ｋ，ｍ
ｄｎ：時間埋め込み次元
ｄｉ，ｄｉ’，ｄｊ，ｄｊ’：空間埋め込み次元
局所ダイナミクスを下記のＧで表現記述し、これをデータに対して最適化する。
【００２８】
Ｘｔ＋ｎ（ｉ，ｊ）＝Ｇ（ＶＸｔ（ｉ，ｊ））      （４）
（３）式が時空間埋め込みの操作である。（４）式を満たすＧを求めれば局所ダイナミクスが得られる。この際各種推定法が用いられる。
つぎに、統計量演算部１３は、これら各種方法によって得られた局所ダイナミクスＦから統計量、Ｌｙａｐｕｎｏｖ指数（λ）、相関次元（ＤＣ  ）等を計算する（Ｓ１０９）。各分割領域（Ｐ）毎に、各種統計量を割り当てて計算することができる。また、統計量演算部１３は、計算された各種統計量を統計量格納部２３に記憶する。
【００２９】
一方、予測演算部１４では、（１）ある初期値による将来予測と、（２）決定論的非線形予測等により予測・演算を行う。
予測演算部１４は、「（１）ある初期値による将来予測」については、以下のように行う。得られた局所ダイナミクスＦを用いて、ある時刻までの領域内のセルの状態（耕地属性）から、それ以降の時刻のセルの状態（耕地属性）を推定して予測することができる。（Ｓ１１０）すなわち、予測演算部１４は、パラメータ格納部２２からパラメータを読み込み、パラメータにより特定される局所ダイナミクスＦを用いて、各セルに任意の初期値を与え、耕地属性の空間分布データの時間変化を予測する。予測演算部１４は、耕地属性データに関する予測結果を耕地空間フレーム（Ｓ）に対応して生成し、予測耕地空間フレームを作成し、その結果を予測格納部２４に出力する。
【００３０】
また、決定論的非線形予測等により予測誤差の各種統計量（相関係数、２乗平均誤差（ｒｍｓ）、相対２乗平均誤差（ｒｒｍｓ）、ＫＳエントロピー）を求めことも可能である。予測演算部１４は、「（２）決定論的非線形予測」については次のように行う。予測演算部１４及び統計量演算部１３の処理に関し、決定論的非線形予測の概要を（１）式又は（２）式のモデルを対象に例示する。
【００３１】
図６に、決定論的非線形予測の説明図を示す。
まず、時刻ｔ位置（ｉ，ｊ）における状態ベクトルは、
ＶＸｔ（ｉ，ｊ）＝｛Ｘｔ（ｉ＋ｍ，ｊ＋ｎ）｜（−１≦ｍ，ｎ≦１）｝      （６）
である。また、既に同定された局所ダイナミクスＦと（１）式より、
Ｘｔ＋１１（ｉ，ｊ）＝Ｆ１（ＶＸｔ（ｉ，ｊ））      （７）
となる。ここで、左辺の上付１は、時間ステップ１先の予測値であることを示す。また、Ｆの上付１は予測を計算を一回行ったことを示す。
【００３２】
（７）式により、位置（ｉ，ｊ）の局所近傍グリッドの各セルの値（状態ベクトルＶＸｔ（ｉ，ｊ））によって、時刻ｔ＋１における（ｉ，ｊ）のセルの状態値を予測できる。時刻ｔ＋１における局所近傍グリッド内のセルに対しても同様の操作を行うことによって、グリッド内全ての状態値を予測する事ができる（注：この操作は領域全体について実施するため、空間フレーム内の全てのセルの状態を予測できる。換言すれば、空間フレーム内の全てのセルの状態を更新することができる。）。これにより、時刻ｔ＋１位置（ｉ，ｊ）における状態ベクトルを、
【００３３】
0003776719.tif 000003
のように構成できる。
【００３４】
つぎに、既に同定された局所ダイナミクスＦと（１）式より、
Ｘｔ＋２１（ｉ，ｊ）＝Ｆ１（ＶＸｔ＋１（ｉ，ｊ））      （９）
となる。このプロセスを繰り返すと、
Ｘｔ＋ｋ１（ｉ，ｊ）＝Ｆ１（ＶＸｔ＋ｋ−１（ｉ，ｊ））      （１０）
となり、任意時間ステップｋ先の予測値が得られる。
【００３５】
つぎに、これをより一般化した方法について以下に説明する。
時刻ｔ位置（ｉ，ｊ）における状態ベクトルは、
ＶＸｔ（ｉ，ｊ）＝｛Ｘｔ（ｉ＋ｍ，ｊ＋ｎ）｜（−１≦ｍ，ｎ≦１）｝      （１１）
である。局所ダイナミクスＦ’を下記のように定義する。
ＶＸｔ＋１（ｉ，ｊ）＝Ｆ’（ＶＸｔ（ｉ，ｊ））
既に同定された局所ダイナミクスＶＦより、
ＶＸｔ＋１１（ｉ，ｊ）＝Ｆ’１（ＶＸｔ（ｉ，ｊ））
ここで、左辺の上付１は、時間ステップ１先の予測値であることを示す。また、Ｆの上付１は予測を計算を一回行ったことを示す。
【００３６】
同様に、
ＶＸｔ＋２１（ｉ，ｊ）＝Ｆ’１（ＶＸｔ＋１（ｉ，ｊ））＝Ｆ’２（ＶＸｔ（ｉ，ｊ））
これを一般化すると、
ＶＸｔ＋ｋｋ（ｉ，ｊ）＝Ｆ’ｋ（ＶＸｔ（ｉ、ｊ））
ただし、Ｘｔ＋ｋｋは、時刻ｔにおける状態ベクトルＶＸｔ（ｉ、ｊ）より、局所ダイナミクスによる予測をｋ回フィードバックを繰返し演算した耕地属性の予測値を示す。
【００３７】
このように、予測演算部１４により、（５）〜（８）式の計算を繰り返すことにより、時刻ｔ＋ｋでの各セルの耕地属性の予測値が計算され、推定される。また、統計量演算部１３により、計算された予測値と実際値との相関係数、２乗平均誤差、相対２乗平均誤差などが予測性能の統計量として計算される。
【００３８】
なお、予測の評価には、自己予測と非自己予測の２者がある。自己予測では、局所ダイナミクスの同定に用いたデータセットに対して予測を実施する。非自己予測では、局所ダイナミクスの同定に用いた
データセットとは異なるデータセットに対して予測を実施する。
フローチャートに戻り、次にマップ生成部１５は、局所ダイナミクス演算部１２により求められた局所ダイナミクスとパラメータ、統計量演算部１３により求められた各種統計値、予測演算部１４により求められた、予測結果等の中間結果及び最終結果について各種マップを生成する。（Ｓ１１１）その際、必要な情報を記憶部２０の各格納部から適宜読み出し、生成された各種マップをマップ格納部２５に記憶する。
【００３９】
マップ生成部１５によりに生成されるマップの項目としては、例えば、次のようなものが挙げられる。
１）耕地属性のｋ時間ステップ後の予測値を各セルに与え、これによって生成される耕地空間フレームのマップ
２）時空間埋め込みによって得られた局所ダイナミクスから計算される統計量（リアプノフ指数、フラクタル次元、情報エントロピーなど）を領域毎に与え、これによって生成される耕地空間フレームのマップ
３）時空間埋め込みによって得られた局所ダイナミクスに基づいて、決定論的非線形予測を行い、その際に計算される統計量（相関係数、２乗平均誤差、相対２乗平均誤差など）を分割領域（Ｐ）毎に与え、これによって生成される耕地空間フレームのマップ
４）モデル依存法において、あらかじめユーザーによって与えられたモデルの式のパラメータを各領域毎にあたえ、これによって生成される耕地空間フレームのマップ
５）得られた局所ダイナミクスを用いて、各領域のもつ統計量を当該領域の特徴量として計算し、これらによって耕地空間フレーム全体を分けた分級マップ
【００４０】
つぎに、マップ格納部２５に格納された各種マップは、出力部４０によって、各種ディスプレイにより表示されたり、ハードコピー、各種媒体、インターネット等へ出力される（Ｓ１１３）。
【００４１】
図７に、予測マップの出力例の説明図を示す。この図は、決定論的非線形予測によって得られたものである。図７（Ａ）は、ある時刻から２０日目の植生分布の実際値、図７（Ｂ）は、１０日目植生分布からの２０日目植生分布の予測値を示す。両者のパターンは、ほぼ一致していることがわかる。
【００４２】
図８に、予測ステップに対応する相対２乗平均誤差及び相関関係について示した図を示す。予測ステップｋが増えるに従って、相関関数ＣＣが増加し、相対２乗平均誤差ＲＲＭＳＥが減少している。
なお、本発明の時空間ダイナミクスの同定方法により同定された局所ダイナミクスを用い、属性データの変動を制御する変動制御方法に、本発明を適用することができる。局所ダイナミクスの同定に影響するシステムパラメータとしては、例えば、気温・降水量・日照・土壌タイプ・バックグラウンド等がある。教師データとして各種の耕地属性データを入力することで局所ダイナミクスを求め、得られた局所ダイナミクスを基に、ある時刻ｔの領域内の各セルの状態（耕地属性）と、将来時刻ｔ＋Ｔの希望のセルの状態（耕地属性）を入力し、各種の時空制御論を利用して、それ以降のセルおよび領域を希望のセルの状態（耕地属性）にするためのこれらシステムパラメータの操作量を計算することもできる。
【００４３】
また、データ入力の供給としては、例えば、衛星画像・航空機画像・近接画像等を用いることができる。また、領域分割としては、例えば、確定分割法、自動分割法を用いることができる。自動分別法には、対称型、非対象型（ローカル最適化、グローバル最適化）がある。
【００４４】
本発明の時空間ダイナミクスの同定方法、変動予測方法及び変動制御方法は、時空間ダイナミクスの同定プログラム、変動予測及び変動制御方法プログラムを記録したコンピュータ読み取り可能な記録媒体、時空間ダイナミクスの同定プログラム、変動予測及び変動制御プログラムを含みコンピュータの内部メモリにロード可能なプログラム又はプログラム製品、時空間ダイナミクスの同定プログラム、変動予測及び変動制御プログラムを含みコンピュータが使用可能な記録媒体にストアされたプログラム又はプログラム製品等により提供されることができる。
【００４５】
【発明の効果】
本発明によると、以上のように、局所ダイナミクスに基づいて、従来では実現できなかった耕地属性の時空間変動ダイナミクスの同定、分級、予測、耕地属性の制御等を高精度で実現することができる。本発明によると、耕地属性のように時間的・空間的変動を引き起こすダイナミクスを局所的に同定し、得られたダイナミクスを用いて、耕地の分級、耕地の状態の将来予測、対象耕地の時空間制御を行うことができる。
【００４６】
本発明によると、力学系の同定に各種の手法を適用し、また均一ではない力学系に対し、領域分割の最適化で対応し、また、力学系の空間的クラスタリングを行うようにした、特に耕地属性についての時空間ダイナミクスの同定方法その変動予測方法及びその制御方法を提供することができる。
【図面の簡単な説明】
【図１】空間データの関係についての説明図。
【図２】簡略化した局所ダイナミクスについての説明図。
【図３】本発明に関する耕地時空間ダイナミクス同定・予測システムの構成図。
【図４】本発明に関する耕地属性時空間ダイナミクスの同定システムおよびそれに基づく耕地分級システム、耕地属性時空間変動予測システムの動作についてのフローチャート。
【図５】ニューラルネットワークによる局所ダイナミクス同定法の説明図。
【図６】決定論的非線形予測の説明図。
【図７】予測マップの出口側の説明図。
【図８】予測ステップに対応する相対２乗平均誤差及び相関関係について示した図。
【符号の説明】
１０  処理部
１１  領域分割部
１２  局所ダイナミクス
１３  統計量演算部
１４  予測演算部
１５  マップ生成部
２０  記憶部
２１  入力データ格納部
２２  パラメータ格納部
２３  統計量格納部
２４  予測格納部
２５  マップ格納部
２６  属性データ格納部
３０  入力部
４０  出力部
【図１】
0003776719.tif 000004
【図２】
0003776719.tif 000005
【図３】
0003776719.tif 000006
【図４】
0003776719.tif 000007
【図５】
0003776719.tif 000008
【図６】
0003776719.tif 000009
【図７】
0003776719.tif 000010
【図８】
0003776719.tif 000011