本發(fā)明涉及雷電定位與電磁場(chǎng)計(jì)算領(lǐng)域，具體的指一種使用Matlab語(yǔ)言計(jì)算衰減因子，與理想場(chǎng)進(jìn)行卷積求得衰減電磁場(chǎng)的方法。

背景技術(shù)：

雷電產(chǎn)生的甚低頻(VLF,3-30kHz)和低頻(LF,30-300kHz)信號(hào)可以傳播很遠(yuǎn)，在數(shù)百公里到上千公里處都可以觀測(cè)到。地波沿著地球表面進(jìn)行傳播，在傳播過(guò)程中主要受地表電導(dǎo)率的影響，發(fā)生衰減與色散。同時(shí)由于傳播距離遠(yuǎn)，要考慮地球曲率造成的影響，會(huì)使觀測(cè)結(jié)果相較于沿平直方向傳播不同。這些影響因子會(huì)干擾觀測(cè)者根據(jù)觀測(cè)數(shù)據(jù)判斷閃電類型、反演電流矩、進(jìn)行地基閃電定位等。

雷電產(chǎn)生的電磁場(chǎng)在傳播過(guò)程中包含水平電場(chǎng)、垂直電場(chǎng)和水平磁場(chǎng)三部分。Delfino^[1],[2]研究了Cooray-Rubinstein(C-R)算法的有效性，此算法用來(lái)計(jì)算雷電產(chǎn)生的水平電場(chǎng)，Delfino指出C-R算法有一個(gè)理想的計(jì)算精度。水平電場(chǎng)包含大量高頻成分，在傳播距離很遠(yuǎn)時(shí)幾乎全部衰減掉，遠(yuǎn)距離觀測(cè)主要采集垂直電場(chǎng)信號(hào)和水平磁場(chǎng)信號(hào)。對(duì)于垂直電場(chǎng)的研究，Cooray^[3],[4]等利用Sommerfeld積分的精確解驗(yàn)證了近似算法的準(zhǔn)確度。Shooray^[5]等利用時(shí)域有限差分法(FDTD)研究了Wait算法的準(zhǔn)確性，并指出在距離閃電通道10-100km范圍內(nèi)，Wait算法計(jì)算出的垂直電場(chǎng)峰值與波形與FDTD的計(jì)算結(jié)果有很好的對(duì)應(yīng)。在信號(hào)傳播距離達(dá)到數(shù)千公里時(shí)，長(zhǎng)波長(zhǎng)的VLF波導(dǎo)傳播模型已被成功開發(fā)應(yīng)用，此模型更適用于窄頻段的信標(biāo)信號(hào)。對(duì)于自然閃電產(chǎn)生的寬帶信號(hào)而言，這種模型濾掉了過(guò)多波形特征，導(dǎo)致結(jié)果與回?fù)舢a(chǎn)生的信號(hào)不同，使接收到的信號(hào)更像一個(gè)脈沖信號(hào)而不是閃電信號(hào)。Cummer^[6]和Hu^[7]提出了一種改進(jìn)的二維柱坐標(biāo)FDTD模型，此模型可用來(lái)計(jì)算地面-電離層電磁場(chǎng)傳播，截止頻率為30kHz。這種算法的結(jié)果與Pappert和Feguson^[8]的波導(dǎo)模型有很好的一致性，與ELF-VLF的閃電觀測(cè)結(jié)果相對(duì)應(yīng)。但相較于Wait算法，此模型計(jì)算時(shí)間較長(zhǎng)，且對(duì)計(jì)算頻率有限制，計(jì)算過(guò)程中存在色散問題。

本發(fā)明主要關(guān)注閃電產(chǎn)生的電磁場(chǎng)信號(hào)傳播至遠(yuǎn)距離處(數(shù)百至上千公里)的特征，這些信號(hào)對(duì)于全球閃電定位、電流矩與電荷矩的反演起到至關(guān)重要的作用。盡管Wait算法^{[9],[10],[11]}已經(jīng)指出遠(yuǎn)距離電磁場(chǎng)傳播的計(jì)算方法，但此算法在計(jì)算過(guò)程中過(guò)于復(fù)雜，本發(fā)明試圖利用一種近似算法對(duì)遠(yuǎn)距離電磁場(chǎng)傳播進(jìn)行計(jì)算，以簡(jiǎn)化計(jì)算過(guò)程。

本發(fā)明的技術(shù)解決方案步驟如下：

步驟1，計(jì)算理想條件(地球曲率半徑無(wú)窮大，電導(dǎo)率無(wú)窮大)下傳播一定距離(數(shù)百千米)的電場(chǎng)E₀和磁場(chǎng)B₀，進(jìn)行傅立葉變換得到頻域場(chǎng)E_w和B_w；

E_w＝fft(E₀)，B_w＝fft (1)

步驟2，根據(jù)W.Sollfery^[12]提出的近似算法，計(jì)算曲率衰減因子W₁：

ts可表示為：

t_s＝e^-iπ/3τ_s＝e^-iπ/3(3πv_s/2)^2/3

步驟3，根據(jù)Wait算法，計(jì)算電導(dǎo)率衰減因子W₂：

Δ為土壤等效表面阻抗，表達(dá)式為

步驟4，將計(jì)算得到的衰減因子W₁、W₂分別與理想頻域場(chǎng)E_w、B_w進(jìn)行卷積得到頻域衰減場(chǎng)E_G、B_G，進(jìn)行傅立葉逆變換即得到時(shí)域的衰減電場(chǎng)E_t與磁場(chǎng)B_t。

E_G＝E_w×W₁×W₂，B_G＝B_w×W₁×W₂ (7)

E_t＝ifft(E_G)，B_t＝ifft(B_G) (8)

本發(fā)明的優(yōu)點(diǎn)是：

1.相比于傳統(tǒng)Wait算法，本發(fā)明提出的近似算法簡(jiǎn)化了計(jì)算過(guò)程，同時(shí)又可以保證計(jì)算結(jié)果在一定范圍內(nèi)的準(zhǔn)確性，為遠(yuǎn)距離電磁場(chǎng)傳播的衰減提供了新的計(jì)算思路。

相比于時(shí)域有限差分法(FDTD)，利用本算法可以節(jié)省計(jì)算時(shí)間，同時(shí)兩種方法得到的結(jié)果相近，優(yōu)化了傳統(tǒng)Wait算法。

附圖說(shuō)明

圖1本發(fā)明的基于簡(jiǎn)化遠(yuǎn)距離電磁場(chǎng)計(jì)算的方法總體流程圖；

圖2-1為300km處理想電場(chǎng)波形圖；

圖2-2為300km處理想磁場(chǎng)波形圖；

圖3-1為300km處衰減電場(chǎng)波形圖；

圖3-2為300km處衰減磁場(chǎng)波形圖。

具體實(shí)施方式

本發(fā)明為基于簡(jiǎn)化遠(yuǎn)距離電磁場(chǎng)計(jì)算的方法。其特征是：首先計(jì)算出理想條件(地球曲率半徑無(wú)窮大，電導(dǎo)率無(wú)窮大)下傳播一定距離的電場(chǎng)和磁場(chǎng)，并進(jìn)行傅立葉變換得到頻域場(chǎng)；分別計(jì)算電導(dǎo)率衰減因子W₁與地球曲率半徑衰減因子W₂，將兩因子與理想場(chǎng)進(jìn)行卷積，得到頻域衰減的電場(chǎng)和磁場(chǎng)；進(jìn)行傅立葉逆變換即得到時(shí)域衰減電場(chǎng)和磁場(chǎng)。由于Cummer提出的時(shí)域有限差分法(FDTD)已在其論文中有了充分論述證明，本專利利用這種方法對(duì)計(jì)算結(jié)果進(jìn)行對(duì)比修正，使結(jié)果更加準(zhǔn)確可信，從而達(dá)到簡(jiǎn)化Wait算法，更加便捷計(jì)算遠(yuǎn)距離電磁場(chǎng)傳播的目的。流程圖如圖1所示。

(1)首先計(jì)算理想電場(chǎng)與磁場(chǎng)，閃電源采用電流矩源，這種電流矩便于加入到FDTD模型中進(jìn)行計(jì)算，表達(dá)式為

式中I₀＝20kA，ν₀＝8×10⁷m/s，γ＝3×10⁴/s，a＝2×10⁴/s,b＝2×10⁵/s。

利用此電流矩源計(jì)算出傳播至一定距離處的理想電場(chǎng)與理想磁場(chǎng)，進(jìn)行傅立葉變換得到頻域的理想場(chǎng)E₀和B₀：

E_w＝fft(E₀)，B_w＝fft(B₀) (2)

(2)根據(jù)W.Sollfery提出的近似算法，計(jì)算曲率衰減因子W₁：

式中x＝(k₀R_e/2)^1/3φ，k₀為波數(shù)，表示為φ代表閃電源與接收端間的弧度，R_e表示地球半徑，閃電源與接收端的球面距離可以表示為ρ＝R_e×φr。

t_s可表示為：

t_s＝e^-iπ/3τ_s＝e^-iπ/3(3πv_s/2)^2/3

s的取值為0,1,2,3...其中s＝0時(shí)誤差為0，s＝1時(shí)誤差為0.0024，s＝2時(shí)誤差小于0.0012，s取值越大，誤差越小，在大于2時(shí)誤差均小于0.0005。把方程(4)中求出的t_s代入方程(3)，求出的W即為地球曲率衰減因子W₁。

(3)根據(jù)Wait算法，計(jì)算電導(dǎo)率衰減因子W₂：

式中“erfc”表示誤差方程，d是閃電通道與觀測(cè)點(diǎn)之間距離，ω表示角頻率，c表示光速，Δ為土壤等效表面阻抗，表達(dá)式為

ω是信號(hào)角頻率，ε_r為地面相對(duì)介電常數(shù)，ε₀為真空介電常數(shù)，μ₀為真空磁導(dǎo)率，σ為地表電導(dǎo)率。

(4)將計(jì)算得到的衰減因子W₁、W₂與理想頻域場(chǎng)E_W和B_W進(jìn)行卷積得到衰減頻域場(chǎng)E_G和B_G，卷積過(guò)程如下：

E_G＝E_w×W₁×W₂，B_G＝B_w×W₁ (8)

(5)對(duì)頻域衰減場(chǎng)E_G和B_G進(jìn)行傅立葉逆變換，即可得到時(shí)域遠(yuǎn)距離衰減場(chǎng)E_t與B_t。傅立葉逆變換過(guò)程如下：

E_t＝ifft(E_G)，B_t＝ifft (9)

具體實(shí)例

見圖1 300km處理想電場(chǎng)與磁場(chǎng)波形圖：

兩種方法得到的理想場(chǎng)峰值有差距，可用FDTD結(jié)果對(duì)本發(fā)明的算法進(jìn)行修正，使兩種方法計(jì)算出的理想電場(chǎng)相近，以便后面衰減結(jié)果分析。

在采用相同理想場(chǎng)后，接下來(lái)我們要對(duì)兩種算法計(jì)算得到的衰減電場(chǎng)和衰減磁場(chǎng)進(jìn)行對(duì)比分析。圖2是300km處衰減電場(chǎng)與磁場(chǎng)波形圖；

對(duì)比可以看出，兩種方法得到的結(jié)果非常相近，說(shuō)明采用本專利提出的新算法可以準(zhǔn)確計(jì)算遠(yuǎn)距離電磁場(chǎng)的衰減。