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

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

0 引言

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

在對(duì)雷電電磁場(chǎng)的初期研究中，國(guó)內(nèi)外學(xué)者利用簡(jiǎn)單近似算法研究其沿地表的傳播情況。Wait 等人[1-3]提出了Wait算法，利用簡(jiǎn)單積分公式計(jì)算遠(yuǎn)距離地閃回?fù)綦姶艌?chǎng)和陸-海混合路徑下的雷電電磁場(chǎng)。Shoory和Cooray等人[4]向Wait算法中加入分層電導(dǎo)率，實(shí)現(xiàn)了在地表水平分層和垂直分層情況下計(jì)算雷電電磁場(chǎng)傳播特性。Cooray和Rubinstein 提出了計(jì)算近距離，平坦有限電導(dǎo)地面下地閃回?fù)綦姶艌?chǎng)精確解的C-R 算法[5-7]。袁飛等[8]，張其林等[9]分別將一維、二維分型方法和Barrick表面阻抗理論加入Wait算法研究了粗糙地表對(duì)地閃回?fù)綦姶艌?chǎng)傳播的影響。

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

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

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

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

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

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

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

1.2 計(jì)算方法

1.2.1 蒙特卡洛法

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

1.2.2 二維球坐標(biāo)時(shí)域有限差分法

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

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

計(jì)算中觀測(cè)點(diǎn)選在各測(cè)站實(shí)際位置，雷擊點(diǎn)到各測(cè)站距離分別為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順序)。計(jì)算網(wǎng)格大小為50 m，考慮到不發(fā)生色散的條件是最高頻率電磁波對(duì)應(yīng)的波長(zhǎng)大于10倍的網(wǎng)格長(zhǎng)度，故計(jì)算中可以模擬的電磁波最高頻率為600 kHz，高于實(shí)際探測(cè)儀器使用的最高頻段。地閃回?fù)舨捎肕TLE (Modified Transmission Line Model with Exponential Current Decay with Height)模型[22]，此模型假定回?fù)綦娏鞣惦S通道高度以指數(shù)形式衰減，在t時(shí)刻通道高度z′處的雷電流可表示為

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

(1)

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

(2)

為驗(yàn)證本研究采用的計(jì)算方法，將我們的計(jì)算結(jié)果與Tran等人計(jì)算的結(jié)果[23]進(jìn)行了比較。由于Tran等沒(méi)有考慮地球的曲率，因而在計(jì)算中設(shè)置了一個(gè)非常大的地球半徑值來(lái)表示平坦地面。圖4是在夜間，平坦地面，完全電導(dǎo)情況下距離200 km地面處計(jì)算出的波形的對(duì)比?？梢钥闯觯瑑煞N方法計(jì)算的波形基本一致，本研究使用的方法是有效的。

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

2 結(jié)果分析

2.1 二維球坐標(biāo)時(shí)域有限差分法結(jié)果

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

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

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

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

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

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

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

2.2 蒙特卡洛結(jié)果

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

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

2.3 定位誤差分析

選擇合理的地閃回?fù)綦姶艌?chǎng)脈沖到達(dá)時(shí)間定義方法所獲得的到達(dá)時(shí)間可以減小定位上的誤差，因此如何選擇是一個(gè)值得探討的問(wèn)題。根據(jù)相關(guān)研究[12]，最常用的計(jì)算脈沖到達(dá)時(shí)間的公式為

(3)

式中：tp對(duì)應(yīng)測(cè)得磁場(chǎng)脈沖峰值Hp所在時(shí)刻；tT對(duì)應(yīng)觸發(fā)閾值Hth所處時(shí)刻。

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

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

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

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

3 結(jié)論

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

1) 真實(shí)地形和地球曲率都會(huì)引起額外的地閃回?fù)舸艌?chǎng)脈沖峰值到達(dá)時(shí)間延遲，各站真實(shí)地形造成的延遲時(shí)間平均為3.76 μs，變化范圍從2.91 μs至4.50 μs，地球曲率導(dǎo)致時(shí)延平均0.37 μs，變化范圍從0.08 μs至0.92 μs，且隨電磁波傳播距離的增大而增大。

2) 研究中，各站處地球曲率均帶來(lái)了地閃回?fù)綦姶艌?chǎng)脈沖強(qiáng)度的衰減，相比該站本身考慮曲率的光滑地面情況，真實(shí)地形增大了地閃回?fù)綦姶艌?chǎng)脈沖強(qiáng)度，其影響大于地球曲率的影響。

3)在考慮真實(shí)地形和地球曲率時(shí)，測(cè)站數(shù)量和布局以及地閃回?fù)舸艌?chǎng)脈沖到達(dá)時(shí)間的不同定義方式，都對(duì)地閃回?fù)舳ㄎ痪却嬖谟绊?，隨著測(cè)站數(shù)量減少，測(cè)站布局對(duì)于定位精度的影響增大，這時(shí)挑選形成更合理布局的測(cè)站就對(duì)定位精度的保證十分重要。

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

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