CN104807534B - Equipment eigentone self study recognition methods based on on-line vibration data - Google Patents
Equipment eigentone self study recognition methods based on on-line vibration data Download PDFInfo
- Publication number
- CN104807534B CN104807534B CN201510024313.5A CN201510024313A CN104807534B CN 104807534 B CN104807534 B CN 104807534B CN 201510024313 A CN201510024313 A CN 201510024313A CN 104807534 B CN104807534 B CN 104807534B
- Authority
- CN
- China
- Prior art keywords
- vibration
- signal
- wavelet packet
- decomposition
- equipment
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Complex Calculations (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention discloses the equipment eigentone self study recognition methods based on on-line vibration data, carries out a denoising to equipment vibrating signal x (n) using wavelet package transforms algorithm, obtains the signal after denoisingUsing singular value decomposition to signalSecond denoising is carried out, obtains the vibration signal after second denoisingUsing adding window discrete fourier algorithm to the vibration signal after denoisingSpectrum analysis is carried out, vibration equipment frequency spectrum is calculated；Train to obtain vibration equipment frequency spectrum using self-learning algorithm, finally obtain the inherent feature frequency and amplitude of equipment.The beneficial effects of the invention are as follows when research object is complex system, can accurate recognition go out characteristic frequency of parts.
Description
Technical Field
The invention belongs to the technical field of online monitoring, and relates to a self-learning identification method for an inherent vibration mode of equipment based on online vibration data.
Background
On-line monitoring technology is the most basic measure for estimating the health condition of equipment, and especially for mechanical equipment (or systems), vibration monitoring analysis is the most practical and effective method. The vibration signal of the equipment is closely related to the mechanical structure of the equipment, and the condition of the internal mechanical structure of the equipment can be directly and effectively reflected from the vibration signal in the operation process of the equipment. The vibration signals of the equipment are processed and analyzed, and the inherent vibration mode of the equipment under the normal operation condition is extracted and used for distinguishing the vibration mode when the equipment is abnormal or has a fault, so that a decision can be made in time when an abnormal vibration model appears, and a corresponding maintenance plan can be made. The invention provides a vibration data-based automatic learning, identifying and extracting method for the natural vibration mode of equipment.
The natural vibration mode (mode) of the device refers to the natural vibration characteristics of the mechanical structure, including the natural frequency, damping ratio and mode shape of the device. If the characteristics of the main modes of the mechanical structure of the device are known, the actual vibration response of the structure under the action of various vibration sources can be predicted. Therefore, modal analysis, i.e. analyzing the modalities to obtain corresponding modal parameters, is an important method for structure dynamic design and fault diagnosis of equipment. Although the vibration behavior of actual equipment (such as wind power generators, transformers, ships, aircrafts and the like) varies and is complex, modal analysis provides an effective way to research the actual structural characteristics of various equipment.
The existing equipment inherent vibration mode identification method is mainly a computational modal analysis method, and firstly, a high-precision three-dimensional entity model of equipment is constructed by utilizing three-dimensional modeling software; then, importing the three-dimensional model of the equipment into finite element analysis software, setting material attribute parameters of the equipment according to the actual condition of the equipment, dividing the three-dimensional model into a finite element network, and adding displacement constraint to obtain a finite element model of the equipment; and finally, carrying out modal analysis on the finite element model to obtain the natural frequency and the mode of the model.
In finite element modal analysis in the prior art, particularly when a research object is a complex system, a calculated vibration mode is often too messy, and the characteristic frequency of a part really concerned is difficult to identify; or when a certain order of vibration is composed of many local modes, it is very difficult to judge the value of the characteristic frequency.
According to the method, a complex model does not need to be established, the modal parameters are obtained by directly acquiring the vibration data of different parts of the equipment through analysis and processing, and the inherent vibration mode of the equipment is identified.
Disclosure of Invention
The invention aims to provide an equipment inherent vibration mode self-learning identification method based on online vibration data, and solves the problems that in finite element modal analysis in the prior art, particularly when a research object is a complex system, a calculated vibration mode is often too messy, the characteristic frequency of a part really concerned is difficult to identify, and the numerical value of the characteristic frequency is difficult to judge.
The technical scheme adopted by the invention is carried out according to the following steps:
step 1: let the vibration signal of the device directly picked up by the sensor be x (n), and its mathematical model be expressed as
x(n)＝f(n)+noise(n)
Wherein f (n) is an original signal, and noise (n) is a noise signal;
step 2: carrying out primary denoising on the equipment vibration signal x (n) by utilizing a wavelet packet transformation algorithm to obtain a denoised signal
And step 3: using singular value decomposition on the signalCarrying out secondary denoising to obtain a vibration signal subjected to secondary denoising
And 4, step 4: de-noised vibration signal using windowed discrete Fourier algorithmCarrying out frequency spectrum analysis, and calculating to obtain a device vibration frequency spectrum;
and 5: and training by utilizing a self-learning algorithm to obtain the vibration frequency spectrum of the equipment, and finally obtaining the inherent characteristic frequency and amplitude of the equipment.
Further, the process of performing primary denoising on the equipment vibration signal x (n) by the wavelet packet transform algorithm in the step 2 is as follows:
(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 an optimal wavelet packet basis for a given entropy standard according to a minimum cost principle;
entropy is defined as:
M＝-∑P jl lgP jl
wherein the content of the first and second substances,and when P is jl When =0, P jl lgP jl ＝0；
(3) For wavelet packet decomposition coefficient C at each decomposition scale 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,for processed wavelet packetsA 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
Further, the singular value decomposition is used for the signals in the step 3The secondary denoising specifically comprises the following steps:
(1) FromExtracting subsequences { x ] in time series 1 ,x 2 ,L,x n As the first vector y of the n-dimensional phase space 1 ；
(2) Move by one step to the right, decimate { x 2 ,x 3 ,L,x n+1 As the first vector y of the n-dimensional phase space 2 ；
(3) By analogy, a group of column vectors y can be obtained 1 ,y 2 ,L,y m }；
(4) Each vector corresponds to a point in the reconstruction phase space, all vectors forming a matrix D of dimensions m x n m ：
D m For embedding an attractor track 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 Of the trajectory matrix, pair 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 1 ,s 2 ,L,s p ,p＝min(m,n),s 1 ≥s 2 ≥L≥s p Not less than 0, wherein s 1 ,s 2 ,L,s p Is a matrix D m The obtained singular value s 1 ,s 2 ,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
Further, in the step 4, a triangular self-convolution window is adopted to perform windowing discrete Fourier algorithm to calculate the vibration frequency spectrum.
Further, the self-learning algorithm in the step 5 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-th 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 Performing the following steps; 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.
The method has the advantage that when the research object is a complex system, the characteristic frequency of the part can be accurately identified.
Drawings
FIG. 1 is a schematic diagram of a wavelet packet 3-level decomposition tree structure according to the present invention;
FIG. 2 is a schematic diagram of the amplitude-frequency response of the triangular window and the triangular self-convolution window of the present invention;
FIG. 3 is a block diagram of an eight-point fast Fourier transform framework of the present invention;
FIG. 4 is a flow chart of the self-learning (SLNM) algorithm of the present invention.
Detailed Description
The present invention will be described in detail with reference to the following embodiments. The invention discloses a self-learning identification method of equipment natural vibration modes based on online vibration data, which comprises the following steps of:
step 1: the noise background of the working site of the equipment is complex, and the collected vibration signals are very easily polluted by the noise, so that the real vibration data in the signals are submerged under strong background noise, and great difficulty is brought to the calculation of the natural frequency of the equipment. Let the device vibration signal directly acquired by the sensor be x (n), and its mathematical model can be expressed as:
x(n)＝f(n)+noise(n)
wherein f (n) is an original signal; 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 signalThe invention adopts a denoising algorithm combining wavelet packet transformation and singular value decomposition to denoise the vibration signal,
wavelet packet transformation: for a given orthogonal scale functionAnd its corresponding wavelet function phi (t), there is a dual-scale equation
In the formula: { h (n) } and { g (n) } are filter coefficients defined in the multiresolution analysis,phi (t) is a scale function and phi (t) is a wavelet function.
Note bookμ 1 (t) = phi (t), mu is defined by recursion relation m (t) is:
family function [ mu ] l (t) | l =0,1,2, L } is relative to an orthogonal scale functionWherein l is a positive integer.
Energy extracted from the wavelet packet constitutes L 2 One set of orthogonal bases of (R) is referred to as L 2 (R) a wavelet basis, the energy extracted from the wavelet bank constituting L 2 (R) a set of orthonormal bases L 2 Wavelet basis of (R). The wavelet packet decomposition of a signal can be realized by adopting various wavelet packet bases, different wavelet packet bases have different decomposition results on the signal, and the results can reflect different degrees of the signal, so that the search of an optimal wavelet base is an important task for most effectively expressing a signal.
The specific steps of wavelet packet denoising are as follows:
(1) wavelet packet decomposition of the signal. Selecting a wavelet basis (namely a wavelet function), determining the layer number N (generally N is 3 or 4) of a 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 unit of a group<N，C j,l And (n) is the nth wavelet packet decomposition coefficient of the jth node of the jth layer in the decomposition tree. k and n both refer to the specific number of wavelet packet decomposition coefficients of a node.
And taking the 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.
For example, a 3-layer wavelet packet decomposition is performed on the signal. The decomposition structure is shown in fig. 1 (wavelet packet decomposition tree is a more intuitive representation of wavelet packet decomposition). In fig. 1, the node (0, 0) represents the original signal, the node (1, 0) represents the low-frequency component of layer 1 of the wavelet packet decomposition, and the nodes (2, 0) and (3, 0) represent the components of the node of layer 2, layer 3, and node 0, respectively, and so on.
(2) Determining an optimal wavelet packet basis: calculating the optimal wavelet packet basis for a given entropy standard (adopting Shannon entropy) according to the minimum cost principle;
with wavelet packet analysis, a better wavelet packet basis must be selected to describe the signal. To select a better wavelet packet basis, a sequence of cost functions is first given, and the basis that minimizes the cost function is found among all wavelet packet bases, and for a given signal, the basis that minimizes Shannon entropy (cost function) is the wavelet packet basis that most effectively represents the signal, and this basis is called the best basis.
The Shannon entropy (cost function) is defined as:
M＝-∑P jl lgP jl
wherein the content of the first and second substances,and when P is jl When =0, P jl lgP jl ＝0。
The step (2) of the invention is carried out on the basis of the step (1), firstly, the step (1) is carried out for a plurality of times of decomposition, and the decomposition coefficients C of the signals under different wavelet bases are calculated j,l The second step uses the decomposition coefficient C obtained in the step (1) j,l And calculating the entropy value of the shonnon to find the minimum entropy value, wherein the corresponding wavelet basis is the optimal wavelet basis.
(3) For wavelet packet decomposition coefficient C at each decomposition scale j,l An appropriate threshold value is selected for processing, and the method adopts 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,is the processed wavelet packet coefficient.
(4) And 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. 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 the wavelet packet decomposition system of node 2l +1 at level j +1 in the decomposition tree, and h (n-2 k) and g (n-2 k) are filter coefficients defined in the multiresolution analysis.
Denoising the equipment vibration signal x (n) by using a wavelet packet algorithm to obtain a vibration signal subjected to primary denoising
And step 3: singular Value Decomposition (SVD) pairCarrying out secondary denoising:
after wavelet packet denoising is carried out on the vibration signal of the equipment, the vibration signal subjected to primary denoising is obtainedThe invention uses singular value decomposition pairsThe secondary noise reduction method is based on the fact that different influences of smooth system signals and random noise signals on singular values of a reconstructed attractor track matrix are subjected to noise reduction, only the noise signals are removed, and the influence on the system signals is almost eliminated. Firstly, the vibration signal is in phase spaceAnd when reconstruction is carried out, the reconstructed attractor can reflect the dynamic characteristics of the system, singular value decomposition is carried out on the orbit matrix representing the attractor, the noise component in the vibration signal is reduced by utilizing the characteristics of a singular spectrum, and a better effect can be achieved.
Vibration signalIs a discrete time sequence, and carries out phase space reconstruction on the discrete time sequence, and the specific steps are as follows:
(1) FromExtracting subsequences { x ] in time series 1 ,x 2 ,L,x n As the first vector y of the n-dimensional phase space 1 ；
(2) Move by one step to the right, decimate { x 2 ,x 3 ,L,x n+1 As the first vector y of the n-dimensional phase space 2 ；
(3) By analogy, a group of column vectors y can be obtained 1 ,y 2 ,L,y m }；
(4) Each vector corresponds to a point in the reconstruction phase space, all vectors forming a matrix D of dimensions m x n m ：
D m To embed dimension m, an attractor track matrix with a 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 The trajectory matrix of (2). For 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 mxn diagonal matrix, diagonal element S 1 ,s 2 ,L,s p ,p＝min(m,n),s 1 ≥s 2 ≥L≥s p Not less than 0, wherein s 1 ,s 2 ,L,s p Matrix D m The singular value of (c). In the singular value decomposition noise reduction process, the selection of the noise reduction order is very critical, and the difference of the selection order can cause the obvious difference of the noise reduction effect. When the selection order is too high, the noise reduction signal still contains more noise information, and a better noise reduction effect cannot be achieved; if the selected order is too low, the information after noise reduction may be incomplete, and even the waveform may be distorted. The singular value s to be obtained here 1 ,s 2 ,L,s p The front k (k is less than or equal to p) items of the matrix are set to zero, and a new diagonal matrix S' is obtained. Reuse of the inverse of SVD decomposition 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 denoised vibration signals
Complete singular value decomposition pairAfter denoising, obtaining a vibration signal after secondary denoising
And 4, step 4: calculating a vibration frequency spectrum by using a windowed discrete Fourier algorithm:
the invention uses windowing discrete Fourier algorithm to denoise vibration signalsAnd carrying out spectrum analysis and calculating to obtain the vibration spectrum of the equipment.
Spectral analysis is an analysis method in which a time domain signal whose abscissa is time is transformed into a frequency domain signal whose abscissa is frequency by fourier transform, and the phase and amplitude of the frequency component of the original time domain signal are obtained. 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 the frequency spectrum, the invention adopts Fast Fourier Transform (FFT) to process the vibration signal.
Discrete Fourier Transform (DFT) calculation:
in the formula (I), the compound is shown in the specification,is a time domain signal of data length n; f (k) is the result of the transformation into the complex domain, the same length as the time domain signal.
The principle of windowing: since the discrete fourier transform is defined for a finite length sequence, it is possible to use an infinite length sequenceWhen calculating the spectrum F (k), x (n) has to be truncated or segmented, which corresponds to dividing the spectrum F (k) into segmentsAnd a rectangular sequence w of amplitude 1 and length N N (N) = u (N) -u (N-N) multiplication, and N-point sequence x after truncation N (n) is:
this rectangular sequence w of amplitude 1 and length N N The (n) is a rectangular window function, the truncation or segmentation of the signal is windowing, and different functions are selected during windowing, so that different influences are brought to spectrum analysis.
Spectrum leakage: the spectral leakage is a phenomenon in which signal energy of a certain frequency spreads to adjacent frequency points when discrete fourier transform is performed. The method is characterized in that discrete Fourier transform is carried out on a measured signal, windowing of a signal time domain is equivalent to frequency domain convolution, the truncation causes errors of a frequency spectrum, the effect is that the frequency spectrum takes an actual frequency value as a center, the shape of a window function frequency spectrum waveform is diffused to two sides, and a leakage effect is generated. The leakage effect adds new frequency components and changes the spectral value size. From the energy perspective, the frequency leakage phenomenon is equivalent to the energy at various frequency components of the original signal permeating to other frequency components, so it is also called power leakage.
Selecting a window function: the frequency spectrum leakage is inherent to the discrete Fourier transform and is closely related to the side lobe of the window function spectrum, and if the height of the side lobe is enabled to approach zero, so that the energy is relatively concentrated in the main lobe, the frequency spectrum which is closer to the true value can be obtained. In consideration of the characteristics of the vibration signal, the invention adopts a triangular self-convolution window.
The normalized amplitude characteristics of the triangular window and the triangular self-convolution window are given in the following fig. 2 (the dotted line in the figure is the triangular window, and the solid line is the triangular self-convolution window).
As can be seen from FIG. 2, the maximum side lobe of the triangular window is 26dB lower than the main lobe, while the maximum side lobe of the triangular self-convolution window is 52dB lower than the main lobe, and the side lobe attenuation rate is obviously accelerated.
FFT calculates the spectrum: in spectral analysis with fast fourier transform, two important parameters need to be determined: sampling rate f s Number of sum frequency domain samples N fft The sampling rate is determined according to the maximum frequency of the signal and the Nyquist sampling theorem, and the number of sampling points is determined according to the frequency resolution Vf.
Then
The operation rule of the FFT time decimation algorithm is summarized as follows:
the fast Fourier algorithm (FFT) mainly comprises two processes, namely a windowed time sequence f N (n) inversion of code position to obtainNew input sequenceThen toButterfly operation is performed to obtain a frequency spectrum F (k). As shown in fig. 3 (taking 8 points as an example), the fast fourier algorithm first performs code bit inversion on the sequences F (0), F (1), F (2), F (3), F (4), F (5), F (6), and F (7) to obtain new sequences F (0), F (1), F (4), F (2), F (1), F (5), F (3), and F (7), and then performs three-stage butterfly operation on the new sequences to obtain F (0), F (1), F (2), F (3), F (4), F (5), F (6), and F (7).
(1) Inversion of code bit
In order to arrange the output sequence F (k) in natural order, the input sequence F is arranged N (n) arranged in the order of code bit inversion, the algorithm is as follows:
let I be a sequential binary number, J be a reverse binary number, I = J =000B. Setting a constant N 2 Is half of a point number, N fft /2 as the midpoint of the sequence, set constant N 1 Count the total number of points N for the sequence fft -1 as sequence length. Firstly, judging that I is larger than or equal to J, if I = J =000B, performing reverse order, and skipping an index switching part. Set variable K = N 2 If K > J is judged to be =100B, J = J + K =000+100=100B, I = I +1=001B. And judging that I is larger than or equal to J, otherwise, entering an indexing exchange part: the values in the memory cells of I and J at this time are swapped. If J = J-K =100-100=000b, K = K/2=010b, and if K > J, J = J + K =000+010=010b, I = I +1=001+1=010b, and if I is less than the total number of sequential points N, J = J-K =100-100=000b, K is more than 010b, and I is more than 010b 1 ＝N fft And (5) returning to the condition that I is larger than or equal to J, and repeating the judgment and calculation until I reaches the total point number, and exiting the program.
(2) Butterfly operation
All operations are decomposed into M stages of operations, each stage of operation comprising N/2 basic butterflies. In the butterfly diagram, the butterfly operations that intersect each other in any one stage are called a group; in the butterfly diagramIn the sequence from top to bottom, the increment of the corresponding element sequence number between the upper and lower two groups is called the group interval. The L-th operation includes N/2 L The group interval of the L-th level is 2 L . W of each group in the same level N The distribution is the same, and the L-th level has 2 L-1 W N Multiplier, W in each group N The multipliers are distributed from top to bottom according to the following rule:
a first stage:
and a second stage:
and an L stage:L，
and an M stage:L，
the L-th basic butterfly operation input sequence interval is 2 L-1 The operation relationship between the input and the output is:
in the formula, A L (. O) represents an input element of the L-th basic butterfly, A L (. Cndot.) represents the output element of the L-th level basic butterfly.
To the vibration signal after the secondary de-noisingAnd carrying out windowed discrete Fourier transform to obtain a vibration signal frequency spectrum F (k) of the equipment, and taking the frequency spectrum F (k) as the basis of a next self-learning algorithm.
The Self-Learning algorithm (SLNM) for the Natural Mode of vibration of the device is shown in FIG. 4. The whole process of calculating the frequency spectrum of the vibration signal is introduced, in order to identify the natural vibration mode of the equipment, the vibration signal of the equipment in normal operation needs to be monitored for a long time, the vibration signal of the equipment in a reasonable time period is taken, denoising and windowing discrete Fourier transform are carried out on the vibration signal, the frequency spectrum of the equipment is calculated, the frequency spectrum in the time period is trained by using a self-learning algorithm, and finally the natural characteristic frequency and the amplitude of the equipment are obtained. 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) denoising, windowing, 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-th calculation period;
step four: to F t (k) Sequencing according to the amplitude values 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 vibration of the equipment and the corresponding amplitude (including the mean value and the variance); and plots each characteristic frequency component.
After the step of self-learning algorithm (SLNM) is completed, the natural frequency and corresponding amplitude of the equipment vibration number are calculated, and thus the natural vibration mode of the equipment is identified.
Example analysis and verification
The vibration sensor is arranged on the gear box of the wind driven generator, so that the vibration data of the gear box of the fan is measured in real time.
The self-learning time is set to be 4 hours, 1 second of vibration data is selected as samples every 5 minutes, 48 samples are counted, denoising and discrete Fourier transform are respectively carried out on the 48 samples, the frequency spectrum is obtained through calculation, then the self-learning (SLNM) algorithm is used for processing, the vibration characteristic parameters of the gearbox are obtained, and therefore the inherent vibration mode of the gearbox is identified. Table 1 shows the vibration frequency and amplitude of the gearbox of the wind driven generator.
TABLE 1
Note: in table 1, f1 to f6 respectively represent 6 vibration characteristic frequencies, and X1 to X6 respectively represent amplitudes corresponding to f1 to f 6; "/" indicates that the frequency components are absent from the spectrum.
As is apparent from Table 1, the frequency spectrum of the wind turbine gearbox is changing all the time and the change is large. In order to identify the natural vibration modes of the gearbox, the 48 frequency spectra were subjected to statistical analysis, and the characteristic parameters were calculated, the results of which 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, and the six characteristic frequency components are respectively: (622, 0.058), (1495, 0.067), (1213, 0.04), (2075, 0.04), (1825, 0.044), (3270, 0.024).
(1) The denoising of the vibration signal adopts a denoising algorithm combining a wavelet packet and singular value decomposition, and the denoising algorithm firstly carries out wavelet packet denoising on the vibration signal and then carries out singular value decomposition secondary denoising;
(2) In order to accurately calculate the frequency spectrum of the vibration signal and reduce the influence of frequency spectrum leakage, the invention firstly carries out windowing processing on the vibration signal and then calculates the frequency spectrum by utilizing fast Fourier transform;
(3) In order to calculate the characteristic parameters of the natural vibration of the device, the invention provides a self-learning (SLNM) algorithm.
The advantages of the invention are mainly divided into three aspects, summarized as follows:
the accuracy is as follows: an improved triangular window function is adopted when the frequency spectrum of the vibration signal is calculated, so that the influence of frequency spectrum leakage can be greatly reduced, and the frequency spectrum closer to a true value is obtained;
stability: the invention uses a denoising algorithm combining wavelet packet and singular value degradation, can effectively remove noise signals and is little interfered by the surrounding environment;
universality: the invention can be used for the natural vibration mode identification of various devices.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not intended to limit the present invention in any way, and all simple modifications, equivalent variations and modifications made to the above embodiments according to the technical spirit of the present invention are within the scope of the present 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 the content of the first and second substances,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 1 ,x 2 ,L,x n As the first vector y of the n-dimensional phase space 1 ；
(2) Move by one step to the right, decimate { x 2 ,x 3 ,L,x n+1 As the first vector y of the n-dimensional phase space 2 ；
(3) By analogy, a group of column vectors y is obtained 1 ,y 2 ,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 1 ,s 2 ,L,s p ,p＝min(m,n),s 1 ≥s 2 ≥L≥s p Is more than or equal to 0, wherein s 1 ,s 2 ,L,s p Is a matrix D m The obtained singular value s 1 ,s 2 ,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.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510024313.5A CN104807534B (en) | 2015-05-21 | 2015-05-21 | Equipment eigentone self study recognition methods based on on-line vibration data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510024313.5A CN104807534B (en) | 2015-05-21 | 2015-05-21 | Equipment eigentone self study recognition methods based on on-line vibration data |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104807534A CN104807534A (en) | 2015-07-29 |
| CN104807534B true CN104807534B (en) | 2018-04-13 |
Family
ID=53692565
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510024313.5A Expired - Fee Related CN104807534B (en) | 2015-05-21 | 2015-05-21 | Equipment eigentone self study recognition methods based on on-line vibration data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104807534B (en) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105203318A (en) * | 2015-10-23 | 2015-12-30 | 珠海格力电器股份有限公司 | Real-time monitoring method and device for running state of range hood |
| CN107515105B (en) * | 2016-06-15 | 2021-03-30 | 中国科学院沈阳自动化研究所 | Electromagnetic valve fault diagnosis method based on plunger vibration signal |
| CN106706119B (en) * | 2016-12-15 | 2019-05-03 | 北方工业大学 | Vibration source identification method and system based on signal frequency domain characteristics |
| CN107065912B (en) * | 2017-05-04 | 2020-08-11 | 厦门衡空科技有限公司 | Method and device for detecting landing of aircraft |
| CN110348491A (en) * | 2019-06-20 | 2019-10-18 | 燕山大学 | Rolling bearing fault recognition methods based on study dictionary and singular value decomposition |
| CN110376437B (en) * | 2019-07-18 | 2020-04-24 | 北京科技大学 | Order analysis method for overcoming non-order frequency component interference |
| CN111144362B (en) * | 2019-12-31 | 2023-07-25 | 上海数深智能科技有限公司 | Periodic optimization algorithm for vibration fault feature library of rotary equipment |
| CN112241599A (en) * | 2020-11-18 | 2021-01-19 | 中国联合网络通信集团有限公司 | Method and device for establishing fault analysis model |
| CN113155973B (en) * | 2021-05-05 | 2023-05-16 | 温州大学 | Beam damage identification method based on self-adaptive singular value decomposition |
| CN117033911B (en) * | 2023-10-07 | 2024-01-30 | 深圳市魔样科技有限公司 | Step counting analysis method based on intelligent glasses data |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5641905A (en) * | 1995-12-14 | 1997-06-24 | Quatro Corporation | Second derivative resonant ultrasound response analyzer |
| CN101852681A (en) * | 2010-03-31 | 2010-10-06 | 桂林电子科技大学 | Crack identification method of main shaft of boring machine |
| CN101858778A (en) * | 2010-05-28 | 2010-10-13 | 浙江大学 | Vibration monitoring-based wind generator set automatic fault diagnosis method |
| CN102520071B (en) * | 2011-12-20 | 2014-10-15 | 江苏方天电力技术有限公司 | Transmission tower modal parameter identification method based on improved subspace algorithm |
| CN103983412B (en) * | 2014-05-30 | 2017-06-06 | 北京航空航天大学 | For vibrating FEM updating avionic device operation mode measuring method |
| CN104165742B (en) * | 2014-07-17 | 2017-03-01 | 浙江工业大学 | A kind of operational modal analysis experimental technique based on mutual spectral function and device |
| CN104376159A (en) * | 2014-11-05 | 2015-02-25 | 汕头大学 | Large horizontal shaft wind turbine transmission chain and flexible design method thereof |
-
2015
- 2015-05-21 CN CN201510024313.5A patent/CN104807534B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN104807534A (en) | 2015-07-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104807534B (en) | Equipment eigentone self study recognition methods based on on-line vibration data | |
| CN109029977B (en) | Planetary gearbox early fault diagnosis method based on VMD-AMCKD | |
| CN109977920B (en) | Water turbine set fault diagnosis method based on time-frequency spectrogram and convolutional neural network | |
| CN105115594B (en) | Gear-box vibration signal fault signature extracting method based on Wavelet Entropy and information fusion | |
| CN106168942B (en) | A kind of fluctuation types dynamic data reconstructing method based on singular boundary method | |
| CN108801630B (en) | Gear fault diagnosis method for single-channel blind source separation | |
| Dong et al. | A repeated single-channel mechanical signal blind separation method based on morphological filtering and singular value decomposition | |
| CN105981025A (en) | Actually-measured marine environment data assimilation method based on sequence recursive filtering three-dimensional variation | |
| CN114782745B (en) | Ocean sound velocity profile classification method and device based on machine learning | |
| CN110782041B (en) | Structural modal parameter identification method based on machine learning | |
| CN114755010A (en) | Rotary machine vibration fault diagnosis method and system | |
| CN106548031A (en) | A kind of Identification of Modal Parameter | |
| CN113435304B (en) | Method, system, device and storage medium for extracting torsional vibration information of torsional vibration signal | |
| CN117200208B (en) | User-level short-term load prediction method and system based on multi-scale component feature learning | |
| Babu et al. | Fault diagnosis in bevel gearbox using coiflet wavelet and fault classification based on ANN including DNN | |
| CN111144230B (en) | Denoising method of time domain load signal based on VMD | |
| CN117473414A (en) | Bearing fault position identification method based on low-noise time-frequency image | |
| CN108267311A (en) | A kind of mechanical multidimensional big data processing method based on tensor resolution | |
| CN116992217A (en) | Mechanical signal noise reduction method based on multi-scale dynamic weighting multi-dimensional residual convolution | |
| CN107389367B (en) | A kind of signal adaptive solution enveloping method and calculating equipment based on optimal signal-to-noise ratio | |
| CN115291103A (en) | Motor fault diagnosis method based on GR-SWPT wavelet packet algorithm embedded with HD-RCF | |
| CN115166514A (en) | Motor fault identification method and system based on self-adaptive spectrum segmentation and denoising | |
| Chen et al. | Gear fault detection analysis method based on fractional wavelet transform and back propagation neural network | |
| RU2595929C2 (en) | Method and apparatus for compressing data depending on time signal | |
| CN111142134A (en) | Coordinate time series processing method and device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180413 Termination date: 20190521 |