CN102298097B - 一种估计云闪雷电脉冲信号到达时间差的方法 - Google Patents

Abstract

本发明公开了一种估计云闪雷电脉冲信号到达时间差的方法，步骤为：雷电脉冲数据块拼接、波形还原、形成雷电脉冲簇；输入时间区间，提取各站的雷电脉冲簇；确定参考探测站；统计A站雷电脉冲簇中的脉冲个数；雷电脉冲关联与TDOA计算；形成用户指定时间区间内的所有同源脉冲对队列和TDOA队列；改变时间区间，计算并形成新的同源脉冲对队列和TDOA队列。本发明专门针对云闪雷电脉冲信号的特点进行设计，给出了波形关联与同源性判定、到达时间差估计的完整过程，可胜任白噪声和色噪声环境下的TDOA估计要求。本发明方法无需遍历两个探测站接收到的所有雷电脉冲，鲁棒性强，简单可靠，计算结果精度高，易于工程实现。

Description

一种估计云闪雷电脉冲信号到达时间差的方法

技术领域

本发明属于雷电监测预警技术研究领域，涉及一种在由多个探测站和一个中心站组成的云闪雷电探测定位系统中，在中心站精确计算雷电脉冲信号到达任意两个探测站的时间差(TDOA)的方法。

背景技术

基于信号来波方向(DOA)对雷电辐射源进行定位，其定位精度与空间分辨力将随着雷电辐射源与探测站之间距离的增大而降低。相比之下，基于信号到达时间差(TimeDifference Of Arrival，TDOA)的雷电辐射源定位技术，其定位精度与空间分辨力则将随着辐射源与探测站间距的增加而得到提高(当定位双曲面相互正交时为最佳)，但当辐射源非常靠近某个探测站时，TDOA定位的空间分辨力也将急剧下降。

获取辐射源发射信号到达两个探测站的时间差有两种方法：(1)将两个接收站得到的信号到达时间(TOA)测量值作差；(2)利用互相关技术，计算一个探测站的接收信号与另一个探测站的接收信号之间的互相关函数，然后进行峰值搜索。前一种方法不仅要求各探测站之间保持严格的时钟同步，而且还须知道发射信号波形的先验知识，因而多用于规则脉冲信号情形。采用互相关技术直接计算信号到达时间差，可以完全摆脱对发射信号波形的限制，可用于对雷电脉冲到达时间差的估计，其应用的前提条件依然是各探测站之间必须保持高精度的时钟同步。

目前常用的时延估计方法有：广义相关时延估计法，如以P.R.Roth命名的ROTH加权广义相关法，以J.C.Hassab和R.E.Boucher命名的HB加权法；用三阶累积量或三阶矩的傅立叶变换的广义双谱时延估计法；利用功率谱与相关函数互为傅立叶变换的广义相位谱时延估计法；以自适应滤波或自适应迭代为特征的自适应时延估计法；基于子空间分解的超分辨时延估计法。在广义相关时延估计法中，在对两个探测站的信号进行相关处理前，先对其进行预滤波，相当于在频域进行加权，有利于增强接收信号的谱分量，提高信噪比，从而获得较高的时延估计精度。这些方法若直接用来求取雷电脉冲的TDOA，其缺点是需要遍历各探测站接收到的所有雷电脉冲信号，运算量非常大，且由于相关峰判决门限不易确定，导致出现许多虚假估计值。

目前，尚未发现对云闪雷电辐射源进行定位时，具备高精度、鲁棒性雷电脉冲信号到达时间差估计的完整方法和实现流程。

发明内容

本发明的目的在于提供一种估计云闪雷电脉冲信号到达时间差的方法，该方法无需遍历两个探测站接收到的所有雷电脉冲，其计算量和虚假估计结果明显减少，估计精度得到提高，且同时也获得了配对的TDOA信息。

本发明提供一种估计云闪雷电脉冲信号到达时间差的方法，该方法包括下述步骤：

第1步对上传到中心站的雷电波形数据块进行拼接和波形还原，得到一个个完整的雷电脉冲，形成雷电脉冲簇；

第2步中心站根据设定的时间区间提取位于该时间区间内的各探测站雷电脉冲簇；

第3步针对用户输入时间区间内、各探测站的雷电脉冲簇，中心站统计分属于每个探测站的雷电脉冲簇中的脉冲个数，由此计算出这些脉冲所对应的原始雷电波形数据块的总个数，并将该总个数最多的那个探测站作为该时间区间内进行到达时间差计算的参考探测站，命名为A站，其它探测站称之为一般探测站，命名为B站、C站、D站、......，设一般探测站的数量为S；

第4步对所述时间区间内参考探测站雷电脉冲簇中的脉冲个数Q进行统计；

第5步令q＝1，s＝1；

第6步取参考探测站雷电脉冲簇中的第q个脉冲，与第s个一般探测站的脉冲簇中的各个脉冲进行粗关联、细关联，得到同源脉冲配对结果，再计算同源脉冲之间的到达时间差值，最后形成完整的同源脉冲配对及其对应的到达时间差队列；

第7步令s＝s+1，判断s＞S是否成立，如果成立，转入第8步，否则进入第6步；

第8步令q＝q+1，判断q＞Q是否成立，如果成立，转入第9步，否则令s＝1，转入第6步；

第9步改变时间区间，重复前述步骤(2)～(8)，得到位于该时间区间内的同源脉冲对队列和到达时间差队列。

本发明方法鲁棒地估计云闪雷电脉冲信号到达任意两个探测站的时间差，它通过若干预处理和判断过程缩小可能关联的雷电脉冲信号范围直至完成雷电脉冲信号的一一配对(得到同源脉冲对)，再利用广义互相关函数或脉冲的包络能量重心来最终估计出TDOA，这样可以大幅度减小TDOA估计的计算量及可能出现的虚假结果，并依据云闪雷电脉冲的持续时间长度选择不同的TDOA具体计算方法。具体而言，本发明具有以下特点：

(1)采用两级关联进行脉冲配对，方法简单，运算量低，遗漏或虚假配对的情况大幅度减小，方法鲁棒性强；

(2)依据云闪雷电脉冲信号的持续时间长短来选择采用不同的TDOA估计方法，兼顾了计算精度、计算量和结果可信度；

(3)通过在互相关函数峰值点附近进行插值，在运算量增加较少的前提下，使TDOA估计的精度突破了整数倍波形采样间隔的限制，实现了连续时差估计；

(4)本发明方法可应用于低信噪比、高信噪比环境，可应用于白噪声、色噪声环境，易于工程实现。

本发明专门针对云闪雷电脉冲信号的特点进行设计，给出了波形关联与同源性判定、到达时间差估计的完整过程，可胜任白噪声和色噪声环境下的TDOA估计要求。本发明所指方法，鲁棒性强，简单可靠，计算结果精度高，易于工程实现。

附图说明

图1是一种估计云闪雷电脉冲到达时间差的计算流程；

图2是雷电发生源与两个探测站的空间几何关系；

图3是雷电脉冲信号的包络能量重心示意图；

图4是当两个同源脉冲的长度相同时，利用包络能量重心计算TDOA的原理示意图(情况一)；

图5是当两个同源脉冲的长度相同时，利用包络能量重心计算TDOA的原理示意图(情况二)；

图6是当两个同源脉冲的长度不相同时，利用包络能量重心计算TDOA的原理示意图。

具体实施方式

以下结合附图对本发明的内容及具体实施方式作进一步说明。

图1给出了一种在由多个探测站和一个中心站组成的云闪雷电定位系统中，形成同源脉冲队列并计算雷电信号到达任意两个探测站的时间差的方法流程，图2给出了在直角坐标系中雷电发生源与两探测站的空间几何关系，图3给出了雷电信号包络能量重心的定义，图4、图5和图6分别给出了依据两个同源脉冲的包络能量重心计算信号到达时间差(TDOA)的原理示意图。

如图1所示，本发明提供一种鲁棒的云闪雷电脉冲到达时间差的估计方法，其步骤包括：

(1)雷电脉冲数据块拼接、波形还原、形成雷电脉冲簇

为了利用各个探测站上传的雷电波形数据块进行TDOA估计，首先需对上传到中心站的雷电波形数据块进行拼接和波形还原，得到一个个完整的雷电脉冲，形成雷电脉冲簇。

以同一个探测站上传到中心站的各个雷电波形数据块为例进行说明。首先，查询该探测站上传到中心站的雷电波形数据块的起始时间、终止时间标记，若前后两个数据块的终止时间、起始时间是连续的，则将这两个数据块拼接为一个较长的数据块。依次类推，完成该探测站上传的所有数据块的起止时间查询和连续数据块的拼接。若前后两个数据块的终止时间、起始时间是不连续的，则停止本个连续数据块的拼接，重新开始下一个新数据块的拼接过程。

对各个探测站上传的所有雷电波形数据块，按上述方式完成各自全部数据块的拼接处理。

对上述拼接得到的每个数据块内的数据，按下式进行波形还原，形成雷电脉冲信号数据块：

式中，n表示脉冲波形采样点的采样时刻，f₀表示接收信号的中心频率，f_s表示采样率，I(n)和Q(n)分别表示雷电信号的同相分量和正交分量，x(n)表示还原出的雷电波形在时刻n时的采样值。

在每个还原出的雷电脉冲信号数据块前部标记以下信息：所属探测站的序号、探测站接收到该脉冲的起止时刻、探测站开始向中心站上传该脉冲数据的时刻、中心站接收完该脉冲数据的时刻。

对应于同一个探测站的全部雷电脉冲，形成该站的雷电脉冲簇。

(2)输入时间区间，提取各站的雷电脉冲簇

中心站以用户输入的时间区间为依据，取出各探测站雷电脉冲簇中脉冲起始时刻位于该时间区间内的雷电脉冲簇数据。

(3)确定参考探测站

针对用户输入时间区间内、各探测站的雷电脉冲簇，中心站统计分属于每个探测站的雷电脉冲簇中的脉冲个数，由此计算出这些脉冲所对应的原始雷电波形数据块的总个数，并将该总个数最多的那个探测站作为该时间区间内进行TDOA计算的参考探测站，将其命名为A站，而其它探测站命名为B站、C站、D站、......，以下计算过程均以B站为例。

注意到一个完整的云闪雷电过程的持续时间从数十微秒至数百毫秒不等，当相邻探测站的间距不超过200公里时，用户输入时间区间的基本原则是：起止时刻差不超过一整秒。

(4)统计A站雷电脉冲簇中的脉冲个数

确定了参考探测站即A站后，对该时间区间内A站雷电脉冲簇中的脉冲个数Q进行统计。

(5)雷电脉冲关联与TDOA计算

依次从A站雷电脉冲簇中取出第一个、第二个、......、第Q个脉冲，与B站、C站、D站、......等各个站的脉冲簇中的各个脉冲进行粗关联、细关联，得到同源脉冲配对结果，再计算同源脉冲之间的TDOA值，最后形成完整的同源脉冲配对及其对应的TDOA队列。以下以A站第q个脉冲与B站的脉冲簇进行同源性配对和TDOA计算过程为例，说明具体实施步骤。

(5.1)雷电脉冲簇的同源性粗关联

以A站的第q个脉冲为基准，将其与B站的脉冲簇进行粗关联的模块。此处所谓的粗关联，是指对两个雷电脉冲的波形相似性进行初步判断。

如图2，假设τ_S时刻在空间S点处发生雷电，其辐射的雷电脉冲分别于τ_A和τ_B时刻抵达A、B两站。则必有：

|(τ_A-τ_S)·V_c-(τ_B-τ_S)·V_c|≤D_AB；

即：|(τ_A-τ_B)·V_c|≤D_AB

或： ${\mathrm{\τ}}_A-\frac{D_{\mathrm{AB}}}{V_c}\≤{\mathrm{\τ}}_B\≤{\mathrm{\τ}}_A+\frac{D_{\mathrm{AB}}}{V_c}$

其中，D_AB为A站和B站之间的直线距离，V_c为雷电波传播速度。

据此可以对探测站A、B接收到的雷电脉冲簇进行同源性粗关联：以A站脉冲簇中的第q个脉冲起始时标τ_A为基点，将B站脉冲簇中脉冲起始时标位于时间区间

内的所有脉冲粗关联为A站该脉冲的潜在同源脉冲，形成B站潜在同源脉冲簇。

为了减少上述粗关联判决引起的遗漏情况、将B站尽可能多的脉冲纳入到粗关联潜在同源脉冲簇中，可将上述粗关联判决公式修正为：

|(τ_A-τ_S)·V_c-(τ_B-τ_S)·V_c|≤D_AB+V_c·max{T_A，T_B}

即： ${\mathrm{\τ}}_A-\frac{D_{\mathrm{AB}}}{V_c}-\max \{T_A,T_B\}\≤{\mathrm{\τ}}_B\≤{\mathrm{\τ}}_A+\frac{D_{\mathrm{AB}}}{V_c}+\max \{T_A,T_B\}$

即：以A站脉冲簇中的第q个脉冲起始时标τ_A为基点，将B站脉冲簇中脉冲起始时标位于时间区间

内的所有脉冲粗关联为A站该脉冲的潜在同源脉冲，形成B站潜在同源脉冲簇。式中，T_A、T_B分别为进行粗关联的A站、B站雷电脉冲的持续时间。

(5.2)雷电脉冲的同源性细关联

A站第q个脉冲与B站的潜在同源脉冲簇进行细关。此处所谓的细关联，是指对两个雷电脉冲的波形相似性进行最终判断。

在完成前述的雷电脉冲粗关联得到潜在同源脉冲簇后，再进行同源性细关联，实现A站与其它探测站雷电脉冲的一一配对(形成同源脉冲对)。以与B站潜在同源脉冲簇的配对为例，说明本发明所指的三种细关联方法：

①从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲之间的包络互相关系数ρ_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预先设定的互相关系数判决门限进行比较或直接按照取大原则，找出与A站该脉冲真正同源的B站脉冲，形成A、B站同源脉冲对。

其中，雷电脉冲信号的包络为：

$y\left(n\right)=\sqrt{I^2\left(n\right)+Q^2\left(n\right)}$

A站与B站雷电脉冲信号的包络互相关系数按下式计算：

${\mathrm{\ρ}}_i=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left\{\right[y_A\left(n\right)-{\mathrm{\μ}}_A]\·[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}\left]\right\}}{\sqrt{\left\{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^2\right\}\·\left\{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\right\}}},$ i＝1，2，3...，I

式中， ${\mathrm{\μ}}_A=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{y}_A\left(n\right),$ ${\mathrm{\μ}}_{\mathrm{Bi}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{y}_{\mathrm{Bi}}\left(n\right),$ I为B站潜在同源脉冲簇队列中的脉冲个数，y_A(n)和y_Bi(n)分别为A、B站的雷电脉冲包络，其采样点数分别为M_A和M_Bi。当M_A≠M_Bi时(即y_A(n)和y_Bi(n)两个脉冲的长度不相等时)，须对较短的脉冲波形补零直至与另一脉冲波形的长度相等。既可在脉冲前端补零，也可在脉冲后端补零。因此，N＝max{M_A，M_Bi}。

为了在不同的信噪比条件下，兼顾细关联的虚警概率和检测概率，可对计算出的全部互相关系数ρ₁，ρ₂，…，ρ_I按大小排序，并按高于判决门限原则或取大原则筛选出同源脉冲。其中，互相关系数判决门限设定为0.5～0.9之间的一个数值。

②从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲的波形包络归一化均方偏差J_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预先设定的均方偏差判决门限进行比较或直接按照取小原则，找出与A站该脉冲真正同源的B站脉冲，形成同源脉冲对。

其中，雷电脉冲的波形包络归一化均方偏差为：

$J_i=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_A\left(n\right)-G_{\mathrm{ABi}}\·y_{\mathrm{Bi}}\left(n\right)]}^2,$ i＝1，2，3...，I

式中，G_ABi为B站脉冲包络相对于A站脉冲包络峰值的归一化系数，其值等于A站脉冲包络的最大值与B站脉冲包络的最大值之比，I为B站潜在同源脉冲簇队列中的脉冲个数。当M_A≠M_Bi时(即y_A(n)和y_Bi(n)两个脉冲的长度不相等时)，须对较短的脉冲波形补零直至与另一脉冲的长度相等，补零的基本原则为：在较短脉冲的前、后两端同时补零。因此，N＝max{M_A，M_Bi}。

为了在不同的信噪比条件下，兼顾细关联的虚警概率和检测概率，可对计算出的全部均方偏差J₁，J₂，…，J_I按大小排序，并按低于判决门限原则或取小原则筛选出同源脉冲。其中，均方偏差判决门限设定为0～1.0之间的一个数值。

③从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲的波形包络偏度偏差S_i或(和)波形包络峰度偏差K_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预先设定的偏度偏差判决门限或(和)峰度偏差判决门限进行比较，找出与A站该脉冲真正同源的B站脉冲，形成同源脉冲对。其中，偏度偏差判决门限和峰度偏差判决门限均可设定为0～1.0之间的一个数值。

其中，雷电脉冲的波形包络偏度偏差和波形包络峰度偏差分别为：

$S_i=\frac{E\left\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^3\right\}}{{\left\{E\right\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^2\left\}\right\}}^{\frac{3}{2}}}-\frac{E\left\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^3\right\}}{{\left\{E\right\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\left\}\right\}}^{\frac{3}{2}}},$ i＝1，2，3...，I

$K_i=\frac{E\left\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^4\right\}}{{\left\{E\right\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^2\left\}\right\}}^2}-\frac{E\left\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^4\right\}}{{\left\{E\right\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\left\}\right\}}^2},$ i＝1，2，3...，I

式中，E(·)为计算均值。I为B站潜在同源脉冲簇队列中的脉冲个数。当M_A≠M_Bi时(即y_A(n)和y_Bi(n)两个脉冲的长度不相等时)，须对较短的脉冲波形补零直至与另一脉冲的长度相等。既可在脉冲前端补零，也可在脉冲后端补零。因此，N＝max{M_A，M_Bi}。

为了在不同的信噪比条件下，兼顾细关联的虚警概率和检测概率，可对计算出的全部偏度偏差和峰度偏差S₁，S₂，…，S_I和K₁，K₂，…，K_I按大小排序，并按取小原则筛选出同源脉冲。

(5.3)对同源脉冲对进行TDOA估计

由于不同时间和空间点发生的雷电脉冲，其持续时间可能相差很大，因此本发明给出了两类计算TDOA的方法。以下以A、B站的同源脉冲对为例具体说明这两类计算方法。

①短脉冲情况下的TDOA计算

当同源脉冲对中的脉冲持续时间(脉冲长度)较短时(例如，N＝200或400)，按下式分别计算同源脉冲对中的A站、B站脉冲的包络能量重心：

${\mathrm{\Δ}}_A=\frac{\underset{n=-N_A}{\overset{n=N_A}{\mathrm{\Σ}}}{\mathrm{ny}}_A^2\left(n\right)}{\underset{n=-N_A}{\overset{n=N_A}{\mathrm{\Σ}}}{y}_A^2\left(n\right)},$ ${\mathrm{\Δ}}_{\mathrm{Bi}}=\frac{\underset{n=-N_{\mathrm{Bi}}}{\overset{n=N_{\mathrm{Bi}}}{\mathrm{\Σ}}}{\mathrm{ny}}_{\mathrm{Bi}}^2\left(n\right)}{\underset{n=-N_{\mathrm{Bi}}}{\overset{n=N_{\mathrm{Bi}}}{\mathrm{\Σ}}}{y}_{\mathrm{Bi}}^2\left(n\right)},$ n≠0

图3说明了本发明所指的脉冲的包络能量重心的含义。将脉冲包络的时间中心点定为参考零点，其前、后段设置为相同采样点数，则能量重心为相对于该参考零点的偏移量。图4、图5和图6表明了根据两个脉冲的起始时刻和重心偏差计算脉冲到达A、B站的时间差的原理。

所以，上述同源脉冲对中的两个脉冲的能量重心实际时刻为：

$t_{\mathrm{AG}}=t_A+(N_A-1)\·T_0+\frac{T_0}{2}+{\mathrm{\Δ}}_A\·T_0,$ $t_{\mathrm{BGi}}=t_{\mathrm{Bi}}+(N_{\mathrm{Bi}}-1)\·T_0+\frac{T_0}{2}+{\mathrm{\Δ}}_{\mathrm{Bi}}\·T_0$

因此，该雷电脉冲到达A站、B站的时间差为：

${\hat{D}}_{A,\mathrm{Bi}}=t_{\mathrm{AG}}-t_{\mathrm{BGi}}=(t_A-t_{\mathrm{Bi}})+(N_A-N_{\mathrm{Bi}})\·T_0+({\mathrm{\Δ}}_A-{\mathrm{\Δ}}_{\mathrm{Bi}})\·T_0$

式中，M_A＝2N_A为同源脉冲对中的A站脉冲长度(采样点数)，M_Bi＝2N_Bi为同源脉冲对中的B站脉冲长度(采样点数)，T₀为采样点间隔(T₀＝1/f_s，f_s为采样频率)，t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻。

②长脉冲情况下的TDOA计算

当同源脉冲对中的脉冲持续时间(脉冲长度)较长时(例如，N＝600或800或1000等情形)，采用脉冲信号的互相关函数峰值点的位置作为雷电信号到达A站和B站的时间差(TDOA)。本发明所指的长脉冲情况下的TDOA估计方法，包括以下三种：

(a)对于同源脉冲对中的A站与B站脉冲，按下式计算它们的互相关函数：

$R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)=E\left\{x_A\left(n\right)x_{\mathrm{Bi}}(n-\mathrm{\τ})\right\}$

由此得到该雷电脉冲到达A站和B站的时间差为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·\arg \left\{\underset{\mathrm{\τ}}{\max }\right|\left|R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)\right|\left|\right\},$ m＝1，2，3...，M

式中，M为该段雷电信号的互相关函数的峰值点个数，对应于多源或多径雷电信号的个数；T₀为采样点间隔(T₀＝1/f_s，f_s为采样频率)；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻；

表示

的绝对值取极大值时所对应的自变量τ之值。

当A站与B站脉冲长度不相等且分别为N_A和N_Bi时，须首先将两个脉冲序列补零至长度均为N_A+N_Bi-1，再计算互相关函数。对脉冲序列补零的原则：在原脉冲序列的后端补零。

(b)对于同源脉冲对中的A站与B站脉冲，按下式分别计算A站和B站脉冲信号的自谱

和互谱

$S_{x_A{x}_A}\left(k\right)=\mathrm{DFT}\left[x_A\left(n\right)\right]$

$S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}\left(k\right)=\mathrm{DFT}\left[x_{\mathrm{Bi}}\left(n\right)\right]$

$S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)=S_{x_A{x}_A}\left(k\right)\·S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}^{*}\left(k\right)$

由此计算最大似然加权函数：

$W\left(k\right)=\frac{1}{\left|S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)\right|}\·\frac{C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)}{[1-C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)]}$

其中， $C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)=\frac{{\left|S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)\right|}^2}{S_{x_A{x}_A}\left(k\right)S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}\left(k\right)}$

进而计算广义互相关函数：

$R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)=\mathrm{IDFT}[W\left(k\right)\·S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)]$

由此得到该雷电脉冲到达A站和B站的时间差为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·\arg \left\{\underset{\mathrm{\τ}}{\max }\right|\left|R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)\right|\left|\right\},$ m＝1，2，3...，M

式中，M为该段雷电信号的互相关函数的峰值点个数，对应于多源或多径雷电脉冲的个数；T₀为采样点间隔(T₀＝1/f_s，f_s为采样频率)；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻；IDFT[.]表示你傅里叶变换；

当A站与B站脉冲长度不相等且分别为N_A和N_Bi时，须首先将两个脉冲序列补零至长度均为N_A+N_Bi-1，再计算互相关函数。对脉冲序列补零的原则：在原脉冲序列的后端补零。

(c)为了进一步提高TDOA估计的精度，在对前述计算出的互相关函数进行峰值搜索时，对得到的峰值点位置进行次极大值点线性插值和相邻三点抛物线插值后再计算新的峰值点位置，并以此作为TDOA估计值输出。

设对互相关函数

进行搜索得到的峰值点为d_m，其互相关函数值为

d_m点的前、后相邻点为d_m-1和d_m+1；d_m点的前、后相邻次极大值点为d_m-k₁和d_m+k₂，其互相关函数值分别为

k₁和k₂为前、后相邻次极大值点到峰值点d_m的距离。

对d_m-k₁点、d_m点、d_m+k₂点进行线性插值，并将插值结果分别赋予d_m-1点和d_m+1点：

${\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m-1)=\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|R_{x_A{x}_{\mathrm{Bi}}}(d_m-k_l)\right|}{k_l}$

${\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m+1)=\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|R_{x_A{x}_{\mathrm{Bi}}}(d_m+k_2)\right|}{k_2}$

然后，对d_m-1点、d_m点、d_m+1点及其对应的互相关函数值

进行三点抛物线插值，最后得到的TDOA值为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·[\frac{2d_m-1}{2}-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|{\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m-1)\right|}{\left|{\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m+1)\right|-2\·\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|+\left|{\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m-1)\right|}],$ m＝1，2，3...，M

式中，M为该段雷电信号中包含的多源或多径雷电信号的个数，T₀为采样点间隔(T₀＝1/f_s，f_s为采样频率)；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻。

(5.4)形成同源脉冲对队列、同源TDOA队列，计算雷电信号到达任意两个探测站的时间差。

按照前述各步骤，得到A站与B站同源脉冲对和TDOA后，将其加入到同源脉冲对队列和TDOA队列中；然后按照相同步骤得到A站与C站的同源脉冲对和TDOA，并将其加入到同源脉冲对队列和TDOA队列中，得到A站与D站的同源脉冲对和TDOA，并将其加入到同源脉冲对队列和TDOA队列中；如此重复，直至最后形成完整的同源脉冲对队列及对应得TDOA队列。

在上述所有站同源脉冲配对队列中，任取两站对应的雷电脉冲，按照前述步骤(5.3)给出的TDOA计算方法，得到雷电信号到达任意两个探测站的时间差。

(6)形成用户指定时间区间内的所有同源脉冲对队列和TDOA队列

重复前述步骤(5)，在用户指定时间区间内，关联出A站的每一个脉冲与B站、C站、D站、......的同源脉冲对及其中任意两个脉冲的到达时间差TDOA，形成用户指定时间区间内的所有同源脉冲对队列和TDOA队列。

(7)改变时间区间，计算并形成新的同源脉冲对队列和TDOA队列

重复前述步骤(2)～(6)，根据用户新输入的时间区间，重新确定参考探测站，然后进行同源脉冲关联和TDOA计算，最后得到该时间区间内的同源脉冲对队列和TDOA队列。

本发明不仅局限于上述具体实施方式，本领域一般技术人员根据本发明公开的内容，可以采用其它多种具体实施方式实施本发明，因此，凡是采用本发明的设计结构和思路，做一些简单的变化或更改的设计，都落入本发明保护的范围。

Claims (8)

1.一种估计云闪雷电脉冲信号到达时间差的方法，该方法包括下述步骤：

第1步对上传到中心站的雷电波形数据块进行拼接和波形还原，得到一个个完整的雷电脉冲，形成雷电脉冲簇；

第2步中心站根据设定的时间区间提取位于该时间区间内的各探测站雷电脉冲簇；

第3步针对用户输入时间区间内、各探测站的雷电脉冲簇，中心站统计分属于每个探测站的雷电脉冲簇中的脉冲个数，由此计算出这些脉冲所对应的原始雷电波形数据块的总个数，并将该总个数最多的那个探测站作为该时间区间内进行到达时间差计算的参考探测站，命名为A站，其它探测站称之为一般探测站，命名为B站、C站、D站、……，设一般探测站的数量为S；

第4步对所述时间区间内参考探测站雷电脉冲簇中的脉冲个数Q进行统计；

第5步令q=1,s=1；

第6步取参考探测站雷电脉冲簇中的第q个脉冲，与第s个一般探测站的脉冲簇中的各个脉冲进行粗关联、细关联，得到同源脉冲配对结果，再计算同源脉冲之间的到达时间差值，最后形成完整的同源脉冲配对及其对应的到达时间差队列；

第7步令s=s+1，判断s＞S是否成立，如果成立，转入第8步，否则进入第6步；

第8步令q=q+1，判断q＞Q是否成立，如果成立，转入第9步，否则令s=1，转入第6步；

第9步改变时间区间，重复前述步骤（2）～（8），得到位于该时间区间内的同源脉冲对队列和到达时间差队列。

2.根据权利要求1所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第1步中，按照下述方式对雷电波形数据块进行拼接：查询每个探测站上传到中心站的雷电波形数据块的起始时间、终止时间标记，若前后两个数据块的终止时间、起始时间是连续的，则将这两个数据块拼接为一个较长的数据块；依次类推，完成该探测站上传的所有数据块的起止时间查询和连续数据块的拼接；若前后两个数据块的终止时间、起始时间是不连续的，则停止对本次连续数据块的拼接，重新开始下一个新数据块的拼接过程。

3.根据权利要求1所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第1步中，按下式进行波形还原，形成雷电脉冲信号数据块：

式中，n表示脉冲波形采样点的采样时刻，f₀表示接收信号的中心频率，f_s表示采样率，I(n)和Q(n)分别表示雷电信号在采样时刻n时的同相分量和正交分量，x(n)表示还原出的雷电波形在采样时刻n时的采样值。

4.根据权利要求1所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第6步中，粗关联是指对两个雷电脉冲的波形相似性进行初步判断，

假设τ_S时刻在空间S点处发生雷电，其辐射的雷电脉冲分别于τ_A和τ_B时刻抵达A、B两站，粗关联判决公式为：

｜(τ_A-τ_S)·V_c-(τ_B-τ_S)·V_c|≤D_AB+V_c·max{T_A,T_B}

即： ${\mathrm{\τ}}_A-\frac{D_{\mathrm{AB}}}{V_c}-\max \{T_A,T_B\}\≤{\mathrm{\τ}}_B\≤{\mathrm{\τ}}_A+\frac{D_{\mathrm{AB}}}{V_c}+\max \{T_A,T_B\}$

其中，D_AB为A站和B站之间的直线距离，V_c为雷电信号传播速度，T_A、T_B分别为进行粗关联的A站、B站雷电脉冲的持续时间；

以A站脉冲簇中的第q个脉冲起始时标τ_A为基点，将B站脉冲簇中脉冲起始时标位于时间区间 $({\mathrm{\τ}}_A-\frac{D_{\mathrm{AB}}}{V_c}-\max \{T_A,T_B\},{\mathrm{\τ}}_A+\frac{D_{\mathrm{AB}}}{V_c}+\max \{T_A,T_B\left\}\right)$ 内的所有脉冲粗关联为A站该脉冲的潜在同源脉冲，形成B站潜在同源脉冲簇。

5.根据权利要求4所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第6步中，细关联是指对两个雷电脉冲的波形相似性进行最终判断，其过程为：

从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲的包络互相关系数ρ_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预先设定的互相关系数判决门限进行比较或直接按照取大原则，找出与A站该脉冲真正同源的B站脉冲，形成A、B站同源脉冲对；

其中，雷电脉冲信号的包络为：

$y\left(n\right)=\sqrt{I^2\left(n\right)+Q^2\left(n\right)}$

I(n)和Q(n)分别表示雷电信号的同相分量和正交分量；

A站与B站雷电脉冲信号的包络互相关系数按下式计算：

${\mathrm{\ρ}}_i=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\{{[y}_A\left(n\right)-{\mathrm{\μ}}_A]\·[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]\}}{\sqrt{\left\{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^2\right\}\·\left\{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\right\}}}i=\mathrm{1,2,3}...,I$

式中， ${\mathrm{\μ}}_A=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{y}_A\left(n\right),$ ${\mathrm{\μ}}_{\mathrm{Bi}}=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{y}_{\mathrm{Bi}}\left(n\right),$ I为B站潜在同源脉冲簇队列中的脉冲个数，y_A(n)和y_Bi(n)为A、B站在采样时刻n时的雷电脉冲包络波形，其采样点数分别为M_A和M_Bi,当M_A≠M_Bi时，须对较短的脉冲波形补零直至与另一脉冲的长度相等,得到N=max{M_A,M_Bi}；对计算出的全部互相关系数ρ₁ ,ρ₂,…,ρ_I按大小排序，并按过门限原则或取大原则筛选出同源脉冲，其中，互相关系数判决门限设定为0.5～0.9之间的一个数值。

6.根据权利要求4所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第6步中，细关联是指对两个雷电脉冲的波形相似性进行最终判断，其过程为：

从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲的波形包络归一化均方偏差J_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预设均方偏差门限比较或取小原则，找出与A站该脉冲真正同源的B站脉冲，形成同源脉冲对；

其中，雷电脉冲的波形包络归一化均方偏差J_i为：

$J_i=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[y_A\left(n\right)-G_{\mathrm{ABi}}\·y_{\mathrm{Bi}}\left(n\right)]}^2,i=\mathrm{1,2,3}...,I$

式中，N为雷电脉冲信号采样点的总个数，y_A(n)和y_Bi(n)分别为A、B站在采样时刻n时的雷电脉冲包络波形，其采样点数分别为M_A和M_Bi，G_ABi为B站脉冲包络相对于A站脉冲包络峰值的归一化系数，其值等于A站脉冲包络的最大值与B站脉冲包络的最大值之比，I为B站潜在同源脉冲簇队列中的脉冲个数；当M_A≠M_Bi时，须对较短的脉冲波形补零直至与另一脉冲的长度相等，得到N=max{M_A,M_Bi}；

对计算出的全部均方偏差J₁,J₂,…,J_I按大小排序，并按低于门限原则或取小原则筛选出同源脉冲，其中，均方偏差判决门限设定为0～1.0之间的一个数值。

7.根据权利要求4所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第6步中，细关联是指对两个雷电脉冲的波形相似性进行最终判断，其过程为：

从B站潜在同源脉冲簇中，按先后顺序逐次取出各个脉冲并计算其与A站该脉冲的波形包络偏度偏差S_i或/和波形包络峰度偏差K_i，该计算遍历B站潜在同源脉冲簇队列中的所有脉冲，以此作为雷电脉冲是否真正同源的判决参数，通过与预设偏度偏差门限或/和峰度偏差门限比较，找出与A站该脉冲真正同源的B站脉冲，形成同源脉冲对，其中，偏度偏差判决门限和峰度偏差判决门限均为0～1.0之间的一个数值；

其中，雷电脉冲的波形包络偏度偏差S_i和波形包络峰度偏差K_i分别为：

$S_i=\frac{E\left\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^3\right\}}{{\left\{E\right\{{[y}_A\left(n\right)-{\mathrm{\μ}}_A{]}^2\left\}\right\}}^{\frac{3}{2}}}-\frac{E\left\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^3\right\}}{{\left\{E\right\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\left\}\right\}}^{\frac{3}{2}}},i=\mathrm{1,2,3}...,I$

$K_i=\frac{E\left\{{[y_A\left(n\right)-{\mathrm{\μ}}_A]}^4\right\}}{{\left\{E\right\{{[y}_A\left(n\right)-{\mathrm{\μ}}_A{]}^2\left\}\right\}}^2}-\frac{E\left\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^4\right\}}{{\left\{E\right\{{[y_{\mathrm{Bi}}\left(n\right)-{\mathrm{\μ}}_{\mathrm{Bi}}]}^2\left\}\right\}}^2},i=\mathrm{1,2,3}...,I$

式中，E(·)为计算均值，y_A(n)和y_Bi(n)为A、B站在采样时刻n时的雷电脉冲包络波形，其采样点数分别为M_A和M_Bi，当M_A≠M_Bi时须对较短的脉冲波形补零直至与另一脉冲的长度相等，得到N=max{M_A,M_Bi}；

I为B站潜在同源脉冲簇队列中的脉冲个数；对计算出的全部偏度偏差和峰度偏差S₁,S₂,…,S_I和K₁,K2,…,K_I按大小排序，并按取小原则筛选出同源脉冲。

8.根据权利要求6、7或8所述的估计云闪雷电脉冲信号到达时间差的方法，其特征在于，第6步中，对同源脉冲对进行到达时间差估计过程为：

当同源脉冲对中的脉冲持续时间较短时，按下式分别计算同源脉冲对中的A站、B站脉冲的包络能量重心Δ_A和Δ_Bi：

${\mathrm{\Δ}}_A=\frac{\underset{n=-N_A}{\overset{n=N_A}{\mathrm{\Σ}}}{\mathrm{ny}}_A^2\left(n\right)}{\underset{n=-N_A}{\overset{n=N_A}{\mathrm{\Σ}}}{y}_A^2\left(n\right)},$ ${\mathrm{\Δ}}_{\mathrm{Bi}}=\frac{\underset{n=-N_{\mathrm{Bi}}}{\overset{n=N_{\mathrm{Bi}}}{\mathrm{\Σ}}}{\mathrm{ny}}_{\mathrm{Bi}}^2\left(n\right)}{\underset{n=-N_{\mathrm{Bi}}}{\overset{n=N_{\mathrm{Bi}}}{\mathrm{\Σ}}}{y}_{\mathrm{Bi}}^2\left(n\right)},n\≠0$

将脉冲包络的时间中心点定为参考零点，其前、后段设置为相同采样点数，则能量重心为相对于该参考零点的偏移量；

上述同源脉冲对中的两个脉冲的能量重心实际时刻t_AG、t_BGi分别为：

$t_{\mathrm{AG}}=t_A+(N_A-1)\·T_0+\frac{T_0}{2}+{\mathrm{\Δ}}_A\·T_0,$ $t_{\mathrm{BGi}}=t_{\mathrm{Bi}}+(N_{\mathrm{Bi}}-1)\·T_0+\frac{T_0}{2}+{\mathrm{\Δ}}_{\mathrm{Bi}}\·T_0$

因此，该雷电脉冲到达A站、B站的时间差为：

${\hat{D}}_{A,\mathrm{Bi}}=t_{\mathrm{AG}}-t_{\mathrm{BGi}}=(t_A-t_{\mathrm{Bi}})+(N_A-N_{\mathrm{Bi}})\·T_0+({\mathrm{\Δ}}_A-{\mathrm{\Δ}}_{\mathrm{Bi}})\·T_0$

式中，M′_A=2N_A为同源脉冲对中A站脉冲的长度，M′_Bi=2N_Bi为同源脉冲对中B站脉冲的长度，T₀为采样点间隔，t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻；

当同源脉冲对中的脉冲持续时间较长时，采用脉冲信号的互相关函数峰值点的位置作为雷电信号到达A站和B站的时间差，采用下述三种方式之一进行长脉冲情况下的到达时间差估计：

第一种方式：对于同源脉冲对中的A站与B站脉冲，按下式计算它们的互相关函数

$R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)=E\left\{x_A\left(n\right)x_{\mathrm{Bi}}(n-\mathrm{\τ})\right\}$

其中，E(·)为取平均值计算，x_A(n)为A站还原出的雷电波形在时刻n时的采样值，x_Bi(n-τ)为B站还原出的雷电波形在时刻n-τ时的采样值；

由此得到该雷电脉冲到达A站和B站的时间差

为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·\arg \left\{\underset{\mathrm{\τ}}{\max }\right|\left|R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)\right|\left|\right\},m=\mathrm{1,2,3}...,M$

式中，M为该段雷电信号的互相关函数的峰值点个数，对应于多源或多径雷电信号的个数；T₀为采样点间隔；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻；表示

的绝对值取极大值时所对应的自变量τ之值；

第二种方式：对于同源脉冲对中的A站与B站脉冲，按下式分别计算A站和B站脉冲信号的自谱

和互谱

$S_{x_A{x}_A}\left(k\right)=\mathrm{DFT}\left[x_A\left(n\right)\right]$

$S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}\left(k\right)=\mathrm{DFT}\left[x_{\mathrm{Bi}}\left(n\right)\right]$

$S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)=S_{x_A{x}_A}\left(k\right)\·S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}^{*}\left(k\right)$

此计算最大似然加权函数：

$W\left(k\right)=\frac{1}{\left|S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)\right|}\·\frac{C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)}{[1-C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)]}$

其中， $C_{x_A{x}_{\mathrm{Bi}}}\left(k\right)=\frac{{\left|S_{x_A{x}_{\mathrm{Bi}}}\right|}^2}{S_{x_A{x}_A}\left(k\right)S_{x_{\mathrm{Bi}}{x}_{\mathrm{Bi}}}\left(k\right)}$

进而计算广义互相关函数：

$R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)=\mathrm{IDFT}[W\left(k\right)\·S_{x_A{x}_{\mathrm{Bi}}}\left(k\right)]$

由此得到该雷电脉冲到达A站和B站的时间差为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·\arg \left\{\underset{\mathrm{\τ}}{\max }\right|\left|R_{x_A{x}_{\mathrm{Bi}}}\left(\mathrm{\τ}\right)\right|\left|\right\},m=\mathrm{1,2,3}...,M$

式中，M为该段雷电信号的互相关函数的峰值点个数，对应于多源或多径雷电脉冲的个数；T₀为采样点间隔；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻；IDFT[.]表示逆傅里叶变换；第三种方式：在对前述计算出的互相关函数

进行峰值搜索时，对得到的峰值点位置进行次极大值点线性插值和相邻三点抛物线插值后再计算新的峰值点位置，并以此作为到达时间差估计值输出；设对互相关函数

进行搜索得到的峰值点为d_m，其互相关函数值为

d_m点的前、后相邻点为d_m-1和d_m+1；d_m点的前、后相邻次极大值点为dx_-k₁和d_m+k₂，其互相关函数值分别为

k₁和k₂为前、后相邻次极大值点到峰值点d_m的距离；对d_m-k₁点、d_m点、d_m+k₂点进行线性插值，并将插值结果分别赋予d_m-1点和d_m+1点：

${\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m-1)=\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|R_{x_A{x}_{\mathrm{Bi}}}(d_m-k_1)\right|}{k_1}$

${\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m+1)=\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|R_{x_A{x}_{\mathrm{Bi}}}(d_m+k_2)\right|}{k_2}$

然后，对d_m-1点、d_m点、d_m+1点及其对应的互相关函数值

进行三点抛物线插值，最后得到的到达时间差值为：

${\hat{D}}_m=(t_A-t_{\mathrm{Bi}})+T_0\·[\frac{{2d}_m-1}{2}-\frac{\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|-\left|\hat{R}{x}_A{x}_{\mathrm{Bi}}(d_m-1)\right|}{\left|{\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m+1)\right|-2\·\left|R_{x_A{x}_{\mathrm{Bi}}}\left(d_m\right)\right|+\left|{\hat{R}}_{x_A{x}_{\mathrm{Bi}}}(d_m-1)\right|}],m=\mathrm{1,2,3},...,M$

式中，M为该段雷电信号中包含的多源或多径雷电信号的个数，T₀为采样点间隔；t_A和t_Bi分别表示同源脉冲对中的A站和B站脉冲的起始时刻。

Images