GPS 搭載機の横方向経路維持誤差の分布モデル
電子航法開発部 ※天井 治、長岡 栄
１．まえがき従来から航空機は洋上空域において、慣性航法装置（Inertial Navigation System：INS）や慣性基準装置（Inertial Reference System：IRS）と飛行管理システム（Flight Management System：FMS）とを結合した航法装置の情報により飛行している。INS や IRS では、加速度の測定値を 2回積分して位置を求めるため、測位誤差は飛行時間と共に増加する。米国連邦航空局の規格では、INS の 95％の飛行誤差は、10 時間以内の飛行では１時間あたり 2 NM（Nautical Mile：海里）以下である(1)。
近年、ボーイング 747-400（B747-400）、B777
エアバス 340-300 など GPS（Global PositioningSystem：全地球的測位システム）を搭載した航空機（GPS 搭載機と呼ぶ）が飛行している。GPSは人工衛星を用いた測位システムで、現在カーナビゲーションや測量などの分野で広く利用されている。GPS 搭載機では、機上の GPS センサで得た位置情報で自機の位置を更新できる。このため、GPS 搭載機は GPS 非搭載機よりも高い航法精度での飛行が期待できる。GPS 搭載機が実際にどの程度の航法精度で飛行しているかは、文献(2)、(3)等で発表されている。しかし、いずれも標本数が少なく、分布モデルを推定することはできなかった。
本稿では、特に横方向の航法精度に着目し、北太平洋の洋上航空路を飛行した GPS 搭載機の横方向経路逸脱量を調べ、最も良く適合する分布モデルを推定した方法と結果を示す。
２．GPS 搭載機(4)本稿でとりあげる GPS 搭載機は、FANS-1 もしくは FANS-A パッケージを搭載した航空機である。FANS-1 は、国際民間航空機関の将来航空航法システム（Future Air Navigation System：FANS）構想を実現するためにボーイング社が開発したアビオニクス・パッケージで、航法、通信 、 監 視 機 能 と し て それ ぞ れ GPS 、CPDLC(Controller Pilot Data Link Communications) 、ADS(Automatic Dependent Surveillance)の機能を備えている。
FANS-1 は 1995 年の 6 月以来、米国連邦航空局などの承認を得ており、現在 B747-400、B777などに搭載されている。また、FANS-1 を搭載した B747-400 は、次の航法機器およびそのセンサを備えている。
1) 飛行管理計算機（FMC）
２基
2) 慣性基準装置（IRS）
３基
3) GPS 受信機
２基
FANS-A はエアバス社の FANS-1 と同様の装置であり、A330 や A340 に搭載されている。
３．横方向経路逸脱量
航空機の進行方向、横方向および高さ方向に対しそれぞれ航法精度が考えられるが、ここでは横方向の航法精度のみに着目する。
航空機の運航時には、予め航空機が飛行すべき経路（割当経路）が航空交通管制機関により割り当てられる。本稿では、割当経路の中心線上を飛行しようとしている航空機が、航法誤差などのために経路の中心線から逸脱した量（横方向経路逸脱量と呼ぶ）を調べる。
本稿では、地球を回転楕円体と見なして横方向経路逸脱量を計算する。図１のように経路を構成する地点 FIX1（緯度ϕ1, 経度λ1, 高さ 0）、FIX2（ϕ2, λ2, 0）、FIX3（ϕ3, λ3, 0）を考える。
ここでϕ3>ϕ1>ϕ2 とする。回転楕円体上の２点間の距離は厳密には測地線となるが、その距離が 100 km より小さいときには法線面と測地線

図１ 横方向経路逸脱量 d の幾何学的関係

との開きの差は通常無視できる(5)。このことを利用し、
FIX1 における法線ベクトルを含み FIX2または FIX3 を通る平面（法線面）を経路の中心
線 と 見 な す 。 こ こ で 地 心 (Earth-centered,Earth-fixed (ECEF)）座標系では、その法線ベク
ト ル は n1(nX1, nY1, nZ1) と 表 せ る 。 同 様 に 、FIXi(i=1,2,3)は(Xi, Yi, Zi)と表せる。
横方向経路逸脱量 d は航空機の位置 P（ϕaa, λa,ha）から上記の法線面に下した垂線の長さとして定義する。d の符号は、経路の北側（法線面に対し北極を含む側）を正とした。例えば
d=+2.0NM は、割当経路の中心線から北側に2.0NM ずれて飛行したことを意味する。
横方向経路逸脱量 d は、レーダデータを用いて、以下の計算式により求めた。

d =C 1 X a + C 2Y a + C 3 Z a + C 4
C 12 + C 2 2 + C 3 2(1)

ここでX a = ( N (ϕ a ) + h a ) cos ϕ a cos λ a

Ya = ( N (ϕa ) + ha ) cos ϕa sin λa
Z a = ( N (ϕ a )(1 − ε 2 ) + h a ) sin ϕ a(2)(3)(4)

ϕ 3 > ϕ a ≥ ϕ1 のとき
C 1 = (Y 3 − Y 1) n Z 1 − ( Z 3 − Z 1) n Y 1
C 2 = ( Z 3 − Z 1) nX 1 − ( X 3 − X 1) n Z 1
C 3 = ( X 3 − X 1) n Y 1 − (Y 3 − Y 1) n X 1(5)(6)(7)

C 4 = (Y 3Z 1 − Y 1Z 3)nX 1 + ( Z 3 X 1 − Z 1 X 3)nY 1+ ( X 3Y 1 − X 1Y 3)nZ 1(8)ϕ1 > ϕ a ≥ ϕ 2 のとき,C 1 = (Y 1 − Y 2 ) n Z 1 − ( Z 1 − Z 2 ) n Y 1C 2 = ( Z 1 − Z 2 ) nX 1 − ( X 1 − X 2 ) nZ 1C 3 = ( X 1 − X 2 ) n Y 1 − (Y 1 − Y 2 ) n X 1
C4 = (Y 1Z 2 − Y 2Z1)nX 1 + (Z1X 2 − Z 2 X 1)nY 1+ ( X 1Y 2 − X 2Y 1)nZ1

(9)
(10)
(11)
(12)

ここでnX 1 = f (ϕ 1) cos 2 ϕ 1 cosλ 1
n Y 1 = f (ϕ 1) cos 2 ϕ 1 s in λ 1
n Z 1 = f (ϕ 1) sin ϕ 1 cos ϕ 1(13)X i = N (ϕ i ) cos ϕ i cos λ i
Y i = N (ϕ i ) cos ϕ i sin λ i (i = 1,.., 3)
Z i = N (ϕ i )(1 − ε 2 ) sin ϕ i(16)である。但し

(14)
(15)

(17)
(18)

f (ϕ ) = − N (ϕ ) 2N (ϕ ) =1− ε 2
1 − ε 2 sin 2 ϕa(19)(20)

1 − ε 2 sin 2 ϕとした。
ここで、a は回転楕円体の長半径（赤道半径）を、ε2 は離心率を表す。離心率は、扁平率 f を用いて、ε2 = (2 - f)f と表すことができる。実際の計算では測地系として WGS-84 測地系を用い、
a = 3443.918467 NM (1 NM = 1852 m)および ε2 =6.69437999x10-3 とした。
(1)式は航空機と法線面の距離を求める式、(2)〜 (4)式は航空機の位置の座標変換式、 (5)〜 (8)および(9)〜(12)は法線面の方程式の係数である。
(2)〜(4)式で用いた航空機の位置 Pa（ϕa, λa, ha）
は、レーダによる測位データ（レーダアンテナからの距離ｒ，レーダ回転角θ，飛行（気圧）高度 h）を用いて文献(2)に示す方法により計算した。
４．使用データ
使用データとデータと観測対象
4.1 使用データ
横方向経路逸脱量を求めるにあたり、次の３種類のデータを用いた。
(1) 釧 路 航 空 路 監 視 レ ー ダ （ Air Route
Surveillance Radar: ARSR）のデータ(2) 飛行計画情報(3) ADS レポート
ARSR に は 、 二 次 監 視 レ ー ダ （ Secondary
Surveillance Radar：SSR）が付いている。SSR では、航空機に搭載されている自動応答装置（トランスポンダ）に対して質問信号を送り、それに対しトランスポンダが応答信号を返すことにより、航空機の位置座標（レーダアンテナからの距離ｒ，レーダ回転角θ）を測定する。また、応答信号から（航空機の気圧高度計で測った）飛行高度（Mode C コード）h および航空機の識別符号（Mode A コード）の情報を得ることができる。
本稿の計算では SSR データのみを計算に用いた。表１に釧路 ARSR に付いている SSR の諸元を示した。得られる航空機の位置情報は ARSRの回転周期に対応し約 10 秒毎の離散データとなる。

表１ 釧路 ARSR の SSR 部の諸元
項目

仕様

アンテナタイプ
アンテナ回転周期
最大覆域
周波数(送信)(受信)
送信先頭出力
送信繰り返し周波数

コフィード型
約 10 秒250 NM
1030 MHz
1090 MHz
1.5 kW
285 pps
測距系: 0.125 NM
測角系: 0.088°

最小量子化単位

図２ 観測空域
飛行計画情報には、便名、機種名、識別符号、飛行経路、飛行経路上の各位置通報点の通過時
刻および飛行高度、機体登録番号など飛行に関わる様々な情報が記載されている。位置通報点（FIX と呼ぶ）は、航空機がその上空を通過する際、管制機関に対して位置通報（通過時刻、飛行高度などの報告）を行うために定められた地点である。レーダデータからは、便名、機種名、飛行経路、機体登録番号などの情報を得ることができないため、飛行計画情報を使用し、それらを得た。
北太平洋航空路では ADS の評価運用が行われており、航空機の何割かは ADS レポート(6)を航空交通管制官に対し送信する。本稿では、ADSレポートを送信した航空機は GPS 搭載機と見なした。これは、ほとんどの航空機が FANS-1 もしくは FANS-A を搭載して GPS センサと同時にADS 機能を装備しているからである。ADS レポートは東京航空交通管制部の DLP （Data LinkProcessing）システムから抽出した。
4.2 観測対象
次の条件を満たす航空機を観測対象機とした。
(i)北太平洋航空路 R220 の位置通報点 NANNO
上を飛行(ii)GPS を搭載(i)は、飛行計画情報に記載された飛行経路の情報および FIX 情報を調べた。(ii)については、ADS レポートの有無および機体登録番号から調べた。図２に観測空域を示す。北太平洋航空路
は、日本と北米を結ぶ北太平洋上の航空路で、5本の航空路がある。R220 はこの内、最も北側に位置する航空路であり、西行き便（一日に約 70便(7)）のみに使用されている。釧路 ARSR のデータからは、R220 を飛行する航空機の一部の空域におけるレーダ情報が得られる。図中の実線の円は釧路 ARSR の最大覆域（250 NM）を示し、破線は八戸 ARSR の最大覆域（200 NM）を示す。
５．解析方法5.1 GPS 搭載機に対する航跡データの作成
解析にあたり、まず ARSR による測位データを飛行高度と時刻により分類し、データ中のすべての航空機の飛行便毎の航跡データを作成した。
次に、航空機を特定するために航跡データと飛行計画情報との照合を行った。この照合には、双方に含まれる航空機の識別符号の情報を利用した。照合できたものに対して、割当経路や機種名、機体登録番号の情報を航跡データに付加した。
更に、ADS レポートに記載されている機体登録番号のみを GPS 搭載機として抽出した。

5.2 横方向経路逸脱量の計算
航跡データは複数個の約 10 秒毎の離散データから構成される。本稿では、NANNO から経路に沿って測った距離が±1.5 NM 以内のデータを抽出し、その中で NANNO に最も近いデータを横方向経路逸脱量の計算に使用した。横方向経路逸脱量は、３章に示した方法で計算した。
計算に当たっては、レーダ自身のランダム誤差の影響を軽減させるため、航跡データの平滑化処理(8)を行った。

６．結果
2001 年 12 月 1 日から 2002 年 5 月 23 日までの 174 日分のデータを用いて解析を行った。レーダの航跡データを精査し、計画的オフセットのような経路の中心線からの故意の逸脱と見られるものは除外した。

6.1 横方向経路逸脱量の分布
図３に横方向経路逸脱量の相対度数分布を示す。横軸は横方向経路逸脱量を、縦軸は対数目盛の相対度数を表す。標本数は 3,150 便（248 機体）、標準偏差は 0.0259 NM であった。横方向経路逸脱量の 95%値は[-0.045 NM, +0.051 NM]であった。
図４には、GPS 非搭載機の横方向逸脱量の分(9)布 を示す。観測空域は GPS 搭載機と全く同じであるが観測期間は 1993 年 5 月 1 日〜1994 年 4月 30 日である。標本数は 17,702 便、標準偏差は 1.17 NM である。
図３と図４を比較すると、GPS 搭載機の横方向経路逸脱量の標準偏差は、GPS 非搭載機のそれの 45 分の１であることが分かる。
6.2 機種の割合
表２に GPS 搭載機の各機種の割合を示した。
B747-400 は全体の約 76%を占める。また、GPSを搭載した B747-200 が２便観測された。
6.3 分布のモデル化
実測分布に最も良く適合する確率密度関数を知ることは、管制間隔の短縮可能性を調べる上で重要である。ここでは、GPS 搭載機の横方向経路逸脱量の分布に最も良く合う確率密度関数を推定する。
6.3.1 当てはめた分布モデル
表３に当てはめに用いた分布モデルを示す。
N は正規（Normal）分布、DE は両側指数（DoubleExponential ）分布である。 N-N 、 N-DE 、 DDEはそれぞれ、２つの正規分布の混合型分布、正規分布と両側指数分布の混合型分布、２つの両側指数分布の混合型分布である。

図３ GPS 搭載機の横方向経路逸脱量の分布

図４ GPS 非搭載機の横方向経路逸脱量の分布（N は標本数、std.dev.は標準偏差を表す）
表２ GPS 搭載機の機種の割合
機種

便数割合

B747-400

2,38075.5 %B777-20063920.3 %A340-300

1294.1 %B747-20020.1 %3,150100 %総計

表３ 当てはめたモデルの確率密度関数
モデルNDE
N-N
N-DEDDE

確率密度関数exp(-y /(2σ2))/(2πσ2)1/2≡N(y|σ)
exp(-|y|/λ)/(2λ)≡DE(y|λ)(1-α)N(y|σ1)+αN(y|σ2)0<α<1(1-α)N(y|σ)+αDE(y|λ)
0<α<1(1-α)DE(y|λ1)+αDE(y|λ2)0<α<1, 0<λ1<λ22

6.3.2 当てはめ方法
各分布のモデルパラメータは最尤推定法により推定した。この方法では、当てはめの良さを示す尤度関数を最大化するパラメータを推定する。これらは、パラメータが解析的に計算できる N と DE 以外は、Newton-Raphson 法により求めた。Newton-Raphson 法の初期値は、区間分割法と呼んでいる次の方法で求めた。
1) 予め決めたパラメータの範囲をいくつかの区間に分割
2) それらの区間の中で尤度が最大となる区間を選択
3) 1)、2)を区間が十分小さくなるまで繰り返す
図３に示した横方向経路逸脱量の分布は中央値で 0.005 NM 偏りがあったが、各データから 0.005 を引いて中央値を 0 にした分布に対し、当てはめを行った。
6.3.3 当てはめの結果
表４に各モデルのパラメータの推定値を示す。ここで、AIC は赤池情報量基準(10)でモデルの当てはめの良さを表す。これは次式で定義される。
AIC=-2L + 2k(21)
ここで、L はモデルの対数最大尤度で、k は推定するモデルパラメータの数である。値が小さい程良い当てはめであることを示す。AIC を使うことにより、モデルパラメータ数によらずに当てはめの良さを比較することができる。
Dn は Kolmogorov-Smirnov 統計量で、次式で定義される。

Dn = sup Fn ( x ) − F ( x )(22)x

ここで、Fn(x) と F(x) はそれぞれ、実測累積分布と当てはめた累積分布である。Dn も値が小さいほど良い当てはめであることを示す。図５には、実測分布と当てはめた分布を示す。横軸は、横方向経路逸脱量の絶対値を、縦軸は対数目盛の（1-累積相対度数）を示す。
表４の AIC と Dn の値を見ると、N-DE モデ
ルが実測分布に最も良く当てはまる分布であることが分かる。パラメータの推定値は、
α=0.0564, σ=0.0232 NM, λ=0.0380 NM であった。

表４ モデルパラメータの推定値
推定値
AICDnモデル-14,087N
σ = 0.0259 NM0.047(229)-14,104DE
λ = 0.0196 NM
0.052(212)α = 0.0239-14,316N-N
σ1 = 0.0234 NM
0.029(0)σ2 = 0.0749 NM
α = 0.0564-14,316N-DE
σ = 0.0232 NM
0.028
(0)
λ = 0.0380 NM
α = 6.18x10-8-14,099
DDE
σ1 = 0.0196 NM
0.052
(217)
σ2 = 0.0202 NM

図５ モデル当てはめの結果（N は標本数を表す）
７．考察
7.1 GPS 搭載機の横方向誤差の標準偏差６章で示した横方向逸脱量の観測データには、測位に使用したレーダ自身の誤差も含まれている。レーダの測位精度は、測角誤差、測距誤差、データ取得時刻の誤差、データの量子化誤差、電波伝搬誤差などに依存する。量子化誤差と測位誤差が主な要因と考える。釧路 SSR の測角系／測距系および飛行高度のデータの最小量子化単位は、それぞれ 0.088°、0.125 NM、100 ft である。しかしながら、これらの誤差の影響は航跡データの平滑化により軽減されていると考える。
ここで、観測された GPS 搭載機の横方向経路逸脱量の標準偏差 σob が次式で計算できると仮定する。

σ ob = σ GPS + FTE 2 + σ m 2(23)ここで、σGPS+FTE は GPS 測位誤差と飛行技術誤差（Flight Technical Error：FTE）を考慮した GPS搭載機の横方向誤差の標準偏差を、 σm はレーダの測位誤差の標準偏差を表す。σGPS+FTE の推定値は 次式で計算できる。

σˆ GPS + FTE = σ ob 2 − σ m 2(24)一方、σGPS+FTE は次式で計算できる。

σ GPS + FTE = σ GPS 2 + σ FTE 2(25)ここで、σGPS は GPS 測位誤差の標準偏差、σFTEは飛行技術誤差 FTE の標準偏差である。 σGPSの推定値は 0.004 NM 程で、 σFTE の推定値はB747-400 の自動操縦（auto-coupled flight）の場合で 0.0625 NM(11)である。これらを(25)式に代入するとσGPS+FTE=0.0626 NM を得る。
ここでσob の推定値は 0.0259 NM である。σmは未知だが、σ̂ GPS + FTE はσob より大きくなることはない。従って、σ̂ GPS + FTE はσGPS+FTE よりも小さい。これは、σFTE の推定値 0.0625 NM が大きめの見積もりであるかもしれないことを示唆する。
８．まとめ
北太平洋航空路の R220 を飛行する GPS 搭載機に対する横方向経路逸脱量を計算し、分布モデルの当てはめを行った。2001 年 12 月 1 日から 2002 年 5 月 23 日までの 174 日分の釧路航空路監視レーダのデータ、飛行計画情報およびADS レポートを解析に使用した。次の結果が得られた。
(1) 観測された GPS 搭載機の横方向経路逸脱量の標本数 3,150 便（248 機体）であり、標準偏差は、0.0259 NM であった。
(2) GPS 搭載機の標準偏差は、同一観測空域のGPS 非搭載機の標準偏差 1.17 NM の 45 分の１である。
(3) GPS 搭載機の横方向経路逸脱量の 95%値は[-0.045 NM, +0.051 NM]であった。
(4) GPS 搭載機の横方向経路逸脱量の実測分布に最も良く当てはまる分布モデルは、正規分布と両側指数分布の混合型の分布（N-DE）であった。モデルパラメータの推定値は、
α=0.0564, σcore=0.0232 NM, λtail=0.0380 NMであった。
今後の課題には、(1) 釧路 ARSR のレーダの測位誤差の推定(2) Bootstrap 法などによるモデルパラメータの区間推定などがある。
［謝 辞］
データ収集の際にご協力頂いた国土交通省
札幌航空交通管制部、東京航空交通管制部および航空局の関係各位に感謝いたします。
[参 考 文 献](1) 加藤 昭英：航法システムの変遷と今後の展
望(1)、航空技術、350、p.8、1984(2) 天井 治、長岡 栄：洋上航空路における GPS装備機の航法精度の解析、日本航海学会論文集、102 号、pp.27-35、March、2000(3) G. Nezer : Empirical Data on Lateral pathkeeping Accuracy of GPS-equipped Aircraft,
ICAO RGCSP-WG/A-WP/18, May, 1999
(4) 長岡 栄：FANS-1 装備機と管制間隔短縮の動向、日本航海学会誌｢NAVIGATION」、130号、pp.1-8、Dec.、1996(5) 檀原 他：「測地・地球物理」､ 共立出版、
1984(6) ARINC Characteristic 745-2 ： AutomaticDependent Surveillance(ADS), June, 1993
(7) 天井 治、長岡 栄：北太平洋航空路の安全性の評価 −ルート構成変更後の近接通過頻度の計算−、電子情報通信学会技術研究
報告、SSS99-25、Oct.、1999(8) 天井 治、長岡 栄：航空管制用レーダの測角／測距誤差分布の解析、電子情報通信学
会技術研究報告、SANE95-60、April、1989(9) 天井 治、長岡 栄：北太平洋ルートにおける航空機対の水平重畳確率の推定、日本航海学会論文集、96、pp.11-19、March、1997(10) H. Akaike ： Information Theory and anExtension of the Maximum Likelihood
Principle, 2nd International Symposium on
Information Theory, Budapest, 1973(11) T. McDavid：Distribution of Lateral Navigation
Errors for FANS-1 Equipped 747-400 aircraft:
Implications for Reduced Separation, ICAO
RGCSP WG/A-WP/20, Brussels, May, 1995