METHOD AND APPARATUS FOR DETERMINING THE HEALTH OF A COMPONENT USING CONDITION INDICATORS

METHODES ET DISPOSITIF PERMETTANT DE DETERMINER LA CONDITION D'UN COMPOSANT AU MOYEN D'INDICATEURS D'ETAT

English Abstract

Disclosed are techniques used in connection with determining a health
indicator (HI) of a component, such as that of an aircraft component. The HI
is determined using condition indicators (CIs) which parameterize
characteristics about a component minimizing possibility of a false alarm.
Different algorithms are disclosed which may be used in determining one or
more CIs. The HI may be determined using a normalized CI value. Techniques are
also described in connection with selecting particular CIs that provide for
maximizing separation between HI classifications. Given a particular HI at a
point in time for a component, techniques are described for predicting a
future state or health of the component using the Kalman filter. Techniques
are described for estimating data values as an alternative to performing data
acquisitions, as may be used when there is no pre-existing data.

French Abstract

Cette invention concerne des techniques s'utilisant pour déterminer un indicateur de condition (ou de <= santé >=)(HI) d'un composant, tel qu'un composant d'aéronef. A cette fin, on utilise des indicateurs d'état (CI) qui paramètrent des caractéristiques sur un composant dans le but de réduire les risques d'alarme erronée. L'invention concerne divers algorithmes servant à déterminer un ou plusieurs CI. L'indicateur de condition HI peut être déterminé au moyen d'une valeur CI normalisée. Sont également décrites des techniques en rapport avec le choix de CI qui permettent de maximiser la séparation entre classifications HI. A partir d'un HI donné à un moment particulier pour un composant, les techniques décrites permettant de prévoir un état futur ou un niveau de santé du composant au moyen d'un filtre Kalman. Sont décrites des techniques permettant, en lieu et place d'acquisition de données, d'estimer des valeurs de données, comme cela se fait en l'absence de données préexistantes.

Claims

What is Claimed is:
1. A method executed in a computer system for determining a health indicator
associated with a component comprising:
determining a plurality of health classifications;
determining at least one condition indicator quantifying a characteristic of
said
component;
determining a probability associated with each of said health classifications,
said
probability being an estimation that said component is of a particular health
classification
given said at least one condition indicator; and
determining, for a given set of observed values, which of said plurality of
health
classifications is associated with said component using said probabilities
associated with said
health classifications.
2. The method of Claim 1, wherein said determining at least one condition
indicator,
said determining a probability, and said determining which of said health
classifications is
associated with said component, are performed upon two of said plurality of
health
classifications, and
wherein, if a determination is made that a first of said health
classifications is not
associated with said component, repeating said determining at least one
condition indicator,
said determining a probability, and said determining which of said health
classifications is
associated with said component using two other health classifications.
79

3. The method of Claim 2, further comprising:
storing previously observed and calculated data associated with said
component; and
using a portion of said previously observed and calculated data to determine
said
probability associated with said each health classification.
4. The method of Claim 2, further comprising:
obtaining data estimates using models of said component; and
using said data estimates to determine said probability associated with said
each
health classification.
5. The method of Claim l, further comprising:
determining a threshold using a ratio of a portion of said probabilities;
calculating a current quantity using said observed values; and
comparing said current quantity to said ratio to determine which of said
health
classifications is associated with said component.
80

6. The method of Claim 5, wherein said calculating said current quantity
further
comprises:
determining a transformation matrix maximizing a distance between two of said
plurality of health classifications;
determining a covariance matrix for each of said health classifications using
said at
least one condition indicator and said observed data; and
determining said current quantity in accordance with said covariance matrices,
said
transformation matrix and said observed data.
7. The method of Claim 6, wherein, for two of said plurality of health
classifications,
said current quantity, h(X), may be represented as:
<IMG>
wherein,
X represents said at least one conditional indicator forming a vector of
individual
conditional indicator values;
Y represents ~-1/2.PHI. T X ;
~ represents (I - K-1) -1, diagonal matrix of eigenvalues of X using a
characteristic
equation of .SIGMA..PHI. = .PHI.~, .PHI.T.PHI. = identity matrix;
.PHI. represents an nxn matrix of eigenvectors, .PHI.1.. .PHI.n, of X using
said characteristic
equation;
I represents A T.SIGMA.1 A;
K represents A T.SIGMA.2A;
L represents A T(M2 - M1);
81

P2 represents an aposteriori probability of a first of said two health
classifications
given X;
P1 represents an aposteriori probability of a second of said two health
classifications
given X;
A represents said transformation matrix maximizing a distance between
distributions
of said two health classifications represented as ~-1/2.PHI.T;
.SIGMA. i represents a covariance matrix of said first health classification
having a vector of
expected values M1; and
.SIGMA.2 represents a covariance matrix of said second health classification
having a vector
of expected values M2.
8. The method of Claim 1, wherein said at least one condition indicator
indicates a
physical state of a portion of said component.
82

9. A method for determining a health status of a component comprising:
selecting a plurality of condition indicators having a value and each having a
corresponding weighting factor, and at least one threshold value defining at
least two
classifications;
determining a contribution to a health indicator for each of said condition
indicators,
wherein said determining further comprises, for each of said plurality of
indicators:
determining which of said at least two classifications said value of said each
indicator belongs; and
determining said contribution to said health indicator by said each condition
indicator in accordance with a selected one of said at least two
classifications and said
weighted value; and
determining said health indicator in accordance with all contributions by each
of said
condition indicator values.
10. The method of Claim 9, wherein three classifications are associated with
each of
said condition indicators formed by two threshold values, and the method
further comprises:
determining that said contribution to said health indicator for a first of
said condition
indicators is zero if said value of said condition indicator is in said first
classification;
determining that said contribution to said health indicator is a first
multiple of said
weighting factor if said value of said condition indicator is in said second
classification; and
determining that said contribution to said health indicator is a second
multiple of said
weighting factor if said value of said condition indicator is in said third
classification.
11. The method of Claim 10, where the two threshold values are alarm level and
warning level.
83

12. A computer program product for determining a health indicator associated
with a
component comprising machine executable code for
determining a plurality of health classifications;
determining at least one condition indicator quantifying a characteristic of
said
component;
determining a probability associated with each of said health classifications,
said
probability being an estimation that said component is of a particular health
classification
given said at least one condition indicator; and
determining, for a given set of observed values, which of said plurality of
health
classifications is associated with said component using said probabilities
associated with said
health classifications.
13. The computer program product of Claim 12, wherein machine executable code
for determining at least one condition indicator, determining a probability,
and determining
which of said health classifications is associated with said component, are
performed upon
two of said plurality of health classifications, and the computer program
product further
includes machine executable code, wherein, if a determination is made that a
first of said
health classifications is not associated with said component, for repeatedly
determining at
least one condition indicator, repeatedly determining a probability, and
repeatedly
determining which of said health classifications is associated with said
component using two
other health classifications.
14. The computer program product of Claim 13, further comprising machine
executable code for:
storing previously observed and calculated data associated with said
component; and
84

using a portion of said previously observed and calculated data to determine
said
probability associated with said each health classification.
15. The computer program product of Claim 13, further comprising machine
executable code for:
obtaining data estimates using models of said component; and
using said data estimates to determine said probability associated with said
each
health classification.
16. The computer program product of Claim 12, further comprising machine
executable code for:
determining a threshold using a ratio of a portion of said probabilities;
calculating a current quantity using said observed values; and
comparing said current quantity to said ratio to determine which of said
health
classifications is associated with said component.
17. The computer program product of Claim 16, wherein said machine executable
code for calculating said current quantity further comprises machine
executable code for:
determining a transformation matrix maximizing a distance between two of said
plurality of health classifications;
determining a covariance matrix for each of said health classifications using
said at
least one condition indicator and said observed data; and
determining said current quantity in accordance with said covariance matrices,
said
transformation matrix and said observed data.
85

18. The computer program product of Claim 17, wherein, for two of said
plurality of
health classifications, said current quantity, h(X), may be represented as:
<IMG>
wherein,
X represents said at least one conditional indicator forming a vector of
individual
conditional indicator values;
Y represents .LAMBDA.-1/2.PHI.T X ;
.LAMBDA. represents (I - K-1)-1, diagonal matrix of eigenvalues of X using a
characteristic
equation of .SIGMA..PHI.=.PHI..LAMBDA., .PHI.T.PHI. = identity matrix;
.PHI. represents an nxn matrix of eigenvectors, .PHI.1.. .PHI.n, of X using
said characteristic
equation;
I represents A T.SIGMA.2A;
K represents A T.SIGMA.2A;
L represents AT(M2 - M1);
P2 represents an aposteriori probability of a first of said two health
classifications
given X;
P1 represents an aposteriori probability of a second of said two health
classifications
given X;
A represents said transformation matrix maximizing a distance between
distributions
of said two health classifications represented as
.LAMBDA.-1/2.PHI.T ;
.SIGMA.1 represents a covariance matrix of said first health classification
having a vector of
expected values M1; and
86

.SIGMA.2 represents a covariance matrix of said second health classification
having a vector
of expected values M2.
19. The computer program product of Claim 12, wherein said at least one
condition
indicator indicates a physical state of a portion of said component.
87

20. A computer program product for determining a health status of a component
comprising machine executable code for:
selecting a plurality of condition indicators having a value and each having a
corresponding weighting factor, and at least one threshold value defining at
least two
classifications;
determining a contribution to a health indicator for each of said condition
indicators,
wherein said machine executable code for determining further comprises, for
each of said
plurality of indicators, machine executable code for:
determining which of said at least two classifications said value of said each
indicator belongs; and
determining said contribution to said health indicator by said each condition
indicator in accordance with a selected one of said at least two
classifications and said
weighted value; and
determining said health indicator in accordance with all contributions by each
of said
condition indicator values.
21. The computer program product of Claim 20, wherein three classifications
are
associated with each of said condition indicators formed by two threshold
values, and the
computer program product further comprises machine executable code for:
determining that said contribution to said health indicator for a first of
said condition
indicators is zero if said value of said condition indicator is in said first
classification;
determining that said contribution to said health indicator is a first
multiple of said
weighting factor if said value of said condition indicator is in said second
classification; and
determining that said contribution to said health indicator is a second
multiple of said
weighting factor if said value of said condition indicator is in said third
classification.
88

22. The computer program product of Claim 21, where the two threshold values
are
alarm level and warning level.
89

23. A method executed in a computer system for determining an health indicator
of a
component at a subsequent time comprising:
determining a first health indicator of said component at a time, n, in
accordance with
at least one corresponding condition indicator; and
using a three state Kalman filter to determine a second health indicator of
said
component at a time subsequent to time n.
24. The method of Claim 23, further comprising:
estimating said first health indicator using a hypothesis determination
technique.
25. The method of Claim 24, wherein the health of a component does not vary by
a
large amount with respect to a change in time.
26. The method of Claim 24, wherein said at least one condition indicator is a
normalized condition indicator.
90

27. The method of Claim 23, wherein
<IMGS>
in which:
.sigma. is a power spectral density,
R is a measurement error,
P is a covariance,
Q is a plant noise,
H is a measurement matrix,
K is a Kalman gain and
.PHI. is state transition matrix, and
X t¦t-1 = .PHI.X t-1¦t-1 ~~(Equation T1)
P t¦t-1 = .PHI.P t-1¦t-1.PHI.T + Q ~(Equation T2)
K = P t¦t-1H T (HP t¦t-1 H T +R) ~~(Equation T3)
P t¦t = (I - KH)P t¦t-1~~ (Equation T4)
X t¦t = X t¦t + K(HI - HX t¦t-1) (Equation T5)
and the method comprising:
using the above Equations T1 through T5 to determine a value Hi_est which is a
best
estimate of a value HI at a time in the future.
91

28. A computer program product for determining an health indicator of a
component
at a subsequent time comprising machine executable code for:
determining a first health indicator of said component at a time, n, in
accordance with
at least one corresponding condition indicator; and
using a three state Kalman filter to determine a second health indicator of
said
component at a time subsequent to time n.
29. The computer program product of Claim 28, further comprising machine
executable code for:
estimating said first health indicator using a hypothesis determination
technique.
30. The computer program product of Claim 29, wherein the health of a
component
does not vary by a large amount with respect to a change in time.
31. The computer program product of Claim 29, wherein said at least one
condition
indicator is a normalized condition indicator.
92

32. The computer program product of Claim 28, wherein
<IMGS>
in which:
.sigma. is a power spectral density,
R is a measurement error,
P is a covariance,
Q is a plant noise,
H is a measurement matrix,
K is a Kalman gain and
.PHI. is state transition matrix, and
X t¦t-1 = .PHI.X t-1¦t-1 ~~(Equation T1)
P t¦t-1 = .PHI.P t-1¦t-1.PHI.T + Q ~(Equation T2)
K = P t¦t-1H T (HP t¦t-1 H T +R) ~~(Equation T3)
P t¦t = (I - KH)P t¦t-1~~ (Equation T4)
X t¦t = X t¦t + K(HI - HX t¦t-1) (Equation T5)
and the computer program product comprises machine executable code for:
determining, using above Equations T1 through T5, to determine a value Hi_est
which is a best estimate of a value HI at a time in the future.
93

PAGE INTENTIONALLY LEFT BLANK
94

33. A method executed in a computer system for ranking condition indicators
used in
determining a health indicator for a component comprising:
determining a first set of a plurality of said condition indicators;
determining a covariance matrix corresponding to said plurality of condition
indicators;
determining a transformation matrix that whitens the covariance matrix;
using said whitening matrix to determine differences between said first
plurality of
condition indicators and expected values for said condition indicators
belonging to a health
class, each health class having a corresponding health indicator; and
selecting a portion of said plurality of condition indicators in accordance
with those
condition indicators have the smallest of said differences.
34. The method of Claim 33, further comprising:
sorting said differences in descending order, each of said differences having
a
corresponding condition indicator.
35. The method of Claim 33, wherein said first plurality of condition
indicators
correspond to an observed data acquisition.
95

36. The method of Claim 33, further comprising:
determining a measure of between class scatter, Sb, represented as:
<IMG>
where M0 is an expected value of all L classes, Mi represents an expected
value of a
particular class, and Pi is the probability of a class i.
37. The method of Claim 36, further comprising:
determining a whitening transformation of Sb as Sbw represented as:
Sbw=A T S b A, for the whitening transformation matrix A.
38. The method of Claim 37, wherein said whitening transformation matrix A is:
.LAMBDA.1/2 .PHI.T with a corresponding eigenvalue matrix .LAMBDA. and a
corresponding eigenvector
matrix .PHI..
96

39. A computer program product for ranking condition indicators used in
determining
a health indicator for a component comprising machine executable code for:
determining a first set of a plurality of said condition indicators;
determining a covariance matrix corresponding to said plurality of condition
indicators;
determining a transformation matrix that whitens the covariance matrix;
using said whitening matrix to determine differences between said first
plurality of
condition indicators and expected values for said condition indicators
belonging to a health
class, each health class having a corresponding health indicator; and
selecting a portion of said plurality of condition indicators in accordance
with those
condition indicators have the smallest of said differences.
40. The computer program product of Claim 39, further comprising machine
executable code for:
sorting said differences in descending order, each of said differences having
a
corresponding condition indicator.
41. The computer program product of Claim 39, wherein said first plurality of
condition indicators correspond to an observed data acquisition.
42. The computer program product of Claim 39, further comprising machine
executable code for:
determining a measure of between class scatter, Sb, represented as:
<IMG>
97

where M0 is an expected value of all L classes, Mi represents an expected
value of a
particular class, and Pi is the probability of a class i.
43. The computer program product of Claim 42, further comprising machine
executable code for:
determining a whitening transformation of Sb as Sbw represented as:
S bw, = A T S b A, for the whitening transformation matrix A.
44. The computer program product of Claim 43, wherein said whitening
transformation matrix A is:
.LAMBDA.-1/2.PHI.T with a corresponding eigenvalue matrix .LAMBDA. and a
corresponding eigenvector
matrix .PHI.
98

45. A method for estimating a conditional indicator value associated with a
gear pair
comprising:
modeling said gear pair as a damped spring model having a contact line between
said
gears;
determining a force, P, at a point of contact along said contact line causing
linear and
torsional response to each of said two gears in said gear pair;
determining a relative movement, d, of said gear pair, in accordance with said
force,
P, as a sum of four responses and a contact deflection, said relative movement
d representing
a gear model having two degrees of freedom; and
using said relative movement, d, in determining said conditional indicator
value for
transmission error associated with said gear pair.
46. The method of Claim 45, further comprising:
varying parameter values in accordance with simulating different fault
conditions; and
determining mean and threshold values for said different fault conditions.
99

47. A method executed in a computer system for estimating a conditional
indicator
value associated with a bearing comprising:
determining a bearing frequency ratio for said bearing;
determining a periodic impulse in accordance with said bearing frequency
ratio;
determining an intensity of an impulse on a bearing surface as a function of
an angle
relative to a bearing fault;
determining a decay of a unit impulse; and
determining a movement of said bearing;
determining a conditional indicator value associated with said bearing in
accordance
with said movement.
48. The method of Claim 47, wherein said bearing frequency ratio corresponds
to one
of: an inner race frequency and an outer bearing race frequency.
49. The method of Claim 47, further comprising:
varying parameter values in determining conditional indicator values
associated with
different health classifications of said bearing.
50. The method of Claim 49, further comprising:
estimating values associated with a gear noise;
estimating values associated with a bearing noise; and
estimating a condition indicator of a gear and bearing noise as a combined
signal
using said values associated with said gear noise and said bearing noise.
51. The method of Claim 50, further comprising:
100

determining said combined signal s(t) as:
s(t)=[d(t)f(t)q(t)a(t)]*e(t)*h(t)
where:
h(t) is a frequency response of a gear case;
d(t) is a signal associated with gear and shaft transmission error;
f(t) is a bearing frequency ratio;
q(t) is an amplitude at a particular time t,
a(t) is a cosine for one of an inner and outer race condition at a location
theta at which
an impulse is applied to said bearing; and
e(t) is a decay rate of unit impulse;
52. The method of Claim 51, wherein h(t) is determined using one of: a linear
predictive coding (LPC) technique and modal hammer analysis.
53. The method of Claim 51, further comprising:
determining a frequency spectrum of signals of s(t) as S(f) in the frequency
domain
by forming:
s(f )=[D(f)*F(f)*Q(f)*A(f)]E(f)H(f)
in the frequency domain.
54. A computer program product for estimating a conditional indicator value
associated with a gear pair comprising machine executable code for:
101

modeling said gear pair as a damped spring model having a contact line between
said
gears;
determining a force, P, at a point of contact along said contact line causing
linear and
torsional response to each of said two gears in said gear pair;
determining a relative movement, d, of said gear pair, in accordance With said
force,
P, as a sum of four responses and a contact deflection, said relative movement
d representing
a gear model having two degrees of freedom; and
using said relative movement, d, in determining said conditional indicator
value for
transmission error associated with said gear pair.
55. The computer program product of Claim 54, further comprising machine
executable code for:
varying parameter values in accordance with simulating different fault
conditions; and
determining mean and threshold values for said different fault conditions.
56. A computer program product for estimating a conditional indicator value
associated with a bearing comprising machine executable code for:
determining a bearing frequency ratio for said bearing;
determining a periodic impulse in accordance with said bearing frequency
ratio;
determining an intensity of an impulse on a bearing surface as a function of
an angle
relative to a bearing fault;
determining a decay of a unit impulse; and
determining a movement of said bearing;
determining a conditional indicator value associated with said bearing in
accordance
with said movement.
102

57. The computer program product of Claim 56, wherein said bearing frequency
ratio
corresponds to one of an inner race frequency and an outer bearing race
frequency.
58. The computer program product of Claim 56, further comprising machine
executable code for:
varying parameter values in determining conditional indicator values
associated with
different health classifications of said bearing.
59. The computer program product of Claim 58, further comprising machine
executable code for:
estimating values associated with a gear noise;
estimating values associated with a bearing noise; and
estimating a condition indicator of a gear and bearing noise as a combined
signal
using said values associated with said gear noise and said bearing noise.
60. The computer program product of Claim 59, further comprising machine
executable code for:
determining said combined signal s(t ) as:
s(t)=[d(t)f(t)q(t)a(t)]8e(t)*h(t)
where:
h(t) is a frequency response of a gear case;
d(t) is a signal associated with gear and shaft transmission error;
f(t) is a bearing frequency ratio;
103

q(t) is an amplitude at a particular time t,
a(t) is a cosine for one of an inner and outer race condition at a location
theta at which
an impulse is applied to said bearing; and
e(t) is a decay rate of unit impulse;
61. The computer program product of Claim 60, wherein h(t) is determined using
one
of: a linear predictive coding technique and modal hammer analysis.
62. The computer program of Claim 60, further comprising machine executable
code
for:
determining a frequency spectrum of signals of s(t) as S(f) in the frequency
domain
represented as
s(f)=[D(f)*F(f)*Q(f)*A(f)]E(f)H(f)
in the frequency domain.
104

63. A method for normalizing a set of at least one observed condition
indicator
comprising:
determining a plurality of conditional indicators and at least one associated
factor in
accordance with previous data acquisitions;
determining a mean and at least one model coefficient corresponding to said at
least
one associated factor; and
adjusting said set of observed condition indicators in accordance with model
coefficients and said at least one associated factor producing a normalized
set of condition
indicators.
64. The method of Claim 63, wherein said relevant factors include at least one
of:
torque, airspeed, and rotor speed.
65. The method of Claim 64, wherein said relevant factors include torque,
airspeed
and rotor speed (Nr), and the method further comprising:
solving for B using the equation form y = B x for the following:
<IMG>
wherein the vector y includes condition indicator values and the matrix x
includes relevant
factors corresponding to one of said condition indicators included in y, and B
includes
coefficients.
105

66. The method of Claim 65, wherein a normalized CI is determined at a later
point in
time by determining:
CInormalized = CIobs - (B * x)
where CIobs is an actual observed value corresponding to a normalized CI value
Clnormalized.
67. The method of Claim 66, further comprising:
storing said values for B and x in a database; and
retrieving said values for B and x in connection with determining
CInormalized.
68. A computer program product for normalizing a set of at least one observed
condition indicator comprising machine executable code for:
determining a plurality of conditional indicators and at least one associated
factor in
accordance with previous data acquisitions;
determining a mean and at least one model coefficient corresponding to said at
least
one associated factor; and
adjusting said set of observed condition indicators in accordance with model
coefficients and said at least one associated factor producing a normalized
set of condition
indicators.
69. The computer program product of Claim 68, wherein said relevant factors
include
at least one of: torque, airspeed, and rotor speed.
106

70. The computer program product of Claim 69, wherein said relevant factors
include
torque, airspeed and rotor speed (Nr), and the computer program product
further comprising
machine executable code for:
solving for B using the equation form y = B x for the following:
<IMG>
wherein the vector y includes condition indicator values and the matrix x
includes relevant
factors corresponding to one of said condition indicators included in y, and B
includes
coefficients.
71. The computer program product of Claim 70, wherein a normalized CI is
determined at a later point in time, and the computer program product further
comprising
machine executable code for determining:
Clnormalized = CIobs - (B * x)
where CIobs is an actual observed value corresponding to a normalized CI value
Clnormalized.
72. The computer program product of Claim 71, further comprising machine
executable code for:
storing said values for B and x in a database; and
retrieving said values for B and x in connection with determining
CInormalized.
107

73. A method executed in a computer system for determining a condition
indicator about a
characteristic of a component comprising:
determining a distribution of observed data associated with said component;
measuring a difference between said distribution and a normal distribution;
and
determining said condition indicator using said difference.
74. The method of Claim 73, further comprising:
determining whether said distribution of observed data is normally distributed
using
said difference using at least normality test that is one of chi-square
goodness of fit test,
Kolmogorov-Smirnof goodness of fit test, Lilliefors test of normality and
Jarque-Bera test of
normality.
75. The method of Claim 74, wherein said normal distribution is one of a
normal
cumulative distribution function and a normal probability distribution
function in accordance
with said at least one normality test.
76. The method of Claim 75, wherein said distribution of observed data
associated
with a component approximates one of a Gaussian distribution if said component
is healthy
and a non-Gaussian distribution otherwise.
108

77. The method of Claim 74, further comprising:
determining a number of differences between said observed data and expected
data,
said expected data being represented by said normal distribution; and
determining a sum using said differences; and
if said number of differences is greater than a critical value, determining
that said
observed data is not normally distributed, said critical value being a
threshold.
78. The method of Claim 77, further comprising:
determining a score being a maximum deviation from said critical value, said
condition indicator being said score.
79. The method of Claim 78, wherein sensitivity of said condition indicator
increases
as a number of observed data values increases.
80. The method of Claim 73, wherein said normal distribution approximates a
distribution of expected values.
109

81. A computer program product for determining a condition indicator about a
characteristic of a component comprising machine executable code for:
determining a distribution of observed data associated with said component;
measuring a difference between said distribution and a normal distribution;
and
determining said condition indicator using said difference.
82. The computer program product of Claim 81, further comprising:
machine executable code for determining whether said distribution of observed
data is
normally distributed using said difference using at least normality test that
is one of: chi-
square goodness of fit test, Kolmogorov-Smirnof goodness of fit test,
Lilliefors test of
normality and Jarque-Bera test of normality.
83. The computer program product of Claim 82, wherein said normal distribution
is
one of a normal cumulative distribution function and a normal probability
distribution
function in accordance with said at least one normality test.
84. The computer program product of Claim 83, wherein said distribution of
observed data associated with a component approximates one of: a Gaussian
distribution if
said component is healthy and a non-Gaussian distribution otherwise.
110

85. The computer program product of Claim 82, further comprising machine
executable code for:
determining a number of differences between said observed data and expected
data,
said expected data being represented by said normal distribution; and
determining a sum using said differences; and
determining that said observed data is not normally distributed, said critical
value
being a threshold if said number of differences is greater than a critical
value.
86. The computer program product of Claim 85, further comprising:
machine executable code for determining a score being a maximum deviation from
said critical value, said condition indicator being said score.
87. The computer program product of Claim 86, wherein sensitivity of said
condition
indicator increases as a number of observed data values increases.
88. The computer program product of Claim 81, Wherein said normal distribution
approximates a distribution of expected values.
111

89. A method executed in a computer system for determining a condition
indicator
associated with a component, the method comprising:
determining a total impulse signal in accordance with configuration data, said
total
impulse signal being a superposition of gear and bearing noise represented as
a convolution
of a gear and bearing signal with a gearbox transfer function; and
determining a condition indicator in accordance with said total impulse
signal.
90. The method of Claim 89, further comprising:
representing a total impulse signal generated by a configuration of associated
with
said component as:
[impulse]~ f(Gear)~ f(Bearing)~ f(Case).ident.[impulse]~
[f(Gear)~f(Bearing)~f(Case)]
in which ~ represents a convolution operation.
91. The method of Claim 90, further comprising:
representing convolution operations in a time domain to equivalent operations
in a
frequency domain.
92. The method of Claim 90, further comprising:
estimating [f(Gear)~ f(Bearing)~ f(Case)] as a transfer function in a
frequency
domain using a linear predictive coding technique to deconvolute a signal into
its base
components.
93. The method of Claim 92, further comprising:
112

.
estimating said transfer function, H, in said frequency domain as a/B, wherein
a =
(al, .., an), each ai representing an ith coefficient for an order p, n=p+1,
as:
y[n] = .alpha.1x[n -1]+.alpha.2x[n -2]+ .alpha.3x[n - 3]+...
and B is an estimate of an error represented as:
B = .SIGMA.all b
in which:
b = (Y-Yhat)2, Y = Y[1,2,...n],
Y hat is an estimated value of y, yhat = ax,
x is a time delayed signal represented as:
<IMG>
where a = (x T x)-1 x T y, values for al .. an.
94. The method of Claim 93, further comprising:
estimating an impulse, IMP, in said frequency domain of said component as:
IMP = exp(log(Y) - log(H)),
in which:
Y = fft(y) and H = fft(h), where fft is the Fourier transform function, y and
h are in a
time domain, Y and H are in said frequency domain.
95. The method of Claim 94, wherein a value associated with H increases as a
fault
increases.
113

96. The method of Claim 94, wherein said condition indicator is said value of
IMP.
97. The method of Claim 94, further comprising:
calculating a power spectral density of said impulse IMP in a frequency
domain; and
determining a value of said power spectral density at a frequency of interest,
said
condition indicator being said value.
98. The method of Claim 97, wherein said frequency of interest is at least one
of: a
bearing passing frequency for a bearing fault, and a mesh frequency for a gear
fault.
99. The method of Claim 98, further comprising:
performing a Fourier transformation to obtain IMP in said frequency domain.
100. The method of Claim 89, further comprising:
detecting a fault in connection with predetermined values of said health
status using
said condition indicator, wherein said fault being detected is one of a pit
and spall on one of:
a gear tooth, inner bearing race, outer bearing race, and bearing roller
element.
114

101. A computer program product for determining a condition indicator
associated
with a component, the computer program product comprising machine executable
code for:
determining a total impulse signal in accordance with configuration data, said
total
impulse signal being a superposition of gear and bearing noise represented as
a convolution
of a gear and bearing signal with a gearbox transfer function; and
determining a condition indicator in accordance with said total impulse
signal.
102. The computer program product of Claim 101, further comprising machine
executable code for:
representing a total impulse signal generated by a configuration of associated
with
said component as:
[impulse]~ f(Gear)~ f(Bearing)~ f(Case) .ident. [impulse]~ f(Gear)~
f(Bearing)~ f(Case)]
in which ~ represents a convolution operation.
103. The computer program of Claim 102, further comprising machine executable
code for:
representing convolution operations in a time domain to equivalent operations
in a
frequency domain.
115

104. The computer program product of Claim 102, further comprising machine
executable code for:
estimating [.function.(Gear)~.function.(Bearing)~.function.(Case)] as a
transfer function in a frequency
domain using a linear predictive coding technique to deconvolute a signal into
its base
components.
105. The computer program product of Claim 104, further comprising machine
executable code for:
estimating said transfer function, H, in said frequency domain as a/B, wherein
a =
(a1, .., an), each ai representing an ith coefficient for an order p, n=p+1,
as:
y[n]=a1x[n-1]+a2x[n-2]+a3x[n-3]+...
and B is an estimate of an error represented as:
B = .SIGMA.all b
in which:
b = (Y-Y hat)2, y = y[1, 2, ... n],
y hat is an estimated value of y, y hat = ax,
x is a time delayed signal represented as:
<IMG>
where a = (x T x)-1 x T y, values for a1 .. an.
106. The computer program product of Claim 105, further comprising machine
executable code for:
116

estimating an impulse, IMP, in said frequency domain of said component as:
IMP = exp(log(Y) - log(H)),
in which:
Y = fft(y) and H = fft(h), where fft is the Fourier transform function, y and
h are in a
time domain, Y and H are in said frequency domain.
107. The computer program product of Claim 106, wherein a value associated
with H
increases as a fault increases.
108. The computer program product of Claim 106, wherein said condition
indicator
is said value of IMP.
109. The computer program product of Claim 106, further comprising machine
executable code for:
calculating a power spectral density of said impulse IMP in a frequency
domain; and
determining a value of said power spectral density at a frequency of interest,
said
condition indicator being said value.
110. The computer program product of Claim 109, wherein said frequency of
interest
is at least one of: a bearing passing frequency for a bearing fault, and a
mesh frequency for a
gear fault.
111. The computer program product of Claim 110, further comprising machine
executable code for
performing a Fourier transformation to obtain IMP in said frequency domain.
117

112. The computer program product of Claim 111, further comprising machine
executable code for:
detecting a fault in connection with predetermined values of said health
status using
said condition indicator, wherein said fault being detected is one of a pit
and spall on one of:
a gear tooth, inner bearing race, outer bearing race, and bearing roller
element.
118

113. A method executed in a computer system for determining a health status of
a
component using at least one condition indicator, the method comprising:
determining said at least one condition indicator using at least one of an
impulse
determination technique and a statistical normality test; and
determining said health indicator in accordance with said at least one
condition
indicator.
114. The method of Claim 113, wherein said statistical normality test is one
of:chi-
square goodness of fit test, Kolmogorov-Smirnof goodness of fit test,
Lilliefors test of
normality and Jarque-Bera test of normality.
115. The method of Claim 113, wherein expected data values approximate a
normal
distribution.
119

116. A computer program product for determining a health status of a component
using at least one condition indicator, the method comprising:
determining said at least one condition indicator using at least one of: an
impulse
determination technique and a statistical normality test; and
determining said health indicator in accordance with said at least one
condition
indicator.
117. The computer program product of Claim 116, wherein said statistical
normality
test is one of chi-square goodness of fit test, Kolmogorov-Smirnof goodness of
fit test,
Lilliefors test of normality and Jarque-Bera test of normality.
118. The computer program product of Claim 116, wherein expected data values
approximate a normal distribution.
120

Description

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
METHOD AND APPARATUS FOR DETERMINING THE HEALTH OF A
COMPONENT USING CONDITION INDICATORS
Background of The Invention
1. Field of the Invention
This application relates to the field of vibration analysis and more
particularly to
performing vibration analysis for the purpose of device monitoring.
l0 2. Descriution of Related Art
The transmission of power to rotors which propel helicopters and other shafts
that
propel devices within the aircraft induce vibrations in the supporting
structure. The
vibrations occur at frequencies that correspond to the shaft rotation rate,
mesh rate, bearing
passing frequency, and harmonics thereof. The vibration is associated with
transmission
error (TE). Increased levels of TE are associated with transmission failure.
Similar types of
vibrations are produced by transmissions in fixed installations as well.
Parts, such as those that may be included in a helicopter transmission, may be
replaced in accordance with a predetermined maintenance and parts replacement
schedule.
These schedules provide for replacement of parts prior to failure. The
replacement schedules
may indicate replacement time intervals that are too aggressive resulting in
needless
replacement of working parts. This may result in incurnng unnecessary costs as
airplane
parts are expensive. Additionally, new equipment may have installed faulty or
defective
parts that may fail prematurely.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Thus it may be desirable to provide for an efficient technique for detecting
part and
device degradation without unnecessarily replacing parts. It may be desirable
that this
technique also provide for problem determination and detection prior to
failure.
2

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Summary Of The Invention
In accordance with one aspect of the invention are a method executed in a
computer
system and a computer program product for determining a health indicator
associated with a
component. A plurality of health classifications are determined. At least one
condition
indicator is determined quantifying a characteristic of the component. A
probability
associated with each of the health classifications is determined. The
probability is an
estimation that the component is of a particular health classification given
the at least one
indicator. A determination is made as to which of said health classifications
is associated
1o with said component using said probabilities associated with said health
classifications for a
given set of observed values.
In accordance with another aspect of the invention are a method executed in a
computer system and a computer program product for determining a health status
of a
15 component. A plurality of condition indicators are selected having a value
and each having a
corresponding weighting factor, and at least one threshold value defining at
least two
classifications. A contribution to a health indicator is determined for each
of said condition
indicators, wherein said determining further comprises, for each of said
plurality of
indicators: determining which of said at least two classifications said value
of said each
2o indicator belongs; and determining said contribution to said health
indicator by said each
condition indicator in accordance with a selected one of said at least two
classifications and
said weighted value. The health indicator is determined in accordance with all
contributions
by each of said condition indicator values.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
In accordance with one aspect of the invention are a method executed in a
computer
system and a computer program product for determining an health indicator of a
component
at a subsequent time. A first health indicator of said component at a time, n,
is determined in
accordance with at least one corresponding condition indicator. A second
health indicator of
the component is determined using a three state Kalman filter at a time
subsequent to time n.
In accordance with one aspect of the invention are a method executed in a
computer
system and a computer program product for ranking condition indicators used in
determining
1o a health indicator for a component. A first set of a plurality of said
condition indicators is
determined. A covariance matrix corresponding to said plurality of condition
indicators is
determined. A transformation matrix that whitens the covariance matrix is
determined.
Differences between said first plurality of condition indicators and expected
values for said
condition indicators belonging to a health class are determined using the
whitening matrix.
15 Each health class has a corresponding health indicator. A portion of said
plurality of
condition indicators is selected in accordance with those condition indicators
have the
smallest of said differences.
In accordance with one aspect of the invention are a method executed in a
computer
2o system and computer program product for estimating a conditional indicator
value associated
with a gear pair. The gear pair is modeled as a damped spring model having a
contact line
between said gears. A force, P, is determined at a point of contact along said
contact line
causing linear and torsional response to each of said two gears in said gear
pair. A relative
movement, d, is determined of said gear pair, in accordance with said force,
P, as a sum of
25 four responses and a contact deflection, said relative movement d
representing a gear model

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
having two degrees of freedom. The relative movement, d, is used in
determining said
conditional indicator value for transmission error associated with said gear
pair.
In accordance with another aspect of the invention are a method executed in a
computer system and a computer program product for estimating a condition
indicator
associated with a bearing. A bearing frequency ratio is determined for the
bearing. A
periodic impulse is determined in accordance with the bearing frequency ratio.
An intensity
of an impulse on a bearing surface as a function of an angle relative to a
bearing fault is
determined. A decay of a unit impulse is determined. A movement of the bearing
is
to determined. A conditional indicator value is determined in accordance with
the movement.
In accordance with one aspect of the invention are a method executed in a
computer
system and computer program product 'for normalizing a set of at least one
observed
condition indicator. A plurality of conditional indicators and at least one
associated factor are
15 determined in accordance with previous data acquisitions. A mean and at
least one model
coefficient corresponding to said at least one associated factor are
determined. The set of
observed condition indicators is adjusted in accordance with model
coefficients and said at
least one associated factor producing a normalized set of condition
indicators.
20 In accordance with one aspect of the invention are a method executed in a
computer
system and computer program product for determining a condition indicator
about a
characteristic of a component. A distribution of observed data associated with
said
component is determined. A difference between said distribution and a normal
distribution is
determined. The condition indicator is determined using the difference.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
In accordance with another aspect of the invention are a method executed in a
computer system and computer program product for determining a condition
indicator
associated with a component. A total impulse signal is determined in
accordance with
configuration data. The total impulse signal is a superposition of gear and
bearing noise
represented as a convolution of a gear and bearing signal with a gearbox
transfer function. A
condition indicator is determined in accordance with the total impulse signal.
In accordance with yet another aspect of the invention is a method executed in
a
computer system and a computer program product for determining a health status
of a
to component using at least one condition indicator. At least one condition
indicator is
determined using at least one of an impulse determination technique and a
statistical
normality test. The health indicator is determined in accordance with the at
least one
indicator.
6

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Brief Description Of Drawings
Features and advantages of the present invention will become more apparent
from the
following detailed description of exemplary embodiments thereof taken in
conjunction with
the accompanying drawings in which:
Figure 1 is an example'of an embodiment of a system that may be used in
performing
vibration analysis and performing associated monitoring functions;
Figure 2 is an example representation of a data structure that includes
aircraft
1o mechanical data;
Figure 3 is an example of parameters that may be included in the type-specific
data
portions when the descriptor type is an indexer;
15 Figure 4 is an example of parameters that may be included in the type-
specific data
portions when the descriptor type is an accelerometer;
Figure 5 is an example of parameters that may be included in the type-specific
data
portions when the descriptor type is a shaft;
Figure 6 is an example of parameters that may be included in the type-specific
data
portions when the descriptor type is for a gear;
Figure 7 is an example of parameters that may be included in the type-specific
data
portions when the descriptor type is a planetary type;

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Figure 8 is an example of parameters that may be included in the type-specific
data
portions when the descriptor type is bearing type;
Figure 9 is an example of a data structure that includes analysis information;
Figure 10 is a more detailed example of an embodiment of a header descriptor
of
Figure 9;
1 o Figure 11 is an example of a descriptor that may be included in the
acquisition
descriptor group of Figure 9;
Figure 12 is an example of a descriptor that may be included in the
accelerometer
group of Figure 9;
Figure 13 is an example of a descriptor that may be included in the shaft
descriptor
group of Figure 9;
Figure 14 is an example of a descriptor that may be included in the signal
average
2o descriptor group of Figure 9;
Figure 15 is an example of a descriptor that may be included in the envelope
descriptor group of Figure 9;
Figure 16 is an example of a planetary gear arrangement;
8

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Figure 17A is an example of an embodiment of a bearing;
Figure 17B is an example of a cut along a line of Figure 17A;
Figure 18A is an example of a representation of data flow in vector
transformations;
Figure 18B is an example of a representation of some of the CI algorithms that
may
be included in an embodiment, and some of the various inputs and outputs of
each;
Figure 19 is an example of a graphical representation of a probability
distribution
function (PDF) of observed data;
Figure 20 is an example of a graphical representation of a cumulative
distribution
function (CDF) observed data following a gamma (5,20) distribution and the
normal CDF;
Figure 21 is an example of a graphical representation of the difference
between the
two CDFs of Figure 20;
2o Figure 22 is an example of a graphical representation of the PDF of
observed data
following a Gamma (5,20) distribution and a PDF of the normal distribution;
Figure 23 is an example of another graphical representation of the two PDFs
from
Figure 22 shown which quantities as intervals rather than continuous lines;

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Figure 24A is an example of a graphical representation of the differences
between the
two PDFs of observed data and the normally distributed PDF;
Figures 24B-24D are examples of a graphical data displays in connection with a
healthy system;
Figures 24E-24G are examples of graphical data displays in connection with a
system
having a fault;
Figure 25 is a flowchart of steps of one embodiment for determining health
indicators
(HIs);
Figure 26 is a graphical illustration of the probability of a false alarm
(PFA) in one
example;
Figure 27 is a graphical illustration of the probability of detection (PD) in
one
example;
Figure 2g is a graphical illustration of the relationship between PD and PFA
and
2o threshold values in one embodiment;
Figure 29 is an graphical illustration of the probability of Ho and threshold
values in
one embodiment;
Figure 30 is an example of an embodiment of a gear model;
to

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Figure 31 is a graphical representation of an estimated signal having an inner
bearing
fault; and
Figure 32 is a graphical representation of the signal of Figure 31 as a
frequency
spectrum.
11

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Detailed Description of the Preferred Embodiments)
Referring now to Figure l, shown is an example of an embodiment of a system 10
that may be used in performing vibration analysis and monitoring of a machine
such as a
portion of an aircraft. The machine being monitored 12 may be a particular
element within an
aircraft. Sensors 14a through 14c are located on the maclune to gather data
from one or more
components of the machine. Data may be collected by the sensors 14a through
14c and sent
to a processor or a VPU16 for data gathering and analysis. The VPU16 analyzes
and gathers
the data from the Sensors 14a through 14c.
The VPU16 may also use other data in performing analysis. For example, the
VPU16
may use collected data 18. One or more of the Algorithms 20 may be used as
input into the
VPU16 in connection with analyzing data such as may be gathered from the
Sensors 14a
through 14c. Additionally, configuration data 22 may be used by the VPUl6 in
connection
with performing an analysis of the data received for example from the Sensors
14a through
14c. Generally, configuration data may include parameters and the like that
may be stored in
a configuration data file. Each of these will be described in more detail in
paragraphs that
follow.
2o The VPU16 may use as input the collected data 18, one or more of the
algorithms 20,
and configuration data 22 to determine one or more condition indicators or
CIs. In turn, these
condition indicators may be used in determining health indicators or HIs that
may be stored
for example in CI and HI storage 28. CIs describe aspects about a particular
component that
may be useful in making a determination about the state or health of a
component as may be
reflected in an HI depending on one or more CIs. Generally, as will be
described in more
12

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
detail in paragraphs that follow, CIs and HIs may be used in connection with
different
techniques in determining an indication about monitored components such as
Machine 12.
As described in more detail elsewhere herein, the configuration data may
include values for
parameters that may vary in accordance with the type of the component being
monitored.
It hould be noted that the collected data 18 may include data collected over a
period
of time from sensors such as 14a through 14c mounted on Machine 12. A user,
such as a Pilot
26, may use a special service processor, such as the PPU24, connected to the
Machine 12 to
obtain different types of data such as the CI and HI values 28.
to
As described in connection with Figure 1, the VPU16 may receive inputs from
Sensors 14a through 14c. These sensors may be different types of data
gathering monitoring
equipment including, for example, high resolution accelerometers and index
sensors
(indexors) or tachometers that may be mounted on a component of Machine 12 at
carefully
15 selected locations throughout an aircraft. Data from these sensors may be
sampled at high
rates, for example, up to 100 kilohertz, in order for the VPU16 to produce the
necessary CI
and HI indicators. Data from these sensors and accelerometers may be acquired
synchronously at precise intervals in measuring vibration and rotational
speeds.
2o Generally, the different types of data gathering equipment such as 14a-14c
may be
sensors or tachometers and accelerometers. Accelerometers may provide
instantaneous
acceleration data along whatever axis on which they are mounted of a
particular device.
Accelerometers may be used in gathering vibration analysis data and
accordingly may be
positioned to optimally monitor vibration generated by one or more mechanical
components
25 such as gears, shafts, bearings or planetary systems. Each component being
monitored may
13

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
generally be monitored using two independent sensors to provide confirmation
of component
faults and to enable detection of sensor faults.
No accelerometer is completely isolated from any other component. Thus, the
component rotational frequencies share as few common divisors as possible in
order to
maximize the effectiveness of the monitoring function being performed. For
example, all
gears being monitored should have differing number of teeth and all bearings
should have
differing numbers and sizes of balls or rollers. This may allow individual
components to be
spectrally isolated from each other to the extent that their rotational
frequencies are unique.
The indexers (index sensors) or tachometers may also be used as a particular
monitoring component 14a through 14c to gather data about a particular
component of
Machine 12. The indexers produce a periodic analog signal whose frequency is
an integer
multiple of the instantaneous rotation frequency of the shaft that they are
monitoring. These
signals may be generated magnetically using one or more evenly spaced metallic
protrusions
on the shaft passing by the fixed sensor. Alternatively, these may be
monitored optically
using a piece of optically reflective material affixed to the shaft. It should
be noted that each
index point should be fixed in time as precisely as possible. In connection
with magnetic
sensors, this may be accomplished for example by interpolating the zero
crossing times of
each index pulse and similarly for optical sensors by locating either rising
or falling edges.
Assuming the minimal play or strain in the drive train when something is under
load, the
relative position and rate of any component may be calculated using a single
index or wave
form.
14

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Because of the high data rates and lengthy processing intervals, diagnostics
may be
performed, for example, on pilot command or on a predetermined flight regime
or time
interval.
Each of the algorithms 20 produces one or more CIs described elsewhere herein
in
more detail. Generally, the CI may yield useful information about the health
of a monitored
component. This condition indicator or CI as well as HI may be used in
determining or
predicting faults of different components.
to It should be noted that the VPUl6 is intended to be used in a wide variety
of
mechanical and electrical environments. As described herein, different
components of an
aircraft may be monitored. However, this is only one example of a type of
environment in
which the system described herein may be used. As known to those skilled in
the art, the
general principles and techniques described herein have much broader and
general
1 5 applicability beyond a specific aircraft environment that may used in an
example here.
In connection with the use of CIs, the VPU16 uses the CIs as input and
portions of the
data such as, for example, used in connection with an algorithm to provide
HIs. These are
described in more detail in paragraphs that follow.
It should be noted that in a particular embodiment, each mechanical part being
monitored may have one or more sensors associated with it where a sensor may
include for
example an accelerometer or a tachometer. Generally, accelerometers may be
used, for
example, to obtain data regarding vibrations and a tachometer may be used, for
example, to

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
gain information and data regarding rotation or speed of a particular object.
Data may be
obtained and converted from the time to the frequency domain.
A particular algorithm may provide one or more CIs. Each of the algorithms may
produce or be associated with a particular CI. One or more CIs may be used in
combination
with a function to produce an HI for a particular part or type. As will be
described in more
detail herein, each of the algorithms may be associated or classified with a
particular part or
type. The CI generally measures vibrations and applies a function as described
in accordance
for each algorithm. Generally, vibration is a function of the rotational
frequency in the
to amount of torque. Using torque and a particular frequency, a CI is
appropriately determined
in accordance with a selected algoritlnn for a part.
The algorithms 20 may be classified into four families or groups in accordance
with
the different types of parts. In this example, the families of algorithms may
include shaft,
15 gears, bearings, and planetary gears. Associated with each particular part
being monitored
may be a number of CIs. Each CI may be the result or output of applying a
different one of
the algorithms for a particular family. For example, in one embodiment, each
gear may have
an associated 27 CIs, each bearing may have 19 CIs, each shaft may have 22
CIs, and each
planetary gear may have two or three CIs. It should be noted that each one of
these numbers
2o represents in this example a maximum number of CIs that may be used or
associated with a
particular type in accordance with the munber of algorithms associated with a
particular class
or family. Generally, the different number of CIs that may be associated with
a particular
type such as a gear try to take into account the many different ways in which
a particular gear
may fail. Thus, a CI reflects a particular aspect or characteristic about a
gear with regard to
25 how it may fail.
16

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Different techniques used in computing CIs are described, for example, in
"Introduction to
Machinery Analysis and Monitoring, Second Edition", 1993, Penn Well Publishing
Company
of Tulsa, OK, ISBN 0-87814-401-3, and "Machinery Vibration: measurement and
analysis",
1991, McGraw-Hill Publishing, ISBN-0-07-071936-5.
Referring now to Figure 2, shown is an example of a data structure 50 that
includes
aircraft mechanical data. Generally, this data structure includes one or more
descriptors 56a
through 56n. In this embodiment there may be one descriptor for each sensor. A
descriptor
associated with a particular sensor includes the parameters relevant to the
particular
l0 component being monitored. Each of the descriptors such as 56a includes
three portions of
data. The field 52 identifies a particular type of descriptor. Each of the
descriptors also
includes a common data portion 54 which includes those data fields common to
all descriptor
types. Also included is a type specific data portion 56 which includes
different data fields,
for example, that may vary in accordance with the descriptor type 52.
Descriptor types may include, for example, an indexer, an accelerometer, a
shaft, a
gear, a planetary gear, or a bearing descriptor type value corresponding to
each of the
different types of descriptors. The common data portion 54 may include, for
example, a
name, part number and identifier. In this example, the identifier in the
common data filed 54
may uniquely identify the component and type.
Referring now to Figures 3 through 8, what will be described are examples of
descriptor type specific parameters or information that may be included in a
descriptor of a
particular type, such as in area 56 o~f the data structure 50.
17

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Referring now to Figure 3, shown is an example of parameters that may be
included
in a descriptor 60 which is an indexer descriptor type. The parameters that
may be included
are a channel 62, a type 64, a shaft identifier 66, a pulses per revolution
parameter 68, a pulse
width parameter 70, and a frequency of interest 72 for this particular type of
descriptor. It
should be noted that the type in this example for the index ox descriptor may
be one of
sinusoidal, pulse such as 1/ rev, or optical. The shaft identifier 66 is that
as may be read or
viewed by the indexer that calculates the shaft rate. The pulse width 70 is in
seconds as the
unit value. Additionally, the frequency of interest 72 for this descriptor
type is a nominal
pulse frequency that is used in computing the data quality signal to noise
ratio. The use of
to these particular data structures and parameters is described in more detail
in paragraphs that
follow.
Referring now to Figure 4, shown is an example of the parameters that may be
included in an accelerometer descriptor type 80. The descriptor for an
accelerometer type
15 may include the channel 82, a type 84, a sensitivity 86 and a frequency of
interest 88. In this
example for the accelerometer descriptor type, the type may be one of normal,
or remote
charge coupled. The frequency of interest may be used in computing the data
quality signal
to noise ratio. The frequency of interest for a gear is the mesh rate which
may be calculated
from the gear shaft rate and the number of teeth of the gear.
Referring now to Figure 5, shown is an example of descriptor type specific
parameters
or data that may be included when a descriptor type is the shaft descriptor. A
shaft descriptor
90 includes path parameter or data 92 and nominal RPM data 94. The path data
is an even
length sequence of gear tooth counts in the mechanical path between the shaft
in question and
a reference shaft. The driving gears alternate with driven gears such that the
expected
18

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
frequency of a gear, shaft, bearing and the like may be determined based on an
input shaft
RPM.
Referring now to Figure 6, shown is an example of data or parameters that may
be
included in a descriptor when the descriptor type is the gear descriptor.
Included in the gear
descriptor 100 is the shaft identifier 102 to which the gear is mounted and a
parameter 104
indicating the number of teeth in the gear.
Referring now to Figure 7, shown is an example of an embodiment of a planetary
l0 descriptor 110 identifying those parameters or data that may be included
when the type is a
planetary descriptor type. The planetary descriptor 110 may include an input
shaft identifier
112, an output shaft identifier 114, a parameter indicating the number of
planet gears 116, a
parameter indicating the number of teeth on the planet gear, a parameter 120
indicating the
number of teeth on the ring gear, and a parameter 122 indicating the number of
teeth on the
15 sun gear. It should be noted that the number of teeth on a planet gear
relates to a planet
carrier that is assumed to be mounted to the output shaft. Additionally, the
ring gear is
described by parameter 120 is assumed to be stationery and the sun gear 122 as
related to
parameter 122 is assumed to be mounted to the input shaft. It should be noted
that the path
between the input and the output shaft may be reduced to using a value S for
the driving path
20 tooth count and R+S as the driven path tooth count where R and S are the
ring and sun tooth
counts respectively. An example of a planetary type geax is described in more
detail
elsewhere herein.
Referring now to Figure 8, shown is an example of a bearing descriptor 130.
The
25 bearing descriptor 130 may include descriptor type specific fields
including a shaft identifier
19

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
132, a cage ratio 134, a ball spin ratio 136, an outer race ratio 138 and an
inner race ratio 140.
An example of a bearing is described in more detail elsewhere herein.
It should be noted that the data structures described in connection with
Figures 2
through 8 are those that may be used in storing data obtained and gathered by
a sensor such
as 14a when monitoring a particular component of a machine 12. Data may be
gathered and
stored in the data structure for a particular descriptor or descriptors and
sent to the VPU 16
for processing. It should be noted that a particular set of data may be
gathered at a particular
instance and time, for example, in connection with the synchronous data
gathering described
1o elsewhere herein. In connection with this, a data set may include multiple
descriptors from
sampling data at a particular point in time which is sent to the VPU 16.
What will now be described are those data structures that may be associated
with an
analysis definition that consists of a specific data acquisition and a
subsequent processing of
15 this data to produce a set of indicators for each of the desired
components.
Referring now to Figure 9, shown is an example of the data structure 150 that
contains
analysis data. Each instance of analysis data 150 as represented in the data
structure includes
a header descriptor 152 and descriptor groups noted as 164. In this example
there are five
2o descriptor groups although the particular number may vary in an embodiment.
Each of the
descriptor groups 154 through 162 as identified by the group identifier 164
includes one or
more descriptors associated with a particular group type. For example,
descriptor group 154
is the acquisition group that includes a descriptor for each sensor to be
acquired. The
accelerometer group 156 consists of a descriptor for each accelerometer to be
processed. The
25 shaft group 158 includes a descriptor for each shaft to be processed. The
signal average

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
group 160 includes a descriptor for each unique parameter set. The envelope
group 162
includes a descriptor for each unique parameter.
Referring now to Figure 10, shown is a more detailed example of a header
descriptor
170. Parameters that may be included in a header descriptor 170 include: an
analysis
identifier 172, acquisition time out parameter 174 and processing time out
parameter 176. In
this example, the acquisition, time out and processing time out parameters are
in seconds.
Referring now to Figure 11, shown is an example of a descriptor that may be
included
to in the acquisition group. A descriptor 180 included in the acquisition
group may include a
sensor identifier 182, a sample rate parameter in Hz 184, a sample duration in
seconds 186, a
gain control setting, such as "auto" or "fixed" 188, an automatic gain control
(AGC)
acquisition time in seconds 190, an automatic gain control (AGC) headroom
factor as a
number of bits 192 and a DC offset compensation enable 194.
Referring now to Figure 12, shown is an example of a descriptor 200 that may
be
included in the accelerometer group. A descriptor in the accelerometer group
may include a
parameter that is an accelerometer acquisition analysis group identifier 202,
a list of
associated planetary identifiers to be processed 204, a list of associated
shaft analysis group
2o identifiers to be processed 206, a processor identifier 208, a transient
detection block size
210, a transient detection RMS factor 212, a power spectrum decimation factor
214 specified
as a power of 2 and a power spectrum block size also specified as a power of
2.
21

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
In one embodiment, the list of associated planetary identifiers 204 also
includes two
signal average analysis group identifiers for each planetary identifier, first
identifier
corresponding to the input shaft and a second corresponding to an output
shaft.
It should be noted that the processor identifier 208 will be used in
connection with
assigning processing to a particular DSP or digital signal processor.
Referring now to Figure 13, shown is an example of an embodiment of a
descriptor
280 that may be included in the shaft group. The descriptor 220 may include a
shaft
1o identifier 222, a signal average analysis group identifier 224, a list of
gear identifiers to be
processed 226, a list of bearing identifiers to be processed 228 and a list of
associated
envelope analysis group identifiers 230.
Referring now to Figure 14, shown is an example of a descriptor 232 that may
be
15 included in the signal average group. It should be noted that the signal
average group
includes a descriptor for each unique parameter set. The signal average
processing group is
run for each accelerometer and shaft combination even if it has the same
parameters as
another combination. Each descriptor 232 may include a number of output points
per
revolution 234 and a number of revolutions to average 236.
Referring now to Figure 15, shown is an example of a descriptor 240 that may
included in the envelope group. It should be noted that the envelope group
includes a
descriptor for each unique parameter. It is not necessary to repeat an
envelope processing for
each bearing if the parameters are the same. Each descriptor 240 may include a
duration
parameter 242 specifying the seconds of raw data to process, an FFT size 244
which is a
22

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
power of 2, a lower bound frequency in Hz 246, and an upper bound frequency,
also in, Hz
248.
Refernng now to Figure 16, shown is an example of an embodiment 300 of a
planetary gear arrangement. Generally, a planetary gear arrangement as
described in
connection with the different types of gears and items to be monitored by the
system 10 of
Figure 1 may include a plurality of gears as configured, for example, in the
embodiment 300.
Included in the arrangement 300 is a ring gear 302 a plurality of planet gears
304a through
304c and of sun gear 306. Generally, the gears that are designated as planets
move around
to the sun gear similar to that as a solar system, hence the name of planet
gear versus sun gear.
The arrangement shown in Figure 16 is a downward view representing the
different types of
gears included in an arrangement 300.
Refernng now to Figure 17A, shown is an example of an embodiment 320 of a
15 bearing. The bearing 320 includes a ring or track having one or more
spherical or cylindrical
elements (rolling elements) 324 moving in the direction of circular rotation
as indicated by
the arrows. Different characteristics about such a structure of a bearing may
be important as
described in connection with this embodiment. One characteristic is an "inner
race" which
represents the circumference of circle 322a of the inner portion of the ring.
Similarly, the
20 "outer race" or circumference 322b representing the outer portion of the
ring may be a
consideration in connection with a bearing.
Refernng now to Figure 17B, shown is an example of a cut along line 17B of
Figure
17A. Generally, this is cut through the ring or track within which a bearing
or bearings 324
23

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
rotate in a circular direction. The ball bearings move in unison with respect
to the shaft
within a cage that follows a track as well as rotate around each of their own
axis.
Referring now to Figure 18A, shown is an example of a representation 550 of
different transformations that may be performed and the associated data flow
and
dependencies for each particular sensor. The output of the transformations are
transformation
vectors and may be used in addition to analysis data or raw data, such as
bearing frequency,
mesh frequency, and the like, by an algorithm in producing a CI.
Referring to the representation 550, an in going arrow represents data flow
input to a
to transformation. For example, the FF or Fast Fourier transform takes as an
input data from the
A1 signal average data transform. A1 has as input the accelerometer data AD.
It should be
noted that other embodiments may produce different vectors and organize data
inputsloutputs
and intermediate calculations in a variety of different ways as known to those
skilled in the
art.
Referring now to Figure 18B, shown is an example of a representation 350
relating
algoritluns, a portion of input data, such as some transformation vectors, and
CIs produced
for each type of component, that may be included in an embodiment. Other
embodiments
may use different data entities in addition to those shown in connection with
Figure 18B. As
2o described elsewhere herein, each type of component in this example is one
of indexer,
accelerometer, shaft, gear, planetary, or bearing. Certain algorithms may be
used in
connection with determining one or more CIs for more than one component type.
It should
be noted that a variety of different algorithms may be used and are known by
one of ordinary
skill in the ant, as described elsewhere herein in more detail. The following
are examples of
some of the different techniques that may be used in producing CIs.
Additionally, Figure
24

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
18B illustrates an example of relationships between some algorithms, a portion
of their
respective inputs and outputs, as well as how the algorithms may be associated
with different
component types. However, it should be noted that this illustration is not all
inclusive of all
algorithms, all respective inputs and outputs, and all component types.
What will now be described are algorithms and the one or more CIs produced
that
may be included in an embodiment. It should be noted that the number and type
of
algorithms included may vary in accordance with an embodiment. Additionally,
it should be
noted that Figure 18B may not include each and every input and output for an
algorithm as
1o described herein and other embodiments of the algorithms described
generally herein may
also vary.
The data quality (DQ) algorithm 356 may be used as a quality assurance tool
for the
DTD CI. DQ performs an assessment of the raw uncalibrated sensor data to
insure that the
entire system is performing nominally. DQ may be used to identify, for
example, bad wiring
connections, faulty sensors, clipping, and other typical data acquisition
problems. The DQ
indicator checks the output of an accelerometer for "bad data". Such "bad
data" causes the SI
to be also be "bad" and should not be used in determining health calculations.
2o What will now be described are the different indicators that may be
included in an
embodiment of the DQ algorithm. ADC Bit Use measures the number of ADC bits
used in
the current acquisition. The ADC board is typically a 16 bit processor. The
log base 2 value
of the maximum raw data bit acquired is rounded up to the next highest
integer. Channels
with inadequate dynamic range typically use less than 6 bits to represent the
entire dynamic
range. ADC Sensor- Range is the maximum range of the raw acquired data. This
range

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
cannot exceed the operational range of the ADC board, and the threshold value
of 32500 is
just below the maximum permissible value of +32767 or -32768 when the absolute
value is
taken. Dynamic Rcznge is similar to the ADC Sensor Range, except the indicator
reports
dynamic channel range as a percent rather than a fixed bit number. Clipping
indicates the
number of observations of clipping in the raw data. For a specific gain value,
the raw ADC
bit values cannot exceed a specific calculated value. Low F~equehcy Slope
(LowFf°eqSlope)
and Low Frequefzcy Ihte~cept (lowFreqlnt) use the first 10 points of the power
spectral
density calculated from the raw data and perform a simple linear regression to
obtain the
intercept and slope in the frequency-amplitude domain. SNR is the signal to
noise ratio
l0 observed in each specific data channel. A power spectral density is
calculated from the raw
uncalibrated vibration data. For each data channel, there are known
frequencies associated
with certain components. Examples include, but are not limited to, gear mesh
frequencies,
shaft rotation rates, and indexer pulse rates. SNR measures the rise of a
known tone
(corrected for operational speed differences) above the typical minimum
baseline levels in a
15 user-defined bandwidth (generally +/- 8 bins).
The Statistics (ST) algorithm 360 is associated with producing a plurality of
statistical
indicators 360a. The Root-Mec~ra-Square (RMS) value of the raw vibration
amplitude
represents the overall energy level of the vibration. The RMS value can be
used to detect
20 major overall changes in the vibration level. The Peak-To-Peak value of the
raw vibrating
amplitude represents the difference between the two vibration extrema. When
failures occur,
the vibration amplitude tends to increase in both upward and downward
directions and thus
the Peak-To-Peak value increases. The Skewness coefficient (which is the third
statistical
moment) measures the asymmetry of the probability density function (p.d.f.) of
the raw
25 vibration amplitude. Since it is generally believed that the p.d.f. is near
Gaussian and has a
26

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Skewness coefficient of zero, any large deviations of this value from zero may
be an
indication of faults. A localized defect in a machine usually results in
impulsive peaks in the
raw vibration signal, which affects the tails of the p.d.f. of the vibration
amplitude. The fourth
moment (Kuf°tosis) of the distribution has the ability to enhance the
sensitivity of such tail
changes. It has a value of 3 (Gaussian distribution) when the machinery is
healthy. Kurtosis
values larger than 3.5 are usually an indication of localized defects.
However, distributed
defects such as wear tend to smooth the distribution and thus decrease the
Kurtosis values.
The ST algorithm may be performed on the following vectors: AD raw
accelerometer
data, Alsignal average data, RS residual data, NB narrow band data, and EV
envelope data
and others, some of which are listed in 360b.
The Tone czndBase Ehe~gy algo~it7zna(TB) 362 uses tone energy and base energy.
Tone Efaergy is calculated as the sum of all the strong tones in the raw
vibration spectrum.
Localized defects tend to increase the energy levels of the strong tones. This
indicator is
designed to provide an overall indication of localized defects. "Strong tones"
are determined
by applying a threshold which is set based on the mean of all the energy
contents in the
spectrum. Any tones that are above this threshold are attributed to this
indicator. The Base
2o Enefgy measures the remaining energy level when all the strong tones are
removed from the
raw vibration spectrum. Certain failures such as wear, do not seem to affect
the strong tones
created by shaft rotation and gear mesh, the energy in the base of the
spectrum could
potentially be a powerful detection indicator for wear-related failures. Note
that the sum of
Tone Energy and Base Energy equals the overall energy level in the spectnun.
27

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
SI are miscellaneous shaft indicators. SOl (Shaft Order 1 in g) is the once-
per-rev
energy in the signal average, and is used to detect shaft imbalance. S02
(Shaft Order 2 in g)
is the twice-per-rev energy in the signal average, and is used to detect shaft
misalignment.
GDF (Gear detector fault) may be an effective detector for distributed gear
faults such as
wear and multiple tooth cracks, and is a complement of the indicator
signalAverageLl (also
known as gearLocalFault).
In addition to the specifically referenced vectors below, the SI algorithm
takes input
from the indexer zero-crossing vector (ZC).
to
The Defnodulation analysis (DM) 370 is designed to further reveal side band
modulation by using the Hilbert transform on either the narrow band signal
(narrow band
demodulation) or the signal average itself (wide band demodulation) to produce
the
Amplitude Modulation (AM) and Phase Modulation (FM) signals. The procedures
involved
15 to obtain such signals are:
Perform Hilbert transform on the narrow band signal (or signal average).
Compute the amplitude of the obtained complex analytic signal to obtain the
AM signal.
Compute the phase angles of the analytic signal to obtain the FM signal.
2o Compute the instantaneous amplitude of the analytic signal to obtain the
dAM
signal.
Compute the instantaneous phase angles of the analytic signal to obtain the
dFM signal.
28

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The DM algorithm is performed on the band passed filtered data at a frequency
of interest by
taking a Hilbert Window function of the frequency domain data and converting
the data back
to the time domain.
The Sideband Modulation (SM) 368 analysis is designed to reveal any sideband
activities that may be the results of certain gear faults such as
eccentricity, misalignment, or
looseness.
CIs included in 368a are DSMn. DSMn is an indicator that characterizes the
Degree of
Sideband Modulation for the nth sideband (n = l, 2, and 3). The DSMn is
calculated as the
sum of both the nth high and low sideband energies around the strongest gear
meshing
1o harmonic. As indicated in 368b, the SM algorithm is performed on the Fast
Fourier transform
vector (FF).
The Planetary Analysis (PL) 364 extracts.the Amplitude Modulation (AM) signal
produced by individual planet gears and compares the "uniformity" of all the
modulation
15 signals.
In general, when each planet gear orbits between the sun and the ring gears,
its vibration
modulates the vibration generated by the two gears. It is believed that when
one of the planet
gears is faulty, the amplitude modulation of that planet gear would behave
differently than
the rest of the planet gears. The procedure to perform this algorithm is to
obtain signal
2o averages for the input, output, and planet shafts. For each signal average:
Locate the strongest gear meshing harmonic.
Bandpass filter the signal average around this frequency, with the bandwidth
equals to twice the number of planet gears.
Hilbert transform the bandpass filtered signal to obtain the AM signal.
25 Find the maximum(MAX) and minimum(MIN) of the AM signal.
29

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Calculate the Planet Gear Fault (PGF) indicator as included in 364a
according to the equation PGF = MAX(AM) / MIN (AM).
The inputs to the PL algorithm are the raw accelerometer data (AD) and the
indexer zero-
crossing data (ZC).
The Zero-Crossing Indicators (ZI) algorithm 354 is performed on the zero-
crossing
vector (ZC). The zero crossing indicators may be determined as follows:
D~ = Ins+, - In ~ , j = O..N - 2 , the stored zero-crossing intervals
pulselntervalMean = Mean(D)
to The Shaft Indicators (S1) algorithm 358 calculates miscellaneous shaft
indicators
included in 358a. SOI (Shaft Order 1 ifa g) is the once-per-rev energy in the
signal average,
and is used to detect shaft imbalance. S02 (Shaft Order 2 in g) is the twice-
per-rev energy in
the signal average, and is used to detect shaft misalignment.
S03 (Shaft Order 3), is the three-per-rev energy in the signal average, and is
used to
detect shaft misalignment. The miscellaneous shaft indicators may also be
included in an
embodiment defined as follows:
p = nufnPathPairs
p-1
shaftPath2;
shaftRatio = p' 1 - driving
driven
shaftPath2~+,
f=O
p-I
indexPath~~
indexRatio = '-° -_ driving
p-' driven
ifadexPath~;+I
=o
driveRatio = ZndexRatio , pulsesUsed
shaftRatio

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
shaftSpeed = 60
pulselnterwalMeara ~ dr°iveRatio
f-esampleRate = shaft~SOpeed . pointsPer~Rev
RS = residual data,
Al = signal average,
sigraalAverageLl = P~p(Al)
Rms(Al)
FF = FFT of the signal average,
shaftOr~der~ = FF~ , j =1..3
gearDistFault = stdev(RS)
Stdev(Al)
As described elsewhere herein, gearDistFault (GDF) is an effective detector
for
distributed gear faults such as wear and multiple tooth cracks, and is a
complement of the
indicator signalAverageLl (also known as gearLocalFault).
In addition to the specifically referenced vectors below, the SI algorithm
takes input
from the indexer zero-crossing vector (ZC) and may also use others and
indicated above.
2o The following definitions for indicators may also be included in an
embodiment in
connection with the SI algorithm:
shaftPath is defined for the shaft descriptor
indexPath is the path of the shaft seen by the indexer used for signal
averaging
nunaPathPaiYS is the number of path pairs defined for shaftPatla and indexPath
pulses Used is the number of pulses used per revolution of the indexer shaft
pulselntervalMearz is the mean of the zero-crossing (ZC) intervals
pointsPer~Rev is the number of output points per revolution in the signal
average,
31

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The Bearing Energy (BE) algorithm 376 performs an analysis to reveal the four
bearing defect frequencies (cage, ball spin, outer race, and inner race
frequencies) that usually
modulate the bearing shaft frequency. As such, these four frequencies are
calculated based on
the measured shaft speed and bearing geometry. Alternatively, the four
frequency ratios may
be obtained from the bearing manufacturers. The energy levels associated with
these four
frequencies and their harmonics are calculated for bearing fault detection.
They are:
- Cage Efaergy: the total energy associated with the bearing cage defect
to frequency and its harmonics. Usually it is detectable only at the later
stage
of a bearing failure, but some studies show that this indicator may increase
before the others.
- Ball Energy: the total energy associated with the bearing ball spin defect
frequency and its harmonics.
15 - Outer Race Energy: the total energy associated with the bearing outer
race
defect frequency and its harmonics.
- Inner Race Enef gy: the total energy associated with the bearing inner race
defect frequency and its harmonics.
The Total Enef gy indicator gives an overall measure of the bearing defect
energies.
In one embodiment, one or more algorithms may be used in determining a CI
representing a score quantifying a difference between observed or actual test
distribution data
and a normal probability distribution function (PDF) or a normal cumulative
distribution
function (CDF). These one or more algorithms may be categorized as belonging
to a class of
algorithms producing CIs using hypothesis tests ("hypothesis testing
algorithms") that
32

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
provide a measure of difference in determining whether a given distribution is
not normally
distributed. These hypothesis testing algorithms produce a score that is used
as a CI. The
score may be described as a sum of differences between an observed or actual
test
distribution function based on observed data and a normal PDF or normal CDF.
An
algorithm may exist, for example, based on each of the following tests: Chi-
Squared
Goodness of fit (CS), Kolmogorov-Smirnov Goodness of fit (KS), Lilliefors test
of
normality, and Jarque-Bera test of normality (JB). Other embodiments may also
include
other algorithms based on other tests for normality, as known to those of
ordinary skill in the
art. The hypothesis tests compare the test distribution to the normal PDF, for
example as
to with CS test, or the normal CDF, for example as with the KS and Lilliefor
tests.
What will now be described is an example in which the CS test is used in
determining
a score with a test distribution of observed actual data. In this example, the
test distribution
of observed data forms a Gamma (5, 20) distribution function, having and alpha
value of 5
15 and a beta value of 20. The mean of this Gamma(5,20) distribution is alpha
* beta having a
variance of alpha * betaz . The Gamma (5,20) distribution function is a tailed
distribution
which graphically is similar to that of a normal distribution.
Referring now to Figure 19, shown is an example of a graphical representation
400 of
20 observed data.
Referring now to Figure 20, shown is an example of a graphical representation
410 of
the normal CDF and the Gamma (5,20) CDF of random data. Referring now Figure
21,
shown is an example of a graphical representation 420 of the difference
between the normal
25 CDF and the Gamma (5,20) CDF.
33

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
In one embodiment, if there are 1000 test samples used in forming a single
CDF, the graphical representation, for example, in Figure 21 represents
differences in 1000
instances where the difference between the expected value (Normal CDF) and the
maximum
deviation of the (in this case defined as the score) observed gamma CDF can
exceed some
critical value. The critical value is that statistic which represents some
predefined alpha error
(the probability that the test indicates the distribution is not normal when
in fact it is normal -
this is typically set at 5%.) If the score exceeds the critical value, the
distribution is said to be
not normal statistic. The score is the maximum deviation from this statistic
or alpha value.
to
It should be noted that the sensitivity or goodness of the test increases as
the number
of samples or instances (degrees of freedom "n") increases approximately as
the square root
of "n". For example, in the case where 1000 instances or samples are used such
that n=1000,
the sensitivity or ability of this CI to be used in detecting gear faults, for
example, is roughly
15 31 times more powerful than kurtosis in identifying a non normal
distribution.
As another example, in the algorithm using the CS test, the normal PDF is
used.
Refernng now to Figure 22, shown is a graphical representation 430 of the
normal PDF and
the PDF of the Gamma (5,20) distribution. The representations of Figure 22 are
drawn as
2o continuous lines rather than discrete intervals.
Referring now to Figure 23, the quantities of the x-axis represented in Figure
22 are
shown in another representation 440 as being divided into discrete bins,
intervals, or
categories. For example, there may be 4 bins or intervals between any two
integer quantities.
25 Between 0 and l, bin 1 includes values between [0,0.25), bin 2 includes
values between
34

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
[0.25, 0.50), bin 3 includes values between [.050,0.75) and bin 4 includes
values between
[0.75, 1.0). For each bin, determine the number of observed and expected
values, and their
difference. Square each of the differences for each bin and then add all the
differences and
divide by the expected value for each bin. The CS test which sums all the
differences for each
category divided by the expected value for each category represented as:
~~fi-ei)z
;_, ei
for k categories or bins, k-1 degrees of freedom, fi is observed data and ei
is expected data
value or number in accordance with a normal distribution.
l0
For each bin, take the difference between the observed and expected
observation.
Square this value and divided by expected number of observation. Sum over all
bins. The
statistic, the critical value is the x2 at k-1 degrees of freedom may be, for
example, 90.72
which is much greater than the .05 alpha value of a x2, which is 54.57 for 39
degrees of
freedom or 40 categories/bins. Thus, the observed data in this example as
indicated by the
statistic is not normally distributed. Figure 24A represents graphically a
difference between
observed and expected values for each bin or interval of Figure 23.
It should be noted that the foregoing algorithms provide a way of measuring
both the
skewness and kurtosis simultaneously by comparing the PDF or CDF of the test
distribution
against the PDF/CDF of a standard normal distribution in which a score is used
as a CI as
described above.
As known to those of ordinary skill in the art, other algorithms belonging to
the
hypothesis testing class may be used in computing CIs. The particular
examples, algorithms,

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
and tests selected for discussion herein are representations of those that may
be included in
the general class.
What will now be described is another algorithm that may be used in
determining a CI
in an embodiment of the system of Figure 1. This may be referred to as an
impulse
determination algorithm that produces a CI indicating an amount of vibration
that may be
used in detecting a type of fault. The impulse determination algorithm takes
into account the
physical model of the system. One type of fault that this technique may be
used to detect is a
pit or span on either: gear tooth, inner bearing race, outer bearing race or
bearing roller
to element. This technique uses a model designed to detect this type of fault
where the model is
based on knowledge of the physical system. For example, if there is a pit or
spall on a
bearing, this may produce a vibration on a first bearing which may further add
vibrations to
other components connected to or coupled to the bearing.
In one embodiment, a model can be determined for a particular configuration by
using
configuration data, for example. In one configuration, for example, a signal
received at a
sensor may be a superposition of gear and bearing noise that may be
represented as a
convolution of gear/bearing noise and a convolution of the Gear/Bearing signal
with the
2o gearbox transfer function. Given this, if one type of fault is a pit or
spall on either a: gear
tooth, inner bearing race, outer bearing race or bearing roller element, a
model that is
designed to look for this type of fault can take advantage of knowledge of the
physical
system.
The impulse determination algorithm uses Linear Predictive Coding (LPC)
techniques. As known to those skilled in the art, LPC may be characterized as
an adaptive
36

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
type of signal processing algorithm used to deconvolute a signal into its base
components. In
the case of a pit/spall fault, the base signal components are an impulse train
generated by the
fault hitting a surface (e.g gear tooth with geartooth, inner race with roller
element, etc) and
the bearing/case transfer function. The bearing, gear and case have there own
transfer
functions. Convolution here is transitive and multiplicative. As such, LPC
techniques may
be used to estimate the total convolution function of the total vibration that
may be produced.
For example, in this arrangement, the total amount of vibration representing
the total
impulse signal generated by a configuration may be represented as:
to
~irrapulse~ D f (Gear) ~ f (Bearing) ~ f (Case) --- ~impulse~ ~ 'f'(Gear) ~ f
(Bearing) ~ f (Case)
in which ~ represents the convolution operation.
It should also be noted that convolution is a homomorphic system such that it
is
monotonically increasing and that logarithmic transformations hold. Thus the
relationship of
c = a*b also holds for Log c = Log a + Log b. A "dual nature" of convolution
is used in
2o
following representations to equate operations using convolution in the time
domain to
equivalent multiplication operation in the frequency domain.
If "y" represents the total response of all elementary responses, and "h"
represents the
response of the system for a series of elementary input impulses "imp" such
that y is the
convolution of imp and h, then this may be represented as:
y=in2p~h
37

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
and then converting "y" and "h" each, respectively, to the frequency domain
represented as "Y" and "H", as may be represented by the following:
Y = ~~Y~~ H = ~~h~
taking the Fourier transform (FFT) of each where H represents the transfer
function.
The convolution in the time domain may be equated to a multiplication in the
frequency
domain represented as:
Y=IMP~H
in which IMP is the Fourier transformation of imp into the frequency domain.
Above, imp is
in the time domain.
The convolution in the time domain is equivalent to multiplication in the
Frequency
Domain. Refernng to the homomorphic property of convolution, it follows that:
log(Y) = log(IMP)+ log(H),
tlae~efone
log(IMP) = log(Y) - log(H),
IMP = exp(log(Y)- log(H))
and fznally
imp = s-' (IMP)
Using the foregoing, the system transfer function "H" may be estimated for the
Gear/Bearing and Case to recover the impulse response allocated with a Gear or
Bearing
2o pit/spall fault. The estimation of this transfer function "H" may be
accomplished using
Linear Predictive Coding (LPC) techniques. LPC assumes that the Transfer
Function is a
FIR filter, and as such, the auto-correlation of the time domain signal may be
used to solve
for the filter coefficients in a minimum sum of square error sense.
38

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Using the LPC model, there is an impulse that is convoluted with a FIR filter,
such
that:
Y[n] - a~xLn -1]+ aZx[n - 2~+ a3x[n - 3]+ . ..
LPC techniques may be used to estimate the coefficients a=(al .. an) for an
orderp in a
minimum sum of square error sense, n~+l. The standard least squares error
estimators may
to be used, wherein y = y[1, 2, ... n], and x is the time delayed signal, in
which:
x[n 1, 2,..n-
- n p]
-
x= x[n-2,n-3,n-p-1]
where a = (xT x)-1 xT y. These values for al .. an may be used with the
following equation:
Yt~ar = ax, b = (y-y,,~t)Z and the estimator of error B is: E~II b.
Y may also be expressed as:
Y = FFT(y[ 1, 2, .. n])
in which y[ l ..n] are values in the time domain expressed in the frequency
domain as a
Fourier transform of the time domain values. Y represents current time vector
measurements
in the frequency domain.
In terms of a and B, the tr ansfer function H may be estimated and represented
as a/B,
(freq. Domain). Note that "a" is a vector of the values al .. . an obtained
above.
The homomorphic property of convolution as described above may be used to
estimate the impulse as represented in:
39

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
IMP = exp(log(Y) - log(H)) IMP Equation
If there is no fault, the impulse, for example, may be characterized
as "white noise". As the fault progresses, the impulse or the value of H
becomes larger. The
CI is the power spectral density at a bearing passing frequency for a bearing
fault, or a mesh
frequency for a gear fault. Other CIs based on the foregoing value may be a
"score" of the
Lilifers test for normality, or other such test.
to
In the foregoing, a pit or spall may cause a vibration or tapping.
Subsequently, other
elements in contact with the ball bearing may also vibrate exhibiting behavior
from this initial
vibration. Thus, the initial vibration of the pit or spall may cause an
impulse spectrum to be
exhibited by such a component having unusual noise or vibration.
The value of IMP as may be determined using the IMP Equation above represents
the
impulse function that may be used as a "raw" value and at a given frequency
and used as an
input into an HI determination technique. For example, the 1MP at a particular
frequency,
since this the spectrum, determined above may be compared to expected values,
such as may
be obtained from the stored historic data and configuration data. An
embodiment may also
take the power spectrum of this raw impulse spectrum prior to being used, for
example, as
input to an HI calculation where the power spectrum is observed at frequencies
of interest,
such as the inner race frequency. For example, if the impulse function is
within some
predetermined threshold amount, it may be concluded that there is no fault.
What is shown in the Figure 24B and Figure 24C are relative to a healthy
system,
such as a main gearbox, for example, such as in connection with a planetary
race fault of an
SH-60B U.S. Navy Helicopter built by Silorsky.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The Figure 24B representation 700 shows an impulse train in the frequency
domain of
the healthy system.
It should be noted that an embodiment may estimate the transfer function H
using
LPC using different techniques. An embodiment may estimate the transfer
function H using
an autocorrelation technique(AutoLPC). An embodiment may also estimate the
transfer
function H using a covariance technique (CovLPC). Use of autoconrelation may
use less
to mathematical operations, but require more data than using the covariance.
Alternatively, use
of the covariance technique may use more mathematical operations but require
less data. As
the amount of available data increases, the autocorrelation LPC result
converges to the
covariance LPC result. In one example, data samples are at 100KHz with 64,000
data points
used with the autocorrelation technique due to the relatively large number of
data points.
Figure 24C representation 710 shows the data of 700 from Figure 24B in the
time
domain rather than the frequency domain.
2o Figure 24D representation 720 shows the power spectral density of the above
figures
as deconvolved time data of frequency v. dB values in a healthy system.
The foregoing Figures 24B-24D represent data in a graphical display in
connection
with a healthy system. Following are three additional graphical displays shown
in Figures
24E-24G in connection with an unhealthy system, such as a starboard ring
channel which
exhibit data that may be expected in connection with a pit or span fault.
41

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Figure 24E, representation 730, illustrates an impulse train as may be
associated with
an unhealthy system in the time domain. Figure 24F, representation 740,
illustrates a
graphical display of the impulse train in the frequency domain.
In Figure 24G, shown is an illustration 740 is a graphical representation of
the power
spectrum of the impulse train represented in connection with the other two
figures for the
unhealthy system identified by a period impulse train associated with an inner
race bearing
fault. In this example, a spike may be viewed in the graphical display as well
as the
to harmonics thereof.
It should be noted that other algorithms and CIs in addition to those
described herein
may be used in producing CIs used in techniques in connection with HIs
elsewhere herein.
What will now be described is one embodiment in which these CIs may be used.
Referring now to Figure 25, shown is a flow chart of steps of one embodiment
for
determining the health of a part as indicated by an HI. At step 502, raw data
acquisition is
performed. This may be, for example, issuing appropriate commands causing the
VPU to
2o perform a data acquisition. At step 504, the raw data may be adjusted, for
example, in
accordance with particular configuration information producing analysis data
as output. It is
at step 504, for example, that an embodiment may make adjustments to a raw
data item
acquired as may be related to the particular arrangement of components. At
step 506, data
transformations may be performed using the analysis data and other data, such
as raw data.
The output of the data transformations includes transformation output vectors.
At step 508,
42

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
CIs are computed using the analysis data and transformation vector data as may
be specified
in accordance with each algorithm. At step 510, one or more CIs may be
selected. Particular
techniques that may be included in an embodiment for selecting particular CIs
is described
elsewhere herein in more detail. At step 512, CIs may be normalized. This step
is described
in more detail elsewhere herein. At step 514, the selected and normalized CIs
are used in
determining HIs. Particular techniques for determining HIs are described in
more detail
elsewhere herein.
In an embodiment, due to the lengthy processing times, for example, in
executing the
different algorithms described herein, HI computations may not be executed in
real time.
Rather, they may be performed, for example, when a command or request is
issued, such as
from a pilot or at predetermined time intervals.
The hardware and/or software included in each embodiment may vary. in one
embodiment, data acquisition and/or computations may be performed by one or
more digital
signal processors (DSPs) running at a particular clock speed, such as 40MHz,
having a
predetermined numerical precision, such as 32 bits. The processors may have
access to
shared memory. In one embodiment, sensors may be multiplexed and data may be
acquired
in groups, such as ~. Other embodiments may vary the number in each group for
data
sampling. The sampling rates and durations within an acquisition group may
also vary in an
embodiment. Data may be placed in the memory accessed by the DSPs on
acquisition. In one
embodiment, the software may be a combination of ADA95 and machine code.
Processors
may include the VPU as described herein as well as a DSP chip.
43

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
What will now be described are techniques for normalizing CIs in connection
with
determining HIs providing more detailed processing of step 512 as described in
connection
with flowchart 500.
Transmission error (T.E.) depends upon torque. Additionally, vibration depends
upon
the frequency response of a gear. As such, the CI, which also depends upon
T.E. and
vibration, is a function (generally linear) of torque and rotor speed (which
is frequency), and
airspeed as this may change the shape of the airframe. Thus, techniques that
may be used in
connection with determining the "health state" or HI of a component may
normalize CIs to
l0 account for the foregoing since HIs are determined using CIs.
For each bearing, shaft and gear within a power train, a number of CIs may be
determined. An embodiment may compare CI values to threshold values, apply a
weighting
factor, and sum the weighted CIs to determine an HI value for a component at a
particular
15 time.
.Because data acquisitions may be made at different torque (e.g. power
setting) values, the
threshold values may be different for each torque value. For example, an
embodiment may
use 4 torque bands, requiring 4 threshold values and weights for each CI.
Additionally, the
coarseness of the torque bands will result in increased, uncontrolled system
variance.
2o Alternatively, rather than use multiple threshold values and have an
uncontrolled variance, an
embodiment may use a normalization technique which normalizes the CI for
torque and rotor
RPM (Nr), and airspeed, expressed as a percentage, for example, in which a
percentage of
100% is perfect. Use of these normalized CIs allows for a reduction of
configuration such
that, for example, only one threshold is used and variance may also be
reduced.
44

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The normalization technique that will now be described in more detail may be
used in
connection with methods of HI generation, such as the non-linear mapping
method and the
hypothesis testing method of HI generation that are also described in more
detail elsewhere
herein.
It should be noted that a deflection in a spring is linearly related to the
force applied to
the spring. The transmission may be similar in certain aspects to a large,
complex spring.
The displacement of a pinion and its corresponding Transmission Error (T.E.)
is proportional
to the torque applied. T.E. is a what causes vibration, while the intensity of
the vibration is a
to function of the frequency response (NY), where frequency is a function of
RPM. Thus,
vibration and the corresponding CI calculated using a data acquisition are
approximately
linearly proportional to torque, NY, (over the operating range of interest)
and/or airspeed
although at times there may be a linear torque*Nr interaction effect. For
example, gear box
manufacturers may design a gearbox to have minimum T.E. under load, and a
graphical
15 representation of T.E. vs. Torque is linear, or at least piece wise linear.
It should be noted
that test data , for example used in connection with a Bell helicopter H-1
loss of Tube test,
shows a relationship between CI and torque suggesting linearity. Additionally,
tests show that
airspeed is also relevant factor. Other embodiments may take into account any
one or more
of these factors as well as apply the techniques described herein to other
factors that may be
20 relevant in a particular embodiment or other application although in this
example, the factors
of torque, airspeed and Nr are taken into account.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
An equation representing a model minimizing the sum of square error of a
measured
CI for a given torque value in a healthy gear box is:
CI = Bo + B, * Torque + BZ Nr + B3 Air speed + T.E... (Equation 1)
The order of the model may be determined by statistical significance of the
coefficients of Equation 1. In the previous equation, the T.E. of a "healthy"
component may
have, for example, a mean of zero (0) with some expected variance. It should
be noted that if
the model fits well for the lower order. Higher order coefficients are not
required and may
actually induce error in some instances. The following example is built as a
first order
model, higher orders may be solved by extension of that explained in the first
order model.
to This model, written in matrix format is: y = B x where
Ch 1 t~ Nn, Airspeedl
y = CI_.. B = ~Bo B~...B~~ ~ arzd x = 1 t... NR... Airspeed...
Cln 1 tn Nan Airspeedn
Each of the CIs included in the vector y is a particular recorded value for a
CI from previous
data acquisitions, for example, as may be stored and retrieved from the
collected data 18.
Also stored with each occurrence of a CI for a data acquisition in an
embodiment may be a
corresponding value for torque (t), Nr, and Airspeed. These values may also be
stored in the
collected data 18.
The model coefficients for B may be estimated by minimizing the sum of square
error
2o between the measured CI and the model or estimated CI using the observed
performance
data. Solving the foregoing for the unbiased estimator of B = (xTx)-' xTy .
The variance of B
46

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
is: T~az~(B) = E(b -B)(b -B)T = ~Z (xTx)-' where b is an unbiased estimator of
B. The
unbiased estimator of 62 is s2: s Z = eT a - (Y -Y)T (Y - Y) _ YT Y -bT xT Y
zz-p-1 zz-p-1 zz-p-1
In the vector B from y=xB, coefficient Bo represents the mean of the data set
for a
particular component which, for example, may be represented as an offset
value. Each of the
other values B1 .. . Bn are coefficients multiplied by the corresponding
factors, such as
airspeed, torque, and Nr.
The foregoing B values or coefficients may be determined at a time other than
in real-
l0 time, for example, when flying a plane, and then subsequently stored, along
with
corresponding X information, for example, in the collected data store 18.
These stored values
may be used in determining a normalized CI value for a particular observed
instance of a
CIobs in determining an HI. The normalized CI may be represented as:
Cljlarntalized ='f.E. = CIobs - (B~ X)
where CIobs represents an instance of a CI being normalized using previously
deternined
and stored B and x values. Threshold values, as may be used, for example, in
HI
determination, may be expressed in terms of multiples of the standard
deviation Warning =
Bo+ 3*a2(xtx)-1, Alarm = Bo+ 662 (xtx)-I. It should be noted that a covariance
that may be
determined as:
2o
E = s2 (xtx)-1 where s2 is calculated as noted above.
As described elsewhere herein, the foregoing techniques are based upon a
healthy
gear characterized as having noise that is stationary and Gaussian in which
the noise
approximates a normal distribution.
47

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
What will now be described are techniques that may be used in determining an
HI
using the normalized CI values as inputs. In particular, two techniques will
be described for
determining an HI. A first technique may be referred to as the non-linear map
technique.
The second technique may be referred to as the hypothesis test method of HI
generation. It
should be noted that CI values other than normalized CI values may be used in
connection
with HI determination techniques described herein.
It~ should be noted that an embodiment may use CI values that are not
normalized in
to connection with the HI determination techniques described herein. In this
instance, multiple
torque bands may be used, one for each CI or group of CIs belonging to
different torque
bands. Additionally, a larger covariance matrix may be used as there may be a
larger
variance causing decrease in separation between classes.
15 For any generic type of analysis (gear, bearing, or shaft), a subset of the
diagnostics indicators or CIs is selected: The CIs which are best suited to
specify the fault
indication may be developed over time through data analysis. Faults may be
calculated at the
component level and an HI may be calculated for a given component. If there is
a
component fault, then there is a sub-assembly.fault, and therefore a drive
train fault.
Following is a description of a non-linear mapping methodology for determining
an
HI. Given a set of component indicators I1, I2, I3, ...IN, choose the desired
subset of K
indicators such that K <= N. For the chosen group of indicators, let WTi
define the weight of
the ith indicator, Wi the warning threshold, and Ai the alarm threshold. Then
apply the
following processing to the set of chosen indicators.
48

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Health Indicator Contribution Description
for XX = 1:K /* cycle through all K indicators in subset */
If I[XX] < Wi /* if less than warning level Wi, assign 0 */
s Hi contribution = 0
elseif Wi * Ii < Ai
Hi contribution =1 * Wi
else
1 o end
end
Hi contribution = 2*Wi
In the foregoing pseudo-code like description, each indicator or CI is
weighted and
contributes a portion to the HI determination. Subsequently all the Hi
contributions for the
15 selected CIs are summed and may be compared to threshold values for
determining one of
two possible outcomes of "healthy" or "not healthy".
Consider the following example table of information for a selected subset of 9
CIs
along with threshold and weight values. It should be noted that in an
embodiment, any one or
2o more of the values for weights, warning and alarm values may be modified.
CI Value Warning Alarm Weight HI contribution
No. Level Level
I2 3.26 3.5 4.0 1.0 0.0
25 I3 3.45 . 3.0 3.5 1.0 1.0
I6 7.5 6.0 8.0 1.4 1.4
I9 0.88 0.5 0.75 0.9 1.8
I14 4.2 3.5 4.5 1.0 1.0
I17 4.7 3.5 4.5 0.9 1.8
3o I22 5.2 2.0 4.0 1.1 2.2
I23 4.4 3.5 4.5 1.2 1.2
I24 18.9 10.0 20.0 1.0 1.0
Using the foregoing example and values, the sum of the HI contributions is
11.4.
49

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Applying the Health Indicator Contribution technique as set forth in the
foregoing pseudo-
code like description, I2, with a value of 3.26, is below the warning
threshold, so the
contribution to the index is 0. Indicator I3 has a value of 3.45, which
contributes a 1 toward
the index since the weight value is also 1. However, Indicator I6 contributes
a 1.4 to the
index because it crosses the warning level (contributing a value of 1 to the
index) while being
weighted by a factor of 1.4.
In the foregoing example, if no indicators were in alarm, the sum of HI
contributions
would be zero and if all indicators were in alarm, the sum would be 19, the
worst fault case
l0 represented by this detector scheme. The HI may be represented as a value
of 1 for healthy
and 0 for not healthy as associated with a component represented by the
foregoing CI values.
The HI may be determined by dividing 11.4/19, the maximum of worst case
outcome
to obtain 0.6. This overall health index output ratio can then be compared to
another final
15 output threshold, where normal components produce HIs, for example, less
than 0.5; values
between 0.5 and 0.75 represent warning levels, and values over 0.75 represent
alarm.
It should be noted that the weights may be determined using a variety of
different
techniques. The weights of each CI may be determined using any one or more of
a variety of
20 techniques. One embodiment may determine weights for the CIs as:
1
eige~c _ values _ of _ the - cov ariahce _ matrix
It should be noted that other threshold values may be used in HI determination
and
may vary with each embodiment.

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
In one embodiment, using the normalized CI described elsewhere herein with the
non-
linear mapping technique, the threshold values may be represented as: Warning
= Bo+
3*62(xtx)-1, Alarm = Bo+ 6*~2(xtx)-1, where Bo may represent a mean or average
coefficient
as included in the B vector being solved for in the equations described in
connection with CI
normalization. In the foregoing example, the Warning threshold is 3 standard
deviations and
the Alarm level is 6 standard deviations. It should be noted that other
threshold values may
be used in and may vary in accordance with each embodiment.
1o What will now be described is a second technique that may be used in
determining
HIs using CIs, in particular, using normalized CIs.
The technique for HI determination may be referred to as Hypothesis testing
technique for HI determination which minimizes the occurrence of a false alarm
rate, or
15 incorrectly diagnosing the health of a part as being included in the alarm
classification when
in fact the part is not in this particular state. In one embodiment, three
classes of health
indication may be used, for example, normal, warning and alarm classifications
with alarm
being the least "healthy" classification. Other embodiments may use the
techniques described
herein with a different number of classes. As described elsewhere herein, the
class of a part
2o indicating the health of the part may be determined based on measured
vibrations associated
with the part. Additionally, the technique described herein may use a
transformation, such as
the whitening transformation to maximize the class distributions or separation
of values thus
decreasing the likelihood or amount of overlap between the classes. In
particular, this
maximization of class separation or distance attempts to minimize the
misclassification of a
25 part. A description of the whitening transformation used in herein in
following paragraphs
51

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
may be found, for example, in "Detection, Estimation and Modulation Theory",
Harry L.
Van Trees, 1968, John Wiley & Sons , New York Library of Congress Catalog Card
Number
67-23331.
Using the Hypothesis Testing method of HI generation, the HI or classification
h(X)
of a vector of normalized CI values denoted as X may be determined in which,
as discussed
elsewhere herein in more detail, X may be normalized.Using the hypothesis
testing
technique, a determination is made as to which class (normal, warning or
alarm) X belongs.
In our instance, there are three classes. However, a first determination using
the hypothesis
l0 testing may be performed using a first class corresponding to normal, and a
second class
corresponding to not normal. If the determination is normal, then testing may
stop.
Otherwise, if determination is made that the testing results are "not normal",
a further or
second determination using the hypothesis testing may be performed to
determine which "not
normal" class (alarm or warning) X belongs. Thus, the hypothesis testing
technique may be
15 performed more than once in accordance with the particular number of
classes of an
embodiment. For three classes, there are two degrees of freedom such that if
the sample X is
not from A or B classes, then it is from Class C.
X may belong to class c~~ or co,, such that: q1 (X) <> qZ (X) (the notation <>
means
~2
20 that if ql(X) is greater than q2(X), choose class 2, coa, or if ql(X) is
less than q2(X), choose
class 1, w1.) In the foregoing, q; is the a posteriori probability of c~i
given X, which can be
computed, using Bayes theorem in which q1= P;pt(X)lp(~, where p(X) is the
mixed density
function.
52

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The mixed density function is the probability function for all cases where q;
is the
unconditional probability of "i" given the probability of "i" conditioned on
the mixed density
function.
Substituting the foregoing representation of each q1 and q2, since p(X) is
common to
(~ ~'' P
both, now: P, p, (X) <>Pz pz (X) or as a Iikelihoodfunction as .~(X) = p' <> p
. The
pz (~ °'z i
likelihood ratio is a quantity in hypothesis test. The value P2/PI is the
threshold value. hi
some instances, it may be easier to calculate the minus log likelihood ratio.
In this case, the
decision rule becomes (e.g. now called the discriminate function):
to h(X) _-ln.~(X) _-Inp,(X)+Inpz(X)<>Inpz
~= P~
Assume that the p~(X)'s are normally distributed with mean or expected values
in
vectors M~- and covariance matrix ~1 . This assumption may be determined
without loss of
generality in that, any non-normal distribution can be whitened, as with the
whitening
transformation described elsewhere herein, with the appropriate power
transform, or by
increasing the sample size to the point where the sample size is very large.
Given this, the
decision rule becomes:
~Z(x) _ -In e(x)
Equation E1
- 1 (x-W )T~i~(X -M,)- 1 (~'-MZ)T~2'(-~'-Mx)'H 1 In ~~'~o In Px
2 2 2 IE 2 I w P,
Recall that maximization of distance between the two classes is desired to
minimize
the chance of a false alarm or misclassification of a part as broken when it
is actually normal.
A function Z is defined as Z = X M, (e. g. a shift where X is the measured CI
data and
M is the mean CI values for a class), so that:
53

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
d (z) = ZT ~-1Z (this distance is the fa dimensional distance between two
distributions).
Note that E represents the covariance. It may be determined that a particular
Z maximizes the
distance function, subj ect to ZTZ = I, the identity matrix.
Using a standard Lagrange multiplier, ~, to find the local extrema (e.g. the
maximum)
a partial derivative is obtained with respect to Z in the following:
/aZ {ZT ~-~Z f~~ZT Z - I)) = 2E-'Z - 2~Z
where E is the eigenvector of X,
1o which may then be set to zero to find the extrema and solving for Z:
E -' z = ~ z or E z = a, z where a, _ ~ ~ ~ ~ order that a non-null Z exits,
,must be chosen to
satisfy the determinant: ~E - ~,z~= o .
Note that ~, is the eigenvalue of X and E is the corresponding eigenvector. E
is a
symmetric ya x n matrix (e.g. a covariance matrix), there are h real
eigenvalues (Aj...l~") and fa
real eigenvectors ~~ ,., ~". the characteristic equation is: E~=~1~, and ~T~ =
I where ~ is an
n x n matrix consisting of n eigenvectors and n is a diagonal matrix of
eigenvalues (e.g. the
eigenvector matrix and eigenvalue matrix, respectively).
Y, representing the coordinated shifted value of X, may be represented as:
2o Y = ~T X,
having a covariance matrix of y, E y = ~ T ~ x ~ _ ~ where EX represents the
covariance of the vector of matrix x . Continuing, the whitening
transformation may be
defined such that:
Y = ~ I/2~T X = (~~ 1/2)TX' ~y = ~ 1I2~T ~x ~~ 1/2 = ~ 1/2 ~ 1/2 = I'
54

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Thus the transformation that maximizes that distance between distribution or
classes
is:
A =Ii n2~T as shown above.
Using this value of A, define
AT E 1 A = I, AT E2 A = K, and AT (M2-Ml) = L and
(E 1-1 EZ 1) -1 transformed to a diagonal matrix 1~ by A that may be
represented as:
EAT~ [A (I-K-i)AT]-' A=(I-K-i)-'
which may be substituted into the discriminate function defined above:
to
h(X)= 1 YTA-'Y-[(-K-'L)T ]Y+[-1 LTK-'L- 1 lnlKl-lnpZ ]
2 2 2 P
Thus, if the above is less than the threshold, for example; In (PZ/P1), then
the
component is a member of the normal or healthy class. Otherwise, the component
is
classified as having an HI in the broken class, such as one of alarm or
warning. In the latter
case, another iteration of the hypothesis testing technique described herein
may be further
performed to determine which "brolcen" classification, such as alarm or
warning in this
instance, characterizes the health of the component under consideration.
2o In the foregoing technique for hypothesis testing, values, such as the a
poste~iofri
probabilities q1 and q2, may be obtained and determined prior to executing the
hypothesis
testing technique on a particular set of CI normalized values represented as X
above. As
known to those of ordinary skill in the art, Bayes theorem may be used in
determining, for
example, how likely a cause is given that an effect has occurred. In this
example, the effect is

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
the particular CI normalized values and it is being determined how likely each
particular
cause, such as a normal or broken part, given the particular effects.
It should be noted that operating characteristics of a system define the
probability of a
false alarm (PFA) and the probability of detection (PD). The transformation
used to
maximize the distance function optimizes the discrimination between classes.
However, the
threshold value selected given a discriminate function may be used in
determining the PD and
PFA. In some embodiments, the cost of a false alarm may be higher than the
cost of a missed
detection. In these instances, the PFA may be set to define threshold values,
and then accept
l0 the PD (e.g., a constant false alarm rate (CFAR) type of process). The
distance function is a
normal density function, based on the conditional covariance of the tested
values under
consideration. Given that, the PFA may be determined as: PF = P(HoHi), which
means the
probability that the sufficient statistic is greater than some threshold is
the integral of the
threshold to infinity of a normal PDF.
z
PFA = ~ p,~Ho (l ~ Ho )d~ _ ,~ 127C exp( x )'fix
2
where
the lower integral limit of
a = 1n( ~ ) l d + d l 2, and, as before d 2 = (M2 - M, )T Ei' (M2 - M, )
z
2o In this example, the threshold may be the In (P2/P1). This integration is
the
incomplete gamma function. Conversely, the probability of a detection (PD) is:
2
PD - J-m pl~H~ (l ~ H~ )dL = ~~ 12~c exp( ( 2) )'lx
but now
a=-ln(PZP)ld+dl2, and d'=(Ml-MZ)TEz'(M,-MZ)
56

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Note, the distance function is relative to the condition (e.g. Ho or H,) being
investigated.
Refernng now to Figure 26, shown is an example of a graphical illustration of
the
probability of a false alarm PFA represented by the shaded region A3 which
designates the
overlap between the distribution of class H0, denoted by the curve formed by
line Al, and
class Hl, denoted by the curve formed by line A2.
Referring now to Figure 27, shown is an example of a graphical illustration of
the
probability of an appropriate detection (PD) represented as area A4 as
belonging to class
to represented by Hl as represented by the curve formed by line A2.
Refernng now to Figure 2~, shown is a graphical illustration of a relationship
in one
embodiment between the PFA and PD and the threshold value. Note that as the
threshold
increases, the PD increases, but also the PFA increases. If the performance is
not acceptable,
15 such as the PFA is too high, an alternative is to increase the
dimensionality of the classifier,
such as by increasing the population sample size, n. Since the variance is
related by
1/sqroot(n), as n increases the variance is decreased and the normalized
distance between the
distributions will increase. This may characterize the performance of the
system. The
likelihood ratio test used herein is a signal to noise ratio such that the
larger the ratio, (e.g.,
20 the larger the distance between the two distributions), the greater the
system performance.
The process of taking an orthonormal transformation may be characterized as
similar to the of
a matched filter maximizing the signal to noise ratio.
Referring now to Figure 29, shown is an example of a graphical illustration of
how
the threshold may vary in accordance with the probability of determining class
Ho.
57

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
It should be noted that false alarm rate and detection rate are two factors
that may
affect selection of particular values, such as thresholds within a particular
system. In the
example embodiment described herein, false alarm rate is a determining factor,
for example,
because of the high cost associated with false alarms and the fact that they
may corrode
confidence when a real fault is detected. It should be noted that other
embodiments and other
applications may have different considerations. Further in this example of the
system of
Figure l, certain factors may be considered. An acceptable false alarm rate,
for example,
such as 1 false alarm per 100 flight hours, is established. An estimate of the
number of
collection opportunities per flight hours may be determined, such as four data
collections. A
to number of HIs may be selected for the system, such as approximately 800. A
confidence
level may be selected, such as that there is a 90% probability that a false
alarm rate is less
than 1 per 100 flight hours.
In this example, it should be noted that each HI is a an independent
classification
event such that the law of total probability may give the system alarm rate
using the
15 foregoing:
System PFA =1/(100 * 4 * 800) = 3.1250 * 10 -6.
It should also be noted that in the foregoing, when the covariance of two
classes is
approximately the same, or for example, unknown for a class, the logarithm
likelihood ratio
test for classification may be simplified in that the model may be reduced to
a linear rather
20 than quadratic problem having the following model:
(M2 _ M~ )T ~-~ ~ + ~ (Mi ~-iM1 _ Mz ~-iM2 )~j'~ pz
~z
58

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
If the covariance is whitened, the model simplifies further (assuming the
appropriate
transformation is made to the means and measured values).
~MZ -M, )T ~Y+ 1 ~ (Mi Ml -MZ MZ )<>~ln Pz
What will now be described are techniques that may be used in comzection with
selecting a subset of CIs, such as selection of normalized CIs, for example,
under
consideration for use in determining a particular HI.
If we have a two or more classes (such as alarm, warning and normal
classifications),
feature extraction, or determining which CIs to use in this embodiment, may
become a
problem of picking those CIs or features that maximize class separability.
Note that
1o reparability is not a distance. As described elsewhere herein, an
eigenvector matrix
transformation may be used in maximizing the distance between two functions or
distribution
classes. However, this same technique may not be applicable when some of the
information
(e.g. dimensionality) is being reduced. For example, in the following test
case, three features,
or CIs, are available, but only two are to be selected and used in determining
HI
classification. The distributions are:
1 -.5 .5 0 1 .7 .7 3
Covl = -.5 2 .~ , Ml = O , Govt = .7 2.5 1 , MZ = -1
.5 .~ 2.5 0 .7 1 2.5 3
When looking at the eigenvalues of the whitening transformation
(1.9311,3.0945, 0.4744),
the maximum distance of the distribution is an axis y (e.g. 2"a dimension, the
distribution was
whitened and the project dimension (e.g. x, y or z) was plotted), but this
axis has the
2o minimum reparability. Using this as one of the two features will result in
higher false alarm
59

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
rates than another feature. This may identify the importance of feature
selection in
maximizing the separability.
The problem of separability may be characterized as a "mixed" problem in that
differences in means may be normalized by different class covariance. If the
mean values are
the same, or the covariance are the same, techniques such as the Bhattacharyya
Distance may
be used to measure class separability. However, same mean or covariance values
may not be
likely and thus such techniques may not be applicable. Statistical tools
developed in
discriminant analysis may be used to estimate class reparability.
A measure of within class scatter may be represented as the weighted average
of the
L
1o class covariance: Sy" _ ~PE; , for each class I, where Pi is the
probability of the occurrence
!=1
of the covariance ~I'for that class. In one embodiment, there may be two
classes, such as
healthy or unhealthy. When considering the unhealthy status, for example, when
performing
a second round of hypothesis testing described herein, there may be alarm and
warning
classes.
A measure of between class scatter, Sb, may be represented as the mixture of
class
means:
L L
Sb =~I'~Mt -MoOM; -Mo)T ~ Mo =~I'M; .
=i
Note that Mo represents the mean or expected value of the classes and Mi- Mo
is a
difference or variation from the expected value for the classes under
consideration. The

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
formulation for a criteria for class reparability may result in values that
are larger when the
between class scatter is larger, or when the within class scatter is smaller.
A typical criteria
for this is J = diag(S",'Sb) , where In general, SW is not diagonal. One
technique takes the
whitening transformation of Sy" where AT S~"A = I, then define the whitening
transformation of
Sb as:
SbW = AT SGA.
Now taking the diagonal of the foregoing Sbw gives a better representation of
the
class reparability of each feature.
In summary, CIs may be selected in accordance with the technique described
above to
l0 obtain and examine the diagonals of the "whitened" Sb, represented as Sbw.
Let X be a
matrix where rows and columns represent different CIs having a covariance
matrix E. An
embodiment may use normalized CIs and select a portion of these for use. An
embodiment
may also use CIs however, those selected should belong to the same torque
band.
As described elsewhere herein, let 1~ represent the corresponding eigenvalue
matrix
and ~ as the corresponding eigenvector matrix for the CI matrix X. Then, A, as
described
elsewhere herein in connection with the whitening transformation, may be
represented as:
A=~1/2~T
where A is the transformation matrix that whitens the covariance E. If Sb is
defined as above
as the between mean covariance of the classes, the whitening matrix A may be
used to
61

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
normalize the differences and give a distance between the mean values of the
different
classes, such that
Sbw = AT SbA
where Sbw represents the "whitened" Sb. The diagonals of Sbw may then be
sorted in
descending order in which each diagonal represents an approximation of the
size of the
separation between features or CIs. Thus, selection of a subset of "n"
features or CIs from a
possible "m" maximum CIs included in X may be determined by selecting the "n"
largest
diagonals of the matrix Sbw. In particular, the diagonal entry 1,1 corresponds
to the first
column of the covariance matrix and the first CI in the vector X, entry 2,2 to
the second
to column of the covariance matrix and the second CI in the vector X being
considered, and so
on.
Once a particular HI is determined at a point in time, it may be desired to
use
techniques in connection with trending or predicting HI values of the
component at future
points in time. Techniques, such as trending, may be used in establishing, for
example, when
15 maintenance or replacement of a component may be expected. As described
elsewhere
herein, techniques may be used in determining an HI in accordance with a
vector of CI values
having expected CI values included in vector M; for a given HI classification,
i, having a
covariance matrix ~;. One technique uses a three state Kalman filter for
predicting or
trending future HI values.
The I~alman filter may be used for various reasons due to the particular
factors taken
into account in the embodiment and uses described herein. It should be noted
that other
systems embodying concepts and techniques described herein may also take into
account
62

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
other noise factors. In one embodiment, the Kalman filter may be preferred in
that it provides
for taking into account the noise of a particular arrangement of components.
There may be
noise corruption, such as indicated, for example, by the covariance matrices
described and
used herein. It may be desirous to filter out such known noise, such as using
the Kalman
filter, providing for smoothing of data values.
The Kalman filter provides a way to take into account other ap~iori knowledge
of the
system described herein. In particular, the health of a component, for
example, may not
change quickly with time. The difference between the health of a component at
a time t, and
to time t+delta may not be large. This technique may also be used in
connection with
determining future HIs of a particular part, for example, where the part is
old. A part may
have reached a particular state of relatively bad heahth, but still a working
and functional part.
The techniques described herein may be used with an older part, for example,
as well as a
newer part.
In the arrangement with the Kalman filter, state reconstruction may be
performed
using the Ricatti equation, as known to those of ordinary skill in the art.
The technique that
is described herein uses a three-state Kalman filter of HI, and the first and
second derivatives
thereof with respect to changes in time, denoted, respectively, dtz and dt3.
The Ricatti
2o equation in this instance uses a [1x3] vector of time values rather than a
single value, for
example, as may be used in connection with a single state Kalman filter.
What will now be described are equations and values that may be used in
determining a
future value of a particular HI. Let:
63

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
H = [1 0 0]
1 dt dt2/
= 0 1 d/ Lt
0 0 1
dt3 2 dt2 ~ t
Q = 262t dtz 2 dt 1
t 1 It
HI est
X = HI
HI
in which:
6 is the power spectral density of the system,
R is the measurement error,
P is the covariance,
Q is the plant noise,
H is the measurement matrix,
K is the Kalman gain and
~ is the state transition matrix.
H may be characterized as the Jacobian matrix. Since the value of a single HI
is
desired, only the first entry in the H vector is 1 with remaining zeroes.
There are n entries in
the n x 1 vector H for the n state Kalman filter. Similarly, the X vector
above is column
vector of 3 HI entries in accordance with the three-state Kalman filter. The
end value being
determined is the vector X, in this instance which represents a series of HI
values, for which
the first entry, HI est in the vector X is the one of interest as a projected
HI value being
determined Within the vector X, HI represents the first derivative of HI est
and HI
2o represents the second derivative of HI est. t represents the average amount
of time
between measurements or updating of the HI value. In other words, if dt
represents a
measurement or delta value in units of time between HI determinations, and
this is performed
for several instances, t represents the average of the delta values
representing time changes.
64

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
What will now be presented are equations representing the relationships
between the
above quantities as may be used in determining a value of X(1) for predicting
or estimating
an HI value at a future point in time given a current HI value.
Xt~~-, _ ~Xt-,~t-, (Equation T1)
P,~t-, =~P-,~l-,~T +Q (Equation T2)
K = P~t_,HT (HP,~,-,HT +R) (Equation T3)
Pit = (I -KH)P,~t-, (Equation T4)
Xt~, = Xt~t-, + K(HI - HXt~t-, ) (Equation TS)
Note that the subscript notation above, for example, such as "tit-1" refers to
determining a value of at a time t conditioned on the measurement at a time of
"t-1 ".
to Similarly, "tit" refers to, for example, determining an estimate at a time
"t" conditioned on a
measurement of time "t".
The current HI determined, for example, using other techniques described
herein, may
be input into Equation TS to obtain a projected value for HI est, the best
estimate of the
current HI. To project the expect HI "n" units of time into the future, input
the number of
units of time "dt" into ~ (as described above), and use the state update
equation (Equation
Tl) where now Equation T1 becomes: X t+at~t= ~ X tit . This allows the best
prediction of
HI est any number of units of time into the future where HI est is desired. It
should be noted
that as set forth above, the linear matrix operation such as ~ X is equivalent
to an integration
2o from t to dt of the state of X, where X represents the vector of HI values
set forth above.
Different values may be selected for initial conditions in accordance with
each
embodiment. For example, an initial value for P representing the covariance
may be

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
(1/mean time value between failures). An embodiment may use any one of a
variety of
different techniques to select an initial value for P. Additionally, since P
converges rapidly
to an appropriate value and the time between data acquisitions is small in
comparison to the
mean failure time, selecting a particularly good initial value for P may not
be as important as
other initial conditions. A value for a may be selecting in accordance with
apf~iori
information, such as manufacturer's data on the mean time between component
parts'
expected failure time. For example, for a transmission, the mean failure time
may be
approximately 20,000 hours. The spectral density may be set to (1/20,000)2. It
should be
noted that the failure rates may be generally characterized as an exponential
type of
to distribution. The mean time between expected failures is a rate, and the
variance is that rate
to the second power. R may also be determined using ap~iori information, such
as
manufacturer's data, for example, an estimated HI variance of manufacturer's
data of a
healthy component. Q may be characterized as the mean time between failures
and dt (delta
change in time between readings). As the value of dt increases, Q increases by
the third
15 power.
Input data used in the foregoing trending equations may be retrieved from
collected
data, for example, as may be stored in the system of Figure 1.
2o In determining HIs, for example, as in connection with the system of Figure
1 for
particular components, HIs may be derived using one or more CIs. In
calculating CIs, data
acquisitions may occur by recording observed data values using sensors that
monitor different
components. There may be a need for estimating data used in connection with CI
calculations, for example, in instances in which there may be too little or no
observed
25 empirical. For example, in connection with a power train, there may be a
need to obtain
66

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
estimated data, for example, for each bearing, shaft and gear within the power
train to
calculate CIs. However, insufficient empirical data may exist in connection
with gear or
bearing related measurements, such as, for example, those in connection with a
gear or
bearing related measurements, such as, for example, those in connection with a
gear or
bearing fault due to the rare occurrence of such events. In such instances,
mean and threshold
values may be derived using other techniques.
A CI may indicate a level of transmission error, for example, in which
transmission
error is a measure of the change in gear rigidity and spacing. Modeling
transmission error
to may allow one to gauge CI sensitivity and derive threshold and mean values
indicative of
gear/bearing failure. This transmission error modeling may be referred to as
dynamic
analysis. What will now be described is a technique that may be used to model
a gears to
obtain such estimated values. By modeling each gear pair as a damped spring
model with the
contact line between the gears, transmission error may be estimated. It should
be noted that
15 this model uses two degrees of freedom or movement. Other systems may use
other models
which may be more complex having more degrees of freedom. However, for the
purposes of
estimating values, this model has proven accurate in obtaining estimates.
Other
embodiments may use other models in estimating values for use in a system such
as that of
Figure 1.
Referring now to Figure 30, shown is an example of an illustration of a pair
of gears
for which a model will now be described. A force P at the contact gives linear
and torsional
response to each of the 2 gears for a total of four responses as indicated in
Figure 30. The
relative movement d at P is the sum of the 4 responses together with the
contact deflection
due to the contact stiffness s~ and the damping coefficient b~. This may be
represented as:
67

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
( jpz
d 1' 'l SP+j~h-yfzp~z + /kp+ jcxlp-Ipwz +
~w+ jtc~w-mt~2 + ~/kw+ je~w-Iw~z + EQUATION G1
~c+ j~c
in which:
sp is the linear stiffness of the pinion;
j is the square root of -1;
to co is the angular rate that may be obtained from the configuration file
(e.g., shaft rpm
* 60 ~= 2~ to obtain radians per second for the pinion driving the wheel);
by is the linear damping coefficient of the pinion;
mp is the mass of the pinion;
rp is the radius of the pinion;
kp is the angular effective stiffness of the pinion;
qp is the angular damping coefficient of the pinion;
Ip is the angular effective mass of the pinion;
2o sw is the linear stiffness of the wheel;
bw is the linear damping coefficient of the wheel;
mw is the mass of the wheel;
rp is the radius of the pinion;
68

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
kw is the angular effective stiffness of the wheel;
qw is the angular damping coefficient of the wheel;
Iw is the angular effective mass of the wheel;
sc is the linear stiffness of the contact patch where the two gears come into
contact;
be is the linear damping coefficient of the contact patch;
It should be noted that values for the above-referenced variables on the
rights hand
side of EQUATION G1 above, except for P (described below), may be obtained
using
to manufacturer's specifications for a particular arrangement used in an
embodiment. An
embodiment may include quantities for the above-referenced variables in units,
for example,
such as stiffness in units of force/distance (e.g., newtons/meter), mass in kg
units, and the
like.
The relative movement, d, is the T.E., so from d, the above-referenced
equation can
be solved for P, the tooth force. Deflection is the force (input torque
divided by the pinion
base radius) ~ the elastic deflection of the shafts, which may be used in
estimating P
represented as:
2o P = ( 1 /kp * rp) + ( 1 /sp) + ( 1 /sw) + ( 1 /kwrw)EQUATION G2
where the variables are as described above in connection with EQUATION Gl.
Using the above estimate for P with EQUATION Gl, the displacement, such as a
vibration
transmitted through the bearing housing and transmission case (which acts an
additional
transfer function), may be determined.
69

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
Referring again to Figure 30, shown is an example of an illustration of the
gear model
and the different variables used in connection with EQUATION Gl and G2. Lp may
represent the longitudinal stiffness of the pinion and Lw may represent the
longitudinal
stiffness of the wheel. It should be noted that these elements may not be
included in an
embodiment using the two degrees of freedom model.
Bearings may also be modeled to obtain estimates of fault conditions in
instances
where there is little or no empirical data available. With bearings, a
periodic impulse is of
interest. The impulse is the result of a bearing rolling over a pit or spall
on the inner or outer
to bearing race. The intensity of the impulse on the bearing surface is a
function of the angle
relative to the fault, which may be represented as, for example, described in
the Stribeck
equation in a book by T.A. Harris, 1966, Rolling Bearing Analysis. New York:
John Wiley p
148 as:
q(9)-q°[l (12E~1-cosA)~ EQUATIONBl
where n = 3/2 for ball bearings and 10/9 for rolling elements bearing, E < .5,
and 8 is less than
~l2 in accordance with values specified in this particular text for the
different bearings used
in the above-referenced Stribeck equation represented as in EQUATION B 1.
An impulse in a solid surface has an exponential decay constant, which may be
taken
2o into account, along with a periodic system due to rotation of the shaft.
The bearing model
may then be represented as a quantity, "s", which is the multiplication of the
impulse, "imp"
below, the impulse intensity, "q(8)" as may be determined above, the period
shaft rotation,
which is "cos(6)" below, all convoluted by the exponential decay of the
material and
represented as:

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
s = imp x q(9) x cos 8~ ~ exp(T l t) EQUATION B2
where T is the exponential decay and t is the time. It should be noted that
"T" varies with the
material of the solid surface. "exp(T/t)" may be obtained, for example, using
a modal
hammer, to generate the decay response experimentally. An embodiment may also
obtain this
value using other information as may be supplied in accordance with
manufacturer's
information. The value of "t" may be a vector of times starting with the first
time sample and
extending to the end of the simulation. T is generally small, so the
expression "exp(T/t)"
approaches zero rapidly even using a high sampling rate.
"imp" is the impulse train that may be represented as the shaft rate =~
bearing frequncy
ratio * sampling rate for the simulation period.
"s" is the simulated signal that may be used in determining a spectrum, "S",
where "S
= fft(s)", the Fourier transform of s into the frequency domain from the time
domain. As
described in more detail in following paragraphs, in determining a CI in
connection with the
bearing model signal "s" having spectrum "S", for example, the Power Spectral
Density of S
at a bearing passing frequency may be used as a CI. Additionally, for example,
other CI
2o values may be obtained, such as in connection with the CI algorithm
comparing the spectrum
"S" to those associated with transmission error in connection with a normal
distribution using
the PDF/ CDF CI algorithms that may be generally described as hypothesis
testing techniques
providing a measure of difference with regard whether the spectmm is normally
distributed.
71

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
It should be noted that, as described elsewhere herein in connection with gear
models,
values may be used in the foregoing equations in connection with simulating
various fault
conditions and severity levels. The particular values may be determined in
accordance with
what small amount of observed data or manufacturer's data may be available.
For example,
in accordance with observed values, an impulse value of 0.02 for the impulse,
"imp", may
correspond to a fairly severe fault condition. Values ranging from 0.001 to
0.03, for example,
may be used to delimit the range of "imp" values used in simulations.
Following is an example of estimated data using the foregoing equations for a
bearing
1o having the following configurations:
Rpm = 287.1
Roller diameter =.25
i5 Pitch diameter = 1.4171
Contact angle = 0
Number of elements =10
Inner race fault
20 Referring now to Figure 31, shown is an example of a graphical
representation of the
signal for the foregoing configuration when there is some type of bearing
fault as estimated
using the foregoing equations EQUATION B 1 and B2. Figure 32 represents the
estimated
spectrum "S" as may be determined using EQUATION B2 above.
72

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
It should be noted that for bearings, there may be three types of faults, for
example,
estimated using the foregoing equations. There may be an inner race fault, an
outer race fault
or a roller element fault. Localized bearing faults induce an excitation which
can be modeled
as an impulse train, expressed as imp in the above equation. This impulse
"imp" corresponds
to the passing of the rolling elements of the fault. Assuming a constant,inner
ring rotation
speed, the impulse train is periodic and the periodicity depends on the fault
location.
For outer race faults, the bearing frequency ratio, f a, or may be represented
as:
N db
1o fd,or = 2 1- dm cos(a) ( f,.r - for ) EQUATION B3
where:
"db" represents the roller diameter,
"dm" represents the pitch diameter,
"a" =2~ * frequency, f~;
" 1,. " is the rotation frequency of the inner race (e.g. shaft rate), and
'for " the rotation frequency of the outer race (if fixed = 0).
2o For inner race faults, the bearing frequency ratio, f d,;r may be
represented as:
fd,ir = N 1 + ~b cos(a) ( f;, - for ) EQUATION B4
2 d",
Replacing a, with 2~f~, the time response is f(t). This substitution may be
performed
as the initial value of a is based on an angle and not a function of time. In
a simulation, there
is a time dependent response as expressed using f(t).
73

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
The radial load applied to the bearing is not constant and results in a load
distribution,
which is a function of angular position. If the defect is on the outer race,
the amplitude of the
impulse is constant because the fault location is not time varying. For an
inner race fault, the
amplitude with respect to angular position. The function is:
t
q(8) = qo (1- 12E (1- cos B)~ foj. lel < emax EQUATION BS
0 elsewhere
q(t) = q(2TCf(6)). This quantity q(t) is amplitude at a particular time, or
q(theta) representing
the amplitude at a particular angle. Amplitude modulation takes into account
the distance
to from the fault to the sensor. For outer race fault, the quantity cos (8) is
constant (1), for inner
race fault, it is the cosine function, noted as "cos (8)" in the above
equation.
For a linear system, the vibrations at a given frequency may be specified by
the
amplitude and phase of the response and the time constant of the exponential
decay. As the
angle, 0 above, changes, the impulse response, h(t), and the transfer function
H(f) also
change due to the changing transmission path and angle of the applied impulse.
It is assumed
that the exponential decays is independent of the angle 8, so that the
response measured at a
transducer due to an impluse applied to the bearing at the location 6 is
characterized by an
amplitude which is a function of 8.
The impulse response function h(t) and the transfer function H(f) may be
replaced by
a function a(8) giving the amplitude and sign of the transfer function H(fJ
at each angle theta and by the exponential decay of a unit impulse, (e(t)).
For an inner race
defect, rotating at the shaft frequency fs, the instantanous amplitude of the
transfer function
between the defect and the transducer as a fimction of time, a(t) may be
obtained by
74

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
substituting 2~*fs~t for theta. Note that a(t) is periodic. At 8 = 0 relative
to the defect and
transducer, a(t) has its maximum value. At 8 = ~c, a(t) should be a minimmn
because the
distance form the defect to the transducer is a minimum. Additionally, the
sign is negative
because the impulse is in the opposite direction. Because of these properties,
the cost) may
be used for the function a(t).
The impulse train is exponentially decaying. The decay of a unit impulse can
be
defined by:
e(t) = exp(-~ ) EQUATION B6
a
for t> 0, where Te is the time constant of decay.
The bearing fault model is then:
v(t) _ ~i~rap(t)q(t)a(t)~ ~ e(t) EQUATION B7
where:
imp(t), which is the impulse over a time t, = 2~ * shaft rate * time * bearing
frequency ratio, as may be determined using EQUATIONS B3 and B4 above;
a(t) is the cos(6) for an inner race, which is 1 for an outer race, where
cos(6) = 0,
where 8 is time varying;
and q(t) and e(t) are as described above.
An embodiment may include a signal associated at the sensor for gear and
bearing noise
combined from the bearing and the gear model may be represented as:

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
s(t) _ ~d (t) f (t)q(t)a(t)~ * e(t) * h(t) EQUATION B8
where:
h(t) is the frequency response of the gear case, as may be determined, for
example,
using an estimate produced with linear predictive coding (LPC) techniques or
with a modal
hammer analysis;
d(t) is the signal associated with gear/shaft T.E. as may be determined using
the gear
model EQUATION Gl;
and other variables are as described elsewhere herein.
l0 The frequency spectrum of signals representing a combined bearing and gear
model
from EQUATION B8 may be represented as:
S(f)=~1~(f)*F(f)*Q(f)*A(f)~E(f)H(f) EQUATION B9
As described elsewhere herein, healthy data, such as may be obtained using
manufacturer's information, may be used in determining different values, such
as those in
connection with stiffnesses for gear simulation, amplitude and exponential
decay for bearing
faults. In terms of generating fault data, since these systems are linear, the
following may be
defined:
~ For gear faults indicative of a crack, a reduction in the stiffness for a
tooth (e.g. 50 and 20
percent of normal) may be used in estimating median and high fault values.
Additionally,
these values may be varied, for example, using the Monte Carlo simulation to
quantify
variance.
~ For shaft misalignment, shaft alignment within the model may be varied to
estimate
mean fault values
76

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
~ Fox gear spalling faults, the "size" of an impulse may be determined through
trial and
error, and by comparing simulation values with any limited observed fault data
previously
collected.
~ For bearing fault models, which are spalling faults, the size of an impulse,
indicative of a
fault, with lcnown bearing faults, may be determined similarly as with gear
spalling faults
Sensitivity analysis may be performed, for example using range of different
input
values for the different parameters, to provide for increasing the
effectiveness of fault
detection techniques, for example, as described and used herein. For example,
an
to embodiment may be better able to simulate a family of bearing faults to
tailor a particular CI
algorithm to be sensitive to that particular fault.
Using the foregoing, the modulated transmission error of a gear mesh, for
example,
which is a signal may be simulated or estimated. This signal may subsequently
be processed
15 using any one or more of a variety of CI algorithms such.that estimates for
the mean and
threshold values can then be derived for fault conditions. (It is assumed that
the stiffness and
torque are known apt°io~i). Parameter values used in the above
equations corresponding to a
healthy gear, for example, as may be specified using manufacturer's data, may
be modified to
estimate parameter values in connection with different types of faults being
simulated. By
2o modifying these parameter values, different output values may be determined
corresponding
to different fault conditions.
For example, known values for stiffness, masses, and the lilce used in
EQUATION G1
may be varied. A cracked gear tooth may be simulated by making the stiffness
time varying.
25 The contact pitch may be varied with time in simulating a shaft alignment
fault. A modulated
77

CA 02439734 2003-08-28
WO 02/095633 PCT/US02/16380
input pulse on d may be used in simulating a spall on a gear tooth. Different
parameter
values may be used in connection with specifying different degrees of fault
severity, such as
alarm levels and warning levels. A particular parameter value, such as a tooth
stiffness of
70% of the normal manufacturer's specified stiffness, may be used in
simulating warning
levels. A value of 20% of the normal manufacturer's specified stiffness may be
used in
simulating alarm levels. The particular values may be determined in accordance
with
comparing calculated values with the characteristics of real CI data on any
few real faults
collected.
l0
While the invention has been disclosed in connection with the preferred
embodiments
shown and described in detail, various modifications and improvements thereon
will become
readily apparent to those skilled in the art. Accordingly, the spirit and
scope of the present
invention is to be limited only by the following claims.
78

Representative Drawing