Detecting Phenology Using MODIS Vegetation Indices and Forest Type Map in South Korea

MODIS 식생지수와 임상도를 활용한 산림 식물계절 분석

Abstract

Despite the continuous development of phenology detection studies using satellite imagery, verification through comparison with the field observed data is insufficient. Especially, in the case of Korean forests patching in various forms, it is difficult to estimate the start of season (SOS) by using only satellite images due to resolution difference. To improve the accuracy of vegetation phenology estimation, this study reconstructed the large scaled forest type map (1:5,000) with MODIS pixel resolution and produced time series vegetation phenology curves from Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) derived from MODIS images. Based on the field observed data, extraction methods for the vegetation indices and SOS for Korean forests were compared and evaluated. We also analyzed the correlation between the composition ratio of forest types in each pixel and phenology extraction from the vegetation indices. When we compared NDVI and EVI with the field observed SOS data from the Korea National Arboretum, EVI was more accurate for Korean forests, and the first derivative was most suitable for extracting SOS in the phenology curve from the vegetation index. When the eight pixels neighboring the pixels of 7 broadleaved trees with field SOS data (center pixel) were compared to field SOS, the forest types of the best pixels with the highest correlation with the field data were deciduous forest by 67.9%, coniferous forest by 14.3%, and mixed forest by 7.7%, and the mean coefficient of determination ($R^2$) was 0.64. The average national SOS extracted from MODIS EVI were DOY 112.9 in 2014 at the earliest and DOY 129.1 in 2010 at the latest, which is about 0.16 days faster since 2003. In future research, it is necessary to expand the analysis of deciduous and mixed forests' SOS into the extraction of coniferous forest's SOS in order to understand the various climate and geomorphic factors. As such, comprehensive study should be carried out considering the diversity of forest ecosystems in Korea.

인공위성 영상을 활용한 식생 계절 추정 연구의 지속적 발전에도 불구하고 현장 자료와의 비교를 통한 검증은 다소 미흡한 실정이다. 특히, 다양한 형태로 패치화 되어 있는 우리나라 산림의 경우, 인공위성 영상과 해상도가 일치하지 않는 한계가 존재하여 위성 영상만으로 개엽일을 추정하는데 어려움이 따른다. 본 연구는 식생 계절 추정의 신뢰도를 향상시키기 위해 대축적 임상도(1:5,000)를 MODIS 픽셀 해상도에 맞게 재구성하여 MODIS의 NDVI와 EVI 두 식생지수들로부터 시계열 식물계절 곡선을 생산하여 현장 자료를 바탕으로 우리나라의 산림에 적합한 식생 지수와 개엽일 추출 방법을 비교 및 평가하였다. 또한 각 픽셀의 임상 구성비율과 식생지수로부터의 식물계절 추출 간의 상관관계를 분석하였다. 국립수목원에서 실측된 개엽일을 바탕으로 NDVI와 EVI를 비교하였을 때 우리나라의 산림에는 EVI가 더 적합하였고, 식생지수 계절 곡선으로는 일차도함수가 개엽일 추출에 가장 적합하였다. 현장 개엽일이 조사된 7개의 활엽수종이 속한 픽셀과(중앙픽셀) 이웃한 8개의 픽셀을 실측 개엽일과 비교하였을 때, 실측과 가장 높은 상관관계를 보이는 베스트 픽셀을 구성하는 임상은 활엽수림 67.9%, 침엽수림 14.3%, 혼효림 7.7% 구성을 보였고 평균 결정계수($R^2$)는 0.64으로 나타나 높은 상관관계를 보였다. MODIS EVI로 추출한 전국 평균 개엽일은 최소 112.9일(2014년)에서 최대 129.1일(2010년)으로 2003년 이래로 약 0.16일 빨라졌다. 향후 연구에서는 활엽수림과 혼효림의 개엽일 분석 외에도 다양한 기후 및 지형 인자들의 영향을 파악하기 위해 침엽수림의 개엽일 추출까지 확장할 필요가 있으며 우리나라 산림 생태계의 다양성을 고려한 종합적인 연구가 진행되어야 할 것이다.

1. 서론

식물계절학은 기후 요인에 따른 개엽, 개화, 성숙, 단풍, 낙엽 등과 같은 식물의 생활사의 변화를 관측 및 연구하는 학문으로 장기(long-term) 식물계절 변화는 각 생태계에 대한 기후변화의 반응을 모니터링 하는 지표가 된다(Hmimina et al., 2013; Ulsig et al., 2017). 특히 식물 계절에서 개엽과 낙엽으로 결정되는 식생 성장 기간은 잠열과 가시열 등 에너지 순환과 기상 조건의 강한 영향을 받는 동시에 탄소와 물 순환의 지리적, 계절적 변화에 직접적인 영향을 준다. 이는 생태계의 연간 생산성 즉 연 광합성량과 관계되는 중요한 지표로서 지역적 규모부터 전세계적 규모까지 다양한 공간 스케일과 생태계에서 연구되어 왔다(Schwartz and Marotz, 1986; White et al., 1997). 나아가 개엽과 낙엽으로 대표되는 식생 생장은 기상 정보를 바탕으로 실행되는 지역적 혹은 전지구적 생태계 예측 모델, 생지화학 모델들의 주요 계수로 사용되었다(Moulin et al., 1997; Zhao and Schwartz, 2003;Beck et al., 2006; Bradley et al., 2007). 그러나 현장에서 관측된 개엽일과 낙엽일 등의 식물계절 자료는 지역적 혹은 전지구적 규모의 변화를 관측 및 예측하기에 인력, 경제적 및 공간적 제약으로 어려움이 있다. 이러한 제한점을 극복하기 위해 인공위성 자료를 통한 개엽 및 낙엽일 추정이 지난 삼십여년간 지속적으로 발전되어 왔다(Pettorelli et al., 2005a). 중저해상도(~250m<) 인공위성자료는 1970년대 후반부터 발달하기 시작하여 1990년대에는 식물계절연구에 보다 활발히 사용되기 시작하였고 대표적인 인공위성 센서로는 Advanced Very High Resolution Radiometer(AVHRR)과 Moderate Resolution Imaging Spectroradiometer(MODIS) 등이 있다(Justice etal., 1985; Huete et al., 1994). 특히 MODIS는 2000년부터 1일~16일 간격으로 조합된 지상 자료를 확보하고 있어 장기 식생 변화에 대한 기초 자료가 되고 있다. Landsat 8(30 m, 2013년 2월 자료 획득 시작), Sentinel-2(10-20m, 2015년 6월 자료 획득 시작) 등 식생 모니터링 연구에 사용 가능한 고해상도(<~30m) 영상 역시 최근에 식물계절 연구에 점차 쓰이기 시작하고 있다(Eklundh et al., 2012). 그러나 아직 고해상도 식생 영상은 수집된 기간이 짧고 자료가 유료이거나 공개된 자료가 제한적인 경우가 있어 장기적인 계절적 변화가 중요한 산림 식생의 경우, 시간 해상도가 높고 관측 기간이 긴 MODIS 자료의 병행 연구가 필요하다.

인공위성영상 자료로부터 식물계절을 추출하기 위해 사용되는 식생지수(Vegetation indices)들은 적색(Red) 파장과 근적외선(ner ifra-red, NIR) 파장의 비율로 엽록소가 적색 파장(620-670 nm)을 흡수하고 엽육세포가 근적외선 파장(841-876 nm)을 반사하는 원리를 이용하였다(Pettorelli et al., 2005a). 식생지수는 식물계절과 이와 관련된 다른 생태학적 요인들 즉 식생의 시간적 공간적 변화, 동식물서식지 감소 및 파편화, 가뭄, 홍수, 산불 등의 기후 변화 재해들과 상관성이 높아 다방면으로 활용되고 있다(Beck et al., 2006). 그 중 Normalized Difference Vegetation Index(NDVI)와 Enhanced Vegetation Index(EVI)는 가장 효과적으로 널리 사용되는 검증된 식생지수들이다. NDVI는 엽록소량을 효과적으로 탐지하기 때문에(Rouse et al., 1974; Myneni et al., 1995a), 엽면적지수(Xiao et al., 2002), 식물계절 모니터링(Lee et al., 2002;Vina et al., 2004), 식생 모델 검증(White and Nemani, 2003;Peng et al., 2012), 식생커버(Jiang et al., 2006), 다양한 식생에서의 생산량(Lee et al., 2002; Zhang et al., 2003; Wang etal., 2005; Sakamoto et al., 2005; Zhang et al., 2005; Soudani et al., 2012) 등과 강한 상관관계를 보인다. EVI는 Huete et al.(1999)이 개발한 지수로 식생의 수관(vegetation canopy) 반사 민감도를 개선하고 대기 에어로졸에 대한 민감도를 높이는 청색 파장(459-479 nm) 밴드와 식생 수관 적응 팩터(canopy background adjustment factor, L)를 추가하여 계산하였다(Huete et al., 1997 & 2002; Wardlow and Egbert, 2010). EVI는 NDVI와 비교해 대기 잔류물에 의한 에러를 줄여주는 지수로 식생이 빽빽한 산림 지역에서 산림식물계절과 프로세스 기반 모델에서 보다 효과적으로 사용되었다(Huete et al., 1994; Zhang et al.,2003; Xiao et al., 2004b).

인공위성영상 자료로부터 나온 식생지수들을 사용하는데 있어 일차적으로 문제가 되는 것은 센서 자체의 오류, 구름이나 대기 에어로졸로 인한 에러 값이 자료에 포함되어 자료의 품질이 저하되는 것이다. 때문에 인공위성영상 자료를 사용하기 전에 반드시 저품질 자료를 제거하고 여전히 남아있는 에러 값을 적절한 평활 곡선 방법을 적용(curve fitting)하여 시계열 자료로 생산하는 과정이 필요하다(Bradley et al., 2007; Hmimina et al.,2013; Lee et al., 2015). 다양한 평활 곡선 방법 기술이 시도되어 왔고 대표적인 방법에는 Savitzky-Golay 필터(Chen et al., 2004; Jönson ad Eklundh, 2004), Gaussian 함수(Jönsson and Eklundh, 2002), wavelet 변환(Sakamoto et al., 2005), logistic 함수(Zhang et al., 2003; Beck et al., 2006), Fourier 변환 (Menenti et al., 1993; Verhoef et al., 1996; Roerink et al., 2000; Moody and Johnson, 2001) 등이 있다. Atkinson et al.(2012)은 인도에서 해상도 ~4.6 km 지상 엽록소 지수 자료를 이용한 NDVI로 침엽수림과 활엽수림에서 식물계절을 추출하였는데, 활엽수림의 경우 double logistic 함수와 Fourier 변환이 식물계절 패턴을 비교적 잘 모사하고 침엽수림에서는 Whittaker 필터(Eilers, 2003)가 잘 모사하였다. Beck et al.(2006)은 MODIS NDVI로 겨울 NDVI 고려한 이중 로지스틱 방법으로 평활 곡선을 그리고 변곡점을 개엽일과 낙엽일로 추출하여 높은 상관성을 보여주었다. MODIS NDVI로 침엽수림에서 식물계절을 연구한 Ulsig et al.(2017)은 Gaussian과 logistic 함수의 한 케이스인 이중 S자형 함수가 식물계절 추이를 잘 모사함을 보여주었다. Sakamoto et al.(2005)은 wavelet 변환을 사용하여 벼 작물의 식물계절을 EVI로 추출하였고 Xiao et al.(2004a)은 SPOT-4와 MODIS EVI로 열대상록수림의 식물계절의 패턴과 최대값을 비교하였다. 이렇듯 수많은 선행 연구들이 진행 되어왔음에도 여전히 인공위성영상의 식물계절과 현장 자료와의 검증은 활발히 이루어 지지 않고 있는데 이는 두 자료 간의 시공간적 스케일이 크게 다르기 때문이다(Liang et al., 2011). 우리나라의 경우, Kim et al.(2014)이 TIMESAT 알고리즘(Jönsson and Eklundh, 2004)을 활용하여 MODIS의 EVI 시계열자료로 남한 지역 전체의 식물 계절을 현장 측정된 개엽일과 비교하고 공간적으로 설명하였다. 그러나 사용된 현장 자료의 수가 제한적이고 산림과 비산림 그리고 임상을 고려하지 않아 산림생태계에서 식물계절 추이가 서로 다른 침엽수림과 활엽수림 및 혼효림이 구분되지 않았다. 산림 생태계의 경우, 임상별(침엽수림, 활엽수림 및 혼효림)로 식물계절 변화, 특히 개엽일과 낙엽일은 서로 다른 계절적 양상을 가진다(Liang and Schwartz, 2009). 따라서 산림의 임상 구분은 산림생태계의 식물계절을 지역적 규모에서 정확하게 파악하는데 있어 직면해야 하는 과정이다. 또한 우리나라의 산림은 복잡하고 작은 넓이에서 다양한 형태로 패치화 되어있어 인공위성 영상에서 제공하는 저해상도(1 km * 1 km) 지형 지도로는 표현되기 어렵다. 본 연구에서 사용된 MODIS의 하나의 픽셀(1 km* 1 km)에도 세 가지 임상이 비슷한 비율로 섞여있거나 50% 미만의 비율로 산림이 존재하기도 하는 등 각각의 픽셀이 하나의 임상으로 대표되기가 어려운 경우가 많았다.

이에 본 연구에서는 1) 우리나라 전국단위의 실측된 대축척임상도(1:5,000)를 MODIS 픽셀 해상도로 재구성하고 2) MODIS의 NDVI와 EVI 두 식생지수들로부터 시계열 식물계절 곡선을 생산하여 3) 현장 조사된 개엽일과 비교 분석하였다. 식생지수의 식물계절 곡선에서는 개엽일 및 낙엽일 모두 추출 가능하나 낙엽일은 단풍 시작 시기, 단풍 절정 시기, 50% 이상 잎이 떨어진 시기, 전체 잎이 떨어진 시기로 구분이 다양하고 수종마다 낙엽의 색깔, 지는 시기 등이 하나의 날짜로 정하기에 모호하여 이 연구에서는 개엽일만 비교 분석에 사용하였다. 또한 4) 각 픽셀의 임상 구성비율이 식생지수에서의 식물계절 추출과 어떤 상관관계가 있는지 분석하고 그 활용 방안에 대해 고찰하였다. 이 연구에서는 NDVI와EVI 두 가지 식생 지수를 모두 사용하여 분석하였고 어느 식생지수가 우리나라 산림에 더 적합한지 대해 토의하였다. 시계열식물계절 곡선은 Beck et al.(2006)의 double logistic 함수를 각 픽셀마다 적용하였고 개엽일과 낙엽일을 추출하는 방법은 사용자 지정 임계치를 주는 방법과 1차 도함수 방법을 적용하여 실측과 보다 일치하는 식물계절 추출 방법을 제안하고자 하였다.

2. 방법

1) 현장조사

국립수목원의 식물계절 조사는 2009년부터 시행되었고 각 지방 공립수목원 중심으로 5개 권역(강원, 경기 북부, 경기충청, 전라, 경상, 제주권)으로 실시되었다(Fig. 1). 개엽일의 기준은 식물의 최소 세 곳에서 완전히 펴진 잎이 나올 때 즉 잎자루가 완전히 보일 때를 측정하여 개엽일로 하였다(Korea Forest Service and Korea National Arboretum, 2010). 식물계절 조사는 산림청과 국립수목원 ‘기후변화 취약 산림식물종 보전·적응 사업’에 의해 2009년부터 시작되었으며 본격적인 자료 수집이 이루어진 2010년부터 자료를 사용하였다. 본 연구에서 이용한 수종들은 국가산림자원정보(NFI, National Forest Inventory)를 바탕으로 선정된 대표 활엽수종 중 국립수목원에 의해 식물계절이 조사된 수종 7종(아까시나무, 신갈나무, 나도밤나무, 서어나무, 상수리나무, 졸참나무, 고로쇠나무)으로 2010년부터 2016년까지의 조사 자료 중 3년 이상 자료가 존재하는 16개 지점의 자료만 사용하였다(Fig. 1).

Fig. 1. (a) Location of Start of Season (SOS) collected by Korea National Arboretum and (b) study period of the selected field.

2) 임상도 재구성

우리 나라 전국토 산림의 공간정보를 자하게 조사한 임상도는 임상, 수종, 경, 영급 등과 같은 정보를 제공하고 있다. 1:5,000 축척의 수치 임상도는 디지털 항공 영상을 활용하여 복잡 다양한 산림을 판독·구획하여 제작되었다(Jo et al., 2012). 이는 근래 주로 사용되는 고해상도 및 중저해상도 토지피복 인공위성 영상 Landsat 7ETM, CCI V2, GlobeLand30 등과 비교해서 정밀도가 높고 복잡한 우리나라의 지형 환경을 잘 반영하고 있다. 그러나 벡터 형태의 데이터로 만들어져 자료의 용량이 크고 복잡해 전국 단위를 분석하는데 기술적 한계가 있어 전국을 여러 단위로 나누어 분석하였다.

먼저 ArcGIS를 사용하여 각각의 폴리곤으로 되어있는 임분구획을 같은 임상별로 통일하여 합치는 과정(dissolve와 merge)을 거친 후, 침엽수림, 활엽수림, 혼효림만 추출하여 사용하고 도시지역과 나지 등 비산림 지역은 제거하였다. 산림 지역의 각기 다른 크기의 폴리곤을 MODIS 픽셀 크기에 맞추어 분할(intersect)을 하였다(Fig. 2). 좌표는 람베르트 정각원추도법(Lambert Conformal Conic(LCC) WGS84)으로 변환하였다. MODIS의 원 좌표는 Sinusoidal 좌표로 이는 정현곡선도법으로 지구 중심부는 비교적 정확하나 주변부로 갈수록 왜곡이 심해져 사용하기가 적합하지 않다. 또한 좌표를 변환하는 과정을 거치게 되면 그 과정에서 픽셀의 왜곡이 불가피한데 이를 최소화 하기 위해 중위도에서 지역왜곡오차가 최소화 되는 도법인 LCC 좌표로 MODIS 좌표와 임상도의 좌표를 일치시켰다(ArcGIS Help center). MODIS 지상 관측 밴드 자료들도 ‘clear-sky’나 구름, 구름 그림자 등의 자료를 만들 때 Sinusoidal 좌표의 문제점을 인식하고 LCC 좌표로 재조합하여 배포하기도 한다(Luo et al., 2008).이와 같은 과정을 거쳐 임상도는 MODIS 1 km * 1 km 픽셀과 같은 크기의 픽셀로 변환하고 1 km* 1 km 해상도 내에서 활엽수림, 혼효림, 침엽수림이 차지하는 비율을 나타내는 지도로 재구성되었다(Fig. 2).

Fig. 2. Re-construction of the forest type map (a) 1:5,000 forest type digital map, (b) dissolved by forest type, (c) merged by forest types, and (d) intersected based on MODIS pixel size.

현장에서 조사된 개엽자료는 포인트 자료이나 포인트가 속한 임상도의 픽셀이 해당 임상도를 대표한다고 가정하고 해당 픽셀과 일치하는 MODIS 픽셀의 식생지수와 비교하였다. 그러나 픽셀 내에서 개엽일이 조사된 포인트의 공간적인 위치에 따라 해당 픽셀의 값을 대표하지 못할 수도 있었다. 재구성된 임상도는 임분구획이 아닌 MODIS 픽셀과 같은 정각형의 픽셀로 나누어져 있기 때문에 조사된 지점이 셀 내에서 가장자리에 위치하거나 이웃한 픽셀의 임분 구역의 연장선상에 소속되어 있는 등의 경우가 발생할 경우 이웃한 픽셀의 특성을 반영하기도 하였다. 때문에 개엽일 조사 지점이 속하는 픽셀과 이웃하는 8개 픽셀의 개엽일을 모두 추출하여 비교하였고 9개의 픽셀 중 현장 개엽일 자료와 상관성이 가장 높은 픽셀을 ‘베스트 픽셀로 선정하여 그의 임상 구성을 알아보았다.

3) 인공위성영상 처리

이 연구에서 사용된 인공위성영상은 NASA의 지구 관측계획의 일환으로 개발된 Terra와 Aqua 위성의 다중 분광센서(Moderate Resolution Imaging Spectroradiometer, MODIS)의 자료이다. MODIS의 산출물 중 해상도 1 km* 1 km, 16일 단위의 MOD13A2(Terra)와 MYDA2(Auqa)식생지수를 사용하였다. MODIS 산출물의 두 가지 버전 중 collection-6를 사용하였는데 이것은 기존의 재합성 알고리즘을 개선하여 기존에 사용되던 collection-5보다 정확도를 높인 자료이다(Jang et al., 2014). Terra는 DOY 1, 17, 33, …, Aqua는 DOY 9, 25, 41, …의 16일 단위로 계산된 산출물이 나오는데 각 위성의 산출물이 서로 8일씩 차이가 나도록 되어있어 두 위성의 영상을 결합하여 8일 단위(DOY 1, 9, 17, 25, …) 자료로 구성하였다. 사용된 식생지수는 NDVI와 EVI는 다음 식과 같이 계산된다.

\(\mathrm{NDVI}=\frac{\rho_{n i r}-\rho_{r e d}}{\rho_{n i r}+\rho_{r e d}}\)       (1)

\(\mathrm{EVI}=2.5 * \frac{\rho_{n i r}-\rho_{r e d}}{\rho_{n i r}+\left(C_{1} * \rho_{r e d}-C_{2} * \rho_{b l u e}\right)+L}\)       (2)

ρ_nir, ρ_red, ρ_blue는 각각 적외선, 적색, 청색 밴드의 반사도, L은 토양 적응 계수, C₁과 C₂는 대기 산란 보정 계수로 일반적으로 L=1, C₁=6.0, C₂=7.5를 사용한다.

인공위성 영상은 필연적으로 기계적인 복잡성과 대기 조건에 따른 한계성을 가지는데 이는 센서의 해상도, 센서 및 기기 에러, 구름 영향, 대 기후, 지상 상태 등에 의한 에러가 존재하기 때문이다. 이 에러를 감소하기 위해 일차적인 전처리 방법으로 ‘수정된 최대값 중첩법(CV-MVC, Modified the Constrained View angle Maximum Value Composite)’을 적용하여 처리된 영상을 사용하였다. 전처리를 하였음에도 존재하는 노이즈를 제거하기 위해 평활 곡선 작업을 거쳐서 시계열 식물계절 곡선을 생산하였다(Fig. 3, curve fitting). 평활 곡선 방법은 double logistic 함수를 사용하여 각각의 픽셀마다 적용하였다(Fig. 4). 함수 VI(t)는 시간(t)에 따라 겨울 식생지수(wVI), 최대 생지수(mVI), 변화율(증감(mS)과 감소(mA)), 시작 시점(S)과 감소 시점(A)을 계수로 사용하며 시즌 중간에 갑작스러운 떨어짐은 반영하지 않고 무시하도록 하여 식생지수의 과대 혹은 과소 평가를 피한다(식 3).

\(\begin{array}{r} V I(t)=w V I+(m V I-w V I) * \left(\frac{1}{1+\exp (-m S *(t-S))}+\right. \left.\frac{1}{1+\exp (m A *(t-A))}-1\right) \end{array}\)       (3)

Fig. 3. The work-flow chart of processing MODIS data, forest type map, and extraction of start of season (SOS).

Fig. 4. The methods used to extract phenological thresholds a) based on a user-defined threshold when the 40% or 50% of green chromatic coordinate (Gcc) reached (b) phenophases based on local extremes in the first derivative f(t).

곡선의 평활화에 사용되는 함수들은 대부분 각값별로 가중치를 주는데, 만약 겨울 기간의 마이너스 값들도 포함하는 곡선에서 가중치를 적용하면 마이너스 값 때문에 연 초 곡선이 과대 평가되고 이 때문에 다음 구간 커브에서 가중치가 적어져 과소평가 되는 경향을 보이게 된다. Double logistic 함수에서는 겨울 기간 동안의 겨울 값을 고정시켜 연 초와 연 말의 과대 및 과소 평가를 피하도록 하고 곡선이 증가하는 시기(성장시즌이 시작되는 시점)에서 발생하는 값들에 가중치를 더 주어서 봄철과 가을철의 큰 폭의 증가 혹은 감소가 반영되도록 하였다. 고정 겨울은 각각의 연도에서 가장 낮은 식생지수 값으로 하거나 만약 그 값이 0보다 낮을 값일 경우에는 0 사용한다. Fourier 함수보다는 전체 성장기간에 걸쳐 더 나은 결과를 보이고 TIMESAT 알고리즘의 Gaussian과 비교 했을 때는 거의 비슷하나 곡선이 막 시작되거나 끝나는 지점에서의 값을 더 잘 모사하였다. 더 자세한 설명은 Beck et al.(2006)의 연구를 참조하기로 한다.

4) 개엽일 추출 및 비교

Filippa et al.(2016)이 제안한 식물계절 추출의 방법 중 사용자 지정 임계치(trs, thresholds in a seasonal green chromatic coordinate trajectory)와 일차도함수(derivatives,local extreme in the first derivative)를 사용하였다. 식생지수의 식생계절곡선으로부터 식물계절을 추출할 때는 각 이미지의 전체 밝기에 대한 녹색 정도(green chromatic coordinate, G_cc)를 수치화 하여 각각의 방법에 맞는 식물 계절 시점(개엽일 혹은 낙엽일)을 추출한다. 사용자 지정 임계치 방법을 적용할 때 임계치는 White et al.(1997)이 제안한 준을 참고로 하였다. White et al.(1997)은 연식생지수의 최대값과 최소값 차이의 비가 40%(개엽시작 DOY 111)에서 50%(DOY 127)까지 변화시키며 그 결과를 비교하였는데 산림지역에서는 평균 50% 이상일 때 개엽이 시작된다고 보았다. 이 연구에서는 우리나라의 개엽일 조사자료를 바탕으로 우리나라 산림의 실정에 맞게 40%와 50% 두 가지 모두 사용하여 개엽일을 추출하여 비교하였다. 일차도함수 방법은 식생지수의 식생계절곡선을 미분하여 변곡점을 추출하는 것이다(Fig. 4). 위에서 설명한 세 가지 방법의 식물계절의 개엽일 추출 과정은 R 프로그램의 패키지 ‘phenopix’ (http://r-forge.r-project.org/projects/phenopix/)를 사용하였다.

MODIS 식생지수에서 추출된 개엽일과 실측 개엽일을 비교할 때 사용된 통계 지수들은 결정계수(Coefficient of determination, R²), 평균 제곱근 오차(Root Mean SquareSrror, RMSE), AIC(Akaike InformationCriterion)등이있다. AIC는 특정 집단 내에서 형성된 통계 모델의 상대적인 적합도를 보여주는 것으로 이 연구에서는 9개의 픽셀(조사지점이 속한 픽셀과 주변 8개)이 한 집단으로 취급될 수 있으므로 집단 내에서의 상대적 평가를 하기 위한 지수 중 하나로 사용하였다. 때문에 다른 집단과의 AIC비교는 의미가 없으나 같은 실측 자료를 공유하는 9개 픽셀 집단 내에서만 유용하므로 참고 자료로만 사용하고 결과에서는 제외하였다.

3. 결과

1) 식생지수 및 개엽일 추출방법 비교

우리나라 산림의 개엽일과 상관성 높은 식생 지수와 개엽일 추출 방법을 선택하기 위해개엽일 조사 지점별로 결정계수(R²)를 비교하였다. Beck et al.(2006) 곡선이 적용된 NDVI와 EVI의 계절 곡선에 세 가지 개엽일 추출 방법(사용자지정 임계치 40% 와 50%, 일차도함수)으로 2010년부터 2016년까지의 개엽일을 추출하고 각 실측 지점별로 실측 개엽일과 식생지수에서 추출된 개엽일 간의 결정계수를 비교하였다(Fig. 5). NDVI를 사용하여 추출한 개엽일의 경우, R²는 사용자지정 임계치 40%일 때 0.19, 50%일 때 0.21, 일차도함수일 때 0.35, EVI사용한 경우 각각 0.24, 0.20, 0.46이었다. 두 식생지수를 비교했을 때 측정지점 ‘정도리(Jeongdori)’를 제외하고 EVI에서 추출된 개엽일과 실측 개엽일과의 상관성이 더 높음을 볼 수 있었다. EVI를 기반으로 개엽일 추출 방법 세 가지를 비교해보면 16개 측정지점 중 6개 즉, 내장산(Naejangsan), 주왕산(Juwangsan), 정도리(Jeongdori), 한라 1100도로(Halla1100), 교래(Gyorae), 가야산(Gayasan)을 제외하고 일차함수가 높은 상관성을 보였다. 이러한 결과를 바탕으로 EVI 계절 곡선에서 일차도함수로 추출한 개엽일이 가장 높은 신뢰도를 가진다고 판단하고 이웃한 8개 픽셀 모두 분석을 수행하였다.

2) ‘베스트 픽셀’ 선정 평가

실제 측정 포인트가 대표할 수 있는 임분 구성을 가진 픽셀(베스트 픽셀)을 선정하기 위해 실측 개엽일 지점이 속한 픽셀(중앙 픽셀)과 이웃 픽셀 8개의 개엽일을 모두 추출하였다. 중앙 픽셀의 EVI에서 추출된 개엽일이 실측 개엽일을 잘 모사하는 경우도 있으나 실측된 개엽일 포인트가 픽셀 내에서 위치에 따라(픽셀의 가장자리 이거나 이웃한 셀의 연장선에 있는 임분 구성에 속한 다거나 하는 경우 등) 주변 픽셀이 좀 더 높은 상관 관계를 보이기도 했다. 따라서 주변 8개 픽셀까지 고려했을 경우 가장 높은 상관관계를 가지는 픽셀(베스트 픽셀)을 각 픽셀별로 모두 시계열 그래프로 그려 결정계수(R²)와 RMSE, AIC를 고려하여 판단하였다(Table 1). 베스트 픽셀을 정할 때 절대값의 차이도 중요하나 장기간 기후변화에 따른 연변화도 고려해야 하므로 1차적으로 R²를 기준으로 두었으나 RMSE나 AIC 값이 결정계수 간 차이보다 과하게 과소평가 혹은 과대평가 되는 경우가 발생하므로 시계열 경향을 함께 고려하였다. R²가 높다해도 RMSE차이가 2주 이상 나는 경우에는 RMSE가 낮은 픽셀을 택하였다. 이 때 기준을 2주로 한 것은 실측자료의 수집 기간이 1-2주 간격이기 때문이다. 예를 들면 월출산의 경우 졸참나무(Quercus serrata)와 나도밤나무(Meliosma myriantha)의 개엽일이 조사되었는데 중앙 픽셀에 이웃한 두 픽셀에서 나도밤나무 개엽일과의 R²는 0.78과 0.65, RMSE는 각각 26.8일과 1.7일로 R2가 약 17% 상승 하였으나 RMSE에서 약 93% 감소하여 그 차이가 크게 나타났다. 두 픽셀에서의 개엽일과 실측 개엽일과의 연간 패턴을 고려해보니 R2가 0.78, RMSE가 26.8일인 픽셀의(High R-squared and RMSE, Fig. 6) 개엽일이 최대 53일까지 차이 나는 것을 확인하였다. 이럴 경우에는 비록 결정계수가 약간 낮으나 RMSE가 유의미한 수준이고 연간 추세를 보았을 때 실측 개엽일의 추세를 더 잘 모사한다고 판단되는 픽셀을 베스트 픽셀로 선택했다.

Fig. 6. Comparison of field observation and neighbor pixels in Wolchulsan. Black solid line with filled diamond: Start of Season (SOS) from field survey, gray solid line with unfilled rectangular: a pixel with high R-squared (0.78) and RMSE (2.8), gray dotted line with unfilled rectangular: a pixel with relatively low R-squared (0.65) and RMSE (1.7).

Table 1. Statistics of correlation analysis of the start of season (SOS) for the best pixels between EVI and field observation. Pixel location indicates the pixel number including measurement point (0: center) and its neighbors pixel number from the center pixel (1001: lower right, 1000: lower, 1001: lower right, 999: lower left, 1: right, -1: left, -999: upper right, and -1000: upper left, respectively). R²: coefficient of determination, RMSE: root mean square error, DBF: deciduous broadleaf forest, CNF:coniferous needle-leaf forest, MF: mixed forest

베스트 픽셀의 평균 R²는 0.64(0.07–0.99), RMSE 평균 6.6일(1.27–15.15일) 전체 비교 지역의 85% 이상에 해당하는 9개 지점은 실측 개엽일과 MODIS EVI에서 추출된 개엽일의 상관성 높음으로 판단되었다. 실측 개엽일이 7년(2010년-2016년) 미만인 왕산, 가야산, 교래 지점을 제외하면 평균 R²는 0.75, 평균 RMSE 5.4일로 향상되었으며 실측 자료가 부족한 곳의 상관성이 비교적 낮음을 볼 수 있다. R²가 0.4 미만인 지점은 5개, 주왕산(Juwangsan), 가야산(Gayasan), 내장산(Naejangsan), 교래 곶자왈(Gyorae), 한라 1100도로(Halla1100)였고, 그 중 상관성이 매우 낮은(R² < 0.2) 지점은 2개(교래 곶자왈과 내장산)였다. 이 연구에서 사용된 16개의 실측 지점 중비산림이나 산림으로 구분되는 면적이 픽셀의 50% 미만인 2개 지점을 제외하고 중앙 픽셀의 평균 임상 구성은 활엽수림 48.0%, 침엽수림 22.3%, 혼효림 7.2%로 구성되어 있었다. 베스트 픽셀의 평균 임상 구성은 활엽수림 67.9%, 침엽수림 14.3%, 혼효림 7.7% 구성을 보였다. 이 중 실측 기간이 7년 미만인 지점을 제외하면 활엽수림 70.9%, 침엽수림 13.8%, 혼효림 5.5% 구성을 보였다. 재구성된 임상도의 경우 벡터형태의 자료를 래스터 형태로 변환시킬 때 임상도에서 임상이 겹쳐지는 부분들이 중복되 계산된 경우가 있어 몇몇 픽셀들의 임상 구성의 합이 100%를 넘는 경우가 있으나 통계 분석에 영향을 줄 정도가 아니므로 이는 무시하였다.

3) 개엽일 지도 산출

베스트 픽셀을 기준으로 MODIS EVI에서 추출된 개엽일은 각 픽셀의 임상 구성이 활엽수림 50% 이상 포함할 경우 높은 상관을 보였다. 활엽수림이 픽셀의 50% 이상인 픽셀만 선택하여 2003년부터 2016년까지의 전국 개엽일을 추출하였다. MODIS EVI로 추출한 전국 평균 개엽일은 최소 112.9일(2014년)에서 최대 129.1일(2010년)으로 2003년 이래로 약 0.16일 빨라졌으나 유의미하게 감소한다고 보기는 어렵다(Fig. 7). MODIS EVI로 추출한 개엽일과 실측 개엽일과는 평균 최소 0.7일(2010년)에서 최대 2.8일(2015년) 차이가 났다.

Fig. 7. Averaged start of season (SOS) using MODIS EVI from 2003 to 2016 at the forest area in South Korea

봄철 기상 조건과 개엽일과의 관계를 보기 위해 기상청에서 제공하는 봄철(3월–5월) 기상 특징을 바탕으로 전년, 후년 혹은 평년 대비 변화가 뚜렷한 해를 중심으로 전국 산림의 개엽일을 Fig. 8에 제시하였다. 2008년과 2014년은 봄철 강수량이 낮았고 2010년 봄철 평균기온은 1973년 이후 가장 낮았다. 반면 2014년과 2016년은 봄철 평균기온이 높았다(Table 2). 해안 지역 및 칩엽수림, 도시 지역을 제외한 산림 지역에서 2008년은 110일부터 140일, 2010년은 130일부터 150일, 2014년은 100일 130일, 2016년은 100일터 130일에 주로 발생하였다. 강수량과 평균 기온만 개엽일과 비교해 보았을 때는 강수량이 적고 평균기온이 높으면 개엽일이 빨라지고 강수량이 높으나 평균기온이 낮을 경우 개엽일이 늦어지는 경향을 보였다. 지역적 차이를 보면 강원 지방과 경북 지방이 태백산맥을 따라 개엽이 늦게 시작되고 산맥의 동쪽 해안을 따라서는 산맥 서쪽보다 개엽이 빠르게 나타났다. 서해 남부 해안 지역과 내륙 지방은 비교적 빠르게 개엽이 시작되나 남부 내륙의 지리산과 덕유산 지역은 비교적 주변 내륙 지방에 비해 개엽이 늦게 나타났다(Fig. 8).

Fig. 8. Start of season (SOS) in deciduous and mixed forests for the years indicated.

Table 2. Mean air temperature and precipitation during the spring period (March to May) from 30 years normal and selected years

4. 토의 및 결론

이 연구에서는 MODIS의 장기 NDVI와 EVI를 이용하여 복잡 다양한 우리나라 산림 지역을 더 잘 모사하는 식생지수와 개엽일을 추출하는 방법을 평가해 보았다. 국립수목원에서 실측된 개엽일을 바탕으로 NDVI와 EVI를 비교해 보았을 때 우리나라의 산림에는 EVI 가 더 적합하다고 판단하였다(ig. 5). 또한 식생지수 계절 곡선을 위해 평활 함수를 적용했을 때 일차도함수가 개엽일 추출에 가장 적합한 것으로 판단되었다. 이 연구에서 사용된 Beck et al.(2006)의 double logistic 평활 알고리즘은 고위도 지방에서 테스트 하였으나 double logistic 알고리즘은 이미 다양한 식생과 지역에서 검증된 방법이다. 이 알고리즘은 항내성을 가질 수 있게 고정된 윈도우 내에서 자료를 평활화 하므로 복잡 지형에서의 식생계절의 특정 이벤트(개엽, 낙엽, 최대 성장점 등)을 추출하는 방법에 있어 높은 정확도를 가진다(Caiet al., 2017). 지금까지 이루어진 선행 연구들을 보면(Kim et al., 2008; Kim et al., 2014) 현장에서 관측된 개엽일 자료가 많지 않으므로 식생지수와 개엽일 추출 방법을 시도해 보더라도 그 적합성을 비교 판단 하기에는 한계가 있었으나 이 연구에서는 비교적 풍부한 현장 자료를 바탕으로 산림 식물 계절에 적합한 방법을 고찰 할 수 있었다.

인공위성 영상의 식생지수를 사용할 때 발생하는 공간 스케일의 불일치는 모든 연구자들에게 피할 수 없는 불확실성을 감수하게 하였는데 우리나라와 같이 산림의 형태가 다양하고 복잡한 지역에서는 그 불확실성이 더 클수 밖에 없다. 이미 30여년간 인공위성의 발전으로 기기 결함, 구름 제거 등의 에러 제거 대한 기술적 발전이 높은 수준으로 올랐으나 복잡 지형에 대한 정확도는 그 지역의 자세한 현장 자료의 유무 의해 결정되므로 개선하기 위해서는 많은 인력과 시간이 필요하였다. 또한 자료가 구축되었다 하더라도 인공위성 영상과의 공간 스케일이 맞지 않아 다운스케일링 하거나 업스케일링 하여 평균된 값을 사용하였는데 이 방법은 해당 연구지역의 대표값을 가지기는 하나 다양한 식생에 대한 값으로는 보기 어려웠다. 전국 단위의 산림 생태계 연구를 수행할 때 인공위성 영상에 맞게 임상이 평균화 되거나 간소화되는 방식이 사용되었다(Kang et al., 2003; Kim et al., 2008). 이 연구는 이러한 불확실성을 줄이고자 인공위성 영상의 픽셀 크기에 맞게 대축척(1:5,000) 임상도를 재구성하되 대표되는 임상으로 평균 되어버린 값이 아닌 그 픽셀을 구성하는 임상들의 구성 비율 정보를 가지게 하였다. 이를 이용하여 산림과 비산림 지역을 구분하고 침엽수림과 활엽수림 및 혼효림의 구성비 정보를 포함한 전국 산림 개엽일 지도를 제작하였다(Fig. 8). 또한 실측 개엽일 포인트의 픽셀 내부 위치에 따라 이웃한 픽셀과 더 상관관계가 높을 수 있음을 가정하여 실측 자료가 위치한 픽셀 및 이웃한 8개 픽셀의 개엽일을 모두 비교해 보았을 때 픽셀의 평균 임상 구성이 활엽수림 67.9%, 침엽수림 14.3%, 혼효림 7.7% 구성을 보였다. 이 중 실측 기간이 7년 미만인 지점을 제외하면 활엽수림 70.9%, 침엽수림 13.8%, 혼효림 5.5% 구성을 보였다. 이를 바탕으로 활엽수림이 50% 이상인 픽셀에서의 EVI 개엽일 값이 보다 유의미하다고 해석 수 있다.

결정계수와 RMSE를 기준으로 판단하였을 때 실측개엽일과 MODIS EVI에서 추출할 개엽일 간의 상관성이 낮은 5개 지점들 즉, 주왕산(Juwangsan), 가야산(Gayasan), 내장산(Naejangsan), 교래 곶자왈(Gyorae), 한라 1100도로(Halla1100)는 앞서 사용자 지정 임계치와 일차도함수로 개엽일 추출의 적합성을 판단할 때 일차도함수 방법보다 사용자지정 임계치 방법의 개엽일 추출이 더 적합한 지역이었다. 따라서 개엽일을 추출하는 방법을 한가지로 하기보다는 지역에 따라 다르게 적용하는 방안을 고려해 보아야 할 것으로 판단된다. 그러나 지역에 따라 다르게 적용하기 위해서는 지역마다 그기준이 명확해야 하며 이 부분은 향후 반드시 해결해야 하는 과제이다. 또한 주왕산, 가야산과 교래 곶자왈의 경우, 개엽일 조사년도가 다른 지점보다 2년 혹은 3년 적었다. 이 연구에 사용된 총 조사연도가 7년이므로 2-3년 자료 누락은 전체의 약 30%로 통계 분석에 영향을 줄 수 있는 것으로 보인다. 내장산의 베스트 픽셀은 활엽 수림 31.0%, 침엽수림 40.0%, 혼효림 3.0%로 침엽수림의 비중이 비교적 높은 지점이었다. 조사된 수종이 모두 활엽수종임을 감안할 때 침엽수 비중이 높은 픽셀에서는 침엽수 잎의 녹색 반사가 개엽일을 추출을 위한 식생 계절 곡선에 영향을 주었을 것으로 판단된다.

산림 식생 계절의 시작 시점은 봄철 기후의 영향이 클 것으로 판단되어 년과 다른 기후 조건을 보이는 연도를 선택하여 다른 해와 비교해 보았을 때 봄철 평균기온이 평년과 비교해 가장 낮고 강수량이 연구 기간 중 가장 높았던 2010년의 전국 평균 개엽일이 매우 늦었고(129.1일), 봄철 평균 기온이 높았던 2014년과 2016년에는 각각 112.9일, 113.8일로 연구 기간 중 가장 빠른 전국 평균 개엽일을 보였다(Fig. 7). 지역적 차이를 보면 강원 지방이 산맥을 따라 개엽이 늦게 시작되고 서해 남부 해안 지역이 비교적 빠르게 개엽이 시작되었다(Fig. 8). 강원 지방의 경우 연평균 기온 6.6-13.1°C로 서해 남부 지역 지역의 연평균(12.6-14.3°C)이 전반적으로 높고 봄철에 주로 남서 방향에서의 따뜻한 공기 유입의 영향이 있어 서해 남부 지방의 개엽이 빠르게 시작되는 것으로 보인다(Korea Meteorological Administration, 2016).

우리나라 산림의 경우에는 기후 변화에 따른 개엽일의 변화를 판단하기에는 아직 미흡하며 실측과 인공위성영상 모두 보다 장기적인 자료가 반영되어야 할 것으로 보인다. 또한 산림에서의 기상 조건의 영향을 판단하려면 산악 기상과 지형 조건을 고려하여야 할 것으로 보인다. 실측 자료의 경우 관찰자의 숙련도와 관찰 시기, 방법 등이 개엽일 값에 크게 영향을 주므로 현장 자료의 오차 역시 고려되어야 할 대상이다. 따라서 향후 MODIS뿐만이 아닌 다른 고해상도 인공위성 영상과 결합하여 보다 정밀한 공간해상도로 검증 필요가 있는 것으로 판단된다. 이 연구에서는 활엽수림 및 혼효림의 개엽일만 고려하였는데 향후 침엽수림의 개엽일 추출까지 확장하여 다양한 기후 인자들 및 지형 인자들과의 관계를 알아보고 우리나라 산림 생태계의 다양한 임상에 대한 종합적인 식물계절 연구가 진행되어야 할 것이다.