企业官网建设流程全解析

简介：这套Python机器学习滑坡易发性制图项目代码，面向地理信息、地质灾害防治及机器学习交叉领域的开发者与研究人员，针对滑坡易发性空间预测建模全流程需求，提供一套模块化、可复现的工程化实现，尤其适合作为科研实验与业务应用的基线参考。压缩包共3个文件，以inscode环境配置、html结果说明及工程辅助文件为主，整体仅9KB，轻量易用。代码覆盖数据准备与特征工程（统一重采样30米、非滑坡样本生成、GLCM纹理特征提取）、随机森林与LightGBM等模型训练、贝叶斯超参数调优、多维度指标评估及易发性分级制图，并内置十折空间块交叉验证与成功率/预测率曲线，避免空间自相关泄露，提升模型外推能力；同时提供SHAP可解释性分析与异常检测机制。已有45人学习下载，适合希望快速搭建滑坡易发性预测基线流程、理解机器学习与GIS结合实战细节的读者。

1. 项目背景与核心思路

1.1 为什么要做“滑坡易发性制图”

滑坡是一种常见的地质灾害，尤其在山区、丘陵地带以及降雨充沛的区域，每年都会造成不小的经济损失和人员伤亡。做滑坡易发性制图，本质上就是回答一个问题：在某个区域内，哪些地方更容易发生滑坡？但这个“容易”不是拍脑袋猜出来的，而是要根据历史滑坡分布、地形地貌、地质条件、水文环境等一系列因素，用数学模型推算出一个概率值，最后画成一张空间分布图。

传统做法一般是专家打分法、层次分析法、信息量法这些统计模型，缺点是主观性强、非线性关系拟合能力弱。到了近些年，机器学习方法在滑坡易发性领域已经成了主流——随机森林、支持向量机、逻辑回归、XGBoost这些算法可以直接从历史数据中学习滑坡与环境因子之间的复杂映射关系，预测精度普遍优于传统模型。

我自己最初接触这个方向时，最大的困惑是：网上资料很多，但零散得很。有讲数据处理的，有讲模型跑的，但很少有把“从原始数据到最终成图”整条链路串起来的完整项目代码。所以这次我把自己实际跑通的一套流程整理出来，从数据准备、特征工程、模型训练到空间制图，每一步都写清楚原理和代码逻辑，希望给正在做相关课题的同学或工程师一条可以直接“抄作业”的路线。

1.2 方法选型：为什么用机器学习

先说我个人的结论：在滑坡易发性这个场景下，机器学习对比传统模型有三个核心优势。

第一，非线性拟合能力强。滑坡发生和环境因子之间绝不是简单的线性关系，比如坡度在某个区间内越陡越危险，但超过某个临界值后，因为岩土类型不同，规律又可能变化。随机森林这类树模型天生擅长处理这种非线性、交互效应比较强的数据。

第二，特征重要性可解释。这一点经常被忽略，但在实际工程和论文里非常重要。机器学习模型可以输出每个因子的重要性排序，比如告诉你“坡度贡献了30%的预测能力，降雨量贡献了20%”，这为后续的工程治理和区域规划提供了客观依据。

第三，泛化能力和可复现性。只要数据规范、流程清晰，同一套代码可以迁移到不同研究区，只需要换数据和参数，不需要重新设计整个算法框架。这对于做区域推广来说是巨大的效率提升。

当然，机器学习也不是万能的。它非常依赖训练样本质量，如果滑坡编目数据不准确，或者影响因子数据精度太低，模型再怎么调参也白搭。所以我的建议是：数据准备花的时间应该占整个项目的一半以上，而不是把重心全放在调模型上。

2. 数据准备与特征工程

2.1 滑坡编目数据怎么处理

滑坡易发性建模的第一个输入，是滑坡编目数据。说白了，就是一份记录了“历史上哪些地方发生过滑坡”的空间数据，通常是一个点图层或者面图层。

我在实际项目里用的是面图层，但面图层要先转成点，才能和栅格数据对齐做样本提取。这里有几个容易踩坑的地方：

1. 滑坡点不能直接用质心。很多人图省事，直接把每个滑坡面的质心当作发生点。但在坡度大的区域，一个滑坡面可能跨越很大的高程范围，质心不一定落在“真正发生滑动”的位置。更合理的做法是用滑坡面的后缘点，或者至少用质心再加一个随机偏移。

2. 训练样本和验证样本要空间分离。滑坡点往往有空间聚集性，如果不做空间分区，直接把所有点随机划分成训练集和验证集，模型会严重过拟合——因为同一个滑坡区域的邻近点几乎共享相同的特征，验证集里其实混入了训练集的“近亲”。我建议先把研究区按格网分成若干块，一部分块做训练，另一部分块做验证，这才符合真实应用场景。

3. 负样本怎么选。滑坡易发性建模是二分类问题，正样本是“发生过滑坡”的栅格，负样本是“从未发生过滑坡”的栅格。负样本通常用随机采样从无滑坡区域提取，数量一般和正样本按1:1到1:3之间平衡。我试过很多比例，最后觉得1:1最常见、最好解释，1:2在某些数据集上精度略高，但也没到显著提升的程度。

2.2 影响因子怎么选、怎么归一化

影响因子是模型的输入特征，也就是用来预测滑坡发生概率的自变量。常用的因子包括：

我实际项目里一开始选了16个因子，但相关性分析后发现“平面曲率”和“剖面曲率”相关系数高达0.85以上，属于高度共线性，最后去掉了平面曲率，保留了剖面曲率。共线性问题在线性模型里影响很大，而在随机森林这类树模型里影响相对较小，但该做的检查还是要做，不要心存侥幸。

归一化方面，逻辑回归和SVM这类基于距离的模型必须做归一化，否则数值范围大的特征（比如高程，单位是几百米）会主导模型训练。随机森林和XGBoost是树模型，不做归一化影响不大。为了统一流程，我在代码里还是保留了标准化操作，用StandardScaler，这样以后切换模型时不用回头补这一步。

from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test)

这里要注意一个细节：只能用训练数据fit，然后用同一个scaler去transform训练集和测试集。如果先对全部数据做标准化再划分训练集和测试集，会引入数据泄露，模型评估结果会偏乐观。

3. 机器学习建模与训练

3.1 数据集构建与划分

数据准备完成后，接下来就是构建机器学习的输入数据集。这一步的逻辑是：把滑坡点和非滑坡点的坐标提取出来，然后用ArcGIS或者Python的rasterio库，从每个影响因子栅格中提取对应位置的像元值，形成一个二维表格。

我用的代码逻辑大致如下：

import rasterio import numpy as np import pandas as pd # 所有影响因子的栅格文件路径 factor_paths = { 'elevation': 'data/elevation.tif', 'slope': 'data/slope.tif', 'aspect': 'data/aspect.tif', 'curvature': 'data/curvature.tif', 'twi': 'data/twi.tif', 'dist_fault': 'data/dist_fault.tif', 'dist_river': 'data/dist_river.tif', 'rainfall': 'data/rainfall.tif', 'ndvi': 'data/ndvi.tif', 'lithology': 'data/lithology.tif', 'dist_road': 'data/dist_road.tif', } def extract_values(coords, factor_paths): """从每个栅格中提取坐标点的像元值""" features = {} with rasterio.open('data/landslide_points.shp') as src: for name, path in factor_paths.items(): with rasterio.open(path) as factor_src: values = [v[0] for v in factor_src.sample(coords)] features[name] = values return pd.DataFrame(features)

这里有一个性能方面的经验可以分享：使用rasterio.sample()按坐标点批量提取值，效率远远高于逐点用read()读取整个栅格再取行列号。如果样本量几万个点，逐点读取会慢到怀疑人生，sample()方法可以直接传入坐标列表，底层是批量操作的。

数据划分这块，我强烈建议用分层抽样。因为正样本和负样本数量不一样，如果不做分层，很可能训练集里全是正样本或者负样本，模型学不到有效特征。train_test_split里有个参数stratify，传y进去就行。

from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.3, random_state=42, stratify=y )

3.2 模型训练与超参数说明

我在项目里重点对比了三个模型：逻辑回归、随机森林、XGBoost。原因很简单，这三个模型分别是“线性基准”“经典集成”“梯度提升”的代表，横向对比既能看出线性模型的不足，也能了解树模型在这个场景下的表现差异。

核心代码逻辑如下：

from sklearn.linear_model import LogisticRegression from sklearn.ensemble import RandomForestClassifier from xgboost import XGBClassifier models = { 'Logistic Regression': LogisticRegression(max_iter=1000), 'Random Forest': RandomForestClassifier(n_estimators=200, random_state=42), 'XGBoost': XGBClassifier(n_estimators=200, learning_rate=0.1, max_depth=6, random_state=42) } for name, model in models.items(): model.fit(X_train_scaled, y_train) train_acc = model.score(X_train_scaled, y_train) test_acc = model.score(X_test_scaled, y_test) print(f"{name}: Train Acc = {train_acc:.4f}, Test Acc = {test_acc:.4f}")

随机森林的几个关键参数，我简单说明一下：

• n_estimators：树的数量。太少了模型欠拟合，太多了计算量增加但收益变小。实际测试中200棵树已经足够，500棵树精度提升不到0.5个百分点，但训练时间翻倍。
• max_depth：树的最大深度。不限制的话树会一直生长，容易过拟合。我一般设定在10-20之间，配合min_samples_split来防止过拟合。
• random_state：随机种子。必须固定，否则每次运行结果都不一样，论文无法复现。

XGBoost的几个超参数也有讲究：

• learning_rate：学习率，一般0.01-0.3之间。学习率越小，模型越稳定，但需要的树更多。
• max_depth：默认6，对于表格数据通常够用了。
• subsample：行采样比例，0.8左右可以增加随机性，防止过拟合。

3.3 模型评估指标怎么看

滑坡易发性建模的任务是二分类，但这里的“正样本”占比通常很低（滑坡区域在整个研究区里面积占比一般不到10%），所以不能只看准确率（Accuracy）。如果模型把所有样本都预测为“非滑坡”，准确率可能有90%以上，但这样的模型没有任何实际意义。

我建议重点关注以下几个指标：

AUC（ROC曲线下面积）：这是最常用的指标，表示模型把正样本排在负样本前面的概率。0.5是随机水平，0.7-0.8是可用，0.8-0.9是优秀，0.9以上可能存在过拟合嫌疑。我用测试集AUC来横向比较模型，通常是XGBoost > 随机森林 > 逻辑回归。

混淆矩阵和F1-score：尤其是F1-score，它同时考虑了精确率（Precision）和召回率（Recall）。在滑坡易发性场景下，我们更关心的是“真正会滑坡的地方有没有被漏掉”，即Recall要尽量高，但同时Precision也不能太低，否则预警范围会大到失去意义。

from sklearn.metrics import roc_auc_score, classification_report, confusion_matrix for name, model in models.items(): y_pred = model.predict(X_test_scaled) y_proba = model.predict_proba(X_test_scaled)[:, 1] auc = roc_auc_score(y_test, y_proba) print(f"\n{name}: AUC = {auc:.4f}") print(classification_report(y_test, y_pred)) print(confusion_matrix(y_test, y_pred))

从我的实际结果来看，XGBoost的AUC大概0.91，随机森林0.89，逻辑回归0.83。但要注意，AUC高并不代表模型输出的概率值在绝对意义上准确，它只是排序能力的度量。如果后续要做风险分级（比如划分高、中、低易发区），可能还需要用校准方法（如Platt Scaling）对概率值做一步校正，否则分级阈值可能不合理。

4. 易发性制图与结果验证

4.1 栅格预测与制图流程

模型训练好之后，最终的目标是把整张研究区栅格输入模型，得到每个像元的滑坡发生概率，然后输出成一张连续的概率栅格图。这个流程说白了就是：把每个栅格像元对应的特征值提取出来，组成一个二维特征矩阵，喂给模型，得到预测概率，再写回栅格文件。

这里有个关键点：不能用循环遍历全部像元。一张1000×1000的栅格就有100万个像元，循环逐个预测会非常慢。正确做法是把栅格展平成一个大数组，一次性批量预测。

import numpy as np from osgeo import gdal def predict_raster(model, scaler, factor_paths, output_path): # 读取一个因子栅格获取投影信息和尺寸 first_src = gdal.Open(factor_paths[list(factor_paths.keys())[0]]) geotransform = first_src.GetGeoTransform() projection = first_src.GetProjection() cols = first_src.RasterXSize rows = first_src.RasterYSize # 批量读取所有因子的全图数据 factor_ds = [] for path in factor_paths.values(): ds = gdal.Open(path) band = ds.GetRasterBand(1).ReadAsArray().astype(np.float32) factor_ds.append(band) # 掩膜掉所有因子都是NoData的像元 valid_mask = np.ones((rows, cols), dtype=bool) for band in factor_ds: valid_mask &= ~np.isnan(band) & (band != -9999) # 提取有效像元的特征矩阵 valid_indices = np.where(valid_mask) X_all = np.column_stack([band[valid_mask] for band in factor_ds]) # 标准化并批量预测 X_all_scaled = scaler.transform(X_all) pred_proba = model.predict_proba(X_all_scaled)[:, 1] # 写回栅格 result = np.full((rows, cols), np.nan, dtype=np.float32) result[valid_mask] = pred_proba driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(geotransform) out_ds.SetProjection(projection) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(result) out_band.SetNoDataValue(np.nan) out_ds.FlushCache() print(f"预测栅格已保存: {output_path}")

这段代码里有几个细节值得注意：

1. NoData掩膜必须做。如果直接把NoData值给模型，模型可能会把它当成一个真实数值来预测，导致边缘区域出现奇怪的预测结果。一定要先用掩膜把无效像元挑掉，只对有效像元预测。

2. 特征顺序必须和训练时一致。训练时你是按factor_paths字典的顺序提取特征的，预测时也必须用完全相同的顺序。这个看似简单，但非常容易出错——一旦字典顺序变了，模型输入的特征就全对不上了，而且模型不会报错，只会默默给你一个完全错误的结果。

3. 内存管理。如果研究区很大，整张栅格一次性读进内存可能会爆内存。我遇到过一次20000×20000的大图，直接把内存吃满了。解决办法是分块读取，每块1000×1000，逐块预测再拼回去。

4.2 结果验证与分级

栅格预测完成后，不能直接拿去用。我习惯用“成功率曲线”和“预测率曲线”来做验证，这两个曲线和ROC曲线类似，但更贴合空间分析的场景。核心区别是：成功率曲线用的是训练集的滑坡点，预测率曲线用的是验证集的滑坡点，两者之间的差距可以反映模型的过拟合程度。

计算逻辑很简单：把预测概率栅格从高到低分成20个等级（每个等级约5%的区域面积），然后统计每个等级里实际滑坡点的占比，画成累积曲线。如果曲线很快攀升到高位，说明模型把高易发区识别得很准。

实际操作里，还可以把预测概率按自然间断点法（Jenks）或者等间隔法分成高、中、低、极低四个易发等级，这样出图时更直观。

5. 常见问题与实用经验

5.1 样本不平衡问题

滑坡易发性建模里，负样本（非滑坡点）数量一般远大于正样本（滑坡点），如果直接硬训练，模型会倾向于把一切预测为“非滑坡”。前面我提到过，用随机采样控制负样本比例到1:1或1:2可以缓解这个问题。另一个思路是调整类别权重，在sklearn里直接用class_weight='balanced'就行：

RandomForestClassifier(n_estimators=200, class_weight='balanced')

我在实际对比中发现：当正负样本比例在1:1时，class_weight的影响不大；当比例到1:4时，class_weight='balanced'能明显提升Recall，代价是Precision略降。具体怎么取舍，取决于你实际业务里更看重哪个指标。

5.2 过拟合与特征共线性

树模型在训练集上的表现通常非常“好看”（接近100%的准确率），但验证集一测就露馅。我判断过拟合有三个标准：

1. 训练集AUC > 0.98，但验证集AUC掉到0.85以下，差距过大。
2. 特征重要性集中在某个单一特征上（比如高程占了60%以上），其他特征几乎没贡献。
3. 成功率和预测率曲线差距超过15%。

缓解过拟合的方法，我和大家分享三个实际管用的：

• 增加数据量是首选，尤其是补充负样本的多样性。
• 限制树深度（max_depth）、增加min_samples_split、增加min_samples_leaf。
• 用GridSearchCV或Optuna做超参数搜索，但要注意用交叉验证而不是单纯切一次训练集/测试集。

5.3 代码结构建议

做这类项目，代码组织一定要清晰，否则后期改参数、换模型时会一团糟。我最终的项目结构是这样的：

landslide_susceptibility/ ├── data/ │ ├── landslide_points.shp │ ├── factor_rasters/ │ └── processed/ ├── src/ │ ├── data_preprocessing.py │ ├── train_model.py │ ├── predict_raster.py │ └── utils.py ├── output/ │ ├── models/ │ └── maps/ ├── requirements.txt └── README.md

每个脚本只负责一个环节，数据和输出路径用配置文件统一管理，模型参数也集中写在config.py里。这样即使隔了几个月再回来，也能快速上手继续做。

还有一个经验：尽量用joblib.dump保存模型和标准化器，不要用pickle。joblib对numpy数组的序列化效率更高，加载也更快。

import joblib joblib.dump(model, 'output/models/rf_model.pkl') joblib.dump(scaler, 'output/models/scaler.pkl')

加载的时候直接：

model = joblib.load('output/models/rf_model.pkl') scaler = joblib.load('output/models/scaler.pkl')

这样后续做预测时就不用重新训练了，尤其是多个区域复用同一套模型参数时，效率提升非常明显。

6. 写在最后的几点体会

这几轮做下来，我最大的感触是：滑坡易发性制图真正难的从来不是“跑通代码”，而是把每个环节的细节想清楚——数据有没有空间自相关、负样本怎么采、特征怎么对齐、验证怎么做。很多初学者拿到公开数据集直接一通fit下去，AUC看着很高，但放到实际区域一验证就露馅。所以我特别建议，做这个方向的读者在完成整套流程后，一定要回到“数据是否合理”这个本源问题上多问几遍。

另外，如果你打算把这套方法用到自己的研究区，我建议前期花一到两周的时间把多源数据统一到同一坐标系、同一分辨率、同一范围。这个步骤看着不起眼，但后续所有环节的顺畅度都取决于它。数据统一好了，整个管线就像流水线一样顺畅；数据不统一，后面每一步都在和坐标系、像元大小较劲。这也是为什么业内经常说“垃圾进、垃圾出”的原因。

如果后续资料和精力允许，我还会继续分享滑坡易发性制图相关的内容，比如怎么把深度学习模型（如CNN或Transformer）引入这个场景，以及怎么结合InSAR地表形变数据做动态易发性评价。希望这篇博客能帮你少走一些弯路，有完整代码需求或者想深入交流的，欢迎留言讨论。

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