CN116008671A - 一种基于时差和聚类的闪电定位方法 - Google Patents

Abstract

本发明公开了一种基于时差和聚类的闪电定位方法，包括：基于到达时间差原理对闪电回击数据进行四站组合定位，得到初始闪电定位数据，构成初始闪电定位数据集，采用k均值聚类算法对初始闪电定位数据集进行聚类分析得到k个聚类簇集合，选取输出k个聚类簇集合中包含闪电定位数据点最多的簇中心，作为最终的闪电定位结果。剔除定位离群点后得到最终精确定位值。本发明可以有效提高闪电定位精度并减小数据粗差、站网布局等因素对闪电定位效果的影响。

Description

一种基于时差和聚类的闪电定位方法

技术领域

本发明属于气象闪电定位技术领域，具体涉及一种基于时差和聚类的闪电定位方法。

背景技术

闪电是大气中的放电现象，常伴随着冰雹、暴雨等强对流天气过程。闪电具有强大电流、电磁辐射以及炙热高温等物理效应因此常常引起雷击灾害，包括对建筑物、电力设备、信息通信设备、油罐储运等造成的巨大破坏。闪电的监测预警及准确定位对减少雷击灾害事故具有重要的意义。

闪电定位***利用闪电回击辐射的电磁场特性来遥测其发生的时间、位置、强度和极性，已应用于气象、电力、航天航空等领域。定位精度是评价闪电定位***的关键技术指标。到达时间差(Time Difference Of Arrival，TDOA)定位方法由于定位精度较高，已成为主流闪电定位方法。该方法是基于各个闪电探测站与闪电辐射源之间的距离差确定定位双曲线，通过求解双曲线方程组来确定闪电辐射源信号相对于各个探测站相对位置。

TDOA定位方法依据时间差定位，由于雷电电磁场在传播过程中会受到地形、地球电导率等因素干扰，用于定位的原始闪电数据受到各类误差因素的影响会含有粗差，使TDOA算法定位曲线往往无法相交于一点，因此该方法存在定位不准确、抗干扰性差等问题。

发明内容

本发明的目的在于针对上述问题，提出一种基于时差和聚类的闪电定位方法，提高闪电定位精度，进一步地减小数据粗差、站网布局等因素对闪电定位效果的影响。

聚类算法可以从大量模糊、含噪声或随机的实际数据中提取出需要的信息，采用TDOA四站定位方法对闪电回击数据进行组合定位得到初步闪电定位结果，之后利用k均值聚类算法对得到的初步定位结果进行分类并剔除离群点得到最终的闪电定位结果，提高定位的准确率以及抗干扰能力。

为实现上述目的，本发明采取以下技术方案：

第一方面，提供一种基于时差和聚类的闪电定位方法，包括：

步骤S1、获取待定位闪电数据信息，所述闪电数据信息包括至少四个闪电探测站接收到所述待定位闪电信号的闪电回击数据；闪电回击数据包括闪电探测站的站点位置信息以及闪电到达时间；

步骤S2、对所述闪电数据信息进行组合，得到多组闪电数据，其中每组闪电数据包括四个闪电探测站接收到待定位闪电信号的闪电回击数据；

步骤S3、对每组闪电数据，采用四站时差定位法分别对闪电辐射源进行定位，得到初始闪电定位数据；所有的初始闪电定位数据构成闪电定位数据集；

步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，得到k个聚类簇集合；

步骤S5、选取输出k个聚类簇集合中包含闪电定位数据点最多的簇中心n_p,1≤p≤k，作为最终的闪电定位结果。

在一些实施例中，步骤S1中，所述站点位置信息为站点经纬高信息，还包括将站点经纬高信息转换成空间直角坐标系下定位的坐标x、y、z；

其中，闪电探测站的经度L、纬度B、高度H，e为椭球的第一偏心率，N为椭球的卯酉圈曲率半径。

在一些实施例中，步骤S2包括：遍历所有闪电探测站，每四站的闪电数据信息进行组合，得到

组闪电数据，其中n为闪电探测站的总个数。

在一些实施例中，步骤S3中，采用四站时差定位法分别对闪电辐射源进行定位，得到闪电初始定位数据，包括：

根据到达时间差TDOA定位原理，距离差方程为：

其中，闪电的发生位置为(x,y,z)，发生时间为t，第i个闪电探测站的坐标位置为(x_i,y_i,z_i)，到达时间为t_i；i＝0时表示主站，i＝1,2,3…n时表示副站；闪电到主站(x₀,y₀,z₀)的距离为r₀，到第i个副站的距离为r_i，Δr_i为r_i与r₀的距离差，c为电磁波信号的传播速率；

探测站个数为4，即i＝0,1,2,3；式(1)变化为：

ri2-r02＝di-d0 (2)

其中d_i＝(x-x_i)²+(y-y_i)²+(z-z_i)²,d₀＝(x-x₀)²+(y-y₀)²+(z-z₀)²；

对式(2)进行移项、平方、整理化简得：

式(3)中i＝1,2,3，是一个关于(x,y,z,t)的非线性方程组，将r₀作为已知量，得到矩阵表达式：

AX＝B (4)

即：

当闪电探测站不部署在同一平面上时，系数矩阵A的秩等于3，得到：

X ＝ (A^TA)^-1A^TB (6)

将X带入式(1)可以得到方程：

ar₀ ²+br₀+c＝0 (7)

解一元二次方程(7)得到r₀的数值解后，将r₀反代入式(6)中求出闪电辐射源的空间直角坐标系下定位的坐标。

当Δ＝b²-4ac＞0时，r₀有两个解r₀₁、r₀₂，即双曲面有两个交点；r₀表示距离必须为正数，r₀₁、r₀₂一正一负时选取正根作为定位解；r₀₁、r₀₂皆为正数时通过增加方位角辅助信息消除定位模糊；

当Δ＝b²-4ac＝0时r₀有唯一解，不存在定位模糊问题。

在一些实施例中，步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，包括：

S41、第一次，在闪电定位数据集X中随机抽取k个对象构成第一个训练子集T₁(T₂∈X)，其中k＜n，利用K均值算法构建一个局部的包含k个簇的模型，得到k个初步的聚类中心；

S42、第二次，在第一次的基础上，从X中随机抽取C₂个除了T₁以外的对象构成第2个训练子集T₂(T₂∈X-T₁)；利用K均值算法将T₂中的对象加入到k个簇中，并更新每个簇的聚类中心；以此类推，重复迭代步骤S41至步骤S42，直至达到预设迭代停止条件，得到最终的k个类簇的集合S＝{S₁，S₂，S₃，…，S_k}。

进一步地，在一些实施例中，步骤S4包括：

在闪电定位数据集X中随机选取k个闪电定位数据样本点作为初始的聚类中心，记为n⁽⁰⁾＝(n₁ ⁽⁰⁾，…，n_l ⁽⁰⁾，…，n_k ⁽⁰⁾)；其中k＜n，n为闪电探测站的总个数；其中m维向量闪电定位数据集X中，样本点x_i∈X，x_i＝(x_1i，x_2i，…，x_mi)^T；

对固定的聚类中心n^(t)＝(n₁ ^(t)，…，n_l ^(t)，…，n_k ^(t))，其中n_l ^(t)为类簇S_l的聚类中心，按照样本点与聚类中心的距离对闪电定位数据样本点进行聚类：

计算每个样本点到类簇中心的距离，根据计算的距离使每个闪电定位数据样本点归属到与其距离最小的类簇中，得到k个类簇的集合S＝{S₁，S₂，S₃，…，S_k}，生成初步的聚类结果C^(t)；

对聚类结果C^(t)计算当前每个类簇的样本均值u_i：

其中x_i为类别l中的闪电定位数据样本点，s⁽ⁱ⁾为x_i所属类别，z为各个类簇中样本点的总数；将均值作为新的聚类中心u^(t+1)＝(u₁ ^(t+1),…,u_l ^(t+1),…,u_k ^(t+1))；

定义闪电定位数据样本点与其所属类中心之间的距离总和为最终损失函数W(C)：

其中

为第l个类簇的聚类中心；

中I(C(i)＝l)表示取值为0或1的指示函数；函数W(C)表示相同类簇中样本点的相似程度；k均值聚类转换为一个优化问题的求解：

如果迭代收敛或者满足迭代停止条件，即损失函数W(C)达到最小则输出最后聚类结果C^*＝C^(t)，否则继续迭代，令迭代次数t＝t+1且返回重新计算损失函数。

在一些实施例中，步骤S5还包括：将闪电定位结果从空间直角坐标系下定位的坐标x、y、z转换成空间大地坐标系下的经纬高信息；WGS-84椭球模型中经度L、纬度B、高度H的求解公式为：

式中，a、b分别为椭球的长短半轴，a＝6378.137km，b＝6356.752km；e为椭球的第一偏心率，N为椭球的卯酉圈曲率半径。

第二方面，本发明提供了一种基于时差和聚类的闪电定位装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据第一方面所述方法的步骤。

第三方面，本发明提供了一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现第一方面所述方法的步骤。

本发明采用TDOA四站定位方法对闪电回击数据进行组合定位得到闪电定位结果，之后利用k均值聚类算法对得到的定位结果进行分类并剔除离群点得到最终的定位结果，解决传统方法存在定位不准确、抗干扰性差等问题。最后通过测试验证了该方法具有较好的定位效果。

与现有技术相比，本发明的有益效果是：

本发明提供的一种基于时差和聚类的闪电定位方法。TDOA定位方法依据时间差定位，由于雷电电磁场在传播过程中会受到地形、地球电导率等因素干扰，用于定位的原始闪电数据受到各类误差因素的影响会含有粗差，使TDOA算法定位曲线往往无法相交于一点，因此该方法存在定位不准确、抗干扰性差等问题。而聚类算法可以从大量模糊、含噪声或随机的实际数据中提取出需要的信息，基于此提出一种基于聚类算法的闪电定位方法。采用TDOA四站定位方法对闪电回击数据进行组合定位得到闪电定位结果，之后利用k均值聚类算法对得到的定位结果进行分类并剔除离群点得到最终的定位结果，提高定位的准确率以及抗干扰能力。

附图说明

图1是本发明实施例的TDOA定位原理图。

图2是本发明实施例的闪电定位***布站示意图。

图3是本发明实施例的闪电定位方法流程图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

在本发明的描述中，若干的含义是一个以上，多个的含义是两个以上，大于、小于、超过等理解为不包括本数，以上、以下、以内等理解为包括本数。如果有描述到第一、第二只是用于区分技术特征为目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量或者隐含指明所指示的技术特征的先后关系。

本发明的描述中，参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

实施例1

一种基于时差和聚类的闪电定位方法，包括：

步骤S1、获取待定位闪电数据信息，所述闪电数据信息包括至少四个闪电探测站接收到所述待定位闪电信号的闪电回击数据；闪电回击数据包括闪电探测站的站点位置信息以及闪电到达时间；

步骤S2、对所述闪电数据信息进行组合，得到多组闪电数据，其中每组闪电数据包括四个闪电探测站接收到待定位闪电信号的闪电回击数据；

步骤S3、对每组闪电数据，采用四站时差定位法分别对闪电辐射源进行定位，得到初始闪电定位数据；所有的初始闪电定位数据构成闪电定位数据集；

步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，得到k个聚类簇集合；

步骤S5、选取输出k个聚类簇集合中包含闪电定位数据点最多的簇中心n_p，1≤p≤k，作为最终的闪电定位结果。

步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，包括：

S41、第一次，在闪电定位数据集X中随机抽取k个对象构成第一个训练子集T₁(T₂∈X)，其中k＜n，利用K均值算法构建一个局部的包含k个簇的模型，得到k个初步的聚类中心；

S42、第二次，在第一次的基础上，从X中随机抽取C₂个除了T₁以外的对象构成第2个训练子集T₂(T₂∈X-T₁)；利用K均值算法将T₂中的对象加入到k个簇中，并更新每个簇的聚类中心；以此类推，重复迭代步骤S41至步骤S42，直至达到预设迭代停止条件，得到最终的k个类簇的集合S＝{S₁,S₂,S₃,…,S_k}。

在一些实施例中，如图1所示为TDOA定位原理图，TDOA是一种无线定位方法，需要三个或三个以上已知位置坐标的闪电探测站进行定位。图1中共有三个闪电探测站A、B、C，S为闪电发生位置。假设闪电辐射源S到探测站A、B、C的距离分别为r₁、r₂和r₃，实际情况中该距离未知。闪电探测站A和B可测量闪电回击辐射源S发出的电磁波信号到达各自探测站的时间，两站存在时间差，将时间差乘以闪电信号的传播速度可确定两站之间的距离差r₂₁＝r₂-r₁，构成一条以探测站A、B为焦点，距离差r₂₁为长轴的双曲线L₂，闪电发生位置在这条双曲线上的某一点。探测站C与闪电探测站B同样可构成另一条定位双曲线L₁，两条双曲线的交点S即为闪电辐射源位置。

如图2所示，为闪电定位***的布站示意图。假设闪电的发生位置为S(x,y,z)，发生时间为t，第i个闪电探测站的坐标位置为(x_i,y_i,z_i)，到达时间为t_i。i＝0时表示主站，i＝1,2,3…n时表示副站。

在一些具体实施例中，如图3所示，图3为闪电定位流程图，包括：

1)闪电探测站接收到同一条闪电信号的闪电回击数据X；

2)将所有闪电探测站每四站进行组合，最多可以得到

种组合；

3)利用四站时差定位法分别对闪电辐射源进行定位，可以得到闪电定位数据集X＝{p₁,p₂,p₃,…,p_m}；

4)采用k均值聚类算法对所有定位数据进行聚类分析，不含噪声或噪声较小的数据点因聚合性较好可以聚成一类，含噪声较大的数据点聚合成另外一类或几类，共k类：s₁、s₂、s₃…、s_k；

5)选取输出k个聚类簇集合{s₁,s₂,s₃,…,s_k}中包含闪电定位数据点最多的簇中心n_p,1≤p≤k，作为最终的闪电定位结果，排除定位结果中的干扰定位点；

6)根据实际定位结果分析算法的定位性能。

根据到达时间差TDOA定位原理，距离差方程为：

其中，闪电的发生位置为(x,y,z)，发生时间为t，第i个闪电探测站的坐标位置为(x_i,y_i,z_i)，到达时间为t_i；i＝0时表示主站，i＝1,2,3…n时表示副站；闪电到主站(x₀,y₀,z₀)的距离为r₀，到第i个副站的距离为r_i，Δr_i为r_i与r₀的距离差，c为电磁波信号的传播速率；

若探测站个数为4，即i＝0,1,2,3；式(1)变化为：

ri2-r02＝di-d0 (2)

其中d_i＝(x-x_i)²+(y-y_i)²+(z-z_i)²,d₀＝(x-x₀)²+(y-y₀)²+(z-z₀)²；

对式(2)进行移项、平方、整理化简得：

式(3)中i＝1,2,3，是一个关于(x,y,z,t)的非线性方程组，将r₀作为已知量，得到矩阵表达式：

AX＝B (4)

当闪电探测站不部署在同一平面上时，系数矩阵A的秩等于3，得到：

X ＝ (A^TA)^-1A^TB (6)

将X带入式(1)可以得到方程：

ar₀ ²+br₀+c＝0 (7)

解一元二次方程(7)得到r₀的数值解后，将r₀反代入式(6)中求出闪电辐射源的空间直角坐标系下定位的坐标。

进一步地，还包括：将求出的闪电辐射源的空间直角坐标系下定位的坐标x、y、z转换成空间大地坐标系下的经纬高信息。

闪电定位***记录的闪电回击数据是每个闪电探测站的经度、纬度、高度信息(空间大地坐标系)，需要将空间大地坐标系与空间直角坐标系进行转换。

由于地球表面真实形状不是完美的规则球形，使用WGS-84椭球模型作为闪电定位参考模型。WGS-84椭球模型中经度L、纬度B、高度H的求解公式为：

式中，a、b分别为椭球的长短半轴，a＝6378.137km，b＝6356.752km。e为椭球的第一偏心率，N为椭球的卯酉圈曲率半径。

相同基准下空间直角坐标系下定位的坐标x、y、z的求解公式为：

综上，闪电定位***记录接收到闪电回击数据的站点经纬高信息以及到达时间后，经式(9)进行坐标转换，通过方程(7)计算得到距离r₀的解析解，将r₀带入式(6)可求出闪电辐射源的空间直角坐标系下定位的坐标，最终经过式(8)得到闪电的经纬高位置。

在一些实施例中，方程(7)可能会出现解的模糊：

当Δ＝b²-4ac＞0时，r₀有两个解r₀₁、r₀₂，即双曲面有两个交点；r₀表示距离必须为正数，r₀₁、r₀₂一正一负时选取正根作为定位解；r₀₁、r₀₂皆为正数时通过增加方位角辅助信息消除定位模糊；

Δ＝b²-4ac＝0时r₀有唯一解，不存在定位模糊问题；

Δ＝b²-4ac＜0时r₀无解，无法定位闪电辐射源。

在一些实施例中，步骤4)包括如下步骤:

首先初始化闪电定位样本质心。m维向量闪电定位数据集X中，样本元素x_i∈X，x_i＝(x_1i,x_2i,…,x_mi)^T。初次迭代时在集合X中随机选取k(k＜n)个闪电定位数据样本点作为初始化的聚类中心点，记为n⁽⁰⁾＝(n₁ ⁽⁰⁾,…,n_l ⁽⁰⁾,…,n_k ⁽⁰⁾)。

对固定的聚类中心n^(t)＝(n₁ ^(t),…,n_l ^(t),…,n_k ^(t))，其中n_l ^(t)为类簇S_l的中心点，按照样本点与初始聚类中心的距离对闪电定位数据样本点进行聚类。对于样本点x_i与聚类中心点n_j之间的闵式距离为:

采用欧式距离计算闪电数据集X中的每个闪电定位数据样本点到簇聚类中心的距离。式(10)中p＝2可表示为欧式距离:

根据式(11)计算每个样本点到类簇中心的距离，依据计算的距离使每个闪电定位数据样本点归属到与其距离最小的类簇中，可得到k个类簇的集合S＝{S₁,S₂,S₃,…,S_k}，生成初步的聚类结果C^(t)。

对聚类结果C^(t)计算当前每个类簇的样本均值：

其中x_i为类别l中的闪电数据样本点，s⁽ⁱ⁾为x_i所属类别，z为各个类簇中样本点的总数。将均值作为新的聚类中心点u^(t+1)＝(u₁ ^(t+1),…,u_l ^(t+1),…,u_k ^(t+1))。定义闪电样本点与其所属类中心之间的距离总和为最终损失函数：

其中

为第l个类的质心(即聚类中心点)。

中I(C(i)＝l)表示取值为0或1的指示函数。函数W(C)表示相同类簇中样本点的相似程度。k均值聚类转换为一个优化问题的求解：

如果迭代收敛或者满足迭代停止条件，即式(13)中的误差平方和W(C)达到最小则输出最后聚类结果C^*＝C^(t)，否则继续迭代，令迭代次数t＝t+1且返回式(14)重新计算。

本发明提供了一种基于时差和聚类的闪电定位方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

实施例2

第二方面，本实施例提供了一种基于时差和聚类的闪电定位装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据实施例1所述方法的步骤。

实施例3

第三方面，本实施例提供了一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现实施例1所述方法的步骤。

本领域内的技术人员应明白，本申请的实施例可提供为方法、***、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种基于时差和聚类的闪电定位方法，其特征在于，包括：

步骤S1、获取待定位闪电数据信息，所述闪电数据信息包括至少四个闪电探测站接收到所述待定位闪电信号的闪电回击数据；闪电回击数据包括闪电探测站的站点位置信息以及闪电到达时间；

步骤S2、对所述闪电数据信息进行组合，得到多组闪电数据，其中每组闪电数据包括四个闪电探测站接收到待定位闪电信号的闪电回击数据；

步骤S3、对每组闪电数据，采用四站时差定位法分别对闪电辐射源进行定位，得到初始闪电定位数据；所有的初始闪电定位数据构成闪电定位数据集；

步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，得到k个聚类簇集合；

步骤S5、选取输出k个聚类簇集合中包含闪电定位数据点最多的簇中心，作为最终的闪电定位结果。

2.根据权利要求1所述的方法，其特征在于，步骤S1中，所述站点位置信息为站点经纬高信息，还包括将站点经纬高信息转换成空间直角坐标系下定位的坐标x、y、z；

其中，闪电探测站的经度L、纬度B、高度H，e为椭球的第一偏心率，N为椭球的卯酉圈曲率半径。

3.根据权利要求1所述的方法，其特征在于，步骤S2包括：遍历所有闪电探测站，每四站的闪电数据信息进行组合，得到

组闪电数据，其中n为闪电探测站的总个数。

4.根据权利要求1所述的方法，其特征在于，步骤S3中，采用四站时差定位法分别对闪电辐射源进行定位，得到闪电初始定位数据，包括：

根据到达时间差TDOA定位原理，距离差方程为：

其中，闪电的发生位置为(x,y,z)，发生时间为t，第i个闪电探测站的坐标位置为(x_i,y_i,z_i)，到达时间为t_i；i＝0时表示主站，i＝1,2,3…n时表示副站；闪电到主站(x₀,y₀,z₀)的距离为r₀，到第i个副站的距离为r_i，Δr_i为r_i与r₀的距离差，c为电磁波信号的传播速率；

探测站个数为4，即i＝0,1,2,3；式(1)变化为：

r_i ²-r₀ ²＝d_i-d₀ (2)

其中d_i＝(x-x_i)²+(y-y_i)²+(z-z_i)²,d₀＝(x-x₀)²+(y-y₀)²+(z-z₀)²；

对式(2)进行移项、平方、整理化简得：

式(3)中i＝1,2,3，是一个关于(x,y,z,t)的非线性方程组，将r₀作为已知量，得到矩阵表达式：

AX＝B (4)

即：

当闪电探测站不部署在同一平面上时，系数矩阵A的秩等于3，得到：

X＝(A^TA)^-1A^TB (6)

将X带入式(1)可以得到方程：

ar₀ ²+br₀+c＝0 (7)

解一元二次方程(7)得到r₀的数值解后，将r₀反代入式(6)中求出闪电辐射源的空间直角坐标系下定位的坐标。

5.根据权利要求4所述的方法，其特征在于，

当Δ＝b²-4ac＞0时，r₀有两个解r₀₁、r₀₂，即双曲面有两个交点；r₀表示距离必须为正数，r₀₁、r₀₂一正一负时选取正根作为定位解；r₀₁、r₀₂皆为正数时通过增加方位角辅助信息消除定位模糊；

当Δ＝b²-4ac＝0时r₀有唯一解，不存在定位模糊问题。

6.根据权利要求1所述的方法，其特征在于，步骤S4、采用k均值聚类算法对闪电定位数据集中的所有闪电定位数据进行聚类分析，包括：

S41、第一次，在闪电定位数据集X中随机抽取k个对象构成第一个训练子集T₁(T₁∈X)，其中k＜n，利用K均值算法构建一个局部的包含k个簇的模型，得到k个初步的聚类中心；

S42、第二次，在第一次的基础上，从X中随机抽取C₂个除了T₁以外的对象构成第2个训练子集T₂(T₂∈X-T₁)；利用K均值算法将T₂中的对象加入到k个簇中，并更新每个簇的聚类中心；以此类推，重复迭代步骤S41至步骤S42，直至达到预设迭代停止条件，得到最终的k个类簇的集合S＝{S₁,S₂,S₃,…,S_k}。

7.根据权利要求6所述的方法，其特征在于，步骤S4包括：

在闪电定位数据集X中随机选取k个闪电定位数据样本点作为初始的聚类中心，记为n⁽⁰⁾＝(n₁ ⁽⁰⁾,…,n_l ⁽⁰⁾,…,n_k ⁽⁰⁾)；其中k＜n，n为闪电探测站的总个数；其中m维向量闪电定位数据集X中，样本点x_i∈X，x_i＝(x_1i,x_2i,…,x_mi)^T；

对固定的聚类中心n^(t)＝(n₁ ^(t),…,n_l ^(t),…,n_k ^(t))，其中n_l ^(t)为类簇S_l的聚类中心，按照样本点与聚类中心的距离对闪电定位数据样本点进行聚类：

计算每个样本点到类簇中心的距离，根据计算的距离使每个闪电定位数据样本点归属到与其距离最小的类簇中，得到k个类簇的集合S＝{S₁,S₂,S₃,…,S_k}，生成初步的聚类结果C^(t)；

对聚类结果C^(t)计算当前每个类簇的样本均值u_l：

其中x_i为类别l中的闪电定位数据样本点，s⁽ⁱ⁾为x_i所属类别，z为各个类簇中样本点的总数；将均值作为新的聚类中心u^(t+1)＝(u₁ ^(t+1),…,u_l ^(t+1),…,u_k ^(t+1))；

定义闪电定位数据样本点与其所属类中心之间的距离总和为最终损失函数W(C)：

其中

为第l个类簇的聚类中心；

中I(C(i)＝l)表示取值为0或1的指示函数；函数W(C)表示相同类簇中样本点的相似程度；k均值聚类转换为一个优化问题的求解：

如果迭代收敛或者满足迭代停止条件，即损失函数W(C)达到最小则输出最后聚类结果C^*＝C^(t)，否则继续迭代，令迭代次数t＝t+1且返回重新计算损失函数。

8.根据权利要求1所述的方法，其特征在于，步骤S5还包括：将闪电定位结果从空间直角坐标系下定位的坐标x、y、z转换成空间大地坐标系下的经纬高信息；WGS-84椭球模型中经度L、纬度B、高度H的求解公式为：

式中，a、b分别为椭球的长短半轴，a＝6378.137km，b＝6356.752km；e为椭球的第一偏心率，N为椭球的卯酉圈曲率半径。

9.一种基于时差和聚类的闪电定位装置，其特征在于，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据权利要求1至8任一项所述方法的步骤。

10.一种存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至8任一项所述方法的步骤。

Images