공간상에서 얻어진 자료는 서로 멀리 떨어져 있을수록 낮은 상관을 가지며 가까이 있을수록 높은 상관을 가지는 특성을 보인다. 이러한 공간자료의 특성은 공간상에 위치하고 있는 각종 계측 기기로부터 얻어지고 있는 다양한 환경자료에서도 볼 수 있다. 특히, 환경부의 수질측정망에서는 다양한 수질오염자료들이 측정되고 데이터베이스화되어 다양한 수질환경정책의 의사결정에 이용되고 있다. 그러나 수질측정망에서 측정된 자료는 측정되는 지점에 대한 값을 의미하기에 측정망 이외의 지점에서의 수질오염물질에 대한 예측이 필요하다. 이처럼 공간에서 관측된 자료를 토대로 미관측 지점에 대한 예측을 위해 공간분석 기법 중 크리깅(Kriging) 기법이 많이 활용되고 있다. 크리깅 기법은 원래 지질통계 분야에서 널리 이용되었던 기법으로 현재는 환경, 교통, 보건 등 다양한 분야에서 폭넓게 활용되고 있다.

Journel (1984)은 환경자료에 대해 공간통계학적 기법을 적용하고자 하였으며, Istok and Cooper (1988)는 토양의 오염농도에 대한 공간적 예측을 위해 공간분석을 이용하였다. Myers (1984)는 여러 가지 오염문질의 공간적 흐름을 평가하기 위해 공간통계학적 기법을 적용하였다. 또한 강우량, 대기오염 등과 같은 환경자료, 교통자료, 범죄자료 등과 같은 공간상에서 얻어지는 자료에 대해 많은 연구자들이 공간분석을 활용하고 있다(Choi et al., 2010; Heo et al., 2007; Heo et al., 2004; Kim and Choi, 2000; Yang and Jin, 2010).

크리깅 기법을 지하수 수위를 감시하기 위한 공간적 설계에 적용하였으며(Yang et al., 2008), Smith et al. (1997)은 수질 감시 자료에 대해 공간분석을 활용하였으며, 강에 용해되어 있는 주요 화학물질에 대한 공간적 추세 연구(Carleton, 2008) 등과 같은 다양한 연구들이 진행되어 왔다.

Cho et al. (2001)은 공간 상관성이 있는 환경자료에 적용하여 공간 통계 방법이 일반 통계 방법보다 공간 예측력이 좋다는 것을 보였으며, 예측방법으로 공간예측 방법인 범용크리깅(ordinary kriging)방법과 비모수방법인 국소회귀분석방법(loess)과 GAM(Generalized Additive Model)을 비교하였다. Kang et al. (2008)은 다변량 공간자료를 분석하는데 있어서 공동크리깅(Cokriging)기법의 활용성에 대해 연구하였다. Yoo and Um (1999)은 강우강도자료를 이용하여 베리오그램(Variogram) 모형의 모수 추정 방법과 함께 이상치 탐지방법에 대하여 연구하였으며, 최소제곱방법을 사용하여 베리오그램 모형의 모수를 추정하였으며 범용크리깅 방법을 사용하여 공간예측을 실시하였다.

본 연구에서는 금강 유역을 대상으로 수질을 대표할 수 있는 생물화학적 산소요구량(Biochemical Oxygen Demand, BOD)에 대해 베리오그램 모형의 모수를 가중최소제곱방법(method of weigted least squares)에 의해 추정하고, 공간예측을 위해 최적의 모형을 선택하고 공간예측을 위해 범용크리깅과 공동크리깅의 두 가지 방법을 적용시켰다. 베리오그램 모형의 선택기준으로 잔차제곱합(residual sum of squares)을 이용하였으며 공간예측력의 성능 평가를 위해 교차검증방법(Cross-Validation)에 의한 예측잔차제곱합을 계산하였다. 자료의 분석도구로서는 통계패키지 R의 다양한 공간분석 패키지를 사용하였다.

금강 유역은 전북 장수군에서 발원하여 대전, 공주, 부여 등을 경유하여 군산만으로 흘러드는 강으로, 그 유역면적은 9,914.01km²로서 우리나라에서 3번째로 큰 유역이다. 기후는 한반도의 중서부에 위치하여 온대북한계에 가까워 유역의 평균기온은 11.0 ~ 12.5°C이고, 연강수량은 1,100 ~ 1,300mm이나, 지형 등의 요인에 의해 국지적인 변화가 크다(Geum River Basin Environmental Office, 2014). 본 연구에서 대상으로 하고 있는 금강 유역의 62개의 본류 및 지류의 수질측정지점을 상･하류로 구분하여 Fig. 1에 나타내었다.

금강은 침식계곡을 따라서 발달하여 있으며, 풍화와 침식작용이 활발하여 낮은 경사의 저구릉지가 발달하였고 분지형 지형에서 나타날 수 있는 수지상 수계가 발달하여 있다. 금강의 상류는 대부분 변성암으로 구성되며, 중류 및 하류는 중생대 화강암류가 주로 분포되어 있다. 금강 유역의 토지지목별 분포는 임야가 61.2%, 농경지가 27.4%, 시가가 4.6%를 차지한다(Kim, 2006). 또한, 유역 상류에 용담댐, 중류에는 대청댐, 하구에는 하구둑 등 주요 시설물이 설치되어 유역 내 유수를 조절하는 중요한 기능을 담당하고 있다(Kim et al., 2007). 본류로 유입되는 주요 지류 중 상류에는 무주남대천, 영동천, 보청천 등이 합류되며, 중류에는 갑천, 미호천, 유구천, 하류에서는 석성천, 논산천이 합류되어 서해로 흘러가게 된다.

관심공간에서 관측된 공간자료들이 거리에 의하여 영향을 받는다면 공간자료들은 서로 위치가 가까울수록 그 자료들의 값이 높은 상관을 보이고 점점 멀어질수록 낮은 상관을 보인다. 따라서 관심공간 내에서 공간자료들이 서로 독립관계에 있지 않고 상관관계에 있을 때 사용할 수 있는 공간모형 및 모형을 통한 공간적 예측의 이론적인 기반이 필요하다(Isaaks and Srivastava, 1989).

관심공간에 대해서 공간자료들이 연속이라는 가정 하에 공간자료에 대한 공간확률과정(spatial random process)은 식 (1)과 같다.

여기서 Z(s) 는 위치 s 에 대한 함수이며, s 는 공간자료의 위치를 나타내며, 일반적으로 위도와 경도와 같은 좌표를 나타낸다. D 는 공간자료들이 수집되는 관심공간이며, R^d는 d차원의 실수 전체 공간을 나타낸다.

본 연구의 공간변수 {Z(s) : s ∈D는 다음과 같은 정상성(stationarity) 조건을 만족한다고 하면, C(h) 는 공분산 함수로서 거리만의 함수이며, 2γ(h) 와 γ(h) 는 각각 Z(s)의 베리오그램(variogram) 또는 세미베리오그램(semivariogram)이라 한다.

1) E(Z(s)) = μ∀s∈D 2) Cov(Z(s)), Z(s+h)) = C(h) < ∞, ∀s, s+h∈D, h≥ 0 3) Var(Z(s)−Z(s+h)) + 2γ(h), ∀s, s+h∈D, h ≥ 0

베리오그램은 확률장(random field)에서 존재하는 공간적 연관성 및 구조 등에 대한 정보를 전달하는 구조적인 도구이다. 공분산 함수와 베리오그램은 공간예측을 위한 크리깅에서 매우 중요한 역할을 하며, 크리깅 예측 방정식은 베리오그램이나 공분산의 함수로 표현될 수 있으며, 공분산 함수 C(h) 와 세미베리오그램 γ(h)사이에는 γ(h) = C(0) − C(h) 관계식이 성립한다(Cressie, 1993).

이와 같이 공간 자료에 대한 베리오그램의 추정, 추정된 베리오그램을 통한 최적 모형 선택, 선택된 베리오그램을 이용한 공간예측의 단계로 이루어져 있으며 공간분석에 대한 순서도는 Fig. 2와 같다.

베리오그램은 공간분석에서 자주 활용되는 공간 연관성의 정도를 나타내는 도구로 활용된다. 베리오그램 분석은 공간분석을 위한 첫 번째 단계라 할 수 있으며, 공간 자료에 대한 베리오그램 분석은 관측된 자료를 통하여 공간적 변동 및 특성을 알아보기 위하여 널리 사용되고 있다. 베리오그램은 관심공간 상에 분포하고 있는 임의의 두 자료의 차이의 분산으로 다음과 같이 정의된다.

식 (2)에서 정의된 베리오그램에 대해 n개의 공간 자료 Z(s_i),i = 1,⋯,n에 대해, Matheron(1963)이 제안한 (세미)베리오그램 적률추정량(method of moment estimator)은 다음과 같다.

여기서 N(h) = {(s_i,s_j) : |s_i−s_j|=h : i, j= 1,⋯,n}는 유클리드 거리로 h만큼 떨어진 자료의 집합을 나타내고, |N(h)|는 원소의 개수를 나타낸다.

대부분의 세미베리오그램 모형은 뭉치(nugget), 문턱(sill)과 범위(range)와 같은 모수들을 가지고 있다. 공간분석에서 자주 사용되고 있는 세미베리오그램 모형에는 가우시안모형(gaussian model), 지수모형(exponential model), 구형모형(spherical model), 선형모형(linear model)과 파워모형(power model) 등이 있다. 본 연구에서는 가우시안모형, 지수모형과 구형모형을 사용하고자 한다.

가우시안모형, 지수모형과 구형모형에 대한 함수식이 각각 식 (4-6)에 있다.

여기서, θ = (c_n,c_s,c_r)^T, c_n ≥ 0, c_s ≥ 0, c_r ≥ 0을 만족하며, c_n는 뭉치효과를, c_n + c_s는 문턱을, c_s는 부분 문턱(partial sill)을, c_r은 범위를 나타낸다.

세미베리오그램 모형의 모수 중 뭉치효과(nugget effect)는 두 공간위치 간의 거리가 0으로 가까이 갈 때 0이 되기 직전의 세미베리오그램 값으로 정의되는데, 보통 측정오차(measurement error)를 나타낸다. 범위(range)는 공간적으로 더 이상 상관이 존재하지 않는 최소거리를 나타낸다. 자료가 범위 이상의 거리만큼 떨어져 있을 때 세미베리오그램은 일정한 값을 가지는데 이 값을 문턱(sill)이라고 하며 이 값은 자료의 분산과 같다. 세미베리오그램 모형의 모수를 추정하는 방법에는 가장 많이 사용되는 보통 최소제곱방법과 최소제곱방법의 일종인 가중 최소제곱방법이 있으며, 우도방법인 최우추정방법과 잔차최우추정방법 등이 있다(Heo et al., 2004).

베리오그램의 모수를 추정하기 위한 가장 보편적인 방법 은 최소제곱방법으로 최소제곱방법을 통한 베리오그램의 모수를 추정하는 단계는 다음과 같다. 관심공간상의 자료로 부터 베리오그램의 표본 추정값들을 식 (3)과 같이 도출하 고, 표본 추정값들을 식 (4-6)과 같은 이론적 베리오그램 모형에 적합화시킨다(Cressie and Hawkins, 1980). 표본 베 리오그램의 값들 (h₁), (h₂),⋯,(h_k)과 이론적 베리오그램 사이의 제곱합을 최소화시키는 모수벡터(뭉치효과, 문턱과 범위) 추정치를 구한다. 즉,

최소제곱방법은 관심공간 내 자료들이 서로 다른 래그(lag)안에서 중복적으로 사용될 수 있기 때문에 베리오그램의 추정값들이 서로 연관 될 수 있는 단점이 있어 이러한 단점을 극복하기 위해 최소제곱추정방법으로서 가중최소제곱법추정량을 이용하는 것이 바람직하며, 가중최소제곱법추정량은 다음과 같이 간략히 설명될 수 있다.

가중최소제곱방법은 일반화 최소제곱방법의 공분산행렬 V를 (h)의 분산 근사값들의 대각행렬로 대치하는 방법이기에 베리오그램 모형에서 모수의 가중최소제곱추정량은 식 (8)을 최소화하여 구한다.

여기서, 베리오그램 추정값의 분산이 근사적으로 이기에 가중치 행렬 W 는 아래와 같이 표현될 수 있다(Cressie, 1993).

공간분석에서 가장 중요하게 다루는 문제 중의 한 가지는 관측된 자료를 통해 관측되지 않은 미관측 지점에 대한 예측값 및 예측값의 분산을 제공하는 것이다. 공간확률과정 {Z(s) : s ∈D, D⊂ R^d 따르는 n개의 공간 자료를 가지고 공간상의 연관성을 이용하여 미지의 s₀지점에서의 Z(s₀)을 예측하는 것을 크리깅(Kriging)이라 한다.

크리깅은 다양한 통계적 가정에 따라 단순크리깅(simple kriging), 범용크리깅(ordinary kriging), 일반크리깅(universal kriging), 지시크리깅(indicator kriging), 공동크리깅(cokriging) 등으로 구분할 수 있다.

단순 크리깅은 정상성 가정에서 E(Z(s)) = μ(s), s∈D를 안다는 조건하에 사용할 수 있는 가장 간단한 크리깅방법이나 실제 공간자료 적용 시 사전에 μ(s)를 알 수 없는 경우가 대부분이기 때문에 현실성이 결여되어 있다고 할 수 있으며 단순 크리깅 추정 방정식은 불편성을 만족하지 않는다는 특성도 가지고 있다. 이와 같은 문제를 보완하기 위해 좀 더 일반적인 경우인 μ(s)가 미지인 경우 범용크리깅과 일반크리깅이 있으며, 범용크리깅에 대한 간략한 설명은 다음과 같다.

범용크리깅은 모형의 평균 (μ(s))는 미지라는 전제 하에 예측하는 방법으로 다음과 같은 관측된 지점의 값들의 선형결합으로 만들 수 있으며, 이라는 조건은 불편 성을 만족하는 조건이다.

본 연구에서는 기존의 일변량 크리깅 외에 하나의 측정지점에서 두 개 이상의 변수가 관측되는 경우 다변량 크리깅인 공동크리깅을 이용하여 수질오염 물질 중 BOD를 예측하였다. 하나의 측정 지점에서 관측되는 두 가지 이상의 공간변수들의 선형조합을 이용하여 미관측 지점을 예측하는 크리깅기법을 공동크리깅이라 하며, 예측하고자 하는 변수를 주변수라 하며, 주변수를 예측하기 위해 사용되는 보조변수로 구성되며, 보조변수는 하나가 아닌 여러 개가 될 수 있다.

공동크리깅에서 미관측 지점 s₀에 대한 예측값z^* (s₀)에 대한 일반식은 다음과 같다.

여기서, z는 주변수, n은 사용된 주변수의 총 자료수, n_s는 사용된 이차변수의 총개수, u_j는 j번째 이차변수, m_j는 j번째 보조변수의 총 자료수, λ는 가중치, 그리고 s 는 각자료의 위치이다. 따라서 주변수 값을 예측하기 위하여 총 (n_s+1)의 변수와(n+n_s × m_j)개의 자료가 공동크리깅에 사용된다.

공동크리깅은 주변수는 자료를 얻기 어렵고 상대적으로 보조변수가 자료를 얻기 쉬울 때 주로 사용되며, 주변수와 보조변수는 반드시 공간적 상호관계가 있어야 한다. 공동크리깅의 장점은 자료 획득이 어려운 주변수의 공간적 예측을 자료 획득이 상대적으로 쉬운 보조변수들을 사용함으로 예측의 정확도를 향상시키는데 있다(Wackernagel, 1998).

본 연구의 관심공간은 Fig. 1과 같이 금강유역을 대상으로 하였으며, 시간적으로 2011년 1월부터 2011년 12월까지 월별로 측정한 수질오염 물질 중 생분해성 유기물질 관리지표인 생물화학적 산소요구량(Biochemical Oxygen Demand, BOD)에 대한 공간적 예측을 위해 공간분석을 실시하였다. 수질자료는 국립환경과학원의 물환경정보시스템(MOE, 2012)을 이용하였다.

월별로 측정한 BOD 자료를 이용하여 베리오그램 분석을 수행하여 적합한 베리오그램을 선택한 후 두 종류의 크리깅 분석방법을 통해 공간예측을 실시하였으며, 예측된 값의 불확실성에 대한 분석도 실시하였다. 월별 BOD에 대한 공간적 예측을 위해 공동크리깅에서의 주변수로 BOD를 사용하였으며, BOD와 상관관계가 높은 화학적 산소요구량(Chemical Oxygen Demand, COD)을 보조변수로 하였으며, 월별 BOD와 COD의 상관계수는 Table 1과 같다.

관심공간에서 관측된 수질자료 중 BOD에 대한 세미베리오그램 모형으로 가우시안모형, 지수모형과 구형모형을 고려하여 이 중 오차가 가장 작은 모형을 선택하고자 한다. 본 연구에서는 모집단 분포에 대한 가정이 없고 자료가 대용량이라도 계산이 용이하다는 장점을 가지고 있는 가중최소제곱방법을 이용하였다.

가중최소제곱방법으로 세 가지 베리오그램을 적합하여 가장 적은 잔차제곱합을 제시하는 베리오그램을 선택하였으며 월별로 베리오그램을 적합한 결과는 Table 2에 나타내었다. 범용크리깅과 공동크리깅에 대해 대부분의 월에서 가장 작은 잔차제곱합을 가지는 모형은 가우시안 모형이며, 2월의 경우 지수모형이, 8월의 경우 구형모형이 동시에 선택되었다.

공간 예측에 대한 성능을 비교하기 위해 한 자료를 차례로 빼나가는 교차검증(leave-one-out cross validation)을 사용하였으며, 예측력 비교기준으로서 다음과 같이 정의되는 예측잔차제곱합(predicted residual sum of squares, PRESS) 값을 사용하였다.

여기서 Z(s_i)는 i번째 위치한 관측 값이며 (s_−i)는 i번째를 제외한 관측 값으로 i번째 위치한 값을 예측한 것으로 표현된다.

Table 3은 두 가지 크리깅 방법에 대한 예측력 비교를 나타낸 것으로 미관측 지점에 크리깅 예측 값이 얼마나 정확한지를 나타내주는 측도이다. 예측잔차제곱합이 낮을수록 정확한 예측을 했다고 할 수 있어 범용크리깅 방법보다 공동크리깅 방법이 상대적으로 더 정확한 공간 예측 값을 가진다고 할 수 있다.

베리오그램 분석을 통해 최적의 베리오그램 모형 선택과 모수 추정 후 두 가지 방법의 크리깅을 통해 공간예측에 대한 성능을 비교한 결과 범용크리깅에 비해 공동크리깅 방법이 공간예측력이 더 우수한 것으로 확인되었다. 따라서 미관측 지점의 예측 값을 지도상에 표출하기 위해 각 월별로 가중최소제곱추정방법으로 적합 시킨 후 최적의 베리오그램을 선택하여 공동크리깅을 실시하였으며, 격자 공간상에서 예측 값과 예측 값의 불확실성을 Fig. 3과 같이 표출하였다. Fig. 3에서 보는 바와 같이 각 월별 BOD의 공간예측 값이 공간적으로 비균질하다는 것을 알 수 있다. 이러한 공간적 비균질성은 수질오염원의 위치, 토지이용 패턴, 수리수문학적 특성, 하천 유형, 강우 특성 등 다양한 요인에 기인할 수 있다(Tu, 2011). 금강 유역 전체적으로는 상류부보다는 하류부에서 높은 BOD값의 공간적 분포를 보이며, 특히 4월 ~ 6월의 경우 금강 하구 근처에서 BOD값이 높게 분포함을 관찰할 수 있었다. Fig. 4는 예측 값의 불확실성을 나타내는 그림으로 예측된 값의 분산을 그림으로 나타내었으며, 각 월별로 예측 값에 대한 불확실성은 매우 비슷한 것으로 판단되었다.

본 연구의 목적은 미관측 지점을 포함한 관심 공간에서의 수질오염물질에 대한 공간적 예측을 위하여 다양한 크리깅 방법을 적용하여 예측성능을 비교, 평가함으로써 우수한 크리깅 방법을 제안하는데 있다. 본 연구에서는 다양한 크리깅 방법 중 범용크리깅과 공동크리깅 기법을 이용하여 수질오염물질 중 BOD값에 대해 관심공간 내 미관측 지점의 BOD의 공간 예측 값과 불확실성을 제공하였다. 본 연구에서는 공동크리깅의 경우 주변수인 BOD를 예측하기 위해 이차변수로 COD를 사용하였으며, 그 결과 범용크리깅보다 공동크리깅 방법이 예측 성능이 우수함을 교차검증을 통해 확인할 수 있었다. 따라서 수질자료 중 관심 대상인 주변수와 이차변수 간의 연관관계가 만족되면 공동크리깅의 방법이 다른 크리깅 방법보다 보다 정확한 예측 값을 제공하는 것을 알 수 있다.

본 연구에서는 2011년도 자료만을 통해 공간예측을 실시하여 해당 년도에 대한 BOD의 공간적 패턴을 확인할 수 있으나 시간적 흐름에 대한 시계열적 패턴을 고려하지 못한 한계를 가지고 있으나 향후 시간과 공간을 포함한 시공간모형(spatial-temporal model)을 이용하여 시간적 패턴과 공간적 패턴을 동시에 반영하여 예측하는 모형에 대한 개발이 필요하다. 또한 본 연구에서 사용한 공간분석은 관측지점 간 유클리디안 거리를 기반으로 공간모형을 설정하여 분석을 실시하여 하천의 흐름에 대한 거리를 이용하지 못하였다는 단점을 가지고 있으나 향후 연구에서는 관측지점 간 실제 물의 흐름에 따른 거리(water distance)의 계산을 통해 유클리디안 거리와의 결과와 비교할 예정이다.