Lightning low-frequency three-dimensional positioning method for complex environment

Abstract

The invention discloses a lightning low-frequency three-dimensional positioning method for a complex environment, and the method comprises the following steps: constructing a signal processing frame which comprises a graph learning module, a signal path optimizer, a signal decision agent module and a parameter updating module to carry out the preprocessing of a lightning signal; performing pulse matching on the preprocessed lightning signal to obtain initial positioning; and fitting the preliminary positioning through an improved particle swarm optimization (PSO) algorithm to obtain a final positioning result. The method focuses on three key technologies of signal processing, pulse matching and initial result optimization in lightning low-frequency positioning, and through deep research and innovation of the three key links, the precision and reliability of lightning low-frequency positioning in a complex environment are improved.

Description

Technical Field

[0001] The invention relates to the technical field of signal processing, and in particular to a low-frequency three-dimensional lightning positioning method for complex environments. Background Art

[0002] The positioning and imaging technology of lightning radiation sources is an important technology in the development and transmission process of lightning channels. It is not restricted by cloud occlusion (which is often encountered by optical methods) and can depict the discharge process in the cloud. In recent years, by upgrading hardware equipment, improving acquisition schemes, and developing signal processing and positioning algorithms, VLF radiation source positioning has also acquired a certain ability to depict the in-cloud process of lightning activity, such as the initial breakdown process, K process, and other channel re-breakdown phenomena, which have been well studied.

[0003] Time Difference of Arrival (TDOA) is an algorithm widely used in lightning location. Its basic principle is to invert the spatial position of the radiation source by analyzing the differences in the arrival times of electromagnetic signal pulses from the same radiation source observed at different stations. Therefore, the positioning accuracy of the TDOA algorithm depends to a large extent on the accuracy of pulse matching between the stations. However, in practical applications, the TDOA algorithm faces multiple challenges: the electromagnetic environment of each station is complex and changeable, introducing various interferences and noises; low-frequency electromagnetic signals have significant attenuation characteristics during long-distance transmission, resulting in weak signal strength received by distant stations; at the same time, multiple discharge breakdown events may occur in the cloud at the same time, generating similar pulse signals. These factors work together to make it difficult for each station to accurately identify and match pulse signals from the same radiation source, thus affecting the positioning accuracy and reliability of the TDOA algorithm.

[0004] In order to improve the accuracy of pulse matching and the signal-to-noise ratio of signals, researchers have developed a variety of advanced signal processing methods. Empirical Mode Decomposition (EMD) has been widely used in the preprocessing stage of lightning signals. EMD technology adaptively decomposes complex non-stationary signals into several intrinsic mode functions (IMFs), realizing multi-scale analysis and feature extraction of lightning signals. However, these advanced methods still face challenges when processing complex multi-frequency and multi-scale lightning signals, such as the time window selection of Hilbert transform, the endpoint effect of EMD processing signals, and the mode mixing problem of IMF decomposition, which need to be carefully considered and further optimized in practical applications. Frequency domain filtering is another commonly used lightning signal processing method. Frequency domain filters are mainly used to remove noise components in lightning signals. By designing appropriate filter parameters, such as cutoff frequency and order, noise within a specific frequency range can be effectively suppressed and the signal quality can be improved. However, lightning signals have complex spectral and time domain characteristics, and the electromagnetic radiation generated by different discharge processes has significant differences in the spectrum and time domain. This complexity results in certain limitations on the effect of frequency domain filters, making it difficult to always maintain optimal filtering performance.

[0005] Regardless of the matching method used, the signal-to-noise ratio of the data is always a key factor affecting the accuracy of lightning location. With the development of machine learning and big data technology, intelligent information processing algorithms based on massive lightning observation data training are becoming increasingly mature. This type of algorithm combines cutting-edge technologies such as signal processing, data mining, and machine learning, and intelligently analyzes and processes complex lightning data through adaptive noise suppression, signal feature enhancement, and multi-sensor information fusion. These emerging technologies have shown significant advantages in improving the accuracy of low-frequency three-dimensional lightning location.

[0006] On the other hand, the selection of lightning location algorithm is also one of the important factors affecting the positioning accuracy. TDOA establishes a nonlinear mathematical model based on the signal arrival time difference to describe the relationship between the signal source coordinates and the sensor position. Due to the nonlinear characteristics of the model, iterative optimization algorithms, such as Newton's method or gradient descent method, are usually used to determine the precise location of the lightning source by minimizing the residual. The CHAN algorithm can improve the computational efficiency by linearizing the TDOA equation. However, since the weight of the equation group is not considered, it is easily affected by large errors of individual measuring stations, resulting in a significant decrease in positioning accuracy. To further optimize the positioning results, researchers introduced the Levenberg-Marquardt (LM) algorithm, which achieves efficient and robust parameter estimation by locally linearizing the nonlinear problem and introducing a damping factor to balance the convergence speed and stability. Many studies have shown that the use of the LM algorithm to optimize the TDOA positioning results can significantly improve the positioning accuracy and reduce the error to the order of tens of meters horizontally to hundreds of meters vertically. However, the LM algorithm is based on the optimization of the local quadratic approximation model, but when the initial point is far away from the global optimal solution, it may converge to a suboptimal solution or a local minimum point. This problem is particularly obvious when the initial error is large or the objective function is non-convex, which limits the global optimization ability of the algorithm. Summary of the invention

[0007] To solve the above problems, the present invention provides a low-frequency three-dimensional lightning positioning method for complex environments. The specific technical solution is as follows:

[0008] The method for low-frequency three-dimensional lightning positioning in complex environments includes the following steps:

[0009] A signal processing framework is constructed to pre-process lightning signals, wherein the signal processing framework includes: a graph learning module, a signal path optimizer, a signal decision agent module and a parameter update module;

[0010] Pulse matching is performed on the preprocessed lightning signal to obtain a preliminary location;

[0011] The final positioning result is obtained by fitting the preliminary positioning through improving the particle swarm optimization (PSO) algorithm.

[0012] Furthermore, the graph learning module is used to capture the spatial characteristics of lightning signals and output comprehensive site, signal and time embeddings. Given the embedding representation of the original signal, the embedding representation of the signal is iteratively updated. The specific formula is:

[0013]

[0014] in, and They represent the embedding of the current site u, the embedding of the event et at the current time t, and the embedding of the predicted event e't+1 at the next time step t+1, respectively. These embeddings are obtained through the graph learning module and reflect the site characteristics, current event characteristics, and predicted event characteristics, respectively. The formula uses these embeddings to calculate the transition probability from the current site to the next site, which is used to optimize the signal propagation path.

[0015] A u , and They represent the adjacency matrix related to the site u, the current event et and the predicted event e't at the lth level, respectively. and Represent the weight matrices related to site u, current event et and predicted event e't respectively, and They represent the corresponding bias terms respectively, and σ is the nonlinear activation function;

[0016] The signal path optimizer selects the optimal signal path. The specific formula is:

[0017]

[0018] Among them, s t =(u,e t ) represents the current state, including site u and event e at time step t t , W P represents the path selection weight matrix, and Respectively represent the event e t and e t+1 Embedding representation in graph learning module;

[0019] The signal path optimizer generates an action t , the specific formula is:

[0020] a t =(e t →e t+1 ) (5)

[0021] Set reward r(s t ,a t ) is used to evaluate the quality of an action. The specific formula is:

[0022]

[0023] Among them, W r represents the reward weight matrix, N represents the number of samples in a 100 microsecond window, and h i represents the feature representation of the i-th sample, h signalThe embedded representation of the signal in the current state is obtained through the graph learning module, which contains the key features of the lightning electric field signal and is used to evaluate the current signal state and guide the decision-making process; represents the expected signal feature embedding, which reflects the expected impact of a specific action on the signal state and is used to evaluate the effect of the action and optimize the signal processing strategy; in reward calculation, h signal and The comparison can evaluate the impact of actions on signal processing, thereby guiding the system to learn the optimal signal processing strategy;

[0024] The signal decision agent module sets the action value function Q to evaluate action a t The expected utility of is:

[0025]

[0026] Among them, θ Q Represents the set of parameters of the action value function, which is used to learn and predict the long-term rewards of different state-action pairs. By adjusting these parameters, the system can optimize its decision-making strategy to obtain better signal processing effects; γ represents the discount factor;

[0027] In the parameter update module, the policy network parameters θ are included E and signal processing model parameters θ SP Update, the specific formula is:

[0028]

[0029] in, represents the policy network parameters at time t, Represents the E The gradient of Representative strategy In a given state s t Next select action a t The logarithmic probability of , α represents the learning rate, which controls the step size of parameter update. It balances learning speed and stability. Larger values ​​speed up learning but may be unstable, while smaller values ​​slow down learning but are more stable; Δθ SP represents the gradient of the signal processing model parameters, η represents the learning rate, It represents the expected utility of the system for the next state, which is used to evaluate the long-term impact of the current decision and guide the optimization of signal processing strategies;

[0030] The signal processing model parameter update rule is refined, and the specific formula is:

[0031]

[0032] in, represents the signal processing model parameters at time t.

[0033] Further, pulse matching is performed on the preprocessed lightning signal, specifically including:

[0034] EEMD is used to process pulse data. The specific steps are as follows:

[0035] Add a signal with a mean of zero and a variance of σ to the original signal x(t) 2 Gaussian white noise n i (t), and get the new signal x i (t) = x(t) + n i (t); for x i (t) Perform EMD decomposition to obtain a series of IMFs and residual terms; repeat the addition of noise and EMD decomposition M times, adding different noise each time; average the IMFs obtained from the M decompositions to obtain the final EEMD decomposition result;

[0036] Feature extraction is performed on the EEMD decomposition results. The specific steps are as follows:

[0037] The energy proportion and correlation of IMF are calculated respectively. The specific formula is:

[0038]

[0039] R min =min{R 1 ,R 2 ,…,R K} (13)

[0040] Among them, E ratio It represents the ratio of the cumulative energy of the first k IMFs to the total energy, E j represents the energy of the jth IMF, E total represents the total energy of all IMFs; R min represents the minimum correlation coefficient between the first k IMFs and the original signal, R 1 ,R 2 ,…,R K Represents the correlation coefficient between each IMF and the original signal;

[0041] Define the K value to indicate that the IMFs with the top K energy proportions and correlations are considered as characteristic IMFs. The specific formula is:

[0042] C(k)=(E ratio ≥E th )∧(R min ≥R th ) (14)

[0043] Among them, C(k) is a Boolean function. When the energy ratio is not less than the threshold Eth And the minimum correlation is not less than the threshold R th True when ; the minimum k value that meets the conditions is the optimal K value;

[0044] The cross-correlation function is used to evaluate the similarity between different IMFs. The specific formula is:

[0045] CCF ij (τ)=∑[imf i (t)·imf j (t+τ)] (15)

[0046] Among them, CCF ij (τ) represents the cross-correlation function between two stations i and j, τ represents the time offset, imf i (t), imf j (t+τ) represents the characteristic intrinsic mode function of the i-th station;

[0047] By calculating CCF ij The maximum value of (τ) determines the best matching position between the two IMFs, thereby achieving accurate time delay estimation. The specific formula is:

[0048] σ tij =arg max {CCF ij (τ)} (16)

[0049] Among them, σ tij represents the time delay between station i and station j.

[0050] Furthermore, the improved particle swarm optimization PSO algorithm specifically adopts a double particle swarm collaborative optimization strategy based on Tent mapping distribution, which specifically includes the following steps:

[0051] After initialization, the particles in Pop2 are sorted from small to large according to fitness;

[0052] Combine them to generate M subgroups, and the number of particles P in each subgroup is:

[0053]

[0054] Where N is the number of particles in Pop2 population;

[0055] The fitness value calculation formula for each subgroup is:

[0056]

[0057] in, represents the fitness of the mth subgroup, t represents the number of iterations, represents the fitness of the i-th particle in the m-th subgroup;

[0058] Different mutation capabilities are set according to different adaptability. The standard deviation of the t-th mutation operator is:

[0059]

[0060] in, represents the standard deviation of the t-th mutation operator, FitX max and FitX min Represent the maximum and minimum fitness values ​​respectively;

[0061] Pop1 population updates its speed according to the speed update formula. The specific formula is:

[0062]

[0063] Pop2 population updates its speed according to the speed update formula. The specific formula is:

[0064]

[0065] Among them, c 1 and c 2 represents the acceleration constant, ω represents the inertia weight, Represents the position vector of the d-th dimension of the optimal particle in the population Pop2. When the speed in Pop2 is less than the set threshold, speed escape will be performed;

[0066] Gaussian mutation is performed on the m subgroup particles in Pop2, and the velocity vector of the particles in the mth subgroup is updated as follows:

[0067]

[0068] Here, randn represents a random number uniformly distributed in the range [0,1].

[0069] According to one aspect of the present invention, a storage medium is provided, wherein instructions are stored in the storage medium. When a computer reads the instructions, the computer executes any one of the above-mentioned methods for low-frequency three-dimensional positioning of lightning in a complex environment.

[0070] According to another aspect of the present invention, there is provided an electronic device, comprising a processor and the above-mentioned storage medium, wherein the processor executes instructions in the storage medium.

[0071] Compared with the prior art, the present invention has the following beneficial effects:

[0072] 1. The present invention utilizes the spatiotemporal correlation of observation data from multiple stations to introduce the GNN algorithm to pre-process the synchronization signal and improve the denoising capability;

[0073] 2. The present invention introduces the EEMD algorithm to achieve the best matching effect by performing cross-correlation on the IMF that characterizes the pulse;

[0074] 3. The present invention designs a D-PSO algorithm to replace the LM algorithm to fit the initial solution, thereby solving the problem of dependence on the accuracy of the initial solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0076] Figure 1 It is a schematic diagram of a flow chart of an embodiment of the present invention;

[0077] Figure 2 A signal analysis system framework diagram of an embodiment of the present invention;

[0078] Figure 3 A schematic diagram of a particle position update strategy of a dual-group PSO algorithm according to an embodiment of the present invention;

[0079] Figure 4 A flow chart of a dual-group PSO algorithm according to an embodiment of the present invention;

[0080] Figure 5 This is a comparison diagram of positioning errors in an embodiment of the present invention;

[0081] Figure 6 A partial lightning signal diagram of an embodiment of the present invention;

[0082] Figure 7 A cross-sectional comparison diagram of three positioning methods according to an embodiment of the present invention;

[0083] Figure 8 It is a 3D comparison diagram of three positioning methods according to an embodiment of the present invention. DETAILED DESCRIPTION

[0084] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0085] The present invention adopts the lightning dual-time base data acquisition mode (ABA), which realizes the high sampling rate trigger acquisition of single-channel large-amplitude lightning signals and the low sampling rate continuous acquisition of complete signals, taking into account the integrity of the data while ensuring the efficiency of data acquisition. In order to improve the stability and reliability of the data acquisition system, the present invention sets a FIFO (First Input First Output) buffer to effectively temporarily store high-speed collected data to avoid data loss, and introduces dual timestamp technology to ensure the time synchronization of data from different measuring stations, providing a reliable time reference for lightning location. These measures provide high-quality data support for subsequent lightning location and physical process analysis.

[0086] Lightning signal preprocessing

[0087] The signal processing framework proposed in this invention consists of four closely cooperating modules: a graph learning module, a signal path optimizer, a signal decision agent, and an update interpretation model.

[0088] Figure 2 The figure shows the framework of the lightning location and signal analysis system based on graph learning and reinforcement learning proposed by the present invention. The system first collects the electric field signals generated by lightning through spatially distributed monitoring stations and transmits them to the graph learning module. In this module, the system converts the monitoring station network into a graph model representing the connectivity relationship of the stations through operations such as site aggregation, site embedding, signal embedding and time embedding. The site aggregation operation is responsible for summarizing the information of multiple monitoring stations; the site embedding operation maps the information of each station to a low-dimensional vector space; the signal embedding operation maps the electric field signal to a low-dimensional vector; and the time embedding operation is used to map the time information to a low-dimensional vector.

[0089] Based on the construction of the graph model, the signal path optimizer uses a two-step strategy to search for the optimal signal propagation path. The first step is to select the next hop node based on the current strategy; the second step is to continue searching according to the updated strategy until the optimal path is found. In order to strike a balance between exploring new paths and leveraging existing knowledge, the system introduces a counterfactual inference unit to evaluate the performance indicators of the current path. The decision-making agent compares the indicator with the preset benchmark value to decide whether to accept the current path. Once the optimal path is found, the system outputs the location coordinates of the lightning signal source and the key characteristic parameters related to it.

[0090] Finally, the signal processing model uses the acquired positioning and feature information to filter and decompose the original electric field signal, extract a more refined feature representation, and use it to update the graph model so that it can adapt to the dynamic changes of the environment. By optimizing the signal propagation path in the sensor network, the algorithm can overcome the blind spots and errors in lightning signal analysis and improve positioning accuracy and feature extraction capabilities.

[0091] The specific algorithm is as follows: The graph learning module is responsible for capturing the spatial characteristics of lightning signals and outputting comprehensive site, signal and time embeddings. Given the embedding representation h of the original signal (0) , the module iteratively updates the embedded representation of the signal through the following graph convolution formula:

[0092]

[0093] in, and They represent the embedding of the current site u, the embedding of the event et at the current time t, and the embedding of the predicted event e't+1 at the next time step t+1, respectively. These embeddings are obtained through the graph learning module and reflect the site characteristics, current event characteristics, and predicted event characteristics, respectively. The formula uses these embeddings to calculate the transition probability from the current site to the next site, which is used to optimize the signal propagation path.

[0094] A u , and They represent the adjacency matrix related to the site u, the current event et and the predicted event e't at the lth level, respectively. and Represent the weight matrices related to site u, current event et and predicted event e't respectively, and They represent the corresponding bias terms respectively, and σ is the nonlinear activation function;

[0095] The signal path optimizer selects the optimal signal path through the following state transition probabilities:

[0096]

[0097] Among them, s t =(u,e t ) represents the current state, including site u and event e at time step t t , W P represents the path selection weight matrix, and Respectively represent the event e t and e t+1 Embedding representation in graph learning module;

[0098] This step generates the action a t It can be expressed by the following formula:

[0099] a t =(e t →e t+1 ) (5)

[0100] The corresponding reward r(s t ,at ) is used to evaluate the quality of an action and is calculated as follows:

[0101]

[0102] Among them, W r represents the reward weight matrix, N represents the number of samples in a 100 microsecond window, and h i represents the feature representation of the i-th sample, h signal The embedded representation of the signal in the current state is obtained through the graph learning module, which contains the key features of the lightning electric field signal and is used to evaluate the current signal state and guide the decision-making process; represents the expected signal feature embedding, which reflects the expected impact of a specific action on the signal state and is used to evaluate the effect of the action and optimize the signal processing strategy; in reward calculation, h signal and The comparison can evaluate the impact of actions on signal processing, thereby guiding the system to learn the optimal signal processing strategy.

[0103] The action value function Q of the signal decision agent is used to evaluate action a t Expected utility:

[0104]

[0105] Among them, θ Q Represents the set of parameters of the action-value function, which is used to learn and predict the long-term rewards of different state-action pairs. By adjusting these parameters, the system is able to optimize its decision-making strategy to obtain better signal processing effects; γ represents the discount factor, which represents the current value of future rewards.

[0106] Policy network parameters θ E The update depends on the action value function and the policy gradient, as follows:

[0107]

[0108] in, represents the policy network parameters at time t, Represents the E The gradient of Representative strategy In a given state s t Next select action a t The logarithmic probability of , α represents the learning rate, which controls the step size of parameter update. It balances learning speed and stability. Larger values ​​speed up learning but may be unstable, while smaller values ​​slow down learning but are more stable;

[0109] The signal processing model adjusts the model parameters to optimize the signal processing flow based on the output of the signal decision agent and the decision of the signal path optimizer. The parameter update of the model depends on the change of the probability distribution, which reflects the change of the signal characteristics under specific decisions and states. The model updates the parameters by calculating the expected improvement of the current decision:

[0110]

[0111] Among them, Δθ SP represents the gradient of the signal processing model parameters, η represents the learning rate, It represents the expected utility of the system for the next state, which is used to evaluate the long-term impact of the current decision and guide the optimization of signal processing strategies.

[0112] In addition, in order to further improve the accuracy and efficiency of signal processing, the present invention uses the following update rules to refine signal features:

[0113]

[0114] Such a parameter updating strategy ensures that the model moves towards improving the overall performance at each time step. Through this iterative optimization process, the signal processing framework of the present invention can adapt to the ever-changing signal characteristics and improve the processing quality of lightning signals.

[0115] signal decomposition

[0116] In order to improve the accuracy and efficiency of pulse matching, the present invention proposes an ensemble empirical mode decomposition (EEMD) method to process pulse data. The traditional EMD method may affect the decomposition quality due to modal aliasing. During the decomposition process, oscillation modes of different scales appear in the same IMF, causing the IMF to lose its physical meaning. The present invention can alleviate the modal aliasing problem by adding noise terms. Adding appropriate noise can disturb the signal so that oscillation modes of different scales can be better separated during the decomposition process. The presence of noise can change the distribution of local extreme points of the signal, making the construction of the envelope more reasonable, thereby suppressing the occurrence of modal aliasing. By adding Gaussian white noise with a mean of zero to the signal and adjusting the noise intensity according to the characteristics of the signal, the quality of EMD decomposition can be effectively improved.

[0117] The steps are as follows:

[0118] EEMD is used to process pulse data. The specific steps are as follows:

[0119] Add a signal with a mean of zero and a variance of σ to the original signal x(t) 2 Gaussian white noise n i(t), and get the new signal x i (t) = x(t) + n i (t); for x i (t) Perform EMD decomposition to obtain a series of IMFs and residual terms; repeat the addition of noise and EMD decomposition M times, adding different noise each time; average the IMFs obtained from the M decompositions to obtain the final EEMD decomposition result.

[0120] Through EEMD decomposition, a more accurate and stable IMF can be obtained, laying the foundation for subsequent feature extraction and pulse matching.

[0121] Feature extraction

[0122] In order to evaluate the importance of IMF, the present invention adopts a criterion based on energy and correlation. First, the proportion of the energy of each IMF to the total energy is calculated. The IMF with larger energy usually contains the main features of the signal. Secondly, the correlation coefficient between each IMF and the original signal is calculated. The IMF with high correlation is more similar to the original signal and contains more useful information. Taking into account the two indicators of energy and correlation, the IMFs whose energy proportion and correlation are in the top K are selected as feature IMFs, where K is a preset parameter determined according to actual needs and experience. A larger K value means retaining more IMFs and containing more comprehensive signal information; a smaller K value means selecting only the most important few IMFs to highlight the key features. The present invention proposes an adaptive method to determine the K value. This method comprehensively considers the energy distribution of the IMF and the correlation with the original signal. Specifically, the present invention defines the K value as follows:

[0123]

[0124] E ratio It represents the ratio of the cumulative energy of the first k IMFs to the total energy, E j represents the energy of the jth IMF, E total represents the total energy of all IMFs.

[0125] R min =min{R 1 ,R 2 ,…,R K} (13)

[0126] R min represents the minimum correlation coefficient between the first k IMFs and the original signal, R 1 ,R 2 ,…,R K Represents the correlation coefficient between each IMF and the original signal.

[0127] C(k)=(E ratio ≥Eth )∧(R min ≥R th ) (14)

[0128] C(k) is a Boolean function. When the energy ratio is not less than the threshold E th And the minimum correlation is not less than the threshold R th True when . The minimum k value that satisfies the condition is the optimal K value. This method aims to find the minimum K value so that the selected IMF contains the main energy of the signal and maintains a high correlation with the original signal.

[0129] Pulse matching

[0130] In the pulse matching stage, the present invention compares the characteristic IMFs obtained at different stations to find the matching pair with the smallest time delay. In order to evaluate the similarity between different IMFs, the present invention uses the cross-correlation function (CCF) for calculation. The cross-correlation function can measure the correlation between two signals, that is, the degree to which one signal lags behind or leads the other signal. Unlike the Euclidean distance, the CCF takes into account the overall shape and change trend of the signal and is robust to time offset.

[0131] The cross-correlation function between two stations i and j is defined as:

[0132] CCF ij (τ)=∑[imf i (t)·imf j (t+τ)] (15)

[0133] Among them, CCF ij (τ) represents the cross-correlation function between two stations i and j, τ represents the time offset, imf i (t), imf j (t+τ) represents the characteristic intrinsic mode function of the i-th station; CCF ij The maximum value of (τ) can determine the best matching position between the two IMFs, thus achieving accurate time delay estimation:

[0134] σ tij =arg max {CCF ij (τ)} (16)

[0135] σ tij represents the time delay between station i and station j.

[0136] The present invention calculates the cross-correlation function for the characteristic IMF of all station pairs and obtains a time delay matrix Δt. By analyzing the structure of the Δt matrix, the matching relationship and time synchronization between different stations can be found. Ideally, the Δt matrix should be an antisymmetric matrix with zero main diagonal elements. Elements that deviate from this characteristic may indicate abnormal conditions, such as station failures or matching pairs with large errors.

[0137] Initial value optimization calculation

[0138] The present invention introduces a particle swarm optimization (PSO) algorithm to overcome the local convergence problem of the LM algorithm and improve the global search capability of lightning location. PSO achieves efficient solution space search by simulating the behavior of biological groups.

[0139] The present invention adopts a dual-particle swarm collaborative optimization strategy based on Tent mapping distribution. Tent mapping can make particles present uniform distribution and good correlation, forming a pseudo-random or uniform distribution. Two particle swarms are set up: Pop1 is responsible for global search and adopts a conventional PSO algorithm; Pop2 is responsible for local search and adopts a multi-scale mutation particle scheme. The particles in Pop2 are divided into 5 groups according to the initial fitness. Each group uses a Gaussian random number to obtain a mutation operator to change the particle speed. Different mutation operators and multiple groups of populations give Pop2 particles diversity to avoid falling into local optimality. The optimal result of Pop2 competes with the optimal result of Pop1 for learning to ensure the reliability of the global optimal value. Figure 3 The following is the iterative update strategy of the double group. Figure 4 As shown in the figure, when the particles converge to the local optimum and the speed is lower than the set threshold, a speed "escape" mechanism based on the multi-scale mutation operator is introduced to make the particles jump out of the local optimum and continue the global search.

[0140] After initialization, the particles in Pop2 are sorted from small to large according to fitness; then they are combined to generate M subgroups, and the number of particles in each subgroup is:

[0141]

[0142] N is the number of particles in the Pop2 population, and the fitness value of each sub-population is calculated as follows:

[0143]

[0144] Where t is the number of iterations. Different mutation scales compete with each other, and different mutation capabilities are set according to different adaptability. Therefore, the standard deviation of the t-th mutation operator is:

[0145]

[0146] Pop1 and Pop2 populations update their velocities according to the velocity update formula.

[0147] Population Pop1 is updated according to the following formula:

[0148]

[0149] Population Pop2 is updated according to the following formula:

[0150]

[0151] in, It is the position vector of the dth dimension of the optimal particle in the population Pop2. When the speed in Pop2 is less than the set threshold, speed escape will be performed.

[0152] Perform Gaussian mutation on the m subgroup particles in Pop2. The velocity vector of the particles in the mth subgroup is updated as follows:

[0153]

[0154] That is, the speed at the iteration time t plus the Gaussian distribution random number of the mutation variance. In this way, through the mutation operator, Pop2 will calculate the optimal value of the M group population in each iteration, and compete with the optimal value of the Pop1 population for learning. The better optimal value is the optimal value of the entire group of Pop1 and Pop2.

[0155] In order to compare the optimization accuracy of the LM algorithm and the DPSO algorithm for the initial solution, the present invention conducted a simulation experiment. Assuming that the ground lightning occurs at the rocket launch site, the present invention calculates the time when each station receives the lightning signal, and on this basis superimposes a random timing error with a mean of 0, a root mean square Δt of 0.2-2μs and a Gaussian distribution. The present invention uses the CHAN algorithm to calculate the initial value, and then uses the LM and DPSO algorithms to optimize the initial solution respectively, and repeats the calculation 1000 times to take the average value.

[0156] like Figure 5 As shown in the figure, the results show that when the timing error is within 1.6μs, both the LM and DPSO algorithms have good fitting accuracy. The horizontal error of the LM algorithm is between 28-211m and the vertical error is between 126-355m; the horizontal error of the DPSO algorithm is between 35-301m and the vertical error is between 97-362m. However, when the timing error exceeds 1.6μs, the calculation error of the LM algorithm increases sharply, the horizontal error exceeds 600m, and even fails to converge. In contrast, the DPSO algorithm can still maintain a high accuracy even under higher timing errors. When the time error of 1.6μs is superimposed, the horizontal error of the DPSO algorithm is only 250-417m and the vertical error is 334-445m.

[0157] Figure 6 The lightning signals measured by the BEIS station at some stages of a lightning process are given. This lightning was observed synchronously by 6 stations, and almost all of them were affected by noise interference. When the noise and lightning signals are superimposed, an envelope similar to an amplitude modulation (AM) wave is formed, which is manifested as the high-frequency lightning signal "riding" on the background noise interference, showing a periodic envelope fluctuation. When looking for valid pulses, it is usually necessary to raise the threshold, which is an effective method, but it will lose small pulses. Since the energy of the lightning electric field signal in the low-frequency band is mainly concentrated in the range of 3kHz to 30kHz, a band filter is usually used to effectively solve the noise interference problem. However, for wide-range sensors, this method is bound to miss signals outside the frequency band, which are essential for the study of lightning physics, such as high-frequency information from the leader head breakdown process and low-frequency signals from processes such as the M component.

[0158] Figure 7 , 8 a uses the positioning method of Ma et al. (2021), namely a. frequency band filtering; b. matching the pulse arrangement and combination of each station within the set window time; c. determining the optimal combination through the chi-square value; d. CHAN algorithm obtains the initial solution, and PSO or LM performs optimization fitting. This method has also been recognized by various teams and can obtain relatively accurate results in most cases. However, in this case, the signals of the synchronous stations are all interfered to varying degrees, resulting in low channel refinement and inability to fully display the branch channels. Many scattered radiation sources can be seen near the channel, which may be a normal physical phenomenon or may be due to pulse matching errors, resulting in the positioning results not truly reflecting the spatiotemporal characteristics of the radiation source. The present invention tends to believe that this is due to pulse matching errors. First, the lightning appeared at the beginning of the thunderstorm. There was no lightning in the area within the first few tens of minutes, indicating that the environmental charge accumulation was insufficient in the early stage and lightning activity had not yet occurred frequently, so there was a high possibility of a large number of false matches. Secondly, through artificial tracing, it was found that the matched pulses were not truly synchronized. When matching a specific pulse detected at a certain station, a problem is often encountered: whether it is matched to other lightning discharge events with similar time or to noise signals, as long as the time of the erroneous pulse is very close to the true observation value, it is very common for the positioning result to pass the test of a small chi-square value. In this case, although the matched pulses do not come from the same lightning process, they can still pass the chi-square test due to the coincidence in time, resulting in deviations in the positioning results. This result leads to poor initial value accuracy, and it is difficult for the LM algorithm to fit a better positioning accuracy.

[0159] Figure 7 ,8 In b, the present invention replaces the LM algorithm in the previous algorithm with the D-PSO algorithm, and performs fitting optimization on the initial solution. This improvement aims to compare whether the accuracy of the positioning calculation solution is improved when the pulse combination is consistent. As can be seen from the figure, after adopting the D-PSO algorithm, the scattered points around the radiation source are significantly reduced, which indicates that some positioning points with poor accuracy have been corrected. However, there is no particularly significant improvement in channel branching and precision, and the overall number of positioning points has not increased significantly. This result shows that simply replacing the optimization algorithm cannot fundamentally solve the problem of lightning positioning accuracy in complex environments.

[0160] Figure 7 , 8 c shows the effect of the lightning location algorithm proposed in the present invention. The algorithm first uses a graph neural network (GNN) to preprocess the signals collected by multiple stations, and then uses ensemble empirical mode decomposition (EEMD) to extract the intrinsic mode functions (IMFs) of different stations for cross-correlation analysis. By calculating the maximum value of the cross-correlation function, the best matching position between the IMFs can be determined to achieve accurate time delay estimation. On this basis, the present invention uses the CHAN algorithm to obtain the initial solution, and uses the differential particle swarm optimization (D-PSO) algorithm for further optimization, and finally obtains the time difference location result of the lightning event.

[0161] from Figure 7 , 8 It can be clearly seen from Figure c that the algorithm proposed in the present invention has achieved significant improvements in the number and precision of lightning location points. Compared with the previous algorithm, the number of location points obtained by the algorithm of the present invention has increased significantly, and the number of newly added points is about 1 / 5 of the original number of points, reaching 2103. This result is of great significance. It shows that in the previous algorithm, due to the influence of various complex factors, such as noise interference, pulse distortion, small pulse amplitude, etc., many pulse combinations that could have been matched were mistakenly ignored or discarded. The algorithm of the present invention effectively overcomes these interference factors through GNN preprocessing and EEMD-CCF, so that these pulse combinations are correctly matched and identified again. It is precisely because the number of correct pulse matching groups has increased significantly that the number of lightning radiation source points has also increased significantly. This is not only reflected in the increase in the number of location points, but more importantly, the fine structure of the lightning channel has been more clearly displayed. As can be seen from the figure, the trunk and branch channels of the lightning channel are clearly displayed, and the richness of details exceeds the previous algorithm results. This result fully demonstrates the superior performance of the algorithm of the present invention in lightning location in complex environments. It can effectively extract and restore the fine structure of the lightning channel, providing valuable data support for in-depth research on the physical process of lightning.

[0162] The examples described in the present invention are merely descriptions of the preferred implementation modes of the present invention, and are not intended to limit the concept and scope of the present invention. Without departing from the design concept of the present invention, various modifications and improvements made to the technical solutions of the present invention by engineers and technicians in this field should all fall within the protection scope of the present invention.

Claims

1. A lightning low-frequency three-dimensional positioning method for complex environments, characterized in that: The steps include: A signal processing framework is constructed to pre-process lightning signals, wherein the signal processing framework includes: a graph learning module, a signal path optimizer, a signal decision agent module and a parameter update module; Pulse matching is performed on the preprocessed lightning signal to obtain a preliminary location; The final positioning result is obtained by fitting the preliminary positioning through improving the particle swarm optimization (PSO) algorithm.

2. The method according to claim 1, characterized in that: The graph learning module is used to capture the spatial characteristics of lightning signals and output comprehensive site, signal and time embeddings. Given the embedding representation of the original signal, the embedded representation of the signal is iteratively updated. The specific formula is: in, and They represent the embedding of the current site u, the embedding of the event et at the current time t, and the embedding of the predicted event e't+1 at the next time step t+1, respectively. These embeddings are obtained through the graph learning module and reflect the site characteristics, current event characteristics, and predicted event characteristics, respectively. The formula uses these embeddings to calculate the transition probability from the current site to the next site, which is used to optimize the signal propagation path. A u , and They represent the adjacency matrix related to the site u, the current event et and the predicted event e't at the lth level, respectively. and Represent the weight matrices related to site u, current event et and predicted event e't respectively, and They represent the corresponding bias terms respectively, and σ is the nonlinear activation function; The signal path optimizer selects the optimal signal path. The specific formula is: Among them, s t =(u,e t ) represents the current state, including site u and event e at time step t t , W P represents the path selection weight matrix, and Respectively represent the event e t and e t+1 Embedding representation in graph learning module; The signal path optimizer generates an action t , the specific formula is: to t =(and t →and t+1 ) (5) Set reward r(s t ,a t ) is used to evaluate the quality of an action. The specific formula is: Among them, W r represents the reward weight matrix, N represents the number of samples in a 100 microsecond window, and h i represents the feature representation of the i-th sample, h signal represents the embedded representation of the signal in the current state, represents the expected signal feature embedding; The signal decision agent module sets the action value function Q to evaluate action a t The expected utility of is: Among them, θ Q Represents the set of parameters of the action value function, which is used to learn and predict the long-term rewards of different state-action pairs. By adjusting these parameters, the system can optimize its decision-making strategy to obtain better signal processing effects; γ represents the discount factor; In the parameter update module, the policy network parameters Δ E and signal processing model parameters θ SP Update, the specific formula is: in, represents the policy network parameters at time t, Represents the E The gradient of Representative strategy In a given state s t Next select action a t The logarithmic probability of , α represents the learning rate; Δθ SP represents the gradient of the signal processing model parameters, η represents the learning rate, Represents the expected utility of the system for the next state, which is used to evaluate the long-term impact of the current decision and guide the optimization of signal processing strategies; The signal processing model parameter update rule is refined, and the specific formula is: in, represents the signal processing model parameters at time t.

3. The method according to claim 1, characterized in that Pulse matching is performed on the preprocessed lightning signal, including: EEMD is used to process pulse data. The specific steps are as follows: Add a signal with a mean of zero and a variance of σ to the original signal x(t) 2 Gaussian white noise n i (t), and get the new signal x i (t) = x(t) + n i (t); for x i (t) Perform EMD decomposition to obtain a series of IMFs and residual terms; repeat the addition of noise and EMD decomposition M times, adding different noise each time; average the IMFs obtained from the M decompositions to obtain the final EEMD decomposition result; Feature extraction is performed on the EEMD decomposition results. The specific steps are as follows: The energy proportion and correlation of IMF are calculated respectively. The specific formula is: R min =min{R1,R2,…,R K } (13) Among them, E ratio It represents the ratio of the cumulative energy of the first k IMFs to the total energy, E j represents the energy of the jth IMF, E total represents the total energy of all IMFs; R min represents the minimum correlation coefficient between the first k IMFs and the original signal, R1, R2, …, R K Represents the correlation coefficient between each IMF and the original signal; Define the K value to indicate that the IMFs with the top K energy proportions and correlations are considered as characteristic IMFs. The specific formula is: C(k)=(E ratio ≥E th )∧(R min ≥R th ) (14) Among them, C(k) is a Boolean function. When the energy ratio is not less than the threshold E th And the minimum correlation is not less than the threshold R th True when ; the minimum k value that meets the conditions is the optimal K value; The cross-correlation function is used to evaluate the similarity between different IMFs. The specific formula is: CCF ij (τ)=∑[imf i (t)·imf j (t+τ)] (15) Among them, CCF ij (τ) represents the cross-correlation function between two stations i and j, τ represents the time offset, imf i (t), imf j (t+τ) represents the characteristic intrinsic mode function of the i-th station; By calculating CCF ij The maximum value of (τ) determines the best matching position between the two IMFs, thereby achieving accurate time delay estimation. The specific formula is: s tij =arg max {CCF ij (t)} (16) Among them, σ tij represents the time delay between station i and station j.

4. The method according to claim 1, characterized in that: The improved particle swarm optimization PSO algorithm specifically adopts a dual particle swarm collaborative optimization strategy based on Tent mapping distribution, which specifically includes the following steps: After initialization, the particles in Pop2 are sorted from small to large according to fitness; Combine them to generate M subgroups, and the number of particles P in each subgroup is: Where N is the number of particles in Pop2 population; The fitness value calculation formula for each subgroup is: in, represents the fitness of the mth subgroup, t represents the number of iterations, represents the fitness of the i-th particle in the m-th subgroup; Different mutation capabilities are set according to different adaptability. The standard deviation of the t-th mutation operator is: in, represents the standard deviation of the t-th mutation operator, FitX max and FitX min Represent the maximum and minimum fitness values ​​respectively; Pop1 population updates its speed according to the speed update formula. The specific formula is: Pop2 population updates its speed according to the speed update formula. The specific formula is: Where c1 and c2 represent acceleration constants, ω represents inertia weight, Represents the position vector of the d-th dimension of the optimal particle in the population Pop2. When the speed in Pop2 is less than the set threshold, speed escape will be performed; Gaussian mutation is performed on the m subgroup particles in Pop2, and the velocity vector of the particles in the mth subgroup is updated as follows: Here, randn represents a random number uniformly distributed in the range [0,1].

5. A storage medium, characterized in that: The storage medium stores instructions, and when a computer reads the instructions, the computer executes the lightning low-frequency three-dimensional positioning method for complex environments as described in any one of claims 1 to 4.

6. An electronic device, characterized in that: The invention comprises a processor and the storage medium as claimed in claim 5, wherein the processor executes instructions in the storage medium.

Images