<figure><img src="image_1.png" /></figure>

HAL
open science

Modeling time and spatial variability of degradation
through gamma processes for structural reliability
assessment

Mestapha Oumouni, Franck Schoefs, Bruno Castanier

► To cite this version:

Mestapha Oumouni, Franck Schoefs, Bruno Castanier. Modeling time and spatial variability of degradation through gamma processes for structural reliability assessment. Structural Safety, 2019, 76, pp.162-173. 10.1016/j.strusafe.2018.09.003. hal-02440519

HAL Id: hal-02440519

https://hal.science/hal-02440519

Submitted on 17 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Modeling time and spatial variability of degradation
through gamma processes for structural reliability
assessment

M. Oumouni<sup>a,*</sup>, F. Schoefs<sup>a</sup>, B. Castanier<sup>b</sup>

<sup>a</sup>University of Bretagne Loire, Université de Nantes, Research Institute in Civil and
Mechanical Engineering (GeM), UMR CNRS 6183, 2 rue de l'Houssinire, BP 92208, 44322
Nantes, France

<sup>b</sup>University of Bretagne Loire, Université d'Angers, Research Institute in Systems
Engineering of Angers (LARIS), Angers, France

Abstract

The objective of this work is to propose a spatio-temporal random field for as-
sessing the effect of spatial variability on the degradation process in structural
reliability assessment. Our model extends the classical Gamma process which
is usually developed for degradation temporal variability concerns to integrate
spatial variability and heterogeneity issues. A log-normal distributed spatial
random scale is then introduced in the Gamma process. Mathematical models
for a structure degradation in space and time and estimation procedures are
developed in this paper. Approximations and simulations are given to evaluate
the failure time distribution and to characterize the residual lifetime after in-
spection. Simulation results are performed using pseudo-random data based on
Monte-Carlo simulations to fit the model and the inference of their parameters.
The two proposed estimation method - Method of Moments and Pseudo Maxi-
mum Likelihood - are numerically compared according to statistical measures.
The accuracy of the methods are also discussed by numerical examples given
the approximations of quantities of interests.

Keywords: Degradation model, Gamma process, Random field, Reliability,

*Corresponding author
Email addresses: mestapha.oumouni@univ-nantes.fr (M. Oumouni),
franck.schoefs@univ-nantes.fr (F. Schoefs), bruno.castanier@univ-angers.fr (B.
Castanier)
# 1. Introduction

In recent years, the scientific material and structural engineering community pays a lot of attention to the elaboration of mathematical degradation models for integrating spatial variability and uncertainty. It is obvious that the deterioration of structures exposed to environmental conditions is spatially and temporally varying. The variation along the space is caused by the inherent uncertainty through the material at several positions on the structure and the physical parameters involved in the deterioration mechanism. This material non-homogeneity problem is well known for steel and concrete structures and has been, e.g., studied in [1, 2, 3, 4, 5] through concrete diffusion property, concrete cover and chloride external concentration.

A common way to deal with such uncertainties is to use a probabilistic framework by modeling the input data with random fields [2, 10, 19, 21, 34]. For instance, a randomized salty environment (de-icing of the sea natural salts) can be integrated in the transport equation to compute the concentration of the chloride in concrete, in order to estimate the time until failure under uncertainties. In [3, 5] the authors consider a two-dimensional Gaussian random field with a Gaussian correlation to compute the likelihood of corrosion-induced cracking in reinforcing steel bar. They provide an estimation of the time to first crack and time to limit crack widths. In [1], authors propose a probabilistic model for steel corrosion in reinforced concrete structures considering crack effect on the corrosion mechanism, in which an empirical model for the crack propagation stage is developed by the standard gamma process and combines corrosion crack width with steel-bar cross sectional loss. The main disadvantage of these approaches is the problem of larger dimension so called "curse of dimensionality" in which the resolution of a large number of deterministic problems is involved.
Meta-models are commonly used to tackle this problem of curse of dimen-
sionality for degradation prediction in structural engineering. The time-dependent
deterioration processes are modeled by a stochastic model where only the varia-
tion in time is studied. In particular a standard Gamma process is an appropri-
ate mathematical model for predicting deterioration encountered in civil engi-
neering [6, 35]; such as corrosion and crack of reinforced concrete. The authors
in [17, 18] developed a probabilistic framework in which interaction between
shocks and gradual process are combined in view to describe the deterioration
process. They proposed a semi-analytic computation model to estimate with
less cost the time to failure. The work in [9] describes the degradation by a
Gamma process and includes other dependent-parameters as covariate in the
shape function. In [26], a Gamma process with a random scale following a
gamma law is considered to model heterogeneity in the degradation data and
obtained analytic results for reliability assessment. Markov chain is another
widely used approach to model cumulative damage. It is seen as a discrete
Markov process where the deterioration is assumed to be a single step function.
The estimation of the one-step transition matrix requires a large number of
transitions to estimate all its elements [16, 28]. The drawback of this approach
is that time variability is difficult to capture.

In [15], the authors propose the construction of the state-dependent degra-
dation model based on the Gamma process, where the cracking of a submerged
concrete structure subjected to corrosion is described by the proposed bi-variate
model with a suitable parameters. In the same time, the authors in [23] intro-
duce a state-dependent Gamma process for the degradation, the dependency
is modelled in the scale function instead of the shape function as proposed in
[15]. The authors in [24] introduce a bi-variate spatio-temporal field to model
the action induced by the walking of a small group of persons. Based on the
spectral and coherence functions of the forces, they proposed an evaluation of
vertical and transversal accelerations at nodes of a finite element.

All models based on temporal variability assume a uniform degradation and
do not integrate the spatial variability through degradation process. Neverthe-
less, recent studies have shown that this spatial correlation has an important and direct impact on the level of structural reliability estimates [32, 3, 4]. Therefore, incorporating these uncertainties in the degradation processes through mathematical modelling improves their prediction and versatility in terms of maintenance and decision. On the other hand, to construct an accurate model of the degradation, a large amount of data using destructive or nondestructive testing is required from a large amount of structures [31, 32]. Therefore, one way of obtaining accurate and reliable information is to embed the spatial variability in the models. This allows to increase the relevance of the Meta-Model approach where the uncertainties are reduced, improving the accuracy of the inference and extending the use of the non-destructive testing.

Therefore, the major contribution of this paper is a new spatio-temporal random model based on Gamma process for predicting a single degradation measure which takes into account both temporal and spatial variability. Under the stationary assumption, the spatial monitoring data of the structure contributes in the parameters estimate to increase the accuracy of the meta-model approach. Our model requires only few parameters and is thus very suitable for inference when only few components are inspected: that is the case for on-site inspection of civil engineering structures or marine structures where the cost of inspection is high.

The spatio-temporal degradation model is assumed to be an observable process in space and time with limited observations (Non Destructive Testing, distributed sensors). However, such kind of database that considers both hazard time and space of the degradation are not available in the literature. In order to validate the proposed inference framework, we construct a synthetic discrete degradation model through Monte Carlo simulations. Therefore, numerical experiments will be conducted and compared for identifying preliminary properties and advantages of our model in terms of statistical inference and computation of quantities of interest for reliability and maintenance.

The article is organized as follows: Section 2 introduces the classical Gamma process for temporal variability, the Gaussian random field with its simulation
method and the construction of degradation model is detailed. Section 3 develops quantities of interest which are useful in the reliability analysis, namely the distribution of the failure time and the distribution of the remaining lifetime of the unit. Section 4 compares approaches for identifying properties of the model in terms of statistical inference. Finally, Section 5 presents a numerical example illustrating the proposed methodology for model validation.

## 2. Degradation Processes and Random field variability

### 2.1. Standard time variant Gamma Process modeling

The standard Gamma process (GP) is an appropriate mathematical model for modeling the degradation evolution in structural engineering, such as corrosion and cracks of materials, which are the common causes of structural failure. The deterioration is supposed to take place gradually over time in a sequence of tiny increments. Consider $\alpha(\cdot)$ to be a non-decreasing, right-continuous, real-valued function for $t \ge 0$ and vanishing at $t = 0$.

#### Definition 2.1.

A stochastic process $(X_t)_{t \ge 0}$ is said to be a GP with shape function $\alpha(t)$ and identical scale parameter $\beta > 0$ if the process satisfies the following properties:

* $X_0 = 0$ with probability one,

* The increment $X_{t+s} - X_t$ has a Gamma distribution $Ga(\alpha(t+s) - \alpha(t), \beta)$,

* $X_t$ has independent positive increments,

where the Gamma distribution $Ga(\alpha, \beta)$ is defined by the density function:

$$f_{Ga}(x) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}, \text{ for each } x > 0, \quad (1)$$

and $\Gamma$ is the classical Gamma function. The process $X_t$ is said to be stationary if $\alpha(\cdot)$ is a linear function and $X_t$ is non-stationary if $\alpha(\cdot)$ is a non-linear function. The mean and variance of $X_t$ are $E[X_t] = \frac{\alpha(t)}{\beta}$, $var[X_t] = \frac{\alpha(t)}{\beta^2}$, respectively. The process $X_t$ satisfies the scaling property,

$$\gamma X_t = Ga(\alpha, \beta/\gamma), \text{ for each } \gamma > 0, \quad (2)$$
and its logarithm $\log(X_t)$ has the following two first moments:

$$ \mathbb{E}[\log(X_t)] = \psi(\alpha(t)) - \log(\beta), \quad (3) $$

$$ \mathrm{var}[\log(X_t)] = \psi_1(\alpha(t)), \quad (4) $$

where the function $\psi$ is digamma function which is defined as the logarithmic derivative of $\Gamma$, and $\psi_1$ is the trigamma function defined as the derivative of $\psi$.

## 2.2. Gaussian Spatial random field modeling

Spatial variability of material properties is classically modeled by a second-order stationary random field (RF) given by a non-linear transformation $\mathcal{T}(\cdot)$ of a Gaussian random field (GRF) $\mathcal{T}(Y)$ [8, 31]. For example, in the transport equation, the diffusivity coefficient is modeled by a log-normal RF [1, 2, 19], where its distribution is obtained as a limit of physical positive quantities.

We consider $Y(z, \omega)$ to be a spatial GRF in the set $D \times \Omega$, where $D$ is a set in $\mathbb{R}^d$ with $d = 1, 2, 3$ and $\Omega$ is an abstract set of events. The GRF $Y$ is assumed to be homogeneous and then totally defined by its mean $\mu \in \mathbb{R}$ and its covariance function $cov(r)$ [7, 36]. The $cov(r)$ function models the correlation between two spatial random variables on any two points separated by the distance $r$.

In order to simulate $Y$ on $\{z_0, z_1, \dots, z_N\} \subset D$ a set of equidistant points, we choose the circulant embedding matrix approach [14] (also named by DSM discrete spectral method in [11]). This method is a very versatile approach for generating GRF, the discretized RF has the same spatial correlation on the grid points. In [30] the authors develop a continuous spectral method to simulate $Y$ where the spectral density is discretized on a uniform grid, then a Discrete Fourier Transform (DFT) is used to generate an approximation of $Y$. The proposed simulated GRF is asymptotically Gaussian where its correlation structure is an approximation of the target correlation and its accuracy is strongly related to the regularity of $Y$ (see [22]). Another widely used approach to generate $Y$ is the Karhunen-Loève expansion. However, it gives only an approximation of $Y$ by a truncation of an infinite series, where its accuracy depends strongly on the smoothness of $Y$ and its correlation length, [20, 21, 31].
The DSM is based on the matrix factorization approach of the convenient positive definite circulant correlation matrix R using DFT. It has a Toeplitz structure and defined by its first row $r = (r_0, r_1, \dots, r_{N-1}, r_N, r_{N-1}, \dots, r_1)$ where each $r_l = cov(z_l)$, for $l = 0, \dots, N$. For instance, in one-dimensional uniform grid $z_0 = 0, \dots, z_N$, a realization of the random field Y is simulated by the following components:

$$Y_N(z_l) = \frac{1}{\sqrt{2N}} \sum_{k=0}^{2N-1} \sqrt{\lambda_k} \left( \zeta_k \cos 2\pi \frac{kl}{2N} + \zeta_k \sin 2\pi \frac{kl}{2N} \right), \quad (5)$$

where $\lambda_k$ are positive and given by the DFT of the circulant vector $r$. The set $\{\zeta_k\}_{k=0}^{2N-1}$ are independent random variables with standard normal distribution $N(0, 1)$. The random vector $(Y_N(z_0), \dots, Y_N(z_N))$ given in (5) has the same correlation structure of Y on the grid $(z_l)_{l=0}^N$. The Fast Fourier Transform is used to construct exact simulation of Y with $O(N \log(N))$ operations.

140

## 2.3. Matérn-Whittle model

The Matérn covariance functions are commonly used for GRF. There are defined by:

$$cov(r) = \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu r}}{l_c} \right)^{\nu} \mathcal{K}_{\nu} \left( \frac{\sqrt{2\nu r}}{l_c} \right) \quad (6)$$

where $\nu$ and $l_c$ are non-negative real numbers, $l_c$ is the correlation length, $r$ is the Euclidean distance between two points, $\mathcal{K}_{\nu}$ denotes the modified Bessel function of the second kind. The regularity parameter $\nu$ tackles different models: a large value of $\nu$ implies that Y is ([$\nu$] - 1)-times differentiable and small value of $\nu$ implies that Y is rough. When $\nu = \frac{1}{2}$, cov coincides with the exponential covariance,

$$c(r) := e^{-r/l_c}, \quad (7)$$

Y are only continuous. When $\nu \to \infty$, it tends to the gaussian model,

$$c(r) := e^{-r^2/(2l_c^2)}, \quad (8)$$

which is an analytic function and so for the samples paths of Y. Both models are the most used in structural engineering applications [1, 3, 19].
## 2.4. Spatial random field scale for Gamma process

We look for modeling a degradation process by a convenient spatio-temporal RF to consider its aleatory evolution in time and space. A separable model is one simple spatio-temporal model obtained through the tensorial product between a merely stochastic process $(X_t)_{t \ge 0}$ and a spatial RF $Z(z)$, where $z$ is the spatial variable. This class of separable RF is extensively used even in situations in which they are not always physically justifiable since separability gives important computational benefits.

We are interested herein to model the deterioration in structural engineering under the presence of spatial variability in the model. The evolution in time models the intrinsic aleatory while the spatial RF models the variability and uncertainty through the structure. Thus, for simplicity of the model, we assume that randomness in time and in space are independent. In the classical approach for modeling the deterioration over time, the GP is ideally suited to model gradual deterioration which monotonically accumulates over time. This process can be extended spatially to obtain a spatio-temporal random field by considering its parameters to be spatial random fields, i.e the shape function $\alpha(\cdot)$ or the scale parameter $\beta$. The scaling property of the GP in (2) motivates to model the scale parameter with a spatial RF and then to obtain a separable spatio-random RF,

$$G_t(z) := X_t \beta(z)^{-1}, \qquad (9)$$

where $X_t$ is the GP with a shape function $\alpha(t)$ and unit scale parameter $Ga(\alpha(t), 1)$. The positive spatial RF $\beta(\cdot)$ is assumed to be independent of $X_t$. The scaling property satisfied by $X_t$ in (2) motivates us to take the spatio-random field $G_t(z)$ as GP with spatial random scale $\beta(\cdot)$.

In practice, it is difficult to find a convenient positive distribution for the spatial RF. However, some several reasons (detailed above) suggest to choose a Log-normal distribution for the random scale coefficient,

$$\beta(z) = e^{\sigma Y(z) + \mu}, \qquad (10)$$
where $Y$ is a GRF with zero mean, unit variance and defined by the correlation function in (6). The constant $e^\mu$ is seen as the deterministic contribution of the scale parameter $\beta(\cdot)$.

The Log-normal distribution appears naturally as a limit law of physical processes. The Central Limit Theorem applied to the product of positive independent random variables (number of measures greater than 30) ensures that the log-normal distribution occurs. Further, a Log-normal RF is completely defined by its covariance. The maximum likelihood method is a commonly used inference procedure for such field. Another remarkable properties is that $e^{\sigma Y}$ and $e^{-\sigma Y}$ have the same finite-dimensional law, which implies that $Ga(\alpha, e^{\sigma Y+\mu})$ and $Ga(\alpha, e^\mu)e^{\sigma Y}$ have the same finite-dimensional law. However, since in this work, we use only the method of moments to perform the inference procedure; we can consider any choice of positive RF.

# 3. Failure time and reliability evaluation

## 3.1. Marginal density approximation

The marginal distribution noted by $f_t(\cdot)$ of the model is approximated in order to compute some quantities of interest which are used in structural reliability and maintenance. The GRF $Y(z)$ considered in this work for the spatial variability is homogeneous, thus $f_t(\cdot)$ does not depend on the position $z$ but only on the time $t$. Let $\xi(y)$ be the density of the standard Gaussian random variable $N(0, 1)$ and setting $\eta = e^\mu$, the marginal density of $G_t|y$ is a gamma distribution with the shape parameter $\alpha(t)$ and a scale parameter $\eta e^{\sigma y}$ for each $t > 0$. Then, the pdf $f_t$ is given for all $v > 0$ by the following form:

$$f_t(v) = \frac{\eta^{\alpha(t)} v^{(\alpha(t)-1)}}{\Gamma(\alpha(t))} \int_{\mathbb{R}} \exp\left(-v\eta e^{\sigma y} + \alpha(t)\sigma y\right) \xi(y) dy. \quad (11)$$

The integral in (11) has a transcendental form, thus we use Gauss-Hermite quadrature formula to approximate this marginal density $f_t$. We consider $m$ roots $\{y_j\}_{j=1}^m$ of the Hermite polynomial and their associated weights $\{w_j\}_{j=1}^m$.
Thus, an approximation of $f_t$ writes:

$$f_t^m(v) := \frac{\eta^{\alpha(t)} v^{(\alpha(t)-1)}}{\Gamma(\alpha(t))} \sum_{j=1}^{m} \exp\left(-v\eta \exp(\sigma y_j) + \alpha(t)\sigma y_j\right). \quad (12)$$

The convergence of the sequence $f_t^m(v)$ for a fixed $v > 0$ is relatively fast. However, the norm of any $m$-derivative of the integrand depends on the value of the parameters $\alpha(t)$, $\sigma$ and $\eta$. A large value of these parameters requires a large order $m$ of the approximation in (12), in particular for large time $t$. Then, the order $m$ is built by the following stop criterion,

$$|f_t^m - f_t^{m-1}| \le \epsilon, \quad (13)$$

where $\epsilon > 0$ is a convenient threshold value.

## 3.2. Failure Time Distribution

The failure time $T_F$ for a component is defined as the time at which the degradation path $G_t$ first crosses a critical level $g_F$ for any spatial location,

$$T_F = \inf\{t; G_t(\cdot) \ge g_F\}, \quad (14)$$

In what follows, the critical level $g_F$ is assumed to be deterministic. For some simple path models, the distribution $F_T(t) := P(T_F < t)$ of $T_F$ can be expressed in a closed form. However, this is not always possible and it can be numerically computed with Monte Carlo simulations by generating multiple paths of $G_t(z)$.

For the spatio-temporal random model considered here, sample paths are monotonic. Thus the failure time cumulative distribution $F_T(T)$ satisfies:

$$\begin{aligned} F_T(t) &= 1 - P(T_F > t) = 1 - P(G_t \le g_F) \\ &= 1 - \int_0^{g_F} f_t(z) dz, \end{aligned} \quad (15)$$

where $f_t$ is the marginal probability density function of the degradation process $G_t$ given in (11). Therefore, we can approximate the distribution $F_T$ by using an approximation $f_t^m$ of the density function $f_t$ with a suitable order $m$:

$$F_T^m(t) = 1 - \int_0^{g_F} f_t^m(v) dv, \quad (16)$$
The integral (16) can be computed accurately by Legendre-Gauss quadrature formula. The derivative of (15) and (16) with respect to the time variable t provides the probability density function of $T_F$ and its approximation, respectively.

The approximation of $f_t^m$ can be inaccurate when the time $t$ or the variance become large. And the computation of the cdf (16) will require a huge cost. In this case, an estimation of $F_T$ is provided using Monte-Carlo (MC) simulations where a sufficiently large number of sample paths of $G_t(z)$ are conducted.

Let $N_t$ be the desired times, $N_z$ the desired locations and $M$ the realizations of the spatio-temporal paths $\{G_{t_i}^m(z_j)\}$ for $i = 1..., N_t$, $j = 1..., N_z$ and $m = 1, ..., M$. The estimate $\tilde{F}_T$ of $F_T$ is given by:

$$ \tilde{F}_T(t_i) := \frac{\sum_{j,m}^{N_z,M} \mathbb{I}\{G_{t_i}^m(x_j) \ge g_F\}}{MN_z}, \quad (17) $$

where $\mathbb{I}_A$ represents the characteristic function of the set $A$, i.e $\mathbb{I}_A(z) = 1$ if $z \in A$ and zero otherwise. Note that since $G_t(z)$ is homogeneous with respect to the spatial variable, the estimate $\tilde{F}_T$ can be computed also using realizations of $G_t(z_p)$ fixed at any position $z_p$,

$$ \tilde{F}_T(t_i) \approx \frac{\sum_m^M \mathbb{I}\{G_{t_i}^m(z_p) \ge g_F\}}{M}. \quad (18) $$

However, estimation (18) requires more MC simulations of $G_t(z)$ than (17) since the spatial average contributes in the convergence of $\tilde{F}_T$ to $F_T$ (ergodic property).

### 3.3. Remaining Lifetime after inspection

In reliability analysis and survival studies, residual lifetime after inspection is a key indicator. In the maintenance decision analysis, the current measured degradation is used to predict the remaining lifetime (RL) of the structure [33]. If $t$ is the current time of inspection, the residual lifetime is defined by the random variable:

$$ RL_t := \inf \{ \tau > 0, G_{t+\tau} \ge g_F | G_t = g_t \}, \quad (19) $$
where $g_F$ is the critical level and $g_t$ is the measured degradation at a given time $t$. We have implicitly $g_t < g_F$. The Markov property of the model ensures that, from any current state, future states can be predicted. However, if we suppose that a component has survived to a given time $t$ and we have no information or measure about the current degradation path $G_t$, then a conditional reliability function gives an evaluation of the remaining lifetime:

195

$$
\begin{align}
R(\tau|T_F > t) &:= P(T_F \ge \tau + t|T_F > t) \nonumber \\
&= \frac{P(G_{\tau+t}(\cdot) \le g_F)}{P(G_t(\cdot) \le g_F)} \nonumber \\
&= \frac{\int_0^{g_F} f_{\tau+t}(y)dy}{\int_0^{g_F} f_t(y)dy}. \tag{20}
\end{align}
$$

When the current degradation measure path of $G_t$ is available. The probability that the unit survival after time $t + \tau$ given its current state $G_t = g_t$ at time $t$ is:

$$
\begin{align}
P(RL_t > \tau) &= P(G_{t+\tau} < g_F | G_t = g_t) \nonumber \\
&= P(\delta_{\tau} G_t < g_F - g_t | G_t = g_t), \tag{21}
\end{align}
$$

200

where $\delta_{\tau}G_t := G_{t+\tau} - G_t$ is the associated degradation increment. Let denote $f_{\delta_{\tau}G_t|G_t}$ the conditional marginal density of the process $\delta_{\tau}G_t$ given the event $\{G_t = g_t\}$. The cdf $F_{RL} := P(RL_t \le \tau)$ of the residual lifetime becomes:

$$
\begin{align}
P(RL_t \le \tau) &= 1 - P(\delta G_t < g_F - g_t | G_t = g_t) \nonumber \\
&= 1 - \int_0^{g_F - g_t} f_{\delta_{\tau} G_t | G_t}(y) dy. \tag{22}
\end{align}
$$

Increments of the model are independent given a realization of the random field $e^{\sigma Y}$. The marginal density of the bivariate variable $(\delta G_t, G_t)$ given $y$ denoted $\pi_y(u, v)$ for any $u, v > 0$ can thus be defined as the product of the density function of $\delta G_t$ given $y$ and that of $G_t$ given $y$:

$$
\pi_y(u, v) = \Psi(u, v) \exp\left(-(u+v)\eta e^{\sigma y} + \alpha(t+\tau)\sigma u + \alpha(t)\sigma v\right), \quad (23)
$$
where $\Psi(u, v) := \frac{u^{\delta_{\tau}\alpha-1}v^{\alpha(t)-1}}{\Gamma(\delta_{\tau}\alpha)\Gamma(\alpha(t))}\eta^{\alpha(t+\tau)}$ the increment $\delta_{\tau}\alpha := \alpha(t+\tau) - \alpha(t)$ and $\eta = e^{\mu}$.

The joint probability density of the bivariate variable $(\delta_{\tau}G_t, G_t)$ is,

$$ \pi(u, v) = \int_{\mathbb{R}} \pi_y(u, v)\xi(y)dy. \qquad (24) $$

Therefore, the conditional density of $\delta_{\tau}G_t$ given the event $\{G_t = g_t\}$ is,

$$ f_{\delta_{\tau}G_t|G_t}(u) = \frac{\pi(u, g_t)}{f_t(g_t)}, \qquad (25) $$

where $f_t$ is the marginal density of $G_t$ given in (11). When $g_t > 0$ and $t > 0$, this conditional probability density is given by the formula:

$$ f_{\delta_{\tau}G_t|G_t}(u) = \frac{u^{\delta_{\tau}\alpha-1}g_t^{\alpha(t)-1}}{B(\delta_{\tau}\alpha, \alpha(t))(u+g_t)^{(\alpha(\tau+t)-1)}} \frac{f_{\tau+t}(u+g_t)}{f_t(g_t)} \qquad (26) $$

where $B(x, y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ is the beta function. Given the approximations of $f_{t+\tau}^m(u+g_t)$ and $f_t^m(g_t)$, $f_{\delta_{\tau}G_t|G_t}^m(u)$ is a good approximation of the conditional probability density in (26). It follows that $F_{RL}(\tau)$ is given by,

$$ F_{RL}(\tau) \approx 1 - \int_{0}^{g_F - g_t} f_{\delta_{\tau}G_t|G_t}^m(u)du. \qquad (27) $$

Note that when $t = 0$ and $g_t = 0$, equation (26) is not definite and the cumulative probability function $F_{RL}$ is the function $F_T$ given in (15). The derivative of (22) with respect to the variable $h$ provides the probability density function of the Residual lifetime RL.

## 4. Parameters inference of the model

The estimation parameters for spatio-time fields is not so thoroughly developed in the literature. The Maximum likelihood Method (MLE) is considered as the best method for this purpose. However, it assumes Gaussian or Log-normal distribution [25] unlike the Method of Least Squares (MLS). The MLE is promoted due to its ease of implementation (see, [13, 29]) where it only requires a preliminary step to obtain a nonparametric estimate of the variogram or covariance by the method of moments.
Our degradation model is a separable random-field given by a tensorial product of Gamma process and a Log-normal RF. The estimation process is in two phases: the space parameters estimation for GRF is performed, and then the temporal parameters of GP.

220

A typical data set consists of inspection points in different increasing times $t_j$ for $j = 0, \dots, N_t$, with the same period $\tau$. For each time $t_j$, the inspection positions are given in uniform positions $z_l$ for $l = 1, \dots, N_z$ with step size $h$.

For sake of simplicity, $\alpha(\cdot)$ is assumed here to be the power law:

$$ \alpha(t) = at^b, \qquad (28) $$

for some unknown $a > 0$ and a known power $b > 0$.

225

## 4.1. Method of Moments (MOM)

### 4.1.1. Spatial parameters estimation

The spatial parameters are the parameters of the second order stationary RF $Y$, i.e. the variance $\sigma^2$, the correlation length $l_c$ and the regularity parameter $\nu$ of the correlation function in (6). The correlation length $l_c$ is the main parameter which quantifies the spatial uncertainty. Recent reliability studies in structural engineering have shown that this length plays an important role in the level of structural reliability [32, 3, 4]. The following nonparametric estimate of the semi-variogram based on the method of moments (MOM) from realizations of $Y$ is proposed:

$$ \hat{\Upsilon}_Y(h_l) = \frac{1}{2N_{h_l}} \sum_{z_i, z_j \in S_{h_l}} (Y(z_i) - Y(z_j))^2, \qquad (29) $$

where $S_{h_l}$ is the set of the points separated with distance $h_l$ and $N_{h_l}$ its cardinal. The variogram is invariant by translation with any random variable. Then, the empirical variogram $\hat{\Upsilon}_Y$ is provided by computing the variogram of the random field $\log(G_t)$. At any fixed time $t$, we get:

230

$$ \begin{aligned} \Upsilon_{\log(G_t)}(h_l) &= \frac{1}{2N_{h_l}} \sum_{z_i, z_j \in S_{h_l}} \left( \log(G_t(z_i)) - \log(G_t(z_j)) \right)^2 \\ &= \frac{1}{2N_{h_l}} \sum_{z_i, z_j \in S_{h_l}} \left( Y(z_j) - Y(z_i) \right)^2 = \hat{\Upsilon}_Y(h_l). \end{aligned} $$
Note that the estimation of $\hat{\Upsilon}_Y(h_l)$ can be improved by the empirical average when $M$ realizations of the model $G_t(z)$ are available. The RF $Y$ is stationary. So the theoretical variogram is given by,

$$ \Upsilon_Y(h_l) = \sigma^2 - \sigma^2 \operatorname{cov}(h_l), \qquad (30) $$

where the correlation function $\operatorname{cov}(h)$ is given by (6). Therefore, the spatial parameters are estimated by minimizing the quadratic error (Least square method) between the theoretical and the experimental variograms. Thus, $\sigma^2$, $lc$ and $\nu$ are deduced by the following minimization problem:

$$ \min_{\sigma,lc,\nu>0} \sum_{l=1}^{N_z} \left( \hat{\Upsilon}_Y(h_l) + \sigma^2 \operatorname{cov}(h_l) - \sigma^2 \right)^2 \qquad (31) $$

The classical least square method is sometimes not efficient, in particular when the length of the space and the number of positions are not large enough, because the set $\{\hat{\Upsilon}_Y(h_l)^2\}_{l=1}^{N_z}$ are not independent and with various variance. The generalized least squares method should be preferred:

$$ \min_{\sigma,lc,\nu>0} \sum_{l=1}^{N_z} \left( \hat{\Upsilon}_Y(h_l) - \Upsilon_Y(h_l) \right) R_{l,k} \left( \hat{\Upsilon}_Y(h_k) - \Upsilon_Y(h_l) \right), \qquad (32) $$

where $\Upsilon_Y(h_l)$ is given in (30) and the matrix $R$ is the correlation matrix of the set $\{\hat{\Upsilon}_Y(h_l)\}_{l=1}^{N_z}$ (see [29], for more details). A simplified approach is the weighted method where $R$ is restricted at the diagonal matrix defined by entries $R_l = \frac{2\Upsilon_Y(h_l)}{N_{h_l}}$ as suggested in [13] assuming a Gaussian law and no correlation among $\{\hat{\Upsilon}_Y(h_l)\}_{l=1}^{N_z}$.

**Comments 4.1.** Note that the space positions will differ from one path to another. In this case, estimation (29) is slightly modified for non-uniform grid to compute the experimental variogram of $Y$ from the logarithm of $G_t(\cdot)$. The spatial average is computed on all pairs of points whose distances are between $h$ and $h+\delta h$. Then, MOM can be extended for non-equidistant spatial inspections to estimate both spatial and temporal parameters.

Moreover, the proposed variogram of such fields depends only on the distance between positions. And finally, the MOM given by the least square method in
(31) can be used for some non-stationary RFs Y, such that a trend-stationary field with a linear trend or with stationary increments.

245

The temporal parameters are the parameters of the process $X_t$, i.e. the shape function $\alpha(t)$ defined in (28) and the deterministic contribution $\eta = e^\mu$ of the scale RF $\beta$. Since we assume that the power $b$ is known, $a$ and $\eta$ can be estimated. Note that $\eta$ contains the contribution of the spatial mean of the random field Y. From now, $\{G_{t_j}^i(z_l)\}$ is assumed to be a sequence of independent and identically distributed (i.i.d.) simulations of the degradation model described in the previous section. Each $i$-th degradation process is observed at time $t_j$ among $N_t$ times and on location $z_l$ among $N_z$ locations.

250

### Case of time-stationary model

255

Let consider the time-stationary case where $\alpha(t) = at$. Let $\zeta_{j,l}^i := \log(\delta G_{t_j}^i(z_l))$ be the logarithm of the increments of $G_t(z)$. Using the scaling property, each variable writes the following equation,

$$ \zeta_{j,l}^i := \log(G_{t_j}^i(z_l) - G_{t_{j-1}}^i(z_l)) = \log(X_{t_j}^i - X_{t_{j-1}}^i) + Y^i(z_l), \quad (33) $$

where $Y^i(z_l)$ is the $i$-th sample of Y at position $z_l$. Since the process $X_t$ and Y are independent, the first two moments of $\zeta_{j,l}^i$ write from equation (3):

$$ m_1 := \mathbb{E}[\zeta_{j,l}^i] = \psi(a\tau) - \log(\eta) \quad (34) $$

$$ m_2 := \mathrm{var}[\zeta_{j,l}^i] = \psi_1(a\tau) + \sigma^2, \quad (35) $$

where $\tau := (t_j - t_{j-1})$ is constant. The increments $(X_{t_j}^i - X_{t_{j-1}}^i)_{j=1}^{N_t}$ are independent and identically distributed. So, from $MN_t$ realizations of these increments combined with a spatial average, $m_1$ and $m_2$ can be estimated by;

260

$$ \hat{m}_1 := (MN_z N_t)^{-1} \sum_{i=1}^{M} \sum_{j=1}^{N_t} \sum_{l=1}^{N_z} \zeta_{j,l}^i, \quad (36) $$

$$ \hat{m}_2 := (MN_z N_t)^{-1} \sum_{i=1}^{M} \sum_{j=1}^{N_t} \sum_{l=1}^{N_z} (\zeta_{j,l}^i)^2 - \hat{m}_1^2 \quad (37) $$
Hence, an estimation of parameters a and η is given by:

$$ \begin{pmatrix} a \\ \eta \end{pmatrix} = f^{-1} \begin{pmatrix} \hat{m}_1 \\ \hat{m}_2 \end{pmatrix}, \qquad (38) $$

where f is given by $f\begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \psi(u\tau) - \log(v) \\ \psi_1(u\tau) + \hat{\sigma}^2 \end{pmatrix}$, with $\tau = (t_j) - (t_{j-1})$ and $\hat{\sigma}^2$ the estimate of $\sigma^2$ given by (31).

### Case of non-stationary model in time

Let consider now $\alpha(t) = at^b$ where $b \neq 1$ is known. The non-stationarity of GP can be transformed to a stationary GP by performing a monotonic transformation of the operational time [12]. However, the transformed inspection times are not equidistant to perform estimates given in (38). Let define the transformed times by $\nu_j = t_j^b - t_{j-1}^b$, and introduce the increments of ith sample of $G_t^i(z_l)$ denoted $D_{j,l}^i := \delta G_{t_j}^i(z_l)$. These increments are conditionally gamma distributed and conditionally independent. Therefore, according to [12], the following estimates hold true conditionally to the parameter $\beta$,

$$ \frac{a}{\beta} \approx \left. \frac{\sum_{j=1}^{N_t} D_{j,l}^i}{\sum_{j=1}^{N_t} \nu_j} \right|_{\beta} \qquad (39) $$

$$ \frac{a}{\beta^2} \left( \sum_{j=1}^{N_t} \nu_j - \frac{\sum_{j=1}^{N_t} \nu_j^2}{\sum_{j=1}^{N_t} \nu_j} \right) \approx \left. \sum_{j=1}^{N_t} \left( D_{j,l}^i - \frac{\sum_{j=1}^{N_t} D_{j,l}^i}{\sum_{j=1}^{N_t} \nu_j} \nu_j \right)^2 \right|_{\beta} \qquad (40) $$

By considering the empirical mean on the right side, knowing that $E[\beta^{-1}] = \eta^{-1}e^{\sigma^2/2}$ and $E[\beta^{-2}] = \eta^{-2}e^{2\sigma^2}$; estimation of a and $\eta$ are provided as solution of the following system,

$$ \frac{\hat{a}}{\hat{\eta}} = \frac{e^{-\hat{\sigma}^2/2}}{MN_z t_{N_t}^b} \sum_{i=1}^{M} \sum_{l=1}^{N_z} \sum_{j=1}^{N_t} D_{j,l}^i \qquad (41) $$

$$ \frac{\hat{a}}{\hat{\eta}^2} = \frac{t_{N_t}^b e^{-2\hat{\sigma}^2}}{MN_z \left( t_{N_t}^{2b} - \sum_{j=1}^{N_t} \nu_j^2 \right)} \sum_{i=1}^{M} \sum_{l=1}^{N_z} \sum_{j=1}^{N_t} \left( D_{j,l}^i - \frac{\sum_{j=1}^{N_t} D_{j,l}^i}{t_{N_t}^b} \nu_j \right)^2 \qquad (42) $$

where $\hat{\sigma}^2$ is the estimate of $\sigma^2$ given by (31).
## 4.2. Maximum Likelihood (MLM) and Pseudo Maximum Likelihood (PML)

The classical Maximum-Likelihood Estimate for the spatio-temporal field is performed when the RF is Gaussian or log-normal distributed. The model $G_t(\cdot)$ is given by the product of two independent process, the Gamma process and the spatial log-normal field. Because it is not obvious to compute the likelihood of any spatio-temporal increments of $G$, we introduce a method which considers only the model at the final time of inspection. It consists of maximizing the likelihood of the observed deterioration ($G_T(z_1), \dots, G_T(z_{N_z})$) spatially distributed at the final observed time $T$ which contains the most information. The probability density function of ($G_T(z_1), \dots, G_T(z_N)$) given a realization of $X_T = x$ is defined by:

$$ \varrho(v|x) = \phi(v) \exp \left\{ -\frac{1}{2\sigma^2} \left(\log(v/x) - \mu\right)' R^{-1} \left(\log(v/x) - \mu\right) \right\}, \quad (43) $$

where $\phi(v) := \prod_{l=1}^{N_z} \frac{1}{v_l^{-1} \sqrt{(2\pi)^{N_z} \sigma^{N_z} |R|}}$, $|R|$ the determinant of $R$ the $R$ is the correlation matrix of the Gaussian vector ($Y(z_1), \dots, Y(z_{N_z})$), $v$ a real vector with length $N_z$ and $\log(v)$ the logarithm of each component of $v$. Then, the joint probability function of ($G_T(z_1), \dots, G_T(z_N)$) writes:

$$ \varrho(v) = \int_{0}^{\infty} \ell(v|x) \frac{x^{\alpha(t)-1} e^{-x}}{\Gamma(\alpha(t))} dx, \quad (44) $$

The likelihood $\varrho(v)$ can be accurately computed through the Gauss-Laguerre quadrature or the Gauss-Hermite quadrature after transforming the integral on $\mathbb{R}$ with the transformation $y = \log(z)$.

Therefore, from $M$ copies of the vector degradation ($G_T^i(z_1), \dots, G_T^i(z_{N_z})$), for $i = 1 \dots, M$, the likelihood of these observed data is given by

$$ L(V) = \prod_{i=1}^{M} \varrho(v^i); \quad (45) $$

<footnote n="280">280 where $V$ is the matrix of $M$ vectors degradation $v^i = (G_T^i(z_1), \dots, G_T^i(z_{N_z}))$. Therefore, the maximum-likelihood estimates of the parameters $a, \nu, l_c, \sigma^2$ and $\mu$ are provided by maximizing the logarithm of the likelihood $\log(L(V))$.</footnote>
The likelihood in (45) depends on both spatial and temporal parameters. An acceptable estimate of those parameters is reached only with large independent copies of degradation samples. A bad estimation of $\alpha$ affects spatial parameters estimates. MOM provides spatial parameters independently of the temporal ones, this makes the method more flexible. However, it is efficient only for the stationary model.

Another approach which can be combined with MOM for non-stationary temporal variability is the pseudo maximum-likelihood method (PML). It consists of maximizing the likelihood of the increments $(\delta_1 G, \dots, \delta_{t_N} G)$ on a given fixed spatial position $z$. These increments $(\delta_{t_j} G = G_{t_j} - G_{t_{j-1}})_{j=1}^{N_t}$ are conditionally independent and their likelihood is given by the product of the marginal density of each increment, which is similarly computed as in (11),

$$ \ell(v) = \prod_{j=1}^{N_t} \frac{\eta^{\delta_{t_j} \alpha} v_j^{(\delta_{t_j} \alpha - 1)}}{\Gamma(\delta_{t_j} \alpha)} \int_{\mathbb{R}} \exp\left(-v_j \eta e^{\sigma y} + \delta_{t_j} \alpha \sigma y\right) \xi(y) dy, \quad (46) $$

where $\delta_{t_j} \alpha := \alpha(t_j) - \alpha(t_{j-1})$ and $v$ is the observed vector of the increment $(\delta_1 G, \dots, \delta_{t_N} G)$ on a given position $z$. When $b$ is unknown, we can use this latter approach of PML given by the maximum of (46) to estimate temporal parameters. The value of $\sigma$ comes from problem (31). Then, by using an appropriate approximation of each integral as given in (12), PML estimators of temporal parameters are given by maximizing $\log(\ell(v))$.

## 5. Numerical simulations

The deterioration model is assumed to be an observable process in space and time with limited information. However, such kinds of database which vary in space and time are not available in the literature. Therefore, a synthetic discrete degradation model is constructed through MC simulations. The purpose of this section is to validate the inference methodology and the approximation of the quantities of interest which depend on the estimated parameters.

We consider $M = 1, 10, 100$ sample paths of the field $(G^j)_{j=1}^M$, each trajectory $G^j$ being simulated at $N_t$ equidistant periods on the interval time $[0, 30]$ (in
years) and $N_z$ equidistant locations in one- and two-dimensional space. The model in time (its temporal increments are time-independent), i.e the shape parameter is linear $\alpha(t) = at$. The deterministic contribution of the scale RF is given by $\eta = e^{\mu}$ where $\mu = 2/3$, $\sigma^2 = 0.6$, $l_c = 1$ (correlation length) and $\nu = 2$ (smoothness parameter of Y).

### 5.1. Simulation of discrete degradation model

In this subsection, the steps for generating a path of the degradation model $G_t(\cdot)$ on a uniform discrete time-space partition is presented. Let consider $N_t$ periods of uniform discrete-time partition $t_0 = 0, t_1, \dots, t_{N_t}$ on the time interval $[0, T]$ and $N_z$ equidistant positions $z_1, \dots, z_{N_z}$ of the spatial interval.

Firstly, let generate the GRF on the grid $z_1, \dots, z_{N_z}$ using the circulant em-
bedding matrix approach defined in Section 2, see ([27]) for more details of the
method for two or three spatial dimensions. The Gaussian vector $(Y(z_1), \dots, Y(z_{N_z}))$
has an exact covariance matrix given by the continuous covariance model of Y.
Secondly, the degradation model is progressively generated in each time for all
positions. For each $z_j$ the increments of the model are independently generated
with identical gamma distribution. This naïve simulation approach is called
Spatial Gamma Sequential Sampling (SGSS). Note that in the one-dimensional
variability, the discrete model $G$ is given by a matrix with size $(N_t + 1) \times N_z$.
A pseudo-code to simulate one path of the model in one-dimensional variability
is given:

$G_0(\cdot) = 0;$

$\tau = \frac{T}{N_t+1};$

For $i = 1$ to $N_z$

For $j = 1$ to $N_t$

Generate $Q \sim Ga(a\tau, e^{\sigma Y(z_i)+\mu})$;

$G_{t_j}(e^{\sigma Y(z_i)+\mu}) = G_{t_{j-1}}(e^{\sigma Y(z_i)+\mu}) + Q;$

Next $j$

Next $i$
This algorithm is running *M*-times to get *M* samples of discrete time-space paths of the degradation process. Similarly, the pseudo-code for a two-dimensional variability case can be performed where the discrete model *G* is given by three-dimensional array.

## 5.2. One-dimensional variability

The Gaussian random field is defined on the interval $[0, L]$ where $L = 100$ and discretized in $N_z$ equidistant spatial positions. From now and for a shake of simplicity, a fixed value of the smoothness parameter $\nu = 2$ is selected (i.e the paths of $G$ are differentiable in quadratic norm). Figure 1 (left) gives one realization of $G_t(z)$. The surface is discontinuous with respect to $z$ because infinite jumps of the Gamma process $X_t$. Further, the surface shows a strong spatial variability since $l_c \ll L$. The graphic compares (right) the experimental and the exact variograms of the GRF $Y$. The experimental variogram is computed using the spatial trajectory of the logarithm $\log(G_t(z))$ at time $t = 30$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Left: example of one realization of $G_t(\cdot)$, right: variogram ($l_c = 1$, $\nu = 2$, $\sigma^2 = 0.6$).</figcaption></figure>

### Method of moments (MOM) Step 1

The first step of the method of moment consists in estimating spatial parameters, variance $\sigma^2$ and correlation length $l_c$. The quality of the estimation is measured by the mean absolute error. This error is given by the average of the absolute differences between the exact parameter and 10 estimated values
calculated across MOM. Table 1 gives estimates of those parameters obtained by MOM where we minimize equation (31) with $\nu = 2$ and three values of $M$. (Mean absolute errors in brackets).

<table><caption>Table 1: Parameters estimation of one dimensional spatial variability by MOM.</caption><thead><tr><th colspan="3">N<sub>z</sub> = 40</th><th colspan="3">N<sub>z</sub> = 100</th></tr><tr><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.714(0.114)</td><td>0.21(0.791)</td><td>1</td><td>0.546(0.053)</td><td>1.186(0.186)</td></tr><tr><td>10</td><td>0.5720 (0.028)</td><td>0.81(0.19)</td><td>10</td><td>0.584(0.015)</td><td>0.934(0.065)</td></tr><tr><td>100</td><td>0.6048(0.0048)</td><td>0.890(0.1024)</td><td>100</td><td>0.6009(0.009)</td><td>1.024(0.024)</td></tr></tbody></table>

From Table 1, the estimation of the correlation length $l_c$ depends strongly on the number of spatial positions $N_z$, obviously because $l_c$ is very small with respect to the length $L = 100$. In contrast to $l_c$, the estimate of the variance $\sigma^2$ is largely acceptable even with small number of positions for one realization of $G_t(\cdot)$ ($M = 1, N_z = 40$). However, the estimation of quantities of interest depends strongly on the variance $\sigma^2$, so an accurate estimate of $\sigma^2$ is needed to forecast reliable predictions.

### Method of moments Step 2

Once spatial parameters are estimated, the second step of the MOM consists in estimating the temporal parameters given by (38). An estimation of $\sigma^2$ is inserted in (38) for each case of $M$ and $N_z$. Table 2 summarizes estimations of temporal parameters $a$ and $\eta$. Results show that their accuracy depends strongly on the total inspection times $N_t$ and on the total number of positions $N_z$. An acceptable accuracy is reached for small number of realizations $M$ ($M = 10$) when $N_t$ and $N_z$ are significantly large ($N_t = 60, N_z = 40$). In particular, estimate of $\eta$ depends significantly on the number of locations $N_z$ since it contains the stochastic contribution of the random field $Y$.

Figures 2 illustrate the convergence of the MOM for temporal parameters $\alpha$ and $\eta$, with a simulated ($M \times N_z \times N_t$) database. The mean absolute error on $a$ and $\eta$ is a function of the product $N_z \times N_t$ where $N_z = 10, 40, 80, 100$
<table><caption>Table 2: Parameters estimation of temporal variability using MOM</caption><thead><tr><th colspan="3">N<sub>z</sub> = 40, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z</sub> = 40, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>0.79(0.21)</td><td>1.487(0.46)</td><td>1</td><td>0.902(0.098)</td><td>1.617(0.33)</td></tr><tr><td>10</td><td>0.947(0.053)</td><td>2.157(0.21)</td><td>10</td><td>1.027(0.027)</td><td>2.087(0.14)</td></tr><tr><td>100</td><td>1.012(0.012)</td><td>1.986(0.038)</td><td>100</td><td>1.01(0.009)</td><td>1.97(0.022)</td></tr></tbody><thead><tr><th colspan="3">N<sub>z</sub> = 100, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z</sub> = 100, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>0.817(0.18)</td><td>2.377(0.43)</td><td>1</td><td>1.01(0.1)</td><td>1.567(0.38)</td></tr><tr><td>10</td><td>1.057(0.057)</td><td>1.837(0.11)</td><td>10</td><td>1.029(0.029)</td><td>1.817(0.13)</td></tr><tr><td>100</td><td>1.012(0.012)</td><td>2.015(0.068)</td><td>100</td><td>1.01(0.009)</td><td>2.973(0.026)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: Mean absolute error on a and η</figcaption></figure>

<footnote n="375">375 and $N_t = 10, 30, 40, 60$. Let remark a strong improvement of the quality of the estimation for $M > 1$, $N_z > 10$ and $N_t > 10$. Estimate accuracy of a is nearly independent of the spatial positions unlike $η$ which it depends on $N_z$ and strongly on $N_t$. Let note that in the range $4800 < N_z \cdot N_t < 6000$ the error is monotonically decreasing since $N_t = 60$ and $N_t \le 40$ outside of this range. An acceptable accuracy is reached with one component $M = 1$ when $N_t$ and $N_z$ are quite large ($N_t = 60, N_z = 40$).</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Errors on *a* and *η* with respect to $N_z$ and $N_t$</figcaption></figure>

Figures 3 sketches the means of absolute error on *a* and *η* with respect to $N_t$ or $N_z$, which is computed from 500 estimates through MC simulations. The temporal estimates are conducted by fixing $N_z = 40$ and $N_t = (5 : 5 : 100)$. Similarly to the spatial estimation, let fix $N_t = 30$ and $N_z = (5 : 5 : 100)$. The figures shows that the accuracy on *η* depends on $N_z$ and strongly on $N_t$, unlike of the parameter *a* which is nearly independent of the spatial positions.

### Pseudo-maximum-likelihood (PML) estimates

Table 3 gives estimates of *a* and *η* using PML which consists of maximizing (46), the logarithm of the marginal density of the increments of $G_t(\tilde{z})$ at an arbitrary fixed location $\tilde{z}$. Here we choose $N_z = 40$ and the increments of $G_t(\tilde{z})$ are considered at the position $\tilde{z} = \frac{L}{2}$. The likelihood in (46) depends on $\sigma^2$, so let choose its estimate according to *M* and $N_z$ in Table 1. For each case, estimations show that the accuracy of PML method depends strongly on *M* independent copies of $N_t$ increments of $G_t$. Further, the PML method provides several estimates where each one depends on the position $\tilde{z}$ of the increments. One may consider the mean of the estimates computed at different positions.

The method of moments looks to be more attractive for inference procedure since it gives the parameter estimates on two separable stages. First, the spatial
<table><caption>Table 3: Parameters estimate of temporal variability by PML</caption><thead><tr><th colspan="3">N<sub>z</sub> = 40, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z</sub> = 40, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.54(0.54)</td><td>1.009(0.94)</td><td>1</td><td>0.84(0.16)</td><td>2.75(0.8)</td></tr><tr><td>10</td><td>1.269(0.269)</td><td>2.41(0.465)</td><td>10</td><td>1.05(0.05)</td><td>2.09(0.14)</td></tr><tr><td>100</td><td>1.163(0.163)</td><td>2.32(0.381)</td><td>100</td><td>0.965(0.034)</td><td>1.794(0.153)</td></tr></tbody></table>

parameter estimates are provided then the temporal ones. Further, by using the spatial data at both stages, the spatial uncertainties are reduced yielding to an improvement of the inference by exploiting inspection results.

<figure><img src="image_6.png" /><figcaption>Figure 4: Left: estimate of the failure time distribution $F_T$, right: marginal density by sampling and quadrature approach.</figcaption></figure>

Once the model parameter estimation is performed, a $10^3$ MC simulations of the discretized $G$ is done at $N_t = 100$ times and $N_z = 100$ positions in order to estimate the marginal cumulative distribution $F_T$ of the failure time $T_F$ by sampling method (17), Figure 4 (left) plots the cdf $F_T$ computed by sampling approach with exact and estimated parameters and by the quadrature rule (16) with $m = 30$ Gaussian knots. Figure 4 (right) compares marginal density of $G_t$ by the quadrature method (12) at several times $t$ with a convergence criterion $\epsilon \approx 10^{-3}$. Let note that the required order of the quadrature rule increases with time $t$. Figure 5 (left) sketches the effect of conditioning the current state on the
<figure><img src="image_7.png" /><figcaption>Figure 5: Left: distribution of predictive time of inspection to failure, right: reliability evaluation and remaining lifetime function.</figcaption></figure>

failure prediction at the current time of inspection $t_i = 10$ for several measured degradation levels from 0.75 to 1.6. The survival function depends on the current state unlike the reliability function ($g_{t_i}$ is unknown, blue curve) which is computed given that the component is in service. The curves highlight that, depending on the current state, the reliability function (blue curve) gives an underestimate or overestimate of the time to failure compared with the survival function. Hence, this latter quantity gives more accurate prediction about time to failure than the reliability function.

Figure 5 (right) draws several failure curves of the remaining lifetime $F_{RL_{t_i}}(\tau)$ with several times of inspection $t_i$ with the same degradation state $g_{t_i} = 1.25$. These figures highlight an obvious result: later is the observation of a given degradation state, longer is the residual lifetime.

<table><caption>Table 4: Parameter estimates by MOM for one sample path $G_t$ and $X_t$ without spatial data and standard deviation of each estimate is in the brackets.</caption><thead><tr><th>parameters</th><th>a</th><th>η</th><th>σ</th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>G<sub>t</sub></td><td>1.32 (0.31)</td><td>2.44 (0.84)</td><td>0.62 (0.21)</td><td>2.37 (1.41)</td></tr><tr><td>X<sub>t</sub></td><td>0.74 (0.36)</td><td>0.34 (0.28)</td><td></td><td></td></tr></tbody></table>

In order to illustrate the impact of dealing with the spatial variability in
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of Remaining lifetime function, $g_{max} = 25$, $g_{t_i} = 6$, $t_i = 2$.</figcaption></figure>

the degradation analysis, let consider one sample path ($M = 1$) of the model $G_t(z)$ discretized on $N_z = 30$ positions and $N_t = 30$ periods. Further, for each instant $t_i$, let $X_{t_i}$ the important observed states from $G_{t_i}(z)$ (maximum value) among all positions. Thus, $(X_{t_i})_{i=1}^N$ can be modeled with a sample path of the standard GP $Ga(\alpha(t), \eta)$. The classical MOM is used to compute parameters of $X_t$ and the MOM for $G_t(z)$. Both are illustrated in Table 4.

Figure 6 compares the cdf of the remaining lifetime for both models $G_t$ and $X_t$ (without spatial variability). This cdf is computed using the estimate parameters given in Table 4, critical level $g_{max} = 25$ and observed degradation $g_{t_i} = 6$ at time $t_i = 2$. Let remind that the cdf $F_{RL_{X_t}}(\tau)$ of the remaining lifetime function for $X_t$ is given by the formula $F_{RL_{X_t}}(\tau) = 1 - P(X_{t+\tau} - X_t < g_{max} - x_t)$, where $x_t$ is the state at time $t$. From the result of Table 4 and Figure 6, we can conclude that considering the spatial variability yields to an accurate model and more predictions for reliability analysis.

### 5.3. Two dimensions variability

The GRF $Y(z)$ is herein defined on the rectangle $[0, L_1] \times [0, L_2]$ where $L_1 = 100, L_2 = 40$ discretized on $N = N_{z_1} \times N_{z_2}$ equidistant spatial positions with $N_{z_1} = 40, 100$ and $N_{z_2} = 20, 40$. The field $Y(z)$ is assumed to be isotropic where $l_c = 1$ is the same in both directions. Figure 7 (left) gives four slices from a realization of $G_t(z)$ at times $t = 0, 10, 20, 30$; and plots (right) the surface of
the experimental variogram of Y which is computed using the spatial trajectory $log(G_t(z))$ at time $t = 30$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Left: one realization of $G_t(\cdot)$, right: experimental variogram of Y.</figcaption></figure>

MOM. Step 1: estimation of spatial parameters

Table 5 summarizes estimates of $\sigma^2$ and $l_c$ obtained by the MOM with the minimization of the standard least square in (31). The smoothness parameter is fixed $\nu = 2$ for three values of $M$ ($M = 1, 10, 100$), (mean absolute errors in brackets).

<table><caption>N<sub>z1</sub> = 40, N<sub>z2</sub> = 20</caption><thead><tr><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.638(0.038)</td><td>0.913(0.087)</td></tr><tr><td>10</td><td>0.583(0.017)</td><td>1.06(0.06)</td></tr><tr><td>100</td><td>0.6008(0.008)</td><td>1.016(0.016)</td></tr></tbody></table>

<table><caption>Table 5: Parameters estimation of two-dimensional spatial variability by MOM.</caption><caption>N<sub>z1</sub> = 100, N<sub>z2</sub> = 40</caption><thead><tr><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.58(0.02)</td><td>0.949(0.051)</td></tr><tr><td>10</td><td>0.606(0.006)</td><td>1.007(0.007)</td></tr><tr><td>100</td><td>0.597(0.003)</td><td>1.0028(0.0028)</td></tr></tbody></table>

In two-dimensional variability, parameters estimation (Table 5) is more accurate than the one-dimensional estimation (even with poor realizations of $G_t$). This accuracy is explained by the use of large spatial locations in the inference $(N_{z1} \times N_{z2})$.
Table 6 provides the estimations of the spatial parameters with total number positions $N = 100$ where $N_{z_1} = 20$, $N_{z_2} = 5$ and $N_{z_1} = 50$, $N_{z_2} = 2$. Results show that the accuracy is closed to the one-dimensional case, in particular estimation using the discretization $N_{z_1} = 50$, $N_{z_2} = 2$. Therefore, when the model is isotropic, estimations using data on a single direction lead to the same accuracy as one-dimensional model.

<table><thead><tr><th colspan="3">N<sub>z1</sub> = 20, N<sub>z2</sub> = 5</th></tr><tr><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.714(0.114)</td><td>1.967(0.967)</td></tr><tr><td>10</td><td>0.669(0.069)</td><td>1.685(0.685)</td></tr><tr><td>100</td><td>0.619(0.019)</td><td>1.518(0.518)</td></tr></tbody></table>

<table><caption>Table 6: Parameters estimation of two-dimensional spatial variability by MOM.</caption><thead><tr><th colspan="3">N<sub>z1</sub> = 50, N<sub>z2</sub> = 2</th></tr><tr><th>M</th><th>σ<sup>2</sup></th><th>l<sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0.662(0.062)</td><td>1.776(0.776)</td></tr><tr><td>10</td><td>0.565(0.035)</td><td>0.759(0.241)</td></tr><tr><td>100</td><td>0.592(0.007)</td><td>0.942(0.058)</td></tr></tbody></table>

MOM. Step 2: estimation of temporal parameters

Table 7 summarizes estimations of $a$ and $\eta$. Results show that accuracy depends more on $N_t$ than $N_{z_1} \times N_{z_2}$. An acceptable accuracy is reached with few realizations $M$ when $N_t$ and $N_{z_1}, N_{z_2}$ are significantly large. Estimation of $\eta$ depends strongly on the number of positions since each position contains a stochastic contribution of the RF $Y$ given by its mean $\mu$.

In Table 8, as in one-dimensional case, we estimate temporal parameters with total number positions $N = 100$ where $N_{z_1} = 20$ and $N_{z_2} = 5$ or ($N_{z_1} = 50$ and $N_{z_2} = 2$). Estimation provides nearly the same accuracy as in the one-dimensional case (Table 2), since it depends strongly on $N_t$ and significantly on the total number of positions. Therefore, these results combined with Table 6 suggest that when the model is isotropic, the spatio-temporal variability can be appropriately modeled by one-dimensional model.

Figure 8 compares the failure distribution $F_T$ which is computed using realizations of the field $G_t(z)$ with exact and estimate parameters. Figure 8 (left) uses the estimation obtained with $N_{z_1} = 40, N_{z_2} = 20$ and $M = 1, 10$ components. The distribution $F_T$ is computed using 1000 MC simulations of the
<table><caption>Table 7: Parameters estimation of the temporal variability using MOM.</caption><thead><tr><th colspan="3">N<sub>z1</sub> = 40, N<sub>z2</sub> = 20, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z1</sub> = 40, N<sub>z2</sub> = 20, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.613(0.613)</td><td>2.660(0.71)</td><td>1</td><td>1.224(0.224)</td><td>2.447(0.49)</td></tr><tr><td>10</td><td>0.934(0.06)</td><td>1.833(0.134)</td><td>10</td><td>1.034(0.034)</td><td>2.084(0.136)</td></tr><tr><td>100</td><td>0.991(0.009)</td><td>1.968(0.02)</td><td>100</td><td>1.011(0.011)</td><td>1.971(0.023)</td></tr></tbody><thead><tr><th colspan="3">N<sub>z1</sub> = 100, N<sub>z2</sub> = 40, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z1</sub> = 100, N<sub>z2</sub> = 40, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>0.801(0.198)</td><td>2.184(0.236)</td><td>1</td><td>0.872(0.127)</td><td>1.722(0.225)</td></tr><tr><td>10</td><td>1.072(0.07)</td><td>2.123(0.175)</td><td>10</td><td>1.018(0.018)</td><td>2.068(0.121)</td></tr><tr><td>100</td><td>1.09(0.009)</td><td>2.024(0.07)</td><td>100</td><td>1.014(0.014)</td><td>2.017(0.069)</td></tr></tbody></table>

<table><caption>Table 8: Parameters estimation of the temporal variability using MOM.</caption><thead><tr><th colspan="3">N<sub>z1</sub> = 20, N<sub>z2</sub> = 5, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z1</sub> = 20, N<sub>z2</sub> = 5, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.861(0.861)</td><td>1.146(0.801)</td><td>1</td><td>1.287(0.287)</td><td>2.503(0.555)</td></tr><tr><td>10</td><td>0.915(0.085)</td><td>1.578(0.369)</td><td>10</td><td>1.087(0.087)</td><td>2.368(0.420)</td></tr><tr><td>100</td><td>0.968(0.031)</td><td>1.858(0.089)</td><td>100</td><td>0.962(0.037)</td><td>2.031(0.084)</td></tr></tbody><thead><tr><th colspan="3">N<sub>z1</sub> = 50, N<sub>z2</sub> = 2, N<sub>t</sub> = 30</th><th colspan="3">N<sub>z1</sub> = 50, N<sub>z2</sub> = 2, N<sub>t</sub> = 60</th></tr><tr><th>M</th><th>a</th><th>η</th><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.651(0.651)</td><td>2.734(0.787)</td><td>1</td><td>0.814(0.186)</td><td>1.563(0.384)</td></tr><tr><td>10</td><td>0.913(0.086)</td><td>1.836(0.112)</td><td>10</td><td>0.971(0.029)</td><td>2.071(0.124)</td></tr><tr><td>100</td><td>0.968(0.031)</td><td>1.871(0.076)</td><td>100</td><td>1.029(0.029)</td><td>2.034(0.086)</td></tr></tbody></table>

field $G_t(z)$. Figure 8 (right) uses 500 MC simulations of $G_t(z)$ with parameters estimate given in the case $N_{z1} = 100, N_{z2} = 40$ and $M = 1, 10$. This latter shows an accurate estimation of the failure time distribution $F_T$, since it uses more positions and accurate parameters estimation.

The PML method consists of maximizing (46), the marginal density of the
<figure><img src="image_10.png" /><figcaption>Figure 8: Failure time distribution $F_T$. Left 1000 MC simulations, $N_{z1} = 40$, $N_{z2} = 20$. Right $M = 500$ MC simulations, $N_{z1} = 100$, $N_{z2} = 40$</figcaption></figure>

increments of $G_t(\tilde{z})$ in an arbitrary position $\tilde{z}$. Table 9 gives estimates of $a$ and $\eta$ by PML, where $G_t(\cdot)$ is simulated with $N_{z1} = 40, N_{z2} = 20$. The increments of $G_t(\tilde{z})$ are computed in the middle position of the physical domain. *sigma* in the likelihood (46) is replaced by its estimations for each case of $M, N_{z1}$ and $N_{z2}$ in (Table 5). Results of these estimates show that the PML estimators depend strongly on $M$ independent copies of $N_t$ increments of $G_t(z)$.

<table><caption>N<sub>z1</sub> = 40, N<sub>z2</sub> = 20, N<sub>t</sub> = 30</caption><thead><tr><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.449(0.449)</td><td>2.411(0.463)</td></tr><tr><td>10</td><td>1.045(0.045)</td><td>2.117(0.169)</td></tr><tr><td>100</td><td>0.972(0.027)</td><td>2.065(0.118)</td></tr></tbody></table>

<table><caption>Table 9: Parameters estimation of temporal variability by PML</caption><caption>N<sub>z1</sub> = 40, N<sub>z2</sub> = 20, N<sub>t</sub> = 60</caption><thead><tr><th>M</th><th>a</th><th>η</th></tr></thead><tbody><tr><td>1</td><td>1.283(0.283)</td><td>2.544(0.597)</td></tr><tr><td>10</td><td>0.968(0.032)</td><td>1.797(0.15)</td></tr><tr><td>100</td><td>1.011(0.011)</td><td>2.007(0.06)</td></tr></tbody></table>

## Discussion

The comparison between estimates given by MOM and by PML method leads to highlight the interest of MOM. The PML provides only the temporal parameters by using the increments of the model in any fixed position, yielding to a set of acceptable estimations. MOM makes the inference procedure on two stages.
First, it determines the spatial parameters and then the temporal ones. Under the property of spatial ergodicity, the use by MOM of the spatial data in both stages leads to improve the inference accuracy by exploiting inspection results. When the temporal variability is characterized by a non-homogeneous GP, MOM requires some additional transformations in the time step to obtain a stationary process. This transformation is not clear when the non-linear shape function is defined by unknown parameters. Thus, MOM and PML can be combined for more practical inference procedure in such case. This combination consists in, first, estimating the spatial parameters with MOM, second, computing the temporal parameters for the non-homogeneous GP with PML.

## 6. Conclusion

We have developed a spatio-temporal degradation model that tackles both the inert spatial uncertainty and heterogeneity across structural component. It is based on the classical gamma process and non-negative spatial random field. The temporal paths of the process are monotonic with conditionally independent increments, the random field scale follows a log-normal distribution. The marginal probability density and cumulative function of the model has found to have a transcendental form and the Gauss quadrature formula was applied to compute these probabilities. The quantities of interest in reliability analysis and maintenance; namely the distribution of failure time and the distribution of remaining useful lifetime are computed. A comparison between the method of moments and a maximum pseudo-likelihood method is carried out to derive the parameters of the model. The advantage of the model is highlighted with Monte Carlo simulations. It shows how incorporating spatial variability both reduces uncertainties and improves accuracy of the inference by exploiting with better manner the spatial data. An interesting extension of the model that incorporates the spatial variability can consider the multi-variate modeling based on the state dependent Gamma process. Such spatio-temporal field is not sep-
arable when the spatial aleatory is considered in the shape function, it seems
computationally more difficult and is still a challenge.

## Acknowledgements

The work reported in this article was supported by the Region Pays de la
Loire de France through the regional project SI3M (Strategies of Identification
of Meta-models for Structural Maintenance). The authors would like to thank
the region for the support.

## REFERENCES

[1] PAPAKONSTANTINOU KG, SHINOZUKA M. *Probabilistic model for steel corrosion in reinforced concrete structures of large dimensions considering crack effects*. J Eng Struct;57:306-326, 2013

[2] PAPAKONSTANTINOU KG, SHINOZUKA M. *Spatial stochastic direct and inverse analysis for the extent of damage in deteriorated RC structures*. J Comput Struct, 128, 286-296, 2013.

[3] VU K, STEWART M. *Predicting the likelihood and extent of reinforced concrete corrosion-induced cracking*. J. Struct. Eng., Volume 131, Issue 11, 1681-1689; November 2005.

[4] STEWART M., MULLARD JA. *Spatial time-dependent reliability analysis of corrosion damage and the timing of first repair for RC structures*. Eng Struct, 29:1457-64, 2007.

[5] QINGHUI SUO, STEWART M. *Corrosion cracking prediction updating of deteriorating RC structures using inspection information*. Reliability Engineering and System Safety 94, 1340-1348, 2009.

[6] ABDEL-HAMEED M. *A Gamma wear process*. IEEE Transactions on Reliability, 24(2), 152-153, 1975.
[7] ADLER ROBERT.J, *The Geometry of Random Fields*, Classics in Applids Mathematics, SIAM 1981.

[8] BONNET S., SCHOEFS F., SALTA M., *Statistical study and probabilistic modelling of error when building chloride profiles for reliability assessment*, European Journal of Environmental and Civil Engineering, pp 1-16; 2017

[9] BAGDONAVIČIUS V, NIKULIN MS. *Estimation in degradation models with explanatory variables*. Lifetime Data Anal, 7(1), 2001. 85-103.

[10] BEHROUZ SHAFEI, ALICE ALIPOUR. *Application of large-scale non-Gaussian stochastic fields for the study of corrosion-induced structural deterioration*, Engineering Structures, 88 262-276, 2015.

[11] CHILÈS, J.-P., AND DELFINER, P., *Geostatistics: Modeling Spatial Uncertainty*, New York: Wiley, 1999.

[12] ÇINLAR E, BAŽANT ZP, OSMAN E. *Stochastic process for extrapolating concrete creep*. J Eng Mech Div, 103(EM6):1069-88; 1977.

[13] CRESSIE, N. *Statistics for Spatial Data*, Wiley, New York, 1993.

[14] DIETRICH, C. ET G. NEWSAM. *Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix*. SIAM J. SCI. COMPUT. Vol. 18, No. 4, pp. 1088-1107, July 1997.

[15] El Hajj, Boutros, Bruno Castanier, Franck Schoefs, and Thomas Yeung. *A Maintenance-Oriented Degradation Model for a Reinforced Concrete Structure Subject to Cracking*. Journal of Risk and Reliability, 1-13, 2015.

[16] EMILIO BASTIDAS-ARTEAGA, FRANCK SCHOEFS. *Stochastic Improvement of Inspection and Maintenance of Corroding Reinforced Concrete Structures Placed in Unsaturated Environments*. Eng. Struct. 41, 50-62, 2012.

[17] JIA, G., TABANDEH, A., & GARDONI, P. *Lifecycle analysis of engineering systems: Modeling deterioration, instantaneous reliability, and resilience*.
In Risk and Reliability Analysis: Theory and Applications (pp. 465-494).
Springer International Publishing, 2017.

[18] KUMAR, R., CLINE, D. B., & GARDONI, P. *A stochastic framework to model deterioration in engineering systems*. Structural Safety, 53, 36-43; 2015.

[20] J. ERHEL, Z. MGHAZLI, M. OUMOUNI *A combined collocation and Monte Carlo method for advection-diffusion equation of a solute in random porous media* ESAIM-Proceedings, Vol. 45, p. 328-337, 2014.

[22] J. ERHEL, G. PICHOT, M. OUMOUNI *About a Generation of a Stationary Gaussian Field by Spectral method*. In preparation.

[24] F. CASCIATI, S. CASCIATI, L. FARAVELLI *A contribution to the modelling of human induced excitation on pedestrian bridges*. Structural Safety; 66: 51-61, 2017.

[25] GNEITING, T., GENTON, M. G., AND GUTTORP, P., *Geostatistical space-time models, stationarity, separability and full symmetry*, in Statistical Methods for Spatio-Temporal Systems, eds. B. Finkenstadt, L. Held, and V. Isham, Boca Raton: Chapman & Hall/CRC, pp. 151-175, 2007.
605

[26] LAWLESS J, CROWDER M. *Covariates and random effects in a gamma process model with application to degradation and failure*. Lifetime Data Anal; 10(3), 213-227, 2004.

610

[27] G. LORD, C. POWELL, T. SHARDLOW. *An Introduction to Computational Stochastic PDEs*. Cambridge Texts in Applied Mathematics. Cambridge University Press, New York, 2014.

615

[28] MADANAT, S., R. MISHANLANI, AND W.H.W. IBRAHIM, *Estimation of infrastructure transition probabilities from condition rating data*. Journal of Infrastructure Systems 19(2), 120-125, 1995.

[29] PARDO-IGUZQUIZA, E., AND DOWD, P, *Variance covariance matrix of the experimental variogram: assessing variogram uncertainty*, Mathematical Geology, 33, 397-419, 2001.

620

[30] SHINOZUKA, M. ET G. DEODATIS. *Simulation of stochastic processes by spectral representation*. App. Mech. Rev., 44(4), 191-204, 1991.

[31] SCHOEFS F., BASTIDAS-ARTEAGA E., TRAN T.V, *Optimal Embedded Sensor Placement for Spatial Variability Assessment of Stationary Random Fields*, Engineering Structures, 152, 35-44, 2017.

[32] SCHOEFS F, CLEMENT A, NOUY A. *Assessment of spatially dependent ROC curves for inspection of random fields of defects*. Structural Safety; 31:409-19, 2009.

625

[33] SI, XIAO SHENG, WENBIN WANG, CHANG HUA HU, AND DONG HUA ZHOU. *Remaining Useful Life Estimation - A Review on the Statistical Data Driven Approaches*. European Journal of Operational Research. Volume 213, Issue 1, 1-14, August 2011.

[34] SUDRET B. *Probabilistic models for the extent of damage in degrading reinforced concrete structure*. Reliab Eng Syst Saf V 93, Issue 3, 410-422, 2008.
[35] VAN NOORTWIJK, J.M. *A Survey of the Application of Gamma Processes in Maintenance*. Reliability Engineering & System Safety 94; 2-21, 2009

[36] A.M. YAGLOM, *Correlation Theory of Stationary and Related Random Functions*, Springer-Verlag, New York, 1987.