戴炳哲, 張其林, 侯文豪, 劉曉東, 宋昊澤, 王曼霏, 姚年鵬

(1.南京信息工程大學氣象災害教育部重點實驗室/氣候與環(huán)境變化國際合作聯(lián)合實驗室/氣象災害預報預警與評估協(xié)同創(chuàng)新中心/中國氣象局氣溶膠與云降水重點開放實驗室，南京 210044;2.內(nèi)蒙古自治區(qū)雷電預警防護中心，呼和浩特 010051)

0 引言

雷電作為常見的瞬時高電壓、強電流自然現(xiàn)象，直接擊中物體易導致火災等災害，產(chǎn)生的電磁場會對架空高壓線、通訊設備、電器電子設備等造成嚴重影響。因而準確定位雷擊點，并明確該雷擊的相關參數(shù)對防護和災害的后續(xù)處理尤為重要。

在對雷電電磁場的初期研究中，國內(nèi)外學者利用簡單近似算法研究其沿地表的傳播情況。Wait 等人[1-3]提出了Wait算法，利用簡單積分公式計算遠距離地閃回擊電磁場和陸-?；旌下窂较碌睦纂婋姶艌?。Shoory和Cooray等人[4]向Wait算法中加入分層電導率，實現(xiàn)了在地表水平分層和垂直分層情況下計算雷電電磁場傳播特性。Cooray和Rubinstein 提出了計算近距離，平坦有限電導地面下地閃回擊電磁場精確解的C-R 算法[5-7]。袁飛等[8]，張其林等[9]分別將一維、二維分型方法和Barrick表面阻抗理論加入Wait算法研究了粗糙地表對地閃回擊電磁場傳播的影響。

隨著計算機設備和技術(shù)的快速發(fā)展，時域有限差分FDTD(Finite Difference Time Domain)方法成了廣泛應用于計算地閃回擊電磁場沿地表傳播情況的方法。楊春山和周璧華[10]驗證了在近距離地閃回擊通道情況下二維FDTD 算法的準確性，并考慮了地面電導率對地閃回擊電磁場的影響。湯宵等[11]利用三維FDTD方法研究了土壤不同電導率分層對雷電感應過電壓的影響。李東帥等[12]采用三維FDTD算法模擬了雷電電磁場沿山體傳播情況，并在后續(xù)工作中[13]采用二維球坐標 FDTD算法研究了雷電電磁場在地-電離層波導中在山區(qū)的傳播效應。黃凱麗等[14]利用二維柱坐標FDTD算法分析了昆明地區(qū)真實地形對閃電定位的影響。上述這些研究揭示了地閃回擊電磁場在傳播過程中會受多種因素，如土壤電導率，真實地形等影響，導致其峰值和峰值時間發(fā)生變化,從而引起定位和雷電參數(shù)反演的誤差。量化這些因素的影響，找出修正誤差的手段就成了值得研究的問題。

由于簡單近似算法計算精確度不高，適用范圍存在局限，同時，有限電導率帶來的影響已在大多數(shù)研究中被探討，而用FDTD對真實地形和地球曲率的影響同時進行分析的很少。因此，筆者采用可以同時將真實地形和地球曲率加入計算的二維球坐標FDTD算法，對比算出波形結(jié)果來體現(xiàn)真實地形與地球曲率的影響，并基于該結(jié)果初步分析了測站數(shù)量和布局以及地閃回擊電磁場脈沖到達時間的不同定義方式對定位精度的影響。

研究結(jié)果對提高閃電定位系統(tǒng)的定位精度和探測效率具有科學意義和參考價值，同時也在使雷擊災害得到及時處理方面有著應用價值。

1 數(shù)據(jù)與方法

1.1 內(nèi)蒙古自治區(qū)地形地貌數(shù)據(jù)

本研究區(qū)域位于內(nèi)蒙古高原，平均海拔1 000 m～1 200 m。采用的地形數(shù)據(jù)來源于SRTM3(Shuttle Radar Topography Mission 3)，由美國航空航天局和國防部國家測繪局聯(lián)合測量得到，相鄰兩采樣點間隔90 m。測站位置為2019年在內(nèi)蒙古自治區(qū)中部區(qū)域建成的測站實際位置，所在區(qū)域真實地形圖見圖1，圖2是雷擊點到各站的地形剖面圖(a-k分別表示東勝，臨河，化德、土左旗、清水河、豐鎮(zhèn)、察右中旗、滿都拉、達茂、大佘太、四子王旗)，可以看出地閃回擊電磁場到各個不同測站除距離不盡相同，傳播路徑也很復雜。

圖2 模擬雷擊點到各觀測站的地形剖面圖Fig.2 Terrain profiles from the simulated lightning strike point to each observation station

1.2 計算方法

1.2.1 蒙特卡洛法

蒙特卡洛法[15,25]通過形成大量隨機樣本進行模擬來了解一個系統(tǒng)的特性。本研究選取38°N-43.5°N，107°E-114.5°E區(qū)域，間隔為0.01°的所有格點。每個格點在其以光速到各站的基礎時間上添加均值為0，方差為200 ns2的正態(tài)分布誤差后利用TOA(Time of Arrival)技術(shù)[16-17]定位，重復進行100次,完成后，算出定位平均誤差。全部完成后，即可了解研究區(qū)域內(nèi)11個測站的定位誤差特征。

1.2.2 二維球坐標時域有限差分法

采用二維球坐標系時域有限差分算法[18-20]，模擬地閃回擊電磁場在地-電離層波導中的傳播。圖3為該方法的整體構(gòu)型，以地球中心為原點建立球坐標，閃電通道位于仿真域左邊緣，上下邊界由卷積完全匹配層CPML (Convolutional Perfect Match Layer)吸收邊界[21]包圍。假定土壤均勻，電導率為σr，相對介電常數(shù)為εr。地表以上的區(qū)域是空氣和電離層。

圖3 二維球坐標時域有限差分法整體構(gòu)型Fig.3 The framework of two-dimensional spherical coordinates FDTD method

計算中觀測點選在各測站實際位置，雷擊點到各測站距離分別為162.7 km,310.6 km,265.1 km,38.6 km,136.9 km,191.4 km,136 km,180.3 km,87.3 km,159.4 km和75.4 km(列舉順序同圖2中a-k順序)。計算網(wǎng)格大小為50 m，考慮到不發(fā)生色散的條件是最高頻率電磁波對應的波長大于10倍的網(wǎng)格長度，故計算中可以模擬的電磁波最高頻率為600 kHz，高于實際探測儀器使用的最高頻段。地閃回擊采用MTLE (Modified Transmission Line Model with Exponential Current Decay with Height)模型[22]，此模型假定回擊電流幅值隨通道高度以指數(shù)形式衰減，在t時刻通道高度z′處的雷電流可表示為

i(z′,t)=e-z′/λi(0,t-z′/v)

(1)

式(1)中：e-z′/λ為指數(shù)衰減系數(shù)，λ為衰減因子，取2 000 m；v為回擊速度，取1.5×108m/s。通道底部的基電流波形采用雙Heidler函數(shù)模型：

(2)

為驗證本研究采用的計算方法，將我們的計算結(jié)果與Tran等人計算的結(jié)果[23]進行了比較。由于Tran等沒有考慮地球的曲率，因而在計算中設置了一個非常大的地球半徑值來表示平坦地面。圖4是在夜間，平坦地面，完全電導情況下距離200 km地面處計算出的波形的對比。可以看出，兩種方法計算的波形基本一致，本研究使用的方法是有效的。

圖4 本文計算結(jié)果與Tran等人計算結(jié)果比較Fig.4 Calculation result in this paper compared with Tran’s

2 結(jié)果分析

2.1 二維球坐標時域有限差分法結(jié)果

圖5是各站處考慮3種不同情況時計算的磁場時域波形。對考慮曲率的真實地形情況的波形進行頻域分析，結(jié)果展現(xiàn)在圖6中。因為頻率的變化隨傳播距離的增加基本一致，所以圖6中僅顯示離雷擊點最遠的臨河站和最近的土左旗站的結(jié)果(其時域波形分別對應圖5(b)和圖5(d))。發(fā)現(xiàn)在真實地形和曲率作用下地閃回擊電磁場整體波形中較高頻部分衰減程度大于更低頻率的部分，隨著距離的增加波形幅值降低，低頻部分逐漸占據(jù)更大比例。進而在圖5(b)和圖5(d)顯示的時域波形里表現(xiàn)為波形整體變寬，波形上升沿時間相較于光滑地表呈現(xiàn)明顯滯后且上升緩慢的趨勢。這一現(xiàn)象在各站間光滑地面情況下的波形對比中并不顯著。

圖5 二維球坐標FDTD方法計算出的各站3種情況下的磁場波形Fig.5 Lightning magnetic field waveform of each station calculated by two-dimensional spherical coordinates FDTD method

圖6 臨河、土左旗站真實地形有曲率情況下波形頻域分析Fig.6 Frequency domain analysis of waves in Linhe station and Tuzuoqi station considering real terrain and earth curvature

將各站本身不同情況下的3種波形峰值時間相互對比，真實地形和地球曲率均造成了波形上升沿時間的額外增加。見圖7，真實地形引起的時間延遲變化范圍從2.91 μs至4.50 μs；對比各站本身考慮曲率時光滑地面與不考慮曲率時光滑地面結(jié)果，當閃電與測站間距小于30 km時曲率造成的時間延遲接近為零，可以忽略，但當間距超過300 km時，曲率帶來的時延就已達到將近1 μs，此時曲率的影響不應忽視。

圖7 地球曲率、真實地形在各測站帶來的時間延遲Fig.7 Time delay caused by the curvature of the earth, the real terrain at each measuring station

波形峰值方面，地球曲率造成了峰值的衰減，且衰減隨著距離增加而增大。各測站本身考慮曲率和不考慮曲率的光滑地面的兩種情況的波形結(jié)果對比顯示，在距離最遠的臨河站處，衰減可達13.78%，顯然此間距下曲率對波形峰值的影響不應忽視。值得注意的是，在各站的結(jié)果中，與其本身光滑地面情況下的結(jié)果相比，考慮真實地形情況時波形的峰值均表現(xiàn)為增大，這可能是地形對電磁場反射作用所致[24]。四子王旗站處增加百分比最大，達到了42.98%，東勝站處最小為24.95%，平均增大約33.08%。很顯然為了準確反演相關參數(shù)，地球曲率和真實地形對峰值的影響也需要在反演中加以訂正。

從上述結(jié)果中可以看出，在間距比較大的測站組網(wǎng)觀測時，除了考慮地面電導率等帶來的影響，真實地形與地球曲率的影響也不容忽視。

2.2 蒙特卡洛結(jié)果

為驗證定位程序的有效性，蒙特卡洛分加誤差和不加誤差兩次進行，結(jié)果分別展現(xiàn)在圖8(a)、圖8 (b)中。圖8顯示在理想到達時間的情況下定位極為準確，故定位程序有效。從圖8中可以看出11站定位結(jié)果在研究區(qū)域內(nèi)左下和右下部分會出現(xiàn)比較大的偏差。同時，對比圖8(a)和圖8 (b)，可以看到較小時間誤差變化就會給給定位帶來巨大影響，故考慮各種因素帶來的影響，采用合理方法減小因此產(chǎn)生的定位誤差十分必要。

圖8 蒙特卡洛結(jié)果Fig.8 Result of Monte Carlo method

2.3 定位誤差分析

選擇合理的地閃回擊電磁場脈沖到達時間定義方法所獲得的到達時間可以減小定位上的誤差，因此如何選擇是一個值得探討的問題。根據(jù)相關研究[12]，最常用的計算脈沖到達時間的公式為

(3)

式中：tp對應測得磁場脈沖峰值Hp所在時刻；tT對應觸發(fā)閾值Hth所處時刻。

從公式(3)可以發(fā)現(xiàn)，由于到達單一測站的脈沖波形是確定的，那么影響計算脈沖到達時間的最主要就是觸發(fā)閾值Hth的選取，選取5種常見的觸發(fā)閾值來計算到達時間進行定位分析：1) 磁場峰值的10%；2) 磁場峰值的20%；3) 磁場峰值的50%； 4) 磁場峰值；5) 磁場一階導數(shù)的峰值。為說明方便將這幾種選取方法簡稱為方法1-5。用11站真實地形情況下的波形進行定位時，方法1誤差最小為42.1 m,方法5次之，為42.6 m。方法1和方法5的誤差較小，可能是由于這兩種定義方法受波形上升沿時間額外增加的影響相對較小。需要注意的是文中的誤差均是指定位點與模擬雷擊點之間的水平誤差，這是因為由于測站基本是處在同一平面的，在低空的垂直誤差本身就比較大[25]。

定位誤差除與地閃回擊電磁場脈沖到達時間的定義方法有關外，還與參與定位測站的數(shù)量和布局有關。用方法1，采用特定的10站進行定位時，水平誤差只有1.41 m，采用3站時也能達到24.2 m，均優(yōu)于直接使用所有11站波形時的結(jié)果。圖9是利用不同數(shù)量測站進行定位后的定位誤差棒圖，圖9(a)是標準差結(jié)果，圖9(b)展現(xiàn)了定位時最大與最小誤差結(jié)果。由于使用4站和3站定位誤差結(jié)果的標準差與最大誤差相比其他過大，圖中沒有將它們的結(jié)果畫出?？梢钥闯?，定位誤差的平均值和標準差隨著測站數(shù)量的減少而增大，誤差的變化范圍也越來越大，但是在采用特定布局測站的情況下，定位最小誤差仍然能夠保持在較低水平。進而可以看出測站數(shù)量越少，測站布局對于定位精度的影響越大，反之則越小，因而在使用較少的測站定位時測站布局的影響更不可忽略。

圖9 不同數(shù)量測站定位誤差棒圖Fig.9 Location error bar graph of different numbers of stations

各測站實際運行中，總存在只有部分測站能接收到地閃回擊電磁脈沖同步波形的情況，這時考慮測站數(shù)量和布局的影響就對精確定位顯得尤為重要。在利用所有接收到同步波形的測站得到雷擊點大致位置后，選擇數(shù)量更少但布局更合理的部分測站可以進一步提升定位的精確度。

3 結(jié)論

采用二維球坐標時域有限差分方法模擬了真實地形和地球曲率對地閃回擊電磁場的影響，并將這些影響和測站數(shù)量與布局以及地閃回擊電磁場脈沖到達時間的定義方法相結(jié)合分析了對閃電定位的影響。研究中發(fā)現(xiàn)：

1) 真實地形和地球曲率都會引起額外的地閃回擊磁場脈沖峰值到達時間延遲，各站真實地形造成的延遲時間平均為3.76 μs，變化范圍從2.91 μs至4.50 μs，地球曲率導致時延平均0.37 μs，變化范圍從0.08 μs至0.92 μs，且隨電磁波傳播距離的增大而增大。

2) 研究中，各站處地球曲率均帶來了地閃回擊電磁場脈沖強度的衰減，相比該站本身考慮曲率的光滑地面情況，真實地形增大了地閃回擊電磁場脈沖強度，其影響大于地球曲率的影響。

3)在考慮真實地形和地球曲率時，測站數(shù)量和布局以及地閃回擊磁場脈沖到達時間的不同定義方式，都對地閃回擊定位精度存在影響，隨著測站數(shù)量減少，測站布局對于定位精度的影響增大，這時挑選形成更合理布局的測站就對定位精度的保證十分重要。

為提高觀測區(qū)域內(nèi)地閃定位和參數(shù)反演的準確性，需要將觀測區(qū)域細化成多個小區(qū)域并選取其中有代表性的點進行模擬，明確各小區(qū)域內(nèi)這些因素的影響，進而找出定位和反演的最優(yōu)方案。另外，在實際定位中選擇數(shù)量更少但布局更合理測站或許可以進一步提升定位的精確度。

不過本研究的時域有限差分計算中尚未考慮傳播路徑中不同土壤類型帶來的影響，也暫時沒有與測站實際觀測獲取的數(shù)據(jù)進行對比。另外，如何利用并行計算方法大幅提升模擬計算的速度，也是目前正在解決的問題，這些不足將在后續(xù)的工作中進一步完善。