CN112697265A - Self-adaptive multi-fractal method and equipment working condition monitoring device - Google Patents
Self-adaptive multi-fractal method and equipment working condition monitoring device Download PDFInfo
- Publication number
- CN112697265A CN112697265A CN202011240325.9A CN202011240325A CN112697265A CN 112697265 A CN112697265 A CN 112697265A CN 202011240325 A CN202011240325 A CN 202011240325A CN 112697265 A CN112697265 A CN 112697265A
- Authority
- CN
- China
- Prior art keywords
- data
- signal
- fractal
- curve
- algorithm
- 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.)
- Withdrawn
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
- G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
- G06F2218/16—Classification; Matching by matching signal segments
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention discloses a self-adaptive multi-fractal method and a device for monitoring the working condition of equipment, which decompose a vibration signal by using an EMD algorithm, remove a noise component and a trend term by using a nonlinear discrimination algorithm, reserve a fractal signal component, interpolate extreme points by using a cubic spline interpolation function, fit the envelope of each fractal signal component by using a least square method, separate a frequency modulation part, estimate an instantaneous frequency by using a TEO algorithm and calculate a corresponding instantaneous scale, automatically determine a vibration signal detrending result according to the analysis scale, calculate a multi-fractal spectrum of a detrended signal, extract singular indexes corresponding to a left endpoint, a right endpoint and the extreme points of the multi-fractal spectrum as characteristic parameters of the running state of the equipment to identify the running state of the equipment, and finally deploy the algorithm to a device for monitoring the running state of the equipment to accurately distinguish the running state of the equipment, the equipment state monitoring device has good flexibility and portability, and is convenient for engineering application.
Description
Technical Field
The invention relates to the field of equipment state monitoring and fault diagnosis, in particular to a self-adaptive multi-fractal method and an equipment working condition monitoring device.
Background
The device vibration signal contains rich fractal features that can describe the operating state of the device. The box dimension, power spectrum analysis and re-standard range method can estimate the single-fractal parameters of stationary signals, and the de-trend fluctuation analysis (DFA) can estimate the single-fractal dimension of non-stationary signals. However, when the device fails, the vibration signal is usually non-stationary and has a multi-fractal characteristic, and the conventional fractal dimension estimation method generates a relatively large error. The multi-fractal detrending fluctuation analysis (MFDF) can estimate multi-fractal parameters of non-stationary signals, but the MFDF method has the problems that the analysis scale needs to be manually determined, the fitting polynomial trend order is difficult to determine, and the data segment is discontinuous. Currently, there is a document that proposes an MFDFA version (MFDFAemd) based on EMD to solve the problem of MFDFA. However, the linear filtering method adopted by mfdfame is easy to destroy the fractal structure of the original signal, and there is a negative frequency phenomenon, and these defects seriously affect the application effect of mfdfame. In summary, in the prior art, it is difficult to accurately extract the multi-fractal features of the device vibration signal, and it is difficult to accurately detect the device operating state.
Disclosure of Invention
Aiming at the defects, the invention provides a self-adaptive multi-fractal method and a device for monitoring the working condition of equipment (the method provided by the invention is abbreviated as MFDFoomed). The method provided by the invention is adopted to analyze the equipment vibration signal, can effectively extract the multi-fractal characteristics of the equipment vibration signal, overcomes the problems that the analysis scale of the MFDF method needs to be manually determined, the fitting polynomial trend order is difficult to determine and the data section is discontinuous, solves the phenomena of original signal fractal structure damage and negative frequency existing in the MFDF method, and has the advantages of high accuracy and precision of analysis results, high accuracy of equipment operation state identification results and the like.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows: an adaptive multi-fractal method is characterized by comprising the following steps:
step 1: measuring a device vibration signal x (k) by using an acceleration sensor at a sampling frequency fs, wherein k =1,2, …, N and N are lengths of the sampling signal;
step 2: the signal x (k) is decomposed into the sum of n components and a trend term, i.e., using an Empirical Mode Decomposition (EMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the EMD algorithmn(k) Represents a trend term derived from the EMD algorithm;
and step 3: eliminating noise components and trend terms from EMD decomposition results by adopting a nonlinear discrimination algorithm, and reserving components c containing fractal featuresf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
m times, j =1,2, …, m, until
To obtain cf(k) Is/are as followsFrequency modulation part FMm(k)，ej(k) Represents cj(k) Envelope of cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
Step 6: calculating FM by using Teager Energy Operator (TEO)m(k) Obtaining the instantaneous frequency of cf(k) Instantaneous frequency instf off(k)，
To obtain cf(k) Instantaneous scale of
；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
；
And 8: will be provided withY s (k) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
and step 9: calculate variance of each piece of data:
step 10: calculating a q-th order function:
step 11: changing the value of s, s = sfF =1,2, …, p, repeating the above steps 3 to 10, resulting in a square with respect to q and sDifference function Fq(s);
Step 12: if it is notx(k) Presence of fractal features, thenF q (s) And sizesThere is a power law relationship between:F q (s)~s H q()h (q) represents the generalized Hurst index of x (k);
when in useqWhen the value is not less than 0, the reaction time is not less than 0,H(0) determined by the logarithmic averaging procedure defined by:
step 13: calculating a standard scale index τ (q) = qH (q) -1 for the signal x (k);
step 14: calculating the singular index α and the multifractal spectrum f (α) of the signal x (k):
α=H(q)+q H’(q),
f (alpha) = q (alpha-H (q)) +1, wherein H’(q) represents the first derivative of h (q);
step 15: and extracting singular indexes corresponding to a left end point, a right end point and an extreme point of the multi-fractal spectrum f (alpha), and describing the running state of the equipment by using the 3 parameters.
Further, the Empirical Mode Decomposition (EMD) algorithm in step 2 comprises the following steps:
1) the first screening process: finding out upper and lower local extreme points of data x (k), fitting the upper and lower local extreme points by cubic spline curve to obtain local maximum envelope and local minimum envelope of signal x (k), averaging values of corresponding points of the two envelopes to obtain an average curve m1；
Then, the original signal x (k) and the average curve m are obtained1A difference of (i) thath 10=x(k)-m 1Ending the first screening process;
2) the second screening process:h 10re-regarded as the original signal, and the above step 1) is repeated to obtain the signalh 11= h 10-m 11This isMean parameterm 11Representsh 10Is repeated j times until the mean value of the curve is 0.2<SD<0.3 the screening process is stopped, here
At this time, the process of the present invention,h j1= h j1(-1)-m j1in this case, it can be considered thath j1Is an Intrinsic Mode Function (IMF), and the 1 st IMF is defined asc 1=h j1；
3) Subtracting from the original signalc 1Is obtained byr 1=x(k)-c 1Then will ber 1When new data is taken, the two steps are repeated, and the 2 nd IMF can be obtained;
4) repeating the operation of step 3) to obtain a series of IMFs ifr n Having become a monotonic curve, the screening process stops, and the original signal is finally decomposed into the following form:
。
further, the step 3 nonlinear discrimination algorithm includes the following steps:
1) performing rearrangement operation and substitution operation on signals c (k), and using c as data obtained by rearrangement operationshuf(k) Indicating that data obtained after the substitution operation are csurr(k) Represents;
2) for c (k), cshuf(k) And csurr(k) Performing multi-fractal Detrended Fluctuation Analysis (MFDF) respectively to obtain a generalized Hurst index curve, wherein the generalized Hurst index curve of c (k) is represented by H (q); c. Cshuf(k) Generalized Hurst exponential curve of (1) using Hshuf(q) represents; c. Csurr(k) Generalized Hurst exponential curve of (1) using Hsurr(q) represents;
3) two parameters e are defined1And e2，
Further, the data rearrangement operation in the step 1) comprises the following steps: randomly randomizing the order of the components c (k).
Further, the data replacement operation in the step 1) comprises the following steps:
1) performing a discrete fourier transform on component c (k) to obtain the phase of component c (k);
2) replacing the original phase of component c (k) with a set of pseudo-independent identically distributed numbers located in the (-pi, pi) interval;
3) performing inverse discrete Fourier transform on the frequency domain data subjected to phase substitution to obtain data cIFFT(k) To obtain data cIFFT(k) The real part of (a).
Further, the MFDFA method in step 2) includes the following steps:
1) contour of construction x (k), k =1,2, …, NY(i)：
x (k) represents c (k) or c in step 2) of claim 3shuf(k) Or csurr（k）；
2) Signal profileY(i) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
3) fitting a polynomial trend of each section of data by using a least square method, and then calculating the variance of each section of data:
y v (i) Is a first of fittingvTrend of the segment data, if the fitted polynomial trend ismOrder, then note the de-trending process as (MF-) DFAm；
4) Calculate the firstqAverage of the order fluctuation function:
5) if it is notx(k) Presence of self-similar features thenqMean value of order fluctuation functionF q (s) And time scalesThere is a power law relationship between:
F q (s)~s H q()；
when in useqIf =0, the formula in step 4) diverges, at which timeH(0) Determined by the logarithmic averaging procedure defined by:
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
Further, the least square method in the step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) A sequence or pair c generated by interpolating the local maxima off(k) The local minimum value of (a) is interpolated to generate a sequence, n represents the length of the interpolated sequence;
1) a set of functions r is selected in advancek(t)，k=1,2,…,m，m<n, constructioning function
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t) in which rk(t) represents a second order polynomial, a third order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve or a complex function curve;
2) calculating a least squares metric
；
3) Let J pair akPartial derivative of
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Wherein R isTRepresenting the transposed matrix of R, R-1An inverse matrix representing R;
4）rk(t) respectively selecting a second-order polynomial, a third-order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve and a composite function curve for calculation, then comparing least square indexes J generated by various curve forms, and selecting a curve form r corresponding to the minimum J as the curve formk(t) form (a).
Further, the Teager energy operator method in the step 6 comprises the following steps:
1) for signal c (k), k =1,2, …, N, c (k) = FMm(k) Building a function
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2) Let d (k) = c (k) -c (k-1), instantaneous frequency instf (k) of signal c (k) be defined as:
based on the above adaptive multi-fractal method, the device for monitoring the working condition of the equipment implementing the method includes the following parts in step 16: the system comprises a data line, an acceleration sensor, a data acquisition card, a case, a notebook computer and signal analysis software, wherein the acceleration sensor is connected with the data acquisition card through the data line, the data acquisition card is installed in the case, the case is connected with the notebook computer through the data line, the signal analysis software is installed on the notebook computer, and the signal analysis software is used for the algorithm.
The method comprises the following steps:
step 1), the following steps: collecting vibration signals by using an acceleration sensor;
step 2), the step of: decomposing an original signal into different component sum forms, wherein some components correspond to noise and trend terms, and some components contain fractal features;
and 3, step 3: removing noise components and trend items in the signal decomposition result by using a nonlinear discrimination algorithm, and only reserving signal components containing fractal features;
4) to 6): separating the frequency modulation part of each fractal signal component, and estimating the instantaneous frequency and the instantaneous scale of each fractal signal component by using TEO;
step 7), the steps of: selecting proper fractal signal components according to the analysis scale, and summing the selected fractal signal components to be used as a signal de-trend result corresponding to the analysis scale;
8) to 14): performing fluctuation analysis on the signal detrending result corresponding to each analysis scale to obtain a multi-fractal spectrum of the original signal;
step 15): extracting singular indexes corresponding to a left end point, a right end point and an extreme point of the multi-fractal spectrum, and taking the three parameters as characteristic parameters of the running state of the equipment;
16) step: the algorithm is deployed on an equipment state monitoring device to monitor the equipment state.
By adopting the technical scheme, compared with the prior art, the invention has the following advantages:
1) the vibration signal is adaptively decomposed by adopting an EMD method, the noise component and the trend term are removed according to a nonlinear filtering method, the fractal structure of the original signal can be protected, and the damage of the linear filtering method to the fractal structure of the original signal is avoided;
2) the frequency modulation part of the separated signal component estimates the instantaneous frequency of the signal component by using TEO, so that the instantaneous frequency can be ensured to keep a positive value, and the negative frequency phenomenon is avoided;
3) calculating corresponding instantaneous scale according to the instantaneous frequency of the signal component, and performing fluctuation analysis according to the instantaneous scale of the signal component, so that the defect of manually setting the scale is avoided;
4) the EMD method is used for automatically determining the type of the signal trend, the continuity of the signal trend is guaranteed, and the defects of the prior art are effectively overcome.
5) The accuracy and precision of the analysis result are high, and the accuracy of the identification result of the running state of the equipment is high.
The invention is further illustrated with reference to the following figures and examples.
Drawings
FIG. 1 is a flow chart of the method of the present invention in an embodiment of the present invention;
FIG. 2 is a schematic diagram of an apparatus state monitoring device according to an embodiment of the present invention;
fig. 3 is a multi-fractal simulation signal generated by a multi-fractal cascade model in the embodiment of the present invention;
fig. 4 is an instantaneous frequency of a multi-fractal simulation signal obtained by using the MFDFAemd method in the embodiment of the present invention, where the number of signal components is 10;
fig. 5 is an instantaneous frequency of a multi-fractal simulation signal obtained by using the mfdfoomd method in the embodiment of the present invention, where the number of signal components is 10;
FIG. 6 is a comparison graph of the multi-fractal simulation signal analysis results respectively using MFDF, MFDFEMd and MFDFEMd methods in the embodiment of the present invention;
FIG. 7 shows the calculation results of two non-linear discriminant parameters in the embodiment of the present invention, where the symbols "circle" and "square" represent e1And e2；
FIG. 8 is a comparison diagram of analysis results of noisy multi-fractal simulation signals respectively using MFDF, MFDFame, and MFDFAEmd methods in the embodiment of the present invention;
FIG. 9 is a diagram illustrating correlation coefficients between each signal component and an original signal obtained by EMD according to an embodiment of the present invention;
FIG. 10 is a comparison diagram of analysis results of noisy multi-fractal simulation signals respectively using MFDFA, MFDFAEmd based on correlation filtering, and MFDFAEmd based on correlation filtering in the embodiment of the present invention;
FIG. 11 shows four vibration signals of the gearbox in the embodiment of the invention, wherein (a) - (d) respectively represent normal, light scratch, heavy scratch and broken tooth gear states;
FIG. 12 is a multi-fractal spectrum of the four gearbox vibration signals obtained using MFDFA in an embodiment of the present invention;
FIG. 13 is a multi-fractal spectrum of the four gearbox vibration signals obtained by using MFDFame in the embodiment of the present invention;
FIG. 14 is a multi-fractal spectrum of the four gearbox vibration signals obtained by using MFDFoomd in the embodiment of the present invention;
fig. 15 is a classification result of singular indexes corresponding to left, right, and extreme points of a multi-fractal spectrum obtained by MFDFA on the four gear box states in the embodiment of the present invention, where the "circle", "square", "plus", and "diamond" symbols represent normal, light, heavy, and broken gear states, respectively;
fig. 16 is a classification result of singular indexes corresponding to left end point, right end point and extreme point of a multi-fractal spectrum obtained by MFDFAemd on the states of the four gearboxes in the embodiment of the present invention, where symbols "circle", "square", "plus" and "diamond" represent normal, light scratch, heavy scratch and broken tooth gear states, respectively;
fig. 17 shows the classification results of the singular indexes corresponding to the left end point, the right end point, and the extreme point of the multi-fractal spectrum obtained from MFDFAoemd according to the embodiment of the present invention on the states of the four gearboxes, where the symbols "circle", "square", "plus", and "diamond" represent the normal, light, heavy, and broken gear states, respectively.
Detailed Description
An embodiment, as shown in fig. 1 and fig. 2, is a device for monitoring operating conditions of a self-adaptive multi-fractal method and equipment, including the following steps:
step 1: measuring a device vibration signal x (k) by using an acceleration sensor at a sampling frequency fs, wherein k =1,2, …, N and N are lengths of the sampling signal;
step 2: the signal x (k) is decomposed into the sum of n components and a trend term, i.e., using an Empirical Mode Decomposition (EMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the EMD algorithmn(k) Represents a trend term derived from the EMD algorithm;
and step 3: eliminating noise components and trend terms from EMD decomposition results by adopting a nonlinear discrimination algorithm, and reserving components c containing fractal featuresf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executeFormula (II)
m times, j =1,2, …, m, until
To obtain cf(k) Frequency modulation part FMm(k)，ej(k) Represents cj(k) Envelope of cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
Step 6: calculating FM by using Teager Energy Operator (TEO)m(k) Obtaining the instantaneous frequency of cf(k) Instantaneous frequency instf off(k)，
To obtain cf(k) Instantaneous scale of
；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
；
And 8: will be provided withY s (k) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
and step 9: calculate variance of each piece of data:
step 10: calculating a q-th order function:
step 11: changing the value of s, s = sfF =1,2, …, p, repeating the above steps 3 to 10, resulting in a variance function F about q and sq(s);
Step 12: if it is notx(k) Presence of fractal features, thenF q (s) And sizesThere is a power law relationship between:F q (s)~s H q()h (q) represents the generalized Hurst index of x (k);
when in useqWhen the value is not less than 0, the reaction time is not less than 0,H(0) determined by the logarithmic averaging procedure defined by:
step 13: calculating a standard scale index τ (q) = qH (q) -1 for the signal x (k);
step 14: calculating the singular index α and the multifractal spectrum f (α) of the signal x (k):
α=H(q)+q H’(q),
f (alpha) = q (alpha-H (q)) +1, wherein H’(q) represents the first derivative of h (q);
step 15: extracting singular indexes corresponding to a left end point, a right end point and an extreme point of the multi-fractal spectrum f (alpha), wherein an Empirical Mode Decomposition (EMD) algorithm in the step 2 comprises the following steps:
1) the first screening process: finding out upper and lower local extreme points of data x (k), fitting the upper and lower local extreme points by cubic spline curve to obtain local maximum envelope and local minimum envelope of signal x (k), averaging values of corresponding points of the two envelopes to obtain an average curve m1；
Then, the original signal x (k) and the average curve m are obtained1A difference of (i) thath 10=x(k)-m 1Ending the first screening process;
2) the second screening process:h 10re-regarded as the original signal, and the above step 1) is repeated to obtain the signalh 11= h 10-m 11Here parameterm 11Representsh 10Is repeated j times until the mean value of the curve is 0.2<SD<0.3 the screening process is stopped, here
At this time, the process of the present invention,h j1= h j1(-1)-m j1in this case, it can be considered thath j1Is an Intrinsic Mode Function (IMF), and the 1 st IMF is defined asc 1=h j1；
3) Subtracting from the original signalc 1Is obtained byr 1=x(k)-c 1Then will ber 1When new data is taken, the two steps are repeated, and the 2 nd IMF can be obtained;
4) repeating the operation of step 3) to obtain a series of IMFs ifr n Having become a monotonic curve, the screening process stops, and the original signal is finally decomposed into the following form:
。
the step 3 of the nonlinear discriminant algorithm comprises the following steps:
1) performing rearrangement operation and substitution operation on signals c (k), and using c as data obtained by rearrangement operationshuf(k) Indicating that data obtained after the substitution operation are csurr(k) Represents;
2) for c (k), cshuf(k) And csurr(k) Performing multi-fractal Detrended Fluctuation Analysis (MFDF) respectively to obtain a generalized Hurst index curve, wherein the generalized Hurst index curve of c (k) is represented by H (q); c. Cshuf(k) Generalized Hurst exponential curve ofBy Hshuf(q) represents; c. Csurr(k) Generalized Hurst exponential curve of (1) using Hsurr(q) represents;
3) two parameters e are defined1And e2，
The data rearrangement operation in the step 1) comprises the following steps: randomly randomizing the order of the components c (k).
The data replacement operation in the step 1) comprises the following steps:
1) performing a discrete fourier transform on component c (k) to obtain the phase of component c (k);
2) replacing the original phase of component c (k) with a set of pseudo-independent identically distributed numbers located in the (-pi, pi) interval;
3) performing inverse discrete Fourier transform on the frequency domain data subjected to phase substitution to obtain data cIFFT(k) To obtain data cIFFT(k) The real part of (a).
The MFDF method in the step 2) comprises the following steps:
1) contour of construction x (k), k =1,2, …, NY(i)：
x (k) represents c (k) or c in step 2) of claim 3shuf(k) Or csurr（k）；
2) Signal profileY(i) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
3) fitting a polynomial trend of each section of data by using a least square method, and then calculating the variance of each section of data:
y v (i) Is a first of fittingvTrend of the segment data, if the fitted polynomial trend ismOrder, then note the de-trending process as (MF-) DFAm；
4) Calculate the firstqAverage of the order fluctuation function:
5) if it is notx(k) Presence of self-similar features thenqMean value of order fluctuation functionF q (s) And time scalesThere is a power law relationship between:
F q (s)~s H q()；
when in useqIf =0, the formula in step 4) diverges, at which timeH(0) Determined by the logarithmic averaging procedure defined by:
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
The least square method in the step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) A sequence or pair c generated by interpolating the local maxima off(k) The local minimum value of (a) is interpolated to generate a sequence, n represents the length of the interpolated sequence;
1) a set of functions r is selected in advancek(t)，k=1,2,…,m，m<n, constructioning function
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t) in which rk(t) represents a second order polynomial, a third order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve or a complex function curve;
2) calculating a least squares metric
；
3) Let J pair akPartial derivative of
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Wherein R isTRepresenting the transposed matrix of R, R-1An inverse matrix representing R;
4）rk(t) respectively selecting a second-order polynomial, a third-order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve and a composite function curve for calculation, then comparing least square indexes J generated by various curve forms, and selecting a curve form r corresponding to the minimum J as the curve formk(t) form (a).
The Teager energy operator method in the step 6 comprises the following steps:
1) for signal c (k), k =1,2, …, N, c (k) = FMm(k) Building a function
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2) Let d (k) = c (k) -c (k-1), instantaneous frequency instf (k) of signal c (k) be defined as:
Firstly, a multi-fractal cascade model is adopted
The generated multi-fractal simulation signal verifies the performance of the MFDFA, the MFDFAemd and the MFDFAemd. In this example, p =0.375 and n =14, the resulting multi-fractal simulation signal is shown in fig. 3. The instantaneous frequency of the multi-fractal simulation signal is calculated by using the MFDFAemd, and the result is shown in fig. 4. As can be seen from fig. 4, the instantaneous frequency calculated by mfdfame has many negative frequencies, and the mfdfame analysis result has a large error because the negative frequencies have no physical significance. The instantaneous frequency of the multi-fractal simulation signal is calculated by using the MFDFAoemd, and the result is shown in fig. 5. As can be seen from fig. 5, the instantaneous frequencies calculated by MFDFAoemd are all positive frequencies, so the MFDFAoemd analysis result is in line with the actual situation. Next, a multi-fractal spectrum of the multi-fractal simulation signal is calculated using MFDFA, MFDFAemd, and MFDFAoemd, respectively, and the result is shown in fig. 6. According to the results shown in fig. 6, it can be found through calculation that the average of the absolute errors between the fractal spectrum obtained from MFDFA and the theoretical value is 0.064, the average of the relative errors is 12.21%, the average of the absolute errors between the fractal spectrum obtained from MFDFAemd and the theoretical value is 0.035, the average of the relative errors is 6.57%, and the average of the absolute errors between the fractal spectrum obtained from MFDFAemd and the theoretical value is flatThe mean value is 0.028, and the mean value of relative error is 5.18%, so that the multi-fractal spectrum obtained by MFDFAEmd is reduced by 56.25% and the mean value of relative error is reduced by 57.58% compared with the mean value of absolute error of multi-fractal spectrum obtained by MFDFA, the multi-fractal spectrum obtained by MFDFAEmd is reduced by 20.00% and the mean value of relative error is reduced by 21.16%. Fig. 7 shows the calculation results of two non-linear discrimination parameters in the embodiment of the present invention, and it can be seen that all signal components include fractal signal components. Then, a noisy signal with a signal-to-noise ratio of 20dB is constructed by a method of adding white Gaussian noise to the multi-fractal simulation signal. The multi-fractal spectrum of the noise-containing multi-fractal simulation signal was calculated by using MFDFA, MFDFAemd, and MFDFAemd, respectively, and the result is shown in fig. 8. According to the results shown in fig. 8, the multi-fractal spectrum obtained from mfdfaamd completely deviates from the theoretical value, the mean absolute error value of the multi-fractal spectrum obtained from MFDFAemd and the theoretical value is 0.084, the mean relative error value is 11.81%, the mean absolute error value of the multi-fractal spectrum obtained from MFDFAemd and the theoretical value is 0.033, and the mean relative error value is 5.38%, so that the multi-fractal spectrum obtained from MFDFAemd is reduced by 60.71% and the mean relative error value is reduced by 54.45% compared with the mean absolute error value of the multi-fractal spectrum obtained from MFDFAemd. It follows that MFDFAoemd has better noise immunity than MFDFA and mfdfaeemd. Fig. 9 shows the correlation coefficient between each signal component and the original signal obtained by EMD in the embodiment of the present invention, and it can be seen that the 7 th component has the weakest correlation with the original signal and should be removed from the original signal. Fig. 10 is a comparison diagram of analysis results of a noisy multi-fractal simulation signal by using MFDFA, mfdfame based on correlation filtering, and mfdfame based on correlation filtering, respectively, in the embodiment of the present invention. As can be seen from fig. 10, the analysis result of the noise-containing multi-fractal simulation signal by the mfdfame based on the correlation filtering and the mfdfame based on the correlation filtering completely deviates from the theoretical value, so that the fractal structure of the original signal is easily damaged by the correlation filtering method.
The invention relates to a gearbox fault simulation experiment table for gearbox vibration data. This experiment simulates the process of a tooth from normal to failure by making different degrees of scoring on the root of a tooth until the tooth is finally completely destroyed. The gear box used in the experiment is in two-stage gear transmission, the number of gear teeth from the input end to the output end is respectively 25, 40, 22 and 55, the fault gear teeth are positioned on the input shaft gear, the rotating speed of the driving motor is 2000RMP, and the vibration signal of the gear box is measured by an acceleration sensor positioned on the shell of the input end. The collected vibration signals comprise four fault states of normal, light scratch, heavy scratch and broken tooth, and the vibration signals represent the process of the gear tooth from normal to failure to a certain extent. The vibration signal sampling frequency was 16384Hz, and 20 pieces of data with a length of 10000 points were collected in each gearbox state, and these four gearbox vibration signals are shown in FIG. 11. Firstly, the four gearbox vibration signals are analyzed by using an MFDF method, and the multi-fractal spectrums corresponding to the four gearbox vibration signals are obtained as shown in FIG. 12, so that the multi-fractal spectrums corresponding to the light scratches and the heavy scratches are seriously overlapped. Then, the four gearbox vibration signals are analyzed by using an MFDFame method, and the multi-fractal spectrums corresponding to the four gearbox vibration signals are obtained as shown in FIG. 13, so that the multi-fractal spectrums corresponding to the four gearbox states are seriously overlapped. Finally, the four gearbox vibration signals are analyzed by using an MFDFoomd method, and the multi-fractal spectrums corresponding to the four gearbox vibration signals are obtained as shown in fig. 14, so that the multi-fractal spectrums of the vibration signals in the tooth breaking state are obviously different from the multi-fractal spectrums corresponding to the states of the other three gearboxes, and the multi-fractal spectrums corresponding to the normal, light-scratch and heavy-scratch gear states can be clearly separated when alpha is less than 0.4. Singular indexes corresponding to a left end point, a right end point and an extreme point of a multi-fractal spectrum obtained by the MFDF, MFDFAemd and MFDFAemd methods are respectively extracted to classify the four states of the gearbox, and the results are respectively shown in FIGS. 15-17. As can be seen from fig. 15, the normal state and the tooth breakage state can be correctly distinguished by using the singular indexes corresponding to the left end point, the right end point, and the extreme point of the multi-fractal spectrum obtained by the MFDFA method, but the light scratch state and the heavy scratch state cannot be distinguished, so that the gear box state recognition rate is 50%. As can be seen from fig. 16, the normal state and the tooth breakage state can be correctly distinguished by using the singular indexes corresponding to the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained by the MFDFAemd method, but the light scratch state and the heavy scratch state cannot be distinguished, so that the gear box state recognition rate is 50%. As can be seen from fig. 17, the four gear box states can be correctly distinguished by using the singular indexes corresponding to the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained by the MFDFAoemd method, so that the gear box state identification rate is 100%. It can be seen that the accuracy of the state identification of the gearbox can be improved by 50% by adopting the MFDFoomed method.
From the test results, it was assumed after analysis that:
1) the vibration signal is adaptively decomposed by adopting an EMD method, the noise component and the trend term are removed according to a nonlinear filtering method, the fractal structure of the original signal can be protected, and the damage of the linear filtering method to the fractal structure of the original signal is avoided;
2) the frequency modulation part of the separated signal component estimates the instantaneous frequency of the signal component by using TEO, so that the instantaneous frequency can be ensured to keep a positive value, and the negative frequency phenomenon is avoided;
3) calculating corresponding instantaneous scale according to the instantaneous frequency of the signal component, and performing fluctuation analysis according to the instantaneous scale of the signal component, so that the defect of manually setting the scale is avoided;
4) the EMD method is used for automatically determining the type of the signal trend, the continuity of the signal trend is ensured, and the defects of the prior art are effectively overcome;
5) the accuracy and precision of the analysis result are high, and the accuracy of the identification result of the running state of the equipment is high.
It should be appreciated by those skilled in the art that the foregoing embodiments are merely exemplary for better understanding of the present invention, and should not be construed as limiting the scope of the present invention as long as the modifications are made according to the technical solution of the present invention.
Claims (9)
1. An adaptive multi-fractal method is characterized by comprising the following steps:
step 1: measuring a device vibration signal x (k) by using an acceleration sensor at a sampling frequency fs, wherein k =1,2, …, N and N are lengths of the sampling signal;
step 2: the signal x (k) is decomposed into the sum of n components and a trend term, i.e., using an Empirical Mode Decomposition (EMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the EMD algorithmn(k) Represents a trend term derived from the EMD algorithm;
and step 3: eliminating noise components and trend terms from EMD decomposition results by adopting a nonlinear discrimination algorithm, and reserving components c containing fractal featuresf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
m times, j =1,2, …, m, until
To obtain cf(k) Frequency modulation part FMm(k)，ej(k) Represents cj(k) Envelope of cj(k)=FM(j-1)(k)，c1(k)= cf（k）；
Step 6: calculating FM by using Teager Energy Operator (TEO)m(k) Obtaining the instantaneous frequency of cf(k) Instantaneous frequency instf off(k)，
To obtain cf(k) Instantaneous scale of
；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
；
And 8: will be provided withY s (k) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
and step 9: calculate variance of each piece of data:
step 10: calculating a q-th order function:
step 11: changing the value of s, s = sfF =1,2, …, p, repeating the above steps 3 to 10, resulting in a variance function F about q and sq(s);
Step 12: if it is notx(k) Presence of fractal features, thenF q (s) And sizesThere is a power law relationship between:F q (s)~s H q()h (q) represents the generalized Hurst index of x (k);
when in useqWhen the value is not less than 0, the reaction time is not less than 0,H(0) determined by the logarithmic averaging procedure defined by:
step 13: calculating a standard scale index τ (q) = qH (q) -1 for the signal x (k);
step 14: calculating the singular index α and the multifractal spectrum f (α) of the signal x (k):
α=H(q)+q H’(q),
f (alpha) = q (alpha-H (q)) +1, wherein H’(q) represents the first derivative of h (q);
step 15: and extracting singular indexes corresponding to a left end point, a right end point and an extreme point of the multi-fractal spectrum f (alpha), and describing the running state of the equipment by using the 3 parameters.
2. The adaptive multi-fractal method according to claim 1, wherein the Empirical Mode Decomposition (EMD) algorithm in step 2 comprises the following steps:
1) the first screening process: finding out upper and lower local extreme points of data x (k), fitting the upper and lower local extreme points by cubic spline curve to obtain local maximum envelope and local minimum envelope of signal x (k), averaging values of corresponding points of the two envelopes to obtain an average curve m1；
Then, the original signal x (k) and the average curve m are obtained1A difference of (i) thath 10=x(k)-m 1Ending the first screening process;
2) the second screening process:h 10is re-regarded as the original signal and,repeating the step 1) to obtain the producth 11= h 10-m 11Here parameterm 11Representsh 10Is repeated j times until the mean value of the curve is 0.2<SD<0.3 the screening process is stopped, here
At this time, the process of the present invention,h j1= h j1(-1)-m j1in this case, it can be considered thath j1Is an Intrinsic Mode Function (IMF), and the 1 st IMF is defined asc 1=h j1；
3) Subtracting from the original signalc 1Is obtained byr 1=x(k)-c 1Then will ber 1When new data is taken, the two steps are repeated, and the 2 nd IMF can be obtained;
4) repeating the operation of step 3) to obtain a series of IMFs ifr n Having become a monotonic curve, the screening process stops, and the original signal is finally decomposed into the following form:
。
3. the adaptive multi-fractal method according to claim 1, wherein the step 3 non-linear discriminant algorithm includes the steps of:
1) performing rearrangement operation and substitution operation on signals c (k), and using c as data obtained by rearrangement operationshuf(k) Indicating that data obtained after the substitution operation are csurr(k) Represents;
2) for c (k), cshuf(k) And csurr(k) Performing multi-fractal Detrended Fluctuation Analysis (MFDF) to obtain generalized Hurst index curve, and using H (q) tableShown in the specification; c. Cshuf(k) Generalized Hurst exponential curve of (1) using Hshuf(q) represents; c. Csurr(k) Generalized Hurst exponential curve of (1) using Hsurr(q) represents;
3) two parameters e are defined1And e2，
4. The adaptive multi-fractal method as claimed in claim 3, wherein the data reordering operation in step 1) comprises the following steps: randomly randomizing the order of the components c (k).
5. The adaptive multi-fractal method of claim 3, characterized in that: the data replacement operation in the step 1) comprises the following steps:
1) performing a discrete fourier transform on component c (k) to obtain the phase of component c (k);
2) replacing the original phase of component c (k) with a set of pseudo-independent identically distributed numbers located in the (-pi, pi) interval;
3) performing inverse discrete Fourier transform on the frequency domain data subjected to phase substitution to obtain data cIFFT(k) To obtain data cIFFT(k) The real part of (a).
6. The adaptive multi-fractal method of claim 3, characterized in that: the MFDF method in the step 2) comprises the following steps:
1) contour of construction x (k), k =1,2, …, NY(i)：
x (k) represents c (k) or c in step 2) of claim 3shuf(k) Or csurr（k）；
2) Signal profileY(i) Divided into non-overlappingN s Length of segment beingsDue to data lengthNUsually cannot be removedsSo that a segment of data remains unavailable; in order to fully utilize the length of the data, the data is segmented with the same length from the reverse direction of the data, so that 2 is obtainedN s Segment data;
3) fitting a polynomial trend of each section of data by using a least square method, and then calculating the variance of each section of data:
y v (i) Is a first of fittingvTrend of the segment data, if the fitted polynomial trend ismOrder, then note the de-trending process as (MF-) DFAm；
4) Calculate the firstqAverage of the order fluctuation function:
5) if it is notx(k) Presence of self-similar features thenqMean value of order fluctuation functionF q (s) And time scalesThere is a power law relationship between:
F q (s)~s H q()；
when in useqIf =0, the formula in step 4) diverges, at which timeH(0) Determined by the logarithmic averaging procedure defined by:
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
7. The adaptive multi-fractal method as claimed in claim 1, wherein the least square method in step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) A sequence or pair c generated by interpolating the local maxima off(k) The local minimum value of (a) is interpolated to generate a sequence, n represents the length of the interpolated sequence;
1) a set of functions r is selected in advancek(t)，k=1,2,…,m，m<n, constructioning function
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t) in which rk(t) represents a second order polynomial, a third order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve or a complex function curve;
2) calculating a least squares metric
；
3) Let J pair akPartial derivative of
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Wherein R isTRepresenting the transposed matrix of R, R-1An inverse matrix representing R;
4）rk(t) respectively selecting a second-order polynomial, a third-order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve and a composite function curve for calculation, then comparing least square indexes J generated by various curve forms, and selecting a curve form r corresponding to the minimum J as the curve formk(t) form (a).
8. The adaptive multi-fractal method of claim 1, characterized in that: the Teager energy operator method in the step 6 comprises the following steps:
1) for signal c (k), k =1,2, …, N, c (k) = FMm(k) Building a function
ψ(c(k))=c2(k)- c(k+1) c(k-1);
2) Let d (k) = c (k) -c (k-1), instantaneous frequency instf (k) of signal c (k) be defined as:
9. an apparatus condition monitoring device for implementing the adaptive multi-fractal method as claimed in any of claims 1-8, characterized in that: the algorithm comprises a data line, an acceleration sensor, a data acquisition card, a case, a notebook computer and signal analysis software, wherein the acceleration sensor is connected with the data acquisition card through the data line, the data acquisition card is installed in the case, the case is connected with the notebook computer through the data line, the signal analysis software is installed on the notebook computer, and the signal analysis software is used for realizing the algorithm.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011240325.9A CN112697265A (en) | 2020-11-09 | 2020-11-09 | Self-adaptive multi-fractal method and equipment working condition monitoring device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011240325.9A CN112697265A (en) | 2020-11-09 | 2020-11-09 | Self-adaptive multi-fractal method and equipment working condition monitoring device |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112697265A true CN112697265A (en) | 2021-04-23 |
Family
ID=75506421
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011240325.9A Withdrawn CN112697265A (en) | 2020-11-09 | 2020-11-09 | Self-adaptive multi-fractal method and equipment working condition monitoring device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112697265A (en) |
-
2020
- 2020-11-09 CN CN202011240325.9A patent/CN112697265A/en not_active Withdrawn
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yang et al. | Vibration feature extraction techniques for fault diagnosis of rotating machinery: a literature survey | |
| CN113375939B (en) | Mechanical part fault diagnosis method based on SVD and VMD | |
| CN112697484A (en) | SSD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697483A (en) | LMD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697266A (en) | ITD and GZC machine state monitoring method and device | |
| CN112683393A (en) | LMD equipment fault diagnosis method and system | |
| CN116750290B (en) | Intelligent monitoring method for running state of bag type packaging machine | |
| CN112697472A (en) | WD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697479A (en) | ITD multi-scale fluctuation analysis state monitoring method and device | |
| CN106053069B (en) | A kind of SSD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN112697473A (en) | EMD and GZC machine state monitoring method and device | |
| CN112697263A (en) | EEMD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697477A (en) | EMD equipment fault diagnosis method and system | |
| CN112697482A (en) | VMD multi-scale fluctuation analysis state monitoring method and device | |
| CN112597958B (en) | Automatic identification method and system for rolling bearing fault | |
| CN105954031A (en) | Envelopment analysis method based on singular spectrum decomposition filtering | |
| CN112697470A (en) | SSD device fault diagnosis method and system | |
| CN112683395A (en) | Method and device for monitoring states of ELMD and GZC machines | |
| CN113435281A (en) | Ripple compensator fault diagnosis method based on improved HHT conversion | |
| CN112881046A (en) | VMD and GZC machine state monitoring method and device | |
| CN106153333B (en) | A kind of envelope Analysis Method based on wavelet decomposition filtering | |
| CN106096199B (en) | A kind of WT of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN112697265A (en) | Self-adaptive multi-fractal method and equipment working condition monitoring device | |
| CN106053060B (en) | A kind of envelope Analysis Method that filtering is decomposed based on nonlinear model | |
| CN112683392A (en) | ELMD multi-scale fluctuation analysis state monitoring method and device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WW01 | Invention patent application withdrawn after publication | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20210423 |