简介：面向地质灾害防治与空间数据分析人员，这份Python机器学习滑坡易发性制图项目代码包，聚焦历史滑坡点与环境因子结合下的随机森林模型训练与概率制图流程。压缩包仅9KB，共3个文件，以inscode配置、html展示页面和gitignore版本控制文件为主，简洁且便于快速查看与集成。已有45人浏览学习。代码完整覆盖统一地理框架搭建、非滑坡样本创建、特征提取与清洗、模型评估（ROC-AUC、精确率、召回率等）、超参数调优以及五级易发性分级制图；同时处理了空间自相关泄露、特征重要性异常等实际难题，并附带成功率与预测率曲线验证，用于检验模型空间外推能力。对正在搭建滑坡易发性评价流程的研究者或工程师而言，这是一份轻量、可复用且工程化程度较高的参考实现。

1. 项目概述与核心需求解析

1.1 滑坡易发性制图到底在解决什么问题

滑坡易发性制图，说白了就是回答一个问题：在某个山区里，未来哪些地方更容易发生滑坡？这不是算命，而是基于历史滑坡点和地形、地质、环境等影响因素，用数学模型去推演整个区域的“危险程度分布图”。

我刚开始接触这个方向时也纳闷，为什么非要绕一大圈用机器学习来做，传统方法不行吗？实际上传统方法（如层次分析法、信息量法）不是不行，而是它们依赖专家打分和线性叠加，主观性强、非线性关系捕捉能力弱。而机器学习方法，比如随机森林、XGBoost、支持向量机，能从历史数据里自动学习各因子与滑坡之间的复杂非线性关系，在多个研究区的对比实验里，AUC（曲线下面积）普遍比传统方法高出一截——这不是个别现象，是学术界的共识。

这个项目的核心价值在于，它是一套完整可跑的 Python 代码，从数据预处理、因子筛选、模型训练、精度评估到最终的易发性制图，链路完整。对刚刚入门机器学习、地质灾害或 GIS 分析方向的同学来说，这份代码比单纯看教材有用得多，因为它把“机器学习”和“空间分析”两个领域的知识真正串在了一起。我自己当年在跑通这套流程时，最大感受就是：原来课本上的随机森林、交叉验证，落到真实地理数据上时，会有这么多细节问题要处理。

1.2 项目的技术栈与适用人群

从技术栈上看，这个项目以 Python 为核心语言，依赖 Geopandas、Shapely、Rasterio、Scikit-learn、Matplotlib 等库，覆盖了空间数据处理、机器学习建模和结果可视化三大块。如果你对这些库不熟，没关系，后面我会逐个说清楚它们各自扮演什么角色。

适用人群方面，我觉得有三类朋友会特别受用：

• 正在做地质灾害相关毕业设计或论文的同学，尤其是地质工程、地理信息系统、测绘工程等专业
• 想学习“机器学习如何落地到实际行业数据”的 Python 开发者，这套代码是一个极好的范例
• 从事自然资源、应急管理、交通选线等领域的技术人员，需要快速产出区域危险性评价图

无论你属于哪一类，我都建议你先不要急着跑代码，而是花十分钟把整个技术流程的骨架理解清楚。接下来我就按我实际做项目时的思考和操作顺序，把每一步拆开讲透。

2. 数据准备与因子体系构建——易发性制图的地基

2.1 核心数据：滑坡编目与影响因子

任何机器学习项目，数据是第一道关口。滑坡易发性制图需要两类基础数据：一是历史滑坡点数据，也就是“正样本”；二是研究区的影响因子数据，也就是模型的输入特征。

滑坡编目数据通常来自野外调查、遥感解译或历史灾害记录，常见的格式是点状矢量数据（Shapefile）。这里的坑在于，很多公开数据集的滑坡点并没有经过严格核实，有些点可能是误判的，有些滑坡体很大但只记录了中心点。我在实际项目中会先做一步“数据清洗”，把明显落在河流中心、道路正中间或者建筑屋顶上的异常点剔除掉，再用高分辨率影像随机抽检一部分点的可靠性。看起来多花了一两天时间，但对模型精度的提升是决定性的。

影响因子方面，学术界没有统一标准，但最常用的因子集中在以下几类：

• 地形因子：高程、坡度、坡向、曲率、地形湿度指数（TWI）、地形起伏度
• 地质因子：岩性、距断层距离
• 水文因子：距河流距离、降雨量
• 人类活动因子：距道路距离、土地利用类型

这里有一个容易犯的错误：因子不是越多越好。如果两个因子之间相关性极高（比如高程和地形起伏度往往强相关），模型容易出现多重共线性问题，虽然树模型对此不敏感，但逻辑回归和 SVM 就会受到很大影响。我的做法是先算一下各因子之间的皮尔逊相关系数，把相关系数高于 0.8 的因子对做一个取舍，保留从物理意义上更合理的那个。

2.2 因子数据的准备与统一坐标

拿到原始数据后，最耗时但最不能跳过的步骤就是因子图层的制备。所有因子栅格必须满足三个条件：同一坐标系、同一范围、同一分辨率。

以我常用的研究区为例，30米分辨率的 DEM 数据是基础，坡度、坡向、曲率、TWI 都从 DEM 派生。距断层、河流、道路的距离，用欧氏距离工具计算。岩性和土地利用这类类别型因子，需要重分类编码成数值。如果这些因子层的坐标系不一致，后续的样本提取环节会直接报错或者提取到错误位置的值，这是新手最容易卡住的地方。

我这里强调一下分辨率的细节：30米分辨率不等于简单地统一像元大小。重采样时我建议用双线性插值处理连续型因子（如高程、坡度），而类别型因子（如岩性、土地利用）要用最近邻法，否则类别边界会出现插值产生的“伪类别”。这个小细节在很多教程里都没人提，但做出来图的差异肉眼可见。

2.3 样本数据的生成：正样本与负样本的博弈

样本制作是滑坡易发性制图里最具争议、也最影响模型上限的环节。正样本是已知的历史滑坡点，负样本呢？没有人能确切告诉你“哪些地方绝对不会滑坡”，所以负样本只能用策略来生成。

常见的负样本生成方法有几种：

• 在研究区内随机生成与正样本等量的点
• 在距滑坡点一定缓冲距离之外的区域生成
• 先做缓冲区排除，再结合坡度极低区域（如平原水体）生成负样本

我个人的实践经验是：正负样本比例 1:1 是多数研究的默认配置，但 1:1 结构在真实场景中可能偏低。实际做滑坡危险性分区时，负样本生成范围建议避开已有滑坡的缓冲区，这个缓冲半径一般取 100 到 200 米比较合理，太大会把潜在易滑区域误判成“安全区”，负样本里混入隐性的真实滑坡点，模型精度虚高。

样本量方面，几百个滑坡点的研究区比较常见。如果滑坡点太少（低于 50 个），机器学习模型很容易过拟合，此时优先选择逻辑回归这类简单模型，或者做 SMOTE 过采样。

3. 模型构建与样本集划分——机器学习核心环节

3.1 编程环境准备与库安装

跑这个项目，我建议直接用 Anaconda 创建独立环境，避免各个库之间的依赖冲突。特别说明一下，很多新手把大量时间耗费在配置环境上，这个项目相对轻量，依赖不算复杂，但你仍需要在 PyPI 源不可用时切换国内镜像源。

安装核心库的命令如下：

conda create -n landslide python=3.9
conda activate landslide
pip install geopandas rasterio shapely scikit-learn pandas numpy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple

解释一下每个库的作用：Geopandas 负责读写矢量数据（滑坡点），Rasterio 负责读写栅格数据（因子图层），Shapely 做几何运算，Scikit-learn 提供机器学习模型，Matplotlib 出图。这套组合是目前 Python 空间分析的主流配置，用熟了之后做其他地理类项目也能直接复用。

很多教程会推荐直接用 QGIS 或 ArcGIS 的 Python 环境，但这类 GIS 软件内置的 Python 环境版本一般比较旧，且包管理权限受限。我试过在 ArcGIS 自带的 Python 里安装最新版的 Scikit-learn，结果版本冲突直接导致软件崩溃，最后老老实实回到 Anaconda 单独建环境。这个坑你提前知道，能省掉半天折腾时间。

3.2 样本提取与数据集划分

因子准备齐了，接下来的操作逻辑是：把滑坡点坐标落到每个因子栅格上，提取对应的因子值，每个点最终得到一条形如 [坡度, 高程, 距断层距离, ..., 标签] 的数据记录。

这里的代码实现是：

import geopandas as gpd
import rasterio
import numpy as np
import pandas as pd
from rasterio.transform import rowcol

# 读取滑坡点数据（正样本）
landslides = gpd.read_file('landslide_points.shp')

# 读取因子栅格文件列表
factor_paths = ['slope.tif', 'elevation.tif', 'dist_fault.tif', 'twi.tif', 'curvature.tif']
factors = [rasterio.open(fp) for fp in factor_paths]

# 统一检查坐标系
for f in factors:
    assert f.crs == factors[0].crs, "坐标系不一致，请先统一"

# 提取每个滑坡点的因子值
rows, cols = rowcol(factors[0].transform, 
                    landslides.geometry.x.values, 
                    landslides.geometry.y.values)

sample_data = []
for i in range(len(landslides)):
    row, col = rows[i], cols[i]
    vals = [factor.read(1)[row, col] for factor in factors]
    sample_data.append(vals + [1])  # 1 表示正样本

# 生成负样本（随机点）
x_min, y_min, x_max, y_max = landslides.total_bounds
np.random.seed(42)
n_neg = len(landslides)
random_x = np.random.uniform(x_min, x_max, n_neg)
random_y = np.random.uniform(y_min, y_max, n_neg)

for i in range(n_neg):
    row, col = rowcol(factors[0].transform, random_x[i], random_y[i])
    vals = [factor.read(1)[row, col] for factor in factors]
    sample_data.append(vals + [0])  # 0 表示负样本

# 转 DataFrame
df = pd.DataFrame(sample_data, columns=['slope', 'elevation', 'dist_fault', 'twi', 'curvature', 'label'])
df = df.replace([np.inf, -np.inf], np.nan).dropna()

代码里有两处容易出错：一是负样本随机点可能落在研究区边界之外或者无数据区域（NoData），提取出来的值是 NaN，所以最后做了一步替换和剔除；二是 rowcol 索引出来的行和列直接用于 factor.read(1)[row, col] ，这个索引顺序和普通的二维数组完全一致，不用转换。

样本集划分通常采用 7:3 比例，但要注意划分方式。滑坡点往往在空间上存在聚集性，如果直接用随机划分，训练集和测试集可能包含距离极近的样本点，导致模型评估结果虚高。学术界更严格的做法是按空间分区划分（如按格网分区），但作为入门项目，随机划分已经足够，空间划分可以作为后续进阶方向。

3.3 机器学习模型的选择与训练

在滑坡易发性研究里，应用最广的模型是随机森林和 XGBoost。逻辑回归作为基准线，SVM 在很多研究里也有不俗表现，但它在样本量较大时训练速度慢，且对超参数敏感。

回归到代码上，Scikit-learn 实现随机森林训练非常简洁：

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, classification_report, confusion_matrix

X = df.drop('label', axis=1)
y = df['label']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 随机森林模型
rf = RandomForestClassifier(
    n_estimators=300,
    max_depth=10,
    min_samples_split=5,
    min_samples_leaf=2,
    random_state=42,
    n_jobs=-1
)
rf.fit(X_train, y_train)

# 预测与评估
y_pred = rf.predict(X_test)
y_proba = rf.predict_proba(X_test)[:, 1]

print("AUC: ", roc_auc_score(y_test, y_proba))
print(classification_report(y_test, y_pred))

核心要理解几个参数的物理意义： n_estimators 是决策树数量，太少模型偏差大，太多训练时间线性增加但精度提升趋于平缓，300 是一个比较稳妥的默认值； max_depth 限制树的深度，防止模型逐样本死记硬背，10 对于因子数量在 5 到 15 个之间的数据集是个合理起点； min_samples_leaf 控制叶子节点的最小样本数，它是我调参时最敏感的开关，建议先在 2 到 5 之间搜索。

超参数调优我推荐用网格搜索加交叉验证，初学者直接用以下代码即可：

from sklearn.model_selection import GridSearchCV

param_grid = {
    'n_estimators': [200, 300, 400],
    'max_depth': [5, 10, 15],
    'min_samples_leaf': [1, 2, 5]
}

grid = GridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid,
    cv=5,
    scoring='roc_auc',
    n_jobs=-1,
    verbose=1
)
grid.fit(X_train, y_train)

print("Best params:", grid.best_params_)
print("Best AUC:", grid.best_score_)

网格搜索是很直观的，它把这个参数表里的所有组合都训练了一遍。不过说实话，如果你只是完成毕业设计或做一个初步评价图，直接使用默认参数加少量手动调整已经能拿到不错的效果，网格搜索更多是让你心里有数，知道模型在当前数据上的上限在哪里。

3.4 模型精度评估与因子重要性分析

滑坡易发性模型评估，行业标配是 ROC 曲线和 AUC 值。AUC 的含义是“随机抽取一个正样本和一个负样本，模型给正样本打更高分的概率”，因此 AUC 越接近 1，模型区分危险区和安全区的能力越强。通常来说，AUC 在 0.75 以上属于可接受，0.85 以上属于优秀。

因子重要性则是随附在树模型上的一个优势。随机森林的 feature_importances_ 属性给出了每个因子在全部树节点分裂中带来的纯度增益总和，归一化后就是重要性占比。

import matplotlib.pyplot as plt

importance = rf.feature_importances_
features = X.columns

plt.figure(figsize=(10, 6))
plt.barh(features, importance, color='steelblue')
plt.xlabel("Importance")
plt.title("Random Forest Feature Importance")
plt.tight_layout()
plt.savefig('feature_importance.png', dpi=300)

从我的多个研究区经验来看，坡度和高程通常稳居重要性前列，距断层距离在某些构造活跃区占比很高，而距道路距离在不同区域差异巨大。这个排序结果本身就是有价值的科学产出，写论文时可以直接用。

4. 易发性制图与结果可视化——把模型输出变成一张图

4.1 全区域预测与栅格图生成

模型训练好了，测试集 AUC 也看了，但真正交付的成果并不是模型文件，而是一整张研究区的易发性分区图。这张图的生成逻辑是：把研究区内每一个像元对应的因子值提取出来，丢进模型得到滑坡发生概率，再把概率值写回栅格。

具体实现是不用像样本提取那样逐个枚举像元的，因为数据量太大（30米分辨率下，一个县可能几百万个像元，逐个循环 Python 会慢到怀疑人生）。更高效的做法是用 numpy 将各因子栅格读成数组，然后展开成二维表，用模型批量预测，再 reshape 回栅格矩阵。

import rasterio
import numpy as np

# 读取所有因子栅格
factor_arrays = []
for fp in factor_paths:
    with rasterio.open(fp) as src:
        factor_arrays.append(src.read(1))

# 堆叠成三维数组
stacked = np.stack(factor_arrays, axis=-1)  # 形状 (height, width, n_factors)
h, w, _ = stacked.shape
flat = stacked.reshape(-1, len(factor_paths))

# 预测概率
proba = rf.predict_proba(flat)[:, 1]
proba_map = proba.reshape(h, w)

# 写为 GeoTIFF
with rasterio.open(factor_paths[0]) as ref:
    profile = ref.profile
    profile.update(dtype=rasterio.float32, count=1)

with rasterio.open('landslide_probability.tif', 'w', **profile) as dst:
    dst.write(proba_map.astype(rasterio.float32), 1)

这里有两个细节很容易出问题。第一，因子栅格堆叠时，如果各因子的边界或 NoData 掩膜不完全一致，把 NoData 值直接输入模型会导致推理结果完全没有意义。我习惯先把所有因子统一填 NoData 为 0 或 -9999，再用掩膜把无效区域覆盖率高的像元单独标记，预测后重新掩掉。第二，reshape 成二维表时，行列顺序必须和栅格数组的行列顺序完全一致，因为后续 reshape 回二维矩阵是按行优先展开，顺序错一位整张图就是错位的，这个问题肉眼很难察觉，需要和专业底图叠加检验。

4.2 易发性等级划分的阈值选择

连续的概率栅格图已经能看，但业务上更常用的是把概率划分为极高易发、高易发、中易发、低易发、极低易发五个等级。等级划分主要有两种方法：

• 等距划分：直接将 0 到 1 的概率区间五等分，简单直观，但往往导致各等级面积严重失衡
• 分位数划分：按概率值排序后等面积划分，每个等级面积基本相同，更便于管理应用

我强烈建议使用分位数法。自然灾害管理视角下，不同等级面积相近时，管理部门更容易制定差异化的防治措施。分位数划分实现极为简单：

percentiles = np.percentile(proba_map, [20, 40, 60, 80])
levels = np.zeros_like(proba_map, dtype=np.int32)
for i in range(4):
    levels[proba_map > percentiles[i]] = i + 1

with rasterio.open('landslide_susceptibility.tif', 'w', **profile) as dst:
    dst.write(levels.astype(np.int32), 1)

写到这里顺便说一句，不同等级的面积占比结果，在论文里通常要画成柱状图或饼图，同时按等级统计历史滑坡点落入数量，计算各等级滑坡频率比。这个频率比如果呈现出随等级升高而单调递增的趋势，说明模型分区结果符合认知，审稿人会很喜欢看到这种结果。

4.3 出图规范与美化技巧

出图是一张易发性图能不能被期刊接收、被业务部门采用的关键。学术出版领域主流的出图方式是将结果在 ArcGIS 或 QGIS 中完成符号化与制图排版，但 Python 也有快速出图方案。

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

# 自定义等级配色
colors = ['#1a9850', '#fee08b', '#fdae61', '#f46d43', '#d73027']
cmap = ListedColormap(colors)

fig, ax = plt.subplots(figsize=(12, 10))
im = ax.imshow(levels, cmap=cmap, extent=[x_min, x_max, y_min, y_max])
cbar = plt.colorbar(im, ax=ax, ticks=[0, 1, 2, 3, 4])
cbar.ax.set_yticklabels(['极低', '低', '中', '高', '极高'])
ax.set_title("Landslide Susceptibility Mapping based on Random Forest")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.tight_layout()
plt.savefig('final_map.png', dpi=300)

关于配色，有一点经验分享：绿到红渐变色是灾害类成图的主流选择，也被多数读者直觉理解，但一定要避开红绿色盲难以区分的配色组合。另外，叠加一些水系、道路、居民点等基础地理信息会让图面更专业，但等值线不要太密集，抢了分区色块的视觉焦点。

5. 常见问题与排查技巧实录

5.1 坐标系统的混乱与修复

这是新手绕不开的第一个坎。你也许会碰上滑坡点文件是地理坐标系（WGS84 经纬度），而 DEM 是投影坐标系（UTM 或高斯克吕格），此时读取点坐标提取栅格值时会得到完全离谱的结果，或者干脆超出栅格范围报错。

排查思路很简单：每次读取数据都打印 CRS 信息，并确认目标因子栅格的 CRS，再看滑坡点 CRS 是否一致。如果不一致就用 to_crs() 转换：

landslides = landslides.to_crs(factors[0].crs)

多花两分钟在这个检查上，后面会省下一堆乱七八糟的问题。

5.2 负样本生成不合理导致精度虚高

我在最初做实验的时候，曾经把负样本限制在坡度极低的河谷平原，结果模型的 AUC 直接冲到 0.95 以上。那时候很高兴，但仔细一分析才意识到这是“作弊”了——低坡度区域和高坡度区域的滑坡特征差异太明显，模型根本不需要学习复杂的因子交互就能将两者分开。这种精度虚高并没有实际指导意义，因为管理部门的决策恰恰需要区分的是坡度为 20 到 40 度之间的中高危险区域。

更合理的负样本策略是在整个研究区随机生成，并剔除掉距已有滑坡点 200 米以内的点。这样一来，模型学习的重点是“在相似地形条件下，哪些因子组合更容易触发滑坡”，结果更具可解释性，也不容易被审稿人质疑。

5.3 内存不足与特征矩阵过大

如果研究区范围很大，比如几个县连在一起，30 米分辨率的因子栅格全部读入内存后，几百兆甚至几个 G 的内存占用是很常见的。如果你的机器内存只有 8G，建议把全区域预测改成区域分块预测：

block_size = 512

for i in range(0, h, block_size):
    for j in range(0, w, block_size):
        block = stacked[i:i+block_size, j:j+block_size]
        proba_block = rf.predict_proba(block.reshape(-1, len(factor_paths)))[:, 1]
        proba_map[i:i+block_size, j:j+block_size] = proba_block.reshape(block.shape[:2])

这种分块计算是一招鲜吃遍天的思路。虽然比一次性算慢一点，但换来了程序不崩的稳定性。我在处理过一个大范围区域时，内存从峰值 6G 降到了不到 2G，再配合数据降采样，一般电脑都能稳稳跑完。

5.4 代码可重复性与工程化整理

最后聊聊项目代码怎么上传到 Git 仓库的问题。很多人一开始把代码和数据全放在同一个目录里，传到仓库后又发现数据文件太大推不上去。我的习惯是建立清晰的项目目录结构：

landslide-susceptibility/
├── data/          # 原始数据（用 .gitignore 忽略大文件）
├── factors/       # 处理后的因子栅格
├── scripts/       # 全部 Python 脚本
│   ├── 01_preprocess.py
│   ├── 02_sample_extract.py
│   ├── 03_model_train.py
│   └── 04_mapping.py
├── output/        # 模型输出、预测栅格、图表
└── README.md

数据文件太大时，可以使用 Git LFS 或把数据放到网盘并在 README 里提供下载链接，只把代码和配置文件传到仓库。这样既保证代码可复现，又不会让仓库膨胀。

6. 项目扩展与经验总结

跑通这套流程之后，整个滑坡易发性制图的核心能力就算掌握了，但距离“做出有价值的研究”还有一段路。我觉得有几个扩展方向很值得尝试：

第一，引入更多影响因子，比如降雨量、地震动峰值加速度、植被覆盖度（NDVI）、土壤类型等，因子数量从 5 个扩到 12 个以上，模型的表达能力会明显提升，但前期的数据收集成本也会成倍增加。第二，做模型对比实验，把随机森林和 XGBoost、SVM、逻辑回归放在同一数据集上比较，用 DeLong 检验判断 AUC 差异是否显著，这会让论文的学术性上一个档次。第三，尝试深度学习方法，比如卷积神经网络直接以多因子栅格作为输入，在空间特征提取上潜力很大，但需要更多数据，对显卡要求也更高。

从工程角度来说，把训练好的模型用 Joblib 保存下来，再用 Flask 或 FastAPI 封装一个推理接口，就能做成一个可以在线的滑坡易发性评估系统。这样研究不仅能发论文，还能落地到业务系统里，投入产出比更高。

说说我对这个项目的体会：滑坡易发性制图技术上完全不复杂，真正难的是对数据的理解和对每一个细节的把控。同样的代码，换了数据源和参数策略，结果可能天壤之别。所以建议每个跑通流程的朋友，都把数据可视化出来多看看，而不是只看几个精度指标——模型只能告诉你概率，但为什么这个区域概率高、那个区域概率低，需要你对研究区的实际情况有判断力，这恰恰是机器学习无法替代的部分。

本文还有配套的精品资源，点击获取