第 539 回 URSI-F 会合資料 2009 年 11 月 20 日

1

レイトレース計算を用いた馬蹄型トンネル内伝搬の解析
Analysis of propagation inside an arched tunnel
using ray tracing simulation
吉敷由起子 1,2
高田潤一 2
Yukiko Kishiki
Jun-ichi Takada
川村 雅彦 2
高尾 肇 3Masahiko Kawamura
Hajime Takao
東京工業大学 1
Tokyo Institute of Technology

1

ナワラット ラーツィリソーポン 2
ギルバート シー チン 2
Navarat Lertsirisopon
Ching, Gilbert Siy
杉原 善洋 3
松永 繁秋 3
上坂 文哉 3
Yoshihiro Sugihara
Shigeaki Matsunaga
Fumiya Uesaka
株式会社 構造計画研究所 2
Kozo Keikaku Engineering Inc.

日揮株式会社 3
JGC CORPORATION

まえがき

新世代移動通信の実現に向けて，あらゆる環境における効率的な伝送研究が盛んに研究されている．最近では伝搬環境によるマルチパスの影響をモデル化できることからレイトレース法が広く利用されつつある [1, 2, 3]．レイト
レース法は FDTD 法 (Finite Diﬀerence Time Domain)や有限要素法よりも広範囲の計算ができ，計算条件の設定も簡易である．またパスごとの到来方向や遅延時間が求まるため，到来方向推定の研究 [4] や ITS(IntelligentTransport Systems) の研究 [5, 6] にも使われている．その他，レイトレース法で計算しにくい散乱体に関しての研究 [7, 8] もされている．
しかしながらレイトレース法には多くの欠点が指摘されている．(1) 幾何光学近似が成り立たない数波長程度以下の物体は扱えない，(2) 曲面や曲線などが扱いづら
い [9]．(1) に関しては幾何光学近似できる範囲を検討し，適用範囲外の散乱については複素散乱断面積を導入することにより解決した [10]．本稿では (2) について検討を行った．レイトレース法では曲面や曲線をそのまま扱うことはできず，平面に分割して計算を行うことが通例である．本研究ではレイトレース法に曲面による反射波と回折波を導入し，円筒面と平面で形成される馬蹄型トンネルにおいて，提案法と実測結果との比較を行った．

円筒面
分割面

円筒面を平面に分割

図1

レイトレース法での円筒面の扱い
yes
円筒面上にパスがない
no

分割面上のパスが全て発生＆ 全パスの長さが同じ

noTx-Rx ラインの分割面にパスがある

円反射なし

yes

円弧の中心点でパスを発生(Caustics パス)
(a)

yes

Tx-Rx ライン上のパスを選択(円反射パス)(b)noTx- 中心、Rx- 中心のラインを生成

２つのエリア内で（Tx- 中心-Rx 扇型、正反対の扇形）それぞれ一番短いパスを選択( 円反射パス) (c)2

円筒面でのパスの軌跡

通常のレイトレース法では構造物のモデル化は平面のみを扱うため，図 1 のような曲率を有する構造物に対して，構造物の構成面を複数の平面板に分割して平面として扱う．また回折波は分割された面の各エッジで発生する．近似精度をあげようとして分割する平面板の大きさを小さくすれば，反射係数を求める際に，”面が波長に対して十分に大きい”というレイトレース法の仮定が成り立たなくなる．
本検討では通常のイメージング法でレイを探索した後に，円筒面内で発生しているレイから反射波と回折波を選び，円筒面の反射係数と回折係数で電界強度を求める．
シミュレータでは円筒面の認識を行うため，円筒面の材質の設定で，”円筒面 材質名”として定義し，円弧の端点と円の中心点の位置情報から円筒面上での反射係数と図2円筒面上での反射波の生成フロー回折係数を求める．

2.1

円筒面上での反射波と回折波

まず円筒面上での反射パスの取り扱いについて述べる．イメージング法でレイを探索すると，円弧上では複数の平板で複数のパスが発生する．フェルマの原理を満足する反射点を探査するため，複数パスから最短距離を選択する．具体的な円筒面上での反射波の生成のフローチャートを図 2 に示す．
簡単に送信点から円筒面上で反射し，受信されるパスを考える．送信点と受信点と円孤の位置により，レイ生成の方法は３種類に分かれる．まず送信点と受信点が円の中心点に配置されている場合には Caustics となる (図rnRxS

Tx
Tx

Rxss'(a) 分割面すべてで反射波発生
Caustic(b) パス長が同じ→Caustics付近
Tx,Rx一直線ライン上のパスを選択

反射面

QR
Rx

円弧(c) それ以外→最短パスを選択

径 ρr1 ，ρr2 を用いて，[11]

図5Tx

図3Rxa) 従来法

図4√

円筒面上での反射波の生成方法

Tx

Txrρr1 ρr2e−jks(ρr1 + s)(ρr2 + s)

(1)

フレネルの反射係数 R̃ はダイアディックで与えられる．
s は QR 点から s 点まで，s ′ は QR 点までの距離である．
平方根の項は曲面による歪みを表す空間減衰係数である．
反射面での主曲率半径 ρr1 ，ρr2 は，入射面での主曲率
半径 ρi1 ，ρi2 と焦点距離 f1 ，f2 を用いて，}{
1 1 111=++(2)ρr1 2 ρi1 ρi2f1}{
1 1 111=++(3)ρr2 2 ρi1 ρi2f2Rxb) 提案法

曲面と平面で形成されるエッジでの回折波球波面で入射するので ρi1 = ρi2 = s ′ となり，焦点距離 f1 ，f2 は以下の式で求められる．
( 2)
11sin θ2 sin2 θ1=+f1,2 cos θia1a2(4)√( 2)21sin θ2 sin2 θ14±+−cos2 θi
a1a2a1 a2sin2 θ1 = cos2 α+sin2 α cos2 θi(5)sin2 θ2 = sin2 α+cos2 α cos2 θi(6)α は入射面と円筒面とのなす角度である．円筒面なので円筒面の曲率半径 a1 ，a2 の一つは無限大で，凹面反射であるためもう一つは負の値となる．ここで a1 を無限大すると，焦点距離は f1 = ∞，f2 < 0 となり，ρr1 ，ρr2 は以下のように簡単になる．

曲率を有する面での散乱係数

点波源から放射された電波は球面状に伝搬するが，曲率を有する面で反射すると，反射波の波形が歪む．曲率半径を用いて，波面の歪みを考慮した円筒面での反射係数と回折係数を導入する．

ρr1 = s ′ ,3.2
3.1iE (S) = E (QR )· R̃

3(a))．円弧上すべての点が反射点，つまり位相が停留することとなり，全分割面において反射波が発生する．
Caustics は幾何光学的な反射係数が発散してしまうた
め，現状は考慮しない．続いて，送信点と受信点と円の中心点が一直線上にある場合には，一直線上付近の円弧上で停留点が停滞しており，複数のパスが発生する (図3(b))．この場合には一直線上のパスを選択する．送信点と受信点の位置が (a)，(b) 以外の時は，送信点と円の中心と受信点から形成される扇型内とその正反対にある扇型内で，２種類のパスが発生する (図 3(c))．各扇型内で述べ距離が最短のパスを円反射パスとする．１つのパスで反射や回折を繰り返す場合には，円弧での反射点を基準としてその前後の散乱点を送信点、受信点に置き換えて考えればよい．
円筒面での回折パスは図 4 のように各エッジで回折波が発生するが，本検討では曲率を有する面と平面で形成されるエッジにて発生する．

3円弧上での反射

円筒面での反射係数

1 1
1= +ρr2 s ′ f2(7)曲面と平面で形成される回折波

回折波は通常平面と平面から成るエッジで発生するが，本検討では円筒面と平面から成るエッジで回折する波を考慮する．円筒面と平面から成る完全導体のエッジでの回折係数 [12] と有限導体のエッジ上での回折係数 [13] か電界 Ei の電波が QR 点で反射した場合（図 5），S 点での電界 Er は，反射点 QR での反射波面での主曲率半2ら，有限導体での円筒面と平面から成るエッジでの回折係数を導出した．他の各パラメータの定義は参考文献に委ねる [12]．

(()
2 sin(π/n)F kLi a(ϕ−ϕ′ )e−jπ/4Ds,h = √2n 2πk sin β0 cos(π/n)−cos ((ϕ−ϕ′ )/n){())π +(ϕ+ϕ′) ( rn +⊥,∥− Rn cot
F kL a (ϕ+ϕ′ )2n)(})′)
π−(ϕ+ϕ) ( r0⊥,∥′+ R0 cot
F kL a(ϕ+ϕ )
2n(8)

5.0m
6.3m

6.0m(a)a(χ) = 2 cos β/2，a (χ) = 2 cos (2πn − β)/2 の値をとる．β0 ，β は入射，出射ベクトルとエッジとのなす角，ϕ′ ，ϕ は入射，出射ベクトルを水平面に射影したベクトルと 0face との角度である．F (x) はフレネル積分，Rn,oはフレネルの反射係数である．入射波に対する距離変数
Li は入射波が球面波である場合には以下のように近似できる．
2Li =+2ss′sin2 β0s + s′

Length: 250m
Position of
Tx & Rx

3.15m 1.0m
5.0m
6.3m
1.8m

6.0m

(9)

(b)

s′ ，s は前回から現在の回折点までの延べ距離と現在から次の回折点までの延べ距離である．反射影境界 (reﬂection
shadow baundaries) での距離変数は
Lr0,n =r0,n+s)ρr0,ns(ρr0,ne
1 ρ2
sin2 β0r0,nr0,n r0,nρe (ρ1 +s)(ρ2 +s)

Length: 300m
Position of
Tx & Rx

(10)

1.0m
5.0m

2.2mr0,nは 0face と nface での出射波の等位相
ここで ρr0,n1 ，ρ2面の主曲率半径である．

12
1r0,n = i − 0,n ′ρ1 a1 (ŝ · n̂0,n )ρ11r0,n =ρ2

1 2(ŝ′ · n̂0,n )−ρi2a0,n2提案法の解析結果

4.1

実測環境

(11)

(c)
図 6 トンネル形状 (a) 実測環境, (b) 馬蹄型トンネルモデル，(c) 方形モデル(12)

反射波が大きく影響を及ぼすため，最大反射回数を何回にするかを検討する必要がある．文献 [15] によると送受信間距離によって有効な最大反射回数が決まる．送受信位置が中央と壁側の両方とも（図 8），反射回数が 4 回以上になるとだんだんと収束に近づいていることがわかった．壁側でも同じような結果となった．反射回数が増えるに従って変動幅は小さくなっている．最大反射回数 2回と最大反射回数 2 回+回折回数 1 回の結果を比較すると回折波を考慮した方が受信レベルが大きくなっている．
この結果から本検討では比較的近距離であるため，最大反射回数が 6 回で十分収束しているといえる．

実験環境は図 6(a) のようなトンネルで行われた．本トンネルの断面サイズ (a×b) は 6.3m×5m で，使用周波数は 5GHz 帯（波長 λ = 0.06m）であり，波長とトンネ
ル断面サイズの関係は a/λ = 103，b/λ = 83 となる．一
般的に 20 以上であれば伝搬形態は幾何学的であるといえる [14]．実験の条件は表 1 に示す．

4.23.0m
6.0m

ρr0,nは出射波の曲率半径で円筒面と平面からなるエッジ
eは直線であるため，ρr0,n= s′ である．
e42.2m4.3

提案法と実測結果との比較

実測結果 [16] と馬蹄型トンネルモデルでの解析（提案法と従来法）を比較した．ここでは馬蹄型トンネルと大きさが等しい方形モデル (図 6(c)) での解析も行った．提案法と従来法の最大反射回数は 6 回で方形モデルも 6 回

提案法の解析結果

馬蹄型トンネルモデル (図 6 (b)) で解析を行った．計算条件は表 2 に示す．トンネル内ではモード伝搬であり，

3表 1 実験条件
周波数4920 MHz(λ = 0.061m)
送受信アンテナ パッチアンテナ
アンテナ高Tx:2.2 m, Rx:2.2 m
送受信間距離0 m 〜 240 m
送受信位置
トンネル断面の中央，壁側 (図 6)
実験回数
中央 3 回，壁側 3 回-50

Path Gain[dB]

-60

Ref:1
Ref:2
Ref:3
Ref:4
Ref:5
Ref:6
Ref:2 Diff:1-90-110
060120Distance[m]180240180240(a)-50である．実測値と各シミュレーションは自由空間損失よりも受信レベルが高くなった．これはトンネル内でモード伝搬になっているためであると考えられる．送受信位置が中央の時，提案法の受信レベルが大きな値を示している．方形モデルのほうが実測値に近い値を示した．送受信位置が中央の場合，円筒面での反射係数が 1 を超えていた．これは送信点と受信点が円の中心位置付近，つまり Caustics 付近となり位相が停留している．そのため受信レベルが非常に大きな値となっていた．方形モデ
ルと実測値の結果がほぼ同じ傾向を示す理由は，方形モ
デルでの反射点が馬蹄型トンネルでの停留点とほぼ同じ位置であるため，同じ傾向を示したと考えられる．
一方，送受信位置が壁側では従来法での受信レベルはだんだんと上がっていき，実測値との隔たりが大きくなっている．提案法は実測値と同じような傾向を示す．
また方形モデルも同じような傾向を示すが，変動幅は方形モデルのほうが大きい．

Path Gain[dB]-60-70-80

Ref:1
Ref:2
Ref:3
Ref:4
Ref:5
Ref:6
Ref:2 Diff:1-90-100-110
060120Distance[m](b)
図 7 トンネル形状での受信レベル 10m 区間の中央値(a) TxRx：中央, (b) TxRx：壁側妥当性を評価した．送受信位置が中央にある場合にはCaustics 現象となり，円筒面上の反射係数が発散してしまい，正しく計算することができなかったが，停留点位置が正しく計算できている方形モデルとは実測値と比較的一致していた．壁側に送受信位置がある場合には，本提案手法が実測値と比較的一致していた．

減衰定数との比較

トンネル内では壁面で反射を繰り返しながら伝搬すると考えられ，減衰定数 α はモード解析から導出された式の近似式より求められている [17]．トンネル断面が幅a，高さ b として，以下の式で求まる．αh と αv は水平偏波，垂直偏波の場合である．
)(1ϵr2√+ √[dB /m] (13)αh = Kh λa3 ϵr − 1 b3 ϵr −1()1ϵr2√√+[dB /m] (14)αv= Kv λa3 ϵr − 1 b3 ϵr −1

参考文献[1] 今井哲朗，藤井輝也，”レイトレース計算を用いた市街地対応移動通信伝搬推定システム，” 信学技報，
RCS97-37, Jan. 1997.
[2] 今井哲朗，藤井輝也，”レイトレースを用いた屋内エリア推定システムの高速アルゴリズム，” 信学論(B), vol.J83-B, No.8, pp.1167-1177, Aug. 2000.
[3] S. Y. Tan, and H. S. Tan, ”A theory for propagation path-loss characteristics in a city-street grid,”
IEEE Trans. on Electromagn. Compat., vol.37,
No.3, pp.333-342, Aug. 1995.

馬蹄型トンネルでは，Kh = 5.13，Kv = 5.09 と報告されており，図 8 には式 (18) から求めた減衰定数の結果を近似直線として示した．減衰定数の近似直線と実測値の傾きは一致している．

5-80-100

表 2 シミュレーション計算条件
最大反射回数
1 回〜6 回
最大回折回数
0 回〜1 回材質
コンクリート
材質の電気的特性 ϵr :6.76,µr :1.0,σ:2.3E-3 [S/m]4.4-70[4] T. Fugen, J. Maurer, W. Wiesbeck, ”Cluster
characterization in urban macrocellular environments with ray-tracing,” IEEE Vehicular Technology Conference (VTC2005-Fall), Dallas, Texas,
Sep. 2005.

まとめ

本研究ではレイトレース計算に円筒面での反射波と回折波を導入し，実測結果と比較することによりその4Pat-90

Measurement (1m median value)
Proposed (1m median value)Proposed (10m median value)
Conventional (1m median value)
Square (1m median value)Free space-100-110
0-5060120
Distance[m]180240Introducing the Complex Radar Cross Section,”International Symposium on Antenna and Propagation (ISAP2008), pp.790-793, Taipei, Oct. 2008.

-50-60-60

Path Gain[dB]-70[11] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley and Sons, New York, 1989.

Path Gain[dB]-70-80-80-90-90-100-100-110
1-110[12] R. G. Kouyoumjian and P.H. Pathak, A Uniform
Geometrical Theory of Diﬀraction for an Edge in aPerfectly. Conducting Surface, Proc. IEEE, vol.
62 , pp. 1448-1461, Nov. 1974.

Measurement (10m median value)
Proposed (10m median value)
Measurement(1m(10mmedianvalue)value)
Conventional
medianProposedvalue)
Square (1m(10mmedianmedianvalue)
Proposed(10m median value)Free space
Conventional(1mmedian value)
Attenuation
coeff.
Square(1m medianvalue)Free space
10100[13] R.Luebbers Finite Conductivity Uniform GTD
versus Knife Edge Diﬀraction in Prediction of
Propagation Path Loss, Trans. IEEE, vol. AP32 , no.1, pp. 70-76, Nov. 1984.

1000

Distance[m]

060

120
Distance[m]180240(a)-50[14] 細矢良雄監修，電波伝搬ハンドブック，リアライズ社，東京，1999.

Path Gain[dB]

-60
-70

[15] 今井哲朗, レイトレース法を用いたトンネル内伝搬特性の推定, 信学技報, vol. J85-B, no. 2, pp.
216-226, Feb. 2002

-80

Measurement (10m median value)
Proposed (10m median value)
Conventional (10m median value)
Square (10m median value)Free space
Attenuation coeff.

-90
-100
-110
1

10100

[16] 高田潤一，高尾肇，杉原善洋，松永繁秋，上坂文哉,前島嘉明, 市塚善己, ”4.9GHz 帯無線 LAN のためのトンネル内伝搬特性の測定,” 第 539 回 URSI-F
会合, Nov. 2009

1000

Distance[m][17] 山口芳雄, 阿部武雄, 関口利男, ”任意断面をもつトンネル内電波減衰定数の近似式について,” 電子通信
学会論文誌 B, 技術談話室, Vol.J67-B, no.3, pp.352- 353 March 1984(b)
図 8 実測との比較結果 10m 区間中央値 (a) TxRx：中央, (b) TxRx：壁側[5] 一之瀬寛智, 尾崎 敏, 津留孝俊, 橋本 修，”複数レーン ETC 料金所における電波環境に関する検討，” 信学技報, ED2005-220, MW2005-174, pp.63-68, Jan.
2006.
[6] K. Ito, N. Itoh, K. Sanda, Y. Karasawa, ”A NovelMIMO-STBC Scheme for Inter-Vehicle Communications at intersection,” IEEE VTC-2006, Sydney
Australia, May 2006.
[7] 中川真弓，藤井泰憲，内田一徳，”2 枚の半無限平行導体板による電磁波散乱問題の厳密解との比較による簡単化レイ・トレース法の精度について，” 信学論 (C)，Vol.J89-C, No.11, pp.951-955, Nov. 2006.
[8] H.Budiarto, K.Horihata, K.Haneda, J.Takada,”Experimental Study of Non-specular Wave Scattering from Building Surface Roughness for the
Mobile Propagation Modeling(Antenna and Propagation), ” IEICE Trans. Com.，Vol.E87-B, No.4,pp.958-966, Apr. 2004.
[9] 安藤真，”開口面アンテナの指向特性解析法 −物
理工学近似と開口面法の正確な理解と使い分け−,”信学論 (B), Vol.J86-B, No.9, pp.1694-1705, Sep.
2003.
[10] Y.Kishiki, J.Takada, ”Improvement of 3D Ray
Tracing Simulation in Microcell Environment by5