Keyless phase order tracking method for rolling bearing fault diagnosis under highly variable speed

Abstract

The present invention discloses a keyless phase order tracking method for rolling bearing fault diagnosis under strongly variable speed, which reduces noise of the vibration signal of the rolling bearing under strongly variable speed working condition, and then adaptively determines the optimal time-varying frequency modulation frequency of each frequency band of the signal by maximizing the local matching kurtosis operator, and then obtains the high-resolution time-frequency distribution of the multi-component signal; redistributes the energy to obtain the target ridge line with minimum loss; according to the target ridge line with minimum loss, extracts the time domain waveform of the high-order harmonics of the instantaneous rotation frequency, and obtains the phase information of the instantaneous speed; uses the phase information of the instantaneous speed to perform angular domain resampling and envelope demodulation analysis on the filtered signal, and obtains the fault characteristic order diagram of the bearing. The present invention can overcome the problems of low aggregation of time-frequency distribution and fuzzy instantaneous speed estimation caused by the strong frequency modulation characteristics of the signal, and realizes accurate estimation of the fundamental frequency and multiple frequency of the rolling bearing under strongly variable speed conditions.

Description

Technical Field

[0001] The invention belongs to the field of rolling bearing fault diagnosis under mechanical vibration signals and strongly variable speed working conditions, and in particular relates to a keyless phase order tracking method for rolling bearing fault diagnosis under strongly variable speed. Background Art

[0002] Rolling bearings are an important part of mechanical equipment. Their operating reliability and working condition directly determine whether the mechanical equipment can operate safely and smoothly for a long time. Most traditional rolling bearing fault diagnosis methods are based on constant speed and weak speed change conditions. Fault diagnosis methods for rolling bearings under strong speed change conditions are less studied. However, in actual engineering, rolling bearings often work under conditions with drastic speed fluctuations. For example, the speed of the bearings in the turbo pump of a rocket engine system accelerates from stop to rated speed in just a few seconds during the start-up and shutdown phases. The operating condition of the bearings in these few seconds directly determines the success of the launch. In addition, since most mechanical equipment has a compact structure and little spare space, it is difficult to directly use speed sensors to measure the speed. Therefore, the research on the keyless phase order tracking method of rolling bearings under strong speed change conditions is of great significance to the safe and reliable operation of equipment.

[0003] Keyless phase order tracking technology is currently the main method to solve the analysis of time-varying speed signals when the speed cannot be directly obtained. The core of keyless phase order tracking technology lies in the extraction of instantaneous rotational frequency of variable speed signals. The quality of instantaneous rotational frequency extraction directly determines the final order tracking result; the commonly used keyless phase order tracking technology is difficult to overcome the problems of too low aggregation of time-frequency distribution and fuzzy instantaneous rotational frequency caused by the strong frequency-amplitude modulation characteristics of strong variable speed signals, which leads to the extracted fault feature order being not accurate enough. Therefore, in view of the deficiency of existing keyless phase order tracking that it is difficult to extract the instantaneous rotational frequency of strong frequency-amplitude modulation signals, how to improve the accuracy of signal time-frequency distribution and accurately extract the instantaneous rotational frequency is a problem to be solved. Summary of the invention

[0004] The purpose of the present invention is to overcome the problems that the existing time-frequency analysis technology has difficulty in realizing high-precision time-frequency analysis of strongly variable speed signals and instantaneous frequency ambiguity, and proposes a keyless phase order tracking method for rolling bearing fault diagnosis under strongly variable speed. This method can accurately extract the fundamental frequency and its multiples of the strongly variable speed signal without a keyless phase signal, and realize rolling bearing fault diagnosis under strongly variable speed conditions.

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

[0006] A keyless phase order tracking method for rolling bearing fault diagnosis under strongly variable speed comprises the following steps:

[0007] (1) De-noise the vibration signal of the rolling bearing under strong variable speed conditions to obtain the low-frequency characteristics of the signal;

[0008] (2) The optimal time-varying frequency modulation frequency of each frequency band of the signal is adaptively determined by maximizing the local matching kurtosis operator based on the low-frequency characteristics of the signal, and then the high-resolution time-frequency distribution of the multi-component signal is obtained by matching Chirplet transform;

[0009] (3) Taking the local maximum value in the frequency direction as the center, the energy is redistributed using the local maximum synchronous compression transform operator;

[0010] (4) Using a ridge extraction algorithm based on a cost function, the target ridge with the minimum loss is obtained according to the global optimal time-frequency distribution after energy redistribution;

[0011] (5) according to the target ridge line with minimum loss obtained in step (4), extract the time domain waveform of the higher harmonics of the instantaneous rotation frequency to obtain the phase information of the instantaneous rotation speed;

[0012] (6) bandpass filtering the vibration signal of the rolling bearing under the strong variable speed condition using the center frequency and the optimal bandwidth to obtain a filtered signal;

[0013] (7) Using the phase information of the instantaneous speed obtained in step (5), the filtered signal in step (6) is resampled in the angular domain and subjected to envelope demodulation analysis to obtain a fault characteristic order diagram of the bearing.

[0014] Furthermore, the specific process of step (2) is as follows:

[0015] According to the variable frequency modulation frequency c, the low-frequency characteristic signal of the signal obtained in step (1) is matched with Chirplet transformation to obtain the time-frequency distribution G(c, t, ω) under different time-varying frequency modulation frequencies c:

[0016]

[0017] Where t′ is time, ω is frequency, g(ut′) is Gaussian window function, s(u) is discrete signal, i is imaginary sign, u is independent variable, -c(ut′) 2 / 2 is the demodulation operator;

[0018] Then, set the discrete number of frequency bands to N c , along the frequency direction, the signal at t c The spectrum range at the moment is divided into N c equal-spaced frequency bands, according to the local matching kurtosis operator LK(c,t c ,N c ), calculate the signal t c The local matching kurtosis LK(c,tc ,N c );

[0019] Finally, by maximizing the signal at t c The local matching kurtosis operator LK(c,t c ,N c ), and obtain t c The local time-varying frequency sequence of Get the global optimal time-frequency distribution

[0020] Furthermore, the signal t c Time i n The local matching kurtosis of a frequency band generated by different time-varying FM frequencies is calculated by the following formula:

[0021]

[0022]

[0023]

[0024]

[0025]

[0026] In the formula, t c The signal at time i n The optimal time-varying FM frequency for each frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the first frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the second frequency band, t c The time signal is at the Nth c The optimal time-varying FM frequency of the frequency band;

[0027] Furthermore, the specific process of step (3) is as follows:

[0028] First, for the high-resolution time-frequency distribution of the multi-component signal obtained in step (2), the local maximum value in the frequency direction is selected as the center, and the spectrum is divided into several frequency bands according to the bandwidth that meets the separation conditions of different instantaneous frequency signals:

[0029] Then, using the local maximum synchronization compression operator ω m (c, t, ω), with the local maximum as the center, the high-resolution time-frequency distribution of the multi-component signal obtained in step (2) is redistributed along the frequency direction to obtain the global optimal time-frequency distribution after energy redistribution.

[0030] Furthermore, the frequency band bandwidth that satisfies the conditions for separating signals of different instantaneous frequencies satisfies the following formula:

[0031]

[0032] Among them, Δ is the frequency bandwidth that satisfies the separation conditions of different instantaneous frequency signals, is the instantaneous frequency of the k+1th component, is the instantaneous frequency of the kth component;

[0033] The global optimal time-frequency distribution after energy redistribution is as follows:

[0034]

[0035]

[0036] Among them, ω m (c, t, ω) is the local maximum synchrosqueezing operator, δ is the delta function, and η represents the frequency domain coordinate after energy redistribution.

[0037] Furthermore, the specific process of step (4) is as follows:

[0038] The point with the largest energy in the global optimal time-frequency distribution after energy redistribution is searched as the initial point of the time-frequency ridge, and the recursive update formula of the cost function is obtained according to the initial point;

[0039] Finally, according to the initial point, the maximum frequency transformation range and the recursive update formula, the minimum cost function point is searched on both sides of the initial point along the time axis to obtain the target ridge line.

[0040] Furthermore, the recursive update formula of the cost function is as follows:

[0041]

[0042]

[0043] Among them, j is the independent variable, t m is the horizontal coordinate of the maximum time-frequency energy point, m is t m The corresponding coordinates, ω m is the ordinate of the maximum time-frequency energy point, ω(t j ) is t j The ordinate corresponding to the maximum time-frequency energy point at time T s =F s / 2 is the frequency domain range of the time-frequency distribution; C1 is the amplitude transformation weight, C2 is the curve gradient weight, Cost is the cost function, is the global optimal time-frequency distribution after energy redistribution.

[0044] Furthermore, the specific process of step (5) is as follows: the target ridge line with the minimum loss is used as the center frequency of the Vold-Kalman filter, the instantaneous bandwidth of the filter is set, and the filtered signal corresponding to the high-order harmonics is obtained; then the Hilbert transform, inverse tangent function and phase deconvolution are used to obtain the phase information of the instantaneous rotation speed.

[0045] Furthermore, the center frequency and the optimal bandwidth in step (6) are determined by the following process:

[0046] The strongly variable speed signal is decomposed by using a filter with phase quasi-analysis effect in the form of double downsampling and iterative filtering, and the spectral kurtosis SK(f) of each sub-band is calculated;

[0047]

[0048] Where S(t,f) represents the time-frequency distribution of signal s;

[0049] The spectral kurtosis of all sub-signals is compared to determine the center frequency and optimal bandwidth corresponding to the maximum spectral kurtosis.

[0050] Furthermore, the specific process of step (7) is as follows:

[0051] First, the least squares fitting method is used to establish the mapping relationship between angle and time.

[0052] Secondly, determine the size of the angle domain sampling interval Δθ and find the time corresponding to the equal angle domain sampling point kΔθ Using the cubic spline interpolation method, according to the time corresponding to the sampling point kΔθ in the equiangular domain resampling the filtered signal;

[0053] Finally, the envelope demodulation of the angular domain resampled signal machine is performed to obtain the order diagram that characterizes the rolling bearing fault characteristics;

[0054] Among them, the angle domain sampling interval and the time corresponding to the equal angle domain sampling point kΔθ are as follows:

[0055]

[0056]

[0057] In the formula, F s is the sampling frequency, f(t) is the minimum rotation frequency of the reference axis, and k is the coefficient.

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

[0059] The present invention first filters and reduces noise on the vibration signal of a rolling bearing under a strongly variable speed condition, and secondly, adaptively determines the optimal time-varying frequency modulation frequency of each frequency band of the signal through a local matching kurtosis operator, and obtains a high-resolution time-frequency distribution using a matching Chirplet transform; subsequently, a local maximum synchronous compression transform operator is used to rearrange energy along the frequency direction; after obtaining the target ridge line through a ridge extraction algorithm based on a cost function, the instantaneous frequency of the speed is extracted; finally, the time domain waveform of the higher-order harmonics of the instantaneous rotation frequency is used to obtain the phase information of the instantaneous speed; and the phase information of the instantaneous speed is used to perform angular domain resampling and envelope demodulation on the filtered signal to extract the fault feature order. The present invention can realize order tracking processing without a key phase trigger control circuit, can effectively overcome the instantaneous frequency ambiguity problem caused by strong speed change, can obtain high-resolution time-frequency distribution, can accurately extract the instantaneous shaft frequency ridge under the starting condition, and realize accurate estimation of the fundamental frequency and multiple frequencies of the rolling bearing fault characteristics under the condition of strong speed change. It has strong robustness under strong background noise and adjacent frequency interference, and can be used in various working conditions such as starting and stopping, rapid speed increase and speed decrease. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0061] Figure 1 A flow chart of a keyless phase order tracking method for rolling bearing fault diagnosis under strongly variable speed conditions proposed by the present invention;

[0062] Figure 2 This is a vibration waveform diagram of a bearing under a strongly variable speed working condition according to an embodiment of the present invention;

[0063] Figure 3 The speed curve diagram of the bearing under the condition of strong variable speed in the embodiment of the present invention;

[0064] Figure 4 It is a high-precision time-frequency distribution diagram of the fourth harmonic of the frequency conversion extracted by the present invention; wherein (a) is a time-frequency distribution diagram within the range where the fourth harmonic of the frequency conversion is located, and (b) is a local enlarged diagram of the high-precision time-frequency distribution;

[0065] Figure 5 This is a diagram of the frequency conversion information estimated by the local matching synchronous compression transformation method proposed in the present invention; wherein (a) is the estimated frequency conversion, and (b) is the estimated error of the frequency conversion;

[0066] Figure 6Figure 2 is a diagram of the frequency conversion information obtained by synchronous compression transformation estimation; (a) is the estimated frequency conversion, and (b) is the estimated error of the frequency conversion;

[0067] Figure 7 A diagram showing the results of the fast spectral kurtosis index analysis implemented by the present invention;

[0068] Figure 8 The keyless phase order tracking result after selecting the sensitive frequency band by the method proposed by the present invention;

[0069] Fig. 9 Results of keyless phase-order tracking using the rotation frequency estimated by the synchrocompressed transform. DETAILED DESCRIPTION

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

[0071] Reference Figure 1 As shown, the present invention provides a keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed, comprising the following steps:

[0072] (1) Extraction of low-frequency signal features: First, the vibration signal of the rolling bearing under strong variable speed conditions is measured by an acceleration sensor;

[0073] Then, with the transfer frequency f c With the center as the center, appropriate bandwidth parameters are selected, and a bandpass filter with zero phase offset characteristics is used to reduce the noise of the collected rolling bearing vibration signal to obtain the low-frequency characteristics of the signal;

[0074] (2) Local kurtosis matching Chirplet transform: The optimal time-varying frequency modulation frequency c of each frequency band of the signal is adaptively determined by maximizing the local matching kurtosis operator for the low-frequency characteristics of the signal, and the high-resolution time-frequency distribution of the multi-component signal is obtained by matching Chirplet transform. The specific process is as follows:

[0075] First, a series of discrete time-varying FM frequencies c are used to describe the time-varying AM FM characteristics in Chirplet transform; let the number of discrete points be N a By discretizing the rotation angle at equal intervals in the range of (-π / 2,π / 2), the relationship between the rotation angle α and the discretized candidate time-varying frequency c is obtained:

[0076]

[0077]

[0078] Among them, F s is the sampling frequency, T sis the signal length, g(ut′) is the Gaussian window function, and c is the time-varying frequency modulation frequency used to match the local time-frequency characteristics of the signal;

[0079] Then, according to a series of candidate time-varying FM frequencies c constructed, the low-frequency characteristic signal of the signal obtained in step (1) is matched with Chirplet transform to obtain the time-frequency distribution G(c, t, ω) under different time-varying FM frequencies c:

[0080]

[0081] Where t′ is time, ω is frequency, g(ut′) is Gaussian window function, s(u) is discrete signal, i is imaginary sign, u is independent variable, -c(ut′) 2 / 2 is the demodulation operator; when the time-varying frequency modulation frequency c is 0, the demodulation operator is set to 0, and the matching Chirplet transform degenerates into the traditional short-time Fourier transform;

[0082] Then, set the discrete number of frequency bands to N c , along the frequency direction, the signal at t c The spectrum range at the moment is divided into N c equal-spaced frequency bands, according to the local matching kurtosis operator LK(c,t c ,N c ), calculate the signal t c The local matching kurtosis LK(c,t c ,N c ), with c as the i-th demodulation factor n The expression of the local matching kurtosis operator of a frequency band is:

[0083]

[0084]

[0085] Among them, i is the independent variable, i n Indicates the i n frequency band.

[0086] Finally, by maximizing the signal at t c The local matching kurtosis operator LK(c,t c ,N c ), you can get t c The local time-varying frequency sequence of So as to obtain the global optimal time-frequency distribution

[0087]

[0088]

[0089]

[0090] In the formula, t c The signal at time i n The optimal time-varying FM frequency for each frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the first frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the second frequency band, t c The time signal is at the Nth c The optimal time-varying FM frequency for the frequency band.

[0091] (3) Local maximum matching synchronous compression transform: With the local maximum value in the frequency direction as the center, the local maximum synchronous compression transform operator is used to redistribute the energy to improve the robustness of the effective features in the time-frequency domain under strong interference; the specific process is as follows: First, in order to redistribute the low-aggregation time-frequency distribution energy along the frequency direction to the theoretical frequency ridge, for the high-resolution time-frequency distribution of the multi-component signal obtained in step (2), the local maximum value in the frequency direction is selected as the center, and the spectrum is divided into several frequency bands with Δ as the frequency band bandwidth:

[0092]

[0093] Among them, Δ is the frequency bandwidth that satisfies the separation conditions of different instantaneous frequency signals, is the instantaneous frequency of the k+1th component, is the instantaneous frequency of the kth component;

[0094] Then, using the local maximum synchronization compression operator ω m (c, t, ω), with its local maximum as the center, redistribute the energy of the time-frequency distribution obtained in step (2) along the frequency direction to obtain the global optimal time-frequency distribution after energy redistribution:

[0095]

[0096]

[0097] (4) Time-frequency ridge extraction based on cost function: A ridge extraction algorithm based on cost function is used to obtain the target ridge with minimum loss according to the global optimal time-frequency distribution after energy redistribution, so as to solve the discontinuity problem of instantaneous frequency conversion ridge caused by noise. The specific process is as follows: First, initialize the parameters, set the amplitude transformation weight C1 and the curve gradient weight C2, construct the cost function Cost, and input the global optimal time-frequency distribution after energy redistribution.

[0098]

[0099] Second, search The point with the largest energy in is taken as the initial point of the time-frequency ridge line IF(t), and the recursive update formula of the cost function Cost is obtained based on the initial point:

[0100]

[0101]

[0102] Among them, j is the independent variable, t m is the horizontal coordinate of the maximum time-frequency energy point, m is t m The corresponding coordinates, ω m is the ordinate of the maximum time-frequency energy point, ω(t j ) is t j The ordinate corresponding to the maximum time-frequency energy point at time T s =F s / 2 is the frequency domain range of time-frequency distribution;

[0103] Finally, according to the initial point (t m ,ω m ), the maximum frequency transformation range Δω and the recursive update formula, search for the minimum cost function point on both sides of the initial point along the time axis, and obtain the complete target ridge IF(t):

[0104] ω R =ω m ,ω L =ω m

[0105]

[0106] (5) Estimation of instantaneous phase of rotation frequency: Based on the target ridge line with minimum loss obtained in step (4), the Vold-Kalman filter is used to accurately extract the time domain waveform of the higher harmonics of the instantaneous rotation frequency and obtain the phase information of the instantaneous rotation speed. The specific process is as follows: the target ridge line with minimum loss is used as the center frequency of the Vold-Kalman filter, and the instantaneous bandwidth B of the filter is set. d , and obtain the filtered signal corresponding to the high-order harmonics; then use Hilbert transform, inverse tangent function and phase deconvolution to obtain the phase information of the instantaneous speed.

[0107] (6) Sensitive frequency band screening based on fast spectral kurtosis: First, a filter with phase quasi-analysis effect is used to perform M-level decomposition of the strongly variable speed signal in the form of double downsampling and iterative filtering, and the spectral kurtosis value of each sub-band is calculated; secondly, the spectral kurtosis of all sub-signals is compared to determine the maximum spectral kurtosis SK max The corresponding center frequency f m and the optimal bandwidth Δf m ; Finally, using the center frequency f m and the optimal bandwidth Δf m Perform bandpass filtering on the collected rolling bearing vibration signal to obtain the filtered signal x f ; The specific process is as follows:

[0108] First, a low-pass filter h1(n) and a high-pass filter h2(n) with a phase quasi-analysis effect are constructed by using FIR filters:

[0109]

[0110] Where h0(n) is the cutoff frequency;

[0111] Then, the constructed filter is used to perform M-level decomposition of the collected rolling bearing vibration signal in the form of double downsampling and iterative filtering. Each level of filtering generates 2 m sub-signals, and generate sub-signals to calculate the spectral kurtosis:

[0112]

[0113] Where S(t,f) represents the time-frequency distribution of signal s;

[0114] Finally, compare the spectral kurtosis of all sub-signals to obtain the maximum spectral kurtosis SK max The corresponding center frequency f m and the optimal bandwidth Δf m , and use f m and Δf m Perform bandpass filtering on the collected rolling bearing vibration signal to obtain the filtered signal x f .

[0115] (7) Fault feature order diagram extraction: The phase information of the instantaneous speed obtained in step (5) is used to extract the bearing fault sensitive frequency band signal, i.e., the filtered signal x, extracted in step (6). f Perform angular domain resampling and envelope demodulation analysis to obtain the fault characteristic order diagram of the bearing. The specific process is as follows: First, use the least squares fitting method to establish the mapping relationship between angle and time Secondly, determine the size of the angle domain sampling interval Δθ and find the time corresponding to the equal angle domain sampling point kΔθ Further using cubic spline interpolation, according to The filtered signal is resampled; finally, the angular domain resampled signal is envelope demodulated to obtain an order diagram representing the rolling bearing fault characteristics;

[0116]

[0117]

[0118] Among them, F s is the sampling frequency, f(t) is the minimum rotation frequency of the reference axis, and k is the coefficient.

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

[0120] In this example, rolling bearings are used as the research object, and a strong variable speed experiment is carried out using a rotating machinery fault implantation test bench. The experiment collects vibration acceleration signals through a PCB acceleration sensor, and the sampling frequency is set to 25600Hz; the speed fluctuates within the range of 1-2400r / min; for the purpose of experimental comparison, a laser tachometer is used to measure the key phase signal of the input shaft.

[0121] The vibration data of the first 7s startup process in the experimental data was intercepted for verification of the embodiment. The signal vibration waveform is as follows: Figure 2 As shown, the speed waveform is as follows Figure 3 shown.

[0122] Using the method of the present invention, the keyless phase order tracking method of rolling bearing under strong variable speed conditions is as follows:

[0123] (1) Signal low-frequency feature extraction

[0124] First, the vibration signal of the rolling bearing under the condition of strong variable speed is measured by the acceleration sensor; then, the rotation frequency f c With the frequency band as the center and 1000 Hz as the bandwidth parameter, a bandpass filter with zero phase shift characteristics is used to reduce the noise of the rolling bearing vibration signal.

[0125] (2) Local matching synchronous compression transformation

[0126] First, set the number of discrete points N a 23, the instantaneous frequency transformation range (-π / 2, π / 2) is discretized to obtain the rotation angle α and the discrete time-varying frequency modulation frequency c; secondly, the strongly variable speed signal is matched with the Chirplet transform according to the candidate time-varying frequency modulation frequency c to obtain the time-frequency distribution G(c, t, ω) under different time-varying frequency modulation frequencies c; then, the discrete number of frequency bands N is set c is 10, and the local matching kurtosis operator LK(c,t c ,Nc );Calculate t c The local matching kurtosis LK(c,t c ,N c ); Finally, according to t c The maximum local matching kurtosis LK(c,t c ,N c ), obtain the local optimal time-varying frequency modulation sequence So as to obtain the global optimal time-frequency distribution

[0127] (3) Local Maximum Synchronous Compression Transform

[0128] First, the parameter Δ is set to 20 according to the shortest distance between the instantaneous frequencies of adjacent components of the multi-component signal. Second, the local maximum synchronous compression operator ω is used. m (c, t, ω), with its local maximum as the center, redistributes the energy of the local kurtosis matching Chirplet result along the frequency direction to obtain the global optimal time-frequency distribution with high energy concentration

[0129] (4) Time-frequency ridge extraction based on cost function

[0130] First, set the amplitude transformation weight C1 to 0.2, the curve gradient weight C2 to 0.1, and construct the cost function Cost; secondly, search according to the cost function The point with the largest energy in the time-frequency ridge is taken as the initial point of the time-frequency ridge IF(t); the recursive update formula of the cost function Cost is further obtained; finally, according to the initial point (t m ,ω m ), the maximum frequency transformation range Δω and the recursive update formula, search for the minimum cost function point on both sides of the initial point along the time axis to obtain the complete target ridge line IF(t);

[0131] (5) Frequency conversion instantaneous phase estimation

[0132] The result of the ridge extraction algorithm based on the cost function is used as the center frequency of the Vold-Kalman filter. In this example, the instantaneous bandwidth B of the filter is set d =2, and the filtered signal corresponding to the high-order harmonics is obtained; then the instantaneous phase in the cumulative phase space is obtained by using Hilbert transform, inverse tangent function and phase deconvolution;

[0133] (6) Fault-sensitive frequency band screening based on fast spectral kurtosis

[0134] Firstly, the strong variable speed signal is decomposed into 10 levels by using a filter with phase quasi-analysis effect in the form of double downsampling and iterative filtering. Each level of filtering produces 2 msub-signals, and generate sub-signals to calculate the spectral kurtosis; secondly, compare the spectral kurtosis calculated by all sub-signals to obtain the maximum spectral kurtosis SK max The corresponding center frequency f m and the optimal bandwidth Δf m ; Finally, using f m and Δf m Bandpass filter the original signal to obtain the filtered signal x f ;

[0135] (7) Fault feature order diagram extraction

[0136] The phase information of the instantaneous rotation speed obtained in step (4) is used to perform angular domain resampling and envelope demodulation analysis on the bearing fault sensitive frequency band signal extracted in step (5) to obtain the bearing fault characteristic order diagram.

[0137] In order to verify the effectiveness of the present invention, the time-frequency distribution of the fourth harmonic extracted by the present invention is shown in the figure below: Figure 4 As shown in (a) and (b) in Figure 4 It can be seen in (a) that the proposed method can obtain the high-precision time-frequency distribution of the 4th harmonic; Figure 4 (b) shows that the proposed method can accurately extract low-frequency signals; the frequency multiplication estimation results extracted by the present invention are shown in Figure 5 As shown in (a) and (b) in Figure 5 As can be seen in (a), the rotation frequency extracted by the proposed method is very close to the actual rotation frequency, with an absolute error within 1 Hz; Figure 5 (b) shows that the relative estimation error of the proposed method for the rotation frequency is less than 2.7%, which can well meet the actual needs of engineering. Figure 6 As shown in (a) and (b) in Figure 6 It can be seen from (a) that the rotation frequency estimated by the synchronous compression transformation method has a large gap with the actual rotation frequency within the initial 0.5s, especially in the range of 0-0.1s, the absolute error is obviously too large; Figure 6 The method shown in (b) has an error of more than 20% in the initial 1s rotation frequency estimation, which is large and affects the effect of subsequent order tracking. Compared with the synchronous compression transform, the local matching synchronous compression transform used in the present invention has a higher rotation speed estimation accuracy and can be applied to the startup phase. The endpoint effect and error are smaller (less than 2.7%). The index results of the fault sensitive frequency band selection method based on fast spectral kurtosis involved in the present invention are as follows: Figure 7 As shown, at level 5, K max =16.4151, optimal window length = 64, center frequency = 2.6kHz, bandwidth = 0.4kHz, Figure 8 The bondless phase order tracking results obtained by the method of the present invention are explained. Fig. 9 This is the result of keyless phase order tracking using the instantaneous rotation frequency extracted by synchronous compression transform; Figure 8 and Fig. 9 It can be seen that the local kurtosis matching Chirplet transform and local maximum synchronous compression transform methods proposed in the present invention can obtain a more accurate fault order (the fault order is 3.94), with obvious fault characteristics and clear fundamental frequency, multiple frequency and sidebands.

[0138] Based on the bearing outer ring fault data under the starting condition of the rotating machinery fault implantation test bench, the present invention realizes the accurate extraction of high-precision time-frequency distribution and instantaneous rotation frequency by using the local matching synchronous compression transformation method and the ridge extraction method based on the cost function; the Vold-Kalman filter is used to accurately estimate the instantaneous phase of the shaft; the signal is adaptively filtered based on the fast spectral kurtosis theory to select the fault sensitive frequency band; the fault sensitive frequency band signal is resampled in the angular domain and envelope demodulated and analyzed according to the phase information of the instantaneous rotation speed, and the fault characteristic order diagram of the bearing under the strongly variable speed condition is obtained. Compared with other traditional methods, the present invention can overcome the problems of low aggregation of time-frequency distribution and fuzzy instantaneous rotation frequency caused by the strong frequency modulation characteristics of the strongly variable speed signal, realize the acquisition of high-precision time-frequency distribution, accurately extract the instantaneous shaft frequency ridge under the starting condition, and realize the accurate estimation of the fundamental frequency and multiple frequency of the rolling bearing fault characteristics under the strongly variable speed condition.

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

Claims

1. A keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed, characterized in that: The steps include: (1) De-noise the vibration signal of the rolling bearing under strong variable speed conditions to obtain the low-frequency characteristics of the signal; (2) The optimal time-varying frequency modulation frequency of each frequency band of the signal is adaptively determined by maximizing the local matching kurtosis operator based on the low-frequency characteristics of the signal, and then the high-resolution time-frequency distribution of the multi-component signal is obtained by matching Chirplet transform; (3) Taking the local maximum value in the frequency direction as the center, the energy is redistributed using the local maximum synchronous compression transform operator; (4) Using a ridge extraction algorithm based on a cost function, the target ridge with the minimum loss is obtained according to the global optimal time-frequency distribution after energy redistribution; (5) according to the target ridge line with minimum loss obtained in step (4), extract the time domain waveform of the higher harmonics of the instantaneous rotation frequency to obtain the phase information of the instantaneous rotation speed; (6) bandpass filtering the vibration signal of the rolling bearing under the strong variable speed condition using the center frequency and the optimal bandwidth to obtain a filtered signal; (7) using the phase information of the instantaneous speed obtained in step (5) to perform angular domain resampling and envelope demodulation analysis on the filtered signal in step (6) to obtain a fault characteristic order diagram of the bearing; The specific process of step (2) is as follows: According to the variable frequency modulation frequency c, the low-frequency characteristic signal of the signal obtained in step (1) is matched with Chirplet transformation to obtain the time-frequency distribution G(c, t, ω) under different time-varying frequency modulation frequencies c: Where t′ is time, ω is frequency, g(ut′) is Gaussian window function, s(u) is discrete signal, i is imaginary sign, u is independent variable, -c(ut′) 2 / 2 is the demodulation operator; Then, set the discrete number of frequency bands to N c , along the frequency direction, the signal at t c The spectrum range at the moment is divided into N c equal-spaced frequency bands, according to the local matching kurtosis operator LK(c,t c ,N c ), calculate the signal t c The local matching kurtosis LK(c,t c ,N c ); Finally, by maximizing the signal at t c The local matching kurtosis operator LK(c,t c ,N c ), and obtain t c The local time-varying frequency sequence of Get the global optimal time-frequency distribution 2. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1 is characterized in that: Signal c Time i n The local matching kurtosis of a frequency band generated by different time-varying FM frequencies is calculated by the following formula: In the formula, t c The signal at time i n The optimal time-varying FM frequency for each frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the first frequency band, t c The optimal time-varying FM frequency of the instantaneous signal in the second frequency band, t c The time signal is at the Nth c The optimal time-varying FM frequency for the frequency band.

3. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1 is characterized in that: The specific process of step (3) is as follows: First, for the high-resolution time-frequency distribution of the multi-component signal obtained in step (2), the local maximum value in the frequency direction is selected as the center, and the spectrum is divided into several frequency bands according to the bandwidth that meets the separation conditions of different instantaneous frequency signals: Then, using the local maximum synchronization compression operator ω m (c, t, ω), with the local maximum as the center, the high-resolution time-frequency distribution of the multi-component signal obtained in step (2) is redistributed along the frequency direction to obtain the global optimal time-frequency distribution after energy redistribution.

4. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 3 is characterized in that: The frequency bandwidth that satisfies the separation conditions of different instantaneous frequency signals satisfies the following formula: Among them, Δ is the frequency bandwidth that satisfies the separation conditions of different instantaneous frequency signals, is the instantaneous frequency of the k+1th component, is the instantaneous frequency of the kth component; The global optimal time-frequency distribution after energy redistribution is as follows: Among them, ω m (c, t, ω) is the local maximum synchrocompression operator, δ is the delta function, and η represents the frequency domain coordinate after energy redistribution.

5. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1 is characterized in that: Step (4) Specific process: The point with the largest energy in the global optimal time-frequency distribution after energy redistribution is searched as the initial point of the time-frequency ridge, and the recursive update formula of the cost function is obtained according to the initial point; Finally, according to the initial point, the maximum frequency transformation range and the recursive update formula, the minimum cost function point is searched on both sides of the initial point along the time axis to obtain the target ridge line.

6. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 5 is characterized in that: The recursive update formula of the cost function is as follows: IF(t m )=ω m Among them, j is the independent variable, t m is the horizontal coordinate of the maximum time-frequency energy point, m is t m The corresponding coordinates, ω m is the ordinate of the maximum time-frequency energy point, ω(t j ) is t j The ordinate corresponding to the maximum time-frequency energy point at time T s =F s / 2 is the frequency domain range of the time-frequency distribution; C1 is the amplitude transformation weight, C2 is the curve gradient weight, Cost is the cost function, is the global optimal time-frequency distribution after energy redistribution, η represents the frequency domain coordinate after energy redistribution, is the global optimal time-frequency distribution after energy redistribution.

7. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1 is characterized in that: The specific process of step (5) is as follows: the target ridge line with the minimum loss is used as the center frequency of the Vold-Kalman filter, the instantaneous bandwidth of the filter is set, and the filtered signal corresponding to the high-order harmonics is obtained; then the Hilbert transform, inverse tangent function and phase deconvolution are used to obtain the phase information of the instantaneous rotation speed.

8. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1 is characterized in that: The center frequency and optimal bandwidth in step (6) are determined by the following process: The strong variable speed signal is decomposed by using a filter with phase quasi-analysis effect in the form of double downsampling and iterative filtering, and the spectral kurtosis SK(ω) of each sub-band is calculated; Among them, S(t,ω) represents the time-frequency distribution of signal s; The spectral kurtosis of all sub-signals is compared to determine the center frequency and optimal bandwidth corresponding to the maximum spectral kurtosis.

9. The keyless phase order tracking method for rolling bearing fault diagnosis under strong variable speed according to claim 1, characterized in that: The specific process of step (7) is as follows: First, the least squares fitting method is used to establish the mapping relationship between angle and time. Secondly, determine the size of the angle domain sampling interval Δθ and find the time corresponding to the equal angle domain sampling point pΔθ Using the cubic spline interpolation method, according to the time corresponding to the sampling point pΔθ in the equiangular domain resampling the filtered signal; Finally, the envelope demodulation of the angular domain resampled signal machine is performed to obtain the order diagram that characterizes the rolling bearing fault characteristics; Among them, the time corresponding to the angle domain sampling interval and the equal angle domain sampling point pΔθ is as follows: In the formula, F s is the sampling frequency, f(t) is the minimum rotation frequency of the reference axis, and p is the coefficient.

Images