1. 서 론

2020년 1월에 공표된 지속가능한 기반시설 관리 기본법(이하, 기반시설관리법)은 기반시설의 체계적인 유지관리와 성능개선을 통하여 국민들이 안전하고 편리하게 기반시설을 활용할 수 있도록 하고 나아가 국가 경제 발전에 기여함을 목적으로 하고 있다.

기반시설관리법 제3조 제1항에 의하면 관리주체(국가 및 지방자치단체, 공공기관 등)는 기반시설의 선제적인 관리를 통해 노후화에 따른 생애주기비용을 최소화하여야 함을 명시하고 있다. 특히, 관리주체는 기반시설의 유지관리와 성능개선에 필요한 예산을 최소유지관리기준 및 성능개선기준에 입각하여 추정하고, 기반시설 관리계획(이하 관리계획)에 반영해야 한다.

비용추정은 이벤트가 발생하지 않은 불확실한 미래를 다루기 때문에 신뢰성 확보가 중요하고, 이를 위해서는 신뢰성 있는 파손예측 모형구축이 중요하다(^{Do, 2011}; ^{Han and Do, 2016}). 국가에서 관리하는 기반시설물들은 유지관리시스템을 구축하여 파손상태의 조사 및 분석을 통해 체계적인 관리가 이루어지고 있다. 하지만, 지자체 소관의 기반시설물은 예산, 인력, 파손모델의 필요성 등에 대한 인식 부족으로 인해 기반시설의 정확한 장래 예산 추정이 어려운 실정이다.

대표적인 국내 기반시설인 도로포장은 일반국도 포장관리시스템(Pavement Management System, PMS)의 도입을 통해 다른 기반시설 대비 체계적인 관리가 이루어져 오고 있다. 하지만 하드웨어인 포장설계 수준과 재료, 이용환경을 대표하는 교통하중을 설명변수로 적용하는 것이 일반적이다(^{Han et al., 2016}; ^{Kim and Kim, 2006}; ^{Prozzi and Madanat, 2004}). 그러나 포장 분야 관계자라면 모두가 영향 인자로 인식하고 있는 온도, 강우, 적설량 등 기후요소에 대한 접근은 쉽게 찾아보기 어렵다. 그 이유는 일반적인 포장관리시스템에서 수집, 관리하는 정보가 아니라는데 있다.

이러한 문제점을 해결하기 위해 ^{Choi et al.(2019)}은 기상관측데이터를 활용하여 지역별로 일반국도 도로포장의 기대수명을 추정한 바 있다. 포장관리시스템에서 취득 가능한 데이터(상태지표, 교통량, 축하중, 포장강도)와 기상청에서 관리하는 종관기상관측장비(Automated Synoptic Observing System, ASOS)의 데이터(기온, 강수량 등)를 활용하여 도로포장의 기대수명을 추정하였다. 하지만, 분석에 활용된 ASOS장비의 기상데이터는 도로포장구간을 기준으로 원거리 위치, 적은 장비 수 등에 의해 영향을 미치는 범위가 상대적으로 작을 수 있다는 문제점을 제시한 바 있다.

한편, ^{Choi et al.(2019)}는 기대수명의 추정을 위해 베이지안 마르코프 혼합 해저드(Bayesian Markov Mixture Hazard) 모델(이하, BMMH 모델)을 활용하였다. BMMH모델은 다양한 환경요인(설명변수)을 고려하여 지역별 기대수명을 추정함에 매우 큰 장점이 있다(^{Kaito et al., 2012}). BMMH모델에는 Markov Chain Monte-Carlo (MCMC)기법과 이질성 요인(Heterogeneity Factor)이 적용되어, 이전 모델에서 발생했던 파라미터의 불완전한 수렴, 데이터의 부족, 샘플의 과분산 문제 등을 해결하였기 때문이다(^{Han et al., 2013}; ^{Han et al., 2014}; ^{Kobayashi et al., 2012}). 이 모델의 중요한 특성은 일반적인 과정(Maximum Likelihood Estimation, MLE) 대신 MCMC기법을 사용한다는 점이다. MCMC기법은 완화된 조건에서 불충분한 샘플 데이터, 과도한 샘플 분산, 불안정하거나 불완전한 수렴 등의 고질적인 문제를 해결함으로써 실질적인 기대수명 추정에 있어 상당한 이점을 갖는다.

따라서, 본 연구에서는 일반국도 모니터링구간과 상대적으로 인접한 방재기상관측장비(Automatic Weather System, AWS)의 기상데이터를 설명변수로 적용하여 아스팔트 도로포장의 균열률에 대한 기대수명을 지역별로 추정하고자 하였다. 여기서 지역( $k$)이란 기후적인 특성에 의한 구분이 아닌, 일반국도 지방국토관리청에 의해 관리되어지는 5개의 지역(서울, 대전, 익산, 원주, 부산)을 의미한다. 분석에 활용한 데이터는 다음과 같다. 먼저, 종속변수는 일반국도 포장관리시스템에 의해 조사(2007~2017년)된 균열률 지표를 대상으로 하였다. 독립변수에는 교통량, 축하중, 포장강도, 기온, 강수량, 제설제 사용량 데이터가 고려되었다. 여기서, 기온, 강수량 데이터를 방재기상관측장비(AWS)에서 취득하여 적용한 것이 선행연구(^{Choi et al., 2019})와의 차별점이라 할 수 있다.

2. BMMH모델

본 연구에서는 도로포장의 파손과정 분석 및 기대수명 추정을 위해 BMMH모델을 활용하였다. BMMH모델을 설명하기 위해서는 1)마르코프 과정(Markov process), 2)다단계지수하자드(Multi State Exponential Hazard) 모델, 3)로컬 혼합 메커니즘(Local Mixing Mechanism) 모델 및 MCMC 시뮬레이션에 대한 포괄적인 이해와 4)Geweke 검정에 대한 이해가 필요하다. 본 장에서는 지면 관계상 모델 전체에 대한 상세한 설명이 어렵기 때문에, 모델에 대한 자세한 이론 및 활용에 대한 정보는 ^{Tsuda et al.(2006)}, ^{Obama et al.(2008)}, ^Geweke(1992), ^{Kobayashi et al.(2012)}를 참조하길 바란다.

2.1 다단계지수하자드모델과 마르코프 연쇄

BMMH모델을 이해하기 위해서는 먼저, 마르코프 연쇄에 대한 이해가 필요하다. 분석에 활용할 도로포장상태지표의 경우 연속적인 숫자로 표현되지만, 마르코프 연쇄에 적용하기 위해서는 $i(i=1, ..., J)$와 같은 이산형태로 표현되어야 한다. 기본적인 마르코프 연쇄모델은 시간간격 ${\tau }_A$와 ${\tau }_B$ 사이의 전이 확률(transition probabilities)인 ${\pi }_{ij}$을 산정하는 과정을 의미한다. 간격은 $z(z=\tau B-\tau A)$로 표시하며, 모델의 기본적인 정의에 따라 마르코프 전이 확률(Markov Transition Probability, MTP) $\Pi$와 ${\pi }_{ij}$는 Eqs. (1) and (2)와 같이 표현된다.

따라서, 확률의 공리를 만족시키기 위해 ${\pi }_{ij}\ge 0$ 와 $\sum _{j=1}^J{\pi }_{ij}=1$ 이 기본적인 전제조건이 된다. 또한, 도로포장의 유지보수는 파손과정을 설명하는데 불필요하기 때문에 ${\pi }_{ij}\ge 0 (i>j)$와 ${\pi }_{JJ}=1$이 추가 전제조건이 된다.

한편, Tsuda et al.(2006)은 기존의 마르코프 모델이 갖는 문제점을 해결하기 위해, 다단계지수하자드모델을 사용하여 MTP 행렬에 대한 추정 방법론을 제안하였다. 기존 모델대비, 개선이 이루어진 사항은 크게 세 가지로 나누어진다. 첫째, 전이 확률을 단위 간격으로 정규화하여 시간 조정 조건 $\prod \left(nZ\right)=[\prod (Z){]}^n$을 만족하게 하였다. 간단히 말해서, 샘플 데이터 세트들 사이에서 z의 시간 동기화를 필요로 하지 않는다. 또한, 샘플드롭문제가 있으며, 이는 최악의 조건에서 샘플 크기가 작은 경우 수렴이 되지 않는 문제를 개선시켰다. 마지막으로, 전이과정에 영향을 미치는 설명변수의 적용을 쉽게 할 수 있게 되었다. 이 모델이 BMMH 모델의 기본 구조가 된다.

2.2 마르코프 혼합 하자드 모델의 설명

각 지역 $k(k=1, ..., K)$ 또는 개별 구간 $s_k(s=1, ...S_k)$은 설명변수의 영향에 의해 서로 다른 파손 특성이 나타나게 된다. 여기서, 각 지역( $k$)이란 지방국토관리청에 의해 관리되어지는 5개의 지역(서울, 대전, 익산, 원주, 부산)을, 개별 구간( $s_k$)이란 지방국토관리청( $k$)에 위치하고 있는 일반국도 도로포장 모니터링 구간을 의미한다.

기본(bench mark)과정에서의 파손속도는 ${\epsilon }^k=1$이 기준이다. 따라서, 지역 $k$ 의 파손속도는 기본과정에서의 상대적인 비율로서 나타내어지며, 항상 양 $(0<{\epsilon }^k<\infty )$의 값이 나타나게 된다. 즉, ${\epsilon }^k$이 1 보다 크면 파손속도가 상대적으로 빠르며, 반대일 경우의 파손속도는 상대적으로 느리다고 해석이 가능하다.

이질성 인자 ${\epsilon }^k$는 각 지역의 다양한 파손특성을 반영하여 확률 변수 또는 함수의 형태로 표현할 수 있다. 여기서, 감마 분포 ${\epsilon }^{k }\varsigma (a,\gamma )$를 기반으로 ${\epsilon }^k$가 도출된다고 가정하였으며 함수 $f({\epsilon }^k;a,\gamma )$는 Eq. (4)와 같다.

여기서, $\Gamma (.)$는 감마 함수, $a$ 와 $\gamma$는 감마 분포의 파라미터이다. 또한, $\alpha \gamma$=1 과 $\alpha {\gamma }^2=1/\Phi$인 경우의 감마 분포의 확률 밀도 함수(Probability Density Function, PDF)는 Eq. (5)로 정리된다.

여기서, 이질성 요인의 값을 ${\overline{\epsilon }}^k$ 로 설정하였다. Eq. (3)을 이용할 경우 특정구간 $s_k$이 상태등급 $i$에서의 기대수명이 $y_i^{s_k}$이상이 될 확률 $R_i\left(y_i^{s_k}\right)$는 Eq. (6)과 같이 변환가능하다.

일반적으로 $\left[\overline{}\right]$는 특정 변수의 평균값을 의미하지만, 본 논문에서의 $\left[\overline{}\right]$는 측정 가능한 값을 의미한다. Eq. (6)은 비슷한 형태의 다단계지수해저드모델에 적용되며, 데이터 조사 간격 $z$에 대해 ${\pi }_{ii}, {\pi }_{i,i+1}, {\pi }_{ij}$ 그리고 ${\pi }_{iJ}$로 분류 가능한 파손경로를 나타낸다.

$\mathrm{where},{\Psi }_{ij}^s\left({\tilde{\lambda }}_i^{s_k}\right)=\prod _{m=i}^{j-1}\frac{{\widetilde{{\lambda }_m}}^{s_k}}{{\tilde{\lambda }}_m^{s_k}-{\tilde{\lambda }}_s^{s_k}} (i=1, ..., J-1; j=i+1, ..., J; k=1, ..., K)$

전제조건이 $\sum _{j=1}^J{\pi }_{ij}=1$인 경우, ${\pi }_{iJ}$는 Eq. (9)로 쉽게 추정할 수 있다.

감마 분포에서의 이질성 요인 ${\epsilon }^k$을 산출하기 위해 Eq. (5)를 사용하였다. 여기서, Eq. (10)을 활용하면 각각의 파손경로 ${\tilde{\pi }}_{ii}\left(z^{s_k}\right)$ 와 ${\tilde{\pi }}_{ij}\left(z^{s_k}\right)$의 확률을 산출할 수 있다.

Eq. (7), Eq. (8), Eq. (9), Eq. (10)의 수식을 활용할 경우, 이질성 요인이 적용되어 지역별 파손특성을 손쉽게 도출할 수 있지만, 환경변수의 특성을 고려하지는 못한다. 해당 문제는 ${\tilde{\lambda }}_i^{ s_k}$ 추정방법을 미지의 파라미터 ${\beta }_i=({\beta }_{i,1}, ..., {\beta }_{i,M})$와 공변량 함수 ${\overline{x}}^{S_k}$를 활용하여 변경하여야 한다. 여기서 $m(m=1,...,M)$은 설명변수의 개수를 의미하며, Eq. (11)과 같이 표현이 가능하다.

Eq. (11)에서 “ $\text{'}$”은 전치 행렬을 의미한다. 산정된 위험률 함수 ${\tilde{\lambda }}_i^{s_k}(i=1, ..., J-1)$를 적용하면, 각 상태등급 $i$에서의 기대수명은 생존 함수 $R_i\left(y_i^{s_{k)}}\right)$의 평균으로 정의할 수 있다(^{Lancaster, 1990}).

즉, 상태등급 $i$로부터 $J$까지의 기대수명은 상태등급별 기대수명의 누적 합으로 정의할 수 있다.

2.3 계층적 베이지안 추정

2.3.1 베이지안 정리 및 모형 모수

베이즈 정리(^{Bayes and Price, 1763})는 약 250년 전에 제안되었다. 대부분 베이즈 정리는 조건부 확률을 계산하는 데 사용된다고 알려져 있다. 하지만, 최근에는 MCMC기법의 개발에 따른 기준화 정수(normalizing constant) 산정과 관련된 문제해결 및 사전분포와 사후분포의 관계를 해결하는 등 광범위하게 적용이 이루어지고 있다. 다수의 문헌들은 1)불충분한 데이터, 2)최대값의 수렴, 3)계산의 복잡성 문제를 해결함에 있어 베이지안 추정기법이 큰 장점을 갖고 있다고 강조하고 있다(^{Train, 2009}). BMMH모델은 단계적인 과정이 아닌, 한 번에 해당 파라미터를 취득할 수 있다는 또 다른 장점이 있으며, 이는 MLE를 적용하는 기존의 혼합하자드모델과 비교할 경우 큰 장점이라고 할 수 있다(^{Obama et al., 2008}).

사후분포 파라미터 $\theta$는 사전분포 파라미터 및 관측 데이터 $\xi$를 통해 산정되며, $\pi (\theta \mid \xi )$는 다양한 $\xi$ 값에 대하여 추정치로서의 가능성을 알아볼 수 있으므로 이를 $\xi$ 의 우도 함수(Likelihood function) $L(\theta \mid \xi )$로 표현가능하며, 사전 분포 $\pi \left(\theta \right)$에 비례하게 된다. 즉,

여기서, $\theta$는 확률 변수이고 $\pi \left(\theta \right)$에 종속되어진다. 따라서, 파라미터 $\theta$의 사후확률밀도함수 $\pi (\theta \mid \xi )$는 베이즈 정리에 의해 Eq. (14)와 같이 표현할 수 있다.

여기서, $L\left(\xi \right)={\int }_{}^{}L(\theta \mid \xi )\pi \left(\theta \right)d\theta$는 기준화 정수이며 한계 확률이다.

앞선 설명을 기반으로 베이지안기법은 크게 1)사전확률분포 $\pi \left(\theta \right)$ 정의, 2)신규 데이터 ξ를 사용하여 우도함수 $L(\theta \mid \xi )$ 정의, 3)베이즈 규칙을 적용하여 사전분포 $\pi \left(\theta \right)$를 수정, 4)사후분포 $\pi (\theta \mid \xi )$ 파라미터 업데이트의 순서로 적용되어진다.

위험률 산정을 위해 필요한 파라미터는 $\theta =(\beta ,\Phi ,\epsilon )$로 정의되며, 필수 데이터는 간단하게 ${\overline{\xi }}^{s_k}=({\overline{\delta }}^{s_k},{\overline{x}}^{s_k},{\overline{z}}^{s_k})$로 표현하였다. 따라서, BMMH 모델은 기초 분석을 위한 파라미터와는 별도로 초매개변수(Hyper parameter) $\Phi$에 대한 사전분포 설정이 필요하다.

여기에서 $\alpha \gamma =1$과 $\alpha {\gamma }^2=1/\Phi$에서 감마 분포를 따르는 것으로 가정하게 되며, 각 모델 파라미터의 조건부결합확률밀도를 각각 계산하게 된다. 파라미터 $\theta$의 기본 형태는 Eq. (15)와 같다.

2.3.2 사후분포의 공식화

먼저, $\theta =({\beta }_{1, ...,} {\beta }_{J-1,} \Phi , {\epsilon }^1, ..., {\epsilon }^k)$에 활용되는 우도 함수에 대한 정의가 필요하며, Eq. (16)에 의해 $L(\theta \mid \xi )$의 공식화가 가능하다(^{Kaito et al., 2012}).

다음으로, 미지의 파라미터 $\theta =(\beta ,\Phi ,\epsilon )$에 대한 사전확률밀도함수 $\pi \left(\theta \right)$가 정의되어야한다. 파라미터 $\beta$에 대해, 이 연구는 $\pi \left(\beta \right)$가 다차원정규분포 ${\beta }_i~N_M({\mu }_i,\sum _i^{})$를 따라 사후확률밀도함수와 유사한 함수로 도출된다고 가정하였다. 여기서, ${\mu }_i$는 기댓값(평균)의 벡터이고 $\sum _i^{}$는 M-다차원정규분포를 따르는 사전 분산-공분산 행렬이다. ${\epsilon }^k$의 밀도함수인 $\pi \left({\epsilon }^k\right)$는 Eq. (5)에 설명된 감마 분포를 따른다. 즉, 계층적 과정 $\pi \left({\epsilon }^k\right)=\pi ({\epsilon }^k:\Phi )$와 $\pi \left(\Phi \right)=h(\Phi :{\alpha }_0,{\gamma }_0)$에 의해 도출된다.

관련 수식을 기준으로 $\pi (\theta \mid \overline{\xi })\propto L(\theta \mid \overline{\xi })\pi \left(\theta \right)$이 도출되어진다. 따라서, 사후결합확률밀도함수(Joint Probability Density Function, JPDF)는 Eq. (19)로 표현된다(^{Kaito et al., 2012}).

2.3.3 사후단계에서의 결합확률밀도함수 추정

파라미터 $\theta =(\beta ,\Phi ,\epsilon )$는 사후결합확률밀도함수(JPDF)를 활용한 Eq. (19)를 통해 산출이 가능하다. 하지만, 직접 샘플링을 통해 산출해 내는 것이 어렵기 때문에, 본 모델에서는 비교적 간단한 Gibbs 샘플링 기법과 Metropolis-Hastings (M-H) 알고리즘을 적용하여 사후결합확률밀도함수(JPDF)를 산출하게 된다.

$\beta$에서 ${\beta }_{e_i}(i=1, ..., J-1)$를 빼는 미지의 파라미터 벡터가 ${\beta }_i^{-e}$로 표시되어 있다고 가정하자. 마찬가지로, $\beta$에서 ${\epsilon }_{e_k}(k=1, ..., K)$를 빼면 ${\epsilon }^{-e_k}$가 된다. Eq. (19)에 의해 CPDF의 ${\beta }_{e_i}$와 $\pi ({\beta }_{e_i}\mid {\beta }^{-e_i},\Phi ,\epsilon ,\overline{\xi })$는 Eq. (20)과 같이 정리되어진다.

여기서, ${\overline{\delta }}_{i{e}_i}^{s_k}$는 더미 변수이며, 이는 샘플 $s_k$의 상태데이터 $h\left({\tau }_A^{s_k}\right)=i$의 상태가 Gibbs 샘플링 $e_i$의 이전 상태등급과 같으면 1이고 그렇지 않으면 0이 된다. ${\epsilon }_{e_k}$에 대한 CPDF는 Eq. (21)과 같다.

추정된 $\beta$와 $\epsilon$을 기반으로, CPDF의 초매개변수 $\Phi ,\pi (\Phi \mid \beta ,\epsilon ,\overline{\xi })$는 다음과 같이 공식화된다.

각 매개 변수에 대한 CPDF는 Eq. (20), Eq. (21), Eq. (22)까지 설명되어 있으며, 이는 Gibbs 샘플링을 사용하여 JPDF를 추정하는 데 사용된다.

2.3.4 조건부확률밀도함수에서의 샘플링

CPDF에서 샘플을 추출한 이후 Eq. (20), Eq. (21), Eq. (22)까지의 적용이 불가능하여, M-H알고리즘을 적용하게 된다. 정확하게는 랜덤워크 M-H가 $\pi ({\beta }_{e_i}^{\left(n\right)}\mid {\beta }^{-e_i(n-1)},{\Phi }^{(n-1)},{\epsilon }^{(n-1)},\overline{\xi })$에 적용되는 반면 $\pi ({\epsilon }_{e_k}^{\left(n\right)}\mid {\beta }^{\left(n\right)},{\Phi }^{(n-1)},{\epsilon }^{-e_k(n-1)},\overline{\xi })$와 $\pi ({\Phi }_{}^{\left(n\right)}\mid {\beta }^{\left(n\right)},{\epsilon }^{\left(n\right)},\overline{\xi })$의 샘플링에는 독립적인 M-H가 필요하게 된다. 여기서 $n(n=1, ..., n, n+1, ..., \overline{n})$는 시행 횟수 또는 샘플 수를 의미한다. $n(n=1, ..., \underset{¯}{n})$의 일부는 최대 수렴에 도달하는 데 필요한 번인(burn-in) 샘플을 의미하며, $n(\underset{¯}{n}+1, ..., \overline{n})$의 다른 부분은 파라미터 샘플을 산정하는데 사용된다. 간단히 말해서, 각 파라미터( $\beta ,\Phi ,\epsilon$)는 파라미터 샘플 $n(\underset{¯}{n}+1, ..., \overline{n})$의 평균을 통해 산정된다.

M-H 알고리즘을 간단하게 설명하면, 다음과 같다. M-H 알고리즘은 ^{Metropolis et al.(1953)} 및 ^{Hastings(1970)}에 의해 생성되고 샘플링에 대략적인 분포(제안 분포 또는 밀도)를 사용하여 대상 분포와의 차이를 감소시키게 된다. 충분히 많은 시행 횟수가 부여될 경우, 마르코프 연쇄는 수렴 근처에서 임의의 작업을 계속 수행하게 된다. M-H 알고리즘의 주요 이슈는 이동 확률 $\alpha =\min (r, 1)$를 계산하고, 확률 $\alpha$로 후보군을 수용하는 것이다. ^{Hastings(1970)}는 임의의 전이확률함수 $q({\tilde{\theta }}^{\left(n\right)}\mid {\theta }^{(n-1)})=\mathrm{Pro}b({\theta }^{n-1}\to {\tilde{\theta }}^{\left(n\right)})$를 사용하고 후보군에 대한 허용 확률을 Eq. (23)과 같이 설정하여 Metropolis 알고리즘을 만들었다.

여기서, 제안분포가 대칭이라고 가정하면, $q\left({\tilde{\theta }}^{\left(n\right)}\right|{\theta }^{(n-1)}=q\left({\theta }^{(n-1)}\right|{\tilde{\theta }}^{\left(n\right)})$는 원래의 Metropolis’s 알고리즘으로 돌아간다.

2.3.5 사후분포에 관한 통계량

베이지안추정을 위한 MCMC 방법에서, 파라미터 $\pi (\theta \mid \xi )$의 기대값은 비모수적 접근법에 의해 얻어지며, 여기서 ${\theta }^{\left(n\right)}=({\beta }^{\left(n\right)},{\Phi }^{\left(n\right)},{\epsilon }^{\left(n\right)})$이다. 먼저, 파라미터 샘플은 번인구간 $(n=1, ..., \underset{¯}{n})$과 샘플링구간 $(n=\underset{¯}{n}+1, ..., \overline{n})$의 두 시퀀스로 분리되어야 한다. 그런 다음 기대값을 추정하기 위한 목적과 진단을 수행하기 위한 목적으로만 사용되어져야 한다.

파라미터샘플 데이터 셋은 $\underset{¯}{n}+1$에서 $\overline{n}$으로 설정된 $M$으로 표시되어지며, 파리미터 $\theta$의 결합확률밀도함수 $I\left(\theta \right)$는 Eq. (24)에 의해 추정되어진다.

여기서, 분자 $⌗({\theta }^{\left(n\right)}\le \theta , n\in M)$는 ${\theta }^{\left(n\right)}\le \theta , n\in M$를 만족하는 샘플의 합을 의미한다. 또한, 파라미터 ${\theta }_i$의 사후 분포에 대한 기대 벡터 $\tilde{\xi }\left(\theta \right)$와 분산-공분산 행렬 $\sum ^{~}\left(\beta \right)$는 다음과 표현할 수 있다.

여기서, ${\theta }_{K_{\theta }}=M(J-1)+K+1$이다. $100(1-2\alpha )\%$의 신뢰구간은 샘플 순서 통계량 $({\underset{¯}{\theta }}_{k_{\theta }}^{\alpha },{\overline{\theta }}_{k_{\theta }}^{\alpha })(k_{\theta }=1,\cdots ,K_{\theta })$에서

Eqs. (26a) and (26b)를 이용하여 ${\underset{¯}{\theta }}_{k_{\theta }}^{\alpha }<{\theta }_{k_{\theta }}<{\overline{\theta }}_{k_{\theta }}^{\alpha }$로 정의된다.

MCMC방법은 초기파라메터 값 ${\theta }^{\left(0\right)}$가 사후분포에서의 표본이라는 가정은 포함하지 않는다. 이러한 문제를 해결하기 위해 해당 모델에서는 ^Geweke(1992)이 제안한 검정기법을 활용하여 가설에 대한 검정을 수행하게 된다. 이는 베이지안 마르코프 하자드 모델을 기반으로 하는 다양한 연구(^{Choi et al., 2019}; ^{Han et al., 2013}; ^{Han et al., 2016}; ^{Kobayashi et al., 2012})에 적용되었으며, 이와 관련된 자세한 사항은 ^Geweke(1992)의 연구를 참조하길 바란다.

3. 실증분석

3.1 데이터 개요

본 연구에서는 BMMH모형을 활용하여 도로포장상태지표 중 균열률의 지역별 기대수명을 추정하였다. 지역별 기대수명을 추정하기 위하여 일반국도포장관리시스템에서 제공되는 도로포장모니터링구간의 데이터를 활용하였다.

모니터링데이터는 크게 구간연장, 위치, 관리사무소, 포장상태자료, 교통량(Annual Average Daily Traffic, AADT), 축하중(Equivalent Single Axle Load, ESAL), 유지보수시기에 대한 정보로 이루어져 있다. 본 연구에서는 2007년부터 2017년까지의 데이터를 활용하였다.

모니터링데이터의 조사구간수는 해당조사구간이 유지보수나 공사와 같이 특정한 사유로 조사가 이루어지지 못하거나, 조사구간의 변경이 이루어짐에 따라 조사년도별 총 조사구간이 2,308~2,445구간으로 차이가 있었다. 하지만 본 연구에서는 포장상태의 시계열적인 데이터를 활용하기 위해 포장상태의 데이터가 누적된 구간만을 분석대상으로 선정하였다. 결과적으로 최종 데이터 구간 수는 총 2,173구간이다.

연도별 도로포장의 평균 상태등급은 Table 1에서 보는 바와 같다. 2007년에는 평균 3.8 %에서 가장 최근 2017년에는 평균 7.3 %의 상태를 보이는 것으로 나타나, 전체적으로 일반국도의 도로포장 상태가 악화되고 있는 것으로 분석된다.

BMMH모델의 적용을 위해서는 각 포장상태지표를 상태등급 자료로 변환해야 한다. 따라서, 본 연구에서는 일반국도 포장관리시스템에서의 유지보수기준(^{Choi et al., 2019})을 준용하여 상태등급을 크게 5등급으로 구분하였다(Table 2).

분석에 활용한 독립변수는 크게 교통량변수, 포장강도변수, 기상변수로 구분되어진다. 먼저 교통량 변수는 모니터링구간데이터에서 취득한 각 구간별 교통량 및 축하중 데이터를 각 구간의 차로수로 나누어 정리하였다. 포장강도변수는 모니터링데이터의 구간별 포장두께자료를 활용하여 포장강도지수를 산정하여 데이터를 구성하였다.

기상변수는 선행연구(^{Choi et al., 2019})와의 비교를 위해 방재기상관측장비(AWS)에서 수집된 기상연보 데이터를 활용하였다. 국내 기상청에서 활용하는 기상관측장비의 종류는 종관기상관측장비(ASOS)와 방재기상관측장비(AWS)의 두 가지로 분류된다. ASOS장비는 전국 96개 지점(20년 10월 기준)에서 운영되며 그 지역의 현재 기상을 실시간 제공 및 기상예보에 활용하고 있는 장비이다. 반면, 무인으로 운영되는 AWS장비의 경우에는 실시간 측정, 연산, 저장, 표출 등의 모든 과정을 자동으로 처리하며 전국 주요 관측지역 및 도서, 산악지역 등 528개 지점에 위험기상 예측을 위해 운영 중이다. 본 연구에서는 두 장비 중 AWS장비를 활용한 관측소에서 취득된 데이터를 활용하였다. Fig. 1에서 각 기상관측장비와 도로포장모니터링구간의 위치를 살펴보면, ASOS장비에 비해 약 5배 많은 AWS장비가 도로포장모니터링 구간에 상대적으로 인접한 것을 확인할 수 있다. ASOS장비와 AWS장비에서 수집되는 기상데이터의 차이점은 아래 Table 3과 같다.

본 연구에서 활용한 기상데이터는 선행연구와의 비교를 위해 동일(평균온도, 최고평균온도, 최저평균온도, 총 강수량 데이터)하게 적용하였다.

모니터링구간에서 가장 가까운 AWS장비의 데이터를 적용하기 위해, 구간과 장비간의 직선거리와 해발고도차이를 고려하여 데이터를 매칭시켰다(Arc GIS S/W 활용). 제설제사용량의 경우, 각 지방 국토관리사무소별로 총 사용량을 제공하고 있어 모니터링구간이 속한 국토관리사무소의 1 km당 평균 제설제사용량 값을 적용하였다.

최초로 구성한 설명변수는 교통량, ESAL, SNP, 평균기온, 최고평균기온, 최저평균기온, 총 강수량, 제설제사용량 데이터로 구성하였다. BMMH 모델을 활용한 분석에 앞서, 변수간 상관성분석과 선행연구(^{Choi et al., 2019})를 참조하여 최종 변수를 선정하였다. 변수선정과정에서 교통량, 평균기온, 최고평균기온변수가 제외되었다. 최종 분석데이터는 각 구간의 포장상태지표인 균열률을 종속변수로 하고, 독립변수는 ESAL, SNP, 총 강수량, 최저평균기온, 제설제사용량을 적용하였다. 여기서, 선행연구와 차이가 존재하는 설명변수는 총 강수량과 최저평균기온이며 이외의 변수는 선행연구와 동일한 값이 적용되었다.

Table 4는 지역별 설명변수의 평균값을 정리한 것이다. 지역별 설명변수는 각 상태등급별로 도출된 파라미터 값과 같이 지역별 기대수명을 추정하는데 활용되어 진다(Eq. (11), Eq. (12a), Eq. (12b)). 지역별 설명변수를 살펴보면 ESAL, SNP, 제설제사용량은 데이터의 전처리 과정에서 발생한 차이 이외에 동일하거나 유사한 값을 보여주고 있다. 하지만, 기상변수(강수량, 최저평균기온)에서는 큰 차이가 발생하는 것으로 나타났다.

Fig. 2에서 보는 바와 같이 강수량과 최저평균기온 모두 ASOS가 AWS보다 수치가 더 높은 것으로 나타났다. 강수량은 AWS와 ASOS 모두 부산이 가장 적은 지역으로 나타났으나 AWS에서는 추가로 대전청이 강수량이 적은 지역으로 분석되었다. 반면 익산청은 두 관측소 모두 강수량 가장 많은 지역으로 분석되었다. 최저평균기온이 높은 지역으로 분석된 지역은 강수량 가장 많은 지역이었던 익산청으로 분석되었으며, 가장 낮은 지역은 원주청으로 분석되었다.

3.2 실증분석

본 연구에서는 BMMH 모형을 활용하여 지역별로 균열률의 기대수명을 추정하였다. 연구의 과정을 간단히 요약하면 1)일반국도모니터링자료의 연도별 취합(2007~2017년), 2)공용역전현상 및 유지보수 수행연도 자료 제거, 3)구간별/연도별 설명변수(ESAL과 SNP)의 계산 및 환경설명변수 매칭, 4)포장상태지표의 변환(상태등급기준으로 이산형 변환), 5)BMMH 분석자료 가공 및 모형분석 6)최종 분석결과 해석 순서와 같다(^{Choi et al., 2019}).

분석을 위해 마련된 데이터는 총 21,340세트이다. 하지만, 데이터 이산화 과정에서 유지보수시행, 데이터 측정오차로 인하여 포장상태가 양호해지는 데이터와 설명변수 자료가 누락된 데이터세트는 모두 제외하였다. 최종적으로 분석에 활용된 데이터 세트는 13,184세트로 기존 연구성과(^{Choi et al., 2019})의 12,005세트에 비해 약 5.5 % 포인트 데이터 세트가 증가한 것으로 나타났다. 해석해보면 ASOS장비에 비해 AWS장비 수가 많아, 기상 데이터 수집율이 상대적으로 높아졌음을 의미한다. Table 5는 BMMH 모델 적용을 위해 AWS의 설명변수 데이터에 대한 표준화를 수행한 결과이며, ASOS의 설명변수에 대한 표준화 결과는 ^{Choi et al.(2019)}를 참조하기 바란다.

파라미터( ${\beta }_i$)의 추정을 위한 반복(Iterations) 횟수는 각 20,000회, MCMC기법의 적용을 위한 초기 랜덤워크(Random walk)는 0.1로 설정하였다. 파라미터의 수렴을 위해 2~3차례 위 과정을 반복하였으며 파라미터의 수렴여부는 Geweke 검정을 통해 수행하였다.

BMMH모형을 통해 추정된 파라미터( ${\beta }_i$)와 이질성 요인( ${\epsilon }^k$)은 Table 6과 같다. 파라미터는 각 설명변수( $x_m$)의 특성에 따라 기대수명 추정에 적합 또는 부적합한 파라미터로 구분된다. ESAL, 강수량, 제설제사용량은 값이 증가(+)하는 경우, 반면에 SNP와 최저평균기온은 값이 감소(-)하는 경우 포장의 기대수명이 줄어드는 특성이 있다. 따라서 Eq. (10)에 의해 각 설명변수의 특성과 상반되는 ESAL, 강수량, 제설제사용량 파라미터가 음수인 경우와 SNP, 최저평균기온 파라미터가 양수인 경우에는 기대수명 추정에 부적합한 파라미터로 판단할 수 있다(^{Choi et al., 2019}).

Table 4에서 확인하였던 각 관리청별 기상변수(강수량, 최저평균기온)의 차이는 파라미터와 이질성 요인에도 영향을 끼친다. Table 6을 살펴보면, 최저평균기온( ${\beta }_4$)의 경우에는 두 시나리오의 1, 2등급 구간에서 동일하게 유효한 변수로 채택되었다. 하지만, AWS 기반 분석에서는 ASOS 기반과 다르게 1등급 구간에서 강수량( ${\beta }_3$)이 추가로 영향을 주는 것으로 분석되었다. 또한 앞선 Table 4에서 분석된 기상변수의 관리청별 차이는 이질성 요인 결과의 차이에서도 도출되었다. 이는 기상변수의 차이가 이질성 요인의 결과에도 영향을 주는 것으로 판단된다.

BMMH 모델을 통한 파라미터 추정 결과 ESAL, SNP, 강수량, 최저기온이 유효한 설명변수로 나타났고 제설제사용량의 경우에는 반복되는 시뮬레이션에서 파라미터의 값이 음수로 나와 유효하지 않은 변수로 분석되었다. ^{Choi et al.(2019)}의 연구와 비교할 경우, 파라미터 및 등급별로 영향을 미치는 설명변수는 유사하게 도출되었다. 하지만, 1등급에서의 기대수명 추정 시에 강수량 변수가 유효한 것으로 나타나, 상대적으로 인접한 지점의 기상장비 데이터를 활용하는 것이 기상에 의한 영향을 고려할 수 있는 것으로 판단된다.

다음으로 BMMH모델의 핵심인 이질성 요인에 대한 분석결과는 다음과 같다. Table 6을 살펴보면 각 관리청별로 값에는 다소 차이가 있으나, 각 지방관리청 간의 전체적인 추세는 동일한 것으로 나타났다. Fig. 3은 MCMC과정에서 얻은 이질성 인자( ${\epsilon }^k$)의 분포를 비교한 것이다. ASOS데이터에 의한 분석에서는 대전과 부산지역의 이질성 인자는 분포 적으로 큰 차이가 없었다. 하지만, AWS데이터를 기반으로 한 분석에서는 대전 지역의 이질성 인자 값이 약 5 % 감소함에 따라, 지역별 파손특성이 ASOS대비 상대적으로 명확히 나타났다. 이것은 AWS장비가 많이 설치되어있어, 각 지역 또는 구간별 기상학적 특징이 세부적으로 구분되었기 때문으로 판단된다.

관리청별 기대수명 값을 추정하기 위해 Eqs. (27) and (28)을 활용하였다. 여기서, 분석에 활용된 변수는 Table 5의 설명변수 값( $x_i$)과 Table 6의 파라미터( ${\beta }_i$) 및 이질성 요인( ${\epsilon }^k$)이다.

먼저, 지역에 상관없이 일반국도 도로포장 전체(bench mark)에 대한 균열률의 기대수명은 14.01년으로 나타났다. 여기서, 전체(bench mark) 기대수명이란 일반국도 도로포장을 대표하는 기대수명으로 해석된다. 지역별 기대수명은 추정된 전체 기대수명을 기준으로, 지역별 이질성요인( ${\epsilon }^k$) 및 설명변수 값( $x_i$)을 활용하여 산정되어지며 최소 13.09년에서 최대 19.61년의 기대수명이 추정되었다(Fig. 4). 대표 기대수명과 비교하여, 익산청의 기대수명이 약 0.92년 정도 낮은 것으로 나타났으며, 원주청은 약 5.6년 정도 기대수명이 높은 것으로 나타났다. ^{Choi et al.(2019)}이 추정한 누적기대수명과 비교할 경우, 누적 기대수명의 차이는 네트워크 전체는 약 0.5년, 각 지역별로 살펴볼 경우 최소 0.4년에서 최대 1.2년 차이가 발생하는 것으로 나타났다.

이러한 기대수명의 차이는 설명변수의 채택 여부에 의해서 발생한 것으로 판단된다. 1등급에서의 기대수명은 ASOS 분석 결과 대비, 평균 157 % 증가한 것을 확인할 수 있다. 반면에, 2~4등급에서의 기대수명은 ASOS대비 97~100 % 수준으로, 감소하거나 유사한 것으로 나타났다(Table 7).

이러한 분석 결과는 Table 6에서의 파라미터 분석 결과를 살펴보면 해석이 가능하다. ASOS 분석 결과에서는 강수량 변수에 해당하는 파라미터( ${\beta }_3$)값이 음수(-)로 산출되어, 강수량 변수가 유효한 변수로 반영되지 않았다. 따라서 AWS 데이터를 활용한 기대수명이 상대적으로 길게 측정된 것으로 판단된다. 반면에, 등급별로 동일한 설명변수가 반영된 2~4등급의 기대수명은 유사한 것을 확인할 수 있다.

본 연구에서 활용한 BMMH모델의 경우 설명변수가 기대수명 추정에 적합한지에 대한 여부는 판별할 수 있으나, 기대수명 추정에 미치는 기여도 등을 판단하기는 어렵다. 따라서 환경변수가 기대수명추정에 미치는 기여도 등을 판별할 수 있는 방안은 후속 연구에서 다루기로 한다.

파손 속도의 개념으로 살펴볼 경우, 전체 네트워크 대비 익산청의 도로포장구간은 파손 속도가 가장 빠르며, 반대로 원주청은 가장 느린 것으로 나타났다. 유지보수 전략 수립 측면에서는 익산청의 도로포장은 균열률에 의한 파손에는 내구성이 약한 것으로 나타나, 균열을 억제하기 위한 특수포장공법의 적용이 바람직할 것으로 판단된다.

4. 결 론

본 연구에서는 BMMH모델을 활용하여, 각 지방국토관리청별 일반국도 아스팔트포장의 균열률에 대한 기대수명을 예측하였다. 분석을 위해 일반국도 포장관리시스템에서 제공(2007~2017년)하는 도로포장 모니터링구간의 균열률 지표를 종속변수로 활용하였다. 독립변수는 크게 교통량변수, 포장강도변수, 기상변수를 적용하였다. BMMH모델의 분석과정에 의해 최종적으로 선정된 독립변수는 축하중, 포장강도지수, 최저평균기온, 총 강수량, 제설제사용량으로 총 5가지 변수이다. 분석결과 각 지방국토관리청별로 균열률에 대한 아스팔트 도로포장의 기대수명을 추정하였으며, 본 연구에서 얻은 성과는 크게 다음과 같다.

(1) 일반국도 아스팔트 도로포장의 균열률에 대한 기대수명에 대한 대푯값은 14.01년으로 나타났다. 각 지방국토관리청별로 살펴보면 대전청 15.09년, 서울청 16.42년, 부산청 14.13년, 익산청 13.09년, 원주청 19.61년으로 나타났다. 대표 기대수명과 비교하여 익산청의 기대수명이 약 0.92년 정도 낮은 것으로 나타났으며, 원주청은 약 5.6년 정도 기대수명이 높은 것으로 나타났다. 이와 같이, 본 연구에서 제안한 이질성요인 계수의 적용을 통해 기대수명의 대푯값을 기준으로 지역적 파손 특성과 차이에 대한 분석이 가능함을 확인하였다.

(2) 동일한 모델을 활용한 선행 연구(^{Choi et al., 2019})와 비교결과, 각 지방국토관리청별 균열률에 대한 기대수명은 유사한 추세를 보이는 것으로 나타났다. 하지만, 선행 연구와 달리, 총 강수량 변수가 기대수명 추정에 영향(도로포장상태등급 1등급 구간의 기대수명 증가)을 주는 것으로 분석되었다. 결과적으로는 선행연구에 적용된 ASOS데이터 대비, 상대적으로 관측소의 수가 많고 모니터링구간과 인접한 AWS 데이터를 활용함으로써 기상에 의한 영향을 고려할 수 있는 것으로 나타났다.

본 연구에서는 설명변수(ASOS와 AWS데이터)의 차이로 인한 기대수명 비교를 목적으로 하였기 때문에, 종속변수를 균열률로 한정하였다는 한계점이 있다. 균열률만을 고려하여 기대수명을 추정할 경우, 다른 분석지표(소성변형, 종단평탄성 등)를 고려하지 않아 기대수명이 실제 수명보다 길게 추정될 수 있다. 또한, 선행연구와의 비교를 목적으로 수행되었기 때문에 적설량, 일조량 등과 같은 다양한 기상변수 고려하지 못하였다는 한계점이 있음을 밝힌다. 따라서, 향후 연구에서 소성변형, 종단평탄성 지표에 대한 기대수명을 추정함에 있어 다양한 기상변수에 대한 고려가 필요할 것으로 판단된다.

본 연구에서 제시한 방법론을 활용할 경우, 각 지역의 특성을 고려한 기대수명 추정이 가능해짐에 따라, 신뢰성 높은 장래유지보수비용 추정을 통한 기반시설관리계획 수립이 가능할 것으로 기대된다.