简介这套Python机器学习滑坡易发性制图项目代码面向地理信息、地质灾害防治及机器学习交叉领域的开发者与研究人员针对滑坡易发性空间预测建模全流程需求提供一套模块化、可复现的工程化实现尤其适合作为科研实验与业务应用的基线参考。压缩包共3个文件以inscode环境配置、html结果说明及工程辅助文件为主整体仅9KB轻量易用。代码覆盖数据准备与特征工程统一重采样30米、非滑坡样本生成、GLCM纹理特征提取、随机森林与LightGBM等模型训练、贝叶斯超参数调优、多维度指标评估及易发性分级制图并内置十折空间块交叉验证与成功率/预测率曲线避免空间自相关泄露提升模型外推能力同时提供SHAP可解释性分析与异常检测机制。已有45人学习下载适合希望快速搭建滑坡易发性预测基线流程、理解机器学习与GIS结合实战细节的读者。1. 项目背景与核心思路1.1 为什么要做“滑坡易发性制图”滑坡是一种常见的地质灾害尤其在山区、丘陵地带以及降雨充沛的区域每年都会造成不小的经济损失和人员伤亡。做滑坡易发性制图本质上就是回答一个问题在某个区域内哪些地方更容易发生滑坡但这个“容易”不是拍脑袋猜出来的而是要根据历史滑坡分布、地形地貌、地质条件、水文环境等一系列因素用数学模型推算出一个概率值最后画成一张空间分布图。传统做法一般是专家打分法、层次分析法、信息量法这些统计模型缺点是主观性强、非线性关系拟合能力弱。到了近些年机器学习方法在滑坡易发性领域已经成了主流——随机森林、支持向量机、逻辑回归、XGBoost这些算法可以直接从历史数据中学习滑坡与环境因子之间的复杂映射关系预测精度普遍优于传统模型。我自己最初接触这个方向时最大的困惑是网上资料很多但零散得很。有讲数据处理的有讲模型跑的但很少有把“从原始数据到最终成图”整条链路串起来的完整项目代码。所以这次我把自己实际跑通的一套流程整理出来从数据准备、特征工程、模型训练到空间制图每一步都写清楚原理和代码逻辑希望给正在做相关课题的同学或工程师一条可以直接“抄作业”的路线。1.2 方法选型为什么用机器学习先说我个人的结论在滑坡易发性这个场景下机器学习对比传统模型有三个核心优势。第一非线性拟合能力强。滑坡发生和环境因子之间绝不是简单的线性关系比如坡度在某个区间内越陡越危险但超过某个临界值后因为岩土类型不同规律又可能变化。随机森林这类树模型天生擅长处理这种非线性、交互效应比较强的数据。第二特征重要性可解释。这一点经常被忽略但在实际工程和论文里非常重要。机器学习模型可以输出每个因子的重要性排序比如告诉你“坡度贡献了30%的预测能力降雨量贡献了20%”这为后续的工程治理和区域规划提供了客观依据。第三泛化能力和可复现性。只要数据规范、流程清晰同一套代码可以迁移到不同研究区只需要换数据和参数不需要重新设计整个算法框架。这对于做区域推广来说是巨大的效率提升。当然机器学习也不是万能的。它非常依赖训练样本质量如果滑坡编目数据不准确或者影响因子数据精度太低模型再怎么调参也白搭。所以我的建议是数据准备花的时间应该占整个项目的一半以上而不是把重心全放在调模型上。2. 数据准备与特征工程2.1 滑坡编目数据怎么处理滑坡易发性建模的第一个输入是滑坡编目数据。说白了就是一份记录了“历史上哪些地方发生过滑坡”的空间数据通常是一个点图层或者面图层。我在实际项目里用的是面图层但面图层要先转成点才能和栅格数据对齐做样本提取。这里有几个容易踩坑的地方滑坡点不能直接用质心。很多人图省事直接把每个滑坡面的质心当作发生点。但在坡度大的区域一个滑坡面可能跨越很大的高程范围质心不一定落在“真正发生滑动”的位置。更合理的做法是用滑坡面的后缘点或者至少用质心再加一个随机偏移。训练样本和验证样本要空间分离。滑坡点往往有空间聚集性如果不做空间分区直接把所有点随机划分成训练集和验证集模型会严重过拟合——因为同一个滑坡区域的邻近点几乎共享相同的特征验证集里其实混入了训练集的“近亲”。我建议先把研究区按格网分成若干块一部分块做训练另一部分块做验证这才符合真实应用场景。负样本怎么选。滑坡易发性建模是二分类问题正样本是“发生过滑坡”的栅格负样本是“从未发生过滑坡”的栅格。负样本通常用随机采样从无滑坡区域提取数量一般和正样本按1:1到1:3之间平衡。我试过很多比例最后觉得1:1最常见、最好解释1:2在某些数据集上精度略高但也没到显著提升的程度。2.2 影响因子怎么选、怎么归一化影响因子是模型的输入特征也就是用来预测滑坡发生概率的自变量。常用的因子包括类别具体因子数据类型地形因子高程、坡度、坡向、曲率、地形湿度指数TWI连续型地质因子岩性、距断层距离类别/连续型水文因子距河流距离、降雨量连续型地表覆盖NDVI、土地利用类型连续/类别型人为因子距道路距离连续型我实际项目里一开始选了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_size0.3, random_state42, stratifyy )3.2 模型训练与超参数说明我在项目里重点对比了三个模型逻辑回归、随机森林、XGBoost。原因很简单这三个模型分别是“线性基准”“经典集成”“梯度提升”的代表横向对比既能看出线性模型的不足也能了解树模型在这个场景下的表现差异。核心代码逻辑如下from sklearn.linear_model import LogisticRegression from sklearn.ensemble import RandomForestClassifier from xgboost import XGBClassifier models { Logistic Regression: LogisticRegression(max_iter1000), Random Forest: RandomForestClassifier(n_estimators200, random_state42), XGBoost: XGBClassifier(n_estimators200, learning_rate0.1, max_depth6, random_state42) } 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%以上但这样的模型没有任何实际意义。我建议重点关注以下几个指标AUCROC曲线下面积这是最常用的指标表示模型把正样本排在负样本前面的概率。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), dtypebool) 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, dtypenp.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})这段代码里有几个细节值得注意NoData掩膜必须做。如果直接把NoData值给模型模型可能会把它当成一个真实数值来预测导致边缘区域出现奇怪的预测结果。一定要先用掩膜把无效像元挑掉只对有效像元预测。特征顺序必须和训练时一致。训练时你是按factor_paths字典的顺序提取特征的预测时也必须用完全相同的顺序。这个看似简单但非常容易出错——一旦字典顺序变了模型输入的特征就全对不上了而且模型不会报错只会默默给你一个完全错误的结果。内存管理。如果研究区很大整张栅格一次性读进内存可能会爆内存。我遇到过一次20000×20000的大图直接把内存吃满了。解决办法是分块读取每块1000×1000逐块预测再拼回去。4.2 结果验证与分级栅格预测完成后不能直接拿去用。我习惯用“成功率曲线”和“预测率曲线”来做验证这两个曲线和ROC曲线类似但更贴合空间分析的场景。核心区别是成功率曲线用的是训练集的滑坡点预测率曲线用的是验证集的滑坡点两者之间的差距可以反映模型的过拟合程度。计算逻辑很简单把预测概率栅格从高到低分成20个等级每个等级约5%的区域面积然后统计每个等级里实际滑坡点的占比画成累积曲线。如果曲线很快攀升到高位说明模型把高易发区识别得很准。实际操作里还可以把预测概率按自然间断点法Jenks或者等间隔法分成高、中、低、极低四个易发等级这样出图时更直观。5. 常见问题与实用经验5.1 样本不平衡问题滑坡易发性建模里负样本非滑坡点数量一般远大于正样本滑坡点如果直接硬训练模型会倾向于把一切预测为“非滑坡”。前面我提到过用随机采样控制负样本比例到1:1或1:2可以缓解这个问题。另一个思路是调整类别权重在sklearn里直接用class_weightbalanced就行RandomForestClassifier(n_estimators200, class_weightbalanced)我在实际对比中发现当正负样本比例在1:1时class_weight的影响不大当比例到1:4时class_weightbalanced能明显提升Recall代价是Precision略降。具体怎么取舍取决于你实际业务里更看重哪个指标。5.2 过拟合与特征共线性树模型在训练集上的表现通常非常“好看”接近100%的准确率但验证集一测就露馅。我判断过拟合有三个标准训练集AUC 0.98但验证集AUC掉到0.85以下差距过大。特征重要性集中在某个单一特征上比如高程占了60%以上其他特征几乎没贡献。成功率和预测率曲线差距超过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地表形变数据做动态易发性评价。希望这篇博客能帮你少走一些弯路有完整代码需求或者想深入交流的欢迎留言讨论。本文还有配套的精品资源点击获取