A Real-time Correction of the Underestimation Noise for GK2A Daily NDVI

GK2A 일단위 NDVI의 과소추정 노이즈 실시간 보정

Abstract

Normalized Difference Vegetation Index (NDVI) is utilized as an indicator to represent the vegetation condition on the land surface in various applications such as land cover, crop yield, agricultural drought, soil moisture, and forest disaster. However, satellite optical sensors for visible and infrared rays cannot see through the clouds, so the NDVI of the cloud pixel is not a valid value for the land surface. This study proposed a real-time correction of the underestimation noise for GEO-KOMPSAT-2A (GK2A) daily NDVI and made sure its feasibility through the quantitative comparisons with Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI and the qualitative interpretation of time-series changes. The underestimation noise was effectively corrected by the procedures such as the time-series correction considering vegetation phenology, the outlier removal using long-term climatology, and the gap filling using rigorous statistical methods. The correlation with MODIS NDVI was higher, and the difference was lower, showing a 32.7% improvement compared to the original NDVI product. The proposed method has an extensibility for use in other satellite products with some modification.

Normalized Difference Vegetation Index (NDVI)는 간단한 수식으로 식생 활력도를 파악할 수 있기 때문에, 토지피복, 곡물수확, 농업가뭄, 토양수분, 산림재해 등의 분야에서 지표면의 식생 상태를 나타내는 지시자로 활용되고 있다. 그러나 가시광선과 적외선을 감지하는 광학위성센서는 구름을 투과하여 센싱(sensing)할 수 없으므로, 구름 화소에 대해서는 지표면의 NDVI가 계산될 수 없다. 본 연구에서는 천리안위성 2A호(GEO-KOMPSAT-2A, GK2A) 일단위 NDVI의 과소추정 노이즈를 개선하기 위한 실시간 보정기법을 개발하였으며, Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI와의 정량적 비교 및 시계열 변화에 대한 정성적 해석을 통해 제안된 기법의 유용성을 확인하였다. 식생의 생장 주기를 반영한 시계열 보정, 장기간 기후치를 이용한 이상치 제거, 엄밀한 통계기반의 결측화소 복원 등의 과정을 거친 보정 NDVI에서는 구름에 의한 과소추정이 효과적으로 개선되었다. MODIS NDVI와의 상관성이 원값보다 높고, 그 차이가 감소한 경향을 보여주었으며, 특히 여름철에는 32.7%의 향상률을 보였다. 본 연구에서 제안한 방법은 일부 수정을 거쳐 타 위성산출물에도 적용될 수 있는 확장성을 가지고 있다.

1. 서론

일반적으로 건강한 식물의 엽록소는 광합성을 위해서 적색광(Red)을 많이 흡수하고 이에 따라 식물 내부 온도 상승을 막기 위해서 근적외선(Near-Infrared, NIR)을 투과 및 반사한다(Jensen, 2016). 정규화 식생지수(Normalized Difference Vegetation Index, NDVI)는 식생의 건강 상태에 따라 적색광과 근적외선을 흡수 또는 반사하는 특성을 이용한 정규화 식을 통해서 계산되며(식(1))(Rouse et al., 1974), 식생 활력도에 따라 0–1 사이의 값으로 나타난다. NDVI는 간단한 수식으로 식생 활력도를 파악할 수 있기 때문에, 토지피복, 곡물수확, 농업 가뭄, 토양수분, 산림재해 등의 분야에서 지표면의 식생 상태를 나타내는 지시자로 활용되고 있다(Peters et al., 2002; Domenikiotis et al, 2002; Han et al., 2010; Huang, et al., 2013; Usman et al., 2015).

\(\begin{aligned}N D V I=\frac{N I R-\operatorname{Red}}{N I R+\operatorname{Red}}\\\end{aligned}\)       (1)

가시광선과 적외선을 감지하는 광학위성센서는 구름을 투과하여 센싱할 수 없으므로, 구름 화소에 대해서는 지표면의 NDVI가 계산될 수 없다(Bradley et al., 2007; Li et al., 2021). 이렇게 미관측된 NDVI를 추정하기 위해서 일정 기간 동안의 이동평균(Lee and Han, 2014) 또는 가장 높은 값을 취하거나(Maximum ValueComposite, MVC) (Holben, 1986; Van Leeuwen et al., 1999), 최소자승법으로 구해진 다항식을 통해 평활화시키는 Savitzky-Golay (SG) 필터(Chen et al., 2004), 3차 다항식을 이용한 시계열내삽(cubic spline) (Vorobiova and Chernov, 2017) 등의 시계열 통계법이 활용되고 있다. 대표적인 NDVI 산출물인 미항공우주국(NationalAeronautics and SpaceAdministration, NASA) Moderate Resolution Imaging Spectroradiometer(MODIS)의 경우, 구름 영향을 최소화한 MVC 합성산출물이 16일 단위 및 월 단위로 제공되고(Kamel et al., 2015), 북미 지역에 한해서는 결측화소를 복원한(gap-filled) 자료(MOD09Q1: https://lpdaac.usgs.gov/products/mod09q1v006/)가 제공되고 있다.

우리나라에서는 2018년 12월에 발사되어 동아시아 영역의 기상관측 임무를 수행하는 정지궤도 위성인 천리안위성 2A호(GEO-KOMPSAT-2A, GK2A)는 Advanced Meteorological Imager (AMI)의 10분 단위 관측을 기반으로 일 단위의 NDVI 산출물을 제공하고 있다(Seong and Han, 2019). 우리나라의 여름은 집중호우나 장마로 인해 여름철에 상당히 오랜 기간 구름이 존재하기 때문에, 농지나 산림임에도 불구하고 0.3 미만의 낮은 NDVI 값을 가지는 화소들이 다수 발견되었다(Fig. 1). 이렇듯 구름 등으로 인한 NDVI 과소추정 노이즈의 개선이 필요하며, GK2A NDVI 자료가 일단위로 생산되기 때문에 이러한 보정처리는 실시간으로 이루어질 필요가 있다.

Fig. 1. GK2A daily NDVI maps between July 19 and August 19, 2020. Many low-value pixels were found in the cropland and forests.

본 연구에서는 GK2A 일단위 NDVI의 과소추정 노이즈를 개선하기 위한 실시간 보정기법을 제안하고자 한다. 이동평균과 MVC에 기초한 시계열 보정, 통계적 이상치의 제거 및 결측화소 복원 등이 과정을 거쳐 통해 GK2A 일단위 NDVI를 개선하였으며, 제안된 기법에 대한 성능 평가를 위해서 MODIS NDVI와의 정량적 비교 및 시계열 변화에 대한 정성적 해석을 수행하였다.

2. 연구 자료 및 방법

1) GK2A NDVI 산출물

GK2A 위성은 16개 채널 가진 고성능 기상센서 AMI를 통해서 0.5–2 km 해상도로 2–10분 단위의 관측을 수행하고, 52종의 다양한 기상산출물을 제공하고 있다. GK2A NDVI 산출물은 2 km 공간 해상도, 1일 간격으로 생산되며, 2020년부터 1월 1일부터 2021년 12월 31일까지 731장의 영상을 사용하여 실험하였다. 연구 지역은 한반도 일대로, 위도 33.0°–43.2°N, 경도 124.0°–130.7°E의 범위이다.

2) 노이즈 보정

GK2A NDVI와 함께 제공되는 Data Quality Flag(DQF)를 사용하여 호수, 저수지 등의 육수(inland water) 부분을 제외하고 지표면 화소만을 대상으로 하였다(Fig. 2(a)). 1년치를 모아서 연말에 시공간 보정을 수행하는 산출물(year-end product)은 타깃 시점에 대한 전·후 시기 영상이 모두 존재하므로 시계열내삽이나 SG 필터 등을 사용할 수 있으나, GK2A NDVI 산출물은 실시간으로 타 산출물의 입력자료로 사용되므로 실시간 보정처리를 해야 하기 때문에, 이전 시기 영상과 NDVI 기후값만을 참조하여 보정처리를 수행해야 한다. 본 연구에서는 시계열 보정, 이상치 제거, Discrete Cosine Transform-based Penalized Least Square Regression (DCTPLS) 결측복원 등의 일련의 과정을 거쳐 GK2A NDVI 개선을 수행하였다(Fig. 2). 먼저, 시계열 보정은 이동평균과 MVC에 기초하며, 식생의 생장 주기를 고려하였다. 식생 생장이 감소하는 NDVI 하강기인 9월부터 3월까지는 이동평균을 적용하였다. NDVI 하강기에 MVC를 적용하면 이전 시기의 높은 값이 타깃 시점에 그대로 반영될 수 있기 때문에, 타깃 시점과 과거 9일치를 합쳐서 총 10일치의 NDVI 값을 이동평균 방식으로 합성하였다(Fig. 2(b)). 식생 생장이 증가하는 NDVI 상승기인 4월부터 8월까지는 이동평균과 MVC가 결합된 방식을 적용하였다. 우리나라에서는 NDVI가 4월부터 상승하기 시작하여 7월이나 8월에 최고치를 나타내는 것이 일반적이지만, 여름철에는 지표면이 아니라 구름 화소의 NDVI가 계산됨으로써 NDVI 과소추정의 경향이 있다. 이때 이동평균을 그대로 취할 경우 과소추정치가 함께 포함되어 실제보다 낮은 NDVI 값을 가질 수 있기 때문에, 이를 방지하기 위하여 NDVI 상승기에는 타깃 일자의 NDVI 원값과 과거 9일치 시계열 보정값의 평균을 구한 뒤, 그 평균값과 원값 중 최대값을 취하였다(Fig. 2(c)).

Fig. 2. Procedure of the real-time correction of the underestimation noise for GK2A daily NDVI.

이러한 시계열 보정을 수행한 후에도 여전히 과소추 정치가 남아있는데, 이는 기후값과의 비교를 통해 보정하였다. 20년동안의(2002–2021년) MODIS월단위NDVI 영상으로부터 각 월별 최소값을 추출하였고, 각 화소별 최소값보다 훨씬 낮은 값이 존재한다면 이를 이상치로 판단하여 제거하였다(Fig. 2(d)). 특히 2020년 7월과 8월은 유례없는 집중호우로 인해(KMA, 2020a; 2020b) GK2A 영상에 구름 화소가 매우 많았고, 이로 인해 NDVI 과소추정 경향이 뚜렷했다(Fig. 3; Fig. 4). MODIS는 1999년과 2002년에 발사된 Terra와 Aqua 위성에 탑재된 다중분광센서로 현재까지 NDVI 산출물을 제공하기 때문에, 본 연구에서는 동일시기 Terra와 Aqua의 NDVI를 평균하여 사용하였다. 월단위 NDVI로부터 각 월별, 각 화소별 20년간 최소값을 도출하였으며, GK2A NDVI 각 화소에 대하여 역대 최소값보다 0.1 이상 더 작을 경우에는 이상치로 간주하여 제거하였다(Fig. 2(d)). 이때, 여러 차례의 민감도 실험을 통해 가장 적절한 임계치로 0.1을 설정하였다(Fig. 5).

Fig. 3. The difference between GK2A and MODIS NDVI in 2020. The values were calculated as GK2A minus MODIS. Many negative values were found in summer.

Fig. 4. The difference between GK2A and MODIS NDVI in 2021. The values were calculated as GK2A minus MODIS. Many negative values were found in summer.

Fig. 5. Ratio of the outlier pixels in land pixels.

제거된 이상치 화소에 대해서는 DCT-PLS 기반의 결측복원을 수행하였다. DCT-PLS 기법은 영상의 모든 화소를 이산코사인변환(DCT)하고, 변환된 데이터의 분포로부터 벌점최소제곱법회귀(PLS)를 수행하여 결측화소를 예측하는 기법으로(Wang et al., 2012), 회귀계수 추정에 사용되는 손실함수에 과적합(overfitting) 방지를 위한 페널티항이 추가된 것이다(Wang et al., 2014). 이 기법은 국지예보모델(Local DataAssimilation and Prediction System, LDAPS)의 상대습도와 풍속, MODIS Land Surface Temperature (LST)와 NDVI 결측복원 실험에서 상관계수 0.95 이상의 원본과의 일치도를 나타내, 격자 자료의 결측복원에 효과적인 것으로 평가된다(Youn et al., 2020). 그러나 전체 영상에서 이상치 화소의 비율의 너무 많을 경우 DCT-PLS의 적용이 어렵기 때문에, 본 연구에서는 육지 화소 대 이상치 화소 비율이 20% 미만인 영상을 DCT-PLS 결측 복원에 사용하였다(Fig. 2(e)). 이상치 화소 비율이 20% 이상인 영상은 731개 중 23개 영상이었으며, 이들 영상은 삭제하고 직전일 보정영상으로 대체하였다(Fig. 2(f)). 이러한 과정을 거쳐 최종적으로 개선된 GK2A NDVI에 대한 평가는 MODIS NDVI와의 비교를 통해 이루어졌으며, 상관계수(Correlation Coefficient, CC)와 평균편의차(Mean Bias Difference, MBD), 평균절대차(Mean Absolute Difference, MAD)를 통해 분석하였다. 이때 시간 단위 통일을 위해 GK2A 일 단위 자료를 월 평균하여 MODIS 월 단위 자료와 비교하였다.

3. 결과 및 토론

일사량이 비교적 충분했던 2021년에 비해, 기록적으로 장마가 길었던 2020년은 여름철 구름으로 인한 GK2A NDVI의 과소추정 경향이 심하게 나타났으나(Fig. 3, Fig. 4), 노이즈가 보정된 GK2A NDVI에서는 문제시되었던 여름철 NDVI의 과소추정 경향이 크게 감소하였다(Fig. 6, Fig. 7). 보정된 GK2A NDVI는 원본 GK2A NDVI보다 MODIS NDVI와의 상관성이 높고(CC=0.732), 차이가 감소하였으며(MBD=–0.030, MAD=0.078), 여름을 비롯한 모든 계절에서 MODIS NDVI와의 일치성이 향상되었다. 특히 여름철에는 32.7%의 향상률을 보였다((0.110–0.074)/0.110)(Table 1). 보정된 GK2A NDVI는 2020년 8월 강원도 산림지역에서 MODIS NDVI에 비해 큰 값을 가지는 경향을 보였는데(Fig. 6), 2020년 8월 집중호우로 인해 MODIS NDVI 역시 산림지역에서 0.6 이하의 값을 보였기 때문이다(Fig. 8(a), Fig. 8(b)). 그러나, 우리나라 산림의 여름철 NDVI가 0.6 이하인 경우는 드물기 때문에, 이는 MODIS NDVI 역시 구름에 의한 노이즈 보정이 완벽하지 못했기 때문으로 사료되며, 오히려 보정된 GK2A NDVI가 더 합리적인 값을 가지는 것으로 해석된다(Fig. 8). 앞서 Fig. 1에서 2020년 7월 19일부터 8월 19일까지 원본 GK2A NDVI에서 나타난 과소추정 경향이 효과적으로 개선되었음을 Fig. 9를 통해 확인할 수 있었다.

Fig. 6. The difference between the corrected GK2A and MODIS NDVI in 2020. The values were calculated as the corrected GK2A minus MODIS and seemed almost around zero in summer.

Fig. 7. The difference between the corrected GK2A and MODIS NDVI in 2021. The values were calculated as the corrected GK2A minus MODIS and seemed almost around zero in summer.

Table 1. Correlation and differences between (original and corrected) GK2A and MODIS NDVI

Fig. 8. MODIS (a and b) and corrected GK2A (c and d) NDVI in August 2020.

Fig. 9. The distribution of corrected GK2A NDVI from July 19 to August 19, 2020.

GK2A NDVI 보정 효과를 시계열적으로 확인하기 위해서, 상당히 균질한 배경화소가 확보되는 농지와 산림 지역을 선정하여, 농지에서는 연백평야와 나주평야 픽셀, 그리고 산림에서 지리산과 대령산 픽셀에 대한 NDVI 시계열 분포를 확인하였다(Fig. 10). GK2A NDVI원값에서 구름에 의해 여름철 NDVI가 급락하는 노이즈 화소가 효과적으로 개선되었음을 확인할 수 있다. 이는 GK2A NDVI 원값이 가지는 패턴과 함께 식생 생장 주기를 반영하는 점진적인 시계열 변화를 적절히 표현하는 것으로 보인다. 연백평야와 지리산 픽셀의 경우, 2020년 8월에 일주일 이상 구름으로 인해 원본값이 매우 낮게 산출되었는데, 이렇듯 장기간 구름이 계속되면 일련의 보정 처리에도 불구하고 NDVI 값이 정상적인 수준으로 완전히 회복되지 못할 수도 있다(Fig. 10(a), Fig. 10(c)). 이는 당일과 이전시기 자료만을 사용해야 하는 실시간 처리의 불가피한 한계이며, 향후 이에 대한 추가적인 알고리즘 개선이 필요할 것으로 보인다.

Fig. 10. Corrected NDVI time series of (a) Yeonbaek Plain, (b) Naju Plain, (c) Jirisan Mountain, and (d) Taeryongsan Mountain.

4. 결론

본 연구에서는 GK2A 일단위 NDVI의 과소추정 노이즈를 개선하기 위한 실시간 보정기법을 개발하였으며, MODIS NDVI와의 정량적 비교 및 시계열 변화에 대한 정성적 해석을 통해 제안된 기법의 유용성을 확인하였다. 식생의생장주기를반영한시계열보정,장기간기 후치를 이용한 이상치 제거, 엄밀한 통계기반의 결측화소 복원 등의 과정을 거친 보정 NDVI에서는 구름에 의한 과소추정이 효과적으로 개선되었다. MODIS NDVI와의 상관성이 원값보다 높고, 그 차이가 감소한 경향을 보여주었으며, 특히 여름철에는 32.7%의 향상률을 보였다. GK2A 일단위 NDVI는 타산출물의 입력자료로도 활용되기 때문에 실시간 산출이 필요하며, NDVI 보정도 실시간으로 이루어져야 하는데, 본 연구에서는 실시간처리의 한계에도 불구하고 적절한 보정 결과를 도출하였다. 그러나 연말 산출물처럼 전·후 시기 자료를 사용하지 못하기 때문에 구름이 장기간 지속될 때에는 일부 과소추정이 남아있는 경우도 존재했다. 이는 향후 다양한 보조자료의 활용을 통한 알고리듬 개선을 통해 해결될 수 있을 것으로 사료되며, 본 연구에서 제안한 방법은 일부 수정을 거쳐 타 위성산출물에도 적용될 수 있는 확장성을 가지고 있다.