A Similarity Comparison Method of Homologous Fault Response Fragments under Variable Rotational Speed

Effective fault detection and diagnosis (FDD) for rotatingmachinery is always a focus issue in improving the prognostic and health management (PHM) of the equipments. +e 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. +irdly, 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.+e experiments proved the effectiveness of the proposed method.


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 timedomain 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 [2][3][4][5]. Since each CW is a sequence, the clustering problem of CW is naturally transformed into the similarity calculation of sequence. e similarity measurement of time series is always the focus and difficulty. e 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. e 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. e 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. ere 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 [10][11][12] 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. e vibration signal will be stretched or compressed (deformed) when it is sampled in the angular domain by the OTmethod [16]. e 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. e 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. e 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 boxscoring 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 phasespace and box-scoring calculation is proposed, which is called the fast-meshed phase portrait (FMPP) method. e 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.

Fault Signal Fragments' Similarity Calculation Method
is section describes the four steps which are used to calculate the fault signal fragments' similarity. e frame is shown in Figure 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. e fault response fragment which contains a complete impulse response reflects the health of the component. e fault response fragment which contains an impulse response reflects the complete health information of the component. e 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 F COO . When the bearing running 1/F COO 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 F CO 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 as where f s is the sampling frequency and Δt � 1/f s . e vibration signal is then sliced at equal angle ∆θ, and ∆θ is an interval angle of two fault responses. en, the ∆θ calculation method is as follows: e vibration signal is represented by X � (x 1 , x 2 , . . . , x N ). When the length of a single response fragment is n, the response fragments group after signal slice is Z, where A j is the jth response in the response fragment group. Among them, Z can be obtained according to where p is a positive integer. Due to the change of the rotational speed, at this time, A 1 , A 2 , . . ., A L are time-domain unequal response fragments, but they are vibration data of the bearing run equal angle (Δθ).

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 A i , 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. e wavelet corresponding to the quality factor should be as similar as possible to the waveform of the response that needs to be extracted. e 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]: e corresponding wavelet center frequencies and bandwidths of different decomposition levels are different.
e center frequency f c and bandwidth B w 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 lowresonance component, they are linear sums of a series of wavelets of different levels. e center frequency of the jth level wavelet can be estimated by using equation (6). e 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.  proposed based on phase-space reconstruction and boxscoring. e aim is to overcome the effects of the lengths and energy significant differences between fault response fragments at variable speeds. e 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.

Signal Fragments
Mapping. e mesh phase portraits method consists of two parts. In the first part, the phasespace 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. en, we calculate the distance between the filled grid cell and the line y � x. e 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 � (a 1 , a 2 , . . ., a K ). 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). e two-dimensional phase space of the fault response fragments is meshed, and the mesh size is 2 M × 2 M . e 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. e distance is used as a weight, and it is assigned to the grid unit where the phasespace track is located. e weight is calculated as follows: where [u(i), v(i)] is the coordinate of the grid unit.

Selection of the Number of Grid
Units. e mesh phase portraits of the fault response fragments are two-dimensional matrices (2 M × 2 M ) with weights. e 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. e 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. e number of grid units is limited by the value of the full-scale voltage V fsr and the peak of fault response V fault . V fault 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. e minimum M value of the experiment involved in this paper is 5:

Delay Parameter Selection and Incremental Sampling.
For phase-space delay, there is no specific requirement and selection method in the embedding theory. e 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. ere 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.
en τ is 1.5 sampling periods. For discrete data, a large error occurs when τ takes 1 or 2 sample periods. is 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 f s . 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 f s .
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. e 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. at 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. e 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.

Advantage of Variable Speed Response Mapping in Mesh Phase Portrait.
e 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. e length of the fault response fragments are not the same in the time domain. e 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. erefore, 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. en, 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. e correspondence between the fault response fragments and the mesh phase portrait is shown in Figure 5. e mesh phase portraits of the fault response fragments at different rotational speeds are shown in Figure 6.

Calculation of the Cosine Similarity of the Mesh Phase Portraits of the Signal Fragments.
e 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. is 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. e cosine distance, also known as the cosine similarity, is similar to the cosine of the angles of the two vectors in space. e 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.

Signal Simulation Validity Verification.
In order to verify the effectiveness of the proposed method, the simulation signal is verified. e impact caused by bearing Considering the slip effect of the bearing rolling elements, the fault bearing vibration signal model is as follows [26,27]: where A m 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 t m 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 t 0 � 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. e 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. e fault characteristic order of outer ring (F COO ) is 3.6. e simulation signal waveform is shown in Figure 7(a), and the rotation speed is as shown in Figure 7(b). e simulated fault bearing vibration signal processing under the speed increase is processed by the mesh phase portrait similarity comparison method proposed in this paper. e simulated signal time is 3 seconds. e 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. e 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. e 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. e experiment still uses the simulated bearing fault vibration data of the increased speed as the object. According to the bearing inner ring (F COI ) 6.1, the fault response fragments are   Shock and Vibration 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. e experimental results are shown as the black curve in Figure 6. e similarity of the nonhomologous vibration source fault response fragments has large fluctuations, and the similarity is low. e nonhomologous average similarity is 0.47, which less than threshold 0.5. is is not consistent with the actual situation, and the effectiveness of this method is verified from the reverse side.

Comparative Analysis.
In order to verify the superiority of the method, it is compared with OT and DTW methods. e 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. e comparison method 2 is DTW. is method can directly calculate the similarity between signals of unequal length, but the time complexity is high. e similarity calculated by the DTW is distance, and the smaller the distance, the more similar it is. e 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. e similarity average of the homologous vibration source response fragments is low, which is inconsistent with the facts. e 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. e platform for running the three methods is Matlab2014a, and the hardware is inkpadE430c notebook computer. e calculation result of the consumption time is shown in Figure 10.
e FMPP took 1.14 s. e 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.

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. e experimental device and the acquired       vibration signal waveforms are shown in Figures 11 and  12. e acceleration sensor is located at the bracket on the upper side of the bearing, and the tachometer is mounted on the shaft end. e acquisition device is the YE6231 acquisition card and its supporting software. e type, the geometric parameters, and the outer characteristic fault characteristic order of fault bearing are shown in Table 3.
e outer ring fault characteristic order F COO and F COI can be calculated according to the geometric parameters of the bearing [18]: where n r is the number of rolling elements and D 1 is the diameter of the rolling element. D 2 is the diameter of the pitch circle, and β is the contact angle. e 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). e number of fault response fragments is 20. e delay time τ is 0.042 ms, and the sampling incrementing multiple is 8.
e mesh phase portraits of the homologous vibration source fault response fragments are shown in Figure 14(a). e similarity between each response fragments and other response fragments is calculated, as shown by the red curve in Figure 13  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. e experiment still uses the bearing fault vibration data of the increased speed as the object. According to the F COI 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. e experimental results are shown as the black curve in Figure 13(a). e similarity of the nonhomologous vibration source fault response fragments has large fluctuations, and the similarity is low. e average similarity is 0.38, which less than threshold 0.5. is experiment shows the validity of FMPP from the reverse side.

Comparative Analysis.
In order to verify the superiority of the method, it is compared with OT and DTW methods. e 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. e 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. e result of the consumption time is shown in Figure 16. FMPP took 1.12 s. e 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.

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) e 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. Data Availability e test data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interests regarding the publication of this paper. Shock and Vibration 11