企业官网建设流程全解析

简介：面向地学与城市规划领域科研人员及工程师的Python机器学习滑坡易发性制图项目代码包，基于历史滑坡点与地形、植被、降雨等多类环境因子训练随机森林等模型，用于生成易发性概率图并支撑灾害风险区划。包内共3个文件，含inscode工程入口与核心配置脚本、HTML交互式结果展示页面以及gitignore版本控制辅助文件，压缩包仅9KB，虽为轻量级代码骨架但目录职责清晰，便于直接阅读与二次开发。目前已有51人学习下载，代码覆盖数据准备、特征工程、模型训练与评估（ROC-AUC、精确率、召回率）、超参数调优及制图可视化全流程，适合快速搭建本区域滑坡易发性预测基线。模块化设计严格，包含统一地理框架、正负样本构建、多尺度特征提取、十折空间交叉验证及SHAP重要性解释，同时给出空间自相关泄露、特征重要性异常等问题的修正方案，并预留时序信息、深度学习与不确定性量化等进阶接口，可迁移至其他区域坡度、降雨、岩性等栅格数据。

1. 滑坡易发性制图：机器学习在区域地质灾害评估里到底怎么落地

滑坡易发性制图，本质是把一片山区按“容易发生滑坡的程度”划分成不同等级的空间底图，Pytho n机器学习在其中解决的是“怎么把多因子评价从打分表变成数据驱动的分类器”的问题。它直接服务于国土空间规划、道路选线、水电选址和地灾隐患排查，不是学术自嗨。我做过几个项目后最大的体会是：卡脖子的从来不是模型，而是样本怎么构建、因子怎么加工——这两步决定了结果上限，换什么算法都救不回来。这篇文章从数据准备讲到出图分级，代码可以直接照着改，适合正做地灾评估、区域规划或遥感制图方向，想用机器学习产出成果图的工程师。

2. 样本数据构建：滑坡点与负样本的采集策略直接决定模型上限

2.1 滑坡点整理的三个步骤：面转点、去重、限距

滑坡易发性建模里，正样本来源一般是野外调查记录、历史灾害台账或遥感解译，常见形态是“滑坡面”而不是“滑坡点”。直接把面转成无数个像元当正样本是第一个坑——滑坡面上的相邻像元高度相似，模型会在局部学到大量重复信息，AUC虚高但实际泛化能力很差。我处理多期解译数据时，习惯每个滑坡面只取一个代表点，再按研究区分辨率做去重。

import geopandas as gpd # 读取滑坡面数据 slides = gpd.read_file('slides.shp') # 每个面只取一个代表点（面内点，不一定是质心） centers = slides.geometry.representative_point() pos = gpd.GeoDataFrame(geometry=centers, crs=slides.crs) # 按30米像元索引去重，防止同一个像元内出现多个正样本 pos['grid_id'] = (pos.geometry.x // 30).astype(int).astype(str) + '_' + (pos.geometry.y // 30).astype(int).astype(str) pos = pos.drop_duplicates(subset='grid_id') print(f'正样本数量: {len(pos)}')

representative_point()和质心的区别在于：质心可能落在面外，尤其对弧形滑坡体，而representative_point()保证返回面内点。grid_id 的玩法需要和因子栅格分辨率对齐，这里是 30 米像元，如果换 10 米或 90 米数据，// 30要同步改。同一个像元里出现多个代表点，说明原始解译面有重叠，去重能有效降低局部样本密度带来的过拟合。

2.2 非滑坡样本怎么取：随机采样、环带采样与分层采样

负样本（非滑坡点）是滑坡易发性建模里最容易翻车的地方。最朴素的做法是在研究区内完全随机采样，但灾害点往往分布在特定岩性和坡度区间，完全随机会让负样本大面积落在平坦河谷区，模型学到的“滑坡”其实是“山区”，而不是真正的地形组合。我常用的三种策略没有绝对优劣，关键看研究区面积和滑坡分布形态。

环带采样是业内相对稳的做法，核心思想是避开“滑坡影响区”和“滑坡体本身”，只在其外一定范围内取负样本，让模型比较“滑坡区”和“邻近的非滑坡区”之间的因子差异，而不是比较山区和平原。

import numpy as np import geopandas as gpd from shapely.geometry import Point from shapely.ops import unary_union # 在某个验证过的研究区范围内生成候选点 x0, x1, y0, y1 = study_area.total_bounds rng = np.random.default_rng(42) candidates = [Point(rng.uniform(x0, x1), rng.uniform(y0, y1)) for _ in range(50000)] cand = gpd.GeoDataFrame(geometry=candidates, crs=study_area.crs) # 缓冲：内缓冲200米剔除滑坡影响区，外缓冲1000米限定采样范围 buffer_inner = unary_union(slides.geometry.buffer(200)) buffer_outer = unary_union(slides.geometry.buffer(1000)) ring = buffer_outer.difference(buffer_inner) # 取落在环形带内的候选点作为负样本 neg = gpd.sjoin(cand, gpd.GeoDataFrame(geometry=[ring], crs=slides.crs), how='inner', predicate='intersects') print(f'负样本候选数量: {len(neg)}')

200 米和 1000 米这两个参数不是固定的，取决于滑坡规模和因子栅格分辨率。滑坡体长几百米时 200 米内缓冲合理；换成高分辨率数据可以缩到 100 米。采样比例上我一般把正负样本控制在 1:1 到 1:5 之间，负样本太多会让模型偏向预测“非滑坡”，太少则训练不稳定。分层取样要再叠一层按坡度分位或岩性类型的筛选，避免环带中平坦区占比过高。

2.3 把样本映射到因子栅格：像元提取与行列索引

样本点建好后，下一步是把每个点的因子值从栅格里提出来。这一步看起来简单，但栅格坐标系、投影不一致、nodata 值处理都会在这里暴露。我习惯用 rasterio 的索引方式提取，而不是逐点读坐标，因为 rasterio 的src.index(x, y)本身会处理坐标转换和越界。

import rasterio import pandas as pd factor_paths = { 'slope': 'factor/slope.tif', 'elevation': 'factor/elevation.tif', 'dist_fault': 'factor/dist_fault.tif', 'ndvi': 'factor/ndvi.tif', } def extract_factors(gdf_points, factor_paths): rows = [] with rasterio.open(list(factor_paths.values())[0]) as src: transform = src.transform for geom in gdf_points.geometry: x, y = geom.x, geom.y row = {} for name, path in factor_paths.items(): with rasterio.open(path) as src: try: r, c = src.index(x, y) val = src.read(1)[r, c] row[name] = float(val) if val != src.nodata else np.nan except Exception: row[name] = np.nan rows.append(row) return pd.DataFrame(rows) X_pos = extract_factors(pos, factor_paths) X_neg = extract_factors(neg, factor_paths)

这段代码里最容易被忽略的是val != src.nodata判断。很多 DEM 和遥感产品在无数据区用 -9999 或 0 填充，如果不处理，模型会把无数据区当成真实的地形值，结果图上经常出现整片异常易发区。提取完因子矩阵后，先打印一下每列的缺失值比例，超过 5% 的样本考虑剔除而不是补值，因为山区边缘的缺失往往对应地形突变区，补值会掩盖真实边界。

3. 评价因子体系：选什么、怎么加工、怎么排查共线性

3.1 必选因子与不能碰的因子：一张因子清单说清楚

滑坡易发性因子没有标准答案，但有几类是公认必选的：地形类、水文距离类、地质类和植被类。我一般会把高程、坡度、坡向分量、曲率、距河流距离、距断层距离、距道路距离、岩性、NDVI 作为基础集合，再根据研究区是否有历史降雨资料决定要不要加降雨因子。每个因子的数据来源不同，加工方式也不同。

有一类因子要坚决避开：任何依赖滑坡点空间位置的因子都不能进模型。比如“距历史滑坡点距离”“历史滑坡密度”“滑坡缓冲区概率”，这类因子本质上是把标签的邻域信息喂给模型，预测时模型只会把已知滑坡周边一圈都标成高易发，图面看起来漂亮，但对无记录区毫无预测能力。判断标准很简单：如果因子栅格的制作过程需要用到滑坡点坐标，它就是泄漏源。

3.2 连续因子的加工细节：坡度、坡向分量与归一化

坡度和坡向从 DEM 生成，不需要自己写窗口计算，GDAL 的 DEMProcessing 足够稳。坡向是环形变量，0 度和 360 度是同一个方向，直接作为数值特征会让模型错误地认为“正北”和“正北偏西”差距巨大。我习惯在做因子时把坡向拆成南北分量和东西分量，保留方向信息同时避免环形断裂。

gdal.DEMProcessing slope.tif dem.tif -computeEdges gdal.DEMProcessing aspect.tif dem.tif -computeEdges gdal.DEMProcessing curvature.tif dem.tif

import numpy as np from osgeo import gdal aspect_ds = gdal.Open('aspect.tif') aspect = np.deg2rad(aspect_ds.ReadAsArray()) north = np.cos(aspect) east = np.sin(aspect) # 保存成与aspect同尺寸的tif profile = {'driver': 'GTiff', 'height': aspect.shape[0], 'width': aspect.shape[1], 'count': 1, 'dtype': 'float32', 'crs': aspect_ds.GetProjection(), 'transform': aspect_ds.GetGeoTransform()} with gdal.GetDriverByName('GTiff').Create('north_component.tif', aspect.shape[1], aspect.shape[0], 1, gdal.GDT_Float32) as ds: ds.SetGeoTransform(profile['transform']) ds.SetProjection(profile['crs']) ds.GetRasterBand(1).WriteArray(north)

坡向拆分量后，后续所有因子在训练前还要做一次标准化。随机森林和 XGBoost 对归一化不敏感，但如果要对比逻辑回归或做特征重要性解释，统一标准化是必要的。我用sklearn.preprocessing.StandardScaler对所有连续因子做标准化，注意拟合时只基于训练样本，避免测试集信息泄漏。距离类因子我习惯先取对数再标准化，因为距河流、断层距离往往长尾分布，直接标准化会被极少数远距离样本拉偏。

3.3 分类因子编码与VIF多重共线性筛查

岩性、土地利用这类分类因子直接塞进模型会有两个问题：一是类别值本身没有顺序语义，随机森林可能把“灰岩=5”和“砂岩=7”的错误比较利用起来；二是分类数过多时会增加过拟合风险。我常用的方案是岩性按频率比编码，用正样本占比除以研究区占比得到一个连续值；土地利用则用 one-hot。

共线性筛查是因子加工的最后一步。滑坡因子之间天然存在相关性，比如坡度和高程经常同向变化，曲率与坡度也有关联。多重共线性不会让随机森林崩溃，但会扭曲特征重要性排序，让你误判哪个因子才是主导。我用 VIF 做快速筛查，阈值取 10，超过的因子剔除或替换，保留业务解释性更强的那一个。

from statsmodels.stats.outliers_influence import variance_inflation_factor from statsmodels.tools.tools import add_constant # X_num 是只包含连续因子的特征矩阵 X_const = add_constant(X_num) vif_df = pd.DataFrame({ 'factor': X_const.columns, 'vif': [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])] }) vif_df = vif_df[vif_df['factor'] != 'const'].sort_values('vif', ascending=False) print(vif_df)

VIF 是线性回归的诊断指标，用在树模型上只能算参考，但它的价值在于暴露因子间的冗余关系。如果坡度 VIF 是 18、高程 VIF 是 14，这两个因子信息高度重叠，我一般保留坡度剔除高程，因为坡度是更直接的滑坡控制因子，解释性更强。VIF 筛查后还要看一眼相关性热力图，避免出现“坡向sin”和“坡向cos”之间意外的高相关。

4. 机器学习模型训练与评估：随机森林起步，XGBoost提精度

4.1 模型选型对比：逻辑回归、随机森林、XGBoost与SVM

滑坡易发性制图的模型选择，业内主流已经收敛到随机森林和 XGBoost。逻辑回归作为基线模型仍然有价值，它的输出概率有明确的物理含义，可以作为结果合理性对照；SVM 在样本量只有几百的小区域还行，一旦特征矩阵超过几万行，训练时间和内存开销都会失控。实际项目中我一般先用随机森林跑通全流程，出图后再对照 XGBoost 结果，两者分歧大的地方往往是值得检查数据质量的信号。

4.2 训练代码与必调参数：从特征矩阵到AUC

把第2章提取的正负样本因子矩阵合并，加上标签，就是标准表格数据。随机森林的训练代码没有太多花样，关键是几个参数怎么设。n_estimators我常用 300，再大收益微弱且拖慢预测速度；min_samples_leaf设到 5 左右，给叶子节点一点平滑，防止单个极端值影响；max_depth不限制，靠叶节点样本数来控制生长；class_weight根据正负样本比例决定要不要平衡。

import pandas as pd import numpy as np from sklearn.model_selection import train_test_split, StratifiedKFold from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import roc_auc_score, accuracy_score # X_pos 和 X_neg 来自第2章的特征提取，合并后打标签 X = pd.concat([X_pos, X_neg], axis=0).reset_index(drop=True) y = np.array([1] * len(X_pos) + [0] * len(X_neg)) # 分层划分：滑坡样本在训练/测试集中的比例保持一致 X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42, stratify=y) model = RandomForestClassifier( n_estimators=300, min_samples_leaf=5, max_features='sqrt', n_jobs=-1, random_state=42, class_weight='balanced' ) model.fit(X_train, y_train) y_prob = model.predict_proba(X_test)[:, 1] y_pred = (y_prob >= 0.5).astype(int) print('AUC:', round(roc_auc_score(y_test, y_prob), 4)) print('Accuracy:', round(accuracy_score(y_test, y_pred), 4))

class_weight='balanced'能自动按样本数量调整权重，但注意它不等于 SMOTE，不会改变样本分布，只是在损失函数里加重少数类的代价。正负样本 1:2 时我通常不加它；超过 1:5 才开。XGBoost 替换的话，核心参数是eta（学习率）、max_depth、subsample和colsample_bytree，eta从 0.05 起步，max_depth控制在 6 以内，subsample0.8 左右，用五折交叉验证挑，而不是直接跑默认值。

4.3 特征重要性与空间独立的验证划分

随机森林训练完后第一件事是看特征重要性排序，它既帮你验证因子加工是否正常，也暴露出是否混入了泄漏因子。重要性最高的是坡度、高程这类地形因子，通常是合理的；如果“距历史滑坡距离”排在首位，回头查因子清单。

importance = pd.Series(model.feature_importances_, index=X.columns) importance.sort_values(ascending=False).plot.bar(figsize=(8, 5))

评估指标只看测试集 AUC 远远不够，因为随机划分会让空间自相关样本同时出现在训练和测试集中。更严格的做法是按空间区块划分：把研究区按子流域或网格分成若干大块，一部分区块做训练，另一部分做验证，模型从没见过验证区块的邻域样本，AUC 才会反映真实泛化能力。随机森林在随机划分下 AUC 0.92，换空间独立划分后掉到 0.78，这是正常现象，后者才是值得写进报告的数字。

5. 避坑：滑坡易发性制图的四个高发翻车点

做滑坡易发性制图这几年，真正把结果推翻重来的几乎都是数据操作问题，而不是模型缺陷。下面四条是我见过并亲自踩过的高发坑，按“现象、原因、解决”记录，方便你对照排查。

5.1 滑坡面转点不彻底导致正样本重叠

现象：模型在训练区的准确率很高，AUC 0.95 以上，但出图后滑坡密集区出现整片高易发“团块”，范围远超实际滑坡边界，且单沟内部看不出强弱之分。

原因：滑坡面转栅格后，面内所有像元都被标记为正样本，几百米长的滑坡体变成几千个正样本，相邻像元特征几乎完全一致。模型学到的是“这个位置的相似地形都是滑坡”，而不是“具备什么组合条件才是滑坡”，本质是样本重度空间自相关。

解决：每个滑坡面只取一个代表点，再按像元尺寸去重。如果一个区域内滑坡面密集，相邻滑坡体之间的距离也小于最小像元尺寸，需要人为合并或剔除部分正样本。判断标准：正样本之间的最小欧氏距离不应小于研究区像元尺寸的 2 倍。

5.2 随机划分验证集导致AUC虚高

现象：随机划分测试集 AUC 0.93，换成按子流域划分后只有 0.76，模型效果看起来“缩水”了。出图后发现某条支沟的预测结果和训练区极其相似。

原因：滑坡在空间上聚簇分布，随机划分时同一滑坡聚簇的相邻样本一部分进训练集、一部分进测试集，模型在测试集上其实是靠记忆邻域完成的预测，而不是靠因子规律。这是空间数据评估中最常见的系统性高估。

解决：改用空间分块验证，按子流域或自然山脊线划分区块，保证验证区与训练区在地理上完全分离。至少报告一个空间验证 AUC，如果不忍心放弃随机划分结果，两者都写但明确标注。真正可信的数字是空间验证 AUC。

5.3 因子里混入“历史滑坡密度”造成数据泄漏

现象：特征重要性列表里“历史滑坡密度”排在第一，图面看高易发区就是历史滑坡点的平滑外扩，像做了缓冲区而非真正的易发性评价，给人“事后诸葛亮”的感觉。

原因：某些教程或工具箱会建议加“历史滑坡密度”或“距最近滑坡距离”因子，这类因子在建模阶段用了滑坡点位置生成，拟合阶段自然表现好，但换到没有灾害记录的区域就完全失效，违背了易发性制图“预判未发生区域”的初衷。

解决：因子清单建立后逐项审视“这个栅格是怎么生成的”——凡是需要输入滑坡点坐标的因子全部剔除。可以将“历史滑坡密度”作为检验变量用于结果验证，但绝不能进训练特征。

5.4 框选研究区边界不当导致负样本污染

现象：非滑坡样本在研究区边缘大量聚集，模型在边界上表现出“锯齿状”的易发性跳跃，出图后沿边界有一条不自然的条带。

原因：随机生成非滑坡候选点时未限制在研究区有效范围内，部分负样本落在斜坡阴影区或DEM边缘的无数据带；或者研究区范围划得过大，把大面积平原也纳入采样范围，模型被迫区分“山区”和“平原”。

解决：采样范围严格使用研究区的有效掩膜，掩膜从 DEM 的 nodata 区域和边界生成；负样本的坡度分布先做统计，与正样本坡度分布对比，如果两者重叠很少，说明采样范围偏了。

5.5 自然断点分级遇到“概率集中在一侧”的问题

现象：模型输出概率多数集中在 0.1~0.3，只有少量像元达到 0.8，分级时自然断点把大部分区域归为中易发，图面灰蒙蒙一片，分不出层次。

原因：这是类别不平衡和数据特性的正常结果，不是模型坏了。概率不等于易发性等级，分级的目的是突出相对风险差异，而不是均匀切分所有像元。

解决：用分位数分级代替自然断点，或用频率比方法重新标定等级边界；先叠加历史滑坡点看分布确定阈值，保证高易发区覆盖绝大多数历史灾害点。我习惯上直接落在这层经验：分级阈值让历史滑坡点在高和极高易发区占比超过 75%，再将低易发区压到最小面积。

6. 从概率到易发性分级图：制图输出的关键技巧

6.1 把模型输出还原成概率栅格

模型预测得到的是每个样本点的概率，要变成一张连续的易发性栅格图，需要把研究区所有像元的因子值构造成特征矩阵，批量喂给模型，再把输出概率写回栅格。研究区大时不要一次性把全部像元读进内存，按分块处理更稳。

import rasterio import numpy as np # factor_stack: 堆叠后的因子栅格数组 (n_features, height, width) # mask: 有效像元掩膜 factor_stack = np.stack([ rasterio.open('factor/slope.tif').read(1), rasterio.open('factor/elevation.tif').read(1), rasterio.open('factor/dist_fault.tif').read(1), ], axis=0) height, width = factor_stack.shape[1], factor_stack.shape[2] valid = ~np.any(np.isnan(factor_stack), axis=0) # 只对有效像元做预测 flat_matrix = factor_stack[:, valid].T # (N, n_features) prob_flat = model.predict_proba(flat_matrix)[:, 1] prob_raster = np.full((height, width), np.nan) prob_raster[valid] = prob_flat with rasterio.open('landslide_prob.tif', 'w', driver='GTiff', height=height, width=width, count=1, dtype='float32', crs=src.crs, transform=src.transform) as dst: dst.write(prob_raster, 1)

写回栅格时必须使用和因子栅格完全相同的投影与变换参数，否则结果图与原始底图对不齐，叠加时就会出现偏移。分块预测的写法是用行号切片分批读取和预测，避免大规模矩阵一次性占满内存。输出后第一件事是叠到山体阴影图上肉眼检查：高易发区应该沿山脊、沟谷和特定岩性条带分布，而不是出现沿轨道或云边缘的条带。

6.2 分级与独立验证：自然断点不是唯一选择

概率栅格不是最终成果，还要分级成低、中、高、极高四类。自然断点是默认选择，但滑坡易发性场景下我更推荐用历史滑坡点的分布来确定阈值，因为分级的目的是让最大比例的灾害点落入高易发区，而不是追求数理统计上的断点最优。实际操作是：先把概率栅格按分位数粗分 10 类，统计每类里历史滑坡点的数量，找到分布拐点，再据此合并成最终等级。

用第2章处理好的历史滑坡点叠加到分级图上，统计落在每个等级内的比例是最直观的验证。合格图件的标志通常是：高和极高易发区面积占比不超过 30%，却覆盖了 75% 以上的历史滑坡点。如果达标，叠加图也自然能说明问题；不达标就回头查采样范围或因子加工，而不是换一个分级方法掩盖问题。我有一次为了追求“好看”的分级效果反复调整阈值，最后发现不过是把概率集中在中间区间的问题，真正原因是负样本采样范围太大，重采后事情就顺了。现在每次出图我都保留分位数分级的结果做底稿，历史滑坡点叠上去先看分布，再定最终阈值——这比任何统计指标都更能说明图对不对。希望帮到你。

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