Forest Damage Detection Using Daily Normal Vegetation Index Based on Time Series LANDSAT Images

시계열 위성영상 기반 평년 식생지수 추정을 통한 산림생태계 피해 탐지 기법

Abstract

Tree growth and vitality in forest shows seasonal changes. So, in order to detect forest damage accurately, we have to use satellite images before and after damages taken at the same season. However, temporal resolution of high or medium resolution images is very low,so it is not easy to acquire satellite images of the same seasons. Therefore, in this study, we estimated spectral information of the same DOY using time-series Landsat images and used the estimates as reference values to assess forest damages. The study site is Hwasun, Jeollanam-do, where forest damage occurred due to hail and drought in 2017. Time-series vegetation index (NDVI, EVI, NDMI) maps were produced using all Landsat 8 images taken in the past 3 years. Daily normal vegetation index maps were produced through cloud removal and data interpolation processes. We analyzed the difference of daily normal vegetation index value before damage event and vegetation index value after event at the same DOY, and applied the criteria of forest damage. Finally, forest damage map based on daily normal vegetation index was produced. Forest damage map based on Landsat images could detect better subtle changes of vegetation vitality than the existing map based on UAV images. In the extreme damage areas, forest damage map based on NDMI using the SWIR band showed similar results to the existing forest damage map. The daily normal vegetation index map can used to detect forest damage more rapidly and accurately.

산림지역은 계절에 따라 생장 및 활력 특성이 변화하기 때문에 산림피해를 정확하게 탐지하기 위해서는 과거 동일한 계절적 시기의 영상정보 확보가 필요하다. 그러나 고해상도 또는 중해상도 영상은 영상촬영주가 높지 않아 동일 시기의 영상 정보들을 확보하는 것은 쉽지 않다. 따라서 본 연구에서는 산림생태계의 피해를 평가하기 위해 시계열 영상정보를 통해 피해발생 이전 과거 동일 시점의 분광정보를 추정하여 산림피해 평가의 기준정보로 활용하는 방법을 연구했다. 연구대상지는 2017년 우박과 가뭄으로 인해 대규모 산림피해가 발생한 전라남도 화순지역이며, 과거 3년간 해당 지역에서 촬영된 모든 Landsat 8 영상의 시계열 식생지수(NDVI, EVI, NDMI) 자료를 구축하고 이를 일별 연속자료로 자료보간을 실시하였다. 그리고 이를 통해 교란 발생 이전의 정상적인 일별 식생지수 추정 지도를 제작하였으며, 동일 날짜의 일별 평년 식생지수와 교란발생 이후의 식생지수의 차이값을 구하고 피해등급 기준을 적용하여 최종적인 위성자료 기반의 피해등급지도가 산출되었다. 위성기반 피해등급지도는 기존의 항공사진 기반 피해등급지도에 비해 미세한 식생활력도 변화를 효과적으로 탐지하였으며, 피해극심지역을 대상으로 비교하였을 때 SWIR 밴드를 이용한 식생지수(NDMI)가 기존의 피해등급평가 결과와 유사한 결과를 산출하여 활용도가 높은 것으로 평가되었다. 결과적으로, 일별 평년식생활력도 지도의 제작을 통해 신속하고 정확한 피해지 탐지가 가능해졌다.

I. 서론

산림은 산불, 산림병해충 등 다양한 요인으로 인해 피해를 받고 있다. 또한 최근에는 이상기상으로 인한 대규모 산림피해도 발생하고 있다. 2017년 5월 말에는 전남과 경북 일대에 대형 우박이 내려 농작물은 물론 산림지역에도 많은 피해를 입혔다. 수목의 가지가 꺾이고,잎이 떨어지며, 수간부에 상처가 발생하였는데 이와 함께 지속적인 가뭄이 겹쳐 6월 중순부터는 대규모 수목고사가 발생하였다. 전남(화순, 담양, 곡성 등)과 경북(봉화)지역을 중심으로 약 4천 ha의 산림이 피해를 입었고, 특히 전남 화순은 생태계 회복이 불가능한 피해 극심지역이 약 185 ha 발생했다(Lim et al., 2017). 우박 피해가 각 지역에서 보고된 후, 피해규모와 원인을 파악하기 위해 정부와 관련 기관이 합동조사를 실시하여 항공촬영과 현장조사를 수행하였으나, 대면적에 대한 신속한 파악에는 한계가 있었다. 대규모 피해사례는 이뿐만이 아니다. 매년 발생하는 산불 피해, 소나무재선충병으로 인한 소나무림 피해, 이상기상으로 인한 소나무림고사, 고산지역 상록침엽수 대규모 고사 등 피해발생 사례가 다양해지고 면적도 넓어지고 있는 상황이다. 이러한 산림의 피해면적을 정확히 파악하고 피해정도를 평가하여 산림복구 및 관리방안을 신속히 마련하는 것은 국가적인 해결 과제이다. 이러한 측면에서 광역적인 산림피해 현황을 주기적이고 신속하게 파악하기 위해서는 지속적으로 생산되고 있는 위성자료의 효과적인 활용 방법을 마련하는 일이 매우 중요하다.

최근 산림피해 모니터링을 위해 활용이 가능한 국내·외 위성정보의 양은 급증하고 있다. 그중전 지구적인 영상정보 제공이 무상으로 이루어지는 위성으로는 Landsat을 대표적으로 들 수 있다. 1970년대 시작된 Landsat 위성은 식생활력을 탐지할 수 있는 다양한 분광밴드를 가진 30 m 해상도 정보를 제공하고 있어 산림부문에서 오랜 기간 동안 가장 활발하게 활용되어 오고 있는 위성이며, 현재 Landsat 8호가 운영되고 있다. Landsat 위성 정보가 수십 년간 축적되면서, 장기적인 위성정보를 이용한 산림변화 추세 분석과 교란 이력 모니터링에 많은 연구가 추진되고 있다(Schroeder et al., 2014; Devries et al., 2015; Wulder et al., 2016; Cohen et al., 2017; Cohen et al., 2018). Huang et al.(2010)은 토지피복별 시계열 분광 정보 변화를 기반으로 산림변화 표준화 지수(Integrated Forest Z-score)를 개발하고, 피복별 변화패턴의 차이를 이용하여 변화 탐지를 수행하였다. Kennedy et al.(2010)는 Landsat 영상의 장기간의 변화 탐지를 통해 교란과 회복의 궤적을 변화 구간별로 파악하는 LandTrendr (Landsat-based detection of Trends in Disturbance and Recovery) 알고리즘을 제안하였다. DeVries et al.(2015)는 하모닉 형태의 모형을 활용한 정상 시계열 모형을 개발하고 범위에서 벗어나는 정도에 따라 산림 변화 형태 및 정도를 구분하기도 하였다.

장기간의 위성정보 기반 산림변화 탐지 연구는 주로 Landsat을 통해 이루어져왔다. 그러나 2015년에 발사되어 자료생산이 시작된 Sentinel 2 위성은 Landsat 보다 정밀한 해상도(10 m 이상)와 다양한 분광정보(12개 밴드)을 제공하고 있어 산림부문 활용에 대한 가능성이 커지고 있다. 최근 산림부문 활용성에 대해 Landsat과 Sentinel에 대한 다양한 비교 연구가 추진되고 있다(Korhonen et al., 2017; Sothe et al., 2017). 국내에서는 현재 고해상도 위성정보를 제공하는 아리랑위성 2, 3, 3A호가 운영되고 있으며, 2019년에는 고해상도 위성정보를 제공하는 차세대중형위성 1, 2호가 추가로 발사될 예정에 있다. 2023년에는 농림부문의 식생변화 탐지에 초점을 맞추어 개발되는 중해상도의 고주기 위성정보 생산을 목적으로 하는 차세대중형위성 발사가 계획되어 있다.

위성정보의 양은 급증하고 있지만 산림변화 분석을 위한 시계열 위성 빅데이터 활용에 대한 국내 연구는 아직 초기단계에 있다. 시계열 위성정보를 이용한 경우는 주로 MODIS를 활용한 연구에 집중되어 왔으며 풍부한 시계열 정보를 기반으로 산림변화탐지, 산림생물계절 (개엽 등) 및 생산성 변화 특성을 탐지하는 연구가 수행되어 왔다(Jung and Jang, 2013; Lee et al., 2018). 그러나 MODIS는 공간해상도가 낮아 MODIS 영상만으로는 산림 변화의 세밀한 지역적 경계를 탐지하고 지도화하는 부분에는 적용이 어려운 문제가 있다. 이를 극복하기 위해 MODIS의 시계열적 정보와 Landsat 영상의 공간해상도 수준을 결합하기 위한 영상융합(FSDAF) 기법을 연구한 사례도 있다. Jin et al.(2017)은 MODIS를 기반으로 한 토지피복별 시계열적 변화를 Landsat 영상에 적용하여 미관측 시기의 분광정보를 예측하는 연구를 수행한 바도 있다.

본 연구에서는 Landsat의 시계열 위성정보를 이용하여 산림피해지의 경계를 보다 정확하게 추출하고 피해정도를 평가하기 위한 기술 개발을 목적으로 한다. 일반적으로 기존 피해지 탐지는 두 가지 방향으로 진행되어왔다. 첫째, 피해 산림과 피해를 입지 않은 인근 산림의 분광특성 차이를 이용하여 피해산림을 분류한다. 둘째, 피해 전과 후의 상태를 비교하여 그 차이를 이용하여 피해산림을 분류하고 피해정도를 구분한다. 첫 번째 접근은 피해발생 이후의 단시기 영상만 확보된 경우에 주로 적용되며, 주변산림과의 비교를 통해 평가가 되므로 피해지 본래의 식생특성에 대한 고려와 피해정도(피해등급) 구분이 어렵다. 두 번째 접근은 피해발생 전후를 비교한다는 측면에서 피해지역의 실질적인 변화를 파악할 수 있다는 장점이 있으나, 이는 피해발생에 따른 변화를 제외한 다른 변화요인은 없어야 한다는 것을 전제로 한다. 그러나 산림은 계절에 따라 변화하고 시기에 따라 생장하는 변화하는 생명체이기 때문에 이러한 내재되어 있는 변화요인을 제거하지 못하면 정확한 변화탐지가 불가능하다. 즉, 유사한 계절적 영상의 확보를 하지 못하면 정확한 분석에 어려움이 발생한다. 따라서 두 가지 경우의 한계를 보완하기 위해 본 연구에서는 피해발생 시기와 동일한 시기의 피해발생 이전 정보를 추정하여 이의 차이를 분석하는 것으로 연구의 방향을 설정하였다.

따라서, 본 연구의 목적은 산림피해 정도를 신속하고 정확하게 파악하기 위해 위성정보 기반 과거 건강한 산림활력도의 시계열 정보 산출 및 활용 방법을 개발하는 것이다. 이를 위해 확인되어야 하는 질문은 다음과 같다. 첫째, 촬영주기가 세밀하지 못한(약 2주) Landsat 위성정보를 이용하여 시계열적으로 연속된 과거 시점의 건강한 상태의 활력도 정보를 생산할 수 있을까, 둘째, 산림피해지를 파악의 정확도를 높이기 위해 어떤 영상 밴드와 지수를 활용할 수 있을까, 셋째, 과거기준 평년 산림활력도 정보를 이용해 산림의 변화를 얼마나 정량적으로 평가할 수 있을까.

이를 대한 해답을 마련하기 위해 본 연구에서는 산림피해 탐지의 기준 마련을 위한 평년 식생활력도라는 개념을 도입하였고, 위성 빅데이터를 이용하여 이를 추정하기 위한 기법을 개발하였다. 또한 평년 식생활력도 기준으로 피해지 탐지와 평가를 수행하고 활용성을 평가하였다.

2. 연구자료 및 방법

1) 연구대상지 및 연구자료

본 연구에서 산림생태계 피해 평가 대상지로 2017년 대규모 우박피해가 발생한 전라남도 화순군 동복면 피해지를 포함한 10 km×9 km 영역을 선정하였다(Fig. 1). 해당 대상지에 대해서는 드론 영상이 촬영되었으며 현장조사와 드론영상을 기반으로 피해지 도면이 기 제작되어 있다(Fig. 2).

Fig. 1. Satellite images in study area damaged by hail and drought.

Fig. 2. Forest damage class map based on airphoto images (Lim et al., 2017).

Fig. 2의 피해현황도는 항공사진을 육안 판독하여 피해지 영역을 산출하고 RGB 화소 분석과 현지조사 대조를 통해 피해 등급을 구분한 결과로 제작되어(Lim et al., 2017), 피해등급 평가, 특히 [심], [중] 지역 구분의 신뢰도가 높다고 할 수 있다. 그러나 드론 항공영상 촬영 대상지가 전체 피해지 중 피해가 약했던 일부 지역(남동부)를 포괄하지 못한 문제가 있었고 피해등급 [경]의 경우는 육안판독으로 피해지 외곽 경계를 명확히 구분하기 어려운 지역이 다수 있었던 점을 고려하면 전체 피해발생지 영역 구분은 정확도가 높지 못하다고 말할 수 있다.

따라서 본 연구에서는 위성영상을 이용하여 전체 피해지를 산출하고 새로운 방법으로 피해등급지도를 산출하고, 항공사진과 현지조사 기반의 기존 피해등급지도의 피해 [심] 지역을 기준으로 정확도 평가를 수행하였다. 그리고 피해 [중], [경] 지역 산출의 평가에 대한 부분은 항공사진과 위성영상 활용의 효과를 상호 비교하였다.

위성영상 자료를 이용하여 산림생태계 피해를 평가하기 위해 장기간에 걸쳐 산림의 활력도 변화를 파악할 수 있는 Landsat 자료를 이용하였다. 피해 이전 해당 지역의 건강한 산림상태를 파악하기 위해 과거 위성영상 자료를 수집하였다. 즉, 산림생태계 피해를 평가하기 위한 기준자료(평년 식생활력도) 생성을 위해 2014년에서 2016년 사이에 촬영된 Landsat Level1 영상(Path/Row: 115/36) 64장을 수집하였다. 또한 피해 이후의 변화를 파악하기 위해 우박피해 이후에 촬영된 영상(2017.6.16. 촬영)을 이용하였다.

2) 연구 방법

본 연구에서는 위성영상 기반 피해현황도 제작을 위한 최적의 방법을 찾기 위해, 먼저 위성영상 전처리를 수행하고 피해 여부를 평가하기 위한 다양한 식생지수영상을 제작하였다. 또한 과거 시계열 영상을 통해 피해발생 이전의 평년 식생활력도를 추정하기 위해 최적의 자료보간 방법을 파악하였다. 최종적으로는 평년 식생활력도와 피해 이후의 영상정보의 차이를 통해 피해 지역을 산출하였다. 평년활력도 기반 피해현황 분석 흐름도는 Fig. 3과 같다.

Fig. 3. Research flow.

(1) 시계열 위성영상 자료 전처리와 식생지수 변환

위성영상을 이용한 시계열 식생지수 정보를 이용하기 위해서는 방사보정, 대기보정, 정사보정이 된 영상을 활용해야 한다. 이를 통해 위치가 명확한 지표면 반사도 정보의 시계열적 변화에 대한 비교가 가능하다. 본 연구에서는 USGS에서 제공하고 있는 대기보정과 정사 보정이 수행된 Landsat 8 Level2 영상을 이용하였다. 대기보정은 LaSRC(Landsat 8 Surface Reflectance Code)에 의해 이루어졌으며, 이 방법은 MODIS 영상으로부터 에어로졸, 수증기, 오존 등 대기 중 산란에 영향을 미치는 요인들을 추정하여 Landsat8 다중분광 영상의 대기보정을 수행하는 방법이다(USGS, 2019).

대기보정된 영상을 다운로드 받은 후에 식생활력도 비교를 위한 식생지수로 변환하였다. 식생지수로는 NDVI(Normalized Difference Vegetation Index), EVI (EnhancedVegetation Index),NDMI(NormalizedDifference Moisture Index)를 이용하였다(USGS, 2017).

\(N D V I=\frac{N I R-R e d}{N I R+R e d}\)       (1)

\(E V I=2.5 *\left(\frac{N I R-R e d}{N I R+6 * R e d-7.5 * B l u e+1}\right)\)       (2)

\(N D M I=\frac{N I R-S W I R 1}{N I R+S W I R 1}\)       (3)

기존에 산림분석에서 산림피복분류를 위해 다양하게 이용되는 영상정보(아리랑 위성영상, RapidEye 등)가 대기보정 없이 이용되고 있는 것을 고려했을 때, 대기보정을 수행하지 않은 영상을 통해서도 유사한 결과를 산출할 수 있을지에 대한 평가를 함께 수행했다. 따라서 대기보정이 수행된 Level2 영상과 대기보정이 수행되지 않고 방사보정만 수행한 Level1 영상의 반사도 정보 및 식생지수의 차이를 비교하고 시계열 분석을 위한 대기보정 효과에 대한 평가를 수행했다. 평가는 연구대상지 내 약 1 km * 1 km 산림지역에 대해 QA Band 상에서 Clear Sky로 평가되어 대기간섭 효과가 적다고 평가된 영상자료를 추출하여 비교·평가하였다.

(2) 일별 평년 식생지수 정보의 추정

교란이 발생한 시기의 피해정도를 평가하기 위해서는 교란이 발생하지 않은 이전시기의 자료와 비교하는 것이 가장 효과적이다. 특히 산림의 경우 시계열적인 변화가 있기 때문에, 교란발생시기와 유사한 시기의 정보를 이용하는 것이 중요하다. 그러나 대부분 동일한 시기의 자료를 취득하기 어려운 것이 현실이다. 따라서 해당 시기의 앞시기와 뒷시기 정보를 이용하여 해당 시기의 정보를 추정하는 방법의 활용이 필요하다. 이를 위해 모형 추정과 자료보간 방법 적용이 가능하다.

모형추정을 위해서는 식생의 시계열적 변화와 유사한 모형의 형태를 가지는 Gaussian 함수 또는 Double logistic 함수를 적용할 수 있다. 자료의 개수가 많지 않을 경우에도 적용이 비교적 용이한 장점이 있으나, 실측 자료의 변화 추세가 함수 형태에 적합하지 않을 경우에는 많은 오류가 발생할 수 있는 단점이 있다. 자료 보간(Interpolation)은 알려진 데이터를 이용하여 해당 범위 내에서 비어있는 지점의 새로운 데이터를 채우는 것을 말한다. Spline 등의 자료 보간 기법과 함께 실측데이터의 세부적인 변이를 제거하고 추세를 파악하기 위한 평활 방법(Saviztky-Golay filter 등)이 함께 적용될 수 있다. 이 방법은 충분한 연속된 시계열 자료가 확보되어야 변화추세를 파악할 수 있다는 한계가 있다.

• Gaussian function

가우스 함수의 그래프는 좌우 대칭의 종 모양의 곡선으로 +/– 극한을 향하면서 급격히 감소하는 특성을 가진다. 여기서 a와 c는 곡선의 형태를 결정, b는 그래프 중심의 x좌표 이동, d는 y절편에 영향을 주는 모형의 파라미터이다. 식생활력도가 연중 증가했다가 감소하는 형태를 보이기 때문에 가우스 함수의 이러한 특징 적용이 가능하다.

\(y=a \exp \left(\frac{-(x-b)^{2}}{2 c^{2}}\right)+a\)       (4)

● Double logistic function

식생의 생장형태와 가장 유사한 곡선은 상한이 있는 로지스틱 곡선(S자 곡선)으로 알려져 있다(Zhang et al., 2003). 식생활력도의 연중변화 역시 이러한 변화와 가장 유사한 형태를 보인다. 식생활력도가 봄철에 급격히 증가하고 여름에는 증가폭이 감소하고 일정수준 유지되다가 감소하는 형태를 보인다. 겨울~여름까지는 증가하는 로지스틱 함수, 여름부터 다시 겨울까지는 감소하는 Double Logistic 함수(Roper, 2000; Beck et al., 2006) 를 적용할 수 있다. 여기서 a는 시계열 변화의 y축의 폭, b는 그래프 중심의 x좌표 이동, c는 최대치 유지의 폭, d는 증가사면의 기울기, e는 감소사면의 기울기, f는 y절편에 영향을 주는 모형의 파라미터이다.

\(\begin{aligned} y=& \frac{a}{1+\exp \left(-\frac{x-b+c / 2}{d}\right)} \times \\ &\left[1-\frac{1}{1+\exp \left(-\frac{x-b+c / 2}{c}\right)}\right]+f \end{aligned}\)       (5)

● Spline + Saviztky-Golay filter

자료 보간의 대표적인 방법으로 관측지점의 정보를 이어주는 Spline 방법이 있다. 1차 스플라인은 해당 지점의 점들을 1차방정식으로 이어주는 방식으로 가장 간단하고 직관적인 보간 방법이다. 선형 Spline에서 관측점을 경계로 발생하는 단절을 완화하기 위해 2차와 3차 방정식을 적용하여 부드럽게 자료를 생산해 낼 수 있다. 그러나 Spline의 원칙은 모든 점을 다 통과해야 하므로 자료에 튀는 값이 존재하거나 결측치 간격이 넓은 경우에 2차 이상의 Spline 적용에는 오류가 발생할 수 있다.

데이터 평활을 위한 일단위 자료 제작을 위해 오류발생이 적은 1차 스플라인을 적용하였고, 이 자료에 Saviztky-Golay filter(Savitzky and Golay, 1964)를 적용하여 자료 평활을 수행했다. Saviztky-Golay filter는 Moving average의 한 종류로, 데이터 경향성에 가까운 데이터는 가중치를 더 주고 멀리 떨어진 데이터는 가중치를 적게 주는 변형함수이다. 선형 보간 자료의 Saviztky-Golay filter를 적용하되, 윈도우의 크기를 31(1개월), 61(2개월), 91(3개월)로 변경하여 산출한 결과는 다음과 같다. 윈도우 크기를 증가할수록 부드러운 그래프가 나오나 실측 자료와는 차이가 났으나, 큰 차이가 발생하지는 않았다. 본 연구에서는 window의 크기를 61로 설정하여 자료 평활을 수행하였다.

(3) 일별 평년 식생지수 지도 제작과 피해지역 평가

각 픽셀별 시계열 정보를 2.2.2에서 선정한 방법을 통해 일별 자료로 보간하여 일별 식생지수 지도를 제작하였다. 또한 피해지역 평가를 위해 피해 전과 후의 자료를 준비했다. 피해 후의 자료는 2017년 6월 16일(DOY: 167) 영상이었고, 피해 전의 자료는 일별 자료를 보간하여 제작한 167일의 평년 식생지수 지도 영상이었다. 두 자료의 차이값을 통해 피해정도를 파악할 수 있었다. 피해정도는 값의 정도에 따라 등급으로 구분될 수 있으며, 등급의 임계치 기준은 실제 피해와 연계하여 구분되어야 한다. 기존의 피해 등급은 실측자료에 기반하여 잎의 손실률에 따라 심, 중경등 3단계로 나눈바 있다 (Table 1). 이를 기반으로 해당픽셀의 식생지수의 변화율이 50% 이상인 지역을 [심], 25~50% 이상인 지역을 [중], 10~25% 이하인 지역을 [경]으로 구분하고 피해 기존의 피해등급 산정 결과와 비교하였다.

Table 1. Definition of forest damage class

잎의 손실율, 즉 식생지수의 감소율은 식2에 따라 산정되었다. 식생지수의 변화폭 산정을 위한 식생지수 최소값은 해당 지역 활엽수림의 잎이 전혀 없는 겨울시기 식생지수의 값을 적용하였다.

\(Reduction\ Ratio\ of\ V I=\frac{V I_{\text {before }}-V I_{\text {after }}}{V I_{\text {before }}-V I_{\text {min }}}\)       (6)

3. 연구 결과

1) 산림교란과 식생지수의 변화

건강한 산림의 일반적인 식생지수 연중 변화 경향을 파악하기 위해 최근 3년간의 자료를 이용하여 식생지수 시계열 변화를 파악하였다. 그 결과, 건강한 활엽수림은 봄철에 식생지수가 급격히 증가했다가 여름과 가을에 걸쳐서 유지되고 늦가을부터 지수가 급격히 감소하는 경향을 보였다. 잎을 연중 달고 있는 침엽수림은 활엽수림보다는 식생지수 변동이 있지는 않지만 겨울철에 약간의 감소추세가 있고 당년 잎이 나는 봄철시기에 약간의 증가추세를 보였다(Fig. 4).

Fig. 4. Change of vegetation index before or after forest damage (a) Deciduous broadleaf forest (b) Evergreen conifer Forest.

그러나 산림피해가 발생할 경우 이러한 패턴은 사라지며 식생지수가 급격히 낮아지는 경향을 보였다. Fig. 4(a)에서는 최근 3년 동안 건강한 활엽수림이었다가 심각한 우박피해를 받은 지점의 사례를 보여주고 있다. 최근 3년간의 일반적인 경향과 비교했을 때 상당히 큰 식생지수 감소를 보이고 있다. 활엽수의 경우 겨울철에 잎이 없기 때문에 잎에 의한 활력도를 0이라 볼 수 있는데,겨울철의 활력도와 비교하면 잎의 손실량을 유추할 수 있다. Fig. 4(b)는 최근 3년 동안 건강한 침엽수림이었다가 심각한 우박피해를 받은 지점의 사례이다. 이 지점 역시 기존의 식생지수와 비교했을 때 상당한 감소를 보였다. 그러나, 해당 지역은 겨울에도 활력도가 높은 침엽수림 지역이기 때문에 해당 지역 정보만으로는 잎이 모두 손실된 최저 활력도 값을 추출할 수 없다. 따라서 주변 활엽수림의 최저 활력도(겨울철 활력도) 자료를 이용하여 잎의 손실량 파악이 가능하다.

2) 시계열 위성정보 분석을 위한 대기보정 영향 평가

3.1의 분석 시 영상의 DN값에서 방사보정과 대기보정을 거친 표면반사도(Surface Reflectance) 값을 이용하였다. 지표반사도의 절대적인 값을 분석의 기준정보로 활용하거나 다양한 종류와 시기의 영상을 복합적으로 활용할 경우에는 대기보정을 통해 실제의 지표반사도자료를 산출해야 한다. 그러나 영상분류를 수행하거나 임계치 절대값을 적용하지 않아도 될 경우, 또는 대기보정에 영향을 받지 않는 분광정보의 경우에는 방사보정만을 거친 대기 상층의 반사도(TOA Reflectance)값을 그대로 이용하는 방법도 적용이 가능하다. 대기보정은 추가적인 영상 처리가 필요한 과정이기 때문에 본 연구의 프로세스에서 대기보정 과정이 꼭 필요한 작업인지, 대기보정을 수행하지 않을 경우 어떤 영향을 받을지에 대한 사전 검토가 필요했다.

이를 위해 대기보정 평가 대상지(전체 연구대상지 내 1 km * 1 km) 내 픽셀들의 TOA Reflectance와 Surface Reflectance를 밴드별로 비교하고 이에 따라 산출된 식생지수도 각각 비교하였다. 그 결과, 가시광선 밴드(Blue, Green, Red)는 대기보정 이후 반사도 값에 크게 변동이 있는 것으로 나타났다. 대기보정 후 가시광선 밴드의 반사도값이 전체적으로 감소되는 경향을 보였고, 파장대가 작을수록 차이는 더 큰 것으로 나타났다. 그러나 적외선 밴드(NIR, SWIR1, SWIR2)는 대기보정 여부에 큰 영향을 받지 않는 것으로 나타났다.

이에 따라, Red 밴드를 이용하는 NDVI의 경우 보정 이후의 지수값이 보정 전의 값보다 크게 산출되었다. Red 밴드와 Blue 밴드를 이용하는 EVI의 경우 보정 이후의 지수값이 보정 전의 값보다 상당히 낮게 산출되었으며, 지수값이 커질수록 그 영향이 커지는 것으로 나타났다. 반면, 적외밴드만을 이용하는 NDMI는 보정 전과 후의 지수값이 크게 변동이 없는 것으로 확인되었다.

Fig. 6. Comparison of atmospheric correction effects (a) NDVI (b) EVI (c) NDMI.

이는 일반적으로 파장이 작을수록 에어로졸에 의한 영향이 커지므로 에어로졸 효과를 제거하는 대기보정에 의해 가시광 대역의 영향이 크게 나타났기 때문으로 볼 수 있다. 반면 수증기량에 영향을 많이 받는 적외선 밴드의 경우 흐린 날이나 구름에 의해서 영향을 많이 받지만, 본 연구에서는 청명한 지역의 픽셀만을 추출해서 분석을 수행하였으므로 이에 대한 효과는 크지 않았던 것으로 판단된다. 즉, 맑은 하늘이 있는 지역(clear sky)만을 추출해서 연구할 경우, 가시광선을 포함한 식생지수 절대값 이용에는 대기보정이 필요하며 적외선 밴드만을 이용하는 식생지수는 대기보정 여부에 상대적으로 적은 의존도를 보였다.

따라서, NDVI와 EVI의 절대값을 이용해야 할 경우에는 대기보정을 수행해야 하며 NDMI를 적용할 경우에는 대기보정을 수행하지 않아도 무방할 것으로 판단된다.

3) 시계열 정보 보간 및 평활 기법 평가

일별 평년 식생지수 추정 기법을 선정하기 위해 곡선추정접근으로Gaussian Function,Double Logistic Function 방법을 적용, 자료보간 접근으로 Spline + Saviztky-Golay filter 방법(1차 Spline으로 일별 자료 보간, SV Filter 이동 평균)을 적용하여 세 가지 식생지수에 대해 그 결과를 비교하였다(Table 2, 3, 4).

Table 2. Comparison of time-series data interpolation techniques in deciduous forest

Table 3. Comparison of time-series data interpolation techniques in coniferous forest

Table 4. Comparison of damaged area in various damage class maps

활엽수림의 경우 식생지수가 연중 증가했다가 감소하는 특성이 뚜렷한데, Double Logistic Function이 식생지수의 연중 변화 경향을 효과적으로 추정하였다. Spline + Saviztky-Golay filter의 방법 역시 식생지수의 연중 변화를 잘 추정하였으나, 관측자료가 부족한 겨울철 기간에는 비현실적인 추정치를 보였다. 침엽수림의 경우 식생지수의 연중 증감 특성이 활엽수보다는 뚜렷하지는 않았는데, 이 경우 Spline + Saviztky-Golay filter가 관측값의 추세를 잘 반영하는 추정치를 보였다. 값의 변동이 크지 않은 NDMI의 경우 Double Logistic Function의 함수 추정에 의해 일부 구간에서 비현실적인 추정 결과가 발견되었다. 이에 따라, 원자료의 양이 자료 보간에 충분한 봄~가을철에 적용할 경우 원자료의 특징을 가장 잘 반영하는 Spline + Saviztky-Golay filter 방법을 가장 효과적인 일별 평년 식생지수 추정 기법으로 선정하였다.

자료 보간에 따른 픽셀별 불확도를 평가하기 위해 3월~11월의 관측값과 예측값의 오차(RMSE)를 산정하였다. 변화 경향에서 벗어나는 값이 많아질수록 RMSE이 거치고 식생지수 추정치의 신뢰도가 낮아진다고 볼 수 있다. 불확도 평가 시 적설에 의해 오차발생 가능성이 높은 겨울철 정보는 제외하였다.

4) 평년 활력도 기반 피해 현황 분석

각 픽셀별 시계열 식생지수 정보는 Spline + SaviztkyGolay filter 방법을 통해 일별 자료보간 과정을 거쳤고, 이에 따라 각 365일에 해당하는 3년 평년 식생지수 지도가 제작되었다. 그리고 피해 발생 이후 시기(2017년 6월 16일) 영상과의 비교를 위해 평년 DOY 167의 식생지수 영상(Fig. 7(a)~(c))을 피해 이전의 기준자료로 이용하였다. 평년 식생지수 영상에서 적외 밴드를 이용하여 수분상태를 파악하는 지수인 NDMI는 식생지역 내에 전체적으로 유사한 값을 보였다. 이는 임상 특성에 따라 수분 보유 정도가 크게 차이가 나지 않는다는 것을 보여준다. NDVI는 NDMI보다는 임상 특성에 따른 산림 내 편차가 약간 있지만 EVI 보다는 일관적인 특성을 보였다. NDVI와 달리 Blue band를 추가로 활용하는 EVI는 임상 특성과 지형 조건에 따른 변동성이 상대적으로 컸다.

Fig. 7. Change analysis using normal vegetation index estimates.

식생피해정도를 파악하기 위해 평년 식생지수 추정자료에서 피해발생 이후 실제 식생지수(Fig. 8(d)~(f)) 값의 차이를 구하였고, 이 중에서 산림지역 내 식생지수가 감소된 지역만 그림으로 표출하였다(Fig. 8(g)~(i)). NDVI와 NDMI의 경우 피해가 없었던 지역은 0에 가까운 값을 가지며 피해발생지역은 NDVI의 차이가 최대 약 0.47, NDMI의 차이는 약 0.55까지 발생했다. EVI의 경우에는 피해가 발생하지 않은 지역에서도 식생지수의 감소현상이 산발적으로 나타났으며, 이에 따라 피해 지역 및 등급의 정확한 추출에 있어서 EVI는 유용한 식생지수로 활용되기 어려울 것으로 판단된다.

Fig. 8. Forest damage class mapping using change analysis of vegetation index.

각 픽셀별 평년 식생지수를 기반으로 한 변화율 산정 결과, NDVI와 NDMI는 피해지역의 영역이 명확히 확정되었고 그 안에서 피해등급별로 구분이 되었다. 그러나 EVI의 경우에는 피해발생지역이 산발적으로 흩어져 있는 것으로 나타나, 우박이동 경로에 따라 집중적으로 피해가 발생한 실제 피해 특성과는 부합하지 않은 결과를 보였다(Fig. 8(a)~(c)).

본 피해등급은 평년 식생지수의 추정을 기반으로 제작되었기 때문에, 추정치의 불확도가 높으면 피해등급정보 역시 불확도가 높다고 볼 수 있다. 3년 동안의 영상을 함께 이용하여 산림의 연간 시계열 평년 자료를 만드는 과정은 해당 지역이 큰 피복의 변화 및 교란을 겪지 않았다는 점을 기본 가정으로 전제하고 있다. 따라서 숲가꾸기, 벌목, 산지개발 등 산림피복의 변화가 있는 지역은 추정치의 오차가 크게 발생할 수밖에 없다. 따라서 피해등급 평가결과를 활용함에 있어서 평년 자료 추정의 불확도가 높은 지역은 제외하는 것이 필요하다(Fig. 8(d)~(f)).

Fig. 8의 피해등급별의 평가정확도를 구하기 위해서는 기존의 무인항공사진 기반 피해등급도(Fig. 2)와의 비교가 필요하다. 그러나 기존 피해등급도 제작을 위한 무인항공사진 촬영 범위가 전체 피해 대상지를 포괄하지 못하는 문제가 있었다. 무인항공사진의 경우 대면적 촬영이 어렵기 때문에 대략적인 현장 확인을 통해 촬영대상지를 미리 확정하게 되는데, 이에 따라 화순 피해지의 우측 하단 부분의 무인기 촬영이 수행되지 못했다 (NIFoS, 2017). 누락지역의 경우 피해 “중”과 “경”지역이 매우 넓게 분포하고 있는 지역이어서 전체 피해대상지 산정에 큰 차이가 발생하게 되는 원인이 되었다. 이러한 측면에서 항공사진에서 수행하지 못하는 대면적탐지가 가능한 부분은 위성영상을 이용한 분석의 주요한 장점으로 볼 수 있다.

위성자료 기반의 피해등급도의 피해등급 산정결과와 기존의 무인항공사진 기반의 피해등급도를 정량적으로 비교·평가하기 위해, 무인항공사진이 촬영된 범위에 한정하여 두 자료를 비교하였다. 그 결과, 피해극심지는 NDMI를 이용한 피해등급 면적이 기존의 평가결과와 가장 유사하였으며 NDVI는 과소추정, EVI는 과대추정을 하는 것으로 나타났다. 피해지 전체 면적은 세 지수의 경우 모두 기존 피해등급지도상 면적보다 2배 이상 크게 산출되었다. 이는 기존의 피해등급지도 제작 시 활용한 RGB 항공사진이 피해 극심지 추출을 용이하지만 피해 “중”, “경”지역의 구분 즉, 약도의 피해여부 탐지에는 한계가 있었기 때문으로 판단된다. 즉, 역으로 보면 위성영상의 활용 시 미세한 산림활력 변화 여부를 항공사진보다 더 효과적으로 탐지할 수 있다는 것을 의미한다.

4. 결론 및 고찰

원격탐사 자료를 이용하여 대규모 산림피해 현황을 파악하기 위해서는 피해 이전의 영상과 피해 이후의 영상의 변화를 파악해야 정확한 변화 특성을 파악할 수 있다. 산림은 계절에 따라서 기본적인 연간 변화를 겪기 때문에, 피해에 따른 변화를 정확히 파악하기 위해서는 계절적 변화에 의한 영향을 제거해야 하며 이를 위해 두 시계열 영상의 비교 대상 시기가 동일해야 한다는 기본 전제가 만족되어야 한다. 그러나 실제 산림교란이 발생했을 때 우리가 필요로 하는 동일한 시기의 과거와 현재 자료를 취득하는 것이 쉬운 일이 아니다.

이러한 문제를 극복하기 위해 본 연구에서는 1년 중 동일시기를 기준으로 피해 이전과 피해 이후를 비교할 수 있는 방법을 개발하였다. 이를 위해 해당 지역의 평년 식생활력도라는 개념을 제시하였다. 이는 기상자료에서 흔히 이야기 하는 평년기후와 유사한 개념으로서, 30년 평년기후라고 하면 30년간의 기상정보를 이용하여 평균적인 기후특성을 도출하는 것을 말한다. 기상 이벤트 또는 기상특성 변화를 평가할 때 30년 평년대비 어느 정도로 변화하였는지 정량적으로 평가한다.

이와 유사하게 평년 식생활력도의 개념은 큰 피해가 발생하기 이전의 해당 지역의 평균적인 식생활력도의 연간 변화 정보를 말한다. 과거의 건강한 상태 정보를 미리 알게 되면, 이후 큰 피해가 발생했을 때 손쉽게 변화의 폭을 파악할 수 있는 장점이 있다. 그러나 식생의 경우 생장을 통해 지속적으로 변화하기 때문에 기상정보와는 달리 30년 등 장기간의 평년 정보를 도출하는 것은 유의성이 떨어지는 차이가 있다. 따라서 시계열 정보를 파악하기에 충분한 자료가 모아져야 하는 동시에 자료의 전체 기간은 최대한 짧은 것이 바람직하다고 볼 수 있다.

본 연구에서는 Landsat 8을 이용하여 피해 발생 이전 시기의 평년 식생활력도를 추정하기 위해서는 약 3년 정도의 자료를 모으면 시계열적인 추정이 가능하다는 것을 확인하였다. Landsat의 재방문주기는 16일이어서 1년에 20장 이상의 영상이 생산되지만, 흐리거나 구름이 있어 이용할 수 없는 정보가 많이 포함되어 있기 때문에 1년의 촬영 자료만으로는 연간 시계열 정보 확보가 어렵기 때문이다. 영상 정보가 충분하지 않기 때문에 영상 일부에 구름이 있다고 해서 영상 전체를 이용하지 않게 되면 일부 포함되어 있는 깨끗한 정보도 이용이 어렵기 때문에 픽셀별로 이용 가능한 정보를 평가하여 가능하면 최대한 많은 정보를 추출해서 이용하는 것이 필요하다. 이를 위해서는 영상 내의 구름 여부와 영상 질 평가를 포함한 QA 밴드의 정확도가 매우 중요하다.

피해지 탐지를 위해 또한 확인되어야 하는 것은 어떤 지수가 효과적으로 피해지 탐지하는데 유용하게 활용될 수 있느냐 하는 부분이다. 본 연구에서는 NDVI, EVI, NDMI를 적용하고 비교하였다. EVI는 NDVI를 개선하기 위해 개발된 지수로서, NDVI 값이 큰 지역에서 활력도의 변별력이 없어지는(값이 포화되는) 문제를 해결하기 위해 Blue 밴드를 추가적으로 이용하였다. 그러나 본 연구의 피해지 탐지에 있어서는 EVI가 효과적인 변화탐지를 수행하지 못하는 것으로 평가되었다. 즉, 피해지가 아닌 지역에도 변화가 크게 나타나는 경향을 보였고 이는 식생의 변화뿐 아니라 다른 변화에 의한 오염이 발생한 것으로 판단된다. 따라서 식생변화 탐지에서는 식생상태 파악에 효과적이라고 평가된 Red, NIR, SWIR 밴드를 이용한 NDVI와 NDMI가 활용이 가능하다고 볼 수 있다.

평가등급 구분에 있어서는 NDMI가 NDVI 보다 기존의 피해등급(피해 심) 기준에 상대적으로 더 유사한 결과를 도출하였으며, 이는 현장에서 파악하는 잎의 손실율과 NDMI의 값의 연계가 정확도가 높다는 것을 말해준다. 특히 NDMI의 경우 맑은 날 영상에서 대기보정에 영향을 거의 받지 않기 때문에 대기보정 수행 영상이 부재할 경우에도 적용이 가능한 큰 장점이 있다. Cohen et al.(2018)은 산림교란 탐지에 있어서 밴드와 지수의 평가 정확도를 비교하였으며, 그 결과 SWIR 밴드와 SWIR를 이용한 지수가 가장 정확도가 높다고 평가한 바도 있다. 추후 국내 교란탐지에 부합하는 SWIR 밴드의 효과적 활용 방안에 대한 연구가 필요하다.

본 연구에서는 영상 전처리와 자료 추출, 식생지수 변환, 미측정 시기의 자료 보간의 과정을 거쳐 일별 평년 식생활력 지도를 제작할 수 있었다. 본 논문에서는 피해 발생 전후 비교대상이 되는 DOY 167일의 평년 식생활력지도만 제시되었지만, 분석 과정을 통해 실제 365장의 일별 평년 식생활력도 지도가 제작이 되었다. 이는 어떤 시기에 이벤트가 발생했다 할지라도 곧바로 바로 변화 비교가 가능해졌다는 것을 의미한다. 이러한 접근방법은 본 연구의 우박피해 사례뿐 아니라, 산불피해 및 산지훼손 등 다양한 변화탐지에 활용이 가능한 효과적인 변화탐지 방법이라 볼 수 있다.

그러나 본 연구의 방법을 전국적으로 적용하기 위해서는 몇 가지 고려해야 할 사항이 있다. 첫째, 시계열 정보 추정을 위해 위성정보의 축적이 충분이 이루어져야 한다. 즉, 식생의 연변화 추정이 가능할 만한 양의 영상이 생산되고 경우는 현재로서는 Landsat 8 밖에는 없는 상황이다. 따라서 Landsat 8 공간해상도(30 m)보다 더 정밀한 변화탐지가 필요한 사례는 적용이 어려운 문제가 있다. 즉, 작은 규모의 산림 피해나 산지훼손 또는 선형으로 발생하는 산사태 피해 등은 정확한 탐지가 어렵다. Sentinel 2는 촬영주기가 1주일 이내이고 해상도가 10 m로 본 연구의 프레임을 적용할 수 있으나, 아직 영상의 축적량이 많지 않아 활용에 제한이 있다. 그러나 향후 영상정보가 지속적으로 축적되면 활용도가 점차 확대될 것으로 판단된다. 또한 향후 2022년에 발사 예정으로 추진되고 있는 농림업중형위성은 한반도를 약 3일 주기로 촬영하는 공간해상도 5 m 영상 생산을 목표로 하고 있어 향후 적용 대상이 점차 확대될 수 있다.

둘째, 산림 피해는 전국 어느 지역에서 발생할지 알 수 없기 때문에, 피해 발생 후 즉각적인 현황 파악을 위해서는 평년 식생활력 상태에 대한 지도가 전국적으로 사전에 구축되어 있어야 한다. 이를 위해서는 전국적인 위성정보를 지속적으로 구축하고 평년 식생활력 상황에 대해 1년 단위로 지도를 생산해 놓는 것이 필요하다. 예를 들면, 2019년을 위해 2016~2018년 3년 자료를 이용한 식생 평년자료를 생산하고, 다음 해가 되면 2020년을 위해 2017~2019년 3년 자료를 이용해 식생 평년자료를 갱신하는 것이다. 이렇게 1년 단위로 갱신된 평년 식생활력 자료는 산림변화분석의 기준 자료가 될 뿐만 아니라 장기적인 산림변화 추세를 이해하는 데에도 도움을 줄 수 있다. 이를 위해, 영상의 분석기술 개발을 넘어 위성영상의 수집과 체계적 관리, 지속적인 자료 갱신 체계에 대한 고민이 병행되어야 한다.

마지막으로, 평년 식생활력도 개념은 근본적인 적용의 한계가 있다. 식생 구조가 안정화된 산림은 평년의 산림활력 정보의 도출이 가능하고 이를 기준으로 교란에 대한 평가가 가능하지만, 한번 교란이 발생하게 되면 해당 지역에 대해서는 평년에 대한 개념의 적용이 가능하지 않다. 또한 교란 이후 산림은 회복의 과정을 거치기 때문에 지속적으로 변화하게 된다. 즉, 최근 시기에 교란이 발생하여 산림이 변화 과정에 있는 지역의 경우에는 수년간 정보를 이용하여 평년 식생활력도 연간 변화를 통계적으로 유의한 수준으로 산출할 수 없는 한계가 있다. 이러한 지역은 교란 이력을 별도 관리를 하고 차별화된 평가 기법을 개발하는 것이 필요하다. 추후, 이런 부분을 보완하기 위해 장기간의 위성정보를 이용한 산림교란 발생과 변화의 추적, 교란 이력의 관리, 산림의 회복 과정 모니터링 방법에 대해 연구가 필요하다.