농작물은 짧은 기간동안 빠른 생장으로 인해 다른 피복에 비해 보다 짧은 위성의 재방문 주기를 요구한다. 이를 위해 2025년 발사 예정인 차세대 위성 4호기인 일명 농림위성은 타 극궤도 위성과 다르게 전 지구를 관측하지 않고 한반도 일대를 3일 주기로 촬영할 수 있게 설계되었다. 이러한 농림위성도 타 지구관측위성처럼 다양한 레벨 및 활용 산출물을 생산하기 위한 발사에 앞서 알고리즘 및 시스템 체계를 구축하고 있다. 지구 관측 광학위성 기반 다양한 육상 산출물 중 시계열 정규식생지수영상(Normalized Difference Vegetation Index, NDVI)은 농작물의 생육 모니터링을 위해 활용되는 대표적인 산출물이다(Houborg and McCabe, 2016; Kasampalis et al., 2018). 시계열 정규식생지수영상의 활용에 있어 광학위성에서 자주 나타나는 구름은 주기적 모니터링에 있어 가장 큰 방해물로 지적되고 있다(Roy et al., 2008; Moreno-Martínez et al., 2020).

이러한 제약점을 해결하기 위해 일정 주기동안 구름의 영향을 받지 않는 화소들을 합성하는 구름 제거 주기합성기법들이 제시되었다. 실제 주기해상도가 높은MODIS와 AVHRR 위성영상을 대상으로 목적에 맞게 다양한 주기로 합성된 구름 제거 식생지수영상이 서비스되었고, 광역의 농경지와 산림 식생의 시계열 모니터링에 활용되었다(van Leeuwen et al., 2006; Albarakat and Lakshmi, 2019; Kumhálová and Matêjková, 2017). 그러나 국내 주기합성된 정규식생지수의 시계열 패턴을 보면 장마철의 경우 장기간 구름으로 주기합성 영상에서도 구름의 영향을 받은 것으로 예상되는 낮은 NDVI 값을 보인다(Kim and Eun, 2022). 특히 장마철은 농작물의 생장이 가장 빠른 시기로 그 모니터링일 중요한데 장기간 구름은 식생 정보의 수집 및 분석에 큰 제약을 준다.

이러한 문제를 해결하고 농작물의 연속적인 시계열 모니터링을 위해 해당 기간의 결측자료를 결측 전과 후 자료를 이용해 추정하는 결측 보정 기법(gap-filling technique)이 적용되었다. MODIS 센서의 경우 결측 보정한 자료를 기존 산출물 이름에 Gap-Filled, Smoothed (GFS)를 붙여 명명하였으며 대표적으로MOD09GFS와MOD15GFS가 서비스 되고 있다(Gao et al., 2008). 단, 이러한 결측 보정은 결측 자료가 전체 자료의 25%를 차지하는 경우 결측 보정 후 결과가 왜곡될 경우를 예상하여 시기적 결측 보정을 수행하지 않는다(Gao et al., 2008). 이 경우 결측 화소 주변에 동일 피복이며, 좋은 품질의 화소를 찾아 참조화소로 정의하고 결측 보정 시 참조화소의 시계열 정보를 사용한다(Gao et al., 2008). Sentinel-2A/B위성의 경우 결측 보정된 산출물을 제공하지 않고 Decomposition and Analysis of Time Series Software (DATimeS)라는 결측 보정 프로그램을 제공하며, 해당 프로그램에서는 시계열 위성영상을 대상으로 다양한 시기적 결측 보정 기법을 제공한다(Belda et al., 2020).

시기적결측보정기법은크게내삽(interpolation)과평활(smoothing) 기법으로 나뉘며, 내삽 기법은 선형 내삽(Linear Interpolation, LI)과 최대값 합성(Maximum Value Composite, MVC) 기법으로 나눌 수 있다(Yu et al., 2021; Garioud et al., 2021). 내삽 기법은 섭동을 줄이는 가장 간단한 방법이나, 식생 성장이나 시계열의 많은 특징들을 고려하지 않는다는 한계점이 있다(Garioud et al., 2021). 평활 기법은 단일 계절 자료에 적합한 비대칭 함수(asymmetric function)와 다중 계절 자료에 적용이 가능한 필터링 기법이 있다(Kandasamy et al., 2013). 필터링 기법은 지역(local)과 전체(overall) 필터링으로 나눌 수 있으며, 해당 기법들은 구름이나 그림자로 인한 노이즈의 영향을 받아 낮은 NDVI를 보인다는 한계점이 있다(Yu et al., 2021). 최근 구름이나 그림자의 영향을 최소화하는 새로운 기법들이 제시되었으나, 대기 변이와 Bidirectional Reflectance Distribution Function (BRDF) 영향을 줄이지 못한다는 문제점도 제시되었다. 최근 국내에서는 이상 기후로 인해 장마뿐만 아니라 국지성 호우가 잦아지며 결측 보정 기법의 개 발이 더욱 필요해졌다. 본 연구에서는 차세대중형위성 4호(CAS500-4) 혹은 농림위성과 가장 비슷한 Sentinel-2A/B 위성의 주기 합성 NDVI 영상을 대상으로 결측 보정 기법을 제시하려 한다. 또한 기존 시계열 결측 보정 기법에 공간적 결측 보정기법이 결합된 시공간적 결측 보정 기법을 제시하고자 한다.

2.1. 연구 자료

본 연구에서는 단일 작물로 넓게 분포한 당진시 석문면에 위치한 교로리 일대를 중심으로 50 km × 50 km (2,500 km²) 면적을 포함하는 논 지역을 선택하였다(Fig. 1). 석문면 교로리 일대는 전부 논으로 구성되어 있으며, 5월 모내기철을 시작으로 10월 추수를 진행한다. 이 지역에 대해 2023년 모내기철인 5월부터 추수가 된 10월까지, 농림위성의 밴드 범위 및 사양이 동일할 것으로 예상되는 Sentinel-2A/B 위성영상의 표면반사율 자료인 L2A를 수집하였다(CopernicusOpen AccessHub, https://scihub.copernicus.eu/dhus/). Sentinel-2A/B위성영상은 10 m급 공간해상도와 두 쌍둥이 위성이 동시 운영되면서 5일 주기로 자료가 제공되고 있으며, 290 km의 촬영폭으로 촬영된다. 위성은 공식적으로 정규식생지수를 제공하지 않고 있어, 대기보정된 반사율인 L2A 자료 중 665 nm중심파장의 Band 4와 833 nm의 Band 8영상을 이용해 일별 정규식생지수영상을 제작하였으며, 해당 자료는 Table 1과 같이 촬영일자가 표시된 파일로 명명되었다. 이 일별 식생지수영상은 구름의 영향을 줄이기 위해 세개의 일별 영상을 최대값 합성기법(Maximum NDVI Composite, MNC)에 적용하여 구름의 영향을 줄인 주기합성영상을 획득하였다(Kim and Eun, 2022). 주기합성 결과는 Fig. 2와 같으며, 주기합성 후에도 구름의 영향을 받은 화소는 주기 합성 후 구름 마스크(mask) 자료를 이용하여 정의한다(Fig. 2e). 이와 같이 제작된 세개의 합성영상과 구름 마스크 정보를 이용하여 구름이 있는 결측화소에 대해 결측보정기법을 적용하였다.

Figure 1. Google color image of the study area (Dangjin rice paddy).

Figure 2. Sentinel-2A/B daily NDVI images: (a) July 21st, (b) July 26th, (c) August 5th, (d) composited NDVI image, and (e) cloud mask extracted from the composite image (white: cloud).

Table 1 . Types and specifications of the dataset used.

Date of daily NDVI data	List of composite data
May 2, 2023 May 22, 2023 June 6, 2023	C4_20230502_20230606_R003_122_L3VP_NDVI.tif C4_20230502_20230606_R003_122_L3VP_Cloud.tif
June 16, 2023 July 1, 2023 July 6, 2023	C4_20230616_20230706_R003_122_L3VP_NDVI.tif C4_20230616_20230706_R003_122_L3VP_Cloud.tif
July 21, 2023 July 26, 2023 August 5, 2023	C4_20230721_20230805_R003_122_L3VP_NDVI.tif C4_20230721_20230805_R003_122_L3VP_Cloud.tif
August 20, 2023 August 25, 2023 September 9, 2023	C4_20230820_20230909_R003_122_L3VP_NDVI.tif C4_20230820_20230909_R003_122_L3VP_Cloud.tif
September 24, 2023 October 4, 2023 October 29, 2023	C4_20230924_20231029_R003_122_L3VP_NDVI.tif C4_20230924_20231029_R003_122_L3VP_Cloud.tif

2.2. 연구 방법

2.2.1. 시공간적 결측 보정 기법

일반적인 시간 결측 보정에 앞서 주기 합성 이후 남아 있는 구름으로 인한 결측 화소의 수를 줄이기 위한 방안으로 공간적 결측 보정 기법을 개발하였다. 공간적 결측 보정 기법은 기본적으로 동일 피복에, 결측 화소 주변에 품질이 좋은 화소 정보를 사용한다. 국내와 같이 상대적으로 좁은 규모의 농경지를 고려할 때, 공간적 결측 보정 기법이 넓은 결측 지역에 적용하는 것은 왜곡을 가져올 수 있다. 이에 따라 본 연구에서는 국내 농경지 규모를 고려하여 5개 화소 이하의 소규모 결측 지역을 대상으로 공간적 결측 보정이 수행되었다. 결측화소가 6개 이상인 상대적으로 큰 면적의 결측 지역은 시간적 결측 보정 기법을 정의한다. 따라서 결측 화소를 중심으로 5 × 5 윈도우 크기안에 가장 가까운 참조 화소들의 평균으로 정의하고 정규식생지수 값을 취한다. 이를 통해 소규모 결측 지역과 결측 지역의 경계지역이 보정되게 된다. 공간적 결측 보정 기법이 적용된 영상을 대상으로 시간적 결측 보정이 수행된다. 시간적 결측 보정을 위해서는 결측 화소가 획득된 날짜를 기준으로 전에 촬영된 영상 2개, 후에 촬영된 영상 1개가 요구된다. 사전 연구에서는 본 연구지역과 동일한 당진 논지역의 2021년 Sentinel-2A/B NDVI영상을 대상으로 가우스 과정 회귀(Gaussian Process Regression, GPR) 모델, 선형 내삽, 최근린 기법을 비교, 적용하였다. 시험 결과 세개의 기법 중 GPR 기법이 가장 낮은 오차(0.173 RootMean Sqaured Error [RMSE], 0.131 MeanAbsolute Error [MAE])를 나타내었으며, 기존 연구 사례에서도 많은 알고리즘 중 GPR 알고리즘이 가장 낮은 RMSE를 보이는 것으로 보고되었다. 이에 따라 본 연구에서도GPR 기법을 선택, 적용하였다(Belda et al., 2020). GPR 모델은 회귀 및 함수 근사를 위한 최첨단 통계 방법으로 광학 위성영상에서 생물리학적 매개변수를 추정하는 데 성공적으로 적용한 모델이다(Belda et al., 2020). 또한 비모수적 모델링을 통해 평균 추정치와 함께 불확성 구간도 제공한다(Belda et al., 2020). 시공간적 결측 보간은 Fig. 3과 같은 절차로 수행되며, 해당 알고리즘은 파이썬(Python)을 이용해 개발된 프로그램(CAS-GPR.exe)으로 구현되었다.

Figure 3. Flowchart of spatio-temporal gap-filling algorithm for Sentinel-2A/B NDVI images.

2.2.2. 검증 방법

공간적 결측 보정의 결과를 분석하기 위해 우선 보정 전과 후의 결측 화소의 수를 비교한다. 시간적 결측 보정 후에는 모든 결측 화소는 보정되기 때문에 보정 후 NDVI 값이 왜곡없이 추정되었는지 확인이 필요하다. 이를 위해서는 구름이 있는 시간과 동일 지점에서 현장 측정이 수행되어야 하나, 구름이 많은 시기에 현장 측정의 정확도도 떨어지기 때문에 현실적으로 어렵다. 또한 NDVI의 현장 측정을 위해서는 분광반사율을 측정하고 NDVI 값을 추출해야 하나, 센서별 특징이 달라 절대적 산출은 어렵다. 일반적으로 시계열 NDVI 패턴에서 결측으로 인해 감소한 값이 회복하는 지를 육안으로 판단하기도 하나, 이는 정량적인 검증이 어렵다는 한계가 있다. 따라서 본 연구에서는 정상 화소 10,000개를 랜덤으로 추출하고, 이를 고의적으로 결측 화소로 수정 후 해당 결측 보정기법을 적용하여, 원 정상 화소의 NDVI 값과 비교하는 방법으로 검증을 수행하였다. 또한 정량적 검증을 위해 상관계수(correlation coefficient, r), 평균 제곱근 오차(RMSE), 편향(bias) 등 세 개의 정량적 척도가 사용되었다.

우선 본 연구의 입력자료인 주기 합성 자료는 Fig. 4(a)와 같이 구름으로 인해 전체 영상의 8.1%의 결측이 발생하였다. 이 주기합성 영상은 본 연구에서 개발한 공간적 결측 보정 프로그램(CAS-GPR.exe)을 통해 결측 비율이 3%가 감소된 5.1%의 결측 비율 결과를 산출할 수 있었으며(Fig. 4b), 결측 화소 경계나 좁은 규모의 결측 지역이 보정된 것을 확인할 수 있다. 해당 영상의 공간적 결측 보정은 대상 2분의 처리 시간이 소요되었으며, RTX 4070 TIGPU, i9-13900K CPU, 128 GB 메모리를 보유한 PC에서 테스트 되었다. 이렇게 공간적 결측 보정된 영상은 여전히 남아 있는 지역에 대해GPR 모델 기반 시간적 결측 보정을 수행하였으며, Fig. 4(c)와 같이 모든 화소에 대해 결측 보정이 된 결과를 육안으로 확인할 수 있다. 이러한 시간적 결측 모듈은 동일 사양의 PC에서 3분 정도 소요되었다.

Figure 4. Cloud-free composited NDVI image. (a) C4_20230924_20231029_R003_122_L3VP_NDVI.tif, (blue: cloud mask). (b) Spatial gap-filling NDVI image (blue: cloud mask) (c) Temporal gap-filled NDVI image.

구름양에 따른 시간적 결측 보정 결과의 정량적 분석을 위해 두 시기 자료를 비교하였다. 우선 상대적으로 구름양이 많은 7%의 구름이 남아있는 전, 후 영상을 사용해 보정한C4_20230721_20230805_R003_122_L3VP영상(Fig. 5a)과 평균 4%의 구름양이 남아있는 전, 후 영상을 사용해 보정한C4_20230820_20230909_R003_122_L3VP영상(Fig. 5b)을 비교하였다. 정상 화소를 결측 화소로 처리 후 시간적 결측 보정한 결과와 원 정상 화소의 NDVI값을 비교한 결과 Fig. 5와 같이 10,000개의 화소를 대상으로 15% 이하의 RMSE를 나타내는 것으로 나타냈다. 특히 구름양이 상대적으로 많은 입력 자료의 경우 NDVI가 과추정 되는 경향(Fig. 5a)을 볼 수 있었으나 전반적으로 일대일 대응선에 밀집하는 경향을 볼 수 있었다.

Figure 5. Scatter plot of NDVI before and after gap-filling. (a) C4_20230721_20230805_R003_122_L3VP. (b) C4_20230820_20230909_R003_122_L3VP (dotted line: 1 to 1 line, solid line: trend line).

2025년에 발사될 농림위성의 처리 알고리즘 중 하나로써 구름의 영향을 최소화한 시계열NDVI를 산출하기 위한 결측 보정 알고리즘이 개발되었다. 이를 위해 농림위성과 가장 비슷한 Sentinel-2A/B 위성영상이 사용되었으며, 자체적으로NDVI 영상을 제작하고 계획된 농림위성영상 크기로 만들어 사용하였다. 1차로 일별 NDVI 영상 3장을 이용해 주기 합성 영상을 제작하고, 남아있는 결측 화소를 공간적으로 보정한 후, 시계열 자료를 이용하여 시간적으로 보정하였다. 이러한 공정은 하나의 프로그램으로 개발되어 5분 내에 자동으로 산출되었으며, 이 프로그램은 농림위성영상을 분석하는 전문가용으로 제공될 예정이다. 다양한 시간적 결측 보정기법이 개발되었으며, 연구자용으로 오픈되어 제공되는 프로그램도 있으나 대부분 시간적 보정기법을 지원한다. 본 연구에서 개발된 알고리즘은 공간적 보정기법도 적용되어 효율성과 처리 속도 측면에서 기존 알고리즘보다 향상되었다. 단, 해당 시공간적 결측 보정 알고리즘의 경우 결측 시기 후에 획득된 영상이 있어야 실행시킬 수 있어 위성영상 획득과 함께 실시간으로 실행될 수 없다는 한계점이 있다. 이와 함께 구름 양에 따라 보정 정확도의 변이가 발생하여 결측 화소의 수에 따른 정확도를 정량적으로 분석할 필요가 있다. 본 연구에서는 2023년 농경지지역을 대상으로 적용, 검증이 수행되었으며, 향후 연구에서는 다양한 피 복의 대상지역과 보다 많은 시계열 자료에 적용할 예정이다. 또한 주기 합성영상뿐만 아니라 일별 식생지수영상이나 타 생물리학적 인 자의 결측 보정에 적용할 계획이다.

본 논문은 농촌진흥청 공동연구사업(RS-2021-RD009991)의 지원으로 이루어졌으며, 이에 감사드립니다.

No potential conflict of interest relevant to this article was reported.