A Study on the Retrieval of River Turbidity Based on KOMPSAT-3/3A Images

KOMPSAT-3/3A 영상 기반 하천의 탁도 산출 연구

Abstract

Turbidity, the measure of the cloudiness of water, is used as an important index for water quality management. The turbidity can vary greatly in small river systems, which affects water quality in national rivers. Therefore, the generation of high-resolution spatial information on turbidity is very important. In this study, a turbidity retrieval model using the Korea Multi-Purpose Satellite-3 and -3A (KOMPSAT-3/3A) images was developed for high-resolution turbidity mapping of Han River system based on eXtreme Gradient Boosting (XGBoost) algorithm. To this end, the top of atmosphere (TOA) spectral reflectance was calculated from a total of 24 KOMPSAT-3/3A images and 150 Landsat-8 images. The Landsat-8 TOA spectral reflectance was cross-calibrated to the KOMPSAT-3/3A bands. The turbidity measured by the National Water Quality Monitoring Network was used as a reference dataset, and as input variables, the TOA spectral reflectance at the locations of in situ turbidity measurement, the spectral indices (the normalized difference vegetation index, normalized difference water index, and normalized difference turbidity index), and the Moderate Resolution Imaging Spectroradiometer (MODIS)-derived atmospheric products(the atmospheric optical thickness, water vapor, and ozone) were used. Furthermore, by analyzing the KOMPSAT-3/3A TOA spectral reflectance of different turbidities, a new spectral index, new normalized difference turbidity index (nNDTI), was proposed, and it was added as an input variable to the turbidity retrieval model. The XGBoost model showed excellent performance for the retrieval of turbidity with a root mean square error (RMSE) of 2.70 NTU and a normalized RMSE (NRMSE) of 14.70% compared to in situ turbidity, in which the nNDTI proposed in this study was used as the most important variable. The developed turbidity retrieval model was applied to the KOMPSAT-3/3A images to map high-resolution river turbidity, and it was possible to analyze the spatiotemporal variations of turbidity. Through this study, we could confirm that the KOMPSAT-3/3A images are very useful for retrieving high-resolution and accurate spatial information on the river turbidity.

탁도는 부유물질에 의한 빛의 산란 또는 흡수로 인한 수체의 흐림을 나타내는 수치로 수질 관리 분야에서 중요 지표로 활용되고 있다. 탁도는 소규모의 하천에서 변동성이 심할 수 있으며, 이는 국가하천의 수질에 직접적으로 영향을 준다. 따라서 고해상도의 탁도 공간정보 산출은 매우 중요하다. 이 연구에서는 Korea Multi-Purpose Satellite-3 및 -3A (KOMPSAT-3/3A) 영상으로부터 한강 수계 하천의 고해상도 탁도 매핑을 위한 eXtreme Gradient Boosting (XGBoost) 알고리즘 기반의 탁도 산출 모델을 개발하였다. 이를 위해 총 24장의 KOMPSAT-3/3A 영상과 150장의 Landsat-8 영상으로부터 계산된 대기 상단(Top Of Atmosphere, TOA) 반사율을 활용하였으며, Landsat-8 TOA 반사율은 KOMPSAT-3/3A의 관측 파장 대역에 적합하도록 교차검보정을 수행하였다. 국가수질자동관측망에서 측정된 탁도를 탁도 산출 모델의 참조자료로 사용하였고, 입력 변수로는 탁도가 실측된 위치에서의 TOA 분광반사율과 탁도 분석에 널리 이용되어 온 분광지수인 정규식생지수, 정규수분지수, 정규탁도지수, 그리고 Moderate Resolution Imaging Spectroradiometer (MODIS)의 대기 산출물(에어로졸 광학 두께, 수증기량, 오존)을 사용하였다. 또한 고탁도와 저탁도에 대한 KOMPSAT-3/3A TOA 분광반사율을 분석하여 탁도를 설명할 수 있는 새로운 정규탁도지수(new normalized difference turbidity index, nNDTI)를 제안하였고, 이를 탁도 산출 모델에 입력 변수로 추가하였다. XGBoost 기반 탁도 산출 모델은 현장관측 탁도와 비교하여 2.70 NTU의 평균 제곱근 오차(root mean square error, RMSE) 및 14.70%의 정규화된 RMSE(normalized RMSE)를 가지는 탁도를 예측하여 우수한 성능을 보였으며, 이 연구에서 새롭게 제안한 nNDTI가 탁도 산출에 있어 가장 중요한 변수로 사용되었다. 개발된 탁도 산출 모델을 KOMPSAT-3/3A 영상에 적용하여 하천 탁도를 고해상도로 매핑하였으며, 탁도의 시공간적 변동에 대한 분석이 가능하였다. 이 연구를 통하여 고해상도의 정확한 탁도 공간정보 산출에 KOMPSAT-3/3A 영상이 매우 유용함을 확인할 수 있었다.

1. 서론

탁도(turbidity)는 수중의 생물 및 비생물 부유 입자에 의해 발생하는 수체의 흐림 정도를 정량적으로 나타내는 수치로, 수질 관리를 위한 모니터링에 널리 이용되고 있는 지표이다(Guttler et al., 2013; Warren et al., 2021). 우리나라에서는 1990년대 후반부터 부유물질의 농도와 탁수 발생 기간의 증가로 인한 하천의 고탁도 현상이 빈번하게 발생하여 주요한 환경 이슈로 대두되기 시작하였다(Kim and Jung, 2007). 하천의 탁도 증가는 수중 생물의 개체 수를 감소시키는 등 수중 생태계에 큰 영향을 미칠 뿐만 아니라, 하천 및 호소수의 오염도를 높여 공공수의 활용에 제약으로 작용한다. 따라서 수체의 탁도 모니터링을 통해 수질 현황과 오염도를 관리할 필요가 있다. 우리나라는 환경정책기본법과 물환경보전법에 근거하여 주요 하천 및 호소에 대해 수질측정소(국가수질자동측정망)를 설치하여 탁도를 측정하고 있다. 국가수질자동측정망에서 측정된 탁도는 1시간 주기로 제공되므로 매우 높은 시간해상도를 가지지만, 측정소 위치에 대한 정보만을 제공하는 점(point) 자료이므로 수계 공간에 대한 탁도 분포를 설명하기엔 한계가 있다. 탁도는 수체의 흐름 특성과 수심, 지류 하천의 유입 등 다양한 요인에 영향을 받아 시공간적으로 변화할수 있기 때문에 탁도에 대한 고해상도의 공간정보 산출이 요구된다.

탁도는 수중의 다양한 부유 입자에 의한 빛의 산란 또는 흡수로 결정된다. 따라서 인공위성 다중분광 영상이 탁도의 공간정보 산출에 효과적으로 이용될 수 있다. 영상 자료가 무상으로 제공되어 접근성이 높은 Landsat 시리즈와 Sentinel-2 위성 영상의 분광 반사율을 이용한 수체의 탁도 산출 및 모니터링 연구가 다수 수행된 바 있으며(Bustamante et al., 2009; Güttler et al., 2013; Braga et al., 2017; Quang et al., 2017; Garg et al., 2020; Peterson et al., 2020; Ma et al., 2021; Warren et al., 2021; El-Alem et al., 2022), 이 중 일부 연구들에서는 탁도 산출 모델 개발을 위해 다양한 종류의 기계학습 기반 회귀분석이 사용되었다. 이 연구들은 가시광선과 근적외선 파장 대역의 반사율이 탁도에 민감하게 반응하며, 기계학습의 적용을 통해 비교적 정확한 탁도 산출이 가능함을 설명하였다. Landsat 시리즈와 Sentinel-2 위성의 다중분광 영상은 각각 30 m와 10–20 m의 공간해상도를 가진다. 이는 폭이 수 백 m 이상으로 넓은 하천이나 대형 호소의 탁도 모니터링에 유용할 수 있으나, 폭이 좁은 소규모 하천에는 적용이 어렵다. 일반적으로 소규모 하천은 수심이 얕고 수계의 유동 특성 변화가 심하다. 따라서 동일한 양의 부유물질이 존재하더라도 소규모 하천에서는 대규모 하천 및 호소에 비해 탁도의 시공간 변동성이 높을 수 있다. 따라서 소규모 하천에 대한 탁도 공간정보 산출은 매우 중요하며, 고해상도의 다중분광 위성 영상이 필수적으로 요구된다.

우리나라의 국가 고유 자원인 Korea Multi-Purpose Satellite-3 (KOMPSAT-3)는 2.8 m의 공간해상도로 가시광선(청색, 녹색, 적색) 및 근적외선 파장 대역의 영상을 제공하며, KOMPSAT-3A는 동일한 다중분광 파장 대역에서 2.2 m 해상도의 영상을 제공한다. KOMPSAT-3/3A 영상의 공간해상도는 소규모 수계의 모니터링에 적용 가능할 정도로 고해상도이다. 그러나 현재까지 KOMPSAT-3/3A 영상을 통해 수체의 탁도를 산출하거나 분석한 연구 사례는 거의 없다. 기존의 다중분광 위성영상을 이용한 연구들에서는 가시광선과 근적외선 파장 대역에서 수체의 반사율이 탁도에 따라 다르게 나타난다고 보고하였으며(Lacaux et al., 2007; Wu et al., 2014; Garg et al., 2017), 이는 KOMPSAT-3/3A가 Landsat 시리즈나 Sentinel-2 위성에 비해 적은 수의 분광 밴드를 가짐에도 탁도 산출 연구에 충분히 활용될 수 있음을 시사한다.

이 연구에서는 한강 유역 하천에 대해 KOMPSAT-3/3A 영상에 적용 가능한 탁도 산출 모델을 개발하였다. KOMPSAT-3/3A의 다중분광 파장 대역에 대한 반사율과 기계학습 기반 회귀분석 기법을 활용하여 탁도 산출 모델을 개발하고 성능을 평가하였으며, 소규모 하천에 대해 분석이 가능한 고해상도의 탁도 공간정보를 생산하였다. 그리고 탁도의 시공간적 변화와 그 원인에 대한 분석도 수행하였다.

2. 연구지역 및 자료

1) 연구지역

이 연구에서는 한강 유역의 하천에 대해 KOMPSAT-3/3A 기반의 탁도 산출 연구를 수행하였다. 한강은 우리나라의 최대 하천으로 크게 남한강, 북한강, 서울을 관통하는 한강으로 구분된다. 한강유역환경청에 의하면, 한강의 유역 면적은 약 26,000 km²로 국토의 약 23%를 차지하고 있으며, 유로연장은 514.8 km, 유역 내 하천의 총 길이는 약 7,250 km에 달한다(http://www.me.go.kr/hg/web/index.do?menuId=1265). 한강은 인구가 밀집된 서울과 경기도를 가로지르고 있어 수계의 탁도는 수도권의 먹는 물 및 농업용수 공급에 영향을 줄 수 있다. 북한강과 남한강 상류는 상대적으로 인구밀도가 낮지만 강우 시 부유토사 유출 등으로 인한 고탁도 현상이 자주 발생하고 있다(Choi, 2006). 한강 본류와 남한강, 북한강이 만나는 지점 부근에 형성된 하천형 호소인 팔당호에서는 식물 플랑크톤의 증감으로 인한 탁도의 변화가 보고된 바 있다(Kim et al., 2002). 한강 수계 하천의 탁도는 시공간적으로 높은 변동성을 가지며, 변화 기작이 다양하여 고해상도 위성 영상을 이용한 공간정보 산출이 요구된다고 할 수 있다.

2) 연구자료

(1) KOMPSAT-3/3A 영상

한강 수계 내의 하천에 대해 2014년 12월부터 2020년 12월 사이에 획득된 KOMPSAT-3 Advanced Earth Imaging Sensor System (AEISS) 영상 20장과 KOMPSAT-3A AEISS-A 영상 4장을 사용하였다(Table 1). 모든 KOMPSAT-3/3A 영상은 운량 조건 30% 이하에서 획득되어 Level 1G로 배포된 영상으로, 대한원격탐사학회 주관 2022 위성영상 활용 새싹연구 지원사업을 통해 구축되었다. 현재 한강 수계에서 운용 중인 하천의 수질자동측정소들은 서로 수–수십 km 떨어져 있다. KOMPSAT-3와 KOMPSAT-3A는 각각 ~13 km와 ~16 km의 관측 폭(swath)을 가지며, 이 연구에서 사용된 영상들에는 각각 한 개의 수질자동측정소만이 포함된다 (Table 1). KOMPSAT-3와 KOMPSAT-3A는 모두 서브미터 급 공간해상도의 전정색(panchromatic band) 영상을 제공한다. 그러나 이 연구에서는 가시광선 및 근적외선 파장 대역의 반사율 정보에 기초하여 탁도 산출 모델을 개발하고자 하였으며, 이에 따라 전정색 영상은 연구에 사용하지 않았다.

Table 1. KOMPSAT-3 and KOMPSAT-3A images used in this study

(2) Landsat-8 OLI 영상

KOMPSAT-3/3A 영상에는 한 개의 수질자동관측소가 포함되므로 탁도 산출 모델 개발을 위한 기계학습에 활용할 수 있는 샘플은 24개에 불과하다. 이는 기계학습을 수행하기에 매우 부족한 개수이다. 샘플 개수의 증대를 위해 KOMPSAT-3/3A와 같이 청색, 녹색, 적색, 근적외선 파장 대역을 모두 가지고 있는 Landsat-8 Operational Land Imager (OLI) 영상 150장을 구축하였다. Landsat-8 OLI 영상은 한강 유역이 포함되는 Path 115와 116, Row 34에서 2013년부터 2021년까지 운량 조건 30% 이하에 대해 촬영된 Level 1 Precision and Terrain Correction (L1TP) 영상으로 구축하였다. Landsat-8 OLI 영상은 30 m의 공간해상도를 가지므로 최소 150 m (5 픽셀) 이상의 하천 폭을 가지는 곳에 수질측정소가 설치되어야 해당 측정소 위치에서 수체에 대한 분광 반사율 정보를 취득할 수 있다. 이와 같은 조건에 해당하는 수질자동측정소는 모두 7개(여주, 능서, 흥천, 가평, 양평, 구리, 의암 측정소)이다(Table 2). 이 연구에서는 Landsat-8 OLI 다중분광 영상에서 KOMPSAT-3/3A와 파장 대역이 유사한 가시광선의 3개 밴드 영상과 근적외선 밴드의 영상만을 사용하였다.

Table 2. Landsat-8 OLI images used in this study

(3) Terra/Aqua MODIS 산출물

다중분광 영상으로부터 탁도와 같은 수질 분석을 위해서는 대기보정을 통해 정확한 지표 반사율을 산출하여야 한다. 이 연구에서 구축된 KOMPSAT-3/3A 및 Landsat-8 영상은 모두 대기보정이 수행되지 않은 자료이다. Landsat-8의 경우 미국지질조사소(United States Geological Survey)에서 대기보정이 수행된 지표 반사율 영상을 배포하고는 있으나 수체 표면에서의 작은 신호대 잡음비와 불완전한 에어로졸 영향 제거로 인해 수체에서 부정확한 정보를 제공할 수 있다(Xu et al., 2020). KOMPSAT-3/3A 영상 또한 대기보정이 적용된 자료를 제공받을 수 있으나 Landsat-8 영상과 같은 이유로 대기보정 결과에 부정확성이 내재될 수 있다. 이 연구에서는 KOMPSAT-3/3A와 Landsat-8 영상에 대기보정을 적용하지 않는 대신 Terra/Aqua 위성의 Moderate Resolution Imaging Spectroradiometer (MODIS) 관측 자료를 통해 제작된 에어로졸 광학두께(Aerosol Optical Thickness, AOT), 수증기량(Water Vapor, WV), 오존(Ozone, OZ) 산출물을 탁도 산출 모델 개발을 위한 기계학습의 입력 변수로 활용하여 KOMPSAT-3/3A 및 Landsat-8 영상에 포함된 대기 영향을 고려하고자 하였다. AOT와 OZ은 MOD/MYD09CMA product, OZ은 MOD/MYD09CMG product에서 제공된다. 연구에 사용된 MODIS 산출물은 하루 동안의 MODIS 영상들을 일간 합성하여 제작된 daily product로 0.05°의 공간해상도를 가지며, 이는 Landsat-8 영상의 대기보정에 널리 활용되고 있다(Vermote et al., 2016; Claverie et al., 2018; Kuhn et al., 2019). KOMPSAT-3/3A와 Landsat-8 위성의 촬영 시각이 각각 지방시(local time)로 오후 2시 30분경, 오전 11시경임을 고려하여 KOMPSAT-3/3A 영상에 대해서는 촬영 시각이 지방시 오후 1시 30분경인 Aqua MODIS, Landsat-8 영상에 대해서는 촬영 시각이 지방시 오전 10시 30분경인 Terra MODIS의 산출물을 활용하였다. 각각의 MODIS 산출물은 30 m 격자 크기로 내삽하여 사용하였다.

(4) 탁도 현장관측 자료

환경부에서 운영하고 국립환경과학원에서 관리하는 국가수질자동측정망에서 측정된 탁도(단위: Nepthelometric Turbidity Unit, NTU)를 탁도 산출 모델의 참조자료(reference dataset)로 사용하였다. 수질자동측정소에서 측정된 탁도는 물환경정보시스템(https://water.nier.go.kr/)을 통해 1시간 간격으로 제공되고 있으며, 연구에 사용된 위성 영상의 촬영 시각과 가장 가까운 시간대의 자료로 구축하였다. KOMPSAT-3/3A 영상에 포함된 수질자동측정소는 총 10개소이며, Landsat-8 영상에서 수체로 식별 가능한 위치에 설치된 측정소는 총 7개소로 각각의 위치는 Fig. 1에 표시한 것과 같다. KOMPSAT-3/3A 영상에 대해 구축된 탁도 관측값은 총 24개로 0.4–18.8 NTU의 범위를 가지며, Landsat-8 영상에 대해 구축된 관측값은 총 704개로 0.3–149.7 NTU의 범위를 가진다. Landsat-8 영상에서는 구름 또는 구름의 그림자로 인해 수질측정소 주변의 수체 식별이 어려운 경우 해당 위치에서의 탁도 현장관측 자료는 구축하지 않았다.

Fig. 1. Landsat-8 OLI true color image obtained on March 23, 2015. The red dots indicate the locations of water quality monitoring stations.

3. 연구방법

1) 탁도 산출 모델에 활용하기 위한 TOA 분광반사율 계산

KOMPSAT-3/3A 다중분광 영상과 Landsat-8 다중분광 영상으로부터 대기 상단(Top Of Atmosphere, TOA) 반사율을 계산하였다. KOMPSAT-3/3A 영상에 대해 식(1)을 이용하여 화소 값(digital number, DN)을 파장(λ)별 TOA 복사휘도(L_λ)로 변환하였다(SIIS, 2018).

L_λ = Gain × (DN) + Offset      (1)

식(1)에서 Gain과 Offset은 영상의 메타데이터에서 제공되는 밴드 별 변환계수이다. L_λ를 계산한 후 식(2)를 이용하여 KOMPSAT-3/3A 영상의 파장 별 TOA 반사율 (ρ_{λ, K})을 계산하였다.

\(\begin{aligned}\rho_{\lambda, K}=\frac{\pi L_{\lambda} d^{2}}{E S U N_{\lambda} \cos \theta_{s}}\\\end{aligned}\)       (2)

식(2)에서 d는 태양-지구 거리, ESUN_λ는 파장 별 외기권에서의 태양복사조도, θ_s는 태양 천정각을 의미한다. Landsat-8 영상은 식(3)을 이용하여 화소 값을 TOA 반사율(ρ_{λ, L})로 변환하였다(USGS, 2019).

\(\begin{aligned}\rho_{\lambda, L}=\frac{\left(M_{p} \times(D N)+A_{p}\right)}{\cos \theta_{s}}\\\end{aligned}\)       (3)

식(3)에서 M_p와 A_p는 밴드 별 변환계수로 영상의 메타데이터에서 제공된다.

KOMPSAT-3/3A AEISS/AEISS-A와 Landsat-8 OLI의 가시광선 및 근적외선 파장 대역은 서로 차이가 있다(Table 1 and 2). 따라서 두 이종 위성의 영상을 탁도 산출 연구에 상호 융합하여 사용하기 위해서는 교차검보정(cross calibration) 과정이 필요하다. KOMPSAT-3/3A와 Landsat-8 사이의 교차검보정에는 분광 밴드 조정치 (Spectral Band Adjustment Factor, SBAF)를 활용하는 방안이 활발히 논의되어 오고 있다(Ahn et al., 2020; Jin and Choi, 2021). 이 연구에서는 KOMPSAT-3/3A 영상 개수의 한계를 보완하면서 탁도 산출을 위한 기계학습에 충분한 데이터를 제공할 목적으로 분광 밴드 별로 KOMPSAT-3/3A와 Landsat-8 영상의 TOA 반사율 간에 교차검보정을 수행하였다. Jin and Choi (2021)에서 제시된 SBAF와 식(4)를 이용하여 KOMPSAT-3/3A의 분광 대역에 교차검보정된 Landsat-8 파장 별 TOA 반사율(ρ_{λ, L*})을 산출하였다.

\(\begin{aligned}\rho_{\lambda, L^{*}}=\frac{\rho_{\lambda, L}}{S B A F}\\\end{aligned}\)       (4)

두 이종 위성 영상 간의 교차검보정 결과는 동일한 날짜에 획득된 KOMPSAT-3 영상과 Landsat-8 영상 중 구름이 없고 AOT 및 WV가 작은 조건에서 획득된 영상을 선정한 후, 다양한 토지피복에서 추출한 ρ_{λ, K} 및 ρ_{λ, L*} 사이의 상관성을 분석하여 평가하였다. 이때두 이종 위성 영상 간의 공간해상도를 고려하여 최소 150 m × 150 m 이상의 영역에서 균일한 특성을 보이는 피복에 대해 ρ_{λ, K} 및 ρ_{λ, L*}을 추출하여 분석하였다.

2) 탁도 산출을 위한 분광지수 계산 및 제안

정규식생지수(Normalized DifferenceVegetation Index, NDVI)와 정규수분지수(Normalized Difference Water Index, NDWI)는 탁도를 분석하는데 이용될 수 있는 대표적인 분광지수이다(Hou et al., 2017). 부유물질에 의해 높아지는 탁도는 NDVI 값의 증가와 NDWI 값의 감소를 야기할 수 있다. 이 연구에서는 KOMPSAT-3/3A 및 교차검보정 된 Landsat-8 TOA 반사율로부터 식(5)와 (6)을 통해 NDVI와 NDWI를 계산하고 탁도 산출 모델의 입력 변수로 고려하였다.

\(\begin{aligned}NDVI=\frac{\rho_{NIR}-\rho_{Red}}{\rho_{NIR}+\rho_{NIR}}\end{aligned}\)       (5)

\(\begin{aligned}NDWI=\frac{\rho_{Green}-\rho_{NIR}}{\rho_{Green}+\rho_{NIR}}\end{aligned}\)       (6)

여기서 ρ_NIR, ρ_Red, ρ_Green는 각각 근적외선, 적색, 녹색 파장 대역의 반사율을 의미한다.

Lacaux et al. (2007)은 호소수의 탁도 분석을 위해 ρ_Red과 ρ_Green을 이용하여 식(7)과 같은 정규탁도지수 (Normalized Difference Turbidity Index, NDTI)를 제안하였다.

\(\begin{aligned}NDTI=\frac{\rho_{Red}-\rho_{Green}}{\rho_{Red}+\rho_{Green}}\end{aligned}\)       (7)

높은 탁도의 호소수는 낮은 ρ_Green와 상대적으로 높은 ρ_Red에 의해 양의 NDTI 값을 가지며, 탁도가 높아질수록 그 값이 커진다고 보고되고 있다(Lacaux et al., 2007; Garg et al., 2017). 식(7)을 통해 계산된 NDTI를 NDVI 및 NDWI와 함께 탁도 산출 모델 개발에 활용하였다.

탁수의 분광반사율은 호소와 하천에서 다르게 나타날 수 있으며, 탁도를 유발시키는 다양한 원인으로 인해 지역에 따라서도 상이할 수 있다. 따라서 Lacaux et al.(2007)에서 제안된 NDTI는 한강 수계 하천의 탁도를 설명하는데 적합하지 않을 가능성이 있다. 이 연구에서는 AOT 0.15 이하, 수증기량 2.0 g/cm² 이하로 분광반사율에 대한 대기의 영향이 작은 KOMPSAT-3/3A 영상에서 고탁도와 저탁도에 따른 분광반사율을 분석하고 한강 수계의 탁도를 설명할 수 있는 새로운 분광지수의 제안을 검토하였다.

3) eXtreme Gradient Boosting 기반 탁도 산출 모델 개발

이 연구에서는 eXtreme Gradient Boosting (XGBoost) 기법을 이용하여 KOMPSAT-3/3A 영상을 위한 탁도 산출 모델을 개발하였다. 수질자동측정소에서 측정된 탁도를 XGBoost 모델의 종속 변수로 사용하였다. 수질측정소 위치에서 추출된 교차검보정된 Landsat-8 TOA 분광반사율과 분광지수들, Terra MODIS의 AOT, WV, OZ을 모델의 훈련 자료(training dataset)로 사용하였고, KOMPSAT-3/3A TOA 분광반사율, 분광지수들 및 Aqua MODIS의 대기 산출물들은 모델의 테스트 자료 (test dataset)로 사용하였다(Fig. 2). TOA 분광반사율과 분광지수들은 측정소가 위치한 픽셀을 중심으로 3 × 3 픽셀 윈도우에 대한 평균값으로 산출하여 사용하였다. 훈련 자료와 테스트 자료는 각각 704개와 24개의 샘플로 구성되었다.

Fig. 2. An overall flowchart of this study.

XGBoost 기법은 의사결정나무 기반 알고리즘의 Boosting 기법을 개선한 Gradient Boosting 기법 중 하나로, 모델의 과적합을 줄이고 분산 및 병렬 연산을 통해 기존 의사결정나무 기반 알고리즘보다 빠르고 효율적인 학습을 가능하게 한다(Chen and Guestrin, 2016). Boosting 기법은 먼저 훈련자료의 일부를 임의로 추출하여 약한 예측 모델을 만들고 정확도를 평가한 후 새로운 약한 예측 모델을 만드는데, 이 때 이전 모델에서 오차가 큰 샘플들에 가중치를 증가시켜 사용한다. 이와 같은 과정을 반복하여 여러 예측 모델을 앙상블함으로써 최종 예측 모델을 구현한다. Gradient Boosting 기법은 경사하강법(gradient descent)을 이용하여 모델의 정확도가 높아지는 새로운 모델을 구성하여 추가하는 Boosting 기법이며, 이를 분산 및 병렬 연산으로 구현되도록 한 것이 XGBoost 기법이다. 이 연구에서는 R 소프트웨어의 xgboost 패키지를 이용하였다. XGBoost 기반 탁도 산출 모델의 성능은 모델에서 예측된 탁도와 현장관측 탁도 사이의 root mean square error (RMSE)와 normalized RMSE (NRMSE)로 평가하였다.

KOMPSAT-3/3A 영상으로부터 가장 정확한 탁도를 산출할 수 있는 최적의 입력 변수 조합을 파악하기 위해 총 4개 case에 대한 탁도 산출 모델을 개발하고 성능을 평가하였다. 첫 번째 case는 TOA 분광반사율만을 입력 변수로 사용하는 경우, 두 번째는 TOA 분광반사율과 MODIS 기반 대기 정보(AOT, WV, OZ)를 사용하는 경우, 세 번째는 탁도 분석에 사용되어 온 분광지수들(NDVI, NDWI, NDTI)을 추가하여 사용하는 경우, 마지막으로 네 번째 case는 이 연구에서 제안 가능한 새로운 분광지수까지 추가하여 사용하는 경우이다. XGBoost 기법은 사용되는 hyperparameter에 따라 다른 모델 성능을 보일 수 있다. 이 연구에서는 각각의 case에 대해 XGBoost의 대표적인 hyperparameter인 subsample(샘플링에 사용되는 비율), gamma(leaf node의 분할을 결정할 최소 손실 감소 값), max_depth(의사결정나무의 최대 깊이), min_child_weight(child node에 필요한 관측치 가중치의 최소 합), learning_rate(학습률)의 최적화를 위해 그리드 탐색(grid search)을 수행하였다. 그리드 탐색에 사용된 각 hyperparameter의 후보 값들은 Table 3과 같다. 최적의 hyperparameter는 훈련 자료의 5-fold cross validation을 통한 모델 검증에서 가장 작은 NRMSE를 보이는 값으로 선택하였다.

Table 3. Hyperparameters of the XGBoost turbidity retrieval models

가장 정확한 탁도를 산출하는 모델에 대해 XGBoost 의 변수중요도 지표 중 하나인 Gain을 이용하여 입력 변수들의 중요도를 해석하였다. Gain은 해당 변수를 사용했을 때 얻는 훈련 성능 상의 이득으로, 그 변수가 모델 예측에 어느 정도 영향을 미쳤는지를 나타내는 지표이다.

4. 연구결과 및 토의

1) KOMPSAT-3/3A와 Landsat-8 영상 간 교차검보정 결과 평가

구름이 전혀 없고 AOT와 OZ이 각각 0.2, 1.2 g/cm² 이하로 대기의 영향이 비교적 적은 조건에서 획득된 2015년 3월 30일의 화천 지역 KOMPSAT-3 영상과 Landsat-8 영상을 이용하여 두 이종 위성 영상 사이의 교차검보정 결과를 평가하였다. 공간해상도가 다른 두 영상에서 공통적으로 식별 가능한 수체, 나지, 초지, 산림, 대형 건물의 지붕을 대상으로 각 영상에서 10개의 TOA 분광반사율 값을 산출한 후 상관분석을 수행하였다(Fig. 3). 분광반사율은 각각의 피복에서 3 × 3 픽셀 윈도우에 대한 평균값으로 계산하였다. KOMPSAT-3 TOA 반사율과 SBAF가 적용된 Landsat-8 TOA 반사율은 4개 파장 대역에서 모두 0.97 이상의 높은 상관성을 보이며, RMSE는 2% 이하로 매우 작았다. 이와 같은 결과는 교차검보정 결과가 신뢰할 만한 수준으로 수행되었음을 보이며, SBAF가 적용된 Landsat-8 TOA 분광반사율이 KOMPSAT-3/3A를 위한 탁도 산출 모델 개발에 활용 될 수 있음을 설명한다. 근적외선 밴드에서 교차검보정에 의한 RMSE는 1.8%로 가시광선 밴드들에 비해 2–3배 컸는데, 이는 대기 중 수증기에 영향을 많이 받는 근적외선 밴드의 특성에 기인하는 것으로 추정된다(Jin and Choi, 2021).

Fig. 3. Scatter plots between KOMPSAT-3/3A TOA reflectance and the Landsat-8 TOA reflectance calibrated with SBAF for the (a) blue, (b) green, (c) red, and (d) near infrared bands.

2) New Normalized Difference Turbidity Index (nNDTI)

Fig. 4는 대기의 영향이 적은 조건에서 획득된 KOMPSAT-3/3A 영상에서 수질측정소 위치에 대해 Lacaux et al. (2007)이 제안한 NDTI와 실측 탁도 사이의 산점도를 나타낸 것이다. NDTI는 현장관측 탁도를 전혀 설명하지 못하는 것으로 나타난다(R²=0.081). 이는 NDTI가 하천이 아닌 대규모의 호소수를 대상으로 제안되었으며, 해당 호소수는 한강 수계 하천과는 다른 탁도 발생 기작을 가질 수 있기 때문인 것으로 추정된다. Fig. 5는 18.8 NTU의 고탁도(2018년 5월 23일, 경안천 측정소)와 1.1 NTU의 저탁도(2015년 10월 29일, 원주 측정소)가 측정된 위치에서 산출한 KOMPSAT-3/3A TOA 분광반사율을 보여준다. 청색 밴드의 TOA 반사율은 탁도에 민감하지 않은 반면, 나머지 3개 밴드의 반사율은 저탁도보다 고탁도에서 높게 나타난다. 이 중에서 녹색밴드의 TOA 반사율이 탁도에 따라 가장 크게 변화하는 것으로 나타났다. 이러한 경향성을 통해 이 연구에서는 식(8)과 같이 새로운 정규탁도지수(new Normalized Difference Turbidity Index, nNDTI)를 제안하였다.

Fig. 4. Scatter plot between turbidity and NDTI calculated from the KOMPSAT-3/3A TOA spectral reflectance.

Fig. 5. The KOMPSAT-3/3A TOA spectral reflectance for low and high turbidities.

\(\begin{aligned}nNDVI=\frac{\rho_{Green}-\rho_{Blue}}{\rho_{Green}+\rho_{Blue}}\end{aligned}\)       (8)

nNDTI는 저탁수에서 음의 값을 가지며, 탁도가 높아질수록 증가하는 경향을 보이게 된다.

Fig. 6은 Fig. 4의 현장관측 탁도와 이에 상응하는 nNDTI 사이의 산점도를 나타낸 것이다. 이 연구에서 제안한 nNDTI는 탁도가 높아질수록 증가하는 경향을 보이며, 기존의 NDTI보다 탁도를 더욱 잘 설명할 수 있음이 확인되었다(R²=0.649). 이 연구에서는 nNDTI를 새로운 분광지수로 탁도 산출 모델의 입력 변수에 추가하여 활용하였다.

Fig. 6. Scatter plot between turbidity and nNDTI calculated from the KOMPSAT-3/3A TOA spectral reflectance.

3) 최적의 탁도 산출 모델 개발 및 평가

총 4개의 입력 변수 조합 case에 대해 최적의 hyperparameter로 선택된 값들은 Table 3에 나열되어 있다. Table 4는 각각의 case 별로 최적의 hyperparameter를 이용하여 모델을 훈련시킨 후 테스트 자료에 대해 예측한 결과를 보여준다. TOA 분광반사율만을 입력 변수로 사용했을 때 모델의 성능이 가장 좋지 않았으며(NRMSE=51.13%, RMSE=9.41 NTU), nNDTI를 포함하여 분광지수들(NDVI, NDWI, NDTI), TOA 분광반사율, MODIS 대기 산출물(AOT, WV, OZ)을 모두 입력 변수로 사용한 모델이 가장 높은 성능(NRMSE=14.70%, RMSE=2.70 NTU)을 보였다. 따라서 모든 입력 변수를 사용한 모델을 KOMPSAT-3/3A 영상으로부터 가장 정확한 탁도를 산출할 수 있는 최적의 모델로 결정하였다. Fig. 7은 테스트 자료에 대해 현장관측 탁도와 최적의 XGBoost 모델로 예측된 탁도 사이의 상관도이다. 관측된 탁도와 예측된 탁도 사이의 상관계수는 0.837로 높았고, 평균절대오차(mean absolute error, MAE)와 평균편향오차(mean bias error, MBE)는 각각 1.77 NTU와 1.39 NTU로 매우 작았다.

Table 4. Performance of the XGBoost turbidity retrieval models

Fig. 7. Comparison of the turbidity values measured by the water quality monitoring system (in situ) and those predicted by the optimal XGBoost model.

다중분광 영상과 기계학습을 통해 개발된 기존의 모델들(El-Alem et al., 2022; Ma et al., 2021)은 실측 탁도와 비교하여 ~0.99의 높은 상관계수를 보이며 RMSE 10–35 NTU의 탁도를 예측하였다. 이 연구에서 개발된 KOMPSAT-3/3A 탁도 산출 모델은 기존 모델들에 비해 실측 탁도와의 상관성이 다소 낮았다. 이는 이 연구에 사용된 KOMPSAT-3/3A 영상에 대기보정이 적용되지 않았고, Landsat 위성 및 Sentinel-2에 비해 적은 분광 정보를 제공하기 때문일 수 있다. 그러나 이 연구에서 개발된 모델은 기존 모델들에 비해 훨씬 작은 RMSE를 가지는 탁도를 예측하여 우수한 성능을 보였다. 이 연구의 테스트 자료에 포함된 실측 탁도 최댓값은 18.8 NTU로 비교적 작다. 따라서 매우 높은 탁도에 대한 예측 성능까지 기존 모델들에 비해 우수하다고 단정하기 어려우며, 고탁도의 수체를 관측한 KOMPSAT-3/3A 영상을 추가적으로 확보하여 분석할 필요가 있다.

최적의 탁도 산출 모델에 대한 변수중요도는 Fig. 8에 도시한 바와 같다. 이 연구에서 제안한 nNDTI가 탁도 산출에 있어 가장 중요한 변수로 사용되었는데, 이는 nNDTI가 다른 변수들에 비해 탁도를 잘 설명할 수 있기 때문으로 판단된다(Fig. 6). 청색 밴드의 TOA 반사율과 AOT가 차순위의 중요 변수로 분석되었다. 청색 밴드의 반사율은 nNDTI 계산에 사용된 변수이다. 그리고 청색 파장 대역은 에어로졸에 매우 큰 영향을 받는데, XGBoost 모델에 AOT를 입력 변수로 고려해 줌으로써 그 영향이 보상되었을 수 있고 그 결과로 변수중요도에서 높은 우선순위를 가질 수 있다.

Fig. 8. Variable importance of the optimal XGBoost model.

4) KOMPSAT-3/3A 기반 고해상도 탁도 매핑

최적의 탁도 산출 모델을 이용하여 KOMPSAT-3/3A 영상에 포함된 모든 수계에 대해 탁도 매핑을 수행하였다. 이를 위해 먼저 NDWI를 기반으로 수계 영역을 제외한 나머지 영역을 마스킹 하였다. 영상 마다 조도, 대기 등의 영향도가 다르기 때문에 수계 영역 추출을 위한 NDWI 임계값은 영상에 따라 다르게 설정해 주었다. Fig. 9는 경기도 구리시의 구리 측정소 주변 한강 본류에 대한 탁도 매핑 영상(2014년 12월 5일)이다. 흰색 화살표는 수계의 흐름 방향을 나타낸다. 유로가 굽은 지점에서 부유물질이 쌓임에 따라 탁도가 높아지는 것이 관측되며, 일부 수계 가장자리에서는 부유물이 축적되어 12 NTU 이상의 고탁도가 관측된다. 영상에서 지형 지물의 그림자가 있는 영역은 0 NTU에 가까운 매우 낮은 탁도가 산출되었는데, 이는 정확한 탁도가 아니라고 할 수 있다.

Fig. 9. The XGBoost turbidity map retrieved from KOMPSAT-3 image obtained on December 5, 2014.

Fig. 10은 2015년 3월 30일, 2015년 10월 29일, 2018년 11월 1일, 2019년 11월 29일에 강원도 원주시의 원주 측정소 주변 하천이 촬영된 KOMPSAT-3/3A 영상에 대한 탁도 매핑 결과이다. Fig. 10(a)의 흰색 화살표는 수계의 흐름 방향을 나타낸다. 수 십 m의 하천 폭을 가지는 지역으로 KOMPSAT-3/3A 영상의 고해상도 관측에 의해 탁도의 공간 분포 파악이 가능하다. 이 지역 하천의 탁도는 가을철(Fig. 10(b)–(d))보다 봄철(Fig. 10(a))에 더 높은 경향을 보임을 알 수 있다. 이는 봄철에 발생하는 황사와 증가하는 수온으로 인한 유기물의 활동성 증대로 인해 탁도가 더 높게 발생하기 때문인 것으로 추측된다. 또한 이 지역 하천의 가을철 탁도는 2015년에서 2019년으로 갈수록 조금씩 높아지는 것을 볼 수 있다. Fig. 10을 통해 KOMPSAT-3/3A 기반 하천 탁도 산출 모델이 소규모 수계의 수질 분석 및 시공간 모니터링에 유용하게 사용될 수 있음을 확인할 수 있다.

Fig. 10. The XGBoost turbidity maps retrieved from KOMPSAT-3/3A images obtained on (a) March 30, 2015, (b) October 29, 2015, (c) November 1, 2018, and (d) November 29, 2019.

북한강과 한강 본류가 만나는 팔당댐 주변 지역에 대한 2018년 5월 23일과 2019년 11월 14일 KOMPSAT-3 영상의 탁도 매핑 결과는 Fig. 11(a)와 (b)에 도시하였다. 원주 측정소 주변 하천과 유사하게 봄철인 5월과 가을철인 11월의 탁도 분포가 확연한 차이를 보이는 것을 볼 수 있다. 5월에는 대규모 하천보다 폭이 좁은 소규모 하천의 탁도가 더 높게 산출되었으며, 이로부터 소규모 하천의 탁도 모니터링이 중요함을 다시 확인할 수 있다. 북쪽의 북한강, 동쪽의 남한강, 남쪽의 경안천이 만나는 지점에서는 여러 수계가 합류하여 부유물질이 혼합되면서 탁도의 공간적인 변화가 발생하는 것도 확인할 수 있다.

Fig. 11. The XGBoost turbidity maps retrieved from KOMPSAT-3/3A images obtained on (a) May 23, 2018 and (b) November 14, 2019. (c) An enlarged image corresponding to the area of red box in (a). (d) The locations of sluice gates in the region of (c).

Fig. 11(a)로부터 일부 수계의 가장자리에서 40 NTU 이상의 고탁도가 산출된 것을 볼 수 있으며, 해당 영역을 확대하여 Fig. 11(c)에 나타냈다. Fig. 11(d)는 하천관리지리정보시스템(http://www.river.go.kr/)에서 제공되는 팔당댐 주변의 배수 통문 지점을 보여준다. 탁도 매핑 결과에서 고탁도가 산출된 위치는 배수 통문의 위치와 매우 유사하며, 이는 하천으로의 배수 시설 주변에 탁수가 집적될 수 있음을 보여준다.

5. 결론

이 연구에서는 KOMPSAT-3/3A 영상으로부터 고해상도의 탁도 공간정보를 산출하기 위해 XGBoost 기법 기반의 탁도 산출 모델을 개발하였다. 탁도에 따른 KOMPSAT-3/3A 영상의 TOA 분광반사율을 분석하여 한강 수계 하천의 탁도를 설명할 수 있는 새로운 분광지수인 nNDTI를 제안하였다. 연구에 사용된 KOMPSAT-3/3A 영상 개수에 제한이 있었으나, Landsat-8 OLI 영상의 TOA 분광반사율에 대한 교차검보정을 통해 기계학습을 위한 훈련 자료를 구축할 수 있었다. 청색, 녹색, 적색 및 근적외선 파장 대역의 TOA 반사율과 분광지수(NDVI, NDWI, NDTI, nNDTI), MODIS 대기 산출물(AOT, WV, OZ)을 입력 변수로 활용한 XGBoost 탁도 산출 모델은 현장관측 탁도와 비교하여 2.70 NTU의 RMSE 및 14.70%의 NRMSE를 가지는 탁도를 예측하여 우수한 성능을 보였다. XGBoost 탁도 산출 모델에는 이 연구에서 새롭게 제안한 nNDTI가 가장 중요한 변수로 사용되었고, 청색 밴드의 TOA 반사율과 AOT가 각각 두 번째와 세 번째 중요 변수로 사용되었다. 개발된 탁도 산출 모델을 KOMPSAT-3/3A 영상에 적용하여 고해상도의 하천 탁도를 매핑할 수 있었고 계절, 하천의 규모 및 지형 지세와 유동 특성, 하천 배수구 위치 등에 따른 탁도의 시공간적 변동 분석이 가능하였다. 이 연구를 통하여 고해상도의 KOMPSAT-3/3A 위성 영상이 소규모 수계에 대한 탁도 공간정보 산출에 효과적으로 활용될 수 있음을 확인하였으며, 우리나라의 수질 관리 모니터링 분야에 크게 기여할 수 있으리라 기대된다.