CN112697475A - ITD equipment fault diagnosis method and system - Google Patents

ITD equipment fault diagnosis method and system Download PDF

Info

Publication number
CN112697475A
CN112697475A CN202011238927.0A CN202011238927A CN112697475A CN 112697475 A CN112697475 A CN 112697475A CN 202011238927 A CN202011238927 A CN 202011238927A CN 112697475 A CN112697475 A CN 112697475A
Authority
CN
China
Prior art keywords
data
signal
itd
curve
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
CN202011238927.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 CN202011238927.0A priority Critical patent/CN112697475A/en
Publication of CN112697475A publication Critical patent/CN112697475A/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)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Acoustics & Sound (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a method and a system for diagnosing the fault of ITD equipment, which decompose the vibration signal of the equipment by using an ITD algorithm, reserve the component of a fractal signal by using nonlinear discrimination algorithm noise and a trend term, interpolate extreme points by using a cubic spline interpolation function, fit an envelope by using a least square method, separate a frequency modulation part, estimate instantaneous frequency by using a DQ algorithm and calculate corresponding instantaneous scale, determine the detrending result of the vibration signal according to the analysis scale, calculate the multi-fractal spectrum of the detrending signal, extract the left end point, the right end point and the extreme point coordinates of the multi-fractal spectrum as the characteristic parameters of the running state of the equipment, identify the running state of the equipment, and deploy the algorithm to an equipment state monitoring system. Is convenient for engineering application.

Description

ITD equipment fault diagnosis method and system
Technical Field
The invention relates to the field of equipment state monitoring and fault diagnosis, in particular to a method and a system for diagnosing faults of an ITD (integrated circuit 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 system for diagnosing faults of an ITD device (the method provided by the invention is abbreviated as MFDFOItd) 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: an ITD equipment fault diagnosis method 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: using eigentime scalesThe solution time-scale Decomposition (ITD) algorithm decomposes the signal x (k) into the sum of n components and a trend term, i.e.
Figure 378425DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the ITD algorithmn(k) Represents the trend term derived by the ITD algorithm, in this example, n = 10;
and step 3: noise components and trend terms are eliminated from the ITD decomposition result by adopting a nonlinear discrimination algorithm, and components c containing fractal features are reservedf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
Figure 305930DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 569421DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 152849DEST_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: FM calculation by Direct integration (DQ)m(k) Is instantaneous angle of, get cf(k) To obtain the instantaneous frequency of cf(k) S of instantaneous dimensionf
And 7: when the dimension is s, then vibrateThe detrending result of the motion signal x (k) is
Figure 467811DEST_PATH_IMAGE006
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 515402DEST_PATH_IMAGE008
Figure 825160DEST_PATH_IMAGE010
step 10: calculating a q-th order function:
Figure 958201DEST_PATH_IMAGE011
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 11608DEST_PATH_IMAGE012
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 ITD algorithm in step 2 includes the following steps:
1) for signal xtT =1,2, …, N, defining an operator
Figure 913705DEST_PATH_IMAGE013
For extracting the low frequency baseline signal, namely:
Figure 456682DEST_PATH_IMAGE014
wherein
Figure 749123DEST_PATH_IMAGE015
Is the baseline signal of the signal from which,
Figure 668537DEST_PATH_IMAGE016
is a natural rotational component, given
Figure 362824DEST_PATH_IMAGE017
Is a real-valued signal that is,
Figure 76702DEST_PATH_IMAGE018
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure 122018DEST_PATH_IMAGE019
. If xtHaving a constant value over a certain interval, we still consider x to be x, taking into account the presence of fluctuations in the adjacent signalstIn this interval, an extreme value is included
Figure 579544DEST_PATH_IMAGE020
Is the right end point of the interval. For convenience, define
Figure 128337DEST_PATH_IMAGE021
Figure 293344DEST_PATH_IMAGE022
Suppose that
Figure 622695DEST_PATH_IMAGE023
And
Figure 821595DEST_PATH_IMAGE024
in that
Figure 287211DEST_PATH_IMAGE025
Above has a definition of xkIn that
Figure 546154DEST_PATH_IMAGE026
Is defined above in the interval
Figure 362800DEST_PATH_IMAGE027
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 162129DEST_PATH_IMAGE013
Namely:
Figure 685514DEST_PATH_IMAGE028
wherein
Figure 177676DEST_PATH_IMAGE029
Where the parameters
Figure 153722DEST_PATH_IMAGE030
Is a linear gain that is a function of the gain,
Figure 756742DEST_PATH_IMAGE031
in this case
Figure 400212DEST_PATH_IMAGE032
2) Defining an intrinsic rotation component extraction operator
Figure 797696DEST_PATH_IMAGE033
Namely:
Figure 591864DEST_PATH_IMAGE034
wherein the content of the first and second substances,
Figure 670678DEST_PATH_IMAGE035
the component obtained by the 1 st iteration decomposition;
Figure 230973DEST_PATH_IMAGE036
baseline signals obtained for the 1 st iterative decomposition; in iteration 1, xtRepresents x (k) in step 2 of claim 1;
3) handle bar
Figure 737040DEST_PATH_IMAGE036
By repeating the above steps as new data, the inherent rotation component whose frequency is successively lowered can be separated
Figure 15575DEST_PATH_IMAGE037
Until the baseline signal becomes monotonic. Thus xtThe whole decomposition process of (a) can be written as:
Figure 898080DEST_PATH_IMAGE038
wherein
Figure 312881DEST_PATH_IMAGE039
Representing the component resulting from the i-th iterative decomposition,
Figure 317746DEST_PATH_IMAGE040
representing the baseline signal resulting from the i-th iterative 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
Figure 755681DEST_PATH_IMAGE042
Figure 238615DEST_PATH_IMAGE044
If e is1And e2All less than 10%, the signal c (k) is discriminated as a noise component or a trend term, and c (k) represents the signal component obtained by the ITD 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 711184DEST_PATH_IMAGE045
Figure 621372DEST_PATH_IMAGE046
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 546602DEST_PATH_IMAGE047
Figure 833227DEST_PATH_IMAGE048
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 160303DEST_PATH_IMAGE049
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 509901DEST_PATH_IMAGE012
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 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 hyperbola, an exponential curve, a pairA number curve or a complex function curve;
2) calculating a least squares metric
Figure 984744DEST_PATH_IMAGE051
3) Let J pair akPartial derivative of
Figure 747164DEST_PATH_IMAGE053
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure 991064DEST_PATH_IMAGE055
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 DQ method in step 6 includes the following steps:
1) let cos (phi (t)) = F (t)), then
Figure 508633DEST_PATH_IMAGE056
Phi (t) stands for FMm(k) T =1,2, …, N, f (t) = FMm(k);
2) Instantaneous angle
Figure 408455DEST_PATH_IMAGE057
Then the instantaneous frequency is
Figure 36883DEST_PATH_IMAGE058
Figure 72972DEST_PATH_IMAGE059
Representing the first derivative of phi (t), the instantaneous scale is defined as
Figure 761442DEST_PATH_IMAGE060
The system based on the ITD equipment fault diagnosis method is characterized in that: the state monitoring system 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 using DQ;
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 system 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 ITD method, noise components and trend terms 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 component estimates the instantaneous frequency of the signal component by using DQ, which can ensure that the instantaneous frequency keeps a positive value and avoid the negative frequency phenomenon;
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 ITD 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 MFDFAoitd 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 MFDFOItd 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 MFDFOItd 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 mfdfaamd based on correlation filtering, mfdfeamed based on correlation filtering, and MFDFAoitd based on correlation filtering according to 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 the five gearbox vibration signals obtained using MFDFOItd 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 the classification of the coordinates of the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained from MFDFAoitd according to the embodiment of the present invention on the states of the five gearboxes, where the symbols "circle", "square", "plus", "diamond" and "left triangle" represent the normal, worn, pitting, broken teeth and worn + pitting gear states, respectively.
Detailed Description
In an embodiment, as shown in fig. 1 and 2, an ITD device fault diagnosis method includes the following steps:
step 1: measuring a device vibration signal x (k) by using an acceleration sensor at a sampling frequency fs, wherein k =1,2, …, N and N are lengths of the sampling signal;
step 2: the signal x (k) is decomposed into the sum of n components and a trend term, i.e., using the Intrinsic time-scale Decomposition (ITD) algorithm
Figure 148561DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the ITD algorithmn(k) Represents the trend term derived by the ITD algorithm, in this example, n = 10;
and step 3: noise components and trend terms are eliminated from the ITD decomposition result by adopting a nonlinear discrimination algorithm, and components c containing fractal features are reservedf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
Figure 315100DEST_PATH_IMAGE061
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 533592DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 395893DEST_PATH_IMAGE004
To obtain cf(k) Frequency modulation ofPart FMm(k),ej(k) Represents cj(k) Envelope of cj(k)=FM(j-1)(k),c1(k)= cf(k);
Step 6: FM calculation by Direct integration (DQ)m(k) Is instantaneous angle of, get cf(k) To obtain the instantaneous frequency of cf(k) S of instantaneous dimensionf
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 4729DEST_PATH_IMAGE062
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 240539DEST_PATH_IMAGE063
Figure 313537DEST_PATH_IMAGE064
step 10: calculating a q-th order function:
Figure 281493DEST_PATH_IMAGE011
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 betweenComprises the following steps: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 502259DEST_PATH_IMAGE012
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 ITD algorithm in the step 2 comprises the following steps:
1) for signal xtT =1,2, …, N, defining an operator
Figure 948283DEST_PATH_IMAGE013
For extracting the low frequency baseline signal, namely:
Figure 141367DEST_PATH_IMAGE065
wherein
Figure 280225DEST_PATH_IMAGE015
Is the baseline signal of the signal from which,
Figure 925970DEST_PATH_IMAGE016
is a natural rotational component, given
Figure 175685DEST_PATH_IMAGE017
Is a real-valued signal that is,
Figure 223276DEST_PATH_IMAGE018
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure 533034DEST_PATH_IMAGE066
. If xtHaving a constant value over a certain interval, we still consider x to be x, taking into account the presence of fluctuations in the adjacent signalstIn this interval, an extreme value is included
Figure 657286DEST_PATH_IMAGE020
Is the right end point of the interval. For convenience, define
Figure 710693DEST_PATH_IMAGE021
Figure 612790DEST_PATH_IMAGE022
Suppose that
Figure 93450DEST_PATH_IMAGE023
And
Figure 448208DEST_PATH_IMAGE024
in that
Figure 305305DEST_PATH_IMAGE025
Above has a definition of xkIn that
Figure 61909DEST_PATH_IMAGE067
Is defined above in the interval
Figure 775787DEST_PATH_IMAGE027
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 821103DEST_PATH_IMAGE013
Namely:
Figure 278630DEST_PATH_IMAGE068
wherein
Figure 827423DEST_PATH_IMAGE029
Where the parameters
Figure 915464DEST_PATH_IMAGE030
Is a linear gain that is a function of the gain,
Figure 244814DEST_PATH_IMAGE031
in this case
Figure 178135DEST_PATH_IMAGE032
2) Defining an intrinsic rotation component extraction operator
Figure 909331DEST_PATH_IMAGE033
Namely:
Figure DEST_PATH_IMAGE069
wherein the content of the first and second substances,
Figure 233521DEST_PATH_IMAGE035
the component obtained by the 1 st iteration decomposition;
Figure 722271DEST_PATH_IMAGE036
baseline signals obtained for the 1 st iterative decomposition; in iteration 1, xtRepresents x (k) in step 2 of claim 1;
3) handle bar
Figure 787179DEST_PATH_IMAGE036
By repeating the above steps as new data, the inherent rotation component whose frequency is successively lowered can be separated
Figure 310564DEST_PATH_IMAGE037
Until the baseline signal becomes monotonic. Thus xtThe whole decomposition process of (a) can be written as:
Figure 537146DEST_PATH_IMAGE070
wherein
Figure DEST_PATH_IMAGE071
Representing the component resulting from the i-th iterative decomposition,
Figure 841088DEST_PATH_IMAGE072
representing the baseline signal resulting from the i-th iterative 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
Figure DEST_PATH_IMAGE073
Figure 444108DEST_PATH_IMAGE074
If e is1And e2Are all less than 10%, the signal c (k) is distinguished as a noise component or a trend term, and c (k) represents the signal component obtained by the ITD 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 822000DEST_PATH_IMAGE045
Figure 422745DEST_PATH_IMAGE046
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 213984DEST_PATH_IMAGE047
Figure 558378DEST_PATH_IMAGE048
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 853093DEST_PATH_IMAGE049
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 359160DEST_PATH_IMAGE012
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 interpolation sequenceThe length of the column is such that,
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 DEST_PATH_IMAGE075
3) Let J pair akPartial derivative of
Figure 637695DEST_PATH_IMAGE076
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure DEST_PATH_IMAGE077
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 DQ method in step 6 includes the steps of:
1) let cos (phi (t)) = F (t)), then
Figure 585447DEST_PATH_IMAGE056
Phi (t) stands for FMm(k) T =1,2, …, N, f (t) = FMm(k);
2) Instantaneous angle
Figure 248DEST_PATH_IMAGE057
Then the instantaneous frequency is
Figure 677217DEST_PATH_IMAGE058
Figure 443047DEST_PATH_IMAGE059
Representing the first derivative of phi (t), the instantaneous scale is defined as
Figure 129244DEST_PATH_IMAGE060
Based on the system of the fault diagnosis method of the ITD equipment, the state monitoring system 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.
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_IMAGE079
The generated multi-fractal simulation signal verifies the performance of the MFDFA, the MFDFAemd and the MFDFOItd. 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 was calculated using MFDFAoitd, and the result is shown in fig. 5. As can be seen from fig. 5, the instantaneous frequencies calculated by MFDFAoitd are all positive frequencies, and thus the MFDFAoitd analysis results are in line with the actual situation. Next, the multi-component of the multi-component simulation signal is calculated by using the MFDFA, the MFDFAemd, and the MFDFAoitd, respectivelyThe spectrum was shaped and the results are 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 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 mfdfaoid and the theoretical value is 0.022, and the average of the relative errors is 5.33%, so that the average of the absolute errors between the multi-fractal spectrum obtained from mfdfaoid is reduced by 65.63%, the average of the relative errors is reduced by 56.35%, the average of the absolute errors between the multi-fractal spectrum obtained from mfdfaoid is reduced by 37.14%, and the average of the relative errors is reduced by 18.87%. 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 noisy multi-fractal simulation signal was calculated using mfdfaamd, MFDFAemd, and MFDFAoitd, respectively, and the result 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 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 MFDFAoitd from the theoretical value is 0.040, and the mean relative error value is 5.56%, so that the multi-fractal spectrum obtained from mfdfaoid is 52.38% smaller than the mean absolute error value of the multi-fractal spectrum obtained from MFDFAemd, and the mean relative error value is 52.92% smaller. It can be seen that MFDFAoitd 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 the mfdfaard, the mfdfeamed based on the correlation filtering, and the MFDFAoitd based on the correlation filtering, respectively, in the embodiment of the present invention. As can be seen from FIG. 10, correlation-based filteringThe analysis result of the MFDFame and the MFDFOItd method based on the related filtering on the noisy multi-fractal simulation signal completely deviates from the theoretical value, so that the related filtering method is easy to damage the fractal structure of the original signal.
Experiment 2 the performance of the algorithm of the invention was verified using gearbox experimental signals.
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 five gearbox vibration signals are analyzed by using an MFDFOItd method, and the multi-fractal spectrums corresponding to the five gearbox vibration signals 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 the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained by the MFDF, MFDFEMd and MFDFOItd methods are respectively extracted to classify the five states of the gearbox, 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 MFDFAoitd method, and therefore the gear box state recognition rate is 100%. It can be seen that the MFDFOItd method can improve the accuracy of the gearbox state identification by 80%.
From the test results, it was assumed after analysis that:
1) the vibration signal is adaptively decomposed by adopting an ITD method, noise components and trend terms 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 component estimates the instantaneous frequency of the signal component by using DQ, which can ensure that the instantaneous frequency keeps a positive value and avoid the negative frequency phenomenon;
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 ITD method is used for automatically determining the type of the signal trend, the continuity of the signal trend is ensured, and the defects of the prior art are effectively overcome;
5) the accuracy and precision of the analysis result are high, and the accuracy of the identification result of the running state of the equipment is high.
It should be appreciated by those skilled in the art that the foregoing embodiments are merely exemplary for better understanding of the present invention, and should not be construed as limiting the scope of the present invention as long as the modifications are made according to the technical solution of the present invention.

Claims (9)

1. An ITD equipment fault diagnosis method 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., using the Intrinsic time-scale Decomposition (ITD) algorithm
Figure 850201DEST_PATH_IMAGE001
Wherein c isi(k) Representing the i-th component, r, obtained by the ITD algorithmn(k) Represents a trend term derived by the ITD algorithm;
and step 3: noise components and trend terms are eliminated from the ITD decomposition result by adopting a nonlinear discrimination algorithm, and components c containing fractal features are reservedf(k) F =1,2, …, p, p represents the number of residual components after filtering;
and 4, step 4: determination of cf(k) Respectively comparing the local maximum value and the local minimum value of c by adopting a cubic spline interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) The upper envelope u (k) and the lower envelope l (k), then cf(k) Is defined as
Figure 767341DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 874975DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure 270184DEST_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: FM calculation by Direct integration (DQ)m(k) Is instantaneous angle of, get cf(k) To obtain the instantaneous frequency of cf(k) S of instantaneous dimensionf
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 146873DEST_PATH_IMAGE006
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 969335DEST_PATH_IMAGE008
Figure 33106DEST_PATH_IMAGE010
step 10: calculating a q-th order function:
Figure 232007DEST_PATH_IMAGE011
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 228781DEST_PATH_IMAGE012
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 ITD device fault diagnosis method of claim 1, wherein: the ITD algorithm in the step 2 comprises the following steps:
1) for signal xtT =1,2, …, N, defining an operator
Figure 222145DEST_PATH_IMAGE013
For extracting the low frequency baseline signal, namely:
Figure 773212DEST_PATH_IMAGE014
wherein
Figure 510224DEST_PATH_IMAGE015
Is the baseline signal of the signal from which,
Figure 364435DEST_PATH_IMAGE016
is a natural rotational component, given
Figure 794279DEST_PATH_IMAGE017
Is a real-valued signal that is,
Figure 35905DEST_PATH_IMAGE018
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure 373345DEST_PATH_IMAGE019
If x istHaving a constant value over a certain interval, we still consider x to be x, taking into account the presence of fluctuations in the adjacent signalstIn this interval, an extreme value is included
Figure 16816DEST_PATH_IMAGE020
Is the right end point of the interval, for convenience, defined
Figure 148720DEST_PATH_IMAGE021
Figure 877642DEST_PATH_IMAGE022
Suppose that
Figure 284353DEST_PATH_IMAGE023
And
Figure 782330DEST_PATH_IMAGE024
in that
Figure 350715DEST_PATH_IMAGE025
Above has a definition of xkIn that
Figure 301353DEST_PATH_IMAGE026
Is defined above in the interval
Figure 246175DEST_PATH_IMAGE027
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 598659DEST_PATH_IMAGE013
Namely:
Figure 869104DEST_PATH_IMAGE028
wherein
Figure 307038DEST_PATH_IMAGE029
Where the parameters
Figure 524393DEST_PATH_IMAGE030
Is a linear gain that is a function of the gain,
Figure 996962DEST_PATH_IMAGE031
in this case
Figure 265457DEST_PATH_IMAGE032
2) Defining an intrinsic rotation component extraction operator
Figure 190687DEST_PATH_IMAGE033
Namely:
Figure 414995DEST_PATH_IMAGE034
wherein the content of the first and second substances,
Figure 273230DEST_PATH_IMAGE035
the component obtained by the 1 st iteration decomposition;
Figure 619897DEST_PATH_IMAGE036
baseline signals obtained for the 1 st iterative decomposition; in iteration 1, xtRepresentative claims1 in step 2 x (k);
3) handle bar
Figure 32424DEST_PATH_IMAGE037
By repeating the above steps as new data, the inherent rotation component whose frequency is successively lowered can be separated
Figure 388319DEST_PATH_IMAGE038
Until the baseline signal becomes monotonic, such that xtThe whole decomposition process of (a) can be written as:
Figure 569902DEST_PATH_IMAGE039
wherein
Figure 556312DEST_PATH_IMAGE040
Representing the component resulting from the i-th iterative decomposition,
Figure 456135DEST_PATH_IMAGE041
representing the baseline signal resulting from the i-th iterative decomposition.
3. The ITD device fault diagnosis method of 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 ofBy Hsurr(q) represents;
3) two parameters e are defined1And e2
Figure 84563DEST_PATH_IMAGE043
Figure 386231DEST_PATH_IMAGE045
If e is1And e2All less than 10%, the signal c (k) is discriminated as a noise component or a trend term, and c (k) represents the signal component obtained by the ITD algorithm.
4. The ITD device fault diagnosis method of 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 ITD device fault diagnosis method of 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 ITD device fault diagnosis method of 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):
Figure 809122DEST_PATH_IMAGE046
Figure 196241DEST_PATH_IMAGE047
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 628360DEST_PATH_IMAGE048
Figure 518955DEST_PATH_IMAGE049
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 4) Calculate the firstqAverage of the order fluctuation function:
Figure 646836DEST_PATH_IMAGE050
5) if it is notx(k) Presence of self-similar features thenqMean value of order fluctuation functionF q (s) And time scalesStore betweenIn the power law relationship:
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 255672DEST_PATH_IMAGE012
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 ITD device fault diagnosis method of 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
Figure 429164DEST_PATH_IMAGE052
3) Let J pair akPartial derivative of
Figure 236583DEST_PATH_IMAGE053
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX, X=(x(1),x(2),…,x(n)) T
Figure DEST_PATH_IMAGE055
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 ITD device fault diagnosis method of claim 1, wherein: the DQ method in step 6 includes the steps of:
1) let cos (phi (t)) = F (t)), then
Figure 266856DEST_PATH_IMAGE056
Phi (t) stands for FMm(k) T =1,2, …, N, f (t) = FMm(k);
2) Instantaneous angle
Figure DEST_PATH_IMAGE057
Then the instantaneous frequency is
Figure 425305DEST_PATH_IMAGE058
Figure DEST_PATH_IMAGE059
Representing the first derivative of phi (t), the instantaneous scale is defined as
Figure DEST_PATH_IMAGE061
9. A system for implementing a method of fault diagnosis of an ITD device according to any of claims 1 to 8, characterized in that: the system 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.
CN202011238927.0A 2020-11-09 2020-11-09 ITD equipment fault diagnosis method and system Withdrawn CN112697475A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011238927.0A CN112697475A (en) 2020-11-09 2020-11-09 ITD equipment fault diagnosis method and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011238927.0A CN112697475A (en) 2020-11-09 2020-11-09 ITD equipment fault diagnosis method and system

Publications (1)

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

Family

ID=75506092

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011238927.0A Withdrawn CN112697475A (en) 2020-11-09 2020-11-09 ITD equipment fault diagnosis method and system

Country Status (1)

Country Link
CN (1) CN112697475A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117909621A (en) * 2024-03-19 2024-04-19 青蛙泵业股份有限公司 Method and system for monitoring running state of sewage and dirt submerged electric pump

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117909621A (en) * 2024-03-19 2024-04-19 青蛙泵业股份有限公司 Method and system for monitoring running state of sewage and dirt submerged electric pump
CN117909621B (en) * 2024-03-19 2024-05-24 青蛙泵业股份有限公司 Method and system for monitoring running state of sewage and dirt submerged electric pump

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
CN113375939B (en) Mechanical part fault diagnosis method based on SVD and VMD
CN111238813B (en) Method for extracting fault features of rolling bearing under strong interference
CN112683393A (en) LMD equipment fault diagnosis method and system
CN112697266A (en) ITD and GZC machine state monitoring method and device
CN112697484A (en) SSD multi-scale fluctuation analysis 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
CN112697475A (en) ITD equipment fault diagnosis method and system
CN112697477A (en) EMD equipment fault diagnosis method and system
CN112697470A (en) SSD device fault diagnosis method and system
CN112697479A (en) ITD multi-scale fluctuation analysis state monitoring method and device
CN112683395A (en) Method and device for monitoring states of ELMD and GZC machines
CN112881046A (en) VMD 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
CN112697473A (en) EMD and GZC machine 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
CN106153333B (en) A kind of envelope Analysis Method based on wavelet decomposition filtering
CN112683394A (en) EEMD equipment fault diagnosis method and system
CN112697478A (en) WD equipment fault diagnosis method and system
CN112697480A (en) NMD equipment fault diagnosis method and system
CN112683390A (en) EITD equipment fault diagnosis method and system

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