A method for extracting rolling bearing fault features under variable speed conditions

Abstract

The invention discloses a rolling bearing fault feature extraction method under variable speed working condition. The method constructs a Savitzky‑Golay (SG) high-pass filter convolution coefficient based on a Savitzky‑Golay (SG) filter, determines the window width of the SG high-pass filter convolution coefficient through an energy mean ratio index, and suppresses the DC component in a signal; secondly, a band-pass filter is constructed through the SG filter and the SG high-pass filter convolution coefficient, and the window width of the band-pass filter is obtained through an IDF index; finally, a band-pass filter is used to perform band-pass filtering on the signal after the DC component is eliminated through a band-pass filter, and the rolling bearing fault feature under variable speed working condition is revealed through an envelope order spectrum; the variable speed working condition fault feature extraction method proposed by the invention can effectively enhance the rolling bearing fault feature component, and has a good suppression effect on the DC component.

Description

Technical Field

[0001] The invention relates to a method for extracting rolling bearing fault features under variable speed working conditions, and belongs to the field of rotating machinery fault diagnosis. Background Art

[0002] There are a large number of variable speed working conditions in engineering practice. The processing of vibration signals in variable speed working conditions requires the calculation of speed information through algorithms such as order tracking, which brings certain difficulties to the fault diagnosis of rotating machinery with complex motion paths such as industrial robots. In recent years, the instantaneous angular speed (IAS) signal obtained by calculating the rotary encoder signal has gradually become a potential tool in the fault diagnosis of rotating machinery. Studies have shown that the IAS signal is not affected by the time-varying transmission path and can be directly obtained from the built-in encoder of the servo motor without the need for an external sensor. Therefore, rolling bearing fault feature extraction based on IAS signal is one of the hot topics in the field of rolling bearing fault diagnosis. There is a strong DC component in the variable speed IAS signal. A high-pass filter is usually used to filter out the interference of the DC component on fault feature extraction, but the parameter setting of the high-pass filter depends on experience; at the same time, the traditional fast spectral kurtosis (Fast Kurtogram, FK) algorithm is easily contaminated by noise, resulting in the failure of fault feature extraction. Summary of the invention

[0003] The present invention provides a method for extracting rolling bearing fault features under variable speed working conditions. The method constructs a Savitzky-Golay (SG) high-pass filter convolution coefficient based on a Savitzky-Golay (SG) filter, determines the window width of the SG high-pass filter convolution coefficient for suppressing a DC component through an energy mean ratio index, and suppresses the DC component in a signal; secondly, constructs a band-pass filter through the SG filter and the SG high-pass filter convolution coefficient, and obtains the window width of the band-pass filter through an improved diagnostic feature (IDF) index; finally, band-pass filtering is performed on the signal after the DC component is eliminated according to the band-pass filter, and the fault features of the rolling bearing under variable speed working conditions are revealed through envelope order spectrum analysis.

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

[0005] Step 1: Obtain the encoder instantaneous angular displacement information and corresponding time information through the PicoScope signal acquisition system, and calculate IAS using the forward difference method i The signal is expressed as:

[0006]

[0007] Where IAS irepresents the instantaneous angular velocity at the i-th moment, i = 1, 2, 3, ..., Δφ = 2π / N, Δt i =t i+1 -t i ; N represents the number of grating grids of the encoder.

[0008] Step 2: Eliminate the influence of DC component on fault feature extraction

[0009] 2-1. Constructing the convolution coefficient of Savitzky-Golay (SG) high-pass filter based on IAS i The signal is convolved to obtain RIAS i , for RIAS i Perform envelope analysis to obtain the envelope order spectrum E(RIAS i );E(RIAS) is characterized by the energy mean ratio index i ) in the DC component, the energy mean ratio index is expressed as:

[0010]

[0011] In the formula, E n is the energy average ratio index, E(RIAS i ) is RIAS i The envelope order spectrum of 1 and n 2 is the number of data points evaluated in the envelope order spectrum, n 1 =floor(ΔO -1 ), n 2 =floor(C / ΔO -1 ), floor(·) means rounding down; u=1,2,…,n 1 , v = n 1 ,n 1 +1,n 1 +2,…,n 2 ; ΔO is the order resolution, C is greater than 1× less than O reb The order of O reb is the theoretical characteristic order of the outer ring of the rolling bearing;

[0012] 2-2. Calculate the filter window width of the SG high-pass filter convolution coefficient corresponding to the maximum value of the energy mean ratio index, and obtain the filter window width of the SG high-pass filter convolution coefficient [2M+1] op , expressed as:

[0013] [2M+1] op = arg max[E n ]

[0014] In the formula, argmax[·] means returning En The corresponding SG high-pass filter convolution coefficient window width at maximum is [2M+1] op .

[0015] 2-3. According to [2M+1] op About IAS i Convolution is performed to obtain a signal with DC component eliminated [RIAS i ] op , the calculation formula is as follows:

[0016] [RIAS i ] op =IAS i -H*IAS i =(EH)*IAS i

[0017] Where, IAS i represents the original instantaneous angular velocity signal, H is the convolution coefficient of the SG filter, E is the unit matrix, EH is the convolution coefficient of the SG high-pass filter constructed by the SG filter, and * represents the convolution operation.

[0018] Step 3: Construct a bandpass filter by the convolution coefficient H of the SG filter and the convolution coefficient EH of the SG high-pass filter, expressed as:

[0019] R B =E([RIAS i ] op )·C oeff =E([RIAS i ] op )·I n (F re (H h ),F re (H)

[0020] H h =EH,H=(A T A) -1 A T ,A=[(B) 0 ,…,(B) P ] T

[0021] In the formula, E([RIAS i ] op ) is [RIAS i ] op The envelope order spectrum, C oeff and I n ((F re (H h )),(F re(H))) is the frequency response coefficient of the bandpass filter, R B is the envelope order spectrum after bandpass filtering, I n (a,b) means a and b intersect, F re (·) represents the frequency response curve calculation operation, H represents the convolution coefficient of the SG filter, and H h =EH is the convolution coefficient of the SG high-pass filter, (A T A) - 1 A T Indicates the calculation of the pseudo-inverse matrix of matrix A, A T is the transpose of matrix A, A -1 is the inverse of matrix A, B 0 =[-(M) 0 ,...,(-1) 0 ,1,1 0 ,...,(M) 0 ], B p =[-(M) p ,...,(-1) p ,0,1 p ,...,(M) p ], M is the half window length of the SG filter, and P is the filter order;

[0022] Step 4: Calculate the filter window width of the bandpass filter corresponding to the maximum value of the improved diagnostic feature (IDF) index, and obtain the optimized filter window width of the bandpass filter [2M h (k)+1,2M l (k)+1] op , expressed as:

[0023]

[0024] In the formula, O T =0.02O reb represents the maximum random slip of the rolling bearing, O b is the background noise width, R B is the envelope order spectrum after bandpass filtering, K represents the harmonic order of rolling bearing fault, argmax{·} represents the optimal window width of the bandpass filter corresponding to the maximum return IDF [2M h (k)+1,2M l (k)+1] op , 2M h (k)+1 is the window width of the SG high-pass filter convolution coefficient, 2M l (k)+1 is the window width of the SG filter convolution coefficient.

[0025] Step 5: Based on the bandpass filter window width determined in step 4, the signal [RIAS] that eliminates the DC component i ] op Perform K times of bandpass filtering and add the results of K times of bandpass filtering to realize fault feature extraction, where the bandpass filtering order P = 3:

[0026]

[0027] Where Y represents the fault feature extraction result.

[0028] The beneficial effects of the present invention are:

[0029] (1) The present invention proposes an energy mean ratio index, which is used to determine the window width of the convolution coefficient of the SG high-pass filter for eliminating the DC component, thereby solving the problem that the selection of the filter window width depends on experience when the high-pass filter suppresses the DC component;

[0030] (2) The present invention proposes to construct a bandpass filter by using the convolution coefficients of the SG filter and the SG high-pass filter, and determine the bandpass filter window width [2M] by using the IDF index. h (k)+1 and 2M l (k)+1] op , realizing the extraction of rolling bearing fault features under variable speed conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 The simulated IAS in Example 1 i Signal and its corresponding envelope order spectrum, where Figure (a) is the simulated IAS i Signal, Figure (b) is the simulated IAS i The envelope order spectrum corresponding to the signal;

[0032] Figure 2 The optimized window length [2M+1] for suppressing the DC component of the SG high-pass filter convolution coefficient calculated by the energy mean ratio index in Example 1 op ;

[0033] Figure 3 The simulated IAS for eliminating the DC component in Example 1 i Signal and envelope analysis results, where (a) is the simulated IAS with DC component removed i Signal, Figure (b) is the simulated IAS i Envelope order spectrum of the signal;

[0034] Figure 4The optimized window width of the bandpass filter calculated based on the IDF index in Example 1 and the feature extraction result of the simulation signal by the method of the present invention, wherein Figures (a)-(c) are the optimized window width of the bandpass filter when k=1, 2, 3, and Figure (d) is the feature extraction result of the simulation signal by the method of the present invention;

[0035] Figure 5 The experimental platform of Example 2;

[0036] Figure 6 is the measured IAS calculated by the forward difference method in Example 2 i Signal and corresponding envelope order spectrum, where (a) is the measured IAS i Signal, Figure (b) is the measured IAS i Envelope order spectrum of the signal;

[0037] Figure 7 The optimized window length [2M+1] for suppressing the DC component of the SG high-pass filter convolution coefficient calculated by the energy mean ratio index in Example 2 op ;

[0038] Figure 8 The measured IAS of Example 2 after eliminating the DC component i Signal and envelope analysis results, where (a) is the measured IAS with the DC component removed i Signal, Figure (b) is the measured IAS i Envelope order spectrum of the signal;

[0039] Fig. 9 In order to verify the effectiveness of the proposed energy mean index through Example 2, the convolution coefficient window width of the SG high-pass filter is set to 1001, and the envelope order spectrum with the DC component eliminated is obtained;

[0040] Fig.10 The optimized window width of the bandpass filter calculated based on the IDF index in Example 2 and the feature extraction result of the measured signal by the method of the present invention, wherein Figures (a)-(c) are the optimized window width of the bandpass filter when k=1, 2, 3, and Figure (d) is the feature extraction result of the simulation signal by the method of the present invention;

[0041] Fig.11 The figure shows the result of analyzing the actual signal of the rolling bearing outer race fault with the DC component eliminated in Example 2 using the Fast Kurtogram (FK) algorithm, wherein Figure (a) is a demodulated frequency band containing rich bearing fault information, and Figure (b) is an order spectrum corresponding to the optimized demodulated frequency band. DETAILED DESCRIPTION

[0042] The following is a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. The methods not specially described in the embodiments are all conventional methods.

[0043] The following description sets forth many specific details for a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0044] Embodiment 1: Simulating the outer ring fault of a rolling bearing, and using the method of the present invention to extract the outer ring fault characteristics of the rolling bearing under variable speed conditions, the specific process includes:

[0045]

[0046]

[0047]

[0048] In the formula, w b The variable speed simulation signal indicating the outer race fault, is the outer race fault component, is the frequency conversion component, A u is the disturbance amplitude of spalling, is the angular position sequence, Θ is the average angular interval between adjacent fault disturbances, sin(·) is the sine function, and the outer ring fault order is O reb =2π / Θ,χ u represents the random sliding angle, O n is the fixed frequency of the system, O r1 is the rotation frequency of the large gear shaft, O r2 is the rotation frequency of the pinion shaft, ξ is the attenuation coefficient, is Gaussian white noise, is the initial phase, U is the number of IAS mutations caused by rolling bearing failure, and the speed change trend is simulated by the first quarter cycle of the sine function sin, and the speed increases from 1rad / s to 4rad / s; r1 =1×;O r2 =2.33×; Θ=2π / O reb ;O reb =3.1; A u =0.1;O n=0.001×;ξ=0.03;θ=2π / N, 4π / N, 6π / N...;

[0049] Step 1: Simulated variable speed outer race fault IAS i Signals such as Figure 1 (a) shows, its corresponding envelope order spectrum is as follows Figure 1 As shown in (b), it can be seen that there is a strong DC component in the envelope order spectrum, and the fault characteristics are completely submerged by the DC component, so the DC component needs to be eliminated.

[0050] Step 2: Eliminate the influence of DC component on fault feature extraction

[0051] 2-1. Constructing the convolution coefficient of Savitzky-Golay (SG) high-pass filter based on IAS i The signal is convolved to obtain RIAS i , for RIAS i Perform envelope analysis to obtain the envelope order spectrum E(RIAS i ); The energy mean ratio index is used to characterize the degree of the DC component in E(RIAS), and the energy mean ratio index is expressed as:

[0052]

[0053] In the formula, E n is the energy average ratio index, E(RIAS i ) indicates RIAS i The envelope order spectrum of 1 and n 2 are the number of data points evaluated in the envelope order spectrum; n 1 =floor(ΔO -1 ),n 2 =floor(C / ΔO -1 ), floor(·) means rounding down, ΔO is the order resolution, where C = 2;

[0054] 2-2. Calculate the filter window width of the SG high-pass filter convolution coefficient corresponding to the maximum value of the energy mean ratio index, and obtain the filter optimization window width of the SG high-pass filter convolution coefficient [2M+1] op , expressed as:

[0055] [2M+1] op = arg max[E n ]

[0056] In the formula, argmax[·] means returning E n The corresponding SG high-pass filter convolution coefficient window width is the maximum; Figure 2As shown, the SG high-pass filter convolution coefficient optimizes the window width [2M+1] op =5037;

[0057] 2-3. Set the filter order P=3 and filter optimization window width [2M+1] of the SG high-pass filter convolution coefficient op =5037 pairs of IAS i The signal is convolved to obtain a signal with DC component eliminated [RIAS i ] op , the calculation formula is as follows:

[0058] [RIAS i ] op =IAS i -H*IAS i =(EH)*IAS i

[0059] Where, IAS i represents the original instantaneous angular velocity signal, H is the convolution coefficient of the SG filter, E is the unit matrix, EH is the convolution coefficient of the SG high-pass filter constructed by the SG filter, and * represents the convolution operation. The signal with the DC component eliminated [RIAS i ] op like Figure 3 (a) shows, the corresponding envelope order spectrum is as follows Figure 3 As shown in (b), it can be seen that the DC component is effectively suppressed.

[0060] Step 3: Construct a bandpass filter by the convolution coefficient H of the SG filter and the convolution coefficient EH of the SG high-pass filter, expressed as:

[0061] R B =E([RIAS i ] op )·C oeff =E([RIAS i ] op )·I n (F re (H h ),F re (H)

[0062] H h =EH,H=(A T A) -1 A T ,A=[(B) 0 ,…,(B) P ] T

[0063] In the formula, E([RIAS i ] op) is [RIAS i ] op The envelope order spectrum, C oeff and I n ((F re (H h )),(F re (H l ))) is the frequency response coefficient of the bandpass filter, R B is the envelope order spectrum after bandpass filtering, I n (a,b) means a and b intersect, F re (·) represents the frequency response curve calculation operation, H represents the convolution coefficient of the SG filter, and H h =EH is the high-pass filter convolution coefficient, (A T A) -1 A T Indicates the calculation of the pseudo-inverse matrix of matrix A, A T is the transpose of matrix A, A -1 is the inverse of matrix A, B 0 =[-(M) 0 ,...,(-1) 0 ,1,1 0 ,...,(M) 0 ], B p =[-(M) p ,...,(-1) p ,0,1 p ,...,(M) p ], M is the half window length of the SG filter, and P is the filter order.

[0064] Step 4: Calculate the filter window width of the bandpass filter corresponding to the maximum value of the improved diagnostic feature (IDF) index, and obtain the optimized filter window width of the bandpass filter [2M h (k)+1,2M l (k)+1] op , expressed as:

[0065]

[0066] In the formula, O T =0.02O reb represents the maximum random slip of the rolling bearing, O b is the background noise width, R B is the envelope order spectrum after bandpass filtering, K represents the harmonic order of rolling bearing fault, argmax{·} represents the optimal window width of the bandpass filter corresponding to the maximum return IDF [2M h (k)+1,2M l (k)+1]op When k = 1, the optimized window width of the bandpass filter is [5055, 3801], such as Figure 4 (a) As shown; when k = 2, the optimized window width of the bandpass filter is [1493, 2351], as shown in Figure 4 (b) As shown; when k = 3, the optimized window width of the bandpass filter is [1059, 1631], as shown in Figure 4 (c) as shown.

[0067] Step 5: Based on the bandpass filter determined in step 4, the optimized window width is used to eliminate the DC component of the signal [RIAS i ] op Perform K times of bandpass filtering and add the results of K times of bandpass filtering to realize fault feature extraction, where the bandpass filtering order P = 3 and the number of bandpass filtering times is K = 3:

[0068]

[0069] Where Y represents the simulated variable speed rolling bearing outer ring fault feature extraction result, such as Figure 4 As shown in (d), the characteristic order of the rolling bearing outer ring fault is O reb , 2O reb and 3O reb Get effective extraction;

[0070] Embodiment 2: This embodiment describes the process of applying the method of the present invention to extract the fault characteristics of the outer ring of a rolling bearing under actual variable speed conditions:

[0071] Use Figure 5 The rolling bearing test bench shown in the figure is used for verification. In this embodiment, the motor speed is accelerated from 0 to 40r / min, the number of teeth of the small gear in the gearbox is 30, the number of teeth of the large gear is 70, the reduction ratio is 2.33, the number of lines of the encoder is 5000, and it is installed on the large gear shaft; the bearing model used in the experiment is 6204 deep groove ball bearing, the rolling element diameter d is 7.54mm, the pitch circle diameter D is 33.5mm, the number of rolling elements Z is 8, and the contact angle α is 0 degree; in order to simulate the outer ring fault of the rolling bearing, the outer ring of the bearing is processed by the electric spark machining method with a width of about 0.2mm, and the outer ring fault order O of the rolling bearing is calculated by the following calculation formula reb =3.1×.

[0072]

[0073] In the formula, O reb It represents the characteristic order of outer ring fault, Z represents the number of rolling elements, d represents the rolling element diameter, D represents the pitch circle diameter, and α represents the contact angle.

[0074] Step 1: Obtain the encoder instantaneous angular displacement information and corresponding time information through the PicoScope signal acquisition system, and calculate IAS using the forward difference method i The signal is expressed as:

[0075]

[0076] Where, IAS i represents the instantaneous angular velocity at the i-th moment, Δt i =t i+1 -t i , N is the number of gratings of the encoder N = 5000. IAS i Signals such as Figure 6 As shown in (a), Figure 6 As shown in (b), there is a strong DC component in the corresponding envelope order spectrum, which has a great impact on the extraction of rolling bearing fault features;

[0077] Step 2: Eliminate the influence of DC component on fault feature extraction.

[0078] 2-1. Constructing the convolution coefficient of Savitzky-Golay (SG) high-pass filter based on IAS i The signal is convolved to obtain RIAS i , for RIAS i Perform envelope analysis to obtain the envelope order spectrum E(RIAS i );E(RIAS) is characterized by the energy mean ratio index i ) in the DC component, the energy mean ratio index is expressed as:

[0079]

[0080] In the formula, E n is the energy average ratio index, E(RIAS i ) indicates RIAS i The envelope order spectrum of 1 and n 2 are the number of data points evaluated in the envelope order spectrum; n 1 =floor(ΔO -1 ),n 2 =floor(C / ΔO -1 ), floor(·) means rounding down, ΔO is the order resolution, where C = 2;

[0081] 2-2. Calculate the filter window width of the SG high-pass filter convolution coefficient corresponding to the maximum value of the energy mean ratio index, and obtain the filter optimization window width of the SG high-pass filter convolution coefficient [2M+1] op, expressed as:

[0082] [2M+1] op = argmax[E n ]

[0083] In the formula, argmax[·] means returning E n The corresponding SG high-pass filter convolution coefficient window width is the maximum, such as Figure 7 As shown, the SG high-pass filter convolution coefficient optimizes the window width [2M+1] op =4628;

[0084] 2-3. Set the filter order P=3 and the optimized window width [2M+1] of the SG high-pass filter convolution coefficient op =4628 pairs of IAS i The signal is convolved to obtain a signal with DC component eliminated [RIAS i ] op , the calculation formula is the same as in Example 1, the signal with DC component eliminated [RIAS i ] op like Figure 8 (a) shows, the corresponding envelope order spectrum is Figure 8 As shown in (b), it can be seen that the DC component is effectively suppressed;

[0085] According to the traditional method, the window width of the SG high-pass filter convolution coefficient is selected as 1001 based on experience, and the signal with DC component eliminated is obtained and its envelope order spectrum is calculated, such as Fig. 9 As shown, it can be seen that the characteristic order of the rolling bearing outer ring fault is O reb , 2O reb and 3O reb It is obviously suppressed, which means that the selection of filter window width based on experience is not only time-consuming and labor-intensive, but also easily leads to the loss of fault characteristic information;

[0086] Step 3: construct a bandpass filter by using the convolution coefficient H of the SG filter and the convolution coefficient EH of the SG high-pass filter, and the construction method is the same as that of Example 1;

[0087] Step 4: Calculate the filter window width of the bandpass filter corresponding to the maximum value of the improved diagnostic feature (IDF) index, and obtain the optimized filter window width of the bandpass filter [2M h (k)+1,2M l (k)+1] op , expressed as:

[0088]

[0089] In the formula, O T =0.02Oreb represents the maximum random slip of the rolling bearing, O b is the background noise width, R B is the envelope order spectrum after bandpass filtering, K represents the harmonic order of rolling bearing fault, argmax{·} represents the optimal window width of the bandpass filter corresponding to the maximum return IDF [2M h (k)+1,2M l (k)+1] op When k = 1, the optimized window width of the bandpass filter is [4401, 3601], such as Fig.10 (a) As shown; when k = 2, the optimized window width of the bandpass filter is [1613, 2451], as shown in Fig.10 (b) As shown; when k = 3, the optimized window width of the bandpass filter is [919, 1451], as shown in Fig.10 (c) as shown.

[0090] Step 5: Based on the bandpass filter determined in step 4, the optimized window width is used to eliminate the DC component of the signal [RIAS i ] op Perform K times of bandpass filtering and add the results of K times of bandpass filtering to realize fault feature extraction, where the bandpass filtering order P = 3 and the number of bandpass filtering times is K = 3:

[0091]

[0092] Where Y represents the actual variable speed rolling bearing outer ring fault feature extraction result; Fig.10 As shown in (d), the characteristic order of the rolling bearing outer ring fault is O reb , 2O reb and 3O reb Effectively extracted.

[0093] In order to further verify the effectiveness of the bandpass filter constructed in the present invention in fault feature extraction, the Fast Kurtogram (FK) algorithm is used to determine the demodulation frequency band containing rich bearing fault information in the above embodiment signal. The results are as follows: Fig.11 As shown in (a), the center frequency is 712.89×, the bandwidth is 19.53×, and the order spectrum corresponding to the demodulation band is as follows Fig.11 As shown in (b), the characteristic order of the rolling bearing outer ring fault is O reb , 2O reb and 3O reb Unable to identify effectively.

Claims

1. A method for extracting rolling bearing fault features under variable speed conditions, characterized in that: Here are the steps: (1) Obtain the instantaneous angular displacement information and corresponding time information of the optical encoder, and use the forward difference method to calculate the instantaneous angular velocity IAS under variable speed conditions i Signal; (2) Constructing the convolution coefficient of SG high-pass filter based on SG filter for IAS i Perform convolution to obtain RIAS i ; for RIAS i Perform envelope analysis to obtain the envelope order spectrum E(RIAS i );E(RIAS) is characterized by the energy mean ratio index i ) and obtain the window width of the SG high-pass filter convolution coefficient [2M+1] op ; Based on [2M+1] op About IAS i Convolution is performed to obtain a signal with DC component eliminated [RIAS i ] op ; (3) Construct a bandpass filter through the convolution coefficient H of the SG filter and the convolution coefficient EH of the SG high-pass filter, calculate the filter window width of the bandpass filter corresponding to the maximum value of the IDF index, and obtain the filter window width of the bandpass filter [2M h (k)+1,2M l (k)+1] op , based on [2M h (k)+1,2M l (k)+1] op To eliminate the DC component of the signal [RIAS i ] op Perform K times of bandpass filtering, where k = 1, 2, ... K, and the bandpass filter is expressed as: R B =AND([RIAS i ] op )·W oeff =AND([RIAS i ] op )·I n (F re (H h ),F re (H)) H h =E-H,H=(A T A) -1 A T ,A=[(B) 0 ,…,(B) P ] T In the formula, E([RIAS i ] op ) is [RIAS i ] op The envelope order spectrum, C oeff and I n ((F re (H h )),(F re (H))) is the frequency response coefficient of the bandpass filter, R B is the envelope order spectrum after bandpass filtering, I n (a,b) means a and b intersect, F re (·) represents the frequency response curve calculation operation, H represents the convolution coefficient of the SG filter, and H h =EH is the convolution coefficient of the SG high-pass filter, (A T A) -1 A T Indicates the calculation of the pseudo-inverse matrix of matrix A, A T is the transpose of matrix A, A -1 is the inverse of matrix A, B 0 =[-(M) 0 ,...,(-1) 0 ,1,1 0 ,...,(M) 0 ], B p =[-(M) p ,...,(-1) p ,0,1 p ,...,(M) p ], M is the half window length of the SG filter, and P is the filter order; (4) Add the results of K times of bandpass filtering to realize rolling bearing fault feature extraction, which is expressed as: Where Y represents the fault feature extraction result.

2. The method for extracting rolling bearing fault features under variable speed conditions according to claim 1 is characterized in that: The energy mean ratio (RIAS) is used to characterize i ) to obtain the window width of the SG high-pass filter convolution coefficient [2M+1] op : In the formula, E n is the energy average ratio index, E(RIAS i ) is RIAS i The envelope order spectrum of , where n1 and n2 are the number of data points evaluated in the envelope order spectrum, n1 = floor (ΔO -1 ), n2 = floor(C / ΔO -1 ), floor(·) means rounding down; u=1,2,…,n1, v=n1,n1+1,n1+2,…,n2; ΔO is the order resolution, C is greater than 1× less than O reb The order of O reb is the theoretical characteristic order of the outer ring of the rolling bearing; argmax[·] means returning E n The window width of the SG high-pass filter convolution coefficient corresponding to the maximum value is [2M+1] op .

3. The method for extracting rolling bearing fault features under variable speed conditions according to claim 1, characterized in that: According to [2M+1] op About IAS i Convolution is performed to obtain a signal with DC component eliminated [RIAS i ] op , the calculation formula is as follows: [RIAS i ] op =IAS i -H*IAS i =(E-H)*IAS i In the formula, IAS i represents the original instantaneous angular velocity signal, H is the convolution coefficient of the SG filter, E is the unit matrix, EH is the convolution coefficient of the SG high-pass filter constructed by the SG filter, and * represents the convolution operation.

4. The method for extracting rolling bearing fault features under variable speed conditions according to claim 1, characterized in that: Calculate the filter window width of the bandpass filter corresponding to the maximum value of the IDF index, and obtain the filter window width of the bandpass filter [2M h (k)+1,2M l (k)+1] op , the calculation formula is as follows: In the formula, O T =0.02O reb represents the maximum random slip of the rolling bearing, O b is the background noise width, O reb is the theoretical characteristic order of the outer ring of the rolling bearing, R B is the envelope order spectrum after bandpass filtering, k = 1, 2, ... K, argmax{·} represents the corresponding bandpass filter window width when the IDF is the largest [2M h (k)+1,2M l (k)+1] op , 2M h (k)+1 is the window width of the SG high-pass filter convolution coefficient, 2M l (k)+1 is the window width of the SG filter convolution coefficient.

Images