Industrial pipeline flow velocity measurement method based on differential array and mvdr algorithm

Abstract

The application discloses an industrial pipeline flow velocity measurement method based on a differential array and an MVDR algorithm, wherein step 1 is that fluid flow signals in a pipeline are collected through a multi-sensor array distributed on the outer wall of the pipeline, and adjacent channel differential preprocessing is carried out to suppress common-mode noise; step 2 is that a frequency domain signal is constructed by using STFT, a time recursive smoothing mechanism is adopted to update a covariance matrix, and dynamic adaptability is enhanced; step 3 is that an improved MVDR algorithm is used to calculate a wave number-frequency (k-f) spectrum diagram in combination with adaptive diagonal loading and SVD regularization technology; step 4 is that background noise suppression and energy enhancement processing are carried out on the k-f two-dimensional power spectrum diagram, and effective turbulent flow characteristics are highlighted; and step 5 is that high-precision real-time measurement of fluid flow velocity is realized through a tilt superposition method in combination with a spline interpolation fine search strategy. The application effectively solves the problems of inaccurate covariance matrix estimation and difficult feature extraction under low signal-to-noise ratio in complex working conditions.

Description

Technical Field

[0001] This invention belongs to the field of pipeline fluid measurement and signal processing technology. Specifically, it relates to a method for non-invasive real-time measurement of fluid velocity and flow rate in industrial pipelines by using a sensor array distributed on the outer wall of a pipeline, combined with differential signal processing and an improved minimum variance distortionless response (MVDR) algorithm. This method is suitable for high-noise industrial environments. Background Technology

[0002] Pipeline flow velocity is a key parameter in industrial process control. Current measurement methods mainly include electromagnetic flow meters, ultrasonic flow meters, and orifice plate flow meters. Among them, the passive sonar measurement method based on the principle of turbulence-induced vibration has attracted much attention due to its advantages such as non-invasiveness, no pressure loss, and convenient installation. This method utilizes the turbulent pressure fluctuations (Tbl) generated when fluid flows through the pipe to convect along the flow direction, and calculates the flow velocity by analyzing the correlation of signals received by a sensor array on the outer wall of the pipe.

[0003] However, in practical industrial applications, flow velocity measurement based on traditional array signal processing methods (such as the traditional MVDR algorithm) has the following problems:

[0004] Severe common-mode environmental noise interference: Industrial pipelines are typically connected to mechanical equipment such as pumps and valves. The mechanical vibrations generated by these devices act as strong common-mode signals, affecting all sensors simultaneously. Traditional methods struggle to distinguish between turbulent signals and overall pipeline vibration noise, resulting in extremely low signal-to-noise ratios and measurement failures.

[0005] Covariance matrix estimation is unstable: Fluid flow is non-stationary, and to ensure real-time performance, the number of available data snapshots is usually limited. The sample covariance matrix is ​​often ill-conditioned or even irreversible, and direct inversion will produce huge numerical errors, causing MVDR beamformers to fail.

[0006] Weak signal feature extraction is difficult: In low-velocity or high-viscosity fluids, the vibration signals caused by turbulence are extremely weak and are often submerged in background noise. Traditional peak search algorithms have difficulty accurately locating velocity ridges in cluttered wavenumber-frequency (kf) spectra.

[0007] Therefore, there is an urgent need for a real-time flow velocity measurement method that can effectively suppress common-mode noise, maintain computational robustness under insufficient snapshot conditions, and enhance weak turbulence characteristics. Summary of the Invention

[0008] The purpose of this invention is to overcome the problems of low accuracy and poor algorithm stability in existing technologies for flow velocity measurement in complex industrial noise environments, and to provide a real-time flow velocity measurement method for industrial pipelines based on differential arrays and an improved MVDR algorithm.

[0009] To achieve the above objectives, the present invention provides the following technical solution: A method for measuring flow velocity in industrial pipelines based on differential array and MVDR algorithm, comprising: Step 1: Acquire vibration signals caused by fluid turbulence using a multi-sensor array; Multiple sensors in a multi-sensor array are arranged at equal intervals along the pipe axis to achieve non-invasive measurement of fluid velocity inside the pipe; The original sensor signal is , where i represents the i-th sensor; Step 2: After performing differential and bandpass filtering preprocessing on the vibration signal, the time domain signal is converted to the frequency domain using short-time Fourier transform (STFT) to obtain the frequency domain signal vector X(f) of the turbulent vibration. In step 2, when converting the time-domain signal to the frequency domain using the Short-Time Fourier Transform (STFT), overlapping frame processing is employed, and a window function is added to suppress spectral leakage. The window function used is the Hanning window, but other window functions are not limited to. The differential steering vector is: ; Among them, a i (k) represents the component of the differential steering vector in the i-th differential channel. The overall differential steering vector a(k) is composed of the steering components a corresponding to each differential channel. i (k) is arranged in channel order; k is the wave number, used to characterize the phase change per unit length when the signal propagates along the axial direction of the pipe, and j is the imaginary unit. The complex exponential steering vector is constructed based on the installation location information of the (i+1)th sensor. The complex exponential steering vector is constructed based on the installation location information of the i-th sensor; Step 3: Construct the covariance matrix of the sensor array and use a recursive weighted smoothing algorithm to update the covariance matrix in real time; The instantaneous sampling covariance matrix of the current time window at frequency f is: ; In the formula, T is the number of frames in the time window, X(f) is the frequency domain signal vector, and H represents the conjugate transpose operation; The recursive weighted smoothing algorithm is as follows: Introducing a smoothing factor ,use and the covariance matrix of the previous time step Recursively update the matrix to obtain the covariance matrix currently used for calculation: This recursive mechanism is used to enhance the robustness of the covariance matrix under dynamic flow velocity environments. Step 4: Use the MVDR algorithm to calculate the power spectrum values ​​corresponding to each frequency point and wavenumber point, and construct a two-dimensional power spectrum diagram of wavenumber-frequency (kf two-dimensional power spectrum diagram). The MVDR algorithm includes adaptive diagonal loading based on signal-to-noise ratio estimation and robust inverse matrix solving based on singular value decomposition (SVD). Adaptive diagonal loading based on signal-to-noise ratio estimation refers to: performing eigenvalue decomposition on the covariance matrix, estimating the signal-to-noise ratio (SNR) at the current frequency point based on the eigenvalue distribution, and dynamically adjusting the diagonal loading amount δ accordingly to obtain a regularized matrix. Forming a regularized matrix ; in, The trace of the matrix, It is the identity matrix. Diagonal loading amount; It is usually taken as 0.1; Robust inverse matrix solving based on singular value decomposition (SVD) refers to: right Decompose to obtain ,in, It is a left singular vector matrix. It is a right singular vector matrix. For a singular value diagonal matrix, for a singular value diagonal matrix The small singular values ​​in the matrix are truncated or regularized to obtain the regularized eigenvalue inverse matrix. The regularization process refers to: excluding tiny singular values ​​smaller than a preset threshold from the inversion operation, or restricting their inversion values, to avoid numerical instability caused by excessively small singular values ​​during inversion, thereby improving the robustness of the inverse matrix calculation; and subsequently reconstructing the robust inverse matrix. Truncating the minute singular values ​​in a singular value diagonal matrix means removing the minute singular values ​​from the singular value diagonal matrix. To ensure the numerical stability of the inverse matrix calculation, singular values ​​smaller than 1 / 1000 of the maximum singular value are considered small singular values. ; Then calculate the MVDR power spectrum: ; The two-dimensional kf power spectrum is a two-dimensional representation of the MVDR power spectrum; In step 4, "MVDR power spectrum" refers to the spectral values ​​calculated using the MVDR algorithm at various wavenumber and frequency points; "kf two-dimensional power spectrum" refers to the two-dimensional spectral distribution map formed by arranging the above MVDR power spectrum values ​​according to two-dimensional coordinates of wavenumber and frequency. Therefore, the kf two-dimensional power spectrum is the two-dimensional representation of the MVDR power spectrum. In the main text, it can be consistently referred to as "MVDR power spectrum" and "kf two-dimensional power spectrum". Step 5: Perform noise reduction and enhancement processing on the kf two-dimensional power spectrum; This includes employing Gaussian smoothing to suppress local random noise, and normalization based on background noise floor to enhance effective turbulence features; Step 6: For the enhanced kf two-dimensional power spectrum, energy integration is performed using the tilt superposition method, and the optimal fluid velocity is calculated by combining coarse search and interpolation fine search.

[0010] The signal-to-noise ratio normalization based on background noise floor mentioned in step 5 refers to: estimating the background noise floor in the non-signal region for the frequency row of the kf two-dimensional power spectrum, and subtracting the floor value from the spectrum intensity to suppress background noise floor and retain effective turbulence features with high signal-to-noise ratio.

[0011] The smoothing process described in step 5 uses two-dimensional Gaussian filtering to eliminate local random noise in the kf two-dimensional power spectrum.

[0012] The energy integration using the tilted superposition method described in step 6 refers to: for a given trial flow velocity u, along the dispersion relation line in the kf two-dimensional power spectrum. Perform energy integration: ; Where E(u) is the cumulative energy corresponding to the flow velocity u. The enhanced two-dimensional energy spectrum; The frequency variable corresponding to the short-time Fourier transform. and These are the lower and upper limits for frequency integration, preferably determined based on the frequency band where the effective turbulence characteristics are located; in this embodiment, and They are 5 Hz and 50 Hz respectively; The optimal flow velocity estimate is determined by finding the maximum value of the velocity energy function E(u).

[0013] By employing a two-stage strategy of "coarse search for peak values ​​+ spline interpolation for fine search," high-precision and high-resolution solutions for fluid velocity are achieved, ultimately yielding the optimal velocity. The specific steps are as follows: First, a coarse search is performed within the global velocity range to obtain a rough velocity estimate based on the peak position of the velocity-energy function E(u), where E(u) is a function obtained by integrating the enhanced two-dimensional energy spectrum along the dispersion relation. Then, a local search interval is constructed near the coarse velocity estimate, preferably within ±10% of the coarse velocity estimate, and spline interpolation is performed on the energy curve, where the energy curve is the curve formed by the velocity-energy function E(u) changing with the trial velocity u. The discrete energy points are fitted into a continuous smooth curve, and the peak value is found on this continuous curve to achieve a fine search, thereby obtaining the final optimal velocity.

[0014] In step 1, a multi-channel sensor array is arranged along the pipe axis to collect vibration signals and differential pairs between adjacent channels are constructed. Physical spatial differentiation is used to cancel the overall common-mode vibration noise of the pipe, extracting pure turbulent flow characteristic signals. These characteristic signals reflect the structural vibration signals generated by the pipe wall under turbulent pulsating pressure. These signals originate from the pressure pulsations in the fluid turbulence acting on the pipe wall, causing minute elastic vibrations. The sensors measure the acceleration or strain response signals of these vibrations at the sensor location. Therefore, the differential signal essentially characterizes the pipe wall vibration response characteristics under turbulent pressure pulsation excitation. Piezoelectric vibration sensors are selected, with a frequency range of 0.1-20kHz and a sensitivity requirement greater than 10 pC / N. In step 2, bandpass filtering and short-time Fourier transform are performed on the turbulent flow characteristic signal (i.e., the differential signal) to obtain the complex spectrum representation of the signal from each sensor channel in the time-frequency domain, which describes the frequency amplitude and phase distribution of the signal within different time windows. . The time spectrum reflects the frequency structure of the pipe wall vibration signal under turbulent pressure pulsation excitation and its time-varying characteristics, and can be used to extract the frequency domain correlation between different channels. The filtering process is mainly used to filter out low-frequency structural vibration, power frequency electromagnetic interference and its harmonic components, as well as non-flow excitation signals such as high-frequency electronic noise. Among them, low-frequency components usually come from the overall pipe structure vibration, equipment operation, or environmental disturbances in the frequency range of 5-10Hz; power frequency and its harmonics are mainly generated by the power supply system and motor equipment; high-frequency components are mostly sensor electronic noise or acquisition system noise, with high-frequency noise above 500Hz.

[0015] In step 3, through calculation The covariance yields the instantaneous covariance matrix at the current moment. (Right now (This is the instantaneous covariance matrix calculated from the STFT results of the current time frame). ; Where T is the number of data segments in the STFT, which is generally taken as 5 to 15, with a preferred value of 10.

[0016] When constructing the frequency domain covariance matrix, a time-recursive forgetting factor mechanism is introduced: ; in, This is the frequency domain covariance matrix updated at the current time. Forgetting factor, satisfying 0 < <1, The covariance matrix of the previous time step The initial value is a matrix of all zeros, and subsequently... The value is typically 0.5.

[0017] Based on this, a frequency domain covariance matrix is ​​constructed using the frequency domain spectral vectors under each time window as samples. A time-recursive forgetting factor mechanism is introduced to apply exponentially decaying weights to historical samples, enabling the covariance matrix to better track the time-varying characteristics of turbulence signals while ensuring statistical stability. The statistical characteristics of the previous time step are used to smooth the current estimate, solving the problem of covariance matrix fluctuation caused by insufficient short-term snapshots.

[0018] In step 4, the improved MVDR algorithm is used to calculate the wavenumber-frequency (kf) spectrum. The frequency domain covariance matrix is ​​updated at the current time. In the inversion process, an adaptive diagonal loading method is used to improve the matrix condition number. Specifically, the trace of the current covariance matrix is ​​calculated, a small loading amount is adaptively determined, and this loading amount is added to the diagonal elements of the matrix. By adding appropriate positive values ​​to the diagonal, the minimum eigenvalue of the matrix can be effectively increased, thereby improving the matrix condition number and reducing the sensitivity of the inversion process to noise and numerical errors. ; in, The trace of the matrix, It is the identity matrix. Typically, the value is set to 0.1. Then, singular value decomposition (SVD) is performed to invert the singular value matrix. ; ; Subsequently, by deleting tiny singularities, we obtained... To ensure the numerical stability of the inverse matrix calculation, small singular values ​​are determined by whether they are less than 1 / 1000 of the largest singular value. Subsequent calculations The inverse matrix is ​​obtained And then calculate ; in, To represent the estimation of the inverse of the covariance matrix after numerical stabilization, This is the right singular vector matrix obtained from singular value decomposition. This is the left singular vector matrix obtained from singular value decomposition. Let be the singular value inverse matrix after small singular value truncation. This is the conjugate transpose.

[0019] The differential steering vector is constructed as follows: ; in, Let be the physical installation position of the i-th sensor along the axial direction of the pipe.

[0020] In step 5, a background noise estimation algorithm is used on the generated original kf two-dimensional power spectrum to subtract the background noise floor from the energy of each frequency row of the spectrum, significantly improving the contrast of turbulent ridges. First, the corresponding noise floor is estimated: ; Median represents the median of the index set.

[0021] Subsequently, the background noise base was extracted from the original spectral values ​​to obtain... ; In the formula, P(i,j) represents the spectral energy value corresponding to the i-th frequency point and the j-th wavenumber point in the original calculated two-dimensional spatial spectral matrix; P en (i, j) represents the enhanced spectral energy value (i.e., enhanced spectral value) that retains the true signal-to-noise ratio intensity after background noise removal; N(i) represents the estimated background noise floor at the i-th frequency point; K outer K is a set of wavenumber indices for the pre-defined background noise region outside the spectrum. outer =0 indicates that there is no available outer region, K outer The value ≠0 indicates the existence of an "outer region," and the energy of the outer region can be used to statistically analyze the background noise.

[0022] Beneficial effects: The advantages and innovations of this invention are as follows: (1) Strong anti-interference capability: The combination of physical difference and mathematical difference guided vector effectively filters out common-mode noise such as pump vibration, and improves availability in low signal-to-noise ratio environments; (2) The algorithm is highly robust: the recursive smooth update combined with SVD regularization inversion completely solves the problems of matrix singularity and computational divergence in traditional MVDR when there are insufficient snapshots or strong signal correlation. (3) High measurement accuracy: Spectrum background suppression and interpolation fine search enable the system to lock in minute flow velocity changes or weak signal features, and achieve sensitive capture of flow field changes; (4) Strong real-time processing capability: Sliding window analysis enables continuous monitoring and meets the requirements of dynamic flow field measurement. Attached Figure Description

[0023] Figure 1 This is a flowchart illustrating the overall processing of the method proposed in this invention.

[0024] Figure 2 This is a schematic diagram of the sensor array arrangement and the differential principle of adjacent channels;

[0025] Figure 3 Comparison of the effects before and after processing the kf two-dimensional power spectrum (Figure) Figure 3 (a) is the traditional MVDR spectrum. Figure 3 (b) is the enhanced spectrum after background noise suppression.

[0026] Figure 4 The image shows the energy-velocity integral curve and velocity result obtained based on the tilt stacking method. Detailed Implementation

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

[0028] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0029] This embodiment takes the flow velocity measurement of a section of industrial pipeline as an example, and uses an 8-channel piezoelectric sensor array with a sensor spacing d=0.0605m and a system sampling frequency Fs=50kHz.

[0030] Step 1: The system synchronously acquires vibration signals from 8 channels. To suppress common-mode interference, adjacent channels are differentially analyzed, and the differences between adjacent channels are calculated to obtain 7 sets of differential signals.

[0031] Step 2: Perform bandpass filtering on the differential signal, setting the passband range to [5,50]Hz, and then perform Fourier transform to convert the time-domain data into a frequency-domain snapshot matrix X(f).

[0032] Step 3: Construct the covariance matrix. ; To adapt to the time-varying characteristics of the flow velocity, a recursive smoothing algorithm is used to update the matrix: ; This step effectively utilizes historical information, allowing the covariance matrix to remain stable even with insufficient snapshots.

[0033] Step 4: Robustness processing is performed. First, adaptive diagonal loading is applied to calculate the eigenvalues ​​of R and estimate the signal-to-noise ratio. To prevent matrix singularities from causing inversion failure, singular value decomposition (SVD) is used. After decomposing R, the small singular values ​​are regularized and their reciprocals are used to reconstruct the inverse matrix. ; Based on the differential geometry, the guiding vector is constructed as follows: ; Finally, the MVDR spectrum was calculated: ; Step 5: After obtaining the original kf spectrum, perform SNR normalization enhancement operation to subtract the background noise value from the frequency row and highlight the energy of the turbulent ridge.

[0034] Finally, the flow velocity search range was set to [0.5, 8] m / s. The cumulative energy at each trial flow velocity was calculated using the tilt stacking method. A coarse search of 400 points was first performed to find the peak value, and then a fine search of 100 points was performed within ±10% of the peak value using spline interpolation. The optimal flow velocity was then finally output.

[0035] Experimental results are as follows Figure 3 As shown, this method can still generate clear "ridges" under low signal-to-noise ratio conditions. The ridges correspond to highly concentrated energy trajectories, reflecting the dominant fluid velocity component in the signal. The clearer the ridges, the more significant the velocity component and the stronger the algorithm's ability to suppress noise.

[0036] The algorithm's calculation results are highly consistent with the measured values ​​of the high-precision electromagnetic flowmeter, as shown in Table 1: ; Simultaneously, cross-calculation was performed using different channel and array combinations to verify the spatial consistency of flow velocity; the variation of flow velocity with flow rate and pump speed was analyzed, and the results conformed to theoretical physical characteristics; recovery tests were conducted under simulated flow velocity signals, showing that the algorithm can accurately reconstruct the flow velocity, comprehensively verifying the reliability of the measurement results from multiple perspectives. The real-time flow velocity tracking curve is smooth and stable, effectively solving the problem of covariance matrix estimation divergence in traditional methods under complex working conditions.

Claims

1. A method for measuring flow velocity in industrial pipelines based on differential array and MVDR algorithm, characterized in that, include: Step 1: Acquire vibration signals caused by fluid turbulence using a multi-sensor array; Multiple sensors in a multi-sensor array are arranged at equal intervals along the pipe axis to achieve non-invasive measurement of fluid velocity inside the pipe; The original sensor signal is , where i represents the i-th sensor; Step 2: After performing differential and bandpass filtering preprocessing on the vibration signal, the time domain signal is converted to the frequency domain using short-time Fourier transform to obtain the frequency domain signal vector X(f) of the turbulent vibration. The differential steering vector is: ； Among them, a i (k) represents the component of the differential steering vector in the i-th differential channel. The overall differential steering vector a(k) is composed of the steering components a corresponding to each differential channel. i (k) is arranged in channel order; k is the wave number, used to characterize the phase change per unit length when the signal propagates along the axial direction of the pipe, and j is the imaginary unit. The complex exponential steering vector is constructed based on the installation location information of the (i+1)th sensor. The complex exponential steering vector is constructed based on the installation location information of the i-th sensor; Step 3: Construct the covariance matrix of the sensor array and use a recursive weighted smoothing algorithm to update the covariance matrix in real time; The instantaneous sampling covariance matrix of the current time window at frequency f is: ； In the formula, T is the number of frames in the time window, X(f) is the frequency domain signal vector, and H represents the conjugate transpose operation; The recursive weighted smoothing algorithm is as follows: Introduce a smoothing factor α, 0 < α < 1, and utilize... and the covariance matrix of the previous time step Recursively update the matrix to obtain the covariance matrix currently used for calculation: ； Step 4: Use the MVDR algorithm to calculate the power spectrum values ​​corresponding to each frequency point and wavenumber point, and construct a wavenumber-frequency two-dimensional power spectrum map kf two-dimensional power spectrum map; The MVDR algorithm includes adaptive diagonal loading based on signal-to-noise ratio estimation and robust inverse matrix solving based on singular value decomposition; Adaptive diagonal loading based on signal-to-noise ratio estimation refers to: forming a regularized matrix. ; in, The trace of the matrix, It is the identity matrix. Diagonal loading amount; Robust inverse matrix solving based on singular value decomposition refers to: right Decompose to obtain ,in, It is a left singular vector matrix. It is a right singular vector matrix. For a singular value diagonal matrix, for a singular value diagonal matrix The small singular values ​​in the matrix are truncated or regularized to obtain the regularized eigenvalue inverse matrix. The regularization process refers to: excluding tiny singular values ​​smaller than a preset threshold from the inversion operation, or restricting their inversion values, to avoid numerical instability caused by excessively small singular values ​​during inversion, thereby improving the robustness of the inverse matrix calculation; and subsequently reconstructing the robust inverse matrix. : ； Then calculate the MVDR power spectrum: ； The two-dimensional kf power spectrum is a two-dimensional representation of the MVDR power spectrum; Step 5: Perform noise reduction and enhancement processing on the kf two-dimensional power spectrum; This includes employing Gaussian smoothing to suppress local random noise, and normalization based on background noise floor to enhance effective turbulence features; Step 6: For the enhanced kf two-dimensional power spectrum, energy integration is performed using the tilt superposition method, and the optimal fluid velocity is calculated by combining coarse search and interpolation fine search.

2. The industrial pipeline flow velocity measurement method based on differential array and MVDR algorithm according to claim 1, characterized in that, The signal-to-noise ratio normalization based on background noise floor mentioned in step 5 refers to: estimating the background noise floor in the non-signal region for the frequency row of the kf two-dimensional power spectrum, and subtracting the floor value from the spectrum intensity to suppress background noise floor and retain effective turbulence features with high signal-to-noise ratio.

3. The industrial pipeline flow velocity measurement method based on differential array and MVDR algorithm according to claim 1, characterized in that, The Gaussian smoothing described in step 5 uses two-dimensional Gaussian filtering to eliminate local random noise in the kf two-dimensional power spectrum.

4. The industrial pipeline flow velocity measurement method based on differential array and MVDR algorithm according to any one of claims 1-3, characterized in that, The energy integration using the tilted superposition method described in step 6 refers to: for a given trial flow velocity u, along the dispersion relation line in the kf two-dimensional power spectrum. Perform energy integration: ； Where E(u) is the cumulative energy corresponding to the flow velocity u. The enhanced two-dimensional energy spectrum; The frequency variable corresponding to the short-time Fourier transform. and These are the lower and upper limits for frequency integration, respectively; The optimal flow velocity estimate is determined by finding the maximum value of the velocity energy function E(u).

5. The industrial pipeline flow velocity measurement method based on differential array and MVDR algorithm according to claim 4, characterized in that, By employing a two-stage strategy of "coarse search peak + spline interpolation fine search", high-precision and high-resolution solution of fluid velocity is achieved, ultimately obtaining the optimal velocity.

Images