APSO-HMM-based pipeline leakage aperture identification method

Abstract

The invention relates to the technical field of pipeline leakage, in particular to a pipeline leakage aperture recognition method based on APSO-HMM. The method comprises the following steps: constructing a plurality of training sets, wherein each training set corresponds to one pipeline leakage aperture; inputting each training set into an HMM (Hidden Markov Model) for training iteration, and optimizing an initial value of an observed value in each state, namely an initial value of B, by applying an APSO (Advanced Persistent Swarm Optimization) algorithm in a training process; the optimal initial value obtained by solving through the APSO algorithm is substituted into the HMM model, iterative calculation is conducted on the HMM model through a Baum-Welch algorithm, and a trained APSO-HMM model is obtained; and extracting a feature vector T of the to-be-identified infrasonic wave original signal leakage signal x (t), respectively inputting the feature vector T into each trained APSO-HMM model, and outputting a pipeline leakage aperture corresponding to the APSO-HMM model with the maximum log-likelihood probability, namely an identification result. The method can avoid improper selection of the initial value of the model, thereby improving the recognition effect of the pipeline leakage aperture.

Description

Technical Field

[0001] The invention relates to the technical field of pipeline leakage, and in particular to a pipeline leakage aperture identification method based on APSO-HMM. Background Art

[0002] Most urban pipeline networks are distributed in a network shape, covering a wide range and numerous pipeline nodes. For tiny leakage signals of urban buried pipelines, the collected signals are mixed with a large amount of non-related signals such as environmental noise and medium noise. Tiny pipeline leaks, especially persistent pinhole leaks, are difficult to detect. Therefore, how to effectively detect and locate tiny leaks in the early stages of pipelines has good economic value and social significance.

[0003] The Hidden Markov Model (HMM) is a time series model and a memoryless non-stationary random process. It has a strong ability to process the statistical information and time series information contained in the sound signal, and can characterize the time-varying signal very well. According to the different signal characteristics of pipeline leakage, the Hidden Markov Model under different conditions is established to classify and identify pipeline leakage. It has a good recognition rate, especially for small continuous leakage, and has a good research effect.

[0004] The Baum-Welch algorithm is an important method for solving the parameters of the hidden Markov model, and is used when the hidden sequence cannot be obtained. However, the existing HMM model has a poor recognition effect on pipelines. Summary of the invention

[0005] The technical problem to be solved by the present invention is to overcome the defects of the prior art and provide a pipeline leakage aperture identification method based on APSO-HMM, which can avoid improper selection of the initial value of the model, thereby improving the recognition effect of the pipeline leakage aperture.

[0006] In order to solve the above technical problems, the technical solution of the present invention is: a pipeline leakage aperture identification method based on APSO-HMM, comprising:

[0007] Construct multiple training sets, each training set corresponds to a pipeline leakage aperture; wherein the construction process of each training set includes: obtaining multiple infrasound original leakage signals x(t) corresponding to the pipeline leakage aperture; extracting the feature vector T of each infrasound original leakage signal x(t), and the set formed by the feature vectors T of the multiple infrasound original leakage signals x(t) is used as the training set of the HMM model for the pipeline leakage aperture;

[0008] Each training set is input into an HMM model for training iteration. During the training process, the APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; the optimal initial value obtained by the APSO algorithm is brought into the HMM model, and the HMM model is iteratively calculated using the Baum-Welch algorithm to obtain a trained APSO-HMM model;

[0009] Extract the feature vector T of the original infrasound signal leakage x(t) to be identified, input the feature vector T into each trained APSO-HMM model, and output the pipeline leakage aperture corresponding to the APSO-HMM model with the largest log-likelihood probability, which is the identification result;

[0010] Among them, the steps of extracting the characteristic vector T of the original infrasound leakage signal x(t) are as follows: the original infrasound leakage signal x(t) is subjected to Gaussian low-pass filtering and denoising to obtain a denoised signal x'(t); the denoised signal x'(t) is decomposed using the EEMD algorithm, and its effective IMF components are selected from the decomposed IMF components based on the correlation coefficient; the energy entropy of the effective IMF components is used as an element to form the characteristic vector T of the original infrasound leakage signal x(t).

[0011] Furthermore, the EEMD algorithm is used to decompose the denoised signal x'(t); specifically:

[0012] Step 11, introduce Gaussian white noise into the denoised signal x'(t), and the calculation formula is:

[0013] e i (t)×x'(t)+n i (t)

[0014] Where e i (t) is the signal obtained after adding the ith noise, n i (t) is the i-th white noise signal;

[0015] Step 12, according to the calculation method in EMD, perform EMD decomposition on the signal after adding white noise to obtain a set of IMF components:

[0016]

[0017] Where m ij (t) is the signal e added with noise after EMD decomposition for the i-th time i (t) The jth IMF component obtained, r ij (t) is the jth residual component obtained by decomposing the signal with noise added for the i-th time;

[0018] Step 13, return to step 11 and repeat N times;

[0019] Step 14, calculate the mean of each IMF component to obtain the final decomposition result. The calculation formula is:

[0020]

[0021] Where m j (t) is the jth IMF component obtained by EEMD decomposition, and r(t) is the residual component obtained by EEMD decomposition.

[0022] Furthermore, the effective IMF components are selected from the decomposed IMF components based on the correlation coefficient; specifically,

[0023] The correlation coefficient formula is used to calculate the correlation coefficient between each IMF component and the denoised signal x'(t);

[0024]

[0025] Where: cov(x,y) is the covariance of x and y; D(x) and D(y) are the variances of x and y respectively; ρ xy The value range is [0,1], x, y represent the IMF component to be calculated and the denoised signal x'(t) respectively;

[0026] The IMF components are sorted from large to small according to the correlation coefficient, and the first J IMF components are taken as valid IMF components, where the sum of the energy values ​​of the first J-1 IMF components is less than a preset percentage of the total energy value of the denoised signal x'(t), and the energy values ​​of the first J IMF components are not less than a preset percentage of the total energy value of the denoised signal x'(t).

[0027] Furthermore, the calculation formula of energy entropy is:

[0028]

[0029] Where m j (t) is the jth effective IMF component, E j is the energy entropy of the jth effective IMF component.

[0030] Furthermore, it also includes normalizing the feature vector T; wherein,

[0031] During the training process, each feature vector T in the training set is normalized and then input into the HMM model training;

[0032] During the recognition process, the feature vector T of the leakage signal x(t) of the original infrasound wave to be identified is normalized and then input into each trained APSO-HMM model;

[0033] T=[E1,E2,...,E J ];

[0034]

[0035] T'=[E1 / E,E2 / E,...,E J / E]

[0036] Where E1, E2, ..., E J represents the energy entropy of the 1st to Jth effective IMF components, and T' represents the normalized vector.

[0037] Furthermore, the APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; specifically,

[0038] Step 21, first define the number of groups as M, each particle searches for the best position in the D-dimensional space at a certain speed, and the optimal position of the nth particle is defined as pbest n ; The optimal position of the group is defined as gbest. In each iteration, each particle and group finds its optimal position and updates its speed and position according to the formula, thus completing an iteration process;

[0039] The speed update formula is:

[0040]

[0041] The position update formula is:

[0042]

[0043] Where: is the velocity and position of the nth particle in the tth iteration in d dimension; is the velocity and position of the nth particle in the t-1th iteration in d dimension; is the optimal position of the nth particle in the t-1th iteration; is the optimal position of the group in the t-1th iteration; w is the inertia weight, c1 and c2 are learning factors, r1, r2∈[0,1] are random numbers; define the following variables:

[0044]

[0045] Where: is the current global optimal value, is the global optimal value of the previous game, α t is the population evolution speed: the ratio of the current global optimal value to the global optimal value of the previous game, α tThe smaller the value of , the faster the population evolves, and w should be increased to improve the global search performance of the population, and vice versa; is the average fitness of the current population, β t is the population aggregation speed factor at the tth iteration, representing the population aggregation degree: the ratio of the global optimum to the average fitness of the current population. The larger its value, the easier it is for the population to fall into the local optimum. w should be increased to make the particles diverge, so that the population can easily jump out of the local area, and vice versa. t is the population premature convergence factor at the tth iteration, representing γ t The smaller the value of , the earlier the population tends to converge prematurely;

[0046] Step 22, when w shows a linear or nonlinear decrease, the algorithm is more efficient, and the mapping relationship between w and γ is designed as

[0047]

[0048] Where: w t is the inertia weight in the tth iteration;

[0049] When the number of iterations is less than the preset value T rec , that is, t≤T rec When , w always maintains a larger value to allow particles to search and expand in a larger range. The expansion method is:

[0050]

[0051] When the preset iteration value is reached, the inertia weight factor should be adaptively reduced, that is, t>T rec hour,

[0052]

[0053] Where: w min and w max are the upper and lower bounds of the inertia weight factor, q1 and q2 are empirical constants, and T max is the maximum number of iterations, τ is the threshold;

[0054] Step 23: To facilitate analysis and expression, simplify the problem space into one dimension and define pbest i =p i ,gbest=g i , then the speed and position update formula is

[0055]

[0056] Through the iterative relationship, we get

[0057]

[0058] Where: is the velocity and position of the ith particle in the tth iteration in d dimensions; is the velocity and position of the ith particle in the d-dimensional space during the t+1th iteration; is the velocity and position of the ith particle in the d-dimensional space during the t+2th iteration;

[0059] Combining the two, we can get the second-order difference equation for particle velocity:

[0060]

[0061] By analyzing the stability of the equation, we can obtain the following adaptive factors:

[0062]

[0063] In the formula, is the value of learning factors c1 and c2 in the tth iteration;

[0064] Step 24, after improving the weight factor in step 22 and the learning factor in step 23, the optimal solution B is obtained.

[0065] Furthermore, the steps for the model to obtain the log-likelihood probability are:

[0066] Step 41, local state initialization:

[0067] δ t (i) = π i b i (O1),i=1,2,3,...N

[0068] In the formula, δ t (i) is the maximum probability of all transition paths of hidden state i at time t.

[0069]

[0070] In the formula, The hidden state of the previous node with the maximum probability among all transfer paths of hidden state i at time t;

[0071] Step 42, perform dynamic programming to recursively calculate the local state at time t (t=1, 2, 3, ..., T);

[0072]

[0073] Step 43, solve δ at time T T (i) Maximum probability of occurrence, and solve for the maximum δ at time T T (i)

[0074] P * =max 1≤j≤N δ T (i)

[0075]

[0076] Step 44, pass Start backtracking, that is

[0077]

[0078] Get the hidden state sequence with the highest probability

[0079]

[0080] The recursive expression of δ is obtained as

[0081]

[0082] estimate The recursive expression of

[0083]

[0084] The maximum probability value P* calculated at the last moment T is the log-likelihood probability.

[0085] After adopting the above technical scheme, the present invention searches for the optimal solution of parameter B through an adaptive particle swarm optimization algorithm, which can well solve the problem of rapid convergence during search, obtain the optimal solution, avoid the problem of premature convergence, and thus improve the recognition effect of pipeline leakage aperture. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 It is the experimental site map of the present invention;

[0087] Figure 2 This is a schematic diagram of a pipeline leakage simulation test of the present invention;

[0088] Figure 3 It is a flow chart of the pipeline leakage aperture identification method based on APSO-HMM of the present invention;

[0089] Figure 4 It is a flow chart of extracting characteristic vector T from infrasound original leakage signal x(t) of the present invention;

[0090] Figure 5 is a time domain waveform diagram of an original infrasound leakage signal x(t);

[0091] Figure 6 for Figure 4The time domain waveform of the original infrasound leakage signal x(t) after denoising;

[0092] Figure 7 The correlation coefficient diagram of each IMF component obtained by EEMD analysis of the denoised signal;

[0093] Figure 8 Adaptive particle swarm optimization (APSO) algorithm is used to optimize the initial values ​​of observations under each state, i.e. the initial value graph of B;

[0094] Fig. 9 Comparison of the model fitness convergence value and convergence speed between the APSO-HMM model and the conventional HMM model with an aperture of 1 mm;

[0095] Fig.10 The difference between the APSO-HMM model and the conventional HMM model in terms of model fitness convergence value and convergence speed under 2mm aperture;

[0096] Fig.11 Comparison of the differences in model fitness convergence value and convergence speed between the APSO-HMM model and the conventional HMM model with an aperture of 3 mm. DETAILED DESCRIPTION

[0097] In order to make the contents of the present invention more clearly understood, the present invention is further described in detail below based on specific embodiments in conjunction with the accompanying drawings.

[0098] like Figure 3 As shown, a pipeline leakage aperture identification method based on APSO-HMM includes:

[0099] The first step is to construct multiple training sets, each training set corresponding to a pipeline leakage aperture; wherein the construction process of each training set includes: obtaining multiple infrasound original leakage signals x(t) corresponding to the pipeline leakage aperture; extracting the feature vector T of each infrasound original leakage signal x(t), and the set formed by the feature vectors T of the multiple infrasound original leakage signals x(t) is used as the training set of the HMM model for the pipeline leakage aperture;

[0100] In the second step, each training set is input into an HMM model for training iteration. During the training process, the APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; the optimal initial value obtained by the APSO algorithm is brought into the HMM model, and the HMM model is iteratively calculated using the Baum-Welch algorithm to obtain a trained APSO-HMM model;

[0101] The third step is to extract the feature vector T of the original infrasound signal leakage x(t) to be identified, input the feature vector T into each trained APSO-HMM model, and output the pipeline leakage aperture corresponding to the APSO-HMM model with the largest log-likelihood probability, which is the identification result;

[0102] Among them, Figure 4 As shown in the figure, the steps of extracting the characteristic vector T of the original infrasound leakage signal x(t) are as follows: the original infrasound leakage signal x(t) is subjected to Gaussian low-pass filtering and denoising to obtain a denoised signal x'(t); since the EEMD algorithm can completely decompose the non-stationary signal into independent IMF signals of different scales, the decomposed signal still has good signal integrity, the EEMD algorithm is used to decompose the denoised signal x'(t); the effective IMF components are selected from the decomposed IMF components based on the correlation coefficient; the energy entropy of the effective IMF components is used as an element to form the characteristic vector T of the original infrasound leakage signal x(t).

[0103] It should be noted that the leakage aperture is continuous, and the leakage aperture interval can be set to several point values ​​as needed. In the process of constructing the training set, the leakage aperture is considered to be the point value to which the actual value of the leakage aperture is closer. After the training is completed, in the recognition process, the point value of the pipeline leakage aperture corresponding to the APSO-HMM model with the largest log-likelihood probability is closer to the actual value of the leakage aperture and can be used as the identified leakage aperture.

[0104] Specifically, the Baum-welch algorithm has a strong dependence on the initial value λ=(π, A, B). So far, the user gives the empirical estimated value of each parameter at the initial moment of the hidden Markov model. Through continuous iteration, each parameter gradually tends to a more reasonable and optimal value. The inventor of this patent found that the parameters π and A have little effect on the selection of the initial value. They can be randomly selected or uniformly selected as long as the random constraint conditions are met, but the selection of B will directly affect the accuracy of the model training. In this embodiment, in order to avoid improper selection of the initial value of the model parameter, the parameter B is optimized by an adaptive particle swarm optimization algorithm. By improving the existing pipeline leakage identification method, the traditional hidden Markov search is improved to quickly converge to obtain the optimal solution, and the premature convergence problem leads to the disadvantage of reduced accuracy of the leakage signal, thereby improving the accuracy of urban pipeline leakage identification, preventing micro-inciple and maintaining urban public safety. Improve the current situation of inaccurate urban pipeline leakage identification, reduce potential energy waste and harm, reduce urban environmental damage and pollution, save energy, and improve economic benefits. Identify pipeline leaks in a timely manner and take measures to control them, thereby eliminating potential accidents, reducing or controlling the occurrence of pipeline accidents, protecting urban public safety, and reducing the loss of life and property of urban residents.

[0105] In one embodiment, the original infrasound leakage signal x(t) is subjected to Gaussian low-pass filtering and denoising, wherein the transfer function of the Gaussian low-pass filter is expressed as follows:

[0106]

[0107] Where G(x,y) is the transfer function of the Gaussian low-pass filter; δ is the standard deviation of the normal distribution. The size of δ will affect the final processing effect. The larger δ is, the more obvious the processing effect on noise will be, but the image will also be more blurred. This value is best selected by comparison in practical applications; (x,y) is the coordinate of the pixel point, which is the point coordinate and can be considered as an integer in image processing.

[0108] In one embodiment, EEMD is an improvement based on the EMD method, which is a method to assist data decomposition by adding noise signals. This method introduces Gaussian white noise into the original signal. Since its mean is zero, it will be eliminated after multiple average calculations, which realizes the automatic distribution of the signal on a suitable time scale, and the final integrated mean is the calculation result, which effectively solves the problem of modal aliasing. The EEMD method shows better feature extraction capabilities in processing non-stationary and nonlinear complex time series data. The superposition effect of the decomposed data is the superposition of the original data and Gaussian white noise. The ensemble average calculation of the obtained IMF components can make the added Gaussian white noise cancel each other. When the number of iterations is large enough, it can be regarded as the decomposition result of the original signal. The EEMD algorithm is used to decompose the denoised signal x'(t); specifically:

[0109] Step 11, introduce Gaussian white noise into the denoised signal x'(t), and the calculation formula is:

[0110] e i (t) = x'(t) + n i (t)(2)

[0111] Where e i (t) is the signal obtained after adding the i-th noise, n i (t) is the i-th white noise signal;

[0112] Step 12, according to the calculation method in EMD, perform EMD decomposition on the signal after adding white noise to obtain a set of IMF components:

[0113]

[0114] Where m ij (t) is the signal e added with noise after EMD decomposition for the i-th time i (t) The jth IMF component obtained, r ij(t) is the jth residual component obtained by decomposing the signal with noise added for the i-th time;

[0115] Step 13, return to step 11 and repeat N times;

[0116] Step 14, in order to offset the influence of Gaussian white noise added in each calculation, the mean of each IMF component is calculated to obtain the final decomposition result. The calculation formula is:

[0117]

[0118] Where m j (t) is the jth IMF component obtained by EEMD decomposition, and r(t) is the residual component obtained by EEMD decomposition.

[0119] After the infrasound signal is decomposed by the EEMD algorithm, several independent IMFs of different scales are obtained. Each IMF has different physical meanings at its own scale. In order to avoid the redundant signals in the signal from misleading the recognition results and reduce the redundancy of the signal, the decomposed IMF signals need to be reasonably screened. In one embodiment, the effective IMF components are selected from the decomposed IMF components based on the correlation coefficient, and the redundant components are eliminated; specifically,

[0120] The correlation coefficient formula is used to calculate the correlation coefficient between each IMF component and the denoised signal x'(t);

[0121]

[0122] Where: cov(x,y) is the covariance of x and y; D(x) and D(y) are the variances of x and y respectively; ρ xy The value range is [0,1]. The correlation coefficient is close to 0, indicating low correlation; the correlation coefficient is close to 1, indicating high correlation; x, y represent the IMF component to be calculated and the denoised signal x'(t) respectively;

[0123] The IMF components are sorted from large to small according to the correlation coefficient, and the first J IMF components are taken as valid IMF components, wherein the sum of the energy values ​​of the first J-1 IMF components is less than a preset percentage of the total energy value of the denoised signal x'(t), and the energy values ​​of the first J IMF components are not less than a preset percentage of the total energy value of the denoised signal x'(t), and the preset percentage may be but is not limited to 90%.

[0124] The calculation formula of energy entropy is:

[0125]

[0126] Where m j (t) is the jth effective IMF component, Ej is the energy entropy of the jth effective IMF component.

[0127] In one embodiment, considering that the energy is sometimes relatively large, in order to facilitate subsequent analysis and processing, T is adjusted by normalizing the features. Therefore, the method further includes normalizing the feature vector T; wherein,

[0128] During the training process, each feature vector T in the training set is normalized and then input into the HMM model training;

[0129] During the recognition process, the feature vector T of the leakage signal x(t) of the original infrasound wave to be identified is normalized and then input into each trained APSO-HMM model;

[0130] T=[E1,E2,...,E J ];

[0131]

[0132] T'=[E1 / E,E2 / E,...,E J / E] (8)

[0133] Where E1, E2, ..., E J represents the energy entropy of the 1st to Jth effective IMF components, and T' represents the normalized vector.

[0134] In one embodiment, the APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; specifically,

[0135] Step 21, first define the number of groups as M, each particle searches for the best position in the D-dimensional space at a certain speed, and the optimal position of the nth particle is defined as pbest n ; The optimal position of the group is defined as gbest. In each iteration, each particle and group finds its optimal position and updates its speed and position according to the formula, thus completing an iteration process;

[0136] The speed update formula is:

[0137]

[0138] The position update formula is:

[0139]

[0140] Where: is the velocity and position of the nth particle in the tth iteration in d dimension; is the velocity and position of the nth particle in the t-1th iteration in d dimension; is the optimal position of the nth particle in the t-1th iteration; is the optimal position of the group in the t-1th iteration; w is the inertia weight, c1 and c2 are learning factors, and r1,r2∈[0,1] are random numbers.

[0141] In practical applications, premature convergence and inappropriate selection of fixed parameters w and c1, c2 will cause the optimization process to fall into the local optimum, affecting the performance of the PSO algorithm. Therefore, an adaptive parameter mechanism will be used to overcome this problem. Define the following variables:

[0142]

[0143] Where: is the current global optimal value, is the global optimal value of the previous game, α t is the population evolution speed: the ratio of the current global optimal value to the global optimal value of the previous game, α t The smaller the value of , the faster the population evolves, and w should be increased to improve the global search performance of the population, and vice versa; is the average fitness of the current population, β t is the population aggregation speed factor at the tth iteration, representing the population aggregation degree: the ratio of the global optimum to the average fitness of the current population. The larger its value, the easier it is for the population to fall into the local optimum. w should be increased to make the particles diverge, so that the population can easily jump out of the local area, and vice versa. t is the population premature convergence factor at the tth iteration, representing γ t The smaller the value of , the earlier the population tends to converge prematurely;

[0144] Step 22, when w shows a linear or nonlinear decrease, the algorithm is more efficient, and the mapping relationship between w and γ is designed as

[0145]

[0146] Where: w t is the inertia weight in the tth iteration;

[0147] When the number of iterations is less than the preset value T rec , that is, t≤T rec When , w always maintains a larger value to allow particles to search and expand in a larger range. The expansion method is:

[0148]

[0149] When the preset iteration value is reached, the inertia weight factor should be adaptively reduced, that is, t>Trec hour,

[0150]

[0151] Where: w min and w max are the upper and lower bounds of the inertia weight factor, q1 and q2 are empirical constants, and T max is the maximum number of iterations, τ is the threshold, τ = 0.002;

[0152] Step 23, in the PSO algorithm, the changes of self and social learning factors are not synchronized. At the beginning of the iteration, c1 increases and c2 decreases. The population has a strong "self-cognition" ability. In the final stage of the iteration, c1 decreases and c2 increases. This embodiment will adopt an adaptive learning factor improvement method. In order to facilitate analysis and expression, the problem space is simplified to one dimension and pbest is defined. i =p i ,gbest=g i , then the speed and position update formula is

[0153]

[0154] Through the iterative relationship, combined with formula (15), we can get

[0155]

[0156] Where: is the velocity and position of the ith particle in the tth iteration in d dimensions; is the velocity and position of the ith particle in the d-dimensional space during the t+1th iteration; is the velocity and position of the ith particle in the d-dimensional space during the t+2th iteration;

[0157] Combining formula (15) and formula (16), we can get the second-order difference equation for particle velocity:

[0158]

[0159] By analyzing the stability of the equation, we can obtain the following adaptive factors:

[0160]

[0161] In the formula, is the value of learning factors c1 and c2 in the tth iteration;

[0162] Step 24, after improving the weight factor in step 22 and the learning factor in step 23, the optimal solution B is obtained.

[0163] In one embodiment, the steps of obtaining the log-likelihood probability of the model are:

[0164] Step 41, local state initialization:

[0165] δ t (i) = π i b i (O1),i=1,2,3,...N (20)

[0166] In the formula, δ t (i) is the maximum probability of all transition paths of hidden state i at time t.

[0167]

[0168] In the formula, The hidden state of the previous node with the maximum probability among all transfer paths of hidden state i at time t;

[0169] Step 42, perform dynamic programming to recursively calculate the local state at time t (t=1, 2, 3, ..., T);

[0170]

[0171] Step 43, solve δ at time T T (i) Maximum probability of occurrence, and solve for the maximum δ at time T T (i)

[0172] P * =max 1≤j≤N δ T (i)

[0173]

[0174] Step 44, pass Start backtracking, that is

[0175]

[0176] Get the hidden state sequence with the highest probability

[0177]

[0178] The recursive expression of δ is obtained as

[0179]

[0180] estimate The recursive expression of

[0181]

[0182] The maximum probability value P* calculated at the last moment T is the log-likelihood probability.

[0183] The pipeline leakage aperture identification method based on APSO-HMM involved in the above embodiment is described in detail below in conjunction with the accompanying drawings of laboratory simulation.

[0184] Figure 1 The experimental site map is shown below. The experimental parameters are as follows: The experimental pipeline system is a U-shaped PE pipe, the pipe is Φ110×10mm, and the total length is 3599cm. The air pressure is 0.1MPa and the flow rate is 16m / s. In the straight section on the left side of the pipeline, an inlet ball valve, an upstream infrasonic sensor, a pressure sensor, and a turbine flowmeter are arranged in sequence from the inlet end of the pipeline. In the straight section on the right side of the pipeline, an outlet ball valve, a downstream infrasonic sensor, a pressure sensor, a turbine flowmeter, and a leakage port B are arranged in sequence from the outlet end of the pipeline. In the curved section (middle horizontal section) of the pipeline, a pressure sensor and a leakage port A are arranged in sequence along the flow direction. Adjustable leakage ports A and B with leakage apertures of 1mm, 2mm, and 3mm are used to simulate small leakages in urban pipelines, and an air compressor is used to simulate the operation process of urban gas pipelines. The upstream and downstream infrasonic sensors are installed at both ends of the pipeline, and connected to the digital network transmission instrument through cables, and then transmitted to the analysis software on the PC through the network cable. See Figure 2 .

[0185] A pipeline leakage aperture identification method based on APSO-HMM includes the following steps:

[0186] S1: When the pressure condition is 0.3MPa, the leaking pipeline is tested. When the leaking pipeline is tested, the leaking aperture size is changed, and the pipeline infrasound signal data under different leaking apertures at 0.3MPa pressure are collected. The leaking aperture sizes are 1mm, 2mm, and 3mm respectively. The waveform of the original infrasound signal data collected under the leaking condition is as follows: Figure 5 shown.

[0187] S2: performing Gaussian low-pass filtering and denoising processing on the original signal of the pipeline leakage infrasound wave to obtain a denoised signal; Figure 5 The waveform in the denoised image is as follows Figure 6 shown.

[0188] S3: The EEMD algorithm is used to completely decompose the denoised signal in step S2 into independent IMF signals of different scales. The size of the energy entropy can reflect the change of the signal energy with the frequency distribution. Therefore, the IMF energy entropy of the infrasound signal is used as the characteristic value to reflect the change of infrasound caused by pipeline leakage.

[0189] In adopting Figure 1In the process of collecting leakage signals with the experimental device, the nine IMF components IMF1 to IMF9 and the residual term IMF10 obtained by the EEMD analysis method are used for the leakage signals under different apertures. The correlation coefficients of the IMF components are calculated according to formula (6). The results are as follows: Figure 7 As shown, the data are shown in Table 1.

[0190] Table 1 shows the data for calculating the correlation coefficients of the IMF components.

[0191]

[0192] According to the calculation results of the correlation coefficient between each component of IMF and the original signal, the conclusions are as follows:

[0193] (1) The correlation coefficients between IMF1~IMF10 and the original signal show a decreasing trend;

[0194] (2) Statistics show that the energy values ​​of the components IMF1 to IMF6 account for more than 90% of the total energy value of the original signal.

[0195] Based on the above two points, refer to formula (7) and take the vector of the root mean square values ​​of IMF1 to IMF6 [E1, E2, E3, E4, E5, E6] as the input key feature value for APSO-HMM pipeline leakage model training and model judgment.

[0196] S4: Use the adaptive particle swarm optimization (APSO) algorithm to optimize the initial value of the observation value under each state, that is, the initial value of B. The comparison before and after optimization is shown in the figure Figure 8 shown.

[0197] S5: The fitness value refers to the value of the objective function. The following is a comparison of the APSO-HMM pipeline leakage model and the conventional HMM pipeline leakage model in terms of model fitness convergence value and convergence speed. Figures 9 to 11 From the results, it can be seen that the three different aperture leakage states can significantly improve the convergence speed of HMM by adding parameter optimization of adaptive particle swarm global optimization and avoid falling into the situation of local minimization.

[0198] Table 2 shows the output log-likelihood probability and recognition results for apertures of 1 mm, 2 mm, and 3 mm

[0199]

[0200] In Table 2, taking operation number 1 (the actual leakage aperture is 1 mm) as an example, its feature vector is input into each trained APSO-HMM model. The likelihood probability is the largest in the model with a leakage aperture of 1 mm, that is, -19.1923>-31.0048>-65.2358. The recognition effect is more obvious, which proves the feasibility of applying APSO-HMM to pipeline leakage identification.

[0201] In order to compare with the conventional HMM pipeline leakage model, the same training samples as the APSO-HMM model were selected for leakage identification using HMM. The labels of the leakage apertures of 1, 2, and 3 mm under the pressure of 0.1 MPa were set to

[123] , where "1" represents the leakage aperture of 1 mm, "2" represents the leakage aperture of 2 mm, and "3" represents the leakage aperture of 3 mm. Accordingly, the data of the pressure of 0.1 MPa and the leakage apertures of 1, 2, and 3 mm were tested in turn, and the result label was

[112] , that is, only the 1 mm aperture leakage was identified, and the accuracy was only 33.33%, which was lower than the recognition result of APSO-HMM.

[0202] Through multiple sample tests, the recognition results of APSO-HMM model and HMM are compared in Table 3. As can be seen from the table, the APSO-HMM model has a good recognition effect.

[0203] Table 3 Comparison of recognition results between APSO-HMM model and HMM model

[0204]

[0205] Based on the above ideal embodiments of the present invention, the relevant staff can make various changes and modifications without departing from the technical concept of the present invention through the above description. The technical scope of the present invention is not limited to the contents of the specification, and its technical scope must be determined according to the scope of the claims.

Claims

1. A pipeline leakage aperture identification method based on APSO-HMM, characterized in that: include: Construct multiple training sets, each training set corresponds to a pipeline leakage aperture; wherein the construction process of each training set includes: obtaining multiple infrasound original leakage signals x(t) corresponding to the pipeline leakage aperture; extracting the feature vector T of each infrasound original leakage signal x(t), and the set formed by the feature vectors T of the multiple infrasound original leakage signals x(t) is used as the training set of the HMM model for the pipeline leakage aperture; Each training set is input into an HMM model for training iteration. During the training process, the APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; the optimal initial value obtained by the APSO algorithm is brought into the HMM model, and the HMM model is iteratively calculated using the Baum-Welch algorithm to obtain a trained APSO-HMM model; Extract the feature vector T of the original infrasound signal leakage x(t) to be identified, input the feature vector T into each trained APSO-HMM model, and output the pipeline leakage aperture corresponding to the APSO-HMM model with the largest log-likelihood probability, which is the identification result; Among them, the steps of extracting the characteristic vector T of the original infrasound leakage signal x(t) are as follows: the original infrasound leakage signal x(t) is subjected to Gaussian low-pass filtering and denoising to obtain a denoised signal x'(t); the denoised signal x'(t) is decomposed using the EEMD algorithm, and its effective IMF components are selected from the decomposed IMF components based on the correlation coefficient; the energy entropy of the effective IMF components is used as an element to form the characteristic vector T of the original infrasound leakage signal x(t).

2. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: The EEMD algorithm is used to decompose the denoised signal x'(t); specifically: Step 11, introduce Gaussian white noise into the denoised signal x'(t), and the calculation formula is: e i (t)=x'(t)+n i (t) Where e i (t) is the signal obtained after adding the ith noise, n i (t) is the i-th white noise signal; Step 12, according to the calculation method in EMD, perform EMD decomposition on the signal after adding white noise to obtain a set of IMF components: Where m ij (t) is the signal e added with noise after EMD decomposition for the ith time i (t) The jth IMF component obtained, r ij (t) is the jth residual component obtained by decomposing the signal with noise added for the i-th time; Step 13, return to step 11 and repeat N times; Step 14, calculate the mean of each IMF component to obtain the final decomposition result. The calculation formula is: Where m j (t) is the jth IMF component obtained by EEMD decomposition, and r(t) is the residual component obtained by EEMD decomposition.

3. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: Based on the correlation coefficient, the effective IMF components are selected from the decomposed IMF components; specifically, The correlation coefficient formula is used to calculate the correlation coefficient between each IMF component and the denoised signal x'(t); Where: cov(x,y) is the covariance of x and y; D(x) and D(y) are the variances of x and y respectively; ρ xy The value range is [0,1], x, y represent the IMF component to be calculated and the denoised signal x'(t) respectively; The IMF components are sorted from large to small according to the correlation coefficient, and the first J IMF components are taken as valid IMF components, where the sum of the energy values ​​of the first J-1 IMF components is less than a preset percentage of the total energy value of the denoised signal x'(t), and the energy values ​​of the first J IMF components are not less than a preset percentage of the total energy value of the denoised signal x'(t).

4. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: The calculation formula of energy entropy is: Where m j (t) is the jth effective IMF component, E j is the energy entropy of the jth effective IMF component.

5. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: It also includes normalizing the feature vector T; wherein, During the training process, each feature vector T in the training set is normalized and then input into the HMM model training; During the recognition process, the feature vector T of the original infrasound signal leakage signal x(t) to be recognized is normalized and then input into each trained APSO-HMM model; T=[E1,E2,...,E J ]; T'=[E1 / E,E2 / E,...,E J / THE] Where E1, E2, ..., E J represents the energy entropy of the 1st to Jth effective IMF components, and T' represents the normalized vector.

6. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: The APSO algorithm is used to optimize the initial value of the observation value under each state, that is, the initial value of B; specifically, Step 21, first define the number of groups as M, each particle searches for the best position in the D-dimensional space at a certain speed, and the optimal position of the nth particle is defined as pbest n The optimal position of the group is defined as gbest. In each iteration, each particle and group finds its optimal position and updates its speed and position according to the formula, thus completing an iteration process. The speed update formula is: The position update formula is: Where: is the velocity and position of the nth particle in the tth iteration in d dimension; is the velocity and position of the nth particle in the t-1th iteration in d dimension; is the optimal position of the nth particle in the t-1th iteration; is the optimal position of the group in the t-1th iteration; w is the inertia weight, c1 and c2 are learning factors, r1, r2∈[0,1] are random numbers; define the following variables: Where: is the current global optimal value, is the global optimal value of the previous game, α t is the population evolution speed: the ratio of the current global optimal value to the global optimal value of the previous game, α t The smaller the value of , the faster the population evolves, and w should be increased to improve the global search performance of the population, and vice versa; is the average fitness of the current population, β t is the population aggregation speed factor at the tth iteration, representing the population aggregation degree: the ratio of the global optimum to the average fitness of the current population. The larger its value, the easier it is for the population to fall into the local optimum. w should be increased to make the particles diverge, so that the population can easily jump out of the local area, and vice versa. t is the population premature convergence factor at the tth iteration, representing γ t The smaller the value of , the earlier the population tends to converge prematurely; Step 22, when w shows a linear or nonlinear decrease, the algorithm is more efficient, and the mapping relationship between w and γ is designed as Where: w t is the inertia weight in the tth iteration; When the number of iterations is less than the preset value T rec , that is, t≤T rec When , w always maintains a larger value to allow particles to search and expand in a larger range. The expansion method is: When the preset iteration value is reached, the inertia weight factor should be adaptively reduced, that is, t>T rec hour, Where: w min and w max are the upper and lower bounds of the inertia weight factor, q1 and q2 are empirical constants, and T max is the maximum number of iterations, τ is the threshold; Step 23: To facilitate analysis and expression, simplify the problem space into one dimension and define pbest i =p i ,gbest=g i , then the speed and position update formula is Through the iterative relationship, we get Where: is the velocity and position of the ith particle in the tth iteration in d dimensions; is the velocity and position of the ith particle in the d-dimensional space during the t+1th iteration; is the velocity and position of the ith particle in the d-dimensional space during the t+2th iteration; Combining the two, we can get the second-order difference equation for particle velocity: By analyzing the stability of the equation, we can obtain the following adaptive factors: In the formula, is the value of learning factors c1 and c2 in the tth iteration; Step 24, after improving the weight factor in step 22 and the learning factor in step 23, the optimal solution B is obtained.

7. The pipeline leakage aperture identification method based on APSO-HMM according to claim 1 is characterized in that: The steps to obtain the log-likelihood probability of the model are: Step 41, local state initialization: δ t (i)=π i b i (O1),i=1,2,3,...N In the formula, δ t (i) is the maximum probability of all transition paths of hidden state i at time t. In the formula, The hidden state of the previous node with the maximum probability among all transfer paths of hidden state i at time t; Step 42, perform dynamic programming to recursively calculate the local state at time t (t=1, 2, 3, ..., T); Step 43, solve δ at time T T (i) Maximum probability of occurrence, and solve for the maximum δ at time T T (i) P * =max 1≤j≤N δ T (i) Step 44, pass Start backtracking, that is Get the hidden state sequence with the highest probability The recursive expression of δ is obtained as estimate The recursive expression of The maximum probability value P* calculated at the last moment T is the log-likelihood probability.

Images