Sage Journals: Discover world-class research

Abstract

In the process of prediction and health management (PHM) about remaining useful life (RUL) of aero-engine, the state evaluation of aero-engine by constructing health index (HI) plays a very important role. At present, in the research of aero-engine state evaluation, there are still many problems: the analysis of the impact between operating conditions and sensor measurement data is insufficient, the standard for sensor data sequence selection is not established, and the HI construction is not reasonable. Aiming at these problems, firstly, a data reconstruction method is used to standardize the sensor data sequences, so as to eliminate the coupling effect between operating conditions and sensor measurement data. Secondly, different from the screening of sensor data sequences by observing the trend manually, a monotonicity standard is proposed. Then, a HI calculation method based on weighted Hellinger distance is proposed, which can construct HI more reasonably and better characterize the health condition of aero-engine. Finally, based on the standard data set of aero-engine degradation process provided by NASA in 2008, a typical HI construction method based on linear regression mapping and a similarity-based RUL prediction method are selected for comparison to verify performance of this method. The results show that the proposed method has better evaluation performance, and the HI constructed by this method is also better in RUL prediction.

Keywords

Aero-engine weighted Hellinger distance state evaluation health index RUL prediction

Introduction

Aero-engine is the main core equipment supporting the normal operation of aircraft. Its reliability, maintainability, safety, supportability and testability (RMSST) have attracted much attention.¹ To improve aero-engine performance and reduce maintenance cost, prediction and health management (PHM) technology plays an important role. It is the core and key of PHM technology that PHM for the remaining useful life (RUL) of the aero-engine and guiding the aero-engine maintenance and management, through the data collected by various sensors. The RUL prediction method predicts the running time of the equipment from the current time to the failure time by analyzing the historical performance degradation trend and real-time state data of the equipment.² Among them, the state evaluation of building health index (HI) based on equipment degradation data and setting the failure threshold according to experience is the prerequisite process for predicting the RUL, and it is also one of the most important contents in PHM technology. Therefore, an appropriate state evaluation method for calculating HI is of great significance to the daily maintenance and management of aero-engines.

State evaluation judges the health state of the system through effective evaluation methods and complete operation data, to understand its health state in real-time. Generally speaking, state evaluation methods can be divided into three categories: model-based methods, data-driven methods and hybrid methods.³ For a complex electromechanical hydromagnetic coupling system such as an aero-engine, there are many performance parameters, so it is difficult to establish an accurate physical model. The data-driven method only pay attention to the data and learn from the data to evaluate the health state of the system. At present, most of the research on aero-engine state evaluation focus on data-driven methods.

Data-driven methods can be divided into univariate methods⁴ and multivariate methods⁵ according to the number of degenerate variables used. In the case of univariate, the output parameters of degenerate models (such as the Wiener model^6–11 and inverse Gaussian model^12–17) are directly used as HI. These methods only use a single parameter to calculate HI, which cannot objectively reflect the overall health state of the aero-engine and have certain disadvantages. However, in the case of multivariable, a comprehensive HI reflecting the overall health state of the aero-engine is constructed by using the multivariable data after reasonable screening. Although the multivariable data-driven method of state evaluation is still in the development stage, there have been many research results.

In literature,¹⁸ after multiple variables are selected by the relief algorithm, HI is constructed based on principal component analysis. In literature,¹⁹ HI is established by clustering and fusing multiple parameters and used for similarity RUL prediction based on the correlation vector machine modeling method. In literature,²⁰ a variety of sensor data are used to calculate the comprehensive HI, and a general HI system is established according to the comprehensive HI. In literature,²¹ a new data fusion method is proposed to establish a one-dimensional HI by automatically selecting and combining multiple sensor signals to better describe the degradation performance. In literature,²² a new signal-to-noise ratio (SNR) index is proposed, which can be adjusted according to the needs of degraded signals. Then, based on the new signal quality measure, a data level fusion model is established to construct HI by fusing multiple sensor data based on degradation. Experiments show that the HI can well describe the health status of the aero-engine, and the performance is improved when it is used to predict RUL.

Although the HI constructed by the above methods can well reflect the health state of aero-engine to a certain extent, there are still some deficiencies. First, the interaction between various working conditions and sensor data is not considered. Second, when screening useful sensor data, it is generally through manual selection, and there is no screening standard. Third, when calculating HI, the linear model is directly used to map the sensor data and HI, without considering the different weight of the degradation characteristics represented by the sensor data. Based on the above problems, the main contributions of this paper can be summarized as follows:

A data reconstruction method based on the standard core state is adopted,²³ which can eliminate the influence of working conditions when standardizing sensor data.

When screening useful sensor data, the overall monotonicity of each sensor data sequence is calculated. Then the sensor data sequences with obvious degradation trends are selected according to the monotonicity, which can better reflect the health state of aero-engine and reduce computational complexity.

A HI calculation method based on weighted Hellinger distance is proposed, which can solve the problem of simple linear mapping between sensor data and HI in the above method. Meanwhile, it also considers the different contribution of each sensor data sequence to the characterization of degradation trend.

The main contents of this paper are as follows: Section I is the introduction, which details the research significance and current situation of aero-engine state evaluation, as well as the structural arrangement of this paper. In Section II, the data reconstruction method based on the standard core state is used to denoise and standardize the data. In Section III, a state evaluation method based on weighted Hellinger distance is proposed, and the HI is calculated by using the sensor data sequence screened according to the overall monotonicity. In Section IV, the aero-engine data set provided by NASA is used for verification. In Section V, the work content and results of this paper are summarized, and the next work direction is put forward.

Data reconstruction based on standard core state

In the process of aero-engine operation, different operating states will cause the degradation trend of the sensor data sequence is not obvious, which will have a certain impact on the screening of sensor data sequence. Therefore, two problems must be solved before HI calculation: one is to solve the problem of sensor data difference caused by multiple operating conditions, the second is to select the characteristic parameters that can fully reflect the degradation performance of the system. For the first problem, this section adopts the data reconstruction method based on the standard core state to establish a data reconstruction mode process, which can be applied to single working conditions and multi-working conditions.

Data denoising

The working environment of the engine is very complex, and it is easy to produce errors in data acquisition and transmission. In particular, the data measured by the sensor generally contain some noise. Therefore, to ensure the effectiveness of the experiment, the collected initial data of the sensor should be denoised in advance. The wavelet analysis method is a commonly used denoising method in engineering. Wavelet transform can be considered as a mathematical tool that converts a signal in time domain into a different domain.

A wavelet function, which is a small wave, possesses oscillating wavelike characteristics and concentrates its energy short in time, is needed to implement the wavelet transform. Its main process is: decomposing the noisy signal into various scales after wavelet transform, removing the wavelet coefficients belonging to noise, retaining and enhancing the wavelet coefficients of the signal, and then recovering through wavelet transform to achieve the purpose of denoising. The wavelet transform of a signal $f (t)$ can be performed through a convolution operation between the signal $f (t)$ and the complex conjugate of a family of wavelets, which is expressed as²⁴:

W T_{f} (a, τ) = \frac{1}{\sqrt{a}} \int_{- \infty}^{\infty} f (t) ψ^{*} (\frac{t - τ}{a}) dt

(1)

where, $ψ (\cdot)$ is wavelet function, also known as basic wavelet or mother wavelet; $ψ^{*} (\cdot)$ represents the conjugate operation of $ψ (\cdot)$ ; $a$ is a scaling factor corresponding to frequency information; $τ$ is a translation factor corresponding to time information.

In this paper, Daubechies wavelet function is selected, and there is no analytical expression for this series of wavelet functions. At the same time, the wavelet soft threshold denoising method based on extreme value threshold estimation is selected. The extreme value threshold estimation formula is as follows.²⁵

λ = {\begin{matrix} 0.3936 + 0.1829 (\frac{\ln N}{\ln 2}), N > 32 \\ 0, N \leq 32 \end{matrix}

(2)

where, $N$ is the number of denoised data.

All sensor data are denoised by the above denoising method. Figure 1 shows the comparison of data measured by the sensor #8 of aero-engine #16 before and after denoising. The black curve is the original data, the red curve is the result of denoising by wavelet analysis, and the blue curve is the result of the typical smoothing denoising method.

Figure 1.

Denoising result of the Sensor #8 in the aero-engine #16.

Figure 1 shows that compared with the typical smoothing denoising method, the denoising processing using the wavelet analysis method has a better smoothing effect. Therefore, before calculating the HI, the wavelet analysis method is used to denoise all the original sensor data in this paper.

Data standardization

Aero-engine data is composed of conditional operating values and sensor measurement data. Different operating conditions make the aero-engine in different operating states, which will cause great differences in the data measured by the same sensor. In order to solve such problems, the data reconstruction method based on the standard core state is adopted to establish a set of data standardization mode process suitable for single working condition and multi-working conditions.

It is assumed that the degradation data of the historical aero-engine go through $E$ time observation points from the beginning to the failure, and the characteristic variables are composed of $M$ dimensional conditional operation values and $N$ dimensional sensor measurement data. The original feature sequence $x_{t_{i}}$ at $t_{i}$ is composed of:

x_{t_{i}} = (O_{t_{i}}; S_{t_{i}}) = (O_{t_{i}}^{1}, O_{t_{i}}^{2}, . . ., O_{t_{i}}^{M}; S_{t_{i}}^{1}, S_{t_{i}}^{2}, . . ., S_{t_{i}}^{N})

(3)

$E$ sets of feature sequences $x_{t_{i}} (i = 1, 2, . . ., E)$ are formed into a complete data set $X$ in chronological order.

X = [\begin{matrix} x_{t_{1}} \\ x_{t_{2}} \\ ⋮ \\ x_{t_{E}} \end{matrix}] = [\begin{matrix} O_{1}^{1} & \dots & O_{1}^{M} & S_{1}^{1} & \dots & S_{1}^{N} \\ O_{2}^{1} & \dots & O_{2}^{M} & S_{2}^{1} & \dots & S_{2}^{N} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ O_{E}^{1} & \dots & O_{E}^{M} & S_{E}^{1} & \dots & S_{E}^{N} \end{matrix}]

(4)

In the original aero-engine data set $X$ , the data obtained by the sensor at each time is always affected by the operating parameters. During operation, the system is constantly under different operating conditions, among which the operating conditions that can significantly distinguish and affect the eigenvalues of the system are called core states.²⁶ For the feature set $X$ , the feature sequence $x_{i}$ at any time $t_{i}$ can be divided into a core state. Because the influence of different operating conditions on the trend of the sensor data sequence is different in each core state, large errors will be generated when using these sensor data to directly calculate HI. Therefore, to reduce the impact of different operating conditions on the characteristic trend, the sensor data under different core states are standardized to obtain the characteristic sequence under the standard core state. The standardized formula is:

{\hat{s}}_{i, j} = \frac{s_{i, j} - μ_{i, j}}{σ_{i, j}}

(5)

where, $σ_{i, j} \neq 0$ , $s_{i, j}$ is the j ^th sensor feature sequence in the i ^th core state, $μ_{i, j}$ and $σ_{i, j}$ are the average value and standard deviation respectively, ${\hat{s}}_{i, j}$ is the value obtained after standardization by $s_{i, j}$ .

Through the above process, the feature set $X$ in the equation (4) is standardized to obtain all new sensor feature sequences $\hat{S} = ({\hat{S}}^{1}, {\hat{S}}^{2}, . . ., {\hat{S}}^{N})$ , where ${\hat{S}}^{j} = ({\hat{S}}_{1}^{j}, {\hat{S}}_{2}^{j}, . . ., {\hat{S}}_{E}^{j})^{'} (j = 1, 2, . . ., N)$ . Therefore, the new feature set $X_{new}$ :

X_{new} = [\begin{matrix} O_{1}^{1} & \dots & O_{1}^{M} & {\hat{S}}_{1}^{1} & \dots & {\hat{S}}_{1}^{N} & c s_{o} \\ O_{2}^{1} & \dots & O_{2}^{M} & {\hat{S}}_{2}^{1} & \dots & {\hat{S}}_{2}^{N} & c s_{o} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ O_{E}^{1} & \dots & O_{E}^{M} & {\hat{S}}_{E}^{1} & \dots & {\hat{S}}_{E}^{N} & c s_{o} \end{matrix}]

(6)

where, $c s_{o}$ indicates that all feature sequences have been converted to the standard core state.

HI calculation based on weighted Hellinger distance

Sensor data sequence selection

In the new feature set $X_{new}$ reconstructed by standard core state processing, not every sensor data sequence has a certain degradation trend. Therefore, in order to reduce the amount of calculation, the sensor data sequences with little or no degradation trend should be further removed. The commonly used method is to manually select useful sensor data sequences by directly observing the changing trend of sensor data, which has a certain subjectivity. Hereto, the overall monotonicity of each sensor data sequence is calculated in this paper, and the sensor data sequence is selected according to the monotonicity. The monotonicity is calculated as follows:

Mo n_{j} = \frac{1}{λ} \sum_{i = 1}^{λ} | \sum_{k = 1}^{L_{i} - 1} \frac{sgn ({\hat{S}}^{j, i} (k + 1) - {\hat{S}}^{j, i} (k))}{L_{i} - 1} |

(7)

where, $0 < Mo n_{j} < 1$ , ${\hat{S}}^{j, i}$ is the feature sequence after the standardization of the j ^th sensor in the i ^th group of historical aero-engine degradation data, $L_{i}$ is the corresponding feature sequence length, and $λ$ is the number of historical aero-engine degradation data groups. The greater the $Mo n_{j}$ , the stronger the monotonicity, indicating that the degradation trend of the j^th sensor data sequence is more obvious.

According to the monotonicity of all sensor data sequences, the sensor data sequences are divided into two categories by K-means, and the sensor data sequences with greater monotonicity are selected. The screened sensor data sequences are:

{\hat{S}}^{*} = ({\hat{S}}^{1 *}, {\hat{S}}^{2 *}, . . ., {\hat{S}}^{R *}) \subset ({\hat{S}}^{1}, {\hat{S}}^{2}, . . ., {\hat{S}}^{N}) (R < N)

(8)

At this time, the sensor data sequences obtained after screening have a certain degradation trend.

HI calculation

In probability and statistics, Hellinger distance is used to measure the similarity between two probability distributions. Assuming that there are probability distributions $P = {p_{i}}_{i \in [n]}$ and $Q = {q_{i}}_{i \in [n]}$ , the Hellinger distance between them is defined as:

h (P, Q) = \frac{1}{\sqrt{2}} \cdot ‖ \sqrt{P} - \sqrt{Q} ‖_{2}

(9)

In this paper, Hellinger distance is used to calculate HI. For the screened sensor data sequences ${\hat{S}}^{j *} = ({\hat{S}}_{t_{1}}^{j *}, {\hat{S}}_{t_{2}}^{j *}, . . ., {\hat{S}}_{t_{E}}^{j *})^{'} (j = 1, 2, . . ., R)$ , a sequence with length $L_{ref}$ from the initial time is selected as the reference sequence, which is recorded as:

{\hat{S}}_{ref}^{j *} = {({\hat{S}}_{t_{1}}^{j *}, {\hat{S}}_{t_{2}}^{j *}, \dots, {\hat{S}}_{t_{L_{ref}}}^{j *})}^{'} (j = 1, 2, \dots, R)

(10)

To calculate the similarity between a certain time and the initial time, a sequence with length $L_{ref}$ need to be selected backwards at a certain time to calculate the Hellinger distance with the reference sequence. The sequence selected at a certain time $t_{i}$ can be expressed as:

{\hat{S}}_{t_{i}}^{j *} = {({\hat{S}}_{t_{i}}^{j *}, {\hat{S}}_{t_{i + 1}}^{j *}, \dots, {\hat{S}}_{t_{i + L_{ref} - 1}}^{j *})}^{'} (j = 1, 2, \dots, R)

(11)

The Hellinger distance $h_{t_{i}}^{j *} (i = 1, 2, . . ., E; j = 1, 2, . . ., R)$ between the reference sequence represented by equation (10) and the sequence represented by equation (11) can be calculated by equation (9):

h_{t_{i}}^{j *} = \frac{1}{\sqrt{2}} ‖ \sqrt{{\hat{S}}_{ref}^{j *}} - \sqrt{{\hat{S}}_{t_{i}}^{j *}} ‖_{2}

(12)

At this time, the reference sequences of each sensor characterize the sensor measurements when the aero-engine is in the initial normal operation state, that is, these reference sequences represent that the aero-engine is in a completely healthy state. Therefore, to obtain the aero-engine health state represented by the sensor data at subsequent times, the sequence at each time can be selected backwards with step 1 and length $L_{ref}$ to calculate the Hellinger distance with the reference sequence.

In addition, the data measured by each sensor represents the health degree of the engine at a certain time to a certain extent. Therefore, at each time, the Hellinger distances calculated by all sensors are weighted according to the equation (13) to obtain a weighted value, which is used as a health index $H I_{t_{i}}$ to comprehensively characterize the aero-engine state at that time.

H I_{t_{i}} = 1 - \frac{\sum_{j = 1}^{R} {(h_{t_{i}}^{j *})}^{2}}{\sum_{j = 1}^{R} h_{t_{i}}^{j *}}

(13)

The equations (12) and (13) contain the core idea of the method based on the Hellinger distance for calculating HI. $h_{t_{i}}^{j *}$ is the Hellinger distance between the data at $t_{i}$ and the data at the initial moment, in the j ^th sensor data sequence. Then, the Hellinger distances calculated at the same moment in different sensor data sequences are weighted so as to fuse the information from multiple sensors. The closer the Hellinger distance is to 0 means the more similar the two distributions are. It is assumed that the aero-engine is normal at the initial moment and the HI is 1. Therefore, the $H I_{t_{i}}$ at $t_{i}$ is obtained by subtracting the weighted value from 1.

Evaluation index

In order to judge the state evaluation method, the overall fitting trend $F^{2}$ (the calculation method is equation (14)) and root mean square error $RMSE$ (the calculation method is equation (15)) are selected as the evaluation indexes.

F^{2} = 1 - \frac{\sum_{i = 1}^{E} {({HI}_{t_{i}}^{true} - H I_{t_{i}})}^{2}}{\sum_{i = 1}^{E} {({HI}_{t_{i}}^{true} - {HI}_{mean}^{true})}^{2}}

(14)

RMSE = \sqrt{\frac{\sum_{i = 1}^{E} {({HI}_{t_{i}}^{true} - H I_{t_{i}})}^{2}}{E}}

(15)

where, ${HI}_{t_{i}}^{true}$ represents the real HI value at the time $t_{i}$ , and ${HI}_{mean}^{true}$ represents the average value of the real HI. The smaller $RMSE$ is, the better the evaluation performance is; otherwise, the worse the evaluation performance is. The closer $F^{2}$ is to 1, the closer the HI constructed by this method is to the real HI.

Experimental verification

Data sources

The data in this paper comes from the C-MAPSS aero-engine simulation model of NASA,²⁷ and the main components of the aero-engine are shown in Figure 2. This data set has been widely used in the field of aero-engine life prediction, and its scientificity and reliability have been fully affirmed in many pieces of literature. Literature^28–30 have used this data set to complete the research of aero-engine RUL. The data set records the performance parameters of aero-engine in 24 dimensions of each flight cycle, including 3 operating condition variables and 21 sensor measurement variables (see Table 1). In each operation cycle, the aero-engine experiences the whole process from different normal states to failure, and 3 operating condition values and 21 sensor measurement values are recorded in the whole process. According to the simulated operating conditions, fault types and quantities, the data set is divided into four categories: FD001, FD002, FD003, and FD004, see Table 2 for details. Each data set is further divided into the training set and test set. At the same time, the real remaining life of the test set is given, which can be used to evaluate the performance of the prediction model. In this paper, the FD001 data set is used for experimental verification.

Figure 2.

Simplified diagram of aero-engine simulated in C-MAPSS.³¹

Table 1.

Sensor measurements of the aircraft engine.²⁷

No.	Symbol	Description	Units
1	T2	Total temperature at fan inlet	°R
2	T24	Total temperature at LPC outlet	°R
3	T30	Total temperature at HPC outlet	°R
4	T50	Total temperature at LPT outlet	°R
5	P2	Pressure at fan inlet	Psia
6	P15	Total pressure in bypass-duct	Psia
7	P30	Total pressure at HPC outlet	Psia
8	Nf	Physical fan speed	Rpm
9	Nc	Physical core speed	Rpm
10	epr	Engine pressure ratio(P50/P2)	-
11	Ps30	Static pressure at HPC outlet	Psia
12	phi	Ratio of fuel flow to Ps30	pps/psi
13	NRf	Corrected fan speed	Rpm
14	NRc	Corrected core speed	Rpm
15	BPR	Bypass Ratio	-
16	farB	Burner fuel-air ratio	-
17	htBleed	Bleed Enthalpy	-
18	Nf_dmd	Demanded fan speed	Rpm
19	PCNfR_dmd	Demanded corrected fan speed	Rpm
20	W31	HPT coolant bleed	lbm/s
21	W32	LPT coolant bleed	lbm/s

Table 2.

Description of four data sets.

Data set	Operating category	Failure mode	Number of training sets	Number of test sets
FD001	One type	One kind of failure	100	100
FD002	Six types	One kind of failure	260	259
FD003	One type	Two kinds of failure	100	100
FD004	Six types	Two kinds of failure	248	249

Experimental results and discussion

Data preprocessing and sensor data sequence choice

After obtaining aero-engine data, due to the noise interference during aero-engine operation and sensor measurement, it is generally necessary to denoise the original data in advance. According to the previous analysis, the wavelet analysis method is used to denoise the original aero-engine data. At the same time, the data is standardized through the equation (5) to eliminate the impact of different working conditions on the data. Figure 3 shows the original data sequences of sensors #7, #11, and #14 in aero-engines #20, #30, and #40. Figure 4 shows the corresponding sensor data sequences after denoising and standardization.

Figure 3.

The original data of sensors #7, #11, and #14 in the aero-engine #20, #30, and #40.

Figure 4.

The results of sensors #7, #11, and #14 in the aero-engine #20, #30, and #40 after denoising.

These figures show that the original data fluctuates greatly, and the unit range of each sensor data sequence is inconsistent. After denoising and standardization, the sensor data sequences show an obvious smoothing trend, and the unit ranges are standardized to a certain range. Moreover, Figure 4 indicates that the monotonicity trend of most sensor data sequences is consistent in different groups, but there are also sensor data sequences with large differences in monotonicity trends in different groups, such as sensor #14. Therefore, it is very necessary to calculate the overall monotonicity of the sensor data sequences according to all groups of data and take it as the standard for screening.

For the standardized data, not every sensor data sequence has an obvious degradation trend. In order to decrease computational complexity, it is necessary to select the sensor data sequence with an obvious degradation trend. The overall monotonicity of 21 sensor data sequences are calculated by equation (7), and the calculation results are shown in Figure 5. On the basis of above, using the two-classification method, the sensor data sequences are divided into two categories according to the monotonicity, and the sensor data sequences with large monotonicity are selected as the calculation data of the comprehensive HI. Selected sensor data sequences: ${\hat{S}}^{*} = ({\hat{S}}^{1 *}, {\hat{S}}^{2 *}, . . ., {\hat{S}}^{10 *}) = ({\hat{S}}^{2}, {\hat{S}}^{7}, {\hat{S}}^{8}, {\hat{S}}^{11}, {\hat{S}}^{12}, {\hat{S}}^{13}, {\hat{S}}^{15}, {\hat{S}}^{17}, {\hat{S}}^{20}, {\hat{S}}^{21})$ . Figure 6 shows the sensor data sequences selected after a group of aero-engine data screening. The selected 10 sensor data sequences have an obvious monotonic trend, which can better reflect the degradation trend of aero-engine.

Figure 5.

Overall monotonicity of all sensor data sequences.

Figure 6.

Sensor data sequences selected from the aero-engine #54.

Comprehensive HI calculation and evaluation

There is no real HI data in the aero-engine data set, but according to the existing literature, the general engine degradation trend is exponential degradation. Therefore, in this paper, it is assumed that HI at all initial times is 1, and it is exponentially degraded to failure time according to a certain degradation slope, that is, HI is 0. In this regard, the HI calculated according to the following equation (16) is used as the real HI in equations (14) and (15) to evaluate the HI calculation method based on the weighted Hellinger distance.

H I_{t} = 1 + a (1 - e^{b \cdot t})

(16)

where, $a$ and $b$ are the constants affecting the degradation trend. In this paper, to simplify the calculation, $a = e^{3.5} - 1$ and $b$ is calculated according to the standard that HI of each group of data decreases to 0 at the last time.

At the same time, the method of directly fusing sensor data based on a linear regression model to calculate HI is selected as the comparison.³² The expression of linear regression model is:

{\hat{h}}_{t_{i}} ({\hat{S}}_{t_{i}}^{j *}, θ) = θ_{0} + \sum_{j = 1}^{R} θ_{j} {\hat{S}}_{t_{i}}^{j *}

(17)

where, ${\hat{S}}_{t_{i}}^{j *}$ is the j ^th sensor data after standardization at $t_{i}$ , ${\hat{h}}_{t_{i}} ({\hat{S}}_{t_{i}}^{j *}, θ)$ is the HI calculated by the method based on linear regression model, $R$ is the number of sensor data sequences after screening. And, weight vector $θ = [θ_{0}, θ_{1}, \dots, θ_{R}]$ can be estimated using an iterative weighted least squares algorithm.³³

Figure 7 is a comparison diagram of HI (red curve) estimated by the method proposed in this paper, HI (black curve) calculated by the method based on linear regression model and real HI (blue curve) calculated according to equation (16) for a group of aero-engine data. It can be seen that the HI calculated by the method proposed in this paper and the HI obtained by the method based on linear regression model can better fit the real HI curve, which verifies the feasibility of the method proposed in this paper. At the same time, we can find that the HI calculated by the method based on linear regression model at the last failure time of the aero-engine is not zero, indicating that the HI calculated by the method at this time does not well reflect that the aero-engine has degenerated to failure, but the method proposed in this paper has some improvement in this aspect.

Figure 7.

The HI curves are calculated with the data of aero-engine #4.

In addition, using all the training data of FD001, the HI is calculated by the method proposed in this paper and the method based on linear regression model respectively, then the overall fitting trend $F^{2}$ and root mean square error $RMSE$ of the two methods are analyzed. Figures 8 and 9 show the curves of $F^{2}$ and $RMSE$ calculated by the method proposed in this paper. In Figure 8, the minimum $F^{2}$ is 0.6556, the maximum $F^{2}$ is 0.9807, and the average value is 0.8921. In Figure 9, the minimum $RMSE$ is 0.0371, the maximum $RMSE$ is 0.1569, and the average value is 0.0832.

Figure 8.

The overall fitting trend $F^{2}$ of HI calculated by the method proposed in this paper (the method 1) and the method based on linear regression model (the method 2).

Figure 9.

The root mean square error $RMSE$ of HI calculated by the method proposed in this paper (the method 1) and the method based on linear regression model (the method 2).

From Figures 8 and 9, in some groups of aero-engine data, the calculated HI is quite different from the real HI at some time, but the rule of overall degradation trend is similar for both HI. At the same time, as a comparison, the method based on linear regression model is used to calculate HI, and the evaluation indexes are shown in the Table 3. On the whole, the state evaluation performance of the method proposed in this paper is better.

Table 3.

Comparison of the results of the two methods.

Evaluation method	Average $RMSE$	Average $F^{2}$
The method proposed in this paper	0.0832	0.8921
The method based on linear regression model	0.1098	0.7973

Further analysis

From the previous experimental results and analysis, it can be seen that the HI calculation method based on weighted Hellinger distance proposed in this paper has good results, and the aero-engine degradation trend can be well fitted by the calculated HI. However, in the PHM of aero-engine, the construction of HI is generally for the purpose of RUL prediction. Therefore, in this section, the performance of HI calculated by the method proposed in this paper is further verified by a set of RUL prediction experiments. In order to test and compare the role of HI calculated by the two methods (the method proposed in this paper and the method based on linear regression model) in the process of RUL prediction, a similarity RUL prediction method based on dynamic time warping (DTW) distance²³ is selected for experiment, and its basic idea is shown in the Figure 10.

Figure 10.

Flow chart of similarity RUL prediction based on DTW distance.

In the DTW-based life prediction, the similarity between the history sample and the test sample is calculated by DTW distance, and the RUL prediction of the test sample is further based on the actual RUL of the history sample. Assuming that the end moments of the test sample and the $l^{th}$ history sample are $t_{o} = k \cdot Δ t$ and $t_{l} = q \cdot Δ t$ , respectively ( $Δ t$ is the data acquisition interval of the aero-engine sensors), the HI used to calculate the similarity in the test sample and the $l^{th}$ historical sample are:

H_{o} (k) = (h_{o} (1 \cdot Δ t), h_{o} (2 \cdot Δ t), . . ., h_{o} (k \cdot Δ t))

(18)

H_{l} (q) = (h_{l} (1 \cdot Δ t), h_{l} (2 \cdot Δ t), . . ., h_{l} (q \cdot Δ t))

(19)

where, $h_{o} (k \cdot Δ t)$ denotes the HI of the test sample at time $k \cdot Δ t$ , $h_{l} (q \cdot Δ t)$ denotes the HI of the $l^{th}$ history sample at time $q \cdot Δ t$ .

To calculate the similarity of two HI sequences, their shortest distance is firstly calculated using the DTW method. The shortest distance $D (k, q)$ between $H_{o} (k)$ and $H_{l} (q)$ can be calculated by equation (20).

D (i, j) = {\begin{matrix} D (i, j) + d_{i, j} \\ D (i, j) + d_{i, j} \\ D (i - 1, j - 2) + d_{i, j} \end{matrix}

(20)

where, $d_{i, j} (i = 1, 2, . . ., k; j = 1, 2, . . ., q)$ can be calculated by the following equation (21), and the shortest path must be searched continuously from the starting point $d_{1, 1}$ to the end point $d_{m, n}$ , that is, when the current distance is $d_{i, j}$ , the next search path can only be $d_{i + 1, j}$ or $d_{i, j + 1}$ or $d_{i + 1, j + 1}$ .

d_{i, j} = \sqrt{{(h_{o} (i \cdot Δ t) - h_{l} (j \cdot Δ t))}^{2}}

(21)

Then, the exponential function based on the DTW distance is chosen as the similarity function to calculate the similarity between the test sample and the $l^{th}$ history sample, and the similarity function is:

S_{o, l} = e^{- {(D (k, q))}^{2}}

(22)

After calculating the similarity between the test sample and the $l^{th}$ history sample according to equation (22), the most similar moment ${\overset{⌢}{t}}_{l} = \overset{⌢}{q} \cdot Δ t$ corresponding to the end point $t_{o} = k \cdot Δ t$ of the test sample in the $l^{th}$ history sample is calculated by equation (23).

\overset{⌢}{q} = \arg S_{o, l}

(23)

Therefore, the actual RUL $AR L_{l} ({\overset{⌢}{t}}_{l})$ of the $l^{th}$ history sample at the time ${\overset{⌢}{t}}_{l}$ can be obtained according to equation (24), which can be used to predict the RUL $PR L_{l}$ of the test sample using the $l^{th}$ history sample.

PR L_{l} = AR L_{l} ({\overset{⌢}{t}}_{l}) = (q - \overset{⌢}{q}) \cdot Δ t

(24)

Finally, weights are calculated based on the similarity between each history sample and the test sample, and the RUL predicted by all history samples are weighted and summed to obtain the final RUL of the test sample.

PRL = \sum_{l = 1}^{L} w_{l} \cdot PR L_{l}

(25)

where, $w_{l}$ is the weight calculated based on the similarity of the $l^{th}$ history sample to the test sample.

In the process of RUL prediction, the HI calculated by the two methods are used as the reference HI curves, and the test data set in FD001 is used for experiment of RUL prediction, then the prediction error is calculated. As shown in the Figure 11, the error curves of RUL prediction based on two calculation methods of HI are shown, in which the method #1 refers to using the method proposed in this paper to calculate HI for RUL prediction, and method #2 refers to using the method based on linear regression model to calculate HI for RUL prediction.

Figure 11.

RUL prediction errors of two methods for calculating HI.

In the Figure 11, the maximum error of the method #1 is −97.2337, the minimum error is 0.9597, and the average value is 0.5357; the maximum error of the method #2 is −116.4621, the minimum error is 0.2828, and the average value is −32.8302. By further calculating the root mean square error of the predicted results, it is obtained that the root mean square error of the method #1 is 27.6250 and the root mean square error of the method #2 is 46.3806.

Figure 11 shows that in the similarity RUL prediction process based on DTW distance, the overall effect of RUL prediction obtained by using the HI calculated by the method proposed in this paper is better than that obtained by using the HI calculated by the method based on linear regression model. However, when the HI calculated by the method proposed in this paper is used in the RUL prediction based on DTW distance similarity, there are some cases with large prediction errors of test samples. A possible explanation for this is that in the process of calculating the HI of historical samples, there are cases where the $F^{2}$ of HI in some historical samples is small and the $RMSE$ is large, see Figures 8 and 9. This situation continues to have an impact in the process of similarity RUL prediction based on DTW distance.

In addition, through the previous analysis, compared with the method based on linear regression model, the method proposed in this paper improves accuracy in failure prediction to a certain extent. Due to the sensor data screening according to monotonicity, this method also reduces the computational complexity and the cost of condition monitoring. Moreover, considering the influence of different working conditions on sensor data, data reconstruction based on standard core state is carried out to improve the applicability of the method to different working condition situations. And, the HI calculation method based on weighted Hellinger distance does not require the construction of various mapping models between sensor data and HI, which improves system reliability.

Conclusion

Aiming at the problem of aero-engine state evaluation, a HI calculation method based on weighted Hellinger distance is proposed in this paper. Firstly, wavelet analysis denoising and data reconstruction based on the standard core state are used to preprocess the original aero-engine data, so as to eliminate the influence of noise and conditional operation value on the sensor measurement data. Secondly, the overall monotonicity of the sensor data sequences are calculated, and the sensor data sequences with relatively obvious degradation trend are selected according to this standard. Then, the HI is calculated based on the weighted Hellinger distance, which can well evaluate the state of the aero-engine. Finally, the aero-engine simulation data set released by NASA in 2008 is used to verify the feasibility and excellent evaluation performance of the method proposed in this paper.

From the previous analysis, the method proposed in this paper can eliminate the influence of different working conditions on sensor data, and screening data according to monotonicity can reduce the computational complexity, then the HI calculated based on weighted Hellinger distance can better characterize the health state of aero-engine without considering various mapping relationships between sensor data and HI. Of course, there are shortcomings here. First, monotonicity is not the only kind of changing trend of sensor data, and other information should be further explored in order to more comprehensively characterize the health state of aero-engine. Second, the method based on the weighted Hellinger distance is affected by the calculating step and the state of aero-engine at the initial moment, which will also be investigated to improve the robustness of the method in the future.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Key R&D Program of China Under Grant 2016YFB1200600.

ORCID iD

Biao Yang

References

Liu

Lei

Pan

, et al. Prediction of remaining useful life of multi-stage aero-engine based on clustering and LSTM fusion. Reliab Eng Syst Saf 2021; 214: 107807.

Huang

Shi

, et al. Research progress of data-driven life prediction and health management technology. Inf Control 2017; 42: 72–82.

Wang

Chen

, et al. Equipment health assessment method based on unlabeled, unbalanced and uncertain initial value data. Control Decis Making 2020; 35: 2687–2695.

Tamilselvan

Wang

. Failure diagnosis using deep belief learning based health state classification. Reliab Eng Syst Saf 2013; 115: 124–135.

Javed

Gouriveau

Zerhouni

. A new multivariate approach for prognostics based on extreme learning machine and fuzzy clustering. IEEE Trans Cybern 2015; 45: 2626–2639.

Shen

Wang

, et al. On modeling bivariate wiener degradation process. IEEE Trans Reliab 2018; 67: 897–906.

Shao

, et al. An adaptive and generalized Wiener process model with a recursive filtering algorithm for remaining useful life estimation. Reliab Eng Syst Saf 2022; 217: 108099.

Kim

, et al. A nonlinear-drift-driven wiener process model for remaining useful life estimation considering three sources of variability. Reliab Eng Syst Saf 2021; 212: 107631.

Lei

Yan

, et al. A Wiener-process-model-based method for remaining useful life prediction considering unit-to-unit variability. IEEE Trans Ind Electron 2019; 66: 2092–2101.

10.

Zhu

Zuo

Cai

. Performance reliability prediction of civil aviation engine based on Wiener process. J Aeronaut Power 2013; 28: 1006–1012.

11.

Peng

Zhou

Sun

, et al. Residual life prediction of products based on degradation and life data fusion. Syst Eng Electron Technol 2011; 33: 1073–1078.

12.

Sun

Wang

, et al. An improved inverse Gaussian process with random effects and measurement errors for RUL prediction of hydraulic piston pump. Measurement 2021; 173: 108604.

13.

Chen

Sun

, et al. Remaining useful life prediction based on an adaptive inverse Gaussian degradation process with measurement errors. IEEE Access 2020; 8: 3498–3510.

14.

Zhao

Sun

Qiu

, et al. Optimal inspection and mission abort policies for systems subject to degradation. Eur J Oper Res 2021; 292: 610–621.

15.

Guida

Pulcini

. The inverse gamma process: A family of continuous stochastic models for describing state-dependent deterioration phenomena. Reliab Eng Syst Saf 2013; 120: 72–79.

16.

Peng

Zhu

Shen

. The transformed inverse Gaussian process as an age- and state-dependent degradation model. Appl Math Model 2019; 75: 837–852.

17.

Sujith

Balakrishna

. Autoregressive inverse Gaussian process and the stochastic volatility modeling. Commun Stat Theory Methods Epub ahead of print 16 September 2021. DOI: 10.1080/03610926.2021.1977324.

18.

Zhang

Wang

, et al. Remaining useful life prediction for aero-engine based on the similarity of degradation characteristics. Syst Eng Electron Technol 2019; 41: 1414–1421.

19.

Shi

Wang

, et al. Health index synthetization and remaining useful life estimation for turbofan engines based on run-to-failure datasets. Eksploatacja Niezawoonosc-Maintenance Relib 2016; 18: 621–631.

20.

Wang

Youn

. A generic probabilistic framework for structural health prognostics and uncertainty management. Mech Syst Signal Process 2012; 28: 622–637.

21.

Kim

Song

Liu

. A generic health index approach for multisensor degradation modeling and sensor selection. IEEE Trans Autom Sci Eng 2019; 16: 1426–1437.

22.

Liu

Chehade

Song

. Optimize the signal quality of the Composite Health ã Index via data fusion for degradation modeling and prognostic analysis. IEEE Trans Autom Sci Eng 2017; 14: 1504–1514.

23.

Wang

. Research on key technologies of fault prediction and health management of suspension system of medium and low speed maglev train. PhD Thesis, National University of Defense Technology, China, 2021.

24.

Yan

Gao

Chen

. Wavelets for fault diagnosis of rotary machines: A review with applications. Signal Process 2014; 96: 1–15.

25.

Zhang

Liu

. Laser doppler signal denoising based on wavelet packet thresholding method. Int J Optics 2019; 2019: 1–11.

26.

Meng

Xie

, et al. Study on health assessment method of a braking system of a mine hoist. Sensors 2019; 19: 769.

27.

Frederick

DeCastro

Litt

. Technical manual, user’s guide for the commercial modular aero-propulsion system simulation (C-MAPSS). Cleveland, OH: National Aeronautics and Space Administration, John H. Glenn Research Center at Lewis Field, 2007.

28.

Wang

Zhu

, et al. A data-driven degradation prognostic strategy for aero-engine under various operational conditions. Neurocomputing 2021; 462: 195–207.

29.

Miao

. A deep domain adaptative network for remaining useful life prediction of machines under different working conditions and fault modes. IEEE Trans Instrum Meas 2021; 70: 1–14.

30.

Zhang

Xin

Liu

, et al. Health status assessment and remaining useful life prediction of aero-engine based on BiGRU and MMoE. Reliab Eng Syst Saf 2022; 220: 108263.

31.

Saxena

Goebel

Simon

, et al. Damage propagation modeling for aircraft engine run-to-failure simulation. In: 2008 International conference on prognostics and health management, Denver, CO, 6–9 October 2008, pp.1–6. New York, NY: IEEE.

32.

Xiong

Cheng

. Remaining useful life prognostics of aircraft engines based on damage propagation modeling and data analysis. In: 2015 8th International symposium on computational intelligence and design, Hangzhou, China, 12–13 December 2015. New York, NY: IEEE.

33.

Tao

Hua

, et al. An iterated reweighting total least squares algorithm formulated by standard least-squares theory. Surv Rev 2021; 53: 454–463.

An aero-engine state evaluation method based on weighted Hellinger distance

Abstract

Keywords

Introduction

Data reconstruction based on standard core state

Data denoising

Data standardization

HI calculation based on weighted Hellinger distance

Sensor data sequence selection

HI calculation

Evaluation index

Experimental verification

Data sources

Experimental results and discussion

Data preprocessing and sensor data sequence choice

Comprehensive HI calculation and evaluation

Further analysis

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References