大地－電離層諧振腔中極低頻電磁波傳播的數(shù)值模擬

譚仕文

（92330部隊(duì) 山東 青島 266102）

極低頻電磁波在地球物理勘探、地震預(yù)報等領(lǐng)域具有重要的用途。由于頻率很低，波長很長，極低頻電磁波的激勵與傳播可以看作是大地-電離層諧振腔中電磁場的空間分布問題[1]。這是一個典型的球面分層介質(zhì)中的邊值問題[2]，分析方法有近似解析方法和數(shù)值方法。目前國外學(xué)者主要采用解析方法來求解該問題，而解析方法在具體的處理技巧上又有所不同。其中一種方法是利用波長遠(yuǎn)大于電離層高度的特殊條件，得到大地-電離層腔體中齊次波動方程的近似解，然后基于這種近似的本征模獲得天線激勵的完備解，這種方法主要被美國學(xué)者所采用。而俄羅斯學(xué)者則是從經(jīng)典的球諧級數(shù)解和散射迭加原理出發(fā)，運(yùn)用Watson變換，將天線輻射場的球諧級數(shù)解變換到復(fù)平面上的圍線積分，然后從中分離出橫向的特征模式解。這兩種處理方法本質(zhì)上是等效的，可以通稱為模式理論。由球面波導(dǎo)模式理論得到的計(jì)算公式中，含有待定本征值的特殊函數(shù)，只有在一定條件下，才能導(dǎo)出本征值及特殊函數(shù)的近似表達(dá)式。

除了解析方法之外，作為一種典型的電磁場數(shù)值算法，時域有限差分法（FDTD）[3-5]也開始應(yīng)用極低頻電磁場的模擬中。FDTD以差分原理為基礎(chǔ)，直接從Maxwell方程出發(fā)，將其轉(zhuǎn)換為差分方程組，在一定體積內(nèi)和一段時間上對連續(xù)電磁場的數(shù)據(jù)取樣。因此，它是對電磁場問題的最原始、最本質(zhì)、最完備的數(shù)值模擬，具有最廣泛的適用性。FDTD算法簡單直觀、容易掌握。從Maxwell出發(fā)，不需任何導(dǎo)出方程，避免了使用更多的數(shù)學(xué)工具，使得它成為所有電磁場計(jì)算方法中最簡單的一種。其次，它可以在計(jì)算機(jī)所能提供的離散數(shù)值時空中仿真再現(xiàn)電磁現(xiàn)象的物理過程，非常直觀，易于可視化顯示。

本文基于積分形式的Maxwell方程，推導(dǎo)了大地—電離層諧振腔中的FDTD迭代方程，計(jì)算了極低頻電磁波的輻射和傳播過程，得到了可視化的結(jié)果。

1 FDTD的剖分和時間步進(jìn)迭代公式

對于ELF頻段，電離層高度遠(yuǎn)小于波長，因此可近似認(rèn)為場沿高度方向沒有變化，這樣可以簡化為二維問題來處理[6-8]。對于地球表面。用經(jīng)緯線進(jìn)行剖分是比較方便的，如圖1所示?？梢娺@種剖分得到的離散網(wǎng)格是非均勻，越靠近南北兩極，沿緯線方向的網(wǎng)格長度越小。

設(shè)沿經(jīng)線方向的剖分單元數(shù)為M，沿緯線方向的剖分單元數(shù)為2×M，沿緯線方向的網(wǎng)格編號為i，沿經(jīng)線方向的網(wǎng)格編號為j，網(wǎng)格編號的參考點(diǎn)位于網(wǎng)格中心。沿緯線方向的網(wǎng)格長度為Δw-e，沿經(jīng)線方向的網(wǎng)格長度為Δs-n，網(wǎng)格面積為S。

圖1 地球表面剖分網(wǎng)格Fig.1 Earth’s surface mesh

在各向同性無源空間中，Maxwell旋度方程為：

安培定律式（2）中的面積分和環(huán)路積分可近似為：

從而可以得到Ez的差分方程為

電磁感應(yīng)定律式（1）中的面積分和環(huán)路積分可以近似寫作：

其中h為電離層高度，于是可以得到Hx的差分方程為：

同樣可以得到Hy的差分方程為

在南北兩極上，四邊形網(wǎng)格退化為三角形，此時差分方程變?yōu)椋?/p>

當(dāng)存在電流源時，安培定律變?yōu)椋?/p>

其差分近似形式為：

在靠近南北兩極時，沿經(jīng)度方向與沿緯度方向的網(wǎng)格長度差距將越來越大，將會影響到數(shù)值計(jì)算的穩(wěn)定性和計(jì)算效率，因此實(shí)際編程中需仔細(xì)處理，合理選擇時間步進(jìn)。

2 計(jì)算示例

2.1 簡諧振蕩電流源在地－電離層波導(dǎo)中激勵的場

假設(shè)在武漢地區(qū)（約北緯30度，東經(jīng)110度）設(shè)置一個簡諧振蕩的垂直電流源，頻率為100 Hz。計(jì)算了不同時刻球面上電磁場的分布狀況，如圖3～5所示。在二維平面示意圖中，下邊緣對應(yīng)北極，上邊緣對應(yīng)南極，左邊緣對應(yīng)經(jīng)度的0度。時間步進(jìn)為 5 μs。

圖2 一個周期時的球面電磁場分布示意圖Fig.2 A cycle of spherical electromagnetic field distribution diagram

圖3 3個周期時的球面電磁場分布示意圖Fig.3 Spherical electromagnetic field distribution diagram of three cycles

圖4 8個周期時的球面電磁場分布示意圖Fig.4 Spherical electromagnetic field distribution diagram of eight cycles

圖2為1個周期時的場分布，可見電磁場以源點(diǎn)為中心，均勻地向四周擴(kuò)散開來。

圖3為3個周期時的電磁場分布圖，此時電磁波已越過北極地區(qū)，可見時間步進(jìn)迭代過程是穩(wěn)定的。

圖4為8個周期時的電磁場分布圖，此時，電磁波已到達(dá)對極點(diǎn)，可以看出經(jīng)球面?zhèn)鞑ズ?，電磁波匯聚到對極點(diǎn)，從而在對極點(diǎn)附近形成一個強(qiáng)信號區(qū)域。

2.2 高斯脈沖電流源在地－電離層波導(dǎo)中激勵的場

雷電等自然極低頻信號源可用高斯脈沖來模擬，將上例中的簡諧源換成高斯脈沖，重新計(jì)算了電磁場的分布。圖5為迭代了50 000步后電磁場的分布圖。此時激勵源信號早已消失，但地-電離層腔體中仍有穩(wěn)定的電磁波存在，這表明地-電離層波導(dǎo)中存在諧振現(xiàn)象。圖6為對極點(diǎn)電場強(qiáng)度隨時間的變化曲線，可以看出兩次正的峰值之間的時間間隔約為0.13 s，由此可以估算出地-電離層腔體的諧振頻率約為7.7 Hz。這個數(shù)值與有關(guān)文獻(xiàn)給出的地球舒曼（Schumem）諧振的基波諧振頻率十分接近。

圖5 高斯脈沖信號激勵的球面電磁場分布示意圖Fig.5 Spherical electromagnetic field distribution diagram with Gauss pulse signal excitation

圖6 對極點(diǎn)電場的時域波形Fig.6 Electric field waveform in time domain on pole

3 結(jié)束語

文中采用FDTD方法計(jì)算了ELF電磁場在地-電離層波導(dǎo)中的分布，并給出了可視化的結(jié)果。從計(jì)算結(jié)果來看，大地-電離層腔體中的ELF電磁場衰減很小，而且球面波導(dǎo)結(jié)構(gòu)可以將電磁波匯聚到天線的對極點(diǎn)，使得對極點(diǎn)附近成為強(qiáng)電磁場區(qū)域。從高斯脈沖源的FDTD計(jì)算結(jié)果來看，大地-電離層腔體在極低頻段存在諧振點(diǎn)，可以觀察到明顯的舒曼（Schumem）諧振效應(yīng)。

由于缺乏全球的大地及電離層參數(shù)，在本文的計(jì)算中，暫時沒有計(jì)入大地和電離層電導(dǎo)率的有限性和非均勻性。但FDTD是一種基于差分方程的算法，適合于處理媒質(zhì)的非均勻性問題。只要建立了各個地區(qū)的大地電導(dǎo)率及電離層電子濃度數(shù)據(jù)庫，應(yīng)用該方法就可以模擬出極低頻電磁場在全球的傳播與分布。

[1]潘威炎.長波、超長波、極長波傳播[M].成都：電子科技大學(xué)出版社，2004.

[2]Wait J R.Electromagnetic wave in stratified media[M].London:Pergamon press，1962.

[3]王秉中.計(jì)算電磁學(xué)[M].北京:科學(xué)出版社，2002.

[4]王長華，王秩雄.基于時域有限差分法的諧振腔分析[J].現(xiàn)代電子技術(shù)，2011（2）：103-105，108.

WANG Chang-hua，WANG Zhi-xiong.Analysis of cavity resonator based on finite difference time domain[J].Modern Electronics Technique，2011（2）：103-105，108

[5]肖健，高愛華.飛秒激光觸發(fā)光電導(dǎo)天線產(chǎn)生太赫茲波的研究[J].電子科技，2010（3）：7-9.

XIAO Jian，GAO Ai-hua.Investigation of generation of terahertz with the photoconductor triggered by the femosecond laser pulse[J].Electronic Science and Technology，2010（3）：7-9.

[6]Burrows M L.ELF Communications antennas[J].Peter Peregrinus Ltd.1978.

[7]Simpson J J，Taflove A.Two-dimensional FDTD model of antipodal ELF propagation and Schumann resonance of the earth[J].IEEE Antennas and Wireless Propagation Letters，2002（1）：53-56.

[8]Simpson J J，Taflove A.Three-dimensional FDTD modeling of impulsive ELF propagation about the earth-sphere[J].IEEE TransactionsonAntennasandPropagation，2004（52）：443-451.