本发明属于气象探测，具体涉及一种基于雷声的闪电通道定位方法和定位系统。
背景技术：
1、雷暴气候是发生在中小尺度对流体统中的常见天气现象，雷暴天气往往伴随着闪电的活动，可引起森林和油库火灾、造成供电及通讯信息系统故障或损坏，对航天航空、矿山及一些重要而敏感的高技术装备等具有重大威胁。八十年代以后，闪电引起的危害显著增加，特别是与高新技术关系密切的领域，如航空航天、国防、通讯、电力、计算机、电子工业等由于广泛应用对闪电电磁干扰极为敏感的大规模及超大规模集成电路致使遭雷击的几率大大增加。所以闪电辐射源定位技术显得尤为重要。
2、闪电辐射源定位技术的发展可以有助于雷暴的预测，也可以为雷电防护设计起到针对性的理论指导。同时，通过闪电辐射源不同时间段的叠加，能够完整的描绘闪电通道《alow frequency 3d lighming mapping network in north china》(《atmosphericresearch》(2022))，可以清晰重构闪电的传输，基于此可以研究闪电的物理特征《characteristics of impulsive currents superimposing on continuous/continuingcurrent of rocket-triggered lightning》(《ieee transactions on electromagneticcompatibi1ity》(2019))，研究雷暴云内电荷极性分布《evolution of the chargestructure and lightning discharge characteristics of a qinghai tibet plateauthunderstorm dominated by negative cloud-to-ground flashes.》(《journal ofgeophysical research》(2020))，为闪电物理研究提供了重要的技术支持。
3、相较于基于闪电电磁辐射源定位，雷声定位的定位范围更远，声音信号也比电磁场信号拥有更强的扛干扰能力，不容易受到复杂地形，电磁环境噪声的影响。而且，闪电电磁场辐射源定位需要更多的观测站布局，或者更复杂的计算方法，而雷声定位则相对较易实现。总之，雷声定位在闪电定位技术中起到不可替代的作用。
4、但目前雷声定位存在一些问题：
5、首先，环境影响与信号传播特性不确定性，雷声定位算法依赖于声波在不同介质中的传播速度来计算闪电的位置。然而，环境因素如大气湿度、温度和风速等会对声波的传播速度产生影响，从而引发定位误差。此外，由于雷电和声波在不同介质中的传播特性具有一定的不确定性，可能导致定位结果的不准确性《传播介质对阵列雷声定位结果的影响模拟研究》(《电声技术》(2021))。其次，多径效应和噪声干扰《基于vtr-plfm扩频调制的移动水声通信方法》(《西北工业大学航海学院》(2021))，雷声在传播过程中容易发生多次反射、绕射和衍射，使得多条路径上的声波同时到达接收器，引起定位误差。最后，定位精度受限，通常采用时间差法(time difference of arrival，toa)进行定位计算，传感器的布局和时差测量精度都会对定位精度产生一定的限制。
6、综上所述，目前亟需要设计一个雷声定位方法以提高定位精度及效率。
技术实现思路
1、解决的技术问题：本发明提出了一种基于雷声的闪电通道定位方法和定位系统，能够精确计算雷声源的位置和处理定位的不确定性，提高定位效率。
2、技术方案：
3、一种基于雷声的闪电通道定位方法，所述闪电通道定位方法包括：
4、s1，采用多个麦克风构成麦克风阵列，收集每个麦克风接收到的雷声信号，对收集的雷声信号进行预处理；
5、s2，使用交叉相关函数方法进行时间差计算，获得第一时间差估计值；再使用互信息方法得到第二时间差估计值；采用监督学习的方法，将计算得到的第一时间差估计值和第二时间差估计值作为训练标签，雷声信号的特征作为输入，对神经网络模型进行训练，得到时间差的预测值；
6、s3，对第一时间差估计值、第二时间差估计值和时间差的预测值进行融合，得到到达时间差；
7、s4，基于雷声源与每个麦克风之间的到达时间差和声速，计算得到雷声源与每个麦克风之间的距离；
8、s5，引入多普勒效应对计算得到的雷声源与每个麦克风之间的距离进行修正；
9、s6，综合修正后的雷声源与每个麦克风的距离，计算得到雷声源的位置坐标。
10、进一步地，步骤s1中，采用5个麦克风构成麦克风阵列，其中一个麦克风位于中心位置，其余麦克风围绕中心位置的麦克风形成一个正方形；以麦克风阵列中心作为坐标系原点，五个麦克风的坐标分别为：(0，0，a)、(a，0，0)、(-a，0，0)、(0，a，0)和(0，-a，0)。
11、进一步地，步骤s1中，对收集的雷声信号进行预处理的过程包括以下步骤：
12、对每个麦克风接收到的雷声信号si(t)进行傅立叶变换，转换到频域：
13、si(f)＝∫[si(t)·e-2πift]dt
14、其中，si(f)是第i个麦克风在时刻t接收到的雷声信号si(t)的傅立叶变换，f是频率，∫表示积分；
15、设计带通滤波器b(f)对转换到频域的雷声信号进行去噪处理，得到滤波后的频域信号si′(f)：
16、si′(f)＝b(f)*si(f)
17、
18、其中，fi是滤波器的低频截止频率，fh是滤波器的高频截止频率；
19、使用逆傅立叶变换得到滤波后的时域信号si′(t)：
20、s′i(t)＝f-1[s′i(f)]
21、其中，f-1[·]表示逆傅立叶变换。
22、进一步地，步骤s2中，采用下述公式计算得到第i个麦克风和第j个麦克风之间的交叉相关函数rij(τ)：
23、rij(τ1)＝∫[s′i(t)*s′j(t+τ1)]dt
24、其中，s′i(t)是第i个麦克风t时刻和第j个麦克风t+τ1时刻滤波后的雷声信号，s′j(t+τ1)是第j个麦克风t+τ1时刻滤波后的雷声信号，τ是时间延迟；通过找到使rij(τ1)达到最大值的τ1来获得第一时间差估计值
25、第一时间差估计值的计算流程包括以下步骤：
26、s201：雷声信号的采集
27、使用两个麦克风分别记录雷声；假设这两个麦克风分别为mi和mj；
28、s202：预处理
29、对每个麦克风接收到的雷声信号si(t)和sj(t)进行预处理；预处理包括傅立叶变换、带通滤波和逆傅立叶变换；
30、s203：交叉相关函数的计算
31、计算两个滤波后的雷声信号si(t)和sj(t)之间的交叉相关函数rij(τ)；交叉相关函数的计算公式为：但在实际操作中，因为信号是离散采样的，所以使用的是离散信号的交叉相关函数的计算公式：其中，k是离散时间延迟的索引，n是离散时间索引。
32、s204：估计时间差
33、通过寻找交叉相关函数rij(τ)的最大值来确定时间差τ；即：找到rij(k)的峰值所对应的k值，这个k值即为第一时间差估计值；
34、采用下述公式计算得到第i个麦克风和第j个麦克风之间的互信息：
35、lij(τ2τ)∑[pij(t，t+τ2)*log 2(pij(t，t+τ2))/(pi(t)*pj(t+τ2))]
36、其中，pij(t，t+τ2)是雷声信号在时间t和t+τ2的联合概率密度函数，pi(t)和pj(t+τ2)分别是雷声信号在时间t和t+τ2的概率密度函数，∑表示求和；通过找到使lij(τ2)达到最大值的τ2来获得第二时间差估计值
37、第二时间差估计值的计算流程包括以下步骤：
38、s211：信号的预处理
39、与计算第一时间差估计值相同，首先采集两个麦克风mi和mj接收到的雷声信号si(t)和sj(t)并进行预处理。
40、s212：计算互信息
41、互信息的计算需要先估计雷声信号的概率密度函数。应用核密度估计来实现。kde是通过将数据点周围的每个点平滑处理来工作的，这涉及到使用核函数，本发明采用高斯核。选择核函数kernel和带宽h。对于每个数据点xi，计算核函数在x的值；将所有核函数值相并归一化以得到概率密度函数的估计：核函数kernel选用高斯函数，它的形式是：其中u是观测值xi到点x的归一化距离除以带宽h：互信息的计算公式为：
42、
43、式中，p(si(t)，sj(t+τ))是在时间t和t+τ下si和sj的联合概率密度函数，p(si(t))和p(sj(t+τ))分别是si和sj在各自时间下的边缘概率密度函数。
44、s213：时间差的估计
45、找到使iij(τ)达到最大值的τ，这个τ即为第二时间差估计值；采用网络搜索法：网格搜索是一种穷举搜索方法，它在预定义的参数范围内系统地遍历多个参数组合；计算步骤具体包括：定义τ的搜索范围，比如-t到t，其中t是可能的最大时间延迟；定义步长δτ，这将决定网格的密度；对于每个t值，计算iij(τ)；找到使iij(τ)最大化的τ值。
46、进一步地，步骤s2中，采用监督学习的方法，将计算得到的第一时间差估计值知第二时间差估计值作为训练标签y，雷声信号的特征x作为输入，训练神经网络模型，得到时间差的预测值
47、假设神经网络模型参数为θ，损失函数为l，训练过程包括：
48、在旧数据上训练神经网络模型，获得模型参数θold，使相应的旧的损失函数lold达到最小值；
49、计算损失函数lo1d关于参数θ的二阶导数，获取fisher信息矩阵f，fisher信息矩阵f是损失函数的曲率，用于决定参数在哪些方向上的变化会对损失函数产生大的影响；
50、当新数据到来时，在新的损失函数lnew上训练模型，并加入一个额外的惩罚项，惩罚项用于抵抗参数在旧任务上有重要影响的方向的改变：
51、
52、
53、其中，λ是超参数，用于决定对旧任务的保护程度；e[]表示期望，是关于参数θ的梯度，p(y|θ)是模型在参数θ下，对旧任务标签y的预测概率。
54、进一步地，步骤s3中，采用下述公式对第一时间差估计值、第二时间差估计值和时间差的预测值进行加权平均，融合得到到达时间差：
55、
56、其中，和分别是第一时间差估计值、第二时间差估计值和时间差的预测值，wccf、wmi和wml是相应的权重，与时间差的计算方法的可靠性相关，且wccf+wmi+wml＝1。
57、进一步地，步骤s4中，采用下述公式计算得到雷声源与每个麦克风之间的距离：
58、
59、式中，是雷声源与第i个麦克风之间的到达时间差，v是声速；di是雷声源与第i个麦克风之间的距离；
60、步骤s5中，采用下述公式修正雷声源与每个麦克风之间的距离di：
61、dcorrected_i＝di-δf*δti/(2*f0)
62、δf＝2*vs*f0/v
63、其中，f0是雷声信号的初始频率，δf是由多普勒效应引起的频率改变，vs是雷声源的速度。进一步地，步骤s6中，使用修正后的距离dcorrected_i，建立以下方程组：
64、dcorrected_i2＝(x-xi)2+(y-yi)2+(z-zi)2
65、其中，(x，y，z)是雷声源的坐标，(xi，yi，zi)是第i个麦克风的坐标；
66、通过最小二乘法求解这个方程组，得到雷声源位置的点估计值。
67、进一步地，所述定位方法还包括：引入贝叶斯推理，计算声源位置的后验概率分布；具体包括以下步骤：
68、s71，假设定位误差服从高斯分布，观测模型p(d|x)为高斯分布，d表示观测到的数据，即修正后的距离dcorrected_i，x表示声源位置；在高斯分布的假设下，有：
69、
70、式中，f(x)表示根据声源位置x计算出的理论值，即声源位置(x，y，z)与麦克风i的距离d，σ表示观测噪声的标准差，通过以下公式计算理论值f(x)：
71、
72、s72，假设先验概率分布p(x)为均匀分布，在均匀分布的假设下，有：
73、p(x)＝1/(b-a)
74、式中，a和b分别代表最小边界值和最大边界值；
75、根据贝叶斯定理，通过以下公式求解后验概率分布p(x|d)：
76、
77、s73，使用步骤s6中得到的声源位置的点估计值(x0，y0，z0)作为初始点，通过马尔科夫链蒙特卡洛算法采样后验概率分布，采样过程包括以下步骤：
78、s731，选择初始声源位置(x0，y0，z0)和步长j，步长决定每次迭代时参数移动的距离；在当前位置的基础上加上一个小的随机扰动，生成一个新的声源位置(x′，y′，z′)：
79、x′＝x+n(0，j)
80、y′＝y+n(0，j)
81、z′＝z+n(0，j)
82、式中，n(0，j)表示均值为0，标准差为j的正态分布：
83、s732，计算接受概率，接受概率是新旧声源位置的后验概率之比：
84、a＝min(1，p(x′|d)/p(x|d))
85、其中p(x|d)和p(x′|d)分别是旧声源位置和新声源位置的后验概率；
86、s733，根据接受概率决定是否接受新的声源位置：从均匀分布u(0，1)中生成一个随机数u，如果u＜a，则接受新的声源位置，将(x′，y′，z′)设置为新的声源位置；否则，保持当前声源位置不变；
87、s734，重复步骤s732至s733，直到满足停止条件；
88、s74，生成n个样本(x1，y1，z1)，(x2，y2，z2)，...，(xn，yn，zn)，通过计算样本的平均值，得到声源位置的最大后验估计：
89、
90、其中，∑表示求和，k＝1，2，...，n。
91、本发明公开了一种基于雷声的闪电通道定位系统，所述闪电通道定位系统包括：
92、麦克风阵列，由多个分布在不同地点的麦克风组成，用于实时收集雷声信号；
93、数据预处理模块，用于收集每个麦克风接收到的雷声信号，对收集的雷声信号进行预处理；
94、数据存储模块，用于存储预处理后的雷声信号；
95、定位计算模块，用于采用如前所述定位方法对数据存储模块中存储的雷声信号进行处理，计算得到雷声源的位置坐标；
96、结果展示和警告模块，用于在地图上标示出计算得到的雷声源的位置坐标，并且对超出预设安全范围的雷声源位置，触发警告信息。
97、有益效果：
98、第一，本发明的基于雷声的闪电通道定位方法和定位系统，集成了交叉相关函数法(ccf)、互信息法(mi)以及基于机器学习的预测模型，三者相互协作提高了定位的准确度。交叉相关函数法(ccf)为提供一个粗略的时间差估计值；互信息法(mi)进行更精细的时间差计算；基于机器学习的预测模型则进一步提高了时间差的计算精度。这种综合使用多种方法的计算到达时间差的策略，可以在不同程度和不同角度上处理信号数据，优化定位精度。
99、第二，本发明的基于雷声的闪电通道定位方法和定位系统，在计算雷声源到麦克风的距离时，考虑到了雷声源的速度对频率的影响，即多普勒效应，这使得计算的距离更为准确，进一步提高了雷声源定位的精度。
100、第三，本发明的基于雷声的闪电通道定位方法和定位系统，引入了贝叶斯推理，利用马尔科夫链蒙特卡洛(mcmc)方法考虑定位的不确定性，通过计算声源位置的后验概率分布得到声源位置的最大后验估计，这使得定位结果更为精确，能够在考虑了实际测量误差和模型误差等因素的前提下，得到更为可靠的定位结果。
1.一种基于雷声的闪电通道定位方法，其特征在于，所述闪电通道定位方法包括：
2.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s1中，采用5个麦克风构成麦克风阵列，其中一个麦克风位于中心位置，其余麦克风围绕中心位置的麦克风形成一个正方形；以麦克风阵列中心作为坐标系原点，五个麦克风的坐标分别为：(0，0，a)、(a，0，0)、(-a，0，0)、(0，a，0)和(0，-a，0)。
3.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s1中，对收集的雷声信号进行预处理的过程包括以下步骤：
4.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s2中，采用下述公式计算得到第i个麦克风和第j个麦克风之间的交叉相关函数rij(τ)：
5.根据权利要求4所述的基于雷声的闪电通道定位方法，其特征在于，步骤s2中，采用监督学习的方法，将计算得到的第一时间差估计值和第二时间差估计值作为训练标签y，雷声信号的特征x作为输入，训练神经网络模型，得到时间差的预测值
6.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s3中，采用下述公式对第一时间差估计值、第二时间差估计值和时间差的预测值进行加权平均，融合得到到达时间差：
7.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s4中，采用下述公式计算得到雷声源与每个麦克风之间的距离：
8.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，步骤s6中，使用修正后的距离dcorrected_i，建立以下方程组：
9.根据权利要求1所述的基于雷声的闪电通道定位方法，其特征在于，所述定位方法还包括：引入贝叶斯推理，计算声源位置的后验概率分布；具体包括以下步骤：
10.一种基于雷声的闪电通道定位系统，其特征在于，所述闪电通道定位系统包括：
