CN104807534A - Equipment natural vibration mode self-learning recognition method based on online vibration data - Google Patents

Abstract

The invention discloses an equipment natural vibration mode self-learning recognition method based on online vibration data, which is characterized by comprising the steps of carrying out primary denoising on equipment vibration signals x(n) by using a wavelet packet transform algorithm so as to acquire denoised signals f1<^>(n); carrying out secondary denoising on the signals f1<^>(n) by using singular value decomposition so as to acquire secondarily denoised vibration signals f2<^>(n); carrying out spectral analysis on the denoised vibration signals f2<^>(n) by using windowed discrete Fourier algorithm, and calculating to acquire equipment vibration spectrum; training by using a self-learning algorithm to acquire the equipment vibration spectrum, and finally acquiring the natural characteristic frequency and the amplitude of equipment. The equipment natural vibration mode self-learning recognition method based on the online vibration data has the beneficial effect that the characteristic frequency of components can be recognized accurately when a research object of a complicated system.

Description

Based on the equipment eigentone self study recognition methods of on-line vibration data

Technical field

The invention belongs to on-line monitoring technique field, relate to the equipment eigentone self study recognition methods based on on-line vibration data.

Background technology

On-line monitoring technique is the most basic measures estimated equipment health status, and especially for plant equipment (or system), Vibration Monitoring, Analyzing is the most practical and effective method.The vibration signal of equipment and its physical construction closely related, in equipment running process, the situation of its internal mechanical structure can directly reflect effectively from vibration signal.By carrying out Treatment Analysis to equipment vibrating signal; the vibration mode that equipment under extraction normal operational condition is intrinsic; for distinguishing the vibration mode when exception or fault appear in equipment, to make a policy in time when abnormal vibrations model occurs, and formulate corresponding maintenance schedule.The present invention proposes a kind of vibration equipment natural mode automatic learning identification based on vibration data and abstracting method, pass through the analyzing and processing to vibration data and the training of one period of reasonable time, calculate the proper vibration parameter of equipment, it can be used as the health characteristics that equipment normally runs, thus provide a feasible evaluation index for the operation health status of equipment, guarantee the healthy operation and maintenance of equipment.

The eigentone (mode) of equipment refers to and comprises natural frequency, damping ratio and Mode Shape that equipment has by the natural vibration characteristic of physical construction.If understood the characteristic of plant machinery structure primary modal, then can foretell the actual vibration response under various vibration source effect of this structure.Thus, namely model analysis is analyzed these mode thus is obtained the important method that corresponding modal parameter is the fault diagnosis of structure dynamic design and equipment.Although the vibration Behavioral change of the equipment in reality (as aerogenerator, transformer, boats and ships, aircraft etc.) is different and complicated, model analysis provides the effective way studying various equipment practical structures characteristic.

Existing equipment eigentone recognition methods mainly computational modal analysis method, first utilizes 3 d modeling software, builds the high-precision three-dimensional solid model of equipment; Then, the three-dimensional model of equipment is imported in finite element analysis software, according to the actual conditions of equipment, the material properties parameter of set device, and mesh of finite element is divided into three-dimensional model, add displacement constraint, obtain the finite element model of equipment; Finally, model analysis is carried out to finite element model, obtain model natural frequency and mode.

Prior art is in finite element modal analysis, and particularly when research object is comparatively complicated system, the mode of oscillation calculated is often too in disorder, is difficult to the characteristic frequency of parts picking out real concern; Or when being made up of the vibration of certain rank many local modes, the numerical value of judging characteristic frequency can be very difficult to.

The present invention is without the need to setting up complicated model, and by the vibration data of direct collecting device different parts, process obtains modal parameter, the eigentone of identification equipment by analysis.

Summary of the invention

The object of the present invention is to provide the equipment eigentone self study recognition methods based on on-line vibration data, solve prior art in finite element modal analysis, particularly when research object is comparatively complicated system, the mode of oscillation calculated is often too in disorder, is difficult to pick out the characteristic frequency of parts of real concern, be difficult to the problem of the numerical value of judging characteristic frequency.

The technical solution adopted in the present invention is carried out according to following steps:

Step 1: establish the equipment vibrating signal directly collected by sensor to be x (n), its mathematical model is expressed as

x(n)＝f(n)+noise(n)

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

Step 2: utilize wavelet package transforms algorithm to carry out a denoising to equipment vibrating signal x (n), obtain the signal after denoising

Step 3: utilize svd to signal carry out second denoising, obtain the vibration signal after second denoising

Step 4: use windowing discrete fourier algorithm to the vibration signal after denoising carry out spectrum analysis, calculate vibration equipment frequency spectrum;

Step 5: utilize self-learning algorithm to train and obtain vibration equipment frequency spectrum, finally obtain inherent feature frequency and the amplitude of equipment.

Further, in described step 2, wavelet package transforms algorithm to the process that equipment vibrating signal x (n) carries out a denoising is:

1. the WAVELET PACKET DECOMPOSITION of signal; Select a wavelet basis, and determine the number of plies N of a WAVELET PACKET DECOMPOSITION, then carry out N layer WAVELET PACKET DECOMPOSITION to signal, decomposition formula is:

$\left\{\begin{array}{c}{C}_{j+\mathrm{1,2}l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)h(n-2k)\\ C_{j+\mathrm{1,2}l+1}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)g(n-2k)\end{array}\right.$

In formula: j<N, C _j,ln () is the n-th WAVELET PACKET DECOMPOSITION coefficient of jth layer l node in decomposition tree, k and n all refers to certain node which WAVELET PACKET DECOMPOSITION coefficient concrete;

Get vibration signal x (n) for signal to be analyzed, carry out the WAVELET PACKET DECOMPOSITION of signal, signal decomposition to different time frequency space;

2. Optimum Wavelet Packet is determined: according to principle of minimum cost, to a given entropy criterion calculation Optimum Wavelet Packet;

Entropy is defined as:

M＝-∑P _jllgP _jl

Wherein, and work as P _jlwhen=0, P _jllgP _jl=0;

3. to the WAVELET PACKET DECOMPOSITION coefficient C under each decomposition scale _j,lselect a threshold values to process, adopt soft threshold method;

Soft threshold values noise-eliminating method is defined as:

$\stackrel{\&OverBar;}{C_{j,l}}=\left\{\begin{array}{cc}{C}_{j,l}-\mathrm{\&lambda;},& C_{j,l}\&GreaterEqual;\mathrm{\&lambda;}\\ 0,& \left|C_{j,l}\right|<\mathrm{\&lambda;}\\ C_{j,l}+\mathrm{\&lambda;},& C_{j,l}\&le;-\mathrm{\&lambda;}\end{array}\right.$

In formula, λ is threshold values, C _j,lfor wavelet packet coefficient, for the wavelet packet coefficient after process;

4. wavelet package reconstruction, carry out wavelet reconstruction according to the WAVELET PACKET DECOMPOSITION low frequency coefficient of n-th layer and quantification treatment coefficient, reconstruction formula is as follows:

$C_{j,l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l}\left(n\right)h(n-2k)+\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l+1}\left(n\right)g(n-2k)$

In formula, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l the node of jth in decomposition tree, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l+1 the node of jth in decomposition tree, h (n-2k) and g (n-2k) filter coefficient for defining in multiresolution analysis, use Wavelet Packet Algorithm to carry out denoising to equipment vibrating signal x (n), obtain the vibration signal after a denoising

Further, utilize svd to signal in described step 3 carry out second denoising concrete steps as follows:

(1) from subsequence { x is extracted in time series ₁, x ₂, L, x _nfirst vectorial y of phase space is tieed up as n ₁;

(2) move right a step-length, extracts { x ₂, x ₃, L, x _n+1first vectorial y of phase space is tieed up as n ₂;

(3) one group of column vector { y can by that analogy, be obtained ₁, y ₂, L, y _m;

(4) point in each vectorial corresponding phase space reconstruction, institute's directed quantity forms the matrix D of m × n dimension _m:

$D_m=\left[\begin{array}{cccc}{x}_1& x_2& L& x_n\\ x_2& x_3& L& x_{n+1}\\ L& L& M& L\\ x_m& x_{m+1}& L& x_{m+n+1}\end{array}\right]$

D _mfor the attractor track matrix that Embedded dimensions m, time delay are 1, if the vibration signal measured contains certain noise, then D _m=D+W, wherein D, W represent smooth signal and D corresponding to noise signal respectively _min track matrix, to matrix D _mmake svd, D _m=USV', U and V are respectively m × n and n × n rank matrix, and UU'=E, VV'=E, E are unit matrix, and S is m × n diagonal matrix, diagonal element s ₁, s ₂, L, s _p, p=min (m, n), s ₁>=s ₂>=L>=s _p>=0, wherein s ₁, s ₂, L, s _pfor matrix D _msingular value, by the singular value s obtained ₁, s ₂, L, s _pfront k (k≤p) item, other zero setting, obtain new diagonal matrix S', recycling SVD decompose inverse process D _m'=US'V' obtains new matrix D _m', matrix D _m' be D _mthe best approach matrix, according to the process of phase space reconfiguration, by D _m' obtain the vibration signal after second denoising

Further, in described step 4, adopt triangular self-convolution window to carry out windowing discrete fourier algorithm and calculate rumble spectrum.

Further, in described step 5, self-learning algorithm step is:

Step one: the time of getting computation period is 1 second, determine that training T.T. is T computation period, get t=1, t represents calculating sequence number;

Step 2: the vibration data x inputting equipment in t computation period _t(n);

Step 3: to x _tn () carries out denoising, windowing, discrete Fourier transformation (FFT), calculate frequency spectrum F _t(k), F _tfrequency spectrum in (k) indication equipment t computation period;

Step 4: to F _tk () is sorted from big to small by amplitude, obtain B _t(k);

Step 5: to B _tk (), setting threshold value ρ, if when amplitude is greater than threshold value ρ, thinks that frequency content corresponding to this amplitude is characteristic frequency, by it stored in eigenmatrix C _tin; If be less than threshold value ρ, cast out;

Step 6: t=t+1, enters next computation period, judges t≤T, if so, then skips to step 2; Otherwise skip to step 7;

Step 7: six steps of above-mentioned algorithm obtain the Faults by Vibrating set C of each monitoring point in the training period, the amplitude of computing equipment vibration natural frequency and correspondence thereof; And draw the curve of each characteristic frequency composition, obtain natural frequency and the respective magnitudes thereof of the equipment number of shaking, thus identify the eigentone of equipment.

The invention has the beneficial effects as follows when research object is comparatively complicated system, characteristic frequency of parts can be gone out by accurate recognition.

Accompanying drawing explanation

Fig. 1 is wavelet packet of the present invention 3 layers of decomposition tree structural representation;

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

Fig. 3 is the Fast Fourier Transform (FFT) frame diagram of the present invention's eight points;

Fig. 4 is self study of the present invention (SLNM) algorithm flow chart.

Embodiment

Below in conjunction with embodiment, the present invention is described in detail.The equipment eigentone self study recognition methods that the present invention is based on on-line vibration data is carried out according to following steps:

Step 1: equipment working site noise background is complicated, the vibration signal gathered is very easy to the pollution being subject to these noises, indeed vibrations data in signal are submerged under powerful ground unrest, thus bring very large difficulty to the calculating of equipment natural frequency.If the equipment vibrating signal directly collected by sensor is x (n), its mathematical model can be expressed as:

x(n)＝f(n)+noise(n)

Wherein, f (n) is original signal; Noise (n) is noise signal.

Step 2: utilize wavelet package transforms algorithm to carry out a denoising to equipment vibrating signal x (n), obtain the signal after denoising the denoise algorithm that the present invention adopts wavelet package transforms to combine with svd carries out denoising to vibration signal,

Wavelet package transforms: for given orthogonal scaling function and wavelet function φ (t) of correspondence, there is Double-scaling equation

In formula: { h (n) } and { g (n) } filter coefficient for defining in multiresolution analysis, for scaling function, φ (t) is wavelet function.

Note μ ₁t ()=φ (t), defines μ by recurrence relation _m(t) be:

${\mathrm{\&mu;}}_{2l}\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\&Sigma;}}h\left(n\right){\mathrm{\&mu;}}_l(2t-n),$

${\mathrm{\&mu;}}_{2l+1}\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\&Sigma;}}g\left(n\right){\mathrm{\&mu;}}_l(2t-n).$

Claim race function { μ _l(t) | l=0,1,2, L} are relative to orthogonal scaling function orthogonal Wavelet Packet, wherein l is positive integer.

That extracts from wavelet packet can form L ²(R) one group of orthogonal basis is just called L ²(R) a wavelet packet basis, that extracts from wavelet library can form L ²(R) one group of orthonormal basis is L ²(R) wavelet packet basis.WAVELET PACKET DECOMPOSITION for a signal can adopt multiple wavelet packet basis to realize, different wavelet packet basiss has different decomposition result to signal, its result can reflect that the degree of signal is also different, therefore seeking the wavelet basis of one group of optimum, is can the vital task of a most effective expression signal.

The concrete steps of wavelet packet denoising:

1. the WAVELET PACKET DECOMPOSITION of signal.Select a wavelet basis (i.e. wavelet function), and determine the number of plies N (General N is 3 or 4) of a WAVELET PACKET DECOMPOSITION, then carry out N layer WAVELET PACKET DECOMPOSITION to signal, decomposition formula is:

$\left\{\begin{array}{c}{C}_{j+\mathrm{1,2}l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)h(n-2k)\\ C_{j+\mathrm{1,2}l+1}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)g(n-2k)\end{array}\right.$

In formula: j<N, C _j,ln () is the n-th WAVELET PACKET DECOMPOSITION coefficient of jth layer l node in decomposition tree.K and n all refers to certain node which WAVELET PACKET DECOMPOSITION coefficient concrete.

Get vibration signal x (n) for signal to be analyzed, carry out the WAVELET PACKET DECOMPOSITION of signal, signal decomposition to different time frequency space.

Such as 3 layers of WAVELET PACKET DECOMPOSITION are carried out to signal.Decomposition texture is (WAVELET PACKET DECOMPOSITION tree is the one method for expressing comparatively intuitively of WAVELET PACKET DECOMPOSITION) as shown in Figure 1.In Fig. 1, (0,0) node represents original signal, (1,0) node represents the 1st layer of low frequency component of WAVELET PACKET DECOMPOSITION, (2,0) and (3,0) node represents the 2nd layer respectively, and the component of the 3rd layer of the 0th node, other by that analogy.

2. determine Optimum Wavelet Packet: according to principle of minimum cost, to a given entropy standard (adopting Shannon entropy), calculate Optimum Wavelet Packet;

Utilize wavelet packet analysis, a good wavelet packet basis must be selected to be used for describing signal.In order to select a good wavelet packet basis, the cost function of a first given sequence, the base making cost function minimum is found in all wavelet packet basiss, for a Setting signal, the base making Shannon entropy (cost function) minimum is exactly the wavelet packet basis the most effectively representing this signal, and this base is just called best base.

Shannon entropy (cost function) is defined as:

M＝-∑P _jllgP _jl

Wherein, and work as P _jlwhen=0, P _jllgP _jl=0.

The present invention the 2. step carries out on the basis 1. walked, first the 1. step repeatedly decompose, calculate the coefficient of dissociation C of signal under different wavelet basis _j,l, second step utilizes 1. to walk the coefficient of dissociation C obtained _j,lcalculate shonnon entropy, find minimum entropy, the wavelet basis of its correspondence is Optimal wavelet bases.

3. to the WAVELET PACKET DECOMPOSITION coefficient C under each decomposition scale _j,lselect a suitable threshold values to process, the present invention adopts soft threshold method.

Soft threshold values noise-eliminating method is defined as:

$\stackrel{\&OverBar;}{C_{j,l}}=\left\{\begin{array}{cc}{C}_{j,l}-\mathrm{\&lambda;},& C_{j,l}\&GreaterEqual;\mathrm{\&lambda;}\\ 0,& \left|C_{j,l}\right|<\mathrm{\&lambda;}\\ C_{j,l}+\mathrm{\&lambda;},& C_{j,l}\&le;-\mathrm{\&lambda;}\end{array}\right.$

In formula, λ is threshold values, C _j,lfor wavelet packet coefficient, for the wavelet packet coefficient after process.

4. wavelet package reconstruction, carries out wavelet reconstruction according to the WAVELET PACKET DECOMPOSITION low frequency coefficient of n-th layer and quantification treatment coefficient.Reconstruction formula is as follows:

$C_{j,l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l}\left(n\right)h(n-2k)+\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l+1}\left(n\right)g(n-2k)$

In formula, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l the node of jth in decomposition tree, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l+1 the node of jth in decomposition tree, h (n-2k) and g (n-2k) filter coefficient for defining in multiresolution analysis.

Use Wavelet Packet Algorithm to carry out denoising to equipment vibrating signal x (n), obtain the vibration signal after a denoising ${\hat{f}}_1\left(n\right)[x_1,x_2,L,x_N].$

Step 3: svd (Single Value Decomposition, SVD) is right carry out second denoising:

After wavelet packet denoising is carried out to the vibration signal of equipment, obtain the vibration signal after a denoising the present invention uses svd pair the method of secondary noise reduction carries out noise reduction based on Smooth Systems signal, the Different Effects of random noise signal to continuous wavelet transform track singular values of a matrix, and the method cancelling noise signal, does not almost affect system signal.First vibration signal is reconstructed in phase space, then continuous wavelet transform can reflect the dynamics of system, for the track matrix characterizing attractor carries out svd, utilize the characteristic of singular spectrum to reduce the noise contribution in vibration signal, good effect can be reached.

Vibration signal be a discrete-time series, carry out phase space reconfiguration to it, its concrete steps are as follows:

(1) from subsequence { x is extracted in time series ₁, x ₂, L, x _nfirst vectorial y of phase space is tieed up as n ₁;

(2) move right a step-length, extracts { x ₂, x ₃, L, x _n+1first vectorial y of phase space is tieed up as n ₂;

(3) one group of column vector { y can by that analogy, be obtained ₁, y ₂, L, y _m;

(4) point in each vectorial corresponding phase space reconstruction, institute's directed quantity forms the matrix D of m × n dimension _m:

$D_m=\left[\begin{array}{cccc}{x}_1& x_2& L& x_n\\ x_2& x_3& L& x_{n+1}\\ L& L& M& L\\ x_m& x_{m+1}& L& x_{m+n+1}\end{array}\right]$

D _mfor Embedded dimensions m, time delay is the attractor track matrix of 1.If the vibration signal measured contains certain noise, then D _m=D+W, wherein D, W represent smooth signal and D corresponding to noise signal respectively _min track matrix.To matrix D _mmake svd, D _m=USV', U and V are respectively m × n and n × n rank matrix, and UU'=E, VV'=E (E is unit matrix).S is m × n diagonal matrix, diagonal element s ₁, s ₂, L, s _p, p=min (m, n), s ₁>=s ₂>=L>=s _p>=0, wherein s ₁, s ₂, L, s _pmatrix D _msingular value.In svd noise reduction process, choosing of noise reduction order is very crucial, selects the difference of order that noise reduction can be caused obviously different.When selecting order too high, can make still to contain more noise information in de-noising signal, good noise reduction cannot be reached; When selected order is too low, the imperfect of the later information of noise reduction can be caused, even cause the distortion of waveform.Here the singular value s will obtained ₁, s ₂, L, s _pfront k (k≤p) item, other zero setting, obtain new diagonal matrix S'.The inverse process D that recycling SVD decomposes _m'=US'V' obtains new matrix D _m', matrix D _m' be D _mthe best approach matrix, according to the process of phase space reconfiguration, by D _m' obtain the vibration signal after denoising

Complete svd pair after carrying out denoising, obtain the vibration signal after second denoising

Step 4: windowing discrete fourier algorithm calculates rumble spectrum:

The present invention uses windowing discrete fourier algorithm to the vibration signal after denoising carry out spectrum analysis, calculate vibration equipment frequency spectrum.

Spectrum analysis is the frequency-region signal of horizontal ordinate being that the time-domain signal of horizontal ordinate become by Fourier transform with frequency with time, and then a kind of analytical approach of the phase place of trying to achieve about former frequency of time domain signal composition and amplitude.The figure drawn by frequency and phase place or amplitude is called phase spectrum and amplitude spectrum, and comparatively conventional in engineering is amplitude spectrum.In order to improve the efficiency calculating frequency spectrum, the present invention adopts Fast Fourier Transform (FFT) (FFT) to process vibration signal.

The calculating formula of discrete Fourier transformation (DFT):

$F\left(k\right)=\mathrm{DFT}\left[{\hat{f}}_2\left(n\right)\right]=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\hat{f}}_2\left(n\right)e^{-j\frac{2\mathrm{\&pi;}}{N}\mathrm{kn}},(k=0~N-1)$

In formula, the time-domain signal of to be data length be n; F (k) is for transforming to the result of complex field, identical with time-domain signal length.

Windowing principle: because discrete Fourier transformation is to finite length sequence definition, therefore, to Infinite Sequences when calculating frequency spectrum F (k), must block or segmentation x (n), this is equivalent to handle be 1 with amplitude, length is the rectangle sequence w of N _nn ()=u (n)-u (n-N) is multiplied, the N point sequence x after blocking _n(n) be:

$f_N\left(n\right)={\hat{f}}_2\left(n\right)w_N\left(n\right)$

This amplitude is 1, and length is the rectangle sequence w of N _nn () is exactly rectangular window function, this to signal block or segmentation is exactly windowing, the function selected during windowing process is different, and the impact brought to spectrum analysis is also different.

Spectrum leakage: so-called spectrum leakage refers to that the signal energy of a certain frequency when carrying out discrete Fourier transformation is diffused into the phenomenon of side frequency point.Discrete Fourier transformation is carried out to measured signal, signal time-domain windowed is equivalent to frequency domain convolution, and this blocking causes frequency spectrum to occur error, and its effect makes frequency spectrum centered by actual frequency values, spread to both sides with the shape of window function spectrum waveform, produce " leakage effect ".Leakage effect can increase new frequency content, and spectrum size is changed.From energy point of view, the energy that frequency leakage phenomenon is equivalent to the various frequency content places of original signal penetrates in other frequency contents, so also known as Power leakage.

Window function is chosen: spectrum leakage be discrete Fourier transformation intrinsic, the secondary lobe that it and window function are composed is closely related, if make the height of secondary lobe go to zero, thus makes energy be relatively concentrated in main lobe, just can obtain comparatively close to the frequency spectrum of actual value.Consider the feature of vibration signal, the present invention adopts triangular self-convolution window.

Lower Fig. 2 gives the normalization amplitude characteristic (figure dot-dashed line is quarter window, and solid line is triangular self-convolution window) of quarter window and triangular self-convolution window.

As can be seen from Figure 2, the secondary lobe 26dB lower than main lobe that quarter window is maximum, and the secondary lobe 52dB lower than main lobe that triangular self-convolution window is maximum, and side lobe attenuation speed obviously accelerates, contrast can find that the energy of triangular self-convolution window mainly concentrates in main lobe, secondary lobe reduces greatly, and triangular self-convolution window can suppress the impact of spectrum leakage effectively.

FFT calculates frequency spectrum: when carrying out spectrum analysis with Fast Fourier Transform (FFT), needs to determine two important parameters: sample rate f _swith frequency domain sample points N _fft, sampling rate, according to the highest frequency of signal, is determined according to nyquist sampling theorem, and sampling number is determined according to frequency resolution Vf.

$\frac{f_s}{N}\&le;\mathrm{Vf},$ Then $N_{\mathrm{fft}}\&GreaterEqual;\frac{f_s}{\mathrm{Vf}}$

The operation law of FFT decimation in time algorithm is summarized as follows:

Fast Fourier algorithm (FFT) mainly comprises two processes, first will add the time series f of window _nn () carries out Binary Code Inversion, obtain new list entries then right carry out butterfly computation, obtain frequency spectrum F (k).The calculating process of fast Fourier algorithm as shown in Figure 3 (for 8 points), first to sequence f (0), f (1), f (2), f (3), f (4), f (5), f (6), f (7) carries out Binary Code Inversion, obtain new sequence f (0), f (1), f (4), f (2), f (1), f (5), f (3), f (7), again three grades of butterfly computations are carried out to new sequence, obtain F (0), F (1), F (2), F (3), F (4), F (5), F (6), F (7).

(1) Binary Code Inversion

In order to make output sequence F (k) by natural order arrangement, by list entries f _nn (), by the arrangement of Binary Code Inversion order, algorithm is as follows:

If I is order binary number, J is inverted order binary number, I=J=000B.If constant N ₂for the half of counting, N _fft/ 2 as the mid point of sequence, if constant N ₁for the total points N of sequence _fft-1, as sequence length.First judge I >=J, I=J=000B meets then to be needed to carry out inverted order, skips index switching part.If variable K=N ₂=100B, judges K > J, is then J=J+K=000+100=100B, I=I+1=001B.Judge I >=J again, otherwise enter index switching part: the value in the storage unit of I and J is now exchanged.Judge K > J now again, otherwise J=J-K=100-100=000B, K=K/2=010B, judge K > J now again, be then J=J+K=000+010=010B, I=I+1=001+1=010B, judge that I is less than the total points N of sequence ₁=N _fft-1, then turn back to I>=J, by that analogy, cycle criterion always and calculating, until I reaches backed off after random program of always counting.

(2) butterfly computation

Whole computing is broken down into the computing of M level and carries out, and every grade of computing all comprises N/2 basic butterfly computation.In butterfly diagram, the butterfly computation mutually intersected in any one-level is called group; By order from top to bottom in butterfly diagram, between upper and lower two groups, the increment of corresponding element sequence number is called group interval.The computing of L level comprises N/2 ^lindividual group, the group of L level is spaced apart 2 ^l.With the W of each group in one-level _ndistribute identical, in L level, have 2 ^l-1individual W _nmultiplier, W in each group _nmultiplier is from top to bottom by following rule distribution:

The first order:

The second level:

L level: l,

M level: l,

The basic butterfly computation list entries of L level is spaced apart 2 ^l-1, between its constrained input, operation relation is:

$\left\{\begin{array}{c}{A}_{L+1}\left(p\right)=A_L\left(p\right)+W_N^r{A}_L\left({p+2}^{L-1}\right)\\ A_{L+1}(p+2^{L-1})=A_L\left(p\right)-W_N^r{A}_L\left({p+2}^{L-1}\right)\end{array}\right.$

In formula, A _l() represents the input element of the basic butterfly computation of L level, A _l() represents the output element of the basic butterfly computation of L level.

To the vibration signal after second denoising carry out windowing DFT algorithm, obtain vibration signals spectrograph F (k) of equipment, as the basis of next step self-learning algorithm.

The self-learning algorithm (Self-Learning Natural Mode, SLNM) of vibration equipment natural mode as shown in Figure 4.Forgoing describe the whole process calculating vibration signals spectrograph, in order to the eigentone of identification equipment, vibration signal when needing long-time monitoring equipment normally to run, get one section of equipment vibrating signal within reasonable time, denoising and windowing DFT algorithm are carried out to vibration signal, the frequency spectrum of computing equipment, utilizes self-learning algorithm to train frequency spectrum in this time period, finally obtains inherent feature frequency and the amplitude of equipment.Self-learning algorithm step is as follows:

Step one: the time of getting computation period is 1 second, determine that training T.T. is T computation period, get t=1, t represents calculating sequence number;

Step 2: the vibration data x inputting equipment in t computation period _t(n);

Step 3: to x _tn () carries out denoising, windowing, discrete Fourier transformation (FFT), calculate frequency spectrum F _t(k), F _tfrequency spectrum in (k) indication equipment t computation period;

Step 4: to F _tk () is sorted from big to small by amplitude, obtain B _t(k);

Step 5: to B _tk (), setting threshold value ρ, if when amplitude is greater than threshold value ρ, thinks that frequency content corresponding to this amplitude is characteristic frequency, by it stored in eigenmatrix C _tin; If be less than threshold value ρ, cast out;

Step 6: t=t+1, enters next computation period, judges t≤T, if so, then skips to step 2; Otherwise skip to step 7;

Step 7: six steps of above-mentioned algorithm obtain the Faults by Vibrating set C of each monitoring point in the training period, the amplitude (comprising: average, variance) of computing equipment vibration natural frequency and correspondence thereof; And draw the curve of each characteristic frequency composition.

After completing the step of self-learning algorithm (SLNM), calculate natural frequency and the respective magnitudes thereof of the equipment number of shaking, thus identified the eigentone of equipment.

Sample calculation analysis is verified

By installing vibration transducer at wind-driven generator wheel-box, measure the vibration data of gear case of blower in real time, the present invention, to this data analysis process, calculates the Faults by Vibrating of wind-driven generator wheel-box.

The time of self study is set as 4 hours, the vibration data of 1 second is chosen as sample every 5 minutes, amount to 48 samples, respectively denoising, discrete Fourier transformation are carried out to 48 samples, calculate frequency spectrum, recycling self study (SLNM) algorithm processes, and obtains the Faults by Vibrating of gear case, thus identifies the eigentone of gear case.Table 1 is wind-driven generator wheel-box vibration frequency and amplitude thereof.

Table 1

Note: what in table 1, f1 ~ f6 represented respectively is 6 vibration performance frequencies, and what X1 ~ X6 represented respectively is the amplitude that f1 ~ f6 is corresponding; "/" to represent in frequency spectrum not this frequency content.

Can obviously find out from table 1, the frequency spectrum of wind-driven generator wheel-box is changing always, and amplitude of variation is larger.In order to identify the eigentone of gear case, carrying out statistical study to 48 frequency spectrums, calculating characteristic parameter, result is as shown in table 2 is the Faults by Vibrating of wind-driven generator wheel-box.

Table 2

Table 2 shows the Faults by Vibrating of gear case, six characteristic frequency compositions are respectively: (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 adopted is the denoise algorithm that wavelet packet combines with svd, first carries out wavelet packet denoising to vibration signal, then carries out svd secondary noise reduction;

(2) in order to calculate vibration signals spectrograph exactly, reduce the impact of spectrum leakage, the present invention first carries out windowing process to vibration signal, and recycling Fast Fourier Transform (FFT) calculates frequency spectrum;

(3) in order to the characteristic parameter of computing equipment proper vibration, the present invention proposes a kind of self study (SLNM) algorithm.

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

Accuracy: what adopt when calculating vibration signals spectrograph is the quarter window function improved, and can affect with greatly reducing spectrum leakage, obtain more close to the frequency spectrum of actual value;

Stability: the present invention uses wavelet packet and singular value to degrade the denoise algorithm combined, and effectively can remove noise signal, the interference by surrounding environment is little;

Versatility: the present invention can be used for the eigentone identification of plurality of devices.

The above is only to better embodiment of the present invention, not any pro forma restriction is done to the present invention, every any simple modification done above embodiment according to technical spirit of the present invention, equivalent variations and modification, all belong in the scope of technical solution of the present invention.

Claims (5)

1., based on the equipment eigentone self study recognition methods of on-line vibration data, it is characterized in that, carry out according to following steps:

Step 1: establish the equipment vibrating signal directly collected by sensor to be x (n), its mathematical model is expressed as

x(n)＝f(n)+noise(n)

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

Step 2: utilize wavelet package transforms algorithm to carry out a denoising to equipment vibrating signal x (n), obtain the signal after denoising

Step 3: utilize svd to signal carry out second denoising, obtain the vibration signal after second denoising

Step 4: use windowing discrete fourier algorithm to the vibration signal after denoising carry out spectrum analysis, calculate vibration equipment frequency spectrum;

Step 5: utilize self-learning algorithm to train and obtain vibration equipment frequency spectrum, finally obtain inherent feature frequency and the amplitude of equipment.

2., according to the equipment eigentone self study recognition methods based on on-line vibration data described in claim 1, it is characterized in that: in described step 2, wavelet package transforms algorithm to the process that equipment vibrating signal x (n) carries out a denoising is:

1. the WAVELET PACKET DECOMPOSITION of signal; Select a wavelet basis, and determine the number of plies N of a WAVELET PACKET DECOMPOSITION, then carry out N layer WAVELET PACKET DECOMPOSITION to signal, decomposition formula is:

$\left\{\begin{array}{c}{C}_{j+\mathrm{1,2}l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)h(n-2k)\\ C_{j+\mathrm{1,2}l+1}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j,1}\left(n\right)g(n-2k)\end{array}\right.$

In formula: j<N, C _j,ln () is the n-th WAVELET PACKET DECOMPOSITION coefficient of jth layer l node in decomposition tree, k and n all refers to certain node which WAVELET PACKET DECOMPOSITION coefficient concrete;

Get vibration signal x (n) for signal to be analyzed, carry out the WAVELET PACKET DECOMPOSITION of signal, signal decomposition to different time frequency space;

2. Optimum Wavelet Packet is determined: according to principle of minimum cost, to a given entropy criterion calculation Optimum Wavelet Packet;

Entropy is defined as:

M＝-ΣP _jllgP _jl

Wherein, and work as P _jlwhen=0, P _jllgP _jl=0;

3. to the WAVELET PACKET DECOMPOSITION coefficient C under each decomposition scale _j,lselect a threshold values to process, adopt soft threshold method;

Soft threshold values noise-eliminating method is defined as:

$\stackrel{\&OverBar;}{C_{j,l}}=\left\{\begin{array}{cc}{C}_{j,l}-\mathrm{\&lambda;},& C_{j,l}\&GreaterEqual;\mathrm{\&lambda;}\\ 0,& \left|C_{j,l}\right|<\mathrm{\&lambda;}\\ C_{j,l}+\mathrm{\&lambda;},& C_{j,l}\&le;-\mathrm{\&lambda;}\end{array}\right.$

In formula, λ is threshold values, C _j,lfor wavelet packet coefficient, for the wavelet packet coefficient after process;

4. wavelet package reconstruction, carry out wavelet reconstruction according to the WAVELET PACKET DECOMPOSITION low frequency coefficient of n-th layer and quantification treatment coefficient, reconstruction formula is as follows:

$C_{j,l}\left(k\right)=\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l}\left(n\right)h(n-2k)+\underset{n}{\mathrm{\&Sigma;}}{C}_{j+\mathrm{1,2}l+1}\left(n\right)g(n-2k)$

In formula, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l the node of jth in decomposition tree, C _j+1,2ln () is the WAVELET PACKET DECOMPOSITION system of+1 layer of 2l+1 the node of jth in decomposition tree, h (n-2k) and g (n-2k) filter coefficient for defining in multiresolution analysis, use Wavelet Packet Algorithm to carry out denoising to equipment vibrating signal x (n), obtain the vibration signal after a denoising ${\hat{f}}_1\left(n\right)=[x_1,x_2,L,x_N].$

3., according to the equipment eigentone self study recognition methods based on on-line vibration data described in claim 1, it is characterized in that:

Utilize svd to signal in described step 3 carry out second denoising concrete steps as follows:

(1) from subsequence { x is extracted in time series ₁, x ₂, L, x _nfirst vectorial y of phase space is tieed up as n ₁;

(2) move right a step-length, extracts { x ₂, x ₃, L, x _n+1first vectorial y of phase space is tieed up as n ₂;

(3) one group of column vector { y by that analogy, is obtained ₁, y ₂, L, y _m;

(4) point in each vectorial corresponding phase space reconstruction, institute's directed quantity forms the matrix D of m × n dimension _m:

$D_m=\left[\begin{array}{cccc}{x}_1& x_2& L& x_n\\ x_2& x_3& L& x_{n+1}\\ L& L& M& L\\ x_m& x_{m+1}& L& x_{m+n+1}\end{array}\right]$

D _mfor the attractor track matrix that Embedded dimensions m, time delay are 1, if the vibration signal measured contains certain noise, then D _m=D+W, wherein D, W represent smooth signal and D corresponding to noise signal respectively _min track matrix, to matrix D _mmake svd, D _m=USV', U and V are respectively m × n and n × n rank matrix, and UU'=E, VV'=E, E are unit matrix, and S is m × n diagonal matrix, diagonal element s ₁, s ₂, L, s _p, p=min (m, n), s ₁>=s ₂>=L>=s _p>=0, wherein s ₁, s ₂, L, s _pfor matrix D _msingular value, by the singular value s obtained ₁, s ₂, L, s _pfront k (k≤p) item, other zero setting, obtain new diagonal matrix S', recycling SVD decompose inverse process D _m'=US'V' obtains new matrix D _m', matrix D _m' be D _mthe best approach matrix, according to the process of phase space reconfiguration, by D _m' obtain the vibration signal after second denoising ${\hat{f}}_2\left(n\right).$

4. according to the equipment eigentone self study recognition methods based on on-line vibration data described in claim 1, it is characterized in that: in described step 4, adopt triangular self-convolution window to carry out windowing discrete fourier algorithm and calculate rumble spectrum.

5., according to the equipment eigentone self study recognition methods based on on-line vibration data described in claim 1, it is characterized in that: in described step 5, self-learning algorithm step is:

Step one: the time of getting computation period is 1 second, determine that training T.T. is T computation period, get t=1, t represents calculating sequence number;

Step 2: the vibration data x inputting equipment in t computation period _t(n);

Step 3: to x _tn () carries out denoising, windowing, discrete Fourier transformation (FFT), calculate frequency spectrum F _t(k), F _tfrequency spectrum in (k) indication equipment t computation period;

Step 4: to F _tk () is sorted from big to small by amplitude, obtain B _t(k);

Step 5: to B _tk (), setting threshold value ρ, if when amplitude is greater than threshold value ρ, thinks that frequency content corresponding to this amplitude is characteristic frequency, by it stored in eigenmatrix C _tin; If be less than threshold value ρ, cast out;

Step 6: t=t+1, enters next computation period, judges t≤T, if so, then skips to step 2; Otherwise skip to step 7;

Step 7: six steps of above-mentioned algorithm obtain the Faults by Vibrating set C of each monitoring point in the training period, the amplitude of computing equipment vibration natural frequency and correspondence thereof; And draw the curve of each characteristic frequency composition, obtain natural frequency and the respective magnitudes thereof of the equipment number of shaking, thus identify the eigentone of equipment.