• The relative importance of conditioning factors in landslide susceptibility analysis was evaluated using the SHAP summary plot.

• The spatial pattern analysis based on SHAP values provided a spatially explicit interpretation of the influence of conditioning factors.

산사태 취약성 분석(landslide susceptibility analysis)은산사태 발생 위치와 산사태 발생에 영향을 미치는 다양한인자 간의 공간적 분포를 분석하여 특정 지역에서 산사태 발생 가능성을 예측하는 방법이다 (Park and Lee, 2022).산사태 취약성 분석 기법 중 데이터 기반 분석(data-driven analysis)은 과거 산사태 유발 조건과 유사한 조건에서미래에 산사태가 발생할 것이라는 전제로, 과거 산사태 발생 위치와 산사태 영향인자 간의 상관성을 분석하는 방법이다 (Zhu et al., 2018; Park and Lee, 2022). 최근 데이터기반 산사태 취약성 분석에서는 지질, 지형, 임상, 토양 등산사태 발생에 영향을 미치는 다양한 인자들의 복잡한패턴을 분석하기 위해 Artificial Intelligence (AI) 모델이널리 활용되고 있다 (Catani et al., 2013; Chen et al., 2017; Huang and Zhao, 2018; Wang et al., 2019; Sahin, 2020). 그중 앙상블 모델은 여러 개의 기본 학습기(base learner)를 결합하고, 기본 학습기의 예측 결과를 통합하여 단일 모델보다 더 우수한 예측 성능을 제공하는 앙상블 학습(ensemble learning) 방법론을 활용하는 모델이다 (Sahin, 2020). 앙상블 학습에서 가장 대표적으로 사용되는 결합 방법은 배깅(bagging)과 부스팅(boosting)이다. 배깅은 부트스트랩샘플링 기법을 활용하여 학습 데이터 세트를 무작위로 생성하고, 각 기본 학습기마다 독립적으로 모델 예측 과정을수행하는 방법이며, 부스팅은 이전 학습기의 학습 결과를기반으로 순차적으로 학습을 수행하는 방법이다 (Hong et al., 2018; Dou et al., 2020). 이러한 결합 구조를 포함하는 앙상블 모델은 비선형 데이터 세트에 대한 일반화 성능이 뛰어나 데이터 기반 산사태 취약성 분석에서 효과적으로 활용되고 있다 (Kutlug Sahin and Colkesen, 2021).

그러나 이러한 모델은 방대한 양의 파라미터를 추정해야하므로 복잡하고 비선형적인 구조를 포함하는 블랙박스(black-box) 모델, 즉 모델의 예측 과정을 명확히 알기어려운 모델로 분류된다 (Li, 2022). 예측 과정을 명확히이해하기 어려운 경우, 예측 결과가 도출된 근거를 파악하기 어려워, 예측 결과에 대한 신뢰성이 떨어질 수 있다.이에 따라 AI 모델의 예측 과정을 설명하고 분석 결과의신뢰성 확보를 위한 연구가 활발히 진행되고 있으며, 특히 XAI 기법이 주목받고 있다 (Bach et al., 2015; Samek et al., 2017; Toms et al., 2020; Ullah et al., 2020; Hsu and Li, 2023). XAI 기법은 AI 모델이 특정 결론에 도달한 과정을 이해할 수 있는 방식으로 설명하여 분석 결과의투명성을 높이고 예측 결과에 대한 신뢰성을 강화한다.특히, 데이터 기반 산사태 취약성 분석에서는 XAI 기법 중SHAP 기법이 주로 활용되고 있다. 기존 연구에서는 SHAP value를 활용하여 산사태 취약성 분석에 사용된 인자들의기여도를 정량화하고, 이를 다양한 plot으로 시각화하여각 인자가 분석 결과에 미치는 영향을 평가하는 데 초점을맞추고 있다 (Collini et al., 2022; Dahal and Lombardo, 2023; Pradhan et al., 2023; Sun et al., 2023; Zhang et al., 2023). 그러나 SHAP value는 산사태 취약성 분석에 대한인자들의 영향력을 숫자 값으로만 표현하며, 산사태 발생위치와 산사태 영향인자 간의 공간 분포에 따른 예측과정을 충분히 설명하지 못하는 한계가 있다. 현재까지산사태 취약성 분석 결과를 공간적으로 해석하는 방법에대해서는 상대적으로 연구된 바가 적어, 이에 따라 공간적관점에서 인자들의 영향력을 분석하는 연구의 필요성이제기되고 있다 (Li, 2022; Wang et al., 2024).

본 연구는 앙상블 모델인 AdaBoost, Random Forest, XGBoost를 적용하여 산사태 취약성 분석을 수행하였으며, SHAP 기법을 적용하여 인자의 중요도를 평가하였다.이를 위하여 2020년 7~8월 집중호우로 다수의 산사태가발생한 충청북도 제천시 봉양읍 일대를 연구 지역으로선정하였다. 또한, 중요도가 높은 인자들의 공간적 영향을 분석하기 위해 SHAP value를 공간 형태로 시각화하는 방법을 적용하였다. 이를 통해 인자들이 산사태 취약성에 미친 영향을 공간적 관점에서 분석하고자 하였다.

2.1. AdaBoost

AdaBoost 모델은 부스팅 알고리즘에 기반한 모델로, 여러개의 약한 분류기(weak classifier)를 순차적으로 학습시키고, 약한 분류기의 결과를 통합하여 예측 성능을 향상시키는 알고리즘이다 (Freund and Schapire, 1997)(Fig. 1a). AdaBoost는 초기 학습 단계에서 모든 데이터에 동일한 가중치를 부여하며, 학습 과정에서 발생한 분류 오류를 보완하기 위해 오분류된 데이터에 더 높은 가중치를 부여하는 방식을 사용한다. 이를 통해 다음 학습 단계에서는 오분류된 데이터를 더 집중적으로 학습하여 성능을 개선한다 (Kadavi et al., 2018). 본 연구에서는 AdaBoost의 약한분류기로 의사결정나무(decision tree)를 사용하였으며, 하이퍼파라미터로 약한 분류기의 개수(n_estimators)와 트리의 최대 깊이(max_depth)를 설정하여 분석을 수행하였다.

Figure 1. Diagram of ensemble learning methods: (a) Adaboost, (b) Random Forest and (c) XGBoost (Le et al., 2024; Salcedo Sanz et al., 2024).

2.2. Random Forest

Random Forest 모델은 배깅 알고리즘에 기반하여 여러개의 의사결정나무를 결합한 앙상블 모델이다 (Breiman, 2001)(Fig. 1b). Random Forest는 데이터 샘플링 단계에서부스트스랩 방식을 활용하여 다양한 데이터 세트를 생성하고, 각 데이터 세트를 학습한 의사결정나무의 결과를통합하여 최종 예측값을 도출한다. 이러한 과정은 단일의사결정나무에서 발생할 수 있는 과적합(overfitting) 문제를 완화하는 데 기여한다 (Wang et al., 2021). 본 연구에서는 Random Forest의 하이퍼파라미터로 트리의 개수(n_estimators), 트리의 최대 깊이(max_depth), 노드를분할하는 데 필요한 최소 샘플 개수(min_samples_split),리프 노드에 있어야 하는 최소 샘플 개수(min_samples_leaf)를 설정하여 분석을 수행하였다.

2.3. XGBoost

XGBoost 모델은 부스팅 알고리즘에 경사 하강법(gradient descent)을 결합하여 개발된 고성능 앙상블 모델이며, 기본학습기로 의사결정나무를 사용한다 (Chen and Guestrin, 2016)(Fig. 1c). XGBoost는 경사 하강법을 활용하여 손실함수(loss function)를 정의하고, 이를 최소화하는 방향으로학습을 진행한다. 이 과정에서 정규화(regularization) 항을 포함한 손실 함수를 사용하여 모델 복잡도를 제어하고과적합을 방지한다 (Shua et al., 2024). 본 연구에서는XGBoost의 하이퍼파라미터로 학습률(eta), 리프 노드의추가 분할 여부를 결정하는 손실 함수의 최소 감소 값(gamma), 트리의 최대 깊이(max_depth), 리프 노드에 포함되는 최소 관측치의 수(min_child_weight), 학습 데이터의 샘플링 비율(subsample), feature 샘플링 비율(colsample_bytree), L1 정규화 항(reg_alpha)을 설정하여 분석을 수행하였다.

2.4. 모델 최적화

본 연구에서는 하이퍼파라미터 최적화를 위해 TPE (Tree-structured Parzen Estimators) 알고리즘을 선택하였다. 기존의 하이퍼파라미터 탐색 방법으로는 그리드 탐색(grid search)과 랜덤 탐색(random search)이 주로 활용된다.그러나 그리드 탐색은 차원이 증가할수록 탐색해야 할조합의 수가 기하급수적으로 증가하는 문제가 있으며, 랜덤탐색은 탐색 공간에서 임의의 샘플을 선택하는 방식이므로높은 정확도를 도출해 내기 어려운 한계가 있다 (Choi et al., 2020; Lee et al., 2021).

반면, TPE 알고리즘은 기존 방식과 달리 성능이 좋은하이퍼파라미터와 그렇지 않은 하이퍼파라미터를 분리하여각각의 분포를 학습한 뒤, 성능이 우수한 영역을 집중적으로 탐색하는 방식이다. 이를 통해 탐색 효율을 높이고,고차원 공간이나 범주형 및 조건부 하이퍼파라미터가 포함된 경우에도 효과적으로 최적의 조합을 찾을 수 있다(Bergstra et al., 2011).

본 연구에서 사용한 모델은 비선형 구조를 가지며, 탐색해야 할 하이퍼파라미터 공간이 고차원이므로, 그리드또는 랜덤 탐색보다 탐색 공간 내에서 중요한 변수들을효과적으로 조정할 방법이 필요하였다. 따라서 본 연구에서는 TPE 알고리즘을 적용하여 주요 하이퍼파라미터의최적 조합을 도출하였다(Table 1).

Table 1 . Hyperparameters of AdaBoost, Random Forest and XGBoost models.

Model	Selected Hyperparameters	Range of Search
AdaBoost	n_estimators = 300 max_depth = 3	[50, 100, 150, 200, 300] [3, 4, 5, 6, 7, 8, 9, 10, 15]
Random Forest	n_estimators = 200 max_depth = 9 min_samples_split = 10 min_samples_leaf = 10	[50, 100, 150, 200, 300] [2, 3, 4, 5, 6, 7, 8, 9, 10] [2, 4, 6, 8, 10] [2, 4, 6, 8, 10]
XGBoost	eta = 0.0407 min_child_weight = 1 gamma = 0.1459 max_depth = 4 subsample = 0.8864 colsample_bytree = 0.4599 reg_alpha = 0	(0.001, 0.5) [1, 2, 3, 4, 5, 6, 7, 8] (0.1, 0.5) [2, 3, 4, 5, 6, 7, 8, 9, 10] (0.0, 0.9) (0.0, 0.9) [0, 0.01, 0.05, 0.1, 0.5]

3.1. SHAP

SHAP은 게임 이론에서 유래된 개념으로, 머신러닝모델의 예측 과정을 이해하기 위해 사용되는 XAI 기법이다 (Lundberg and Lee, 2017; Iban and Bilgilioglu, 2023). SHAP은 각 입력 변수가 모델의 예측에 기여하는정도를 정량화하여, 모델의 예측 결과를 각 입력 변수의SHAP value의 합으로 설명한다 (Zhang et al., 2023). 변수i의 SHAP value(φ_i)는 다음 식 (1)과 같이 정의된다.

$${\varphi }_i=\sum _{S\subseteq N\backslash \left\{i\right\}}\frac{\left|S\right|!\left(N-\left|S\right|-1\right)!}{N!}\left[v\left(S\cup \left\{i\right\}\right)-v\left(S\right)\right]$$(1)

여기서 N은 예측에 사용된 입력 변수의 전체 집합을, S는 N에서 변수 i를 제외한 모든 부분 집합을 나타낸다. v(S)는 S에 포함된 변수들로 계산된 기여도를, $v\left(S\cup \left\{i\right\}\right)$는 S에 변수 i를 추가했을 때의 기여도를 의미한다. 따라서 변수i의 SHAP value(φ_i)는 모델 예측값에서 변수 i가 포함되었을 때와 포함되지 않았을 때의 기여도 차이를 모든 변수조합에 대해 평균 내어 계산한 값이다.

SHAP과 함께 XAI 기법으로 주로 활용되는 LIME (Local Interpretable Model-agnostic Explanations)은 개별(local)예측값에서 변수 기여도를 양(+)과 음(−)의 방향으로 설명할 수 있으나, 전역적(global) 변수 중요도를 분석하는데 한계가 있다 (Ribeiro et al., 2016). 또한, Permutation Importance 기법은 변수의 전역적 중요도를 평가하는 데유용하지만, 개별 예측값에서 특정 변수가 예측 결과에미치는 영향의 방향성을 분석하는 데 한계가 있다 (Altmann et al., 2010). 이에 비해 SHAP은 국소적 및 전역적 설명을모두 제공할 수 있어, 개별 예측값뿐만 아니라 전체데이터 세트를 대상으로 포괄적인 해석이 가능하다. 또한,특정 변수가 예측값에 미치는 영향 방향성을 정량적으로 분석할 수 있다. 따라서 본 연구에서는 이러한 SHAP의특성을 활용하여 산사태 발생 영향인자의 중요도를 평가하는 것에 더하여, SHAP value를 공간적으로 시각화함으로써 연구 지역 내에서 특정 변수가 산사태 취약성을증가 또는 감소시키는 데 미친 영향을 분석하고, 이를 통해공간적 패턴을 분석하는 데 중점을 두었다.

SHAP은 분석에 사용하는 모델 종류에 따라 다양한 알고리즘이 존재하며, 대표적으로 KernelSHAP, TreeSHAP, DeepSHAP 등이 있다. 본 연구에서는 트리 기반 모델과호환성이 좋은 TreeSHAP 알고리즘을 적용하였다 (Lundberg et al., 2018).

4.1. 연구 지역

본 연구에서는 2020년 7~8월 집중호우로 인해 211개의산사태가 발생한 충청북도 제천시 봉양읍 일대를 연구지역으로 선정하여 산사태 취약성 분석을 수행하였다(Figs. 2a & b). 연구 지역은 지리 좌표상 위도 37° 03′54″ ~ 37° 09′ 37″, 경도 128° 03′ 11″ ~ 128° 10′ 31″에위치하고, 가로 10.87 km, 세로 10.59 km의 범위를 포함하며, 연구 지역의 전체 면적은 72.41 km²이다. 해발고도는175 m에서 688 m 사이에 분포하며, 평균 해발고도는 296 m로 구릉지형의 특성을 보인다. 지질학적으로는 연구 지역의중앙부에 화강암과 제4기 충적층이 분포하며, 북쪽 지역에는 편마암이 주로 분포한다. 특히, 화강암은 쥐라기의흑운모 화강암과 백악기의 반상 화강암으로 구성되어 있다(Kim et al., 1967). 데이터 기반 산사태 취약성 분석은 과거산사태 발생 위치 정보를 기반으로 수행되므로, 산사태현황도(inventory map)의 구축이 필수적이다. 이에 따라본 연구에서는 산사태 발생 위치 정보를 획득하기 위하여Google Earth를 활용해 산사태 발생 전후 시기인 2019년6월에서 2021년 2월까지 Airbus와 Maxar Technologies에서 제공하는 위성사진을 비교함으로써 산사태 발생 위치를확인하였으며, 이를 바탕으로 점 형태의 산사태 현황도를구축하였다(Fig. 2b).

Figure 2. Study area: (a) location of Jecheon area and (b) shaded relief map of Jecheon area with landslide locations.

4.2. 공간 데이터베이스 및 Mapping Unit

산사태는 영향인자들의 복잡한 상호작용에 의해 발생하기 때문에, 산사태 취약성 분석을 수행하기 위해 연구지역으로부터 산사태 발생과 관련된 영향인자들을 획득하고 이를 기반으로 공간 데이터베이스를 구축하는 과정이필수적이다. 본 연구에서는 국토지리정보원에서 제공하는1:5,000 축척의 수치지형도를 활용하여 DEM (Digital Elevation Model)을 생성하였으며, DEM을 통해 고도(altitude), 사면경사(slope angle), 사면방향(aspect), 평면곡률(planform curvature), 측면곡률(profile curvature),표준곡률(standard curvature), TWI (Topographic Wetness Index), SPI (Stream Power Index)의 지형 인자를 획득하였다. 획득한 지형 인자 중, 평면곡률은 최대 경사방향에수직인 방향의 곡률을 의미하며, 측면곡률은 최대 경사방향의 곡률을 의미한다 (Ko et al., 2014). 또한, 표준곡률은 평면곡률과 측면곡률을 통합한 개념이다. TWI는 지형형태에 따른 물의 흐름과 집적을 나타내는 인자이며, SPI는흐르는 물의 침식력을 나타내는 인자이다 (Guo et al., 2015; Kalantar et al., 2018). 토양 인자는 산림청에서 제공하는 1:5,000 축척의 산림입지도를 통해 토양심도(soil depth) 및 토양속성(soil texture) 인자를 획득하였고, 임상 인자는 산림청에서 제공하는 1:5,000 축척의 수치임상도를기반으로 경급(timber diameter), 영급(timber age), 임상종류(forest type), 임상밀도(forest density), 수목높이(forest height)를 획득하였다. 지질 인자는 한국지질자원연구원에서 제공하는 1:50,000 축척의 수치지질도를 통해 암종(lithology)을 획득하였으며, 토지이용 인자는 환경부에서제공하는 1:5,000 축척의 토지피복지도를 활용하였다.

한편, 산사태 취약성 분석에서는 지형 공간을 구분하는기준인 mapping unit의 선정이 필수적이다. 대표적인 mapping unit으로는 grid unit과 slope unit이 있다. Grid unit은 지형 공간을 동일한 크기의 격자로 분할하여 각격자에서 산사태 발생 가능성을 예측하는 방식으로, 상대적으로 모델링이 단순하다는 장점이 있다. 그러나 고정된 격자 구조로 인해 지형 경계가 인위적으로 나뉘어 실제사면의 형태를 충분히 반영하지 못하는 한계가 존재한다.반면, slope unit은 배수선(drainage line)과 분수선(divide line)을 기준으로 지형 공간을 분할하여, 지형학적 및 수문학적 특성을 보다 효과적으로 반영할 수 있다 (Schlögel et al., 2018). 그러나 복잡한 지역에서는 slope unit의경계가 명확하지 않을 수 있으며, grid unit과는 달리 slope unit의 크기는 지형 조건에 따라 달라질 수 있다는 점을고려해야 한다.

따라서 본 연구에서는 실제 지형 공간의 형태를 고려하여 분석을 수행하기 위해 slope unit을 기반으로 공간데이터베이스를 구축하였다 (Erener and Düzgün, 2012)(Fig. 3). 래스터 형식으로 구축된 데이터베이스를 slope unit 기반의 데이터베이스로 변환하는 과정에서 하나의slope unit 내 여러 인자 속성값이 포함되므로, slope unit의 인자 속성값을 결정하기 위한 기준이 필요하다 (Chang et al., 2023). 연속형 인자의 경우, 기존 연구에서는 slope unit에 인자의 평균값을 부여하는 방식을 주로 사용하였다 (Ba et al., 2018; Erener and Düzgün, 2012). 그러나이러한 방식은 slope unit 내 인자의 다양한 속성값을 충분히 반영하지 못해 지나치게 단순화되는 문제가 발생할수 있다. 이러한 한계를 보완하기 위해, 최근 연구에서는slope unit 내 변동성을 고려하는 방법으로 평균값과 함께 표준편차 값을 사용하는 방법이 제안되고 있다 (Chang et al., 2023). 예를 들어 평면곡률은 지형에 따라 물의 흐름이 수렴하는지 발산하는지를 나타내며, 형태가 볼록한경우 양수, 오목한 경우 음수, 평평한 경우 0의 값을 갖는다 (Guo et al., 2015). 따라서 평면곡률 표준편차가 큰경우, 하나의 slope unit 내에서 곡률의 변동성이 크다는의미이며, 해당 지역이 계곡 또는 능선일 가능성이 높음을의미한다. 본 연구에서도 이러한 접근법을 적용하여 각slope unit에 연속형 인자인 고도(Figs. 4a & b), 사면경사(Figs. 4c & d), 평면곡률(Figs. 4e & f), 측면곡률(Figs. 4g & h), 표준곡률(Figs. 4i & j), TWI (Figs. 4k & l), SPI (Figs. 4m & n)의 평균값과 표준편차 값을 부여하였다. 한편, 범주형 인자의 경우, 평균과 표준편차를 계산할 수 없으므로이러한 방법을 적용할 수 없다. 이에 따라 본 연구에서는 기존 연구에서 활용된 방식인 최빈값을 이용하여 각slope unit에 인자 속성값을 부여하였다 (Chang et al., 2023).범주형 인자로는 사면방향(Fig. 4o), 토양심도(Fig. 4p), 토양속성(Fig. 4q), 경급(Fig. 4r), 영급(Fig. 4s), 임상종류(Fig. 4t), 임상밀도(Fig. 4u), 수목높이(Fig. 4v), 암종(Fig. 4w),토지이용(Fig. 4x)을 사용하였다. 이를 바탕으로 총 24개의 공간 데이터베이스를 구축하여 산사태 취약성 분석에사용하였다.

Figure 3. Shaded relief map of Jecheon area with slope units.

Figure 4. Thematic maps of continuous and categorical landslide conditioning factors: (a) mean of altitude, (b) standard deviation of altitude, (c) mean of slope, (d) standard deviation of slope, (e) mean of planform curvature, (f) standard deviation of planform curvature, (g) mean of profile curvature, (h) standard deviation of profile curvature, (i) mean of standard curvature, (j) standard deviation of standard curvature, (k) mean of SPI, (l) standard deviation of SPI, (m) mean of TWI, (n) standard deviation of TWI, (o) aspect, (p) soil depth, (q) soil texture, (r) timber diameter, (s) timber age, (t) forest type, (u) forest density, (v) forest height, (w) lithology and (x) landuse.

마지막으로, 산사태 발생 위치는 slope unit 내에 산사태가 하나 이상 발생한 경우를 산사태 발생 unit으로정의하였으며, 총 116개의 산사태 발생 unit을 획득하였다(Fig. 5).

Figure 5. Landslide inventory map of slope units.

4.3. 분석 데이터 세트 구축

본 연구에서는 데이터 기반 산사태 취약성 분석을 수행하기 위해, 앞서 획득한 데이터를 학습 데이터(train data)와 검증 데이터(test data)로 분할하여 분석 데이터세트를 구축하였다. 연구 지역에서 획득한 데이터는 산사태 발생부 116개와 산사태 미발생부 2,431개로 산사태미발생부의 개수가 더 많은 불균형 데이터이다. 이러한불균형 데이터의 특성상, 전체 데이터를 무작위로 분할할 경우, 학습 데이터가 산사태 미발생부에 편향될 가능성이 높다 (Gupta and Shukla, 2023). 따라서 본 연구에서는 산사태 발생부와 미발생부의 비율이 1:1이 되도록무작위로 샘플링하는 기존 연구의 방식을 적용하여 각각116개씩 샘플링한 뒤, 이 데이터를 70:30의 비율로 학습데이터(70%)와 검증 데이터(30%)로 분할하였다 (Sun et al., 2023; Zhang et al., 2023). 또한, 샘플링되지 않은 산사태 미발생부 데이터는 모두 검증 데이터에 포함하였다.최종적으로 학습 데이터는 산사태 발생부 81개와 미발생부81개로 총 162개로 구성되었으며, 검증 데이터는 산사태발생부 35개와 미발생부 2,350개로 총 2,385개로 구성되었다.

본 연구에서는 slope unit을 기반으로 데이터 세트를 구축하였으며, 산사태 취약성 분석에 앙상블 학습 기반의AdaBoost, Random Forest, XGBoost 모델을 적용하여 산사태 취약지수(Landslide susceptibility index)를 0에서 1사이의 연속적인 값으로 획득하였다. 획득한 산사태 취약지수를 0.2 단위로 재분류하여 very low (0~0.2), low (0.2~0.4), moderate (0.4~0.6), high (0.6~0.8), very high (0.8~1.0)의 다섯 단계로 구분한 산사태 취약성도를 구축하였다(Fig. 6).

Figure 6. Landslide susceptibility maps: (a) AdaBoost, (b) Random Forest and (c) XGBoost (Zhou et al., 2021; Ahn et al., 2023).

AdaBoost 모델을 활용하여 구축된 산사태 취약성도의경우, very low와 very high 등급의 면적 비율이 0%로 나타나 산사태 취약성이 low~high 등급으로만 평가되었다(Table 2; Fig. 6a). 이는 AdaBoost 모델이 산사태 위험이매우 낮거나 매우 높은 지역에 대해 예측하지 못했음을 의미하며, 결과적으로 산사태 취약성을 효과적으로 분류하는 데 한계가 있었던 것으로 판단된다. 반면, Random Forest 모델과 XGBoost 모델은 very low부터 very high등급까지 모두 예측하였다. 또한, 두 모델은 등급별 면적비율이 유사한 값을 보였는데, 이는 두 모델이 트리 기반 알고리즘의 특성을 바탕으로 산사태 발생부와의 관계를 유사하게 분석한 것으로 판단된다(Table 2). 또한, Random Forest 모델과 XGBoost 모델을 활용하여 구축된 산사태 취약성도에서는 실제 산사태가 발생한 연구지역의 남동부 구릉 지역을 중심으로, 산사태 취약성이high와 very high 등급인 지역이 집중적으로 분포하였다(Figs. 6b & c). 이러한 결과는, Random Forest와 XGBoost모델이 과거 산사태 발생부와 산사태 영향인자 간의 상관성을 효과적으로 분석했음을 보여준다.

Table 2 . Percentage of area for landslide susceptibility classes.

	AdaBoost		Random Forest		XGBoost
	Number of slope units	Percentage	Number of slope units	Percentage	Number of slope units	Percentage
Very high	0	0	147	5.77	169	6.64
High	322	12.64	243	9.54	225	8.83
Moderate	1286	50.49	607	23.83	530	20.81
Low	939	36.87	965	37.89	916	35.96
Very low	0	0	585	22.97	707	27.76

산사태 취약성 분석 모델의 성능을 비교하기 위해 모델 검증에 널리 사용되는 confusion matrix 기반의 ROC (Receiver Operating Characteristics) 그래프의 하단 면적인 AUC (Area Under Curve)를 활용하였다 (Broeckx et al., 2018). Confusion matrix는 위험한 것으로 예측된 지역에서 실제 산사태가 발생한 경우를 나타내는 TP (True Positive), 실제 산사태가 발생하지 않았지만, 위험한 것으로 예측된 경우를 나타내는 FP (False Positive), 안전한 것으로 예측된 지역에서 실제 산사태가 발생한 경우를 나타내는 FN (False Negative), 그리고 실제 산사태가발생하지 않았고, 안전하다고 예측된 경우를 나타내는TN (True Negative)로 구성된다. Confusion matrix를 기반으로 TPR (True Positive Rate)과 FPR (False Positive Rate)을 계산하여 두 지표 간의 관계를 시각화한 ROC 그래프를 작성하였다. TPR(식 (2))은 모델의 예측 성능을나타내는 지표이며, FPR(식 (3))은 모델의 예측 오류를나타내는 지표이다.

$$TPR=\frac{TP}{TP+FN}$$(2)

$$FPR=\frac{FP}{FP+TN}$$(3)

모델의 성능은 ROC 곡선 하단 면적인 AUC를 통해 평가되었으며, AUC 값이 클수록 모델의 예측 성능이 우수함을 의미한다. 본 연구에서 AdaBoost 모델의 AUC는87.51%, Random Forest 모델의 AUC는 85.69%, XGBoost모델의 AUC는 89.92%로 계산되었으며, 이를 통해 XGBoost 모델이 가장 우수한 예측 성능을 보였음을 확인하였다(Fig. 7).

Figure 7. Receiver operating characteristic curves of AdaBoost, Random Forest and XGBoost models.

6.1. 인자 중요도

앞서 제시된 ROC 기반 AUC와 같은 모델 평가지표는모델의 성능을 정량적으로 평가하는 데 유용하지만, 예측값을 도출하는 과정에서 각 인자가 모델의 예측 과정에 어떻게 기여했는지와 같은 정보를 파악하기 어려운한계가 있다. 이에 본 연구에서는 SHAP 기법을 활용하여 각 인자가 산사태 취약성에 기여한 정도를 정량화하고, 이를 바탕으로 인자 중요도(feature importance)를 평가하였다. 또한, Random Forest와 XGBoost 모델이 산사태 발생부와의 관계를 유사하게 분석한 것으로 판단됨에따라, 두 모델이 모델 예측 과정에서 각 인자를 어떻게판단하였는지에 중점을 두어 분석하기 위해, AdaBoost모델을 제외한 Random Forest와 XGBoost 모델을 대상으로 인자의 중요도 평가를 수행하였다.

6.1.1. SHAP summary plot

SHAP summary plot에서 가로축은 SHAP value를,세로축은 산사태 취약성 분석에 사용된 인자를 의미한다.각 점은 인자의 개별 속성을 나타내며, 점의 위치는 해당 속성이 갖는 SHAP value를 나타낸다(Fig. 8). 점의 색상은 인자의 속성값을 나타내며, 파란색 점은 속성값이작음을, 빨간색 점은 속성값이 큼을 의미한다. 예를 들어, Fig. 8b에서 mean slope의 경우, 빨간색 점은 경사 값이크다는 것을, 파란색 점은 경사 값이 작음을 의미한다.또한, SHAP value는 음수에서 양수까지의 범위를 갖는데, SHAP value가 양수인 경우 해당 인자가 산사태 발생 가능성을 높이는 데 기여함을, 음수인 경우 산사태 발생 가능성을 낮추는 데 기여함을 의미한다. 예를 들어, Fig. 8b에서 mean slope는 경사 값이 클수록 양의 SHAP value를 보여, 산사태 발생 가능성을 높이는 데 기여하고,경사 값이 작을수록 음의 SHAP value를 보여 산사태 발생 가능성을 낮추는 데 기여한다. 따라서 summary plot을 통해 인자의 속성값에 따라 산사태 발생 가능성을 증가시키는지, 또는 감소시키는지를 확인할 수 있다.

Figure 8. SHAP summary plots: (a) Random Forest and (b) XGBoost.

Random Forest 모델과 XGBoost 모델의 SHAP summary plot을 비교한 결과, 중요도가 높은 인자는 summary plot의 상위에 배치되었으며, 두 모델에서 중요도가 높게 평가된 상위 4개 인자 그룹(A, B)은 모두 수목높이, 평면곡률 표준편차, 표준곡률 표준편차, 고도로 확인되었다. 또한, 상위 4개 인자의 SHAP value의 분포 양상 또한 유사한 경향을 보였다(Fig. 8). 먼저, 수목높이는 상위 4개인자 중 가장 높은 중요도를 보였다. 낮은 수목높이(파란색 점)는 양의 SHAP value를 나타내어 특정 지역에서 수목높이가 낮을수록 산사태 발생 가능성을 높이는 데 기여하였다. 다음으로, 평면곡률 표준편차와 표준곡률 표준편차는 값이 클 때(빨간색 점) 모두 양의 SHAP value를나타내며, 산사태 발생 가능성을 높이는 데 기여하였다.특히, 평면곡률 표준편차는 값이 증가함에 따라 SHAP value가 음수에서 양수로 전환되는 경향을 보여, 평면곡률표준편차가 커질수록 산사태 발생 가능성을 높이는 데기여하는 인자로 해석되었다. 고도의 경우, 낮은 고도(파란색 점)에서는 음의 SHAP value를 나타내어 산사태 발생가능성을 낮추는 데 기여하였다. 반면, 중간 고도(보라색점)에서는 양의 SHAP value를 보여 산사태 발생 가능성을높이는 데 기여하였으며, 고도가 높은 지점(빨간색 점)에서는 주로 음에서 0 근처의 SHAP value를 보여, 산사태발생 가능성이 크게 증가하지 않는 것으로 분석되었다.

이처럼 두 모델에서 동일한 4개의 인자가 상위 중요도를 보였으며, 해당 인자들의 SHAP value 분포 양상 또한 유사하게 나타났다. 그러나 ROC 기반 AUC 분석 결과에서는 XGBoost 모델이 Random Forest 모델보다 더 높은예측 성능을 보였다. 이를 SHAP summary plot을 통해분석해 보면, Random Forest 모델은 상위 8개 인자를 제외한 나머지 인자들의 영향력이 거의 0에 수렴하는 반면, XGBoost 모델은 상위 18개 인자까지 일정한 영향력이관찰되었다. 즉, XGBoost 모델이 특정 인자에 편향되지않고, 다양한 인자의 영향을 균형 있게 반영하여 예측 정확도를 향상시킨 것으로 판단된다.

6.2. Spatial Pattern Analysis

최근 산사태 취약성 분석에서는 다양한 SHAP plot을통해 모델의 예측 결과를 설명하려는 연구가 활발히 진행되고 있다. 그러나 기존 연구는 SHAP value를 단순히값의 크기 관점에서만 분석하며, 인자들의 영향력을 공간적 관점에서 분석하지 못했다는 한계를 가지고 있다.이에 본 연구에서는 SHAP value를 공간적으로 시각화하여, 특정 인자가 연구 지역 내에서 산사태 취약성에 어떤 영향을 미쳤는지를 분석하고자 하였다. Summary plot의분석 결과를 통해 다양한 인자의 영향을 가장 효과적으로반영하여 분석을 수행한 것으로 판단된 XGBoost 모델을 기반으로 중요도가 높은 인자의 공간 패턴 분석을 수행하였다. 앞선 XGBoost 모델의 summary plot 분석 결과,수목높이, 평면곡률 표준편차, 표준곡률 표준편차, 고도가 중요도가 높은 인자로 확인되었다. 이 중 평면곡률과표준곡률은 모두 곡률의 특성을 나타내기 때문에, 더 높은중요도를 보인 평면곡률 표준편차에 대해서만 공간 패턴분석을 수행하였다. 따라서 수목높이, 평면곡률 표준편차,고도에 대한 공간 패턴 분석을 진행하였다(Figs. 9a, 10a & 11a). 또한, 인자의 개별 속성값이 모델의 예측 결과에 미친영향을 구체적으로 파악하기 위해 single-factor dependence plot을 함께 활용하였다(Figs. 9b, 10b & 11b). Single-factor dependence plot은 특정 인자의 속성값과 해당 인자의SHAP value 간의 관계를 시각화한 그래프로, 가로축은특정 인자의 속성값을, 세로축은 해당 인자의 SHAP value를 나타낸다. 이를 통해 각 인자의 특성이 산사태발생 가능성에 미친 영향을 보다 정량적으로 분석하였다.

Figure 9. Forest height via SHAP values: (a) spatial effects and (b) single-factor dependence plot.

Figure 10. The standard deviation of planform curvature via SHAP values: (a) spatial effects and (b) single-factor dependence plot.

Figure 11. The mean of altitude via SHAP values: (a) spatial effects and (b) single-factor dependence plot.

6.2.1. 수목높이

수목높이는 해당 지역의 평균적인 수목의 높이를 나타내는 인자이다. 수목높이 SHAP value의 공간 분포를 나타낸 Fig. 9a에서 빨간색은 양의 SHAP value를 나타내며 산사태 발생 가능성을 높이는 데 기여함을, 파란색은음의 SHAP value를 나타내며 산사태 발생 가능성을 낮추는 데 기여함을 의미한다. 빨간색으로 표시된 지역은summary plot에서 확인된 수목높이가 낮은 지역과 일치하며, 이러한 지역은 과거 산사태 이력과 높은 연관성을보였다. 추가로 single-factor dependence plot을 통해 수목높이와 SHAP value 간의 관계를 분석한 결과, 수목높이가 0~1 m인 지역에서만 모두 양의 SHAP value를보였다(Fig. 9b). 즉, 빨간색으로 나타나는 지역의 수목높이가 0~1 m인 지역임을 의미하며, 해당 지역이 산사태발생 가능성을 높이는 데 기여함을 보여준다. 이는 일반적으로 수목의 높이가 높을수록 나무의 크기가 크고 뿌리가많이 발달되어 물의 흡수 능력이 좋고 토양을 고정시켜산사태 발생 가능성을 낮추는 데 기여하는 반면, 수목높이가 낮은 지역에서는 이러한 능력이 부족하여 산사태발생 가능성을 높이는 데 기여한 것으로 분석된다.

6.2.2. 평면곡률 표준편차

평면곡률 표준편차의 SHAP value를 공간적으로 분석한결과, 연구 지역 중 산악지형이 분포하는 북서부와 남동부에서는 주로 양의 SHAP value가 나타난 반면, 비교적평탄한 지형은 대체로 음의 SHAP value를 보이는 것으로확인되었다(Fig. 10a). 이를 통해 앞서 summary plot에서양의 SHAP value를 보였던 지역이 실제로 산악지역에해당함을 확인하였다. 추가로 single-factor dependence plot을 활용하여 평면곡률 표준편차와 SHAP value 간의관계를 분석한 결과, 평면곡률 표준편차가 약 1.1 이하인경우 SHAP value는 음수를, 약 1.1에서 1.7 사이인 경우0에 근접한 값을, 약 1.7 이상인 경우에는 양의 SHAP value를 나타냈다(Fig. 10b). 이를 통해 Fig. 10a에서 SHAP value가 빨간색으로 나타나는 지역은 평면곡률 표준편차가1.7 이상인 산악지역임을 확인할 수 있었다. 이러한 결과는평면곡률 표준편차가 큰 계곡 지형에서 집중호우 시강우가 집수되어 지반에 흡수되는 물의 양이 증가하기때문에 산사태 발생 가능성을 높이는 데 기여한 것으로분석된다.

6.2.3. 고도

고도의 SHAP value를 공간적으로 분석한 결과, 연구지역 내 산악지형에서는 전반적으로 양의 SHAP value가나타났으나, 북서부 산악지형에서는 양의 SHAP value와음의 SHAP value가 동시에 나타나는 것으로 확인되었다(Fig. 11a). 이를 통해 앞서 summary plot에서 고도가 높아음의 SHAP value를 나타내는 지역이 주로 북서부 산악지역에 해당함을 확인하였다. 추가로 single-factor dependence plot을 활용하여 고도와 SHAP value 간의 관계를 분석한결과, 고도가 265 m 이하이거나 400 m 이상인 지역에서는음의 SHAP value를 나타내어 산사태 발생 가능성을낮추는 데 기여하는 것으로 나타났다. 반면, 265 m에서400 m 사이에 해당하는 지역에서는 양의 SHAP value를보이며 산사태 발생 가능성을 높이는 데 기여하는 것으로분석되었다(Fig. 11b). 따라서, 중간 고도는 산사태 발생가능성을 높이는 데 기여하며, 실제로 과거 산사태가 발생한 위치도 주로 중간 고도의 지역에 분포하고 있음을확인하였다. 이러한 경향은 높은 고도에서는 기반암의 노출이 많아 산사태의 대상이 되는 토층이 얕게 분포하여산사태 발생 가능성을 낮추는 데 기여하는 반면, 중간 고도에서는 산사태 발생에 취약한 토양층이 주로 분포하기때문으로 판단된다 (Kim et al., 2007).

본 연구에서는 제천시 봉양읍 일대를 대상으로 AdaBoost, Random Forest, XGBoost 모델을 활용하여 산사태 취약성분석을 수행하였다. 분석 결과를 바탕으로 산사태 취약성도를 구축하였으며, 이를 통해 과거 산사태 발생부와산사태 영향인자 간의 상관성을 유사하게 분석한 것으로판단된 Random Forest와 XGBoost 모델에 SHAP 기법을적용하여 summary plot을 통해 산사태 영향인자들의 중요도를 평가하였다. summary plot 분석 결과, 두 모델 모두에서 수목높이, 평면곡률 표준편차, 표준곡률 표준편차,고도가 산사태 취약성에 가장 큰 영향을 미치는 주요인자로 확인되었다. 또한, summary plot의 분석 결과를바탕으로 다양한 인자의 영향을 가장 효과적으로 반영하여분석을 수행한 것으로 판단된 XGBoost 모델을 대상으로, 중요도가 높게 평가된 인자의 SHAP value를 공간적으로변환하여 지도 형태로 시각화하였다. 이를 통해 인자들의 공간 분포에 따른 기여도를 직관적으로 확인할 수 있었으며, 각 인자가 연구 지역 내 산사태 취약성에 미친영향을 공간적으로 분석하였다. 공간 패턴 분석 결과, 수목높이는 평균 높이가 0~1 m 범위인 지역에서 산사태취약성을 높이는 데 기여하였으며, 이는 과거 산사태 이력과 높은 연관성을 보였다. 평면곡률 표준편차는 약 1.7이상의 곡률의 변동성이 큰 산악지역에서 양의 SHAP value를 보여 산사태 취약성을 높이는 데 기여한 것으로분석되었다. 고도의 경우, 265 m에서 400 m의 중간 고도를 보이는 지역에서 산사태 취약성을 높이는 데 기여하였으며, 실제로 산사태가 발생한 위치 또한 주로 중간고도의 지역에 분포함을 확인하였다.

본 연구에서는 비선형 데이터 세트에 대한 일반화 성능이 뛰어난 앙상블 모델을 활용하여 산사태 취약성 분석을 수행하고, SHAP 기법을 적용하여 산사태 발생 영향인자의 기여도를 정량적으로 평가하였다. 특히, 산사태발생 가능성에 대한 각 인자의 영향을 공간적 관점에서분석함으로써 기존 연구의 한계를 보완하고자 하였다.이러한 연구는 산사태 취약성 평가 결과의 해석 가능성을높이고, 위험성이 높은 지역을 정밀하게 식별하는 데기여할 수 있다. 이를 바탕으로 신뢰성 높은 산사태 위험지도를 제작하고, 효과적인 예방 정책을 수립하는 데활용될 수 있다. 또한, 정량적 평가 결과를 반영하여 객관적인 의사결정을 지원할 수 있을 것으로 판단된다.

이 연구는 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원(RS-2023-00222563)을 받아 수행된 연구입니다.