Predicting Site Quality by Partial Least Squares Regression Using Site and Soil Attributes in Quercus mongolica Stands

신갈나무 임분의 입지 및 토양 속성을 이용한 부분최소제곱 회귀의 지위추정 모형

Abstract

Predicting forest productivity is essential to evaluate sustainable forest management or to enhance forest ecosystem services. Ordinary least squares (OLS) and partial least squares (PLS) regression models were used to develop predictive models for forest productivity (site index) from the site characteristics and soil profile, along with soil physical and chemical properties, of 112 Quercus mongolica stands. The adjusted coefficients of determination (adjusted R²) in the regression models were higher for the site characteristics and soil profile of B horizon (R²=0.32) and of A horizon (R²=0.29) than for the soil physical and chemical properties of B horizon (R²=0.21) and A horizon (R²=0.09). The PLS models (R²=0.20-0.32) were better predictors of site index than the OLS models (R²=0.09-0.31). These results suggest that the regression models for Q. mongolica can be applied to predict the forest productivity, but new variables may need to be developed to enhance the explanatory power of regression models.

산림생산력의 예측은 지속가능한 산림경영이나 산림생태계서비스 증진을 위한 산림관리에 필수적인 것으로 알려져 있다. 본 연구는 전국 112개 신갈나무 임분을 대상으로 입지 및 토양 특성의 속성변수를 이용하여 지위 추정 모형을 개발하였다. 신갈나무 임분의 지위지수는 입지 및 토양 특성을 독립변수 한 일반최소제곱(Ordinary Least Squares, OLS) 및 부분최소제곱(Partial Least Squares, PLS) 회귀모형을 이용하여 유의적인 추정(P<0.05)이 가능하였다. 지위추정 회귀모형의 수정 결정계수(adjusted R²)는 입지 및 토양단면 속성변수의 회귀모형(A층: R²=0.29; B: R²=0.32)이, 토양 물리·화학적특성의 속성변수(A층: R²=0.09; B층: R²=0.21)보다 높게 나타났다. 한편, PLS 회귀모형(R²=0.20~0.32)은 OLS 회귀모형(R²=0.09~0.31)에 비해 지위지수 추정식의 설명력이 높았다. 본 연구로부터 신갈나무 임분의 입지 및 토양 특성을 이용한 지위 추정 회귀 모형이 개발되었으나, 결정계수 값이 낮아 회귀모형의 설명력을 향상시킬 수 있는 새로운 변수 개발이 필요할 것으로 사료되었다.

서론

지위지수(Site Index)는 산림 생태계의 지위(site quality)를 판정하기 위해 일반적으로 사용하는 측정 방법으로 강수량, 온도, 임목 생육기간, 임분 밀도, 경쟁식생, 시비(Hagglund, 1981) 등에 의해 영향을 받기도 하지만 토양의 물리·화학적 특성(Lee et al., 1996; Woolery et al., 2002)에 의해 결정되기도 한다. 지위지수에 의한 산림 생산능력 파악은 산림경영적인 관점에서 중요하며 임분 내 우세목과 준우세목의 수고를 이용하여 지위지수 곡선으로부터 지위를 판정하는 방법이 일반적으로 사용되고 있다(Hagglund, 1981; Sturtevant and Seagle, 2004; Lee et al., 2007; Subedi and Fox, 2016). 그러나 공간적으로 넓은 범위에 분포하는 수종의 경우 지위 추정에 적합한 임분이나 임목을 찾는 것이 쉽지 않고, 또한 천연림은 다양한 영급과 수종의 혼효로 인하여 지위지수 곡선으로부터 지위 판정이 상당히 어려운 실정이다. 그러한 경우 토양, 지형, 기후 등과 같은 입지환경 및 토양 특성과 지위지수와의 관계를 조사하기도 한다(Woolery et al., 2002; Hamel et al., 2004; Pietrzykowski et al., 2015).

국내에서 참나무류의 입지환경 및 토양 특성을 이용한 산림생산력 추정관련 연구로 굴참나무(Lee et al., 1996)와 상수리나무(Jung et al., 2017) 임분의 토양 물리·화학적 특성을 이용한 지위 추정 모형, 1:25,000 산림입지도의 입지환경 및 토양단면의 속성 정보를 이용한 상수리나무나 신갈나무(Lee et al., 2007), 상수리나무(Kim et al., 2013), 신갈나무(Lee et al., 2014) 등의 적지 판정 및 지위 추정 연구가 실시된 바 있다. 그러나 이들 연구의 대부분은 지위추정에 일반최소제곱(Ordinary Least Squares, OLS) 회귀모형을 제시하고 있으나, 토양 특성을 결정하는 변수들 사이에는 강한 상관관계가 존재하여 다중공선성(multicollinearity) 문제가 발생할 수 있다(Subedi and Fox, 2016). 이러한 문제점을 제거하기 위해서는 독립변수 상호 간 분산팽창인자(Variance Inflation Factor, VIF) 값으로부터 다중공선성이 있는 독립변수를 제거하고 회귀분석을 실시하기도 한다(Subedi and Fox, 2016). 한편, 부분최소제곱(Partial Least Squares, PLS) 회귀는 각 독립변수와 종속변수 간의 다중공선성 문제를 제거할 수 있는 통계 모형으로 알려져 있다(Carrascal et al., 2009). 산림 분야에서 수행된 PLS 연구로는 미국 테다소나무 임분의 토양 물리·화학적 특성을 이용한 지위 추정모형이 개발된 바 있으며(Subedi and Fox, 2016), 국내에서는 산불 발생 요인에 대한 기상 및 수문학적 요인 분석(Kim et al., 2021)이나, 미래 육상식생의 생산성 분석 예측(Choi et al., 2017) 등에 PLS를 적용한 연구가 보고된 바 있다.

국내에서 개발된 참나무류의 토양 물리·화학적 특성을 이용한 지위추정 모형은 굴참나무 임분의 단계적 회귀모형(Lee et al., 1996), 상수리나무 임분의 단계적 회귀모형과 주성분 회귀(Principle Component Regressions)를 이용한 연구가 수행된 바 있다(Jung et al., 2017). 그러나 이전에 수행된 굴참나무나 상수리나무 등의 연구는 변수 선발에 다중공선성을 고려하지 않거나, 회귀식에 10개 이상의 설명변수가 포함되어 자료수집에 많은 노력이 요구되는 문제점이 있었다.

우리나라 온대 낙엽활엽수류 중 대표적 우점종인 참나무류는 분포 면적이 약 980천ha 정도이며(Park et al., 2020), 이 중 신갈나무의 분포 면적은 414천ha로 가장 넓은 면적을 차지하고 있다(Lee et al., 2018). 주로 순림의 상태로 존재하기 때문에 집약적인 산림경영을 위한 수종으로 고려될 수 있으며, 탄소저장능력도 높은 것으로 보고되고 있다(Lee et al., 2022). 한편, 신갈나무는 지리적으로 광범위한 지역에 분포하기 때문에 공간적으로 넓은 범위에서 조사된 입지 및 토양 환경 정보를 이용한 지위판정이 요구된다. 따라서 본 연구는 전국에 분포하는 신갈나무 임분을 대상으로 산림생산력 판정이나 적지선정을 위한 기초자료 제공을 목적으로 입지환경 및 토양단면의 속성변수와 토양 물리·화학적 특성을 이용하여 지위지수를 추정하는 OLS나 PLS 회귀모형을 제시하고자 수행하였다.

재료 및 방법

1. 조사지 선정 및 조사방법

본 연구를 위한 조사는 2016년에서 2021년 사이 전국 경제림 육성단지를 중심으로 기존 문헌 분석 및 임상도(1:5,000)를 이용하여 신갈나무 임분의 분포 현황을 파악하였고, 그 중 신갈나무의 점유비율이 75% 이상인 임분을 조사구로 선정하였다. 조사한 표본수는 강원도가 68개소로 가장 많았으며, 충청도 19개소, 경상도 9개소, 경기도 9개소, 전라도 8개소 등 총 112개소였다.

조사구는 반경 11.3 m의 원형 표본 조사구(400 m²)를 사면경사와 방위 등을 고려하여 설치하고, 지황(지형, 경사, 방위 등)과 임황(수종, 흉고직경, 수고 등)을 조사하였다. 각 조사구 내 임목들은 직경테이프를 이용하여 흉고직경(1.2 m)을 측정하였으며, 임분의 지위 추정을 위해 각 조사구 당 4∼5본의 우세목 및 준우세목의 수고를 height meter (Suunto PM-5/1520)로 측정하였다. 또한, 각 조사구로부터 대표적인 표본목을 1본 선정하고 생장추를 이용하여 흉고 부위에 목편(increment core)을 채취한 후 임분 연령을 추정하였다. 임목 축적은 산림청과 국립산림과학원에서 개발한 임목수간재적표(Korea Forest Service and National Institute of Forest Science, 2021)를 이용하여 계산하였으며, 지위지수는 국가산림자원조사자료를 바탕으로 개발된 참나무류(Park et al., 2020)의 지위지수 추정식을 이용하여 계산하였다. 조사한 신갈나무 임분의 평균 임령은 62년, 임분밀도는 742본 ha^-1, 흉고단면적은 30.2 m² ha^-1, 임목축적은 213 m³ ha^-1였다(Table 1).

Table 1. General stand characteristics in Quercus mongolica stands.

2. 토양 특성 분석

토양 특성 분석을 위해 가로 80 cm, 세로 80 cm에 모재층 깊이까지 토양단면을 제작하고 A층과 B층에 대한 토양 단면조사를 실시한 후 토양 물리·화학적 특성 분석용 시료를 채취하였다. 채취된 토양 시료는 실내에서 2주 이상 음건 후 2 mm 체로 선별하여 분석용 시료를 조제하였다. 토양 유기물 분석은 Walkley-Black 방법(Kalra and Maynard, 1991), 유효인산은 Lancaster침출법, 교환성양이온(K⁺, Ca²⁺, Mg²⁺, Na⁺)은 1N-NH₄OAc (pH 7.0)로 침출 후 ICP를 이용하여 분석하였다. 양이온교환용량(Cation Exchange Capacity)은 1N-NH₄OAc (pH 7.0)로 침출 후 NH₄⁺-N 증류법, 토양 pH는 토양:증류수(1:5) 현탁액을 pH 메터로 측정하였다. 전질소(Total nitrogen)는 Kjeldahl 황산분해증류법, 모래, 미사, 점토 함량은 피펫법을 이용하였다(Kalra and Maynard, 1991).

3. 자료분석

지위추정 회귀모형의 종속변수는 지위지수이며, 독립변수는 입지 및 토양 단면 조사로부터 얻어진 7개 속성변수와 A층과 B층의 토양단면에 대한 물리·화학적 특성 분석으로부터 얻어진 13개 속성변수였다. 입지 및 토양단면 속성변수는 해발고, 색상, 명도, 채도, 경사, 방위, 토심 등이며, 토양 물리·화학적 특성의 속성변수는 모래, 미사, 점토, 토양 pH, 토양유기물, 전질소, C/N비(탄질률), 유효인산, 교환성 칼슘(Ca²⁺), 교환성 마그네슘(Mg²⁺), 교환성 포타슘(K⁺), 교환성 소듐(Na⁺)과 양이온교환용량(CEC) 등이다. OLS 회귀 모형의 변수 선발은 다중공선성 문제점을 해결하기 위해 독립변수 상호 간 분산팽창인자(Variance Inflation Factor, VIF) 값을 점검하고 일반적으로 사용하는 변수선택 기준 값인 VIF 10 이상 변수를 제외한 후(O’brien, 2007), SAS의 변수증감법(Stepwise Multiple Linear Regression, SMLR)으로 유의 수준 P<0.15에서 최적 변수들이 선택될 수 있도록 하였다(Statistical Analysis System Inc., 2003). PLS 회귀모형은 변수중요척도(Variable Importance in the Projection, VIP)를 XLSTAT에 의해 계산하고 VIP 0.8 이상의 변수를 회귀모형에 포함하였다(Subedi and Fox, 2016). 회귀모형의 설명력은 수정결정계수(adjusted R²)로 비교하였다.

결과

1. 입지 및 토양단면과 토양 물리·화학적 특성

신갈나무 임분의 평균 지위지수는 15.5였으며 해발고는 166 m에서 1,268 m까지 분포하여 수직적 분포범위가 크게 나타났다(Table 2). 평균토심은 A층 16 cm, B층 40 cm였으며, 토색 중 채도와 명도는 B층이 A층에 비해 높게 나타났다. 토양 물리·화학적 특성 중 모래, 미사, 점토 함량은 두 층위 간 큰 차이는 없었으나, 토양 pH는 A층이 pH 4.2로 B층의 pH 4.5에 비해 산성화되어 있었다(Table 3). 유효인산과 포타슘, 칼슘, 마그네슘 같은 교환성 양이온은 A층이 B층에 비해 높은 값을 보였다.

Table 2. Selected variables of site or soil profile property in Quercus mongolica stands.

Table 3. Selected variables of soil physical and chemical properties in Quercus mongolica stands.

Note: OM: organic matter content; TN: total nitrogen; CEC: cation exchange capacity

2. 입지 및 토양 속성 사이 상관관계

입지환경 속성인 해발고(r=0.24; r=0.34) 및 경사(r=0.20; r=0.18)는 지위지수와 정의 상관, 토양단면 속성인 색상(r=-0.27; r=-0.28)은 지위지수와 부의 상관이 있었다. 입지 및 토양단면 속성 중 해발고는 토색이나 토심과 부의 상관, 색상, 채도, 명도 등은 변수 사이에 강한 정의 상관을 보였다(Table 4). A층의 경우 물리·화학적 특성의 속성변수와 지위는 상관관계가 나타나지 않았으나, B층의 전질소와 지위는 정의 상관(r=0.17), 토양 pH(r=-0.25) 및 교환성 K⁺(r=-0.24)와 지위는 부의 상관이 있었다(Table 5). 한편, 토양 A층과 B층 모두 토양 pH, 유기물, 전질소, 탄질률 등의 변수 사이에 강한 상관관계를 보였다.

Table 4. Pearson’s correlation of site index and site or soil profile properties of the soil horizons in Quercus mongolica stands.

Note: Red or blue color indicates negative or positive values, respectively. Intensity of the color represents the magnitude of the values. The bold values indicate a significance between two variables at P<0.05.

Table 5. Pearson’s correlation of site index and soil physical and chemical properties of the soil horizons in Quercus mongolica stands.

Note: Red or blue color indicates negative or positive values, respectively. Intensity of the color represents the magnitude of the values. The bold values indicate a significance between two variables at P<0.05.

3. 변수 중요척도

입지 및 토양단면 속성변수로부터 얻어진 VIP에 따르면 해발고(A층: 1.3, B층 1.8), 색상(A층: 1.2, B층 1.5), 채도(A층: 0.9, B층 1.2) 등이 PLS 회귀에 포함될 수 있는 중요 변수였다(Figure 1). 토양 물리·화학적 특성은 A층과 B층의 경우 유효인산(A층: 1.6, B층: 1.2), 토양 pH(A층: 1.4, B층: 1.4), 전질소(A층: 1.1, B층: 1.3), 탄질률(A층: 1.2, B층: 0.9), 교환성 마그네슘(A층: 0.9, B층: 1.3) 등이 회귀모형에 포함될 변수로 선발되었다(Figure 1). 입지 및 토양단면과 토양 물리·화학적 특성을 모두 포함하는 20개의 속성변수는 색상(A층: 2.1, B층: 1.8), 해발고(A층: 1.9, B층: 2.1), 경사(A층: 1.6, B층: 1.2), pH(A층: 1.3, B층: 1.6) 등이 회귀모형에 포함될 변수였다(Figure 1).

Figure 1. Variable Importance in the Projection (VIP) values for the candidate soil variables in the Partial Least Squares (PLS) regression models for predicting site index of Quercus mongolica stands.

4. 일반최소제곱과 부분최소제곱 회귀에 의한 지위추정 회귀 모형

지위추정을 위한 OLS 회귀모형의 수정결정계수(adjusted R²)는 입지 및 토양단면 속성변수 회귀모형의 경우 R²=0.29와 R²=0.31이었으나, 토양 물리·화학적 특성의 속성변수는 R²=0.09와 R²=0.21로 회귀모형의 설명력은 입지 및 토양단면 속성변수가 높게 나타났다(Table 6). 유사한 결과로 PLS 회귀모형의 경우도 입지 및 토양단면 속성 변수는 R²=0.25와 R²=0.32로 토양 물리·화학적 특성 속성 변수 R²=0.05와 R²=0.20보다 설명력이 높았다. 입지 및 토양단면과 토양 물리·화학적 특성의 총 20개 속성변수를 이용한 OLS나 PLS 회귀모형은 입지 및 토양단면의 7개 속성변수나 토양 물리·화학적 특성의 13개 속성변수로 부터 개발된 각각의 회귀모형 결정계수 값과 큰 차이가 없었다.

Table 6. Ordinary Least Squares (OLS) and Partial Least Squares (PLS) regressions to predict site index using the attributes of site and soil profile, and the attributes of soil physical and chemical properties in Quercus mongolica stands.

Note: SI: site index; TN: total nitrogen

고찰

토양은 임분의 지위지수를 결정하는 요인 중의 하나이다. 신갈나무 임분의 입지 및 토양단면의 속성변수나 토양 물리·화학적 특성의 속성변수는 변수 사이에 상관관계가 강하게 나타나 OLS 회귀모형에 의한 지위추정은 다중공선성이 발생할 수 있는 것으로 나타났다. 이러한 문제점을 제거하기 위해 VIF 값이 높은 속성변수를 제외한 회귀 모형은 A층의 토양 물리·화학적 특성의 속성변수를 이용한 회귀식을 제외하고, 지위지수와 유의적인 회귀 관계를 보였다. VIP로부터 얻어진 변수를 이용한 PLS 회귀모형도 입지 및 토양단면의 속성변수나 토양 물리·화학적 특성의 속성변수를 이용하여 유의적인 지위 추정이 가능하였다. 한편, 수정결정계수는 OLS 회귀모형보다 PLS 회귀모형이 높아 지위 추정 관련 설명력은 PLS 회귀모형이 다소 크게 나타났다. 이는 OLS 회귀모형의 경우 독립변수가 2∼3개 이내였던 반면 PLS 회귀모형의 독립변수는 3∼8개로 독립변수 증가에 따라 잔차제곱합(Error Sum of Squares, SSE)이 감소하여 수정결정계수 값이 증가하였기 때문으로 사료된다(Carrascal et al., 2009). 타 연구에서도 VIP를 이용한 PLS 회귀모형이 OLS 회귀모형의 설명력에 비해 우수하다고 보고된 바 있다(Carrascal et al., 2009; Subedi and Fox, 2016). 현편, 입지환경 및 토양단면과 토양 물리·화학적 특성의 속성변수 모두를 이용한 OLS와 PLS 회귀 모형 모두 지위지수와 유의적인 회귀관계를 보였다. 그러나 회귀모형에 포함되는 변수가 증가하였음에도 결정계수 값은 각각의 속성변수를 이용한 회귀모형과 큰 차이는 나지 않아, 변수 수집의 시간과 경비를 고려할 때 실용적이지는 않을 것으로 사료된다.

본 연구의 입지 및 토양 특성의 속성변수를 이용한 회귀 모형의 설명력은 타 연구와 유사한 수준이었다. 국내에서 조사된 1:25,000 산림입지도의 입지환경 및 토양단면 속성정보를 이용한 신갈나무 임분 지위추정 회귀모형의 결정계수는 온대 북부산림대는 R²=0.30, 온대 중부 산림대는 R²=0.25 정도로 보고된 바 있다(Lee et al., 2007). 한편 굴참나무 임분의 토양 물리·화학적 특성의 속성변수를 이용한 단계적회귀분석 지위추정 모형의 결정계수는 A층 R²=0.32, B층 R²=0.54로 보고된 바 있다(Lee et al., 1996). 상수리나무 임분의 토양 물리·화학적 특성의 속성변수로 부터 주성분분석 회귀모형을 이용한 지위추정 결정계수는 R²=0.40(Jung et al., 2017)로 본 연구 결과보다 약간 높거나 유사한 값을 보이고 있다. 국외에서 조사된 토양 특성의 속성변수를 이용한 참나무류의 지위추정 모형의 경우 포르투갈 co rk o ak (Quercus suber L.) 임분은 R²=0.41 (Paulo et al., 2015), 미국에서 조사된 white oak (Q. alba L.)과 northern red oak (Q. rubra L.)의 지위추정 모형은 각각 R²=0.51과 R²=0.70(Woolery et al., 2002)로 본 연구결과보다 높게 나타났다. 본 연구에서 지위추정 회귀모형의 설명력이 낮게 나타났지만, 이는 입지 및 토양단면이나 토양 물리·화학적 특성의 속성변수가 지위추정 회귀 모형에 영향을 주지 않는다는 의미라기보다는 본 연구가 광범위한 공간적 범위에서 조사되어 해발고, 기후, 지형 등의 차이에 따른 입지 및 토양 특성의 변이 폭이 크게 나타났기 때문으로 사료된다. 타 연구에서도 공간적으로 넓은 범위에서 조사된 토양 특성 및 입지 인자에 의한 지위추정 회귀식의 결정계수는 좁은 범위나 서로 유사한 토양형에서 실시된 지위 추정 모형보다 낮은 값을 보이는 것으로 알려져 있다(Bergès et al., 2005). 또한, 지위지수 추정을 위한 입지 및 토양 특성 분석이 하나의 점(point) 단위인 토양 단면 속성을 기반으로 조사되고 분석되어 임분 단위로 추정되는 지위지수 변이를 올바로 반영하지 못할 가능성도 있다(Landsberg et al., 2003).

신갈나무 임분의 PLS 회귀모형에 포함되는 중요변수로 해발고, 색상, 채도 등과 전질소, 탄질률, P₂O₅, Mg²⁺ 등이 지위 추정 모형에 포함되었다. 입지 및 토양단면 속성 중 해발고는 A층의 깊이나 지위지수 등과 강한 정의 상관관계를 보이고 있으며(Table 4), A층은 식물의 양분 흡수나 수분 흡수와 밀접한 관련이 있는 세근이 많이 분포하여 산림생산력에 영향을 끼칠 수 있기 때문으로 사료된다. 한편 PLS 회귀모형의 중요변수인 토색의 색상(Hue) 회귀계수는 마이너스 값으로, 이는 토색 중 색상은 주로 교환성 양이온함량과 관련이 있으며, 5YR과 같이 적색에 가까운 색상 값이 황색에 가까운 10YR 같은 높은 색상 값보다 교환성 K⁺ 함량이 높게 나타나는(Koné et al., 2016) 등 양분함량과 관계되기 때문으로 사료된다. 토양 물리·화학적 특성의 속성변수 중 회귀계수에 플러스 값을 보인 전질소와 유효 인은 국내 산림토양의 산림생산력을 제한하는 대표적인 양분으로서 알려져 있다(Jeong et al., 2002). 그러나 마이너스 회귀계수를 보인 탄질률은 토양 내 유기물과 결합된 유기 질소를 무기질소로 변화시키는 질소무기화나 미생물에 의한 질소부동화와 관련이 있으며, 높은 탄질률을 보이는 토양은 낮은 탄질률을 보이는 토양에 비해 순 질소무기화량이 낮고, 높은 질소부동화율을 보였다(Janssen, 1996; Gan et al., 2010). 한편, Mooshammer et al.(2014)은 높은 토양 탄질율은 미생물에 의해 부동화된 질소로 인하여 미생물 질소이용효율(microbial nitrogen use efficiency)이 높아 식생이 이용할 수 있는 무기질소량이 제한될 수 있음을 보고한 바 있다.

결론

국내 분포하는 신갈나무 임분은 입지 및 토양단면과 토양 물리·화학적 특성의 속성변수를 이용하여 OLS와 PLS 회귀 모형에 의해 산림생산력 추정이 가능하였으며, OLS 보다 PLS 회귀모형의 설명력이 높게 나타났다. 신갈나무 임분의 지위 추정은 B층의 입지 및 토양단면의 속성변수를 이용한 PLS 회귀모형이 가장 설명력이 높아 토양채취 후 물리·화학적 특성 분석이 필요한 회귀모형과는 다르게, 현장에서 지위지수 추정이 가능할 것으로 사료되었다. 그러나 본 연구에서 제시된 OLS와 PLS 회귀모형은 결정계수 값이 낮아 산림생산력의 정확한 예측에 한계가 있을 수 있어, 현장에서 지위지수 추정에 활용할 수 있으면서도 회귀모형의 설명력을 향상시킬 수 있는 새로운 변수 개발이나, 모암, 토양형, 기후대 등과 같이 공간적인 범위의 축소를 통한 새로운 회귀모형의 개발도 필요할 것으로 사료되었다.