海南地區山體地形對地閃回擊定位精度的影響

鄒 斌， 崔海華， 張其林*， 李 靜， 戴炳哲，姚年鵬， 蔣 堯， 楊 敢， 羅 華

(1.南京信息工程大學氣象災害教育部重點實驗室， 南京 210044； 2.南京信息工程大學氣候與環境變化國際合作 聯合實驗室， 南京 210044； 3.南京信息工程大學氣象災害預報預警與評估協同創新中心， 南京 210044； 4.南京信息工程大學中國氣象局氣溶膠與云降水重點開放實驗室， 南京 210044； 5.河北省氣象行政技術服務中心， 石家莊 050000)

雷電是一種常見的災害性天氣，每年造成大量的損失[1]。因此，對雷電進行深入的研究，提高閃電探測能力和定位精度，對降低雷電災害損失，保障人類生產活動的有序進行和人員的生命財產安全具有積極的意義[2]。然而現實中除了引雷實驗和極少數地方有引雷塔，可以對現實中雷電定位算法的定位誤差進行評估，對于大部分的雷電很難得知其定位誤差，而隨著計算機技術的發展，模擬算法為評估雷電定位誤差提供了有效途徑。雷電電磁波在近地面傳播過程中受土壤色散、土壤電導率以及高低起伏的不規則地形等因素的影響較大，會導致上升沿時間、電磁波波形峰值等發生改變，從而降低了雷擊點的定位精度[3]。因此，利用二維時域有限差分法(finite difference time domain，FDTD)研究雷電電磁波傳播特性并進行定位誤差分析有重要意義。

在雷電回擊電磁場的研究中，主要包括精確解析計算、近似計算和數值計算3種求解方法。Wait[4]利用等效表面阻抗定義的衰減函數，給出了垂直電場的近似計算公式。Cooray[5]和Rubinstein[6]給出了水平電場的近似計算公式。此后，經過學者們不斷地改進，近似算法被不斷推廣，可以應用到分層土壤、粗糙地表和海陸交接等多種地形環境中[7-10]。然而，解析法在山體等非規則地形中的電磁場計算仍受到較多限制。

中外學者已利用FDTD算法，圍繞土壤電導率[11]、土壤各向異性[12-13]、土壤色散效應[14]和地形[15-16]等多個因素進行了大量的分析討論。張明霞等[17]利用矩量法，討論了矩形山體模型對電磁輻射磁場的影響。Li等[18]在研究地形對雷電定位的影響中發現當傳播路徑中有錐形山體存在時，發現電磁場波形和峰值均會受到嚴重影響，并且信號的時延隨山體高度的增加而增大，隨地面電導率的減小而減小。Li 等[19]通過結合雷擊山頂鐵塔的實測電場波形和FDTD數值模擬，討論了真實不平坦地形的影響，證實了山體地形對電場(15 km觀測點)具有增強作用，結果表明忽略地形起伏的平地假設導致對反演估算的雷電流幅值達到真實值的1.8倍。陶玉郎等[20]利用3D數值模擬研究發現土壤的特征對雷電流沖擊有很大影響。Hou等[21]建立2D FDTD錐形山體模型比較了0.3～50 km處電磁場的變化規律，發現對于垂直電場，近距離受山體的屏蔽作用影響被削弱，而遠距離處的幅值則由于山腳的反射而增強。張金波等[22]利用2D FDTD算法對土壤電導率進行研究，結果表明山體地形坡度對雷電過電壓有顯著影響，尤其在土壤電導率為有限值的情況下，過電壓幅值隨地形坡度的升高而明顯增大。

總的來說，以上研究大部分都局限在孤立山體，涉及實際山體地形的雷電電磁波傳播特性較少。雷電電磁場沿復雜起伏地形傳播時，會以繞射和反射的方式經過山體，會引起電磁場到達測站時間的增加，另外電磁場也會出現強度的變化，從而使根據雷電電磁場波形進行的閃電定位產生誤差。因此如何將地形地貌帶來的誤差進行修正，提高閃電的探測效率和定位精度，是個亟需解決的問題。

為此，在二維柱坐標系下FDTD模型加入了真實地形，研究了真實地形對切向磁場波形以及閃電定位的影響。針對現有的TOA(time of arrival)時差定位算法的缺陷，即忽略了地形對電磁波的到達時間延遲作用，提出了不同步長的地形包絡法，對定位精度進行優化。

1 研究內容和方法

1.1 測站和雷擊點位置

海南地區位于中國南端，主要以山地為主，地勢起伏比較嚴重，山地面積較大，緯度低四面環海，雷暴活動相對比較頻繁，是中國雷電災害最為嚴重的省份之一。以實際安裝在海南省的7個測站位置為基礎，將它們應用于FDTD模式中，7個測站分別為瓊中站、樂東站、三亞站、昌江站、萬寧站、文昌站、澄邁站，具體測站位置如表1所示。在閃電探測網區域選取了3個雷擊點L1、L2、L3，具體位置如表2所示。

表1 測站坐標Table 1 Station coordinates

表2 雷擊點坐標Table 2 Lightning coordinates

圖1為海南地區測站和雷擊點位置，在閃電探測網區域選取了3個雷擊點L1、L2、L3,其中雷擊點L1在中間山地區域，雷擊點L2在東北平原地區，雷擊點L3在西南山地附近的盆地區域，選取這3個雷擊點具有代表性。

1.2 FDTD模式簡介

FDTD是由麥克斯韋旋度方程的微分方程出發的，通過二階中心差分方法，把原來的偏微分方程運算形式轉換成差分方程，從而在一定的時間上和一定體積內，對連續的電磁場數據可以進行抽樣壓縮。

Maxwell方程組的旋度方程組形式為

(1)

式(1)中：E為電場強度；H為磁場強度；D為電位移矢量；ρ為電荷密度；B為磁感應強度；ε為電介質常數；σ為電導率；μ為磁導系數；t為時間。直角坐標系中，式(1)可表示為

(2)

式(2)中：x、y、z分別為直角坐標系的3個方向。

采用磁場波形進行閃電定位，則切向磁場的表達式為[23]

(3)

式(3)中：Δta為時間步長；Δr為水平方向網格長度；Δz為豎直方向格點長度；i為水平方向格點數；j為豎直方向格點數；Hφ為切向磁場；Ez為垂直電場；Er為水平電場；n為第n個計算單元；μ0為真空中的磁導率(即真空中的磁導系數)。

圖2為FDTD模式示意圖，模式中網格寬度Δr×Δz=10 m×10 m，滿足數值穩定性條件，同時也與插值后的地形網格數據相匹配。時間步長Δt=16.7 ns，符合Courant穩定性條件。土壤電導率σg=0.001 S/m，相對介電常數εr=10。空氣電導率σair=0 S/m，相對介電常數εr=1。邊界采用卷積完全匹配層(convolution perfectly matched layer，CPML)吸收邊界條件。閃電回擊通道放置于計算區域的最左側山頂的上面，通道高度H=7 500 m。

z為垂直方向；φ為切向；r為水平方向圖2 模型示意圖Fig.2 Model diagram

1.3 地閃回擊模型

采用MTLE(modified transmission line)源回擊模型，即回擊通道中的回擊電流以指數的形式沿著通道高度在衰減，回擊通道中t時刻在高度z′處的電流分布可表示為

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

(4)

式(4)中：e-z′/λ為數值衰減系數；i為電流強度；z′為高度；λ為衰減因子，取2 000 m；ν為回擊速度，取1.5×108m/s。

使用Heidler雙指數函數[24]來模擬繼后回擊下雷電流通道底部的基電流，如式(5)所示。典型的繼后回擊的電流參數如表3所示。

表3 繼后回擊電流參數典型值Table 3 Typical value of subsequent return stroke current parameters

(5)

(6)

(7)

式中：i01和i02分別為擊穿電流和電暈電流峰值；η1、η2分別為擊穿電流和電暈電流峰值的修正因子；τ11和τ21分別為擊穿電流和電暈電流的波形上升時間；τ12和τ22分別為擊穿電流和電暈電流的波形下降時間。

1.4 波形到達時間計算方法

時差法是利用電磁信號到達不同測站的時間差進行定位，結合高精度的全球定位系統(global positioning system，GPS)同步時間測量，能夠實現對輻射源點的準確定位。在利用TOA 算法計算輻射原點的過程中，首先對信號進行波形識別篩選出不同類型的閃電事件，同時對不同的同步脈沖進行分組，再計算出同步信號到達不同測站的時間差，然后假定信號沿直線距離以光速向測站傳播，求解非線性方程組獲得輻射源點位置的初值，并借助非線性最小二乘擬合算法對初解進行優化，最終得到輻射源點的位置。圖3為到達時間示意圖。

Upeak為閃電峰值電壓；tt為閃電上升沿開始時間； tpeak為波形峰值時間，即逐峰法定義的時間圖3 到達時間示意圖Fig.3 Diagram of arrival time

為了求解出地閃回擊點的空間坐標，需要先得到信號到達不同測站的時間差。在時間差的計算中，可以通過互相關方法直接計算出兩兩測站之間的時間差，也可以以電磁波峰值的到達時間為準，先計算出單站信號的到達時間，再與參考站點的到達時間相減獲得兩個測站之間的時間差。不同方法計算得到時間差可能存在一定的差異，為了對比不同方法的定位效果，選取互相關，逐峰法這2種方法的定位結果進行比較。

1.5 地形包絡法

雖然對閃電定位普遍采用的是TOA時差法，但顯然這種算法并沒有考慮實際地形對閃電電磁波產生的路徑延長問題，即信號到達每個測站的時間都會比理想情況下略長，因此采用不同步長的地形包絡法來近似估計地形引起的延長時間。

如圖4所示，地形包絡即用連續線段將地形包圍起來。考慮到所有測站和雷擊點的距離大小，每隔5、10、20、40 km步長距離取每段地形的最高點并連線形成包絡線。將這些地形包絡線長度作為電磁波在傳播中的近似傳播距離，反映出信號傳播路徑隨山體高度的變化。然后用包絡線長度減去水平距離后除以電磁波的傳播速度作為電磁波到達時間的誤差訂正值。假設真實路徑下用包絡法得到測站至雷擊點路徑長度為l1，水平路徑為l2，可以得到誤差訂正時間Δtb的計算公式為

(8)

式(8)中：c為光速。

將原各個測站的信號到達時間分別減去對應的誤差訂正時間，重新進行雷擊點定位。

1.6 時間補償法

忽略測站GPS 計時精度的誤差，信號從雷擊點向測站傳播的時間t，會因為地形、電導率等的影響而變化，與同等水平傳播距離的理想平坦地表時的傳播時間t存在一定的偏差，即

t=tflat+Δtc

(9)

式(9)中：tflat為信號以光速傳播過平坦地表的時間；Δtc為真實崎嶇地形情況與理想平坦地表情況的傳播時間偏差。

進行時間補償時，將信號以光速沿紅色折線和黑色直線傳播的時間差作為對式(8)中Δt的估算，減少傳播時間t與理想情況下的tflat的偏差。有必要時，考慮進行多次補償，提高最終定位結果的準確性。

1.7 雷擊點到測站真實地形

圖5～圖7為L1、L2、L3雷擊點到測站的真實地形高程圖，可以看出，雷擊點L1雷電電磁波傳播路徑的總體趨勢為由高到低，雷擊點在山地上，傳播路徑總體上比較復雜，地形起伏較大；雷擊點L2在平原地區，傳播的總體趨勢為由低到高，總體上受山地影響較小；雷擊點L3在萬寧測站附近，雷擊點海拔高度低，但其四周地勢都比雷擊點高，在山地的低洼區域，受山地影響更明顯。

圖5 雷擊點L1到測站地形圖Fig.5 Topographic map of lightning L1 to the station

圖7 雷擊點L3到測站地形圖Fig.7 Topographic map of lightning L3 to the station

2 模擬結果與討論

2.1 平地與真實的波形對比

圖8～圖10為雷擊點L1、L2、L3到每個測站平地與真實地形的波形對比。藍色線和紅色線分別為平地和真實地形的模擬結果。山地地形對波形的上升沿時間和波形峰值影響較大。上升沿時間增加范圍在1～5 μs，主要是電磁波在傳播中受土壤有限電導率和色散效應影響，高頻電磁場存在更大的衰減。同樣峰值時間也都是滯后于平地的，但峰值的大小變化卻不盡相同。

圖9 雷擊點L2到測站磁場波形Fig.9 Magnetic field waveform from lightning L2 to the station

圖10 雷擊點L3到測站磁場波形Fig.10 Magnetic field waveform from lightning L3 to the station

對于雷擊點L1，瓊中站真實地形下峰值明顯比平地情況大很多，峰值時間較平地情況也提前，而其他6個測站真實地形下的峰值對比與平地都略有減少，真實地形下峰值時間明顯滯后于平地傳播的峰值，這主要受山地影響，雷電高頻電磁波衰減較大，而瓊中站距雷擊點距離較近，受山地影響較小，這也說明了地形對電磁波影響的復雜性。

對于雷擊點L2，7個測站真實地形下的峰值都略小于平地情況，峰值時間也都滯后，三亞站和昌江站滯后比較明顯，峰值減小也較明顯，其他5個測站真實地形下的電磁波波形及到達時間基本上一致，這主要受傳播路徑影響，閃電從雷擊點L2到三亞站和昌江站受山地影響較大，而從雷擊點L2到另外5個測站地形相對比較平坦。

對于雷擊點L3, 6個測站真實地形下的峰值都明顯小于平地情況，峰值時間除樂東站外也都明顯滯后，主要是樂東站從雷擊點L3到測站傳播路徑比較平坦，和雷擊點L2一致。

2.2 真實地形下濾波前后波形對比

考慮到雷電探測傳感器探測都具有特定的帶寬，實際接收的波形已被過濾。在此，使用帶通濾波器，對2-D FDTD模型的模擬磁場波形進行過濾，以分析帶寬對雷電磁場的影響。圖11～圖13分別為雷擊點 L1～L3 在真實地形下，各測站接收的磁場波形過濾前后對比，頻率范圍為10～500 kHz。與過濾前相比，磁場峰值減小，波形變陡，上升沿時間變短，而峰值到達時間與過濾前相比幾乎相等，這樣對于逐峰法定位來說，誤差幾乎相似。

圖11 雷擊點L1真實地形下到測站濾波前后磁場波形Fig.11 Waveforms of magnetic field before and after filtering at lightning point L1 in real terrain

圖12 雷擊點L2真實地形下到測站濾波前后磁場波形Fig.12 Waveforms of magnetic field before and after filtering at lightning point L2 in real terrain

圖13 雷擊點L3真實地形下到測站濾波前后磁場波形Fig.13 Waveforms of magnetic field before and after filtering at lightning point L3 in real terrain

2.3 定位點的偏差及修正

圖14為在逐峰法波形到達時間算法下3個雷擊點在平地、真實地形、以及不同步長地形包絡法修正后的定位誤差對比圖。可以看出，雷擊點L2的誤差較小，主要是雷擊點L2到測站的傳播受地形影響較小，且用包絡法進行修正，隨著包絡空間步長的增加，定位誤差在減小；對于雷擊點L1和L3，定位誤差明顯大于雷擊點L2的定位誤差，對于平坦地表，定位誤差只有十幾米，定位效果較好，而使用包絡法和時間補償法進行修正，定位誤差明顯減小，且包絡空間步長越小，定位誤差越小。

圖14 逐峰法定位誤差Fig.14 Location error by peak method

圖15為在互相關法波形到達時間算法下3個雷擊點在平地、真實地形、以及不同步長地形包絡法修正后的定位誤差對比。可以看出，雷擊點L1的定位誤差相比較其他兩個雷擊點要小很多，雷擊點L3使用六站進行定位相比較雷擊點L1的7站定位，定位誤差明顯大了很多；對于3個雷擊點，使用包絡法和時間補償法進行修正，定位誤差明顯減小，且包絡空間步長越小，定位誤差越小。

圖15 互相關法定位誤差Fig.15 Location error by cross-correlation method

對比3個雷擊點在兩種定位算法下的定位誤差，從圖14、圖15可以看出：對于平坦地表，雷擊點逐峰法下的定位效果更好，定位誤差都只有幾十米，而對于雷擊點L1，因為傳播受地形影響較大，使用互相關算法定位誤差較好，而雷擊點L2、L3，因為傳播路徑都相對來說比較平坦，使用逐峰法定位效果較好，因此，在實際的閃電定位中，因為測站接收電磁波受各種因素影響，使用互相關效果應該更好。

3 結論

(1)山地地形會明顯改變切向磁場的波形峰值大小和上升沿時間，相比較平地，切向磁場波形峰值變化最大的減小了38%，上升沿時間增加范圍在1～3 μs，峰值到達時間也滯后于平地情況，而經過帶通濾波器進行濾波處理，與過濾前相比，磁場峰值減小，波形變陡，上升沿時間變短，峰值到達時間與過濾前相比幾乎相等不變。

(2)對于逐峰法，雷擊點在真實地形路徑下定位誤差分別為351、31、261 m，相比平地明顯增大很多；而對于互相關算法，3個點的定位誤差明顯比逐峰法大很多，真實地形路徑下定位誤差比平地情況下大，但差距明顯有所減小。

(3)使用一種可以修正閃電定位誤差的算法，即不同空間步長的地形包絡和時間補償法，在一定程度上減少地形帶來的閃電定位誤差。在考慮海南地區的實際地形情況下，雷擊點L3最大修正了135 m，且不同空間步長在相同地形包絡法下修正效果不同。對于逐峰法，閃電L1和L3使用5 km步長的包絡修訂較好，而閃電L2使用40 km步長的包絡修訂效果較好，因為閃電L1和L3傳播受地形影響較大；使用互相關法，閃電L1使用10 km包絡修訂效果較好，而閃電L2和L3使用5 km包絡修訂效果較好，定位誤差都在百米量級內，相比較其他算法閃電定位精度有很大提高。