計算アルゴリズム — Python ソースコード

本ページは、本書第10章〜第12章の 1/3 オクターブ受音点レベル L_p(Z)L_p(A)L_p(C)L_p(G) および 重み付け総合値を算出した Python ソースコード一式を整理したものです。各ファイルは Python 標準ライブラリのみを使用しており、NumPy・SciPy 等の外部依存はありません。

計算は次の 3 系統からなります。

動作環境・実行方法

Python 3 のみで動作します。python/ ディレクトリに移動して各メインスクリプトを直接実行してください。

cd python
python3 iso_compute.py             # ISO 9613-2 要約、9 ケース
python3 iso_compute.py --per-band  # 1/3 オクターブ詳細
python3 compute_weather.py         # ケース1〜4 気象補正、9 × 4 = 36 総合値
python3 compute_weather.py --per-band   # 1/3 オクターブ詳細
python3 compute_cylindrical.py     # 円筒拡散による上限値、9 ケース
python3 iso_emit_tables.py         # 付録の標準参照テーブル(α、A/C/G/Z 重み)を LaTeX 形式で出力

ファイル一覧

各ファイルの目的と本書中の位置づけを以下に示す。ソースコードのリンクから生のスクリプトを参照できる。

標準実装層

iso9613_1.py

ISO 9613-1:1993 の大気吸収係数 α(f, T, RH, p_a) を解析公式で計算する。第6.2節の式 (3)〜(5) をそのまま実装。

iso9613_2.py

ISO 9613-2 の点音源屋外伝搬モデル。幾何拡散 A_div、大気吸収 A_atm、地表効果 A_gr を算出する。地表項は海域・陸域の2 セグメントに分けて処理する。

iso_weightings.py

1/3 オクターブバンド中心周波数における周波数重み付けカーブ:IEC 61672-1:2013 の A 特性・C 特性、ISO 7196:1995 の G 特性。

iso_layouts.py

3 配置案(A・B・C)の風車座標と受音点座標を定義する。本書第6.1節・第12.1.1節に対応。

iso_vestas_spectrum.py

Vestas V236-15 MW の 1/3 オクターブ音響パワースペクトル(Z・A)。エストニア Hiiumaa EIA(Ramboll、2022)で採用された安全マージン +2 dB を含む(L_WA = 120 dB(A), L_WZ = 134 dB)。

メインスクリプト・補助

iso_compute.py

ISO 9613-2 ベースラインのメインスクリプト。9 ケース(3 配置 × 3 受音点)の要約と、オプション --per-band による 1/3 オクターブバンド詳細を出力する。表 12.2 を再現。

iso_emit_tables.py

本書付録の標準参照テーブル(10 °C / 70 % および 20 °C / 80 % における α、ならびに A・C・G・Z 重み付け値)を LaTeX 形式で出力する。

weather_corrections.py

ケース1〜4 の帯域別気象補正を実装する。屈折集束 ΔL_refrac(表 11.2.4)、ダクト内拡散 ΔL_duct(実効拡散指数 n_i、表 11.2.3・表 11.2.6)、および海面反射補正 ΔL_sea(表 11.2.5)の 3 補正からなる。

compute_weather.py

ケース1〜4 のメインスクリプト。9 × 4 = 36 通りの気象補正後総合値を算出する。表 12.3 および付録 C を再現。

cylindrical_spreading.py

第11.3節の円筒拡散モデル。ISO 9613-2 の球面拡散 20·log10(r) を、ダクト高 h = 150 m の円筒拡散 10·log10(2π r h) に置き換える。大気吸収(ISO 9613-1)と海面反射 +3 dB はそのまま用いる。

compute_cylindrical.py

円筒拡散モデルのメインスクリプト。9 ケースの要約を出力し、表 12.4 の物理的上限参照行を再現する。

英文の総合 README(インストール条件・モデル定式化・検証結果):python/README.md

計算モデルの要点

ISO 9613-2 ベースライン

各風車 k および各 1/3 オクターブバンド i について、受音点における音圧レベル L_p,i,k を次式で算出する。

L_p,i,k = L_W,i − A_div(r_k) − A_atm,i(r_k) − A_gr,i(r_k)

全風車の寄与を帯域ごとにエネルギー和で集約し、最後に重み付けを加える。

L_p,i  = 10·log10 Σ_k 10^(L_p,i,k / 10)
L_p(W) = 10·log10 Σ_i 10^((L_p,i + W_i) / 10)

重み W_i は標準値:A・C は IEC 61672-1:2013 表 3(0.1 dB 単位)、G は ISO 7196:1995 表 2、Z は全帯域 0 dB。

気象補正(ケース1〜4)

ISO 9613-2 §7.3.1 注 10 では、海域上の気温逆転条件は適用外と明記されている。本書第11.2章では、ISO ベースライン値の上に 3 つの帯域別補正を重ねることで、この問題に対応する。

  1. 屈折集束 ΔL_refrac(表 11.2.4)— 風速勾配および温度勾配により下方に屈折する音響エネルギー。距離非依存・帯域別・ケース依存。
  2. ダクト内拡散補正 ΔL_duct — 表面ダクトに捕捉された音は r^nn < 2)で拡散する。実効拡散指数 n_i を表 11.2.3(海域)・表 11.2.6(陸域)から取り、海域・陸域セグメントを別々に計算したのち合算する。距離依存。
  3. 海面反射補正 ΔL_sea(表 11.2.5)— 強風時(ケース 1, 2)は海面散乱で +3 dB 反射が中高域でやや減少、無風時(ケース 3, 4)は鏡面反射に近づく。

4 ケースは独立した 2 つの気象軸の組み合わせとして定義される。

3 補正はすべて風車ごとのエネルギー和の内側で適用されるため、ダクト項は風車ごとの実際の経路長を反映する。

L_p,i,k_corr = L_p,i,k_ISO + ΔL_refrac,i + ΔL_duct,i,k + ΔL_sea,i
L_p,i_corr   = 10·log10 Σ_k 10^(L_p,i,k_corr / 10)

円筒拡散による上限値(第11.3節)

ISO 9613-2 の球面拡散 20·log10(r) を、表面ダクト高 h = 150 m の円筒拡散 10·log10(2π r h) に置き換えた物理的上限値を算出する。大気吸収(ISO 9613-1)と海面反射 +3 dB は維持される。

本書記載値との対応

本書本文および付録に掲載されている計算結果は、すべて以下のスクリプト群の出力を 0.1 dB に丸めたものである。本ページから取得したコードを再実行することで、本書記載値と次の精度範囲で同一の数値が得られる。

条件設定

すべての計算で次の条件を用いる(ISO 9613-2 参照条件)。

← 主要参考資料に戻る本書のページに戻る