GPR 偏波解析による月地下媒質不均一性の計測―

シミュレーションによる手法の検討

小林 敬生（独）科学技術振興機構

―

佐藤 源之
東北アジア研究センター

1. はじめに月の表面はレゴリスと呼ばれる粒径～10-6mオーダーのパウダー状の物質に覆われている。レゴリス層の厚さは数mから数十mで、深さとともに粒径が増し、粒径の大きな礫状物質で構成されるスーパーレゴリス層に移行して基岩盤につながっていると考えられている[1]。しかし、この月地下構造モデルを実証するデータはまだない。
地質学的にミクロなこの月地下構造情報は月面地形の進化および風化の歴史を理解するうえで重要である。この情報を得るためには、表面で掘削調査をするのが理想的だが費用対効果の観点からは疑問が残る。この問題を軽減するために、多地点における観測が可能で、かつ必要な地下構造情報を得る手段として我々は月着陸機に搭載できる地中レーダ（Ground Penetrating Radar : GPR）の検討を始めた。
地中レーダはその名のとおり地中に向けて電磁波パルスを送信し、反射波の受信時間を計測することにより地中の電波反射体の位置を決定するレーダである。その基本原理は雨雲の位置を調べる降雨レーダや飛行機の位置を調べる航空レーダのものと変わるところはない。我々が提唱する手法は、GPR 受信波の偏波情報を解析することにより地下媒質（レゴリス～スーパーレゴリス層）の不均一性の空間スケール（粒径）の変化を推定するものである。周波数が数百 MHz（波長 1m）から 10GHz（波長数 cm）程度まで、帯域約 10GHz の超広帯域で計測し、電波が波長と同スケールの物体によって強く散乱される性質を利用して地下に分布するスーパーレゴリス礫の粒径の変化を検出する。検討に当たって、まず、GPR 受信波の偏波比が深度方向の礫粒径変化を反映して期待するような時間変化をするのかどうかを 3 次元 FDTD[2]シミュレーションによって確認することにした。本稿は、そのシミュレーションの報告である。

2. モデル2.1. シミュレーション空間
700×700×2200 メッシュのシミュレーション空間の中に、送受信アンテナと月面（および地下媒質）をモデル化した。メッシュ寸法は 1.5×1.5×1.5mm である。モデルの実空間スケールはおよそ 1m×1m×3m である。今後行なう予定の実験室実験の評価・解釈に役立てるよう、実空間における計測をそのスケールを保存したまま、シミュレーション空間に写した。シミュレーション空間全体は 10 層厚の PML 吸収境界[3]に囲まれている。一度の計算に必要なメモリ容量は、電場、磁場のベクトル３成分と媒質の電磁特性をメッシュ毎に定義する係数すべてを含めて約 230GB である。理想を述べれば、水平方向も鉛直方向と同じ程度の大きさ、すなわち 3m×3m の広がりがほしいところである（次回の計算機システム更新に期待をしているところである）。

2.2. アンテナ
超広帯域の計測を行なうためにビバルディアンテナ[4]を使う。本稿でモデル化したビバ
ルディアンテナを図１に示す。アンテナは厚さ 1.5mm の誘電体板（比誘電率 4.5）と完全導体箔で構成されている。図の黒い部分は表面の金属箔アンテナエレメント、灰色の部分は裏面のアンテナエレメントを示す。アンテナ
給電はギャップ給電モデル[5]でモデル化し給電線と信号源は Maloney[6]の 1 次元給電線モデルを適用した。図中、給電点は上端辺中点に位置する。このアンテナ４枚を使って偏波面が直交する二組の送受信アンテナ対（図２）を構成し、地中散乱波の平行偏波と交差偏波を受信する。
図３にモデル化されたアンテナの S11 計測（のシミュレーション）をして得られたリターンロスを示す。アンテナのリターンロスはアンテナ入力電力のうち電磁波として放射されずに送信機側に戻ってきた電力の大きさとして定義される。図３をみると、500MHz よりも低い周波数では入力電力がほとんど放射されずに送信機側に戻ってくる（0dB）が、2GHzよりも高い周波数では入力電力の 90％以上（リターンロス-10dB 以下）がアンテナから放射されており、このアンテナが広い帯域で効率よく働くことがわかる。
図４と図５、図６にアンテナ端より 1.5m 離れた完全導体平面からの反射波測定によるアンテナの送信パルス波形とそのスペクトルを示す。図４は平行偏波成分（E 面成分）のパルス波形、図５は交差偏波成分（H 面成分）のパルス波形である。図６より、送信波の偏波比は約 20dB であることが分かり、このアンテナによる有意な偏波比測定のダイナミックレンジは 20dB であることがわかる。
送受信アンテナ（図２）は、アンテナ下端が月表面平均面より 9cm 上に位置するように設置した。レーダーパルスを地面に向けて発射し地中からの散乱波を受信する。

100mm100mm
図１．ビバルディアンテナ

Y 受信アンテナ

X 受信アンテナ

Y 送信アンテナ
X 送信アンテナ
図２．送受信アンテナ

図３．アンテナリターンロス2.3. 地下媒質
地下媒質モデルとして３種類を検討した。
まず、月地下は月物質が稠密に存在するのではなく、体積占有率が 30％の比較的疎な状態であると仮定した。３つのモデルの違いは、平均粒径分布にある。図７にそれぞれのモデ
ル月面の垂直断面図を示す。黒い部分は月物質が充填されている部分で、白く抜けている部分は真空である。
モデル１では、月物質の平均粒径は深さによらず 20mm で一定であると仮定した。モデル２では、月物質の平均粒径が 20mm から 300mまで深さの 1 次関数として変化すると仮定し、モデル３では同様に 40mm から 600mm（モデル２の 2 倍の大きさ）まで変化すると仮定した。
すべてのモデルにおいて粒形は 3 軸不等楕円体で、3 軸の軸長は Monte Carlo 法でランダムに与える。軸長の分布は正規分布でその標準偏差はそれぞれの深さでの平均値の 30％の値としている。また、月物質の誘電率は 1.3（深度０）から 5.3（深度 3m）まで、月表面からの深さの 1 次関数として変化すると仮定した。導電率は 0 としている。

図４．送信パルス波形（平行偏波成分）

図５．送信パルス波形（交差偏波成分）

3. シミュレーションおよび解析
3.1. 計測
シミュレーションでは、この GPR システムが月面を移動可能なプラットフォームに装荷されたと仮定し、月面上を移動しながら 30cmの間隔をおいて計測を行い直線軌道上の 10 地点で計測を行なった。
各計測点では、送信アンテナを Gaussian パルスで励振して平行偏波、交差偏波の受信計測を同時に行なう。Gaussian パルスで励振さ
図６．送信パルススペクトル．
れたアンテナは図４（アンテナ面平行成分）と図５（アンテナ面垂直成分）に示す波形の Co-Pol は平行偏波 Cx-Pol は交差偏波を示す．
成分を含む電磁波パルスを送信する。
図６に示すように、このパルスは非常に広い帯域の成分を含んでいるため（時間的に）大変短いパルスとなっている。
2 組の送受信アンテナが直交するように設置されるので、2 組の独立した偏波受信データが得られる。便宜的にこれら 2 組の送受信アンテナを X アンテナ、Y アンテナと呼ぶことにすれば、平行偏波成分としては XX と YY のアンテナの組による計測、交差偏波成分としては XY、YX のアンテナの組による計測の４個の受信データが１計測点で得られることになる。
得られた受信波データは時間長 6.144nsec（4096 点）の時間窓を連続的に動かしながら短時間フーリエ解析（STFT）を施し、１ＧＨz から９ＧＨz までの 50 の周波数レンジの偏波比の時間変化を調べる。時々刻々の偏波比変化はバラつきが非常に大きいので、本稿では各偏波成分のパワーの観測値平均の比をとって偏波比の代表値とした。
すなわち、周波数 f、時刻 t における偏波比 R(f,t)を∑ ⎡⎢ PXY , n ( f , t ) + PYX , n ( f , t )⎤⎥⎣⎦R( f , t ) = n（１）∑ ⎡⎢ PXX , n ( f , t ) + PYY , n ( f , t )⎤⎥⎦n⎣のように定義する。ここで、P の添字の X,Y は送受信アンテナの組を示す。たとえば、添字が YX の場合は送信アンテナが Y アンテナ、受信アンテナは X アンテナであることを示す。n は観測の番号で、１から 10 までである。

（１）モデル１
3.2. 受信波波形（２）モデル２
図７．地下媒質構造モデル垂直断面図．

（３）モデル３

図７にモデル２の計測によって得られたデータの１例を示す。図 8(a)は X アンテナ送信、X アンテナ受信による平行偏波の受信波形、図(a)７(b)は X アンテナ送信、Y アンテナ受信による交差偏波の受信波形である。両者を見比べると、両者とも 3nsec あたりで受信が始まる。
月表面付近の反射波はパルシブな波形で広帯域であることをうかがわせるが、時間の経過とともに波形の山と谷の区別が容易になり、高周波成分が減衰していることがわかる。さらに、表面近傍での反射波振幅は平行偏波成 (b)分のそれが交差偏波成分に比べ非常に大きいが時間の経過とともにその比は小さくなり、
25nsec 以降（1.5m 以深からの反射）ではもはやその違いは明らかではない。これは送信波が地下媒質の不均一性により散乱を受けてそのエネルギーの一部が交差偏波成分のエネルギーに変換されている様子を示している。

3.3. 偏波比スペクトル時間変化
図９に各モデルの計測から得られた偏波比スペクトルの時間変化を示す。偏波スペクト
ル強度の比をカラーコードで表している。時刻は平均月面からの反射波の受信開始時刻を基準に示している。
図中描かれた白線は偏波比-15dB のおよその境界線（等高線）である。図６から分かる
ように、-15dB の偏波比は送信パルスの偏波比が保存されている状態を示す。図 9 では、この-15dB 線より左側、つまり受信開始後の早い時間には偏波に変化はなく、地下深くへ進むにつれて高周波（短波長）成分から散乱の影響が出始め、順次低周波側に変化が伝搬することが容易にわかる。これは、直感的にも容易に予想されることである。
この図９を地下媒質構造モデル図（図７）と比較することにより、散乱の影響の現れ方が地下媒質の不均一性の分布を反映していることも容易にうかがえる。たとえば、粒径の小さな粒が比較的均一に存在する場合（モデル１）には、散乱の影響を受け始めるまでにかかる時間が他の２つのモデルの場合にくらべ長い点があげられる。

図８．受信波波形．
(a)平行偏波成分と(b)交差偏波成分

図９．偏波比スペクトル比時間変化．
上から、モデル１、モデル２、モデル３での計測結果を示す。

4. おわりに超広帯域 GPR の偏波情報を利用した地下媒質不均一性の計測のシミュレーションをFDTD 法により行なった。送受信アンテナと地下媒質を同時にモデル化し、実寸法による大規模なシミュレーションを行なった。
シミュレーションによって、地下媒質散乱波の偏波比スペクトルの時間変化を示すことができた。本稿のモデルでは、計算時間を節約することと現象の変化を強調することを目的として真空・誘電体の２媒質のみで地下媒質をモデル化した。
地下の不均一媒質をモデル化するに当たり、本稿では滑らかな表面を持ち、一様な内部構造をもつ３軸不等楕円体を地下空間に分布させたが、これは現実を著しく単純化したモデルである。この単純モデルの有効性については、これから予定されている実験による検証で確認する必要がある。モデルで今後改善すべき点はまず、地下媒質誘電体の滑らかな表面を粗くすること、そして本稿では真空として単純化した地下空間を誘電体で埋めることである。
本稿では、偏波情報の解析を偏波比スペクトルの時間変化を調べることで行なった。
合成開口レーダ（SAR）の世界では偏波情報の解析を行なうことで計測ターゲットの識別を試みる手法が確立されている[7]。今後はこうしたリモートセンシングにおける偏波解析の手法を応用して地下の構造を推定する手法の研究を進める必要がある。
FDTD 法の主要部は前年度までに完成させていたので、本年度の課題は解くべき問題のモデル化の部分であった。地下不均一媒質のモデル化は Monte Carlo 法を適用してシミュレーション空間の中に誘電体媒質をランダムに充填する方法を取った。計算時間を現実的な時間に抑えるために誘電体の占める空間体積比率を 30％に抑えたが、それでもモデル１のケースでは、予想を大きく超える時間がかかってしまった。体積比率を確認しながらランダムに順次誘電体を充填する方法はベクトル化にも並列化にも向かず計算時間を長くする要因である。モデル３では、礫が大きいおかげで個数が少なくなり計算全体としての効率は比較的よいものであった。32CPU による計算で、並列化率は 82.5%で、ベクトル化率は 99.7％、総計算時間は 396308 秒であった。しかし、モデル１の場合には並列化率が 20%まで低下した。乱数発生により礫粒子を１個ずつ発生させていたためである。ベクトル化率は 99.0%で、総計算時間は 438715秒であった。２CPU 以上の並列化率は 85%なので、本シミュレーションにおける律速段階はこの礫生成にあることがわかる。
大規模な計算を行なうときは、モニター用のデータ出力を随時行なうようにしておかないと、計算の進み具合は一般ユーザにはわからない。しかし、データの出力は計算の効率に大きな影響を与える。今回の計算では情報シナジーセンターのジョブモニ
ターのおかげで計算効率を著しく阻害していた小さなバグを早期に検出することができた。この場を借りて感謝したい。

謝辞本研究におけるシミュレーションコードの開発は東北大学情報シナジーセンターとの共同開発で行なわれた。また本研究の一部は科学研究費補助金基盤研究(S)14102024 によるものである。

参考文献[1] Heiken, G., D. Vaniman, and B. M. French, Lunar Source Book, Cambridge UniversityPress, Cambridge,[2] Yee, K. S., Numerical solution of initial boundary value problems involving Maxwell’sequations in isotropic media, IEEE Trans. Antennas and Propagation, Vol.14, p.302, 1966
[3] Berenger, J. P., Three-dimensional perfectly matched layer for the absorption ofelectromagnetic waves, J. Computational Physics, Vol.127, p.363, 1996
[4] Langley, L. D., P. S. Hall, and P. Newham, Balanced antipodal Vivaldi antenna for widebandwidth phased arrays, IEE Proc. Microwave, Antennas and Propagation, Vol.143, p.97,1996[5] Watanabe, S. and M. Taki, An improved FDTD model for the feeding gap of a thin wire
antenna, IEEE Trans. of Microwave Theory and Techniques, Vol.46, p.2131, 1998
[6] Maloney , J. G., K. K. Shlager, and G. S. Smith, A simple FDTD model for transient
excitation of antennas by transmission lines, IEEE Trans. Antennas and Propagation,Vol.42, p. 289, 1994,[7] Cloude, S. R. and E. Pottier, A Review of Target Decomposition Theorems in RadarPolarimetry, IEEE Trans. on Geoscience. and Remote Sensing, Vol. 34, p.498 1996.