(58)【調査した分野】(Int.Cl.,DB名)
【発明を実施するための形態】
【0010】
以下、図面を参照しつつ、本発明に係る種々の実施の形態について詳細に説明する。なお、図面全体において同一符号を付された構成要素は、同一構成及び同一機能を有するものとする。
【0011】
実施の形態1.
図1は、本発明に係る実施の形態1である目標監視システム1の構成を概略的に示す図である。
図1に示されるように目標監視システム1は、監視エリアSA内の電波発信源である目標Tgtから複数の通信経路を経て到来した電波信号をそれぞれ受信して受信信号S1,S2,S3を出力する受信アンテナ(受信局)Rx1,Rx2,Rx3と、受信信号S1,S2,S3に基づいて目標Tgtの状態を推定する目標監視装置2とを備えて構成されている。本実施の形態では、受信アンテナRx1,Rx2,Rx3と目標Tgtとの間の3本の通信経路にそれぞれ3基の人工衛星St1,St2,St3(以下、単に「衛星St1,St2,St3」という。)が存在する。衛星St1,St2,St3に搭載されたトランスポンダ(中継器)は、目標Tgtから発信された電波信号を受信し、当該受信された電波信号を増幅してダウンリンク信号を生成し、これらダウンリンク信号を目標監視システム1の受信アンテナRx1,Rx2,Rx3にそれぞれ送信する。衛星St1,St2,St2としては、静止軌道に投入された静止衛星が想定されている。なお、3基の衛星St1,St2,St3に限定されずに2基またはN基の衛星(Nは4以上の整数)からそれぞれダウンリンク信号を受信可能なように本実施の形態の目標監視システム1の構成が変更されてもよい。
【0012】
目標Tgtは、移動自在な電波発信源(船舶、航空機または車両などの移動体)であり、特定の周波数帯域の電波信号を衛星に向けて発信する機能を有する。しかしながら、目標Tgtから発信された電波信号が、地上の無線通信網または衛星通信網の送信電波と干渉して通信障害を引き起こすことがある。このような場合に、干渉波となる電波信号を発信する目標Tgtの位置を特定したいというニーズがある。また、目標Tgtから発信された捜索救助用のビーコン信号に基づいて目標Tgtの位置を特定したいというニーズもある。本実施の形態の目標監視装置2は、受信信号S1〜S3に基づいて、目標Tgtの動きの有無、すなわち目標Tgtの状態が静止状態または移動状態のいずれであるかを検出することができる。目標Tgtが静止状態にあるとき、目標監視装置2は、目標Tgtの推定位置(静止位置)を示す状態推定値を逐次算出することができ、目標Tgtが移動状態にあるときは、目標監視装置2は、目標Tgtの推定位置(追尾位置)および移動状態(たとえば速度)を示す状態推定値を逐次算出することができる。
【0013】
目標監視装置2は、
図1に示されるように、高周波帯域の受信信号S1,S2,S3をディジタル受信信号D1,D2,D3に変換する受信器10と、ディジタル受信信号D1,D2,D3に基づいて目標Tgtの状態を示す状態推定値を逐次算出する監視処理部11Aと、当該算出された状態推定値に基づく画像情報を表示する表示器12とを備える。
【0014】
図2は、実施の形態1における受信器10の構成例を示すブロック図である。
図2に示される受信器10は、周波数変換用の局部発振信号を出力する局部発振器20と、高周波帯域の受信信号S1,S2,S3にそれぞれフィルタ処理を施すバンドパスフィルタ31,41,51と、局部発振信号を用いてバンドパスフィルタ31,41,51の高周波出力F1,F2,F3にそれぞれ周波数変換を施すダウンコンバータ32,42,52と、A/D変換器(ADC)33,43,53とを有している。ダウンコンバータ32,42,52は、高周波出力F1,F2,F3を局部発振信号と混合して、より低い周波数帯域のアナログ信号C1,C2,C3を生成する周波数変換器である。ADC33,43,53は、アナログ信号C1,C2,C3をそれぞれディジタル受信信号D1,D2,D3に変換し、当該ディジタル受信信号D1,D2,D3を監視処理部11Aに出力する。ディジタル受信信号D1,D2,D3の各々は、振幅成分と位相成分とを含む複素信号である。
【0015】
図3は、実施の形態1における監視処理部11Aの概略構成を示すブロック図である。
図3に示される監視処理部11Aは、ディジタル受信信号D1,D2,D3に基づいて電波の到来時間差および到来周波数差の計測値を逐次算出する計測部61と、これら計測値のいずれかに基づいて目標Tgtの動きの有無を検出する状態検出部62Aと、目標Tgtの動きが無いことが検出されたときは、当該計測値を用いた測位演算を実行して目標Tgtの推定位置を示す状態推定値を逐次算出する測位演算部63Aと、目標Tgtの動きが有ることが検出されたきは、当該計測値を用いた追尾演算を実行して目標Tgtの推定位置および移動状態(たとえば推定速度)の双方を示す状態推定値を逐次算出する追尾演算部64Aと、測位演算部63Aおよび追尾演算部64Aの一方から他方へ演算結果のデータを移行させるデータ移行部65と、状態検出部62Aでの処理に必要なデータ(衛星軌道情報OD,監視エリア情報ADおよび目標揺動情報SD)を格納するデータ記憶部67と、外部サーバ(図示せず。)から当該データを取得するデータ取得部68とを有する。ここで、kは時刻t
kを示す整数である。
【0016】
また、監視処理部11Aは、測位演算部63Aおよび追尾演算部64Aで算出された状態推定値を基に画像情報DDを生成する出力制御部66を備える。出力制御部66は、当該画像情報DDを液晶ディスプレイまたは有機ELディスプレイなどの表示器12に出力する。画像情報DDとしては、たとえば、目標Tgtの推定位置の座標値を表す英数字情報、その推定位置を視覚的に表す地図情報、目標Tgtの推定速度の値を示す英数字情報、および、その推定速度の遷移を視覚的に表すグラフが挙げられる。ユーザは、表示器12に表示された画像情報DDを視ることで目標Tgtの状態を把握することができる。
【0017】
計測部61は、ディジタル受信信号D1,D2,D3に基づき、受信信号S1,S2間の到来時間差TDOA
12(k)、受信信号S2,S3間の到来時間差TDOA
23(k)、受信信号S3,S1間の到来時間差TDOA
31(k)、受信信号S1,S2間の到来周波数差FDOA
12(k)、受信信号S2,S3間の到来周波数差FDOA
23(k)、および、受信信号S3,S1間の到来周波数差FDOA
31(k)の計測値を逐次算出する機能を有する。計測値を算出する具体的な方法としては、公知の方法が採用されればよく、特に限定されるものではない。たとえば、時間方向および周波数方向についてディジタル受信信号Di,Dj間の2次元の相互相関を算出し、当該算出された相互相関に現れるピークを検出し、当該ピークに対応する到来時間差TDOA
ij(k)および到来周波数FDOA
ij(k)の組を得るという算出法を採用することができる。ここで、i,jは、1〜3の範囲内の整数である(i≠j)。
【0018】
また、計測部61は、到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)および到来周波数差FDOA
12(k),FDOA
23(k),FDOA
31(k)の計測値を、状態検出部62A、測位演算部63Aおよび追尾演算部64Aに供給する。
【0019】
状態検出部62Aは、計測部61から逐次供給される計測値のうちの少なくとも1つを監視対象とし、この監視対象の計測値の時間変動を常時監視して目標Tgtの動きの有無を検出することができ、その検出結果を示す検出信号ESを測位演算部63A,追尾演算部64Aおよびデータ移行部65に供給する。
【0020】
図4は、実施の形態1の状態検出部62Aの構成例を概略的に示すブロック図である。最初に、到来時間差TDOA
ij(k)を監視対象として選択した場合の状態検出部62Aの構成および動作について説明する。
【0021】
図4に示される状態検出部62Aは、監視対象の計測値の時系列データにフィルタ処理を施して平滑値<TDOA>
kを算出する平滑値算出部81と、閾値TH1,TH2の組を設定する閾値設定部82Aと、平滑値<TDOA>
kを、閾値TH1を下限としかつ閾値TH2を上限とする数値範囲Δ(TH1,TH2)と比較して目標Tgtの状態が静止状態または移動状態のいずれであるかを判定する状態判定部83とを有する。ここで、kは時刻t
kを示す整数である。状態判定部83は、目標Tgtの状態が静止状態または移動状態のいずれであるかを示す判定結果ES1を判定出力部89に与える。たとえば、N個の到来時間差TDOA
ij(k−N+1),TDOA
ij(k−N+2),…,TDOA
ij(k−1),TDOA
ij(k)の計測値にフィルタ処理(たとえば、平均化)を施すことで、時刻t
kでの平滑値<TDOA>
kを算出することが可能である。到来時間差TDOA
ij(k)の値には、計測誤差に起因するバラツキが存在する。平滑値<TDOA>
kを使用することにより、そのようなバラツキの影響を低減させることが可能となる。
【0022】
状態判定部83は、平滑値<TDOA>
kの時間変動を常時監視し、平滑値<TDOA>
kが一定期間T1の間、数値範囲Δ(TH1,TH2)内に継続して収まるときに、目標Tgtは静止状態にあると判定する。
図5Aは、平滑値<TDOA>
kの時間変動の一例を表すグラフである。
図5Aのグラフに示されるように、目標Tgtが静止状態にあるとき(t<t
mのとき)、平滑値<TDOA>
kは、閾値TH1,TH2間の数値範囲内に収まるように変動する。
【0023】
平滑値<TDOA>
kが数値範囲Δ(TH1,TH2)内から数値範囲Δ(TH1,TH2)外へ変動したとき、状態判定部83は、目標Tgtの静止状態から移動状態への状態変化が生じたと判定する。
図5Aの例では、平滑値<TDOA>
kは、時刻t
mで上限の閾値TH2を超えるように変動しているため、状態判定部83は、時刻t
mで目標Tgtの状態変化が生じたと判定することができる。
【0024】
閾値TH1,TH2については、目標Tgtの静止状態から移動状態への状態変化による平滑値<TDOA>
kの変動を速やかに検知することができるように数値範囲Δ(TH1,TH2)の幅を極力小さくする一方で、その状態変化以外の要因による平滑値<TDOA>
kの変動を検知しないように数値範囲Δ(TH1,TH2)の幅を大きくすることが望ましい。その状態変化以外の要因としては、主に、衛星Sti,Stjの動きが挙げられる。したがって、目標Tgtの状態変化による平滑値<TDOA>
kの変動と、衛星Sti,Stjの動きによる平滑値<TDOA>
kの変動とを考慮して、適切な閾値TH1,TH2を設定することが望ましい。
【0025】
閾値設定部82Aは、データ記憶部67から供給された衛星軌道情報ODおよび監視エリア情報ADに基づき、次の状態判定式(1),(2)に示される閾値TH1,TH2の組を設定することができる。
【0026】
ここで、TDOA
minは、i番目の衛星Stiおよびj番目の衛星Stj(i≠j)の動きを考慮して設定された最小値であり、TDOA
maxは、衛星Sti,Stjの動きを考慮して設定された最大値である。また、σ
TDOAは、到来時間差の計測誤差の標準偏差である。3σ
TDOAは、到来時間差の計測誤差を考慮して標準偏差の3倍相当のマージンを与えるものである。
【0027】
状態判定部83は、状態判定式(1),(2)のいずれか一方が満たされたことを検知したときに、目標Tgtの静止状態から移動状態への状態変化が生じたと判定することができる。最小値TDOA
minおよび最大値TDOA
maxは、衛星軌道情報ODおよび監視エリア情報ADを用いて算出可能である。監視エリア情報ADは、
図1に示した監視エリアSAを定める座標情報である。衛星軌道情報ODは、衛星Sti,Stjの軌道計算に必要な軌道要素情報、または、軌道要素情報を用いた軌道計算アルゴリズムにより算出された軌道計算値であればよい。
【0028】
たとえば、閾値設定部82Aは、あらかじめ、監視エリア情報ADで定められた監視エリアSA内の複数の地点Ψ
1,Ψ
2,…,Ψ
L(Lは正整数)に電波発信源が存在するとの仮定の下で、衛星Sti,Stjの最新の軌道計算値を用いて地点Ψ
1,Ψ
2,…,Ψ
Lにそれぞれ対応する到来時間差TDOA(Ψ
1,k),TDOA
k(Ψ
2,k),…,TDOA(Ψ
L,k)を予測する(kは時刻t
kを示す整数)。閾値設定部82Aは、予測された到来時間差TDOA(Ψ
1,k)〜TDOA(Ψ
L,k)のうちの最小値および最大値をそれぞれ最小値TDOA
minおよび最大値TDOA
maxとして決定することができる。衛星Sti,Stjの軌道は24時間周期を持つため、閾値設定部82Aは、最新の軌道計算値を用いて24時間分の最小値TDOA
min(k)および最大値TDOA
max(k)をあらかじめ計算すれば十分である。また、閾値設定部82Aは、予測された到来時間差TDOA(Ψ
1,k)〜TDOA(Ψ
L,k)の中から、平滑値<TDOA>
kに近い到来時間差を選択し、当該選択された到来時間差に基づいて最小値TDOA
minおよび最大値TDOA
maxを決定してもよい。
【0029】
次に、目標Tgtの移動状態から静止状態への状態変化が生じた場合、平滑値<TDOA>
kの変動幅が一定範囲内に収まる。この場合、状態判定部83は、平滑値<TDOA>
kが一定期間T1の間、数値範囲Δ(TH1,TH2)内に継続して収まることを検知したときに、目標Tgtの移動状態から静止状態への状態変化が生じたと判定する。一定期間T1の時間長は、パラメータ設定値である。
図5Bは、平滑値<TDOA>
kの時間変動の他の例を表すグラフである。
図5Bのグラフに示されるように、平滑値<TDOA>
kの変動幅が時刻t
d以後小さくなり、時刻t
d〜時刻t
sの間、平滑値<TDOA>
kは、数値範囲Δ(TH1,TH2)内に継続して収まるように推移している。このとき、状態判定部83は、時刻t
sで目標Tgtの移動状態から静止状態への状態変化が生じたと判定することができる。
【0030】
より具体的には、状態判定部83は、一定期間T1の間、次の状態判定式(3),(4)の双方が満たされることを継続して検知したときに、目標Tgtの移動状態から静止状態への状態変化が生じたと判定することができる。
【0032】
次に、
図4を参照しつつ、到来周波数差FDOA
ij(k)を監視対象として選択した場合の状態検出部62Aの構成および動作について説明する。
【0033】
図4に示される状態検出部62Aは、監視対象である計測値の時系列データにフィルタ処理を施して平滑値<FDOA>
kを算出する平滑値算出部85と、閾値TH3,TH4の組を設定する閾値設定部86Aと、平滑値<FDOA>
kを、閾値TH3を下限としかつ閾値TH4を上限とする数値範囲Δ(TH3,TH4)と比較して目標Tgtの状態が静止状態または移動状態のいずれであるかを判定する状態判定部87とを有する。ここで、kは時刻t
kを示す整数である。たとえば、N個の到来周波数差FDOA
ij(k−N+1),FDOA
ij(k−N+2),…,FDOA
ij(k−1),FDOA
ij(k)の計測値にフィルタ処理(たとえば、平均化)を施すことで、時刻t
kでの平滑値<FDOA>
kを算出することが可能である。到来周波数差FDOA
ij(k)の値には、計測誤差に起因するバラツキが存在する。平滑値<FDOA>
kを使用することにより、そのようなバラツキの影響を低減させることが可能となる。
【0034】
状態判定部87は、平滑値<FDOA>
kの時間変動を常時監視し、平滑値<FDOA>
kが一定期間T2の間、数値範囲Δ(TH3,TH4)内に継続して収まるときに、目標Tgtは静止状態にあると判定することができる。また、平滑値<FDOA>
kが数値範囲Δ(TH3,TH4)内から数値範囲Δ(TH3,TH4)外へ変動したとき、状態判定部87は、目標Tgtの静止状態から移動状態への状態変化が生じたと判定することができる。
【0035】
閾値TH3,TH4については、目標Tgtの静止状態から移動状態への状態変化による平滑値<FDOA>
kの変動を速やかに検知することができるように数値範囲Δ(TH3,TH4)の幅を極力小さくする一方で、その状態変化以外の要因による平滑値<FDOA>
kの変動を検知しないように数値範囲Δ(TH3,TH4)の幅を大きくすることが望ましい。その状態変化以外の要因としては、主に、衛星Sti,Stjの動きと、目標Tgtの姿勢変化などの動揺によるドップラー効果とが挙げられる。たとえば、目標Tgtが船舶の場合、波または風の影響による船舶の動揺(たとえば、縦揺れおよび横揺れ)が生じ、これに起因するドップラー効果が発生することがある。したがって、目標Tgtの状態変化による平滑値<FDOA>
kの変動と、人工衛星Sti,Stjの動きによる平滑値<FDOA>
kの変動と、目標Tgtの動揺による平滑値<FDOA>
kの変動とを考慮して、適切な閾値TH3,TH4を設定することが望ましい。
【0036】
閾値設定部86Aは、データ記憶部67から供給された衛星軌道情報OD,監視エリア情報ADおよび目標揺動情報SDに基づき、次の状態判定式(5),(6)に示される閾値TH3,TH4の組を設定することができる。
【0037】
ここで、FDOA
minは、i番目の衛星Stiおよびj番目の衛星Stj(i≠j)の動きを考慮して設定された最小値であり、FDOA
maxは、衛星Sti,Stjの動きを考慮して設定された最大値である。また、σ
FDOAは、到来周波数差の計測誤差の標準偏差である。3σ
FDOAは、到来周波数差の計測誤差を考慮して標準偏差の3倍相当のマージンを与えるものである。さらに、σ
motionは、目標揺動情報SDに基づき、目標Tgtの動揺によるドップラー効果を考慮して設定された値であり、当該ドップラー効果に起因する到来周波数差の揺らぎを示す値である。たとえば、目標Tgtが船舶の場合、海洋上の多数の地点における波に関するデータが提供されている。データ取得部68は、外部サーバ(図示せず)から、このような波に関するデータを目標揺動情報SDとして取得することができるので、閾値設定部86Aは、目標揺動情報SDに基づき、任意の地点における船舶の動揺の見積もりを実行して到来周波数差の揺らぎを示す値σ
motionを算出することができる。
【0038】
閾値設定部86Aは、状態判定式(5),(6)のいずれか一方が満たされたことを検知したときに、目標Tgtの静止状態から移動状態への状態変化が生じたと判定することができる。最小値TDOA
minおよび最大値TDOA
maxを算出する場合と同様に、最小値FDOA
minおよび最大値FDOA
maxは、衛星軌道情報ODおよび監視エリア情報ADを用いて算出可能である。また、揺らぎσ
motionは、目標Tgtの位置と衛星Sti,Stjの位置とが不変であるとの仮定の下、目標Tgtの動揺速度の分散値を含む目標揺動情報SDに基づいて算出可能である。目標Tgtの動揺速度の分散値が与えられれば、マージン3σ
motionを算出することができる。
【0039】
一方、目標Tgtの移動状態から静止状態への状態変化が生じる場合、平滑値<FDOA>
kの変動幅が一定範囲内に収まる。この場合、状態判定部87は、平滑値<FDOA>
kが一定期間T2の間、数値範囲Δ(TH3,TH4)内に継続して収まることを検知したときに、目標Tgtの移動状態から静止状態への状態変化が生じたと判定する。一定期間T2の時間長は、パラメータ設定値である。より具体的には、状態判定部87は、一定期間T2の間、次の状態判定式(7),(8)の双方が満たされることを継続して検知したときに、目標Tgtの移動状態から静止状態への状態変化が生じたと判定することができる。
【0041】
以上に説明したように、状態判定部83は、目標Tgtの動きの有無、すなわち、目標Tgtの状態が静止状態または移動状態のいずれであるかを判定し、その判定結果ES1を判定出力部89に与える。判定出力部89は、
図3に示されるように、判定結果ES1を示す検出信号ESを測位演算部63A,追尾演算部64Aおよびデータ移行部65に供給する。
【0042】
次に、測位演算部63A、追尾演算部64Aおよび追尾演算部64Aを詳細に説明する前に、実施の形態1および後述の実施の形態2〜4で使用されるカルマンフィルタおよび非線形最小二乗法について説明する。カルマンフィルタは、フィルタリング理論の一種である。
【0043】
カルマンフィルタの目的は、時刻t
kにおいて雑音の混入した観測値系列z
1,…,z
k−1,z
kが得られたとき、当該観測値系列と状態空間モデルとを用いて目標Tgtの静止状態および移動状態を示す物理量を推定することである。状態空間モデルは、目標Tgtの運動モデル(状態方程式)と観測モデル(観測方程式)とで構成される。運動モデルは次式(9)で与えられる。
【0044】
ここで、x
kは、時刻t
kにおけるn行1列の状態ベクトルであり、観測不能な真の状態を表している。nは2以上の整数である。w
kは、時刻t
kにおけるn行1列の雑音ベクトルであり、F
kは、時刻t
kにおけるn行n列の状態遷移行列である。後述するように、状態ベクトルx
kの成分は、測位演算に適用される場合と追尾演算で適用される場合とで異なる。雑音ベクトルw
kの平均は零ベクトルとし、雑音ベクトルw
kの共分散行列はQ
kで表現されるものとする。
【0045】
一方、観測モデルは次式(10)で与えられる。
【0046】
ここで、z
kは、時刻t
kにおけるm行1列の観測ベクトルであり、v
kは、時刻t
kにおけるm行1列の雑音ベクトルであり、mは正整数である。雑音ベクトルv
kの平均は零ベクトルとし、雑音ベクトルv
kの共分散行列はR
kで表現されるものとする。また、h[x
k]は、状態ベクトルx
kに関する観測関数であり、m行1列のベクトルである。m=1のとき、z
k,v
k,h[x
k]は、スカラー量となる。
【0047】
カルマンフィルタによる推定処理(反復推定処理)は、予測処理、平滑化処理(更新処理)および外れ値判定処理を反復して実行することにより実現される。以下、x(k|p)は、時刻t
pの時点での時刻t
kの状態推定値からなるn行1列のベクトルを表すものとする。
【0048】
予測処理では、次式(11),(12)に従い、事前状態推定ベクトルx(k|k−1)および事前誤差共分散行列P(k|k−1)が計算される。
【0049】
ここで、F
kTは、状態遷移行列F
kに対する転置行列である。
【0050】
平滑化処理(更新処理)では、次式(13),(14)に従い、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)が算出される。
【0051】
ここで、e
kは、観測値と状態推定値との間の残差からなる観測残差ベクトルである。また、S
kは観測残差共分散行列、K
kはカルマンゲイン、H
kは観測行列である。観測残差ベクトルe
k、観測残差共分散行列S
kおよびカルマンゲインK
kは、それぞれ、次式(15),(16),(17)で与えられる。
【0053】
観測行列H
kは、次式(18)で示されるヤコビアン行列(Jacobian matrix)によって定義される。
【0055】
また、上記予測処理の後に、観測値が外れ値(outlier values)か否かを判定する外れ値判定処理が実行される。外れ値判定処理では、まず、次式(19)に基づいて、次式(20)で定義される観測残差共分散行列S
kを用いた重み付き残差平方和D
kが算出される。
【0057】
式(19)は、上式(15)で示される観測残差ベクトルe
kを用いれば、次式(21)として表現可能である。
【0059】
その後、次の判定式(22)が成立するか否かが判定される。
【0060】
ここで、γは、カイ二乗分布に従うゲートサイズパラメータである。
【0061】
判定式(22)が成立する場合には、観測値は外れ値ではないと判定され、平滑化処理(更新処理)が実行される。一方、判定式(22)が成立しない場合には、観測値は外れ値であると判定されるので、平滑化処理は実行されない。この場合、上式(13),(14)の代わりに次式(23)に従って、状態推定値x(k|k)および事後誤差共分散行列P(k|k)が算出される。
【0063】
図6は、カルマンフィルタによる反復推定処理の手順の一例を概略的に示すフローチャートである。まず、ステップST11では、初期化が実行される。すなわち、時刻t
0(k=0)における状態推定ベクトルx(0|0)および事後誤差共分散行列P(0|0)が与えられる。
【0064】
次いで、時刻t
1(k=1)における観測値の組すなわち観測ベクトルz
kが得られるまで待機し(ステップST12のNO)。観測ベクトルz
kが得られると(ステップST12のYES)、時刻t
k(k=1)について上記予測処理が実行される(ステップST13)。この結果、事前状態推定ベクトルx(k|k−1)および事前誤差共分散行列P(k|k−1)が得られる。続けて、上記外れ値判定処理が実行される(ステップST14)。その結果、観測ベクトルz
kにおける観測値が外れ値ではないと判定されたときは(ステップST15のNO)、上記平滑化処理が実行される(ステップST20)。これにより、時刻t
k(k=1)における状態推定値x(k|k)および事後誤差共分散行列P(k|k)が計算される。
【0065】
平滑化処理(ステップST20)の後、収束条件が満たされるか否かが判定される(ステップST21)。収束条件は、各時刻t
kにおける状態推定値x(k|k)および事後誤差共分散行列P(k|k)が十分に収束していると推定される場合に満たされる条件である。収束条件が満たされない場合(ステップST21のNO)、事前状態推定ベクトルx(k|k−1)および事前誤差共分散行列P(k|k−1)が、ステップST20で算出された状態推定値x(k|k)および事後誤差共分散行列P(k|k)でそれぞれ置き換えられる(ステップST22)。次いで、平滑化処理(ステップST20)が再度実行される。時刻t
kについて、最終的に収束条件が満たされた場合に(ステップST21のYES)、次のステップS24が実行される。このように収束条件が満たされるまで、平滑化処理を繰り返し実行すること(巡回処理)によって、精度の高い状態推定値x(k|k)および事後誤差共分散行列P(k|k)を得ることができる。
【0066】
ステップST21の収束条件は、たとえば、各時刻t
kについて平滑化処理があらかじめ設定された回数だけ繰り返し実行された場合、または、各時刻t
kについて、直近のj回目に算出された状態推定値x
j(k|k)とj−1回目に算出された状態推定値x
j−1(k|k)との間のノルムが閾値以下となった場合に満たされるものとすることができる(jは、巡回処理における平滑化処理の繰り返し回数)。あるいは、直近のj回目に算出された事後誤差共分散行列P
j(k|k)とj−1回目に算出された事後誤差共分散行列P
j−1(k|k)との間のノルムが閾値以下となり、かつ、状態推定値x
j(k|k),x
j−1(k|k)間のノルムが閾値以下となった場合にステップST21の収束条件が満たされてもよい。
【0067】
一方、観測値が外れ値であると判定されたときは(ステップST15のYES)、ステップST20の平滑化処理は実行されず、上式(23)に従って、状態推定値x(k|k)および事前誤差共分散行列P(k|k)が算出される(ステップST16)。
【0068】
ステップST16またはステップST21の後、あらかじめ定められた反復条件が満たされると判定されたとき(ステップST24のYES)、次の時刻t
k(k=2)について上記処理手順が反復して実行される(ステップST25,ST12以後)。最終的に、反復条件が満たされないと判定されたときに(ステップST24のNO)、反復推定処理は終了する。
【0069】
次に、非線形最小二乗法について説明する。
【0070】
まず、次式(24)の観測モデル(観測方程式)が成立すると仮定する。
【0071】
ここで、y
kは、時刻t
kにおけるp行1列の観測ベクトルであり、η
kは、時刻t
kにおけるp行1列の雑音ベクトルであり、pは2以上の整数である。雑音ベクトルη
kの平均は零ベクトルとし、雑音ベクトルη
kの共分散行列(観測雑音共分散行列)はΩ
kで表現されるものとする。また、φ
k[χ]は、状態ベクトルχに関する観測関数であり、p行1列のベクトルである。
【0072】
K回の計測結果(Kは2以上の整数)に対して、観測残差の大きさを評価する評価関数J(χ)が次式(25)により定義される。
【0074】
評価関数J(χ)は、雑音環境下での観測ベクトルy
kとベクトルφ
k[χ]との間の距離(「マハラノビス距離(Mahalanobis distance)」と呼ばれる。)の二乗の総和を表している。
【0075】
非線形最小二乗法の目的は、次式(26)に示されるように評価関数J(χ)を最小にするχの近似解χ
minを求めることである。
【0077】
仮に、観測関数φ
k[χ]が、次式(27)に示されるように線形でかつ行列Γ
kを用いて表現可能であるとする。
【0079】
このとき、非線形最小二乗法による近似解χ
minは、次式(28)で与えられる。
【0081】
一方、観測関数φ
k[χ]が非線形である場合には、再帰最小二乗法(Recursive Least−Squares,RLS)アルゴリズムにより、近似解χ
minを求めることができる。たとえば、観測関数φ
k[χ]が、近似的に、次式(29)で示されるヤコビアン行列(関数行列)を用いて表現可能であるとする。
【0083】
RLSアルゴリズムでは、たとえば、次の収束演算式(30),(31),(32)を使用することができる(mは、再帰処理における繰り返し回数)。
【0085】
RLSアルゴリズムに基づく再帰処理の手順は以下のとおりである。まず、初期化が実行される。すなわち、解の候補の初期値χ
0(m=0)が設定される。次に、繰り返し回数m=1,2,…について、上式(30),(31),(32)に基づいた演算が反復して実行される。このとき、次式(33)で示される反復終了条件が満たされたときに再帰処理が正常に終了する。
【0086】
ここで、εは、微少な正の実数であり、Mは、繰り返し回数の上限値である。最終的に、反復終了条件を満たす解χ
mが、近似解χ
minとして与えられる。
【0087】
次に、測位演算部63Aの構成および動作について説明する。
【0088】
測位演算部63Aは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)の計測値の中から2以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた測位演算を実行することにより目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。測位演算部63Aは、互いに同期して得られる到来時間差TDOA
12(k),TDOA
23(k)の計測値の組を選択し、この組を観測ベクトルとして用いた測位演算を実行することができる。
【0089】
あるいは、測位演算部63Aは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)および到来周波数差FDOA
12(k),FDOA
23(k),FDOA
31(k)の計測値の中から、少なくとも1つの到来時間差および少なくとも1つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた測位演算を実行することにより目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を逐次算出する機能を有している。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、測位演算部63Aは、到来時間差TDOA
12(k)と到来周波数差FDOA
12(k)の計測値の組を選択し、この組を観測ベクトルとして用いた測位演算を実行することができる。
【0090】
図7は、実施の形態1における測位演算部63A、追尾演算部64Aおよびデータ移行部65からなる演算部の構成を概略的に示すブロック図である。
【0091】
図7に示されるように、測位演算部63Aは、単測位部91Aと反復推定部92Aとを有する。単測位部91Aは、当該選択された計測値の組を観測ベクトルとして用いた非線形最小二乗法による測位演算を実行して目標Tgtの位置を示す測位ベクトルPOS
kを算出し、当該測位ベクトルPOS
kを反復推定部92Aに出力する。反復推定部92Aは、測位ベクトルPOS
kを用いたカルマンフィルタによる反復推定処理を実行して目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を算出する。
【0092】
より具体的には、単測位部91Aは、上記した再帰最小二乗法(RLS)アルゴリズムに従い、当該選択された計測値の組を観測ベクトルy
kとして用いた再帰処理を実行することができる。たとえば、TDOA
12(k),TDOA
23(k)の計測値の組が観測ベクトルy
kとして選択されている場合、観測ベクトルy
kは、次式(34)で与えられる。
【0094】
このとき、時刻t
kにおける観測関数φ
k[χ=p
tgt]は、次式(35)で与えられる。
【0095】
ここで、p
tgtは、地球の重心を原点とする、目標Tgtの位置座標を示す位置ベクトルであり、p
1(k)は、時刻t
kでの衛星St1の位置座標を示す位置ベクトルであり、p
2(k)は、時刻t
kでの衛星St2の位置座標を示す位置ベクトルである。
【0096】
観測雑音共分散行列Ω
kは、次式(36)で与えられるものとする。
【0098】
単測位部91Aは、複数回の計測結果に対して上記したRLSアルゴリズムによる再帰処理を実行することにより、評価関数J(χ=p
tgt)を最小にする近似解χ
minを測位ベクトルPOS
kとして算出することができる。単測位部91Aは、測位ベクトルPOS
kを反復推定部92Aに供給する。
【0099】
単測位演算における理論的な測位誤差を表す誤差共分散行列Λは、次式(36A)で与えられる。
ここで、Φ
k[χ]は、上式(29)で示されるヤコビアン行列である。
【0100】
次に、反復推定部92Aは、上式(9),(10)で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を実行する。反復推定部92Aに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(37),(38),(39)で与えられる。
【0101】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度である。目標Tgtは静止状態にあるので、状態遷移行列F
kは単位行列である。また、目標Tgtの運動に不確かさがなく、雑音ベクトルw
kは零ベクトルとなる。よって、共分散行列Q
kは零行列である。
【0102】
また、反復推定部92Aに適用される観測モデル(観測方程式)では、観測ベクトルz
k、観測行列H
kおよび共分散行列R
kは、たとえば、次式(40),(41),(42)で与えられる。
【0103】
ここで、式(40)におけるOLN
k,OLT
kは、目標Tgtの観測位置を示す経度および緯度であり、式(42)におけるDOPは、精度劣化度(Dilution of Precision)と呼ばれ、本実施の形態では、単測位演算における理論的な測位誤差式によって与えられる。具体的には、DOPは、上式(36A)の誤差共分散行列Λに設定されればよい。
【0104】
反復推定部92Aは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部92Aは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)として出力する。上記のとおり、単測位部91Aは、複数回の計測結果(累積計測結果)を用いた非線形最小二乗法による測位演算を実行して高精度な測位ベクトルPOS
kを算出することができるので、反復推定部92Aは、高精度な測位ベクトルPOS
kに基づき、精度の高い測位演算結果を生成することができる。
【0105】
次に、追尾演算部64Aの構成および動作について説明する。
【0106】
追尾演算部64Aは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)の計測値の中から2以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Tgtの推定位置および移動状態(たとえば推定速度)を示す状態推定ベクトルx
trk(k)を逐次算出する機能を有する。あるいは、追尾演算部64Aは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)および到来周波数差FDOA
12(k),FDOA
23(k),FDOA
31(k)の計測値の中から、少なくとも1つの到来時間差および少なくとも1つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)を逐次算出する機能を有している。
【0107】
より具体的には、追尾演算部64Aは、
図7に示されるように単測位部93Aと反復推定部94Aとを有している。単測位部93Aの構成および機能は、単測位部91Aのそれらと同じである。反復推定部94Aは、単測位部93Aで生成された測位ベクトルPOS
kを用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)を算出することができる。
【0108】
反復推定部94Aは、上式(9),(10)で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部94Aに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(43),(44),(45)で与えられる。
【0109】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度、VLN
kは経度方向の速度、LT
kは緯度方向の速度であり、T=t
k−t
k−1である。また、I
2は2行2列の単位行列であり、qはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。目標Tgtは移動状態にあるので、状態遷移行列F
kは、等速直進運動を想定して設定されている。
【0110】
また、反復推定部94Aに適用される観測モデル(観測方程式)では、観測ベクトルz
k、観測行列H
kおよび共分散行列R
kは、たとえば、次式(46),(47),(48)で与えられる。
【0111】
ここで、OLN
k,OLT
kは、目標Tgtの観測位置を示す経度および緯度であり、式(48)におけるDOPは、精度劣化度(Dilution of Precision)と呼ばれ、本実施の形態では、単測位演算における理論的な測位誤差式によって与えられる。具体的には、DOPは、上式(36A)の誤差共分散行列Λに設定されればよい。
【0112】
反復推定部94Aは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部94Aは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)として出力する。
【0113】
次に、状態検出部62Aが目標Tgtの状態変化を検出した場合の測位演算部63A、追尾演算部64Aおよびデータ移行部65の動作について説明する。
【0114】
目標Tgtの静止状態から移動状態への状態変化が検出されたとき、測位演算部63Aの反復推定部92Aは、検出信号ESに応じて、直近の測位演算結果である状態推定ベクトルx
po(k)(=x(k|k))および事後誤差共分散行列P
po(k)(=P(k|k))をデータ移行部65を介して追尾演算部64Aに供給する。
【0115】
このとき、追尾演算部64Aの反復推定部94Aは、測位演算部63Aから供給された状態推定ベクトルx
po(k)および事後誤差共分散行列P
po(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部94Aは、それら初期値x(0|0),P(0|0)と測位ベクトルPOS
kとを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
trk(k)を算出することができる。
【0116】
この場合に、反復推定部94Aで使用される初期値x(0|0),P(0|0)は、たとえば、次式(49),(50)に示すように設定可能である。
【0117】
ここで、0
2×1は2行1列の零行列、0
2×2は2行2列の零行列である。式(49)の状態推定ベクトルx(0|0)は、目標Tgtの初期速度が零となるように設定されている。また、式(50)の初期値P(0|0)は、推定速度成分に誤差はないとの仮定に基づいて設定されている。
【0118】
一方、目標Tgtの移動状態から静止状態への状態変化が検出されたとき、追尾演算部64Aの反復推定部94Aは、検出信号ESに応じて、直近の追尾演算結果である状態推定ベクトルx
trk(k)(=x(k|k))および事後誤差共分散行列P
trk(k)(=P(k|k))をデータ移行部65を介して測位演算部63Aに供給する。
【0119】
このとき、測位演算部63Aにおける反復推定部92Aは、追尾演算部64Aから供給された状態推定ベクトルx
trk(k)および事後誤差共分散行列P
trk(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部92Aは、それら初期値x(0|0),P(0|0)と測位ベクトルPOS
kとを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
po(k)を算出することができる。
【0120】
この場合に、反復推定部92Aで使用される初期値x(0|0),P(0|0)は、たとえば、次式(51),(52)に示すように設定可能である。
【0121】
ここで、[x
trk(k)]
1,1は、状態推定ベクトルx
trk(k)の1行1列目成分(1番目成分)、[x
trk(k)]
2,1は、状態推定ベクトルx
trk(k)の2行1列目成分(2番目成分)である。また、[P
trk(k)]
1,1は、事後誤差共分散行列P
trk(k)の1行1列目成分、[P
trk(k)]
2,1は、事後誤差共分散行列P
trk(k)の2行1列目成分、[P
trk(k)]
1,2は、事後誤差共分散行列P
trk(k)の1行2列目成分、[P
trk(k)]
2,2は、事後誤差共分散行列P
trk(k)の2行2列目成分である。
【0122】
以上に説明したように、目標Tgtの静止状態から移動状態への状態変化が検出されたとき、追尾演算部64Aの反復推定部94Aは、直近の測位演算結果x
po(k),P
po(k)を初期値として用いかつ計測値を用いた追尾演算を実行することができる。逆に、目標Tgtの移動状態から静止状態への状態変化が検出されたときは、測位演算部63Aの反復推定部92Aは、直近の追尾演算結果x
trk(k),P
trk(k)を初期値として用いかつ計測値を用いた測位演算を実行することができる。
【0123】
上記した監視処理部11Aのハードウェア構成は、たとえば、DSP(Digital Signal Processor),ASIC(Application Specific Integrated Circuit)またはFPGA(Field−Programmable Gate Array)などの半導体集積回路を有する単数または複数のプロセッサで実現されればよい。あるいは、監視処理部11Aのハードウェア構成は、不揮発性メモリから読み出されたソフトウェアまたはファームウェアのプログラムコードを実行する、CPU(Central Processing Unit)またはGPU(Graphics Processing Unit)などの演算装置を含む単数または複数のプロセッサで実現されてもよい。DSPなどの半導体集積回路とCPUなどの演算装置との組み合わせを含む単数または複数のプロセッサで監視処理部11Aのハードウェア構成が実現されてもよい。
【0124】
図8は、監視処理部11Aの機能を実現するハードウェア構成例である信号処理装置70の概略構成を示すブロック図である。
図8に示される信号処理装置70は、プロセッサ71、メモリ72、入力インタフェース部73、出力インタフェース部74、通信回路75および信号路76を含んで構成されている。信号路76は、プロセッサ71、メモリ72、入力インタフェース部73、出力インタフェース部74および通信回路75を相互に接続するためのバスである。入力インタフェース部73は、外部から入力されたディジタル受信信号D1,D2,D3を信号路76を介してプロセッサ71に転送する機能を有する。プロセッサ71は、転送されたディジタル受信信号D1,D2,D3に基づいて測位演算および追尾演算を実行することで目標Tgtの静止状態および移動状態を示す状態推定値を算出することができる。また、プロセッサ71は、算出された状態推定値を基に画像情報DDを生成し、この画像情報DDを信号路75および出力インタフェース部74を介して外部に出力することができる。
【0125】
ここで、通信回路75は、外部サーバ(図示せず)と通信して状態検出部62Aでの処理に必要なデータを受信する回路である。また、メモリ72は、プロセッサ71がディジタル信号処理を実行する際に使用されるデータ記憶領域である。プロセッサ71がCPUなどの演算装置を内蔵する場合には、メモリ72は、プロセッサ71により実行されるソフトウェアまたはファームウェアのプログラムコードを記憶するデータ記憶領域を有する。メモリ72としては、たとえば、ROM(Read Only Memory)およびSDRAM(Synchronous Dynamic Random Access Memory)などの半導体メモリを使用することが可能である。
【0126】
以上に説明したように実施の形態1の目標監視装置2においては、状態検出部62Aは、複数の計測値のうちの少なくとも1つの計測値の時間変動を監視して目標Tgtの動きの有無を検出し、その検出結果を示す検出信号ESを測位演算部63A、追尾演算部64Aおよびデータ移行部65に供給するので、測位演算部63Aおよび追尾演算部64Aは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、監視処理部11Aは、目標Tgtの静止状態から移動状態への状態変化、あるいは目標Tgtの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。
【0127】
また、目標Tgtの静止状態から移動状態への状態変化が検出されたとき、追尾演算部64Aの反復推定部94Aは、直近の測位演算結果x
po(k),P
po(k)を初期値として用いかつ計測値の組を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Tgtに対する高い追尾精度を確保することができる。
【0128】
実施の形態2.
次に、本発明に係る実施の形態2について説明する。
図9は、本発明に係る実施の形態2における測位演算部63B、追尾演算部64Bおよびデータ移行部65からなる演算部の構成を概略的に示すブロック図である。本実施の形態の目標監視システムおよび目標監視装置(監視処理部を含む。)の構成は、
図3の測位演算部63Aおよび追尾演算部64Aに代えて、
図9の測位演算部63Bおよび追尾演算部64Bを有する点を除いて、上記実施の形態1の目標監視システム1および目標監視装置2の構成と同じである。
【0129】
本実施の形態の測位演算部63Bは、
図9に示されるように反復推定部92Bを有する。反復推定部92Bは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)の計測値の中から2以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、反復推定部92Bは、互いに同期して得られる到来時間差TDOA
12(k),TDOA
23(k)の計測値を選択し、これら計測値の組を観測ベクトルとして用いた測位演算を実行することができる。
【0130】
あるいは、反復推定部92Bは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)および到来周波数差FDOA
12(k),FDOA
23(k),FDOA
31(k)の計測値の中から、少なくとも1つの到来時間差および少なくとも1つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を逐次算出する機能を有する。これら計測値は、互いに同期して、または互いにほぼ同期して得られる値である。たとえば、反復推定部92Bは、互いに同期して得られる到来時間差TDOA
12(k)および到来周波数差FDOA
12(k)の計測値を選択し、これら計測値の組を観測ベクトルとして用いた測位演算を実行することができる。
【0131】
より具体的には、測位演算部63Bの反復推定部92Bは、上式(9),(10)で示される状態空間モデルを用いてカルマンフィルタによる反復推定処理を測位演算として実行することができる。反復推定部92Bに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(53),(54),(55)で与えられる。
【0132】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度である。目標Tgtは静止状態にあるので、状態遷移行列F
kは単位行列である。また、目標Tgtの運動に不確かさがなく、雑音ベクトルw
kは零ベクトルとなる。よって、共分散行列Q
kは零行列である。
【0133】
また、反復推定部92Bに適用される観測モデル(観測方程式)では、観測ベクトルz
k、観測関数h[x
k]および共分散行列R
kは、たとえば、次式(56),(57),(58)で与えられる。
【0134】
ここで、観測関数h[x
k]は、上式(18)に従い、ヤコビアン行列H
kに変換される。
【0135】
反復推定部92Bは、上記状態空間モデルを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部92Bは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)として出力する。
【0136】
一方、追尾演算部64Bは、
図9に示されるように反復推定部94Bを有する。反復推定部94Bは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)の計測値の中から2以上の計測値の組を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Tgtの推定位置および移動状態(たとえば推定速度)を示す状態推定ベクトルx
trk(k)を逐次算出する機能を有する。当該選択された計測値は、互いに同期して、または互いにほぼ同期して得られる値である。あるいは、追尾演算部64Bは、計測部61から供給された到来時間差TDOA
12(k),TDOA
23(k),TDOA
31(k)および到来周波数差FDOA
12(k),FDOA
23(k),FDOA
31(k)の計測値の中から、少なくとも1つの到来時間差および少なくとも1つの到来周波数差の計測値を選択し、当該選択された計測値の組を観測ベクトルとして用いた追尾演算を実行することにより目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)を逐次算出する機能を有する。当該選択された計測値は、互いに同期して、または互いにほぼ同期して得られる値である。
【0137】
より具体的には、反復推定部94Bは、上式(9),(10)で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部94Bに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(59),(60),(61)で与えられる。
【0138】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度、VLN
kは経度方向の速度、LT
kは緯度方向の速度であり、T=t
k−t
k−1である。また、I
2は2行2列の単位行列であり、qはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。
【0139】
また、反復推定部92Bに適用される観測モデル(観測方程式)では、観測ベクトルz
k、観測関数h[x
k]および共分散行列R
kは、たとえば、次式(62),(63),(64)で与えられる。
【0140】
ここで、観測関数h[x
k]は、上式(18)に従い、ヤコビアン行列H
kに変換される。
【0141】
反復推定部94Bは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部94Bは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)として出力する。
【0142】
次に、目標Tgtの状態変化が検出された場合の測位演算部63B、追尾演算部64Bおよびデータ移行部65の動作について説明する。
【0143】
目標Tgtの静止状態から移動状態への状態変化が検出されたとき、測位演算部63Bの反復推定部92Bは、検出信号ESに応じて、直近の測位演算結果である状態推定ベクトルx
po(k)(=x(k|k))および事後誤差共分散行列P
po(k)(=P(k|k))をデータ移行部65を介して追尾演算部64Bに供給する。
【0144】
このとき、追尾演算部64Bの反復推定部94Bは、測位演算部63Bから供給された状態推定ベクトルx
po(k)および事後誤差共分散行列P
po(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部94Bは、それら初期値x(0|0),P(0|0)と観測ベクトルとを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
trk(k)を算出することができる。この場合に、反復推定部92Bで使用される初期値x(0|0),P(0|0)は、たとえば、上式(49),(50)に示すように設定可能である。
【0145】
一方、目標Tgtの移動状態から静止状態への状態変化が検出されたとき、追尾演算部64Bの反復推定部94Bは、検出信号ESに応じて、直近の追尾演算結果である状態推定ベクトルx
trk(k)(=x(k|k))および事後誤差共分散行列P
trk(k)(=P(k|k))をデータ移行部65を介して測位演算部63Bに供給する。
【0146】
このとき、測位演算部63Bの反復推定部92Bは、追尾演算部64Bから供給された状態推定ベクトルx
trk(k)および事後誤差共分散行列P
trk(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部92Bは、それら初期値x(0|0),P(0|0)と観測ベクトルとを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
po(k)を算出することができる。この場合に、反復推定部92Bで使用される初期値x(0|0),P(0|0)は、たとえば、次式(51),(52)に示すように設定可能である。
【0147】
以上に説明したように実施の形態2の目標監視装置は、実施の形態1の目標監視装置2と同様に状態検出部62Aを有するので、測位演算部63Bおよび追尾演算部64Bは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態2の監視処理部は、目標Tgtの静止状態から移動状態への状態変化、あるいは目標Tgtの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Tgtの静止状態から移動状態への状態変化が検出されたとき、追尾演算部64Bの反復推定部94Bは、直近の測位演算結果x
po(k),P
po(k)を初期値として用い、かつ計測値の組を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Tgtに対する高い追尾精度を確保することができる。
【0148】
実施の形態3.
次に、本発明に係る実施の形態3について説明する。上記実施の形態1,2では、時刻t
k毎の観測ベクトルの成分として使用される複数の計測値(たとえば、到来時間差TDOA
12(k),TDOA
23(k)の計測値)は、互いに同期して、または互いにほぼ同期して得られる観測値である。これに対し、本実施の形態では、複数の計測値は、互いに非同期に得られる観測値、言い換えれば、同時刻に同時に得られない観測値である。
【0149】
たとえば、受信アンテナ(受信局)が3つある場合には、到来時間差TDOA
12(k),TDOA
23(k)の計測値を互いに同期する形で同時に得ることができる。これに対し、使用可能な受信アンテナ(受信局)が2つのみである場合には、衛星St1,St2間の計測値と衛星St2,St3間の計測値との組を得るために、アンテナ指向方向の切り替えを行う必要がある。しかしながら、その切り替えに時間ずれが生じると、それら2種類の計測値が互いに非同期に得られることがある。また、到来時間差TDOA
12(k),TDOA
23(k)の計測値が同時に得られたとしても、これら計測値の一方が外れ値であると判定された場合など何らかの理由で欠落する可能性がある。本実施の形態の監視処理部は、非同期に得られた複数の計測値を用いて測位演算および追尾演算を実行することができる。
【0150】
図10は、本発明に係る実施の形態3における監視処理部11Cの概略構成を示すブロック図である。本実施の形態の目標監視システムの構成は、
図3の監視処理部11Aに代えて
図10の監視処理部11Cを有する点を除いて、上記実施の形態1の目標監視システム1の構成と同じである。監視処理部11Cは、ディジタル受信信号D1,D2,D3に基づいて電波の到来時間差TDOA
ij(k)および到来周波数差FDOA
ij(k)の計測値を逐次算出する計測部61Cを備える。到来時間差TDOA
ij(k)は、到来時間差TDOA
12(k)またはTDOA
23(k)のいずれか一方であり、到来周波数差FDOA
ij(k)は、到来周波数差FDOA
12(k)またはFDOA
23(k)のいずれか一方である。たとえば、時刻t
k1に到来時間差TDOA
12(k1)の計測値が得られ、時刻t
k2(≠t
k1)に到来時間差TDOA
23(k2)の計測値が得られる場合、到来時間差TDOA
12(k1),TDOA
23(k2)の計測値は、互いに非同期に得られた観測値である。
【0151】
監視処理部11Cは、さらに、到来時間差TDOA
ij(k)および到来周波数差FDOA
ij(k)の計測値のうちのいずれかの計測値の時間変動を常時監視して目標Tgtの動きの有無を検出する状態検出部62Aと、目標Tgtの動きが無いことが検出されたときは、当該計測値を用いた測位演算を実行して目標Tgtの推定位置を示す状態推定値を逐次算出する測位演算部63Cと、目標Tgtの動きが有ることが検出されたきは、当該計測値を用いた追尾演算を実行して目標Tgtの推定位置および移動状態(たとえば推定速度)の双方を示す状態推定値を算出する追尾演算部64Cと、測位演算部63Cおよび追尾演算部64Cの一方から他方へ演算結果のデータを移行させるデータ移行部65と、状態検出部62Aでの処理に必要なデータを格納するデータ記憶部67とを有する。
【0152】
図11は、実施の形態3における測位演算部63C、追尾演算部64Cおよびデータ移行部65からなる演算部の構成を概略的に示すブロック図である。本実施の形態の測位演算部63Cは反復推定部92Cを有する。反復推定部92Cは、計測部61で非同期に得られた到来時間差TDOA
12(k1),TDOA
23(k2)の計測値(k1≠k2)の組を用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を算出する機能を有する。
【0153】
あるいは、反復推定部92Cは、計測部61で非同期に得られた到来時間差TDOA
12(k1)と到来周波数差FDOA
23(k2)の計測値(k1≠k2)の組を用いたカルマンフィルタによる反復推定処理を測位演算として実行することにより、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を算出してもよい。
【0154】
より具体的には、測位演算部63Cの反復推定部92Cは、上式(9),(10)で示される状態空間モデルを用いてカルマンフィルタによる反復推定処理を測位演算として実行することができる。反復推定部92Cに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(65),(66),(67)で与えられる。
【0155】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度である。目標Tgtは静止状態にあるので、状態遷移行列F
kは単位行列である。また、目標Tgtの運動に不確かさがなく、雑音ベクトルw
kは零ベクトルとなる。よって、共分散行列Q
kは零行列である。
【0156】
また、反復推定部92Cに適用される観測モデル(観測方程式)では、観測値z
k、観測関数h[x
k]および共分散R
kは、たとえば、次式(68),(69),(70)で与えられる。
【0157】
ここで、式(68),(69)における(i,j)は、各時刻t
kで(1,2)または(2,3)のいずれかである。たとえば、時刻t
k1では、到来時間差TDOA
12(k1)の観測値(計測値)が得られ(i=1;j=2)、時刻t
k2ではTDOA
23(k2)の観測値(計測値)が得られる(i=2;j=3)。
【0158】
反復推定部92Cは、上記状態空間モデルを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部92Cは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)として出力する。
【0159】
一方、追尾演算部64Cは、
図11に示されるように反復推定部94Cを有する。反復推定部94Cは、計測部61で非同期に得られた到来時間差TDOA
12(k1),TDOA
23(k2)の計測値(k1≠k2)の組を用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Tgtの推定位置および移動状態(たとえば推定速度)を示す状態推定ベクトルx
trk(k)を逐次算出する機能を有する。
【0160】
あるいは、追尾演算部64Cは、計測部61で非同期に得られた到来時間差TDOA
12(k1)および到来周波数差FDOA
23(k2)の計測値(k1≠k2)の組を用いたカルマンフィルタによる反復推定処理を追尾演算として実行することにより、目標Tgtの推定位置および移動状態(たとえば推定速度)を示す状態推定ベクトルx
trk(k)を算出してもよい。
【0161】
より具体的には、反復推定部94Cは、上式(9),(10)で示される状態空間モデルを用いたカルマンフィルタによる反復推定処理を追尾演算として実行する。反復推定部94Cに適用される運動モデル(状態方程式)では、状態ベクトルx
k、状態遷移行列F
kおよび共分散行列Q
kは、たとえば、次式(71),(72),(73)で与えられる。
【0162】
ここで、LN
k,LT
kは、目標Tgtの真の位置を示す経度および緯度、VLN
kは経度方向の速度、LT
kは緯度方向の速度であり、T=t
k−t
k−1である。また、I
2は2行2列の単位行列であり、qはパワースペクトラム密度と呼ばれ、運動モデルの不確かさを表す設定パラメータ値である。
【0163】
また、反復推定部94Cに適用される観測モデル(観測方程式)では、観測値z
k、観測関数h[x
k]および共分散R
kは、たとえば、次式(74),(75),(76)で与えられる。
【0164】
ここで、式(74),(75)における(i,j)は、各時刻t
kで(1,2)または(2,3)のいずれかである。たとえば、時刻t
k1では、到来時間差TDOA
12(k1)の計測値が得られ(i=1;j=2)、時刻t
k2ではTDOA
23(k2)の計測値が得られる(i=2;j=3)。
【0165】
反復推定部94Cは、上記状態空間モデルを用いてカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)を算出することができる。反復推定部94Cは、当該状態推定ベクトルx(k|k)を、目標Tgtの推定位置および移動状態を示す状態推定ベクトルx
trk(k)として出力する。
【0166】
次に、目標Tgtの状態変化が検出された場合の測位演算部63C、追尾演算部64Cおよびデータ移行部65の動作について説明する。
【0167】
目標Tgtの静止状態から移動状態への状態変化が検出されたとき、測位演算部63Cの反復推定部92Cは、検出信号ESに応じて、直近の測位演算結果である状態推定ベクトルx
po(k)(=x(k|k))および事後誤差共分散行列P
po(k)(=P(k|k))をデータ移行部65を介して追尾演算部64Cに供給する。
【0168】
このとき、追尾演算部64Cの反復推定部94Cは、測位演算部63Cから供給された状態推定ベクトルx
po(k)および事後誤差共分散行列P
po(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部94Cは、それら初期値x(0|0),P(0|0)と観測値(計測値)とを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
trk(k)を算出することができる。この場合に、反復推定部92Bで使用される初期値x(0|0),P(0|0)は、たとえば、上式(49),(50)に示すように設定可能である。
【0169】
一方、目標Tgtの移動状態から静止状態への状態変化が検出されたとき、追尾演算部64Cの反復推定部94Cは、検出信号ESに応じて、直近の追尾演算結果である状態推定ベクトルx
trk(k)(=x(k|k))および事後誤差共分散行列P
trk(k)(=P(k|k))をデータ移行部65を介して測位演算部63Cに供給する。
【0170】
このとき、測位演算部63Cの反復推定部92Cは、追尾演算部64Cから供給された状態推定ベクトルx
trk(k)および事後誤差共分散行列P
trk(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部92Cは、それら初期値x(0|0),P(0|0)と観測値(計測値)とを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
po(k)を算出することができる。この場合に、反復推定部92Bで使用される初期値x(0|0),P(0|0)は、たとえば、次式(51),(52)に示すように設定可能である。
【0171】
以上に説明したように実施の形態3の目標監視装置は、実施の形態1の目標監視装置2と同様に状態検出部62Aを有するので、測位演算部63Cおよび追尾演算部64Cは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態3の監視処理部11Cは、目標Tgtの静止状態から移動状態への状態変化、あるいは目標Tgtの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Tgtの静止状態から移動状態への状態変化が検出されたとき、追尾演算部64Cの反復推定部94Cは、直近の測位演算結果x
po(k),P
po(k)を初期値として用い、かつ計測値を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Tgtに対する高い追尾精度を確保することができる。
【0172】
実施の形態4.
次に、本発明に係る実施の形態4について説明する。上記実施の形態3では、測位演算部63Cは、各時刻t
kに計測値が入力される度に測位演算を逐次的に実行する。これに対し、実施の形態4の測位演算は、複数回分の計測値を一括して処理するバッチ型の測位演算である。
【0173】
図12は、本発明に係る実施の形態4における測位演算部63D、追尾演算部64Cおよびデータ移行部65からなる演算部の構成を概略的に示すブロック図である。本実施の形態の目標監視システムおよび目標監視装置の構成は、
図10の測位演算部63Cに代えて
図10の測位演算部63Dを有する点を除いて、上記実施の形態3の目標監視システムおよび目標監視装置の構成と同じである。
【0174】
本実施の形態の測位演算部63Dには、計測部61C(
図10)から到来時間差TDOA
ij(k)および到来周波数差FDOA
ij(k)の計測値が入力される。到来時間差TDOA
ij(k)は、到来時間差TDOA
12(k)またはTDOA
23(k)のいずれか一方であり、到来周波数差FDOA
ij(k)は、到来周波数差FDOA
12(k)またはFDOA
23(k)のいずれか一方である。たとえば、時刻t
k1に到来時間差TDOA
12(k1)の計測値が得られ、時刻t
k2(≠t
k1)に到来時間差TDOA
23(k2)の計測値が得られる場合、到来時間差TDOA
12(k1),TDOA
23(k2)の計測値は、互いに非同期に得られた観測値である。
【0175】
本実施の形態の測位演算部63Dは、
図12に示されるようにバッファメモリ100とバッチ型反復推定部92Dとを有している。バッファメモリ100は、少なくとも、時刻t
k−K+1から時刻t
k−1までのK−1回分の過去の計測値を記憶する容量を有している。バッチ型反復推定部92Dは、時刻t
kにおける現在の計測値が得られる度に、バッファメモリ100からK−1回分の過去の計測値を読み出す。そして、バッチ型反復推定部92Dは、当該過去の計測値と当該現在の計測値とを一括して用いた非線形最小二乗法によるバッチ処理を測位演算として実行することにより、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)を算出することができる。
【0176】
より具体的には、バッチ型反復推定部92Dは、上記した再帰最小二乗法(RLS)アルゴリズムに従い、次式(77)に示される観測値y
kを用いた再帰処理を実行することができる。
【0178】
このとき、時刻t
kにおける観測関数φ
k[χ=p
tgt]は、次式(78)で与えられる。
【0179】
ここで、p
tgtは、地球の重心を原点とする、目標Tgtの位置座標を示す位置ベクトルであり、p
i(k)は、時刻t
kでの衛星Stiの位置座標を示す位置ベクトルであり、p
j(k)は、時刻t
kでの衛星Stjの位置座標を示す位置ベクトルである。
【0180】
観測雑音共分散Ω
kは、次式(79)で与えられるものとする。
【0182】
バッチ型反復推定部92Dは、K回の計測結果に対して上記したRLSアルゴリズムによる再帰処理を実行することにより、評価関数J(χ=p
tgt)を最小にする近似解χ
minを、目標Tgtの推定位置を示す状態推定ベクトルx
po(k)として算出することができる。
【0183】
次に、目標Tgtの状態変化が検出された場合の測位演算部63D、追尾演算部64Cおよびデータ移行部65の動作について説明する。
【0184】
目標Tgtの静止状態から移動状態への状態変化が検出されたとき、測位演算部63Dの反復推定部92Dは、検出信号ESに応じて、直近の測位演算結果である状態推定ベクトルx
po(k)(=x(k|k))および事後誤差共分散行列P
po(k)(=P(k|k))をデータ移行部65を介して追尾演算部64Cに供給する。
【0185】
このとき、追尾演算部64Cの反復推定部94Cは、測位演算部63Dから供給された状態推定ベクトルx
po(k)および事後誤差共分散行列P
po(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部94Cは、それら初期値x(0|0),P(0|0)と観測値(計測値)とを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
trk(k)を算出することができる。
【0186】
一方、目標Tgtの移動状態から静止状態への状態変化が検出されたとき、追尾演算部64Cの反復推定部94Cは、検出信号ESに応じて、直近の追尾演算結果である状態推定ベクトルx
trk(k)(=x(k|k))および事後誤差共分散行列P
trk(k)(=P(k|k))をデータ移行部65を介して測位演算部63Dに供給する。
【0187】
このとき、測位演算部63Dの反復推定部92Dは、追尾演算部64Cから供給された状態推定ベクトルx
trk(k)および事後誤差共分散行列P
trk(k)に基づいて、状態推定ベクトルx(k|k)および事後誤差共分散行列P(k|k)の初期値x(0|0),P(0|0)を生成する。そして、反復推定部92Dは、それら初期値x(0|0),P(0|0)と観測値(計測値)とを用いたカルマンフィルタによる反復推定処理(
図6)を実行することにより、時刻t
k毎の状態推定ベクトルx
po(k)を算出することができる。この場合に、反復推定部92Dで使用される初期値x(0|0),P(0|0)は、たとえば、次式(51),(52)に示すように設定可能である。
【0188】
以上に説明したように実施の形態4の目標監視装置は、実施の形態1の目標監視装置2と同様に状態検出部62Aを有するので、測位演算部63Dおよび追尾演算部64Cは、測位演算を追尾演算に切り替えるべきタイミング、あるいは、追尾演算を測位演算に切り替えるべきタイミングを高い確度で決定することができる。したがって、実施の形態4の監視処理部は、目標Tgtの静止状態から移動状態への状態変化、あるいは目標Tgtの移動状態から静止状態への状態変化に応じて精度の高い演算を実行することができる。また、目標Tgtの静止状態から移動状態への状態変化が検出されたとき、追尾演算部64Cの反復推定部94Cは、直近の測位演算結果x
po(k),P
po(k)を初期値として用い、かつ計測値を用いた追尾演算を実行することができる。したがって、測位演算が追尾演算に切り替えられたとしても、目標Tgtに対する高い追尾精度を確保することができる。
【0189】
実施の形態5.
次に、本発明に係る実施の形態5について説明する。実施の形態1の状態検出部62A(
図4)においては、閾値設定部82Aは、監視エリアSAを定める監視エリア情報ADを用いて閾値TH1,TH2を設定し、閾値設定部86Aも、監視エリア情報ADを用いて閾値TH3,TH4を設定している。実施の形態5,6では、監視エリア情報ADに代えて、直近の測位演算結果または直近の追尾演算結果を用いて閾値TH1,TH2,TH3,TH4が設定される。
【0190】
図13は、実施の形態5における監視処理部11Eの概略構成を示すブロック図である。監視処理部11Eの構成は、
図3の状態検出部62Aに代えて
図13の状態検出部62Eおよび遅延素子78,79を有する点を除いて、実施の形態1の監視処理部11Aの構成と同じである。
図13に示されるように、一方の遅延素子78は、現在の時刻t
kよりも1回前の時刻t
k−1での測位演算結果である状態推定ベクトルx
po(k−1)を状態検出部62Eに供給し、他方の遅延素子79は、時刻t
k−1での追尾演算結果である状態推定ベクトルx
trk(k−1)を状態検出部62Eに供給する。
【0191】
図14は、実施の形態5における状態検出部62Eの概略構成を示すブロック図である。状態検出部62Eの構成は、
図4の閾値設定部82A,86Aに代えて
図14の閾値設定部82E,86Eを有する点を除いて、上記状態検出部62Aの構成と同じである。閾値設定部82Eは、測位演算結果x
po(k−1)または追尾演算結果x
trk(k−1)のいずれか一方と衛星軌道情報ODとに基づき、上式(1),(2)に示される閾値TH1,TH2の組を設定することができる。また、閾値設定部86Eは、測位演算結果x
po(k−1)または追尾演算結果x
trk(k−1)のいずれか一方と衛星軌道情報ODと目標揺動情報SDとに基づき、上式(5),(6)に示される閾値TH3,TH4の組を設定することができる。
【0192】
以上に説明したように実施の形態5では、目標Tgtに関する状態推定ベクトルx
po(k−1),x
trk(k−1)を使用して閾値TH1,TH2,TH3,TH4が設定されるので、目標Tgtの状態変化を検知するための信頼性の高い数値範囲Δ(TH1,TH2),Δ(TH3,TH4)を決定することが可能となる。
【0193】
実施の形態6.
次に、本発明に係る実施の形態6について説明する。
図15は、実施の形態6における監視処理部11Fの概略構成を示すブロック図である。監視処理部11Fの構成は、
図3の状態検出部62Aに代えて
図15の状態検出部62Fおよび遅延素子78,79,91,92を有する点を除いて、実施の形態1の監視処理部11Aの構成と同じである。
図15に示されるように、遅延素子78,91は、現在の時刻t
kよりも1回前の時刻t
k−1での測位演算結果である状態推定ベクトルx
po(k−1)および事後誤差共分散行列P
po(k−1)を状態検出部62Eに供給し、遅延素子79,92は、時刻t
k−1での追尾演算結果である状態推定ベクトルx
trk(k−1)および事後誤差共分散行列P
trk(k−1)を状態検出部62Eに供給する。
【0194】
図16は、実施の形態6における状態検出部62Fの概略構成を示すブロック図である。状態検出部62Fの構成は、
図4の閾値設定部82A,86Aに代えて
図16の閾値設定部82E,86Eおよび誤差楕円算出部84,88を有する点を除いて、上記状態検出部62Aの構成と同じである。
【0195】
誤差楕円算出部84は、測位演算結果x
po(k−1),P
po(k−1)に基づいて、次式(80)に示される誤差楕円で定まる領域データを閾値設定部82F,86Fに供給する。
【0196】
ここで、ζは、設定パラメータである。誤差楕円とは、測位結果の誤差範囲を示す確率的な範囲のことであり、測位位置を中心とする楕円状の領域内に目標Tgtの真の位置があることを意味する。
【0197】
一方、誤差楕円算出部88は、追尾演算結果x
trk(k−1),P
trk(k−1)に基づいて、次式(81)に示される誤差楕円で定まる領域データを閾値設定部82F,86Fに供給する。
【0198】
ここで、ζ’は、設定パラメータである。
【0199】
閾値設定部82Fは、式(80)または(81)で示される誤差楕円で定まる領域データと衛星軌道情報ODとに基づき、上式(1),(2)に示される閾値TH1,TH2の組を設定することができる。また、閾値設定部86Fは、式(80)または(81)で示される誤差楕円で定まる領域データと衛星軌道情報ODと目標揺動情報SDとに基づき、上式(5),(6)に示される閾値TH3,TH4の組を設定することができる。
【0200】
以上に説明したように実施の形態6では、誤差楕円で定まる領域データを使用して閾値TH1,TH2,TH3,TH4が設定されるので、目標Tgtの状態変化を検知するための信頼性の高い数値範囲Δ(TH1,TH2),Δ(TH3,TH4)を決定することが可能となる。また、領域データを人為的に設定する必要がないという利点がある。
【0201】
実施の形態1〜6の変形例.
以上、図面を参照して本発明に係る種々の実施の形態について述べたが、上記実施の形態は本発明の例示であり、これら実施の形態以外の様々な形態を採用することもできる。
【0202】
たとえば、測位演算部63A〜63Dのうちから選択された任意の測位演算部と、追尾演算部64A〜63Cのうちから選択された任意の追尾演算部と、データ移行部65との組み合わせを有する監視処理部を構成することが可能である。
【0203】
また、上記実施の形態2〜5の各監視処理部のハードウェア構成は、実施の形態1の場合と同様に、たとえば、DSP,ASICまたはFPGAなどの半導体集積回路を有する単数または複数のプロセッサで実現されればよい。あるいは、監視処理部11Aのハードウェア構成は、不揮発性メモリから読み出されたソフトウェアまたはファームウェアのプログラムコードを実行する、CPUまたはGPUなどの演算装置を含む単数または複数のプロセッサで実現されてもよい。DSPなどの半導体集積回路とCPUなどの演算装置との組み合わせを含む単数または複数のプロセッサで各監視処理部のハードウェア構成が実現されてもよい。さらには、実施の形態2〜5の各監視処理部のハードウェア構成は、
図8に示した信号処理装置70により実現されてもよい。
【0204】
本発明の範囲内において、上記実施の形態1〜6の自由な組み合わせ、各実施の形態の任意の構成要素の変形、または各実施の形態の任意の構成要素の省略が可能である。
目標監視装置は、目標から複数の通信経路を経て到来した電波信号をそれぞれ受信する複数の受信アンテナで得られた複数の受信信号に基づいて当該目標の状態を推定する。目標監視装置は、複数の受信信号に基づいて複数の計測値を逐次算出する計測部(61)、複数の計測値に含まれる到来時間差または到来周波数差を監視して目標の動きの有無を検出する状態検出部(62A)と、測位演算および追尾演算を実行する演算部(63A,64A,65)とを備える。演算部(63A,64A,65)は、目標の動きが無いことが検出されたときは、測位演算を実行して目標の推定位置を示す状態推定値を算出し、目標の静止状態から移動状態への変化が検出されたときは、追尾演算を実行して当該目標の推定位置および移動状態の双方を示す状態推定値を逐次算出する。