Harmonic mode detection and removal method based on blind source separation

Abstract

A kind of harmonic modal detection and removal method based on blind source separation, rotating machinery generates harmonic excitation to structure in running state, based on the vibration response acceleration signal obtained by measurement, harmonic modal will appear in modal parameter identification, cause modal parameter identification error, low identification accuracy and other problems, harmonic modal should be detected and removed.In this method, based on the high similarity between modal expansion and blind separation principle, the vibration signal is analyzed under the condition of underdetermination sparse component, the energy peak scatter is clustered to estimate the mixing matrix by using the maximum energy method, and then the source signal is recovered based on the principle of minimum L1 norm;The probability density curve of each source signal is calculated, and the position of the harmonic modal is preliminarily determined;For the source signal containing harmonic, a spectral kurtosis curve diagram is drawn, so as to detect the frequency value corresponding to the harmonic modal;The corresponding harmonic source is removed, and the signal without harmonic interference is reconstructed.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of vibration signal analysis of mechanical structure, in particular to a vibration response signal processing method of rotating mechanical structure under operation state. BACKGROUND

[0002] Rotating mechanical structure is usually excited by periodic force during operation. The structure vibration of gas turbine under operation is mainly caused by multi-harmonic excitation of unbalance of rotor and critical speed, which can induce harmonic modal interference in the structure modal. The structure vibration during operation is the combination of response to random disturbance and harmonic excitation caused by rotating parts, and the complex periodic aerodynamic effect during operation of gas turbine can also introduce harmonic interference. These harmonics usually appear as false resonance peaks in the response spectrum to disturb the identification of modal parameters of real structure. Therefore, harmonic modal detection and removal of mechanical structure is the prerequisite for correct modal parameter identification under operation state.

[0003] The traditional harmonic detection method mainly includes time domain and frequency domain methods. The time domain method needs to perform single frequency filtering on the frequency component to be detected, which is low in detection and calculation efficiency and needs subjective judgment of human. Although the frequency domain method has the characteristics of high calculation efficiency and synchronous detection of full-band components without using digital filtering method to extract frequency components, when processing response signals containing harmonics, the frequency spectrum curve of harmonic detection will appear near-frequency interference, which is not easy to identify, and in the process of harmonic removal, the modal information of the original signal may be lost.

[0004] Therefore, a harmonic modal detection and removal method based on blind source separation is studied. The method does not need prior information of harmonic frequency, can automatically separate harmonic modal and system modal signals, and the blind source separation method used can still effectively separate signals under underdetermined condition, i.e. the number of measuring points is less than the number of modes. The harmonic detection is performed on each separated source signal, which avoids the process of subjective selection of frequency components to be detected and filtering, and has high detection efficiency. The process of blind separation actually decomposes complex signals into simple signals, which is more clear and accurate in frequency domain detection, avoiding the phenomenon of near-frequency interference. In the process of reconstruction of source signals, the complete information of system modal is retained, which is beneficial to subsequent modal parameter identification.

[0005] The above information disclosed in the background section is only used to enhance the understanding of the background of the present application, and therefore can contain information that is not prior art known to those of ordinary skill in the art. SUMMARY

[0006] In view of the problems in the prior art, the application provides a harmonic modal detection and removal method based on blind source separation, which is based on the high similarity between modal expansion and blind source separation principles to overcome the problems of low detection efficiency and the need for subjective human judgment in the harmonic modal detection and removal.

[0007] The application aims to realize the technical scheme, and the harmonic modal detection and removal method based on blind source separation comprises the following steps.

[0008] Step 1, collecting vibration response acceleration signals of each measuring point under the running state of a rotating machine;

[0009] Step 2, performing polynomial trend item elimination and smoothing processing on the vibration response acceleration signals to obtain a signal sequence;

[0010] Step 3, obtaining a time-frequency domain via short-time Fourier transform on the signal sequence;

[0011] Step 4, clustering each time-frequency scatter point of the time-frequency domain based on an energy maximum method to obtain an estimated value of a mixing matrix;

[0012] Step 5, converting a source recovery problem into an optimization problem based on an L1 norm minimum method to realize source signal recovery;

[0013] Step 6, performing probability density curve calculation on the source signals to find source signals with harmonic characteristics;

[0014] Step 7, performing spectral kurtosis detection on the source signals with harmonic characteristics to obtain harmonic modal frequencies;

[0015] Step 8, directly removing the source signals with obvious harmonic characteristics in the probability density curve, and filtering and removing the source signals with non-obvious harmonic characteristics based on the harmonic modal frequencies, wherein the harmonic characteristics are obvious when there are two peaks in the probability density curve, and the harmonic characteristics are non-obvious when only one peak can be seen;

[0016] Step 9, combining other source signals containing system characteristics and the source signals after removing the harmonic frequencies together, combining the mixing matrix, and realizing signal reconstruction, wherein the reconstructed signal is the signal after removing the harmonics;

[0017] Step 10, comparing the vibration response acceleration signals with the signal after removing the harmonics, and the difference before and after removing the harmonics can be clearly and intuitively seen, and the signal after removing the harmonics is the output result.

[0018] The harmonic modal detection and removal method based on blind source separation, in step 1, the number of the arrangement of the measuring point of the rotating machinery is less than the sum of the number of the harmonic modes and the number of the system modes, wherein the number of the harmonic modes is obtained by harmonic detection, and the number of the system modes is obtained by modal analysis.

[0019] The harmonic modal detection and removal method based on blind source separation, in step 2, in the polynomial trend item elimination, the sampled data of the measured vibration response acceleration signal is x k (k=1, 2, 3, …, n), n represents the number of sampling points, and the polynomial function is set as:

[0020]

[0021] In the formula, a j (j=0, 1, …, m) is each undetermined coefficient of the polynomial, m is the highest order of the polynomial, m=1-6 is taken, each undetermined coefficient of the function is determined, so that the error sum of squares of the function and x k is minimum, that is:

[0022]

[0023] E represents the error sum of squares of two functions, in order to minimize E, there is an extreme condition as follows:

[0024]

[0025] The partial derivative of E with respect to a is taken in turn, and an m+1 element linear equation group is generated:

[0026]

[0027] The equation group is solved to obtain the undetermined coefficient;

[0028] The formula for eliminating the polynomial trend item is:

[0029]

[0030] The five-point multiple smoothing processing includes:

[0031]

[0032]

[0033]

[0034]

[0035] In the above formula, x iThe measured vibration response acceleration signal sequence is (i=1, …, n), n represents the number of sampling points, y i The smoothed signal sequence is (i=1, …, n).

[0036] In the short-time Fourier transform, the formula of the frequency energy function is defined as follows:

[0037]

[0038] In the formula, E(f) is the frequency energy value, X(t, f) is the time-frequency transformed signal, t represents time, f represents frequency, and ∞ represents infinity; the energy distribution of a single channel in the frequency domain is calculated first in the time-frequency domain, then the energy of multiple channels at the same frequency point is added, that is:

[0039]

[0040] In the formula, E s (f) is the energy sum of signals received by all sensors in the frequency domain; Real(X i (t, f)) is the real part of the time-frequency transformed signal of the i-th sensor; Img(X i (t, f)) is the imaginary part of the time-frequency transformed signal of the i-th sensor; m is the number of measurement channels, t represents time, and f represents frequency; the peak value detection method is used to extract the peak frequency points.

[0041] In the step 5 of the method, in the L1 norm minimum method, all points in the sparse time-frequency domain are:

[0042]

[0043] In the formula, s i (t, f) is each source signal; A is a mixing matrix, x(t, f) is a measurement channel signal; s.t. is a constraint condition; t represents time, and f represents frequency; the L1 norm of the solution is required to be minimum, wherein s(t, f) is the optimal estimation of the source signal.

[0044] In the step 6 of the method, it is assumed that the mixed signal x(n) is composed of a Gaussian white noise signal and a harmonic signal, that is:

[0045] x(n) = g(n) + h(n),

[0046] Wherein, n=1, 2, …, N, N is the number of sampling points, g(n) is a Gaussian signal sequence and h(n) is a harmonic signal sequence, the sparse component analysis method is used to separate g(n) and h(n) signals, and the ideal Gaussian white noise signal meets the normal distribution, and the probability density is:

[0047]

[0048] The probability density of the two random variables X and T is f X (x) and f T (t), and x=h(t), according to the joint probability density function theory, the following can be obtained:

[0049] f X (x)=f t [h(x)]|h′(x)|, α<x<β,

[0050] Where h'(x) is the first derivative of h(x), and α and β are the value range of x, for the harmonic signal, it is considered that:

[0051] x=h(t)=a sin(2πft),

[0052] Where the time t belongs to the uniform distribution, and the probability density function is a constant f t , then the probability density function of the harmonic response is:

[0053]

[0054] For any frequency f, when the harmonic has a constant amplitude a, the above formula is always true, when the mean value of the harmonic component is 0, that is, x→a or x→-a, the probability density function tends to infinity, and there are two obvious peaks in the probability density curve.

[0055] In the harmonic modal detection and removal method based on blind source separation, the spectral kurtosis in step 7 can be defined as:

[0056]

[0057] In the formula: E{|X(f)| 4} is the fourth-order cumulant; E[|X(f)| 2 ] is the second-order cumulant,

[0058] For the mixed signal of Gaussian white noise signal and harmonic signal, the Fourier transform result is:

[0059] X(f)=G(f)+H(f),

[0060] The spectral kurtosis value of the mixed signal is:

[0061]

[0062] In the formula, f0 is a harmonic frequency value, and the spectrum kurtosis value at the corresponding frequency is -1, so that the harmonic component is determined.

[0063] The rotating machinery includes a gas turbine or a wind turbine.

[0064] In the method, multiple accelerator sensors are used to collect vibration response acceleration signals of each measuring point under the operating state of the rotating machinery.

[0065] Compared with the prior art, the method has the following advantages: the method combines sparse component analysis, probability density and spectrum kurtosis to realize detection and removal of harmonic modes, thereby improving operating modal parameter identification precision. BRIEF DESCRIPTION OF DRAWINGS

[0066] Various other advantages and benefits of the present application will become apparent to those of ordinary skill in the art, upon reading the following detailed description of the preferred embodiments. The accompanying drawings are included to provide a description of preferred embodiments of the application and are not intended to limit the scope of the application. It should be apparent to those of ordinary skill in the art that the drawings described below are merely illustrative of embodiments of the application and should not be construed to limit the scope of the application. It should be apparent to those of ordinary skill in the art that other embodiments can be employed without departing from the scope of the application. Throughout the drawings, like reference numerals will be used to refer to like components.

[0067] In the drawings:

[0068] Figure 1 is a flowchart of a method for detecting and removing harmonic modes based on blind source separation according to an embodiment of the present application;

[0069] Figure 2 is a schematic diagram of a three-degree-of-freedom mass-spring-damper system according to an embodiment of the present application;

[0070] Figure 3 is a three-dimensional time-domain scatter plot of a collected vibration response signal according to an embodiment of the present application;

[0071] Figure 4 is a three-dimensional time-frequency-domain scatter plot after sparse transformation according to an embodiment of the present application;

[0072] Figure 5is a frequency energy three-dimensional scatter diagram provided by one embodiment of the present application;

[0073] Fig. 6(a), Fig. 6(b) are probability density curve harmonic detection diagrams provided by one embodiment of the present application;

[0074] Figure 7 is a spectral kurtosis curve harmonic detection diagram of source signal 1 provided by one embodiment of the present application;

[0075] Figure 8 is a spectral kurtosis curve harmonic detection diagram of source signal 3 provided by one embodiment of the present application;

[0076] Figure 9 is an energy peak diagram containing harmonic modal provided by one embodiment of the present application;

[0077] Figure 10 is an energy peak diagram after removing harmonic modal provided by one embodiment of the present application.

[0078] The present application will be further explained in conjunction with the accompanying drawings and embodiments. DETAILED DESCRIPTION

[0079] The present application will be further explained in conjunction with the accompanying drawings and embodiments. Figures 1 to 10 The specific embodiments of the present application will be described in further detail below. Although specific embodiments of the present application are shown in the drawings, it should be understood that the present application can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided so that this application will be thoroughly and completely understood, and will fully convey the scope of the application to those skilled in the art.

[0080] It should be noted that certain terms are used throughout the present specification and claims which have particular meanings as set forth below. Those skilled in the art will understand that not all embodiments of the present application will include all of the terms as described below. Additionally, some embodiments of the present application will include terms that are not listed below. Accordingly, the use of particular terms in particular embodiments of the present application should not be used to restrict the scope of the present application or that particular embodiment.

[0081] For the purpose of promoting an understanding of the principles of the application, reference will now be made to the embodiments illustrated in the drawings and specific language will be used to describe the same. It will, nevertheless, be understood that no limitation of the scope of the application is intended by this specific disclosure. Alterations, modifications, and improvements to such specific embodiments are contemplated and, where applicable, some substitutions apart from those specifically mentioned are possible without departing from the spirit and scope of the application.

[0082] As shown in the figure, the harmonic modal detection and removal method based on blind source separation includes: Figures 1 to 10 As shown in the figure, the harmonic modal detection and removal method based on blind source separation includes:

[0083] Step 1, collect the vibration response acceleration signal of each measuring point under the running state of the rotating machinery;

[0084] Step 2, polynomial trend term elimination and smoothing processing are performed on the vibration signal;

[0085] Step 3, time-frequency transformation is performed on the time-domain vibration signal to convert to a time-frequency domain with better sparsity;

[0086] Step 4, based on the energy maximum method, each time-frequency scatter point is clustered to obtain the estimated value of the mixing matrix;

[0087] Step 5, based on the L1 norm minimum method, the source recovery problem is converted into an optimization problem to realize source signal recovery;

[0088] Step 6, the probability density curve of each source signal is calculated to find the source signal with harmonic characteristics;

[0089] Step 7, the spectral kurtosis of the source signal with harmonic characteristics is detected to obtain the harmonic modal frequency;

[0090] Step 8, for the source signal with obvious harmonic characteristics in the probability density curve, direct removal is performed, and for the source signal without obvious harmonic characteristics, filtering removal should be performed according to the harmonic modal frequency;

[0091] Step 9, other source signals containing system characteristics and the source signals after removing the harmonic frequency are combined together to realize signal reconstruction in combination with the mixing matrix, and the reconstructed signal is the signal after removing the harmonic;

[0092] Step 10, the original signal is compared with the signal after removing the harmonic.

[0093] In the harmonic modal detection and removal method based on blind source separation, the underdetermined blind source separation situation can be processed, the number of measuring points is usually less than the sum of the number of harmonic modes and the number of system modes, and the method is suitable for the case of few measuring points.

[0094] In the harmonic modal detection and removal method based on blind source separation, in step 2, in the polynomial trend term elimination, the sampled data of the measured vibration response acceleration signal is {x k}(k=1,2,3,…,n), n represents the number of sampling points, and the polynomial function is set as:

[0095]

[0096] In the formula, a j (j=0,1,…,m) is the undetermined coefficient of the polynomial, m is the highest order of the polynomial, m is taken as 1-6, and the function each undetermined coefficient, so that the function and the error square sum of x k is minimum, that is:

[0097]

[0098] E represents the error square sum of two functions, in order to minimize E, the extreme value condition should be met:

[0099]

[0100] Take the partial derivative of E to a in turn, and generate an m+1 linear equation group:

[0101]

[0102] Solve the equation group to obtain the undetermined coefficient;

[0103] The formula for eliminating the polynomial trend term is:

[0104]

[0105] Five-point multiple smoothing processing includes:

[0106]

[0107]

[0108]

[0109]

[0110] In the above formula, x i (i=1,..., n) is the original vibration response acceleration signal sequence, n represents the number of sampling points, y i (i=1,..., n) is the signal sequence after smoothing processing.

[0111] In the harmonic modal detection and removal method based on blind source separation, in step 3, selecting a suitable sparse transform method is a prerequisite for performing sparse component analysis method, and the short-time Fourier transform in time-frequency transform is the simplest and most effective method to realize sparsification, which can be written as the following formula:

[0112]

[0113] In the formula, w(τ-t) represents a window function moving with time t, x(τ) is a truncated signal, and X(t, f) is a signal after time-frequency transform, and the short-time Fourier transform considering the phase difference is:

[0114]

[0115] According to the Parseval theorem, this can be written as:

[0116]

[0117] where, is the Fourier transform of the window function, which has a frequency support in the range [-Δ, Δ], is the Fourier transform of the initial signal. For two harmonic signals of different frequencies:

[0118]

[0119]

[0120] where f1and f2represent the respective harmonic frequencies and A1and A2represent the respective amplitudes. The corresponding short-time Fourier transform can be written as:

[0121]

[0122]

[0123] Two source signals are linearly instantaneous mixed and measured by two measurement channels:

[0124] y2(t) = a 11 · x1(t) + a 12 · x2(t)

[0125] y2(t) = a 21 · x1(t) + a 22 · x2(t)

[0126] where y1(t), y2(t) represent the mixed measurement channel signals, a 11 , a 12 , a 21 , a 22 represent the mixing coefficients. Based on the linearity of the short-time Fourier transform, the short-time Fourier transform of the two measurement channels can be written as:

[0127] Y1(t, f) = a 11 · X1(t, f) + a 12 · X2(t, f)

[0128] Y2(t, f) = a 21 · X1(t, f) + a 22 · X2(t, f)

[0129] It should be noted that the interval of harmonic frequency should be greater than the frequency maintaining range of window function, that is, |f1-f2|>2Δ, so the window function with longer window length is selected to ensure that the frequency maintaining range is small enough. The two harmonic signals do not overlap in the time-frequency domain, and for f∈[f1-Δ, f1+Δ], the following formula can be constructed:

[0130]

[0131] For f∈[f2-Δ, f2+Δ], the following formula can be constructed:

[0132]

[0133] In the formula, Re represents the real part of a complex number. According to the above formula, when drawing the real part scatter plot of the short-time Fourier transform of the measurement channel signal in this case, for the frequency point, it is actually equivalent to drawing a series of special case Lissajous figures, which in this example are straight lines, and the direction of the straight line is the direction of the column vector of the mixing matrix, as shown in Or These cluster straight lines appear in the scatter plot, which is the theoretical basis for estimating the mixing matrix.

[0134] In the harmonic modal detection and removal method based on blind source separation, the mixing matrix is estimated based on the energy maximum method in step 4. For estimating the mixing matrix of the source signal, there are mainly potential function method and clustering method. The potential function method is to expand the angle of the cluster straight line to the polar coordinate axis, and determine the angle between the cluster straight line and the coordinate axis by calculating the peak point, so as to estimate the column vectors of the mixing matrix, but it can only be used for two-channel mixed signals, and has certain limitations. While estimating the center of the cluster straight line to determine the mixing matrix, it is not limited by the number of channels, so it is widely used. Whether at the harmonic frequency or at the structural resonance frequency, the energy amplitude corresponding to the frequency is very large, so it is reasonable to separate the harmonic modal and the system modal by clustering through the energy maximum method. From the perspective of signal energy, the signal energy in the frequency band is estimated. Based on the short-time Fourier transform of the signal, the formula of the frequency energy function is defined as:

[0135]

[0136] In the formula, E(f) is the frequency energy value; X(t, f) is the time-frequency transformed signal, t represents time, f represents frequency, and ∞ represents infinity; the energy distribution of a single channel in the frequency domain is calculated in the time-frequency domain first, and then the energy of multiple channels at the same frequency point is added, that is:

[0137]

[0138] In the formula, E s(f) the energy sum of all received signals in the frequency domain; Real(X i (t, f)) is the real part of the time-frequency transform of the i-th sensor signal; Img(X i (t, f)) is the imaginary part of the time-frequency transform of the i-th sensor signal; m is the number of measurement channels, t represents time, and f represents frequency. The peak value detection method is used to extract the peak frequency points. In the peak value detection method, if the commonly used zero derivative method is used, the first derivative will accidentally pass through zero due to the noise in the measured signal, resulting in false detection. If some low-pass filter is used to smooth the curve, the effective information of the signal will be lost to some extent. Therefore, the method of finding the peak value and the valley value is used to define a peak value with lower points around it. A peak value threshold e is set, which requires that the peak value and the points around it have at least a threshold level difference, so as to define it as a peak value.

[0139] In the method, in step 5, the L1 norm minimum method is used to recover the source signal. Since the number of unknown sources is greater than the number of equations, the equation does not have a unique solution. A constraint condition is introduced, and the L1 norm of the solution is minimized. Then, the source recovery is converted into an optimization problem. For all points in the sparse domain:

[0140]

[0141] In the formula, s i (t, f) is each source signal; A is a mixing matrix, x(t, f) is a measurement channel signal, t represents time, f represents frequency, and s.t. is a constraint condition; the formula converts the source recovery into an optimization problem. Since the number of unknown sources is greater than the number of equations, the equation does not have a unique solution. A constraint condition is introduced, and the L1 norm of the solution is minimized. The optimal estimate of the source signal is s(t, f). In theory, the equation uses the L0 norm, which represents the number of non-zero elements in a vector. Naturally, the sparse solution with the fewest non-zero elements can be found among all feasible solutions. However, the constraint equation is non-convex and highly discrete, and it is difficult to solve by numerical calculation. However, the sparse representation theory points out that if the solution is indeed sparse, the L1 norm, as the optimal convex approximation of the L0 norm, can obtain the same sparse solution as the L0 norm and is easier to optimize and solve.

[0142] Taking the two measurement channel signals as an example, the constraints show that x(t) must be linearly decomposed in the direction of any two column vectors in the mixing matrix to find the shortest path from O to the observed signal x(t). a1a2a3 are the three column vectors of the mixing matrix, and point E is any point of the observed signal. At this time, the mixing matrix has no inverse matrix, and the solution with the smallest L1 norm is the solution that satisfies the shortest distance from point E to the origin. Reduce the dimension of the mixing matrix to generate three 2×2 sub-matrices (that is, three combinations of the column vectors of the mixing matrix). Assuming that point E is between vectors a2a3, the sub-matrix composed of a1a2 is the optimal solution matrix, and OBOC is the optimal solution of point Z in the direction of a2a3. The inverse of the sub-matrix is ​​multiplied by the observation point to obtain an estimate of the source signal. The specific steps of the L1 norm minimization algorithm are as follows

[0143] (1) Find the mixing matrix There are m×m dimensional submatrices, let B k ,

[0144] (2) Find all possible solutions for a point X in the sparse domain,

[0145] (3) Find the L1 norm of each solution and take the one with the smallest norm As the optimal estimate of the source signal,

[0146] (4) Repeat the above steps to find the optimal solution for all points in the sparse domain.

[0147] (5) Perform sparse inverse transform on the source signal to obtain the time domain estimation of the source signal.

[0148] In the harmonic mode detection and removal method based on blind source separation, in step 6, a probability density function is calculated based on the probability density curve of each source signal. The probability density function represents the probability that the signal amplitude falls within a specified area. The detection method is based on the estimation of the probability density curves of the harmonic signal and the Gaussian white noise signal. The core idea is to estimate the probability density curve of the separated single-frequency time domain signal and distinguish the harmonic response from the random response based on the different statistical characteristics.

[0149] Assume that the mixed signal x(n) consists of Gaussian white noise signal and harmonic signal, that is:

[0150] x(n)=g(n)+h(n)

[0151] Wherein, n=1, 2, …, N, N is the number of sampling points; g(n) is a Gaussian signal sequence and h(n) is a harmonic signal sequence, the aforementioned sparse component analysis method can be used to separate g(n) and h(n) signals for research respectively. The ideal Gaussian white noise signal conforms to normal distribution and its probability density is:

[0152]

[0153] Although the structure in actual engineering is excited by multiple Gaussian white signals with different mean values and different variances and is independent of each other, and the probability density curve of the signal is not an ideal Gaussian white noise signal, the probability density curve is still approximately subject to normal distribution, and even if the signal is subjected to a series of signal processing, the distribution form of the probability density function will not be changed.

[0154] Suppose that the probability densities of two random variables X and T are f x (x) and f T (t) respectively, and x=h(t), according to the joint probability density function theory, the following equation can be obtained:

[0155] f X (x)=f t [h(x)]|h′(x)|, α<x<β

[0156] Where h'(x) is the first derivative of h(x), and α and β are the value ranges of x. For a harmonic signal, the following equation is considered:

[0157] x=h(t)=asin(2πft)

[0158] Where the time t belongs to uniform distribution, and its probability density function is a constant f t , then the probability density function of the harmonic response is:

[0159]

[0160] For any frequency f, when the harmonic has a constant amplitude a, the above equation is always true. When the mean value of the harmonic component is 0, i.e. x→a or x→-a, the probability density function tends to infinity, and there are two obvious peaks in the probability density curve.

[0161] In the harmonic modal detection and removal method based on blind source separation, the spectrum kurtosis curve of each source signal is calculated in step 7, and the harmonic frequency is detected.

[0162] The spectrum kurtosis can be defined as:

[0163]

[0164] In the formula: E{|X(f)| 4- Fourth order cumulant; E[|X(f) 2 - Second order cumulant.

[0165] The time series is divided into M segments, and Fourier transform X i (f) (i = 1, 2, …, M) is substituted into the above formula to calculate the unbiased estimation formula of the fourth order cumulant:

[0166]

[0167] Similarly, the unbiased estimation formula of the second order cumulant is obtained:

[0168]

[0169] Finally, the unbiased estimation formula of the spectral kurtosis is obtained:

[0170]

[0171] For a linear system, the response will not change the signal characteristics of the input. For a Gaussian white noise signal g(n), after Fourier transform, it is G(f). The spectral kurtosis value SK g (f) = 0 is calculated by the above formula.

[0172] Let the harmonic signal with frequency f0 and phase be:

[0173]

[0174] After Fourier transform, it is:

[0175]

[0176] where δ(f0) is the Dirac function. The spectral kurtosis value SK h (f0) = -1 is calculated by the formula.

[0177] The linear property of Fourier transform, for a mixed signal of Gaussian white noise and harmonic signal, its Fourier transform result is:

[0178] X(f) = G(f) + H(f)

[0179] Therefore, the spectral kurtosis value of the mixed signal is:

[0180]

[0181] Therefore, according to the spectral kurtosis value at the corresponding frequency being -1, it can be determined that it is a harmonic component.

[0182] The one kind based on blind source separation's harmonic mode detection and removal method, in step 9, the harmonic signal of detection is removed, if probability density indicates that entire source signal is harmonic signal, then directly to entire source signal removal.For example, the i-th source signal is harmonic signal S i (t,f), S i (t,f) in source signal is removed to obtain S R (t,f), and the corresponding a i in mixing matrix is removed to obtain A R , and the reconstructed signal X R (t,f) = A R S R (t,f) is obtained, and then the signal is converted to time domain to obtain x R (t).If harmonic energy is too small in entire source signal, probability density fails to clearly indicate that it is harmonic, and spectral kurtosis method indicates that the harmonic frequency in source signal, then entire source signal does not need to be removed, and stop-band filtering is applied to remove corresponding harmonic frequency component, and then participate in subsequent signal reconstruction.

[0183] In step 10, the initial signal is compared with the signal after removing the harmonic, and more intuitive results are obtained.

[0184] In one embodiment, as Figures 1 to 10 shown, a harmonic mode detection and removal method based on blind source separation includes,

[0185] 1) Measure the vibration response acceleration signal of all measuring points. For more clear and concise description of the process and results of the method, a three-degree-of-freedom mass-spring-damper simulation system is established, as shown in the accompanying Figure 2 , white noise excitation is applied to it, and the vibration response acceleration signal is measured, and the method is described based on the simulation data. The system parameters are M1=M2=M3=1 kg, K1=K2=1×10 6 N / m, K3=1.5×10 6 N / m, C=M+10 -6 K, then the modal parameters theoretical value can be obtained by modal analysis theory, as shown in Table 1.

[0186] Table 1: Theoretical values of natural frequency and damping ratio of three-degree-of-freedom system

[0187]

[0188] 2) select a certain period of signal, the signal is eliminated and five-point polynomial trend term multiple smoothing processing; the purpose of eliminating the polynomial trend term is that the signal often deviates from the baseline due to the zero drift caused by temperature change, the unstable low frequency performance outside the frequency range of the sensor and the interference of the environment around the sensor, and even the size of the baseline deviation will change with time. The whole process is called the trend term of the signal, which directly affects the correctness of the signal and should be removed;

[0189] The specific steps of polynomial trend elimination are:

[0190] The measured vibration signal sampling data is {x k}(k = 1, 2, 3, …, n), set a polynomial function as:

[0191]

[0192] Generally, m = 1 ~ 6 can be taken according to the situation;

[0193] Determine the function Each undetermined coefficient, so that the function The error sum of squares of x k is minimum, that is:

[0194]

[0195] The extreme condition of E is satisfied:

[0196]

[0197] Take the partial derivative of E to a in turn, which can produce an m+1 linear equation group:

[0198]

[0199] Solve the equation group to obtain the undetermined coefficient;

[0200] The formula for eliminating the polynomial trend term is:

[0201]

[0202] The purpose of smoothing processing is that the collected vibration signal often superimposes noise, in addition to some periodic interference signals, there are also some irregular random interference signals. Because the frequency band of the interference signal is wide, sometimes the proportion of high frequency component is very large, so that the collected signal presents many burrs and is very rough. In order to weaken the influence of interference signal and improve the smoothness of vibration curve, the signal needs to be smoothed;

[0203] The specific steps of five-point multiple smoothing processing are:

[0204]

[0205]

[0206]

[0207]

[0208] x (i) = x (i) - x (i) (i = 1, 2,..., n) (1) where x (i) is the original vibration response acceleration signal sequence, n represents the number of sampling points, x (i) is the smoothed signal sequence, and x (i) is the mean value of the signal sequence. i (i = 1,..., n) is the original vibration response acceleration signal sequence, n represents the number of sampling points, y i (i = 1,..., n) is the smoothed signal sequence.

[0209] 3) Convert the time-domain signal to the time-frequency domain by short-time Fourier transform. The short-time Fourier transform in time-frequency transform is the simplest and most effective method to achieve sparsity, which can be written as follows:

[0210]

[0211] In the formula, w (τ-t) represents a window function moving with time t, x (τ) is a truncated signal, and X (t, ff) is a signal after time-frequency transform. The response signals of the three measurement channels in the example are selected, and a three-dimensional scatter plot in the time domain is drawn. It can be found that the scatter signals are mixed together in the time domain, and the sparsity is poor, as shown in Figure 3 . The time-domain signal is converted to the time-frequency domain by using the above-mentioned sparse transform method based on short-time Fourier transform, and a three-dimensional scatter plot in the time-frequency domain is drawn. It can be found that the three-dimensional plot presents obvious clustering straight lines, which meets the prerequisite assumption of good sparsity, as shown in Figure 4 .

[0212] 4) Estimate the mixing matrix based on the energy maximum method. First, define the formula of the frequency energy function as follows:

[0213]

[0214] In the formula: E (f) - frequency energy value; X (t, f) - signal after time-frequency transform; t represents time, f represents frequency, and ∞ represents infinity. For the modal response signal, the energy is concentrated at certain frequency points, and the clustering direction at the local maximum frequency point represents the clustering direction of the source signal, and the mixing matrix can be estimated by only calculating the clustering direction at these frequency points. The maximum value of each clustering line is located at the end of the line, and the cosine distance of other points to the end is the same. The specific process is as follows: first, calculate the energy distribution of a single channel in the frequency domain in the time-frequency domain, and then add the energies of multiple channels at the same frequency point, that is:

[0215]

[0216] In the formula: Es (f) - Energy sum of all sensor received signals in frequency domain; Real(X i (t, f)) - Real part of the i-th sensor signal after time-frequency transform; Img(X i (t, f)) - Imaginary part of the i-th sensor signal after time-frequency transform; t represents time, f represents frequency, and m represents the number of measurement channels. The peak value detection method is used to extract each peak frequency point. If the commonly used zero derivative method is used in the peak value detection method, the first derivative will accidentally pass zero due to the noise in the measured signal, resulting in false detection. If some low-pass filter is used to smooth the curve, the effective information of the signal will be lost to some extent. Therefore, the method of finding the peak value and the valley value is adopted here, that is, a peak value is defined as having lower points around it. A peak value threshold is set, which requires that the peak value and the surrounding points have at least a difference of the threshold level, so as to define it as a peak value. Each point in the clustering line corresponds to the frequency band of the source in the time-frequency domain representation, as shown in Figure 5 .

[0217] 5) Source signal recovery based on L1 norm minimization. Since the number of unknown sources is greater than the number of equations, the equation does not have a unique solution. A constraint condition is introduced, which requires the L1 norm of the solution to be minimized. Then the source recovery is converted into an optimization problem. For all points in the sparse time-frequency domain:

[0218]

[0219] In the formula: s i (t, f) is each source signal; A is the mixing matrix, x(t, f) is the measurement channel signal, t represents time, f represents frequency, and s.t. represents the constraint condition; the formula converts the source recovery into an optimization problem. Since the number of unknown sources is greater than the number of equations, the equation does not have a unique solution. A constraint condition is introduced, which requires the L1 norm of the solution to be minimized, where s(t, f) is the optimal estimate of the source signal. In theory, this equation uses L0 norm, which represents the number of non-zero elements in a vector. Naturally, the sparse solution with the least non-zero elements can be found among all feasible solutions. However, this constraint equation is non-convex and highly discrete, making it difficult to solve numerically. However, the sparse representation theory points out that if the solution is indeed sparse, the L1 norm, as the optimal convex approximation of the L0 norm, can obtain the same sparse solution as the L0 norm and is easier to optimize and solve. In this embodiment, a total of 5 source signals are separated by this method.

[0220] 6) Calculate the probability density curve of the 5 source signals to make preliminary harmonic detection. As shown in FIG. 6(a) and FIG. 6(b), the horizontal axis is amplitude, and the vertical axis is the probability density function value, which represents the probability of the signal amplitude falling in the specified area. This detection method is based on the estimation of the probability density curve of the harmonic signal and the Gaussian white noise signal. The core idea is to estimate the probability density curve of the separated single-frequency time-domain signal, and distinguish the harmonic response and the random response according to the different statistical characteristics. The harmonic signal will show two obvious peaks, while the white noise signal will have only one peak. In this embodiment, it can be seen that the probability density curves of source signal 1 and source signal 3 show two obvious peaks, so they are determined as harmonic signals.

[0221] 7) Calculate the spectral kurtosis curve of source signal 1 and source signal 3 and locate the frequency corresponding to the harmonic. The spectral kurtosis can be defined as:

[0222]

[0223] wherein: E{|X(f)| 4}-fourth-order cumulant; E[|X(f) 2 ]-second-order cumulant.

[0224] Divide the time series into M segments, and perform Fourier transform X i (f) on each segment (i=1, 2, …, M), substitute the above formula to calculate the unbiased estimation formula of the fourth-order cumulant:

[0225]

[0226] Similarly, the unbiased estimation formula of the second-order cumulant is obtained:

[0227]

[0228] Finally, the unbiased estimation formula of the spectral kurtosis is obtained:

[0229]

[0230] For a linear system, the response will not change the signal characteristics of the input. For a Gaussian white noise signal g(n), the Fourier transform is G(f), and the spectral kurtosis value SK g (f) calculated by the above formula is 0.

[0231] Let the harmonic signal with frequency f0 and phase be:

[0232]

[0233] After Fourier transform, it is:

[0234]

[0235] where δ(f0) is the Dirac function, the spectral kurtosis value SK at f0 is calculated by the formula h (f0) = -1.

[0236] The linear property of Fourier transform, for the mixed signal of Gaussian white noise signal and harmonic signal, the Fourier transform result is:

[0237] X(f) = G(f) + H(f)

[0238] Therefore, the spectral kurtosis value at the corresponding frequency is -1, which can be determined as a harmonic component. As shown in

[0239]

[0240] Therefore, the spectral kurtosis value at the corresponding frequency is -1, which can be determined as a harmonic component. As shown in Figure 7 and Figure 8 The spectral graph and the corresponding spectral kurtosis graph of the source signal 1 and the source signal 3 are drawn, and it can be found that the spectral kurtosis value of the source signal 1 at 50Hz is close to -1, and the spectral kurtosis value of the source signal 3 at 150Hz is close to -1, so the harmonic frequencies are 50Hz and 150Hz respectively.

[0241] 8) Harmonic modal removal and signal reconstruction. The detected harmonic signal is removed, and if the probability density indicates that the entire source signal is a harmonic signal, the entire source signal is directly removed. For example, the i-th source signal is a harmonic signal S i (t, f), S i (t, f) in the source signal is removed to obtain S R (t, f), and the corresponding a i in the mixing matrix is removed to obtain A R , and the reconstructed signal x R (t, f) = A R S R (t, f) is obtained, and then the signal is converted to time domain to obtain x R (t). If the harmonic energy is too small in the entire source signal, the probability density cannot clearly indicate that it is a harmonic, and the spectral kurtosis method indicates the harmonic frequency in the source signal, then the entire source signal does not need to be removed, and the corresponding harmonic frequency component is removed by applying a band-stop filter, and then participates in the subsequent signal reconstruction.

[0242] In this embodiment, the probability density can clearly indicate that the source signal 1 and the source signal 3 are harmonics, so the two source signals are directly removed, and the corresponding mixing coefficients in the mixing matrix are removed, and then the signal is reconstructed based on the reconstruction formula.

[0243] 9) The initial signal is compared with the deharmonized reconstructed signal. The frequency energy graphs of the initial signal and the deharmonized reconstructed signal are drawn respectively, and the corresponding peak values at the harmonic frequencies are observed. It can be found that the peak values corresponding to 50Hz and 150Hz in the initial signal are very obvious, which are harmonic modes, and the peak values at 50Hz and 150Hz in the deharmonized reconstructed signal are not obvious, and the remaining three frequency peak values are all peak values of system modes, and finally the detection and removal of the harmonic modes are realized.

[0244] Although the embodiments of the present application are described above with reference to the drawings, the present application is not limited to the above-described specific embodiments and application fields, and the above-described specific embodiments are merely illustrative and instructive, but not restrictive. A person of ordinary skill in the art can make many forms under the guidance of the present specification and without departing from the scope protected by the claims of the present application, and these all belong to the protection of the present application.

Claims

1. A method for harmonic modal detection and removal based on blind source separation, characterized in that, It comprises the following steps, Step 1, collecting vibration response acceleration signals of each measuring point of the rotating machinery in the running state; Step 2, performing polynomial trend term elimination and smoothing processing on the vibration response acceleration signals to obtain a signal sequence; Step 3, obtaining a time-frequency domain via short-time Fourier transform on the signal sequence; Step 4, clustering each time-frequency scatter point of the time-frequency domain based on the energy maximum method to obtain an estimated value of a mixing matrix; Step 5, converting a source recovery problem into an optimization problem based on the L1 norm minimum method to realize source signal recovery; Step 6, calculating a probability density curve of the source signal to find a source signal with harmonic characteristics; Step 7, performing spectral kurtosis detection on the source signal with harmonic characteristics to obtain a harmonic modal frequency; Step 8, directly removing the source signal with obvious harmonic characteristics in the probability density curve, and filtering and removing the source signal with non-obvious harmonic characteristics based on the harmonic modal frequency, wherein the harmonic characteristics are obvious when there are two peaks in the probability density curve, and the harmonic characteristics are non-obvious when there is only one peak; Step 9, combining other source signals containing system characteristics and the source signal after removing the harmonic frequency together to realize signal reconstruction in combination with the mixing matrix, wherein the reconstructed signal is the signal after removing the harmonic; Step 10, the signal after removing the harmonic is the output result.

2. The method of claim 1, wherein, In step 1, the number of measuring points of the rotating machinery is less than the sum of the number of harmonic modes and the number of system modes, wherein the number of harmonic modes is obtained through harmonic detection, and the number of system modes is obtained through modal analysis.

3. The method of claim 1, wherein, In the polynomial trend term elimination in Step 2, the measured vibration response acceleration signal sample data is , represents the sample point number, and the polynomial function is set as: ， In the formula, are the undetermined coefficients of the polynomial, is the highest order of the polynomial, and are determined, the undetermined coefficients of the function are determined, and the error square sum of the function and is minimum, specifically: ， denotes the sum of the squared errors of two functions, and for minimizing, there is an extreme condition: ， Take in order To Partial derivatives, resulting in a Linear equations: ， Solve the equation set to obtain the undetermined coefficients; The formula for eliminating the polynomial trend term is: ， Five-point multiple smoothing processing includes: ， ， ， ， In the above formula is a measured vibration response acceleration signal sequence, represents the number of sampling points, is a smoothed signal sequence.

4. The method of claim 1, wherein, In the short-time Fourier transform, first define the formula of the frequency energy function as: ， In the formula: is the frequency energy value; is the signal after time-frequency transformation; represents time, represents frequency, represents infinity; in the time-frequency domain, the energy distribution of a single channel in the frequency domain is calculated first, and then the energy of multiple channels at the same frequency point is added, specifically: ， wherein: is the energy in the frequency domain of the signal received by all sensors; is the real part of the time-frequency transform of the signal of the th sensor; is the imaginary part of the time-frequency transform of the signal of the th sensor; m is the number of measurement channels, representative time, representative frequency, the peak value frequency point is extracted by using the peak value detection method.

5. The method of claim 1, wherein, In the step 5, based on the L1 norm minimum method, all points in the sparse time-frequency domain are: ， wherein: are the individual source signals; A is a mixing matrix, is the measured channel signal; is the constraint; represents time, represents frequency; the L1 norm of the solution of the equation is minimized, where is the optimal estimate of the source signal.

6. The method of claim 1, wherein, In the step 6, assume that the mixed signal x(n) is composed of a Gaussian white noise signal and a harmonic signal, i.e., ， wherein, N is the number of sampling points; is a Gaussian signal sequence and is a harmonic signal sequence, the sparse component analysis method is used to separate and signals, and the ideal Gaussian white noise signal conforms to the normal distribution, and the probability density is: ， The probability densities of the two random variables X and T are and , and According to the joint probability density function theory, we have: ， wherein is the first derivative of and is the range of values of x, for a harmonic signal, considered to be: ， where time t belongs to a uniform distribution with a constant probability density function The probability density function of the harmonic response is then given by ， For any frequency f, the above holds constant when the harmonic has a constant amplitude a, and when the mean of the harmonic component is zero, i.e. or then the probability density function tends to infinity and there are two distinct peaks in the probability density curve.

7. The method of claim 6, wherein, In the step 7, the spectral kurtosis can be defined as: ， wherein: is the fourth order cumulant; is the second order cumulant, For the mixed signal of the Gaussian white noise signal and the harmonic signal, the Fourier transform result is: ， The spectral kurtosis value of the mixed signal is: ， In the formula, is a harmonic frequency value, and the spectrum kurtosis value at the corresponding frequency is -1, so it is determined that it is a harmonic component.

8. The method of claim 1, wherein, The rotating machinery includes a gas turbine or a wind turbine.

9. The method of claim 1, wherein, A plurality of accelerator sensors are used to collect vibration response acceleration signals of each measuring point of the rotating machinery in the running state. A plurality of accelerator sensors are used to collect vibration response acceleration signals of each measuring point of the rotating machinery in the running state.

Images