A bearing fault diagnosis method, system, device and medium

Abstract

The present invention discloses a bearing fault diagnosis method, system, equipment and medium, including obtaining a vibration acceleration signal of a bearing with strong interference and weak impact characteristics, processing it, and constructing a candidate embedding dimension sequence; starting the search from the starting dimension of the candidate embedding dimension sequence to obtain a symplectic geometric component, and repeating the search until all the candidate embedding dimension sequences are searched, obtaining the optimal embedding dimension corresponding to the minimum distribution consistency index, and performing symplectic geometric modal decomposition, calculating the improved sideband power value; taking the symplectic geometric component exceeding the upper threshold as a fault sensitive feature, performing envelope demodulation analysis, and determining the fault type. Symplectic geometric modal decomposition based on Wasserstein distance not only has powerful signal decomposition energy, but also can retain the essential characteristics of the vibration signal during the decomposition process, and while achieving fault feature decoupling and weak feature extraction, ensures that the decomposed information can fully reflect the bearing operation status.

Description

Technical Field

[0001] The invention belongs to the field of composite fault decoupling and weak fault feature extraction, and specifically relates to a bearing fault diagnosis method, system, equipment and medium. Background Art

[0002] With the growth of industrial production capacity, high speed, precision, high efficiency and intelligence have gradually become the development trend of mechanical equipment. Ensuring that these equipment can work efficiently and reliably for a long time and reducing the number of accidental downtime can bring huge economic benefits and competitive advantages to production companies. Rotating parts (such as bearings, gears, etc.) are important components of mechanical equipment, and their performance directly determines whether the mechanical equipment can operate safely and smoothly for a long time.

[0003] Major national equipment, represented by rocket engines, is often a complex nonlinear system composed of multiple bearings with strong coupling of mechanical, combustion, flow and other processes, and operates under extreme conditions of high pressure, high speed, flammability and low temperature. Therefore, the collected vibration signals often contain time-varying components with different amplitudes caused by multi-source vibrations such as unstable combustion, fluid excitation, and turbine pump vibration, which makes the vibration analysis of the equipment extremely difficult. Therefore, it is of great value to carry out fault feature decoupling and weak fault feature extraction technology under strong noise background.

[0004] In recent years, methods represented by variational mode decomposition (VMD) and its variants have been increasingly widely used in fault feature decoupling and feature extraction. However, these methods require the determination of the number of decompositions K and the quadratic penalty factor α based on the analyst's expert experience, and the selection of parameters has a significant impact on the decomposition effect. The parameters determined by beginners often have poor decomposition accuracy and are difficult to effectively achieve feature extraction. Therefore, for signals with low signal-to-noise ratio and coupled fault features, how to achieve fault feature decoupling and extract effective weak fault features to achieve bearing fault diagnosis is a problem to be solved. Summary of the invention

[0005] In order to overcome the problem that the feature extraction methods in the prior art are difficult to achieve feature decoupling and weak feature extraction from low signal-to-noise ratio and feature coupling signals, the purpose of the present invention is to propose a bearing fault diagnosis method, system, equipment and medium, which can separate a single fault feature from the feature coupling signal caused by multi-source vibration, and achieve fault feature decoupling and weak feature extraction.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A bearing fault diagnosis method comprises the following steps:

[0008] Obtain vibration acceleration signals with strong interference and weak impact characteristics of different bearing fault states;

[0009] The vibration acceleration signal of the bearing during its entire life cycle is transformed into a power spectrum by Fourier transform, and the initial value of the embedding dimension is determined based on the power spectrum, and a candidate embedding dimension sequence is constructed.

[0010] Starting the search from the starting dimension of the candidate embedding dimension sequence, obtaining the symplectic geometric component and the initial single signal component, and repeating the search until the entire candidate embedding dimension sequence is searched, obtaining the optimal embedding dimension corresponding to the minimum distribution consistency index;

[0011] According to the optimal embedding dimension corresponding to the minimum distribution consistency index, the symplectic geometric mode decomposition based on Wasserstein distance is performed to obtain the optimal symplectic geometric component;

[0012] According to the optimal symplectic geometric component, an improved sideband power value of the symplectic geometric component is calculated;

[0013] The symplectic geometric components whose improved sideband power values ​​exceed the upper threshold are taken as fault sensitive features, and envelope demodulation analysis is performed on the fault sensitive features to determine the fault type.

[0014] Furthermore, the vibration acceleration signal of the bearing throughout its life cycle is Fourier transformed to obtain a power spectrum, and the initial value of the embedding dimension is determined according to the power spectrum, and a candidate embedding dimension sequence is constructed, including the following steps:

[0015] The bearing vibration signal is transformed by Fourier transform to obtain the power spectrum, and the fundamental frequency and rotation frequency F of the signal are determined according to the power spectrum. r , according to f∈[F r -Δf,F r +Δf] with the maximum power spectral density is used as the initial embedding dimension, and a candidate embedding dimension sequence is constructed according to the initial embedding dimension, where f is the frequency.

[0016] Furthermore, the frequency f and the initial embedding dimension d p The relationship is:

[0017]

[0018] Among them, F s is the sampling frequency, f max For f∈[F r -Δf,F r +Δf] has the frequency corresponding to the maximum power spectral density, Δf is used to constrain the frequency band range, and n is the signal length;

[0019] Candidate embedding dimension sequence D w for:

[0020] D w =[dp -delta,…,d p +delta]

[0021] delta=max(15,c×d p )

[0022] Among them, delta is the parameter and c is the weight coefficient.

[0023] Furthermore, the symplectic geometric component and the initial single signal component are obtained by the following process:

[0024] The initial value of the embedding dimension is taken as a variable parameter to construct the trajectory matrix with the set embedding dimension. The trajectory matrix is ​​decomposed using the symplectic geometric mode decomposition method based on Wasserstein distance to obtain the symplectic geometric components and the initial single signal components.

[0025] Furthermore, the symplectic geometry components are obtained by the following process:

[0026] According to Takens embedding theorem, the original analysis signal sig=x1,x2,…,x n , n is the signal length, which is expanded into a trajectory matrix using topological methods;

[0027] Perform autocorrelation analysis on the trajectory matrix X to obtain the covariance matrix; and establish the Hamilton matrix based on the covariance matrix; and use the Hamilton matrix to construct a symplectic geometry matrix that meets the requirements of symplectic geometry analysis:

[0028] The symplectic geometry matrix that meets the requirements of symplectic geometry analysis is solved by QR to obtain the eigenvector matrix, and the initial single-component matrix is ​​constructed according to the eigenvector matrix and the trajectory matrix;

[0029] The multi-dimensional initial single component matrix is ​​converted into a one-dimensional vector by using diagonal averaging to obtain the initial single signal component;

[0030] Assuming that d initial single signal components are obtained by decomposition, the first initial single signal component is taken as the benchmark, and the Wasserstein distance between it and the remaining d-1 initial single signal components is evaluated. The initial single signal components whose Wasserstein distance exceeds the threshold are reorganized to obtain the first symplectic geometric component; then, the initial single signal component with the largest energy proportion is selected from the remaining initial single signal components as the benchmark, and the reorganization is continuously performed to obtain the symplectic geometric component.

[0031] Furthermore, the distribution consistency index cross_correlation is calculated by the following formula:

[0032]

[0033] Where a is the first coefficient, b is the second coefficient, kstest(C,D) represents the probability that the first sequence C and the second sequence D obey the same distribution obtained by performing KS hypothesis test on the first sequence C and the second sequence D; d i is the number of initial single-component signals, L is the number of symplectic geometric components, SGC i is the i-th symplectic geometric component, SGC j is the jth symplectic geometric component, inti_SGCd i is the embedding dimension d i The initial single signal component, sig is the original analysis signal, d i is the embedding dimension.

[0034] Furthermore, the improved sideband power value SPI of the symplectic geometric component is calculated by the following formula:

[0035]

[0036] Among them, P max is the local maximum value of the characteristic frequency sideband range, mean(P f ) is the average power in the sideband range, P f,j is the power amplitude in the sideband range, P T is the total energy in the frequency domain, and K is the number of power spectra in the sideband.

[0037] A bearing fault diagnosis system, comprising:

[0038] The signal acquisition module is used to obtain the vibration signal of the bearing throughout its life cycle from a healthy state to a damaged state, and obtain the vibration acceleration signal of different fault states;

[0039] A sequence construction module is used to perform Fourier transform on the vibration acceleration signal of the bearing throughout its life cycle to obtain a power spectrum, determine the initial value of the embedding dimension based on the power spectrum, and construct a candidate embedding dimension sequence;

[0040] A search module is used to start searching from the starting dimension of the candidate embedding dimension sequence, obtain the symplectic geometric component and the initial single signal component, and repeat the search until all the candidate embedding dimension sequences are searched to obtain the optimal embedding dimension corresponding to the minimum distribution consistency index;

[0041] A decomposition module is used to perform symplectic geometric mode decomposition based on Wasserstein distance according to the optimal embedding dimension corresponding to the minimum distribution consistency index to obtain the optimal symplectic geometric component;

[0042] An improved sideband power value calculation module is used to calculate the improved sideband power value of the symplectic geometric component according to the optimal symplectic geometric component;

[0043] The parsing module takes the symplectic geometric components whose improved sideband power values ​​exceed the upper threshold as fault sensitive features, performs envelope demodulation analysis on the fault sensitive features, and determines the fault type.

[0044] An electronic device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the bearing fault diagnosis method is implemented when the processor executes the computer program.

[0045] A computer-readable storage medium stores a computer program, and the computer program implements the bearing fault diagnosis method when executed by a processor.

[0046] Compared with the prior art, the present invention has the following beneficial technical effects:

[0047] The present invention realizes the diagnosis of bearing faults based on improved symplectic geometric modal decomposition and distribution consistency index. Compared with the VMD decomposition method of selecting parameters according to artificial experience, the present invention can adaptively select key parameters, and solves the problem of poor decomposition accuracy caused by selecting parameters according to artificial experience; and the symplectic geometric modal decomposition based on Wasserstein distance not only has powerful signal decomposition energy, but also can retain the essential characteristics of the vibration signal during the decomposition process, while realizing fault feature decoupling and weak feature extraction, it ensures that the decomposed information can fully reflect the operating status of the equipment, and overcomes the problem that the existing bearing feature extraction method is difficult to realize feature decoupling and weak feature extraction from low signal-to-noise ratio and feature coupling signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0049] Figure 1 A flow chart of a bearing fault diagnosis method based on improved symplectic geometric modal decomposition and distribution consistency index proposed by the present invention;

[0050] Figure 2 This is a vibration waveform diagram of an early composite fault of a bearing according to an embodiment of the present invention;

[0051] Figure 3 The envelope spectrum of the early composite fault of the bearing according to the embodiment of the present invention;

[0052] Figure 4This is a graph of the objective function results in the optimal parameter determination method based on distribution consistency proposed by the present invention;

[0053] Figure 5 This is a result diagram of the improved sideband power index proposed by the present invention;

[0054] Figure 6 Hilbert envelope spectra of fault-sensitive components screened out by the present invention; wherein (a) is the envelope spectrum of SGC6, and (b) is the envelope spectrum of SGC12;

[0055] Figure 7 Hilbert envelope spectrum of CEEMDAN decomposition results; (a) is the envelope spectrum of IMF3, (b) is the envelope spectrum of IMF5;

[0056] Figure 8 It is a schematic diagram of the bearing fault diagnosis system of the present invention. DETAILED DESCRIPTION

[0057] The present invention is further described in detail below in conjunction with specific embodiments, which are intended to explain the present invention rather than to limit it.

[0058] Symplectic geometry mode decomposition (SGMD) is a method that can decompose nonlinear signals, reflect nonlinear structures, and is suitable for analyzing nonlinear signals. Its core idea is to obtain eigenvalues ​​and eigenvectors in the symplectic geometry framework by performing symplectic geometry similarity transformation on the trajectory matrix (constructed by the analysis signal), and to reconstruct the single-component signal of the shock by selecting the corresponding eigenvectors through correlation indicators.

[0059] Reference Figure 1As shown, the bearing fault diagnosis method based on improved symplectic geometric modal decomposition and distribution consistency index of the present invention determines the initial optimal embedding dimension based on the power spectrum of the strong interference weak impact signal and constructs the embedding dimension sequence; secondly, the symplectic geometric modal decomposition based on the Wasserstein distance is performed under different embedding dimension states to obtain the decomposed symplectic geometric components (SGCs) and the initial single signal components (Inti_SGCs); then, the optimal embedding dimension of the symplectic geometric modal decomposition is determined based on the distribution consistency index; then, the trajectory matrix is ​​constructed based on the optimal embedding dimension, and the symplectic geometric modal decomposition is performed on it to obtain the symplectic geometric components; further, the fault-sensitive components are screened out from the numerous symplectic geometric components based on the improved sideband power index. Finally, the fault-sensitive components are demodulated by Hilbert envelope to realize the bearing fault mode discrimination. The present invention can overcome the problem that the existing bearing feature extraction method is difficult to realize feature decoupling and weak feature extraction from low signal-to-noise ratio and feature coupling signals, and realize fault feature decoupling and weak feature extraction. Specifically includes the following steps:

[0060] (1) Vibration signal acquisition: The acceleration sensor is used to measure the full life cycle vibration signal from the healthy state to the damaged state with strong interference and weak impact characteristics under high speed and heavy load conditions, and the vibration acceleration signals of different fault states are obtained.

[0061] (2) Determine the initial embedding dimension and embedding dimension search sequence based on the power spectrum PSD: Perform Fourier transform on the collected vibration acceleration signals of different fault states to obtain their power spectrum PSD, and determine the embedding dimension d based on the power spectrum. p The initial value of the embedding dimension d is constructed to construct the search sequence of the candidate embedding dimension sequence; the specific process is as follows:

[0062] First, the bearing vibration signal is transformed by Fourier transform to obtain its power spectrum PSD, and the fundamental frequency and rotation frequency F of the signal are determined according to the power spectrum. r , according to f∈[F r -Δf,F r The embedding dimension determined by the frequency corresponding to the maximum power spectral density within [Δf] is taken as the initial value d p , frequency f and initial embedding dimension d p Relationship:

[0063]

[0064] Among them, F s is the sampling frequency, f max For f∈[F r -Δf,F r +Δf] has the frequency corresponding to the maximum power spectral density. Δf is used to constrain the frequency band range. (0.5~0.8)F is recommended. r, n is the signal length;

[0065] Secondly, according to the initial embedding dimension d p and parameter delta to construct a candidate embedding dimension sequence D w :

[0066] D w =[d p -delta,…,d p +delta]

[0067] delta=max(15,c×d p )

[0068] Among them, c is the weight coefficient, which is used to constrain the length of the embedding dimension sequence.

[0069] (3) Based on the symplectic geometric mode decomposition method of Wasserstein distance, d i ∈D w As a variable parameter, construct an embedding dimension of d i The trajectory matrix of the case is decomposed by the symplectic geometric mode decomposition method based on Wasserstein distance to obtain the corresponding symplectic geometric components SGCs and the initial single signal components inti_SGCs; the specific process is as follows:

[0070] First, according to Takens embedding theorem, the original analysis signal sig = x1, x2, ..., x n (n is the signal length) Using topological methods, it can be expanded into a multidimensional time series matrix (usually called a trajectory matrix):

[0071]

[0072] Where d is the embedding dimension, τ represents the delay time, and there are m=n-(d-1)τ rows.

[0073] Secondly, the covariance matrix A can be obtained by performing autocorrelation analysis on the trajectory matrix X. Based on the covariance matrix A, the Hamilton matrix W can be further established. The Hamilton matrix W can be used to construct the symplectic geometry matrix G that meets the requirements of symplectic geometry analysis:

[0074] A=X T X

[0075]

[0076] W=M 2

[0077]

[0078] Among them, M is a diagonal matrix.

[0079] Then, according to W = M 2 Convert it to Without loss of generality, assume that the eigenvalues of the covariance matrix A are And λ1, λ2, …, λ d is the first eigenvalue, the second eigenvalue, …, the d-th eigenvalue of matrix B; Since the symplectic geometric similarity transformation can ensure that the time series still maintains its essential characteristics after transformation, according to the square relationship between the eigenvalues of covariance matrices A and B, the problem of solving matrix B can be simplified to the problem of solving the eigenvalues and eigenmatrices of matrix A; Use QR decomposition to decompose matrix A 2 to obtain the eigenvector matrix Q. Construct the coefficient matrix S and the reconstruction matrix Z based on the eigenvector matrix Q and the original trajectory matrix X:

[0080] S i = Q i T X

[0081] Z i = Q i T S i

[0082] Z = Z1 + Z2 +... + Z d

[0083] Among them, Q i (i = 1, 2, …, d) is the eigenvector corresponding to each eigenvalue of matrix A; X is the trajectory matrix; Z i (i = 1, 2, …, d) is called the initial single-component reconstruction matrix; S i is the element of the coefficient matrix S.

[0084] Subsequently, for the convenience of subsequent analysis and decomposition effect description, use diagonal averaging to convert the m×d-dimensional initial single-component matrix into a one-dimensional vector to obtain the initial single-signal component (inti_SGCs). The elements in Z i are z ij (1 ≤ i ≤ m, 1 ≤ j ≤ d), the minimum value of the number of rows and the embedding dimension d* = min(m, d), the maximum value of the number of rows and the embedding dimension m* = max(m, d), and 1 ≤ i ≤ m, 1 ≤ j ≤ d, where i is the number of rows and j is the number of columns. If m < d, then there is a coordinate transformation value z ij * = z ij , otherwise z ij * = z ji ;

[0085] Each element y of the initial single signal component matrix k for:

[0086]

[0087] Finally, assuming that d initial single signal components are decomposed, the first initial single signal component is used as the benchmark to evaluate the Wasserstein distance between it and the remaining d-1 vectors, and the initial single signal components exceeding the threshold Wasserstein distance are reorganized to obtain the first symplectic geometric component (SGC1); then, the initial single signal component with the largest energy proportion is selected from the remaining initial decomposition results as the benchmark, and the above process is repeated, and the symplectic geometric component (SGC s );

[0088] The Wasserstein distance W |p,q| The calculation formula is as follows:

[0089]

[0090] Among them, (p,q) --The set of all possible joint distributions of distribution p and distribution q. Sample x and sample y are obtained from the joint distribution (x, y) ~ γ, and the expected value E of the sample distance is calculated under the joint distribution γ. (x,y)～γ ||xy||

[0091] (4) Distribution consistency index calculation: Based on the symplectic geometric components SGCs obtained by symplectic geometric mode decomposition and the initial single signal components inti_SGCs, the distribution consistency index is used as the objective function to calculate the embedding dimension d i The objective function value when ; the specific process is as follows:

[0092] Firstly, based on the Kolmogorov-Smirnov (KS) hypothesis test, the probability that the initial single signal component inti_SGCs and the original analysis signal (sig) obey the same distribution and the probability that each symplectic geometric component (SGCs) obeys the same distribution are evaluated respectively.

[0093] Then, according to the distribution consistency calculation formula, d = d i In this case, the value of the distribution consistency index is:

[0094]

[0095] The coefficients a and b are 0.5 and 1 respectively; kstest(A,B) represents the probability that A and B follow the same distribution after the KS hypothesis test on the series A and B, d iis the number of initial single-component signals, L is the number of symplectic geometric components, SGC i is the i-th symplectic geometric component, SGC j is the jth symplectic geometric component.

[0096] (5) Determine the optimal parameters based on distribution consistency: Based on the candidate embedding dimension sequence, starting from the starting dimension of the sequence, i.e. the first point, repeat steps 3 and 4 until the entire sequence is traversed and the optimal embedding dimension corresponding to the minimum objective function is obtained. The optimal embedding dimension Substitute it into the symplectic geometric mode decomposition based on Wasserstein distance to obtain the optimal symplectic geometric components SGCs;

[0097] (6) Determine the fault sensitive component based on the improved sideband power index: According to the symplectic geometric mode decomposition result based on the optimal embedding dimension, calculate the improved sideband power value of the symplectic geometric component SGCs, set the threshold according to the 3σ criterion, and take the SGCs whose improved sideband power index exceeds the upper limit of the threshold as the fault sensitive feature; perform envelope demodulation analysis on the extracted fault sensitive features to determine the fault type; the specific process is as follows:

[0098] Firstly, the symplectic geometric components SGCs obtained by SGMD are subjected to Hilbert transform and Fourier transform to obtain the envelope spectrum of SGCs;

[0099] Secondly, the sideband power index is calculated according to the envelope spectrum of each SGCs, and the threshold is set according to the 3σ criterion. The SGCs whose improved sideband power index exceeds the upper limit of the threshold are taken as fault sensitive features. The rolling bearing fault characteristic frequency extracted from the fault sensitive feature envelope spectrum is compared with the theoretical bearing fault characteristic frequency to determine the bearing fault type.

[0100]

[0101] Among them, P max is the local maximum value of the characteristic frequency sideband range, mean(P f ) is the average power in the sideband range; P max / mean(P f ) represents the prominence of the envelope spectrum fault characteristics, while The energy proportion that characterizes the local characteristics of the envelope spectrum, P f,j is the power amplitude in the sideband range, P T is the total energy in the frequency domain, and K is the number of power spectra in the sideband.

[0102] See also Figure 8 One embodiment of the present invention provides a bearing fault diagnosis system, comprising:

[0103] The signal acquisition module is used to obtain the full life cycle vibration signal of the bearing from the healthy state to the damaged state with the characteristics of strong interference and weak impact, and obtain the vibration acceleration signal of different fault states;

[0104] A sequence construction module is used to perform Fourier transform on the vibration acceleration signal of the bearing throughout its life cycle to obtain a power spectrum, determine the initial value of the embedding dimension based on the power spectrum, and construct a candidate embedding dimension sequence;

[0105] A search module is used to start searching from the starting dimension of the candidate embedding dimension sequence, obtain the symplectic geometric component and the initial single signal component, and repeat the search until all the candidate embedding dimension sequences are searched to obtain the optimal embedding dimension corresponding to the minimum distribution consistency index;

[0106] A decomposition module is used to perform symplectic geometric mode decomposition based on Wasserstein distance according to the optimal embedding dimension corresponding to the minimum distribution consistency index to obtain the optimal symplectic geometric component;

[0107] An improved sideband power value calculation module is used to calculate the improved sideband power value of the symplectic geometric component according to the optimal symplectic geometric component;

[0108] The parsing module takes the symplectic geometric components whose improved sideband power values ​​exceed the upper threshold as fault sensitive features, performs envelope demodulation analysis on the fault sensitive features, and determines the fault type.

[0109] An embodiment of the present invention provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the bearing fault diagnosis method is implemented when the processor executes the computer program.

[0110] An embodiment of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the bearing fault diagnosis method is implemented.

[0111] A specific application example process is given below, and the effectiveness of the present invention in engineering applications is verified.

[0112] In this example, rolling bearings are taken as the research object, and data are collected using the XJTU-SY rolling bearing test bench; the test bench consists of an AC motor, a motor speed controller, a support shaft, two support bearings (heavy-loaded roller bearings) and a hydraulic loading system. The motor speed controller is connected to the AC motor, and the shaft of the AC motor passes through the inner ring of the test bearing and is connected to the test bearing. The support bearing is used to support the shaft. A support bearing is set at each end of the test bearing. The test bench is used to perform accelerated degradation tests on the test bearing under different working conditions to obtain complete run-to-failure data. The radial force is generated by the hydraulic loading system and applied to the housing of the test bearing, and the speed is set and maintained by the motor speed controller of the AC motor. In order to collect the vibration signal of the test bearing, two PCB 352C33 accelerometers are placed at the position of the housing of the test bearing, that is, one is placed on the horizontal axis and the other is placed on the vertical axis. The sampling frequency is set to 25.6kHz, and 32,768 data points (i.e., 1.28s) are recorded every 1 minute. At a speed of 2100rpm, the detailed information of the five sets of data tested is shown in Table 1.

[0113] Table 1XJTU-SY 2100rpm bearing data set information

[0114]

[0115] Taking the Bearing 1_5 dataset as the analysis object, the vibration data of the 25th second is intercepted as an example verification. The time domain diagram of the signal is as follows: Figure 2 As shown, the envelope spectrum is Figure 3 Shown

[0116] Using the method of the present invention, the bearing fault diagnosis method based on improved symplectic geometric mode decomposition and distribution consistency index is as follows:

[0117] (1) Vibration signal acquisition:

[0118] The acceleration sensor is used to measure the vibration signal of the entire life cycle from the healthy state to the damaged state under high speed and heavy load conditions with the characteristics of strong interference and weak impact, and the vibration acceleration signals of different fault states are obtained.

[0119] (2) Determine the initial embedding dimension and embedding dimension search sequence based on PSD

[0120] First, the Fourier transform is performed on the analyzed object to obtain its power spectrum PSD, and the fundamental frequency and the conversion frequency F of the signal are determined according to the power spectrum. r is 34 Hz, according to f∈F r -Δf,F r +Δf] with the maximum power spectral density determines the embedding dimension as the initial value d p is 84, f and the initial embedding dimension dp Relationship:

[0121]

[0122] Among them, the sampling frequency F s is 25600Hz, f max For f∈[F r -Δf,F r +Δf] has the frequency corresponding to the maximum power spectral density, the specific value is 403, Δf is used to constrain the frequency band range, in this case it is set to 0.5F r , n is the signal length;

[0123] Secondly, according to the initial embedding dimension d p And parameter delta constructs the candidate embedding dimension sequence D w :

[0124] D w =[d p -delta,…,d p +delta]

[0125] delta=max(15,0.3*d p )

[0126] (3) Symplectic geometric mode decomposition method based on Wasserstein distance; i ∈D w As a variable parameter, construct an embedding dimension of d i The trajectory matrix of the case is decomposed by the symplectic geometric mode decomposition method based on Wasserstein distance to obtain the corresponding symplectic geometric components SGCs and the initial single signal components inti_SGCs; the specific process is as follows:

[0127] First, according to Takens embedding theorem, we transform the one-dimensional time series x=x1,x2,…,x n (n is the signal length) Using topological methods, it can be expanded into a multidimensional time series matrix (usually called a trajectory matrix):

[0128]

[0129] Where d is the embedding dimension, which changes with iterations; and τ represents the delay time, which is set to 1, and m=n-(d-1)τ.

[0130] Secondly, the covariance matrix A can be obtained by performing autocorrelation analysis on the trajectory matrix X. Based on the covariance matrix A, the Hamilton matrices M and W can be further established. The Hamilton matrices M and W can be used to construct the symplectic geometry matrix G that meets the requirements of symplectic geometry analysis:

[0131] A=X T X

[0132]

[0133] W=M 2

[0134]

[0135] Then, according to W = M 2 Convert it to Let's assume that the moment

[0136] The eigenvalues ​​of the matrix A are And λ1,λ2,…,λ d is the eigenvalue of matrix B; since symplectic geometric similarity transformation can ensure that the time series still maintains its essential characteristics after transformation, according to the square relationship between the eigenvalues ​​of matrices A and B, the problem of solving matrix B can be simplified to the problem of solving the eigenvalue and characteristic matrix of matrix A; using QR decomposition to solve matrix A 2 Decompose and obtain the eigenvector matrix Q. Construct the coefficient matrix S and reconstruction matrix Z based on the eigenvector matrix Q and the original trajectory matrix X:

[0137] S i =Q i T X

[0138] Z i =Q i T S i

[0139] Z=Z1+Z2+...+Z d

[0140] Among them, Q i (i=1,2,…,d) is the eigenvector corresponding to each eigenvalue of matrix A; X is the trajectory matrix; Z i (i=1,2,…,d) is called the initial single-component reconstruction matrix;

[0141] Subsequently, in order to facilitate subsequent analysis and decomposition effect description, the initial single component matrix of m×d dimensions is converted into a one-dimensional vector by diagonal averaging to obtain the initial single signal component (inti_SGCs). i The element in is z ij(1 ≤ i ≤ m, 1 ≤ j ≤ d), d* = min(m, d), m* = max(m, d) and 1 ≤ i ≤ m, 1 ≤ j ≤ d. If m < d, then there is z ij * = z ij , otherwise z ij * = z ji :

[0142]

[0143] Finally, set the threshold to 0.005. Assuming that d initial single-signal components are obtained by decomposition, first, taking the first initial single-signal component as a reference, evaluate the Wasserstein distance between it and the remaining d - 1 vectors, and recombine the initial single-signal components with a Wasserstein distance exceeding the threshold to obtain the first symplectic geometric component (SGC1); subsequently, select the initial single-signal component with the largest energy proportion from the remaining initial decomposition results as a reference, repeat the above process, and continuously perform recombination to obtain the symplectic geometric components (SGC s );

[0144]

[0145] where, ∏ (p,q) -- the set of all possible joint distributions combining distribution p and distribution q. Samples x and y are sampled from the joint distribution (x, y) ~ γ, and the expected value E of the sample distance is calculated under the joint distribution γ (x,y)～γ ||x - y||

[0146] (4) Calculation of distribution consistency index: Based on the symplectic geometric components SGCs and the initial single-signal components inti_SGCs obtained by symplectic geometric mode decomposition, with the distribution consistency index as the objective function, calculate the value of the objective function when the embedding dimension is d i ; the specific process is as follows:

[0147] First, based on the Kolmogorov - Smirnov (K - S) hypothesis test, evaluate the probability that the initial single-signal components inti_SGCs and the original analysis signal (sig) follow the same distribution, and the probability that the symplectic geometric components (SGCs) follow the same distribution respectively.

[0148] Subsequently, according to the distribution consistency calculation formula, calculate the value of the distribution consistency index when d = d i :

[0149]

[0150] Among them, the coefficients a and b are 0.5 and 1 respectively; kstest(A,B) represents the probability that A and B obey the same distribution obtained by performing KS hypothesis test on series A and B.

[0151] (5) Determine the optimal parameters based on distribution consistency: Based on the candidate embedding dimension sequence, starting from the first point in the sequence, repeat steps 3 and 4 until the entire sequence is traversed; obtain the optimal embedding dimension corresponding to the minimum objective function Determine the optimal embedding dimension to be 81; Substitute it into the symplectic geometric mode decomposition based on Wasserstein distance to obtain the optimal decomposition result. The objective function values ​​corresponding to different embedding dimensions are as follows: Figure 4 As shown;

[0152] (6) According to the symplectic geometric mode decomposition results based on the optimal embedding dimension, the improved sideband power values ​​of the symplectic geometric components SGCs are calculated, and the threshold is set according to the 3σ criterion. The SGCs exceeding the upper threshold are taken as fault sensitive features; the extracted fault sensitive features are subjected to envelope demodulation analysis to determine the fault type. The specific process is as follows:

[0153] Firstly, the symplectic geometric components SGCs obtained by SGMD are subjected to Hilbert transform and Fourier transform to obtain the envelope spectrum of SGCs;

[0154] Secondly, the improved sideband power index is calculated based on the envelope spectrum of each SGCs, as follows: Figure 5 As shown; the upper limit of the improved sideband power index is determined according to the 3σ criterion, and Hilbert envelope demodulation is performed on SGC6 and SGC12 that exceed the upper limit. The bearing fault characteristics extracted from the envelope spectrum are the bearing fault characteristic frequency and its sideband. The rolling bearing fault characteristic frequency extracted from the envelope spectrum is compared with the theoretical bearing fault characteristic frequency to determine the bearing fault type.

[0155]

[0156] Among them, P max is the local maximum value of the characteristic frequency sideband range, mean(P f ) is the average power in the sideband range; P max / mean(P f ) represents the prominence of the envelope spectrum fault characteristics, while The energy proportion that characterizes the local characteristics of the envelope spectrum; P f,j is the power amplitude in the sideband range, P T is the total energy in the frequency domain, and K is the number of power spectra in the sideband.

[0157] In order to verify the effectiveness of the present invention, the complete ensemble empirical mode decomposition of adaptive noise (CEEMDAN) which is widely used in signal decomposition and feature extraction is selected for comparison. The envelope spectra of the screening components SGC6 and SGC12 of the proposed method are shown in Figure 2. Figure 6 As shown in (a) and (b), Figure 6 In the envelope spectrum of SGC6 shown in (a), the fundamental frequency, fundamental frequency multiples, outer race fault characteristic frequency and its sidebands can be clearly identified. Figure 6 In the envelope spectrum of SGC12 shown in (b), the fundamental frequency, fundamental frequency multiples, inner ring fault frequency and its sidebands can be clearly identified; the envelope spectrum of the components extracted by CEEMDAN is shown in Figure 7 As shown in (a) and (b), the analysis Figure 7 In (a), it is found that the envelope spectrum of IMF3 contains not only the characteristic frequency of the inner circle fault and its sideband, but also the characteristic frequency of the outer circle fault. Figure 7 Although IMF5 reflected in (b) can identify the characteristic frequency of the outer ring fault, the amplitude is low, which shows that CEEMDAN has not achieved complete separation of the inner and outer ring composite faults. Figure 6 and Figure 7 It can be seen that the method proposed in the present invention can realize the decoupling and separation of composite faults and can extract weak fault features from strong interference and weak impact signals.

[0158] Based on the composite fault data in the XJTU-SY data set, the present invention realizes feature decoupling by using the symplectic geometric modal decomposition method based on the Wasserstein distance, and determines the optimal parameters based on the distribution consistency, thereby getting rid of the influence of manual selection on the decomposition accuracy, reducing noise interference to the greatest extent, and realizing the decoupling function of the coupled signal; in view of the problem that there are too many symplectic geometric modal decomposition components and manual screening is time-consuming and labor-intensive, a fault-sensitive frequency band screening method based on the improved sideband power index is proposed, and the adaptive screening of fault-sensitive features is realized; compared with other traditional methods, the present invention can overcome the problem that the existing feature extraction methods are difficult to realize feature decoupling and weak feature extraction from low signal-to-noise ratio and feature coupling signals, retain the essential characteristics of the vibration signal in the decomposition process, and ensure that the decomposed information can fully reflect the operating status of the equipment while realizing fault feature decoupling and weak feature extraction.

[0159] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A bearing fault diagnosis method, characterized in that: The steps include: Obtain vibration acceleration signals with strong interference and weak impact characteristics of different bearing fault states; The vibration acceleration signal of the bearing during its entire life cycle is transformed into a power spectrum by Fourier transform, and the initial value of the embedding dimension is determined based on the power spectrum, and a candidate embedding dimension sequence is constructed. Starting the search from the starting dimension of the candidate embedding dimension sequence, obtaining the symplectic geometric component and the initial single signal component, and repeating the search until the entire candidate embedding dimension sequence is searched, obtaining the optimal embedding dimension corresponding to the minimum distribution consistency index; According to the optimal embedding dimension corresponding to the minimum distribution consistency index, the symplectic geometric mode decomposition based on Wasserstein distance is performed to obtain the optimal symplectic geometric component; According to the optimal symplectic geometric component, an improved sideband power value of the symplectic geometric component is calculated; The symplectic geometric components whose improved sideband power values ​​exceed the upper threshold are taken as fault sensitive features, and the fault sensitive features are subjected to envelope demodulation analysis to determine the fault type; The distribution consistency indicator cross_correlation is calculated by the following formula: Where a is the first coefficient, b is the second coefficient, kstest(C,D) represents the probability that the first sequence C and the second sequence D obey the same distribution obtained by performing KS hypothesis test on the first sequence C and the second sequence D; d i is the number of initial single-component signals, L is the number of symplectic geometric components, SGC i is the i-th symplectic geometric component, SGC j is the jth symplectic geometric component, inti_SGCd i is the embedding dimension d i The initial single signal component, sig is the original analysis signal, d i is the embedding dimension; The improved sideband power value SPI of the symplectic geometric component is calculated by the following formula: Among them, P max is the local maximum value of the characteristic frequency sideband range, mean(P f ) is the average power in the sideband range, P f,j is the power amplitude in the sideband range, P T is the total energy in the frequency domain, and K is the number of power spectra in the sideband.

2. The bearing fault diagnosis method according to claim 1, characterized in that: The power spectrum is obtained by Fourier transforming the vibration acceleration signal of the bearing throughout its life cycle, and the initial value of the embedding dimension is determined according to the power spectrum, and a candidate embedding dimension sequence is constructed, including the following steps: The bearing vibration signal is transformed by Fourier transform to obtain the power spectrum, and the fundamental frequency and rotation frequency F of the signal are determined according to the power spectrum. r , according to f∈[F r -Δf,F r +Δf] with the maximum power spectral density is used as the initial embedding dimension, and a candidate embedding dimension sequence is constructed according to the initial embedding dimension, where f is the frequency.

3. The bearing fault diagnosis method according to claim 2, characterized in that: Frequency f and initial embedding dimension d p The relationship is: Among them, F s is the sampling frequency, f max For f∈[F r -Δf,F r +Δf] has the frequency corresponding to the maximum power spectral density, Δf is used to constrain the frequency band range, and n is the signal length; Candidate embedding dimension sequence D w for: D w =[d p -delta,…,d p +delta] delta=max(15,c×d p ) Among them, delta is the parameter and c is the weight coefficient.

4. The bearing fault diagnosis method according to claim 1, characterized in that: The symplectic geometric components and the initial single signal components are obtained by the following process: The initial value of the embedding dimension is taken as a variable parameter to construct the trajectory matrix with the set embedding dimension. The trajectory matrix is ​​decomposed using the symplectic geometric mode decomposition method based on Wasserstein distance to obtain the symplectic geometric components and the initial single signal components.

5. The bearing fault diagnosis method according to claim 1 or 4, characterized in that: The symplectic geometry components are obtained by the following process: According to Takens embedding theorem, the original analysis signal sig=x1,x2,…,x n , n is the signal length, which is expanded into a trajectory matrix using topological methods; Perform autocorrelation analysis on the trajectory matrix X to obtain the covariance matrix; And based on the covariance matrix, a Hamilton matrix is ​​established; and the Hamilton matrix is ​​used to construct a symplectic geometry matrix that meets the requirements of symplectic geometry analysis: The symplectic geometry matrix that meets the requirements of symplectic geometry analysis is solved by QR to obtain the eigenvector matrix, and the initial single-component matrix is ​​constructed according to the eigenvector matrix and the trajectory matrix; The multi-dimensional initial single component matrix is ​​converted into a one-dimensional vector by using diagonal averaging to obtain the initial single signal component; Assuming that d initial single signal components are obtained by decomposition, firstly, the first initial single signal component is used as a benchmark to evaluate the Wasserstein distance between it and the remaining d-1 initial single signal components, and the initial single signal components whose Wasserstein distance exceeds the threshold are reorganized to obtain the first symplectic geometric component; Then, the initial single signal component with the largest energy proportion is selected from the remaining initial single signal components as a reference, and the components are continuously reorganized to obtain the symplectic geometric components.

6. A bearing fault diagnosis system, characterized in that: include: The signal acquisition module is used to obtain the vibration signal of the bearing throughout its life cycle from a healthy state to a damaged state, and obtain the vibration acceleration signal of different fault states; A sequence construction module is used to perform Fourier transform on the vibration acceleration signal of the bearing throughout its life cycle to obtain a power spectrum, determine the initial value of the embedding dimension based on the power spectrum, and construct a candidate embedding dimension sequence; A search module is used to start searching from the starting dimension of the candidate embedding dimension sequence, obtain the symplectic geometric component and the initial single signal component, and repeat the search until all the candidate embedding dimension sequences are searched to obtain the optimal embedding dimension corresponding to the minimum distribution consistency index; A decomposition module is used to perform symplectic geometric mode decomposition based on Wasserstein distance according to the optimal embedding dimension corresponding to the minimum distribution consistency index to obtain the optimal symplectic geometric component; An improved sideband power value calculation module is used to calculate the improved sideband power value of the symplectic geometric component according to the optimal symplectic geometric component; The analysis module takes the symplectic geometric component whose improved sideband power value exceeds the upper threshold as the fault sensitive feature, performs envelope demodulation analysis on the fault sensitive feature, and determines the fault type; The distribution consistency indicator cross_correlation is calculated by the following formula: Where a is the first coefficient, b is the second coefficient, kstest(C,D) represents the probability that the first sequence C and the second sequence D obey the same distribution obtained by performing KS hypothesis test on the first sequence C and the second sequence D; d i is the number of initial single-component signals, L is the number of symplectic geometric components, SGC i is the i-th symplectic geometric component, SGC j is the jth symplectic geometric component, inti_SGCd i is the embedding dimension d i The initial single signal component, sig is the original analysis signal, d i is the embedding dimension; The improved sideband power value SPI of the symplectic geometric component is calculated by the following formula: Among them, P max is the local maximum value of the characteristic frequency sideband range, mean(P f ) is the average power in the sideband range, P f,j is the power amplitude in the sideband range, P T is the total energy in the frequency domain, and K is the number of power spectra in the sideband.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the bearing fault diagnosis method as described in any one of claims 1 to 5 is implemented.

8. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the bearing fault diagnosis method according to any one of claims 1 to 5 is implemented.

Images