Econ. Environ. Geol. 2022; 55(4): 353-366
Published online August 30, 2022
https://doi.org/10.9719/EEG.2022.55.4.353
© THE KOREAN SOCIETY OF ECONOMIC AND ENVIRONMENTAL GEOLOGY
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided original work is properly cited.
Spatial estimation of geoscience data (geo-data) is challenging due to spatial heterogeneity, data scarcity, and high dimensionality. A novel spatial estimation method is needed to consider the characteristics of geo-data. In this study, we proposed the application of Gaussian Mixture Model (GMM) among machine learning algorithms with multivariate data for robust spatial predictions. The performance of the proposed approach was tested through soil chemical concentration data from a former smelting area. The concentrations of As and Pb determined by ex-situ ICP-AES were the primary variables to be interpolated, while the other metal concentrations by ICP-AES and all data determined by in-situ portable X-ray fluorescence (PXRF) were used as auxiliary variables in GMM and ordinary cokriging (OCK). Among the multidimensional auxiliary variables, important variables were selected using a variable selection method based on the random forest. The results of GMM with important multivariate auxiliary data decreased the root mean-squared error (RMSE) down to 0.11 for As and 0.33 for Pb and increased the correlations (r) up to 0.31 for As and 0.46 for Pb compared to those from ordinary kriging and OCK using univariate or bivariate data. The use of GMM improved the performance of spatial interpretation of anthropogenic metals in soil. The multivariate spatial approach can be applied to understand complex and heterogeneous geological and geochemical features.
Keywords Gaussian Mixture Model (GMM), multivariate, geoscience data (geo-data), machine learning, soil contamination
김호림1 · 유순영2 · 윤성택2 · 김경호3 · 이군택4 · 이정호1 · 허철호1 · 류동우1*
1한국지질자원연구원
2고려대학교
3한국환경연구원
4서울대학교 NICEM
지구과학 데이터(지오데이터)의 공간 이질성, 희소성 및 고차원성으로 인해 공간 분포 추정에 어려움이 있다. 따라서 지구과학의 많은 응용 분야에서 지오데이터의 고유 특성을 고려할 수 있는 공간 추정 기법이 필요하다. 본 연구에서는 기계 학습 알고리즘 중 하나인 가우시안 혼합 모델(Gaussian Mixture Model; GMM)을 이용하여 공간 예측 방법을 제공하고자 하였다. 제안된 기법의 성능을 검증하기 위해, 옛 제련소 부지에서 휴대용 X선 형광분석기(PXRF) 및 유도결합플라즈마-원자방출분광법(ICPAES)을 이용하여 분석된 토양 농도 자료를 활용하였다. ICP-AES를 이용해 분석된 As와 Pb를 주변수로 하고, 나머지 자료는 보조변수로 활용하였다. 다차원의 보조변수 중 중요 변수를 선별하기 위해 랜덤포레스트 기반의 변수선택법을 적용하였다. ICPAES 및 PXRF를 통해 구축된 다변량 데이터를 사용한 GMM의 결과를 단변량 및 이변량 데이터를 사용한 정규 크리깅(Ordinary Kriging; OK) 및 정규 공동크리깅(Ordinary Co-Kriging; OCK)의 결과와 비교하였다. GMM의 결과는 OK 및 OCK의 결과보다 낮은 평균 제곱근 편차(RMSE; 비소는 최대 0.11 및 납은 0.33까지 향상)와 높은 상관관계(r; 비소는 최대 0.31 및 납은 0.46까지 향상)를 제공하였다. 이는 GMM을 사용할 경우 토양 오염의 범위 해석의 성능을 향상시킬 수 있음을 지시한다. 본 연구는 다변량 공간추정 접근법이복잡하고 이질적인 지질 및 지구 화학자료의 특징을 이해하는 데 효과적으로적용될 수 있음을증명하였다.
주요어 가우시안 혼합모형, 다변량, 지구과학데이터(지오데이터), 기계학습, 토양오염
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
예를 들면, 지질자원 분야의 관심변수(예: 대형재난, 광물자원 매장량 등)는 시공간적으로 드물게 발생할 뿐만 아니라 지하 수 미터에서 수 킬로미터에 다다르는 심부공간에서 발생한다. 다양한 제약조건(예: 막대한 시추 비용 등)에 따른 하드 데이터(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
본 연구에서는 다변량의 확률분포를 추정하고 이를 반영한 공간예측결과가 기존 지구통계기법 중 단변량을 활용한 정규 크리깅(Ordinary Kriging; OK)과 이변량을 활용한 정규 공동크리깅(Ordinary Co-Kriging; OCK))에 비해 정확성 및 정밀성을 얼마나 향상시키는 지를 평가하였다. 성능 평가를 위해 충남 서천군 옛 장항제련소 인근 토양오염 정밀조사 부지를 테스트 베드로 선정하였다. Kim
다변량 변수와 GMM을 이용한 공간분포 추정 방법의 성능을 평가하기 위해 충남 서천군 옛 장항제련소 인근 토양오염부지 정밀조사 결과를 활용하였다. 장항제련소(1936-1989년)는 운영 과정에서 발생한 오염물질로 인해 토양오염 및 농작물 피해 등 다양한 환경문제를 유발하였다.
Moon
토양정밀조사지침(KMOE, 2009b)에 따르면 600~1,000m2당 1개 지점을 조사해야 하며, 따라서 연구부지 면적(3,300m2)에서는 최대 6개 지점을 조사하면 되었다. 그러나 본 연구에서는 정밀한 토양오염 면적 산정과 GMM 기법의 신뢰성을 검토하기 위해 정밀조사지침 기준보다 많은 156개 지점에서 토양시료를 조사하였다. 기준필지를 대상으로 정방형의 5m 간격으로 시료를 채취하였다. 현장시료채취 및 분석은 서울대학교 NICEM(National Instrumentation Center for Environmental Management)에서 수행하였다.
실험실 분석용 시료(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
PXRF는 고체물질에 X선을 흡수(혹은 투과), 회절(혹은 산란)시키고, 다른 색의 X선을 만들어 주사하여 화학물질의 전 함량을 분석한다. 예를 들면, 특정원소에 X선을 주사하면 K 준위에 있던 전자들은 X선으로부터 에너지를 받아 들뜬 상태에 이르며, 들뜬 상태의 원자가 원래의 에너지 준위로 돌아올 때 특정 에너지 파장을 방출한다. 비소의 경우, L에서 K로 갈 때 10.53 KeV를, M에서 K로 갈 때 11.73 KeV를 방출한다. 이 방출파장은 각 원소마다 특정 값을 지니기 때문에 각 peak들의 에너지 값과 면적을 통해 원소의 종류와 함량을 평가할 수 있다.
그러나 현장 PXRF 분석은 토양 특성에 따른 여러 가지 간섭 효과에 의해 정밀도가 떨어진다는 문제가 있다. 예를 들어 토양 사이의 공극이나 물질들의 물리적 배열상태(matrix effect)에 따라 오차가 발생하며, 토양 내 수분함량이 높은 경우에도 간섭이 발생할 수 있어, 이를 극복하기 위해 많은 연구가 진행 중에 있다(Ryan
본 연구에서 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
다변량의 보조변수들을 이용하여 주변수인 납과 비소의 공간분포를 산정하기에 앞서, 예측 모형의 과적합(overfitting)을 방지하고, 모델 복잡도를 최소화하는 한편, 예측 속도를 향상시키기 위해 변수 선택 기법을 적용하였다. 변수 선택을 위해 랜덤포레스트 기반의 Boruta를 활용하였으며(Fig. 1 Step 1-②), Boruta는 다음과 같은 순서로 변수를 선정한다(Kursa and Rudnicki, 2010).
① 모든 변수의 복사본 생성(adding copies of all variables)
② 그림자 변수(shadow attributes)2) 생성: 변수 복사본을 무작위로 섞어서 변수 간 상관성을 없앰. 이들은 무작위성을 통해 만들어진 새로운 변수임
③ 그림자 변수를 포함한 표본들을 대상으로 랜덤포레스트 모형을 구동하고, 각 변수에 관한 정확도 손실(accuracy loss)을 표준화 점수(Z-score)로 계산
④ 그림자 변수들 중에서 최대 표준화 점수를 갖는 변수(maximum Z score among shadow attributes; MZSA)를 확인
⑤ MZSA를 기준으로, 입력 변수의 표준화 점수가 높으면 중요(important), 낮으면 중요하지 않다(unimportant)고 간주
⑥ 모든 변수의 중요도(importance)3)가 정해질 때까지 ②-⑤ 과정을 반복
일반적으로 확률분포의 추정은 모집단(population)의 모수(parameter; 예: 가우시안 분포의 경우 1차 모멘트인 표본평균 μi와 2차 모멘트인 표본분산 σi2)를 최대우도추정(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
GMM은 개별밀도함수를 전체 확률밀도함수의 성분(component)인 커널(kernel)로 간주한다. 다중 가우시안 분포를 이용할 수 있으므로 여러 개의 중심점을 갖는 다차원 데이터에 대해 견고한 모델링이 가능하다.
GMM을 통한 확률밀도함수 추정은 복수(m)의 가우시안 확률밀도함수에 대한 선형 결합으로 표현되며, 이를 식으로 표현하면 아래와 같다(Reynolds, 2009):
위 식에서
가우시안 확률밀도함수의 모수(parameter; θ)는 평균 μi, 분산 σi2, 그리고 가중치 α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-②).
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
본 연구에서 조사된 비소 및 중금속의 농도는 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
각 변수별 상관성은 Fig. 2와 같다. 도표의 대각선 방향은 각 원소의 분포도(히스토그램)와 적합(fitting) 결과(붉은색 선)를 함께 도시한 것이다. 대각선 상단은 변수별 상관계수를 표시하였으며, 숫자 위의 붉은 별의 개수는 신뢰수준을 의미한다. 대각선 하단은 변수별 이변량도(bivariate scatter plots)와 적합 결과를 함께 도시한 결과이다. 일부 항목에서(예: ICP-AES로 분석된 구리와 니켈과의 상관성 및 적합 결과) 이변량 분포 간 비선형성 및 이분산적(heteroscedastic) 특징이 관찰되며, 따라서 본 연구 데이터는 다변량 확률분포 모델링과 같은 접근방식이 필요하다.
ICP-AES 분석 결과 간의 상관성을 비교해보면, 토양 내 비소와 납(0.91), 비소와 구리(0.92), 납과 구리(0.96)에서 유의한 신뢰성과 높은 상관성을 보인다. 앞서 기초통계의 평균 및 환경부 기준 등을 고려하였을 때, 제련소에 의한 토양 내 비소와 납 오염은 파악하였으나 구리는 대부분 시료에서 우려기준을 초과하지 않았기 때문에 구리 농도에 대한 제련소의 영향이 없는 것으로 판단하였다. 그러나 상관성을 고려할 시 비소와 납뿐만 아니라 구리도 제련활동의 영향을 받는 것으로 간주될 수 있다. 실제 연구 지역의 구리 농도는 국내 토양의 배경 농도(15.3 mg/kg; KMOE, 2011; Yang
구리는 자연계에서 자연동으로 산출되며, 주요광석은 황동석(CuFeS2), 휘동석(Cu2S), 공작석(Cu2CO3(OH)), 적동석(Cu2O) 등이 있다. 지각 내 함유량은 55~65 mg/kg으로 추정되며, 화성암은 55 mg/kg, 퇴적암은 5~45 mg/kg, 토양은 20 ppm로 추정하고 있다(Koljonen, 1992; Galán
비소와 납 농도의 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를 이용하여 분석한 비소, 납, 구리 농도는 각각 비소, 납, 구리의 공간분포 예측을 위한 보조변수로 충분히 활용할 수 있을 것으로 보인다.
랜덤포레스트 기반의 Boruta를 이용하여 보조변수들이 비소와 납 농도에 미치는 중요도를 계산하였다(Fig. 1; Step 1-②). 먼저 비소 예측을 위한 변수 중요도 계산 결과, 비소를 제외한 9개 설명변수 중에서 PXRF 니켈과 아연은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었다. 확정된 보조변수들 중에서 ICP-AES 납과 구리가 상위 영향 인자(상대적으로 높은 중요도)로 파악된다. 이는 앞선 상관성 분석 결과와도 일치하며, 제련과정에서 발생한 비소, 납, 구리의 영향이 함께 반영된 것으로 판단할 수 있다.
납에서도 마찬가지로 상위 중요도로 분류된 변수는 ICP-AES 비소와 구리이다. PXRF 니켈은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었으며, PXRF 아연은 최대 그림자 변수(shadowMax)와 유사한 중요도를 보여 잠정(tentative) 변수로 분류되었으나, 중요도 이력 정보 및 이상치 분포 등을 종합적으로 고려하여 납의 공간 예측 시 PXRF 아연은 활용하지 않았다.
따라서 GMM을 이용한 비소 농도의 공간 분포에는 ICP-AES 납, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하고, 납 농도의 공간 분포 추정에는 ICP-AES 비소, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하였다.
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은 이질성을 내포한 데이터 모집단의 통계 구조를 모사·반영하고 확률밀도 추정 시 강건한 모델링을 제공할 수 있음을 지시한다.
기법별 공간분포 예측결과 중 비소의 공간분포 추정 결과를 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)와 비교하였다.
OK, OCK, GMM 분석결과 모두 서쪽에 분포한 고농도 비소 오염지역(핫스폿; hot spot)을 효과적으로 추정하였다. 그러나 OK는 서쪽 중앙부에 위치한 저농도 지점(콜드스폿; cold spot)을 예측하는 데 실패하였다. 이는 콜드스폿에 대한 주변수(ICP-AES의 As 49개) 정보가 빠져있기 때문이다.
반면, OCK나 GMM은 주요 변수(ICE-AES로 분석한 As 농도)가 저농도 지점에 대한 정보를 갖지 않음에도 불구하고 보조변수의 공간적 상관성을 이용하여 원래 자료의 콜드스폿을 효율적으로 찾을 수 있었다. 특히 GMM은 비소 대책기준 미만 농도에 준하는 경계(콜드스폿)를 예측하였으며 대조군 자료와도 가장 유사한 결과를 보였다(Fig. 4c-d).
환경재해예측이나 자원탐사 분야에서는 핫스폿 혹은 콜드스폿 등 공간이상치(outlier)를 파악하는 것이 중요하다. 위와 같이 원 변수가 중요지점에 대한 정보를 모르는 경우에도 다변량 변수 간의 공간적인 관계성을 토대로 한 OCK나 GMM의 기법들이 원 변수 공간분포 추정에 매우 효과적일 수 있음을 시사한다.
단변량(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 |
*rGMM ‒ rgeost.: performance comparison (r) between GMM and geostatistical approach (OK or OCK);
**RMSEGMM ‒RMSEgeost: performance comparison (RMSE) between GMM and geostatistical approach (OK or OCK).
동일 샘플링 밀도 기준으로 정확도와 정밀도를 비교해 보았다. 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)의 상관계수 값이 증가하였다.
상기 결과는 다변량을 활용한 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)의 지원으로 수행되었으며, 이에 감사드립니다.
1) 지오데이터(Geo-data)는 지질자원분야에서 발생할 수 있는 모든 자료의 집합체로서, 자연현상으로부터 계측되는 자료, 가상모델로부터 해석을 위해 생성되는 자료 등 지질자원분야와 관련된 모든 시공간 자료를 일컫는다.
2) 그림자 변수(shadow attributes)는 변수 중요도의 유의성을 평가하는 데 활용한다. 그림자 변수는 모든 변수의 복사본을 섞어서 변수 간 상관성을 제거한 새로운 생성 변수를 지시한다. 무작위 변동에서 발생한 변수인 그림자 변수보다 낮은 중요도를 가진 변수들의 경우, 모델에서 중요하지 않다고 판단한다.
3) Boruta의 중요도는 mean decrease accuracy에 기반을 둔다. mean decrease accuracy는 각 변수를 제거할 때 모형이 얼마나 많은 정확도를 상실하는지를 점수화한다. 즉 변수 제거에 따른 모형의 정확도가 크게 낮아지면, 해당 변수는 예측에 중요한 인자로 간주할 수 있다.
Econ. Environ. Geol. 2022; 55(4): 353-366
Published online August 30, 2022 https://doi.org/10.9719/EEG.2022.55.4.353
Copyright © THE KOREAN SOCIETY OF ECONOMIC AND ENVIRONMENTAL GEOLOGY.
Ho-Rim Kim1, Soonyoung Yu2, Seong-Taek Yun2, Kyoung-Ho Kim3, Goon-Taek Lee4, Jeong-Ho Lee1, Chul-Ho Heo1, Dong-Woo Ryu1,*
1Korea Institute of Geoscience and Mineral Resources, Republic of Korea
2Korea University, Republic of Korea
3Korea Environment Institute, Republic of Korea
4National Instrumentation Center for Environmental Management, Seoul National University, Republic of Korea
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided original work is properly cited.
Spatial estimation of geoscience data (geo-data) is challenging due to spatial heterogeneity, data scarcity, and high dimensionality. A novel spatial estimation method is needed to consider the characteristics of geo-data. In this study, we proposed the application of Gaussian Mixture Model (GMM) among machine learning algorithms with multivariate data for robust spatial predictions. The performance of the proposed approach was tested through soil chemical concentration data from a former smelting area. The concentrations of As and Pb determined by ex-situ ICP-AES were the primary variables to be interpolated, while the other metal concentrations by ICP-AES and all data determined by in-situ portable X-ray fluorescence (PXRF) were used as auxiliary variables in GMM and ordinary cokriging (OCK). Among the multidimensional auxiliary variables, important variables were selected using a variable selection method based on the random forest. The results of GMM with important multivariate auxiliary data decreased the root mean-squared error (RMSE) down to 0.11 for As and 0.33 for Pb and increased the correlations (r) up to 0.31 for As and 0.46 for Pb compared to those from ordinary kriging and OCK using univariate or bivariate data. The use of GMM improved the performance of spatial interpretation of anthropogenic metals in soil. The multivariate spatial approach can be applied to understand complex and heterogeneous geological and geochemical features.
Keywords Gaussian Mixture Model (GMM), multivariate, geoscience data (geo-data), machine learning, soil contamination
김호림1 · 유순영2 · 윤성택2 · 김경호3 · 이군택4 · 이정호1 · 허철호1 · 류동우1*
1한국지질자원연구원
2고려대학교
3한국환경연구원
4서울대학교 NICEM
지구과학 데이터(지오데이터)의 공간 이질성, 희소성 및 고차원성으로 인해 공간 분포 추정에 어려움이 있다. 따라서 지구과학의 많은 응용 분야에서 지오데이터의 고유 특성을 고려할 수 있는 공간 추정 기법이 필요하다. 본 연구에서는 기계 학습 알고리즘 중 하나인 가우시안 혼합 모델(Gaussian Mixture Model; GMM)을 이용하여 공간 예측 방법을 제공하고자 하였다. 제안된 기법의 성능을 검증하기 위해, 옛 제련소 부지에서 휴대용 X선 형광분석기(PXRF) 및 유도결합플라즈마-원자방출분광법(ICPAES)을 이용하여 분석된 토양 농도 자료를 활용하였다. ICP-AES를 이용해 분석된 As와 Pb를 주변수로 하고, 나머지 자료는 보조변수로 활용하였다. 다차원의 보조변수 중 중요 변수를 선별하기 위해 랜덤포레스트 기반의 변수선택법을 적용하였다. ICPAES 및 PXRF를 통해 구축된 다변량 데이터를 사용한 GMM의 결과를 단변량 및 이변량 데이터를 사용한 정규 크리깅(Ordinary Kriging; OK) 및 정규 공동크리깅(Ordinary Co-Kriging; OCK)의 결과와 비교하였다. GMM의 결과는 OK 및 OCK의 결과보다 낮은 평균 제곱근 편차(RMSE; 비소는 최대 0.11 및 납은 0.33까지 향상)와 높은 상관관계(r; 비소는 최대 0.31 및 납은 0.46까지 향상)를 제공하였다. 이는 GMM을 사용할 경우 토양 오염의 범위 해석의 성능을 향상시킬 수 있음을 지시한다. 본 연구는 다변량 공간추정 접근법이복잡하고 이질적인 지질 및 지구 화학자료의 특징을 이해하는 데 효과적으로적용될 수 있음을증명하였다.
주요어 가우시안 혼합모형, 다변량, 지구과학데이터(지오데이터), 기계학습, 토양오염
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
예를 들면, 지질자원 분야의 관심변수(예: 대형재난, 광물자원 매장량 등)는 시공간적으로 드물게 발생할 뿐만 아니라 지하 수 미터에서 수 킬로미터에 다다르는 심부공간에서 발생한다. 다양한 제약조건(예: 막대한 시추 비용 등)에 따른 하드 데이터(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
본 연구에서는 다변량의 확률분포를 추정하고 이를 반영한 공간예측결과가 기존 지구통계기법 중 단변량을 활용한 정규 크리깅(Ordinary Kriging; OK)과 이변량을 활용한 정규 공동크리깅(Ordinary Co-Kriging; OCK))에 비해 정확성 및 정밀성을 얼마나 향상시키는 지를 평가하였다. 성능 평가를 위해 충남 서천군 옛 장항제련소 인근 토양오염 정밀조사 부지를 테스트 베드로 선정하였다. Kim
다변량 변수와 GMM을 이용한 공간분포 추정 방법의 성능을 평가하기 위해 충남 서천군 옛 장항제련소 인근 토양오염부지 정밀조사 결과를 활용하였다. 장항제련소(1936-1989년)는 운영 과정에서 발생한 오염물질로 인해 토양오염 및 농작물 피해 등 다양한 환경문제를 유발하였다.
Moon
토양정밀조사지침(KMOE, 2009b)에 따르면 600~1,000m2당 1개 지점을 조사해야 하며, 따라서 연구부지 면적(3,300m2)에서는 최대 6개 지점을 조사하면 되었다. 그러나 본 연구에서는 정밀한 토양오염 면적 산정과 GMM 기법의 신뢰성을 검토하기 위해 정밀조사지침 기준보다 많은 156개 지점에서 토양시료를 조사하였다. 기준필지를 대상으로 정방형의 5m 간격으로 시료를 채취하였다. 현장시료채취 및 분석은 서울대학교 NICEM(National Instrumentation Center for Environmental Management)에서 수행하였다.
실험실 분석용 시료(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
PXRF는 고체물질에 X선을 흡수(혹은 투과), 회절(혹은 산란)시키고, 다른 색의 X선을 만들어 주사하여 화학물질의 전 함량을 분석한다. 예를 들면, 특정원소에 X선을 주사하면 K 준위에 있던 전자들은 X선으로부터 에너지를 받아 들뜬 상태에 이르며, 들뜬 상태의 원자가 원래의 에너지 준위로 돌아올 때 특정 에너지 파장을 방출한다. 비소의 경우, L에서 K로 갈 때 10.53 KeV를, M에서 K로 갈 때 11.73 KeV를 방출한다. 이 방출파장은 각 원소마다 특정 값을 지니기 때문에 각 peak들의 에너지 값과 면적을 통해 원소의 종류와 함량을 평가할 수 있다.
그러나 현장 PXRF 분석은 토양 특성에 따른 여러 가지 간섭 효과에 의해 정밀도가 떨어진다는 문제가 있다. 예를 들어 토양 사이의 공극이나 물질들의 물리적 배열상태(matrix effect)에 따라 오차가 발생하며, 토양 내 수분함량이 높은 경우에도 간섭이 발생할 수 있어, 이를 극복하기 위해 많은 연구가 진행 중에 있다(Ryan
본 연구에서 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
다변량의 보조변수들을 이용하여 주변수인 납과 비소의 공간분포를 산정하기에 앞서, 예측 모형의 과적합(overfitting)을 방지하고, 모델 복잡도를 최소화하는 한편, 예측 속도를 향상시키기 위해 변수 선택 기법을 적용하였다. 변수 선택을 위해 랜덤포레스트 기반의 Boruta를 활용하였으며(Fig. 1 Step 1-②), Boruta는 다음과 같은 순서로 변수를 선정한다(Kursa and Rudnicki, 2010).
① 모든 변수의 복사본 생성(adding copies of all variables)
② 그림자 변수(shadow attributes)2) 생성: 변수 복사본을 무작위로 섞어서 변수 간 상관성을 없앰. 이들은 무작위성을 통해 만들어진 새로운 변수임
③ 그림자 변수를 포함한 표본들을 대상으로 랜덤포레스트 모형을 구동하고, 각 변수에 관한 정확도 손실(accuracy loss)을 표준화 점수(Z-score)로 계산
④ 그림자 변수들 중에서 최대 표준화 점수를 갖는 변수(maximum Z score among shadow attributes; MZSA)를 확인
⑤ MZSA를 기준으로, 입력 변수의 표준화 점수가 높으면 중요(important), 낮으면 중요하지 않다(unimportant)고 간주
⑥ 모든 변수의 중요도(importance)3)가 정해질 때까지 ②-⑤ 과정을 반복
일반적으로 확률분포의 추정은 모집단(population)의 모수(parameter; 예: 가우시안 분포의 경우 1차 모멘트인 표본평균 μi와 2차 모멘트인 표본분산 σi2)를 최대우도추정(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
GMM은 개별밀도함수를 전체 확률밀도함수의 성분(component)인 커널(kernel)로 간주한다. 다중 가우시안 분포를 이용할 수 있으므로 여러 개의 중심점을 갖는 다차원 데이터에 대해 견고한 모델링이 가능하다.
GMM을 통한 확률밀도함수 추정은 복수(m)의 가우시안 확률밀도함수에 대한 선형 결합으로 표현되며, 이를 식으로 표현하면 아래와 같다(Reynolds, 2009):
위 식에서
가우시안 확률밀도함수의 모수(parameter; θ)는 평균 μi, 분산 σi2, 그리고 가중치 α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-②).
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
본 연구에서 조사된 비소 및 중금속의 농도는 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
각 변수별 상관성은 Fig. 2와 같다. 도표의 대각선 방향은 각 원소의 분포도(히스토그램)와 적합(fitting) 결과(붉은색 선)를 함께 도시한 것이다. 대각선 상단은 변수별 상관계수를 표시하였으며, 숫자 위의 붉은 별의 개수는 신뢰수준을 의미한다. 대각선 하단은 변수별 이변량도(bivariate scatter plots)와 적합 결과를 함께 도시한 결과이다. 일부 항목에서(예: ICP-AES로 분석된 구리와 니켈과의 상관성 및 적합 결과) 이변량 분포 간 비선형성 및 이분산적(heteroscedastic) 특징이 관찰되며, 따라서 본 연구 데이터는 다변량 확률분포 모델링과 같은 접근방식이 필요하다.
ICP-AES 분석 결과 간의 상관성을 비교해보면, 토양 내 비소와 납(0.91), 비소와 구리(0.92), 납과 구리(0.96)에서 유의한 신뢰성과 높은 상관성을 보인다. 앞서 기초통계의 평균 및 환경부 기준 등을 고려하였을 때, 제련소에 의한 토양 내 비소와 납 오염은 파악하였으나 구리는 대부분 시료에서 우려기준을 초과하지 않았기 때문에 구리 농도에 대한 제련소의 영향이 없는 것으로 판단하였다. 그러나 상관성을 고려할 시 비소와 납뿐만 아니라 구리도 제련활동의 영향을 받는 것으로 간주될 수 있다. 실제 연구 지역의 구리 농도는 국내 토양의 배경 농도(15.3 mg/kg; KMOE, 2011; Yang
구리는 자연계에서 자연동으로 산출되며, 주요광석은 황동석(CuFeS2), 휘동석(Cu2S), 공작석(Cu2CO3(OH)), 적동석(Cu2O) 등이 있다. 지각 내 함유량은 55~65 mg/kg으로 추정되며, 화성암은 55 mg/kg, 퇴적암은 5~45 mg/kg, 토양은 20 ppm로 추정하고 있다(Koljonen, 1992; Galán
비소와 납 농도의 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를 이용하여 분석한 비소, 납, 구리 농도는 각각 비소, 납, 구리의 공간분포 예측을 위한 보조변수로 충분히 활용할 수 있을 것으로 보인다.
랜덤포레스트 기반의 Boruta를 이용하여 보조변수들이 비소와 납 농도에 미치는 중요도를 계산하였다(Fig. 1; Step 1-②). 먼저 비소 예측을 위한 변수 중요도 계산 결과, 비소를 제외한 9개 설명변수 중에서 PXRF 니켈과 아연은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었다. 확정된 보조변수들 중에서 ICP-AES 납과 구리가 상위 영향 인자(상대적으로 높은 중요도)로 파악된다. 이는 앞선 상관성 분석 결과와도 일치하며, 제련과정에서 발생한 비소, 납, 구리의 영향이 함께 반영된 것으로 판단할 수 있다.
납에서도 마찬가지로 상위 중요도로 분류된 변수는 ICP-AES 비소와 구리이다. PXRF 니켈은 그림자 변수보다 낮은 중요도를 보여 변수 선정에서 기각되었으며, PXRF 아연은 최대 그림자 변수(shadowMax)와 유사한 중요도를 보여 잠정(tentative) 변수로 분류되었으나, 중요도 이력 정보 및 이상치 분포 등을 종합적으로 고려하여 납의 공간 예측 시 PXRF 아연은 활용하지 않았다.
따라서 GMM을 이용한 비소 농도의 공간 분포에는 ICP-AES 납, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하고, 납 농도의 공간 분포 추정에는 ICP-AES 비소, 구리, 니켈, 아연 및 PXRF 비소, 납, 구리를 보조변수로 활용하였다.
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은 이질성을 내포한 데이터 모집단의 통계 구조를 모사·반영하고 확률밀도 추정 시 강건한 모델링을 제공할 수 있음을 지시한다.
기법별 공간분포 예측결과 중 비소의 공간분포 추정 결과를 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)와 비교하였다.
OK, OCK, GMM 분석결과 모두 서쪽에 분포한 고농도 비소 오염지역(핫스폿; hot spot)을 효과적으로 추정하였다. 그러나 OK는 서쪽 중앙부에 위치한 저농도 지점(콜드스폿; cold spot)을 예측하는 데 실패하였다. 이는 콜드스폿에 대한 주변수(ICP-AES의 As 49개) 정보가 빠져있기 때문이다.
반면, OCK나 GMM은 주요 변수(ICE-AES로 분석한 As 농도)가 저농도 지점에 대한 정보를 갖지 않음에도 불구하고 보조변수의 공간적 상관성을 이용하여 원래 자료의 콜드스폿을 효율적으로 찾을 수 있었다. 특히 GMM은 비소 대책기준 미만 농도에 준하는 경계(콜드스폿)를 예측하였으며 대조군 자료와도 가장 유사한 결과를 보였다(Fig. 4c-d).
환경재해예측이나 자원탐사 분야에서는 핫스폿 혹은 콜드스폿 등 공간이상치(outlier)를 파악하는 것이 중요하다. 위와 같이 원 변수가 중요지점에 대한 정보를 모르는 경우에도 다변량 변수 간의 공간적인 관계성을 토대로 한 OCK나 GMM의 기법들이 원 변수 공간분포 추정에 매우 효과적일 수 있음을 시사한다.
단변량(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 |
*rGMM ‒ rgeost.: performance comparison (r) between GMM and geostatistical approach (OK or OCK);.
**RMSEGMM ‒RMSEgeost: performance comparison (RMSE) between GMM and geostatistical approach (OK or OCK)..
동일 샘플링 밀도 기준으로 정확도와 정밀도를 비교해 보았다. 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)의 상관계수 값이 증가하였다.
상기 결과는 다변량을 활용한 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)의 지원으로 수행되었으며, 이에 감사드립니다.
1) 지오데이터(Geo-data)는 지질자원분야에서 발생할 수 있는 모든 자료의 집합체로서, 자연현상으로부터 계측되는 자료, 가상모델로부터 해석을 위해 생성되는 자료 등 지질자원분야와 관련된 모든 시공간 자료를 일컫는다.
2) 그림자 변수(shadow attributes)는 변수 중요도의 유의성을 평가하는 데 활용한다. 그림자 변수는 모든 변수의 복사본을 섞어서 변수 간 상관성을 제거한 새로운 생성 변수를 지시한다. 무작위 변동에서 발생한 변수인 그림자 변수보다 낮은 중요도를 가진 변수들의 경우, 모델에서 중요하지 않다고 판단한다.
3) Boruta의 중요도는 mean decrease accuracy에 기반을 둔다. mean decrease accuracy는 각 변수를 제거할 때 모형이 얼마나 많은 정확도를 상실하는지를 점수화한다. 즉 변수 제거에 따른 모형의 정확도가 크게 낮아지면, 해당 변수는 예측에 중요한 인자로 간주할 수 있다.
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.
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 |
*rGMM ‒ rgeost.: performance comparison (r) between GMM and geostatistical approach (OK or OCK);.
**RMSEGMM ‒RMSEgeost: performance comparison (RMSE) between GMM and geostatistical approach (OK or OCK)..
Ju Young Park, Sun Young Park, Jiyoung Choi, Sungil Kim, Yuri Kim, Bo Yeon Yi, Kyungbook Lee
Econ. Environ. Geol. 2024; 57(5): 529-537Kyoungeun Lee, Jaehyung Yu, Chanhyeok Park, Trung Hieu Pham
Econ. Environ. Geol. 2024; 57(4): 353-362Kalaivanan K, Vellingiri J
Econ. Environ. Geol. 2024; 57(3): 329-342