서 론

임도는 산림의 경영 및 관리를 위하여 설치된 도로로서, 산림이라는 공간에 도로를 매개체로 시장 또는 생산 ‧ 생활공간을 연결해주는 기반 시설이다(KFS, 2000). Fig. 1은 임도의 전형적인 횡단구조를 보여주는 것으로서 임도 개설 시 일반적으로 임도 상 ‧ 하부에 절 ‧ 성토면이 생성된다(KFS, 2020). 산림에 임도가 개설되면 유수의 흐름이 바뀌게 되고, 변화된 유수의 흐름과 유량은 임도, 임도비탈면, 인접한 산림의 안정성을 저해시킬 수 있다(Wemple et al., 1996; Choi et al., 2011). 2004년 충청북도 증평군 율리마을의 산사태, 2011년 경상남도 밀양시 양지마을의 산사태, 2013년 강원도 춘천시와 홍천군, 경상남도 밀양시 일대의 산사태, 2014년 전라북도 고창군 성송면의 산사태 등이 모두 임도로부터 유발된 산사태로 보고되어 있다(Lee et al., 2004; Park and Park, 2015; Baek et al., 2016).

산사태는 인간활동(도로건설, 광산개발, 토지이용 등), 물리적 특성, 역학적 특성, 지형적 특성 등에 의해 복합적으로 발생하기도 하지만, 강우나 지진과 같이 하나의 유발요인(triggering factor)에 의해서 발생할 수도 있다(Wieczorek, 1996; Clague and Stead, 2012). 국내에서는 집중호우 시기인 여름철에 산사태가 자주 발생하고 있어서, 산사태 피해 저감을 위해 산사태 발생지점의 위치 및 해당 지점의 물리 ‧ 역학적, 지형적 특성들에 대한 정보를 수집하여 산사태 취약지도를 작성하고, 이를 활용하여 방재대책을 마련할 필요가 있다. 하지만 다수의 연구에서는 광역적인 연구지역에 대한 현장조사의 어려움으로 수치지형도로부터 취득한 지형적 데이터만을 이용하고 있다(Lee et al., 2006; Jo et al., 2007; Quan et al., 2008; Kim and Park, 2013; Lee and Oh, 2019; Urmi et al., 2020). 또한, 물리 ‧ 역학적 특성이 고려된 경우에도 강우 데이터 수집 지점들 간의 거리가 수십 km 이상 이격되어 있어 강우조건이 상이한 조건 하에서 연구된 사례가 많다(Chae et al., 2004; Quan et al., 2011).

본 연구에서는 동일 시간대 반경 2 km 이내 좁은 지역에서 다수의 산사태가 발생하여 강우 조건 및 식생 조건이 동일한 것으로 가정할 수 있는 충주시 산척면 상산마을 일대를 연구지역으로 설정하였다. 산사태가 집중적으로 발생한 연구지역 주변 임도를 대상으로 물리적, 역학적, 지형적 데이터를 수집하였으며, 이를 대상으로 구조방정식 모형과 로지스틱 회귀분석을 실시하여 산사태 주요 영향인자를 선정하고자 한다.

Fig. 1.

Cross section of a typical forest road (modified from KFS, 2020).

연구지역 현황 및 분석자료 수집

강우량

연구지역에서는 2020년 8월 2일에 집중 강우로 다수의 산사태가 발생하였다. 기상청의 공공기관 기상관측자료를 이용하여 연구지역과 4 km 이격된 관측소(송강(3201))의 강우 자료를 확인한 결과, 8월 2일의 일강우량은 343 mm, 최대시우량 74 mm로 과거 산림청의 산사태 예보제 기준인 누적 강우량 200 mm, 일강수량 150 mm, 최대시우량 30 mm를 모두 초과하는 것으로 파악되었다(KMA, 2022; Fig. 2).

Fig. 2.

Hourly and accumulated rainfall around the study area on August 2, 2020.

산사태 발생 현황

연구지역에서 관찰된 산사태 발생지점은 총 40개소이며, 임도비탈면에서 발생한 산사태 36개소, 자연비탈면에서 발생한 산사태 4개소로 전체 산사태의 90%가 임도 비탈면을 따라 발생하였다. 산사태는 대부분 Fig. 3a와 같이 두부 급사면(head scarp)에서 원호파괴 형태로 시작되며, 일부 토층 두께가 얇은 지역에서는 Fig. 3b와 같이 전이형 산사태(translational landslide)가 관찰되기도 한다.

Fig. 3.

Photographs of typical types of landslide along the forest road. (a) Rotational landslide in the head scarp of the lower slope and (b) translational landslide in the upper slope.

분석자료 수집

산사태 영향인자 선정을 위하여 산사태 발생지점 37개소(40개소 중 시료채취가 불가능한 3개소 제외)와 미발생 지점 45개소를 대상으로 토층의 물리 ‧ 역학적 특성과 비탈면의 지형 데이터를 분석하였다(Fig. 4). 토층의 물리 ‧ 역학적 특성은 시료채취 및 실내시험을 수행하여 획득하였으며, 현장에서 채취한 불교란 시료를 대상으로 함수비시험(ASTM D2216-10, 2010), 비중시험(ASTM D854-10, 2010), 직접전단시험(ASTM D3080-98, 1998)을 수행하여 건조단위중량, 포화단위중량, 포화함수비, 공극률, 비중, 내부 마찰각, 점착력을 획득하였다. 지형 데이터는 현장조사와 지형 공간분석을 통해 획득하였으며, 현장조사 시 전토층(regolith)을 측정하거나 콘관입시험기(Solsolution社 Panda 2)를 이용하여 토층의 두께를 측정하였고, 토층을 제외한 지형적 특성은 국토지리정보원의 1:5,000 수치지도를 SAGA GIS를 활용하여 분석함으로써 5 m × 5 m 해상도의 사면의 경사, 측면곡률(profile curvature), 평면곡률(plan curvature), 지형습윤지수(topographic wetness index, TWI)를 획득하였다.

Fig. 4.

Locations of undisturbed sampling and investigation of regolith thickness in the study area.

산사태 영향인자 선정

데이터 전처리

통계 분석을 위한 데이터 전처리는 최소-최대 정규화(min-max normalization), 다중공선성(multicollinearity) 진단 순으로 진행하였다. 수집된 12개의 독립변수(건조단위중량, 포화단위중량, 비중, 공극률, 포화함수비, 마찰각, 점착력, 토층 두께, 경사, 측면곡률, 평면곡률, TWI)를 대상으로 상자-수염 그림을 도시하여 수염범위를 벗어난 이상치를 제거한 다음 최소-최대 정규화를 수행하였다(식 (1)).

상기 식에서 Xn은 정규화가 적용된 값, X는 관측값, Xmin은 관측값의 최솟값, XMax는 관측값의 최댓값을 의미한다. 최소-최대 정규화를 수행하면 모든 데이터의 분포범위가 최솟값 0, 최댓값 1로 정규화되기 때문에 각각의 독립변수들이 종속변수(산사태 발생)에 미치는 영향을 직접적으로 비교할 수 있다.

공선성이란 회귀분석에 사용되는 독립변수 중 일부 독립변수 간의 상관정도가 높을 때 회귀모형 변수의 과대평가, 회귀결과 신뢰도 저하 등 부정적인 영향이 나타날 수 있는 현상이다. 따라서, 회귀분석 이전에 공선성을 필수적으로 확인해야 하며, 공선성 진단은 식 (2)의 분산팽창계수(variation inflation factor, VIF)를 이용하고, 분산팽창계수가 10 이상일 경우 다중공선성이 존재하는 것으로 진단한다(Kutner et al., 2004).

상기 식에서 VIF는 분산팽창계수, R2은 결정계수이다. Table 1은 수집된 12개의 독립변수를 대상으로 다중공선성을 진단한 결과로, 분산팽창계수가 10 이상인 독립변수는 건조단위중량, 포화단위중량, 비중, 공극률, 포화함수비이다. 상관성이 높은 5개의 독립변수는 식 (3), (4), (5)와 같이 독립변수들 간의 상관성을 설명할 수 있다.

여기서, Υd는 건조단위중량, Υsat은 포화단위중량, Gs는 비중, e는 공극률, w는 함수비, S는 포화도이다. 독립변수의 수에 따라 공선성이 달라지므로, 다중공선성이 높은 인자들을 조합하여 모든 독립변수의 분산팽창계수가 10 이하를 만족하도록 하였다. 다중공선성 진단을 통해 최종적으로 건조단위중량, 비중, 포화함수비가 제거되었고, 9개의 독립변수를 이용하여 다중공선성 진단을 수행한 결과는 Table 2와 같다.

구조방정식 모형분석

구조방정식 모형분석(structural equation model analysis, SEM)은 Wright(1921)가 제안한 통계분석 기법으로 복잡한 상호 관련성을 가지고 있는 여러 독립변수들 간의 인과관계를 표현하고, 상관성을 설명하는 분석방법이다. 변수 간의 상관관계를 그림으로 표현하여 직관적으로 이해할 수 있다는 장점이 있지만(Hox and Bechger, 1999; Yung, 2008; Ullman and Bentler, 2013), 신뢰성 있는 결과를 위해서는 분석대상에 대한 사전지식이 필요하고, 추정식을 미제공하므로 변수들의 인과관계만 파악할 수 있다는 단점이 있다.

구조방정식 모형은 측정변수(observed variables), 잠재변수(latent variable), 오차항(error)으로 구성된다. 측정변수는 직접 관측한 변수로 현장조사 및 실내시험 등으로 획득한 데이터가 이에 해당한다. 잠재변수는 변수 간의 인과관계에 대한 이해를 돕기 위해 가상으로 결정되는 변수로, 잠재변수 결정은 분석자의 전문적인 지식이 반영된다. 오차항은 측정변수에 붙어 있는 측정오차(measurement error)와 잠재변수에 붙어 있는 구조오차(structural error)로 구분할 수 있으며, 측정변수로 잠재변수를 100% 설명할 수 없으므로 오차항이 필요하다. 여기서, 다른 변수를 설명하는 변수는 외생변수(exogenous variable)로 명명되며, 다른 변수에 의해 설명이 되어지는 변수는 내생변수(endogenous variable)라고 한다. 외부모형(external model)은 측정변수, 잠재변수, 측정오차로 구성되고, 내부모형(inner model)은 잠재변수와 구조오차로 구성된다. 전체 모형은 외부모형과 내부모형을 모두 포함하는 형태로 이루어져 있다.

구조방정식 모형은 Fig. 5와 같이 설정하였으며, 측정변수는 현장조사, 실내시험, 지형공간분석으로 획득한 자료이고, 잠재변수는 데이터의 특성을 고려하여 물리적 특성, 역학적 특성, 지형적 특성, 산사태 발생여부로 설정하였다. 물리적 특성에 해당하는 측정변수는 포화단위중량, 공극률이며, 역학적 특성은 마찰각, 점착력, 그리고 지형적 특성은 토층 두께, 경사, TWI, 측면곡률, 평면곡률로 설정하였다. Fig. 5에서 각 화살표의 숫자는 구조방정식 모형분석을 수행한 결과로 통계적 가중치를 나타내고, 부호(+, -)는 각각 비례, 반비례를 의미한다. 전체 모형의 신뢰성 평가는 P-value와 전체 모형의 접합도 지수(goodness of fit index, GFI)로 평가한다. 전체 모형의 P-value는 0.001 이하로 P < 0.05 (신뢰수준 95%)를 만족하고, 모형의 적합도 지수는 0.763으로 Sanchez(2013)의 기준 적용 시 신뢰도 높음(high)에 해당하여 측정변수와 잠재변수 관계에 대한 모형은 잘 구성된 것으로 평가된다. Table 3은 산사태 발생에 영향을 미치는 각 측정변수의 영향도이며 외부모형과 내부모형의 계수의 곱으로 계산된다. 토층 두께의 영향도가 0.55로 가장 높고, 공극률, 포화단위중량, 경사, 마찰각, 평면곡률, TWI, 측면곡률, 점착력 순으로 영향도가 높은 것으로 나타난다.

Fig. 5.

Results of SEM analysis. Each arrow’s numerical label indicates its statistical weight; the sign indicates a proportional (positive) or inversely proportional (negative) relationship.

로지스틱 회귀분석

로지스틱 회귀분석은 하나의 종속변수(dependent variable)와 이를 설명하는 다수의 독립변수(independent variable) 간의 상관성 분석으로 사건의 발생 확률을 계산하여 0과 1사이의 값으로 표현하고, 0.5를 기준으로 0.5보다 크면 1, 작으면 0으로 판단한다(Fig. 6). 로지스틱 회귀분석을 통해 사건이 발생할 확률(Pz)은 식 (6), (7)과 같으며, 여기서 Z는 로지스틱 회귀식, α는 상수항, βn은 독립변수(Xn)의 회귀계수를 의미한다.

Fig. 6.

Distribution of a dependent variable after logistic regression.

로지스틱 회귀분석 결과에 대한 신뢰도를 분석하기 위하여 Nagelkerke R-제곱, Hosmor-Lemeshow 검정, 오차행렬(confusion matrix) 검증방법을 이용하였다. Nagelkerke R제곱은 독립변수가 종속변수를 설명할 수 있는 정도를 의미하는 것으로서 0.2 이상인 경우 설명력이 있다고 판단한다. Hosmor-Lemeshow 검정은 회귀모형의 전체적인 적합도를 판단하는 검정으로 유의수준 0.05 보다 큰 경우 설명력이 있는 모형으로 간주된다. 오차행렬은 예측 결과가 실제 결과와 일치하는지를 평가하는 방법으로 오차행렬의 정확도가 높을수록 예측 결과와 실제 결과가 일치도가 높음을 의미하며, 식 (8) 및 Table 4와 같이 표현된다. 여기서, TP는 참 긍정(true positive)으로 실제 긍정값이 긍정으로 예측된 횟수, TN은 참 부정(true negative)으로 실제 부정값이 부정으로 예측된 횟수, FP는 거짓 긍정(false positive)으로 부정값을 긍정으로 잘못 예측한 횟수, FN은 거짓 부정(false negative)으로 긍정값을 부정으로 잘못 예측한 횟수를 의미한다.

전처리 과정으로 선별된 9개 독립변수를 대상으로 로지스틱 회귀분석을 실시하였으며, 독립변수가 다수이기 때문에 독립변수가 2개, 3개, 4개 … 9개인 모든 경우의 모델을 분석하였다. 분석 결과, 7개의 독립변수(공극률, 포화단위중량, 점착력, 토층두께, 사면경사, TWI, 측면곡률)로 구성된 회귀모델의 정확도가 83.9%로 가장 높았으며, Nagelkerke R-제곱이 0.543, Hosmor-Lemeshow 검정 유의 확률이 0.123으로 나타나 유의수준 및 신뢰도 기준도 만족하는 것으로 분석되었다. 상기 회귀모델의 회귀계수와 산사태 발생에 영향을 미치는 각 독립변수의 영향도는 Table 5와 같으며, 공극률의 영향도가 38.71%로 가장 높고, 토층 두께, 포화단위중량, 측면곡률, 점착력, 경사, TWI 순으로 영향도가 높게 나타난다.

주요 영향인자 선정

주요 영향인자를 비교하기 위해서는 동일한 변수들을 비교해야지만, 9개의 변수를 모두 적용한 로지스틱 회귀분석의 오차행렬의 정확도가 낮아 통계적으로 유의미하지 않은 결과를 보였다. 따라서 오차행렬의 정확도가 83.9%로 가장 높은 7개의 변수를 적용한 로지스틱 회귀분석과 구조방정식 모형분석에 사용된 공통된 7개의 변수(공극률, 포화단위중량, 점착력, 토층 두께, 사면 경사, TWI, 평면곡률)를 대상으로 산사태 주요 영향인자를 선정하였다(Table 6). 토층 두께, 공극률, 포화단위중량은 두 분석기법에서 모두 상위권에 자리매김하였으며, 토층 두께와 공극률은 분석방법에 따라 1순위와 2순위가 바뀌어 나타나기도 하지만 상기 3개 영향인자의 영향도 합이 전체의 71~83%로 큰 비중을 차지하여 영향력이 큰 것으로 분석된다. 나머지 독립변수(점착력, 사면 경사, TWI, 평면곡률)의 경우 분석방법에 따라 영향도 크기 및 순위가 다르게 나타난다.

Table 6의 맨 우측 열은 구조방정식 모형분석의 영향도와 로지스틱 회귀분석의 영향도를 산술 평균한 값으로, 영향도 순위는 토층 두께 36.09%, 공극률 29.09%, 포화단위중량 12.18%, 사면경사 6.58%, 측면곡률 3.77%, 점착력 3.09%, TWI 1.54% 순이다. 상위 3개 변수들의 영향도는 77.36%로, 해당 인자들이 연구지역의 산사태 발생에 주된 영향을 미치는 것으로 분석된다. 토층두께가 가장 중요한 영향인자로 분류된 이유는 다른인자에 비해 산사태 발생지점과 미발생 지점의 데이터 차이가 뚜렷하기 때문인 것으로 판단된다. 토층을 제외한 인자들은 산사태 발생지점과 미발생 지점의 평균값이 근소한 차이를 나타내는 반면, 산사태 발생 지점의 평균 토층두께는 2.43 m인 반면 미발생 지점은 1.07 m로 측정되어 발생지점의 토층이 2배 이상 두껍게 분포하는 것으로 분석된다. 또한, 지질공학적 관점에서는 토층의 역학적 특성(마찰각, 점착력)이 산사태 발생에 크게 영향을 미치는 것으로 알려져 있지만(Skempton and Delory, 1957; Kim et al., 2005; Park et al., 2010), 본 연구지역에서는 모든 조사지점의 역학적 특성이 유사한 토질(SW-SP)로 인해 유의한 통계적 차이를 보이지 않기 때문에 영향도가 낮은 것으로 판단된다.

결 론

본 연구에서는 충주시 산척면 상산마을 주변 임도를 대상으로 산사태 발생 지점과 미발생 지점에서 토층의 물리 ‧ 역학적 특성 및 지형적 특성에 관련된 독립변수를 수집하고, 구조방정식 모형분석과 로지스틱 회귀분석을 실시하여 산사태의 주요 영향인자를 선정하였다. 연구지역은 집중호우 기간에 반경 2 km 이내에서 집중적으로 다수의 산사태가 발생한 지역으로서, 강우 조건과 식생 조건이 유사한 것으로 평가되어 그 영향도를 배제할 수 있어 토층의 물리 ‧ 역학적 특성 및 지형적 요인과 관련된 인자들만의 영향도를 파악할 수 있다.

산사태 발생지점 37개소와 미발생 지점 45개소를 대상으로 현장조사, 실내시험, 지형 공간분석 등을 통해 12개의 독립변수(건조단위중량, 포화단위중량, 비중, 공극률, 포화함수비, 마찰각, 점착력, 토층 두께, 사면경사, 측면곡률, 평면곡률, TWI)데이터를 수집하였고, 모든 독립변수의 값을 대상으로 0에서 1 사이에 분포하도록 최소-최대 정규화를 수행하였다. 그 후 다중공선성 진단을 통하여 최종적으로 9개의 독립변수(포화단위중량, 공극률, 마찰각, 점착력, 토층 두께, 사면경사, 측면곡룰, 평면곡률, TWI)를 선정하였다.

구조방정식 모형분석에는 9개의 독립변수가 사용되었으며, 토층 두께의 영향도가 35.37%로 가장 높고, 공극률 19.46%, 포화단위중량 16.98%, 경사 10.82%, 마찰각 10.36%, 평면곡률 5.01%, TWI 1.31%, 측면곡률 0.67%, 점착력 0.02%로 분석되었다. 로지스틱 회귀분석에는 7개의 독립변수가 사용되었으며, 공극률의 영향도가 38.71%로 가장 높고, 토층 두께 36.81%, 포화단위중량 7.37%, 측면곡률 6.86%, 점착력 6.16%, 사면 경사 2.33%, TWI 1.76%로 분석되었다. 공통된 독립변수의 영향도를 산술 평균한 결과, 영향도 순위는 토층 두께 36.09%, 공극률 29.09%, 포화단위중량 12.18%, 사면경사 6.58%, 측면곡률 3.77%, 점착력 3.09%, TWI 1.54% 순이다. 토층 두께, 공극률, 포화단위중량 등 상위 3개 변수의 영향도는 77.36%로서 연구지역의 산사태 발생에 큰 영향을 미치는 것으로 분석된다.