简介一份演示机器学习在成矿预测中应用的轻量级源码包面向地质科研人员、地学数据工程师及相关专业学生系统梳理从证据权重法到随机森林的建模路径重点展示机器学习如何突破传统方法条件独立性假设带来的局限。资源包内共3个文件核心为html页面与inscode项目配置另附gitignore文件整包约9KB便于查看源码结构与页面实现。目前已有59人学习/浏览。读者可借助该页面快速对比传统统计方法与机器学习在地学大数据场景下的建模差异理解随机森林如何捕捉变量间非线性关系与交互作用也能以此为参考骨架进一步扩展高维地质数据预测实验尤其适合初学者快速把握地学大数据驱动建模的思路与落地实现。1. 这份源码包能跑通什么成矿预测不是玄学手里有已知矿点的 shp、1:5 万水系沉积物化探数据和航磁数据想快速圈出一批找矿靶区放在十年前标准流程是证据权法或多因子叠加构造权重全靠经验现在这套流程的底层逻辑变了——把矿点当标签、把多源地质特征当特征交给机器学习模型去拟合。这份机器学习成矿预测源码包就是把整条流水线串起来的工程代码从多源数据栅格化、正负样本构建、随机森林与 XGBoost 训练、AUC 评估到靶区概率栅格输出全部有可跑的 Python 代码。适合地质工程师、GIS 开发以及想把手头化探物探数据快速变成靶区图的从业者拿到手按顺序执行就行不用自己从零造轮子。2. 数据与样本构建从地质图到特征矩阵的三步2.1 多源数据栅格化分辨率与坐标系统一成矿预测的第一步永远不是建模而是把地质图、化探、物探、遥感这些格式各异的数据压到同一张网格上。常见做法是统一到一个坐标系比如 CGCS2000 或 WGS84 UTM 投影再统一分辨率——我一般用 50m 或 100m 像元矿田尺度用 100m 足够矿区精细评价才加密到 50m。地层、岩体、断裂这些矢量数据要转成栅格特征地层按年代编码成整数栅格断裂和岩体先栅格化成 0/1 二值图再做欧氏距离变换得到到断裂距离这类连续特征。化探数据因为是不规则采样点要先插值成规则网格常用的有反距离加权IDW或普通克里金。import rasterio import numpy as np from scipy.ndimage import distance_transform_edt def build_feature_stack(base_path): # 读入已经配准并裁剪到同一范围的化探栅格 with rasterio.open(f{base_path}/geochem_Cu.tif) as src: cu src.read(1) profile src.profile # 统一用这个profile写输出栅格 # 断裂距离场先读二值栅格值为1的像元表示有断裂穿过 with rasterio.open(f{base_path}/fault_binary.tif) as src: fault_bin (src.read(1) 0).astype(np.uint8) # 欧氏距离变换得到每个像元到最近断裂的像素距离 fault_dist distance_transform_edt(1 - fault_bin) * profile[transform].a # 把所有特征沿通道方向堆叠形状为 (通道数, 高, 宽) stack np.stack([cu, fault_dist], axis0) return stack, profile, fault_distdistance_transform_edt计算的是像素距离乘以transform.a是因为 transform 矩阵里记录了每个像元的实际地面尺寸比如 100m这一步把像素距离换算成米。堆叠后的stack就是后面训练用的特征立方体profile保存了坐标系和变换参数最后把概率预测结果写回栅格时直接复用。2.2 正负样本采集随机采样与缓冲区采样差异样本构建是成矿预测里最容易被低估的一步。正样本是已知矿点直接取矿点所在像元的特征值即可负样本则不能随便撒——如果你在距离矿点几十米的地方取负样本模型学到的是距离矿近这个假象而不是真正的地质控矿规律。我一般会用缓冲区剔除以已知矿点为中心做 500m 缓冲区缓冲区外随机采样作为负样本。正负样本比例建议控制在 1:3 到 1:5太高的负样本比例会把模型推成全部预测为负的保守形态。import geopandas as gpd import numpy as np from shapely.geometry import Point def generate_negative_samples(ore_points, study_area, exclude_dist500, n_samples3000): # 矿点周围500m范围内不采负样本 buffer gpd.GeoDataFrame(geometryore_points.geometry.buffer(exclude_dist), crsore_points.crs) merged_buffer buffer.geometry.union_all() x_min, y_min, x_max, y_max study_area.total_bounds candidates [] while len(candidates) n_samples: xs np.random.uniform(x_min, x_max, 10000) ys np.random.uniform(y_min, y_max, 10000) pts [Point(x, y) for x, y in zip(xs, ys)] # 剔除落在研究区边界外的点和缓冲区内的点 for p in pts: if study_area.contains(p) and not merged_buffer.contains(p): candidates.append(p) if len(candidates) n_samples: break return gpd.GeoDataFrame(geometrycandidates, crsstudy_area.crs)exclude_dist是最关键的参数它控制负样本与矿点的最小距离。100m 网格下我建议至少 500m也就是 5 个像元如果你研究的矿体本身很大比如斑岩型铜矿缓冲区要放大到 1–2km否则负样本里混入矿化蚀变带的像元模型会显著低估靶区面积。union_all把矿点的缓冲区合并成一个整体避免逐个圆判断造成采样点落在两个缓冲区的缝隙里。2.3 样本特征对齐与归一化样本采集完成之后要把每个样本点的多源特征对齐。规则是每个样本点正样本是矿点负样本是随机点通过其坐标反查特征立方体中对应位置的像元值抽取出一行特征向量。这个过程我列一个清单用rasterio.index(x, y)将地理坐标转为行列号注意坐标系必须一致特征包括化探元素含量、断裂距离、岩性编码、航磁异常值、遥感蚀变指数岩性这类类别特征用 one-hot 编码不要让模型把地层 3和地层 7当成有序关系数值特征全部过StandardScaler很多树模型虽然不要求归一化但在特征量纲差异极大时归一化能加快调参时的收敛速度from sklearn.preprocessing import StandardScaler # samples 为DataFrame包含正负样本点的特征列 feature_cols [Cu, Pb, fault_dist, litho_code, mag_anomaly] X samples[feature_cols].values y samples[label].values # 归一化后再训练注意scaler只fit训练集避免信息泄漏 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_val_scaled scaler.transform(X_val)这里的核心坑是scaler一定不能在全量数据上 fit。正确做法是先切分样本按空间切分下一章细说再在训练集上fit_transform验证集和测试集只做transform。否则验证集的信息在训练前就已经参与了均值和方差的估计评估结果会偏乐观。特征列的选择上我建议先全量跑一版随机森林看特征重要性再剔除重要性极低且与其他特征高相关的列比如 Au 与 As 常高度相关保留一个即可。3. 模型训练与调参随机森林与 XGBoost 怎么选3.1 基线模型逻辑回归帮你看数据下限不要一上来就上树模型集成。先用逻辑回归跑一遍它会强迫你审视特征本身是否具备线性可分性。逻辑回归在成矿预测里往往 AUC 在 0.75 左右这个数值作为后续模型的对照基线非常有用——如果随机森林只比逻辑回归高 0.05那说明数据里新增的非线性信息有限问题可能出在特征工程而不是模型复杂度上。from sklearn.linear_model import LogisticRegression from sklearn.model_selection import cross_val_score # C1.0 为正则化强度小数据量上默认即可 lr LogisticRegression(C1.0, max_iter1000, random_state42) scores cross_val_score(lr, X_train_scaled, y_train, cv5, scoringroc_auc) print(fLR CV-AUC: {scores.mean():.3f} ± {scores.std():.3f})cv5是五折交叉验证每个折的 AUC 差异如果超过 0.1说明样本空间分布很不均匀后续要用分组交叉验证。逻辑回归的系数也可以输出正系数最大的特征就是与成矿正相关的因子负系数最大的则相反这一步能帮你发现解释性强的控矿因素比如到断裂距离系数为负表示离断裂越近成矿概率越高这在地质上完全说得通。3.2 随机森林为什么它是成矿预测的默认选择随机森林在这个场景下的地位几乎不可替代。它对特征量纲不敏感能处理特征交互天然支持并行最重要的是能输出特征重要性地质学家要解释哪个因素在控矿时这是刚需。我见过的成矿预测项目里随机森林的 AUC 通常比逻辑回归高 0.1 左右达到 0.85–0.9再往上去就依赖更细致的特征工程。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import GroupKFold # 按空间块划分样本同一矿化带的样本必须分到同一折 block_id samples[block_id].values # 每个空间区块一个编号 gkf GroupKFold(n_splits5) rf RandomForestClassifier( n_estimators300, max_depth8, min_samples_leaf10, max_featuressqrt, n_jobs-1, random_state42 ) for train_idx, val_idx in gkf.split(X_train_scaled, y_train, groupsblock_id): rf.fit(X_train_scaled[train_idx], y_train[train_idx]) print(val AUC:, roc_auc_score(y_train[val_idx], rf.predict_proba(X_train_scaled[val_idx])[:, 1]))GroupKFold是整个训练环节最重要的细节。普通KFold是随机划分但如果同一个矿化带延伸几公里矿点落在训练集和验证集各一半模型相当于开卷考试。GroupKFold要求你预先给每个样本分一个空间区块号按区块划分确保验证集的样本与训练集在空间上完全不相邻。max_depth8和min_samples_leaf10是两个最有效的正则化参数树太深对空间噪声拟合过强建议max_depth从 8 起步调到 12 如果验证集 AUC 没有明显收益就守住 8。3.3 XGBoost 调参学习率、树深与早停XGBoost 在成矿预测里的表现通常略优于随机森林但差距不是玄学主要来自它对特征交互的建模更精细。代价是需要多设几个正则化参数。我常用的配置是learning_rate0.03max_depth5subsample0.8colsample_bytree0.8配合early_stopping_rounds50。这个组合在多个项目里都能稳定收敛。import xgboost as xgb xgb_model xgb.XGBClassifier( n_estimators1000, learning_rate0.03, # 学习率越小越稳但需要的树更多 max_depth5, # 树深超过6容易过拟合空间噪声 subsample0.8, # 每棵树随机用80%样本 colsample_bytree0.8, # 每棵树随机抽80%特征 reg_lambda1.0, # L2正则控制权重 eval_metricauc, early_stopping_rounds50, random_state42 ) xgb_model.fit( X_train_scaled, y_train, eval_set[(X_val_scaled, y_val)], verboseFalse # 关闭训练日志只看最终结果 )early_stopping_rounds50的意思是连续 50 轮验证集 AUC 没有提升就停止训练返回最优模型。这里必须搭配一个稳定的验证集——如果你用随机划分早停选出的最优轮次可能恰好是过拟合的某一个中间态搭配上一节的GroupKFold划分才能保证早停选出的模型具备真实泛化能力。我一般会跑两遍第一遍用默认参数跑全量看 AUC 和特征重要性第二遍把learning_rate从 0.03 微调到 0.01同时把n_estimators上限提到 1500通常 AUC 能再涨 0.02–0.05。4. 评估与靶区输出AUC 之外还要看空间表现4.1 混淆矩阵和命中率的正确打开方式成矿预测的评估指标与一般分类问题有一个关键差异你真正关心的不是把所有负样本全找对而是预测的靶区里能不能尽量多地把已知矿点圈住。术语叫命中率或召回率。一个 AUC 0.9 的模型如果输出的高概率区只覆盖 30% 的已知矿点对找矿没有实际价值。因此我评估时固定看三件事AUC 是否大于 0.8、在 50% 面积约束下能命中多少已知矿点、高概率区在空间上是否聚集在断裂与岩体接触带附近。from sklearn.metrics import roc_auc_score, classification_report y_prob rf.predict_proba(X_val_scaled)[:, 1] auc roc_auc_score(y_val, y_prob) # 假设靶区面积占总面积的30%取概率最高的30%像元作为靶区 area_ratio 0.3 threshold np.percentile(y_prob, (1 - area_ratio) * 100) hits (y_val 1) (y_prob threshold) hit_rate hits.sum() / (y_val 1).sum() print(fAUC{auc:.3f}, 30%面积命中率{hit_rate:.2%})np.percentile(y_prob, 70)的含义是概率最高的 30% 像元进入靶区然后检查这些栅格里包含了多少已知矿点。这个指标比 AUC 更贴近地质决策——你圈出的靶区面积固定能捕到多少矿才是核心诉求。我见过 AUC 0.92 但命中率只有 50% 的模型原因就是高概率区集中在某一岩体边缘而忽略了另一片热液活动的区域这种模型在业务上不如 AUC 0.85 但命中率 75% 的模型。4.2 概率阈值选择别用 0.5用召回率反推默认的 0.5 分类阈值在成矿预测里几乎永远是错的。因为正负样本比例是人工设定的并非自然分布0.5 对应的地质含义是成矿概率 50%——这在找矿场景里过于苛刻。常见做法是把阈值与面积预算绑定你打算投入钻探验证的面积占总面积的 10%就让阈值取概率分布的第 90 百分位如果想保守一点按 85% 命中率反推阈值。# 按目标命中率反推阈值 target_recall 0.85 idx_sorted np.argsort(y_prob)[::-1] cum_ore_hits np.cumsum(y_val[idx_sorted]) total_ore y_val.sum() cut np.searchsorted(cum_ore_hits, target_recall * total_ore) threshold_recall y_prob[idx_sorted[min(cut, len(y_val) - 1)]] print(f85%命中率对应阈值{threshold_recall:.4f})np.argsort(y_prob)[::-1]把概率从高到低排序cum_ore_hits是累加命中的矿点数。searchsorted找到累计命中数达到总矿点数 85% 的位置该位置对应概率值就是你要的阈值。这个方法保证按这个阈值圈靶区理论上能覆盖 85% 的已知矿点同时把无效面积控制在最小。4.3 概率栅格输出与在 GIS 里叠加验证模型训练完后把整幅研究区内每个像元的特征矩阵都推算一遍得到一副连续的成矿概率栅格这才是交付给地质人员使用的最终成果。输出时直接把概率数组写入 GeoTIFF注意复用之前特征栅格的profile保证坐标系和分辨率完全一致。import rasterio import numpy as np # all_features 为全图特征立方体形状 (特征数, 高, 宽) h, w all_features.shape[1], all_features.shape[2] X_all all_features.reshape(all_features.shape[0], -1).T # (像元数, 特征数) X_all_scaled scaler.transform(X_all) prob_flat rf.predict_proba(X_all_scaled)[:, 1] prob_grid prob_flat.reshape(h, w).astype(float32) # 裁剪掉屏蔽区域的概率值 prob_grid[nodata_mask] np.nan with rasterio.open(output/prob_rf.tif, w, **profile) as dst: dst.write(prob_grid, 1)这里最容易翻车的动作是X_all的特征排列顺序必须与训练时的X_train列顺序完全一致少一列或列顺序不同会导致概率值彻底错乱。我自己的习惯是训练前把特征列名列表存成一份feature_cols.json预测时读取该文件名列表再构造X_all从根上杜绝顺序不一致的问题。输出后的 GeoTIFF 直接拖进 QGIS与地质图叠加看高概率区是否落在已知成矿构造带上这一步是模型与经验判断的节点对齐不可跳过。5. 避坑排查成矿预测里五个翻车场景5.1 现象AUC 高达 0.95但靶区分布完全不合理原因空间自相关导致数据泄漏。同一矿化带内的像元特征高度相似随机划分训练集和测试集时测试样本相当于训练样本的近亲模型靠记忆邻近像元而非识别地质规律得分。解决改用GroupKFold或者缓冲区切分。将研究区划分为 2km×2km 的网格块每个块内样本作为整体放入同一折保证同一区块不会同时出现在训练集和验证集。切分后 AUC 通常会掉 0.05–0.1这是真实的泛化水平不要恐慌。5.2 现象训练集 AUC 高验证集 AUC 低原因树模型过拟合空间细节。成矿数据本身噪声大化探元素含量分布受采样密度影响树深度过深时会把局部异常当规律。解决优先调max_depth和min_samples_leaf。把max_depth从默认 10 以上收紧到 6–8min_samples_leaf提到 10–20随机森林和 XGBoost 的性能下降微乎其微泛化稳定性明显提升。如果做 XGBoost 记得用早停不要靠固定n_estimators硬跑。5.3 现象正负样本比例 1:1 时 AUC 好看样本比例一变结果剧烈波动原因负样本采样方式太随意。每次随机采的负样本点特征分布不同导致模型评估结果不稳定。解决负样本采完存成固定文件每次试验用同一批样本。更稳的做法是采用多次随机负采样的平均结果跑 5 次每次重新采负样本报告 AUC 的均值与标准差标准差超过 0.03 说明负样本空间分布影响过大需要增大负样本数量或调整缓冲区。5.4 现象特征排序与地质认识矛盾比如距断裂距离重要性极低原因断裂距离场可能没有正确生成。常见问题包括断裂线 shp 的坐标系与研究区栅格不一致、距离变换前没有将二值栅格翻转导致距离方向颠倒、缓冲区内负采样把断裂附近样本全部剔除从而掩盖了该特征的作用。解决逐一排查。坐标系不一致可以在rasterio.open读取后检查crs与profile[crs]是否一致。距离方向问题用distance_transform_edt(1 - fault_bin)确认二值图中断裂处为 0、背景为 1。缓冲区的负样本剔除范围不能覆盖整条断裂适当缩小exclude_dist或在断裂附近保留部分负样本点。5.5 现象输出概率栅格有带状或放射状条纹原因插值参数不当。化探数据做 IDW 插值时幂指数过大局部高值点向外扩散衰减太剧烈产生以采样点为圆心的放射状纹理克吕金变差函数拟合不当也会产生类似效果。解决IDW 的幂指数默认 2 建议降到 1.5 试跑对比条纹是否减弱更推荐使用普通克里金并调整块金值块金值太小会导致插值面过度贴近训练点。插值完成后做一遍低通滤波或中值滤波比如 3×3 窗口能有效压掉网格化噪声代价是靶区边界略变平滑。6. 往新研究区迁移改参数而不是改代码拿到这份源码包之后换到另一个研究区时你不需要动核心算法把配置文件改清楚即可。我建议先把以下参数填好研究区边界 shp、目标坐标系、栅格分辨率、特征列表、已知矿点 shp、负样本采样距离与比例。下面这个配置是我惯用的结构# config/study_area.yaml study_area: shp_path: D:/geodata/new_area_boundary.shp target_crs: EPSG:32650 grid_resolution: 100 features: geochem: [Cu, Pb, Zn, Au, As] geophysics: [mag, gravity] structure: [fault_dist, granite_dist] litho: [litho_code] samples: positive_shp: ore_points.shp negative_ratio: 3 exclude_buffer_m: 500 block_size_m: 2000换研究区时我会强制自己走一遍流程先加载配置画一张研究区范围图确认边界和坐标系再输出特征栅格叠加上已知矿点目视检查是否存在矿点落在栅格空值区的对齐问题接着跑一个快速随机森林打印特征重要性并和地质背景比对一遍比如花岗岩距离特征优先于地层编码通常说明该区域岩浆热液作用主导成矿逻辑上成立。从换区到拿到第一版靶区图我控制的正常周期是一天。实际工程里最容易拖慢我的不是模型精度而是数据坐标系的来回折腾。我现在每接手一个新研究区第一件事就是把所有矢量数据统一to_crs(target_crs)所有栅格用rasterio重采样到同一分辨率并裁剪到同一范围这一步做完后面的一切都顺畅。曾经有一次跳过统一坐标系直接训练AUC 虽然不低但靶区图与研究区边界错位了几公里整个上午都在排查特征对齐问题。从那以后我配置里首先写死target_crs并强制每个数据源都做坐标系检查不再信任任何默认投影。希望这份源码和这套流程能帮你在自己的数据上少走同样的弯路。本文还有配套的精品资源点击获取