閃電輻射源多脈沖匹配定位方法

王彥輝， 閔應昌， 劉亞櫟， 劉銀萍

(1.南京信息工程大學氣象災害預報預警與評估協同創新中心， 南京 210044； 2.南京信息工程大學大氣物理學院， 南京 210044； 3.中國科學院西北生態環境資源研究院高原大氣物理研究室， 蘭州 730000)

現代信號處理技術和高速數據采集技術的應用提升了雷電探測的能力和水平，閃電甚高頻(very high frequency, VHF)輻射源定位技術能夠實現閃電通道發展的時空演變過程分析，對閃電物理機理研究及雷電電磁脈沖防護具有重要的意義[1-4]。

在整個閃電定位系統中最核心的技術是閃電定位算法[5]。定位過程主要包括波形預處理、尋找脈沖、匹配脈沖以及到達時差方法解算最優解[6]。利用閃電輻射源的脈沖信號到達不同測站的到達時間差(time difference of arrival, TDOA)進行閃電通道定位的關鍵就是準確地找到輻射源脈沖，并將不同測站的同一個放電事件產生的脈沖匹配在一起，錯誤的脈沖匹配會得到不準確甚至錯誤的定位結果[6]。因此，將閃電產生的相同輻射源脈沖準確匹配獲得脈沖的到達時間差對定位結果至關重要。Shao等[7]將原始波形經過Hilbert變換后獲取歸一化的功率波形，再從歸一化的功率波形中提取脈沖的峰值信息,并進一步判斷兩兩脈沖峰的時間差是否小于兩站之間的光程差(距離/光速)，如果兩站之間的時間差小于光程差，那么說明這兩個脈沖的匹配是合理的。Shi等[8]在該方法基礎上做了改進，在波形匹配后，首先，判斷兩兩脈沖峰的時間差是否小于兩站之間的光程差，其次，為保證兩個測站的脈沖匹配唯一性，其脈沖到該段波形最大幅值的高度差應最接近。基于原始波形轉換的歸一化功率波形進行尋峰運算，獲得滿足條件的脈沖峰值幅度和峰值時間；基于兩個子站間脈沖的單脈沖峰值時間差限制條件和幅度相近條件，獲得匹配的脈沖峰值時間，并結合子站位置信息，組成到達時間非線性方程組，使用最優化算法獲得最優解。該方法由于僅僅采用單脈沖峰值時間和幅度進行匹配，在脈沖信號比較豐富的情況下，會出現較多的錯誤匹配，導致較低的定位計算效率；同時，采用單脈沖相似幅度的匹配，也會導致部分放電事件不能有效匹配的情況，需要進一步改進算法。

中外一些定位系統在利用閃電低頻/甚低頻(low frequency/very low frequency, LF/VLF)電場脈沖信號的定位研究中采用了互相關(cross-correlation)算法來實現閃電信號到達各個測站的時間估計[9-11]。

劉波等[12]將頻域時間反轉(time reversal, TR)聚焦算法應用到閃電輻射源定位中，李書磊等[2]將正交傳播算子(orthogonal propagator method, OPM)方法應用于閃電寬帶輻射源定位。

為了增加脈沖信號的匹配成功率和提高脈沖峰值時間的精度，從而提高定位精度和精細化程度，Fan等[13]將經驗模態分解(empirical mode decomposition, EMD)信號處理技術應用到閃電波形處理中，實現多站閃電電場波形的低頻濾波和高頻噪聲降低，并基于處理后的波形，通過逐步縮減窗口進行波形相關匹配后轉化形成匹配的功率波形，尋找脈沖獲得匹配的脈沖峰值時間，并采用到達時差方法解算最優解。但由于EMD固有的模式混疊現象，采用EMD方法處理閃電電場信號后，利用波形互相關進行脈沖匹配時仍存在多站波形匹配不正確，導致脈沖信息提取不準確[14]。

基于上述原因，采用EMD方法對原始信號進行分解后再部分合成，并對合成后的波形信號進行歸一化處理，獲得較干凈和平穩的閃電電場信號波形。提出了四脈沖序列的高度比例特征和間隔特征作為脈沖匹配條件，對相同閃電輻射源脈沖進行匹配，并基于到達時間差對閃電輻射源進行三維定位，描繪出閃電通道的發展路徑并計算出各個輻射源發生的時間。

1 閃電輻射源定位方法

利用到達時間差對閃電輻射源進行定位，4個及4個以上測站的同一輻射源的到達時間差定位出一個輻射源的空間位置及其發生的時間，多個輻射源點能較好地描繪出閃電通道的基本形狀和放電通道的發展過程。對于同一個閃電輻射源，(xi,yi,zi)和ti表示第i個測站的三維坐標和接收到輻射源信號的時間，(x,y,z)和t表示輻射源空間位置和發生時間，c為光速，則第i個測站接收到閃電輻射源的時間滿足關系式為

(1)

每個測站都可得到式(1)所示的一個關系式，4個及4個以上測站得到一個方程組，解方程組求得輻射源的空間位置(x,y,z)和輻射源發生的時間t。

2 數據預處理方法

經驗模態分解(EMD)是一種新的處理非平穩信號的方法[15]，能將復雜信號分解為有限個本征模函數(intrinsic mode function, IMF)。EMD適合于非線性、非平穩信號的分析，也適合于線性、平穩信號的分析。經驗模態分解方法的實質是通過特征時間尺度來識別信號中所內含的所有振動模態，為了從原始信號中分解出本征模函數，經驗模態分解方法過程如下：找到測站i的原始信號Xi(t)所有的極值點，用3次樣條曲線分別擬合出上下極值點的包絡線，原始信號減去兩條包絡線的平均值得到c1，根據預設判據判斷c1是否為IMF，如果不是，則以c1代替Xi(t)，重復以上步驟直到c1滿足判據，則c1就是需要提取的IMF分量，每得到一階IMF分量，就從原信號中扣除該分量，重復以上步驟，直到信號最后剩余部分rn滿足預設判據，剩余部分rn通常是直流分量、單調分量和低頻周期分量或3種分量的疊加。

利用EMD方法將原始信號Xi(t)分解成一系列IMF分量以及剩余部分rn的線性疊加，可以表述為

(2)

根據EMD的信號分解原理，EMD方法將組成原始信號的各尺度分量不斷從高頻到低頻進行提取，其中殘差分量rn通常是單調序列或者常值序列或者低頻周期序列或三者的疊加，也稱為趨勢分量，利用EMD方法對原始信號進行分解之后，再由去除趨勢分量之外的其他分量合成得到關系式

(3)

的信號，獲得較干凈的信號xi(t)，從而獲得非線性和非平穩信號中的輻射源脈沖。

閃電輻射源的寬帶電場原始信號經過EMD方法分解后再合成，由于不同的測站對信號的放大能力不完全相同，以及測站與輻射源的距離也會影響測站接收到的輻射源信號強度，因此需要對EMD方法處理之后的閃電電場信號進行歸一化處理，以使得不同測站的同一輻射源脈沖信號強度更趨于一致，并根據波形對應的GPS(global positioning system)時間實現信號波形的整體對齊。

3 脈沖匹配方法

本文中使用的數據對應的采樣率為20 MS/s，單次閃電的數據記錄時長為1.2 s，即每個測站的一次閃電信號由2 400萬個數據點組成，相鄰兩個數據點的時間間隔為50 ns。

信號的局部極大值與其左右極小值的高度差為該脈沖相對高度，以預設的脈沖最小相對高度和脈沖最小間隔為條件提取各個測站的脈沖。以明德測站為主站，明德測站的脈沖序列為樣板。由于本文中使用的數據其探測站間最遠距離不超過15 km，即相同輻射源脈沖到達各個測站的光程差在±0.05 ms范圍內，本文以各個測站0.15 ms信號波段的GPS時間相同確定各測站滑動窗口的位置，滑動窗口的數據長度為3 000個數據點，窗口時間長度為0.15 ms，即在滑動窗口內匹配相同輻射源的脈沖信號。各測站完成當前窗口的脈沖匹配，滑動窗口移動到主站下一個脈沖所在位置，以主站的該脈沖時間為中心在各測站截取GPS時間相同的0.15 ms信號波段以獲得新的滑動窗口。

對于同一個輻射源的脈沖，歸一化處理之后各測站脈沖的幅度仍然存在差異，對所有輻射源脈沖又不能明確這種差異在某一特定的范圍。本文中采取的方法是：在各測站的滑動窗口內，采用四個脈沖序列的高度比例特征和間隔特征作為脈沖匹配的條件，每一次匹配出四組相同輻射源脈沖。

如圖1(a)所示，主站的第j、k、m、n4個脈沖在主站的當前滑動窗口內；如圖1(b)所示，第i個測站的第ij、ik、im、in4個脈沖在第i個測站的當前滑動窗口內。主站的第j個脈沖高度為基線，第k個脈沖到基線的高度為hk，第m個脈沖到基線的高度為hm，第n個脈沖到基線的高度為hn；同樣的，第i個測站的第ij個脈沖高度為基線，第ik個脈沖到基線的高度為Hk，第im個脈沖到基線的高度為Hm，第in個脈沖到基線的高度為Hn，hk與Hk的比例為λk，hm與Hm的比例為λm，hn與Hn的比例為λn，當λk、λm、λn的兩兩之差在一定的容差范圍內，即滿足4個脈沖序列的高度比例特征；同時，主站的第j個脈沖與第k個脈沖的距離為xjk，第k個脈沖與第m個脈沖的距離為xkm，第m個脈沖與第n個脈沖的距離為xmn，第i個測站的第ij個脈沖與第ik個脈沖的距離為Xijik，第ik個脈沖與第im個脈沖的距離為Xikim，第im個脈沖與第in個脈沖的距離為Ximin，當xjk與Xijik之差在一定容差范圍內，xkm與Xikim之差在一定容差范圍內，xmn與Ximin之差在一定容差范圍內，即滿足4個脈沖序列的間隔特征。

圖1 脈沖分布示意圖Fig.1 Schematic diagram of pulses distribution

滿足上述的脈沖匹配條件，則第i個測站的第ij、ik、im、in4個脈沖分別與主站的第j、k、m、n4個脈沖匹配成功；下一個測站與主站在當前滑動窗口內進行脈沖匹配，直到最后一個測站完成當前窗口內的脈沖匹配，滑動窗口向后滑動，進行下一組脈沖的匹配，重復上述步驟，直到滑動窗口移動到信號末尾。

如圖1所示，當hk、hm、hn、Hk、Hm、Hn均為非零時，設

(4)

設α為高度比例差的容差范圍，β為脈沖間隔差的容差范圍，本文中設定α=50%，β=0.75 μs。若λk、λm、λn滿足關系式

(5)

則滿足4個脈沖序列的高度比例特征；若xjk、xkm、xmn和Xijik、Xikim、Ximin滿足

(6)

則滿足4個脈沖序列的間隔特征。同時滿足4個脈沖序列的高度比例特征和間隔特征，則完成第i個測站的第ij、ik、im、in4個脈沖分別與主站的第j、k、m、n4個脈沖的匹配。當這4個脈沖序列的某一組脈沖相對于第一組脈沖高度為零或趨近于零時，例如當hm或Hm為零或者趨近于零時，則不計算該組脈沖的高度比例λm，對于該組脈沖只需滿足四脈沖序列的間隔特征即可認為該組輻射源脈沖匹配完成。

本文中將利用一次探測到的閃電寬帶電場信號驗證所提出的多脈沖匹配方法的效果，用于進一步補充說明本文的多脈沖匹配方法。

4 實驗分析

4.1 數據采集

閃電甚高頻(very high frequency,VHF)輻射源定位系統[3]主要由VHF天線、寬帶電場變化測量天線、帶通濾波器、對數放大器、高速A/D數據采集卡、高精度時鐘(時間精度50 ns)、處理器及數據無線傳輸模塊組成。VHF接收系統中心頻率270 MHz，帶寬為6 MHz，高速A/D數據采集卡的采樣率為20 MS/s，寬帶電場系統帶寬0～10 M。為把接收到的閃電寬帶電場信號波形不失真的放大，該系統采用了對數功放電路。閃電輻射的脈沖信號幅度峰值和峰值時間由數字化模塊記錄緩存，記錄長度根據基線長度及處理數據的分辨率而定，本文采用的單次數據記錄時長為1.2 s,系統觸發及數據記錄的時間與GPS高精度時鐘同步。信號的噪聲電平由同步觸發閾值電路控制，處理器接收、處理和記錄輻射源脈沖峰值絕對時間信息，并實時通過無線數據傳輸系統將時間信息發送到中心站進行時差計算及實時處理顯示。

采用和分析的數據來自于在青海省大通縣進行的雷電綜合觀測實驗，閃電發生在2011年8月15日0：2：41。大通縣位于青海省東北部，平均海拔約2 600 m，是一個多山地區，地形復雜、氣候垂直差異明顯，由于地處內陸高原地區，海拔較高，受高原天氣系統與西風帶天氣系統相互作用，加之河谷地形和下墊面性質的復雜性，造成短時強降水較多，并常伴有雷暴和冰雹等強對流天氣。

如圖2所示，雷電綜合觀測實驗由7個觀測站組成，其中以明德站為中心，各測站分布半徑約為8 km。基于時差法的閃電VHF輻射源三維定位系統中,觀測試驗的7個測站的寬帶電場變化測量系統通過無線寬帶通信系統連接組成同步觀測網絡，采用主站控制方式，并用GPS同步的高精度時鐘(50 ns)記錄觸發時間，從各測站記錄的寬帶電場變化波形中提取并匹配相同輻射源的脈沖，計算各測站相同輻射源的脈沖信號到達時間并利用時差法定位技術實現閃電輻射源的三維定位，精確描繪閃電放電通道，其中，明德、藥草、新莊、苗圃、極樂、斜溝、良教這7個測站的寬帶電場變化測量系統數據采樣率為20 MS/s,單次閃電數據記錄時長為1.2 s。

根據圖2所示各測站的經緯度和海拔高度，以明德站為中心站，確定各測站相對于中心站的坐標，進行閃電輻射源定位的輻射源空間位置也是以明德站為三維坐標原點，各個測站的相對坐標分別為明德(0，0，0)，藥草(3 178.44，2 921.02，131.538)，新莊(-3 355.01，5 572.89，64.680 7)，苗圃(-6 013.87，2 442.51，44.694 9)，極樂(-8 325.07，-2 473.84，196.083)，斜溝(-4 165.98，-4 407.31，242.115)，良教(883.798，-4 003.41，-9.318 57)，其中，各測站的相對坐標以米為單位。

圖2 青海省大通縣測站地理位置Fig.2 Location of stations in Datong County, Qinghai Province

4.2 數據預處理結果

如圖3(a)所示，7個測站單次閃電的數據記錄時長為1.2 s，數據采樣率為20 MS/s，其中，明德站的觸發前數據占比為28%，藥草、新莊、苗圃、良教4個測站的觸發前數據占比為30%，極樂、斜溝兩個測站的觸發前數據占比為35%。

使用EMD方法對原始閃電數據分解再部分合成，從原始信號中提取出圖3(b)所示的信號分量，去除了圖3(c)所示的信號分量，去除的信號分量包括直流分量、單調分量和低頻周期分量。每個測站原始信號經過EMD方法分解再部分合成，再對合成的信號進行歸一化處理，得到的信號最大值為1,最小值為-1，并根據閃電信號波形對應的GPS時間完成信號波形的整體對齊，如圖3(d)所示。

圖3 數據預處理結果Fig.3 Results of data preprocessing

4.3 輻射源脈沖匹配結果

如圖4所示為該次閃電的其中一個四脈沖序列匹配過程，其中，圖4(a)所示為7個閃電探測站在相同GPS時間采集到的閃電電場波形；圖4(b)是圖4(a)中標出信號的幅度放大波形，明德(主站)的四脈沖序列與藥草(子站)的四脈沖序列進行匹配，分別用紅虛線和藍虛線標出明德測站與藥草測站的四脈沖序列；圖4(c)為該四脈沖序列的最終匹配結果，可以看到，7個測站中有6個測站實現了脈沖匹配，實現6站同步，可以計算出相應的閃電輻射源位置。

圖4 閃電四脈沖序列匹配過程Fig.4 Four pulses sequence matching process of lightning

如圖5(a)所示為稀疏的大脈沖的匹配結果，如圖5(b)所示為密集的小脈沖的匹配結果，對于稀疏的大脈沖和密集的小脈沖，使用本文方法能獲得較好的脈沖匹配效果。以明德站的脈沖序列為樣板，根據滑動窗口內4個脈沖序列的高度比例特征和間隔特征，一次匹配出四組相同輻射源的脈沖，滑動窗口再次向后移動到主站的下一個脈沖所在位置，完成其他輻射源脈沖的匹配，并計算出各個測站相同輻射源脈沖的到達時間。如圖5(c)所示為第884組輻射源脈沖匹配時滑動窗口的信號波形，可以看到，各個測站的脈沖特征差異較大，多個輻射源脈沖信號在時間上交錯，可能是由于閃電的多個分支通道同時發展，使得輻射源天線在較短時間內同時接收來自多個閃電通道的輻射源脈沖信號。如圖5(d)所示為圖5(c)所示滑動窗口內的脈沖匹配結果，其中極樂測站未匹配出輻射源脈沖的原因是該處接收到的輻射源信號較弱，干擾噪聲大于輻射源脈沖信號，而良教測站未匹配出輻射源脈沖的原因是四脈沖序列中存在可能其他分支通道的輻射源脈沖，且不同分支通道相距較遠。對于2011年8月15日0:2:41的該次閃電共匹配出1 897組相同輻射源脈沖。

4.4 定位結果對比與閃電過程分析

如圖6所示，圖6(a)和圖6(b)分別是本文中多脈沖匹配方法與傳統的單脈沖匹配方法的定位結果。圖6中閃電VHF輻射源3D定位系統的顯示界面分為5個視圖，從上到下、從左到右依次為：輻射源高度隨時間變化圖、東西方向的垂直投影圖、輻射源數量隨高度變化圖、水平投影圖、南北方向的垂直投影圖。其中，輻射源數量隨高度變化圖的縱坐標表示輻射源高度，橫坐標表示相應高度的輻射源數量，從圖6可以看出閃電輻射源的高度分布特征；垂直投影和水平投影的坐標是以明德測站為原點，表示輻射源相對于明德測站的三維位置。本文中的閃電信號預處理方法及脈沖匹配方法實現了較多相同輻射源的脈沖匹配，利用輻射源脈沖的到達時間差定位方法計算輻射源三維坐標及其發生時間，獲得了較清晰的閃電發展通道，定位出了一次閃電的多個放電發展過程，其中一些放電過程是單脈沖匹配方法未定位出的，定位結果如圖6(a)所示，根據顏色區分不同放電過程發生的時間及通道長度。

單脈沖匹配方法未將原始信號中的直流分量、單調分量和低頻周期分量去除，并且是基于多個測站的單脈沖匹配，依據脈沖到達時的GPS時間最接近匹配相同輻射源的脈沖，其定位結果如圖6(b)所示，部分通道定位效果不明顯，甚至不能定位出閃電的部分放電通道。從整體定位效果來看，本文提出的輻射源脈沖匹配定位方法效果更好。

如圖6(a)所示，閃電起始區域的東南角的閃電通道的輻射源數目較少，該時段的閃電信號波形如圖5(c)所示，部分測站接收到的輻射源脈沖信號較弱，小于該處的噪聲干擾，或是四脈沖序列中存在可能其他分支通道的輻射源脈沖，不同放電通道距離較遠，且輻射源信號天線同時接收到來自多個放電通道的輻射源脈沖信號，使得多個放電通道的多個輻射源脈沖信號在時間上交錯，在這種情況下四脈沖序列同時匹配的匹配效果不夠理想，而傳統的單脈沖匹配方法則具有一些優勢。

圖5 脈沖匹配結果Fig.5 Results of pulses matching

閃電發生于2011年8月15日0：2：41，利用相同輻射源到達時間差對閃電輻射源進行定位，顏色顯示先后順序為藍、青、綠、黃、紅、洋紅，顏色代表閃電輻射源發生的時間先后，如圖6(a)所示。閃電起始位置在明德站的西北方向，相對坐標大概(-2 km,8 km)處，第一個閃電輻射源發生在41.715 000 s。從41.840 s開始，閃電放電通道從閃電起始區域向西南方向發展，該通道在41.939 s結束，通道發展時間約100 ms，通道長度約7.8 km，通道發展速度約7.8×104m/s。該通道從42.169 35 s開始出現第一次再放電過程，通道長度約6 km，該次再放電過程發展較快，持續時間約7 ms，速度約8.6×105m/s，與曹冬杰等[16]的研究結果在同一個量級。

圖6 閃電輻射源定位結果對比Fig.6 Results of lightning radiation localization

從閃電起始區域向西北方向發展的閃電通道中輻射源脈沖數量較少。從42.326 6 s開始該通道再次放電，從輻射源定位結果來看，該通道再放電過程是沿著原來的通道快速向閃電起始區域發展，其發展方向與該通道第一次放電過程的發展方向相反，其通道長度約3.6 km，持續時間約3.5 ms，速度約為1×106m/s，曹冬杰等[16]對云閃的反沖流光速度估計為106～107m/s，與這里的結果基本一致。

在閃電起始區域西北方向的通道再次放電結束后間隔約1.58 ms，從閃電起始區域向西南方向發展的閃電通道出現第二次再放電過程，沿著第一次放電的通道向西南方向發展，該通道第二次再放電過程的放電通道長度約4 km，持續時間約16 ms，通道發展二維速度約2.5×105m/s，與Liu等[17]和張泉等[18]的研究結果基本一致。最后一個閃電輻射源發生在42.347 675 s，整個閃電過程持續時間約632.675 ms。

5 結論

采用經驗模態分解(EMD)方法對原始信號進行分解后再部分合成，并對合成后的波形信號進行歸一化處理，獲得較干凈和平穩的閃電電場信號波形。在采用EMD方法處理閃電寬帶電場信號之后，不采用波形互相關算法匹配脈沖，而是在單脈沖峰值時間差限制條件和幅度相近條件匹配閃電輻射脈沖的基礎上，提出了四脈沖序列的高度比例特征和間隔特征作為脈沖匹配條件，對相同閃電輻射源脈沖進行匹配。并且，利用一次探測到的閃電寬帶電場信號，以本文方法與單脈沖峰值時間差限制條件和幅度相近條件匹配閃電輻射脈沖的方法進行了脈沖匹配效果的比較，所提出的多脈沖匹配方法具有較好的自適應性和準確性，實現了更多組相同輻射源的脈沖匹配。對于本文中使用的這一次閃電數據，本文中所提出的多脈沖匹配方法最終定位出了1 897個閃電輻射源，而單脈沖匹配方法最終定位出了727個閃電輻射源。

通過上述閃電輻射源多脈沖匹配定位方法的提出和驗證，得到以下結論。

(1)通過EMD方法對閃電輻射信號進行分解后再部分合成，去除了原始信號中的單調分量、直流分量和低頻周期分量，可以使信號波形更干凈和平穩，便于后續的脈沖提取和匹配。

(2)對分解再合成的信號進行歸一化處理，可以使不同測站的同一輻射源信號強度更趨于一致。

(3)以滑動窗口內4個脈沖序列的高度比例特征和間隔特征作為各測站間相同輻射源的脈沖匹配條件，每一組脈沖匹配時四脈沖序列的高度比例特征和間隔特征不同，該方法具有較好的自適應性和準確性，能實現較多組相同輻射源的脈沖匹配。

(4)從整體定位效果來看，本文中提出的閃電輻射源多脈沖匹配定位方法效果較好，獲得了較清晰的閃電發展通道，定位出了單脈沖匹配定位方法未定位出的一次閃電的多個放電發展過程及閃電的多個放電通道形狀。

閃電輻射信號的預處理及同一輻射源的脈沖匹配不僅對閃電三維定位具有重要意義，也能夠幫助研究同一個輻射源發出的電磁脈沖到達不同測站的幅度和寬度差異，對雷電電磁脈沖防護具有一定的意義。在閃電發展過程中，多個雷電分支通道同時發展使閃電輻射源天線接收信號發生錯亂，滑動窗口內不同分支通道的輻射源脈沖雜糅在一起，單個輻射源脈沖幅度較弱或者多個輻射源脈沖信號在時間上交錯，這種情況下的相同輻射源脈沖匹配效果不夠理想，該問題值得進一步研究，這將是下一步重點研究的內容。