第 480 回電波研連Ｆ分科会
2003 年 12 月 12 日

模擬降水粒子によるミリ波後方強調散乱の計算機シミュレーション

小口 知宏 1

井原俊夫 21

2

関東学院大学工学総合研究所
関東学院大学工学部電気・電子工学科

１．はじめに衛星搭載降雨レーダの雨量測定精度向上などのために，マイクロ波だけでなくミリ波の様な高い周波数帯をも用いる計画がある[1]． この様な高い周波数になると，雨滴間の多重散乱が無視出来なくなり，特にレーダのように送受信アンテナがほぼ同一位置にある場合，後方強調散乱と言われる現象が生じている可能性がある[2,3]．多重散乱とその強調が起こっている場合に通常の雨量推定アルゴリズムを用いると，推定された
雨量に誤差を生ずる事になる．強調が本当に起こるのか，また強調の程度はどのくらいになるかを調べるために，模擬降
水粒子を用いた室内モデル実験の計画が進んでいる[4]． この計画と平行して実験モデルに対応する計算機シミュレーシ
ョンも行っている．シミュレーションに関しては，平面スラブ領
域での Rayleigh 散乱体に対するモンテカルロシミュレーションなども試みられているが［5］，ここでは出来るだけ実験の条件を正しく再現するようにモデル化し，とりあえず第一段階として最も単純な，正直な方法でどこまで問題に対処出来るかを調べてみた．本稿では今までに得られた結果と問題点について報告する．

二つの経路は存在しない．

散乱体積
受信アンテナ3
2
1

送信アンテナ

水球(a) 従来の配置

送信アンテナ
散乱体積
受信アンテナ3
2
1
水球

２．モデル実験
後方強調散乱を調べるためのモデル実験シムテムを図１に示す．散乱体は水入りポリスチロール球で，これをポリフォー
ム厚板中に埋め込んで積み重ねたものが散乱体積となる．散乱計はネットワークアナライザを用いる３０GHｚ帯のもので，直線偏波，円偏波の切り替えとポラリメトリック測定が可能である[4,6]．従来の研究結果から強調の現象は散乱角が１ﾟ以下の所で起こることが予測されるため[7]，送受信アンテナの間隔を狭くしてその存在を調べるのであるが，従来の配置(a)では散乱角を１ﾟ以下にすることが物理的に不可能であった．現在計
画中の配置(b)では，この点を改善するために受信アンテナと散乱体の間にガラス板を置き，ガラス板の反射・透過特性を利用して等価的に送受信アンテナの作る散乱角が１ﾟ以下に出来る様にしている．

2003/11/8

ガラス板

5(b) 計画中の配置

図１ モデル実験システム

１回散乱

1

２回散乱

2
1
3

３回散乱

2
ケース１

３．散乱の過程と基本式
２回以上の散乱では図２に示すように１→２，２→１の様に互いに逆周りの経路が存在する．送受信アンテナが同一地点にあると，これらの経路の波は同位相で加算されるので強め合い（強調散乱），二つのアンテナが離れると位相差が生じてこの様な加算は起こらなくなる．また，３回散乱では３の散乱体が１の散乱体そのものになる場合もあり，このときは上のような

ケース２

1
3
1

2

2003/11/10

図２ 散乱過程の模式図
なお，通常のレーダ方程式で解析の対象となるのは，一番上の経路である１回散乱のみである．

1

（ｒ）はｚ＝０の平面上任意の位置（xr, yr, ０）にあるとする．ただし，受信アンテナのビームセンターは常に散乱体積の中心を向いている．これらの記号と送受信アンテナ，および散乱体の座標ベクトルとの関係は付録 A に示す．S(1,n), S(n,1) はこの様な波の干渉を式で表すと以下のようになる．水平および垂直偏波の電界を，Eh および Ev とし， この二つを纏めて

x
R(t,1)t１

θ(t,1)⊿R(t,1)

Beam centern2

θ(n,r)r

S(1,n) = S(n)S(n−1) LS( 2)S(1)⊿R(1,2)R0(t,1)(1)

S ( n ,1) = S S

z

(n,r)
0

R(n,r)

図３

LS( n −1)

S

(n)

（6）

であり，行きと帰りの経路についてｎ個の散乱マトリクスを全て掛け合わせた物である．S の上付のバーは帰りの（時間反転）経路であることを表す．また，ｎ＝１の場合は⊿R(n,r)RR 0(n,r)

( 2)

n−1

∏{L} =1（7）  θ ( i, j ) 2
 log e 2Fe(i , j ) = exp − 2  θ b 


（8）

(

（9）

i=1

記号の説明図と置く．

E 
E =  h
Ev 

(1)

F p(i , j ) = Fe( i , j )の様に表し，コヒーレンシー・マトリクスの要素を縦並びにしたものをＤで表すと，Ｄは

D = E ⊗ E*と書ける．記号 ⊗ はクロネッカー積である． 経路１→ｎ，およびｎ→１の散乱を繰り返した後の受信電界を E(1,n)，E(n,1)と書くと，
Ｄは

(

) (

(
+ (E

(1, n )⊗ E ( n ,1)* + E ( n ,1) ⊗ E (1,n )*

T

[

)
)

Q = diag (−1,1)（3）

E (1, n ) ⊗ E (1, n )*

[

（4）

[

)

]

n −1

]

i =1

（5）

[

 exp 2 Im(ke )∆R(i ,i +1)

∏ 


(∆R

(12)

]

[

[

]

(

)

 exp − ik ( R0(t ,1) − R0(1, r ) + R0( n, r ) − R0(t , n ) ) =R(t ,1) R(1, r ) R( n, r ) R(t , n )
exp − ike (∆R(t ,1) + ∆R( n, r ) ) exp ike* (∆R(1, r ) + ∆R(t , n) )と書ける．E(n,1)についても同様に[(

E(1, n ) ⊗ E( n,1)*

 exp(−ike ∆R(i ,i +1) )  ( n, r ) (t ,1) (1, n ) 0Fe Fe S E∏∆R (i ,i +1)
i =1 

n −1

 exp − ik ( R 0(1, r ) + R 0( t , n ) ) E ( n ,1) = R (1, r ) R ( t , n )(1, r )
(t ,n )
+ ∆R
exp − ik e ( ∆ R)] exp 2 Im(k e )(∆R ( t ,1) + ∆R ( n , r ) ) 
=

2


R ( t ,1) R ( n , r )( i , i +1) n −1 
 exp 2 Im(k e )∆R∏( i , i +1) 2
i =1 
∆R( n,r )( t ,1)(1, n )
(1, n )*
F p F p (S⊗S)(E 0 ⊗ E 0* )

[

]

（11）

で与えられる対角マトリクスである．なお、本報告では常に波の進行方向を基準にとる座標表示（FSA 表示）を用いる．これ
らの関係を（３）式に適用すると]

 exp − ik ( R0( n, r ) + R0(t ,1) ) 
E =R ( n, r ) R (t ,1)( n,r )( t ,1)+ ∆R )
exp − ike (∆R

[

（10）

の様な reciprocity の関係が成り立つ．ここでＱはと書ける．（３）式最後の上下の括弧が，散乱問題で ladderterm，cyclical term と云われるものに相当する．
アンテナ出力ポートにおける等価的な電界 E0 と E(1,n)の関係
は，図３に示す各記号を用いて(1, n )2

S ( n ,1) = QS (1,n ) Q)= E (1,n ) ⊗ E (1,n )* + E ( n ,1) ⊗ E ( n ,1)`2はアンテナの電界指向性および電力指向性を表す．
一方、散乱マトリクスについては，何回かの散乱を繰り返しても，時間反転経路と元の経路に対する散乱マトリクスの間には（2）

D = E(1,n) + E(n,1) ⊗ E(1,n)* + E(n,1)*))( i , i +1) 2] []F F]( n , r ) ( t ,1) (1, r ) ( t , n )
e
e
e
e

F

F

T*

(S(1, n ) ⊗ QS(1, n ) Q)(E0 ⊗ E0* ) exp( − ik e ∆ R ( i ,i +1) )  (1, r ) ( t , n ) ( n ,1) 0
Fe Fe S E∆ R ( i , i +1)
i =1 


n −1

(13)

∏

と書ける． ただし，ｎ＝１の場合はと書ける．ｋは自由空間の伝搬定数，ke は散乱体積内部の平均伝搬定数である．送信点（ｔ）は座標の原点にあり，受信点

2５．計算結果および考察
１回から３回散乱において，散乱角と受信電力の関係を入射波と同一の偏波と交差偏波について計算したものを，図４および図５に示す．入射波は電界が単位の大きさの H 偏波とする．
記号 l, c はそれぞれ ladder term, cyclical term を表す．真後ろ（散乱角 0 ﾟ）への１回散乱の強度 Pr はn−1∏{L} =1(14)

i=1と置く．ここでの番号付け 1, 2, …,n は散乱体に固定された
番号ではなく，単に１回目，２回目，・・・，ｎ回目に散乱したことを表す物であることに注意されたい．

Pr =４．シミュレーション計算
４．１ モデル化
図３に示すようにモデル化では出来るだけ実験の配置，パラメータに合うようにした．送信波は散乱体積までは自由空間の，散乱体積内では内部の平均伝搬定数を持つ球面波で伝搬する．送信アンテナを固定し，受信アンテナの位置は可変とする．散乱の計算では各散乱回数ごとに D の加算を行い，異なった散乱回数のものとは経路長の相関が無いと考えられる
事から独立であるとしている．
４．２ 計算の手順
球体の散乱マトリクスを計算するためには，球体への入射角，散乱角を知る必要がある．ここで各球体の中心を原点とす
る局所的な座標系（ｘ’，ｙ’，ｚ’）を導入する．ｚ’は常に上向き（散乱計システムの座標では ｙ 軸方向）にとり，入射波の進行
方向単位ベクトルがｘ’ｚ’平面上に乗るようにｘ’の方向を定め
る．ｘ’，ｙ’，ｚ’を基準とする球座標での入射角α，散乱角θ，φは前後の散乱体の位置ベクトルから求められる（付録 B 参照）．同様の方法で各散乱体について入射角，散乱角を次々に求めることが出来る．
球体の場合には，入射波と散乱波の作る散乱平面上での散乱マトリクスを計算しておけば，上述の球座標を基準とした散乱マトリクスに変換することは容易である（付録 C 参照 ）．
従って前もって計算し，ファイルに保存する散乱マトリクス・データは，この球の散乱平面上の多くの散乱角について求めた物である．
次に乱数を発生，散乱体の３次元的にランダムな配置を与える．散乱体のｎ個の組み合わせについて，散乱マトリクス・データを参照しｎ回散乱を計算，これを全てのｎ個の組み合わせについて加え合わせる．ここで加算は（３）式の上下括弧中の第一項のみについて行う．第二項は逆の散乱経路の場合に計算される．ついで新しい乱数を発生，散乱体の配置を変えて上の計算を再び実行，これを３２回行う．この３２回の平均をとった物が最終データとなる．この手順を必要な散乱角について繰り返し行う．
４．３ 計算のパラメータ
今回の計算に用いたパラメータを表１に示す．このパラメー
タは以前に行った散乱実験の際のパラメータと同じ物にしている［6］．
表１

Pt Gt 0 Gr 0 λ2(D)1(4π ) 2(15)で与えられる．ここで Pt は送信電力，Gt0，Gr0 はそれぞれ送受信アンテナのビームセンターでの利得，λは波長である．計算された D の値を代入すると受信電力は-59.7dBm, これに対し従来のレーダ方程式による値は-59.6dBm でほぼ完全に一致している．

Co–polarization
SingleEhxEh*

0.0006

Double (l+c)

0.0004

Triple (l+c)

Double (l)Triple (l)Triple (c)Double (c)0.0002

0

–2

0

2

ANGLE (deg.)

図４

散乱角と受信電力の関係（同一偏波）–5[×10 ]
8

Cross polarization
Double (l+c)

EvxEv*

6

4
Double (l)Triple (l+c)

2Triple (l)Double (c)Triple (c)

0

–2

0

2

ANGLE (deg.)

図５

散乱角と受信電力の関係（交差偏波）

２回，３回散乱では散乱角０ﾟ付近で強調が起こっている事が分かるが，特に交差偏波で顕著である．２回散乱では散乱角０ﾟで ladder term と cyclical term の値は同じとなるが，３回散乱では cyclical term の値は ladder term の値に比べかな
り小さくなる．これは図２に示す様に３回散乱では１→２→１の経路があるためであるが，予想より減少の割合は大きい．強調の起こる散乱角の範囲は，波長を散乱の平均自由行程で割った程度であると云われているが[7]，今回の散乱パラメータから求めた値は約 0.8 ﾟであり，オーダとしては合致している．
４回以上の散乱は現実的な計算時間をはるかに超えるため，今回は計算出来なかった．しかし，２回散乱に比べ３回散乱がシミュレーション計算のためのパラメータ

周波数：30 GHz
散乱体積の大きさ：90cm×90cm×180cm
散乱球の大きさと数：半径 1.25cm, 1736 個
送信アンテナから散乱体積前面までの距離：1m
送信および受信アンテナのビーム幅：25 ﾟ
送受アンテナの作る散乱角の範囲：0 ﾟ～3 ﾟ
送信電力：0 dBm
送受アンテナの利得：17.2 dB

3

あまり小さくなっていない事から，４回以上の散乱の評価が多重散乱全体の大きさの定量的評価に必要な事が分かる．完
全に同じモデルではないが，今回の計算に用いたのと同じ散乱体を同じ密度で埋め込んだ厚さ 1.8m の平面スラブに平面波が入射した場合，放射伝達方程式から求めた ladder term
全体の値はほぼ１回散乱の値と同程度になる．Mishchenko によると同一偏波に対する強調の割合（enhancement factor）ζ
はς=

G + 2G
G1 + G L
1

L

表２
散乱回数

ランダム数起動

１回

３２回

２分

球の散乱マトリクスはサ

２回

３２回

１６分

球の散乱マトリクスは散乱平面上 0.001 ﾟおきに

３回

EhxEh*

１２時間４０分
球の散乱マトリクスは散乱平面上 0.1 ﾟおきに計算したデータによる．

４回

１回～916 日同上

１回 （～22 時間）スーパーコン，或いは PC
クラスターで 1000 倍の速度が得られたとした場合＊計算は CPU Pentium 4，クロック 2GHz の PC による乱までの計算結果を得た．これから散乱角１ﾟ以内で強調散乱が起こることが確認出来た．特に，交差偏波で強調の効果は顕著である．強調の起こり方は入射波が平面波か球面波かには無関係の様である．また，４回以上の多重回散乱の影響が無視出来ない程度に大きそうなことが分かった．今回の直接的方法で，現実的な時間内で４回以上の多重回散乱の計算を行うにはどうすればよいか，不正直（？）な方法による計算とその正確さの検討などが今後の課題である．

0.00040.00021020
ランダム配置回数

参考文献[1]http://www.eorc.nasda.go.jp/TRMM/gpm/index_j.htm[2] Yu. N. Barabanenkov, Yu. A. Kravtsov, V. D. Ozrin, A. I.
Saichev, “Enhanced backscattering: The universal wave
Phenomena”, Proc. IEEE, Vol.79, No.10, pp.1367-1370 (1991).
[3] 岩井俊昭，岡本卓，朝倉利光，“光散乱現象研究の展開―単一散乱から多重散乱まで―”，応用物理，第６３巻，第１号，pp．14-22(1994).
[4] 井原俊夫，小口知宏，田崎民生，金子和正，“模擬降水粒子によるミリ波後方強調多重散乱の測定法の検討”，2002
年 度 関 東 学 院 大 学 工 学 部 研 究 発 表 会 , 75, pp.149-150(2002).
[5] M.P.van Albada, A.Lagendijk, “Vector character of light inweak localization: Spatial anisotropy in coherent backscatteringfrom a random medium”, Physical Review B, Vol.36, No.4,
pp.2353-2356 (1987).
[6] T.Tazaki, H.Tabuchi, K. Ikeda, T.Oguchi, S.Ito, “Laboratory
measurements of polarimetric radar signatures of randomly
distributed spherical and spheroidal scatterers at 30 GHz”, IEE
Proc.-Mocrow. Antennas Propag., Vol.147, No.1, pp.8-12(2000).
[7] Y.Kuga, L.Tsang, A.Ishimaru, “Depolarization effects of theenhanced retroreflectance from a dense distribution of sphericalparticles”, JOSA Communications, Vol.2, No.4, pp.616-618(1985).
[8] M. I. Mishchenko, “Enhanced backscattering of polarized
light from discrete random media: calculations in exactly thebackscattering direction”, J. Opt. Soc. Am. A/ Vol.9, No.6,
pp.978-982 (1992).
[9] S. Chandrasekhar, “Radiative Transfer”, Dover (1960).

30

ランダム配置ごとの受信電力の変動（１回散乱）
２回散乱（Ladder term）

EhxEh*

0.00020.0001102030

ランダム配置回数

図７

１回

(16)

0.000600

備考

計算したデータによる．

１回散乱

図６

計算時間

ブ・ルーチンによる．

で与えられる[8]．ここで G1，GL はそれぞれ１回散乱，ladderterm の電力に比例する量である．上述の放射伝達方程式による ladder term の値が正しいとすれば，今回のモデルに対する enhancement factor は 1.5 となる．

00

計算に要する時間＊

ランダム配置ごとの受信電力の変動（２回散乱；ladder term）

図 6, 7 にランダム配置ごとの受信電力の変動を示す．
散乱回数を一定にした条件のもとでの変動であるから受信電力全体の変動とは云えないが，実験で行う 32 回の平均操作はほぼ十分であると云えそうである．表 2 に計算に要する時間を纏めた．散乱回数を増やすことは上限が散乱体数（1736）の Do loop を増やす事になるので，散乱回数が増えるとすぐに非現実的な時間が掛かるようになる．
6．まとめ
直接的な計算機シミュレーション法を用いて１回から３回散

4

r r

付録Ａ
図３の各記号と座標ベクトル等の関係を纏めておく．送信点
t から散乱体 i までの長さは自由空間部分と散乱体積内部の部分の和でR (t ,i ) = ri = R0(t ,i ) + ∆R (t ,i )

θ = cos −1 (j • k i +1 )r

z 0 ri(B-3)

rr − ri −1
ki = iri − ri −1

(B-4)

rr
r
(r − ri −1 ) − (ri − ri −1 ) • j j
li = irr(ri − ri −1 ) − (ri − ri −1 ) • j j

(B-5)

ただし

(A-1)

(A-2)

k • riz ∆R ( t ,i ) = ri 1 − 0  k•r 
i 


r

(A-3)

である．ただし， j は散乱実験の座標系でのｙ軸方向単位ベクトルである．

と書ける． k はｚ軸方向の単位ベクトル，z0 は散乱体積前面のｚ座標である．同様に散乱体 i から受信点 r までの長さもR ( r ,i ) = ri − rr = R0( r ,i ) + ∆R ( r ,i )R0( r ,i ) =z 0 ri − rr

付録 C
球体の散乱では，入射波と散乱波の作る散乱平面を偏波基準とする散乱マトリクスから球座標基準のものに変換することは容易である［9］．

(A-4)

(A-5)

k • ( ri − rr )



z0

∆R ( r ,i ) = ri − rr 1 −
 k • (r − r ) ir Scattered wave (θ,φ)∆R= ri +1 − ri

Incident wave (θ’,φ’)

z’

(A-6)

となる．また，二つの隣り合う散乱体間の長さは( i ,i +1)r

φ = cos −1 (−li • li +1 )ただし

R0( t , i ) =

α2
(A-7)

α1

θs

である．送受信アンテナのビームセンターと散乱体 i のなす角は，それぞれθ ( t ,i ) = cos −1  k •θ( r ,i )

(B-2)

ri 
ri 

y’

(A-8)


r − rr = cos  k r • i
ri − rr −1x’
図 C-1 球体の散乱平面(A-9)

図 C-1 で座標原点と入射波・散乱波の作る球の散乱平面上での散乱マトリクスを S’(θｓ)と書くと，球座標上の散乱マトリクス S(θ，φ；θ’，φ’)との関係はで与えられる．ただし， k r は受信アンテナのビームセンター
方向の単位ベクトルである．

S(θ , φ;θ ' , φ ' ) = R(α 2 + π / 2)S' (θ s )R(α1 + π / 2)
付録 B
散乱体 i の散乱マトリクスを計算するための入射角α，散乱角（θ，φ）は散乱体 i と前後の散乱体 （i－1, i+1） の位置
ベクトルを用いて以下のように書ける（図 B-1 参照）．図 B-1 の座標（ｘ’、ｙ’、ｘ’）は各散乱体に付随した局所的な座標系で
ある．
z’
α

i-1

(C-1)

と書ける．ただし，

θ s = cos −1 [cos θ ' cos θ + sin θ ' sin θ cos(φ '−φ )]

(C-2)

 1

[cosθ sin θ '− sin θ cosθ ' cos(φ '−φ )]
 sin θ s

α 1 = cos −1 

i+1




1
[cos θ ' sin θ − sin θ ' cos θ cos( φ '−φ ) ]
sin
θ
s



α 2 = cos −1 

θ

(C-3)

i

y’

であり，また R はφ cos βR(β ) = − sin βx’
図 B-1 入射角と散乱角r r

α = cos −1 (−j • k i )sin β cos β で与えられる偏波の回転マトリクスである．
(B-1)

5

(C-4)