JP6625295B1 - 目標監視装置および目標監視システム - Google Patents

Abstract

目標監視装置は、目標から複数の通信経路を経て到来した電波信号をそれぞれ受信する複数の受信アンテナで得られた複数の受信信号に基づいて当該目標の状態を推定する。目標監視装置は、複数の受信信号に基づいて複数の計測値を逐次算出する計測部（６１）、複数の計測値に含まれる到来時間差または到来周波数差を監視して目標の動きの有無を検出する状態検出部（６２Ａ）と、測位演算および追尾演算を実行する演算部（６３Ａ，６４Ａ，６５）とを備える。演算部（６３Ａ，６４Ａ，６５）は、目標の動きが無いことが検出されたときは、測位演算を実行して目標の推定位置を示す状態推定値を算出し、目標の静止状態から移動状態への変化が検出されたときは、追尾演算を実行して当該目標の推定位置および移動状態の双方を示す状態推定値を逐次算出する。

Description

本発明は、電波発信源から複数の衛星を経由して受信された電波に基づき、当該電波発信源の位置などの状態を推定する技術に関するものである。

複数の衛星を用いて電波発信源の位置を推定する技術が知られている。たとえば、下記の非特許文献１には、電波発信源から複数の衛星を経由して受信された電波間の到来時間差（Ｔｉｍｅ−Ｄｉｆｆｅｒｅｎｃｅｓ−Ｏｆ−Ａｒｒｉｖａｌ，ＴＤＯＡｓ）に基づく測位方法と、当該受信された電波間のドップラー周波数差を示す到来周波数差（Ｆｒｅｑｕｅｎｃｙ−Ｄｉｆｆｅｒｅｎｃｅｓ−Ｏｆ−Ａｒｒｉｖａｌ，ＦＤＯＡｓ）および到来時間差（ＴＤＯＡｓ）に基づく測位方法とが開示されている。

このような測位方法を採用する目標監視システムでは、電波発信源である目標から発信された電波は、複数の衛星にそれぞれ搭載されたトランスポンダ（Ｔｒａｎｓｐｏｎｄｅｒｓ）を経由した後に、地上の受信局における複数の受信アンテナで受信される。測位演算器は、当該複数の受信アンテナで得られた受信電波間の到来時間差（ＴＤＯＡｓ）を測定し、あるいは、当該受信電波間の到来周波数差（ＦＤＯＡｓ）および到来時間差（ＴＤＯＡｓ）を計測する。そして、測位演算器は、計測された到来時間差、あるいは、計測された到来周波数差および到来時間差に基づき、電波発信源が地球の地表面に存在するとの仮定の下で、電波発信源の位置を推定することができる。当該複数の衛星の軌道位置は互いに異なるので、当該受信電波間で伝搬時間の差が生じる。これにより、当該受信電波間の到来時間差（ＴＤＯＡｓ）の計測が可能となる。また、当該複数の衛星は、たとえ静止衛星であったとしても、地表面に対して若干運動しているため、当該受信電波間でドップラー周波数の差が生じる。これにより、当該受信電波間の到来周波数差（ＦＤＯＡｓ）の計測が可能となる。

K. C. Ho and Y. T. Chan, "Geolocation of a known altitude object from TDOA and FDOA measurements," IEEE Transactions on Aerospace and Electronic Systems, Vol. 33, Issue 3, pp. 770-783, 1997.

移動自在な電波発信源（たとえば、船舶）の状態を監視しようとする場合に、従来の目標監視システムでは、目標となる電波発信源が静止状態（船舶の場合は、当該船舶が停泊している状態）にあるときから監視が開始されることがある。このとき、目標が静止状態にあるときは当該目標の瞬時位置を推定する測位演算が実行され、当該目標が移動を開始した後は当該目標を追尾する追尾（ｔｒａｃｋｉｎｇ）演算が実行される。移動目標の追尾時には、移動目標の、時々刻々と変化する位置を推定することのみならず、時刻毎の移動速度を推定することが求められる。目標の状態が静止状態から移動状態へ変化すると、測位演算を追尾演算に切り替える必要がある。しかしながら、従来の目標監視システムでは、その変化に応じて測位演算を追尾演算へ切り替えるべきタイミングを正確に決定することができず、移動目標の追尾に失敗するおそれがあった。

上記に鑑みて本発明の目的は、電波発信源である目標の状態変化に応じて測位演算を追尾演算に切り替えるべきタイミングを正確に決定し、当該目標に対する高い追尾精度を確保することができる目標監視装置および目標監視システムを提供する点にある。

本発明の一態様による目標監視装置は、電波発信源である目標から、複数の衛星を含む３つ以上の通信経路を経て到来した電波信号をそれぞれ受信する複数の受信アンテナで得られた複数の受信信号に基づいて当該目標の状態を推定する目標監視装置であって、前記複数の受信信号に基づき、前記受信信号間の到来時間差、または前記受信信号間の到来時間差および到来周波数差を示す複数の計測値を逐次算出する計測部と、前記複数の計測値のうちの少なくとも１つの計測値の時間変動を監視して当該目標の動きの有無を検出する状態検出部と、前記状態検出部により前記目標の動きが無いことが検出されたときは、前記複数の計測値を用いた測位演算を実行して前記目標の推定位置を示す第１の状態推定値を算出し、前記状態検出部により前記目標の静止状態から移動状態への状態変化が検出されたときは、追尾演算を実行して当該目標の推定位置および移動状態の双方を示す第２の状態推定値を算出する演算部とを備え、前記状態検出部は、監視エリアを定める座標情報、前記測位演算もしくは前記追尾演算により算出された当該目標の推定位置、または誤差楕円のいずれか１つと、前記複数の衛星それぞれの軌道計算値とに基づいて数値範囲を設定し、前記状態検出部は、前記少なくとも１つの計測値にフィルタ処理を施して平滑値を算出し、当該平滑値が、前記数値範囲内から当該数値範囲外へ変動したときに、前記目標の静止状態から移動状態への状態変化が生じたと判定することを特徴とする。

本発明によれば、目標の静止状態から移動状態への状態変化に応じて測位演算を追尾演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、当該目標に対する高い追尾精度を確保することができる。

本発明に係る実施の形態１である目標監視システムの構成を概略的に示す図である。 実施の形態１における受信器の構成例を示すブロック図である。 実施の形態１における監視処理部の概略構成を示すブロック図である。 実施の形態１の状態検出部の構成例を概略的に示すブロック図である。 図５Ａ，図５Ｂは、平滑値の時間変動の例を表すグラフである。 カルマンフィルタによる反復推定処理の手順の一例を概略的に示すフローチャートである。 実施の形態１における測位演算部、追尾演算部およびデータ移行部からなる演算部の構成を概略的に示すブロック図である。 監視処理部の機能を実現するハードウェア構成例を示すブロック図である。 本発明に係る実施の形態２における測位演算部、追尾演算部およびデータ移行部からなる演算部の構成を概略的に示すブロック図である。 本発明に係る実施の形態３における監視処理部の概略構成を示すブロック図である。 実施の形態３における測位演算部、追尾演算部およびデータ移行部からなる演算部の構成を概略的に示すブロック図である。 本発明に係る実施の形態４における測位演算部、追尾演算部およびデータ移行部からなる演算部の構成を概略的に示すブロック図である。 本発明に係る実施の形態５における監視処理部の概略構成を示すブロック図である。 実施の形態５における状態検出部の概略構成を示すブロック図である。 本発明に係る実施の形態６における監視処理部の概略構成を示すブロック図である。 実施の形態６における状態検出部の概略構成を示すブロック図である。

以下、図面を参照しつつ、本発明に係る種々の実施の形態について詳細に説明する。なお、図面全体において同一符号を付された構成要素は、同一構成及び同一機能を有するものとする。

実施の形態１．
図１は、本発明に係る実施の形態１である目標監視システム１の構成を概略的に示す図である。図１に示されるように目標監視システム１は、監視エリアＳＡ内の電波発信源である目標Ｔｇｔから複数の通信経路を経て到来した電波信号をそれぞれ受信して受信信号Ｓ１，Ｓ２，Ｓ３を出力する受信アンテナ（受信局）Ｒｘ１，Ｒｘ２，Ｒｘ３と、受信信号Ｓ１，Ｓ２，Ｓ３に基づいて目標Ｔｇｔの状態を推定する目標監視装置２とを備えて構成されている。本実施の形態では、受信アンテナＲｘ１，Ｒｘ２，Ｒｘ３と目標Ｔｇｔとの間の３本の通信経路にそれぞれ３基の人工衛星Ｓｔ１，Ｓｔ２，Ｓｔ３（以下、単に「衛星Ｓｔ１，Ｓｔ２，Ｓｔ３」という。）が存在する。衛星Ｓｔ１，Ｓｔ２，Ｓｔ３に搭載されたトランスポンダ（中継器）は、目標Ｔｇｔから発信された電波信号を受信し、当該受信された電波信号を増幅してダウンリンク信号を生成し、これらダウンリンク信号を目標監視システム１の受信アンテナＲｘ１，Ｒｘ２，Ｒｘ３にそれぞれ送信する。衛星Ｓｔ１，Ｓｔ２，Ｓｔ２としては、静止軌道に投入された静止衛星が想定されている。なお、３基の衛星Ｓｔ１，Ｓｔ２，Ｓｔ３に限定されずに２基またはＮ基の衛星（Ｎは４以上の整数）からそれぞれダウンリンク信号を受信可能なように本実施の形態の目標監視システム１の構成が変更されてもよい。

目標Ｔｇｔは、移動自在な電波発信源（船舶、航空機または車両などの移動体）であり、特定の周波数帯域の電波信号を衛星に向けて発信する機能を有する。しかしながら、目標Ｔｇｔから発信された電波信号が、地上の無線通信網または衛星通信網の送信電波と干渉して通信障害を引き起こすことがある。このような場合に、干渉波となる電波信号を発信する目標Ｔｇｔの位置を特定したいというニーズがある。また、目標Ｔｇｔから発信された捜索救助用のビーコン信号に基づいて目標Ｔｇｔの位置を特定したいというニーズもある。本実施の形態の目標監視装置２は、受信信号Ｓ１〜Ｓ３に基づいて、目標Ｔｇｔの動きの有無、すなわち目標Ｔｇｔの状態が静止状態または移動状態のいずれであるかを検出することができる。目標Ｔｇｔが静止状態にあるとき、目標監視装置２は、目標Ｔｇｔの推定位置（静止位置）を示す状態推定値を逐次算出することができ、目標Ｔｇｔが移動状態にあるときは、目標監視装置２は、目標Ｔｇｔの推定位置（追尾位置）および移動状態（たとえば速度）を示す状態推定値を逐次算出することができる。

目標監視装置２は、図１に示されるように、高周波帯域の受信信号Ｓ１，Ｓ２，Ｓ３をディジタル受信信号Ｄ１，Ｄ２，Ｄ３に変換する受信器１０と、ディジタル受信信号Ｄ１，Ｄ２，Ｄ３に基づいて目標Ｔｇｔの状態を示す状態推定値を逐次算出する監視処理部１１Ａと、当該算出された状態推定値に基づく画像情報を表示する表示器１２とを備える。

図２は、実施の形態１における受信器１０の構成例を示すブロック図である。図２に示される受信器１０は、周波数変換用の局部発振信号を出力する局部発振器２０と、高周波帯域の受信信号Ｓ１，Ｓ２，Ｓ３にそれぞれフィルタ処理を施すバンドパスフィルタ３１，４１，５１と、局部発振信号を用いてバンドパスフィルタ３１，４１，５１の高周波出力Ｆ１，Ｆ２，Ｆ３にそれぞれ周波数変換を施すダウンコンバータ３２，４２，５２と、Ａ／Ｄ変換器（ＡＤＣ）３３，４３，５３とを有している。ダウンコンバータ３２，４２，５２は、高周波出力Ｆ１，Ｆ２，Ｆ３を局部発振信号と混合して、より低い周波数帯域のアナログ信号Ｃ１，Ｃ２，Ｃ３を生成する周波数変換器である。ＡＤＣ３３，４３，５３は、アナログ信号Ｃ１，Ｃ２，Ｃ３をそれぞれディジタル受信信号Ｄ１，Ｄ２，Ｄ３に変換し、当該ディジタル受信信号Ｄ１，Ｄ２，Ｄ３を監視処理部１１Ａに出力する。ディジタル受信信号Ｄ１，Ｄ２，Ｄ３の各々は、振幅成分と位相成分とを含む複素信号である。

図３は、実施の形態１における監視処理部１１Ａの概略構成を示すブロック図である。図３に示される監視処理部１１Ａは、ディジタル受信信号Ｄ１，Ｄ２，Ｄ３に基づいて電波の到来時間差および到来周波数差の計測値を逐次算出する計測部６１と、これら計測値のいずれかに基づいて目標Ｔｇｔの動きの有無を検出する状態検出部６２Ａと、目標Ｔｇｔの動きが無いことが検出されたときは、当該計測値を用いた測位演算を実行して目標Ｔｇｔの推定位置を示す状態推定値を逐次算出する測位演算部６３Ａと、目標Ｔｇｔの動きが有ることが検出されたきは、当該計測値を用いた追尾演算を実行して目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）の双方を示す状態推定値を逐次算出する追尾演算部６４Ａと、測位演算部６３Ａおよび追尾演算部６４Ａの一方から他方へ演算結果のデータを移行させるデータ移行部６５と、状態検出部６２Ａでの処理に必要なデータ（衛星軌道情報ＯＤ，監視エリア情報ＡＤおよび目標揺動情報ＳＤ）を格納するデータ記憶部６７と、外部サーバ（図示せず。）から当該データを取得するデータ取得部６８とを有する。ここで、ｋは時刻ｔ_ｋを示す整数である。

また、監視処理部１１Ａは、測位演算部６３Ａおよび追尾演算部６４Ａで算出された状態推定値を基に画像情報ＤＤを生成する出力制御部６６を備える。出力制御部６６は、当該画像情報ＤＤを液晶ディスプレイまたは有機ＥＬディスプレイなどの表示器１２に出力する。画像情報ＤＤとしては、たとえば、目標Ｔｇｔの推定位置の座標値を表す英数字情報、その推定位置を視覚的に表す地図情報、目標Ｔｇｔの推定速度の値を示す英数字情報、および、その推定速度の遷移を視覚的に表すグラフが挙げられる。ユーザは、表示器１２に表示された画像情報ＤＤを視ることで目標Ｔｇｔの状態を把握することができる。

計測部６１は、ディジタル受信信号Ｄ１，Ｄ２，Ｄ３に基づき、受信信号Ｓ１，Ｓ２間の到来時間差ＴＤＯＡ_１２（ｋ）、受信信号Ｓ２，Ｓ３間の到来時間差ＴＤＯＡ_２３（ｋ）、受信信号Ｓ３，Ｓ１間の到来時間差ＴＤＯＡ_３１（ｋ）、受信信号Ｓ１，Ｓ２間の到来周波数差ＦＤＯＡ_１２（ｋ）、受信信号Ｓ２，Ｓ３間の到来周波数差ＦＤＯＡ_２３（ｋ）、および、受信信号Ｓ３，Ｓ１間の到来周波数差ＦＤＯＡ_３１（ｋ）の計測値を逐次算出する機能を有する。計測値を算出する具体的な方法としては、公知の方法が採用されればよく、特に限定されるものではない。たとえば、時間方向および周波数方向についてディジタル受信信号Ｄｉ，Ｄｊ間の２次元の相互相関を算出し、当該算出された相互相関に現れるピークを検出し、当該ピークに対応する到来時間差ＴＤＯＡ_ｉｊ（ｋ）および到来周波数ＦＤＯＡ_ｉｊ（ｋ）の組を得るという算出法を採用することができる。ここで、ｉ，ｊは、１〜３の範囲内の整数である（ｉ≠ｊ）。

また、計測部６１は、到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ），ＦＤＯＡ_２３（ｋ），ＦＤＯＡ_３１（ｋ）の計測値を、状態検出部６２Ａ、測位演算部６３Ａおよび追尾演算部６４Ａに供給する。

状態検出部６２Ａは、計測部６１から逐次供給される計測値のうちの少なくとも１つを監視対象とし、この監視対象の計測値の時間変動を常時監視して目標Ｔｇｔの動きの有無を検出することができ、その検出結果を示す検出信号ＥＳを測位演算部６３Ａ，追尾演算部６４Ａおよびデータ移行部６５に供給する。

図４は、実施の形態１の状態検出部６２Ａの構成例を概略的に示すブロック図である。最初に、到来時間差ＴＤＯＡ_ｉｊ（ｋ）を監視対象として選択した場合の状態検出部６２Ａの構成および動作について説明する。

図４に示される状態検出部６２Ａは、監視対象の計測値の時系列データにフィルタ処理を施して平滑値＜ＴＤＯＡ＞_ｋを算出する平滑値算出部８１と、閾値ＴＨ１，ＴＨ２の組を設定する閾値設定部８２Ａと、平滑値＜ＴＤＯＡ＞_ｋを、閾値ＴＨ１を下限としかつ閾値ＴＨ２を上限とする数値範囲Δ（ＴＨ１，ＴＨ２）と比較して目標Ｔｇｔの状態が静止状態または移動状態のいずれであるかを判定する状態判定部８３とを有する。ここで、ｋは時刻ｔ_ｋを示す整数である。状態判定部８３は、目標Ｔｇｔの状態が静止状態または移動状態のいずれであるかを示す判定結果ＥＳ１を判定出力部８９に与える。たとえば、Ｎ個の到来時間差ＴＤＯＡ_ｉｊ（ｋ−Ｎ＋１），ＴＤＯＡ_ｉｊ（ｋ−Ｎ＋２），…，ＴＤＯＡ_ｉｊ（ｋ−１），ＴＤＯＡ_ｉｊ（ｋ）の計測値にフィルタ処理（たとえば、平均化）を施すことで、時刻ｔ_ｋでの平滑値＜ＴＤＯＡ＞_ｋを算出することが可能である。到来時間差ＴＤＯＡ_ｉｊ（ｋ）の値には、計測誤差に起因するバラツキが存在する。平滑値＜ＴＤＯＡ＞_ｋを使用することにより、そのようなバラツキの影響を低減させることが可能となる。

状態判定部８３は、平滑値＜ＴＤＯＡ＞_ｋの時間変動を常時監視し、平滑値＜ＴＤＯＡ＞_ｋが一定期間Ｔ１の間、数値範囲Δ（ＴＨ１，ＴＨ２）内に継続して収まるときに、目標Ｔｇｔは静止状態にあると判定する。図５Ａは、平滑値＜ＴＤＯＡ＞_ｋの時間変動の一例を表すグラフである。図５Ａのグラフに示されるように、目標Ｔｇｔが静止状態にあるとき（ｔ＜ｔ_ｍのとき）、平滑値＜ＴＤＯＡ＞_ｋは、閾値ＴＨ１，ＴＨ２間の数値範囲内に収まるように変動する。

平滑値＜ＴＤＯＡ＞_ｋが数値範囲Δ（ＴＨ１，ＴＨ２）内から数値範囲Δ（ＴＨ１，ＴＨ２）外へ変動したとき、状態判定部８３は、目標Ｔｇｔの静止状態から移動状態への状態変化が生じたと判定する。図５Ａの例では、平滑値＜ＴＤＯＡ＞_ｋは、時刻ｔ_ｍで上限の閾値ＴＨ２を超えるように変動しているため、状態判定部８３は、時刻ｔ_ｍで目標Ｔｇｔの状態変化が生じたと判定することができる。

閾値ＴＨ１，ＴＨ２については、目標Ｔｇｔの静止状態から移動状態への状態変化による平滑値＜ＴＤＯＡ＞_ｋの変動を速やかに検知することができるように数値範囲Δ（ＴＨ１，ＴＨ２）の幅を極力小さくする一方で、その状態変化以外の要因による平滑値＜ＴＤＯＡ＞_ｋの変動を検知しないように数値範囲Δ（ＴＨ１，ＴＨ２）の幅を大きくすることが望ましい。その状態変化以外の要因としては、主に、衛星Ｓｔｉ，Ｓｔｊの動きが挙げられる。したがって、目標Ｔｇｔの状態変化による平滑値＜ＴＤＯＡ＞_ｋの変動と、衛星Ｓｔｉ，Ｓｔｊの動きによる平滑値＜ＴＤＯＡ＞_ｋの変動とを考慮して、適切な閾値ＴＨ１，ＴＨ２を設定することが望ましい。

閾値設定部８２Ａは、データ記憶部６７から供給された衛星軌道情報ＯＤおよび監視エリア情報ＡＤに基づき、次の状態判定式（１），（２）に示される閾値ＴＨ１，ＴＨ２の組を設定することができる。

ここで、ＴＤＯＡ_ｍｉｎは、ｉ番目の衛星Ｓｔｉおよびｊ番目の衛星Ｓｔｊ（ｉ≠ｊ）の動きを考慮して設定された最小値であり、ＴＤＯＡ_ｍａｘは、衛星Ｓｔｉ，Ｓｔｊの動きを考慮して設定された最大値である。また、σ_ＴＤＯＡは、到来時間差の計測誤差の標準偏差である。３σ_ＴＤＯＡは、到来時間差の計測誤差を考慮して標準偏差の３倍相当のマージンを与えるものである。

状態判定部８３は、状態判定式（１），（２）のいずれか一方が満たされたことを検知したときに、目標Ｔｇｔの静止状態から移動状態への状態変化が生じたと判定することができる。最小値ＴＤＯＡ_ｍｉｎおよび最大値ＴＤＯＡ_ｍａｘは、衛星軌道情報ＯＤおよび監視エリア情報ＡＤを用いて算出可能である。監視エリア情報ＡＤは、図１に示した監視エリアＳＡを定める座標情報である。衛星軌道情報ＯＤは、衛星Ｓｔｉ，Ｓｔｊの軌道計算に必要な軌道要素情報、または、軌道要素情報を用いた軌道計算アルゴリズムにより算出された軌道計算値であればよい。

たとえば、閾値設定部８２Ａは、あらかじめ、監視エリア情報ＡＤで定められた監視エリアＳＡ内の複数の地点Ψ_１，Ψ_２，…，Ψ_Ｌ（Ｌは正整数）に電波発信源が存在するとの仮定の下で、衛星Ｓｔｉ，Ｓｔｊの最新の軌道計算値を用いて地点Ψ_１，Ψ_２，…，Ψ_Ｌにそれぞれ対応する到来時間差ＴＤＯＡ（Ψ_１，ｋ），ＴＤＯＡ_ｋ（Ψ_２，ｋ），…，ＴＤＯＡ（Ψ_Ｌ，ｋ）を予測する（ｋは時刻ｔ_ｋを示す整数）。閾値設定部８２Ａは、予測された到来時間差ＴＤＯＡ（Ψ_１，ｋ）〜ＴＤＯＡ（Ψ_Ｌ，ｋ）のうちの最小値および最大値をそれぞれ最小値ＴＤＯＡ_ｍｉｎおよび最大値ＴＤＯＡ_ｍａｘとして決定することができる。衛星Ｓｔｉ，Ｓｔｊの軌道は２４時間周期を持つため、閾値設定部８２Ａは、最新の軌道計算値を用いて２４時間分の最小値ＴＤＯＡ_ｍｉｎ（ｋ）および最大値ＴＤＯＡ_ｍａｘ（ｋ）をあらかじめ計算すれば十分である。また、閾値設定部８２Ａは、予測された到来時間差ＴＤＯＡ（Ψ_１，ｋ）〜ＴＤＯＡ（Ψ_Ｌ，ｋ）の中から、平滑値＜ＴＤＯＡ＞_ｋに近い到来時間差を選択し、当該選択された到来時間差に基づいて最小値ＴＤＯＡ_ｍｉｎおよび最大値ＴＤＯＡ_ｍａｘを決定してもよい。

次に、目標Ｔｇｔの移動状態から静止状態への状態変化が生じた場合、平滑値＜ＴＤＯＡ＞_ｋの変動幅が一定範囲内に収まる。この場合、状態判定部８３は、平滑値＜ＴＤＯＡ＞_ｋが一定期間Ｔ１の間、数値範囲Δ（ＴＨ１，ＴＨ２）内に継続して収まることを検知したときに、目標Ｔｇｔの移動状態から静止状態への状態変化が生じたと判定する。一定期間Ｔ１の時間長は、パラメータ設定値である。図５Ｂは、平滑値＜ＴＤＯＡ＞_ｋの時間変動の他の例を表すグラフである。図５Ｂのグラフに示されるように、平滑値＜ＴＤＯＡ＞_ｋの変動幅が時刻ｔ_ｄ以後小さくなり、時刻ｔ_ｄ〜時刻ｔ_ｓの間、平滑値＜ＴＤＯＡ＞_ｋは、数値範囲Δ（ＴＨ１，ＴＨ２）内に継続して収まるように推移している。このとき、状態判定部８３は、時刻ｔ_ｓで目標Ｔｇｔの移動状態から静止状態への状態変化が生じたと判定することができる。

より具体的には、状態判定部８３は、一定期間Ｔ１の間、次の状態判定式（３），（４）の双方が満たされることを継続して検知したときに、目標Ｔｇｔの移動状態から静止状態への状態変化が生じたと判定することができる。

次に、図４を参照しつつ、到来周波数差ＦＤＯＡ_ｉｊ（ｋ）を監視対象として選択した場合の状態検出部６２Ａの構成および動作について説明する。

図４に示される状態検出部６２Ａは、監視対象である計測値の時系列データにフィルタ処理を施して平滑値＜ＦＤＯＡ＞_ｋを算出する平滑値算出部８５と、閾値ＴＨ３，ＴＨ４の組を設定する閾値設定部８６Ａと、平滑値＜ＦＤＯＡ＞_ｋを、閾値ＴＨ３を下限としかつ閾値ＴＨ４を上限とする数値範囲Δ（ＴＨ３，ＴＨ４）と比較して目標Ｔｇｔの状態が静止状態または移動状態のいずれであるかを判定する状態判定部８７とを有する。ここで、ｋは時刻ｔ_ｋを示す整数である。たとえば、Ｎ個の到来周波数差ＦＤＯＡ_ｉｊ（ｋ−Ｎ＋１），ＦＤＯＡ_ｉｊ（ｋ−Ｎ＋２），…，ＦＤＯＡ_ｉｊ（ｋ−１），ＦＤＯＡ_ｉｊ（ｋ）の計測値にフィルタ処理（たとえば、平均化）を施すことで、時刻ｔ_ｋでの平滑値＜ＦＤＯＡ＞_ｋを算出することが可能である。到来周波数差ＦＤＯＡ_ｉｊ（ｋ）の値には、計測誤差に起因するバラツキが存在する。平滑値＜ＦＤＯＡ＞_ｋを使用することにより、そのようなバラツキの影響を低減させることが可能となる。

状態判定部８７は、平滑値＜ＦＤＯＡ＞_ｋの時間変動を常時監視し、平滑値＜ＦＤＯＡ＞_ｋが一定期間Ｔ２の間、数値範囲Δ（ＴＨ３，ＴＨ４）内に継続して収まるときに、目標Ｔｇｔは静止状態にあると判定することができる。また、平滑値＜ＦＤＯＡ＞_ｋが数値範囲Δ（ＴＨ３，ＴＨ４）内から数値範囲Δ（ＴＨ３，ＴＨ４）外へ変動したとき、状態判定部８７は、目標Ｔｇｔの静止状態から移動状態への状態変化が生じたと判定することができる。

閾値ＴＨ３，ＴＨ４については、目標Ｔｇｔの静止状態から移動状態への状態変化による平滑値＜ＦＤＯＡ＞_ｋの変動を速やかに検知することができるように数値範囲Δ（ＴＨ３，ＴＨ４）の幅を極力小さくする一方で、その状態変化以外の要因による平滑値＜ＦＤＯＡ＞_ｋの変動を検知しないように数値範囲Δ（ＴＨ３，ＴＨ４）の幅を大きくすることが望ましい。その状態変化以外の要因としては、主に、衛星Ｓｔｉ，Ｓｔｊの動きと、目標Ｔｇｔの姿勢変化などの動揺によるドップラー効果とが挙げられる。たとえば、目標Ｔｇｔが船舶の場合、波または風の影響による船舶の動揺（たとえば、縦揺れおよび横揺れ）が生じ、これに起因するドップラー効果が発生することがある。したがって、目標Ｔｇｔの状態変化による平滑値＜ＦＤＯＡ＞_ｋの変動と、人工衛星Ｓｔｉ，Ｓｔｊの動きによる平滑値＜ＦＤＯＡ＞_ｋの変動と、目標Ｔｇｔの動揺による平滑値＜ＦＤＯＡ＞_ｋの変動とを考慮して、適切な閾値ＴＨ３，ＴＨ４を設定することが望ましい。

閾値設定部８６Ａは、データ記憶部６７から供給された衛星軌道情報ＯＤ，監視エリア情報ＡＤおよび目標揺動情報ＳＤに基づき、次の状態判定式（５），（６）に示される閾値ＴＨ３，ＴＨ４の組を設定することができる。

ここで、ＦＤＯＡ_ｍｉｎは、ｉ番目の衛星Ｓｔｉおよびｊ番目の衛星Ｓｔｊ（ｉ≠ｊ）の動きを考慮して設定された最小値であり、ＦＤＯＡ_ｍａｘは、衛星Ｓｔｉ，Ｓｔｊの動きを考慮して設定された最大値である。また、σ_ＦＤＯＡは、到来周波数差の計測誤差の標準偏差である。３σ_ＦＤＯＡは、到来周波数差の計測誤差を考慮して標準偏差の３倍相当のマージンを与えるものである。さらに、σ_{ｍｏｔｉｏｎ}は、目標揺動情報ＳＤに基づき、目標Ｔｇｔの動揺によるドップラー効果を考慮して設定された値であり、当該ドップラー効果に起因する到来周波数差の揺らぎを示す値である。たとえば、目標Ｔｇｔが船舶の場合、海洋上の多数の地点における波に関するデータが提供されている。データ取得部６８は、外部サーバ（図示せず）から、このような波に関するデータを目標揺動情報ＳＤとして取得することができるので、閾値設定部８６Ａは、目標揺動情報ＳＤに基づき、任意の地点における船舶の動揺の見積もりを実行して到来周波数差の揺らぎを示す値σ_{ｍｏｔｉｏｎ}を算出することができる。

閾値設定部８６Ａは、状態判定式（５），（６）のいずれか一方が満たされたことを検知したときに、目標Ｔｇｔの静止状態から移動状態への状態変化が生じたと判定することができる。最小値ＴＤＯＡ_ｍｉｎおよび最大値ＴＤＯＡ_ｍａｘを算出する場合と同様に、最小値ＦＤＯＡ_ｍｉｎおよび最大値ＦＤＯＡ_ｍａｘは、衛星軌道情報ＯＤおよび監視エリア情報ＡＤを用いて算出可能である。また、揺らぎσ_{ｍｏｔｉｏｎ}は、目標Ｔｇｔの位置と衛星Ｓｔｉ，Ｓｔｊの位置とが不変であるとの仮定の下、目標Ｔｇｔの動揺速度の分散値を含む目標揺動情報ＳＤに基づいて算出可能である。目標Ｔｇｔの動揺速度の分散値が与えられれば、マージン３σ_{ｍｏｔｉｏｎ}を算出することができる。

一方、目標Ｔｇｔの移動状態から静止状態への状態変化が生じる場合、平滑値＜ＦＤＯＡ＞_ｋの変動幅が一定範囲内に収まる。この場合、状態判定部８７は、平滑値＜ＦＤＯＡ＞_ｋが一定期間Ｔ２の間、数値範囲Δ（ＴＨ３，ＴＨ４）内に継続して収まることを検知したときに、目標Ｔｇｔの移動状態から静止状態への状態変化が生じたと判定する。一定期間Ｔ２の時間長は、パラメータ設定値である。より具体的には、状態判定部８７は、一定期間Ｔ２の間、次の状態判定式（７），（８）の双方が満たされることを継続して検知したときに、目標Ｔｇｔの移動状態から静止状態への状態変化が生じたと判定することができる。

以上に説明したように、状態判定部８３は、目標Ｔｇｔの動きの有無、すなわち、目標Ｔｇｔの状態が静止状態または移動状態のいずれであるかを判定し、その判定結果ＥＳ１を判定出力部８９に与える。判定出力部８９は、図３に示されるように、判定結果ＥＳ１を示す検出信号ＥＳを測位演算部６３Ａ，追尾演算部６４Ａおよびデータ移行部６５に供給する。

次に、測位演算部６３Ａ、追尾演算部６４Ａおよび追尾演算部６４Ａを詳細に説明する前に、実施の形態１および後述の実施の形態２〜４で使用されるカルマンフィルタおよび非線形最小二乗法について説明する。カルマンフィルタは、フィルタリング理論の一種である。

カルマンフィルタの目的は、時刻ｔ_ｋにおいて雑音の混入した観測値系列ｚ_１，…，ｚ_ｋ−１，ｚ_ｋが得られたとき、当該観測値系列と状態空間モデルとを用いて目標Ｔｇｔの静止状態および移動状態を示す物理量を推定することである。状態空間モデルは、目標Ｔｇｔの運動モデル（状態方程式）と観測モデル（観測方程式）とで構成される。運動モデルは次式（９）で与えられる。

ここで、ｘ_ｋは、時刻ｔ_ｋにおけるｎ行１列の状態ベクトルであり、観測不能な真の状態を表している。ｎは２以上の整数である。ｗ_ｋは、時刻ｔ_ｋにおけるｎ行１列の雑音ベクトルであり、Ｆ_ｋは、時刻ｔ_ｋにおけるｎ行ｎ列の状態遷移行列である。後述するように、状態ベクトルｘ_ｋの成分は、測位演算に適用される場合と追尾演算で適用される場合とで異なる。雑音ベクトルｗ_ｋの平均は零ベクトルとし、雑音ベクトルｗ_ｋの共分散行列はＱ_ｋで表現されるものとする。

一方、観測モデルは次式（１０）で与えられる。

ここで、ｚ_ｋは、時刻ｔ_ｋにおけるｍ行１列の観測ベクトルであり、ｖ_ｋは、時刻ｔ_ｋにおけるｍ行１列の雑音ベクトルであり、ｍは正整数である。雑音ベクトルｖ_ｋの平均は零ベクトルとし、雑音ベクトルｖ_ｋの共分散行列はＲ_ｋで表現されるものとする。また、ｈ［ｘ_ｋ］は、状態ベクトルｘ_ｋに関する観測関数であり、ｍ行１列のベクトルである。ｍ＝１のとき、ｚ_ｋ，ｖ_ｋ，ｈ［ｘ_ｋ］は、スカラー量となる。

カルマンフィルタによる推定処理（反復推定処理）は、予測処理、平滑化処理（更新処理）および外れ値判定処理を反復して実行することにより実現される。以下、ｘ（ｋ｜ｐ）は、時刻ｔ_ｐの時点での時刻ｔ_ｋの状態推定値からなるｎ行１列のベクトルを表すものとする。

予測処理では、次式（１１），（１２）に従い、事前状態推定ベクトルｘ（ｋ｜ｋ−１）および事前誤差共分散行列Ｐ（ｋ｜ｋ−１）が計算される。

ここで、Ｆ_ｋ ^Ｔは、状態遷移行列Ｆ_ｋに対する転置行列である。

平滑化処理（更新処理）では、次式（１３），（１４）に従い、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）が算出される。

ここで、ｅ_ｋは、観測値と状態推定値との間の残差からなる観測残差ベクトルである。また、Ｓ_ｋは観測残差共分散行列、Ｋ_ｋはカルマンゲイン、Ｈ_ｋは観測行列である。観測残差ベクトルｅ_ｋ、観測残差共分散行列Ｓ_ｋおよびカルマンゲインＫ_ｋは、それぞれ、次式（１５），（１６），（１７）で与えられる。

観測行列Ｈ_ｋは、次式（１８）で示されるヤコビアン行列（Ｊａｃｏｂｉａｎ ｍａｔｒｉｘ）によって定義される。

また、上記予測処理の後に、観測値が外れ値（ｏｕｔｌｉｅｒ ｖａｌｕｅｓ）か否かを判定する外れ値判定処理が実行される。外れ値判定処理では、まず、次式（１９）に基づいて、次式（２０）で定義される観測残差共分散行列Ｓ_ｋを用いた重み付き残差平方和Ｄ_ｋが算出される。

式（１９）は、上式（１５）で示される観測残差ベクトルｅ_ｋを用いれば、次式（２１）として表現可能である。

その後、次の判定式（２２）が成立するか否かが判定される。

ここで、γは、カイ二乗分布に従うゲートサイズパラメータである。

判定式（２２）が成立する場合には、観測値は外れ値ではないと判定され、平滑化処理（更新処理）が実行される。一方、判定式（２２）が成立しない場合には、観測値は外れ値であると判定されるので、平滑化処理は実行されない。この場合、上式（１３），（１４）の代わりに次式（２３）に従って、状態推定値ｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）が算出される。

図６は、カルマンフィルタによる反復推定処理の手順の一例を概略的に示すフローチャートである。まず、ステップＳＴ１１では、初期化が実行される。すなわち、時刻ｔ_０（ｋ＝０）における状態推定ベクトルｘ（０｜０）および事後誤差共分散行列Ｐ（０｜０）が与えられる。

次いで、時刻ｔ_１（ｋ＝１）における観測値の組すなわち観測ベクトルｚ_ｋが得られるまで待機し（ステップＳＴ１２のＮＯ）。観測ベクトルｚ_ｋが得られると（ステップＳＴ１２のＹＥＳ）、時刻ｔ_ｋ（ｋ＝１）について上記予測処理が実行される（ステップＳＴ１３）。この結果、事前状態推定ベクトルｘ（ｋ｜ｋ−１）および事前誤差共分散行列Ｐ（ｋ｜ｋ−１）が得られる。続けて、上記外れ値判定処理が実行される（ステップＳＴ１４）。その結果、観測ベクトルｚ_ｋにおける観測値が外れ値ではないと判定されたときは（ステップＳＴ１５のＮＯ）、上記平滑化処理が実行される（ステップＳＴ２０）。これにより、時刻ｔ_ｋ（ｋ＝１）における状態推定値ｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）が計算される。

平滑化処理（ステップＳＴ２０）の後、収束条件が満たされるか否かが判定される（ステップＳＴ２１）。収束条件は、各時刻ｔ_ｋにおける状態推定値ｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）が十分に収束していると推定される場合に満たされる条件である。収束条件が満たされない場合（ステップＳＴ２１のＮＯ）、事前状態推定ベクトルｘ（ｋ｜ｋ−１）および事前誤差共分散行列Ｐ（ｋ｜ｋ−１）が、ステップＳＴ２０で算出された状態推定値ｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）でそれぞれ置き換えられる（ステップＳＴ２２）。次いで、平滑化処理（ステップＳＴ２０）が再度実行される。時刻ｔ_ｋについて、最終的に収束条件が満たされた場合に（ステップＳＴ２１のＹＥＳ）、次のステップＳ２４が実行される。このように収束条件が満たされるまで、平滑化処理を繰り返し実行すること（巡回処理）によって、精度の高い状態推定値ｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を得ることができる。

ステップＳＴ２１の収束条件は、たとえば、各時刻ｔ_ｋについて平滑化処理があらかじめ設定された回数だけ繰り返し実行された場合、または、各時刻ｔ_ｋについて、直近のｊ回目に算出された状態推定値ｘ_ｊ（ｋ｜ｋ）とｊ−１回目に算出された状態推定値ｘ_ｊ−１（ｋ｜ｋ）との間のノルムが閾値以下となった場合に満たされるものとすることができる（ｊは、巡回処理における平滑化処理の繰り返し回数）。あるいは、直近のｊ回目に算出された事後誤差共分散行列Ｐ_ｊ（ｋ｜ｋ）とｊ−１回目に算出された事後誤差共分散行列Ｐ_ｊ−１（ｋ｜ｋ）との間のノルムが閾値以下となり、かつ、状態推定値ｘ_ｊ（ｋ｜ｋ），ｘ_ｊ−１（ｋ｜ｋ）間のノルムが閾値以下となった場合にステップＳＴ２１の収束条件が満たされてもよい。

一方、観測値が外れ値であると判定されたときは（ステップＳＴ１５のＹＥＳ）、ステップＳＴ２０の平滑化処理は実行されず、上式（２３）に従って、状態推定値ｘ（ｋ｜ｋ）および事前誤差共分散行列Ｐ（ｋ｜ｋ）が算出される（ステップＳＴ１６）。

ステップＳＴ１６またはステップＳＴ２１の後、あらかじめ定められた反復条件が満たされると判定されたとき（ステップＳＴ２４のＹＥＳ）、次の時刻ｔ_ｋ（ｋ＝２）について上記処理手順が反復して実行される（ステップＳＴ２５，ＳＴ１２以後）。最終的に、反復条件が満たされないと判定されたときに（ステップＳＴ２４のＮＯ）、反復推定処理は終了する。

次に、非線形最小二乗法について説明する。

まず、次式（２４）の観測モデル（観測方程式）が成立すると仮定する。

ここで、ｙ_ｋは、時刻ｔ_ｋにおけるｐ行１列の観測ベクトルであり、η_ｋは、時刻ｔ_ｋにおけるｐ行１列の雑音ベクトルであり、ｐは２以上の整数である。雑音ベクトルη_ｋの平均は零ベクトルとし、雑音ベクトルη_ｋの共分散行列（観測雑音共分散行列）はΩ_ｋで表現されるものとする。また、φ_ｋ［χ］は、状態ベクトルχに関する観測関数であり、ｐ行１列のベクトルである。

Ｋ回の計測結果（Ｋは２以上の整数）に対して、観測残差の大きさを評価する評価関数Ｊ（χ）が次式（２５）により定義される。

評価関数Ｊ（χ）は、雑音環境下での観測ベクトルｙ_ｋとベクトルφ_ｋ［χ］との間の距離（「マハラノビス距離（Ｍａｈａｌａｎｏｂｉｓ ｄｉｓｔａｎｃｅ）」と呼ばれる。）の二乗の総和を表している。

非線形最小二乗法の目的は、次式（２６）に示されるように評価関数Ｊ（χ）を最小にするχの近似解χ_ｍｉｎを求めることである。

仮に、観測関数φ_ｋ［χ］が、次式（２７）に示されるように線形でかつ行列Γ_ｋを用いて表現可能であるとする。

このとき、非線形最小二乗法による近似解χ_ｍｉｎは、次式（２８）で与えられる。

一方、観測関数φ_ｋ［χ］が非線形である場合には、再帰最小二乗法（Ｒｅｃｕｒｓｉｖｅ Ｌｅａｓｔ−Ｓｑｕａｒｅｓ，ＲＬＳ）アルゴリズムにより、近似解χ_ｍｉｎを求めることができる。たとえば、観測関数φ_ｋ［χ］が、近似的に、次式（２９）で示されるヤコビアン行列（関数行列）を用いて表現可能であるとする。

ＲＬＳアルゴリズムでは、たとえば、次の収束演算式（３０），（３１），（３２）を使用することができる（ｍは、再帰処理における繰り返し回数）。

ＲＬＳアルゴリズムに基づく再帰処理の手順は以下のとおりである。まず、初期化が実行される。すなわち、解の候補の初期値χ_０（ｍ＝０）が設定される。次に、繰り返し回数ｍ＝１，２，…について、上式（３０），（３１），（３２）に基づいた演算が反復して実行される。このとき、次式（３３）で示される反復終了条件が満たされたときに再帰処理が正常に終了する。

ここで、εは、微少な正の実数であり、Ｍは、繰り返し回数の上限値である。最終的に、反復終了条件を満たす解χ_ｍが、近似解χ_ｍｉｎとして与えられる。

次に、測位演算部６３Ａの構成および動作について説明する。

測位演算部６３Ａは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）の計測値の中から２以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた測位演算を実行することにより目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。測位演算部６３Ａは、互いに同期して得られる到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値の組を選択し、この組を観測ベクトルとして用いた測位演算を実行することができる。

あるいは、測位演算部６３Ａは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ），ＦＤＯＡ_２３（ｋ），ＦＤＯＡ_３１（ｋ）の計測値の中から、少なくとも１つの到来時間差および少なくとも１つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた測位演算を実行することにより目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を逐次算出する機能を有している。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、測位演算部６３Ａは、到来時間差ＴＤＯＡ_１２（ｋ）と到来周波数差ＦＤＯＡ_１２（ｋ）の計測値の組を選択し、この組を観測ベクトルとして用いた測位演算を実行することができる。

図７は、実施の形態１における測位演算部６３Ａ、追尾演算部６４Ａおよびデータ移行部６５からなる演算部の構成を概略的に示すブロック図である。

図７に示されるように、測位演算部６３Ａは、単測位部９１Ａと反復推定部９２Ａとを有する。単測位部９１Ａは、当該選択された計測値の組を観測ベクトルとして用いた非線形最小二乗法による測位演算を実行して目標Ｔｇｔの位置を示す測位ベクトルＰＯＳ_ｋを算出し、当該測位ベクトルＰＯＳ_ｋを反復推定部９２Ａに出力する。反復推定部９２Ａは、測位ベクトルＰＯＳ_ｋを用いたカルマンフィルタによる反復推定処理を実行して目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を算出する。

より具体的には、単測位部９１Ａは、上記した再帰最小二乗法（ＲＬＳ）アルゴリズムに従い、当該選択された計測値の組を観測ベクトルｙ_ｋとして用いた再帰処理を実行することができる。たとえば、ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値の組が観測ベクトルｙ_ｋとして選択されている場合、観測ベクトルｙ_ｋは、次式（３４）で与えられる。

このとき、時刻ｔ_ｋにおける観測関数φ_ｋ［χ＝ｐ_ｔｇｔ］は、次式（３５）で与えられる。

ここで、ｐ_ｔｇｔは、地球の重心を原点とする、目標Ｔｇｔの位置座標を示す位置ベクトルであり、ｐ_１（ｋ）は、時刻ｔ_ｋでの衛星Ｓｔ１の位置座標を示す位置ベクトルであり、ｐ_２（ｋ）は、時刻ｔ_ｋでの衛星Ｓｔ２の位置座標を示す位置ベクトルである。

観測雑音共分散行列Ω_ｋは、次式（３６）で与えられるものとする。

単測位部９１Ａは、複数回の計測結果に対して上記したＲＬＳアルゴリズムによる再帰処理を実行することにより、評価関数Ｊ（χ＝ｐ_ｔｇｔ）を最小にする近似解χ_ｍｉｎを測位ベクトルＰＯＳ_ｋとして算出することができる。単測位部９１Ａは、測位ベクトルＰＯＳ_ｋを反復推定部９２Ａに供給する。

単測位演算における理論的な測位誤差を表す誤差共分散行列Λは、次式（３６Ａ）で与えられる。

ここで、Φ_ｋ［χ］は、上式（２９）で示されるヤコビアン行列である。

次に、反復推定部９２Ａは、上式（９），（１０）で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を実行する。反復推定部９２Ａに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（３７），（３８），（３９）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度である。目標Ｔｇｔは静止状態にあるので、状態遷移行列Ｆ_ｋは単位行列である。また、目標Ｔｇｔの運動に不確かさがなく、雑音ベクトルｗ_ｋは零ベクトルとなる。よって、共分散行列Ｑ_ｋは零行列である。

また、反復推定部９２Ａに適用される観測モデル（観測方程式）では、観測ベクトルｚ_ｋ、観測行列Ｈ_ｋおよび共分散行列Ｒ_ｋは、たとえば、次式（４０），（４１），（４２）で与えられる。

ここで、式（４０）におけるＯＬＮ_ｋ，ＯＬＴ_ｋは、目標Ｔｇｔの観測位置を示す経度および緯度であり、式（４２）におけるＤＯＰは、精度劣化度（Ｄｉｌｕｔｉｏｎ ｏｆ Ｐｒｅｃｉｓｉｏｎ）と呼ばれ、本実施の形態では、単測位演算における理論的な測位誤差式によって与えられる。具体的には、ＤＯＰは、上式（３６Ａ）の誤差共分散行列Λに設定されればよい。

反復推定部９２Ａは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９２Ａは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）として出力する。上記のとおり、単測位部９１Ａは、複数回の計測結果（累積計測結果）を用いた非線形最小二乗法による測位演算を実行して高精度な測位ベクトルＰＯＳ_ｋを算出することができるので、反復推定部９２Ａは、高精度な測位ベクトルＰＯＳ_ｋに基づき、精度の高い測位演算結果を生成することができる。

次に、追尾演算部６４Ａの構成および動作について説明する。

追尾演算部６４Ａは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）の計測値の中から２以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を逐次算出する機能を有する。あるいは、追尾演算部６４Ａは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ），ＦＤＯＡ_２３（ｋ），ＦＤＯＡ_３１（ｋ）の計測値の中から、少なくとも１つの到来時間差および少なくとも１つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を逐次算出する機能を有している。

より具体的には、追尾演算部６４Ａは、図７に示されるように単測位部９３Ａと反復推定部９４Ａとを有している。単測位部９３Ａの構成および機能は、単測位部９１Ａのそれらと同じである。反復推定部９４Ａは、単測位部９３Ａで生成された測位ベクトルＰＯＳ_ｋを用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出することができる。

反復推定部９４Ａは、上式（９），（１０）で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部９４Ａに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（４３），（４４），（４５）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度、ＶＬＮ_ｋは経度方向の速度、ＬＴ_ｋは緯度方向の速度であり、Ｔ＝ｔ_ｋ−ｔ_ｋ−１である。また、Ｉ_２は２行２列の単位行列であり、ｑはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。目標Ｔｇｔは移動状態にあるので、状態遷移行列Ｆ_ｋは、等速直進運動を想定して設定されている。

また、反復推定部９４Ａに適用される観測モデル（観測方程式）では、観測ベクトルｚ_ｋ、観測行列Ｈ_ｋおよび共分散行列Ｒ_ｋは、たとえば、次式（４６），（４７），（４８）で与えられる。

ここで、ＯＬＮ_ｋ，ＯＬＴ_ｋは、目標Ｔｇｔの観測位置を示す経度および緯度であり、式（４８）におけるＤＯＰは、精度劣化度（Ｄｉｌｕｔｉｏｎ ｏｆ Ｐｒｅｃｉｓｉｏｎ）と呼ばれ、本実施の形態では、単測位演算における理論的な測位誤差式によって与えられる。具体的には、ＤＯＰは、上式（３６Ａ）の誤差共分散行列Λに設定されればよい。

反復推定部９４Ａは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９４Ａは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）として出力する。

次に、状態検出部６２Ａが目標Ｔｇｔの状態変化を検出した場合の測位演算部６３Ａ、追尾演算部６４Ａおよびデータ移行部６５の動作について説明する。

目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、測位演算部６３Ａの反復推定部９２Ａは、検出信号ＥＳに応じて、直近の測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｐｏ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して追尾演算部６４Ａに供給する。

このとき、追尾演算部６４Ａの反復推定部９４Ａは、測位演算部６３Ａから供給された状態推定ベクトルｘ_ｐｏ（ｋ）および事後誤差共分散行列Ｐ_ｐｏ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９４Ａは、それら初期値ｘ（０｜０），Ｐ（０｜０）と測位ベクトルＰＯＳ_ｋとを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出することができる。

この場合に、反復推定部９４Ａで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、次式（４９），（５０）に示すように設定可能である。

ここで、０_２×１は２行１列の零行列、０_２×２は２行２列の零行列である。式（４９）の状態推定ベクトルｘ（０｜０）は、目標Ｔｇｔの初期速度が零となるように設定されている。また、式（５０）の初期値Ｐ（０｜０）は、推定速度成分に誤差はないとの仮定に基づいて設定されている。

一方、目標Ｔｇｔの移動状態から静止状態への状態変化が検出されたとき、追尾演算部６４Ａの反復推定部９４Ａは、検出信号ＥＳに応じて、直近の追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して測位演算部６３Ａに供給する。

このとき、測位演算部６３Ａにおける反復推定部９２Ａは、追尾演算部６４Ａから供給された状態推定ベクトルｘ_ｔｒｋ（ｋ）および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９２Ａは、それら初期値ｘ（０｜０），Ｐ（０｜０）と測位ベクトルＰＯＳ_ｋとを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｐｏ（ｋ）を算出することができる。

この場合に、反復推定部９２Ａで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、次式（５１），（５２）に示すように設定可能である。

ここで、［ｘ_ｔｒｋ（ｋ）］_１，１は、状態推定ベクトルｘ_ｔｒｋ（ｋ）の１行１列目成分（１番目成分）、［ｘ_ｔｒｋ（ｋ）］_２，１は、状態推定ベクトルｘ_ｔｒｋ（ｋ）の２行１列目成分（２番目成分）である。また、［Ｐ_ｔｒｋ（ｋ）］_１，１は、事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）の１行１列目成分、［Ｐ_ｔｒｋ（ｋ）］_２，１は、事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）の２行１列目成分、［Ｐ_ｔｒｋ（ｋ）］_１，２は、事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）の１行２列目成分、［Ｐ_ｔｒｋ（ｋ）］_２，２は、事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）の２行２列目成分である。

以上に説明したように、目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、追尾演算部６４Ａの反復推定部９４Ａは、直近の測位演算結果ｘ_ｐｏ（ｋ），Ｐ_ｐｏ（ｋ）を初期値として用いかつ計測値を用いた追尾演算を実行することができる。逆に、目標Ｔｇｔの移動状態から静止状態への状態変化が検出されたときは、測位演算部６３Ａの反復推定部９２Ａは、直近の追尾演算結果ｘ_ｔｒｋ（ｋ），Ｐ_ｔｒｋ（ｋ）を初期値として用いかつ計測値を用いた測位演算を実行することができる。

上記した監視処理部１１Ａのハードウェア構成は、たとえば、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ），ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）またはＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）などの半導体集積回路を有する単数または複数のプロセッサで実現されればよい。あるいは、監視処理部１１Ａのハードウェア構成は、不揮発性メモリから読み出されたソフトウェアまたはファームウェアのプログラムコードを実行する、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）またはＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）などの演算装置を含む単数または複数のプロセッサで実現されてもよい。ＤＳＰなどの半導体集積回路とＣＰＵなどの演算装置との組み合わせを含む単数または複数のプロセッサで監視処理部１１Ａのハードウェア構成が実現されてもよい。

図８は、監視処理部１１Ａの機能を実現するハードウェア構成例である信号処理装置７０の概略構成を示すブロック図である。図８に示される信号処理装置７０は、プロセッサ７１、メモリ７２、入力インタフェース部７３、出力インタフェース部７４、通信回路７５および信号路７６を含んで構成されている。信号路７６は、プロセッサ７１、メモリ７２、入力インタフェース部７３、出力インタフェース部７４および通信回路７５を相互に接続するためのバスである。入力インタフェース部７３は、外部から入力されたディジタル受信信号Ｄ１，Ｄ２，Ｄ３を信号路７６を介してプロセッサ７１に転送する機能を有する。プロセッサ７１は、転送されたディジタル受信信号Ｄ１，Ｄ２，Ｄ３に基づいて測位演算および追尾演算を実行することで目標Ｔｇｔの静止状態および移動状態を示す状態推定値を算出することができる。また、プロセッサ７１は、算出された状態推定値を基に画像情報ＤＤを生成し、この画像情報ＤＤを信号路７５および出力インタフェース部７４を介して外部に出力することができる。

ここで、通信回路７５は、外部サーバ（図示せず）と通信して状態検出部６２Ａでの処理に必要なデータを受信する回路である。また、メモリ７２は、プロセッサ７１がディジタル信号処理を実行する際に使用されるデータ記憶領域である。プロセッサ７１がＣＰＵなどの演算装置を内蔵する場合には、メモリ７２は、プロセッサ７１により実行されるソフトウェアまたはファームウェアのプログラムコードを記憶するデータ記憶領域を有する。メモリ７２としては、たとえば、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）およびＳＤＲＡＭ（Ｓｙｎｃｈｒｏｎｏｕｓ Ｄｙｎａｍｉｃ Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）などの半導体メモリを使用することが可能である。

以上に説明したように実施の形態１の目標監視装置２においては、状態検出部６２Ａは、複数の計測値のうちの少なくとも１つの計測値の時間変動を監視して目標Ｔｇｔの動きの有無を検出し、その検出結果を示す検出信号ＥＳを測位演算部６３Ａ、追尾演算部６４Ａおよびデータ移行部６５に供給するので、測位演算部６３Ａおよび追尾演算部６４Ａは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、監視処理部１１Ａは、目標Ｔｇｔの静止状態から移動状態への状態変化、あるいは目標Ｔｇｔの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。

また、目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、追尾演算部６４Ａの反復推定部９４Ａは、直近の測位演算結果ｘ_ｐｏ（ｋ），Ｐ_ｐｏ（ｋ）を初期値として用いかつ計測値の組を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Ｔｇｔに対する高い追尾精度を確保することができる。

実施の形態２．
次に、本発明に係る実施の形態２について説明する。図９は、本発明に係る実施の形態２における測位演算部６３Ｂ、追尾演算部６４Ｂおよびデータ移行部６５からなる演算部の構成を概略的に示すブロック図である。本実施の形態の目標監視システムおよび目標監視装置（監視処理部を含む。）の構成は、図３の測位演算部６３Ａおよび追尾演算部６４Ａに代えて、図９の測位演算部６３Ｂおよび追尾演算部６４Ｂを有する点を除いて、上記実施の形態１の目標監視システム１および目標監視装置２の構成と同じである。

本実施の形態の測位演算部６３Ｂは、図９に示されるように反復推定部９２Ｂを有する。反復推定部９２Ｂは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）の計測値の中から２以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、反復推定部９２Ｂは、互いに同期して得られる到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値を選択し、これら計測値の組を観測ベクトルとして用いた測位演算を実行することができる。

あるいは、反復推定部９２Ｂは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ），ＦＤＯＡ_２３（ｋ），ＦＤＯＡ_３１（ｋ）の計測値の中から、少なくとも１つの到来時間差および少なくとも１つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、反復推定部９２Ｂは、互いに同期して得られる到来時間差ＴＤＯＡ_１２（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ）の計測値を選択し、これら計測値の組を観測ベクトルとして用いた測位演算を実行することができる。

より具体的には、測位演算部６３Ｂの反復推定部９２Ｂは、上式（９），（１０）で示される状態空間モデルを用いてカルマンフィルタによる反復推定処理を測位演算として実行することができる。反復推定部９２Ｂに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（５３），（５４），（５５）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度である。目標Ｔｇｔは静止状態にあるので、状態遷移行列Ｆ_ｋは単位行列である。また、目標Ｔｇｔの運動に不確かさがなく、雑音ベクトルｗ_ｋは零ベクトルとなる。よって、共分散行列Ｑ_ｋは零行列である。

また、反復推定部９２Ｂに適用される観測モデル（観測方程式）では、観測ベクトルｚ_ｋ、観測関数ｈ［ｘ_ｋ］および共分散行列Ｒ_ｋは、たとえば、次式（５６），（５７），（５８）で与えられる。

ここで、観測関数ｈ［ｘ_ｋ］は、上式（１８）に従い、ヤコビアン行列Ｈ_ｋに変換される。

反復推定部９２Ｂは、上記状態空間モデルを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９２Ｂは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）として出力する。

一方、追尾演算部６４Ｂは、図９に示されるように反復推定部９４Ｂを有する。反復推定部９４Ｂは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）の計測値の中から２以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を逐次算出する機能を有する。当該選択された計測値は、互いに同期して、または互いにほぼ同期して得られる値である。あるいは、追尾演算部６４Ｂは、計測部６１から供給された到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ），ＴＤＯＡ_３１（ｋ）および到来周波数差ＦＤＯＡ_１２（ｋ），ＦＤＯＡ_２３（ｋ），ＦＤＯＡ_３１（ｋ）の計測値の中から、少なくとも１つの到来時間差および少なくとも１つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を逐次算出する機能を有する。当該選択された計測値は、互いに同期して、または互いにほぼ同期して得られる値である。

より具体的には、反復推定部９４Ｂは、上式（９），（１０）で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部９４Ｂに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（５９），（６０），（６１）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度、ＶＬＮ_ｋは経度方向の速度、ＬＴ_ｋは緯度方向の速度であり、Ｔ＝ｔ_ｋ−ｔ_ｋ−１である。また、Ｉ_２は２行２列の単位行列であり、ｑはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。

また、反復推定部９２Ｂに適用される観測モデル（観測方程式）では、観測ベクトルｚ_ｋ、観測関数ｈ［ｘ_ｋ］および共分散行列Ｒ_ｋは、たとえば、次式（６２），（６３），（６４）で与えられる。

ここで、観測関数ｈ［ｘ_ｋ］は、上式（１８）に従い、ヤコビアン行列Ｈ_ｋに変換される。

反復推定部９４Ｂは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９４Ｂは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）として出力する。

次に、目標Ｔｇｔの状態変化が検出された場合の測位演算部６３Ｂ、追尾演算部６４Ｂおよびデータ移行部６５の動作について説明する。

目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、測位演算部６３Ｂの反復推定部９２Ｂは、検出信号ＥＳに応じて、直近の測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｐｏ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して追尾演算部６４Ｂに供給する。

このとき、追尾演算部６４Ｂの反復推定部９４Ｂは、測位演算部６３Ｂから供給された状態推定ベクトルｘ_ｐｏ（ｋ）および事後誤差共分散行列Ｐ_ｐｏ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９４Ｂは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測ベクトルとを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出することができる。この場合に、反復推定部９２Ｂで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、上式（４９），（５０）に示すように設定可能である。

一方、目標Ｔｇｔの移動状態から静止状態への状態変化が検出されたとき、追尾演算部６４Ｂの反復推定部９４Ｂは、検出信号ＥＳに応じて、直近の追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して測位演算部６３Ｂに供給する。

このとき、測位演算部６３Ｂの反復推定部９２Ｂは、追尾演算部６４Ｂから供給された状態推定ベクトルｘ_ｔｒｋ（ｋ）および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９２Ｂは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測ベクトルとを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｐｏ（ｋ）を算出することができる。この場合に、反復推定部９２Ｂで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、次式（５１），（５２）に示すように設定可能である。

以上に説明したように実施の形態２の目標監視装置は、実施の形態１の目標監視装置２と同様に状態検出部６２Ａを有するので、測位演算部６３Ｂおよび追尾演算部６４Ｂは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態２の監視処理部は、目標Ｔｇｔの静止状態から移動状態への状態変化、あるいは目標Ｔｇｔの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、追尾演算部６４Ｂの反復推定部９４Ｂは、直近の測位演算結果ｘ_ｐｏ（ｋ），Ｐ_ｐｏ（ｋ）を初期値として用い、かつ計測値の組を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Ｔｇｔに対する高い追尾精度を確保することができる。

実施の形態３．
次に、本発明に係る実施の形態３について説明する。上記実施の形態１，２では、時刻ｔ_ｋ毎の観測ベクトルの成分として使用される複数の計測値（たとえば、到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値）は、互いに同期して、または互いにほぼ同期して得られる観測値である。これに対し、本実施の形態では、複数の計測値は、互いに非同期に得られる観測値、言い換えれば、同時刻に同時に得られない観測値である。

たとえば、受信アンテナ（受信局）が３つある場合には、到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値を互いに同期する形で同時に得ることができる。これに対し、使用可能な受信アンテナ（受信局）が２つのみである場合には、衛星Ｓｔ１，Ｓｔ２間の計測値と衛星Ｓｔ２，Ｓｔ３間の計測値との組を得るために、アンテナ指向方向の切り替えを行う必要がある。しかしながら、その切り替えに時間ずれが生じると、それら２種類の計測値が互いに非同期に得られることがある。また、到来時間差ＴＤＯＡ_１２（ｋ），ＴＤＯＡ_２３（ｋ）の計測値が同時に得られたとしても、これら計測値の一方が外れ値であると判定された場合など何らかの理由で欠落する可能性がある。本実施の形態の監視処理部は、非同期に得られた複数の計測値を用いて測位演算および追尾演算を実行することができる。

図１０は、本発明に係る実施の形態３における監視処理部１１Ｃの概略構成を示すブロック図である。本実施の形態の目標監視システムの構成は、図３の監視処理部１１Ａに代えて図１０の監視処理部１１Ｃを有する点を除いて、上記実施の形態１の目標監視システム１の構成と同じである。監視処理部１１Ｃは、ディジタル受信信号Ｄ１，Ｄ２，Ｄ３に基づいて電波の到来時間差ＴＤＯＡ_ｉｊ（ｋ）および到来周波数差ＦＤＯＡ_ｉｊ（ｋ）の計測値を逐次算出する計測部６１Ｃを備える。到来時間差ＴＤＯＡ_ｉｊ（ｋ）は、到来時間差ＴＤＯＡ_１２（ｋ）またはＴＤＯＡ_２３（ｋ）のいずれか一方であり、到来周波数差ＦＤＯＡ_ｉｊ（ｋ）は、到来周波数差ＦＤＯＡ_１２（ｋ）またはＦＤＯＡ_２３（ｋ）のいずれか一方である。たとえば、時刻ｔ_ｋ１に到来時間差ＴＤＯＡ_１２（ｋ１）の計測値が得られ、時刻ｔ_ｋ２（≠ｔ_ｋ１）に到来時間差ＴＤＯＡ_２３（ｋ２）の計測値が得られる場合、到来時間差ＴＤＯＡ_１２（ｋ１），ＴＤＯＡ_２３（ｋ２）の計測値は、互いに非同期に得られた観測値である。

監視処理部１１Ｃは、さらに、到来時間差ＴＤＯＡ_ｉｊ（ｋ）および到来周波数差ＦＤＯＡ_ｉｊ（ｋ）の計測値のうちのいずれかの計測値の時間変動を常時監視して目標Ｔｇｔの動きの有無を検出する状態検出部６２Ａと、目標Ｔｇｔの動きが無いことが検出されたときは、当該計測値を用いた測位演算を実行して目標Ｔｇｔの推定位置を示す状態推定値を逐次算出する測位演算部６３Ｃと、目標Ｔｇｔの動きが有ることが検出されたきは、当該計測値を用いた追尾演算を実行して目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）の双方を示す状態推定値を算出する追尾演算部６４Ｃと、測位演算部６３Ｃおよび追尾演算部６４Ｃの一方から他方へ演算結果のデータを移行させるデータ移行部６５と、状態検出部６２Ａでの処理に必要なデータを格納するデータ記憶部６７とを有する。

図１１は、実施の形態３における測位演算部６３Ｃ、追尾演算部６４Ｃおよびデータ移行部６５からなる演算部の構成を概略的に示すブロック図である。本実施の形態の測位演算部６３Ｃは反復推定部９２Ｃを有する。反復推定部９２Ｃは、計測部６１で非同期に得られた到来時間差ＴＤＯＡ_１２（ｋ１），ＴＤＯＡ_２３（ｋ２）の計測値（ｋ１≠ｋ２）の組を用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を算出する機能を有する。

あるいは、反復推定部９２Ｃは、計測部６１で非同期に得られた到来時間差ＴＤＯＡ_１２（ｋ１）と到来周波数差ＦＤＯＡ_２３（ｋ２）の計測値（ｋ１≠ｋ２）の組を用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を算出してもよい。

より具体的には、測位演算部６３Ｃの反復推定部９２Ｃは、上式（９），（１０）で示される状態空間モデルを用いてカルマンフィルタによる反復推定処理を測位演算として実行することができる。反復推定部９２Ｃに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（６５），（６６），（６７）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度である。目標Ｔｇｔは静止状態にあるので、状態遷移行列Ｆ_ｋは単位行列である。また、目標Ｔｇｔの運動に不確かさがなく、雑音ベクトルｗ_ｋは零ベクトルとなる。よって、共分散行列Ｑ_ｋは零行列である。

また、反復推定部９２Ｃに適用される観測モデル（観測方程式）では、観測値ｚ_ｋ、観測関数ｈ［ｘ_ｋ］および共分散Ｒ_ｋは、たとえば、次式（６８），（６９），（７０）で与えられる。

ここで、式（６８），（６９）における（ｉ，ｊ）は、各時刻ｔ_ｋで（１，２）または（２，３）のいずれかである。たとえば、時刻ｔ_ｋ１では、到来時間差ＴＤＯＡ_１２（ｋ１）の観測値（計測値）が得られ（ｉ＝１；ｊ＝２）、時刻ｔ_ｋ２ではＴＤＯＡ_２３（ｋ２）の観測値（計測値）が得られる（ｉ＝２；ｊ＝３）。

反復推定部９２Ｃは、上記状態空間モデルを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９２Ｃは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）として出力する。

一方、追尾演算部６４Ｃは、図１１に示されるように反復推定部９４Ｃを有する。反復推定部９４Ｃは、計測部６１で非同期に得られた到来時間差ＴＤＯＡ_１２（ｋ１），ＴＤＯＡ_２３（ｋ２）の計測値（ｋ１≠ｋ２）の組を用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を逐次算出する機能を有する。

あるいは、追尾演算部６４Ｃは、計測部６１で非同期に得られた到来時間差ＴＤＯＡ_１２（ｋ１）および到来周波数差ＦＤＯＡ_２３（ｋ２）の計測値（ｋ１≠ｋ２）の組を用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Ｔｇｔの推定位置および移動状態（たとえば推定速度）を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出してもよい。

より具体的には、反復推定部９４Ｃは、上式（９），（１０）で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部９４Ｃに適用される運動モデル（状態方程式）では、状態ベクトルｘ_ｋ、状態遷移行列Ｆ_ｋおよび共分散行列Ｑ_ｋは、たとえば、次式（７１），（７２），（７３）で与えられる。

ここで、ＬＮ_ｋ，ＬＴ_ｋは、目標Ｔｇｔの真の位置を示す経度および緯度、ＶＬＮ_ｋは経度方向の速度、ＬＴ_ｋは緯度方向の速度であり、Ｔ＝ｔ_ｋ−ｔ_ｋ−１である。また、Ｉ_２は２行２列の単位行列であり、ｑはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。

また、反復推定部９４Ｃに適用される観測モデル（観測方程式）では、観測値ｚ_ｋ、観測関数ｈ［ｘ_ｋ］および共分散Ｒ_ｋは、たとえば、次式（７４），（７５），（７６）で与えられる。

ここで、式（７４），（７５）における（ｉ，ｊ）は、各時刻ｔ_ｋで（１，２）または（２，３）のいずれかである。たとえば、時刻ｔ_ｋ１では、到来時間差ＴＤＯＡ_１２（ｋ１）の計測値が得られ（ｉ＝１；ｊ＝２）、時刻ｔ_ｋ２ではＴＤＯＡ_２３（ｋ２）の計測値が得られる（ｉ＝２；ｊ＝３）。

反復推定部９４Ｃは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）を算出することができる。反復推定部９４Ｃは、当該状態推定ベクトルｘ（ｋ｜ｋ）を、目標Ｔｇｔの推定位置および移動状態を示す状態推定ベクトルｘ_ｔｒｋ（ｋ）として出力する。

次に、目標Ｔｇｔの状態変化が検出された場合の測位演算部６３Ｃ、追尾演算部６４Ｃおよびデータ移行部６５の動作について説明する。

目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、測位演算部６３Ｃの反復推定部９２Ｃは、検出信号ＥＳに応じて、直近の測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｐｏ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して追尾演算部６４Ｃに供給する。

このとき、追尾演算部６４Ｃの反復推定部９４Ｃは、測位演算部６３Ｃから供給された状態推定ベクトルｘ_ｐｏ（ｋ）および事後誤差共分散行列Ｐ_ｐｏ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９４Ｃは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測値（計測値）とを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出することができる。この場合に、反復推定部９２Ｂで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、上式（４９），（５０）に示すように設定可能である。

一方、目標Ｔｇｔの移動状態から静止状態への状態変化が検出されたとき、追尾演算部６４Ｃの反復推定部９４Ｃは、検出信号ＥＳに応じて、直近の追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して測位演算部６３Ｃに供給する。

このとき、測位演算部６３Ｃの反復推定部９２Ｃは、追尾演算部６４Ｃから供給された状態推定ベクトルｘ_ｔｒｋ（ｋ）および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９２Ｃは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測値（計測値）とを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｐｏ（ｋ）を算出することができる。この場合に、反復推定部９２Ｂで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、次式（５１），（５２）に示すように設定可能である。

以上に説明したように実施の形態３の目標監視装置は、実施の形態１の目標監視装置２と同様に状態検出部６２Ａを有するので、測位演算部６３Ｃおよび追尾演算部６４Ｃは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態３の監視処理部１１Ｃは、目標Ｔｇｔの静止状態から移動状態への状態変化、あるいは目標Ｔｇｔの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、追尾演算部６４Ｃの反復推定部９４Ｃは、直近の測位演算結果ｘ_ｐｏ（ｋ），Ｐ_ｐｏ（ｋ）を初期値として用い、かつ計測値を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Ｔｇｔに対する高い追尾精度を確保することができる。

実施の形態４．
次に、本発明に係る実施の形態４について説明する。上記実施の形態３では、測位演算部６３Ｃは、各時刻ｔ_ｋに計測値が入力される度に測位演算を逐次的に実行する。これに対し、実施の形態４の測位演算は、複数回分の計測値を一括して処理するバッチ型の測位演算である。

図１２は、本発明に係る実施の形態４における測位演算部６３Ｄ、追尾演算部６４Ｃおよびデータ移行部６５からなる演算部の構成を概略的に示すブロック図である。本実施の形態の目標監視システムおよび目標監視装置の構成は、図１０の測位演算部６３Ｃに代えて図１０の測位演算部６３Ｄを有する点を除いて、上記実施の形態３の目標監視システムおよび目標監視装置の構成と同じである。

本実施の形態の測位演算部６３Ｄには、計測部６１Ｃ（図１０）から到来時間差ＴＤＯＡ_ｉｊ（ｋ）および到来周波数差ＦＤＯＡ_ｉｊ（ｋ）の計測値が入力される。到来時間差ＴＤＯＡ_ｉｊ（ｋ）は、到来時間差ＴＤＯＡ_１２（ｋ）またはＴＤＯＡ_２３（ｋ）のいずれか一方であり、到来周波数差ＦＤＯＡ_ｉｊ（ｋ）は、到来周波数差ＦＤＯＡ_１２（ｋ）またはＦＤＯＡ_２３（ｋ）のいずれか一方である。たとえば、時刻ｔ_ｋ１に到来時間差ＴＤＯＡ_１２（ｋ１）の計測値が得られ、時刻ｔ_ｋ２（≠ｔ_ｋ１）に到来時間差ＴＤＯＡ_２３（ｋ２）の計測値が得られる場合、到来時間差ＴＤＯＡ_１２（ｋ１），ＴＤＯＡ_２３（ｋ２）の計測値は、互いに非同期に得られた観測値である。

本実施の形態の測位演算部６３Ｄは、図１２に示されるようにバッファメモリ１００とバッチ型反復推定部９２Ｄとを有している。バッファメモリ１００は、少なくとも、時刻ｔ_{ｋ−Ｋ＋１}から時刻ｔ_ｋ−１までのＫ−１回分の過去の計測値を記憶する容量を有している。バッチ型反復推定部９２Ｄは、時刻ｔ_ｋにおける現在の計測値が得られる度に、バッファメモリ１００からＫ−１回分の過去の計測値を読み出す。そして、バッチ型反復推定部９２Ｄは、当該過去の計測値と当該現在の計測値とを一括して用いた非線形最小二乗法によるバッチ処理を測位演算として実行することにより、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）を算出することができる。

より具体的には、バッチ型反復推定部９２Ｄは、上記した再帰最小二乗法（ＲＬＳ）アルゴリズムに従い、次式（７７）に示される観測値ｙ_ｋを用いた再帰処理を実行することができる。

このとき、時刻ｔ_ｋにおける観測関数φ_ｋ［χ＝ｐ_ｔｇｔ］は、次式（７８）で与えられる。

ここで、ｐ_ｔｇｔは、地球の重心を原点とする、目標Ｔｇｔの位置座標を示す位置ベクトルであり、ｐ_ｉ（ｋ）は、時刻ｔ_ｋでの衛星Ｓｔｉの位置座標を示す位置ベクトルであり、ｐ_ｊ（ｋ）は、時刻ｔ_ｋでの衛星Ｓｔｊの位置座標を示す位置ベクトルである。

観測雑音共分散Ω_ｋは、次式（７９）で与えられるものとする。

バッチ型反復推定部９２Ｄは、Ｋ回の計測結果に対して上記したＲＬＳアルゴリズムによる再帰処理を実行することにより、評価関数Ｊ（χ＝ｐ_ｔｇｔ）を最小にする近似解χ_ｍｉｎを、目標Ｔｇｔの推定位置を示す状態推定ベクトルｘ_ｐｏ（ｋ）として算出することができる。

次に、目標Ｔｇｔの状態変化が検出された場合の測位演算部６３Ｄ、追尾演算部６４Ｃおよびデータ移行部６５の動作について説明する。

目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、測位演算部６３Ｄの反復推定部９２Ｄは、検出信号ＥＳに応じて、直近の測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｐｏ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して追尾演算部６４Ｃに供給する。

このとき、追尾演算部６４Ｃの反復推定部９４Ｃは、測位演算部６３Ｄから供給された状態推定ベクトルｘ_ｐｏ（ｋ）および事後誤差共分散行列Ｐ_ｐｏ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９４Ｃは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測値（計測値）とを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｔｒｋ（ｋ）を算出することができる。

一方、目標Ｔｇｔの移動状態から静止状態への状態変化が検出されたとき、追尾演算部６４Ｃの反復推定部９４Ｃは、検出信号ＥＳに応じて、直近の追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ）（＝ｘ（ｋ｜ｋ））および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）（＝Ｐ（ｋ｜ｋ））をデータ移行部６５を介して測位演算部６３Ｄに供給する。

このとき、測位演算部６３Ｄの反復推定部９２Ｄは、追尾演算部６４Ｃから供給された状態推定ベクトルｘ_ｔｒｋ（ｋ）および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ）に基づいて、状態推定ベクトルｘ（ｋ｜ｋ）および事後誤差共分散行列Ｐ（ｋ｜ｋ）の初期値ｘ（０｜０），Ｐ（０｜０）を生成する。そして、反復推定部９２Ｄは、それら初期値ｘ（０｜０），Ｐ（０｜０）と観測値（計測値）とを用いたカルマンフィルタによる反復推定処理（図６）を実行することにより、時刻ｔ_ｋ毎の状態推定ベクトルｘ_ｐｏ（ｋ）を算出することができる。この場合に、反復推定部９２Ｄで使用される初期値ｘ（０｜０），Ｐ（０｜０）は、たとえば、次式（５１），（５２）に示すように設定可能である。

以上に説明したように実施の形態４の目標監視装置は、実施の形態１の目標監視装置２と同様に状態検出部６２Ａを有するので、測位演算部６３Ｄおよび追尾演算部６４Ｃは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態４の監視処理部は、目標Ｔｇｔの静止状態から移動状態への状態変化、あるいは目標Ｔｇｔの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Ｔｇｔの静止状態から移動状態への状態変化が検出されたとき、追尾演算部６４Ｃの反復推定部９４Ｃは、直近の測位演算結果ｘ_ｐｏ（ｋ），Ｐ_ｐｏ（ｋ）を初期値として用い、かつ計測値を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Ｔｇｔに対する高い追尾精度を確保することができる。

実施の形態５．
次に、本発明に係る実施の形態５について説明する。実施の形態１の状態検出部６２Ａ（図４）においては、閾値設定部８２Ａは、監視エリアＳＡを定める監視エリア情報ＡＤを用いて閾値ＴＨ１，ＴＨ２を設定し、閾値設定部８６Ａも、監視エリア情報ＡＤを用いて閾値ＴＨ３，ＴＨ４を設定している。実施の形態５，６では、監視エリア情報ＡＤに代えて、直近の測位演算結果または直近の追尾演算結果を用いて閾値ＴＨ１，ＴＨ２，ＴＨ３，ＴＨ４が設定される。

図１３は、実施の形態５における監視処理部１１Ｅの概略構成を示すブロック図である。監視処理部１１Ｅの構成は、図３の状態検出部６２Ａに代えて図１３の状態検出部６２Ｅおよび遅延素子７８，７９を有する点を除いて、実施の形態１の監視処理部１１Ａの構成と同じである。図１３に示されるように、一方の遅延素子７８は、現在の時刻ｔ_ｋよりも１回前の時刻ｔ_ｋ−１での測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ−１）を状態検出部６２Ｅに供給し、他方の遅延素子７９は、時刻ｔ_ｋ−１での追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ−１）を状態検出部６２Ｅに供給する。

図１４は、実施の形態５における状態検出部６２Ｅの概略構成を示すブロック図である。状態検出部６２Ｅの構成は、図４の閾値設定部８２Ａ，８６Ａに代えて図１４の閾値設定部８２Ｅ，８６Ｅを有する点を除いて、上記状態検出部６２Ａの構成と同じである。閾値設定部８２Ｅは、測位演算結果ｘ_ｐｏ（ｋ−１）または追尾演算結果ｘ_ｔｒｋ（ｋ−１）のいずれか一方と衛星軌道情報ＯＤとに基づき、上式（１），（２）に示される閾値ＴＨ１，ＴＨ２の組を設定することができる。また、閾値設定部８６Ｅは、測位演算結果ｘ_ｐｏ（ｋ−１）または追尾演算結果ｘ_ｔｒｋ（ｋ−１）のいずれか一方と衛星軌道情報ＯＤと目標揺動情報ＳＤとに基づき、上式（５），（６）に示される閾値ＴＨ３，ＴＨ４の組を設定することができる。

以上に説明したように実施の形態５では、目標Ｔｇｔに関する状態推定ベクトルｘ_ｐｏ（ｋ−１），ｘ_ｔｒｋ（ｋ−１）を使用して閾値ＴＨ１，ＴＨ２，ＴＨ３，ＴＨ４が設定されるので、目標Ｔｇｔの状態変化を検知するための信頼性の高い数値範囲Δ（ＴＨ１，ＴＨ２），Δ（ＴＨ３，ＴＨ４）を決定することが可能となる。

実施の形態６．
次に、本発明に係る実施の形態６について説明する。図１５は、実施の形態６における監視処理部１１Ｆの概略構成を示すブロック図である。監視処理部１１Ｆの構成は、図３の状態検出部６２Ａに代えて図１５の状態検出部６２Ｆおよび遅延素子７８，７９，９１，９２を有する点を除いて、実施の形態１の監視処理部１１Ａの構成と同じである。図１５に示されるように、遅延素子７８，９１は、現在の時刻ｔ_ｋよりも１回前の時刻ｔ_ｋ−１での測位演算結果である状態推定ベクトルｘ_ｐｏ（ｋ−１）および事後誤差共分散行列Ｐ_ｐｏ（ｋ−１）を状態検出部６２Ｅに供給し、遅延素子７９，９２は、時刻ｔ_ｋ−１での追尾演算結果である状態推定ベクトルｘ_ｔｒｋ（ｋ−１）および事後誤差共分散行列Ｐ_ｔｒｋ（ｋ−１）を状態検出部６２Ｅに供給する。

図１６は、実施の形態６における状態検出部６２Ｆの概略構成を示すブロック図である。状態検出部６２Ｆの構成は、図４の閾値設定部８２Ａ，８６Ａに代えて図１６の閾値設定部８２Ｅ，８６Ｅおよび誤差楕円算出部８４，８８を有する点を除いて、上記状態検出部６２Ａの構成と同じである。

誤差楕円算出部８４は、測位演算結果ｘ_ｐｏ（ｋ−１），Ｐ_ｐｏ（ｋ−１）に基づいて、次式（８０）に示される誤差楕円で定まる領域データを閾値設定部８２Ｆ，８６Ｆに供給する。

ここで、ζは、設定パラメータである。誤差楕円とは、測位結果の誤差範囲を示す確率的な範囲のことであり、測位位置を中心とする楕円状の領域内に目標Ｔｇｔの真の位置があることを意味する。

一方、誤差楕円算出部８８は、追尾演算結果ｘ_ｔｒｋ（ｋ−１），Ｐ_ｔｒｋ（ｋ−１）に基づいて、次式（８１）に示される誤差楕円で定まる領域データを閾値設定部８２Ｆ，８６Ｆに供給する。

ここで、ζ’は、設定パラメータである。

閾値設定部８２Ｆは、式（８０）または（８１）で示される誤差楕円で定まる領域データと衛星軌道情報ＯＤとに基づき、上式（１），（２）に示される閾値ＴＨ１，ＴＨ２の組を設定することができる。また、閾値設定部８６Ｆは、式（８０）または（８１）で示される誤差楕円で定まる領域データと衛星軌道情報ＯＤと目標揺動情報ＳＤとに基づき、上式（５），（６）に示される閾値ＴＨ３，ＴＨ４の組を設定することができる。

以上に説明したように実施の形態６では、誤差楕円で定まる領域データを使用して閾値ＴＨ１，ＴＨ２，ＴＨ３，ＴＨ４が設定されるので、目標Ｔｇｔの状態変化を検知するための信頼性の高い数値範囲Δ（ＴＨ１，ＴＨ２），Δ（ＴＨ３，ＴＨ４）を決定することが可能となる。また、領域データを人為的に設定する必要がないという利点がある。

実施の形態１〜６の変形例．
以上、図面を参照して本発明に係る種々の実施の形態について述べたが、上記実施の形態は本発明の例示であり、これら実施の形態以外の様々な形態を採用することもできる。

たとえば、測位演算部６３Ａ〜６３Ｄのうちから選択された任意の測位演算部と、追尾演算部６４Ａ〜６３Ｃのうちから選択された任意の追尾演算部と、データ移行部６５との組み合わせを有する監視処理部を構成することが可能である。

また、上記実施の形態２〜５の各監視処理部のハードウェア構成は、実施の形態１の場合と同様に、たとえば、ＤＳＰ，ＡＳＩＣまたはＦＰＧＡなどの半導体集積回路を有する単数または複数のプロセッサで実現されればよい。あるいは、監視処理部１１Ａのハードウェア構成は、不揮発性メモリから読み出されたソフトウェアまたはファームウェアのプログラムコードを実行する、ＣＰＵまたはＧＰＵなどの演算装置を含む単数または複数のプロセッサで実現されてもよい。ＤＳＰなどの半導体集積回路とＣＰＵなどの演算装置との組み合わせを含む単数または複数のプロセッサで各監視処理部のハードウェア構成が実現されてもよい。さらには、実施の形態２〜５の各監視処理部のハードウェア構成は、図８に示した信号処理装置７０により実現されてもよい。

本発明の範囲内において、上記実施の形態１〜６の自由な組み合わせ、各実施の形態の任意の構成要素の変形、または各実施の形態の任意の構成要素の省略が可能である。

本発明に係る目標監視装置および目標監視システムは、複数の衛星を用いて、船舶、航空機または車両などの移動自在な電波発信源である目標の位置および移動速度などの状態を高い精度で推定することができるので、たとえば、目標追尾システムおよび捜索救助システムへの使用に適している。

Ｓｔ１〜Ｓｔ３ 人工衛星、Ｒｘ１〜Ｒｘ２ 受信アンテナ、ＳＡ 監視エリア、Ｔｇｔ 目標、１ 目標監視システム、２ 目標監視装置、１０ 受信器、１１Ａ，１１Ｃ，１１Ｅ，１１Ｆ 監視処理部、１２ 表示器、２０ 局部発振器、３１，４１，５１ バンドパスフィルタ、３２，４２，５２ ダウンコンバータ、３３，４３，５３ Ａ／Ｄ変換器（ＡＤＣ）、６１，６１Ｃ 計測部、６２Ａ，６２Ｅ 状態検出部、６３Ａ，６３Ｂ，６３Ｃ，６３Ｄ 測位演算部、６４Ａ，６４Ｂ，６４Ｃ 追尾演算部、６５ データ移行部、６６ 出力制御部、６７ データ記憶部、７０ 信号処理装置、７１ プロセッサ、７２ メモリ、７３ 入力インタフェース部、７４ 出力インタフェース部、７５ 通信回路、７６ 信号路、７８，７９ 遅延素子、８１ 平滑値算出部、８２Ａ，８２Ｂ，８２Ｈ，８２Ｅ，８２Ｆ 閾値設定部、８３ 状態判定部、８４ 誤差楕円算出部、８５ 平滑値算出部、８６Ａ，８６Ｂ，８６Ｈ，８６Ｅ，８６Ｆ 閾値設定部、８７ 状態判定部、８８ 誤差楕円算出部、８９ 判定出力部、９１Ａ 単測位部、９２Ａ 反復推定部、９２Ｂ 反復推定部、９２Ｃ 反復推定部、９２Ｄ バッチ型反復推定部、９３Ａ 単測位部、９４Ａ 反復推定部、９４Ｂ 反復推定部、９４Ｃ 反復推定部、１００ バッファメモリ、１０１ バッファメモリ。

Claims (17)

1. 電波発信源である目標から、複数の衛星を含む３つ以上の通信経路を経て到来した電波信号をそれぞれ受信する複数の受信アンテナで得られた複数の受信信号に基づいて当該目標の状態を推定する目標監視装置であって、
   前記複数の受信信号に基づき、前記受信信号間の到来時間差、または前記受信信号間の到来時間差および到来周波数差を示す複数の計測値を逐次算出する計測部と、
   前記複数の計測値のうちの少なくとも１つの計測値の時間変動を監視して当該目標の動きの有無を検出する状態検出部と、
   前記状態検出部により前記目標の動きが無いことが検出されたときは、前記複数の計測値を用いた測位演算を実行して前記目標の推定位置を示す第１の状態推定値を算出し、前記状態検出部により前記目標の静止状態から移動状態への状態変化が検出されたときは、追尾演算を実行して当該目標の推定位置および移動状態の双方を示す第２の状態推定値を算出する演算部と
   を備え、
   前記状態検出部は、監視エリアを定める座標情報、前記測位演算もしくは前記追尾演算により算出された当該目標の推定位置、または誤差楕円のいずれか１つと、前記複数の衛星それぞれの軌道計算値とに基づいて数値範囲を設定し、
   前記状態検出部は、前記少なくとも１つの計測値にフィルタ処理を施して平滑値を算出し、当該平滑値が、前記数値範囲内から当該数値範囲外へ変動したときに、前記目標の静止状態から移動状態への状態変化が生じたと判定する、
   ことを特徴とする目標監視装置。
2. 請求項１記載の目標監視装置であって、前記状態検出部は、前記平滑値が前記数値範囲外から当該数値範囲内へ変動した後に、前記平滑値が当該数値範囲内に一定期間の間、継続して収まるときに、前記目標の移動状態から静止状態への状態変化が生じたと判定することを特徴とする目標監視装置。
3. 請求項１記載の目標監視装置であって、前記演算部は、前記状態検出部により前記目標の静止状態から移動状態への状態変化が検出されたときは、前記測位演算で得られた直近の演算結果を初期値として用い、かつ前記計測部で算出された複数の計測値を用いて前記追尾演算を実行することを特徴とする目標監視装置。
4. 請求項２記載の目標監視装置であって、前記演算部は、前記状態検出部により前記目標の移動状態から静止状態への状態変化が検出されたときは、前記追尾演算により得られた直近の演算結果を初期値として用い、かつ前記計測部で算出された複数の計測値を用いて前記測位演算を実行することを特徴とする目標監視装置。
5. 請求項１記載の目標監視装置であって、
   前記演算部は、
   前記複数の計測値の組を観測ベクトルとして用いた非線形最小二乗法による演算を実行して測位ベクトルを算出する単測位部と、
   前記測位ベクトルを用いたカルマンフィルタによる反復推定処理を実行して前記第１の状態推定値を算出する反復推定部と
   を含むことを特徴とする目標監視装置。
6. 請求項１記載の目標監視装置であって、
   前記複数の計測値は、時刻毎に互いに同期して、または互いにほぼ同期して得られる観測値であり、
   前記演算部は、前記複数の計測値の組を観測ベクトルとして用いたカルマンフィルタによる反復推定処理を前記測位演算として実行することを特徴とする目標監視装置。
7. 請求項１記載の目標監視装置であって、
   前記複数の計測値は、互いに非同期に得られる観測値であり、
   前記演算部は、前記複数の計測値を用いたカルマンフィルタによる反復推定処理を前記測位演算として実行することを特徴とする目標監視装置。
8. 請求項１記載の目標監視装置であって、
   前記複数の計測値は、互いに非同期に得られる観測値であり、
   前記演算部は、前記複数の計測値の複数回分を一括して用いた非線形最小二乗法によるバッチ処理を前記測位演算として実行することを特徴とする目標監視装置。
9. 請求項５記載の目標監視装置であって、前記カルマンフィルタによる反復推定処理は、前記測位ベクトルが外れ値か否かを判定する外れ値判定処理を含むことを特徴とする目標監視装置。
10. 請求項６記載の目標監視装置であって、前記カルマンフィルタによる反復推定処理は、前記観測ベクトルが外れ値か否かを判定する外れ値判定処理を含むことを特徴とする目標監視装置。
11. 請求項１記載の目標監視装置であって、
    前記演算部は、
    前記複数の計測値の組を観測ベクトルとして用いた非線形最小二乗法による演算を実行して測位ベクトルを算出する単測位部と、
    前記測位ベクトルを用いたカルマンフィルタによる反復推定処理を実行して前記第２の状態推定値を算出する反復推定部と
    を含むことを特徴とする目標監視装置。
12. 請求項１記載の目標監視装置であって、
    前記複数の計測値は、時刻毎に互いに同期して、または互いにほぼ同期して得られる観測値であり、
    前記演算部は、前記複数の計測値を用いたカルマンフィルタによる反復推定処理を前記追尾演算として実行することを特徴とする目標監視装置。
13. 請求項１記載の目標監視装置であって、
    前記複数の計測値は、互いに非同期に得られる観測値であり、
    前記演算部は、前記複数の計測値を用いたカルマンフィルタによる反復推定処理を前記追尾演算として実行することを特徴とする目標監視装置。
14. 請求項１１記載の目標監視装置であって、前記カルマンフィルタによる反復推定処理は、前記測位ベクトルが外れ値か否かを判定する外れ値判定処理を含むことを特徴とする目標監視装置。
15. 請求項１２記載の目標監視装置であって、前記カルマンフィルタによる反復推定処理は、前記複数の計測値が外れ値か否かを判定する外れ値判定処理を含むことを特徴とする目標監視装置。
16. 請求項１３記載の目標監視装置であって、前記カルマンフィルタによる反復推定処理は、前記複数の計測値が外れ値か否かを判定する外れ値判定処理を含むことを特徴とする目標監視装置。
17. 請求項１から請求項１６のうちのいずれか１項記載の目標監視装置と、
    前記複数の受信信号を出力する複数の受信アンテナと
    を備えることを特徴とする目標監視システム。

Images