下部電離圏における電子密度分布同定の高解像度化に関する研究
代表研究者

安 藤 芳 晃電気通信大学 大学院情報理工学研究科 准教授

1 はじめに
電離圏は電磁波技術の登場とともに発見され、これまで電磁波によるサウンディングや、ロケットなどによる in situ 観測など、古くから様々な手法によって調べられてきている[1]。しかし、電子密度が低い高度80 km 程度の D 層と呼ばれる領域は有力な観測手法がなく、VLF/LF 帯送信局の電波を用いた低周波の電波サウンディングが唯一の有効手段と言える[2]。しかしながら、VLF/LF 帯は波長が 5 - 15 km 程度と長く、電子密度の解像度が粗くなってしまうという欠点がある。
VLF/LF 帯の電波サウンディングによる D 層の電子密度同定は、地表面の多点で観測したデータから、電磁界シミュレーションと最適化によって電子密度分布を推測する手法といえる[3]。そのため、その高解像度化には、電磁界シミュレーションの高精度化と、観測点数を増やすための小型観測デバイスの開発が必要となる。
そこで本研究では VLF/LF 帯大地-電離圏導波管伝搬の数値解析手法を確立するため、そのもっとも重要な要素となる吸収境界条件として、高性能である Perfectly Matched Layer(PML)を開発する。これまでは簡易的な吸収境界条件を用いてきたが、それでは不要な反射による誤差が生じたり、Late Time Instability と呼ばれる長時間計算した際の不安定性が生じるものであった。従って、安定かつ高性能な吸収境界条件の開発が、本技術のボトルネックになっていたと言ってよい。本研究では、複素座標変換に基づく PML を球座標系の一部の領域に適用できるよう拡張し、その性能を評価する。また、小型観測デバイスとして、長時間安定して所望の周波数の観測を行える回路を開発し、JJY 送信局（40 kHz、福島県大鷹鳥谷山）の連続観測を行う。

2 高精度な VLF/LF 帯電磁波の大地-電離圏間伝搬の数値解析手法の確立
2-1 高精度解析のための吸収境界条件
地表面で放射された VLF/LF(Very low frequency/Low frequency)帯の電磁波は大地と電離圏下部の間を反射しながら伝搬してゆく。放射源として雷放電といった自然由来のものと、旧来は航行用として、現在は水中通信や標準電波などで用いられる人工起源の送信局によるものがある。VLF/LF 帯電磁波の伝搬経路の形状は球殻状であることから、モデル化においては球座標系を用いるのが精度よく解析するために必要である。
また、電離圏や大地のもつ複雑性を考慮するには有限差分時間領域(Finite Difference Time Domain; FDTD)
法[4]を用いるのが適当と言える。
VLF/LF 帯大地-電離圏導波管伝搬を数値シミュレーションする際、高度方向は電離圏と大地に囲まれており大半は反射するので問題ないが、地表面に対して平行な方向は反射体はなく、途中で解析領域を打ち切る必要が生じる。その際に、解析領域端で電磁波を反射させないよう「吸収境界条件」を設ける必要が生じる。
吸収境界条件は、解析的な一方向伝搬を想定した条件を課すものと、仮想的吸収媒質を設けるものがあり、現在では後者のうち J. P. Berenger が開発した Perfectly Matched layer (PML)[5]が主流である。この PMLはその後に複素座標変換の観点で定式化され[6]、その変換を用いて円筒座標系や球座標系の FDTD 法に適用できるものが開発された[7]。
しかし、これまで開発された球座標系 FDTD 法の PML は、いわば「全球」を取り扱うものであり、r 方向にのみ吸収できるものであった。一方で、VLF/LF 帯大地電離圏導波管伝搬では、伝搬パスを中心に解析すればよく、これまでの PML は適用できない。
そこで本研究では VLF/LF 大地-電離圏導波管伝搬を解析するために、球座標系において方向にも吸収し、いわば「球座標系で方向にも終端させた解析領域」に適用できるPML を開発する。
2-2 定式化
複素座標変換を用いて、球座標系の各座標を次の式で変数変換する。

1  ~   s ( )d    0    j0~(1)ここで、 a  r ,  ,  であり、 s  1  j  ( ) /  とする。 ,  は長さの単位を持つものではないが、それぞれ の 微 小 な 増 分 に 対 す る 位 置 の 変 化 rd , r sin  d を 座 標 変 換 し て い る と 解 釈 す る 。 ま た 、
~     ( )d  である。r についても座標変換を施すが、たとえば r d は ~r s d とせず、 rs d と近0似する。このとき、変換された周波数領域での Maxwell の方程式は以下となる。

1 H  cot 
1 1 H H  jErrs rr sin  s ~
1 1 H r 1 1  (r H  ) jE
r sin  s  ~r sr rr H  ) 1 1 H r
1 1  (~ jE
~r sr rr s 1 E cot 
1 1 EE   jH r
rs rr sin  s ~
1 1 Er 1 1  ( r E )  jH r sin  s  ~r sr rr E ) 1 1 Er
1 1  (~  jH ~r sr rr s (2)
(3)

(4)
(5)
(6)
(7)

ここで、式(2), (5)において左辺第一、二項目はの微分をした結果に対して座標変換を行っているが、第二項については方向の伝搬には影響しないので座標変換を行わない。この簡略化によって計算回数とメモリが抑制できる。
上記の式を時間領域に変換する。それぞれ関連する場を成分毎に分けることで更新式を得る。以下の電界の更新式のみを載せる。

Ecot E1 H H    r* ,  r    Er ,rtr tE
1 H r   E ,r sin  
t~~
Er
H ~
  r Er ,
r
t~
E
Er r  ~r Er  r ,ttEr  Er  Er *  Er ,E  E  Er ,~EE r   ~r E ,
r
t

E r
1 H 
   E r ,r sin  t~~H 
E~  r   r Er ,tr
E1 H r      E ,r t~Er ~
Err,
  r Er t
t


~
E
r

E  Er  E ,

r

E
t

 ~r E .

(8)(9)(10)(11)(12)(13)

磁界の更新式もほぼ同様である。この更新式を通常の有限差分によって離散化すればよい。ただし、式(8)の第二式の H  については FDTD 法で用いる Yee セル上に割りあたらないので、近傍の値を用いて平均によって求める。

2

2-3 数値計算による吸収性能の検証
数値計算を用いて開発した吸収性能を検証する。
図 1 は開発した PML を用いて、球の一部分の計算領域中央に電流源を置いた時の電界の強度分布を、時刻の経過とともに示した図である。解析領域は、内側の半径を ri = 21.2, 外側の半径を ro = 26.8 とし、

  / 4     / 4 , 0     / 8 として、140×400×200 セルで分割した。波源を r = 24,  = /2,  = /4 の位置におき、電流源を次式とする。

J  V   {t  t 0 }2 t  t0exp S2 S 2 (14)ここでV は波源における 1 セルの体積であり、S = 12t、t0 = 6S とした。PML の減衰を表すパラメタは( M  1)c0 ln R0    ( )  2 L L M(15)としている。ここで、R0 は反射係数で R0  10-6、c0 は光速、L が層数で M はのプロファイルを決定するパラメタであり、PML の反射が最少となるよう調整する。図 1 は    / 2 における面の電界強度をいくつかの時刻で濃淡で表した図である。L  10, M  3 とした。n は FDTD 法の反復回数を表している。見てわかる
とおり、中央で発生した電磁波は n = 180 ステップあたりで r 方向に垂直な PML に到達して吸収されている。
また、n = 250 ステップあたりでφ方向に垂直な PML に到達してこれも吸収されている。n = 320 ステップではすべて吸収されて、反射波は見えなくなっていることがわかる。ここから、本手法で開発した球座標系 PMLが高い性能を持っていることがわかる。

(a) n = 110(b) n = 180(c) n = 250
図 1. 開発した PML の吸収の様子(d) n = 320

図 2 は一点において観測された電界の波形である。解析領域は、内側の半径を ri = 3.2, 外側の半径を ro = 7.49と し 、 ( 12  32)    ( 12  329 ) ,0    34  と し て 、 210 × 225 × 300 セ ル で 分 割 し た 。 波 源 を(rs , s ,  s )  (5, 2 , 83  ) におき、観測点に (ro , o , o )  (4.5, 2 , 83  ) おいて観測した。パラメタは各 PML からの反射波が異なる時刻に観測されるよう調整した。層数は L = 4, 8, 16 とし、次数 M は調整した。波源から観測点に直接到達する波は約 4 ns あたりに観測されている。また、13, 20, 27, 34 ns あたりに観測されている波は、それぞれ r(内側), r(外側), , 方向に垂直な PML からの反射波である。L = 4 ではまだ反射波大きいものの、L = 8, 16 とすることで急速に反射波小さくなり、本研究で開発した PML が良い性能を持っていることがわかる。
図 3 は、開発した PML の反射誤差の入射角度依存性である。PML の層数として L = 8 および 16 層について示している。それぞれのラベルは「その方向に垂直な PML からの反射」を意味している。見てわかるとおり、
16 層 PML であれば 100dB 以上に反射が抑制されており、十分な吸収性能を持っていると言える。r 方向に垂直な PML（内側を ri, 外側を ro で示している）は、入射角度が大きくなると吸収性能が悪くなるが、これを適用する VLF/LF 大地-電離圏導波管伝搬においては、大地表面または電離圏上部となり、使用しないので大3

きな問題とはならない。

図 2. 層数を変化させたときの観測波形

図 3. PML からの反射の入射角度特性

3 小形観測デバイスの開発
VLF/LF 帯電磁波の受信は、現在の PC であればサウンドカード等で十分なサンプリングがとれるため、直接受信した電磁波をサウンドカードに入力して録音することで観測システムが構成できる。ただし、電子密度同定のためには小型観測デバイスを多数配置する必要があるため、マイコンを用いた観測機器を開発する必要がある。また、アンテナ近傍のスイッチング電源による雑音等の影響を防ぐためのフィルタが必要となる。
また、VLF/LF 帯電波の日変化パタンはこれまでよく研究されており、この日変化パタンは電子密度同定において強力な要素となる。日変化パタンを受信するには連続的に長時間観測する必要があり、安定して受信できることが要件となる。
そこで本研究では中波用バーアンテナで受信し、増幅およびフィルタ回路を用いて、試験観測として標準電波 JJY（40 kHz）を東京都調布市で観測した。サンプリング周波数 fs = 96 kHz、高速フーリエ変換(FFT)
点数 214 = 16384 点、Hann の窓関数を用いている。
図 4 に 38kHz-42kHz の周波数帯域におけるフーリエ変換したデータの時間変化をカラーマップで示す[9]。
横軸は周波数、縦軸は時刻であり、図の上側が新しい時刻となっている。受信強度が強い周波数は白い色であらわされており、40kHz の電波が受信できていることがわかる。40 kHz(JJY)は電波時計の時刻データを送信する標準電波である。0.8 秒、0.5 秒の２つの符号で送信され、時刻データが構成される。また、10 秒ごとに 0.2 秒のマーカーが現れるが、こちらも正確に 10 秒ごとに受信できていることがわかる。また、その他の周波数は 40 kHz に比べると強度が低くなっており、受信電波の同調もよくできていることが分かる。
図 5 は連続データを取得した結果であり、横軸が現地時刻、縦軸が振幅である。40kHz における信号の 1秒間の最大値を 5 分毎に平均をとったものを 2018 年 6 月のある 2 日間を表示している。同時に、39.8kHz お
よび 40.2kHz のデータもとり、それらが 40kHz より十分小さく、広帯域なノイズが混入していないことを確認している。
図から見てわかるとおり、太陽活動によって D 層が現れる昼間は低い振幅が安定して受信されており、1 dB以内の変動となっている。また、それが消失する夜間は激しい変動が現れれていることがわかる。昼夜の間に一旦振幅が増加しているが、これは日出・日没線（ターミネータ）が伝搬パス上を通過しているときに現れる特徴である。JJY 送信局（福島県おおたかどや）と調布は直線距離が 225km と短い伝搬距離であり、電離圏の影響があまり現れずに地上波と呼ばれる空間波が支配的となり、また数 dB の変化しかないのでノイズも目立つが、それでも電離圏の一般的な性質が見て取れることがわかる。

4

図 4. 38-42kHz のスペクトルの時間変化

図 5. JJY-調布の日変化パタン

4 まとめ
本研究では、VLF/LF 帯大地-電離圏間導波管伝搬を用いた電子密度同定手法の高解像度化として、電磁界シミュレーションの高精度化、および小型観測デバイスの開発を行った。
電磁界シミュレーションの高精度化として、ボトルネックとなっていた吸収境界条件を開発した。複素座標に変換することで伝搬する電磁波を反射させずに吸収する手法を、VLF/LF 帯導波管伝搬へと適用した。開発した PML は r, , 方向に垂直な PML 全てで吸収することができた。これによって、これまで球全体の領域しか適用できなかった PML を、球座標の一部の領域のみにも適用できることとなり、VLF/LF 帯導波管伝搬で安定して動作する吸収境界条件が開発できたと言える。開発した PML は 16 層程度を設けることで  100dB以上の吸収性能を持ち、反射をほぼ完全に抑制できていることが分かった。
また、開発した小型観測デバイスを用いて、JJY の電波について、東京都調布市にて連続観測を行った。
サウンドカードで録音した波形を高速フーリエ変換し、JJY の信号が受信できていることを確認した。また、40kHz の信号のみを取り出し、平均化を用いて 5 分ごとのデータとして日変化パタンを観測した。観測結果
から、短距離の伝搬路（約 225km）でありながら、昼夜の電離圏の状態が反映されているものが観測できた。
今後の課題として、高精度シミュレーションにおいては、地形などの影響の考慮が挙げられる。大地の山岳などは大きくて半波長程度の寸法を持つため、伝搬への影響は無視できない可能性がある。また観測デバ
イスとしては、たとえば西日本の JJY (60 kHz)の受信が挙げられる。小型化をしているため、最大のサンプリング周波数が 96 kHz となっていることから、直接の信号処理ができない。そこで、ミキサ回路を加えて周波数を落として観測する必要があり、その回路の開発が今後の課題となる。

【参考文献】[1] K. Davies, Ionosphere Radio, Peter Peregrinus Ltd., London, UK, 1990.
[2] G. Bainbridge and U. S. Inan, “Ionospheric D region electron density profiles derived from the measured
interference pattern of VLF waveguide modes,” Radio Sci., vol. 38, no. 4, doi: 10.1029/2002RS002686, 2003.
[3] 安藤, 伊藤, “VLF 帯大地-電離圏導波管伝搬を用いた電子密度同定問題における適切なパラメタ選択に関する検討,” 電気学会論文誌 A, vol. 136, no. 5, pp. 304-309, 2016.
[4] A. Taflove and S. Hagness, Computational Electrodynamics: the finite-difference time-domain method, 2nd Ed.,
Artech House, MA, 2000.

5

[5] J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phys., vol.
114, pp. 185-200, 1994.
[6] W. C. Chew, J. M. Jin, and E. Michielssen, “Complex coordinate stretching as a generalized absorbing
boundary condition,” Micro. Opt. Tech. Lett., vol. 15, no. 6, pp. 363-369.
[7] F. L. Teixeira and W. C. Chew, “PML-FDTD in cylindrical and spherical grids,” IEEE Micro. Guided Wave
Lett., vol. 7, no. 9, pp. 285-287, 1997.
[8] F. L. Teixeira and W. C. Chew, “Systematic derivation of anisotropic PML absorbing media in cylindrical andspherical coordinates,” IEEE Micro. Guided Wave Lett., vol. 7, no. 11, pp. 371-373, 1997.
[9] Spectrum Lab, http://www.qsl.net/dl4yhf/spectra1.html〈発題

名An Approximate PML Applied toCylindrical and Spherical Coordinate
Sectors
VLF/LF 帯 電 波 の 多 点 観 測 及 び 3 次 元
FDTD 法を用いた金環日食に伴う電波伝搬特性の解明に関する研究

表資料〉

掲載誌・学会名等IEEE Microwave andComponent Letteres

発表年月
Wireless

Japan Geoscience Union Meeting
20186

採録決定

2018 年 5 月