3. Spatial-temporal process
3.2 Bayesian estimation considering spatial correlation structure
In order to improve the estimation and to enable local estimation, utilization of Bayesian estimation including spatial correlation is proposed in this research. This approach uses prior information of the unknown parameters which characterize soil behavior, e.g. model parameters or soil properties, and the observation data, e.g. observed settlement or movement, from all observation points to search for the best estimates of the unknown parameters. The formulation consists of two statistical components, namely, the observation model and the prior information model. These two models will then be combined by Bayes’ theorem to obtain the solution.
3.2.1 Observation model
This model relates the observation data to the unknown parameters which are defined in a multivariate stochastic Gaussian field Ζ(x) = [z1(x), z2(x), ... , zP(x)]T where x is a spatial vector coordinate; P is total number of the unknown parameters; and zi(x) (for i = 1, 2, … , P) is a random function of an unknown parameter (e.g. soil or model parameter) at any location x in a specific domain. This paper proposes a method to identify the best estimator of Z for a discrete spatial point field, x1, x2, ... , xn, xn+1, ... , xm, which is defined as
[
1 2]
ˆ ˆ , ˆ , , ˆP T
Z = z z L z (3-1)
where
[
1 2]
ˆi ˆi( ) , ˆi( ) , , ˆi( )m ; 1, 2,...,
z = z x z x L z x i= P (3-2)
Suppose that a set of observations Yk (e.g. ground settlement) at the discrete time step k, i.e. k
= 1, 2, ... , K, has been obtained at n observation points x1, x2, ... , xn. Yk is defined as
[
k k k n]
Tk y x y x y x
Y = ( 1) , ( 2) , L , ( ) (3-3)
It should be noted that xn+1, xn+2, ... , xm are defined as any arbitrary points at which the unknown parameters are to be estimated, i.e. m-n interpolation points.
The general formulations of the observation model have been presented in several literatures, e.g. Hoshiya and Yoshida (1996), Honjo and Kashiwagi (1999), etc. Here it is assumed that the observation Yk is expressed as a linear function of Z with observation error of ε as follows:
k k
Y =M Z+ε (3-4)
where ε is the Gaussian observation error vector which is assumed to follow N(0, Vε). Vε is defined as a covariance matrix of ε where Vε = σε2·In,n. σε2 is the variance of the observation error and In,n is an n × n unit matrix. This implies that the observation errors are assumed to be spatially independent.
Mk is the n × (P·m) coefficient matrix, which is defined as
1 2
, 0 , , 0 , P, 0 ,
k n n n m n n n n m n n n n m n
M =⎡⎣M − M − L M − ⎤⎦ (3-5)
where Min,n denotes n × n coefficient matrix, relating zi to Yk; 0n,m-n denotes n × (m-n) zero matrix, attaching to each Min,n to eliminate the unknown parameters at m-n arbitrary points (i.e. xn+1, xn+2, ... , xm) from the observation model.
Given Z and σε2, the predicted settlement distribution at any time step k can be represented by the following multivariate normal distribution
(
k , 2) ( )
2 n/ 2 1/ 2 exp 12(
k k)
T 1(
k k)
p Y Z σε = π − Vε − ⋅ ⎡⎢⎣− Y −M Z Vε− Y −M Z ⎤⎥⎦ (3-6)
3.2.2 Prior information model
It is assumed that the prior information of the unknown parameters has the following structure
Z =Z0+δ (3-7)
where Z0 is the prior mean vector (P·m dimension) at points x1, x2, ... , xm. It can be defined as
0 0,1, 0,2 , , 0,P T
Z =⎡⎣z z L z ⎤⎦ (3-8)
where
0,i 0,i( ) ,1 0,i( ) ,2 , 0,i( )m ; 1, 2,...,
z =⎡⎣z x z x L z x ⎤⎦ i= P (3-9)
z0,i(xj) can be generally defined as z0,i(xj) = Ci1+ Ci2·x’j + Ci3·y’j + Ci4·x’j y’j + …, depending on the shape of the trend components considered to be suitable for the specific model parameters. Note that x’j and y’j denote spatial coordinates at point xj, while Ci1, Ci2, Ci3, … represent the constant coefficients of the trend for the corresponding unknown parameters zi (i = 1, 2, … , P). These parameters can be either deterministic or unknown, depending on the assumption made. For the later case, these coefficients can be estimated as one of the hyperparameters based on ABIC which will be presented in Section 3.3.
δ represents the uncertainty of the prior mean of the unknown parameters which is assumed to follow N(0, VZ) where VZ is a covariance matrix. By introducing the spatial correlation structure in the formulation of VZ, we have
2 1
2 2
2
0
0
z C
z C Z
zP C
V V V
V σ
σ
σ
⎡ ⎤
⎢ ⎥
⎢ ⎥
=⎢ ⎥
⎢ ⎥
⎢ ⎥
⎣ ⎦
O
(3-10)
where σ2z1, σ2z2, ... , σ2zP represent the prior variance of the unknown parameters z1, z2, ... , zP, respectively. These variances also can be assumed to be either deterministic or unknown and, in the same way with the prior means, they can be estimated as one of the hyperparameters based on ABIC (see Section 3.3 for detail). VC is the auto-covariance matrix which is defined as
( ) ( )
( ) ( )
⎥⎥⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−
−
−
=
m m m
m C
x x x
x
x x x
x V
ρ ρ
ρ ρ
L
M O
M
L
1
1 1
1 (3-11)
ρ(|xi - xj|) denotes the auto-correlation function where xi, xj = spatial vector coordinate. Several analytical expressions have been proposed for the auto-correlation function but, in fact, none of them can claim any fundamental basis (Vanmarcke 1977a). The exponential type auto-correlation function is chosen for the current study because it is commonly used in geotechnical applications (e.g.
Vanmarcke 1977a, Fenton and Griffiths 2002, Griffiths and Fenton 2004 etc.). The function is given as
( ) [
η]
ρ xi −xj =exp− xi −xj / (3-12)
where η = auto-correlation distance. To emphasize, this parameter is assumed to be constant at any directions in the horizontal plane. This implies that the anisotropy of soil is not considered in this case. In addition, it should be kept in mind that this type of auto-correlation function is, in fact, chosen only as an example for an application of the proposed method. In practice, several types of autocorrelation functions may be tested and the one which fits most to the observation should be used.
From the above definitions, it is clear that the spatial correlation structure is included in the form of the spatial correlation of unknown parameters, which relate to soil properties, instead of that of ground behavior. The authors believe that this is the most suitable way to introduce the spatial correlation structure to the geotechnical model due to the fact that the physical correlation of the observed ground behavior actually results from the spatial correlation of soil properties.
It should also be noted that, for the sake of simplification, there are two important assumptions about the correlation structure for formulating the above covariance matrix (VZ). Firstly, the unknown parameters, z1, z2, ... , zP, are assumed to be independent of each other. Secondly, the correlation structures of these parameters are identical, meaning that they share the same auto-correlation distance. In fact, these assumptions can be released without major change of the
formulation, if the observation data is available in the amount that the detail specification of the spatial correlation is possible.
Given η, prior means, and prior variances of the unknown parameters, the prior distribution of the unknown parameters is also assumed as a multivariate normal distribution of the following form
( ) ( )
2 (P m)2 Z 1/ 2exp 12(
0)
T Z1(
0)
p Zη = π − ⋅ V − ⎡⎢⎣− Z−Z V− Z−Z ⎤⎥⎦ (3-13)
3.2.3 Bayesian estimation
Suppose that the set of observations Yk at the discrete time step k = 1, 2, ... , K has already been obtained. By employing Bayes’ theorem, the posterior distribution of the state vector Z can be formulated as
(
, 2,) ( )
K1(
k , 2)
k
p Z Y σ ηε c p Zη p Y Z σε
=
= ⋅
∏
(3-14)where Y denotes the set of all observation data, i.e. Y = (Y1, Y2, ... , YK), and c denotes the normalizing constant. By substituting Eq. (3-6) and (3-13) into the above equation, we have
(
, 2,) ( )
2 [P m K n]2 Z 1/ 2 K/ 2p Z Y σ ηε = ⋅c π − ⋅ + ⋅ V − Vε −
(
0)
1(
0) ( )
1( )
1
exp 1 2
T K T
Z k k k k
k
Z Z V− Z Z Y M Z Vε− Y M Z
=
⎧ ⎡ ⎤⎫
⋅ ⎨⎩− ⎢⎣ − − +
∑
− − ⎥⎦⎬⎭ (3-15)The Bayesian estimator of Z, i.e. ˆZ, is the one that maximizes the above function. Therefore, it is equivalent to minimizing the following objective function
(
2,) (
0)
T Z1(
0)
K1(
k k)
T 1(
k k)
k
J Z σ ηε Z Z V− Z Z Y M Z Vε− Y M Z
=
= − − +
∑
− − (3-16)It should be noted that σε2 and η are assumed to be given in this case. The Bayesian method, however, does not provide the rational way to determine these values. In order to choose the most
appropriate values of σε2 and η based on the information in hand, Akaike’s Bayesian Information Criterion (ABIC) is introduced and presented in the next section.