机器学习复现造山型金矿黄铁矿微量元素分析：从数据预处理到SHAP解释

简介：面向地质学与数据科学交叉领域研究者的一份复现论文资源，聚焦造山型金矿床中黄铁矿微量元素变化规律，可辅助理解金矿化阶段判别与成矿温度预测。文档基于Python完整演示数据清洗与预处理、KNN插补和中心对数比转换、PCA与PLS-DA降维判别、随机森林分类与回归建模，并通过网格搜索优化模型参数；同时配有各步骤代码注释、结果解读和多种可视化图表，能帮助读者掌握将大数据和机器学习引入矿床学研究的实际分析流程，并根据自己的数据集灵活调整参数。压缩包仅含1个docx文档，约20KB，篇幅精炼但覆盖数据预处理、统计分析、模型评估与解释全链路。目前已有52人学习，适合矿物学、地球化学和机器学习方向科研工作者参考借鉴。

1. 造山型金矿床的黄铁矿微量元素：为什么大数据分析与机器学习值得复现这篇论文

手上有几千个黄铁矿测点的 LA-ICP-MS 数据，想厘清哪个期次的黄铁矿真正载金、哪些微量元素能作为找矿标志——这是造山型金矿床研究中很常见的处境。传统做法是画 Au-As 二元图解、做相关性矩阵，一旦数据量超过两三千个测点，纸面图件就开始失效：不同期次的黄铁矿在元素组成上高度重叠，肉眼和二维散点都分不开。这篇论文的价值，是把大数据分析与机器学习算法搬进矿床地球化学，让模型去抓多元素联动的非线性指纹。复现它的意义不只是跑通代码，而是把一套可迁移的分析流程拿到自己的数据集上，量化判断哪一期黄铁矿的 Au、As、Te 组合真正指示富矿体。适合的人群是矿床学研究生、地质大数据工程师和有 LA-ICP-MS 数据但不会建模的从业者。

2. 复现前的数据准备：把 LA-ICP-MS 测点整理成机器学习能吃的表格

2.1 原始数据长什么样：一行一个测点，一列一种元素

造山型金矿中的黄铁矿通常被划分为三期：沉积变质期 Py1、主成矿期 Py2、晚期热液期 Py3。每期黄铁矿的微量元素组成有差异，但这种差异不是单一元素能刻画的——Py2 的 As、Au、Te 偏高，Py3 的 Sb、Pb 偏高，Py1 的 Co、Ni 偏高且 Au 普遍低于检出限。论文里的大数据分析起点，就是把这些测点数据组织成结构化表格。

实际从 ICP-MS 仪器导出、或者从论文附件下载的原始表，长这样：

常见的坑是：每列元素单位不一致。Au 用 ppb，其他主微量元素用 ppm，到来的数值差几百到几千倍。正规做法是进模型之前先做量纲统一，要么全部对数变换、再做标准化，要么把单位换算成同一量纲后在行内做归一化。论文里的处理路径更倾向"对数变换 + 标准化"，因为黄铁矿微量元素跨数量级分布，直接用原始数值会让高含量元素支配模型权重。

另一个需要盯紧的字段是"期次"列。这是你的监督标签，也是后续分类模型的 y。它来自岩相学和矿相学观察——颗粒形貌、背散射图像灰度、与脉体的穿插关系，由地质人员在点测时人工判定。论文复现时，如果原始数据文件里没有这一列，就得根据样品的产状信息手动补标，这一步的工作量通常是数据预处理里最大的。

2.2 构建训练集：样品分组、标签编码与数据划分的先后顺序

拿到原始表后，第一件事不是标准化，而是先检查样本量分布。Py2 往往占一半以上，Py3 可能只有不到一百个测点——这种不平衡是造山型金矿数据的常态。处理不平衡的最佳顺序是：先分组、再划分训练集和测试集、最后做标准化。顺序错了，后面模型评估全是虚高。

import numpy as np import pandas as pd from sklearn.model_selection import GroupKFold from sklearn.preprocessing import StandardScaler df = pd.read_csv("pyrite_trace_element.csv") print(df["phase"].value_counts()) # 只保留元素列作为特征 element_cols = ["Au_ppb", "As_ppm", "Sb_ppm", "Te_ppm", "Bi_ppm", "Ag_ppm", "Cu_ppm", "Pb_ppm", "Zn_ppm", "Co_ppm", "Ni_ppm", "Se_ppm"] X = df[element_cols].copy() y = df["phase"].astype(str).copy() groups = df["sample_id"].copy() # 样品ID # 对数变换：压低数量级差异，把乘法关系变成加法关系 for col in X.columns: X[col] = np.log1p(X[col].clip(lower=0))

这段代码做了两件事：把所有元素做 log1p 变换，以及保留样品 ID 作为分组依据。先说 log1p 的必要性——黄铁矿中 Au 可以从小于 1 ppb 到几千 ppb，跨三个数量级；如果不做对数变换，树模型虽然对单调变换不敏感，但后续要做的标准化、聚类、SHAP 依赖图都会倾向于高值点。clip 到 0 是为了把负检出值提前兜住，这一步在第 5 章的避坑清单里会详细展开。

分组这个变量现在不用，但它决定了交叉验证的切分方式。同一个黄铁矿颗粒上打的多个激光剥蚀坑，微量元素高度自相关；如果按测点纯随机切分，训练集里出现颗粒 A 的 3 个点、测试集里出现颗粒 A 的另外 2 个点，模型就等于带着答案考试。GroupKFold 按样品 ID 切分，能保证同一个样品的点不会同时出现在训练集和测试集里，这是复现此论文分类结果的关键细节。

2.3 检出限、缺失值与异常值的处理：低于检出限的数值是灾难源

LA-ICP-MS 数据里，低于检出限的测点很常见。很多论文的数据表里，未检出点被写成"0"或空值。如果直接丢给机器学习模型，会造成两个问题：一是 0 值在对数变换后变成极大的负数，被模型当作一类特殊信号；二是一旦测试集里未检出比例和训练集不同，模型立即失效。

常见的处理顺序是这样的：先把 < DL 的数值替换为 DL/2（检出限的一半），然后加一列标志变量区分"实测值"和"未检出值"，再做对数变换。比如 Au 的检出限是 0.5 ppb，那么所有低于 0.5 的测量值填 0.25，同时新建一列 Au_belowDL，值为 1 表示未检出，0 表示实测。这样模型既能看到数值大小，也能区分数据质量。

dl_dict = {"Au_ppb": 0.5, "As_ppm": 2, "Sb_ppm": 1, "Te_ppm": 0.5, "Bi_ppm": 0.2, "Ag_ppm": 0.3, "Cu_ppm": 5, "Pb_ppm": 2, "Zn_ppm": 5, "Co_ppm": 2, "Ni_ppm": 5, "Se_ppm": 1} for col in element_cols: dl = dl_dict[col] below_col = f"{col}_belowDL" df[below_col] = (df[col] < dl).astype(int) df[col] = np.where(df[col] < dl, dl / 2.0, df[col])

替换为 DL/2 而不是 0，是为了让未检出值保留在含量区间的下端，而不是凭空造出一个人为的零值。加标志列是论文中常被忽略但很实用的一步——很多时候"这个元素到底有没有被检出"本身就是重要的地质信息，Py1 黄铁矿的 Au 和 Te 普遍未检出，这个模式对分类器来说和低含量同样有区分度。

做完这一步，再把 data 和 belowDL 标志列合并，加上 phase、sample_id，才是一个能进模型的完整数据集。标准化放在划分训练集之后做，第 5 章会讲为什么。

3. 用随机森林分类器区分黄铁矿期次：特征重要性比散点图诚实

3.1 为什么选随机森林而不是深度网络：地质数据不需要花哨模型

数据集经过预处理后，特征维度大概在 12 个元素列加上 12 个检出标志列，样本量是两三千个测点。这种规模下，深度神经网络没有优势；样本少、特征维度中等、噪声大，正是随机森林和梯度提升树的舒适区。论文选择随机森林而不是神经网络的理由很实在：一是特征重要性可以直接输出，对应到具体元素做地质解释；二是对类别不平衡和异常值有天然韧性；三是训练时间短，不需要调学习率、网络宽度这些玄学参数。

梯度提升树在这个任务里也很常用，但 LightGBM 对噪声敏感，在 LA-ICP-MS 这种带强噪声的数据上容易过拟合。随机森林的每棵树独立采样、结果取平均，对单个异常测点的抗性更好。如果数据量超过一万个测点，再考虑换成 LightGBM 调参，现在这个规模用随机森林更稳。

3.2 用 GroupKFold 做交叉验证：核心代码与参数说明

先划分折再标准化是硬性要求。组内划分保证同一个黄铁矿颗粒的测点不会跨训练和测试集；划分完成后，StandardScaler 只在训练折上 fit，再转换测试折。

from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import classification_report, confusion_matrix from sklearn.pipeline import Pipeline from scipy.stats import randint from sklearn.model_selection import RandomizedSearchCV feature_cols = element_cols + [f"{c}_belowDL" for c in element_cols] X = df[feature_cols].copy() outer_cv = GroupKFold(n_splits=5) train_idx, test_idx = next(iter(outer_cv.split(X, y, groups))) X_train, X_test = X.iloc[train_idx], X.iloc[test_idx] y_train, y_test = y.iloc[train_idx], y.iloc[test_idx] scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test) clf = RandomForestClassifier( n_estimators=500, min_samples_leaf=3, max_depth=None, class_weight="balanced", random_state=42, n_jobs=-1 ) clf.fit(X_train_scaled, y_train) y_pred = clf.predict(X_test_scaled) print(classification_report(y_test, y_pred, digits=3))

n_estimators=500 是收敛的稳妥值，树太少模型抖动大，树太多纯耗时。min_samples_leaf=3 是这里最关键的参数——叶节点至少 3 个样本，能有效抑制模型去死记单个异常测点的含量组合。class_weight="balanced" 自动给数量少的 Py3 类更大的惩罚权重，解决数据不平衡。max_depth 不限制，靠 min_samples_leaf 控制复杂度即可。

所有特征在进入树前都做了标准化。严格说，树模型不要求特征标准化，但下面要对比特征重要性和做聚合分析，标准化后各特征的尺度一致，输出的重要性分数更适合跨元素比较。特征全谱包括 12 个元素列和 12 个 belowDL 标志列，一共 24 维。

3.3 读特征重要性：哪些元素在真正区分期次

训练完成后立刻输出特征重要性排序：

import matplotlib.pyplot as plt importance = pd.DataFrame({ "feature": feature_cols, "importance": clf.feature_importances_ }).sort_values("importance", ascending=False) print(importance.head(12))

一个典型的复现结果里，As_ppm、Au_ppb、Te_ppm、Sb_ppm 会排在前四。这和造山型金矿的地球化学认知一致：主成矿期 Py2 富集 As 和 Au，伴随 Te、Bi、Ag；Py3 晚期则 Sb、Pb 抬升。真正的价值在于，特征重要性会告诉你 Py1 和 Py2 的区分主力是哪几个元素——如果 Au 的重要性极低，说明这个矿区的黄铁矿载金性和期次关系不紧密，需要换找矿指标。

看重要性图的时候要留个心眼：随机森林的 feature_importances_ 偏向高基数特征，对相关性强的元素组合会分散权重。比如 Au 和 As 高度相关时，重要性会在这两个元素之间分摊，单看排名可能低估 Au 的贡献。所以第 4 章要用 SHAP 验证。

4. 用回归模型量化 Au 含量：从微量元素到品位约束

4.1 目标变量怎么定：直接预测 Au 还是预测 Au/As 比值

分类任务回答了"这粒黄铁矿属于哪一期"，回归任务回答的是"这粒黄铁矿到底带了多少 Au"。论文里的机器学习约束，有一部分就是把 Au 作为目标变量，用其他元素做特征做回归预测——本质上是在模拟黄铁矿的载金能力，为圈定富矿段提供指标。

这里有个关键选择：直接预测 Au_ppb 还是预测 log(Au)？Au 含量跨 0.1 到 5000 ppb，直接回归会让模型拼命拟合高值点，低值点几乎被忽略。常见做法是对目标变量做 log1p 变换，让误差在数量级意义上均匀分布。

另一个更稳健的做法是预测 Au/As 比值。造山型金矿中 As 与 Au 常呈强正相关，Au/As 比值可以消除"热液强度"这一共性变量的影响，突出不同成矿阶段的元素分异。但这个选择要看你复现的论文具体做了什么约束——如果原论文聚焦品位估算，用 log(Au)；如果聚焦元素共生分异，用 Au/As 更适合。

4.2 LightGBM 回归与 SHAP 拆解：代码与参数落地

分类用随机森林，回归可以尝试 LightGBM 加 SHAP。LightGBM 对 log1p 变换后的连续目标拟合能力强，且 SHAP 能解释每一个测点预测值里各个元素的贡献方向。先把数据切成训练测试，再做训练：

import lightgbm as lgb import shap from sklearn.metrics import mean_absolute_error, r2_score target_col = "Au_ppb" feature_cols_r = element_cols + [f"{c}_belowDL" for c in element_cols] # 注意：特征里不要包含 Au 自身，但 As、Te 都要保留 X_r = df[feature_cols_r].copy() y_r = np.log1p(df[target_col].clip(lower=0)) Xr_train, Xr_test, yr_train, yr_test = train_test_split( X_r, y_r, test_size=0.25, random_state=42, stratify=df["phase"] ) model = lgb.LGBMRegressor( n_estimators=800, learning_rate=0.05, num_leaves=31, min_child_samples=10, subsample=0.8, colsample_bytree=0.8, random_state=42 ) model.fit(Xr_train, yr_train, eval_set=[(Xr_test, yr_test)], callbacks=[lgb.early_stopping(50)]) pred_log = model.predict(Xr_test) print("R2:", r2_score(yr_test, pred_log)) print("MAE(log):", mean_absolute_error(yr_test, pred_log)) # 反变换回原始 Au 单位，便于和地质解释对应 pred_au = np.expm1(pred_log) true_au = np.expm1(yr_test)

n_estimators=800 配合 early_stopping，模型会在验证集 50 轮不提升时自动截断，防止过拟合。learning_rate=0.05 偏低，配合大树数量是 LightGBM 的常见组合。subsample=0.8 和 colsample_bytree=0.8 做行和列采样，增加模型鲁棒性——LA-ICP-MS 数据里有设备漂移造成的系统偏差，行列采样能削弱这种干扰。

SHAP 解释放在模型训练后：

explainer = shap.TreeExplainer(model) shap_values = explainer.shap_values(Xr_test) shap.summary_plot(shap_values, Xr_test, max_display=15)

SHAP 散点图里每行是一个特征，横轴是 SHAP 值——正值表示该特征把 Au 预测值推高。造山型金矿数据上最常见的图形是：As ppm 和 Te ppm 的 SHAP 值随含量增大而显著右移，说明这两个元素对 Au 的贡献是强正向的；而 Co、Ni 的 SHAP 点集中分布在 0 附近，说明它们对 Au 预测几乎没有信息量。这比随机森林的特征重要性多了一个方向维度，能直接回答"As 高一定对应 Au 高吗"——不一定，要看 As 是在哪个期次、和什么元素组合在一起的。

4.3 回归模型的边界：预测不是外推，别拿模型当勘探承诺

模型 R² 在测试集上达到 0.85 甚至 0.9，容易让人激动，但要清楚边界：训练数据来自同一矿区的有限样品，模型学到的是这个矿区黄铁矿的元素共生规律。换一个矿区，热液流体性质变了，As-Au 相关性可能完全失效。参考论文的做法，应该报告模型在矿区内的适用区间——比如 Au 预测值在 10 到 2000 ppb 范围内可靠，低于 1 ppb 时模型倾向于高估。

还有一层地质约束：回归模型预测的是黄铁矿颗粒本身的 Au 含量，不是矿石品位。矿石品位还取决于黄铁矿的丰度、粒度、赋存状态。复现论文时，你要把模型的输出定义为"单颗粒载金能力指标"，它可以作为找矿评价的辅助因子，但要落到资源量估算上，必须乘以矿相学统计的黄铁矿体积分数。这是模型边界，不是模型缺陷。

5. 复现论文的避坑清单：数据、标签、模型的三类翻车现场

5.1 低于检出限的值被填成 0，模型学到的是"缺失模式"而不是丰度梯度

现象：模型在分类任务上验证精度极高，达到 0.97，但画出的特征重要性排序里，Au 排到了倒数几位，而 belowDL 标志列排名第一。整个结果和地质认知完全冲突。

原因：未检出值填 0 后，log1p(0)=0，与真实低含量值形成两簇分布——一簇是已检出的低值，另一簇是未检出填零值。树模型抓到了这个人为的断崖，"是否存在实测值"成了最省力的分类特征，真实的含量梯度反而被忽略。

解决：按第 2 节的处理流程，未检出值统一替换为 DL/2，同时保留标志列。填 DL/2 之后，未检出值落入含量分布的低尾，模型被迫去学习真实含量模式。如果跑完特征重要性发现 belowDL 列仍然碾压元素列，说明检出率过低，这时应该考虑删掉该元素列，或者把任务从定量分析降级为检出/未检出二分类。

5.2 标签泄漏：把人工编录信息混进特征，特征重要性全乱

现象：模型精度高得离谱，但特征重要性榜首出现了一个不是微量元素的字段——比如"颗粒形态""反射率等级"这类编录信息。

原因：原始数据表里除了 LA-ICP-MS 元素之外，往往还有岩相学描述列。颗粒形态、颜色、蚀变强度这些东西本身就是地质人员判定期次的依据，等于把自变量从元素偷偷换成了接近标签的变量。模型学到的是"长得像 Py2 的就是 Py2"，而不是微量元素指纹。

解决：特征矩阵里只保留仪器输出的元素含量列和对应检出标志列，所有人工编录字段一律不放。检查方式很简单：把特征重要性输出里出现任何非元素字段的情况视为数据错误，立即重做预处理。这个坑在复现早期最容易被忽略，因为多余列通常排在数据表末尾。

5.3 同一个黄铁矿颗粒的测点被分成训练和测试两份，验证分数虚高

现象：用普通 KFold 交叉验证时 F1 分数 0.93，换成 GroupKFold 后掉到 0.81，差距明显。很多人此时怀疑分组出了问题，其实分组才是对的。

原因：同一个颗粒上的剥蚀点，微量元素空间自相关极强——同一颗粒打 5 个点，Au 含量相差不到 20%。纯随机划分时，同一个颗粒的一部分点进了训练集，另一部分进了测试集，模型等于考试时看到了同题答案。KFold 分数虚高的幅度，取决于每个颗粒的平均测点数。

解决：划分交叉验证时用 GroupKFold，groups 参数传样品 ID 或颗粒 ID。更严格的做法是按钻孔或按勘探线分组——如果目的是把模型推广到未知钻孔，就要以钻孔为组，否则训练集和测试集共享了钻孔级别的系统性偏差。选哪个粒度取决于你的业务目标，地层对比用样品分组，勘探预测用钻孔分组。

5.4 类别不平衡被忽略，模型把所有测点都判成 Py2

现象：Py1、Py2、Py3 三类样本比例大约是 3:5:2，随机森林预测结果里 Py2 的召回率 0.95，Py1 和 Py3 召回率不到 0.3，模型输出几乎全是 Py2。

原因：没有处理不平衡时，树分裂点的信息增益被多数类主导，少数类即使能被识别也被整体淹没。这在造山型金矿数据里几乎必然发生，Py3 通常只在晚期构造裂隙里发育，测点天然少。

解决：随机森林里加 class_weight="balanced"，等价于给少数类样本乘以一个较大的权重；如果比例悬殊到 10:1，先用 SMOTE 做过采样再训练。但要注意，SMOTE 合成样本是插值出来的，在微量元素高度非线性的数据上可能造出不符合物理化学规律的组合。论文复现时优先用 class_weight 和 LightGBM 的 is_unbalance=True，少用 SMOTE。

5.5 标准化顺序出错：用全数据拟合 scaler 造成信息泄漏

现象：模型评估指标很漂亮，但换到新矿区的数据上性能崩塌，R² 下降一半。

原因：预处理时先在全数据集上调用 scaler.fit_transform，然后才切分训练测试集。均值和方差从测试集里"偷看"了信息，导致测试集的标准化结果不纯粹。这在早期复现中非常常见，因为很多教程都是先清洗后切分，但带监督的建模流程里"先切分再预处理"才是更稳妥的顺序。

解决：把标准化写进 Pipeline，或者严格执行 fit_transform 只在训练集，transform 只在测试集。值得多说一句：如果是无监督的聚类和可视化，全数据集标准化没问题；只要涉及模型评估，就一定要先分组切分、再 fit 标准化器。

6. 把模型结果落回地质解释：SHAP 依赖图与元素共生组合

模型跑通只是复现的中点，终点是拿结果推进地质认识。有两个技巧我习惯在收尾时用：一是画 SHAP 依赖图检查成矿元素对 Au 的协同作用；二是做特征聚类，把微量元素划分成几个共生组合，锁定找矿指示标志。

shap.dependence_plot("As_ppm", shap_values, Xr_test, interaction_index="Te_ppm")

SHAP 依赖图横轴是 As_ppm 值，纵轴是 As 对预测的 SHAP 贡献值，点的颜色是 Te_ppm。造山型金矿数据上常见的模式：As 低于 1000 ppm 时 SHAP 值在 0 附近波动，超过 5000 ppm 后 SHAP 值急剧拉升且颜色变红——说明只有 As 和 Te 同步升高时，Au 预测才被显著推高。这个模式直接用相关矩阵是看不出的，因为相关矩阵只看线性关系，而黄铁矿微量元素的共生是分段的。

特征聚类用层次聚类看哪些元素天然走在一起：

from scipy.cluster.hierarchy import linkage, fcluster from scipy.spatial.distance import squareform corr = df[element_cols].corr().values Z = linkage(squareform(1 - corr), method="average") labels = fcluster(Z, t=0.6, criterion="distance")

阈值 0.6 大致对应距离 0.6 处的聚类数，实际运行时可观察树状图调整。典型的聚类结果会分成三簇：Au-As-Te-Bi-Ag 一簇，反映主成矿热液中的硫化物共沉淀组合；Sb-Pb-Zn 一簇，反映晚期低温脉体；Co-Ni-Cu 一簇，反映黄铁矿的变质沉积成因。这个分组和第 3 节随机森林的特征重要性互为印证，使模型输出具备可解释性。

我的习惯是做完 SHAP 再做特征聚类，两个结果对照着看。如果某元素在两个方法里地位不一致——比如 SHAP 显示重要但聚类归属模糊——回到原始数据做单元素分布检查，通常会找到一批异常测点作祟。复现论文不是终点，把数据、模型和地质解释拧成闭环才是。这篇复现里踩过的坑——检出限填充、分组泄漏、标准化顺序，几乎每个都会在换个数据集后重新出现，把这些细节守住，模型约束的结论才经得起推敲。希望帮到你。

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