简介本资源为造山型金矿床黄铁矿微量元素大数据与机器学习研究的复现资料面向地质学、地球化学及数据科学领域的研究人员与研究生帮助其掌握从数据清洗到模型评估的完整分析流程。包内仅含1个docx文档约20KB以文字与代码块形式呈现涵盖数据预处理、PCA与PLS-DA降维分类、随机森林分类与回归建模及网格搜索调参等核心环节。文档基于Python示例具体演示了KNN插补、中心对数比转换、碎石图绘制、金矿化阶段预测与温度回归等操作并附有逐段解释便于读者理解成分数据处理与模型优化的思路。目前已有52人学习适合希望将机器学习方法应用于矿物学与矿床学问题、需要可运行代码参考的读者实际使用时可根据自身数据情况调整参数与流程。1. 从一块黄铁矿说起这套代码到底能帮你跑出什么结果如果你手上有造山型金矿床的 LA-ICP-MS 黄铁矿微量元素数据或者正在做矿床地球化学的数据挖掘这套复现代码值得你花一个下午跑一遍。它做的事情很具体把黄铁矿的微量元素组成当作高维数据矩阵用 PCA 降维看样品分群用随机森林做元素重要性排序和成因判别再用 PLS-DA 做有监督的分类建模最终回答一个矿床学核心问题——黄铁矿的微量元素变化到底记录了什么样的成矿信息。整套流程从数据预处理到模型评估全部可运行不是伪代码也不是只给框架让你自己填。我第一次拿到这类数据的时候最大的困惑不是算法本身而是微量元素有几十个维度检出限以下的值怎么处理不同样品的信号计数怎么归一化PCA 跑出来 PC1 解释了 40% 的方差到底代表什么地质意义这些问题论文里往往一笔带过但代码里必须给出明确答案。这套复现资源的价值就在于它把这些“论文里不写但你必须做”的步骤全部展开了。适合谁看做矿床地球化学的研究生、需要快速上手机器学习的地质从业者、以及想用 Python 把地球化学数据跑出可解释结果的工程师。不适合谁已经有成熟机器学习流水线、只想找一个现成分类器调包的人——这里的重点不是模型多先进而是每一步为什么这么做。2. 数据预处理从原始信号到可建模矩阵的五个关键决策2.1 为什么预处理决定了后面所有结论的可信度黄铁矿微量元素数据有几个天然特点元素种类多通常 2050 个、检出限以下的值大量存在、不同元素的浓度跨度可达 56 个数量级、部分元素之间存在强烈的共生或替代关系。如果你直接把原始浓度矩阵丢进 PCA 或随机森林结果大概率是不可用的——要么被高浓度元素主导要么因为缺失值报错要么分类器学到的是采样偏差而不是地质信号。常见做法是分五步走缺失值处理、数据变换、标准化、异常值筛查、特征筛选。每一步都有多种选择选错了不会报错但结论会偏。下面我按代码实际执行的顺序拆开讲。2.2 缺失值与检出限处理不是简单填个零LA-ICP-MS 数据里低于检出限的值通常标记为 bdlbelow detection limit或负数。直接填零会引入虚假的“低值聚集”直接删样本又会损失大量信息。代码里采用的是替换为检出限的一半1/2 LOD这是地球化学领域的常见做法。import numpy as np import pandas as pd # 读取原始数据假设第一列是样品编号其余列是元素浓度 df pd.read_csv(pyrite_te_raw.csv, index_col0) # 查看缺失和 bdl 情况 print(原始数据形状:, df.shape) print(各列缺失值数量:) print(df.isnull().sum()[df.isnull().sum() 0]) # 将 bdl 标记假设用 -999 表示替换为 NaN df df.replace(-999, np.nan) # 对每个元素用该元素检出限的一半填充 # 这里假设检出限存在一个字典里实际使用时从论文或仪器报告获取 lod_dict { Co: 0.01, Ni: 0.02, As: 0.05, Se: 0.03, Ag: 0.01, Sb: 0.02, Te: 0.01, Au: 0.005, Cu: 0.05, Zn: 0.1, Pb: 0.05, Bi: 0.01 } for elem, lod in lod_dict.items(): if elem in df.columns: df[elem] df[elem].fillna(lod / 2) # 剩余未指定检出限的列用该列最小值的一半填充 df df.fillna(df.min() / 2) print(填充后缺失值总数:, df.isnull().sum().sum())这段代码的逻辑是先识别 bdl 标记再按元素逐个用 1/2 LOD 填充。参数方面lod_dict需要你根据实际仪器的检出限报告来填不同实验室的 LOD 差异很大不能照搬。如果某个元素没有指定 LOD用该列最小值的一半兜底虽然粗糙但比填零合理。注意填充之后一定要检查各元素的分布是否出现异常的尖峰——如果某个元素大量样本都是 bdl填充后会在低值端形成一个假峰这种元素后续要考虑剔除。2.3 对数变换与标准化让不同量级的元素公平竞争微量元素浓度跨度极大Au 可能在 0.001 ppm 级别Fe 可能在百分含量级别。如果不做变换PCA 的第一主成分几乎必然被高浓度元素占据而真正有成因指示意义的低浓度元素如 Au、Ag、Te会被淹没。代码里采用 log10 变换后再做 z-score 标准化。from sklearn.preprocessing import StandardScaler # 对数变换加一个小常数避免 log(0) df_log np.log10(df 1e-6) # 检查变换后是否有无穷值 print(无穷值数量:, np.isinf(df_log).sum().sum()) df_log df_log.replace([np.inf, -np.inf], np.nan).fillna(df_log.min()) # z-score 标准化 scaler StandardScaler() df_scaled pd.DataFrame( scaler.fit_transform(df_log), indexdf_log.index, columnsdf_log.columns ) print(标准化后均值应接近0:, df_scaled.mean().abs().max()) print(标准化后标准差应接近1:, df_scaled.std().mean())log10 变换把乘性关系转为加性关系更适合后续的线性方法。加1e-6是为了避免零值取对数报错这个常数要远小于你的最小非零浓度。标准化让每个元素的权重相同避免量级差异主导。注意如果你后续要用 PLS-DA标准化是必须的如果只用随机森林标准化不影响树的分裂点但做了也没坏处保持流程统一。2.4 异常值筛查别让一个坏点毁掉整个 PCA地球化学数据里经常混入不同期次或不同成因的样品这些样品在 PCA 图上会表现为远离主群的离群点。如果直接纳入建模会严重拉扯主成分方向。代码里用马氏距离Mahalanobis distance做筛查。from scipy.spatial.distance import mahalanobis from numpy.linalg import inv # 计算马氏距离 cov_matrix np.cov(df_scaled.T) inv_cov inv(cov_matrix) mean_vec df_scaled.mean().values mahal_dist df_scaled.apply( lambda row: mahalanobis(row.values, mean_vec, inv_cov), axis1 ) # 以卡方分布 97.5% 分位数为阈值自由度为特征数 from scipy.stats import chi2 threshold chi2.ppf(0.975, dfdf_scaled.shape[1]) outliers mahal_dist[mahal_dist threshold] print(f异常值阈值: {threshold:.2f}) print(f检出异常值数量: {len(outliers)}) print(异常值样品编号:, outliers.index.tolist()) # 标记但不自动删除后续对比分析 df_scaled[is_outlier] (mahal_dist threshold).astype(int)马氏距离考虑了特征之间的协方差比欧氏距离更适合高维相关数据。阈值用卡方分布的 97.5% 分位数自由度等于特征数。这里我建议标记而不是直接删除——先看看这些异常值是不是有地质意义比如不同成矿期次的样品如果确认是分析误差或采样污染再删。参数方面如果你希望更保守可以把分位数调到 99%如果数据本身 heterogeneity 很强调到 95% 也可以。2.5 特征筛选去掉冗余元素保留成因指示信息几十个微量元素里有些是强相关的比如 Co 和 Ni 在黄铁矿中经常协同变化有些则几乎不变。保留所有特征不仅增加计算量还会稀释真正有判别力的信号。代码里用方差膨胀因子VIF和相关性矩阵双重筛选。from statsmodels.stats.outliers_influence import variance_inflation_factor # 计算 VIF feature_cols [c for c in df_scaled.columns if c ! is_outlier] X df_scaled[feature_cols].values vif_data pd.DataFrame() vif_data[feature] feature_cols vif_data[VIF] [ variance_inflation_factor(X, i) for i in range(X.shape[1]) ] # 筛选 VIF 10 的特征 high_vif vif_data[vif_data[VIF] 10][feature].tolist() print(高 VIF 特征考虑剔除:, high_vif) # 相关性矩阵剔除相关系数 0.85 的冗余对中的后者 corr_matrix df_scaled[feature_cols].corr().abs() upper corr_matrix.where( np.triu(np.ones(corr_matrix.shape), k1).astype(bool) ) to_drop [col for col in upper.columns if any(upper[col] 0.85)] print(高相关冗余特征:, to_drop) # 最终保留的特征 final_features [c for c in feature_cols if c not in set(high_vif to_drop)] print(f最终保留特征数: {len(final_features)}) print(保留特征:, final_features)VIF 大于 10 说明该特征与其他特征存在严重共线性考虑剔除。相关系数阈值 0.85 是经验值你可以根据数据调整——如果两个元素地质意义不同但统计相关保留哪个需要结合矿床知识判断。这一步的输出直接决定后续 PCA 和随机森林的输入维度建议把筛选前后的结果都跑一遍对比。3. PCA 与随机森林降维看分群、建模看重要性3.1 PCA 实战从碎石图到双标图的地质解读PCA 在这套流程里的角色是探索性分析——先看看样品自然分不分群分群跟哪些元素有关。代码里用 sklearn 的 PCA但关键不在跑通而在解读。from sklearn.decomposition import PCA import matplotlib.pyplot as plt # 用筛选后的特征做 PCA X_final df_scaled[final_features].values pca PCA(n_componentsmin(len(final_features), 10)) scores pca.fit_transform(X_final) # 碎石图看前几个主成分解释了多少方差 explained pca.explained_variance_ratio_ cumulative np.cumsum(explained) fig, axes plt.subplots(1, 2, figsize(12, 5)) axes[0].bar(range(1, len(explained)1), explained * 100, alpha0.7) axes[0].plot(range(1, len(explained)1), cumulative * 100, ro-) axes[0].set_xlabel(主成分) axes[0].set_ylabel(解释方差 (%)) axes[0].axhline(y80, colorgray, linestyle--, alpha0.5) # 双标图PC1 vs PC2 的样品得分和元素载荷 axes[1].scatter(scores[:, 0], scores[:, 1], csteelblue, alpha0.6) for i, feat in enumerate(final_features): axes[1].arrow(0, 0, pca.components_[0, i] * 3, pca.components_[1, i] * 3, colorcoral, alpha0.7, head_width0.1) axes[1].text(pca.components_[0, i] * 3.2, pca.components_[1, i] * 3.2, feat, fontsize8, colorcoral) axes[1].set_xlabel(fPC1 ({explained[0]*100:.1f}%)) axes[1].set_ylabel(fPC2 ({explained[1]*100:.1f}%)) axes[1].axhline(0, colorgray, linewidth0.5) axes[1].axvline(0, colorgray, linewidth0.5) plt.tight_layout() plt.savefig(pca_biplot.png, dpi300) plt.show() # 输出载荷矩阵看每个主成分主要受哪些元素控制 loadings pd.DataFrame( pca.components_[:3].T, indexfinal_features, columns[PC1, PC2, PC3] ) print(前三个主成分的载荷绝对值排序:) for pc in [PC1, PC2, PC3]: print(f\n{pc} 主要载荷:) print(loadings[pc].abs().sort_values(ascendingFalse).head(5))碎石图告诉你保留几个主成分——通常累计解释方差到 70%80% 即可。双标图是地质解读的核心箭头方向代表元素载荷箭头越长该元素对该主成分贡献越大样品点聚集区域对应的元素箭头方向就是该群的特征元素组合。比如 PC1 正方向如果 As、Au、Ag 载荷高那 PC1 正半轴的样品可能代表富 As-Au 的成矿期次。参数方面n_components先设大一点看碎石图再定箭头缩放系数3是可视化用的不影响实际载荷值。3.2 随机森林元素重要性排序与分类判别随机森林在这套流程里干两件事一是输出元素重要性排序告诉你哪些微量元素对分类贡献最大二是做交叉验证的分类判别评估模型能不能区分不同成因类型的黄铁矿。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import StratifiedKFold, cross_val_score from sklearn.metrics import classification_report, confusion_matrix # 假设有一个成因标签列 genesis_type0沉积期1变质期2岩浆期 # 实际使用时根据你的样品标签替换 labels df_scaled[genesis_type].values # 这里需要你提前准备好标签 X_rf df_scaled[final_features].values # 随机森林建模 rf RandomForestClassifier( n_estimators500, # 树的数量500 通常足够稳定 max_featuressqrt, # 每次分裂考虑 sqrt(特征数) 个特征 min_samples_leaf2, # 叶节点最少样本数防止过拟合 random_state42, class_weightbalanced # 类别不平衡时自动加权 ) # 交叉验证评估 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(rf, X_rf, labels, cvcv, scoringf1_macro) print(f5折交叉验证 F1-macro: {scores.mean():.3f} ± {scores.std():.3f}) # 全量训练后输出特征重要性 rf.fit(X_rf, labels) importance pd.DataFrame({ feature: final_features, importance: rf.feature_importances_ }).sort_values(importance, ascendingFalse) print(\n元素重要性排序前10:) print(importance.head(10).to_string(indexFalse)) # 混淆矩阵 from sklearn.model_selection import cross_val_predict y_pred cross_val_predict(rf, X_rf, labels, cvcv) print(\n混淆矩阵:) print(confusion_matrix(labels, y_pred)) print(\n分类报告:) print(classification_report(labels, y_pred))参数方面n_estimators500是稳定性和计算成本的折中再多提升有限max_featuressqrt是分类任务的默认推荐min_samples_leaf2防止每棵树记住噪声class_weightbalanced在类别不均衡时很重要——地质样品里某一成因类型的样品往往远少于其他类型。交叉验证用 StratifiedKFold 保证每折的类别比例一致。特征重要性输出后前几个元素通常就是判别成因的关键指示元素结合矿床知识判断是否合理。3.3 PLS-DA有监督降维与分类边界可视化PLS-DA 和 PCA 的区别在于它利用了类别标签信息来寻找最大区分度的投影方向。当 PCA 分群不明显但你知道类别标签时PLS-DA 往往能给出更清晰的分类边界。from sklearn.cross_decomposition import PLSRegression from sklearn.preprocessing import LabelBinarizer # 将标签做 one-hot 编码 lb LabelBinarizer() Y_onehot lb.fit_transform(labels) # PLS-DA本质是对 one-hot 标签做 PLS 回归 plsda PLSRegression(n_components3, scaleFalse) # 数据已标准化不再重复 plsda.fit(X_rf, Y_onehot) # 得分图 scores_pls plsda.x_scores_ fig, ax plt.subplots(figsize(8, 6)) scatter ax.scatter(scores_pls[:, 0], scores_pls[:, 1], clabels, cmapviridis, alpha0.7, edgecolorsk) ax.set_xlabel(PLS Component 1) ax.set_ylabel(PLS Component 2) plt.colorbar(scatter, label成因类型) plt.tight_layout() plt.savefig(plsda_scores.png, dpi300) plt.show() # VIP 值Variable Importance in Projection t plsda.x_scores_ w plsda.x_weights_ q plsda.y_loadings_ p X_rf.shape[1] vip np.sqrt(p * np.sum((w**2) * np.sum(t**2, axis0) * np.sum(q**2, axis0), axis1) / np.sum(np.sum(t**2, axis0) * np.sum(q**2, axis0))) vip_df pd.DataFrame({ feature: final_features, VIP: vip }).sort_values(VIP, ascendingFalse) print(VIP 1 的元素对分类有显著贡献:) print(vip_df[vip_df[VIP] 1].to_string(indexFalse))n_components3是先用默认值试再看交叉验证的分类准确率调整。VIP 值大于 1 通常认为该元素对分类有显著贡献这和随机森林的重要性排序可以互相验证——如果两种方法都排在前面的元素可信度更高。注意PLS-DA 容易过拟合样本量小于 30 时慎用或者必须用严格的留一交叉验证来评估。4. 避坑与排查跑这套代码时最容易翻车的五个地方4.1 坑一检出限填充后分布出现假峰现象某个元素在直方图上低值端出现一个异常尖峰log 变换后仍然存在。原因该元素大量样本低于检出限用 1/2 LOD 填充后这些样本全部集中在同一个值上形成人工峰。解决统计每个元素的 bdl 比例如果超过 30% 的样本都是 bdl考虑剔除该元素而不是填充。代码里可以加一行bdl_ratio (df_raw -999).mean()来筛查。4.2 坑二PCA 双标图箭头重叠看不清现象几十个元素的箭头挤在原点附近完全无法解读。原因标准化后所有元素的载荷值都在 -1 到 1 之间低载荷元素的箭头太短。解决只标注载荷绝对值排前 810 的元素其余不画箭头或者用载荷热图替代双标图。代码里把for i, feat in enumerate(final_features)改成只遍历 top 元素即可。4.3 坑三随机森林交叉验证 F1 很高但混淆矩阵很差现象交叉验证 F1-macro 有 0.85但混淆矩阵显示某一类几乎全部被错分。原因类别不均衡导致 F1-macro 被多数类拉高少数类的召回率很低。解决看classification_report里每个类的 precision/recall不要只看总体 F1。如果少数类召回率低于 0.5需要增加该类样本、用 SMOTE 过采样、或者调整class_weight参数。4.4 坑四PLS-DA 的 VIP 值和随机森林重要性排序矛盾现象随机森林说 As 最重要PLS-DA 的 VIP 却把 Co 排第一。原因两种方法的数学原理不同——随机森林基于分裂不纯度减少PLS-DA 基于投影方向上的协方差最大化。它们捕捉的信号角度不一样。解决不要强行统一把两种排序都列出来取交集作为高置信度的指示元素。如果某个元素在两种方法里都排前 5那它大概率是真正有判别力的。4.5 坑五换一批数据后代码报形状错误现象用新数据跑代码时在标准化或 PCA 步骤报ValueError: shapes not aligned。原因新数据的元素列数和列顺序与训练时不一致或者某些元素全部缺失被 pandas 自动丢弃。解决在预处理开头强制对齐列——df df.reindex(columnsexpected_columns, fill_valuenp.nan)然后再做填充。把训练时的final_features列表保存下来预测新数据时按这个列表取列。5. 进阶技巧用置换检验验证模型显著性别让随机森林骗了你随机森林的特征重要性有一个容易被忽略的问题即使标签是随机打乱的它也会给出一个“重要性排序”。换句话说你看到的元素重要性可能只是噪声。要判断模型是否真的学到了信号置换检验permutation test是必须走一遍的。做法很直接把标签随机打乱 1000 次每次重新训练随机森林并记录交叉验证 F1得到一个“随机基线分布”。然后看真实标签的 F1 在这个分布里的位置——如果真实 F1 远高于 95% 分位数说明模型显著如果只是略高那你的分类结果可能不可靠。from sklearn.utils import shuffle n_permutations 1000 perm_scores [] for i in range(n_permutations): labels_shuffled shuffle(labels, random_statei) score cross_val_score( rf, X_rf, labels_shuffled, cvcv, scoringf1_macro ).mean() perm_scores.append(score) perm_scores np.array(perm_scores) real_score cross_val_score( rf, X_rf, labels, cvcv, scoringf1_macro ).mean() p_value (perm_scores real_score).mean() print(f真实 F1: {real_score:.3f}) print(f置换检验 p 值: {p_value:.4f}) print(f随机基线 95% 分位数: {np.percentile(perm_scores, 95):.3f}) if p_value 0.05: print(模型显著分类结果可信) else: print(模型不显著需要检查特征或标签质量)1000 次置换大概需要几分钟到十几分钟取决于样本量和树的数量。如果计算资源有限200 次也能给出粗略判断。p 值小于 0.05 是常规阈值但地质数据样本量往往不大我一般会放宽到 0.1 并结合地质意义判断。另一个实用技巧是把随机森林的feature_importances_和 PLS-DA 的 VIP 值做秩相关分析。如果两种方法的排序显著相关Spearman 相关系数 0.6说明元素的重要性排序是稳健的如果完全不相关那你的结论可能依赖于方法选择需要谨慎。from scipy.stats import spearmanr # 假设 importance 和 vip_df 已经按特征名对齐 merged importance.merge(vip_df, onfeature) rho, p_spearman spearmanr(merged[importance], merged[VIP]) print(f随机森林重要性 vs PLS-DA VIP 的 Spearman 相关系数: {rho:.3f}) print(fp 值: {p_spearman:.4f})从那以后我每次跑完分类模型都强制走一遍置换检验和秩相关分析——这两个步骤花不了多少时间但能帮你避开“模型看起来很好但其实什么都没学到”的坑。希望帮到你。本文还有配套的精品资源点击获取