CN110514444B - Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis - Google Patents

Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis Download PDF

Info

Publication number
CN110514444B
CN110514444B CN201910504140.5A CN201910504140A CN110514444B CN 110514444 B CN110514444 B CN 110514444B CN 201910504140 A CN201910504140 A CN 201910504140A CN 110514444 B CN110514444 B CN 110514444B
Authority
CN
China
Prior art keywords
phase space
component
signal
kurtosis
matrix
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.)
Active
Application number
CN201910504140.5A
Other languages
Chinese (zh)
Other versions
CN110514444A (en
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.)
Southeast University
Original Assignee
Southeast University
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 Southeast University filed Critical Southeast University
Priority to CN201910504140.5A priority Critical patent/CN110514444B/en
Publication of CN110514444A publication Critical patent/CN110514444A/en
Application granted granted Critical
Publication of CN110514444B publication Critical patent/CN110514444B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00—Testing of machine parts
    • G01M13/04—Bearings
    • G01M13/045—Acoustic or vibration analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention relates to a rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis, which comprises the following steps of: acquiring a vibration signal of a rolling bearing by using a vibration sensor to obtain a weak fault signal x (t); setting a punishment factor alpha and a decomposition layer number K, and carrying out variation mode decomposition on the vibration signal to obtain a corresponding intrinsic mode function component; selecting the intrinsic mode function component with the highest kurtosis as the optimal component according to a kurtosis criterion; performing phase space parallel factor analysis on the optimal component to obtain at least one independent component in a high-dimensional phase space; selecting the independent component with the highest kurtosis from the at least one independent component obtained in the fourth step to carry out envelope spectrum analysis, and obtaining fault characteristic frequency; the method effectively extracts the fault characteristic frequency and realizes accurate identification of the fault type of the rolling bearing.

Description

Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis
Technical Field
The invention relates to a rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis, and belongs to the technical field of fault diagnosis of mechanical equipment.
Background
Rolling bearings are one of the most common and critical components of rotating mechanical equipment, the operating conditions of which directly affect the working performance of the whole equipment. By incomplete statistics, about 30% of mechanical failures are caused by bearing failure. This means that local defect detection of rolling bearings is essential for the field of practical engineering; therefore, how to identify the rolling bearing failure as early as possible and prevent the occurrence of accidents is an important issue of current research.
In the early failure stage of the rolling bearing, the periodic pulse failure signal is susceptible to the tiny size of local defects and strong background noise, and the failure characteristic information of the periodic pulse failure signal is difficult to accurately detect by using a traditional signal processing algorithm. Some experts and scholars propose a lot of effective weak fault diagnosis methods. However, there are some limitations, such as the Stochastic Resonance (SR) algorithm requires more parameters; the robustness of the minimum entropy deconvolution is poor; empirical mode decomposition has the problems of mode aliasing and the like.
Disclosure of Invention
The invention provides a rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis, which can effectively extract fault feature frequency and realize accurate identification of the fault type of a rolling bearing.
The technical scheme adopted by the invention for solving the technical problems is as follows:
a rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis comprises the following steps:
the first step is as follows: acquiring a vibration signal of a rolling bearing by using a vibration sensor to obtain a weak fault signal x (t);
the second step is that: setting a punishment factor alpha and a decomposition layer number K, and carrying out variation mode decomposition on the vibration signal to obtain a corresponding intrinsic mode function component;
the third step: selecting the intrinsic mode function component with the highest kurtosis as the optimal component according to a kurtosis criterion;
the fourth step: performing phase space parallel factor analysis on the optimal component to obtain at least one independent component in a high-dimensional phase space;
the fifth step: selecting the independent component with the highest kurtosis from the at least one independent component obtained in the fourth step to carry out envelope spectrum analysis, and obtaining fault characteristic frequency;
as a further preferred aspect of the present invention, the variation modal decomposition in the second step can be regarded as a constraint variation optimization problem:
Figure GDA0002238500470000021
the above is the expression of the constraint variation model of the obtained original weak fault signal x, where | | · |. luminance |2Is a 2 norm, delta (t) is a pulse function,
Figure GDA0002238500470000022
for the gradient of the function, j is the unit of an imaginary number, let uk(t) is the kth eigenmode function, let wkIs ukThe center frequency of (d); budgeting a secondary penalty factor alpha and a Lagrange multiplier lambda, and changing the constraint variable problem into an unconstrained variable problem:
Figure GDA0002238500470000023
solving the above formula by using an alternating direction multiplier, and updating each eigenmode function component u for K equal to 1,2, …, K, wherein K is the maximum value of all eigenmode functions, K and all frequencies ω ≧ 0k(t) and its center frequency ωk：
Figure GDA0002238500470000024
Figure GDA0002238500470000025
Wherein u isiTo remove ukThe components of other intrinsic mode functions, alpha is a secondary penalty factor, x (t) is a weak fault signal, omegakIs the center frequency; for all frequencies ω ≧ 0, update λ:
Figure GDA0002238500470000026
wherein, lambda is Lagrange multiplication operator, theta is updated parameter of lambda, ukIs the k-th intrinsic mode function; the above process is repeatedly updated until the following convergence condition is satisfied:
Figure GDA0002238500470000031
wherein epsilon is generally set to 10-6Thereby obtaining a corresponding eigenmode function component uk(t)，k＝1,2,...,K；
As a further preferred embodiment of the present invention, the kurtosis criterion in the third step is
Figure GDA0002238500470000032
Wherein
Figure GDA0002238500470000033
T is the time length of the signal,
Figure GDA0002238500470000034
is ukMean of (t), s.t.: subject to indicates that the condition u is satisfiedkFor the function of the k-th intrinsic mode,
Figure GDA0002238500470000035
is the optimal component;
as a further preferred aspect of the present invention, the phase space parallelism analysis described in the fourth step comprises the main steps of
Let y be the optimal component selected in the third step, whose higher dimensional phase space can be written as
Figure GDA0002238500470000036
Wherein the content of the first and second substances,
Figure GDA0002238500470000037
for the reconstructed phase space signal, L ═ N- (M-1) τ is the number of points in the phase space, N is the signal length, M is the embedding dimension, τ is the delay time, where M and τ values are predetermined, other parameters are known, and the high-dimensional phase space signal Y can be expressed as Y ═ a · S
Wherein the content of the first and second substances,
Figure GDA0002238500470000038
is a mixing matrix, S is an independent component; dividing the high-dimensional phase space signal Y into D non-overlapping data blocks, wherein each data block contains P ═ L/D data points, and L is the point number of the phase space; thus, the high-dimensional phase space signal Y can also be expressed as
Y＝[Y1 Y2 … YD]＝A·[S1 S2 … SD]
Wherein
Figure GDA0002238500470000041
M high-dimensional phase space signals for the d data block;
Figure GDA0002238500470000042
m independent components for the d data block; y isdCan be expressed as
Figure GDA0002238500470000043
Wherein the content of the first and second substances,
Figure GDA0002238500470000044
Figure GDA0002238500470000045
is SdThe covariance matrix of (a);
to simplify the mathematical notation of the subsequent derivation, the notation is annotated as follows:
Figure GDA0002238500470000046
wherein the content of the first and second substances,
Figure GDA0002238500470000047
for m in the d data block1A high-dimensional phase space signal, m1＝1,2,…,M；
Figure GDA0002238500470000048
For m in the d data block2A high-dimensional phase space signal, m1＝1,2,…,M； m 21,2, …, M; r is a covariance matrix;
Figure GDA0002238500470000049
wherein the content of the first and second substances,
Figure GDA00022385004700000410
for m in the d data block1An independent component, m1＝1，2，…，M；
Figure GDA00022385004700000411
For m in the d data block2An independent component, m1＝1,2,…,M； m 21,2, …, M; according to
Figure GDA00022385004700000412
m1≠m2Is obtained by
Figure GDA00022385004700000413
Thus, the device is provided with
Figure GDA00022385004700000414
Can be expressed as
Figure GDA0002238500470000051
Wherein, ammIs an element in the matrix a, M is 1, 2. Defining a matrix
Figure GDA0002238500470000052
The expression is
Figure GDA0002238500470000053
Each covariance matrix is divided into
Figure GDA0002238500470000054
D is superimposed as a third order tensor form
Figure GDA0002238500470000055
The elements of which can be represented as
Figure GDA0002238500470000056
Wherein the content of the first and second substances,
Figure GDA0002238500470000057
is an element in the tensor R, ammFor the elements in the matrix a, it is,
Figure GDA0002238500470000058
is a matrix RsThe upper formula of the element in (1) is a PARAFAC model; sectioning the third order tensor R along the z-axisOperating to obtain an equivalent two-dimensional matrix RMD×M：
Figure GDA0002238500470000059
Defining the following objective function
Figure GDA00022385004700000510
Where off (-) is the sum of the squares of the matrix off-diagonal elements and W is the separation matrix; performing minimum optimization on the objective function by using Givens rotation, and estimating a separation matrix W, so as to obtain an independent component S which is W.Y, wherein Y is a reconstructed phase space signal;
for the selection of the delay time and the embedding dimension, wherein M represents the embedding dimension, and tau represents the delay time, in order to avoid the blindness of artificial selection, the global kurtosis maximum criterion of independent components is provided; the basic steps of the criterion are as follows: 1) setting the variation range of the delay time and the embedding dimension; 2) under different parameters, obtaining corresponding independent components by utilizing a phase space parallel factor analysis algorithm, and calculating corresponding kurtosis values; 3) obtaining a relation between the two parameters and the maximum kurtosis value of the independent component; 4) and parameters corresponding to the global maximum kurtosis value are the optimal delay time and the embedding dimension. Its cost function can be expressed as
Figure GDA0002238500470000061
Wherein, kurt(s)m) M is 1,2 …, M is the M-th independent component smThe kurtosis value of (a) is,
Figure GDA0002238500470000062
is the maximum value of the embedding dimension;
Figure GDA0002238500470000063
is the maximum value of the delay time, N is the signal length, M is the embeddingDimension, MoFor the optimal embedding dimension, τoTo an optimal delay time.
Through the technical scheme, compared with the prior art, the invention has the following beneficial effects:
1. the invention decomposes the fault signal into a plurality of intrinsic mode components by using variation mode decomposition, and selects the component with the largest kurtosis for further analysis by using a kurtosis rule, and aims to eliminate the noise interference of irrelevant characteristic frequency bands;
2. performing phase space parallel factor analysis on the optimal component, and selecting an independent component with the highest global kurtosis, wherein the purpose is to suppress noise interference embedded in a characteristic frequency band;
3. envelope spectrum analysis is carried out on the independent component with the maximum global kurtosis, and fault characteristic frequency is effectively extracted, so that the fault type of the rolling bearing is accurately identified;
4. in order to avoid blindness of artificial parameter selection, the method provides a criterion of maximum global kurtosis of the independent component, and selects the optimal embedding dimension and delay time parameter.
Drawings
The invention is further illustrated with reference to the following figures and examples.
FIG. 1 is a weak fault feature extraction flow diagram of the preferred embodiment of the present invention;
fig. 2 is various data of weak fault signals of a rolling bearing collected according to a preferred embodiment of the present invention, where 2a is a time domain waveform, 2b is an FFT spectrum, and 2c is an envelope spectrum;
FIG. 3 is the data of four eigenmode function components obtained by eigenmode decomposition of the vibration signal according to the preferred embodiment of the present invention, where 3a is the time domain waveform and 3b is the FFT spectrum;
FIG. 4 is a kurtosis value of an eigenmode-function component of a preferred embodiment of the present invention.
FIG. 5 is a relationship between the maximum kurtosis value of an independent component and two parameters (delay time, embedding dimension) according to the preferred embodiment of the present invention;
FIG. 6 is a diagram of the independent components in the high dimensional phase space of the preferred embodiment of the present invention;
FIG. 7 is a kurtosis value corresponding to an independent component in accordance with a preferred embodiment of the present invention;
fig. 8 is an envelope spectrum obtained by envelope demodulation analysis of the independent component with the maximum kurtosis according to the preferred embodiment of the present invention, where 8a is the envelope spectrum analysis and 8b is the obtained fault feature frequency.
Detailed Description
The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic views illustrating only the basic structure of the present invention in a schematic manner, and thus show only the constitution related to the present invention.
As shown in fig. 1, the rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis of the present invention includes the following steps:
the first step is as follows: acquiring a vibration signal of a rolling bearing by using a vibration sensor to obtain a weak fault signal x (t);
the second step is that: setting a punishment factor alpha and a decomposition layer number K, and carrying out variation mode decomposition on the vibration signal to obtain a corresponding intrinsic mode function component;
the third step: selecting the intrinsic mode function component with the highest kurtosis as the optimal component according to a kurtosis criterion;
the fourth step: performing phase space parallel factor analysis on the optimal component to obtain at least one independent component in a high-dimensional phase space;
the fifth step: selecting the independent component with the highest kurtosis from the at least one independent component obtained in the fourth step to carry out envelope spectrum analysis, and obtaining fault characteristic frequency;
specifically, the method comprises the following steps of acquiring a rolling bearing vibration signal by using a vibration sensor to obtain a weak fault signal x (t), wherein various data of the weak fault signal are shown in fig. 2, wherein 2a is a time domain waveform, 2b is an FFT spectrum, and 2c is an envelope spectrum;
example 1:
setting a punishment factor alpha as 4000 and a decomposition layer number K as 4, and carrying out variation mode decomposition on the vibration signal to obtain a series of intrinsic mode function components;
the variation modal decomposition can be regarded as a constraint variation optimization problem:
Figure GDA0002238500470000081
the above is the expression of the constraint variation model of the obtained original weak fault signal x, where | | · |. luminance |2Is a 2 norm, delta (t) is a pulse function,
Figure GDA0002238500470000082
for the gradient of the function, j is the unit of an imaginary number, let uk(t) is the kth eigenmode function, let wkIs ukThe center frequency of (d); budgeting a secondary penalty factor alpha and a Lagrange multiplier lambda, and changing the constraint variable problem into an unconstrained variable problem:
Figure GDA0002238500470000083
solving the above formula by using an alternating direction multiplier, and updating each eigenmode function component u for K equal to 1,2, …, K, wherein K is the maximum value of all eigenmode functions, K and all frequencies ω ≧ 0k(t) and its center frequency ωk：
Figure GDA0002238500470000084
Figure GDA0002238500470000085
Wherein u isiTo remove ukThe components of other intrinsic mode functions, alpha is a secondary penalty factor, x (t) is a weak fault signal, omegakIs the center frequency; for all frequencies ω ≧ 0, update λ:
Figure GDA0002238500470000086
wherein, lambda is Lagrange multiplication operator, theta is updated parameter of lambda, ukIs the k-th intrinsic mode function; the above process is repeatedly updated until the following convergence condition is satisfied:
Figure GDA0002238500470000091
wherein epsilon is generally set to 10-6The vibration signal is decomposed into four eigenmode components through variation mode, the time domain waveform of the vibration signal is shown as 3a in fig. 3, and the FFT spectrum is shown as 3b in fig. 3.
Example 2:
FIG. 4 shows that the kurtosis value of the eigenmode is calculated according to the kurtosis criterion
Figure GDA0002238500470000092
Wherein
Figure GDA0002238500470000093
T is the time length of the signal,
Figure GDA0002238500470000094
is ukMean of (t), s.t.: subject to indicates that the condition u is satisfiedkSelecting the intrinsic mode component 2 with the maximum kurtosis as the optimal component for further analysis;
example 3:
the optimal component is subjected to phase space parallel factor analysis to obtain a plurality of independent components in a high phase space, and the method mainly comprises the following steps
Let y be the optimal component selected in the third step, whose higher dimensional phase space can be written as
Figure GDA0002238500470000095
Wherein the content of the first and second substances,
Figure GDA0002238500470000096
for the reconstructed phase space signal, L ═ N- (M-1) τ is the number of points in the phase space, N is the signal length, M is the embedding dimension, τ is the delay time, where M and τ values are predetermined, other parameters are known, and the high-dimensional phase space signal Y can be expressed as Y ═ a · S
Wherein the content of the first and second substances,
Figure GDA0002238500470000101
is a mixing matrix, S is an independent component; dividing the high-dimensional phase space signal Y into D non-overlapping data blocks, wherein each data block contains P ═ L/D data points, and L is the point number of the phase space; thus, the high-dimensional phase space signal Y can also be expressed as
Y＝[Y1 Y2 … YD]＝A·[S1 S2 … SD]
Wherein
Figure GDA0002238500470000102
M high-dimensional phase space signals for the d data block;
Figure GDA0002238500470000103
m independent components for the d data block; y isdCan be expressed as
Figure GDA0002238500470000104
Wherein the content of the first and second substances,
Figure GDA0002238500470000105
Figure GDA0002238500470000106
is SdThe covariance matrix of (a);
to simplify the mathematical notation of the subsequent derivation, the notation is annotated as follows:
Figure GDA0002238500470000107
wherein the content of the first and second substances,
Figure GDA0002238500470000108
for m in the d data block1A high-dimensional phase space signal, m1＝1,2,…,M；
Figure GDA0002238500470000109
For m in the d data block2A high-dimensional phase space signal, m1＝1,2,…,M； m 21,2, …, M; r is a covariance matrix;
Figure GDA00022385004700001010
wherein the content of the first and second substances,
Figure GDA00022385004700001011
for m in the d data block1An independent component, m1＝1，2，…，M；
Figure GDA00022385004700001012
For m in the d data block2An independent component, m1＝1,2,…,M； m 21,2, …, M; according to
Figure GDA0002238500470000111
m1≠m2Is obtained by
Figure GDA0002238500470000112
Thus, the device is provided with
Figure GDA0002238500470000113
Can be expressed as
Figure GDA0002238500470000114
Wherein amm is an element in the matrix a, and M is 1, 2. Defining a matrix
Figure GDA0002238500470000115
The expression is
Figure GDA0002238500470000116
Each covariance matrix is divided into
Figure GDA0002238500470000117
D is superimposed as a third order tensor form
Figure GDA0002238500470000118
The elements of which can be represented as
Figure GDA0002238500470000119
Wherein the content of the first and second substances,
Figure GDA00022385004700001110
being the elements in tensor R, amm being the elements in matrix a,
Figure GDA00022385004700001111
is a matrix RsThe upper formula of the element in (1) is a PARAFAC model; performing section operation on the third-order tensor R along the z axis to obtain an equivalent two-dimensional matrix RMD×M：
Figure GDA0002238500470000121
Defining the following objective function
Figure GDA0002238500470000122
Where off (-) is the sum of the squares of the matrix off-diagonal elements and W is the separation matrix; performing minimum optimization on the objective function by using Givens rotation, and estimating a separation matrix W, so as to obtain an independent component S which is W.Y, wherein Y is a reconstructed phase space signal;
example 4:
for the selection of the delay time and the embedding dimension, M represents the embedding dimension, the variation range of the embedding dimension is set to be [1,200], T represents the delay time, the variation range of the delay time is set to be [1,11], under different parameters, a phase space parallel factor analysis algorithm is utilized to obtain corresponding independent components, and corresponding kurtosis values are calculated; thus, a relationship between the two parameters and the maximum kurtosis values of the independent components is obtained, as shown in FIG. 5; as can be seen from fig. 5, the optimal delay time and embedding dimension for the global maximum kurtosis value are 10 and 163, respectively, fig. 6 is an independent separation in the high-dimensional phase space under the optimal parameters, and the corresponding kurtosis value is calculated, as shown in fig. 7;
example 5:
as shown in 8a of fig. 8, the independent component with the highest kurtosis is selected for envelope spectrum analysis, and the fault characteristic frequency is obtained, and as shown in 8b of fig. 8, the fault characteristic frequency and 2 multiples thereof can be seen.
It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the prior art and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.
The meaning of "and/or" as used herein is intended to include both the individual components or both.
The term "connected" as used herein may mean either a direct connection between components or an indirect connection between components via other components.
In light of the foregoing description of the preferred embodiment of the present invention, many modifications and variations will be apparent to those skilled in the art without departing from the spirit and scope of the invention. The technical scope of the present invention is not limited to the content of the specification, and must be determined according to the scope of the claims.

Claims (3)

1. A rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis is characterized by comprising the following steps: the method comprises the following steps:
the first step is as follows: acquiring a vibration signal of a rolling bearing by using a vibration sensor to obtain a weak fault signal x (t);
the second step is that: setting a punishment factor alpha and a decomposition layer number K, and carrying out variation mode decomposition on the vibration signal to obtain a corresponding intrinsic mode function component;
the third step: selecting the intrinsic mode function component with the highest kurtosis as the optimal component according to a kurtosis criterion;
the fourth step: performing phase space parallel factor analysis on the optimal component to obtain at least one independent component in a high-dimensional phase space;
the fifth step: selecting the independent component with the highest kurtosis from the at least one independent component obtained in the fourth step to carry out envelope spectrum analysis, and obtaining fault characteristic frequency;
the fourth step is the analysis of the phase space parallelism factor, which comprises the main steps of
Let y be the optimal component selected in the third step, whose higher dimensional phase space can be written as
Figure FDA0002834250700000011
Wherein the content of the first and second substances,
Figure FDA0002834250700000012
for the reconstructed phase space signal, L ═ N- (M-1) τ is the number of points in the phase space, N is the signal length, M is the embedding dimension, τ is the delay time, where M, τ values can be preset, other parameters are known, and the high-dimensional phase space signal Y can be expressed as
Y＝A·S
Wherein the content of the first and second substances,
Figure FDA0002834250700000013
is a mixing matrix, S is an independent component; dividing the high-dimensional phase space signal Y into D non-overlapping data blocks, wherein each data block contains P ═ L/D data points, and L is the point number of the phase space; thus, the high-dimensional phase space signal Y can also be expressed as
Y＝[Y1 Y2…YD]＝A·[S1 S2…SD]
Wherein
Figure FDA0002834250700000021
M high-dimensional phase space signals for the d data block;
Figure FDA0002834250700000022
m independent components for the d data block; y isdCan be expressed as
Figure FDA0002834250700000023
Wherein the content of the first and second substances,
Figure FDA0002834250700000024
Figure FDA0002834250700000025
is SdThe covariance matrix of (a);
to simplify the mathematical notation of the subsequent derivation, the notation is annotated as follows:
Figure FDA0002834250700000026
wherein the content of the first and second substances,
Figure FDA0002834250700000027
for m in the d data block1A high-dimensional phase space signal, m1＝1,2,…,M；
Figure FDA0002834250700000028
For m in the d data block2A high-dimensional phase space signal, m1＝1,2,…,M；m21,2, …, M; r is a covariance matrix;
Figure FDA0002834250700000029
wherein the content of the first and second substances,
Figure FDA00028342507000000210
for m in the d data block1An independent component, m1＝1,2,…,M；
Figure FDA00028342507000000211
For m in the d data block2An independent component, m1＝1,2,…,M；m21,2, …, M; according to
Figure FDA00028342507000000212
m1≠m2Is obtained by
Figure FDA00028342507000000213
Thus, the device is provided with
Figure FDA00028342507000000214
Can be expressed as
Figure FDA0002834250700000031
Wherein, ammIs an element in the matrix A, M is 1,2, …, M; defining a matrix
Figure FDA0002834250700000032
The expression is
Figure FDA0002834250700000033
Each covariance matrix is divided into
Figure FDA0002834250700000034
Superimposed as third order tensor form
Figure FDA0002834250700000035
The elements of which can be represented as
Figure FDA0002834250700000036
Wherein the content of the first and second substances,
Figure FDA0002834250700000037
is an element in the tensor R, ammFor the elements in the matrix a, it is,
Figure FDA0002834250700000038
is a matrix RsThe upper formula of the element in (1) is a PARAFAC model; performing section operation on the third-order tensor R along the z axis to obtain an equivalent two-dimensional matrix RMD×M：
Figure FDA0002834250700000039
Defining the following objective function
Figure FDA0002834250700000041
Where off (-) is the sum of the squares of the matrix off-diagonal elements and W is the separation matrix; performing minimum optimization on the objective function by using Givens rotation, and estimating a separation matrix W, so as to obtain an independent component S which is W.Y, wherein Y is a reconstructed phase space signal;
for the selection of the delay time and the embedding dimension, wherein M represents the embedding dimension, and tau represents the delay time, in order to avoid the blindness of artificial selection, the global kurtosis maximum criterion of independent components is provided; the basic steps of the criterion are as follows: 1) setting the variation range of the delay time and the embedding dimension; 2) under different parameters, obtaining corresponding independent components by utilizing a phase space parallel factor analysis algorithm, and calculating corresponding kurtosis values; 3) obtaining a relation between the two parameters and the maximum kurtosis value of the independent component; 4) the parameters corresponding to the global maximum kurtosis value are the optimal delay time and the embedding dimension; its cost function can be expressed as
Figure FDA0002834250700000042
Wherein, kurt(s)m) M is 1,2 …, M is the M-th independent component smThe kurtosis value of (a) is,
Figure FDA0002834250700000043
is the maximum value of the embedding dimension;
Figure FDA0002834250700000044
is the maximum value of the delay time, N is the signal length and M is the embedding dimension.
2. The rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis according to claim 1, characterized in that: the variation modal decomposition in the second step can be regarded as a constraint variation optimization problem:
Figure FDA0002834250700000051
the above is the expression of the constraint variation model of the obtained original weak fault signal x, where | | · |. luminance |2Is a 2 norm, delta (t) is a pulse function,
Figure FDA0002834250700000052
for the gradient of the function, j is the unit of an imaginary number, let uk(t) is the kth eigenmode function, let wkIs ukThe center frequency of (d); budgeting a secondary penalty factor alpha and a Lagrange multiplier lambda, and changing the constraint variable problem into an unconstrained variable problem:
Figure FDA0002834250700000053
solving the above formula by using an alternating direction multiplier, and updating each eigenmode function component u for K equal to 1,2, …, K, wherein K is the maximum value of all eigenmode functions, K and all frequencies ω ≧ 0k(t) and its center frequency ωk：
Figure FDA0002834250700000054
Figure FDA0002834250700000055
Wherein u isiTo remove ukThe components of other intrinsic mode functions, alpha is a secondary penalty factor, x (t) is a weak fault signal, omegakIs the center frequency; for all frequencies ω ≧ 0, update λ:
Figure FDA0002834250700000056
wherein, lambda is Lagrange multiplication operator, theta is updated parameter of lambda, ukIs the k-th intrinsic mode function; the above process is repeatedly updated until the following convergence condition is satisfied:
Figure FDA0002834250700000061
wherein epsilon is generally set to 10-6Thereby obtaining a corresponding eigenmode function component uk(t)，k＝1,2,…,K。
3. The method for extracting weak fault features of a rolling bearing according to claim 1, wherein: the kurtosis criterion in the third step is
Figure FDA0002834250700000062
Wherein
Figure FDA0002834250700000063
T is the time length of the signal,
Figure FDA0002834250700000064
is ukMean of (t), s.t.: subject to indicates that the condition u is satisfiedkFor the function of the k-th intrinsic mode,
Figure FDA0002834250700000065
is the optimal component.
CN201910504140.5A 2019-06-12 2019-06-12 Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis Active CN110514444B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910504140.5A CN110514444B (en) 2019-06-12 2019-06-12 Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910504140.5A CN110514444B (en) 2019-06-12 2019-06-12 Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis

Publications (2)

Publication Number Publication Date
CN110514444A CN110514444A (en) 2019-11-29
CN110514444B true CN110514444B (en) 2021-04-06

Family

ID=68622413

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910504140.5A Active CN110514444B (en) 2019-06-12 2019-06-12 Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis

Country Status (1)

Country Link
CN (1) CN110514444B (en)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110916716A (en) * 2019-12-30 2020-03-27 龙岩学院 Wearable heart sound monitoring facilities
CN111178318B (en) * 2020-01-06 2023-07-11 东南大学 Rolling bearing early compound fault feature extraction method based on progressive VMD
CN112113766B (en) * 2020-09-01 2021-11-09 兰州理工大学 Characteristic extraction method for early damage state of rolling bearing
CN112431753B (en) * 2021-01-25 2021-04-16 赛腾机电科技（常州）有限公司 Multiple quantitative diagnosis method for shoe loosening fault of axial plunger pump
CN113204051B (en) * 2021-06-10 2022-04-15 成都理工大学 Low-rank tensor seismic data denoising method based on variational modal decomposition

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109187024A (en) * 2018-09-04 2019-01-11 温州大学激光与光电智能制造研究院 A kind of piston type air compressor crankcase Fault Diagnosis of Roller Bearings
CN109632310A (en) * 2019-01-18 2019-04-16 北京化工大学 A kind of Method for Bearing Fault Diagnosis based on feature enhancing
CN109781412A (en) * 2019-02-26 2019-05-21 东南大学 A kind of rolling bearing adaptive resonance demodulation method based on EEMD

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109187024A (en) * 2018-09-04 2019-01-11 温州大学激光与光电智能制造研究院 A kind of piston type air compressor crankcase Fault Diagnosis of Roller Bearings
CN109632310A (en) * 2019-01-18 2019-04-16 北京化工大学 A kind of Method for Bearing Fault Diagnosis based on feature enhancing
CN109781412A (en) * 2019-02-26 2019-05-21 东南大学 A kind of rolling bearing adaptive resonance demodulation method based on EEMD

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
"Decentralized modal identification of structures using parallel factor decomposition and sparse blind source separation";A. Sadhu等;《Mechanical Systems and Signal Processing》;20130930;第41卷;第396-419页 *
"基于VMD和谱峭度的滚动轴承早期故障诊断方法";唐贵基等;《中国测试》;20170930;第43卷(第9期);第112-117页 *
"基于变分模态分解和排列熵的滚动轴承故障诊断";郑小霞等;《振动与冲击》;20171130;第36卷(第22期);第22-28页 *
"基于相空间独立分量分析及峭度贡献系数的早期故障分析方法研究";陈建国等;《振动与冲击》;20160630;第35卷(第12期);第155-159页 *
"基于相空间重构的独立分量分析及其工程应用";赵长生等;《轴承》;20130131(第1期);第51-54页 *

Also Published As

Publication number Publication date
CN110514444A (en) 2019-11-29

Similar Documents

Publication Publication Date Title
CN110514444B (en) Rolling bearing weak fault feature extraction method based on variational modal decomposition and phase space parallel factor analysis
CN109829402B (en) GS-SVM-based bearing damage degree diagnosis method under different working conditions
CN110303380A (en) A kind of cutting tool for CNC machine method for predicting residual useful life
CN109632963B (en) Structural damage four-dimensional imaging method constructed based on time-invariant characteristic signals
CN112686096A (en) Rolling bearing fault diagnosis method based on multi-scale diffusion entropy and VPMCD
CN114296075A (en) Ground penetrating radar image artificial intelligence identification method and device
CN109918417B (en) Time sequence data self-adaptive segmentation, dimension reduction and characterization method based on wavelet transformation and application
CN107358588A (en) Phase goes to roll up pleat method, MRI scan method and MR imaging apparatus
Saravanakumar et al. Hierarchical symbolic analysis and particle swarm optimization based fault diagnosis model for rotating machineries with deep neural networks
CN104779144A (en) Signature analytics for improving lithographic process of manufacturing semiconductor devices
Shi et al. The VMD-scale space based hoyergram and its application in rolling bearing fault diagnosis
CN110426191A (en) A kind of method for diagnosing faults of anti-interference rotating machinery
CN111007390B (en) Analog circuit fault diagnosis model based on algebraic method
CN112733612A (en) Cross-domain rotating machinery fault diagnosis model establishing method and application thereof
CN108985357A (en) The hyperspectral image classification method of set empirical mode decomposition based on characteristics of image
CN109782158B (en) Analog circuit diagnosis method based on multi-stage classification
Safeiee et al. Optimization using Mahalanobis-Taguchi System for inductor component
CN111079591B (en) Bad data restoration method and system based on improved multi-scale principal component analysis
CN110335287A (en) The extracting method and device of Architectural drawing data
CN105809177A (en) Method used for actuating remote sensing image classification
CN107389367B (en) A kind of signal adaptive solution enveloping method and calculating equipment based on optimal signal-to-noise ratio
CN108664936A (en) A kind of diagnostic method and system based on mechanical disorder
CN115014748A (en) Fault diagnosis method for seed cotton sorting spray valve
Wu et al. Hilbert ID considering multi-window feature extraction for transformer deep vision fault positioning
CN109272054B (en) Vibration signal denoising method and system based on independence

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
GR01 Patent grant
GR01 Patent grant