A method for measuring real-time flow of industrial pipeline based on MUSIC algorithm

Abstract

The application discloses a kind of industrial pipeline real-time flow measurement methods based on MUSIC algorithm, step 1, through the multiple sensor array of distribution in the outer wall of pipeline, fluid flow signal in pipeline is collected;Step 2, the collected signal is preprocessed, including window function processing and band-pass filtering;Step 3, the robust covariance matrix of array measurement signal is constructed;Step 4, the robust covariance matrix is processed, and the noise subspace characteristic vector matrix Un is obtained;Step 5, calculate the MUSIC spectrum, construct the power spectrum function, calculate k-f spectrum diagram;Step 6, after drawing energy curve, i.e.wave number-frequency diagram, the optimal flow rate is searched in k-f spectrum diagram by tilting superposition method, to calculate and output real-time flow rate.The application greatly improves the signal source estimation accuracy and measurement precision under low signal-to-noise ratio.

Description

Technical Field

[0001] This invention belongs to the field of pipeline flow monitoring for industrial users, and relates to a method for real-time flow measurement of industrial pipelines based on an external clamp-on sensor array and an array processing algorithm. Background Technology

[0002] Flow rate is a key parameter characterizing the operating state of fluids and the level of energy transfer in a system. Since the acceleration of industrialization in the 21st century, the scale and complexity of fluid transport in energy production, industrial manufacturing, water resource utilization, and environmental governance have been continuously increasing. Currently, flow measurement is widely used in important fields such as industrial process control, energy trade settlement, water conservancy scheduling, safety production, and environmental monitoring, and its measurement results are directly related to production efficiency, economic benefits, and operational safety. As my country's industrial system develops towards large-scale, continuous, and intelligent operations, the requirements for the accuracy, reliability, and long-term stability of flow measurement are increasing.

[0003] Sonar flowmeters, as a non-invasive flow measurement technology, demonstrate unique advantages in large-diameter pipes and complex operating conditions. Internationally, sonar flow measurement methods based on multi-sensor arrays utilize theories such as beamforming, tomography, and inverse problem solving to passively acquire turbulent vibration signals using an array of sensors placed outside the pipe, and then estimate flow velocity using covariance matrix analysis and eigenvalue decomposition. Driven by this, China has also made some progress in the engineering application of non-invasive flowmeters.

[0004] However, in real-world industrial scenarios, existing sensor array-based flow velocity measurement methods (especially traditional MUSIC algorithms) face the following serious technical bottlenecks when practically implemented: (1) Ill-conditioned and time-non-stationary covariance matrix: Due to factors such as limited sampling time, high correlation of sensor channels and concentrated energy of turbulent signals, the constructed covariance matrix often exhibits ill-conditioned characteristics or has a rank close to 1; at the same time, the fluid flow state changes rapidly over time, resulting in extremely poor numerical stability of the eigenvalue decomposition and inversion process. (2) Strong mechanical noise interference: Industrial sites are often accompanied by strong low-frequency mechanical noise such as pump and valve vibration. These strong mechanical interferences usually mask weak fluid turbulence signals, causing subsequent spatial spectrum estimation to fail. (3) Failure of signal source quantity estimation under low signal-to-noise ratio: Traditional MUSIC algorithms rely heavily on accurate signal source quantity estimation. However, under harsh conditions of low signal-to-noise ratio (SNR) or small sample size, traditional information theory-based estimation algorithms (such as AIC or MDL criteria) often make serious misjudgments (overestimation or underestimation), resulting in the mixing of signal subspace and noise subspace, a sharp decline in algorithm performance, and even the generation of false flow velocity peaks; (4) Insufficient real-time processing capability: In order to improve robustness, existing methods usually introduce complex regularization or multi-dimensional matrix global operations, which significantly increases the computational load of embedded platforms such as DSP or MCU, making it difficult to meet the high real-time and low power consumption requirements of dynamic flow field monitoring.

[0005] Therefore, how to achieve efficient, stable, and accurate effective signal source estimation under complex flow conditions with strong interference and ill-conditioned covariance matrix, and ultimately complete real-time flow measurement suitable for embedded engineering implementation, remains a key technical problem that urgently needs to be solved in the field of industrial flow monitoring. Summary of the Invention

[0006] This invention provides a real-time flow measurement method for industrial pipelines based on the MUSIC (Multiple Signal Classification) algorithm to solve the problems of poor flow velocity stability, large computational load, large error, and low accuracy in current array flow measurement. MUSIC is a high-precision algorithm used in array signal processing to extract parameters from multiple signal sources. Its core is based on the eigenspace decomposition of the covariance matrix of the array received data.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for measuring fluid velocity in industrial pipelines based on the MUSIC algorithm includes the following steps: Step 1: Collect flow signals caused by fluid turbulence in the pipe by using a multi-sensor array that is equally spaced along the pipe axis; Step 2: Preprocess the acquired signal, including bandpass filtering based on window function processing; Step 3: Use Fast Fourier Transform to convert the preprocessed data to the frequency domain to obtain X(n), and combine it with the time domain to construct the robust covariance matrix of the array measurement signal through first-order exponential smoothing update; ; In the formula, For a robust covariance matrix, X is the frequency domain signal matrix, n is the frequency index of the signal, and a piezoelectric vibration sensor is selected with a frequency range of 0.1-20kHz and a sensitivity requirement of greater than 10pC / N. Step 4, for Explicit stripping is performed to obtain the covariance matrix R after explicit stripping. new Then for R new The eigenvalue decomposition of the matrix yields the eigenvector matrix of the noise subspace. ; Step 5: Calculate the MUSIC spatial power spectrum value based on the noise subspace eigenvector matrix and the array steering vector. ; In the formula, The power spectrum value of the MUSIC space, a(ω) i ,k j ) is the array steering vector. Because Pmusic is a two-dimensional variable, ω is needed. i Provide indexes; Where x is a vector consisting of the installation locations of multiple sensors; X represents the frequency domain signal, and x represents the installation location. ω i k is the angular frequency. j Wave number; k j The value range is [-103.85, 103.85] rad / m, ω i The value range is [-0, 314] rad / s. The steering vector reflects the phase delay relationship when a signal with a specific wavenumber arrives at each sensor array in space. The wavenumber is used to characterize the phase change per unit length when the signal propagates along the axial direction of the pipe. Its specific range is determined by the sensor strip spacing d, which is [ , ]; Based on the MUSIC spatial power spectrum values, a two-dimensional spectrum of wavenumber-frequency (kf) can be plotted; w and k are the independent variables (coordinate axes) for constructing the two-dimensional spectrum. Step 6: Plot the energy curve, i.e., the kf spectrum (wavenumber-frequency diagram). By using the tilt stacking method, perform energy integration along the fluid dispersion relation line in the kf spectrum to search for the optimal flow velocity, thereby calculating and continuously outputting the real-time flow velocity.

[0008] The multi-sensor array in step 1 is a non-invasive piezoelectric sensor array, with the mounting section of each sensing unit perpendicular to the fluid flow direction within the pipe. Preferably, it is a PVDF piezoelectric array containing at least four sensing units, but it is not limited to any specific piezoelectric sensor array.

[0009] In step 2, the Hamming window is used as the window function, and the FIR filter is used for the bandpass filter; the allowed frequency range for the bandpass filter is [10, 500] Hz.

[0010] In step 4, the following formula is used to... Perform explicit stripping: ; In the formula, R newThe covariance matrix is ​​the result of explicit noise stripping (i.e., stripping away strong interference); lam1 is the largest mechanical noise eigenvalue, and v1 is the largest mechanical noise eigenvector; the power method is used to quickly calculate the largest mechanical noise eigenvalue lam1 and vector v1, which together form the covariance matrix. The maximum eigenvalue and its corresponding principal eigenvector are obtained by using the power method, a mature and classic algorithm for calculating the principal eigenvalue and eigenvector of a matrix. Its advantages include low computational cost and the ability to extract the strongest mechanical noise component without full decomposition. Then, for R new Eigenvalues ​​are obtained using Matrix Eigenvalue Decomposition (EVD) and sorted in descending order. Adaptive fusion estimation is performed using the AIC or MDL information theory criterion combined with the eigenvalue ratio threshold method (the threshold is preferably 0.01) to determine the number of effective signal sources K, where K is the number of signal sources exceeding the threshold, and the value of K is generally 5. The eigenvectors corresponding to the top K largest eigenvalues ​​are classified into the "signal subspace", and the eigenvectors corresponding to the remaining MK small eigenvalues ​​are classified into the "noise subspace", thus separating the eigenvector matrix of the signal subspace and the eigenvector matrix of the noise subspace, where M is the total number of vibration sensors. ; In the formula, Let U be the covariance matrix after explicit stripping, U be the eigenvector matrix, Σ be a diagonal matrix whose diagonal elements are eigenvalues, Us be the eigenvector matrix of the signal subspace, and Un be the eigenvector matrix of the noise subspace.

[0011] In step 6, for a given trial flow velocity u, with a value ranging from [0.5, 7] m / s, along the dispersion relation line in the k−f spectrum... Perform energy integration: ; Here, E(u) is the cumulative energy corresponding to the flow velocity u. The optimal flow velocity estimate is determined by finding the maximum value of E(u) using a peak-finding algorithm, and the minimum frequency value f is obtained. min The frequency is 10Hz, and the maximum frequency f is... max 500Hz; The peak-finding algorithm is used to find the u value corresponding to the highest energy point on the energy curve (i.e., the flow velocity u corresponding to the maximum value of E(u)), which is the best estimated flow velocity at the current moment. The peak-finding algorithm adopts a two-stage peak-finding mode of coarse search + fine search. The coarse search calculates E(u) with a fixed step size in the range of [0.5, 7] m / s to find the discrete point with the maximum energy. The fine search uses cubic spline interpolation or parabola to find the extreme point within ±10% of the maximum point in order to obtain the sub-step velocity resolution.

[0012] In step 4, QR decomposition is performed on the stripped matrix: ; In the formula, R new Q is the covariance matrix after explicit stripping. S Q is a submatrix formed by the first K columns of an orthogonal matrix Q. N R is the submatrix formed by the remaining MK columns of the orthogonal matrix Q. 11 R 12 R 22 These are the upper triangular block matrices divided according to their respective dimensions.

[0013] Beneficial effects: The advantages and innovations of this invention are as follows: 1. The improved MUSIC algorithm proposed in this invention greatly enhances the accuracy of signal source estimation and measurement precision under low signal-to-noise ratio: In response to the fatal flaw of weak fluid turbulence signals and low signal-to-noise ratio in actual pipeline measurements, which causes the traditional AIC / MDL criteria to fail, this invention innovatively integrates eigenvalue spectrum morphology (proportional threshold) and information theory criteria, realizing adaptive and accurate estimation of the number of effective signal sources under harsh working conditions, avoiding spectral peak aliasing and false flow velocity peaks caused by overestimation or underestimation of the number of sources; 2. This algorithm achieves efficient removal of strong mechanical interference while reducing its complexity: by introducing power-law iteration, it quickly and accurately captures and explicitly removes the main mechanical noise features, such as pump and valve vibrations. This operation eliminates the need for complete full-space feature decomposition of the original large matrix, significantly reducing embedded computational overhead while filtering out strong interference. Combined with a piezoelectric array sensor from a passive sonar system and employing sliding window analysis, this algorithm can continuously monitor flow velocity changes and measure pipeline flow velocity in real time. 3. The algorithm enhances the matrix robustness in dynamic non-stationary flow fields: by combining the first-order exponential smoothing mechanism in the time domain with the adaptive diagonal loading strategy, it effectively overcomes the ill-conditioned inversion problem of the covariance matrix under finite sampling snapshots and small samples, and ensures the numerical stability of the eigenvalue decomposition. 4. Image-level enhancement of weak ridges and high-precision optimization: After conventional MUSIC spatial spectrum calculation, a two-dimensional Gaussian smoothing kernel is introduced for the first time to perform spatial energy denoising on the kf spectrum, filtering out shot noise. Combined with a two-step method of coarse search and spline interpolation fine search, weak convective energy ridges are highlighted, thereby achieving flow velocity measurement accuracy far exceeding that of traditional grid search. Attached Figure Description

[0014] Figure 1 This is an overall flowchart of the flow velocity measurement method based on the improved MUSIC algorithm described in this invention.

[0015] Figure 2This is a schematic diagram of the multi-sensor array arrangement in an embodiment of the present invention.

[0016] Figure 3 This is a schematic diagram of the overall experimental platform in an embodiment of the present invention.

[0017] Figure 4 The wavenumber-frequency (kf) two-dimensional spatial spectrum is calculated and smoothed using the embodiments of the present invention.

[0018] Figure 5 The diagram shows the energy-velocity integral curve and the optimal velocity result obtained by the tilt superposition method and interpolation fine optimization in this embodiment of the invention.

[0019] Figure 6 This is a graph showing the long-term measurement results of the present invention under constant flow rate. Detailed Implementation

[0020] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Example 1

[0021] Taking the measurement of flow velocity in an industrial pipeline using an 8-channel PVDF piezoelectric sensor array as an example, such as Figure 2 As shown, eight sensor patches are arranged at equal intervals along the pipe axis (e.g., a spacing of 0.0605m). The sensor strips are directly attached to the outer wall of the pipe and fixed by threaded holes at multiple locations to ensure phase consistency of the signals in each channel. The overall experimental platform is as follows. Figure 3 As shown, the overall system is driven by a variable frequency water pump to achieve circulating flow between the water tank and the pipeline, providing a stable flow velocity range of 0–3 m / s. An electromagnetic flowmeter is installed to compare and verify the algorithm calculation results. Multiple sets of data from 0.8–3 m / s were collected, with the vibration signal sampling rate at 50000 Hz and each data set lasting 60 seconds. Implementing the measurement method proposed in this invention specifically includes the following steps: Step 1: Acquire time series data of pipe flow velocity measured by array sensor; Using the aforementioned 8-channel piezoelectric sensor array, at a set sampling rate (e.g., 50000Hz), the minute acoustic / vibration signals caused by the turbulent vortex flow of fluid in the pipe passing over the sensor surface are simultaneously acquired. Step 2, signal preprocessing, including window function processing and bandpass filtering; A sliding data window (e.g., a window length of 15 seconds and a step size of 2 seconds) is used to truncate the time series signal, and window function processing is performed sequentially to reduce spectral leakage. Subsequently, a bandpass filter (passband range set to [5,50]Hz) is used to extract effective analysis frequency band data containing the main convection information; Step 3: Construct a robust covariance matrix; The preprocessed data is converted to the frequency domain using Fast Fourier Transform (SFFT) to obtain X(n), and then the robust covariance matrix of the array measurement signal is constructed by first-order exponential smoothing update in the time domain. ; To overcome the data nonstationarity caused by rapid changes in the actual flow field, this invention employs a first-order exponential smoothing mechanism in the time domain to update the covariance matrix: ; In the formula, This is the smoothed robust covariance matrix at the current time. This is the covariance matrix initially calculated based on the frequency domain signal matrix at the current moment. The robust covariance matrix of the previous time step. This is a time smoothing factor (e.g., a value of 0.5). Step 4: Perform diagonal loading to explicitly remove strong mechanical interference, and finally perform matrix decomposition. First, calculate the eigenvalues ​​of the current covariance matrix and evaluate the signal-to-noise ratio (SNR) of the current snapshot signal. When the SNR is lower than the preset severe working condition threshold, adaptively calculate the diagonal loading amount and inject it into the main diagonal of the covariance matrix to solve the ill-conditioned matrix inversion problem. Secondly, explicit removal of strong mechanical interference is performed. For the low-frequency strong vibration of pumps and valves that may exist in the pipeline system, the power method is used to quickly calculate the maximum mechanical noise eigenvalue (lam1) and vector (v1), and then it is removed. ; In the formula, Let be the covariance matrix, lam1 be the largest mechanical noise eigenvalue, and v1 be the largest mechanical noise eigenvector; This operation avoids directly performing full-dimensional feature decomposition on large matrices containing strong noise, reducing the computational complexity of embedded hardware such as DSPs / MCUs. Finally, the stripped matrix is ​​decomposed using QR-EVD. To address the shortcomings of traditional AIC or MDL information theory criteria, which are prone to failure at low signal-to-noise ratios, this embodiment adopts a fusion strategy: first, the first estimated number of sources is selected based on the eigenvalue spectrum shape (e.g., setting the eigenvalue energy ratio threshold to 0.01); then, the second estimated number of sources is calculated using the AIC criterion; the two are then weighted and averaged or fused using the median to finally output an extremely stable number of effective signal sources, which is used as the basis for dividing a high-purity noise subspace. ; In the formula, Let U be the covariance matrix after explicit stripping, U be the eigenvector matrix, Σ be a diagonal matrix whose diagonal elements are eigenvalues, Us be the eigenvector matrix of the signal subspace, and Un be the eigenvector matrix of the noise subspace. Step 5: Calculate the MUSIC spectrum, construct the power spectral function, and calculate the kf spectrum; Calculate the MUSIC spectrum and construct the power spectrum function according to the formula. ; In the formula, U represents the spatial power spectrum value of the MUSIC system. n Let a(ω) be the eigenvector matrix of the noise subspace. i ,k j ) is the array steering vector, ω i k is the angular frequency. j Wave number; After calculating the basic two-dimensional wavenumber-frequency (kf) spectrum, as follows: Figure 4 As shown, this embodiment further introduces a two-dimensional Gaussian smoothing filter kernel from the field of computer vision to perform two-dimensional convolution processing on the two-dimensional matrix. This effectively smooths out discrete "spiking" noise on the spatial spectrum, making the energy ridges representing fluid convection velocities more continuous and prominent, as shown in the figure. Figure 4 As shown, this method can still generate clear "ridges" under low signal-to-noise ratio conditions; Step 6: Use tilted superposition integral search to find the optimal flow velocity; After plotting the energy curve (wavenumber-frequency plot), the optimal flow velocity is searched in the k-f spectrum using the tilt stacking method to calculate the real-time flow velocity. For a given trial flow velocity u, the velocity is calculated along the dispersion relation line in the k-f spectrum. Perform energy integration: ; Here, E(u) is the cumulative energy corresponding to the flow velocity u, and the optimal flow velocity estimate is determined by finding the maximum value of E(u).

[0022] like Figure 5 As shown, the flow velocity search range is set to [0.5, 7] m / s. The cumulative energy at each test flow velocity is calculated using the tilt stacking method. First, a coarse search of 400 points is performed to find the peak value, and then a fine search of 100 points is performed within ±10% of the peak value using spline interpolation. Finally, the highest point on the interpolation curve is output as the optimal flow velocity at the current moment (e.g., ...). Figure 5 (The location is marked by the middle circle).

[0023] Furthermore, when this method is run on embedded hardware (such as DSP / MCU), it only performs operations on extremely low-dimensional (8*8) matrices and avoids full-dimensional feature decomposition by explicitly stripping away strong interference through power-law iteration. Therefore, a single flow rate calculation takes less than one second. Combined with a sliding window with a set step size (e.g., 2 seconds), the system can achieve continuous, delay-free flow rate refresh, fully possessing industrial-grade real-time processing capabilities.

[0024] To verify the accuracy of the real-time measurement method proposed in this invention, the measurement output value of this method was synchronously compared with the calculation results of a high-precision electromagnetic flowmeter and the algorithm under the same pipe diameter and operating conditions. Multiple sets of different test flow rates were set, and the measurement data after the system stabilized were recorded. The relative error between the flow rate measured by the algorithm of this invention and the actual flow rate of the electromagnetic flowmeter was within ±2%.

[0025] ; Table 1 is a comparison table of electromagnetic flowmeter readings and flow velocity measurements using this algorithm. Table 1 shows that the improved MUSIC algorithm of this invention can achieve high-precision, high-real-time pipeline flow monitoring.

[0026] To evaluate the dynamic response characteristics and long-term stability of the algorithm under time-varying conditions, a sliding window processing was applied to continuous stable flow velocity data for 60 s. Considering the non-stationarity of the flow process and the need to establish statistical characteristics, the sliding step was 2 s. The flow velocity evolution over time and its statistical characteristics obtained by the two algorithms are shown in Figure 5. The measurement results remained stable throughout the observation period, with an average flow velocity of 2.874 m / s and a standard deviation of 0.040 m / s. No obvious divergence was observed, indicating that the algorithm can achieve long-term stable measurement and monitoring of flow velocity.

Claims

1. A method for real-time flow measurement of industrial pipelines based on the MUSIC algorithm, characterized in that, Includes the following steps: Step 1: Collect flow signals caused by fluid turbulence in the pipe by using a multi-sensor array that is equally spaced along the pipe axis; Step 2: Preprocess the acquired signal, including bandpass filtering based on window function processing; Step 3: Use Fast Fourier Transform to convert the preprocessed data to the frequency domain to obtain X(n), and combine it with the time domain to construct the robust covariance matrix of the array measurement signal through first-order exponential smoothing update; ； In the formula, Let X be the robust covariance matrix, X be the frequency domain signal matrix, and n be the frequency index of the signal. Step 4, for Explicit stripping is performed to obtain the covariance matrix R after explicit stripping. new Then for R new The eigenvalue decomposition of the matrix yields the eigenvector matrix of the noise subspace. ; Step 5: Calculate the MUSIC spatial power spectrum value based on the noise subspace eigenvector matrix and the array steering vector. ； In the formula, The power spectral density of the MUSIC space is a(ω) i ,k j ) is the array steering vector. ; Where x is a vector consisting of the installation locations of multiple sensors; ω i k is the angular frequency. j Wave number; Step 6: Plot the energy curve, i.e., the kf spectrum. Integrate the energy along the fluid dispersion relation line in the kf spectrum using the tilt stacking method to search for the optimal flow velocity, thereby calculating and continuously outputting the real-time flow velocity. In step 4, the following formula is used to... Perform explicit stripping: ； In the formula, R new is the covariance matrix after explicit stripping; lam1 is the largest mechanical noise eigenvalue, and v1 is the largest mechanical noise eigenvector; Then, for R new The eigenvalues ​​are obtained by matrix eigenvalue decomposition and sorted in descending order. Adaptive fusion estimation is performed by combining the AIC or MDL information theory criteria with the eigenvalue ratio threshold method to determine the number of effective signal sources K, where K is the number of signal sources exceeding the threshold. The eigenvectors corresponding to the top K largest eigenvalues ​​are classified into the "signal subspace", and the eigenvectors corresponding to the remaining MK small eigenvalues ​​are classified into the "noise subspace", thus separating the eigenvector matrix of the signal subspace and the eigenvector matrix of the noise subspace, where M is the total number of vibration sensors. ； In the formula, Let U be the covariance matrix after explicit stripping, U be the eigenvector matrix, Σ be a diagonal matrix whose diagonal elements are eigenvalues, Us be the eigenvector matrix of the signal subspace, and Un be the eigenvector matrix of the noise subspace.

2. The method for real-time flow measurement of industrial pipelines based on the MUSIC algorithm according to claim 1, characterized in that, The multi-sensor array in step 1 is a non-invasive piezoelectric sensor array, and the installation section of each sensing unit is perpendicular to the direction of fluid flow in the pipe.

3. The method for real-time flow measurement of industrial pipelines based on the MUSIC algorithm according to claim 1, characterized in that, In step 2, the Hamming window is used as the window function, and the FIR filter is used as the bandpass filter.

4. The method for real-time flow measurement of industrial pipelines based on the MUSIC algorithm according to any one of claims 1-3, characterized in that, In step 6, for a given trial flow velocity u, with a value ranging from [0.5, 7] m / s, along the dispersion relation line in the k−f spectrum... Perform energy integration: ； Here, E(u) is the cumulative energy corresponding to the flow velocity u. The optimal flow velocity estimate is determined by finding the maximum value of E(u) using a peak-finding algorithm, and the minimum frequency value f is obtained. min The frequency is 10Hz, and the maximum frequency f is... max 500Hz; The peak-finding algorithm is used to find the u-value corresponding to the highest energy point on the energy curve, which is the best estimated flow velocity at the current moment.

5. The method for real-time flow measurement of industrial pipelines based on the MUSIC algorithm according to claim 3, characterized in that, The peak-finding algorithm adopts a two-stage peak-finding mode of coarse search + fine search. The coarse search calculates E(u) with a fixed step size in the range of [0.5, 7] m / s to find the discrete point with the maximum energy. The fine search uses cubic spline interpolation or parabola to find the extreme point within ±10% of the maximum point in order to obtain the sub-step velocity resolution.

Images