Blind Deconvolution Denoising for Helicopter Vibration Signals
558
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
Blind Deconvolution Denoising for Helicopter Vibration Signals Bin Zhang, Senior Member, IEEE, Taimoor Khawaja, Romano Patrick, and George Vachtsevanos, Senior Member, IEEE
Abstract—Critical aircraft assets are required to be available when needed, while exhibiting attributes of reliability, robustness, and high confidence under a variety of flight regimes, and maintained on the basis of their current condition rather than on the basis of scheduled maintenance practices. New and innovative technologies must be developed and implemented to address these concerns. Condition-based maintenance requires that the health of critical components/systems be monitored and diagnostic/prognostic strategies be developed to detect and identify incipient failures and predict the failing component’s remaining useful life. Typically, vibration and other key indicators onboard an aircraft are severely corrupted by noise, thus curtailing the ability to accurately diagnose and predict failures. This paper introduces a novel blind deconvolution denoising scheme that employs a vibration model in the frequency domain and attempts to arrive at the true vibration signal through an iterative optimization process. Performance indexes are defined and data from a helicopter are used to demonstrate the effectiveness of the proposed approach.
Fig. 1.
Crack of planetary gear carrier plate of the UH-60A helicopter.
Index Terms—Blind deconvolution, planetary gear train, vibration signal denoising.
I. INTRODUCTION AULT DETECTION, isolation, and identification, especially for safety critical components and subsystems, have recently drawn increasing interest in the condition-based maintenance (CBM) community [1]–[4]. Epicyclic, or planetary, gear trains are widely used in the main transmission of many systems, such as the helicopter and other aircraft [5], [7]. This kind of gear system consists of an inner sun gear, two or more rotating planet gears, a stationary outer ring gear, and a planetary carrier. The planet gears, which surround the sun gear, are riding on the planetary carrier and also rotating inside the outer ring gear. In the operation, torque is transmitted through the sun gear to the planets and planetary carrier. The carrier, in turn, transmits torque to the main rotor shaft and blades [6]–[9]. The gearbox is a critical component that directly impacts the safety and performance of the aircraft.
F
Manuscript received October 8, 2006; revised April 25, 2008. Current version published October 8, 2008. Recommended by Technical Editor M.-Y. Chow. This work was supported in part by the Defense Advanced Research Projects Agency (DARPA) and in part by Northrop Grumman under DARPA Prognosis Program Contract HR0011-04-C-003. B. Zhang and T. Khawaja are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail: [email protected]). R. Patrick is with Impact Technologies, LLC, Rochester, NY 14623 USA (e-mail: [email protected]). G. Vachtsevanosis is with Georgia Institute of Technology, Atlanta, GA 30332 USA, and also with Impact Technologies, LLC, Rochester, NY 14623 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMECH.2008.2002324
The UH-60A Blackhawk and Seahawk main transmission employs a five-planet-gear epicyclic gear system [7]. Recently, a crack in the planetary carrier was found, as shown in Fig. 1. This resulted in flight restrictions on many helicopters. Manual inspection of all transmissions is not only costly, but also time prohibitive. A CBM based on vibration signal analysis is a cost-effective solution [16]. To measure the vibration signals, the transducer is mounted at a fixed point on the frame of the gearbox. If a fault occurs on the gear system, it is expected that the vibration signal will exhibit a characteristic signature that reveals the presence, severity, and location of the fault. In fault diagnosis and failure prognosis, gear health is typically described in terms of feature or condition indicator (CI) values, which are extracted from the vibration signals [7]. The success of diagnosis and prediction of the remaining useful life of a gear system highly depends on the quality of these features. However, during the course of a flight, the vibration signals derived from accelerometers often have multiple excitation sources, and in most cases, information about the noise characteristics is not available. The measured output signal is, therefore, the combined effect from these different sources. Noise signals will degrade the quality of the features, and consequently, the performance of diagnosis and prognosis, such as the detection threshold, the accuracy of remaining useful life of components, etc. For example, when a fault is at its early stage, the indication in the vibration signals and the features is often masked by noise. It is, therefore, important to develop an efficient and reliable denoising algorithm to remove noise signals as much as possible and make the fault characteristics perceptible in the extracted features.
1083-4435/$25.00 © 2008 IEEE Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
Fig. 2.
559
Blind deconvolution denoising scheme.
The widely used denoising technique in planetary gears is time-synchronous averaging (TSA) [7]–[9]. To realize TSA, a pulse signal that is synchronized to the rotation of a gear indicating the start of an individual revolution is assumed to be available. The sampling data collected from different revolutions are interpolated to the same length, and then, ensemble averaging is performed. The interpolation operation is often computationally intensive. A TSA technique in the frequency domain has been proposed to reduce the computation drastically [17]. Some features are successfully extracted from the TSA data and used in the fault diagnosis and prognosis [8], [16], [17]. The TSA technique enhances the components at the frequencies that are multiples of the shaft frequency, which are often related to the meshing of the gear teeth [8], [16]. At the same time, it tends to average out the external random disturbances and noise that are asynchronous with the rotation of the gear. However, in the frequency domain, it is possible that some noise components appear at frequencies that are multiples of the shaft frequency and enhanced by TSA processing. Hence, it is necessary to further denoise the TSA data to achieve better SNR and improve the quality of the features extracted from the TSA data. The vibration signal collected from the transducer on the frame of the gearbox is amplitude modulated and corrupted by additive noise. Blind deconvolution has shown very promising results for similarly formulated problems in image and speech processing [10], [11], but it is rarely applied to vibration signals [12]–[15]. In this paper, a blind deconvolution scheme, which consists of vibration analysis, nonlinear projection, and cost function optimization, is employed to the denoising of the vibration signals from a UH-60A helicopter gearbox, whose planetary carrier has a seeded crack fault growing with the evolving operation of the system. The analysis of the system, the theoretical basis of the blind deconvolution algorithm, and the proposed performance indexes are addressed in detail. Experimental results show the effectiveness of the proposed denoising scheme. II. DENOISING SCHEME ARCHITECTURE The architecture of the proposed blind deconvolution denoising scheme is shown in Fig. 2. In this scheme, nonlinear projection, which is based on vibration analysis in the frequency domain, and cost function minimization are critical components, which will be described in later sections. The proposed denoising approach first transforms the measured vibration signal s(t) ¯ ), which is an estimate of to S(f ). Then, an inverse filter Z(f
Fig. 3.
Configuration of an epicyclic gear system.
the inverse of the modulating signal in the frequency domain, is convoluted with the measured vibration signal S(f ) to demodulate S(f ) and give a rough estimate of the noise-free vibration ¯ ). The signal B(f ¯ ) passes through the nonlinear signal B(f ¯ ) to a subspace that contains only projection, which maps B(f known characteristics of the vibration signal, to yield Bn l (f ). ¯ ) and Bn l (f ) is denoted as E(f ). The difference between B(f ¯ ) iteratively to minimize the E(f ) and when By adjusting Z(f ¯ ) → B(f ) can be E(f ) reaches a minimal value, the signal B(f regarded as the denoised vibration signal. Through an inverse Fourier transform, the denoised vibration signal in the time domain can be obtained as well.
III. SYSTEM DESCRIPTION The vibration signals are derived from the main transmission gearbox of the UH-60A helicopter. The gearbox is an epicyclic gear system with five planet gears, whose configuration is shown in Fig. 3. The following analysis assumes an ideal system. The accelerometer is mounted at a fixed point on the annulus gear at the position θ = 0. Since the vibration signal is generated from the meshing of gear teeth and the planetary gears are rotating inside the angular gear, the vibration signal is amplitude modulated to the static accelerometer. The observed vibration amplitude will be large when the planetary gear is close to the accelerometer, and it will be small when the planetary gear is far. Suppose that there is only one planetary gear, then the vibration observed by the transducer should have the largest amplitude when the planetary gear is at θ = 0, 2π, 4π, . . .. Similarly, the vibration should have the smallest amplitude at θ = π, 3π, . . ., i.e., the teeth meshing vibration signal is amplitude modulated as the planetary carrier rotates. In the ideal case, the Np = 5 planetary gears are evenly spaced and suppose that the planetary carrier has a rotation frequency of f s . Then, the planetary gear p at time instant t has a phase of p−1 θp = 2π f s t + . Np
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
(1)
560
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
In this case, the amplitude-modulating signal for planetary gear p can be written in the time domain as ap (t) =
N
αn cos(nθp )
(2)
n =−N
where N is the number of sidebands about the harmonics under consideration and αn is the amplitude of the component of the modulating signal at frequency nf s . The “sidebands” are defined as the frequency components that appear as a harmonically spaced series [6]. In the ideal case, it is natural to assume that all mesh vibrations generated from different planetary gears are of the same amplitude but of different phases. The annulus gear has Nt = 228 teeth. Since the speed at which teeth meshing takes place is proportional to the angular velocity of the planetary carrier, the meshing vibration appears at frequencies Nt f s [6]. In addition, suppose that, in the frequency domain, the meshing vibration signal has amplitudes βm at its harmonics, then the vibration signal generated from planet gear p can be written in the form of bp (t) =
M
βm sin(mNt θp )
Fig. 4. Superposition of vibration from different gears with N p = 5. (a) mN t + n = kN p . (b) mN t + n = kN p .
(3)
m =1
where M is the number of harmonics under consideration. Then, to the static accelerometer, the vibration signal of planetary gear p appears as the product of the meshing vibration signal and the amplitude-modulating signal. It is denoted as yp (t) and is given as N M 1 αn βm sin((mNt + n)θp ). 2 m =1 n =−N (4) When there are more than one, say Np , planetary gears, the vibration signal observed by the accelerometer is the superposition of the Np vibration signals generated from Np different planetary gears. This superposition vibration signal has the form
Fig. 5.
Spectrum of combined signal with N t = 228 and N p = 5.
yp (t) = ap (t)bp (t) =
Np M N 1 αn βm sin((mNt + n)θp ) y(t) = 2 p=1 m =1 n =−N
N 1 mNt + n = αn βm sin 2π(p − 1) (5) 2 p=1 m =1 Np Np
M
n =−N
where (1) and the fact that sin(2kπ + θ) = sin(θ) for any integer k are used. Since the planetary gears are evenly spaced, the phase angel of the sidebands will be evenly spaced along 2π [6]. From (5), it is obvious that if sideband mNt + n is not a multiple of Np and (mNt + n)/Np has a remainder of γ, the vibration components from different gears are evenly spaced by the angle 2γπ/Np . In this case, when the vibrations generated from different planetary gears are combined, these sidebands add destructively and become zero, as illustrated in Fig. 4(a), in which φp,m ,n with 1 ≤ p ≤ 5 indicates the frequency components of gear p. Therefore, the frequency components appear at sidebands where mNt + n = kNp are termed as nonregular
meshing components (NonRMC). On the contrary, if sideband mNt + n is a multiple of Np , the remainder of (mNt + n)/Np will be zero. In this case, the vibration components from different gears do not have a phase difference. When the vibration signals from different planetary gears are combined, these sidebands add constructively and are reinforced, as illustrated in Fig. 4(b). These frequency components that appear at sidebands where mNt + n = kNp are termed as regular meshing components (RMC) or apparent sidebands. This process of frequency components adding destrutively/constructively finally generates asymmetrical sidebands. A partial frequency spectrum of the sidebands about the firstorder harmonic that illustrates this asymmetry is shown in Fig. 5, where Np = 5 and Nt = 228. Note that the peak of the spectrum does not appear at n = 0 (order 228) but at n = −3 (order 225), n = 2 (order 230), etc. The largest spectral amplitude (also known as dominant sideband) appears at the frequency closest to the gear meshing frequency. The vibration analysis in the frequency domain and previous research results [6], [7] suggest that, for an ideal system, only terms at frequencies that are multiples of the number of planetary gears survive. Then, the Fourier transform of the vibration data can be written as Y (f ) = fn l (γm ,n (mNt + n)f s )
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
(6)
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
where γm ,n is the magnitude of the spectral amplitude at (mNt + n)f s and fn l is the nonlinear projection 1, if mNt + n is a multiple of Np fn l = (7) 0, otherwise. From the previous analysis of the system vibration behavior, the following assumptions are made and used in the blinddeconvolution scheme. 1) The sidebands of the vibration spectra exist only at frequencies (with a resolution of the planetary carrier rotation frequency f s ) that are multiples of the number of planetary gears. 2) The phase of each sideband depends on the geometrical angular position of the planet and is given by θp (mNt + n) for sideband (m, n). 3) The amplitude of the modulating signal a(t) decreases monotonically on either side of the maximum until it reaches the minimum. Note that this assumption is somewhat restrictive. IV. BLIND DECONVOLUTION DENOISING SCHEME
(8)
where s(t) is the measured vibration signal, b(t) is the noise-free unmodulated vibration signal, a(t) is the modulating signal, and n(t) is the cumulative additive noise. Note that the modulating signal a(t) is itself affected by noise in the system. Let a ˆ(t) denote the ideal or noise-free modulating signal and na (t) the noise introduced in this signal. Consequently, a(t) can be represented as a(t) = a ˆ(t) + na (t).
(9)
Thus, (8) can be rewritten as follows: s(t) = (ˆ a(t) + na (t))b(t) + n(t) =a ˆ(t)b(t) + n ˆ (t)
(10)
where n ˆ (t) = na (t)b(t) + n(t) contains the total additive noise in the system. On the other hand, the factor a ˆ(t) describes the multiplicative noise in the system. The goal for a denoising scheme, such as the one described here, is to recover the unknown vibration signal b(t) from the observed signal s(t) given partial information about the noise sources and characteristics of the vibration signal. A typical approach would be to find the inverse of a ˆ(t) zˆ(t) =
1 a ˆ(t)
(11)
such that
¯b(t) = s(t)¯ z (t).
(13)
If partial knowledge about how the plate system is influenced by the modulating signal a ˆ(t) and a reasonable understanding of the true vibration signal are available, then the ideal characteristics of the vibration signal can be obtained by projecting this estimated signal ¯b(t) into a subspace with only the known ideal characteristics of the vibration signal to yield a refined signal bn l (t). Since this nonlinear projection, as the subscript signifies, removes all uncharacterisitic components that exist in the rough estimate ¯b(t), it is necessary to stress the importance of a good understanding of the underlying process. An iterative scheme then refines these results by minimizing the error between the two signals ¯b(t) and bn l (t), i.e., (14)
Previous research results detail the spectral characteristics of vibration signals for rotating equipment [6], [7]. It is, therefore, appropriate to investigate the measured noisy vibration in the frequency domain. The convolution theorem states that the product of two signals in the time domain is equivalent to their convolution in the frequency domain. Thus, model (10) can be written in the frequency domain as ˆ ) ∗ B(f ) + N ˆ (f ) S(f ) = A(f
(15)
ˆ ), B(f ), with ∗ being the convolution operator and S(f ), A(f ˆ and N (f ) are the Fourier transforms of s(t), a ˆ(t), b(t), and n ˆ (t), respectively. Then, the goal in the frequency domain is to recover B(f ). Writing (13) in the frequency domain, we have ¯ ) = S(f ) ∗ Z(f ¯ ) B(f
(16)
¯ ) and Z(f ¯ ) being the Fourier transforms of ¯b(t) and with B(f ¯ ) through the nonlinear projecz¯(t), respectively. Passing B(f tion, it yields Bn l (f ). Then, in the frequency domain, we will ¯ ) minimize the difference between Bn l (f ) and B(f ¯ ) − Bn l (f ). min E(f ) = min B(f
(17)
¯ ) to minimize (17). When it The iterative process refines Z(f ¯ reaches the minimal value, Z(f ) converges to Z(f ). Then, with ¯ ) in (16), a good estimate for B(f ) is this Z(f ) replacing Z(f obtained as B(f ) = S(f ) ∗ Z(f ).
(18)
Last, the estimate is transformed back into the time domain to recover the noise-free vibration signal b(t) ∞ ei2π f t B(f ) df . (19) b(t) = F −1 B(f ) = −∞
b(t) = (s(t) − n ˆ (t))ˆ z (t) = s(t)ˆ z (t) − n ˆ (t)ˆ z (t).
Note, however, that little can be assumed about n ˆ (t), a ˆ(t) is not available, and zˆ(t) is not applicable. To solve this problem, rather than using zˆ(t), we propose an iterative denoising scheme that starts with z¯(t), a very rough initial estimate of the inverse of the modulating signal, which demodulates the observed signal s(t) to give a rough noise-free vibration signal ¯b(t)
min e(t) = min ¯b(t) − bn l (t).
From the vibration analysis of the gearbox, we know that the vibration signals collected from the transducer are amplitude modulated [6]. Multiple sources of noise may further corrupt the signal. A simplified model for such a complex signal may be defined as s(t) = a(t)b(t) + n(t)
561
(12)
To solve this problem, the following additional assumptions are made (continued from the previous section).
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
562
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
TABLE I GROUND TRUTH DATA OF CRACK LENGTH (IN INCHES)
1) Z(f ) exists and Z(f ) < ∞. 2) Since the modulating signal a ˆ(t) is always positive, its inverse should be positive too. Hence, the Fourier transform of its inverse Z(f ) contains a dc component. With the previous two assumptions being taken into consideration, the cost function is defined as J=
[B(f ) − Bn l (f )]2 +
2 Z(f ) − 1
Fig. 6.
Growth of crack length.
Fig. 7.
(a) SNR at 40% torque. (b) SNR at 100% torque.
(20)
f ∈D s u p
where Dsup is the frequency range that contains the main vibration information. Because of the periodic fade of signal spectra between harmonics, a window centered at harmonic frequencies is used to define critical frequencies. All these windows form the support Dsup . In (20), assumptions 1) and 2) are used to arrive at the second term to avoid an all-zero inverse filter Z(f ), which leads to the trivial solution for error minimization. Moreover, an iterative optimization routine is required to implement this scheme. The iterative conjugate gradient method is called upon to address the optimization problem. This method has faster convergence rate in general as compared to the steepest descent method [11]. V. EXPERIMENTAL RESULTS The planetary carrier has a seeded crack with an initial length of 1.344 in. The crack, as shown in Fig. 1, grows with the evolving operation of the gearbox. The gearbox is installed on a test bed and operates over a large number of ground-air-ground (GAG) cycles at different torque levels. This way, vibration data are acquired at different torque levels and different crack lengths. In each GAG cycle, the torque increases from 20% to 40%, and finally to 100%, and then, decreases to 20% for the next cycle. A. Actual Crack Growth The ground truth crack length data at discrete GAG cycles are available and tabulated in Table I. With these data, the crack length for GAG cycles from 1 to 1000 can be obtained via interpolation, which results in the crack length growth curve shown in Fig. 6. B. Signal-to-Noise Ratio The SNR, before and after denoising, is investigated, and the results at 40% and 100% torque levels are shown in Fig. 7(a) and (b), respectively. The improvement in the SNR value is significant.
C. Feature Performance Indexes Although the blind deconvolution routine shows a significant improvement in the SNR, it is desirable that it improves also the quality of the features or condition indicators. The accuracy and precision of mappings between the evolution of features and the actual crack growth have an important impact on the performance of diagnostic and prognostic algorithms and the performance of the CBM system overall. To evaluate the quality of the features, two performance indexes are introduced. The first performance index is an accuracy measure defined as the linear correlation between the feature values and the crack length growth along the GAG cycle axis. Suppose that x is the
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
feature vector and y the crack growth curve with x ¯ and y¯ being their means, respectively. The correlation coefficient between them is ss2xy (21) CCR(x, y) = ssxx ssy y where ss ¯)(yi − y¯), ssxx = (xi − x ¯)2 , and i −x xy = (x 2 ssy y = (yi − y¯) , respectively. The extracted feature will be used to map and predict the crack growth. Hence, a high correlation coefficient is expected to generate an accurate estimate of actual crack length and is preferred. Then, for a good feature, the value of the correlation coefficient should be near to 1. Because of changes in operating conditions and disturbances, the feature values along a GAG cycle are often very noisy and need to be smoothed through a low-pass filtering operation. This filtering will not change the characteristics of the feature curve. In this case, the linear correlation is calculated based on the smoothed feature curve. The second performance index is a precision measure corresponding to a normalized measure of the signal dispersion. It is termed as percent mean deviation (PMD) and defined as n x i −˜x i PMD(x, x ˜) =
i=1
lx
x ˜i
× 100
Fig. 8.
563
Groups and D su p at the first harmonic with N t = 228 and N p = 5.
(22)
˜ where lx is the number of entities in the feature vector x and x is the smoothed feature vector. In the fault detection and prognosis algorithms, the extracted feature values are used as a measurement input. Therefore, this performance index is closely related to the detection threshold and precision of the prognosis or prediction of the remaining useful life. From a precise feature, the fault detection and prognosis algorithms can detect the incipient failure and/or predict the fault growth with a high confidence. PMD for a precise feature should have a very small value close to 0. D. Experimental Results As the previous analysis shows, only RMC, which include the dominant and apparent components appearing at frequencies that are multiples of Np , will survive in the vibration signal. However, the actual system is not ideal, particularly under faulty conditions. Then, the NonRMC, which appear at frequencies that are not multiples of Np , will not be zero. Hence, the frequency spectra are divided into two categories: RMC and NonRMC. To simplify the description and relevant computations, Dsup is defined in terms of groups of frequency spectra. From the previous spectra analysis, for a gearbox with Np planetary gears, the RMCs are evenly spaced by Np − 1 sidebands. Then, it is natural to define the group size as Np . For the five-planetarygear gearbox in this research, the group size is defined as 5. Then, Dsup is given as the number of groups on each side of the dominant component. The concepts of “group” and “Dsup ” are illustrated in Fig. 8, where only the groups and part of Dsup at the first harmonic are shown. The entire Dsup is the combination of these groups at all the harmonics under consideration.
Fig. 9. Feature SBR(2) versus GAG before and after denoising. (a) SBR(2) at 40% torque. (b) SBR(2) at 100% torque.
The relative size of the NonRMC sidebands to the RMC sidebands in Dsup or part of Dsup may indicate the presence of a fault on the planetary carrier. As the crack grows, the operating condition deviates from the ideal one. In this case, the constructive superposition of vibration signals from different planetary gears is attenuated while the destructive superposition is reinforced. This phenomenon results in smaller RMC and larger NonRMC sidebands. This indicates that the ratio between RMC and NonRMC sidebands can be used as a characteristic feature. Based on this notion, features or condition indicators are
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
564
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
TABLE II CORRELATION COEFFICIENTS OF THE FEATURES
TABLE III PERCENT MEAN DEVIATION OF THE FEATURES
respectively. In the tables, D-N simply means de-noised. It is clear that for most cases, substantial feature improvements are achieved via the application of the denoising routine. It is worth noting that the influence of crack on vibration spectra at low torque level 20% is not as remarkable as those at high torque levels. Since the nonlinear projection in the developed denoising routine depends on the analysis of vibration spectra, the application of this routine shows more reliable performance on high torque levels. VI. CONCLUSION
Fig. 10. Feature relSX(6,3) versus GAG before and after denoising. (a) relSX(6,3) at 40% torque. (b) relSX(6,3) at 100% torque.
successfully extracted in the frequency domain. We will use two such features, the sideband ratio (SBR) and the relative size of sidebands (relSX), to demonstrate the effectiveness of the denoising routine. SBR: Let X denote the number of groups on each side of the dominant component, the SBR is then defined as the ratio between the energy of the RMC and NonRMC components in Dsup formulated as m X h=1 g =−X NonRMC SBR(X) = m X . (23) h=1 g =−X NonRMC + RMC The results of SBR(2) at 40% and 100% torque levels are shown in Fig. 9(a) and (b), respectively. relSX: Let m denote the index of harmonics and X the number of groups on each side of the dominant component, the relSX is then defined as X avg(NonRM C) relSX(m, X) =
g =−X NonRM C+RM C
2X + 1
.
(24)
The results of relSX(6,3) at 40% and 100% torque levels are shown in Fig. 10(a) and (b), respectively. The correlation coefficient and PMD for these two features at different torque levels are summarized in Tables II and III,
The quality of the vibration signals and the features derived from such data is closely related to the success of fault diagnosis and failure prognosis. A blind deconvolution denoising scheme was developed for vibration signals acquired from the gearbox of a helicopter test bed. The denoising scheme uses a recursive filtering technique to minimize a given cost function to improve the SNR of the vibration signals. Then, the qualities of the features extracted from the de-noised data and those from TSA data were compared in terms of accuracy and precision indexes. Experimental results regarding the SNR and performance indexes demonstrate the efficiency of the blind deconvolution denoising routine. Such signal processing techniques promise to significantly lower the fault detection thresholds and improve the performance of diagnositic/prognoistic algorithms, particularly when implemented onboard an aircraft. The blind deconvolution de-nosing scheme developed in this research is a generic one; it can be applied to various systems if nonlinear projections for these systems can be developed. ACKNOWLEDGMENT The authors gratefully acknowledge the support and assistance received from our sponsors. REFERENCES ¨ [1] I. Y. Onel and M. Benbouzid, “Induction motor bearing failure detection and diagnosis: Park and concordia transform approaches comparative study,” IEEE/ASME Trans. Mechatronics, vol. 13, no. 2, pp. 257–262, Apr. 2008. [2] M. Djeziri, R. Merzouki, B. Bouamama, and G. Dauphin-Tanguy, “Robust fault diagnosis by using bond graph approach,” IEEE/ASME Trans. Mechatronics, vol. 12, no. 6, pp. 599–611, Dec. 2007.
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
[3] I. Cade, P. Keogh, and M. Sahinkaya, “Fault identification in rotor/magnetic bearing systems using discrete time wavelet coefficients,” IEEE/ASME Trans. Mechatronics, vol. 10, no. 6, pp. 648–657, Dec. 2005. [4] Y. Izumikawa, K. Yubai, and J. Hirai, “Fault-tolerant control system of flexible arm for sensor fault by using reaction force observer,” IEEE/ASME Trans. Mechatronics, vol. 10, no. 4, pp. 391–396, Aug. 2005. [5] H.-S. Yan and Y.-C. Wu, “A novel design of a brushless dc motor integrated with an embedded planetary gear train,” IEEE/ASME Trans. Mechatronics, vol. 11, no. 5, pp. 551–557, Oct. 2006. [6] P. D. McFadden and J. D. Smith, “An explanation for the asymmetry of the modulation sidebands about the tooth meshing frequency in epicyclic gear vibration,” Proc. Inst. Mech. Eng.—Mech. Eng. Sci., vol. 199, no. 1, pp. 65–70, 1985. [7] J. Keller and P. Grabill, “Vibration monitoring of a UH-60A main transmission planetary carrier fault,” in Proc. Amer. Helicopter Soc. 59th Annu. Forum, Phoenix, AZ, May 2003, pp. 1–11. [8] B. Wu, A. Saxena, T. Khawaja, R. Patrick, G. Vachtsevanos, and R. Sparis, “An approach to fault diagnosis of helicopter planetary gears,” in Proc. AUTOTESTCON, San Antonio, TX, Sep. 2004, pp. 475–481. [9] G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its applications,” Opt. Lett., vol. 13, no. 7, pp. 547–549, Jul. 1988. [10] D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Process. Mag., vol. 13, no. 3, pp. 43–64, May 1996. [11] R. Prost and R. Goutte, “Discrete constrained iterative deconvolution with optimized rate of convergence,” Signal Process., vol. 7, pp. 209–230, Dec. 1984. [12] J. Antoni, “Blind separation of vibration components: Principles and demonstrations,” Mech. Syst. Signal Process., vol. 19, pp. 1166–1180, 2005. [13] A. K. Nandi, D. Mampel, and B. Roscher, “Blind deconvolution of ultrasonic signals in nondestructive testing applications,” IEEE Trans. Signal Process., vol. 45, no. 5, pp. 1382–1390, May 1997. [14] R. Peled, S. Braun, and M. Zacksenhouse, “A blind deconvolution separation of multiple sources, with application to bearing diagnostics,” Mech. Syst. Signal Process., vol. 19, pp. 1181–1195, 2005. [15] G. Gelle, M. Colas, and G. Delaunay, “Blind sources separation applied to rotating machines monitoring by acoustical and vibrations analysis,” Mech. Syst. Signal Process., vol. 14, no. 3, pp. 427–442, 2000. [16] A. Saxena, B. Wu, and G. Vachtsevanos, “A methodology for analyzing vibration data from planetary gear system using complex Morelt wavelets,” in Proc. Amer. Control Conf., vol. 7, Portland, OR, Jun. 2005, pp. 4730– 4735. [17] B. Wu, A. Saxena, R. Romano, and G. Vachtesevanos, “Vibration monitoring for fault diagnosis of helicopter planetary gears,” in Proc. 16th IFAC World Cong., Prague, Czech Republic, Jul. 2005, pp. 1–6.
Bin Zhang (M’04–SM’06) received the B.E. and M.S.E. degrees in mechanical engineering from Nanjing University of Science and Technology, Nanjing, China, in 1993 and 1999, respectively, and the Ph.D. degree in electrical engineering from Nanyang Technological University, Singapore, in 2007. Since 2005, he has been a Postdoctoral Researcher in the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta. He is the author or coauthor of more than 60 technical papers. His current research interests include fault diagnosis and failure prognosis, systems and control, digital signal processing, learning control, intelligent systems and their applications to robotics, power electronics, and various mechanical systems.
565
Taimoor Khawaja is currently working toward the Ph.D. degree in electrical engineering in the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta. His current research interests include design and implementation of real-time architectures for prognostics and health management (PHM) systems, blind deconvolution theory for designing efficient methods for signal denoising and exploring data-driven methods based on support vector machines and Bayesian inference for online fault detection and identification (FDI), and failure prognosis.
Romano Patrick received the degrees in electrical engineering from the University of Texas, Arlington, and the University of Panamericana, Guadalajara, Mexico, and the M.B.A. degree and the Ph.D. degree in electrical and computer engineering from Georgia Institute of Technology, Atlanta, in 2007. He is currently a Project Manager with Impact Technologies, LLC, Rochester, NY. His current research interests include interdisciplinary integration of hardware, software, and techniques to support costeffective design and implementation of engineering systems. He was involved in electronic, mechanical, and software design, stateof-the-art process automation, and novel machine health monitoring for a variety of industrial and government sponsors. He was a Graduate Professor and a Program Coordinator at the University of Panamericana.
George Vachtsevanos (S’62–M’63–SM’89) received the B.E.E. degree in electrical engineering from the City College of New York, New York, NY, in 1962, the M.E.E. degree in electrical engineering from New York University, New York, in 1963, and the Ph.D. degree in electrical engineering from the City University of New York, New York, in 1970. He is currently a Professor Emeritus of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, where he directs the Intelligent Control Systems Laboratory. His work is funded by government agencies and industry. He is the author or coauthor of more than 240 technical papers. Prof. Vachtsevanos was the recipient of the IEEE Control Systems Magazine Outstanding Paper Award for the years 2002 and 2003. He was also awarded the 2002 and 2003 Georgia Tech School of Electrical and Computer Engineering Distinguished Professor Award and the 2003 and 2004 Georgia Institute of Technology Outstanding Interdisciplinary Activities Award.
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
Blind Deconvolution Denoising for Helicopter Vibration Signals Bin Zhang, Senior Member, IEEE, Taimoor Khawaja, Romano Patrick, and George Vachtsevanos, Senior Member, IEEE
Abstract—Critical aircraft assets are required to be available when needed, while exhibiting attributes of reliability, robustness, and high confidence under a variety of flight regimes, and maintained on the basis of their current condition rather than on the basis of scheduled maintenance practices. New and innovative technologies must be developed and implemented to address these concerns. Condition-based maintenance requires that the health of critical components/systems be monitored and diagnostic/prognostic strategies be developed to detect and identify incipient failures and predict the failing component’s remaining useful life. Typically, vibration and other key indicators onboard an aircraft are severely corrupted by noise, thus curtailing the ability to accurately diagnose and predict failures. This paper introduces a novel blind deconvolution denoising scheme that employs a vibration model in the frequency domain and attempts to arrive at the true vibration signal through an iterative optimization process. Performance indexes are defined and data from a helicopter are used to demonstrate the effectiveness of the proposed approach.
Fig. 1.
Crack of planetary gear carrier plate of the UH-60A helicopter.
Index Terms—Blind deconvolution, planetary gear train, vibration signal denoising.
I. INTRODUCTION AULT DETECTION, isolation, and identification, especially for safety critical components and subsystems, have recently drawn increasing interest in the condition-based maintenance (CBM) community [1]–[4]. Epicyclic, or planetary, gear trains are widely used in the main transmission of many systems, such as the helicopter and other aircraft [5], [7]. This kind of gear system consists of an inner sun gear, two or more rotating planet gears, a stationary outer ring gear, and a planetary carrier. The planet gears, which surround the sun gear, are riding on the planetary carrier and also rotating inside the outer ring gear. In the operation, torque is transmitted through the sun gear to the planets and planetary carrier. The carrier, in turn, transmits torque to the main rotor shaft and blades [6]–[9]. The gearbox is a critical component that directly impacts the safety and performance of the aircraft.
F
Manuscript received October 8, 2006; revised April 25, 2008. Current version published October 8, 2008. Recommended by Technical Editor M.-Y. Chow. This work was supported in part by the Defense Advanced Research Projects Agency (DARPA) and in part by Northrop Grumman under DARPA Prognosis Program Contract HR0011-04-C-003. B. Zhang and T. Khawaja are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail: [email protected]). R. Patrick is with Impact Technologies, LLC, Rochester, NY 14623 USA (e-mail: [email protected]). G. Vachtsevanosis is with Georgia Institute of Technology, Atlanta, GA 30332 USA, and also with Impact Technologies, LLC, Rochester, NY 14623 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMECH.2008.2002324
The UH-60A Blackhawk and Seahawk main transmission employs a five-planet-gear epicyclic gear system [7]. Recently, a crack in the planetary carrier was found, as shown in Fig. 1. This resulted in flight restrictions on many helicopters. Manual inspection of all transmissions is not only costly, but also time prohibitive. A CBM based on vibration signal analysis is a cost-effective solution [16]. To measure the vibration signals, the transducer is mounted at a fixed point on the frame of the gearbox. If a fault occurs on the gear system, it is expected that the vibration signal will exhibit a characteristic signature that reveals the presence, severity, and location of the fault. In fault diagnosis and failure prognosis, gear health is typically described in terms of feature or condition indicator (CI) values, which are extracted from the vibration signals [7]. The success of diagnosis and prediction of the remaining useful life of a gear system highly depends on the quality of these features. However, during the course of a flight, the vibration signals derived from accelerometers often have multiple excitation sources, and in most cases, information about the noise characteristics is not available. The measured output signal is, therefore, the combined effect from these different sources. Noise signals will degrade the quality of the features, and consequently, the performance of diagnosis and prognosis, such as the detection threshold, the accuracy of remaining useful life of components, etc. For example, when a fault is at its early stage, the indication in the vibration signals and the features is often masked by noise. It is, therefore, important to develop an efficient and reliable denoising algorithm to remove noise signals as much as possible and make the fault characteristics perceptible in the extracted features.
1083-4435/$25.00 © 2008 IEEE Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
Fig. 2.
559
Blind deconvolution denoising scheme.
The widely used denoising technique in planetary gears is time-synchronous averaging (TSA) [7]–[9]. To realize TSA, a pulse signal that is synchronized to the rotation of a gear indicating the start of an individual revolution is assumed to be available. The sampling data collected from different revolutions are interpolated to the same length, and then, ensemble averaging is performed. The interpolation operation is often computationally intensive. A TSA technique in the frequency domain has been proposed to reduce the computation drastically [17]. Some features are successfully extracted from the TSA data and used in the fault diagnosis and prognosis [8], [16], [17]. The TSA technique enhances the components at the frequencies that are multiples of the shaft frequency, which are often related to the meshing of the gear teeth [8], [16]. At the same time, it tends to average out the external random disturbances and noise that are asynchronous with the rotation of the gear. However, in the frequency domain, it is possible that some noise components appear at frequencies that are multiples of the shaft frequency and enhanced by TSA processing. Hence, it is necessary to further denoise the TSA data to achieve better SNR and improve the quality of the features extracted from the TSA data. The vibration signal collected from the transducer on the frame of the gearbox is amplitude modulated and corrupted by additive noise. Blind deconvolution has shown very promising results for similarly formulated problems in image and speech processing [10], [11], but it is rarely applied to vibration signals [12]–[15]. In this paper, a blind deconvolution scheme, which consists of vibration analysis, nonlinear projection, and cost function optimization, is employed to the denoising of the vibration signals from a UH-60A helicopter gearbox, whose planetary carrier has a seeded crack fault growing with the evolving operation of the system. The analysis of the system, the theoretical basis of the blind deconvolution algorithm, and the proposed performance indexes are addressed in detail. Experimental results show the effectiveness of the proposed denoising scheme. II. DENOISING SCHEME ARCHITECTURE The architecture of the proposed blind deconvolution denoising scheme is shown in Fig. 2. In this scheme, nonlinear projection, which is based on vibration analysis in the frequency domain, and cost function minimization are critical components, which will be described in later sections. The proposed denoising approach first transforms the measured vibration signal s(t) ¯ ), which is an estimate of to S(f ). Then, an inverse filter Z(f
Fig. 3.
Configuration of an epicyclic gear system.
the inverse of the modulating signal in the frequency domain, is convoluted with the measured vibration signal S(f ) to demodulate S(f ) and give a rough estimate of the noise-free vibration ¯ ). The signal B(f ¯ ) passes through the nonlinear signal B(f ¯ ) to a subspace that contains only projection, which maps B(f known characteristics of the vibration signal, to yield Bn l (f ). ¯ ) and Bn l (f ) is denoted as E(f ). The difference between B(f ¯ ) iteratively to minimize the E(f ) and when By adjusting Z(f ¯ ) → B(f ) can be E(f ) reaches a minimal value, the signal B(f regarded as the denoised vibration signal. Through an inverse Fourier transform, the denoised vibration signal in the time domain can be obtained as well.
III. SYSTEM DESCRIPTION The vibration signals are derived from the main transmission gearbox of the UH-60A helicopter. The gearbox is an epicyclic gear system with five planet gears, whose configuration is shown in Fig. 3. The following analysis assumes an ideal system. The accelerometer is mounted at a fixed point on the annulus gear at the position θ = 0. Since the vibration signal is generated from the meshing of gear teeth and the planetary gears are rotating inside the angular gear, the vibration signal is amplitude modulated to the static accelerometer. The observed vibration amplitude will be large when the planetary gear is close to the accelerometer, and it will be small when the planetary gear is far. Suppose that there is only one planetary gear, then the vibration observed by the transducer should have the largest amplitude when the planetary gear is at θ = 0, 2π, 4π, . . .. Similarly, the vibration should have the smallest amplitude at θ = π, 3π, . . ., i.e., the teeth meshing vibration signal is amplitude modulated as the planetary carrier rotates. In the ideal case, the Np = 5 planetary gears are evenly spaced and suppose that the planetary carrier has a rotation frequency of f s . Then, the planetary gear p at time instant t has a phase of p−1 θp = 2π f s t + . Np
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
(1)
560
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
In this case, the amplitude-modulating signal for planetary gear p can be written in the time domain as ap (t) =
N
αn cos(nθp )
(2)
n =−N
where N is the number of sidebands about the harmonics under consideration and αn is the amplitude of the component of the modulating signal at frequency nf s . The “sidebands” are defined as the frequency components that appear as a harmonically spaced series [6]. In the ideal case, it is natural to assume that all mesh vibrations generated from different planetary gears are of the same amplitude but of different phases. The annulus gear has Nt = 228 teeth. Since the speed at which teeth meshing takes place is proportional to the angular velocity of the planetary carrier, the meshing vibration appears at frequencies Nt f s [6]. In addition, suppose that, in the frequency domain, the meshing vibration signal has amplitudes βm at its harmonics, then the vibration signal generated from planet gear p can be written in the form of bp (t) =
M
βm sin(mNt θp )
Fig. 4. Superposition of vibration from different gears with N p = 5. (a) mN t + n = kN p . (b) mN t + n = kN p .
(3)
m =1
where M is the number of harmonics under consideration. Then, to the static accelerometer, the vibration signal of planetary gear p appears as the product of the meshing vibration signal and the amplitude-modulating signal. It is denoted as yp (t) and is given as N M 1 αn βm sin((mNt + n)θp ). 2 m =1 n =−N (4) When there are more than one, say Np , planetary gears, the vibration signal observed by the accelerometer is the superposition of the Np vibration signals generated from Np different planetary gears. This superposition vibration signal has the form
Fig. 5.
Spectrum of combined signal with N t = 228 and N p = 5.
yp (t) = ap (t)bp (t) =
Np M N 1 αn βm sin((mNt + n)θp ) y(t) = 2 p=1 m =1 n =−N
N 1 mNt + n = αn βm sin 2π(p − 1) (5) 2 p=1 m =1 Np Np
M
n =−N
where (1) and the fact that sin(2kπ + θ) = sin(θ) for any integer k are used. Since the planetary gears are evenly spaced, the phase angel of the sidebands will be evenly spaced along 2π [6]. From (5), it is obvious that if sideband mNt + n is not a multiple of Np and (mNt + n)/Np has a remainder of γ, the vibration components from different gears are evenly spaced by the angle 2γπ/Np . In this case, when the vibrations generated from different planetary gears are combined, these sidebands add destructively and become zero, as illustrated in Fig. 4(a), in which φp,m ,n with 1 ≤ p ≤ 5 indicates the frequency components of gear p. Therefore, the frequency components appear at sidebands where mNt + n = kNp are termed as nonregular
meshing components (NonRMC). On the contrary, if sideband mNt + n is a multiple of Np , the remainder of (mNt + n)/Np will be zero. In this case, the vibration components from different gears do not have a phase difference. When the vibration signals from different planetary gears are combined, these sidebands add constructively and are reinforced, as illustrated in Fig. 4(b). These frequency components that appear at sidebands where mNt + n = kNp are termed as regular meshing components (RMC) or apparent sidebands. This process of frequency components adding destrutively/constructively finally generates asymmetrical sidebands. A partial frequency spectrum of the sidebands about the firstorder harmonic that illustrates this asymmetry is shown in Fig. 5, where Np = 5 and Nt = 228. Note that the peak of the spectrum does not appear at n = 0 (order 228) but at n = −3 (order 225), n = 2 (order 230), etc. The largest spectral amplitude (also known as dominant sideband) appears at the frequency closest to the gear meshing frequency. The vibration analysis in the frequency domain and previous research results [6], [7] suggest that, for an ideal system, only terms at frequencies that are multiples of the number of planetary gears survive. Then, the Fourier transform of the vibration data can be written as Y (f ) = fn l (γm ,n (mNt + n)f s )
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
(6)
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
where γm ,n is the magnitude of the spectral amplitude at (mNt + n)f s and fn l is the nonlinear projection 1, if mNt + n is a multiple of Np fn l = (7) 0, otherwise. From the previous analysis of the system vibration behavior, the following assumptions are made and used in the blinddeconvolution scheme. 1) The sidebands of the vibration spectra exist only at frequencies (with a resolution of the planetary carrier rotation frequency f s ) that are multiples of the number of planetary gears. 2) The phase of each sideband depends on the geometrical angular position of the planet and is given by θp (mNt + n) for sideband (m, n). 3) The amplitude of the modulating signal a(t) decreases monotonically on either side of the maximum until it reaches the minimum. Note that this assumption is somewhat restrictive. IV. BLIND DECONVOLUTION DENOISING SCHEME
(8)
where s(t) is the measured vibration signal, b(t) is the noise-free unmodulated vibration signal, a(t) is the modulating signal, and n(t) is the cumulative additive noise. Note that the modulating signal a(t) is itself affected by noise in the system. Let a ˆ(t) denote the ideal or noise-free modulating signal and na (t) the noise introduced in this signal. Consequently, a(t) can be represented as a(t) = a ˆ(t) + na (t).
(9)
Thus, (8) can be rewritten as follows: s(t) = (ˆ a(t) + na (t))b(t) + n(t) =a ˆ(t)b(t) + n ˆ (t)
(10)
where n ˆ (t) = na (t)b(t) + n(t) contains the total additive noise in the system. On the other hand, the factor a ˆ(t) describes the multiplicative noise in the system. The goal for a denoising scheme, such as the one described here, is to recover the unknown vibration signal b(t) from the observed signal s(t) given partial information about the noise sources and characteristics of the vibration signal. A typical approach would be to find the inverse of a ˆ(t) zˆ(t) =
1 a ˆ(t)
(11)
such that
¯b(t) = s(t)¯ z (t).
(13)
If partial knowledge about how the plate system is influenced by the modulating signal a ˆ(t) and a reasonable understanding of the true vibration signal are available, then the ideal characteristics of the vibration signal can be obtained by projecting this estimated signal ¯b(t) into a subspace with only the known ideal characteristics of the vibration signal to yield a refined signal bn l (t). Since this nonlinear projection, as the subscript signifies, removes all uncharacterisitic components that exist in the rough estimate ¯b(t), it is necessary to stress the importance of a good understanding of the underlying process. An iterative scheme then refines these results by minimizing the error between the two signals ¯b(t) and bn l (t), i.e., (14)
Previous research results detail the spectral characteristics of vibration signals for rotating equipment [6], [7]. It is, therefore, appropriate to investigate the measured noisy vibration in the frequency domain. The convolution theorem states that the product of two signals in the time domain is equivalent to their convolution in the frequency domain. Thus, model (10) can be written in the frequency domain as ˆ ) ∗ B(f ) + N ˆ (f ) S(f ) = A(f
(15)
ˆ ), B(f ), with ∗ being the convolution operator and S(f ), A(f ˆ and N (f ) are the Fourier transforms of s(t), a ˆ(t), b(t), and n ˆ (t), respectively. Then, the goal in the frequency domain is to recover B(f ). Writing (13) in the frequency domain, we have ¯ ) = S(f ) ∗ Z(f ¯ ) B(f
(16)
¯ ) and Z(f ¯ ) being the Fourier transforms of ¯b(t) and with B(f ¯ ) through the nonlinear projecz¯(t), respectively. Passing B(f tion, it yields Bn l (f ). Then, in the frequency domain, we will ¯ ) minimize the difference between Bn l (f ) and B(f ¯ ) − Bn l (f ). min E(f ) = min B(f
(17)
¯ ) to minimize (17). When it The iterative process refines Z(f ¯ reaches the minimal value, Z(f ) converges to Z(f ). Then, with ¯ ) in (16), a good estimate for B(f ) is this Z(f ) replacing Z(f obtained as B(f ) = S(f ) ∗ Z(f ).
(18)
Last, the estimate is transformed back into the time domain to recover the noise-free vibration signal b(t) ∞ ei2π f t B(f ) df . (19) b(t) = F −1 B(f ) = −∞
b(t) = (s(t) − n ˆ (t))ˆ z (t) = s(t)ˆ z (t) − n ˆ (t)ˆ z (t).
Note, however, that little can be assumed about n ˆ (t), a ˆ(t) is not available, and zˆ(t) is not applicable. To solve this problem, rather than using zˆ(t), we propose an iterative denoising scheme that starts with z¯(t), a very rough initial estimate of the inverse of the modulating signal, which demodulates the observed signal s(t) to give a rough noise-free vibration signal ¯b(t)
min e(t) = min ¯b(t) − bn l (t).
From the vibration analysis of the gearbox, we know that the vibration signals collected from the transducer are amplitude modulated [6]. Multiple sources of noise may further corrupt the signal. A simplified model for such a complex signal may be defined as s(t) = a(t)b(t) + n(t)
561
(12)
To solve this problem, the following additional assumptions are made (continued from the previous section).
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
562
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
TABLE I GROUND TRUTH DATA OF CRACK LENGTH (IN INCHES)
1) Z(f ) exists and Z(f ) < ∞. 2) Since the modulating signal a ˆ(t) is always positive, its inverse should be positive too. Hence, the Fourier transform of its inverse Z(f ) contains a dc component. With the previous two assumptions being taken into consideration, the cost function is defined as J=
[B(f ) − Bn l (f )]2 +
2 Z(f ) − 1
Fig. 6.
Growth of crack length.
Fig. 7.
(a) SNR at 40% torque. (b) SNR at 100% torque.
(20)
f ∈D s u p
where Dsup is the frequency range that contains the main vibration information. Because of the periodic fade of signal spectra between harmonics, a window centered at harmonic frequencies is used to define critical frequencies. All these windows form the support Dsup . In (20), assumptions 1) and 2) are used to arrive at the second term to avoid an all-zero inverse filter Z(f ), which leads to the trivial solution for error minimization. Moreover, an iterative optimization routine is required to implement this scheme. The iterative conjugate gradient method is called upon to address the optimization problem. This method has faster convergence rate in general as compared to the steepest descent method [11]. V. EXPERIMENTAL RESULTS The planetary carrier has a seeded crack with an initial length of 1.344 in. The crack, as shown in Fig. 1, grows with the evolving operation of the gearbox. The gearbox is installed on a test bed and operates over a large number of ground-air-ground (GAG) cycles at different torque levels. This way, vibration data are acquired at different torque levels and different crack lengths. In each GAG cycle, the torque increases from 20% to 40%, and finally to 100%, and then, decreases to 20% for the next cycle. A. Actual Crack Growth The ground truth crack length data at discrete GAG cycles are available and tabulated in Table I. With these data, the crack length for GAG cycles from 1 to 1000 can be obtained via interpolation, which results in the crack length growth curve shown in Fig. 6. B. Signal-to-Noise Ratio The SNR, before and after denoising, is investigated, and the results at 40% and 100% torque levels are shown in Fig. 7(a) and (b), respectively. The improvement in the SNR value is significant.
C. Feature Performance Indexes Although the blind deconvolution routine shows a significant improvement in the SNR, it is desirable that it improves also the quality of the features or condition indicators. The accuracy and precision of mappings between the evolution of features and the actual crack growth have an important impact on the performance of diagnostic and prognostic algorithms and the performance of the CBM system overall. To evaluate the quality of the features, two performance indexes are introduced. The first performance index is an accuracy measure defined as the linear correlation between the feature values and the crack length growth along the GAG cycle axis. Suppose that x is the
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
feature vector and y the crack growth curve with x ¯ and y¯ being their means, respectively. The correlation coefficient between them is ss2xy (21) CCR(x, y) = ssxx ssy y where ss ¯)(yi − y¯), ssxx = (xi − x ¯)2 , and i −x xy = (x 2 ssy y = (yi − y¯) , respectively. The extracted feature will be used to map and predict the crack growth. Hence, a high correlation coefficient is expected to generate an accurate estimate of actual crack length and is preferred. Then, for a good feature, the value of the correlation coefficient should be near to 1. Because of changes in operating conditions and disturbances, the feature values along a GAG cycle are often very noisy and need to be smoothed through a low-pass filtering operation. This filtering will not change the characteristics of the feature curve. In this case, the linear correlation is calculated based on the smoothed feature curve. The second performance index is a precision measure corresponding to a normalized measure of the signal dispersion. It is termed as percent mean deviation (PMD) and defined as n x i −˜x i PMD(x, x ˜) =
i=1
lx
x ˜i
× 100
Fig. 8.
563
Groups and D su p at the first harmonic with N t = 228 and N p = 5.
(22)
˜ where lx is the number of entities in the feature vector x and x is the smoothed feature vector. In the fault detection and prognosis algorithms, the extracted feature values are used as a measurement input. Therefore, this performance index is closely related to the detection threshold and precision of the prognosis or prediction of the remaining useful life. From a precise feature, the fault detection and prognosis algorithms can detect the incipient failure and/or predict the fault growth with a high confidence. PMD for a precise feature should have a very small value close to 0. D. Experimental Results As the previous analysis shows, only RMC, which include the dominant and apparent components appearing at frequencies that are multiples of Np , will survive in the vibration signal. However, the actual system is not ideal, particularly under faulty conditions. Then, the NonRMC, which appear at frequencies that are not multiples of Np , will not be zero. Hence, the frequency spectra are divided into two categories: RMC and NonRMC. To simplify the description and relevant computations, Dsup is defined in terms of groups of frequency spectra. From the previous spectra analysis, for a gearbox with Np planetary gears, the RMCs are evenly spaced by Np − 1 sidebands. Then, it is natural to define the group size as Np . For the five-planetarygear gearbox in this research, the group size is defined as 5. Then, Dsup is given as the number of groups on each side of the dominant component. The concepts of “group” and “Dsup ” are illustrated in Fig. 8, where only the groups and part of Dsup at the first harmonic are shown. The entire Dsup is the combination of these groups at all the harmonics under consideration.
Fig. 9. Feature SBR(2) versus GAG before and after denoising. (a) SBR(2) at 40% torque. (b) SBR(2) at 100% torque.
The relative size of the NonRMC sidebands to the RMC sidebands in Dsup or part of Dsup may indicate the presence of a fault on the planetary carrier. As the crack grows, the operating condition deviates from the ideal one. In this case, the constructive superposition of vibration signals from different planetary gears is attenuated while the destructive superposition is reinforced. This phenomenon results in smaller RMC and larger NonRMC sidebands. This indicates that the ratio between RMC and NonRMC sidebands can be used as a characteristic feature. Based on this notion, features or condition indicators are
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
564
IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 13, NO. 5, OCTOBER 2008
TABLE II CORRELATION COEFFICIENTS OF THE FEATURES
TABLE III PERCENT MEAN DEVIATION OF THE FEATURES
respectively. In the tables, D-N simply means de-noised. It is clear that for most cases, substantial feature improvements are achieved via the application of the denoising routine. It is worth noting that the influence of crack on vibration spectra at low torque level 20% is not as remarkable as those at high torque levels. Since the nonlinear projection in the developed denoising routine depends on the analysis of vibration spectra, the application of this routine shows more reliable performance on high torque levels. VI. CONCLUSION
Fig. 10. Feature relSX(6,3) versus GAG before and after denoising. (a) relSX(6,3) at 40% torque. (b) relSX(6,3) at 100% torque.
successfully extracted in the frequency domain. We will use two such features, the sideband ratio (SBR) and the relative size of sidebands (relSX), to demonstrate the effectiveness of the denoising routine. SBR: Let X denote the number of groups on each side of the dominant component, the SBR is then defined as the ratio between the energy of the RMC and NonRMC components in Dsup formulated as m X h=1 g =−X NonRMC SBR(X) = m X . (23) h=1 g =−X NonRMC + RMC The results of SBR(2) at 40% and 100% torque levels are shown in Fig. 9(a) and (b), respectively. relSX: Let m denote the index of harmonics and X the number of groups on each side of the dominant component, the relSX is then defined as X avg(NonRM C) relSX(m, X) =
g =−X NonRM C+RM C
2X + 1
.
(24)
The results of relSX(6,3) at 40% and 100% torque levels are shown in Fig. 10(a) and (b), respectively. The correlation coefficient and PMD for these two features at different torque levels are summarized in Tables II and III,
The quality of the vibration signals and the features derived from such data is closely related to the success of fault diagnosis and failure prognosis. A blind deconvolution denoising scheme was developed for vibration signals acquired from the gearbox of a helicopter test bed. The denoising scheme uses a recursive filtering technique to minimize a given cost function to improve the SNR of the vibration signals. Then, the qualities of the features extracted from the de-noised data and those from TSA data were compared in terms of accuracy and precision indexes. Experimental results regarding the SNR and performance indexes demonstrate the efficiency of the blind deconvolution denoising routine. Such signal processing techniques promise to significantly lower the fault detection thresholds and improve the performance of diagnositic/prognoistic algorithms, particularly when implemented onboard an aircraft. The blind deconvolution de-nosing scheme developed in this research is a generic one; it can be applied to various systems if nonlinear projections for these systems can be developed. ACKNOWLEDGMENT The authors gratefully acknowledge the support and assistance received from our sponsors. REFERENCES ¨ [1] I. Y. Onel and M. Benbouzid, “Induction motor bearing failure detection and diagnosis: Park and concordia transform approaches comparative study,” IEEE/ASME Trans. Mechatronics, vol. 13, no. 2, pp. 257–262, Apr. 2008. [2] M. Djeziri, R. Merzouki, B. Bouamama, and G. Dauphin-Tanguy, “Robust fault diagnosis by using bond graph approach,” IEEE/ASME Trans. Mechatronics, vol. 12, no. 6, pp. 599–611, Dec. 2007.
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
ZHANG et al.: BLIND DECONVOLUTION DENOISING FOR HELICOPTER VIBRATION SIGNALS
[3] I. Cade, P. Keogh, and M. Sahinkaya, “Fault identification in rotor/magnetic bearing systems using discrete time wavelet coefficients,” IEEE/ASME Trans. Mechatronics, vol. 10, no. 6, pp. 648–657, Dec. 2005. [4] Y. Izumikawa, K. Yubai, and J. Hirai, “Fault-tolerant control system of flexible arm for sensor fault by using reaction force observer,” IEEE/ASME Trans. Mechatronics, vol. 10, no. 4, pp. 391–396, Aug. 2005. [5] H.-S. Yan and Y.-C. Wu, “A novel design of a brushless dc motor integrated with an embedded planetary gear train,” IEEE/ASME Trans. Mechatronics, vol. 11, no. 5, pp. 551–557, Oct. 2006. [6] P. D. McFadden and J. D. Smith, “An explanation for the asymmetry of the modulation sidebands about the tooth meshing frequency in epicyclic gear vibration,” Proc. Inst. Mech. Eng.—Mech. Eng. Sci., vol. 199, no. 1, pp. 65–70, 1985. [7] J. Keller and P. Grabill, “Vibration monitoring of a UH-60A main transmission planetary carrier fault,” in Proc. Amer. Helicopter Soc. 59th Annu. Forum, Phoenix, AZ, May 2003, pp. 1–11. [8] B. Wu, A. Saxena, T. Khawaja, R. Patrick, G. Vachtsevanos, and R. Sparis, “An approach to fault diagnosis of helicopter planetary gears,” in Proc. AUTOTESTCON, San Antonio, TX, Sep. 2004, pp. 475–481. [9] G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its applications,” Opt. Lett., vol. 13, no. 7, pp. 547–549, Jul. 1988. [10] D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Process. Mag., vol. 13, no. 3, pp. 43–64, May 1996. [11] R. Prost and R. Goutte, “Discrete constrained iterative deconvolution with optimized rate of convergence,” Signal Process., vol. 7, pp. 209–230, Dec. 1984. [12] J. Antoni, “Blind separation of vibration components: Principles and demonstrations,” Mech. Syst. Signal Process., vol. 19, pp. 1166–1180, 2005. [13] A. K. Nandi, D. Mampel, and B. Roscher, “Blind deconvolution of ultrasonic signals in nondestructive testing applications,” IEEE Trans. Signal Process., vol. 45, no. 5, pp. 1382–1390, May 1997. [14] R. Peled, S. Braun, and M. Zacksenhouse, “A blind deconvolution separation of multiple sources, with application to bearing diagnostics,” Mech. Syst. Signal Process., vol. 19, pp. 1181–1195, 2005. [15] G. Gelle, M. Colas, and G. Delaunay, “Blind sources separation applied to rotating machines monitoring by acoustical and vibrations analysis,” Mech. Syst. Signal Process., vol. 14, no. 3, pp. 427–442, 2000. [16] A. Saxena, B. Wu, and G. Vachtsevanos, “A methodology for analyzing vibration data from planetary gear system using complex Morelt wavelets,” in Proc. Amer. Control Conf., vol. 7, Portland, OR, Jun. 2005, pp. 4730– 4735. [17] B. Wu, A. Saxena, R. Romano, and G. Vachtesevanos, “Vibration monitoring for fault diagnosis of helicopter planetary gears,” in Proc. 16th IFAC World Cong., Prague, Czech Republic, Jul. 2005, pp. 1–6.
Bin Zhang (M’04–SM’06) received the B.E. and M.S.E. degrees in mechanical engineering from Nanjing University of Science and Technology, Nanjing, China, in 1993 and 1999, respectively, and the Ph.D. degree in electrical engineering from Nanyang Technological University, Singapore, in 2007. Since 2005, he has been a Postdoctoral Researcher in the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta. He is the author or coauthor of more than 60 technical papers. His current research interests include fault diagnosis and failure prognosis, systems and control, digital signal processing, learning control, intelligent systems and their applications to robotics, power electronics, and various mechanical systems.
565
Taimoor Khawaja is currently working toward the Ph.D. degree in electrical engineering in the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta. His current research interests include design and implementation of real-time architectures for prognostics and health management (PHM) systems, blind deconvolution theory for designing efficient methods for signal denoising and exploring data-driven methods based on support vector machines and Bayesian inference for online fault detection and identification (FDI), and failure prognosis.
Romano Patrick received the degrees in electrical engineering from the University of Texas, Arlington, and the University of Panamericana, Guadalajara, Mexico, and the M.B.A. degree and the Ph.D. degree in electrical and computer engineering from Georgia Institute of Technology, Atlanta, in 2007. He is currently a Project Manager with Impact Technologies, LLC, Rochester, NY. His current research interests include interdisciplinary integration of hardware, software, and techniques to support costeffective design and implementation of engineering systems. He was involved in electronic, mechanical, and software design, stateof-the-art process automation, and novel machine health monitoring for a variety of industrial and government sponsors. He was a Graduate Professor and a Program Coordinator at the University of Panamericana.
George Vachtsevanos (S’62–M’63–SM’89) received the B.E.E. degree in electrical engineering from the City College of New York, New York, NY, in 1962, the M.E.E. degree in electrical engineering from New York University, New York, in 1963, and the Ph.D. degree in electrical engineering from the City University of New York, New York, in 1970. He is currently a Professor Emeritus of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, where he directs the Intelligent Control Systems Laboratory. His work is funded by government agencies and industry. He is the author or coauthor of more than 240 technical papers. Prof. Vachtsevanos was the recipient of the IEEE Control Systems Magazine Outstanding Paper Award for the years 2002 and 2003. He was also awarded the 2002 and 2003 Georgia Tech School of Electrical and Computer Engineering Distinguished Professor Award and the 2003 and 2004 Georgia Institute of Technology Outstanding Interdisciplinary Activities Award.
Authorized licensed use limited to: IEEE Xplore. Downloaded on October 29, 2008 at 12:33 from IEEE Xplore. Restrictions apply.
