简介这份资源围绕造山型金矿床黄铁矿微量元素变化这一课题复现了论文《Trace element variations of pyrite in orogenic gold deposits》中的大数据与机器学习分析流程面向地质学、地球化学及数据科学领域的研究人员与机器学习爱好者。内容涵盖数据预处理、PCA与PLS-DA多元统计降维分类、随机森林分类与回归建模以及网格搜索调参与模型评估完整呈现从数据清洗到结果解释的技术链路。资源包为1个docx文档约20KB以文字与Python代码示例为主便于对照学习。已有52人学习下载。读者可借此掌握KNN插补、中心对数比转换、碎石图绘制、随机森林预测金矿化阶段与温度等具体操作并获得可迁移至地质数据分析的建模与排错思路。1. 造山型金矿与黄铁矿微量元素为什么大数据分析和机器学习成了新标配造山型金矿的成因研究里黄铁矿是绕不开的矿物。它几乎贯穿了从早期成矿前到主成矿期再到后期改造的整个演化序列微量元素Au、As、Co、Ni、Sb、Te、Bi 等在黄铁矿中的分布直接记录了流体来源、温度、压力和硫逸度等信息。传统做法是挑几十个点做 LA-ICP-MS 或 EPMA然后画几张散点图、算几个比值靠经验判断成矿阶段。问题是一个矿床动辄几百上千个分析点多矿床对比时数据量轻松上万人工处理根本扛不住而且不同分析批次之间的系统偏差很容易被忽略。大数据分析加机器学习约束解决的正是这个痛点用聚类和分类算法从高维微量元素数据中自动识别成矿相关群体用特征重要性排序反推哪些元素组合最具判别力再用回归模型定量约束温度或流体演化路径。这套方法适合两类人一是做矿床地球化学研究、手头有大量 EPMA/LA-ICP-MS 数据但不知道怎么深挖的研究生和青年学者二是想用机器学习方法处理地质数据的从业者需要一套能跑通、能复现的代码框架。下面我从数据准备、特征工程、模型选型到避坑把整条链路拆开讲。2. 数据准备与特征工程从原始微量元素表到可训练矩阵2.1 黄铁矿微量元素数据的常见来源与清洗规则做这个方向数据来源无非三类自己测的 LA-ICP-MS 时间分辨数据、已发表论文的补充材料表格、以及公开数据库如 GEOROC 中部分硫化物数据。自己测的数据最干净但要注意时间分辨信号中包裹体尖峰的剔除文献数据最大的坑是不同实验室的检出限和标准化流程不同直接合并会导致系统偏差。清洗规则我一般按这几步走第一剔除总和小于 95% 或大于 105% 的 EPMA 数据点这通常是分析误差或表面污染第二对 LA-ICP-MS 数据用 Si 或 Fe 作为内标归一化后再剔除信号中明显异常尖峰通常对应微米级包裹体第三对每个元素若某分析点低于检出限用该元素检出限的一半替代而不是直接删掉——删掉会丢失成矿信息替代虽然保守但保留样本量。import pandas as pd import numpy as np # 读取原始数据假设列名为元素名行为分析点 df pd.read_csv(pyrite_trace_elements.csv) # 剔除总和异常点EPMA 数据常用 element_cols [Au, As, Co, Ni, Sb, Te, Bi, Se, Ag, Cu] df[total] df[element_cols].sum(axis1) df df[(df[total] 95) (df[total] 105)] # 对低于检出限的值用 1/2 DL 替代 detection_limits {Au: 0.01, As: 5, Co: 0.5, Ni: 0.5, Sb: 0.05, Te: 0.02, Bi: 0.02, Se: 1, Ag: 0.05, Cu: 0.5} for elem, dl in detection_limits.items(): df[elem] df[elem].apply(lambda x: dl/2 if x dl else x) # 对数变换微量元素常呈对数正态分布 for elem in element_cols: df[flog_{elem}] np.log10(df[elem] 1e-6) print(df[[log_Au, log_As, log_Co, log_Ni]].describe())这段代码的逻辑很直接先按总和过滤掉明显不可靠的分析点再用检出限的一半替代低于检出限的值最后做对数变换。参数上总和阈值 95~105 适用于 EPMA 数据如果是 LA-ICP-MS 数据这个范围可以放宽到 90~110因为 LA-ICP-MS 的绝对含量精度本身不如 EPMA。对数变换的底数用 10 还是自然对数不影响后续聚类结果但用 10 更符合地球化学惯例方便和文献对比。2.2 特征构造比值、组合与降维的取舍原始元素浓度直接扔进模型不是不行但地球化学里很多信息藏在比值里。比如 Co/Ni 比值是经典的黄铁矿成因判别指标Au/As 比值能反映金的赋存状态Sb/Te 比值对温度敏感。我一般会构造三类特征第一类是元素浓度的对数变换值第二类是成对元素的比值对数比如 log(Co/Ni)、log(Au/As)、log(Sb/Te)第三类是组合特征比如 AsAuSb 的总和代表成矿元素富集程度。构造完特征后面临降维问题。如果特征数超过 20而样本量只有几百直接上聚类会过拟合。常见做法是用主成分分析PCA先降到 5~8 个主成分再输入聚类算法。但 PCA 的缺点是解释性差主成分的地质意义不直观。我的折中方案是先用随机森林做特征重要性排序保留重要性前 10 的特征再输入聚类。这样既降了维又保留了地质可解释性。from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler # 构造比值特征 df[log_Co_Ni] df[log_Co] - df[log_Ni] df[log_Au_As] df[log_Au] - df[log_As] df[log_Sb_Te] df[log_Sb] - df[log_Te] df[log_As_Au_Sb] np.log10(df[As] df[Au] df[Sb] 1e-6) feature_cols [flog_{e} for e in element_cols] [log_Co_Ni, log_Au_As, log_Sb_Te, log_As_Au_Sb] X df[feature_cols].values # 标准化 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 用随机森林评估特征重要性这里假设有成矿阶段标签若无标签可用无监督方法 # 若无标签可跳过此步直接用全部特征做聚类 rf RandomForestClassifier(n_estimators200, random_state42) # 假设 df[stage] 是成矿阶段标签0 为成矿前1 为主成矿期2 为成矿后 rf.fit(X_scaled, df[stage]) importances pd.Series(rf.feature_importances_, indexfeature_cols).sort_values(ascendingFalse) print(importances.head(10))这里随机森林的n_estimators200是经验值特征数在 10~30 之间时200 棵树足够稳定。random_state42是为了复现实际研究中应该跑多次不同随机种子取平均重要性。如果没有成矿阶段标签这一步可以换成无监督的方差过滤或互信息法但解释性会差一些。标准化用StandardScaler而不是归一化是因为微量元素对数变换后近似正态分布标准化更合适。3. 机器学习模型选型与训练聚类、分类与回归的分工3.1 无监督聚类识别成矿群体KMeans 还是 DBSCAN黄铁矿微量元素数据的一个特点是成矿相关群体往往不是均匀分布的而是几个密集簇加上一些散点。KMeans 假设簇是球形且大小相近对黄铁矿数据来说这个假设经常不成立——主成矿期的黄铁矿可能形成一个密集簇而成矿前和成矿后的点更分散。DBSCAN 不需要预设簇数能识别任意形状的簇还能把低密度区域标为噪声更适合这种场景。但 DBSCAN 的调参比 KMeans 麻烦。两个关键参数eps邻域半径和min_samples核心点的最小邻域样本数。我的经验是先对标准化后的数据做 k-距离图找拐点确定eps的初始值然后min_samples设为特征数的 2 倍左右。如果数据量超过 1000min_samples可以适当调大到 10~15避免把噪声当成簇。from sklearn.cluster import DBSCAN, KMeans from sklearn.neighbors import NearestNeighbors import matplotlib.pyplot as plt # 用 k-距离图确定 eps k 5 # 对应 min_samples nbrs NearestNeighbors(n_neighborsk).fit(X_scaled) distances, indices nbrs.kneighbors(X_scaled) distances np.sort(distances[:, k-1], axis0) plt.plot(distances) plt.xlabel(Points sorted by distance) plt.ylabel(f{k}-th nearest neighbor distance) plt.title(k-distance graph for DBSCAN eps selection) plt.show() # 根据拐点选择 eps假设拐点在 1.2 左右 dbscan DBSCAN(eps1.2, min_samples10) labels_db dbscan.fit_predict(X_scaled) df[cluster_db] labels_db # 对比 KMeans kmeans KMeans(n_clusters3, random_state42, n_init10) labels_km kmeans.fit_predict(X_scaled) df[cluster_km] labels_km print(DBSCAN 簇分布, pd.Series(labels_db).value_counts().to_dict()) print(KMeans 簇分布, pd.Series(labels_km).value_counts().to_dict())这段代码先画 k-距离图辅助选eps然后分别跑 DBSCAN 和 KMeans。n_init10是 KMeans 的默认推荐值表示用 10 次不同初始化取最优结果。实际对比时如果 DBSCAN 把大量点标为 -1噪声说明eps太小或min_samples太大如果所有点都在一个簇里说明eps太大。我一般会尝试 3~5 组参数结合地质背景判断哪个结果最合理——比如主成矿期样本应该形成一个独立簇如果聚类结果把主成矿期和成矿前混在一起那参数就需要调整。3.2 有监督分类与回归用随机森林和 XGBoost 约束成矿阶段如果有明确的成矿阶段标签比如通过岩相学观察确定就可以做有监督分类。随机森林和 XGBoost 是地质数据里最常用的两个模型。随机森林对超参数不敏感解释性好能输出特征重要性XGBoost 在样本量较大时精度更高但调参更复杂。我一般先用随机森林做基线如果精度不够再上 XGBoost。回归任务主要是预测温度或流体盐度。黄铁矿微量元素中As 含量和温度有一定相关性Co/Ni 比值也常用来估算温度。但要注意这种回归的 R² 通常不会太高0.5~0.7 就算不错因为微量元素还受其他因素影响。如果看到 R² 超过 0.9大概率是数据泄漏——比如用了温度直接相关的特征或者训练集和测试集有重叠样本。from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split, cross_val_score from sklearn.metrics import classification_report, r2_score import xgboost as xgb # 分类任务预测成矿阶段 X_train, X_test, y_train, y_test train_test_split( X_scaled, df[stage], test_size0.3, random_state42, stratifydf[stage] ) rf_clf RandomForestClassifier(n_estimators300, max_depth8, random_state42) rf_clf.fit(X_train, y_train) y_pred rf_clf.predict(X_test) print(classification_report(y_test, y_pred)) # 回归任务预测温度假设 df[temperature] 存在 X_train_r, X_test_r, y_train_r, y_test_r train_test_split( X_scaled, df[temperature], test_size0.3, random_state42 ) xgb_reg xgb.XGBRegressor(n_estimators300, learning_rate0.05, max_depth5, random_state42) xgb_reg.fit(X_train_r, y_train_r) y_pred_r xgb_reg.predict(X_test_r) print(fXGBoost R²: {r2_score(y_test_r, y_pred_r):.3f}) # 交叉验证 cv_scores cross_val_score(rf_clf, X_scaled, df[stage], cv5, scoringf1_macro) print(fRandom Forest 5-fold CV F1: {cv_scores.mean():.3f} ± {cv_scores.std():.3f})分类部分max_depth8是防止过拟合的常用值如果特征数少可以降到 5~6。stratifydf[stage]保证训练集和测试集的类别比例一致这对不平衡数据很重要。回归部分XGBoost 的learning_rate0.05配合n_estimators300是比较保守的组合不容易过拟合。交叉验证用f1_macro而不是准确率是因为成矿阶段样本往往不平衡准确率会被多数类主导。3.3 模型可解释性SHAP 值如何反推地质意义机器学习模型跑出高精度不是终点地质研究需要知道“为什么”。SHAPSHapley Additive exPlanations是目前最常用的可解释性工具它能给出每个特征对每个样本预测结果的贡献值。对黄铁矿数据来说如果 SHAP 显示 As 和 Co/Ni 比值对主成矿期判别贡献最大那就和传统地质认识一致说明模型学到了真实信号如果 SHAP 显示某个无关元素贡献最大那就要检查数据是否有问题。import shap # 用随机森林模型做 SHAP 分析 explainer shap.TreeExplainer(rf_clf) shap_values explainer.shap_values(X_test) # 对主成矿期假设类别 1画 summary plot shap.summary_plot(shap_values[1], X_test, feature_namesfeature_cols) # 输出前 5 个重要特征 shap_importance np.abs(shap_values[1]).mean(axis0) top_features pd.Series(shap_importance, indexfeature_cols).sort_values(ascendingFalse) print(top_features.head())SHAP 的计算量随样本数增加而增大如果样本超过 5000建议用shap.sample(X_test, 1000)做子采样。TreeExplainer对随机森林和 XGBoost 都适用计算速度比KernelExplainer快很多。解读 SHAP 图时注意看特征值的颜色分布——如果高 As 含量对应正 SHAP 值推动预测为主成矿期那就符合地质预期。4. 避坑与排查黄铁矿微量元素机器学习里最容易翻车的五个地方4.1 数据泄漏标准化在划分训练集之前做了现象模型在测试集上准确率 95% 以上但交叉验证结果波动很大换一批数据就崩。原因先用整个数据集做了标准化或 PCA再划分训练集和测试集导致测试集的信息泄漏到了训练过程中。解决标准化、降维、特征选择都必须在训练集上拟合然后应用到测试集。用sklearn的Pipeline可以避免这个错误。from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.ensemble import RandomForestClassifier pipe Pipeline([ (scaler, StandardScaler()), (clf, RandomForestClassifier(n_estimators300, random_state42)) ]) # 交叉验证时Pipeline 会自动在每折训练集上拟合 scaler cv_scores cross_val_score(pipe, X, df[stage], cv5, scoringf1_macro) print(fPipeline CV F1: {cv_scores.mean():.3f})4.2 类别不平衡主成矿期样本太少导致模型偏向多数类现象分类报告里主成矿期的召回率很低模型几乎把所有样本都预测成成矿前。原因主成矿期样本通常只占 10%~20%模型为了最大化整体准确率倾向于预测多数类。解决用class_weightbalanced让模型自动调整权重或者用 SMOTE 过采样。但 SMOTE 对地球化学数据要慎用因为它会在特征空间插值可能生成地质上不合理的组合。4.3 检出限替代引入的偏差用 1/2 DL 替代后聚类结果偏移现象DBSCAN 把大量低含量样本聚成一簇但这些样本其实分属不同成矿阶段。原因低于检出限的值用 1/2 DL 替代后这些样本在特征空间里被压缩到同一个区域失去了区分度。解决对低于检出限比例超过 30% 的元素考虑剔除该元素或者用 Tobit 回归等专门处理截断数据的方法。如果必须保留可以在聚类时降低这些元素的权重。4.4 过拟合特征数远大于样本量时随机森林也会翻车现象训练集准确率 100%测试集不到 60%。原因特征数 30样本量只有 200随机森林的每棵树都能完美拟合训练集。解决先做特征选择保留重要性前 10 的特征或者用正则化更强的模型如带 L1 正则的逻辑回归增加样本量永远是最有效的但在无法增加的情况下用 5 折交叉验证的 F1 而不是准确率来评估。4.5 地质意义缺失模型精度高但无法解释成矿机制现象审稿人问“你的模型学到了什么地质信息”答不上来。原因只关注了精度指标没有做可解释性分析。解决SHAP 分析是必须的另外可以把聚类结果投影回主成分空间看不同簇在 PC1-PC2 图上的分布是否和成矿阶段对应。如果模型用的特征和传统地质指标Co/Ni、Au/As一致那解释起来就有底气。5. 从复现到进阶用 SHAP 依赖图约束温度区间的实操技巧跑通前面的流程后最有价值的进阶用法是用 SHAP 依赖图做定量约束。具体做法是对温度回归模型画 As 含量和 Co/Ni 比值的 SHAP 依赖图看这两个特征在什么区间内对温度预测的贡献最大。如果 SHAP 显示 As 含量在 100~1000 ppm 区间内对高温预测贡献显著那就可以把这个区间作为高温成矿的判别指标。import shap import matplotlib.pyplot as plt # 用 XGBoost 回归模型做 SHAP 依赖图 explainer_reg shap.TreeExplainer(xgb_reg) shap_values_reg explainer_reg.shap_values(X_test_r) # 画 As 含量的 SHAP 依赖图 shap.dependence_plot(log_As, shap_values_reg, X_test_r, feature_namesfeature_cols, interaction_indexlog_Co_Ni) # 画 Co/Ni 比值的 SHAP 依赖图 shap.dependence_plot(log_Co_Ni, shap_values_reg, X_test_r, feature_namesfeature_cols, interaction_indexlog_As)interaction_index参数指定用哪个特征着色这样能同时看到两个特征的交互效应。如果 As 和 Co/Ni 的交互在图上表现为明显的分群说明这两个指标联合约束温度比单独用更可靠。我一般会把 SHAP 依赖图和传统温度计如 As 温度计的标定曲线叠在一起看如果趋势一致说明机器学习模型没有偏离地质认识如果趋势相反那就要检查数据里是否有混入不同矿床的样本。还有一个实用技巧用 SHAP 值做样本筛选。如果某个样本的 SHAP 值在所有特征上都接近 0说明模型对这个样本的预测没有把握可能是分析误差较大的点可以考虑剔除后重新训练。这个操作要谨慎剔除比例不要超过 5%否则会引入选择偏差。我自己在这个方向踩过最大的坑是早期直接用原始浓度做聚类结果 Co 和 Ni 因为含量高、方差大主导了整个聚类结果As 和 Au 的信息完全被淹没。后来改成先做对数变换再标准化聚类结果才和岩相学观察对上了。另一个习惯是每次跑完模型都会把 SHAP 图和传统地球化学判别图如 Co/Ni vs As并排看如果两者矛盾优先信传统图因为那是几十年积累的经验。希望帮到你。本文还有配套的精品资源点击获取 造山型金矿黄铁矿微量元素：大数据与机器学习分析实战
简介这份资源围绕造山型金矿床黄铁矿微量元素变化这一课题复现了论文《Trace element variations of pyrite in orogenic gold deposits》中的大数据与机器学习分析流程面向地质学、地球化学及数据科学领域的研究人员与机器学习爱好者。内容涵盖数据预处理、PCA与PLS-DA多元统计降维分类、随机森林分类与回归建模以及网格搜索调参与模型评估完整呈现从数据清洗到结果解释的技术链路。资源包为1个docx文档约20KB以文字与Python代码示例为主便于对照学习。已有52人学习下载。读者可借此掌握KNN插补、中心对数比转换、碎石图绘制、随机森林预测金矿化阶段与温度等具体操作并获得可迁移至地质数据分析的建模与排错思路。1. 造山型金矿与黄铁矿微量元素为什么大数据分析和机器学习成了新标配造山型金矿的成因研究里黄铁矿是绕不开的矿物。它几乎贯穿了从早期成矿前到主成矿期再到后期改造的整个演化序列微量元素Au、As、Co、Ni、Sb、Te、Bi 等在黄铁矿中的分布直接记录了流体来源、温度、压力和硫逸度等信息。传统做法是挑几十个点做 LA-ICP-MS 或 EPMA然后画几张散点图、算几个比值靠经验判断成矿阶段。问题是一个矿床动辄几百上千个分析点多矿床对比时数据量轻松上万人工处理根本扛不住而且不同分析批次之间的系统偏差很容易被忽略。大数据分析加机器学习约束解决的正是这个痛点用聚类和分类算法从高维微量元素数据中自动识别成矿相关群体用特征重要性排序反推哪些元素组合最具判别力再用回归模型定量约束温度或流体演化路径。这套方法适合两类人一是做矿床地球化学研究、手头有大量 EPMA/LA-ICP-MS 数据但不知道怎么深挖的研究生和青年学者二是想用机器学习方法处理地质数据的从业者需要一套能跑通、能复现的代码框架。下面我从数据准备、特征工程、模型选型到避坑把整条链路拆开讲。2. 数据准备与特征工程从原始微量元素表到可训练矩阵2.1 黄铁矿微量元素数据的常见来源与清洗规则做这个方向数据来源无非三类自己测的 LA-ICP-MS 时间分辨数据、已发表论文的补充材料表格、以及公开数据库如 GEOROC 中部分硫化物数据。自己测的数据最干净但要注意时间分辨信号中包裹体尖峰的剔除文献数据最大的坑是不同实验室的检出限和标准化流程不同直接合并会导致系统偏差。清洗规则我一般按这几步走第一剔除总和小于 95% 或大于 105% 的 EPMA 数据点这通常是分析误差或表面污染第二对 LA-ICP-MS 数据用 Si 或 Fe 作为内标归一化后再剔除信号中明显异常尖峰通常对应微米级包裹体第三对每个元素若某分析点低于检出限用该元素检出限的一半替代而不是直接删掉——删掉会丢失成矿信息替代虽然保守但保留样本量。import pandas as pd import numpy as np # 读取原始数据假设列名为元素名行为分析点 df pd.read_csv(pyrite_trace_elements.csv) # 剔除总和异常点EPMA 数据常用 element_cols [Au, As, Co, Ni, Sb, Te, Bi, Se, Ag, Cu] df[total] df[element_cols].sum(axis1) df df[(df[total] 95) (df[total] 105)] # 对低于检出限的值用 1/2 DL 替代 detection_limits {Au: 0.01, As: 5, Co: 0.5, Ni: 0.5, Sb: 0.05, Te: 0.02, Bi: 0.02, Se: 1, Ag: 0.05, Cu: 0.5} for elem, dl in detection_limits.items(): df[elem] df[elem].apply(lambda x: dl/2 if x dl else x) # 对数变换微量元素常呈对数正态分布 for elem in element_cols: df[flog_{elem}] np.log10(df[elem] 1e-6) print(df[[log_Au, log_As, log_Co, log_Ni]].describe())这段代码的逻辑很直接先按总和过滤掉明显不可靠的分析点再用检出限的一半替代低于检出限的值最后做对数变换。参数上总和阈值 95~105 适用于 EPMA 数据如果是 LA-ICP-MS 数据这个范围可以放宽到 90~110因为 LA-ICP-MS 的绝对含量精度本身不如 EPMA。对数变换的底数用 10 还是自然对数不影响后续聚类结果但用 10 更符合地球化学惯例方便和文献对比。2.2 特征构造比值、组合与降维的取舍原始元素浓度直接扔进模型不是不行但地球化学里很多信息藏在比值里。比如 Co/Ni 比值是经典的黄铁矿成因判别指标Au/As 比值能反映金的赋存状态Sb/Te 比值对温度敏感。我一般会构造三类特征第一类是元素浓度的对数变换值第二类是成对元素的比值对数比如 log(Co/Ni)、log(Au/As)、log(Sb/Te)第三类是组合特征比如 AsAuSb 的总和代表成矿元素富集程度。构造完特征后面临降维问题。如果特征数超过 20而样本量只有几百直接上聚类会过拟合。常见做法是用主成分分析PCA先降到 5~8 个主成分再输入聚类算法。但 PCA 的缺点是解释性差主成分的地质意义不直观。我的折中方案是先用随机森林做特征重要性排序保留重要性前 10 的特征再输入聚类。这样既降了维又保留了地质可解释性。from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler # 构造比值特征 df[log_Co_Ni] df[log_Co] - df[log_Ni] df[log_Au_As] df[log_Au] - df[log_As] df[log_Sb_Te] df[log_Sb] - df[log_Te] df[log_As_Au_Sb] np.log10(df[As] df[Au] df[Sb] 1e-6) feature_cols [flog_{e} for e in element_cols] [log_Co_Ni, log_Au_As, log_Sb_Te, log_As_Au_Sb] X df[feature_cols].values # 标准化 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 用随机森林评估特征重要性这里假设有成矿阶段标签若无标签可用无监督方法 # 若无标签可跳过此步直接用全部特征做聚类 rf RandomForestClassifier(n_estimators200, random_state42) # 假设 df[stage] 是成矿阶段标签0 为成矿前1 为主成矿期2 为成矿后 rf.fit(X_scaled, df[stage]) importances pd.Series(rf.feature_importances_, indexfeature_cols).sort_values(ascendingFalse) print(importances.head(10))这里随机森林的n_estimators200是经验值特征数在 10~30 之间时200 棵树足够稳定。random_state42是为了复现实际研究中应该跑多次不同随机种子取平均重要性。如果没有成矿阶段标签这一步可以换成无监督的方差过滤或互信息法但解释性会差一些。标准化用StandardScaler而不是归一化是因为微量元素对数变换后近似正态分布标准化更合适。3. 机器学习模型选型与训练聚类、分类与回归的分工3.1 无监督聚类识别成矿群体KMeans 还是 DBSCAN黄铁矿微量元素数据的一个特点是成矿相关群体往往不是均匀分布的而是几个密集簇加上一些散点。KMeans 假设簇是球形且大小相近对黄铁矿数据来说这个假设经常不成立——主成矿期的黄铁矿可能形成一个密集簇而成矿前和成矿后的点更分散。DBSCAN 不需要预设簇数能识别任意形状的簇还能把低密度区域标为噪声更适合这种场景。但 DBSCAN 的调参比 KMeans 麻烦。两个关键参数eps邻域半径和min_samples核心点的最小邻域样本数。我的经验是先对标准化后的数据做 k-距离图找拐点确定eps的初始值然后min_samples设为特征数的 2 倍左右。如果数据量超过 1000min_samples可以适当调大到 10~15避免把噪声当成簇。from sklearn.cluster import DBSCAN, KMeans from sklearn.neighbors import NearestNeighbors import matplotlib.pyplot as plt # 用 k-距离图确定 eps k 5 # 对应 min_samples nbrs NearestNeighbors(n_neighborsk).fit(X_scaled) distances, indices nbrs.kneighbors(X_scaled) distances np.sort(distances[:, k-1], axis0) plt.plot(distances) plt.xlabel(Points sorted by distance) plt.ylabel(f{k}-th nearest neighbor distance) plt.title(k-distance graph for DBSCAN eps selection) plt.show() # 根据拐点选择 eps假设拐点在 1.2 左右 dbscan DBSCAN(eps1.2, min_samples10) labels_db dbscan.fit_predict(X_scaled) df[cluster_db] labels_db # 对比 KMeans kmeans KMeans(n_clusters3, random_state42, n_init10) labels_km kmeans.fit_predict(X_scaled) df[cluster_km] labels_km print(DBSCAN 簇分布, pd.Series(labels_db).value_counts().to_dict()) print(KMeans 簇分布, pd.Series(labels_km).value_counts().to_dict())这段代码先画 k-距离图辅助选eps然后分别跑 DBSCAN 和 KMeans。n_init10是 KMeans 的默认推荐值表示用 10 次不同初始化取最优结果。实际对比时如果 DBSCAN 把大量点标为 -1噪声说明eps太小或min_samples太大如果所有点都在一个簇里说明eps太大。我一般会尝试 3~5 组参数结合地质背景判断哪个结果最合理——比如主成矿期样本应该形成一个独立簇如果聚类结果把主成矿期和成矿前混在一起那参数就需要调整。3.2 有监督分类与回归用随机森林和 XGBoost 约束成矿阶段如果有明确的成矿阶段标签比如通过岩相学观察确定就可以做有监督分类。随机森林和 XGBoost 是地质数据里最常用的两个模型。随机森林对超参数不敏感解释性好能输出特征重要性XGBoost 在样本量较大时精度更高但调参更复杂。我一般先用随机森林做基线如果精度不够再上 XGBoost。回归任务主要是预测温度或流体盐度。黄铁矿微量元素中As 含量和温度有一定相关性Co/Ni 比值也常用来估算温度。但要注意这种回归的 R² 通常不会太高0.5~0.7 就算不错因为微量元素还受其他因素影响。如果看到 R² 超过 0.9大概率是数据泄漏——比如用了温度直接相关的特征或者训练集和测试集有重叠样本。from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split, cross_val_score from sklearn.metrics import classification_report, r2_score import xgboost as xgb # 分类任务预测成矿阶段 X_train, X_test, y_train, y_test train_test_split( X_scaled, df[stage], test_size0.3, random_state42, stratifydf[stage] ) rf_clf RandomForestClassifier(n_estimators300, max_depth8, random_state42) rf_clf.fit(X_train, y_train) y_pred rf_clf.predict(X_test) print(classification_report(y_test, y_pred)) # 回归任务预测温度假设 df[temperature] 存在 X_train_r, X_test_r, y_train_r, y_test_r train_test_split( X_scaled, df[temperature], test_size0.3, random_state42 ) xgb_reg xgb.XGBRegressor(n_estimators300, learning_rate0.05, max_depth5, random_state42) xgb_reg.fit(X_train_r, y_train_r) y_pred_r xgb_reg.predict(X_test_r) print(fXGBoost R²: {r2_score(y_test_r, y_pred_r):.3f}) # 交叉验证 cv_scores cross_val_score(rf_clf, X_scaled, df[stage], cv5, scoringf1_macro) print(fRandom Forest 5-fold CV F1: {cv_scores.mean():.3f} ± {cv_scores.std():.3f})分类部分max_depth8是防止过拟合的常用值如果特征数少可以降到 5~6。stratifydf[stage]保证训练集和测试集的类别比例一致这对不平衡数据很重要。回归部分XGBoost 的learning_rate0.05配合n_estimators300是比较保守的组合不容易过拟合。交叉验证用f1_macro而不是准确率是因为成矿阶段样本往往不平衡准确率会被多数类主导。3.3 模型可解释性SHAP 值如何反推地质意义机器学习模型跑出高精度不是终点地质研究需要知道“为什么”。SHAPSHapley Additive exPlanations是目前最常用的可解释性工具它能给出每个特征对每个样本预测结果的贡献值。对黄铁矿数据来说如果 SHAP 显示 As 和 Co/Ni 比值对主成矿期判别贡献最大那就和传统地质认识一致说明模型学到了真实信号如果 SHAP 显示某个无关元素贡献最大那就要检查数据是否有问题。import shap # 用随机森林模型做 SHAP 分析 explainer shap.TreeExplainer(rf_clf) shap_values explainer.shap_values(X_test) # 对主成矿期假设类别 1画 summary plot shap.summary_plot(shap_values[1], X_test, feature_namesfeature_cols) # 输出前 5 个重要特征 shap_importance np.abs(shap_values[1]).mean(axis0) top_features pd.Series(shap_importance, indexfeature_cols).sort_values(ascendingFalse) print(top_features.head())SHAP 的计算量随样本数增加而增大如果样本超过 5000建议用shap.sample(X_test, 1000)做子采样。TreeExplainer对随机森林和 XGBoost 都适用计算速度比KernelExplainer快很多。解读 SHAP 图时注意看特征值的颜色分布——如果高 As 含量对应正 SHAP 值推动预测为主成矿期那就符合地质预期。4. 避坑与排查黄铁矿微量元素机器学习里最容易翻车的五个地方4.1 数据泄漏标准化在划分训练集之前做了现象模型在测试集上准确率 95% 以上但交叉验证结果波动很大换一批数据就崩。原因先用整个数据集做了标准化或 PCA再划分训练集和测试集导致测试集的信息泄漏到了训练过程中。解决标准化、降维、特征选择都必须在训练集上拟合然后应用到测试集。用sklearn的Pipeline可以避免这个错误。from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.ensemble import RandomForestClassifier pipe Pipeline([ (scaler, StandardScaler()), (clf, RandomForestClassifier(n_estimators300, random_state42)) ]) # 交叉验证时Pipeline 会自动在每折训练集上拟合 scaler cv_scores cross_val_score(pipe, X, df[stage], cv5, scoringf1_macro) print(fPipeline CV F1: {cv_scores.mean():.3f})4.2 类别不平衡主成矿期样本太少导致模型偏向多数类现象分类报告里主成矿期的召回率很低模型几乎把所有样本都预测成成矿前。原因主成矿期样本通常只占 10%~20%模型为了最大化整体准确率倾向于预测多数类。解决用class_weightbalanced让模型自动调整权重或者用 SMOTE 过采样。但 SMOTE 对地球化学数据要慎用因为它会在特征空间插值可能生成地质上不合理的组合。4.3 检出限替代引入的偏差用 1/2 DL 替代后聚类结果偏移现象DBSCAN 把大量低含量样本聚成一簇但这些样本其实分属不同成矿阶段。原因低于检出限的值用 1/2 DL 替代后这些样本在特征空间里被压缩到同一个区域失去了区分度。解决对低于检出限比例超过 30% 的元素考虑剔除该元素或者用 Tobit 回归等专门处理截断数据的方法。如果必须保留可以在聚类时降低这些元素的权重。4.4 过拟合特征数远大于样本量时随机森林也会翻车现象训练集准确率 100%测试集不到 60%。原因特征数 30样本量只有 200随机森林的每棵树都能完美拟合训练集。解决先做特征选择保留重要性前 10 的特征或者用正则化更强的模型如带 L1 正则的逻辑回归增加样本量永远是最有效的但在无法增加的情况下用 5 折交叉验证的 F1 而不是准确率来评估。4.5 地质意义缺失模型精度高但无法解释成矿机制现象审稿人问“你的模型学到了什么地质信息”答不上来。原因只关注了精度指标没有做可解释性分析。解决SHAP 分析是必须的另外可以把聚类结果投影回主成分空间看不同簇在 PC1-PC2 图上的分布是否和成矿阶段对应。如果模型用的特征和传统地质指标Co/Ni、Au/As一致那解释起来就有底气。5. 从复现到进阶用 SHAP 依赖图约束温度区间的实操技巧跑通前面的流程后最有价值的进阶用法是用 SHAP 依赖图做定量约束。具体做法是对温度回归模型画 As 含量和 Co/Ni 比值的 SHAP 依赖图看这两个特征在什么区间内对温度预测的贡献最大。如果 SHAP 显示 As 含量在 100~1000 ppm 区间内对高温预测贡献显著那就可以把这个区间作为高温成矿的判别指标。import shap import matplotlib.pyplot as plt # 用 XGBoost 回归模型做 SHAP 依赖图 explainer_reg shap.TreeExplainer(xgb_reg) shap_values_reg explainer_reg.shap_values(X_test_r) # 画 As 含量的 SHAP 依赖图 shap.dependence_plot(log_As, shap_values_reg, X_test_r, feature_namesfeature_cols, interaction_indexlog_Co_Ni) # 画 Co/Ni 比值的 SHAP 依赖图 shap.dependence_plot(log_Co_Ni, shap_values_reg, X_test_r, feature_namesfeature_cols, interaction_indexlog_As)interaction_index参数指定用哪个特征着色这样能同时看到两个特征的交互效应。如果 As 和 Co/Ni 的交互在图上表现为明显的分群说明这两个指标联合约束温度比单独用更可靠。我一般会把 SHAP 依赖图和传统温度计如 As 温度计的标定曲线叠在一起看如果趋势一致说明机器学习模型没有偏离地质认识如果趋势相反那就要检查数据里是否有混入不同矿床的样本。还有一个实用技巧用 SHAP 值做样本筛选。如果某个样本的 SHAP 值在所有特征上都接近 0说明模型对这个样本的预测没有把握可能是分析误差较大的点可以考虑剔除后重新训练。这个操作要谨慎剔除比例不要超过 5%否则会引入选择偏差。我自己在这个方向踩过最大的坑是早期直接用原始浓度做聚类结果 Co 和 Ni 因为含量高、方差大主导了整个聚类结果As 和 Au 的信息完全被淹没。后来改成先做对数变换再标准化聚类结果才和岩相学观察对上了。另一个习惯是每次跑完模型都会把 SHAP 图和传统地球化学判别图如 Co/Ni vs As并排看如果两者矛盾优先信传统图因为那是几十年积累的经验。希望帮到你。本文还有配套的精品资源点击获取