Rolling bearing fault diagnosis method and system based on small sample data

Abstract

The invention relates to a rolling bearing fault diagnosis method and system based on small sample data, and the method comprises the steps: taking a variable working condition vibration signal as an input sample, carrying out the optimization of the number of modal components and penalty factors of the variational modal decomposition of the input sample through a sparrow search algorithm, and obtaining an optimal parameter combination; performing variational mode decomposition on the input sample based on the optimal parameter combination to obtain a plurality of intrinsic mode function components, and obtaining an envelope order spectrum of each intrinsic mode function component through calculation order tracking and envelope demodulation; according to the sparseness of the envelope spectrum, screening out the intrinsic mode function component with a better sparse effect to carry out signal reconstruction; and converting the reconstructed signal into a two-dimensional image, and performing model training on the multi-scale subtraction network to obtain a diagnosis model. According to the method, only limited data samples under part of working conditions are needed, and the bearing state can be accurately identified under small samples.

Description

Technical Field

[0001] The present invention relates to the technical field of rolling bearing fault diagnosis, and in particular to a rolling bearing fault diagnosis method and system based on small sample data. Background Art

[0002] Wind energy is intermittent and random, which causes the unit to operate under variable speed and load conditions for a long time. Changes in operating conditions will affect the vibration characteristics of the rolling bearings of wind-driven motors, increasing the difficulty of rolling bearing fault diagnosis in wind turbines. The fault characteristics of rolling bearings will change with changes in operating conditions. How to extract fault characteristics that are less affected by changes in operating conditions becomes the key to rolling bearing fault diagnosis under variable conditions. Although existing diagnostic methods have achieved good diagnostic results, their models require a large number of fault samples under different operating conditions during the training process. In reality, it is difficult to obtain enough variable operating condition bearing fault samples. Therefore, variable operating condition fault diagnosis of rolling bearings in wind turbines often encounters the problem of insufficient sample size. Summary of the invention

[0003] In view of the deficiencies in the prior art, the present invention provides a rolling bearing fault diagnosis method and system based on small sample data, the purpose of which is to require only limited data samples under certain working conditions, that is, to accurately identify the bearing status under small samples.

[0004] The technical solution adopted by the present invention is as follows:

[0005] The present invention provides a rolling bearing fault diagnosis method based on small sample data, comprising:

[0006] Taking the variable working condition vibration signal as the input sample, the number of modal components k and the penalty factor α of the variational modal decomposition of the input sample are optimized by the sparrow search algorithm to obtain the optimal parameter combination [k, α];

[0007] Performing variational mode decomposition on the input sample based on the optimal parameter combination to obtain k intrinsic mode function components, and obtaining the envelope order spectrum of each intrinsic mode function component by calculating order tracking and envelope demodulation;

[0008] According to the sparsity of the envelope spectrum, the intrinsic mode function component with better sparsity effect is selected to perform signal reconstruction;

[0009] Converting the reconstructed signal into a two-dimensional image and dividing it into a training set and a test set; inputting the training set into a multi-scale subtraction network for model training to obtain a diagnostic model;

[0010] The test set is input into the diagnostic model to perform bearing fault diagnosis.

[0011] Further technical solutions are:

[0012] The optimizing the number k of modal components and the penalty factor α of the variational modal decomposition of the input sample by using the sparrow search algorithm includes:

[0013] Set a sparrow population X, initialize a task sequence for each sparrow, and the task sequence includes the variable values ​​of the number of modal components k and the penalty factor α;

[0014]

[0015] d is the number of variable values ​​to be optimized, and n is the number of sparrows in the initial population;

[0016] Using envelope entropy as the fitness function value, calculate the fitness F of the initial population X and find the best and worst individuals;

[0017]

[0018] f() is the fitness function;

[0019] Generate alarm values;

[0020] Calculate and discover new positions of explorers, followers, and sentinels;

[0021] If the new position is better than the current position, the position is updated until the iteration stop requirement is met. The global optimal solution in the updated sparrow population is determined as the optimal solution of the variational problem mathematical model of variational mode decomposition, and the optimal parameter combination and fitness value are output.

[0022] The variational mode decomposition of the input sample comprises:

[0023] Constructing mathematical model of variational problem:

[0024]

[0025] In the formula, u k is the kth modal component obtained by decomposition; ω k for u k The center frequency of; δ(t) is the pulse function; k is the number of modal components; f is the spectrum of the original vibration signal;

[0026] The alternating direction multiplier method is used to continuously update each modal component, center frequency and Lagrange multiplier and find the optimal solution. The updated expression is as follows:

[0027]

[0028] In the formula, α is the penalty factor, λ is the Lagrange multiplier, i represents the mode number; n represents the number of updates; ω represents the frequency coordinate in the spectrum; represents the estimated value of the Lagrange multiplier; represents the kth modal component u of the n+1th iteration k An estimated value of Represents the center frequency of the n+1th iteration The estimated value of .

[0029] The explorer's location update includes:

[0030]

[0031] Where t is the current iteration number, j = 1, 2…, d; iter max is the maximum number of iterations, which is a constant; X i,j is the position information of the i-th sparrow in the j-th dimension; random number ξ∈[0,1]; warning value R 2 ∈[0,1]; safety value ST∈[0.5,1]; q is a random number that follows a normal distribution; L is a 1×d matrix whose elements are all 1;

[0032] The follower's location update includes:

[0033]

[0034] Where, X p The best position for the current explorer; X worst is the worst position in the current population; A is a 1×d matrix column whose elements are randomly 1 or -1, and A + =A T (AA T ) -1 .

[0035] The calculation order tracking comprises:

[0036] In the first step, the input sample vibration signal and the corresponding speed signal are sampled at a constant time interval Δt, and are recorded as x(n) and s(n) respectively; wherein s(n) is a sinusoidal voltage signal whose frequency changes with the speed, and each zero crossing point thereof corresponds to the time value of the rotating shaft rotating one circle; n represents a sampling point in the signal;

[0037] The second step is to search for a point in s(n) that satisfies s(n)≤0 and s(n+1)≥0, and calculate the time value t corresponding to each rotation of the rotation axis by the following formula: zi :

[0038]

[0039] Among them, s(n) and s(n+1) correspond to time t s(n) and t s(n+1)The waveform between them is linear;

[0040] Calculate the time series t corresponding to the equal angle increment Δθ through s(n) ik , i=1,2,3…, for t ik Calculate the difference to obtain the time T required for the rotating axis to rotate once i ,include:

[0041] The initial velocity v of the first two rotations is calculated by the following formula: 10 and angular acceleration a 1 :

[0042]

[0043] The time value t corresponding to each sampling point in the first two rotations is calculated by the following formula: ik , instantaneous speed v ik and the initial velocity v of the third cycle 30 :

[0044]

[0045] In the above formula, i is the number of rotations, i=1,2; k is the kth sampling point in each rotation; M is the number of equal-angle sampling points for each rotation of the rotating shaft, and 2kπ / M is the angle of rotation of the rotating shaft corresponding to the kth sampling point;

[0046] Further deducing the angular acceleration a for i ≥ 3 cycles i and initial velocity v i+1,0 , and the corresponding t ik and v ik ;

[0047] The third step is to calculate the ik , i = 1, 2, 3…, interpolate x(n) and calculate the value of x(n) at t ik The corresponding amplitude is obtained to obtain the equal-angle resampled signal x(t ik ).

[0048] The multi-scale subtraction network includes a feature extraction module, a multi-scale subtraction feature extraction module, a lightweight module, an attention mechanism, a global average pooling layer and a fully connected layer;

[0049] The feature extraction module uses Res2Net-50 to extract features at four levels and reduces the feature channels to 64 using 3×3 convolution;

[0050] The multi-scale subtraction feature extraction module takes the features of the four levels as input, obtains four complementary enhancement features, and outputs the four complementary enhancement features after fusion;

[0051] The lightweight module uses the fused features output by the multi-scale subtraction feature extraction module and performs lightweight extraction;

[0052] The channel attention mechanism gives more weight to important features in the lightweight extracted features;

[0053] The global average pooling layer and the fully connected layer output diagnostic classification results.

[0054] The present invention also provides a rolling bearing fault diagnosis system based on small sample data, comprising:

[0055] An optimization module, which uses a variable working condition vibration signal as an input sample, optimizes the number of modal components k and the penalty factor α of the variational modal decomposition of the input sample by a sparrow search algorithm, and obtains an optimal parameter combination [k, α];

[0056] A variational mode decomposition module, which performs variational mode decomposition on the input sample based on an optimal parameter combination to obtain k intrinsic mode function components, and obtains the envelope order spectrum of each intrinsic mode function component by calculating order tracking and envelope demodulation;

[0057] A signal reconstruction module, which selects the intrinsic mode function component with better sparse effect according to the envelope spectrum sparsity to perform signal reconstruction;

[0058] An image diagnosis model building module converts the reconstructed signal into a two-dimensional image and divides it into a training set and a test set; the training set is input into a multi-scale subtraction network for model training to obtain a diagnosis model;

[0059] A diagnosis module inputs the test set into the diagnosis model to perform bearing fault diagnosis.

[0060] Further technical solutions are:

[0061] The optimizing the number k of modal components and the penalty factor α of the variational modal decomposition of the input sample by using the sparrow search algorithm includes:

[0062] Set a sparrow population X, initialize a task sequence for each sparrow, and the task sequence includes the variable values ​​of the number of modal components k and the penalty factor α;

[0063]

[0064] d is the number of variable values ​​to be optimized, and n is the number of sparrows in the initial population;

[0065] Using envelope entropy as the fitness function value, calculate the fitness F of the initial population X and find the best and worst individuals;

[0066]

[0067] f() is the fitness function;

[0068] Generate alarm values;

[0069] Calculate and discover new positions of explorers, followers, and sentinels;

[0070] If the new position is better than the current position, the position is updated until the iteration stop requirement is met, and the global optimal solution in the updated sparrow population is determined as the optimal solution of the variational problem mathematical model of variational mode decomposition, and the optimal parameter combination and fitness value are output;

[0071] The variational mode decomposition of the input sample comprises:

[0072] Constructing mathematical model of variational problem:

[0073]

[0074] In the formula, u k is the kth modal component obtained by decomposition; ω k for u k The center frequency of; δ(t) is the pulse function; k is the number of modal components; f is the spectrum of the original vibration signal;

[0075] The alternating direction multiplier method is used to continuously update each modal component, center frequency and Lagrange multiplier and find the optimal solution. The updated expression is as follows:

[0076]

[0077] In the formula, α is the penalty factor, λ is the Lagrange multiplier, i represents the mode number; n represents the number of updates; ω represents the frequency coordinate in the spectrum; represents the estimated value of the Lagrange multiplier λ; represents the kth modal component u of the n+1th iteration k An estimated value of Represents the center frequency of the n+1th iteration The estimated value of .

[0078] The calculation order tracking comprises:

[0079] In the first step, the input sample vibration signal and the corresponding speed signal are sampled at a constant time interval Δt, and are recorded as x(n) and s(n) respectively; wherein s(n) is a sinusoidal voltage signal whose frequency changes with the speed, and each zero crossing point thereof corresponds to the time value of the rotating shaft rotating one circle; n represents a sampling point in the signal;

[0080] The second step is to search for a point in s(n) that satisfies s(n)≤0 and s(n+1)≥0, and calculate the time value t corresponding to each rotation of the rotation axis by the following formula: zi :

[0081]

[0082] Among them, s(n) and s(n+1) correspond to time t s(n) and t s(n+1) The waveform between them is linear;

[0083] Calculate the time series t corresponding to the equal angle increment Δθ through s(n) ik , i=1,2,3…, for t ik Calculate the difference to obtain the time T required for the rotating axis to rotate once i ,include:

[0084] The initial velocity v of the first two rotations is calculated by the following formula: 10 and angular acceleration a 1 :

[0085]

[0086] The time value t corresponding to each sampling point in the first two rotations is calculated by the following formula: ik , instantaneous speed v ik and the initial velocity v of the third cycle 30 :

[0087]

[0088] In the above formula, i is the number of rotations, i=1,2; k is the kth sampling point in each rotation; M is the number of equal-angle sampling points for each rotation of the rotating shaft, and 2kπ / M is the angle of rotation of the rotating shaft corresponding to the kth sampling point;

[0089] Further deducing the angular acceleration a for i ≥ 3 cycles i and initial velocity v i+1,0 , and the corresponding t ik and v ik ;

[0090] The third step is to calculate the ik , i=1,2,3…, interpolate x(n) and calculate the value of x(n) at t ik The corresponding amplitude is obtained to obtain the equal-angle resampled signal x(t ik ).

[0091] The multi-scale subtraction network includes a feature extraction module, a multi-scale subtraction feature extraction module, a lightweight module, an attention mechanism, a global average pooling layer and a fully connected layer;

[0092] The feature extraction module uses Res2Net-50 to extract features at four levels and reduces the feature channels to 64 using 3×3 convolution;

[0093] The multi-scale subtraction feature extraction module takes the features of the four levels as input, obtains four complementary enhancement features, and outputs the four complementary enhancement features after fusion;

[0094] The lightweight module uses the fused features output by the multi-scale subtraction feature extraction module and performs lightweight extraction;

[0095] The channel attention mechanism gives more weight to important features in the lightweight extracted features;

[0096] The global average pooling layer and the fully connected layer output diagnostic classification results.

[0097] The beneficial effects of the present invention are as follows:

[0098] The present invention combines signal processing with deep learning and proposes a rolling bearing fault diagnosis method based on optimized variational mode decomposition and improved multi-scale subtraction network. It has higher diagnostic accuracy than commonly used models and has good recognition performance under small sample conditions and variable working conditions. The specific advantages include the following:

[0099] The deep learning neural network used in the present invention can achieve high-precision diagnosis of rolling bearing fault types under variable conditions by only training samples under some working conditions. The diagnostic model of the present invention is less dependent on the number of training samples. According to experiments, a higher diagnostic accuracy can be obtained when the number of training samples is 25, which has great diagnostic superiority in the case of small samples.

[0100] The present invention optimizes variational modal decomposition parameters through a sparrow search algorithm, calculates order tracking and envelope demodulation to obtain the envelope order spectrum of the intrinsic mode function components, realizes the stable angular domain conversion of non-stationary time domain signals, and reduces the influence of operating condition changes on fault feature extraction.

[0101] Other features and advantages of the present invention will be set forth in the following description or may be learned by practicing the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] Figure 1 Schematic diagram of the process of the embodiment of the present invention.

[0103] Figure 2 It is a structural schematic diagram of a multi-scale subtraction network diagnosis model (SVMD-COT-MSSNN) according to an embodiment of the present invention.

[0104] Figure 3It is the diagnostic accuracy of different models under variable load conditions of the embodiment of the present invention.

[0105] Figure 4 It is the diagnostic accuracy of the multi-scale subtraction network diagnostic model under small sample conditions in the embodiment of the present invention. DETAILED DESCRIPTION

[0106] The specific implementation of the present invention is described below with reference to the accompanying drawings.

[0107] See also Figure 1 , a rolling bearing fault diagnosis method based on small sample data of this embodiment includes:

[0108] Taking the variable working condition vibration signal as the input sample, the number of modal components k and the penalty factor α of the variational mode decomposition (VMD) of the input sample are optimized by the Sparrow Search Algorithm (SSA) to obtain the optimal parameter combination [k, α];

[0109] Based on the optimal parameter combination, the input sample is subjected to variational mode decomposition to obtain k intrinsic mode function (IMF) components, and the envelope order spectrum of each intrinsic mode function component is obtained by computed order tracking (COT) and envelope demodulation;

[0110] According to the sparsity of the envelope spectrum, the intrinsic mode function component with better sparsity effect is selected to perform signal reconstruction;

[0111] The reconstructed signal is converted into a two-dimensional image, and divided into a training set and a test set; the training set is input into a multi-scale subtraction neural network (MSSNN) for model training to obtain a diagnostic model;

[0112] The test set is input into the diagnostic model to perform bearing fault diagnosis.

[0113] The optimizing the number k of modal components and the penalty factor α of the variational modal decomposition of the input sample by using the sparrow search algorithm includes:

[0114] Set a sparrow population X, initialize a task sequence for each sparrow, and the task sequence includes the variable values ​​of the number of modal components k and the penalty factor α;

[0115]

[0116] d is the number of variable values ​​to be optimized, and n is the number of sparrows in the initial population;

[0117] Using envelope entropy as the fitness function value, calculate the fitness F of the initial population X and find the best and worst individuals;

[0118]

[0119] f() is the fitness function;

[0120] Generate alarm values;

[0121] Calculate and discover new positions of explorers, followers, and sentinels;

[0122] If the new position is better than the current position, the position is updated until the iteration stop requirement is met. The global optimal solution in the updated sparrow population is determined as the optimal solution of the variational problem mathematical model of variational mode decomposition, and the optimal parameter combination and fitness value are output.

[0123] Specifically, in this embodiment, the number of sparrow populations is set to 30, the number of iterations is set to 20, the envelope entropy is used as the fitness function value, the optimization range of k is [3,9], the optimization range of α is [10,5000]; the iteration accuracy of VMD is set to 10 -7 .

[0124] Specifically, variational mode decomposition is a signal processing technology that is further developed on the basis of classical Wiener filtering and Hilbert-Huang transform. This method believes that real signals are composed of multiple sine waves or cosine waves, and the amplitude, frequency and number of decomposition waves of different signals are different. Variational mode decomposition constructs a variational problem to solve each modal component. It must be ensured that each decomposed eigenmode component has a finite bandwidth with a center frequency, and the sum of the estimated bandwidths of each modal component is always required to be kept to a minimum. Therefore, its powerful adaptive ability can determine the number of modal decompositions, finite bandwidth and optimal center frequency according to actual needs. The mathematical model of its variational problem is as follows:

[0125]

[0126] In the formula, u k is the kth modal component obtained by decomposition; ω k for u k The center frequency of; δ(t) is the pulse function; k is the number of modal components; f is the spectrum of the original vibration signal;

[0127] The alternating direction multiplier method is used to continuously update each modal component Center frequency and the Lagrange multiplier λ n+1And find the optimal solution, the update expression is as follows:

[0128]

[0129] In the formula, α is the penalty factor, λ is the Lagrange multiplier, i represents the mode number; n represents the number of updates; ω represents the frequency coordinate in the spectrum; represents the estimated value of the Lagrange multiplier λ; represents the kth modal component u of the n+1th iteration k An estimated value of Represents the center frequency of the n+1th iteration The estimated value of .

[0130] Specifically, the fitness value of the explorer is relatively high, and they have priority in discovering and obtaining food. At the same time, their search range is much larger than that of the follower. The position update of the explorer includes:

[0131]

[0132] Where t is the current iteration number, j = 1, 2…, d; iter max is the maximum number of iterations, which is a constant; X i,j is the position information of the i-th sparrow in the j-th dimension; random number ξ∈[0,1]; warning value R 2 ∈[0,1]; safety value ST∈[0.5,1]; q is a random number that follows a normal distribution; L is a 1×d matrix whose elements are all 1;

[0133] Followers will closely follow the explorer's steps. Once they find that the explorer has found better food, they will immediately grab the food. If successful, they can get the food immediately. Otherwise, they can only continue to follow. The position update of the follower includes:

[0134]

[0135] Where, X p The best position for the current explorer; X worst is the worst position in the current population; A is a 1×d matrix column whose elements are randomly 1 or -1, and A + =A T (AA T ) -1 .

[0136] The calculation order tracking comprises:

[0137] In the first step, the input sample vibration signal and the corresponding speed signal are sampled at a constant time interval Δt, and are recorded as x(n) and s(n) respectively; wherein s(n) is a sinusoidal voltage signal whose frequency changes with the speed, and each zero crossing point thereof corresponds to the time value of the rotating shaft rotating one circle; n represents a sampling point in the signal;

[0138] The second step is to search for a point in s(n) that satisfies s(n)≤0 and s(n+1)≥0, and calculate the time value t corresponding to each rotation of the rotation axis by the following formula: zi :

[0139]

[0140] Among them, s(n) and s(n+1) correspond to time t s(n) and t s(n+1) The waveform between them is linear;

[0141] Calculate the time series t corresponding to the equal angle increment Δθ through s(n) ik , i=1,2,3…, for t ik Calculate the difference to obtain the time T required for the rotating axis to rotate once i ,include:

[0142] The initial velocity v of the first two rotations is calculated by the following formula: 10 and angular acceleration a 1 :

[0143]

[0144] Specifically, the rotating axis performs uniformly accelerated motion during one rotation, and maintains the same angular acceleration during the first two rotations;

[0145] The time value t corresponding to each sampling point in the first two rotations is calculated by the following formula: ik , instantaneous speed v ik and the initial velocity v of the third cycle 30 :

[0146]

[0147] In the above formula, i is the number of rotations, i=1,2; k is the kth sampling point in each rotation; M is the number of equal-angle sampling points for each rotation of the rotating shaft, and 2kπ / M is the angle of rotation of the rotating shaft corresponding to the kth sampling point;

[0148] Further deducing the angular acceleration a for i ≥ 3 cycles i and initial velocity v i+1,0 , and the corresponding t ik and v ik ;

[0149] The third step is to calculate the ik , i=1,2,3…, interpolate x(n) and calculate the value of x(n) at t ik The corresponding amplitude is obtained to obtain the equal-angle resampled signal x(t ik ). At this time, the equal-angle resampled signal x(t ik ) is already a stable angular domain signal and can be processed by traditional frequency domain analysis and other means.

[0150] The envelope spectrum sparsity can measure the strength of the impulse signal in the frequency domain. If there are some peaks in the impulse signal and the area between the peaks is relatively flat, then this impulse signal can be called a sparse signal. When the signal amplitude distribution is relatively uniform and the difference is not large, the sparsity of the signal is weak and the sparsity index is small. When a small number of sharp pulses appear in the signal and the amplitude of individual positions is more prominent, the signal shows a strong sparse characteristic, and the coefficient index also increases accordingly.

[0151] See also Figure 2 , the multi-scale subtraction network includes a feature extraction module, a multi-scale subtraction feature extraction module, a lightweight module, an attention mechanism, a global average pooling layer and a fully connected layer;

[0152] The feature extraction module uses Res2Net-50 to extract features at four levels and reduces the feature channels to 64 using 3×3 convolution;

[0153] The multi-scale subtraction feature extraction module takes the features of the four levels as input, obtains four complementary enhancement features, and outputs the four complementary enhancement features after fusion;

[0154] The lightweight module uses the fused features output by the multi-scale subtraction feature extraction module and performs lightweight extraction;

[0155] The channel attention mechanism gives more weight to important features in the lightweight extracted features;

[0156] The global average pooling layer and the fully connected layer output diagnostic classification results.

[0157] Specifically, the construction and testing of the model of this embodiment are completed in the Tensorflow2.8 deep learning framework. During training, the input of the multi-scale subtraction network is a 32×32 two-dimensional image sample converted from a one-dimensional vibration signal with noise added, and the output is the bearing status 1 to 10. The model training uses the Adan optimizer with a learning rate of 0.001 and a decay rate of 0.005. The cross entropy error is selected as the minimization objective function, the mini-batch is set to 32, and the number of iterations is set to 50.

[0158] In order to further verify the effectiveness of this embodiment, the following experimental analysis is performed by comparing the four conventional models with the model "optimized variational mode decomposition, computational order tracking and improved multi-scale subtraction network model (SVMD-COT-MSSNN)" adopted in this embodiment. The four comparison models are: (1) Multi-branch Convolution Neural Network (MBCNN); (2) Residual Neural Network (ResNet); (3) Improved Multi-scale Subtraction Network Model (MSSNN); Computational Order Tracking and Improved Multi-scale Subtraction Network Model (COT-MSSNN).

[0159] 1.1 Experimental Setup and Dataset

[0160] The two most common variable operating conditions, variable load and variable speed, were selected as the research background. Rolling bearing fault diagnosis experiments under variable load conditions, variable speed conditions and small sample conditions were set up. Two types of data sets were used: Dataset I uses the rolling bearing vibration acceleration data set of the Bearing Data Center of Case Western Reserve University, and Dataset II is the bearing variable speed experimental data of the rolling bearing simulation test bench of the University of Ottawa, Canada.

[0161] (1) Dataset I

[0162] Rolling bearing fault diagnosis under simulated variable load conditions. Four conditions with similar speeds and different loads are selected: 1797rpm / 0, 1772rpm / 1hp, 1750rpm / 2hp, and 1730rpm / 3hp. Among them, 1797rpm represents the speed, 0hp represents no load, the load size corresponding to 1hp is 0.735kw, the load size corresponding to 2hp is 1.471kw, and the load size corresponding to 3hp is 2.206kw. Each condition includes four bearing states: inner ring fault, outer ring fault, rolling element fault, and normal. The sampling frequency is 12kHz, and 1024 data points are selected as a sampling sample. A total of 450 samples are collected under each bearing state. See Table 1 for details.

[0163] Table 1 Dataset I Samples of variable load conditions

[0164]

[0165] 2) Dataset II

[0166] The rolling bearing simulation experiment of the University of Ottawa in Canada uses ER16K rolling bearings, selects four test speed ranges of 900-1080rpm, 1080-1260rpm, 1260-1440rpm, and 1440-1620rpm, and sets the motor load to 0. Experimental bearings with inner ring faults, outer ring faults, and normal conditions are selected for simulation experiments. The bearing vibration data is collected by an acceleration sensor with a sampling frequency of 200kHz. 1024 data points are selected as a sampling sample, and 125 samples are collected under each bearing state. The samples are detailed in Table 2.

[0167] Table 2 Data set II variable speed condition samples

[0168]

[0169] The variable speed working condition diagnosis tasks are divided according to Table 3. Taking the diagnosis task A in data set II as an example, 75 samples of three types of bearing conditions, inner ring fault, outer ring fault and normal, under working condition 1: 900-1080rpm are selected as training samples for this variable speed working condition diagnosis task, and the remaining 300 samples under this working condition are all used as test set samples to test the small sample diagnosis accuracy of the model under variable speed working conditions.

[0170] Table 3. Division of diagnostic tasks for variable speed condition in dataset II

[0171]

[0172]

[0173] 1.2 Fault diagnosis under variable load conditions

[0174] In order to verify the superiority of SVMD-COT-MSSNN diagnostic performance under variable load conditions, it is compared with MBCNN, ResNet, MSSNN, and COT-MSSNN in the variable load condition diagnosis tasks A to D of dataset I. The comparison results are shown in the figure. Figure 3 As shown. Figure 3 It can be seen that when the training sample conditions are completely different from the test sample conditions, the SVMD-COT-MSSNN used in this application can still achieve a test set accuracy of 99.58%, which is much higher than other models, indicating that only training samples under some working conditions are needed to accurately identify the bearing status under other working conditions, verifying the feasibility of this application.

[0175] 1.3 Fault diagnosis of variable speed conditions

[0176] In order to verify the superiority of SVMD-COT-MSSNN diagnostic performance under variable speed conditions, five models, namely MBCNN, ResNet, MSSNN, COT-MSSNN, and SVMD-COT-MSSNN, are compared in variable speed condition diagnosis tasks A to D in dataset II. The comparison results are shown in Table 4.

[0177] Table 4. Diagnostic accuracy of different models under variable speed condition in dataset II

[0178]

[0179] As can be seen from Table 4, the test set accuracy of SVMD-COT-MSSNN and COT-MSSNN is significantly higher than that of MSSNN, ResNet, and MBCNN, indicating that the use of computational order tracking to process non-stationary time domain signals under variable speed can solve the problem of data distribution differences of rolling bearings under variable speed conditions, verifying the effectiveness of computational order tracking. The test set accuracy of SVMD-COT-MSSNN can reach 99.12%, the highest among the five models, indicating that the use of optimized VMD to decompose the signal can improve the performance of rolling bearing fault identification, reflecting the advantages of optimized variational mode decomposition.

[0180] 1.4 Fault diagnosis in small sample situations

[0181] In order to test the diagnostic performance of SVMD-COT-MSSNN under small sample conditions, the number of training samples is changed to observe the changing trend of the diagnostic accuracy of SVMD-COT-MSSNN under variable load and variable speed conditions as the number of training samples changes. Figure 4 It can be seen that SVMD-COT-MSSNN only needs more than 25 training samples to obtain a test accuracy of more than 97%, which verifies the feasibility of this method for rolling bearing fault diagnosis with small samples.

[0182] In summary, SVMD-COT-MSSNN can achieve high-precision diagnosis of rolling bearing fault types under variable working conditions by training samples under only some working conditions: the average diagnostic accuracy under variable load conditions can reach 99.58%, and the average diagnostic accuracy under variable speed conditions can reach 99.12%; the introduction of SVMD-COT realizes the smooth angular domain conversion of non-stationary time domain signals, reducing the impact of working condition changes on fault feature extraction; SVMD-COT-MSSNN has little dependence on the number of training samples, and a higher diagnostic accuracy can be obtained when the number of training samples is 25, which verifies its diagnostic superiority in small sample conditions.

[0183] This embodiment also provides a rolling bearing fault diagnosis system based on small sample data, including:

[0184] An optimization module, which uses a variable working condition vibration signal as an input sample, optimizes the number of modal components k and the penalty factor α of the variational modal decomposition of the input sample by a sparrow search algorithm, and obtains an optimal parameter combination [k, α];

[0185] A variational mode decomposition module, which performs variational mode decomposition on the input sample based on an optimal parameter combination to obtain k intrinsic mode function components, and obtains the envelope order spectrum of each intrinsic mode function component by calculating order tracking and envelope demodulation;

[0186] A signal reconstruction module, which selects the intrinsic mode function component with better sparse effect according to the envelope spectrum sparsity to perform signal reconstruction;

[0187] An image diagnosis model building module converts the reconstructed signal into a two-dimensional image and divides it into a training set and a test set; the training set is input into a multi-scale subtraction network for model training to obtain a diagnosis model;

[0188] A diagnosis module inputs the test set into the diagnosis model to perform bearing fault diagnosis.

[0189] The rolling bearing fault diagnosis method described is used to execute the small sample data.

[0190] Specifically, the diagnostic system includes a memory storing a computer program executable on a processor, and a processor that executes the method when the computer program is executed.

[0191] Those skilled in the art can understand that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention is described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A rolling bearing fault diagnosis method based on small sample data, characterized in that: include: Taking the variable working condition vibration signal as the input sample, the number of modal components k and the penalty factor α of the variational modal decomposition of the input sample are optimized by the sparrow search algorithm to obtain the optimal parameter combination [k, α]; Performing variational mode decomposition on the input sample based on the optimal parameter combination to obtain k intrinsic mode function components, and obtaining the envelope order spectrum of each intrinsic mode function component by calculating order tracking and envelope demodulation; According to the sparsity of the envelope spectrum, the intrinsic mode function component with better sparsity effect is selected to perform signal reconstruction; Converting the reconstructed signal into a two-dimensional image and dividing it into a training set and a test set; inputting the training set into a multi-scale subtraction network for model training to obtain a diagnostic model; The test set is input into the diagnostic model to perform bearing fault diagnosis.

2. The method according to claim 1, characterized in that The optimizing the number k of modal components and the penalty factor α of the variational modal decomposition of the input sample by using the sparrow search algorithm includes: Set a sparrow population X, initialize a task sequence for each sparrow, and the task sequence includes the variable values ​​of the number of modal components k and the penalty factor α; d is the number of variable values ​​to be optimized, and n is the number of sparrows in the initial population; Using envelope entropy as the fitness function value, calculate the fitness F of the initial population X and find the best and worst individuals; f() is the fitness function; Generate alarm values; Calculate and discover new positions of explorers, followers, and sentinels; If the new position is better than the current position, the position is updated until the iteration stop requirement is met. The global optimal solution in the updated sparrow population is determined as the optimal solution of the variational problem mathematical model of variational mode decomposition, and the optimal parameter combination and fitness value are output.

3. The method according to claim 2, characterized in that The variational mode decomposition of the input sample comprises: Constructing mathematical model of variational problem: In the formula, u k is the kth modal component obtained by decomposition; ω k for u k The center frequency of; δ(t) is the pulse function; k is the number of modal components; f is the spectrum of the original vibration signal; The alternating direction multiplier method is used to continuously update each modal component, center frequency and Lagrange multiplier and find the optimal solution. The updated expression is as follows: In the formula, α is the penalty factor, λ is the Lagrange multiplier, i represents the mode number; n represents the number of updates; ω represents the frequency coordinate in the spectrum; represents the estimated value of the Lagrange multiplier λ; represents the kth modal component u of the n+1th iteration k An estimated value of Represents the center frequency of the n+1th iteration The estimated value of .

4. The method according to claim 2, characterized in that: The explorer's location update includes: Where t is the current iteration number, j = 1, 2…, d; iter max is the maximum number of iterations, which is a constant; X i,j is the position information of the i-th sparrow in the j-th dimension; random number ξ∈[0,1]; warning value R2∈[0,1]; safety value ST∈[0.5,1]; q is a random number that follows a normal distribution; L is a 1×d matrix with all elements being 1; The follower's location update includes: Where, X p The best position for the current explorer; X worst is the worst position in the current population; A is a 1×d matrix column whose elements are randomly 1 or -1, and A + =a T (AA T ) -1 .

5. The method according to claim 1, characterized in that The calculation order tracking comprises: In the first step, the input sample vibration signal and the corresponding speed signal are sampled at a constant time interval Δt, and are recorded as x(n) and s(n) respectively; wherein s(n) is a sinusoidal voltage signal whose frequency changes with the speed, and each zero crossing point thereof corresponds to the time value of the rotating shaft rotating one circle; n represents a sampling point in the signal; The second step is to search for a point in s(n) that satisfies s(n)≤0 and s(n+1)≥0, and calculate the time value t corresponding to each rotation of the rotation axis by the following formula: zi : Among them, s(n) and s(n+1) correspond to time t s(n) and t s(n+1) The waveform between them is linear; Calculate the time series t corresponding to the equal angle increment Δθ through s(n) ik , i=1,2,3…, for t ik Calculate the difference to obtain the time T required for the rotating axis to rotate once i ,include: The initial velocity v of the first two rotations is calculated by the following formula: 10 And the angular acceleration a1: The time value t corresponding to each sampling point in the first two rotations is calculated by the following formula: ik , instantaneous speed v ik and the initial velocity v of the third cycle 30 : In the above formula, i is the number of rotations, i=1,2; k is the kth sampling point in each rotation; M is the number of equal-angle sampling points for each rotation of the rotating shaft, and 2kπ / M is the angle of rotation of the rotating shaft corresponding to the kth sampling point; Further deducing the angular acceleration a for i ≥ 3 cycles i and initial velocity v i+1,0 , and the corresponding t ik and v ik ; The third step is to calculate the ik , i = 1, 2, 3…, interpolate x(n) and calculate the value of x(n) at t ik The corresponding amplitude is obtained to obtain the equal-angle resampled signal x(t ik ).

6. The method according to claim 1, characterized in that The multi-scale subtraction network includes a feature extraction module, a multi-scale subtraction feature extraction module, a lightweight module, an attention mechanism, a global average pooling layer and a fully connected layer; The feature extraction module uses Res2Net-50 to extract features at four levels and reduces the feature channels to 64 using 3×3 convolution; The multi-scale subtraction feature extraction module takes the features of the four levels as input, obtains four complementary enhancement features, and outputs the four complementary enhancement features after fusion; The lightweight module uses the fused features output by the multi-scale subtraction feature extraction module and performs lightweight extraction; The channel attention mechanism gives more weight to important features in the lightweight extracted features; The global average pooling layer and the fully connected layer output diagnostic classification results.

7. A rolling bearing fault diagnosis system based on small sample data, characterized in that: include: An optimization module, which uses a variable working condition vibration signal as an input sample, optimizes the number of modal components k and the penalty factor α of the variational modal decomposition of the input sample by a sparrow search algorithm, and obtains an optimal parameter combination [k, α]; A variational mode decomposition module, which performs variational mode decomposition on the input sample based on an optimal parameter combination to obtain k intrinsic mode function components, and obtains the envelope order spectrum of each intrinsic mode function component by calculating order tracking and envelope demodulation; A signal reconstruction module, which selects the intrinsic mode function component with better sparse effect according to the envelope spectrum sparsity to perform signal reconstruction; An image diagnosis model building module converts the reconstructed signal into a two-dimensional image and divides it into a training set and a test set; the training set is input into a multi-scale subtraction network for model training to obtain a diagnosis model; A diagnosis module inputs the test set into the diagnosis model to perform bearing fault diagnosis.

8. The diagnostic system according to claim 7, characterized in that The optimizing the number k of modal components and the penalty factor α of the variational modal decomposition of the input sample by using the sparrow search algorithm includes: Set a sparrow population X, initialize a task sequence for each sparrow, and the task sequence includes the variable values ​​of the number of modal components k and the penalty factor α; d is the number of variable values ​​to be optimized, and n is the number of sparrows in the initial population; Using envelope entropy as the fitness function value, calculate the fitness F of the initial population X and find the best and worst individuals; f() is the fitness function; Generate alarm values; Calculate and discover new positions of explorers, followers, and sentinels; If the new position is better than the current position, the position is updated until the iteration stop requirement is met, and the global optimal solution in the updated sparrow population is determined as the optimal solution of the variational problem mathematical model of variational mode decomposition, and the optimal parameter combination and fitness value are output; The variational mode decomposition of the input sample comprises: Constructing mathematical model of variational problem: In the formula, u k is the u-th modal component obtained by decomposition; ω k for u k The center frequency of; δ(t) is the pulse function; k is the number of modal components; f is the spectrum of the original vibration signal; The alternating direction multiplier method is used to continuously update each modal component, center frequency and Lagrange multiplier and find the optimal solution. The updated expression is as follows: Where a is the penalty factor, λ is the Lagrange multiplier, i represents the mode number, n represents the number of updates, and ω represents the frequency coordinate in the spectrum. represents the estimated value of the Lagrange multiplier λ; represents the kth modal component u of the n+1th iteration k An estimated value of Represents the center frequency of the n+1th iteration The estimated value of .

9. The diagnostic system according to claim 7, characterized in that: The calculation order tracking comprises: In the first step, the input sample vibration signal and the corresponding speed signal are sampled at a constant time interval Δt, and are recorded as x(n) and s(n) respectively; wherein s(n) is a sinusoidal voltage signal whose frequency changes with the speed, and each zero crossing point thereof corresponds to the time value of the rotating shaft rotating one circle; n represents a sampling point in the signal; The second step is to search for a point in s(n) that satisfies s(n)≤0 and s(n+1)≥0, and calculate the time value t corresponding to each rotation of the rotation axis by the following formula: zi : Among them, s(n) and s(n+1) correspond to time t s(n) and t s(n+1) The waveform between them is linear; Calculate the time series t corresponding to the equal angle increment Δθ through s(n) ik , i=1,2,3…, for t ik Calculate the difference to obtain the time T required for the rotating axis to rotate once i ,include: The initial velocity v of the first two rotations is calculated by the following formula: 10 And the angular acceleration a1: The time value t corresponding to each sampling point in the first two rotations is calculated by the following formula: ik , instantaneous speed v ik and the initial velocity v of the third cycle 30 : In the above formula, i is the number of rotations, i=1,2; k is the kth sampling point in each rotation; M is the number of equal-angle sampling points for each rotation of the rotating shaft, and 2kπ / M is the angle of rotation of the rotating shaft corresponding to the kth sampling point; Further deducing the angular acceleration a for i ≥ 3 cycles i and initial velocity v i+1,0 , and the corresponding r ik and v ik ; The third step is to calculate the ik , i = 1, 2, 3…, interpolate x(n) and calculate the value of x(n) at t ik The corresponding amplitude is obtained to obtain the equal-angle resampled signal x(t ik ).

10. The diagnostic system according to claim 7, characterized in that The multi-scale subtraction network includes a feature extraction module, a multi-scale subtraction feature extraction module, a lightweight module, an attention mechanism, a global average pooling layer and a fully connected layer; The feature extraction module uses Res2Net-50 to extract features at four levels and reduces the feature channels to 64 using 3×3 convolution; The multi-scale subtraction feature extraction module takes the features of the four levels as input, obtains four complementary enhancement features, and outputs the four complementary enhancement features after fusion; The lightweight module uses the fused features output by the multi-scale subtraction feature extraction module and performs lightweight extraction; The channel attention mechanism gives more weight to important features in the lightweight extracted features; The global average pooling layer and the fully connected layer output diagnostic classification results.

Images