A method for identifying flow patterns of oil-gas-water three-phase flow by multi-mode ultrasonic

Abstract

The application provides a horizontal pipeline oil-gas-water three-phase flow pattern recognition method based on multi-mode ultrasonic waves, which combines a single-crystal pulse wave ultrasonic sensor with a double-crystal continuous wave ultrasonic sensor to respectively obtain echo intensity data reflecting the distribution characteristics of a phase interface and Doppler frequency shift data reflecting the fluctuation characteristics of a flow velocity; based on the response characteristics of test data under different flow patterns, the data fluctuation characteristics are extracted from multiple angles of time domain, frequency domain and space domain, and the equidistant feature mapping method is used for feature dimension reduction, so that the multi-domain joint representation and quantification of flow characteristics are completed; and a nonlinear support vector machine is used to build a classifier, so that the horizontal pipeline oil-gas three-phase flow pattern is accurately recognized.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of multiphase flow measurement, and relates to a method for realizing non-invasive identification of flow patterns of oil-gas-water three-phase flow in a horizontal pipeline by using an ultrasonic sensor in two different working modes of pulse wave and continuous wave, and through effective feature extraction and dimension reduction. BACKGROUND

[0002] Oil-gas-water three-phase flow is a common form of output and raw material in the oil and gas industry, often appearing in related pipelines in the energy and chemical industries, and is related to production processes such as mining, transportation, storage and refining. The mixing of three-phase media with different densities, viscosities and other dynamic properties results in a complex non-uniform spatiotemporal distribution structure of the phase interface in oil-gas-water three-phase flow, which is called flow pattern. Flow pattern has an important influence on the pressure drop characteristics, mass transfer / heat transfer characteristics and the like of the fluid, and is not only related to the efficiency improvement and safety guarantee in industrial production processes, but also has important guiding significance for the accurate measurement of multiphase flow process parameters and the establishment of related mechanism models. Therefore, accurate identification of the flow pattern of oil-gas-water three-phase flow is an indispensable basic problem in the field of multiphase flow.

[0003] The accurate recognition of flow patterns often depends on high-precision sensing technology, effective feature extraction method and appropriate classification algorithm. Among them, high-precision sensing technology is used to obtain the process fluctuation signals of the fluid, which are usually closely related to the distribution or fluctuation characteristics of the flow pattern, such as flow rate, concentration, pressure, etc. At present, many researchers at home and abroad have developed various fluid measurement sensors based on different sensitive principles according to the specific characteristics of the measured object, including electrical sensors (probe, wire mesh, planar electrode, cross-section distributed electrode, etc.), optical sensors (laser Doppler velocimeter, particle image velocimeter, etc.), radiation sensors (X-ray, gamma ray, etc.), microwave sensors (antenna, resonant cavity, etc.), ultrasonic sensors, etc. Different sensing technologies have their own characteristics and application limitations. For example, electrical probe and wire mesh sensor are invasive measurement, which will interfere with the flow field and bring potential threat of pressure loss and blockage; optical sensor has strict requirements on the light transmittance of fluid and pipeline, and the test system is complex and costly; radiation sensor has radioactivity, and strict radiation isolation protection is required to ensure personnel and environmental safety; microwave sensor is easily affected by fluid mineralization and flow pattern, and the design of hardware electronic circuit and sensor is difficult; ultrasonic sensor is easily affected by temperature fluctuation. In actual application, the measurement sensor should be reasonably selected according to the specific use scene. Effective feature extraction method is to extract objective characteristic parameters from test signals to realize effective representation of flow state according to the response characteristics of test signals to flow pattern on the basis of sensing technology. At present, time domain statistical characteristics (mean, variance, maximum, minimum, etc.), frequency domain characteristics (maximum frequency, average frequency, etc.), time-frequency domain characteristics (band correlation coefficient, energy proportion, etc.), nonlinear characteristics (entropy, fractal dimension, etc.) have been used to represent flow characteristics. Due to the complexity of the flow process, sometimes multiple features need to be combined to achieve a comprehensive description of the flow characteristics, but the redundancy between features needs to be considered. Finally, the extracted flow characteristics are sent to the classifier to realize the recognition of flow pattern. Compared with traditional classification methods, machine learning-based methods such as support vector machine, decision tree, random forest, artificial neural network, etc. have obvious advantages in solving the nonlinear and random problem of multiphase flow pattern recognition, and have obtained many successful applications.

[0004] In the past few decades, many domestic and foreign researchers have proposed various flow pattern recognition methods, but the objects are mostly concentrated on the relatively simple gas-water or oil-water two-phase flow. Compared with two-phase flow, oil-gas-water three-phase flow has more complex phase distribution structure and phase interaction, and the accurate recognition of its flow pattern needs to rely on more abundant test information, more effective feature expression and more superior classification algorithm. At present, the multi-sensor fusion strategy has become an effective solution to improve the accuracy of flow pattern recognition. However, in actual engineering application, the increase of sensing modalities can greatly improve the complexity, installation difficulty and application cost of the test system, and the most ideal flow pattern recognition solution should be accurate and simple. Among various sensing technologies, ultrasonic method can form different measurement modes by using different acoustic phenomena such as transmission, reflection, scattering, diffraction and Doppler emitted by sound waves in multiphase flow, which can be used to obtain various flow parameter information reflecting holdup, interface position and flow velocity, and has the advantages of rich information acquisition, simple sensor structure, no interference with flow state, easy installation and maintenance, safety and no radiation, low cost and the ability to measure non-transparent fluids in non-transparent pipelines. Therefore, the present application proposes a three-phase flow pattern recognition method based on multi-mode ultrasonic, which involves the combination of continuous wave and pulse wave acoustic testing means, multi-dimensional feature extraction and dimension reduction method and support vector machine classification algorithm. SUMMARY

[0005] The purpose of the present application is to use pulse wave and continuous wave ultrasonic sensors with two different working modes to provide an accurate, non-invasive and low-cost horizontal pipeline oil-gas-water three-phase flow pattern recognition method. The technical solution is as follows:

[0006] A horizontal pipeline oil-gas-water three-phase flow pattern recognition method based on multi-mode ultrasonic combines single-crystal pulse wave ultrasonic sensors and double-crystal continuous wave ultrasonic sensors to respectively obtain echo intensity data reflecting the distribution characteristics of the phase interface and Doppler frequency shift data reflecting the flow velocity fluctuation characteristics; based on the response characteristics of the test data under different flow patterns, the data fluctuation characteristics are extracted from multiple angles of time domain, frequency domain and spatial domain, and the equidistant feature mapping method is used for feature dimension reduction to complete the multi-domain joint representation and quantification of flow characteristics; a nonlinear support vector machine is used to build a classifier to accurately identify the horizontal pipeline oil-gas three-phase flow; including the following steps:

[0007] (1) synchronously collecting horizontal pipeline oil-gas three-phase flow information under different working conditions by using pulse wave ultrasonic sensors and continuous wave ultrasonic sensors: the pulse wave ultrasonic sensor comprises a piezoelectric ceramic wafer installed vertically to the pipe wall, which is used to obtain echo intensity profile data reflecting the distribution characteristics of the fluid phase interface; the continuous wave ultrasonic sensor comprises two piezoelectric ceramic wafers installed obliquely, which is used to obtain Doppler frequency shift data reflecting the flow velocity fluctuation characteristics of the fluid;

[0008] (2) In the echo intensity profile data, the radial position of the maximum echo amplitude in the profile obtained in each pulse repetition period is located, and the probability distribution thereof is counted; based on the statistical result, the kurtosis, skewness, weighted average value and position corresponding to the maximum probability are extracted in the probability distribution graph, which are respectively used to quantify the flatness, symmetry, fluctuation distribution and extreme value of the probability distribution, as the spatial domain features of the flow pattern representation;

[0009] (3) The time-frequency spectrum is obtained by performing short-time Fourier transform on the Doppler shift data, and the average frequency f(τ) of the spectrum at τ time and the corresponding fluid average flow velocity u(τ) are calculated; based on the calculation result, the mean, variance and average zero-crossing rate in the average flow velocity fluctuation time sequence are extracted, which are respectively used to quantify the concentration tendency, fluctuation degree and morphological change degree of the data, as the time domain features of the flow pattern representation;

[0010] (4) The ensemble empirical mode decomposition method is used to decompose the Doppler shift data into a sum of several intrinsic mode functions and residual terms; based on the decomposition result, the correlation coefficients of each level intrinsic mode function and the original Doppler shift data are calculated, and the first several intrinsic mode function components with larger correlation coefficients are selected as the principal components representing the significant features of the Doppler shift data, and the energy proportion corresponding to the principal components is taken as the multi-scale domain feature of the flow pattern representation;

[0011] (5) A multi-dimensional feature vector based on pulsed wave and continuous wave ultrasonic test data is formed;

[0012] (6) The multi-dimensional feature vector is used as an input data set, and the isometric feature mapping method is used to reduce the dimension of the feature vector to obtain a reduced feature vector;

[0013] (7) The oil-gas-water three-phase flow pattern label is added to the reduced feature vector, and the feature vectors under different working conditions are formed into a data set;

[0014] (8) The data set is divided into a training data set and a test data set according to a certain proportion, the training data set is used for training the flow pattern classifier model, and the test data set is used to evaluate the classification error of the trained classifier model;

[0015] (9) A nonlinear support vector machine flow pattern classifier for outputting the flow pattern of an unknown sample is obtained.

[0016] Further, in step (1), the pulsed wave ultrasonic sensor intermittently emits ultrasonic pulses into the fluid at a certain repetition frequency when working, and receives echoes during the pulse interval. The echo intensity profile data A(y, t PRF ) reflecting the distribution characteristics of the fluid phase interface is obtained by amplitude demodulation of the received echoes, where t PRF ∈(0,nT PRF) is the measurement time, T PRF is the pulse repetition period, n is the number of pulse excitations within the sampling time, y ∈ (0, D) is the distance from the pipe bottom, and D is the pipe diameter.

[0017] Further, in step (1), the continuous wave ultrasonic sensor comprises two piezoelectric ceramic wafers installed at an angle θ (45° ≤ θ ≤ 60°) with the horizontal direction, which are respectively responsible for continuously transmitting and receiving ultrasonic waves in the fluid during operation; the Doppler shift data f d (t), where t ∈ (0, T) is the time and T is the sampling time.

[0018] Further, the method of step (2) is:

[0019] In the echo intensity profile data A(y, t PRF ), t PRF ∈ (0, nT PRF ), the position where the maximum echo amplitude value occurs is located in each profile obtained within each pulse repetition period, and the probability distribution thereof is counted:

[0020]

[0021] where P(y) represents the probability of the maximum echo amplitude value A max occurring in the position interval [y, y+Δy] within the sampling time, and Num(y≤A max ≤y+Δy) represents the number of times the maximum echo amplitude value A max occurs in the position interval [y, y+Δy] within the sampling time.

[0022] Based on the statistical results, the kurtosis P kur , skewness P ske , weighted average P wei , and position P lo corresponding to the maximum probability are extracted from the probability distribution graph, which are respectively used to quantify the flatness, symmetry, fluctuation distribution, and extreme value of the probability distribution as spatial domain features of the flow pattern representation:

[0023]

[0024] where E(*) represents the mean value operation, P max is the maximum probability value in the probability distribution, P mean and P var are the mean and variance of the probability distribution statistical values, respectively.

[0025] Further, the method of step (3) is:

[0026] Doppler shift data f d (t),t∈(0,T) is subjected to short-time Fourier transform to obtain time-frequency spectrum S d (τ,f), and the average frequency of the spectrum at time τ is calculated by applying the following formula to each column of S d (τ,f) and the corresponding fluid average flow velocity

[0027]

[0028] where c is the sound speed in the acoustic wedge, f0 is the frequency of the emitted sound wave, and θ is the angle between the central axis of the piezoelectric ceramic wafer and the horizontal direction.

[0029] Based on the calculation results, the mean value , variance , and average zero-crossing rate of the average flow velocity fluctuation time series are extracted, which are used to quantify the central tendency, fluctuation degree, and morphological change intensity of the data, respectively, as time-domain features of the flow pattern representation:

[0030]

[0031] where sgn(*) is the sign function, is the data obtained by subtracting the mean value of the average flow velocity fluctuation time series from itself, and M is the total number of data points in the average flow velocity fluctuation time series.

[0032] Further, the method of step (4) is:

[0033] The Doppler shift data f d (t),t∈(0,T) is decomposed into a sum of intrinsic mode functions IMF j (t)(j=1,2,...,N) and a residual term r N (t) using the ensemble empirical mode decomposition method:

[0034]

[0035] where N is the number of intrinsic mode functions obtained by decomposition;

[0036] Based on the decomposition results, the correlation coefficients of each level of intrinsic mode function and the original Doppler shift data are calculated, and the first b intrinsic mode function components IMF i ,i∈{j1,j2,j3,...,j b} with larger correlation coefficients are selected as the principal components representing the significant features of the Doppler shift data, and the energy proportion R i corresponding to the principal components is taken as a multi-scale domain feature of the flow pattern representation:

[0037]

[0038] Further, in step (5), the multi-dimensional feature vector S is:

[0039]

[0040] Further, in step (6), the multi-dimensional feature vector S is taken as the input data set X = {x i | x1, x2, ···, x b+7}, where x i represents each feature parameter in the multi-dimensional feature vector S, and the isometric feature mapping method is used to reduce its dimension to m-dimensional space (m < b + 7) to obtain the reduced feature vector, which is realized by the following sub-steps:

[0041] (6.1) The k-nearest neighbors of each data point x i in the data set are determined using the K-nearest neighbor search algorithm, and its neighborhood space N i = {x i1 , x i2 , ···, x ik} is obtained, i ∈ (1, b + 7);

[0042] (6.2) The neighborhood connection graph G = [V, E] is constructed, where V represents the nodes in the neighborhood connection graph, which are all data points in the data set X, and E represents the edges in the neighborhood connection graph, whose length is measured by the distance d G between each data point in the data set X: when the data x i and x j belong to the same neighborhood space, d G is set as the Euclidean distance between the two points, i.e. d G = d2(x i , x j ) = ||x i - x j ||2, and when the data x i and x j do not belong to the same neighborhood space, d G = ∞, indicating no connection;

[0043] (6.3) The shortest distance between any two points in the neighborhood connection graph G is determined using the shortest path algorithm, and the geodesic matrix D G = {dist(x i , x j )}, x i , x j ∈ X is obtained, where

[0044]

[0045] (6.4) Constructing low-dimensional embedding: Adopting Multidimensional Scaling (MDS) algorithm to construct kernel matrix where S G = [dist(x i ,x j )] 2 , H = I - J / (b+7), I is a unit diagonal matrix, J is a full 1 matrix, and singular value decomposition is performed on τ(D G ), and the first m eigenvalues λ1≤λ2≤...≤λ m and the corresponding eigenvectors v1,v2,...,v m are taken.

[0046]

[0047] where Y is the feature vector after dimensionality reduction, Λ m = diag(λ1,λ2,...,λ m ), V m = [v1,v2,…,v m ]

[0048] Further, the shortest path algorithm is Dijkstra or Floyds algorithm.

[0049] The substantial features of the present application are: combining ultrasonic sensors of two different working modes of pulse wave and continuous wave to realize non-invasive synchronous measurement of oil-gas-water three-phase flow, respectively acquiring echo intensity profile data reflecting the distribution characteristics of fluid phase interface and Doppler frequency shift data reflecting the fluctuation characteristics of fluid velocity; deeply mining ultrasonic test data, extracting feature parameters from multiple angles of spatial domain, time domain and multi-scale domain to realize objective and comprehensive characterization of the flow process; realizing comprehensive dimensionality reduction of nonlinear global feature relationship in high-dimensional feature space through equidistant feature mapping; adopting nonlinear support vector machine method suitable for small sample size to construct flow pattern classifier, and realizing accurate identification of oil-gas-water three-phase flow pattern.

[0050] 1. The method adopts a multi-mode ultrasonic combined test method, and has the advantages of non-invasiveness, simple structure and low cost.

[0051] 2. The method extracts data features from spatial domain, time domain and multi-scale domain, and can realize comprehensive and objective characterization of the flow process.

[0052] 3. The method uses equidistant feature mapping method to realize feature dimension reduction, which helps to reduce the calculation complexity, improve the model generalization, and reduce the overfitting in the flow pattern classifier construction process. BRIEF DESCRIPTION OF DRAWINGS

[0053] The following drawings and tables describe selected embodiments of the present application, which are all exemplary drawings and tables rather than exhaustive or limiting, wherein:

[0054] Figure 1 The schematic diagram of the pulsed wave and continuous wave ultrasonic testing system in the measurement method of the present application;

[0055] Figure 2 The schematic diagram of the pulsed wave and continuous wave ultrasonic sensor structure in the measurement method of the present application;

[0056] Figure 3 The flow chart of the short-time Fourier transform method in the measurement method of the present application;

[0057] Figure 4 The flow chart of the ensemble empirical mode decomposition method in the measurement method of the present application;

[0058] Figure 5 The algorithm flow chart of the nonlinear support vector machine in the measurement method of the present application;

[0059] Figure 6 The flow chart of the oil-gas-water three-phase flow pattern recognition in the measurement method of the present application. DETAILED DESCRIPTION

[0060] The embodiments of the present application will be described in detail below in combination with the drawings and examples of the specification.

[0061] Figure 1 The schematic diagram of the pulsed wave and continuous wave ultrasonic testing system in the measurement method of the present application. The testing system is composed of three modules, namely, an excitation unit I, a sensing unit II, and an acquisition unit III. The sensing unit II includes a single-crystal pulsed wave ultrasonic sensor 0 and a double-crystal continuous wave ultrasonic sensor 1, which are installed at a certain interval at the bottom of the pipeline and work simultaneously, and the installation sequence is not limited. The excitation unit I includes a signal generator 2 and a power amplifier 3, which are used to generate continuous wave and pulsed wave excitation signals required for the operation of the ultrasonic sensor. The acquisition unit III includes a data processing card 4, a data acquisition card 5, and a computer 6, wherein the data processing card 4 is used to process the fluid-modulated acoustic wave signals received by the ultrasonic sensor and output the echo intensity data related to the phase interface distribution characteristics and the Doppler frequency shift data related to the flow velocity fluctuation characteristics, the data acquisition card 5 synchronously acquires and uploads the output data of the data processing card 4 to the computer 6, and completes the real-time storage of the data.

[0062] Figure 2The figure is a schematic diagram of the pulse wave and continuous wave ultrasonic sensor structure in the measurement method of the present application. The pulse wave ultrasonic sensor 0 and the continuous wave ultrasonic sensor 1 used in the present application are installed at a certain interval at the bottom of the pipeline 7, ensuring that the measurement spaces 8 of the two sensors do not interfere with each other. The pulse wave ultrasonic sensor 0 contains a piezoelectric ceramic wafer 901 with a center frequency of 1 MHz attached to an acoustic wedge 101. When in operation, the piezoelectric ceramic wafer intermittently emits ultrasonic pulses into the fluid at a certain repetition frequency and receives echoes during the pulse interval. The amplitude information of the received acoustic waves is related to the distribution characteristics of the phase interface of the measured fluid. The continuous wave ultrasonic sensor 1 contains two piezoelectric ceramic wafers 902 and 903 with a center frequency of 1 MHz attached to acoustic wedges 102 and 103. When in operation, the two piezoelectric ceramic wafers are responsible for continuously emitting and receiving ultrasonic waves into the fluid, respectively. The frequency offset between the emitted and received ultrasonic waves is related to the velocity of the measured fluid. The acoustic wedges 102 and 103 are cut into a fixed shape to ensure that the angle between the normal of the piezoelectric ceramic wafer and the horizontal direction is θ (45° ≤ θ ≤ 60°), and the beam path direction is against the flow direction 12 of the fluid. In addition, the back of the piezoelectric ceramic wafer in the pulse wave ultrasonic sensor 0 and the continuous wave ultrasonic sensor 1 is filled with sound-absorbing material 11 to effectively reduce the reflection and scattering of acoustic waves inside the sensor, improving the signal-to-noise ratio and resolution of the measurement.

[0063] Figure 3 The figure is a flowchart of the short-time Fourier transform method in the measurement method of the present application. Short-time Fourier transform (STFT) is an important time-frequency analysis method. The idea is to use a shorter window function to slide along the time axis and calculate the frequency characteristics of the signal at different times in different windows, which can effectively capture the transient changes and time-varying characteristics of the signal. Mathematically, the STFT of a time-domain signal s(t) is defined as: where h(*) is the window function, which is used to define the time range of the analysis. Commonly used window functions include rectangular window, Hamming window, Hanning window, Blackman window, etc. Implementing STFT on a given time-domain signal s(t) mainly includes the following steps: 1) setting the window function: select the window function and its length according to the application requirements, and set the frame shift, which determines the sliding step of the window function on the time axis. Common frame shifts are 50% or 75% of the window length; 2) framing: use the window function to divide the signal into multiple overlapping short-time frames. Each frame is multiplied by the window function to obtain the windowed signal segment. The length of each signal segment is the length of the window function; 3) perform Fourier transform on each signal segment to obtain its frequency domain representation; 4) arrange the Fourier transform results of all signal segments in time order to form a time-frequency matrix S(τ, f).

[0064] Figure 4 is a flow chart of the ensemble empirical mode decomposition method in the measurement method of the present application. Ensemble Empirical Mode Decomposition (EEMD) is a data adaptive multi-resolution technique, which can decompose the fluctuation components of different frequencies contained in a signal into intrinsic mode functions of different scales. EEMD is proposed on the basis of the traditional Empirical Mode Decomposition (EMD) method. It utilizes the frequency uniform distribution and zero mean characteristics of white noise, and through adding multiple groups of white noise to the signal, it can overcome the mode mixing problem in the EMD method and improve the accuracy of decomposition, so the EEMD method is also called white noise assisted analysis method. For a given signal x(t), EEMD decomposition mainly includes the following steps: 1) add a Gaussian white noise ω i (t) with a mean of 0 and a standard deviation of 0.1-0.4 times the standard deviation of x(t) to x(t) to construct a new time series X i (t) = x(t) + ω i (t), i = 1, 2, …, g; 2) perform EMD decomposition on X i (t) to obtain where IMF ij is the jth intrinsic mode function of the ith time series, and r in is the decomposition residual term of the ith time series; 3) average the decomposition results of g times to remove the influence of white noise, to obtain the intrinsic mode functions IMF j of different frequencies representing the fluctuation characteristics in the signal and the residual term r N representing the trend characteristics in the signal: Finally, x(t) can be expressed as:

[0065] Figure 5 is a flow chart of the algorithm of the nonlinear support vector machine in the measurement method of the present application. Support Vector Machine (SVM) is a semi-supervised machine learning method based on statistical learning theory and structural risk minimization principle, which is widely used to solve small sample, nonlinear and high-dimensional pattern recognition and regression problems. For nonlinearly separable classification problems, SVM is essentially a binary classifier, and its basic idea is to use kernel functions to map data from low-dimensional linearly inseparable feature space to high-dimensional Hilbert space, and find a separation hyperplane in the high-dimensional space that can separate different classes of samples in the data set and has the maximum geometric interval, which can be realized by convex quadratic programming optimization algorithm. The specific steps are as follows: 1) input the training data set X = {(x1, y1), (x2, y2), …, (x n,y n )}, where x i ∈R D , i=1,2,…,n is n samples in the D-dimensional feature space, y i ∈{-1,+1} is the sample category; 2) Select an appropriate kernel function L(*) and penalty factor C, construct and solve the convex quadratic programming problem: Get the optimal solution Select α * A positive component of Calculation parameters Then the separating hyperplane can be expressed as Use the separating hyperplane to construct the classification decision function: Where sgn(*) represents the sign function. For samples of unknown category, the output of the classification decision function can be used to determine its category.

[0066] Figure 6 This is a flow chart for identifying the oil-gas-water three-phase flow pattern in the measurement method of the present invention.

[0067] The oil-gas-water three-phase flow pattern identification method proposed in this patent involves an acoustic testing method combining continuous waves and pulse waves, a multidimensional feature extraction and dimensionality reduction method, and a support vector machine classification algorithm. The specific implementation steps are as follows:

[0068] Step 1: Exploit Figure 1 、 Figure 2 The ultrasonic testing system and sensor shown synchronously obtain pulse wave and continuous wave ultrasonic measurement data under different flow patterns of oil, gas and water three-phase flow in a horizontal pipeline.

[0069] 1) The pulse wave ultrasonic sensor works in the self-transmitting and self-receiving mode, that is, the piezoelectric ceramic chip is measured at a repetition frequency f PRF Ultrasonic pulses are intermittently emitted into the fluid, and the echoes are received by the same piezoelectric ceramic chip during the pulse interval. By demodulating the amplitude of the received echoes, the echo intensity profile data A(y,t) reflecting the distribution characteristics of the fluid interface can be obtained. PRF ), where t PRF ∈(0,nT PRF ) is the measurement time, T PRF =1 / f PRF is the pulse repetition period, n is the number of pulse excitations within the sampling time, y∈(0,D) is the distance from the bottom of the pipe, and D is the pipe diameter.

[0070] 2) The continuous wave ultrasonic sensor operates in a transmit-receive mode, i.e., during measurement, one piezoelectric ceramic chip continuously transmits ultrasonic waves into the fluid, while the other piezoelectric ceramic chip continuously receives the echoes. Based on the Doppler effect, by frequency demodulating the received echoes, the Doppler frequency shift data f that reflects the fluid velocity fluctuation characteristics can be obtained. d (t), where t∈(0,T) is time and T is the sampling time.

[0071] Step 2: Process the echo intensity data and Doppler shift data, and perform feature extraction from multiple angles to achieve flow pattern characterization.

[0072] 1) In the echo intensity profile data A(y,t PRF ),t PRF ∈(0,nT PRF ), locate the position where the maximum echo amplitude occurs within each pulse repetition period, and calculate its probability distribution:

[0073]

[0074] Among them, P(y) represents the probability that the maximum echo amplitude within the sampling time appears in the position interval [y,y+Δy], Num(y≤A max ≤y+Δy) represents the number of times the maximum echo amplitude appears in the position interval [y,y+Δy] within the sampling time.

[0075] The shape of the probability distribution graph is modulated by the spatial distribution characteristics of the fluid phase interface at the radial position of the pipeline. Therefore, the kurtosis P is extracted from the probability distribution graph. kur , skewness P ske , weighted average P wei The position P corresponding to the maximum probability lo , which are used to quantify the flatness, symmetry, fluctuation distribution and extreme value of the probability distribution, as the spatial domain characteristics of the flow pattern:

[0076]

[0077] Among them, E(*) represents the average operation, P mean and P var The mean and variance of the probability distribution statistics respectively.

[0078] 2) Utilize Figure 3 The short-time Fourier transform algorithm process shown in the figure processes the Doppler frequency shift data f d (t), t∈(0,T), obtain the time spectrum S d (τ,f). S d(τ,f) is a two-dimensional matrix, and each column vector represents the spectrum obtained by Fourier transform after windowing the signal at different time τ. d For each column of (τ,f), the average frequency of the spectrum is calculated using the following formula:

[0079]

[0080] Furthermore, based on the Doppler effect, the average flow velocity u of the fluid at time τ is calculated using the following formula:

[0081]

[0082] Where c0 is the speed of sound in the fluid, f0 is the frequency of the emitted sound wave, and θ0 is the angle between the sound wave in the fluid and the horizontal direction.

[0083] Based on the calculation results, extract the average flow velocity fluctuation time series The mean variance Average zero-crossing rate They are used to quantify the central tendency, fluctuation degree and morphological change intensity of the data, respectively, as the time domain characteristics of flow pattern characterization:

[0084]

[0085] Among them, sgn(*) is the sign function, is the data after subtracting its own mean from the average flow velocity fluctuation time series, and M is the total number of data points in the average flow velocity fluctuation time series.

[0086] 3) Utilize Figure 4 The algorithm flow of the ensemble empirical mode decomposition method shown in the figure is to transform the Doppler frequency shift data f d (t),t∈(0,T) is decomposed into several intrinsic mode functions IMF j (t)(j=1,2,...,N) and the residual term r N The sum of (t):

[0087]

[0088] Where N is the number of intrinsic mode functions obtained by decomposition. The eigenmode functions of each level obtained by ensemble empirical mode decomposition reflect the different frequency components in the original signal, while the residual term represents the trend information in the original signal.

[0089] Based on the decomposition results, calculate the energy proportion R of each level of eigenmode function j :

[0090]

[0091] Further, the correlation coefficient CC of each level intrinsic mode function and the original Doppler shift data is calculated by using the following formula IMFj :

[0092]

[0093] wherein E[*] represents the mean value operation. Generally, the larger the correlation coefficient is, the more the corresponding intrinsic mode function can reflect the most significant characteristics of the original signal.

[0094] The intrinsic mode function components are sorted according to the correlation coefficient The first 10 intrinsic mode function components with larger correlation coefficients are selected as the principal components representing the significant characteristics of the Doppler shift data, and the corresponding energy proportion is taken as the multi-scale domain feature of the flow pattern representation:

[0095] R j ,j={j1,j2,j3,...,j b} (14)

[0096] 4) Form a multi-dimensional feature vector S based on the pulsed wave and continuous wave ultrasonic test data:

[0097]

[0098] The dimension of the multi-dimensional feature vector S is 17.

[0099] Step 3: Use Isometric Feature Mapping method to reduce the dimension of the multi-dimensional feature vector S.

[0100] Isometric Feature Mapping (Isomap) is a popular learning algorithm that can approximate the geodesic distance matrix on the manifold by constructing an adjacency graph and calculating the shortest path, thereby mapping high-dimensional data points to low-dimensional space to achieve data dimension reduction, with the advantages of high computational efficiency, global optimization and good gradual convergence effect, and is widely used in comprehensive dimension reduction of high-dimensional space nonlinear global feature relationship. The multi-dimensional feature vector S is taken as the input data set X = {x i |x1,x2,···,x 17}, wherein x i represents each feature parameter in the multi-dimensional feature vector, and the Isomap method is used to reduce it to m-dimensional space (m < 17), mainly including the following sub-steps:

[0101] 1) Use the K-nearest neighbor search algorithm to determine the k-nearest neighbors of each data point x i in the data set, and obtain its neighborhood space N i = {x i1 , x i2 , ···, x ik}, i e (1, 17);

[0102] 2) Construct a neighborhood connection graph G = [V, E], where V represents the nodes in the neighborhood connection graph, V = {x1, x2,..., x17} is all the data points in the data set X, E represents the edges in the neighborhood connection graph, and the length of each edge is the distance d between the data points in the data set X. G Measure: when data x i belongs to the same neighborhood space as x j , d G is set as the Euclidean distance between the two points, i.e., d G = d2(x i , x j ) = ||x i - x j ||2, and when data x i does not belong to the same neighborhood space as x j , d G = ∞, indicating no connection.

[0103] 3) Determine the shortest distance between any two points in the neighborhood connection graph G using the shortest path algorithm (Dijkstra or Floyd's algorithm) to obtain the geodesic matrix: D G = {dist(x i , x j )}, x i , x j ∈ X, where

[0104] 4) Construct a low-dimensional embedding: use the Multidimensional Scaling (MDS) algorithm to construct a kernel matrix where H = I - J / 17, I is a unit diagonal matrix, and J is a full 1 matrix, and perform singular value decomposition on τ(D G ), and take the first m eigenvalues λ1≤ λ2≤...≤ λ m and the corresponding eigenvectors v1, v2,..., v m , then the new coordinates of the original data set X formed by the multidimensional eigenvectors S in the m-dimensional space of low dimension can be represented as:

[0105]

[0106] where Y is the feature vector after dimensionality reduction, Λ m = diag(λ1, λ2,..., λ m ), and V m = [v1, v2,..., v m ].

[0107] Step 4: A flow pattern classifier is constructed by using the reduced feature vectors by a nonlinear support vector machine method to realize oil-gas-water three-phase flow pattern recognition, including the following sub-steps:

[0108] 1) Adding labels to the data set: according to the actual observation, adding flow pattern labels to the reduced feature vector data set to form a data set W. In the embodiment of the present application, a total of 610 working conditions of oil-gas-water three-phase flow ultrasonic measurement data are collected, including 8 flow pattern labels, which are water-based-oil-water separation-layer flow, water-based-oil-water dispersion-layer flow, water-based-oil-water separation-wave flow, water-based-oil-water dispersion-wave flow, water-based-oil-water dispersion-plug flow, water-based-oil-water dispersion-elastic flow, water-based-oil-water dispersion-penetrating elastic flow and water-based-oil-water dispersion-ring flow.

[0109] 2) Using the leave-one-out method to divide the data set: dividing the data set into two subsets Q and Z according to the ratio of 7:3, satisfying Q∪Z=W and Wherein the subset Q will be used as training data set for training the flow pattern classifier model, and the subset Z will be used as test data set for evaluating the classification error of the trained classifier model.

[0110] 3) Constructing a flow pattern classifier based on a nonlinear support vector machine method and testing: the classical SVM algorithm is designed for solving two-classification problems, and in solving multi-classification problems, a "one-to-one" or "one-to-many" multi-classifier can be constructed. In the embodiment of the present application, a "one-to-one" strategy is used to construct a multi-classifier for oil-gas-water three-phase flow pattern recognition. Specifically, the algorithm shown in the formula is used to train a decision function for any two flow pattern sample data sets, and a total of (8×7) / 2=28 decision functions can be obtained under 8 flow patterns, corresponding to 28 SVM two-classifiers, and each SVM two-classifier only classifies between two flow patterns. When the test sample is flow pattern predicted, the prediction results of the 28 SVM two-classifiers are voted and counted, and the result with the highest number of votes is the flow pattern prediction result of the test sample. Figure 5

[0111] Table 1 Oil-gas-water three-phase flow pattern recognition results

[0112]

[0113]

[0114] ​Table 1 is the identification result of the oil-gas-water three-phase flow pattern in the measuring method of the present application. From the table, it can be seen that the identification rates of the eight flow patterns of water-based-oil-water separation-laminar flow, water-based-oil-water dispersion-laminar flow, water-based-oil-water separation-wavy flow, water-based-oil-water dispersion-wavy flow, water-based-oil-water dispersion-plug flow, water-based-oil-water dispersion-elastic flow, water-based-oil-water dispersion-penetrating elastic flow and water-based-oil-water dispersion-ring flow are 93.33%, 100%, 88.89%, 94.11%, 90.00%, 93.75%, 93.33% and 100% respectively, the average identification rate of all flow patterns is 93.44%, indicating that the method is effective.

[0115] Finally, it should be noted that: the above examples are used to illustrate the technical solutions of the present application, but not to limit it. Although the present application has been described in detail with reference to the above examples, those skilled in the art should understand that the specific embodiments of the present application can be modified or replaced by equivalents without departing from the spirit and scope of the present application. Any modification or equivalent replacement without departing from the spirit and scope of the present application should be covered within the protection scope of the claims of the present application.

Claims

1. A multi-mode ultrasound-based method for identifying flow patterns in horizontal pipeline oil-gas-water three-phase flows. This method combines a single-crystal pulsed-wave ultrasonic sensor with a dual-crystal continuous-wave ultrasonic sensor to acquire echo intensity data reflecting the distribution characteristics of the phase interface and Doppler frequency shift data reflecting the flow velocity fluctuation characteristics, respectively. Based on the response characteristics of test data under different flow patterns, data fluctuation characteristics are extracted from multiple perspectives in the time, frequency, and spatial domains. Feature dimensionality reduction is performed using an equidistant feature mapping method to achieve multi-domain joint characterization and quantification of flow characteristics. A classifier is constructed based on a nonlinear support vector machine to accurately identify flow patterns in horizontal pipeline oil-gas three-phase flows. The method comprises the following steps: (1) Using a pulse wave ultrasonic sensor and a continuous wave ultrasonic sensor to synchronously collect flow process information of oil and gas three-phase flow in a horizontal pipeline under different working conditions: the pulse wave ultrasonic sensor includes a piezoelectric ceramic chip installed perpendicular to the pipe wall, which is used to obtain echo intensity profile data reflecting the distribution characteristics of the fluid phase interface; the continuous wave ultrasonic sensor includes two obliquely installed piezoelectric ceramic chips, which are used to obtain Doppler frequency shift data reflecting the fluid velocity fluctuation characteristics; (2) In the echo intensity profile data, the radial position of the pipeline where the maximum echo amplitude occurs is located in the profile obtained within each pulse repetition period, and its probability distribution is statistically analyzed; based on the statistical results, the kurtosis, skewness, weighted average value, and the position corresponding to the maximum probability are extracted from the probability distribution diagram, which are used to quantify the flatness, symmetry, fluctuation distribution, and extreme value of the probability distribution, respectively, as spatial domain features for flow pattern characterization; (3) Perform short-time Fourier transform on the Doppler frequency shift data to obtain the time spectrum, and calculate the average frequency of the spectrum at time τ And the corresponding average fluid velocity Based on the calculation results, the mean, variance, and average zero-crossing rate are extracted from the average velocity fluctuation time series to quantify the central tendency, fluctuation degree, and morphological change intensity of the data, respectively, as time domain features for flow pattern characterization. (4) The Doppler shift data is decomposed into the sum of several intrinsic mode functions and residual terms using the ensemble empirical mode decomposition method; based on the decomposition results, the correlation coefficients between the eigenmode functions at each level and the original Doppler shift data are calculated, and the first several eigenmode function components with the largest correlation coefficients are selected as the principal components representing the significant characteristics of the Doppler shift data, and the energy proportion corresponding to the principal components is used as the multi-scale domain feature of the flow pattern characterization; (5) forming a multidimensional feature vector based on pulse wave and continuous wave ultrasonic test data; (6) Taking the multidimensional feature vector as the input data set, the dimensionality is reduced using the isometric feature mapping method to obtain the reduced-dimensional feature vector; (7) Add the oil-gas-water three-phase flow pattern label to the feature vector after dimensionality reduction, and form a data set with the feature vectors under different working conditions; (8) Divide the data set into a training data set and a test data set according to a certain ratio. The training data set is used to train the flow type classifier model, and the test data set is used to evaluate the classification error of the trained classifier model. (9) Obtain a nonlinear support vector machine flow type classifier for outputting the flow type of unknown samples.

2. The method for identifying the flow pattern of oil-gas-water three-phase flow in a horizontal pipeline according to claim 1 is characterized in that: In step (1), the pulse wave ultrasonic sensor transmits ultrasonic pulses to the fluid intermittently at a certain repetition frequency, and receives echoes during the pulse interval. By demodulating the amplitude of the received echoes, the echo intensity profile data A(y,t) reflecting the distribution characteristics of the fluid interface is obtained. PRF ), where t PRF ∈(0,nT PRF ) is the measurement time, T PRF is the pulse repetition period, n is the number of pulse excitations within the sampling time, y∈(0,D) is the distance from the bottom of the pipe, and D is the pipe diameter.

3. The method for identifying the flow pattern of oil-gas-water three-phase flow in a horizontal pipeline according to claim 2, characterized in that: In step (1), the continuous wave ultrasonic sensor comprises two piezoelectric ceramic chips installed at an angle θ to the horizontal direction, which are responsible for continuously transmitting and receiving ultrasonic waves into the fluid during operation; by frequency demodulating the received ultrasonic waves, the Doppler frequency shift data f reflecting the fluid velocity fluctuation characteristics is obtained. d (t), where t∈(0,T) is time and T is the sampling time.

4. The method for identifying the flow pattern of oil-gas-water three-phase flow in a horizontal pipeline according to claim 3 is characterized in that: 45°≤θ≤60°.

5. The method for identifying flow patterns of oil-gas-water three-phase flow in a horizontal pipeline according to claim 3, characterized in that: The method of step (2) is: In the echo intensity profile data A(y,t PRF ),t PRF ∈(0,nT PRF ), locate the position where the maximum echo amplitude occurs in the profile obtained within each pulse repetition period, and calculate its probability distribution: Where P(y) represents the maximum echo amplitude A within the sampling time. max The probability of appearing in the position interval [y,y+Δy], Num(y≤A max ≤y+Δy) represents the maximum echo amplitude A within the sampling time max The number of occurrences in the position interval [y,y+Δy]; Based on the statistical results, the kurtosis P is extracted from the probability distribution graph. kur , skewness P ske , weighted average P wei The position P corresponding to the maximum probability lo , which are used to quantify the flatness, symmetry, fluctuation distribution and extreme value of the probability distribution, as the spatial domain characteristics of the flow pattern: P lo =y, when P=P max Among them, E(*) represents the average operation, P max is the maximum probability value in the probability distribution, P mean and P var are the mean and variance of the probability distribution statistics, respectively.

6. The method for identifying the flow pattern of oil-gas-water three-phase flow in a horizontal pipeline according to claim 5, characterized in that: The method of step (3) is: Doppler frequency shift data f d (t),t∈(0,T) is short-time Fourier transform to obtain the time-frequency spectrum S d (τ,f), and for S d Each column of (τ,f) uses the following formula to calculate the average frequency of the spectrum at time τ And the corresponding average fluid velocity Where c is the speed of sound in the acoustic wedge, f0 is the frequency of the transmitted sound wave, and θ is the angle between the central axis of the piezoelectric ceramic chip and the horizontal direction; Based on the calculation results, the average flow velocity fluctuation time series Extract the mean variance Average zero-crossing rate They are used to quantify the central tendency, fluctuation degree and morphological change intensity of the data, respectively, as the time domain characteristics of flow pattern characterization: Among them, sgn(*) is the sign function, is the data after subtracting its own mean from the average flow velocity fluctuation time series, and M is the total number of data points in the average flow velocity fluctuation time series.

7. The method for identifying flow patterns of oil-gas-water three-phase flow in a horizontal pipeline according to claim 6, characterized in that: The method of step (4) is: The Doppler frequency shift data f is transformed into d (t),t∈(0,T) is decomposed into several intrinsic mode functions IMF j (t),j=1,2,...,N, and the residual term r N The sum of (t): Where N is the number of eigenmode functions obtained by decomposition; Based on the decomposition results, the correlation coefficients between the intrinsic mode functions of each level and the original Doppler frequency shift data are calculated, and the first b intrinsic mode function components IMF with the largest correlation coefficient are selected. i ,i∈{j1,j2,j3,...,j b } as the principal component to characterize the significant characteristics of Doppler frequency shift data, and the energy proportion R corresponding to the principal component i Multi-scale domain features as flow pattern representation:

8. The method for identifying flow patterns of oil-gas-water three-phase flow in a horizontal pipeline according to claim 7, characterized in that: In step (5), the multidimensional feature vector S is:

9. The method for identifying the flow pattern of oil-gas-water three-phase flow in a horizontal pipeline according to claim 8, characterized in that: In step (6), the multidimensional feature vector S is used as the input data set X = {x i |x1,x2,···,x b+7 }, where x i Represent each feature parameter in the multidimensional feature vector S, and use the isometric feature mapping method to reduce its dimension to m-dimensional space, m < b + 7, to obtain the feature vector after dimension reduction, which is achieved by the following sub-steps: (6.1) Use the K nearest neighbor search algorithm to determine the value of each data point x in the data set. i The k nearest neighbors of the , get its neighborhood space N i ={x i1 ,x i2 ,···,x ik },i∈(1,b+7); (6.2) Construct a neighborhood connection graph G = [V, E], where V represents the nodes in the neighborhood connection graph, which are all data points in the dataset X, and E represents the edges in the neighborhood connection graph, whose length is the distance d between each data point in the dataset X. G Measure: When data x i with x j When they belong to the same neighborhood space, d G Let the Euclidean distance between two points be d G =d2(x i ,x j )=||x i -x j ||2, when data x i with x j When they do not belong to the same neighborhood space, d G =∞, indicating no connection; (6.3) Use the shortest path algorithm to determine the shortest distance between any two points in the neighborhood connection graph G and obtain the geodesic matrix: D G ={dist(x i ,x j )},x i ,x j ∈X, where (6.4) Constructing low-dimensional embedding: Using the Multidimensional Scaling (MDS) algorithm, construct the kernel matrix Among them S G =[dist(x i ,x j )] 2 , H=IJ / (b+7), I is the unit diagonal matrix, J is the full 1 matrix, and τ(D G ) performs singular value decomposition and takes the first m eigenvalues ​​λ1≤λ2≤...≤λ m and the corresponding eigenvectors v1,v2,...,v m , then the new coordinates of the original data set X formed by the multidimensional feature vector S in the low-dimensional m-dimensional space are expressed as: Among them, Y is the eigenvector after dimensionality reduction, Λ m =diag(λ1,λ2,...,λ m ), V m =[v1,v2,...,v m ].

10. The method for identifying flow patterns of oil-gas-water three-phase flow in a horizontal pipeline according to claim 9, characterized in that: The shortest path algorithm is Dijkstra's or Floyds' algorithm.

Images