Abstract

Effective fault detection and diagnosis (FDD) for rotating machinery is always a focus issue in improving the prognostic and health management (PHM) of the equipments. The existing usage of similarity measurement has been widely spread in searching the homologous fault responses from vibration signals, but most of them are just suitable for stable speed and cannot be applied in all variable speed conditions. In order to improve measurement performance, a fast-meshed phase portrait (FMPP) frame combining the phase-space technique and box-scoring calculation is proposed. Firstly, the variable-speed signal is divided into multiple undetermined fragments according to fault characteristic orders (FCOs). Secondly, the undetermined fragments are reconstructed into corresponding phase-space trajectories to overcome the time-delay matching inconsistency of the variable speed fragments. Thirdly, the phase-space trajectories are mapped into meshed phase portraits via box-scoring calculation. Such decisive calculation can effectively transforms the diverse unequal fragments into the phase diagrams with same size, which saves time for the subsequent similarity measurement. Finally, the proposed FMPP is tested both for accuracy and timeliness on a self-built bearing bench, where the order tracking (OT) and dynamic time warping (DTW) methods are used for the comparisons. The experiments proved the effectiveness of the proposed method.

1. Introduction

Time-domain waveform clearly and intuitively reflects the inherent characteristics of machine vibration, which is a key technique widely used in mining nonobvious fault information from vibration signals. In addition, the characteristic waveform (CW), as the small fragments cut from time-domain waveform, is intuitive, clear physical representation and comprehensive reflection of the inherent characteristics [1] when compared with the characteristic statistics and the spectral analysis method. By clustering a large number of CW, the possible fault types can be mined so as to provide samples for the establishment of fault diagnosis and identification system, which applies machine learning or deep neural network and other methods [25]. Since each CW is a sequence, the clustering problem of CW is naturally transformed into the similarity calculation of sequence.

The similarity measurement of time series is always the focus and difficulty. The commonly used measurement methods of time series similarity are cosine similarity, symbolic aggregate approximation (SAX), Euclidean distance, distance with real penalty, etc., [6]. In 2006, Dong et al. [7] used the cosine similarity to measure the matching affinity between two vectors. The results showed that the algorithm based on the cosine similarity has a good application prospect in online signal monitoring. Georgoulas et al. [8] introduced a novel bearing fault detection method that uses the symbolic aggregate approximation (SAX) framework and associated related intelligent to represent the diagnostic representation of the extracted vibration record. The fault diagnosis is considered as a classification problem, the vibration signal and its feature vector are classified, and the classification accuracy is significantly improved. Lines and Anthony [9] had compared the classification accuracy of nine distance similarity measures on 38 data sets from different scientific fields. One of the main conclusions is that although the latest similarity measure is theoretically attractive, in most cases, its effectiveness is lower than the widely used measure. Among them, these methods have been well applied in fault diagnosis under stationary conditions.

However, in engineering applications, mechanical equipment operating in nonstationary conditions (especially variable speed) accounts a large proportion. There are two main difficulties in the similarity measurement of the rotational mechanical vibration signal under variable speed conditions. On the one hand, the length of the time series is inconsistent. On the other hand, the energy of the vibration signal is different, which is reflected in the difference in signal amplitude due to the SAX and cosine similarity methods cannot measure time series similarities of different lengths. Consequently, they are limited under nonstationary conditions. At present, the common methods to overcome the influence of rotational speed change are the OT [1012] and the similarity calculation method of DTW. Firstly, the essence of OT is to resample the original vibration signal at a constant angle increment and transform the nonstationary signal in the time domain into the stationary signal in the angle domain [13, 14]. In 2006, Saavedra and Rodriguez [15] carried out detailed analysis and verification of the algorithm for calculating order ratio analysis. The vibration signal will be stretched or compressed (deformed) when it is sampled in the angular domain by the OT method [16]. The deformation of the vibration signal seriously affects the accuracy of its similarity estimation. Secondly, the DTW method can calculate the similarity of unequal long time series. In 2017, Han et al. [1] proposed a rolling bearing fault identification method based on multiscale dynamic time warping (MDTW) which is applied to measure the similarity of the normalized time-domain characteristic waveform, and the results are used to classify different types of bearing faults. Different types of bearing faults are classified according to the calculation results. In 2013, Zhen et al. [17] analyzed the stator current signal of the motor driver fault using the DTW method. In 2019, Entezami and Shariatmadar [18] proposed a correlation-based DTW method to detect damage by using randomly high-dimensional multivariate features. The experimental results suggest that this method can detect and locate the damage under the nonstationary signals. However, the high time complexity of DTW is difficult to meet the needs of high-frequency computing for the era of big data.

For this case, it is necessary to find a low time complexity way to calculating vibration signals’ similarity, which not only overcomes the influence of variable speed but also maintains the essential characteristics of fault response. The phase space [19] can maintain the geometric invariance of the intrinsic structure of the source dynamical system corresponding to the time-domain signal. Hence, the phase in space can be used to represent the intrinsic link of each response that the fault state excited in the signal. And box-scoring calculation [20] can map time series of different lengths into equal phase portraits to overcome the influence of variable speed.

In view of above analysis, a new method based on phase-space and box-scoring calculation is proposed, which is called the fast-meshed phase portrait (FMPP) method. The method has four steps: (a) acquiring the signal fragments of fault response, (b) demodulating the signal fragments (the resonance sparse decomposition method used in the paper), (c) mapping the signal fragments into equal-sized mesh phase portraits, and (d) calculating the cosine similarity of the mesh phase portraits of the signal fragments.

2. Fault Signal Fragments’ Similarity Calculation Method

This section describes the four steps which are used to calculate the fault signal fragments’ similarity. The frame is shown in Figure 1.

2.1. Acquiring the Signal Fragments of Fault Response at Equal Angles

A fault response fragment in the vibration signal is the collision of the damaged component and other component. In other papers, fault response fragments are called CW. The fault response fragment which contains a complete impulse response reflects the health of the component. The fault response fragment which contains an impulse response reflects the complete health information of the component. The following describes how to get the response fragments.

In the acquired time-domain vibration data at variable speeds, the period of the fault response is not stable. Due to the different geometrical characteristics of mechanical components, different types of faults have different fault characteristic orders (FCOs) such as each bearing running one circle, and the number of times the rolling element passes the outer ring defect position is the fault characteristic order of outer ring FCOO. When the bearing running 1/FCOO circle, a fault response will be triggered, and the vibration data collected during this process are fault response fragments. Hence, we can obtain the fault response fragments according to cutting the signal at equal angle based on FCO in the angular domain.

First, the angle change curve is calculated according to the following equation based on the integral of the rotational speed at time t:

For discrete signals, equation (1) can be rewritten aswhere fs is the sampling frequency and . The vibration signal is then sliced at equal angle ∆θ, and ∆θ is an interval angle of two fault responses. Then, the ∆θ calculation method is as follows:

The vibration signal is represented by . When the length of a single response fragment is n, the response fragments group after signal slice is Z, where is the jth response in the response fragment group. Among them, Z can be obtained according towhere p is a positive integer. Due to the change of the rotational speed, at this time, A1, A2, …, AL are time-domain unequal response fragments, but they are vibration data of the bearing run equal angle (Δθ).

2.2. Demodulating the Signal Fragments

In this paper, taking the resonance sparse decomposition as an example, the fault mode signal is demodulated from the fault response fragments Ai, and the noise is eliminated. Since Selesnick [21] proposed a resonance-based signal sparse decomposition method in 2011, it has been widely used in fault diagnosis. Because of it, the band overlap problem can be solved. In signal resonance sparse decomposition using two high- and low-quality factor wavelets to fit the signal, one can find the components of the vibration signal that exhibit different waveform characteristics. It consists of two parts: the quality factor adjustable wavelet change and the morphological component analysis, which can decompose a signal into a high-resonance component and a low-resonance component.

Among them, the setting of high and low quality factors is very important for the ability to fit well. The wavelet corresponding to the quality factor should be as similar as possible to the waveform of the response that needs to be extracted. The quality factor can be estimated by Q = center frequency/bandwidth. When the quality factor Q and the redundancy r are set, the low-pass scale factor α and the high-pass scale factor β are calculated as follows [22]:

The corresponding wavelet center frequencies and bandwidths of different decomposition levels are different. The center frequency fc and bandwidth of each subband can be expressed as follows:where j is the number of levels in which the wavelet is located.

Whether it is a high-resonance component or a low-resonance component, they are linear sums of a series of wavelets of different levels. The center frequency of the jth level wavelet can be estimated by using equation (6). The larger j, the smaller the center frequency of the wavelet (because α < 1). It is necessary to ensure that the wavelet function library has wavelets whose center frequency is close to the oscillation frequency of the component to be extracted. However, if j is too large, the amount of calculation will increase, and for a signal of length N, the maximum number of decomposition stages is determined by the following equation:where represents the largest integer that does not exceed the number.

2.3. Mapping the Signal Fragments into Equal-Sized Mesh Phase Portraits

The mesh phase portraits method is proposed based on phase-space reconstruction and box-scoring. The aim is to overcome the effects of the lengths and energy significant differences between fault response fragments at variable speeds. The mesh phase portraits method is to map a fault response fragments into a two-dimensional mesh phase portrait to realize representation of the fault response.

2.3.1. Signal Fragments Mapping

The mesh phase portraits method consists of two parts. In the first part, the phase-space reconstruction is used to reconstruct the two-dimensional phase trajectory portraits of the fault response fragments. In the second part, the phase portraits are meshed, and the fault response trajectory map fills the part of the grid units. Then, we calculate the distance between the filled grid cell and the line y = x. The distance is filled as a weight value into the corresponding grid unit. Finally, mesh phase portraits will be obtained.

Step 1 (the fault response reconstruction using the phase space). For convenience of presentation, the fault response fragments is represented as A = (a1, a2, …, aK). According to the embedding method, the response fragment A is transformed into an m-dimensional phase-space vector by delay time coordinates, m is an embedding dimension, and τ is the delay time. According to the study of Hou et al. [20], the mesh phase portraits are developed in two-dimensional phase space, so m = 2.

Step 2 (meshing the phase space). The two-dimensional phase space of the fault response fragments is meshed, and the mesh size is 2M × 2M. The value of M is related to the sampling accuracy, and the value setting method will be explained below.

Step 3 (weight calculation). With the attenuation of the fault signal amplitude (such as impact-type fault signals), the distance from the trajectory position of the phase space to the straight line (y = x) becomes smaller and smaller. Hence, the distance between the trajectory and line y = x can be used to represent the change in energy. The distance is used as a weight, and it is assigned to the grid unit where the phase-space track is located. The weight is calculated as follows:where is the coordinate of the grid unit.

2.3.2. Selection of the Number of Grid Units

The mesh phase portraits of the fault response fragments are two-dimensional matrices (2M × 2M) with weights. The larger the number of grid units, the more the signal details will be retained. But too large grid units will increase the unnecessary calculation amount. The value range of M must satisfy the mesh phase portraits to clearly express the fault response fragments and reduce the calculation amount as much as possible. The number of grid units is limited by the value of the full-scale voltage Vfsr and the peak of fault response Vfault. Vfault needs to be much larger than the resolution of the mesh phase portraits. According to equation (9), the minimum value of M can be calculated. The minimum M value of the experiment involved in this paper is 5:

2.3.3. Delay Parameter Selection and Incremental Sampling

For phase-space delay, there is no specific requirement and selection method in the embedding theory. The embedding theory is applicable to noise-free and infinite-length signals [22]. Because, in the process of actually collecting signals, it is inevitable that noise will be mixed in and the length of the processed data is not infinite, it is difficult to meet the requirements. Different delays τ affect the amount of information contained in mesh phase portraits. If τ is too large, the mutual information of the delayed coordinate elements on the mesh phase portraits will be lost, and the signal traces will appear to be folded. If τ is too small, the redundancy is large, and the information of the original attractor is small in the mesh phase portraits. In the mesh phase portraits, the trajectory shrinks toward the main diagonal.

There are many methods of τ setting, such as the mutual information method [14] and the average autocorrelation method (AUD) [23]. In order to obtain the value of τ, the mutual information method is used. We find that as τ varies from one to n − 1 sample periods, the cross-correlation value does not show a local minimum, so the optimal τ cannot be determined. τ calculated by AUD is 0.075 ms, and the signal sampling period of the fault response fragments is 0.05 ms. Then τ is 1.5 sampling periods. For discrete data, a large error occurs when τ takes 1 or 2 sample periods. This indicates that the fault response fragment sampling frequency is insufficient. Hence, the sampling frequency is increased by interpolation. As shown in Figure 2, they are the fault response fragments mesh phase portraits of 1–10 times fs. As can be seen from Figure 2, a clear reconstruction trajectory can be observed at a sampling frequency of 5 times or more. After increasing the sampling in this paper, the sampling frequency is >6 times fs.

With τ taking different values, the mesh phase portraits of the fault response fragments are shown in Figure 3. It can be known that when τ = 0.075 ms, the mesh phase-space track area accounts for the largest proportion in the portrait. The value of τ is also 0.075 ms that is calculated by the AUD method and is exactly equal to 1/4 of the resonance period in the fault response fragments. That is, τ = 1/(4 × fr), where fr is the resonant frequency of the fault response.

In order to further explore the relationship between the optimal value of τ and the resonance period, the τ value of the different resonance frequency fault response fragments acquired using the AUD method is compared. The calculation results are shown in Table 1. And the fault response fragment mesh phase portraits are shown in Figure 4. We can know from Table 1 and Figure 5 that the AUD method yields a value of τ equal to about 1/4 times the resonance period. Hence, the value can be quickly determined by the resonance period, thereby reducing the amount of calculation.

2.3.4. Advantage of Variable Speed Response Mapping in Mesh Phase Portrait

The differences of device fault response fragment between the variable speed and the constant speed are as follows: (a) the period of the fault response is changed when the speed is changed. The length of the fault response fragments are not the same in the time domain. The faster the speed, the shorter the length of the fragments; (b) the energy of the fault response is different when the speed is changed. And the faster the speed, the greater the energy. It is expressed by the greater amplitude of the vibration. In the mesh phase portrait method, the geometrical characteristics of the intrinsic structure of the motive system are unchanged. Therefore, the reconstructed trajectory in the mesh phase portrait is not related to the data length. It only increases or decreases the content of the mesh phase portrait. And it does not change the characteristics of the mesh phase portrait according to the length of the data. Then, the mesh phase portrait method maps response fragments of different energy into an equal number of grid units, thereby eliminating the energy interference of the rotational speed. It can be known that the mesh phase portrait method can well eliminate the influence of the variable speed. The correspondence between the fault response fragments and the mesh phase portrait is shown in Figure 5. The mesh phase portraits of the fault response fragments at different rotational speeds are shown in Figure 6.

2.4. Calculation of the Cosine Similarity of the Mesh Phase Portraits of the Signal Fragments

The similarity comparison of the fault response fragments is achieved by calculating the mesh phase portraits and determining whether they belong to the corresponding type of fault. Since the response fragment waveforms of the same fault are not exactly the same, there are individual differences. This results in similar trajectory shapes and weights in mesh phase portraits, but not equal. In order to eliminate the influence of the individual difference of the fault source on the similarity calculation, a 3 × 3 Gaussian filtering is performed on the mesh phase portrait. The cosine distance, also known as the cosine similarity, is similar to the cosine of the angles of the two vectors in space. The cosine value is closer to 1, the angle is closer to 0 degrees, and the two samples are more similar. Here, the two-dimensional mesh phase portrait is vectorized, and then the cosine value is calculated. After a comparison between 100 homologous fragments and 100 nonhomologous ones, the similarity threshold can be set as 0.5 which will be used in the subsequent experiments.

3. Simulated Analysis

3.1. Signal Simulation Validity Verification

In order to verify the effectiveness of the proposed method, the simulation signal is verified. The impact caused by bearing failure can be modeled by the single-degree-of-freedom mass-spring-damper impact response [24, 25]:

Considering the slip effect of the bearing rolling elements, the fault bearing vibration signal model is as follows [26, 27]:where Am is the amplitude of the m-th fault impulse response, u(t) is the unit step function, β is the structural attenuation coefficient, ωr is the resonance frequency excited by the bearing fault, and tm is the time when the m-th impact occurs and can be calculated as follows:where τd represents the error between the fault impact intervals caused by the sliding of the rolling elements, and its value is generally 0.01∼0.02. And t0 = 0. n represents the number of impact responses produced by each revolution of the shaft. f(t) represents the frequency of the different times, and f(t) = 10t + 15. The parameters of the simulation signal are shown in Table 2.

According to the above model, the vibration signal of the bearing outer ring fault is simulated. The fault characteristic order of outer ring (FCOO) is 3.6. The simulation signal waveform is shown in Figure 7(a), and the rotation speed is as shown in Figure 7(b).

The simulated fault bearing vibration signal processing under the speed increase is processed by the mesh phase portrait similarity comparison method proposed in this paper. The simulated signal time is 3 seconds. The speed range is 15–45 r/s. According to the fast-delay determination method, the delay time τ is 0.068 ms, and the sampling incrementing multiple is 8. The number of fault response fragments in this experiment was 20. Secondly, the selected fault response fragments are demodulated by the resonance sparse decomposition method and then transformed into mesh phase portraits, as shown in Figure 8(a). Finally, the similarity between each response fragments and other response fragments is calculated, as shown by the red curve in Figure 6. The mean similarity of each fault response with other responses is 0.64, which is greater than threshold 0.5. Simulation results show the effectiveness of the proposed method.

In order to further verify whether the method has a misjudgment, the similarity of the fault response fragments of the nonhomologous vibration source is also calculated. The experiment still uses the simulated bearing fault vibration data of the increased speed as the object. According to the bearing inner ring (FCOI) 6.1, the fault response fragments are intercepted. Because there is only an outer ring fault in these data, there is no inner ring fault, so the inner ring fault response fragments are nonhomologous. The experimental results are shown as the black curve in Figure 6. The similarity of the nonhomologous vibration source fault response fragments has large fluctuations, and the similarity is low. The nonhomologous average similarity is 0.47, which less than threshold 0.5. This is not consistent with the actual situation, and the effectiveness of this method is verified from the reverse side.

3.2. Comparative Analysis

In order to verify the superiority of the method, it is compared with OT and DTW methods. The comparison method 1 is OT, which eliminates the effect of the rotational speed by resampling the vibration signal in the angular domain to keep the length of the fault response fragments consistent. The comparison method 2 is DTW. This method can directly calculate the similarity between signals of unequal length, but the time complexity is high. The similarity calculated by the DTW is distance, and the smaller the distance, the more similar it is.

The experimental results are shown in Figure 9. It can be seen from Figure 9(a) that there are more overlapping regions with the similarity between the homologous vibration source response fragments and the nonhomologous vibration source response fragments by OT, and the discrimination effect is not obvious. The similarity average of the homologous vibration source response fragments is low, which is inconsistent with the facts. The DTW also cannot distinguish the homologous and nonhomologous vibration source response fragments.

In order to verify the time effectiveness of method, the comparison methods 1 and 2 were calculated. The platform for running the three methods is Matlab2014a, and the hardware is ThinkpadE430c notebook computer. The calculation result of the consumption time is shown in Figure 10. The FMPP took 1.14 s. The calculation time was reduced by 5% relative to the OT. Compared with DTW, the calculation time is reduced by 93%. In general, FMPP has low time complexity compared to the other two comparison methods.

4. Experiment Analysis

4.1. Validity Verification

In order to verify the effectiveness of the proposed algorithm, this experiment collected the signal of the outer ring fault crack bearing during increase in the speed, and the sampling frequency is 24 kHz. The experimental device and the acquired vibration signal waveforms are shown in Figures 11 and 12. The acceleration sensor is located at the bracket on the upper side of the bearing, and the tachometer is mounted on the shaft end. The acquisition device is the YE6231 acquisition card and its supporting software. The type, the geometric parameters, and the outer characteristic fault characteristic order of fault bearing are shown in Table 3. The outer ring fault characteristic order FCOO and FCOI can be calculated according to the geometric parameters of the bearing [18]:where nr is the number of rolling elements and D1 is the diameter of the rolling element. D2 is the diameter of the pitch circle, and β is the contact angle.

The fault bearing vibration signal processing under increase in the speed is processed by the mesh phase portrait similarity comparison method proposed in this paper. In this experiment, the fault response fragments’ speed range is 50 r/s–76 r/s, as shown in Figure 13(b). The number of fault response fragments is 20. The delay time τ is 0.042 ms, and the sampling incrementing multiple is 8. The mesh phase portraits of the homologous vibration source fault response fragments are shown in Figure 14(a). The similarity between each response fragments and other response fragments is calculated, as shown by the red curve in Figure 13(a). The average similarity of each fault response with other responses is 0.57, which is greater than threshold 0.5. The test results show that the FMPP method is effective in measuring the similarity for fault response under variable speed.

In order to further verify whether the method has a misjudgment, the similarity of the fault response fragments of the nonhomologous vibration source is also calculated. The experiment still uses the bearing fault vibration data of the increased speed as the object. According to the FCOI value of 4.452, the fault response fragments are intercepted. Because there is only an outer ring fault, and there is no inner ring fault, the inner ring fault response fragments is nonhomologous response fragments. The experimental results are shown as the black curve in Figure 13(a). The similarity of the nonhomologous vibration source fault response fragments has large fluctuations, and the similarity is low. The average similarity is 0.38, which less than threshold 0.5. This experiment shows the validity of FMPP from the reverse side.

4.2. Comparative Analysis

In order to verify the superiority of the method, it is compared with OT and DTW methods. The experimental results are shown in Figure 15. It can be seen from Figure 15(a) that there are more overlapping regions with the similarity between the homologous vibration source response fragments and the nonhomologous vibration source response fragments by OT, and the discrimination effect is not obvious. The similarity average of the homologous vibration source response fragments is low, which is inconsistent with the facts. It can be seen from Figure 15(a) that the DTW also cannot distinguish the homologous and nonhomologous vibration source response fragments.

In order to verify the time effectiveness of the method, the comparison methods 1 and 2 were calculated. The result of the consumption time is shown in Figure 16. FMPP took 1.12 s. The calculation time was reduced by 84% with respect to the DTW. Compared with the OT method, the calculation time is reduced by 19%. In general, FMPP has lower time complexity compared to OT and DTW methods.

5. Conclusion

Based on phase-space technology and box-scoring calculation, a similarity measurement method FMPP for fault response fragments at variable speed is proposed for the first time. And the following conclusions can be obtained by the experiments:(1)The homogenous bearing fault response fragments at variable speed have high geometric invariance in mesh phase portrait. And it can be used to measure the similarity for response fragments with different length and energy.(2)The simulation analysis and test analysis show that the FMPP method can accurately measure similarity to fault response segments under variable speed. Compared with OT and DTW, FMPP has a better ability to distinguish homologous and nonhomologous response fragments. And FMPP has less time-consumption.(3)FMPP can quickly set the value of delay τ by 1/4 resonance period, thus reducing the calculation amount.

Data Availability

The test data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interests regarding the publication of this paper.

Acknowledgments

This study was supported by the Beijing Natural Science Foundation (3192025). The support is greatly appreciated.