2.1 연구 대상 지역 및 데이터 수집

본 연구에서는 4대강(한강, 금강, 낙동강, 영산강) 본류 및 지류의 생물측정망 지점에서 2016년부터 2021년까지 수집된 수질, 수문, 기상 자료 및 저서성 대형무척추동물 자료를 활용하여 분석을 수행했다(Fig. 1). 한강, 금강, 낙동강, 영산강 본류의 유로 연장은 각각 480 km, 395 km, 510 km, 115 km이며, 유역 면적은 각각 22,298 km², 9,903 km², 23,690 km², 3,473 km²이다^{(Ji et al., 2012; Yu et al., 2018)}. 온난한 계절성 기후의 영향을 받으며, 연간 강수량의 54% (약 1,300 mm)가 여름철인 7월부터 9월까지 발생한다. 기온의 경우 가장 따뜻한 8월에는 월평균 기온이 22.1°C에서 31.0°C까지, 가장 춥고 건조한 1월에는 –7.2°C에서 3.6°C까지의 기온 범위가 나타난다^{(KMA, 2021)}.

Fig. 1. Study sites.

현재 저서성 대형무척추동물의 경우 생물측정망 지점별로 연 2회, 봄(4∼5월)과 가을(9∼10월)에 한 번씩 조사를 실시하고 있다^{(NIER, 2021)}. 생물측정망 지점의 경우 2008년부터 540개에서 시작해 2016년까지 3,035개로 지속적으로 증가해왔다(Table 1). 2015년까지는 모든 생물측정망 지점에서 매년 조사를 실시했고, 2016년부터 일반 지점 2,814지점은 3년(매년 약 938개)에 한 번씩, 상시지점(221개)은 매년 조사를 실시한다. 본 연구는 현행 BMI 등급체계가 시행된 시점인 2016년부터 2021년까지의 생물측정망 지점에서 수집된 총 13,912건의 표본 단위의 모니터링 자료를 수집했다.

모니터링 자료 내 저서성 대형무척추동물 자료의 경우, 출현 종의 개체 수와 종별 개체밀도(개체/m²) 자료를 기반으로 계산된 BMI 자료를 수집했다. 통계분석에 활용할 환경인자는 모니터링 자료 내에 저서성 대형무척추동물의 영향인자 중 가용한 변수로 총 8가지이다. 이는 pH, 5일 BOD(5-day biochemical oxygen demand, BOD5)(mg/L), 총질소(Total nitrogen, T-N)(mg/L), 총인(Total phosphorus, T-P)(mg/L)의 4가지 수질 변수와 수심(cm), 유속(cm/s) 같은 수문 변수, 그리고 대한민국 기상청에서 수집한 평균 기온(℃)과 일일 강수량(mm)의 2가지 기상 변수로 구성된다. 현행 저서성 대형무척추동물 건강성 평가 등급 기준 산정에 활용된 총부유물질(Total suspended solids, TSS)(mg/L)의 경우 생물측정망 자료 내에서의 결측 비율이 전체 자료의 94.94%로 매우 높아 활용 변수에서 제외했다.

MRT 기반의 다변량 통계분석에 앞서 데이터의 구조적 특성을 고려해 데이터 전처리를 수행했다. 먼저 각 변수에서 결측치가 포함된 행을 제거했다. 또한, 마할라노비스 거리(Mahalanobis distance)를 기반으로 이상치를 탐지하고 제거했다. 마할라노비스 거리는 다변량 데이터에서 각 변수의 분산을 고려하여 거리를 계산하는 방법으로, 이상치를 식별하는 데 적합한 방법론이다^{(De Maesschalck et al., 2000)}. 각 데이터셋의 변수들을 표준화(standardization)한 후, 데이터 표본($x$)과 평균 벡터($\mu$)의 차이를 공분산 행렬의 역행렬을 반영하여 마할라노비스 거리를 계산하였다(식 (1)). 마할라노비스 거리기반 이상치를 판별하기 위한 적절한 임계값을 설정하기 위해, 각 항목의 데이터 값이 평균($\mu$)에서 얼마나 떨어져 있는지를 표준편차 단위($\sigma$)로 나타내는 Z-score를 확인하였다(식 (2)). Z-score가 3 이상인 데이터를 이상치로 간주하였으며, 이는 데이터 값의 평균에서 3 표준편차 이상 떨어진 값으로 전체 데이터의 약 0.3%에 해당하는 매우 드문 데이터 값임을 의미한다^{(Anusha et al., 2019)}. 분석 결과, pH를 제외한 다른 항목들에서 약 100∼200개의 데이터를 이상치로 판별하였으며, 평균적으로 약 139개의 이상치가 확인되었다(Table 2). 이를 바탕으로 마할라노비스 거리 기반 이상치 판별을 위해 전체 데이터의 약 1%를 이상치로 제거하도록 임계값을 7.5로 설정하였다(Fig. 2). 최종적으로, 1,020개의 결측치와 104개의 이상치를 제거하였으며, 총 12,788개의 데이터 표본을 MRT 구축에 활용하였다.

Fig. 2. Distribution of Mahalanobis distance in the dataset.

$ x: Data \,point \\ \mu : Mean \,vector \,of \,data\\ \sum : Covariance $

$ Z :z-score\\ X : Data \,point \\ \mu : Mean \,of \,data\\ \sigma : Standard \,deviation \,of \,data $

2.2 Multivariate Regression Tree (MRT)

본 연구에서는 MRT를 활용하여 BMI 등급 체계와 환경변수 간의 관계를 분석하고, BMI의 최적 분할점을 체계적으로 탐색하였다. MRT는 기존의 Classification and Regression Tree (CART) 방법론을 확장한 분석 기법으로, 독립변수를 기준으로 데이터를 이진 트리 구조로 분할하며, 각 노드에서 종속변수의 편차를 최소화하는 방향으로 작동한다^{(Breiman et al., 1984)}. 트리 생성 과정은 뿌리 노드에서 시작하여 데이터를 반복적으로 이진 분할하며, greedy search를 통해 각 단계에서 종속변수 편차를 최소화하는 독립변수의 최적 분할 기준을 탐색한다. 이를 통해 데이터의 복잡한 구조와 종속변수 간의 관계를 반영한 최적화된 트리를 구축할 수 있다^{(De’Ath, 2002; Larsen and Speckman, 2004)}.

본 연구에서는 BMI를 독립변수로 설정하고, 환경변수를 종속변수로 사용하여 다변량 데이터의 구조와 변수 간의 상호작용을 반영한 트리를 구축하였다. 생성된 트리는 5개의 노드를 가지는 구조를 기본으로 설정하여 BMI 등급 분할점을 결정하였으며, 이를 통해 BMI와 환경변수 간의 다변량 관계를 정량적으로 평가하였다^{(Segal, 1992; Zhang, 1998)}. 추가적으로, 생성된 트리를 2∼4개의 노드를 가지는 구조로 축소하여 각 트리에서 도출된 분할점과 BMI 등급 간의 차이를 비교하였다. 이러한 분석은 BMI 등급 수의 조정이 분류 성능과 해석에 미치는 영향을 평가하기 위해 수행되었다.

2.3 후속 분석

MRT 기반 새로운 BMI 등급체계의 유용성을 확인하기 위해 후속 분석을 수행하였다. 먼저 변수중요도 지수(variable importance index)를 계산하여 MRT 기반 BMI 등급체계와 현행 BMI 등급체계에 대하여 환경변수들이 등급화에 미친 영향의 상대적 중요도를 분석하고, 주요 환경변수에 대한 BMI 등급 간 자료 분포 변별력을 평가하였다. 또한, 랜덤포레스트(random forest) 기반 인공지능 모델을 활용하여 MRT 기반 BMI 등급체계와 현행 BMI 등급체계를 각각 예측하고 재현성을 비교하였으며, MRT 기반 BMI 분류 시 등급 수 조정에 따른 예측력 변화를 추가적으로 분석하였다. 연구 흐름도는 Fig. 3과 같다.

Fig. 3. Research flow chart.

3.1 MRT 결과 분석

MRT를 수행한 결과, ‘매우 좋음’ 등급은 BMI가 91.45 이상일 때, ‘좋음’ 등급은 82.25 이상 91.45 미만, ‘보통’ 등급은 70.35 이상 82.25 미만, ‘나쁨’ 등급은 51.95 이상 70.35 미만, 마지막으로 ‘매우 나쁨’ 등급은 51.95 미만으로 분할되었다(Fig. 4). 이는 현행 BMI 건강성 평가 기준과 비교하였을 때, MRT에서 도출된 분할점이 현행 기준의 분할점보다 전반적으로 더 높은 경향이 나타났다(Fig. 5). 예를 들어, 현행 기준에서 ‘매우 좋음’ 등급은 BMI가 80 이상으로 설정되어 있으나, MRT 결과에서는 91.45 이상으로 상승하였다. 특히 ‘매우 좋음’, ‘좋음’ 같은 상위 등급의 경우, 분할점이 상향 조정됨과 동시에 해당 등급의 범위가 좁아져 현행 등급체계보다 엄격하게 설정되었다.

Fig. 4. Primary splits derived from multivariate regression tree models. The tree is constrained to a size of 5 nodes. A tree model is created based on the split point of the independent variable, BMI. Different colors represent the BMI grade associated with each final split node. N denotes the number of observations at each node.

Fig. 5. Comparison of the distribution of a) BMI grades based on MRT-based grade system b) BMI grades based on current grade system.

MRT 기반 BMI 등급체계에서의 데이터 분포를 분석한 결과, ‘매우 나쁨’ 등급에 해당하는 데이터가 3,756개로 가장 많았으며, 이어서 ‘나쁨’ 등급이 2,912개, ‘보통’ 등급이 2,980개, ‘좋음’ 등급이 2,443개, ‘매우 좋음’ 등급이 697개로 나타났다. 반면, 현행 BMI 등급체계에서는 ‘매우 좋음’ 등급이 3,774개로 가장 많았고, 이어서 ‘좋음’ 등급이 3,367개, ‘보통’ 등급이 2,157개, ‘나쁨’ 등급이 1,968개, ‘매우 나쁨’ 등급이 1,522개로 확인되었다. MRT 기반 등급체계에서는 하위 등급(‘매우 나쁨’ 및 ‘나쁨’)에서 데이터가 집중되었으나, 상위 등급(‘매우 좋음’ 및 ‘좋음’)의 비율이 대폭 감소하였다. 이는 MRT가 다양한 환경변수를 종합적으로 고려하여 BMI 등급을 분류함으로써 상위 등급에 보다 엄격한 기준을 적용했음을 의미한다. 반면, 현행 BMI 등급체계는 상위 등급 데이터의 비중이 과도하게 높아 실제 환경 조건을 과대평가했을 가능성이 있음을 확인했다. 이러한 결과는 환경변수와 BMI 간의 다변량 관계를 통합적으로 분석함으로써 실제로 환경이 우수한 지역을 명확히 구분하고, 보다 신뢰도 높은 평가 결과를 제공할 수 있음을 보여준다.

3.2 변수중요도 분석 결과

변수중요도 지수 산정을 통해 다양한 환경변수가 MRT 기반 BMI 등급화에 미치는 영향을 분석하였다(Fig. 6). 그 결과, 유속의 변수중요도 지수가 53.38%로 가장 높은 BMI 등급 구분에 가장 큰 영향을 미치는 것으로 나타났다. 이는 유속이 저서성 대형무척추동물의 서식지 안정성, 퇴적물 이동, 먹이 공급 등과 같은 서식 환경에 직접적으로 영향을 미치는 특성으로서 저서성 대형무척추동물 군집구조 및 조성에 영향을 미쳤기 때문인 것으로 사료된다(Table 3).

Fig. 6. Visualization of variable importance index results.

다음으로 기온이 12.92%로 높은 변수중요도를 보였는데, 이는 기온 변동에 따라 저서성 대형무척추동물의 서식지가 변화하면서 군집구조에 영향을 미치기 때문으로 판단된다. 이는 우리나라의 뚜렷한 사계절 특성을 고려할 때, 저서성 대형무척추동물의 건강성이 계절적 변동성을 보인다는 것을 의미한다.

BOD, T-N의 변수중요도 지수도 각각 12.36%, 12.08%로 기온과 유사한 중요도를 보여 저서성 대형무척추동물이 수질 변화에도 민감하게 반응함을 보여준다. 수심은 5.54%의 변수중요도를 기록하였고, T-P, 수심, 일일 강수량은 각각 2.68%, 1.93%, 0.77%로 낮은 변수중요도를 보여 하천 환경의 일부 특성만이 BMI 등급에 큰 영향을 미치는 것으로 나타났다.

이러한 결과는 저서성 대형무척추동물의 건강성이 하천의 물리적 특성과 화학적 특성 간의 상호작용에 의해 크게 영향을 받는다는 것을 의미한다. 특히 유속은 단일 변수로 가장 높은 중요도를 기록한 만큼, 수생태계 건강성 평가에서 우선적으로 고려해야 할 요소임을 보여준다. 또한, 기온도 저서성 대형무척추동물이 계절적 기온 변화에 민감하게 반응하도록 하며, 군집구조와 서식지 안정성에 중요한 영향을 미치는 것으로 나타났다. BOD, T-N과 같은 수질 변수도 생물학적 건강성과 밀접한 연관성을 가지며, 이를 포함한 종합적인 접근이 건강성 평가의 신뢰성을 높일 수 있다. 이러한 분석은 저서성 대형무척추동물의 건강성 평가를 통해 유속 변화와 계절적 기온 변동성, 주요 수질 변수의 영향을 구분해야 함을 시사한다.

3.3 변수들의 등급별 데이터 분포 비교 결과

현행 BMI 등급체계에서 환경변수의 등급별 분포 변화를 분석한 결과, BOD와 TP의 경우 E등급에서 A등급으로 갈수록 데이터의 평균값이 점진적으로 증가하는 경향이 관찰되어, 등급 간 분포 차이를 확인했다(Table 4). 그러나 유속의 경우 A, B, C 등급 간 변별력은 나타났으나, D와 E 등급 간 변별력은 떨어지는 것으로 나타났다. 등급 간 데이터 평균을 비교한 결과 D등급에서의 유속 평균은 18.2 m/s, E등급에서의 유속 평균은 18.4 m/s로, E등급에서 오히려 높은 유속이 확인되었다. 이는 현행 BMI 등급체계에서 유속과 같은 물리적 특성를 충분히 반영하지 못해 하위 등급 간의 구분이 뚜렷하지 못한 것으로 판단된다.

MRT 기반으로 재분류된 BMI 등급별 데이터 분포를 확인한 결과, BOD와 TP 뿐만 아니라 유속에 대해서도 등급 간 데이터 분포의 변별력이 나타났다(Fig. 7a). 유속의 등급 간 데이터 평균을 비교하였을 때, A등급은 48.6 cm/s, E등급은 18.4 cm/s로 현행 BMI의 등급체계와 비교하여 등급 간 변별력이 더 뛰어나다는 사실을 확인하였다(Table 4a). 두 번째로 높은 변수중요도를 보인 기온의 경우, 현행 BMI 등급체계에 비해 A등급과 E등급 간의 평균 차이가 크게 나타났다. 반면, BOD와 T-P에서는 A등급과 E등급 간 평균 차이가 작았으나, 평균값이 A > B > C > D > E 순으로 순차적으로 분포하였다. 이는 MRT 기반 BMI 등급체계가 현행 BMI 등급체계에서 보다 하천의 유속, 기온과 같은 물리적 특성을 효과적으로 반영했음을 보여준다.

Fig. 7. Boxplot comparison of a) the distribution of variables by grade based on MRT-based grade system and b) the distribution of variables by grade based on current grade system.

종합적으로, MRT 기반 BMI 등급의 재분류를 통해 현행 체계에서 등급 간 변별력이 낮았던 환경변수의 변별력이 상당히 향상되었음을 확인하였다. 이러한 결과는 MRT를 통해 BOD와 T-P뿐만 아니라, 저서성 대형무척추동물과 상관성이 높은 유속이나 기온 같은 변수에서도 BMI 등급체계를 적용하여 등급 간 데이터 분포의 변별력이 새롭게 드러났음을 의미한다. 이는 MRT가 다양한 환경변수 간의 상호작용을 통합적으로 분석함으로써, 환경변수와 BMI 간의 관계를 보다 정밀하게 반영한 결과이다. 특히, 다양한 환경요인에 영향을 받는 저서성 대형무척추동물의 건강성이 MRT 기반 체계에서 유속과 같은 물리적 특성을 등급 간 차이를 효과적으로 반영함으로써, 저서성 대형무척추동물의 건강성을 보다 효과적으로 평가할 수 있는 근거를 제공한다. 또한, MRT가 현행 BMI 등급체계가 시행된 2016년 이후의 생물측정망 데이터를 활용하여 환경적 변화를 반영하는 분류 기준을 설정했기 때문에, 현행 체계보다 현재 환경 상태를 더 정확히 진단할 수 있을 것으로 판단된다. 따라서, 본 연구의 접근 방식은 저서성 대형무척추동물의 건강성 평가에서 환경변수를 종합적으로 고려하고 변별력이 높은 기준을 제공하여, 더 정확하고 신뢰성 있는 평가를 가능케 할 수 있음을 시사한다.

3.4 저서성 대형무척추동물의 건강성 예측결과

MRT 기반 BMI 등급의 재분류가 저서성 대형무척추동물의 건강성 등급 분류모델의 예측 성능 향상에 얼마나 기여하였는지 확인하기 위해 분류모델링 후 성능을 비교하였다. 기존의 5등급 체계에서의 분류모델의 성능뿐만 아니라 2, 3, 4개 분할의 등급체계에서의 분류모델의 예측 성능도 확인하였다. 각 모델링의 조건에서 하이퍼파라미터 최적화를 수행한 결과는 Table 5와 같다.

현행 BMI 등급체계를 기반으로 한 분류모델의 성능을 확인해보았을 때, accuracy는 0.457이며, Cohen’s kappa 계수는 0.283으로 나타났다(Table 6). MRT 결과 기반 분류모델 성능의 경우 accuracy는 0.472이고, Cohen’s kappa 계수는 0.301이었다. 또한, 등급 수가 적어질수록 분류모델의 성능도 같이 증가하였는데, 마지막 2개 등급 분할에서의 accuracy는 0.780, Cohen’s kappa 계수는 0.560까지 향상되었다.

5개 등급의 랜덤포레스트 분류모델에서의 변수중요도를 MDA와 MDG 지표를 통해 평가한 결과, 앞선 변수중요도 지수의 결과와 유사하게 유속, BOD가 모든 지표에서 높은 중요도를 나타내며, 등급 간 차이를 설명하는 데 중요한 변수로 작용하였다(Fig. 8). 특히, 유속은 두 지표 모두에서 가장 높은 중요도를 기록하며 저서성 대형무척추동물의 서식지 특성을 결정짓는 주요 요인임을 확인하였다. 반면, 강수량, 기온과 같은 일부 변수는 낮은 중요도를 보여, 등급화 과정에서의 기여도가 상대적으로 적었다. 특히, 기온의 경우에는 MRT 결과에 대해서는 높은 변수중요도가 나타났으나 분류모델에서의 변수중요도는 낮았다. 이는 기온이 MRT 기반 군집 간 등급 차별성에는 큰 기여를 했으나, 분류모델에서는 다른 변수들과의 상관관계 및 모델이 예측에 활용한 비선형적 상호작용 속에서 상대적 중요도가 낮아진 것으로 판단된다. 이러한 결과는 MRT 기반 BMI 등급체계가 다양한 환경변수 간의 상호작용을 잘 반영하여 기존 체계보다 개선된 평가 기준으로 활용될 수 있음을 보여준다.

MRT 결과 기반의 분류모델이 다양한 환경변수들을 종합적으로 고려하여 더 나은 예측 성능을 보였다. 이는 모델 내 높은 변수중요도를 보인 유속과 같은 물리적 특성을 MRT를 통해 변수에 대한 등급 간 구분을 명확히 한 것이 성능 향상에 기여한 것으로 판단된다. 그러나 여전히 일부 등급 간 데이터 분포의 변별력이 부족한 환경변수가 존재하여, 분류모델의 성능 향상에 한계가 발생한다. 활용되는 변수의 증가로 모델이 복잡해지면서 과적합의 가능성이 높아져, 예측 성능의 향상이 제한된 것으로 판단된다. 하지만 이러한 결과는 현행 BMI 등급체계 기반의 분류모델도 일정 수준 이상의 예측 성능을 보이지만, 성능 향상의 여지가 충분히 존재한다는 것을 의미한다.

또한, 등급 수를 조정한 결과, 4개 등급에서는 accuracy가 0.517, Cohen’s kappa 계수가 0.347, 3개 등급에서는 accuracy가 0.648, Cohen’s kappa 계수가 0.391, 2개 등급에서는 accuracy가 0.780, Cohen’s kappa 계수가 0.560으로 나타나, 3개 이하의 분할에서 예측 성능이 크게 향상됨을 확인하였다. 분류모델의 특성상 분할의 수가 적어지면 예측 성능이 향상되는 것이 일반적이지만, 실제 하천 수생태계의 훼손 여부를 판단하기 위한 지표인 BMI에 다양한 등급체계를 적용해 개별 등급 간의 상대적 건강성 차이를 더 잘 반영함과 동시에 예측력을 향상시킬 수 있음을 시사한다.

Fig. 8. Comparison of Mean Decrease Accuracy and Mean Decrease Gini for variables in a) MRT-Based grade system and b) current grade system.