JP5424338B2 - 衛星測位システムの異常値検出装置、異常値検出方法及び異常値検出プログラム - Google Patents

Description

本発明は衛星測位システムの異常値検出装置、異常値検出方法及び異常値検出プログラムに係り、特に測位信号から算出される擬似距離に含まれる異常値を検出するＧＢＡＳ等の衛星測位システムの異常値検出装置、異常値検出方法及び異常値検出プログラムに関する。

人工衛星が送信する測位信号を受信機で受信し、その測位信号から受信機自身の位置を知ることができる全地球測位システム（ＧＰＳ：Global Positioning System）などの衛星測位システムは、車両、航空機、または船舶などの航法支援などに広く利用されている。この衛星測位システムにおける受信機は、測位信号から人工衛星との距離を算出し、その距離を基に自身の位置を算出する。この測位信号から算出される受信機と人工衛星との間の距離は、擬似距離と呼ばれる。

上記の衛星測位システムでは、各人工衛星と受信機との間の擬似距離に基づいて測位計算を行うので、この擬似距離の精度は、測位精度に大きな影響を与える。この擬似距離には、時計誤差、電離層による誤差、対流圏による誤差、ノイズ等の様々な原因により生じた異常値が含まれている。このため、この異常値を適切に取り除くことが、衛星測位システムの測位精度を高め、利用時の安定性を向上させるうえで重要である。

そこで、衛星測位システムにおいては、通常、逐次測定される擬似距離等の時系列データにおける異常値の発生状態を監視する異常値検出装置（ＩＭ：Integrity Monitor）が用いられることが多い。特に、地上型補強システム（ＧＢＡＳ：Ground Based Augmentation System）に代表されるような、衛星測位システムを航空機の航法に利用するシステムでは、高い安全性が要求されるため、異常値検出装置は不可欠である。

非特許文献１に記載された衛星測位システムの異常値検出装置は、被モニタ事象のモデルを前日以前の過去データを用いてモデル化している。この異常値検出装置は、モデルとして１０度の仰角幅毎にガウス分布を仮定し、この１０度の仰角幅毎に求めた分散の値を用いて各データを正規化する。そして、この異常値検出装置は、正規化したデータの分布が標準ガウス分布に従うと仮定し、この標準ガウス分布に従わない値を異常値として検出する。

また、特許文献１に記載された衛星測位システムの異常値検出装置は、所定の有効期間内における衛星軌道位置の値同士の差の最悪値と閾値とを比較し、最悪値が閾値を超えていた場合に、最悪値の元となった衛星軌道位置の値を異常値とする。

特開２００９−６８９２７号公報

Gang Xie,"OPTIMAL ON-AIRPORT MONITORING OF THE INTERGRITY OF GPS-BASE LANDING SYSTEMS",Ph.D. Dissertation, Stanford University,pp.26-32,March 2004.

しかし、非特許文献１や特許文献１に記載された異常値検出装置では、精度良く異常値を検出できない場合がある。

すなわち、非特許文献１や特許文献１に記載された異常値検出装置では、過去のデータに基づく不変のモデルを仮定しているため、観測環境が変化する場合、不変のモデルでは対応できないことがある。このような場合は、計算結果が信用できないものとなり得る。

また、これらの異常値検出装置では、モデル化において、各データ間に時系列的な相関を仮定していないため、時系列的な相関関係がある場合、その関係を捉えることができない。この場合も、計算結果に誤りが生じ得る。

本発明は以上の点に鑑みなされたもので、観測環境に変化が生じたり、各データ間に時系列的な相関関係がある場合であっても、精度良く異常値を検出し得るＧＢＡＳ等の衛星測位システムの異常値検出装置、異常値検出方法及び異常値検出プログラムを提供することを目的とする。

上記の目的を達成するため、本発明の衛星測位システムの異常値検出装置は、衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１確率分布取得手段と、第１の期間内の単位時刻である各時刻ｔについて第１確率分布取得手段により取得された第１の確率分布に基づき、時刻ｔの直前の時刻ｔ−１における第１の確率変数の選択情報量を不確実性指標として算出する不確実性指標算出手段と、時刻ｔ以前の、第１の期間より短い第２の期間内で、不確実性指標算出手段により算出された不確実性指標の平均値を算出し、その平均値を第２の確率変数として、第２の確率変数の重み付け統計量に基づき、第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第２確率分布取得手段と、第２確率分布取得手段により各時刻ｔ毎に取得された第２の確率分布に基づき、第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の第２の期間内の第２の確率分布の平均情報量を変化点指標として算出する変化点指標算出手段と、変化点指標算出手段により算出された変化点指標を予め設定された閾値と比較し、閾値より大きな値の変化点指標に対応する時系列値を異常値として検出する異常値検出手段とを有することを特徴とする。

また、上記の目的を達成するため、本発明の衛星測位システムの異常値検出方法は、衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、第１の期間内の単位時刻である各時刻ｔについて第１のステップにより取得された第１の確率分布に基づき、時刻ｔの直前の時刻ｔ−１における第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、時刻ｔ以前の、第１の期間より短い第２の期間内で、第２のステップにより算出された不確実性指標の平均値を算出する第３のステップと、第３のステップで算出された平均値を第２の確率変数として、第２の確率変数の重み付け統計量に基づき、第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、第４のステップにより各時刻ｔ毎に取得された第２の確率分布に基づき、第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の第２の期間内の第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、第５のステップにより算出された変化点指標を予め設定された閾値と比較し、閾値より大きな値の変化点指標に対応する時系列値を異常値として検出する第６のステップとを含むことを特徴とする。

更に、上記の目的を達成するため、本発明の衛星測位システムの異常値検出プログラムは、コンピュータに、
衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、第１の期間内の単位時刻である各時刻ｔについて第１のステップにより取得された第１の確率分布に基づき、時刻ｔの直前の時刻ｔ−１における第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、時刻ｔ以前の、第１の期間より短い第２の期間内で、第２のステップにより算出された不確実性指標の平均値を算出する第３のステップと、第３のステップで算出された平均値を第２の確率変数として、第２の確率変数の重み付け統計量に基づき、第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、第４のステップにより各時刻ｔ毎に取得された第２の確率分布に基づき、第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の第２の期間内の第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、第５のステップにより算出された変化点指標を予め設定された閾値と比較し、閾値より大きな値の変化点指標に対応する時系列値を異常値として検出する第６のステップとを実行させることを特徴とする。

本発明によれば、時系列値を確率変数として第１の確率密度関数を取得し、その第１の確率密度関数から算出した選択情報量である不確実性指標の平均値を確率変数として第２の確率密度関数を算出し、更にその第２の確率密度関数の平均情報量である変化点指標を算出し、閾値より大きい変化点指標に対応する時系列値を異常値として検出するため、各データ間に時系列的な相関関係のある場合や、観測環境の変化が生じた場合であっても、精度の高い異常値の検出ができる。

本発明の衛星測位システムの異常値検出装置の一実施形態のブロック図である。 図１中の動的モデル構成部の一実施形態のブロック図である。 本発明の衛星測位システムの異常値検出方法の一実施形態のフローチャートである。 時系列で観測される加速度値の時系列データに、意図的に混入される異常値の値と、その異常値を混入する時間の一例を示す図である。 図１の異常値検出装置による、図４に示した異常値が意図的に混入された加速度値の時系列データの異常値検出結果を示す図である。 本発明の衛星測位システムの異常検出装置が適用される地上型補強システム（ＧＢＡＳ）の一実施例のシステム構成図である。

次に、本発明の実施形態について図面を参照して説明する。

図１は、本発明になる衛星測位システムの異常値検出装置の一実施形態のブロック図を示す。同図に示すように、本実施形態の異常値検出装置１０は、異常値指標算出部１１、動的モデル構成部１２及び異常値検出部１３より構成され、入力される時系列データを監視し、衛星測位システムの異常値を検出する。

異常値指標算出部１１は、時系列データが入力され、その時系列データの各時刻における異常値指標を算出する。上記の時系列データは、衛星測位システムの各人工衛星と受信機との間の擬似距離や受信測位信号の位相などの、受信機において人工衛星からの測位信号に基づいて時系列に取得されたデータである。異常値指標とは、異常値の検出対象となる時系列値である。ここでは、異常値指標は、擬似距離又は位相の加速度等であり、検出したい異常に鋭敏に反応する値である。この異常値指標は、衛星測位システムの異常値の検出対象である本発明の時系列値に相当する。

動的モデル構成部１２は、一定期間において、異常値指標算出部１１により時系列に求められた異常値指標から動的にモデルを構成する。また、動的モデル構成部１２は、構成した動的モデルに基づいて、時系列の異常値指標から変化点指標を算出する。変化点指標は、突発的に増減した時系列値が存在する場合、その時系列値が単発的な動的モデルからの外れ値であるのか、それとも入力データの動的モデル自体が変化しているのかを判断する指標である。本実施形態において、この変化点指標は、動的モデル自体が変化している可能性が高いほど、高くなる値とされ、ある閾値より大きい変化点指標が算出された場合、その時刻を変化点として扱う。

図２は、動的モデル構成部１２の一実施形態のブロック図を示す。同図に示すように、動的モデル構成部１２は、第１確率分布取得部１２１、不確実性指標算出部１２２、第２確率分布取得部１２３、及び変化点指標算出部１２４より構成される。

第１確率分布取得部１２１は、衛星測位システムの異常値の検出対象である時系列値（すなわち、異常値指標）を第１の確率変数として、その異常値指標の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して、それを第１の確率分布として取得する。

不確実性指標算出部１２２は、上記第１の期間内の単位時刻である各時刻ｔについて、第１確率分布取得部１２１により取得された第１の確率分布に基づき、時刻ｔ−１に対応する第１の確率変数の選択情報量を不確実性指標として算出する。

第２確率分布取得部１２３は、各時刻ｔについて、上記第１の期間よりも短い期間を第２の期間として、時刻ｔ以前の第２の期間内で不確実性指標算出部１２２により算出された不確実性指標の平均値を算出し、その平均値を第２の確率変数として、第２の確率変数の重み付け統計量に基づき、上記第１の期間における第２の確率密度関数を算出して、それを第２の確率分布として取得する。

変化点指標算出部１２４は、第２確率分布取得部１２３により各時刻ｔ毎に取得された第２の確率分布に基づき、第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の第２の期間内の第２の確率分布（すなわち、第２の確率密度関数）の平均情報量を変化点指標として算出する。この変化点指標の算出において、動的モデル構成部１２は、その計算のための母集団に時系列モデルを導入している。このため、動的モデル構成部１２は、変化点指標の算出において、異常値指標の時系列な挙動を反映させることができる。

図１に示す異常値検出部１３は、動的モデル構成部１２により各時刻ｔ毎に算出された変化点指標を予め設定した閾値と比較し、変化点指標の値が閾値よりも大きければ、その時刻の変化点指標に対応する異常値指標（時系列値）を異常値として検出する。

次に、図１及び図２に示した本実施形態の異常値検出装置１０の動作について、図３のフローチャート等を併せ参照して詳細に説明する。

まず、異常値検出装置１０は、異常値指標算出部１１により、入力された時系列データから各時刻の異常値指標を算出する（ステップＳ１）。続いて、異常値検出装置１０は第１確率分布取得部１２１により、ステップＳ１で算出された異常値指標を第１の確率変数として、現時刻以前の所定の第１の期間Ｔ１の期間内における、その第１の確率変数の第１の確率密度関数ｐ(ｘ_t)を次式により求める（ステップＳ２）。

（１）式において、ｔは各異常値指標（第１の確率変数）が算出された時刻であり、ｘ_tは時刻ｔにおける異常値指標（第１の確率変数）の値である。また、（１）式において、ｗ１_t-1は、期間Ｔ１のうち、時刻ｔ−１以前の期間における、異常値指標（第１の確率変数）ｘ_tの重み付け平均値である。この重み付け平均値ｗ１_t-1は、現時刻に近いほど大きな値になるように重み付けを行った異常値指標（第１の確率変数）ｘ_tの平均値である。更に、（１）式において、σ１_t-1は、期間Ｔ１のうち、時刻ｔ−１以前の期間における、異常値指標（第１の確率変数）ｘ_tの重み付け分散値である。この重み付け分散値σ１_t-1は、現時刻に近いほど大きな値になるように重み付けを行った異常値指標（第１の確率変数）ｘ_tの分散値である。

続いて、異常値検出装置１０は、不確実性指標算出部１２２により、第１確率分布取得部１２１から入力される時刻ｔ−１における第１の確率密度関数に基づいて、次式によりｔ＝２以降の各時刻ｔにおける不確実性指標ＵＬＩ(ｔ)を算出する（ステップＳ３）。

時刻ｔの不確実性指標ＵＬＩ(ｔ)は、（２）式から時刻ｔの直前の時刻ｔ−１における異常値指標の選択情報量である。上記（２）式の右辺のｐ_t-1(ｘ_t|ｘ^t-1)は、時刻ｔ−１において（１）式から算出された第１の確率密度関数である。

すなわち、このｐ_t-1(ｘ_t|ｘ^t-1)は、ｘ^t-1(＝ｘ₁,ｘ₂,・・・,ｘ_t-1；時刻１からｔ−１までのデータ)に基づいて推定したパラメータθ_t-1で表される確率密度分布である。つまり、ｐ_t-1(ｘ_t|ｘ^t-1)＝ｐ_t-1(ｘ_t|θ_t-1)である。因みに、θ_t-1＝（ｗ1_t-1,σ1_t-1）である。この（２）式の右辺は、時刻ｔ＝２以降の各時刻ｔにおいて、時刻ｔ−１における確率密度関数ｐ_t-1(ｘ_t|ｘ^t-1)の対数のマイナス値、すなわち選択情報量を算出することを意味する。

次に、異常値検出装置１０は、第２確率分布取得部１２３により、不確実性指標算出部１２２により算出された不確実性指標ＵＬＩ(ｔ)に基づき、次式により各時刻ｔについて、時刻ｔ以前の期間Ｔ２における不確実性指標の平均値ｙ_tを算出する（ステップＳ４）。ここで、期間Ｔ２の長さは、期間Ｔ１より短く設定される。

（３）式において、ＵＬＩ(ｉ)は、時刻ｉに対応する不確実性指標を示す。

第２確率分布取得部１２３は、続いて、上記の不確実性指標の平均値ｙ_tを第２の確率変数として、次式に基づき時刻ｔ以前の期間Ｔ３における確率変数の第２の確率密度関数ｑ(ｙ_t)を算出する（ステップＳ５）。この第２の確率密度関数ｑ(ｙ_t)は第２の確率分布として取得される。ここで、期間Ｔ３の長さは、期間Ｔ１以下に設定される。

（４）式において、ｗ２_t-1は期間Ｔ３のうち、時刻ｔ−１以前の期間における第２の確率変数（不確実性指標の平均値）ｙ_tの重み付け平均値である。この重み付け平均値ｗ２_t-1は、現時刻に近いほど値が大きくなるように重み付けを行ったｙ_tの平均値である。また、（４）式において、σ２_t-1は期間Ｔ３のうち、時刻ｔ−１以前の期間における第２の確率変数（不確実性指標の平均値）ｙ_tの重み付け分散値である。この重み付け分散値σ２_t-1は、現時刻に近いほど値が大きくなるように重み付けを行ったｙ_tの分散値である。

次に、異常値検出装置１０は、変化点指標算出部１２４により、第２確率分布取得部１２３により算出された第２の確率密度関数に基づき、次式により各時刻ｔにおける変化点指標ＣＰＩ(ｔ)を算出する（ステップＳ６）。

（５）式において、Ｔ４は時刻ｔ−１以前の期間であり、また、ｑ_i-1(ｙ_i|ｙ^i-1)は時刻ｉにおける第２の確率密度関数（第２の確率分布）である。（５）式の右辺は、各時刻ｔにおいて、時刻ｔ−１以前の期間Ｔ４の第２の確率密度関数ｑ_i-1(ｙ_i|ｙ^i-1)の対数のマイナス値の平均値、すなわち平均情報量を算出することを意味する。ここで、期間Ｔ４の長さは、期間Ｔ１より短く設定される。また、本実施形態では期間Ｔ４の長さは、期間Ｔ２以下に設定される。

異常値検出装置１０は、期間Ｔ２及びＴ４の長さが短いほど演算量が小さくなり、変化点指標を速く算出することができるが、異常値検出における誤検出率が高くなる。一方、期間Ｔ２及びＴ４の長さが長いほど演算量が多くなり、変化点指標の算出が遅くなるが、異常値検出装置１０は、高い精度で異常値を検出することができる。従って、上記の期間Ｔ２及びＴ４の長さは、変化点指標の算出速度（これは異常値の検出速度に対応する）と異常値の検出精度とのバランスを考慮して定められる。

最後に、異常値検出装置１０は、異常値検出部１３により、各時刻について変化点指標算出部１２４において算出された変化点指標と予め設定した閾値とを比較し、閾値より変化点指標の値が大きければその時刻の異常値指標を異常値として検出する（ステップＳ７）。

なお、上記の期間Ｔ１及びＴ３は、本発明における第１の期間に相当し、上記の期間Ｔ２及びＴ４は、本発明における第２の期間に相当する。

次に、図４及び図５を参照して、異常値検出装置１０が異常値を検出した結果の一例について説明する。図４及び図５では、入力された位相の加速度を異常値指標とする。ここで、人工衛星の軌道運動による自然な加速度の変化分は事前に排除している。

図４は、時系列で観測される加速度値の時系列データに、意図的に混入される異常値の値と、その異常値を混入する時間の一例を示す。同図において、縦軸は、意図的に混入される異常値の値を示し、横軸は、性能シミュレーションの開始時刻から経過した時間を示す。同図に斜線で示すように、異常値は、基準時から２７時間経過した時点から２８時間経過した時点までの期間において、徐々に大きな値となるように意図的に混入される。

図５は、本実施形態の異常値検出装置１０による、異常値が図４に示したように意図的に混入された加速度値の時系列データの異常値検出結果を示す。図５は、左側の縦軸が異常値指標を示し、右側の縦軸が変化点指標を示し、横軸が、加速度値の時系列データを取得した時間を示す。また、図５において、Iは本実施形態の異常値検出装置１０による異常値変換指標の時間変化、IIは変化点指標の時間変化を示す。

図５に示すように、本実施形態の異常値検出装置１０は、丸で囲んだ期間IIIにおいて、閾値（例えば「１.５」）以上の比較的高い値の変化点指標を算出している。この期間IIIは基準時から２７時間経過した時点であり、図４に示したように加速度値の時系列データに意図的に混入された異常値の混入時間と一致する。一方、図４に示したように、基準時から２７時間経過するまでの期間と、２８時間経過した後では、異常値は時系列データに混入されていないが、この期間において本実施形態の異常値検出装置１０により得られる変化点指標は図５にIIで示すように閾値以下であり、異常値は検出されない。従って、本実施形態の異常値検出装置１０によれば、時系列データに混入された異常値だけを精度良く検出できることが分かる。

以上説明したように、本実施形態の異常値検出装置１０は、時系列値を確率変数として、それぞれ現時刻に近いほど大きくなるように重み付けを行った重み付け平均値ｗ１_t-1及び重み付け分散値σ１_t-1を用いることで、直近の観測環境に即した第１の確率密度関数を算出し、その第１の確率密度関数から選択情報量である不確実性指標を算出する。そして、異常値検出装置１０は、更にその不確実性指標の平均値を確率変数として、それぞれ現時刻に近いほど大きくなるように重み付けを行った重み付け平均値ｗ２_t-1及び重み付け分散値σ２_t-1を用いることで、直近の観測環境に即した第２の確率密度関数を算出し、その第２の確率密度関数から平均情報量である変化点指標を算出し、閾値より大きい変化点指標に対応する時系列値を異常値として検出する。

これにより、異常値検出装置１０は、第１の確率密度関数の算出では捉えられない時系列的な相関関係を、第２の確率密度関数の算出により捉えることができ、また、第１の確率密度関数の算出では捉えられない観測環境の変化を、第２の確率密度関数の算出により捉えることができる。このため、異常値検出装置１０は、時系列的な相関関係のある場合や、観測環境の変化が生じた場合であっても、精度の高い異常値の検出ができる。

図６は、本発明になる異常検出装置が適用される地上型補強システム（ＧＢＡＳ）の一実施例のシステム構成図を示す。同図において、衛星測位システムである本実施例のＧＢＡＳは、データ処理装置１と、ＶＨＦ送信機３０と、基準局４１〜４４と、ＧＰＳを構成する人工衛星（ＧＰＳ衛星）５１〜５４と、ＧＢＡＳユーザ（航空機）６０とから構成されている。

データ処理装置１とＶＨＦ送信機３０と基準局４１〜４４とはＧＢＡＳ地上局を構成している。データ処理装置１は、上記の実施形態の異常値検出装置１０を有し、またディファレンシャル補正処理部２０を有する。

基準局４１〜４４は、ＧＰＳ衛星５１〜５４等から送信された測位信号を受信し、その受信信号であるＧＰＳ観測データをデータ処理装置１内のディファレンシャル補正処理部２０へそれぞれ出力する。ディファレンシャル補正処理部２０は、ＧＰＳ観測データを基に異常値検出装置１０の入力用に加工した統計情報量である時系列データを生成して異常値検出装置１０へ供給する。

異常値検出装置１０は、前述した実施形態の動作を行い、異常値を検出した時には検出した異常値をディファレンシャル補正処理部２０へ出力する。ディファレンシャル補正処理部２０は、異常値が入力されると該当するＧＰＳ衛星又は基準局を使用禁止にする。異常値が入力されない場合は、すべてのＧＰＳ衛星５１〜５４や基準局４１〜４４の組み合わせにより、ＧＰＳによる航法の精度や安全性を向上させるためのＧＢＡＳ補正データを公知の方法により生成してＶＨＳ受信機３０へ供給する。ＶＨＳ送信機３０は、入力されたＧＢＡＳ補正データで変調されたＶＨＳ信号を送信する。

ＧＢＡＳユーザである航空機６０は、ＧＰＳ衛星５１〜５４から送信された測位信号を受信して位置情報を取得するＧＰＳ受信機と、ＶＨＳ送信機３０から送信されたＧＢＡＳ補正データで変調されたＶＨＳ信号を受信する受信機（いずれも図示せず）を搭載している。航空機６０は、ＧＰＳ受信機の受信信号から取得した所定の情報と、ＶＨＳ信号の受信機の受信信号から復調したＧＢＡＳ補正データとを、公知の方法により着陸誘導等に用いることで、安全な航法を確保する。

なお、本発明は以上の実施形態及び実施例に限定されるものではなく、例えば、上記実施形態の異常値検出装置１０は、（１）式及び（４）式により計２回確率分布（確率密度関数）を求めているが、（４）式で算出された値を確率変数として、その確率変数の重み付け統計量に基づき更に確率密度関数を求め、その確率密度関数から変曲点指標を算出し、その変曲点指標と閾値との比較結果に基づいて、異常値を検知する構成としてもよい。

また、本発明は図３のフローチャートの全部又は一部のステップをコンピュータである中央処理装置（ＣＰＵ:Central Processing Unit）に実行させる異常値検出プログラムも包含する。この場合、異常値検出装置１０は、処理装置と、主記憶装置と、本発明の異常値検出プログラムを記憶した記憶装置とを備え、処理装置と主記憶装置により上記の異常値検出プログラムを実行することにより、図３のフローチャートの全部又は一部のステップを実現する構成とする。ここで、上記の処理装置には異常値検出プログラムに従って処理を実行するＣＰＵが設けられる。また、上記の主記憶装置は、図３のフローチャートで説明した処理に必要なデータ、処理装置による演算処理の過程で算出されるデータ及び演算処理の結果を保持する。

更に、本発明は人工衛星からの測位信号を利用する衛星測位システムであれば、ＧＢＡＳ以外の公知の他の衛星測位システムにも適用可能である。

以上の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。

（付記１）
衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１確率分布取得手段と、
前記第１の期間内の単位時刻である各時刻ｔについて前記第１確率分布取得手段により取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する不確実性指標算出手段と、
前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記不確実性指標算出手段により算出された前記不確実性指標の平均値を算出し、その平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第２確率分布取得手段と、
前記第２確率分布取得手段により各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する変化点指標算出手段と、
前記変化点指標算出手段により算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する異常値検出手段と
を有することを特徴とする衛星測位システムの異常値検出装置。

（付記２）
前記第２確率分布取得手段は、前記第２の期間内で前記不確実性指標の平均値を算出し、その平均値を前記第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間の長さ以下の期間である第３の期間における前記第２の確率密度関数を算出して前記第２の確率分布として取得し、
前記変化点指標算出手段は、前記時刻ｔ−１に対応する前記第２の確率分布に基づき、前記時刻ｔ−１以前で、前記第２の期間以下の長さである第４の期間内の前記第２の確率分布の平均情報量を前記変化点指標として算出する
ことを特徴とする付記１記載の衛星測位システムの異常値検出装置。

（付記３）
前記時系列値は、衛星測位システムの各人工衛星と受信機との間の擬似距離や測位信号の位相などの、前記受信機において受信した前記人工衛星からの測位信号に基づいて時系列に取得されたデータであることを特徴とする付記１又は２記載の衛星測位システムの異常値検出装置。

（付記４）
前記第１確率分布取得手段は、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け平均値をｗ１_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、新しい時系列値ほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け分散値をσ１_t-1としたとき、次式

により算出される前記第１の確率密度関数ｐ(ｘ_t)を前記第１の確率分布として取得することを特徴とする付記１乃至３のうちいずれか一項記載の衛星測位システムの異常値検出装置。

（付記５）
前記不確実性指標算出手段は、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記時刻ｔ−１における前記第１の確率密度関数をｐ_t-1(ｘ_t|ｘ^t-1)としたとき、次式

により算出されるＵＬＩ(ｔ)を前記時刻ｔにおける前記不確実性指標として算出することを特徴とする付記１乃至４のうちいずれか一項記載の衛星測位システムの異常値検出装置。

（付記６）
前記第２確率分布取得手段は、前記第２の期間をＴ２、時刻ｉに対応する前記不確実性指標をＵＬＩ(ｉ)としたとき、次式

により、前記時刻ｔ以前の前記第２の期間における前記不確実性指標の平均値である前記第２の確率変数ｙ_tを算出し、更に前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け平均値をｗ２_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け分散値をσ２_t-1としたとき、次式

により算出される前記第２の確率密度関数ｑ(ｙ_t)を前記第２の確率分布として取得することを特徴とする付記１乃至５のうちいずれか一項記載の衛星測位システムの異常値検出装置。

（付記７）
前記変化点指標算出手段は、前記時刻ｔ−１以前の前記第４の期間をＴ４、前記時刻ｉにおける前記第２の確率密度関数をｑ_i-1(ｙ_i|ｙ^i-1)としたとき、次式

により算出されるＣＰＩ(t)を前記時刻ｔにおける前記変化点指標として算出することを特徴とする付記２又は６記載の衛星測位システムの異常値検出装置。

（付記８）
前記第２確率分布取得手段により各時刻ｔ毎に取得された前記第２の確率分布を確率変数として、その確率変数の重み付け統計量に基づき、第３の確率密度関数を算出する算出手段を更に有し、
前記変化点指標算出手段は、前記第３の確率密度関数から変曲点指標を算出し、前記異常値検出手段は、前記変曲点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変曲点指標に対応する前記時系列値を異常値として検出することを特徴とする付記１記載の衛星測位システムの異常値検出装置。

（付記９）
衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、
前記第１の期間内の単位時刻である各時刻ｔについて前記第１のステップにより取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、
前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記第２のステップにより算出された前記不確実性指標の平均値を算出する第３のステップと、
前記第３のステップで算出された前記平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、
前記第４のステップにより各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、
前記第５のステップにより算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する第６のステップと
を含むことを特徴とする衛星測位システムの異常値検出方法。

（付記１０）
前記第４のステップは、前記平均値を前記第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間の長さ以下の期間である第３の期間における前記第２の確率密度関数を算出して前記第２の確率分布として取得し、
前記第５のステップは、前記時刻ｔ−１に対応する前記第２の確率分布に基づき、前記時刻ｔ−１以前で、前記第２の期間以下の長さである第４の期間内の前記第２の確率分布の平均情報量を前記変化点指標として算出することを特徴とする付記９記載の衛星測位システムの異常値検出方法。

（付記１１）
前記時系列値は、衛星測位システムの各人工衛星と受信機との間の擬似距離や測位信号の位相などの、前記受信機において受信した前記人工衛星からの測位信号に基づいて時系列に取得されたデータであることを特徴とする付記９又は１０記載の衛星測位システムの異常値検出方法。

（付記１２）
前記第１のステップは、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け平均値をｗ１_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、新しい時系列値ほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け分散値をσ１_t-1としたとき、次式

により算出される前記第１の確率密度関数ｐ(ｘ_t)を前記第１の確率分布として取得することを特徴とする付記９乃至１１のうちいずれか一項記載の衛星測位システムの異常値検出方法。

（付記１３）
前記第２のステップは、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記時刻ｔ−１における前記第１の確率密度関数をｐ_t-1(ｘ_t|ｘ^t-1)としたとき、次式

により算出されるＵＬＩ(ｔ)を前記時刻ｔにおける前記不確実性指標として算出することを特徴とする付記９乃至１２のうちいずれか一項記載の衛星測位システムの異常値検出方法。

（付記１４）
前記第３のステップは、前記第２の期間をＴ２、時刻ｉに対応する前記不確実性指標をＵＬＩ(ｉ)としたとき、次式

により、前記時刻ｔ以前の前記第２の期間における前記不確実性指標の平均値である前記第２の確率変数ｙ_tを算出することを特徴とする付記９乃至１３のうちいずれか一項記載の衛星測位システムの異常値検出方法。

（付記１５）
前記第４のステップは、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け平均値をｗ２_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け分散値をσ２_t-1としたとき、次式

により算出される前記第２の確率密度関数ｑ(ｙ_t)を前記第２の確率分布として取得することを特徴とする付記１４記載の衛星測位システムの異常値検出方法。

（付記１６）
前記第５のステップは、前記時刻ｔ−１以前の前記第４の期間をＴ４、前記時刻ｉにおける前記第２の確率密度関数をｑ_i-1(ｙ_i|ｙ^i-1)としたとき、次式

により算出されるＣＰＩ(t)を前記時刻ｔにおける前記変化点指標として算出することを特徴とする付記１０又は１５記載の衛星測位システムの異常値検出方法。

（付記１７）
コンピュータに、
衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、
前記第１の期間内の単位時刻である各時刻ｔについて前記第１のステップにより取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、
前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記第２のステップにより算出された前記不確実性指標の平均値を算出する第３のステップと、
前記第３のステップで算出された前記平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、
前記第４のステップにより各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、
前記第５のステップにより算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する第６のステップと
を実行させることを特徴とする衛星測位システムの異常値検出プログラム。

（付記１８）
前記第４のステップは、前記平均値を前記第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間の長さ以下の期間である第３の期間における前記第２の確率密度関数を算出して前記第２の確率分布として取得し、
前記第５のステップは、前記時刻ｔ−１に対応する前記第２の確率分布に基づき、前記時刻ｔ−１以前で、前記第２の期間以下の長さである第４の期間内の前記第２の確率分布の平均情報量を前記変化点指標として算出することを特徴とする付記１７記載の衛星測位システムの異常値検出プログラム。

（付記１９）
前記時系列値は、衛星測位システムの各人工衛星と受信機との間の擬似距離や測位信号の位相などの、前記受信機において受信した前記人工衛星からの測位信号に基づいて時系列に取得されたデータであることを特徴とする付記１７又は１８記載の衛星測位システムの異常値検出プログラム。

（付記２０）
前記第１のステップは、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け平均値をｗ１_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、新しい時系列値ほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け分散値をσ１_t-1としたとき、次式

により算出される前記第１の確率密度関数ｐ(ｘ_t)を前記第１の確率分布として取得することを特徴とする付記１７乃至１９のうちいずれか一項記載の衛星測位システムの異常値検出プログラム。

（付記２１）
前記第２のステップは、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記時刻ｔ−１における前記第１の確率密度関数をｐ_t-1(ｘ_t|ｘ^t-1)としたとき、次式

により算出されるＵＬＩ(ｔ)を前記時刻ｔにおける前記不確実性指標として算出することを特徴とする付記１７乃至２０のうちいずれか一項記載の衛星測位システムの異常値検出プログラム。

（付記２２）
前記第３のステップは、前記第２の期間をＴ２、時刻ｉに対応する前記不確実性指標をＵＬＩ(ｉ)としたとき、次式

により、前記時刻ｔ以前の前記第２の期間における前記不確実性指標の平均値である前記第２の確率変数ｙ_tを算出することを特徴とする付記１７乃至２１のうちいずれか一項記載の衛星測位システムの異常値検出プログラム。

（付記２３）
前記第４のステップは、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け平均値をｗ２_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け分散値をσ２_t-1としたとき、次式

により算出される前記第２の確率密度関数ｑ(ｙ_t)を前記第２の確率分布として取得することを特徴とする付記２２記載の衛星測位システムの異常値検出プログラム。

（付記２４）
前記第５のステップは、前記時刻ｔ−１以前の前記第４の期間をＴ４、前記時刻ｉにおける前記第２の確率密度関数をｑ_i-1(ｙ_i|ｙ^i-1)としたとき、次式

により算出されるＣＰＩ(t)を前記時刻ｔにおける前記変化点指標として算出することを特徴とする付記１８又は２３記載の衛星測位システムの異常値検出プログラム。

１０ 異常値検出装置
１１ 異常値指標算出部
１２ 動的モデル構成部
１３ 異常値検出部
１２１ 第１確率分布取得部
１２２ 不確実性指標算出部
１２３ 第２確率分布取得部
１２４ 変化点指標算出部

Claims (10)

1. 衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１確率分布取得手段と、
   前記第１の期間内の単位時刻である各時刻ｔについて前記第１確率分布取得手段により取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する不確実性指標算出手段と、
   前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記不確実性指標算出手段により算出された前記不確実性指標の平均値を算出し、その平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第２確率分布取得手段と、
   前記第２確率分布取得手段により各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する変化点指標算出手段と、
   前記変化点指標算出手段により算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する異常値検出手段と
   を有することを特徴とする衛星測位システムの異常値検出装置。
2. 前記第２確率分布取得手段は、前記第２の期間内で前記不確実性指標の平均値を算出し、その平均値を前記第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間の長さ以下の期間である第３の期間における前記第２の確率密度関数を算出して前記第２の確率分布として取得し、
   前記変化点指標算出手段は、前記時刻ｔ−１に対応する前記第２の確率分布に基づき、前記時刻ｔ−１以前で、前記第２の期間以下の長さである第４の期間内の前記第２の確率分布の平均情報量を前記変化点指標として算出する
   ことを特徴とする請求項１記載の衛星測位システムの異常値検出装置。
3. 前記時系列値は、衛星測位システムの各人工衛星と受信機との間の擬似距離や測位信号の位相などの、前記受信機において受信した前記人工衛星からの測位信号に基づいて時系列に取得されたデータであることを特徴とする請求項１又は２記載の衛星測位システムの異常値検出装置。
4. 前記第１確率分布取得手段は、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け平均値をｗ１_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、新しい時系列値ほど大きな値になるように重み付けを行った前記第１の確率変数ｘ_tの重み付け分散値をσ１_t-1としたとき、次式
   

   

   により算出される前記第１の確率密度関数ｐ(ｘ_t)を前記第１の確率分布として取得することを特徴とする請求項１乃至３のうちいずれか一項記載の衛星測位システムの異常値検出装置。
5. 前記不確実性指標算出手段は、前記時刻ｔにおける前記第１の確率変数をｘ_t、前記時刻ｔ−１における前記第１の確率密度関数をｐ_t-1(ｘ_t|ｘ^t-1)としたとき、次式
   

   

   により算出されるＵＬＩ(ｔ)を前記時刻ｔにおける前記不確実性指標として算出することを特徴とする請求項１乃至４のうちいずれか一項記載の衛星測位システムの異常値検出装置。
6. 前記第２確率分布取得手段は、前記第２の期間をＴ２、時刻ｉに対応する前記不確実性指標をＵＬＩ(ｉ)としたとき、次式
   

   

   により、前記時刻ｔ以前の前記第２の期間における前記不確実性指標の平均値である前記第２の確率変数ｙ_tを算出し、更に前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け平均値をｗ２_t-1、前記第１の期間のうち、前記時刻ｔ−１以前の期間における、前記時系列値が新しいほど大きな値になるように重み付けを行った前記第２の確率変数ｙ_tの重み付け分散値をσ２_t-1としたとき、次式
   

   

   により算出される前記第２の確率密度関数ｑ(ｙ_t)を前記第２の確率分布として取得することを特徴とする請求項１乃至５のうちいずれか一項記載の衛星測位システムの異常値検出装置。
7. 前記変化点指標算出手段は、前記時刻ｔ−１以前の前記第４の期間をＴ４、前記時刻ｉにおける前記第２の確率密度関数をｑ_i-1(ｙ_i|ｙ^i-1)としたとき、次式
   

   

   により算出されるＣＰＩ(t)を前記時刻ｔにおける前記変化点指標として算出することを特徴とする請求項２記載の衛星測位システムの異常値検出装置。
8. 衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、
   前記第１の期間内の単位時刻である各時刻ｔについて前記第１のステップにより取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、
   前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記第２のステップにより算出された前記不確実性指標の平均値を算出する第３のステップと、
   前記第３のステップで算出された前記平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、
   前記第４のステップにより各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、
   前記第５のステップにより算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する第６のステップと
   を含むことを特徴とする衛星測位システムの異常値検出方法。
9. 前記第４のステップは、前記平均値を前記第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間の長さ以下の期間である第３の期間における前記第２の確率密度関数を算出して前記第２の確率分布として取得し、
   前記第５のステップは、前記時刻ｔ−１に対応する前記第２の確率分布に基づき、前記時刻ｔ−１以前で、前記第２の期間以下の長さである第４の期間内の前記第２の確率分布の平均情報量を前記変化点指標として算出することを特徴とする請求項８記載の衛星測位システムの異常値検出方法。
10. コンピュータに、
    衛星測位システムの異常値の検出対象である時系列値を第１の確率変数として、前記時系列値の重み付け統計量に基づき、所定の時刻以前の第１の期間における第１の確率密度関数を算出して第１の確率分布として取得する第１のステップと、
    前記第１の期間内の単位時刻である各時刻ｔについて前記第１のステップにより取得された前記第１の確率分布に基づき、前記時刻ｔの直前の時刻ｔ−１における前記第１の確率変数の選択情報量を不確実性指標として算出する第２のステップと、
    前記時刻ｔ以前の、前記第１の期間より短い第２の期間内で、前記第２のステップにより算出された前記不確実性指標の平均値を算出する第３のステップと、
    前記第３のステップで算出された前記平均値を第２の確率変数として、前記第２の確率変数の重み付け統計量に基づき、前記第１の期間における第２の確率密度関数を算出して第２の確率分布として取得する第４のステップと、
    前記第４のステップにより各時刻ｔ毎に取得された前記第２の確率分布に基づき、前記第２の確率分布の取得時刻ｔの直前の時刻ｔ−１以前の前記第２の期間内の前記第２の確率分布の平均情報量を変化点指標として算出する第５のステップと、
    前記第５のステップにより算出された前記変化点指標を予め設定された閾値と比較し、前記閾値より大きな値の前記変化点指標に対応する前記時系列値を前記異常値として検出する第６のステップと
    を実行させることを特徴とする衛星測位システムの異常値検出プログラム。

Images