Method for monitoring rotary mechanical system of aircraft

Abstract

The invention relates to a method for monitoring a rotary mechanical system of an aircraft. The method comprises, inter alia, the steps of determining (205) an estimated power spectral density (V) of the noise (v) from a signal (x) representative of vibrations of the rotary mechanical system, determining (207) an equalized power spectral density (Q) by dividing the estimated power spectral density (X) of the signal (x) by the estimated power spectral density (V) of the noise (v), extracting (215) a peak (Q-detection) from the equalized power spectral density (Q), and identifying (217) a fault in the rotating mechanical system from the extracted peaks.

Description

Technical Field

[0001] The present invention relates to the field of monitoring the state of a mechanical system constituting an aircraft turbomachine by analyzing vibration signals. In particular, the present invention relates to a method for monitoring a rotating mechanical system of an aircraft. Background Art

[0002] The prior art includes in particular the documents WO-A1-2010 / 094915, WO-A1-2014 / 184657 and a scientific paper entitled “Bayesian extreme value statistics for novelty detection in gas-turbine engines” published by David A. Clifton et al. at the 2008 IEEE Aerospace Conference.

[0003] In aircraft, rotating machine elements, in addition to normal operation, generate vibrations which may be caused by, for example, imbalances, eccentric masses, worn rolling bearings, bent shafts or misaligned shafts.

[0004] In this case, the vibrations generated are periodic in nature, with the period of the vibrations generated depending on the component and the type of damage suffered by the component.

[0005] For example, an unbalance fault in a rotor produces a pure sinusoidal component whose instantaneous frequency is equal to the rotational frequency and whose amplitude is roughly proportional to the square of the rotational speed.

[0006] Another example is parallel misalignment of the shaft, which results in sinusoidal vibrations with a frequency twice that of the shaft.

[0007] These failures may result in a shortened service life of the associated mechanical component or even in a premature shutdown of the mechanical system into which the mechanical component is integrated. In addition, the failure of one mechanical component may affect another mechanical component belonging to the same mechanical system. For example, misalignment of a shaft causes premature wear of bearings, seals, and couplings in the drive train to which the shaft belongs.

[0008] Spectral analysis is a well-known technique for analyzing vibration signals. The periodicity of the vibration results in a signature of the vibration consisting of a collection of peaks in the spectral domain (or Fourier domain).

[0009] Thus, based on a priori knowledge of the kinematics of the monitored mechanical system, the operator can infer the state of said mechanical system, or even the state of a single mechanical element of the system, by monitoring the evolution of the amplitude of this feature. However, in practice, a priori knowledge of the monitored mechanical system cannot be guaranteed.

[0010] Furthermore, in environments such as aircraft turbines, blind detection of spectral peaks is complicated by the presence of high random noise.Vibration signals from aero-mechanical systems are known to have very low signal-to-noise ratios and to exhibit colored noise (as opposed to white noise).

[0011] In most prior art methods, the noise is usually considered to be white Gaussian noise, but in reality this is not the case.

[0012] Furthermore, for a rotating aircraft mechanical system, the frequency structure of the vibration signal spectrum is very complex, making it impossible to know the number of peaks. Summary of the invention

[0013] The present invention proposes a solution to these disadvantages.

[0014] To this end, according to a first aspect, the invention relates to a method for monitoring a rotating mechanical system of an aircraft, the method comprising the following steps:

[0015] - acquiring a signal representative of vibrations in a rotating mechanical system;

[0016] - determining an estimated power spectral density of the signal;

[0017] - determining an estimated power spectral density of the noise from the estimated power spectral density of the signal;

[0018] The monitoring method is characterized in that the method further comprises the following steps:

[0019] - determining an equalized power spectral density, the equalized power spectral density being equal to the estimated power spectral density of the signal divided by the estimated power spectral density of the noise;

[0020] - determining a cut-off power spectral density, the cut-off power spectral density being equal to the set of values ​​of the equalized power spectral density that are below a first threshold;

[0021] -Adjust the truncated power spectral density by truncating the gamma distribution;

[0022] - determining a second detection threshold based on a truncated gamma distribution;

[0023] - extracting peaks above a detection threshold from the equalized power spectral density; and

[0024] - Identify faults in rotating machinery systems by comparing the extracted peaks with a reference power spectral density.

[0025] The method according to the invention may include one or more of the following features, taken alone or in combination with one another:

[0026] The signal representative of the vibration is a vibration signal or a sound signal, the vibration signal is preferably obtained from an accelerometer or the sound signal is preferably obtained from a microphone.

[0027] -The estimated power spectral density of the signal is determined by calculating the average periodogram of the signal.

[0028] - The estimated power spectral density of the noise is determined by a regression applied to the estimated power spectral density of the signal.

[0029] -The regression is applied to the natural logarithm of the estimated power spectral density of the signal.

[0030] - The regression function used is the B-spline function.

[0031] - The first threshold is set such that a given percentage of the value of the equalized power spectral density is retained in the truncated power spectral density.

[0032] - The conditioning of the truncated power spectral density by a truncated gamma distribution is performed by maximizing a likelihood function of the probability density of said truncated gamma distribution.

[0033] The detection threshold is determined so that the probability density of the truncated gamma distribution is equal to a predetermined value called the false alarm probability.

[0034] The reference power spectral density corresponds to the equilibrium power spectral density obtained for the rotating mechanical system of a fault-free aircraft turbomachine. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] The invention will be better understood and other details, features and advantages of the invention will become more apparent from the following description made by way of non-limiting example and with reference to the accompanying drawings, in which:

[0036] Figure 1 is a schematic diagram of the steps in a monitoring method according to an embodiment of the present invention;

[0037] Figure 2 is a schematic diagram of steps in a monitoring method according to an embodiment of the present invention;

[0038] Figure 3 is an example of a power spectral density of a signal representing vibrations of a rotating mechanical system according to the present invention;

[0039] Figure 4 is an example of an equalized power spectral density of a signal representing vibrations of a rotating mechanical system according to the present invention;

[0040] Figure 5 is based on Figure 4 An example of a histogram obtained by equalizing the power spectral density is shown;

[0041] Figure 6 is based on Figure 4 An example of a truncated histogram obtained by equalizing the power spectral density is shown;

[0042] Figure 7 is an example of equalized power spectral density and detection threshold according to the present invention; and

[0043] Figure 8 is an example of a peak value detected from a signal representing vibrations of a rotating mechanical system according to the present invention.

[0044] Elements having the same function in different exemplary embodiments are provided with the same reference numerals in the figures. DETAILED DESCRIPTION

[0045] Reference Figure 1 and Figure 2 , we will now describe an example of implementation of a method according to the invention for monitoring a rotating mechanical system of an aircraft.

[0046] In the following, a rotating machine system is an assembly of rotating machine elements that have a single function and are spatially grouped together. For example, an aircraft turbine consists of multiple rotating machine systems, including an auxiliary gearbox (AGB). Various machine elements generate vibrations during operation, and the frequency and amplitude of the vibrations are affected by the degradation of the condition of the machine elements.

[0047] In the example shown, step 201 comprises acquiring a signal x representing a vibration of the rotating mechanical system. In a non-limiting manner, the signal in question may be a vibration signal obtained, for example, from an accelerometer, or an acoustic signal obtained, for example, from a microphone.

[0048] Furthermore, assume that the signal x is represented as follows:

[0049] [Mathematical formula 1]

[0050] x(n)=s(n)+v(n)

[0051] where n is the sample index, s(n) is the sum of the real numbers of the sinusoids characterized by the power spectral density (DSP) with finite spectral quality at discrete frequencies, and v(n) is the number with variance The random stationary broadband component of (also real).

[0052] For example, in the case of acoustic signals, s(n) represents the periodic elements produced by the rotating machine and v(n) represents the effects of flow noise, turbulence, and transient events.

[0053] Step 203 includes determining an estimated power spectral density of signal x In the described embodiment, the DSP of a signal x is determined by computing the average periodogram of the signal. DSP Calculated according to the following equation:

[0054] [Mathematical formula 2]

[0055]

[0056] Where k is the frequency index, H is the overlap size, L is the signal size x, n is the segment size, is the number of available averages and w(n) is the apodization window used.

[0057] Furthermore, those skilled in the art will appreciate that in alternative embodiments of the monitoring method of the present invention, step 203 may include determining a time transformation of the signal x.

[0058] Step 205 comprises determining an estimated power spectral density of the noise v from the estimated power spectral density of the signal x.

[0059] This step is based on the definition of the DSP of the noise (when x(n)=v(n)), such that:

[0060] [Mathematical formula 3]

[0061]

[0062] If the signal x(n) contains a deterministic part s(n) in addition to v(n), then is not available, but can be The trend of including peak values).

[0063] So, for example, the DSP of the noise v The estimated power spectral density applied to the signal x can be obtained from This operation is expressed as the following equation:

[0064] [Formula 4]

[0065]

[0066] Where f is the smooth regression function and ∈ is the error.

[0067] In the described non-limiting example, if the model of the power spectral density of the noise is unknown (for example, when the noise is called pink noise), then the DSP Apply global regression, but along the DSP The sliding window is used for local regression.

[0068] In other words, the regression function is approximated by a base de splines and its coefficients are estimated.

[0069] These splines are continuous and have continuous derivatives called knots. The number of knots determines the bandwidth of the spline, and their conservation properties (such as the absence of moments and edge effects of the data) make them relevant for this regression. In addition, splines have the advantage of simple curvature control, and in particular, B-Splines are numerically stable.

[0070] In the example described, the regression function used is a B-spline function. Therefore, the regression first consists in constructing a B-spline basis of order o and number of knots O, and then using this basis to apply The return of.

[0071] Furthermore, in the described example, regression is applied to This is because the use of the natural logarithm makes it possible to even out the amplitude of the DSP, and because The distribution ratio of the natural logarithm of The distribution of is more symmetrical, thus making the peak robust.

[0072] Specifically, the regression is expressed as follows:

[0073] [Mathematical formula 5]

[0074]

[0075] Where B is the B-spline basis, the coefficients {c j ; j=1,...,O} are the coefficients of the B-spline.

[0076] The distribution of the error ∈(k) makes the regression more robust. The distribution used is called a heavy-tailed distribution to take into account the The presence of a peak in .

[0077] In particular, in the example shown, the distribution used is the Student distribution, but the distribution used may be another heavy-tailed distribution. More precisely, the maximum likelihood value is calculated from this distribution using an iteratively reweighted least squares (IRLS) algorithm. This algorithm consists in estimating the error weights by optimizing an objective function, wherein the error itself depends on the weighting weights. Given the estimated value {c 最小，j ; j = 1, ..., o}, so The estimated value of is finally obtained by the following calculation:

[0078] [Mathematical formula 6]

[0079]

[0080] Furthermore, in each iteration of the IRLS algorithm, samples with errors larger than a specified threshold (in this case a constant times the median of the absolute values ​​of the errors) may be truncated to obtain better robustness against peaks.

[0081] Figure 3 The estimated power spectral density of signal x is shown (in black) , the signal x represents the vibration of the same rotating mechanical system. In curve 301, the rotating mechanical system is healthy in the sense that it has not suffered any degradation (compared to the nominal state of the rotating mechanical system), while in curve 303, the rotating mechanical system has suffered degradation that has caused a change in the vibration behavior of the rotating mechanical system. In addition, the white curve represents the estimated power spectral density of the noise v obtained by the method described above.

[0082] Step 207 includes determining the equalized power spectral density Q. This is considered equalized because the power spectral density can be assimilated by the DSP of white noise. In practice, the Q-equalized DSP is obtained by the following calculation:

[0083] [Formula 7]

[0084]

[0085] The equalization can be referred to as the local oscillator frequency because the way the DSP has estimated the noise can be likened to a local sliding polynomial regression.

[0086] Advantageously, since the equalization DSP Q is comparable to the DSP of white noise, statistics designed for white noise (or for uniform DSP) can be applied to the equalization DSP. To this end, in the following steps, various operations are applied to the Q-equalization DSP instead of the estimated DSP of the signal x.

[0087] Figure 4 (in black) shows the Figure 3 DSP shown Two examples of obtained Q-balanced power spectral densities. Curve 401 corresponds to a healthy rotating machine system and curve 403 corresponds to a damaged rotating machine system.

[0088] Advantageously, the combination of the steps described so far enables a robust neutralization of the colored noise to be obtained, ie the influence of the colored noise on the signal is avoided. Furthermore, in the method, the DSP of the noise is estimated in a way that is robust to the presence of peaks.

[0089] Then step 209 includes determining a truncated power spectral density R, the truncated power spectral density being equal to the equalized power spectral density Q below the first threshold Q r A collection of values.

[0090] Since Q contains peaks in addition to spectral noise, and only the model of the probability distribution of the spectral noise is known (if x(n) = v(n), then the model of the probability distribution of the spectral noise corresponds to the probability distribution of Q(k)), while the number of peaks and their probability distribution are unknown, it is useful to truncate Q to eliminate the highest values ​​of Q, thereby retaining only the part that can be considered as spectral noise.

[0091] This allows us to use a model with a known noise distribution and thus be robust to spikes.

[0092] Therefore, R corresponds to R = {Q(k) | Q r}, that is, below the threshold Q r A collection of spectral data.

[0093] In addition, if Figure 5 As shown, the threshold Q r can be set to equalize a given percentage P of the power spectral density Q 裁剪 The value of is preserved in the truncated power spectral density R. In particular, Figure 5 In the example shown, the same percentage of 70% is applied to both the healthy case 501 and the impaired case 503. In other words, the top 30% of the Q values ​​are eliminated to obtain R.

[0094] Step 211 includes obtaining a truncated gamma distribution g​r To adjust the truncated power spectral density R.

[0095] The values ​​of Q are assumed to follow a gamma distribution when (n) = v(n). The values ​​of R are assumed to be distributed according to a right-truncated gamma distribution defined by the following equation:

[0096] [Mathematical formula 8]

[0097]

[0098] Where g(u|a,b)=Γ(a) -1 b -a e -u / b u a-1 is the full gamma distribution, Γ(.) is the full gamma function, and c is the normalization constant expressed as follows:

[0099] [Mathematical formula 9]

[0100]

[0101] in, is the upper incomplete gamma function.

[0102] The constant c can therefore be expressed as:

[0103] [Formula 10]

[0104]

[0105] c corresponds to the percentile Q assigned to the full gamma distribution r H0 assumes that there is only spectral noise (x(n) = v(n)).

[0106] The next goal is to use the known R to right-truncate the gamma distribution g r The probability density of infers the values ​​of a and b.

[0107] To this end, the approach adopted is to find the values ​​of a and b that enable the histogram of R, p(R) to be optimally adjusted in the direction of the cost function.

[0108] In the described example, this is achieved by using the Kullback-Leibler divergence as a similarity measure. This divergence is defined as follows:

[0109] [Mathematical formula 11]

[0110]

[0111] The difference between the logarithms of the probability density is weighted by p(R). This means that high probabilities are considered where the density of occurrences in the histogram is high, and in particular the optimization problem described by the following equation is solved:

[0112] [Mathematical formula 12]

[0113]

[0114] This can be simplified to:

[0115] [Mathematical formula 13]

[0116]

[0117] Finally, when Q is distributed as p(R), the cost function can be approximated using the Monte Carlo principle, so that the cost function is expressed as:

[0118] [Formula 14]

[0119]

[0120] In the described example, by making the truncated gamma distribution g r The likelihood function of the probability density is maximized using the truncated gamma distribution g r The truncated power spectral density R is adjusted.

[0121] Such a cost function is non-convex and non-differentiable. For example, this can be computed using the Nelder-Mead algorithm.

[0122] Figure 6 Shown (in grey) are histograms of R obtained from Q truncated on the right using the method described above for the healthy case 601 and the degraded case 603. Furthermore, the black curves show the adjustment of these histograms by the probability density of a right-truncated gamma distribution, also using the method described above.

[0123] Step 213 includes extracting from the truncated gamma distribution g r Determine the second detection threshold V 检测 .

[0124] Once the estimated value is obtained The method that follows first involves setting a value called the false alarm probability α. The false alarm probability represents the probability of detecting a sample as a peak when it is actually noise. The false alarm probability corresponds to the peak detection threshold, for which the higher the value of α, the greater the number of peaks detected and the greater the risk of detecting spectral noise. α is typically set between 0.01 and 0.001.

[0125] Then the threshold V consistent with the false alarm probability α is obtained by the following equation 检测 :

[0126] [Mathematical formula 15]

[0127]

[0128] Thus, once the equalization step has been performed, the detection threshold V 检测 It has nothing to do with the frequency k.

[0129] Step 215 includes extracting the detection threshold V from the equalized power spectrum density Q 检测 Above the peak Q detection.

[0130] Once you have obtained V 检测 , a statistical test is performed to decide which samples are considered peaks based on strong evidence against the H0 hypothesis. This is done by the following hypothesis test:

[0131] [Formula 16]

[0132]

[0133] Here, we assume that H0 corresponds only to The spectral noise in .

[0134] Therefore, the detection sample of Q is the peak Q 检测 .

[0135] For the healthy case 703 and the impaired case 701, the results of the hypothesis test are as follows Figure 7 As shown. Threshold V 检测 Represented as a grey line on the power spectral density (left) or probability density (right).

[0136] Finally, step 217 includes converting the extracted peak Q 检测 With reference power spectral density Q 参考 Comparisons are performed to identify faults in rotating machinery systems.

[0137] Since the proposed method enables detection of mechanical faults in a blind manner, it is useful to compare the peak values ​​detected in healthy cases with the peak values ​​detected in faulty (ie degraded) cases.

[0138] In this case, in the example described, the reference power spectral density Q 参考 corresponds to the equilibrium power spectral density obtained for the rotating mechanical system in a fault-free aircraft, and both ( Figure 8 The black Q 检测 and Figure 8 The gray Q 参考) enables fault identification.

[0139] Figure 8 Signal 801 in Figure 8 shows this comparison. In the case of a fault, the d 检测 The peaks detected in (black with “+”) contain new families of cycles that are not present in the healthy reference measurement Q 参考 The set of peaks detected in (grey with "o") does not exist in .

[0140] As a non-limiting example, one can 检测 minus Q 参考 , to highlight only anomalies / defects. Therefore, the energy criterion (ie the sum of the squared amplitudes) of the harmonics detected in this set can be used as an indication of the presence of one or more faults.

[0141] Advantageously, this enables detection of anomalies even without knowing the kinematics of the system and, if necessary, triggering of fault alarms to the extent of the entire harmonic family characterized by repetition (harmonics) and modulation (sidebands).

[0142] Advantageously, the monitoring method according to the invention is blind, ie the monitoring method does not require any a priori knowledge of the rotating mechanical system being monitored (ie the kinematics of the rotating mechanical system).

[0143] Furthermore, advantageously, the described monitoring method is robust to the use of signals with unknown DSPs with colored noise.

[0144] Advantageously, the detection threshold obtained is independent of frequency.

[0145] Furthermore, the monitoring method according to the method means that not all peaks may be monitored, but only frequencies of interest may be targeted.

[0146] Furthermore, monitoring is based on the set of frequencies that occur, unlike conventional methods that consider amplitude and sometimes phase.

[0147] Finally, the monitoring method according to the invention is also robust to application to rotating mechanical systems where the number of mechanical elements generating vibrations is unknown.

Claims

1. A method for monitoring a rotating mechanical system (101) of an aircraft, the method comprising the following steps: - acquiring (201) a signal (x) representing a vibration of the rotating mechanical system (101); - determining (203) an estimated power spectral density ( ); - From the estimated power spectral density ( ) in determining (205) the estimated power spectral density ( ); The monitoring method is characterized in that the method further comprises the following steps: - determining (207) an equalized power spectral density (Q), said equalized power spectral density being equal to the estimated power spectral density ( ) divided by the estimated power spectral density ( ); - determining (209) a truncated power spectral density (R), said truncated power spectral density being equal to the power spectral density (Q) below the first threshold (Q r )’s value set; - by truncating the gamma distribution (g r ) to adjust (211) the truncated power spectral density (R); -Based on the truncated gamma distribution (g r ) to determine (213) a second detection threshold (V 检测 ); - extracting (215) from the equalized power spectral density (Q) the 检测 ) above the peak value (Q 检测 );as well as - By extracting the peak value (Q 检测 ) and the reference power spectral density (Q 参考 ) are compared to identify (217) a fault in the rotating mechanical system (101).

2. The method according to claim 1, wherein: The signal (x) representing the vibration is a vibration signal, preferably obtained from an accelerometer, or a sound signal, preferably obtained from a microphone.

3. The method according to claim 2, wherein: The estimated power spectral density ( ) is determined by calculating the average periodogram of the signal (x).

4. A method according to any one of the preceding claims, wherein: The estimated power spectral density ( ) by applying the estimated power spectral density ( ) is determined by regression.

5. The method according to claim 4, wherein: The regression is applied to the estimated power spectral density ( ) is the natural logarithm of .

6. The method according to claim 4 or 5, wherein: The regression function used is a B-spline type function.

7. A method according to any one of the preceding claims, wherein: The first threshold (Q r ) is set so that a certain percentage (P 裁剪 ) is retained in the truncated power spectral density (R).

8. A method according to any one of the preceding claims, wherein: By truncating the gamma distribution (g r ) The adjustment of the truncated power spectral density (R) is based on making the truncated gamma distribution (g r ) is performed by maximizing the likelihood function of the probability density of .

9. A method according to any one of the preceding claims, wherein: The detection threshold (V 检测 ) is determined so that the truncated gamma distribution (g r ) is equal to a predetermined value (α) called the false alarm probability.

10. A method according to any one of the preceding claims, wherein: The reference power spectral density (Q 参考 ) corresponds to the equilibrium power spectral density obtained for the rotating mechanical system of a fault-free aircraft.

Images