Methods for improved forewarning of critical events across multiple data channels

[0001] This is a continuation-in-part of Hively, U.S. patent application Ser. No. 10/195,626 filed Jul. 12, 2002 now U.S. Pat. No. 7,139,677.

[0002] This invention was made with assistance under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The Government has certain rights in this invention.

[0003] The field of the invention is computer methods for analyzing nonlinear processes to forewarn of critical events in nonlinear processes

[0004] Examples of critical events are mechanical or electrical failures in machines; epileptic seizures, ventricular fibrillations, fainting (syncope), breathing difficulty, and sepsis in human medical patients; and other physical processes. Further examples of nonlinear processes include brain waves, heart waves, chest sounds, transients in power systems, airflow over automobiles and airplanes, weather and climate dynamics, water flow around submarines, machine tool-part interaction (e.g., tool chatter), nuclear reactor instabilities, fusion plasma instabilities, earthquakes, turbulent flow in pipes, planetary/satellite motion.

[0005] Engineering, medical, and research applications frequently must distinguish a difference between two apparently similar but actually different states in a nonlinear process. Examples include various data for machine failure, pre-seizure versus non-seizure brain waves, pre-fibrillation versus fibrillation heart waves, pre-syncope versus syncope heart waves, pre-sepsis versus sepsis heart waves, and normal versus abnormal chest sounds as an indicator of breathing difficulty. The electrical/mechanical community calls this problem, “condition monitoring.” The equivalent term in the medical community is diagnostic, medical, or health monitoring. In the computer/networking world, the term is “security monitoring.”

[0006] Hively et al., U.S. Pat. Nos. 5,743,860 and 5,857,978 disclose methods for detecting and predicting epileptic seizures by acquiring brain wave data from a patient, and analyzing the data with traditional nonlinear methods.

[0007] Hively et al., U.S. Pat. No. 5,815,413 discloses the applicability of nonlinear techniques to monitor machine conditions, such as the condition of a drill bit or the performance of an electrical motor driving a pump.

[0008] Clapp et al., U.S. Pat. No. 5,626,145 discloses the removal of artifacts representing eye blinks from EEG data using a zero-phase quadratic filter. As one normally skilled in the art can appreciate, the same method can remove other artifacts from data, such as quasi-periodic variations from three-phase electrical current, voltage, or power; and quasi-periodic oscillations from one or more channels of acceleration; and breathing artifacts from heart waves and chest sounds;.

[0009] Hively et al., U.S. Pat. No. 6,484,132 discloses a plurality of nonlinear techniques for using and enhancing phase space dissimilarity measures (PSDM) to forewarn of machine failures, as well as impending epileptic events from scalp EEG in ambulatory settings. PSDM yield superior performance over traditional nonlinear indicators, such as Kolmogorov entropy, Lyapunov exponents, and correlation dimension.

[0010] The present invention improves on prior disclosed methods involving the use of phase space dissimilarity measures (PSDM) to provide forewarning indications of critical events. Such events include failures in critical mechanical or electrical equipment, in medical treatment of human patients, and in many other applications.

[0011] The invention provides for better selection of process-indicative data for such analysis, such that one data channel can be used in place of multiple data channels without sacrificing consistent forewarning of critical events. Examples of such a process-indicative data are vibration power for mechanical systems, electrical power for electrically-driven machines, and differences between adjacent scalp EEG (electroencephalogram) channels.

[0012] The invention also provides a novel method to check the quality of the data, such that batches of inadequate quality data can be identified and replaced by adequate quality data.

[0013] The invention also provides a novel method including at least two separate filtering operations to remove both high-frequency and low-frequency artifacts using a zero-phase quadratic filter, which previously has been used only removing low-frequency artifacts in the prior art.

[0014] The invention also gives a novel method to combine several channels of time-serial data into a dC-dimensional multi-channel phase-space (PS) vector, which has the form: s_i=[s(1)_i, s(1)_i+λ, . . . , s(1)_i+(d−1)λ, . . . , s(C)_i, s(C)_iλ, . . . s(C)_i+(d-1)λ]. Here, s_i(j) denotes the symbolized data for j-th channel at time t_i. C denotes the total number of data channels in the multi-channel phase-space reconstruction.

[0015] The invention also provides an end-of-life forewarning indicator (G), which is computed from a composite measure of dissimilarity (C_i). Here C_i is the sum of the four renormalized PSDMs for the i-th subset:
C_i=U(X²)+U(X_c²)+U(L)+U(L_c)  (1)
The end-of-life indication (G) is quantified from this composite measure as follows. The appropriate process-indicative data is selected and acquired (for example) over 10-second snapshots at periodic intervals (such as once per minute), and then divided into subsets of N data points, which typically has a value in the range, 10,000≦N≦100,000. The analysis then proceeds as follows: (a) each subset is artifact filtered; (b) the data in each subset is symbolized; (c) each subset is converted into a sequence of multi-channel phase-space vectors; (d) the phase-space vectors are tabulated into non-connected and connected phase-space distribution functions; (e) the distribution functions for each sequential test state (subset) are compared to the nominal state distribution functions via the four measures of dissimilarity in Eqs. (7)–(10); (f) each measure of dissimilarity is renormalized; (g) these four renormalized measures of dissimilarity are used to construct sequential values for the composite measure of condition change C_i via Eq. (1); and (h) a standard least-squares method is used to fit several sequential values of C_i to a straight line:
y_i=ai+b  (2)
We choose a window for this straight-line fit with a length of m=10 values of C_i (and y_i below), that is consistent with the number of data subsets in each data snapshot. Other values of m give inferior indication of condition change. Next, the variance, σ₁², measures the variability of the C_i values about this straight-line fit:
σ₁²=Σ_i(y_i−C_i)²/(m−1)  (3)
Finally, the end-of-life statistical indicator, G, measures the variability of next m values of C_i about an extrapolation of the straight-line fit as follows:
G=Σ_i(y_i−C_i)²/σ₁²  (4)
The index, i, runs over the m values of C_i and y_i in Eqs. (2)–(4). G has the form of a conventional chi-squared statistic, but that notation is not used to avoid confusion with the two chi-squared PSDMs, U(X²) and U(X_c²). Moreover, this disclosure does not use the G measure for a statistical test of null hypothesis. Rather in this invention, G is a relative measure, providing end-of-life forewarning, plus an indication of the failure onset.

[0016] Other objects and advantages of the invention, besides those discussed above, will be apparent to those of ordinary skill in the art from the description of the preferred embodiments, which follows. In the description reference is made to the accompanying drawings, which form a part hereof, and which illustrate examples of the invention. These examples, however are not exhaustive for the various embodiments of the invention, and therefore reference is made to the claims, which follow the description for determining the scope of the invention.

[0026] The following description incorporates the methods disclosed in U.S. patent application Ser. No. 10/195,626 filed Jul. 12, 2002. To the extent that the methods of the present invention build upon the methods disclosed there, the disclosure of that application is hereby incorporated by reference.

[0027] The methods of this disclosure can be applied to electric motor predictive maintenance, other machinery and physical processes, as well as biomedical data for diagnosis and treatment of human patients. In one example, data sets were recorded in snapshots of 1.5 seconds, sampled at 40 kHz (60,000 total time-serial samples), including three-phase voltages and currents from an electric motor. The subsequent description describes analysis of an induced fault, referred to as a seeded fault.

[0028] In such a motor, the components of the three-phase voltages (V_i) and currents (I_i), (see FIGS. 1A–1F) can be converted into an instantaneous motor power, P=Σ_iI_iV_i, where the sum runs over the three phases. The electrical power, P, (FIG. 1G) is advantageous in forewarning of critical events, because it reduces the number of data channels that must be analyzed to detect the condition change. The PSDM were computed for each of the three-phase currents and voltages separately, and compared with the condition change for components of the three-phase power. FIGS. 2A–2L show examples of the four dissimilarity measures for each component (I_iV_i) of three-phase power. All four dissimilarity measures of I₂V₂ (middle column) show an almost linear increase, and all four dissimilarity measures of I₃V₃ increase monotonically. All four dissimilarity measures of I₁V₁ increase from test #1 (no fault) to test #2 (first seeded fault), then remain constant from test #2 (first seeded fault) to test #3 (second seeded fault in addition to the first seeded fault). Thus, it was determined that the phase space dissimilarity measures (PSDM) method described above can extract condition change from one (or more) components of three-phase electrical power as well as from individual phase currents and voltages. The advantage of selecting electrical power as the process-indicative signal is the reduction in the number of channels of data for analysis.

[0029] The data for this same test also included tri-axial acceleration data from inboard (IB) and outboard (OB) motor locations. Vibration (mechanical) power can be determined from tri-axial acceleration, as follows. Acceleration, a, is a three-dimensional vector that can be integrated once in time to give velocity (vector), v=∫adt. Mass, m, times acceleration (vector) is force (vector), F=ma. The dot-product of force and velocity converts these vector quantities into the mechanical power (scalar), P=F·v=ma·∫a dt.

[0030] FIGS. 3A–3G show the three components of inboard acceleration (FIGS. 3A–3C), as well as the corresponding three components of inboard velocity (FIGS. 3D–3F). The noteworthy feature of these plots is the clear periodicity with additional complex nonlinear features, which also appear in the resulting vibration power (FIG. 3G).

[0031] FIGS. 4A–4D show a nearly linear rise in all four PSDM as computed from the vibration power of the previous paragraph. To obtain this result, a sequence of single parameter searches was executed: varying S with λ fixed, then varying λ with S set to the value from the first search. This approach, yielded S=85 and λ=42 for d=3. Other dimensions (d=2 and 4) give inferior results. Subsequently, an exhaustive search was made over a range of S and λ for d=3 to obtain best condition change at S=84 and λ=47.

[0032] The above examples show that electrical power has renormalized values of PSDM in the range of 5–55 standard deviations from the mean; two of the PSDM rise monotonically (not linearly), thus yielding less desirable change discrimination. It was found that vibration power has PSDM values on the order of unity; all four PSDM increase linearly, thus yielding good change discrimination. A search to find one (or more) data channels for consistent event forewarning requires the analysis of all three components of acceleration, or all three components of velocity, or all three components of three-phase electrical power, or all six components of three-phase voltage and current. However, only one channel of electrical power or one of vibration (mechanical) power needs to be analyzed to obtain consistent failure forewarning. As one normally skilled in the art can readily appreciate, the savings in computation effort is an important and novel improvement to the art.

[0033] Another example uses multi-channel EEG data that is combined into a multi-channel PS reconstruction for critical-event forewarning. FIG. 5 shows an array of monopolar scalp EEG electrodes, each of which measures a local voltage with respect to some reference electrode. The voltage difference between adjacent monopolar electrodes creates a “bipolar montage.” One bipolar channel is (for example) the difference, F4−C4, on the right hemisphere; another is F3−C3 on the left hemisphere. Each bipolar channel extracts a time-serial signal that potentially allows removal of most locally-generated artifacts, while retaining a differential measure of the brain activity. Our analysis of EEG data has determined that data from the bipolar montage without artifact filtering provides a reasonable alternative to artifact filtering, yielding better forewarning than the monopolar montage with (or without) artifact filtering (and not much different from bipolar forewarning with artifact filtering). Moreover, two (or more) of these bipolar channels can be combined via the multi-channel phase-space vector, to provide forewarning of epileptic seizures that is comparable to the analysis of all nineteen of the individual monopolar channels. As one normally skilled in the art can readily appreciate, the savings in computational effort (one two-channel bipolar PS reconstruction vs. nineteen individual PS reconstructions) is an important and novel improvement.

[0034] After choosing the appropriate process-indicative data in a manner that reduces the analysis effort, as discussed above, the present invention next provides an improved method to identify inadequate quality data. Raw process-indicative data can come from a variety of sensors, such as electroencephalogram data, electrocardiogram data, chest sounds, single- or multi-axial accelerometer(s), single- or multi-phase motor current(s) and voltage(s). This data is sampled at equal time intervals, τ, starting at time t₀, yielding a sequence of time serial data, x_i=x(t₀+iτ), i=1, 2, . . . , N. These signals are digitized and recorded as either integer or decimal numbers, which have a range between from some minimum value (x_min) to maximum value (x_max).

[0035] Determination of the data quality from the time serial data, x_i, involves the following steps: (a) sorting the x_i values into ascending order from x_min to x_max; (b) determining the number of unique signal values (n); (c) obtaining the occurrence frequency (F_k) for each unique signal value (v_k); (d) estimating the occurrence probability for a single outlier in a normal distribution by solving: 1/n=½ erfc(z/√2). Here, n is the number of unique signal values from step (b); erfc is the complementary error function; z is the number of standard deviations from the mean for such a single outlier; and the mean occurrence frequency is estimated as N/n. Next, (e) plotting F_k versus v_k in each BIN, as shown (for example) in FIGS. 6A–6C; (f) reject any outlier bin with [v_k−(N/n)]/σ₃>z, where σ₃ is the standard deviation of the occurrence frequency over all of the unique signal values from step (b). In this example, FIGS. 6A–6C show frequency distributions from real test data. Huge spikes in log₁₀FREQUENCY occur above the underlying baseline (top two subplots for channel 1 (C1) and channel 2 (C2), indicating an accumulation of data values in one (or a few) unique signal value(s). These spikes far exceed the above limit. Another problem is that the c1- and c2-signals have a narrow range of values (−0.099 to +0.089), indicating that the signal is improperly scaled for digitization. Moreover, the c1- and c2-distribution functions have large variability, which further indicates excessive noise.

[0036] Table 1 summarizes the results for samples of tri-axial accelerometer data that was acquired from a test motor. The data were sampled for 10 seconds at 128 kHz for a total of 1,280,000 points from a shaft crack test. The table below shows the minimum (min) and maximum (max) values by dataset for each of three accelerometer channels (c1–c3).

[0037]

TABLE 1
Tri-axial Acceleration Data
DATASET	min (c1)	max (c1)	min (c2)	max (c2)	min (c3)	max (c3)	#BINS (c3)
FSF40c3	−.08755	+.07760	−.08748	+.06009	−.00323	−.00254	11
FSF20c5	−.08899	+.07375	−.08247	+.05514	−.00322	−.00261	10
FSF20c4	−.09895	+.07684	−.08893	+.05906	−.00316	−.00261	9
UNCRc1	−.09257	+.08611	−.07705	+.05892	−.00316	0.0	11
UNCRc3	−.08549	+.07725	−.08364	+.05572	−.00309	−.00247	10

These data have discrete values with five significant figures (x.xxxx) for all three channels. The right-most column of Table 1 (with the header “#BINS”) shows the number of different discrete c3-values, which varies between nine and eleven over a very narrow range of signal values (−0.0032 to +0.0026). This small number of discrete signal values is also seen in FIG. 6C, corresponding to a data precision of log₂(11)=3.45. This low data precision is well below that attainable by modern data acquisition systems, which can acquire data with >10 bits of precision, and more typically 16 (or more) bits. Thus, the c3-signal is due to bit chatter (noise), rather than a real signal.

[0038] Rejection of this inadequate quality data from further analysis avoids meaningless or nonsensical results. In the prior art represented by U.S. Pat. No. 5,815,413, the data range was partitioned into an equal number of non-overlapping uniform-size bins to obtain the occurrence frequency. The improvement of this disclosure is the use of the unique signal values in the raw data. As one of ordinary skill in the art can appreciate, this novel improvement uses the explicit features in the data to determine its quality, rather than artificial constructs of the prior art (such as analyst-chosen bins to tabulate the occurrence distribution function).

[0039] Data of the types discussed above, may also have artifacts that would otherwise confound the analysis, and therefore must be removed. In Clapp et al., U.S. Pat. No. 5,626,145, “Method and Apparatus for Extraction of Low-Frequency Artifacts from Brain Waves for Alertness Detection,” a method to remove low-frequency artifacts (such as eye-blinks in EEG data or breathing artifacts from ECG data) via a zero-phase quadratic filter is disclosed. This filter uses a moving window of process-indicative data points, e_i, with the same number of data points, w, on either side of a central point. A quadratic curve is fitted in a least-squares sense to these 2w+1 data points. The central point of this fit is taken as the best estimate of the low-frequency artifact, f_i. The residual signal, g_i=e_i−f_i, has essentially no low-frequency artifact activity.

[0040] However, data may have both low-frequency artifacts (such as eye-blinks in EEG, or breathing in electrocardiogram, ECG, data), and high-frequency artifacts (such as 60-Hz noise in ECG). The novel improvement taught in this disclosure is application of the original zero-phase quadratic filter first to extract the low-frequency f_i-data of the previous paragraph using an appropriate (small) value for the filtering-window width, w, which effectively removes the high-frequency artifacts. Subsequently, the same filter is applied to this f_i-data with an appropriate (large) filtering window width for w, which removes the low-frequency artifacts from the resultant second-stage g_i-data of the previous paragraph. The use of the zero-phase quadratic filter in two stages as applied to PSDM data analysis has not been seen in the prior art and improves forewarning analysis, for example in ECG data with both low-frequency breathing artifacts and high-frequency 60-Hz noise. As one normally skilled in the art can appreciate, the removal of both low- and high-frequency artifacts while preserving the nonlinear data features is an important and novel component of the analysis.

[0041] Next, each artifact-filtered value, g_i, is converted into a symbolized form, s_i, that is, one of S different integers, 0,1, . . . , s−1:
0≦s_i=INT[S(g_i−g_min)/(g_max−g_min)]≦S−1  (5)
Here, the function INT converts a decimal number to the closest lower integer; g_min and g_max denote the minimum and maximum values of artifact-filtered data, g_i, respectively, over the baseline (nominal state) data. To maintain S distinct symbols, the following expression holds, namely s_i=S−1 when g_i=g_max. Expression (5) creates symbols that are uniformly spaced between the minimum and maximum in signal amplitude (uniform symbols). Alternatively, one can use equiprobable symbols, by ordering all N base case data points from the smallest to the largest value. The first N/S of these ordered data values correspond to the first symbol, 0. Ordered data values (N/S)+1 through 2N/S correspond to the second symbol, 1, and so on up to the last symbol, S−1. By definition, equiprobable symbols have non-uniform partitions in signal amplitude and present the advantage that dynamical structure arises only from the phase-space reconstruction. Large negative or positive values of g_i have little effect on equiprobable symbolization, but significantly change the partitions for uniform symbols. Moreover, the mutual information function is a smooth function of the reconstruction parameters for equiprobable symbols, but is a noisy function of these same parameters for uniform symbols. Thus, equiprobable symbols provide better indication of condition change than uniform symbols when constructing a (non)connected phase space.

[0042] Phase-space (PS) construction uses time-delay vectors, s_i=[s_i, s_i+λ, . . . , s_{i+(d-1) λ}] to unfold the underlying dynamics. Critical parameters in this approach are the time delay, λ, and system dimensionality, d, and the type and number of symbols, S. However, information exchange between (for example) complex machinery components or between different brain locations connect the local processes. This disclosure provides a novel method to combine several channels of time-serial data into a multi-channel phase-space (PS) reconstruction, which has the form of a dC-dimensional vector: s_i=[s(1)_i, s(1)_i+λ, . . . , s(1)_i+(d-1)λ, . . . , s(C)_i, s(C)_i+λ, . . . s(C)_{i+(d-1)λ]. Here, si}(j) denotes the symbolized data for j-th channel at time t_i. C denotes the total number of data channels in the multi-channel phase-space reconstruction.

[0043] Conversion of the multi-channel time serial data to a sequence of phase-space vectors via the method of the previous paragraph divides the phase space into S^Cd discrete bins, which are uniquely labeled via base-S arithmetic:

[0044] $\begin{array}{cc}0\le y\left(i\right)=y\left({\underset{\_}{s}}_i\right)=\sum _{j=1}^c\sum _{k=1}^{d-1}{s\left(j\right)}_{i+k\lambda }{S}^{k+d\left(j-1\right)}\le S^{\mathrm{Cd}}-1& \left(6\right)\end{array}$
The resulting distribution function (DF) is a discretized density, which is obtained by counting the number of points that occur in each phase space bin. The population of the ith bin of the distribution function, is denoted Q_i, for the base case, and R_i for a test case, respectively.

[0045] The test case is compared to the base case via dissimilarity measures between Q_i with R_i that are defined by:

[0046] $\begin{array}{cc}{\chi }^2=\sum _i^{}{\left(Q_i-R_i\right)}^2/\left(Q_i+R_i\right),& \left(7\right)\\ L=\sum _i^{}Q_i-R_i,& \left(8\right)\end{array}$
Here, the summations run over all of the populated phase space cells. These measures account for changes in the geometry and visitation frequency of the resultant geometric object, called an “attractor.” The X² measure is one of the most powerful, robust, and widely used statistics for comparison between observed and expected frequencies. In this context, X² is a relative measure of dissimilarity, rather than an unbiased statistic for accepting or rejecting a null statistical hypothesis. The L distance is the natural metric for distribution functions by its direct relation to the total invariant measure on the phase-space attractor. Consistent calculation of these measures obviously requires the same number of points in both the base case and the test case distribution functions, identically sampled; otherwise, the distribution functions must be properly rescaled.

[0047] The process dynamics yield successive PS points, as prescribed by s_i→s_i+l. Thus, a discrete representation of the process flow is obtained in the form of a 2dC-dimensional multi-channel PS vector, S_i=[s_i, s_i+1], that is formed by adjoining two successive vectors from the dC-dimensional reconstructed PS. S_i is a vector for the multi-channel connected phase space (CPS). As before, Q and R denote the CPS distribution functions for the base case and test case, respectively. The measure of dissimilarity between the two distribution functions for the CPS, signified by the “c” subscript are thus defined as follows:

[0048] $\begin{array}{cc}{\chi }_c^2=\sum _{\mathrm{ij}}^{}{\left(Q_{\mathrm{ij}}-R_{\mathrm{ij}}\right)}^2/\left(Q_{\mathrm{ij}}+R_{\mathrm{ij}}\right)& \left(9\right)\\ L_c=\sum _{\mathrm{ij}}^{}Q_{\mathrm{ij}}-R_{\mathrm{ij}}.& \left(10\right)\end{array}$
The first subscript in Eqs. (9)–(10) denotes the initial PS point, and the second subscript denotes the sequel PS point in the CPS pair. These CPS measures have higher discriminating power than unconnected PS measures of dissimilarity. Indeed the measures defined in Eqs. (7)–(10) satisfy the following inequalities: X²≦L, X_c²≦L_c, L≦L_c, and X²≦X_c², where X_c² and L_c are dissimilarity measures for connected phase space and X² and L are dissimilarity measures for non-connected PS.

[0049] To assure robustness, the construction of the base case data requires careful statistics to eliminate outlier base case subsets. For a description of this part of the methodology, reference is made to Hively, U.S. patent application Ser. No. 10/195,626 filed Jul. 12, 2002.

[0050] A consistent method of comparison is needed to interpret the disparate range and variability of the dissimilarity measures. To this end, the dissimilarity measures are renormalized, as described next. The B non-outlier base case subsets from the previous paragraph are compared to each test case subset, to obtain the corresponding average dissimilarity value, V_i, of the ith subset for each dissimilarity measure. Here, V denotes each dissimilarity measure from the set, V={L, L_c, X², and X_c²}. The mean value, V, and the corresponding standard deviation, σ₂, of the dissimilarity measure V are calculated by comparison of all of the unique pairs of base case subsets, after the outliers have been eliminated. The renormalized dissimilarity is the number of standard deviations that the test case deviates from the base case mean: U(V)=|V_i−V|/σ₂.

[0051] Once the renormalized measures for the test and base cases have been obtained, a threshold, U_c, is selected for each renormalized measure U to distinguish between normal (base) and possibly abnormal (test) regimes. The choice of a reasonable threshold is critical for obtaining robust, accurate, and timely results. A forewarning indication is obtained when a renormalized measure of dissimilarity exceeds the threshold, U≧U_c, for a specified number, N_OCC, of sequential occurrences.

[0052] The remainder of the PSDM analysis is described by Hively, U.S. patent application Ser. No. 10/195,626 filed Jul. 12, 2002, cited above. As one normally skilled in the art can readily appreciate, the novel improvement to the prior art involves combining multiple-channel time-serial data into a one multi-channel time-delay phase-space (PS) reconstruction, rather than performing the PS reconstruction for each of several channels.

[0053] An alternative embodiment of this invention is illustrated by the following example. The vibration power data, as discussed above, was acquired in 10-second snapshots at 50 kHz at fifteen-minute intervals; see Table 2 for the interval (Δτ) between each snapshot. Each 500,000-point snapshot was divided into ten non-overlapping subsets of 50,000-points each for subsequent analysis. This data characterized an accelerated-failure test of an over-loaded gearbox. The test protocol involves a break-in period at the nominal (1×) load (530 ft-lbs) for one hour, followed by twice (2×) or three times (3×) the normal load, as shown in Table 2 The end-of-life failure included pinion gear damage, broken teeth, and a sheared shaft. The failure occurred after 162.5 hours, corresponding to 650 snapshots. FIGS. 7A–7L show the very complex, nonlinear features in the three components of the tri-axial accelerometer data (a). The individual snapshots were combined into one long dataset (12.7 GB). Acceleration was subsequently converted to a long stream (4.1 GB) of vibration power (FIGS. 8A–8D), using the previously described method to convert tri-axial acceleration into power via time-integration to velocity (v=∫adt) with a subsequent vector dot-product to produce power (P˜a·v). FIGS. 9A–9D show the PSDM from this data. The phase-space parameters are N=50,000, S=274, d=2, λ=1. All four dissimilarity measures rise systematically (FIGS. 9A–9D) to provide forewarning of the failure. However, a more robust and specific end-of-life (EOL) indicator is needed. We observe that all four of the PSDMs have similar tends, suggesting the definition of a composite measure, C_i, as the sum of the four renormalized PSDMs for the i-th dataset (FIG. 9E):
C_i=U(X²)+U(x_c²)+U(L)+U(L_c)  (11)
This composite measure is expected to be more robust than any one of the individual measures of phase space dissimilarity, while accurately indicating condition change. Subsequent description presents another novel improvement of this disclosure to provide end-of-life forewarning, and indication of failure onset.

[0054] The end-of-life indication from the composite measure of the previous paragraph is quantified as follows. We use contiguous, non-overlapping windows of C_i to obtain the best straight-line fit in a least-squares sense:
y_i=ai+b  (12)
The window length of m=10 values of C_i (and y_i below) is chosen consistent with the number of subsets in each snapshot. Other values of m give inferior indication of condition change. Next, the variance, σ₁², measures the variability of the C_i values about this straight-line fit:
σ₁²=Σ_i(y_i−C_i)²/(m−1)  (13)

[0055] Other fits (such as, quadratic, cubic, and quartic) give inferior condition change indication, due to the poor extrapolation outside the fitting window. Finally, G measures the variability of next m values of C_i about an extrapolation of this straight-line:
G=Σ_i(y_i−C_i)²/σ₁²  (14)
The index, i, in Eqs. (12)–(14) runs over the m values of C_i and y_i. Note that G has the form of a conventional chi-squared statistic, but we do not use that notation to avoid confusion with the two X² PSDMs, U(X²) and U(X_c²). A statistical test for G would involve (for example) the null hypothesis that deviations from the straight-line fit are normally distributed. Standard chi-squared statistical tables give the corresponding value of G≦28.5 for m=10 degrees of freedom with a probability of one out of the number (650) of snapshots (1/650˜1.5×10⁻³). However, we observe many instances with G>28.5 throughout this test sequence. These outliers occur because the underlying tri-axial acceleration has dynamical correlations, thus violating the requirement for independent, identically distributed samples. Instead, we use G (solid curve in FIG. 9F) as a relative measure of end-of-life. As time moves forward, we obtain the running maximum of G, as the larger of G and G_max (dashed curve in FIG. 9F). We neglect the first six G-values (for example) to avoid startup transients. This running maximum rises in modest increments to 376 over the first 159.75 hours of the test, while intermediate values of G fall well below G_max. The chain curve (-.-) in FIG. 9F is the ratio, R=(G_max)_k/(G_max)_k-1, of the current running maximum in G, (G_max)_k, to the previous maximum in G, (G_max)_k-1. G_max rises to 2,493 at 160 hours, with a corresponding R=6.62, which is well above the largest non-end-of-life value, R=2.22, at 2 hours. Thus, G provides clear end-of-life forewarning, plus indication of the failure (G=244,655). Table 2 gives: (a) the largest non-EOL value of R (R_NEOL) and the corresponding value of G (G_NEOL); (b) values of R (R_EOL) and G (G_EOL) that indicate the end of life, and the matching time (T_EOL/T_FAIL); (c) the value of G at failure onset (G_ONSET) and the corresponding time (T_ONSET/T_FAIL); and (d) the failure-endpoint time (T_FAIL).

[0056]

TABLE 2
Summary of Gearbox Overload Results
	Over-	Δτ					TEOL		TONSET	TFAIL
Run	load	min.	RNEOL	GNEOL	REOL	GEOL	TFAIL	GONSET	TFAIL	hr.
36	2X	15	2.22	376	6.62	2,493	.985	244,655	0.998	162.50
37	3X	1	1.79	333	8.07	2,690	.956	16,284	0.996	8.55
38	3X	1	6.20	374	11.71	13,486	.938	48,379	0.990	4.02
39	2X	1	2.32	853	3.89	5,231	.980	5,231	0.980	8.60
39	3X	1	2.88	1,151	29.03	33,415	.972	44,552	0.994	8.60

Runs #36–38 have largest non-EOL values: R_NEOL=6.20 and G_NEOL=376. The smallest EOL values are: R_EOL=6.62 and G_EOL=2,493. Thus, limits (for example) of R>6.4 and G>1,800 provide EOL forewarning. Moreover, we find that the largest EOL value of G_EOL=13,486, while the smallest failure-onset value is G_ONSET=16,284. Thus, an intermediate value (for example) of G>15,000 distinguishes the EOL from failure onset forewarning This approach gives quantitative limits for transitions from nominal operation (green-light in a traffic signal metaphor), to forewarning of failure (yellow light for caution), and finally to failure onset (red-light for stop).

[0057] Run #39 involves a different test protocol: a one-hour break-in period at nominal load (1×), followed by 2× load for two hours, after which the load alternates between 3× and 2× loads for ten and five minutes, respectively. Run #39 seeks failure forewarning in the presence of load changes. Table 2 shows that the above limits for G and R also distinguish between the non-EOL (green) and EOL (yellow) states for the 3×-portion of this test, because the higher overload drives the failure. These limits do not apply to the 2× test, due to the reduced damage at the lower overload. Unsurprisingly, a different limit of G>38,000 (for example) distinguishes between the EOL and failure onset forewarnings, due to the change in test protocol. The green-yellow-red traffic-signal metaphor still applies for this test.

[0058] As one of ordinary skill in the art can appreciate, these examples give quantitative limits for transitions from normal operation (normal condition), to forewarning of failure (caution condition), and finally to failure onset (critical event), as a novel improvement in the art.

[0059] This has been a description of detailed examples of the invention. These examples illustrate the technical improvements, as taught in the present invention. It will be apparent to those of ordinary skill in the art that certain modifications might be made without departing from the scope of the invention, which is defined by the following claims.