A multi-dimensional data compression method for rolling bearing fault diagnosis

Abstract

The application discloses a kind of multidimensional data compression methods for rolling bearing fault diagnosis, it is compressed to the sample length by carrying out short-time Fourier transformation to signal first, while extracting the time-frequency information of signal and storing in short-time Fourier coefficient matrix;Subsequently, the time-frequency information is modeled using hidden Markov model, and the hidden two state information is mined, and the state sequence obtained can be expressed as binary sequence, so that the bit number dimension compression is realized;Finally, the wavelet decomposition of state sequence is used for fault frequency resolution positioning, and the signal is down-sampled based on the corresponding fault frequency resolution, the process of approximate wavelet compression is realized, and the data in component dimension is compressed.Said above, by multidimensional compression, while realizing super-high compression ratio, the fault characteristics in original data are effectively retained.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of fault diagnosis, and more particularly relates to a multi-dimensional data compression method for rolling bearing fault diagnosis. BACKGROUND

[0002] In the process of rolling bearing fault diagnosis, vibration monitoring data is the most commonly used data type, because based on vibration data not only good diagnosis effect can be achieved, but also clear machine interpretation can be obtained. By collecting high-frequency vibration data, the fault characteristics are effectively retained, but at the same time, the big data problem is also brought. With the increase of state monitoring points, this problem is becoming more and more prominent.

[0003] At the same time, real industrial equipment often experiences a long running time and a slow degradation process, which inevitably aggravates the big data problem. Massive data not only increases the pressure of data storage and processing, but also consumes a lot of energy in the data transmission process. Therefore, many data compression (DC) methods are used to solve this problem.

[0004] The most core parameter in the DC method is the compression ratio (CR), which is defined as the ratio of the bit number of the original data to the compressed data. Undoubtedly, the goal of the DC method is to achieve as high CR as possible while ensuring effective information retention. Although many methods have been developed to achieve DC under specific conditions, such as using coding techniques in the field of speech compression to achieve compression, using the sparsity of signals to filter out invalid components, and other mechanism-based methods. However, in order not to affect the diagnosis effect of the data, it is still very difficult to achieve high CR. SUMMARY

[0005] The purpose of the present application is to overcome the shortcomings of the prior art and provide a multi-dimensional data compression method for rolling bearing fault diagnosis. This method introduces bit number dimension compression while avoiding mutual constraints between different dimension compressions. The constructed multi-dimensional data compression method retains fault characteristic information while achieving ultra-high compression ratio, thereby overcoming the bottleneck problem of further improving CR in current fault diagnosis DC methods.

[0006] To achieve the above-mentioned application purposes, the present application provides a multi-dimensional data compression method for rolling bearing fault diagnosis, characterized in that it comprises the following steps:

[0007] (1) time-frequency transformation of vibration data;

[0008] (1.1) collecting vibration data X={x1, x2, …, x t ,…,xn},x t represents the vibration data value of the t-th sampling point collected, and n is the total number of sampling points;

[0009] (1.2), performing short-time Fourier transform on the vibration data X, so as to convert the time-domain vibration data into time-frequency domain data;

[0010] (1.3), storing the short-time Fourier transform coefficient after time-frequency conversion in a matrix X stft , wherein the element x stft (i,f) in the i-th column and the f-th row in the matrix X stft is calculated by the following formula:

[0011]

[0012] wherein w[] represents an analysis window function; x[] refers to a sequence divided by the window function, u represents the serial number of the sample in the sequence; N w represents the window width; R is the time shift amount of the window; i represents the serial number of the window, i = 1, 2, …, L, L = floor((n+R-N w ) / R); f represents the serial number of the frequency in the spectrum distribution, f = 1, 2, …, N f , N f = N w +1, and floor() is a floor function;

[0013] (2), time-frequency domain state modeling;

[0014] (2.1), defining a Hidden Markov Model (HMM);

[0015] The HMM is represented by a parameter set λ = {π, A, C}, wherein π is an initial probability vector, A and C are respectively a transition probability matrix and a covariance matrix of an observation distribution;

[0016] (2.2), initialization of parameters of the HMM;

[0017] The parameters of the HMM are initialized by a sampling random method, which needs to satisfy the following constraints:

[0018]

[0019] wherein π I represents the probability of the initial state being I, and a IJ represents the transition probability of the state I to the state J;

[0020] (2.3), taking the column vectors in the matrix X stft as the observation values of the HMM, so that the matrix X stftRe-expressing as observation sequence O = {o i}, i = 1, 2, …, L, o i represents the i-th observation value, and the state sequence corresponding to the observation sequence O is represented as S = {s i}, s i represents the state corresponding to the i-th observation value, s i ∈ {0, 1};

[0021] (2.4), using the expectation-maximization algorithm to estimate the parameters of the HMM;

[0022] (2.4.1), setting the maximum number of iterations Z, initializing the current iteration number z = 1, z = 1, 2, …, Z; setting the threshold value ε;

[0023] (2.4.2), defining the forward probability and the backward probability of the i-th observation value as α i (J) and β i (I) respectively:

[0024]

[0025] where s i (J) represents the state value of the i-th observation value in state J;

[0026] (2.4.3), at the z-th iteration, the forward and backward probabilities of the i-th observation value are solved;

[0027]

[0028] where α i-1 (I) represents the forward probability of the i-1-th observation value;

[0029] (2.4.4), similarly, by traversing each observation value o i , the forward and backward probability values of all observation values at the z-th iteration are calculated according to formula (4);

[0030] (2.4.5), calculating the state occupancy probability and the state transition probability;

[0031] Given the parameter λ and the observation value o i , the occupancy probability γ i (I) of state I of the i-th observation value and the state transition probability ξ i (I, J) are calculated;

[0032]

[0033] (2.4.6), estimating the model parameters;

[0034]

[0035] wherein: represents o i , the conjugate transpose of is an estimate of the model parameters;

[0036] (2.4.7), calculating the likelihood probability value p z (O) of the model after the zth iteration:

[0037] p z (O) = a L (0) + a L (1) (7)

[0038] wherein a L (0) represents the forward probability of the Lth observation being in state 0, and a L (1) represents the forward probability of the Lth observation being in state 1;

[0039] (2.4.8), if the difference between the probability values of the parameter estimates after two adjacent iterations is less than a specified threshold value e or z > Z, then go to step (3); otherwise, update the elements of the model parameters in the model parameter set l using the model parameters estimated after this round of iteration, then let z = z + 1, and return to step (2.4.3);

[0040] (3), state recognition and conversion;

[0041] (3.1), defining a path optimization mathematical model;

[0042]

[0043] (3.2), using the Viterbi algorithm to iteratively solve the path optimization mathematical model to obtain the joint probability d i (I) of the path with the maximum transition probability of the state I of the ith observation, and record the previous state x i (I) of the optimal path of the current ith observation;

[0044]

[0045] wherein d i (I) represents the joint probability of the path with the maximum transition probability of the state I of the ith observation; p(o i | s i = I) represents the probability of the observation o i when the state s i is in state I;

[0046] (3.3), obtaining the final state;

[0047]

[0048] (3.4) calculating the optimal path;

[0049] The states are backtracked in turn:

[0050]

[0051] After backtracking, the states are rearranged in chronological order to obtain the final state sequence s * = {s i * , i = 1, 2, …, L};

[0052] (4) downsampling based on component reduction;

[0053] (4.1) wavelet decomposition of the state sequence s * ;

[0054]

[0055] wherein, is the wavelet coefficient obtained by wavelet orthogonal transformation, ψ m,k (t) is the scale transformation and translation of the analysis wavelet, and m and k represent scale and translation parameters;

[0056] (4.2) after wavelet decomposition of the state sequence using different m and k parameter values, small values in the corresponding wavelet coefficients are set to zero by thresholding, and then the state sequence s * is restored, and after restoration, spectral analysis is performed to determine the optimal scale parameter m * ;

[0057] (4.3) equidistant downsampling of the state sequence s * at a downsampling rate of to obtain the final compressed data.

[0058] The invention aims to achieve the following:

[0059] The application is a multi-dimensional data compression method for rolling bearing fault diagnosis, which firstly realizes sample length compression by short-time Fourier transformation of the signal, extracts time-frequency information of the signal and stores it in a short-time Fourier coefficient matrix; then, hidden state information is mined by modeling the time-frequency information with a hidden Markov model, and the obtained state sequence can be expressed as a binary sequence, thereby realizing bit number dimension compression; finally, wavelet decomposition of the state sequence is used for fault frequency resolution positioning, and the signal is down-sampled based on the corresponding fault frequency resolution, which is similar to the wavelet compression process, thereby realizing component dimension compression. In summary, through multi-dimensional compression, ultra-high compression ratio is realized while effectively retaining the fault features in the original data.

[0060] Meanwhile, the multi-dimensional data compression method for rolling bearing fault diagnosis has the following beneficial effects:

[0061] The application is a data compression method for rolling bearing fault diagnosis, which aims to further improve the compression ratio in current research, expands the data compression range by introducing bit number dimension compression in mechanism-based methods, improves the bottleneck problem that the compression of this type of method is only limited to sample length and component dimensions, and avoids the mutual restriction problem in different dimension data compression by down-sampling to approximate the component compression process. BRIEF DESCRIPTION OF DRAWINGS

[0062] Figure 1 is a flowchart of the multi-dimensional data compression method for rolling bearing fault diagnosis of the application;

[0063] Figure 2 is a test bench for collecting the Case Western Reserve rolling bearing data set;

[0064] Figure 3 is a time-frequency domain graph of the inner ring fault bearing before and after data compression; DETAILED DESCRIPTION

[0065] The specific embodiments of the application are described below with reference to the accompanying drawings, so that those skilled in the art can better understand the application. It should be particularly noted that in the following description, when the detailed description of known functions and designs may obscure the main content of the application, these descriptions will be omitted here.

[0066] EMBODIMENT

[0067] Figure 1 is a flowchart of the multi-dimensional data compression method for rolling bearing fault diagnosis of the application.

[0068] In this embodiment, as Figure 1As shown, the application is a multi-dimensional data compression method for rolling bearing fault diagnosis, and a specific embodiment is used to test and verify the application.

[0069] S1, time-frequency transformation of vibration data;

[0070] S1.1, the embodiment uses the rolling bearing data set published by Case Western Reserve University for testing, and the data collection is carried out on the test bench at Figure 2 . The test bearing is driven by a 2-horsepower motor, and the load is 0 to 3 horsepower. Two sampling frequencies are used to collect vibration data: 12 kHz and 48 kHz. In order to approach the big data problem we are interested in, we select 48 kHz data for the experiment. The experimental data of the inner race bearing fault with a speed of 1750 rpm and a defect size of 0.1778 mm is represented as X = {x1, x2, …, x t ,…,x 485643}, x t represents the vibration data value of the tthsample point collected;

[0071] S1.2, short-time Fourier transform is performed on the vibration data X, the selected analysis window function w[] is a periodic Hann window, the window width N w is 256, and the window time shift R is 32, so that the time domain vibration data is transformed into time-frequency domain data;

[0072] S1.3, the short-time Fourier transform coefficient after time-frequency transformation is stored in matrix X stft , wherein the element x stft (i,f) in the ithcolumn and the fthrow of matrix X stft is calculated as follows:

[0073]

[0074] Here, x[] refers to the sequence divided by the window function, u represents the sequence number of the sample in the sequence; i represents the sequence number of the window, i = 1, 2, …, 15169; f represents the sequence number of the frequency in the spectrum distribution, f = 1, 2, …, 257;

[0075] In this way, the length dimension of the vibration data is reduced through time-frequency transformation, i.e. the sample length is reduced from 485643 to 15169. Although additional spectral information is introduced, these spectral information will be used in the subsequent time-frequency domain state modeling and state transformation process, and then they will be removed.

[0076] S2, time-frequency domain state modeling;

[0077] S2.1, according to the pulse characteristics of the signal, the time-varying data is divided into two states: random steady state and repeated transient state. The HMM is used to model the potential state sequence, and the HMM is defined first, and the parameter set λ={π, A, C} is used to represent, wherein π is the initial probability vector, A and C are the transition probability matrix and the covariance matrix of the observation distribution respectively;

[0078] S2.2, parameter initialization of HMM;

[0079] The sampling random method is used to initialize the parameters of the HMM, which needs to meet the following constraints:

[0080]

[0081] Wherein, π I represents the probability of the initial state being I, a IJ represents the transition probability of state I to state J;

[0082] S2.3, the column vector in the matrix X stft is taken as the observation value of the HMM, so that the matrix X stft is re-expressed as the observation sequence O={o i}, i=1, 2, …, 15169, o i represents the ith observation value, and the state sequence corresponding to the observation sequence O is represented as S={s i}, s i represents the state corresponding to the ith observation value, s i ∈{0, 1};

[0083] S2.4, the expectation maximum algorithm is used to estimate the parameters of the HMM;

[0084] S2.4.1, set the maximum iteration number Z=50, initialize the current iteration number z=1, z=1, 2, …, 50; set the threshold value ε=e 0.00001 ;

[0085] S2.4.2, define the forward probability and the backward probability of the ith observation value as α i (J) and β i (I) respectively:

[0086]

[0087] Wherein, s i (J) represents the state value of the ith observation value in state J;

[0088] S2.4.3, at the zth iteration, the forward and backward probabilities of the ith observation value are solved;

[0089]

[0090] where α i-1 (I) denotes the forward probability of the i-1th observation;

[0091] S2.4.4. Similarly, the forward and backward probability values of all observations at the zth iteration are calculated according to formula (4) by traversing each observation o i

[0092] S2.4.5. Calculate the state occupancy probability and state transition probability;

[0093] Given the parameter λ and the observation o i , the occupancy probability γ i (I) of the state I of the i-th observation and the state transition probability ξ i (I, J) are calculated.

[0094]

[0095] S2.4.6. Estimate the model parameters;

[0096]

[0097] where: denotes the conjugate transpose of o i , and is the estimated value of the model parameter;

[0098] S2.4.7. Calculate the likelihood probability value p z (O) of the model after the zth iteration:

[0099] p z (O) = α 15169 (0) + α 15169 (1) (7)

[0100] where α 15169 (0) denotes the forward probability of the 15169th observation being in state 0, and α 15169 (1) denotes the forward probability of the 15169th observation being in state 1.

[0101] S2.4.8. If the difference between the probability values of the parameter estimates after two adjacent iterations is less than a specified threshold e 0.00001 or z ≥ 50, go to step S3; otherwise, use the model parameters estimated after this round of iteration to update the elements of the parameters in the model parameter set λ, and then let z = z + 1, and return to step S2.4.3.

[0102] S3. State recognition and conversion;

[0103] S3.1. Define a path optimization mathematical model;​

[0104]

[0105] S3.2. Use the Viterbi algorithm to iteratively solve the mathematical model of path optimization, and obtain the joint probability δ of the path with the highest transition probability of the i-th observation being in state I. i (I) Then record the previous state χ of the optimal path for the current i-th observation. i (I);

[0106]

[0107] Where, δ i (I) represents the joint probability of the path with the highest transition probability when the state of the i-th observation is I; p(o i |s i =I) represents state s i The observed value is o when in state I. i The probability of;

[0108] S3.3 Obtain the final state;

[0109]

[0110] S3.4 Calculate the optimal path;

[0111] Backtrack through the states sequentially:

[0112]

[0113] After backtracking is complete, the states are rearranged in chronological order to obtain the final state sequence s. * ={s i * ,i=1,2,…,15169};Since s * Since there are only two states, binary data can be stored in a single logical variable, thereby reducing the number of bits per sample and significantly decreasing the data size.

[0114] S4. Downsampling based on component reduction;

[0115] S4.1, For the state sequence s * Perform wavelet decomposition and select the analytic wavelet ψ(t) as the Halkyne function;

[0116]

[0117] in, ψ represents the wavelet coefficients obtained from the wavelet orthogonal transform. m,k(t) is the scale transform and translation of the analysis wavelet, m, k represent scale and translation parameters;

[0118] S4.2, after the state sequence is decomposed using different m, k parameter values, small values in the corresponding wavelet coefficients are set to zero using threshold method, and then the state sequence s * is restored, and after restoration, frequency spectrum analysis is performed to obtain an optimal scale parameter m* = 3; since the wavelet compressed data loses the advantage of binary state sequence, down-sampling is used to approximate the wavelet compression process to avoid the mutual restriction problem between the single sample bit dimension and the component dimension in the data compression process.

[0119] S4.3, according to the optimal scale parameter, the state sequence s * is subjected to equidistant down-sampling at a down-sampling rate of 4 times to obtain final compressed data with a length of 3793.

[0120] S5, the diagnosis effect performance of the compressed data is verified using the fault diagnosis result;

[0121] S5.1, a rolling bearing fault characteristic frequency: ball pass frequency on inner race (BPFI).

[0122]

[0123] wherein f r is the rotation frequency of the shaft, d and n are the diameter and number of the rolling elements, and D and are the pitch diameter of the bearing and the contact angle of the rolling elements.

[0124] S5.2, time-frequency conversion is performed on the data before and after compression using fast Fourier transform, and the results are shown in Figure 3 . Wherein, a) and b) are the original vibration data and its spectrum, and c) and d) are the compressed data and its spectrum. From the time domain, the data length is obviously compressed, and the binary of the vibration data is realized, which reflects that the size of the data is effectively compressed. The compression ratio of the compressed experimental data is 8194, and the data compression effect is much better than the method in the existing literature.

[0125] S5.3, based on the BPFI = 157.94 Hz calculated in S5.1, the peak value at the corresponding position in the spectrum is searched to perform fault diagnosis. From the spectrum obtained in S5.2, the data before and after compression both have obvious peak values at BPFI, the inner ring fault of the rolling bearing can be accurately identified from the compressed data, and the fault characteristic frequency is more prominent in the spectrum of the compressed data.

[0126] In summary, the experimental results verify the effectiveness of the proposed multi-dimensional data compression method, which can greatly improve the compression effect on the basis of preserving fault feature information.

[0127] Although the above describes the illustrative embodiments of the present application in order to assist those skilled in the art to understand the present application, it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, as long as various changes are within the spirit and scope of the present application defined and determined by the appended claims, all the inventions utilizing the concept of the present application are within the protection of the present application.

Claims

1. A method for multi-dimensional data compression for bearing fault diagnosis, characterized in that, The method comprises the following steps: (1) time-frequency transformation of vibration data; (1.1) Collect vibration data of the bearing under test X={x1,x2,…,x t ,…,x n }, x t This represents the vibration data value collected at the t-th sampling point, where n is the total number of sampling points; (1.2) performing short-time Fourier transform on the vibration data X, so as to transform the time-domain vibration data into time-frequency domain data; (1.3) storing the short-time Fourier transform coefficients after time-frequency conversion in a matrix X stft wherein the element x stft (i,f) of the i-th column and f-th row of the matrix X stft is calculated by wherein w[] represents an analysis window function; x[] refers to a sequence divided by the window function, u represents the serial number of samples in the sequence; N w represents the window width; R is the time shift amount of the window; i represents the serial number of the window, i = 1, 2, …, L, L = floor((n + R - N w ) / R); f represents the serial number of the frequency in the frequency spectrum distribution, f = 1, 2, …, N f , N f = N w +1, and floor() is a floor function. (2) state modeling in time-frequency domain; (2.1) defining a Hidden Markov Model (HMM); The HMM is represented by a parameter set λ = {π, A, C}, wherein π is an initial probability vector, A and C are respectively a transition probability matrix and a covariance matrix of an observation distribution; (2.2) initialization of parameters of the HMM; The parameters of the HMM are initialized by a sampling random method, and the following constraints need to be met: wherein, π I represents the probability of the initial state being I, a IJ represents the transition probability from state I to state J; (2.3), the column vectors in matrix X stft are taken as the observation values of HMM, so that matrix X stft is re-expressed as an observation sequence O = {o i}, i = 1, 2, …, L, o i represents the i-th observation value, and the state sequence corresponding to the observation sequence O is represented as S = {s i}, s i represents the state corresponding to the i-th observation value, s i ∈ {0, 1}. (2.4) parameter estimation of the HMM by using an expectation maximization algorithm; (2.4.1) setting a maximum iteration number Z, initializing a current iteration number z = 1, z = 1, 2, …, Z, and setting a threshold value ε; (2.4.2), the forward and backward probabilities of defining the ith observation are α i (J) and β i (I): where s i (J) denotes the state value of the ith observation being in state J. (2.4.3) in the zth iteration, solving the forward and backward probabilities of the ith observation value; wherein a i-1 (I) denotes the forward probability of the i-1th observation; (2.4.4) Similarly, the forward and backward probabilities for all observations are computed at the zth iteration according to equation (4) by iterating through each observation o i , the forward and backward probabilities for all observations are computed at the zth iteration according to equation (4) by iterating through each observation o (2.4.5) calculating state occupation probability and state transition probability; The occupation probability γ i (I) of state I for the i-th observation and the state transition probability ξ i (I, J) are calculated when the parameter λ and the observation o i (I, J) are given. (2.4.6) estimating model parameters; wherein: represents o i the conjugate transpose of o, is an estimate of the model parameter; (2.4.7) calculating the likelihood probability value p of the model after the zth iteration z (O): p z (O) = a L (0) + a L (1) (7) where α L (0) denotes the forward probability that the Lth observation is in state 0, and α L (1) denotes the forward probability that the Lth observation is in state 1. (2.4.8) if the difference between the probability values of the parameter estimates after two adjacent iterations is less than the specified threshold value ε or z ≥ Z, going to step (3); otherwise, updating the elements of the parameters in the model parameter set λ by using the model parameters estimated after the current iteration, then letting z = z + 1, and returning to step (2.4.3); (3) state recognition and transformation; (3.1) defining a path optimization mathematical model; (3.2), the path optimization mathematical model is solved by using the Viterbi algorithm to obtain the joint probability δ of the path with the maximum transition probability of the state I of the i-th observation value i (I), and the last state χ of the optimal path of the current i-th observation value is recorded i (I). where δ i (I) denotes the joint probability of the path with the largest transition probability of state I for the i-th observation; p(o i |s i ) denotes the probability of the observation o i given that the state s i is in state I. (3.3) obtaining a final state; (3.4) calculating an optimal path; backtracking the states in sequence: After the backtracking is completed, the final state sequence s is obtained by rearranging the state sequence s in chronological order * = {s i * , i = 1, 2, …, L} (4) down-sampling based on component reduction; (4.1) for a state sequence s * wavelet decomposition; wherein ψ are wavelet coefficients obtained from a wavelet orthogonal transform m,k (t) are scale and translation parameters of the analysis wavelet (4.2), after wavelet decomposition of the state sequence with different m, k parameter values, the small values in the corresponding wavelet coefficients are set to zero by threshold method, then the state sequence s * is restored, and after restoration, frequency spectrum analysis is performed to determine the optimal scale parameter m * ; (4.3), with downsampling the state sequence s * at a constant rate to obtain the final compressed data.

Images