2.1 연구대상지 및 분석 데이터

영산강은 전라남도 담양에서 발원하여 서해로 흐르는 강으로, 광주광역시 전역과 전라남도 7개 시⋅군을 지난다. 2012년 본류에 설치된 승촌보와 죽산보는 영산강의 중류인 광주광역시 남구 승촌동, 전라남도 나주시 다시면에 위치하고 있다(Fig. 1). 영산강의 중권역 대표 지점으로는 광주, 지석천, 나주 등이 있으며, 보 대표 지점으로는 승촌보에 가장 가깝게 위치한 광산(영본B, 용봉), 죽산보에 가장 가깝게 위치한 죽산(나주) 지점이 있다.

Fig. 1. Map of Yeongsan River Study Sites

승촌보 상류 구간에 위치한 공공하수처리장에서는 일평균 72만 톤의 방류수가 영산강으로 유입되고 있으며 이로 인해 승촌보는 방류수의 직접적인 영향을 받고 있다(^{Seo et al., 2018}). 따라서, 본 연구에서는 1989년에 전라남도 나주시 노안면 학산리, 승촌보 상류 약 500m 지점에 설치된 광산 수질측정망 데이터를 활용하여 승촌보 개방 전⋅후 장기간에 걸친 수질 변화를 분석하였다.

현재 광산 수질측정망은 영산강물환경연구소에서 관리되고 있으며, 수질 데이터는 환경부에서 운영하는 물환경정보시스템(Water Environment Information System)의 일별 데이터 중 승촌보가 설치된 후 개방 전 2013년부터 2017년까지 데이터와 보 완전 개방 이후 2018년부터 2021년까지의 데이터를 분석에 이용하였다(^{ME, 2022}). 본 연구에서는 대상 수질 지표 항목으로 Water Temperature, Dissolved Oxygen (DO), Biochemical Oxygen Demand (BOD) Chemical Oxygen Demand (COD), Suspended Solids (SS), Total Nitrogen (TN), Nitrate (NO₃-N), Ammonia (NH₃-N), Total Phosphorus (TP), Phosphate (PO₄-P), potential of Hydrogen (pH) 및 Chl-a를 사용하였으며, 기상 환경이 Chl-a 농도에 미치는 영향을 고려하기 위해 기상청에서 운영하는 기상자료개방포털(Open MET Data Portal) 데이터 중 2013년부터 2021년까지의 광주 지점의 평균 기온, 강우량, 최대 풍속, 최대 일사량을 일별로 측정한 데이터를 분석 항목에 포함하였다(^{KMA, 2022}). 또한, 2018년 승촌보가 완전 개방되면서 수위 및 저수량의 변화가 관측되어 수위 변화에 따른 Chl-a 농도 변화를 살펴보기 위해 수자원공사에서 운영하는 물정보포털(My Water)의 승촌보 수위 데이터를 분석에 활용하였다(^{K-water, 2022}).

2.2 분석 방법

2.2.3 RF

RF는 결정 트리(decision tree)의 과적합(overfitting) 문제를 보완하기 위해 ^{Breiman (2001)}에 의해 소개된 앙상블(ensemble) 알고리즘의 하나로 분류 및 회귀 분석 등 다양한 문제에 폭넓게 사용되고 있다(^{Kim, Lee et al., 2021}). 앙상블이란 주어진 자료를 통해 다수의 예측 모델을 만들어 학습한 뒤, 모델들을 하나로 결합하여 더 강한 최종 예측 모델을 만들어내는 기법이다(^{Kang and Park, 2019}). 앙상블 방법은 크게 배깅(Bootstrap Aggregating)과 부스팅(Boosting)으로 나눌 수 있다. RF는 이 중 배깅을 사용하는 대표적인 알고리즘으로 Fig. 2의 과정을 통해 설계되어 다음 네 단계로 요약될 수 있다(^{Dittman et al., 2015}; ^{Kang and Park, 2019}).

Fig. 2. Random Forest Algorithm

• (1) 입력 데이터로부터 부트스트랩(Bootstrap) 샘플링을 통해 임의의 부트스트랩 표본을 추출한다. 이때, 부트스트랩이란 입력 데이터로부터 중복을 허용하여 입력 데이터와 같은 크기의 새로운 표본 집합을 생성하는 기법을 말한다.

• (2) 부트스트랩 표본과 함께 전체 특성(feature)중 랜덤하게 m개를 선택하여 최적의 분할을 만드는 결정 트리를 생성한다.

• (3) 위 (1), (2)번 과정을 N번 반복하여 총 N개의 결정 트리를 생성한다.

• (4) N개의 결정 트리를 하나의 분류기로 결합하여(Aggregate) RF 모델을 만들고, 각 결정 트리의 결과를 평균하여 최종 결과를 도출한다.

위 과정을 통해 생성된 RF는 중복을 허용하여 표본을 추출하기 때문에 각 결정 트리가 예측한 결과가 비상관화(decorrelation)되어 일반화된 결과를 도출할 뿐 아니라 노이즈가 포함된 데이터에 대해서도 강한 특성을 갖게 된다(^{Kim and Seo, 2020}; ^{Lee, Park et al., 2020}). 이로 인해 RF는 범주형 변수뿐 아니라 연속형 변수 예측에도 높은 예측력을 보여 다양한 분야에서 활용되고 있다(^{Kim, Lee et al., 2021}).

RF의 주요 하이퍼파라미터는 무작위성을 제어하기 위한 하이퍼파라미터와 복잡도를 제어하기 위한 하이퍼파라미터로 구분된다(^{Müller and Guido, 2016}). 무작위성을 제어하기 위한 하이퍼파라미터에는 결정 트리의 개수(n_estimators)와 각 노드를 가장 잘 분류할 수 있는 최대 특성의 개수(max_features)가 있으며, 복잡도를 제어하기 위한 하이퍼파라미터로는 트리의 최대 깊이(max_depth), 결정 트리의 말단인 리프 노드의 최대 개수(max_leaf_nodes), 노드에 넣을 수 있는 최소 샘플의 개수(min_samples_split), 리프 노드까지의 최소한의 샘플 개수(min_samples_leaf) 등이 있다(^{Tyralis et al., 2019}).

2.3 데이터 전처리 및 스케일링

데이터 전처리(Data pre-processing)는 머신러닝 모델을 통해 학습을 진행하기 전, 분석에 적합한 형태로 데이터를 변형하는 과정으로 머신러닝 모델의 성능을 높이기 위해서는 초기 데이터를 가공하는 과정이 필요하다(^{Lee and Sun, 2020}). 광산 수질측정망의 일별 데이터는 측정 날짜와 간격이 일정하지 않고 월 4~5회 측정되어, 일별 자료가 있는 기상, 수위 데이터와의 시간적 해상도를 맞추기 위해 수질 데이터의 시간적 해상도를 높이는 보간 방법을 진행하였다. ^{Gnauck (2004)}은 불규칙하게 측정된 수질 데이터를 4가지 보간 방법을 통해 일일 데이터로 보간하여 수질 모델링 분석에 적합한 규칙적인 데이터 셋을 구성하였다. ^{Gnauck (2004)}의 연구 결과 선형 보간(linear interpolation)이 수질 데이터 사이의 간격을 채우기에 가장 적합함을 확인할 수 있었으며, 여러 선행 연구에서도 데이터 보간 시 선형 보간 방법을 이용하여, 본 연구에서도 수질 데이터는 초기 데이터를 일별 데이터로 선형 보간한 후 분석에 사용하였다(^{Junninen et al., 2004}; ^{Karace and baleanu, 2020}; ^{Lepot et al., 2017}). 기상 데이터 인자인 평균 기온, 강우량, 최대 풍속, 최대 일사량과 수위 데이터의 경우 일별 측정 데이터가 존재하여 선형 보간하지 않고 초기 데이터를 그대로 사용하였다. 전처리를 마친 2013년부터 2021년까지 총 3287개의 일별 데이터를 전체 데이터로 구성하였으며, 이 중 목적 변수인 Chl-a 농도는 출력 데이터(output data)로, 설명 변수인 나머지 항목들을 입력 데이터(input data)로 사용하였다(Table 1). 또한, 수질측정 데이터에서 필연적으로 발생하는 이상치, 결측값 및 노이즈 등을 제거하기 위하여 Savitzky-Golay 필터를 이용하여 평활화(smoothing) 작업을 진행하였다. Savitzky-Golay 필터는 n차 다항식에 최소제곱법(method of least squares)을 이용하여 데이터를 피팅(fitting)하는 과정을 통해 데이터의 손상을 최소화하여 노이즈를 제거하는 방법이다(^{Savitzky and Golay, 1964}). 12개의 수질 데이터는 Savitzky-Golay 필터를 이용하여 노이즈를 제거하였으나, 기상 데이터 중 강우량의 경우 평활화 작업을 진행하면 강우량이 없는 무강우일에도 소량의 강우가 있는 것으로 보정되는 등 데이터의 중요 의미가 훼손될 가능성이 있다. 이에 기상 데이터와 수위 데이터는 전체 데이터 중 상위 0.05% 이상의 백분위에 해당되는 값을 가진 데이터를 0.1%의 백분위 값으로 대체하는 보정을 진행하였다(Fig. 3, Fig. 4). Fig. 3, Fig. 4는 Savitzky-Golay 필터를 이용하여 노이즈를 제거하기 전⋅후를 나타내며 이때, before는 노이즈를 제거하기 전 데이터, after는 노이즈 제거 후 데이터, drop은 노이즈 제거 전⋅후 차이를 의미한다. 최종적으로 데이터의 단위를 0~1사이로 정규화(Normalization)하는 데이터 스케일링을 통해 데이터 간 단위 차이로 인한 분석 결과의 오차가 생기지 않도록 정규화 된 데이터를 분석 데이터로 이용하였다. 이후 분석 데이터는 모델 학습을 위해 Scikit-Learn의 model_selection 패키지 내 train_test_split 모듈을 이용하여 Train data와 Test data로 구분하였다. 2013년부터 2021년까지의 총 3287개의 일별 데이터 중 약 70%인 임의의 2301개의 데이터를 Train data로, 나머지 30%인 986개의 데이터를 Test data로 나누어 DNN, RF 및 XGBoost 모델 학습을 진행하였다.

Fig. 3. Plots of input variables recorded at the study area. (before: cyan line and raw data, after: black line and smoothed data, drop: the difference between before and after)

Fig. 4. Plot of output variable recorded at the study area. (before: cyan line and raw data, after: black line and smoothed data, drop: the difference between before and after)

3.1 상관성 분석

SPSS (ver. 25.0)을 이용하여 Shapiro-Wilk과 Kolmogorov- Smirnov검증을 통해 분석 데이터의 정규성을 확인 후, 본 연구에는 스피어만 상관분석이 적합하다고 판단하여 스피어만 상관분석을 실시하였다. Chl-a와 수질, 기상, 수량 등 다양한 환경인자와의 상관 관계를 파악하기 위하여 실시한 상관분석 결과 중 p-value값이 0.05 미만인 값을 Fig. 5의 히트맵(heatmap) 그래프를 통해 시각화하였다.

Fig. 5. Correlation analysis of input and output data characteristics using heat map plot (p<0.05)

Chl-a와 가장 높은 양의 상관성을 가진 입력 변수는 BOD, COD, SS, pH 순으로 각각 0.62, 0.60, 0.42, 0.41의 강한 상관관계를 보였다. 유기물질을 나타내는 수질 지표인 BOD, COD는 조류가 성장함에 따라 수체 내부에서 발생하는 자생 오염물질(autochthonous matter)과의 관련성을 추정할 수 있으며(^{Kim et al., 2013}; ^{Lee, Park et al., 2020}), ^{Jung and Kim (2017b)}의 연구에서는 조류가 증식함에 따라 광합성에 의해 수체 내 탄산염 및 중탄산염이 흡수되어 pH가 증가한 것으로 보고하였다.

PO₄-P는 –0.42로 강한 음의 상관관계를 보였고, ^{Jung and Kim. (2017a)}과 ^{Park et al. (2008)}의 연구에 따르면, 식물성 플랑크톤 현존량이 높은 시기에 PO₄-P는 조류 증식과정에서 영양물질로 소모되는 것으로 보고하였다. 이에 반해 Water Temperature, NO₃-N, Average Temperature는 상관계수가 0.1 미만으로 Chl-a와 유의미한 상관관계가 나타나지 않았다.

3.2 DNN, RF 및 XGBoost 모델 최적화 과정

본 연구에서 Chl-a 농도를 예측에 사용한 DNN 모델은 Python 3.9.12 환경에서 TensorFlow 2.8.0과 Keras 2.8.0을 기반으로 설계되었으며, RF와 XGBoost 모델은 각각 Python 3.9.12 환경에서 Scikit-Learn 라이브러리의 RandomForestRegressor와 XGBRegressor를 이용하였다.

모델을 구축하는 과정에서 높은 예측 성능을 얻기 위해서 최적의 하이퍼파라미터를 설정하는 작업은 필수적이다(^{Oh et al., 2019}). ^{Müller and Guido (2016)}에 의하면 모델의 성능을 높이기 위해서 다수의 하이퍼파라미터를 활용하는 것은 오히려 모델을 복잡하게 만들 수 있다. 이에 본 연구에서는 모든 모델에서 4개의 하이퍼파라미터를 이용하였으며, Train data를 훈련시키는 과정에서 grid search를 통한 4-fold 교차검증으로 최적의 파라미터 조합을 선정하였다(Table 2). grid search는 연구자가 정의한 범위 내에서 파라미터 간의 다양한 조합을 설정한 뒤, K-fold 교차검증(K-fold cross validation)을 통해 예측 성능이 가장 우수한 파라미터 조합을 찾아내는 방법이다(^{Kim et al., 2018}). 이때 K-fold 교차검증은 학습 데이터를 k개의 그룹으로 나눈 뒤, 1개의 fold는 Validation data로 나머지 k-1개의 fold는 Train data로 사용하여 다음 과정을 k번 반복하여 학습을 진행하는 방법이다(^{Ha et al., 2017}).

학습을 통해 설정된 최적의 하이퍼파라미터 조합을 이용하여 각 모델을 학습하고 모델 예측 정확도를 비교하기 위하여 Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), R-Square (R²)와 Adjusted R-Square (Adj R²)를 비교한 결과를 Table 3에 제시하였으며, Train data와 Test data의 예측 정확도를 나타내는 그림을 Fig. 6에 제시하였다. MAE는 실제값과 예측값 차이 절대값의 평균, MSE는 실제값과 예측값 차이 제곱의 평균, RMSE는 MSE의 제곱근 값으로 세 값 모두 0에 가까울수록 예측 성능이 우수하다고 평가한다. R²는 실제값과 예측값의 선형회귀분석을 통해 계산하며 값이 1에 가까울수록 모델 재현성이 높다고 판단한다. 이때, R²는 독립변수의 개수가 증가할수록 값이 증가하는 경향이 있어 이를 보완하기 위해 Adj R²를 함께 평가하였다.

Fig. 6. Comparison of the actual Chl-a concentration and prediction results obtained using the three ML models: DNN, RF, and XGBoost

본 연구에서 Test data를 기준으로 예측 성능을 비교하였을 때, 예측 성능이 우수한 모델은 XGBoost, DNN, RF 순으로 각각 R²값이 0.9493, 0.9348, 0.9163으로 모든 모델에서 R²값이 0.9 이상인 매우 높은 Chl-a 농도 예측 정확도를 보였다. 특히, XGBoost의 경우, Train data 학습 성능은 0.9995, Test data는 0.9493으로 매우 높은 정확도를 보였다.

3.3 Feature Importance 분석

Feature Importance는 학습된 모델에서 목적 변수에 영향을 미치는 설명 변수의 중요도를 의미한다(^{Kim, Park et al., 2021}). 머신러닝의 학습 데이터는 다양한 설명 변수들의 특징값들로 이루어지는데, Feature Importance는 특정 특징값이 다른 값으로 대체될 경우 발생하는 오차의 크기를 가늠하여 계산되며, Feature Importance를 통해 입력된 특징값이 결과에 미치는 영향을 파악할 수 있어 머신러닝 모델 해석에 필수적인 과정이다(^{Lee, Oh et al., 2020}). DNN의 경우 TensorFlow의 GradientTape을 이용하였으며, RF와 XGBoost는 Scikit- Learn 라이브러리에서 제공되는 feature_importances_를 이용하여, 각 모델에 대한 Feature Importance 결과를 도출하였다(Fig. 7).

Fig. 7. Feature Importance of the prediction results obtained using the three ML models: DNN, RF, and XGBoost

앙상블 알고리즘 계열에 속하는 RF와 XGBoost는 상위 25%인 4순위까지의 Feature Importance가 BOD, COD, SS, PO₄-P로 동일하게 나타났으며, 수리 인자인 Elevation에 비해 Precipitaion, Average Temperature, Max Windspeed, Max Radiation의 기상 인자들은 상대적으로 낮은 중요도를 보였다. 각 모델을 통해 상위 50%에 도출된 인자는 DO, BOD, COD, SS, pH, PO₄-P, NH₃-N, NO₃-N, TN, TP, Water Temperature, Elevation이며 이 중 COD, SS, PO₄-P, NH₃-N는 모든 모델에서 Chl-a 농도 예측 시 높은 중요도를 가지는 주요 인자로 도출되었다.

3.4 Discussion

본 연구에서는 영산강 Chl-a 최적 관리 방안 마련을 위해 DNN, RF 및 XGboost 모델을 이용하여 모델별 예측 정확성을 확인하고, Feature Importance를 도출하여 Chl-a농도 예측 시 영향을 미치는 주요 요인을 파악하였다. R²을 통해 평가한 예측 정확성은 모든 모델에서 0.9 이상으로 매우 높은 것을 확인할 수 있었다. 특히, XGBoost의 경우 R² 값이 Train data에서 0.9995, Test data에서 0.9493으로 매우 높은 예측 정확성을 보였다. XGBoost는 부스팅 기법을 통해 높은 정확도를 얻을 수 있어, 머신러닝 기반의 알고리즘 중에서도 정확도가 높은 것으로 알려져 있다(^{Han et al., 2019}). 또한 최근에는 데이터 기반 온라인 경진대회 플랫폼인 Kaggle에서도 널리 사용되면서 모델의 우수성이 입증되고 있다(^{Lee, Oh et al., 2020}).

상관성 분석 결과 BOD, COD, SS, pH, PO₄-P는 Chl-a와 유의미한 상관관계를 가졌으며 BOD, COD, SS 및 PO₄-P는 모든 예측 모델에서의 Feature Importance 중요도가 높게 나타났다. 반면, pH는 DNN과 XGBoost을 이용한 Feature Importance 결과에서 상위 50% 이내 주요 인자에 포함되지 않았다. 또한, DO, NH₃-N, NO₃-N, TN, TP, Water Temperature 및 Elevation은 상관성 분석 결과에서는 유의미한 상관관계를 보이지 않았지만, NH₃-N은 모든 모델에서, TN과 TP는 DNN과 XGBoost에서, NO₃-N, Water Temperature와 DO, Elevation은 각각 DNN, RF, XGBoost에서 변수 중요도가 상위 50% 이내인 주요 인자로 도출되었다. 일반적으로 상관성 분석은 목적 변수와 설명 변수 간의 단순 선형 관계는 파악이 가능하나, 변수 간 상대적 영향도 및 비선형적 관계를 갖는 데이터에 대해서는 정확한 분석이 제한적이다(^{Jeong et al., 2021}). 이에 반해, 머신러닝을 통해 Feature Importance를 도출하면 선형 및 비선형 관계 모두에서 상대적 중요도를 파악할 수 있다(^{Jeong et al., 2021}; ^{Kim and Jun, 2019}). 따라서 Chl-a 농도 예측 시 영향을 미치는 주요 인자를 파악하기 위해서는 상관성 분석뿐 아니라 머신러닝을 이용한 Feature Importance 결과도 함께 분석할 필요가 있다.

낙동강 하류 물금지점에서 Chl-a와 여러 수질, 기상 및 수리 인자와 상관성 분석을 시행한 ^{Lim et al. (2015)}의 연구 결과에서 Water Temperature, pH, DO, BOD, COD, TN, NO₃-N 및 PO₄-P가 Chl-a와 유의미한 상관관계를 나타내었고, 기상 및 수리 인자는 상대적으로 유의미한 상관관계를 나타내지 않았다. 이는 본 연구에서 진행한 상관성 분석과 Feature Importance 결과와도 유사하였는데, 본 연구에서도 기상 및 수리 인자보다는 수질 인자가 Chl-a와 유의미한 상관성 및 모델 예측 시 높은 중요도를 나타냈다. 특히, PO₄-P는 Chl-a와 높은 음의 상관성을 보였으며, Feature Importance 결과에서도 Chl-a 농도 예측 시 중요한 변수로 도출되었다. 전 세계적으로 점오염원에 의한 인(Phosphorus)의 수층 유입은 감소하는 추세지만, 비점오염원에 의한 인의 유입은 여전히 조류 발생을 야기하는 원인이 되고 있다(^{Lee et al., 2013}). 본 연구 대상 지역인 승촌보 상류에는 영산강 유역에서 가장 큰 광주하수처리장이 위치하고 있으며, 총인 처리시설이 추가로 도입되어 PO₄-P와 TP 발생량이 감소하였다(^{Shin et al., 2015}). 그러나 ^{Kim et al. (2020)}의 연구에서는 승촌보 상류 지역에 다수의 배수문이 있어 도시, 산업지역 및 농경지로부터 오염물질이 강우 시 또는 항시 유입될 가능성이 존재한다고 하였다. 특히, 승촌보 상류에 위치한 신촌과 벽진 배수문은 타 배수문에 비해 오염부하량이 높았으며, 신촌의 경우 TP 부하량이 크고, 벽진의 경우 농경지의 비점오염물질이 혼재되어 있었다(^{Kim et al., 2020}). 일반적으로 수생태계에서 인의 작용은 매우 복잡한 것으로 알려져 있다. 낮은 분자량을 가진 불안정한 인 화합물은 콜로이드 형성을 위해 결합하거나, PO₄-P로 빠르게 가수분해 된다(^{Wetzel and Likens, 2013}). 이러한 인의 화학적 특징과 승촌보 상류의 높은 TP 부하량으로 인해 PO₄-P가 Chl-a 농도를 예측하는 데 주요 인자가 된 것으로 판단된다(^{Kim and Ahn, 2022}; ^{Shin et al., 2019}).

NH₃-N은 본 연구에서 Chl-a 농도 예측에 활용된 모든 모델에서 주요 인자로 도출되었다. 일반적으로 NH₃-N은 pH에 따라 NH₃, NH₄⁺ 형태로 나타나는데 자연수계에서 대부분 NH₄⁺ 형태로 존재하고 있다(^{Seo et al., 2018}). 이때 자정능력이 탁월한 하천의 경우에는 NH₃-N가 유입된 후 물리적, 화학적 및 생물학적 자정작용에 의해 질산성 질소로 산화되거나 생물체로 동화되어 농도가 낮아지지만, 지속적인 오염원 유입이 있는 경우 하천에서의 NH₃-N 농도가 높게 유지될 수 있다(^{Han et al., 2015}). 공공하수처리장의 방류수는 NH₃-N 농도를 증가시키는 주요 원인으로 알려져 있으며 본 연구 지점인 영산강 광산 지점의 상류에는 영산강에서 승촌보로 유입되는 광주하수처리장이 위치하고있다. 현재 광주 제1 하수처리장의 가동률은 100%에 육박하여 증설 요구가 지속되고 있으며, 방류수의 높은 총질소량을 낮추기 위해 질소 처리장치의 확충을 계획하고 있다(^{Kim et al., 2020}).

따라서 효과적인 Chl-a관리를 위해서는 Chl-a 농도 예측 시 주요 요인으로 도출된 PO₄-P와 NH₃-N에 대한 연구가 더욱 필요할 것으로 보이며, NH₃-N는 광주하수처리장의 총질소처리 장치 확충을 통해, PO₄-P는 점오염원인 광주하수처리장보단 영산강 상류 유역의 비점오염원 관리를 통해 오염물질 관리가 필요할 것으로 사료된다.