Markov Chain Monte Carlo Simulation to Estimate Material Properties of a Layered Half-space

Jin Ho Lee; Hieu Van Nguyen; Se Hyeok Lee

doi:10.7734/COSEIK.2023.36.3.203

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 30 June 2023. 203-211
https://doi.org/10.7734/COSEIK.2023.36.3.203

Markov Chain Monte Carlo Simulation to Estimate Material Properties of a Layered Half-space

층상 반무한 지반의 물성치 추정을 위한 마르코프 연쇄 몬테카를로 모사 기법

Jin Ho Lee²^*

Hieu Van Nguyen¹

Se Hyeok Lee³

이 진호²^*

Nguyen Van Hieu¹

이 세혁³

¹Associate Professor, Department of Ocean Engineering, Pukyong National University, Busan, 48513, Korea

²Graduate Student, Department of Ocean Engineering, Pukyong National University, Busan, 48513, Korea

³Associate Researcher, Department of Structural Engineering Research, Korea Institute of Civil Engineering and Building Technology, Goyang, 10223, Korea

¹부경대학교 해양공학과 부교수

²부경대학교 해양공학과 박사과정

³한국건설기술연구원 구조연구본부 수석연구원

^{*Corresponding Author}

ABSTRACT

A Markov chain Monte Carlo (MCMC) simulation is proposed for probabilistic full waveform inversion (FWI) in a layered half-space. Dynamic responses on the half-space surface are estimated using the thin-layer method when a harmonic vertical force is applied. Subsequently, a posterior probability distribution function and the corresponding objective function are formulated to minimize the difference between estimations and observed data as well as that of model parameters from prior information. Based on the gradient of the objective function, a proposal distribution and an acceptance probability for MCMC samples are proposed. The proposed MCMC simulation is applied to several layered half-space examples. It is demonstrated that the proposed MCMC simulation for probabilistic FWI can estimate probabilistic material properties such as the shear-wave velocities of a layered half-space.

Keywords

markov chain monte carlo simulation

full waveform inversion

inverse problem

layered half-space

thin-layer method

층상 반무한체에서의 확률론적 완전파형역산을 위한 Markov chain Monte Carlo (MCMC) 모사 기법을 정식화한다. Thin-layer method를 사용하여 조화 수직 하중이 작용하는 층상 반무한체의 지표면에서 추정된 동적 응답과 관측 데이터와의 차이 및 모델 변수의 사전 정보와의 차이를 최소화하도록 목적함수와 모델 변수의 사후 확률밀도함수를 정의한다. 목적함수의 기울기에 기반하여 MCMC 표본을 제안하기 위한 분포함수와 이를 수락 또는 거절할지 결정하는 수락함수를 결정한다. 기본 진동모드 뿐만이 아니라 고차 진동모드가 우세한 경우를 포함하여 다양한 층상 반무한체의 전단파 속도 추정에 제안된 MCMC 모사 기법을 적용하고 그 정확성을 검증한다. 제안된 확률론적 완전파형역산을 위한 MCMC 모사 기법은 층상 반무한체의 전단파 속도와 같은 재료 특성의 확률적 특성을 추정하는 데 적합함을 확인할 수 있다.

키워드

마르코프 연쇄 몬테카를로 모사

완전파형역산

역해석

층상 반무한체

thin-layer method

MAIN

1. 서 론
2. 층상 반무한체의 동적 응답
3. 마르코프 연쇄 몬테카를로 시뮬레이션
4. 적용 예제
5. 결 론

1. 서 론

완전파형역산(full wave inversion) 기법을 이용하여 무한 매질의 물성치를 추정하는 기술은 지구물리학적 지하구조 영상화, 지반공학 부지 조사 및 비파괴 검사 등 관련 분야의 다양한 문제에 적용되어 왔다(Tarantola, 1987; Virieux and Operto, 2009). 완전파형역산 기법은 일반적인 불규칙 구조를 가진 매질에 적용 가능하지만, 파전파 현상과 이와 관련된 역해석 문제가 층상 매질에서 발생하는 경우가 많기 때문에 층상 매질의 특성을 고려한 완전파형역산 기법의 개발이 지속적으로 연구되어 왔다. 특히, thin-layer method를 사용하여 층상 매질에서의 역해석 문제에 대한 해를 구하는 방법이 개발되기도 하였다(Lee and Lee, 2021a; 2021b; Lee et al., 2022; Mashayekh et al., 2018; Nguyen et al., 2022). 이 방법은 수평 방향 또는 반경 방향으로의 푸리에 변환을 해석적으로 수행하고, 층상 구조를 이루고 있는 수직 방향으로는 유한요소의 개념을 도입하기 때문에 층상 매질에서의 파전파 문제에 대한 해를 정확하고 효율적으로 얻을 수 있다(Kausel, 1981; 1996).

역해석 문제의 대상이 되는 물리계(physical system)는 모델 공간(model space)에서 정의되는 모델 변수(model parameter)에 의해 특징 지어진다. 관측 데이터와 추정 데이터의 차이를 최소화하여 얻어지는 모델 변수는 ill-posed problem의 해를 구하여 얻어지고, 모델의 불확실성, 사전 정보의 불확실성, 관측 데이터의 잡음 등으로 인해 불확실성을 가지게 된다. 이와 같은 불확실성은 모델 공간에서 정의된 확률분포함수를 사용하여 표현 가능하므로, 이를 산정하기 위한 확률론적 역해석 기법의 개발이 필요하다. 많은 연구에서는 모델 변수의 사전 정보를 관측 데이터 및 모델의 불확실성과 결합한 베이즈 추론 방법이 확률론적 역해석을 위해 많이 사용되고 있다(Tarantola, 1987). 이때, 일반적으로 비선형성을 보이는 정해석(forward analysis) 문제를 선형화하고, 사후 확률밀도함수는 정규분포를 maximum a posterior model에 국부적으로 근사시켜 추정한다. 그러므로 역해석 문제의 비선형성 등을 고려할 때, 이와 같이 추정한 역해석 문제의 해는 가능한 해를 완전히 나타내지 않을 수도 있다. 이 문제를 해결하기 위한 방법의 하나로서 목적함수의 기울기(gradient)에 기반한 마르코프 연쇄 몬테카를로 (Markov chain Monte Carlo, MCMC) 모사 기법을 사용할 수 있다(Martin et al., 2012; Qi and Minka, 2002; Zhao and Sen, 2021). 이 기법에서는 제안함수에 따라 모델 공간을 직접 탐색하여 사후 확률분포함수를 추정하므로, 비선형 정해석 문제의 선형화와 사후 확률밀도함수의 국부 최적화로 인한 문제를 해결할 수 있다.

이 연구에서는 층상 반무한체에서의 확률론적 완전파형역산을 위한 MCMC 모사 기법을 정식화하고자 한다. 단, 반무한 영역에서의 파전파 현상을 정확하고 효율적으로 모사하기 위하여 진동수 영역 완전파형역산을 수행할 것이다(Eslaminia et al., 2022; Virieux and Operto, 2009). Thin-layer method를 사용하여 조화 수직 하중이 작용하는 층상 반무한체의 지표면에서 추정된 동적 응답과 관측 데이터와의 차이 및 모델 변수의 사전 정보와의 차이를 최소화하도록 목적함수와 모델 변수의 사후 확률밀도함수를 정의한다. 이 목적함수의 기울기에 기반하여 MCMC 표본을 제안하기 위한 분포함수와 이를 수락 또는 거절할지 결정하는 수락함수를 결정한다. 기본 진동모드 뿐만이 아니라 고차 진동모드가 우세한 경우를 포함하여 다양한 층상 반무한체의 전단파 속도 추정에 제안된 MCMC 모사 기법을 적용하여 그 정확성을 검증한다.

2. 층상 반무한체의 동적 응답

이 연구에서는 지표면에 작용하는 조화 수직 원형 분포하중에 대한 층상 반무한체의 수직 변위 응답(Fig. 1)을 구하는 것을 정해석 문제로 정의한다. 이 문제에 대한 해는 유한요소법과 같은 수치적 방법을 사용하여 구할 수 있지만, 이 연구에서는 층상 매질의 특성을 고려할 수 있는 thin-layer method를 사용하여 그 해를 구할 것이다(Kausel, 1981; 1996). 그 자세한 과정은 Kausel(1981)에 제시되어 있으므로, 여기에서는 그 결과만을 정리하여 제시하고자 한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F1.jpg

Fig. 1.

Layered half-space

표면에 조화 수직 원형 분포하중이 작용하는 층상 반무한체의 표면에서의 수직 변위 $u_{z} (r, ω)$ 는 다음과 같이 얻어진다(Kausel, 1981).

(1a)

u_{z} (r, ω) = q R \sum_{l = 1}^{2 N} (ϕ_{l, z 1})^{2} I_{l}

(1b)

I_{l} = \frac{π}{2 i k_{l}} J_{0} (k_{l} r) H_{1}^{(2)} (k_{l} R) - \frac{1}{R k_{l}^{2}} for 0 \leq r \leq R

(1c)

I_{l} = \frac{π}{2 i k_{l}} J_{1} (k_{l} R) H_{0}^{(2)} (k_{l} r) for r \geq R

여기서 r은 원통형 좌표계에서의 반경방향 좌표, 𝜔는 가진 진동수, q와 R은 각각 수직 원형 분포하중의 강도와 반경, $k_{l}$ 은 l번째 모드의 고유값, $ϕ_{l, z 1}$ 은 l번째 고유벡터 $ϕ_{l}$ 의 요소 중 표면에서의 수직 변위에 해당하는 요소, $J_{0}$ 와 $J_{1}$ 은 각각 제1종 0차 또는 1차 Bessel 함수, $H_{0}^{(2)}$ 와 $H_{1}^{(2)}$ 는 각각 제2종 0차 또는 1차 Hankel 함수, N은 thin-layer element의 수, $i = \sqrt{- 1}$ 는 허수 단위이다. 고유값 $k_{l}$ 과 고유벡터 $ϕ_{l}$ 은 층상 반무한체의 다음과 같은 고유값 문제로부터 얻어진다.

(2)

[k_{l}^{2} A + k_{l} B + G - ω^{2} M] ϕ_{l} = [k_{l}^{2} A + k_{l} B + C] ϕ_{l} = 0

여기서, $C = G - ω^{2} M$ 이고, 식 (2)의 행렬 A, B, G, M은 각각의 유한요소 행렬로부터 구성된다(Kausel, 1981). Thin-layer method에 대한 최초의 연구에서는 기존의 선형 유한요소를 사용하여 층상 매질을 나타내었지만, 최근에는 mid-point integrated finite element 또는 complex-length finite element를 사용하기도 한다(Astaneh and Guddati, 2016a; 2016b; Guddati et al., 2016). 이 수치모형을 사용하면 이산화의 영향이 제거된 층상 반무한체의 분산 방정식을 효과적으로 얻을 수 있기 때문에, 이 연구에서는 mid-point integrated finite element를 사용하여 층상 매질을 나타낼 것이다. 또한, 기저를 이루는 반무한체에서 발생하는 에너지 방사를 고려하기 위해 perfectly matched discrete layer(PMDL)을 사용할 것이다(Lee and Tassoulas, 2011). Mid- point integrated finite element와 PMDL로 구성된 N개의 유한요소를 사용하여 층상 매질을 나타내면, 2N개의 Rayleigh파 모드와 N개의 Love파 모드가 식 (2)로부터 얻어지는데, 수직 원형 분포하중에 대한 응답을 얻어야 하므로 이 연구에서는 Rayleigh파 모드만 고려할 것이다. l번째 Rayleigh 고유모드는 다음을 만족하도록 정규화한다.

(3)

\frac{1}{2} ϕ_{l}^{T} (2 k_{l} A + B) ϕ_{l} = k_{l}

식 (1)은 앞에서 정의한 정해석 문제의 해를 나타내고, 이 연구에서는 이에 기반하여 마르코프 연쇄 몬테카를로 시뮬레이션을 정식화할 것이다.

3. 마르코프 연쇄 몬테카를로 시뮬레이션

확률론적 역해석을 위해 MCMC 모사 기법을 정식화한다. 이를 위해서는 표본을 제안하기 위한 분포함수(proposal distribution)와 이를 수락 또는 거절할지 결정하는 함수(acceptance function)를 결정해야 한다. 이 연구에서는 목적함수의 기울기에 기반하여 제안함수와 수락함수를 결정하고자 한다.

관측 데이터 $d_{obs}$ 가 주어졌을 때, 모델 변수 m의 사후 확률밀도함수 $p_{post} ({m d}_{obs})$ 는 베이즈 정리에 의해 다음과 같이 주어진다(Tarantola, 1987; Virieux and Operto, 2009).

(4a)

p_{post} ({m | d}_{obs}) \propto p_{like} (d_{obs} | m) p_{prior} (m) \propto \exp [- \frac{1}{2} (d_{cal} (m) - d_{obs})^{T} C_{obs}^{- 1} (d_{cal} (m) - d_{obs})] ∙ \exp [- \frac{1}{2} (m - m_{prior})^{T} C_{prior}^{- 1} (m - m_{prior})] \propto \exp [- E (m)]

(4b)

E (m) = - \frac{1}{2} (d_{cal} (m) - d_{obs})^{T} C_{obs}^{- 1} (d_{cal} (m) - d_{obs}) - \frac{1}{2} (m - m_{prior})^{T} C_{prior}^{- 1} (m - m_{prior})

여기서 $p_{like} (d_{obs} | m)$ 는 정규분포함수로 근사한 가능도(likelihood) 함수, $p_{prior} (m)$ 는 정규분포함수로 근사한 매개변수에 대한 사전 확률밀도함수, $d_{cal} (m)$ 는 정해석 문제로부터 얻어지는 예측 데이터, $C_{obs}$ 는 관측 데이터의 공분산 행렬, $m_{prior}$ 와 $C_{prior}$ 는 각각 매개변수의 사전 평균과 공분산 행렬, E(m)은 목적함수이다. $p_{post} (m | d_{obs})$ 를 최대로 하거나 E(m)을 최소로 하는 최적해는 Newton 방법을 사용한 반복 계산에 의해 얻을 수 있다. 최적해의 k번째 근사 $m_{k}$ 에 대해 식 (4a)를 다음과 같이 표현할 수 있다(Tarantola, 1987; Virieux and Operto, 2009).

(5a)

p_{post} (m | d_{obs}) \propto \sqrt{\det H_{k}} \exp (- \frac{1}{2} {\{m - (m_{k} - H_{k}^{- 1} g_{k})\}}^{T} H_{k} \{m - (m_{k} - H_{k}^{- 1} g_{k})\})

(5b)

g_{k} = J_{k}^{T} C_{o b s}^{- 1} \{d_{c a l} (m_{k}) - d_{obs}\} + C_{prior}^{- 1} (m_{k} - m_{prior})

(5c)

H_{k} = J_{k}^{T} C_{obs}^{- 1} J_{k} + R_{k} + C_{prior}^{- 1} \approx J_{k}^{T} C_{obs}^{- 1} J_{k} + C_{prior}^{- 1}

(5d)

J_{k} = \frac{\partial d_{cal}}{\partial m} |_{{m = m}_{k}}

(5e)

R_{k} = \frac{\partial J_{k}^{T}}{\partial m^{T}} C_{obs}^{- 1} \{d_{cal} (m_{k}) - d_{obs}\}

여기서 $g_{k}$ 와 $H_{k}$ 는 각각 목적함수 E(m)의 기울기와 Hessian matrix, $J_{k}$ 는 Jacobian matrix, sensitivity 또는 Fréchet derivative matrix, $R_{k}$ 는 Hessian matrix의 2차항으로 최적해 근방에서는 무시할 수 있다. 식 (5c)에서는 고차 미분항을 무시하는 Gauss- Newton 방법을 적용하여 $H_{k}$ 를 근사하였다. 식 (1)에 주어진 정해석 문제의 해는 $g_{k}$ 와 $H_{k}$ 를 산정하기 위하여 $d_{cal} (m_{k})$ 와 $J_{k}$ 를 계산할 때 활용된다(Lee and Lee, 2021a).

식 (5)의 $p_{post} ({m | d}_{obs})$ 에 기반하여 제안함수 $g ({m | m}_{k})$ 를 다음과 같이 정의한다(Qi and Minka, 2002).

(6)

g (m | m_{k}) \propto \sqrt{\frac{\det H_{k}}{β^{2}}} \exp (- \frac{1}{2} {\{m - (m_{k} - α H_{k}^{- 1} g_{k})\}}^{T} \frac{H_{k}}{β^{2}} \{m - (m_{k} - α H_{k}^{- 1} g_{k})\})

여기서 𝛼는 preconditioned gradient-descent update와 random walk의 상대적 중요도를 조절하는 변수이고, 𝛽는 random walk의 크기를 결정하는 변수이다(Zhao and Sen, 2021). 식 (6)의 분포함수를 사용하여 MCMC 표본 $m'$ 을 다변수 정규분포 $N (m_{k} - α H_{k}^{- 1} g_{k}, β^{2} H_{k}^{- 1})$ 로부터 추출한다.

추출된 표본 $m'$ 을 수락 또는 거절할지 결정하는 수락함수 $γ (m', m_{k})$ 는 식 (7)과 같이 정의한다(Qi and Minka, 2002).

(7)

γ (m', m_{k}) = \min (1, \frac{p_{post} (m' | d_{obs}) g (m_{k} | m')}{p_{post} (m_{k} | d_{obs}) g (m' | m_{k})})

식 (6)과 (7)의 제안함수와 수락함수를 사용하여 목적함수의 기울기에 기반한 MCMC 모사 기법을 Table 1과 같이 정식화할 수 있다.

Table 1.

Algorithm of a gradient-based Markov chain Monte Carlo simulation

Given: Prior mean and covariance

m_{prior}

and

C_{prior}

for parameters m.
Given: Covariance

C_{obs}

representing the uncertainties in observations

d_{obs}

.
Choose

m_{0}

, 𝛼, and 𝛽.
For

k = 0, 1, \dots, N - 1

Jacobian matrix

J_{k} = \frac{\partial d_{cal}}{\partial m} |_{m = m_{k}}

.
Gradient of a objective function

g_{k} = J_{k}^{T} C_{o b s}^{- 1} \{d_{c a l} (m_{k}) - d_{obs}\} + C_{prior}^{- 1} (m_{k} - m_{prior})

.
Hessian matrix

H_{k} = J_{k}^{T} C_{obs}^{- 1} J_{k} + R_{k} + C_{prior}^{- 1} \approx J_{k}^{T} C_{obs}^{- 1} J_{k} + C_{prior}^{- 1}

where

R_{k} = \frac{\partial J_{k}^{T}}{\partial m^{T}} C_{obs}^{- 1} \{d_{cal} (m_{k}) - d_{obs}\}

.
Posterior distribution

p_{post} (m | d_{obs}) \propto \sqrt{\det H_{k}} \exp (- \frac{1}{2} {\{m - (m_{k} - H_{k}^{- 1} g_{k})\}}^{T} H_{k} \{m - (m_{k} - H_{k}^{- 1} g_{k})\})

.
Proposal distribution

g (m | m_{k}) \propto \sqrt{\frac{\det H_{k}}{β^{2}}} \exp (- \frac{1}{2} {\{m - (m_{k} - α H_{k}^{- 1} g_{k})\}}^{T} \frac{H_{k}}{β^{2}} \{m - (m_{k} - α H_{k}^{- 1} g_{k})\})

.
Draw sample

m'

from a normal distribution

N (m_{k} - α H_{k}^{- 1} g_{k}, β^{2} H_{k}^{- 1})

γ (m', m_{k}) = \min (1, \frac{p_{post} (m' | d_{obs}) g (m_{k} | m')}{p_{post} (m_{k} | d_{obs}) g (m' | m_{k})})

Draw a random number u from a uniform distribution on [0,1].
If

u \leq γ (m', m_{k})

then accept:

m_{k + 1} = m'

else reject:

m_{k + 1} = m_{k}

end.
end.

Table 2.

Structures of layered half-spaces

Layer	Shear-wave velocity (m/s)			Thickness (m)	Poisson’s ratio	Density (kg/m³)	Damping ratio
Layer	Model 1	Model 2	Model 3	Thickness (m)	Poisson’s ratio	Density (kg/m³)	Damping ratio
1	200	300	300	6	0.35	1800	0.02
2	300	200	400	4
3	400	400	200	8
Half-space	500	500	500	∞

4. 적용 예제

3장의 MCMC 모사 기법을 사용하여 예제 층상 반무한체(Table 2)의 전단파 속도 $C_{S, m} (m = 1, 2, \dots, N_{p})$ 를 추정한다. 모델 1에서는 최저차의 기본 모드가 반무한체의 동적 거동에 주로 영향을 주지만, 모델 2와 3에서는 전단파 속도가 큰 두 개의 층 사이에 전단파 속도가 작은 층이 위치하고 있어 다양한 고차 모드가 시스템의 동적 거동에 영향을 미친다. 전단파 속도를 20m 깊이까지 결정하기 위하여 두께가 0.5m인 40개의 mid- point integrated finite element와 4개의 PMDL로 반무한체의 수치 모형을 구성한다. 즉, N=44이다. 기저를 이루고 있는 PMDL은 하나의 모델 변수를 공유하기 때문에 추정해야 하는 모델 변수 m의 수는 41, 즉 $N_{p}$ = 41이다. 반무한체의 영향을 고려하기 위해 PMDL을 사용하기 때문에 반무한체의 전단파 속도가 균일(homogeneous)한 경우에만 적용 가능하지만, 일반적으로 대부분의 지반공학적 문제에서 반무한 지반의 물성치는 일정 깊이부터는 균일하다고 가정하기 때문에 이 연구에서 제안하는 역해석 기법의 적용성이 제한되지는 않는다.

이 에제에서 관측 데이터 $d_{obs}$ 및 예측 데이터 $d_{cal}$ 는 반무한체 지표면에서의 수직 변위로 정의되는데, 각각 식 (1)의 정해석 기법을 사용하여 합성되거나 계산한다. 관측 및 예측 데이터의 계산을 위해 반경 R 및 강도 q가 각각 0.1m이고 $1 / π R^{2}$ = 31.83N/m²인 수직 원형 분포하중에 의해 발생하는 반무한체 지표면에서의 수직 변위를 r = 5,6,…,19m의 위치에서 $f = ω / 2 π$ = 5,10,…50Hz의 진동수에 대하여 계산한다. 단, inverse crime을 방지하기 위하여 13dB의 신호 대 잡음비(signal-to- noise ratio, SNR)(잡음의 진폭비 5%)를 가지는 백색 가우스 잡음을 합성 관측 데이터에 추가한다. 여기서, inverse crime이란 합성 데이터를 생성한 수치 모형과 동일한 모형을 사용하여 역해석을 수행했을 때 역해석의 알고리즘이 옳지 않더라도 목표치와 잘 일치하는 추정치를 산출하는 현상이다.

역해석 문제로부터 얻어진 해의 매끈함(smoothness)를 보장하기 위해 예측 데이터와 관측 데이터에 각각 $C_{S, m + 1} - C_{S, m} (m = 1, 2, \dots, N_{p} - 1)$ 과 ( $N_{p} - 1$ )개의 0을 추가한다. 이는 연속하는 두 층의 전단파 속도에 대한 차분이 0를 중심으로 정규 분포하도록 하여 전단파 속도 분포의 매끈함을 보장한다. 단, 지표면에서의 수직 변위와 전단파 속도는 크기의 차이가 있으므로, 수직 변위와 전단파 속도를 각각 지표면에서 발생하는 수직 변위의 최고값과 100 m/s로 정규화하여 유사한 크기를 갖도록 한다. 모델 1~3에 대하여 잡음이 없을 때 지표면에서 발생하는 수직 변위의 최고값 $\max | u_{z}^{obs} |$ 은 각각 4.3097 × 10^-10m, 2.3644 × 10^-10m, 1.9954 × 10^-10m이다.

3장의 MCMC 모사 기법을 사용하여 전단파 속도의 수직 분포를 추정한다. 전단파 속도의 사전 평균과 표준표차는 각각 350m/s와 70m/s이고, 하한 및 상한은 각각 100m/s 및 1000m/s 이다. 관측 수직 변위의 표준편차는 지표면에서 관측된 수직 변위 최고값 $\max | u_{z}^{obs} |$ 의 1%로 가정한다. 전단파 속도 차이 $C_{S, m + 1} - C_{S, m}$ 의 평균과 표준편차는 각각 0와 $0.01 / \sqrt{2 C_{R}}$ 로 가정하는데 모델 1~3에 대하여 $C_{R}$ 은 각각 0.1, 0.01, 0.01이다. 𝛼 = 0.1, 𝛽 = 1의 값을 사용하여 총 10000개의 MCMC 표본을 생성하였는데, burn-in 현상을 고려하여 초기 1000개의 표본을 제외한 9000개의 표본을 결과로서 생성된 표본으로 간주하였다. Burn-in 구간은 마르코프 연쇄가 목표 확률밀도함수를 대표하지 못하는 초기 표본에서 이를 잘 대표할 수 있는 표본으로 발전하는 예비 구간으로 일반적으로 이 구간의 표본을 제외하고 최종 표본을 제시한다.

Fig. 2, 3, 4는 모델 1~3에 대하여 생성된 MCMC 표본의 히스토그램을 보여주고 있고, Fig. 5는 평균 전단파 속도의 수직 분포를 보여준다. Fig. 6, 7, 8은 Fig. 5의 평균 전단파 속도 수직 분포로부터 계산된 지표면 변위와 합성된 관측 데이터를 비교하여 보여주고 있는데, 잡음이 포함되지 않은 관측 데이터도 함께 보여주고 있다. Fig. 5, 6, 7, 8의 비교로부터 MCMC 모사기법을 사용하여 층상 반무한체의 전단파 속도를 만족스럽게 추정할 수 있음을 확인할 수 있다. 역해석 결과의 정확성을 조사하기 위해 전단파 속도에 대한 상대 오차를 다음과 같이 계산한다.

(8)

E = \sqrt{\frac{\sum_{m = 1}^{N_{p}} (C_{S, m}^{target} - C_{S, m})^{2}}{\sum_{m = 1}^{N_{p}} (C_{S, m}^{target})^{2}}}

여기서 $C_{s, m}^{target}$ 은 m번째 요소의 목표 전단파 속도이다. 식 (8)로부터 계산된 모델 1, 2, 3의 상대오차가 Table 3에 기존 연구의 상대오차와 함께 제시되어 있다. Lee 등(2022)을 제외하고, 이 연구의 결과는 기존 연구의 결과보다 우수한 것을 확인할 수 있다. Lee 등(2022)에서는 불연속적인 전단파 속도 수직 분포를 고려할 수 있기 때문에, 연속적인 분포를 고려하는 다른 연구보다 우수한 결과를 얻을 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F2.jpg

Fig. 2.

Shear-wave velocities of Model 1

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F3.jpg

Fig. 3.

Shear-wave velocities of Model 2

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F4.jpg

Fig. 4.

Shear-wave velocities of Model 3

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F5.jpg

Fig. 5.

Estimated profiles of shear-wave velocities in layered half-spaces (black: target, red: estimated)

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F6.jpg

Fig. 6.

Vertical displacements on the surface of Model 1

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F7.jpg

Fig. 7.

Vertical displacements on the surface of Model 2

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-03/N0040360307/images/Figure_jcoseik_36_03_07_F8.jpg

Fig. 8.

Vertical displacements on the surface of Model 3

Table 3.

Relative errors (%)

	This study	Lee and Lee (2021a)	Lee and Lee (2021b)	Lee et al. (2022)
Model 1	7.27	9.77	8.81	5.66
Model 2	6.34	9.31	10.10	8.28
Model 3	15.70	19.60	21.22	11.41

이 예제를 통해서 MCMC 표본을 활용한 확률론적 FWI 기법을 사용하여 다양한 층상 반무한체의 전단파 속도 등과 같은 물성치의 확률분포함수를 성공적으로 추정할 수 있는 것을 확인할 수 있다.

5. 결 론

이 연구에서는 목적함수의 기울기에 기반하여 층상 반무한체에서의 확률론적 완전파형역산을 위한 MCMC 모사 기법을 정식화하였다. 이를 위해 MCMC 표본을 제안하기 위한 분포함수와 이를 수락 또는 거절할지 결정하는 수락함수를 목적함수의 기울기에 기반하여 결정하였다. 조화 수직 하중이 작용하는 층상 반무한체의 지표면에서 추정된 동적 응답과 관측 데이터와의 차이 및 모델 변수의 사전 정보와의 차이를 최소화하도록 목적함수와 모델 변수의 사후 확률밀도함수가 정의하였다. 이때, 층상 반무한체의 동적응답은 mid-point integrated finite element와 PMDL을 사용하여 구성된 thin-layer model의 정해석 문제의 해로부터 얻어졌다. 모델 변수의 사후 확률밀도함수를 기반으로 제안함수와 수락함수를 결정하고, 이로부터 MCMC 표본을 생성하였다.

제안된 기법을 다양한 층상 반무한체의 전단파 속도 추정에 적용하였는데, 기본 진동모드 뿐만이 아니라 고차 진동모드가 우세한 경우도 포함하였다. 적용 예제를 통해 제안된 확률론적 완전파형역산을 위한 MCMC 모사 기법이 층상 반무한체의 전단파 속도와 같은 재료 특성의 확률적 특성을 추정하는 데 적합함을 확인할 수 있었다. 향후 연구에서는 불규칙한 재료 특성과 기하 형상을 가진 매질에 대한 확장 및 적용이 다루어져야 할 것이다.

Acknowledgements

이 논문(또는 저서)은 부경대학교 자율창의학술연구비(2021년)에 의하여 연구되었음.

References

Astaneh, A.V., Guddati, M.N. (2016a) Improved Inversion Algorithms for Near-Surface Characterization, Geophys. J. Int., 206(2), pp.1410~1423. 10.1093/gji/ggw192

Astaneh, A.V., Guddati, M.N. (2016b) Efficient Computation of Dispersion Curves for Multilayered Waveguides and Half- spaces, Comput. Methods Appl. Mech. & Eng., 300, pp. 27~46. 10.1016/j.cma.2015.11.019

Eslaminia, M, Elmeliegy, A.M., Guddati, M.N. (2022) Full Waveform Inversion through Double-Sweeping Solver, J. Comput. Phys., 453, p.110914. 10.1016/j.jcp.2021.110914

Guddati, M.N., Druskin, V., Astaneh, A.V. (2016) Exponential Convergence through Linear Finite Element Discretization of Stratified Subdomains, J. Comput. Phys., 322, pp.492~447. 10.1016/j.jcp.2016.06.045

Kausel, E. (1981) An Explicit Solution for the Green's Functions for Dynamic Loads in Layered Media, MIT Research Report R81-13, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA.

Kausel, E. (1996) Discussion: Dynamic Response of a Multi- Layered Poroelastic Medium by Rajapakse RKND, Senjuntichai T. Earthq. Eng. & Struct. Dyn., 25(10), pp.1165~1167. 10.1002/(SICI)1096-9845(199610)25:10<1165::AID-EQE600>3.0.CO;2-I

Lee, J.H., Lee, S.-H. (2021a) Full Waveform Inversion to Estimate the Material Properties of a Layered Half-Space, Soil Dyn. & Earthq. Eng., 151, p.106956. 10.1016/j.soildyn.2021.106956

Lee, J.H., Lee, S.-H. (2021b) Estimation of Shear-Wave Velocities of a Layered Half-space Using a Full Waveform Inversion with a Genetic Algorithm, Computational Methods in Applied Mechanics and Engineering, 34(4), pp.221~230. 10.7734/COSEIK.2021.34.4.221

Lee, J.H., Tassoulas, J.L. (2011) Consistent Transmitting Boundary with Continued-Fraction Absorbing Boundary Conditions for Analysis of Soil-Structure Interaction in a Layered Half- Space, Comput. Methods Appl. Mech. & Eng., 200, pp. 1509~1525. 10.1016/j.cma.2011.01.004

Lee, S.-H., Han, S.-W., Lee, J.H. (2022) Probabilistic Full- Waveform Inversion Considering the Discontinuities of the Material Properties of a Layered Half-Space, Int. J. Geomech., 22(11), p.04022207. 10.1061/(ASCE)GM.1943-5622.0002583

Martin, J., Wilcox, L.C., Burstedd, C., Ghattas, O. (2012) A Stochastic Newton MCMC Method for Large-scale Statistical Inverse Problems with Application to Seismic Inversion, SIAM J. Sci. Comput., 34(3), pp.A1460~A1487. 10.1137/110845598

Mashayekh, H., Kallivokas, L.F., Tassoulas, J.L. (2018) Parameter Estimation in Layered Media Using Dispersion-Constrained Inversion, J. Eng. Mech., 144(11), p.04018099. 10.1061/(ASCE)EM.1943-7889.0001527

Nguyen, V.H., Lee, S.-H., Lee, J.H. (2022) Application of the Unscented Kalman Filter to Estimate the Material Properties of a Layered Half-Space, Int. J. Appl. Mech., 15(2), p.2350005. 10.1142/S1758825123500059

Qi, Y., Minka, T.P. (2002) Hessian-based Markov Chain Monte- Carlo Algorithms, Proceedings of the First Cape Cod Workshop on Monte Carlo Methods, Massachusetts, USA: Cape Cod

Tarantola, A. (1987) Inverse Problem Theory: Method for Data Fitting and Model Parameter Estimation, Society for Industrial and Applied Mathematics, Philadelphia, USA.

Virieux, J., Operto, S. (2009) An Overview of Full-Waveform Inversion in Exploration Geophysics, Geophys., 74(6), pp. WCC127~WCC152. 10.1190/1.3238367

Zhao, Z., Sen, M.K. (2021) A Gradient-Based Markov Chain Monte Carlo Method for Full-Waveform Inversion and Uncertainty Analysis, Geophys., 86(1), pp.R15~R30. 10.1190/geo2019-0585.1

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Markov Chain Monte Carlo Simulation to Estimate Material Properties of a Layered Half-space

ABSTRACT

MAIN

Fig. 1.

Layered half-space

(1a)

(1b)

(1c)

(2)

(3)

(4a)

(4b)

(5a)

(5b)

(5c)

(5d)

(5e)

(6)

(7)

Table 1.

Algorithm of a gradient-based Markov chain Monte Carlo simulation

Table 2.

Structures of layered half-spaces

(8)

Fig. 2.

Shear-wave velocities of Model 1

Fig. 3.

Shear-wave velocities of Model 2

Fig. 4.

Shear-wave velocities of Model 3

Fig. 5.

Estimated profiles of shear-wave velocities in layered half-spaces (black: target, red: estimated)

Fig. 6.

Vertical displacements on the surface of Model 1

Fig. 7.

Vertical displacements on the surface of Model 2

Fig. 8.

Vertical displacements on the surface of Model 3

Table 3.

Relative errors (%)

Acknowledgements

References