構造可変な飛翔体のテレメトリデータ解析装置及びテレメトリデータ解析プログラム
【課題】飛翔体のテレメトリデータからカオス力学系に基づく異常疑のある動作を検知するテレメトリデータ解析装置を提供する。
【解決手段】飛翔体のテレメトリデータ群から時系列データ群を切り出すデータ切出部と、時系列データ群の埋め込み次元を算出するデータ前処理部と、時系列データ群からＲＳ及びＦＴサロゲートデータ群を作成するサロゲートデータ作成部と、時系列データ群とＲＳ及びＦＴサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、非線形統計量に関する有意性検定を行い、時系列データ群がカオス的であるか否かを判定するデータ検定部と、カオス的であると判定された時系列データ群の非線形統計量と過去の時系列データ群の非線形統計量との差から、飛翔体の異常疑のある動作を検出する異常疑検出部と、データ前処理部に対して、飛翔体の構造が変化する度に、時系列データ群の埋め込み次元を算出させる制御部とを備える。
【解決手段】飛翔体のテレメトリデータ群から時系列データ群を切り出すデータ切出部と、時系列データ群の埋め込み次元を算出するデータ前処理部と、時系列データ群からＲＳ及びＦＴサロゲートデータ群を作成するサロゲートデータ作成部と、時系列データ群とＲＳ及びＦＴサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、非線形統計量に関する有意性検定を行い、時系列データ群がカオス的であるか否かを判定するデータ検定部と、カオス的であると判定された時系列データ群の非線形統計量と過去の時系列データ群の非線形統計量との差から、飛翔体の異常疑のある動作を検出する異常疑検出部と、データ前処理部に対して、飛翔体の構造が変化する度に、時系列データ群の埋め込み次元を算出させる制御部とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ロケットなどの構造可変な飛翔体から送信されるテレメトリデータを解析して、飛翔体の異常疑のある動作を検出するテレメトリデータ解析装置及びテレメトリデータ解析プログラムに関する。
【背景技術】
【０００２】
飛行中のロケット、軌道上の人工衛星、無人宇宙船や有人宇宙船の状態は、これらの飛翔体から送信されるテレメトリデータを通して、地上の観測基地において常時観測されている。テレメトリデータとは、飛翔体の状態を表す各種の掲載機器から送信される時系列データのことであり、主なものとしては、各種搭載機器の電圧や電流、動作のオンオフ状態、飛翔体の各部の温度や圧力、送受信装置の送受信電力、飛翔体の外部環境に関する情報などのデータがある。これらのデータは各々、飛翔体に搭載されたテレメータ送信機のエンコーダにて収集され、一定のフォーマットにて地上の観測基地に送信される。図１０に示すように、観測基地Ｂにおいては、飛翔体（ロケットＲ）から送信された各種テレメトリデータＳ１〜ＳＮをアンテナ等の受信装置Ａ（Ａ１〜ＡＮ）にて受信した後、信号号復調変換装置Ｄ（Ｄ１〜ＤＮ）にて復調し、Ａ／Ｄ変換や各種の物理量変換を行う。そして、テレメトリデータ解析装置Ｚにおいて、変換された時系列データを統計解析や周波数解析したり、データの可視化を行い表示装置Ｐに表示させたりして、その結果を過去のデータ、設計値、又は予測値などと比較する。このような時系列データ解析により、飛翔体が計画通りに飛行できたことを確認したり、飛翔体の異常な動作を検知したり、また、飛翔体に何らかの事故が発生した場合にはその原因を調査して搭載機器の改良に役立てたりする（特許文献１）。
【０００３】
このような従来のテレメトリデータ解析（時系列データ解析）では、解析対象となるテレメトリデータを過去のデータ、設計値、又は予測値と比較したときに差異があまり大きくない場合には、その差異は通常ランダムノイズとして認識される。
【０００４】
しかし、近年数理物理学の分野において、一見するとランダムに見える時系列データがカオス力学系における決定論的な法則に従っている場合があることが明らかになってきた。カオス力学系の特徴は、そのダイナミクスがある種の非線形な因果関係に基づいているため、その時系列データがランダムノイズに似た複雑な振る舞いを示す点にある。ある時系列データがカオス力学系に基づく時系列データであるか否かを判定するためには、カオス力学系の特徴量であるフラクタル次元、リャプノフスペクトルなどの非線形統計量を算出したり、カオス力学系に特有なアトラクタ軌道（ストレンジアトラクタ）を相空間にて可視化したりするといったデータ解析を行う必要がある。このような時系列データ解析を特にカオス時系列データ解析といい、現状では主に生体内の各種生理信号、プラント等の状態信号など、各種工学システムの計測データに実用されている（非特許文献１）。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】米国特許公開２００９０１０５８９２号公報
【非特許文献】
【０００６】
【非特許文献１】合原 一幸、池口 徹、山田 泰司、小室 元政著、カオス時系列解析の基礎と応用、産業図書、２０００。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
従って、ロケット等の飛翔体のテレメトリデータに対する従来のデータ解析では、テレメトリデータが一見するとランダムに見えるような時系列データであった場合には、その中に潜む非線形性に基づく決定論的な要因、つまり飛翔体の動作における何らかのカオス力学系、の存在を見逃してしまう可能性があった。このことは、従来の飛翔体のテレメトリデータ解析では、カオス力学系に基づく飛翔体の異常疑のある動作を見逃してしまうことを意味する。
【０００８】
本発明は、上記の実情を鑑みて為されたものであり、ロケットなどの構造可変な飛翔体から送信されるテレメトリデータからカオス力学系に基づく飛翔体の異常疑な動作を検知することを目的とする。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、請求項１に記載の発明は、構造可変な飛翔体のテレメトリデータ解析装置であって、構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出部と、切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理部と、算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理部と、切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成部と、リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定部と、カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出部と、前記データ前処理部に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御部と、を備えたことを特徴とする。
【００１０】
また、請求項２に記載の発明は、コンピュータを、構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出手段と、切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理手段と、算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理手段と、切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成手段と、リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定手段と、飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定手段と、カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出手段と、前記データ前処理手段に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御手段と、して機能させることを特徴とする構造可変な飛翔体のテレメトリデータ解析プログラムであることを特徴とする。
【発明の効果】
【００１１】
本発明によれば、構造可変な飛翔体からの一見するとランダムのように見えるテレメトリデータの中からカオス力学系に基づく決定論的な要因を系統だって見つけ出すことが可能となるので、従来見過ごされていた可能性のある飛翔体の異常疑のある動作を検出することができるようになる。さらに、この装置を飛行後解析に適用した場合には、既に発生してしまった飛翔体の故障の原因を効率良く特定して、将来発生する可能性のある故障を予測することが可能となるので、次号機以降の飛翔体の設計の改善を行うことができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の一実施の形態に係るテレメトリデータ解析装置において解析対象となる多段式ロケットからのテレメトリデータを示す図である。
【図２】図１に示した多段式ロケットから観測基地に送られてくるテレメトリデータの特徴を示す図である。
【図３】本発明の一実施の形態に係るテレメトリデータ解析装置のブロック図である。
【図４】図３に示したテレメトリデータ解析装置において行われるテレメトリデータ解析の全体的な流れを示すフローチャートである。
【図５】図４に示したデータの前処理の詳細を示すフローチャートである。
【図６】図４に示したデータの可視化処理の詳細を示すフローチャートである。
【図７】図４に示した非線形統計量の推定処理の詳細を示すフローチャートである。
【図８】図４に示したサロゲートデータの作成処理の詳細を示すフローチャートである。
【図９】図４に示した非線形統計量の推定処理の詳細を示すフローチャートである。
【図１０】ロケットからのテレメトリデータを受信する観測基地の主要構成要素を示すブロック図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照しながら、本発明の一実施の形態に係る構造可変な飛翔体のテレメトリデータ解析装置及びテレメトリデータ解析プログラムについて説明する。図１は、本発明の一実施の形態に係るテレメトリデータ解析装置１（図３）において解析対象となる多段式ロケットＲからのテレメトリデータを示す図である。
【００１４】
［解析対象となるテレメトリデータ］
本発明において解析対象となるのは、図１に示すように、多段式ロケットのような構造可変な飛翔体Ｒから送信されるテレメトリデータである。多段式ロケットの一例として二段式ロケットについて説明する。本ロケットは、各段Ｒ１，Ｒ２に液体酸素と液体水素によるロケットエンジンを備えた２段式の液体燃料ロケットであり、１段目Ｒ１には推力増強のための大型固体ロケットブースタが２本装備されており、２段目Ｒ２の先端にペイロードフェアリングで覆われた人工衛星Ｒ３を装備している。また、さらに推力を増強するために、１段目に固体補助ロケットが複数本装備される場合もある。
【００１５】
これらのパーツＲ１〜Ｒ３にはそれぞれ、飛行中のロケットＲ（の搭載機器）の状態をチェックするためにテレメータ送信機が装備されており、これらのテレメータ送信機は、パーツＲ１〜Ｒ３に装備された各種の搭載機器の状態を表すテレメトリデータ群ＳＧ１〜ＳＧ３をそれぞれ地上の観測基地Ｂへとリアルタイムに送信する。テレメトリデータ群ＳＧ１〜ＳＧ３はそれぞれ、Ｎ１個のテレメトリデータＳ１１〜Ｓ１Ｎ１，Ｎ２個のテレメトリデータＳ２１〜Ｓ２Ｎ２，Ｎ３個のテレメトリデータＳ３１〜Ｓ３Ｎ３を含む。これらテレメトリデータ群ＳＧ１〜ＳＧ３中のテレメトリデータの個数は、一般にはパーツごと異なっている。なお、搭載機器には非常に様々なものがあるが、通信系、計測系、制御系、タイマ・点火系、電源・管制系、配線系などに関する搭載機器はロケットＲの可動に必須である。
【００１６】
パーツＲ１，Ｒ２は、ロケットＲの飛行中に下から順に作動し、ロケットＲが必要な速度に達するための推力を出し終わったら順次本体から切り離されていき、最終的には、人工衛星Ｒ３のみが地球の周回軌道上に投入される。このようなロケットＲからのパーツＲ１，Ｒ２の切り離しに応じて、観測基地Ｂに送信されてくるテレメトリデータは、図２に示すように時間的に変化する。図２は、図１に示した多段式ロケットＲから観測基地Ｂに送られてくるテレメトリデータの特徴を示す図である。時刻ｔ０にてロケットＲが打ち上げられ、時刻ｔ１にてパーツＲ１が品体から切り離されると、時刻ｔ１以降、パーツＲ１から観測基地Ｂに送信されていたテレメトリデータ群ＳＧ１は送信されなくなる。次いで、時刻ｔ２にてパーツＲ２が本体から切り離されると、時刻ｔ２以降、パーツＲ２から観測基地Ｂに送信されていたテレメトリデータ群ＳＧ２は送信されなくなる。そして、最終的には、それ以降（現時刻ｔ３）、周回軌道に投入された人工衛星Ｒ０からのテレメトリデータ群ＳＧ３のみが観測基地Ｂに送信され続けることになる。つまり、多段式ロケットＲのような構造可変な飛翔体から送信されるテレメトリデータは、観測基地Ｂから送信されるコマンドに応じて当該飛翔体の構造が変化する度に、観測基地Ｂにて受信されるテレメトリデータ数が変化することを特徴としている。
【００１７】
［解析装置の構成］
次に、本発明の一実施の形態に係るテレメトリデータ解析装置１の構成について説明する。図３は、本発明の一実施の形態に係るテレメトリデータ解析装置１のブロック図である。テレメトリデータ解析装置１は、データ受信部１０と、データ切出部２０と、データ前処理部３０と、データ可視化処理部４０と、非線形統計量推定部５０と、サロゲートデータ作成部６０と、データ検定部７０と、異常疑検出部８０と、データ送信部９０と、入力部１００と、記憶部１１０と、制御部１２０とを備えている。
【００１８】
データ受信部１０は、観測基地Ｂ（アンテナＡ１〜ＡＮ）にて受信され、信号復調変換装置Ｄ（Ｄ１〜ＤＮ）にて復調されＡ／Ｄ変換等が行われたテレメトリデータ群を受信する。
【００１９】
データ切出部２０は、データ受信部１０にて受信したテレメトリデータ群からテレメトリデータ解析を実施する時間領域の時系列データ群を切り出す処理を行う。
【００２０】
データ前処理部３０は、データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群であった場合には、誤り近傍点法によって当該時系列データ群が描く軌道を埋め込み可能とする相空間Ｈの次元（埋め込み次元）ｍを算出する処理を行い、切り出された時系列データ群が多次元であった場合には、その次元ｎを当該時系列データ群の埋め込み次元として設定する処理を行う。
【００２１】
データ可視化処理部４０は、データ前処理部３０にて算出された埋め込み次元を有する高次元相空間Ｈ中の２点間の距離を規定する所定のノルムに基づき時系列データ群を二値化して視覚化するリカレンスプロットと、高次元相空間Ｈ中の２点における方向ベクトル間のばらつき具合を評価するノルムに基づき時系列データ群を二値可して視覚化する同方向性リカレンスプロットと、これら２つのプロットの結果に基づき同方向的近傍プロットを行い、同方向的近傍プロットの結果に基づき時系列データ群が決定論的な法則に従っているか否かを判断する処理を行う。
【００２２】
非線形統計量推定部５０は、切り出された時系列データ群の非線形統計量として、ＧＰ（Ｇｒａｓｓｂｅｒｇｅｒ−Ｐｒｏｃａｃｃｉａ）法によって当該の時系列データ群からフラクタル次元の１つである相関次元Ｄｃを推定し、また、この時系列データ群からリャプノフスペクトル及びリャプノフ指数を推定する処理を行う。
【００２３】
サロゲートデータ作成部６０は、本データ解析がリアルタイム解析ではなく飛行後解析である場合に、切り出された時系列データ群からＲＳ（Ｒａｎｄｏｍ Ｓｈｕｆｆｌｅ）法にて１次統計量を保存するＲＳサロゲートデータを作成すると共に、ＦＴ（Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）法にて２次統計量を保存するＦＴサロゲートデータを作成する処理を行う。
【００２４】
なお、非線形統計量推定部５０は、サロゲートデータの非線形統計量として、サロゲートデータ作成部６０にて作成されたＦＴサロゲートデータとＲＳサロゲートデータに対しても、ＧＰ法によってフラクタル次元の１つである相関次元Ｄｃを推定し、また、これらのサロゲートデータからリャプノフスペクトル及びリャプノフ指数を推定する処理も行う。
【００２５】
データ検定部７０は、切り出された時系列データ群から推定された非線形統計量と各サロゲートデータから推定された非線形統計量とをモンテカルロ有意性検定により比較して、当該の時系列データ群がカオス的な決定論に従うか否かを判断する処理を行う。
【００２６】
異常疑検出部８０は、データ検定部７０にてカオス的な決定論に従うと判断された時系列データ群を、後述する記憶部１１０から読み出された過去の時系列データ群の非線形統計量（フラクタル次元、リャプノフスペクトルなど）や可視化データ（リカレンスプロット、相空間中のアトラクタ軌道など）と比較して、その差を求めることにより、ロケットＲの異常疑のある動作を検出する処理を行う。
【００２７】
データ送信部９０は、データ切出部２０、データ前処理部３０、データ可視化処理部４０、非線形統計量推定部５０、サロゲートデータ作成部６０、データ検定部７０、及び異常疑検出部８０にて得られたデータを記憶部１１０へと送信したり、また、必要に応じて、これらのデータを観測基地Ｂの表示装置Ｐに送信したりする。
【００２８】
なお、データ受信部１０、データ切出部２０、データ前処理部３０、データ可視化処理部４０、非線形統計量推定部５０、サロゲートデータ作成部６０、データ検定部７０、及び異常疑検出部８０、データ送信部９０は、本実施の形態に係るテレメトリデータ解析装置１の主要部であるテレメトリデータ解析部１Ａを構成している。
【００２９】
また、入力部１００は、テレメトリデータ解析部１Ａにおける上記の処理を行う上で必要となるパラメータ（詳細は後述）の入力をユーザーに促す。
【００３０】
記憶部１１０は、テレメトリデータ解析部１Ａにおける上記の処理にて得られた各種データの格納を行うサブメモリと、上記の処理の各ステップ（後述）を実行するためのコンピュータに読み取り可能なプログラムを格納するメインメモリとから構成される。記憶手段１１０は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓｓ Ｍｅｍｏｒｙ）やＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）などから構成される。さらに、記憶部１１０のサブメモリとメインメモリとを別体として構成してもよい。
【００３１】
制御部１２０は、記憶部１１０から読み出したプログラムに従って、テレメトリデータ解析部１Ａの各動作を制御するＣＰＵ(Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ)を備える。特に、制御部１２０は、飛翔体Ｒの構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させるようにデータ前処理部３０を制御する点に特徴を有する。
【００３２】
本実施の形態では、テレメトリデータ解析装置１を、テレメトリデータ解析部１Ａと、入力部１００、記憶部１１０、及び制御部１２０とを一体化した構成としたが、記憶部１１０を独立した記憶装置としてテレメトリデータ解析部１Ａ、入力部１００、制御部１２０から切り離した構成としてもよい。また、テレメトリデータ解析部１Ａからデータ受送信のインターフェースであるデータ受信部１０及びデータ送信部９０を取り除いた構成をもってテレメトリデータ解析部１Ｂとしてもよい。
【００３３】
ここで、コンピュータとは、構造化された入力を所定の規則に従って処理し、処理した結果を構造化して出力する装置のことを指し、例えば、汎用コンピュータ、スーパーコンピュータ、メインフレーム、ワークステーション、マイクロコンピュータ、サーバ等が含まれる。また、通信ネットワーク（例えば、イントラネット、ローカルエリアネットワーク（LAN）、ワイドエリアネットワーク（WAN）、及びこれらの組み合わせから成る通信ネットワーク）を介して接続された２つ以上のコンピュータから成る構成（例えば、分散コンピュータシステム）であってもよい。
【００３４】
［データ解析方法］
次に、テレメトリデータ解析装置１において行われるテレメトリデータ解析の全体的な流れについて説明する。図４は、図３に示したテレメトリデータ解析装置において行われるテレメトリデータ解析の全体的な流れを示すフローチャートである。
【００３５】
データ受信部１０が観測基地Ｂの信号復調変換装置Ｄ（Ｄ１〜ＤＮ）からテレメトリデータ群を受け取ると、制御部１２０は、記憶部１１０に格納されたプログラムに従い、テレメトリデータ解析部１Ａの各部が以下に示す処理を行うように各部を制御する。
【００３６】
ステップＳ１０において、データ切出部２０は、データ受信部１０にて受信したテレメトリデータ群を読み込み、ステップＳ２０において、読み込んだテレメトリデータ群からデータ解析を実施する時間領域の時系列データ群を切り出す処理を行う。その際、制御部１２０は、ユーザーに対して、入力部１００からこの時間領域を入力するように促す。
【００３７】
ステップＳ３０において、制御部１２０は、データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群か否かを判断し、一次元の時系列データ群であると判断した場合には、データ前処理部３０に対して、誤り近傍点法によって当該時系列データ群を埋め込み可能とする相空間の次元（埋め込み次元）ｍを算出する処理を行わせ、切り出された時系列データ群が多次元であった場合には、その次元ｎを当該時系列データ群の埋め込み次元として設定する処理を行わせる。ここで、切り出された時系列データ群が何次元の時系列データ群であるかは飛翔体Ｒの構造の変化に依存する。従って、制御部１２０は、多段式ロケットにおいて各段が切り離されるなど飛翔体Ｒの構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる（又は、切替える）ようにデータ前処理部３０を制御する。
【００３８】
図５は、この埋め込み次元の設定処理の詳細を示すフローチャートである。データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群ｘ（ｔ）であった場合、ステップＳ４１にて、制御部１２０は、ユーザーに入力部１００から遅れ時間τの入力を促し、データ前処理部３０に対して、入力された遅れ時間τに対してｍ次元相空間Ｈを構成させる。このとき、切り出された時系列データ群の軌跡上の点は、構成されたｍ次元相空間Ｈにおいて、
【数１】
【００３９】
と表される。また、制御部１２０は、ｍ次元相空間Ｈにおける時系列データ群の軌跡ｙｍ（ｔ）を別途、観測基地Ｂの表示装置Ｐに表示させる。
【００４０】
このとき、適当な遅れ時間τを設定しないと、ｍ次元相空間Ｈにおける時系列データ群の軌跡がアトラクタに漸近するのかどうかの判断がつかないことになる。具体的には、遅れ時間τが小さすぎると相関が大きくなりすぎて、時系列データ群の軌跡がｍ次元相空間内の各象限の対角線上につぶれた形となる一方、遅れ時間τが大きすぎると相関が小さくなりすぎて、時系列データ群の軌跡はｍ次元相空間内においてランダムな軌跡を描くことになる。それゆえ、制御部１２０は、ユーザーに対して、表示装置Ｐに表示されたｍ次元相空間中の時系列データ群の軌跡ｙｍ（ｔ）の画像をチェックさせて、入力部１００から適当な遅れ時間τを入力するように促す。
【００４１】
次に、ステップＳ４２において、データ前処理部３０は、誤り近傍点法により最も低い埋め込み次元を求めて、それを当該相空間Ｈの次元ｍとして設定する。まず、データ前処理部３０は、式（１）で表されるｍ次元相空間Ｈにおける切り出された時系列データ群の軌跡上の点ｙｍ（ｔ）の最近傍点
【数２】
【００４２】
をヒープソートによって算出させる。ここで、ｎ（ｔ，ｍ）は最近傍点のインデックスである。そして、データ前処理部３０は、ｍ次元相空間Ｈでの２点間処理とｍ＋１次元相空間での２点間距離との比
【数３】
【００４３】
を計算する。ここで、分子、分母共に距離はユークリッドノルムで規定している。そして、データ前処理部３０は、全データ点数Ｎ（ｔ＝０，１，・・・，Ｎ）に対するａの平均値
【数４】
【００４４】
を算出し、相空間の次元をｍ次元からｍ＋１次元に変化させたときの変化率の指標
【数５】
【００４５】
を算出する。
【００４６】
次に、ステップＳ４３において、データ前処理部３０は、次元ｍを増加させて式（５）の値を計算させ、その値が収束するか否かを判定させる。切り出された時系列データ群が何らかの決定論的な法則に従っている場合には、式（５）の値は収束するので、データ前処理部３０は、閾値内まで式（５）の値が変化しなくなったときのｍの値を埋め込み可能な最低次元として設定する。この閾値は、制御部１２０が、ユーザーに対して、入力部１００から設定するように促す。
【００４７】
そして、ステップＳ４４において、データ前処理部３０は、設定された埋め込み次元ｍの相空間Ｈに切り出された時系列データ群を埋め込む処理を行う。
【００４８】
また、ステップＳ３０において、制御部１２０は、データ切出部２０にて切り出された時系列データ群が多次元の時系列データ群であると判断した場合には、データ前処理部３０に対して、その次元ｎを当該時系列データ群の埋め込み次元として設定させる。
【００４９】
また、ステップＳ４３において、収束解が見つからなかった場合には、ステップＳ４５において、制御部１２０は、ユーザーに入力部１００から適当な埋め込み次元ｍを設定するように促した後に、ステップＳ４４の処理を行わせる。
【００５０】
次に、ステップＳ５０において、データ可視化処理部４０は、切り出された時系列データ群を表示装置Ｐに表示させることを可能とするための可視化処理を行わせる。
【００５１】
図６は、この可視化処理の詳細を示すフローチャートである。まず、ステップＳ５１において、データ可視化処理部４０は、データ前処理部３０にて算出された埋め込み次元（ｍ又はｎ：以降では、仮にｍとする）を有する高次元相空間Ｈにおける切り出された時系列データ群の軌跡上のの２点ｙｍ（ｉ），ｙｍ（ｊ）間の距離
【数６】
【００５２】
をｉ，ｊの全ての組（ｉ，ｊ）に対して算出する。ここで、２点間距離はユークリッドノルムにて規定している。また、ｉ及びｊは共に、切り出された時系列データ群における時間のインデックスである。
【００５３】
ここで、制御部１２０は、ユーザーに対して、入力部１００から、切り出された時系列データ群を二値化してリカレンスプロットにて視覚化するための閾値Ｔｈ１を入力するように促す。次に、データ可視化処理部４０は、式（６）で算出されたＤｉｊを入力された閾値Ｔｈ１に基づき二値化する。その際、Ｄｉｊが閾値以下の点（ｉ，ｊ）をＴｒｕｅ，閾値以上の点（ｉ，ｊ）をＦａｌｓｅとする。また、制御部１２０は、データ可視化処理部４０にて二値化された時系列データ群をリカレンスプロットとして表示装置Ｐに表示させるために、式（６）で算出されたＤｉｊ（ｉ，ｊ＝０，１，・・・，Ｎ−１）をデータ送信部９０を介して表示装置Ｐへと送信することもできる。
【００５４】
次に、ステップＳ５２において、データ可視化処理部４０は、ｍ次元相空間Ｈにおける切り出された時系列データ群の軌跡上の２点ｙｍ（ｉ），ｙｍ（ｊ）における方向ベクトル
【数７】
【００５５】
間のばらつき具合を評価するために、
【数８】
【００５６】
を算出する。
【００５７】
ここで、制御部１２０は、ユーザーに対して、入力部１００から、切り出された時系列データ群を二値化して同方向リカレンスプロットにて視覚化するための閾値Ｔｈ２及び式（７），（８）におけるパラメータ（微小な時間差）Δを入力するように促す。次に、データ可視化処理部４０は、式（８）で算出されたＤＶｉｊを入力された閾値Ｔｈ２に基づき二値化する。その際、ＤＶｉｊが閾値以下の点（ｉ，ｊ）をＴｒｕｅ，閾値以上の点（ｉ，ｊ）をＦａｌｓｅとする。また、制御部１２０は、データ可視化処理部４０にて二値化された時系列データ群を同方向リカレンスプロットとして表示装置Ｐに表示させるために、式（８）で算出されたＤＶｉｊ（ｉ，ｊ＝０，１，・・・，Ｎ−１）をデータ送信部９０を介して表示装置Ｐへと送信することもできる。
【００５８】
次に、ステップＳ５３において、データ可視化処理部４０は、式（６）で算出されたＤｉｊ＝（Ｔｒｕｅ、又はＦａｌｓｅ）と式（８）で算出されたＤＶｉｊ＝（Ｔｒｕｅ、又はＦａｌｓｅ）との論理積
【数９】
【００５９】
を計算する。このＩＤＮＰも、Ｔｒｕｅ又はＦａｌｓｅの２値を取る。
【００６０】
次に、ステップＳ６０において、データ可視化処理部４０は、
【数１０】
【００６１】
を算出し、算出されたＲの値が１に近い場合には、切り出された時系列データ群は決定論的な法則に従うものと判断し、ステップＳ７０の処理に移り、そうでない場合には、切り出された時系列データ群は非決定論的であると判断して、本データ解析を終了する。ここで、制御部１２０は、ユーザーに対して、入力部１００から、Ｒの値が１に近いかどうかの判断基準となる閾値を入力するように促す。
【００６２】
次に、ステップＳ７０において、制御部１２０は、本データ解析がリアルタイム解析であるか否かを判断し、リアルタイム解析であると判断した場合には、ステップＳ８０において、制御部１２０は、非線形統計量推定部５０に対して、切り出された時系列データ群の非線形統計量を推定させる。
【００６３】
図７は、この非線形統計量の算出処理の詳細を示すフローチャートである。まず、ステップＳ８１において、非線形統計量推定部５０は、ＧＰ法によって切り出された時系列データ群からフラクタル次元の１つである相関次元Ｄｃを推定する。具体的には、非線形量推定部５０は、
【数１１】
【００６４】
を算出する。ここで、Ｉはヘビサイド関数
【数１２】
【００６５】
である。Ｃｍ（ｒ）は相関積分と呼ばれる量であり、定性的には、ｍ次元相空間Ｈにおける切り取られた時系列データ群の軌跡上の点ｙｍ（ｉ）を中心とした半径ｒのｍ次元超球内に入る点ｙｍ（ｊ）の数をカウントして、その全データ点数Ｎに対する割合を計算するということを軌跡上の全点ｙｍ（ｉ）に対して行い、その平均を取ったものである。
【００６６】
次に、非線形統計量推定部５０は、半径ｒの適当な範囲に対して、
【数１３】
【００６７】
を満たす次元ｄ（ｍ）を埋め込み次元ｍごとに算出する。具体的には、式（１３）を対数グラフにプロットした際にその傾斜が一定となるｌｏｇｒの領域におけるグラフの傾きを算出する。ここで、傾斜が一定であるとは、１ディケード（対数周波数軸で１０倍になる区間）以上に亘って傾斜の変動がある閾値以内である場合をいう。そこで、この際に、制御部１２０は、ユーザーに対して、入力部１００から、この閾値を入力するように促す。そして、非線形統計量推定部５０は、埋め込み次元ｍの変化に伴う収束先の傾斜をＤｃとして算出し、このＤｃが、条件
【数１４】
【００６８】
を満たすか否かを判定する。この条件を満たす場合、非線形統計量推定部５０は、Ｄｃを相関次元として推定する。一方、この条件を満たさない場合には、制御部１２０は、ユーザーに表示装置Ｐを介して、データ点数を増やし、本データ解析を再度実行するように通知する。
【００６９】
次に、ステップＳ８２において、非線形統計量推定部５０は、切り出された時系列データ群からリャプノフスペクトル及びリャプノフ次元を推定する処理を行う。
【００７０】
まず、非線形統計量推定部５０は、ｍ次元相空間Ｈにおける切り出された時系列データ群の軌道上にＮ個あるデータ点数の各点ｙｍ（ｔ）（ｔ＝０，１，・・・，Ｎ−１）を中心とするＭ個の近傍点ｙｍ（ｋｉ）（ｉ＝１，２，・・・，Ｍ）を含むｍ次元超球を考え、次の微小変位ベクトル
【数１５】
【００７１】
が
【数１６】
【００７２】
を満たすｍ×ｍ行列Ｊｔを算出する処理を行う。ここで、ｓは単位ステップ時間を表す。そのため、まず、制御部１２０は、ユーザーに対して、入力部１００から、超球内に含まれる軌道上の近傍点ｙｍ（ｋｉ）の個数Ｍを入力するように促す。なお、近傍点ｙｍ（ｋｉ）（ｉ＝１，２，・・・，Ｍ）について、制御部１２０は、ステップＳ６１にて算出されたＤｉｊをヒープソートアルゴリズムにて昇順にソーティングして、インデックスを呼び出す処理を行う。
【００７３】
次に、非線形統計量推定部５０は、時刻ｔ＝０に対して、
【数１７】
【００７４】
を算出する。ここで、εｉｋ及びεｉｌはそれぞれ、式（１５）の微小変位ベクトルεｉ（ｉ＝１，２，・・・，Ｍ）の第ｋ及び第ｌ成分である。また、ｗｋｌ及びｃｋｌはそれぞれ、
【数１８】
【００７５】
を満たすｍ×ｍ行列の成分であり、Ｗ及びＣはそれぞれ、分散行列及び共分散行列である。ちなみに、式（１８）は、
【数１９】
【００７６】
を最小にするための極小条件
【数２０】
【００７７】
から得られる関係式である。
【００７８】
そして、非線形統計量推定部５０は、算出された式（１７）を用いて式（１８）の右辺をＬＵ分解法により解き、軌道上の点ｙｍ（ｔ）（ｔ＝０）におけるＪｔを近似的なヤコビ行列として算出する。このとき、Ｍ≧ｍで縮退がなければ、ヤコビ行列Ｊｔを一意に算出することができる。そして、非線形統計量推定部５０は、この処理を全部のデータ点列に対して行う（すなわち、Ｎ回繰り返す）ことによって、Ｎ個のヤコビ行列Ｊｔ（ｔ＝０，１，・・・，Ｎ−１）を算出する。
【００７９】
次に、非線形統計量推定部５０は、予め用意されたｍ次元の正規直交基底
【数２１】
【００８０】
に対して上記処理にて算出された時刻ｔ＝０でのヤコビ行列
【数２２】
【００８１】
を乗じて、
【数２３】
【００８２】
へと写像する。ここで、
【数２４】
【００８３】
である。
【００８４】
そして、非線形統計量推定部５０は、式（２３）のｚ０ｉの各列ベクトルをグラム−シュミットの直交化法にて再度直交化して、
【数２５】
【００８５】
を算出する。このとき、ｉ番目の直交ベクトルのノルム‖ｖ１ｉ‖はヤコビ行列Ｊのｉ番目の固有値に等しく、その指数部分はリャプノフスペクトルを近似している。
【００８６】
次いで、非線形統計量推定部５０は、Ｖ１を再度正規化してＷ１を算出して、上記の処理をＮ回繰り返すことにより、Ｖ１，Ｖ２，・・・，ＶＮを算出して、Ｎ個の列ベクトルｖｊｉ（ｊ＝１，２，・・・，Ｎ）のノルムを平均して、リャプノフスペクトル
【数２６】
【００８７】
を推定する。
【００８８】
続いて、非線形統計量推定部５０は、式（２６）にて推定したリャプノフスペクトルλｉの和
【数２７】
【００８９】
を満たす最大の整数Ｐに対して、次の定義
【数２８】
【００９０】
に従い、リャプノフ次元ＤＬを算出させる。
【００９１】
次に、ステップＳ９０にて、制御部１２０は、非線形統計量推定部５０にて算出されたリャプノフスペクトルλｉの内の少なくとも最大リャプノフ指数λ１が正であれば、切り出しされた時系列データ群がカオス的な法則に従っていると判断して、ステップＳ１００にて、以上の処理過程において得られたデータを記憶部１１０に記憶させる。一方、最大リャプノフ指数λ１が負である場合には、切り出した時系列データ群はカオス的ではない決定論的な法則に従っているとして本データ解析を終了する。
【００９２】
ここで、ステップＳ７０の処理に戻って説明を続ける。ステップＳ７０において、制御部１２０は、本データ解析が、リアルタイム解析ではなく飛行後解析であると判断した場合には、ステップＳ１１０において、サロゲートデータ作成部６０に対して、切り出された時系列データ群からＲＳ（Ｒａｎｄｏｍ Ｓｈｕｆｆｌｅ）法にてＲＳサロゲートデータを作成させると共に、ＦＴ（Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）法にてＦＴサロゲートデータを作成させる。
【００９３】
図８は、このサロゲートデータの作成処理の詳細を示すフローチャートである。まず、ステップＳ１１１において、制御部１２０は、データ前処理部３０にて埋め込まれた時系列データ群が元々は一次元の時系列データ群であったか否かを判断し、一次元の時系列データ群であると判断した場合には、制御部１２０は、ステップＳ１１２において、サロゲートデータ作成部６０に対して、元の一次元の時系列データ群におけるデータ点列の並びをランダムに変化させたＲＳサロゲートデータを所定の組数作成させる。このとき作成されたＲＳサロゲートデータにおいては、平均値、分散、中央値などの１次統計量が保存されている。
【００９４】
次に、ステップＳ１１３において、サロゲートデータ作成部６０は、元の一次元の時系列データ群に高速フーリエ変換を施すことによって、当該一次元の時系列データ群を実空間から周波数空間へと変換する。そして、サロゲートデータ作成部６０は、高速フーリエ変換された周波数空間における各波の振幅を保存したままで位相だけランダム化した後に、逆高速フーリエ変換して周波数空間から実空間へと逆変換したＦＴサロゲートデータを所定の組数作成する。このとき作成されたＦＴサロゲートデータにおいては、パワースペクトル密度などの２次統計量が保存されている。
【００９５】
一方、ステップＳ１１１において、制御部１２０が、データ前処理部３０にて埋め込まれた時系列データ群が元々は多次元の時系列データ群であると判断した場合には、制御部１２０は、ステップＳ１１４において、サロゲートデータ作成部６０に対して、それぞれの（一次元の）時系列データ群に対して、ステップＳ１１２及びＳ１１３において、ＲＳ及びＴＦサロゲートデータを所定の組数作成させる。
【００９６】
なお、ステップＳ１１２において、制御部１２０は、ユーザーに対して、入力部１００から、ステップＳ１１２において作成するＲＳサロゲートデータの組数Ｓ１とステップＳ１１３において作成するＦＴサロゲートデータの組数Ｓ２を入力するように促す。以降では、サロゲートデータの総組数をＳ（＝Ｓ１＋Ｓ２）とする。
【００９７】
次に、ステップＳ１２０において、非線形統計量推定部５０は、切り出された時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータ（Ｓ１組のＲＳサロゲートデータとＳ２組のＦＴサロゲートデータ）のそれぞれに対して、非線形統計量として、相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）、及びリャプノフ次元ＤＬを推定する。
【００９８】
図９は、この非線形統計量の推定処理の詳細を示すフローチャートである。まず、ステップＳ１２１において、制御部１２０は、切り出された元の時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータのそれぞれが元々一次元の時系列データであったと判断した場合には、ステップＳ１２２において、データ前処理部３０にて算出された埋め込み次元ｍの相空間Ｈへのデータの埋め込みを行う。
【００９９】
次に、ステップＳ１２３において、非線形統計量推定部５０は、ｍ次元相空間Ｈに埋め込まれた元の時系列データ群と各サロゲートデータに対して、ステップＳ８１にて行った処理と同様の処理を施すことによって、フラクタル次元の１つである相関次元Ｄｃを推定する。
【０１００】
次に、ステップＳ１２４において、非線形統計量推定部５０は、ｍ次元相空間Ｈに埋め込まれた元の時系列データ群と各サロゲートデータに対して、ステップＳ８２にて行った処理と同様の処理を施すことによって、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）（のノルムと最大リャプノフ指数）、及びリャプノフ次元ＤＬを推定する。
【０１０１】
一方、ステップＳ１２１において、制御部１２０は、元の時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータのそれぞれが元々多次元の時系列データであったと判断した場合には、ステップＳ１２５において、その次元そのものを埋め込み次元とした相空間Ｈへのデータの埋め込みを行い、ステップＳ１２３へと処理を勧める。
【０１０２】
次に、ステップＳ１３０において、データ検定部７０は、非線形統計量推定部５０によって元の時系列データ群から推定された非線形統計量Ｑ０と各サロゲートデータから推定された非線形統計量Ｑｌ（ｌ＝１，２，・・・，Ｓ）とをモンテカルロ有意性検定により比較する。その際、制御部１２０は、ユーザーに入力部１００から、モンテカルロ有意性検定における有意水準α（例えば、α＝０．０５）を入力するように促す。ここで、Ｑ０及びＱｌは、相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）（のノルムと最大リャプノフ指数）、及びリャプノフ次元ＤＬである。
【０１０３】
次に、ステップＳ１４０において、制御部１２０は、Ｑ０がＱｌの最も大きい方から（Ｂｓ＋１）・α／２番目、或いは、最も小さい方から（Ｂｓ＋１）・α／２番目の間に挟まれているか否かを判定し、Ｑ０がこの間に挟まれていない場合には、元の時系列データ群はカオス的な法則に従っていると判断して、ステップＳ１５０へと処理を進める。
【０１０４】
一方、ステップＳ１４０において、Ｑ０がこの間に挟まれていると判断された場合には、元の時系列データ群はカオス的な法則に従ってはいないとして、本データ処理を終了する。
【０１０５】
次に、ステップＳ１５０において、異常疑検出部８０は、データ検定部７０にてカオス的な決定論に従うと判断された元の時系列データ群の非線形統計量Ｑ０（相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）、及びリャプノフ次元ＤＬ）や可視化データ（リカレンスプロット、同方向性リカレンスプロット、同方向近傍プロット、相空間中のアトラクタ軌道）を、記憶部１１０から読み出された過去の時系列データ群の非線形統計量、及びリャプノフ次元ＤＬ）や可視化データと比較して、その差を検出する。そして、制御部１２０は、異常疑検出部８０にて検出された差を表示装置Ｐなどによりユーザーに提示する。これによって、ユーザーは、解析対象である時系列データ群を送信した構造可変な飛翔体Ｒのカオス力学系に基づく異常疑のある動作を検出することができる。
【０１０６】
なお、本発明は、上述した実施の形態に限定されるものではなく、その要旨を逸脱しない範囲でその他の構成にても具現化することができる。
【０１０７】
例えば、その他の実施の形態として、３層ニューラルネットワークなどの学習機構を備えたテレメトリデータ解析部１Ｃを適用することもできる。テレメトリデータ解析部１Ｃでは、この学習機構を用いて、例えばバックプロパゲーション学習則に基づき、切り出された時系列データ群からカオス力学系のモデルを直接構築することによって、時系列データ群からカオス力学系を同定してしまう。そして、構築されたカオス力学系のモデルを用いて、解析対象となる時系列データ群からリャプノフスペクトルなどの非線形統計量を推定し、それを記憶部１１０に記憶された過去のデータと比較することによって、異常疑検出部８０にて、カオス力学系のモデルに基づく飛翔体の異常疑のある動作を検出することができる。
【符号の説明】
【０１０８】
１ テレメトリデータ解析装置
１Ａ テレメトリデータ解析部
１０ データ受信部
２０ データ切出部
３０ データ前処理部
４０ データ可視化処理部
５０ 非線形統計量推定部
６０ サロゲートデータ作成部
７０ データ検定部
８０ 異常疑検出部
９０ データ送信部
１００ 入力部
１１０ 記憶部
１２０ 制御部
【技術分野】
【０００１】
本発明は、ロケットなどの構造可変な飛翔体から送信されるテレメトリデータを解析して、飛翔体の異常疑のある動作を検出するテレメトリデータ解析装置及びテレメトリデータ解析プログラムに関する。
【背景技術】
【０００２】
飛行中のロケット、軌道上の人工衛星、無人宇宙船や有人宇宙船の状態は、これらの飛翔体から送信されるテレメトリデータを通して、地上の観測基地において常時観測されている。テレメトリデータとは、飛翔体の状態を表す各種の掲載機器から送信される時系列データのことであり、主なものとしては、各種搭載機器の電圧や電流、動作のオンオフ状態、飛翔体の各部の温度や圧力、送受信装置の送受信電力、飛翔体の外部環境に関する情報などのデータがある。これらのデータは各々、飛翔体に搭載されたテレメータ送信機のエンコーダにて収集され、一定のフォーマットにて地上の観測基地に送信される。図１０に示すように、観測基地Ｂにおいては、飛翔体（ロケットＲ）から送信された各種テレメトリデータＳ１〜ＳＮをアンテナ等の受信装置Ａ（Ａ１〜ＡＮ）にて受信した後、信号号復調変換装置Ｄ（Ｄ１〜ＤＮ）にて復調し、Ａ／Ｄ変換や各種の物理量変換を行う。そして、テレメトリデータ解析装置Ｚにおいて、変換された時系列データを統計解析や周波数解析したり、データの可視化を行い表示装置Ｐに表示させたりして、その結果を過去のデータ、設計値、又は予測値などと比較する。このような時系列データ解析により、飛翔体が計画通りに飛行できたことを確認したり、飛翔体の異常な動作を検知したり、また、飛翔体に何らかの事故が発生した場合にはその原因を調査して搭載機器の改良に役立てたりする（特許文献１）。
【０００３】
このような従来のテレメトリデータ解析（時系列データ解析）では、解析対象となるテレメトリデータを過去のデータ、設計値、又は予測値と比較したときに差異があまり大きくない場合には、その差異は通常ランダムノイズとして認識される。
【０００４】
しかし、近年数理物理学の分野において、一見するとランダムに見える時系列データがカオス力学系における決定論的な法則に従っている場合があることが明らかになってきた。カオス力学系の特徴は、そのダイナミクスがある種の非線形な因果関係に基づいているため、その時系列データがランダムノイズに似た複雑な振る舞いを示す点にある。ある時系列データがカオス力学系に基づく時系列データであるか否かを判定するためには、カオス力学系の特徴量であるフラクタル次元、リャプノフスペクトルなどの非線形統計量を算出したり、カオス力学系に特有なアトラクタ軌道（ストレンジアトラクタ）を相空間にて可視化したりするといったデータ解析を行う必要がある。このような時系列データ解析を特にカオス時系列データ解析といい、現状では主に生体内の各種生理信号、プラント等の状態信号など、各種工学システムの計測データに実用されている（非特許文献１）。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】米国特許公開２００９０１０５８９２号公報
【非特許文献】
【０００６】
【非特許文献１】合原 一幸、池口 徹、山田 泰司、小室 元政著、カオス時系列解析の基礎と応用、産業図書、２０００。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
従って、ロケット等の飛翔体のテレメトリデータに対する従来のデータ解析では、テレメトリデータが一見するとランダムに見えるような時系列データであった場合には、その中に潜む非線形性に基づく決定論的な要因、つまり飛翔体の動作における何らかのカオス力学系、の存在を見逃してしまう可能性があった。このことは、従来の飛翔体のテレメトリデータ解析では、カオス力学系に基づく飛翔体の異常疑のある動作を見逃してしまうことを意味する。
【０００８】
本発明は、上記の実情を鑑みて為されたものであり、ロケットなどの構造可変な飛翔体から送信されるテレメトリデータからカオス力学系に基づく飛翔体の異常疑な動作を検知することを目的とする。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、請求項１に記載の発明は、構造可変な飛翔体のテレメトリデータ解析装置であって、構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出部と、切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理部と、算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理部と、切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成部と、リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定部と、カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出部と、前記データ前処理部に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御部と、を備えたことを特徴とする。
【００１０】
また、請求項２に記載の発明は、コンピュータを、構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出手段と、切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理手段と、算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理手段と、切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成手段と、リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定手段と、飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定手段と、カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出手段と、前記データ前処理手段に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御手段と、して機能させることを特徴とする構造可変な飛翔体のテレメトリデータ解析プログラムであることを特徴とする。
【発明の効果】
【００１１】
本発明によれば、構造可変な飛翔体からの一見するとランダムのように見えるテレメトリデータの中からカオス力学系に基づく決定論的な要因を系統だって見つけ出すことが可能となるので、従来見過ごされていた可能性のある飛翔体の異常疑のある動作を検出することができるようになる。さらに、この装置を飛行後解析に適用した場合には、既に発生してしまった飛翔体の故障の原因を効率良く特定して、将来発生する可能性のある故障を予測することが可能となるので、次号機以降の飛翔体の設計の改善を行うことができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の一実施の形態に係るテレメトリデータ解析装置において解析対象となる多段式ロケットからのテレメトリデータを示す図である。
【図２】図１に示した多段式ロケットから観測基地に送られてくるテレメトリデータの特徴を示す図である。
【図３】本発明の一実施の形態に係るテレメトリデータ解析装置のブロック図である。
【図４】図３に示したテレメトリデータ解析装置において行われるテレメトリデータ解析の全体的な流れを示すフローチャートである。
【図５】図４に示したデータの前処理の詳細を示すフローチャートである。
【図６】図４に示したデータの可視化処理の詳細を示すフローチャートである。
【図７】図４に示した非線形統計量の推定処理の詳細を示すフローチャートである。
【図８】図４に示したサロゲートデータの作成処理の詳細を示すフローチャートである。
【図９】図４に示した非線形統計量の推定処理の詳細を示すフローチャートである。
【図１０】ロケットからのテレメトリデータを受信する観測基地の主要構成要素を示すブロック図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照しながら、本発明の一実施の形態に係る構造可変な飛翔体のテレメトリデータ解析装置及びテレメトリデータ解析プログラムについて説明する。図１は、本発明の一実施の形態に係るテレメトリデータ解析装置１（図３）において解析対象となる多段式ロケットＲからのテレメトリデータを示す図である。
【００１４】
［解析対象となるテレメトリデータ］
本発明において解析対象となるのは、図１に示すように、多段式ロケットのような構造可変な飛翔体Ｒから送信されるテレメトリデータである。多段式ロケットの一例として二段式ロケットについて説明する。本ロケットは、各段Ｒ１，Ｒ２に液体酸素と液体水素によるロケットエンジンを備えた２段式の液体燃料ロケットであり、１段目Ｒ１には推力増強のための大型固体ロケットブースタが２本装備されており、２段目Ｒ２の先端にペイロードフェアリングで覆われた人工衛星Ｒ３を装備している。また、さらに推力を増強するために、１段目に固体補助ロケットが複数本装備される場合もある。
【００１５】
これらのパーツＲ１〜Ｒ３にはそれぞれ、飛行中のロケットＲ（の搭載機器）の状態をチェックするためにテレメータ送信機が装備されており、これらのテレメータ送信機は、パーツＲ１〜Ｒ３に装備された各種の搭載機器の状態を表すテレメトリデータ群ＳＧ１〜ＳＧ３をそれぞれ地上の観測基地Ｂへとリアルタイムに送信する。テレメトリデータ群ＳＧ１〜ＳＧ３はそれぞれ、Ｎ１個のテレメトリデータＳ１１〜Ｓ１Ｎ１，Ｎ２個のテレメトリデータＳ２１〜Ｓ２Ｎ２，Ｎ３個のテレメトリデータＳ３１〜Ｓ３Ｎ３を含む。これらテレメトリデータ群ＳＧ１〜ＳＧ３中のテレメトリデータの個数は、一般にはパーツごと異なっている。なお、搭載機器には非常に様々なものがあるが、通信系、計測系、制御系、タイマ・点火系、電源・管制系、配線系などに関する搭載機器はロケットＲの可動に必須である。
【００１６】
パーツＲ１，Ｒ２は、ロケットＲの飛行中に下から順に作動し、ロケットＲが必要な速度に達するための推力を出し終わったら順次本体から切り離されていき、最終的には、人工衛星Ｒ３のみが地球の周回軌道上に投入される。このようなロケットＲからのパーツＲ１，Ｒ２の切り離しに応じて、観測基地Ｂに送信されてくるテレメトリデータは、図２に示すように時間的に変化する。図２は、図１に示した多段式ロケットＲから観測基地Ｂに送られてくるテレメトリデータの特徴を示す図である。時刻ｔ０にてロケットＲが打ち上げられ、時刻ｔ１にてパーツＲ１が品体から切り離されると、時刻ｔ１以降、パーツＲ１から観測基地Ｂに送信されていたテレメトリデータ群ＳＧ１は送信されなくなる。次いで、時刻ｔ２にてパーツＲ２が本体から切り離されると、時刻ｔ２以降、パーツＲ２から観測基地Ｂに送信されていたテレメトリデータ群ＳＧ２は送信されなくなる。そして、最終的には、それ以降（現時刻ｔ３）、周回軌道に投入された人工衛星Ｒ０からのテレメトリデータ群ＳＧ３のみが観測基地Ｂに送信され続けることになる。つまり、多段式ロケットＲのような構造可変な飛翔体から送信されるテレメトリデータは、観測基地Ｂから送信されるコマンドに応じて当該飛翔体の構造が変化する度に、観測基地Ｂにて受信されるテレメトリデータ数が変化することを特徴としている。
【００１７】
［解析装置の構成］
次に、本発明の一実施の形態に係るテレメトリデータ解析装置１の構成について説明する。図３は、本発明の一実施の形態に係るテレメトリデータ解析装置１のブロック図である。テレメトリデータ解析装置１は、データ受信部１０と、データ切出部２０と、データ前処理部３０と、データ可視化処理部４０と、非線形統計量推定部５０と、サロゲートデータ作成部６０と、データ検定部７０と、異常疑検出部８０と、データ送信部９０と、入力部１００と、記憶部１１０と、制御部１２０とを備えている。
【００１８】
データ受信部１０は、観測基地Ｂ（アンテナＡ１〜ＡＮ）にて受信され、信号復調変換装置Ｄ（Ｄ１〜ＤＮ）にて復調されＡ／Ｄ変換等が行われたテレメトリデータ群を受信する。
【００１９】
データ切出部２０は、データ受信部１０にて受信したテレメトリデータ群からテレメトリデータ解析を実施する時間領域の時系列データ群を切り出す処理を行う。
【００２０】
データ前処理部３０は、データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群であった場合には、誤り近傍点法によって当該時系列データ群が描く軌道を埋め込み可能とする相空間Ｈの次元（埋め込み次元）ｍを算出する処理を行い、切り出された時系列データ群が多次元であった場合には、その次元ｎを当該時系列データ群の埋め込み次元として設定する処理を行う。
【００２１】
データ可視化処理部４０は、データ前処理部３０にて算出された埋め込み次元を有する高次元相空間Ｈ中の２点間の距離を規定する所定のノルムに基づき時系列データ群を二値化して視覚化するリカレンスプロットと、高次元相空間Ｈ中の２点における方向ベクトル間のばらつき具合を評価するノルムに基づき時系列データ群を二値可して視覚化する同方向性リカレンスプロットと、これら２つのプロットの結果に基づき同方向的近傍プロットを行い、同方向的近傍プロットの結果に基づき時系列データ群が決定論的な法則に従っているか否かを判断する処理を行う。
【００２２】
非線形統計量推定部５０は、切り出された時系列データ群の非線形統計量として、ＧＰ（Ｇｒａｓｓｂｅｒｇｅｒ−Ｐｒｏｃａｃｃｉａ）法によって当該の時系列データ群からフラクタル次元の１つである相関次元Ｄｃを推定し、また、この時系列データ群からリャプノフスペクトル及びリャプノフ指数を推定する処理を行う。
【００２３】
サロゲートデータ作成部６０は、本データ解析がリアルタイム解析ではなく飛行後解析である場合に、切り出された時系列データ群からＲＳ（Ｒａｎｄｏｍ Ｓｈｕｆｆｌｅ）法にて１次統計量を保存するＲＳサロゲートデータを作成すると共に、ＦＴ（Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）法にて２次統計量を保存するＦＴサロゲートデータを作成する処理を行う。
【００２４】
なお、非線形統計量推定部５０は、サロゲートデータの非線形統計量として、サロゲートデータ作成部６０にて作成されたＦＴサロゲートデータとＲＳサロゲートデータに対しても、ＧＰ法によってフラクタル次元の１つである相関次元Ｄｃを推定し、また、これらのサロゲートデータからリャプノフスペクトル及びリャプノフ指数を推定する処理も行う。
【００２５】
データ検定部７０は、切り出された時系列データ群から推定された非線形統計量と各サロゲートデータから推定された非線形統計量とをモンテカルロ有意性検定により比較して、当該の時系列データ群がカオス的な決定論に従うか否かを判断する処理を行う。
【００２６】
異常疑検出部８０は、データ検定部７０にてカオス的な決定論に従うと判断された時系列データ群を、後述する記憶部１１０から読み出された過去の時系列データ群の非線形統計量（フラクタル次元、リャプノフスペクトルなど）や可視化データ（リカレンスプロット、相空間中のアトラクタ軌道など）と比較して、その差を求めることにより、ロケットＲの異常疑のある動作を検出する処理を行う。
【００２７】
データ送信部９０は、データ切出部２０、データ前処理部３０、データ可視化処理部４０、非線形統計量推定部５０、サロゲートデータ作成部６０、データ検定部７０、及び異常疑検出部８０にて得られたデータを記憶部１１０へと送信したり、また、必要に応じて、これらのデータを観測基地Ｂの表示装置Ｐに送信したりする。
【００２８】
なお、データ受信部１０、データ切出部２０、データ前処理部３０、データ可視化処理部４０、非線形統計量推定部５０、サロゲートデータ作成部６０、データ検定部７０、及び異常疑検出部８０、データ送信部９０は、本実施の形態に係るテレメトリデータ解析装置１の主要部であるテレメトリデータ解析部１Ａを構成している。
【００２９】
また、入力部１００は、テレメトリデータ解析部１Ａにおける上記の処理を行う上で必要となるパラメータ（詳細は後述）の入力をユーザーに促す。
【００３０】
記憶部１１０は、テレメトリデータ解析部１Ａにおける上記の処理にて得られた各種データの格納を行うサブメモリと、上記の処理の各ステップ（後述）を実行するためのコンピュータに読み取り可能なプログラムを格納するメインメモリとから構成される。記憶手段１１０は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓｓ Ｍｅｍｏｒｙ）やＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）などから構成される。さらに、記憶部１１０のサブメモリとメインメモリとを別体として構成してもよい。
【００３１】
制御部１２０は、記憶部１１０から読み出したプログラムに従って、テレメトリデータ解析部１Ａの各動作を制御するＣＰＵ(Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ)を備える。特に、制御部１２０は、飛翔体Ｒの構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させるようにデータ前処理部３０を制御する点に特徴を有する。
【００３２】
本実施の形態では、テレメトリデータ解析装置１を、テレメトリデータ解析部１Ａと、入力部１００、記憶部１１０、及び制御部１２０とを一体化した構成としたが、記憶部１１０を独立した記憶装置としてテレメトリデータ解析部１Ａ、入力部１００、制御部１２０から切り離した構成としてもよい。また、テレメトリデータ解析部１Ａからデータ受送信のインターフェースであるデータ受信部１０及びデータ送信部９０を取り除いた構成をもってテレメトリデータ解析部１Ｂとしてもよい。
【００３３】
ここで、コンピュータとは、構造化された入力を所定の規則に従って処理し、処理した結果を構造化して出力する装置のことを指し、例えば、汎用コンピュータ、スーパーコンピュータ、メインフレーム、ワークステーション、マイクロコンピュータ、サーバ等が含まれる。また、通信ネットワーク（例えば、イントラネット、ローカルエリアネットワーク（LAN）、ワイドエリアネットワーク（WAN）、及びこれらの組み合わせから成る通信ネットワーク）を介して接続された２つ以上のコンピュータから成る構成（例えば、分散コンピュータシステム）であってもよい。
【００３４】
［データ解析方法］
次に、テレメトリデータ解析装置１において行われるテレメトリデータ解析の全体的な流れについて説明する。図４は、図３に示したテレメトリデータ解析装置において行われるテレメトリデータ解析の全体的な流れを示すフローチャートである。
【００３５】
データ受信部１０が観測基地Ｂの信号復調変換装置Ｄ（Ｄ１〜ＤＮ）からテレメトリデータ群を受け取ると、制御部１２０は、記憶部１１０に格納されたプログラムに従い、テレメトリデータ解析部１Ａの各部が以下に示す処理を行うように各部を制御する。
【００３６】
ステップＳ１０において、データ切出部２０は、データ受信部１０にて受信したテレメトリデータ群を読み込み、ステップＳ２０において、読み込んだテレメトリデータ群からデータ解析を実施する時間領域の時系列データ群を切り出す処理を行う。その際、制御部１２０は、ユーザーに対して、入力部１００からこの時間領域を入力するように促す。
【００３７】
ステップＳ３０において、制御部１２０は、データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群か否かを判断し、一次元の時系列データ群であると判断した場合には、データ前処理部３０に対して、誤り近傍点法によって当該時系列データ群を埋め込み可能とする相空間の次元（埋め込み次元）ｍを算出する処理を行わせ、切り出された時系列データ群が多次元であった場合には、その次元ｎを当該時系列データ群の埋め込み次元として設定する処理を行わせる。ここで、切り出された時系列データ群が何次元の時系列データ群であるかは飛翔体Ｒの構造の変化に依存する。従って、制御部１２０は、多段式ロケットにおいて各段が切り離されるなど飛翔体Ｒの構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる（又は、切替える）ようにデータ前処理部３０を制御する。
【００３８】
図５は、この埋め込み次元の設定処理の詳細を示すフローチャートである。データ切出部２０にて切り出された時系列データ群が一次元の時系列データ群ｘ（ｔ）であった場合、ステップＳ４１にて、制御部１２０は、ユーザーに入力部１００から遅れ時間τの入力を促し、データ前処理部３０に対して、入力された遅れ時間τに対してｍ次元相空間Ｈを構成させる。このとき、切り出された時系列データ群の軌跡上の点は、構成されたｍ次元相空間Ｈにおいて、
【数１】
【００３９】
と表される。また、制御部１２０は、ｍ次元相空間Ｈにおける時系列データ群の軌跡ｙｍ（ｔ）を別途、観測基地Ｂの表示装置Ｐに表示させる。
【００４０】
このとき、適当な遅れ時間τを設定しないと、ｍ次元相空間Ｈにおける時系列データ群の軌跡がアトラクタに漸近するのかどうかの判断がつかないことになる。具体的には、遅れ時間τが小さすぎると相関が大きくなりすぎて、時系列データ群の軌跡がｍ次元相空間内の各象限の対角線上につぶれた形となる一方、遅れ時間τが大きすぎると相関が小さくなりすぎて、時系列データ群の軌跡はｍ次元相空間内においてランダムな軌跡を描くことになる。それゆえ、制御部１２０は、ユーザーに対して、表示装置Ｐに表示されたｍ次元相空間中の時系列データ群の軌跡ｙｍ（ｔ）の画像をチェックさせて、入力部１００から適当な遅れ時間τを入力するように促す。
【００４１】
次に、ステップＳ４２において、データ前処理部３０は、誤り近傍点法により最も低い埋め込み次元を求めて、それを当該相空間Ｈの次元ｍとして設定する。まず、データ前処理部３０は、式（１）で表されるｍ次元相空間Ｈにおける切り出された時系列データ群の軌跡上の点ｙｍ（ｔ）の最近傍点
【数２】
【００４２】
をヒープソートによって算出させる。ここで、ｎ（ｔ，ｍ）は最近傍点のインデックスである。そして、データ前処理部３０は、ｍ次元相空間Ｈでの２点間処理とｍ＋１次元相空間での２点間距離との比
【数３】
【００４３】
を計算する。ここで、分子、分母共に距離はユークリッドノルムで規定している。そして、データ前処理部３０は、全データ点数Ｎ（ｔ＝０，１，・・・，Ｎ）に対するａの平均値
【数４】
【００４４】
を算出し、相空間の次元をｍ次元からｍ＋１次元に変化させたときの変化率の指標
【数５】
【００４５】
を算出する。
【００４６】
次に、ステップＳ４３において、データ前処理部３０は、次元ｍを増加させて式（５）の値を計算させ、その値が収束するか否かを判定させる。切り出された時系列データ群が何らかの決定論的な法則に従っている場合には、式（５）の値は収束するので、データ前処理部３０は、閾値内まで式（５）の値が変化しなくなったときのｍの値を埋め込み可能な最低次元として設定する。この閾値は、制御部１２０が、ユーザーに対して、入力部１００から設定するように促す。
【００４７】
そして、ステップＳ４４において、データ前処理部３０は、設定された埋め込み次元ｍの相空間Ｈに切り出された時系列データ群を埋め込む処理を行う。
【００４８】
また、ステップＳ３０において、制御部１２０は、データ切出部２０にて切り出された時系列データ群が多次元の時系列データ群であると判断した場合には、データ前処理部３０に対して、その次元ｎを当該時系列データ群の埋め込み次元として設定させる。
【００４９】
また、ステップＳ４３において、収束解が見つからなかった場合には、ステップＳ４５において、制御部１２０は、ユーザーに入力部１００から適当な埋め込み次元ｍを設定するように促した後に、ステップＳ４４の処理を行わせる。
【００５０】
次に、ステップＳ５０において、データ可視化処理部４０は、切り出された時系列データ群を表示装置Ｐに表示させることを可能とするための可視化処理を行わせる。
【００５１】
図６は、この可視化処理の詳細を示すフローチャートである。まず、ステップＳ５１において、データ可視化処理部４０は、データ前処理部３０にて算出された埋め込み次元（ｍ又はｎ：以降では、仮にｍとする）を有する高次元相空間Ｈにおける切り出された時系列データ群の軌跡上のの２点ｙｍ（ｉ），ｙｍ（ｊ）間の距離
【数６】
【００５２】
をｉ，ｊの全ての組（ｉ，ｊ）に対して算出する。ここで、２点間距離はユークリッドノルムにて規定している。また、ｉ及びｊは共に、切り出された時系列データ群における時間のインデックスである。
【００５３】
ここで、制御部１２０は、ユーザーに対して、入力部１００から、切り出された時系列データ群を二値化してリカレンスプロットにて視覚化するための閾値Ｔｈ１を入力するように促す。次に、データ可視化処理部４０は、式（６）で算出されたＤｉｊを入力された閾値Ｔｈ１に基づき二値化する。その際、Ｄｉｊが閾値以下の点（ｉ，ｊ）をＴｒｕｅ，閾値以上の点（ｉ，ｊ）をＦａｌｓｅとする。また、制御部１２０は、データ可視化処理部４０にて二値化された時系列データ群をリカレンスプロットとして表示装置Ｐに表示させるために、式（６）で算出されたＤｉｊ（ｉ，ｊ＝０，１，・・・，Ｎ−１）をデータ送信部９０を介して表示装置Ｐへと送信することもできる。
【００５４】
次に、ステップＳ５２において、データ可視化処理部４０は、ｍ次元相空間Ｈにおける切り出された時系列データ群の軌跡上の２点ｙｍ（ｉ），ｙｍ（ｊ）における方向ベクトル
【数７】
【００５５】
間のばらつき具合を評価するために、
【数８】
【００５６】
を算出する。
【００５７】
ここで、制御部１２０は、ユーザーに対して、入力部１００から、切り出された時系列データ群を二値化して同方向リカレンスプロットにて視覚化するための閾値Ｔｈ２及び式（７），（８）におけるパラメータ（微小な時間差）Δを入力するように促す。次に、データ可視化処理部４０は、式（８）で算出されたＤＶｉｊを入力された閾値Ｔｈ２に基づき二値化する。その際、ＤＶｉｊが閾値以下の点（ｉ，ｊ）をＴｒｕｅ，閾値以上の点（ｉ，ｊ）をＦａｌｓｅとする。また、制御部１２０は、データ可視化処理部４０にて二値化された時系列データ群を同方向リカレンスプロットとして表示装置Ｐに表示させるために、式（８）で算出されたＤＶｉｊ（ｉ，ｊ＝０，１，・・・，Ｎ−１）をデータ送信部９０を介して表示装置Ｐへと送信することもできる。
【００５８】
次に、ステップＳ５３において、データ可視化処理部４０は、式（６）で算出されたＤｉｊ＝（Ｔｒｕｅ、又はＦａｌｓｅ）と式（８）で算出されたＤＶｉｊ＝（Ｔｒｕｅ、又はＦａｌｓｅ）との論理積
【数９】
【００５９】
を計算する。このＩＤＮＰも、Ｔｒｕｅ又はＦａｌｓｅの２値を取る。
【００６０】
次に、ステップＳ６０において、データ可視化処理部４０は、
【数１０】
【００６１】
を算出し、算出されたＲの値が１に近い場合には、切り出された時系列データ群は決定論的な法則に従うものと判断し、ステップＳ７０の処理に移り、そうでない場合には、切り出された時系列データ群は非決定論的であると判断して、本データ解析を終了する。ここで、制御部１２０は、ユーザーに対して、入力部１００から、Ｒの値が１に近いかどうかの判断基準となる閾値を入力するように促す。
【００６２】
次に、ステップＳ７０において、制御部１２０は、本データ解析がリアルタイム解析であるか否かを判断し、リアルタイム解析であると判断した場合には、ステップＳ８０において、制御部１２０は、非線形統計量推定部５０に対して、切り出された時系列データ群の非線形統計量を推定させる。
【００６３】
図７は、この非線形統計量の算出処理の詳細を示すフローチャートである。まず、ステップＳ８１において、非線形統計量推定部５０は、ＧＰ法によって切り出された時系列データ群からフラクタル次元の１つである相関次元Ｄｃを推定する。具体的には、非線形量推定部５０は、
【数１１】
【００６４】
を算出する。ここで、Ｉはヘビサイド関数
【数１２】
【００６５】
である。Ｃｍ（ｒ）は相関積分と呼ばれる量であり、定性的には、ｍ次元相空間Ｈにおける切り取られた時系列データ群の軌跡上の点ｙｍ（ｉ）を中心とした半径ｒのｍ次元超球内に入る点ｙｍ（ｊ）の数をカウントして、その全データ点数Ｎに対する割合を計算するということを軌跡上の全点ｙｍ（ｉ）に対して行い、その平均を取ったものである。
【００６６】
次に、非線形統計量推定部５０は、半径ｒの適当な範囲に対して、
【数１３】
【００６７】
を満たす次元ｄ（ｍ）を埋め込み次元ｍごとに算出する。具体的には、式（１３）を対数グラフにプロットした際にその傾斜が一定となるｌｏｇｒの領域におけるグラフの傾きを算出する。ここで、傾斜が一定であるとは、１ディケード（対数周波数軸で１０倍になる区間）以上に亘って傾斜の変動がある閾値以内である場合をいう。そこで、この際に、制御部１２０は、ユーザーに対して、入力部１００から、この閾値を入力するように促す。そして、非線形統計量推定部５０は、埋め込み次元ｍの変化に伴う収束先の傾斜をＤｃとして算出し、このＤｃが、条件
【数１４】
【００６８】
を満たすか否かを判定する。この条件を満たす場合、非線形統計量推定部５０は、Ｄｃを相関次元として推定する。一方、この条件を満たさない場合には、制御部１２０は、ユーザーに表示装置Ｐを介して、データ点数を増やし、本データ解析を再度実行するように通知する。
【００６９】
次に、ステップＳ８２において、非線形統計量推定部５０は、切り出された時系列データ群からリャプノフスペクトル及びリャプノフ次元を推定する処理を行う。
【００７０】
まず、非線形統計量推定部５０は、ｍ次元相空間Ｈにおける切り出された時系列データ群の軌道上にＮ個あるデータ点数の各点ｙｍ（ｔ）（ｔ＝０，１，・・・，Ｎ−１）を中心とするＭ個の近傍点ｙｍ（ｋｉ）（ｉ＝１，２，・・・，Ｍ）を含むｍ次元超球を考え、次の微小変位ベクトル
【数１５】
【００７１】
が
【数１６】
【００７２】
を満たすｍ×ｍ行列Ｊｔを算出する処理を行う。ここで、ｓは単位ステップ時間を表す。そのため、まず、制御部１２０は、ユーザーに対して、入力部１００から、超球内に含まれる軌道上の近傍点ｙｍ（ｋｉ）の個数Ｍを入力するように促す。なお、近傍点ｙｍ（ｋｉ）（ｉ＝１，２，・・・，Ｍ）について、制御部１２０は、ステップＳ６１にて算出されたＤｉｊをヒープソートアルゴリズムにて昇順にソーティングして、インデックスを呼び出す処理を行う。
【００７３】
次に、非線形統計量推定部５０は、時刻ｔ＝０に対して、
【数１７】
【００７４】
を算出する。ここで、εｉｋ及びεｉｌはそれぞれ、式（１５）の微小変位ベクトルεｉ（ｉ＝１，２，・・・，Ｍ）の第ｋ及び第ｌ成分である。また、ｗｋｌ及びｃｋｌはそれぞれ、
【数１８】
【００７５】
を満たすｍ×ｍ行列の成分であり、Ｗ及びＣはそれぞれ、分散行列及び共分散行列である。ちなみに、式（１８）は、
【数１９】
【００７６】
を最小にするための極小条件
【数２０】
【００７７】
から得られる関係式である。
【００７８】
そして、非線形統計量推定部５０は、算出された式（１７）を用いて式（１８）の右辺をＬＵ分解法により解き、軌道上の点ｙｍ（ｔ）（ｔ＝０）におけるＪｔを近似的なヤコビ行列として算出する。このとき、Ｍ≧ｍで縮退がなければ、ヤコビ行列Ｊｔを一意に算出することができる。そして、非線形統計量推定部５０は、この処理を全部のデータ点列に対して行う（すなわち、Ｎ回繰り返す）ことによって、Ｎ個のヤコビ行列Ｊｔ（ｔ＝０，１，・・・，Ｎ−１）を算出する。
【００７９】
次に、非線形統計量推定部５０は、予め用意されたｍ次元の正規直交基底
【数２１】
【００８０】
に対して上記処理にて算出された時刻ｔ＝０でのヤコビ行列
【数２２】
【００８１】
を乗じて、
【数２３】
【００８２】
へと写像する。ここで、
【数２４】
【００８３】
である。
【００８４】
そして、非線形統計量推定部５０は、式（２３）のｚ０ｉの各列ベクトルをグラム−シュミットの直交化法にて再度直交化して、
【数２５】
【００８５】
を算出する。このとき、ｉ番目の直交ベクトルのノルム‖ｖ１ｉ‖はヤコビ行列Ｊのｉ番目の固有値に等しく、その指数部分はリャプノフスペクトルを近似している。
【００８６】
次いで、非線形統計量推定部５０は、Ｖ１を再度正規化してＷ１を算出して、上記の処理をＮ回繰り返すことにより、Ｖ１，Ｖ２，・・・，ＶＮを算出して、Ｎ個の列ベクトルｖｊｉ（ｊ＝１，２，・・・，Ｎ）のノルムを平均して、リャプノフスペクトル
【数２６】
【００８７】
を推定する。
【００８８】
続いて、非線形統計量推定部５０は、式（２６）にて推定したリャプノフスペクトルλｉの和
【数２７】
【００８９】
を満たす最大の整数Ｐに対して、次の定義
【数２８】
【００９０】
に従い、リャプノフ次元ＤＬを算出させる。
【００９１】
次に、ステップＳ９０にて、制御部１２０は、非線形統計量推定部５０にて算出されたリャプノフスペクトルλｉの内の少なくとも最大リャプノフ指数λ１が正であれば、切り出しされた時系列データ群がカオス的な法則に従っていると判断して、ステップＳ１００にて、以上の処理過程において得られたデータを記憶部１１０に記憶させる。一方、最大リャプノフ指数λ１が負である場合には、切り出した時系列データ群はカオス的ではない決定論的な法則に従っているとして本データ解析を終了する。
【００９２】
ここで、ステップＳ７０の処理に戻って説明を続ける。ステップＳ７０において、制御部１２０は、本データ解析が、リアルタイム解析ではなく飛行後解析であると判断した場合には、ステップＳ１１０において、サロゲートデータ作成部６０に対して、切り出された時系列データ群からＲＳ（Ｒａｎｄｏｍ Ｓｈｕｆｆｌｅ）法にてＲＳサロゲートデータを作成させると共に、ＦＴ（Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）法にてＦＴサロゲートデータを作成させる。
【００９３】
図８は、このサロゲートデータの作成処理の詳細を示すフローチャートである。まず、ステップＳ１１１において、制御部１２０は、データ前処理部３０にて埋め込まれた時系列データ群が元々は一次元の時系列データ群であったか否かを判断し、一次元の時系列データ群であると判断した場合には、制御部１２０は、ステップＳ１１２において、サロゲートデータ作成部６０に対して、元の一次元の時系列データ群におけるデータ点列の並びをランダムに変化させたＲＳサロゲートデータを所定の組数作成させる。このとき作成されたＲＳサロゲートデータにおいては、平均値、分散、中央値などの１次統計量が保存されている。
【００９４】
次に、ステップＳ１１３において、サロゲートデータ作成部６０は、元の一次元の時系列データ群に高速フーリエ変換を施すことによって、当該一次元の時系列データ群を実空間から周波数空間へと変換する。そして、サロゲートデータ作成部６０は、高速フーリエ変換された周波数空間における各波の振幅を保存したままで位相だけランダム化した後に、逆高速フーリエ変換して周波数空間から実空間へと逆変換したＦＴサロゲートデータを所定の組数作成する。このとき作成されたＦＴサロゲートデータにおいては、パワースペクトル密度などの２次統計量が保存されている。
【００９５】
一方、ステップＳ１１１において、制御部１２０が、データ前処理部３０にて埋め込まれた時系列データ群が元々は多次元の時系列データ群であると判断した場合には、制御部１２０は、ステップＳ１１４において、サロゲートデータ作成部６０に対して、それぞれの（一次元の）時系列データ群に対して、ステップＳ１１２及びＳ１１３において、ＲＳ及びＴＦサロゲートデータを所定の組数作成させる。
【００９６】
なお、ステップＳ１１２において、制御部１２０は、ユーザーに対して、入力部１００から、ステップＳ１１２において作成するＲＳサロゲートデータの組数Ｓ１とステップＳ１１３において作成するＦＴサロゲートデータの組数Ｓ２を入力するように促す。以降では、サロゲートデータの総組数をＳ（＝Ｓ１＋Ｓ２）とする。
【００９７】
次に、ステップＳ１２０において、非線形統計量推定部５０は、切り出された時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータ（Ｓ１組のＲＳサロゲートデータとＳ２組のＦＴサロゲートデータ）のそれぞれに対して、非線形統計量として、相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）、及びリャプノフ次元ＤＬを推定する。
【００９８】
図９は、この非線形統計量の推定処理の詳細を示すフローチャートである。まず、ステップＳ１２１において、制御部１２０は、切り出された元の時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータのそれぞれが元々一次元の時系列データであったと判断した場合には、ステップＳ１２２において、データ前処理部３０にて算出された埋め込み次元ｍの相空間Ｈへのデータの埋め込みを行う。
【００９９】
次に、ステップＳ１２３において、非線形統計量推定部５０は、ｍ次元相空間Ｈに埋め込まれた元の時系列データ群と各サロゲートデータに対して、ステップＳ８１にて行った処理と同様の処理を施すことによって、フラクタル次元の１つである相関次元Ｄｃを推定する。
【０１００】
次に、ステップＳ１２４において、非線形統計量推定部５０は、ｍ次元相空間Ｈに埋め込まれた元の時系列データ群と各サロゲートデータに対して、ステップＳ８２にて行った処理と同様の処理を施すことによって、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）（のノルムと最大リャプノフ指数）、及びリャプノフ次元ＤＬを推定する。
【０１０１】
一方、ステップＳ１２１において、制御部１２０は、元の時系列データ群と、サロゲートデータ作成部６０にて作成されたＳ組のサロゲートデータのそれぞれが元々多次元の時系列データであったと判断した場合には、ステップＳ１２５において、その次元そのものを埋め込み次元とした相空間Ｈへのデータの埋め込みを行い、ステップＳ１２３へと処理を勧める。
【０１０２】
次に、ステップＳ１３０において、データ検定部７０は、非線形統計量推定部５０によって元の時系列データ群から推定された非線形統計量Ｑ０と各サロゲートデータから推定された非線形統計量Ｑｌ（ｌ＝１，２，・・・，Ｓ）とをモンテカルロ有意性検定により比較する。その際、制御部１２０は、ユーザーに入力部１００から、モンテカルロ有意性検定における有意水準α（例えば、α＝０．０５）を入力するように促す。ここで、Ｑ０及びＱｌは、相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）（のノルムと最大リャプノフ指数）、及びリャプノフ次元ＤＬである。
【０１０３】
次に、ステップＳ１４０において、制御部１２０は、Ｑ０がＱｌの最も大きい方から（Ｂｓ＋１）・α／２番目、或いは、最も小さい方から（Ｂｓ＋１）・α／２番目の間に挟まれているか否かを判定し、Ｑ０がこの間に挟まれていない場合には、元の時系列データ群はカオス的な法則に従っていると判断して、ステップＳ１５０へと処理を進める。
【０１０４】
一方、ステップＳ１４０において、Ｑ０がこの間に挟まれていると判断された場合には、元の時系列データ群はカオス的な法則に従ってはいないとして、本データ処理を終了する。
【０１０５】
次に、ステップＳ１５０において、異常疑検出部８０は、データ検定部７０にてカオス的な決定論に従うと判断された元の時系列データ群の非線形統計量Ｑ０（相関次元Ｄｃ、リャプノフスペクトルλｉ（ｉ＝１，２，・・・，Ｎ）、及びリャプノフ次元ＤＬ）や可視化データ（リカレンスプロット、同方向性リカレンスプロット、同方向近傍プロット、相空間中のアトラクタ軌道）を、記憶部１１０から読み出された過去の時系列データ群の非線形統計量、及びリャプノフ次元ＤＬ）や可視化データと比較して、その差を検出する。そして、制御部１２０は、異常疑検出部８０にて検出された差を表示装置Ｐなどによりユーザーに提示する。これによって、ユーザーは、解析対象である時系列データ群を送信した構造可変な飛翔体Ｒのカオス力学系に基づく異常疑のある動作を検出することができる。
【０１０６】
なお、本発明は、上述した実施の形態に限定されるものではなく、その要旨を逸脱しない範囲でその他の構成にても具現化することができる。
【０１０７】
例えば、その他の実施の形態として、３層ニューラルネットワークなどの学習機構を備えたテレメトリデータ解析部１Ｃを適用することもできる。テレメトリデータ解析部１Ｃでは、この学習機構を用いて、例えばバックプロパゲーション学習則に基づき、切り出された時系列データ群からカオス力学系のモデルを直接構築することによって、時系列データ群からカオス力学系を同定してしまう。そして、構築されたカオス力学系のモデルを用いて、解析対象となる時系列データ群からリャプノフスペクトルなどの非線形統計量を推定し、それを記憶部１１０に記憶された過去のデータと比較することによって、異常疑検出部８０にて、カオス力学系のモデルに基づく飛翔体の異常疑のある動作を検出することができる。
【符号の説明】
【０１０８】
１ テレメトリデータ解析装置
１Ａ テレメトリデータ解析部
１０ データ受信部
２０ データ切出部
３０ データ前処理部
４０ データ可視化処理部
５０ 非線形統計量推定部
６０ サロゲートデータ作成部
７０ データ検定部
８０ 異常疑検出部
９０ データ送信部
１００ 入力部
１１０ 記憶部
１２０ 制御部
【特許請求の範囲】
【請求項１】
構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出部と、
切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理部と、
算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理部と、
切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成部と、
リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、
飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定部と、
カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出部と、
前記データ前処理部に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御部と、
を備えたことを特徴とする構造可変な飛翔体のテレメトリデータ解析装置。
【請求項２】
コンピュータを、
構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出手段と、
切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理手段と、
算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理手段と、
切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成手段と、
リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定手段と、
飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定手段と、
カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出手段と、
前記データ前処理手段に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御手段と、
して機能させることを特徴とする構造可変な飛翔体のテレメトリデータ解析プログラム。
【請求項１】
構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出部と、
切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理部と、
算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理部と、
切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成部と、
リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定部と、
飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定部と、
カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出部と、
前記データ前処理部に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御部と、
を備えたことを特徴とする構造可変な飛翔体のテレメトリデータ解析装置。
【請求項２】
コンピュータを、
構造可変な飛翔体から送信されるテレメトリデータ群から解析対象となる時系列データ群を切り出すデータ切出手段と、
切り出された時系列データ群を埋め込み可能とする相空間の埋め込み次元を算出するデータ前処理手段と、
算出された次元を有する相空間に埋め込まれた時系列データ群から可視化情報を算出するデータ可視化処理手段と、
切り出された時系列データ群から１次統計量を保存する第１のサロゲートデータ群と２次統計量を保存する第２のサロゲートデータ群とを作成するサロゲートデータ作成手段と、
リアルタイム解析時には、可視化された時系列データ群から非線形統計量を推定し、飛行後解析時には、可視化された時系列データ群と第１及び第２のサロゲートデータ群とから非線形統計量を推定する非線形統計量推定手段と、
飛行後解析時に、第１及び第２のサロゲートデータ群から推定された非線形統計量を母集団として、可視化された時系列データ群から推定された非線形統計量の有意性検定を行い、切り出された時系列データ群がカオス力学系に従うか否かを判定するデータ検定手段と、
カオス力学系に従うと判定された時系列データ群の非線形統計量及び視覚化情報と過去の時系列データ群の非線形統計量及び視覚化情報との差を算出し、カオス力学系に基づく飛翔体の異常疑のある動作を検出する異常疑検出手段と、
前記データ前処理手段に対して、飛翔体の構造が変化する度に、切り出された時系列データ群の埋め込み次元を算出させる制御手段と、
して機能させることを特徴とする構造可変な飛翔体のテレメトリデータ解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１１−６８３２２（Ｐ２０１１−６８３２２Ａ）
【公開日】平成２３年４月７日（２０１１．４．７）
【国際特許分類】
【出願番号】特願２００９−２２３２７４（Ｐ２００９−２２３２７４）
【出願日】平成２１年９月２８日（２００９．９．２８）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２１年度独立行政法人新エネルギー・産業技術総合開発機構「次世代輸送系システム設計基盤技術開発 ミッション対応設計高度化技術」委託研究、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（０００００００９９）株式会社ＩＨＩ (5,014)
【公開日】平成２３年４月７日（２０１１．４．７）
【国際特許分類】
【出願日】平成２１年９月２８日（２００９．９．２８）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２１年度独立行政法人新エネルギー・産業技術総合開発機構「次世代輸送系システム設計基盤技術開発 ミッション対応設計高度化技術」委託研究、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（０００００００９９）株式会社ＩＨＩ (5,014)
[ Back to top ]