JP5012168B2 - 目標状態量推定方法 - Google Patents

Description

本発明は、目標から放射される音を、運動可能な艦船等の移動体に取り付けた受波器センサアレイで受信し、雑音に乱された観測量から、移動目標の位置と速度に関する内部状態量の推定を行う目標状態量推定方法に関するものである。

従来、『非線形方程式を解いて目標の運動解析を行なう際の解の収束を速める。』技術として、『目標の放射音に対応する受信信号が、方位情報算出部２２と周波数情報算出部２３に入力され、観測方位角θｍ（ｔ）と観測周波数成分ν_km（ｔ）が算出され、初期値算出部２４及びマニューバ検出部２５に渡される。マニューバ検出部２５は、マニューバ発生時点での内部状態量の推定値ＶＸ₀ ^*の速度に関する成分をゼロとした内部状態量の推定値ＶＸ_d0*と推定情報行列ΣＶＸ_d0* ^-1を、初期値算出部２４に送る。初期値算出部は、推定値ＶＸ_d0*と行列ΣＶＸ_d0* ^-1を先見情報として、その情報と目標に対する観測量とで構成された擬似線形方程式を解き、解析初期値ＶＸ₀₀ ^*を求めて内部状態量推定部２６に送る。内部状態量推定部２６は、初期値ＶＸ₀₀ ^*を適用したと評価関数で目標の運動解析を行う。』というものが提案されている（特許文献１）。

特開平９−１３３７５０号公報（要約）

上記特許文献１に記載の目標運動解析方法では、目標はマニューバ検出時刻ｔ_dまで等速直線運動を続け、瞬時にマニューバを終了し、マニューバ検出時刻ｔ_dから新たな等速直線運動を開始することを仮定している。
このため、目標がマニューバを開始してから終了するまでに時間が掛かるような場合、マニューバ検出時刻ｔ_dまでの観測値に基づいて推定した時刻ｔ_dにおける目標状態量の推定値の位置成分を次ステップの先見情報に用いても、実際の目標位置は変化してしまうため、次ステップにおける推定の誤差を大きくするという課題があった。

そのため、目標がマニューバを開始してから終了するまでの時間がかかるような場合でも、マニューバ終了後の次ステップにおける推定誤差を低減することのできる目標状態量推定方法が望まれていた。

本発明に係る目標状態量推定方法は、移動する目標から放射される音を、移動体に取り付けた受波器で受信し、その目標の位置と速度に関する状態量を推定する方法であって、前記受波器で受信した信号に基づき前記目標の観測方位角とその誤差分散を算出する方位情報算出ステップと、前記方位情報算出ステップの結果に基づき前記目標の状態量を推定する状態量推定ステップと、前記方位情報算出ステップの結果に基づき前記目標のマニューバを検出するマニューバ検出ステップと、を有し、前記マニューバ検出ステップにおいてマニューバが検出された際には、マニューバ開始時刻およびマニューバ終了時刻を検出するステップと、マニューバ開始時刻における前記目標の推定状態量とその推定情報行列を、マニューバ終了時刻における前記目標の推定状態量とその推定情報行列に変換する変換ステップと、前記変換ステップで変換した前記目標の推定状態量とその推定情報行列に基づき構成した評価関数を最小化する状態量を求めることにより、前記目標の状態量を推定するステップと、を実行するものである。

本発明に係る目標状態量推定方法によれば、マニューバ開始時刻における目標の推定状態量とその推定情報行列を、マニューバ終了時刻における目標の推定状態量とその推定情報行列に変換し、変換後の各値を次ステップの先見情報として使用しているので、マニューバ中の目標位置の変化による推定への悪影響を抑制することができ、推定精度の向上を実現することが可能となる。

図１は、目標状態量推定方法における観測系及び運動系を示す幾何学的説明図である。
図１において、（Ｘ、Ｙ）は原点ｏの固定座標系、１は航走体、２は状態量を推定する対象（以下、「目標」）である。各記号の意味は以下の通りである。
（１）ｘ ₁（ｔ）：時刻ｔにおける航走体１の位置ベクトル
（２）ｘ ₂（ｔ）：時刻ｔにおける目標２の位置ベクトル
（３）ｒ（ｔ）：時刻ｔにおける航走体１と目標２間の距離‖ｘ ₂（ｔ）−ｘ ₁（ｔ）‖
（‖‖はベクトルのノルム）
（４）θ（ｔ）：時刻ｔにおける航走体１から見たＹ軸を基準とする目標２の方位角

目標状態量推定は、目標２が等速直線運動を行っていると仮定し、雑音に乱された目標音源の方位角の観測時系列から、目標２の位置と速度に関する状態量を推定することを目的とする。

ここで、本発明の理解を容易にするため、まず目標状態量推定方法の原理について説明し、次にその原理を適用した従来の目標状態量推定装置の構成を図３で説明する。

時刻ｔにおける方位角θ（ｔ）の観測値、および音源の第ｋ固有周波数成分ｆ_kのドップラシフト周波数ν_k（ｔ）の観測値を（式１）のように表す。

次に、目標に対する観測者から見た相対位置ベクトルを次の（式２）のように表す。

また、時刻ｔ₀における目標状態量として、次の（式３）のＸ ₂（ｔ₀）を定義する。

目標２が等速直線運動を行うと仮定すると、時刻ｔにおける目標状態量Ｘ ₂（ｔ）とＸ ₂（ｔ₀）の関係は、次の（式４）の状態遷移式で与えられる。

一方、時刻ｔにおける航走体１の状態量Ｘ ₁（ｔ）は既知であり、次の（式５）で表されるものとする。

このとき、目標状態量Ｘ ₂（ｔ₀）を推定することは、バッチ処理区間ｔ＝ｔ₁・・・ｔ_nにおける観測値の組から構成される、次の（式６）の評価関数Ｌ（Ｘ ₂（ｔ₀））を最小とするＸ ₂（ｔ₀）を求めることと等価である。

ここで、航走体１および目標２の運動速度は音速に比べ遥かに小さいと仮定しており、観測雑音の誤差共分散行列Σを含む諸量は、次の（式７）で表される。

評価関数Ｌ（Ｘ ₂（ｔ₀））の最小化は、例えばガウス−ニュートン法等の非線形最小二乗法を用いて行うことができる。

観測量時系列時刻ｔ＝ｔ₁・・・ｔ_nの区間データに対する目標状態量Ｘ ₂（ｔ₀）の推定値の推定誤差共分散行列は、次の（式８）で定義される情報行列の逆行列で与えられる。

（式８）の情報行列は、目標状態量の真値Ｘ ₂（ｔ₀）に対して定義されており、実際の目標状態量推定では使用することができない。
そこでまず、目標状態量Ｘ ₂（ｔ₀）の推定値を用いて、その推定値に対する情報推定行列を次の（式９）のように定義する。

（式９）の推定情報行列の逆行列を、目標状態量Ｘ ₂（ｔ₀）の推定値に対する推定誤差共分散行列として扱う。

時刻ｔ＝ｔ₀の推定状態量に基づく時刻ｔの推定状態量は、（式１０）で定義される。

また、（式９）で得られる推定情報行列に基づく時刻ｔの推定情報行列は、（式１１）で定義される。

ところが、目標２が変針あるいは変速等のマニューバを行った場合には、目標状態量推定の仮定条件が成立しなくなる。
例えば、時刻ｔ＝ｔ_dでマニューバが検出された場合、時刻ｔ_dにおける目標２の位置および固有周波数に関する推定状態量とその推定情報行列は、時刻ｔ_d以降の状態量推定における先見情報として考えられる。

よって、（式１０）でｔ＝ｔ_dを代入して算出した推定状態量と、（式１１）でｔ＝ｔ_dを代入して算出した推定情報行列のそれぞれについて、速度に関する成分をゼロとしたものを求め、これらの値を初期値として、以後の状態量推定を行う。
即ち、マニューバ検出時刻ｔ_dを改めてｔ＝ｔ₀とし、（式１２）の評価関数を最小化するＸ ₂（ｔ₀）を求め、目標状態量Ｘ ₂（ｔ）を推定する。

以上が目標状態量推定の原理である。次に、この原理を用いた目標状態量推定装置の構成と動作について図３で説明する。

図３は、従来の目標状態量推定装置の構成と動作を説明するものである。
図３において、１１は受波器センサからの信号入力端子、１２は受信信号から目標信号の到来方位の算出を行う方位情報算出部、１３は受信信号から目標信号の周波数成分の算出を行う周波数情報算出部、１４は状態量推定の推定初期値を設定する初期値設定端子、１５は観測量時系列から目標のマニューバ検出を行うとともに、推定状態量及び推定情報行列を算出するマニューバ検出部、１６は目標の状態量を推定する目標状態量推定部、１７は目標状態量推定結果を出力する出力端子である。

次に、この目標状態量推定装置の動作についてステップ毎に説明する。

（１）信号入力端子１１から入力される受信センサアレイで受信された信号源からの受信信号は、それぞれ方位情報算出部１２および周波数情報算出部１３に送られる。

（２）受信信号が方位情報算出部１２および周波数情報算出部１３に入力されると、方位情報算出部１２は観測方位角θｍ（ｔ）およびその誤差分散σ_θ ²（ｔ）を算出し、周波数情報算出部１３はドップラシフト観測周波数成分ν_km（ｔ）およびその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）を算出する。
算出結果は、マニューバ検出部１５と目標状態量推定部１６に出力される。

（３）目標状態量推定部１６は、初期値設定端子１４から設定された推定初期値Ｘ ₀₀を状態量推定の初期値とし、（式６）の評価関数を最小化する目標状態量推定を行って目標状態量の推定値を求め、マニューバ検出部１５に送るとともに、その推定値を用いて、（式９）で推定情報行列を求めて、マニューバ検出部１５に送る。
さらに、（式４）を用いて目標状態量Ｘ ₂（ｔ）の推定値を求め、出力端子１７を介して出力する。

（４）次ステップにおいて、各時刻ｔ＝ｔ_n+1・・・ｔ_m（ｔ_m＞ｔ_n）の（ｍ−ｎ）組の観測方位角θｍ（ｔ）とその誤差分散σ_θ ²（ｔ）、および、ドップラシフト観測周波数成分νｋｍ（ｔ）とその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）が、マニューバ検出部１５に入力されると、マニューバ検出部１５では、次の（式１３）で定義される推定残差εの統計分布を調査し、その調査結果でマニューバ検出を行う。

（５）マニューバ検出部１５でマニューバが検出された場合、そのマニューバ検出部１５は、（式１０）および（式１１）を用いて、マニューバ検出時点ｔ＝ｔ_dの推定状態量および推定情報行列を算出し、速度に関する成分をゼロと置いた目標状態量の推定値およびその推定情報行列を、それぞれ目標状態量推定部１６に送る。

（６）目標状態量推定部１６は、初期値設定端子１４から設定された推定初期値Ｘ ₀₀を状態量推定の初期値とし、（式１２）の評価関数を最小化する目標状態量推定を行う。
そして、目標状態量推定部１６は、（式４）を用いて目標状態量Ｘ₂（ｔ）の推定値を出力端子１７を介して出力する。

（７）マニューバ検出部１５でマニューバが検出されなかった場合、目標状態量推定部１６は、推定初期値Ｘ ₀₀の代わりに、ステップ（３）で推定した時刻ｔ＝ｔ₀における推定状態量を状態量推定の初期値とする。
次に、各時刻ｔ＝ｔ₁・・・ｔ_mの（ｍ）組の観測方位角θ_m（ｔ）とその誤差分散σ_θ ²（ｔ）、およびドップラシフト観測周波数成分ν_km（ｔ）とその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）に対して（式６）の評価関数を最小化する目標状態量推定を行って目標状態量の推定値を求め、マニューバ検出部１５に送る。
また、（式４）を用いて目標状態量Ｘ ₂（ｔ）の推定値を求め、出力端子１７を介して出力する。

以上、従来の目標状態量推定装置の構成と動作について説明した。次に、本発明に係る目標状態量推定方法の原理について説明する。

目標２がマニューバすることによる針路、速力は未知であるので、目標がマニューバする前までに推定していた目標状態量の推定値に対して、目標がマニューバすることで針路Ｃ、速力Ｍの誤差が増えたと近似する。
図１に示すように、目標速度ベクトルは針路Ｃ、速力Ｍの関数となるので、針路、速力の誤差の誤差分散をそれぞれσ_C ²、σ_M ²とすると、目標速度ベクトルの誤差共分散行列は、（式１４）で算出できる。

目標のマニューバ開始時刻をｔ_s、マニューバ終了時刻をｔ_eとする。目標がマニューバすることによって、時刻ｔ_sの目標状態量の推定誤差共分散行列の速度に関する成分が、上述の（式１４）で求められる誤差共分散行列の分だけ増えたと考えられる。
このように速度に関する誤差を増やした推定誤差共分散行列は、次の（式１５）で算出する。

推定誤差共分散行列が上記（式１５）で求められる目標状態量に基づいて、時刻ｔ_eの目標状態量を推定すると、その推定誤差共分散行列は次の（式１６）で算出できる。

また、目標マニューバ前後で目標距離の変化は小さいと仮定して、時刻ｔ_eの目標距離を時刻ｔ_sの目標距離で近似すると、時刻ｔ_eの目標位置は（式１７）で算出できる。

ただし、

である。

よって、（式１６）で算出される推定誤差共分散行列の速度に関する成分をゼロとおいたものを初期値として、以後の状態量推定を行う。
即ち、マニューバ終了時刻ｔ_eを改めてｔ＝ｔ₀とし、（式１９）の評価関数を最小化するＸ ₂（ｔ₀）を求めることにより、目標２がマニューバする前の情報を活かして、目標状態量Ｘ ₂（ｔ）を推定する。

以上が本発明に係る目標状態量推定の原理である。次に、この原理を用いた目標状態量推定装置の構成と動作について図２で説明する。

図２は、本発明に係る目標状態量推定装置の構成と動作を説明するものである。
図２において、２１は受波器センサからの信号入力端子、２２は受信信号から目標信号の到来方位および方位の誤差分散の算出を行なう方位情報算出部、２３は受信信号から目標信号の周波数成分および周波数の誤差分散の算出を行う周波数情報算出部、２４は状態量推定の推定初期値を設定する初期値設定端子、２５は観測量時系列から目標のマニューバ検出を行うマニューバ検出部、２６は目標の状態量を推定する目標状態量推定部、２７は目標状態量推定結果を出力する出力端子、２８は推定状態量及び推定情報行列の算出及び変換を行う情報量変換部、２９は情報量変換部２８に対して、針路誤差分散、速力誤差分散を設定する設定端子である。

次に、本発明に係る目標状態量推定装置の動作についてステップ毎に説明する。

（１）信号入力端子２１から入力される受信センサアレイで受信された信号源からの受信信号は、それぞれ方位情報算出部２２および周波数情報算出部２３に送られる。

（２）受信信号が方位情報算出部２２および周波数情報算出部２３に入力されると、方位情報算出部２２は観測方位角θ_m（ｔ）およびその誤差分散σ_θ ²（ｔ）を算出し、周波数情報算出部２３はドップラシフト観測周波数成分ν_km（ｔ）およびその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）を算出する。
算出結果は、マニューバ検出部２５と目標状態量推定部２６に出力される。

（３）目標状態量推定部２６は、初期値設定端子２４から設定された推定初期値Ｘ ₀₀を状態量推定の初期値とし、（式６）の評価関数を最小化する目標状態量推定を行って目標状態量の推定値を求め、マニューバ検出部２５に送る。
さらに、（式４）を用いて目標状態量Ｘ ₂（ｔ）の推定値を、出力端子２７を介して出力する。

（４）次ステップにおいて、各時刻ｔ＝ｔ_n+1・・・ｔ_m（ｔ_m＞ｔ_n）の（ｍ−ｎ）組の観測方位角θ_m（ｔ）とその誤差分散σ_θ ²（ｔ）、および、ドップラシフト観測周波数成分ν_km（ｔ）とその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）が、マニューバ検出部２５に入力されると、マニューバ検出部２５では、（式１３）で定義される推定残差εの統計分布を調査し、その調査結果でマニューバ開始時刻ｔ_s、マニューバ終了時刻ｔ_eの検出を行う。

（５）マニューバ検出部２５でマニューバが検出された場合、マニューバ検出部２５は、時刻ｔ＝ｔ₀の目標状態量の推定値およびその推定情報行列から、（式１０）および（式１１）を用いて、マニューバ開始時刻ｔ＝ｔ_sの推定状態量および推定情報行列を算出し、情報量変換部２８に送る。

（６）情報量変換部２８は、マニューバ検出部２５から送られた推定状態量およびその推定情報行列と、設定端子２９から設定された針路誤差分散σ_C ²、速力誤差分散σ_M ²を用いて、（式１５）により、速度に関する誤差を増やした推定誤差共分散行列を算出する。
なお、針路誤差分散σ_C ²および速力誤差分散σ_M ²は、物理的に想定される範囲で適切に設定すれば良い。

（７）情報量変換部２８は、さらに、（式１６）を用いて、マニューバ終了時刻ｔ＝ｔｅの推定誤差共分散行列を算出した後、その逆行列を算出し、速度に関する成分をゼロとおいた推定情報行列を、推定情報行列の新たな初期値とする。
（８）一方、情報量変換部２８は、マニューバ開始時刻ｔ＝ｔ_sの推定状態量およびマニューバ終了時刻ｔ＝ｔ_eの観測方位角θ_m（ｔ_e）から、（式１７）を用いて時刻ｔ_eの目標位置を算出する。
（９）情報量変換部２８は、マニューバ開始時刻ｔ＝ｔ_sの推定状態量の成分のうち、位置に関する成分をステップ（８）で算出した値、速度に関する成分をゼロとおいた推定状態量を算出する。
（１０）情報量変換部２８は、ステップ（９）で算出した新たな推定状態量と、ステップ（７）で算出した新たな推定情報行列を、目標状態量推定部２６に送る。

（１１）目標状態量推定部２６は、初期値設定端子２４から設定された推定初期値Ｘ ₀₀を状態量推定の初期値とし、（式１９）の評価関数を最小化する目標状態量推定を行う。そして、目標状態量推定部２６は、（式４）を用いてＸ ₂（ｔ）の推定値を求め、出力端子２７を介して出力する。

（１２）マニューバ検出部２５でマニューバが検出されなかった場合、目標状態量推定部２６は、推定初期値Ｘ ₀₀の代わりに、ステップ（３）で推定した時刻ｔ＝ｔ₀における推定状態量を状態量推定の初期値とする。
次に、各時刻ｔ＝ｔ１・・・ｔ_mの（ｍ）組の観測方位角θ_m（ｔ）とその誤差分散σ_θ ²（ｔ）、およびドップラシフト観測周波数成分ν_km（ｔ）とその誤差分散σ_νk ²（ｔ）（ｋ＝１・・・ｐ）に対して、（式６）の評価関数を最小化する目標状態量推定を行って目標状態量の推定値を求め、マニューバ検出部２５に送る。
また、（式４）を用いて目標状態量Ｘ ₂（ｔ）の推定値を求め、出力端子２７を介して出力する。

以上、本発明に係る目標状態量推定装置の構成と動作について説明した。
なお、本発明は上記の実施の形態に限定されず、種々の変形が可能である。例えば、上述の説明では、推定する目標状態量を、目標２の位置、速度及び固有音源周波数としているが、目標２の位置と速度だけを推定する場合にも適用可能である。
この場合には、（式３）、（式４）、（式６）、（式９）、（式１２）から固有音源周波数及びドップラシフト周波数を取り除けばよい。

なお、以上の説明では、航走体１を例に説明したが、これに限られるものではなく、任意の移動体としてよい。

以上のように、本実施の形態によれば、目標のマニューバが検出されると、マニューバ開始時刻における目標状態量の推定値の位置に関する成分は、マニューバ終了時刻に置ける観測方位に基づいて位置を変換し、推定情報行列は目標針路、目標速力の誤差が増加したとして、マニューバ終了時刻に置ける推定情報行列に変換した後に先見情報として使用するようにしている。
このため、目標マニューバ中の目標位置の変化による推定への悪影響を抑制することができ、推定精度の向上を実現することが可能となる。

目標状態量推定における観測系及び運動系を示す幾何学的説明図である。 本発明に係る目標状態量推定装置の構成と動作を説明するものである。 従来の目標状態量推定装置の構成と動作を説明するものである。

符号の説明

１１ 信号入力端子、１２ 方位情報算出部、１３ 周波数情報算出部、１４ 初期値設定端子、１５ マニューバ検出部、１６ 目標状態量推定部、１７ 出力端子、２１ 信号入力端子、２２ 方位情報算出部、２３ 周波数情報算出部、２４ 初期値設定端子、２５ マニューバ検出部、２６ 目標状態量推定部、２７ 出力端子、２８ 情報量変換部、２９ 設定端子。

Claims (4)

1. 移動する目標から放射される音を、移動体に取り付けた受波器で受信し、その目標の位置と速度に関する状態量を推定する方法であって、
   前記受波器で受信した信号に基づき前記目標の観測方位角とその誤差分散を算出する方位情報算出ステップと、
   前記方位情報算出ステップの結果に基づき前記目標の状態量を推定する状態量推定ステップと、
   前記方位情報算出ステップの結果に基づき前記目標のマニューバを検出するマニューバ検出ステップと、
   を有し、
   前記マニューバ検出ステップにおいてマニューバが検出された際には、
   マニューバ開始時刻およびマニューバ終了時刻を検出するステップと、
   マニューバ開始時刻における前記目標の推定状態量とその推定情報行列を、マニューバ終了時刻における前記目標の推定状態量とその推定情報行列に変換する変換ステップと、
   前記変換ステップで変換した前記目標の推定状態量とその推定情報行列に基づき構成した評価関数を最小化する状態量を求めることにより、前記目標の状態量を推定するステップと、
   を実行することを特徴とする目標状態量推定方法。
2. 前記変換ステップにおいて、
   前記目標がマニューバする前の推定状態に対して、針路と速力に関する誤差が増加したものと仮定し、その増加分を反映するように前記変換を実行する
   ことを特徴とする請求項１に記載の目標状態量推定方法。
3. 前記目標のマニューバによる針路誤差分散と速力誤差分散を入力するステップを有し、
   前記変換ステップでは、
   前記針路誤差分散と前記速力誤差分散を用いて前記変換を実行する
   ことを特徴とする請求項２に記載の目標状態量推定方法。
4. 前記受波器で受信した信号に基づき前記目標のドップラシフト観測周波数とその誤差分散を算出する周波数情報算出ステップ、
   を有し、
   前記状態量推定ステップでは、
   前記方位情報算出ステップの結果と前記周波数情報算出ステップの結果に基づき前記目標の位置、速度、および固有音源周波数に関する状態量を推定し、
   前記マニューバ検出ステップでは、
   前記方位情報算出ステップと前記周波数情報算出ステップの結果に基づき前記目標のマニューバを検出する
   ことを特徴とする請求項１ないし請求項３のいずれかに記載の目標状態量推定方法。

Images