CN112683395A - Method and device for monitoring states of ELMD and GZC machines - Google Patents
Method and device for monitoring states of ELMD and GZC machines Download PDFInfo
- Publication number
- CN112683395A CN112683395A CN202011240587.5A CN202011240587A CN112683395A CN 112683395 A CN112683395 A CN 112683395A CN 202011240587 A CN202011240587 A CN 202011240587A CN 112683395 A CN112683395 A CN 112683395A
- Authority
- CN
- China
- Prior art keywords
- data
- elmd
- signal
- curve
- local
- 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
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention discloses a monitoring method and a device for ELMD and GZC machine states, which decompose equipment vibration signals by using an ELMD algorithm, remove decomposed noise components and trend terms by using a nonlinear discrimination algorithm, reserve fractal signal components, interpolate extreme points by using a rational spline interpolation function, fit envelopes by using a least square method, separate a frequency modulation part, estimate instantaneous frequency by using a GZC algorithm and calculate corresponding instantaneous scales, automatically determine a de-trend result of the vibration signals according to the analysis scales, calculate a multi-fractal spectrum of de-trend signals, extract left end points, right end points and extreme point coordinates of the multi-fractal spectrum as characteristic parameters of equipment running states, identify the equipment running states, deploy the algorithm to an equipment state monitoring device, can accurately distinguish the equipment running states, and the equipment state monitoring device has good flexibility and portability, is convenient for engineering application.
Description
Technical Field
The invention relates to the field of equipment state monitoring and fault diagnosis, in particular to an ELMD and GZC machine state monitoring method and 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
The invention provides a method and a device for monitoring the states of ELMD and GZC machines (the method provided by the invention is abbreviated as MFDFooelmd) aiming at the defects. 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: a method for monitoring the state of ELMD and GZC machines is characterized in that: the method comprises 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. the sum of n components and a trend term, by using an Ensemble Local Mean Decomposition (ELMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the ELMD algorithmn(k) Represents the trend term obtained by the ELMD algorithm, in this example, n = 10;
and step 3: eliminating noise components and trend terms from ELMD 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 rational 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: FM was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
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), in this case q is taken in the range (-5, 5);
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), and describing the running state of the equipment by using the 3 parameters;
step 16: the method in the steps is deployed on a state monitoring device to monitor the state of equipment.
Further, the ELMD algorithm in step 2 includes the following steps:
1) to data x0(k) Adding a white noise sequence produces a new data xj(k) ：
Std[x0(k)]Representative data x0(k) Standard deviation of (Wn)j(k) Representative wnjThe kth data, wn injRepresents the jth randomly generated white noise sequence, wnjThe amplitude is 1, j is more than or equal to 1 and less than or equal to K; x is the number of0(k) Represents x (k) in step 2 of claim 1; in this example, K = 100;
2) for xj(k) Performing local mean decomposition to obtain n components and a trend term
cij(k) Represents a pair xj(k) The i-th component, r, resulting from performing a local mean decompositionnj(k) Represents a pair xj(k) Performing a trend term obtained by local mean decomposition;
3) calculating the average value of the decomposition results of K times
ci(k) Represents a pair x0(k) The ith component, r, obtained by performing a local mean decomposition of the setn(k) Represents a pair x0(k) And performing a trend term obtained by the set local mean decomposition.
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(ii) a In this example, m = 1;
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 minima of the 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) representsA 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 GZC method in step 6 includes the following steps:
1) determining local maximum, minimum and zero points of the signal c (k), k =1,2, …, N, c (k) = FMm(k)；
2) Counting the time intervals T between two continuous local maximum points, two continuous local minimum points, two continuous rising zeros and two continuous falling zeros according to the time sequence1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, counting neighbouring offices in time orderTime intervals T between partial maximum point and zero-point, adjacent zero-point and local minimum point, adjacent local minimum point and zero-point, and adjacent zero-point and local maximum point3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
The device based on the ELMD and GZC machine state monitoring method is characterized in that: the state monitoring device in the step 16 comprises the following parts: 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 realizing the algorithm.
The method comprises the following steps:
step 1), the following steps: collecting a vibration signal;
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 utilizing the GZC;
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 the vertical coordinates of the left end point, the right end point and the extreme point of the multi-fractal spectrum, and taking the three parameters as the 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 ELMD 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) a frequency modulation part for separating the signal components, and estimating the instantaneous frequency of the signal components by using the GZC, 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 ELMD method is utilized to automatically determine the type of the signal trend, ensure the continuity of the signal trend and effectively solve the defects of the prior art;
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 MFDFAoelmd 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 MFDFaeelmd methods in the embodiment of the present invention;
FIG. 7 is a diagram illustrating the calculation results of two non-linear discriminant parameters, wherein the symbols "circle" and "square" represent e1 and e2, respectively;
FIG. 8 is a comparison diagram of analysis results of noisy multi-fractal simulation signals respectively using MFDF, MFDFEMd, and MFDFaeelmd 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 ELMD 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 MFDFAElmd based on correlation filtering in the embodiment of the present invention;
FIG. 11 is a vibration signal of five gear boxes in the embodiment of the invention, wherein (a) - (e) respectively represent normal, abrasion, pitting, tooth breakage and abrasion + pitting gear states;
FIG. 12 is a multi-fractal spectrum of the five gearbox vibration signals obtained using MFDFA in an embodiment of the present invention;
FIG. 13 is a multi-fractal spectrum of vibration signals of the five gearboxes obtained by using MFDFame in the embodiment of the present invention;
FIG. 14 is a multi-fractal spectrum of vibration signals of the five gearboxes obtained using MFDFaeelmd in an embodiment of the present invention;
FIG. 15 shows the results of the classification of the left, right and extreme coordinates of the multi-fractal spectrum obtained from MFDFA into the five gear box states in the example of the present invention, where the "circle", "square", "plus", "diamond" and "left triangle" symbols represent the normal, wear, pitting, tooth breakage and wear + pitting gear states, respectively;
FIG. 16 shows the results of the classification of the left, right, and extreme coordinates of the multi-fractal spectrum obtained from MFDFame on the states of the five gearboxes according to the present invention, where the "circle", "square", "plus", "diamond", and "left triangle" symbols represent the normal, wear, pitting, tooth breakage, and wear + pitting gear states, respectively;
fig. 17 shows the results of classifying the five gear box states by the coordinates of the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained from MFDFAoelmd in the embodiment of the present invention, and the symbols "circle", "square", "plus", "diamond" and "left triangle" represent the normal, worn, pitting, broken tooth and worn + pitting gear states, respectively.
Detailed Description
Embodiment, as shown in fig. 1 and fig. 2, a method for monitoring states of ELMD and GZC machines is characterized in that: the method comprises 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. the sum of n components and a trend term, by using an Ensemble Local Mean Decomposition (ELMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the ELMD algorithmn(k) Represents the trend term obtained by the ELMD algorithm, in this example, n = 10;
and step 3: eliminating noise components and trend terms from ELMD 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 rational 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 | denotes the absolute of xFor the value;
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: FM was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
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), in this case q is taken in the range (-5, 5);
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), and describing the running state of the equipment by using the 3 parameters;
step 16: the method in the steps is deployed on a state monitoring device to monitor the state of equipment.
The ELMD algorithm in the step 2 comprises the following steps:
1) to data x0(k) Adding a white noise sequence produces a new data xj(k) ：
Std[x0(k)]Representative data x0(k) Standard deviation of (Wn)j(k) Representative wnjThe kth data, wn injRepresents the jth randomly generated white noise sequence, wnjThe amplitude is 1, j is more than or equal to 1 and less than or equal to K; x is the number of0(k) Represents x (k) in step 2 of claim 1; in this example, K = 100;
2) for xj(k) Performing local mean decomposition to obtain n components and a trend term
cij(k) Represents a pair xj(k) The i-th component, r, resulting from performing a local mean decompositionnj(k) Represents a pair xj(k) Performing a trend term obtained by local mean decomposition;
3) calculating the average value of the decomposition results of K times
ci(k) Represents a pair x0(k) The ith component, r, obtained by performing a local mean decomposition of the setn(k) Represents a pair x0(k) And performing a trend term obtained by the set local mean decomposition.
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 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，
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(ii) a In this example, m = 1;
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 minima of the 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 GZC method in the step 6 comprises the following steps:
1) determining local maxima, local maxima of signals c (k), k =1,2, …, NMinimum point and zero point, c (k) = FMm(k)；
2) Counting the time intervals T between two continuous local maximum points, two continuous local minimum points, two continuous rising zeros and two continuous falling zeros according to the time sequence1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and descending zero points, adjacent descending zero points and local minimum points, adjacent local minimum points and ascending zero points, and adjacent ascending zero points and local maximum points according to the time sequence3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
Based on the above-mentioned ELMD and GZC machine status monitoring method, the status monitoring device in step 16 includes the following parts: 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 realizing the algorithm.
Firstly, a multi-fractal cascade model is adopted
The generated multi-fractal simulation signal verifies the performance of MFDFA, MFDFAemd and MFDFAelmd. 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 MFDFAoelmd, and the result is shown in fig. 5. As can be seen from fig. 5, the instantaneous frequencies calculated by MFDFAoelmd are all positive frequencies, so the MFDFAoelmd analysis result conforms to the actual situation. Next, a multi-fractal spectrum of the multi-fractal simulation signal is calculated using MFDFA, MFDFAemd, and mfdfaeelmd, respectively, and the result is shown in fig. 6. According to the results shown in fig. 6, it is 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%, the average of the absolute errors between the fractal spectrum obtained from mfdfaeelmd and the theoretical value is 0.025, and the average of the relative errors is 5.29%, so that the average of the absolute errors between the fractal spectrum obtained from mfdfaeelmd is 60.94% smaller than the average of the absolute errors between the fractal spectrum obtained from MFDFAemd, the average of the relative errors is 56.67%, the average of the fractal spectrum obtained from mfdfaeelmd is 28.57% smaller than the average of the absolute errors between the fractal spectrum obtained from MFDFAemd, and the average of the relative errors is 19.48%. 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 mfdfaeelmd, respectively, and the result is shown in fig. 8. According to the results shown in FIG. 8, the fractal spectrum obtained from MFDFA completely deviates from the theoretical value, the mean absolute error value of the fractal spectrum obtained from MFDFAEmd from 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 MFDFAoelmd and the theoretical value is 0.036, and the mean relative error value is 5.52%, so that the multi-fractal spectrum obtained from MFDFAoelmd is reduced by 57.14% and the mean relative error value is reduced by 53.26% compared with the mean absolute error value of the multi-fractal spectrum obtained from MFDFAemd. It follows that MFDFAoelmd has better noise immunity than MFDFA and MFDFAemd. 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 mfdfaeelmd based on correlation filtering, respectively, in the embodiment of the present invention. As can be seen from fig. 10, the analysis result of the MFDFAemd based on the correlation filtering and the mfdfaeelmd based on the correlation filtering on the noisy multi-fractal simulation signal completely deviates from the theoretical value, so that the fractal structure of the original signal is easily damaged by the correlation filtering method.
The gearbox vibration data used in the invention is from a gearbox fault simulation experiment table. This experiment simulates three single point gear failures: wear, pitting and tooth breakage, and one compound gear failure: wear + pitting. The collected vibration signals comprise five gear running states of normal, abrasion, pitting corrosion, tooth breakage and abrasion + pitting corrosion. The motor speed was 1500RPM, the vibration signal sampling frequency was 5120Hz, and 5 segments of data with a length of 10000 points were collected at each gear state, and the five gearbox vibration signals are shown in fig. 11. Firstly, the five gearbox vibration signals are analyzed by adopting an MFDF method, and the multi-fractal spectrums corresponding to the five gearbox vibration signals are obtained as shown in figure 12, so that the multi-fractal spectrums corresponding to normal, abrasion, broken tooth and compound fault states are seriously overlapped, and the multi-fractal spectrum corresponding to a pitting state is abnormal in shape. Then, the five gearbox vibration signals are analyzed by using an MFDFame method, and the multi-fractal spectrums corresponding to the five gearbox vibration signals are obtained as shown in FIG. 13, so that the multi-fractal spectrums corresponding to normal, abrasion, broken teeth and compound fault states are seriously overlapped. Finally, the vibration signals of the five gearboxes are analyzed by using an MFMFMFDOelmd method, and the multi-fractal spectrums corresponding to the vibration signals of the five gearboxes are obtained as shown in fig. 14, so that the multi-fractal spectrums corresponding to the five gear states can be separated. 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, MFDFEMd and MFDFOelmd methods are respectively extracted to classify the states of the five gear boxes, and the results are respectively shown in FIGS. 15-17. As can be seen from fig. 15, the pitting gear 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 MFDFA method, but the remaining four gear states cannot be distinguished, so the gear box state recognition rate is 20%. As can be seen from fig. 16, the left end point, the right end point, and the singular index corresponding to the extreme point of the multi-fractal spectrum obtained by the MFDFAemd method can correctly distinguish the states of the pitting gears, but cannot distinguish the remaining four gear states, so the gear box state recognition rate is 20%. As can be seen from fig. 17, the five gear 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 MFDFAoelmd method, and therefore the gear box state recognition rate is 100%. It can be seen that the MFDFooelmd method can improve the accuracy of the state identification of the gearbox by 80%.
From the test results, it was assumed after analysis that:
1) the vibration signal is adaptively decomposed by adopting an ELMD 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 for separating the signal components estimates the instantaneous scale of the signal components by utilizing the GZC, so that the instantaneous frequency can be ensured to keep a positive value, and the negative frequency phenomenon is avoided;
3) the fluctuation analysis is carried out according to the instantaneous scale of the signal component, so that the defect of manually setting the scale is avoided;
4) the ELMD method is utilized to automatically determine the type of the signal trend, ensure the continuity of the signal trend and effectively solve the defects of the prior art;
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. A method for monitoring the state of ELMD and GZC machines is characterized in that: the method comprises 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. the sum of n components and a trend term, by using an Ensemble Local Mean Decomposition (ELMD) algorithm
Wherein c isi(k) Representing the i-th component, r, obtained by the ELMD algorithmn(k) Represents a trend term derived from the ELMD algorithm;
and step 3: eliminating noise components and trend terms from ELMD 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 rational 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: FM was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
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 ELMD and GZC machine state monitoring method according to claim 1, wherein: the ELMD algorithm in the step 2 comprises the following steps:
1) to data x0(k) Adding a white noise sequence produces a new data xj(k) ：
Std[x0(k)]Representative data x0(k) Standard deviation of (Wn)j(k) Representative wnjThe kth data, wn injRepresents the jth randomly generated white noise sequence, wnjThe amplitude is 1, j is more than or equal to 1 and less than or equal to K; x is the number of0(k) Represents x (k) in step 2 of claim 1;
2) for xj(k) Performing local mean decomposition to obtain n components and a trend term
cij(k) Represents a pair xj(k) The i-th component, r, resulting from performing a local mean decompositionnj(k) Represents a pair xj(k) Performing a trend term obtained by local mean decomposition;
3) calculating the average value of the decomposition results of K times
ci(k) Represents a pair x0(k) The ith component, r, obtained by performing a local mean decomposition of the setn(k) Represents a pair x0(k) And performing a trend term obtained by the set local mean decomposition.
3. The ELMD and GZC machine state monitoring method according to claim 1, wherein: 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 multiple fractal separatelyRemoving trend Fluctuation Analysis (MFDF) 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，
4. An ELMD and GZC machine state monitoring device according to claim 3, wherein: the data rearrangement operation in the step 1) comprises the following steps: randomly randomizing the order of the components c (k).
5. The ELMD and GZC machine state monitoring method according to claim 3, wherein: 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 ELMD and GZC machine state monitoring method according to claim 3, wherein: 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 ELMD and GZC machine state monitoring method according to claim 1, wherein: 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 minima of the 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 ELMD and GZC machine state monitoring method according to claim 1, wherein: the GZC method in the step 6 comprises the following steps:
1) determining local maximum, minimum and zero points of the signal c (k), k =1,2, …, N, c (k) = FMm(k)；
2) Counting the time intervals T between two continuous local maximum points, two continuous local minimum points, two continuous rising zeros and two continuous falling zeros according to the time sequence1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, counting adjacent local maximum points and dips in time orderZero point, adjacent falling zero point and local minimum point, adjacent local minimum point and rising zero point, and time interval T between adjacent rising zero point and local maximum point3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
9. Apparatus for implementing an ELMD and GZC machine condition monitoring method according to any of claims 1 to 8, characterized in that: the device comprises the following parts: 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 realizing the algorithm.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011240587.5A CN112683395A (en) | 2020-11-09 | 2020-11-09 | Method and device for monitoring states of ELMD and GZC machines |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011240587.5A CN112683395A (en) | 2020-11-09 | 2020-11-09 | Method and device for monitoring states of ELMD and GZC machines |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112683395A true CN112683395A (en) | 2021-04-20 |
Family
ID=75445937
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011240587.5A Withdrawn CN112683395A (en) | 2020-11-09 | 2020-11-09 | Method and device for monitoring states of ELMD and GZC machines |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112683395A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116992393A (en) * | 2023-09-27 | 2023-11-03 | 联通(江苏)产业互联网有限公司 | Safety production monitoring method based on industrial Internet of things |
-
2020
- 2020-11-09 CN CN202011240587.5A patent/CN112683395A/en not_active Withdrawn
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116992393A (en) * | 2023-09-27 | 2023-11-03 | 联通(江苏)产业互联网有限公司 | Safety production monitoring method based on industrial Internet of things |
| CN116992393B (en) * | 2023-09-27 | 2023-12-22 | 联通(江苏)产业互联网有限公司 | Safety production monitoring method based on industrial Internet of things |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106198015B (en) | A kind of VMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN113375939B (en) | Mechanical part fault diagnosis method based on SVD and VMD | |
| CN108801630B (en) | Gear fault diagnosis method for single-channel blind source separation | |
| CN110046476B (en) | Ternary binary fractal wavelet sparse diagnosis method for rolling bearing faults | |
| CN112683393A (en) | LMD equipment fault diagnosis method and system | |
| CN112697484A (en) | SSD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697266A (en) | ITD and GZC machine state monitoring method and device | |
| CN112697483A (en) | LMD multi-scale fluctuation analysis state monitoring method and device | |
| Laval et al. | Amplitude and phase interaction in Hilbert demodulation of vibration signals: Natural gear wear modeling and time tracking for condition monitoring | |
| Zhang et al. | Improved local cepstrum and its applications for gearbox and rolling bearing fault detection | |
| CN112683395A (en) | Method and device for monitoring states of ELMD and GZC machines | |
| CN112881046A (en) | VMD and GZC machine state monitoring method and device | |
| CN106053069B (en) | A kind of SSD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN112697477A (en) | EMD equipment fault diagnosis method and system | |
| CN112697470A (en) | SSD device fault diagnosis method and system | |
| CN112697473A (en) | EMD and GZC machine state monitoring method and device | |
| CN112697472A (en) | WD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697482A (en) | VMD multi-scale fluctuation analysis state monitoring method and device | |
| CN105954031A (en) | Envelopment analysis method based on singular spectrum decomposition filtering | |
| CN112697479A (en) | ITD multi-scale fluctuation analysis state monitoring method and device | |
| CN112697263A (en) | EEMD multi-scale fluctuation analysis state monitoring method and device | |
| CN106096199B (en) | A kind of WT of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN112683391A (en) | EEMD and GZC machine state monitoring method and device | |
| CN112697481A (en) | NMD and GZC machine state monitoring method and device | |
| CN112629848A (en) | EFDM and GZC machine 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: 20210420 |