CN112697266A - ITD and GZC machine state monitoring method and device - Google Patents

ITD and GZC machine state monitoring method and device Download PDF

Info

Publication number
CN112697266A
CN112697266A CN202011240613.4A CN202011240613A CN112697266A CN 112697266 A CN112697266 A CN 112697266A CN 202011240613 A CN202011240613 A CN 202011240613A CN 112697266 A CN112697266 A CN 112697266A
Authority
CN
China
Prior art keywords
signal
data
itd
points
curve
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
CN202011240613.4A
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 CN202011240613.4A priority Critical patent/CN112697266A/en
Publication of CN112697266A publication Critical patent/CN112697266A/en
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses an ITD and GZC machine state monitoring method and a device, which decompose a vibration signal of equipment by using an ITD algorithm, remove a noise component and a trend term by using a nonlinear discrimination algorithm, reserve a fractal signal component, interpolate extreme points by using a piecewise linear interpolation function, fit an envelope by using a least square method, separate a frequency modulation part, estimate an instantaneous frequency by using a GZC algorithm and calculate a corresponding instantaneous scale, determine a vibration signal detrending result according to an analysis scale, calculate a multi-fractal spectrum of a detrended signal, extract left end points, right end points and extreme point coordinates of the multi-fractal spectrum as characteristic parameters of the running state of the equipment, identify the running state of the equipment, deploy the algorithm to the equipment state monitoring device, can accurately distinguish the running state of the equipment, and the equipment state monitoring device has good flexibility and portability, is convenient for engineering application.

Description

ITD and GZC machine state monitoring method and device
Technical Field
The invention relates to the field of equipment state monitoring and fault diagnosis, in particular to a method and a device for monitoring states of an ITD (integrated circuit) machine and a GZC (GZC) machine.
Background
The device vibration signal contains rich fractal features that can describe the operating state of the device. The box dimension, power spectrum analysis and re-standard range method can estimate the single-fractal parameters of stationary signals, and the de-trend fluctuation analysis (DFA) can estimate the single-fractal dimension of non-stationary signals. However, when the device fails, the vibration signal is usually non-stationary and has a multi-fractal characteristic, and the conventional fractal dimension estimation method generates a relatively large error. The multi-fractal detrending fluctuation analysis (MFDF) can estimate multi-fractal parameters of non-stationary signals, but the MFDF method has the problems that the analysis scale needs to be manually determined, the fitting polynomial trend order is difficult to determine, and the data segment is discontinuous. Currently, there is a document that proposes an MFDFA version (MFDFAemd) based on EMD to solve the problem of MFDFA. However, the linear filtering method adopted by mfdfame is easy to destroy the fractal structure of the original signal, and there is a negative frequency phenomenon, and these defects seriously affect the application effect of mfdfame. In summary, in the prior art, it is difficult to accurately extract the multi-fractal features of the device vibration signal, and it is difficult to accurately detect the device operating state.
Disclosure of Invention
The invention provides a method and a device for monitoring the states of an ITD (integrated circuit) machine and a GZC machine (the method is abbreviated as MFDFOItd). 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 and GZC machine state monitoring 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 594204DEST_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) The local maximum value and the local minimum value of c are respectively compared by adopting a piecewise linear 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 DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure DEST_PATH_IMAGE003
m times, j =1,2, …, m, until
Figure DEST_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 was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 42503DEST_PATH_IMAGE005
；
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 DEST_PATH_IMAGE006
Figure 166448DEST_PATH_IMAGE007
；
step 10: calculating a q-th order function:
Figure DEST_PATH_IMAGE008
；
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 DEST_PATH_IMAGE009
；
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(α) = q (α -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 xt， t=1, 2, …,N，xtRepresenting the signal x (k) in step 2 of claim 1, defining an operator
Figure DEST_PATH_IMAGE010
For extracting the low frequency baseline signal, namely:
Figure 669979DEST_PATH_IMAGE011
wherein
Figure DEST_PATH_IMAGE012
Is the baseline signal of the signal from which,
Figure 218772DEST_PATH_IMAGE013
is a natural rotational component, given
Figure DEST_PATH_IMAGE014
Is a real-valued signal that is,
Figure 837973DEST_PATH_IMAGE015
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure DEST_PATH_IMAGE016
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 714793DEST_PATH_IMAGE017
Is the right end point of the interval, for convenience, defined
Figure DEST_PATH_IMAGE018
，
Figure 976010DEST_PATH_IMAGE019
Suppose that
Figure DEST_PATH_IMAGE020
And
Figure 685034DEST_PATH_IMAGE021
in that
Figure DEST_PATH_IMAGE022
Above has a definition of xkIn that
Figure 412819DEST_PATH_IMAGE023
Is defined above in the interval
Figure DEST_PATH_IMAGE024
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 963886DEST_PATH_IMAGE010
Namely:
Figure 966477DEST_PATH_IMAGE025
wherein
Figure DEST_PATH_IMAGE026
Where the parameters
Figure 37332DEST_PATH_IMAGE027
Is a linear gain that is a function of the gain,
Figure DEST_PATH_IMAGE028
in this case
Figure 529494DEST_PATH_IMAGE029
；
2) Defining an intrinsic rotation component extraction operator
Figure DEST_PATH_IMAGE030
Namely:
Figure 613862DEST_PATH_IMAGE031
wherein the content of the first and second substances,
Figure DEST_PATH_IMAGE032
the component obtained by the 1 st iteration decomposition;
Figure 623406DEST_PATH_IMAGE033
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 860353DEST_PATH_IMAGE033
By repeating the above steps as new data, the inherent rotation component whose frequency is successively lowered can be separated
Figure DEST_PATH_IMAGE034
Until the baseline signal becomes monotonic, such that xtThe whole decomposition process of (a) can be written as:
Figure 742989DEST_PATH_IMAGE035
wherein
Figure DEST_PATH_IMAGE036
Representing the component resulting from the i-th iterative decomposition,
Figure 206332DEST_PATH_IMAGE037
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 DEST_PATH_IMAGE038
，
Figure 81884DEST_PATH_IMAGE039
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 discrete Fourier inversion on frequency domain data subjected to phase substitutionGet 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 DEST_PATH_IMAGE040
Figure 688183DEST_PATH_IMAGE041
；
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 DEST_PATH_IMAGE042
Figure 928672DEST_PATH_IMAGE043
；
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 DEST_PATH_IMAGE044
；
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 676048DEST_PATH_IMAGE009
；
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c ，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
Further, the least square method in the step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) A sequence or pair c generated by interpolating the local maxima off(k) The local minima of the sequence, n represents the length of the interpolated sequence,
1) a set of functions r is selected in advancek(t)，k=1,2,…,m，m<n, constructioning function
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t) in which rk(t) 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 902761DEST_PATH_IMAGE045
；
3) Let J pair akPartial derivative of
Figure DEST_PATH_IMAGE046
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 989666DEST_PATH_IMAGE047
Wherein R isTRepresenting the transposed matrix of R, R-1An inverse matrix representing R;
4）rk(t) respectively selecting a second-order polynomial, a third-order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve and a composite function curve for calculation, then comparing least square indexes J generated by various curve forms, and selecting a curve form r corresponding to the minimum J as the curve formk(t) form (a).
Further, the GZC method in step 6 includes the following steps:
1) determining local maximum, minimum and zero points of the signal c (k), k =1,2, …, N, c (k) = FMm(k)；
2) Counting the time intervals T between two continuous local maximum points, two continuous local minimum points, two continuous rising zeros and two continuous falling zeros according to the time sequence1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and descending zero points, adjacent descending zero points and local minimum points, adjacent local minimum points and ascending zero points, and adjacent ascending zero points and local maximum points according to the time sequence3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
Figure DEST_PATH_IMAGE048
The temporal scale of the signal c (k) is sf=T。
Based on the ITD and GZC machine state monitoring method, the device for realizing the method is characterized in that: the state monitoring device in the step 16 comprises the following parts: the system comprises a data line, an acceleration sensor, a data acquisition card, a case, a notebook computer and signal analysis software, wherein the acceleration sensor is connected with the data acquisition card through the data line, the data acquisition card is installed in the case, the case is connected with the notebook computer through the data line, the signal analysis software is installed on the notebook computer, and the signal analysis software is used for realizing the algorithm.
The method comprises the following steps:
step 1), the following steps: collecting a vibration signal;
step 2), the step of: decomposing an original signal into different component sum forms, wherein some components correspond to noise and trend terms, and some components contain fractal features;
and 3, step 3: removing noise components and trend items in the signal decomposition result by using a nonlinear discrimination algorithm, and only reserving signal components containing fractal features;
4) to 6): separating the frequency modulation part of each fractal signal component, and estimating the instantaneous frequency and the instantaneous scale of each fractal signal component by utilizing the GZC;
step 7), the steps of: selecting proper fractal signal components according to the analysis scale, and summing the selected fractal signal components to be used as a signal de-trend result corresponding to the analysis scale;
8) to 14): performing fluctuation analysis on the signal detrending result corresponding to each analysis scale to obtain a multi-fractal spectrum of the original signal;
step 15): extracting the vertical coordinates of the left end point, the right end point and the extreme point of the multi-fractal spectrum, and taking the three parameters as the characteristic parameters of the running state of the equipment;
16) step: the algorithm is deployed on an equipment state monitoring device to monitor the equipment state.
By adopting the technical scheme, compared with the prior art, the invention has the following advantages:
1) the vibration signal is adaptively decomposed by adopting an 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) a frequency modulation part for separating the signal components, and estimating the instantaneous frequency of the signal components by using the GZC, so that the instantaneous frequency can be ensured to keep a positive value, and the negative frequency phenomenon is avoided;
3) calculating corresponding instantaneous scale according to the instantaneous frequency of the signal component, and performing fluctuation analysis according to the instantaneous scale of the signal component, so that the defect of manually setting the scale is avoided;
4) the 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,
the number of signal components is 10;
figure 5 is a graph showing the instantaneous frequency of a multi-fractal simulation signal obtained by the MFDFAoitd method according to an embodiment of the present invention,
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 ITD 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 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 using MFDFOItd in an embodiment of the present invention;
FIG. 15 is a diagram showing the classification results of the left end point, the right end point and the extreme point coordinates of the multi-fractal spectrum obtained from MFDFA on the four gear box states in the example of the present invention, where the symbols "circle", "square", "plus" and "diamond" represent the normal, light scratch, heavy scratch and broken tooth gear states, respectively;
fig. 16 is a classification result of coordinates of a left end point, a right end point, and an extreme point of a multi-fractal spectrum obtained from MFDFAemd on 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 results of classifying the four gear box states by the coordinates of the left end point, the right end point and the extreme point of the multi-fractal spectrum obtained from MFDFAoitd according to 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
In an embodiment, as shown in fig. 1 and fig. 2, an ITD and GZC machine state monitoring 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 525689DEST_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) The local maximum value and the local minimum value of c are respectively compared by adopting a piecewise linear 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 432465DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 430246DEST_PATH_IMAGE049
m times, j =1,2, …, m, until
Figure 168395DEST_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 was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 750686DEST_PATH_IMAGE005
；
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 3813DEST_PATH_IMAGE006
Figure 493700DEST_PATH_IMAGE007
；
step 10: calculating a q-th order function:
Figure DEST_PATH_IMAGE050
；
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 sizesHave power betweenLaw relation: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 102667DEST_PATH_IMAGE009
；
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 xt， t=1, 2, …,N，xtRepresenting the signal x (k) in step 2 of claim 1, defining an operator
Figure 652597DEST_PATH_IMAGE010
For extracting the low frequency baseline signal, namely:
Figure 533965DEST_PATH_IMAGE011
wherein
Figure 889860DEST_PATH_IMAGE012
Is the baseline signal of the signal from which,
Figure 337022DEST_PATH_IMAGE013
is a natural rotational component, given
Figure 995537DEST_PATH_IMAGE014
Is a real-valued signal that is,
Figure 478383DEST_PATH_IMAGE015
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure 513335DEST_PATH_IMAGE016
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 80582DEST_PATH_IMAGE017
Is the right end point of the interval, for convenience, defined
Figure 769053DEST_PATH_IMAGE018
，
Figure 625013DEST_PATH_IMAGE019
Suppose that
Figure 260394DEST_PATH_IMAGE020
And
Figure 495197DEST_PATH_IMAGE021
in that
Figure 495514DEST_PATH_IMAGE022
Above has a definition of xkIn that
Figure 697825DEST_PATH_IMAGE023
Is defined above in the interval
Figure 136897DEST_PATH_IMAGE024
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 616420DEST_PATH_IMAGE010
Namely:
Figure 161540DEST_PATH_IMAGE025
wherein
Figure 726513DEST_PATH_IMAGE026
Where the parameters
Figure 703696DEST_PATH_IMAGE027
Is a linear gain that is a function of the gain,
Figure 896780DEST_PATH_IMAGE028
in this case
Figure 504479DEST_PATH_IMAGE051
；
2) Defining an intrinsic rotation component extraction operator
Figure 353487DEST_PATH_IMAGE030
Namely:
Figure 681831DEST_PATH_IMAGE031
wherein the content of the first and second substances,
Figure 932684DEST_PATH_IMAGE032
the component obtained by the 1 st iteration decomposition;
Figure 711284DEST_PATH_IMAGE033
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 109904DEST_PATH_IMAGE033
As new dataBy repeating the above steps, the inherent rotation component with the sequentially decreased frequency can be separated
Figure 428890DEST_PATH_IMAGE034
Until the baseline signal becomes monotonic, such that xtThe whole decomposition process of (a) can be written as:
Figure 471932DEST_PATH_IMAGE035
wherein
Figure 218171DEST_PATH_IMAGE036
Representing the component resulting from the i-th iterative decomposition,
Figure 353355DEST_PATH_IMAGE037
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 413715DEST_PATH_IMAGE038
，
Figure 373581DEST_PATH_IMAGE039
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.
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 618618DEST_PATH_IMAGE040
Figure 867196DEST_PATH_IMAGE041
；
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 527985DEST_PATH_IMAGE042
Figure 155406DEST_PATH_IMAGE043
；
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 712290DEST_PATH_IMAGE044
；
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 244902DEST_PATH_IMAGE009
；
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c ，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
The least square method in the step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) Is interpolated to produce a sequence of local maximaOr to cf(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 771698DEST_PATH_IMAGE045
；
3) Let J pair akPartial derivative of
Figure 643840DEST_PATH_IMAGE046
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 168362DEST_PATH_IMAGE047
Wherein R isTRepresenting the transposed matrix of R, R-1An inverse matrix representing R;
4）rk(t) respectively selecting a second-order polynomial, a third-order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve and a composite function curve for calculation, then comparing least square indexes J generated by various curve forms, and selecting a curve form r corresponding to the minimum J as the curve formk(t) form (a).
The GZC method in the step 6 comprises the following steps:
1) determining local maximum, minimum and zero points of the signal c (k), k =1,2, …, N, c (k) = FMm(k)；
2) Counting the time intervals T between two continuous local maximum points, two continuous local minimum points, two continuous rising zeros and two continuous falling zeros according to the time sequence1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and descending zero points, adjacent descending zero points and local minimum points, adjacent local minimum points and ascending zero points, and adjacent ascending zero points and local maximum points according to the time sequence3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
Figure 499855DEST_PATH_IMAGE048
The temporal scale of the signal c (k) is sf=T。
Based on the above-mentioned ITD and GZC machine status monitoring method, the device for implementing the method, wherein the status monitoring device in 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_IMAGE052
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 the view of figure 4,the instantaneous frequency calculated by the MFDFame has many negative frequencies, and the MFDFame analysis result has large errors 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, a multi-fractal spectrum of the multi-fractal simulation signal is calculated using mfdfaard, mfdfeaitd, and MFDFAoitd, respectively, and the result is shown in fig. 6. According to the results shown in fig. 6, it is found through calculation that the average of the absolute errors between the 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.027, and the average of the relative errors is 5.24%, so that the average of the absolute errors between the multi-fractal spectrum obtained from mfdfaoid is reduced by 57.81%, the average of the relative errors is reduced by 57.08%, the average of the absolute errors between the multi-fractal spectrum obtained from mfdfaoid is reduced by 22.86%, and the average of the relative errors is reduced by 20.24%. 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 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 mfdfaoidtd and the theoretical value is 0.031, and the mean relative error value is 5.47%, so that the multi-fractal spectrum obtained from mfdfaoidtd is reduced by 63.10% and the mean relative error value is reduced by 53.68%. It can be seen that MFDF and MFDMFDFOItd has better noise immunity than FAemd. 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, the analysis result of the MFDFAemd based on the correlation filtering and the MFDFAoitd based on the correlation filtering on the noisy multi-fractal simulation signal completely deviates from the theoretical value, so that the correlation filtering method easily destroys the fractal structure of the original signal.
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 MFDFOItd 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 MFDFOItd 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 gearbox 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 gearbox state recognition rate is 100%. It can be seen that the MFDFOItd 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 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 components estimates the instantaneous scale of the signal components by utilizing the GZC, so that the instantaneous frequency can be ensured to keep a positive value, and the negative frequency phenomenon is avoided;
3) the fluctuation analysis is carried out according to the instantaneous scale of the signal component, so that the defect of manually setting the scale is avoided;
4) the 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 and GZC machine state monitoring 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 598912DEST_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) The local maximum value and the local minimum value of c are respectively compared by adopting a piecewise linear interpolation functionf(k) The local maximum and local minimum are interpolated, and c is fitted by least square methodf(k) Upper part ofEnvelope u (k) and lower envelope l (k), then cf(k) Is defined as
Figure 139615DEST_PATH_IMAGE002
The symbol | x | represents taking the absolute value of x;
and 5: repeatedly executing formula
Figure 232479DEST_PATH_IMAGE002
m times, j =1,2, …, m, until
Figure 55936DEST_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 was calculated using the Generalized zero-crossing method (GZC)m(k) Average period of (2) to obtain cf(k) S of instantaneous dimensionf；
And 7: when the scale is s, the detrending result of the vibration signal x (k) is
Figure 722540DEST_PATH_IMAGE005
；
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 801355DEST_PATH_IMAGE006
Figure 237015DEST_PATH_IMAGE007
；
step 10: calculating a q-th order function:
Figure 24986DEST_PATH_IMAGE004
；
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 395519DEST_PATH_IMAGE009
；
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 and GZC machine state monitoring method according to claim 1, wherein: the ITD algorithm in the step 2 comprises the following steps:
1) for signal xt， t=1, 2, …,N，xtRepresenting the signal x (k) in step 2 of claim 1, defining an operator
Figure 278024DEST_PATH_IMAGE010
For extracting the low frequency baseline signal, namely:
Figure 568191DEST_PATH_IMAGE011
wherein
Figure 510739DEST_PATH_IMAGE012
Is the baseline signal of the signal from which,
Figure 886357DEST_PATH_IMAGE013
is a natural rotational component, given
Figure 306974DEST_PATH_IMAGE014
Is a real-valued signal that is,
Figure 717227DEST_PATH_IMAGE015
represents xtThe time corresponding to the local extremum of (c) is defined for convenience
Figure 565097DEST_PATH_IMAGE016
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 926546DEST_PATH_IMAGE017
Is the right end point of the interval, for convenience, defined
Figure 150854DEST_PATH_IMAGE018
，
Figure 415613DEST_PATH_IMAGE019
Suppose that
Figure 699964DEST_PATH_IMAGE020
And
Figure 50174DEST_PATH_IMAGE021
in that
Figure 812594DEST_PATH_IMAGE022
Above has a definition of xkIn that
Figure 931859DEST_PATH_IMAGE023
Is defined above in the interval
Figure 387112DEST_PATH_IMAGE024
Defining a piecewise linear baseline signal extraction operator between successive extreme points of (A) to (B)
Figure 723153DEST_PATH_IMAGE010
Namely:
Figure 289263DEST_PATH_IMAGE025
wherein
Figure 263035DEST_PATH_IMAGE026
Where the parameters
Figure 889189DEST_PATH_IMAGE027
Is a linear gain that is a function of the gain,
Figure 213991DEST_PATH_IMAGE028
；
2) defining an intrinsic rotation component extraction operator
Figure 318213DEST_PATH_IMAGE029
Namely:
Figure 412071DEST_PATH_IMAGE030
wherein the content of the first and second substances,
Figure 209126DEST_PATH_IMAGE031
the component obtained by the 1 st iteration decomposition;
Figure 260039DEST_PATH_IMAGE032
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 433532DEST_PATH_IMAGE032
By repeating the above steps as new data, the inherent rotation component whose frequency is successively lowered can be separated
Figure 381896DEST_PATH_IMAGE033
Until the baseline signal becomes monotonic, such that xtThe whole decomposition process of (a) can be written as:
Figure 349852DEST_PATH_IMAGE034
wherein
Figure 383667DEST_PATH_IMAGE035
Representing the component resulting from the i-th iterative decomposition,
Figure 829692DEST_PATH_IMAGE036
representing the baseline signal resulting from the i-th iterative decomposition.
3. The ITD and GZC machine state monitoring method according to claim 1, wherein: the step 3 of the nonlinear discriminant algorithm comprises the following steps:
1) performing rearrangement operation and substitution operation on signals c (k), and using c as data obtained by rearrangement operationshuf(k) Indicating that data obtained after the substitution operation are csurr(k) Represents;
2) for c (k), cshuf(k) And csurr(k) Performing 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 614843DEST_PATH_IMAGE006
，
Figure 10052DEST_PATH_IMAGE008
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 and GZC machine state monitoring method according to claim 3, wherein: the data rearrangement operation in the step 1) comprises the following steps: randomly randomizing the order of the components c (k).
5. The ITD and GZC machine state monitoring method according to claim 3, wherein: the data replacement operation in the step 1) comprises the following steps:
1) performing a discrete fourier transform on component c (k) to obtain the phase of component c (k);
2) replacing the original phase of component c (k) with a set of pseudo-independent identically distributed numbers located in the (-pi, pi) interval;
3) performing inverse discrete Fourier transform on the frequency domain data subjected to phase substitution to obtain data cIFFT(k) To obtain data cIFFT(k) The real part of (a).
6. The ITD and GZC machine state monitoring method according to claim 3, wherein: the MFDF method in the step 2) comprises the following steps:
1) contour of construction x (k), k =1,2, …, NY(i)：
Figure 56646DEST_PATH_IMAGE039
Figure 306362DEST_PATH_IMAGE040
；
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 291635DEST_PATH_IMAGE041
Figure 539077DEST_PATH_IMAGE042
；
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 609801DEST_PATH_IMAGE043
；
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 600891DEST_PATH_IMAGE009
；
6) taking logarithm of both sides of the formula in step 5) to obtain ln [ 2 ]F q (s)]=H(q)ln(s)+c ，cIs constant, whereby the slope of the straight line can be obtainedH(q)。
7. The ITD and GZC machine state monitoring method according to claim 1, wherein: the least square method in the step 4 comprises the following steps: for x (t), t =1,2, …, n, x (t) represents the pair c in step 4f(k) A sequence or pair c generated by interpolating the local maxima off(k) The local minima of the sequence, n represents the length of the interpolated sequence,
1) a set of functions r is selected in advancek(t)，k=1,2,…,m，m<n, constructioning function
f(t)=a1r1(t)+a2r2(t)+…+ amrm(t) in which rk(t) represents a second order polynomial, a third order polynomial, a hyperbolic curve, an exponential curve, a logarithmic curve or a complex function curve;
2) calculating a least squares metric
Figure 440671DEST_PATH_IMAGE044
；
3) Let J pair akPartial derivative of
Figure 859014DEST_PATH_IMAGE045
K =1,2, …, m, when (a)1,a2,…,am)T =(RTR)-1RTX， X=(x(1),x(2),…,x(n)) T，
Figure 151455DEST_PATH_IMAGE046
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 and GZC machine state monitoring method according to claim 1, wherein: the GZC method in the step 6 comprises the following steps:
1) determining local maximum, minimum and zero points of the signal c (k), k =1,2, …, N, c (k) = FMm(k)；
2) According to the time sequence, two continuous local maximum value points, two continuous local minimum value points, two continuous rising zero points and two continuous falling points are countedTime interval T between zero points1j，j=1,2,…,N1，N1Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and local minimum points, adjacent local minimum points and local maximum points, adjacent rising zero points and falling zero points, and adjacent falling zero points and rising zero points according to the time sequence2j，j=1,2,…,N2，N2Representing the number of time intervals, and counting the time intervals T between adjacent local maximum points and descending zero points, adjacent descending zero points and local minimum points, adjacent local minimum points and ascending zero points, and adjacent ascending zero points and local maximum points according to the time sequence3j，j=1,2,…,N3，N3Calculating the average period T of the signal c (k) representing the number of time intervals
Figure 67832DEST_PATH_IMAGE010
The temporal scale of the signal c (k) is sf=T。
9. Apparatus for implementing an ITD and GZC machine condition monitoring method according to any of claims 1 to 8, characterized in that: the device comprises the following parts: the system comprises a data line, an acceleration sensor, a data acquisition card, a case, a notebook computer and signal analysis software, wherein the acceleration sensor is connected with the data acquisition card through the data line, the data acquisition card is installed in the case, the case is connected with the notebook computer through the data line, the signal analysis software is installed on the notebook computer, and the signal analysis software is used for realizing the algorithm.
CN202011240613.4A 2020-11-09 2020-11-09 ITD and GZC machine state monitoring method and device Withdrawn CN112697266A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011240613.4A CN112697266A (en) 2020-11-09 2020-11-09 ITD and GZC machine state monitoring method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011240613.4A CN112697266A (en) 2020-11-09 2020-11-09 ITD and GZC machine state monitoring method and device

Publications (1)

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

Family

ID=75506414

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011240613.4A Withdrawn CN112697266A (en) 2020-11-09 2020-11-09 ITD and GZC machine state monitoring method and device

Country Status (1)

Country Link
CN (1) CN112697266A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114969631A (en) * 2022-05-26 2022-08-30 Oppo广东移动通信有限公司 Baseband chip, channel estimation method, data processing method and equipment
CN116992393A (en) * 2023-09-27 2023-11-03 联通(江苏)产业互联网有限公司 Safety production monitoring method based on industrial Internet of things

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114969631A (en) * 2022-05-26 2022-08-30 Oppo广东移动通信有限公司 Baseband chip, channel estimation method, data processing method and equipment
CN114969631B (en) * 2022-05-26 2024-05-10 Oppo广东移动通信有限公司 Baseband chip, channel estimation method, data processing method and equipment
CN116992393A (en) * 2023-09-27 2023-11-03 联通(江苏)产业互联网有限公司 Safety production monitoring method based on industrial Internet of things
CN116992393B (en) * 2023-09-27 2023-12-22 联通(江苏)产业互联网有限公司 Safety production monitoring method based on industrial Internet of things

Similar Documents

Publication Publication Date Title
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
CN108845250B (en) On-load tap-changer fault identification method based on vibration signal feature extraction
CN108801630B (en) Gear fault diagnosis method for single-channel blind source separation
CN112697266A (en) ITD and GZC machine state monitoring method and device
CN110991055B (en) Residual life prediction system for rotary mechanical products
CN112697484A (en) SSD multi-scale fluctuation analysis state monitoring method and device
CN112697483A (en) LMD multi-scale fluctuation analysis state monitoring method and device
CN113375939A (en) Mechanical part fault diagnosis method based on SVD and VMD
CN112683393A (en) LMD equipment fault diagnosis method and system
CN110717472B (en) Fault diagnosis method and system based on improved wavelet threshold denoising
CN116304570B (en) Hydraulic turbine fault signal denoising method based on EEMD combined Chebyshev filtering
Zhang et al. Improved local cepstrum and its applications for gearbox and rolling bearing fault detection
CN112697479A (en) ITD multi-scale fluctuation analysis state monitoring method and device
CN112697473A (en) EMD and GZC machine state monitoring method and device
CN112697472A (en) WD multi-scale fluctuation analysis state monitoring method and device
CN112697475A (en) ITD equipment fault diagnosis method and system
CN112697482A (en) VMD 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
CN112697470A (en) SSD device fault diagnosis method and system
CN112697477A (en) EMD equipment fault diagnosis method and system
CN112697263A (en) EEMD multi-scale fluctuation analysis state monitoring method and device
CN112597958B (en) Automatic identification method and system for rolling bearing fault
CN113435281A (en) Ripple compensator fault diagnosis method based on improved HHT conversion
CN112697469A (en) LMD 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