Forest Fire Severity Classification Using Probability Density Function and KOMPSAT-3A

확률밀도함수와 KOMPSAT-3A를 활용한 산불피해강도 분류

Abstract

This research deals with algorithm for forest fire severity classification using multi-temporal KOMPSAT-3A image to mapping forest fire areas. The recent satellite of the KOMPSAT series, KOMPSAT-3A, demonstrates high resolution and multi-spectral imagery with infrared and high resolution electro-optical bands. However, there is a lack of research to classify forest fire severity using KOMPSAT-3A. Therefore, the purpose of this study is to analyze forest fire severity using KOMPSAT-3A images. In addition, this research used pre-fire and post-fire Sentinel-2 with differenced Normalized Burn Ratio (dNBR) to taking for burn severity distribution map. To test the effectiveness of the proposed procedure on April 4, 2019, Gangneung wildfires were considered as a case study. This research used the probability density function for the classification of forest fire damage severity based on R software, a free software environment of statistical computing and graphics. The burn severities were estimated by changing NDVI before and after forest fire. Furthermore, standard deviation of probability density function was used to calculate the size of each class interval. A total of five distribution of forest fire severity were effectively classified.

본 연구는 산불 전후 KOMPSAT-3A 영상을 사용하여 산불피해지역을 분석하는 것을 목적으로 한다. KOMPSAT 시리즈 중 KOMPSAT-3A는 적외선 및 고해상도의 멀티 스펙트럼 밴드를 가진 VHR위성이다. 하지만, KOMPSAT-3A를 활용하여 산불피해강도를 분류하는 연구는 부족한 실정이다. 따라서 본 연구에서는 KOMPSAT-3A의 산불 피해강도를 분류하기 위한 새로운 알고리즘을 제시하는 것을 목표로 한다. 또한, 본 연구에서는 산불 피해지역에 대한 참조자료로 Sentinel-2로 생성한 dNBR을 사용하였다. 본 연구의 연구 지역은 2019년 4월 4일 강릉에서 발생한 산불 피해지역으로 선정하였다. 본 연구에서는 산불피해구간을 산정하기 위한 알고리즘으로 오픈 소스 통계 프로그램인 R software의 확률분포함수를 사용하였다. KOMPSAT-3A에서 산불 피해지역은 산불 전, 후 NDVI의 변화에 따라 생성되었다. 산불피해강도는 분포 함수의 표준 편차를 사용하여 각 등급 크기를 산정하였다. 총 5개 구간에 따른 산불 피해 강도가 효과적으로 분류되었다.

1. 서론

우리나라 산림에서 발생하는 산불은 그 피해 면적이 넓을 뿐만 아니라 봄철에 집중되어 발생하기 때문에 피해 지역의 면적을 신속하게 분석하여 대책을 마련할 기반 자료 구축이 필요하다. 산불에서 발생한 화재는 항공사진을 통해 분석되기도 한다(Harris et al., 2011). 하지만, 특정 임분에서의 산불위험도는 잠재된 연료량이나, 풍속, 지형 등에 따라 다양한 형태로 산불피해강도에 영향을 미치기 때문에 이에 대한 피해 특성을 규명하기 위해 위성영상자료는 필수적으로 활용된다 (Won et al., 2014). 산불피해강도는 산림의 피해 면적, 산불 전후의 피해 특성, 훼손된 산림의 종류를 파악하는데 중요한 기초자료로 활용될 수 있기 때문에 정확한 평가가 필요하다(Lentile et al., 2006).

인공위성 영상자료를 이용하는 것은 도시 전체나 국가 규모와 같은 넓은 지역을 효과적으로 관리하고 분석할 수 있기 때문에 산림피해조사와 관리 및 식생 변화 모니터링 분야에서 사용되고 있다(Kim et al., 2010; Jung et al., 2012). 인공위성의 적외선 파장대는 서로 다른 스펙트럼 지수를 통해 화재가 발생한 지역의 토지 식별에 사용될 수 있다(Chuvieco et al., 2002).

일반적으로 산불은 정규화된 식생 지수를 통해 파악할 수 있지만(Tucker, 1979), 산림뿐만이 아니라 모든 식생 지역에 대한 정규화 지수이기 때문에 산불의 특성을 고려한 산불 지수가 필요하다. 산불 피해는 일반적으로 NBR과 같이 Short wave infrared(SWIR)밴드를 포함한 분광밴드의 조합에 따라 나타낸다(Key and Benson, 2006; Wu et al., 2015). 다시기 위성영상을 활용한 NDVI는 식생 변화 지역을 탐지할 수 있기 때문에 산림 수종분류를 통해 산불에 영향을 미치는 요소를 분석할 수 있다(Lee, 2009). 산림의 종류는 침엽수, 활엽수 혹은 수종에 따라 나타낼 수 있으며 이를 분류하기 위해 기계학습 분류 기법을 적용한 연구가 있다(Wessel et al., 2018). 또한, 영상분류기법이 고도화됨에 따라 최근의 연구들은 산불 탐지를 위해 전통적인 감독분류기법뿐만 아니라 객체 기반 분류 기법과 PCA 등을 활용하기도 한다 (Katagis, 2014; Choi, 2015; Sertel and Alganci, 2016). Sentinel-2와 같이 다중분광밴드로 획득한 위성영상은 새로운 산불지수 알고리즘을 개발하여 산불 탐지에 활용되었다(Filipponi, 2019; Park et al., 2019). 산불탐지를 위한 고해상도 위성영상 중에는 R, G, B NIR 4개 밴드만 제공하는 Very-High Resolution(VHR)위성영상을 활용한 사례가 있다. SWIR밴드가 없는 위성영상은 정규화된 산불 지수가 없기 때문에 다양한 정규식생지수와 산불 피해지역의 상관분석을 통해 연관성을 나타내는 연구도 진행되고 있으나(Meng et al., 2017), 이러한 분석결과는 참고하는 수준에 머무르고 있는 실정이다. 본 연구는 국내에서 개발한 위성 중 KOMPSAT-3A를 활용하여 VHR 위성영상을 활용한 산불 피해 강도를 분석하는 데 연구 목적이 있다. 또한, KOMPSAT-3A의 공간해상력과 분광해상력을 활용하여 정밀한 산불피해지도를 제작하는데 본 연구의 성과가 기여할 수 있는 기초연구를 수행하였다.

2. 연구 자료 및 방법

본 연구에서는 KOMPSAT-3A의 위성영상을 바탕으로 산불피해강도를 나타내기 위해 2019년 4월 4일 산불이 발생한 강원도 강릉시 일대를 연구대상지로 선정하였다(Fig. 1).

Fig. 1. Geographical location and false color image of KOMPSAT-3A, Sentinel-2 in study site - Gangneung city, Gangwon province, South Korea. (a) K3A March 8th, 2019 (b) K3A April 12th, 2019 (c) Sentinel-2 March 1st, 2019 (d) Sentinel-2 April 15th, 2019.

우리나라는 공간해상력이 흑백 0.55 m, 컬러 2.2 m급으로 2015년 3월에 발사된 고해상도 위성영상 아리랑 3A호(KOMPSAT-3A)를 보유하고 있다. 특히, 가시광선 밴드를 비롯하여 적외선(Near-infrared)밴드를 제공하고 있어 원격탐사 분야 전반에 걸쳐서 활용되고 있다 (Yong et al., 2016). 따라서 본 연구에서는 산불 피해를 나타내기 위해 Fig. 1과 같이 산불 피해 전후에 촬영된  KOMPSAT-3A 위성영상을 활용하였다.

본 연구는 KOMPSAT-3A를 활용하여 산불피해강도를 재분류 하는 것을 목적으로 하며 전체적인 흐름은 Fig. 2에 따른다. 본 연구에서는 산불 피해 지역에 대한 참조 자료를 생성하기 위해 Sentinel-2 위성영상을 사용하였다. Table 1과 같이 Sentinel-2 위성은 다중분광센서를 탑재하고 있으며 특히, 20 m의 공간해상력을 가진 SWIR밴드와 10 m의 공간해상력을 가진 Red, Green, Blue, NIR 밴드를 제공하고 있다.

Table 1. Data source and list of spatial information

Fig. 2. Flowchart showing the methodology of the study.

Fig. 3은 Sentinel-2 위성영상의 NIR 및 SWIR 밴드를 사용하여 정규탄화지수(Normalized Burn Ratio: NBR)를 생성하였고, 산불 피해 전후 영상에서 dNBR(Differenced NBR)을 통해 산불피해 규모를 나타냈다(Table 2). 미국 지질조사국(USGS)에서 제공하는 dNBR 수치에 따른 산불 피해 기준은 High, Moderate-high, Moderate-low, Low, Unburned까지 총 5단계로 나뉜다. 이때 각 피해영역에 대해 산불 이후 KOMPSAT-3A 영상도 함께 나타내었다. 본 연구에서는 dNBR을 통한 산불피해지역과 비교하기 위해 정규식생지수(Normalized Difference Vegetation Index: NDVI)와 dNDVI(Differenced NDVI)를 사용하였다. NDVI는 녹색 식물의 특성상 높은 반사율을 나타내는 NIR 밴드와 낮은 반사율을 나타내는 Red 밴드의 조합으로 만든 전통적인 식생지수이다.

Fig. 3. Definitions of five burn severity levels in study area. (a) the false color image is Sentinel-2 scene in April 15, 2019. (b) burn severity distribution map of dNBR between March 1st, 2019 and April 15th, 2019 for Sentinel-2 images.

Table 2. Spectral indices evaluated in this study

또한, dNDVI는 산불발생 이전 시점의 NDVI값과 산불 발생 이후 시점의 NDVI값의 차를 의미하기 때문에 산불 피해영역을 나타내는 데 유용하다고 판단하였다. 일반적으로 NBR이 NDVI보다 산불피해지의 공간적 이질성과 복잡성을 더 효과적으로 표현할 수 있으며, 가을철에 촬영된 영상보다 봄철 생장기에 촬영된 영상이 피해 강도를 더 효과적으로 평가할 수 있다(Key and Bensen, 2006). 하지만 Fig. 3의 false color 영상은 산불피해가 나타나지 않았지만, dNBR에서는 산불피해가 발생한 지역이 있어 오탐지 되는 영역이 존재한다. 또한, 국내기술로 제작한 VHR 위성영상 중에는 SWIR 밴드를 가진 위성이 없어 NBR을 생성하기 힘들기 때문에 KOMPSAT-3A 영상의 NDVI를 통해 산불피해 구간을 산정하였다. 본 연구에서는 NDVI에 따른 산불 피해 구간을 나타내기 위해 R software의 확률밀도함수를 활용하였다. 확률밀도함수는 확률론에서 확률 변수의 분포를 나타내는 함수로 특정 구간에서 입력된 변수의 크기를 나타내는 데 활용된다.

\(P(a \leq x \leq b)=\int_{a}^{b} f(x) d x\)       (1)

식 (1)은 확률밀도함수의 수식을 나타낸 것이며 확률변수 X의 특정 구간인 a, b에서 넓이를 구하는 것이다. 본 연구에서는 산불 발생 시기가 4월이기 때문에 활엽수림의 식생이 증가하는 것을 고려하기 위해 중분류 토지피복지도를 사용하였다. 산림 수종이 산불 전후에 NDVI가 변화하는 구간을 확률밀도함수를 통해 나타내고 산불 피해가 발생한 영역과 산불 피해가 발생하지 않은 영역으로 구분하였다. 본 연구에서는 산불 피해가 발생한 영역에 대해 육안 검사에 의한 3가지 패턴으로 나누어 피해 구간을 산정하였다. 산불 피해지역은 완전 연소된 지역, 일부 피해가 받은 식생이 남아 있는 지역, 온전한 식생과 연소된 식생이 공존하는 지역으로 나누어 각각의 지역에서 확률밀도함수를 나타내고 표준편차를 산정하였다. 표준편차 결과들은 모두 합산하여 평균치를 구하고 산불 피해 구간의 계급 크기로 사용하여 최종적으로 dNDVI에 따른 산불 피해도를 작성하였다.

3. 연구 결과

Fig. 4는 연구대상지의 산불 피해지역을 나타내기 위해 NDVI를 사용하여 KOMPSAT-3A 위성영상을 나타낸 것이다. Fig. 4(a)는 산불이 발생하기 이전인 3월 8일 NDVI를 나타낸 것이고, Fig. 4(b)는 산불이 발생한 이후인 4월 12일 NDVI를 나타낸 것이다. 또한, Fig. 4(c)은 산불 피해 전후 NDVI 변화를 나타낸 dNDVI를 나타낸 것이다. dNDVI의 결과는 산불피해 지역에서 식생지수의 변화를 명확히 확인할 수 있다.

Fig. 4. Delta-Normalized Difference Vegetation Index with KOMPSAT-3A image (a) NDVI March 8, 2019 (b) NDVI April 12, 2019 (c) dNDVI between March 8 and April 12, 2019.

Fig. 5(a)은 연구 지역에서 산림 영역만 추출하여 Sentinel-2 영상의 dNBR 영역과 KOMPSAT-3A 영상의 dNDVI영역을 중첩한 것이다. dNBR은 다섯 단계로 산불 피해지역을 나타낼 수 있고 피해강도의 등급을 피해상(High), 피해중상(Moderate-high), 피해중하(Moderatelow), 피해경(Low), 피해무(Unburn)으로 구분하였다. 또한, 각 피해 영역마다 KOMPSAT-3A와 중첩하여 화소들의 dNDVI값과 빈도수를 계산하여 확률밀도함수로 나타냈다. Fig. 5(b)는 dNBR을 통해 나타낸 피해영역에 따라 dNDVI값과 빈도 수를 히스토그램과 확률밀도함수를 통해 구간에 따른 차이를 나타낸 것이다. 또한, Fig. 5(c)는 각 피해 영역마다 Sentinel-2 영상의 dNDVI를 확률밀도함수로 나타낸 것이다. Fig. 5(b)는 각 구간에서 최빈도수를 나타낸 피해도를 중심으로 dNDVI값에 따른 산불 피해도를 나타냈다. 5개의 확률밀도함수가 중첩되는 지점은 각각 0.04, 0.11, 0.24, 0.39로 나타났다. Fig. 5(c)에서 산불피해지역과 피해무(Unburn)지역의 확률밀도함수가 중첩되는 지점이 상대적으로 Fig. 5(b)보다 낮은 곳에서 나타난다. 이는 KOMPSAT-3A가 Sentinel-2보다 높은 공간해상력을 가졌기 때문에 나타난 현상으로 판단된다.

Fig. 5. Map showing the forest fire severity distribution generated with KOMPSAT-3A and Sentinel-2 images of forest areas. The histogram is representation of the dNDVI according to the dNBR burn severity distribution area. (a) Burn severity distribution overlay map of dNBR and dNDVI. (b) dNDVI probability density function for each severity area generated with KOMPSAT-3A histogram. (c) dNDVI probability density function for each severity area generated by Sentinel-2 histogram.

Fig. 6은 산림 종류에 따라 dNDVI의 확률밀도함수를 나타낸 것이다. 본 연구의 전체 산림 지역은 산불피해가 발생하지 않은 영역이 다수 분포하고 있기 때문에 dNDVI 값이 0에 가까운 화소들이 높게 분포하는 것을 확인할 수 있다.

Fig. 6. NDVI probability density function of forest species.

Fig. 7(a)는 산불 전후 산림 수종에 따른 NDVI 확률밀도함수를 표현한 것이다. 확률밀도함수에서 활엽수림은 일부 식생이 감소한 지역이 있으나, 식생이 증가한 지역도 있다는 것을 확인할 수 있다. 이는 본 연구의 연구지역에서 활엽수림의 산불 피해가 다른 수종에 비해 크지 않고, 봄철 시기상 1개월이 경과함에 따라 식생이 증가하게 됐다는 것을 나타낸다. 활엽수림과 침엽수림에서는 산불 피해 이후 NDVI 확률밀도함수가 0 ~ 0.2구간에서 급증한 것을 확인할 수 있다. 따라서 이 구간을 산불 피해가 발생한 영역의 NDVI 계급값으로 선정하였다. Fig. 7(b)은 NDVI 커널 밀집도로서 산불 전후 NDVI 변화를 나타냈다. 활엽수림은 확률밀도함수에서도 확인했듯이 NDVI가 증가한 곳과 감소한 곳이 공존하나 변화 크기가 상대적으로 작다. 하지만 침엽수림에서는 산불 이후 NDVI 값이 감소한 곳을 확인할 수 있다.

Fig. 7. Probability density function and kernel density plots of NDVI according to tree species. In density plots, the x-axis is KOMPSAT-3A NDVI, March 8th, and y-axis is KOMPSAT-3A NDVI, April 12th. (a) Probability density function of NDVI (b) Kernel density plots.

Fig. 8은 산불 피해가 발생한 지역내에서 육안검사를 통해 3가지 경우에 따른 NDVI 확률밀도함수를 나타낸 것이다. Area 1은 식생이 손상된 산림과 완전 연소된 산림이 존재하는 지역이다. 확률밀도함수에서도 0 이상부터 상승하다가 0.4 ~ 0.5에서 가장 높은 확률밀도함수를 보였다. Area 2는 완전 연소된 산림만 존재하는 지역으로, 0.5 에서 가장 높은 확률밀도함수를 보였으며, Area 1보다 상승 폭이 가파르다. Area 3은 온전한 식생이 있는 산림과 연소된 산림이 공존하는 지역이다. 확률밀도함수에서도 dNDVI 산불 이후 상승한 영역과 감소한 영역이 모두 있다는 것을 확인할 수 있다. 산불 피해 구간은 Fig. 8의 3개 영역에서 나타낸 확률밀도함수의 표준편차를 통해 산정하였다.

Fig. 8. Probability density function according to three patterns of forest fire areas. Area 1 is moderate distributed forest fire area. Area 2 is highly damaged area and Area 3 is low damaged area.

Area 1보다 상승 폭이 가파르다. Area 3은 온전한 식생이 있는 산림과 연소된 산림이 공존하는 지역이다. 확률밀도함수에서도 dNDVI 산불 이후 상승한 영역과 감소한 영역이 모두 있다는 것을 확인할 수 있다. 산불 피해 구간은 Fig. 8의 3개 영역에서 나타낸 확률밀도함수의 표준편차를 통해 산정하였다.

Area 1의 표준편차는 0.17, Area 2의 표준편차는 0.22, Area 3의 표준편차는 0.16으로 나타났다. 각 표준편차는 모두 합산하여 평균치를 산정하여 0.18로 산정되었으며, 이 중 양의 편차와 음의 편차가 중복될 수가 있으므로 절반인 0.09를 산불 피해지역의 계급 구간으로 선정하였다. Fig. 9는 Fig. 8에서 선정한 계급 구간을 이용하여 KOMPSAT-3A의 dNDVI를 재분류한 결과이다. 총 5개 영역으로 나타냈으며 피해무(Unburn)등급에 대한 구분값은 Fig. 7에서 확률밀도함수가 중첩되는 지점을 이용하였다. 재분류 결과 Fig. 3의 dNBR에서 산불이 발생하지 않은 영역에서 산불이 발생했다고 하는 오류가 개선된 것을 확인할 수 있다.

Fig. 9. Reclassified dNDVI map with KOPMSAT-3A in forest fire area.

4. 결론 및 고찰

본 연구에서는 2019년 강릉 옥계면에서 발생한 산불의 피해 강도를 분류하기 위해 KOMPSAT-3A를 사용하였다. 산불 피해지역에 대한 참조자료는 Sentinel-2에서 생성한 dNBR을 사용하였고, dNBR 등급에 따른 5개 벡터 영역에 포함된 KOMPSAT-3A의 dNDVI를 확률밀도함수로 나타냈다. dNDVI의 계급 구간은 산불 피해지역에서 육안 판독을 통해 다른 형태를 보인 3개 영역을 추출하고, 각 영역의 확률밀도함수가 중첩되는 지점을 통해 계급 구간을 산정하였다. 총 5개 구간으로 산불피해도를 나누었으며, 이로부터 얻어낸 결과는 다음과 같다

첫째, 위성영상을 이용하여 산불 피해도를 산정하는 방법은 SWIR 밴드를 활용한 dNBR을 활용하는 것이 일반적이지만, 본 연구에서는 R, G, B, NIR 4개 밴드만을 가진 고해상도 위성영상으로 dNDVI를 통해 산불 피해지수를 나타냈으며, 효과적으로 산불 피해도 분류를 수행하였다.

둘째, KOMPSAT-3A 영상의 산불 전후 dNDVI를 활용하여 산불 영역을 dNBR로 나타냈을 때 산불이 발생하지 않은 지역에서 산불이 발생했다고 나타난 오류를 개선하였다.

일반적으로 산불 피해 분류는 SWIR을 포함한 dNBR을 이용하는 선행 연구가 진행되었다(Won et al., 2014). 하지만 본 연구에서는 산불 피해 강도를 분류하기 위해 SWIR이 없는 VHR 위성인 KOMPSAT-3A 위성을 이용했다는 데에 의의가 있다. 하지만 연구 지역의 산림 특성상 음영지역이 존재하여 사면에 따라 다른 분광값을 가질 수 있고, 강릉시 지역만을 대상으로 하였기에 일반화하기에 한계점이 존재한다. 따라서 향후 추가 연구를 통해 다른 연구 지역을 대상으로 지형,임상도를 적용한 산불피해의 정량화, 위성촬영시점, 경사면 등을 고려하여 KOMPSAT-3A 영상에서의 정규화된 산불 지수를 나타내고 범용적으로 활용할 수 있는 표준 모델을 제시하는 것을 향후 연구 목표로 삼고자 한다.