CN112697483A - LMD multi-scale fluctuation analysis state monitoring method and device - Google Patents

LMD multi-scale fluctuation analysis state monitoring method and device Download PDF

Info

Publication number
CN112697483A
CN112697483A CN202011240586.0A CN202011240586A CN112697483A CN 112697483 A CN112697483 A CN 112697483A CN 202011240586 A CN202011240586 A CN 202011240586A CN 112697483 A CN112697483 A CN 112697483A
Authority
CN
China
Prior art keywords
data
signal
lmd
function
fractal
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
Application number
CN202011240586.0A
Other languages
Chinese (zh)
Inventor
豆春玲
寇兴磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Kerishen Intelligent Technology Co ltd
Original Assignee
Shandong Kerishen Intelligent Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shandong Kerishen Intelligent Technology Co ltd filed Critical Shandong Kerishen Intelligent Technology Co ltd
Priority to CN202011240586.0A priority Critical patent/CN112697483A/en
Publication of CN112697483A publication Critical patent/CN112697483A/en
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00Subject matter not provided for in other groups of this subclass
    • G01M99/004Testing the effects of speed or acceleration
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00Testing of machine parts
    • G01M13/02Gearings; Transmission mechanisms
    • G01M13/021Gearings
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00Testing of machine parts
    • G01M13/02Gearings; Transmission mechanisms
    • G01M13/028Acoustic or vibration analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/15Correlation function computation including computation of convolution operations

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Acoustics & Sound (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a method and a device for monitoring an LMD multi-scale fluctuation analysis state. Decomposing the equipment vibration signal by using an LMD algorithm, removing noise components and trend terms by using a nonlinear discrimination algorithm, retaining fractal signal components, interpolating extreme points of each fractal signal component by using a Newton interpolation function, fitting an envelope by using a least square method, estimating instantaneous frequency by using a TEO algorithm, calculating corresponding instantaneous scale, determining a vibration signal detrending result according to an analysis scale, calculating a multi-fractal spectrum, extracting coordinates of a left end point, a right end point and the extreme points of the multi-fractal spectrum as characteristic parameters of the equipment running state, identifying the equipment running state, deploying the algorithm to an equipment state monitoring device, and being capable of accurately distinguishing the equipment running state.

Description

LMD multi-scale fluctuation analysis state monitoring method and device
Technical Field
The invention relates to the field of equipment state monitoring and fault diagnosis, in particular to an LMD multi-scale fluctuation analysis 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 an LMD multi-scale fluctuation analysis state (the method provided by the invention is abbreviated as MFDFOLmd) 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 an LMD multi-scale fluctuation analysis state is characterized by comprising the following steps: 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: office of adoptionThe partial Mean Decomposition (LMD) algorithm decomposes the signal x (k) into the sum of n components and a trend term, i.e.
Figure 398916DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the LMD algorithmn(k) Represents the trend term derived by the LMD algorithm, in this example, n = 10;
and step 3: eliminating noise components and trend terms from LMD 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 using a Newton interpolation function to the local maximum value and the local minimum value of cf(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
Figure 683267DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 158111DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 248427DEST_PATH_IMAGE004
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
Figure 430009DEST_PATH_IMAGE006
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 681999DEST_PATH_IMAGE007
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:
Figure 581822DEST_PATH_IMAGE008
Figure 210249DEST_PATH_IMAGE009
step 10: calculating a q-th order function:
Figure 511918DEST_PATH_IMAGE010
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:
Figure 934809DEST_PATH_IMAGE011
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 LMD algorithm in step 2 includes the following steps:
1) calculating a local mean function: determining all local extreme points n of original signals x (k)iCalculating two adjacent extreme points niAnd ni+1Average value m ofiI.e. by
Figure 384245DEST_PATH_IMAGE012
Average value m of all adjacent two extreme pointsiConnected by broken lines, and then smoothed by a moving average method to obtain a local mean function m11(k) (ii) a In this example, the smoothing step size in the moving average method is set to 5;
2) estimating the envelope value of the signal: using local extreme points niCalculating an envelope estimate ai
Figure 3314DEST_PATH_IMAGE013
Similarly, all two adjacent envelope estimates aiConnected by a broken line, and then smoothed by a moving average method to obtain an envelope estimateFunction a11(k);
3) The local mean function m11(k) Separated from the original signal x (k) to obtain
Figure 962086DEST_PATH_IMAGE014
4) By using h11(k) Division by an envelope estimation function a11(k) Thus to h11(k) Demodulating to obtain
Figure 759140DEST_PATH_IMAGE015
Ideally, s11(k) Is a pure frequency-modulated signal, i.e. its envelope estimation function a12(k) Satisfies a12(k) =1, if s11(k) If this condition is not satisfied, s is11(k) Repeating the above iteration process m times as new data until a pure frequency-modulated signal s is obtained1m(k) I.e. s1m(k) Satisfies-1. ltoreq. s1m(k) 1 or less, its envelope estimation function a1(m+1)(k) Satisfies a1(m+1)(k) =1, therefore, there are
Figure 430293DEST_PATH_IMAGE016
In the formula
Figure 603786DEST_PATH_IMAGE017
The condition for the iteration to terminate is
Figure 676784DEST_PATH_IMAGE018
In practical application, a variation delta can be set when 1-delta ≦ a is satisfied1m(k) When the sum of the delta and the delta is less than or equal to 1, the iteration is terminated; the amount of variation Δ =0.01 in this example;
5) the envelope signal is obtained by multiplying all envelope estimation functions generated in the iterative process
Figure 707057DEST_PATH_IMAGE019
6) Envelope signal a1(k) And a pure frequency-modulated signal s1m(k) The multiplication can obtain the 1 st component of the original signal
Figure 865506DEST_PATH_IMAGE020
7) The 1 st component c1(k) Separating from original signal x (k) to obtain a new signal r1(k) R is to1(k) Repeating the above steps as new data, and circulating n times until rn(k) Is a monotonic function
Figure 577110DEST_PATH_IMAGE021
To this end, the original signal x (k) is decomposed into n components and a monotonic function rn(k) To sum, i.e.
Figure 504614DEST_PATH_IMAGE022
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 ofExponential curve Hsurr(q) represents;
3) two parameters e are defined1And e2
Figure 768105DEST_PATH_IMAGE023
Figure 413850DEST_PATH_IMAGE024
If e is1And e2All less than 10%, the signal c (k) is discriminated as a noise component or a trend term, c (k) represents the signal component obtained by the LMD algorithm.
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):
Figure 663566DEST_PATH_IMAGE025
Figure 773473DEST_PATH_IMAGE026
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:
Figure 83232DEST_PATH_IMAGE027
Figure 216273DEST_PATH_IMAGE028
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:
Figure 269680DEST_PATH_IMAGE029
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:
Figure 171777DEST_PATH_IMAGE011
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+ccIs 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
Figure 652437DEST_PATH_IMAGE030
3) Let J pair akPartial derivative of
Figure 7195DEST_PATH_IMAGE031
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure 926609DEST_PATH_IMAGE032
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 the corresponding curve when J is minimumThe line has the form of rk(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:
Figure 620896DEST_PATH_IMAGE033
based on the above method for monitoring the state of the LMD multi-scale fluctuation analysis, the apparatus for 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 realizing 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 LMD method is adopted to decompose the vibration signal, noise components and trend items are removed according to the 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 LMD 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.
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 mfdfaold 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 by respectively using the MFDF, MFDFEMd and MFDFOLmd 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, MFDFEMd and MFDFOLmd 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, MFDFAMd based on correlation filtering, and MFDFAMd 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 MFDFOLmd 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 mfdfaold for the four gear box states in the embodiment of the present invention, and the symbols "circle", "square", "plus", and "diamond" represent the normal, light scratch, heavy scratch, and broken tooth gear states, respectively.
Detailed Description
An embodiment, as shown in fig. 1 and fig. 2, is a method for monitoring an LMD multi-scale fluctuation analysis state, wherein: 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 Local Mean Decomposition (LMD) algorithm is used to decompose the signal x (k) into the sum of n components and a trend term, i.e. the sum
Figure 334774DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the LMD algorithmn(k) Represents the trend term derived by the LMD algorithm, in this example, n = 10;
and step 3: eliminating noise components and trend terms from LMD 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 using a Newton interpolation function to the local maximum value and the local minimum value of cf(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
Figure 380090DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 831757DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 380550DEST_PATH_IMAGE004
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
Figure 530908DEST_PATH_IMAGE034
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 797942DEST_PATH_IMAGE007
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:
Figure 59159DEST_PATH_IMAGE008
Figure 462458DEST_PATH_IMAGE035
step 10: calculating a q-th order function:
Figure 783718DEST_PATH_IMAGE010
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:
Figure 538048DEST_PATH_IMAGE011
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 LMD algorithm in the step 2 comprises the following steps:
1) calculating a local mean function: determining all local extreme points n of original signals x (k)iCalculating two adjacent extreme points niAnd ni+1Average value m ofiI.e. by
Figure 337376DEST_PATH_IMAGE012
Average value m of all adjacent two extreme pointsiConnected by broken lines, and then smoothed by a moving average method to obtain a local mean function m11(k) (ii) a In this example, the smoothing step size in the moving average method is set to 5;
2) estimating the envelope value of the signal: using local extreme points niCalculating an envelope estimate ai
Figure 923078DEST_PATH_IMAGE013
Similarly, all two adjacent envelope estimates aiConnected by a broken line, and then smoothed by a moving average method to obtain an envelope estimation function a11(k);
3) The local mean function m11(k) Separated from the original signal x (k) to obtain
Figure 415240DEST_PATH_IMAGE036
4) By using h11(k) Division by an envelope estimation function a11(k) Thus to h11(k) Demodulating to obtain
Figure 391286DEST_PATH_IMAGE037
Ideally, s11(k) Is a pure frequency-modulated signal, i.e. its envelope estimation function a12(k) Satisfies a12(k) =1, if s11(k) If this condition is not satisfied, s is11(k) Repeating the above iteration process m times as new data until a pure frequency-modulated signal s is obtained1m(k) I.e. s1m(k) Satisfies-1. ltoreq. s1m(k) 1 or less, its envelope estimation function a1(m+1)(k) Satisfies a1(m+1)(k) =1, therefore, there are
Figure 994306DEST_PATH_IMAGE038
In the formula
Figure 637777DEST_PATH_IMAGE039
The condition for the iteration to terminate is
Figure 35260DEST_PATH_IMAGE018
In practical application, a variation delta can be set when 1-delta ≦ a is satisfied1m(k) When the sum of the delta and the delta is less than or equal to 1, the iteration is terminated; the amount of variation Δ =0.01 in this example;
5) the envelope signal is obtained by multiplying all envelope estimation functions generated in the iterative process
Figure 764181DEST_PATH_IMAGE019
6) Envelope signal a1(k) And a pure frequency-modulated signal s1m(k) The multiplication can obtain the 1 st component of the original signal
Figure 842996DEST_PATH_IMAGE040
7) The 1 st component c1(k) Separating from original signal x (k) to obtain a new signal r1(k) R is to1(k) Repeating the above steps as new data, and circulating n times until rn(k) Is a monotonic function
Figure 403290DEST_PATH_IMAGE021
To this end, the original signal x (k) is decomposed into n components and a singleTransfer function rn(k) To sum, i.e.
Figure 971675DEST_PATH_IMAGE022
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
Figure 187892DEST_PATH_IMAGE041
Figure 132715DEST_PATH_IMAGE042
If e is1And e2All less than 10%, the signal c (k) is discriminated as a noise component or a trend term, c (k) represents the signal component obtained by the LMD algorithm.
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):
Figure 485199DEST_PATH_IMAGE025
Figure 490064DEST_PATH_IMAGE026
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:
Figure 927998DEST_PATH_IMAGE027
Figure 410932DEST_PATH_IMAGE028
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:
Figure 883502DEST_PATH_IMAGE029
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:
Figure 793689DEST_PATH_IMAGE011
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+ccIs 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
Figure 781237DEST_PATH_IMAGE030
3) Let J pair akPartial derivative of
Figure 5545DEST_PATH_IMAGE031
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure 400797DEST_PATH_IMAGE032
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:
Figure 685148DEST_PATH_IMAGE033
based on the above method for monitoring the state of the LMD multi-scale fluctuation analysis, the apparatus for 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 realizing the algorithm.
Experiment 1 the performance of the algorithm of the present invention was verified using a multi-fractal simulation signal generated by a multi-fractal cascade model.
Firstly, a multi-fractal cascade model is adopted
Figure DEST_PATH_IMAGE043
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 mfdfaold, and the result is shown in fig. 5. As can be seen from fig. 5, the instantaneous frequencies calculated by the mfdfaold are all positive frequencies, so the mfdfaold analysis result conforms to the actual situation. Next, a multi-fractal spectrum of the multi-fractal simulation signal is calculated using MFDFA, MFDFAemd, and mfdfaold, 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 multi-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 multi-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 multi-fractal spectrum obtained from MFDFAolmd and the theoretical value is 0.023, and the average of the relative errors is 5.13%, so that the average of the absolute errors between the multi-fractal spectrum obtained from MFDFAolmd is reduced by 64.06%, the average of the relative errors is reduced by 57.99%, the average of the absolute errors between the multi-fractal spectrum obtained from MFDFAolmd is reduced by 34.29%, and the average of the relative errors is reduced by 21.92%. 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. Respectively adopting MFDFA, MFDFAemd and MFDFAemd to calculate noiseThe result of the multi-fractal spectrum of the multi-fractal simulation signal is shown in fig. 8. According to the result 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 mfdfaamd 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 MFDFAolmd and the theoretical value is 0.031, and the mean relative error value is 5.36%, so that the mean absolute error value of the multi-fractal spectrum obtained from MFDFAolmd is reduced by 63.10%, and the mean relative error value is reduced by 54.61%. It follows that MFDFAolmd 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 noisy multi-fractal simulation signals by using the MFDFA, the mfdfame based on the correlation filtering, and the mfdfaold based on the correlation filtering 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 MFDFAolmd 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.
Experiment 2 the performance of the algorithm of the invention was verified using gearbox experimental signals.
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 MFDFOLmd 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 other three gearbox states, 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, MFDFEMd and MFDFOLmd 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 mfdfaold method, so that the gear box state identification rate is 100%. It can be seen that the MFDFOLmd method can improve the accuracy of the state identification of the gearbox by 50%.
From the test results, it was assumed after analysis that:
1) the vibration signal is adaptively decomposed by adopting an LMD (local mean decomposition) method, the noise component and the trend item 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 LMD 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. A method for monitoring an LMD multi-scale fluctuation analysis state is characterized by comprising the following steps: 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 Local Mean Decomposition (LMD) algorithm is used to decompose the signal x (k) into the sum of n components and a trend term, i.e. the sum
Figure 955166DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the LMD algorithmn(k) Represents a trend term derived by the LMD algorithm;
and step 3: eliminating noise components and trend terms from LMD 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 using a Newton interpolation function to the local maximum value and the local minimum value of cf(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
Figure 934623DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 776677DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 171887DEST_PATH_IMAGE004
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
Figure 48576DEST_PATH_IMAGE006
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 871038DEST_PATH_IMAGE007
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:
Figure 217967DEST_PATH_IMAGE008
Figure 416867DEST_PATH_IMAGE009
step 10: calculating a q-th order function:
Figure 148062DEST_PATH_IMAGE010
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:
Figure 141426DEST_PATH_IMAGE011
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 LMD multi-scale fluctuation analysis state monitoring method according to claim 1, wherein: the LMD algorithm in the step 2 comprises the following steps:
1) calculating a local mean function: determining all local extreme points n of original signals x (k)iCalculating two adjacent extreme points niAnd ni+1Average value m ofiI.e. by
Figure 958072DEST_PATH_IMAGE012
Average value m of all adjacent two extreme pointsiConnected by broken lines, and then smoothed by a moving average method to obtain a local mean function m11(k);
2) Estimating the envelope value of the signal: using local extreme points niCalculating an envelope estimate ai
Figure 695084DEST_PATH_IMAGE013
Similarly, all two adjacent envelope estimates aiConnected by a broken line, and then smoothed by a moving average method to obtain an envelope estimation function a11(k);
3) The local mean function m11(k) Separated from the original signal x (k) to obtain
Figure 280786DEST_PATH_IMAGE014
4) By using h11(k) Division by an envelope estimation function a11(k) Thus to h11(k) Demodulating to obtain
Figure 710631DEST_PATH_IMAGE015
Ideally, s11(k) Is a pure frequency-modulated signal, i.e. its envelope estimation function a12(k) Satisfies a12(k) =1, if s11(k) If this condition is not satisfied, s is11(k) Repeating the above iteration process m times as new data until a pure frequency-modulated signal s is obtained1m(k) I.e. s1m(k) Satisfies-1. ltoreq. s1m(k) 1 or less, its envelope estimation function a1(m+1)(k) Satisfies a1(m+1)(k) =1, therefore, there are
Figure 14573DEST_PATH_IMAGE016
In the formula
Figure 289697DEST_PATH_IMAGE017
The condition for the iteration to terminate is
Figure 995485DEST_PATH_IMAGE018
In practical application, a variation delta can be set when 1-delta ≦ a is satisfied1m(k) When the sum of the delta and the delta is less than or equal to 1, the iteration is terminated;
5) the envelope signal is obtained by multiplying all envelope estimation functions generated in the iterative process
Figure 330651DEST_PATH_IMAGE019
6) Envelope signal a1(k) And a pure frequency-modulated signal s1m(k) The multiplication can obtain the 1 st component of the original signal
Figure 121889DEST_PATH_IMAGE020
7) The 1 st component c1(k) Separating from original signal x (k) to obtain a new signal r1(k) R is to1(k) Repeating the above steps as new data, and circulating n times until rn(k) Is a monotonic function
Figure 263021DEST_PATH_IMAGE021
To this end, the original signal x (k) is decomposed into n components and a monotonic function rn(k) To sum, i.e.
Figure 760998DEST_PATH_IMAGE022
3. The LMD multi-scale fluctuation analysis 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 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
Figure 594962DEST_PATH_IMAGE023
Figure 545600DEST_PATH_IMAGE024
If e is1And e2All less than 10%, the signal c (k) is discriminated as a noise component or a trend term, c (k) represents the signal component obtained by the LMD algorithm.
4. A LMD multi-scale fluctuation analysis state monitoring method according to claim 3, characterized in that: the data rearrangement operation in the step 1) comprises the following steps: randomly randomizing the order of the components c (k).
5. A LMD multi-scale fluctuation analysis state monitoring method according to 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. A LMD multi-scale fluctuation analysis state monitoring method according to claim 3, characterized in that: the MFDF method in the step 2) comprises the following steps:
1) structure of the devicex (k), k =1,2, …, N, profileY(i):
Figure 490423DEST_PATH_IMAGE025
Figure 842907DEST_PATH_IMAGE026
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:
Figure 847772DEST_PATH_IMAGE027
Figure 348023DEST_PATH_IMAGE028
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:
Figure 768640DEST_PATH_IMAGE029
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:
Figure 303527DEST_PATH_IMAGE011
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+ccIs constant, whereby the slope of the straight line can be obtainedH(q)。
7. The LMD multi-scale fluctuation analysis 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 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
Figure 416976DEST_PATH_IMAGE030
3) Let J pair akPartial derivative of
Figure 404524DEST_PATH_IMAGE031
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure 691149DEST_PATH_IMAGE032
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 LMD multi-scale fluctuation analysis state monitoring method according to claim 1, wherein: 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:
Figure DEST_PATH_IMAGE033
9. apparatus for implementing a method for LMD multi-scale fluctuation analysis status monitoring 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, and the signal analysis software is installed on the notebook computer.
CN202011240586.0A 2020-11-09 2020-11-09 LMD multi-scale fluctuation analysis state monitoring method and device Withdrawn CN112697483A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011240586.0A CN112697483A (en) 2020-11-09 2020-11-09 LMD multi-scale fluctuation analysis state monitoring method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011240586.0A CN112697483A (en) 2020-11-09 2020-11-09 LMD multi-scale fluctuation analysis state monitoring method and device

Publications (1)

Publication Number Publication Date
CN112697483A true CN112697483A (en) 2021-04-23

Family

ID=75506413

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011240586.0A Withdrawn CN112697483A (en) 2020-11-09 2020-11-09 LMD multi-scale fluctuation analysis state monitoring method and device

Country Status (1)

Country Link
CN (1) CN112697483A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115015682A (en) * 2022-08-09 2022-09-06 南京佑友软件技术有限公司 Real-time online monitoring method for power quality
CN116776168A (en) * 2023-08-22 2023-09-19 惠州帝恩科技有限公司 Intelligent analysis method and system for production data of reagent tubes

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115015682A (en) * 2022-08-09 2022-09-06 南京佑友软件技术有限公司 Real-time online monitoring method for power quality
CN115015682B (en) * 2022-08-09 2022-11-08 南京佑友软件技术有限公司 Real-time online monitoring method for power quality
CN116776168A (en) * 2023-08-22 2023-09-19 惠州帝恩科技有限公司 Intelligent analysis method and system for production data of reagent tubes
CN116776168B (en) * 2023-08-22 2023-11-21 惠州帝恩科技有限公司 Intelligent analysis method and system for production data of reagent tubes

Similar Documents

Publication Publication Date Title
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
CN112697483A (en) LMD multi-scale fluctuation analysis state monitoring method and device
CN113375939B (en) Mechanical part fault diagnosis method based on SVD and VMD
CN108801630B (en) Gear fault diagnosis method for single-channel blind source separation
CN110991055B (en) Residual life prediction system for rotary mechanical products
CN112697484A (en) SSD multi-scale fluctuation analysis state monitoring method and device
CN112683393A (en) LMD equipment fault diagnosis method and system
CN112697266A (en) ITD and GZC machine state monitoring method and device
Jiang et al. A novel method for self-adaptive feature extraction using scaling crossover characteristics of signals and combining with LS-SVM for multi-fault diagnosis of gearbox
CN112697472A (en) WD multi-scale fluctuation analysis state monitoring method and device
CN112697479A (en) ITD multi-scale fluctuation analysis state monitoring method and device
CN112697482A (en) VMD multi-scale fluctuation analysis state monitoring method and device
CN112697473A (en) EMD and GZC machine state monitoring method and device
CN112697263A (en) EEMD multi-scale fluctuation analysis state monitoring method and device
CN112881046A (en) VMD and GZC machine state monitoring method and device
CN112697477A (en) EMD equipment fault diagnosis method and system
CN112697470A (en) SSD device fault diagnosis method and system
CN112683395A (en) Method and device for monitoring states of ELMD and GZC machines
CN106096199B (en) A kind of WT of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method
CN112683392A (en) ELMD multi-scale fluctuation analysis state monitoring method and device
CN106153333B (en) A kind of envelope Analysis Method based on wavelet decomposition filtering
CN112697474A (en) NMD multi-scale fluctuation analysis state monitoring method and device
CN112697265A (en) Self-adaptive multi-fractal method and equipment working condition monitoring device
CN112629849A (en) FDM multi-scale fluctuation analysis state monitoring method and device
CN112697471A (en) SSD 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: 20210423