ADS 縦方向予測誤差分布の特徴とモデル化
電子航法開発部 ※河道貴宏*

長岡栄

天井

治

*前 実 習 生 （ 日 本 大 学 大 学 院 ， 現 ㈱ ナ ナ オ ）1. まえがき
我が国では，自動従属監視（Automatic DependentSurveillance: ADS）を用いた距離間隔基準に基づく洋上管制方式の導入が検討されている[1]．この導入により航空交通容量の増加が期待されている．現在は評価運用を行っている．本格運用を行うに当たり，ADS を用いた管制方式の安全性評価が重要である．
洋上航空路の安全性評価では，Reich の衝突危険度モデル[2]が衝突危険度（単位時間内に発生する空中衝突事故件数の期待値）の推定にしばしば用いられる．Anderson はこのモデルを発展させ，ADS 環境下の衝突危険度モデルを提案した[3]．このモデルでは，衝突危険度を計算するためのパラメータの一つに外挿速度誤差がある．ところが，過去に実データから外挿速度誤差を計算した例は見られない．本研究では，実データを利用し，この外挿速度誤差に相当する ADS による現在位置と予測位置の差異（予測誤差）の分布を調査した．本稿では，進行方向（縦方向）に注目して ADS 予測誤差の計算方法とその分布を示す．次に分布の特徴を示す．さらに衝突危険度の計算で取扱い易いように予測誤差分布の数学的モデルを提案する．そのモデル化の方法と結果についても述べる．

2. 使用データ
北太平洋航空路（NOPAC）の洋上航空路を飛行する航空機から送られた周期レポート（現在位置，次の通報予定時の予測位置など）をもとに ADS 縦方向予測誤差分布を求めた．調査期間は 2000 年 4 月 6日から 2003 年 9 月 3 日である．ADS レポートの通報周期はすべて 16 分に設定されていた．

3. 予測誤差の定義と計算方法
予測誤差は一つ前の通報による予測位置と現在の通報位置との差異とする．いま，ある航空機の i 番目の通報時刻を ti，位置を Pi とする．ここで，Piは地球の中心を原点とする位置ベクトルとする．時刻 ti における次の通報の予測時刻を ti+1,p，位置をPi+1,p とする．しかし， i + 1 番目の通報時刻 ti+1 は予測時刻 ti+1,p とは数秒異なる場合が少なくない．このため，予測誤差を求めるに際し，予測時刻 ti+1,pにおける推定位置 Pi+1,e を次式により求め，これを予測位置に対応する通報位置とした．
Pi+1,e = Pi+1 + vi∆t
ただし
vi+1 = (Pi+1 – Pi) / (ti+1 – ti)

(1)
(2)

(3)
∆t = ti+1 – ti+1,pである．なお，vi+1 は速度ベクトルの推定値である．
予測誤差ベクトルを次式で定義する．
∆P ≡ Pi+1,p – Pi+1,e(4)このベクトルは，座標変換により，Pi+1,e を原点，進行（縦）方向を X 軸，それに直交する軸を Y 軸とする座標系で表せる． ∆P は次式で与えられる．
(5)∆P ≡ (∆PX, ∆PY, ∆Ph)T
縦方向予測誤差∆PX および横方向予測誤差∆PY の位置関係を図 1 に示す．
Pi-1 (t = ti-1)

Pi (t = ti)vi+1∆t∆PPi+1 (t = ti+1)
Pi+1,p(t = ti+1,p)

Pi+1,e (t = ti+1,p)∆PX∆PY
図 1 各位置の関係

4. ADS 縦方向予測誤差分布
4.1 全体の分布
本稿では，縦方向衝突危険度と関る ADS 縦方向予測誤差（∆PX）のみを取り扱う．図 2 に∆PX の分布および基本統計量を示す．∆PX の区間幅は 1NM で，[∆PX-0.5, ∆PX+0.5)NM のように設定した．縦軸は相対度数を対数で表示した．

1.0E+00

1.0E+00

Log10(Relative Frequency)

1.0E-01
1.0E-02

Mean=-1.339NMMedian=-1.340NM
SD=2.725NM
Skewness=0.639
Kurtosis=3.246

Log10(Relative Frequency)

Sample size=66336
Min.=-27.42 NM
Max=27.98 NM

1.0E-03
1.0E-04
1.0E-05

Sample size=12289
Min.=-8.662NM
Max=12.30NM
Median=0.163NM
Mean=0.334NM
SD=1.826NM
Skewness=0.766
Kurtosis=2.623

1.0E-01
1.0E-02
1.0E-03

B777-200

1.0E-04

1.0E-06-30-20-1001020

1.0E-0530

Prediction Error ∆PX [NM]-30-20

図 2 ADS 縦方向予測誤差分布-1001020
Prediction Error ∆PX[NM]30

図 5 B777-200 の分布

Log10(Relative Frequency)

1.0E+00
1.0E-01
1.0E-02
1.0E-03

Sample size=1480
Min.=-7.441NM
Max=27.98NM
Median=-1.340NM
Mean=-1.339NM
SD=5.475NM
Skewness=0.766
Kurtosis=0.375

A340-300

4.3 運航者（社）別の分布
同一機種（B747-400 型機）において運航者（航空会社）別に分類した結果，運航者毎でそれぞれ異なる分布形状を得た．本稿では全 16 社のうち，例として C 社と E 社の分布をそれぞれ図 6，図 7 に示す．
1.0E+00
Log10(Relative Frequency)

4.2 機種別の分布
ADS の情報は機上装置に依存している．そこで，機種別の∆PX 分布を求めた．機種の判別に，東京航空交通管制部の飛行計画情報処理システム（FlightData Processing system: FDP）から得た飛行計画情報を利用した．A340-300，B747-400，B777-200 の分布をそれぞれ図 3 から図 5 に示す．その他の機種は標本数が少ないため割愛する．その結果，3 種の機種でそれぞれ異なる形状の分布を得た．

1.0E-01

C社

1.0E-02
1.0E-03
1.0E-04
1.0E-05

1.0E-04-30

1.0E-05-20-1001020
Prediction Error ∆PX [NM]1.0E-01
1.0E-02

Sample size=52455
Min.=-27.42NM
Max=12.22NM
Median=-1.827NM
Mean=-1.901NM
SD=2.400NM-1001020
Prediction Error ∆PX[NM]30

1.0E+00

図 3 A340-300 の分布
1.0E+00-20

図 6 C 社の分布30

B747-400
Skewness=-0.329
Kurtosis=3.081

1.0E-03
1.0E-04

Log10(Relative Frequency)-30

Log10（Relative Frequency）

Sample Size = 6769
Min. = -26.89 NM
Max. = 7.984 NM
Median = -2.182 NM
Mean = -2.258NM
SD = 2.332 NM
Skewness = -0.359
Kurtosis = 3.056

1.0E-01
1.0E-02

Sample Size = 5526
Min. = -11.69 NM
Max. = 9.147 NMMedian = -1.096 NM
Mean = -1.115 NM
SD = 2.017 NM
Skewness = 0.029
Kurtosis = 1.537

E社

1.0E-03
1.0E-041.0E-05-301.0E-05-30-20-10010Prediction Error ∆PX[NM]図 4 B747-400 の分布

2030-20-1001020Prediction Error ∆PX[NM]図 7 E 社の分布

305. 分布のモデル化
5.1 分布モデルの候補
正規分布，両側指数（DE）分布，非対称 DE（ADE）
分布，正規分布と DE 分布の混合分布（N-DE 分布），正規分布と ADE 分布の混合分布（N-ADE 分布）を予測誤差分布に当てはめた．それぞれの確率密度関数を式(6)から式(10)に示す．
Ⅰ）正規分布f (x) =−12π σe( x −m)22σ 2(6)ただし，m，σ(>0)はパラメータ．
Ⅱ）両側指数(DE)分布

5.2 モデル化の方法
各モデルのパラメータは，最尤推定法[4]で推定した．正規分布と DE 分布は解析的に最大尤度を求めた．
ADE 分布の当てはめでは，パラメータ m を定数とし，中央値を用いた．m は既知として，λ1，λ2 は最尤推定法で解析的に推定した．
N-DE 及び N-ADE 分布では，尤度方程式が複雑であるため，パラメータの最尤推定量を解析的に計算するのは容易でない．このため，実際の計算では，分布のヒストグラムをみて各パラメータを変化させる範囲を決め，その中で区間を分割し近似的に最適な値を探す区間分割法を用いて尤度を計算した．

| x − m|f (x) =1 − λ1e2λ 1(7)

ただし，m， λ 1 (>0)はパラメータ．

Ⅲ）非対称 DE（ADE）分布

5.3 モデル当てはめの結果各分布の当てはめ結果を図 9 から図 13 に示す．また，当てはめの良さを比較するため，AIC（赤池情報量基準）[4]を計算した．AIC は次式で定義される．
AIC= –2{ln(最大尤度) + (自由パラメータ)} 1efor x > m
 2λf ( x ) =  1 | x − m|
 1 − λ2for x ≤ m
 2λ e 2| x − m|−λ1(11)AIC の値はそれぞれ図 9 から図 13 に示した．
一般に，AIC の値が小さいほどモデルの当てはまりがよいと言われている[4]．表 1 に N-DE 分布の AIC で規格化した各モデルの AIC をまとめた．図と表 1 から
N-ADE 分布がよく当てはまることが分かる．

(8)

ただし，ｍ，λ1(>0)，λ2(>0)はパラメータ．
Ⅳ）N-DE 分布

2π σ−e(x − m )22σ2x −mα −+ 1 e λ1λ1(9)ただし，m,σ(>0)，λ 1 (>0)，0<α<1 はパラメータ．

Ⅴ）N-ADE 分布x−m( x − m )2−− (1 − α1 − α 2 ) e 2σ 2 + α1 e λ1
λ12π σf (x) = x −m( x − m )2
 (1 − α1 − α 2 ) − 2α −e 2σ + 2 e λ 2λ22π σfor x > m(10)for x ≤ m

ただし，m, σ，λi (>0), 0 < αi < 1，0 < (α1 + α2) < 1 (i=1, 2) はパラメータ．
一般に，分布の裾部が長い分布を用いて衝突危険度を推定すると，推定値は大きめな見積もりとなる[2]．その代表的なものとして DE および N-DE 分布を当てはめた．分布の非対称性も考慮して，非対称型分布の当てはめも試みた．

1.0E+00

Log10(Relative Frequency)

f (x) =(1 − α 1 )m=-1.340 NM
m=-1.339
σ=2.725σ=2.725 NM
AIC=321224
AIC=3.27e-10

1.0E-01

Fitted

1.0E-02

1.0E-03

1.0E-04

1.0E-05-30-20-10010
20
Prediction Error ∆PX[NM]

図 9 正規分布の当てはめ

30

1.0E+00

m=-1.339NM
m=-1.340
λ=1.967
λ=1.967NM
1AIC=314392
AIC=314383.7

1.0E-01

FittedLog10(Relative Frequency)

Log10(Relative Frequency)

1.0E+00

1.0E-02

1.0E-03

1.0E-01

1.0E-02m= -1.577 NM
σ=2.161 NM
α 1=1.623e-01
α 2=8.825e-02
λ1=2.846 NM
λ2=2.584 NM
AIC=312207

Fitted

1.0E-03

1.0E-04

1.0E-04
1.0E-05

1.0E-05-30-30-20-100102030

Prediction Error ∆PX[NM]-20-1001020
Prediction Error ∆PX[NM]30

図 13 N-ADE 分布の当てはめ

図 10 DE 分布の当てはめ
表 1 AIC の比較モデル
正規DE
ADE
N-DE
N- ADE

Log10(Relative Frequency)

1.0E+00
m=-1.339
m=-1.340 NM=1.968 NM
λλ11=1.968=1.965 NM
λλ22=1.965
AIC=314392
AIC=181685.2

1.0E-01

Fitted

1.0E-02

6. まとめ1.0E-031.0E-04
1.0E-05-30-20-1001020
Prediction Error ∆PX[NM]30

図 11 ADE 分布の当てはめ
1.0E+00m=-1.387 NM
α=0.2761σ=2.124 NM
NM
λ=2.700
1AIC=312373

1.0E-01
Log10(Relative Frequency)

規格化 AIC8851
2019
20190-166

1.0E-02

Fitted

参考文献

1.0E-03

1.0E-04

1.0E-05-30-20

約 3 年 5 ヶ月のデータを用いて，NOPAC における ADS の縦方向予測誤差の分布を調査した．分布の形状は機種や運航者ごとに異なることが分かった．
さらに，正規分布，DE 分布，ADE 分布，N-DE 分
布，N-ADE 分布を，それぞれ全体の分布に当てはめた．AIC の比較では，N-ADE 分布が最もよいモデルであると言える．衝突危険度を計算する場合，空域の特性に応じたモデルとパラメータの選択が重要と考える．今後は衝突危険度モデルに適用され，航空路の安全性評価に役立つものと期待する．

-1001020Prediction Error ∆PX[NM]図 12 N-DE 分布の当てはめ

30[1] 今若善紀：運輸多目的衛星（MTSAT）を利用した航空管制，信学技報，No. SANE 2002-6，
pp.51-54，2002.6．
[2] Reich, P. G.: Analysis of Long-Range Air Traffic
Systems Separation Standards-I, II & III, The
Journal of the Institute of Navigation Vol. 19, 1966.
[3] Anderson, D.: A General Collision Risk Model for
Distance-Based Separation on Intersecting andCoincident Tracks, ICAO SASP-WG/A-WP/21,
Madrid, May 2001.
[4] 広津千尋：実験データの解析−分散分析を超えて−，共立出版，1992.