航空管制における高信頼度目標追跡法に関する研究

全文

(1)九州大学学術情報リポジトリ Kyushu University Institutional Repository. 航空管制における高信頼度目標追跡法に関する研究高林, 佑樹. https://doi.org/10.15017/1543998 出版情報：Kyushu University, 2015, 博士（工学）, 課程博士バージョン：権利関係：Fulltext available..

(2) 航空管制における高信頼度目標追跡法に関する研究. 2015 年. 高林佑樹.

(3) 目. 次. 第 1 章序論 ....................................................................................................... 1 1.1 本研究の背景 ............................................................................................................. 2 1.1.1 TDOA 測位システムの概要 ........................................................................... 3 1.1.2 モノパルスレーダによる目標高度推定の概要 .......................................... 6 1.2 本研究の目的 ............................................................................................................. 8 1.3 本論文の構成 ............................................................................................................. 9. 第 2 章目標追尾技術 .................................................................................... 10 2.1 緒言 ........................................................................................................................... 10 2.2 追尾技術の概要 ....................................................................................................... 10 2.2.1 追尾フィルタ.................................................................................................. 11 2.2.2 相関処理 .......................................................................................................... 15 2.2.3 初期化処理 ...................................................................................................... 17. 第 3 章非同期測位方式 ................................................................................ 18 3.1 緒言 ........................................................................................................................... 18 3.2 従来の TDOA 測位 .................................................................................................. 19 3.3 ASREX ...................................................................................................................... 20 3.3.1 TDOA 観測時の処理式詳細 ......................................................................... 26 3.3.2 TDOA と FDOA 観測時の処理式詳細 ........................................................ 27 3.4 計算機シミュレーションによる性能評価 ......................................................... 30 3.4.1 シミュレーション条件 ................................................................................. 30 3.4.2 シミュレーション結果 ................................................................................. 34 3.4.3 初期値妥当性.................................................................................................. 42 i.

(4) 3.5 結言 ........................................................................................................................... 45. 第 4 章非同期追尾方式 ................................................................................ 46 4.1 緒言 ........................................................................................................................... 46 4.2 従来手法と課題 ....................................................................................................... 47 4.3 ASTRA ...................................................................................................................... 49 4.3.1 目標の運動モデル ......................................................................................... 49 4.3.2 センサの観測モデル ..................................................................................... 50 4.3.3 予測と平滑 ...................................................................................................... 53 4.4 計算機シミュレーションによる性能評価 ......................................................... 55 4.4.1 シミュレーション条件 ................................................................................. 55 4.4.2 シミュレーション結果 ................................................................................. 58 4.5 結言 ........................................................................................................................... 70 4.6 付録 ........................................................................................................................... 71 4.6.1 安定性の評価.................................................................................................. 71 4.6.2 観測精度 .......................................................................................................... 73. 第 5 章マルチパス環境下における高度推定方式 ................................... 74 5.1 緒言 ........................................................................................................................... 74 5.2 マルチパス伝搬モデル .......................................................................................... 75 5.3 複数固定仮高度型推定方式 .................................................................................. 79 5.3.1 追尾処理 .......................................................................................................... 80 5.3.2 マルチパス高度計算 ..................................................................................... 81 5.3.3 仮高度信頼度計算 ......................................................................................... 82 5.3.4 最適仮高度選択処理 ..................................................................................... 83 5.4 粒子フィルタ型推定方式 ...................................................................................... 84 5.4.1 仮高度重み計算 ............................................................................................. 86 5.4.2 リサンプリング ............................................................................................. 86 5.4.3 高度推定処理.................................................................................................. 87 5.5 複数固定仮高度型推定方式の性能評価 ............................................................. 88 ii.

(5) 5.5.1 シミュレーション条件 ................................................................................. 89 5.5.2 シミュレーション結果 ................................................................................. 90 5.6 粒子フィルタ型推定方式の性能評価.................................................................. 96 5.6.1 シミュレーション条件 ................................................................................. 96 5.6.2 シミュレーション結果 ................................................................................. 97 5.7 結言 ......................................................................................................................... 100. 第 6 章結論 ................................................................................................... 101 謝辞 ................................................................................................................ 103 参考文献 ........................................................................................................ 104. iii.

(6) 第 1 章序論. 第1章序論. 航空交通管制システムは，航空交通量の増加とともに発展してきた．航空機が発明された当時の航空交通量は多くなかったため，パイロットは現在と比べ自由に飛行することが可能であった．しかし航空機の運航数が増加すると，航空機同士の衝突を防ぎつつ航空機を効率的に飛行させるための交通整理として，航空交通管制の必要性が生じた．国土交通省航空局の統計によると，航空機の運航便数は 1998 年の 88 万便に対して 2012 年では 130 万便であり，過去 15 年間で約 1.5 倍まで増大している．さらに 2025 年頃までには 166 万便まで増加すると予測されている[1]．今後，交通量の増大に伴い，航空機離着陸の遅延が起こり，最悪ダイヤ通りに運航ができないことによって欠航となるケースが発生する状況が想定されている．斯かる交通量増大に対する施策の一例として，混雑空港における同時並行離着陸方式の導入が検討されている[1][2]．本運用方式は，従来は１機ごとの離着陸に対して２機同時に離着陸をすることによって運航の効率性が向上するものである．また，他の施策として，航空路上の航空機同士の飛行間隔を短縮し，効率的に運航させる広域航法と呼ばれる航法も検討されている[3]．ただし，これらの運用を実現するためには安全性の確保が最も重要であり，航空機を高精度かつ高頻度に監視することが必要となる．航空交通管制システムにおいて，航空管制官は航空機の安全運航や定時運航を目的として，航空交通の監視や判断，パイロットへの指示など主要な機能を担っている[4]．航空管制官は航空機同士の衝突を避けるため，将来の航空交通状況を予測し，各航空機に対して適切な指示を与える必要がある．そのため，航空機の位置のみならず，速度や加速度などの動態情報を推定する追尾技術が用いられている．さらに，追尾技術は，航空機などの移動体（目標）の運動緒元（位置，速度，加速度など）を推定するだけでなく，航空機以外の雲などの不要信号を除去する機能を有しており，航空管制システムにおいて重要な技術の一つである．本論文は，航空交通管制における目標追尾手法に関する研究を纏めたものである．本章では，まず研究の背景を示した後，本研究の目的と課題について述べ，最後に本論文の構成を示す．. 1.

(7) 第 1 章序論. 1.1 本研究の背景航空管制システムにおいて，これまでレーダによる航空機監視が行われてきた．レーダとは，自ら電波を空間に放射し，航空機等の目標からの反射波を探知することで目標位置を観測するアクティブセンサである．レーダは装置規模が大きく，高価であるものの，昼夜，天候にかかわらず，遠距離の目標を観測できる特長を持つ．主に航空管制で用いられるレーダは，開口アンテナによって，ビーム幅が上下に広く，方位面内で狭い，いわゆるファンビームを形成する２次元レーダであり，機械式に 360°回転することで目標の距離，方位角を観測することができる．また，航空管制用レーダには，空港周辺の航空機を監視する空港監視レーダ(ASR: Airport Surveillance radar)と航空路上の航空機を監視する航空路監視レーダ (ARSR: Air Route Surveillance Radar)等がある．ASR，ARSR は，一般的に１次監視レーダ (PSR: Primary Surveillance Radar)と２次監視レーダ(SSR: Secondary Surveillance Radar)が組み合わされて構成される．PSR は回転する指向性の PSR アンテナから電波を送受して目標の距離，方位角を観測するレーダである．SSR は PSR アンテナと共に回転する指向性アンテナから質問信号を送信して，航空機に搭載されたトランスポンダ（レーダ信号応答装置）から送り返されてくる応答信号を処理して航空機の気圧高度及び識別情報を得るレーダである．ただし，現行の航空管制用レーダである ASR は 4 秒に 1 回，ARSR は 10 秒に 1 回の間隔で回転するため，航空機位置の更新頻度はそれぞれ 4 秒，10 秒と長く，レーダの性質として遠距離ほど方位方向に対しての精度が劣化するため，同時並行離着陸方式や管制間隔短縮化等の高度な運用方式を導入する際の課題となっている．また，現行のレーダでは，トランスポンダが搭載されていない，もしくは故障した航空機からは気圧高度情報が取得できないため，正確な高度が不明であることも課題となっている．特に，テロ機などの非協調型の航空機について正確な位置が取得できないことが問題視されている．斯かる課題に対して，近年，上記の高度な運用方式を導入するために，広域マルチラテレーション（WAM: Wide Area Multilateration）による到来時間差（TDOA: Time Deference of Arrival）測位が検討されている[5][6]．また，トランスポンダ未搭載の航空機高度情報を取得するために，モノパルス測角方式を採用した３次元レーダ（以下，モノパルスレーダ）の導入が検討されている[7]．. 2.

(8) 第 1 章序論. そこで，本論文は上記２つの方式に着目し，TDOA 測位における追尾性能改善およびモノパルスレーダにおける高度推定精度の改善により，高精度かつ高信頼性を備えたレーダ管制を実現するものである．以下，上記２つの方式の概要及び課題を示す．. 1.1.1. TDOA 測位システムの概要. 近年，新しい監視システムとして，空港面ではマルチラテレーション（ MLAT: MultiLATeration），空港周辺及び航空路では広域マルチラテレーション（WAM: Wide Area Multilateration）の導入が検討されている．MLAT，WAM は航空機トランスポンダが送信する信号を地上に配置された複数の受信局（パッシブセンサ）で検出して，その検出時間差から航空機の位置を算出するシステムである．航空機トランスポンダからの応答信号を利用することで，更新頻度が平均 1 秒と航空管制用レーダの 4 秒もしくは 10 秒周期に比べて短い特長を持つ．また，MLAT，WAM は到来時間差（TDOA: Time Deference of Arrival）測位の原理を利用することで，単一レーダによる測位に比べて精度が測角精度に影響されず遠距離での精度が高いといった特長を持つ．. 1.1.1.1. TDOA 測位の原理. TDOA 測位の原理について説明する．TDOA 測位では，航空機に搭載されたトランスポンダが送信する SSR の応答信号を複数の受信局で検出して電波到来時刻を測定する．次に，測定した電波到来時刻から受信局間の到来時間差を求めて，航空機と各受信局との距離差に変換する．そして，「距離差が一定」との条件からなる双曲面同士の交点を求めることで航空機の位置を算出する（図 1.1）．. 3.

(9) 第 1 章序論. f34 Sensor 1. f12. Sensor 4. Target. TDOA f23. Sensor 2. Sensor 3. 図 1.1. 1.1.1.2. TDOA 測位の概念. 従来研究の概要と課題. 従来の TDOA 測位（図 1.2）では，GPS 測位[28]と同様の原理で各受信局の位置と速度を既知とした上で，目標の位置や速度を未知とする測位方程式を逐次最小 2 乗法(ILS: Iterated Least Squares) [22]-[24][26][28]を用いて解くことにより，位置と速度からなる測位結果を得る．ただし，測位結果を得るためには，最低でも未知数＋1 個のセンサ数が必要である．なぜならば，TDOA 観測値を用いた測位は，図 1.1 に示すように各センサ組と目標間の TDOA から得られる双曲面の交点を求めることに相当するため，例えばセンサ位置を既知として，3 次元の目標位置[x,y,z]を測位する場合，最低 4 つのセンサを要することになる．仮にセンサ数を 3 つとした場合には，センサ 1 と 2，センサ 1 と 3 による最大 2 つの TDOA 観測値のみ入力されることとなり，測位方程式より 3 次元位置の解を得ることはできない（センサ 2 と 3 による TDOA はセンサ 1 と 2 及びセンサ 1 と 3 による TDOA の差で表現されるため，1 次従属となる）．ゆえに，一時的に観測値が欠落するセンサが発生する場合，同時刻に測位方程式の数が揃わず，測位結果が得られなくなる問題がある（図 1.3）．センサ観測値が欠落する要因の一つとして，目標からの電波の反射以外に構造物等からの反射など２つ以上の電波伝播経路が発生する（マルチパスと呼ぶ）場合に起こるマルチパスフェージング現象（干渉性フェージングとも呼ばれる）によって，電波強度が低下することが挙げられる．マルチ. 4.

(10) 第 1 章序論. パスフェージング現象とは，直接波と間接波が干渉し合い，電波を強めあったり弱めあったりする現象である．また，他の要因として，構造物などで電波が遮蔽される可能性もある．先行研究ではこのうちマルチパス対策として，二階層にわたって粒子フィルタによる異常値除去と屋内移動体の状態量推定を行っている[17]．ただし，屋内移動体の速度が小さく，距離情報の時間変化も緩やかである仮定のもとで近似して位置を測位しているため，自由空間で高速に移動する目標に対しては近似が成立せずに性能に限界がある．一方，他の先行研究では追尾フィルタで得られる TDOA の予測値からマルチパス干渉下にあり，精度が劣化したセンサを推定し，マルチパス干渉下でないセンサの TDOA 観測値を用いて測位する方式について報告されている[18]．しかし，マルチパス干渉下にないセンサ数が測位に必要な数に満たない場合については，追尾フィルタは動作しないため性能が劣化する．また，他の先行研究[20][21]も文献[18]と同様，同時刻に測位に必要なセンサ数が得られることが前提であり，必要なセンサ数が得られない場合の検討はなされていない．つまり，従来の TDOA 測位[17][18][20][21]では，同一時刻に得られる各センサ組の観測値が測位に必要な数に満たない場合，測位不能となる．特に，後述するように追尾フィルタとしてカルマンフィルタを用いた場合には，追尾開始時に位置や速度等の状態量に適切な初期値を設定する必要がある．しかし，測位が不能となる場合は，初期値を設定できず，追尾開始が遅れることとなる．. Sensor1. TDOA/ FDOA. Sensor2. TDOA/ FDOA. Sensor3. ILS. Residual computation. Gradient computation. TDOA/ FDOA. ・・・. ・・・. Least squares. 図 1.2. 従来測位方式のブロック図. 5. Estimated state vector.

(11) 第 1 章序論. Positioning is impossible f34. Sensor 1. f12. Sensor 4. Target. TDOA f23. Sensor 2. Sensor 3. 図 1.3. 1.1.2. TDOA 観測値欠落時の問題. モノパルスレーダによる目標高度推定の概要. 近年，航空管制の分野において，テロ機等の非協調型の航空機を監視することが求められている．特に，欧州では欧州各国が横断的に参画する研究プロジェクトが立ち上がり，従来の航空管制とは異なるコンセプトを打ち出している[7]．その中で，従来，航空管制で利用されてきた目標の距離，方位角を観測する２次元レーダではなく，３次元レーダを用いて目標高度を含む３次元位置を測位することが示唆されている．３次元レーダでは，高度を観測するためにモノパルス測角方式[38]を想定しており，目標の距離，方位角に加え，仰角を観測することができる．上記の３次元レーダは防衛システムでは既に自国の領空を侵犯する諸外国の航空機を監視する目的で利用されてきた．しかし，海面を低高度で飛行する目標に対してはレーダと目標間の直接波と海面反射による間接波によってマルチパス環境となり，マルチパスフェージング現象が発生することがある．この時，モノパルスアンテナでの仰角観測値（高度観測値）に大きな誤差が生じることが知られている．しかも，上記誤差はレーダ受信機雑音のようなランダム性の誤差とは異なり，レーダと目標の相対位置関係に依存するバイアス性の誤差である．そのため，マルチパスが問題となるような低高度目標追尾時に，正しく目標高度を推定することは，レーダ追尾において重要な課題となっている．. 6.

(12) 第 1 章序論. 1.1.2.1. マルチパス伝搬モデル. 上記のマルチパス環境における仰角観測値の誤差を表現するモデルとして，マルチパス伝搬モデルが知られている．図 1.4 にマルチパス伝搬モデルの概念図を示す．図 1.4 のように，本論文ではマルチパス波は高度方向にのみ存在し，レーダは目標からの直接波および海面を反射する間接波を受信する．マルチパスフェージング特性は間接波の海面反射点における反射係数，及び直接波と間接波の経路差に依存する[39]．海面反射には Specular（鏡面）反射と Diffuse（拡散）反射があるが，特に海面が穏やかな場合には Specular 反射が支配的となる．この Specular 反射は，直接波と相関が強く，複数の伝搬経路で受信した電波が強め合う，もしくは弱め合うことでマルチパスフェージングが生じ，高度観測値の精度が著しく劣化する問題がある．例えば，このような状況では，レーダと目標との距離に応じて高度観測値がスパイク状に変化する，もしくは海面付近に張り付くように観測されることが知られている.. Target Radar θ0 ψg. ψg 図 1.4. 1.1.2.2. Surface. マルチパス伝搬の概念図. 従来研究の概要と課題. マルチパスによる高度観測値への影響に対し，一般的な追尾フィルタの観測モデル[23]では，真値にランダム誤差が付加されるモデルを前提としており，1.1.2 で説明したようなマルチパスフェージングによって発生するバイアス誤差は考慮していない．そのため，実際の高度観測値と，この観測モデルとの差異によって，追尾精度が劣化する[41][42][43]．従来，このようなマルチパスフェージング環境下において，高度推定値の安定化を図る追尾方式が提案されてきた[44][45]．文献[44]では，マルチパス誤差を駆動雑音が付加された時刻間に相関がある誤差と仮定することで，マルチパスの影響を軽減する追尾アルゴリズムが提案されている．. 7.

(13) 第 1 章序論. 文献[45]では，マルチパスフェージングの影響が周波数に依存する特徴を利用して，2 周波を時分割で切り替え，各周波数で計算した追尾航跡を統合することにより，追尾精度の劣化を抑圧する方法が提案されている．しかし，スパイク状の観測誤差については改善されているものの，海面付近に張り付くように観測される高度に生じる誤差については改善が不十分であった．また，文献[50]では，アレーアンテナをサブアレー構成かつ周波数ホッピングによるフェージング対策を施し，高度を正確に推定する方法が提案されている．しかし，レーダの装置構成を大幅に変更する必要があり，本論文で想定しているレーダとは異なる．. 1.2. 本研究の目的. 本研究の目的は，TDOA 測位における追尾性能改善およびモノパルスレーダの目標高度推定精度の改善である．これらを実現するためには下記のような課題がある．. 課題１. TDOA 測位における同時刻に必要なセンサが揃わない場合の追尾開始早期化. 課題２. TDOA 測位における追尾維持中に測位結果が得られない場合の追尾安定化. 課題３. モノパルスレーダによるマルチパス環境下での目標高度バイアス誤差の除去. 本論文では，課題１を解決するために，目標が等速直線運動するとの仮定の下に，異なる時刻に得られた各センサ組の観測値を用いて逐次最小 2 乗法により，追尾開始に必要な初期状態ベクトルを推定する非同期測位法を提案する．ここで，非同期とは異なる時刻に得られた観測値を利用することを意味する．本方法によれば，同一時刻に得られる各センサ組の観測値が測位に必要な数に満たない場合でも，異なる時刻かつ推定時刻より前に得られた観測値も時間を遡って使用することにより測位可能であるため，早期追尾開始が期待できる．また，課題２を解決するために，非同期に得られる各センサの TDOA 観測値から測位処理を介さずに目標の位置，速度を直接推定する追尾方式を提案する．本方式が実現できれば，同時刻に測位に必要なセンサ数が得られない場合でも，残りの有効なセンサの観測値を使用してデータ更新が可能であるため，追尾性能の向上が期待できる． 8.

(14) 第 1 章序論. ここで，将来的にドップラー周波数差(FDOA：Frequency Difference of Arrival)も抽出することを想定して，課題１及び２に対しては，TDOA 観測値と FDOA 観測値を併用した方法も検討する．さらに，課題３を解決するために，モノパルス測角を前提としたマルチパス伝搬モデルに基づく高度推定方式を提案する．事前に想定される目標高度の範囲を仮高度として保持しておき，マルチパス伝搬モデルに基づくバイアス誤差を考慮した仮高度の信頼度を計算し，仮高度の信頼度より高度推定値を算出することで正確な高度推定を図る．なお，事前に想定する仮高度を固定とする複数固定仮高度による高度推定方式，仮高度を固定とせずに Particle Filter を利用して仮高度を時刻ごとにリサンプリングする粒子フィルタ型高度推定方式の２手法を提案する．. 1.3. 本論文の構成. 本論文は，第１章「序論」，第２章「目標追尾技術」，第３章「非同期測位方式」，第４章「非同期追尾方式」，第５章「マルチパス環境下における高度推定方式」，及び第６章「結論」から構成されている．以下に各章の概要を述べる．第２章では，本論文の基本技術である目標追尾技術について概要を述べる．第３章では，複数の受信センサによる測位システム及び従来の測位方式の課題を述べる．さらに斯かる課題に対して非同期測位方式を提案し，性能評価結果を述べる．第４章では複数の受信センサによる測位システムを利用した場合に従来の追尾維持で起こる問題を述べる．さらに斯かる課題に対して非同期追尾方式を提案し，性能評価結果を述べる．第５章では，マルチパス環境下における高度推定の課題を述べる．また，複数固定仮高度による高度推定方式及び粒子フィルタ型の高度推定方式を提案し，性能評価結果を述べる．第６章では，本論文で得られた結果を総括する．. 9.

(15) 第 2 章目標追尾技術. 第2章目標追尾技術. 2.1. 緒言. 目標追尾は，レーダなどのセンサの観測結果に基づき，航空機などの移動体（目標）の運動緒元（位置，速度，加速度など）を推定する技術である．目標追尾は，車両，船舶，航空機の交通路の安全を確保するための交通管制システム，自国の領海，領空を侵犯する諸外国の船舶及び航空機の動きを監視する防衛システムなどにおいて重要な役割を担っている [4][8]．本章では，本論文で扱う目標追尾技術の概要を説明する．. 2.2. 追尾技術の概要. 図 2.1 に目標追尾の処理ブロック図を示す．本節では，追尾フィルタについて述べ，次に相関処理，最後に初期化処理の順に説明を述べる． Initialization. Initial state vector. Tracking filter. Sensor. Radar measurements. Coordinates conversion. Correlation. Smoothing. Prediction. 図 2.1. 目標追尾の処理ブロック. 10. Estimated state vector.

(16) 第 2 章目標追尾技術. 図 2.1 において，センサから目標位置に関する観測値が得られると，追尾処理部は，必要に応じて適当な座標系の位置データに変換する．次に「予測・平滑処理」において，前サンプリング時刻に得られた各追尾目標の位置及び速度などの推定値をもとに，現時刻における予測位置を算出する．この予測値とセンサからの観測値との相対的な位置関係から，各目標にいずれの観測値が対応するかの割り当てを決定する処理が「相関処理」である．割り当て済みの目標については，観測値の情報をもとに，「予測・平滑処理」において，現時刻における新たな推定値が算出される．さらに「相関処理」において目標に割り当てられなかった観測値を新規目標として登録し，「予測・平滑処理」の目標初期値を設定する処理が「初期化処理」である．図 2.1 の「予測・平滑処理」の機能は，センサ観測値に含まれる雑音成分を除去する低域通過フィルタに相当するため，追尾フィルタと呼ばれる．本論文では，追尾フィルタが動作している時間帯を追尾維持中もしくは追尾維持フェーズと定義する．「相関処理」はサンプリング時刻ごとに得られた複数の観測値を追尾目標ごとの信号や不要信号に弁別する機能を担っている．特に，陸地，海面，もしくは雲等から観測されるクラッタのような不要信号が多く得られた場合に重要となる．「相関処理」が失敗すると，不要信号にとらえられて目標を見失うこととなる．「初期化処理」は追尾フィルタに与える初期値（位置，速度など）を計算する機能である．例えば，センサ観測値が目標位置のみの場合，2 サンプル分の観測位置の時間差分から初期速度を算出し，最新のサンプリング時刻における位置を初期位置として出力する．本論文では，この初期値計算が終了した時刻を追尾開始もしくは追尾開始フェーズと定義する．. 2.2.1. 追尾フィルタ. センサから取得した目標位置などの観測値には，受信機雑音に起因する様々な誤差が加わるため，真の目標位置は得られない．追尾フィルタは相関処理によって目標に割り当てられた観測値から誤差を抑圧し，位置のみならず，速度等も含めた目標運動緒元を推定する．追尾フィルタの代表例は，カルマンフィルタ[12][22][23]，αβフィルタ，αβγフィルタ[22][23][46]-[48]を適用したものである．カルマンフィルタは，目標の運動モデルやセンサの観測モデルの設定に基づき，目標運動緒元の最小分散推定値を与える最適フィルタである．また，カルマンフィルタは，例えば，直交座標で定義される x,y,z 軸の各軸成分を一括して扱える干渉型のフィルタである．一方， 11.

(17) 第 2 章目標追尾技術. αβフィルタ，αβγフィルタは，各軸独立に構成される非干渉型のフィルタであり，カルマンフィルタの誤差共分散行列，フィルタゲインの逐次計算を省略し，フィルタゲインに固定値を使用する簡易型のフィルタである．αβフィルタは，目標の位置及び速度を推定し， αβγフィルタは加速度まで推定するフィルタである．これらのフィルタはカルマンフィルタに比べ演算負荷が軽いため，計算機の能力が不十分であった 1980 年代まで多くのシステムに実装されてきた．しかし，その後の計算機性能の向上に伴い，現在ではカルマンフィルタの適用が一般的となっている．以降，目標運動が近似的に等速直線運動とみなせる場合を対象とし，追尾フィルタにカルマンフィルタを適用する際の座標系，目標の運動モデル，センサの観測モデル，予測処理，平滑処理について，一般的な方法を示す[22][23]．なお，ここでは一例として，3 次元レーダを使用する場合を想定して説明する．3 次元レーダとは，目標の距離，仰角，方位角を観測できるレーダである．（１）座標系目標の運動緒元を記述するための座標系として，図 2.2 のような北基準直交座標系が良く知られている．北基準直交座標系は空間に固定した基準点 O を原点とし，東を x 軸の正，北を y 軸の正，鉛直上方を z 軸の正に取った座標系である．北基準直交座標系におけるセンサの位置を“Radar”として，図 2.2 に示す．. Z (Up) z. Target R. Radar  xr , yr , zr . x. E Az. y. 0(origin). X（East）図 2.2. 北基準直交座標 12. Y (North).

(18) 第 2 章目標追尾技術. また，図 2.2 に示すように，センサと目標間の距離を R，水平面より目標までの仰角を E，水平面内で北方向より目標までの方位角を Az とした座標を極座標として定義する．北基準直交座標と極座標には式(2.1)あるいは式(2.2)の関係がある．ここで，北基準直交座標におけるセンサの位置ベクトルを式(2.3)，目標位置ベクトルを式(2.4)で定義する．T は行列の転置を表す．   2 2 2     x  xr    y  yr    z  zr    R     z  zr    1     E    tan  2 2   Az   x  xr    y  yr            1  x  xr  tan      y  yr   . (2.1).  R cos E sin By    x   R cos E cos By   xr   R sin E  . (2.2). xr   xr. (2.3). x  x. zr . T. yr z. T. y. (2.4). （２）運動モデル一般的に航空機などの移動体は，全体の飛行時間に対して，ほぼ一定速度で直進している時間が大半を占める．このため，等速直線運動を線形式で記述できる北基準直交座標系で座標系を定義することが有利とされている．サンプリング時刻 tk における目標の状態ベクトル（運動緒元）を式(2.5)，目標の運動モデルを式(2.6)で定義する． xk   xk. yk. zk. xk. yk. zk . T. (2.5). xk  Φk 1 xk 1  Γ1,k 1wk 1. (2.6). ここで， Φk 1 は時刻 tk 1 より tk への状態ベクトルの既知の推移行列，状態ベクトル xk は北基準直交座標における 3 次元空間の目標位置・速度から成る 6 次元ベクトル， wk は目標が等速直線運動から逸脱した際の誤差を表すための駆動雑音ベクトルで，平均 0，共分散行列. Qk の白色ガウス雑音， Γ1,k は駆動雑音ベクトルの既知の変換行列である．また， Φk 1 は式 (2.7)， Γ1, k 1 は式(2.8)で定義する．. 13.

(19) 第 2 章目標追尾技術. (tk  tk 1 ) I33   I33 . I Φk 1   33 0nn.  (t  t )2 Γ1,k 1   k k 1 I33 2 . (2.7).  (tk  tk 1 ) I33  . T. (2.8). ここで， I nn ， 0nn はそれぞれ n 行 n 列の単位行列，零行列を意味する．また，共分散行列 Qk は目標の機動性に応じて設定されるパラメータである．共分散行列. Qk の対角項の値が大きいほど等速直線運動から目標運動が逸脱した際の追従性能が向上するものの，観測雑音の抑圧効果が下がり，等速直線運動に対する追尾精度が低下する．一方，値が小さいほど等速直線運動に対する追尾精度が向上するものの，運動の変化に対して追従性能が低下する．すなわち，フィルタの帯域幅を決定するパラメータに相当している．. （３）観測モデルサンプリング時刻 tk における極座標系から北基準直交座標系に変換されたセンサの観測ベクトルを式(2.9)で表す．この時，観測モデルは式(2.10)で記述される． zk   xok. zk  H. yok. x. p,k. zok . (2.9).  xr   vk. (2.10). T. ここで，H は状態ベクトルより，位置ベクトルのみを抽出するための既知の観測行列であり，次式で表す．. H   I3 x3. 03 x3 . (2.11). また， v k は北基準直交座標系におけるセンサの観測雑音ベクトルで，平均 0，共分散行列. Bk の白色ガウス雑音である．センサがレーダである場合，極座標系で目標を観測するため，極座標系での観測雑音は座標軸間の相関がないと仮定できるものの，北基準直交座標系では座標軸間に相関が発生する．極座標系の観測雑音ベクトルを v k とすると，北基準直交座標系の観測雑音ベクトル v k は式(2.2)をサンプリング時刻 tk における極座標系の位置ベクトル. pk で線形 1 次近似することにより，式(2.12)で表す．また， Γ 2,k は観測雑音ベクトルの変換行列であり，式(2.13)で表す．極座標系の位置ベクトル pk は式(2.14)で表す．. vk  Γ 2,k vk. Γ 2, k . x p. p  pk. (2.12)  cos Ek sin Byk  cos Ek cos Byk  sin Ek.  Rk sin Ek sin Byk  Rk sin Ek cos Byk Rk cos Ek. 14. Rk cos Ek cos Byk   Rk cos Ek sin Byk   0. (2.13).

(20) 第 2 章目標追尾技術. pk   Rk. Ek. Byk . T. (2.14). 上記式(2.12)，(2.13)の関係より，北基準直交座標系の観測雑音ベクトルの共分散行列 Bk は近似的に次式で表すことができる．. Bk  Γ 2, k.  R2  0 0 . 0. . 2 E. 0. 0   0  Γ 2,T k 2   By . (2.15). ここで，  R ，  E ，  By はそれぞれ距離，仰角，方位角の観測雑音標準偏差である．これらパラメータには通常センサの性能使用値を与える．. （４）予測処理サンプリング時刻 tk における目標の予測ベクトル x p , k ，予測誤差共分散行列 Pp , k は 1 サンプル前の時刻における平滑ベクトル xs , k 1 ，平滑誤差共分散行列 Ps , k 1 より，下記式にしたがって算出される．. x p,k  Φk 1 xs ,k 1. (2.16). Pp,k  Φk 1 Ps ,k 1ΦkT1  Γk 1Qk 1 ΓkT1. (2.17). （５）平滑処理サンプリング時刻 tk における目標の平滑ベクトル xs , k ，平滑誤差共分散行列 Ps , k 及びゲイン行列 K k は予測ベクトル x p , k ，予測誤差共分散行列 Pp , k より，下記式にしたがって算出される．. . x s , k  x p , k  K k zk  H. x. p,k.  xr . . (2.18). Ps ,k   I 66  K k H   Pp,k. K k  Pp,k H T  HPp,k H T  Bk . 2.2.2. (2.19) 1. (2.20). 相関処理. 図 2.1 の「相関処理」に使用される基本的な手法は追尾ゲート内外判定である．追尾ゲートとは，各目標に相関する観測値を限定するために，各目標の予測位置を中心として設定される領域である．追尾ゲート内外判定では，追尾ゲート内に入った観測値のみが当該目標に 15.

(21) 第 2 章目標追尾技術. 対応する観測値候補とされる．一般的に追尾フィルタにカルマンフィルタを適用した場合は，式(2.21)によって追尾ゲート内外判定を実施する．式(2.21)は目標予測位置と観測値のマハラノビス平方距離であり，カイ二乗検定に基づき，観測値が次式を満たす場合に，追尾ゲート内に存在すると判定する．. z. k. H. x. p,k.  xr . . T. . Sk1 zk  H. x. p ,k. .  xr   d. (2.21). ここで S k は残差共分散行列，d はゲートサイズパラメータである[23]．なお，センサが目標探知に失敗し，追尾ゲート内に 1 つも観測値が得られない場合，カルマンフィルタによる前サンプリング時刻における状態ベクトル予測値を現サンプリング時刻まで外挿して推定値を算出し，追尾処理を継続する．この操作はメモリトラック処理と呼ばれる．一方，クラッタなどの不要信号の発生頻度が高い領域や，多くの目標が密集する領域では，追尾ゲート内に複数の観測値が得られて，相関の判定が困難となる場合がある．この場合に対処する方法として，NN（Nearest Neighbor），GNN（Global Nearest Neighbor），JPDA （Joint Probabilistic Data Association），MHT（Multiple Hypothesis Tracking）といった相関処理がある[52][53]．NN は各目標に対して，式(2.21)の左辺が最小となる観測値を選択する方法であり，処理が簡単で実装が容易である反面，複数目標が近接する場合に同一の観測値を割り当てる可能性がある．GNN は観測値の重複割り当てを防ぐために，重複割り当てがないという制約条件の下で，目標の予測位置と観測値のマハラノビス平方距離の総計が最も小さくなるように割り当てを行う．JPDA と MHT は，目標探知失敗の可能性，不要信号発生の可能性，近接目標存在の可能性を考慮した上で，各目標にいずれかの観測値が対応するかについて複数の仮説をたて，各仮説が正しい確率を評価する．さらに，JPDA は各サンプリング時刻において，複数の仮説の影響を重みづけ統合することにより状態推定値を一意に決定する方法であるのに対して，MHT はいくつかの有力な仮説が現れた際に，各仮説を並行して計算し，次サンプリング時刻以降に正しい仮説を選択しようとする方法である．不要信号環境下や近接目標環境下において追尾性能が向上することが確認されている．航空管制の場合は，航空機から識別番号が送られてくることでクラッタとの弁別が容易であり，管制間隔を空けて飛行することで目標が密集することが無いため，本論文の追尾処理には処理負荷が最も軽い NN を用いることとする．. 16.

(22) 第 2 章目標追尾技術. 2.2.3. 初期化処理. 図 2.1 の「初期化処理」では，予測処理，平滑処理が動作する前の平滑ベクトル及び平滑誤差共分散行列の初期値を算出する．ここでは，センサ観測値が観測位置のみである場合について説明する．まず，目標初探知の平滑ベクトルを式(2.22)のように平滑位置を観測位置，平滑速度を零ベクトルとして定義する．なお，初探知では平滑誤差共分散行列は定義しない． xs ,1   z1T. 013 . T. (2.22). 次の観測位置が取得されるサンプリング時刻において，平滑ベクトルを式(2.23)，平滑誤差共分散行列を式(2.24)で定義する．.  xs ,2   z2T . z2T  z1T   t2  t1 .   B2 Ps ,2    B2   t2  t1. T. (2.23).     1 B  B2  Q pp ,1   2  1  t2  t1   B2 t2  t1. (2.24). ここで，上記式の Bk は北基準直交座標系の観測雑音ベクトルの共分散行列（式(2.15)）であり， Q pp , k は下記式(2.25)の Qk の位置に係わる相関項である． Q pp ,k Qk   T Q pv, k. Q pv ,k   Γ1,k Qk Γ1,Tk Qvv, k . (2.25). すなわち，センサ観測値が目標観測位置のみである場合は，少なくとも２探知分の観測値が必要となる．. 17.

(23) 第 3 章非同期測位方式. 第3章非同期測位方式. 3.1. 緒言. 複数の受信センサによる TDOA 観測値を用いた測位では，単一レーダによる測位に比べて精度が測角精度に影響されない特長を持つ反面，マルチパス波や電波遮蔽などの影響により一時的に著しく精度が劣化したセンサが発生することがある[17][18]．ここで，文献[18] では追尾フィルタで得られる TDOA の予測値からマルチパス干渉下にあり，精度が劣化したセンサを推定し，マルチパス干渉下でないセンサの TDOA 観測値を用いて測位する方式について報告されている．しかし，マルチパス干渉下にないセンサの数が測位に必要な数に満たない場合については，追尾フィルタでは推定不可能であるため追尾精度が劣化する．また，文献[20][21]も文献[18]と同様，同時刻に測位に必要なセンサ数が得られることが前提であり，必要なセンサ数が得られない場合の検討はなされていない．つまり，従来の TDOA 測位方式[17]-[21]では，同一時刻に得られる各センサ組の観測値が測位に必要な数に満たない場合，測位不能となる．その際，2.2.3 項で述べた追尾初期化が行われず，追尾開始が不可となる問題がある．この追尾初期化の課題を解決するため，目標が等速直線運動するとの仮定の下に，異なる時刻に得られた各センサ組の観測値を用いて逐次最小 2 乗法により，追尾開始に必要な初期状態ベクトルを推定する非同期測位方式を提案する．ここで，非同期とは異なる時刻に得られた観測値を利用することを意味する．本方法によれば，同一時刻に得られる各センサ組の観測値が測位に必要な数に満たない場合でも，異なる時刻かつ推定時刻より前に得られた観測値も時間を遡って使用することにより測位可能であるため，早期追尾開始が期待できる．また，将来的な拡張を考慮して，提案する非同期測位法に FDOA も利用する方法についても検討する． 18.

(24) 第 3 章非同期測位方式. 以下，本章では，まず初めに従来の測位方式の概要を述べる．次に，提案手法である異なる時刻に得られた TDOA/FDOA 観測値を利用する非同期測位方式（ASREX; ASynchronous TDOA/FDOA localization by Reverse Extrapolation of target state vector）の概要と処理詳細について述べる．最後に，計算機シミュレーションにより，追尾開始性能を評価する．. 3.2. 従来の TDOA 測位. 従来の TDOA 測位（図 1.1）では，GPS 測位[29]と同様の原理で各受信局の位置と速度を既知とした上で，目標の位置や速度を未知とする非線形な測位方程式を逐次最小 2 乗法 (ILS: Iterated Least Squares) [22][23][24][28]を用いて解くことにより，位置と速度からなる測位結果を得る．ここで，従来方式の例として，4 つのセンサで得られる TDOA 観測値のみを用いて，ある時刻における 3 次元の目標位置ベクトル x を測位する場合の測位方程式の解き方を以下に示す．まず，目標とセンサ i，j との距離差観測値を rij，距離差真値を hij  x  ，観測誤差を  ij と定義すると，以下の非線形方程式が成立する．. rij  hij  x    ij , i, j 1,. , 4. (3.1). ここで，目標位置ベクトル x を 3 次元位置として下記式で定義する．また，添え字 T は転置を表す．. x  x. y z. T. (3.2). 目標推定位置ベクトル x̂ とセンサ i，j との距離差 hij ( xˆ ) と観測値 rij との残差は以下のように求められる．. hij  rij  hij ( xˆ ), i, j 1, , 4. (3.3). 例えばセンサ 1 を基準とすると，残差ベクトル δh は以下の式で与えられる．.  h   h21 h31 h41 . T. (3.4). この残差ベクトルを用いて目標推定位置ベクトル x̂ を修正する．修正ベクトルと残差ベクトルの関係は以下の式となる．.  h  G  x. (3.5). ここで，G はヤコビ行列，δx は修正ベクトルである. 19.

(25) 第 3 章非同期測位方式.  h h31 h41  G   21   x x x . T. (3.6) x  xˆ.  x  x  y z . T. (3.7). ILS では，残差ベクトル δh の各成分の 2 乗和が最小となるように，式(3.8)を用いて修正ベクトルを求める．そして，式(3.9)のように目標推定位置ベクトル x̂ を更新する．.  x  GT R1G  GT R1 h. (3.8). xˆ (itr 1)  xˆ (itr )   x. (3.9). 1. ここで，添え字 itr は ILS の反復回数を意味し，観測誤差共分散行列 R は以下のように与える．E[]は期待値を表す． T R  E  21  31  41   21  31  41   . (3.10). 式(3.9)の修正ベクトル δx が十分小さくなるまで反復することにより，目標推定位置ベクトル x̂ が求められる．ただし，前述の測位方程式を解くには，最低でも未知数＋1 個のセンサ数が必要である．なぜならば，TDOA 観測値を用いた測位は，図 1.1 に示すように各センサ組と目標間の TDOA から得られる双曲面の交点を求めることに相当するため，例えばセンサ位置を既知として，3 次元の目標位置[x,y,z]を測位する場合，最低 4 つのセンサを要することになる．仮にセンサ数を 3 つとした場合には，センサ 1 と 2，センサ 1 と 3，センサ 2 と 3 による TDOA が得られるものの，センサ 2 と 3 による TDOA はセンサ 1 と 2 及びセンサ 1 と 3 による TDOA の差で表現されるため，1 次従属であり，測位方程式の数は最大 2 つとなるため， 3 次元位置の解を得ることはできない．ゆえに，構造物等によるマルチパスの発生や電波の遮蔽などにより，一時的に観測値が欠落するセンサが発生する場合，同時刻に測位方程式の数が揃わず，測位結果が得られなくなる問題がある．. 3.3. ASREX. 3.2 節で述べたように従来の測位方式では，同時刻に必要な TDOA/FDOA 観測値が得られない場合，測位不可能である．本節では，異なる時刻で得られた各センサ組の. 20.

(26) 第 3 章非同期測位方式. TDOA/FDOA 観測値を用いて追尾開始に必要な初期状態ベクトルを推定する非同期測位法 ASREX を提案する[25][26]．センサと目標の位置関係を図 3.1 に示す．図 3.1 のセンサ i の位置ベクトルを si，速度ベクトルを si ，目標の位置ベクトルを u とする． z. Target. ri. u. Sensor i. rj. si sj y Sensor j. x. 図 3.1. 目標とセンサの配置. 提案方式では図 3.2 のような非同期で TDOA/FDOA 観測値が得られる状況を想定し，目標が等速直進するという仮定のもとに初期状態ベクトルを推定する．ここで，航空管制において本方式を適用するケースは，遠方から WAM 覆域に入る航空機の追尾開始時である．この場合，航空機は概ね等速直進で接近するため，等速直進運動の仮定を置くことは妥当といえる．. 21.

(27) 第 3 章非同期測位方式. TDOA/FDOA between sensor 1 and 2. TDOA/FDOA between sensor 1 and 3 TDOA/FDOA between sensor 1 and 4. Sensor 3. Sensor 2. t1  t  t1 t2  t  t2. Reverse extrapolation Sensor 4. 図 3.2. t. Sensor 1. ASREX 処理概要. ここで， xˆ (t ) を基準（最新）時刻 t における初期状態ベクトル推定値，その成分 uˆ (t ) ， uˆ (t ) はそれぞれ初期位置ベクトル推定値，初期速度ベクトル推定値とする．また， xˆ k は基. 準時刻 t より過去の時刻 tk(tk＜t)における目標の状態ベクトル推定値， uˆ k は時刻 tk における目標の位置ベクトル推定値， uˆ k は時刻 tk における目標の速度ベクトル推定値， tk を基準時刻 t と時刻 tk の時刻差とすると， xˆ (t ) ， xˆ k ， tk は以下のように定義できる． uˆ (t )T  . xˆ (t )  uˆ (t )T . xˆ k   uˆkT. uˆkT . T. (3.11). T. (3.12). tk  t  tk. (3.13). 目標が 3 次元空間で等速直進すると仮定すると，. xˆ k  xˆ (t  tk )  Φk xˆ (t ) I Φk   33 0nn. (3.14). tk I33  I33 . (3.15). となる．ここで， Φk は状態遷移行列，In×n は n 行 n 列の単位行列を表す．式(3.14)のように，過去の時刻 tk における目標の状態ベクトル推定値 xˆ k を基準時刻 t における初期状態ベク 22.

(28) 第 3 章非同期測位方式. トル推定値 xˆ (t ) で表現することにより，過去の時刻 tk で得られた TDOA/FDOA 観測値も用いて初期状態ベクトル推定値 xˆ (t ) を算出することが可能となる．図 3.2 のように，過去の時刻 t1 にセンサ 1 と 2，過去の時刻 t2 にセンサ 1 と 3，基準時刻 t （=t3）でセンサ 1 と 4 より TDOA と FDOA 観測値が得られる場合を例として，非同期に得られた観測値を用いて目標の初期状態ベクトル推定値 xˆ (t ) を推定する方法について説明する．まず，現時刻 t に対する過去の時刻 tk における目標の状態ベクトル推定値 xˆ k とセンサ i，j との距離差と距離差観測値との残差を Δhij,k，ドップラ速度差とドップラ速度差観測値の残差を Δhij , k とすると，図 3.2 の場合，残差ベクトル δh は以下のように求められる．  h   h21,1 h31,2. h41,3 h21,1 h31,2. h41,3 . T. (3.16). また，ヤコビ行列 G は以下の式となる．  h G   21,1  xˆ (t ). h31,2 xˆ (t ). h41,3 xˆ (t ). h21,1 xˆ (t ). h31,2 xˆ (t ). h41,3   xˆ (t ) . T. (3.17). 従来方式では残差ベクトルの内，基準時刻 t（=t3）と異なる時刻における残差と基準時刻における修正ベクトルには対応関係がなく，式(3.17)の. h21,1 h31,2 h21,1 h31,2 , , , を算出するこ xˆ (t ) xˆ (t ) xˆ (t ) xˆ (t ). とはできない．しかし，提案方式は式(3.14)の運動モデルを用いることで，式(3.17)の全要素を算出できる．例えば，時刻 t2 においてセンサ 1 と 3 で観測される距離差の残差を目標の初期状態ベクトル推定値 xˆ (t ) で偏微分する成分を下記のように表す． h31,2 h31,2 xˆ 2  xˆ (t ) xˆ 2 xˆ (t ). ここで，. (3.18). xˆ 2 は，式(3.14)より求まる．他の要素も同様にして求めることができる． xˆ (t ). 図 3.3 に ASREX の処理ブロック図を示す．従来の測位法が同一時刻で各センサ組の TDOA/FDOA 観測値が必要数得られなければ測位できないのに対し，本手法は同一時刻で TDOA/FDOA 観測値が必要数揃わなくとも，目標の運動状態を考慮することで異なる時刻で得られた各センサ組の TDOA/FDOA 観測値を用いて測位することが可能となる．そのため早期の追尾開始が期待できる．. 23.

(29) 第 3 章非同期測位方式. Proposed method Reverse extrapolation TDOA/FDOA memory. Residual computation. Gradient computation. TDOA/FDOA memory. Residual computation. Gradient computation. TDOA/FDOA memory. Residual computation. TDOA/ FDOA. Sensor3. ・・・. Gradient computation. ・・・. ・・・. ・・・. TDOA/ FDOA. ・・・. Sensor2. TDOA/ FDOA. ・・・. Sensor1. Least squares. 図 3.3. Estimated state vector. ASREX 処理ブロック図. また，ASREX の処理フローチャートを図 3.4 に示す．図 3.4 を用いて本方式の処理の流れを説明する．まず，センサ組毎に観測値を蓄積していき，式(3.17)で与えられるヤコビ行列 G がフルランクとなるまで時刻 tk を更新する（ブロック(2)(3)）．ここで，具体的には，目標の 3 次元位置及び速度を含む 6 次元ベクトルの解を得るためには，ヤコビ行列 G のランクが 6 となるまで更新することになる（この場合，ブロック(2)では Nrank=6 とする）．次に観測値が必要数揃った場合，現時刻を基準時刻 t として設定し，適当な初期値 x(t)init を設定する(ブロック(4),(5))．そして，式(3.14)より過去の時刻 tk に遡って状態ベクトル xk を求める(ブロック(6))．式(3.16)，(3.17)を用いて，状態ベクトルに応じた残差ベクトル δh とヤコビ行列 G を求め，さらに式(3.8)，(3.9)を用いて状態ベクトルを更新する(ブロック(7)～(10))．以上の処理を繰り返し，  x が事前に決めた閾値εより小さくなった場合あるいは，逐次最小 2 乗法の反復回数”itr”が事前に決めた閾値”Lim”に達した場合に，計算を中断し初期状態ベクトルの推定値を得る(ブロック(11),(12))．ここで，δx は初期状態ベクトルの補正量であり，. はユークリッドノルムを表す．. 24.

(30) 第 3 章非同期測位方式. (1)START. (3)Updated time tk k=k+1. (2)Full rank judgment rank(G)=Nrank. No. Yes (4)Setting reference time t t= tk (5)Initial state vector of x(t) x(t)=x(t)init (6)Reverse extrapolation (Eq.(3.14)) (13)Updated iteration itr=itr+1. (7)Jacobian matrix G (Eq.(3.17)) (8)Estimate error vector δh (Eq.(3.16)) (9)Correction vector δx (Eq.(3.8)) (10)Update state vector x(t) (Eq.(3.9)). (11)Correction vector No judgment No ||δx||<ε Yes. (12)Iteration judgment itr=Lim. No. Yes. (14)END. 図 3.4. ASREX のフローチャート. 次に，3 次元空間上における初期状態ベクトルの補正量  x の導出方法を具体的に示す． 3.3.1 項では TDOA 観測値しか観測されない場合，3.3.2 項では TDOA と FDOA 観測値が観測される場合について補正量を導出する. まず，時刻 tk におけるセンサ i の位置ベクトル si,k，目標の初期状態ベクトル推定値 xˆ (t ) ，時刻 tk における目標の状態ベクトル推定値 xˆ k は以下の式(3.19)-(3.21)になる． si ,k   six,k. siy ,k. siz ,k . T. (3.19) 25.

(31) 第 3 章非同期測位方式. xˆ (t )  uˆ (t )T . xˆ k  uˆkT. T. uˆ (t )T    xˆ(t ) yˆ (t ) zˆ(t ) xˆ(t ) yˆ (t ) zˆ(t )     T. uˆkT    xˆk. yˆk. zˆk. xˆk. yˆk. zˆ k . T. (3.20). T. (3.21). 目標が等速直進すると仮定すると，式(3.14)より ˆ  xˆk   xˆ (t )  x(t )  tk    uˆ k   yˆ k    yˆ (t )  yˆ (t )  tk    zˆ    k   zˆ (t )  zˆ(t )  tk . (3.22).  xˆk   xˆ (t )      uˆ k   yˆ k    yˆ (t )       zˆk   zˆ (t )     . (3.23). が成立する．. 3.3.1. TDOA 観測時の処理式詳細. まず，式(3.14)で定義される時刻 tk における目標の状態ベクトル推定値 xˆ k とセンサ i，センサ j との距離差 hij ( xˆ k ) は以下の式となる． hij ( xˆ k ) .  uˆ. k.  si ,k . T.  uˆ. k.  uˆ.  si ,k  . k.  s j ,k.   uˆ T. k.  s j ,k. . (3.24). 時刻 tk において観測される目標とセンサ i，センサ j との距離差を rij,k とすると，距離差の残差 Δhij,k は以下の式で定義される．. hij ,k  rij ,k  hij ( xˆ k ). (3.25). 例えばサンプリング時刻毎にセンサ一組の TDOA 観測値しか得られない状況では 6 サンプリング分の観測値を使用することになる.この場合，残差ベクトル δh は以下の式のようになる．.  h  h21,1 h31,2. h41,6 . T. (3.26). 次にヤコビ行列 G は. 26.

(32) 第 3 章非同期測位方式.  h21,1 xˆ1  ˆ  x1 xˆ (t )   G   h31,5 xˆ5  xˆ xˆ (t )  5  h41,6   xˆ (t ). h21,1 yˆ1 yˆ1 yˆ (t ). h21,1 zˆ1 zˆ1 zˆ(t ). h21,1 xˆ1 xˆ1 xˆ (t ). h21,1 yˆ1 yˆ1 yˆ (t ). h31,5 yˆ5 yˆ5 yˆ (t ). h31,5 zˆ5 zˆ5 zˆ(t ). h41,6 yˆ (t ). h41,6 zˆ(t ). h31,5 xˆ5 xˆ5 xˆ (t ) h41,6 xˆ (t ). h31,5 yˆ5 yˆ5 yˆ (t ) h41,6 yˆ (t ). h21,1 zˆ1   zˆ1 zˆ (t )    h31,5 zˆ5  zˆ5 zˆ (t )   h41,6   zˆ (t ) . (3.27). となる．以上より，式(3.8)に式(3.26)，(3.27)を代入して補正量  x を求める．ここで，式(3.27)の各要素は以下の式で表される． hij , k  xˆk hij , k  yˆ k. hij ,k zˆk. .  xˆ. k.  six , k .  uˆk  si,k . T.  yˆ. k.  uˆ. k.  siy , k.  si , k . T.  zˆ. k.  uˆk  si,k . .  uˆ. k.  si ,k .  siz , k .  uˆk  si,k   uˆk  si,k  T. . .  xˆ . uˆ k  s j , k.  yˆ. k.  uˆ. k. .  s jx, k. k. . uˆ k  s j ,k.  s jy , k.   uˆ.  s j ,k.  zˆ.  s jz , k. . T. k. uˆ k  s j ,k.  T. . .  s j ,k k. .  T. uˆ k  s j ,k. . . (3.28). (3.29). (3.30). xˆk yˆ k zˆk   1 xˆ (t ) yˆ (t ) zˆ(t ). (3.31). xˆk yˆ k zˆ   k  tk xˆ (t ) yˆ (t ) zˆ(t ). (3.32). xˆ (t ) yˆ (t ) zˆ(t )   0 xˆ (t ) yˆ (t ) zˆ(t ). (3.33). hij , k h h ， ij , k ， ij , k は式(3.28)-(3.30)の xˆk ， yˆ k ， zˆk ， uˆ k をそれぞれ xˆ (t ) ， yˆ (t ) ， zˆ(t ) ， xˆ (t ) yˆ (t ) zˆ (t ). uˆ (t ) に置き換えれば求まる．. 3.3.2. TDOA と FDOA 観測時の処理式詳細. 3.3.1 節の TDOA 観測値に加え，FDOA 観測値も考慮すると，式(3.14)で定義される目標の状態ベクトル推定値 xˆ k とセンサ i，センサ j とのドップラ速度差 hij ( xˆ k ) は以下の式となる．. 27.

(33) 第 3 章非同期測位方式. hij ( xˆ k ) .  uˆ.  si ,k.  uˆ.  si ,k   uˆ k  si ,k . k. k.   uˆ T. k.  si ,k . T. .  uˆ.  s j ,k.  uˆ.  s j ,k   uˆ k  s j ,k . k. k.   uˆ T. k.  s j ,k . (3.34). T. 時刻 tk において観測される目標とセンサ i，センサ j とのドップラ速度差を rij ,k とすると，ドップラ速度差の残差 hij ,k は以下の式で定義される．. hij ,k  rij ,k  hij ( xˆ k ). (3.35). TDOA と FDOA 観測値が一組で得られるため，例えばサンプリング時刻毎に TDOA と FDOA 観測値の組が一つしか得られない状況では 3 サンプリング分の観測値を使用することになる(図 3.2 参照)．この場合，残差ベクトル δh は以下の式となる．.  h  h21,1 h31,2 h41,3 h21,1 h31,2 h41,3 . T. (3.36). 次にヤコビ行列 G は.  h21,1 xˆ1  ˆ  x1 xˆ (t )  h31,2 xˆ 2  ˆ ˆ  x  x ( t)  2  h 41,3   xˆ (t ) G   h21,1 xˆ1   xˆ1 xˆ (t )   h31,2 xˆ2  xˆ2 xˆ (t )   h41,3   xˆ (t ). h21,1 yˆ1 yˆ1 yˆ (t ). h21,1 zˆ1 zˆ1 zˆ(t ). h21,1 xˆ1 xˆ1 xˆ (t ). h41,3 zˆ(t ). h41,3 xˆ (t ) xˆ (t ) xˆ (t ) h21,1 xˆ. h21,1 zˆ1 zˆ1 zˆ (t ). h41,3 zˆ (t ). 1. xˆ1 xˆ (t ). h41,3 xˆ (t ). h21,1 yˆ1 yˆ1 yˆ (t ). h21,1 zˆ1   zˆ1 zˆ (t )  h31,2 zˆ2   zˆ2 zˆ (t )  h41,3 zˆ(t )  zˆ (t ) zˆ (t )   h21,1 zˆ1   zˆ1 zˆ (t )   h31,2 zˆ2  zˆ2 zˆ (t )   h41,3   zˆ (t ) . (3.37). となる．以上より，式(3.8)に式(3.36)，(3.37)を代入して補正量  x を求める．ここで，式(3.37)の各要素は以下の式で表される．ここで，Rk,i ，Rk,j を以下の式(3.38)， (3.39)のように定義した．. Rk ,i  (uˆk  si ,k )T (uˆk  si ,k ). (3.38). Rk , j  (uˆk  s j ,k )T (uˆk  s j ,k ). (3.39). 28.

(34) 第 3 章非同期測位方式.  (uˆ  s )T (uˆ  s )  (uˆk  si ,k )T (uˆk  si ,k )  k  i ,k k i ,k ˆ ˆ ( x  s )  ( x  s )  k ix , k k ix , k  T T (uˆk  si ,k ) (uˆk  si ,k )  (uˆ k  si ,k ) (uˆ k  si ,k )  T   (uˆ k  s j , k )T (uˆ k  s j ,k ) 1  (uˆ k  s j , k ) (uˆ k  s j ,k ) ˆ  ˆ  ( x  s )  ( x  s )  k jx , k k jx , k  T Rk , j 2  (uˆ k  s j , k )T (uˆ k  s j ,k ) (uˆ k  s j , k ) (uˆ k  s j ,k )  . (3.40).  (uˆ  s )T (uˆ  s )  (uˆk  si ,k )T (uˆk  si ,k )  k  i,k k i ,k ˆ ˆ ( y  s )  ( y  s )  k iy , k k iy , k  T T (uˆk  si ,k ) (uˆk  si ,k )  (uˆ k  si ,k ) (uˆ k  si ,k )  T   (uˆ k  s j , k )T (uˆ k  s j ,k ) 1  (uˆ k  s j , k ) (uˆ k  s j ,k ) ˆ  ˆ  ( y  s )  ( y  s )  k jy , k k jy , k  2 T T Rk , j  (uˆ k  s j , k ) (uˆ k  s j ,k ) (uˆ k  s j , k ) (uˆ k  s j ,k )  . (3.41).  (uˆ  s )T (uˆ  s )  (uˆk  si ,k )T (uˆk  si ,k )  k  i ,k k i ,k ( zˆk  siz , k )  ( zˆk  siz ,k )   T T (uˆk  si ,k ) (uˆk  si ,k )  (uˆ k  si ,k ) (uˆ k  si ,k )  T   (uˆ k  s j , k )T (uˆ k  s j ,k ) 1  (uˆ k  s j , k ) (uˆ k  s j ,k ) ˆ  ˆ  ( z  s )  ( z  s )   k jz , k k jz , k 2 T T Rk , j  (uˆ k  s j , k ) (uˆ k  s j ,k ) ˆ k  s j , k ) (uˆ k  s j ,k ) ( u  . (3.42). hij , k 1  ˆ xk Rk ,i 2. hij , k 1  yˆ k Rk ,i 2. hij , k zˆk. . hij ,k  xˆ k. hij ,k  yˆ k. hij ,k  zˆ k. 1 Rk ,i 2.  xˆ. k.  six ,k .  uˆk  si,k   uˆk  si,k  T.  yˆ. k.  uˆ. k.  siy , k. .  si ,k   uˆ k  si ,k  T.  zˆ. k.  siz , k .  uˆk  si,k   uˆk  si,k  T. . .  xˆ. k. . uˆ k  s j ,k.  yˆ. k.  uˆ. k. . k.  T. . uˆ k  s j ,k.   uˆ T. k.  s jz , k. uˆ k  s j ,k.  T.  s j ,k. . (3.44). . (3.45). . uˆ k  s j ,k. xˆk yˆ k zˆ   k 1 xˆ (t ) yˆ (t ) zˆ(t ). また，. . (3.43). .  s jy , k.  s j ,k.  zˆ .  s jx, k. (3.46). hij , k h h ， ij , k ， ij , k は式(3.40)-(3.42)の xˆk ， yˆ k ， zˆk ， uˆ k をそれぞれ xˆ (t ) ， yˆ (t ) ， xˆ (t ) yˆ (t ) zˆ (t ). zˆ(t ) ， uˆ (t ) に置き換えれば求まる．さらに. hij , k hij , k hij , k ，，についても同様に式(3.43)zˆ (t ) xˆ (t ) yˆ (t ). 29.

(35) 第 3 章非同期測位方式. (3.45)の xˆk ， yˆ k ， zˆk ， uˆ k をそれぞれ xˆ (t ) ， yˆ (t ) ， zˆ(t ) ， uˆ (t ) に置き換えれば求まる．他の要素については，3.3.1 と同様である（式(3.28)-(3.33)参照）．. 3.4. 計算機シミュレーションによる性能評価. 本章では追尾開始時における目標の初期状態ベクトルを推定する方法として従来測位方式と 3.3 節で定式化した提案方式 ASREX の性能比較を行う．. 3.4.1. 3.4.1.1. シミュレーション条件. シナリオ. シミュレーションシナリオを図 3.5 に示す．受信局（受信専用センサ）の数を 4 つ，送信局（送信専用レーダ）の数を 1 つとした．受信局と送信局の位置は既知であり，受信局 1 を. [0km,-20km,0km], 受信局 2 を [-40km,0km,0km], 受信局 3 を [40km,0km,0km], 受信局 4 を [0km,0km,10km] 送信局を[0km,40km,20km]に設置した．図 3.5 において，四角い印がレーダとセンサを表す．観測に用いる受信局の組合せは 1 と 2，1 と 3，1 と 4 とした．ここで， 3.2 節で述べたように他の受信局の組み合わせは上記の組み合わせで表現されるため，上記 3 つの TDOA が同時に観測されない場合は測位されないシナリオとした．目標は[-20km,70km,10km]を出発して 300m/s で 200 秒間等速直進し，目標の終点を[20km,10km,10km]とした．シナリオとして，送信局が電波を目標に照射して，その目標からの反射波を各受信局が受信する場合を想定した．ここで，本シナリオでは図 3.5 のようにセンサを分散配置させている．これは， TDOA/FDOA を用いた測位ではセンサが目標から見て一列に並んでいる場合には逐次最小 2 乗法が解けない，もしくは精度が劣化する可能性があるからである．. 30.

(36) 第 3 章非同期測位方式. 80. Target (Altitude:10km). 60 Y[km]. 40 20 Sensor2 (Altitude:0km). -50. Radar (Altitude:20km). -30. 0 -10 -20. Sensor4 (Altitude:10km). 10. Sensor3 (Altitude:0km). 30. 50. Sensor1(Altitude:0km). X[km] 図 3.5. 3.4.1.2. シミュレーションシナリオ. 観測条件. 本章では，センサを移動させることも想定している．移動させる場合は，センサの時刻同期誤差とセンサの位置計測誤差が生じる．そのため，時刻同期誤差および位置計測誤差が発生しない理想状態を観測条件Ⅰ，時刻同期誤差および位置計測誤差が発生する状態を観測条件Ⅱとして 2 ケースのシミュレーションを実施した．また，観測誤差及びサンプリング間隔，データ検出率，センサ時刻同期誤差および位置計測誤差についての観測条件を以下に述べる．. （１）観測誤差及びサンプリング間隔距離差観測値に加わる誤差の標準偏差を 18.7m，ドップラー速度差観測値に加わる誤差の標準偏差を 0.187m/s とした．ここで，上記の観測誤差は 4.6.2 節で後述する受信 S/N 理論値より設定した．また，サンプリング間隔は 2s かつ一定とした．観測条件Ⅰ，Ⅱともに以上の条件を用いる．. 31.

(37) 第 3 章非同期測位方式. （２）データ検出率マルチパスや電波の遮蔽等の影響によりデータが欠落する場合を簡易的に模擬するために，各センサで信号が検出される確率を表すパラメータとして，データ検出率 Pd を定義した．データ検出率については実際に空港内で測定した結果，最も低くて 0.62 程度まで検出率が落ちる事例が報告されている[19]．そこで，更に最悪となるケースを想定し，最低検出率を 0.5 として設定し，0.62 より高い設定値として 0.75（0.5 と 1 の中間値），欠落なしの設定値として 1 の計 3 通りのデータ検出率について検証する．なお，観測条件Ⅰ，Ⅱともに以上の条件を用いることとする．. （３）センサ時刻同期誤差および位置計測誤差. センサには GPS(Global Positioning System)が搭載されているとする．文献[28][30]によれば， GPS を用いるとセンサの時刻同期誤差は数 ns, 位置計測誤差の標準偏差は各軸 1m であるとされている．ここで，センサの時計がずれていることによるセンサ間の時刻同期誤差を数 ns とすると TDOA 観測値にそのまま数 ns 分バイアス誤差が含まれることになる．本章では上記誤差の影響を確かめるために，文献[28][30]より大きい誤差を与えることにした．つまり，基準センサ 1 に対してセンサ 2,3,4 の時刻同期誤差を 60ns（3.4.1.2（１）の距離差観測誤差 18.7m 相当），センサ位置計測誤差の標準偏差を各軸 2m とした．観測条件Ⅱのみ，以上の条件を用いる．. 3.4.1.3. パラメータ. 各種パラメータ設定値を以下に示す．（１）観測精度式(3.10)の観測誤差共分散行列 R は，式(3.47)のように非対角成分の値を 0 とし，対角成分もセンサ i,j の組合せに依らず一定とした．  2  R  0 . 0     2 . (3.47). 表 3.1 より，TDOA を利用する場合は式(3.47)の σ として距離差精度 σTDOA，FDOA を利用する場合はドップラ速度差精度 σFDOA を参照した（設定した観測誤差の標準偏差と同じ値とした）． 32.