A Method for Health Indicator Evaluation for Condition Monitoring of Industrial Robot Gears

: Condition monitoring of industrial robots has the potential to decrease downtimes in highly automated production systems. In this context, we propose a new method to evaluate health indicators for this application and suggest a new health indicator (HI) based on vibration data measurements, Short-time Fourier transform and Z-scores. By executing the method, we ﬁnd that the proposed health indicator can detect varying faults better, has lower temperature sensitivity and works better in instationary velocity regimes compared to several state-of-the-art HIs. A discussion of the validity of the results concludes our contribution.


Introduction
Industrial robots are a fundamental part of highly automated production systems, which can be found in the automotive or electronics industry [1]. Since they operate in complex production cells and as a part of linear production lines, robot malfunctions lead to long downtimes for repair or replacement and, hence, to increased costs. In particular, robot gear faults are responsible for the longest downtimes because they often require the replacement of the whole robot [2]. The condition monitoring (CM) of these gears offers the potential to resolve this issue. CM is the monitoring of an asset's health using sensor data. The health state represents a wear reserve before a failure occurs. This health state is quantified with a health indicator (HI). A significant monitored change in this health indicator can be used as a decision-making aid in the planning of maintenance actions [3].

State-of-the-Art
In recent years, different HIs based on vibration data for several industrial robot components, such as bearings, gears and motors, and their specific faults have been investigated. Furthermore, several approaches to cope with instationary signals in CM have been presented. The next two sections give a short overview of these topics followed by a section stating the contribution of our publication.

Vibration-Based Robot Condition Monitoring
A fault detection method was developed in [4], which first uses a novel phase-based, time-domain averaging method to remove the deterministic part of the vibration signal. Subsequently, the root mean square (RMS) and power spectrum entropy of the remaining residual signal are calculated as health indicators. A vibration signal based CM system for SCARA robots was implemented in [5], which in the first step uses statistical HIs of the time-domain signal to detect the occurrence of a defect and in the second step uses an artificial neural network to diagnose the fault type. A three-layer architecture for remote fault diagnosis of industrial robot gearboxes was proposed using vibration signals in [6]. In the diagnosis layer, the authors present a performance evaluation approach using a support vector machine (SVM), a remaining useful life prediction by a Markov changing robot axes' velocities, changing temperatures of the gears due to unbalanced robot utilization and unknown robot gear fault types. This is why we present a new HI for robot gear condition monitoring, which potentially copes with these characteristics. Furthermore, we propose a method to evaluate the suitability of HIs for the task of robot gear condition monitoring. We apply this method on the newly formulated HI and several HIs from the state-of-the-art.

Materials and Methods
This section is divided in two parts. First, the newly developed HI is presented. Afterwards, the methodologies to evaluate the HI's performance and data sets used in this context are explained.

Time-Frequency-Domain-Based Z-Score
The concept of the newly designed HI is based on two cornerstones. To deal with instationary velocity regimes, which are found in robot applications due to the typical movement patterns of a robot, the HI is based on time-frequency-domain data. Simultaneously, the HI must take into account a certain variance of this data due to environmental changes such as temperature fluctuations. This is realized by the concept of Z-scores, a common similarity measure from statistics [34]. The process to calculate the new HI is depicted in Figure 1. In detail, the new HI is based on high-frequency sampled acceleration sensor data. Data from one measurement are transformed to a time-frequency-spectrogram by usage of the STFT, which is calculated according to Equation (1). Here, τ and ω are time and frequency indices, x(n) is the time series signal of the vibration signal at timestep n and w is a windowing function with the length R.
To set up the HI, a certain number of vibration signal spectrograms must be collected for the robot to capture its signal signature in a healthy state with its stochastic variations. This takes place in an initialization phase. For this, initially, two measurements must be collected. In this context, a measurement is defined as the collection of vibration data over one single movement. Based on this data, the two spectrograms are calculated. To determine whether this reference quantity of two spectrograms captures the stochastic variation of the signal, the overall mean (Equation (2)) and standard deviation (Equation (3)) of the spectrograms are calculated.
In these formulas, k describes the number of measurements in the reference quantity. T is the time length of each measurement, F is the sampling frequency and spec(τ, ω) avg is the average value of spec(τ, ω) over measurements 0 to k. Afterwards, one measurement is added to the reference quantity at a time, and again avg spec,overall and std spec,overall are calculated. Plotting these standard deviations over the number of measurements in the reference quantity usually first shows an increase in std spec,overall and then a saturation as can be seen in Figure 2. If this saturation is reached, the reference quantity can sufficiently represent the stochastic behavior of the signal signature. In the shown example, this saturation is reached after 5 measurements.  After the initialization, an HI can be determined based on a newly collected measurement. For this, the measurement's spectrogram overall Z-score is determined according to Equation (4) .
In this context, spec(τ, ω) avg,re f and spec(τ, ω) std,re f are the mean value and the standard deviation of spec(τ, ω) for all measurements in the reference quantity. In Figure 3, the STFT and Z-score spectrograms of exemplary vibration measurements from a healthy and a faulty robot gear are depicted. The Z-score-based spectrogram of the faulty measurement shows more prominent changes compared to the STFT-based spectrogram.

Hi Evaluation Method
To compare the ability of the newly designed HI to cope with industrial robot application characteristics, we followed a three step approach. First of all, we investigated how well the designed HI can detect different kinds of faults in comparison to HIs from the state-of-the-art. Second, we investigated the temperature sensitivity of HIs from the stateof-the-art meeting this criterion and our HI. Third, we investigated the trend behavior of HIs showing a low temperature sensitivity on data from two accelerated wear tests. These three steps are now described more precisely. The overall process of our investigations is also described in Figure 4.

Varying Fault Detection Analysis
We used the FEMTO data set, which is described in detail in [35], to select HIs capable of detecting different faults. The data set is available in [36]. This data set provides run to failure vibration data from 16 identical bearings and for different faults and working conditions defined by the applied load and the rotational speed. The acceleration sensor sampled data with 25.6 kHz, one measurement has a length of 0.1 s and measurements were taken in equidistant timesteps of 10 s for all bearings. The test run for one bearing ended when the signal from the acceleration sensor exceeded 20 g. Therefore, different numbers of measurements are available per bearing ranging from 230 to 2803. We calculated the HIs summarized in Table 1 for all measurements of one sensor. These HIs were derived from several review papers regarding gearbox and bearing CM [37][38][39][40] and the publications mentioned in Section 1. Therefore, the HI calculation was based either on the raw acceleration signal, an enveloped signal as described in [41] or the residual signal as suggested by [4]. Additionally, the newly designed HI presented in Section 2 was calculated for the measurements based on the raw signals. To detect whether these HIs are sensitive to multiple faults, different techniques can be applied. In addition to filter techniques, ensemble, wrapper and embedded methods exist [42]. However, the latter three techniques combine classification or regression models with HIs for their evaluation. Hence, this evaluation is always dependent on the used models. Thus, we chose to use filter methods for the evaluation. Here, different figures of merit for regression and classification tasks can be applied, such as trendability, robustness, monotony or discriminance [42]. To combine these different performance indicators, we fitted different basic functions on the HIs calculated for the last 20 percent of measurements per bearing. These functions were first and second degree polynoms, exponential and sigmoid functions. For each of the fits, we calculated the R² value. This means that we received four R² values per HI and bearing. High R² values of these fits correlate with a high trendability, monotony, robustness and discriminance, which are desirable for HIs. To evaluate whether an HI can detect several damages, we considered only the best R² value per HI and bearing. We plotted the statistics of these 16 remaining R² values per HI as a boxplot. Suitable HIs should show high R² values with low variance.

Temperature Sensitivity Analysis
HIs showing this behavior were analyzed regarding their temperature sensitivity. For this purpose, we acquired vibration data from an industrial robot test rig. This test rig consists of a KUKA KR510 industrial robot with an attached load of 365 kg. We attached acceleration sensors close to the gearboxes as shown in Figure 5 on the right side. These sensors have a sampling rate of 26 kHz. The acceleration direction of the sensors was orthogonal to their contact area. For data acquisition, the robot performed a trajectory where each joint was moved individually at different speeds in an angle area of 10°, as described in Figure 6, and for different gear temperatures in the range of 25°C and 60°C and 5°C steps. One measurement per axis lasted 8 s. The gear temperature was measured at the gearbox cap with an infrared thermometer. For each temperature step, four measurements were made. For each measurement at each temperature step, the remaining HIs were calculated. To determine the temperature sensitivity, we divided the average HI values calculated from measurements at the highest gear temperatures by the values calculated from measurements at the lowest temperature. HIs with a high sensitivity were eliminated for the last step.

Accelerated Wear Test Analysis
Here, we calculated the remaining HIs for measurements from two data sets from accelerated robot wear tests to see how these HIs perform in a more industry like setting and how they cope with instationary velocity behavior. The first data set was collected during a time range of approximately one year with an ABB robot of type RB 6600-255/2.55. During the data acquisition, the robot performed an isolated movement of the second axis in an angle area of 150°for each measurement. Vibration data were only acquired with a sensor attached axially at the robot axis 2 gearbox. At the end of the experiment, the gearbox was dismantled and faults on the bearings and the shafts of the gear were found. A total of 2290 measurements, equally distributed over time, were taken for our analysis from this data set. One measurement lasted 1.6 s and the sampling rate was 10 kHz. More detailed information about this experiment can be found in [7,8]. The second data set was derived from another experiment. Here, the second axis of an ABB IRB 7600-340/2.8 was moved in an angle area of 80°continuously over the time frame of three months. The vibration sensor attached to the gearbox cap of axis 2 sampled with 20 kHz and one measurement lasted 2.15 s. The measurement setup is presented on the left side in Figure 5. The experiment ended after a roller element of a bearing had cracked and had blocked the gear. In this time range, 920 vibration measurements were taken in total in equidistant time steps. The faults, which occurred in both experiments, can be seen in Figure 7. In both experiments, environmental conditions such as load and trajectory were kept constant. Fluctuations of the temperature were kept at a minimum due to the constant movements of the robots. In this way, signal changes are likely to be correlated to increasing wear.

Results
This section is divided in three parts. First, the results from the varying fault detection experiments are shown. Secondly, the results from the temperature sensitivity analysis are presented. Finally, the application of the HIs on the two accelerated wear tests is described.

Varying Fault Detection Analysis
From the 16 bearing experiments, the HIs presented in Table 1 were calculated. We used the first 100 measurements per bearing as the reference quantity for the Z-score-HI and set R to 128. Figure 8 shows the R² values for a selection of different HIs as a box plot. The R² statistics for all HIs can be found in Appendix A. The abbreviations of the HIs are explained in Table 1. The PtP-, Peak-, RMS-, Std-and Z-score-HI show the highest R² values on average. They also show the lowest variance between the different bearings. This means that these HIs detect different faults most reliably. Other HIs show also high trend values but only for some of the bearings. HIs derived from the frequency-domain (DomF, SpC, SpE, SpF, SpRO) perform worse compared to HIs from the time-domain. The preprocessing steps of enveloping the signal or calculating the residual signal do not affect the HI trend behavior significantly, which can be seen in Tables A1-A3. The TDI-, and  DWTRMS-HI for specific frequency bands also show high average values with changing  variance (see Table A4). If these HIs would be used for robot gear condition monitoring, the progress of all frequency band specific HIs would have to be tracked as different faults stimulate changes in different frequency bands.

Temperature Sensitivity Analysis
Based on this result, we conducted the temperature sensitivity analysis for the PtP-, Peak-, RMS-, Std-, TDI-, DWTRMS-and Z-score-HI. Here, we used one measurement per temperature step as the reference quantity for the Z-score-HI and set R to 128. Figure 9 shows the change of the HIs per axis in percent for the PtP-, Peak-, RMS-, Std-and Z-score-HI. The RMS-and Z-score-HI show the lowest temperature sensitivity overall. Figure 10 shows the results for the DWTRMS-HIs. Here, high sensitivities for different detail coefficient DWTRMS-HIs exist. Figure A1 shows the temperature sensitivity of the TDI-HIs of different frequency bands. Here, a similar result can be seen compared to the DWTRMS-HIs. The data of Figures 10 and A1 can also be found in Tables A5 and A6. In general, the data from axis 4 show the highest temperature sensitivity for all HIs. The comparably higher sensitivity of the HI values derived from data at axis 4 can be related to the robot trajectory. During the trajectory, the robot arm was stretched out, which leads to greater elasticity at the position of the sensor at axis 4. This can cause increased vibrations, which are magnified under changing temperature influences. Given the results of the temperature sensitivity analysis, we analyzed the data sets from the accelerated wear tests with only the RMS-and the Z-score-HI. The other HIs were excluded due to their high temperature sensitivity. Even though some frequency band specific DWTRMS-HIs and TDI-HIs show low sensitivity, they were excluded as robot gear faults do not have to stimulate these frequency bands with low sensitivity.

Accelerated Wear Tests Analysis
In this analysis, we used the first 100 measurements as the reference quantity for the Z-score-HI and set R to 256. For smoothing, we applied a rolling average with a window length of 15 on both HI series. The progress of the HIs in the accelerated wear test of the ABB IRB 7600 is shown in Figure 11. Both HIs show a plateau with increased values at the end of the experiment. It can be assumed that, at this point in time, faults have already been present. Here, the increased HI values over a longer time period could have been used as a decision criterion for maintenance actions. The measurements at the very end show decreased values again. We assume that this decrease is correlated to a part of the bearing roller. In the end of the experiment, one of the roller elements showed a large pit. During the measurements showing the higher HI values this detached part of the roller element could have been still slightly fixed at the roller element and thus could have caused high vibration. After full detachment, this noise level decreased again. For the measurements before the plateau, the RMS-HI shows higher fluctuations compared to the Z-score-HI. For instance, the RMS-HI shows a first high peak around measurement 100. Such peaks could lead to false alarms in a condition monitoring scenario and should be avoided.
The progress of the HIs in the other accelerated wear test performed with the ABB IRB 6600 is shown in Figure 12. Here, the Z-score-HI shows a trending behavior and the RMS shows a stationary progress. Both HIs show a high increase during the last measurements. In this experiment, the trending behavior of the Z-score could have been a criterion to execute maintenance actions. This information is not present in the RMSprogress. Based on the fact that the Z-score showed a better trend behavior in the ABB IRB 6600 experiment and less noisy behavior in the ABB IRB 7600 experiment, we suggest the use of the Z-score-HI for the condition monitoring of robot gears.

Discussion
The discussion is divided in four parts. First, some remarks regarding our designed HI are given. Afterwards, three parts make up the Results subsections.
To derive the spectrograms required for the Z-score-HI, the length of the window function must be defined. High values for R result in a high frequency resolution and low values in a high time resolution. For the individual experiments, we chose window lengths that lead to a good compromise between time and frequency resolution by inspecting spectrograms created with different window lengths. We chose window lengths that lead to spectrograms appearing the least noisy in a visual inspection. In an industrial setting, an automated approach should be developed for this dependent on the robot's trajectory and the used sensor.
The motivation to use the FEMTO data set to investigate HI performance was to assess HIs' capability to detect multiple faults. Within a robot gearbox, which are mostly RV reducers, not only bearings but also the gear teeth can have faults. Such faults are not taken into account by our analysis explicitly. However, the bearing faults present in the FEMTO data set, e.g., pitting, are similar to typical gear teeth or shaft damage from a signal analysis point of view. Damage from all components modulate the acceleration signals at a specific frequency and its sidebands. Exactly this capability to track such changes in the signal was investigated in our analysis. There also exist HIs that track energy changes at the specific component fault frequencies. Such HIs were excluded from our analysis because expert knowledge about the geometric characteristics of the gears, e.g., the bearing diameters or the number of roller elements, is required to calculate these HIs. This expert knowledge is usually not available to industrial robot users. We also excluded HIs that could be derived automatically from machine learning models, such as autoencoders, as the physical interpretation of these HIs is difficult and hence a transferability between different robot systems is questionable from our point of view.
Regarding the results of the temperature sensitivity analysis, it must be pointed out that the results are valid only for the chosen robot trajectory. As the dynamic behavior of the robot changes within its working space, this analysis should be performed individually for trajectories and robot systems. However, from a theoretical point of view, the Z-score-HI possesses the ability to cope with these temperature fluctuations independently of the trajectory. Temperature variations lead to variance in the STFT spectrograms. This variance is taken into account in the spec(τ, ω) avg,re f and spec(τ, ω) std,re f during the initialization phase. Hence, Z-score-HIs derived from measurements from functional robot gears and different temperatures will show only little differences in the Z-score-HI value. This becomes more clear considering Figure 13. Here, the STFT and Z-score spectrograms from two vibration measurements of the temperature sensitivity experiment are shown. On the left side, the spectrograms from a cold gear measurement are depicted. On the right side, the spectrograms from a warm gear measurement are shown. Differences are visible in the STFT spectrograms around seconds 1 and 2. No differences are visible in the Z-score spectrograms. The scales of the STFT spectrograms reach from −5 to 0 and the scales of the Z-score spectrograms from 0 to 1.5. Hence, the relative changes of the STFT spectrograms are bigger compared to the Z-score spectrograms. In this example, the total relative change in energy in the STFT spectrogram is 9.15 percent, whereas the total relative change in the Z-score spectrogram is just 1.63 percent.
Finally, the results from the accelerated wear tests show noisy progress over time. This hinders a simple or automated detection of faults in a condition monitoring behavior. To establish an automated CM system, a suitable trend detection in combination with an outlier detection system must be set up. A trend detection system could identify HI progress shown as in Figure 12, whereas an outlier detection system could detect progress as depicted in Figure 11. The development of such a system also marks the outlook of our future work.

Conclusions
Condition monitoring of robot gears has the potential to decrease production system downtimes. The state-of-the-art provides many health indicators to track the health state of gears. We analyzed these health indicators regarding specific requirements rising from typical industrial robot applications. These requirements are the ability to detect different faults, low temperature sensitivity and the capability to deal with instationary velocity behavior. Additionally, we suggested a new health indicator based on STFT spectrograms and Z-scores that can cope with these requirements. Our analysis showed that the RMS health indicator and our suggested health indicator meet the defined requirements the best. Data from accelerated wear tests show that for an automatic condition monitoring system a combination of a trend detection and an outlier detection system that can deal with a noisy signal is required. Funding: We express our gratitude to the Bavarian Ministry of Economic Affairs, Regional Development, and Energy for the funding of our research. The formulated outlook will be investigated as part of the research project "KIVI" (grant number IUK-1809-0008 IUK597/003) and will be further developed and implemented.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality reasons.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.  TDI0  TDI2  TDI4  TDI6  TDI8  TDI10  TDI12  TDI14  TDI16  TDI18  TDI20  TDI22  TDI24  TDI26  TDI28  TDI30  TDI32  TDI34  TDI36  TDI38  TDI40  TDI42  TDI44  TDI46  TDI48  TDI50  TDI52  TDI54  TDI56  TDI58  TDI60  TDI62  TDI64   0