Landsat 8 OLI영상의 NDVI를 이용한 식생피복지수 분석

Analysis of Vegetation Cover Fraction on Landsat OLI using NDVI

Abstract

대기의 에너지를 측정하거나 지표면유출을 예측하는 기상 및 수문모델에서 지표면특성(식생피복)을 파악하는 것은 매우 중요한 요소이다. 1978년 Deardorff가 식생피복을 정량적으로 파악하기 위하여 식생피복지수(Vegetation Cover Fraction)를 제안한 후 식생피복지수에 관한 연구가 활발해졌다. 그러나 선행연구에서는 AVHRR, MODIS 그리고 KOMPSAT-2영상과 같은 고 저해상도 위성영상을 이용한 많은 연구가 있었으나, 중해상도 영상인 Landsat에 대한 연구는 미비한 실정이다. 따라서 본 연구는 Landsat OLI영상을 이용하여 식생피복지수 산정방법을 연구하였다. 정확하고 효율적인 식생피복지수 산정방법을 연구하기 위하여, 본 연구에서 제안된 방법과 선행연구방법을 비교평가 하였다. 실험결과 NDVI와 식생피복지수는 많은 연관성을 지니는 것으로 분석되었으며, 본 연구에서 제안된 방법을 이용한 식생피복지수가 특이점을 제외한 RMSE 7.3%로 전체 방법 중에서 가장 높은 정확도를 보였다.

The Vegetation cover is a significant factor to comprehend characteristics of the ground surface for meterological and hydrological models, which measure energy in the atmosphere or predict the runoff of ground surface. Deardorff introduced vegetation cover fraction to quantitatively comprehend the vegetation cover in 1978. After Deardorff, most of previous researches were conducted on low-resolution or high-resolution images, but only few researches on Landsat that are in medium-resolution images. Therefore, this study aims to investigate a way of calculating the vegetation cover fraction by using NDVI of Landsat images, which were hardly handled previously. For accurate vegetation cover fraction, we compared the evaluated parameters from this study with past vegetation cover fraction parameters that have been calculated for using NDVI of Landsat OLI images. The result of research was shown that NDVI is quite correlated with the vegetation fraction cover in the previous researches. In fact, RMSE of vegetation cover fraction values that obtained through the suggested parameters on this study showed the highest accuracy of 7.3% among all the cases.

1. 서 론

지표면의 특성을 파악하는 것은 대기의 에너지를 측정하거나 지표면의 유출을 예측하는 수치모델을 구성하는데 있어 매우 중요한 요소이다. 기존 연구에서는 1978년 Deardorff가 이를 정량적으로 파악하기 위하여 식생피복지수(Vegetation Cover Fraction)를 제안하였다(Deardorff, 1978). 그러나 1990년대 이전에는 식생피복지수를 이용하여 수치모델을 구성하는데 있어 제약이 많아 식생피복지수의 활용은 크지 않았다.

이후 다양한 원격탐사 자료에 대한 취득이 가능해지고 처리가 용이해지면서 MODIS 및 AVHRR등의 위성영상을 이용 하여 식생피복지수를 수치모델에 적용하는 연구가 활발하게 진행되어 왔다(Gutman et al., 1997; Xubin et al ., 1999; Gallo et al., 2001; LI et al., 2003; Rahmani et al., 1997). 위성영상을 이용하는 식생피복지수는 지표면의 반사율을 측정하여 산정하며, 수치적으로 표현되는 산정 값을 이용하여 다양한 모델에서 식생의 수평 밀도를 측정하는 자료로 활용된다(Gutman et al., 1997).

뿐만 아니라, 정확한 식생피복지수는 잠열유량(latent heat flux)을 계산하는데 있어 매우 중요한 역할을 한다. 잠열유량은 기상상태, 토양재질, 잎기공저항(leaf stomatal resistance)들과 같은 수많은 특성을 모델링 하는데 중요한 변수이며, 여러 연구결과에서 식생피복지수는 잠열유량 값에 크게 영향을 미치는 것으로 분석되었다(Gallo et al., 2001; LI et al., 2003).

국내에서는 KOMPSAT-2위성의 NDVI와 현장식생피복상태를 비교분석결과 NDVI와 현장식생피복상태의 상관성이 큰 것으로 나타났다(Kim et al., 2011). 또한 국립재난안전연구원에서 SPOT5 NDVI와 MODIS영상을 이용하여 식생피복 지수를 구하는 식을 제안했지만, 수치모델에서 활용하고 있는 식생피복지수 산정방법인 GI모델에서는 LANDSAT을 이용한 식생피복지수산정에 관한 연구는 미비한 편이다. 따라서 본 연구는 Landsat OLI영상을 이용하여 본 연구에서 제안한 방법과 기존의 식생피복지수 방법과 비교하고, 이를 실제 현장데이터와 비교평가함으로써 국내 환경 모델링에 최적화된 식생피복지수 산정방법을 제안하였다.

 

2. 활용데이터 및 처리방법

2.1 Landsat OLI data

본 연구에서 매개변수를 산정하기 위한 위성영상은 2013년 4월부터 12월에 취득된 Landsat OLI영상(WRS-2: Path 114, Row 35)이다. Landsat OLI의 특징으로는 Fig. 1과 같이 기존의 Landsat ETM+에서 2개의 밴드(Ultra blue, Cirrus)가 추가 되었으며, Table 1과 같은 QA밴드를 지닌 BQA영상을 같이 제공하고 있다. 영상의 상태를 파악하기 위한 QA밴드는 위성 전체 픽셀의 각각의 센서상태, 대기상태, 지표면상태 등을 시각화하여 표현한다. 따라서, 이러한 특징을 가진 BQA영상을 활용한다면 노이즈 픽셀의 마스킹이 용이하며, 특정 상태의 영상 추출이 가능하다.

Table 1.Level 1 QA band of Landsat OLI

Fig. 1.Bandpass wavelengths for Landsat 8 OLI and TIRS sensor, compared to Landsat 7 ETM+ sensor(USGS)

2.2 전처리

다중분광영상 자료들을 효과적으로 활용하기 위해서는 보다 정확한 복사 및 대기보정 등이 필수적으로 요구되는데, 기존연구에서 대기 및 태양의 흡수, 굴절 등으로 인하여 Landsat TM영상의 band1,2,3 값은 일반적으로 20% 적게 band4는 20%이상 크게 나타나는 것으로 판명되었다 (Chavez et al., 1996). 따라서 이를 보정하기 위한 대표적인 알고리즘으로 QUAC 또는 MODTRAN을 이용한다. MODTRAN의 결과는 QUAC의 결과보다 향상된 결과가 도출되나 MODTRAN의 경우 센서높이, 비행시간, 대상지역의 특성 등을 완벽하게 알고 있어야 한다는 제약이 있다(Agrawal et al, 1996). 본 연구에서는 USGS에서 제공된 특성값을 적용하여 FLAASH 대기보정을 수행하였으며, (1)식과 같다.

ρ: pixel surface reflectance ρe: average surface reflectance for the surrounding region S: spherical albedo of the atmosphere La: radiance backscattered by the atmosphere without reaching the surface A,B: surface independent coefficients that vary with atmospheric and geometric conditions.

대기보정이 끝난 위성영상의 노이즈를 제거하기 위하여 Landsat OLI에서 제공되는 BQA영상을 이용한다. QA밴드를 이용하여 노이즈를 마스킹한 영상을 추출하였고, 추출된 영상은 동일 지역의 다중시기 영상이다. 이와 같이 동일 장소에 대하여 여러 시기에 촬영된 영상의 NDVI 값은 다양한 요소에 의하여 불안정한 값을 포함하므로 이를 제거하기 위하여 최대값 합성 알고리즘(Maximum Value Compositing Algorithm)을 이용한다(Xubin et al., 1999). 따라서 본 연구 에서도 최대값 합성 알고리즘을 이용하여 하나의 NDVI 영상을 제작하였다.

2.3 식생피복지수

식생피복지수 산정을 위한 방법으로는 크게 1차 선형과 2차 선형방법이 있으나 본 연구에서는 1차 선형모델이 NOAH Land-Surface Model과 NAM Eta model을 제작함에 있어 공통적으로 사용되는 모델이기 때문에 1차 선형모델을 기준으로 대상지역을 분석하였다(Xubin et al., 1999; Gallo et al., 2001; LI et al., 2003). 1차 선형모델인 Gutman and Ignatov(GI)모델에서의 식생피복지수 결정은 (2)식과 같다.

NDVIPIXEL : pixel NDVI value NDVIo: bare soil NDVINDVI∞: live vegetation NDVI

NDVIo과 NDVI∞의 대표적인 산출방법으로 1)GI모델에서 제안하는 NDVIo=년중 비식생지역에서의 최저값, NDVI∞=년중 식생지역에서의 최대값, 2)Zeng이 제안하는 각 피복별 퍼센트타일을 이용, 3)Gallo이 제안하는 특정지역, 특정월의 5년 평균값의 퍼센트타일을 이용 4)Matsui이 제안하는 20년 최대, 최소 값의 퍼센트 타일을 이용하는 방법 등 이 있다.

2.4 개선된 NDVIo과 NDVI∞ 산정방법

GI모델을 적용한 식생피복지수를 산정하기 위하여 국외의 경우 기존에 제안된 매개변수를 활용하거나 International Geospere-Biosphere Program(IGBP)를 이용하였고, 국내의 경우 항공사진을 이용하여 비식생지역과 식생지역을 디지타이징하여 이를 구하였다. 그러나 국외의 방법은 Table 2와 같이 매개변수 산정을 위하여 저해상도 영상을 사용하였고, 산정된 매개변수 또한 국내와 환경이 다른 지역을 대상으로 연구되었다. 국내의 경우 국외의 방법보다 정확한 값을 산 정할 수 있으나 대상지역을 구분하는데 있어 많은 시간과 비용이 들어 대규모 지역의 수치모델을 구성하는데 있어 어려 움이 있다.

Table 2.Study area and dataset

본 연구에서는 이러한 문제점을 해결하고자 2가지 방법을 적용하였다. 첫째, 토지피복분류를 빠르고 정확하게 하기 위하여 환경부 정밀토지피복도를 이용하였다. 분류시 대상지역에서 수계와 같은 비식생지역을 제외하고 특이점을 제거하기 위하여 각 중분류별 95% 신뢰구간을 이용하여 재추출하였다. 둘째, GI모델에서 사용한 클러스터분석(Gutman et al., 1997)의 정확도를 높이기 위하여 EM클러스터를 이용하여 분석을 수행하고 분석된 결과값을 조정하기 위하여 각 중분류별 포인트 수를 이용한 가중치 분석을 수행하였다(Gallo et al., 2001).이와 같이 분석된 값을 이용하여 매개변수 삼각망을 구성하였다.

 

3. 실험 및 분석

본 연구에서는 Landsat OLI영상과 정밀토지피복도의 중분류지도를 이용하여 연구하였고, 흐름도는 Fig. 2와 같다.

Fig. 2.Flow chart of vegetation fraction map of study area

3.1 연구대상지역 및 데이터

Landsat 8 OLI는 2013년 4월부터 영상을 제공하고 있어, 아직 많은 영상을 확보할 수가 없다. 따라서 매개변수를 산정하기 위한 대상지역으로 양호한 영상을 많이 확보할 수 있는 경상북도(35°25'47.82"-37°02'13.34"N, 128°02'44.77"- 130° 34'30.54"E) 대구를 중심으로 90×90km를 선정하였으며, 사용 영상의 취득날짜와 구름의 상태는 Table 3과 같다. 또한 산정된 NDVI0과 NDVI∞를 이용한 식생피복지수를 비교하고자 114개의 현장데이터와 경기도 OLI 영상을 이용하였다. 아래 Fig. 3의 (a)는 매개변수 산정을 위한 경상북도의 LANDSAT OLI 영상이며, (b)는 현장데이터와 비교평가를 하기 위한 경기도 지역의(WRS-2:Path 116, Row 34, 2013- MAY-11) LANDSAT OLI 영상이다. 현장데이터를 구하기 위하여 항공사진의 식생픽셀을 계산하고(Purevdorj et al., 1998) 이를 보완하고자 감독자 검수도 같이 적용하였다.

Fig. 3.Scenes of Landsat OLI

Table 3.Landsat OLI scenes used to map in the study area

3.2 전처리

첫 번째로 동일 지역의 다중시기에 획득된 전체 영상의 NDVI를 얻기 위하여 각 BQA영상을 이용하여 노이즈를 마스킹하여 추출하고 이후 최대값 합성 알고리즘을 이용, 하나의 NDVI값을 획득하였다. 그 결과로 얻어진 마스킹 영상과 NDVI영상은 Fig. 4와 같다.

Fig. 4.Result of NDVI of Landsat OLI : (a) RGB image, (b) BQA image, (c) Masking result(yellow is no data), (d) NDVI result with MVC.

위에서 구해진 NDVI영상을 이용한 식생피복지수 산정은 기존의 연구에서 활용된 IGBP의 토지피복분류를 참고하여 환경부 정밀토지피복도를 사용하였다. 정밀토지피복지도를 활용한 광역적 범위의 피복 분류는 기존 연구에서 항공사진이나 위성영상을 이용하여 특정 피복분류를 추출하는 것보다 식생피복지수를 구하는데 있어 보다 효과적인 것으로 나타났다. 토지피복분류 결과 연구대상지역은 숲이 전체에서 66%를 차지하고 있으며, 가장 두드러진 특성으로 숲을 제외한 논이 11%로 가장 큰 부분을 차지하고 있는 것으로 나타났다. 이는 국내의 지역적인 특성을 반영하고 있는 것으로 판단되며, 분류결과는 Table 4와 같다.

Table 4.Pixel numbers [relative to the total pixel number of 8.7×108over land (I.E, categories 110-620)

본 연구에서 제안하기 위한 식생피복지수 매개변수산정과 대상지역의 최대 NDVI영상과 최소 NDVI영상을 제작하기 위하여 4월에서 12월까지 NDVI영상을 구하고, 각 영상을 최대값 합성 알고리즘을 이용하여 Fig. 5와 같은 하나의 NDVI 영상으로 제작하였다. NDVI영상에서의 특이점으로는 Table 5에서 보는 것과 같이 논의 경우 4월은 습지보다 높지만 나지와 비슷한 값을 모내기가 시작되는 5월은 습지보다 낮은 값이 나오는 것을 알 수가 있다.

Table 5.Comparison of NDVI in APRIL and MAY

Fig. 5.Results of minimum and maximum NDVI

이는 외국과 다르게 국내는 논의 비율이 전체에서 많이 차지하고 있기 때문에 특정 계절에서 수치모델을 적용할 경우 (예, 대기모형, 유사유출모형 등)기존 피복도별 절대값을 이용한 방법의 경우 많은 오차를 발생할 수 있다는 것을 확인할 수 있다.

3.3 식생피복지수

ZENG의 1999년 연구에서 AVHRR영상과 LANDSAT TM영상을 비교한 결과 초지와 관목지에서는 AVHRR이 크게 나오나 산림과 사막지역에서는 작게 나오는 양상을 보였다. 또한 1,3,5m급 항공사진과 비교에서 LANDSAT TM데이터 간의 상관성이 크다는 것을 알 수 있었다. 또한, 기존의 선행연구에서 정해진 지수를 사용하는 것이 아닌 히스토그램을 이용한 지수선정방법이 보다 효과적인 것으로 나타났다 (Xubin et al., 1999). 국내 연구에서 위성영상의 NDVI와 현장 식생피복이 높은 상관성을 가지는 것으로 나타났다(Kim et al., 2011).따라서 본 연구에서는 기존의 연구방법 중 대표적인 방법을 통하여 산정된 매개변수와 LANDSAT 8 OLI NDVI를 이용하여 제안한 매개변수를 적용하여 대상지역을 분석 하였다.

기존 연구방법은 9개의 케이스 별로 분류하여 적용하였으며, 분류된 내용은 아래 Table 6과 같다. 식생피복지수를 구하기 위하여 첫 번째로 Table 6에서 매개변수 결정방법을 각 중분류별 토지피복에 적용하여 NDVI0와 NDVI∞을 산정하였으며, 각 케이스별 대표 중분류(110, 300, 600번대)의 값은 Table 7과 같다. 그 결과 NDVI0는0.04~0.297의 값을 가지며, 케이스 6의 경우 다른 케이스와 다르게 전체 영역에서의 3th 퍼센트타일을 적용하여 가장 큰 값을 가지게 된 것으로 판단된다.

Table 6.Method of FG in recent publications

Table 7.Results of parameters of vegetation cover fraction

기존의 매개변수 산정방법 외에 본 연구에서는 국내실험에 적합한 매개변수를 산정하기 위하여 대상지역에서 수계와 같은 비식생지역을 제외한 8,782,088개의 점을 추출한 후, 이를 활용하여 클러스터분석을 수행하였다. 클러스터 분석 방법으로는 각 픽셀의 최대값과 최소값을 선택하여 전체 픽셀에 대하여 EM클러스터분석(7개의 클러스터)을 실행하였다. 분석결과값을 조정하기 위하여 각 중분류별 포인트 수를 이용한 가중치 분석을 수행하였다(Gallo et al., 2001). 그 결과 비식생지역의 클러스터는 0.14

이와 같은 7개의 클러스터 결과값을 도식화 하여 삼각형으로 표현한 그림은 Fig. 6과 같다. 이것은 대상지역에서 계절 변화지역의 NDVI는 단기간에 NDVI의 변화 즉, 식생의 변화가 크다는 것을 보여주고 있으며, 단 기간에 변하는 NDVI는 매개변수로 정하기에 식생피복지수의 결과오차가 매우 클 수 밖에 없다. 또한, 100% 식생지역과 순수한 토양지역을 찾는 것이 어려우므로 변화의 변동 폭이 적은 비식생지역의 최소값을 NDVI0 (0.14)로 선정하고, 식생지역의 최대값을 NDVI∞ (0.86)로 선정하였다.

Fig. 6.Results of cluster analysis using NDVI’s maximum and minimum

제안한 방법을 10번째 케이스로 적용하여 식생피복지수를 산정한 결과 케이스 6의 식생지수가 가장 낮게 나왔으며, 이는 의 값이 0.269로 가장 크게 산정되었기 때문으로 판단된다. 또한 고정 매개변수의 값의 경우 대체적으로 과다하게 식생피복지수가 산정된 것으로 보인다. 이는 Landsat 영상의 값이 아닌 MODIS영상과 같은 다른 특성을 가진 영상의 매개변수 값을 이용했기 때문인 것으로 나타났으며, 전체 중분류별 값은 Fig. 7과 같다.

Fig. 7.Mean of vegetation cover fraction

위의 결과 값에서 Fig. 8과 같이 300번대의 고정 매개변수의 경우 식생피복지수가 과다하게 산정된 것을 볼 수 있으며, 전체 페센트타일을 쓴 경우 과소하게 산정된 것을 알 수 있다.

Fig. 8.Mean of vegetation cover fraction in 300

Fig. 9.Mean of vegetation cover fraction in 100

위의 300번대의 식생피복지수의 결과보다 비식생지역인 100번대의 경우 그 결과의 차이가 크게 나타나며, 전반적으로 히스토그램을 이용한 값과 기존 연구 매개변수를 이용한 값과의 차이가 2배 이상 나는 것으로 분석되었다.

위의 결과 값에서 비식생지역 중 140번대의 경우 다른 지역에 비하여 높은 식생피복지수를 보이고 있다. 이는 140번대의 경우 문화/체육/휴양지역을 나타내는 피복이며, 대부분 운동장으로 식생이 충분히 자랄 수 있었던 환경으로 분석 되었다.

3.4 케이스별 식생피복지수 비교

식생피복지수는 각 픽셀에 대하여 독립적으로 취득이 되기 때문에 항공사진에서 감독자가 분류한 식생지수와 비교평가가 가능하다. 이를 위하여 가장 비슷한 시기의 정사항공사진을 취득할 수 있는 지역을 선정하였으며, 경기도 일대의 0.52m급 정사항공사진 9도엽을 취득할 수 있었다. 이를 활용하여 위성영상에서 114개의 테스트 필드를 선정하고 이를 항공사진에서 30×30m 그리드 격자를 생성 항공사진의 식생픽셀을 계산하고(Purevdorj et al., 1998) 이를 보완하고자 감독자의 검수에 의하여 식생피복을 구하였다. 114개의 테스트 필드는 식생피복지수를 정확하게 비교하기 위하여 경계지역을 위주로 선정하였으며, 식생피복구분이 난해한 일부 산림지역과 물이 차있는 논지역은 제외하였다. 물이 차 있는 논의 경우 식생피복지수 값이 수계 혹은 습지와 유사하게 나오므로 선행연구에서와 같이 제외하고, 식생피복지수가 85%가 넘는 산림의 경우 항공사진을 이용한 식생피복지수 결정이 난해하여 제외하였다.

Fig. 10.Results of fraction vegetation and aerial photograph(0.52m)

테스트필드 114개 중 특이점 오차가 12개 발생하였고, 이는 2012년에 촬영된 항공사진과 위성영상의 시기적 차이에 의하여 건물의 신축 혹은 식생생장에 의하여 발생한 것으로 분석되었다. 각각의 케이스별 특이점 제거 전·후 RMSE값은 Table 7과 같다.

Table 7은 보는 바와 같이 특이점을 제거하기 전에는 제안매개변수를 활용한 방법과(케이스 10) Oleson제안식의 RMSE가 11이하로 가장 좋은 결과를 나타냈으며, 특이점을 제거한 후에는 본 연구에서 제안한 RMSE가 7.3으로 가장 좋은 결과를 나타내고 있다. 이와 같이 제안된 매개변수를 이용한 방법은 기존의 제안 방법보다 국내 지형에 보다 적합하게 적용하여 수치모델의 정확도를 높일 수 있을 것으로 분석되었다.

Table 8.RMSE of each methods and used to exclude outliner

 

4. 결 론

본 연구에서 제안된 방법을 GI모델에서 적용한 결과 다음과 같은 결론을 얻을 수 있었다.

첫째, 기존의 연구에서 제안된 고정된 식생지수를 사용하는 경우보다 히스토그램을 이용한 방법(케이스3~4, 케이스 7~8)의 정확도가 10%가 높았으며, 본 연구에서 제안된 방법의 경우 GI모델 고정매개변수보다 3배 이상 정확도가 높은 것으로 나타났다.

둘째, 정밀토지피복도를 이용하여 연구지역의 피복을 분류하는 방법은 광범위한 지역에서의 피복분류를 실행할 경우 기존의 항공사진 혹은 위성영상을 통하여 비식생지역과 식생지역을 구분하는 방법보다 효율적이며, 논에 물이 차거나 토지이용도가 변한 부분에 대해서는 기존 정밀토지피복도를 이용하여 피복을 분류한 경우 파악을 하기가 어려운 것으로 나타났다.

셋째, 연구 대상지의 특정 지역에서는 항공사진과 비교하여 식생피복지수가 다소 높거나 낮게 나오는 지역이 있었으나 그런 지역은 기존의 연구(Kim et al., 2011)에서 나온 결과와 같이 위성 센서에서 수집되는 키 작은 식물군에 의하여 근적외선 반사치량이 적거나 나무군락의 잎에 의하여 식생피복의 양에 비하여 과다하게 산정되는 것으로 보여졌다. 하지만, 위와 같은 오차에도 불구하고 위성영상을 이용하여 얻은 식생피복지수와 항공사진을 이용하여 얻은 식생피복지수는 높은 상관성을 가지는 것으로 판단되었다.

넷째, 본 연구에서 제안된 방법을 이용하여 수치모델에 적용한다면 특정 계절에 많은 오차를 포함하는 지역을 가진 국내 환경에 맞는 정확하고 정밀한 기상 및 수문 모델링이 가능 할 것이다.