1. 서 론

강풍의 위력은 내습할 때마다 다르고 예측이 어려워서 공학과 기술이 발전한 지금에도 이에 의한 사고와 재해는 건설 관련 엔지니어들에게 큰 고심을 가져오게 한다. 엔지니어들은 구조물의 안전을 확보하기 위하여 설계 풍하중을 산정하여 내풍설계를 수행한다. 일반적으로 설계 풍속은 시공될 구조물 인근에서 관측된 장기 풍속자료의 극치분석을 통하여 추정한다.

극치풍속을 추정하는 연구는 연구자들이 해결하고자 하는 목적에 맞추어 다양하게 진행되어 왔다. 이러한 극치자료를 추정하는 방법에는 크게 3가지의 방법으로 설명할 수가 있다. 첫 번째는 자료의 기간이 충분할 때 Generalized Extreme Value (GEV) 분포 즉, Gumbel, Frechet, Weibull 분포의 형태 중 하나를 이용하는 방법이 있다(^{Jenkinson, 1955}). 그리고 풍속자료의 극치분석은 주로 Gumbel 분포를 이용하여 추정해왔다(^{Simiu and Scanlan, 1996}). 그러나 자료의 기간이 충분하지 않은 경우에는 The Method of Independent Storms (MIS)와 Peak Over Threshold method (POT) 방법을 이용한다. 이러한 방법들은 일정한 임계값(Threshold) 이상의 자료를 모두 취하여 빈도해석에 이용함으로써 극치분포의 우측꼬리에 추가적으로 정보를 제공하게 되고 이에 따라 매개변수와 빈도에 따른 값을 보다 적절하게 추정할 수 있다. MIS 방법은 각각 다른 태풍을 독립적으로 나누어 그 때의 최대풍속을 극치풍속자료로 사용하여 Gumbel 분포로 추정한다(^{Cook, 1982}). 그리고 POT 방법은 지정한 임계값을 초과하는 극치자료에 대해서 Pareto 분포 또는 Generalized Pareto 분포를 이용하는 방법으로 수자원공학 분야에서 많이 활용하고 있다(^{Palutikof et al., 1999}).

이전 대부분의 연구들은 시간에 따라 변화가 없는 정상성을 가정하여 수행하였다. 하지만 지구 온난화로 인하여 기상변동성이 증가하고 이에 따라 지표풍속, 강수량, 온도와 같은 기후 요소들이 변화하고 있다는 보고나 연구들이 발표되고 있다(^{Zwiers and Kharin, 1998}; ^{Yan et al., 2006}).

최근 누적된 최대풍속 관측자료에 경향이나 변화가 있는지 분석하는 연구들이 진행되었다. ^{Yan et al.(2006)}은 1958~1998년 동안의 NCEP 재해석자료와 Gamma 분포와 Weibull 분포 모델을 이용하여 북서 유럽의 일 최대풍속 시계열을 모의하였다. ^{Hundecha et al.(2008)}은 북아메리카 지역의 NARR 재해석자료와 기상 관측자료로부터 연 최대풍속의 비정상성 극치분석을 통하여 캐나다의 St. Lawrence만 주변의 풍속 변화를 조사하였다. ^{Wang and Li(2016)}는 미국 Joint Typhoon Warning Center에서 제공하는 과거 71년 동안의 북서태평양 열대저기압 자료를 이용하여 기후 변화로 인한 열대저기압의 강도가 증가한다고 확인하였다. 그리고 남중국해와 북서태평양 지역의 극치풍속을 추정하기 위한 비정상성 모델을 제안하였다. 하지만 국내에는 비정상성을 가지는 극치풍속의 추정에 대한 연구는 상대적으로 미비하다.

이에 본 연구에서는 경향성을 반영할 수 있는 비정상성 빈도해석을 수행하기 위하여 Bayesian 모형의 매개변수들에 대한 선형회귀의 모형식을 수립하였다. 그리고 매개변수들의 추정에 있어 Markov Chain Monte Carlo (MCMC) 모의를 수행하였으며 제주공항의 연 최대풍속 자료를 대상으로 비정상성 빈도해석의 적합성을 평가하였다.

2. 극치풍속의 추정 방법

본 연구에서는 경향성을 반영할 수 있는 비정상성 빈도해석을 수행하기 위하여 Bayesian 모형에 매개변수들의 경향성을 반영할 수 있는 선형회귀의 모형식을 수립하였고, MCMC 모의를 통해서 Bayesian 모형의 매개변수들을 추정하였다.

2.1 Bayesian 방법

Bayesian 확률은 사전확률을 기준으로 자료가 증가되어 정보가 갱신되고, 최종에는 조건부로 바뀌는 사후확률을 추정하는 과정이다. Bayes 정리(Bayes’s rule)는 두 확률 변수의 조건부(Conditional) 확률분포와 주변부(Marginal) 확률분포를 연관 짓는 확률이론을 말한다. 그리고 확률이론에서의 Bayesian 해석에 따르면 Bayes 정리는 사전확률과 추가적인 정보를 통해 사후확률을 구할 수 있다(^{Gelman et al., 2003}).

확률변수 $\theta$와 $y$가 주어졌을 때, 이 확률변수들의 결합확률분포를 생성하는 모형에 대해 결합확률밀도함수는 사전분포 $p\left(\theta \right)$와 자료분포 $p\left(y\overline{)\theta }\right)$의 곱으로 Eq. (1)처럼 나타낼 수 있다. Bayes 정리에서 확률사건 $\theta$와 $y$의 조건부 확률과 경계확률의 관계는 Eq. (2)와 같다. 여기서, $p\left(\theta \overline{)y}\right)$는 $y$가 주어졌을 때 $\theta$의 조건부 확률(사후확률)이라 하고 $l\left(\theta \overline{)y}\right)$는 $y$가 고정되었을 때 $\theta$의 우도(likelihood)이다. 즉, Bayesian 추론에서 사후확률의 개념은 사전확률과 우도의 곱에 비례한다는 점이다. 이것을 앞에서 언급한 Bayes 정리의 각 항을 이용하여 나타내면 Eq. (3)과 같이 표현할 수 있다.

$$p(\theta,y)=p(\theta)p(y\vert\theta)$$

(1)

$$p(\theta\vert y)=\frac{p(\theta,y)}{p(y)}=\frac{p(\theta)p(y\vert\theta)}{p(y)}\propto l(\theta\vert y)p(\theta)$$

(2)

$$p(y)= \int _{} ^{} {p(y, \theta )d \theta =} \int _{} ^{} {p( \theta )} p(y | \theta )d \theta$$

(3)

2.2 Markov Chain Monte Carlo 모의

Monte Carlo 기법은 수리적으로 계산이 불가능하거나 복잡한 적분, 추정 등의 문제에 활용되고 응용통계학 분야에서 매개변수의 추정 및 불확실성을 평가하는 수단으로 이용되고 있다. 이 기법이 확률변수들 간의 독립성을 가정으로 이루어지는 샘플링 방법이라면, MCMC (Markov Chain Monte Carlo) 기법은 확률변수들 간의 종속성을 기준으로 조건부 샘플링이 가능한 방법이다. 난수를 발생시킨다고 하여 관심이 되는 확률분포를 정확하게 따르지는 않지만, 일정 시간동안 반복 후 얻는 난수들은 추정하려는 분포에 수렴하게 된다. 따라서 MCMC 기법은 복잡한 확률변수들 간의 확률분포 및 매개변수의 추정을 요하는 문제에서 주로 사용된다(^{Kwon et al., 2008}).

MCMC 기법의 대표적인 방법으로 메트로폴리스-헤스팅 알고리즘(Metropolis-Hastings algorithm)과 깁스표본법(Gibbs sampling)이 있다. 본 연구에서는 깁스표본법을 이용하였다. 깁스표본법은 원하는 다변량의 확률분포에서 독립적이며 동일한 분포(iid)를 따르는 표본추출이 복잡하거나 어려운 경우에 이용할 수 있는 방법이다. 2개의 변수를 갖는 다변량의 확률분포를 이용하여 설명하면 다음과 같다. 2개의 변수 $\mu$와 $\delta$를 갖는 다변량의 확률밀도함수를 $f(\mu ,\sigma )$라고 정의한다. 깁스표본법은 확률밀도함수로부터 직접 표본을 추출할 수는 없으나 각각의 변수들에 대해서 다른 두 변수들이 주어졌을 때의 조건부 확률분포를 알고 이로부터의 표본추출이 가능한 경우에 사용할 수 있다(^{Kwon et al., 2012}). 알고리즘을 정리하면 다음과 같다(^{Lee, 2010}).

[1] Gumbel 분포의 매개변수들에 대한 초기값 $({\mu }^{\left(0\right)},{\sigma }^{\left(0\right)})$을 부여한다.

[2] $i$번째 난수 벡터 $({\mu }^{\left(i\right)},{\sigma }^{\left(i\right)})$가 주어졌을 때 $i+1$번째 난수 벡터를 다음과 같은 조건부 분포에서 추출한다.

(1) ${\mu }^{i+1}~f\left(\mu \right|{\sigma }^{\left(i\right)})$

(2) ${\sigma }^{i+1}~f\left(\sigma \right|{\mu }^{(i+1)})$

[3] 위의 과정을 충분히 반복한 후 초기의 일정부분 난수를 제거한 이후의 난수들을 이용한다. 이러한 제거과정을 모의담금질이라고 하며 Bayesian해석에서 일반적으로 요구되는 단계이다(^{George and McCulloch, 1993}).

3. 비정상성 빈도해석

3.1 대상지점 및 자료

강우 자료와 달리 풍속 자료는 주변 지형의 영향을 많이 받는다. 즉, 주변의 지표 조도가 바뀌면 풍속 또한 달라진다. 가장 대표적인 경우가 도시화로 인하여 기상관측소 주변에 높은 건물들이 들어서면서 병풍을 두른 효과가 있다. 그리고 기상관측소가 이전한 경우에도 지표 조도가 크게 바뀌게 된다. 따라서 경향성을 분석하기 위해서는 자료의 연속성을 보장할 수 있는 풍속자료를 수집해야 한다.

공항을 둘러싼 주변 고도는 항공기의 안전한 이륙 및 착륙을 위해서 제한된다. 그리고 안전한 항공기의 이착륙을 위하여 공항기상관측장비(ASOS)을 사용하여 실시간으로 기상 정보를 제공하며 풍속을 관측한다. 따라서 공항에서 관측한 풍속자료는 지표 조도의 변화가 거의 없어서 풍속의 비정상성을 분석하기에 적절하다.

본 연구에서는 제주도 해안에 위치하고 있어 오랜 기간 동안 지표 조도의 변화를 주는 요소가 없는 제주공항의 1984년 1월에서 2017년 12월까지 관측된 10분 평균풍속 자료를 이용하였다. 또한 하루 동안에 관측된 자료들 중 최댓값을 일 최대풍속으로 정하고 같은 방식으로 연, 월에 해당하는 최대풍속을 추출하였다. 이 지역의 최대풍속을 발생시키는 요인은 여름철의 태풍과 겨울철의 계절풍이 있는데 제주 지역은 빈번하게 태풍의 경로에 놓여 있다. 따라서 연 최대풍속에 지배적인 요인은 태풍이다. 그리고 제주공항 지점의 풍속계에 대한 변경 이력은 기상청에서 제공하는 기상월보를 통해 확인하였고 Table 1과 같다. 풍속계 높이가 변경된 이력이 있는 기간의 풍속은 보정을 하였다. 지표면 가까이에서는 조도에 의한 마찰력이 풍속에 가장 지배적이므로 태풍 풍속의 고도 보정에도 지수법칙의 사용이 가능하다. 보정 방법으로 Eq. (4)와 같은 지수법칙을 이용하였다.

Table 1. Changed History of Wind Gauge

Period	Height
1984.01 ~ 1985.12	10 m
1986.01 ~ 2003.01	15 m
2003.02 ~ 2017.12	10 m

$$V(z)=V(z_r)\lbrack\frac z{z_r}\rbrack^\alpha$$

(4)

여기서, $V\left(z\right)$는 높이 $Z$에서 평균풍속(m/s), $V\left(z_r\right)$은 기준 높이 $z_r$에서 평균풍속(m/s)이고 $z_r$은 10 m를 기준으로 하는 경우가 많다. a는 지표면의 조도에 따라 결정되는 지수이며, 개활지에 해당하는 0.16을 사용하였다.

3.2 경향성을 고려한 빈도해석 모형

Bayesian 기법을 통한 매개변수 추정은 기존의 최우도법이나 모멘트법과는 다르게 매개변수를 하나의 확률변수로 취급한다. 이는 매개변수가 하나의 값이 아닌 확률분포의 형태로 부여되며 결국 매개변수의 사후분포를 추정하는데 목적을 두고, Eq. (3)에 기반한다.

수집된 풍속자료를 임의의 확률분포로 가정할 경우에는 그 확률분포의 유용성을 검토하기 위한 적합도(Goodness of fit, GOF) 검정이 필요하다. 적합도를 검정하기 위하여 χ² 검정(Chi-Square test)와 K-S 검정(Kolmogorov-Smirnov test)을 5 %의 유의수준에서 수행하였다. 그리고 빈도해석 모형으로 더 알맞은 확률분포를 찾기 위하여 부정로그우도(Negative Log Likelihood, NLL) 함수와 BIC (Bayesian Information Criterion) 값을 이용하였다. Table 2에는 각각의 적합도 검정 평가 결과를 나타내었다. χ² 검정과 K-S 검정 결과는 값이 0인 경우에 적합을, 1인 경우에 부적합을 의미한다. 값의 대부분이 0인 것으로 보아 가정한 분포들은 적합하다고 할 수 있다. 그렇지만 χ² 검정에서 정규 분포는 부적합하다고 판정되었다. NLL과 BIC 값의 크기가 가장 작은 확률 분포는 Gumbel 분포이므로 이 분포를 빈도해석의 확률분포 모형으로 선택하였다.

$p\left(\theta \right|y)$는 사후분포로 제주공항 풍속자료의 모집단이며 Gumbel 분포를 따른다고 가정한다. 정상성 모형 기준으로 $\theta =[\mu ,\sigma ]$는 Gumbel 분포에 대한 매개변수들의 집합을 나타낸다. $p\left(y\right)$는 극치풍속자료 $y$의 주변분포(Marginal), $p\left(\theta \right)$는 매개변수들의 사전분포를, $p\left(y\right|\theta )$는 극치풍속자료 $y$의 우도함수를 의미하고 Eq. (5)와 같이 나타낼 수 있다. 여기서, N은 연 최대풍속 자료의 개수를 나타내며, 위치매개변수와 규모매개변수는 사전분포로 각각 정규분포(N)와 감마분포( $\Gamma$)로 가정하였다.

Table 2. Goodness of Fit Tests and Bayesian Information Criterion Value for Probability Distributions

Distributions	Gumbel	GEV	Gamma	Normal	Weibull
χ2 test	0	0	0	1	0
K-S test	0	0	0	0	0
NLL	-82	-82	-85	-87	-90
BIC	172	175	177	181	187

$$p(\theta\vert y)=\Pi_{t=1}^NGumbel(y_t\vert\mu,\sigma)\bullet N(0.1,0.1)\bullet\Gamma(1,0.1)$$

(5)

$$\begin{array}{l}p(\theta\vert y)\propto\Pi_{t=1}^NGumbel(y_t\vert\mu(t),\sigma(t))\bullet N(\mu(t)\vert a+bt,\sigma_\mu)\\\bullet\Gamma(\sigma(t)\vert c+dt,\sigma_\sigma)\end{array}$$

(6)

본 연구에서 제시하는 비정상성 빈도해석은 분포의 두 매개변수의 선형경향성을 고려하였고, Eq. (6)처럼 확장할 수 있다. Eq. (6)에서 a와 b는 위치매개변수와 c와 d는 규모매개변수의 선형경향성을 고려하기 위한 회귀계수를 의미한다. 해석적으로 모든 매개변수들에 대한 사후분포를 직접적으로 추정하는 것은 불가능하다. 따라서 본 연구에서는 MCMC방법 중 깁스표본법을 이용하여 매개변수들의 사후분포를 추정하였다.

3.3 비정상성 빈도해석 결과

Gumbel 확률지에 제주공항의 극치풍속을 도시한 결과는 Fig. 1(a)와 같고, 좌·우측꼬리 부분의 거동을 잘 묘사하고 있다. Fig. 1(b)에서는 경험적인 누가확률밀도함수(CDF)와 Gumbel 분포의 누가확률밀도함수가 유사한 거동을 하고 있음을 확인할 수 있다.

Fig. 2는 극치풍속자료를 Normal 분포로 가정한 경향선(사각형)과 Gumbel 분포(원형)로 가정한 경향선을 비교한 결과이다. 경향선의 기울기가 Gumbel 분포보다 Normal 분포가 큰 것을 확인할 수 있다. 중앙값 외에 5 %, 25 %, 75 % 그리고 95 %의 불확실성 구간(Credible Interval; CI)을 사후분포로부터 추정하여 Fig. 2에 표시하였다. Fig. 2에서 보는바와 같이 불확실성 구간이 중앙값을 기준으로 비대칭하게, 우측꼬리(upper tail)가 좌측꼬리(lower tail)보다 두껍게 나타내는 것을 알 수 있으며 이는 Gumbel 분포로 극치자료의 경향성을 나타낸 결과라고 할 수 있다.

도로교 설계기준(^{Korea Road & Transportation, 2010})과 건축 구조기준(^{Architectural Institute of Korea, 2016})에서는 교량이나 건축물의 설계 풍하중을 산정하는데 있어서 재현기간 100년에 해당하는 기본풍속을 주로 이용하기 때문에 Eq. (6)을 통해서 추정된 분포의 매개변수를 활용하여 재현기간 100년에 대한 기본풍속을 추정한다. 그리고 결과는 Fig. 3과 같다. Fig. 3은 경향성을 고려하여 추정된 100년 빈도 풍속의 불확실성 구간을 사후분포로부터 추정하여 나타내었다. 경향성이 존재한다고 가정하면 1998년 이후부터 현재의 기본풍속을 초과하고 있으며, 결과적으로 정상성을 기준으로 추정되는 기본풍속은 시간이 경과됨에 따라 과소추정 될 수 있음을 의미한다.

Bayesian 모형은 모든 매개변수의 사후분포로 통계적인 추론이 가능하여 추정된 매개변수의 유의성 여부를 정량적으로 파악할 수 있다. Fig. 4는 극치풍속자료가 경향성을 가지고 Gumbel 확률분포를 따를 때 위치 및 규모매개변수에 대한 기울기 b와 d의 확률분포이다. 이는 기울기 b와 d의 값이 각각 나타날 확률을 의미하고, 기울기의 값이 음수이면 감소하는 경향을 양수이면 증가하는 경향을 갖는다. Fig. 4(b)에서 기울기 b가 0.12 m/s 근처에서의 확률이 큰 것으로 보아, 극치풍속이 해마다 0.12 m/s씩 증가할 확률이 높다는 것을 의미한다. 또한 수치적분을 통하여 Fig. 4(b)의 면적을 구함으로써 기울기의 구간에 대한 확률을 산출할 수 있다. Fig. 4(b)에서 양수 구간에 대해 수치적분을 하면, 증가하는 경향에 대한 확률을 추정할 수 있고 b는 0.989를 d는 0.999의 값을 얻었다. 기울기에 양수가 나타나지 않을 확률과 그 확률을 허용할 수 있는 기준인 유의수준을 비교하여 통계적으로 유의한지 아닌지를 검정할 수 있다. 검증 결과를 보면, 기울기 b와 d가 음수일 확률은 각각 0.011과 0.001이고, 95 % 신뢰도 기준의 유의수준인 0.05보다 작다. 즉, 극치풍속자료의 경향성이 증가하는 것은 통계적으로 유의하다고 판단할 수 있다.

Fig. 5(a)는 경향성이 고려된 재현기간 동안의 기본풍속을 5년 간격으로 추정하여 나타낸 빈도곡선이고 빨간색 점선은 기존 방법과 같이 정상성을 가정하여 추정한 기본풍속이다. 1998년 이후로 비정상성 빈도해석 결과값을 초과하는 것을 확인할 수 있다. Fig. 5(b)는 정상성 빈도해석과 비정상성 빈도해석을 이용하여 추정한 기본풍속과 신뢰구간이다. 재현기간 100년일 때 기본풍속의 차이는 약 1 m/s로 비정상성 빈도해석 결과로 추정한 값이 더 크다. 이는 극치풍속자료에 경향성이 존재하는 경우에는 기본풍속이 과소추정 될 가능성이 있다는 것을 의미한다.

정상성 빈도해석 과정에서 추정된 Gumbel 확률밀도함수와 시간에 따라 연속적으로 추정된 확률밀도함수를 비교한 것은 Fig. 6과 같다. 현재 극치풍속의 확률밀도함수뿐만 아니라 현재의 경향성이 미래에도 계속된다는 가정으로 확률밀도함수를 추정하였다. 현재(파랑선)와 미래(검정파선)의 분포를 보면 우측꼬리가 경향성 증가로 인해 두꺼워지는 것을 확인할 수 있다. 이는 우측꼬리 부분에 속하는 극치풍속의 발생확률이 미래에는 증가할 수 있다는 것을 의미하며, ^{Choi and Kim(2007)}이 우리나라에 영향을 미치는 태풍의 특성과 이동경로에 따른 빈도를 분석하여 태풍의 빈도가 증가한다고 제시한 결과와 비슷하다. ^{Choi and Moon(2008)}은 56년 동안의 기상관측 지점의 풍속자료를 이용하여 56년간 연 최대순간풍속이 15 m/s정도 증가함을 확인하였다. 증가하는 경향을 보이는 본 연구결과와 비슷하다. 규모매개변수의 불확실성 구간이 넓어지는 것을 보여준다. 경향성을 고려한 비정상성해석결과와 매개변수들의 사후분포에 대한 정보는 Table 3에서 확인할 수 있다. Table 3에서 2.5 %, 97.5 %는 유의수준을 5 %로 고려했을 때, 사후분포의 불확실성구간을 표현하기 위한 경계값이다.

Table 3. Posterior Density of Parameters and Uncertainty Bounds

Parameter	Mean	SD	2.5 %	50 %	97.5 %
Stationary $\mu$	18.11	0.43	17.22	18.12	18.92
Nonstationary $\mu$	18.41	0.67	16.35	18.41	20.46
Stationary $\sigma$	2.35	0.34	1.79	2.32	3.11
Nonstationary $\sigma$	2.54	0.49	1.64	2.54	3.44
$\mu$ intercept - a	16.18	0.88	14.43	16.24	17.80
$\mu$ slope - b	0.13	0.05	0.03	0.12	0.22
$\sigma$ intercept - c	1.57	0.51	0.67	1.49	2.69
$\sigma$ slope - d	0.06	0.03	0.01	0.05	0.13

4. 결 론

본 연구에서는 수공학 분야에서 활용되고 있는 비정상성 해석 기법(^{Lee et al., 2010)}을 극치풍속 자료에 적용하여 선형 경향성을 고려한 비정상성 빈도해석을 수행하고 극치풍속의 경향성을 평가하였다. 제주공항 연 최대풍속을 대상으로 빈도분석을 수행한 결과를 보면, 제주 지역의 극치풍속은 해마다 0.12 m/s씩 증가할 확률이 가장 크고, 재현기간 100년의 기본풍속은 기존의 정상성 빈도해석 결과보다 약 1 m/s정도 크다. 제주공항의 결과에서 보듯이 극치풍속 자료에 경향성이 존재하는 경우에는 과소추정 될 가능성이 있으므로, 통계적으로 유의한 경향성을 가지는 지역의 경우에는 비정상성을 고려한 극치풍속 추정방안의 도입이 요구된다.