闪电辐射源多脉冲匹配定位方法

王彦辉， 闵应昌， 刘亚栎， 刘银萍

(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)从整体定位效果来看，本文中提出的闪电辐射源多脉冲匹配定位方法效果较好，获得了较清晰的闪电发展通道，定位出了单脉冲匹配定位方法未定位出的一次闪电的多个放电发展过程及闪电的多个放电通道形状。

闪电辐射信号的预处理及同一辐射源的脉冲匹配不仅对闪电三维定位具有重要意义，也能够帮助研究同一个辐射源发出的电磁脉冲到达不同测站的幅度和宽度差异，对雷电电磁脉冲防护具有一定的意义。在闪电发展过程中，多个雷电分支通道同时发展使闪电辐射源天线接收信号发生错乱，滑动窗口内不同分支通道的辐射源脉冲杂糅在一起，单个辐射源脉冲幅度较弱或者多个辐射源脉冲信号在时间上交错，这种情况下的相同辐射源脉冲匹配效果不够理想，该问题值得进一步研究，这将是下一步重点研究的内容。