(57)【要約】      （修正有）【課題】飛翔体の位置に関する計算の結果が収束するまでに要する時間を短縮する。【解決手段】本位置推定方法は、飛翔体の観測情報を取得し、共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを観測情報に適用して、飛翔体の推定位置を算出する処理を含む。【選択図】図５
2018151300.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、
飛翔体の観測情報を取得し、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記観測情報に適用して、前記飛翔体の推定位置を算出する、
処理を実行させる位置推定プログラム。
【請求項２】
前記第１の所定値は、全特異点の留数の総和に基づく値であり、
前記第２の所定値は、前記観測情報において発生しうる最大の誤差に基づく値である、
請求項１記載の位置推定プログラム。
【請求項３】
前記非線形のベイズフィルタは、拡張カルマンフィルタである、
請求項１又は２記載の位置推定プログラム。
【請求項４】
前記推定位置を算出する処理において、
前記観測情報または前回算出された推定値に基づき、非慣性系での運動方程式に従って前記飛翔体の予測位置を算出し、
前記予測位置と前記観測情報とに基づき、カルマンゲインを算出し、
前記予測位置と前記カルマンゲインとに基づき、前記推定位置を算出する、
請求項３記載の位置推定プログラム。
【請求項５】
前記予測位置を算出する処理において、
前記予測位置の座標系を回転座標系から局所座標系に変換し、
前記推定位置を算出する処理において、
前記局所座標系での前記予測位置と前記カルマンゲインとに基づき、前記局所座標系での前記推定位置を算出する、
請求項４記載の位置推定プログラム。
【請求項６】
前記観測情報は、赤外線検出装置により取得された、前記飛翔体を含む二次元画像の情報を含み、
前記推定位置を算出する処理において、
少なくとも前記二次元画像の情報に基づき、前記飛翔体の方位を算出し、
前記方位に基づき、前記飛翔体の観測位置を算出する、
請求項４記載の位置推定プログラム。
【請求項７】
前記推定位置を算出する処理において、
前記観測情報から前記飛翔体の推定速度をさらに算出し、
前記コンピュータに、
前記推定位置および前記推定速度から、前記飛翔体の落下位置の範囲を算出し、
前記落下位置の範囲の情報を出力する、
処理をさらに実行させる請求項１記載の位置推定プログラム。
【請求項８】
前記システム誤差の共分散は、三次元位置の各成分および三次元速度の各成分についての要素を含む６行６列の行列で表される、
請求項２記載の位置推定プログラム。
【請求項９】
前記状態方程式は、地球の中心を原点とした直交座標系での位置および速度についての項を含み、
前記観測方程式は、局所座標系での角度についての項を含む、
請求項１記載の位置推定プログラム。
【請求項１０】
コンピュータが、
飛翔体の観測情報を取得し、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記観測情報に適用して、前記飛翔体の推定位置を算出する、
処理を実行する位置推定方法。
【請求項１１】
飛翔体の観測情報を取得する取得部と、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記取得部により取得された前記観測情報に適用して、前記飛翔体の推定位置を算出する算出部と、
を有する位置推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、飛翔体の位置の推定技術に関する。
【背景技術】
【０００２】
飛翔体（弾道弾、ロケット、或いは人工衛星等）の位置を推定する技術が存在する。例えば、或る文献は、飛翔している目標の観測情報を用いて、最小二乗法によるバッチフィルタに基づき目標位置を推定する技術を開示する。
【０００３】
上記のような技術においては、計算結果が収束した場合あるいは計算結果が収束したとみなすことができる時間が経過した場合に、飛翔体の推定位置および推定速度等を用いて飛翔体の落下位置を予測する計算が開始される。従って、計算結果が収束するまでの時間を短縮することができれば、より早い段階から飛翔体の落下の予測位置を知ることができる。これにより、飛翔体の落下に備えた対策をより早い段階で講じることができる。
【０００４】
しかし、最小二乗法に基づく計算を実行する場合、観測値のサンプリングの間隔を短くすることでサンプル数を増加させると、計算量が増大する。結果として、計算結果が収束するまでの時間が長くなる。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００２−２０２３６７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明の目的は、１つの側面では、飛翔体の位置に関する計算の結果が収束するまでに要する時間を短縮するための技術を提供することである。
【課題を解決するための手段】
【０００７】
一態様に係る位置推定方法は、飛翔体の観測情報を取得し、共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを観測情報に適用して、飛翔体の推定位置を算出する処理を含む。
【発明の効果】
【０００８】
１つの側面では、飛翔体の位置に関する計算の結果が収束するまでに要する時間を短縮できるようになる。
【図面の簡単な説明】
【０００９】
【図１】図１は、飛翔体である弾道弾の位置の算出について説明するための図である。
【図２】図２は、運動方程式に基づき弾道弾の位置を予測する処理について説明するための図である。
【図３】図３は、標定処理に基づく着弾位置の予測について説明するための図である。
【図４】図４は、本実施の形態に係るシステムの構成図である。
【図５】図５は、情報処理装置の機能ブロック図である。
【図６】図６は、センサ信号処理部が実行する処理の処理フローを示す図である。
【図７】図７は、標定処理部および着弾予測部が実行する処理の処理フローを示す図である。
【図８】図８は、局所座標系での原点を示す図である。
【図９】図９は、局所座標系での原点と弾道弾の位置との関係を示す図である。
【図１０】図１０は、ＥＣＥＦ座標系とＮＥＤ座標系との関係を示す図である。
【図１１】図１１は、ＮＥＤ座標系と機体座標系と視軸座標系との関係を示す図である。
【図１２】図１２は、標定処理の処理フローを示す図である。
【図１３】図１３は、弾道弾の位置と観測機の位置との関係を示す図である。
【図１４】図１４は、着弾予測処理の処理フローを示す図である。
【図１５】図１５は、予測着弾範囲の算出について説明するための図である。
【図１６】図１６は、シミュレーションにおける弾道弾を示す図である。
【図１７】図１７は、飛翔シナリオを示す図である。
【図１８】図１８は、シミュレーション条件を示す図である。
【図１９】図１９は、シミュレーション条件を示す図である。
【図２０】図２０は、シミュレーションの処理フローを示す図である。
【図２１】図２１は、計算結果が収束するまでの時間を示す図である。
【図２２】図２２は、サンプリング間隔毎に標定誤差を示す図である。
【図２３】図２３は、標定誤差と観測時間との関係を示す図である。
【図２４】図２４は、システム誤差の共分散の設定値を変化させた場合における標定誤差を示す図である。
【図２５】図２５は、第２の設定値を使用した場合における標定誤差を示す図である。
【図２６】図２６は、着弾位置の予測誤差を示す図である。
【図２７】図２７は、着弾位置の予測誤差と観測時間との関係を示す図である。
【発明を実施するための形態】
【００１０】
図１は、飛翔体である弾道弾３の位置の算出について説明するための図である。本実施の形態の主要な処理を実行する観測機１は、弾道弾３の位置情報を二次元カメラ座標（Ｈ，Ｖ）（すなわち画像上の画素位置）として取得する。観測機１のＥＣＥＦ（Earth Center Earth Fixed）座標系での位置（ａ，ｂ，ｃ）、機体座標系での視軸方位（θe，Φe，εe）およびＮＥＤ（North-East-Down）座標系での姿勢角（α，β，γ）は、観測機１において計測可能である。
【００１１】
上記の二次元カメラ座標、姿勢角および視軸方位を用いて、弾道弾３の局所座標系での目標方位（θ，Φ，ε）を算出することはできる。しかし、観測機１から弾道弾３までの距離が判明しなければ、弾道弾３の位置（ｘ，ｙ，ｚ）を算出することはできない。
【００１２】
そこで本実施の形態においては、観測機１から弾道弾３までの距離が所定距離であることを仮定して弾道弾３の位置（ｘ，ｙ，ｚ）が計算される。
【００１３】
図２は、運動方程式に基づき弾道弾３の位置を予測する処理について説明するための図である。本実施の形態のように目標が弾道弾３である場合、発射台から発射された後にブーストフェーズが終了すると弾道弾３は自由落下するので、弾道弾３の移動を１つの運動方程式で表すことができる。具体的には、観測機１から弾道弾３までの距離が所定距離であることを仮定して計算された位置（ｘt1，ｙt1，ｚt1）を起点として運動方程式に基づき計算を行えば、単位時間後の位置（ｘt2，ｙt2，ｚt2）を予測することができる。
【００１４】
図３は、標定処理に基づく着弾位置の予測について説明するための図である。本実施の形態においては、弾道弾３のブーストフェーズが終了すると標定処理が開始される。標定処理においては、サンプリング間隔ごとに計算された目標方位とその予測値との誤差が小さくなるように弾道弾３の推定位置および推定速度が修正される。本実施の形態においては、計算結果が収束したとみなすことができる時間が経過した場合（例えば数十秒）、最終的に算出された推定位置および推定速度を用いて、運動方程式に従って予測の着弾位置（より具体的には、着弾の予測範囲）が計算される。
【００１５】
図４は、本実施の形態に係るシステムの構成図である。本実施の形態における観測機１の機外には、ＧＰＳ（Global Positioning System）アンテナ１０００と、捜索センサ１００１と、追尾センサ１００２と、ジンバル１００３と、位置方位基準器１００４とが設置される。観測機１の機内には、情報処理装置１０が設置される。
【００１６】
ＧＰＳアンテナ１０００により、観測機１のＥＣＥＦ座標系での位置の情報が取得される。捜索センサ１００１および追尾センサ１００２は、パッシブ二次元撮像装置である赤外線センサである。捜索センサ１００１は、弾道弾３を探知し、弾道弾３の二次元カメラ座標を取得する。追尾センサ１００２は、捜索センサ１００１により検知された弾道弾３の飛翔に合わせて弾道弾３の二次元カメラ座標を取得する。ジンバル１００３は、捜索センサ１００１および追尾センサ１００２の向きを制御する。位置方位基準器１００４は、機体情報（例えば、姿勢角および視軸方位等の情報）を取得する。
【００１７】
情報処理装置１０は、１又は複数のＣＰＵ（Central Processing Unit）１５１と、１又は複数のメモリ１５３（例えばＤＲＡＭ（Dynamic Random Access Memory））と、１又は複数のＨＤＤ（Hard Disk Drive）１５５とを有する。
【００１８】
情報処理装置１０は、取得されたカメラ座標、位置情報、姿勢角および視軸方位等に基づき処理を実行する。また、情報処理装置１０は、処理結果に基づきジンバル１００３を制御する。
【００１９】
このように、本実施の形態に係るシステムはＩＲＳＴ（Infra-Red Search and Track）システムの機能を有する。
【００２０】
図５は、情報処理装置１０の機能ブロック図である。情報処理装置１０は、ジンバル制御部１０１と、センサ信号処理部１０３と、標定処理部１０５と、着弾予測部１０７と、標定結果格納部１１１と、出力データ格納部１１３とを有する。
【００２１】
ジンバル制御部１０１、センサ信号処理部１０３、標定処理部１０５および着弾予測部１０７は、例えば、図４に示したメモリ１５３にロードされたプログラムがＣＰＵ１５１に実行されることで実現される。標定結果格納部１１１および出力データ格納部１１３は、例えば、メモリ１５３又はＨＤＤ１５５に設けられる。
【００２２】
ジンバル制御部１０１は、後述の標定処理の結果に基づきジンバル１００３の動作を制御する。センサ信号処理部１０３は、弾道弾３の方位（以下、目標方位と呼ぶ）を算出し、算出した目標方位の情報を標定処理部１０５に出力する。標定処理部１０５は、センサ信号処理部１０３の処理の結果およびその他の情報（例えば設定値等）に基づき、後述の標定処理を実行し、標定処理の結果を標定結果格納部１１１に格納する。着弾予測部１０７は、標定結果格納部１１１に格納されているデータに基づき処理を実行し、処理結果を出力データ格納部１１３に格納する。
【００２３】
次に、図６乃至図１５を用いて、情報処理装置１０の動作について説明する。
【００２４】
図６は、センサ信号処理部１０３が実行する処理の処理フローを示す図である。本処理は、サンプリング間隔（例えば５００ミリ秒）ごとに実行される。
【００２５】
センサ信号処理部１０３は、センサデータ（具体的には、弾道弾３の二次元カメラ座標）を捜索センサ１００１および追尾センサ１００２から取得する。また、センサ信号処理部１０３は、機体情報（例えば、姿勢角および視軸方位等の情報）をジンバル制御部１０１を介して位置方位基準器１００４から取得する。さらに、センサ信号処理部１０３は、観測機１の位置情報をＧＰＳアンテナ１０００から取得する（図６：ステップＳ１）。
【００２６】
センサ信号処理部１０３は、ステップＳ１において取得した情報に基づき弾道弾３の目標方位を算出し（ステップＳ３）、メモリ１５３等の記憶装置に格納する。具体的には、アジマス角ＡＺL（ε＝０，Φ）及びエレベーション角ＥＬL（ε＝０，θ）が算出される。なお、座標系は局所座標系であり、視軸方位と同じ座標系であってもよいし異なる座標系であってもよい。そして処理は終了する。
【００２７】
以上のような処理により算出された弾道弾３の目標方位を用いれば、後の処理において弾道弾３の位置を算出できるようになる。
【００２８】
図７は、標定処理部１０５および着弾予測部１０７が実行する処理の処理フローを示す図である。本処理は、捜索センサ１００１により弾道弾３が探知され且つ追尾センサ１００２による追尾が開始された後に実行される。
【００２９】
まず、本実施の形態のシステムモデルを示す。本実施の形態のシステムモデルは拡張カルマンフィルタに基づいており、状態方程式は以下のように表される。
【００３０】
【数１】
2018151300.tif 000003
【００３１】
ここで、ｖ（ｔ）はシステム誤差であり、共分散行列がＱ（ｔ）であり且つ０を平均とする分布に従う。
【００３２】
【数２】
2018151300.tif 000004
【００３３】
本実施の形態における観測方程式は、以下の観測方程式であるとする。
【００３４】
【数３】
2018151300.tif 000005
【００３５】
ここで、ｗ（ｔ）は観測残差であり、共分散行列がＲ（ｔ）であり且つ０を平均とする分布に従う。
【００３６】
【数４】
2018151300.tif 000006
【００３７】
なお、以下では、共分散行列のことを単に共分散と呼ぶ場合が有る。
【００３８】
時間遷移モデルＦ（ｔ）は、６行６列の行列として以下のように表される。
【００３９】
【数５】
2018151300.tif 000007
【００４０】
なお、ハットの記号は予測値を表し、バーの記号は推定値を表す。
【００４１】
状態変数ｘ（ｔ）は、以下のように表される。
【００４２】
【数６】
2018151300.tif 000008
【００４３】
観測モデルＨ（ｔ）は、２行６列の行列として以下のように表される。
【００４４】
【数７】
2018151300.tif 000009
【００４５】
そして、標定処理部１０５は、弾道弾３のブーストフェーズが終了したか判定する（図７：ステップＳ１１）。ブーストフェーズが終了したか否かは、例えば、二次元カメラ座標の遷移から推定される加速度等に基づいて判定される。
【００４６】
弾道弾３のブーストフェーズが終了していない場合（ステップＳ１１：Ｎｏルート）、処理はステップＳ１１に戻る。
【００４７】
一方、弾道弾３のブーストフェーズが終了した場合（ステップＳ１１：Ｙｅｓルート）、標定処理部１０５は、弾道弾３の局所座標系での初期位置を算出する（ステップＳ１３）。ステップＳ１３においては、例えば以下の式に従って初期位置が算出される。
【００４８】
【数８】
2018151300.tif 000010
【００４９】
Ｒは、観測機１から弾道弾３が初めて探知された位置までの距離の設定値を表す。本実施の形態においては、例えば、Ｒ＝１２００（ｋｍ）としてＲが設定される。
【００５０】
図８は、局所座標系での原点を示す図である。ここでは、局所座標系は三次元直交座標系であり、捜索センサ１００１および追尾センサ１００２の位置（すなわち観測機１の位置）が原点Ｏである。Ｘ軸の正方向は北方向であり、Ｙ軸の正方向は東方向であり、Ｚ軸の正方向は天頂方向である。
【００５１】
図９は、局所座標系での原点と弾道弾３の位置との関係を示す図である。Ｏは原点を表し、観測機１の位置に相当する。原点からの距離がＲである位置９０１が弾道弾３の位置に相当する。弾道弾３の位置９０１は、原点を基準として、アジマス角Φおよびエレベーション角θにより特定される。
【００５２】
なお、標定処理部１０５は、弾道弾３の初期位置の座標を、最終的にＥＣＥＦ座標系での座標に変換する。
【００５３】
ここで、座標系の変換について説明する。
【００５４】
図１０は、ＥＣＥＦ座標系とＮＥＤ座標系との関係を示す図である。ＥＣＥＦ座標系は直交座標系であり、原点は地球の中心に有る。一方、ＮＥＤ座標系はセンサ（本実施の形態においては観測機１）の位置（図１０においては、経度がλであり且つ緯度がδである位置）を原点とする直交座標系である。上で述べたように、Ｘ軸の正方向は北方向であり、Ｙ軸の正方向は東方向であり、Ｚ軸の正方向は天頂方向である。なお、ＲECEF_NEDが、ＥＣＥＦ座標系を観測機１の位置を原点とする局所座標系の座標に変換する行列式であるとすると、ＲNED_ECEF＝ＲECEF_NED-1が成立する。
【００５５】
図１１は、ＮＥＤ座標系と機体座標系と視軸座標系との関係を示す図である。図１１において、ＲNED_ACは姿勢角（ロール，ピッチ，ヨー）を用いてＮＥＤ座標系を機体座標系に変換する行列式であり、ＲAC_VIは視軸方位を用いて機体座標系を視軸座標系に変換する行列式である。この場合、ＲAC_NED＝ＲNED_AC-1が成立し、また、ＲVI_AC＝ＲAC_VI-1が成立する。なお、本実施の形態においては、視軸座標系でローテーションθが０である場合、ε＝ＡＺLとし、Φ＝ＥＬLとしている。
【００５６】
また、ＥＣＥＦ座標系と測地座標系との関係は以下のとおりである。
【００５７】
【数９】
2018151300.tif 000011
【００５８】
【数１０】
2018151300.tif 000012
【００５９】
【数１１】
2018151300.tif 000013
【００６０】
ここで、Ｘ、ＹおよびＺはそれぞれＥＣＥＦ座標系の座標を表す。λ、Φ及びｈはそれぞれ緯度、経度および高度を表す。ａは赤道半径を表し、６，３７８，１３７（ｍ）である。ｆは扁平率を表し、１／２９８．２５７２２３５６３である。
【００６１】
ＥＣＩ（Earth-Centered Inertial）座標系およびＥＣＥＦ座標系には、地球自転によるＺ軸まわりの回転角についてΔΦの違いがあり、ＥＣＩ座標系とＥＣＥＦ座標系との関係は以下のとおりである。
【００６２】
【数１２】
2018151300.tif 000014
【００６３】
ここで、ｘ、ｙおよびｚはそれぞれＥＣＩ座標系の座標を表す。ΔΦはＥＣＩ座標系に対する、ＥＣＥＦ座標系のｚ軸（Ｚ軸）軸周りの回転角を表す。
【００６４】
測地座標系からＥＣＩ座標系またはＥＣＥＦ座標系への変換は、上で示した式への代入により行うことができる。
【００６５】
ＥＣＩ座標系またはＥＣＥＦ座標系から測地座標系への座標変換は、以下の式のとおりである。
【００６６】
【数１３】
2018151300.tif 000015
【００６７】
ここで、ａは地球長半径を表し、６３７８１３７．０（ｍ）である。ｅは離心率を表し、０．０８１８１９１９１０４２８１５８である。ｆは扁平率を表し、０．００３３５２８１１である。
【００６８】
図７の説明に戻り、標定処理部１０５は、弾道弾３の初期速度を算出する（ステップＳ１５）。弾道弾３の初期速度は、例えば、弾道弾３の初期位置と運動モデルとから算出される。但し、その他の方法で算出してもよいし、予め初期速度を設定してもよい。
【００６９】
標定処理部１０５は、システム誤差の初期共分散および観測残差の初期共分散を設定する（ステップＳ１７）。システム誤差の初期共分散は、例えば、微小時間において弾道弾３が等加速に移動するという仮定した場合における全特異点の留数の総和に基づく。観測残差の初期共分散は、例えば、観測時に発生しうる最大誤差に基づく。但し、微小時間において弾道弾３が等加速に移動するという仮定してシステム誤差の初期共分散および観測残差の初期共分散を設定してもよい。
【００７０】
標定処理部１０５は、弾道弾３の高度が０より小さいか判定する（ステップＳ１９）。
【００７１】
弾道弾３の高度が０より小さい場合（ステップＳ１９：Ｙｅｓルート）、処理は終了する。弾道弾３の高度が０より小さくない場合（ステップＳ１９：Ｎｏルート）、標定処理部１０５は、標定処理を実行する（ステップＳ２１）。標定処理については、図１２及び図１３を用いて説明する。
【００７２】
まず、標定処理部１０５は、座標系をＥＣＥＦ座標系から回転座標系（すなわち非慣性座標系）に変換する（図１２：ステップＳ３１）。
【００７３】
標定処理部１０５は、弾道弾３の予測位置および予測速度を算出する（ステップＳ３３）。ステップＳ３３においては、以下の式に基づき予測位置が算出される。
【００７４】
【数１４】
2018151300.tif 000016
【００７５】
ここで、ｘハット（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）におけるｘ座標の予測値であり、ｙハット（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）におけるｙ座標の予測値であり、ｚハット（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）におけるｚ座標の予測値である。ｘバー（ｔ）は時刻ｔにおけるｘ座標の推定値であり、ｙバー（ｔ）は時刻ｔにおけるｙ座標の推定値であり、ｚバー（ｔ）は時刻ｔにおけるｚ座標の推定値であるが、ｔ＝０の場合には弾道弾３の初期位置である。Ｖｘバー（ｔ）は時刻ｔにおけるｘ方向の速度の推定値であり、Ｖｙバー（ｔ）は時刻ｔにおけるｙ方向の速度の推定値であり、Ｖｚバー（ｔ）は時刻ｔにおけるｚ方向の速度の推定値であるが、ｔ＝０の場合には弾道弾３の初期速度である。Δｔはサンプリング間隔（秒）である。
【００７６】
なお、弾道弾３の予測速度は、例えば、時刻ｔにおける弾道弾３の速度の推定値と、重力加速度（ｍ／ｓ2）と、コリオリ力と、遠心力とに基づき算出される。
【００７７】
標定処理部１０５は、システム誤差の共分散の予測値を算出する（ステップＳ３５）。システム誤差の共分散の予測値は、以下のように算出される。
【００７８】
【数１５】
2018151300.tif 000017
【００７９】
ここで、Ｐハット（ｔ＋Δｔ｜ｔ）は、時刻ｔにおけるシステム誤差の共分散の予測値を表す。Ｐバー（ｔ｜ｔ）は、時刻ｔにおけるシステム誤差の共分散の推定値を表す。
【００８０】
標定処理部１０５は、座標系を回転座標系から局所座標系（すなわち慣性座標系）に変換する（ステップＳ３７）。ステップＳ３７においては、時刻（ｔ＋Δｔ）における観測機１の位置に基づいて変換が行われる。
【００８１】
標定処理部１０５は、時刻（ｔ＋Δｔ）における位置を基準として目標方位の予測値を算出する（ステップＳ３９）。目標方位の予測値は、以下のように算出される。
【００８２】
【数１６】
2018151300.tif 000018
【００８３】
ここで、ＥＬL（ｔ｜ｔ）およびＡＺL（ｔ｜ｔ）は、第２象限および第４象限の場合、以下のとおりである。
【００８４】
【数１７】
2018151300.tif 000019
【００８５】
図１３は、弾道弾３の位置と観測機１の位置との関係を示す図である。Ｏは原点を表し、観測機１の位置に相当する。弾道弾３の位置１３０１は、原点を基準として、アジマス角Φおよびエレベーション角θにより特定される。Ｖｘ、ＶｙおよびＶｚはそれぞれ速度ベクトルを表す。Ｖｘ、ＶｙおよびＶｚの合成ベクトルの方向は、弾道弾３の移動方向に相当する。
【００８６】
標定処理部１０５は、観測残差の共分散の予測値を算出する（ステップＳ４１）。観測残差の共分散の予測値は、以下のように算出される。
【００８７】
【数１８】
2018151300.tif 000020
【００８８】
ここで、Ｈ（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）における観測行列を表す。
【００８９】
標定処理部１０５は、予測値の修正量（すなわちカルマンゲイン）を算出する（ステップＳ４３）。カルマンゲインは、以下のように算出される。
【００９０】
【数１９】
2018151300.tif 000021
【００９１】
ここで、Ｋ（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）におけるカルマンゲインを表す。Ｓ（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）における観測残差の共分散を表す。
【００９２】
標定処理部１０５は、推定位置および推定速度を算出し（ステップＳ４５）、推定位置および推定速度を標定結果格納部１１１に格納する。推定位置および推定速度は、以下のように算出される。
【００９３】
【数２０】
2018151300.tif 000022
【００９４】
ここで、ｚ（ｔ＋Δｔ）は時刻（ｔ＋Δｔ）における観測値（目標方位の情報を含む）を表し、センサ信号処理部１０３から取得される。
【００９５】
標定処理部１０５は、システム誤差の共分散の推定値を算出する（ステップＳ４７）。そして処理は呼び出し元に戻る。システム誤差の共分散の推定値は、以下のように算出される。
【００９６】
【数２１】
2018151300.tif 000023
【００９７】
ここで、Ｐバー（ｔ＋Δｔ｜ｔ＋Δｔ）は、時刻（ｔ＋Δｔ）におけるシステム誤差の共分散の推定値を表す。
【００９８】
このように、上で述べたような状態方程式および観測方程式に基づく拡張カルマンフィルタ（拡張カルマンフィルタはベイズフィルタの一種である）を使用すれば計算量が少なくなるので、サンプリング間隔が短くなりサンプル数が増えたとしても、計算量が増大することを抑制できるようになる。なお、拡張カルマンフィルタは、特に物理法則に従った運動に対して適用することで収束までの時間を短くすることができる技術である。
【００９９】
図７の説明に戻り、標定処理部１０５は、ブーストフェーズが終了してから所定時間（例えば２００秒）が経過したか判定する（ステップＳ２３）。
【０１００】
所定時間が経過していない場合（ステップＳ２３：Ｎｏルート）、処理はステップＳ１９に戻る。所定時間が経過した場合（ステップＳ２３：Ｙｅｓルート）、標定処理部１０５は、標定処理が終了したことを着弾予測部１０７に通知する。これに応じ、着弾予測部１０７は、着弾予測処理を実行する（ステップＳ２５）。着弾予測処理については、図１４及び図１５を用いて説明する。
【０１０１】
まず、着弾予測部１０７は、標定処理の結果（すなわち、標定処理において最後に算出された推定位置および推定速度）を標定結果格納部１１１から読み出す（図１４：ステップＳ５１）。
【０１０２】
着弾予測部１０７は、ステップＳ５１において読み出した推定位置および推定速度を用いて、弾道弾３の予測位置および予測速度を算出する（ステップＳ５３）。ステップＳ５３においては、式（１４）に基づき予測位置および予測速度が算出される。但し、ｔは、標定結果格納部１１１から読み出された推定位置および推定位置に対応する時刻（以下、時刻ｔs）である。なお、Δｔの初期値は所定の単位時間であり、ステップＳ５３の処理が実行される度に２Δｔ、３Δｔ、４Δｔ、・・・、と更新される。
【０１０３】
着弾予測部１０７は、弾道弾３の予測位置の高度が０より小さいか判定する（ステップＳ５５）。弾道弾３の予測位置の高度が０より小さくない場合（ステップＳ５５：Ｎｏルート）、着弾予測部１０７はΔｔを更新する。そして処理はステップＳ５３に戻る。
【０１０４】
一方、弾道弾３の予測位置の高度が０より小さい場合（ステップＳ５５：Ｙｅｓルート）、着弾予測部１０７は、予測着弾範囲を算出する（ステップＳ５７）。
【０１０５】
着弾予測部１０７は、ステップＳ５７において算出した予測着弾範囲のデータを出力データ格納部１１３に格納する（ステップＳ５９）。そして処理は呼び出し元に戻り終了する。
【０１０６】
なお、予測着弾範囲のデータは、例えば地上に有るコンピュータに送信されてその表示装置等に表示され、管理者により確認される。
【０１０７】
図１５は、予測着弾範囲の算出について説明するための図である。予測位置の高度が０になった場合における弾道弾３の位置が予測着弾位置である。予測着弾範囲は、予測着弾位置を中心とする誤差楕円体を地表面に対して速度ベクトル方向に射影した場合における楕円に相当する。なお、時刻ｔfは予測着弾時刻である。
【０１０８】
次に、図１６乃至図２７を用いて、本実施の形態の方法に基づくシミュレーションについて説明する。
【０１０９】
図１６は、本シミュレーションにおける弾道弾３を示す図である。弾道弾３は１３００ｋｍ級の中距離弾道ミサイルである。弾道弾３の全長は１６．０（ｍ）であり、幅は１．３２（ｍ）である。断面積は１３７００（ｃｍ2）であり、抵抗係数は０．１である。ペイロードの質量は１０００（ｋｇ）であり、発射時の総質量は１６２５０（ｋｇ）である。一点鎖線で表された部分はプルームを表す。
【０１１０】
図１７は、飛翔シナリオを示す図である。本シミュレーションにおいては、３つのシナリオが試行される。観測機１ａについてのシナリオをシナリオＡとし、観測機１ｂについてのシナリオをシナリオＢとし、観測機１ｃについてのシナリオをシナリオＣとする。各観測機の飛翔高度は２５０００フィートである。観測機１ａ乃至１ｃは、弾道弾３が初めて検知された位置からの距離が８００（ｋｍ）の位置に滞空する。但し、観測機１ａは弾道弾３の軌道の直下に滞空し、観測機１ｂは軌道の直下から２００（ｋｍ）離れた位置に滞空し、観測機１ｃは軌道の直下から４００（ｋｍ）離れた位置に滞空する。
【０１１１】
シミュレーション条件は、図１８及び図１９に示すとおりである。観測機１ａ乃至１ｃから弾道弾３が初めて検知された位置までの距離は、上で述べたように８００（ｋｍ）であるが、図１９に示すように４００（ｋｍ）の相対距離誤差が重畳されるので、距離Ｒは１２００（ｋｍ）に設定される。
【０１１２】
図２０は、本シミュレーションの処理フローを示す図である。ここでは、観測機１ａの動作を例として説明するが、観測機１ｂおよび観測機１ｃも同様に動作する。
【０１１３】
観測機１ａは、飛行経路を生成し（ステップＳ１９０１）、また、目標である弾道弾３の軌道（以下、目標軌道と呼ぶ）を生成する（ステップＳ１９０２）。
【０１１４】
観測機１ａは、生成された飛行経路および目標軌道に基づき、観測角（観測角は、上で述べた目標方位に相当する）を計算する（ステップＳ１９０３）。ステップＳ１９０３において計算された観測角は真値として取り扱われる。
【０１１５】
観測機１ａは、生成された飛行経路から観測機１ａの位置誤差を生成する（ステップＳ１９０４）。また、観測機１ａは、真値である観測角から測角誤差を生成する（ステップＳ１９０５）。
【０１１６】
観測機１ａは、観測機１ａの観測機位置（誤差が含まれる）および観測角（誤差が含まれる）に基づき、標定処理を実行する（ステップＳ１９０６）。
【０１１７】
ステップＳ１９０６の標定処理により生成された標定結果と、目標軌道に基づく目標位置の真値とから、標定誤差（本実施の形態においては、位置のずれ）が算出される。以上が本シミュレーションの処理フローである。
【０１１８】
図２１は、各シナリオについて計算結果が収束するまでの時間を示す図である。サンプリング間隔は５００（ミリ秒）である。表における上段は特許文献１に記載の最小二乗法のアルゴリズムに基づき同一のシミュレーション条件でシミュレーションを実行した場合の時間を示し、下段は本実施の形態の方法に基づくシミュレーションを実行した場合の時間を示す。図２１に示すように、本実施の形態においては全シナリオにおいて時間が短縮されており、特にシナリオＢおよびシナリオＣについては約５割時間が短縮された。
【０１１９】
また、１サイクルの計算量（ここでは、標定処理を１回実行した場合における計算量）については、本実施の形態の方法の実行した場合の計算量は特許文献１に記載の最小二乗法のアルゴリズムを使用した場合の計算量と比較して大幅に少なくなることが確認された。
【０１２０】
図２２は、サンプリング間隔毎に収束までの時間を示す図である。上段はサンプリング間隔が５００（ミリ秒）である場合の時間を示し、下段はサンプリング間隔が１００（ミリ秒）である場合の時間を示す。本実施の形態の方法により計算量が少なくなるので、サンプリング間隔を短くすることができる。そして、サンプリング間隔を短くすることによって、収束するまでの時間が約１割短くなることが確認された。
【０１２１】
図２３は、標定誤差と観測時間との関係を示す図である。図２３はシナリオＡについてのグラフを示しており、サンプリング間隔は１００（ミリ秒）である。縦軸は標定誤差を表し、横軸は観測時間を表す。図２３に示すように、観測時間が１００秒より長くなると急激に標定誤差が小さくなり、観測時間が約１７０秒になった場合に標定誤差が１０（ｋｍ）以下になることが確認された。
【０１２２】
図２４は、システム誤差の初期共分散の設定値を変化させた場合における収束までの時間を示す図である。サンプリング間隔は１００（ミリ秒）である。図２４の表の上段は第１の設定値でシミュレーションを実行した場合の時間を示し、図２４の表の下段は第１の設定値より大きい第２の設定値でシミュレーションを実行した場合の時間を示す。設定値を大きくすることで、実際の観測値の分散との差分が大きくなるため、位置および速度の修正量が大きくなる。そのため、シナリオＡについては、第２の設定値を使用することで収束までの時間が短くなることが確認された。但し、シナリオＢ及びシナリオＣについては、結果はほとんど同じであった。
【０１２３】
図２５は、第２の設定値を使用した場合における標定誤差を示す図である。図２５はシナリオＡについてのグラフを示しており、サンプリング間隔は１００（ミリ秒）である。縦軸は標定誤差を表し、横軸は観測時間を表す。第２の設定値を使用した場合には、観測時間が１００（秒）になる前において標定誤差が急激に小さくなることが確認された。
【０１２４】
図２６は、着弾位置の予測誤差を示す図である。サンプリンク間隔は１００（ミリ秒）である。図２６の表の上段は特許文献１に記載の最小二乗法のアルゴリズムに基づき同一のシミュレーション条件でシミュレーションを実行した場合の予測誤差を示し、下段は本実施の形態の方法に基づくシミュレーションを実行した場合の予測誤差を示す。各シナリオにおいて、顕著な精度向上は確認されなかった。
【０１２５】
図２７は、着弾位置の予測誤差と観測時間との関係を示す図である。図２７はシナリオＡについてのグラフを示しており、サンプリング間隔は１００（ミリ秒）である。縦軸は着弾位置の予測誤差を表し、横軸は観測時間を表す。図２７に示すように、観測時間が約１３０秒になった場合に着弾位置の予測誤差が１０（ｋｍＣＥＰ）以下になることが確認された。
【０１２６】
なお、ベイズ推定では、真値が確率変数として取り扱われ、サンプルをもとに真値の分布が更新される。そのため、サンプルの分布を予め想定できるのであれば、設計時に盛り込むことによって、収束を早くすることを期待できる。本実施の形態では、システム誤差の初期共分散および観測残差の初期共分散を予め適切に設定することにより、計算量を少なくすることができる。また、最小二乗法に基づくアルゴリズムを使用した場合にはサンプル数に比例して１サイクルの計算量が多くなることがあるが、本実施の形態の方法であれば、サンプル数にかかわらず計算量を一定にすることができる。
【０１２７】
以上より、サンプリング間隔を短くすることができるので、同一の観測期間内に取得できるサンプルの数を増やすことができるようになる。これにより、より早い段階で弾道弾３の着弾位置を予測することができるようになる。
【０１２８】
以上本発明の一実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、上で説明した情報処理装置１０の機能ブロック構成は実際のプログラムモジュール構成に一致しない場合もある。
【０１２９】
また、処理フローにおいても、処理結果が変わらなければ処理の順番を入れ替えることも可能である。さらに、並列に実行させるようにしても良い。
【０１３０】
なお、弾道弾３以外の飛翔体に対しても本実施の形態を適用することが可能である。また、慣性系の座標ではなく非慣性系の座標で計算を行ってもよい。
【０１３１】
以上述べた本発明の実施の形態をまとめると、以下のようになる。
【０１３２】
本実施の形態の第１の態様に係る位置推定方法は、（Ａ）飛翔体の観測情報を取得し、（Ｂ）共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを観測情報に適用して、飛翔体の推定位置を算出する処理を含む。
【０１３３】
上記のような非線形のベイズフィルタを用いているので、飛翔体の位置に関する計算の結果が収束するまでに要する時間を短縮することができるようになる。
【０１３４】
また、第１の所定値は、全特異点の留数の総和に基づく値であってもよく、第２の所定値は、観測情報において発生しうる最大の誤差に基づく値であってもよい。
【０１３５】
妥当な共分散が与えられることにより収束までの時間が短縮されるようになる。
【０１３６】
また、非線形のベイズフィルタは、拡張カルマンフィルタであってもよい。
【０１３７】
また、推定位置を算出する処理において、（ｂ１）観測情報または前回算出された推定値に基づき、非慣性系での運動方程式に従って飛翔体の予測位置を算出し、（ｂ２）予測位置と観測情報とに基づき、カルマンゲインを算出し、（ｂ３）予測位置とカルマンゲインとに基づき、推定位置を算出してもよい。
【０１３８】
例えばコリオリ力および遠心力等を推定位置に反映できるようになる。
【０１３９】
また、予測位置を算出する処理において、（ｂ１１）予測位置の座標系を回転座標系から局所座標系に変換し、推定位置を算出する処理において、（ｂ３１）局所座標系での予測位置とカルマンゲインとに基づき、局所座標系での推定位置を算出してもよい。
【０１４０】
位置座標の値が適切に算出されるようになる。
【０１４１】
また、観測情報は、赤外線検出装置により取得された、飛翔体を含む二次元画像の情報を含んでもよい。そして、推定位置を算出する処理において、（ｂ３１１）少なくとも二次元画像の情報に基づき、飛翔体の方位を算出し、方位に基づき、飛翔体の観測位置を算出してもよい。
【０１４２】
飛翔体が遠方に位置する場合であっても、赤外線装置により飛翔体を捉えることができる場合がある。従って、赤外線検出装置により取得された二次元画像の情報（例えば、飛翔体の画素位置の情報）を用いることで、遠方の飛翔体の位置を推定できるようになる。
【０１４３】
また、推定位置を算出する処理において、（ｂ４）観測情報から飛翔体の推定速度をさらに算出してもよい。そして、本位置推定方法は、（Ｃ）推定位置および推定速度から、飛翔体の落下位置の範囲を算出し、（Ｄ）落下位置の範囲の情報を出力する処理をさらに含んでもよい。
【０１４４】
飛翔体の落下への対策を講じることができるようになる。
【０１４５】
また、システム誤差の共分散は、三次元位置の各成分および三次元速度の各成分についての要素を含む６行６列の行列で表されてもよい。
【０１４６】
また、状態方程式は、地球の中心を原点とした直交座標系での位置および速度についての項を含んでもよく、観測方程式は、局所座標系での角度についての項を含んでもよい。
【０１４７】
本実施の形態の第２の態様に係る位置推定装置は、（Ｅ）飛翔体の観測情報を取得する取得部（実施の形態の捜索センサ１００１および追尾センサ１００２は、取得部の一例である）と、（Ｆ）共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを取得部により取得された観測情報に適用して、飛翔体の推定位置を算出する算出部（実施の形態の標定処理部１０５は、算出部の一例である）とを有する。
【０１４８】
なお、上記方法による処理をコンピュータに行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
【０１４９】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【０１５０】
（付記１）
コンピュータに、
飛翔体の観測情報を取得し、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記観測情報に適用して、前記飛翔体の推定位置を算出する、
処理を実行させる位置推定プログラム。
【０１５１】
（付記２）
前記第１の所定値は、全特異点の留数の総和に基づく値であり、
前記第２の所定値は、前記観測情報において発生しうる最大の誤差に基づく値である、
付記１記載の位置推定プログラム。
【０１５２】
（付記３）
前記非線形のベイズフィルタは、拡張カルマンフィルタである、
付記１又は２記載の位置推定プログラム。
【０１５３】
（付記４）
前記推定位置を算出する処理において、
前記観測情報または前回算出された推定値に基づき、非慣性系での運動方程式に従って前記飛翔体の予測位置を算出し、
前記予測位置と前記観測情報とに基づき、カルマンゲインを算出し、
前記予測位置と前記カルマンゲインとに基づき、前記推定位置を算出する、
付記３記載の位置推定プログラム。
【０１５４】
（付記５）
前記予測位置を算出する処理において、
前記予測位置の座標系を回転座標系から局所座標系に変換し、
前記推定位置を算出する処理において、
前記局所座標系での前記予測位置と前記カルマンゲインとに基づき、前記局所座標系での前記推定位置を算出する、
付記４記載の位置推定プログラム。
【０１５５】
（付記６）
前記観測情報は、赤外線検出装置により取得された、前記飛翔体を含む二次元画像の情報を含み、
前記推定位置を算出する処理において、
少なくとも前記二次元画像の情報に基づき、前記飛翔体の方位を算出し、
前記方位に基づき、前記飛翔体の観測位置を算出する、
付記４記載の位置推定プログラム。
【０１５６】
（付記７）
前記推定位置を算出する処理において、
前記観測情報から前記飛翔体の推定速度をさらに算出し、
前記コンピュータに、
前記推定位置および前記推定速度から、前記飛翔体の落下位置の範囲を算出し、
前記落下位置の範囲の情報を出力する、
処理をさらに実行させる付記１記載の位置推定プログラム。
【０１５７】
（付記８）
前記システム誤差の共分散は、三次元位置の各成分および三次元速度の各成分についての要素を含む６行６列の行列で表される、
付記２記載の位置推定プログラム。
【０１５８】
（付記９）
前記状態方程式は、地球の中心を原点とした直交座標系での位置および速度についての項を含み、
前記観測方程式は、局所座標系での角度についての項を含む、
付記１記載の位置推定プログラム。
【０１５９】
（付記１０）
コンピュータが、
飛翔体の観測情報を取得し、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記観測情報に適用して、前記飛翔体の推定位置を算出する、
処理を実行する位置推定方法。
【０１６０】
（付記１１）
飛翔体の観測情報を取得する取得部と、
共分散の初期値が第１の所定値に設定されたシステム誤差の項を含む状態方程式と、共分散の初期値が第２の所定値に設定された観測残差の項を含む観測方程式とに基づく非線形のベイズフィルタを前記取得部により取得された前記観測情報に適用して、前記飛翔体の推定位置を算出する算出部と、
を有する位置推定装置。
【符号の説明】
【０１６１】
１  観測機    １０  情報処理装置
１５１  ＣＰＵ    １５３  メモリ
１５５  ＨＤＤ    １０００  ＧＰＳアンテナ
１００１  捜索センサ    １００２  追尾センサ
１００３  ジンバル    １００４  位置方位基準器
１０１  ジンバル制御部    １０３  センサ信号処理部
１０５  標定処理部    １０７  着弾予測部
１１１  標定結果格納部    １１３  出力データ格納部
【図１】
2018151300.tif 000024
【図２】
2018151300.tif 000025
【図３】
2018151300.tif 000026
【図４】
2018151300.tif 000027
【図５】
2018151300.tif 000028
【図６】
2018151300.tif 000029
【図７】
2018151300.tif 000030
【図８】
2018151300.tif 000031
【図９】
2018151300.tif 000032
【図１０】
2018151300.tif 000033
【図１１】
2018151300.tif 000034
【図１２】
2018151300.tif 000035
【図１３】
2018151300.tif 000036
【図１４】
2018151300.tif 000037
【図１５】
2018151300.tif 000038
【図１６】
2018151300.tif 000039
【図１７】
2018151300.tif 000040
【図１８】
2018151300.tif 000041
【図１９】
2018151300.tif 000042
【図２０】
2018151300.tif 000043
【図２１】
2018151300.tif 000044
【図２２】
2018151300.tif 000045
【図２３】
2018151300.tif 000046
【図２４】
2018151300.tif 000047
【図２５】
2018151300.tif 000048
【図２６】
2018151300.tif 000049
【図２７】
2018151300.tif 000050