JP2009058367A - 目標運動解析方法及び装置 - Google Patents

Abstract

【課題】目標体の運動モデルを等速直進運動とすると、目標体が変針変速することにより、変針変速前後の目標の運動を正しく解析できないという問題があった。
【解決手段】本発明の目標運動解析方法は、目標体の変針変速を検出した上で、時系列観測データを用いて変針変速時刻を推定し、さらに時系列観測データを用いて、基準時刻における目標体位置座標及び、推定した変針変速時刻前後の等速直進区間毎の目標体の針路、速力を推定パラメータとする非線形最適化問題を解くことである。
【選択図】図8

Description

本発明は、目標体から放射される音波を時系列データとして観測し、目標の位置、針路、速力等を推定する目標運動解析方法及び装置に関する。

目標運動解析方法は移動可能な観測体に設けられた音響センサにより、目標体から放射される音波の到来方位及び周波数を時系列観測データとして蓄積し、前記時系列観測データと時系列推定データの残差が最小となる目標体の運動モデルのパラメータを推定するものである。

これまでの目標運動解析方法は目標体の運動モデルを等速直進運動として扱っており、推定するパラメータは基準時刻における目標体の位置座標、針路、速力、固有周波数であった(特許文献１)。

前述した時系列観測データと時系列推定データの残差を最小化し、目標体の運動ベクトルを推定する問題は、非線形最適化問題として知られており、カルマンフィルタによる方法、最小二乗法による方法、遺伝的アルゴリズムによる方法等があった。

また、等速直進運動ではない目標体、すなわち変針変速を行うような目標体については、変針変速を検知し、オペレータに提示する方法があった (特許文献2)。

特開平10−62511公報 特開平10−62507公報

前述した技術では、目標体の運動モデルを等速直進運動として扱っているため、目標体が変針変速を行うと、変針変速前後の目標の運動を正しく解析できないという問題があった。

図１１は前述した従来技術による目標運動解析のシミュレーション結果である。1101は観測体の航跡、1102は目標体の航跡、1103は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1104は観測体と目標体の相対距離、1105は目標運動解析結果による観測体と目標体の推定相対距離、1106は目標体の速力、1107は目標運動解析結果による目標体の推定速力、1108は目標体の針路、1109は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体変針以降においては目標運動解析を正しく実施できていないことが分かる。

また、目標体の変針変速を何らかの検定手段により看破し、変針変速後の時系列データを用いて、解析を行う方法においては、変針変速直後は解析が不能となり、解析結果が不安定となるという問題があった。これは、変針変速後は解析に使用するデータ数が減少するためである。

図１２は何らかの検定手段により、変針変速を検出し、変針変速検出後は変針変速時刻以降の観測データを用いて目標の運動を解析する技術による目標運動解析のシミュレーション結果である。1201は観測体の航跡、1202は目標体の航跡、1203は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1204は観測体と目標体の相対距離、1205は目標運動解析結果による観測体と目標体の推定相対距離、1206は目標体の速力、1207は目標運動解析結果による目標体の推定速力、1208は目標体の針路、1209は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体の変針以降は一度解析不能状態となり、約600秒間は解析結果が不安定な状態となっていることが分かる。

前述の課題を解決する方法として、目標体の運動モデルとして等速直進運動ではなく、変針変速を考慮して、2区間以上の等速直進区間の組み合わせとして目標体の運動をモデル化する方法が考えられる（以降、このモデルを等速直進モデルに対して、変針変速モデルと呼ぶ）。しかし、変針変速モデルの時系列推定方位関数及び時系列周波数関数は不連続なものとなってしまい、カルマンフィルタあるいは最小二乗法のように評価値の導関数を使用する最適化手法を適用することができない。遺伝的アルゴリズムへの変針変速モデルの適用は比較的容易であるが、遺伝的アルゴリズムは推定パラメータが増加すると、解析を収束させるための計算量は飛躍的に増加することが知られている。

そのため、計算量の大幅な増加を伴わずに変針変速前後の目標体の運動ベクトルを安定して出力する目標運動解析方法及び装置が求められている。

本発明は、目標体から放射される音波の到来方位及び周波数を音響センサにより観測した、時系列観測方位データ及び時系列観測周波数データから、目標体の位置及び針路、速力を推定する目標運動解析方法である。その目標運動解析方法は、目標体の変針変速を検出した上で変針変速時刻を推定し、時系列観測データを用いて、基準時刻における目標体位置座標及び、等速直進区間毎の目標体の針路、速力を推定パラメータとする非線形最適化問題を解くことである。

より具体的な態様は、まず目標体の変針変速の時刻を仮定し、仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、その等速直進運動区間毎に、時系列観測方位データを基にした時系列推定方位データと時系列観測周波数データを基にした時系列推定周波数データとを求め、時系列観測方位データと時系列推定方位データとの残差及び時系列観測周波数データと時系列推定周波数データとの残差の和を評価値として、変針変速の時刻を推定する。この推定した変針変速の時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、この等速直進運動区間毎に、時系列観測方位データを基にした時系列推定方位データと時系列観測周波数データを基にした時系列推定周波数データとを求め、時系列観測方位データと時系列推定方位データとの残差及び時系列観測周波数データと時系列推定周波数データとの残差の和を評価値として、目標体の運動ベクトルを推定する。

本発明の目標運動解析方法及び装置により、変針変速を繰り返すジグザグ運動をする目標体であっても、目標体の運動ベクトルを安定して推定することができる。

以下、本発明による目標運動解析方法及び装置の実施の形態について図面を参照して説明する。

図１は、目標運動解析方法を実現するための機能ブロック図である。変針変速検出部102は、観測情報データベース101と変針変速データベース103を参照し、新たな変針変速があったかを検出すると共に、変針変速時刻を推定する。目標運動解析部104は観測情報データベース101と変針変速データベース103を参照し、基準時刻における目標体位置と等速直進毎の針路、速力を推定する解析を行い、解析履歴データベース105に解析結果を格納する。変針変速データベース103上に変針変速が登録されていない場合には、基準時刻における目標位置と速力を推定する解析を行い、解析結果データベース105に解析結果を格納する。

図２は観測情報データベース部101を説明する図である。図に示すように、観測方位データベースと観測周波数データベースは独立に構成する。[]は時系列にデータを蓄積するために配列構造を示しており、蓄積されているデータ数を格納データ数として格納しておく。なお、観測方位誤差標準偏差及び観測周波数誤差標準偏差は観測した情報の正確さを示すものであり、必ずしも必要な情報ではないが、これを用いて解析することにより、より高精度の解析が可能となる。

図３は変針変速データベース103の構造を説明する図である。変針変速情報として推定変針変速時刻を格納するレコードを持ち、変針変速回数を保持するレコードを持つ構造となっている。

図４は変針変速検出部102の処理を説明するフローチャートある。401の処理では、変針変速データベース103から最新の変針変速時刻を取得し、観測情報データベース101中の最新変速時刻以降のデータ、すなわち最新の区間のデータを取得し（以降、解析用データと記述する）、以降の処理を実施する。ここで、最新の変針変速時刻とは、以前に推定した変針変速時刻の中で、最も新しい時刻データである。したがって、最新の区間のデータとは、最新の変針変速時刻から最も新しい観測時刻までの時系列の観測情報（解析用データ）である。

402の変針変速検出処理は解析用データから変針変速の有無を統計的な手段で検出する処理であり、多数の特許が出願(たとえば、特許文献2)されているため、ここでの説明を省略する。変針変速検出処理の結果を403にて判定し、変針変速が検出された場合には404以降の処理により変針変速時刻を推定する。

考え方は最新の観測データ時刻から最新の変針変速時刻まで、変針変速時刻を仮定値として2レグ解析を実施して、最も評価値の良い仮定値を推定変針変速時刻とするものである。2レグ解析の詳細については後述するが、２つの等速直進区間の時系列観測データを用いて、基準時刻における目標体の位置座標、変針変速前後の速力、固有周波数を推定し、評価値として正規化方位残差と正規化周波数残差の和を出力するものである。

まず、処理404において最良の評価値Φに初期値Φ_initを設定する。評価値は正規化方位残差と正規化周波数残差であるため、小さい方が良い評価値である。したがって、Φ_initは十分に大きな値とする。次に処理405において変針変速時刻の仮定値t_zigに最新の解析用データの観測時刻t_newを設定する。処理406の2レグ解析は変針変速時刻を仮定値t_zigとして基準時刻における目標体の位置座標、変針変速前後の速力、固有周波数を推定し、評価値φを出力する。次に分岐407において、最良の評価値Φと２レグ解析にて計算した評価値φを比較し、2レグ解析にて計算した評価値φの方が良ければ、処理408及び処理409にて最良の評価値Φと、最良の評価値となる推定変針変速時刻T_zigを更新する。処理410にて変針変速時刻の仮定値をdtだけ過去に遡り、分岐411にて変針変速最古時刻T_maxより過去でなければ、処理406以降を繰り返す。最終的に変針変速時刻はT_zigとなり、処理412にて変針変速データベース103に格納する。

変針変速最古時刻T_maxは、最新の区間の中の最も古い時刻の意であり、前述の最新の変針変速時刻である。図４の処理が時間を遡るように繰り返すために、変針変速時刻の仮定値をdtだけ過去に遡るようにしているのであって、図４の処理が変針変速最古時刻T_maxから最新の観測時刻に向けて処理するならば、時間の推移は逆になることは自明であろう。またdtは、推定する変針変速時刻に要求される分解能より小さな値であればよい。

図５は解析履歴データベース105の構造を説明する図である。解析結果は固有周波数とレグと呼ぶ等速直進区間の運動ベクトルで管理する。解析履歴データベースに格納しているレグ数は格納レグ数に保持する。レグには解析基準時刻、解析基準時刻における目標位置X座標、解析基準時刻におけるY座標、速力X成分、速力Y成分を格納する構造となっている。

図６は目標運動解析部104の処理を説明するフローチャートである。まず、分岐601において、変針変速データベース103を参照して、変針変速回数が０の場合は処理602の１レグ解析を実行し、変針変速回数が１以上の場合は603の２レグ解析を実行する。

ここで、1レグ解析は、前述の最新の変針変速時刻から最新の観測時刻に至る区間の、目標体の運動モデルとして等速直進モデルを適用するものであり、解析基準時刻における目標体の位置座標、速力、固有周波数を推定するものである。２レグ解析は目標体の運動モデルとして、前述の最新の変針変速時刻から推定した変針変速時刻までの区間と推定した変針変速時刻から最新の観測時刻までの区間とのそれぞれの区間に等速直進モデルを適用するものであり、解析基準時刻における目標体の位置座標、固有周波数、等速直進区間毎の針路、速力を推定するものである。

最後に、処理604にて解析結果を解析履歴データベース105に格納する。

図７は１レグ解析において推定する運動パラメータを説明する観測体と目標体の幾何学的関係を示した図である。701は観測体の運動軌跡、702は目標体の運動軌跡を表している。703〜705は目標体から放射される音波を観測した観測体の位置、705が最新の観測データを取得した観測体の位置を表しており、最新観測データを取得した時刻を解析基準時刻とする。１レグ解析において推定するパラメータは解析基準時刻における目標体位置X座標707、解析基準時刻における目標体位置Y座標708、速力X成分710、速力Y成分711及び目標体固有周波数である。

これらのパラメータは最小二乗法により式(1)を最小化することにより推定する。ここで、tは解析基準時刻を0として正規化した時刻を表し、Cは音速、co,moは観測体の針路、速力を表す。また、y_i、z_iは時刻t_iにおける観測方位及び観測周波数であり、式(2)(3)はそれぞれ推定方位、推定周波数である。また、xo_i、yo_iは時刻t_iにおける観測体の位置（xo_i，yo_i）である。

特に、式(2)(3)から構成される評価値(1)はパラメータx₁,x₂,x₃,x₄,x₅に対して非線形であるため、反復により解く必要があり、例えば、Levenberg-Marquardt法ではΔXが十分小さくなるまで式(7)(8)の計算繰り返し実施する。なお、λは定数、Ｉは単位行列である。またbは、ヤコビアン行列と残差ベクトルの積であり、ここでは、方位ヤコビアン行列の転置行列と方位残差ベクトルの積、および周波数ヤコビアン行列の転置行列と周波数残差ベクトルの積の二つの積の和である。

式(8)のA、Bは方位情報及び周波数情報のヤコビアン行列であり、次式となる。

図８はｎレグ解析において等速直進区間が２区間の場合に、推定する運動パラメータを説明する観測体と目標体の幾何学的関係を示した図である。801は観測体の運動軌跡、802は目標体の運動軌跡を表している。803は最古の観測体の位置座標、804は推定変針変速時刻における観測体の位置座標、805は最新の観測体の位置座標を表している。また806は、推定した変針変速時刻における目標体の位置であり、そのX座標を807で、そのY座標を808で表している。

解析基準時刻は最新の観測時刻とし、２レグ解析において推定するパラメータは解析基準時刻における目標体位置812のX座標813、解析基準時刻における目標体位置812のY座標814、変針変速前の速力809のX成分810、変針変速前の速力809のY成分811、変針変速後の速力815のX成分816、変針変速後の速力815のY成分817及び固有周波数である。

これらのパラメータは１レグ解析と同様に式(1)を最小化することにより推定する。なお、図４のフローチャートでは、２レグ解析における式(1)のS(X)をΦで表している。ただし、推定方位及び推定周波数は推定変針変速時刻をt_z（図４のフローチャートでは、T_zigと表記している。）として、式(21)〜(24)となる。

Levenberg-Marquardt法等により計算するためのヤコビアンは式(30)〜(56)となる。

実施例として、本発明による目標運動解析装置の構成例を図９に示す。目標運動解析装置901は図９に示すように、CPU904、メモリ905、出力装置906、I/F(インターフェイス)907および908、入力装置909から構成されるような電子計算機で実現可能であり、音響センサ902及び航海センサ903に接続される。目標運動解析装置のサイズや演算速度に対する制約により、専用装置として構成されても良い。

音響センサ902は目標体から放射される音波を受波して、音波の到来方位及び周波数を目標運動解析装置901に入力する。航海センサ903は電磁ログ、ジャイロ、GPS等により、観測体の位置、速力、針路を目標運動解析装置901に入力する。入力された音波の到来方位、周波数、観測体の位置、速力、針路はメモリ905の観測情報データベースに蓄積される。変針変速検出部及び目標運動解析部は前述した方法に基づく処理を行い、その結果を変針変速データベース及び解析履歴データベースに格納する。

メモリ905上で保有する観測情報データベース、変針変速データベース、解析履歴データベースの内容をCRT等の出力装置906に表示し、オペレータに提示する。オペレータは出力装置906に出力された観測情報データベース、変針変速データベース、解析履歴データベースの内容を参照することにより、目標の変針変速を看破した場合には入力装置909により変針変速時刻を変針変速データベースに格納することも可能である。

例えば、図１０のような横軸に時刻、縦軸に観測方位とした図を出力装置906に表示する。図１０は目標の運動パターンを表１のとおりとした場合に音響センサで観測される音波の到来方位であり、図１０からも700秒に観測方位の傾向が変化していることがわかる。

しかしながら、実際に観測される方位は音波伝搬中の屈折の影響や音響センサにおける計測誤差もあり、図１０のように鮮明な方位の傾向変化を読み取る事は難しいため、本装置は統計的な手段により変針変速を看破する変針変速検出部を備えている。

図１３は本実施例による目標運動解析方法のシミュレーション結果であり、1301は観測体の航跡、1302は目標体の航跡、1303は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1304は観測体と目標体の相対距離、1305は目標運動解析結果による観測体と目標体の推定相対距離、1306は目標体の速力、1307は目標運動解析結果による目標体の推定速力、1308は目標体の針路、1309は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体の変針以降も解析結果を安定して出力している。なお、本シミュレーションの条件は図１１、１２と同様であり、本実施例により目標体が変針変速している場合でも、安定した解析を実施することができていることが分かる。

目標運動解析を実施するためのブロック図である。 観測情報データベース部の構造を説明する図である。 変針変速データベース部の構造を説明する図である。 変針変速検出部のフローチャート図である。 解析履歴データベース部の構造を説明する図である。 目標運動解析部のフローチャート図である。 等速直進モデルの目標運動解析における観測体と目標体の幾何学的関係を示した図である。 変針変速モデルの目標運動解析における観測体と目標体の幾何学的関係を示した図である。 目標運動解析装置の構成例である。 目標が変針した場合の音波の到来方位を時系列に示したグラフである。 変針する目標体に対して、従来技術である等速直進モデルによる目標運動解析を適用した場合のシミュレーション結果である。 変針する目標体に対して、従来技術である変針変速検出後の観測データを用いて目標運動解析を適用した場合のシミュレーション結果である。 変針する目標体に対して、目標運動解析を適用した場合のシミュレーション結果である。

符号の説明

101…観測情報データベース部、102…変針変速検出部、103…変針変速データベース部、104…目標運動解析部、105…解析履歴データベース部、901…目標運動解析装置、902…音響センサ、903…航海センサ、904…ＣＰＵ、905…メモリ、906…出力装置、907〜908…インターフェース、909…入力装置。

Claims (6)

1. 目標体から放射される音波の到来方位及び周波数を音響センサにより時系列的に観測し、蓄積した時系列観測方位データ及び時系列観測周波数データを用いて、前記目標体の位置、針路及び速力を含む運動ベクトルを推定する目標運動解析方法において、
   (a)目標体の変針変速の検出に応じて、前記変針変速の時刻を仮定し、
   (b)前記仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、前記区間毎に、前記時系列観測方位データを基にした時系列推定方位データと前記時系列観測周波数データを基にした時系列推定周波数データとを求め、
   (c)前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、前記変針変速の時刻を推定し、
   (d)前記推定した変針変速の時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、前記区間毎に、前記時系列観測方位データを基にした時系列推定方位データと前記時系列観測周波数データを基にした時系列推定周波数データとを求め、
   (e)前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、前記目標体の前記運動ベクトルを推定する目標運動解析方法。
2. 前記仮定した変針変速の時刻を変化させて、前記(b)及び(c)の手順を繰り返す請求項１記載の目標運動解析方法。
3. 前記仮定した変針変速の時刻を変化させて、前記(b)及び(c)の手順を繰り返した結果、前記評価値が最小となる変針変速時刻を出力する請求項２載の目標運動解析方法。
4. 前記時系列観測方位データの出力に応じてオペレータによって入力される変針変速時刻を前記仮定した変針変速時刻とする請求項１記載の目標運動解析方法。
5. 目標体から放射される音波の到来方位及び周波数を時系列的に観測する音響センサ、
   前記音響観測船差により観測された前記音波の到来方位及び周波数を、時系列観測方位データ及び時系列観測周波数データとして格納するデータベース、及び
   (a)目標体の変針変速の検出に応じて、前記変針変速の時刻を仮定する処理、(b)前記仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、前記区間毎に、前記時系列観測方位データを基にした時系列推定方位データと前記時系列観測周波数データを基にした時系列推定周波数データを求める処理、(c)前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、前記変針変速の時刻を推定する処理、(d)前記推定した変針変速の時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、前記区間毎に、前記時系列観測方位データを基にした時系列推定方位データと前記時系列観測周波数データを基にした時系列推定周波数データとを求める処理、および(e)前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、前記目標体の位置、針路及び速力を含む運動ベクトルを推定する処理を実行するプロセッサを有する目標運動解析装置。
6. 前記プロセッサは、前記仮定した変針変速の時刻を変化させて、前記(b)及び(c)の処理を繰り返した結果、前記評価値が最小となる変針変速時刻を出力する請求項５載の目標運動解析装置。

Images