A lightning multi-frequency signal data fusion positioning method

Abstract

The present application discloses a lightning radiation source positioning method based on multi-station synchronous observation, and belongs to the technical field of atmospheric science detection. The method comprises the following steps: collecting lightning broadband electromagnetic signals by multiple stations synchronously and performing time-frequency analysis; intelligently selecting high-quality frequency points to suppress noise and interference by calculating cross-station coherence, interference level and other indicators; reconstructing a narrowband signal for each high-quality frequency point and performing high-precision time difference (TDOA) estimation by using an improved generalized cross-correlation method; obtaining the three-dimensional coordinates of lightning by weighting and fusing the TDOA measurement values of multiple high-quality frequency points and optimizing the solution by using a Grid Search algorithm. The present application intelligently selects high-quality frequency points, accurately estimates the time difference and optimizes the positioning algorithm, and effectively overcomes the shortcomings of traditional single-frequency positioning methods, such as being easily disturbed and limited precision, by multi-band information fusion and intelligent processing, thereby significantly improving the precision, anti-interference ability and reliability of lightning positioning, and being suitable for lightning positioning requirements in various environments.

Description

Technical Field

[0001] This invention relates to the field of atmospheric science detection technology, and in particular to a lightning radiation source location method based on multi-station synchronous observation. Specifically, it relates to a lightning location method that improves location accuracy by intelligently selecting high-quality frequency points and fusing multi-frequency band data. Background Technology

[0002] Lightning is a powerful atmospheric discharge phenomenon that poses a serious threat to aerospace, power transmission, forest fire prevention, and the safety of people's lives and property. Accurate detection and location of the location, intensity, and development process of lightning are crucial for lightning early warning and scientific research. The Time Difference of Arrival (TDOA) method is a widely used technique for lightning location. Its basic principle is to measure the time difference between the arrival of the electromagnetic signal radiated by lightning at multiple different monitoring stations, construct a hyperbolic equation system, and ultimately solve for the spatial location of the lightning radiation source. However, existing single-band TDOA location methods face many challenges in practical applications, resulting in limited location accuracy and reliability: 1) Limited frequency selection and interference sensitivity: Lightning radiation signals have a wide spectrum, but traditional methods often use a fixed frequency band (such as low frequency) for location. The actual electromagnetic environment is complex; this frequency band may be subject to strong communication broadcasting or industrial noise interference, leading to signal distortion, increased TDOA measurement errors, and even location failure. 2) Gross errors in measurements: Due to multipath effects and noise interference during electromagnetic wave propagation, the TDOA values ​​measured at certain stations or at certain times may be outliers (gross errors) that significantly deviate from the true values. Traditional least squares algorithms are very sensitive to gross errors, and even a few outliers can cause the positioning results to deviate significantly from the true location. 3) Model errors and solution stability: The TDOA positioning equation is nonlinear, and the solution process usually involves linearization approximation. When the station geometry is poor or the initial values ​​are not chosen properly, it is easy to get trapped in local optima or cause the solution to fail. Although the Chan algorithm can provide a closed-form solution, its accuracy decreases when the error is large; the grid optimization method has high accuracy but a large computational load. 4) Weak signal processing: Traditional cross-correlation algorithms directly process the original signal without fully considering noise and baseline drift in the signal, which affects the accuracy of time difference estimation.

[0003] To address these issues, researchers have explored various improvement methods, such as using multiple frequency bands for independent positioning followed by averaging, or employing more complex signal processing algorithms. However, these methods either fail to fundamentally achieve the organic fusion of inter-frequency information or suffer from low computational efficiency, making it difficult to meet the demands of real-time, precise positioning. Therefore, there is an urgent need for a high-precision lightning positioning method that can intelligently adapt to complex electromagnetic environments, effectively fuse multi-frequency band information, and possess strong anti-interference capabilities. Summary of the Invention

[0004] The purpose of this invention is to provide a lightning multi-frequency signal data fusion and localization method, comprising: synchronously acquiring broadband electromagnetic signals of lightning from multiple stations and performing time-frequency analysis; intelligently selecting high-quality frequency points to suppress noise and interference by calculating indicators such as cross-station coherence and interference level; reconstructing narrowband signals for each high-quality frequency point and using an improved generalized cross-correlation method to perform high-precision time difference of origin (TDOA) estimation; and finally obtaining the three-dimensional coordinates of the lightning by weighted fusion of TDOA measurements from multiple high-quality frequency points and optimizing the solution using the Grid Search algorithm. This invention, through multi-band information fusion and intelligent processing, effectively overcomes the shortcomings of traditional single-band localization methods, such as susceptibility to interference and limited accuracy, and significantly improves the accuracy, anti-interference capability, and reliability of lightning localization.

[0005] This invention is achieved through the following technical solution: a lightning multi-frequency signal data fusion positioning method, comprising the following steps:

[0006] S1: Broadband Signal Acquisition and Time-Frequency Analysis: Broadband electromagnetic signals from lightning radiation are simultaneously acquired at multiple stations. Short-Time Fourier Transform (STFT) is performed on the signals from each station to obtain the time-frequency matrix X. i (f,t).

[0007] S2: Intelligent Selection of High-Quality Frequency Points: Based on the time-frequency matrix, the cross-site coherence index and interference level index of each frequency point are calculated, and a comprehensive frequency point scoring function is constructed.

[0008] ,

[0009] Where M is the total number of features considered simultaneously in this multivariate analysis; λ represents the squared amplitude coherence between station i and station j; max (f) is the largest eigenvalue of the cross-spectral matrix C(f) at frequency point f; λ k (f) is the k-th eigenvalue of a matrix associated with feature f; P int (f) represents the estimated interference power at frequency point f; ε is a small normal number. The K frequency points with the highest overall scores are selected to form a high-quality frequency point set F. sel .

[0010] S3, Multi-frequency TDOA Measurement: For each high-quality frequency point f k Perform the following operations:

[0011] S31. Signal Extraction: Extract the narrowband signal near the frequency point from the Short Time Fourier Transform (STFT) result;

[0012] S32. Time Difference Estimation: The time difference of signal arrival at each station is calculated using the generalized cross-correlation method.

[0013] ,

[0014] in, R is the time difference between the arrival of the signal at each station. ij (τ) is the cross-correlation function between stations i and j;

[0015] S33. Optimal TDOA Calculation: For each pair of stations (i,j), calculate an optimal TDOA estimate by integrating information from all high-quality frequency points. Using the quality score of each high-quality frequency point as a weight, calculate the weighted average as the optimal TDOA estimate for that station pair. The mathematical expression is as follows:

[0016] .

[0017] S4. Multi-frequency TDOA fusion positioning: The TDOA measurement values ​​of each high-quality frequency point are weighted and fused, and the lightning location r is solved through an optimization algorithm.

[0018] ,

[0019] Where c is the speed of light, r i , r j Let i be the position of station i and station j.

[0020] S5: Positioning Result Optimization: The Grid Search algorithm is used to optimize the preliminary positioning results and output the final three-dimensional coordinates.

[0021] Furthermore, in step S2, the selection of the high-quality frequency point set adopts an adaptive threshold determination method:

[0022] Calculate the statistical characteristics of the comprehensive score for all frequency points, including the mean μ and standard deviation σ;

[0023] Dynamically set the selection threshold: Threshold = μ + α·σ, where α is an adjustable parameter;

[0024] Statistical features are recalculated at specific time windows to achieve adaptive updates of the threshold.

[0025] Furthermore, the generalized cross-correlation method described in step S32 employs the PHAT algorithm for weighting:

[0026] ,

[0027] Among them, X i (f) represents the signal spectrum of station i; X *j (f) represents the complex conjugate of the Fourier transform of signal j; τ represents the time delay; It is a complex exponential function, where h is the imaginary unit.

[0028] Furthermore, a two-stage localization strategy is adopted in step S4:

[0029] Phase 1: Use the Chan algorithm to solve the closed-form solution and obtain the initial position estimate;

[0030] The second stage: using the Chan algorithm result as the initial value, the Grid Search method is used for iterative optimization.

[0031] Further, in step S31, the extracted narrowband signal undergoes envelope alignment preprocessing:

[0032] First, the envelope of the signal at each station is extracted by Hilbert transform to obtain the time-domain waveform that reflects the change in signal energy;

[0033] Next, the cross-correlation method is used to calculate the signal envelope of each station, and the rough time difference is estimated based on a large time window to achieve preliminary signal alignment.

[0034] Finally, the original signal within a small time window after initial alignment is extracted, and high-precision cross-correlation analysis is performed using the GCC-PHAT method. The time delay estimation is further refined using phase information to achieve accurate alignment at the sub-sampling interval level.

[0035] Furthermore, in step S4, a statistical identification method based on the mean and standard deviation is used to detect and remove outliers in the TDOA measurements:

[0036] First, calculate the mean and standard deviation of all TDOA measurements;

[0037] Measurements that deviate from the mean by more than a preset multiple of the standard deviation are identified as outliers and removed.

[0038] The mean is recalculated using the remaining valid measurements, and the above process is iterated until no new outliers appear.

[0039] Finally, lightning location calculation is performed based on the converged set of valid TDOA measurements.

[0040] Furthermore, in the weighted fusion process of step S4:

[0041] Assign a confidence weight to the TDOA measurement value for each high-quality frequency point;

[0042] The credibility weight is determined by the overall score, signal-to-noise ratio, and cross-site consistency of the frequency point.

[0043] The beneficial effects of this invention are as follows: By intelligently selecting high-quality frequency points, the invention ensures that the frequency points participating in positioning contribute the most to the positioning, thus improving positioning accuracy. The use of generalized cross-correlation and the PHAT algorithm for weighted calculation effectively suppresses the influence of noise and interference, improving the accuracy of time difference estimation. The adoption of a two-stage positioning strategy and the Grid Search algorithm improves the accuracy and reliability of the positioning results. The adaptive threshold update method dynamically adjusts the frequency point selection threshold according to environmental changes, enhancing the method's adaptability. Envelope alignment preprocessing and high-precision cross-correlation analysis further improve the signal alignment accuracy, providing a foundation for accurate time difference estimation. Attached Figure Description

[0044] Figure 1 This is a flowchart of the positioning method of the present invention;

[0045] Figure 2 A schematic diagram of the intelligent selection process for high-quality frequency points;

[0046] Figure 3 This is a schematic diagram of the multi-frequency TDOA measurement and fusion positioning process;

[0047] Figure 4 A schematic diagram illustrating the optimization process of the Grid Search algorithm. Detailed Implementation

[0048] To enable those skilled in the art to better understand the present invention, in conjunction with Figures 1-4 Further explanation of this application is provided, and the content mentioned in the embodiments is not intended to limit the invention.

[0049] Implementation environment: Deploy at least 5 lightning signal stations with broadband reception capability and high-precision GNSS clock synchronization, with the distance between the stations ranging from several kilometers to hundreds of kilometers.

[0050] The lightning multi-frequency signal data fusion localization method of this application is illustrated in the flowchart below. Figure 1 This includes the following steps:

[0051] S1. Broadband Signal Acquisition and Time-Frequency Analysis: Each station uses a high-speed ADC acquisition card to synchronously acquire broadband electromagnetic signals from lightning radiation. Multiple stations perform Short-Time Fourier Transform (STFT) on the signals from each station. Through STFT processing, the one-dimensional time-domain signal is converted into a two-dimensional time-frequency matrix X. i (f,t) is used to obtain the energy distribution of the signal at different times and frequencies, providing a data foundation for subsequent selection of high-quality frequency points and signal processing.

[0052] S2. Intelligent Selection of High-Quality Frequency Points: Based on indicators such as cross-station coherence and interference level, a comprehensive scoring function is constructed to adaptively select a set of high-quality frequency points with high signal-to-noise ratio and strong anti-interference capability from a wide spectrum, ensuring that subsequent processing is based on the most reliable data and providing a reliable data foundation for subsequent processing. The selection of the high-quality frequency point set adopts an adaptive threshold determination method, calculating the statistical characteristics of the comprehensive score of all frequency points, including the mean μ and standard deviation σ; dynamically setting the selection threshold: Threshold = μ + α·σ, where α is an adjustable parameter; recalculating the statistical characteristics at specific time windows to achieve adaptive updating of the threshold. Based on the time-frequency matrix, the cross-station coherence index and interference level index of each frequency point are calculated, and the maximum eigenvalue λ of the cross-spectral matrix C(f) composed of signals from all stations is calculated. max (f). Interference power

[0053] P int The calculation of (f) takes a neighborhood centered at frequency point f, calculates the minimum power value within this neighborhood after excluding possible main lobe signal regions, and constructs a frequency point comprehensive scoring function:

[0054] ,

[0055] Where M is the total number of features (or signal channels, number of sensors) considered simultaneously in this multivariate analysis. It is the dimensional basis of the entire formula and determines the number of feature pairs for coherence calculation and the number of feature values ​​obtained after feature decomposition. λ represents the squared amplitude coherence between station i and station j; max (f) is the largest eigenvalue of the cross-spectral matrix C(f) at frequency point f; λ k (f) is the k-th eigenvalue of a matrix associated with feature f, which quantifies the strength of a certain pattern contained in the matrix; P int (f) represents the estimated interference power at frequency point f; ε is a small normal number (e.g., 0.001 or 0.00001) introduced to ensure the mathematical and numerical stability of the formula. Its main function is to avoid the interference power P from being too high. int (f) indicates situations where the value is divided by zero or becomes unstable when the value is too small. The K frequency points with the highest overall scores are selected to form a high-quality frequency point set F. sel .

[0056] A schematic diagram of the intelligent selection process for high-quality frequencies is shown below. Figure 2 A detailed explanation of the intelligent selection process for high-quality frequencies, which includes seven core steps:

[0057] Step 1: Input the time-frequency matrix. The starting point of the process is the time-frequency matrix Xi(f, t) obtained by multiple stations through short-time Fourier transform (STFT). This matrix fully represents the energy distribution of the signal at different frequencies f and times t, and is the data basis for all subsequent analyses.

[0058] Step 2: Calculate the cross-station coherence index, used to assess the waveform similarity or correlation of the same frequency signal across different stations. Higher coherence indicates a greater likelihood that the signal originates from the same source (lightning), less susceptibility to local noise interference, and higher reliability for timing and location. Two key indices are calculated: 1) the maximum eigenvalue λ of the cross-spectral matrix. max (f), for each frequency point f, calculate the cross-spectral matrix formed by the signals from all stations, and find its largest eigenvalue λ. max The magnitude of (f) reflects the spatial coherence of the signal at that frequency point; a larger value indicates better coherence. 2) Amplitude Squared Coherence MSC ij (f) Calculate the amplitude squared coherence of every two stations i and j at the same frequency point f, and then average it for all station pairs as a comprehensive measure of coherence at that frequency point.

[0059] Step 3: Calculate the interference level index, which quantifies the severity of background noise and man-made interference at each frequency point. Calculate the interference power estimate P. int (f). The usual practice is to calculate the minimum power spectrum of a small neighborhood centered at frequency point f (excluding the part that may contain the main signal), and use this minimum value as an estimate of the interference level of that frequency point. The smaller the value, the "cleaner" the frequency point is.

[0060] Step 4: Construct a frequency point comprehensive scoring function to integrate the multiple indicators calculated in Steps 2 and 3 into a unified, quantifiable comprehensive score for sorting and filtering all frequency points. This function aims to find frequency points with both "high coherence" and "low interference level." The coherence index is used as the numerator, and the interference level is used as the denominator, ensuring that frequency points that simultaneously meet both conditions receive the highest comprehensive score.

[0061] Step 5: Determine the dynamic threshold. Set an automated screening criterion, rather than a fixed threshold, to allow the system to adapt to changing electromagnetic environments at different times and locations. Calculate the dynamic threshold based on the statistical characteristics (such as the global mean μ and standard deviation σ) of the comprehensive scores of all frequency points. Recalculate the statistical characteristics and update the threshold at specific time windows. A common method is: Threshold = μ + α·σ, where α is an adjustable parameter (e.g., set to 1.5). Frequency points with scores higher than this threshold are considered "high-quality".

[0062] Step Six: Select high-quality frequency points. Choose the top K frequency points with a comprehensive score higher than the threshold, sort them from highest to lowest, and compare the comprehensive score of each frequency point with the dynamic threshold. Select all frequency points with a comprehensive score higher than the dynamic threshold. Sometimes, to further control the number, the top K frequency points with the highest scores will be selected again to form the final set of high-quality frequency points F. sel .

[0063] Step 7: Output the set of high-quality frequency points. The selected set of high-quality frequency points F... sel The output serves as input for the next process (i.e., "multi-frequency TDOA measurement"). Subsequent processing will only be performed on these selected frequencies, thereby significantly improving data quality and processing efficiency.

[0064] S3, Multi-frequency TDOA Measurement: For each high-quality frequency point f k Perform the following operations:

[0065] S31. Signal Extraction: Extract the narrowband signal near the specified frequency from the Short Time Fourier Transform (STFT) results. The extracted narrowband signal undergoes envelope alignment preprocessing: First, the envelope of each station's signal is extracted using Hilbert transform to obtain a time-domain waveform reflecting signal energy changes. Next, the cross-correlation method is used to calculate the envelope of each station's signal, estimating a coarse time difference based on a large time window to achieve initial signal alignment. Finally, the original signal within a small time window after initial alignment is truncated, and high-precision cross-correlation analysis is performed using the GCC-PHAT method. Phase information is then used to further refine the time delay estimation, achieving precise alignment at the sub-sampling interval level.

[0066] S32. Time Difference Estimation: The GCC-PHAT method is used to calculate the cross-correlation function between all station pairs. By detecting the peak position of the cross-correlation function, the TDOA measurement value τ of the signal arriving at each station is calculated. ij :

[0067] ,

[0068] in, R is the time difference between the arrival of the signal at each station. ij (τ) is the cross-correlation function between stations i and j; The symbol "^" indicates that this is an estimate; This time delay estimate is for a specific frequency component f. k Performed; This means finding the value of the independent variable τ that makes a certain function reach its maximum value; simply put, it means "finding the x-coordinate value corresponding to the peak point".

[0069] The generalized cross-correlation method uses the PHAT algorithm for weighting:

[0070] ,

[0071] Among them, X i (f) represents the signal spectrum of station i; X * j (f) is the complex conjugate of the Fourier transform of signal j; in the frequency domain, finding the complex conjugate is equivalent to performing a time reversal operation in the time domain; τ is the time delay, which represents the time delay of signal j relative to signal i when calculating coherence. By changing the value of τ, we can analyze how the correlation between the two signals changes with the time delay. It is a complex exponential function, the core of the inverse Fourier transform, where h is the imaginary unit (hi). 2 =-1), which is responsible for mapping the phase information in the frequency domain to the time delay τ. Molecular X i (f)X * j (f) is the kernel of the cross-power spectral density of the two signals in the frequency domain, which contains information about the amplitude relationship and phase difference between the two signals at different frequencies. The denominator |X i (f)X * j (f) is the modulus (or absolute value) of the molecule above. Its function is to normalize. Dividing by it eliminates the amplitude information of each frequency component in the molecule, leaving only the phase information.

[0072] S33. Optimal TDOA Calculation: Perform the calculation independently for each pair of stations (i,j), and combine the information from all high-quality frequency points to calculate an optimal TDOA estimate.

[0073] For a specific station pair (i, j), the input is a set of TDOA measurements obtained from step S32: τ ij (f1), τ ij (f2), ..., τ ij (f k ), and the corresponding quality score (f) for each frequency point. k To improve robustness, a simple statistical screening was first performed on the TDOA values ​​of this group, calculating all τ values. ij (f k The mean μ τ and standard deviation σ τ Measurements deviating from the mean by more than a preset number of standard deviations (e.g., N=2) are considered outliers and temporarily removed, forming a set of valid measurements. Using the quality score of each high-quality frequency point as a weight, a weighted average is calculated as the optimal TDOA estimate for that station pair. Its mathematical expression is:

[0074] .

[0075] S4. Multi-Frequency TDOA Fusion Positioning: A credibility weight is assigned to the TDOA measurement value of each high-quality frequency point. The credibility weight is determined by the comprehensive score, signal-to-noise ratio, and cross-site consistency of that frequency point. TDOA values ​​measured from all high-quality frequency points are combined into a measurement set. An outlier detection and removal method based on mean and standard deviation is used: First, the mean and standard deviation of all TDOA measurements are calculated; measurements deviating from the mean by more than a preset multiple of the standard deviation are identified as outliers and removed; the mean is recalculated using the remaining valid measurements, and the above process is iterated until no new outliers appear; a two-stage positioning strategy is adopted.

[0076] The first stage (coarse localization) uses the Chan algorithm to approximate the linearization of the problem and quickly calculates a closed-form solution as the initial position estimate (X0, Y0, Z0). This initial value effectively narrows the search range of the global optimal solution and avoids grid optimization from getting stuck in local minima, thus obtaining the initial position estimate.

[0077] The second stage (fine localization): Using the preliminary localization results (X0, Y0, Z0) obtained by the Chan algorithm as the center and reference, a three-dimensional search grid is established. The grid range and step size can be adaptively adjusted according to the prior error. The Grid Search algorithm is used to perform a fine iterative search in the local three-dimensional space. For each point within the grid, the weighted sum of squared residuals between the theoretical TDOA value and the actual measured TDOA value is calculated as the cost function. By continuously narrowing the search range and step size (e.g., from ±5 km, step size 500 m to ±2.5 km, step size 250 m), the optimal solution (X0, Y0, Z0) is gradually approached and the global weighted sum of squared residuals is minimized. opt Y opt Z opt Finally, based on the converged set of valid TDOA measurements, lightning location is calculated, and the lightning position r is determined using an optimized algorithm.

[0078] ,

[0079] Where c is the speed of light, r i , r j Let i be the location of station i and station j; ||r - rᵢ|| represents the distance from the lightning target source r to the i-th station, and ||r - rᵢ|| represents the distance from the lightning target source r to the i-th station. j ||This represents the distance from the lightning target source r to the j-th measuring station;||r - rᵢ||-|r - r j || represents the distance difference between the lightning strikes at stations i and j;

[0080] This step outputs the preliminary positioning result, namely the mathematically optimal three-dimensional coordinates (X). opt Y opt Z opt The result includes the corresponding residuals, which may contain jitter due to measurement noise or model error.

[0081] See the schematic diagram of the multi-frequency TDOA measurement and fusion positioning process. Figure 3 The multi-frequency TDOA measurement and fusion positioning process is a complete computational flow from raw frequency data to preliminary positioning results. Its core logic is "parallel processing and serial fusion." The entire process can be divided into three main steps: multi-frequency parallel measurement, data centralization and weighting, and hybrid algorithm positioning solution.

[0082] Step 1: Parallel measurement of TDOA at multiple frequencies. In this stage, multiple high-quality frequencies are processed in parallel and independently to generate massive amounts of raw time difference measurement data.

[0083] The process begins with the set of high-quality frequency points F. sel This set contains K high-quality frequency points (f1, f2, ..., f3) with high signal-to-noise ratio and strong anti-interference capability, which were obtained through intelligent screening in the early stage. k Each frequency point f in the set. k Each signal is processed through a separate, identically structured parallel branch. Each branch executes the following two sub-steps sequentially: 1) Narrowband signal extraction: Extract the corresponding frequency point f from the original wideband signal. k 1) For narrowband time-domain signals, specific frequency components are separated to reduce interference from other frequency bands. 2) The PHAT algorithm calculates the TDOA value. The generalized cross-correlation-phase transform (GCC-PHAT) algorithm is used to perform cross-correlation calculation on the narrowband signal. The PHAT algorithm weights and sharpens the correlation peaks, thereby estimating the time difference (TDOA) of the signal arriving at different stations with high accuracy, denoted as τ. ij (f k );

[0084] The TDOA values ​​calculated from all K parallel branches are ultimately merged into a "multi-frequency TDOA measurement result pool". The pool contains the original TDOA measurement values ​​of all frequencies and all station pairs, forming the data foundation for subsequent fusion positioning.

[0085] Step 2: Data Centralization and Weighting. This step performs serial quality control on the raw data generated in the first stage, providing high-quality, high-reliability input for localization. Selected steps: 1) Perform statistical outlier detection on all TDOA values ​​in the "Measurement Result Pool". Specifically, calculate the mean (μ) and standard deviation (σ) of all values, and identify and remove measurements that deviate from the mean by more than N times the standard deviation (e.g., N=2). The output is the "cleaned TDOA dataset". 2) Calculate credibility weights. Assign a credibility weight coefficient to each TDOA measurement value in the "cleaned TDOA dataset". This weight is calculated based on multiple indicators: Frequency Point Score (f... k ): Quality score of the frequency point to which this TDOA value belongs; the higher the score, the greater the weight. Signal-to-noise ratio (SNR): The signal-to-noise ratio of the signal used to calculate this TDOA value. Cross-site consistency: The degree of agreement between this TDOA value and other relevant measurements. The output is a "weighted TDOA dataset," with each data point assigned its weight, ready for optimization calculations.

[0086] Step 3: Hybrid algorithm for localization and solution. This stage utilizes a weighted TDOA dataset and a two-stage hybrid algorithm to solve for the lightning location:

[0087] The first stage involves using the Chan algorithm to solve a closed-form solution for coarse localization. Weighted TDOA values ​​are input into the Chan algorithm, which approximates the nonlinear equations through mathematical transformations, quickly calculating a closed-form solution as the initial estimate of the lightning location (X0, Y0, Z0). This step is computationally efficient, providing a good starting point for the next fine-grained search and avoiding the massive computational burden of a global search.

[0088] The second stage: Using the Chan algorithm result as the initial value, the Grid Search algorithm is used for iterative optimization and fine localization. A search grid is constructed in a three-dimensional space centered on the initial position (X0, Y0, Z0) given by the Chan algorithm, and the Grid Search algorithm is used for iterative optimization. The algorithm traverses each point in the grid, calculating the sum of squared residuals between the theoretical TDOA value and the weighted actual measured TDOA value for that point. By continuously narrowing the search range and step size, the point that minimizes the global weighted sum of squared residuals is iteratively found. This step is a fine correction of the Chan algorithm result, effectively approximating the optimal solution of the nonlinear optimization problem.

[0089] The preliminary positioning result is output. The optimal spatial point output by the grid optimization algorithm is the preliminary positioning result, which is the estimated three-dimensional coordinate value of the lightning.

[0090] A schematic diagram of the Grid Search algorithm optimization process can be found here. Figure 4 The entire process forms a closed-loop iterative system until the convergence condition is met. The detailed optimization process of the hybrid algorithm is as follows:

[0091] Step 1: Input and Initialization. The process begins by receiving the output results from the previous steps, namely the "weighted TDOA dataset" and the "spatial coordinates of the station." This is all the data necessary for the location calculation. Based on this input data, the lightning location coordinates that minimize the objective function (weighted sum of squared residuals) are calculated.

[0092] Step Two: The Chan algorithm quickly obtains an initial estimate. This algorithm is an analytical method that transforms the nonlinear TDOA equations into a linear model by introducing intermediate variables and two least-squares estimates, thus directly obtaining the closed-form solution. The purpose of this step is to perform fast coarse localization. It is computationally efficient, outputting an initial position estimate (X0, Y0, Z0) within milliseconds. This solution may have some bias due to measurement noise, but it provides a crucial, near-true initial point for the subsequent time-consuming fine-grained search, avoiding the blind global search of Grid Search.

[0093] Step 3: Mesh Parameter Settings. Using the initial estimate (X0, Y0, Z0) given by the Chan algorithm in the previous step as the center, construct a 3D cubic search mesh. Two key parameters need to be set here: 1) Search range: For example, offset ±5km in each of the X, Y, and Z directions. 2) Initial step size: For example, 500m, which is the interval between each point in the mesh.

[0094] Step 4: Grid Search Iteration. The algorithm traverses each grid point (X) within the defined 3D grid. i Y j Z k For each point, its cost function value is calculated (usually the weighted sum of squared residuals of all TDOA measurements). The smaller the cost function value, the closer the grid point is to the actual lightning location.

[0095] Step 5: Grid optimization. After traversing all points within the grid, find the grid point with the smallest cost function value in this round and mark its coordinates as (X). best Y best Z best This is the optimal solution for this round of search.

[0096] Step Six: Convergence Assessment. This is a decision point to determine whether the calculation results meet the accuracy requirements. Check if the current grid step size is less than or equal to the preset convergence threshold ε (e.g., ε = 10 meters). Branch One (No): If the step size is still greater than the threshold, it means the required accuracy has not yet been achieved. The process will proceed to the "narrowing the search range and reducing the step size" step. For example, using the currently found optimal point (X... best Y best Z best Using the new center, the search range is halved (e.g., from ±5km to ±2.5km), and the step size is also halved (e.g., from 500m to 250m). Then, the process returns to the "Grid Parameter Settings" step with an arrow to begin a new round of finer grid search. Branch Two (Yes): If the step size is less than or equal to the convergence threshold, the result is accurate enough, and iterative optimization can terminate.

[0097] Step 7: Output the final result. Once the convergence condition is met, the current optimal grid point (X) is... best Y best Z best This is the final solution. The process ends, and the final optimized 3D coordinates (X, Y, Z) of the lightning are output.

[0098] S5: Positioning Result Optimization: This step optimizes and constrains the preliminary mathematical results output from S4 in a physical sense, inputting the preliminary positioning result (X) from step S4. opt Y opt Z opt The system performs time-series filtering on the location results of continuous time series, using moving averages for smoothing to filter out unreasonable positional jumps, making the lightning channel path smoother and more reasonable, improving its smoothness and actual reliability, and outputs the results. Based on factors such as residuals calculated by S4, the system estimates observation noise in real time and outputs the accuracy evaluation indicators of this location (such as error ellipse and confidence level). Finally, it outputs the optimized three-dimensional coordinates of the lightning occurrence site and the accuracy evaluation indicators such as the standard deviation of the residuals of this location.

[0099] The specification and drawings of this application are merely one specific embodiment and are not restrictive. Those skilled in the art can make many other modifications based on the teachings of this application without departing from the spirit and scope of the claims, all of which are within the scope of protection of this application.

Claims

1. A method for fusion and localization of lightning multi-frequency signal data, characterized in that, Includes the following steps: S1. Broadband Signal Acquisition and Time-Frequency Analysis: Broadband electromagnetic signals from lightning radiation are simultaneously acquired at multiple stations. Short-Time Fourier Transform (STFT) is performed on the signals from each station to obtain the time-frequency matrix. X i ( f , t ); S2. Intelligent Selection of High-Quality Frequency Points: Based on the time-frequency matrix, the cross-site coherence index and interference level index of each frequency point are calculated, and a comprehensive frequency point scoring function is constructed. ， in, M This refers to the total number of features considered simultaneously in this multivariate analysis; Inter-station amplitude squared coherence; λ max ( f (Frequency point) f Cross-spectral matrix C( f The largest eigenvalue; λ k ( f ) for features f The related matrix of the first k One eigenvalue; P int ( f (Frequency point) f The estimated interference power; ε It is a tiny normal number; the one with the highest overall score is selected. K Each frequency point constitutes a set of high-quality frequency points F sel ; S3, Multi-frequency TDOA Measurement: For each high-quality frequency point f k Perform the following operations: S31. Signal Extraction: Extract the narrowband signal near the frequency point from the Short Time Fourier Transform (STFT) result; S32. Time Difference Estimation: The time difference of signal arrival at each station is calculated using the generalized cross-correlation method. ， in, For specific frequency components f k The time difference of signal arrival at each station R ij (τ) represents the station. i and j The cross-correlation function; S33, Optimal TDOA Calculation: For each pair of stations ( i,j The process is conducted independently, integrating information from all high-quality frequency points to calculate an optimal TDOA estimate. Using the quality score of each high-quality frequency point as a weight, a weighted average is calculated as the optimal TDOA estimate for that station pair. The mathematical expression for this is: ， S4. Multi-frequency TDOA fusion positioning: The TDOA measurement values ​​of each high-quality frequency point are weighted and fused, and the lightning location r is solved through an optimization algorithm. ， Where c is the speed of light, r i , r j For the station i and monitoring station j Location; S5: Positioning Result Optimization: The Grid Search algorithm is used to optimize the preliminary positioning results and output the final three-dimensional coordinates.

2. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, In step S2, the selection of the high-quality frequency point set adopts an adaptive threshold determination method: Calculate the statistical characteristics of the overall score for all frequency points, including the mean. μ and standard deviation σ ; Dynamically set the selection threshold: Threshold = μ + α · σ ,in, α These are adjustable parameters; Statistical features are recalculated at specific time windows to achieve adaptive updates of the threshold.

3. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, The generalized cross-correlation method described in step S32 uses the PHAT algorithm for weighting: ， in, X i ( f ) for measuring station i The signal spectrum; X * j ( f ) is a signal j The complex conjugate of the Fourier transform; τ For time delay; It is a complex exponential function, where h It is the imaginary unit.

4. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, Step S4 employs a two-stage localization strategy: Phase 1: Use the Chan algorithm to solve the closed-form solution and obtain the initial position estimate; The second stage: using the Chan algorithm result as the initial value, the Grid Search method is used for iterative optimization.

5. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, In step S31, the extracted narrowband signal undergoes envelope alignment preprocessing: First, the envelope of the signal at each station is extracted by Hilbert transform to obtain the time-domain waveform that reflects the change in signal energy; Next, the cross-correlation method is used to calculate the signal envelope of each station, and the rough time difference is estimated based on a large time window to achieve preliminary signal alignment. Finally, the original signal within a small time window after initial alignment is extracted, and high-precision cross-correlation analysis is performed using the GCC-PHAT method. The time delay estimation is further refined using phase information to achieve accurate alignment at the sub-sampling interval level.

6. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, In step S4, a statistical identification method based on the mean and standard deviation is used to detect and remove outliers in the TDOA measurements: First, calculate the mean and standard deviation of all TDOA measurements; Measurements that deviate from the mean by more than a preset multiple of the standard deviation are identified as outliers and removed. The mean is recalculated using the remaining valid measurements, and the above process is iterated until no new outliers appear. Finally, lightning location calculation is performed based on the converged set of valid TDOA measurements.

7. The lightning multi-frequency signal data fusion positioning method according to claim 1, characterized in that, During the weighted fusion process in step S4: Assign a confidence weight to the TDOA measurement value for each high-quality frequency point; The credibility weight is determined by the overall score, signal-to-noise ratio, and cross-site consistency of the frequency point.

Images