鄒 斌, 崔海華, 張其林*, 李 靜, 戴炳哲,姚年鵬, 蔣 堯, 楊 敢, 羅 華
(1.南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室, 南京 210044; 2.南京信息工程大學(xué)氣候與環(huán)境變化國際合作 聯(lián)合實(shí)驗(yàn)室, 南京 210044; 3.南京信息工程大學(xué)氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心, 南京 210044; 4.南京信息工程大學(xué)中國氣象局氣溶膠與云降水重點(diǎn)開放實(shí)驗(yàn)室, 南京 210044; 5.河北省氣象行政技術(shù)服務(wù)中心, 石家莊 050000)
雷電是一種常見的災(zāi)害性天氣,每年造成大量的損失[1]。因此,對(duì)雷電進(jìn)行深入的研究,提高閃電探測能力和定位精度,對(duì)降低雷電災(zāi)害損失,保障人類生產(chǎn)活動(dòng)的有序進(jìn)行和人員的生命財(cái)產(chǎn)安全具有積極的意義[2]。然而現(xiàn)實(shí)中除了引雷實(shí)驗(yàn)和極少數(shù)地方有引雷塔,可以對(duì)現(xiàn)實(shí)中雷電定位算法的定位誤差進(jìn)行評(píng)估,對(duì)于大部分的雷電很難得知其定位誤差,而隨著計(jì)算機(jī)技術(shù)的發(fā)展,模擬算法為評(píng)估雷電定位誤差提供了有效途徑。雷電電磁波在近地面?zhèn)鞑ミ^程中受土壤色散、土壤電導(dǎo)率以及高低起伏的不規(guī)則地形等因素的影響較大,會(huì)導(dǎo)致上升沿時(shí)間、電磁波波形峰值等發(fā)生改變,從而降低了雷擊點(diǎn)的定位精度[3]。因此,利用二維時(shí)域有限差分法(finite difference time domain,F(xiàn)DTD)研究雷電電磁波傳播特性并進(jìn)行定位誤差分析有重要意義。
在雷電回?fù)綦姶艌龅难芯恐?,主要包括精確解析計(jì)算、近似計(jì)算和數(shù)值計(jì)算3種求解方法。Wait[4]利用等效表面阻抗定義的衰減函數(shù),給出了垂直電場的近似計(jì)算公式。Cooray[5]和Rubinstein[6]給出了水平電場的近似計(jì)算公式。此后,經(jīng)過學(xué)者們不斷地改進(jìn),近似算法被不斷推廣,可以應(yīng)用到分層土壤、粗糙地表和海陸交接等多種地形環(huán)境中[7-10]。然而,解析法在山體等非規(guī)則地形中的電磁場計(jì)算仍受到較多限制。
中外學(xué)者已利用FDTD算法,圍繞土壤電導(dǎo)率[11]、土壤各向異性[12-13]、土壤色散效應(yīng)[14]和地形[15-16]等多個(gè)因素進(jìn)行了大量的分析討論。張明霞等[17]利用矩量法,討論了矩形山體模型對(duì)電磁輻射磁場的影響。Li等[18]在研究地形對(duì)雷電定位的影響中發(fā)現(xiàn)當(dāng)傳播路徑中有錐形山體存在時(shí),發(fā)現(xiàn)電磁場波形和峰值均會(huì)受到嚴(yán)重影響,并且信號(hào)的時(shí)延隨山體高度的增加而增大,隨地面電導(dǎo)率的減小而減小。Li 等[19]通過結(jié)合雷擊山頂鐵塔的實(shí)測電場波形和FDTD數(shù)值模擬,討論了真實(shí)不平坦地形的影響,證實(shí)了山體地形對(duì)電場(15 km觀測點(diǎn))具有增強(qiáng)作用,結(jié)果表明忽略地形起伏的平地假設(shè)導(dǎo)致對(duì)反演估算的雷電流幅值達(dá)到真實(shí)值的1.8倍。陶玉郎等[20]利用3D數(shù)值模擬研究發(fā)現(xiàn)土壤的特征對(duì)雷電流沖擊有很大影響。Hou等[21]建立2D FDTD錐形山體模型比較了0.3~50 km處電磁場的變化規(guī)律,發(fā)現(xiàn)對(duì)于垂直電場,近距離受山體的屏蔽作用影響被削弱,而遠(yuǎn)距離處的幅值則由于山腳的反射而增強(qiáng)。張金波等[22]利用2D FDTD算法對(duì)土壤電導(dǎo)率進(jìn)行研究,結(jié)果表明山體地形坡度對(duì)雷電過電壓有顯著影響,尤其在土壤電導(dǎo)率為有限值的情況下,過電壓幅值隨地形坡度的升高而明顯增大。
總的來說,以上研究大部分都局限在孤立山體,涉及實(shí)際山體地形的雷電電磁波傳播特性較少。雷電電磁場沿復(fù)雜起伏地形傳播時(shí),會(huì)以繞射和反射的方式經(jīng)過山體,會(huì)引起電磁場到達(dá)測站時(shí)間的增加,另外電磁場也會(huì)出現(xiàn)強(qiáng)度的變化,從而使根據(jù)雷電電磁場波形進(jìn)行的閃電定位產(chǎn)生誤差。因此如何將地形地貌帶來的誤差進(jìn)行修正,提高閃電的探測效率和定位精度,是個(gè)亟需解決的問題。
為此,在二維柱坐標(biāo)系下FDTD模型加入了真實(shí)地形,研究了真實(shí)地形對(duì)切向磁場波形以及閃電定位的影響。針對(duì)現(xiàn)有的TOA(time of arrival)時(shí)差定位算法的缺陷,即忽略了地形對(duì)電磁波的到達(dá)時(shí)間延遲作用,提出了不同步長的地形包絡(luò)法,對(duì)定位精度進(jìn)行優(yōu)化。
海南地區(qū)位于中國南端,主要以山地為主,地勢起伏比較嚴(yán)重,山地面積較大,緯度低四面環(huán)海,雷暴活動(dòng)相對(duì)比較頻繁,是中國雷電災(zāi)害最為嚴(yán)重的省份之一。以實(shí)際安裝在海南省的7個(gè)測站位置為基礎(chǔ),將它們應(yīng)用于FDTD模式中,7個(gè)測站分別為瓊中站、樂東站、三亞站、昌江站、萬寧站、文昌站、澄邁站,具體測站位置如表1所示。在閃電探測網(wǎng)區(qū)域選取了3個(gè)雷擊點(diǎn)L1、L2、L3,具體位置如表2所示。
表1 測站坐標(biāo)Table 1 Station coordinates
表2 雷擊點(diǎn)坐標(biāo)Table 2 Lightning coordinates
圖1為海南地區(qū)測站和雷擊點(diǎn)位置,在閃電探測網(wǎng)區(qū)域選取了3個(gè)雷擊點(diǎn)L1、L2、L3,其中雷擊點(diǎn)L1在中間山地區(qū)域,雷擊點(diǎn)L2在東北平原地區(qū),雷擊點(diǎn)L3在西南山地附近的盆地區(qū)域,選取這3個(gè)雷擊點(diǎn)具有代表性。
FDTD是由麥克斯韋旋度方程的微分方程出發(fā)的,通過二階中心差分方法,把原來的偏微分方程運(yùn)算形式轉(zhuǎn)換成差分方程,從而在一定的時(shí)間上和一定體積內(nèi),對(duì)連續(xù)的電磁場數(shù)據(jù)可以進(jìn)行抽樣壓縮。
Maxwell方程組的旋度方程組形式為
(1)
式(1)中:E為電場強(qiáng)度;H為磁場強(qiáng)度;D為電位移矢量;ρ為電荷密度;B為磁感應(yīng)強(qiáng)度;ε為電介質(zhì)常數(shù);σ為電導(dǎo)率;μ為磁導(dǎo)系數(shù);t為時(shí)間。直角坐標(biāo)系中,式(1)可表示為
(2)
式(2)中:x、y、z分別為直角坐標(biāo)系的3個(gè)方向。
采用磁場波形進(jìn)行閃電定位,則切向磁場的表達(dá)式為[23]
(3)
式(3)中:Δta為時(shí)間步長;Δr為水平方向網(wǎng)格長度;Δz為豎直方向格點(diǎn)長度;i為水平方向格點(diǎn)數(shù);j為豎直方向格點(diǎn)數(shù);Hφ為切向磁場;Ez為垂直電場;Er為水平電場;n為第n個(gè)計(jì)算單元;μ0為真空中的磁導(dǎo)率(即真空中的磁導(dǎo)系數(shù))。
圖2為FDTD模式示意圖,模式中網(wǎng)格寬度Δr×Δz=10 m×10 m,滿足數(shù)值穩(wěn)定性條件,同時(shí)也與插值后的地形網(wǎng)格數(shù)據(jù)相匹配。時(shí)間步長Δt=16.7 ns,符合Courant穩(wěn)定性條件。土壤電導(dǎo)率σg=0.001 S/m,相對(duì)介電常數(shù)εr=10??諝怆妼?dǎo)率σair=0 S/m,相對(duì)介電常數(shù)εr=1。邊界采用卷積完全匹配層(convolution perfectly matched layer,CPML)吸收邊界條件。閃電回?fù)敉ǖ婪胖糜谟?jì)算區(qū)域的最左側(cè)山頂?shù)纳厦?,通道高度H=7 500 m。
z為垂直方向;φ為切向;r為水平方向圖2 模型示意圖Fig.2 Model diagram
采用MTLE(modified transmission line)源回?fù)裟P停椿負(fù)敉ǖ乐械幕負(fù)綦娏饕灾笖?shù)的形式沿著通道高度在衰減,回?fù)敉ǖ乐衪時(shí)刻在高度z′處的電流分布可表示為
i(z′,t)=e-z′/λi(0,t-z′/ν)
(4)
式(4)中:e-z′/λ為數(shù)值衰減系數(shù);i為電流強(qiáng)度;z′為高度;λ為衰減因子,取2 000 m;ν為回?fù)羲俣?,?.5×108m/s。
使用Heidler雙指數(shù)函數(shù)[24]來模擬繼后回?fù)粝吕纂娏魍ǖ赖撞康幕娏?,如?5)所示。典型的繼后回?fù)舻碾娏鲄?shù)如表3所示。
表3 繼后回?fù)綦娏鲄?shù)典型值Table 3 Typical value of subsequent return stroke current parameters
(5)
(6)
(7)
式中:i01和i02分別為擊穿電流和電暈電流峰值;η1、η2分別為擊穿電流和電暈電流峰值的修正因子;τ11和τ21分別為擊穿電流和電暈電流的波形上升時(shí)間;τ12和τ22分別為擊穿電流和電暈電流的波形下降時(shí)間。
時(shí)差法是利用電磁信號(hào)到達(dá)不同測站的時(shí)間差進(jìn)行定位,結(jié)合高精度的全球定位系統(tǒng)(global positioning system,GPS)同步時(shí)間測量,能夠?qū)崿F(xiàn)對(duì)輻射源點(diǎn)的準(zhǔn)確定位。在利用TOA 算法計(jì)算輻射原點(diǎn)的過程中,首先對(duì)信號(hào)進(jìn)行波形識(shí)別篩選出不同類型的閃電事件,同時(shí)對(duì)不同的同步脈沖進(jìn)行分組,再計(jì)算出同步信號(hào)到達(dá)不同測站的時(shí)間差,然后假定信號(hào)沿直線距離以光速向測站傳播,求解非線性方程組獲得輻射源點(diǎn)位置的初值,并借助非線性最小二乘擬合算法對(duì)初解進(jìn)行優(yōu)化,最終得到輻射源點(diǎn)的位置。圖3為到達(dá)時(shí)間示意圖。
Upeak為閃電峰值電壓;tt為閃電上升沿開始時(shí)間; tpeak為波形峰值時(shí)間,即逐峰法定義的時(shí)間圖3 到達(dá)時(shí)間示意圖Fig.3 Diagram of arrival time
為了求解出地閃回?fù)酎c(diǎn)的空間坐標(biāo),需要先得到信號(hào)到達(dá)不同測站的時(shí)間差。在時(shí)間差的計(jì)算中,可以通過互相關(guān)方法直接計(jì)算出兩兩測站之間的時(shí)間差,也可以以電磁波峰值的到達(dá)時(shí)間為準(zhǔn),先計(jì)算出單站信號(hào)的到達(dá)時(shí)間,再與參考站點(diǎn)的到達(dá)時(shí)間相減獲得兩個(gè)測站之間的時(shí)間差。不同方法計(jì)算得到時(shí)間差可能存在一定的差異,為了對(duì)比不同方法的定位效果,選取互相關(guān),逐峰法這2種方法的定位結(jié)果進(jìn)行比較。
雖然對(duì)閃電定位普遍采用的是TOA時(shí)差法,但顯然這種算法并沒有考慮實(shí)際地形對(duì)閃電電磁波產(chǎn)生的路徑延長問題,即信號(hào)到達(dá)每個(gè)測站的時(shí)間都會(huì)比理想情況下略長,因此采用不同步長的地形包絡(luò)法來近似估計(jì)地形引起的延長時(shí)間。
如圖4所示,地形包絡(luò)即用連續(xù)線段將地形包圍起來??紤]到所有測站和雷擊點(diǎn)的距離大小,每隔5、10、20、40 km步長距離取每段地形的最高點(diǎn)并連線形成包絡(luò)線。將這些地形包絡(luò)線長度作為電磁波在傳播中的近似傳播距離,反映出信號(hào)傳播路徑隨山體高度的變化。然后用包絡(luò)線長度減去水平距離后除以電磁波的傳播速度作為電磁波到達(dá)時(shí)間的誤差訂正值。假設(shè)真實(shí)路徑下用包絡(luò)法得到測站至雷擊點(diǎn)路徑長度為l1,水平路徑為l2,可以得到誤差訂正時(shí)間Δtb的計(jì)算公式為
(8)
式(8)中:c為光速。
將原各個(gè)測站的信號(hào)到達(dá)時(shí)間分別減去對(duì)應(yīng)的誤差訂正時(shí)間,重新進(jìn)行雷擊點(diǎn)定位。
忽略測站GPS 計(jì)時(shí)精度的誤差,信號(hào)從雷擊點(diǎn)向測站傳播的時(shí)間t,會(huì)因?yàn)榈匦巍㈦妼?dǎo)率等的影響而變化,與同等水平傳播距離的理想平坦地表時(shí)的傳播時(shí)間t存在一定的偏差,即
t=tflat+Δtc
(9)
式(9)中:tflat為信號(hào)以光速傳播過平坦地表的時(shí)間;Δtc為真實(shí)崎嶇地形情況與理想平坦地表情況的傳播時(shí)間偏差。
進(jìn)行時(shí)間補(bǔ)償時(shí),將信號(hào)以光速沿紅色折線和黑色直線傳播的時(shí)間差作為對(duì)式(8)中Δt的估算,減少傳播時(shí)間t與理想情況下的tflat的偏差。有必要時(shí),考慮進(jìn)行多次補(bǔ)償,提高最終定位結(jié)果的準(zhǔn)確性。
圖5~圖7為L1、L2、L3雷擊點(diǎn)到測站的真實(shí)地形高程圖,可以看出,雷擊點(diǎn)L1雷電電磁波傳播路徑的總體趨勢為由高到低,雷擊點(diǎn)在山地上,傳播路徑總體上比較復(fù)雜,地形起伏較大;雷擊點(diǎn)L2在平原地區(qū),傳播的總體趨勢為由低到高,總體上受山地影響較??;雷擊點(diǎn)L3在萬寧測站附近,雷擊點(diǎn)海拔高度低,但其四周地勢都比雷擊點(diǎn)高,在山地的低洼區(qū)域,受山地影響更明顯。
圖5 雷擊點(diǎn)L1到測站地形圖Fig.5 Topographic map of lightning L1 to the station
圖7 雷擊點(diǎn)L3到測站地形圖Fig.7 Topographic map of lightning L3 to the station
圖8~圖10為雷擊點(diǎn)L1、L2、L3到每個(gè)測站平地與真實(shí)地形的波形對(duì)比。藍(lán)色線和紅色線分別為平地和真實(shí)地形的模擬結(jié)果。山地地形對(duì)波形的上升沿時(shí)間和波形峰值影響較大。上升沿時(shí)間增加范圍在1~5 μs,主要是電磁波在傳播中受土壤有限電導(dǎo)率和色散效應(yīng)影響,高頻電磁場存在更大的衰減。同樣峰值時(shí)間也都是滯后于平地的,但峰值的大小變化卻不盡相同。
圖9 雷擊點(diǎn)L2到測站磁場波形Fig.9 Magnetic field waveform from lightning L2 to the station
圖10 雷擊點(diǎn)L3到測站磁場波形Fig.10 Magnetic field waveform from lightning L3 to the station
對(duì)于雷擊點(diǎn)L1,瓊中站真實(shí)地形下峰值明顯比平地情況大很多,峰值時(shí)間較平地情況也提前,而其他6個(gè)測站真實(shí)地形下的峰值對(duì)比與平地都略有減少,真實(shí)地形下峰值時(shí)間明顯滯后于平地傳播的峰值,這主要受山地影響,雷電高頻電磁波衰減較大,而瓊中站距雷擊點(diǎn)距離較近,受山地影響較小,這也說明了地形對(duì)電磁波影響的復(fù)雜性。
對(duì)于雷擊點(diǎn)L2,7個(gè)測站真實(shí)地形下的峰值都略小于平地情況,峰值時(shí)間也都滯后,三亞站和昌江站滯后比較明顯,峰值減小也較明顯,其他5個(gè)測站真實(shí)地形下的電磁波波形及到達(dá)時(shí)間基本上一致,這主要受傳播路徑影響,閃電從雷擊點(diǎn)L2到三亞站和昌江站受山地影響較大,而從雷擊點(diǎn)L2到另外5個(gè)測站地形相對(duì)比較平坦。
對(duì)于雷擊點(diǎn)L3, 6個(gè)測站真實(shí)地形下的峰值都明顯小于平地情況,峰值時(shí)間除樂東站外也都明顯滯后,主要是樂東站從雷擊點(diǎn)L3到測站傳播路徑比較平坦,和雷擊點(diǎn)L2一致。
考慮到雷電探測傳感器探測都具有特定的帶寬,實(shí)際接收的波形已被過濾。在此,使用帶通濾波器,對(duì)2-D FDTD模型的模擬磁場波形進(jìn)行過濾,以分析帶寬對(duì)雷電磁場的影響。圖11~圖13分別為雷擊點(diǎn) L1~L3 在真實(shí)地形下,各測站接收的磁場波形過濾前后對(duì)比,頻率范圍為10~500 kHz。與過濾前相比,磁場峰值減小,波形變陡,上升沿時(shí)間變短,而峰值到達(dá)時(shí)間與過濾前相比幾乎相等,這樣對(duì)于逐峰法定位來說,誤差幾乎相似。
圖11 雷擊點(diǎn)L1真實(shí)地形下到測站濾波前后磁場波形Fig.11 Waveforms of magnetic field before and after filtering at lightning point L1 in real terrain
圖12 雷擊點(diǎn)L2真實(shí)地形下到測站濾波前后磁場波形Fig.12 Waveforms of magnetic field before and after filtering at lightning point L2 in real terrain
圖13 雷擊點(diǎn)L3真實(shí)地形下到測站濾波前后磁場波形Fig.13 Waveforms of magnetic field before and after filtering at lightning point L3 in real terrain
圖14為在逐峰法波形到達(dá)時(shí)間算法下3個(gè)雷擊點(diǎn)在平地、真實(shí)地形、以及不同步長地形包絡(luò)法修正后的定位誤差對(duì)比圖??梢钥闯觯讚酎c(diǎn)L2的誤差較小,主要是雷擊點(diǎn)L2到測站的傳播受地形影響較小,且用包絡(luò)法進(jìn)行修正,隨著包絡(luò)空間步長的增加,定位誤差在減??;對(duì)于雷擊點(diǎn)L1和L3,定位誤差明顯大于雷擊點(diǎn)L2的定位誤差,對(duì)于平坦地表,定位誤差只有十幾米,定位效果較好,而使用包絡(luò)法和時(shí)間補(bǔ)償法進(jìn)行修正,定位誤差明顯減小,且包絡(luò)空間步長越小,定位誤差越小。
圖14 逐峰法定位誤差Fig.14 Location error by peak method
圖15為在互相關(guān)法波形到達(dá)時(shí)間算法下3個(gè)雷擊點(diǎn)在平地、真實(shí)地形、以及不同步長地形包絡(luò)法修正后的定位誤差對(duì)比??梢钥闯?,雷擊點(diǎn)L1的定位誤差相比較其他兩個(gè)雷擊點(diǎn)要小很多,雷擊點(diǎn)L3使用六站進(jìn)行定位相比較雷擊點(diǎn)L1的7站定位,定位誤差明顯大了很多;對(duì)于3個(gè)雷擊點(diǎn),使用包絡(luò)法和時(shí)間補(bǔ)償法進(jìn)行修正,定位誤差明顯減小,且包絡(luò)空間步長越小,定位誤差越小。
圖15 互相關(guān)法定位誤差Fig.15 Location error by cross-correlation method
對(duì)比3個(gè)雷擊點(diǎn)在兩種定位算法下的定位誤差,從圖14、圖15可以看出:對(duì)于平坦地表,雷擊點(diǎn)逐峰法下的定位效果更好,定位誤差都只有幾十米,而對(duì)于雷擊點(diǎn)L1,因?yàn)閭鞑ナ艿匦斡绊戄^大,使用互相關(guān)算法定位誤差較好,而雷擊點(diǎn)L2、L3,因?yàn)閭鞑ヂ窂蕉枷鄬?duì)來說比較平坦,使用逐峰法定位效果較好,因此,在實(shí)際的閃電定位中,因?yàn)闇y站接收電磁波受各種因素影響,使用互相關(guān)效果應(yīng)該更好。
(1)山地地形會(huì)明顯改變切向磁場的波形峰值大小和上升沿時(shí)間,相比較平地,切向磁場波形峰值變化最大的減小了38%,上升沿時(shí)間增加范圍在1~3 μs,峰值到達(dá)時(shí)間也滯后于平地情況,而經(jīng)過帶通濾波器進(jìn)行濾波處理,與過濾前相比,磁場峰值減小,波形變陡,上升沿時(shí)間變短,峰值到達(dá)時(shí)間與過濾前相比幾乎相等不變。
(2)對(duì)于逐峰法,雷擊點(diǎn)在真實(shí)地形路徑下定位誤差分別為351、31、261 m,相比平地明顯增大很多;而對(duì)于互相關(guān)算法,3個(gè)點(diǎn)的定位誤差明顯比逐峰法大很多,真實(shí)地形路徑下定位誤差比平地情況下大,但差距明顯有所減小。
(3)使用一種可以修正閃電定位誤差的算法,即不同空間步長的地形包絡(luò)和時(shí)間補(bǔ)償法,在一定程度上減少地形帶來的閃電定位誤差。在考慮海南地區(qū)的實(shí)際地形情況下,雷擊點(diǎn)L3最大修正了135 m,且不同空間步長在相同地形包絡(luò)法下修正效果不同。對(duì)于逐峰法,閃電L1和L3使用5 km步長的包絡(luò)修訂較好,而閃電L2使用40 km步長的包絡(luò)修訂效果較好,因?yàn)殚W電L1和L3傳播受地形影響較大;使用互相關(guān)法,閃電L1使用10 km包絡(luò)修訂效果較好,而閃電L2和L3使用5 km包絡(luò)修訂效果較好,定位誤差都在百米量級(jí)內(nèi),相比較其他算法閃電定位精度有很大提高。