# Bearing Diagnostics of Hydro Power Plants Using Wavelet Packet Transform and a Hidden Markov Model with Orbit Curves

Authors: Gabriel Pino; José Roberto Ribas; Luciana Fernandes Guimarães
Journal: Shock and Vibration
 (2018)

---

## Abstract

The contribution of a medium-sized hydro power plant to the power grid can be either at base load or at peak load. When the latter is the most common operation mode, it increases the start and stop frequency, intensifying the hydro turbine components’ degradation, such as the guide bearings. This happens due to more frequent operation in transient states, which means being outside the service point of the machines’ nominal condition, consisting of speed, flow, and gross head. Such transient state operation increases the runner bearings’ mechanical vibration. The readings are acquired during the runner start-ups and filtered by a DC component mean value and a wavelet packet transform. The filtered series are used to estimate the relationship between the maximum orbit curve displacement and the accumulated operating hours. The estimated equation associated with the ISO 7919-5 vibration standards establishes the sojourn times of the degradation states, sufficient to obtain the transition probability distribution. Thereafter, a triangular probability function is used to determine the observation probability distribution in each state. Both matrices are inputs required by a hidden Markov model aiming to simulate the equipment deterioration process, given a sequence of maximum orbit curve displacements.


---

## Body

## 1. Introduction

Hydro power plants (HPPs) are generally expected to run either at base load or at peak load. The former condition means they should operate uninterruptedly, irrespective of the interconnected grid’s instantaneous demand variation, due to their lower marginal costs. Most electricity flowing in national interconnected grid (NIG) is generated by large HPPs. In addition, peak demand is usually supplied by thermal power plants (TPPs) and small HPPs, due to their higher marginal costs and lower energy efficiency. In theory, HPPs should only perform a normal shutdown during planned maintenance and for reservoir water level control. The continuous operation of large HPPs near the best efficiency point (BEP) demands low intermittency, which in turn depends on the reservoir gross head constancy, regular rainfall periods, and rigorous maintenance plans [1]. The national system operator (NSO) forces all Brazilian HPPs to start, shut down, or change their generated power to attend the grid’s instantaneous requirements. During a system demand variation, runners can assume different operation states: static, start-up, no-load, low-load, nominal-load, and shut down. Start-up is one of the most damaging states regarding runner operating condition [2, 3].Depending on the HPP’s operation regime, its power unit components can be subject to greater degradation. A stressful start-stop cycle makes the runner operate in transient state many times, increasing degeneration of mechanical elements. Off-design conditions, such as unit start-stop, switching between operation regimes, load rejection, and out-of-phase synchronization, can influence rotors’ transient operations [4, 5]. The aging generated by electrical, thermal, and environmental stress due to excessive starts and stops can be compared with the years of operation and the estimated annual equivalent cost of reduced remaining lifetime [6]. There are four main failure modes related to volatile HPP conditions: cavitation, erosion, material defects, and fatigue [7]. For instance, pressure fluctuation due to rotor-stator interaction and occurrence of vortex rope in the draft tube at partial load operation are hydrodynamic effects addressed by researchers [8].Orbit curve diagrams have been used in condition-based maintenance (CBM) studies. They are included among the instruments used for vibration measurement and signal processing techniques for condition monitoring of machine tools in manufacturing operations [9, 10]. For instance, orbit curves were used in the evaluation of the dynamic behavior of power units affected by steam-whirl instability, a common phenomenon in steam turbines caused by a rise of the steam flow and pressure [11, 12].It is important to predict degradation of HPP’s components before they reach failure limits. It is also valuable to know how the system operation and the degradation process are related. This paper focuses on one of the most important rotating machinery components, the turbine guide bearing [13]. During start-ups, they are submitted to hydraulic transient state [14]. This phenomenon is associated with mechanical vibrations that lead to a higher shaft surface displacement and a likely guide bearing trimming process. Continuous start-stop cycles cause progressive degradation of rotating elements. Shaft vibration measurement is essential to understand and quantify the entire process, so vibration monitoring systems are commonly installed in power units’ bearings. These systems usually have limited resources regarding digital-to-analogue sampling frequency, memory buffer register and offset error induction. This lowers the quality measurements due to the presence of noise in a significant part of the information. More details are given in Section 3.2 and described in [15–18].The acquired transient vibration waveforms are nonstationary and present low signal-to-noise ratio (SNR). These original noisy signals are unsuitable to analyze and design orbit curves, mainly for two reasons: high DC levels distort the orbit curve formation from the origin of the Cartesian plane; and low SNR signals do not properly represent the axis displacement due to vibration sparks and undesired harmonic content. Section3.1 shows the difference of applying filtering techniques to better interpret the transient vibration behavior. The nonstationary waveforms require sophisticated filtering tools to handle time and frequency dependent data [19, 20]. There are several mathematical methods to deal with vibration signals. They have chronologically evolved from frequency-domain analysis regarding digital filters [21] to wavelet transform denoising [22]. The wavelet transform method, applied in this paper, has the mathematical mechanisms to properly denoise and process the primary data acquired.The filtered vibration database is accurate to develop a statistical method to predict and track the bearing degradation process. The signal processing applied to the database results in higher SNR and better represents the axis displacements as described in Section4.1. The stochastic phenomenon evolves from “good as new” to “imminent failure.” The state progression moves through adjacent states, in such way that the next state only depends on the current one, characterizing a Markov process. Markovian models are suitable for dividing equipment conditions into any meaningful state, facilitating general understanding of the specification of the process model [23]. Furthermore, they can characterize the stochastic relationship between the features extracted from condition monitoring data and the actual health states of the equipment [24].Nevertheless, most of the time the bearing degradation state is not detectable because it is not possible to perform inspections during the HPP operation. So, the degradation process cannot be followed directly, meaning evolution is hidden. To overcome this limitation, vibration measurements are used as an observable variable indicating the degradation states of the bearing. The entire process can be statistically better represented as a hidden Markov model (HMM) [25–28]. Besides the proposed purpose, HMM has other successful applications in speech recognition [29–31], monitoring of machine conditions, and fault diagnosis [32, 33].This paper develops a procedure that does not require extensive historical vibration data to perform statistical assessment. The model also incorporates previous technical experiences acquired during preventive and corrective maintenance performed at the power plants and these parameters determine the decision model for maintenance based on the condition-based maintenance (CBM) of the equipment. It is possible to ascertain the guide bearing deterioration condition, for instance, based on available measurements and the proportional vibration limits established in ISO 7919-5 [34, 35]. After the theoretical background presented in Section 2 and the technical development in Section 3 regarding vibration and stochastic methods, Section 4 discusses and applies the proposed method to measurements obtained from generating unit 1 (GU-1) of Corumbá IV HPP, located in Brazil, and finally Section 5 presents the conclusion.


## 2. Theoretical Background

The use of techniques for prognosis of bearing degradation should be divided into two parts: data acquisition and the data processing [36]. When considering data acquisition, this study performs a nonstationary signal analysis for bearing fault diagnosis. Feng et al. [37] classified the various time-frequency analysis methods applied to machinery fault diagnosis into four categories: linear and bilinear time-frequency representations (e.g., wavelet transform); adaptive parametric time-frequency analysis; adaptive nonparametric time-frequency analysis (e.g., Hilbert-Huang transform) [38]; and time varying higher order spectra. Various methods adapted to specific conditions have been used since then. The nuisance attribute projection (NAP) is introduced in the bearing performance degradation assessment to mitigate the influence of problems irrelevant to the degradation state introduced by operation conditions [25]. The Rényi entropy based features exploit the idea that a progressing fault implies rising dissimilarity in the distribution of energies across the vibrational spectral band sensitive to the bearing faults [39]. Envelope analysis has been applied to overcome the constraint of constant operating speed of the rolling element bearing. The squared envelope spectrum has been extended to cases in which small speed fluctuations occur [40, 41]. The spectral kurtosis technique adopts the concept of kurtosis to capture the impulsiveness of a signal. It uses a combination of short-time Fourier transform-based SK, kurtogram, adaptive SK, and protrugram [42]. However, despite such specific situations, the wavelet transform has been the most popular denoising technique for the extraction of the defect vibratory signature from the measured signal in which the random noise and other parameters of the bearing are immersed [43–45]. With respect to data processing, the stochastic process inherent to bearing wear can be classified as a single component which depends on the nature of the degradation state: discrete or continuous [46]. When the technical condition of a bearing cannot be described by a continuous, measurable value, its deterioration is often assessed by visual inspections or by other methods that only lead to qualitative results. As an alternative to continuous measures, several discrete deterioration states can be defined [47]. Markov processes are used in such cases when measurement of discrete degradation states is not precise.


## 3. Materials and Methods

The reasons why a guide bearing degradation increases are strongly related to the shaft’s axial displacement. This leads to possible relationship between the shaft vibration intensity and the deterioration caused by mechanical contact of rotating and static elements. Vibration measurements and degradation status are so closely related that many power plants have installed vibration monitoring systems to detect its evolution [48]. This type of system records a substantial increase of vibration level during hydro power units’ start-up, which takes approximately five minutes, depending on speed governor parameterization [49]. The vibration data acquired during these events are crucial to understand its consequences to the whole rotating group, including the evolution of the guide bearings’ degradation.


### 3.1. Data Acquisition

The power unit shaft’s mechanical vibration usually presents two signal sources, orthogonally installed, as shown in Figure1. The configuration of both signals permits parameterizing them so as to describe the external profile of the shaft. The graph of the signal parameterization and composition that suppresses the time dependency is generally called a Lissajous curve [50], and specifically for shaft surface displacement it is called an orbit curve [51, 52]. This curve is fundamental to infer how much the shaft perimeter is displaced from its center line. This particular measurement, called Smax, is the maximum axis vector distance and its measurement is also used as an indicator for technical standards [34]. It is calculated as in (1)Smax=max⁡xt2+yt2,0≤t≤1nrpm,where nrpm is the rated turbine speed in revolutions per minute. This article considers Smax a collusion variable of the entire vibration phenomenon per machine revolution in transient and steady states operation, so it is stored in a database for stochastic analysis.Figure 1
Sensors displacement and measurements.Vibration monitoring systems installed in power plants usually use inductive proximity sensors to measure the distance between probe tips and the shaft surface [53]. Inductive sensors present low-pass filter response for displacement measurements; thus they physically attenuate high-frequency signals that could affect the vibration waveform analysis [54]. Besides this, most commercial electronic drives responsible for signal processing and analogue-to-digital converting operate with maximum frequency bandwidth up to 20 kHz [55, 56]. Another technical limitation of commercial systems is the memory buffers required to acquire long waveforms during five minutes and store the data. Lower frequency sampling implies worse information quality regarding vibration phenomenon due to the inherent presence of flicker noise [57]. Associated with this noise, HPP turbines usually present low rated speed, which implies their first harmonic is up to 10 Hz, meaning that the fundamental signal frequency component is polluted by flicker and distributed thermal noise [58].Technical issues related to data acquisition require signal processing using filtering techniques for better understanding of transient and steady state vibration waveforms.


### 3.2. Data Filtering

The technical limitations of vibration monitoring systems directly impact the precision and accuracy of orbit curve calculation. Hence, proper filtering tools must be applied to acquired waveforms. The signal processing depends on the correct mathematical approach given to each data characteristic, operation state, sensor accuracy, associated errors of analogue-to-digital converters, frequency sampling, SNR, noise distribution, and window lengths. Figure2 illustrates the orbit curves formation regarding three different inputs: original vibration signals from the monitoring system, signals after DC filtering, and signals after DC and AC filtering. The first image does not properly represent the axis displacement and its center line is shifted from the origin. The second one is aligned to the center line but is still very noisy. The last one, after signal processing of the vibration inputs, better represents the axis displacement.Figure 2
Orbit curve formation for different inputs quality.


#### 3.2.1. Frequency Spectrum

Frequency spectrum analysis is applied before choosing filtering tools in order to better understand the signal behavior and how its processing will be properly performed. This preliminary mathematical feature indicates the signal content and which harmonic components are the most important [59].Specifically, the transient vibration signals acquired during machine start-ups present a nonstationary characteristic; that is, the frequency spectrum is not constant in the whole window analysis. This condition demands more general time and frequency-domain analysis such as short-time Fourier transform (STFT) [60]. Equation (2) shows the integral kernel dependency on time and frequency and (3) the Haning window function with period τ.(2)SFxt=Xf,τ=∫-∞+∞xte-j2πftWtdt,(3)Wt=121-cos⁡2πtτ.Using the time and frequency excursion of the STFT it is possible to infer a diverse frequency content of the AC component and also a DC level. The DC level content indicates the sensor setup before the data acquisition. The intrinsic differences between the original AC and DC vibration components require different mathematical tools to properly filter them.


#### 3.2.2. DC Component Filtering

Orbit curve calculation necessarily uses peak-to-peak vibration values. Otherwise they would be shifted from the origin of the Cartesian plane and the waveform would be distorted, as can be inferred by (1) and seen in Figure 1. So the DC component must be filtered. The spectrum calculated by (2) possibly indicates a DC content during transient and steady states, which can be interpreted as the stand still distance between the shaft surface and the inductive sensor probe [55].There are many techniques that can filter the DC component such as digital FIR or IIR high pass filtering [61]. However, these techniques present amplitude and phase frequency response. The behavior of a very low cutoff frequency filter implies nonlinear phase response near the inflection, which is not desirable since the fundamental frequency is usually up to 10 Hz. In order to avoid phase shifting, this paper uses the mean vibration value as the DC reference.


#### 3.2.3. AC Component Filtering

The vibration measurement during the runner start-up tends to assume nonstationary behavior, as can also be inferred using (2). This type of signal is not properly filtered by applying conventional techniques such as impulse response filtering. Wavelet transform is a very efficient denoising tool able to adapt the original waveform characteristics by means of frequency and amplitude modulation [62]. Equation (4) denotes the wavelet transform of a continuous signal x(t). Although the acquired signals are discrete, the continuous expression better exemplifies the advantages of using the wavelet packet.(4)Wxt=Xa,b=1a∫-∞+∞xtψt-badt,where ψ=ψ(t) is the mother wavelet function, a is the frequency dilatation parameter and b is the shifting parameter related to signal excursion in time-domain. The kernel composition of time and frequency parameters enables the integral transform to adapt them to the original signal characteristics, which is desirable for nonstationary measurements.After the signal processing, the vibration data and the shaft replacement are ready to be statistically analyzed.


### 3.3. Degradation Prediction Model

#### 3.3.1. Theoretical Limits forSmax

According to ISO 7919-5, there is a relationship betweenSmax at the nominal machine operation and the guide bearing deterioration stage, depending on machine nominal speed. The vibration magnitudes can be related to four zones, which are defined as follows. The vibration of newly commissioned machines would normally fall within zone A. In zone B, machines are usually considered acceptable for unrestricted long-term operation. Machines with vibration within zone C are ordinarily considered unsatisfactory for long-term continuous operation. Generally, the machine can be operated for a limited period in this condition until a suitable opportunity arises for maintenance. In zone D, vibration values are normally considered to be of sufficient severity to cause damage to the runner. These limits define the degradation states and the observable vibration limits in the Hidden Markov Model.


#### 3.3.2. Triangular Probability Function

The association of a discrete nominal scale represented by linguistic expressions like “small deterioration” with a continuous degradation process is a fertile field for errors of judgment. SinceSmax is a discrete observed variable associated with each state emission, it imbeds an imprecision related to both the nominal classification and the quality of the observed variable acquisition and filtering. A simple way to mitigate such imprecision is to work with intervals, in which case a symmetric triangular probability distribution is proposed. In this representation, the observed variable is defined by a set of three numbers (L, M, U). These numbers are determined by M, a discrete value reflecting the calculated Smax, and the interval, a sensitive value used to obtain the lower (L) and upper (U) bounds. The higher the degree, the higher the perceived imprecision regarding the observed variable. The triangular-density function is customarily used when function parameters are subjective. A typical example is for estimating activity durations by identifying a minimum, most likely, and maximum values, representing uncertain activity durations [63]. It has been shown that the cutoff points of such triangular function contain the confidence intervals of any symmetric probability distribution with the same mode and support [64].


#### 3.3.3. Hidden Markov Model

Markov chains are stochastic processes in which the progression to the next state depends exclusively on the previous one [29]. The Hidden Markov Model (HMM) is a Markov progression generally observed indirectly, as a function of the transitions among states of the Markov chain governing process, and in which successive state transitions are determined by the likelihood model transition [29, 31–33, 65–68]. Since inspections are not carried out continuously, unlike continuous and progressive equipment deterioration, the states cannot be observed, and as a result a Hidden Markov Model applies. Figure 3 demonstrates a hidden Markov process. In this case, one or more states elements S1,S2,…,SM represent the degradation states and they are not directly observable. However, they emit observations O1,O2,…,OM that can be identified and associated with the states.Figure 3
Emitted and observable states in HMM.The elementsaij refer to the transition probabilities between states. For instance, a11 is the probability that the process at time T+1 remains in state S1 knowing that at time T it was in state S1. For machine degradation, the transitions between states do not allow the return to the previous state or the passing to the next state. This means the degradation of the equipment from one moment T to the next T+1 will only stay in the home state or it will switch to the next. Thus, for three states, the transition matrix assumes (5)T=P11P1200P22P23001.It can be observed thatP33=1; that is, once the last degradation state is reached, a potential failure condition occurs and corrective maintenance of the equipment is necessary. So the process is reset and the state is reestablished to “good as new.” Based on the sojourn times, it is possible to estimate the transition matrix by the Bayesian method. The elements bij show the emission probabilities. These are the probabilities with which the signals of the observable sequence are produced in each of the states. For example, state S1 can emit several values for the observable signal: the probability of issuing O1 is represented by b11, of issuing O2 by b22 and so on. In addition, the model has as initial distribution of states πi, where πi=P(S1=i) for 1≤i≤M. The HMM process provides ways to perform the task of learning the parameters; that is, given a sequence of results, the algorithms estimate by maximum likelihood the transition and emission probabilities that best represent them. Figure 4 summarizes the proposed methodology.Figure 4
Methodology flowchart.


## 3.1. Data Acquisition

The power unit shaft’s mechanical vibration usually presents two signal sources, orthogonally installed, as shown in Figure1. The configuration of both signals permits parameterizing them so as to describe the external profile of the shaft. The graph of the signal parameterization and composition that suppresses the time dependency is generally called a Lissajous curve [50], and specifically for shaft surface displacement it is called an orbit curve [51, 52]. This curve is fundamental to infer how much the shaft perimeter is displaced from its center line. This particular measurement, called Smax, is the maximum axis vector distance and its measurement is also used as an indicator for technical standards [34]. It is calculated as in (1)Smax=max⁡xt2+yt2,0≤t≤1nrpm,where nrpm is the rated turbine speed in revolutions per minute. This article considers Smax a collusion variable of the entire vibration phenomenon per machine revolution in transient and steady states operation, so it is stored in a database for stochastic analysis.Figure 1
Sensors displacement and measurements.Vibration monitoring systems installed in power plants usually use inductive proximity sensors to measure the distance between probe tips and the shaft surface [53]. Inductive sensors present low-pass filter response for displacement measurements; thus they physically attenuate high-frequency signals that could affect the vibration waveform analysis [54]. Besides this, most commercial electronic drives responsible for signal processing and analogue-to-digital converting operate with maximum frequency bandwidth up to 20 kHz [55, 56]. Another technical limitation of commercial systems is the memory buffers required to acquire long waveforms during five minutes and store the data. Lower frequency sampling implies worse information quality regarding vibration phenomenon due to the inherent presence of flicker noise [57]. Associated with this noise, HPP turbines usually present low rated speed, which implies their first harmonic is up to 10 Hz, meaning that the fundamental signal frequency component is polluted by flicker and distributed thermal noise [58].Technical issues related to data acquisition require signal processing using filtering techniques for better understanding of transient and steady state vibration waveforms.


## 3.2. Data Filtering

The technical limitations of vibration monitoring systems directly impact the precision and accuracy of orbit curve calculation. Hence, proper filtering tools must be applied to acquired waveforms. The signal processing depends on the correct mathematical approach given to each data characteristic, operation state, sensor accuracy, associated errors of analogue-to-digital converters, frequency sampling, SNR, noise distribution, and window lengths. Figure2 illustrates the orbit curves formation regarding three different inputs: original vibration signals from the monitoring system, signals after DC filtering, and signals after DC and AC filtering. The first image does not properly represent the axis displacement and its center line is shifted from the origin. The second one is aligned to the center line but is still very noisy. The last one, after signal processing of the vibration inputs, better represents the axis displacement.Figure 2
Orbit curve formation for different inputs quality.


### 3.2.1. Frequency Spectrum

Frequency spectrum analysis is applied before choosing filtering tools in order to better understand the signal behavior and how its processing will be properly performed. This preliminary mathematical feature indicates the signal content and which harmonic components are the most important [59].Specifically, the transient vibration signals acquired during machine start-ups present a nonstationary characteristic; that is, the frequency spectrum is not constant in the whole window analysis. This condition demands more general time and frequency-domain analysis such as short-time Fourier transform (STFT) [60]. Equation (2) shows the integral kernel dependency on time and frequency and (3) the Haning window function with period τ.(2)SFxt=Xf,τ=∫-∞+∞xte-j2πftWtdt,(3)Wt=121-cos⁡2πtτ.Using the time and frequency excursion of the STFT it is possible to infer a diverse frequency content of the AC component and also a DC level. The DC level content indicates the sensor setup before the data acquisition. The intrinsic differences between the original AC and DC vibration components require different mathematical tools to properly filter them.


### 3.2.2. DC Component Filtering

Orbit curve calculation necessarily uses peak-to-peak vibration values. Otherwise they would be shifted from the origin of the Cartesian plane and the waveform would be distorted, as can be inferred by (1) and seen in Figure 1. So the DC component must be filtered. The spectrum calculated by (2) possibly indicates a DC content during transient and steady states, which can be interpreted as the stand still distance between the shaft surface and the inductive sensor probe [55].There are many techniques that can filter the DC component such as digital FIR or IIR high pass filtering [61]. However, these techniques present amplitude and phase frequency response. The behavior of a very low cutoff frequency filter implies nonlinear phase response near the inflection, which is not desirable since the fundamental frequency is usually up to 10 Hz. In order to avoid phase shifting, this paper uses the mean vibration value as the DC reference.


### 3.2.3. AC Component Filtering

The vibration measurement during the runner start-up tends to assume nonstationary behavior, as can also be inferred using (2). This type of signal is not properly filtered by applying conventional techniques such as impulse response filtering. Wavelet transform is a very efficient denoising tool able to adapt the original waveform characteristics by means of frequency and amplitude modulation [62]. Equation (4) denotes the wavelet transform of a continuous signal x(t). Although the acquired signals are discrete, the continuous expression better exemplifies the advantages of using the wavelet packet.(4)Wxt=Xa,b=1a∫-∞+∞xtψt-badt,where ψ=ψ(t) is the mother wavelet function, a is the frequency dilatation parameter and b is the shifting parameter related to signal excursion in time-domain. The kernel composition of time and frequency parameters enables the integral transform to adapt them to the original signal characteristics, which is desirable for nonstationary measurements.After the signal processing, the vibration data and the shaft replacement are ready to be statistically analyzed.


## 3.2.1. Frequency Spectrum

Frequency spectrum analysis is applied before choosing filtering tools in order to better understand the signal behavior and how its processing will be properly performed. This preliminary mathematical feature indicates the signal content and which harmonic components are the most important [59].Specifically, the transient vibration signals acquired during machine start-ups present a nonstationary characteristic; that is, the frequency spectrum is not constant in the whole window analysis. This condition demands more general time and frequency-domain analysis such as short-time Fourier transform (STFT) [60]. Equation (2) shows the integral kernel dependency on time and frequency and (3) the Haning window function with period τ.(2)SFxt=Xf,τ=∫-∞+∞xte-j2πftWtdt,(3)Wt=121-cos⁡2πtτ.Using the time and frequency excursion of the STFT it is possible to infer a diverse frequency content of the AC component and also a DC level. The DC level content indicates the sensor setup before the data acquisition. The intrinsic differences between the original AC and DC vibration components require different mathematical tools to properly filter them.


## 3.2.2. DC Component Filtering

Orbit curve calculation necessarily uses peak-to-peak vibration values. Otherwise they would be shifted from the origin of the Cartesian plane and the waveform would be distorted, as can be inferred by (1) and seen in Figure 1. So the DC component must be filtered. The spectrum calculated by (2) possibly indicates a DC content during transient and steady states, which can be interpreted as the stand still distance between the shaft surface and the inductive sensor probe [55].There are many techniques that can filter the DC component such as digital FIR or IIR high pass filtering [61]. However, these techniques present amplitude and phase frequency response. The behavior of a very low cutoff frequency filter implies nonlinear phase response near the inflection, which is not desirable since the fundamental frequency is usually up to 10 Hz. In order to avoid phase shifting, this paper uses the mean vibration value as the DC reference.


## 3.2.3. AC Component Filtering

The vibration measurement during the runner start-up tends to assume nonstationary behavior, as can also be inferred using (2). This type of signal is not properly filtered by applying conventional techniques such as impulse response filtering. Wavelet transform is a very efficient denoising tool able to adapt the original waveform characteristics by means of frequency and amplitude modulation [62]. Equation (4) denotes the wavelet transform of a continuous signal x(t). Although the acquired signals are discrete, the continuous expression better exemplifies the advantages of using the wavelet packet.(4)Wxt=Xa,b=1a∫-∞+∞xtψt-badt,where ψ=ψ(t) is the mother wavelet function, a is the frequency dilatation parameter and b is the shifting parameter related to signal excursion in time-domain. The kernel composition of time and frequency parameters enables the integral transform to adapt them to the original signal characteristics, which is desirable for nonstationary measurements.After the signal processing, the vibration data and the shaft replacement are ready to be statistically analyzed.


## 3.3. Degradation Prediction Model

### 3.3.1. Theoretical Limits forSmax

According to ISO 7919-5, there is a relationship betweenSmax at the nominal machine operation and the guide bearing deterioration stage, depending on machine nominal speed. The vibration magnitudes can be related to four zones, which are defined as follows. The vibration of newly commissioned machines would normally fall within zone A. In zone B, machines are usually considered acceptable for unrestricted long-term operation. Machines with vibration within zone C are ordinarily considered unsatisfactory for long-term continuous operation. Generally, the machine can be operated for a limited period in this condition until a suitable opportunity arises for maintenance. In zone D, vibration values are normally considered to be of sufficient severity to cause damage to the runner. These limits define the degradation states and the observable vibration limits in the Hidden Markov Model.


### 3.3.2. Triangular Probability Function

The association of a discrete nominal scale represented by linguistic expressions like “small deterioration” with a continuous degradation process is a fertile field for errors of judgment. SinceSmax is a discrete observed variable associated with each state emission, it imbeds an imprecision related to both the nominal classification and the quality of the observed variable acquisition and filtering. A simple way to mitigate such imprecision is to work with intervals, in which case a symmetric triangular probability distribution is proposed. In this representation, the observed variable is defined by a set of three numbers (L, M, U). These numbers are determined by M, a discrete value reflecting the calculated Smax, and the interval, a sensitive value used to obtain the lower (L) and upper (U) bounds. The higher the degree, the higher the perceived imprecision regarding the observed variable. The triangular-density function is customarily used when function parameters are subjective. A typical example is for estimating activity durations by identifying a minimum, most likely, and maximum values, representing uncertain activity durations [63]. It has been shown that the cutoff points of such triangular function contain the confidence intervals of any symmetric probability distribution with the same mode and support [64].


### 3.3.3. Hidden Markov Model

Markov chains are stochastic processes in which the progression to the next state depends exclusively on the previous one [29]. The Hidden Markov Model (HMM) is a Markov progression generally observed indirectly, as a function of the transitions among states of the Markov chain governing process, and in which successive state transitions are determined by the likelihood model transition [29, 31–33, 65–68]. Since inspections are not carried out continuously, unlike continuous and progressive equipment deterioration, the states cannot be observed, and as a result a Hidden Markov Model applies. Figure 3 demonstrates a hidden Markov process. In this case, one or more states elements S1,S2,…,SM represent the degradation states and they are not directly observable. However, they emit observations O1,O2,…,OM that can be identified and associated with the states.Figure 3
Emitted and observable states in HMM.The elementsaij refer to the transition probabilities between states. For instance, a11 is the probability that the process at time T+1 remains in state S1 knowing that at time T it was in state S1. For machine degradation, the transitions between states do not allow the return to the previous state or the passing to the next state. This means the degradation of the equipment from one moment T to the next T+1 will only stay in the home state or it will switch to the next. Thus, for three states, the transition matrix assumes (5)T=P11P1200P22P23001.It can be observed thatP33=1; that is, once the last degradation state is reached, a potential failure condition occurs and corrective maintenance of the equipment is necessary. So the process is reset and the state is reestablished to “good as new.” Based on the sojourn times, it is possible to estimate the transition matrix by the Bayesian method. The elements bij show the emission probabilities. These are the probabilities with which the signals of the observable sequence are produced in each of the states. For example, state S1 can emit several values for the observable signal: the probability of issuing O1 is represented by b11, of issuing O2 by b22 and so on. In addition, the model has as initial distribution of states πi, where πi=P(S1=i) for 1≤i≤M. The HMM process provides ways to perform the task of learning the parameters; that is, given a sequence of results, the algorithms estimate by maximum likelihood the transition and emission probabilities that best represent them. Figure 4 summarizes the proposed methodology.Figure 4
Methodology flowchart.


## 3.3.1. Theoretical Limits forSmax

According to ISO 7919-5, there is a relationship betweenSmax at the nominal machine operation and the guide bearing deterioration stage, depending on machine nominal speed. The vibration magnitudes can be related to four zones, which are defined as follows. The vibration of newly commissioned machines would normally fall within zone A. In zone B, machines are usually considered acceptable for unrestricted long-term operation. Machines with vibration within zone C are ordinarily considered unsatisfactory for long-term continuous operation. Generally, the machine can be operated for a limited period in this condition until a suitable opportunity arises for maintenance. In zone D, vibration values are normally considered to be of sufficient severity to cause damage to the runner. These limits define the degradation states and the observable vibration limits in the Hidden Markov Model.


## 3.3.2. Triangular Probability Function

The association of a discrete nominal scale represented by linguistic expressions like “small deterioration” with a continuous degradation process is a fertile field for errors of judgment. SinceSmax is a discrete observed variable associated with each state emission, it imbeds an imprecision related to both the nominal classification and the quality of the observed variable acquisition and filtering. A simple way to mitigate such imprecision is to work with intervals, in which case a symmetric triangular probability distribution is proposed. In this representation, the observed variable is defined by a set of three numbers (L, M, U). These numbers are determined by M, a discrete value reflecting the calculated Smax, and the interval, a sensitive value used to obtain the lower (L) and upper (U) bounds. The higher the degree, the higher the perceived imprecision regarding the observed variable. The triangular-density function is customarily used when function parameters are subjective. A typical example is for estimating activity durations by identifying a minimum, most likely, and maximum values, representing uncertain activity durations [63]. It has been shown that the cutoff points of such triangular function contain the confidence intervals of any symmetric probability distribution with the same mode and support [64].


## 3.3.3. Hidden Markov Model

Markov chains are stochastic processes in which the progression to the next state depends exclusively on the previous one [29]. The Hidden Markov Model (HMM) is a Markov progression generally observed indirectly, as a function of the transitions among states of the Markov chain governing process, and in which successive state transitions are determined by the likelihood model transition [29, 31–33, 65–68]. Since inspections are not carried out continuously, unlike continuous and progressive equipment deterioration, the states cannot be observed, and as a result a Hidden Markov Model applies. Figure 3 demonstrates a hidden Markov process. In this case, one or more states elements S1,S2,…,SM represent the degradation states and they are not directly observable. However, they emit observations O1,O2,…,OM that can be identified and associated with the states.Figure 3
Emitted and observable states in HMM.The elementsaij refer to the transition probabilities between states. For instance, a11 is the probability that the process at time T+1 remains in state S1 knowing that at time T it was in state S1. For machine degradation, the transitions between states do not allow the return to the previous state or the passing to the next state. This means the degradation of the equipment from one moment T to the next T+1 will only stay in the home state or it will switch to the next. Thus, for three states, the transition matrix assumes (5)T=P11P1200P22P23001.It can be observed thatP33=1; that is, once the last degradation state is reached, a potential failure condition occurs and corrective maintenance of the equipment is necessary. So the process is reset and the state is reestablished to “good as new.” Based on the sojourn times, it is possible to estimate the transition matrix by the Bayesian method. The elements bij show the emission probabilities. These are the probabilities with which the signals of the observable sequence are produced in each of the states. For example, state S1 can emit several values for the observable signal: the probability of issuing O1 is represented by b11, of issuing O2 by b22 and so on. In addition, the model has as initial distribution of states πi, where πi=P(S1=i) for 1≤i≤M. The HMM process provides ways to perform the task of learning the parameters; that is, given a sequence of results, the algorithms estimate by maximum likelihood the transition and emission probabilities that best represent them. Figure 4 summarizes the proposed methodology.Figure 4
Methodology flowchart.


## 4. Results and Discussion

The method presented in the previous section is applied as a case study to a power unit of Corumbá IV, which is a medium-sized hydro power plant located in Brazil. Its main technical characteristics are reported in Table1. The runner of Corumbá IV was chosen for the case study due to previous technical issues regarding the accelerated guide bearing degeneration process, which required corrective maintenance on August 16, 2015.Table 1
Corumbá IV HPP main technical characteristics.
Number of power units 2 Rated turbine power 64.8 MW Turbine diameter 3287 mm Rated turbine speed 200 RPM Net head 65 m Reservoir capacity 3.8 billion m3From August 2016 to May 2017, 14 measurements were acquired by the vibration monitoring system during the runner start-up. The method depicted in Figure4 is applied to these acquisitions.


### 4.1. Data Processing

It is strongly recommended to acquire the start-up vibration measurement from stand still until the no-load speed condition, to assure complete acquisition of vibration data. The process takes approximately five minutes. According to the technical data sheet of the vibration monitoring system [55], the maximum sampling rate is 128 Hz regarding a five-minute memory buffer. Obeying the Nyquist criteria, the information is completely determined up to 64 Hz, as can be inferred in the next subsection.


#### 4.1.1. Frequency Spectrum

Applying (2), for a 2 Hz Haning window, the result is pictured in Figure 5. The surface relates the normalized vibration power density V=V(f,t) to the excursion in time t and frequency f, given in minutes and hertz, respectively. The vibration function domain is DV={(f,t)∈R2∣0≤f≤64,0≤t≤5}.Figure 5
Short term Fourier transform of a transient state vibration.The normalized absolute value given by its image presents clear nonstationary behavior due to its nonconstant frequency arrangement. The hill over low frequencies presents decay up to 10 Hz, suggesting an intense flicker noise content due to the signal’s low sampling frequency [57]. Besides this, there are some random noise clouds over various frequencies of the AC component, even over the fundamental harmonic, which indicate distributed thermal noise pollution. There is also a constant DC level even after steady state is achieved, where the harmonic frequencies all act in the same direction, in approximately 250 s.


#### 4.1.2. DC Component Filtering

The DC filtering techniques discussed in Section3.2.2 are compared using a FIR filter and the mean value. Table 2 shows the amplitude of the DC components filtered by each method and their relative difference. The comparison infers that their relative differences are not substantial, so the mean value technique is adopted as an advantage, by not distorting the phase frequency response.Table 2
DC filtering techniques result.
Acquisition Mean value FIR filter Difference 16/04/2017 1089.76 [μm] 1087.34 [μm] 0.22 [%] 02/05/2017 1023.05 [μm] 1022.54 [μm] 0.05 [%] 20/03/2017 1018.53 [μm] 1018.09 [μm] 0.04 [%] 23/05/2017 1035.73 [μm] 1036.67 [μm] 0.09 [%] 09/09/2016 1063.90 [μm] 1063.50 [μm] 0.04 [%] 17/04/2017 1064.40 [μm] 1063.52 [μm] 0.08 [%]


#### 4.1.3. AC Component Filtering

Figure5 shows random noise over the entire surface domain of the unfiltered signal. The prominent flicker noise is more concentrated in regions up to 10 Hz and the thermal noise power density is well distributed in the whole frequency excursion. The wavelet packet used for denoising presents the Symlet 8 mother wavelet due to its good performance, with quasisinusoidal waveforms, and the denoising process considers soft thresholding at seven frequency levels. Table 3 shows a SNR comparison before and after the signal processing for some frequency harmonics [69, 70]. It can be noted that the signal processing considerably increased the relative signal power component over the noisy one. Figure 6 compares the original and filtered signals. Zooming reveals a large difference between each waveform, thus good adherence and filter behavior.Table 3
SNR comparison.
Harmonic Frequency [Hz] SNR before SNR after 0 0 37.39 317.62 1 3.33 35.68 304.38 3 10 25.99 205.45 6 20 17.77 191.3 8 26.66 15.16 132.9Figure 6
Comparison between original and filtered signals.


#### 4.1.4. Sensitivity Analysis betweenSmax and Operation Time

The maximum displacement evolution over the operating time fits an exponential function nicely, as depicted in Figure7. In a smoothing process, when the mean extracted from a subset containing the k greatest numbers of a sample is its statistic, the exponential adherence among the statistics of the 14 samples is validated. Such confirmation has an important consequence on practical grounds, since the exponential relation is preserved even when two or more greater numbers are considered, mainly from the presence of outliers. This paper considers the relationship function for the maximum absolute value, k=1, over operation time, as observed in Table 4.Table 4
S m a x values over operation time.
Acquisition k = 1 k = 10 k = 20 09/09/2016 149.92 147.71 146.59 19/12/2016 166.45 163.65 158.91 06/01/2017 169.22 161.93 156.87 19/01/2017 171.97 170.36 166.24 20/02/2017 178.09 175.69 170.64 21/02/2017 176.28 174.07 170.95 20/03/2017 183.35 180.22 174.52 05/04/2017 186.39 183.62 179.97 16/04/2017 188.74 181.79 177.00 17/04/2017 185.51 182.72 177.15 24/04/2017 190.84 183.91 178.70 02/05/2017 191.62 183.55 176.31 23/05/2017 196.98 191.45 184.93 30/05/2017 197.68 189.71 180.82Figure 7
Exponential adherence betweenSmax and operation time.


### 4.2. Degradation Model

#### 4.2.1. Degradation Model Inputs

The runner nominal speed is 200 RPM, and according to [34], the Smax theoretical limits for a 200 RPM machine are shown in Table 5. In this case study, the B and C bands of minimum and maximum degradation, respectively, are split into 5 and 15 levels. This implies increments of 10.8 μm for section B and 9 μm for section C. The total number of states is divided into 22 subdivisions, 20 for B and C bands plus A and D. This aims to refine the control bands, as well as improve the estimation of the proposed model probabilities. So, the degradation prediction model was developed based on readings of the vibration sensors. The cumulative hours of operation since the last preventive maintenance, which had reestablished good-as-new condition, were also recorded. The bearing condition before the preventive maintenance is shown in Figure 8.Table 5
S m a x limits for 200 RPM.
State Maximum displacement [μm] A S m a x < 82 B 82 ≤Smax < 136 C 136 ≤Smax < 271 D S m a x ≥ 271Figure 8
Guide bearing overview.The deterioration process often accelerates as the operating time increases until reaching the remaining useful life. The expectation of the distribution of states in degradation prediction models presents features similar to those of exponential distributions [71]. The maximum displacement measurements of the Corumbá IV HPP by sensors follow this exponential behavior, which was used to calculate the sojourn times for the transition probabilities in the Markov Model.


#### 4.2.2. Transition Probabilities

To simulate how the transitions states will evolve, an exponential curve is plotted, supported by the exponential behavior observed through the measurements. This exponential curve, shown in (6), relates the operating time to Smax and allows calculation of the sojourn times in each state according to the bands defined in the previous section.(6)Smaxt=e3.1395924+0.00006944t.After calculating the sojourn times and knowing that on average Corumbá IV HPP operates 80.48 hours per start-stop cycle, it is possible to determine how many start-stop cycles, here calledn, are associated with each state. This information is used to estimate the probability of remaining in the home condition. Since there is only one transition between each state, the probability of going to the next stage is 1/n and of remaining in the home state is 1-1/n.


#### 4.2.3. Emission Probabilities

In this case study, the estimated values forSmax are adopted as reference values for the observed variable. These values identify the change between states, and therefore must be linked to a probability distribution that binds them to the different states. Since the emissions are discrete and have a certain imprecision, due to the subjective nature of the degradation state, such boundaries are treated as probabilities. In this case, these sets are represented as a symmetric triangular probability function. Thus, for a given value of Smax treated as median, the lower and upper limits are determined according to this criterion. The link of ISO 7919-5 absolute limits with the corresponding linguistic terms of the degradation states may result in a certain level of vagueness, with the interval range having attenuation to such weakness. The decision about a fair value in μm for such range is based on sensitivity. When the value is too low, the emission probabilities are concentrated around few corresponding states and the ability to predict chances of being in some other near state is affected. When the value is too high, the emission probabilities are widely distributed around the most probable state causing insecurity to the decision maker. Therefore, different interval range values should be tested to define which value fairly reflects such level of imprecision.


#### 4.2.4. Model Validation through Inspection and Maintenance Data

The plant’s maintenance records showed that when the power unit operating time reached 24,000 hours, it underwent preventive maintenance. On this occasion, no anomaly was identified in the guide bearing and therefore this was judged subjectively by the technicians as if it still remained in operating condition. However, after operating 36,000 hours, the power unit underwent new preventive maintenance, when severe scratching of the coating metal was detected. As a result, the guide bearing was identified as being in a state of high deterioration and likely to compromise the regular operating condition of the power unit. The guide bearing was thus replaced and the equipment was classified as good as new again. This information about the guide bearing degradation identified in the 24,000 and 36,000 hour inspections was used to validate the model parameterization.


#### 4.2.5. HMM Parameters Inference

Due to the characteristic ofSmax progression, as shown in Figure 7, the left-right HMM [29, 72] is considered. This means that the degradation process is irreversible and once the equipment is in one specific state condition, it cannot move to the previous state anymore, remaining in the same state or advancing to the next one. The usage of HMM for diagnostics applies two complementary steps: estimation of the parameters and posterior probabilities. The former comprises assessment of both the transition and emission matrices. The Smax values obtained from the 14 vibration samples and the sojourn time of the states are used to estimate the transition probabilities. On the other hand, the emission probabilities feature the thresholds as stated by ISO 7919-5 and a symmetric triangular probability function. Both emission and transition probabilities are inputs of the inferred HMM model and are calculated by means of the forward-backward algorithm [73, 74]. In the second phase, the HMM parameters are used to estimate the posterior probabilities given a sequence of two hypothetical Smax values. Following the forward and backward probabilities [29], each value represents the guide bearing degradation state at a specific time. In this case study, the observed Smax values are set hypothetically as 253, 262, and 271 μm, coincident with the boundaries of the 19th, 20th, and 21st degradation states, in a scale where the 1st state is “good as new” and the 22nd state is “imminent failure.” The basic idea is to evaluate how responsive is the proposed model in estimating the degradation condition of the guide bearing, which is essential for the CBM.Table6 shows the posterior probabilities of being in states 19, 20, 21, or 22 at times of T1 and T2, relative to the vibration measurements onwards. These states are referred to as “excessive deterioration,” “prefailure level 1,” “prefailure level 2,” and “imminent failure,” respectively. The values for the operating times T1 and T2 depend on the state determined by Smax⁡, where the relationship is defined by (6). In the present simulation, the upper and lower limits of the triangular function sufficient to calculate the emission probabilities were set at 9 μm relative to Smax⁡. In a first round, when two consecutive readings were 253 and 262 μm, the probability of having the degradation referred by state 20 after T1 equal to 521.6 hours was 76% and after T2 equal to 1,025 hours was 93%. In a second round, since the next Smax⁡ reading increased to 271 μm, the probability of being in state 21 after T1 equal to 503.4 hours was 73% and after T2 equal to 989.8 hours was 93%. In a third round, when the two consecutive readings were equal to 271 μm, the probabilities of having the degradation referred by states 21 and 22 after T1 equal to 486.4 hours were 79% and 21%, respectively. In this round, the probabilities of being in the same states 21 and 22 after T2 equal to 972.8 hours changed to 71% and 29%, respectively.Table 6
Posteriori probabilities.
S m a x [μm] T 1 T 2 State 19 20 21 22 19 20 21 22 253–262 24% 76% 0 0 0 93% 7% 0 262–271 0 27% 73% 0 0 0 93% 7% ≥271 0 0 79% 21% 0 0 71% 29%


## 4.1. Data Processing

It is strongly recommended to acquire the start-up vibration measurement from stand still until the no-load speed condition, to assure complete acquisition of vibration data. The process takes approximately five minutes. According to the technical data sheet of the vibration monitoring system [55], the maximum sampling rate is 128 Hz regarding a five-minute memory buffer. Obeying the Nyquist criteria, the information is completely determined up to 64 Hz, as can be inferred in the next subsection.


### 4.1.1. Frequency Spectrum

Applying (2), for a 2 Hz Haning window, the result is pictured in Figure 5. The surface relates the normalized vibration power density V=V(f,t) to the excursion in time t and frequency f, given in minutes and hertz, respectively. The vibration function domain is DV={(f,t)∈R2∣0≤f≤64,0≤t≤5}.Figure 5
Short term Fourier transform of a transient state vibration.The normalized absolute value given by its image presents clear nonstationary behavior due to its nonconstant frequency arrangement. The hill over low frequencies presents decay up to 10 Hz, suggesting an intense flicker noise content due to the signal’s low sampling frequency [57]. Besides this, there are some random noise clouds over various frequencies of the AC component, even over the fundamental harmonic, which indicate distributed thermal noise pollution. There is also a constant DC level even after steady state is achieved, where the harmonic frequencies all act in the same direction, in approximately 250 s.


### 4.1.2. DC Component Filtering

The DC filtering techniques discussed in Section3.2.2 are compared using a FIR filter and the mean value. Table 2 shows the amplitude of the DC components filtered by each method and their relative difference. The comparison infers that their relative differences are not substantial, so the mean value technique is adopted as an advantage, by not distorting the phase frequency response.Table 2
DC filtering techniques result.
Acquisition Mean value FIR filter Difference 16/04/2017 1089.76 [μm] 1087.34 [μm] 0.22 [%] 02/05/2017 1023.05 [μm] 1022.54 [μm] 0.05 [%] 20/03/2017 1018.53 [μm] 1018.09 [μm] 0.04 [%] 23/05/2017 1035.73 [μm] 1036.67 [μm] 0.09 [%] 09/09/2016 1063.90 [μm] 1063.50 [μm] 0.04 [%] 17/04/2017 1064.40 [μm] 1063.52 [μm] 0.08 [%]


### 4.1.3. AC Component Filtering

Figure5 shows random noise over the entire surface domain of the unfiltered signal. The prominent flicker noise is more concentrated in regions up to 10 Hz and the thermal noise power density is well distributed in the whole frequency excursion. The wavelet packet used for denoising presents the Symlet 8 mother wavelet due to its good performance, with quasisinusoidal waveforms, and the denoising process considers soft thresholding at seven frequency levels. Table 3 shows a SNR comparison before and after the signal processing for some frequency harmonics [69, 70]. It can be noted that the signal processing considerably increased the relative signal power component over the noisy one. Figure 6 compares the original and filtered signals. Zooming reveals a large difference between each waveform, thus good adherence and filter behavior.Table 3
SNR comparison.
Harmonic Frequency [Hz] SNR before SNR after 0 0 37.39 317.62 1 3.33 35.68 304.38 3 10 25.99 205.45 6 20 17.77 191.3 8 26.66 15.16 132.9Figure 6
Comparison between original and filtered signals.


### 4.1.4. Sensitivity Analysis betweenSmax and Operation Time

The maximum displacement evolution over the operating time fits an exponential function nicely, as depicted in Figure7. In a smoothing process, when the mean extracted from a subset containing the k greatest numbers of a sample is its statistic, the exponential adherence among the statistics of the 14 samples is validated. Such confirmation has an important consequence on practical grounds, since the exponential relation is preserved even when two or more greater numbers are considered, mainly from the presence of outliers. This paper considers the relationship function for the maximum absolute value, k=1, over operation time, as observed in Table 4.Table 4
S m a x values over operation time.
Acquisition k = 1 k = 10 k = 20 09/09/2016 149.92 147.71 146.59 19/12/2016 166.45 163.65 158.91 06/01/2017 169.22 161.93 156.87 19/01/2017 171.97 170.36 166.24 20/02/2017 178.09 175.69 170.64 21/02/2017 176.28 174.07 170.95 20/03/2017 183.35 180.22 174.52 05/04/2017 186.39 183.62 179.97 16/04/2017 188.74 181.79 177.00 17/04/2017 185.51 182.72 177.15 24/04/2017 190.84 183.91 178.70 02/05/2017 191.62 183.55 176.31 23/05/2017 196.98 191.45 184.93 30/05/2017 197.68 189.71 180.82Figure 7
Exponential adherence betweenSmax and operation time.


## 4.1.1. Frequency Spectrum

Applying (2), for a 2 Hz Haning window, the result is pictured in Figure 5. The surface relates the normalized vibration power density V=V(f,t) to the excursion in time t and frequency f, given in minutes and hertz, respectively. The vibration function domain is DV={(f,t)∈R2∣0≤f≤64,0≤t≤5}.Figure 5
Short term Fourier transform of a transient state vibration.The normalized absolute value given by its image presents clear nonstationary behavior due to its nonconstant frequency arrangement. The hill over low frequencies presents decay up to 10 Hz, suggesting an intense flicker noise content due to the signal’s low sampling frequency [57]. Besides this, there are some random noise clouds over various frequencies of the AC component, even over the fundamental harmonic, which indicate distributed thermal noise pollution. There is also a constant DC level even after steady state is achieved, where the harmonic frequencies all act in the same direction, in approximately 250 s.


## 4.1.2. DC Component Filtering

The DC filtering techniques discussed in Section3.2.2 are compared using a FIR filter and the mean value. Table 2 shows the amplitude of the DC components filtered by each method and their relative difference. The comparison infers that their relative differences are not substantial, so the mean value technique is adopted as an advantage, by not distorting the phase frequency response.Table 2
DC filtering techniques result.
Acquisition Mean value FIR filter Difference 16/04/2017 1089.76 [μm] 1087.34 [μm] 0.22 [%] 02/05/2017 1023.05 [μm] 1022.54 [μm] 0.05 [%] 20/03/2017 1018.53 [μm] 1018.09 [μm] 0.04 [%] 23/05/2017 1035.73 [μm] 1036.67 [μm] 0.09 [%] 09/09/2016 1063.90 [μm] 1063.50 [μm] 0.04 [%] 17/04/2017 1064.40 [μm] 1063.52 [μm] 0.08 [%]


## 4.1.3. AC Component Filtering

Figure5 shows random noise over the entire surface domain of the unfiltered signal. The prominent flicker noise is more concentrated in regions up to 10 Hz and the thermal noise power density is well distributed in the whole frequency excursion. The wavelet packet used for denoising presents the Symlet 8 mother wavelet due to its good performance, with quasisinusoidal waveforms, and the denoising process considers soft thresholding at seven frequency levels. Table 3 shows a SNR comparison before and after the signal processing for some frequency harmonics [69, 70]. It can be noted that the signal processing considerably increased the relative signal power component over the noisy one. Figure 6 compares the original and filtered signals. Zooming reveals a large difference between each waveform, thus good adherence and filter behavior.Table 3
SNR comparison.
Harmonic Frequency [Hz] SNR before SNR after 0 0 37.39 317.62 1 3.33 35.68 304.38 3 10 25.99 205.45 6 20 17.77 191.3 8 26.66 15.16 132.9Figure 6
Comparison between original and filtered signals.


## 4.1.4. Sensitivity Analysis betweenSmax and Operation Time

The maximum displacement evolution over the operating time fits an exponential function nicely, as depicted in Figure7. In a smoothing process, when the mean extracted from a subset containing the k greatest numbers of a sample is its statistic, the exponential adherence among the statistics of the 14 samples is validated. Such confirmation has an important consequence on practical grounds, since the exponential relation is preserved even when two or more greater numbers are considered, mainly from the presence of outliers. This paper considers the relationship function for the maximum absolute value, k=1, over operation time, as observed in Table 4.Table 4
S m a x values over operation time.
Acquisition k = 1 k = 10 k = 20 09/09/2016 149.92 147.71 146.59 19/12/2016 166.45 163.65 158.91 06/01/2017 169.22 161.93 156.87 19/01/2017 171.97 170.36 166.24 20/02/2017 178.09 175.69 170.64 21/02/2017 176.28 174.07 170.95 20/03/2017 183.35 180.22 174.52 05/04/2017 186.39 183.62 179.97 16/04/2017 188.74 181.79 177.00 17/04/2017 185.51 182.72 177.15 24/04/2017 190.84 183.91 178.70 02/05/2017 191.62 183.55 176.31 23/05/2017 196.98 191.45 184.93 30/05/2017 197.68 189.71 180.82Figure 7
Exponential adherence betweenSmax and operation time.


## 4.2. Degradation Model

### 4.2.1. Degradation Model Inputs

The runner nominal speed is 200 RPM, and according to [34], the Smax theoretical limits for a 200 RPM machine are shown in Table 5. In this case study, the B and C bands of minimum and maximum degradation, respectively, are split into 5 and 15 levels. This implies increments of 10.8 μm for section B and 9 μm for section C. The total number of states is divided into 22 subdivisions, 20 for B and C bands plus A and D. This aims to refine the control bands, as well as improve the estimation of the proposed model probabilities. So, the degradation prediction model was developed based on readings of the vibration sensors. The cumulative hours of operation since the last preventive maintenance, which had reestablished good-as-new condition, were also recorded. The bearing condition before the preventive maintenance is shown in Figure 8.Table 5
S m a x limits for 200 RPM.
State Maximum displacement [μm] A S m a x < 82 B 82 ≤Smax < 136 C 136 ≤Smax < 271 D S m a x ≥ 271Figure 8
Guide bearing overview.The deterioration process often accelerates as the operating time increases until reaching the remaining useful life. The expectation of the distribution of states in degradation prediction models presents features similar to those of exponential distributions [71]. The maximum displacement measurements of the Corumbá IV HPP by sensors follow this exponential behavior, which was used to calculate the sojourn times for the transition probabilities in the Markov Model.


### 4.2.2. Transition Probabilities

To simulate how the transitions states will evolve, an exponential curve is plotted, supported by the exponential behavior observed through the measurements. This exponential curve, shown in (6), relates the operating time to Smax and allows calculation of the sojourn times in each state according to the bands defined in the previous section.(6)Smaxt=e3.1395924+0.00006944t.After calculating the sojourn times and knowing that on average Corumbá IV HPP operates 80.48 hours per start-stop cycle, it is possible to determine how many start-stop cycles, here calledn, are associated with each state. This information is used to estimate the probability of remaining in the home condition. Since there is only one transition between each state, the probability of going to the next stage is 1/n and of remaining in the home state is 1-1/n.


### 4.2.3. Emission Probabilities

In this case study, the estimated values forSmax are adopted as reference values for the observed variable. These values identify the change between states, and therefore must be linked to a probability distribution that binds them to the different states. Since the emissions are discrete and have a certain imprecision, due to the subjective nature of the degradation state, such boundaries are treated as probabilities. In this case, these sets are represented as a symmetric triangular probability function. Thus, for a given value of Smax treated as median, the lower and upper limits are determined according to this criterion. The link of ISO 7919-5 absolute limits with the corresponding linguistic terms of the degradation states may result in a certain level of vagueness, with the interval range having attenuation to such weakness. The decision about a fair value in μm for such range is based on sensitivity. When the value is too low, the emission probabilities are concentrated around few corresponding states and the ability to predict chances of being in some other near state is affected. When the value is too high, the emission probabilities are widely distributed around the most probable state causing insecurity to the decision maker. Therefore, different interval range values should be tested to define which value fairly reflects such level of imprecision.


### 4.2.4. Model Validation through Inspection and Maintenance Data

The plant’s maintenance records showed that when the power unit operating time reached 24,000 hours, it underwent preventive maintenance. On this occasion, no anomaly was identified in the guide bearing and therefore this was judged subjectively by the technicians as if it still remained in operating condition. However, after operating 36,000 hours, the power unit underwent new preventive maintenance, when severe scratching of the coating metal was detected. As a result, the guide bearing was identified as being in a state of high deterioration and likely to compromise the regular operating condition of the power unit. The guide bearing was thus replaced and the equipment was classified as good as new again. This information about the guide bearing degradation identified in the 24,000 and 36,000 hour inspections was used to validate the model parameterization.


### 4.2.5. HMM Parameters Inference

Due to the characteristic ofSmax progression, as shown in Figure 7, the left-right HMM [29, 72] is considered. This means that the degradation process is irreversible and once the equipment is in one specific state condition, it cannot move to the previous state anymore, remaining in the same state or advancing to the next one. The usage of HMM for diagnostics applies two complementary steps: estimation of the parameters and posterior probabilities. The former comprises assessment of both the transition and emission matrices. The Smax values obtained from the 14 vibration samples and the sojourn time of the states are used to estimate the transition probabilities. On the other hand, the emission probabilities feature the thresholds as stated by ISO 7919-5 and a symmetric triangular probability function. Both emission and transition probabilities are inputs of the inferred HMM model and are calculated by means of the forward-backward algorithm [73, 74]. In the second phase, the HMM parameters are used to estimate the posterior probabilities given a sequence of two hypothetical Smax values. Following the forward and backward probabilities [29], each value represents the guide bearing degradation state at a specific time. In this case study, the observed Smax values are set hypothetically as 253, 262, and 271 μm, coincident with the boundaries of the 19th, 20th, and 21st degradation states, in a scale where the 1st state is “good as new” and the 22nd state is “imminent failure.” The basic idea is to evaluate how responsive is the proposed model in estimating the degradation condition of the guide bearing, which is essential for the CBM.Table6 shows the posterior probabilities of being in states 19, 20, 21, or 22 at times of T1 and T2, relative to the vibration measurements onwards. These states are referred to as “excessive deterioration,” “prefailure level 1,” “prefailure level 2,” and “imminent failure,” respectively. The values for the operating times T1 and T2 depend on the state determined by Smax⁡, where the relationship is defined by (6). In the present simulation, the upper and lower limits of the triangular function sufficient to calculate the emission probabilities were set at 9 μm relative to Smax⁡. In a first round, when two consecutive readings were 253 and 262 μm, the probability of having the degradation referred by state 20 after T1 equal to 521.6 hours was 76% and after T2 equal to 1,025 hours was 93%. In a second round, since the next Smax⁡ reading increased to 271 μm, the probability of being in state 21 after T1 equal to 503.4 hours was 73% and after T2 equal to 989.8 hours was 93%. In a third round, when the two consecutive readings were equal to 271 μm, the probabilities of having the degradation referred by states 21 and 22 after T1 equal to 486.4 hours were 79% and 21%, respectively. In this round, the probabilities of being in the same states 21 and 22 after T2 equal to 972.8 hours changed to 71% and 29%, respectively.Table 6
Posteriori probabilities.
S m a x [μm] T 1 T 2 State 19 20 21 22 19 20 21 22 253–262 24% 76% 0 0 0 93% 7% 0 262–271 0 27% 73% 0 0 0 93% 7% ≥271 0 0 79% 21% 0 0 71% 29%


## 4.2.1. Degradation Model Inputs

The runner nominal speed is 200 RPM, and according to [34], the Smax theoretical limits for a 200 RPM machine are shown in Table 5. In this case study, the B and C bands of minimum and maximum degradation, respectively, are split into 5 and 15 levels. This implies increments of 10.8 μm for section B and 9 μm for section C. The total number of states is divided into 22 subdivisions, 20 for B and C bands plus A and D. This aims to refine the control bands, as well as improve the estimation of the proposed model probabilities. So, the degradation prediction model was developed based on readings of the vibration sensors. The cumulative hours of operation since the last preventive maintenance, which had reestablished good-as-new condition, were also recorded. The bearing condition before the preventive maintenance is shown in Figure 8.Table 5
S m a x limits for 200 RPM.
State Maximum displacement [μm] A S m a x < 82 B 82 ≤Smax < 136 C 136 ≤Smax < 271 D S m a x ≥ 271Figure 8
Guide bearing overview.The deterioration process often accelerates as the operating time increases until reaching the remaining useful life. The expectation of the distribution of states in degradation prediction models presents features similar to those of exponential distributions [71]. The maximum displacement measurements of the Corumbá IV HPP by sensors follow this exponential behavior, which was used to calculate the sojourn times for the transition probabilities in the Markov Model.


## 4.2.2. Transition Probabilities

To simulate how the transitions states will evolve, an exponential curve is plotted, supported by the exponential behavior observed through the measurements. This exponential curve, shown in (6), relates the operating time to Smax and allows calculation of the sojourn times in each state according to the bands defined in the previous section.(6)Smaxt=e3.1395924+0.00006944t.After calculating the sojourn times and knowing that on average Corumbá IV HPP operates 80.48 hours per start-stop cycle, it is possible to determine how many start-stop cycles, here calledn, are associated with each state. This information is used to estimate the probability of remaining in the home condition. Since there is only one transition between each state, the probability of going to the next stage is 1/n and of remaining in the home state is 1-1/n.


## 4.2.3. Emission Probabilities

In this case study, the estimated values forSmax are adopted as reference values for the observed variable. These values identify the change between states, and therefore must be linked to a probability distribution that binds them to the different states. Since the emissions are discrete and have a certain imprecision, due to the subjective nature of the degradation state, such boundaries are treated as probabilities. In this case, these sets are represented as a symmetric triangular probability function. Thus, for a given value of Smax treated as median, the lower and upper limits are determined according to this criterion. The link of ISO 7919-5 absolute limits with the corresponding linguistic terms of the degradation states may result in a certain level of vagueness, with the interval range having attenuation to such weakness. The decision about a fair value in μm for such range is based on sensitivity. When the value is too low, the emission probabilities are concentrated around few corresponding states and the ability to predict chances of being in some other near state is affected. When the value is too high, the emission probabilities are widely distributed around the most probable state causing insecurity to the decision maker. Therefore, different interval range values should be tested to define which value fairly reflects such level of imprecision.


## 4.2.4. Model Validation through Inspection and Maintenance Data

The plant’s maintenance records showed that when the power unit operating time reached 24,000 hours, it underwent preventive maintenance. On this occasion, no anomaly was identified in the guide bearing and therefore this was judged subjectively by the technicians as if it still remained in operating condition. However, after operating 36,000 hours, the power unit underwent new preventive maintenance, when severe scratching of the coating metal was detected. As a result, the guide bearing was identified as being in a state of high deterioration and likely to compromise the regular operating condition of the power unit. The guide bearing was thus replaced and the equipment was classified as good as new again. This information about the guide bearing degradation identified in the 24,000 and 36,000 hour inspections was used to validate the model parameterization.


## 4.2.5. HMM Parameters Inference

Due to the characteristic ofSmax progression, as shown in Figure 7, the left-right HMM [29, 72] is considered. This means that the degradation process is irreversible and once the equipment is in one specific state condition, it cannot move to the previous state anymore, remaining in the same state or advancing to the next one. The usage of HMM for diagnostics applies two complementary steps: estimation of the parameters and posterior probabilities. The former comprises assessment of both the transition and emission matrices. The Smax values obtained from the 14 vibration samples and the sojourn time of the states are used to estimate the transition probabilities. On the other hand, the emission probabilities feature the thresholds as stated by ISO 7919-5 and a symmetric triangular probability function. Both emission and transition probabilities are inputs of the inferred HMM model and are calculated by means of the forward-backward algorithm [73, 74]. In the second phase, the HMM parameters are used to estimate the posterior probabilities given a sequence of two hypothetical Smax values. Following the forward and backward probabilities [29], each value represents the guide bearing degradation state at a specific time. In this case study, the observed Smax values are set hypothetically as 253, 262, and 271 μm, coincident with the boundaries of the 19th, 20th, and 21st degradation states, in a scale where the 1st state is “good as new” and the 22nd state is “imminent failure.” The basic idea is to evaluate how responsive is the proposed model in estimating the degradation condition of the guide bearing, which is essential for the CBM.Table6 shows the posterior probabilities of being in states 19, 20, 21, or 22 at times of T1 and T2, relative to the vibration measurements onwards. These states are referred to as “excessive deterioration,” “prefailure level 1,” “prefailure level 2,” and “imminent failure,” respectively. The values for the operating times T1 and T2 depend on the state determined by Smax⁡, where the relationship is defined by (6). In the present simulation, the upper and lower limits of the triangular function sufficient to calculate the emission probabilities were set at 9 μm relative to Smax⁡. In a first round, when two consecutive readings were 253 and 262 μm, the probability of having the degradation referred by state 20 after T1 equal to 521.6 hours was 76% and after T2 equal to 1,025 hours was 93%. In a second round, since the next Smax⁡ reading increased to 271 μm, the probability of being in state 21 after T1 equal to 503.4 hours was 73% and after T2 equal to 989.8 hours was 93%. In a third round, when the two consecutive readings were equal to 271 μm, the probabilities of having the degradation referred by states 21 and 22 after T1 equal to 486.4 hours were 79% and 21%, respectively. In this round, the probabilities of being in the same states 21 and 22 after T2 equal to 972.8 hours changed to 71% and 29%, respectively.Table 6
Posteriori probabilities.
S m a x [μm] T 1 T 2 State 19 20 21 22 19 20 21 22 253–262 24% 76% 0 0 0 93% 7% 0 262–271 0 27% 73% 0 0 0 93% 7% ≥271 0 0 79% 21% 0 0 71% 29%


## 5. Conclusion

The signal processing methods applied in this article, such as mean DC level and wavelet packet, presented fair performance in handling nonstationary data. Flicker and thermal noise were attenuated and the obtained signals presented high SNR. The resulting vibration values and the maximum axis displacement were consistent for use as inputs for the proposed stochastic model. The proposed method has some important aspects. First, the vibration readings are acquired in a real operational environment, which means application of standard vibration monitoring systems. While on the one hand this method has the downside of reading imprecision, on the other it enables the operator himself to collect the signals without external assistance. The adoption of HMM has two aspects, both important features of the proposed method. First, it is related to the main aspect of a Markov technique that allows the operator to carry out the guide bearing diagnosis without the need of putting the runner on hold to perform inspections. The second is concerned with the possibility of modeling guide bearing diagnosis by means of corrective maintenance records, therefore based on past experience. In the case of Corumbá IV HPP, the filtering and HMM parameters are adapted to a routine diagnosis system. In this case, the plant operator simply acquires the axial vibration signals during the runner start-up and feeds the model with such raw data. The automated routine calculates the probabilities of the guide bearing being at certain stages of degradation. Therefore, the proposed method has substantial advantages by reducing failure risk and maintenance costs.



---

*Source: 5981089-2018-01-02.xml*