• We proposed a novel approach for geoscience data mapping using the Gaussian Mixture Model (GMM).

• GMM with multivariate data was better for assessing the spatial distribution of soil contamination than OK and OCK.

• The use of multivariate data with GMM enhanced the understanding of spatial distribution.

정보통신기술(Information & Communications Technology; ICT) 및 센서의 발달로 방대한 양의 데이터 전송·수집이 가능해지면서, 지질자원 분야에서도 다량의 데이터 구축과 그 활용에 대한 관심이 증가하고 있다(Zhang et al., 2018; Bergen et al., 2019; Reichstein et al., 2019). 그러나 지질자원 분야에서 다루는 기작 및 지오데이터1)의 특수성으로 인해, 기계학습과 같은 대용량 자료 해석 기법이 제한적으로 적용되거나 혹은 강건하지 않은 추정 성능을 제공해왔다(Ryu, 2019).

예를 들면, 지질자원 분야의 관심변수(예: 대형재난, 광물자원 매장량 등)는 시공간적으로 드물게 발생할 뿐만 아니라 지하 수 미터에서 수 킬로미터에 다다르는 심부공간에서 발생한다. 다양한 제약조건(예: 막대한 시추 비용 등)에 따른 하드 데이터(hard data)의 희소성을 극복하기 위해, 비저항탐사 등의 간접 관측 결과나 주관적 사후 평가가 반영된 소프트 데이터(soft data; 예: 토지피복 등)가 타 분야에 비해 빈번히 보조변수로 활용된다. 또한, 지오데이터는 공간적 이질성 및 불확실성이 높은 특징을 지니는데, 이는 지질 매체에 내재된 불균질성, 한정된 시료 개수에 따른 통계적 결함, 혹은 공간적 이질성에 따른 변동성 차이 등에 의해 유발된다.

상기 언급한 지오데이터의 고유 특징으로 인해 가중선형조합(weighted linear combination)을 이용하는 지구통계 접근방식은 지오데이터를 해석하는데 유효하지 않을 수 있으며, 자연계 현상(실제 모집단)을 반영하지 못하는 불량한 결과를 제공할 수 있다. 또한 다변량의 지오데이터는 기존의 지리통계적 접근법(geostatistical approach)으로 재현하기 어려운 다변량 복잡성(multivariate complexities)을 내재하여, 이를 극복할 수 있는 새로운 통계적 기법을 요구한다(Barnett and Deutsch, 2015; Silva and Deutsch, 2018).

이 연구는 공간분포를 추정하는 전통적인 선형기반의 지구통계기법에서 탈피하여 다변량 및 공간 이질성을 지닌 지오데이터의 특성을 반영하기 위해 다중분포모사에 효과적인 가우시안 혼합 모델(Gaussian Mixture Model, GMM)을 적용하여 강건한 공간 예측 결과를 도출하고자 하였다. 지질자원 분야에서 GMM은 데이터 군집화(Kim et al., 2014; Herms et al., 2021), 분류(Grana et al., 2017; Kim et al., 2019), 확률밀도추정(Chen et al., 2006)에 활용된 바 있다. 그러나 데이터의 공간분포 추정울 위해 GMM을 적용한 사례는 찾아보기 어렵다. 최근 Qu and Deutsch (2018)은 정상성(stationarity) 가정을 위배한 원시 데이터를 대상으로 GMM을 통한 변환과정을 거쳐 자료의 비정상성(non-stationarity)을 제거하였으며, 이를 통해 지구통계 시뮬레이션의 공간추정 성능 개선을 보고한 바 있다.

본 연구에서는 다변량의 확률분포를 추정하고 이를 반영한 공간예측결과가 기존 지구통계기법 중 단변량을 활용한 정규 크리깅(Ordinary Kriging; OK)과 이변량을 활용한 정규 공동크리깅(Ordinary Co-Kriging; OCK))에 비해 정확성 및 정밀성을 얼마나 향상시키는 지를 평가하였다. 성능 평가를 위해 충남 서천군 옛 장항제련소 인근 토양오염 정밀조사 부지를 테스트 베드로 선정하였다. Kim et al., (2019)은 동일지역을 대상으로 휴대용 X선형광분석기(Portable X-Ray Fluorescence Spectroscopy; PXRF)를 보조변수로 활용하는 OCK 기법을 제안한 바 있다. OCK는 보조자료를 공간추정과정에 활용하는 확장된 크리깅 방법으로 다양한 지구과학 분야에서 활용되었으며(e.g., Goovaerts, 2000), 이외에도 다변량 변수를 고려하기 위해 다중변수를 단일 병합 보조변수(single merged secondary dataset)로 변환한 후 OCK를 활용한 사례(Badak and Deutsch, 2009)나 다중가우시안 크리깅(multi Gaussian kriging)과 다변량회귀분석(multivariate regression)을 활용한 연구도 있다(Goovaerts et al., 2005). 본 연구는 다변량 및 다봉분포(multimodal distribution; 서로 다른 두 개 이상의 최빈값을 갖는 연속확률분포)의 구조를 고려한 GMM을 공간추정모델링을 활용한다는 측면에서 기존 연구들과는 차별점이 있다.

2.1. 연구지역 개황

다변량 변수와 GMM을 이용한 공간분포 추정 방법의 성능을 평가하기 위해 충남 서천군 옛 장항제련소 인근 토양오염부지 정밀조사 결과를 활용하였다. 장항제련소(1936-1989년)는 운영 과정에서 발생한 오염물질로 인해 토양오염 및 농작물 피해 등 다양한 환경문제를 유발하였다.

Moon et al. (2011) 및 Kim et al. (2016)에 따르면 장항제련소 인근 토양의 중금속 오염은 원광석에 포함된 비소, 납, 구리 등이 제련과정에서 굴뚝을 통하여 비산되거나 원광석의 운반과정 중 분진이 비산되어 유발된 것으로 판단된다. 환경부는 비소 및 중금속 오염토양이 제련소 굴뚝을 중심으로 반경 2km까지, 수직적으로 지표하 최대 1m까지 분포되어 있는 것으로 보고하였다(KMOE, 2009a). 또한, 제련소 굴뚝으로부터 거리가 멀어질수록 오염물질의 종류 및 농도가 점차 감소되는 경향을 보였다. 이는 굴뚝으로부터 비산되는 오염물질에 의한 표층오염과 낙하거리가 상관성이 있음을 보여준다.

2.2. 환경부 정밀조사지침 및 시료개수

토양정밀조사지침(KMOE, 2009b)에 따르면 600~1,000m²당 1개 지점을 조사해야 하며, 따라서 연구부지 면적(3,300m²)에서는 최대 6개 지점을 조사하면 되었다. 그러나 본 연구에서는 정밀한 토양오염 면적 산정과 GMM 기법의 신뢰성을 검토하기 위해 정밀조사지침 기준보다 많은 156개 지점에서 토양시료를 조사하였다. 기준필지를 대상으로 정방형의 5m 간격으로 시료를 채취하였다. 현장시료채취 및 분석은 서울대학교 NICEM(National Instrumentation Center for Environmental Management)에서 수행하였다.

2.3. 토양분석방법

2.3.1. 실험실 분석

실험실 분석용 시료(ex-situ; n=153개)는 5-10개의 토양시료를 채취한 후 혼합하여 복합시료 형태로 채취하였다. 채취한 시료는 토양분석실로 이송한 후, PE 용기에 균일한 두께로 펼친 후 약 10일간 풍건시켰다. 건조된 토양을 10-mesh 체(<2 mm)나 막자사발에 넣어 분쇄한 후 100-mesh 체(<125 μm)를 이용하여 체질 후 토양 중 비소, 납, 구리, 니켈, 아연 용출 시험에 이용하였다. ISO 11466에 따라 왕수에 의한 용출실험을 수행하고(ISO, 1995), 유도결합플라즈마-원자발광분광법(Inductively Coupled Plasma Atomic Emission Spectroscopy, ICP-AES)으로 토양 내 중금속(납, 구리, 니켈, 아연) 및 비소 농도를 측정하였다.

토양오염공정시험기준법에서는 ICP-AES를 사용하여 토양에 함유된 화학물질의 농도를 분석하도록 정하고 있다. ICP-AES는 여러 개의 화학원소를 동시에 분석가능하며 검량선의 직선영역이 넓은 장점을 갖지만, 추출 과정에서 다량의 산을 사용하여야 하고 지각 물질(예: Al이나 Fe 등)의 간섭을 받을 수 있다(Lee et al., 2015). 또한, 복잡한 토양 매질에 의한 영향 및 분광학적 간섭으로 인해 측정하는 파장에 따라 실제 농도보다 과대 혹은 과소평가될 우려가 있다.

2.3.2. 현장 분석

PXRF는 고체물질에 X선을 흡수(혹은 투과), 회절(혹은 산란)시키고, 다른 색의 X선을 만들어 주사하여 화학물질의 전 함량을 분석한다. 예를 들면, 특정원소에 X선을 주사하면 K 준위에 있던 전자들은 X선으로부터 에너지를 받아 들뜬 상태에 이르며, 들뜬 상태의 원자가 원래의 에너지 준위로 돌아올 때 특정 에너지 파장을 방출한다. 비소의 경우, L에서 K로 갈 때 10.53 KeV를, M에서 K로 갈 때 11.73 KeV를 방출한다. 이 방출파장은 각 원소마다 특정 값을 지니기 때문에 각 peak들의 에너지 값과 면적을 통해 원소의 종류와 함량을 평가할 수 있다.

그러나 현장 PXRF 분석은 토양 특성에 따른 여러 가지 간섭 효과에 의해 정밀도가 떨어진다는 문제가 있다. 예를 들어 토양 사이의 공극이나 물질들의 물리적 배열상태(matrix effect)에 따라 오차가 발생하며, 토양 내 수분함량이 높은 경우에도 간섭이 발생할 수 있어, 이를 극복하기 위해 많은 연구가 진행 중에 있다(Ryan et al., 2017; Lemiere, 2018).

본 연구에서 PXRF을 통한 현장분석(in-situ; n=156개)은 US EPA (2007)의 PXRF 비파괴분석법을 준용하였다. 비파괴 분석방법의 일종인 PXRF를 이용한 측정법은 미국 Resource Conservation and Recovery Act의 SW486에 포함된 EPA method 6200에 명시되어 토양, 하천 침전물 내 중금속 분석 방법으로 활용된다.

EPA method 6200은 분석방법을 크게 두 가지로 구분하고 있다: in-situ 현장 스크리닝은 짧은 시간 안에 여러 지점에 걸쳐 많은 시료를 측정하며, ex-situ 조사는 토양시료를 전처리하여 분석하는 ex-situ 정밀조사 방법이다. PXRF를 이용한 ex-situ 정밀분석은 토양오염공정시험법과 비슷하게 토양시료를 건조 및 125㎛ 미만으로 체질하여 토양수분이나 matrix effect에 의한 간섭효과를 최소화한다.

본 연구는 in-situ 현장 스크리닝 분석방법을 활용하여 정확도는 다소 낮지만 경제적이고 편리한 방법으로 토양화학자료를 구축하였다. 현장에서 자갈 등 조립질 입자와 유기물(예: 낙엽, 잔가지 등)을 제거한 후 PXRF를 통해 토양 중 비소, 납, 구리, 니켈, 아연 농도를 측정하였다.

연구지역에서 비소와 납의 공간분포를 파악하기 위해 ICP-AES로 측정된 비소와 납을 주변수로, ICP-AES로 측정된 그 외 물질(구리, 아연, 니켈)과 PXRF로 측정된 모든 항목을 보조변수로 활용하였다. 지구통계적 기법에 대한 설명은 Kim et al. (2019)에 자세히 기재하였으며, Fig. 1을 통해 통계절차 흐름도를 도식화하였다.

Figure 1. Flow chart of statistical procedures for estimating soil contamination areas.

3.1. 중요 변수 선택(variable selection) 기법

다변량의 보조변수들을 이용하여 주변수인 납과 비소의 공간분포를 산정하기에 앞서, 예측 모형의 과적합(overfitting)을 방지하고, 모델 복잡도를 최소화하는 한편, 예측 속도를 향상시키기 위해 변수 선택 기법을 적용하였다. 변수 선택을 위해 랜덤포레스트 기반의 Boruta를 활용하였으며(Fig. 1 Step 1-②), Boruta는 다음과 같은 순서로 변수를 선정한다(Kursa and Rudnicki, 2010).

① 모든 변수의 복사본 생성(adding copies of all variables)

② 그림자 변수(shadow attributes)²⁾ 생성: 변수 복사본을 무작위로 섞어서 변수 간 상관성을 없앰. 이들은 무작위성을 통해 만들어진 새로운 변수임

③ 그림자 변수를 포함한 표본들을 대상으로 랜덤포레스트 모형을 구동하고, 각 변수에 관한 정확도 손실(accuracy loss)을 표준화 점수(Z-score)로 계산

④ 그림자 변수들 중에서 최대 표준화 점수를 갖는 변수(maximum Z score among shadow attributes; MZSA)를 확인

⑤ MZSA를 기준으로, 입력 변수의 표준화 점수가 높으면 중요(important), 낮으면 중요하지 않다(unimportant)고 간주

⑥ 모든 변수의 중요도(importance)3)가 정해질 때까지 ②-⑤ 과정을 반복

3.2. GMM을 통한 확률분포추정

3.2.1.기존 비모수 밀도추정을 통한 확률분포 추정법의 한계 및 대안

일반적으로 확률분포의 추정은 모집단(population)의 모수(parameter; 예: 가우시안 분포의 경우 1차 모멘트인 표본평균 μ_i와 2차 모멘트인 표본분산 σ_i²)를 최대우도추정(Maximum Likelihood Estimation, MLE)에 기초하여 추정한다. 이때 확률밀도함수는 단봉분포(univariate distribution) 형태를 가정한다. 그러나 실제 데이터(특히 지오데이터)는 대부분 복잡한 다변량 분포특성을 지니며, 다봉분포(multimodal distribution) 형태인 경우가 대부분이다(Silverman, 1998).

상기 문제를 극복하기 위해 지오데이터의 확률밀도함수 추정 시 비모수적(non-parametric) 밀도 추정법을 활용해왔다. 그러나 비모수적 방법(예: Kernel density estimation)은 확률분포 계산에 드는 시간과 비용이 크게 증가하는 문제를 유발하며(Qu and Deutsch, 2017; Silva and Deutsch, 2018; Silverman, 1998), 또한 지질자원 분야와 같이 방대한 보조변수를 활용해야 하는 경우에는 비모수 밀도추정기법의 적용이 어렵다.

GMM과 같은 반모수적(semi-parametric) 방법이 대안으로 사용될 수 있다(Scrucca et al., 2016). GMM을 통한 확률분포추정은 기존의 비모수적 기법과 동일하게 밀도함수의 분포 형태를 가정하고 모수를 추정(예: 가우시안의 경우 평균과 분산)하는데, 해당 밀도함수가 복수(다봉)의 확률분포로 구성된다는 특징을 갖는다.

3.2.2. 가우시안 혼합모형(GMM)을 통한 복수의 확률밀도함수 추정

GMM은 개별밀도함수를 전체 확률밀도함수의 성분(component)인 커널(kernel)로 간주한다. 다중 가우시안 분포를 이용할 수 있으므로 여러 개의 중심점을 갖는 다차원 데이터에 대해 견고한 모델링이 가능하다.

GMM을 통한 확률밀도함수 추정은 복수(m)의 가우시안 확률밀도함수에 대한 선형 결합으로 표현되며, 이를 식으로 표현하면 아래와 같다(Reynolds, 2009):

(1)$$p\left(x|\theta \right)=\sum _{i=1}^{m}p\left(x|{\omega }_i,{\theta }_i\right){\alpha }_i$$

위 식에서 $p\left(x|{\omega }_i,{\theta }_i\right)$는 i 번째 가우시안 확률밀도함수를 의미한다. αi는 혼합가중치로 각 확률밀도함수(즉, m 개 성분)의 상대적인 중요도를 의미하며, 아래와 같은 제약조건을 따른다(Han, 2009):

(2)$$0\le a\le 1and\sum _{i=1}^m\alpha =1$$

가우시안 확률밀도함수의 모수(parameter; θ)는 평균 μ_i, 분산 σ_i², 그리고 가중치 αi로 구성된다. 전체 모델을 이루는 가우시안 분포들은 대각(diagonal), 구형(spherical), 혹은 타원체(ellipsoidal) 형태의 공분산 행렬을 갖는다. GMM은 EM 알고리즘(Expectation-Maximization algorithm; 기댓값 최대화 알고리즘)을 활용하여 데이터의 로그우도(p[x|w])를 최대로 하는 가우시안 성분들의 모수, 즉 식(1)의 평균 μ_i, 분산 θ_i, 가중치 α_i를 추정한다. EM 알고리즘은 모집단의 모수 추정값(예: 최대우도 혹은 최대사후확률)이 최적해에 수렴할 때까지 계산을 반복하는 방법이다(Reynolds, 2009).

GMM에는 다중 확률밀도함수를 추정하는 다양한 모델을 제공한다. 각 성분 가우시안 분포는 평균벡터 μ를 중심으로 하는 타원체이며, 서로 다른 기하학적 특징(예: 부피, 형상 및 방향)에 의한 공분산 행렬로 결정된다.

이 연구에서는 혼합 가우시안 분포의 성분 개수(m)와 mclust에 탑재된 14가지 기하학적 조합을 가진 다변량 혼합 모형(multivariate mixture)을 고려하였다. 각 혼합모형의 MLE가 최대가 되도록 EM 알고리즘을 수행하였다. EM 알고리즘을 통해 도출된 각 혼합 가우시안 분포의 성분 개수(m) 및 기하학적 모형들의 조합 중에서 BIC(Bayes Information Criterion)가 가장 큰 모형을 최종 가우시안 혼합모형으로 선정하였다(Fig. 1 step 2-②).

3.3. 통계적용절차

ICP-AES로 측정된 비소와 납 데이터 중 123개(약 80%) 시료를 6개의 하위 군집으로 재구성(30, 49, 61, 76, 91, 107개)하여 GMM 학습에 활용하였다. 다변량자료와 GMM을 적용하여 토양화학자료의 공간분포를 산정한 결과의 성능을 비교하기 위하여 단변량자료를 이용한 OK 및 이변량자료를 이용한 OCK도 수행하였다. OK 및 OCK를 다변량 공간 추정기법의 대조군으로 활용하여, 사용변수 증가에 따른 모델 예측력 향상 정도를 함께 고찰하고자한다. 총 10개의 화학분석결과 중 각 기법별 사용변수는 아래와 같다:

○ 단변량자료를 이용한 OK: ICP-AES의 납 또는 비소농도

○ 이변량자료를 이용한 OCK: 주변수로 ICP-AES의 납 또는 비소 농도를 사용, 보조변수로는 PXRF의 납 또는 비소 결과를 사용

○ 다변량자료를 이용한 GMM: ICP-AES의 납 혹은 비소를 주변수로 사용, 나머지 in-situ와 ex-situ 8개 변수들(ICP-AES로 측정한 구리, 니켈, 아연 농도 및 PXRF 분석 결과)을 보조변수로 활용

ICP-AES로 측정된 비소와 납 데이터 중 나머지 시료(30개; 약 20%)는 모델 검증에 활용하였다. ICP-AES 분석값과 모델 예측값 간 차이의 평균제곱근오차(Root Mean Square Error, RMSE) 및 피어슨 상관계수(Pearson correlation coefficient, r)를 평가하고 이를 성능평가지표로 활용하였다. RMSE가 작을수록, 상관계수(r)가 클수록 모델 예측력이 개선되었음을 지시한다. 이외에도 모든 ICP-AES 분석 결과(153개)를 활용하여 OK로 추정한 공간 분포를 대조군(control group)으로 활용하였다.

OK 및 OCK는 비대칭도(strong skewness)를 보이는 비정규(non-Gaussian)분포에 대한 공간 예측에 한계를 갖는다. 정규성 확보를 위해 토양화학 10개 변수에 로그변환을 수행하였다.

GMM(mclust v5.4.10)(Scrucca et al., 2016) 및 Boruta(Boruta v7.0.0)(Kursa and Rudnicki, 2010)을 포함한 통계패키지는 R-3.6.1(R Development Core Team)을 통해 수행되었으며, OK 및 OCK 등 지구통계기법은 SGeMS(Remy et al., 2009)를 통해 수행되었다.

4.1. 토양시료의 중(준)금속 오염특성

본 연구에서 조사된 비소 및 중금속의 농도는 Table 1과 같다. ICP-AES를 이용하여 분석된 각 원소별 평균농도는 74.48(비소), 196.93(납), 56.02(구리), 11.45(니켈), 43.25(아연) mg/kg이다. 각 원소별 평균농도와 국내 토양오염 우려/대책기준(25/75(비소), 200/600(납), 150/450(구리), 100/300(니켈), 300/900(아연) mg/kg)을 비교하였을 때, 비소의 경우 우려기준을 상회하며, 대책기준에 준하였다. 납의 경우에도 평균농도가 우려농도에 준하나, 나머지 원소들의 경우는 우려 및 대책기준에 비해 낮은 농도특성을 보인다. 우려기준을 대상으로 시료별 초과 개수(초과율)는 비소, 납, 구리, 니켈, 아연이 각각 122개(80%), 63개(41%), 1개(1%,), 0개(0%), 0개(0%)이며, 대책기준을 대상으로 살펴보면, 67개(44%, 비소), 2개(1%, 납), 0개(0%, 구리), 0개(0%, 니켈), 0개(0%, 아연)가 초과하였다. 따라서 연구지역의 토양은 특히 비소와 납의 오염이 심각한 것으로 드러났다. 반면 토양 중 구리, 니켈, 아연 등의 미량 중금속 원소의 경우 기준치 이하로 나타난바, 제련 활동으로 인한 영향은 미비한 것으로 판단할 수 있다.

Table 1 . The descriptive statistics of metal (loid) contents in soil samples by ex-situ (ICP-AES) and in-situ (Portable XRF) measurements.

Unit: mg kg-1	Laboratory analysis using ICP-AES (n=153)					Portable XRF (PXRF) measurements in the field (n=156)
	As	Pb	Cu	Ni	Zn	As	Pb	Cu	Ni	Zn
Minimum	4.41	13.71	6.21	3.87	16.00	0.50	12.00	12.00	14.60	137.00
Maximum	236.70	961.67	167.70	33.70	112.83	143.00	430.00	145.00	25.00	205.00
Range	232.29	947.96	161.49	29.83	96.83	142.50	418.00	133.00	10.40	68.00
Median	66.33	160.47	53.23	11.47	43.47	22.00	60.50	34.00	20.30	163.00
Mean	74.48	196.93	56.02	11.45	43.25	27.95	74.33	38.05	20.30	164.40
SE.mean*	4.17	13.42	2.67	0.30	1.14	2.04	4.77	1.46	0.17	0.79
CI.mean**	8.24	26.51	5.27	0.60	2.26	4.04	9.41	2.89	0.33	1.57
Var.***	2662.33	27540.86	1087.04	13.89	200.45	651.63	3542.20	333.70	4.49	98.59
Std.dev.****	51.60	165.95	32.97	3.73	14.16	25.53	59.52	18.27	2.12	9.93
Coef.var*****	0.69	0.84	0.59	0.33	0.33	0.91	0.80	0.48	0.10	0.06

* SE.mean: the standard error of the mean; ** CI.mean: the confidence interval of the mean at the p level of 0.95; *** Var: the variance; **** Std.dev: the standard deviation; ***** Coef.var: the variation coefficient defined as the standard deviation divided by the mean.

표준편차와 분산의 경우에도 평균대비 높은 수준의 변동계수를 보인다. ICP-AES 및 PXRF 분석 결과의 변동계수(coefficient of variation; CV)를 살펴보면, ICP-AES의 변동계수는 59(최소; 아연)~84(최대; 납과 구리)%를 보인다. 변동계수의 크기 차이는 곧 공간적 이질성의 차이가 큼을 지시한다. PXRF의 변동계수는 6(최소; 아연)-48(최대; 구리)%로 ICP-AES 분석 결과에 비해 낮은 변동계수를 보인다.

Kim et al. (2016)은 SM&T(Standards, Measurements, and Testing program of the European Commission) 방법을 통해 이온교환(F1), 산화물에 결합(F2), 유기물 및 황화물에 결합(F3), 잔류상 형태(F4)로 구분하여 비소 및 중금속의 농도를 측정하였다. 4단계 연속추출을 수행한 결과, 본 연구지역에 가깝게 위치한(제련소 반경 1 km 이내) 토양의 비소가 주로 산화물과 결합한 형태(F2; 60%)인 것으로 확인되었으며, 고농도 비소와 납의 원인은 제련활동에 의한 인위적 영향으로 추정하였다. Lee et al. (2019)는 제련소에서 4-7 km 반경 북동쪽에 위치한 토양 및 암석 시료를 대상으로 지화학 및 납 동위원소를 분석하고, 광화작용(자연기원)의 의한 비소 오염을 확인한 바있다. 장항제련소 내 토양 오염은 오염원과의 거리 및 토양광물 성상에 따라 그 기원의 차이가 있음을 시사한다.

4.2. 토양환경변수 간 상관성 고찰

각 변수별 상관성은 Fig. 2와 같다. 도표의 대각선 방향은 각 원소의 분포도(히스토그램)와 적합(fitting) 결과(붉은색 선)를 함께 도시한 것이다. 대각선 상단은 변수별 상관계수를 표시하였으며, 숫자 위의 붉은 별의 개수는 신뢰수준을 의미한다. 대각선 하단은 변수별 이변량도(bivariate scatter plots)와 적합 결과를 함께 도시한 결과이다. 일부 항목에서(예: ICP-AES로 분석된 구리와 니켈과의 상관성 및 적합 결과) 이변량 분포 간 비선형성 및 이분산적(heteroscedastic) 특징이 관찰되며, 따라서 본 연구 데이터는 다변량 확률분포 모델링과 같은 접근방식이 필요하다.

Figure 2. Correlation chart of multivariate soil data. The histogram of each variable is shown on the diagonal with primary variables (As and Pb determined using ICP-AES) in blue and auxiliary variables (Cu, Ni, and Zn using ICP-AES and all PXRF data) in green. The values on the upper right of the diagonal indicate the correlation coefficients. The stars indicate the significance levels (*** p < 0.001, ** p < 0.01, * p < 0.05). On the bottom left of the diagonal, the bivariate scatter plots are displayed with a fitted red line. Values on the boundaries indicate the concentrations on a log scale.

ICP-AES 분석 결과 간의 상관성을 비교해보면, 토양 내 비소와 납(0.91), 비소와 구리(0.92), 납과 구리(0.96)에서 유의한 신뢰성과 높은 상관성을 보인다. 앞서 기초통계의 평균 및 환경부 기준 등을 고려하였을 때, 제련소에 의한 토양 내 비소와 납 오염은 파악하였으나 구리는 대부분 시료에서 우려기준을 초과하지 않았기 때문에 구리 농도에 대한 제련소의 영향이 없는 것으로 판단하였다. 그러나 상관성을 고려할 시 비소와 납뿐만 아니라 구리도 제련활동의 영향을 받는 것으로 간주될 수 있다. 실제 연구 지역의 구리 농도는 국내 토양의 배경 농도(15.3 mg/kg; KMOE, 2011; Yang et al., 2014)과 비교했을 때 높은 수준이다.

구리는 자연계에서 자연동으로 산출되며, 주요광석은 황동석(CuFeS₂), 휘동석(Cu₂S), 공작석(Cu₂CO₃(OH)), 적동석(Cu₂O) 등이 있다. 지각 내 함유량은 55~65 mg/kg으로 추정되며, 화성암은 55 mg/kg, 퇴적암은 5~45 mg/kg, 토양은 20 ppm로 추정하고 있다(Koljonen, 1992; Galán et al., 2008). 인위적인 구리의 주요 발생원은 (1) 제련공정, (2) 가공공정(합금 등), (3) 화합물제조공정을 들 수 있다. 상기의 정보들과 기초통계 및 상관성 등을 종합해보면 연구 지역 토양 내 구리의 농도 분포는 단순히 자연계에 분포하는 지각이나 토양에서 기원한 것이 아니라 연구지역의 주요한 인위적 인자인 제련활동의 영향을 받은 것으로 보인다. 정확한 원인 파악을 위해서는 중금속 동위원소분석 등 환경정밀분석이 필요하다.

4.3. PXRF와 ICP-AES 분석 결과 비교

비소와 납 농도의 ICP-AES와 PXRF 분석 결과 간의 상관계수는 0.8(p=0.00)로 서로 유의한 상관성을 보인다. 구리도 ICP-AES와 PXRF 분석 농도 간의 상관계수는 0.6으로 상관성을 보인다. 반면, 니켈과 아연의 ICP-AES와 PXRF 분석 결과 간의 상관계수는 각각 0.1과 0.38로 선형적으로 유의하지 않은 관계를 나타내고 있다.

결정계수가 유의미한 비소, 납, 구리의 농도를 대상으로 관계식을 도출해보면 다음과 같다:

ICP-AES와 PXRF의 비소 분석 결과: ICP As = 1.72 × PXRF As + 29.66

ICP-AES와 PXRF의 납 분석 결과: ICP Pb = 1.89 × PXRF Pb + 58.34

ICP-AES와 PXRF의 구리 분석 결과: ICP Cu = 1.02 × PXRF Cu + 17.60

특히 비소와 납의 관계식을 보면 PXRF 분석 결과 대비 ICP-AES 분석 결과가 1.7-1.9배까지 차이를 보여, PXRF 분석 결과의 정확도가 다소 결여되어 있음을 알 수 있다. 위와 같은 차이는 (1) PXRF의 현장 측정에 따른 간섭효과와 (2) 제련 광종의 생성 기원에 따른 차이가 함께 영향을 준 것으로 판단된다. 그러나 현장에서 전처리를 거치지 않고 바로 분석한 PXRF 분석 결과가 전함량법에 의한 ICP-AES 분석 결과와 강한 선형 상관성을 보여(비소와 납 모두 r=0.8로 높은 수준의 정밀도를 보임) PXRF를 이용하여 분석한 비소, 납, 구리 농도는 각각 비소, 납, 구리의 공간분포 예측을 위한 보조변수로 충분히 활용할 수 있을 것으로 보인다.

4.4. 변수 선택법 적용 결과

랜덤포레스트 기반의 Boruta를 이용하여 보조변수들이 비소와 납 농도에 미치는 중요도를 계산하였다(Fig. 1; Step 1-②). 먼저 비소 예측을 위한 변수 중요도 계산 결과, 비소를 제외한 9개 설명변수 중에서 PXRF 니켈과 아연은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었다. 확정된 보조변수들 중에서 ICP-AES 납과 구리가 상위 영향 인자(상대적으로 높은 중요도)로 파악된다. 이는 앞선 상관성 분석 결과와도 일치하며, 제련과정에서 발생한 비소, 납, 구리의 영향이 함께 반영된 것으로 판단할 수 있다.

납에서도 마찬가지로 상위 중요도로 분류된 변수는 ICP-AES 비소와 구리이다. PXRF 니켈은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었으며, PXRF 아연은 최대 그림자 변수(shadowMax)와 유사한 중요도를 보여 잠정(tentative) 변수로 분류되었으나, 중요도 이력 정보 및 이상치 분포 등을 종합적으로 고려하여 납의 공간 예측 시 PXRF 아연은 활용하지 않았다.

따라서 GMM을 이용한 비소 농도의 공간 분포에는 ICP-AES 납, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하고, 납 농도의 공간 분포 추정에는 ICP-AES 비소, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하였다.

4.5. 다변량 토양화학자료기반 공간분포 예측 결과

4.5.1. GMM을 이용한 다변량 통계분포 추정

ICP-AES의 납 또는 비소를 주변수로 나머지 변수(PXRF 아연과 니켈 제외)들을 보조변수로 활용하여 GMM을 수행하였다. m개 성분의 가우시안 분포의 모수인 μ_i, θ_i과 성분별 가중치 α_i를 구하고 개별 데이터(시료)들이 어느 성분(m)에 속하는지를 평가하기 위해 모든 기하학적 모델의 조합을 대상으로 EM 알고리즘을 수행하였다. 각 모델별 BIC를 비교한 결과, 가우시안 확률분포의 성분 개수(m)는 2개에서, 14가지 기하학적 모델 중에서는 EVE(가우시안 분포의 기하학적 특징이 ellipsoidal, equal volume and orientation인 모형)가 가장 큰 BIC 값을 제공했다. 이를 통해 본 연구의 토양화학분석결과는 2개의 가우시안 혼합모형(쌍봉분포; bimodal distribution)이 모집단의 다변량 확률분포 예측에 최적인 것으로 판단할 수 있었다.

Fig. 3은 EVE 모형 및 2개 성분의 가우시안 혼합분포를 적합시킨 결과이며, 직관적 이해를 위해 다차원 공간(8-d)의 다중(쌍봉)확률밀도함수를 2차원 단면으로 잘라서 이변량도에 도시하였다. 즉, 8차원(주변수 및 보조변수 7개)의 다변량 분포를 일반적인 정규분포의 형태가 아닌 쌍봉 형태의 확률밀도함수로 표현한 것이다. GMM을 통한 다변량 확률밀도 함수 추정 결과, 선형관계(예: Fig. 3 내 ICP-AES를 이용하여 분석한 비소와 납 농도) 뿐만 아니라 비선형적 관계(예: Fig. 3 내 ICP-AES를 이용하여 분석한 납과 구리 농도)에 대한 확률분포 추정이 가능함을 알 수 있다. 특히 쌍봉의 가우시안 확률밀도함수의 적합을 통해, ICP-AES의 구리, 니켈, 아연의 비선형적 특징을 반영한 통계분포 모델링이 수행되었음을 확인할 수 있다. 한편, Fig. 3에 ICP-AES 비소와 납의 2차원 단면은 단봉분포 및 선형적 관계로 보이지만, 이는 자료의 내재적 특징에 기인함을 확인할 수 있다(Fig. 2의 비소 납의 선형관계 참조). 상기 결과를 종합해보면, GMM은 이질성을 내포한 데이터 모집단의 통계 구조를 모사·반영하고 확률밀도 추정 시 강건한 모델링을 제공할 수 있음을 지시한다.

Figure 3. Pairwise probability density functions fitted to the Gaussian Mixture Model (GMM) and EVE (ellipsoidal, equal volume and orientation) model. Values on the boundaries indicate probability density.

4.5.2. 기법별 공간분포 예측결과

기법별 공간분포 예측결과 중 비소의 공간분포 추정 결과를 Fig. 4에 도시하였다. 주변수를 모두 사용한 단변량의 OK 결과(Fig. 4d; ICP-AES 비소 153개 사용)를 49개의 비소 농도만을 이용한 OK(Fig. 4a), PXRF 비소 농도(n=156)를 보조변수로 활용한 이변량의 OCK(Fig. 4b), 그리고 변수 선택법을 통해 선정된 다변량의 보조변수(ICPAES를 이용하여 분석한 납, 구리, 니켈, 아연 농도(n=153) 및 PXRF의 비소, 납, 구리 분석 결과(n=156))를 활용한 GMM 공간분포(Fig. 4c)와 비교하였다.

Figure 4. Spatial distributions of As concentrations in the study area: (a) the results of ordinary kriging (OK) using 49 training data; (b) and (c) the results of ordinary co-kriging (OCK) and Gaussian Mixture Model (GMM), respectively, using 49 training data and 156 auxiliary data (As determined by PXRF for OCK and ICP-AES Pb, Cu, Ni, and Zn and PXRF As, Pb, and Cu for GMM); (d) the most realistic distribution of contamination that was estimated by OK using all the ICP-AES data (n=153). The black line indicates the regulatory level (25 mg/kg) by the Soil Environmental Conservation Act of the Republic of Korea.

OK, OCK, GMM 분석결과 모두 서쪽에 분포한 고농도 비소 오염지역(핫스폿; hot spot)을 효과적으로 추정하였다. 그러나 OK는 서쪽 중앙부에 위치한 저농도 지점(콜드스폿; cold spot)을 예측하는 데 실패하였다. 이는 콜드스폿에 대한 주변수(ICP-AES의 As 49개) 정보가 빠져있기 때문이다.

반면, OCK나 GMM은 주요 변수(ICE-AES로 분석한 As 농도)가 저농도 지점에 대한 정보를 갖지 않음에도 불구하고 보조변수의 공간적 상관성을 이용하여 원래 자료의 콜드스폿을 효율적으로 찾을 수 있었다. 특히 GMM은 비소 대책기준 미만 농도에 준하는 경계(콜드스폿)를 예측하였으며 대조군 자료와도 가장 유사한 결과를 보였다(Fig. 4c-d).

환경재해예측이나 자원탐사 분야에서는 핫스폿 혹은 콜드스폿 등 공간이상치(outlier)를 파악하는 것이 중요하다. 위와 같이 원 변수가 중요지점에 대한 정보를 모르는 경우에도 다변량 변수 간의 공간적인 관계성을 토대로 한 OCK나 GMM의 기법들이 원 변수 공간분포 추정에 매우 효과적일 수 있음을 시사한다.

4.5.3. 기법별 예측결과 통계량 비교

단변량(OK), 이변량(OCK) 및 다변량(GMM)을 통한 기법들의 예측성능을 RMSE와 상관계수(r)를 활용하여 비교 평가하였다(Table 2; Fig. 5&6). 세 기법 모두 샘플링 밀도 증가에 따라 전반적으로 RMSE는 낮아지며 상관계수 값은 높아지는 경향을 보였다(Fig. 5&6). 이는 샘플링 밀도 증가에 따라 예측의 정확도/정밀도가 향상된다는 것을 의미한다.

Table 2 . Correlation coefficient (r) and root mean-squared error (RMSE) between the measured and predicted values using the validation data (n=30 determined by ICE-AES) at different sampling densities for training (n=30 to 107 determined by ICP-AES) by each model: OK (ordinary kriging), OCK (ordinary co-kriging), GMM (Gaussian mixture model).

As						Pb
Sampling density	Prediction method	r	RMSE	rGMM ‒ rgeost.*	RMSEGMM‒RMSEgeost.**	Sampling density	Prediction method	r	RMSE	rGMM ‒ rgeost.*	RMSEGMM‒RMSEgeost.**
30	OK	0.59	0.28	0.27	-0.05	30	OK	0.46	0.37	0.46	-0.2
	OCK	0.64	0.26	0.22	-0.03		OCK	0.6	0.5	0.32	-0.33
	GMM	0.86	0.23				GMM	0.92	0.17
49	OK	0.61	0.28	0.31	-0.11	49	OK	0.52	0.34	0.39	-0.17
	OCK	0.66	0.26	0.26	-0.09		OCK	0.61	0.31	0.3	-0.14
	GMM	0.92	0.17				GMM	0.91	0.17
61	OK	0.67	0.25	0.26	-0.09	61	OK	0.54	0.33	0.43	-0.25
	OCK	0.72	0.24	0.21	-0.08		OCK	0.62	0.31	0.35	-0.23
	GMM	0.93	0.16				GMM	0.97	0.08
76	OK	0.73	0.24	0.11	-0.01	76	OK	0.58	0.32	0.38	-0.2
	OCK	0.76	0.22	0.08	0.01		OCK	0.66	0.29	0.3	-0.17
	GMM	0.84	0.23				GMM	0.96	0.12
91	OK	0.73	0.24	0.22	-0.1	91	OK	0.61	0.31	0.37	-0.23
	OCK	0.77	0.22	0.18	-0.08		OCK	0.68	0.29	0.3	-0.21
	GMM	0.95	0.14				GMM	0.98	0.08
107	OK	0.75	0.23	0.21	-0.1	107	OK	0.62	0.31	0.35	-0.2
	OCK	0.79	0.21	0.17	-0.08		OCK	0.68	0.28	0.29	-0.17
	GMM	0.96	0.13				GMM	0.97	0.11

*r_GMM ‒ r_geost.: performance comparison (r) between GMM and geostatistical approach (OK or OCK);.

**RMSE_GMM ‒RMSE_geost: performance comparison (RMSE) between GMM and geostatistical approach (OK or OCK)..

Figure 5. Root mean-square errors (RMSE) between the measured and predicted values using the validation data (n=30) at different sampling densities for training (n=30 to 107 in Table 2). (a) As, (b) Pb.

동일 샘플링 밀도 기준으로 정확도와 정밀도를 비교해 보았다. RMSE는 일부 훈련 데이터셋(76개의 비소)을 제외하고 OK > OCK > GMM 순으로 높았다(Table 2& Fig. 5). 비소의 GMM은 OK 대비 최소 0.01(n=76)에서 최대 0.11(n=49)의 RMSE 감소를 보였고, OCK 대비 최대 0.09(n=49)의 RMSE 감소를 보였다. 납의 GMM 결과는 OK 대비 최소 0.17(n=49)에서 최대 0.25(n=61)의 RMSE 감소가 있었고, OCK 대비 최소 0.14(n=49)에서 최대 0.33(n=30)의 RMSE 감소를 보였다.

상관계수 값은 GMM > OCK > OK 순으로 높았다(Table 2& Fig. 6). 비소의 GMM은 OK 대비 최소 0.11(n=76)에서 최대 0.31(n=49)의 상관계수 값을 증가시켰으며, OCK대비 최소 0.08(n=76)에서 최대 0.26(n=49)의 상관계수 값을 증가시켰다. 납의 GMM은 OK 대비 최소 0.35(n=107)에서 최대 0.46(n=30)의 상관계수 값이 증가하였으며, OCK 대비 최소 0.29(n=107)에서 최대 0.35(n=61)의 상관계수 값이 증가하였다.

Figure 6. Pearson correlation coefficients (r) between the measured and predicted values using the validation data (n=30) at different sampling densities for training (n=30 to 107 in Table 2): (a) As, (b) Pb.

상기 결과는 다변량을 활용한 GMM이 단변량 혹은 이변량을 사용한 지구통계기법보다 예측 정확도(RMSE) 및 상관계수(r)가 높음을 지시한다. 더군다나, 상관성이 높은 변수들을 함께 활용하는 것이 모델 예측력 향상에 도움을 줄 수 있다는 경험법칙(rule of thumb)이 본 연구지역과 같이 실제 오염대상 부지에서의 공간분포 추정에도 적용될 수 있음을 시사한다. 위를 종합하였을 때, GMM을 통한 다변량 기법이 단변량(OK) 혹은 이변량을 이용한 OCK에 비해 우수한 공간분포 추정 결과를 제공할 수 있음을 시사한다. 다만, OCK 및 GMM의 상호비교에 있어 추가 보조변수의 차이(OCK의 경우 2차원, GMM은 8차원의 데이터 활용)에 따른 공간추정 성능 차이를 유발할 수 있음을 유념해야 한다.

본 연구는 충남 서천군 옛 장항제련소 인근 토양오염부지 정밀조사 결과를 활용하여 가우시안 혼합 모형(GMM) 기반 다변량 지구과학 데이터의 공간추정 가능성을 평가하였다. GMM은 기존 지구통계 추론기법들인 OK나 OCK 대비 비소에서 최소 0.01에서 0.11까지, 납에서는 0.14에서 0.33까지 RMSE를 감소시켰으며, 상관계수 값은 비소에서 0.08에서 0.31까지, 납에서는 0.29에서 0.46까지 개선되었다. 이는 GMM이 다변량 공간구조의 이해를 바탕으로 토양오염분포(지오데이터의 공간모집단)를 예측하여 정확성과 정밀성을 높이고 강건한 예측결과를 제공할 수 있음을 시사한다.

위 결과들을 종합적으로 고려해보면, 제안된 GMM이 자연계(real fields) 상의 이질적이고 복잡한 자연현상들을 모사하는 데 효과적인 방법임을 알 수 있다. 향후 연구에서는 다변량 변수와 더불어 다중 해상도의 소프트 데이터(예: DEM, 물리탐사 데이터)를 GMM과 함께 활용하여 다양한 제약조건(예: 잡음원, 분석의 간섭효과 등)의 영향을 받는 지오데이터의 공간 분포 예측 성능을 향상시키는 연구를 수행하고자 한다.

본 논문은 한국지질자원연구원의 주요사업(22-3415, 22-3117, 22-3412-2) 및 환경산업기술원의 지중환경오염·위해관리기술개발사업(과제번호: 2018002440002)의 지원으로 수행되었으며, 이에 감사드립니다.

¹⁾ 지오데이터(Geo-data)는 지질자원분야에서 발생할 수 있는 모든 자료의 집합체로서, 자연현상으로부터 계측되는 자료, 가상모델로부터 해석을 위해 생성되는 자료 등 지질자원분야와 관련된 모든 시공간 자료를 일컫는다.

²⁾ 그림자 변수(shadow attributes)는 변수 중요도의 유의성을 평가하는 데 활용한다. 그림자 변수는 모든 변수의 복사본을 섞어서 변수 간 상관성을 제거한 새로운 생성 변수를 지시한다. 무작위 변동에서 발생한 변수인 그림자 변수보다 낮은 중요도를 가진 변수들의 경우, 모델에서 중요하지 않다고 판단한다.

³⁾ Boruta의 중요도는 mean decrease accuracy에 기반을 둔다. mean decrease accuracy는 각 변수를 제거할 때 모형이 얼마나 많은 정확도를 상실하는지를 점수화한다. 즉 변수 제거에 따른 모형의 정확도가 크게 낮아지면, 해당 변수는 예측에 중요한 인자로 간주할 수 있다.