引言

雷电定位系统能够实时监测回击发生的时间、位置等信息，便于气象、电力、航天、森林防火等领域进行灾害预警及灾后处理^[1-2]。目前，大多定位系统通过磁向(direction finding，DF)法和到达时间差(time difference of arrival，TDOA)法两种定位算法进行定位^[3]。常用定位算法包括迭代法和非迭代法两类^[4]，由于雷电波沿地球椭球面传播，使用迭代法计算更精确。传统的迭代算法包括最小二乘法^[5]和泰勒级数展开法^[6]，迭代算法需要一个与最终解接近的初始值，且迭代易于发散。为了克服这一缺点，提出网格搜索方法^[7]去寻找最优定位，通过不断缩小雷击点与探测点之间的误差最终获得雷击点目标，在雷电定位中引入粒子群优化(particle swarm optimization，PSO)^[8]和改进的果蝇优化(modified fruit fly optimization algorithm，MFOA)^[9]等算法后，提高了定位精度。该类优化算法易于陷入局部最优且计算量较大，实际应用受限于计算机运算速度，抗误差干扰能力差。刘达新等^[10]提出一种能抵抗粗差的定位算法。本文提出一种新的雷电定位算法DG-LLA(DBSCAN and grid-search lightning location algorithm)，在雷电定位计算中引入DBSC- AN(density-based spatial clustering of applications with noise)方法^[11]与网格搜索方法，对时差法的粗定位结果^[12-13]进行聚类分析^[14]，并结合定位结果的数值特征改进DBSCAN方法，TDOA法对时间误差异常敏感，而原始数据的测量误差无法避免，定位精度低，引入DBSCAN方法后，能有效识别出误差数据，精确稳定地计算出雷击点发生的位置，增强了抗误差干扰的能力；网格搜索方法采用多级搜索策略，提高了局部搜索和全局搜索的能力，克服了传统迭代算法易于发散和陷入局部最优等缺点，另外，在初始值接近目标位置的条件下，优化算法收敛快，精度高。

1 雷电定位原理

TDOA法是一种利用时间差进行定位的方法。通过测量闪电产生的电磁波信号到达不同探测站的时间差，得到以探测站为焦点，距离差为长轴的时差双曲线，两条时差双曲线的交点即为闪电位置。假设闪电位置为(x，y，z)，第i个探测站的位置和闪电到达时间表示为(x_i，y_i，z_i，t_i)，通过探测站与各个探测站的距离建立时差方程组，求解得到闪电位置^[15]。

2 DG-LLA算法设计

TDOA法通过建立时差方程组求解得到雷电的时间和位置信息，计算过程较为简单，能得到基于测量数据的解析解，但该方法只能利用3站或4站进行定位，不能很好利用冗余信息，增大了测量值选取不当带来的误差。本文提出一种新的定位算法DG-LLA，算法思想是当测量数据较少时，通过TDOA法定位得到初定位结果，当测量数据较多时，利用TDOA法对所有测量数据进行3站和4站定位，对得到的若干个定位结果进行DBSCAN聚类分析，将聚合性较好的类中心作为初定位结果，最后通过网格搜索法对初定位结果进行优化，算法流程如图 1所示。

2.1 DBSCAN聚类算法

传统DBSCAN方法需要人工选择邻域参数e和m，e为邻域半径大小，m为邻域内最小样本量。本文需对雷电定位结果进行聚类分析，而实际定位结果间聚合程度各不相同，邻域参数难以唯一确定，很难实现聚类的最优化。本文提出一种改进自适应聚类算法，根据定位结果数据集的统计特性判断邻域参数的取值，实现参数的自适应判别。

2.1.1 e的自适应判别

为了实现领域半径e的自适应判别，首先需要计算数据集的距离分布矩阵，对每一个样本与数据集中其他样本的距离进行升序排列，得到数据集的k近邻分布^[16]。从数据集的k近邻分布可以看到，随着k的增加，排序后的距离曲线大致沿纵轴向上平移，且任何一条曲线都是在平稳变化后急剧上升，因此, 可以判断数据集值大部分值分布较为密集, 将分布曲线峰值点所对应的最近邻距离值作为e。试验发现通过逆高斯分布统计模型进行拟合求解分布曲线峰值效果最好, 文献[16]给出了逆高斯分布的概率密度公式，对逆高斯分布的概率密度公式求一阶导数，在导数为0时获得极值，得到邻域半径e的计算公式。

2.1.2 m的自适应判别

确定邻域半径后，根据核密度估计理论，以对象x为核心，e为半径的空间内存在对象数量即为m值，本文计算数据集中每个点邻域内对象数量，最后求数学期望得到m。文献[11]定义了核函数，通过计算核函数得到以对象x为中心的邻域内的数据量，遍历数据集中所有数据并求数学期望得到对应的m。

2.1.3 基于自适应DBSCAN方法的地闪定位

假设共接收到同一次地闪的数据n条，理论上可以得到s=C_n⁴+C_n³个定位结果，对这s个定位结果进行聚类分析，不含噪声或噪声较小的数据得到的定位结果间聚合性较好，将被聚成一类，噪声较大的数据产生的定位结果将被分成另外一类或几类。算法描述如下：①将n条原始数据进行3站和4站时差定位得到解析解R_i=(B_i，L_i)，其中，L_i，B_i为定位结果的经度和纬度，定位结果集合V={R₁，R₂，…，R_s}；②对集合V中的定位结果进行聚类分析，分成o类(V₁，V₂，…，V_o)；③选出聚合性最好的集合V_c；④求出集合V_c的簇中心R_c=(B_c，L_c)。

2.2 网格搜索法

网格搜索法通过将探测网的覆盖区域划分成若干个网格，将每个网格的位置作为初始位置代入适应函数χ²计算，以χ²最小的格点位置作为最优解。为了提高搜索效率，采取多级网格划分的策略，在搜索过程中不断更新解空间区域与迭代步长，使搜索在初始阶段具有较好的全局搜索能力，在后期搜索阶段具有较好的局部搜索能力。搜索采用方位角和时差法联合误差建立适应函数：

(1)

式(1)中，t_i和α_i为第i个站点接收到闪电的时间与方位信息，t_u为当前计算得到的闪电发生时间，d_i和θ_i是通过高斯平均引数法大地主题反算得到的大地线长度和方位角，w₁和w₂为权重系数。网格搜索方法具体步骤如下：①确定初始区域搜索空间，取椭球定位和聚类定位得到的初定位结果邻近区域作为搜索空间；②在第1级搜索时，将初始搜索区域分成若干个网格，计算每个格点的经纬度坐标；③计算相邻格点经纬度的差值；④计算各个格点的适应函数，取最小适应函数值对应的格点作为当前最优位置，以该格点为中心构造下一级搜索区域；⑤重复上述过程，当达到收敛条件后停止搜索得到最优雷电位置。

3 仿真分析 3.1 实例仿真 3.1.1 雷击实例

文献[8]和文献[9]给出了某地雷电监测网中的某雷击实例，包括各探测站的地理坐标及其探测到的到达时间，利用本文算法对该实例进行定位。任取6条回击数据中的3条或4条数据进行定位，对得到的35个定位结果进行聚类分析，聚合性较好的簇中心对应的纬度和经度分别为39.339°N和114.380°E，聚类效果如图 2所示，再通过网格搜索法得到最终定位结果纬度和经度分别为39.339°N和114.382°E，实际雷击点坐标的纬度和经度分别为39.340°N和114.385°E，两者相距84.09 m，验证了本文算法实现雷电定位的可行性。

从传统方法以及文献[8]提出的PSO方法与文献[9]提出的MFOA方法的定位结果看，TDOA法的定位结果误差最大，达到了1314 m；Taylor级数展开法作为经典的迭代算法，误差为881 m；PSO方法和MFOA方法均较大提升了定位精度，误差均小于500 m，分别为395 m和362 m，本文提出的定位算法DG-LLA误差为84 m, 验证了该算法实现雷击定位的有效性。

3.1.2 人工引雷实例

从2006年开始, 中国气象科学研究院每年夏季在广州从化进行人工引雷试验，这为检验闪电定位系统的精度提供了很好的对比观测数据。本文选取了3次人工触发闪电回击作为研究对象，根据国家业务闪电定位网各探测子站的原始记录重新进行定位计算，通过对比触发闪电实际回击位置获得定位误差的评估。人工引雷试验场位于广州从化北部光联村，引流杆(即实际回击位置)的纬度和经度分别为23.639°N, 113.595°E。表 1给出了3次人工触发闪电回击定位结果对比，根据保密法相关规定，本文对实例中的到达时间进行了技术处理。

定位结果与2018年从化3次人工触发闪电回击实际位置对比表明，本文算法DG-LLA定位结果比国家雷电监测网更精确。国家雷电监测网定位平均误差为153.47 m，采用DG-LLA定位后，定位误差降低至104.18 m，比原来降低32.2%。本文计算中所用的实例样本量较少，下一步将结合区域仿真进行算法验证。

3.2 区域仿真 3.2.1 仿真方法

为了定量评估新算法DG-LLA的区域定位误差，选取国家雷电监测网中湖北省雷电监测站网作为仿真区域，采用蒙特卡罗法对仿真区域进行计算机模拟^[17]，站网共包含13个探测站。在模拟过程中，在研究区域内任意给定一个闪电位置，其到各个探测站的准确时间可知，增加平均值为0、均方差为1 μs的正态分布随机误差。将研究区域等分成若干个网格，每个网格中模拟发生100个闪电，100个闪电定位误差的平均值认为是该网格点的定位误差。定位误差计算主要包括以下步骤：①计算闪电从网格中心到精度范围内探测站的到达时间；②在到达时间上增加随机时间误差，模拟闪电的真实传输情况；③利用定位算法对闪电进行定位；④计算定位结果与真实位置之差，即闪电的定位误差。

3.2.2 仿真结果

图 3为湖北区域定位误差分布对比图。由图 3可知，理论定位误差与距探测网远近呈正相关，TDOA法和Taylor级数展开法误差较大，最大误差分别达到了10.5 km和6.2 km，区域均方根误差分别为982 m和668 m；在定位计算中加入DBSCAN方法定位误差明显减小，最大定位误差为2.6 km，区域均方根误差为406 m；新算法DG-LLA的区域均方根误差减小为349 m。

图 4为时间误差标准差从0.1 μs增大至1 μs时区域定位均方根误差变化。由图 4可以看到，区域内定位误差随所加时间误差的增大而增大；当所加时间误差较小时，定位效果改善不明显，随着时间误差增大，本文定位算法DG-LLA优于传统定位算法，抗时间误差干扰能力与定位精度均有明显提高。

为了评估定位算法的效率，对不同探测站数量时各个算法模块的运行时间做了统计分析。结果显示：Taylor级数展开法平均耗时约41 ms；DBSCAN法效率较高，运行时间随探测站数量增多而增多，平均耗时约12 ms；DG-LLA运行时间平均耗时约32 ms，相较于Taylor级数展开法运行效率得到了提升，且定位更精确。

4 效果评估 4.1 定位数据及质量控制

国家雷电监测网采用ADTD闪电定位系统，采用时差法和定向时差联合法进行定位，仅能探测地闪，闪电监测网单站的探测范围为150 km，其探测率可达80%~90%，定位精度在500 m范围内^[18]。回击数据要素包括地闪回击时间、位置、极性、峰值强度等。本研究选取2018年8月1日的国家雷电监测网数据。

一般需要先对定位数据进行质量控制，然后进行对比分析^[19-21]。地闪定位系统直接探测回击，一次地闪通常包含有多次回击，需要将一定的时间和空间内的多次回击归为一次地闪。国家雷电监测网通过类似美国国家闪电监测网(NLDN)对闪击归类的方式进行归闪：首先确定首次回击，若继后回击距首次回击500 km内且时间间隔在3 ms内则被认为是一次闪电。如果一次回击同时可以归在多个闪电过程中，则归到距离最近的那个。

观测表明：正地闪可通过延伸很长的云闪分支开始，云闪中具有小峰值的电流符合正地闪标准而被闪电定位系统误作正地闪^[22]。本文仅统计10~500 kA的正地闪以及0~500 kA的负地闪，国家雷电监测网中有很多电流大于500 kA的闪电，一般被认为误差较大^[23-25], 对比时不予统计。

4.2 闪电次数及时间分布对比

2018年8月1日国家雷电监测网共接收到回击数据1276870条，对这些回击数据进行归闪并采用DG-LLA进行重新定位，与国家雷电监测网定位结果进行对比分析。结果显示：DG-LLA得到的正地闪次数(5966)、负地闪次数(158410)均比国家雷电监测网多(正地闪4982次，负地闪133624次)。两者得到的地闪极性接近，国家雷电监测网正地闪占总地闪百分比(3.5%)稍低于DG-LLA(3.6%)。国家雷电监测网正地闪电流(62.92 kA)大于DG-LLA正地闪电流(48.43 kA)，负地闪电流(38.61 kA，取绝对值，全文同)大于本文算法负地闪电流(36.64 kA)。

图 5为两系统2018年8月1日地闪次数分布对比。图 5a正地闪分布对比显示，两系统在13:00—14:00(北京时，下同)以及16:00—17:00探测到的次数较多，不同地方在于DG-LLA得到的正地闪次数较国家雷电监测网更多。由图 5b负地闪分布可知，12:00—18:00探测到的负地闪均较多，该时段为负地闪的高发时段。不同的是DG-LLA得到的负地闪次数较国家雷电监测网更多，且在15:00—17:00差异较大。初步分析两系统出现差异的原因主要有两个：①DG-LLA对回击数据进行重新归闪，归闪方式存在差异，国家雷电监测网通过实时归闪定位计算，存在部分数据未能有效利用，而该算法对原始所有回击数据重新排序归闪，数据利用率增加；②国家雷电监测网归闪成功的回击数据未能成功定位而DG-LLA成功定位，因此，增加了部分定位结果。

4.3 闪电空间分布对比

图 6为2018年8月1日正地闪空间分布对比，由图 6可以看到，该日地闪主要分布在中南部区域，国家雷电监测网监测到4982次正地闪，在东部和北部地区出现了部分偏离探测网的地闪，分布较为分散；DG-LLA得到正地闪5966次，离散点变少，分布更为集中。

图 7为2018年8月1日负地闪空间分布对比，由图 7可以看到，该日监测到的负地闪数量较多，分布范围较广，国家雷电监测网监测到133624次负地闪，采用DG-LLA定位后，得到正地闪158410次，离散点变少，分布更为集中。

不考虑两站定位，不同探测站数量的定位效果对比如表 2所示，包括定位数量，参与定位回击数据数量及定位结果与测站间的平均距离。

由表 2可以看到，在不同探测站数时，DG-LLA的定位数量均比国家雷电监测网偏多，参与定位的回击数据量也增多，回击数据利用率从41.5%提高至54.6%；从定位结果与探测站间的平均距离看，DG-LLA的定位结果距站网更近，定位结果更真实有效。总体而言，国家雷电监测网在回击数据量较少或回击数据存在较大观测误差时，定位算法的鲁棒性和稳定性较差，得到的定位结果分布较为离散，部分定位结果距探测网较远；DG-LLA能有效避免探测站数量和观测误差对定位计算的干扰，提高了定位计算的稳定性和准确性。

4.4 区域空间分布对比

2018年8月1日09:00—16:00湖北区域定位分布如图 8所示。在该时间段内，湖北区域监测地闪次数较多。国家雷电监测网共监测到21265次地闪，略少于DG-LLA的23664次。

由图 8可以看到，两系统总体分布情况大致相同，DG-LLA新增了一些国家雷电监测网未检测到的地闪，分布范围更广。为了进一步验证DG-LLA定位结果空间位置的可靠性，通过雷达回波和地闪定位结果的叠加对比图进行效果验证(如图 9所示)。

雷达组合反射率因子^[26]以10 min为单位生成文件，雷达时间为扫描结束时间，而地闪定位结果时间精度精确到小数点后7位，以整点时刻去匹配定位结果并不可靠，所以通过叠加当前时刻前后1 min内的地闪验证定位效果。2018年8月1日09:10国家雷电监测网探测到地闪8次，部分闪电周围回波较弱，组合反射率因子为10~20 dBZ，DG-LLA得到地闪15次，有4次地闪组合反射率因子超过30 dBZ; 14:20探测到的地闪数量较多，国家闪电监测网32次，DG-LLA得到117次，在东南区域出现了较多新增定位结果且区域回波较强，组合反射率因子最大达到60 dBZ；15:00—16:00探测到的地闪集中在西南区域，国家雷电监测网探测到地闪66次，DG-LLA得到172次，从15:10和15:40叠加图(图 9)可以看到，DG-LLA定位周围区域有回波出现，组合反射率因子为20~30 dBZ。通过雷达组合反射率因子叠加对比可知，DG-LLA部分新增定位结果具有一定的真实性和准确性。

叠加图能粗略验证定位结果空间位置的可靠性，为了对定位结果进行定量评估对比，选取国家闪电监测网与DG-LLA差异最为明显的区域统计雷达组合反射率因子与发生雷电数量。对比结果如表 3所示。

从选取的6个时间点的雷达组合反射率因子与定位叠加对比效果(表 3)可以看到，DG-LLA新增定位结果区域均在该闪电发生时刻出现较强的回波。

2018年8月1日14:20在29°N~30°N, 113°E~114°E区域内雷达组合反射率因子平均值为42.7 dBZ，国家雷电监测网探测到12次闪电，DG-LLA得到45次闪电，新增定位真实性较高。

5 结论

本文提出一种基于DBSCAN和网格搜索的雷电定位算法DG-LLA，通过仿真分析和实际定位效果对比，得到以下结论：

1) 雷击实例仿真结果说明基于DBSCAN的雷电定位算法能有效识别噪声数据，增强了抗误差干扰的能力。

2) 区域仿真结果表明：基于DBSCAN和网格搜索的雷电定位算法提高了迭代计算在解空间的局部搜索能力和全局搜索能力，克服了传统迭代算法易发散、优化算法陷入局部最优等缺点，稳定并精确求解出雷击点。

3) 通过与国家雷电监测网定位结果对比分析，新算法DG-LLA取得了更好的定位效果。在相同回击数量的条件下，DG-LLA得到了更多的定位，提高了回击数据的利用率，新增定位结果周围雷达回波特征较强，定位效果较好。

致谢: 感谢中国气象科学研究院灾害天气国家重点实验室雷电科研团队提供相关的人工触发闪电数据。