INM による Lden コンターマップ作成マニュアル

INM による
Lden コンターマップ
作成マニュアル
北海道立総合研究機構
環境・地質研究本部
環境科学研究センター
地球・大気環境グループ
2015 年 3 月
目次
１ INM による Lden コンターマップ作成マニュアルについて ......................................... 2
１-１はじめに ......................................................................................................... 2
１-２手順の流れ...................................................................................................... 2
１-３実行に必要なソフトウェア ............................................................................. 2
１-４本マニュアルの表記について ......................................................................... 3
２飛行航跡データの収集 ............................................................................................. 4
２-１ ADS-B 信号受信システムの作成 ....................................................................... 4
２-２航跡データの収集 ........................................................................................... 5
２-３データ処理...................................................................................................... 6
２-４飛行高度の気圧補正 ....................................................................................... 9
３コース選択予測モデルの作成 ................................................................................ 10
４ INM によるシミュレーション手順 .......................................................................... 17
４-１対象空港の条件設定 ..................................................................................... 17
４-２離発着コースデータの作成 ........................................................................... 20
４-３プロファイルデータの作成 ........................................................................... 26
４-４機種・コース別運行本数の作成 .................................................................... 33
４-５騒音評価方法の設定 ..................................................................................... 37
４-６周辺地形の設定 ............................................................................................ 38
４-７シミュレーションの実行 .............................................................................. 39
４-８結果の表示・保存 ......................................................................................... 41
1
１INM による Lden コンターマップ作成マニュアルについて
１-１はじめに
H25 年度より航空機騒音の評価指標が Lden へと移行し、今後環境基準類型の地域の指定範
囲の見直しや新規指定を検討するためには Lden コンターマップが必要です。INM はアメリカ
合衆国連邦航空宇宙庁が制作した航空機騒音シミュレーションソフトウェアで安価（＄３
００）で入手でき、航空機騒音の発生・減衰データが豊富に収録されています。本マニュア
ルは、INM を用いて Lden コンターマップを作成するためのデータ入手と設定の手順を取りま
とめたものであり、安価に Lden コンターマップを作成することができます。
１-２手順の流れ
ADS-B 信号受信システム
想定運航スケジュール
ベイズ統計を用い
飛行航跡データ
気象条件
コース選択予測モデル
たパターン解析
機種・コース別平均運行数
気象データ
シミュレーション設定
速度プロファイル
離着陸コース
高度プロファイル
ＩＮＭによるシミュレーション
⇒Ｌｄｅｎコンターマップ作成
１-３実行に必要なソフトウェア
・INM は FAA より購入できます。
HP アドレス：
https://www.faa.gov/about/office_org/headquarters_offices/apl/research
/models/inm_model/
・本手順では大規模 GIS データを扱うための Python と関連するライブラリおよび
PostgreSQL、地理空間情報を扱うための QGIS と関連するライブラリおよび PostGIS、統
2
計モデル解析のための R および OpenBUGS を使用しています。これらはすべてオープンソ
ースソフトウェアであり、無料で入手可能です。
・必要なソフトウェア
INM, Python, Python ライブラリ(NumPy,Pandas,psycopg2,guiqwt),
PsotgreSQL, PostGIS, QGIS, R, R ライブラリ（R2OpenBUGS）
ベイズ推計ソフトウェア（OpenBUGS）
・INM の詳しい操作方法については、ソフトウェア付属のマニュアルを参照してください
１-４本マニュアルの表記について
本マニュアルでは枠で囲ってある部分は実際の作業を示しています。使用するソフトウ
ェアは枠の最初の行に記載してあります。
3
２飛行航跡データの収集
２-１ADS-B 信号受信システムの作成
① シングルボードコンピューター（RaspberryPi）を用意し、”Raspbian”をインストー
ル、初期設定を行う。
参考 HP：http://qiita.com/makoto_kw/items/393e098f214f81449c9f
② R820T チップを搭載する USB ワンセグチューナーを用意する。
③ RspberryPi に必要なソフトウェアをインストールする。
#/***RaspberryPi コンソール******/
#ライブラリのインストール
$pi>sudo apt-get update
$pi>sudo apt-get upgrade
$pi>sudo apt-get install autoconf automake libtool
$pi>sudo apt-get install git libusb-1.0.0-dev
# rtl-sdr のインストール
$pi>git clone git://git.osmocom.org/rtl-sdr.git
$pi>cd rtl-sdr
$pi>sutoreconf –i
$pi>./configure
$pi>make
$pi>sudo make install
$pi>sudo ldconfig
$pi>sudo make install-udev-rules
# rtl-sdr のインストール
$pi>git clone git://github.com/MalcolmRobb/dump1090.git
$pi>cd dump1090
$pi>make
$pi>cd /usr/local/bin/
$pi>ln –s ~/dump1090/dump1090
4
④ テキストエディタを用いて pi のホームディレクトリに dump1090.sh を作成。
#******RaspberryPi テキストエディタ********
#!bin/si
echo –n “DateTime Set Y/M/D H:M:S
----- “
read ans
sudo date –s “$ans”
./dump1090/dump1090 --net --quiet &
telnet localhost 30003 | tee `date + %Y%m%d%H%M`.log
２-２航跡データの収集
① USB ワンセグチューナーとアンテナを三脚に固定し設置。
設置例
② Raspberry Pi に接続し、起動。
③ ログインし、dump1090.sh を実行。
#/***RaspberryPi コンソール******/
$pi>sh dump1090.sh
5
２-３データ処理
① Python コンソールを起動し、log ファイルを読み込む。
#/***Python コンソール******/
import numpy as np
import pandas as pd
df = pd.read_csv('????.log',header=None,skiprows=3)
df.rename(columns={0:"MSG type",1:"Transtype",2:"SessionID",\
3:"AircraftID",4:"HexIdent",5:"FlightID",\6:"date gen",\
7:"time gen",8:"date log",9:"time log",10:"callsign",\
11:"altitude",12:"grnd speed", 13:"track",14:"latitude",\
15:"longitude",16:"vert rate", 17:"squawk",18:"alert",\
19:"emergency",20:"SPI",21:"on ground"}, inplace=True)
② 不必要なレコードを削除。
#/***Python コンソール******/
df = df[df["Trans type"]!=2]
df = df[df["Trans type"]!=5]
df = df[df["Trans type"]!=6]
df = df[df["Trans type"]!=7]
df = df[df["Trans type"]!=8]
df = df[df["Trans type"].notnull()]
df = df[df.HexIdent!="000000"]
df = df[df.altitude!=0]
hex_alt_min = df.altitude.groupby(df.HexIdent).min()
df = df[np.in1d(df.HexIdent,hex_alt_min.index[hex_alt_min<3000].values)]
③ 使用機材 ID と受信時刻から離着陸を分割し、コールサインを割り振る。
#/***Python コンソール******/
#コールサインを抽出する関数
def res_cs_date(sub_df):
call_ser = pd.Series(sub_df.callsign.dropna().unique())
call_ser = call_ser[np.logical_not(call_ser.str.startswith('?'))]
date_str = sub_df.timestamp.iloc[0].date().strftime("%y%m%d")
return call_ser[0].strip() + "_" + date_str
6
#/***Python コンソール******/
#日付時刻フォーマットデータを作成
df["timestamp"] = pd.to_datetime(df['date log'] + " " +df['time log'])
#HexIdent 毎に処理
hex_ser = df.HexIdent.unique()
for h_hex in hex_ser:
tmp_df = df.ix[df.HexIdent==h_hex]
ts_df = tmp_df[tmp_df["Trans type"] == 3]
peri_start_ts = ts_df[ts_df.timestamp.diff().gt(10*60*1E9)].timestamp
peri_start_arr = pd.Series(ts_df.iloc[0].timestamp).append(peri_start_ts)
peri_end_arr = peri_start_ts.append(pd.Series(ts_df.iloc[-1].timestamp
+ np.timedelta64(1E9)))
for i in range(peri_start_arr.count()):
sub_df = tmp_df[np.logical_and(
tmp_df.timestamp >= peri_start_arr.iloc[i],
tmp_df.timestamp < peri_end_arr.iloc[i])]\
if sub_df.altitude.min() < 1500:
df.ix[np.logical_and(df.HexIdent == h_hex,
np.logical_and(df.timestamp >= peri_start_arr.iloc[i],
df.timestamp < peri_end_arr.iloc[i])),
"cs_date"] = res_cs_date(sub_df)
df = df[np.logical_not(df.cs_date.isnull())]
7
④ 緯度・経度を保有するレコードの速度、機首方向を前後のレコードから補完する。
#/***Python コンソール******/
#cs_date 毎に処理
cs_date_ser = df.cs_date.unique()
for h_cs in cs_date_ser:
df.ix[df.cs_date==h_cs,"track"]\
=df.ix[df.cs_date==h_cs,"track"].fillna(method='ffill')
df.ix[df.cs_date==h_cs,"track"]\
=df.ix[df.cs_date==h_cs,"track"].fillna(method='bfill')
df.ix[df.cs_date==h_cs,"grnd speed"]\
=df.ix[df.cs_date==h_cs,"grnd speed"].fillna(method='ffill')
df.ix[df.cs_date==h_cs,"grnd speed"]\
=df.ix[df.cs_date==h_cs,"grnd speed"].fillna(method='bfill')
⑤ PostGIS にインポートする。
#/***Python コンソール******/
import psycopg2
conn = psycopg2.connect("dbname=???? user=???? host=localhost")
cur = conn.cursor()
tmp_df = df[df.latitude.notnull()]
cur.execute("CREATE TABLE segm_points(id SERIAL, cs_date VARCHAR,\
latitude DOUBLE PRECISION, longitude DOUBLE PRECISION,\
altitude DOUBLE PRECISION,heading DOUBLE PRECISION,\
speed DOUBLE PRECISION, T_log TIMESTAMP)")
cur.executemany("INSERT INTO segm_points(\
cs_date,latitude,longitude,altitude,heading,speed,t_log\
) VALUES(%s,%s,%s,%s,%s,%s,%s)",\
[x[1:] for x in tmp_df[\
["cs_date","latitude","longitude","altitude",\
"track","grnd speed","timestamp"]].itertuples()])
cur.execute("ALTER TABLE segm_points \
ADD COLUMN the_geom geometry(POINT,4326)")
cur.execute("UPDATE segm_points SET the_geom \
= ST_SetSRID(ST_MakePoint(longitude,latitude),4326)")
conn.commit()
8
２-４飛行高度の気圧補正
① 気象庁・過去の気象データ・ダウンロードページより、対象空港に近い官署の海面気圧
時別値を入手。
② Python コンソールを起動し、気圧データを postgreSQL にインポートする。
#/***Python コンソール******/
import numpy as np
import pandas as pd
import psycopg2
conn = psycopg2.connect("dbname=adsb user=postgres host=localhost")
cur = conn.cursor()
ame_df = pd.read_csv("data.csv",encoding = 'SHIFT-JIS',skiprows=4)
ame_df.columns = ['date','pressure','d_flg','no']
ame_df["datetime"] = pd.to_datetime(ame_df.date)
cur.execute("CREATE TABLE ame_pres(
id SERIAL, T_log TIMESTAMP, pressure DOUBLE PRECISION)")
conn.commit()
cur.executemany("INSERT INTO ame_pres(T_log, pressure) VALUES( %s, %s)",
[x[1:] for x in
ame_df[["datetime","pressure"]].itertuples()])
conn.commit()
③ 飛行高度の気圧補正を行う。
#/***Python コンソール******/
cur.execute("ALTER TABLE segm_points ADD COLUMN qnh_alt \
DOUBLE PRECISION")
cur.execute("UPDATE segm_points a SET qnh_alt = \
(SELECT (b.pressure - 1013.25 )/ 1013.25 * 8434.515630756852 /0.3048 + \
a.altitude FROM ame_pres b WHERE date_trunc('hour',a.t_log) = b.t_log)")
conn.commit()
9
３コース選択予測モデルの作成
① 電子航空路誌（eAIP）から Standard Departure Chart 及び Instrument Approach Chart
を入手する。
HP アドレス
https://aisjapan.mlit.go.jp/Login.do
離陸コース
着陸コース
RWY34,RWY16
LOC Z, LOC Y, RNAV
② Python コンソールを起動し、フライトリストを PostgreSQL にインポートする。
#/***Python コンソール******/
import numpy as np
import pandas as pd
#高度が最小値のレコードを抽出する関数
def f_head(df):
return df.sort("altitude").head(1)
#cs_date でグループ化して処理し、csv に出力
df = pd.read_picke(“adsb.df”)
tmp_df = df[df.latitude.notnull()]
apdp_df = tmp_df.groupby("cs_date").apply(f_head)
apdp_df[["cs_date","HexIdent","timestamp"]].to_csv(\
"apdp_list.csv",index=False,date_format='%Y/%m/%d %H:%M:%S')
10
③ フライト毎に飛行コースを判別する。
#/*******QGIS＆Excel******/
Ⅰ．航跡データを QGIS に読み込む
レイヤ→レイヤの追加
→PostGIS レイヤの追加
→segm_points を選択
Ⅱ．cs_date で分類して表示
レイヤのプロパティ→
→スタイルタブ
→”分類された”を選択
カラム：”cs_date”を選択
Ⅲ．cs_date 別に飛行コースを判別
apdp_list.csv に couse 列を追加
QGIS 上で cs_date 別の航跡を
表示し、飛行コースを判別。
11
④ 離発着時刻における気象データを入力する。
#/*******QGIS＆Excel******/
apdp_list.csv に風速、風向、視程列を追
加
気象庁・過去のアメ
ダスデータより、対
象空港に最も近い観
測地点の離発着時刻
のデータを入力
⑤ 気象条件を説明変数、離発着方向を応答変数とした一般化線形モデルを作成する。
#旭川空港におけるモデル
・ILS が R34 方向にしか設置されていないため、視程距離が短い場合、風向に関係な
く R34 に着陸する。そのため、ゼロインフレーションベルヌーイ分布を用いたモデ
ルを作成した。
0.
with prbability 𝑞𝑖
𝑌𝑖 ~ {
𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝑝𝑖 ) with probability (1 − 𝑞𝑖 )
𝑝𝑖 =
𝑞𝑖 =
ここで、𝑌𝑖 : R16 方向の生起
1
1 + exp(−(𝛼1 + 𝛽1 𝑤𝑠𝑖 ))
1
1 + exp(−(𝛼2 + 𝛽2 𝑣𝑖𝑠𝑖 ))
ws𝑖 : R34 方向の風速成分 vis𝑖 : 視程距離
12
⑥ ベイズ推計ソフトウェア OpenBUGS コードによりモデルを記述する。
#/*******Text Editor*******/
#BUGSmodel の作成 → model.txt に保存
model{
for( i in 1:N_dat){
Y[i] ~ dbern(p1[i])
p1[i] <- q1[i] * p2[i]
p2[i] ~ dbern(q2[i])
logit(q1[i]) <- logit_q1[i]
logit_q1[i] <- a1 + b1 * ws[i]
logit(q2[i]) <- logit_q2[i]
logit_q2[i] <- a2 + b2 * vis[i]
}
a1 ~ dnorm(0,0.00001)
a2 ~ dnorm(0,0.00001)
b1 ~ dnorm(0,0.00001)
b2 ~ dnorm(0,0.00001)
}
⑦ R を用いてデータの前処理を行う。
#/*******R コンソール******/
apdp.df <- read.csv("apdp_list.csv",fileEncoding="shift-jis")
apdp.df[apdp.df$couse == "RWY34","op_dir"] = 0
apdp.df[apdp.df$couse == "RWY16",”op_dir"] = 1
apdp.df[apdp.df$couse == "LOCY","op_dir"] = 0
apdp.df[apdp.df$couse == "LOCZ",”op_dir"] = 0
apdp.df[apdp.df$couse == "RNAV",”op_dir"] = 1
apdp.df[apdp.df$w_direction == "北","w_dir_deg"] = 0.
apdp.df[apdp.df$w_direction == "北北東","w_dir_deg"] = 22.5 * 1
apdp.df[apdp.df$w_direction == "北東","w_dir_deg"] = 22.5 * 2
apdp.df[apdp.df$w_direction == "東北東","w_dir_deg"] = 22.5 * 3
13
apdp.df[apdp.df$w_direction == "東","w_dir_deg"] = 22.5 * 4
apdp.df[apdp.df$w_direction == "東南東","w_dir_deg"] = 22.5*5
apdp.df[apdp.df$w_direction == "南東","w_dir_deg"] = 22.5*6
apdp.df[apdp.df$w_direction == "南南東","w_dir_deg"] = 22.5*7
apdp.df[apdp.df$w_direction == "南","w_dir_deg"] = 22.5*8
apdp.df[apdp.df$w_direction == "南南西","w_dir_deg"] = 22.5*9
apdp.df[apdp.df$w_direction == "南西","w_dir_deg"] = 22.5*10
apdp.df[apdp.df$w_direction == "西南西","w_dir_deg"] = 22.5*11
apdp.df[apdp.df$w_direction == "西","w_dir_deg"] = 22.5*12
apdp.df[apdp.df$w_direction == "西北西","w_dir_deg"] = 22.5*13
apdp.df[apdp.df$w_direction == "北西","w_dir_deg"] = 22.5*14
apdp.df[apdp.df$w_direction == "北北西","w_dir_deg"] = 22.5*15
apdp.df["w_cos"] = cos((apdp.df$w_dir_deg + 180
- 334.17)/180*pi)
apdp.df["w_comp"] = apdp.df["w_cos"] * apdp.df["w_speed"]
⑧ OpenBUGS を用いて、パラメータ推定を行う。
#/*******R コンソール******/
library(R2OpenBUGS)
Y <- apdp.df$op_dir
vis <- apdp.df$visibility
ws <- apdp.df$w_comp
N_dat <- dim(apdp.df)[1]
source <- list(“Y”,”vis”,”ws”,”N_dat”)
parameters <- list(“a1”,”a2”,”b1”,”b2”)
inits <- function(){
list(a1 = rnorm(1,0,0.001), a2 = rnorm(1,0,0.001),
b1 = rnorm(1,0,0.001), b2 = rnorm(1,0,0.001))
}
sim <- bugs( data=source, inits, parameters, model.file=”model.txt”,
n.chains = 3, n.iter = 15000)
14
#結果の表示
sim
Inference for Bugs model at "model.txt",
Current: 3 chains, each with 15000 iterations (first 10000 discarded)
Cumulative: n.sims = 15000 iterations saved
mean
sd
2.5%
25%
50%
75% 97.5% Rhat n.eff
a1
-1.4 0.2
-1.8 -1.5 -1.4 -1.2 -1.0
1
430
a2
-1.3 0.7
-2.8 -1.6 -1.2 -0.8 -0.1
1
210
b1
0.9 0.1
0.7
0.8
0.9
0.9
1.2
1
160
b2
0.5 0.3
0.2
0.3
0.4
0.5
1.3
1
93
deviance 312.3 31.5 243.5 291.1 315.6 337.9 359.4
1
85
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 485.3 and DIC = 797.6
DIC is an estimate of expected predictive error (lower deviance is better).
#旭川空港では、R34 の着陸コースが LOCY と LOCZ の 2 種類存在するが、出発地との
関係性がみられないため、航跡データの割合で単純分配する。
p_LOCY = sum(apdp.df$couse == "LOCY")
/sum(apdp.df$couse %in% c("LOCY","LOCZ")) #0.56
15
⑨ 得られたパラメータの値を使用して、コース選択予測モデルを作成する。
#コース選択予測モデル
q(𝑣𝑖𝑠𝑖 ) = 1 −
p(𝑤𝑠𝑖 ) =
1
1 + 𝑒𝑥𝑝(−(−1.3 + 0.5 ∗ 𝑣𝑖𝑠𝑖 ))
1
1 + 𝑒𝑥𝑝(−(−1.4 + 0.9 ∗ 𝑤𝑠𝑖 ))
Pr(𝑅𝑊𝑌34|𝑤𝑠, 𝑣𝑖𝑠) = q(𝑣𝑖𝑠𝑖 ) + (1 − q(𝑣𝑖𝑠𝑖 )) ∗ (1 − p(𝑤𝑠𝑖 ))
Pr(𝑅𝑊𝑌16|𝑤𝑠, 𝑣𝑖𝑠) = (1 − q(𝑣𝑖𝑠𝑖 )) ∗ p(𝑤𝑠𝑖 )
Pr(𝐿𝑂𝐶𝑌|𝑤𝑠, 𝑣𝑖𝑠) = (q(𝑣𝑖𝑠𝑖 ) + (1 − q(𝑣𝑖𝑠𝑖 )) ∗ (1 − p(𝑤𝑠𝑖 ))) ∗ 0.56
Pr(𝐿𝑂𝐶𝑍|𝑤𝑠, 𝑣𝑖𝑠) = (q(𝑣𝑖𝑠𝑖 ) + (1 − q(𝑣𝑖𝑠𝑖 )) ∗ (1 − p(𝑤𝑠𝑖 ))) ∗ 0.44
Pr(𝑅𝑁𝐴𝑉|𝑤𝑠, 𝑣𝑖𝑠) = (1 − q(𝑣𝑖𝑠𝑖 )) ∗ p(𝑤𝑠𝑖 )
16
４INM によるシミュレーション手順
４-１対象空港の条件設定
① AISJapan の公開している AIP より対象空港の情報を入手する。
#AISJapan ページ： https://aisjapan.mlit.go.jp/Login.do
#旭川空港 RJEC Asahikawa(Civil)
ARP coordinates : 434015N/1422651E
Elevation : 690ft
Demendions of RUNWAY : 2500x60(M)
RUNWAY PYSICAL CHARACTERISTICS
RUMWAY NR
TRUE BRG
THR coordinates
THR elevation
16
154.17°
434051.04N/1422626.96E
660ft
34
334.17°
433938.12N/1422715.59
721ft
② INM の Study を作成し、中心緯度経度および高度を設定する。
#/*******INM******/
Ⅰ．File→NewStudy を選択し、
Study を作成。
Ⅱ．使用単位を Engrish System に
設定。
Ⅲ．原点と空港名を設定
Latitude =43+40/60+15/3600
= 43.6708333
Longitude =142+26/60+51/3600
= 142.447500
Elevation=690ft
空港名：AKJ
17
③ シミュレーションを行う機種を登録。
#/*******INM******/
Ⅰ．Setup→Civil Airplane を選択し、シミュレーショ
ンする機種を登録する。
④ 滑走路端の位置情報を設定。
#/*******INM******/
Ⅰ．Tracks→Runway and Helipad
Identifer を選択する。
Ⅱ．Edit→Add Record でレコードを
追加し、滑走路端の識別名を設定
End #1: 34
End #2: 16
Width：60ｍ/0.3048=196.85ft
Ⅲ．Edit→Commit
Record で確定
Ⅳ．Tracks→Runway Ends and Helipads
を選択する。
18
Ⅴ．Runway Ends 毎に位置情報を設定。
Runway:16
Cordinates:Lat/Long
Lat:43+40/60+51.04/3600
= 43.68084444
Long:142+26/60+26.96/3600
= 142.440822
Elavation MSL(ft): 660
Ⅵ．Edit→Commit Record で確定。
Ⅶ．Runway:34 も同様に設定する。
⑤ ケースの作成。
#/*******INM******/
Ⅰ．Setup→Cases を選択する。
Ⅱ．Edit→Add Record を選択し、
Case 名を入力。
Ⅲ．Edit→Commit Record で確定
⑥ シナリオの作成。
#/*******INM******/
Ⅰ．Setup→Scenarios を選択する。
Ⅱ．Edit→Add Record を選択し、
Senario 名を入力、Case を登録。
Ⅲ．Edit→Commit Record で確定
19
４-２離発着コースデータの作成
① Python コンソールを起動し、フライトリストを PostGIS にインポートする。
import numpy as np
import pandas as pd
import psycopg2
conn = psycopg2.connect("dbname=adsb user=postgres host=localhost")
cur = conn.cursor()
apdp_df = pd.read_csv("apdp_list.csv",encoding="shift-jis")
cur.execute("CREATE TABLE apdp_list(id SERIAL, cs_date VARCHAR,
couse VARCHAR)")
conn.commit()
cur.executemany("INSERT INTO apdp_list(cs_date,couse)
VALUES (%s,%s)",
[x[1:] for x in apdp_df[["cs_date","couse"]].itertuples()])
conn.commit()
cur.execute("CREATE VIEW segm_view AS SELECT a.*,b.couse
FROM segm_points a, apdp_list b WHERE a.cs_date = b.cs_date")
conn.commit()
20
② 滑走路端のポイントデータを作成。
#/*******QGIS＆Excel******/
Ⅰ．Excel を用いて滑走路端の位置情報シート
を作成し、CSV カンマ区切り形式で保存。
Ⅱ．滑走路端の位置情報を QGIS に読み込む。
レイヤ→レイヤの追加
→デリミティッドテキストレイヤの追加
インポート設定
ファイル名：RunwayEnds.csv
ファイル形式：CSV
ジオメトリ定義：ポイント座標
X フィールド：Long
Y フィールド：Lat
空間参照システム選択
フィルター：4326
空間参照システム：WGS 84
21
③ QGIS に対象とするコースの航跡データをロードする。
#/*******QGIS＆Excel******/
Ⅰ．航跡データを QGIS に読み込む
レイヤ→レイヤの追加
→PostGIS レイヤの追加
→segm_view を選択
Ⅱ．フィルタの設定
フィルタの設定をクリック→
→フィルタ式を入力
“couse” = 作成するコース名 AND qnh_alt <= 6000（着陸コース）
<= 10000(離陸コース)
→設定ウィンドウを閉じ、追加をクリック。
Ⅲ．表示の変更
レイヤのプロパティ→
→スタイルタブ
→シンプルマーカー
塗りつぶし：#f5f500
境界線：透過
地物混合モード：乗算
22
④ GIS 上で複数の航跡の中心を通るラインベクターを作成。
#/*******QGIS＆Excel******/
Ⅰ．ラインベクターの新規作成
レイヤ→レイヤの作成
→新規シェープファイルレイヤ
レイヤの設定
タイプ：ライン
座標参照システム：WGS 84
レイヤに名前を付けて保存
コース名_line.shp
Ⅱ．ラインの作成
レイヤ→編集モード切替
編集→地物の追加
航跡データの中心を通るよう
に、左クリックでノードを追加
し、ラインを作成してゆく。
滑走路端で右クリックし、入力
を終了（着陸コース）
。
id を設定する。
※離陸コースの場合は、滑走路
端から入力を始める。
Ⅱ．編集を終了
レイヤ→編集モード切替
保存をクリック
23
⑤ ノードの緯度経度を計算する。
#/*******QGIS＆Excel******/
Ⅰ．ラインのノードを展開する
ベクタ→ジオメトリツール
→ノードを展開する
入力ベクタレイヤ：LOCZ_line
出力ファイル：LOCZ_point.shp
Ⅱ．地物 id を追加する
属性テーブル→編集モード切替
→フィールド計算機
地物 id の追加
フィールド名：node_id
フィールドタイプ：整数値
フィールド幅：10
式：$id
Ⅲ．編集を終了
レイヤ→編集モード切替
保存をクリック
Ⅳ．CSV にエクスポート
名前を付けて保存する
→形式：カンマで区切られた値
パス：LOCY_node
レイヤオプション：
GEOMETRY=AS_XY
24
⑥ INM に離着陸 Track を設定。
#/*******INM******/
Ⅰ．Tracks→Imput Graphics を選択。
Ⅱ．Track を入力
Edit→Add Track で Track
入力モードに入る。
右下の緯度経度情報を参
考に大まかな位置で左ク
リックし、始点を作成
右クリックで緯度経度入力
ウィンドウを起動し、ノー
ド毎に緯度経度を入力。
滑走路端でダブルクリッ
クし、Track 識別名を入力。
Edit→Delete Point で
始点を削除。
25
４-３プロファイルデータの作成
① 航跡データのそれぞれのポイントについて、滑走路端までの航行距離を計算する。
#/***Python コンソール******/
import numpy as np
import pandas as pd
import psycopg2
conn = psycopg2.connect("dbname=adsb user=postgres host=localhost")
cur = conn.cursor()
#滑走路端の位置情報を postGIS にインポート
R_end_df = pd.read_csv("RunwayEnds.csv",encoding="shift-jis")
cur.execute("CREATE TABLE runway_ends( \
id SERIAL, name VARCHAR, lat DOUBLE PRECISION, \
long DOUBLE PRECISION, the_geom GEOMETRY(POINT,4326))")
conn.commit()
cur.executemany("INSERT INTO runway_ends(name,lat,long)
VALUES(%s,%s,%s)", \
[(str(x[1]),)+x[2:] for x in
R_end_df[["Runway","Lat","Long"]].itertuples()])
conn.commit()
cur.execute("UPDATE runway_ends SET the_geom = ST_SetSRID( \
ST_MakePoint(long,lat), 4326)")
conn.commit()
#航跡のラインデータを作成する
cur.execute("ALTER TABLE apdp_list ADD COLUMN the_geom \
GEOMETRY(LINESTRING,2454)")
conn.commit()
26
#離陸コースの処理
cur.execute("SELECT cs_date FROM apdp_list WHERE couse = 'RWY34'")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS ( SELECT ST_Transform(ST_MakeLine(\
ARRAY(SELECT the_geom FROM segm_points WHERE cs_date = %s \
ORDER BY t_log)),2454) AS the_geom ) UPDATE apdp_list SET the_geom\
= ST_LineSubstring((SELECT the_geom FROM raw_line),\
ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
(SELECT ST_Transform(the_geom,2454) FROM runway_ends \
WHERE name = '16')),1) WHERE cs_date = %s",[x+x for x in h_list])
conn.commit()
cur.execute("SELECT cs_date FROM apdp_list WHERE couse = 'RWY16'")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS ( SELECT ST_Transform(ST_MakeLine(\
ARRAY(SELECT the_geom FROM segm_points WHERE cs_date = %s \
ORDER BY t_log)),2454) AS the_geom ) UPDATE apdp_list SET the_geom\
= ST_LineSubstring((SELECT the_geom FROM raw_line),\
ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
(SELECT ST_Transform(the_geom,2454) FROM runway_ends \
WHERE name = '34')),1) WHERE cs_date = %s",[x+x for x in h_list])
conn.commit()
#着陸コースの処理（RNAV）
cur.execute("SELECT cs_date FROM apdp_list WHERE couse = 'RNAV'")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS ( SELECT ST_Transform(ST_MakeLine(￥
ARRAY(SELECT the_geom FROM segm_points WHERE cs_date = %s \
ORDER BY t_log)),2454) AS the_geom ) UPDATE apdp_list SET the_geom \
= ST_LineSubstring((SELECT the_geom FROM raw_line),0.0,\
ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
(SELECT ST_Transform(the_geom,2454) FROM runway_ends \
WHERE name = '16'))) WHERE cs_date = %s",[x+x for x in h_list])
conn.commit()
27
#着陸コースの処理（LOCZ,LOCY）
cur.execute("SELECT cs_date FROM apdp_list \
WHERE couse IN ('LOCZ','LOCY')")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS ( SELECT ST_Transform(ST_MakeLine(\
ARRAY(SELECT the_geom FROM segm_points WHERE cs_date = %s \
ORDER BY t_log)),2454) AS the_geom ) UPDATE apdp_list SET the_geom \
= ST_LineSubstring((SELECT the_geom FROM raw_line),0.0,\
ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
(SELECT ST_Transform(the_geom,2454) FROM runway_ends \
WHERE name = '34'))) WHERE cs_date = %s",[x+x for x in h_list])
conn.commit()
#滑走路端からの飛行距離の算出
cur.execute("ALTER TABLE segm_points \
ADD COLUMN dist_km DOUBLE PRECISION")
conn.commit()
cur.execute("SELECT cs_date FROM apdp_list \
WHERE couse IN ('LOCZ','LOCY','RNAV')")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS (\
SELECT the_geom,ST_Length(the_geom) AS len \
FROM apdp_list WHERE cs_date = %s) \
UPDATE segm_points SET dist_km = (\
1-ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
ST_Transform(the_geom,2454)))*(SELECT len FROM raw_line)/1000 \
WHERE cs_date = %s",[x+x for x in h_list])
conn.commit()
cur.execute("SELECT cs_date FROM apdp_list \
WHERE couse IN ('RWY34','RWY16')")
h_list = cur.fetchall()
cur.executemany("WITH raw_line AS (\
SELECT the_geom,ST_Length(the_geom) AS len \
FROM apdp_list WHERE cs_date = %s) \
UPDATE segm_points SET dist_km28
= \
ST_LineLocatePoint((SELECT the_geom FROM raw_line),\
ST_Transform(the_geom,2454))*(SELECT len FROM raw_line)/1000 \
WHERE cs_date = %s",[x+x for x in h_list])
② 着陸について、コース別に滑走路端までの航行距離と高度の散布図を表示し、平均的な
高度プロファイルを作成。
#/***Python コンソール******/
#データの読み込みと処理
import numpy as np
import pandas as pd
import psycopg2
conn = psycopg2.connect("dbname=adsb user=postgres host=localhost")
cur = conn.cursor()
cur.execute("SELECT cs_date,speed,qnh_alt,dist_km FROM segm_points \
WHERE cs_date IN (SELECT cs_date FROM apdp_list \
WHERE couse = 'LOCY') AND qnh_alt < 8000")
df = pd.DataFrame(cur.fetchall())
df.columns = ["cs_date","speed","qnh_alt","dist_km"]
Xa= (df.dist_km/0.3048).values
Ya = (df.qnh_alt - 721).values
#散布図の作成
import guidata
from guiqwt.plot import CurveDialog
from guiqwt.builder import make
from guiqwt.tools import MultiLineTool
app = guidata.qapplication()
win = CurveDialog(edit=True, toolbar=True, wintitle="ScatterPlot",\
options=dict(title="Dist-Alt Scatter",
ylabel="Altitude(ft)",\
xlabel="TrackDistance(1000ft)"))
plot = win.get_plot()
s_plt = make.curve(Xa,Ya,marker="+",linestyle="NoPen")
plot.add_item(s_plt)
win.add_tool(MultiLineTool)
win.show()
app.exec_()
29
Ⅱ．平均的な高度プロファイルを作成
Plyline ツールを用い
て平均的なプロファイ
ルを描画
FreeMarkers ツールを用
いて変曲点の高度・距離
を取得
CSV（カンマ区切り）で保存
③ Procedural Profile のステップを作成。
#/***Python コンソール******/
#高度プロファイルデータの読み込み
prof_df = pd.read_csv("locy_profile.csv",encoding="shift_jis")
#DecendStep 高度の TrackDistance を補完
prof_df = prof_df.append(pd.DataFrame(\
[1000,1500,3000,6000],columns=["Altitude"]))
prof_df.index = prof_df.Altitude
prof_df.sort(inplace=True)
prof_df = prof_df.interpolate(method='values')
#平均速度を求める
prof_df.index = prof_df.Distance
for h_dist in prof_df.Distance:
tmp_df = res_df[np.logical_and(\
res_df.dist_km > h_dist*0.3048 - 0.5,\
res_df.dist_km < h_dist*0.3048 + 0.5)]
prof_df.ix[prof_df.Distance == h_dist,"speed"] = tmp_df.speed.mean()
30
#降下角度を求める
diff_df = prof_df.diff()
prof_df["deg"] =\
np.degrees(np.arctan2(diff_df.Altitude,diff_df.Distance*1000))
#csv に保存
prof_df.to_csv("locy_steps.csv")
④ プロファイルの登録。
#/*******INM******/
Ⅰ.Civil→Profile Identifers を選択
し、機種を選択。APP-STANDARD1 を表示
し、Aircraft Weight をメモしておく。
Ⅱ．Edit→Add Record を選択し、
Operation：APP を選択。Profile
Group 名と Aircraft Weight を入力
して Commit する。
Ⅲ．Civil→Procedual Profiles を選択
し、機種と Profile を選択。
Ⅳ．Edit→Add Record を選択
し、APP-STANDARD１を参考に、高
度 6000ft から順に Step データを
入力。
31
Ⅳ．Civil→Profile Graph を選択
し、作成プロファイルを確認。
32
４-４機種・コース別運行本数の作成
① 一週間単位の運航スケジュールを作成。
#/*******Excel******/
Ⅰ.運航スケジュールを作成し、CSV で保存。
発着地
週運航日数
発着地までの
飛行マイル
発着時間帯の区分
day
使用機材
7:00-19:00
evening 19:00-23:00
night
23:00-7:00
② 気象条件を作成。
#/*******Excel******/
Ⅰ. Excel を用いて風向風速階級別頻度表を作成し、CSV で保存。
列：16 方位（北から右回り）
行：階級幅 1m/s の風速
Ⅱ. 視程距離の階級別頻度表を作成し、CSV で保存。
階級幅 1ｋｍの視程距離
10km 以上は合算
33
③ コース選択予測モデルを用いてコース選択確率を計算する。
#/***Python コンソール******/
import numpy as np
import pandas as pd
#視程の影響の計算
vis_freq_df = pd.read_csv("vis_freq.csv")
vis_mean_ser = vis_freq_df.vis.str.split("-")
vis_freq_df["mean"] = np.nan
vis_freq_df.iloc[:10]["mean"] \
= [np.array(x).astype(np.float).mean() for x in vis_mean_ser[:10]]
vis_freq_df["q34"] = 1 - 1/(1+np.exp(-(-1.3+0.5*vis_freq_df["mean"])))
vis_freq_df.ix[vis_freq_df.vis=="10<","q34"] = 0
q_vis = (vis_freq_df.freq * vis_freq_df.q34).sum() /
vis_freq_df.freq.sum()
#コース選択予測確率の計算
ws_freq_df = pd.read_csv("wind_freq.csv")
ws_freq_df.rename(columns={'Unnamed: 0':'ws'},inplace=True)
ws_freq1d_df =
pd.DataFrame(ws_freq_df.iloc[:,1:].as_matrix().flatten(),columns=["freq"])
ws_mat = np.matrix(np.arange(0.5,12.5,1)).T
cos_mat = np.matrix(np.cos(np.radians(22.5*np.arange(0,16,1) + 180334.17)))
comp_mat = ws_mat * cos_mat
ws_freq1d_df["comp"] = comp_mat.getA().flatten()
ws_freq1d_df["gp"] = pd.cut( ws_freq1d_df.comp,\
np.arange(-12,13,1),labels=np.arange(-11.5,12.5,1))
ws_p_df = pd.DataFrame(ws_freq1d_df.freq.groupby(ws_freq1d_df.gp).sum(),\
columns=["freq"])
ws_p_df["comp"] = ws_p_df.index.values
ws_p_df["p16"] = 1/(1+np.exp(-(-1.4+0.9*ws_p_df.comp)))
p_ws = (ws_p_df.freq * ws_p_df.p16).sum() / ws_p_df.freq.sum()
34
#選択確率の dict を作成
couse_p_dict = {
"RWY34": q_vis + (1-q_vis)*(1-p_ws),
"RWY16": (1-q_vis)*p_ws,
"LOCY": (q_vis + (1-q_vis)*(1-p_ws)) * 0.56,
"LOCZ": (q_vis + (1-q_vis)*(1-p_ws)) * 0.44,
"RNAV": (1-q_vis)*p_ws}
④ 機種・コース別の日平均運航本数を計算する。
#/***Python コンソール******/
#機材毎の離着陸本数を集計
sche_df = pd.read_csv("flight_schedule.csv",encoding="shift_jis")
sche_df["P_Stage"] = pd.cut(sche_df[u"飛行距離（マイル）"],\
np.r_[np.arange(0,2000,500),np.arange(2500,7500,1000),\
[10000]],labels=np.arange(1,10,1))
app_gp = sche_df[u'週運航日数'].groupby(\
[sche_df[u'使用機材'],sche_df[u'到着時間帯']]).sum()
locz_freq =
app_gp.unstack()*couse_p_dict["LOCZ"]/7
locy_freq =
app_gp.unstack()*couse_p_dict["LOCY"]/7
rnav_freq = app_gp.unstack()*couse_p_dict["RNAV"]/7
app_df = pd.concat(dict(￥
LOCZ=locz_freq,LOCY=locy_freq,RNAV=rnav_freq),axis=1)
app_df.to_csv("app_num.csv",encoding="shift-jis")
dep_gp = sche_df[u'週運航日数'].groupby(\
[sche_df[u'使用機材'],sche_df["P_Stage"],￥
sche_df[u'出発時間帯']]).sum()
rwy34_freq = dep_gp.unstack(1).unstack()* couse_p_dict["RWY34"]/7
rwy16_freq = dep_gp.unstack(1).unstack()* couse_p_dict["RWY16"]/7
dep_df = pd.concat(dict(RWY34=rwy34_freq,RWY16=rwy16_freq),axis=1)
#CSV に保存
dep_df.to_csv("dep_num.csv",encoding="shift-jis")
35
⑤ オペレーションを入力。
#/*******INM******/
Ⅰ.Operations→Civil Flights を選択し、
計算する Case を選択。
Ⅱ．Edit→Add Record を選択し、
日平均運航本数データを基に、
Flight Operation データを入力
し、Edit→Commit Record で確定。
36
４-５騒音評価方法の設定
① Noise Metrics に Lden を追加する。
#/*******INM******/
Ⅰ.Setup→Noise Metrics を選択し、
Edit→Add Record で新規 Metrics を作成。
Ⅱ.以下のパラメータを設定し、
Edit→Commit Record で確定。
Metric ID: LDEN
Family: A-weighted
Type: Exposure
Day Multiplier: 1
Evning Multiplier: 3.16
Night Multiplier: 10
10 log(T): 49.37
② 騒音評価地点の設定。
#/*******INM******/
Ⅰ.Setup→Population Points を選択。
Edit→Add Record で新規 Point を作成。
Ⅱ.地点名、緯度経度を入力し、Edit→
Commit Record で確定。
37
４-６周辺地形の設定
① 数値標高モデルから、対象空港周辺の 3TX ファイルを作成。
② INM に読み込む。
#/*******INM******/
Ⅰ.File→Import Data into Study
→Terrain Files を選択。
Ⅱ.Terrian Format：3CD/3TX を
選択し、ファイルの保存されて
いるディレクトリを指定。
38
４-７シミュレーションの実行
① グリッドの設定。
#/*******INM******/
Ⅰ.Run→Grid Setup を選択し、Scinario
を選択。
Ⅱ. Edit→Add Record で新規レコードを
追加。GridType：contor を選択し、範囲
を指定、Edit→Commit Record で確定。
Ⅲ. Edit→Add Record で新規レコードを
追加。GridType：contor を選択し、範囲
を指定、Edit→Commit Record で確定。
② Run Option の設定。
#/*******INM******/
Ⅰ.Run→Run Options を選択。
Run Type: Single Metric
Noise Metrics：LDEN
実行するシナリオ
を選択
Do Terrain をチェック
Terrain Format: 3CD/3TX
Do Line-of-Sight Blokage をチェック
Do Contor をチェック
Do Population Points をチェッ
Ⅱ．Edit→Commit Record で確定
39
③ シミュレーションの実行。
#/*******INM******/
Ⅰ.Run→Run Start を選択。
Ⅱ.実行するシナリオを登録し、
Multisleded をチェックして実行。
Ⅲ.処理状況が表示される。
40
４-８結果の表示・保存
① Output の設定。
#/*******INM******/
Ⅰ.Output→Output Setup を選択。
Edit→Add Record で新規 Output を作成。
Ⅱ.Output 名を入力。Metric に LDEN を選
択し、表示する Senario を選択。
Ⅱ．Edit→Commit Record で確定
② コンター表示とシェープファイルへのエクスポート。
#/*******INM******/
Ⅰ.Output→OutputGraphics を選択し、表
示する Output を選択。
ⅡFile→Export as Shapefile を選択し、
Export Unit を設定して実行。
41
③ 騒音評価地点の結果の表示。
#/*******INM******/
Ⅰ.Output→Noise at PoP Pints を選択し、
表示する Scinario を選択。
42

Download Report