A method for diagnosing bearing faults in petrochemical units

Abstract

The invention discloses a method for diagnosing bearing faults of a petrochemical unit, comprising the following steps: extracting vibration signals of rolling bearings of petrochemical units; improving the hyperbolic sine-cosine optimization algorithm SCHO by using random inertia weights and Levy flight strategies to obtain an ISCHO algorithm; using ISCHO to optimize the parameters of time-varying filter empirical mode decomposition TVFEMD, and using the optimized TVFEMD to adaptively decompose the bearing fault signal; obtaining intrinsic mode parameters IMF and extracting features thereof to construct feature vectors; using random forests to screen feature vectors; establishing and training a MSCNN‑LTransformer fault diagnosis model; and finally verifying the fault diagnosis accuracy of the model. The invention can effectively extract the fault information of rolling bearings of petrochemical units, analyze the causes of the faults, and improve the accuracy of bearing fault diagnosis of petrochemical units.

Description

Technical Field

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

[0002] The petrochemical industry is an important basic industry in my country, and plays an important role in improving the national economy, ensuring the supply of industrial foundations and promoting agricultural development. With the increasing development of the petrochemical industry, petrochemical units have become more and more important in production practice. Once the equipment fails, it may cause damage to the machine, and in severe cases, it may cause casualties of workers, resulting in serious consequences. Bearings are one of the most important components in petrochemical units. Because bearings have the advantages of convenient loading and unloading and lubrication, they are widely used in petrochemical units. When a bearing fails, it will cause the running machine to fail. Therefore, fault diagnosis of petrochemical unit bearings can detect faults at an early stage and avoid subsequent greater damage by timely maintenance.

[0003] Since petrochemical units are mostly in a complex working environment of high temperature and high pressure, there are strong noise and other interference in the extracted bearing vibration signal. In the past, most of the petrochemical unit bearing fault diagnosis methods used traditional machine learning for fault diagnosis, but the fault diagnosis model using traditional machine learning cannot well separate the noise and other interference information in the signal. The fault diagnosis model using traditional machine learning also has problems such as long running time and large amount of calculation.

[0004] Compared with the prior art, the petrochemical unit bearing fault diagnosis method designed by the present invention uses the time-varying filter empirical mode decomposition TVFEMD optimized by the improved hyperbolic sine-cosine optimization algorithm ISCHO to decompose the vibration signal, and selects the intrinsic mode parameter IMF with more fault feature information for calculation, which can effectively remove the noise information in the signal and speed up the signal decomposition. Considering that MSCNN can extract features from time series with different spans, LTransformer can capture the global features of time series. Therefore, the MSCNN-LTransformer fault diagnosis model is established to diagnose the bearing fault of the petrochemical unit, which can effectively extract the fault information in the vibration signal, analyze the cause of the fault, realize the fault diagnosis of the petrochemical unit bearing, and improve the fault diagnosis accuracy of the petrochemical unit. Summary of the invention

[0005] Purpose of the invention: The present invention aims to provide a petrochemical unit bearing fault diagnosis method with high fault diagnosis accuracy and fast running speed, by using the time-varying filter empirical mode decomposition TVFEMD optimized by the improved hyperbolic sine-cosine optimization algorithm ISCHO combined with the multi-scale neural network MSCNN and LinearTransformer (LTransformer), to realize the fault diagnosis of petrochemical unit bearings.

[0006] Technical solution: The present invention discloses a method for diagnosing bearing faults of a petrochemical unit, comprising the following steps:

[0007] Step 1: Collect vibration signals of four different rolling bearing conditions: normal condition, inner ring fault, outer ring fault, and rolling element fault at the same speed.

[0008] Step 2: Use random inertia weights and Levy flight strategy to improve the hyperbolic sine-cosine optimization algorithm SCHO, speed up the algorithm, avoid local optimal situations in the algorithm, and obtain the improved hyperbolic sine-cosine optimization algorithm ISCHO.

[0009] Step 3: Use ISCHO to optimize the parameters of the time-varying filter empirical mode decomposition TVFEMD, which include the bandwidth threshold ξ and the B-spline order. Use the optimized TVFEMD to decompose the collected bearing vibration signal to obtain multiple intrinsic mode function components IMF, and further select the IMF with more fault characteristics.

[0010] Step 4: Use the sliding window method to perform overlapping sampling on the selected IMF and divide it into multiple sample sets. Calculate the time domain features, frequency domain features, energy entropy features, sample entropy features and singular value features of the vibration signal to construct a feature vector.

[0011] Step 5: Use random forest to screen the feature vectors, obtain the optimal feature vector, and divide the sample set into training set and test set in a ratio of 7:3.

[0012] Step 6: Construct and train the MSCNN-LTransformer fault diagnosis model by fusing the multi-scale convolutional neural network MSCNN and LinearTransformer (LTransformer). MSCNN extracts features from time series of different spans, and LTransformer captures the global features of the time series.

[0013] Step 7: Use the test set data to test the trained MSCNN-LTransformer fault diagnosis model.

[0014] Furthermore, the improved hyperbolic sine-cosine optimization algorithm ISCHO in step 2 comprises the following steps:

[0015] Step 2.1: Since the second stage of the hyperbolic sine-cosine optimization algorithm SCHO exploration is to perform an undirected search at the current position, a random inertia weight strategy is introduced to avoid the local optimal situation. The improved formula is as follows:

[0016]

[0017] w=W max +r 14 (W max -W min )+σ·R

[0018] Where t is the current iteration number; is the jth position of the ith solution in the current iteration, i = 1, 2, ... N, N is the number of candidate solutions, j = 1, 2, ... dim, dim is the dimension of the problem; ε represents a very small positive number, which is 0.003 in this paper; is the jth position of the current optimal solution; W2 is The weight coefficient; r is a random number between the interval [0,1]; w is the random inertia weight; W max is the maximum value of random inertia weight, which is taken as 0.9 in this paper; W min is the minimum value of the random inertia weight, which is taken as 0.4 in this paper; σ is the deviation between the random inertia weight and the mathematical expectation, which is taken as 0.3 in this paper; R is a random number of normal distribution.

[0019] Step 2.2: Since the second stage of SCHO development is to conduct in-depth development around the optimal solution, the Levy flight strategy is introduced to speed up the algorithm. The improved formula is as follows:

[0020]

[0021] Where l is the step size scaling factor; levy(λ) is the random path of Levy flight.

[0022] Further, the optimization of TVFEMD using ISCHO in step 3 includes the following steps:

[0023] Step 3.1: Initialize the relevant parameters of ISCHO, including population size, problem dimension dim, number of iterations t, and the switching value T between the first and second stages of exploration and development. The calculation formula of T is as follows:

[0024]

[0025] Among them, floor is the round-down function; Max_iteration represents the maximum value of iteration; ct is the coefficient of two-stage switching, which is 3.6 in this paper.

[0026] Step 3.2: Randomly initialize a set of candidate solutions. The candidate solution calculation formula is as follows:

[0027] X=rand(N,dim)×(ub-lb)+lb

[0028] Among them, X is a set of random candidate solutions; rand is a random number in [0,1]; ub is the upper bound of the variable; lb is the lower bound of the variable.

[0029] Step 3.3: Calculate the objective function value of the candidate solution and determine the current optimal solution. The objective function used in this paper is the minimum envelope entropy. The calculation formula of the minimum envelope entropy Ep is:

[0030]

[0031] Where a(h) is the envelope signal, h = 1, 2, ... N; p h is the normalized form of a(h).

[0032] Step 3.4: Calculate the value of A, which is the value converted by the Sinh and Cosh conversion mechanisms in the exploration and development stages; when A is greater than 1, it enters the exploration stage, and when A is less than 1, it enters the development stage. The calculation formula of A is as follows:

[0033]

[0034] Among them, p is the balance coefficient of the exploration stage in the iterative process, which is taken as 10 in this paper; q is the balance coefficient of the development stage in the iterative process, which is taken as 9 in this paper; r is a random number in the interval [0,1].

[0035] If t = BS k When , find out the current position of the second solution and update BS k and search space, redistributing candidate solutions.

[0036] BS k The calculation formula is as follows:

[0037]

[0038] Among them, BS k is the number of iterations of the current bounded search strategy; α is the sensitivity coefficient value, which is taken as 4.6 in this paper.

[0039] The calculation formula of the search space is as follows:

[0040]

[0041]

[0042] Among them, ub kis the upper bound of the search space; lb k is the lower bound of the search space; is the jth position of the current optimal solution; represents the jth position of the suboptimal solution.

[0043] Step 3.5: When A>1, enter the exploration stage. Calculate the values ​​of W1 and W2. The calculation formulas of W1 and W2 are as follows:

[0044] W1=r3×a1×(coshr4+u×sinhr4-1)

[0045] W2=r6×a2

[0046] Among them, W1 is The weight coefficient of The weight coefficient of ; a1, a2 are monotonically decreasing functions; u is the sensitivity coefficient for controlling the exploration accuracy in the first stage, which is taken as 0.388 in this paper.

[0047] When t≤T, the first stage of exploration begins. The position update formula of the candidate solution in the first stage of exploration is as follows:

[0048]

[0049] When t ≥ T, the second stage of exploration begins. The position update formula of the candidate solution in the second stage of exploration is as follows:

[0050]

[0051] w=W max +r 14 (W max -W min )+σ·R

[0052] When A<1, it enters the development stage and calculates the weight coefficient W3 of the first stage of development. The formula is as follows:

[0053] W3=r9×a1×(coshr 10 +u×sinhr 10 )

[0054] When t≤T, the first stage of development begins. The position update formula of the candidate solution in the first stage of development is as follows:

[0055]

[0056] When t≥T, the second stage of development begins. The position update formula of the candidate solution in the second stage of development is as follows:

[0057]

[0058] Among them, W2 controls the extent of the second stage of development.

[0059] Step 3.6: Determine whether the termination condition is met. If it is met, output the optimal solution. If not, repeat steps 3.3 to 3.5.

[0060] Step 3.7: Use TVFEMD to decompose the vibration signal. The calculation formula is:

[0061]

[0062] Among them, θ(n) is the cutoff criterion; B Loughlin is the Loughlin instantaneous bandwidth; is the weighted mean instantaneous frequency.

[0063] Furthermore, the step 4 of calculating the time domain features, frequency domain features, energy entropy features, sample entropy features and singular value features of the vibration signal includes the following steps:

[0064] Step 4.1: Calculating the time domain characteristics of the vibration signal includes calculating the peak value Y1, average value Y2, and root mean square value Y3 of the vibration signal. The calculation formulas for the peak value Y1, average value Y2, and root mean square value Y3 are:

[0065] Y1=max(x(n))

[0066]

[0067]

[0068] Wherein, x(n) is the original signal, n=1,2,3,…N, and N is the data length.

[0069] Step 4.2: Calculating the frequency domain characteristics of the vibration signal includes calculating the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 of the vibration signal. The calculation formulas for the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 are:

[0070]

[0071]

[0072]

[0073] Where, the spectrum sequence of s(k) is k=1,2,3,…K, where K is the spectrum line; f k is the frequency value of the kth spectrum line.

[0074] Step 4.3: Energy entropy feature EF i The calculation formula is:

[0075]

[0076] Among them, ES i is the energy value of the IMF component, and E represents the sum of the IMF components.

[0077] The calculation formula of sample entropy feature SampEn is:

[0078]

[0079] Where m is the pattern dimension; r is the similarity tolerance threshold; B m (r) is the probability that two sequences match m points under r, A m (r) is the probability that the two sequences match m+1 points.

[0080] Singular Value Feature SVF i The calculation formula is

[0081]

[0082] Among them, A represents the matrix constructed by IMF components, U m×m and V n×n are the left and right singular matrices of A, Λ is an m×n diagonal matrix, λ1,λ2,…,λ i are the singular values ​​of the matrix A.

[0083] Furthermore, the MSCNN-LTransformer fault diagnosis model in step 6 includes the following steps:

[0084] Step 6.1: Input the training set data into the MSCNN layer. The MSCNN layer is composed of convolution kernels of different sizes. Convolution kernels of different sizes can obtain different feature sizes, so MSCNN can obtain richer feature information. This article uses convolution kernel sizes of 1, 3, 7, and 9. The operation of CNN can be expressed as:

[0085] S1=C1(x)

[0086] S2=C3(S1)

[0087] S3=C7(S2)

[0088] S4=C9(S3)

[0089] S=Concat(S2,S3,S4)

[0090] Among them, C n represents a one-dimensional convolutional element with a convolution kernel size of n, where n = 1, 3, 7, 9; x represents the input; S represents the output state of the MCNN layer; Concat(·) represents the concatenation operation.

[0091] Step 6.2: Encode the local feature information extracted from step 6.1 and input it into LTransformer. LTransformer is an improved model of Transformer, which mainly consists of input layer, feature encoder layer, feature decoder layer and output layer. The feature encoder and feature decoder are mainly composed of multi-head attention module MSA and feedforward neural network module FFN. Each module uses residual connection Add and then performs layer normalization Norm.

[0092] LTransformer transforms the traditional softmax attention into dot product attention based on feature mapping, which reduces the computational cost and complexity of the algorithm. i 'The calculation formula is as follows:

[0093]

[0094] Among them, sim(·) is the similarity function; Q is the query vector; K is the key vector; V is the value vector.

[0095] MSA can learn features from different feature subspaces and has good global information extraction capabilities. The calculation formula of MSA is as follows:

[0096] M(Q,K,V)=Concat(H1,H2,…H h )·W

[0097] H i =V i '(Q i ,K i ,V i ),i=1,2…h

[0098] Among them, M(·) is the multi-head attention mechanism operation, H i is the i-th head of the attention mechanism, W is the multi-head attention weight matrix, V i 'For Q i The cumulative value of .

[0099] FFN consists of two fully connected layers. The calculation formula of FFN is as follows:

[0100] FFN(x)=max(0,xW1+b1)W2+b2

[0101] Among them, x is the normalized output matrix of the attention layer; W1, W2 are weight vectors; and b is the bias term.

[0102] A residual connection Add is used between each submodule of the feature encoder and feature decoder to prevent network degradation, and then the layer normalization Norm is used to normalize the activation value of each layer to avoid overfitting and speed up the training. The calculation formula is:

[0103] X B =LN(M(X A )+X A )

[0104] X C =LN(FFN(X B )+X B )

[0105] Among them, X A is the input matrix of the encoding layer; X B Input matrix for FFN layer; X C is the output matrix of the coding layer; LN is the layer normalization process.

[0106] Step 6.3: Input the feature encoder output result into the linear layer and the Softmax classifier layer and then output the fault diagnosis result.

[0107] Beneficial effects:

[0108] The present invention uses TVFEMD optimized by ISCHO to decompose the vibration signal, effectively removes the noise information in the vibration signal, and reduces the noise interference in the signal. An MSCNN-LTransformer fault diagnosis model is established, in which MSCNN extracts features from time series of different spans, and LTransformer captures the global features of the time series. The fault information in the vibration signal is effectively extracted, the cause of the fault is analyzed, the fault diagnosis of the petrochemical unit bearing is realized, and the accuracy of the fault diagnosis of the petrochemical unit is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0109] Figure 1 is a flow chart of the present invention;

[0110] Figure 2 To improve the hyperbolic sine-cosine optimization algorithm flow chart;

[0111] Figure 3 A schematic diagram of a bearing signal decomposed by TVFEMD provided by the present invention;

[0112] Figure 4 This is the MSCNN-LTransformer fault diagnosis model diagram;

[0113] Figure 5 This is a comparison chart of the fault diagnosis accuracy of the present invention and other models. DETAILED DESCRIPTION

[0114] The present invention is described in detail below with reference to specific examples.

[0115] The present invention discloses a method for diagnosing bearing faults of a petrochemical unit. Figures 1 to 4 As shown, the specific steps include:

[0116] Step 1: Collect vibration signals of four different rolling bearing conditions: normal condition, inner ring fault, outer ring fault, and rolling element fault at the same speed.

[0117] Step 2: Use random inertia weights and Levy flight strategy to improve the hyperbolic sine-cosine optimization algorithm SCHO, speed up the algorithm, avoid local optimal situations in the algorithm, and obtain the improved hyperbolic sine-cosine optimization algorithm ISCHO.

[0118] Step 2 specifically includes the following steps:

[0119] Step 2.1: Since the second stage of SCHO exploration is to perform an undirected search at the current position, a random inertia weight w is introduced to avoid the local optimal situation. The improved formula is as follows:

[0120]

[0121] w=W max +r 14 (W max -W min )+σ·R

[0122] Where t is the current iteration number; is the jth position of the ith solution in the current iteration, i = 1, 2, ... N, N is the number of candidate solutions, j = 1, 2, ... dim, dim is the dimension of the problem; ε represents a very small positive number, which is 0.003 in this paper; is the jth position of the current optimal solution; W2 is The weight coefficient; r is a random number between the interval [0,1]; w is the random inertia weight; W max is the maximum value of random inertia weight, which is taken as 0.9 in this paper; W min is the minimum value of the random inertia weight, which is taken as 0.4 in this paper; σ is the deviation between the random inertia weight and the mathematical expectation, which is taken as 0.3 in this paper; R is a random number of normal distribution.

[0123] Step 2.2: Since the second stage of SCHO development is to conduct in-depth development around the optimal solution, the Levy flight strategy is introduced to speed up the algorithm. The improved formula is as follows:

[0124]

[0125] Where l is the step size scaling factor; levy(λ) is the random path of Levy flight.

[0126] Step 3: Use ISCHO to optimize the parameters of the time-varying filter empirical mode decomposition TVFEMD, which include the bandwidth threshold ξ and the B-spline order. Use the optimized TVFEMD to decompose the collected bearing vibration signal to obtain multiple intrinsic mode function components IMF, and further select the IMF with more fault characteristics.

[0127] Step 3 specifically includes the following steps:

[0128] Step 3.1: Initialize the relevant parameters of ISCHO, including population size, problem dimension dim, switching coefficient ct between the two stages, sensitivity coefficient α, balance coefficient p in the exploration stage, balance coefficient q in the development stage, number of iterations t, and switching value T between the first and second stages of exploration and development. The calculation formula of T is as follows:

[0129]

[0130] Among them, floor is the round-down function; Max_iteration represents the maximum value of iteration; ct is the coefficient of two-stage switching, which is 3.6 in this paper.

[0131] Step 3.2: Randomly initialize a set of candidate solutions. The candidate solution calculation formula is as follows:

[0132] X=rand(N,dim)×(ub-lb)+lb

[0133] Among them, X is a set of random candidate solutions; rand is a random number in [0,1]; ub is the upper bound of the variable; lb is the lower bound of the variable.

[0134] Step 3.3: Calculate the objective function value of the candidate solution and determine the current optimal solution. The objective function used in this paper is the minimum envelope entropy. The calculation formula of the minimum envelope entropy Ep is:

[0135]

[0136] Where a(h) is the envelope signal, h = 1, 2, ... N; p h is the normalized form of a(h).

[0137] Step 3.4: Calculate the value of A, which is the value converted by the Sinh and Cosh conversion mechanisms in the exploration and development stages; when A is greater than 1, it enters the exploration stage, and when A is less than 1, it enters the development stage. The calculation formula of A is as follows:

[0138]

[0139] Among them, p is the balance coefficient of the exploration stage in the iterative process, which is taken as 10 in this paper; q is the balance coefficient of the development stage in the iterative process, which is taken as 9 in this paper; r is a random number in the interval [0,1].

[0140] If t = BS k When , find out the current position of the second solution and update BS k and search space, redistributing candidate solutions.

[0141] BS k The calculation formula is as follows:

[0142]

[0143] Among them, BS k is the number of iterations of the current bounded search strategy; α is the sensitivity coefficient value, which is taken as 4.6 in this paper.

[0144] The calculation formula of the search space is as follows:

[0145]

[0146]

[0147] Among them, ub k is the upper bound of the search space; lb k is the lower bound of the search space; is the jth position of the current optimal solution; represents the jth position of the suboptimal solution.

[0148] Step 3.5: When A>1, enter the exploration stage. Calculate the values ​​of W1 and W2. The calculation formulas of W1 and W2 are as follows:

[0149] W1=r3×a1×(coshr4+u×sinhr4-1)

[0150]

[0151] Among them, W1 is The weight coefficient of ; u is the sensitivity coefficient for controlling the exploration accuracy in the first stage, which is taken as 0.388 in this paper; a1 is a monotonically decreasing function; m is the sensitivity coefficient for controlling the exploration accuracy, which is taken as 0.45 in this paper.

[0152] W2=r6×a2

[0153]

[0154] Among them, W2 is The weight coefficient of ; a2 is a monotonically decreasing function; n is the sensitivity coefficient for controlling the exploration accuracy in the second stage, which is taken as 0.5 in this paper.

[0155] When t≤T, the first stage of exploration begins. The position update formula of the candidate solution in the first stage of exploration is as follows:

[0156]

[0157] When t ≥ T, the second stage of exploration begins. The position update formula of the candidate solution in the second stage of exploration is as follows:

[0158]

[0159] w=W max +r 14 (W max -W min )+σ·R

[0160] When A<1, it enters the development stage and calculates the weight coefficient W3 of the first stage of development. The formula is as follows:

[0161] W3=r9×a1×(coshr 10 +u×sinhr 10 )

[0162] When t≤T, the first stage of development begins. The position update formula of the candidate solution in the first stage of development is as follows:

[0163]

[0164] When t≥T, the second stage of development begins. The position update formula of the candidate solution in the second stage of development is as follows:

[0165]

[0166] Among them, W2 controls the extent of the second stage of development.

[0167] Step 3.6: Determine whether the termination condition is met. If it is met, output the optimal solution. If not, repeat steps 3.3 to 3.5.

[0168] Step 3.7: Use TVFEMD to decompose the vibration signal. The calculation formula is:

[0169]

[0170] Among them, θ(n) is the cutoff criterion; B Loughlin is the Loughlin instantaneous bandwidth; is the weighted mean instantaneous frequency.

[0171] Step 4: Use the sliding window method to perform overlapping sampling on the selected IMF and divide it into multiple sample sets. Calculate the time domain features, frequency domain features, energy entropy features, sample entropy features and singular value features of the vibration signal to construct a feature vector.

[0172] Step 4 specifically includes the following steps:

[0173] Step 4.1: Calculating the time domain characteristics of the vibration signal includes calculating the peak value Y1, average value Y2, and root mean square value Y3 of the vibration signal. The calculation formulas for the peak value Y1, average value Y2, and root mean square value Y3 are:

[0174] Y1=max(x(n))

[0175]

[0176]

[0177] Wherein, x(n) is the original signal, n=1,2,3,…N, and N is the data length.

[0178] Step 4.2: Calculating the frequency domain characteristics of the vibration signal includes calculating the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 of the vibration signal. The calculation formulas for the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 are:

[0179]

[0180]

[0181]

[0182] Where, the spectrum sequence of s(k) is k=1,2,3,…K, where K is the spectrum line; f k is the frequency value of the kth spectrum line.

[0183] Step 4.3: Energy entropy feature EF i The calculation formula is:

[0184]

[0185] Among them, ES i is the energy value of the IMF component, and E represents the sum of the IMF components.

[0186] The calculation formula of sample entropy feature SampEn is:

[0187]

[0188] Where m is the pattern dimension; r is the similarity tolerance threshold; B m(r) is the probability that two sequences match m points under r, A m (r) is the probability that the two sequences match m+1 points.

[0189] Singular Value Feature SVF i The calculation formula is

[0190]

[0191] Among them, A represents the matrix constructed by IMF components, U m×m and V n×n are the left and right singular matrices of A, Λ is an m×n diagonal matrix, λ1,λ2,…,λ i are the singular values ​​of the matrix A.

[0192] Step 5: Use random forest to screen the feature vectors, obtain the optimal feature vector, and divide the sample set into training set and test set in a ratio of 7:3. The classification accuracy CA of random forest is:

[0193]

[0194] Among them, CM(i,j) is the number of times the i-th sample in the mixing matrix is ​​classified into type j, and when i=j, it indicates correct classification.

[0195] Step 6: Construct and train the MSCNN-LTransformer fault diagnosis model by fusing the multi-scale convolutional neural network MSCNN and LinearTransformer (LTransformer). MSCNN extracts features from time series of different spans, and LTransformer captures the global features of the time series.

[0196] Step 6 specifically includes the following steps:

[0197] Step 6.1: Input the training set data into the MSCNN layer. The MSCNN layer is composed of convolution kernels of different sizes. Convolution kernels of different sizes can obtain different feature sizes, so MSCNN can obtain richer feature information. This article uses convolution kernel sizes of 1, 3, 7, and 9. The operation of MSCNN can be expressed as:

[0198] S1=C1(x)

[0199] S2=C3(S1)

[0200] S3=C7(S2)

[0201] S4=C9(S3)

[0202] S=Concat(S2,S3,S4)

[0203] Among them, C n represents a one-dimensional convolutional element with a convolution kernel size of n, where n = 1, 3, 7, 9; x represents the input; S represents the output state of the MCNN layer; Concat(·) represents the concatenation operation.

[0204] Step 6.2: Encode the local feature information extracted from step 6.1 and input it into LTransformer. LTransformer is an improved model of Transformer, which mainly consists of input layer, feature encoder layer, feature decoder layer and output layer. The feature encoder and feature decoder are mainly composed of multi-head attention module MSA and feedforward neural network module FFN. Each module uses residual connection Add and then performs layer normalization Norm.

[0205] LTransformer transforms the traditional softmax attention into dot product attention based on feature mapping, which reduces the computational cost and complexity of the algorithm. i 'The calculation formula is as follows:

[0206]

[0207] Among them, sim(·) is the similarity function; Q is the query vector; K is the key vector; V is the value vector.

[0208] MSA can learn features from different feature subspaces and has good global information extraction capabilities. The calculation formula of MSA is as follows:

[0209] M(Q,K,V)=Concat(H1,H2,…H h )·W

[0210] H i =V i '(Q i ,K i ,V i ),i=1,2…h

[0211] Among them, M(·) is the multi-head attention mechanism operation, H i is the i-th head of the attention mechanism, W is the multi-head attention weight matrix, V i 'For Q i The cumulative value of .

[0212] FFN consists of two fully connected layers. The calculation formula of FFN is as follows:

[0213] FFN(x)=max(0,xW1+b1)W2+b2

[0214] Among them, x is the normalized output matrix of the attention layer; W1, W2 are weight vectors; and b is the bias term.

[0215] A residual connection Add is used between each submodule of the feature encoder and feature decoder to prevent network degradation, and then the layer normalization Norm is used to normalize the activation value of each layer to avoid overfitting and speed up the training. The calculation formula is:

[0216] X B =LN(M(X A )+X A )

[0217] X C =LN(FFN(X B )+X B )

[0218] Among them, X A is the input matrix of the encoding layer; X B Input matrix for FFN layer; X C is the output matrix of the coding layer; LN is the layer normalization process.

[0219] Step 6.3: Input the feature encoder output result into the linear layer and the Softmax classifier layer and then output the fault diagnosis result.

[0220] Step 7: Use the test set data to test the trained MSCNN-LTransformer fault diagnosis model.

[0221] In order to compare the fault diagnosis accuracy of a petrochemical unit bearing fault diagnosis method (ISCHO-TVFEMD-MSCNN-LTransformer) proposed in this invention, four fault diagnosis models, TVFEMD-MSCNN, TVFEMD-MSCNN-LTransformer, SCHO-TVFEMD-MSCNN-LTransformer and ISCHO-TVFEMD-MSCNN-LTransformer, are compared at the same time. The fault diagnosis accuracy of different fault diagnosis models is shown in Table 1 and Figure 5 As shown in Table 1 and Figure 5 It can be clearly seen that the fault diagnosis accuracy of ISCHO-TVFEMD-MSCNN-LTransformer is significantly better than other models.

[0222] Table 1 Comparison of fault diagnosis accuracy results

[0223] Fault diagnosis model accuracy TVFEMD-MSCNN 80.33％ TVFEMD-MSCNN-LTransformer 87.33％ SCHO-TVFEMD-MSCNN-LTransformer 94.67％ ISCHO-TVFEMD-MSCNN-LTransformer 99.33％

[0224] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments. Those skilled in the art can make various deformations or modifications within the scope of the claims, which will not affect the essential content of the present invention.

Claims

1. A method for diagnosing bearing faults of a petrochemical unit, characterized in that: The steps include: Step 1: Collect vibration signals of four different rolling bearing conditions: normal condition, inner ring fault, outer ring fault, and rolling element fault at the same speed; Step 2: Use random inertia weight and Levy flight strategy to improve the hyperbolic sine-cosine optimization algorithm SCHO to obtain the improved hyperbolic sine-cosine optimization algorithm ISCHO; Step 2.1: Since the second stage of the hyperbolic sine-cosine optimization algorithm SCHO exploration is to perform an undirected search at the current position, a random inertia weight strategy is introduced to avoid the local optimal situation. The improved formula is as follows: w=W max +r 14 (W max -W min )+σ·R Where t is the current iteration number; is the jth position of the ith solution in the current iteration, i = 1, 2, ... N, N is the number of candidate solutions, j = 1, 2, ... dim, dim is the dimension of the problem; ε represents a very small positive number; is the jth position of the current optimal solution; W2 is The weight coefficient; r is a random number between the interval [0,1]; w is the random inertia weight; W max is the maximum value of random inertia weight; W min is the minimum value of random inertia weight; σ is the deviation between random inertia weight and mathematical expectation; R is a random number of normal distribution; Step 2.2: Since the second stage of SCHO development is to conduct in-depth development around the optimal solution, the Levy flight strategy is introduced to speed up the algorithm. The improved formula is as follows: Where l is the step size scaling factor; levy(λ) is the random path of Levy flight; Step 3: Use ISCHO to optimize the parameters of the time-varying filter empirical mode decomposition TVFEMD, which include the bandwidth threshold ξ and the B-spline order. Use the optimized TVFEMD to decompose the collected bearing vibration signal to obtain multiple intrinsic mode function components IMF, and further select the IMF containing more fault characteristics; Step 4: Use the sliding window method to perform overlapping sampling on the selected IMFs, divide them into multiple sample sets, calculate the time domain features, frequency domain features, energy entropy features, sample entropy features and singular value features of the vibration signal, and construct the feature vector; Step 5: Use random forest to screen the feature vectors, obtain the optimal feature vector, and divide the sample set into training set and test set in a ratio of 7:3; Step 6: Construct and train the MSCNN-LTransformer fault diagnosis model by fusing the multi-scale convolutional neural network MSCNN and Linear Transformer (LTransformer). MSCNN extracts features from time series of different spans, and LTransformer captures the global features of the time series. Step 7: Use the test set data to test the trained MSCNN-LTransformer fault diagnosis model.

2. A petrochemical unit bearing fault diagnosis method according to claim 1, characterized in that: The step 3 of optimizing TVFEMD using ISCHO includes the following steps: Step 3.1: Initialize the relevant parameters of ISCHO, including population size, problem dimension dim, number of iterations t, and the switching value T between the first and second stages of exploration and development, where T is calculated as follows: Among them, floor is the rounding down function; Max_iteration represents the maximum value of iteration; ct is the coefficient of the two-stage switching; Step 3.2: Randomly initialize a set of candidate solutions. The candidate solution calculation formula is as follows: X=rand(N,dim)×(ub-lb)+lb Among them, X is a set of random candidate solutions; rand is a random number in [0,1]; ub is the upper bound of the variable; lb is the lower bound of the variable; Step 3.3: Calculate the objective function value of the candidate solution and determine the current optimal solution. The objective function used is the minimum envelope entropy. The calculation formula of the minimum envelope entropy Ep is: Where a(h) is the envelope signal, h = 1, 2, ... N; p h is the normalized form of a(h); Step 3.4: Calculate the value of A, which is the value converted by the Sinh and Cosh conversion mechanisms in the exploration and development stages. When A is greater than 1, it enters the exploration stage, and when A is less than 1, it enters the development stage. The calculation formula of A is as follows: Among them, p is the balance coefficient of the exploration phase in the iterative process; q is the balance coefficient of the development phase in the iterative process; r is a random number between the interval [0,1]; If t = BS k When , find the current position of the second solution and update BS k and search space, reallocate candidate solutions, BS k The calculation formula is as follows: Among them, BS k is the number of iterations of the current bounded search strategy; α is the sensitivity coefficient value; The calculation formula of the search space is as follows: Among them, ub k is the upper bound of the search space; lb k is the lower bound of the search space; is the jth position of the current optimal solution; represents the jth position of the suboptimal solution; Step 3.5: When A>1, enter the exploration stage and calculate W1 and W2. The calculation formulas of W1 and W2 are as follows: W1=r3×a1×(coshr4+u×sinhr4-1) W2=r6×a2 Among them, W1 is The weight coefficient of The weight coefficient of a1 and a2 are monotonically decreasing functions; u is the sensitivity coefficient for controlling the exploration accuracy in the first stage; When t≤T, the first stage of exploration begins. The position update formula of the candidate solution in the first stage of exploration is as follows: When t ≥ T, the second stage of exploration begins. The position update formula of the candidate solution in the second stage of exploration is as follows: w=W max +r 14 (W max -W min )+σ·R When A<1, it enters the development stage and calculates the weight coefficient W3 of the first stage of development. The calculation formula is as follows: <h2 style=";text-align:left;direction:ltr">W3 = r9×a1×(coshr<h2 style=";text-align:left;direction:ltr"> 10 <h2 style=";text-align:left;direction:ltr"> +u×sinhr<h2 style=";text-align:left;direction:ltr"> 10 <h2 style=";text-align:left;direction:ltr"> ) When t≤T, the first stage of development begins. The position update formula of the candidate solution in the first stage of development is as follows: When t≥T, the second stage of development begins. The position update formula of the candidate solution in the second stage of development is as follows: Among them, W2 controls the extent of the second stage of development; Step 3.6: Determine whether the termination condition is met. If it is met, output the optimal solution. If not, repeat steps 3.3 to 3.

5.

3. A petrochemical unit bearing fault diagnosis method according to claim 1, characterized in that: The optimized TVFEMD is used to decompose the collected bearing vibration signal, and the calculation formula is: Among them, θ(n) is the cutoff criterion; B Loughlin is the Loughlin instantaneous bandwidth; is the weighted mean instantaneous frequency.

4. A petrochemical unit bearing fault diagnosis method according to claim 1, characterized in that: The step 4 of calculating the time domain features, frequency domain features, energy entropy features, sample entropy features and singular value features of the vibration signal includes the following steps: Step 4.1: Calculating the time domain characteristics of the vibration signal includes calculating the peak value Y1, average value Y2, and root mean square value Y3 of the vibration signal. The calculation formulas for the peak value Y1, average value Y2, and root mean square value Y3 are: Y1=max(x(n)) Where x(n) is the original signal, n=1,2,3,…N, N is the data length; Step 4.2: Calculating the frequency domain characteristics of the vibration signal includes calculating the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 of the vibration signal. The calculation formulas for the frequency domain amplitude mean value Y4, the frequency centroid value Y5, and the frequency amplitude root mean square value Y6 are: Where, the spectrum sequence of s(k) is k=1,2,3,…K, where K is the spectrum line; f k is the frequency value of the kth spectrum line; Step 4.3: Energy entropy feature EF i The calculation formula is: Among them, ES i is the energy value of the IMF component, E represents the sum of the IMF components; The calculation formula of sample entropy feature SampEn is: Where m is the pattern dimension; r is the similarity tolerance threshold; B m (r) is the probability that two sequences match m points under r, A m (r) is the probability that two sequences match m+1 points; Singular Value Feature SVF i The calculation formula is Among them, A represents the matrix constructed by IMF components, U m×m and V n×n are the left and right singular matrices of A, Λ is an m×n diagonal matrix, λ1,λ2,…,λ i are the singular values ​​of the matrix A.

5. A petrochemical unit bearing fault diagnosis method according to claim 1, characterized in that: The MSCNN-LTransformer fault diagnosis model in step 6 includes the following steps: Step 6.1: Input the training set data into the MSCNN layer. The MSCNN layer is composed of convolution kernels of different sizes. Convolution kernels of different sizes can obtain different feature sizes, so MSCNN can obtain richer feature information. Using convolution kernel sizes of 1, 3, 7, and 9, the operation of MSCNN can be expressed as: S1=C1(x) S2=C3(S1) S3=C7(S2) S4=C9(S3) S=Concat(S2,S3,S4) Among them, C n represents a one-dimensional convolutional element with a convolution kernel size of n, where n = 1, 3, 7, 9; x represents the input; S represents the output state of the MCNN layer; Concat(·) represents the concatenation operation; Step 6.2: Encode the local feature information extracted from step 6.1 and input it into LTransformer. LTransformer consists of input layer, feature encoder layer, feature decoder layer and output layer. The feature encoder and feature decoder consist of multi-head attention module MSA and feed-forward neural network module FFN. Each module uses residual connection Add and then performs layer normalization Norm. LTransformer transforms the traditional softmax attention into dot product attention based on feature map. i 'The calculation formula is as follows: Where sim(·) is the similarity function; Q is the query vector; K is the key vector; V is the value vector; The multi-head attention module MSA learns features from different feature subspaces. The calculation formula of MSA is as follows: M(Q,K,V)=Concat(H1,H2,…H h )·W H i =V i '(Q i ,K i ,V i ),i=1,2…h Among them, M(·) is the multi-head attention mechanism operation, H i is the i-th head of the attention mechanism, W is the multi-head attention weight matrix, V i 'For Q i The cumulative value of FFN consists of two fully connected layers. The calculation formula of FFN is as follows: FFN(x)=max(0,xW1+b1)W2+b2 Among them, x is the normalized output matrix of the attention layer; W1, W2 are weight vectors; b is the bias term; A residual connection Add is used between each submodule of the feature encoder and feature decoder to prevent network degradation, and then the layer normalization Norm is used to normalize the activation value of each layer to avoid overfitting and speed up the training. The calculation formula is: X B =LN(M(X A )+X A ) X C =LN(FFN(X B )+X B ) Among them, X A is the input matrix of the encoding layer; X B Input matrix for FFN layer; X C is the output matrix of the coding layer; LN is the layer normalization process; Step 6.3: Input the feature encoder output result into the linear layer and the Softmax classifier layer and then output the fault diagnosis result.

Images