서론
산불은 다양한 교란 유형 중 빈도와 피해규모 면에서 현재 전 세계 산림생태계의 주요 교란 중 하나이다(Chandler et al., 1983; Whelan, 1995). 특히 산불로 인한 나무의 고사는 탄소 저장, 생물 다양성 보전, 수자원, 경제 및 사회적 서비스 등에 큰 영향을 미친다(Bowman et al., 2009). 고사한 나무는 강풍에 의해 도복되면서 주변의 도로나 전신주와 같은 인프라를 파괴하는 2차 피해가 발생할 우려가 있다. 또한 소나무 재선충병은 산불피해지에서 가까울수록 확산 가능성이 크며(Shim, 2019), 산불피해가 심하고 시간이 지날수록 매개충의 서식 밀도가 증가한다(Jung et al., 2020). 이렇듯 기후변화로 인해 산불이 점차 대형화되고 발생 빈도가 증가하면서 산불피해목의 처리를 결정하기 위한 과학적 근거가 더욱 필요해지고 있다.
기존에는 산불피해목 존치 및 벌채 여부를 산불 직후 현장 조사를 통한 교목들의 수관 피해 여부를 기준으로 결정하였다(Korea Forest Service, 2006). 최근에는 위성영상이나 무인기를 이용한 원격탐사로 얻어진 분광 지수들로 산불피해 등급을 심･중･경으로 분류하여 산불피해지 복원을 위한 근거로 활용하고 있다(Lee et al., 2012; Lee and Jung, 2019; Lee et al., 2023). 산불피해 등급이 ‘경’으로 분류되는 지역은 지표화 피해지로, 수관층의 피해가 거의 없는 지역을 말한다(Won et al., 2013). 일반적인 산불 피해지 복원에서는 산불피해 등급이 ‘경’인 지역의 피해목은 존치하고 있다(Lee et al., 2006).
소나무는 활엽수림에 비해서 수지분이 많고 발열량이 크기 때문에 낮은 온도에서도 산불발화의 위험성이 높다(Lee and Lee, 2006). 또한 수관층에 피해가 없고 수간에 그을음이 남을 정도의 산불피해에도 고사가 발생한다(Lee et al., 2008). 이러한 소나무의 특성을 고려하면, 복원계획 수립에서 산불피해도 ‘경’인 소나무에 대하여 일괄적인 존치 결정이 아닌 잠재적인 고사 위험성을 고려한 피해목 처리의 기준을 마련할 필요성이 있다.
산불 후 산림의 고사는 수관층 및 수간의 직접적인 연소에 의한 손상으로 인해 단기간에 발생하기도 하고(Dickinson and Johnson, 2001), 토양 및 식생의 변화와 병해충 피해, 기후에 의한 생리적 스트레스에 의해서 중장기적으로 발생하기도 한다(Van Mantgem et al., 2011; Das et al., 2016; Kane et al., 2017). 또한 고도, 경사, 사면 방향 등의 지형적인 요인도 임목의 고사에 영향을 미친다(Maringer et al., 2016). 산불피해에 따른 수목의 생존 연구는 다양한 수종들에 대해서 주로 수관 피해(Crown scorch), 수간 그을림(Bark char), 뿌리 열해(Heat to roots)와 같은 산불피해 정도와 수피 두께(Bark thickness), 흉고직경, 수고 등의 피해 수용 능력과 같은 다양한 인자들을 이용하여 이루어졌다(Michaletz and Johnson, 2007; Lawes et al., 2011; Brando et al., 2012). 또한 최근에는 원격탐사로 얻어진 여러 광학지수들을 이용한 수목의 고사예측 연구도 많아지고 있다(Rao et al., 2019; Furniss et al., 2020; Pascual et al., 2022).
생존분석은 환자의 사망과 같이 관심 있는 사건(Event)이 발생할 때까지 걸리는 생존시간(survival time)을 추정하는 통계적인 추론 방법이다. 생존분석의 대표적인 방법인 콕스 비례위험모형(Cox proportional hazard model)은 개별 대상의 생존시간에 영향을 미치는 공변량(covariate)들을 이용하여 그 관계성을 확인할 수 있는 통계 모형으로, 위험비(Hazard ratio)를 계산함으로서 공변량들의 영향력을 보여줄 수 있다(Cox, 1972). 콕스 비례위험모형은 연구 기간 동안 관찰된 대상에 대해 “위험비는 시간에 따라 변하지 않으며 항상 일정하다.”라는 비례위험 가정(proportional hazard assumstion)을 전제로 한다. 대형산불의 경우 피해면적 대비 한정된 인력과 자원을 이용하여 다년간에 걸친 복원계획을 수립하기 때문에, 피해목 벌채와 복원사업의 우선순위를 정하는 것이 필요하다. 생존분석은 산불피해 소나무의 생존여부 판별뿐만 아니라 고사발생에 소요되는 시간을 예측할 수 있다.
본 연구에서는 콕스 비례위험모형을 이용하여 산불피해 소나무의 고사에 영향을 미치는 인자를 밝히고, 시간에 따른 생존율 변화를 비교함으로써 산불피해목 고사에 미치는 영향력을 분석하고자 하였다.
재료 및 방법
1. 연구 대상지 및 대상목
강원도 삼척시의 기후 평년값(1991~2020)은 기온 13.3°C, 강수량 1509.7 mm, 풍속 2.4 m/s이며(Korea Meteorolgical Administration, 2024), 연구 대상지가 위치한 도계읍 일대(37°15′52.60″N/ 129°1′18.77″E)는 2017년 5월 6일에 산불이 발생하여 765 ha의 산림이 연소되었다(Korea Forest Service, 2018). 산불피해 지역 중 소나무 단순림이며 산불 피해 등급이 ‘경’으로 분류된 지역 3개소를 연구 대상지로 선정하였다(Figure 1). 대상지 내에 위치한 소나무들 중에서 수관층에 직접적인 연소로 인한 그을림이나 열해로 인한 잎의 갈변이 없으며, 지표화 피해로 인해 수간의 수피에 그을음이 남아 산불피해를 확인할 수 있는 소나무 340 그루를 대상목으로 선정하였다. 대상목들은 6영급의 소나무로 평균 수고는 17.7 m, 평균 흉고직경은 33.6 cm로 중경목에서 대경목에 해당 되었다. 2017년부터 2023년까지 매년 1회씩 대상목들의 고사 발생 여부를 확인하고 기록하여 모니터링하였다.
Figure 1. Location of the 3 research sites(left) in Samsheok, Korea and photograph of the burned forest(right).
2. 고사 영향 인자 자료 수집
1) 현장 조사 자료
현장 조사를 통해 대상목의 흉고직경과 수간의 수피에 남겨진 그을음의 높이와 상대적인 면적을 측정하였다. 그을음의 측정은 대상목을 사방면(동, 서, 남, 북)에서 각각 바라보고 해당 면에서 가장 높은 그을음의 높이와 그 높이를 가지는 가상의 사각형을 그렸을 때 그을음이 차지하는 비율을 기록하였다. 조사를 통해 얻어진 그을음의 높이는 평균을 내어 수피 그을음 높이(Bark Scorch Height, BSH)를 계산했으며, 그을음의 상대적인 면적은 각 면의 그을음 높이와 비율을 곱한 후 합한 수피 그을음 지수(Bark Scorch Index, BSI)를 산출하여 분석에 이용하였다(Kwon et al., 2021)(Figure 2).
Figure 2. Measurment and calculation of BSH and BSI.
H: Bark scorsh area height, R: Bark scorch area rate, BSH: Average of bark scorch height, BSI: bark scorch index
2) GIS 및 원격탐사 자료
각 대상목의 위치좌표와 대상지의 수치표고모델(Digital Elevation Model, DEM) 자료를 이용하여 각 대상목이 위치한 지형을 분석하였다. 수치표고모델은 국토지리정보 플랫폼에서 제공하는 1:5,000 수치지형도의 등고선을 이용하여 생성하였다. 결과적으로 대상지의 고도, 경사, 사면방향을 지형 인자로 선정하여 분석에 활용하였다.
산불 발생 전・후의 수관층 및 환경적 변화를 확인하기 위하여 유럽우주국(European Space Agency, ESA)에서 운용중인 Sentinel-2의 위성영상을 이용하였다(Table 1). 산불 발생 전・후 1개월 이내의 정규식생지수(Normalized Difference Vegetation Index, NDVI), 정규탄화지수(Normalized Burn Ratio, NBR), 수분스트레스지수(Moisture Stress Index, MSI)변화를 계산하였다(Table 2).
Table 1. Image data information.
Image source: Sentinel-2(https://apps.sentinel-hub.com/eo-browser/)
Table 2. Fomula of multispectral index.
NDVI는 식생의 활력을 정량화하고 모니터링하기 위해 가장 일반적으로 이용되는 다중 분광 지수 중 하나로 –1에서 1 사이의 값을 가지며 값이 클수록 건강한 식생을 나타낸다. NBR은 산불로 인한 피해를 정량화하기 위해 이용되는 지수로 식생에 민감하게 반응하는 NIR(Near InfraRed)과 토양 수분에 민감하게 반응하는 SWIR(Short-Wave InfraRed)를 이용하여 계산된다. 산불 전후 차이로 계산되는 dNDVI와 dNBR은 산불피해로 인한 나무의 주변 환경과 수관층의 변화를 간접적으로 나타낸다(Sparks et al., 2016; Furniss et al., 2020; Dindaroglu et al., 2021). MSI는 수관층의 수분함량에 민감하게 반응하는 지수로 값이 클수록 수관층이 큰 수분스트레스와 낮은 수분함량을 의미한다. 본 연구에서는 산불 전후의 각 지수들을 산출한 후 그 차이를 계산한 dNDVI(delta Normalized Difference Vegetation Index), dNBR(delta Normalized Burn Ratio), dMSI(delta Moisture Stress Index)을 산불피해 소나무의 생존분석의 인자로써 활용하였다. DEM자료에 기반한 지형 분석 및 Sentinel-2 위성영상을 이용한 다중분광지수의 계산은 오픈소스 프로그램인 QGIS (Quantum Geographical Information Software) 3.28.11를 이용하여 수행하였다. 조사 및 분석을 통해 얻어진 변수들은 각각 산불 피해목 인자, 지형 인자, 환경변화 인자로 분류하였다(Table 3).
Table 3. Factors and variables used for survival analysis.
3. 생존 분석
1) 콕스 비례위험모형
산불피해 소나무의 생존분석을 위하여 콕스 비례위험모형을 이용하였다. 생존분석에서는 생존시간을 이용하여 관찰 대상이 특정한 시간 t 이상 생존했을 때, t로부터 매우 단시간에 사건이 발생할 확률을 의미하는 위험 함수(hazard function)를 식 1과 같이 나타낼 수 있다.
\(\begin{align}h(t)=\lim _{\Delta t \rightarrow 0} \frac{P(t<T<t+\Delta t \mid T>t)}{\Delta t}\end{align}\) (1)
h(t) : hazard funtion
t : survival time
P : probability of occurrence of event
T : variable for survival time
콕스 비례위험모형에서는 위험함수 h(t)를 대상의 여러 특성(공변량, covariate)을 이용한 회귀식이 포함된 식 2와 같이 나타낸다.
h(t) = h0(t)e(β1x1 + β2x2 + ... + βkxk) (2)
h0(t) : basline hazard funtion at t
β : regression coefficient
x : covariate
모형의 식에서 h0(t)는 시간 t에서 모든 공변량이 0일 때, 사건에 대한 위험함수를 말하며 이를 기저위험함수라 부른다. 위험비(hazard ratio, HR)는 서로 공변량의 값이 다른 어떤 대상의 위험을 다른 대상의 위험으로 나눈 값이며 식 3과 같이 나타낼 수 있다.
\(\begin{align}\begin{aligned} H R & =\frac{h_{i}(t)}{h_{j}(t)}=\frac{h_{0}(t) e^{\left(\beta_{1} x_{i 1}+\ldots+\beta_{k} x_{i k}\right)}}{h_{0}(t) e^{\left(\beta_{1} x_{j 1}+\ldots+\beta_{k} x_{j k}\right)}} \\ & =e^{\left(\beta_{1}\left(x_{i 1}-x_{j 1}\right)+\ldots+\beta_{k}\left(x_{i k}-x_{j k}\right)\right)}\end{aligned}\end{align}\) (3)
HR: Hazard Ratio
결과적으로 위험비는 시간에 따라 변하는 기저위험함수가 사라지고 변수인 공변량과 회귀계수만 남게된다. 따라서 위험비를 통해 각 공변량이 사건의 발생에 미치는 영향을 분석하였다.
2) 비례위험 가정 검증
콕스 비례위험모형의 결과에 대한 신뢰성을 확보하기 위해 기본 가정인 “위험비는 시간에 따라 변하지 않으며 항상 일정하다.”라는 비례위험 가정(proportional hazard assumstion)에 대한 검증이 필요하다. 비례위험 가정의 검증에는 martingale 잔차 분석을 통해 모형에 사용된 변수(공변량)의 적합도를 평가하였다(Therneau et al., 1990).
관측치 중 이상치 확인을 위하여 삭제사례 편차(case deletion residuals) 분석을 실시하였다(Zhu et al., 2015). 삭제사례 편차는 모든 관측값을 이용하여 계산된 회귀계수와 특정 관측값을 제외하고 계산된 회귀계수의 차이를 계산하여 산출한 후 정규화를 실시한 dfbetas를 이용하였다. 삭제사례 편차 분석을 통해 모형의 결과에 큰 영향을 준 관측값을 확인하였다. 그 후 모든 변수들에 대하여 schoenfeld 잔차 그래프를 그리고 통계적인 분석을 통해 시간에 따라 위험비가 변하지 않는다는 비례위험 가정을 검증하였다(Grambsch and Therneau, 1994).
3) 생존곡선을 이용한 변수 영향력 분석
변수마다의 기준으로 여러 집단으로 나누어준 후 시간의 흐름에 따른 생존율 변화를 확인할 수 있는 생존곡선을 작성하였다(Table 4). 작성된 생존곡선을 분석한 결과를 기반으로 각 변수들이 산불피해 소나무의 생존율에 미치는 영향을 평가하였다.
Table 4. Group criteria for each variable for survival curves.
결과 및 고찰
1. 산불피해 소나무의 콕스 비례위험모형 위험비
산불피해 소나무의 콕스 비례위험모형의 분석 결과 각 변수별 위험비가 제시되었다(Figure 3). 산불피해목 인자의 경우, 흉고직경의 위험비는 1 미만인 0.93으로 나타나 소나무의 흉고직경이 클수록 고사 위험성이 감소하였다. 반면 BSH와 BSI의 위험비는 각각 1.01과 1.02로, 산불에 그을린 흔적이 클수록 고사 위험성이 증가하였다. 이러한 결과는 산불 이후 나무에 남겨진 그을음이 고사 발생과 상관성을 가진다는 선행 연구들과 일치하는 결과이다(Ordóñez et al., 2005; Regelbrugge and Conard, 1993; McHugh and Kolb, 2003; Furniss et al., 2018).
Figure 3. Hazard ratio of Individual variables.
C.I.: Confidence Interval.
**Correlation is significant at the 0.01 level(2-tailed).
***Correlation is significant at the 0.001 level(2-tailed).
지형 인자에서 경사의 위험비는 1.06으로 경사가 증가할수록 고사 위험성이 증가하는 것으로 분석되었다. 이러한 경향성은 경사가 낮을수록 산불 심각도가 낮아지고(Lentile et al., 2006), 경사가 높을수록 도복이 발생할 위험성이 커지기 때문으로 보여진다(Gibbons et al., 2008). 반면 고도의 경우 0.97의 위험비를 보여 높은 고도에 위치한 나무일수록 고사 위험성이 감소하는 것으로 나타났다. 그리고 사면방향은 위험비가 1을 나타내어 사면방향에 따른 고사위험성 차이는 명확하지 않았Furniss다.
환경변화 인자에서 dNDVI와 dNBR의 위험비는 각각 5.21과 4.05로 나타나 산불에 의한 식생 활력의 감소가 크고, 산불 심각도가 증가할수록 소나무의 고사 위험성이 증가하였다. 반면 dMSI의 위험비는 0.38로, 산불에 의한 수관층의 수분 스트레스 증가량이 클수록(dMSI의 값이 작을수록) 고사위험성이 증가하는 것으로 분석되었다. dNDVI와 dNBR은 산불피해로 인한 나무의 주변 환경과수관층의 변화를 나타내며(Sparks et al., 2016; Furniss et al., 2020; Dindaroglu et al., 2021), MSI는 수관층의 수분스트레스를 나타냄과 동시에 산불로 인한 교란의 정도와 높은 상관성을 지닌다(Avetisyan et al., 2022). 따라서 산불 전･후의 다중분광지수 변화가 산불피해 소나무의 활력 및 환경 변화를 나타내고 이를 통해 고사 발생을 예측하는 유의미한 변수로서 도출된 것으로 보인다.
2. 비례위험 가정 검증
1) martingale 잔차를 이용한 변수 적합도 평가
각각의 고사 영향 인자들의 변화에 대해 위험비가 일정하게 유지되는지를 검정하기 위해 martingale 잔차 그래프를 작성하였다(Figure 4). 흉고직경, 고도, 경사, dNDVI, dNBR, dMSI의 경우 martingale 잔차가 0을 중심으로 일정한 패턴을 보이지 않고 무작위로 분포되어 파란색의 회귀 곡선이 큰 기울어짐 없이 직선의 형태를 그렸다. 이는 해당 변수들의 값이 변화하더라도 위험비는 항상 일정하게 유지되어 해당 변수가 콕스 비례위험모형에 적합함을 의미한다. 범주형 변수인 사면방향은 방위별로 상자수염도를 그려 martingale 잔차의 분포를 확인하였다. 사방위 모두 잔차가 양이나 음의 값이 아닌 0을 중심으로 잔차가 분포되어 있어 사면방향이 고사예측에 적합한 변수임을 확인하였다. BSH와 BSI의 경우 각각 600 cm와 200 이하에서는 다른 변수들처럼 잔차가 0을 기준으로 특정한 패턴이 없이 무작위로 분포되어있어 콕스 비례위험모형에 적합한 것으로 나타났다. 그러나 그 이상에서는 martingale 잔차의 회귀곡선이 크게 음의 값으로 기울어지는 경향을 보였다. 이는 매우 높은 그을음 높이(BSH)와 면적(BSI)를 가지는 1개의 이상치가 원인임을 확인하였다.
Figure 4. The martingale residual plots of variables.
2) 삭제사례편차를 이용한 이상치 평가
BSH와 BSI의 변수적합도에 영향을 미친 이상치를 확인하기 위해 삭제사례편차(Case deletion residuals)의 한 종류인 dfbetas 그래프를 작성하였다(Figure 5). 데이터 중에서 절대값이 0.2 이상으로 비교적 큰 영향력을 나타내는 관측치는 BSH에서 7건(34, 215, 254, 264, 268, 280, 295), BSI에서는 6건(34, 215, 254, 264, 268, 295)을 확인하였다. martingale 잔차 그래프에서 BSH와 BSI에서 이상치로 지목된 관측치는 34번 소나무로 BSH가 765cm, BSI가 241.2로 산불로 인한 피해를 나타내는 그을음 특성이 다른 대상목들에 비해 매우 커 빠른 고사가 예상되는 나무임에도 불구하고 산불 발생 7년 차까지 생존하면서 이상치로 나타났다. 반면 나머지 이상치들은 큰 영향력을 나타내기는 하지만 그을음 자체가 다른 관측치들과 비슷하기 때문에 martingale 잔차 그래프에서는 이상치로 나타나지 않았다. 따라서 본 연구에서는 높은 그을음 높이와 면적을 가지는 대상목들이 적어 BSH와 BSI가 각각 600 cm와 200 이하일 때 고사예측에 적합한 변수임을 확인하였다.
Figure 5. Case deletion residuals plots of BSH(upper) and BSI(lower). The number in the boxes are the sample lables of the ouliers
3) 시간의 변화에 따른 비례위험 가정 검정
각각의 고사 영향 인자들의 위험비가 시간의 변화에도 항상 일정하게 유지되는지를 확인하기 위해 schoenfeld 잔차를 계산하여 시간에 대한 그래프로 나타내었다(Figure 6). schoenfeld 잔차로 그려진 회귀곡선(파란색)을 보면 위험비가 시간에 따라 항상 일정하지 않고 변화하는 것으로 보인다. 이는 관측된 산불피해 소나무의 고사가 연구기간 동안 항상 일정하게 발생한 것이 아니라 산불 발생 후 초기(1~3년 차)에 많은 고사가 발생하고 이후 시간이 지남에 따라 점차 고사 발생이 감소했기 때문이다. 반면 schoenfeld 잔차의 회귀직선(빨간색)을 보면 대부분의 변수가 0을 중심으로 작은 기울기를 보여 비례위험 가정이 만족되는 것으로 보인다. 이러한 회귀직선을 이용한 비례위험 가정의 검정이 통계적으로도 유의한지 확인하기 위해 잔차의 회귀직선과 시간에 따른 위험비의 변화가 없음을 의미하는 y=0의 직선에 대하여 카이제곱 검정(Chi-Squared Test)을 시행하였다. 결과적으로 모든 변수의 Schoenfeld 회귀직선이 y=0직선과 유의미한 차이가 없는 것으로 나타났다(p>0.05). 그러나 사면방향의 경우 회귀직선이 일정구간에서 회귀곡선의 95% 신뢰구간 범위(회색)에서 벗어나 있어 회귀직선이 회귀곡선을 대변한다고 할 수 없다. 따라서 사면방향을 제외한 모든 변수가 비레위험 가정을 만족시키는 것으로 나타났다.
Figure 6. schoenfeld residual plot of variables. blue curve: regression curve of schoenfeldresiduals. red line: regression line of schoenfeldresiduals. Schoenfeld Test p = Results of chi-square test of residual regression straight and y=0 straight lines p-values.
3. 변수별 산불피해 소나무의 생존곡선
1) 산불피해 소나무 고사 모니터링 결과
산불피해 소나무의 생존곡선을 작성한 결과(Figure 7), 산불 발생 5년 차까지 지속적으로 고사가 발생했으며 그중에서도 3년 차에 가장 큰 생존율의 감소가 확인되었다. 반면 산불 발생 6년 차부터는 거의 고사가 발생하지 않아 7년차에는 생존율이 80.29%로 안정되는 결과를 나타내었다.
Figure 7. Survival curve of forest fire damaged Pinus densiflora.
2) 산불피해목 인자(DBH, BSH, BSI)
산불피해목 인자가 소나무의 고사에 미치는 영향을 확인하기 위해 변수마다 집단을 나누어 생존곡선을 그려 분석하였다(Figure 8). 흉고직경의 경우 대경목(Large, DBH ≥ 30cm)의 생존곡선이 소경목(Small, DBH < 30cm)의 생존곡선보다 7년동안 일관되게 높은 생존율을 유지하였다. 결과적으로 대경목이 소경목보다 13.59% 더 높은 생존율을 보여 큰 흉고직경을 가지는 소나무 임분의 경우 산불 이후 생존 가능성이 높은 것으로 나타났다.
Figure 8. Survival curves by variable groups.
BSH도 생존곡선들이 교차하지 않아 그을음 높이에 따른 생존율을 차이를 확인할 수 있었다. 낮음(Low, BSH < 150cm)의 경우 고사가 거의 발생하지 않아 96.28%의 높은 생존율을 나타내는 반면, 중간(Middle, 150cm ≤ BSH < 300cm)과 높음(High, BSH ≥ 300cm)은 생존율이 각각 64.41%와 50.00%로 나타나 수간에 남겨진 그을음의 높이가 높은 소나무는 낮은 소나무에 비해 고사 확률이 높았다.
BSI의 경우 이러한 경향이 BSH에 비해 매우 명확하게 나타났다. 낮음(Low, BSI < 150)의 경우 96.26%의 높은 생존율을 보이는 반면 중간(Middle, 150 ≤ BSI < 300)은 59.60%, 높음(High, BSI ≥ 300)은 33.33%로 현저하게 낮은 생존율을 보였다.
BSI에 따른 생존율 차이가 BSH에 따른 생존율 차이보다 더 크게 나타낸 것은 BSH가 그을음의 평균적인 높이를 나타내는데 비해 BSI는 그을음의 높이 대비 비율을 합으로 계산되어 그을음의 전체적인 면적을 반영하기 때문이다. 따라서 BSI가 BSH보다 더 명확한 고사 예측 인자로 보여진다. 그러나 높이 대비 그을음 비율의 측정은 조사자에 따라 달라질 수 있으며 높이만 측정하는 것에 비해 많은 시간이 소요되기 때문에 대면적의 조사가 필요한 경우에는 간단한 BSH의 측정이 유용할 것으로 생각된다.
3) 지형 인자(고도, 경사, 사면방향)
고도의 경우 고고도(High, Elevation ≥ 800m)에 위치한 소나무들이 저고도(Low, Elevation < 800m) 대비 27.05% 더 높은 생존율을 보여, 높은 고도에 위치한 산불피해 소나무일수록 높은 생존 가능성을 가졌다. 그러나 이러한 차이는 비교적 저고도에 위치한 연구 대상지1과 2에 비해 고고도에 위치한 연구대상지 3의 소나무는 비교적 산불피해(그을음 등)가 약하고 큰 흉고직경을 가지기 때문에 온전히 고도에 의한 차이로 보기에는 한계가 있는 것으로 보인다.
경사의 경우 급경사(Steep, Slope ≥ 20°)가 완경사(Gentle, Slope < 20°)에 비해 15.32% 더 낮은 생존율을 나타내었다. 이는 연구 초기에 생존이 확인되었던 소나무가 다음 해에 도복으로 인해 고사한 사례가 드물게 관측되었기 때문으로, 산불 이후 경사가 클수록 도복 위험성이 커진다는 Gallaway et al.(2009)의 연구 결과와 일치한다. 그러나 급경사와 완경사 집단 둘 모두가 산불 발생 5년 차까지 꾸준히 고사가 발생하고 비교적 작은 생존율 차이를 나타내는 것으로 보아 경사가 산불피해 소나무의 고사에 미치는 영향은 산불피해목 인자에 비해 작은 것으로 보인다.
사면방향의 경우 각 방위별로 생존율의 차이가 발생하기는 하나 동사면에 위치한 소나무가 다른 사면에 비해 매우 적어 실질적으로 사면방향에 따른 생존율 차이를 확인하기에는 한계가 있었다(Figure 9). 또한 북사면(N)과 서사면(W)의 생존곡선이 교차하는 것으로 나타나 비례위험 가정을 만족하지 못함을 다시 한번 더 확인할 수 있었다.
Figure. 9 Histogram of the number of forest fire damaged Pinus densiflora by aspect.
4) 환경변화 인자(dNDVI, dNBR, dMSI)
변수별 산불피해 소나무 집단의 생존 곡선 dNDVI에서 산불 이후 식생 활력의 변화가 없음을 나타내는 변화 없음(Unaltered, dNDVI < 0.00)의 경우 모든 소나무가 연구 기간 동안 생존한 반면, 조금 감소(Slight decrease, 0.00 ≤ dNDVI < 0.05)와 심각하게 감소(Severe decrease, dNDVI ≥ 0.05)의 경우 각각 88.68%와 68.59%의 생존율을 보이며 산불 이후 식생활력이 감소할수록 고사위험성이 커졌다. 심각하게 감소한 집단은 연구 기간 동안 꾸준히 고사가 발생하는 경향을 보인 반면 조금 감소한 집단은 산불 발생 1년차에는 고사가 발생하지 않고 2~4년차에만 고사가 발생한 후 생존이 안정되는 차이를 나타내었다. 이러한 dNDVI에 따른 생존율 차이는 수관층에 직접적으로 영향을 미치지 않는 지표화 피해를 입더라도 산불피해가 수목의 전체에 영향을 미치며, 더 나아가 고사 발생을 일으킬 수 있음을 시사한다.
산불 심각도를 나타내는 dNBR은 피해 없음(Unbured, dNBR < 0.100)에서 고사가 발생하지 않아 dNDVI와 비슷한 경향을 보였다. 반면 매우 낮은 심각도(Very low severity, 0.100 ≤ dNBR < 0.185)와 낮은 심각도(Low severity, 0.185 ≤ dNBR < 0.270)는 각각 생존율이 79.49%, 71.43%로 둘 모두 산불 발생 5년 차까지 지속적으로 고사가 발생하였다. 그룹을 나누는 기준값인 0.100과 0.270은 미국지질조사국(Unuted States Geological Survey, USGS)에 의해 제시된 수치이며 0.185는 그 중간값으로 추가한 값이다. 이렇듯 산불 피해도 “경”에 해당하는 Low severity라는 1개의 등급 내에서 dNBR을 두 집단으로나눴음에도 생존율의 차이가 발생한 것으로 보아 산불피해 소나무의 처리 기준 마련을 위해서는 기존의 피해 등급 기준보다는 새로운 기준의 제시가 필요할 것으로 생각된다.
dMSI 역시 수분 스트레스 증가가 없는 안정적인 집단(Stable, dMSI < -0.15)은 고사가 발생하지 않았다. 중간정도의 스트레스(Moderate stess, -0.15 ≤ dMSI < 0.00)의 경우 연구기간 동안 지속적으로 고사가 발생하긴 했으나 산불 발생 3년차에 가장 많은 고사가 발생하여 7년 차에는 78.11%의 생존율을 나타낸 반면, 심각한 스트레스(Severe stress, dMSI ≥ 0.00)는 산불 발생 5년 차까지 지속적인 고사가 발생한 점은 동일하나 산불 발생 1년 차와 4년 차에 고사가 집중적으로 발생하는 차이를 보였으며 7년 차에는 60.27%의 생존율을 보였다. 이러한 결과는 산불피해가 수목의 관다발계에 공동화(cavitation)를 일으켜 수분 스트레스가 증가하고 고사까지 발생시킨다는 Kavanagh et al. (2010)의 연구결과와 일치한다.
환경변화 인자로 이용된 dNDVI와 dNBR, dMSI는 위성 영상 기반의 영향인자로 대면적의 산불피해지나 사람이 접근하기 어려운 지역의 산불피해 소나무의 처리에 대하여 의사결정이 필요한 경우 활용가능할 것으로 보인다.
결론
본 연구는 산불피해 소나무의 고사 모니터링을 바탕으로 콕스 비례위험모형을 이용하여 고사 발생에 미치는 영향 인자를 도출하고 생존율 비교를 통해 인자별 영향력을 분석하고자 하였다. 사면방향을 제외한 산불피해목 인자(DBH, BSH, BSI), 지형인자(고도, 경사), 환경변화 인자(dNDVI, dNBR, dMSI)들이 비례위험 가정을 만족하여 산불피해 소나무의 고사에 영향을 미치는 것으로 나타났다. 단 BSH와 BSI의 경우 매우 높은 그을음 높이와 면적을 가지는 대상목들이 적어 본 연구의 결과를 산불 피해도가 큰 지표화 피해지에 적용하기에는 한계가 있을 것으로 보인다. 생존율 비교를 통해 인자별 영향력을 분석한 결과, 인자들 중 나무의 수피에 남겨진 그을음과 관련된 변수인 BSH, BSI와 산불 전후 수관층의 수분스트레스 변화를 나타내는 dMSI가 산불피해 소나무의 생존 예측에 큰 영향력을 가지는 것으로 나타났다. 현장 조사를 통해 얻어지는 BSH와 BSI는 고사예측을 통한 피해 산정이 필요한 경우에 활용될수 있을 것으로 보이며, 원격탐사로 얻어진 dMSI는 기존의 NDVI와 dNBR을 이용한 대면적의 산불 피해 산정에 추가적인 분석 자료로써 이용할 수 있을 것으로 보여진다. 본 연구는 산불피해 소나무 고사 발생에 영향을 미치는 인자를 규명하던 이전의 연구에서 더 나아가 그 영향력을 평가함으로써, 산불복원 계획 수립에 필요한 산불피해 소나무의 처리 의사결정 및 산불피해지 조사인자 선택에 유용할 것으로 기대된다. 또한 향후 기후인자를 고려하는 등의 다양한 산불피해목 생존 진단 연구 확장에 기초자료로 활용될 것이다.
감사의 글
본 연구는 국립산림과학원의 ‘산불피해지 복원 프로세스 및 내화숲 기능증진 연구’(과제번호 FE0100-2022-02-2024)의 지원으로 수행되었습니다.
References
- Avetisyan, D., Velizarova, E. and Filchev, L. 2022. Post-Fire Forest Vegetation State Monitoring through Satellite Remote Sensing and In Situ Data. Remote Sensing 14(24): 6266.
- Bowman, D.M. et al. 2009. Fire in the Earth system. Science 324(5926): 481-484.
- Brando, P.M., Nepstad, D.C., Balch, J.K., Bolker, B., Christman, M.C., Coe, M. and Putz, F.E. 2012. Fire-induced tree mortality in a neotropical forest: the roles of bark traits, tree size, wood density and fire behavior. Global Change Biology 18(2): 630-641.
- Chandler, C., Cheney, P., Thomas, P., Trabaud, L. and Williams, D. 1983. Fire in forestry. Volume 1. Forest fire behavior and effects. New York: Wiley.
- Cox, D.R. 1972. Regression Models and Life-Tables. Journal of the Royal Statistical Society. Series B 34 (2): 187-220.
- Das, A.J., Stephenson, N.L. and Davis, K.P. 2016. Why do trees die? Characterizing the drivers of background tree mortality. Ecology 97(10): 2616-2627.
- Dickinson, M. and Johnson, E. 2001. Fire effects on trees Forest Fires: Behavior and Ecological Effects. E A Johnson and K Miyanishi. New York: Academic
- Dindaroglu, T., Babur, E., Yakupoglu, T., Rodrigo-Comino, J. and Cerda, A. 2021. Evaluation of geomorphometric characteristics and soil properties after a wildfire using Sentinel-2 MSI imagery for future fire-safe forest. Fire Safety Journal 122: 103318.
- Furniss, T.J., Kane, V.R., Larson, A.J. and Lutz, J.A. 2018. Multi-scale assessment of post-fire tree mortality models. International Journal of Wildland Fire 28(1): 46-61.
- Furniss, T.J., Kane, V.R., Larson, A.J. and Lutz, J.A. 2020. Detecting tree mortality with Landsat-derived spectral indices: Improving ecological accuracy by examining uncertainty. Remote Sensing of Environment 237: 111497.
- Gallaway, J.M., Martin, Y.E. and Johnson, E.A. Johnson. 2009. Sediment transport due to tree root throw: integrating tree population dynamics, wildfire and geomorphic response. Earth Surface Processes and Landforms 34(9): 1255-1269.
- Gibbons, P., Cunningham, R.B. and Lindenmayer, D.B. 2008. What factors influence the collapse of trees retained on logged sites?: a case-control study. Forest Ecology and Management 255(1): 62-67.
- Grambsch, P.M. and Therneau., T.M. 1994. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika 81(3): 515-526.
- Jung, J.K., Kim, M., Nam, Y. and Koh, S.H. 2020. Changes in spatial and temporal distributions of Monochamus beetles along the fire severity in burned Pinus densiflora forests. Journal of Asia-Pacific Entomology 23(2): 404-410.
- Kane, J.M., Varner, J.M., Metz, M.R. and van Mantgem, P.J. 2017. Characterizing interactions between fire and other disturbances and their impacts on tree mortality in western US Forests. Forest Ecology and Management 405: 188-99
- Kavanagh, K.L., Dickinson, M.B. and Bova, A.S. 2010. A way forward for fire-caused tree mortality prediction: modeling a physiological consequence of fire. Fire Ecology 6: 80-94.
- Korea Forest Service. 2019. Forest Fire Statistics Annual Report for 2018. Korea Forest Service.
- Korea Meteorolgical Administration. Korea's 30-year average climate data. 2024. https://data.kma.go.kr/climate/average30Years/selectAverage30YearsList.do?pgmNo=113. (2024.02.28.)
- Kwon, S.M., Kim, S.H., Kim, J.H., Kang, W.S., Park, K.H., Kim, C.B., and Girona, M.M. 2021. Predicting post-fire tree mortality in a temperate pine forest, Korea. Sustainability 13(2): 569.
- Lawes, M.J., Richards, A., Dathe, J., and Midgley, J.J. 2011. Bark thickness determines fire resistance of selected tree species from fire-prone tropical savanna in north Australia. Plant Ecology 212: 2057-2069.
- Lee, H.N., Yun, K.H. and Kim, G.H. 2023. Mapping burned forests using a k-Nearest neighbors classifier in complex land cover. Korean Society of Civil Engineers 43(5): 883-896.
- Lee, H.J., Lee, J.M., Won, M.S. and Lee, S.W. 2012. Development and Verification of Korean Composit Burn Index (KCBI). Journal of Korean Forest Society 101(1): 163-170.
- Lee, S.M. and Jung, J.C. 2019. Classification of forest fire damage intensity using probability density functions and KOMPSAT-3A. Korean Journal of Remote Sensing 35(6-4): 1341-1350.
- Lee, S.Y. and Lee, H.P. 2006. Analysis of forest fire occurrence in Korea. Journal of Korean Institute of Fire Science & Engineering 20(2): 54-63.
- Lee, S.Y., Jeon, G.W., Lee, M.W. and Jeon, G.W. 2008. Mortality in Pine Stand and Vegetation Recovery after Forest Fire. Journal of the Korean Society of Hazard Mitigation 8(1): 71-79.
- Lee, Y.G., Gu, K.S., Jo, B.S., Kang, Y.H., Lee, M.B., Jeon, B.K., Im, J.H., Lee, K.S. and Lee, Y.S. 2006. Manual for the Restoration of Fire-Damaged Sites for Sustainable Forest Recovery and Ecological Restoration. National Institute of Forest Science.
- Lentile, L.B., Smith, F.W. and Shepperd, W.D. 2006. Influence of topography and forest structure on patterns of mixed severity fire in ponderosa pine forests of the South Dakota Black Hills, USA. International Journal of Wildland Fire 15(4): 557-566.
- Maringer, J., Ascoli, D., Kuffer, N., Schmidtlein, S. and Conedera, M. 2016. What drives European beech (Fagus sylvatica L.) mortality after forest fires of varying severity?. Forest Ecology and Management 368: 81-93.
- McHugh, C.W. and Kolb, T.E. 2003. Ponderosa pine mortality following fire in northern Arizona. International Journal of Wildland Fire 12(1): 7-22.
- Michaletz, S.T. and Johnson, E.A. 2007. How forest fires kill trees: a review of the fundamental biophysical processes. Scandinavian Journal of Forest Research 22(6): 500-515.
- Ordonez, J.L., Retana, J. and Espelta, J.M. 2005. Effects of tree size, crown damage, and tree location on post-fire survival and cone production of Pinus nigra trees. Forest Ecology and Management 206(1-3): 109-117.
- Pascual, A., Tupinamba-Simoes, F., Guerra-Hernandez, J. and Bravo, F. 2022. High-resolution planet satellite imagery and multi-temporal surveys to predict risk of tree mortality in tropical eucalypt forestry. Journal of Environmental Management 310: 114804.
- Rao, K., Anderegg, W.R., Sala, A., Martinez-Vilalta, J. and Konings, A.G. 2019. Satellite-based vegetation optical depth as an indicator of drought-driven tree mortality. Remote Sensing of Environment 227: 125-136.
- Regelbrugge, J.C. and Conard S.G. 1993. Modeling tree mortality following wildfire in Pinus ponderosa forests in the central Sierra-Nevada of California. International Journal of Wildland Fire 3(3): 139-148.
- Shim, Sang Taek. 2018. Empirical Analysis of the Proliferation Factors of Forest Disaster Using Spatial Big Data: Focusing on Pine Wilt Disease. Kongju. Kongju National University.
- Sparks, A.M., Kolden, C.A., Talhelm, A.F., Smith, A.M., Apostol, K.G., Johnson, D.M. and Boschetti, L. 2016. Spectral indices accurately quantify changes in seedling physiology following fire: towards mechanistic assessments of post-fire carbon cycling. Remote Sensing 8(7): 572.
- Therneau, T.M., Grambsch, P.M. and Fleming, T.R. 1990) Martingale-based residuals for survival models. Biometrika 77(1): 147-160.
- Van Mantgem, P.J., Stephenson, N.L., Knapp, E., Battles, J. and Keeley, J.E. 2011. Long-term effects of prescribed fire on mixed conifer forest structure in the Sierra Nevada, California. Forest Ecology and Management 261: 989-94.
- Won, M.S., Kim, K.H. and Lee, S.W. 2013. Development of Quantitative Evaluation Techniques for Forest Fire Damage Intensity and Identification of Damage Characteristics. National Institute of Forest Science.
- Whelan, R.J. 1995. The ecology of fire. Cambridge university press.
- Zhu, H., Ibrahim, J.G. and Chen, M.H. 2015. Diagnostic measures for the Cox regression model with missing covariates. Biometrika 102(4): 907-923.