JP2005017032A - 目標運動解析方法及び装置 - Google Patents
目標運動解析方法及び装置 Download PDFInfo
- Publication number
- JP2005017032A JP2005017032A JP2003179895A JP2003179895A JP2005017032A JP 2005017032 A JP2005017032 A JP 2005017032A JP 2003179895 A JP2003179895 A JP 2003179895A JP 2003179895 A JP2003179895 A JP 2003179895A JP 2005017032 A JP2005017032 A JP 2005017032A
- Authority
- JP
- Japan
- Prior art keywords
- target
- motion
- observation
- observation ship
- analysis
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Landscapes
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
【課題】水中及び水上を運動する目標物の運動解析収束時間を最小化することが可能に観測船の運動パターンを決定する。
【解決手段】水中及び水上を運動する音源から放射される音波を音響センサ２０１により観測し、得られる到来音波の方位時系列データ及び周波数時系列データを、観測船の位置の時系列データと共に記憶器２０５に格納する。目標運動解析装置２１０は、最適運動演算器２０８を備え、目標及び観測船の運動パターンに関する予測計算を実施し、音響センサの受信可能範囲や観測船運動性能に関する制約条件を考慮してフィッシャー情報行列が最大化される観測船の運動パターンを計算してリコメンドし、観測船を操舵させる。これにより、解析収束時間を加速することができる。
【選択図】 図２
【解決手段】水中及び水上を運動する音源から放射される音波を音響センサ２０１により観測し、得られる到来音波の方位時系列データ及び周波数時系列データを、観測船の位置の時系列データと共に記憶器２０５に格納する。目標運動解析装置２１０は、最適運動演算器２０８を備え、目標及び観測船の運動パターンに関する予測計算を実施し、音響センサの受信可能範囲や観測船運動性能に関する制約条件を考慮してフィッシャー情報行列が最大化される観測船の運動パターンを計算してリコメンドし、観測船を操舵させる。これにより、解析収束時間を加速することができる。
【選択図】 図２
Description
【０００１】
【発明の属する技術分野】
本発明は、目標運動解析方法及び装置に係り、特に、水中及び水上を運動する目標が有する音源から放射される音波を音響センサにより観測し、そこから得られる到来音波の方位及び周波数を用いて目標の運動諸元を解析する目標運動解析方法及び装置に関する。
【０００２】
【従来の技術】
一般に、目標運動の解析は、観測船に設けられた音響センサを使用して、目標が放射する音波を受信することにより、到来音波の方位及び周波数を観測し、観測された方位及び周波数、観測時刻及び観測位置の時系列データに基づいて目標の状態ベクトル（距離、方位、速力、針路、放射音波の周波数）を解析するという方法により行われている。
【０００３】
前述した解析方法は、状態ベクトルから推定される各観測時刻の解析方位及び解析周波数と、観測方位及び周波数との残差二乗を観測誤差分散により正規化し、これらを時系列的に合計した評価関数が最小となるような状態ベクトルを解析する非線形最適化問題に帰着する。
【０００４】
前述のような非線形最適化問題を解析する方法としては、これまで、カルマンフィルタや最小二乗法によるものが知られており、カルマンフィルタを用いた目標運動解析については、非特許文献１等に記載されたものが、また、最小二乗法を用いた目標運動解析については、非特許文献２等に記載されたものが知られている。
【０００５】
【非特許文献１】
ＩＥＥＥ会報１９８５年 Ｖｏｌ ＡＣ−３０ Ｎｏ １０ 「Ａ Ｓｔｏｃｈａｓｔｉｃ Ａｎａｌｙｓｉｓ ｏｆ ａ Ｍｏｄｉｆｉｅｄ Ｇａｉｎ Ｅｘｔｅｎｄｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ ｗｉｔｈ Ａｐｐｌｉｃａｔｉｏｎ ｔｏ Ｅｓｔｉｍａｔｉｏｎ ｗｉｔｈ Ｂｅａｒｉｎｇ Ｏｎｌｙ Ｍｅａｓｕｒｅｍｅｎｔｓ」
【０００６】
【非特許文献２】
ＩＥＥＥ会報１９８４年 Ｖｏｌ ＡＣ−２９ Ｎｏ ９ 「Ｆｕｎｄａｍｅｎｔａｌ Ｐｒｏｐｅｒｔｉｅｓ ａｎｄ Ｐｅｒｆｏｒｍａｎｃｅ ｏｆ Ｃｏｎｖｅｎｔｉｏｎａｌ Ｂｅａｒｉｎｇ−ｏｎｌｙ Ｔａｒｇｅｔ Ｍｏｔｉｏｎ Ａｎａｌｙｓｉｓ」
【０００７】
【発明が解決しようとする課題】
前述した従来技術は、何れも、既に観測された到来音波の方位時系列データ及び周波数時系列データを用いて目標運動諸元を解析する技術であるが、目標が遠距離で時系列データの変化率が小さく、時系列データ誤差が大きい場合で、時系列データが十分蓄積されていない場合に目標運動諸元を一義的に決定する条件である可観測性が成立せず、解析が不能となるという問題点を有している。
【０００８】
また、前述した従来技術は、何れも、観測船が運動を実施しても適切な運動でなければ解析不能な状態から脱することができないという問題点を有している。
【０００９】
図９は目標が遠距離で時系列データの変化率が小さく、時系列データの誤差が大きい場合で、時系列データが十分蓄積されていない場合に解析が不能となるという問題点を説明する図、図１０は観測船の運動が適切でない場合に解析不能な状態から脱することができないという問題点について説明する図であり、これらの図を参照して、前述した問題点について説明する。
【００１０】
図９において、音響センサの受信可能範囲が９０２として示すような範囲であり、観測船９０１が停止状態にあるとし、さらに、周波数時系列データが得られないものとする。このような場合、方位時系列データの観測誤差９０３が、図示のように大きい場合、観測結果から得られる目標の運動ベクトルが誤差の範囲となってしまい、目標の運動ベクトル（真値）９０５と目標の運動ベクトル解析値（誤値）９０４、９０６との区別がつかないことになる。
【００１１】
また、観測船が運動を実施しても適切な運動でない場合、例えば、図１０に示すように、観測船１００１が目標と並行運動して時系列データを蓄積しても、目標運動ベクトル（真値）１００４と目標運動ベクトル解析値（誤値）１００３、１００５との区別がつかない。
【００１２】
さらに、従来から知られている目標運動解析装置は、パッシブ音響センサで受信した音波の到来方位及び周波数が海中雑音による誤差やセンサ特有の測定誤差を有していたため、解析精度を一定値以下に収束させるのに長時間を要するという問題点を有している。また、従来、方位変化率を最大化する観測船の運動パターンが解析収束時間を加速させると考えられていたが、音響センサの受信可能範囲や観測船の運動性能の制約から観測船をどのように操舵するのが最適かを指示する指針は存在していない。
【００１３】
本発明の目的は、前述した従来技術の問題点を解決し、目標物の運動解析収束時間を最小化することが可能に観測船の運動パターンを決定して、目標物の運動解析を行う目標運動解析方法及び装置を提供することにある。
【００１４】
【課題を解決するための手段】
本発明によれば前記目的は、水中及び水上を運動する目標が有する音源から放射される音波を観測船に設けた音響センサにより観測し、そこから得られる到来音波の方位時系列データ及び周波数時系列データを用いて目標の距離、方位、速力、針路を含む運動諸元を解析する目標運動解析方法において、目標及び観測船の運動パターンに関する予測計算を実施し、算出された目標及び観測船の予測運動パターンについて、ヤコビアン行列の可観測性を評価することにより、解析収束時間を最小化することができる観測船の運動パターンを指示することにより達成される。
【００１５】
前述した手段を備える本発明によれば、その結果としての観測船の運動パターンを観測船の操舵装置へ送出することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。
【００１６】
【発明の実施の形態】
以下、本発明による目標運動解析方法及び装置の実施形態を図面により詳細に説明する。
【００１７】
図１は本発明に適用する目標運動解析の原理を説明する図である。図１において、１０１は観測船、１０２は目標である。
【００１８】
観測船１０１は、自身も移動しながら音響センサを使用して、水上または水中を移動する船舶等の目標１０２が放射する音波を受信することにより、到来音波の方位Ｂ〜ｙ（なお、〜 は、本来直前の符号の上に付けて表記すべきものである。以下の説明で、＾ と記述するものもあるが、これも同様である）及び周波数ｆ〜 を観測する。そして、観測された方位及び周波数、そして観測時刻及び観測位置の時系列データに基づいて、目標１０２の状態ベクトルＸ（距離ｒ、方位Ｂｙ 、速力Ｍｔ 、針路Ｃｔ 、放射音波の周波数ｆ０）が解析される。図１０には、観測船１０１が移動しながら到来音波の方位Ｂ〜ｙ及び周波数ｆ〜 を観測し、それらの観測結果の時系列データから目標１０２の状態ベクトルＸ（距離ｒ、方位Ｂｙ 、速力Ｍｔ 、針路Ｃｔ 、放射音波の周波数ｆ）を得ることができることを示している。なお、状態ベクトルを示す各要素の時刻の順を各要素の符号の後に１、２、……、ｉ、ｌを付けて示している。
【００１９】
前述した解析は、状態ベクトルＸから推定される各観測時刻の解析方位Ｂｙ 及び解析周波数ｆと、観測方位Ｂ〜ｙ及び周波数ｆ〜 との残差二乗を観測誤差分散σＢ〜ｙ２ 及びσｆ〜 ２により正規化し、これらを時系列的に合計した評価関数Φ（Ｘ） が最小となるような状態ベクトルＸを解析する非線形最適化問題として行うことができ、図１１に示す式（１）、式（２）を解くことにより解析することができる。
【００２０】
図２は本発明の一実施形態による目標運動解析装置の構成例を示すブロック図である。図２において、２０１は音響センサ、２０２は航海センサ、２０３は操舵装置、２０４は位置演算器、２０５は観測データベース記憶器、２０６は状態ベクトル推定器、２０７は解析値表示器、２０８は最適運動演算器、２０９は最適運動表示器、２１０は目標運動解析装置、２１１はオペレータコンソールである。
【００２１】
図２に示す本発明の実施形態による目標運動解析装置２１０は、音響センサ２０１で受信した目標音源からの音波到来方位及び受信周波数と、自船のログ（船速）及びジャイロ（方位）を計測する航海センサ２０２からの計測値に基づき位置演算器２０４で算出された観測船位置とを観測データベース記録器２０５に、時刻を整合させた時系列の観測データベースとして蓄積する。状態ベクトル推定器２０６は、ここで蓄積した観測データベースに対して最小二乗法や適応型フィルタにより目標運動諸元に関する非線形解析を実施し、その解析過程を解析値表示器２０７に表示すると共に、現在時刻における目標運動ベクトル解析値を最適運動演算器２０８に送出する。最適運動演算器２０８は、状態ベクトル推定器２０６で解析された目標運動ベクトルと観測船の複数の運動パターンの候補に基づいて将来の予測計算を実施し、図１３に示す式（１６）を使用してフィッシャー情報行列値の増加状況を比較評価する。そして、最適運動演算器２０８は、式（１６）に示すフィッシャー情報行列値が最大となる観測船の運動パターンを最適運動表示器２０９に表示する。なお、図１３における式（１７）〜式（２０）は、式（１６）の各項を示す式である。
【００２２】
通常、観測船の操舵装置２０３は、オペレータコンソール２１１からの指示により作動するが、最適運動演算器２０８により算出された運動パターンを操舵装置２０３に送出することにより、状態ベクトル推定器２０６における可観測性が成立し、目標の運動解析が可能な状態を生成して、解析収束速度を加速させるような観測船の運動を実施させることが可能となる。
【００２３】
図３は観測船の最適運動演算の原理を説明する図である。図３において、観測船３０８の現在位置３０５とその時点における目標運動ベクトルの解析値３０３及び観測方位誤差３０９との関係から、目標の運動ベクトルが現在の解析値３０３を維持して運動することを前提に予測計算を行うことにより、予測運動ベクトル３０４を得ることができる。そして、観測船の運動パターンに関する複数の候補３０６について、それぞれの運動を実施した場合における図１３に示すフィッシャー情報行列値の式（１６）を算出して比較評価し、フィッシャー情報行列値の式（１６）が最大化される観測船の運動パターンを決定して操舵することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。また、目標の運動ベクトルの現在の解析値３０３と観測船の運動パターンに関する複数の候補３０６との関係から予測計算を行うことにより、予測運動ベクトル３０４を得る場合、音響センサの受信可能範囲３０７や誤差特性及び観測船３０８の運動性能を考慮することにより、より正確に観測船３０８の運動パターンを決定することができる。
【００２４】
なお、初回観測時点から目標の運動ベクトルの現在の解析値３０３が得られるまでは、音響センサの受信感度から目標距離を推定すると共に、図１１に示している式（３）を用いて、音波到来方位の変化方向に目標の運動ベクトル３０１を仮定し、前述した最適運動演算処理を実施する。
【００２５】
図４は最適運動演算器２０８での処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図４には２つのフローを示しているが、ステップ４０１〜４０３を持つフローは、もう一方のフローにおけるステップ４１１の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００２６】
（１）処理が開始されると、まず、図１３に示す式（１６）の各観測船の運動パターン番号毎のフィッシャー情報行列値を初期化（ＦＩＭ（０）＝０に初期化）し、観測船の運動パターン番号を１に初期化（現在の運動パターン）する（ステップ４０４〜４０７）。
【００２７】
（２）次に、観測船の運動パターン番号ｎと目標運動ベクトルの解析値とに基づいて将来の仮想位置の予測計算を実施し、音響センサの受信可能範囲３０７で受信可能な仮想観測データを図１３に示す式（１８）を用いて生成し、さらに、生成した仮想観測データを仮想観測データベースに蓄積してデータベースを更新する（ステップ４０８〜４１０）。
【００２８】
（３）ここで、蓄積された仮想観測データベースの情報に基づいて、図１３に示す式（１６）フィッシャー情報行列値の計算を実施する（ステップ４１１）。
【００２９】
（４）フィッシャー情報行列値の計算の処理において、まず、図１３に式（１７）として示す予測ヤコビアン行列及び図１３に式（２０）として示す予測重み行列を計算し、さらに、図１３に示す式（１６）のフィッシャー情報行列値を算出する（ステップ４０１〜４０３）。
【００３０】
（５）次に、観測船の運動パターン番号ｎに関するフィッシャー情報行列値ＦＩＭ（ｎ）が観測船の運動パターン番号ｎ−１に関するフィッシャー情報行列値ＦＩＭ（ｎ−１）よりも大きくなるか否かを判定し、大きくなる場合、最適運動パターン番号をｎに置換える（ステップ４１２、４１３）。
【００３１】
（６）ステップ４１２の判定で、観測船の運動パターン番号ｎに関するフィッシャー情報行列値が観測船の運動パターン番号ｎ−１に関するフィッシャー情報行列値よりも大きくならない場合、また、ステップ４１３の処理後、観測船の全ての運動パターン番号に関する比較評価が完了しているか否か、すなわち、運動パターンｎが最大値を越えたか否かを判定し、完了していれば、ここでの処理を終了し、完了していなければ、運動パターン番号をインクリメントして、ステップ４０８からの処理に戻って、処理を繰り返す（ステップ４１４、４１５）。
【００３２】
図５は Ｍａｒｑｕａｒｄｔ法を利用した最小二乗法による状態ベクトル推定器２０６の処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図５には２つのフローを示しているが、ステップ５０１〜５０４を持つフローは、もう一方のフローにおけるステップ５０６の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００３３】
（１）疑似線形解析の処理の実行が開始されると、まず、データベース内の全観測時刻のデータを使用して、初期値算出の処理が実行される（ステップ５０５、５０６）。
【００３４】
（２）初期値算出の処理において、まず、観測データベースに蓄積されている全観測時刻のデータを使用し、観測誤差が微少であることを仮定した疑似線形解析の処理を行って初期解析値を算出する。この解析値は、ステップ５０５での非線形解析処理により反復修正されて最終的な解析値として算出される（ステップ５０１）。
【００３５】
（３）前回時刻における非線形解析処理５０５による解析値が「解有り」として存在するか否かを判定し、存在した場合、これを今回の初期解析値とし、疑似線形解析処理をスキップする（ステップ５０２）。
【００３６】
（４）ステップ５０２の判定で、前回の処理での解析値が存在していなかった場合、目標運動諸元の状態ベクトルを図１１に示す式（４）のように定義し、係数行列Ｆ及び定数行列Ｇのそれぞれを、図１１に示す式（６）、式（７）により計算する線形解析準備の処理を実行する。その後、図１１に示す式（５）により初期解析値を算出する（ステップ５０３、５０４）。
【００３７】
（５）次に、非線形解析処理で、目標運動諸元の状態ベクトルを、図１２に示す式（８）のように定義し、前述した疑似線形解析処理で算出した初期解析値とこれに基づく目標音源の放射周波数を、図１２に示す式（９）を用いて算出して設定する。
【００３８】
（６）次に、Ｍａｒｑｕａｒｄｔ 法解析における反復修正を安定的行うことを目的とした縮小因子λの初期値を設定する縮小因子λ設定の処理を実行し、さらに、反復修正回数の初期化（Ｌｏｏｐ１＝０）を行う（ステップ５０７、５０８）。
【００３９】
（７）次に、図１１に示している式（１）により、各時刻における観測情報と解析情報との残差の二乗和Φを算出する残差二乗和算出処理を実行し、さらに、図１２に示す式（１１）を用いて、方位情報及び周波数情報に関する偏微分係数を算出するヤコビアン行列形成の処理を実行する（ステップ５０９、５１０）。
【００４０】
（８）ステップ５０９の処理で算出した残差二乗行列とステップ５１０の処理で算出したヤコビアン行列とを用い、図１２に示す式（１０）により状態ベクトル修正量Δｘの算出を行う修正量Δｘ算出処理を実行する。なお、図１２における式（１１）〜式（１５）は、前述の式（１０）の各項を示す式である（ステップ５１１）。
【００４１】
（９）反復修正回数Ｌｏｏｐ１ が最大値を越えたか否か、あるいは、前回の残差の二乗和Φが今回の残差の二乗和Φより大きいか否かを判定し、判定結果がＮＧであった場合、縮小因子λを２倍にしてステップ５０９からの処理に戻って処理を繰り返す（ステップ５１２、５１８）。
【００４２】
（１０）ステップ５１２の判定で、判定結果がＯＫであり、判定条件を満足していた場合、状態ベクトルｘを、ｘ＋Δｘに修正し、反復修正回数Ｌｏｏｐ１ が１より大きいか否かを判定し、大きかった場合、縮小因子λを半分にして次回の反復修正処理に備える（ステップ５１３〜５１５）。
【００４３】
（１１）ステップ５１４の判定で、反復修正回数Ｌｏｏｐ１ が１より大きくなかった場合、あるいは、ステップ５１５の処理後、状態ベクトルの修正量Δｘが極小値ｍｉｎ より小さくなったか否か、反復修正回数Ｌｏｏｐ２ が上限値に達した否かを判定する。判定の結果がＮＧであった場合、ステップ５０９からの処理に戻って処理を繰り返す（ステップ５１６）。
【００４４】
（１２）ステップ５１６の判定で、状態ベクトルの修正量が極小値ｍｉｎ より小さくなるか、反復修正回数Ｌｏｏｐ２ が上限値に達した場合、反復を打ち切り、図１２に示す式（８）で求められた状態ベクトルを目標運動諸元（距離ｒ，方位Ｂｙ ，速力Ｍｔ ，針路Ｃｔ ，放射周波数ｆ０ ）に変換して処理を終了する（ステップ５１７）。
【００４５】
前述した処理において、ステップ５１６の判定で、反復修正回数Ｌｏｏｐ２ が上限値に達して終了した場合「解なし」として扱う。
【００４６】
なお、誤差を伴って観測される情報から内部の状態モデルを推定する最小二乗法及びＭａｒｑｕａｒｄｔ 法の詳細については以下の文献に記載されている。
【００４７】
「最小二乗法による実験データ解析」中川徹 他 東京大学出版会
図６はカルマンフィルタを利用した適応型フィルタによる状態ベクトル推定器２０６の処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図６には２つのフローを示しているが、ステップ６０１〜６０３を持つフローは、もう一方のフローにおけるステップ６０５の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００４８】
ここでの処理は、目標運動諸元の状態ベクトルを、図１４に式（２１）として示すように定義し、観測データベースに蓄積されている各観測時刻毎のデータを１つ１つ使用して状態ベクトル及び共分散行列の初期値を逐次更新することにより解析値を収束させる処理である。
【００４９】
（１）逐次推定の処理の実行が開始されると、初期値設定の処理を実行する（ステップ６０４、６０５）。
【００５０】
（２）初期値設定の処理において、まず、前回の解析値が存在しない場合、図１４に示す式（２２）を用いて、状態ベクトルの初期値を設定する状態ベクトル初期化の処理を実行し、図１４に示す式（２３）を用いて、共分散行列の初期値を設定する共分散行列初期化の処理を実行する（ステップ６０１〜６０３）。
【００５１】
（３）次に、状態ベクトルの初期値を、観測データべースから取り出したデータの観測時刻まで、図１４に示す式（２４）により時間遷移させる状態ベクトル状態遷移処理を実行し、共分散行列状態遷移では、共分散行列の初期値を、観測データべースから取り出したデータの観測時刻まで、図１４に示す式（２５）により時間遷移させる共分散行列状態遷移の処理を実行する（ステップ６０６、６０７）。
【００５２】
（４）次に、図１４に示す式（２６）により、観測データベースから取り出した観測データに基づき状態ベクトルを修正更新する状態ベクトルフィルタリングの処理を実行し、図１４に示す式（２７）により、共分散行列を修正更新する共分散ベクトルフィルタリングの処理を実行する（ステップ６０８、６０９）。
【００５３】
（５）さらに、図１４に示す式（２８）により状態ベクトルを目標運動諸元に変換し、図１４に示す式（２９）により共分散行列を目標運動諸元解析値の解析誤差標準偏差に変換する目標運動諸元変換の処理を実行する（ステップ６１０）。
【００５４】
（６）最後に、観測データベースファイルがＥＯＦ（終了）となっているか否かを判定し、ＥＯＦとなっていた場合、ここでの処理を終了し、そうでない場合、ステップ６０６からの処理に戻り、観測データベースから取り出すデータがなくなるまで処理を繰り返す（ステップ６１１）。
【００５５】
なお、誤差を伴って観測される情報から内部の状態モデルを推定するカルマンフィルタの詳細については以下の文献に記載されている。
【００５６】
「応用カルマンフィルタ」片山徹 朝倉書店
前述した本発明の実施形態における各処理は、処理プログラムとして構成することができ、この処理プログラムは、ＨＤ、ＤＡＴ、ＦＤ、ＭＯ、ＤＶＤ−ＲＯＭ、ＣＤ−ＲＯＭ等の記録媒体に格納して提供することができる。
【００５７】
図７は観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に時系列方位データの変化率が最大化される方向へ針路を変更した場合における状態ベクトル推定器２０６の解析結果の状況を示す図である。
【００５８】
図７（ａ）には、観測船７０９の航跡と目標運動ベクトル（真値）７１１と解析目標位置７１０とを表示した平面図を示している。また、図７（ｂ）には、グラフ７０２〜７０６を示している。グラフ７０２には、横軸を時間軸とし、縦軸を距離とした目標距離（真値）７１２と解析距離７１３との時系列推移が表示されており、グラフ７０３には、横軸を時間軸とし、縦軸を速力とした目標速力（真値）７１４と解析速力７１５との時系列推移が表示されている。また、グラフ７０４には、横軸を時間軸とし、縦軸を針路とした目標針路（真値）７１６と解析針路７１７との時系列推移が表示されており、グラフ７０５には、横軸を時間軸とし、縦軸を方位とした観測船針路７０７と最適運動リコメンド７０８との時系列推移が表示されている。さらに、グラフ７０６には、横軸を時間軸としたフィッシャー情報行列値７１８の時系列推移が表示されている。
【００５９】
図８は観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に最適運動演算器の計算結果に従って針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【００６０】
図８（ａ）には、観測船８０９の航跡と目標運動ベクトル（真値）８１１と解析目標位置８１０とを表示した平面図を示している。また、図８（ｂ）には、グラフ８０２〜８０６を示している。グラフ８０２には、横軸を時間軸とし、縦軸を距離とした目標距離（真値）８１２と解析距離８１３との時系列推移が表示されており、グラフ８０３には、横軸を時間軸とし、縦軸を速力とした目標速力（真値）８１４と解析速力８１５との時系列推移が表示されている。また、グラフ８０４には、横軸を時間軸とし、縦軸を針路とした目標針路（真値）８１６と解析針路８１７との時系列推移が表示されており、グラフ８０５には、横軸を時間軸とし、縦軸を方位とした観測船針路８０７と最適運動リコメンド８０８との時系列推移が表示されている。さらに、グラフ８０６には、横軸を時間軸としたフィッシャー情報行列値８１８の時系列推移が表示されている。
【００６１】
図７と図８とに示す解析結果を比較してみると、観測船の運動パターンが図８（ａ）に示すように、フィッシャー情報行列を最大化する運動の場合、図７における方位変化率最大化の場合に比較してフッシャー情報行列値８１８が初期段階から増加傾向を示しており、目標距離の解析結果８１３及び目標速力の解析結果８１５の収束時間が加速されていることが判る。
【００６２】
前述した本発明の実施形態によれば、観測船の現在位置とその時点における目標運動ベクトルの解析値及び観測方位誤差との関係から、目標の運動ベクトルが現在の解析値を維持して運動することを前提に予測計算を行い、観測船の運動パターンに関する複数の候補について、それぞれの運動を実施した場合におけるフィッシャー情報行列値を算出して比較評価し、フィッシャー情報行列値を最大化することのできる観測船の運動パターンを決定することができる。そして、この運動パターンにより観測船を操舵することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。
【００６３】
また、前述した本発明の実施形態によれば、目標の運動ベクトルの現在の解析値と観測船の運動パターンに関する複数の候補との関係を用いて予測計算を行う場合に、音響センサの受信可能範囲や誤差特性及び観測船の運動性能を考慮することにより、より正確な観測船の運動パターンを決定することができる。
【００６４】
【発明の効果】
以上説明したように本発明によれば、目標物の運動解析収束時間を最小化することが可能に観測船の運動パターンを決定して、観測船を運動させることができる。
【図面の簡単な説明】
【図１】本発明に適用する目標運動解析の原理を説明する図である。
【図２】本発明の一実施形態による目標運動解析装置の構成例を示すブロック図である。
【図３】観測船の最適運動演算の原理を説明する図である。
【図４】最適運動演算器での処理動作を説明するフローチャートである。
【図５】Ｍａｒｑｕａｒｄｔ 法を利用した最小二乗法による状態ベクトル推定器の処理動作を説明するフローチャートである。
【図６】カルマンフィルタを利用した適応型フィルタによる状態ベクトル推定器の処理動作を説明するフローチャートである。
【図７】観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に時系列方位データの変化率が最大化される方向へ針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【図８】観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に最適運動演算器の計算結果に従って針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【図９】目標が遠距離で時系列データの変化率が小さく、時系列データの誤差が大きい場合で、時系列データが十分蓄積されていない場合に解析が不能となるという問題点を説明する図である。
【図１０】観測船の運動が適切でない場合に解析不能な状態から脱することができないという問題点について説明する図である。
【図１１】説明の中で使用される数式を示す図（その１）である。
【図１２】説明の中で使用される数式を示す図（その２）である。
【図１３】説明の中で使用される数式を示す図（その３）である。
【図１４】説明の中で使用される数式を示す図（その４）である。
【符号の説明】
１０１ 観測船
１０２ 目標
２０１ 音響センサ
２０２ 航海センサ
２０３ 操舵装置
２０４ 位置演算器
２０５ 観測データベース記憶器
２０６ 状態ベクトル推定器
２０７ 解析値表示器
２０８ 最適運動演算器
２０９ 最適運動表示器
２１０ 目標運動解析装置
２１１ オペレータコンソール
【発明の属する技術分野】
本発明は、目標運動解析方法及び装置に係り、特に、水中及び水上を運動する目標が有する音源から放射される音波を音響センサにより観測し、そこから得られる到来音波の方位及び周波数を用いて目標の運動諸元を解析する目標運動解析方法及び装置に関する。
【０００２】
【従来の技術】
一般に、目標運動の解析は、観測船に設けられた音響センサを使用して、目標が放射する音波を受信することにより、到来音波の方位及び周波数を観測し、観測された方位及び周波数、観測時刻及び観測位置の時系列データに基づいて目標の状態ベクトル（距離、方位、速力、針路、放射音波の周波数）を解析するという方法により行われている。
【０００３】
前述した解析方法は、状態ベクトルから推定される各観測時刻の解析方位及び解析周波数と、観測方位及び周波数との残差二乗を観測誤差分散により正規化し、これらを時系列的に合計した評価関数が最小となるような状態ベクトルを解析する非線形最適化問題に帰着する。
【０００４】
前述のような非線形最適化問題を解析する方法としては、これまで、カルマンフィルタや最小二乗法によるものが知られており、カルマンフィルタを用いた目標運動解析については、非特許文献１等に記載されたものが、また、最小二乗法を用いた目標運動解析については、非特許文献２等に記載されたものが知られている。
【０００５】
【非特許文献１】
ＩＥＥＥ会報１９８５年 Ｖｏｌ ＡＣ−３０ Ｎｏ １０ 「Ａ Ｓｔｏｃｈａｓｔｉｃ Ａｎａｌｙｓｉｓ ｏｆ ａ Ｍｏｄｉｆｉｅｄ Ｇａｉｎ Ｅｘｔｅｎｄｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ ｗｉｔｈ Ａｐｐｌｉｃａｔｉｏｎ ｔｏ Ｅｓｔｉｍａｔｉｏｎ ｗｉｔｈ Ｂｅａｒｉｎｇ Ｏｎｌｙ Ｍｅａｓｕｒｅｍｅｎｔｓ」
【０００６】
【非特許文献２】
ＩＥＥＥ会報１９８４年 Ｖｏｌ ＡＣ−２９ Ｎｏ ９ 「Ｆｕｎｄａｍｅｎｔａｌ Ｐｒｏｐｅｒｔｉｅｓ ａｎｄ Ｐｅｒｆｏｒｍａｎｃｅ ｏｆ Ｃｏｎｖｅｎｔｉｏｎａｌ Ｂｅａｒｉｎｇ−ｏｎｌｙ Ｔａｒｇｅｔ Ｍｏｔｉｏｎ Ａｎａｌｙｓｉｓ」
【０００７】
【発明が解決しようとする課題】
前述した従来技術は、何れも、既に観測された到来音波の方位時系列データ及び周波数時系列データを用いて目標運動諸元を解析する技術であるが、目標が遠距離で時系列データの変化率が小さく、時系列データ誤差が大きい場合で、時系列データが十分蓄積されていない場合に目標運動諸元を一義的に決定する条件である可観測性が成立せず、解析が不能となるという問題点を有している。
【０００８】
また、前述した従来技術は、何れも、観測船が運動を実施しても適切な運動でなければ解析不能な状態から脱することができないという問題点を有している。
【０００９】
図９は目標が遠距離で時系列データの変化率が小さく、時系列データの誤差が大きい場合で、時系列データが十分蓄積されていない場合に解析が不能となるという問題点を説明する図、図１０は観測船の運動が適切でない場合に解析不能な状態から脱することができないという問題点について説明する図であり、これらの図を参照して、前述した問題点について説明する。
【００１０】
図９において、音響センサの受信可能範囲が９０２として示すような範囲であり、観測船９０１が停止状態にあるとし、さらに、周波数時系列データが得られないものとする。このような場合、方位時系列データの観測誤差９０３が、図示のように大きい場合、観測結果から得られる目標の運動ベクトルが誤差の範囲となってしまい、目標の運動ベクトル（真値）９０５と目標の運動ベクトル解析値（誤値）９０４、９０６との区別がつかないことになる。
【００１１】
また、観測船が運動を実施しても適切な運動でない場合、例えば、図１０に示すように、観測船１００１が目標と並行運動して時系列データを蓄積しても、目標運動ベクトル（真値）１００４と目標運動ベクトル解析値（誤値）１００３、１００５との区別がつかない。
【００１２】
さらに、従来から知られている目標運動解析装置は、パッシブ音響センサで受信した音波の到来方位及び周波数が海中雑音による誤差やセンサ特有の測定誤差を有していたため、解析精度を一定値以下に収束させるのに長時間を要するという問題点を有している。また、従来、方位変化率を最大化する観測船の運動パターンが解析収束時間を加速させると考えられていたが、音響センサの受信可能範囲や観測船の運動性能の制約から観測船をどのように操舵するのが最適かを指示する指針は存在していない。
【００１３】
本発明の目的は、前述した従来技術の問題点を解決し、目標物の運動解析収束時間を最小化することが可能に観測船の運動パターンを決定して、目標物の運動解析を行う目標運動解析方法及び装置を提供することにある。
【００１４】
【課題を解決するための手段】
本発明によれば前記目的は、水中及び水上を運動する目標が有する音源から放射される音波を観測船に設けた音響センサにより観測し、そこから得られる到来音波の方位時系列データ及び周波数時系列データを用いて目標の距離、方位、速力、針路を含む運動諸元を解析する目標運動解析方法において、目標及び観測船の運動パターンに関する予測計算を実施し、算出された目標及び観測船の予測運動パターンについて、ヤコビアン行列の可観測性を評価することにより、解析収束時間を最小化することができる観測船の運動パターンを指示することにより達成される。
【００１５】
前述した手段を備える本発明によれば、その結果としての観測船の運動パターンを観測船の操舵装置へ送出することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。
【００１６】
【発明の実施の形態】
以下、本発明による目標運動解析方法及び装置の実施形態を図面により詳細に説明する。
【００１７】
図１は本発明に適用する目標運動解析の原理を説明する図である。図１において、１０１は観測船、１０２は目標である。
【００１８】
観測船１０１は、自身も移動しながら音響センサを使用して、水上または水中を移動する船舶等の目標１０２が放射する音波を受信することにより、到来音波の方位Ｂ〜ｙ（なお、〜 は、本来直前の符号の上に付けて表記すべきものである。以下の説明で、＾ と記述するものもあるが、これも同様である）及び周波数ｆ〜 を観測する。そして、観測された方位及び周波数、そして観測時刻及び観測位置の時系列データに基づいて、目標１０２の状態ベクトルＸ（距離ｒ、方位Ｂｙ 、速力Ｍｔ 、針路Ｃｔ 、放射音波の周波数ｆ０）が解析される。図１０には、観測船１０１が移動しながら到来音波の方位Ｂ〜ｙ及び周波数ｆ〜 を観測し、それらの観測結果の時系列データから目標１０２の状態ベクトルＸ（距離ｒ、方位Ｂｙ 、速力Ｍｔ 、針路Ｃｔ 、放射音波の周波数ｆ）を得ることができることを示している。なお、状態ベクトルを示す各要素の時刻の順を各要素の符号の後に１、２、……、ｉ、ｌを付けて示している。
【００１９】
前述した解析は、状態ベクトルＸから推定される各観測時刻の解析方位Ｂｙ 及び解析周波数ｆと、観測方位Ｂ〜ｙ及び周波数ｆ〜 との残差二乗を観測誤差分散σＢ〜ｙ２ 及びσｆ〜 ２により正規化し、これらを時系列的に合計した評価関数Φ（Ｘ） が最小となるような状態ベクトルＸを解析する非線形最適化問題として行うことができ、図１１に示す式（１）、式（２）を解くことにより解析することができる。
【００２０】
図２は本発明の一実施形態による目標運動解析装置の構成例を示すブロック図である。図２において、２０１は音響センサ、２０２は航海センサ、２０３は操舵装置、２０４は位置演算器、２０５は観測データベース記憶器、２０６は状態ベクトル推定器、２０７は解析値表示器、２０８は最適運動演算器、２０９は最適運動表示器、２１０は目標運動解析装置、２１１はオペレータコンソールである。
【００２１】
図２に示す本発明の実施形態による目標運動解析装置２１０は、音響センサ２０１で受信した目標音源からの音波到来方位及び受信周波数と、自船のログ（船速）及びジャイロ（方位）を計測する航海センサ２０２からの計測値に基づき位置演算器２０４で算出された観測船位置とを観測データベース記録器２０５に、時刻を整合させた時系列の観測データベースとして蓄積する。状態ベクトル推定器２０６は、ここで蓄積した観測データベースに対して最小二乗法や適応型フィルタにより目標運動諸元に関する非線形解析を実施し、その解析過程を解析値表示器２０７に表示すると共に、現在時刻における目標運動ベクトル解析値を最適運動演算器２０８に送出する。最適運動演算器２０８は、状態ベクトル推定器２０６で解析された目標運動ベクトルと観測船の複数の運動パターンの候補に基づいて将来の予測計算を実施し、図１３に示す式（１６）を使用してフィッシャー情報行列値の増加状況を比較評価する。そして、最適運動演算器２０８は、式（１６）に示すフィッシャー情報行列値が最大となる観測船の運動パターンを最適運動表示器２０９に表示する。なお、図１３における式（１７）〜式（２０）は、式（１６）の各項を示す式である。
【００２２】
通常、観測船の操舵装置２０３は、オペレータコンソール２１１からの指示により作動するが、最適運動演算器２０８により算出された運動パターンを操舵装置２０３に送出することにより、状態ベクトル推定器２０６における可観測性が成立し、目標の運動解析が可能な状態を生成して、解析収束速度を加速させるような観測船の運動を実施させることが可能となる。
【００２３】
図３は観測船の最適運動演算の原理を説明する図である。図３において、観測船３０８の現在位置３０５とその時点における目標運動ベクトルの解析値３０３及び観測方位誤差３０９との関係から、目標の運動ベクトルが現在の解析値３０３を維持して運動することを前提に予測計算を行うことにより、予測運動ベクトル３０４を得ることができる。そして、観測船の運動パターンに関する複数の候補３０６について、それぞれの運動を実施した場合における図１３に示すフィッシャー情報行列値の式（１６）を算出して比較評価し、フィッシャー情報行列値の式（１６）が最大化される観測船の運動パターンを決定して操舵することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。また、目標の運動ベクトルの現在の解析値３０３と観測船の運動パターンに関する複数の候補３０６との関係から予測計算を行うことにより、予測運動ベクトル３０４を得る場合、音響センサの受信可能範囲３０７や誤差特性及び観測船３０８の運動性能を考慮することにより、より正確に観測船３０８の運動パターンを決定することができる。
【００２４】
なお、初回観測時点から目標の運動ベクトルの現在の解析値３０３が得られるまでは、音響センサの受信感度から目標距離を推定すると共に、図１１に示している式（３）を用いて、音波到来方位の変化方向に目標の運動ベクトル３０１を仮定し、前述した最適運動演算処理を実施する。
【００２５】
図４は最適運動演算器２０８での処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図４には２つのフローを示しているが、ステップ４０１〜４０３を持つフローは、もう一方のフローにおけるステップ４１１の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００２６】
（１）処理が開始されると、まず、図１３に示す式（１６）の各観測船の運動パターン番号毎のフィッシャー情報行列値を初期化（ＦＩＭ（０）＝０に初期化）し、観測船の運動パターン番号を１に初期化（現在の運動パターン）する（ステップ４０４〜４０７）。
【００２７】
（２）次に、観測船の運動パターン番号ｎと目標運動ベクトルの解析値とに基づいて将来の仮想位置の予測計算を実施し、音響センサの受信可能範囲３０７で受信可能な仮想観測データを図１３に示す式（１８）を用いて生成し、さらに、生成した仮想観測データを仮想観測データベースに蓄積してデータベースを更新する（ステップ４０８〜４１０）。
【００２８】
（３）ここで、蓄積された仮想観測データベースの情報に基づいて、図１３に示す式（１６）フィッシャー情報行列値の計算を実施する（ステップ４１１）。
【００２９】
（４）フィッシャー情報行列値の計算の処理において、まず、図１３に式（１７）として示す予測ヤコビアン行列及び図１３に式（２０）として示す予測重み行列を計算し、さらに、図１３に示す式（１６）のフィッシャー情報行列値を算出する（ステップ４０１〜４０３）。
【００３０】
（５）次に、観測船の運動パターン番号ｎに関するフィッシャー情報行列値ＦＩＭ（ｎ）が観測船の運動パターン番号ｎ−１に関するフィッシャー情報行列値ＦＩＭ（ｎ−１）よりも大きくなるか否かを判定し、大きくなる場合、最適運動パターン番号をｎに置換える（ステップ４１２、４１３）。
【００３１】
（６）ステップ４１２の判定で、観測船の運動パターン番号ｎに関するフィッシャー情報行列値が観測船の運動パターン番号ｎ−１に関するフィッシャー情報行列値よりも大きくならない場合、また、ステップ４１３の処理後、観測船の全ての運動パターン番号に関する比較評価が完了しているか否か、すなわち、運動パターンｎが最大値を越えたか否かを判定し、完了していれば、ここでの処理を終了し、完了していなければ、運動パターン番号をインクリメントして、ステップ４０８からの処理に戻って、処理を繰り返す（ステップ４１４、４１５）。
【００３２】
図５は Ｍａｒｑｕａｒｄｔ法を利用した最小二乗法による状態ベクトル推定器２０６の処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図５には２つのフローを示しているが、ステップ５０１〜５０４を持つフローは、もう一方のフローにおけるステップ５０６の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００３３】
（１）疑似線形解析の処理の実行が開始されると、まず、データベース内の全観測時刻のデータを使用して、初期値算出の処理が実行される（ステップ５０５、５０６）。
【００３４】
（２）初期値算出の処理において、まず、観測データベースに蓄積されている全観測時刻のデータを使用し、観測誤差が微少であることを仮定した疑似線形解析の処理を行って初期解析値を算出する。この解析値は、ステップ５０５での非線形解析処理により反復修正されて最終的な解析値として算出される（ステップ５０１）。
【００３５】
（３）前回時刻における非線形解析処理５０５による解析値が「解有り」として存在するか否かを判定し、存在した場合、これを今回の初期解析値とし、疑似線形解析処理をスキップする（ステップ５０２）。
【００３６】
（４）ステップ５０２の判定で、前回の処理での解析値が存在していなかった場合、目標運動諸元の状態ベクトルを図１１に示す式（４）のように定義し、係数行列Ｆ及び定数行列Ｇのそれぞれを、図１１に示す式（６）、式（７）により計算する線形解析準備の処理を実行する。その後、図１１に示す式（５）により初期解析値を算出する（ステップ５０３、５０４）。
【００３７】
（５）次に、非線形解析処理で、目標運動諸元の状態ベクトルを、図１２に示す式（８）のように定義し、前述した疑似線形解析処理で算出した初期解析値とこれに基づく目標音源の放射周波数を、図１２に示す式（９）を用いて算出して設定する。
【００３８】
（６）次に、Ｍａｒｑｕａｒｄｔ 法解析における反復修正を安定的行うことを目的とした縮小因子λの初期値を設定する縮小因子λ設定の処理を実行し、さらに、反復修正回数の初期化（Ｌｏｏｐ１＝０）を行う（ステップ５０７、５０８）。
【００３９】
（７）次に、図１１に示している式（１）により、各時刻における観測情報と解析情報との残差の二乗和Φを算出する残差二乗和算出処理を実行し、さらに、図１２に示す式（１１）を用いて、方位情報及び周波数情報に関する偏微分係数を算出するヤコビアン行列形成の処理を実行する（ステップ５０９、５１０）。
【００４０】
（８）ステップ５０９の処理で算出した残差二乗行列とステップ５１０の処理で算出したヤコビアン行列とを用い、図１２に示す式（１０）により状態ベクトル修正量Δｘの算出を行う修正量Δｘ算出処理を実行する。なお、図１２における式（１１）〜式（１５）は、前述の式（１０）の各項を示す式である（ステップ５１１）。
【００４１】
（９）反復修正回数Ｌｏｏｐ１ が最大値を越えたか否か、あるいは、前回の残差の二乗和Φが今回の残差の二乗和Φより大きいか否かを判定し、判定結果がＮＧであった場合、縮小因子λを２倍にしてステップ５０９からの処理に戻って処理を繰り返す（ステップ５１２、５１８）。
【００４２】
（１０）ステップ５１２の判定で、判定結果がＯＫであり、判定条件を満足していた場合、状態ベクトルｘを、ｘ＋Δｘに修正し、反復修正回数Ｌｏｏｐ１ が１より大きいか否かを判定し、大きかった場合、縮小因子λを半分にして次回の反復修正処理に備える（ステップ５１３〜５１５）。
【００４３】
（１１）ステップ５１４の判定で、反復修正回数Ｌｏｏｐ１ が１より大きくなかった場合、あるいは、ステップ５１５の処理後、状態ベクトルの修正量Δｘが極小値ｍｉｎ より小さくなったか否か、反復修正回数Ｌｏｏｐ２ が上限値に達した否かを判定する。判定の結果がＮＧであった場合、ステップ５０９からの処理に戻って処理を繰り返す（ステップ５１６）。
【００４４】
（１２）ステップ５１６の判定で、状態ベクトルの修正量が極小値ｍｉｎ より小さくなるか、反復修正回数Ｌｏｏｐ２ が上限値に達した場合、反復を打ち切り、図１２に示す式（８）で求められた状態ベクトルを目標運動諸元（距離ｒ，方位Ｂｙ ，速力Ｍｔ ，針路Ｃｔ ，放射周波数ｆ０ ）に変換して処理を終了する（ステップ５１７）。
【００４５】
前述した処理において、ステップ５１６の判定で、反復修正回数Ｌｏｏｐ２ が上限値に達して終了した場合「解なし」として扱う。
【００４６】
なお、誤差を伴って観測される情報から内部の状態モデルを推定する最小二乗法及びＭａｒｑｕａｒｄｔ 法の詳細については以下の文献に記載されている。
【００４７】
「最小二乗法による実験データ解析」中川徹 他 東京大学出版会
図６はカルマンフィルタを利用した適応型フィルタによる状態ベクトル推定器２０６の処理動作を説明するフローチャートであり、次に、これについて説明する。なお、図６には２つのフローを示しているが、ステップ６０１〜６０３を持つフローは、もう一方のフローにおけるステップ６０５の処理の詳細を示すものであり、以下では、２つのフローを一連のものとして説明する。
【００４８】
ここでの処理は、目標運動諸元の状態ベクトルを、図１４に式（２１）として示すように定義し、観測データベースに蓄積されている各観測時刻毎のデータを１つ１つ使用して状態ベクトル及び共分散行列の初期値を逐次更新することにより解析値を収束させる処理である。
【００４９】
（１）逐次推定の処理の実行が開始されると、初期値設定の処理を実行する（ステップ６０４、６０５）。
【００５０】
（２）初期値設定の処理において、まず、前回の解析値が存在しない場合、図１４に示す式（２２）を用いて、状態ベクトルの初期値を設定する状態ベクトル初期化の処理を実行し、図１４に示す式（２３）を用いて、共分散行列の初期値を設定する共分散行列初期化の処理を実行する（ステップ６０１〜６０３）。
【００５１】
（３）次に、状態ベクトルの初期値を、観測データべースから取り出したデータの観測時刻まで、図１４に示す式（２４）により時間遷移させる状態ベクトル状態遷移処理を実行し、共分散行列状態遷移では、共分散行列の初期値を、観測データべースから取り出したデータの観測時刻まで、図１４に示す式（２５）により時間遷移させる共分散行列状態遷移の処理を実行する（ステップ６０６、６０７）。
【００５２】
（４）次に、図１４に示す式（２６）により、観測データベースから取り出した観測データに基づき状態ベクトルを修正更新する状態ベクトルフィルタリングの処理を実行し、図１４に示す式（２７）により、共分散行列を修正更新する共分散ベクトルフィルタリングの処理を実行する（ステップ６０８、６０９）。
【００５３】
（５）さらに、図１４に示す式（２８）により状態ベクトルを目標運動諸元に変換し、図１４に示す式（２９）により共分散行列を目標運動諸元解析値の解析誤差標準偏差に変換する目標運動諸元変換の処理を実行する（ステップ６１０）。
【００５４】
（６）最後に、観測データベースファイルがＥＯＦ（終了）となっているか否かを判定し、ＥＯＦとなっていた場合、ここでの処理を終了し、そうでない場合、ステップ６０６からの処理に戻り、観測データベースから取り出すデータがなくなるまで処理を繰り返す（ステップ６１１）。
【００５５】
なお、誤差を伴って観測される情報から内部の状態モデルを推定するカルマンフィルタの詳細については以下の文献に記載されている。
【００５６】
「応用カルマンフィルタ」片山徹 朝倉書店
前述した本発明の実施形態における各処理は、処理プログラムとして構成することができ、この処理プログラムは、ＨＤ、ＤＡＴ、ＦＤ、ＭＯ、ＤＶＤ−ＲＯＭ、ＣＤ−ＲＯＭ等の記録媒体に格納して提供することができる。
【００５７】
図７は観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に時系列方位データの変化率が最大化される方向へ針路を変更した場合における状態ベクトル推定器２０６の解析結果の状況を示す図である。
【００５８】
図７（ａ）には、観測船７０９の航跡と目標運動ベクトル（真値）７１１と解析目標位置７１０とを表示した平面図を示している。また、図７（ｂ）には、グラフ７０２〜７０６を示している。グラフ７０２には、横軸を時間軸とし、縦軸を距離とした目標距離（真値）７１２と解析距離７１３との時系列推移が表示されており、グラフ７０３には、横軸を時間軸とし、縦軸を速力とした目標速力（真値）７１４と解析速力７１５との時系列推移が表示されている。また、グラフ７０４には、横軸を時間軸とし、縦軸を針路とした目標針路（真値）７１６と解析針路７１７との時系列推移が表示されており、グラフ７０５には、横軸を時間軸とし、縦軸を方位とした観測船針路７０７と最適運動リコメンド７０８との時系列推移が表示されている。さらに、グラフ７０６には、横軸を時間軸としたフィッシャー情報行列値７１８の時系列推移が表示されている。
【００５９】
図８は観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に最適運動演算器の計算結果に従って針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【００６０】
図８（ａ）には、観測船８０９の航跡と目標運動ベクトル（真値）８１１と解析目標位置８１０とを表示した平面図を示している。また、図８（ｂ）には、グラフ８０２〜８０６を示している。グラフ８０２には、横軸を時間軸とし、縦軸を距離とした目標距離（真値）８１２と解析距離８１３との時系列推移が表示されており、グラフ８０３には、横軸を時間軸とし、縦軸を速力とした目標速力（真値）８１４と解析速力８１５との時系列推移が表示されている。また、グラフ８０４には、横軸を時間軸とし、縦軸を針路とした目標針路（真値）８１６と解析針路８１７との時系列推移が表示されており、グラフ８０５には、横軸を時間軸とし、縦軸を方位とした観測船針路８０７と最適運動リコメンド８０８との時系列推移が表示されている。さらに、グラフ８０６には、横軸を時間軸としたフィッシャー情報行列値８１８の時系列推移が表示されている。
【００６１】
図７と図８とに示す解析結果を比較してみると、観測船の運動パターンが図８（ａ）に示すように、フィッシャー情報行列を最大化する運動の場合、図７における方位変化率最大化の場合に比較してフッシャー情報行列値８１８が初期段階から増加傾向を示しており、目標距離の解析結果８１３及び目標速力の解析結果８１５の収束時間が加速されていることが判る。
【００６２】
前述した本発明の実施形態によれば、観測船の現在位置とその時点における目標運動ベクトルの解析値及び観測方位誤差との関係から、目標の運動ベクトルが現在の解析値を維持して運動することを前提に予測計算を行い、観測船の運動パターンに関する複数の候補について、それぞれの運動を実施した場合におけるフィッシャー情報行列値を算出して比較評価し、フィッシャー情報行列値を最大化することのできる観測船の運動パターンを決定することができる。そして、この運動パターンにより観測船を操舵することにより、可観測性が成立し解析可能な状態を生成して解析収束速度を加速させることが可能となる。
【００６３】
また、前述した本発明の実施形態によれば、目標の運動ベクトルの現在の解析値と観測船の運動パターンに関する複数の候補との関係を用いて予測計算を行う場合に、音響センサの受信可能範囲や誤差特性及び観測船の運動性能を考慮することにより、より正確な観測船の運動パターンを決定することができる。
【００６４】
【発明の効果】
以上説明したように本発明によれば、目標物の運動解析収束時間を最小化することが可能に観測船の運動パターンを決定して、観測船を運動させることができる。
【図面の簡単な説明】
【図１】本発明に適用する目標運動解析の原理を説明する図である。
【図２】本発明の一実施形態による目標運動解析装置の構成例を示すブロック図である。
【図３】観測船の最適運動演算の原理を説明する図である。
【図４】最適運動演算器での処理動作を説明するフローチャートである。
【図５】Ｍａｒｑｕａｒｄｔ 法を利用した最小二乗法による状態ベクトル推定器の処理動作を説明するフローチャートである。
【図６】カルマンフィルタを利用した適応型フィルタによる状態ベクトル推定器の処理動作を説明するフローチャートである。
【図７】観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に時系列方位データの変化率が最大化される方向へ針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【図８】観測船の運動パターンとして、初回受信した音波到来方位の方向へ一定時間運動した後に最適運動演算器の計算結果に従って針路を変更した場合における状態ベクトル推定器の解析結果の状況を示す図である。
【図９】目標が遠距離で時系列データの変化率が小さく、時系列データの誤差が大きい場合で、時系列データが十分蓄積されていない場合に解析が不能となるという問題点を説明する図である。
【図１０】観測船の運動が適切でない場合に解析不能な状態から脱することができないという問題点について説明する図である。
【図１１】説明の中で使用される数式を示す図（その１）である。
【図１２】説明の中で使用される数式を示す図（その２）である。
【図１３】説明の中で使用される数式を示す図（その３）である。
【図１４】説明の中で使用される数式を示す図（その４）である。
【符号の説明】
１０１ 観測船
１０２ 目標
２０１ 音響センサ
２０２ 航海センサ
２０３ 操舵装置
２０４ 位置演算器
２０５ 観測データベース記憶器
２０６ 状態ベクトル推定器
２０７ 解析値表示器
２０８ 最適運動演算器
２０９ 最適運動表示器
２１０ 目標運動解析装置
２１１ オペレータコンソール
Claims (5)
- 水中及び水上を運動する目標が有する音源から放射される音波を観測船に設けた音響センサにより観測し、そこから得られる到来音波の方位時系列データ及び周波数時系列データを用いて目標の距離、方位、速力、針路を含む運動諸元を解析する目標運動解析方法において、目標及び観測船の運動パターンに関する予測計算を実施し、算出された目標及び観測船の予測運動パターンについて、ヤコビアン行列の可観測性を評価することにより、解析収束時間を最小化することができる観測船の運動パターンを指示することを特徴とする目標運動解析方法。
- 前記ヤコビアン行列の可観測性を評価は、音響センサの受信可能範囲及び誤差特性、観測船の運動特性を考慮しておこなわれることを特徴とする請求項１記載の目標運動解析方法。
- 前記ヤコビアン行列の可観測性を評価は、複数の観測船の運動パターンの候補の中からフィッシャー情報行列値を最大化する観測船の運動パターンを算出するように行われることを特徴とする請求項１または２記載の目標運動解析方法。
- 算出された観測船の最適運動パターンをオペレータの判断により操舵装置へ送出することを特徴とする請求項１、２または３記載の目標運動解析方法。
- 水中及び水上を運動する目標が有する音源から放射される音波を観測船に設けた音響センサにより観測し、そこから得られる到来音波の方位時系列データ及び周波数時系列データを用いて目標の距離、方位、速力、針路を含む運動諸元を解析する目標運動解析装置において、目標及び観測船の運動パターンに関する予測計算を実施し、算出された目標及び観測船の予測運動パターンについて、ヤコビアン行列の可観測性を評価することにより、解析収束時間を最小化することができる観測船の運動パターンを指示する最適運動演算器を備えることを特徴とする目標運動解析装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003179895A JP4192046B2 (ja) | 2003-06-24 | 2003-06-24 | 目標運動解析方法及び装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003179895A JP4192046B2 (ja) | 2003-06-24 | 2003-06-24 | 目標運動解析方法及び装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2005017032A true JP2005017032A (ja) | 2005-01-20 |
| JP4192046B2 JP4192046B2 (ja) | 2008-12-03 |
Family
ID=34181097
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2003179895A Expired - Lifetime JP4192046B2 (ja) | 2003-06-24 | 2003-06-24 | 目標運動解析方法及び装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4192046B2 (ja) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010249593A (ja) * | 2009-04-14 | 2010-11-04 | Hitachi Ltd | 目標運動解析方法、目標運動解析装置およびプログラム |
| JP2011141125A (ja) * | 2010-01-05 | 2011-07-21 | Nec Corp | 目標運動解析方法及び目標運動解析装置 |
| JP2013250277A (ja) * | 2013-08-14 | 2013-12-12 | Nec Corp | 目標運動解析方法及び目標運動解析装置 |
| JP2014025786A (ja) * | 2012-07-26 | 2014-02-06 | Hitachi Ltd | 目標運動解析方法及び目標運動解析装置 |
| CN110319838A (zh) * | 2019-07-09 | 2019-10-11 | 哈尔滨工程大学 | 一种自适应的运动姿态参考系统升沉测量方法 |
| WO2020050084A1 (ja) * | 2018-09-07 | 2020-03-12 | 国立研究開発法人宇宙航空研究開発機構 | 誘導制御プログラムを格納した記憶媒体 |
| CN115616602A (zh) * | 2022-10-14 | 2023-01-17 | 哈尔滨工程大学 | 基于被动声纳纯方位定位检测前跟踪算法的观测者最优机动策略 |
-
2003
- 2003-06-24 JP JP2003179895A patent/JP4192046B2/ja not_active Expired - Lifetime
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010249593A (ja) * | 2009-04-14 | 2010-11-04 | Hitachi Ltd | 目標運動解析方法、目標運動解析装置およびプログラム |
| JP2011141125A (ja) * | 2010-01-05 | 2011-07-21 | Nec Corp | 目標運動解析方法及び目標運動解析装置 |
| JP2014025786A (ja) * | 2012-07-26 | 2014-02-06 | Hitachi Ltd | 目標運動解析方法及び目標運動解析装置 |
| JP2013250277A (ja) * | 2013-08-14 | 2013-12-12 | Nec Corp | 目標運動解析方法及び目標運動解析装置 |
| WO2020050084A1 (ja) * | 2018-09-07 | 2020-03-12 | 国立研究開発法人宇宙航空研究開発機構 | 誘導制御プログラムを格納した記憶媒体 |
| CN110319838A (zh) * | 2019-07-09 | 2019-10-11 | 哈尔滨工程大学 | 一种自适应的运动姿态参考系统升沉测量方法 |
| CN115616602A (zh) * | 2022-10-14 | 2023-01-17 | 哈尔滨工程大学 | 基于被动声纳纯方位定位检测前跟踪算法的观测者最优机动策略 |
| CN115616602B (zh) * | 2022-10-14 | 2023-08-18 | 哈尔滨工程大学 | 基于被动声纳纯方位定位检测前跟踪算法的观测者最优机动策略的航向确定方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| JP4192046B2 (ja) | 2008-12-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chatzi et al. | The unscented Kalman filter and particle filter methods for nonlinear structural system identification with non‐collocated heterogeneous sensing | |
| CN111985093B (zh) | 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法 | |
| CN104155650B (zh) | 一种基于熵权值法点迹质量评估的目标跟踪方法 | |
| CN104596514B (zh) | 加速度计和陀螺仪的实时降噪系统和方法 | |
| CN110286357B (zh) | 一种基于水声探测的目标运动定位方法 | |
| KR20180091372A (ko) | 레이더의 표적 위치 추적 방법 | |
| JP4192046B2 (ja) | 目標運動解析方法及び装置 | |
| CN106200377A (zh) | 一种飞行器控制参数的估计方法 | |
| CN115964603B (zh) | 一种基于改进Kalman滤波的机动目标跟踪方法 | |
| CN108319570A (zh) | 一种异步多传感器空时偏差联合估计与补偿方法及装置 | |
| KR101303417B1 (ko) | 정보 처리 장치, 정보 처리 방법 및 기록 매체 | |
| KR20070036009A (ko) | 관성 항법 장치를 이용한 실시간 위치 검사 방법 및 시스템 | |
| WO2019163701A1 (ja) | システム同定装置、システム同定方法及び記録媒体 | |
| KR20170111582A (ko) | 도플러 비컨을 이용한 수중 이동체의 위치 및 속도 추정 방법 및 시스템 | |
| CN105701292B (zh) | 一种机动目标转弯角速度的解析辨识方法 | |
| US20240045001A1 (en) | System and method for fast magnetometer calibration using gyroscope | |
| Goutaudier et al. | Impulse identification technique by estimating specific modal ponderations from vibration measurements | |
| JP4931145B2 (ja) | 磁気式３次元位置検出装置、磁気式３次元位置検出方法、プログラム及び記録媒体 | |
| Kang et al. | A novel recursive modal parameter estimator for operational time-varying structural dynamic systems based on least squares support vector machine and time series model | |
| CN109274107B (zh) | 一种计及奇异值的低频振荡信号参数辨识方法 | |
| Yang et al. | Doppler-bearing tracking in the presence of observer location error | |
| Usayiwevu et al. | Information driven self-calibration for LiDAR-inertial systems | |
| JP3881967B2 (ja) | オペレータ介在型目標運動解析方法及び装置 | |
| US20210281361A1 (en) | Methods, apparatuses and systems for time delay estimation | |
| CN105133472B (zh) | 路面平整度测量方法及惯性式检测仪 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20051108 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20080222 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080226 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20080422 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20080902 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20080919 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110926 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4192046 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120926 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120926 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130926 Year of fee payment: 5 |
|
| EXPY | Cancellation because of completion of term |