こんにちは．九大航空修士2年，ITストラテジストの岩附陽太です．今回はBADAデータセットを用いたフライトエンベロープの可視化方法について解説します．自分は深層強化学習による航空管制の自動化を専門としていて，航空機特性を考慮する際にこのBADAデータセットを用いています．

BADAデータセットとは

欧州航空航法安全機構，EUROCONTROL(ユーロコントロールと呼びます)が提供している，様々な旅客機の特性を考慮できるデータセットです[1]．これを用いることで，航空機の失速速度や，質量，燃料消費など様々な要件を考慮することができます．使用のためにはEUROCONTROL の OneSky Online に登録し，ライセンスのリクエストが必要です．本記事ではフライトエンベロープ描画のためのプログラムは公開しますが，ライセンスの関係上データそのものは公開しないように注意しますので，必要な人は上記手順でデータを入手してください．

フライトエンベロープとは

フライトエンベロープは飛行包絡線とも呼ばれ，航空機の水平飛行可能な真大気速度と高度の条件を可視化した図になります．これは論文[2]などでもよく用いられています．航空機は揚力を得て飛ぶため，失速速度以下では飛ぶことができません．また，高度が上がりすぎると大気の密度が低くなるため飛べません．また，旅客機は通常マッハ数を超えては飛行しませんし，速度を上げすぎると動圧限界を超えてしまいいます（動圧限界は構造上安全とされる航空機にかかる圧力のことです．静圧ではなく動圧なので，厳密ではないですが簡単に言うと風から受ける圧力なども考慮した限界です．なので速度が上がりすぎるとダメということですね）．

真対気速度（TAS:True Air Speed）とは

航空機にはGS,IAS,CAS,EAS,TASなどたくさん速度表現があります．これを説明してると長くなるので，説明はほかの記事に譲ります．TASは簡単にいうと，ピトー管から図ったIASを誤差，飛行高度による空気密度を補正した大気に対する速度です[3]．これに風速を考慮すると，GS（対地速度）になり，これは地面に対してどれだけ動いたかなので，レーダーデータやGPSから取得する水平位置と同じ意味になります．

SR比とは

SRとはSpecific rangeの略です．SR比は単位燃料流量で進む，水平距離を表す値であり，下の式で定義されます．
$$SR=\frac{V_{TAS}}{FuelFlowRate}$$
結局トータルの燃料消費量は経路によって大きく変わりますし，降下時は位置エネルギーを速度に効率よく変換できる方が効率的なので単純比較はできませんが，エンルート（全く厳密でないですが，一般の方向けにわかりやすく言えば上昇でも下降でもない，空港間を飛行する際の大部分）の水平飛行の場合はSRが高いと燃料使用効率が良いということになります．

実際に得られる図

こんな図が出せます．

こればB772のフライトエンベロープです．単位がknotとfeetであることに注意してください．高高度ほど速い速度で飛行できます．SRはヒートマップとなっており，高高度，高速のほうが良い傾向があります．

必ずしも高高度高速が良いのか

これだけ見ると航空機は高高度高速の方がよいように見えます．これはいつなんどきも正しいのでしょうか？下図はフライトエンベロープ内の単位時間あたりの燃料流量をヒートマップにした図です．

ご覧のように燃料流量は同高度なら低速の方が小さい，つまり単位時間あたりの燃料消費を抑えたい場合は低速側に振ることが効率的になる可能性がああります．このようなシチュエーションについて説明します．

福岡空港の航空交通流の混雑について

福岡空港は日本で４番目に着陸機数が多く，滑走路１本あたりの着陸機は日本１で世界有数の混雑空港です（航空法でも，IATAのWSGのレベル３にも指定されてます）．下図は国交省の空港別順位から作りました．

実際に空域が混雑するとどうなるのでしょうか？下図は2019年12月15日に福岡空港へ着陸した航空機の航空交通流を可視化したものです．国土交通省航空局のCARATSオープンデータを使っているので，民間機の定期旅客便が対象です．

青色の機体はそれほど経路を延伸していませんが，赤色の機体は大きく経路を延伸していることがわかります．この経路延伸はレーダーベクターと呼ばれて，管制官が航空機同士の安全な距離を保つため，あるいは混雑時に着陸時間を遅らせるために指示する経路延伸になります．この日は青色がRWY16で着陸していて，赤色がRWY34のVisual approachをしている機体が多そうなので，おそらくRWY変更をきっかけに混雑が発生し連鎖したのでしょう．（ちなみに私の専門の深層強化学習による航空交通流の混雑解消はこのレーダーベクターをいかに減らすかという研究です．）

待っているときは低高度低速が理想的

エンルートの際は高高度高速の方がSR比が高く効率的であることは述べました．

さて問題となるのは，この時経路延伸させられて「待っている」航空機はどんな速度で飛ぶべきでしょうか．これは高度を落としてできるだけ燃料消費を減らし，空港に近いところで待っていてほしいので，フライトエンベロープ上では左下に近い部分を飛ぶことになります．SR比最大だけ考えていた場合とは真逆の結果になりますね．この低高度低速で待つ，という概念は今年度に入って福岡にも導入されたポイントマージ上でも同様の管制が行われます．

また，迎角０の際の揚抗比をプロットしてみると以下の図のようになります．やはり低速のほうが揚抗比がよい傾向があり，燃料の節約には利いてそうです．

ただし実際の問題はより複雑であることに注意してください．BADAにおける揚力係数Clおよび抗力係数Cdは迎角0の水平飛行を前提としています．そして揚力係数，抗力係数，揚抗比は迎角によって変化します．

また実際は航空機は高度を落とせば速度が上がります．また，高度によって吹いている風も違います．理想的な管制を行うためにはこれらをすべて航空機の通る軌道を未来予測して（つまり位置の３次元だけでなく時間の次元のことまで考えた）複数機体の４次元軌道考える必要があります．これは人間の管制官には無理です．だからこそ管制官の機能をサポート，代替するＡＩ管制官や管制手法をさまざまな研究者が研究しています．

図を表示するプログラム

長くなるので，BADAのcsv化は別記事にします．BADAのcsv化ができたとして，以下のようなプログラムを用いるとフライトエンベロープの可視化が可能です．


import scipy.io
import csv
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import seaborn as sns

#定数パラメータの設定
Rkitai=287.05287#気体定数
g0=9.80665#重力加速度
earth=6378.137*1000#地球半径(m)
R0=6378.137#地球半径のkm版
pi=math.pi#円周率
kap=1.4#比熱比
Tsea=288.15#海面標準大気状態の温度のケルビン表示
ft=3.2808#1m=3.2808ft
kt=1.9438#1m/s=1.9438kt

TargetKisyu='B772'

Bada=pd.read_csv('BadaOpfIntegrationList.csv',index_col=0)
TargetBada=Bada[Bada['BADAKisyu']==TargetKisyu]
TargetBada=TargetBada.reset_index(drop=True)

mas=float(TargetBada['Reference'][0])*1000
S=float(TargetBada['Surf'][0])#翼面積(m^2)

CD0CR=float(TargetBada['Cd01'][0])
CD2CR=float(TargetBada['Cd21'][0])

CTc1=float(TargetBada['Ctc1'][0])
CTc2=float(TargetBada['Ctc2'][0])
CTc3=float(TargetBada['Ctc3'][0])
CTc4=float(TargetBada['Ctc4'][0])
CTc5=float(TargetBada['Ctc5'][0])
CTcr=0.95

Cf1=float(TargetBada['Cfc1'][0])
Cf2=float(TargetBada['Cfc2'][0])
Cfcr=float(TargetBada['CruiseCorr'][0])
Cf3=float(TargetBada['DescFuel'][0])#fuel decent coefficientのこと
Cf4=float(TargetBada['FlowCoefficient'][0])

VMO=float(TargetBada['VmoKcas'][0])
MMO=float(TargetBada['Mmo'][0])
hMO=float(TargetBada['MaxAlt'][0])
Hmax=float(TargetBada['Hmax'][0])
Vstall0=float(TargetBada['Vstall1'][0])

ro0=1.2249
p0=101325
#ヘクトパスカルでなくパスカルなのに注意
kappa=1.4
mu=(kappa-1)/kappa
CVmin=1.3

#knotとfeetのつもり
vtas=list(range(0,600,1))
alt=list(range(0,60000,1))

lv=len(vtas)
la=len(alt)

geo=[0 for i in range(la)]
geokm=[0 for i in range(la)]

for i in range(0,la):
    print('Phase1:',i)
    geo[i]=earth*alt[i]/ft/(earth+alt[i]/ft)
    geokm[i]=geo[i]/1000

d=[0 for i in range(la)]
Vstall=[0 for i in range(la)]
p=[0 for i in range(la)]
Vpls=[0 for i in range(la)]
T=[0 for i in range(la)]
VMach=[0 for i in range(la)]
Time=[0 for i in range(la)]

for i in range(0,la):
    print('Phase2:',i)
    if geokm[i]<=11:
        #T[i]=11-6.5*kikakm[i]
        T[i]=15-6.5*geokm[i]
        p[i]=p0*(Tsea/(T[i]+273.15))**(-5.256)
    else:
        T[i]=-56.5
        p[i]=22632.064*np.exp(-0.1577*(geokm[i]-11))
    
    d[i]=0.0034837*p[i]/(T[i]+273.15)

    Time[i]=i
    
    Vstall[i]=CVmin * (( (2*p[i]) / (mu*d[i]) ) * ( (( 1+(p0/p[i])*((1+(mu*ro0)*((Vstall0/kt)**2)/(2*p0))**(1/mu)-1))**mu)-1))**(1/2)*kt
    Vpls[i]=(( (2*p[i]) / (mu*d[i]) ) * ( (( 1+(p0/p[i])*((1+(mu*ro0)*((VMO/kt)**2)/(2*p0))**(1/mu)-1))**mu)-1))**(1/2)*kt
    VMach[i]=MMO*np.sqrt(kappa*Rkitai*(T[i]+273.15))*kt

hmax=[0 for i in range(lv)]
for i in range(0,lv):
    print('Phase3:',i)
    hmax[i]=hMO

Cl=[[0 for i in range(lv)] for j in range(la)]
Cd=[[0 for i in range(lv)] for j in range(la)]
Drag=[[0 for i in range(lv)] for j in range(la)]
Thr=[[0 for i in range(lv)] for j in range(la)]
Hisuiryoku=[[0 for i in range(lv)] for j in range(la)]
fcr2=[[0 for i in range(lv)] for j in range(la)]
fcr=[[0 for i in range(lv)] for j in range(la)]
SR=[[0 for i in range(lv)] for j in range(la)]

for i in range(0,la):
    print('Phase4:',i)
    for j in range(0,lv):
        if vtas[j]==0:
            vtas[j]=0.1
        Cl[i][j]=2*mas*(g0*earth**2/((earth+alt[i]/ft)**2))/(float(d[i])*(float(vtas[j])/kt)**2*S)

        Cd[i][j]=CD0CR+CD2CR*float(Cl[i][j])**2
        #単位N(ニュートン)
        Drag[i][j]=(float(Cd[i][j])*float(d[i])*(float(vtas[j])/kt)**2)*S/2

        Thr[i][j]=float(Drag[i][j])
        #単位kg/min-kN
        Hisuiryoku[i][j]=Cf1*(1+float(vtas[j])/Cf2)

        #単位kg/min
        fcr2[i][j]=float(Hisuiryoku[i][j])*(float(Thr[i][j])/1000)*Cfcr
        #単位kg/s
        fcr[i][j]=float(fcr2[i][j])/60
        SR[i][j]=float(vtas[j])/kt/float(fcr[i][j])

temp=la*lv
temp2=[0 for i in range(temp)]
nddata=np.array([temp2,temp2,temp2],dtype=float)


for i in range(0,la):
    for j in range(0,lv):
        count=lv*i+j
        nddata[0][count]=vtas[j]
        nddata[1][count]=alt[i]
        nddata[2][count]=SR[i][j]

nddata2=np.transpose(nddata)

df=pd.DataFrame(nddata2,columns=['vtas','alt','SR'])
df_pivot=pd.pivot_table(data=df,values='SR',columns='vtas',index='alt')

for i in range(la):
    for j in range(lv):
        if j<Vstall[i]:
            df_pivot.iat[i,j]=0
        elif j>Vpls[i]:
            df_pivot.iat[i,j]=0
        elif j>VMach[i]:
            df_pivot.iat[i,j]=0
        elif i>hmax[j]>0:
            df_pivot.iat[i,j]=0
    
fig1=plt.figure()
ax1=fig1.add_subplot(1,1,1)

ax1.plot(Vstall,alt,'-b',zorder=1)
ax1.plot(Vpls,alt,'-r',zorder=1)
ax1.plot(VMach,alt,'-y',zorder=1)
ax1.plot(vtas,hmax,'-g',zorder=1)

im=ax1.imshow(df_pivot,aspect='auto')
fig1.colorbar(im, ax=ax1)
ax1.invert_yaxis()
ax1.set_xlabel('真対気速度[knot]',fontname="MS Gothic")
ax1.set_ylabel('幾何高度[feet]',fontname="MS Gothic")
ax1.legend(["stall","VMO","MMO","hMO"])
ax1.set_xlim(0,600)
ax1.set_ylim(0,60000)
plt.savefig('FlightEnvelopeResult'+TargetKisyu+'.png')
plt.show()

プログラムの説明

まず，このプログラムの使い方はTargetKisyuに求めたい機体の機種名を指定してpythonで実行するだけです．この際にBadaOpfIntegrationList.csvというファイルが必要です．この作り方は別記事にします．

次に，Gridを作ります．Gridを作っているのは下の部分で，600×60000のGridを作っています．こんなにいらないのですが，１刻みにしないとヒートマップと曲線がずれるのでこうしています．計算時間がかかるのは面倒ですが，研究上多くのフライトエンベロープの図は不要で，改修にかかる労力をペイしないと思ったので放置してます．

#knotとfeetのつもり
vtas=list(range(0,600,1))
alt=list(range(0,60000,1))

Phase1では幾何高度をジオポテンシャル高度に直しkm単位にしてます．これはPhase2で標準大気における温度,気圧や大気密度,がジオポテンシャル高度の関数になるためです標準大気．下の部分ですね．

for i in range(0,la):
    print('Phase2:',i)
    if geokm[i]<=11:
        #T[i]=11-6.5*kikakm[i]
        T[i]=15-6.5*geokm[i]
        p[i]=p0*(Tsea/(T[i]+273.15))**(-5.256)
    else:
        T[i]=-56.5
        p[i]=22632.064*np.exp(-0.1577*(geokm[i]-11))
    
    d[i]=0.0034837*p[i]/(T[i]+273.15)

そして，フライトエンベロープの曲線を作ります．失速速度Vstallと動圧限界Vplsは，それぞれ機種ごとのBADAでCASとして与えられていて，Vstall0,VMOがそれになります．これをCASをTASに変換する式を使い変換します．このへんはBADAのUser Manualに全部書いてあります．最新版は3.15とかなのですが，3.10ならネットに転がってますね...．

$$V_{TAS}=[\frac{2p}{\mu \rho} \lbrace (1+\frac{p_{0}}{p}[(1+\frac{\mu \rho_{0}}{2p_{0}}V_{CAS}^{2})^{\frac{1}{\mu}}-1])^{\mu}-1 \rbrace]^{\frac{1}{2}}$$

また，マッハ数の限界はBADAではMMOというマッハ数として与えられるので，音速を計算してknotになおして，マッハ数を出します．この際，ジオポテンシャル高度に対する温度分布は高度11kmで折れ曲がるので，マッハ数のフライトエンベロープも高度11kmで折れ曲がります．それらはこの辺の部分です．

    Vstall[i]=CVmin * (( (2*p[i]) / (mu*d[i]) ) * ( (( 1+(p0/p[i])*((1+(mu*ro0)*((Vstall0/kt)**2)/(2*p0))**(1/mu)-1))**mu)-1))**(1/2)*kt
    Vpls[i]=(( (2*p[i]) / (mu*d[i]) ) * ( (( 1+(p0/p[i])*((1+(mu*ro0)*((VMO/kt)**2)/(2*p0))**(1/mu)-1))**mu)-1))**(1/2)*kt
    VMach[i]=MMO*np.sqrt(kappa*Rkitai*(T[i]+273.15))*kt

Phase3では高度限界を書きます．これは純粋にhMOというfeetでBADAから機種ごとに与えられます．

Phase4ではまず航空機の飛行状態を下図のようにとらえます．経路角＝０です．

一定高度のため重力＝揚力の釣り合いから揚力係数を出し，BADAのフレームワークで揚力係数から抗力係数を導出します．BADAでは巡航，降下，着陸で異なる式を使いますが，巡航を使っています．一定速度を仮定するため，推力＝抗力からエンジン推力を出します．$\rho$は大気密度，$S$はBADAから得られる翼面積です．そしてエンジン推力と比推力$\eta$から燃料流量を出します．BADAではエンジンタイプによって式が異なりますが，今回はジェットエンジンの式を用いています．

$$C_{D}=C_{D0,CR}+C_{D2,CR}C_{L}^{2}$$

$$\eta=C_{f1}\times(1+\frac{V_{TAS}}{C_{f2}})$$

$$ FuelFlowRate=\eta Thrf_{cr} $$

あとは冒頭の定義におけるSR比を，単位を揃えて計算すればOKです．次にヒートマップ描画のためのピボットテーブルを準備します．この時，SRでなく燃料流量を指定すれば燃料流量のヒートマップを作れますし，$C_{L}/C_{D}$のリストをつくれば揚抗比が図示できます．

for i in range(0,la):
    for j in range(0,lv):
        count=lv*i+j
        nddata[0][count]=vtas[j]
        nddata[1][count]=alt[i]
        nddata[2][count]=SR[i][j]

nddata2=np.transpose(nddata)

df=pd.DataFrame(nddata2,columns=['vtas','alt','SR'])
df_pivot=pd.pivot_table(data=df,values='SR',columns='vtas',index='alt')

次にフライトエンベロープの外にある部分のヒートマップの値を0にします．その真大気速度と高度では飛べないことを示すためです．

for i in range(la):
    for j in range(lv):
        if j<Vstall[i]:
            df_pivot.iat[i,j]=0
        elif j>Vpls[i]:
            df_pivot.iat[i,j]=0
        elif j>VMach[i]:
            df_pivot.iat[i,j]=0
        elif i>hmax[j]>0:
            df_pivot.iat[i,j]=0

最後に可視化して終了です．zorderは上書きというか描画レイヤの順番を表しているので，これを書いておかないと塗りつぶされます．ただ曲線から書かないと，ヒートマップが上手くいかなかった気がします．

あとがき

ぼくはもともとMRJと防衛研究がやりたくて2018年に九大に入学しただけに，最近，航空宇宙工学の人気が情報学に負けがちなのを，なんだかなあと思っています．「高校生みんなそんなに情報だけやりたいですか？航空宇宙というか理系分野どこでもデータサイエンスも機械学習もやりますよ．データサイエンスはドメイン知識があってナンボなんで，もっとドメイン側に注目したら―？」と勝手におもってます．そもそも高校生あんまり航空×情報を見たことがないのが原因なのかなと思って，こういったネタを投稿してみました．なかなか核心部分のネタをつかっちゃうと研究上不都合があったりするので，まず周辺分野のネタを投稿してみましたが，これからも自分の研究に支障がないネタは投稿していきたいなと思います．九大の航空宇宙自体，日本で２番目の歴史を持ち，東北大，名古屋と並んで航空宇宙に強い大学として知られていました（東大は別格ですが）．航空管制に関しても，前任の宮沢先生が日本で先駆的に研究を始められて，伝統ある研究テーマとして知られています．先日の航空機事故もあり（事故に関しては本当に心を痛めています．）まだまだ研究し社会実装していかなければならないことが多くある分野です．ぜひ多くの人が航空宇宙や航空管制に興味を持ち，ひいては九州大学を選んでくれたらと思ってます（九大が学研や周船寺に住む分にはそこまで田舎じゃないって説明とかもできたらと思います，愛知から九大に行きましたが，いい大学ですよ！）．