セルオートマトン法を用いた騒音伝播モデルの開発
建設設計工学研究室 富樫 孝介
指導教官 宮木 康幸

1はじめに騒音は日常生活に密接な関係があり，関心の高い環境問題の一つである．その中でも，建設業界に起因する騒音問題は，平成 14 年度騒音規制法施行状況調査（環境省）によると，苦情件数 15,464件に対して，建設作業 27.0%，航空機 4.7%，自動車 2.5%，鉄道 0.6%の割合で発生しており，合計34.8%と大きな割合を占めている (図–1 参照)．

図–1 騒音に対する苦情 平成 14 年度騒音規制法施行状況調査 (環境省)

このような状況の中で，特定の騒音に対する騒音伝播解析や騒音制御技術の開発は盛んに行われている．しかし，汎用性を持った騒音伝播予測技術はあまり開発されていない．そこで，本研究は汎用性をもった 2 次元騒音伝播予測モデルの開発を目的とする．
騒音伝播解析には，有限要素法や差分法など，現象を支配方程式によりモデル化し解析する手法が主に用いられてきた．しかし，音の伝播は非常に複雑であり，結果を得ることは非常に難しい．そこで本研究では，近年，機械・電気の分野で注目されているセルオートマトン法 (Cellular Automata：以下 CA) を用いて解析を行うことにした．
本モデルの妥当性の検討に際しては，1 次元に関して，1 次元音響管を用いて開管と閉管の 2 つの場合について，波動方程式の解と CA による計算解を比較する．2 次元については，自由音場における距離減衰の検討として，理論値と比較を行い，また，半無限障壁の回折減衰について Fresnel–Kirchhoff の光の回折理論と前川の実験式の 2 つと比較を行う．

2セルオートマトン法

CA は解析対象をセルと称する均一な大きさの区分領域に分割し，各セル上に定義された離散的な状態量をごく近傍のセル間に設けた局所近傍則によって，離散的時間を追って状態量を推移させる，現象のモデル化手法である．
CA は群集流解析や粒状体解析，流体解析などさまざまな分野への適用が試みられている．また，従来から方程式によるモデル化が可能だった物理現象に対しても本手法は適用されており，波動・拡散現象など数学モデルとして偏微分方程式で記述されてきた現象について解析された例がいくつかある．
局所近傍則によって，系全体としての現象を表現可能となるため，系全体について解析を行うことなく一部分を抽出することが可能となる．この局所近傍則は任意に設定することができ，解析対象が物理現象を表すものであれば物理現象を近傍則として扱うことも可能である．しかし，その一般的な導出方法は確立されていない．局所近傍則は支配方程式に代わるものであるため，解析対象によって十分に吟味する必要がある．また，セルに与える状態量に変位や速度などの物理量を定義することが出来るが，2 次元 CA の場合，解析空間は正方形もしくは正三角形のセルで分割するため，状態量に方向性が与えられる．そこで，空間全体についてある種の平均化操作が必要となる場
合がある．

3
3.1

理論
1 次元音響管の支配方程式

媒質に生ずる圧力の振動は音として観測され，音場は以下の波動方程式により記述される．

∂p(x, t)= c2 ∇2 p(x, t),∂t2sc=Kρ(1)

ここで，p(x, t) は音圧，c は音速を表し，K, ρ はそれぞれ媒質の体積弾性率および密度である．
調和振動する音源によって形成される音場は，調和振動解を仮定することにより一般的に次式で求まり，平面波の場合，次式で表される．
no(2)p(x, t) = jωρ Aej(ωt−kx) Bej(ωt+kx)
ここで，A, B は境界条件より定まる係数であり，ω は音源の角周波数，k は波数 (k = ω/c) である．式 (2) の第 1 項は進行波，第 2 項は後進波を表す．一端に加振端を持つ音響管の音圧に関する調和振動解は，加振端の境界条件として粒子速度u̇(0, t) = u̇0 ejωt を与え，さらに他端が開放端の場
合，p(l, t) = 0 の条件より，p(x, t) = jρcu̇0

sin k(l − x) jωtecos kl(3)と求まる．ここで l は管の長さを表す．また，閉端の場合も同様に，境界条件 u̇(l, t) = 0 を考慮して，p(x, t) = −jρcu̇0cos k(l − x) jωtesin kl(4)と求まる．管の共鳴周波数は，各々の境界条件において，f = mc/2l および f = mc/4l (m =1, 2, · · · ) である．
以上のように，比較的境界条件が単純な場合や媒質の性質が一様である場合については波動方程式を解析的に解くことが出来る．
3.2 自由音場における空気音の距離減衰
3.2.1 点音源からの距離減衰
3 次元の自由空間中にある音響出力 W の点音源から，距離 d の点の音の強さ I は，半径 d，球面の表面積 4πd2 を単位時間に通過するエネルギの総和が W であることより，

I=W4πd2(5)となり，距離の 2 乗に比例する．この点における音の強さのレベルは，音響出力をパワーレベル Lwで表して，

L = Lw − 10 log10 4π − 10 log10 d2= Lw − 11 − 20 log10 d[dB](6)

音源に指向性があって，受音点の方向に対して指向係数 Q であればL = Lw − 11 − 20 log10 d + 10 log10 Q [dB] (7)となる．この場合もひとつの方向に対しては，逆2 乗則が成立するから，距離 d1 の点で音圧レベル
L1 [dB] であれば，距離 d2 = nd1 の点ではL2 = L1 − 20 log10

d2 2= L1 − 20 log10 nd1 2[dB](8)となり，d が 2 倍になるごとに 6dB 減衰する．
3.2.2 線音源からの距離減衰
自動車道路のように，音源が多数並んでいるような場合は，無限長線音源として扱うことができる．位相がランダムで波動性が無視できるような音源が，無数に連続していると仮定し，その単位長さ当たりの音響出力を W とすると，音波はその線音源を軸とする円筒状に広がるから，距離 dの点の音の強さは，半径 d の円筒面を考え
W
I=(9)2πdとなり，距離の 1 乗に反比例する．その強さのレベルは単位長さ当たりの音源出力をパワーレベル
Lw と表しL = Lw − 8 − 10 log10 d[dB](10)となって，距離が 2 倍になるごとに 3dB の減衰をする（3dB/DD と表記することができる）．
有限の長さの線音源の場合には，点音源が x1 から x2 まで連続していると考え，垂直距離 d にある受音点 P の受けるエネルギ密度はZ x2W dx
E=2x1 4πr cZ x2Wdx=24πc x1 (d + x2 )W 1 ³ −1 x2x1 ´=·tan− tan−14πc dddW ψ=·(11)4πc dとなって，P から線音源 x1 · x2 を見る角 ψ に比例し，距離 d に反比例する．
音源の長さ l，x1 = −x2 = l/2 のときの P 点の音圧レベルの距離 d による変化は，式 (11) より，µ¶W 1−1 lE=·tan2πc d2dと求めることができる．このとき，
¶µπl
W 1−1=d ¿ l;tan∴ Enear =·2d24cd¶µllWl 1 
=d À l;tan−1∴ E f ar =·
2d2d4πc d2(12)この両者を等しいとおけば，d = l/x，これは，dが l/π より近い範囲は，倍距離 3dB，遠い範囲は倍距離 6dB の減衰で近似できることを示す．
3.3 半無限障壁による回折減衰
野外において音源と受音点の間に十分大きい塀を立てると，音源側は壁による反射で音圧が上がるが，受音点側では影となって遮音効果が期待できる．この影の中の音圧レベルは，光の回折に関する
Fresnel–Kirchhoff の近似理論やそれをもとに前川が実験的に修正した図（通称前川チャート，図–2）によって求めることができる．

図–2 自由空間の薄い半無限障壁による減衰値

光の回折に関する Fresnel–Kirchhoff の近似理論
において半無限障壁をとりまく空間は完全な自由音場を仮定して，2 次元平面上で点音源を仮定している．これは，3 次元空間で 1 つの線音源を仮定し，それを単位長さで平面に切り取ったものに等しい．前川の実験は，3 次元空間に点音源を配置したものであり，実験の制約上，完全な自由音場ではないと思われる．そのため地面からの反射も実験値に含まれている．
前川チャートより得られた近似曲線を定式化したものを以下に示す．
µ ¶N− 20.0 (1.0 ≤ N )∆L= −10.0 log105∆L= −4.97 log10 (N ) − 13.01(0.1 ≤ N < 1.0)
∆L= −2.09 log10 (N ) − 10.12(0.01 ≤ N < 0.1)
∆L= −77.0 · N − 5.17(|N | < 0.01) !
µ¶ Ã0.33−0.3 < N∆L= 10.0 · log10 √3N ≤ −0.01
N + 1.0∆L= 0.0(N ≤ −0.3)(13)

ここで，N はフレネル数を表し，N = 2δ/λ で求めることができる．また，λ は波長，δ は経路差であり，δ = A + B − d で表される．経路差の模式
図は図–2 の右下部に示す．

4CA による音場のモデル化

1 次元および 2 次元 CA を用いて音場のモデル化を試みる．
4.1 解析空間の設定
1 次元音響管を図–3(a) に，2 次元自由音場を図–4(a) に示すように領域を正方形のセルで分割した．
各セル上に状態量として，セルの状態を壁・媒質(空気)・解析境界の 3 種に区別する整数値 S を定義した．さらに壁でないセルに対して音圧 P およびその時間変化量 G を実数値で定義した．考慮すべき近傍セルとして，1 次元においては図–3(b) に示すように両側に存在する 2 セルを設定した．2 次元の場合は，図–4(b) のように近傍セルを上下左右にとるノイマン近傍を採用した．CA の計算においては空間および時間は離散的な整数値として扱うが，波動方程式との対応づけを行う必要がある．ここでは，単位セル幅を dx = 0.01m と決め，波動の媒体を空気として音速を c = 344m/s と定義した．
時間ステップ dt は dx を用いて，dt = ca · dx/c で求めることが出来る．表–1 に実現象における諸物理量と CA 空間におけるそれとの対応を示す．

図–3 1 次元音響管の要素分割と近傍セル

図–4 2 次元自由音場の要素分割と近傍セル

5

表–1 諸物理量の対応
音速c = 344 m/s
ca = 1 cell/step√ca = 1/ 2 cell/step

実現象
1 次元 CA
2 次元 CA

時間ステップdt = 1/344 sec
dt = 1 step
dt = 1 step

単位セル幅dx = 0.01 mdx = 1 celldx = 1 cell

4.2 局所近傍則
各セルに与える状態量としては 4.1 節に示すよ
うに，音圧 P と時間変化量 G を定義した．これらの状態量を計算ステップを追うごとに更新するための局所近傍則は，(1) 音圧値 P の次の時間ステップにおける時間変化量 G は，注目セルと近傍セルとの音圧値の差から直接決定される．
(2) 次の時間ステップにおける音圧値 P は，そのセルの現在の音圧値に (1) で求めた時間変化量 G を加算することで決定される．
(3) 近傍セルが壁の場合の音圧値 Pwall は，注目セルと同じ音圧値 Pnotice を用いることで完全反射を考慮する．

解析結果および考察

5.1 1 次元音響管
前節で示した局所近傍則を用いて CA により解析した音場モデルと波動方程式の解を比較することで，本手法の妥当性を検討する．
図–5 に示すように，音響管の境界条件として左側を加振端，右側を固定端および開放端の 2 種類を考え，強制加振時の管内の音圧分布を求める．波動方程式による解は調和加振であれば，式 (3),(4)より求めることが出来る．ここで，管の長さ l =2.0m，管路の分割は 200 セルとした．また，波動は管の断面積によらず平面波が生じるものとした．
CA を用いた解析において，セルの状態量として粒子速度を定義していないため，音源に直接音圧値を与えた．この音源の音圧値は，入力する音圧レベルを 100dB とし，音源の最大音圧 Pmax =100
P0 × 10 20 で求める．ここで，P0 は基準音圧でP0 = 2.0×10−5 Pa である．調和加振は余弦波 (cos)を仮定しているので，P (0, t) = Pmax × cos(ωt) で表現した．ここで ω は角周波数であり，周波数 fと時間 t を用いて ω = 2πf で求められる．

具体的に局所近傍則の定式化を行うと以下のようになる．まず時間変化量 G は以下のようになる．

G(x, t + 1) = G(x, t)()
X+ma P (x + dxa , t) − mP (x, t) ,図–5 音響管の形状

aX(14)
ma = ma

さらに，音圧の更新は次式に従う．

P (x, t + 1) = P (x, t) + G(x, t + 1)

(15)

波動が等方的に伝播するとして式 (14) において，ma = m/N(16)とする必要がある．N は近傍数である．m は差分方程式における差分スキームのようなものであり，以下のような式で表すことが出来る．

m = N · c2a(17)ca は伝播速度に相当するパラメータである．設定
可能な伝播速度は，1 次元において ca 5 1，2 次√
元においては ca 5 1/ 2 である．

閉管の共振周波数におけるモード形状の比較を図–6 に示す．比較を行ったモードは開管・閉管ともに 1∼3 次である．実線が波動方程式の解より得られた解析解であり，破線が CA により得られた計算解である．同様に開管の共振周波数におけるモード形状の比較を図–7 に示す．開管に関しては，左の縦軸が波動方程式の音圧を示しており，右の縦軸が CA の音圧値を示している．
閉管に関しては，おおむね正確なモード形状と音圧値が得られている．開管に関しては音圧値のオーダーが異なるものの，モード形状に関しては十分な精度が得られている．閉管でオーダーが異なる結果が得られた原因として，境界条件を壁に当たる直前のセル P (l/dx) と壁セル P ((l/dx)+1) の音圧値を同じにすることで，時間変化量 G(l/dx) = 0を表現しようと試みた点に問題があると考えられる．また，気柱の共鳴の実験においては，音源は管の入り口に数秒間かざすだけで共振が得られる

図–6 閉管の共振周波数におけるモード形状の比較

図–8 点音源を解析領域中央に配置したときの伝播状況

図–7 開管の共振周波数におけるモード形状の比較

が，本解析モデルでは，音源は常に音を出し続けるため，共鳴が連続的に発生することで，音が発散に向かったのではないかと考えられる．

5.22 次元自由音場

2 次元の自由音場に点音源が 1 つ存在する状態を想定して音場のモデル化を行った．これは 3 次元空間に 1 つの線音源を仮定し，それを単位長さで平面に切り取った状態に相当している．
まず，解析領域の中央部に音源を配置し，そこから出される余弦波の伝播形状の妥当性の評価を行った．領域は縦・横 2m で切り取り，領域の分割は縦・横方向ともに 100 セルとした．そして，1セルあたりの幅を dx, dy = 0.02m として解析を行った．このときの音源の扱いは 1 次元音響管の場合と同じ方法を用いた．音源の周波数を 500Hz，計算ステップ数 900step としたときの伝播形状を図–8 に示す．同心円状に音が伝播しており，よい結果が得られている．
次に距離減衰の検討を行った．解析領域を縦・横1000 セル，分割幅 dx, dy = 0.02m とし，音源をその中央に配置した．受音点は音源からの距離をl=120[cell],2l=240[cell] と設定して解析を行った．

図–9 距離 l と 2l の間の距離減衰量

音源の周波数を 100Hz から 2400Hz まで 100Hz 刻
みで l 地点と 2l 地点の距離減衰による音圧レベル差をプロットしたものが図–9 である．検討を行った 24 種の周波数において，平均で 2.97dB の減衰が見られた．本解析モデルで得られた距離減衰は，空間 3 次元に線音源を配置したときの距離減衰量3dB/DD が得られていることがわかった．本解析
モデルは，空間 2 次元に点音源を配置したものであるが，これを空間 3 次元に拡張した場合，奥行き方向に 2 次元モデルを並べたものとなるため，空間3 次元に線音源を配置したモデルと等しい解が得られたと考えられる．理論値との差は最大 0.17dBと小さく，2 次元自由音場における音の距離減衰を CA を用いて表現できると言える．
最後に，半無限障壁の回折減衰の検討を行った．
前川チャートによる回折減衰量は，壁の有無によって変化する受音点の音圧レベルの差を表している．
これと同じ状況で CA を用いて解析を行った．解析領域は，余弦波の伝播形状の妥当性の評価と同じく

Fresnel–Kirchhoff の光の回折近似理論と前川の実験式との間に挟まれる形の結果が出ている．全体的に Fresnel–Kirchhoff の光の回折近似理論に近い値を示しているのは，本研究の境界条件が Fresnel–Kirchhoff の光の回折近似理論と同じく解析領域の外周を無限境界で考えており，壁も下方向に無限（完全な半無限障壁）にとっている．それに対して前川の実験は，現実の空間で行われた実験であり，壁が有限の長さを持っており，かつ，地面による反射の影響を受けている．この違いが，解析結果に表れたのではないかと推測できる．

6
図–10 壁がある場合の伝播状況

CA を用いることで 2 次元空間における騒音レベルの伝播予測ができるということがわかった．

7

図–11 回折減衰量の比較設定した．解析領域の右下部（(x,y)=(25,75)[cell]）に音源セルを配置し，受音点は壁セルを軸に音源と線対称（(x,y)=(75,75)[cell]）に配置した．近傍セルに壁セルが存在する場合，完全反射をするという仮定の下，その壁セルは注目セルと同じ音圧値をとるという形で壁セルに関する境界条件を決定した．解析は始めに壁セルを配置しないパターン（自由音場）で行った．次に解析領域の中央部（(x,y)=(50,50)[cell]）から下方（(x,y)=(50,100)[cell]）に壁セルを配置したパターン（半無限障壁）で解析を行った．音源の周波数を 500Hz，計算ステップ数 900step としたときの伝播形状を図–10 に示す．そして 2 つの解析より得られた受音点の音圧レベルの差をとって得られた減衰量と前川の実験式および Kirchhoff の近似理論より算出された音圧レベルの差を音源の周波数を 100Hz から 2400Hz まで 100Hz 刻みで比較したのが図–11 である．CA による計算解はおおむね

結論

今後の課題

本研究では，2 次元自由音場における解析境界の表現がうまく定式化できなかった．無限境界の定式化が出来ることで，解析精度の向上・解析領域の拡大・計算時間の短縮などが期待できる．また，壁については完全反射のみ考慮したモデル化を行ったが，吸収や分散などの減衰効果を定式化することで，汎用性が向上すると考えられる．そして，本モデルを実現象や他の解析手法と比較することで，本解析モデルの信頼性を確認出来ればよいと考えている．最後に，計算効率の最適化は行っていないので，この点を改良することで，解析領域の拡大・計算速度の向上が図れればよいと思う．
参考文献
1) 小松崎 俊彦ら:セルラオートマトンによる波動伝播解析, Transactions of JSCES, No.19990017, 1999.
2) 加藤 恭義ら:セルオートマトン法 —複雑系の自己組織化と超並列処理—, 森北出版, 1998.
3) 佐藤 清雄:振動と波動, 初等物理シリーズ 2, 培風館,
pp114–149, 1993.
4) 高見 穎郎, 河村 哲也:偏微分方程式の差分解法, 東京
大学基礎工学双書, 東京大学出版会, pp1–74, 1994.
5) 中野 有朋:低騒音化技術, 技術書院, 1993.
6) 前川 純一, 森本 政之, 阪上 公博:建築 · 環境音響学
第 2 版, 共立出版, 2000.
7) 前川 純一:障壁 (塀) の遮音設計に関する実験的研究,
日本音響学会誌, Vol.18 No.4 pp.187–196, 1962.
8) 前川 純一:有限障壁 (衝立) による騒音の回折—障壁(塀) の遮音設計に関する実験的研究 第 2 報, 日本音響学会誌, Vol.21 No.1 pp.1–7, 1965.