A Wavelet Bicoherence-Based Quadratic Nonlinearity Feature for Translational Axis Condition Monitoring

The translational axis is one of the most important subsystems in modern machine tools, as its degradation may result in the loss of the product qualification and lower the control precision. Condition-based maintenance (CBM) has been considered as one of the advanced maintenance schemes to achieve effective, reliable and cost-effective operation of machine systems, however, current vibration-based maintenance schemes cannot be employed directly in the translational axis system, due to its complex structure and the inefficiency of commonly used condition monitoring features. In this paper, a wavelet bicoherence-based quadratic nonlinearity feature is proposed for translational axis condition monitoring by using the torque signature of the drive servomotor. Firstly, the quadratic nonlinearity of the servomotor torque signature is discussed, and then, a biphase randomization wavelet bicoherence is introduced for its quadratic nonlinear detection. On this basis, a quadratic nonlinearity feature is proposed for condition monitoring of the translational axis. The properties of the proposed quadratic nonlinearity feature are investigated by simulations. Subsequently, this feature is applied to the real-world servomotor torque data collected from the X-axis on a high precision vertical machining centre. All the results show that the performance of the proposed feature is much better than that of original condition monitoring features.


Introduction
The translational axis is one of the most important subsystems in modern machine tools, which has been widely used for transforming rotational motion into linear motion, due to its advantages of low sensitivity to inertia variations, good positional accuracy and high driving speeds [1]. Starting from a new condition, each part of the translational axis system is properly manufactured and installed into the system without any errors, and is operated at a particular level of performance. As their service life progresses, faults such as material fatigue, abrasion or adhesion in the components always leads to a performance degradation, which means the parts can no longer meet their original service requirements [2]. This performance degradation of a translational axis system may cause some problems. First of all, it may produce distance errors between the tool tip and the workpiece in the machining process and thus lead to the loss of the product qualification [3,4]. Secondly, any faults occuring in components will lower the control precision [1,5]. It is, therefore, very essential to investigate the technology for translational axis system condition monitoring. Condition-based maintenance (CBM) has been considered as one of the advanced maintenance schemes to achieve effective, reliable and cost-effective operation of machines and systems [6]. In a CBM system a machine's condition before it fails can be monitored and predicted, and optimal maintenance actions can be scheduled to improve reliability and reduce the possibility of breakdowns [7]. In the past few years, condition monitoring methods for machine tools have been widely studied as the key research topics in CBM [8][9][10]. Ottewill and Orkisz [11] used synchronously averaged electric motor signals to monitor gearboxes, showing that this method was extremely adept at locating gear tooth defects. Hwang and Lee [12] studied a condition monitoring method based on the change likelihood of a stochastic model, which was successfully applied to the monitoring of machine condition and weld condition. Sinha and Elbhbah [13] focused on the monitoring of rotating machines by enhancing the computational effort in signal processing to reduce the number of sensors per bearing pedestal. Subsequently, Hassan and Bayoumi et al. [7] proposed a Fourier bicoherence-based index for health condition monitoring of helicopter drive trains. Yu [14] proposed a Generative Topographic Mapping (GTM) and contribution analysis-based method to perform machine health degradation assessment and monitoring. Generally, most of the current methods are based on identifying some significant change in vibrations of the system. However, current vibration-based maintenance schemes cannot be employed in the translational axis system directly, due to its complex structure and the inefficiency of the commonly used condition monitoring features [6,15]. In this paper, a wavelet bicoherence-based quadratic nonlinearity feature is proposed for translational axis condition monitoring by using the torque signature of the drive servomotor. Firstly, the quadratic nonlinearity of the servomotor torque signature is discussed, and then, a biphase randomization wavelet bicoherence is introduced for its quadratic nonlinear detection. On this basis, a quadratic nonlinearity feature is proposed for condition monitoring of the translational axis. The key idea is that as the performance of the translational axis degrades, the torque signature will tend to be more nonlinear, which will result in the generation of new frequencies [16,17]. Those new frequencies are phase coupled to the original characteristic frequencies of the translational axis, and the wavelet bicoherence is a sensitive detector for identifying such phase coupling [18]. The remaining parts of this paper are organized as follows: in Section 2, a theoretical basis is set up to show the quadratic nonlinearity of motor torque signal, under the fault condition of translational axis system. After that, a biphase randomization wavelet bicoherence technique is introduced in Section 3, and a quadratic nonlinearity feature is proposed for condition monitoring. The performance of the proposed quadratic nonlinearity feature is investigated through a simulation in Section 4. In Section 5, this feature is applied to the translational axis condition monitoring of a vertical machining centre in an in-service environment. Finally, conclusions are given in Section 6.

Quadratic Nonlinear Model of Servomotor Torque
As shown in Figure 1, the translational axis system is often composed of an AC servomotor, a ball screw, a table, a guideway and other rotating components (such as support bearing and reducer). These parts work cooperatively to guarantee the table is constrained to move axially on guideways once the servomotor torque overcomes the disturbance and table inertia. In order to calculate the torque of the AC servomotor, the three-phase AC currents of the servomotor are measured and converted into DC current by using the D-Q transformation as follows [19,20] where u I , v I and w I are the u-, v-and w-phase AC current respectively, p n is the number of poles,  is the relative angle between the u and Q-axes. The D-Q transformation method transforms rotating coordinates to the D-axis and Q-axis, in fixed coordinates, and the converted DC current q I is directly related to the motor torque, d I is maintained at zero. During the cutting process, the measured disturbance torque is composed of frictional torque and the torque caused by the cutting force (as shown in Figure 1), which is proportional to the motor DC current q I according to the ratio of torque constant T  . Hence, the total motor torque can be modeled as: where q I is the equivalent DC current of the drive motor, T  is the motor torque constant, J is the equivalent total moment of inertia of the translational axis system,  is the rotational speed of rotor, f T is the friction torque, and c T is the cutting torque. When the table is fed steadily without cutting, the measured disturbance torque only represents the frictional torque of guideways, ball screw and other rotating components (such as support bearing and reducer). Then, 0 c T  , and Equation (2) is modified as: During the assembly of machine tools, a preload must be exerted to minimize the clearance of the ball screw and support bearing to improve the dynamical operational behavior [21]. The increase of the preload will raise the stiffness of the ball screw and support bearing [22]. However, increasing the preload will also increase the friction on the ball screw and support bearing. During the serving cycle, various frictional components such as rolling, bore and sliding friction will lead to different faults such as material fatigue, abrasion or adhesion, and ultimately potentially result in additional torque. Hence, as shown in Equation (3), the torque under no load condition can be used as an attribute of the condition monitoring of the translational axis system [22]. Commonly, with the degradation of the translational axis system, the servomotor torque signature will tend to become more nonlinear. The nonlinearity mechanism is often very complex. In this article, the quadratic nonlinear (a strong indicator of nonlinear signal) of the torque signature is considered [23,24]. The quadratic nonlinear system response, x(t), for an input y(t) is: where a and b are two constants that define the linear and nonlinear components of the system input and output relationship, n(t) is additive Gaussian white noise. The degradation of the part in the translational axis system will lead to an additional friction. Suppose that the additional friction y(t) is a cosinusoid wave, and due to nonlinearity and harmonic production, a suitable input is two mixed cosinusoids: To avoid the additional friction, an additional torque should be provided by the servomotor. Substituting Equation (5) into Equation (4) and applying trigonometric identities, the additional torque response is given as: In the frequency domain, the continuous wavelet transform (CWT) spectrum of Equation (6) can be obtained as: As shown in Equation (7), the nonlinearity signals at 1 2 f and 2 2 f are denoted as the harmonics, and those at 21 ff  and 21 ff  are the nonlinear interaction components. In addition, the phase of the signals caused by the nonlinearity satisfies the relationship, 1   is the phase of the new interaction component. This is the quadratic phase coupling (QPC) due to system nonlinear behavior, which provides a way to distinguish signals at a particular frequency. Because only the components contributed to the QPC are of interest, all other terms are neglected, and Equation (6) is modified as [18,24]: The simplified QPC model is employed widely for condition monitoring [18]. By checking the phase coupling frequency components, various faults can be recognized with high probability. Bispectrum and bicoherence are effective in detecting the QPC in this type of nonlinear signal [25,26], but, they are not appropriate for the signal associated with short-time duration nonlinear interactions. Recently, a wavelet bicoherence (WB) method combining benefits of both wavelet transform and bicoherence analysis is proposed for nonlinear and non-stationary signal analysis [27][28][29]. Compared with the Fourier bicoherence, WB has its advantage in preservation of temporal information, thus it is more suitable for the instantaneous servomotor torque analysis.

Biphase Randomization Wavelet Bicoherence
The WB is firstly introduced by van Milligen [30]. Given a signal x(t), the continuous wavelet transform (CWT) of x(t) is defined as the convolution of x(t) with the scaled and normalized wavelet. We write: where the * indicates the complex conjugate,  () t is an admissible mother wavelet, s is the scale variable, t is the time shift variable, f is the equivalent Fourier frequency.
In an analogous definition to the usual Fourier bispectrum, the wavelet bispectrum is defined as follows [30]: WT where the * indicates the complex conjugate, T is the finite time interval of the signal, frequency values f 1 , f 2 and f 3 satisfy the relationship f 3 = f 1 + f 2 . The wavelet bispectrum is a complex value. So, it can be also expressed by its magnitude A(f 1 , f 2 ) and biphase  (f 1 , f 2 ) as follows: 12 ( , ) , (11) where the biphase  (f 1 , f 2 ) can be calculated by, are phases uniformly distributed within ( , ]   , obtained by CWT. Conventional wavelet spectrum based WB is incapable of eliminating the spurious peaks coming from long coherence time waves and non-QPC waves [23,31]. As a result, direct application may cause incorrect results. To overcome this problem, a biphase randomization wavelet bispectrum is introduced as [32,33]: where  12 ( , ) iR f f e is the biphase randomization term,  (f 1 , f 2 ) is the biphase obtained by Equation (11) f At the same time, a BRWB-based quadratic nonlinearity feature, namely the maximum eigenvalue of the BRWB (MEBRWB feature) [29,34], is proposed to make condition monitoring of the translational axis system more convenient. Since the BRWB is a symmetrical matrix to the main diagonal, denoted as , it can be decomposed as: V is the eigenvetor corresponding to k  . The value of eigenvalue is proportional to the amount of correlation in the direction of their associated eigenvectors, and the maximum eigenvalue linearly increases with the increase of the correlation strength in a cluster. Therefore, the single-value MEBRWB feature can be defined as a global phase coupling feature of the BRWB.

Estimation of BRWB and Its Quadratic Nonlinearity Features
The estimation of the BRWB, the SBRWB and MEBRWB features is detailed in the following steps as illustrated in Figure 2. Compute the CWT of each epoch by (9) Calculate the wavelet bispectrum by (12), and the result is retained as a sample for ensemble average Estimate the biphase randomization wavelet bispectrum of each epoch. Estimate the BRWB by (13) Estimate the SBRWB and MEBRWB feature by (14) and (15) k times Generate a random variable R, and multiply of both R and to generate a new biphase 1 2 ( , ) f f 


Step 1: Divide the servomotor torque signal into n epochs.
Generally, in computing of BRWB, the signal with long data record length is firstly divided into a series of epochs, and then, the reliable bicoherence estimation of the signal is obtained by averaging the BRWB of all epochs. In this study, the BRWB of torque signal is calculated based on smallest meaningful signal segment, due to its non-stationary property. Because that the signal segment of 1/ PFL F ( PFL F is the passing frequency of ball screw lead) length often covers all characteristic frequencies of translational axis system, the torque signal is divided into n epochs with each period equal to 1/ PFL F , with the overlap of 75% [29].
Step 2: Compute the wavelet transform ( , ) W f t  of each epoch by Equation (9).
For the torque signals, impulse components often correspond to the degradation of the translational axis system, thus the basic wavelet used for feature extraction should be similar to an impulse [35]. For this reason, the Morlet wavelet is chosen for feature extraction in this study. Since the wavelet bicoherence calculation involves three scales, the choice of wavelet parameter should tradeoff between the time width of analysis wavelet function and best suited individual scale. Here, takes  0 to 6 to satisfy the admissibility condition, and takes  within 3 to 8 to obtain the appropriate half power time-width of the mother wavelet.
Step 4: Generate a random variable R within ( 12 ( , ) ff. At the same time, the biphase randomization term is also used to eliminate the spurious bicoherence coming from non-QPC waves.
In this step, the biphase randomization term  Step 6: Estimate the biphase randomization wavelet bispectrum of each signal epoch.
In this step, the biphase randomization wavelet bispectrum of this signal epoch is estimated by ensemble average operation of the 100 data samples obtained in Step 5. The spurious bispectrum coming from non-QPC waves can be eliminated by the ensemble averaging operation. After normalization, the BRWB of this signal epoch can be obtained, and the result is retained as the samples for reliable BRWB estimation. From this step, n BRWB data samples can be obtained for reliable BRWB estimation.
Since n signal epochs are used for BRWB estimation, here, the reliable BRWB is estimated by averaging n BRWB data samples obtained in Step 6. The biphase component of each BRWB sample has been made independent over each other in previous steps. Thus, by the average operation, the spurious bicoherence coming from long coherence time waves can be eliminated automatically. Commonly, 10-20 signal epochs are sufficient to obtain a reliable BRWB estimation [32].
The SBRWB feature is calculated by Equation (14), and an eigenvalue series can be obtained from the Equation (15). By extracting the maximum value of this series, a single-value MEBRWB feature is obtained for global QPC measurement of the torque signal.

Simulations
According to the QPC model introduced in Equation (8), the simulation signal model is as follows [24,35] (2  ) (cos(2 ( ) ( ))) (1 )(cos(2 ( ) ))] where f 1 and f 2 are the coupled frequency,  1 and  2 are the initial phase distributed within   ( , ] uniformly, n(t) is white Gaussian noise with zero-mean and unit variance, T is the period of the signal and A is the coupling coefficient. The simulation signal in Figure 3a is configurated as f 1 = 80 Hz, f 2 = 270 Hz, A = 1 and T = 0.25 s with the sampling rate 1,000 Hz, and its wavelet scalogram is shown in Figure 3b.  It can be seen from Figure 4a that the estimated MEWB feature is increasing as the SNR increases, however, it doesn't change with the increase of the coupling coefficient A. Thus, the MEWB feature is not able to detect the coupling components at each SNR levels, when phase-coupled and non-phase-coupled are mixed together. In contrast, as shown in Figure 4b, the estimated MEBRWB feature increases with the increase of the coupling coefficient A, although it decreases with the decreasing SNR levels. In addition, the trend of the MEBRWB feature approximates to linear. Thus the proposed quadratic nonlinearity feature is able to capture changes in global QPC of the signal efficiently, particularly when the SNR is greater than 3 dB.

Description of the Experimental System
The proposed method is applied to monitor the translational axis of a high precision vertical machining centre. Experiments are conducted on the X-axis of this machining centre in an in-service environment, and the experimental time is the whole maintenance period of the X-axis. The experiment system is shown in Figure 5 and the illustration of the X-axis is shown in Figure 6. As shown in Figure 6, the translational axis system (X-axis) is composed of an AC servomotor, a reducer, a precision ball screw and a table. The actuation of this system is provided through the AC servomotor which is attached to the reducer using a diaphragm type coupling. The reducer is a three-stage gearbox attached to the precision ball screw also using a diaphragm type coupling (the teeth number of the each gear is shown in Figure 6). The precision ball screw with 16 mm pitch and 40 mm diameter is supported by two bearings, which drives a table supported on a guideway. Experiments are implemented during the whole maintenance period of the X-axis under the condition that the feed rate is 550 mm/min, when the table is fed steadily without cutting. The three-phase AC current signals of the servomotor are measured synchronously by three current sensors, and then, the motor torque is obtained by Equations (1) and (2). The total experimental time is 192 days, and the data is collected with an interval of approximate 26 days. Each data is sampled with sampling frequency 1,000 Hz. Actually, there are eight torque data samples, and each data sample is a collected data series containing 60,000 data points. The characteristic frequencies of the X-axis are given in Table 1. Figure 6. Illustration of the X-axis.  Figure 7 shows the calculated servomotor torque (left panels) of X-axis and its spectrum (right panels), at 192 days, 62 days 34 days and 2 days before maintenance, respectively. It can be seen from the left panels of Figure 7 that the servomotor torque fluctuates for each data sample. During the last 34 days before maintenance, some impulses appear in the last half of the ball screw travel. As shown in the right panels of Figure 7, the spectrum plot of each data sample shows that the same dominating spectrum peaks appear at the servomotor's rotary frequency and the first, second and third harmonics of the reducer meshing frequency, with very slight changes in the peaks' value. At the same time, a resonance appears at the frequencies band between 330-370 Hz, during the last 62 days before maintenance. Similar results occur in the spectrum of almost every data sample. Thus, it is not easy to describe the degradation progress by directly using the change of the characteristic frequencies. To describe the degradation progress, the proposed BRWB is used to detect the quadratic nonlinear phase coupling frequencies band of the torque signature. Here, the bandwidth from 0-500 Hz is selected for the bicoherence calculation. Figure 8a-d showS the BRWB of the torque data samples at 192 days, 62 days, 34 days and 2 days before maintenance, respectively. As shown in Figure 8a, it is found that the dominant phase coupling peaks appear at the bifrequency (32 Hz, 64-350 Hz), which illustrates that, at the beginning of the service life progress, the quadratic nonlinear interaction only occurs between the 1th, 2th and 3the harmonics of the reducer's meshing frequency. As the service life progresses, new phase coupling presents at the bifrequency (64 Hz, 270 Hz) and (140 Hz, 230 Hz) as shown in Figure 8c,d, which illustrates that quadratic nonlinear interaction between the reducer's meshing frequency and the resonance frequency band results from the degradation of the X-axis. In addition, the peaks value at the new bifrequency (64 Hz, 270 Hz) and (140 Hz, 230 Hz) are increasing as the degradation progresses. To further investigate the interactions, the phase coupling distribution of each data sample is estimated by using the SBRWB feature, as shown in Figure 9a, and the global phase coupling of each data sample is estimated by the MEBRWB feature as shown in Figure 9b. The differences among the BRWB of each data sample can be directly revealed by the SBRWB feature. The SBRWB features of the resonance frequency band, servomotor rotary frequency harmonics and the second harmonic of reducer meshing frequency are increasing significantly during the service cycle. At the same time, the global phase coupling of each data sample is also increasing significantly during the service cycle, as shown in Figure 9b. For convenience, the single-value MEBRWB feature of the eight torque data samples is calculated to monitor the condition of X-axis. The MEBRWB feature is compared with other commonly used feature during the whole maintenance period of the X-axis, as shown in Figure 10. Figure 10a-e shows the kurtosis, root mean square, summation of the rotary frequency harmonics, summation of the meshing frequency harmonics and peak values of the resonance frequency band of the eight measured torque data samples, respectively. Results show that the commonly used features don't perform well in describing the progress of X-axis degradation. By contrast, as shown in Figure 10f, the value of MEBRWB feature starts to increase at 143 day before maintenance. It keeps increasing until the translational axis system maintenance because of the degradation of the components, which physically can be interpreted as increased additional torque quadratic nonlinear interaction due to the performance degradation of the X-axis. Thus, this trend can be used as precocious indication of performance degradation of the translational axis of the machine tools.

Conclusions
A BRWB based quadratic nonlinearity feature is established for translational axis condition monitoring of the machine tools, which shows that it is possible to perform condition monitoring by using servomotor torque signature. Some conclusions are drawn as follows: (1) A BRWB is established to overcome the problem of current WB, which can eliminate the spurious peaks coming from long coherence time waves and non-QPC waves efficiently. Based on the proposed BRWB, a quadratic nonlinearity MEBRWB feature is proposed for translational axis condition monitoring. Numerical example results show that the global QPC of the simulation signal can be tracked approximately linearly by the proposed feature, especially when the SNR is greater than 3 dB.
(2) The proposed quadratic nonlinearity MEBRWB feature is also used to study the torque data collected from a high precision vertical machining centre. Experimental results illustrate the robust experimental performance of the proposed feature, compared with commonly used features. The advantage of MEBRWB feature is that it can exploit the true global QPC of the signal at different frequencies, which contains useful additional information for detection of quadratic nonlinear phenomena induced by mechanical faults. Potentially, this single-valued feature can be used in prognostic models for remanding useful life estimation of mechanical components.