Bivariate-Dependent Reliability Estimation Model Based on Inverse Gaussian Processes and Copulas Fusing Multisource Information

: Reliability estimation for key components of a mechanical system is of great importance in prognosis and health management in aviation industry. Both degradation data and failure time data contain abundant reliability information from different sources. Considering multiple variable-dependent degradation performance indicators for mechanical components is also an effective approach to improve the accuracy of reliability estimation. This study develops a bivariate-dependent reliability estimation model based on inverse Gaussian process and copulas fusing degradation data and failure time data within one computation framework. The inverse Gaussian process model is used to describe the degradation process of each performance indicator. Copula functions are used to capture the dependent relationship between the two performance indicators. In order to improve the reliability estimation accuracy, both degradation data and failure time data are used simultaneously to estimate the unknown parameters in the degradation model based on the likelihood function transformed using the zeros-ones trick. A simulation study and a real application in the reliability estimation of mechanical seal used in airborne hydraulic pump are conducted to validate the effectiveness and accuracy of the proposed model compared with existing reliability models.


Introduction
Accurate reliability estimation for critical components of a mechanical system is crucially important in the aviation industry [1,2]. Reliability information can be obtained from various sources, including degradation data and failure lifetime data [3]. Failure lifetime data represent the reliability information for mechanical component on the time scale. However, sufficient lifetime data are very difficult to obtain only by life testing or even accelerated life testing for those mechanical components with the characteristics of long life and small sample [4]. Degradation data describe the entire degradation process from start to failure [5]. Degradation data can be monitored and collected by sensors. Making full use of available reliability information of mechanical components, including degradation data and failure lifetime data, can improve the accuracy of the reliability estimation and guide optimal maintenance [6].
According to Zhang et al. [7], three main types of reliability estimation models have been developed in recent years, namely knowledge-based, physical-based and data-driven models. A data-driven model is much more flexible to implement and perform the reliability estimation of mechanical components [8]. Regarding degradation data-based models, a stochastic process degradation model is commonly used to build the degradation model studies aimed to solve the problem of evaluating the reliability of system components or systems. However, most of them only consider degradation or failure time data or consider only one performance indicator or two independent performance indicators. Few of them contribute to the model considering both dependent performance indicators and combining different types of data source.
To solve this problem and improve the accuracy of the reliability evaluation and life prediction for lifecycle management, this study develops a reliability evaluation method by integrating likelihood functions with both degradation data and failure lifetime data considering two dependent performance indicators. First, the degradation trajectories of the two performance indicators of the components are both described by the IG process, and copula functions are used to describe the dependency relationship between the two performance indicators. Different from the existing bivariate copula reliability models, both the degradation data and failure lifetime data are integrated in order to improve the accuracy of the reliability estimation and remaining useful life (RUL) prediction. The zerosones trick and Bayesian Markov chain Monte Carlo (MCMC) method are used to estimate the unknown parameters. The Akaike information criterion (AIC) and BIC are used to select the best copula function linking two marginal distributions of the two performance indicators. Finally, the reliability estimation and RUL prediction are performed based on the proposed model and it is validated by the simulation study and mechanical seal degradation experiments.
The rest of this paper is organized as follows. In Section 2, the reliability estimation model based on the IG process and copula function is presented. Section 3 shows the likelihood functions built for two cases, where the first one only considers degradation data and the second one considers both degradation data and failure lifetime data. Parameter estimation procedures are performed based on different likelihood functions. In Section 4, a simulation study is described to illustrate the accuracy of the proposed model. In Section 5, the degradation data and failure lifetime data from mechanical seal are used to validate the proposed model. In Section 6, the overall conclusions of this study are summarized.

The IG Process Model
IG process is a monotonically increasing process and it has been widely applied in the monotonic degradation process for mechanical components [36,37]. In our study, two degradation performance indicators are assumed to be modeled by the IG process. The IG process {X(t), t ≥ 0} has the following properties: (1) the initial value X(0) ≡ 0 with probability one; (2) the degradation increments are independent for separate time intervals.
(3) the degradation increment ∆X(t) = X(t + ∆t) − X(t) follows the IG distribution: 2 (1) where λ, η > 0, Λ(t) is the monotonic time scale function representing the degradation trajectory, and ∆Λ(t) = Λ(t + ∆t) − Λ(t). It can be expressed as a power-law function, exponential function, obtained by the failure mechanism model or fitted by the degradation data [38,39]. In particular, Λ(t) = t q is commonly used. When q = 1, it means the degradation trajectory is linear. When q > 1, it means that the degradation trajectory is concave. When q < 1, it means that the degradation trajectory is convex. Since Λ(t) = t q can describe different types of degradation paths, we will adopt it without the loss of generality. The cumulative distribution function (CDF) and probability destination function (PDF) of ∆X(t). ∆X(t) can be obtained by: where Φ(·) is the standard normal CDF. The failure time T denotes the time at which the degradation path first reaches a predefined threshold d, and it can be expressed as: Considering that X(t) ∼ IG λΛ(t), η(Λ(t)) 2 , the CDF and PDF of T can be obtained by If ηΛ(t) is large, i.e., when t is large enough, X(t) is approximately normally distributed with mean λΛ(t) and variance λ 3 Λ(t) η [40]. The CDF and PDF of the failure time T can be approximately expressed as: The RUL L j (j = 1, 2, · · · , M) denotes the time t j that the component can operate until it fails. It also means that if there is X t j < d, the component continues functioning before time t j . M is the total measurement times. The RUL can be expressed as: Based on the concept of failure time T, RUL L j at time t j can be considered as the failure time of the IG process X l j , l j > 0 crossing the threshold d − X t j , where X l j = X l j + t j − X t j . Based on Equations (7) and (8), the CDF and PDF of the RUL can be, respectively, expressed as: The IG process-based degradation model for each performance degradation indicator will be used as marginal input function in copulas. The following sections will show how to analyze the dependent relationship between two performance indicators using copulas.

Copula Theory and Bivariate Degradation Model
The bivariate copula function is used to construct the joint function of two degradation performance indicators taking their dependency into account. According to Sklar's theory [41], a bivariate copula function can be expressed as: where F 1 (∆X 1 ) and F 2 (∆X 2 ) are the marginal distributions that are uniform on [0, 1], and θ is the parameter in the copula function to evaluate the dependency strength of the two performance indicators. The density function can be expressed as: , the bivariate copula can be considered as the mapping from [0, 1] 2 to [0, 1] Five copula functions commonly used in the reliability estimation are shown in Table 1. The Gaussian copula is an elliptical copula function, which can describe the symmetry dependence relationship between two performance indicators. The Clayton, Frank and Gumbel copulas are typical Archimedean copulas. There is only one parameter governing the degree of dependence in Archimedean copulas, which makes them very flexible for reliability modeling. The Farlie-Gumbel-Morgenstern (FGM) is based on generalizations of the FGM copula.

C(u,v|θ)
The Domain of θ Assuming that the component fails when either of the two performance indicators X 1 (t) and X 2 (t) cross the threshold d 1 and d 2 , respectively, the reliability function considering two performance indicators can be given as follows: The RUL L j (j = 1, 2, · · · , M) of the mechanical component considering two performance indicators can be expressed as: Aerospace 2022, 9, 392 6 of 25 L j = inf l j : X 1 t j + l j ≥ d 1 and X 2 t j + l j ≥ d 2 (15) If two performance indicators are dependent, and the dependency structure is described by copula functions, Equation (14) can be expressed as [42,43]: Furthermore, the PDF function can be obtained by: Regarding the different copula functions, the f s (t) function is different.
The mean time to failure (MTTF) of the mechanical component can be expressed as: The AIC and BIC are used in our study to choose a suitable copula function to describe the dependent relationship between two performance indicators. Both the AIC and BIC deal with the balance between the goodness of fit of the model and the model complexity. They can avoid the risk of overfitting or underfitting problems during model selection. The AIC and BIC can be, respectively, given by: where L denotes the likelihood function of the model, k denotes the number of parameters in the model, and n denotes the number of observations or the sample size. Based on discussion above, we complete the bivariate-dependent reliability analysis model for mechanical components using IG process model and copulas. However, there are unknown parameters which need to be estimated. In addition, the estimation of unknown parameters will affect the accuracy of reliability and RUL analysis. In the following sections, the estimation of unknown parameters fusing degradation data and failure lifetime data will be introduced.

Motivation
Degradation data provide the reliability information for mechanical component during the whole degradation process, and failure lifetime data provide the reliability information from a time-scale perspective. The IG process can describe the degradation process based on degradation data and failure time expression based on IG process given predefined threshold can describe the time-related reliability information. To make full use of degradation data and lifetime data, we integrate the likelihood functions for two types of reliability information under Bayesian information fusion framework. Previous studies have shown the effectiveness of reliability information fusion framework [30,31]. However, to the best of our knowledge, few studies focus on the consideration of both dependent multiple performance indicators and reliability information fusion. As we have discussed above, copula functions link degradation performance indicators considering dependent relationships among them. It is meaningful to combine both degradation data and failure lifetime data with multiple performance indicators' reliability information to obtain accurate reliability estimation. The posterior distributions of unknown parameters are obtained Aerospace 2022, 9, 392 7 of 25 by the integration likelihood function using the Bayesian MCMC method. We will show the basic assumptions, parameter settings and the whole information fusion framework in detail in the following sections.
Assuming that there are N test samples in our test and M i (i = 1, 2, · · · , N) measurements during the predefined test time. The number of failed samples is N f , and the number of censored samples is N c . There is N = N f + N c , evidently. Let X k t ij denote the degradation data for the k(k = 1, 2) th performance indicator of the i(i = 1, 2, · · · , N) th sample at the j(j = 1, 2, · · · , M i ) th measurement. Then, the degradation increment ∆X k t ij = X k t ij − X k t i,j−1 . Let Θ = Θ IG k , θ denote all the unknown parameters in the degradation model, where Θ IG k = (λ k , η k , q k ) denotes the unknown parameters in the IG process model and θ denotes the dependence structure parameter in the copula function. Let T = (T 1 , T 2 , · · · , T N ) denote the failure time data or censoring time data for each test sample. Let D = (D 1 , D 2 ) denote the degradation data, where D 1 denotes the degradation data for performance indicator 1 and D 2 denotes the degradation data for performance indicator 2. The likelihood function based on degradation data is shown in Section 3.2, and the likelihood function based on two types of data, namely degradation data and failure lifetime data, is shown in Section 3.3.

Likelihood Function Based on Degradation Data
The likelihood function is constructed by the IG process model using available degradation data. We develop the likelihood function from the following three aspects: (1) Considering one performance indicator If only one degradation performance indicator is considered, the log-likelihood function based on the degradation data of one performance indicator can be given by Equation (22).
Considering that the degradation process of the performance indicator is described by the IG process, then Equation (21) can be expressed as: (2) Considering two independent performance indicators If two independent performance indicators are considered, the log-likelihood function based on the degradation data of both performance indicators can be given by Equation (24).
(3) Considering two dependent performance indicators If two dependent performance indicators are considered and the dependency structure is described by the copula function, the log-likelihood function based on the degradation data of both performance indicators can be given by Equation (25).

Likelihood Function Fusing Degradation Data and Failure Lifetime Data
As we have discussed above, it is reasonable to integrate both degradation data and failure lifetime data to obtain accurate reliability estimation. The likelihood function based on degradation data is developed by the IG process model, and the likelihood function based on failure lifetime data is developed by the failure lifetime expression based on IG process model. Similarly, we also present the likelihood function fusing degradation data and failure lifetime data from the following aspects: (1) Considering one performance indicator If only one performance indicator is considered, the log-likelihood function based on both the degradation data and failure lifetime data of one performance indicator can be given by Equation (26).
If the degradation process of each performance indicator is described by the IG process, then Equation (25) can be expressed as: (2) Considering two independent performance indicators If two independent performance indicators are considered, the log-likelihood function based on both the degradation data and failure lifetime data of the two performance indicators can be given by Equation (28).
(3) Considering two dependent performance indicators If two dependent performance indicators are considered and the dependency structure is described by the copula function, the log-likelihood function based on both the degradation data and failure lifetime data of two performance indicators can be given by Equation (29).
In Sections 3.2 and 3.3, the likelihood functions based on degradation data and two types of data are constructed, respectively. In the following section, the estimation methods of unknown parameters based on likelihood functions will be shown.

Unknown Parameter Estimation Methods
Two unknown parameter estimation methods are widely used, namely the maximum likelihood estimation (MLE) method and Bayesian MCMC method. The MLE method is widely used in parameter estimation. The aim of MLE is to find the values of unknown parameters that maximize the likelihood function or log-likelihood function over the parameter space based on the observed data Y. Taking the derivative of the log-likelihood function and setting it to 0, the estimation results can be obtained by If there is no closed-form solution to the first order derivative of the log-likelihood function, the MLE can be found by numerical methods.
(2) The Bayesian MCMC method The Bayesian MCMC method is used to estimate the unknown parameters in the model based on the likelihood function or log-likelihood function. The Bayesian MCMC method constructs a Markov chain with the corresponding prior information and generates a sample of the posterior distribution by recording the states from the chain. Enough steps need to be performed to obtain the posterior distribution which is closer to the actual posterior distribution. In this study, the whole Bayesian MCMC process is performed in the OpenBUGS software using the Gibbs sampling algorithm. The Bayesian estimation of unknown parameters can be expressed by where π(Θ) and π(Θ|Y) denote the prior distribution and posterior distribution of unknown parameters Θ, respectively. Prior distributions describe the uncertainty about the unknown parameters in the model. A prior distribution might be a non-informative prior distribution or an informative prior distribution depending on the real applications. The non-informative prior distribution provides the general information about the unknown parameters in the model. If available information about the parameters in the model has been obtained, the informative prior distribution might be under consideration. The Gelman-Rubin ratio is used to evaluate the MCMC convergence by analyzing the difference among Markov chains. The convergence is checked by comparing the between-chain and within-chain variances for each unknown parameter.
Based on the different types of likelihood functions given above, different strategies are performed to separately estimate the unknown parameters.

•
Regarding the likelihood functions (21) and (23), the Bayesian MCMC method is also used to estimate the unknown parameters in the likelihood function. Because the estimation results for one performance indicator have no effect on the estimation results for the other performance indicator, the unknown parameters for each performance indicator can be estimated separately. In addition, Wang et al. [44] has reported the explicit estimation method and results for the IG process.

•
As for the likelihood function (24), two stage estimation strategies, which are called the inference functions for the margins (IFM) method, are used to separately estimate the unknown parameters in the IG process and copula function. The likelihood corresponding to the IG process ln L IG Θ IG k D is estimated first. The margin distributions of the two performance indicators are calculated based on the estimation results, and then used to estimate the unknown parameter in the copula function based on the likelihood function. Details of this procedure have been previously reported [19].
Regarding the likelihood functions (25), (27) and (28), the zeros-ones trick is used to estimate unknown parameters since two different types of reliability information cannot be estimated together. The zeros-ones trick can be used to remodel the likelihood function fusing multiple sources of data using the Bernoulli distribution. That means we treat both lifetime data and degradation data as Bernoulli data in the transformed integration likelihood function. Taking the log-likelihood function Equation (28) as an example, it can be transformed as: M i denotes the number of degradation data samples; N f denotes the number of failed samples; N c denotes the number of censored samples; f B (1; exp(l i ), 1) denotes the binomial probability function with success probability exp(l i ) and N = 1. l i denotes the transformed likelihood function which can be given by: Therefore, the likelihood function in the model can be regarded as the product of the densities of new pseudorandom variables with all observed values that are equally set as one. The new random variables follow the Bernoulli distribution with success probability exp(l i ). To ensure that the probability is lower than one, exp(−C) is multiplied to each likelihood term, where C is a positive large number. Then, the likelihood function (33) can be given by: Due to the calculation complexity of those likelihood functions, the second estimation method, i.e., the Bayesian MCMC method, is performed to estimate the unknown parameters. In this study, the MLE estimation resultsΘ IG = λ ,η,q based on the IG process model and copula function can be obtained as discussed above, and the prior distribution π(Θ) of unknown parameters Θ could be chosen as the informative distribution. Informative prior distributions in reliability estimation under a Bayesian framework have been used in several studies [45,46]. Liu et al. [47] used the different priors for parameters in the stochastic model, which is also applied in our study. As for the parameters θ in the copula function, normal distributions are selected as prior distributions.
After all the unknown parameters are obtained, the RUL, reliability estimation value and MTTF can be calculated. We also present a flow chart to show the whole reliability analysis procedure, as shown in Figure 1.

Simulation Study
In this section, a simulation study is conducted to verify the accuracy of the proposed model.

Simulation Study
In this section, a simulation study is conducted to verify the accuracy of the proposed model.

Data Generation
The degradation of two performance indicators is described by the IG process with parameters λ 1 = 3, η 1 = 24, q 1 = 1.2 and λ 2 = 2, η 2 = 15, q 2 = 1.4. The failure threshold for each performance indicator is d 1 = 15 and d 2 = 12, respectively. The time interval is assumed to be ∆t = 0.2. The dependence structure is described by the Frank copula, and the dependent parameter θ = 12. A total of 21 simulation samples with random degradation data and failure lifetime data fit for the proposed model are generated by the Monte Carlo method. The degradation path for each degradation performance indicator and corresponding failure lifetime for five of them are shown in Figure 2.

Data Generation
The degradation of two performance indicators is described by the IG process with parameters 1 3 . The dependence structure is described by the Frank copula, and the dependent parameter 12   . A total of 21 simulation samples with random degradation data and failure lifetime data fit for the proposed model are generated by the Monte Carlo method. The degradation path for each degradation performance indicator and corresponding failure lifetime for five of them are shown in Figure 2.
(a) (b) The random data generation procedure is as follows: Step 1. Set the initial degradation value: Step 2. Generate the degradation increment for each performance indicator: (i) Generate the random number 1 (iv) The degradation increment can be obtained by through the inverse CDF function.
Step 3. Compute the degradation data and failure time data.
For comparison, we list five models in our study. These models have been investigated in previous studies and the results turned out to be very good in some applications [31,42,48]. Model 5 (M5) is the one proposed in this study, as discussed before, and the The random data generation procedure is as follows: Step 1. Set the initial degradation value: X 1 (0) = 0 and X 2 (0) = 0.
(iv) The degradation increment can be obtained by F 1 ∆X 1 t ij and F 2 ∆X 2 t ij through the inverse CDF function.
Step 3. Compute the degradation data and failure time data.
(i) The degradation data X k t ij = X k t i,j−1 + ∆X k t ij (k = 1, 2; i = 1, · · · , 21; j = 1, · · · , M; M = 20). If X k t ij ≥ d k , set the failure time T ik and terminate the simulation process. (ii) The failure time of each test sample: For comparison, we list five models in our study. These models have been investigated in previous studies and the results turned out to be very good in some applications [31,42,48]. Model 5 (M5) is the one proposed in this study, as discussed before, and the results obtained with the other four models, from Model 1 (M1) to Model 4 (M4), are compared with those obtained with M5.

•
Model 1 (M1): Only considers PI1 and unknown parameters are estimated using the degradation data (The likelihood function can be expressed as Equation (21)). • Model 2 (M2): Only considers PI1 and unknown parameters are estimated using both the degradation data and failure lifetime data (The likelihood function can be expressed as Equation (25) and the zeros-ones trick is applied to transform the likelihood function).
• Model 3 (M3): Considers two independent performance indicators and unknown parameters are estimated using the degradation data (The likelihood function can be expressed as Equation (23)). • Model 4 (M4): Considers two dependent performance indicators and unknown parameters are estimated using the degradation data (The likelihood function can be expressed as Equation (24)). • Model 5 (M5): Considers two dependent performance indicators and unknown parameters are estimated using both the degradation data and failure lifetime data (The likelihood function can be expressed as Equation (28) and the zeros-ones trick is applied to transform the likelihood function).

Parameter Estimation
The parameter estimation process is based on the discussion in Section 3.4. The unknown parameters are estimated by the Bayesian MCMC method using the OpenBUGS software. The MLE method results are also used to determine the prior distributions in the Bayesian MCMC procedure. The first 5000 samples are burn in samples. The Brooks-Gleman-Rubin (BGR) ratio is used to assess the convergence of Markov chains. In order to perform the subsequent reliability estimation and RUL estimation, 50,000 converged samples are obtained from the posterior distribution of each unknown parameter in different models. The converged samples are used to determine the mean value as well as the 2.5% and 97.5% quantile value for each parameter.
The AIC and BIC are used to choose the best copula function to describe the dependence structures of the two performance indicators. The AIC and BIC values, which are shown in Table 2, reveal that the Frank copula has the lowest AIC and BIC values, indicating that the Frank copula is the best copula function to describe the dependence structure, and M5 could be the best model to fit the simulated data. This is consistent with what we initially assumed. Thus, M5 is used to conduct the reliability and RUL estimation in the following sections. The parameter estimation results are shown in Table 3. As for M4 and M5, Table 3 only shows the estimation results of the parameter in the Frank copula. The BGR convergence results and posterior distribution for parameters in M5 are shown in Figures 3 and 4, which reveal that the MCMC became stable and converged after 5426 iterations.

Reliability Estimaton and RUL Prediction
In Section 4.2, the unknown parameters in five models (M1~M5) were obtained. The reliability estimation and RUL prediction are performed based on our discussion in Sections 2.1 and 2.2. The reliability plot for each model and the simulation model (Ms) displayed in Figure 5 reveal obvious differences in the reliability curve among each model. Only the curve of M5 is lower than that of the Ms, the curves of the other models are all higher than that of the Ms. M1 is the one with the largest deviations from the Ms. Choosing different models to fit the degradation data or failure lifetime data can produce different reliability estimation results. Aerospace 2022, 9,

Reliability Estimaton and RUL Prediction
In Section 4.2, the unknown parameters in five models (M1~M5) were obtained. The reliability estimation and RUL prediction are performed based on our discussion in Sections 2.2 and 2.3. The reliability plot for each model and the simulation model (Ms) displayed in Figure 5 reveal obvious differences in the reliability curve among each model. Only the curve of M5 is lower than that of the Ms, the curves of the other models are all higher than that of the Ms. M1 is the one with the largest deviations from the Ms. Choosing different models to fit the degradation data or failure lifetime data can produce different reliability estimation results.
The 95% confidence intervals for M4 and M5, shown in Figure 6, reveal that the confidence interval of M5 is narrower than that of M4, which indicates that M5 has higher prediction precision of the reliability estimation than M4. Therefore, models that consider two types of reliability information in the reliability estimation can achieve higher confidence compared with the models that only consider degradation data.

Reliability Estimaton and RUL Prediction
In Section 4.2, the unknown parameters in five models (M1~M5) were obtained. The reliability estimation and RUL prediction are performed based on our discussion in Sections 2.2 and 2.3. The reliability plot for each model and the simulation model (Ms) displayed in Figure 5 reveal obvious differences in the reliability curve among each model. Only the curve of M5 is lower than that of the Ms, the curves of the other models are all higher than that of the Ms. M1 is the one with the largest deviations from the Ms. Choosing different models to fit the degradation data or failure lifetime data can produce different reliability estimation results.
The 95% confidence intervals for M4 and M5, shown in Figure 6, reveal that the confidence interval of M5 is narrower than that of M4, which indicates that M5 has higher prediction precision of the reliability estimation than M4. Therefore, models that consider two types of reliability information in the reliability estimation can achieve higher confidence compared with the models that only consider degradation data.  The 95% confidence intervals for M4 and M5, shown in Figure 6, reveal that the confidence interval of M5 is narrower than that of M4, which indicates that M5 has higher prediction precision of the reliability estimation than M4. Therefore, models that consider two types of reliability information in the reliability estimation can achieve higher confidence compared with the models that only consider degradation data.
The kernel failure PDFs for M3, M4 and M5 at the time t = 1.4 are shown in Figure 7. The true failure value for the Ms is shown by the red vertical line. It can be seen that the PDF curve of M4 is close to the true value, and the PDF curves of M3 and M4 both deviate from the true value, which demonstrates that M5 has higher prediction accuracy than the other bivariate models. Aerospace 2022, 9, x FOR PEER REVIEW 17 of 27 The kernel failure PDFs for M3, M4 and M5 at the time 1.4 t  are shown in Figure  7. The true failure value for the Ms is shown by the red vertical line. It can be seen that the PDF curve of M4 is close to the true value, and the PDF curves of M3 and M4 both deviate from the true value, which demonstrates that M5 has higher prediction accuracy than the other bivariate models.  Table  4, reveal that M5 has the lowest relative errors compared to other models.
Based on the above discussion, it can be concluded that considering two types of reliability data, namely degradation data and failure lifetime data, in the reliability estimation can improve the accuracy of estimation. Furthermore, it is necessary to consider two performance indicators in the reliability estimation process.  Figure 7. The true failure value for the Ms is shown by the red vertical line. It can be seen that the PDF curve of M4 is close to the true value, and the PDF curves of M3 and M4 both deviate from the true value, which demonstrates that M5 has higher prediction accuracy than the other bivariate models. The kernel failure PDFs for M3, M4 and M5 at the time 1.4 t  are shown in Figure  7. The true failure value for the Ms is shown by the red vertical line. It can be seen that the PDF curve of M4 is close to the true value, and the PDF curves of M3 and M4 both deviate from the true value, which demonstrates that M5 has higher prediction accuracy than the other bivariate models.  Table  4, reveal that M5 has the lowest relative errors compared to other models.
Based on the above discussion, it can be concluded that considering two types of reliability data, namely degradation data and failure lifetime data, in the reliability estimation can improve the accuracy of estimation. Furthermore, it is necessary to consider two performance indicators in the reliability estimation process.  Table 4, reveal that M5 has the lowest relative errors compared to other models. Based on the above discussion, it can be concluded that considering two types of reliability data, namely degradation data and failure lifetime data, in the reliability estimation can improve the accuracy of estimation. Furthermore, it is necessary to consider two performance indicators in the reliability estimation process.
We used the leave-one-out method to verify the model proposed in this study, which has been applied in Guo et al. [32,33]. Twenty samples were used to perform parameter estimation and reliability estimation and the last sample was used to perform the RUL estimation and verify the estimation results. The RUL PDF curves at different observation time points for M4 and M5 are shown in Figure 8. For the purpose of comparing the differences between two models, the RUL PDF curves at t = 1.4 for M4 and M5 are shown in Figure 9, where the red solid vertical line is the true RUL at t = 1.4. These curves indicate that the PDF of M5 is around the true value of the RUL, which also indicates that the prediction accuracy of M5 is high. Thus, considering both the degradation data and failure lifetime data can lead to a more accurate prediction of the RUL.  We used the leave-one-out method to verify the model proposed in this study, which has been applied in Guo et al. [32,33]. Twenty samples were used to perform parameter estimation and reliability estimation and the last sample was used to perform the RUL estimation and verify the estimation results. The RUL PDF curves at different observation time points for M4 and M5 are shown in Figure 8. For the purpose of comparing the differences between two models, the RUL PDF curves at 1.4 t  for M4 and M5 are shown in Figure 9, where the red solid vertical line is the true RUL at 1.4 t  . These curves indicate that the PDF of M5 is around the true value of the RUL, which also indicates that the prediction accuracy of M5 is high. Thus, considering both the degradation data and failure lifetime data can lead to a more accurate prediction of the RUL.

Introduction to Mechanical Seal
In this section, the proposed model is demonstrated using a real dataset with degradation data and failure lifetime data for mechanical seals. Mechanical seals can seal the pump shaft and they also allow the pump shaft to rotate, contain pressure, prevent leak-

Introduction to Mechanical Seal
In this section, the proposed model is demonstrated using a real dataset with degradation data and failure lifetime data for mechanical seals. Mechanical seals can seal the pump shaft and they also allow the pump shaft to rotate, contain pressure, prevent leakage and exclude contamination. A schematic of mechanical seals is shown in Figure 10. Mechanical seals are widely used in the aviation industry. It is one of the most important mechanical components in an aviation pump system. Two performance indicators, namely friction torque and leakage rate, are commonly used to evaluate the performance of mechanical seals [49,50]. The distribution of the film between the rotating ring and static ring can be affected by the loading pressure or thermal effects [51]. It can also result in the performance degradation of the two performance indicators. Therefore, it is important to consider the bivariate-dependent relationship between friction torque and leakage rate when conducting reliability estimation.

Introduction to Mechanical Seal
In this section, the proposed model is demonstrated using a real dataset with degradation data and failure lifetime data for mechanical seals. Mechanical seals can seal the pump shaft and they also allow the pump shaft to rotate, contain pressure, prevent leakage and exclude contamination. A schematic of mechanical seals is shown in Figure 10. Mechanical seals are widely used in the aviation industry. It is one of the most important mechanical components in an aviation pump system. Two performance indicators, namely friction torque and leakage rate, are commonly used to evaluate the performance of mechanical seals [49,50]. The distribution of the film between the rotating ring and static ring can be affected by the loading pressure or thermal effects [51]. It can also result in the performance degradation of the two performance indicators. Therefore, it is important to consider the bivariate-dependent relationship between friction torque and leakage rate when conducting reliability estimation.  The degradation test of the mechanical seals was conducted to monitor the two performance indicators of the mechanical seals, namely friction torque (PI1) and leakage rate (PI2), as shown in Figure 11. The operating conditions of the test are listed in Table 5. The degradation test of the mechanical seals was conducted to monitor the two performance indicators of the mechanical seals, namely friction torque (PI1) and leakage rate (PI2), as shown in Figure 11. The operating conditions of the test are listed in Table 5.

Operating Conditions
Value Rotation speed 5400 r/min Temperature Room temperature (23 °C) Pressure in cylinder chamber 0.15 MPa Figure 11. Degradation test rig used for the degradation test of mechanical seals.

Model Application
The degradation paths of these two performance indicators are shown in Figure 12. The degradation data from 0 to 450 h are not considered due to almost negligible leakage rate. The degradation increment data are transformed into absolute values. Before performing the unknown parameters estimation, the degradation data need to be normalized. The time interval when conducting the reliability estimation is set to be 0.02 t   . Different time interval settings will affect the estimation of unknown parameters. However, it will not affect the reliability, MTTF and RUL estimation results. Similar to the simulation study, we also use six models that are commonly used in other studies to compare

Model Application
The degradation paths of these two performance indicators are shown in Figure 12. The degradation data from 0 to 450 h are not considered due to almost negligible leakage rate. The degradation increment data are transformed into absolute values. Before performing the unknown parameters estimation, the degradation data need to be normalized. The time interval when conducting the reliability estimation is set to be ∆t = 0.02. Different time interval settings will affect the estimation of unknown parameters. However, it will not affect the reliability, MTTF and RUL estimation results. Similar to the simulation study, we also use six models that are commonly used in other studies to compare the goodness of fit [31,42,48], as shown below: Unknown parameters in different models are estimated according to the discussion in Section III. The parameter estimation results for the different models are shown in Table  6. These results reveal that the values of the degradation rate  , time scale function q and dependency parameter  are affected when the failure lifetime data are considered. The degradation rate  affects the velocity of the degradation process and q in the time scale function affects the shape of the degradation path. The dependency parameter  affects the dependent strength between friction torque and leakage rate. All these parameters have an influence on the accuracy of the reliability estimation for a mechanical seal. The AIC and BIC values of different copulas describing the bivariate-dependent structures between two performance indicators are shown in Table 7. These results reveal that Frank copula has the lowest AIC and BIC values, which indicates that the Frank cop- The six test seals' degradation paths of two performance indicators and corresponding failure lifetime are shown in Figure 12.
Unknown parameters in different models are estimated according to the discussion in Section III. The parameter estimation results for the different models are shown in Table 6. These results reveal that the values of the degradation rate λ, time scale function q and dependency parameter θ are affected when the failure lifetime data are considered. The degradation rate λ affects the velocity of the degradation process and q in the time scale function affects the shape of the degradation path. The dependency parameter θ affects the dependent strength between friction torque and leakage rate. All these parameters have an influence on the accuracy of the reliability estimation for a mechanical seal. The AIC and BIC values of different copulas describing the bivariate-dependent structures between two performance indicators are shown in Table 7. These results reveal that Frank copula has the lowest AIC and BIC values, which indicates that the Frank copula is the best to describe the bivariate-dependent relationship between the leakage rate and friction torque for mechanical seals. Based on the above estimations, the reliability and RUL estimations can be obtained. The reliability curves for the six different models in Figure 13a-c show the partial magnification of the curves in the periods of 775-950 h and 900-1050 h, respectively. These results show that the dark red curve for M1 and green curve for M2 are higher than those for the other models, which indicates that only considering one performance indicator results in estimation deviations in the reliability estimation. The purple dashed curve for M4 is initially higher than the blue dashed curve for M5, and then decreases, which indicates that the bivariate-dependent relationship also affects the reliability estimation results. The blue dashed curve for M5 is higher than the pink curve for M6, which indicates that considering the failure lifetime data affects the accuracy of the reliability estimation. The corresponding failure probability density of M5 and M6, which are also shown in Figure 14, reveal that the PDF curve of M6 is much slimmer than that of M5, and this means that M6 is more precise than M5. The MTTF estimation results for these six models are shown in Table 8, where the differences among those models are obvious. The MTTF given by the manufacturer is 944 h. M6 has the lowest estimation deviation. It also indicates that using both two types of reliability information can improve the reliability estimation accuracy.    (b) (c)    Similarly, we also used the hold out method to compare these models mentioned above. The eleven seal samples under examination were used to conduct parameter estimation and reliability estimation, and the last one was used to verify the effectiveness of the proposed model. The RUL estimation values based on M5 and M6 for the twelfth test seal are shown in Figure 15, where the red line represents the actual RUL value and the blue dashed line and red line represent the PDF of the RUL for M5 and M6 at different measurement times, respectively. We can also find that the PDF curve of M6 is much slimmer than that of M5, and this means that M6 is more precise than M5. Therefore, it is necessary to fuse the degradation data and failure lifetime data when conducting reliability estimation. Similarly, we also used the hold out method to compare these models mentioned above. The eleven seal samples under examination were used to conduct parameter estimation and reliability estimation, and the last one was used to verify the effectiveness of the proposed model. The RUL estimation values based on M5 and M6 for the twelfth test seal are shown in Figure 15, where the red line represents the actual RUL value and the blue dashed line and red line represent the PDF of the RUL for M5 and M6 at different measurement times, respectively. We can also find that the PDF curve of M6 is much slimmer than that of M5, and this means that M6 is more precise than M5. Therefore, it is necessary to fuse the degradation data and failure lifetime data when conducting reliability estimation.

Conclusions and Future Work
This study evaluated the reliability estimation model based on the IG process model and bivariate dependence analysis by copula functions, fusing two sources of reliability information, namely the degradation data and failure lifetime data. First, the IG process model for each performance indicator was established to describe the degradation process. The bivariate-dependent model was developed using copula functions to describe the dependent structures between two performance indicators. Then, the likelihood functions were constructed considering different types of reliability information, including only degradation data, only failure lifetime data and both degradation data and failure lifetime data. The MLE and Bayesian MCMC methods were used to estimate the unknown

Conclusions and Future Work
This study evaluated the reliability estimation model based on the IG process model and bivariate dependence analysis by copula functions, fusing two sources of reliability information, namely the degradation data and failure lifetime data. First, the IG process model for each performance indicator was established to describe the degradation process. The bivariate-dependent model was developed using copula functions to describe the dependent structures between two performance indicators. Then, the likelihood functions were constructed considering different types of reliability information, including only degradation data, only failure lifetime data and both degradation data and failure lifetime data. The MLE and Bayesian MCMC methods were used to estimate the unknown parameters in the degradation model considering the various types of likelihood functions. Ultimately, the reliability and RUL values for the mechanical components were obtained from the unknown parameter estimation results and Monte Carlo method. Comparative analysis of the simulation study and real case study showed the effectiveness of the proposed model.
The main contributions provided by this study are as follows: (1) The development of a reliability estimation model considering dependent degradation performance indicators using the IG process model and copula functions. (2) To fully use the reliability information obtained, the likelihood functions are constructed using the degradation data and failure lifetime data based on the developed reliability model. Since the likelihood function contains information from different sources, unknown parameters in the model are estimated by transforming the likelihood function using the zeros-ones trick.
Future research will be focused on the improvement of parameters for real-time updating based on the offline reliability information and online new reliability information.
Combining the physical of failure model and data-driven model in dependent analysis with multiple performance indicators will also be considered in future studies.