Method and system for monitoring gas leaks in underground pipe networks

Abstract

The application discloses a gas leakage monitoring method and system for underground pipe network and relates to the technical field of gas pipeline leakage monitoring. The method comprises the following steps: arranging a sound wave sensing array along the underground pipe network, synchronously collecting pipeline wall vibration signals and pipeline internal medium pressure pulsation signals; constructing a time domain response matrix and performing intrinsic mode decomposition, extracting intrinsic mode component sequences of each monitoring node, and calculating mutual information entropy between adjacent node intrinsic mode component sequences; when the instantaneous change rate of the mutual information entropy exceeds a preset threshold, intercepting vibration signal segments and pressure signal segments of the node and adjacent nodes within a preset time window; arranging the segments to generate a signal matrix and calculating a joint covariance matrix, extracting a main eigenvector corresponding to a maximum eigenvalue through eigenvalue decomposition, and determining the relative position of a leakage sound source according to the amplitude distribution of each node component in the main eigenvector. The method can realize leakage anomaly recognition and high-precision sound source positioning through multi-physical field signal fusion.

Description

Technical Field

[0001] This invention relates to the field of gas pipeline leakage monitoring technology, specifically to a method and system for monitoring gas leaks in underground pipeline networks. Background Technology

[0002] Gas leak monitoring in underground pipe networks primarily relies on deploying sensors along the pipeline to capture the acoustic or vibration signals caused by the leak. Current technologies typically use a single pipe wall vibration signal or internal pipe pressure signal for leak identification. These signals are easily affected by environmental noise, fluid flow, and mechanical operation, making it difficult to reliably extract weak early-stage leak characteristics. For location methods, time-of-arrival (TOA) or beamforming techniques are commonly used, requiring a high signal-to-noise ratio from the leak source signal and that the sensor array meets far-field or plane wave assumptions. However, urban underground pipe networks are often space-constrained, with dense bends and branches, making it difficult to fully meet ideal array conditions for sensor placement. This results in insufficient location accuracy and frequent false alarms or missed detections.

[0003] Vibration signal analysis is sensitive to stress waves propagating through pipe walls, but it is easily affected by external impacts such as construction and traffic. While internal pressure pulsation signals can directly reflect changes in medium flow caused by leaks, they are significantly affected by pressure fluctuations from internal turbulence and valve operation. Both types of signals carry specific information about the leak source; processing them separately fails to fully utilize their complementary time-frequency response characteristics to the same leak event, thus losing crucial correlation dimensions that could enhance detection reliability and refine spatial resolution.

[0004] In recent years, most multi-sensor data fusion solutions have remained at the feature-level or decision-level fusion level, lacking joint spatiotemporal modeling of the original waveform signals. This makes it impossible to capture the subtle abrupt changes in the cross-covariance structure between vibration and pressure at the moment of a leak. Faced with the complex multipath effects of buried pipeline leak sound waves propagating along both the pipe wall and the fluid inside the pipe, existing methods struggle to separate the leak-dominant propagation mode and its energy gradient distribution between adjacent nodes from contaminated multi-channel data. This results in the inability to accurately determine the relative location and orientation of the leak source in situations with extremely low signal-to-noise ratios. Summary of the Invention

[0005] This invention provides a method and system for monitoring gas leaks in underground pipeline networks. By synchronously acquiring the dual physical quantities of vibration and pressure and jointly decomposing the intrinsic mode, a dynamic characterization of the mutual information entropy between monitoring nodes is established, enabling reliable triggering of leak anomalies. Furthermore, the spatial energy distribution of the leak sound source is reconstructed by decomposing the principal eigenvectors of the joint covariance matrix, thereby improving the sensitivity and directional accuracy of underground pipeline network leak monitoring.

[0006] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a method and system for monitoring gas leaks in underground pipelines. It synchronously collects vibration signals of the gas pipeline wall and pressure pulsation signals of the medium inside the pipeline through an acoustic wave sensor array. Based on the instantaneous change characteristics of intrinsic mode decomposition and mutual information entropy, it triggers leak event identification. Furthermore, it utilizes the multi-channel synthetic response intensity distribution obtained by joint covariance matrix eigenvalue decomposition to accurately determine the relative location of the leak sound source.

[0007] This method deploys acoustic wave sensor arrays along preset monitoring nodes of the underground pipeline network to simultaneously collect gas pipeline wall vibration signals and pipeline medium pressure pulsation signals at each monitoring node. Based on the pipeline wall vibration signals and pipeline medium pressure pulsation signals at each monitoring node, a time-domain response matrix is ​​constructed with the monitoring node number as the row index and the sampling time point as the column index. By performing intrinsic mode decomposition on the time-domain response matrix, the intrinsic mode component sequences of each monitoring node are extracted, and the mutual information entropy between the intrinsic mode component sequences of each monitoring node and the intrinsic mode component sequences of adjacent monitoring nodes is calculated. This effectively enhances the weak anomaly characteristics caused by leakage and reflects the changes in the information transmission relationship between signals from different nodes.

[0008] When the instantaneous rate of change of mutual information entropy exceeds a preset threshold, it indicates an abnormal abrupt change in the coupling structure of vibration and pressure pulsation signals in that local area of ​​the pipeline network. Therefore, taking this monitoring node as the center, segments of pipeline wall vibration signals and pipeline medium pressure pulsation signals are extracted from this monitoring node and its adjacent monitoring nodes within a preset time window. The extracted signal segments are arranged according to the monitoring node number, generating a first signal matrix and a second signal matrix. The joint covariance matrix of the two is calculated, and then eigenvalue decomposition is performed on the joint covariance matrix to extract the principal eigenvector corresponding to the largest eigenvalue. This principal eigenvector integrates the responses of vibration and pressure to the same leakage source. Based on the component amplitude distribution of the principal eigenvector at each monitoring node, the synthetic response intensity of each monitoring node can be constructed. The monitoring node with the largest synthetic response intensity is the nearest neighbor node to the leakage source. Furthermore, by combining the ratio of the synthetic response intensity of adjacent nodes with the spatial direction vector, the azimuth bias of the leakage source can be further determined, enabling convenient location without prior knowledge of the precise sound velocity model of the pipeline or the leakage aperture.

[0009] As a preferred technical solution of the present invention, at each preset monitoring node, a piezoelectric vibration sensor is coupled and fixed to the outer wall of the pipe, and a pressure transmitter is connected to the internal medium of the pipe through a pressure tapping pipe. Using the same global clock source as a synchronization reference, the piezoelectric vibration sensor and pressure transmitter of each monitoring node are controlled to start collecting data at the same sampling frequency, so that the pipe wall vibration signal and the pressure pulsation signal of the internal medium of each monitoring node have a unified sampling time axis, thereby ensuring that the subsequent mutual information entropy calculation and joint covariance matrix construction have a strict temporal correspondence.

[0010] As a further preferred embodiment of the present invention, empirical mode decomposition is performed on the signal sequence of each row in the time-domain vibration response matrix and the time-domain pressure response matrix to obtain multiple intrinsic mode function components for each monitoring node. The first three intrinsic mode function components of each monitoring node are selected and superimposed to form a reconstructed component sequence. The mutual information value between the reconstructed component sequence of the monitoring node and the reconstructed component sequences of adjacent monitoring nodes within a preset time delay range is calculated, and the maximum value of the mutual information value is taken as the mutual information entropy. This not only preserves the multi-scale fluctuation characteristics related to leakage, but also captures the dynamic evolution of the signal dependence relationship between adjacent nodes through the mutual information entropy.

[0011] As one embodiment of the present invention, when the instantaneous rate of change of the entropy value of a certain monitoring node is greater than a preset threshold, the monitoring node is marked as a trigger node. Taking the triggering time of the trigger node as the endpoint, a pipe wall vibration signal segment and a pipe medium pressure pulsation signal segment of a preset time window length are extracted. Simultaneously, signal segments of all adjacent monitoring nodes of the trigger node within the same preset time window are extracted, so that the analysis window exactly covers the effective arrival time period of the leakage sound source to each sensor, thereby improving the signal-to-noise ratio contributed by the leakage source in the joint covariance matrix.

[0012] When generating the joint covariance matrix, the pipe wall vibration signal segments from all monitoring nodes are first arranged from top to bottom according to the monitoring node number to form the first signal matrix. The pipe internal medium pressure pulsation signal segments from all monitoring nodes are arranged in the same order to form the second signal matrix. Both are then subjected to mean removal processing. The two matrices are then vertically concatenated to form the joint signal matrix. The joint signal matrix is ​​multiplied by its transpose and divided by the number of time points of the signal segments minus one to obtain the joint covariance matrix. This matrix simultaneously contains the independent statistical characteristics of the vibration and pressure signals and their cross-correlation information, allowing the principal eigenvectors to naturally integrate the directional information of the two sensing modes.

[0013] To determine the relative location of the leakage source, a power-law iterative eigenvalue decomposition is performed on the joint covariance matrix to extract the principal eigenvector corresponding to the largest eigenvalue. The dimension of this principal eigenvector is twice the number of monitoring nodes; the first half of the dimension corresponds to the pipe wall vibration signal channel of each monitoring node, and the second half corresponds to the pipe internal medium pressure pulsation signal channel of each monitoring node. The sum of the squares of the amplitudes of the first and second half of the components corresponding to each monitoring node is calculated to obtain the synthetic response intensity of that monitoring node. The monitoring node with the largest synthetic response intensity is identified as the nearest neighbor node to the leakage source. Then, using this nearest neighbor node as the origin, the synthetic response intensities of all its adjacent nodes are weighted and summed with the spatial direction vector of each adjacent node relative to the nearest neighbor node. The direction of the resulting synthetic vector is the azimuth bias of the leakage source. This method does not rely on the precise calibration of propagation velocity and attenuation coefficient in the analytical model, and the location results have high robustness.

[0014] In the above process, the instantaneous rate of change of entropy can be calculated using a five-point difference scheme, which divides the weighted difference mean of the mutual information entropy at the current moment and the mutual information entropy at the previous four moments by the time interval, thereby achieving a balance between sensitivity and noise resistance and timely capturing the inflection point of change in the mutual information structure.

[0015] Corresponding to the method described above, this invention also provides an underground pipeline gas leak monitoring system, including a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the aforementioned underground pipeline gas leak monitoring method. This system can automatically complete the synchronous acquisition of multi-node acoustic data, construction of the time-domain response matrix, intrinsic mode decomposition and mutual information entropy monitoring, leak event triggering discrimination, joint covariance matrix generation and eigenvalue decomposition, and output of the relative position of the leak sound source. While achieving early detection of weak leaks, it provides reliable leak location information, significantly improving the operational safety of underground pipeline gas lines.

[0016] The technical effects and advantages provided by the present invention in the above technical solution are as follows: Simultaneously, vibration signals from the pipe wall and pressure pulsations within the pipe are acquired. Empirical mode decomposition (EMD) is performed on the time-domain response matrices constructed from these two types of signals to extract the first three intrinsic mode function (IMF) components of each monitoring node and reconstruct them into an IMF component sequence. The maximum mutual information value between each monitoring node and its adjacent nodes within a preset time delay range is calculated as the mutual information entropy. The mutual information entropy is then arranged along the time axis and subjected to sliding difference calculation to obtain the instantaneous rate of change of the entropy value. When the instantaneous rate of change exceeds a preset threshold, it indicates a significant jump in the information coupling strength between vibration and pressure within the corresponding monitoring node region. This jump originates from the structural deviation caused by leakage-induced media impact and pipe wall resonance in the mutual information space. By simultaneously associating the structural acoustic propagation characteristics of vibration and the fluid acoustic propagation characteristics of pressure, the dual-physical-quantity fusion measurement of mutual information entropy can effectively suppress environmental noise and flow disturbances that are difficult to avoid when acting on a single physical quantity. This significantly amplifies the weak characteristics induced by leakage at the multi-dimensional information interaction level, thereby capturing early leakage events with a sensitivity far exceeding that of traditional single-signal threshold detection, avoiding misjudging transient interference as leakage.

[0017] Vibration and pressure signal segments from the trigger node and its adjacent nodes within a preset time window are extracted and arranged into a first signal matrix and a second signal matrix according to node number. After mean removal, they are vertically concatenated to form a joint signal matrix, and the joint covariance matrix is ​​calculated. Power-law iterative eigenvalue decomposition is performed on the joint covariance matrix to extract the principal eigenvector corresponding to the largest eigenvalue. The first half of this vector corresponds to the vibration channel of each node, and the second half corresponds to the pressure channel of each node. The sum of the squares of the vibration channel component amplitude and the pressure channel component amplitude corresponding to each node is taken as the synthetic response intensity of that node. Since the joint covariance matrix contains the second-order statistical correlation structure of the vibration and pressure fields under the same leakage source excitation, the principal eigenvector reflects the dominant spatial propagation mode of these two types of physical fields on the sensor array. The covariant information of vibration and pressure is compressed into the synchronously resolved synthetic response intensity. The node with the largest synthetic response intensity can be directly identified as the nearest node to the leakage source. Then, using the spatial coordinates of this node as a reference, the weighted sum of the synthetic response intensities of its adjacent nodes and the relative spatial direction vector is obtained to obtain the azimuth bias of the leakage source. This method of using vibration and pressure combined spatial spectrum decomposition for localization eliminates the dependence of traditional time-of-flight localization methods on array geometry and wavefront flatness assumptions. It can provide the spatial deviation angle information of the leakage source relative to the neighboring nodes in complex structures such as pipe bends and branches, significantly improving the spatial resolution and directional output capability of leakage points within the pipeline network. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0019] Figure 1 This is a flowchart of a method for monitoring gas leaks in underground pipeline networks; Figure 2 This is a flowchart of the mutual information feature extraction process for pipeline leak monitoring signals; Figure 3 This is a flowchart of a method for locating leakage sound sources based on the combined analysis of pipe wall vibration and pressure pulsation. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] See Figure 1 This invention provides a method for monitoring gas leaks in underground pipeline networks, comprising: deploying an acoustic wave sensor array along preset monitoring nodes of the underground pipeline network to simultaneously collect gas pipeline wall vibration signals and pipeline medium pressure pulsation signals at each monitoring node; constructing a time-domain response matrix with monitoring node numbers as row indices and sampling time points as column indices based on the pipeline wall vibration signals and pipeline medium pressure pulsation signals at each monitoring node; performing intrinsic mode decomposition on the time-domain response matrix to extract the intrinsic mode component sequences of each monitoring node, and calculating the mutual information entropy between the intrinsic mode component sequences of each monitoring node and the intrinsic mode component sequences of adjacent monitoring nodes; when the mutual information entropy is calculated... When the instantaneous rate of change of information entropy exceeds a preset threshold, the pipe wall vibration signal segment and the pipe medium pressure pulsation signal segment within a preset time window are extracted from the monitoring node and its adjacent monitoring nodes, centered on the monitoring node. The pipe wall vibration signal segment and the pipe medium pressure pulsation signal segment are arranged according to the monitoring node number to generate a first signal matrix and a second signal matrix. The joint covariance matrix of the first signal matrix and the second signal matrix is ​​calculated. The joint covariance matrix is ​​subjected to eigenvalue decomposition to extract the principal eigenvector corresponding to the largest eigenvalue. Based on the component amplitude distribution of the principal eigenvector at each monitoring node, the relative position of the leakage sound source is determined.

[0022] Example 1: In practical implementation, acoustic sensor arrays are deployed along pre-set monitoring nodes of the underground pipeline network to synchronously collect gas pipeline wall vibration signals and internal medium pressure pulsation signals at each monitoring node. This is achieved by installing sensor groups at each pre-set monitoring node and realizing synchronous acquisition across all nodes. The pre-set monitoring nodes are sensor placement locations pre-selected based on the pipeline route, burial depth, pipe diameter changes, and valve well locations. The pipeline length between adjacent monitoring nodes is determined as follows: Based on the pipeline material, diameter, wall thickness, medium pressure, and burial conditions, an attenuation model for the propagation of leakage acoustic waves along the pipe wall is established, with attenuation coefficients... according to Calculation, where The scattering attenuation coefficient is... The absorption attenuation coefficient, This is a comprehensive correction factor that includes dispersion and leakage losses from adjacent media; the minimum detectable signal amplitude is determined based on the signal-to-noise ratio parameters of the piezoelectric vibration sensor and pressure transmitter. ; Set the spacing between adjacent monitoring nodes satisfy: ,in The distance of the leaked sound source The amplitude of the signal after that, for The signal-to-noise ratio margin coefficient; under the premise of satisfying the above inequalities, select the maximum value that minimizes the total number of monitoring nodes. The value is used as a fixed spacing. As a preferred implementation, when precise attenuation model parameters are lacking, the spacing between adjacent monitoring nodes is 50 meters; when the pipe diameter is greater than DN300 and the medium pressure is higher than 1.6MPa, the spacing can be increased to 100 meters; when the pipe diameter is less than DN100 and the medium pressure is lower than 0.4MPa, the spacing should be reduced to 30 meters.

[0023] For each preset monitoring node, the piezoelectric vibration sensor is fixed to the outer wall of the pipe using mechanical clamps and a coupling layer. The mechanical clamps use stainless steel cable ties or bow-shaped clamps to press the piezoelectric vibration sensor base firmly against the metal surface of the pipe's outer wall. Silicone grease coupling medium is filled between the contact surface of the piezoelectric vibration sensor base and the pipe's outer wall to eliminate air gaps between the sensor and the pipe wall and to achieve acoustic impedance matching. The sensitive axis of the piezoelectric vibration sensor is parallel to the pipe wall surface to maximize the reception of guided wave signals propagating along the pipe wall. The piezoelectric vibration sensor incorporates a pre-charge amplifier, converting the charge signal generated by the piezoelectric ceramic under mechanical vibration into a low-impedance voltage signal output.

[0024] At the same monitoring node, the pressure transmitter is connected to the internal medium of the pipeline via a pressure tapping pipe. One end of the pressure tapping pipe is welded or threaded to a pre-installed pressure tapping port on the pipeline. This port is pre-installed on the pipeline sidewall during pipeline construction and connected to the pressure tapping pipe via a shut-off valve. The other end of the pressure tapping pipe is connected to the pressure transmitter's diaphragm capsule, which is filled with a separating fluid that transmits the pressure of the internal medium to the silicon piezoresistive or capacitive sensing element inside the pressure transmitter. The pressure transmitter outputs a voltage or current signal that is linearly related to the pressure pulsation of the internal medium in the pipeline.

[0025] Each monitoring node is equipped with a data acquisition unit, which performs synchronous analog-to-digital conversion on the voltage signal output by the piezoelectric vibration sensor and the electrical signal output by the pressure transmitter. The data acquisition unit contains a multi-channel synchronous sampling analog-to-digital converter, and all analog input channels share the same sampling clock and trigger signal to ensure strict time alignment between the pipe wall vibration signal and the pressure pulsation signal of the medium in the pipe within the same monitoring node.

[0026] To achieve synchronized data acquisition across all monitoring nodes, a common global clock source is used as the synchronization reference. This global clock source is generated by a satellite timing receiver module deployed at the monitoring center. This module receives timing signals from the Global Navigation Satellite System and outputs a standard time code and a second pulse signal. The standard time code and second pulse signal are distributed to the data acquisition unit of each monitoring node via a fiber optic time synchronization network. Specifically, a time server is deployed at the monitoring center. The time server locks onto the second pulse signal and standard time code output by the satellite timing receiver module and transmits the time synchronization message and synchronization pulse to the data acquisition unit of each monitoring node via a fiber optic Ethernet switch supporting the IEEE 1588 precision time protocol. Each monitoring node's data acquisition unit has a built-in slave clock module supporting the IEEE 1588 protocol. The slave clock module synchronizes its local clock with the global clock source in phase and frequency based on the received synchronization message and hardware timestamp, ensuring that the time deviation between the local clock and the global clock source is controlled within the sub-microsecond range.

[0027] After the synchronization clock is established, the monitoring center broadcasts a data acquisition start command to the data acquisition units of all monitoring nodes. Upon receiving the start command, each data acquisition unit simultaneously initiates analog-to-digital conversion at the rising edge of the next second pulse, continuously sampling at the same sampling frequency. The sampling frequency is determined based on the frequency range of the pipe wall vibration signal and the pressure pulsation signal of the medium inside the pipe, and is set to a value that satisfies the Nyquist theorem, for example, 5000 Hz. Each data acquisition unit buffers the sampled pipe wall vibration signal data and the pressure pulsation signal data of the medium inside the pipe according to the sampling point sequence, and appends the absolute timestamp corresponding to each sampling point to the data. The absolute timestamp is provided by the clock module and records the instantaneous absolute moment of sampling. In this way, the pipe wall vibration signal and the pressure pulsation signal of the medium inside the pipe of all monitoring nodes are given an absolute time mark based on the same time base, thus having a unified sampling time axis, providing multi-node signal data with precise time alignment for the subsequent construction of the time domain response matrix.

[0028] Example 2: In specific implementation, please refer to Figure 2 When constructing the time-domain response matrix based on the pipe wall vibration signals and pipe internal medium pressure pulsation signals from each monitoring node, the pipe wall vibration signals and pipe internal medium pressure pulsation signals for each monitoring node are arranged into row vectors according to their sampling time points. At each monitoring node, the data acquisition unit collects the pipe wall vibration signal sequence and the pipe internal medium pressure pulsation signal sequence within an observation period of N sampling points, following a unified sampling time axis. For monitoring node m, the pipe wall vibration signals from all N sampling points collected within the observation period are arranged in chronological order, forming a row vector containing N elements. The nth element in the row vector represents the amplitude of the pipe wall vibration signal at monitoring node m at the nth sampling time. Similarly, the pipe internal medium pressure pulsation signals from all N sampling points collected by monitoring node m within the same observation period are arranged in chronological order, forming another row vector containing N elements. The nth element in the row vector represents the amplitude of the pipe internal medium pressure pulsation signal at monitoring node m at the nth sampling time. The time sequence of the sampling points within each row vector strictly corresponds to the unified sampling time axis.

[0029] After obtaining the row vectors of all monitoring nodes, the row vectors of the pipe wall vibration signals of all M monitoring nodes are stacked vertically from 1 to M according to the monitoring node numbers to form a time-domain vibration response matrix. The time-domain vibration response matrix is ​​an M-row, N-column real number matrix, where the m-th row corresponds to the pipe wall vibration signal row vector of monitoring node m, and the n-th column corresponds to the pipe wall vibration signal value of all monitoring nodes at the n-th sampling time. Similarly, the row vectors of the pipe internal medium pressure pulsation signals of all M monitoring nodes are stacked vertically from 1 to M according to the monitoring node numbers to form a time-domain pressure response matrix. The time-domain pressure response matrix is ​​also an M-row, N-column real number matrix, where the m-th row corresponds to the pipe internal medium pressure pulsation signal row vector of monitoring node m, and the n-th column corresponds to the pipe internal medium pressure pulsation signal value of all monitoring nodes at the n-th sampling time. The column indices of both the time-domain vibration response matrix and the time-domain pressure response matrix are the sampling time point numbers, and the row indices are the monitoring node numbers; both have the same row and column dimensions.

[0030] When performing intrinsic mode decomposition (IMD) on the time-domain response matrix and extracting the IMD component sequences of each monitoring node, empirical mode decomposition (EMD) is performed on the signal sequence of each row in the time-domain vibration response matrix. EMD decomposition decomposes a row of the original signal sequence into a set of IMD components and a residual term using a recursive selection process. The recursive selection process involves: identifying all local maxima and local minima of the current signal sequence to be decomposed; using each local maxima and local minima as interpolation nodes; constructing a smooth upper envelope passing through all local maxima and a smooth lower envelope passing through all local minima using cubic spline interpolation; calculating the mean curves of the upper and lower envelopes; subtracting the mean curves from the current signal sequence to be decomposed to obtain an intermediate signal sequence; and checking whether the intermediate signal sequence satisfies the IMD conditions. The IMD conditions include that the number of zeros is equal to or differs from the number of extrema by no more than one, and that the upper envelope obtained by interpolation of local maxima and the lower envelope obtained by interpolation of local minima are symmetrical about the zero line. If the intermediate signal sequence does not satisfy the intrinsic mode function (IMF) condition, the intermediate signal sequence is treated as a new signal sequence to be decomposed, and the above recursive screening is repeated. If the intermediate signal sequence satisfies the IMF condition, the intermediate signal sequence is treated as an IMF component, and this IMF component is subtracted from the current signal sequence to be decomposed to obtain the remaining signal sequence. The remaining signal sequence is then treated as a new signal sequence to be decomposed, and the recursive screening is repeated until all IMF components and the final residual terms that cannot be further decomposed are extracted. For the signal sequence in each row of the time-domain pressure response matrix, the same empirical mode decomposition steps as those for the time-domain vibration response signal sequence are used to obtain multiple IMF components of the pressure pulsation signal in the pipeline at each monitoring node.

[0031] For each monitoring node, the first three intrinsic mode function (IMF) components obtained from empirical mode decomposition are selected and superimposed to form a reconstructed component sequence. The IMF components are arranged in the order they were extracted, with the first extracted IMF component being the highest frequency component, the second being the second highest frequency component, and the third being the third highest frequency component. Leakage-induced signals from pipe wall vibration and medium pressure pulsation are concentrated in the high-frequency region in the frequency domain. The first three IMF components contain the high-frequency oscillation modes in the signal; superimposing these three components can preserve leakage characteristic information and suppress low-frequency flow noise and environmental interference. The values ​​of the first three IMF components at each sampling time are summed to obtain the reconstructed component sequence, which has a length of N sampling points.

[0032] When calculating the mutual information entropy between the reconstructed component sequences of each monitoring node and the reconstructed component sequences of adjacent monitoring nodes, adjacent monitoring nodes are defined as nodes with adjacent monitoring node numbers, that is, upstream and downstream nodes that are immediately adjacent to the monitoring node on the pipeline path. For a monitoring node numbered m, its reconstructed component sequence is denoted as... ,in This represents the nth sampling time, where n ranges from 1 to N; the reconstructed component sequences of its adjacent monitoring nodes numbered m-1 and m+1 are respectively denoted as... and The calculation of mutual information entropy is based on the mutual information values ​​of the reconstructed component sequences within a preset time delay range. The preset time delay range is determined by the minimum delay index. and maximum delay number A limited set of integers, where -20, The value is 20. Assume the distance between adjacent monitoring nodes in the pipeline network is 50 meters, the sound velocity of the gas in the pipeline is approximately 430 m / s, the vibration propagation speed in the pipe wall is approximately 1200 m / s, and the maximum time difference between the arrival of the leak sound source at adjacent nodes does not exceed 50 / 430≈0.116 seconds, corresponding to approximately 580 sampling points at a sampling frequency of 5000 Hz. However, the leak-related components reflected in the reconstructed component sequences exhibit strong similarity between different nodes. Setting a smaller time delay range can reduce the computational load and avoid the introduction of spurious correlations. Therefore, the preset time delay range is set to an interval of -20 to 20 sampling points, corresponding to a time delay window of ±4 milliseconds, which can cover the effects of slight asynchrony between adjacent nodes and changes in the propagation speed of the medium.

[0033] When calculating the mutual information value between monitoring node m and its neighboring monitoring node p, p takes the value of m-1 or m+1. For each integer delay d within the preset time delay range, the... Sequence and When aligning sequences in time, As a reference sequence, with The sequence is a delayed sequence. The mutual information between the two reconstructed component sequences is calculated using the following formula:

[0034] in, Represents the reconstructed component sequence With reconstructed component sequences Mutual information value under delay d Represents the reconstructed component sequence The moderate amplitude value falls within the range The representative value within, Represents the reconstructed component sequence The moderate amplitude value falls within the range The representative value within, and These represent the reconstructed component sequences respectively. and The amplitude range is divided into equal-width bins to obtain a set of discrete amplitude values. The number of bins is uniformly set to 32. The minimum to maximum value of the signal amplitude is divided into 32 intervals. The amplitude value of all sample points in each interval is represented by the center value of the interval, a or b. Represents the reconstructed component sequence The marginal probability distribution of the value 'a' representing the moderate amplitude is derived from... The proportion of sample points within the range to the total number of sample points is estimated. Represents the reconstructed component sequence The marginal probability distribution of the occurrence of the moderate amplitude representative value b is given by The proportion of sample points within the range to the total number of sample points is estimated. This indicates the reconstructed component sequence under delay d. The intermediate amplitude representative value 'a' and the reconstructed component sequence The joint probability distribution of the simultaneous occurrence of the intermediate amplitude representation value b is calculated by aligning the two sequences in time. It falls within the interval of the amplitude representation value a and The proportion of sample point pairs falling within the interval of the magnitude representative value b is estimated to be the total number of valid aligned sample point pairs. The base of the logarithmic function is taken as the natural constant e.

[0035] For the preset time delay range For each integer delay d within the range, calculate the mutual information value. A set of mutual information values ​​is obtained, and the maximum value in this set is taken as the mutual information entropy between monitoring node m and its neighboring monitoring node p. For an internal node with two neighboring monitoring nodes, two mutual information entropies are calculated with its upstream and downstream neighboring nodes, respectively, and the larger of the two mutual information entropies is taken as the mutual information entropy of monitoring node m. For a monitoring node at the end of the pipeline with only one neighboring monitoring node, the mutual information entropy calculated with the only neighboring monitoring node is taken as the mutual information entropy of that monitoring node. Finally, each monitoring node obtains a mutual information entropy value, reflecting the strength of the signal correlation between that node and its neighboring nodes within the leakage characteristic frequency band.

[0036] After the above calculation process, for each monitoring node, two sets of mutual information entropy results will be obtained. One set is calculated from the reconstructed component sequence of the pipe wall vibration signal, and is denoted as the vibration channel mutual information entropy. The other set is calculated from the reconstructed component sequence of the pressure pulsation signal of the medium inside the pipe, and is denoted as the pressure channel mutual information entropy.

[0037] For monitoring nodes within the pipeline network with two adjacent upstream and downstream nodes, the mutual information entropy of the vibration channel is taken as the larger of the two mutual information entropies calculated from the vibration reconstruction component sequences of the node and its upstream and downstream adjacent nodes; the mutual information entropy of the pressure channel is taken as the larger of the two mutual information entropies calculated from the pressure reconstruction component sequences of the node and its upstream and downstream adjacent nodes. For monitoring nodes at the end of the pipeline network with only a single adjacent node, the mutual information entropy of the vibration channel and the mutual information entropy of the pressure channel are directly taken as the mutual information entropy results calculated from the channel reconstruction component sequences corresponding to the node and its only adjacent node.

[0038] To achieve the fusion of two physical quantities at the detection level, a fusion mutual information entropy of monitoring nodes is constructed. The calculation expression for the fusion mutual information entropy is as follows:

[0039] Representative number is The fusion and mutual information entropy of monitoring nodes; Representative number is The mutual information entropy of the vibration channels of the monitoring nodes; Representative number is The mutual information entropy of the pressure channels of the monitoring nodes; Represents the vibration channel weighting coefficient; This represents the pressure channel weighting coefficient.

[0040] The weighting coefficient is configured according to the operating conditions of the pipeline network. When the external environmental vibration interference of the pipeline network is strong, the weighting coefficient is increased. The value of is reduced The value of is determined by the degree of turbulence within the pipe network; when the internal fluid turbulence is severe, the value of is increased. The value of is reduced The value of . Under normal pipeline network operating conditions without special disturbances, and Take the equal values.

[0041] This fused mutual information entropy serves as the node mutual information entropy used for subsequent trigger judgment, achieving mathematical fusion of vibration and pressure physical quantities at the detection trigger level. It eliminates the need to independently use vibration channel mutual information entropy or pressure channel mutual information entropy for leakage trigger judgment.

[0042] Example 3: In the specific implementation, after obtaining the mutual information entropy of each monitoring node, the mutual information entropy of each monitoring node is arranged into an entropy value time series in chronological order. Let the set of monitoring node numbers be {1,2,…,M}. For each monitoring node, when it updates and calculates the mutual information entropy each time, the timestamp of the calculation time and the corresponding mutual information entropy value are recorded. The time interval between two adjacent mutual information entropy calculations is kept at a fixed step size ΔT, which is determined according to the update time window of the time domain response matrix, and is set to 0.1 seconds. When the mutual information entropy is calculated continuously with a period of ΔT, the mutual information entropy value obtained by monitoring node m at time t=t0+k·ΔT is recorded as Hm[k], where k is the time sequence number, and k increases sequentially from 0. According to the time sequence number k from 0 to K, the mutual information entropy values ​​of monitoring node m are arranged into a one-dimensional sequence Hm[0],Hm[1],…,Hm[K], which is the entropy value time series of monitoring node m.

[0043] In practice, the entropy time series of each monitoring node is calculated using a moving difference method to obtain the instantaneous rate of change of entropy. The moving difference calculation employs a five-point difference scheme, whereby for each time index k≥4 in the entropy time series, the instantaneous rate of change of entropy is calculated using a weighted difference method based on the mutual information entropy at that moment and the mutual information entropies of the four preceding times. The specific formula for calculating the instantaneous rate of change of entropy is as follows:

[0044] in, This represents the instantaneous rate of change of the entropy value of monitoring node m at time index k. It is the time interval between adjacent mutual information entropy values ​​in an entropy time series, in seconds. This represents the mutual information entropy value of monitoring node m at time index k. This represents the mutual information entropy value of monitoring node m at time number k-1. This represents the mutual information entropy value of monitoring node m at time number k-2. This represents the mutual information entropy value of monitoring node m at time number k-3. This represents the mutual information entropy value of monitoring node m at time number k-4. Weighting coefficients. These are the five-point backward difference coefficients, and their values ​​are as follows: , , , , This set of weighting coefficients is based on numerical differential theory, using the function values ​​at five equidistant nodes to perform a backward difference approximation of the first derivative, by setting the node spacing to... The mutual information entropy at five time points is expanded using a Taylor series and low-order error terms other than the first derivative term are eliminated. This method can approximate the first-order time derivative with fourth-order precision, thereby accurately capturing the instantaneous changes in mutual information entropy.

[0045] In practice, the instantaneous rate of change of entropy at each monitoring node is compared with a preset threshold to determine whether a gas leak has caused a sudden change in mutual information entropy. The preset threshold is a constant pre-set based on the range of mutual information entropy fluctuations in the pipeline network under normal operating conditions. Specifically, during a stable operating period without pipeline leaks, the time series of mutual information entropy at each monitoring node is continuously collected, and the root mean square value of the instantaneous rate of change of entropy at each monitoring node during the stable period is calculated. The obtained root mean square value is multiplied by a coefficient of 6 to obtain a preset threshold, i.e. The coefficient 6 was chosen based on the following: Under the assumption of a normal distribution, the false alarm probability corresponding to 6 times the root mean square value is approximately... (i.e., more than in a one-sided test) The sample point ratio is extremely low, which ensures sufficient detection sensitivity for entropy value changes caused by leaks while keeping the false alarm rate of occasional entropy value fluctuations caused by environmental noise, fluid turbulence, etc., under normal operating conditions as leak events to an extremely low level. If the coefficient value is too small (e.g., 34), normal fluctuations are likely to exceed the threshold, leading to frequent false alarms; if the coefficient value is too large (e.g., 810), entropy value changes caused by weak leaks may not trigger an alarm, resulting in missed detections. The root mean square (RMS) value is an optimal value for balancing detection sensitivity and false alarm resistance in engineering practice. Those skilled in the art can make adaptive adjustments based on the actual operating environment of the pipeline network. Five-point differential sliding calculations are performed on the time series of mutual information entropy in the vibration channel, pressure channel, and fused mutual information entropy, respectively, to obtain the instantaneous change rate of the vibration channel entropy, pressure channel entropy, and fused mutual information entropy. This method uses the instantaneous change rate of the fused mutual information entropy as the primary trigger criterion, while configuring the instantaneous change rates of the vibration channel entropy and pressure channel entropy as auxiliary verification conditions. The primary trigger criterion is that the instantaneous change rate of the fused mutual information entropy exceeds a preset primary threshold. When the primary trigger criterion is met, the auxiliary conditions are then verified: if at least one of the instantaneous change rates of the vibration channel entropy and pressure channel entropy exceeds the corresponding auxiliary threshold, a leakage anomaly is determined, and the corresponding monitoring node is marked as a trigger node. If only the instantaneous change rate of the vibration channel entropy exceeds the auxiliary threshold, but the instantaneous change rate of the fused mutual information entropy does not reach the main threshold, it is determined to be external vibration interference, and no leakage event is triggered. If only the instantaneous change rate of the pressure channel entropy exceeds the auxiliary threshold, but the instantaneous change rate of the fused mutual information entropy does not reach the main threshold, it is determined to be fluid disturbance within the pipe, and no leakage event is triggered. This triggering logic avoids the false alarm and missed alarm problems caused by single-channel threshold judgment. The change characteristics of the two physical quantities, vibration and pressure, are comprehensively measured by fused mutual information entropy, and then combined with single-channel auxiliary verification to complete event identification.

[0046] The instantaneous change rate of the entropy value of monitoring node m at time number k When the value exceeds a preset threshold, the monitoring node m is marked as a trigger node, and the time corresponding to the time sequence k is taken as the trigger time. Simultaneously with marking the trigger node, the forward time boundary of the trigger time is recorded, and the judgment of other monitoring nodes at the same time sequence is stopped; the signal segment extraction process then proceeds directly.

[0047] In practice, the vibration signal segments of the pipe wall and the pressure pulsation signal segments of the medium inside the pipe are extracted from the trigger node and its adjacent monitoring nodes within a preset time window. The preset time window length is... Number of sampling points. Window length. The selection follows these principles: The larger the value, the higher the frequency resolution (frequency resolution) , (Sampling frequency), but the time positioning accuracy decreases; The smaller the value, the higher the time positioning accuracy, but the lower the frequency resolution. To achieve both complete capture of the time-domain waveform of the leakage signal and effective resolution of its frequency-domain characteristics, a balance needs to be struck between time and frequency resolution. The main energy of the leakage signal waveform is typically concentrated within the first 400 milliseconds, a duration sufficient to contain the main wave information of pipe wall vibration and pressure pulsations caused by the leak; simultaneously, at a sampling frequency of 5000 Hz, The corresponding frequency resolution is 2.5 Hz ( This is sufficient to distinguish the main frequency components in the leaked signal. Therefore, take As a preferred value that balances temporal integrity and frequency domain resolution, those skilled in the art can adaptively adjust this value based on the signal characteristics of the actual pipeline network. Taking the triggering time of the triggering node as the endpoint (i.e., the sampling point number corresponding to the triggering time), L sampling points are extracted forward to obtain the pipe wall vibration signal segment and the pipeline medium pressure pulsation signal segment of the triggering node. Both signal segments are one-dimensional signal sequences of length L. Simultaneously, all adjacent monitoring nodes of the triggering node are determined. Adjacent monitoring nodes refer to monitoring nodes directly connected to the triggering node via pipe segments in the pipeline network topology. For each adjacent monitoring node, the pipe wall vibration signal segment and the pipeline medium pressure pulsation signal segment within the same preset time window range as the triggering node are taken on the same unified sampling time axis. That is, the signal segment within the time window range of L sampling points forward from the triggering time is also taken as the pipe wall vibration signal segment and the pipeline medium pressure pulsation signal segment of the adjacent monitoring node. All extracted signal segments retain absolute timestamp information to verify the consistency of the time windows between nodes. The above method yields a set of signal segments consisting of the trigger node and all its adjacent monitoring nodes, which are then used to construct the first signal matrix and the second signal matrix.

[0048] Example 4: In specific implementation, please refer to Figure 3The pipe wall vibration signal segments from all the monitoring nodes are arranged from top to bottom according to their monitoring node numbers to form the first signal matrix. The set of monitoring nodes participating in the signal segment extraction includes the trigger node and all its adjacent monitoring nodes. Let there be P monitoring nodes in the set, where P equals 1 plus the number of adjacent monitoring nodes of the trigger node. These P monitoring nodes are sorted in ascending order according to their monitoring node numbers. After sorting, the pipe wall vibration signal segment of each monitoring node is sequentially used as a row of the first signal matrix. The pipe wall vibration signal segment of each monitoring node is a one-dimensional signal sequence consisting of W sampling points collected within a preset time window, where W is 2000. The pipe wall vibration signal segments of the P monitoring nodes are stacked vertically in numerical order to form a P-row, W-column real number matrix, which is the first signal matrix.

[0049] In practice, the captured pipeline medium pressure pulsation signal segments from all monitoring nodes are arranged from top to bottom according to the monitoring node number to form a second signal matrix. Corresponding to the P monitoring nodes as arranged above, the pipeline medium pressure pulsation signal segment of each monitoring node is taken as a row of the second signal matrix. Each monitoring node's pipeline medium pressure pulsation signal segment is also a one-dimensional signal sequence composed of W sampling points within the same preset time window. The pipeline medium pressure pulsation signal segments of the P monitoring nodes are stacked vertically in numerical order to form a P-row, W-column real number matrix, which is the second signal matrix. The row indices of both the first and second signal matrices correspond to the monitoring node numbering, and the column indices correspond to the sampling point sequence numbers within the preset time window.

[0050] In the specific implementation, the first signal matrix and the second signal matrix are subjected to mean-reduction processing respectively. For the first signal matrix, the average value of the amplitudes of all W sampling points in each row of the first signal matrix is ​​calculated. The average value of the row is then subtracted from the amplitude of each sampling point in that row to obtain the mean-reduced first signal matrix. The mean of each row of the mean-reduced first signal matrix is ​​zero. For the second signal matrix, the same method is used: the average value of the amplitudes of all W sampling points in each row of the second signal matrix is ​​calculated. The average value of the row is then subtracted from the amplitude of each sampling point in that row to obtain the mean-reduced second signal matrix.

[0051] In practical implementation, the first signal matrix after mean-reduction processing and the second signal matrix after mean-reduction processing are vertically concatenated to form a joint signal matrix. The vertical concatenation method is as follows: all rows of the first signal matrix after mean-reduction processing are kept in the upper half of the joint signal matrix, and all rows of the second signal matrix after mean-reduction processing are kept in the lower half of the joint signal matrix. The joint signal matrix has 2P rows and W columns. The first P rows of the joint signal matrix correspond to the mean-reduction pipe wall vibration signals of P monitoring nodes, and the last P rows correspond to the mean-reduction pipe internal medium pressure pulsation signals of P monitoring nodes. The construction of the joint covariance matrix relies on the signal segments obtained after the aforementioned detection layer completes the dual physical quantity feature verification. The fusion mutual information entropy of the detection layer is only used for event trigger identification and does not participate in the numerical calculation of the joint covariance matrix. The joint covariance matrix fully preserves the second-order statistical relationship of the original vibration signal segments and the original pressure signal segments, including both the signal variance of each channel and the cross-channel covariance term between the vibration channel and the pressure channel. The cross-channel covariance term characterizes the statistical correlation between the pipe wall vibration response and the internal pressure pulsation response under the same leakage excitation. This correlation cannot be represented by the fused mutual information entropy of the detection layer. The fused mutual information entropy of the detection layer determines whether an event has occurred, while the joint covariance matrix directly uses the original dual-physical signal fragments to analyze the sound source location after the event is triggered. The two layers have different functions. The detection layer completes the logical judgment of whether the event has been triggered, while the localization layer uses the complete statistical information of the original dual-physical signals to solve for the location of the leakage sound source.

[0052] In practical implementation, the joint covariance matrix is ​​calculated using the joint signal matrix. Let the joint signal matrix be denoted as... , The transpose of the matrix is ​​denoted as Joint covariance matrix The calculation formula is:

[0053] in, Denotes the joint covariance matrix. The dimension is 2P rows and 2P columns; This indicates the number of sampling points within a preset time window, i.e., the number of time points in the signal segment. The value is 2000; This represents the joint signal matrix after mean removal, with dimensions of 2P rows and W columns; Represents the joint signal matrix The transpose of the matrix has dimensions W rows and 2P columns. Representation matrix and The product of these two elements results in a 2P x 2P square matrix, where each element represents the covariance of the corresponding two rows of signals in the joint signal matrix. times; divided by back, The element in the i-th row and j-th column of the joint signal matrix is ​​the estimated covariance of the i-th row signal and the j-th row signal in the sampling point dimension.

[0054] In practice, the joint covariance matrix is ​​subjected to power-law iterative eigendecomposition to extract the principal eigenvector corresponding to the largest eigenvalue. The specific process of power-law iterative eigendecomposition is as follows: First, a non-zero initial vector of dimension 2P is randomly initialized. This initial vector is then compared with the joint covariance matrix... Multiply the vectors to obtain the updated vector; calculate the Euclidean norm of the updated vector, and normalize it by dividing the updated vector by its norm to obtain the iterative vector; in the next iteration, multiply the iterative vector obtained in the previous iteration with the joint covariance matrix again and normalize it; repeat the multiplication and normalization operations until the cosine similarity between the iterative vectors obtained from two adjacent iterations is greater than 0.9999, at which point the iteration stops, and the iterative vector obtained in the last iteration is taken as the principal feature vector. The principal feature vector is a real column vector of dimension 2P, where the first P components correspond to the monitoring node channels of the first signal matrix after mean removal processing, and the last P components correspond to the monitoring node channels of the second signal matrix after mean removal processing.

[0055] In practice, the sum of the squares of the amplitudes of the first half-dimensional component and the second half-dimensional component corresponding to each monitoring node in the main feature vector is taken to obtain the synthetic response intensity of that monitoring node. For the r-th monitoring node after sorting, where r ranges from 1 to P, its corresponding first half-dimensional component amplitude is the value of the r-th element in the main feature vector, denoted as [equation missing]. The magnitude of the second-half-dimensional component is the value of the (P+r)th element in the principal eigenvector, denoted as... The synthetic response intensity of monitoring node r The calculation method is as follows Calculate the synthetic response intensity for each of the P monitoring nodes, and determine the monitoring node with the largest synthetic response intensity value as the nearest neighbor node to the leakage sound source.

[0056] In practice, after identifying the nearest node to the leakage sound source, the azimuth of the leakage sound source is determined based on the spatial coordinates of this nearest node in the pipeline network and the ratio of the synthesized response intensity of its adjacent nodes. The spatial coordinates of the nearest node in the pipeline network are denoted as follows: The spatial coordinates of each neighboring node of the nearest node are denoted as . The combined response intensity of adjacent nodes is denoted as The orientation bias is determined by taking the nearest neighbor node as the origin and calculating the combined response intensity of its neighboring nodes. The spatial direction vector of the adjacent node relative to its nearest neighbor is weighted and summed. The composite vector obtained by weighted summation is represented as follows: The direction of the synthesized vector is the azimuth offset of the leakage sound source, that is, the spatial direction from the nearest node to the direction of the synthesized vector is the deviation of the leakage sound source from the location of the nearest node.

[0057] Example 5: In practice, the underground pipeline gas leak monitoring system includes at least one monitoring server, multiple field data acquisition terminals, a time synchronization network, and a data communication network connecting the monitoring server and the field data acquisition terminals. The monitoring server, as the core computing and control unit of the system, contains a memory, a processor, and a computer program stored in the memory and running on the processor. Field data acquisition terminals are deployed at each pre-set monitoring node, and each field data acquisition terminal integrates a piezoelectric vibration sensor, a pressure transmitter, and a multi-channel synchronous acquisition module. The time synchronization network provides a unified time reference for all field data acquisition terminals, and the data communication network transmits the signals collected by each field data acquisition terminal to the monitoring server in real time.

[0058] The memory is configured as a combination of volatile random access memory (RAM) and non-volatile flash memory. The volatile RAM temporarily stores intermediate data generated during processor operation, including the time-domain response matrix, intrinsic mode function components, mutual information entropy sequence, signal segment matrix, and joint covariance matrix. The non-volatile flash memory persistently stores the computer program itself, historical monitoring data, preset thresholds and preset time window lengths, as well as the number, spatial coordinates, and pipeline topology connections of each monitoring node. The processor employs a multi-core central processing unit, supporting parallel data stream processing and multi-threaded task scheduling. The processor connects to the memory and communication interface via the system bus. The data communication network uses a ring topology composed of industrial Ethernet switches. Each field data acquisition terminal connects to the switch via an Ethernet interface. The switch connects to the monitoring server's Ethernet port via a fiber optic uplink port. The monitoring server receives pipeline wall vibration signal data and pipeline medium pressure pulsation signal data from all field data acquisition terminals via the Ethernet port.

[0059] The time synchronization network consists of a satellite time receiver module located in the monitoring server room, a time server connected to the satellite time receiver module, fiber optic Ethernet switches supporting the precise time protocol, and clock modules inside the field data acquisition terminals. The satellite time receiver module continuously receives time signals from the Global Navigation Satellite System (GNSS) and outputs a second pulse signal and a standard time code containing Coordinated Universal Time (UTC). The time server locks the second pulse signal, aligns its local clock with UTC, and sends time synchronization messages conforming to the IEEE 1588 precise time protocol to each switch via the fiber optic Ethernet interface. Each switch performs transparent clock processing on the precise time protocol messages and forwards them to each field data acquisition terminal. The clock module inside the field data acquisition terminal receives the precise time protocol messages, parses the master clock timestamp and link delay from the messages, and adjusts the local clock phase and frequency to synchronize the local clock with the global clock source.

[0060] The computer program running in the processor adopts a modular architecture design, including a signal acquisition synchronization module, a time-domain response matrix construction module, an intrinsic mode decomposition and mutual information entropy calculation module, an instantaneous rate of change detection and signal interception module, a joint covariance matrix and principal eigenvector analysis module, and a leakage sound source localization and transmission module. The signal acquisition synchronization module is responsible for issuing synchronization acquisition commands to each field data acquisition terminal. These commands include a unified sampling frequency and the absolute time of acquisition initiation. The field data acquisition terminals synchronously start sampling at this absolute time and generate signal sampling points with absolute timestamps at a unified sampling frequency. The signal acquisition synchronization module continuously receives pipeline wall vibration signal data and pipeline medium pressure pulsation signal data uploaded by each field data acquisition terminal, aligns them with timestamps, and stores them in a circular buffer of the memory.

[0061] The time-domain response matrix construction module retrieves data from all monitoring nodes within the most recently acquired observation period from the memory buffer. After frame alignment of the data within that observation period, it constructs a time-domain vibration response matrix and a time-domain pressure response matrix according to the monitoring node numbers. In the time-domain vibration response matrix, each row represents the pipeline wall vibration signal vector of a single monitoring node, and each column represents the vibration amplitude of all monitoring nodes at the same sampling time. Similarly, in the time-domain pressure response matrix, each row represents the pipeline medium pressure pulsation signal vector of a single monitoring node, and each column represents the pressure pulsation amplitude of all monitoring nodes at the same sampling time. After the construction of both matrices is complete, the matrix data is written to memory for the next module to read.

[0062] The intrinsic mode decomposition (IMD) and mutual information entropy calculation module reads the time-domain vibration response matrix and the time-domain pressure response matrix, and performs empirical mode decomposition on each row of the signal sequence in both matrices. The IMD uses a recursive selection process, implemented as an iterative subroutine in the computer program. The subroutine identifies all extreme points of the signal sequence, constructs upper and lower envelopes, calculates the envelope mean, extracts intermediate sequences, checks whether the intermediate sequences satisfy the intrinsic mode function (IMF) conditions, and outputs them as IMF components if they do. The remaining sequences are then decomposed further until all IMF components are extracted. For each monitoring node, the module retains the first three IMF components and accumulates them according to sampling points to form a reconstructed component sequence. Based on the pipeline topology table stored in memory, the module determines the adjacent monitoring nodes for each monitoring node and calculates the mutual information value between the reconstructed component sequence of that monitoring node and the reconstructed component sequences of its adjacent monitoring nodes within a preset time delay range. When calculating mutual information, the reconstructed component sequences are first discretized into equal-width bins to obtain marginal probability distributions. Then, the joint probability distributions under different delays are calculated, and the mutual information values ​​corresponding to each delay are obtained by substituting them into the mutual information formula. The maximum value is taken as the mutual information entropy. The mutual information entropies of all monitoring nodes are organized into a time series and stored in memory.

[0063] The instantaneous change rate detection and signal interception module periodically reads the entropy value time series of each monitoring node from the memory and calculates the instantaneous change rate of entropy value using a five-point backward difference formula. The module compares the calculated instantaneous change rate of entropy value with a preset threshold node by node. The preset threshold is pre-stored in non-volatile memory and is obtained by multiplying the root mean square value of the instantaneous change rate of entropy value obtained during a leak-free quiet period in the pipeline network by 6. When the instantaneous change rate of entropy value of any monitoring node exceeds the preset threshold, the instantaneous change rate detection and signal interception module immediately marks that monitoring node as a trigger node, determines the trigger time, and intercepts a fixed-length time window. The module queries all adjacent monitoring nodes of the trigger node according to the stored topology relationship. For the trigger node and each adjacent monitoring node, it retrieves the pipe wall vibration signal segment and the pipe medium pressure pulsation signal segment within the corresponding time window from the memory, organizes them into a first signal matrix and a second signal matrix, and passes them to the next module.

[0064] The joint covariance matrix and principal eigenvector analysis module receives the first signal matrix and the second signal matrix, performs row-wise mean removal on both matrices, and then concatenates the mean-removed first signal matrix and the mean-removed second signal matrix vertically to form the joint signal matrix. The module calculates the joint covariance matrix using the formula: the product of the joint signal matrix and its transpose divided by the number of time points in the signal segment minus one. The module performs power-law iterative eigendecomposition on the joint covariance matrix, randomly initializing the iterative vector and repeatedly performing matrix multiplication and normalization until the iterative vector converges, obtaining the principal eigenvector corresponding to the largest eigenvalue. The module sums the squares of the first and second half-dimensional components corresponding to each monitoring node in the principal eigenvector to generate the synthetic response intensity of each monitoring node, and selects the monitoring node with the largest synthetic response intensity as the nearest neighbor node to the leakage source.

[0065] The leak sound source localization and transmission module uses the spatial coordinates of the nearest node and the synthetic response intensity of adjacent nodes as the origin. It then performs a weighted summation of the synthetic response intensity of adjacent nodes and the spatial direction vector pointing from adjacent nodes to the nearest node. The direction of the resulting synthetic vector is used as the azimuth bias of the leak sound source. The module packages the nearest node's number, spatial coordinates, azimuth bias, and corresponding timestamp into leak alarm information, which is then transmitted via the data communication network to the monitoring center's display screen for display. Simultaneously, it stores the information in the historical alarm record database in the memory.

[0066] When the processor executes the computer program, the above six modules run sequentially in a pipeline manner. After each round of full-process calculation is completed, a leak judgment result is output, and the next round of calculation is immediately started, realizing continuous online monitoring of gas leaks in underground pipeline networks.

[0067] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for monitoring gas leaks in underground pipeline networks, characterized in that, The method includes: An acoustic wave sensor array is deployed along the preset monitoring nodes of the underground pipeline network to simultaneously collect the vibration signals of the gas pipeline wall and the pressure pulsation signals of the medium inside the pipeline at each monitoring node. Based on the pipe wall vibration signal and the pressure pulsation signal of the medium inside the pipe at each monitoring node, a time-domain response matrix is ​​constructed with the monitoring node number as the row index and the sampling time point as the column index. The time-domain response matrix is ​​subjected to intrinsic mode decomposition to extract the intrinsic mode component sequence of each monitoring node, and the mutual information entropy between the intrinsic mode component sequence of each monitoring node and the intrinsic mode component sequence of adjacent monitoring nodes is calculated. When the instantaneous change rate of the mutual information entropy exceeds a preset threshold, the monitoring node and its adjacent monitoring nodes are captured within a preset time window as the center, along with segments of the pipe wall vibration signal and the pipe medium pressure pulsation signal. The pipe wall vibration signal segments and the pipe medium pressure pulsation signal segments are arranged according to the monitoring node number to generate a first signal matrix and a second signal matrix. The joint covariance matrix of the first signal matrix and the second signal matrix is ​​calculated. The joint covariance matrix is ​​subjected to eigenvalue decomposition to extract the principal eigenvector corresponding to the largest eigenvalue. Based on the component amplitude distribution of the principal eigenvector at each monitoring node, the relative position of the leakage sound source is determined.

2. The method for monitoring gas leaks in underground pipelines according to claim 1, characterized in that, The pre-set monitoring nodes along the underground pipeline network are equipped with acoustic wave sensor arrays to synchronously collect the gas pipeline wall vibration signals and the pressure pulsation signals of the medium inside the pipeline at each monitoring node, including: At each preset monitoring node, a piezoelectric vibration sensor is coupled and fixed to the outer wall of the pipe, and a pressure transmitter is connected to the internal medium of the pipe through a pressure tapping pipe. Using the same global clock source as the synchronization reference, the piezoelectric vibration sensors and pressure transmitters of each monitoring node are controlled to start collecting data at the same sampling frequency, so that the pipe wall vibration signal and the pressure pulsation signal of the medium in the pipe of each monitoring node have a unified sampling time axis.

3. The method for monitoring gas leaks in underground pipelines according to claim 1, characterized in that, The step of constructing a time-domain response matrix based on the pipe wall vibration signal and the pressure pulsation signal of the medium inside the pipe at each monitoring node, with the monitoring node number as the row index and the sampling time point as the column index, includes: The pipe wall vibration signal and the pressure pulsation signal of the medium inside the pipe at each monitoring node are arranged into row vectors according to the sampling time points; The pipeline wall vibration signal row vectors of all monitoring nodes are stacked vertically according to the monitoring node number to form a time-domain vibration response matrix. The pipeline internal medium pressure pulsation signal row vectors of all monitoring nodes are stacked vertically according to the monitoring node number to form a time-domain pressure response matrix.

4. The method for monitoring gas leaks in underground pipelines according to claim 3, characterized in that, The step of performing intrinsic mode decomposition on the time-domain response matrix, extracting the intrinsic mode component sequences of each monitoring node, and calculating the mutual information entropy between the intrinsic mode component sequences of each monitoring node and the intrinsic mode component sequences of adjacent monitoring nodes includes: Empirical mode decomposition is performed on the signal sequence of each row in the time-domain vibration response matrix and the time-domain pressure response matrix to obtain multiple intrinsic mode function components of each monitoring node; The first three intrinsic mode function components of each monitoring node are selected and superimposed to form a reconstructed component sequence. The mutual information value between the reconstructed component sequence of the monitoring node and the reconstructed component sequences of adjacent monitoring nodes within a preset time delay range is calculated, and the maximum value of the mutual information value is taken as the mutual information entropy.

5. The method for monitoring gas leaks in underground pipelines according to claim 1, characterized in that, When the instantaneous change rate of the mutual information entropy exceeds a preset threshold, taking the monitoring node as the center, the pipe wall vibration signal segment and the pipe medium pressure pulsation signal segment of the monitoring node and its adjacent monitoring nodes within a preset time window are extracted, including: The mutual information entropy of each monitoring node is arranged in chronological order to form an entropy value time series. The sliding difference is then used to calculate the instantaneous rate of change of the entropy value. When the instantaneous rate of change of the entropy value of a certain monitoring node exceeds a preset threshold, the monitoring node is marked as a trigger node. Taking the triggering time of the trigger node as the endpoint, pipe wall vibration signal segments and pipe medium pressure pulsation signal segments of a preset time window length are extracted. At the same time, pipe wall vibration signal segments and pipe medium pressure pulsation signal segments of all adjacent monitoring nodes of the trigger node within the same preset time window are extracted.

6. The method for monitoring gas leaks in underground pipelines according to claim 1, characterized in that, The step of arranging the pipe wall vibration signal segments and the pipe internal medium pressure pulsation signal segments according to the monitoring node number to generate a first signal matrix and a second signal matrix, and calculating the joint covariance matrix of the first signal matrix and the second signal matrix includes: The pipe wall vibration signal segments of all the monitoring nodes are arranged from top to bottom according to the monitoring node number to form the first signal matrix, and the pipe medium pressure pulsation signal segments of all the monitoring nodes are arranged from top to bottom according to the monitoring node number to form the second signal matrix. The first signal matrix and the second signal matrix are respectively subjected to mean removal processing. The mean-removed first signal matrix and the mean-removed second signal matrix are vertically concatenated to form a joint signal matrix. The joint signal matrix is ​​multiplied by its own transpose and divided by the number of time points of the signal segment minus one to obtain the joint covariance matrix.

7. The method for monitoring gas leaks in underground pipelines according to claim 1, characterized in that, The step of performing eigenvalue decomposition on the joint covariance matrix, extracting the principal eigenvector corresponding to the largest eigenvalue, and determining the relative location of the leakage sound source based on the component amplitude distribution of the principal eigenvector at each monitoring node includes: The joint covariance matrix is ​​subjected to power-law iterative eigenvalue decomposition to extract the principal eigenvector corresponding to the largest eigenvalue. The dimension of the principal eigenvector is twice the number of monitoring nodes. The first half of the dimension corresponds to the pipe wall vibration signal channel of each monitoring node, and the second half of the dimension corresponds to the pipe medium pressure pulsation signal channel of each monitoring node. The sum of the squares of the amplitudes of the first half-dimensional component and the second half-dimensional component corresponding to each monitoring node in the main feature vector is obtained to obtain the synthetic response intensity of the monitoring node. The monitoring node with the largest synthetic response intensity is taken as the nearest node of the leakage sound source. The spatial coordinates of the nearest node in the pipeline network are used as a reference, and the azimuth bias of the leakage sound source is determined by combining the ratio of the synthetic response intensity of the adjacent nodes of the nearest node.

8. The method for monitoring gas leaks in underground pipelines according to claim 5, characterized in that, The sliding difference calculation adopts a five-point difference format, and the instantaneous change rate of the entropy value is the ratio of the weighted difference mean of the mutual information entropy at the current moment and the mutual information entropy at the previous four moments to the time interval.

9. The method for monitoring gas leaks in underground pipelines according to claim 7, characterized in that, The orientation bias is determined as follows: taking the nearest node as the origin, the composite response intensity of each of its adjacent nodes is weighted and summed with the spatial direction vector of the adjacent node relative to the nearest node, and the direction of the resulting composite vector is the orientation bias of the leakage sound source.

10. An underground pipeline gas leak monitoring system, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the underground pipeline gas leakage monitoring method according to any one of claims 1 to 9.

Images