CN104807534B - Equipment eigentone self study recognition methods based on on-line vibration data - Google Patents

Abstract

The invention discloses a self-learning identification method for the natural vibration mode of equipment based on online vibration data, and uses a wavelet packet transformation algorithm to perform denoising on the equipment vibration signal x(n) once to obtain the denoised signal Using singular value decomposition to analyze the signal Perform secondary denoising to obtain the vibration signal after secondary denoising Denoised vibration signal using windowed discrete Fourier algorithm Perform spectrum analysis to calculate the vibration spectrum of the equipment; use the self-learning algorithm to train the vibration spectrum of the equipment, and finally obtain the natural characteristic frequency and amplitude of the equipment. The beneficial effect of the invention is that when the research object is a relatively complex system, the characteristic frequency of the components can be accurately identified.

Description

Self-learning identification method of equipment natural vibration mode based on online vibration data

technical field

The invention belongs to the technical field of on-line monitoring, and relates to a self-learning identification method of a natural vibration mode of equipment based on online vibration data.

Background technique

Online monitoring technology is the most basic measure to estimate the health status of equipment, especially for mechanical equipment (or systems), vibration monitoring and analysis is the most practical and effective method. The vibration signal of the equipment is closely related to its mechanical structure, and the condition of its internal mechanical structure during the operation of the equipment can be directly and effectively reflected from the vibration signal. By processing and analyzing the vibration signal of the equipment, the inherent vibration mode of the equipment under normal operating conditions is extracted, which is used to distinguish the vibration mode when the equipment is abnormal or faulty, so as to make timely decisions when the abnormal vibration model appears, and formulate corresponding maintenance plan. The invention proposes a method for automatic learning, identification and extraction of equipment vibration natural modes based on vibration data. After analyzing and processing the vibration data and training for a reasonable period of time, the natural vibration parameters of the equipment are calculated and used as the health of the normal operation of the equipment. Features, so as to provide a feasible evaluation index for the operation health of the equipment, and ensure the healthy operation and maintenance of the equipment.

The natural vibration mode (mode) of the equipment refers to the natural vibration characteristics of the mechanical structure, including the natural frequency, damping ratio and mode shape of the equipment. If the characteristics of the main modes of the mechanical structure of the equipment are known, the actual vibration response of the structure under various vibration sources can be predicted. Therefore, modal analysis is to analyze these modes to obtain the corresponding modal parameters, which is an important method for structural dynamic design and equipment fault diagnosis. Although the vibration behavior of actual equipment (such as wind turbines, transformers, ships, aircraft, etc.) varies and is complex, modal analysis provides an effective way to study the actual structural characteristics of various equipment.

The existing equipment natural vibration mode identification method is mainly the computational modal analysis method. First, use the 3D modeling software to construct a high-precision 3D solid model of the equipment; then, import the 3D model of the equipment into the finite element analysis software, According to the actual situation, set the material property parameters of the equipment, divide the three-dimensional model into finite element networks, add displacement constraints, and obtain the finite element model of the equipment; finally, perform modal analysis on the finite element model to obtain the natural frequency and mode of the model state.

In the prior art, in the finite element modal analysis, especially when the research object is a relatively complex system, the calculated vibration modes are often too messy, and it is difficult to identify the eigenfrequencies of the parts that are really concerned; or when many local When the mode shape constitutes a certain order of vibration, it will be very difficult to judge the value of the eigenfrequency.

The present invention does not need to establish complicated models, directly collects the vibration data of different parts of the equipment, obtains the modal parameters through analysis and processing, and identifies the natural vibration mode of the equipment.

Contents of the invention

The purpose of the present invention is to provide a self-learning identification method for the natural vibration mode of equipment based on online vibration data, which solves the problem of the calculated vibration in the finite element modal analysis of the prior art, especially when the research object is a relatively complex system. The modes are often too messy, and it is difficult to identify the eigenfrequency of the component that is really concerned, and it is difficult to judge the value of the eigenfrequency.

The technical scheme adopted in the present invention is to carry out according to the following steps:

Step 1: Assume that the equipment vibration signal directly collected by the sensor is x(n), and its mathematical model is expressed as

x(n)=f(n)+noise(n)

Among them, f(n) is the original signal, noise(n) is the noise signal;

Step 2: Use the wavelet packet transform algorithm to denoise the equipment vibration signal x(n) once to obtain the denoised signal

Step 3: Use singular value decomposition to signal Perform secondary denoising to obtain the vibration signal after secondary denoising

Step 4: Use the windowed discrete Fourier algorithm to denoise the vibration signal Carry out spectrum analysis and calculate the equipment vibration spectrum;

Step 5: Use the self-learning algorithm to train to obtain the vibration spectrum of the equipment, and finally obtain the natural characteristic frequency and amplitude of the equipment.

Further, in the step 2, the wavelet packet transform algorithm performs a denoising process on the equipment vibration signal x(n) as follows:

①Wavelet packet decomposition of the signal; select a wavelet base, and determine the number of layers N of a wavelet packet decomposition, and then perform N-layer wavelet packet decomposition on the signal. The decomposition formula is:

In the formula: j<N, C _j,l (n) is the nth wavelet packet decomposition coefficient of the lth node in the jth layer in the decomposition tree, k and n both refer to the specific number of wavelet packets of a certain node decomposition factor;

Take the vibration signal x(n) as the signal to be analyzed, perform wavelet packet decomposition of the signal, and decompose the signal into different time-frequency spaces;

② Determining the best wavelet packet base: According to the minimum cost principle, calculate the best wavelet packet base for a given entropy standard;

Entropy is defined as:

M＝-∑P _jl lgP _jl

in, And when P _jl =0, P _jl lgP _jl =0;

③ Select a threshold value for processing the wavelet packet decomposition coefficients C _j,l under each decomposition scale, using the soft threshold method;

The soft threshold denoising method is defined as:

In the formula, λ is the threshold value, C _j,l is the wavelet packet coefficient, is the processed wavelet packet coefficient;

④Wavelet packet reconstruction: perform wavelet reconstruction according to the wavelet packet decomposition low-frequency coefficients and quantization processing coefficients of the Nth layer. The reconstruction formula is as follows:

In the formula, C _j+1,2l (n) is the wavelet packet decomposition system of the 2l node in the j+1th layer in the decomposition tree, and C _j+1,2l (n) is the j+1th layer in the decomposition tree The wavelet packet decomposition system of the 2l+1 node, h(n-2k) and g(n-2k) are the filter coefficients defined in the multi-resolution analysis, using the wavelet packet algorithm to analyze the equipment vibration signal x(n) Perform denoising to obtain a vibration signal after denoising

Further, in the step 3, the singular value decomposition is used to signal The specific steps for secondary denoising are as follows:

(1) from Extract the subsequence {x ₁ ,x ₂ ,L,x _n } from the time series as the first vector y ₁ of the n-dimensional phase space;

(2) Move one step to the right, extract {x ₂ ,x ₃ ,L,x _n+1 } as the first vector y ₂ of the n-dimensional phase space;

(3) By analogy, a set of column vectors {y ₁ , y ₂ , L, y _m } can be obtained;

(4) Each vector corresponds to a point in the reconstructed phase space, and all vectors form an m×n-dimensional matrix D _m :

D _m is the attractor orbital matrix with embedding dimension m and time delay of 1. If the measured vibration signal contains certain noise, then D _m =D+W, where D and W represent the smooth signal and the noise signal respectively. For the trajectory matrix in D _m , perform singular value decomposition on matrix D _m , D _m =USV', U and V are m×n and n×n order matrices respectively, and UU'=E, VV'=E, E is Identity matrix, S is an m×n diagonal matrix, diagonal elements s ₁ , s ₂ , L, s _p , p=min(m, n), s ₁ ≥ s ₂ ≥ L ≥ s _p ≥ 0, where s ₁ , s ₂ , L, s _p are the singular values of the matrix D _m , the first k (k≤p) items of the obtained singular values s ₁ , s ₂ , L, s _p are set to zero, and a new Diagonal matrix S', and then use the inverse process of SVD decomposition D _m '=US'V' to get a new matrix D _m ', matrix D _m 'is the best approximation matrix of D _m , according to the process of phase space reconstruction, Vibration signal after secondary denoising obtained from D _m '

Further, in the step 4, a triangular self-convolution window is used to perform windowed discrete Fourier algorithm to calculate the vibration spectrum.

Further, the self-learning algorithm steps in the step 5 are:

Step 1: Take the calculation cycle time as 1 second, determine the total training time as T calculation cycles, take t=1, and t represents the calculation sequence number;

Step 2: Input the vibration data x _t (n) of the equipment in the tth calculation period;

Step 3: Perform denoising, windowing, and discrete Fourier transform (FFT) on x _t (n), and calculate the frequency spectrum F _t (k), where F _t (k) represents the frequency spectrum in the tth calculation cycle of the device;

Step 4: sort F _t (k) from large to small in magnitude to obtain B _t (k);

Step 5: For B _t (k), set the threshold ρ, if the amplitude is greater than the threshold ρ, consider the frequency component corresponding to the amplitude to be the characteristic frequency, and store it in the characteristic matrix C _t ; if it is less than the threshold The value ρ is discarded;

Step 6: t=t+1, enter the next calculation cycle, judge t≤T, if yes, skip to step 2; otherwise, skip to step 7;

Step 7: The six steps of the above algorithm obtain the vibration characteristic parameter set C of each monitoring point during the training period, calculate the natural frequency of equipment vibration and its corresponding amplitude; draw the curve of each characteristic frequency component, and obtain the vibration signal of the equipment The natural frequency and its corresponding amplitude, thereby identifying the natural vibration mode of the device.

The beneficial effect of the invention is that when the research object is a relatively complex system, the characteristic frequency of the components can be accurately identified.

Description of drawings

Fig. 1 is the 3-layer decomposition tree structure schematic diagram of wavelet packet of the present invention;

Fig. 2 is the magnitude-frequency response schematic diagram of triangular window and triangular self-convolution window of the present invention;

Fig. 3 is the fast Fourier transform frame diagram of eight points of the present invention;

Fig. 4 is a flowchart of the self-learning (SLNM) algorithm of the present invention.

Detailed ways

The present invention will be described in detail below in combination with specific embodiments. The self-learning identification method of the natural vibration mode of the equipment based on the online vibration data of the present invention is carried out according to the following steps:

Step 1: The noise background at the working site of the equipment is complex, and the collected vibration signals are easily polluted by these noises, so that the real vibration data in the signal are submerged under the strong background noise, which brings great problems to the calculation of the natural frequency of the equipment. Difficulties. Assuming that the equipment vibration signal directly collected by the sensor is x(n), its mathematical model can be expressed as:

x(n)=f(n)+noise(n)

Among them, f(n) is the original signal; noise(n) is the noise signal.

Step 2: Use the wavelet packet transform algorithm to denoise the equipment vibration signal x(n) once to obtain the denoised signal The present invention adopts the denoising algorithm combining wavelet packet transform and singular value decomposition to denoise the vibration signal,

Wavelet packet transform: For a given orthogonal scaling function And its corresponding wavelet function φ(t), there is a double-scale equation

where: {h(n)} and {g(n)} are the filter coefficients defined in the multiresolution analysis, is the scaling function, and φ(t) is the wavelet function.

remember μ ₁ (t)＝φ(t), and μ _m (t) is defined by the recurrence relation as:

The family function {μ _l (t)|l＝0,1,2,L} is called relative to the orthogonal scaling function Orthogonal wavelet packets of , where l is a positive integer.

A group of orthogonal bases extracted from the wavelet packet that can form L ² (R) is called a wavelet packet base of L ² (R), and a set of norms extracted from the wavelet library that can form L ² (R) The orthogonal basis is the wavelet packet basis of L ² (R). The wavelet packet decomposition of a signal can be realized by using a variety of wavelet packet bases. Different wavelet packet bases have different decomposition results for the signal, and the results can reflect different degrees of the signal. Therefore, a set of optimal wavelet packets is sought. The base is an important task to express a signal most efficiently.

The specific steps of wavelet packet denoising:

①Wavelet packet decomposition of the signal. Select a wavelet base (namely wavelet function), and determine the number of layers N of a wavelet packet decomposition (generally N is 3 or 4), and then perform N-layer wavelet packet decomposition on the signal, and the decomposition formula is:

In the formula: j<N, C _j,l (n) is the nth wavelet packet decomposition coefficient of the lth node in the jth layer in the decomposition tree. Both k and n refer to the specific number of wavelet packet decomposition coefficients of a certain node.

The vibration signal x(n) is taken as the signal to be analyzed, and the wavelet packet decomposition of the signal is carried out to decompose the signal into different time-frequency spaces.

For example, a 3-layer wavelet packet decomposition is performed on the signal. The decomposition structure is shown in Figure 1 (the wavelet packet decomposition tree is a relatively intuitive representation of wavelet packet decomposition). In Figure 1, the (0,0) node represents the original signal, the (1,0) node represents the low-frequency component of the first layer of wavelet packet decomposition, and the (2,0) and (3,0) nodes represent the second layer, the component of the 0th node in layer 3, and so on.

② Determine the best wavelet packet base: according to the minimum cost principle, calculate the best wavelet packet base for a given entropy standard (using Shannon entropy);

Using wavelet packet analysis, a better wavelet packet base must be selected to describe the signal. In order to choose a better wavelet packet base, firstly given a sequence of cost functions, find the base that minimizes the cost function among all wavelet packet bases, and minimize the Shannon entropy (cost function) for a given signal The basis of is the most effective wavelet packet basis to represent the signal, and this basis is called the best basis.

Shannon entropy (cost function) is defined as:

M＝-∑P _jl lgP _jl

in, And when P _jl =0, P _jl lgP _jl =0.

Step ② of the present invention is carried out on the basis of step ①. First, step ① performs multiple decompositions to calculate the decomposition coefficient C _j,l of the signal under different wavelet bases. The second step utilizes the decomposition obtained in step ①. The coefficient C _j,l calculates the shonnon entropy value, finds the minimum entropy value, and its corresponding wavelet basis is the optimal wavelet basis.

③ An appropriate threshold is selected for processing the wavelet packet decomposition coefficients C _j,l under each decomposition scale, and the present invention adopts a soft threshold method.

The soft threshold denoising method is defined as:

In the formula, λ is the threshold value, C _j,l is the wavelet packet coefficient, is the processed wavelet packet coefficient.

④Wavelet packet reconstruction, according to the wavelet packet decomposition of the Nth layer low-frequency coefficients and quantized processing coefficients for wavelet reconstruction. The reconstruction formula is as follows:

In the formula, C _j+1,2l (n) is the wavelet packet decomposition system of the 2l node in the j+1th layer in the decomposition tree, and C _j+1,2l (n) is the j+1th layer in the decomposition tree The wavelet packet decomposition system of the 2l+1 node, h(n-2k) and g(n-2k) are the filter coefficients defined in the multi-resolution analysis.

Use the wavelet packet algorithm to denoise the vibration signal x(n) of the equipment to obtain a denoised vibration signal

Step 3: Singular Value Decomposition (Single Value Decomposition, SVD) pair Perform secondary denoising:

After wavelet packet denoising is performed on the vibration signal of the equipment, a denoised vibration signal is obtained The present invention uses the singular value decomposition pair The secondary noise reduction method is based on the different effects of the smooth system signal and the random noise signal on the singular value of the reconstructed attractor orbital matrix. This method only removes the noise signal and has almost no effect on the system signal. Firstly, the vibration signal is reconstructed in the phase space, and the reconstructed attractor can reflect the dynamic characteristics of the system. In order to characterize the orbital matrix of the attractor, the singular value decomposition is performed, and the characteristics of the singular spectrum are used to reduce the noise component in the vibration signal. Can achieve better results.

vibration signal is a discrete time series, the phase space reconstruction is performed on it, and the specific steps are as follows:

(1) from Extract the subsequence {x ₁ ,x ₂ ,L,x _n } from the time series as the first vector y ₁ of the n-dimensional phase space;

(2) Move one step to the right, extract {x ₂ ,x ₃ ,L,x _n+1 } as the first vector y ₂ of the n-dimensional phase space;

(3) By analogy, a set of column vectors {y ₁ , y ₂ , L, y _m } can be obtained;

(4) Each vector corresponds to a point in the reconstructed phase space, and all vectors form an m×n-dimensional matrix D _m :

D _m is the attractor track matrix with embedding dimension m and time delay of 1. If the measured vibration signal contains certain noise, then D _m =D+W, where D and W respectively denote the trajectory matrix in D _m corresponding to the smooth signal and the noise signal. Singular value decomposition is performed on the matrix D _m , D _m =USV', U and V are matrices of order m×n and n×n respectively, and UU'=E, VV'=E (E is the unit matrix). S is an m×n diagonal matrix, diagonal elements s ₁ , s ₂ , L, s _p , p=min(m, n), s ₁ ≥ s ₂ ≥ L ≥ s _p ≥ 0, where s ₁ , s ₂ , L, sp the _singular value of matrix D _m . In the singular value decomposition noise reduction process, the selection of the noise reduction order is very critical, and the difference in the selected order will cause the noise reduction effect to be significantly different. When the selected order is too high, the noise reduction signal will still contain more noise information, which cannot achieve a good noise reduction effect; when the selected order is too low, the information after noise reduction will be incomplete, or even cause waveform distortion. Here, the first _k (k≤p) items of the obtained singular values s ₁ , s ₂ , L, sp are set to zero, and a new diagonal matrix S' is obtained. Then use the inverse process of SVD decomposition D _m '= _US'V ' to get a new matrix D _m ', the matrix D _m 'is the best approximation matrix of D _m , according to the process of phase space reconstruction, get the de Vibration signal after noise

Complete the singular value decomposition pair After denoising, the vibration signal after the second denoising is obtained

Step 4: Windowed discrete Fourier algorithm calculates the vibration spectrum:

The present invention uses the windowed discrete Fourier algorithm to denoise the vibration signal Spectrum analysis is carried out to calculate the vibration spectrum of the equipment.

Spectrum analysis is an analysis that converts the time-domain signal with time as the abscissa into a frequency-domain signal with frequency as the abscissa through Fourier transform, and then obtains the phase and amplitude of the frequency components of the original time-domain signal method. The graphs drawn by frequency and phase or amplitude are called phase spectrum and amplitude spectrum respectively, and the amplitude spectrum is more commonly used in engineering. In order to improve the efficiency of calculating frequency spectrum, the present invention uses Fast Fourier Transform (FFT) to process the vibration signal.

The calculation formula of discrete Fourier transform (DFT):

In the formula, is a time-domain signal with a data length of n; F(k) is the result of transformation into the complex domain, which is the same length as the time-domain signal.

Windowing principle: Since the discrete Fourier transform is defined for finite length sequences, therefore, for infinite length sequences When calculating the spectrum F(k), x(n) must be truncated or segmented, which is equivalent to Multiplying with the rectangular sequence w _N (n)=u(n)-u(nN) whose amplitude is 1 and whose length is N, the truncated N-point sequence x _N (n) is:

This rectangular sequence w _N (n) with an amplitude of 1 and a length of N is the rectangular window function. This kind of truncation or segmentation of the signal is windowing. The function selected during windowing processing is different, which affects the spectrum analysis. Also different.

Spectrum leakage: The so-called spectrum leakage refers to the phenomenon that the signal energy of a certain frequency spreads to adjacent frequency points during discrete Fourier transform. Discrete Fourier transform is performed on the measured signal, and time-domain windowing of the signal is equivalent to frequency-domain convolution. This truncation causes errors in the spectrum, and its effect makes the spectrum centered on the actual frequency value, and the shape of the window function spectrum waveform Diffusion to both sides, resulting in a "leakage effect". Leakage effects add new frequency components and change the magnitude of the spectral values. From an energy point of view, the frequency leakage phenomenon is equivalent to the energy at various frequency components of the original signal penetrating into other frequency components, so it is also called power leakage.

Window function selection: Spectrum leakage is inherent in the discrete Fourier transform, and it is closely related to the side lobes of the window function spectrum. If the height of the side lobes tends to zero, so that the energy is relatively concentrated on the main lobe, a relatively spectrum close to the true value. Considering the characteristics of the vibration signal, the present invention uses a triangular self-convolution window.

Figure 2 below shows the normalized amplitude characteristics of the triangular window and the triangular self-convolution window (the dotted line in the figure is the triangular window, and the solid line is the triangular self-convolution window).

It can be seen from Figure 2 that the largest side lobe of the triangular window is 26dB lower than the main lobe, while the largest side lobe of the triangular self-convolution window is 52dB lower than the main lobe, and the attenuation rate of the side lobe is significantly accelerated. The comparison shows that the triangular self-convolution The energy of the product window is mainly concentrated in the main lobe, and the side lobes are greatly reduced. The triangular self-convolution window can effectively suppress the influence of spectrum leakage.

FFT calculation spectrum: When using fast Fourier transform for spectrum analysis, two important parameters need to be determined: sampling rate f _s and frequency domain sampling points N _fft , the sampling rate is based on the highest frequency of the signal, according to the Nyquist sampling theorem To determine, the number of sampling points is determined according to the frequency resolution Vf.

but

The operation rules of the FFT time extraction algorithm are summarized as follows:

The Fast Fourier Algorithm (FFT) mainly includes two processes. First, the windowed time sequence f _N (n) is code-bit inverted to obtain a new input sequence then to Perform butterfly operation to obtain spectrum F(k). The operation process of the Fast Fourier Algorithm is shown in Figure 3 (take 8 points as an example). Firstly, the sequences f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7) perform code bit inversion to obtain new sequences f(0), f(1), f(4), f(2), f(1), f( 5), f(3), f(7), and then perform three-level butterfly operation on the new sequence to obtain F(0), F(1), F(2), F(3), F(4) , F(5), F(6), F(7).

(1) code point inversion

In order to arrange the output sequence F(k) in a natural order, the input sequence f _N (n) must be arranged in the code position inversion order, the algorithm is as follows:

Let I be a sequential binary number, J be a reverse binary number, I=J=000B. Let the constant N ₂ be half of the number of points, N _fft /2 as the midpoint of the sequence, and set the constant N ₁ as the total number of points of the sequence N _fft -1 as the length of the sequence. First judge I≥J, if I=J=000B is satisfied, reverse order is required, and the index exchange part is skipped. Let the variable K=N ₂ =100B, judge K>J, then J=J+K=000+100=100B, I=I+1=001B. Then judge I≥J, otherwise enter the index exchange part: exchange the values in the storage units of I and J at this time. Then judge K>J at this time, otherwise J=JK=100-100=000B, K=K/2=010B, then judge K>J at this time, then J=J+K=000+010=010B , I=I+1=001+1=010B, judging that I has not reached the total number of sequence points N ₁ =N _fft -1, then return to I≥J, and so on, cyclic judgment and calculation, until I reaches the total number of points Then exit the program.

(2) Butterfly operation

All operations are decomposed into M-level operations, and each level of operations includes N/2 basic butterfly operations. In the butterfly diagram, the butterfly operations that intersect each other at any level are called groups; in the butterfly diagram, in the order from top to bottom, the increment of the corresponding element numbers between the upper and lower groups is called group interval. The L-th stage operation includes N/2 ^L groups, and the L-th stage group interval is 2 ^L . The W _N distribution of each group in the same level is the same, there are 2 ^L-1 W _N multipliers in the L level, and the W _N multipliers in each group are distributed according to the following rules from top to bottom:

First level:

second level:

Level L: L,

Class M: L,

The input sequence interval of the L-level basic butterfly operation is 2 ^L-1 , and the operational relationship between its input and output is:

In the formula, AL (·) represents the input element of the _L -th basic butterfly operation, and _AL (·) represents the output element of the L-th basic butterfly operation.

For the vibration signal after secondary denoising Perform windowed discrete Fourier transform to obtain the vibration signal spectrum F(k) of the equipment, which is used as the basis for the next step of the self-learning algorithm.

The self-learning algorithm (Self-Learning Natural Mode, SLNM) of the equipment vibration natural mode is shown in Figure 4. The whole process of calculating the frequency spectrum of the vibration signal was introduced above. In order to identify the natural vibration mode of the equipment, it is necessary to monitor the vibration signal of the equipment during normal operation for a long time, take the vibration signal of the equipment within a reasonable period of time, and denoise and add a window to the vibration signal Discrete Fourier transform, calculate the frequency spectrum of the equipment, use the self-learning algorithm to train the frequency spectrum in this time period, and finally obtain the natural characteristic frequency and amplitude of the equipment. The steps of the self-learning algorithm are as follows:

Step 1: Take the calculation cycle time as 1 second, determine the total training time as T calculation cycles, take t=1, and t represents the calculation sequence number;

Step 2: Input the vibration data x _t (n) of the equipment in the tth calculation period;

Step 3: Perform denoising, windowing, and discrete Fourier transform (FFT) on x _t (n), and calculate the frequency spectrum F _t (k), where F _t (k) represents the frequency spectrum in the tth calculation period of the device;

Step 4: sort F _t (k) from large to small in magnitude to obtain B _t (k);

Step 5: For B _t (k), set the threshold ρ, if the amplitude is greater than the threshold ρ, consider the frequency component corresponding to the amplitude to be the characteristic frequency, and store it in the characteristic matrix C _t ; if it is less than the threshold The value ρ is discarded;

Step 6: t=t+1, enter the next calculation cycle, judge t≤T, if yes, skip to step 2; otherwise, skip to step 7;

Step 7: The six steps of the above algorithm obtain the vibration characteristic parameter set C of each monitoring point during the training period, calculate the natural frequency of equipment vibration and its corresponding amplitude (including: mean value, variance); and draw the curve of each characteristic frequency component .

After completing the steps of the self-learning algorithm (SLNM), the natural frequency of the vibration signal of the equipment and its corresponding amplitude are calculated, thereby identifying the natural vibration mode of the equipment.

Calculation analysis verification

By installing a vibration sensor on the wind turbine gearbox, the vibration data of the wind turbine gearbox is measured in real time, and the invention analyzes and processes the data to calculate the vibration characteristic parameters of the wind turbine gearbox.

The self-learning time is set to 4 hours, and the vibration data of 1 second is selected as a sample every 5 minutes. There are 48 samples in total, and the 48 samples are denoised and discrete Fourier transformed to calculate the frequency spectrum, and then use the self-learning The learning (SLNM) algorithm is used to process the vibration characteristic parameters of the gearbox, so as to identify the natural vibration mode of the gearbox. Table 1 shows the vibration frequency and amplitude of wind turbine gearbox.

Table 1

Note: f1~f6 in Table 1 represent the 6 vibration characteristic frequencies, and X1~X6 represent the corresponding amplitudes of f1~f6 respectively; "/" means that there is no such frequency component in the spectrum.

It can be clearly seen from Table 1 that the frequency spectrum of the wind turbine gearbox has been changing, and the change range is relatively large. In order to identify the natural vibration mode of the gearbox, 48 frequency spectra were statistically analyzed and the characteristic parameters were calculated. The results are shown in Table 2 as the vibration characteristic parameters of the wind turbine gearbox.

Table 2

Table 2 shows the vibration characteristic parameters of the gearbox. The six characteristic frequency components are: (622,0.058), (1495,0.067), (1213,0.04), (2075,0.04), (1825,0.044), ( 3270,0.024).

(1) What vibration signal denoising of the present invention adopts is the denoising algorithm that wavelet packet combines with singular value decomposition, first carries out wavelet packet denoising to vibration signal, then carries out singular value decomposition secondary noise reduction;

(2) In order to accurately calculate the frequency spectrum of the vibration signal and reduce the impact of spectrum leakage, the present invention first carries out windowing processing on the vibration signal, and then utilizes fast Fourier transform to calculate the frequency spectrum;

(3) In order to calculate the characteristic parameters of the natural vibration of the equipment, the present invention proposes a self-learning (SLNM) algorithm.

Advantage of the present invention mainly is divided into three aspects, is summarized as follows:

Accuracy: The improved triangular window function is used to calculate the vibration signal spectrum, which can greatly reduce the influence of spectrum leakage and obtain a spectrum closer to the real value;

Stability: the present invention uses a denoising algorithm combining wavelet packets and singular value degradation, which can effectively remove noise signals and is less disturbed by the surrounding environment;

Versatility: the present invention can be used for natural vibration pattern recognition of various devices.

The above description is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Any simple modifications made to the above embodiments according to the technical essence of the present invention, equivalent changes and modifications, all belong to this invention. within the scope of the technical solution of the invention.

Claims (1)

1. The device inherent vibration mode self-learning identification method based on the online vibration data is characterized by comprising the following steps of:

step 1: let x (n) be the vibration signal of the device directly collected by the sensor, and its mathematical model is expressed as x (n) = f (n) + noise (n)

Wherein f (n) is an original signal, and noise (n) is a noise signal;

and 2, step: carrying out primary denoising on the equipment vibration signal x (n) by utilizing a wavelet packet transformation algorithm to obtain a denoised signal

The process of carrying out primary denoising on the equipment vibration signal x (n) by the wavelet packet transformation algorithm comprises the following steps:

(1) wavelet packet decomposition of the signal; selecting a wavelet base, determining the number N of layers of wavelet packet decomposition, and then performing N-layer wavelet packet decomposition on the signal, wherein the decomposition formula is as follows:

in the formula: j is a function of<N，C _j,l (n) is the nth wavelet packet decomposition coefficient of the jth node of the jth layer in the decomposition tree, and both k and n refer to the specific several wavelet packet decomposition coefficients of a certain node;

taking a vibration signal x (n) as a signal to be analyzed, carrying out wavelet packet decomposition on the signal, and decomposing the signal into different time-frequency spaces;

(2) determining an optimal wavelet packet basis: calculating the optimal wavelet packet basis for a given entropy standard according to the minimum cost principle;

entropy is defined as:

M＝-∑P _jl lgP _jl

wherein,and when P is _jl When =0, P _jl lgP _jl ＝0；

(3) For each decomposition scaleWavelet packet decomposition coefficient C _j,l Selecting a threshold value for processing, and adopting a soft threshold value method;

the soft threshold noise cancellation method is defined as:

wherein λ is a threshold value, C _j,l Is the coefficient of the wavelet packet and is,the processed wavelet packet coefficient;

(4) wavelet packet reconstruction, namely performing wavelet reconstruction according to the wavelet packet decomposition low-frequency coefficient and the quantization processing coefficient of the Nth layer, wherein the reconstruction formula is as follows:

in the formula, C _j+1,2l (n) wavelet packet decomposition system for the 2l node at the j +1 th level in the decomposition tree, C _j+1,2l (n) is a wavelet packet decomposition system of the 2l +1 node of the j +1 layer in the decomposition tree, h (n-2 k) and g (n-2 k) are filter coefficients defined in multi-resolution analysis, and a device vibration signal x (n) is denoised by using a wavelet packet algorithm to obtain a vibration signal after primary denoising

And 3, step 3: using singular value decomposition on the signalCarrying out secondary denoising to obtain a vibration signal subjected to secondary denoising

(1) FromExtracting subsequences { x ] in time series ₁ ,x ₂ ,L,x _n As the first vector y of the n-dimensional phase space ₁ ；

(2) Move by one step to the right, decimate { x ₂ ,x ₃ ,L,x _n+1 As the first vector y of the n-dimensional phase space ₂ ；

(3) By analogy, a group of column vectors y is obtained ₁ ,y ₂ ,L,y _m }；

(4) Each vector corresponds to a point in the reconstruction phase space, and all vectors form a matrix D with dimensions of m × n _m ：

D _m For embedding attractor orbit matrix with dimension m and time delay of 1, if the measured vibration signal contains certain noise, D _m = D + W, where D, W represent D for smooth and noisy signals, respectively _m Track matrix of (1), to matrix D _m By singular value decomposition, D _m (= USV '), U and V are m × n and n × n order matrices, respectively, and UU ' = E, VV ' = E, E is an identity matrix, S is an m × n diagonal matrix, a diagonal element S ₁ ,s ₂ ,L,s _p ,p＝min(m,n),s ₁ ≥s ₂ ≥L≥s _p Is more than or equal to 0, wherein s ₁ ,s ₂ ,L,s _p Is a matrix D _m The obtained singular value s ₁ ,s ₂ ,L,s _p The first k (k is less than or equal to p) item, and other items are set to zero to obtain a new diagonal matrix S', and then the inverse process D of SVD decomposition is utilized _m ' = US ' V ' resulting in a new matrix D _m ', matrix D _m Is' a D _m According to the process of phase space reconstruction, from D _m ' obtaining vibration signal after secondary denoising

And 4, step 4: using windowed discrete Fourier algorithm to removePost-noise vibration signalCarrying out frequency spectrum analysis, and calculating to obtain a device vibration frequency spectrum;

and 5: training by using a self-learning algorithm to obtain a vibration frequency spectrum of the equipment, and finally obtaining the inherent characteristic frequency and amplitude of the equipment;

the self-learning algorithm comprises the following steps:

the method comprises the following steps: taking the time of a calculation period as 1 second, determining the total training time as T calculation periods, and taking T =1,t to represent a calculation sequence number;

step two: inputting vibration data x of the equipment in the t-th calculation period _t (n)；

Step three: for x _t (n) performing denoising, windowing and discrete Fourier transform (FFT), and calculating to obtain a frequency spectrum F _t (k)，F _t (k) Representing the frequency spectrum of the device in the t calculation period;

step four: to F _t (k) Sequencing according to the amplitude from large to small to obtain B _t (k)；

Step five: to B _t (k) Setting a threshold value rho, if the amplitude is greater than the threshold value rho, considering the frequency component corresponding to the amplitude as the characteristic frequency, and storing the characteristic frequency into a characteristic matrix C _t The preparation method comprises the following steps of (1) performing; if the value is smaller than the threshold value rho, the value is discarded;

step six: t = T +1, entering the next calculation period, judging that T is less than or equal to T, and if so, jumping to the step two; otherwise, jumping to the seventh step;

step seven: the six steps of the algorithm obtain a vibration characteristic parameter set C of each monitoring point in a training period, and calculate the natural frequency of the equipment vibration and the corresponding amplitude thereof; and drawing curves of all characteristic frequency components to obtain the natural frequency and the corresponding amplitude of the equipment vibration number, thereby identifying the natural vibration mode of the equipment.