(57)【特許請求の範囲】
【請求項１】
観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データを入力し、前記観測値の観測時刻の時間間隔を算出する第１の計算部と、
観測誤差標準偏差、駆動雑音標準偏差および観測時刻の時間間隔の関数として表されるトラッキングインデックスを用いて、複数の前記観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように、各観測時刻における前記駆動雑音標準偏差を算出する第２の計算部と、
前記第１の計算部によって算出された観測時刻の時間間隔と、前記第２の計算部によって算出された各観測時刻における前記駆動雑音標準偏差と、前記観測対象の状態の平滑ベクトルおよび平滑誤差共分散行列とを用いて、前記観測対象の状態の予測ベクトルおよび予測誤差共分散行列を算出する予測部と、
各観測時刻における前記観測誤差標準偏差と、前記予測部によって算出された前記予測ベクトルおよび前記予測誤差共分散行列とを用いて、前記平滑ベクトルおよび前記平滑誤差共分散行列を算出する平滑部と、
を備えたことを特徴とする平滑フィルタ装置。
【請求項２】
複数の前記観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値になるように設定された前記トラッキングインデックスの初期値を算出し、前記トラッキングインデックスの初期値を前記第２の計算部に設定する初期値設定部を備え、
前記第２の計算部は、前記トラッキングインデックスの初期値、前記観測誤差標準偏差および前記第１の計算部によって算出された観測時刻の時間間隔を用いて、前記駆動雑音標準偏差を算出すること
を特徴とする請求項１記載の平滑フィルタ装置。
【請求項３】
前記初期値設定部は、雑音抑圧比の平方根と前記観測誤差標準偏差の積である前記平滑誤差標準偏差を算出し、複数の前記観測方法にそれぞれ対応する前記平滑誤差標準偏差が同じ値となるように、前記トラッキングインデックスの初期値を前記観測方法ごとに算出すること
を特徴とする請求項２記載の平滑フィルタ装置。
【請求項４】
第１の計算部が、観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データを入力し、前記観測値の観測時刻の時間間隔を算出するステップと、
第２の計算部が、観測誤差標準偏差、駆動雑音標準偏差および観測時刻の時間間隔の関数として表されるトラッキングインデックスを用いて、複数の前記観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように、各観測時刻における前記駆動雑音標準偏差を算出するステップと、
予測部が、前記第１の計算部によって算出された観測時刻の時間間隔と、前記第２の計算部によって算出された各観測時刻における前記駆動雑音標準偏差と、前記観測対象の状態の平滑ベクトルおよび平滑誤差共分散行列とを用いて、前記観測対象の状態の予測ベクトルおよび予測誤差共分散行列を算出するステップと、
平滑部が、各観測時刻における前記観測誤差標準偏差と、前記予測部によって算出された前記予測ベクトルおよび前記予測誤差共分散行列とを用いて、前記平滑ベクトルおよび前記平滑誤差共分散行列を算出するステップと、
を備えたことを特徴とする平滑化方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、観測値の時系列データの平滑処理を行う平滑フィルタ装置および平滑化方法に関する。
【背景技術】
【０００２】
従来から、センサによって観測された観測値の時系列データのばらつきを抑えることを目的として、平滑フィルタによって、観測値の時系列データの平滑処理が行われている。例えば、非特許文献１には、平滑フィルタの一つであるカルマンフィルタが記載されている。平滑処理には、平滑フィルタに設定された観測誤差共分散行列および駆動雑音共分散行列の情報が用いられる。【０００３】
観測誤差共分散行列は、観測値の真値からの誤差の共分散行列であって、観測誤差標準偏差を用いて導出される。観測誤差標準偏差は、観測値を得る観測方法の観測精度を示すパラメータである。例えば、観測精度が異なる観測方法によって得られた観測値が混在した時系列データを平滑処理する場合、平滑フィルタには、観測方法ごとの観測誤差共分散行列が設定される。【０００４】
駆動雑音共分散行列は、観測対象の運動の曖昧さを表す駆動雑音ベクトルの共分散行列であり、駆動雑音標準偏差を用いて導出される。駆動雑音共分散行列は、観測対象の運動モデルの曖昧さの大きさを表すフィルタ設定パラメータである。平滑フィルタの平滑性能は、観測誤差共分散行列を導出するための観測誤差標準偏差の設定値と、駆動雑音共分散行列を導出するための駆動雑音標準偏差の設定値に依存する。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Ｆｒａｎｋ  Ｌ．  Ｌｅｗｉｓ，  “Ｏｐｔｉｍａｌ  Ｅｓｔｉｍａｔｉｏｎ：  Ｗｉｔｈ  ａｎ  Ｉｎｔｒｏｄｕｃｔｉｏｎ  ｔｏ  Ｓｔｏｃｈａｓｔｉｃ  Ｃｏｎｔｒｏｌ  Ｔｈｅｏｒｙ，”  Ｗｉｌｅｙ−Ｉｎｔｅｒｓｃｉｅｎｃｅ，  １９８６．
【発明の概要】
【発明が解決しようとする課題】
【０００６】
従来の平滑フィルタ装置は、観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データを平滑処理する場合、観測方法の違いに応じた駆動雑音共分散行列の調整を行わない。このため、平滑フィルタの平滑性能が観測方法ごとに異なって均一な平滑化が行えない可能性があるという課題があった。【０００７】
本発明は上記課題を解決するものであり、観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データの平滑処理において、複数の観測方法にそれぞれ対応する平滑性能を一致させることができる平滑フィルタ装置および平滑化方法を得ることを目的とする。
【課題を解決するための手段】
【０００８】
本発明に係る平滑フィルタ装置は、観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データを入力し、観測値の観測時刻の時間間隔を算出する第１の計算部と、観測誤差標準偏差、駆動雑音標準偏差および観測時刻の時間間隔の関数として表されるトラッキングインデックスを用いて、複数の観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように、各観測時刻における駆動雑音標準偏差を算出する第２の計算部と、第１の計算部によって算出された観測時刻の時間間隔と、第２の計算部によって算出された各観測時刻における駆動雑音標準偏差と、観測対象の状態の平滑ベクトルおよび平滑誤差共分散行列とを用いて、観測対象の状態の予測ベクトルおよび予測誤差共分散行列を算出する予測部と、各観測時刻における観測誤差標準偏差と、予測部によって算出された予測ベクトルおよび予測誤差共分散行列とを用いて、平滑ベクトルおよび平滑誤差共分散行列を算出する平滑部を備える。
【発明の効果】
【０００９】
本発明によれば、トラッキングインデックスを用いて、複数の観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように、各観測時刻における駆動雑音標準偏差が算出され、算出された駆動雑音標準偏差を用いて観測値の時系列データの平滑処理が行われる。観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データの平滑処理において、複数の観測方法にそれぞれ対応する平滑性能を一致させることができる。
【図面の簡単な説明】
【００１０】
【図１】実施の形態１に係る平滑フィルタ装置の構成を示すブロック図である。
【図２】実施の形態１に係る平滑フィルタ装置による平滑処理を示すイメージ図である。
【図３】ランダム雑音に対する平滑処理を示すイメージ図である。
【図４】トラッキングインデックスの初期値の設定処理を示すフローチャートである。
【図５】実施の形態１に係る平滑化方法を示すフローチャートである。
【図６】図６Ａは、実施の形態１に係る平滑フィルタ装置の機能を実現するハードウェア構成を示すブロック図であり、図６Ｂは、実施の形態１に係る平滑フィルタ装置の機能を実現するソフトウェアを実行するハードウェア構成を示すブロック図である。
【図７】航空監視システムにおける航空機の測位の概要を示すイメージ図である。
【図８】航空監視システムにおける受信局の時刻管理の概要を示すイメージ図である。
【図９】クロックオフセットの観測時刻と観測値の平滑誤差絶対値との関係を示すグラフである。
【図１０】クロックオフセットの観測値の平滑処理のモンテカルロ試行した結果を示すグラフである。
【発明を実施するための形態】
【００１１】
実施の形態１．
図１は、実施の形態１に係る平滑フィルタ装置１の構成を示すブロック図である。平滑フィルタ装置１は、例えば、運動モデルおよび観測モデルを用いて観測値の時系列データを平滑化するカルマンフィルタである。観測対象の運動モデルには、等速直線運動モデルがあり、等速直線運動モデルは、下記式（１）から（６）を用いて定義される。下記式において、ｔｋ（ｋ＝０，１，２，・・・）は、観測対象の観測時刻（サンプリング時刻）である。ｘｋは、時刻ｔｋにおける観測対象の状態（観測対象の位置および速度）を表す状態ベクトルである。Φｋは、時刻ｔｋから時刻ｔｋ＋１への状態ベクトルｘｋの遷移行列である。Ｔｋは、一つ前の時刻ｔｋ−１と時刻ｔｋとの差である。すなわち、Ｔｋは、観測時刻の時間間隔である。ｗｋは、駆動雑音ベクトルであり、観測対象の運動の曖昧さを示す。すなわち、駆動雑音ベクトルｗｋは、観測対象が等速直線運動から逸脱したときの誤差を表している。ｓｋは、観測対象の状態の真値であり、ｓｋドットは、真値ｓｋの時間変化率である。［Ａ］Ｔは、ベクトルＡの転置ベクトルである。Ｑｋは駆動雑音共分散行列であり、駆動雑音ベクトルｗｋの共分散行列である。駆動雑音共分散行列Ｑｋは、観測対象の運動モデルの曖昧さの大きさを表している。Ｅは平均を表す記号であり、ｑｋは、時刻ｔｋにおける駆動雑音のパワースペクトラム密度である。σｗ，ｋは、時刻ｔｋにおける駆動雑音標準偏差である。
0006952944.tif 000002
【００１２】
観測対象を観測するセンサの観測モデルは、例えば、下記式（７）から（９）で表すことができる。下記式において、ｚｋは、時刻ｔｋにおける観測値であり、時刻ｔｋにおける観測対象の状態の観測ベクトルである。Ｈは、状態ベクトルｘｋから位置ベクトルのみを抽出するための既知の観測行列である。ｖｋは、時刻ｔｋにおける観測値の観測雑音（観測誤差とも呼ばれる）である。Ｒｋは、時刻ｔｋにおける観測誤差共分散行列であり、スカラー値で表される。σｖ，ｋ２は、時刻ｔｋにおける観測誤差分散であり、σｖ，ｋは、時刻ｔｋにおける観測誤差標準偏差である。
0006952944.tif 000003
【００１３】
平滑フィルタ装置１は、図１に示すように、初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４を備える。予測部１３は、観測時刻の時間間隔Ｔｋおよび時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを用いて、上記式（２）および（４）から、時刻ｔｋにおける駆動雑音共分散行列Ｑｋおよび遷移行列Φｋを算出する。さらに、予測部１３は、時刻ｔｋにおける駆動雑音共分散行列Ｑｋおよび遷移行列Φｋと、時刻ｔｋにおける観測対象の平滑ベクトルｘｋハットおよび平滑誤差共分散行列Ｐｋハットとを用いて、下記式（１０）に従い、次時刻ｔｋ＋１における観測対象の予測ベクトルｘｋ＋１バーおよび予測誤差共分散行列Ｐｋ＋１バーを算出する。なお、ΦｋＴは、遷移行列Φｋの転置行列である。
0006952944.tif 000004
【００１４】
平滑部１４は、上記式（９）に次時刻ｔｋ＋１における観測誤差標準偏差σｖ，ｋ＋１を用いて、次時刻ｔｋ＋１における観測誤差共分散行列Ｒｋ＋１を算出する。平滑部１４は、次時刻ｔｋ＋１における観測誤差共分散行列Ｒｋ＋１と、次時刻ｔｋ＋１における観測行列と、次時刻ｔｋ＋１における観測対象の予測ベクトルｘｋ＋１バーおよび次時刻ｔｋ＋１における予測誤差共分散行列Ｐｋ＋１バーを用いて、下記式（１１）に従い、次時刻ｔｋ＋１における観測対象の平滑ベクトルｘｋ＋１ハットおよび平滑誤差共分散行列Ｐｋ＋１ハットを算出する。Ｈｋ＋１Ｔは、次時刻ｔｋ＋１における観測行列Ｈｋ＋１の転置行列である。
0006952944.tif 000005
【００１５】
平滑フィルタ装置１による平滑処理には、上記式（１０）および上記式（１１）に示すように駆動雑音共分散行列Ｑｋおよび観測誤差共分散行列Ｒｋが必要である。すなわち、平滑処理には、これらを導出するための駆動雑音標準偏差σｗ，ｋおよび観測誤差標準偏差σｖ，ｋの情報が必要である。平滑フィルタ装置１における平滑性能は、駆動雑音標準偏差σｗ，ｋおよび観測誤差標準偏差σｖ，ｋの設定値に依存する。
【００１６】
また、カルマンフィルタは、ゲインが予測誤差共分散行列Ｐｋ＋１バーと、観測誤差標準偏差σｖ，ｋとから決定されるα−βフィルタであると言える。このため、カルマンフィルタにおける定常ゲインＫと、カルマンフィルタに初期設定されるフィルタ設定パラメータとの間には、下記式（１２）および（１３）に示す関係がある。下記式において、αは、α−βフィルタにおける位置のゲインであり、βは、速度のゲインである。Ｔは、観測時刻の時間間隔（サンプリング間隔）である。
0006952944.tif 000006
【００１７】
トラッキングインデックスΓは、観測誤差と駆動雑音の挙動を表すパラメータであり、下記式（１４）に示すように、観測誤差標準偏差σｖ、駆動雑音標準偏差σｗおよび観測時刻の時間間隔Ｔの関数として表される。さらに、トラッキングインデックスΓは、αおよびβを用いた下記式（１４）に示す関係を有している。ただし、ｑは、駆動雑音のパワースペクトラム密度である。上記式（１３）および下記式（１４）から、トラッキングインデックスΓが決定されれば、αおよびβが一意に決定されることが分かる。
0006952944.tif 000007
【００１８】
図２は、平滑フィルタ装置１による平滑処理を示すイメージ図である。観測誤差標準偏差σｖ，ｋは、観測方法の観測精度に相当するパラメータである。また、観測タイミングは、観測時刻の時間間隔（サンプリング間隔）Ｔｋである。例えば、観測方法（１）における観測誤差標準偏差がσｖ，ｋ（１）であり、サンプリング間隔がＴ（１）である。観測方法（２）における観測誤差標準偏差がσｖ，ｋ（２）であり、サンプリング間隔がＴ（２）である。σｖ，ｋ（１）とσｖ，ｋ（２）は互いに異なり、Ｔ（１）とＴ（２）は互いに異なる。また、図２において、観測方法（１）が観測方法（２）よりも観測精度が高いものとする。
【００１９】
平滑フィルタ装置１に対して、観測方法（１）と観測方法（２）とによって共通の観測対象が観測された観測値が混在した時系列データが入力された場合に、図２の左側に示すように、観測方法（１）は、観測方法（２）よりも観測精度が高いため、丸印のプロットで示す観測値が、真値に近い値となっている。一方、観測精度が低い観測方法（２）によって観測された観測値（四角印のプロットで示す）は、真値から外れる度合いが大きくなっている。観測方法（１）と観測方法（２）では、観測タイミングが異なるため、観測値間の時間間隔も時々刻々とばらばらになっている。
【００２０】
観測方法（１）と観測方法（２）によって観測対象が観測された観測値が混在した時系列データを平滑処理する場合、従来の平滑フィルタには、観測誤差標準偏差（観測方法の観測精度に相当する）は、観測方法ごとに設定されるが、観測対象の運動モデルに関連した駆動雑音標準偏差は、観測方法によらずに設定される。このため、従来の平滑フィルタでは、観測方法（１）と観測方法（２）とによって共通の観測対象が観測された観測値が混在した時系列データを平滑処理するときに、観測方法の違いに応じた駆動雑音共分散行列Ｑｋの調整が行われず、平滑性能が観測方法ごとに異なって観測値を均一に平滑化できない可能性があった。
【００２１】
図３は、ランダム雑音に対する平滑処理を示すイメージ図である。図３において、Ｎ１は、平滑処理前のランダム雑音であり、Ｎ２は、平滑フィルタ装置１による平滑処理後のランダム雑音である。図３に示すようなランダム雑音Ｎ１からランダム雑音Ｎ２へ平滑化する平滑フィルタの平滑性能の指標として、定常状態における平滑値の雑音抑圧比Ｓ（α，β）がある。雑音抑圧比Ｓ（α，β）は、下記式（１５）から算出することができる。
0006952944.tif 000008
【００２２】
上記式（１５）に示すように、雑音抑圧比Ｓ（α，β）は、αおよびβのみの関数である。よって、トラッキングインデックスΓが決定されれば、雑音抑圧比Ｓ（α，β）は、一意に決定される。さらに、雑音抑圧比Ｓ（α，β）が決定されると、平滑誤差標準偏差σ平滑誤差は、雑音抑圧比Ｓ（α，β）の平方根を用いて下記式（１６）から算出することができる。ここで、σｖは、観測誤差標準偏差である。
0006952944.tif 000009
【００２３】
初期値設定部１０は、トラッキングインデックスΓについて観測方法（１）と観測方法（２）で平滑誤差標準偏差σ平滑誤差が同じ値となるように、初期値Γ１，０およびΓ２，０を算出し、算出した初期値Γ１，０およびΓ２，０を第２の計算部１２に設定する。Γ１，０は、観測方法（１）におけるトラッキングインデックスの初期値であり、Γ２，０は、観測方法（２）におけるトラッキングインデックスの初期値である。なお、初期値設定部１０は、平滑フィルタ装置１とは別に設けられた外部装置が備える構成要素とすることができる。この場合、平滑フィルタ装置１が備える第２の計算部１２は、外部装置が備える初期値設定部１０から、トラッキングインデックスの初期値が設定される。
【００２４】
第１の計算部１１は、観測方法（１）および観測方法（２）によって共通の観測対象が観測された観測値が混在した時系列データを入力し、観測値の観測時刻の時間間隔Ｔｋを算出する。第２の計算部１２は、トラッキングインデックスΓを用いて、観測方法（１）と観測方法（２）の平滑誤差標準偏差σ平滑誤差が同じ値となるように、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを算出する。
【００２５】
第２の計算部１２は、時刻ｔｋにおけるトラッキングインデックスΓｋと初期値Γ１，０またはΓ２，０が一致するように、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを算出する。平滑フィルタ装置１は、第２の計算部１２によって算出された駆動雑音標準偏差σｗ，ｋを用いることで、観測方法（１）と観測方法（２）によって共通の観測対象が観測された観測値が混在した時系列データの平滑処理において、観測方法（１）による観測値に対する平滑性能と観測方法（２）による観測値に対する平滑性能とを一致させることができる。
【００２６】
次に、初期値Γ１，０およびΓ２，０の設定処理について詳細に説明する。
図４は、初期値Γ１，０およびΓ２，０の設定処理を示すフローチャートである。図４において、平滑フィルタ装置１には、観測精度および観測タイミングが互いに異なる観測方法（１）と観測方法（２）によって共通の観測対象が観測された観測値が混在した時系列データが入力される。観測方法（１）および観測方法（２）の観測時刻の時間間隔Ｔ０をそれぞれ１．０（秒）とする。この場合、観測方法（１）と観測方法（２）によって観測対象が観測された観測値が混在した時系列データにおいては、観測時刻の時間間隔Ｔｋ（＝ｔｋ−ｔｋ−１）は、０（秒）から２（秒）までの間で変動する。
【００２７】
まず、初期値設定部１０は、観測方法（２）よりも観測精度が高い観測方法（１）に関するトラッキングインデックスΓの初期値Γ１，０を計算する（ステップＳＴ１）。例えば、初期値設定部１０には、予め決定された一定の平滑誤差標準偏差σ平滑誤差を導出するための、観測方法（１）に対応する既知の観測誤差標準偏差σｖ１と、観測方法（２）に対応する既知の観測誤差標準偏差σｖ２とが設定されている。
【００２８】
初期値設定部１０は、観測方法（１）に対応する駆動雑音標準偏差σｗ１，０を決定する。例えば、初期値設定部１０は、下記式（１７）で表される平滑誤差標準偏差を用いて、上記式（１５）および（１４）から駆動雑音標準偏差の初期値σｗ１，０を決定する。
0006952944.tif 000010
【００２９】
また、観測誤差標準偏差σｖ１と時間間隔Ｔ０＝１．０（秒）を用いた平滑処理のシミュレーションを行って平滑誤差標準偏差を評価することにより、適切な平滑誤差標準偏差に対応する駆動雑音標準偏差σｗ１，０を決定してもよい。平滑誤差標準偏差が観測誤差標準偏差σｖ１に比べて小さければ、平滑効果を与える適切な値であると言える。ただし、平滑誤差標準偏差が観測誤差標準偏差σｖ１に比べて著しく小さい場合は、平滑フィルタ装置１における観測対象への追従性能が劣化する。このため、許容可能な範囲内で平滑誤差標準偏差を決定する必要がある。
【００３０】
なお、実際のゲインα１およびβ１は解析的に得られない。このため、平滑フィルタ装置１において、上記式（１０）から（１１）までに示したカルマンフィルタの更新式を、上記式（１２）に示した定常ゲインＫの値が収束するまで繰り返し実行し、収束した定常ゲインＫの値から、ゲインα１およびβ１を決定する。
【００３１】
初期値設定部１０は、決定した駆動雑音標準偏差の初期値σｗ１，０と、既知の観測誤差標準偏差σｖ１と時間間隔Ｔ０（＝１．０（秒））を用いて、下記式（１８）に従い、初期値Γ１，０を算出する。
0006952944.tif 000011
【００３２】
次に、初期値設定部１０は、観測方法（１）に対応するＳ（α１，β１）σｖ１２を算出する（ステップＳＴ２）。観測方法（１）と観測方法（２）にそれぞれ対応する平滑誤差標準偏差が同じ値となるためには、下記式（１９）が成立する必要がある。
0006952944.tif 000012
【００３３】
初期値設定部１０は、下記式（２０）および（２１）に示す最適化問題を解くことで、駆動雑音標準偏差σｗ２，０を算出する。なお、下記の最適化問題は、解析的には解けないため、最小値σｗ２，０，ｍｉｎから最大値σｗ２，０，ｍａｘまでのサーチ範囲を一定値のステップΔσｗ２，０ごとにサーチすることで、観測方法（２）に対応する駆動雑音標準偏差σｗ２，０を決定する。
0006952944.tif 000013
【００３４】
まず、初期値設定部１０は、駆動雑音標準誤差σｗ２，０に最小値σｗ２，０，ｍｉｎを設定する（ステップＳＴ３）。続いて、初期値設定部１０は、ゲインα１およびβ１の場合と同様の手順で、定常ゲインＫの値からゲインα２およびβ２を決定し、ゲインα２およびβ２と既知の観測誤差標準偏差σｖ２を用いて、上記式（１５）および（１７）に従い、観測方法（２）に対応するＳ（α２，β２）σｖ２，０２を算出する（ステップＳＴ４）。次に、初期値設定部１０は、算出した値で上記式（２１）を満たすか否かを確認する（ステップＳＴ５）。
【００３５】
上記式（１９）を満たさない場合（ステップＳＴ５；ＮＯ）、初期値設定部１０は、現在のσｗ２，０にステップΔσｗ２，０を加算した値を、次の駆動雑音標準偏差σｗ２，０に設定する（ステップＳＴ６）。このとき設定された駆動雑音標準偏差σｗ２，０を用いて、ステップＳＴ４からの処理が実行される。これらの処理は上記サーチ範囲内で繰り返される。
【００３６】
上記式（１９）を満たす場合（ステップＳＴ５；ＹＥＳ）、初期値設定部１０は、このときの駆動雑音標準偏差σｗ２，０を採用する。この駆動雑音標準偏差σｗ２，０では、上記式（１９）が成立するので、上記式（２１）から、Ｊ（σｗ２，０）は、ほぼ零になる。
【００３７】
初期値設定部１０は、決定した駆動雑音標準偏差の初期値σｗ２，０と、観測誤差標準偏差σｖ２と、時間間隔Ｔ０＝１．０（秒）を用いて、下記式（２２）に従い、初期値Γ２，０を算出する（ステップＳＴ７）。初期値Γ１，０と初期値Γ２，０では上記式（１９）が成立するので、観測方法（１）と観測方法（２）では平滑誤差標準偏差が同じ値である。初期値設定部１０は、初期値Γ１，０および初期値Γ２，０を、第２の計算部１２に設定する。
0006952944.tif 000014
【００３８】
次に、実施の形態１に係る平滑化方法について説明する。
図５は、実施の形態１に係る平滑化方法を示すフローチャートである。図５において、平滑フィルタ装置１には、観測方法（１）と観測方法（２）とによって共通の観測対象が観測された観測値が混在した時系列データが入力されて、第２の計算部１２には、初期値設定部１０によって算出された初期値Γ１，０および初期値Γ２，０が設定されている。
【００３９】
第１の計算部１１は、観測方法（１）と観測方法（２）によって共通の観測対象が観測された観測値が混在した時系列データが入力されると、観測時刻の時間間隔Ｔｋを順次算出する（ステップＳＴ１ａ）。なお、時間間隔Ｔｋ（＝ｔｋ−ｔｋ−１）は、時刻ｔｋと時刻ｔｋ−１とで観測方法が異なっていてもよい。
【００４０】
第２の計算部１２は、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを更新する（ステップＳＴ２ａ）。第２の計算部１２は、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを、下記式（２３）または下記式（２４）に従って算出する。下記式（２３）には、時刻ｔｋにおける観測誤差標準偏差σｖ，ｋとして、既知のσｖ１が設定され、下記式（２４）には、時刻ｔｋにおける観測誤差標準偏差σｖ，ｋとして、既知のσｖ２が設定されている。
0006952944.tif 000015
【００４１】
第２の計算部１２によって下記式（２３）に従って算出された駆動雑音標準偏差σｗ，ｋおよび上記式（１４）を用いて、下記式（２５）の関係が導出される。下記式（２５）において、トラッキングインデックスΓ１，ｋは、時刻ｔｋにおける、観測方法（１）に対応したトラッキングインデックスである。下記式（２５）が示すように、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋは、トラッキングインデックスΓ１，ｋと初期値Γ１，０が一致するように算出されている。また、下記式（２５）と同様な手順で、下記式（２６）の関係が導出される。下記式（２６）において、トラッキングインデックスΓ２，ｋは、時刻ｔｋにおける、観測方法（２）に対応したトラッキングインデックスである。これにより、平滑フィルタ装置１において、各時刻ｔｋにおける平滑誤差標準偏差が一定に保たれる。
0006952944.tif 000016
【００４２】
予測部１３は予測処理を行う（ステップＳＴ３ａ）。例えば、予測部１３は、上記式（１０）に従い、次時刻ｔｋ＋１における観測対象の予測ベクトルｘｋ＋１バーおよび予測誤差共分散行列Ｐｋ＋１バーを算出する。続いて、平滑部１４は平滑処理を行う（ステップＳＴ４ａ）。例えば、平滑部１４は、上記式（１１）に従い、次時刻ｔｋ＋１における観測対象の平滑ベクトルｘｋ＋１ハットおよび平滑誤差共分散行列Ｐｋ＋１ハットを算出する。平滑部１４は、平滑ベクトルｘｋ＋１ハットおよび平滑誤差共分散行列Ｐｋ＋１ハットを、観測値の平滑結果として出力するとともに、予測部１３にフィードバックする。
【００４３】
次に、平滑フィルタ装置１の機能を実現するハードウェア構成について説明する。
平滑フィルタ装置１における初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の各機能は、処理回路により実現される。すなわち、平滑フィルタ装置１は、図５に示すステップＳＴ１ａからステップＳＴ４ａまでの処理を実行するための処理回路を備える。処理回路は、専用のハードウェアであってもよいが、メモリに記憶されたプログラムを実行するＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）であってもよい。
【００４４】
図６Ａは、平滑フィルタ装置１の機能を実現するハードウェア構成を示すブロック図であり、図６Ｂは、平滑フィルタ装置１の機能を実現するソフトウェアを実行するハードウェア構成を示すブロック図である。図６Ａおよび図６Ｂにおいて、入力インタフェース１００は、センサから平滑フィルタ装置１へ入力される観測値の時系列データを中継するインタフェースである。出力インタフェース１０１は、平滑フィルタ装置１から外部へ出力される平滑結果の情報を中継するインタフェースである。
【００４５】
処理回路が図６Ａに示す専用のハードウェアの処理回路１０２である場合、処理回路１０２は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。平滑フィルタ装置１における初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の機能を、別々の処理回路で実現してもよく、これらの機能をまとめて１つの処理回路で実現してもよい。
【００４６】
処理回路が図６Ｂに示すプロセッサ１０３である場合、平滑フィルタ装置１における初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の機能は、ソフトウェア、ファームウェアまたはソフトウェアとファームウェアとの組み合わせにより実現される。なお、ソフトウェアまたはファームウェアは、プログラムとして記述されてメモリ１０４に記憶される。
【００４７】
プロセッサ１０３は、メモリ１０４に記憶されたプログラムを読み出して実行することで、平滑フィルタ装置１における初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の機能を実現する。例えば、平滑フィルタ装置１は、プロセッサ１０３によって実行されるときに、図５に示したフローチャートにおけるステップＳＴ１ａからステップＳＴ４ａまでの処理が結果的に実行されるプログラムを記憶するためのメモリ１０４を備える。これらのプログラムは、初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の手順または方法をコンピュータに実行させる。メモリ１０４は、コンピュータを、初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４として機能させるためのプログラムが記憶されたコンピュータ可読記憶媒体であってもよい。
【００４８】
メモリ１０４は、例えば、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ−ＥＰＲＯＭ）などの不揮発性または揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤなどが該当する。
【００４９】
平滑フィルタ装置１における初期値設定部１０、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４の機能の一部を専用のハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現してもよい。例えば、初期値設定部１０は、専用のハードウェアである処理回路１０２により機能を実現し、第１の計算部１１、第２の計算部１２、予測部１３および平滑部１４は、プロセッサ１０３がメモリ１０４に記憶されたプログラムを読み出して実行することにより機能を実現する。このように、処理回路は、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせによって、上記機能を実現することができる。
【００５０】
次に、航空監視システム（Ｍｕｌｔｉｌａｔｅｌａｔｉｏｎ；ＭＬＡＴ）に適用された平滑フィルタ装置１について説明する。ＭＬＡＴは、空港面に設置された複数のセンサによって航空機から送信された信号を同期受信し、複数のセンサによって受信された信号間の到来時間差（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ；ＴＤＯＡ）に基づいて、空港面に存在する航空機の測位を行う。ただし、実際のセンサ間には時刻同期誤差（以下、クロックオフセットと呼ぶ）が存在するために、クロックオフセットの値を予め観測して推定しておく必要がある。クロックオフセットの観測方法には、基準局方法と、ＧＰＳ（Ｇｌｏｂａｌ  Ｐｏｓｉｔｉｏｎｉｎｇ  Ｓｙｓｔｅｍ）方法が存在する。
【００５１】
基準局方法は、空港面に分散配置された、位置が既知の基準局から送信された信号を、位置が既知の複数の受信局（センサに相当する）によって受信し、受信局間の受信時刻の相対的なずれを、クロックオフセットとして観測する方法である。また、ＧＰＳ方法は、複数の受信局（センサに相当する）のそれぞれによってＧＰＳ衛星からＧＰＳパルス信号を受信し、１ＰＰＳ（Ｐｕｌｓｅ  Ｐｅｒ  Ｓｅｃｏｎｄ）ごとに得られるＧＰＳ時刻に基づいて、受信局間のＧＰＳ時刻の相対的なずれを、クロックオフセットとして観測する方法である。
【００５２】
図７は、ＭＬＡＴにおける航空機３０の測位の概要を示すイメージ図である。空港面２０に設置された受信局２１〜２４は、航空機３０から送信された信号を同期受信し、受信局間の受信タイミングの時刻差を用いた下記式（２７）を解くことにより、航空機３０の測位が行われる。なお、空港面２０において、受信局２１〜２４の各位置は既知である。下記式（２７）において、ＴＤＯＡｉｊは、第ｉ番目の受信局に到来した信号と第ｊ番目の受信局に到来した信号の到来時間差である。ｘは、測位対象の航空機３０の位置ベクトルである。ｐｉは、第ｉ番目の受信局の位置ベクトルであり、既知の値である。例えば、三次元測位を行う場合には、少なくとも３つのＴＤＯＡが必要であることから、４局以上の受信局が必要となる。
0006952944.tif 000017
【００５３】
ＴＤＯＡの観測には、各受信局によって信号が受信された時刻における各受信局の時計（クロック）の値が参照される。しかしながら、受信局の時計は実際には同期していないため、何らかの方法で非同期成分を推定し、推定した非同期成分を用いてＴＤＯＡを補償する必要がある。
【００５４】
図８は、ＭＬＡＴにおける受信局の時刻管理の概要を示すイメージ図である。例えば、基準局方法では、基準局２５から送信された同一の信号（電波）を、受信局２６〜２８で受信する。受信局２６〜２８がそれぞれ備える時計が同期していない場合、下記式（２８）に示すように、受信局が備える時計を用いて計測された上記信号の到来時間差と、基準局２５から各受信局までの距離から見積もられた電波の到来時間差との間には、差Δが生じる。基準局方法は、この差Δをクロックオフセットとして観測する方法である。下記式において、Δｉｊ，基準局は、基準局方法における観測対象である、第ｉ番目の受信局と第ｊ番目の受信局との間のクロックオフセット観測値であり、Ｔｉ，基準局は、第ｉ番目の受信局が備える時計によって計測された上記信号の受信時刻である。
0006952944.tif 000018
【００５５】
ＧＰＳ方法では、下記式（２９）に示すように、１ＰＰＳごとのＧＰＳパルス信号によって計測されるＧＰＳ時刻と、受信局が備える時計によって計測された時刻との差Δが、クロックオフセットとして観測される。下記式において、Δｉｊ，ＧＰＳは、ＧＰＳ方法における観測対象である、第ｉ番目の受信局と第ｊ番目の受信局との間のクロックオフセット観測値であり、Ｔｉ，ＧＰＳは、第ｉ番目の受信局が備える時計によって計測された上記信号の受信時刻である。
0006952944.tif 000019
【００５６】
Δｉｊ，基準局およびΔｉｊ，ＧＰＳは、観測精度および観測タイミングが互いに異なる基準局方法とＧＰＳ方法とによって観測された観測値である。平滑フィルタ装置１には、これらの観測値が混在した時系列データが入力される。平滑フィルタ装置１は、トラッキングインデックスΓを用いて、基準局方法に対応する平滑誤差標準偏差とＧＰＳ方法に対応する平滑誤差標準偏差とが同じ値となるように、時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを算出する。平滑フィルタ装置１は、算出した時刻ｔｋにおける駆動雑音標準偏差σｗ，ｋを用いて、観測値Δｉｊ，基準局およびΔｉｊ，ＧＰＳが混在した時系列データの平滑処理において、基準局方法による観測値に対する平滑性能とＧＰＳ方法による観測値に対する平滑性能とを一致させることができる。
【００５７】
図９は、クロックオフセットの観測時刻と観測値の平滑誤差絶対値との関係を示すグラフである。図９に示す関係は、基準局方法とＧＰＳ方法によるクロックオフセット観測をそれぞれシミュレーションし、これらのシミュレーションによって得られた観測値が混在した時系列データを平滑処理して得たものである。平滑処理は、平滑フィルタ装置１および従来の平滑フィルタによって行っている。クロックオフセット観測のシミュレーションにおいて、基準局方法における観測精度、すなわち観測誤差標準偏差σｖ１を１（ナノ秒）とし、ＧＰＳ方法における観測誤差標準偏差σｖ２を１５（ナノ秒）としており、ＧＰＳ方法の観測精度は、基準局方法よりも劣っている。
【００５８】
図９において、四角のプロットによって平滑誤差絶対値を示す関係Ａが、従来の平滑フィルタ装置による平滑処理の結果である。丸形のプロットによって平滑誤差絶対値を示す関係Ｂは、平滑フィルタ装置１による平滑処理の結果である。従来の平滑フィルタでは、基準局方法とＧＰＳ方法とで平滑性能を一致させる処理が行われない。このため、従来の平滑フィルタの平滑処理は、平滑フィルタ装置１の平滑処理よりも平滑誤差が大きくなっている。
【００５９】
また、図１０は、クロックオフセットの観測値の平滑処理をモンテカルロ試行した結果を示すグラフであり、図９に示したクロックオフセット観測のシミュレーションによって得られた観測値に対する平滑処理についてモンテカルロシミュレーションを１００回試行したときの各試行における平均二乗誤差（ＲＭＳＥ）を示している。図１０において、四角のプロットによってＲＭＳＥを示す関係Ｃが、従来の平滑フィルタ装置による平滑処理の結果である。丸形のプロットによってＲＭＳＥを示す関係Ｄが、平滑フィルタ装置１による平滑処理の結果である。図１０に示すように、平滑フィルタ装置１の平滑処理は、従来の平滑フィルタの平滑処理よりもＲＭＳＥが小さいと言える。
【００６０】
次に、バイスタティックレーダシステムに平滑フィルタ装置１を適用した場合について説明する。バイスタティックレーダシステムは、レーダ波の送信機と受信機の間を離して配置し、両者を通信で結んだシステムであり、送信機と受信機との間の通信は高度に同期が取られていることが求められる。バイスタティックレーダシステムでは、送信機と受信機との通信の同期を取るため、一般的に、原子時計（例えば、ルビジウム時計）を用いた絶対時刻の精度管理が行われる。原子時計における時刻計測は、ＧＰＳにおける時刻計測よりも計測精度が低く、原子時計によって計測される時刻は、徐々にドリフトして真値から外れることがある。このため、原子時計によって計測される時刻は、ＧＰＳ時刻を用いて修正をかける必要がある。
【００６１】
バイスタティックレーダシステムにおいて、平滑フィルタ装置１は、トラッキングインデックスを用いて、原子時計による時刻計測とＧＰＳにおける時刻計測とで平滑誤差標準偏差とが同じ値となるように駆動雑音標準偏差σｗ，ｋを算出し、算出した駆動雑音標準偏差σｗ，ｋを用いて平滑処理を行う。これにより、原子時計によって計測された時刻の値とＧＰＳにおいて計測された時刻の値とが混在した時系列データの平滑処理において、原子時計によって計測された時刻に対する平滑性能とＧＰＳによって計測られた時刻に対する平滑性能とを一致させることができる。
【００６２】
これまで、観測精度および観測タイミングが互いに異なる２つの観測方法によって観測された観測値が混在した時系列データについて説明したが、平滑フィルタ装置１は、３つ以上の観測方法によって観測された観測値が混在した時系列データを平滑処理することができる。すなわち、平滑フィルタ装置１には、３つ以上の観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように算出されたトラッキングインデックスΓの初期値が設定される。平滑フィルタ装置１は、この初期値を用いて３つ以上の観測方法に対応する平滑性能を一致させることができる。
【００６３】
以上のように、実施の形態１に係る平滑フィルタ装置１において、トラッキングインデックスを用いて、複数の観測方法にそれぞれ対応する平滑誤差標準偏差が同じ値となるように、各観測時刻における駆動雑音標準偏差が算出され、算出された駆動雑音標準偏差を用いて観測値の時系列データの平滑処理が行われる。観測精度および観測タイミングが互いに異なる複数の観測方法によって共通の観測対象が観測された観測値が混在した時系列データの平滑処理において、複数の観測方法にそれぞれ対応する平滑性能を一致させることができる。
【００６４】
なお、本発明は上記実施の形態に限定されるものではなく、本発明の範囲内において、実施の形態の任意の構成要素の変形もしくは実施の形態の任意の構成要素の省略が可能である。
【産業上の利用可能性】
【００６５】
本発明に係る平滑フィルタ装置は、例えば、ＭＬＡＴに利用可能である。
【符号の説明】
【００６６】
１  平滑フィルタ装置、１０  初期値設定部、１１  第１の計算部、１２  第２の計算部、１３  予測部、１４  平滑部、２０  空港面、２１  受信局、２２，２３，２４，２６，２７，２８  受信局、２５  基準局、３０  航空機、１００  入力インタフェース、１０１  出力インタフェース、１０２  処理回路、１０３  プロセッサ、１０４  メモリ。
【図１】
0006952944.tif 000020
【図２】
0006952944.tif 000021
【図３】
0006952944.tif 000022
【図４】
0006952944.tif 000023
【図５】
0006952944.tif 000024
【図６】
0006952944.tif 000025
【図７】
0006952944.tif 000026
【図８】
0006952944.tif 000027
【図９】
0006952944.tif 000028
【図１０】
0006952944.tif 000029