Bearing fault diagnosis method

Abstract

The invention discloses a bearing fault diagnosis method, which comprises the following steps of: firstly, acquiring vibration signals of a bearing at different positions and directions under health and various fault states by adopting multiple sensors, denoising original signals of the bearing by utilizing an SG filtering method, and carrying out feature fusion to improve the adaptability and reliability of the signals; secondly, decomposing the denoised bearing vibration signal by using a regenerative phase shift sine-assisted empirical mode decomposition method to obtain a series of vibration signal component data IMF, and selecting an optimal TE-IMF by using transfer entropy to obtain recombined signal data; secondly, improving an initialization stage and a population search stage of a non-monopoly search algorithm NO by using sparse initialization and a probability disturbance strategy to obtain a CNO; and finally, optimizing hyper-parameters in the bearing fault diagnosis model based on SageFormer by using CNO, and carrying out bearing fault diagnosis by using the trained model. According to the invention, the accuracy of the bearing fault diagnosis result can be effectively improved.

Description

Technical Field

[0001] The invention belongs to the technical field of bearing fault diagnosis, and in particular relates to a bearing fault diagnosis method. Background Art

[0002] Bearing fault diagnosis is an important part of industrial equipment maintenance. It involves real-time monitoring and analysis of the bearing operating status in order to promptly detect and handle potential faults. As an important transmission component in mechanical equipment, the health of the bearing directly affects the performance and operating life of the entire equipment. Therefore, the research and application of bearing fault diagnosis technology has important engineering significance.

[0003] With the improvement of industrial automation and intelligence, higher requirements are placed on the reliability and maintenance cost-effectiveness of bearings. Fault diagnosis technology can help timely detect and warn problems in the operation of bearings, reduce downtime, and avoid major accidents that cause economic losses. Operating bearings involve a variety of equipment and complex physical processes, generating a large amount of data. How to extract useful information from this data and establish an accurate fault diagnosis model is a technical challenge. Machine learning technology, especially deep learning, data decomposition and feature extraction technology, makes it possible to process complex data sets and establish efficient fault diagnosis. In order to improve the accuracy of fault diagnosis, it is necessary to improve and optimize the existing machine learning algorithms to better adapt to specific application scenarios.

[0004] At present, the solutions to this problem can be mainly divided into three categories: sound signal processing methods, physical models and data-driven methods. For example, signal processing technology uses wavelet packet decomposition, Hilbert-Huang transform (HHT), short-time Fourier transform (STFT), etc. to analyze vibration signals and extract features that are helpful for fault diagnosis. However, traditional signal processing methods may ignore some important fault features, especially when processing complex signals, and it is difficult to extract accurate fault information. Support vector machine (SVM), random forest (RF), neural network, etc. are used to realize pattern recognition and classification of bearing faults. However, this method requires a large amount of labeled data for training, and dimensional disaster problems may occur when processing high-dimensional data. Deep learning, reinforcement learning and other algorithms are used in bearing fault diagnosis to automatically learn and extract features from a large amount of data to improve the accuracy and efficiency of diagnosis. However, the algorithm is highly complex, consumes a lot of computing resources, and some algorithms require a large amount of data for training. Using data-driven methods, this type of method usually does not require an in-depth understanding of the physical model of the bearing, but relies on a large amount of historical data, and for new fault types, it may be necessary to collect a large amount of data for training. Therefore, a bearing fault diagnosis method that can overcome the above technical problems is urgently needed. Summary of the invention

[0005] Purpose of the invention: The present invention proposes a bearing fault diagnosis method, which overcomes the problem that traditional signal processing methods ignore some important fault characteristics, and uses CNO to optimize the hyperparameters in the SageFormer model to improve the accuracy of fault detection through data-driven analysis.

[0006] Technical solution: A bearing fault diagnosis method according to the present invention comprises the following steps:

[0007] (1) SG filtering is used to denoise the vibration signals of the bearing at different positions under healthy and various fault conditions acquired in advance; and feature fusion is performed on the denoised multidimensional data;

[0008] (2) The denoised bearing vibration signal is decomposed using the regenerative phase-shifted sine-assisted empirical mode decomposition algorithm RPSEMD to obtain a series of vibration signal component data IMFs; the optimal T-IMF is selected using the transfer entropy TE to obtain the reorganized signal data; the reorganized signal data is normalized to construct the corresponding training set and test set;

[0009] (3) Construct a bearing fault diagnosis model based on SageFormer;

[0010] (4) The sparse initialization and probability perturbation strategies are used to improve the initialization phase and population search phase of the non-monopolistic search algorithm NO, and the improved non-monopolistic search algorithm CNO is obtained;

[0011] (5) Using the improved non-monopolistic search algorithm CNO to optimize the hyperparameters of the bearing fault diagnosis model based on SageFormer;

[0012] (6) The training set is used to train the bearing fault diagnosis model, the optimal parameters of the model are obtained, and the loss function error is minimized; the bearing fault diagnosis model after training and optimization is used to perform online fault diagnosis on the bearing vibration data to obtain the corresponding diagnosis results.

[0013] Furthermore, the implementation process of using SG filtering to perform denoising on the vibration signals of the pre-acquired healthy and various fault states of different positions of the bearing in step (1) is as follows:

[0014] First, a polynomial fit is performed on a set of 2H+1 sample center discrete points centered at m=0:

[0015]

[0016] Among them, T is the polynomial order, k represents the order of the discrete sequence value at this point, and z k is the polynomial coefficient, m is the discrete sequence value of the point, and y(m) is the fitting value of the point;

[0017] Secondly, calculate the sum E of the least - square fitting residuals between all the fitting data and the original data within the filtering window N :

[0018]

[0019] where x p (m) is the p - th group of samples of the original bearing vibration signal, and H represents the boundary of a group of bearing vibration signal samples;

[0020] When E N is the smallest, the fitting effect is the best; take the partial derivative of Equation (2) to obtain the polynomial fitting coefficient z k , calculate the fitting result and assign it to the mid - point of this window as the fitting result output, that is:

[0021] f[0] = y[0] (3)

[0022] Translate the filter one step along the time axis and repeat the above steps until the end of the data to achieve Savitzky - Golay (S - G) moving smoothing; use S - G convolution smoothing to replace the moving smoothing algorithm to expand the noise reduction dimension of the traditional S - G filter and improve the calculation efficiency of the fitting result. The specific implementation process is as shown in Equation (4):

[0023]

[0024] f i [m] represents the vibration signal collected by the i - th (i < w) sensor after noise reduction; after filtering and noise reduction for the w groups of local vibration signals collected by w sensors, we get:

[0025]

[0026] Furthermore, the process of feature fusion for the denoised multi - dimensional data in step (1) is as follows:

[0027] The correlation between any two signals f α [m], f β [m]α, α ∈ [1, w] is defined as: Define as:

[0028]

[0029] where t = 0, 1, 2, 3, …, N - 1 is the time series of the signal; the total correlation energy between the signal f α [m] and other signals f β [m] is:

[0030]

[0031] The energy of n signals is normalized separately, that is:

[0032]

[0033] After energy normalization, the signal The N discrete data points are Variance contribution rate K il Defined as:

[0034]

[0035] In the formula, For signal The first data point in The distribution coefficient k at the first data point il for:

[0036]

[0037] Based on the obtained distribution coefficient, w signals are fused into one signal, and the fused signal R(l) is:

[0038]

[0039] Furthermore, the implementation process of decomposing the denoised bearing vibration signal using the regenerative phase-shifted sine-assisted empirical mode decomposition algorithm RPSEMD in step (2) is as follows:

[0040] S1: Use EMD to process the signal x(t) and determine a based on the obtained intrinsic mode components k and f k θ ki Sampling between [0, 2π] with a phase shift of i (1≤i≤I), thus generating s(t|a k ,f k ,θ ki );

[0041] S2: Get the first intrinsic mode component and average all I to get the last intrinsic mode component c k (t);

[0042] S3: Remove c from x(t) k (t), x(t)←x(t)-c k (t), let k = k + 1;

[0043] S4: Repeat S1 to S3 until all the intrinsic modal components are obtained, and the final x(t) is used as the residual r(t); the fault sequence component IMF obtained by the RPSEMD algorithm is expressed as follows:

[0044] IMF×[imf1 imf2 … imf J ] (13)

[0045] Wherein, J represents the number of components decomposed by the RPSEMD algorithm.

[0046] Furthermore, the specific process of using the transfer entropy TE to select the optimal T-IMF in step (2) to obtain the reorganized signal data is as follows:

[0047] The probability histogram is used to replace the kernel probability density estimation method, and the probability density is defined as:

[0048]

[0049] In the formula, E represents r (r = 1, 2, 3) dimensional variables, and the interval of the variable e of each dimension [e min ,e max ] are divided into D equal parts, then the entire space is divided into D r The volume of each space is expressed as value δ represents the interval length of the δth variable; M represents the total number of samples; b α Indicates the number of samples in each small space;

[0050] Transfer entropy is a direct measure of the dynamic information transferred from one variable to another. Suppose that given two independent variables and They are s-order and k-order Markov processes respectively, and the transfer entropy is passed through the variable and variables Estimating variable x under historical data n+1 To evaluate the change in prediction effect; at time n+1, the entropy rate can be expressed as:

[0051]

[0052] In the formula, x n+1 Represents the sample of variable x at time n+1; Represents x n+1 and Probability of simultaneous occurrence; Indicates a given Under the condition x n+1 Probability of occurrence;

[0053] The joint conditional entropy rate is expressed as:

[0054]

[0055] In the formula, represents the joint probability density function; Indicated in and Under the condition x n+1 The probability density of ; According to equations (15) and (16), the information flow from variable y to variable x can be expressed as:

[0056]

[0057] In the formula, yes Under the condition x n+1 The probability density of ; Since the transfer entropy is directional, information transfer from variable x to variable y is also required; in summary, the transfer entropy from y to x is shown as:

[0058] TE y→x =T y→x -T x→y (18)

[0059] The fault sequence component T-IMF obtained after data reorganization after processing by transfer entropy is expressed as follows:

[0060] T-IMF = [imf1imf2…imf j ](j<J) (19)

[0061] Where j represents the number of fault sequence components obtained after reorganization by transfer entropy.

[0062] Furthermore, the implementation process of step (3) is as follows:

[0063] Before inputting the bearing vibration sequence, learnable global tokens are added to each sequence, which are used to capture the global information of each sequence; the input multivariate time series data is first converted into overlapping blocks, and each block is mapped to a unified latent space through a trainable linear projection. The global token embedding process of the bearing vibration sequence is as follows:

[0064] G(0) c ,:,:=[g1;...;g M ;x 1c E proj ;...;x Nc E proj ]+E pos (20)

[0065] Where G(0) is the embedding of the input bearing target sequence MTS sample, g i is the global token of the bearing sample vibration sequence, x ic is the i-th time step in the sampled vibration sequence c, E proj is a trainable linear projection of the bearing sample vibration sequence, E posis the vibration sequence position embedding;

[0066] The bearing vibration sequence adjacency matrix is ​​learned in an end-to-end manner to capture the relationship between each sampled bearing vibration sequence. Randomly initialized node embeddings are used and transformed through training parameters to learn unidirectional dependencies between sequences. Through a specific operation, the nearest neighbors of each node are determined and the weights of unconnected nodes are set to zero to generate the final sparse adjacency matrix:

[0067]

[0068] Where E is the node embedding of the bearing vibration sequence; Θ1 and Θ2 are trainable parameters; act is the nonlinear activation function of the SageFormer encoding layer; M1 and M2 represent the node representations obtained by transforming the activation function and parameter matrix used to operate the vibration sequence, which is used to capture the relationship between the sequences; ReLU(·) represents the intrinsic activation function of the SageFormer model, which is used to calculate the relationship between M1 and M2; A′ is the learned adjacency matrix between bearing vibration sequences, which is used to represent the preliminary relationship between sequences;

[0069] The embedded tokens are processed through a SageFormer encoding layer, which collects the global tokens of all sequences and processes them through a graph neural network component to fuse the information of each bearing vibration sequence and its neighbors:

[0070]

[0071] Among them, G i is the global token of the i-th layer bearing vibration sequence; is the graph Laplace matrix of the bearing vibration sequence; W d is the graph convolution filter of the bearing vibration sequence, which is used to filter and aggregate the neighbor features of the node during the graph aggregation process; D represents the number of hops that the information propagates in the graph during the graph aggregation stage of the SageFormer model;

[0072] After passing through multiple encoding layers, the output of the model is passed through a decoding head to produce the final prediction; the decoding head uses the information received from the encoder to predict future time steps:

[0073]

[0074] In the formula, It is the characteristic representation of the bearing vibration sequence after being processed by the encoder, and Y is the output of the bearing vibration sequence diagnosis result.

[0075] Furthermore, the implementation process of step (4) is as follows:

[0076] Exploration phase: The sparse initialization method is introduced into the exploration phase of the monopoly search algorithm NO. At the beginning of the algorithm, the initial position of the population is determined using random numbers. The algorithm explores the solution space based on the initial parameters obtained by the sparse initialization method to find new potential optimal solutions; the exploration phase is executed when the iteration counter (t) is less than or equal to half of the total number of iterations (T). The formula is as follows:

[0077] I i =rand(1,D)⊙(I max -I min )+I min (25)

[0078] Anew(j)=rand*A(RP)+I i (26)

[0079] Among them, I i is the position of the two individuals of the i-th SageFormer model learning rate lr and batch size hu, I max and I min are the upper and lower bounds of the learning rate lr and batch size hu in the search space, rand(1, D) is a vector of dimension D of uniformly distributed random numbers in the interval 0, 1, Anew(j) is the jth dimension in the newly generated solution of learning rate lr and batch size hu, rand is a random value between 0 and 1; j is the dimension index defined between 1 and the maximum dimension of the search space, A is the current solution, A(RP) is the selected random dimension of the current solution of learning rate lr and batch size hu; RP is a random number between 1 and D;

[0080] Development stage: Introduce the probability perturbation strategy into the population update process of the NO algorithm to help the algorithm escape from the local optimum; the specific implementation process is as follows:

[0081] Anew(j)=A(j)-(A(SRP)*rand)*eps-(A(j)-NO)(27)

[0082]

[0083] Where: P r is the disturbance probability; D a is the dimension of the algorithm optimization problem and D a =2; T a is the maximum number of iterations; in the early stage of algorithm iteration, the perturbation probability is small, so that the algorithm can quickly find the global optimal solution; in the later stage of algorithm iteration, the perturbation probability increases, which enhances the diversity of the population; T a is a large constant. When the algorithm iterates to the later stage t a →T a , disturbance probability:

[0084]

[0085] After adding the probability perturbation, the formula for updating the position of the algorithm particles is optimized as follows:

[0086] Anew(j)=A(j)-(A(SRP)*rand1)*eps-(A(j)-NO); rand1<P r (30)

[0087] where A(j) is the j-th position of the learning rate lr and batch size hu in the current solution, rand1 is a random value between [0, 1], A(SRP) is the dimension of the learning rate lr and batch size hu randomly selected from the current solution dimension, eps is a very small value, SRP is a solution random dimension, a random number between 1 and the dimension size; NO is a tuning parameter used to adjust the search process of the proposed method.

[0088] Furthermore, the implementation process of step (5) is as follows:

[0089] The matching function of the diagnosis results is used as the objective function of the SageFormer fault diagnosis model;

[0090] R = sum(S = S P ) / length(S)(31)

[0091] Among them, R represents the accuracy of bearing fault diagnosis, S represents the actual result of the bearing status, and S P Results indicating the diagnosed bearing condition;

[0092] Determine the learning rate lr and batch size hu of the population type SageFormer model, as follows:

[0093]

[0094] Among them, X represents the population of NO search algorithm, d represents the number of particles of one type in the population;

[0095] The upper and lower limits of the search space are determined by the range of the learning rate lr and the batch size hu; the specific expression is as follows:

[0096]

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

[0098] 1. Collecting vibration signals of bearings at different positions and directions through multiple sensors can more comprehensively capture the health status and fault characteristics of bearings, providing a richer data source for diagnosis;

[0099] 2. Using the Savitzky-Golay filtering method to denoise the original bearing signal and realize feature fusion can improve the signal quality and enhance the reliability of diagnosis; using the regenerative phase-shifted sine-assisted empirical mode decomposition method to decompose the denoised signal, and selecting the optimal decomposition component through transfer entropy, it is helpful to extract more accurate fault features;

[0100] 3. The monopoly search algorithm has problems such as slow convergence speed and easy to fall into local optimality during the optimization process. Therefore, the present invention adds sparse initialization and probability perturbation strategies to the population update process of the monopoly search algorithm to help the algorithm jump out of the local optimality, constructs the CNO algorithm, enhances the search ability of the original monopoly search algorithm, and improves the search graph efficiency;

[0101] 4. In view of the problems of low operating efficiency and large input data requirement in traditional bearing fault information diagnosis methods, the present invention adopts the SageFormer model to effectively simplify the complexity of the model and improve efficiency; and uses CNO to optimize the hyperparameters in the SageFormer model to achieve accurate diagnosis and identification of bearing faults. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] Figure 1 A schematic diagram of the SageFormer model structure constructed for the present invention;

[0103] Figure 2 A CNO algorithm flow chart is provided for the present invention. DETAILED DESCRIPTION

[0104] The present invention will be further described in detail below in conjunction with the accompanying drawings.

[0105] The present invention provides a bearing fault diagnosis method, comprising the following steps:

[0106] Step 1: Signal acquisition: Use multiple sensors to collect the vibration signals of the bearing at different positions under healthy and various fault conditions; Use the SG filtering method to denoise the vibration signals of the bearing under healthy and various types of fault conditions collected by multiple sensors, and perform feature fusion.

[0107] The sensors are reasonably installed at various important positions of the gears to sample the vibration information and obtain multi-dimensional vibration signal data. The multi-dimensional vibration signal data is filtered and denoised using Savitzky-Golay (SG). The specific implementation process is as follows:

[0108] First, a polynomial fit is performed on a set of 2H+1 sample center discrete points centered at m=0:

[0109]

[0110] Where: T is the polynomial order, k represents the order of the discrete sequence value at this point, z k is the polynomial coefficient, m is the discrete sequence value at this point, and y(m) is the fitting value at this point.

[0111] Secondly, find the sum of the least squares fitting residuals E of all the fitting data and the original data within the filtering window N :

[0112]

[0113] Where, x p (m) (p = 1, …, n) is the p-th group of samples of the original bearing vibration signal, and H represents the boundary of a group of bearing vibration signal samples.

[0114] When E N is the smallest, the fitting effect is the best. Take the partial derivative of Equation (2) to obtain the polynomial fitting coefficient z k , calculate the fitting result and assign it to the midpoint of this window as the fitting result output, that is:

[0115] f[0] = y[0] (3)

[0116] Translate the filter one step along the time axis and repeat the above steps until the end of the data to achieve S-G moving smoothing. Use S-G convolution smoothing to replace the moving smoothing algorithm to expand the noise reduction dimension of the traditional S-G filter and improve the calculation efficiency of the fitting result. The specific implementation process is as shown in Equation (4):

[0117]

[0118] f i [m] represents the vibration signal collected by the i-th (i < w) group of sensors after noise reduction.

[0119] After filtering and noise reduction of the w groups of local vibration signals collected by w sensors, we get:

[0120] F[m] = [f1[m] f2[m] … f w [m]] (5)

[0121] Perform feature fusion on the denoised multi-dimensional data. The specific implementation process is as follows:

[0122] The correlation between any two signals f α [m], f β [m] where α, β ∈ [1, w] can be defined as: It can be defined as:

[0123]

[0124] Where: t = 0, 1, 2, 3, ..., N-1 is the time series of the signal. Signal f α [m] and other signals f β The total correlation energy of [m] (α = 1, 2, 3, ... n, α ≠ β) is:

[0125]

[0126] Taking into account the dimensions and magnitudes of the signals, the energy of the n signals is normalized, namely:

[0127]

[0128] After energy normalization, the signal The N discrete data points are Variance contribution rate K il Defined as:

[0129]

[0130] In the formula, For signal The first data point in the signal The distribution coefficient k at the first data point il for:

[0131]

[0132] Based on the obtained distribution coefficient, w signals are fused into one signal, and the fused signal R(l) is:

[0133]

[0134] Step 2: Decompose the denoised bearing vibration signal using the regenerative phase-shifted sine-assisted empirical mode decomposition (RPSEMD) method to obtain a series of vibration signal component data IMFs, and select the optimal T-IMF using the transfer entropy (TE) method to obtain reorganized signal data; normalize the reorganized signal data to construct corresponding training sets and test sets.

[0135] The RPSEMD algorithm is used to iteratively generate sinusoids with different scales to handle all possible modal aliasing problems in different intrinsic modal components (IMFs). The general form of the sinusoid can be used as an auxiliary signal, which is expressed as:

[0136] s k (t|a k ,f k ,θ k )=ak cos(2πf k t+θ k ) (12)

[0137] Where: k is the number of IMFs obtained by decomposition; a k , f k ,θ k are amplitude, frequency and phase respectively. k (t) Designed according to the natural mode components selected after cluster analysis, it can ensure that RPSEMD can effectively solve the mode aliasing problem. k Moves k Changing the extreme point of (t) can not only help preserve more details of the separated mode aliasing problem, but also ensure that s k (t) is completely offset in the final result. The specific steps of the RPSEMD algorithm are as follows:

[0138] 1) Use EMD to process the signal x(t) and determine a based on the obtained intrinsic modal components k and f k θ ki Sampling between [0, 2π] with a phase shift of i (1≤i≤I), thus generating s(t|a k ,f k ,θ ki ).

[0139] 2) Obtain the first intrinsic mode component and average all I to obtain the last intrinsic mode component c k (t).

[0140] 3) Remove c from x(t) k (t), x(t)←x(t)-c k (t), let k=k+1.

[0141] 4) Repeat steps 1) to 3) until all the natural modal components are obtained, and the final x(t) is taken as the residual r(t).

[0142] The fault sequence component IMF obtained by RPSEMD algorithm decomposition is expressed as follows:

[0143] IMF=[imf1imf2…imf J ](13)

[0144] Wherein, J represents the number of components decomposed by the RPSEMD algorithm.

[0145] The transfer entropy (TE) method is used to reconstruct the characteristics of the fault sequence components IMF obtained by decomposition. The specific process is as follows:

[0146] Considering that the probability distribution p(x) in the original transfer entropy (TE) method is often realized by the kernel probability density estimation method, which often leads to a significant increase in the computational burden and thus greatly increases the computational time due to the increase in the analysis dimension, the present invention improves the kernel probability density estimation method and adopts the probability histogram method to realize this process. Its probability density can be defined as:

[0147]

[0148] In the formula, E represents r (r = 1, 2, 3) dimensional variables, and the interval of the variable e of each dimension [e min ,e max ] are divided into D equal parts, then the entire space is divided into D r The volume of each space is expressed as value δ represents the interval length of the δth variable; M represents the total number of samples; b α Indicates the number of samples in each small space.

[0149] Transfer entropy is a direct measure of the dynamic information transferred from one variable to another. Suppose that given two independent variables and They are s-order and k-order Markov processes respectively, and the transfer entropy is passed through the variable and variables Estimating variable x under historical data n+1 To evaluate the change in prediction effect. At time n+1, the entropy rate can be expressed as:

[0150]

[0151] In the formula, x n+1 Represents the sample of variable x at time n+1; Represents x n+1 and Probability of simultaneous occurrence; Indicates a given Under the condition x n+1 Probability of occurrence.

[0152] The joint conditional entropy rate is expressed as:

[0153]

[0154] In the formula, represents the joint probability density function; Indicated in and Under the condition x n+1According to equations (14) and (15), the information flow from variable y to variable x can be expressed as:

[0155]

[0156] In the formula, yes Under the condition x n+1 The probability density of . Since the transfer entropy is directional, the information transfer from variable x to variable y is also required. In summary, the transfer entropy from y to x is shown as:

[0157] TE y→x =T y→x -T x→y (18)

[0158] The fault sequence component T-IMF obtained after data reorganization after processing by transfer entropy is expressed as follows:

[0159] T-IMF = [imf1imf2…imf j ](j<J)(19)

[0160] Where j represents the number of fault sequence components obtained after reorganization by transfer entropy.

[0161] The obtained j fault sequence components are divided into training set and diagnosis test set in the ratio of 8:2, and the divided j data sets are used as input data of the subsequent diagnosis model.

[0162] Step 3: Figure 1 As shown in the figure, a bearing fault diagnosis model based on SageFormer is constructed. The specific implementation process is as follows:

[0163] Before inputting the bearing vibration sequence, learnable global tokens are added to each sequence, which are used to capture the global information of each sequence. The input multivariate time series (MTS) data is first converted into overlapping patches, and each patch is mapped to a unified latent space through a trainable linear projection. The global token embedding process of the bearing vibration sequence is as follows:

[0164] G(0) c ,:,:=[g1;...;g M ;x 1c E proj ;...;x Nc E proj ]+E pos (20)

[0165] Where G(0) is the embedding of the input bearing target sequence MTS sample, gi is the global token of the bearing sample vibration sequence, x ic is the i-th time step in the sampled vibration sequence c, E proj is a trainable linear projection of the bearing sample vibration sequence, E pos is the vibration sequence position embedding.

[0166] The bearing vibration sequence adjacency matrix is ​​learned in an end-to-end manner to capture the relationship between each sampled bearing vibration sequence, using randomly initialized node embeddings and transforming them through training parameters to learn unidirectional dependencies between sequences. Through a specific operation, the nearest neighbors of each node are determined and the weights of unconnected nodes are set to zero to generate the final sparse adjacency matrix.

[0167]

[0168] where E is the node embedding of the bearing vibration sequence; Θ1 and Θ2 are trainable parameters; act is the nonlinear activation function of the SageFormer encoding layer; M1 and M2 represent the node representations obtained by transforming the activation function and parameter matrix used to operate the vibration sequence, which are used to capture the relationship between the sequences; ReLU(·) represents the intrinsic activation function of the SageFormer model, which is used to calculate the relationship between M1 and M2; A′ is the learned adjacency matrix between bearing vibration sequences, which is used to represent the preliminary relationship between sequences.

[0169] The embedded tokens are processed through a SageFormer encoding layer, the global tokens of all sequences are collected and processed through a graph neural network (GNN) component to fuse the information of each bearing vibration sequence and its neighbors:

[0170]

[0171] Among them, G i is the global token of the i-th layer bearing vibration sequence; is the graph Laplace matrix of the bearing vibration sequence; W d is the graph convolution filter of the bearing vibration sequence, which is used to filter and aggregate the neighbor features of the node during the graph aggregation process; D represents the number of hops that the information propagates in the graph during the graph aggregation stage of the SageFormer model.

[0172] After passing through multiple encoding layers, the output of the model is passed through a decoding head to produce the final prediction. The decoding head uses the information received from the encoder to predict future time steps:

[0173]

[0174] In the formula, It is the characteristic representation of the bearing vibration sequence after being processed by the encoder, and Y is the output of the bearing vibration sequence diagnosis result.

[0175] Step 4: Figure 2 As shown, the sparse initialization and probability perturbation strategies are used to improve the initialization phase and population search phase of the non-monopolistic search algorithm NO, and an improved non-monopolistic search algorithm CNO is obtained.

[0176] The Monopoly Search (NO) algorithm is a new type of optimization algorithm based on local search, which aims to solve global optimization problems. The NO algorithm is characterized by its simplicity and its ability to balance exploration and exploitation. Through randomness, the algorithm is able to avoid falling into local optimal solutions in the search space and has a chance to find the global optimal solution. In addition, the NO algorithm is a single-solution algorithm, which means that it only processes one candidate solution in each iteration. The NO search algorithm has two main phases, namely the exploration phase and the exploitation phase. In these phases, the NO algorithm modifies the candidate solutions in each dimension (decision variable).

[0177] Exploration phase: The sparse initialization method is introduced into the exploration phase of the algorithm. At the beginning of the algorithm, the initial position of the population is determined using random numbers. The algorithm explores the solution space based on the initial parameters obtained by the sparse initialization method to find new potential optimal solutions. The exploration phase is executed when the iteration counter (t) is less than or equal to half of the total number of iterations (T). The formula is as follows:

[0178] I i =rand(1,D)⊙(I max -I min )+I min (25)

[0179] Anew(j)=rand*A(RP)+I i (26)

[0180] Among them, I i is the position of the two individuals of the i-th SageFormer model learning rate lr and batch size hu, I max and I min are the upper and lower bounds of the learning rate lr and batch size hu in the search space, rand(1, D) is a vector of dimension D of uniformly distributed random numbers in the interval 0, 1, Anew(j) is the jth dimension in the newly generated solution of learning rate lr and batch size hu, rand is a random value between 0 and 1. j is the dimension index defined between 1 and the maximum dimension of the search space, A is the current solution, A(RP) is the selected random dimension of the current solution of learning rate lr and batch size hu. RP is a random number between 1 and D.

[0181] Development phase: The algorithm performs local search in the neighborhood of the current solution to refine the quality of the solution. The development phase is executed when the iteration count is higher than half of T. This phase aims to utilize the neighborhood around the current solution. The current solution is randomly perturbed to search for a better position around it. When the algorithm gradually converges, it will be difficult to jump out once it is trapped in the local optimum. Therefore, the present invention introduces a probabilistic perturbation strategy into the algorithm population update process to help the algorithm jump out of the local optimum. The specific implementation process is as follows:

[0182] Anew(j)=A(j)-(A(SRP)*rand)*eps-(A(j)-NO) (27)

[0183]

[0184] Among them, P r is the disturbance probability; D a is the dimension of the algorithm optimization problem and D a =2; T a is the maximum number of iterations. In the early stages of algorithm iteration, the perturbation probability is small, so that the algorithm can quickly find the global optimal solution; in the later stages of algorithm iteration, the perturbation probability increases, which enhances population diversity. Generally speaking, T a is a large constant. When the algorithm iterates to the later stage t a →T a , disturbance probability:

[0185]

[0186] After adding the probability perturbation, the formula for updating the position of the algorithm particles is optimized as follows:

[0187] Anew(j)=A(j)-(A(SRP)*rand1)*eps-(A(j)-NO); rand1<P r (30)

[0188] where A(j) is the jth position in the current solution, rand1 is a random value between [0, 1], A(SRP) is a dimension randomly selected from the current solution dimension, eps is a very small value, SRP is a solution random dimension, a random number between 1 and the dimension size; NO is a tuning parameter used to adjust the search process of the proposed method.

[0189] Step 5: Use the improved non-monopoly search algorithm to optimize the hyperparameters of the bearing fault diagnosis model based on SageFormer, where the hyperparameters are the learning rate lr and the batch size hu. The specific implementation process is as follows:

[0190] The matching function of the diagnosis results is used as the objective function of the SageFormer fault diagnosis model;

[0191] R = sum(S = S P ) / length(S)(31)

[0192] Among them, R represents the accuracy of bearing fault diagnosis, S represents the actual result of the bearing status, and S P Indicates the results of the diagnosis of the bearing condition.

[0193] Determine the learning rate lr and batch size hu of the population type SageFormer model, as follows:

[0194]

[0195] Among them, X represents the population of the NO search algorithm, and d represents the number of particles of one type in the population.

[0196] The upper and lower limits of the search space are determined by the range of the learning rate lr and the batch size hu; the specific expression is as follows:

[0197]

[0198] Step 6. Use the bearing vibration time series data processed in step 2 and the improved non-monopolistic search algorithm to train the bearing fault diagnosis model established by SG-RPSEMD-TE-CNO-SageFormer, find the optimal parameters of the network, and minimize the loss function error; use the trained and optimized bearing fault diagnosis model to perform online fault diagnosis on the bearing vibration data and obtain the corresponding diagnosis results.

[0199] The training set and test set data constructed after feature reorganization and normalization obtained in step 2 are used as the input of the SG-RPSEMD-TE-CVNO-SageFormer bearing fault diagnosis model, the model is trained, and the trained SG-RPSEMD-TE-CVNO-SageFormer bearing fault diagnosis model is obtained. After being processed in step 2, the various sensor data collected in real time are sent to the trained SG-RPSEMD-TE-CVNO-SageFormer fault diagnosis model for fault type diagnosis and identification.

[0200] The above embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable people familiar with the technology to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent transformation or modification made according to the spirit of the present invention should be included in the protection scope of the present invention.

Claims

1. A bearing fault diagnosis method, characterized in that: The following steps are involved: (1) SG filtering is used to denoise the vibration signals of the bearing at different positions under healthy and various fault conditions acquired in advance; and feature fusion is performed on the denoised multidimensional data; (2) The denoised bearing vibration signal is decomposed using the regenerative phase-shifted sine-assisted empirical mode decomposition algorithm RPSEMD to obtain a series of vibration signal component data IMFs; the optimal T-IMF is selected using the transfer entropy TE to obtain the reorganized signal data; the reorganized signal data is normalized to construct the corresponding training set and test set; (3) Construct a bearing fault diagnosis model based on SageFormer; (4) The sparse initialization and probability perturbation strategies are used to improve the initialization phase and population search phase of the non-monopolistic search algorithm NO, and the improved non-monopolistic search algorithm CNO is obtained; (5) Using the improved non-monopolistic search algorithm CNO to optimize the hyperparameters of the bearing fault diagnosis model based on SageFormer; (6) The training set is used to train the bearing fault diagnosis model, the optimal parameters of the model are obtained, and the loss function error is minimized; the bearing fault diagnosis model after training and optimization is used to perform online fault diagnosis on the bearing vibration data to obtain the corresponding diagnosis results.

2. A bearing fault diagnosis method according to claim 1, characterized in that: The process of using SG filtering to denoise the vibration signals of the pre-acquired healthy and various fault states of different positions of the bearing in step (1) is as follows: First, a polynomial fit is performed on a set of 2H+1 sample center discrete points centered at m=0: Among them, T is the polynomial order, k represents the order of the discrete sequence value at this point, and z k is the polynomial coefficient, m is the discrete sequence value of the point, and y(m) is the fitting value of the point; Secondly, calculate the least squares fitting residuals and E of all fitted data and original data in the filter window N : Among them, x p (m) is the pth group of samples of the original bearing vibration signal, and H represents a group of bearing vibration signal sample boundaries; When E N When it is minimum, the fitting effect is the best; take the partial derivative of formula (2) and get the polynomial fitting coefficient z k , calculate the fitting result and assign it to the midpoint of the window as the fitting result output, that is: f[0]=y[0] (3) Shift the filter along the time axis in a single step, repeat the above steps until the end of the data, and realize SG moving smoothing; use SG convolution smoothing to replace the moving smoothing algorithm to expand the denoising dimension of the traditional SG filter, improve the calculation efficiency and fitting results, and the specific implementation process is as shown in formula (4): f i [m] represents the vibration signal collected by the i-th (i < w) group of sensors after noise reduction; after filtering and noise reduction of the w groups of local vibration signals collected by w sensors, we get: F[m]=[f1[m] f2[m] … f w [m]] (5)。 3. A bearing fault diagnosis method according to claim 1, characterized in that: The process of implementing feature fusion on the denoised multidimensional data in step (1) is as follows: Any two signals f α [m], f β [m] the correlation between α, β∈[1,w] Defined as: Where t=0, 1, 2, 3, ..., N-1 is the time series of the signal; signal f α [m] and other signals f β The total correlation energy of [m] is: The energy of n signals is normalized separately, that is: After energy normalization, the signal The N discrete data points are Variance contribution rate K il Defined as: In the formula, For signal The first data point in The distribution coefficient k at the first data point il for: Based on the obtained distribution coefficient, w signals are fused into one signal, and the fused signal R(l) is:

4. A bearing fault diagnosis method according to claim 1, characterized in that: The implementation process of decomposing the denoised bearing vibration signal using the regenerative phase-shifted sine-assisted empirical mode decomposition algorithm RPSEMD in step (2) is as follows: S1: Use EMD to process the signal x(t) and determine a based on the obtained intrinsic mode components k and f k θ ki Sampling between [0, 2π] with a phase shift of i (1≤i≤I), thus generating s(t|a k ,f k ,θ ki ); S2: Get the first intrinsic mode component and average all I to get the last intrinsic mode component c k (t); S3: Remove c from x(t) k (t), x(t)←x(t)-c k (t), let k = k + 1; S4: Repeat S1 to S3 until all the intrinsic modal components are obtained, and the final x(t) is used as the residual r(t); the fault sequence component IMF obtained by the RPSEMD algorithm is expressed as follows: IMF=[imf1 imf2…imf J ](13)Wherein, J represents the number of components decomposed by the RPSEMD algorithm.

5. A bearing fault diagnosis method according to claim 1, characterized in that: The specific process of using the transfer entropy TE to select the optimal T-IMF in step (2) to obtain the reorganized signal data is as follows: The probability histogram is used to replace the kernel probability density estimation method, and the probability density is defined as: In the formula, E represents r (r = 1, 2, 3) dimensional variables, and the interval of the variable e of each dimension [e min ,e max ] are divided into D equal parts, then the entire space is divided into D r The volume of each space is expressed as value δ represents the interval length of the δth variable; M represents the total number of samples; b α Indicates the number of samples in each small space; Transfer entropy is a direct measure of the dynamic information transferred from one variable to another. Suppose that given two independent variables and They are s-order and k-order Markov processes respectively, and the transfer entropy is passed through the variable and variables Estimating variable x under historical data n+1 To evaluate the change in prediction effect; at time n+1, the entropy rate can be expressed as: In the formula, x n+1 Represents the sample of variable x at time n+1; Represents x n+1 and Probability of simultaneous occurrence; Indicates a given Under the condition x n+1 Probability of occurrence; The joint conditional entropy rate is expressed as: In the formula, represents the joint probability density function; Indicated in and Under the condition x n+1 The probability density of ; According to equations (15) and (16), the information flow from variable y to variable x can be expressed as: In the formula, yes Under the condition x n+1 The probability density of ; Since the transfer entropy is directional, information transfer from variable x to variable y is also required; in summary, the transfer entropy from y to x is shown as: TE y→x =T y→x -T x→y (18) After the data is reorganized by the transfer entropy, the fault sequence component T-IMF is expressed as follows: T-IMF = [imf1 imf2…imf j ](j<J)(19)wherein j represents the number of fault sequence components obtained after reorganization by transfer entropy.

6. A bearing fault diagnosis method according to claim 1, characterized in that: The implementation process of step (3) is as follows: Before inputting the bearing vibration sequence, learnable global tokens are added to each sequence, which are used to capture the global information of each sequence; the input multivariate time series data is first converted into overlapping blocks, and each block is mapped to a unified latent space through a trainable linear projection. The global token embedding process of the bearing vibration sequence is as follows: G(0) c ,:,:=[g1;...;g M ;x 1c E proj ;...;x Nc E proj ]+E pos (20) where G(0) is the embedding of the input bearing target sequence MTS sample, g i is the global token of the bearing sample vibration sequence, x ic is the i-th time step in the sampled vibration sequence c, E proj is a trainable linear projection of the bearing sample vibration sequence, E pos is the vibration sequence position embedding; The bearing vibration sequence adjacency matrix is ​​learned in an end-to-end manner to capture the relationship between each sampled bearing vibration sequence. Randomly initialized node embeddings are used and transformed through training parameters to learn unidirectional dependencies between sequences. Through a specific operation, the nearest neighbors of each node are determined and the weights of unconnected nodes are set to zero to generate the final sparse adjacency matrix: Where E is the node embedding of the bearing vibration sequence; Θ1 and Θ2 are trainable parameters; act is the nonlinear activation function of the SageFormer encoding layer; M1 and M2 represent the node representations obtained by transforming the activation function and parameter matrix used to operate the vibration sequence, which is used to capture the relationship between the sequences; ReLU(·) represents the intrinsic activation function of the SageFormer model, which is used to calculate the relationship between M1 and M2; A′ is the learned adjacency matrix between bearing vibration sequences, which is used to represent the preliminary relationship between sequences; The embedded tokens are processed through a SageFormer encoding layer, which collects the global tokens of all sequences and processes them through a graph neural network component to fuse the information of each bearing vibration sequence and its neighbors: Among them, G i is the global token of the i-th bearing vibration sequence; is the graph Laplace matrix of the bearing vibration sequence; W d is the graph convolution filter of the bearing vibration sequence, which is used to filter and aggregate the neighbor features of the node during the graph aggregation process; D represents the number of hops that the information propagates in the graph during the graph aggregation stage of the SageFormer model; After passing through multiple encoding layers, the output of the model is passed through a decoding head to produce the final prediction; the decoding head uses the information received from the encoder to predict future time steps: In the formula, It is the characteristic representation of the bearing vibration sequence after being processed by the encoder, and Y is the output of the bearing vibration sequence diagnosis result.

7. A bearing fault diagnosis method according to claim 1, characterized in that: The implementation process of step (4) is as follows: Exploration phase: The sparse initialization method is introduced into the exploration phase of the monopoly search algorithm NO. At the beginning of the algorithm, the initial position of the population is determined using random numbers. The algorithm explores the solution space based on the initial parameters obtained by the sparse initialization method to find new potential optimal solutions; the exploration phase is executed when the iteration counter (t) is less than or equal to half of the total number of iterations (T). The formula is as follows: I i =rand(1,D)⊙(I max -I min )+I min (25) Anew(j)=rand*A(RP)+I i (26) Among them, I i is the position of the two individuals of the i-th SageFormer model learning rate lr and batch size hu, I max and I min are the upper and lower bounds of the learning rate lr and batch size hu in the search space, rand(1, D) is a vector of dimension D of uniformly distributed random numbers in the interval 0, 1, Anew(j) is the jth dimension in the newly generated solution of learning rate lr and batch size hu, and rand is a random value between 0 and 1; j is a dimension index defined between 1 and the maximum dimension of the search space, A is the current solution, A(RP) is the chosen random dimension of the solution with the current learning rate lr and batch size hu; RP is a random number between 1 and D; Development stage: Introduce the probability perturbation strategy into the population update process of the NO algorithm to help the algorithm escape from the local optimum; the specific implementation process is as follows: Anew(j)=A(j)-(A(SRP)*rand)*eps-(A(j)-NO)(27) Where: P r is the disturbance probability; D a is the dimension of the algorithm optimization problem and D a =2; T a is the maximum number of iterations; in the early stage of algorithm iteration, the perturbation probability is small, so that the algorithm can quickly find the global optimal solution; in the later stage of algorithm iteration, the perturbation probability increases, which enhances the diversity of the population; T a is a large constant. When the algorithm iterates to the later stage t a →T a , disturbance probability: After adding the probability perturbation, the formula for updating the position of the algorithm particles is optimized as follows: Anew(j)=A(j)-(A(SRP)*rand1)*eps-(A(j)-NO); rand1<P r (30) Among them, A(j) is the j-th position of the learning rate lr and batch size hu in the current solution, rand1 is a random value between [0, 1], A(SRP) is the dimension of the learning rate lr and batch size hu randomly selected from the current solution dimension, eps is a very small value, SRP is a solution random dimension, a random number between 1 and the dimension size; NO is a tuning parameter used to adjust the search process of the proposed method.

8. A bearing fault diagnosis method according to claim 1, characterized in that: The implementation process of step (5) is as follows: The matching function of the diagnosis results is used as the objective function of the SageFormer fault diagnosis model; R = sum(S = S P ) / length(S)(31) Where R represents the accuracy of bearing fault diagnosis, S represents the actual result of the bearing status, and S P Results indicating the diagnosed bearing condition; Determine the learning rate lr and batch size hu of the population type SageFormer model, as follows: Among them, X represents the population of NO search algorithm, and d represents the number of particles of one type in the population; The upper and lower limits of the search space are determined by the range of the learning rate lr and the batch size hu; the specific expression is as follows:

Images