Claudio M Rocco

Enrico Zio

Global sensitivity analysis in a Multi-state Physics Model of Component Degradation based on a hybrid State-Space Enrichment and Polynomial Chaos Expansion approach

Keywords: Degradation, Global Sensitivity Analysis, Inhomogeneous Continuous Time Markov Chain, Multi-state physics model, Parameter Uncertainty, Polynomial Chaos Expansion, State-space Enrichment

This paper extends previous works related to the assessment of component degradation, through Markov multi-state physic models. The extension includes the evaluation of the effects of uncertain parameters in the model and the definition of their importance with respect to their influence on the output of the model. Global Sensitivity Analysis (GSA) is selected as the technique because able to: 1) consider the simultaneous effects of parameters variations and 2) to define importance indexes that allow a ranking of the components. GSA requires a large number of evaluations for specific points, identified by an appropriate design of experiment. To avoid the many costly evaluations, a meta-model is built based on polynomial chaos expansion (PCE). A PCE is a multidimensional polynomial approximation of the model with coefficients determined by evaluating the model in a reduced set of predetermined points. Importance indexes values are then derived directly from the PCE. Since in the problem considered, the model provides the timedependent behavior of the state probabilities, the importance indexes are also functions of time. An application is presented, related to the cracking process in an Alloy 82/182 dissimilar metal weld in the primary coolant system of a nuclear power plant.

I. INTRODUCTION

Multi-state models (MSM) [START_REF] Kuo | Optimal Reliability Modeling: Principles and Applications[END_REF]- [START_REF] Levitin | Gregory Optimal Structure of Multi-State Systems With Uncovered Failures[END_REF] have been proposed to assess the degradation processes of components or systems [START_REF] Park | Stochastic degradation models with several accelerating variables[END_REF]- [START_REF] Gebraeel | Residual Life Predictions in the Absence of Prior Degradation Knowledge[END_REF]. In general, a state-space diagram is used to define the stochastic process of transition among system states which represent a range of levels from perfect functioning to complete failure.

Then, a model is developed and solved to represent the dynamical behavior of the probability of being in each state [START_REF] Billinton | Reliability Evaluation of Engineering Systems[END_REF]- [START_REF] Zio | Computational Methods For Reliability And Risk Analysis[END_REF]. In general, transition rates are estimated from field data, but "in practice, it can be difficult or even impossible to collect relevant data especially for the highly reliable devices" [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF].

Recently, Unwin et al. [START_REF] Unwin | Multistate physics models of aging passive components in probabilistic risk assessment[END_REF] proposed a multi-state physics modeling approach to assess the degradation processes (e.g. crack growth) of components or systems. The main idea is to describe the transition rates among states by physic functions rather than estimated from service data. As a result transition rates are time-dependent thus defining a non-Markovian model. Several approaches have been devised to solve the non-Markovian problem, such as state-space enrichment approach [START_REF] Unwin | Multistate physics models of aging passive components in probabilistic risk assessment[END_REF], stochastic Petri net (SPN) [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF]- [START_REF] Schneeweiss | Tutorial: Petri nets as a graphical description medium for many reliability scenarios[END_REF] and Monte Carlo simulation [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF].

In general, the physic functions used to describe transition rates depend on a set of parameters or factors. For example the micro-crack initiation transition rate  1 , of the application considered in Section 5, from initial state S to micro-crack state M is defined as:

   1 1   b t b   
, where  and b are fitting parameters estimated from observed data, and possibly subject to uncertainties.

The existence of uncertainties poses the typical question: What happens to the probability of a specific state if parameters have uncertainty? Moreover, which parameters affect the output of the model the most? In general, these questions are often addressed by analyzing the effect of the variation of a single parameter at a time (also known as Once At a Time evaluation) or using Taylor-series based expansions [START_REF] Havekort | Sensitivity and uncertainty analysis of Markovreward models[END_REF]. However, to account for the possible effects of interactions among parameters (for example between b and  in  1 ), global sensitivity analysis (GSA) methods are required. GSA is known to consider the effect of a factor while all others are varying as well, thus allowing the exploration of the multi-dimensional input space [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF] and is able to produce indexes that assess the importance of the factors, taking into account interactions among them.

In general, importance indexes are calculated by evaluating the model output (i.e. the state probabilities for the case of interest here) at a set of specific points derived using random or special sampling techniques from the joint distribution of the space of input factors [START_REF] Saltelli | Sensitivity Analysis[END_REF].

GSA methods based on variance decomposition decompose the output variance in contributions of each factor (main effects) and possible interactions. GSA has been used in many applications (e.g., [START_REF] Tarantola | Sensitivity analysis using contribution to sample variance plot: Application to a water hammer model[END_REF]- [START_REF] Xu | The effect of parameter uncertainty on achieved safety, integrity of safety system[END_REF]). Recently, Rocco and Zio [START_REF] Rocco | Global sensitivity analysis in Markov reliability models of power system components[END_REF] applied GSA to assess the importance of transition rates, in homogeneous continuous Markov models of power system components.

The number of model evaluations required for a GSA depends on the number of factors in the model and the sample size. The evaluation of the indexes depends on the complexity of the model and could require very large computational times. In these cases, it is possible to use surrogate or meta-models that approximate the original model performance. Several techniques have been suggested to build such metamodels, like regression models, artificial neural networks or support vector machines, among others.

In this paper, a particular technique named Polynomial Chaos Expansion (PCE) is selected to create a meta-model represented by a polynomial of the factors. If the factor uncertainties are modeled as random variables with known probability density functions, then the coefficients of the resulting polynomial contain "the complete response of the model" [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]. In this way, factor importance indexes can be assessed directly from the coefficients of the polynomial decomposition.

PCE has been used in many areas including GSA [START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF]- [START_REF] Tootkaboni | Topology optimization of continuum structures under uncertainty -A Polynomial Chaos approach[END_REF]. To our best knowledge it has not been used to assess importance indexes in multi-state physics models of component degradation.

In this paper, we extend the use of GSA in Markov models [START_REF] Rocco | Global sensitivity analysis in Markov reliability models of power system components[END_REF] to multi-state physics models whose solutions are not known analytically and thus obtained by numerical methods.

Note that multi-state physics models assess the behavior of the system under analysis as a function of time. For this reason, in this paper, the PCE of a dynamic system, recently suggested in [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF] is selected. As a result, a PCE is used at each time in order to derive the dynamical behavior of the importance indexes of the factors.

The presentation of the work is structured in the paper as follows: Section II defines the Markov model setting. Section III briefly reviews some SA concepts, while

Section IV provides a brief introduction to PCE. Section V presents results of experimentation and Section 6 presents conclusions and future work.

Acronyms, Notations and Assumptions

Acronyms

II. MARKOV MODELING

The evolution in time of a system which changes state stochastically as its components change modes of operation at random times (due to degradation, reduced production, failure, repair, etc.) can be mathematically described by a discrete-state, continuous-time Markov chain (CTMC) Z={z(t), t  0}, with finite state space E={1,2,…, m}.

Let Q(t) be the m x m infinitesimal generator matrix whose generic element q ij (t)

represents the transition rate from state i to state j (i,jE), q ii (t)=  ij q ij (t), and P i (t) be the probability that the system is in state i at time t. The m-dimensional probability vector P(t) is obtained by solving:

) ( ) ( ) ( t Q t P t P   (1)
given initial conditions for P(0).

In general, transition rates can be constant (homogenous Markov Chain) or functions of time (non-homogeneous Markov Chain). General techniques for solving Ordinary Differential Equations, like Euler or Runge-Kutta methods [START_REF] Tombuyses | Solving Markovian systems of O.D.E. for availability and reliability calculations[END_REF], the special implicit method [START_REF] Hosea | Analysis and implementation of TR-BDF2[END_REF] in case of stiff models or uniformization-based algorithms can be used for solving (1) [START_REF] Van Moorsely | Numerical Solution of Non-Homogeneous Markov Processes through Uniformization[END_REF].

In some models, as the one of Section 5, some transition rates are functions of the holding times in each state. In these cases, the previous general methods can not be used directly. As suggested in [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF]- [START_REF] Unwin | Multistate physics models of aging passive components in probabilistic risk assessment[END_REF] one of the possible techniques that can be used is the method of state enrichment. Basically the method is based on the "discretization of the component lifetime into equally sized time intervals, during each of which the transition rate remains constant". The new enriched state-space is described by the original states and the vector of discretized times at each state. For example, given a component with 6 original states and a lifetime of 100 years discretized in time steps of 0.5 years, the approximation method would generate a transition matrix with (6x(100/0.5)) 2 entries.

In any case, no matter the solution technique used for solving (1), the numerical values of P i (t) are non-linear functions of the transition rates q ij (t) with possible interactions among them. To evaluate the influence of transition rate uncertainties on P i (t) would mean to consider the variability in their values and the corresponding variability in the model state probabilities output. It is also of interest to identify, through a sensitivity analysis, those transition rates most influential in terms of the variations that they induce on the state probabilities when their variability is considered.

III. SENSITIVITY ANALYSIS

Sensitivity analysis (SA) is "the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input factors" X=(x 1 , x 2 , …., x n ), n being the number of factors in the model [START_REF] Saltelli | Sensitivity Analysis[END_REF].

Global sensitivity analysis (GSA) considers the exploration of the multidimensional input space [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF] and is able to determine the effect of a factor while all others are varying as well. Some of the methods developed, such as Fourier Amplitude Sensitivity Technique (FAST), extended FAST [START_REF] Saltelli | Sensitivity Analysis[END_REF] provide a numerical assessment of the importance of each the factor and are able to consider and quantify the interaction among them.

In this paper, variance-based sensitivity indexes are evaluated. The main idea is to decompose the variance of an output as a function of the main effects of each factor and possible interactions.

Let Y = f(x 1 , x 2 , ..., x n ), be a scalar output and x i the n orthogonal input factors treated as random variables (i.e. uncertainty is characterized in terms of a probability density function (pdf)). The variance V(Y) of the output Y can be written as [START_REF] Saltelli | Sensitivity Analysis[END_REF]:

          i i j j l n ijl i i j ij i i V V V V Y V ... 12 .... ) ( (2)
where:

  i i X Y E V V   is the main effect (or first-order term) due to x i ,     j i j i ij V V X X Y E V V   
, is the two-way interaction between x i and x j , and so on.

Main and total order sensitivity indexes can then be defined as [START_REF] Saltelli | Sensitivity Analysis[END_REF], [START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF]:

    Y i i i i V X Y E V S   (3)         Y i i i Y i i i Y Ti V Y V E V X Y E V V S       X (4)
where X -i = {x 1 , x 2 ,…x i-1 , x i+1 , ...x n ) and E -i (Y│x i ) is the expected value of Y conditioned on x i , thus being a function of factor x i alone. The final expression in ( 4)

is derived from the fact that, given two random variables A and B:

      B A E V B A E A V   ) ( [14].
The main index S i is the fraction of variance of the output V(Y) that can be attributed to x i alone, while S Ti corresponds to the fraction of V(Y) that can be attributed to x i including all its interactions with other factors [START_REF] Saltelli | Sensitivity Analysis[END_REF]. For example, in the case of 3 factors x 1 , x 2 and x 3 : S T1 = S 1 + S 12 + S [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF] .

That means that "the computation of total sensitivity indices is a shortcut to measure the overall effect of "main effect" and all higher order interactions" [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF].

The main index S i is the measure to be used to determine which subset of parameters accounts for most of the output uncertainty [START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF], while S Ti must be used as a screening tool for identifying the subset of non-influential factors, i.e. factors that when fixed to any value within their uncertainty range, do not significantly reduce the output variance [START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF]. The difference (S Ti -S i ) assesses the interaction effects of a factor.

The estimation of S i and S Ti are approximated by assuming independence among factors and using random or special sampling techniques (e.g., Monte Carlo simulation based on Latin Hypercube Sampling (LHS) approach from the joint distribution of the space of input factors). Techniques such as Sobol or extended FAST could be used to determine S i and S Ti [START_REF] Saltelli | Sensitivity Analysis[END_REF]. One of the main disadvantages of this approach is that it requires several model evaluations and, thus, can be computationally time consuming. For example, the "extended-FAST" method [START_REF] Saltelli | A quantitative, model independent method for global sensitivity analysis of model output[END_REF] used to estimate first-order and total indices for all the n factors of a model, has a total computational cost of (ns x n) model evaluations, where ns is an integer (that can vary "from a few hundreds to a few thousands" [START_REF] Saltelli | Global Sensitivity Analysis: The Primer[END_REF]) representing the sample size selected and depends on the convergence properties of the sensitivity indexes. Thus, depending on the number of factors and the computational complexity of the model, the technique could require very large computational times.

IV. POLYNOMIAL CHAOS EXPANSION

A. Introduction

The basic idea of PCE is to approximate the outputs of a system (e.g., state probabilities in our case) through an orthonormal polynomial basis in the uncertain parameter spaces (e.g., transition rates). Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] was the first to propose the use of Polynomial Chaos (PC), when parameters follow a Gaussian distribution. In this case, the basis is defined through Hermite polynomials. However, this limitation was extended to generalized polynomial chaos (gPC), which is able to consider more distributions [START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic Differential equations[END_REF]. Recently [START_REF] Oladyshkin | A concept for data-driven uncertainty quantification and its application to carbon dioxide storage in geological formations[END_REF] defines the arbitrary polynomial chaos (aPC), which generalizes PC and is able to accommodate arbitrary distributions or distributions known with a few statistical moments.

Let =[ 1 ,…,  n ] define the vector of n random input factors and ()=f() a secondorder with finite variance model output. PC theory [START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic Differential equations[END_REF] shows that the spectral expansion the model output () can be expanded as:

()= (5)     0 1 ) ( i n i i ,..., ξ ξ a
where: a i are the coefficients of the expansion  i are random variables with a pdf defined  i () are polynomials of the selected basis.

The polynomial basis is selected based on the pdf of  i . As an example, if  i are standard normal random variables, then the basis is defined as multivariate Hermite polynomials. Table 1 shows the type of distribution and the associated polynomials as suggested in the Wiener-Askey scheme [START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic Differential equations[END_REF].

Table 1: Relation between probability distribution and polynomial basis based on the Wiener-Askey scheme [42]

For example, if  is a standard normal random variable, the Hermite polynomials up to third order are:

k H k () 0 1 1  2 -1+ 2 3 -3+ 3
For practical purpose, the infinite series in ( 2) is limited to a finite number of terms M+1 which depends on n, the number of random variables of the model and on the selected polynomial degree d: [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]. For example, if n=2 and d=3, then M+1=10.

! ! )! ( 1 d n d n M   
The multivariate polynomial  i () is given by the tensor product of the corresponding one dimensional selected polynomial [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]. For example, if  1 and  2 are standard normal random variables, the 10 terms of the PCE expansion for ( 1 , 2 ) up to third order is given by [START_REF] Liang | An EXCEL add-in implementation for collocation-based stochastic response surface method[END_REF]:

) ( ) ( ) 3 ( ) 3 ( ) 1 ( ) 1 ( ) , ( 2 2 1 2 9 1 2 2 1 8 2 3 2 7 1 3 1 6 2 1 5 2 2 4 2 1 3 2 2 1 1 0 2 1                                    a a a a a a a a a a
The terms a i are unknown and must be determined. There are two approaches for their evaluation: the first is based on the modification of the model () and its solution (intrusive approach) while the second does not require any modification of the model equations (non-intrusive approach). In the latter approach, the terms are estimated from the evaluation of the model under study in a set of selected samples of the random vector . In our case, the set of differential equations ( 1) is solved, based on a set of transition rate deviates.

Two techniques are described in the literature for the non-intrusive approach. The first (regression method) is based on the optimization of the mean-square difference between the model and the approximation [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. The optimization is performed on a set of N=2(M+1) samples  j and its corresponding evaluation ( j ). The second technique is based on projection methods which exploit the orthogonality properties of the PC basis. In fact, by evaluating the inner product of the PC expansion (equation ( 2)) with  l ():

(),  l ()  ) ( ), ( 1 1 0 n l n i M i i ,..., ξ ξ ,..., ξ ξ a     (6)
where , is the inner product in the Hilbert space, and taking into account the orthogonality of the basis, the coefficients can be obtained as [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]:

) ( ), ( ) ( ), ( 1 1 1 1 n i n i n i n i ,..., ξ ξ ,..., ξ ξ ,..., ξ ξ ,..., ξ ξ a      (7)
The numerator in (4) is computed as:

      d n n n i n n i n ,..., ξ ξ d ,..., ξ ξ h ,..., ξ ξ ,..., ξ ξ ,..., ξ ξ ,..., ξ ξ ) ( ) ( ) ( ) ( ) ( ), (
As mentioned in [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF], the coefficients could be obtained by solving up to (M+1) multidimensional integrals [START_REF] Park | Stochastic degradation models with several accelerating variables[END_REF]. However, the integrals are numerically approximated by evaluating ( j ) in a set of selected ni integration points  j (j=1,..,ni), using the following general expression [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]:

) (9) ( 1 ) ( ) ( 1 ) ( ) ( 1 1 1 1 1 ) ( ) ( ) ( ) ( ) ( ) ( j ni j j n j i j n j d n n n i n w ,..., ξ ξ ,..., ξ ξ ,..., ξ ξ d ,..., ξ ξ h ,..., ξ ξ ,..., ξ ξ        
where w (j) are specific weights. This means that the computational cost for determining the polynomial coefficients are proportional to the number of integration points ni.

There are two approaches to solve the multidimensional integrals: Deterministic and Probabilistic algorithms [START_REF] Panayirci | On the Capabilities of the Polynomial Chaos Expansion Method within SFE Analysis-An Overview[END_REF]. The first ones are based on high order interpolation schemes. For example, the Smolyak algorithm [START_REF] Petras | Smolyak cubature of given polynomial degree with few nodes for increasing dimension[END_REF] of order oq uses an interpolation approach by evaluating polynomials from lower dimensions to higher order at selected points. Figure 1 shows the set of points generated by the algorithm for d = 2 and oq=4 (33 points) and oq=6 (65 points). Note that the set of points for oq =4 is included in the set of points for oq =6.

The second, like Monte Carlo methods, uses random samples. In both approaches, the convergence of the integrals depends on the smoothness of the model under analysis [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF].

As suggested in [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF], the Smolyak approach is, in general, a convenient method when the number of factors is less than 15, while Monte Carlo methods are "best suited for higher dimensional problems". In this paper, the Petras-technique [START_REF] Petras | Smolyak cubature of given polynomial degree with few nodes for increasing dimension[END_REF], a special Smolyak-based method is used, as implemented in the Scilab toolbox NISP [START_REF] Baudin | Polynômes de chaos sous Scilab via librairie NISP[END_REF].

After determining the coefficients, the polynomial meta-model is already defined. As shown in [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] and [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] sensitivity indexes for all factors, including S i and S Ti are evaluated analytically from the polynomial coefficients. So the computational cost for deriving them corresponds basically to the cost required to evaluate numerically a set of multidimensional integrals.

B. PCE for dynamical systems

The PCE of dynamical systems has been recently proposed in [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]. Since the output of the model is a function of time, the PCE is based on coefficients that also vary as a function of time. This means that at each time instant a set of parameters are required.

If =f(,t) then [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF]:

(,t)   (10)   M i i i t a 0 n 1 ) ξ ,..., ξ ( ) (
As in the previous subsection, each time the variance decomposition is performed analytically from the polynomial coefficients, all sensitivity indexes are determined for each parameter in the model.

Based on the two-step procedure proposed in [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF], the following approach is defined:

1. Based on the parameter distribution, select the corresponding polynomial type 2. Define the polynomial degree d 3. Calculate M+1 4. For t=0 until a specific time limit t max { 5. Estimate the PC coefficients 6. Calculate sensitivity indexes for each parameter }

In step 2, the selection of the appropriate degree d could be done iteratively, that is, trying different values for d and evaluating the average relative error between the model output and the PCE approximation [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF].

V. EXAMPLE: A NUCLEAR COMPONENT UNDERGOING STRESS CORROSION CRACKING

This example refers to the cracking process in an Alloy 82/182 dissimilar metal weld in a primary coolant system of a nuclear power plant, previously analyzed by [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF]- [START_REF] Unwin | Multistate physics models of aging passive components in probabilistic risk assessment[END_REF] using a multi-state physics model. Figure 2 shows the six-state transition diagram model. Transition rates among states are defined as:  1 is the micro-crack initiation transition rate,  2 ,  3 and  4 are crack growth transition rates,  5 and  6 are transition rates that lead to ruptures,  1 ,  2 ,  3 and  4 are transition rates associated with detection/repair at the different states. In [START_REF] Unwin | Multistate physics models of aging passive components in probabilistic risk assessment[END_REF], the transition rates  1 ,  2  3 and  4 are time-dependent and/or stochastic, whereas the others are assumed constant.

An interesting fact in this example is that some transition rates depend on the time the system is in a particular state. For example, the transition rate  2 (between state M and D) describes the transitions from micro-crack state M to radial-crack state D and is defined as a function of u, the time of crack initiation.

The expressions for transition rates  1 ,  2  3 and  4 are presented in Table 2. As shown in Table 2, the transition rates are function of a set of 14 factors (i.e., n=14)

whose base values are shown in Table 3.

Table 2: Expressions for transitions rates  1 ,  2  3 and  4

Table 3: Parameter definition and base value

Figure 3 shows the time-dependent behavior of P i (t) for each state, computed using the state enrichment method for a component lifetime t max =80 years and a time step of 1 year [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF].

Figure 2: Transient behavior for P i (t), base case [9]

As observed in [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF], "1) there is an early, rapid transition from the Initial state to the Micro-crack state; 2) there is a monotonic increase in the probability of the Rupture state". In the time span 5-20 years it is possible to detect a transient behavior.

For the importance indexes assessment, a  10 % variation on base values is assumed described by uniform distributions. Importance sensitivity indexes are derived using the PCE approach, based on a Legendre polynomial basis as suggested from Table 1.

Note that each factor must be mapped on the interval [-1,1].

In this paper, a non-intrusive approach is selected: State probabilities P i (t) are evaluated using the state enrichment method as the "black-box" model. The PCE as well as the sensitivity indexes are evaluated using the Scilab toolbox NISP [START_REF] Baudin | Polynômes de chaos sous Scilab via librairie NISP[END_REF] and a Smolyak integration technique is selected based on Petras integration points [START_REF] Petras | Smolyak cubature of given polynomial degree with few nodes for increasing dimension[END_REF].

As suggested in [START_REF] Haro | Sensitivity study of dynamic systems using polynomial Chaos[END_REF], different values for d were selected and the average relative error between the model output and the PCE approximation were evaluated. For d=3 the PCE is considered sufficient since the maximum error is less than 6.2%, as shown in Figure 4. In this case, the number of models evaluations is equal to 3361 [START_REF] Heiss | Likelihood approximation by numerical integration on sparse grids[END_REF] and the PCE expansion has 681 terms.

Figures 567show the time-dependent behavior of S i , S Ti along with the value of (S Ti -S i ) for each factor considered and for each state probability. The title in each figure is coded as: sx PETRAS yy z-ww, where:

x is the state considered yy is the type of index presented z is the polynomial degree ww is the number of model evaluations (PETRAS stands for the option selected in NISP [START_REF] Baudin | Polynômes de chaos sous Scilab via librairie NISP[END_REF] for the determination of the ni integration points)

Let us consider the time-dependent behavior of the sensitivity indexes related to the initial state S (state 1). As mentioned in Section C, S i is used to determine which parameters account for most of the output uncertainty, S Ti is used for identifying the non-influential factors, and (S Ti -S i ) for quantifying interaction effects. To perform the analysis the following criterion was used: if S i or S Ti or (S Ti -S i ) < 0.1 then the corresponding index is not considered.

The analysis reveals: 7: The index (S Ti -S i ) reflects the interaction effects between a factor and all the others. In general the index is small (i.e., S Ti S i ) during the component lifetime. However in the time span 5-10, there are interactions between  and the rest of the parameters.

A similar analysis can performed for the other states. In general, the importance of each factor is dependent not only on the state considered but also at what time is evaluated.

The index S i  0.1 is used to determine which parameters account for most of the output uncertainty. Table 3 shows the important factors for each probability state during the time-span simulated.

CONCLUSIONS

The output of a model depends on the input parameters. In many situations, inputs are not precisely known because subject to uncertainty. It is important to know how uncertainty in the parameters affects the results of the model. Such assessment is possible using global sensitivity analysis, an approach able to analyze the variance of the output and determine the individual effects of each parameter (factor) as well as the possible interactions among them. In practice, global sensitivity analysis requires many evaluations of the model which may pose a problem of computation time. In order to reduce the computation burden, a surrogate model can be constructed. This paper presents a hybrid approach for evaluating the importance of uncertainty factors in multi-state physics models of component degradation.

The original contributions of the paper are : a) the assessment of the effects of uncertain parameters in a Markov multi-state physics model of component degradation and of their importance with respect to the model output by using GSA; and b) the use of Polynomial Chaos Expansion to derive an efficient meta-model for coping with the computational cost of the analysis.

An example is presented related to a cracking process described by a six-state Markov model with 14 uncertain factors, The meta-model is based on a PCE of third order whose coefficients are determined using only 3341 model evaluations. The metamodel is able to mimic the dynamic behavior of the degradation state probabilities, with a maximum average relative error less than 6.2 %. The sensitivity indexes for each state reveal that the ranking of importance is time-varying and each state has its own set of important factors.

Future works aim to test some approaches for: 1) reducing the number of evaluations without affecting the accuracy of the expansion, for example, through an initial screening approach, as the one suggested in [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]; 2) evaluating the effects of different probability distributions associated to the uncertain factors; 3) considering the robustness of the importance indexes when modeling dependencies among factors; and 4) The convenience of using an integrated importance index able to consider simultaneously the effects of one factor in several outputs. Figure 3: Time-dependent behavior for Pi(t), base case [START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF] Figure 4: Error between the model output and the PCE approximation, using d=3 Figure 5: Time-dependent behavior for S i , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture Figure 6: Time-dependent behavior for S Ti , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture Figure 7: Time-dependent behavior for S Ti -S i , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture

Figure 2 :

 2 Figure 2: Transition diagram of the multi-state physics model of crack

Figure 1 :Figure 2 :

 12 Figure 1: 2D Smolyak points: a) oq = 4; b) oq =6 Figure 2: Transition diagram of the multi-state physics model of crack development in Alloy 82/182 dissimilar metal welds [9]Figure3: Time-dependent behavior for Pi(t), base case[START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF] Figure4: Error between the model output and the PCE approximation, using d=3 Figure5: Time-dependent behavior for S i , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture Figure6: Time-dependent behavior for S Ti , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture Figure7: Time-dependent behavior for S Ti -S i , for states 1:initial; 2:microcrack; 3:circumferential; 4:radial; 5:leak;6:rupture

Figure 1 :Figure 2 :

 12 Figure 1: 2D Smolyak points: a) oq = 4; b) oq =6

Figure 3 :

 3 Figure3: Time-dependent behavior for P i (t), base case[START_REF] Li | A Multi-state Physics Model of Component Degradation based on Stochastic Petri Nets and Simulation[END_REF]

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: Error between the model output and the PCE approximation, using d=3

  1) Top left panel in Figures 5: During the component lifetime, parameter  accounts for most of the output uncertainty. However in the time span 8-10 parameter a C accounts for most of the output. For t > 20,  and  1 are the most important parameters. 2) Top left panel in Figures 6: Several parameters have low S Ti values, indicating that they are non-influential factors 3) Top left panel in Figure

Table 1 :

 1 Relation between probability distribution and polynomial basis based on the Wiener-Askey scheme[START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic Differential equations[END_REF]

Table 2 :

 2 Expressions for transitions rates  1 ,  2  3 and  4

Table 3 :

 3 Parameter definition and base value Table4Important Parameters using a PCE with d=3

Table 1 :

 1 Relation between probability distribution and polynomial basis based on the Wiener-Askey scheme[START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic Differential equations[END_REF]

	Distribution	Polynomials	Support
	Gaussian	Hermite	
	Uniform	Legendre	[-1,1]
	Gamma	Laguerre	[0,[
	Beta	Jacobi	[-1,1]
	Poisson	Charlier	{0,1,…}
	Negative Binomial	Miexner	{0,1,…}
	Binomial	Krawtchouk {0,1,…N}
	Hypergeometric	Hahn	{0,1,…N}

Table 2 :

 2 Expressions for transitions rates  1 ,  2 ,  3 and  4

	Transition rate	Expression

Table 3 :

 3 Parameter definition and base value

	Parameter	Definition	Base value
	b	Weibull shape parameter for crack initiation model	2.0
		Weibull scale parameter for crack initiation model	4 years
	a D	Crack length threshold for radial macro-crack	10 mm
	P D	Probability that micro-crack evolves as radial crack	0.009
	a M	Maximum credible crack growth rate	9.46 mm/yr
	a C	Crack length threshold for circumferential macro-crack	10 mm
	P C	crack Probability that micro-crack evolves as circumferential	0.001
	a L  1  2  3  4  5  6	Crack length threshold for leak Repair transition rate from micro-crack Repair transition rate from radial macro-crack Repair transition rate from circumferential macro-crack 2 x10-2 /yr 20 mm 1 x10-3 /yr 2 x10-2 /yr Repair transition rate from leak 8 x10-1 /yr Leak to rupture transition rate 2x10-2 /yr Macro-crack to rupture transition rate 1x10-5 /yr

Table 4

 4 Important Parameters using a PCE with d=3

1 1 1 1 1 (8)where h( 1 ,…,  n ) is the joint probability density function of .


Zifeng Wang

Gui-Song Xia

Caiming Xiong

Liangpei Zhang

SPECTRAL ACTIVE CLUSTERING OF REMOTE SENSING IMAGES

Keywords: Information mining, remote sensing image clustering, active clustering

Mining useful information from remote sensing images is a longstanding and challenging problem in earth observation, among which images clustering is used to discover meaningful scene information, by grouping similar image pixels into clusters. The main difficulty of image clustering, however, lies in the fact that imperfect similarity measure between images usually leads to bad clustering results. Supervised classification with labeled training samples can partially solve this problem, but the collection of such labeled data is usually time-consuming and sometimes impossible in many real problems. This paper presents an active remote sensing image clustering algorithm by integrating simple human queries into the clustering process. More precisely, we propose a spectral active clustering method that can actively query the oracle (such as human) to improve the image clustering performance. We first construct a k-nearest neighbor (k-NN) graph of the remote sensing images. We then iteratively select the most informative pairwise constraints and purify the k-NN graph, by removing the edges between images from different classes. The final clustering on the purified k-NN graph leads to more accurate result. The proposed method has been evaluated on three high-resolution remote sensing image datasets. It achieves the state-of-the-art performance and demonstrates high potentials in practical remote sensing applications.

INTRODUCTION

Nowadays, remote sensing images capture in details wide surfaces and result in extremely large volumes of data with high spectral and spatial resolution. However, the large amount of remote sensing images are at present weakly exploited due to their large sizes and time-consuming visual analysis [1]. This demands to develop more efficient methods for mining information from these large scale remote sensing images.

In this paper, we are interested in remote sensing image classification. As we know, most of the classification methods require labeled training data from experts, while the annotation of which is usually expensive and sometimes impossible to acquire in many real problems. Alternatively, clustering (or unsupervised classification) approaches are more suitable when no labeled data is available [START_REF] Ng | Efficient and effective clustering methods for spatial data mining[END_REF]. Survey of different clustering methods has been presented in [START_REF] Jain | Algorithms for Clustering Data[END_REF]. One main difficulty of these methods, however, lies in the fact that their performances heavily depend on the similarity matrices between images, which are usually far from perfect on real problem. One possible solution is to integrate acceptable and weak priori knowledge, that are easy to achieve without expertise, into the clustering process, and to build semi-supervised [START_REF] Li | Constrained clustering by spectral kernel learning[END_REF] or active clustering algorithms [START_REF] Eriksson | Active clustering: Robust and efficient hierarchical clustering using adaptively selected similarities[END_REF]6].

This paper introduces a new active clustering algorithm for remote sensing images with weak human queries. More precisely, we propose a spectral active clustering method that can actively query the oracle (such as human) to purify the similarity matrix and finally improve the images clustering performance. In particular, given remote sensing images, we first construct a k-nearest neighbor (k-NN) graph, on which the clustering will perform. To get near-perfect similarity matrices, we propose to interactively purify the k-NN graph, by actively selecting the most informative pairwise constraints. (In our work, the pairwise image constraints are "similar" or "dissimilar" given by an oracle, which is much easier than that of collecting labeled data.) The purified similarity matrices lead to a updated k-NN graph and are finally used to achieve the clustering of images. Unlike previous active clustering method [START_REF] Eriksson | Active clustering: Robust and efficient hierarchical clustering using adaptively selected similarities[END_REF][START_REF] Wang | Active spectral clustering[END_REF], we use vertex uncertainty instead of pairwise uncertainty for selecting informative vertexes on the k-NN graph. To our knowledge, our work is the first time to develop an active clustering method for remote sensing images, which has high potentials in real problem of mining information from earth observation images. The proposed method demonstrates the state-of-the-art classification and annotation results on three high-resolution satellite image datasets.

The rest of this paper is organized as follows. Section 2 recalls the basics of spectral image clustering. Section 3 introduces the proposed spectral active clustering framework for remote sensing images. Section 4 presents experimental results and Section 5 gives some conclusion remarks.

SPECTRAL IMAGE CLUSTERING

Spectral clustering, which makes use of the spectrum of the similarity matrix of the data to perform clustering, is one of the most popular clustering algorithms, due to its efficiency and the simplicity of implementation [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF].

Constructing a k-NN graph of images : Given a set of images (or regions of an image) I = {I 1 , I 2 , • • • , I n }, described by some visual features f i , denote the similarity matrix as S, with each elements s ij representing a measure of the similarity between f i , f j of I i , I j . The similarity graph is thus defined by G = (V, E), where each vertex v i ∈ V represents an image I i and two vertices I i , I j are connected if s ij is larger than a certain threshold, and the edge e ij ∈ E is weighted by s ij . In our case, we refer the graph G as a k-NN graph, which implies that for a vertex v i only the vertices with the first k largest similarity measure are connected to it.

Spectral clustering algorithm : Spectral image clustering [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] then works with the eigenvectors of the normalized Laplacian matrix L of G, using k-means algorithm. The normalized Laplacian matrix L is computed as

L = I -D -1/2 SD -1/2 ,
where I is the identity matrix and D is the diagonal matrix as D ii = n j=1 s ij . In order to deal with large scale remote sensing images, some large scale spectral clustering algorithm can be used to compute L more efficiently.

SPECTRAL ACTIVE CLUSTERING OF REMOTE SENSING IMAGES

This section presents the proposed spectral active clustering approach for remote sensing images. The flowchart of the algorithm is depicted in Fig. 1. Given a set of image (or image regions) as inputs, we first construct the k-NN graph followed by a spectral clustering algorithm, as described in Section 2. We iteratively obtain new constraints by actively selecting the most informative image and querying the oracle. Active learning helps us to find out the most informative image, which is also the most uncertain one, according to the current clustering result and the k-NN graph. With new constraints, the k-NN graph is purified and spectral clustering performs again on the updated k-NN graph. The algorithm iterates this process until the oracle is satisfied with the result or the k-NN graph is fully purified.

In what follows, we precise the criterions for active constraints selection and the k-NN graph purification.

Active Constraints Selection Ideally, in a perfect k-NN graph, each image should belong to the same class as all its neighbors do. Consequently, vertices connected in the k-NN graph will be assigned to the same cluster. In order to look for new constraints to purify the k-NN graph, we try to find out the image whose neighbors contain the "worst edge" (edge linked images belonging to different classes). For the class label c of vertices, we use cluster labels instead of the real labels which are still unavailable. In order to measure the disordered level of an image I i on the graph G, we use an entropy criterion defined as

H(I i ) = - P Ii ( ) log P Ii ( ), (1)
where

P Ii ( ) = #{cv= ,v∈N I i } k
is the probability of the neighbors that are assigned to the cluster , with N Ii as the neighborhoods of I i on G. We thus choose the most disordered image I i * , i.e. I i * = arg max Ij ∈G H(I j ), as the vertex on the graph for querying. Consequently, the active edges are thus {(I i * , I j ); e i * j = 1, I j ∈ G}, which link the image I i * and its neighbors. We query the oracle along these edges to obtain new constraints.

k-NN Graph Purification As mentioned before, there are two kind of pairwise constraints: similar and dissimilar, called "must-link" and "cannot-link" on the graph which respectively link images to the same and different classes. A simple constraint augmentation process is described in Fig. 2. It enables us to obtain more constraints from known ones. By purifying the k-NN graph, all the "cannot-link" edges in the graph will be removed, while the similarity values of "must-link" edges are set to 1.

EXPERIMENTAL RESULTS

We evaluates the proposed spectral active clustering method on remote sensing image classification and annotation.

The three remote sensing image sets we used are: (1) A GeoEye-1 image in Beijing as in Figure 4: the image is of size 4000 × 4000 and is partitioned into 1600 Fig. 3. The evaluation of the similarity matrix with the iterative purification of the k-NN graph. From left to right: the similarity matrices at 1, 75, 150 and 225 iterations. From top to bottom: the similarity matrices of the GeoEye-1 image in Beijing, the WhuSPL datset [START_REF] Xia | Structural high-resolution satellite image indexing[END_REF], and the UCM datset [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF]. Observe that with the iterations the similarity matrices become more and more discriminative.

tiles with size of 100 × 100 pixels. We try to annotate the tiled images as in [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF], who also provides its human-labeled groundtruth. (2) WhuSPL dataset [START_REF] Xia | Structural high-resolution satellite image indexing[END_REF], containing 20-class high-resolution satellite images with about 50 samples per each categories. (3) UCM dataset [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF], which contains 21class high-resolution satellite images with 100 samples per each categories.

In our experiments, to characterize each image or image tile, we calculate the bag of dense SIFT descriptors [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and bag of color descriptors in each image. Thus, the image I i is described by a histogram-like vector h i of length M . The pairwise similarities matrix S is calculated by histogram intersection kernel as

s ij = M m min(h i [m], h j [m]).
We compare our algorithm with three currently proposed methods: spectral clustering (SC) algorithm [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], semi-supervised spectral clustering (S 3 C) algorithm [START_REF] Li | Constrained clustering by spectral kernel learning[END_REF] and the M 3 DA-RF algorithm [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF], a recently proposed fullsupervised method on our task. Remark that, instead of collecting pairwise constraints actively, S 3 C [4] performs spectral clustering by using pairwise constraints provided before clustering. The M 3 DA-RF [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF] is full-supervised and requires labeled data from experts as training sample. In current experiments, our method exits the active clustering loop when the oracle is satisfied with the results. For a fair comparison, we use the same number of pairwise constraints for the S 3 C [4] algorithm. In order to evaluate the performance of clustering algorithms (as that for supervised algorithms, e.g. M 3 DA-RF), the clustering accuracy is computed by assigning the label of each cluster as its closest class in the groundtruth. Fig. 3 illustrates the evaluation of the similarity matrices with the iterative purification of the k-NN graph on each data set. Note that, on all the three dataset, with more active iterations (i.e. more queries from oracle) are used the similarity matrices become more and more discriminative. This verifies the efficiency and the necessariness of the active selection procedure of our proposed method.

We compare the clustering performances of SC, S 3 C and our method on WhuSPL dataset [START_REF] Xia | Structural high-resolution satellite image indexing[END_REF] and UCM dataset [START_REF] Yang | Bag-of-visual-words and spatial extensions for land-use classification[END_REF]. The accuracies are respectively 39.2%, 58.5% and 98.8% on WhuSPL dataset and 35.7%, 47.9% and 96.3% on UCM dataset of the three algorithms. Fig. 4 displays the comparison of our method with the three approaches mentioned above, SC [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], S 3 C [4] and M 3 DA-RF [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF], on the annotation of a big high-resolution satellite image, a GeoEye-1 image in Beijing. The M 3 DA-RF [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF] takes half of the labeled data as training samples and uses multi-level max-margin discriminative random field [START_REF] Hu | Tile-level annotation of satellite images using multi-level max-margin discriminative random field[END_REF] for annotation. Notice that M 3 DA-RF uses actually spatial constraints via conditional random field to improve its performance. From Fig. 4, one can see that S 3 C gives much better result (61.8%) than that of SC (35.8%), which implies that pairwise constraints is very helpful for this task. The M 3 DA-RF produces better result (91.6%) than S 3 C does. This is mainly because it uses fully labeled training samples, the prior information of which is much stronger that that of weak pairwise constraints. Our method outperforms all the three methods and shows near-perfect annotation result (99.2%). This finally verifies our intuition on the proposed method, that the weak prior knowledge from the actively selected pairwise constraints can help the clustering algorithm a lot.

CONCLUSION

This paper addresses the problem of remote sensing image clustering with weak priors. Our method actively queries the oracle to get weak pairwise constraints. Instead of pair uncertainty, the proposed method uses node uncertainty as active select criterion, which can more accurately select useful queries. With these acceptable a few pairwise constraints, the clustering result shows notable improvements. It is very powerful to classify remote sensing images, when there is no available labeled training data or it is hard to acquire such training data. From experiment results, we can see that the proposed method outperforms three state-of-the-art methods and shows high potentials in mining meaningful information from remote sensing images.

Fig. 1 .

 1 Fig. 1. Flowchart for spectral active clustering. Refer to the text for more details.

Fig. 2 .

 2 Fig. 2. Our constraint augmentation process.

Fig. 4 .

 4 Fig. 4. Comparison of different methods on the annotation of a GeoEye-1 image in Beijing. Refer to the text for more details.