Infrasonic wave signal positioning analysis method and system

Abstract

The invention relates to the technical field of safety monitoring and signal processing, in particular to an infrasonic wave signal positioning analysis method and system, and the method comprises the steps: obtaining multi-channel infrasonic wave signals, extracting multipath propagation features, generating propagation compensation parameters, completing the signal correction, and obtaining an infrasonic wave signal set; calculating a time difference of arrival set according to the infrasonic wave signal set, estimating an initial direction of arrival, updating array geometric parameters based on the initial direction of arrival and the time difference of arrival set to estimate the direction of arrival, and performing weighted fusion according to consistency evaluation to obtain a fused direction; inputting the fused direction into a direction filter, updating noise parameters in combination with a historical direction-of-arrival reference sequence, and outputting a refined direction-of-arrival; a time difference of arrival residual is formed based on the refined direction of arrival and is corrected to obtain a final direction of arrival, the final direction of arrival is fused with pipeline geographic information to determine a leakage position, and the positioning stability and the high-precision positioning capability are improved.

Description

Technical Field

[0001] This invention relates to the field of safety monitoring and signal processing technology, specifically to a method and system for infrasound signal localization and analysis. Background Technology

[0002] Oil and gas pipelines are laid over long distances and across regions. Once a leak occurs, it can easily cause environmental pollution, fires, explosions, and losses due to transportation disruptions. The industry has higher requirements for early detection and rapid location, and is gradually incorporating leak monitoring into the fault prediction and health management system to achieve full life cycle risk control.

[0003] Existing leak location methods often face engineering challenges such as strong on-site noise, significant multipath propagation, limited on-site deployment density of smart sensors and power and communication conditions, and deviations in the installation location of smart sensors. These problems lead to unstable time difference of arrival, easy jitter or jump in azimuth estimation, and the location results are sensitive to single observations and geometric errors. At the same time, the fusion of azimuth results and pipeline geographic information lacks residual constraints and self-correction mechanisms, making it difficult to form a verifiable and consistent link. This affects high-precision positioning and reliable alarms in complex pipeline scenarios, and thus restricts the closed-loop implementation of fault prediction and health management. Summary of the Invention

[0004] This invention provides a method and system for infrasound signal localization analysis, which can at least solve the problem that the leakage orientation is unstable and difficult to achieve verifiable high-precision localization under the conditions of multipath noise and array geometric errors.

[0005] In a first aspect, the present invention provides a method for locating and analyzing infrasound signals, comprising the following: The multi-channel infrasound signals of oil and gas pipelines are acquired, multipath propagation characteristics are extracted and propagation compensation parameters are determined, and the multi-channel infrasound signals are corrected based on the propagation compensation parameters to obtain an infrasound signal set. The time difference of arrival (TDA) set is calculated based on the infrasound signal set. The initial azimuth of arrival (AHA) is estimated based on the TDA set. The array geometric parameters are updated and the AHA is estimated based on the initial AHA and the TDA set. The fused AHA is obtained by weighted fusion. The azimuth input azimuth filter is fused, and the azimuth filter noise parameters are updated based on the historical azimuth reference sequence to output a refined azimuth. The arrival time difference residual is calculated based on the refined wave arrival position. The refined wave arrival position is then corrected based on the arrival time difference residual to obtain the final wave arrival position. The leak location is determined based on the final wave arrival position and pipeline geographical information, and the leak location result is output.

[0006] According to some embodiments of the present invention, extracting multipath propagation features includes: The cross-correlation peak width, cross-correlation peak to cross-correlation secondary peak energy ratio, and cross-channel coherence are calculated based on multi-channel infrasound signals to form multipath propagation characteristics. The propagation compensation parameters are determined, including determining the arrival time difference offset correction and channel weighting coefficients based on the multipath propagation characteristics. The multi-channel infrasound signals are corrected based on the propagation compensation parameters, including performing arrival time difference offset correction based on the arrival time difference offset correction and performing channel weighting based on the channel weighting coefficients to obtain the infrasound signal set.

[0007] According to some embodiments of the present invention, calculating the time difference of arrival set based on the infrasound signal set includes: Perform cross-correlation operation on any two channels of infrasound signals in the infrasound signal set, determine the arrival time difference corresponding to the main cross-correlation peak, and collect all arrival time differences to form an arrival time difference set.

[0008] According to some embodiments of the present invention, estimating the initial azimuth of arrival based on the set of time differences of arrival includes: A candidate set of arrival directions is constructed based on array geometric parameters. The theoretical time difference of arrival (TDOA) set is calculated for the arrival directions in the candidate set. The arrival direction with the smallest difference between the TDOA set and the theoretical TDOA set is taken as the initial arrival direction.

[0009] According to some embodiments of the present invention, updating array geometric parameters based on an initial set of azimuth and time difference of arrival includes: Construct an array configuration set, which consists of multiple candidate combinations of array geometric parameters generated based on sensor installation location constraints. Calculate the theoretical arrival time difference set for each array configuration in the array configuration set, and calculate the difference between the arrival time difference set and the corresponding theoretical arrival time difference set to obtain a consistency evaluation value. Select the array configuration with the smallest consistency evaluation value as the updated array geometric parameters.

[0010] According to some embodiments of the present invention, obtaining the fused orientation through weighted fusion includes: Based on the array configuration set, the azimuth of arrival (ADI) of each array configuration in the array configuration set is estimated. Based on the consistency evaluation value of each array configuration, the configuration weight is determined. Based on the configuration weight, each ADI is weighted and fused to obtain the fused ADI.

[0011] According to some embodiments of the present invention, the azimuth filter is a Kalman filter, and the noise parameters of the azimuth filter include process noise parameters and measurement noise parameters.

[0012] According to some embodiments of the present invention, updating the azimuth filter noise parameters based on historical azimuth-arrival reference sequences includes: The difference between the fused azimuth and the azimuth predicted by the azimuth filter is used to form an observation residual sequence. The process noise parameters and measurement noise parameters are updated based on the observation residual sequence to output a refined azimuth of arrival.

[0013] According to some embodiments of the present invention, calculating the time difference of arrival residual based on refined azimuth includes: The theoretical time difference of arrival (TDOA) set is calculated based on the refined azimuth and updated array geometry parameters. The TDOA residual is formed by the difference between the TDOA set and the theoretical TDOA set. The refined azimuth is then corrected based on the TDOA residual by: determining the azimuth correction amount based on the TDOA residual using the least squares method; and correcting the refined azimuth based on the azimuth correction amount to obtain the final azimuth. Determining the leak location based on the final azimuth of arrival and pipeline geographic information involves using the sensor installation location corresponding to the multi-channel infrasound signal as the monitoring point, generating an azimuth line along the final azimuth of arrival, and determining the intersection of the azimuth line and the pipeline centerline represented by the pipeline geographic information as the leak location.

[0014] In a second aspect, the present invention provides an infrasound signal localization and analysis system for implementing an infrasound signal localization and analysis method, the system comprising: The signal correction module is used to acquire multi-channel infrasound signals of oil and gas pipelines, extract multipath propagation characteristics and determine propagation compensation parameters, and correct the multi-channel infrasound signals based on the propagation compensation parameters to obtain an infrasound signal set. The azimuth fusion module is used to calculate the time difference of arrival set based on the infrasound signal set, estimate the initial azimuth of arrival based on the time difference of arrival set, update the array geometric parameters and estimate the azimuth of arrival based on the initial azimuth of arrival and the time difference of arrival set, and obtain the fused azimuth through weighted fusion. The filter update module is used to input the fused azimuth into the azimuth filter, update the noise parameters of the azimuth filter based on the historical azimuth-arrival reference sequence, and output the refined azimuth-arrival. The positioning output module is used to calculate the arrival time difference residual based on the refined wave arrival position, correct the refined wave arrival position according to the arrival time difference residual to obtain the final wave arrival position, determine the leak location based on the final wave arrival position and pipeline geographical information, and output the leak location result.

[0015] Compared with the prior art, the advantages and beneficial effects of the present invention are as follows: By extracting multipath propagation features and generating propagation compensation parameters, the system achieves correction for time-of-arrival (TOA) bias and channel quality differences. Through initial azimuth search driven by TOA and array geometry parameter updates, it achieves self-correction for installation deviations and geometric uncertainties. By determining configuration weights through consistency evaluation values ​​and performing weighted fusion, it achieves robust convergence of multiple azimuth solutions. By updating noise parameters using an azimuth filter and combining historical TOA reference sequences, it achieves suppression of azimuth jitter and anomalies. Finally, by azimuth correction based on TOA residual constraints and intersection with pipeline geographic information, it achieves verifiable high-precision leak location. Attached Figure Description

[0016] Figure 1 This is a flowchart of the infrasound signal localization and analysis method according to an embodiment of the present invention; Figure 2 This is a functional block diagram of the infrasound signal localization and analysis system according to an embodiment of the present invention. Detailed Implementation

[0017] The embodiments of the present disclosure will now be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the disclosure. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of the present disclosure for ease of explanation. However, it will be apparent that one or more embodiments may be practiced without these specific details. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concepts of the present disclosure.

[0018] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit this disclosure. The terms “comprising,” “including,” etc., as used herein indicate the presence of the stated features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0019] All terms used herein (including technical and scientific terms) have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein are to be interpreted in a manner consistent with the context of this specification, and not in an idealized or overly rigid way.

[0020] Infrasound signals typically refer to low-frequency acoustic signals with frequencies below the lower limit of human audibility. They possess propagation characteristics such as long wavelengths, strong diffraction capabilities, and relatively slow attenuation in surface and underground media. In oil and gas pipeline leak scenarios, leak jets and turbulent disturbances excite low-frequency acoustic responses in the pipe and surrounding media. These responses can be collected by smart sensors deployed along the pipeline in a multi-channel manner to extract location information such as time difference of arrival (TDOA) and incident azimuth, and can serve as anomaly evidence inputs for fault prediction and health management. However, infrasound is susceptible to wind noise, mechanical vibration, and terrain reflection in real-world environments, and its propagation path may exhibit multipath superposition, leading to decreased waveform consistency for the same event across different channels. This, in turn, affects the stable extraction of TDOA and the continuity of azimuth estimation. To achieve high-precision leak location, it is necessary to perform propagation compensation, geometric self-correction, and residual constraint azimuth refinement on the multi-channel infrasound signals collected by smart sensors. This ensures that the output azimuth forms a verifiable location link with the pipeline's geographical information, providing traceable and verifiable spatial location data for fault prediction and health management.

[0021] like Figure 1 As shown, a method for locating and analyzing infrasound signals includes the following: The multi-channel infrasound signals of oil and gas pipelines are acquired, multipath propagation characteristics are extracted and propagation compensation parameters are determined, and the multi-channel infrasound signals are corrected based on the propagation compensation parameters to obtain an infrasound signal set. In one embodiment, several infrasound sensors are deployed along the oil and gas pipeline. The outputs of each sensor are uniformly sampled and time-aligned by a data acquisition unit to form a multi-channel infrasound signal. The multi-channel infrasound signal is segmented by time window, and the cross-correlation function is calculated for any two channels. The width of the main cross-correlation peak, the energy ratio of the main peak to the secondary peak, and the cross-channel coherence are extracted as multipath propagation characteristics to characterize the false peak interference and time delay instability caused by multipath propagation. Propagation compensation parameters are determined based on the multipath propagation characteristics. These parameters include the time difference of arrival offset correction and the channel weighting coefficients. The time difference of arrival offset correction is used to perform time correction on the multi-channel infrasound signal, and each channel is weighted according to the channel weighting coefficients to obtain an infrasound signal set. The infrasound signal set contains the corrected signal segments, timestamps, channel weights, and quality markers for subsequent azimuth-of-arrival estimation.

[0022] Extracting multipath propagation features includes: calculating the cross-correlation main peak width, the energy ratio of the cross-correlation main peak to the cross-correlation secondary peak, and the cross-channel coherence based on the multi-channel infrasound signals to form multipath propagation features; determining propagation compensation parameters includes determining the arrival time difference offset correction amount and channel weighting coefficients based on the multipath propagation features; and correcting the multi-channel infrasound signals based on the propagation compensation parameters includes performing arrival time difference offset correction based on the arrival time difference offset correction amount and performing channel weighting based on the channel weighting coefficients to obtain the infrasound signal set.

[0023] In one embodiment, multiple infrasound sensors are deployed along the oil and gas pipeline. The acquisition unit uniformly samples the infrasound signals from each channel and performs time alignment, extracting segments of the multi-channel infrasound signals according to fixed time windows. For each time window, a cross-correlation sequence is calculated for any two channels of infrasound signals, and the main cross-correlation peak and secondary cross-correlation peak are located within the cross-correlation sequence. The width of the main cross-correlation peak is used to characterize the broadening degree of the main cross-correlation peak, which can be obtained by using the time span corresponding to the main cross-correlation peak amplitude decreasing to half its amplitude. The energy ratio of the main cross-correlation peak to the secondary cross-correlation peak is used to distinguish between the main peak and spurious peaks, and the energy ratio can be calculated as the ratio of the energy in the neighborhood of the main peak to the energy in the neighborhood of the secondary peak. in, The ratio of the energy of the main cross-correlation peak to the energy of the secondary cross-correlation peak. For the energy of the neighborhood of the main peak of cross-correlation, This represents the neighborhood energy of the second peak of the cross-correlation. Cross-channel coherence is used to characterize the similarity between two channels in the low-frequency band and can be calculated using the normalized coefficient of the amplitude of the main peak of the cross-correlation. in, For cross-channel coherence, The amplitude of the main cross-correlation peak, The autocorrelation delay value of the first channel is zero. This represents the zero-delay value of the autocorrelation in the second channel. The width of the main cross-correlation peak, the energy ratio between the main cross-correlation peak and the secondary cross-correlation peak, and the cross-channel coherence are summarized in the sensor dimension to form multipath propagation characteristics, which are used to reflect the broadening of the main cross-correlation peak, the enhancement of spurious peaks, and the decrease in channel coherence caused by multipath propagation.

[0024] The propagation compensation parameters include the time difference of arrival (TDOA) offset correction and channel weighting coefficients. The TDOA offset correction is established for each sensor pair. It utilizes the TDOA drift sequence obtained from the changes in the position of the main cross-correlation peak between adjacent time windows. Outliers in the drift sequence are removed, and robust statistics are used as the TDOA offset correction to eliminate systematic offsets introduced by multipath propagation. The channel weighting coefficients are established for each channel. Channel quality is scored based on cross-channel coherence and the energy ratio of the main cross-correlation peak to the secondary cross-correlation peak. Channels with higher scores are assigned larger channel weighting coefficients, and channels with lower scores are assigned smaller channel weighting coefficients. The channel weighting coefficients are then normalized to maintain amplitude scale stability.

[0025] When correcting multi-channel infrasound signals based on propagation compensation parameters, the corresponding channel segments are first time-corrected according to the time difference of arrival (TDOA) offset correction. Time correction can be achieved through interpolation resampling or fractional delay filtering to support fine-grained correction of non-integer sampling points. Then, the multi-channel infrasound signals are weighted according to channel weighting coefficients, outputting an infrasound signal set. The infrasound signal set includes at least the corrected multi-channel infrasound signal segments, channel weighting coefficients, and time window identifiers, used for subsequent TDOA calculation and azimuth-of-arrival (AOA) estimation. This process reduces the impact of false peak interference and time delay drift on subsequent positioning, improving the consistency and stability of subsequent AOA estimation.

[0026] The time difference of arrival (TDA) set is calculated based on the infrasound signal set. The initial azimuth of arrival (AHA) is estimated based on the TDA set. The array geometric parameters are updated and the AHA is estimated based on the initial AHA and the TDA set. The fused AHA is obtained by weighted fusion. In one embodiment, after inputting a set of infrasound signals, cross-correlation sequences are calculated pairwise for each sensor pair. The time delay corresponding to the main peak of the cross-correlation sequence is used as the arrival time difference, and the sets of arrival time differences are collected. The arrival time difference is obtained by the following formula: in, For the time difference of arrival, It is a two-channel cross-correlation sequence. This represents the delay amount.

[0027] Based on the set of time differences of arrival (TDAs), the theoretical TDAs are calculated one by one from the preset candidate azimuth of arrival (ADI) set, and the ADI with the smallest difference is selected as the initial ADI. Subsequently, a candidate set of array geometric parameters is constructed around the nominal values ​​of the array geometric parameters. ADI estimation is repeatedly performed on each candidate set of array geometric parameters to obtain multiple ADIs and corresponding difference values. The array geometric parameter with the smallest difference value is selected as the updated array geometric parameter. Based on the updated array geometric parameters, multiple ADIs are obtained. Weights are determined according to the difference values, and the ADIs are converted into unit direction vectors, weighted, summed, and back-calculated to obtain the fused ADI, thereby reducing the impact of multipath residues and geometric biases on ADI estimation.

[0028] The calculation of the arrival time difference set based on the infrasound signal set includes: performing cross-correlation operation on any two channels of infrasound signals in the infrasound signal set, determining the arrival time difference corresponding to the main peak of cross-correlation, and collecting all arrival time differences to form the arrival time difference set.

[0029] In one embodiment, the infrasound signal set consists of multi-channel infrasound signal segments, time window identifiers, and channel quality markers. To obtain stable time difference of arrival (TDOA) information, the infrasound signal segments of each channel in the infrasound signal set are first traversed according to the time window. For any two channels, a sensor pair is formed, and cross-correlation operations are performed on the sensor pairs to obtain a cross-correlation sequence. The cross-correlation operation is implemented using a sliding delay matching method: taking the first channel infrasound signal segment as a reference, the second channel infrasound signal segment is shifted point by point within a preset delay search range, and the correlation value at each shift position is calculated, forming a cross-correlation sequence that varies with delay. To reduce the impact of amplitude differences, the cross-correlation sequence is normalized to ensure that the amplitude of the cross-correlation sequence falls within a uniform scale.

[0030] The main cross-correlation peak is used to characterize the optimal alignment delay between the two channels within a given time window. The location of the main cross-correlation peak is determined by the global maximum value of the cross-correlation sequence, and a refined search is performed within the neighborhood of this maximum value to obtain a more stable peak position. When multiple local peaks exist in the cross-correlation sequence, the main peak is confirmed in conjunction with channel quality marking: priority is given to peaks whose energy ratio between the main and secondary peaks meets a threshold condition, and the peak neighborhood is smoothed to suppress spike noise. If the main cross-correlation peak does not meet the consistency requirements within a certain time window, the arrival time difference corresponding to that time window is marked as invalid to avoid introducing occasional interference into subsequent azimuth estimation.

[0031] The arrival time difference is obtained in the same way as the formula defined above, where the arrival time difference is obtained from the delay corresponding to the main cross-correlation peak. To support the arrival time difference of non-integer sampling points, interpolation refinement is performed in the neighborhood of the main cross-correlation peak. The interpolation method adopts parabolic fitting or cubic interpolation to obtain a sub-sampling accuracy delay estimate.

[0032] The cross-correlation calculation and peak localization process are repeated for all sensor pairs within the same time window to obtain a set of arrival time differences (OTDs). The OTDs from multiple time windows are then aggregated to form an OTD set. Each element in the OTD set contains a sensor pair identifier, a time window identifier, an OTD value, and a validity marker, used for subsequent initial azimuth-of-arrival (AOA) estimation and array geometry parameter updates. By limiting the extraction of OTDs to cross-correlation peak localization and introducing a validity marker, mismatches can be reduced in the presence of multipath and noise, improving the stability of subsequent AOA estimation.

[0033] Estimating the initial azimuth of arrival based on the time difference of arrival set includes: constructing a candidate set of azimuths of arrival based on array geometric parameters; calculating the theoretical time difference of arrival set for the azimuths of arrival in the candidate set; and taking the azimuth of arrival with the smallest difference between the time difference of arrival set and the theoretical time difference of arrival set as the initial azimuth of arrival.

[0034] In one embodiment, the time difference of arrival set already includes the time difference of arrival values ​​and validity markers for multiple sensor pairs. The array geometry parameters are determined by the sensor installation locations, which can be provided by construction survey coordinates or field calibration data, and stored in the system in two-dimensional planar coordinates for subsequent calculation of the theoretical time difference of arrival for each sensor pair under different incident directions.

[0035] To estimate the initial azimuth of arrival (AHA), a candidate AHA set is first constructed. This candidate AHA set can be generated discretely at equal azimuth intervals, covering the entire target monitoring area. When the pipeline direction is known, the candidate set can also be generated only in the azimuth interval near the pipeline's normal to reduce computational load. Subsequently, for each AHA in the candidate AHA set, the theoretical time difference of arrival (TDOA) is calculated based on the array's geometric parameters. The calculation of the theoretical TDOA is based on the plane wave approximation, which is sufficient for engineering applications when the distance to the leak source is much larger than the array aperture. The propagation velocity is a preset propagation velocity, which can be obtained through field calibration or calculated from environmental parameters.

[0036] To achieve the "minimum difference" criterion, the following difference evaluation function is used to score each candidate azimuth of arrival: in, This is the evaluation value for the difference in azimuth of arrival. For candidate wave direction, For an effective sensor pair set, For sensor weights, For the arrival time difference in the set of arrival time differences, The bearing of the port is The theoretical time difference of arrival. Sensor pair weights can be determined by channel quality labels, with higher-quality sensor pairs assigned greater weights; sensor pairs with invalid labels are not included in the set of valid sensor pairs. To suppress occasional interference, consistency screening of the time differences of arrival can be performed before calculating the difference evaluation value, for example, removing time differences of arrival that deviate too much from the median.

[0037] The difference evaluation value is calculated for each candidate azimuth of arrival (ADI) and the ADI corresponding to the smallest difference evaluation value is selected as the initial ADI. The initial ADI can provide a stable initial ADI value without relying on a complex model, providing a reliable starting point for subsequent array geometry parameter updates and ADI refinement, thereby reducing the risk of ADI jumps caused by multipath residues and geometric biases.

[0038] Updating array geometric parameters based on the initial azimuth and time difference of arrival set includes: constructing an array configuration set, which consists of multiple candidate combinations of array geometric parameters generated based on sensor installation location constraints; calculating the theoretical time difference of arrival set for each array configuration in the array configuration set; calculating the difference between the arrival time difference set and the corresponding theoretical time difference of arrival set to obtain a consistency evaluation value; and selecting the array geometric parameters corresponding to the array configuration with the smallest consistency evaluation value as the updated array geometric parameters.

[0039] In one embodiment, array geometry parameters are used to describe the installation positions of each sensor and the relative geometric relationships between the sensors. Since on-site construction errors, ground subsidence, or equipment movement may cause deviations between the installation positions and nominal values, directly using the nominal array geometry parameters will introduce inconsistencies between the theoretical and actual time differences of arrival (TDAs), thus affecting the stability of the azimuth of arrival (HAA). Therefore, after obtaining the initial HAA, the array geometry parameters are self-calibrated and updated using the initial HAA and the set of TDAs.

[0040] An array configuration set represents a group of candidate array geometry parameters. The array configuration set is constructed centered on the nominal array geometry parameters, introducing small-scale positional perturbations at each sensor installation location to form several candidate configurations. Positional perturbations can be generated using grid enumeration or random sampling. When there are many sensors, sensors with higher channel quality ratings can be kept fixed, while only sensors with lower channel quality ratings or greater uncertainties in the field are perturbed, thus reducing the number of combinations. Each candidate configuration contains the candidate coordinates of all sensor installation locations, thus corresponding to a set of candidate array geometry parameters.

[0041] For each set of candidate array geometric parameters in the array configuration set, a theoretical time difference of arrival (TDOA) set is calculated based on the initial azimuth of arrival (ADI). The calculation method for the theoretical TDOA set is consistent with that used in the initial ADI estimation, only replacing the array geometric parameters with those of the current candidate arrays. Subsequently, the TDOA set is compared with the current theoretical TDOA set to obtain a consistency evaluation value. The consistency evaluation value can be calculated using the difference evaluation function described earlier, and summarized across the set of effective sensor pairs to ensure that the evaluation value is only affected by the reliable TDOA. To enhance robustness, weights can be introduced into the difference calculation, determined by the channel quality label of the sensor pair, allowing higher-quality sensor pairs to contribute more to the consistency evaluation value.

[0042] A consistency evaluation value is calculated for each array configuration set, and the candidate array geometric parameters corresponding to the smallest consistency evaluation value are selected as the updated array geometric parameters. To avoid unreasonable jumps in the update results, a smoothing constraint is applied to the updated array geometric parameters, meaning that the offset of the updated sensor installation position relative to the nominal position does not exceed a preset upper limit, and candidate configurations exceeding the limit are directly eliminated. After the update is completed, the updated array geometric parameters are written into the processing context of the current time window for subsequent re-estimation and weighted fusion of azimuth and arrival. Through the above update process, array installation errors can be compensated without increasing the workload of on-site calibration, improving the consistency between the theoretical time difference of arrival and the actual time difference of arrival, thereby reducing the sensitivity of azimuth estimation to array geometric deviations.

[0043] The weighted fusion method for obtaining the fused azimuth involves: estimating the azimuth of arrival (ADR) of each array configuration in the array configuration set; determining the configuration weight based on the consistency evaluation value of each array configuration; and performing weighted fusion of each ADR based on the configuration weight to obtain the fused azimuth.

[0044] In one embodiment, after updating the array geometry parameters, to reduce the impact of single geometry parameter estimation errors on the azimuth of arrival (ADI), an array configuration set is determined based on the updated array geometry parameters. The array configuration set is generated by introducing small-scale positional perturbations centered on the updated sensor installation locations, resulting in multiple candidate installation location combinations. Each candidate installation location combination corresponds to a set of array geometry parameters. Positional perturbations can be implemented using grid enumeration or random sampling, and priority can be given to keeping sensors with higher channel quality labels unchanged, while introducing perturbations only for sensors with greater uncertainty, in order to control the number of candidates and ensure computational availability.

[0045] For each set of array geometric parameters in the array configuration set, the azimuth estimation process is repeatedly performed using the same set of time differences of arrival (TDAs) to obtain a set of candidate azimuth results. For each set of array geometric parameters, the difference between the set of TDAs and the theoretical set of TDAs corresponding to that set of array geometric parameters is calculated simultaneously to obtain a consistency evaluation value. The consistency evaluation value is calculated in the form of the aforementioned difference evaluation function and is only summarized for valid sensor pairs to avoid interference from invalid TDAs on the evaluation results.

[0046] Configuration weights are determined by consistency evaluation values. Arrays with smaller consistency evaluation values ​​are assigned larger configuration weights. Configuration weights are obtained using reciprocal normalization. in, For the first The configuration weights of the group array configuration. For the first Consistency evaluation value of group array configuration. To prevent positive constants with a denominator of zero, Configure the array with sequence number The consistency evaluation value corresponding to the array configuration. The total number of array configurations in the array configuration set is determined. After calculating the configuration weights, the azimuths of arrival (AHAs) of each group are converted into planar unit direction vectors. These vectors are then weighted and summed according to their configuration weights to obtain a weighted composite direction vector. This weighted composite direction vector is then converted into a fused azimuth. Through this weighted fusion method, the fused azimuth can average and suppress residual errors in array geometric parameters and fluctuations in time difference of arrival (TDOA), thereby improving the stability and consistency of azimuth estimation and providing a more reliable input for subsequent filtering updates and leak localization.

[0047] The azimuth input azimuth filter is fused, and the azimuth filter noise parameters are updated based on the historical azimuth reference sequence to output a refined azimuth. In one embodiment, the fused azimuth is generated continuously according to time windows, and may exhibit jitter or jumps due to transient interference and residual multipath effects. To obtain a stable azimuth usable for positioning, the fused azimuth is input to an azimuth filter for temporal smoothing. The azimuth filter is implemented using a Kalman filter, with the fused azimuth of the current time window as input and the refined azimuth as output. The historical azimuth reference sequence consists of refined azimuths from previous time windows, used to characterize the normal fluctuation level of azimuth changes.

[0048] The azimuth filter noise parameters include process noise parameters and measurement noise parameters. Process noise parameters describe the uncertainty of azimuth changes over time, while measurement noise parameters describe the measurement fluctuations in the fused azimuth. When updating the noise parameters, the observation residuals between the fused azimuth and the filter-predicted azimuth are first calculated to form an observation residual sequence. Then, the measurement noise parameters are adjusted based on the statistics of the observation residual sequence. Simultaneously, the process noise parameters are adjusted based on the azimuth increment statistics of the historical azimuth-arrival reference sequence. This enhances the smoothing effect of the filter when the azimuth is stable and maintains its tracking capability when the azimuth changes rapidly. Through this updating method, the influence of occasional abnormal azimuths on the results can be suppressed, resulting in a refined azimuth-arrival output with better continuity.

[0049] The azimuth filter is a Kalman filter, and the noise parameters of the azimuth filter include process noise parameters and measurement noise parameters.

[0050] In one embodiment, a Kalman filter is used as the azimuth filter to transform the fused azimuth into a continuous and stable azimuth output. The Kalman filter operates in discrete time, with each time window corresponding to a filter update. The input is the fused azimuth of the current time window, and the output is the refined azimuth of arrival. For ease of engineering implementation, the filter state is composed of the azimuth angle and its change. In the prediction phase, the current azimuth is predicted based on the state of the previous time window, and in the update phase, the prediction result is corrected based on the fused azimuth, thereby achieving a balance between smoothing jitter and maintaining response speed.

[0051] Azimuth filter noise parameters include process noise parameters and measurement noise parameters. Process noise parameters describe the uncertainty of natural azimuth changes between adjacent time windows, reflecting the magnitude of azimuth changes caused by leakage sources, environmental disturbances, or changes in propagation conditions. Measurement noise parameters describe the measurement uncertainty of the fused azimuth, reflecting the impact of time-of-arrival fluctuations, multipath propagation, and residual errors in array geometry parameters on the fused azimuth. A larger process noise parameter results in a stronger ability to track rapid azimuth changes but a weaker smoothing effect; a larger measurement noise parameter results in a lower level of confidence in the fused azimuth but a stronger ability to suppress outliers.

[0052] In the system implementation, both process noise parameters and measurement noise parameters are stored in a configurable format and support dynamic updates by time window. To ensure availability in the initial stage, preset initial values ​​can be used for initialization when the filter starts. After accumulating a certain length of historical azimuth-of-arrival (AHA) reference sequence, the noise parameters are adaptively adjusted based on the fluctuation characteristics of the historical AHA reference sequence and the fused azimuth, enabling the filter to adapt to different pipe sections, different construction conditions, and different background noise levels. By dividing the noise parameters into process noise parameters and measurement noise parameters, the filter can suppress random jitter in the fused azimuth and maintain reasonable tracking capability when the leak azimuth changes, thereby outputting a more continuous and refined AHA.

[0053] Updating the azimuth filter noise parameters based on the historical azimuth reference sequence includes: forming an observation residual sequence by the difference between the fused azimuth and the azimuth predicted by the azimuth filter, and updating the process noise parameters and measurement noise parameters based on the observation residual sequence to output a refined azimuth.

[0054] In one embodiment, the Kalman filter includes a prediction phase and an update phase in each time window. The prediction phase outputs a predicted azimuth based on the filtering state of the previous time window. The predicted azimuth reflects a priori estimate of the current azimuth without incorporating the current fused azimuth. To enable the filter to adaptively adjust the smoothing intensity as environmental and signal quality changes, a historical azimuth-arrival reference sequence is introduced to update the noise parameters online. The historical azimuth-arrival reference sequence consists of refined azimuths from the previous several time windows and includes azimuth increment information for the corresponding time windows, reflecting the normal fluctuation range of the azimuth on the time axis.

[0055] Noise parameter updates are based on the observation residual sequence. The observation residual is determined by the difference between the fused azimuth and the predicted azimuth, and is continuously recorded along the time window to form the observation residual sequence. The observation residual sequence can directly reflect the degree of deviation of the fused azimuth from the predicted azimuth: when the overall observation residual sequence is small and the fluctuation is stable, it indicates that the fused azimuth is reliable and the azimuth change is slow; when the observation residual sequence shows a sudden increase or a continuous increase, it indicates that the fused azimuth is more affected by noise or multipath residue, or that the azimuth has indeed changed rapidly, and the noise parameters need to be adjusted accordingly.

[0056] The updating of measurement noise parameters is based on the statistics of the observation residual sequence. The mean and dispersion of the observation residual sequence are calculated within a sliding window, and outliers with excessive dispersion are removed to avoid drastic fluctuations in noise parameters caused by a single anomaly. After outlier removal, the dispersion of the observation residual sequence is used as the basis for updating the measurement noise parameters, allowing them to adaptively increase or decrease with the level of fused azimuth fluctuations. The updating of process noise parameters is based on the azimuth increment statistics of the historical wave arrival-azimuth reference sequence. A azimuth increment sequence is formed from the difference between adjacent refined wave arrival-azimuth, and the dispersion of the azimuth increment sequence characterizes the uncertainty of azimuth change, thereby updating the process noise parameters to reflect the rate of change of the true azimuth.

[0057] After updating the noise parameters, the filter uses the fused azimuth to correct the predicted azimuth during the update phase, outputting a refined azimuth of arrival. By updating the measurement noise parameters and process noise parameters in conjunction with the observation residual sequence, the filter automatically enhances its suppression of abnormal inputs when the fused azimuth fluctuations increase, and automatically enhances its tracking capability when the actual azimuth changes rapidly. Thus, it can output a refined azimuth of arrival with better continuity and less jitter under different operating conditions, providing a stable input for subsequent time-of-arrival residual correction and leakage location.

[0058] The arrival time difference residual is calculated based on the refined wave arrival position. The refined wave arrival position is then corrected based on the arrival time difference residual to obtain the final wave arrival position. The leak location is determined based on the final wave arrival position and pipeline geographical information, and the leak location result is output.

[0059] In one embodiment, after obtaining the refined azimuth of arrival (RAA), the theoretical time difference of arrival (TDOA) set is calculated based on the updated array geometry parameters. The TDOA set is then paired with the theoretical TDOA set to obtain the TDOA residuals. These residuals characterize the degree of theoretical matching at the current azimuth. A systematic deviation in the TDOA residuals indicates that the refined azimuth of arrival still has correctable biases. The azimuth of arrival correction is determined based on the TDOA residuals. This correction can be obtained using least-squares fitting, and the correction is then superimposed on the refined azimuth of arrival to obtain the final azimuth of arrival. The final azimuth of arrival is used for leak location: the coordinates of the monitoring point and the final azimuth of arrival are converted into azimuth lines in a geographic coordinate system, and their intersection with the pipeline centerline in the pipeline geographic information is determined. The intersection point is used as the leak location. When multiple intersection points exist, they are selected based on the consistency between the distance from the intersection point to the monitoring point and the azimuth residual, thus determining the leak location and outputting the leak location result.

[0060] The calculation of the arrival time difference residual based on the refined azimuth of arrival includes: calculating the theoretical arrival time difference set based on the refined azimuth of arrival and the updated array geometry parameters, and forming the arrival time difference residual by the difference between the arrival time difference set and the theoretical arrival time difference set. Correcting the refined azimuth of arrival based on the arrival time difference residual includes: determining the azimuth of arrival correction amount based on the arrival time difference residual using the least squares method; and correcting the refined azimuth of arrival based on the azimuth of arrival correction amount to obtain the final azimuth of arrival. Determining the leak location based on the final azimuth of arrival and pipeline geographic information involves using the sensor installation location corresponding to the multi-channel infrasound signal as the monitoring point, generating an azimuth line along the final azimuth of arrival, and determining the intersection of the azimuth line and the pipeline centerline represented by the pipeline geographic information as the leak location.

[0061] In one embodiment, the refined azimuth of arrival (RAA) is used to provide a stable estimate of the current incident direction of the leaking sound source. However, when residual biases exist in the array geometry parameters and fluctuations occur in the time difference of arrival (TDA), the RAA may still carry a small systematic error. To further improve the consistency between the azimuth and the TDA, a theoretical TDA set is calculated using the RAA and the updated array geometry parameters. The theoretical TDA set is based on the plane wave approximation, with the updated sensor installation position and the RAA as inputs, and the output as the theoretical TDA for each effective sensor pair at that azimuth. Subsequently, the TDA set is subtracted from the theoretical TDA set pairwise, and the difference forms the TDA residual. The TDA residual reflects the fitting error distribution under the current RAA. When the residuals show an overall unidirectional shift or exhibit a regular change with the sensor pair, it indicates that there is a correctable bias in the RAA.

[0062] When refining the bearing of arrival (BHA) based on the time difference of arrival (TDOA) residuals, the BHA correction is used as the estimated variable. A linear relationship is established between the residuals and the BHA correction, and the least squares method is used to solve for the BHA correction, minimizing the difference between the corrected theoretical TDOA set and the actual TDOA set. For ease of engineering implementation, the least squares solution can be iterative: using the refined BHA as the initial value, the TDOA residual is calculated once, and the BHA correction is solved. After updating the BHA, the residual is calculated again, until the BHA correction is less than a preset threshold or the number of iterations reaches the upper limit. The obtained BHA correction is then superimposed on the refined BHA to obtain the final BHA. Through residual-driven BHA correction, the TDOA information can be reconstrained into the BHA estimation, reducing the accumulation of small deviations caused by filtering and smoothing, and improving the interpretability and verifiability of the BHA solution.

[0063] In the leak location phase, the sensor installation locations corresponding to the multi-channel infrasound signals are used as monitoring points. These monitoring points are represented using geographic coordinates or projected coordinates, sourced from field measurement records or station locations exported from a geographic information system. An azimuth line is generated at each monitoring point based on the final azimuth of arrival, extending in the geographic coordinate system along the direction of the final azimuth. Pipeline geographic information includes at least a vector representation of the pipeline centerline, formed by connecting several geographic nodes, and may include mileage information for result presentation. The intersection of the azimuth line and the pipeline centerline is calculated, and the intersection point is taken as the leak location, outputting the leak location result. When multiple intersection points exist between the azimuth line and the pipeline centerline, the distance between the monitoring point and the intersection point, the consistency of the arrival time difference residual, and the continuity of the pipeline direction can be considered to filter the points, selecting the intersection point with better consistency as the leak location. This ensures stable location output even in complex pipe networks or broken pipe sections.

[0064] like Figure 2 As shown, an infrasound signal localization analysis system is used to implement the above-described infrasound signal localization analysis method. The system includes: The signal correction module is used to acquire multi-channel infrasound signals from oil and gas pipelines, extract multipath propagation characteristics, determine propagation compensation parameters, and correct the multi-channel infrasound signals based on the propagation compensation parameters to obtain an infrasound signal set. The signal correction module hardware includes a multi-channel infrasound sensor array, a pre-amplifier analog conditioning circuit, a synchronous sampling analog-to-digital converter circuit, and an edge computing unit. The multi-channel infrasound sensors can be micro-pressure sensors or low-frequency accelerometers, deployed along the pipeline or around valve chambers and stations. The pre-amplifier analog conditioning circuit includes low-noise amplification, switchable gain, band-limited filtering, and overvoltage protection to improve the availability of weak signals and suppress power frequency and mechanical vibration interference. The synchronous sampling analog-to-digital converter circuit uses a shared clock or the same sampling trigger line to achieve channel synchronization, ensuring that the arrival time difference can be calculated. The edge computing unit has a built-in high-precision clock and buffer to locally extract multipath propagation characteristics and generate propagation compensation parameters, perform time delay correction and channel weighting on the sampled data, and output an infrasound signal set for subsequent processing.

[0065] The azimuth fusion module is used to calculate the time difference of arrival (TDOA) set based on the infrasound signal set, estimate the initial azimuth of arrival (AHA) based on the TDOA set, update the array geometry parameters and estimate the azimuth of arrival (AHA) based on the initial AHA and the TDOA set, and obtain the fused azimuth through weighted fusion. The azimuth fusion module includes a high-performance processor or embedded computing board, optional hardware acceleration units, and inter-channel clock consistency maintenance circuitry. The processor performs cross-correlation calculations, generates the TDOA set, searches for the azimuth of arrival (AHA), and updates the array geometry parameters. When the number of channels is large or the time window is short, on-chip vector instructions, digital signal processors, or graphics processors can be used to achieve parallel acceleration of cross-correlation and consistency evaluation. This module works with a unified time reference circuit to obtain the channel sampling timestamps, combines the stored sensor installation locations or nominal coordinates to generate a candidate set of array geometry parameters, generates configuration weights based on the consistency evaluation values, and performs weighted fusion of multiple sets of azimuths of arrival (AHA) to output the fused azimuth, ensuring computational throughput and real-time performance at the hardware level.

[0066] The filter update module is used to input the fused azimuth data into the azimuth filter, update the noise parameters of the azimuth filter based on the historical azimuth-arrival reference sequence, and output a refined azimuth-arrival. The filter update module includes a real-time operating environment, non-volatile memory, and a reliable time window scheduling unit. The real-time operating environment triggers the prediction and update calculations of the azimuth filter according to fixed time windows. The non-volatile memory stores the historical azimuth-arrival reference sequence, initial noise parameter values, and updated noise parameters, ensuring a rapid recovery of stable output after a power outage and restart. The time window scheduling unit is bound to the system clock to ensure that the fused azimuth input and historical sequence updates occur on the same time reference. When hardware resources permit, this module can be configured with a double-buffered structure, receiving the fused azimuth data stream while simultaneously performing observation residual statistics and noise parameter updates, outputting continuous refined azimuth-arrival and reducing the impact of computational jitter on real-time performance.

[0067] The positioning output module is used to calculate the time difference of arrival residual based on the refined azimuth, correct the refined azimuth to obtain the final azimuth, determine the leak location based on the final azimuth and pipeline geographic information, and output the leak location result. The positioning output module includes a geographic information storage and interface unit, a positioning calculation unit, a communication output unit, and a human-machine interface or alarm interface. The geographic information storage and interface unit stores pipeline centerline vector data, mileage marker information, and monitoring point coordinates, and supports synchronous updates with the pipeline geographic information system via Ethernet, cellular communication, or industrial bus. The positioning calculation unit calculates the time difference of arrival residual based on the refined azimuth and completes the azimuth correction, then generates an azimuth line in the geographic coordinate system and finds its intersection to obtain the leak location. The communication output unit reports the leak location result, confidence information, and timestamp to the monitoring platform and supports interface with valve control and linkage alarm systems. The human-machine interface or alarm interface may include local indicator lights, buzzers, displays, or relay outputs for rapid alarm and on-site handling support at the station.

[0068] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects.

[0069] The above are merely embodiments of the present invention and are not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principle of the present invention should be included within the scope of the claims of the present invention.

Claims

1. A method for localizing and analyzing infrasound signals, characterized in that, Including the following: Multi-channel infrasound signals of oil and gas pipelines are acquired, multipath propagation characteristics are extracted and propagation compensation parameters are determined, and the multi-channel infrasound signals are corrected based on the propagation compensation parameters to obtain an infrasound signal set. Based on the infrasound signal set, calculate the time difference of arrival set, estimate the initial azimuth of arrival based on the time difference of arrival set, update the array geometric parameters and estimate the azimuth of arrival based on the initial azimuth of arrival and the time difference of arrival set, and obtain the fused azimuth through weighted fusion. The fused azimuth input is used to the azimuth filter, and the noise parameters of the azimuth filter are updated based on the historical azimuth reference sequence to output a refined azimuth. The arrival time difference residual is calculated based on the refined wave arrival position. The refined wave arrival position is then corrected based on the arrival time difference residual to obtain the final wave arrival position. The leak location is determined based on the final wave arrival position and pipeline geographical information, and the leak location result is output.

2. The infrasound signal localization and analysis method according to claim 1, characterized in that, Extracting the multipath propagation features includes: The multipath propagation characteristics are formed by calculating the cross-correlation peak width, the energy ratio of the cross-correlation peak to the cross-correlation secondary peak, and the cross-channel coherence based on the multipath propagation characteristics. The propagation compensation parameters are determined by determining the arrival time difference offset correction and the channel weighting coefficients based on the multipath propagation characteristics. The multipath propagation compensation parameters are used to correct the multipath propagation characteristics by performing arrival time difference offset correction based on the arrival time difference offset correction and performing channel weighting based on the channel weighting coefficients to obtain the infrasound signal set.

3. The infrasound signal localization and analysis method according to claim 1, characterized in that, The calculation of the time difference of arrival set based on the infrasound signal set includes: Perform cross-correlation operation on any two channels of infrasound signals in the infrasound signal set, determine the arrival time difference corresponding to the main cross-correlation peak, and collect the arrival time differences to form the arrival time difference set.

4. The infrasound signal localization analysis method according to claim 1, characterized in that, Estimating the initial azimuth of arrival based on the set of time differences of arrival includes: Based on the array geometry parameters, a candidate set of arrival directions is constructed. A theoretical time difference of arrival (TDOA) set is calculated for the arrival directions in the candidate set. The arrival direction with the smallest difference between the TDOA set and the theoretical TDOA set is taken as the initial arrival direction.

5. The infrasound signal localization analysis method according to claim 1, characterized in that, Updating the array geometric parameters based on the initial azimuth of arrival and the time difference of arrival set includes: A set of array configurations is constructed, which consists of multiple candidate combinations of array geometric parameters generated based on sensor installation location constraints. For each array configuration in the set of array configurations, a set of theoretical arrival time differences is calculated. The difference between the set of arrival time differences and the corresponding theoretical arrival time difference set is calculated to obtain a consistency evaluation value. The array geometric parameters corresponding to the array configuration with the smallest consistency evaluation value are selected as the updated array geometric parameters.

6. The infrasound signal localization analysis method according to claim 5, characterized in that, The fused orientation obtained through weighted fusion includes: Based on the array configuration set, the azimuth of arrival (ADR) corresponding to each array configuration in the array configuration set is estimated respectively. Based on the consistency evaluation value corresponding to each array configuration, the configuration weight is determined. Based on the configuration weight, each ADR is weighted and fused to obtain the fused ADR.

7. The infrasound signal localization analysis method according to claim 1, characterized in that, The azimuth filter is a Kalman filter, and the noise parameters of the azimuth filter include process noise parameters and measurement noise parameters.

8. The infrasound signal localization analysis method according to claim 7, characterized in that, Updating the azimuth filter noise parameters based on the historical azimuth reference sequence includes: The difference between the fused azimuth and the azimuth predicted by the azimuth filter is used to form an observation residual sequence, and the process noise parameter and the measurement noise parameter are updated based on the observation residual sequence to output the refined azimuth of arrival.

9. The infrasound signal localization analysis method according to claim 1, characterized in that, The calculation of the arrival time difference residual based on the refined wave azimuth includes: Based on the refined azimuth of arrival (RAA) and the updated array geometry parameters, a theoretical time difference of arrival (TDOA) set is calculated. The TDOA residual is formed by the difference between the RAOA set and the theoretical RAOA set. Correcting the refined RAA based on the TDOA residual includes: determining the RAA correction amount using the least squares method based on the TDOA residual; and correcting the refined RAA based on the RAA correction amount to obtain the final RAA. Determining the leak location based on the final azimuth of arrival and pipeline geographic information includes using the sensor installation location corresponding to the multi-channel infrasound signal as the monitoring point, generating an azimuth line along the final azimuth of arrival, and determining the intersection of the azimuth line and the pipeline centerline represented by the pipeline geographic information as the leak location.

10. A system for locating and analyzing infrasound signals, used to implement the method for locating and analyzing infrasound signals according to any one of claims 1 to 9, characterized in that, The system includes: The signal correction module is used to acquire multi-channel infrasound signals of oil and gas pipelines, extract multipath propagation characteristics and determine propagation compensation parameters, and correct the multi-channel infrasound signals based on the propagation compensation parameters to obtain an infrasound signal set. The azimuth fusion module is used to calculate the time difference of arrival set based on the infrasound signal set, estimate the initial azimuth of arrival based on the time difference of arrival set, update the array geometric parameters and estimate the azimuth of arrival based on the initial azimuth of arrival and the time difference of arrival set, and obtain the fused azimuth through weighted fusion. The filtering update module is used to input the fused azimuth into the azimuth filter, update the noise parameters of the azimuth filter based on the historical azimuth reference sequence, and output the refined azimuth. The positioning output module is used to calculate the time difference of arrival residual based on the refined wave arrival position, correct the refined wave arrival position according to the time difference of arrival residual to obtain the final wave arrival position, determine the leak location based on the final wave arrival position and pipeline geographical information, and output the leak location result.

Images