JPH1062508A

JPH1062508A - 目標運動解析方法

Info

Publication number: JPH1062508A
Application number: JP8222980A
Authority: JP
Inventors: Osamu Fujimoto; 治藤本; Yoshio Okita; 芳雄沖田
Original assignee: Oki Electric Industry Co Ltd
Current assignee: Oki Electric Industry Co Ltd
Priority date: 1996-08-26
Filing date: 1996-08-26
Publication date: 1998-03-06

Abstract

(57)【要約】【課題】目標からの音響又は電磁波信号を観測し、そ
の観測時刻が異なる観測誤差には相関性が有るとした処
理方法により、従来よりも状態量推定精度の向上した目
標運動解析方法。【解決手段】目標から放射される音響又は電磁波信号
を逐次受信、この受信信号から前記目標の方位角及びド
ップラ周波数を逐次算出し（１〜５）、この方位角及び
ドップラ周波数の時系列観測値に基づき前記目標の状態
量を推定する目標運動解析方法において、まず前記方位
角及びドップラ周波数の時系列観測値の誤差分散をそれ
ぞれ推定し、次に、前記方位角及びドップラ周波数の時
系列観測値及びその誤差分散をそれぞれ白色化する処理
を行い（６）、次に、前記方位角及びドップラ周波数の
それぞれの白色化処理後の時系列観測値及びその誤差分
散を用いて前記目標の状態量を推定する（７〜１１）方
法。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、例えば移動する目
標から放射される信号を観測位置の移動可能な観測体に
設けられた受波センサアレイで受信し、前記目標の運動
解析を行う方法に関するものである。

【０００２】

【従来の技術】移動する目標体から放射された信号を観
測位置の移動可能な観測体に取り付けた受波器センサア
レイで受信し、上記目標体の運動解析を行う目標運動解
析に関する公知文献としては、例えば「ＩＥＥＥ，ＩＣ
ＡＳＳＰ，ＣＨ2673-2，U2.5,1989，J,M.PASSERIEUX,D.
PILLON,P.BLANC-BENON,C.JAUFFRET"TARGET MOTION ANAL
YSIS WITH BEARINGS AND FREQUENCIES MEASUREMENTS VI
A INSTRUMENTAL VARIAVLE ESITMATOR" PP.2645-2648 」
がある。

【０００３】図２は目標運動解析方法が適用される観測
系及び運動系の幾何学的説明図であり、図において、
（Ｘ，Ｙ）は原点ｏの固定座標系、２１は観測体、２２
は目標体、ｘ₁(t) は時刻ｔにおける観測体２１の位置
ベクトル、ｘ₂(t) は時刻ｔにおける目標体２２の位置
ベクトル、ｒ(t) は時刻ｔにおける観測体２１と目標体
２２間の距離‖ｘ₂(t) −ｘ₁(t) ‖である。なお、‖
‖はベクトルのノルム（ベクトルの長さの概念）を表
す。また、θ(t) は時刻ｔにおける観測体２１から見た
Ｙ軸を基準とする目標体２２の方位角、ν_k(t) は時刻
ｔにおける目標体２２の持つ信号源のｋ番目の周波数ｆ
_k（以下、この周波数を固有周波数と呼ぶ）を観測体２
１で観測したときのドップラシフトした周波数（以下、
この周波数をドップラ周波数と呼ぶ）である。（なお、
観測体２１においては、ドップラ周波数しか観測でき
ず、固有周波数ｆ_kは目標運動解析によって推定され
る。）目標運動解析は、目標体２２が等速直線運動を行ってい
ると仮定し、雑音に乱された目標の方位角あるいは方位
角とドップラ周波数の観測時系列から、目標体２２の状
態量である位置及び速度を推定するものである。

【０００４】次に、従来の目標運動解析装置が採用して
いる目標運動解析方法の原理について説明する。観測値
は真値に観測誤差が加わったものであり、そのため、時
刻ｔにおける方位角の観測値β(t) 及び第ｋドップラ周
波数の観測値γ_k(t) は、次の（１），（２）式のよう
に表される。 β(t) ＝θ(t) ＋ｎ_β(t) …（１） γ_k(t) ＝ν_k(t) ＋ｎ_γk(t) …（２）ただし、ｎ_β(t) は時刻ｔにおける平均０、分散σ_β ²
(t) なる方位角の観測誤差、ｎ_γk(t) は時刻ｔにおけ
る平均０、分散（σ_γk）²(t) なる第ｋドップラ周波
数の観測誤差である。ここで、時刻ｔにおける目標体２
２に対する観測体２１から見た相対位置ベクトルｘ(t)
を（３）式で表し、観測体２１の状態量Ｘ(t) を、
（４）式に示すように、相対位置ベクトルの要素、相対
速度ベクトルの要素、第ｋドップラ周波数の逆数の要素
で定義する。

【０００５】

【数２】

【０００６】なお、相対位置ベクトルｘ(t) を、各座標
軸方向の成分で表記すると、状態量Ｘ(t) は、（５）式
のように定義することもできる。

【０００７】

【数３】

【０００８】ただし、Ｔは転置、′は時間での１階微分
を表すものである。このような状態量Ｘ(t) を、方位角
及びドップラ周波数の観測値時系列β及びγ_k（ｋ＝１
〜Ｋ）を用いて、推定して求めることを通じて、最終的
に出力すべき値（例えば相対位置ベクトルｘ(t) ）を得
る。なお、方位角の観測値時系列ベクトルβ及び第ｋド
ップラ周波数の観測値時系列γ_kは、それぞれ以下の
（６）式及び（７）式で表される。

【０００９】

【数４】

【００１０】時刻ｔにおける状態量Ｘ(t) の推定値は、
時刻ｔ＝ｔ₀における状態量Ｘ₀＝Ｘ（ｔ₀）（以下、
初期状態量と呼ぶ。）と定義すると、（８）式で定義さ
れる評価関数Ｌ（Ｘ₀）を最小とする初期状態量Ｘ₀を
求めた後、（９）式を適用することにより求める。な
お、下記における″は時間での２階微分を表すものであ
る。

【００１１】

【数５】

【００１２】ここで、ベクトルθ（Ｘ₀）、ベクトルν
_k（Ｘ₀）及び行列Φ（ｔ，ｔ₀）はそれぞれ、次の
（１０）式、（１１）式、（１２）式で表される。

【００１３】

【数６】

【００１４】なお、ベクトルθ（Ｘ₀）の各要素θ（ｔ
_N，Ｘ₀）は、時刻ｔ＝ｔ₀における初期状態量Ｘ
₀と、それからの時間経過で定まる時刻ｔ_Nでの相対位
置ベクトルｘ（ｔ_N，Ｘ₀）の要素ｒ_x（ｔ_N，Ｘ₀）
及びｒ_y（ｔ_N，Ｘ₀）を用いて、（１３）式で表さ
れ、ベクトルν_k（Ｘ₀）の各要素ν_k（ｔ_N，Ｘ₀）
は、時刻ｔ＝ｔ₀における初期状態量Ｘ₀と、それから
の時間経過で定まる時刻ｔ_Nでの相対速度ベクトルｘ′
（ｔ_N，Ｘ₀）及び方位角θ（ｔ_N）＝θ（ｔ_N，
Ｘ₀）を用いて、（１４）式で表される。

【００１５】

【数７】

【００１６】（１４）式におけるｃは観測体２１及び目
標体２２間での媒質における音速である。また、（８）
式におけるΣ_β及びΣ_γkはそれぞれ、方位角及び第ｋ
ドップラ周波数の観測値についての誤差共分散行列であ
り、方位角及び第ｋドップラ周波数の観測値についての
観測誤差ベクトルＮ_β及びＮ_γkをそれぞれ、（１５）
式及び（１６）式で表すと、（１７）式及び（１８）式
で定義されるものである。なお、Ｅ［］は期待値を表
している。

【００１７】

【数８】

【００１８】上記（８）式で表される評価関数Ｌ
（Ｘ₀）を最小にする初期状態量Ｘ₀の推定には、オン
ライン処理のカルマンフィルタやバッチ処理の最尤推定
法が用いられる。

【００１９】図３は従来の目標運動解析を実施するため
の目標運動解析装置の構成を示す図である。図におい
て、１は受波器センサアレイからの出力信号が入力され
る信号入力端子、２は周波数分析およびビームフォーミ
ングを行う周波数分析／ビームフォーミング処理部、３
は入力信号の積分を行う積分処理部、４はオペレータに
よる積分時間及び目標探知操作情報が入力される設定情
報入力端子、５は周波数分析／ビームフォーミングされ
た信号から目標の方位と周波数情報を得る方位／周波数
情報算出部、１０は推定結果の出力端子、１１は自船情
報入力部、１２は観測値とその誤差分散を用いて目標の
状態量を推定する状態量推定部である。

【００２０】図３の動作を説明する。信号入力端子１か
ら入力される受波器センサアレイで受信された目標信号
源からの受信信号は、周波数分析／ビームフォーミング
処理部２に送られる。周波数分析／ビームフォーミング
処理部２は、入力された信号に対して、観測体２１の全
周（３６０°）にＭ本のビームを形成するようにビーム
フォーミングし、その出力を２乗してビーム出力信号Ｓ
_θm(t) （ｍ＝１〜Ｍ）として積分処理部３に送ると同
時に、各ビーム出力信号Ｓ_θm(t) をＬビンに周波数分
析し、実数部と虚数部の２乗和をとり、各ビーム方向に
それぞれＬ個の周波数出力信号Ｓ_fml(t) （ｍ＝１〜
Ｍ，ｌ＝１〜Ｌ）として積分処理部３に送る。観測体２
１及び目標体２２が位置する自然環境では背景雑音が高
いため、ビーム出力信号Ｓ_θm(t) 及び周波数出力信号
Ｓ_fml(t) の中らか目標信号を取り出すためには積分を
行って、雑音のレベルを下げる必要がある。そこで、積
分処理部３では、設定情報入力端子４から入力された積
分時間等の情報に従い、入力された信号Ｓ_θm(t) ，Ｓ
_fml(t) を積分して、積分後の信号Ｓ’_θm(t) ，Ｓ’
_fml(t) を方位／周波数情報算出部５に送る。

【００２１】以下、積分処理部３での積分方法を説明す
る。なお、ここでは、ビーム出力信号Ｓ_θm(t) と周波
数出力信号Ｓ_fml(t) を区別することなく一般的な表記
Ｓ(t) によって説明する。ここで、積分前の信号レベル
のピークの位置が積分後の信号レベルのピークの位置に
近似的に等しいとみなす。時刻ｔ＝ｔ_jにおける積分前
の信号レベル（以下、瞬時信号レベルと呼ぶ）をＳ（ｔ
_j）、積分前の信号レベル真値（以下、瞬時信号レベル
真値と呼ぶ）をＳ₀（ｔ_j）、積分前の信号レベルの誤
差（以下、瞬時信号レベル誤差と呼ぶ）をｍ_s（ｔ_j）
とする。瞬時信号レベルは、瞬時信号レベル真値に瞬時
信号レベル誤差が加算されたものとなり、次の（１９）
式で表される。Ｓ（ｔ_j）＝Ｓ₀（ｔ_j）＋ｍ_s（ｔ_j） …（１９）指数積分は、積分後の信号レベル（以下、平滑化信号レ
ベルと呼ぶ）をＳ’（ｔ_j）とすると、（２０）式に示
すように、瞬時信号レベルＳ（ｔ_j）と、その直前時刻
までの平滑化信号レベルＳ’（ｔ_j-1）とを、重みづけ
加算して、その時刻の平滑化信号レベルＳ’（ｔ_j）を
得る方法である。ここで、αは、積分時間をＴ、積分処
理部への入力周期をΔｔとすると、（２１）式で示され
る。

【００２２】

【数９】

【００２３】なお、平滑化信号レベルＳ’（ｔ_j）は、
積分後のビーム出力信号Ｓ’_θm(t) あるいは、積分後
の周波数出力信号Ｓ’_fml(t) のいずれかに相当するも
のである。

【００２４】方位／周波数情報算出部５は入力信号Ｓ’
_θm(t) から信号を探知し、探知したビームとその近傍
のビームの信号レベルから補間により方位を算出し、目
標体２２の方位角の観測値β(t) とし、入力信号Ｓ’
_fml(t) から信号を探知し、探知したビンとその近傍の
ビンの信号レベルから補間により周波数を算出し、目標
音のドップラ周波数の観測値γ_k(t) （ｋ＝１〜Ｋ）と
して状態量推定部１２に送る。ところで、ビーム出力信
号及び周波数出力信号は（１９）式あるいは（２０）式
に示したように、信号レベルの真値に誤差が加わったも
のとなっている。よって、誤差を持つ信号レベルから算
出した方位角またはドップラ周波数の観測値は、それぞ
れの真値に対して誤差を持つ。

【００２５】そこで、さらに方位／周波数情報算出部５
は、方位角の観測値の誤差分散σ_β ²(t) 及びドップラ
周波数の観測値の誤差分散（σ_γk）²(t) を推定す
る。誤差分散の推定方法として、例えば、方位角の観測
値β(t) を平滑化し、ｔ＝ｔ_jにおける方位角の観測値
の平滑値β’（ｔ_j）を、誤差分散を求めるための基準
値として用い、（２２）式に示すように、所定時間内の
サンプル数Ｎの方位角の観測値β（ｔ_j-n）（ｎ＝０…
Ｎ−１）の前記基準値β’（ｔ_j）からの差の二乗平均
を適用することができる。同様に、第ｋドップラ周波数
の観測値γ_k(t) を平均化し、ｔ＝ｔ_jにおける第ｋド
ップラ周波数の観測値の平滑値γ’_k（ｔ_j）を、誤差
分散を求めるための基準値として用い、（２３）式に示
すように、所定時間内のサンプル数Ｎの第ｋドップラ周
波数の観測値γ_k（ｔ_j-n）（ｎ＝０…Ｎ−１）の前記
基準値γ’_k（ｔ_j）からの差の二乗平均を適用するこ
とができる。

【００２６】

【数１０】

【００２７】なお、所定時間（ｔ_j-(N-1)〜ｔ_j）内で
観測値の誤差は定常とみなすことができると仮定してお
り、また観測値を平滑化するときの積分時間は、積分処
理部３で行われる積分処理時間に比べて十分長くとる必
要がある。

【００２８】自船情報入力部１１は、該当観測体２１の
位置ベクトルｘ₁(t) 、速度ベクトルｘ₁′(t) 及び加
速度ベクトルｘ″₁(t) を得て、状態量推定部１２に与
える。状態量推定部１２では、方位／周波数情報算出部
５から入力された目標体２２の方位角の観測値β(t) と
その誤差分散σ_β ²(t) 及び目標音のドップラ周波数の
観測値γ_k(t) とその誤差分散（σ_γk）²（ｋ＝１〜
Ｋ）と、自船情報入力部１１から入力された、該当観測
体２１の位置ベクトルｘ₁(t) 、速度ベクトルｘ₁′
(t) 及び加速度ベクトルｘ₁″(t) とから、（８）式の
評価関数Ｌ（Ｘ₀）を最小とする初期状態量Ｘ₀を求
め、さらに目標体２２の状態量ｘ₂(t) の推定量を求
め、その推定結果を出力端子１０に与える。

【００２９】ここで、評価関数Ｌ（Ｘ₀）を最小とする
初期状態量Ｘ₀の推定には、例えばガウス・ニュートン
法が用いられる。ガウス・ニュートン法は、探索初期値
Ｘ₀ ⁽⁰⁾を出発点として初期状態量を修正して行き、Ｌ
（Ｘ₀）を最小とする初期状態量Ｘ₀を推定する。ま
た、探索初期値Ｘ₀ ⁽⁰⁾の算出には、例えば疑似線形法
が用いられる。疑似線形法は観測誤差が微小であると仮
定して、前記（１３）式及び（１４）式を初期状態量に
関して線形近似し、観測値時系列から初期状態量を推定
する方法である（前記文献を参照）。前記（８）式で表
される評価関数Ｌ（Ｘ₀）を形成するには、方位角及び
ドップラ周波数の観測値だけでなく、方位角及びドップ
ラ周波数の観測値についての誤差共分散行列Σ_β及びΣ
_γkが必要となる。ここで、前述のように積分処理はビ
ーム出力信号及び周波数出力信号に対して行っており、
方位角の観測値あるいはドップラ周波数の観測値を直接
積分するものではない。しかし、ビーム出力信号あるい
は周波数出力信号に対して行う積分は、等価的に、方位
角の観測値あるいはドップラ周波数の観測値に対して積
分を行っていると考えることができる。

【００３０】図４は積分された観測値のサンプル周期Δ
ｔと積分時間Ｔとの関係を示した図である。図４の
（Ａ）に示すように、積分時間Ｔがサンプル周期Δｔよ
り短いならば、各サンプル時刻ｔ₀，ｔ₁，…ｔ_nの積
分後の観測値は完全に切り分けられるので、サンプル時
刻が異なる観測誤差間は無相関となる。また図４の
（Ｂ）に示すように、積分時間Ｔがサンプル周期Δｔよ
り長いならば、各サンプル時刻ｔ₀，ｔ₁，…ｔ_nの積
分後の観測値は重複し、サンプル時刻が異なる観測誤差
間には相関性がある。従来の目標運動解析方法において
は、積分時間は観測値のサンプル周期にくらべて短く、
サンプル時刻が異なる観測誤差間は無相関であるか、あ
るいは積分時間は観測値のサンプル周期に比べて長い場
合でも、サンプル時刻が異なる観測誤差間の相関を無視
して無相関であるとしていた。このため、方位角及びド
ップラ周波数の観測値についての誤差共分散行列Σ_β及
びΣ_γkとして、対角行列を想定していた。すなわち、
下記の（２４）式及び（２５）式に示すように、非対角
要素（共分散要素）が０で、対角要素（分散要素）だけ
に各サンプル時刻の観測誤差の分散σ_β ²(t) ，（σ
_γk）²(t) を代入して共分散行列を求めていた。

【００３１】

【数１１】

【００３２】

【発明が解決しようとする課題】一般に、目標信号の信
号レベルが低い場合には、積分時間を長くしなければ信
号を探知できない。そして上述したように、積分時間Ｔ
がサンプル周期Δｔより長くなると、サンプル時刻が異
なる観測誤差は互に相関を有する。しかしながら、従来
の目標運動解析方法に従い、このような相関を有する観
測誤差を無相関とみなした（２４）式及び（２５）式に
示す誤差共分散行列Σ_β及びΣ_γkを適用すると、目標
運動の解析結果は最小分散推定にはならず、最良な推定
結果が得られない。目標の状態量を推定する手段とし
て、カルマンフィルタを用いる場合、カルマンフィルタ
では、異なる時刻間の観測誤差は互に無相関であること
を前提としているため、観測誤差が互に相関を有する場
合には、前述のように最良な推定結果が得られないとい
う問題があった。

【００３３】一方、目標の状態量を推定する手段とし
て、最尤推定法を用いる場合、（８）式における誤差共
分散行列Σ_β及びΣ_γkとして非対角行列を用いれば、
最良な推定結果を得ることは可能である。（８）式から
明らかなように、誤差共分散行列Σ_β及びΣ_γkの逆行
列Σ_β ^-1及びΣ_γk ^-1が（８）式には適用されている。
誤差共分散行列Σ_β及びΣ_γkが対角行列であれば、対
角要素の逆数を対角要素とすることで、逆行列Σ_β ^-1及
びΣ_γk ^-1を容易に得ることができる。しかし、誤差共
分散行列Σ_β及びΣ_γkが非対角行列となる場合、その
逆行列Σ_β ^-1及びΣ_γk ^-1を求める計算を必要とすると
いう問題があった。

【００３４】

【課題を解決するための手段】本発明に係る目標運動解
析方法は、目標から放射される音響又は電磁波信号を観
測位置の移動可能な観測体に設けられた受波センサアレ
イにより逐次受信し、この受信信号から前記目標の方位
角及びドップラ周波数を逐次算出し、この逐次算出され
た方位角及びドップラ周波数の時系列観測値に基づき前
記目標の状態量を推定する目標運動解析方法において、
まず前記方位角及びドップラ周波数の時系列観測値の誤
差分散をそれぞれ推定し、次に、前記方位角及びドップ
ラ周波数の時系列観測値及びその誤差分散をそれぞれ白
色化する処理を行い、次に、前記方位角及びドップラ周
波数のそれぞれの白色化処理後の時系列観測値及びその
誤差分散を用いて前記目標の状態量を推定するものであ
る。その結果、従来よりも精度の高い推定結果が得られ
ると共に、誤差共分散行列を対角要素のみにし、演算量
を削減できる。また観測値の積分時間を変更しても観測
値のサンプル周期を変更する必要がない。

【００３５】

【発明の実施の形態】本実施形態の説明に先立ち、観測
値が指数積分されている場合の、白色化処理の構成方法
について説明する。なお、ここでは、方位角及びドップ
ラ周波数を区別することなく一般的な表記ｙ(t) によっ
て説明する。時刻ｔ＝ｔ_jにおける積分後の観測値（以
下、平滑化観測値と呼ぶ）をｙ（ｔ_j）、積分後の観測
値の真値（以下、平滑化観測真値と呼ぶ）をｙ
₀（ｔ_j）、積分後の観測誤差（以下、平滑化観測誤差
と呼ぶ）をｎ_y（ｔ_j）とする。先に述べたように積分
処理部で行われる積分は、等価的に方位あるいはドップ
ラ周波数の観測値に対して積分を行っているとみなすこ
とができる。よって仮想的に、観測値の積分される前の
値（以下、瞬時観測値と呼ぶ）ｙ”（ｔ_j）、観測値の
真値の積分される前の値（以下、瞬時観測真値と呼ぶ）
ｙ”₀（ｔ_j）、観測誤差の積分される前の値（以下、
瞬時観測誤差と呼ぶ）ｎ”_y（ｔ_j）を考える。なお、
平滑化観測値ｙ（ｔ_j）は、方位角の観測値β（ｔ_j）
又は、ドップラ周波数の観測値γ_k（ｔ_j）に相当し、
平滑化観測誤差ｎ_y（ｔ_j）は、方位角の観測誤差ｎ_β
（ｔ_j）又は、第ｋドップラ周波数の観測誤差ｎ
_γk（ｔ_j）に相当するものである。

【００３６】想定した瞬時観測値は次の（２６）式で表
されるように瞬時観測真値と瞬時観測誤差の和となる。ｙ”（ｔ_j）＝ｙ”₀（ｔ_j）＋ｎ”_y（ｔ_j） …（２６）平滑化観測値ｙ（ｔ_j）は瞬時観測値を指数積分したも
のであるとして、次の（２７）式に示すように表され
る。

【００３７】

【数１２】

【００３８】このように積分処理された平滑化観測値に
たいして、白色化処理は次の（２８）式に示すように、
その時刻の平滑化観測値ｙ（ｔ_j）と、その直前時刻の
平滑化観測値ｙ（ｔ_j-1）とを、それぞれ重みづけ加算
する処理である。

【００３９】

【数１３】

【００４０】すなわち、平滑化観測値を２時刻分蓄えて
おき、これらから、瞬時観測値に変換する処理である。
よって、（２８）式で得られる、誤差が白色化された観
測値は前記（２６）式で示したように、瞬時観測真値に
瞬時観測誤差が加算されたものになっている。次に、指
数積分された観測値の誤差分散から、白色化後の誤差分
散を算出する原理を説明する。前記（２７）式から、平
滑化観測誤差ｎ_y（ｔ_j）は次の（２９）式で表され
る。

【００４１】

【数１４】

【００４２】ここで、瞬時観測誤差ｎ”_y（ｔ_j）の分
散を（σ_n"）²（ｔ_j）とすると、平滑化観測誤差ｎ_y
（ｔ_j）の分散σ_n ²（ｔ_j）は次の（３０）式で表さ
れる。

【００４３】

【数１５】

【００４４】よって、（３０）式から、白色化後の誤差
分散（σ_n"）²（ｔ_j）は次の（３１）式で表される。

【００４５】

【数１６】

【００４６】なお、以下の実施形態１及び実施形態２で
は、観測項目が方位角とドップラ周波数である場合を示
すが、観測項目が方位角だけである場合にも本発明を適
用することできる。また、以下の実施形態１及び実施形
態２では、状態量として、目標体２２の相対位置、相対
速度及び固有周波数の逆数を用いる場合を示すが、前記
（１３）式、（１４）式あるいは、後述する（４７）
式、（４８）式でモデル化されるように状態量をとれ
ば、固定座標系に本発明を適用することができる。

【００４７】実施形態１．実施形態１では、バッチ処理
を行うことによって目標の状態量を推定する方法とし
て、最尤推定法を用いた場合について説明する。図１は
本発明の実施形態１に係る目標運動解析装置の構成を示
す図である。図１において、１は受波器センサアレイか
らの出力信号が入力される信号入力端子、２は周波数分
析およびビームフォーミングを行う周波数分析／ビーム
フォーミング処理部、３は入力信号の積分を行う積分処
理部、４はオペレータによる積分時間及び目標探知操作
情報が入力される設定情報入力端子、５は周波数分析／
ビームフォーミングされた信号から目標の方位と周波数
情報を得る方位／周波数情報算出部、６は観測誤差の白
色化を行う白色化処理部、７は誤差が白色化された観測
値とその誤差分散を蓄える観測情報蓄積部、８は最尤推
定の探索初期値を設定する探索初期値算出部、９は観測
値時系列及び観測誤差分散時系列を用いて目標の状態量
を推定する最尤推定部、１０は推定結果の出力端子、１
１は自船情報入力部である。

【００４８】図１の動作を説明する。受波器センサアレ
イから信号入力端子１を介して入力される信号に対し
て、周波数分析／ビームフォーミング処理部２及び積分
処理部３が行う信号処理は、従来技術で述べた図３の目
標運動解析装置における処理と同一である。方位／周波
数情報算出部５は、積分処理器３からの入力信号Ｓ’
_θm(t) ，Ｓ’_fml(t) から、目標体２２の方位角の観
測値β(t) と該方位角の観測値の誤差分散σ_β ²(t) 及
び目標音のドップラ周波数の観測値γ_k(t) と該ドップ
ラ周波数の観測値の誤差分散（σ_γk）²(t) （ｋ＝１
〜Ｋ）を、従来技術の目標運動解析装置と同様の方法で
算出し、白色化処理部６に送る。

【００４９】白色化処理部６は入力された、時刻ｔ＝ｔ
_jとｔ＝ｔ_j-1の目標体２２の方位角の平滑化観測値β
（ｔ_j），β（ｔ_j-1）及び、時刻ｔ＝ｔ_jとｔ＝ｔ
_j-1の目標体２２のドップラ周波数の平滑化観測値γ_k
（ｔ_j），γ_k（ｔ_j-1）と、設定情報入力端子４から
入力された積分時間等の設定情報とから、観測値の誤差
を白色化する処理を行い、目標体２２の誤差を白色化し
た方位角の観測値β”（ｔ_j）及び、誤差を白色化した
ドップラ周波数の観測値γ”_k（ｔ_j）を算出し、観測
情報蓄積部７に与える。白色化処理は、前記（２８）式
で説明した処理である。すなわち、時刻ｔ＝ｔ_jの目標
体２２の観測誤差を白色化した方位角の観測値β”（ｔ
_j）を、時刻ｔ＝ｔ_jとｔ＝ｔ_j-1の目標体２２の方位
角の平滑化観測値β（ｔ_j），β（ｔ_j- ₁）から次の
（３２）式で算出する。同様に、時刻ｔ＝ｔ_jの目標体
２２の観測誤差を白色化したドップラ周波数の観測値
γ”_k（ｔ_j）を、時刻ｔ＝ｔ_jとｔ＝ｔ_j-1の目標体
２２のドップラ周波数の平滑化観測値γ_k（ｔ_j），γ
_k（ｔ_j-1）から次の（３３）式で算出する。

【００５０】

【数１７】

【００５１】さらに、白色化処理部６は、入力された時
刻ｔ＝ｔ_jとｔ＝ｔ_j-1の目標体２２の方位角の平滑化
観測値の誤差分散σ_β ²（ｔ_j），σ_β ²（ｔ_j-1）と
時刻ｔ＝ｔ_jとｔ＝ｔ_j-1の目標体２２のドップラ周波
数の平滑化観測値の誤差分散（σ_γk）²（ｔ_j），
（σ_γk）²（ｔ_j-1）から、前述の（３１）式に従
い、誤差を白色化した後の、目標体２２の方位角の観測
値の誤差分散（σ_β”）²を（３４）式で、ドップラ周
波数の観測値の誤差分散（σ_γk"）²（ｔ_j）を（３
５）式で算出し、観測情報蓄積部７に与える。

【００５２】

【数１８】

【００５３】ただし、ｔ＝ｔ₁の誤差分散については、
ｔ＝ｔ₀での平滑化観測値の誤差分散が得られないの
で、ｔ＝ｔ₁以前の瞬時観測誤差の分散は一定であった
と仮定して求めている。

【００５４】観測情報蓄積部７は入力された誤差を白色
化した方位角の観測値β”（ｔ_j）と、誤差を白色化し
たドップラ周波数の観測値γ”_k（ｔ_j）及び、誤差を
白色化した方位角の観測値の誤差分散（σ_β”）²（ｔ
_j）と誤差を白色化したドップラ周波数の観測値の誤差
分散（σ_γk"）²（ｔ_j）を蓄え、ｔ_j＝ｔ₁〜ｔ
_n（現時刻）の時系列データとして最尤推定部９に送る
とともに、観測値を用いて探索初期値を設定する必要が
ある場合には、探索初期値算出部８にも送る。探索初期
値算出部８は最尤推定法によって初期状態量Ｘ₀を推定
する際の、探索初期値Ｘ₀ ⁽⁰⁾を算出して、最尤推定部
９に送る。探索初期値の算出方法として例えば、観測情
報蓄積部７から観測値時系列を得て、疑似線形法によっ
て近似的な初期状態量を推定し、これを探索初期値とす
る。疑似線形法については従来技術及び前記公知文献で
説明されているので、ここでの説明は省略する。自船情
報入力部１１は、当該観測体１の位置ベクトルｘ₁(t)
、速度ベクトルｘ₁′(t) 及び加速度ベクトルｘ₁″
(t) を得て、最尤推定部９に与える。

【００５５】最尤推定部９は、観測情報蓄積部７から与
えられた、誤差を白色化した後の目標体２２の方位角の
観測値β”（ｔ_j）とその誤差分散（σ_β”）
²（ｔ_j）の時系列（ｔ_j＝ｔ₁〜ｔ_n）及びドップラ
周波数の観測値γ”_k(ｔ_j）とその誤差分散
（σ_γk"）²（ｔ_j）の時系列（ｔ_j＝ｔ₁〜ｔ_n）を
用いて、探索初期値算出部８から与えられた探索初期値
Ｘ₀ ⁽⁰⁾を初期値として、上述した（８）式に示す評価
関数Ｌ（Ｘ₀）をたて、これが最小となる、初期時刻ｔ
＝ｔ₀における初期状態量Ｘ₀を求め、これを初期状態
量の推定値Ｘ^* ₀とし、得られた初期状態量の推定値Ｘ
^* ₀と、解析結果を得たい時刻ｔの時間差ｔ−ｔ₀と、
自船情報入力部１１からの加速度ベクトルｘ₁″(t) と
を、上記（９）式を適用して時刻ｔにおける状態量Ｘ
(t) を求める。そして、状態量Ｘ(t) に含まれている、
又は、その要素から簡単に求めた出力形式の解析結果を
解析結果出力端子１０に与える。

【００５６】例えば、出力形式が時刻ｔ＝ｔ_nにおける
観測体２１に対する目標体２２の相対位置ベクトルｘ
(t) であれば、状態量Ｘ（ｔ_n）の該当要素をそのまま
出力する。また、出力形式が時刻ｔ＝ｔ_nにおける目標
体２２の絶対位置ベクトルｘ₂(t) であれば、状態量Ｘ
（ｔ_n）に含まれている相対位置ベクトルｘ（ｔ_n）に
加速度ベクトルｘ₁″(t) から求めた観測体２１の絶対
位置ベクトルｘ₁（ｔ_n）を加算して出力する。

【００５７】以下、最尤推定部９による処理の具体例を
説明する。観測情報蓄積部７から与えられた観測値は、
誤差が白色化されているので、その誤差共分散行列の非
対角要素は０となる。よって、上記（８）式は下記の
（３６）式のように変形できる。

【００５８】

【数１９】

【００５９】よって、本実施形態の最尤推定部９では上
記（３６）式の評価関数Ｌ′（Ｘ₀）が最小となるＸ₀
を求めることになる。評価関数Ｌ′（Ｘ₀）の最小化に
際しては（３７）式を解くことになる。

【００６０】

【数２０】

【００６１】上記（３７）式の非線形方程式を解く方法
として、ガウス・ニュートン法等の反復修正により、
（３６）式を最小化する初期状態量Ｘ₀を探索すること
で解く方法がある。すなわち、ガウス・ニュートン法で
は、探索初期値Ｘ₀ ⁽⁰⁾を出発点とし、（３８）式に従
って初期状態量Ｘ₀を探す。（３８）式において、μは
縮小因子であって０＜μ≦１であり、ΔＸ₀は（３９）
式を満たすような方向に向いた修正ベクトルである。

【００６２】

【数２１】

【００６３】次に、修正ベクトルΔＸ₀の決定方法につ
いて説明する。なおここでは、（４０）式に示すよう
に、方位角の観測値及びドップラ周波数の観測値をまと
めた観測値ベクトルＺ_mによって説明する。また、各観
測値に対応する、前記（１３）式、（１４）式のモデル
で計算される値のベクトルＺ（Ｘ₀）は（４１）式で表
す。また、共分散行列の逆行列Ｚ^-1は（４２）式で表さ
れる。なお、これらは以下の説明の都合から上記原理の
説明と別な表現を用いただけであり、その本質に変わり
はない。

【００６４】

【数２２】

【００６５】修正ベクトルΔＸ₀は、（４５）式に示す
正規方程式を解くことによって求められる。（４５）式
における行列Ｂはヤコビアン行列であり、このヤコビア
ン行列Ｂは、（１３）式、（１４）式で表されるモデル
の、時刻ｔ＝ｔ₀における初期状態量Ｘ₀での偏微分で
あり、（４６）式で定義される。

【００６６】

【数２３】

【００６７】すなわち、（４５）式に従って修正ベクト
ルΔＸ₀を求め、この修正ベクトルΔＸ₀を用いて（３
８）式で初期状態量Ｘ₀を更新しながら、評価関数Ｌ′
（Ｘ₀）を最小とする状態量を探せば良い。

【００６８】以上説明したように、本実施形態１によれ
ば、指数積分された観測値に対して、予め白色化する処
理を行い、観測誤差を無相関にしているため、重みとし
て用いる観測誤差の共分散行列が、前記（４２）式から
（４４）式に示したように対角行列となり、その逆行列
が各時刻の分散の逆数をとるだけで求められ、非対角行
列を求めるのに比べ、演算量の削減ができる。

【００６９】実施形態２．本実施形態２では、オンライ
ン処理を行うことによって目標の状態量を推定する方法
として、カルマンフィルタを用いた場合で説明する。説
明に先立ち、実施形態２での目標体２２の状態量をカル
マンフィルタを用いて推定する原理について説明する。
状態量Ｘ(t) を前述の（４）式のように定めたとき、状
態量の遷移関係は（４７）式で表され、状態量と観測値
の関係は（４８）式で表される。

【００７０】

【数２４】

【００７１】ここで、ｎ（ｔ_j）は時刻ｔ＝ｔ_jにおけ
る方位角及びドップラ周波数の観測誤差であり、（４
９）式で表す。

【００７２】

【数２５】

【００７３】カルマンフィルタは、この観測誤差が平均
０のガウス性白色雑音であることを前提としている。す
なわち、ｎ（ｔ_j）及びｎ（ｔ_j）ｎ^T（ｔ_j）の期待
値は式（５０）及び（５１）となる。

【００７４】

【数２６】

【００７５】前記（４７）式、（４８）式における関数
ｇ_tj，ｈ_tjは非線形関数であるので、ｇ_tj，ｈ_tjをそれ
ぞれ推定値Ｘ^* _tj/tj及びＸ^* _tj/tj-1のまわりに線形
化した、拡張カルマンフィルタを用いる。Ｘ^* _tj/tj-1
は時刻ｔ＝ｔ_j-1までの観測値から時刻ｔ＝ｔ_jの状態
量Ｘ（ｔ_j）を予測した値、Ｘ^* _tj/tjは時刻ｔ＝ｔ_j
までの観測値からのＸ（ｔ_j）の推定状態量である。こ
こで、（５２）式、（５３）式で表されるＧ_tj及びＨ_tj
を定義する。

【００７６】

【数２７】

【００７７】拡張カルマンフィルタのアルゴリズムは
（５４）式から（５８）式で与えられる。

【００７８】

【数２８】

【００７９】ここで、Ｘ^* _tj/tjが、時刻ｔ＝ｔ_jにお
ける状態量Ｘ（ｔ_j）の推定結果となる。

【００８０】図５は本発明の実施形態２に係る目標運動
解析装置の構成を示す図である。図５において、１は受
波器センサアレイからの信号が入力される信号入力端
子、２は周波数分析およびビームフォーミングを行う周
波数分析／ビームフォーミング処理部、３は入力信号の
積分を行う積分処理部、４はオペレータによる積分時間
及び目標探知操作情報が入力される設定情報入力端子、
５は周波数分析／ビームフォーミングされた信号から目
標の方位と周波数情報を得る方位／周波数情報算出部、
６は観測誤差の白色化を行う白色化処理部、１０は推定
結果の出力端子、１１は自船情報入力部、１３は状態量
を推定する初期値を設定する初期値設定端子、１４は状
態量を推定するカルマンフィルタ部である。なお、実施
形態１と同じ動作をするものについては、同一の番号を
付与している。

【００８１】図５の動作を説明する。図５の受波器セン
サアレイ１、周波数分析／ビームフォーミング処理部
２、積分処理部３、設定情報入力端子４、方位／周波数
情報算出部５、白色化処理部６は、図１の場合と同様の
動作を行うので説明を省略する。初期値設定端子１３は
時刻ｔ＝ｔ₁のとき状態量の初期値Ｘ^* _t0/t0及び推定
誤差共分散行列の初期値Ｐ_t0/t0を設定する。カルマン
フィルタ部１４は、白色化処理部６から与えられた、時
刻ｔ＝ｔ_jにおける、誤差を白色化した後の目標体２２
の方位角の観測値β”（ｔ_j）とその誤差分散
（σ_β”）²（ｔ_j）及びドップラ周波数の観測値γ”
_k（ｔ_j）とその誤差分散（σ_rk”）²（ｔ_j）とか
ら、観測値ベクトルＺ（ｔ_j）を（５９）式に示すよう
に、観測誤差共分散行列Ｒ_tjを（６０）式に示すように
生成し、この観測値ベクトルＺ（ｔ_j）と、観測誤差共
分散行列Ｒ_tj及び時刻ｔ＝ｔ_j-1における推定結果Ｘ^*
_tj-1/tj-1（ただし、ｔ＝ｔ₁のときは初期値設定端子
１３から設定されたＸ^* _t0/t0）と推定誤差共分散行列
Ｐ_tj-1/tj-1（だたし、ｔ＝ｔ₁のときは初期値設定端
子１３から設定されたＰ_t0/t0）から、上述の（５４）
式から（５８）式のアルゴリズムに基づいて時刻ｔ＝ｔ
_jにおける状態量Ｘ（ｔ_j）を推定するとともに、次の
時刻での推定のために、推定結果Ｘ^* _tj/tjと推定誤差
共分散行列Ｐ_tj/tjを蓄積する。

【００８２】

【数２９】

【００８３】そして、推定状態量Ｘ^* _tj/tjに含まれて
いる、又は、その要素から簡単に求めた出力形式の推定
結果を推定結果出力端子１０に与える。出力形式につい
ては、実施形態１と同様である。以上のように、本実施
形態２によれば、積分された観測値に対して、予め白色
化処理を施すことによって、観測誤差を無相関にしてい
るため、白色ガウス雑音を仮定しているカルマンフィル
タを用いて、オンラインで目標運動解析を行った場合、
最良な推定結果が得られる。

【００８４】なお、観測誤差が互に相関性を持つ場合
に、本発明（実施形態１又は実施形態２）を適用しなく
ても、図４の（Ｂ）に示したような観測値の重複がなく
なるように、観測値のサンプル周期を長くすればよい。
しかし、このようにすれば、積分時間に応じて状態量を
推定する処理周期（最尤推定部９あるいはカルマンフィ
ルタ部１４の処理周期）を変える必要が生じる。これに
対し、本発明（実施形態１又は実施形態２）を適用すれ
ば積分時間にかかわらず、状態量を推定する処理周期を
一定にすることができる。

【００８５】

【発明の効果】以上にように本発明によれば、目標から
放射される音響又は電磁波信号を観測位置の移動可能な
観測体に設けられた受波センサアレイにより逐次受信
し、この受信信号から前記目標の方位角及びドップラ周
波数を逐次算出し、この逐次算出された方位角及びドッ
プラ周波数の時系列観測値に基づき前記目標の状態量を
推定する目標運動解析方法において、まず前記方位角及
びドップラ周波数の時系列観測値の誤差分散をそれぞれ
推定し、次に、前記方位角及びドップラ周波数の時系列
観測値及びその誤差分散をそれぞれ白色化する処理を行
い、次に、前記方位角及びドップラ周波数のそれぞれの
白色化処理後の時系列観測値及びその誤差分散を用いて
前記目標の状態量を推定するようにしたので、その結
果、従来よりも精度の高い推定結果が得られると共に、
誤差共分散行列を対角要素のみにし、演算量を削減でき
る。また観測値の積分時間を変更しても観測値のサンプ
ル周期を変更する必要がないという効果が得られる。

【図面の簡単な説明】

【図１】本発明の実施形態１に係る目標運動解析装置の
構成を示す図である。

【図２】目標運動解析方法が適用される観測系及び運動
系の説明図である。

【図３】従来の目標運動解析装置の構成を示す図であ
る。

【図４】観測値のサンプル周期と積分時間との関係を示
す図である。

【図５】本発明の実施形態２に係る目標運動解析装置の
構成を示す図である。

【符号の説明】

１信号入力端子２周波数分析／ビームフォーミング処理部３積分処理器４設定情報入力端子５方位／周波数情報算出部６白色化処理部７観測情報蓄積部８探索初期値算出部９最尤推定部１０推定結果出力端子１１自船情報入力部１３初期値設定端子１４カルマンフィルタ部

───────────────────────────────────────────────────── フロントページの続き (51)Int.Cl.⁶ 識別記号庁内整理番号ＦＩ技術表示箇所Ｇ０１Ｓ 17/02 Ｇ０１Ｓ 17/02 Ｚ

Claims

【特許請求の範囲】

【請求項１】目標から放射される音響又は電磁波信号
を観測位置の移動可能な観測体に設けられた受波センサ
アレイにより逐次受信し、この受信信号から前記目標の
方位角及びドップラ周波数を逐次算出し、この逐次算出
された方位角及びドップラ周波数の時系列観測値に基づ
き前記目標の状態量を推定する目標運動解析方法におい
て、まず前記方位角及びドップラ周波数の時系列観測値の誤
差分散をそれぞれ推定し、次に、前記方位角及びドップラ周波数の時系列観測値及
びその誤差分散をそれぞれ白色化する処理を行い、次に、前記方位角及びドップラ周波数のそれぞれの白色
化処理後の時系列観測値及びその誤差分散を用いて前記
目標の状態量を推定することを特徴とする目標運動解析
方法。
【請求項２】前記時系列観測値及びその誤差分散に対
する白色化処理において、前記時系列観測値の積分方法として、現時刻の瞬時観測
値と、一入力周期前までの積分値にそれぞれ固定の重み
係数α及び（１−α）を乗算して加算することにより求
める指数積分が施されている場合に、誤差を白色化した
観測値ｙ”（ｔ_j）を、下記式で算出することを特徴と
する請求項１記載の目標運動解析方法。【数１】
【請求項３】前記目標の状態量を推定する方法におい
て、前記目標の方位角度及びドップラ周波数の時系列観測値
を蓄積した時系列と、状態量から算出される観測値時系
列の各時刻に対応する目標の方位角及びドップラ周波数
の推定値の時系列との差の重み付き２乗和を最小とする
ように、前記蓄積した観測値を一括処理して推定を行う
ことを特徴とする請求項１記載の目標運動解析方法。
【請求項４】前記目標の状態量を推定する方法におい
て、各入力周期の前記目標の方位角及びドップラ周波数の観
測値が入力される毎に、入力された観測値と、一入力周
期前の推定状態量から予測される現時刻の目標の方位角
及びドップラ周波数の予測値との差から求めた予測誤差
が最小となるように、一入力周期前の推定状態量を修正
することにより状態量の推定を行うことを特徴とする請
求項１記載の目標運動解析方法。