简介复现论文《Trace element variations of pyrite in orogenic gold deposits》的 Python 实现资料面向地质学家、数据科学家及机器学习研究人员旨在借助大数据分析与机器学习手段揭示造山型金矿床中黄铁矿微量元素与金矿化阶段、温度之间的内在关联。内容涵盖数据清洗与 KNN 插补、中心对数比clr转换、主成分分析PCA与偏最小二乘判别分析PLS-DA、随机森林分类与回归、网格搜索参数优化以及模型评估图表与指标解读每一部分均配有可运行的 Python 代码及逐行解释。资源包仅含 1 个 docx 文档约 20KB以整理好的文档形式呈现完整分析流程与代码注释方便直接对照学习。读者按照文档操作即可掌握从原始数据预处理到机器学习建模评估的完整链路并能结合实际数据灵活调整有效提升复现论文实验的效率与可信度。该资源已有 52 人学习适合需要系统掌握地质大数据分析流程的科研工作者参考使用。1. 造山型金矿黄铁矿微量元素的大数据分析这套数据到底能看出什么黄铁矿是造山型金矿里最普适的载金矿物它的微量元素组合As-Sb-Tl-Au 与 Co-Ni 的相对变化记录了从成矿流体冷却到围岩混染的完整过程。过去我们习惯拿 Au-As 散点图和判别三角图讲故事一张图配一个解释但样品一多、元素一全手工作图就撑不住了。这份资源是一篇论文的完整复现核心思路是把大数据分析里的机器学习方法直接用到几百个 LA-ICP-MS 测点上——数据预处理、PCA、随机森林、PLS-DA 每一步都有可运行代码和逐行解释。适合做矿床学研究的硕博生、想从散点图切换到脚本化分析的地质工程师以及所有被高维微量元素数据卡住的人。2. 数据预处理把 LA-ICP-MS 原始数据变成能喂进模型的特征矩阵2.1 原始数据长什么样元素列、样品列、检出限的天然坑电子探针和 LA-ICP-MS 导出的数据表通常是一个测点一行、一个元素一列单位是 ppm 或 wt%旁边还挂着一列样品编号和矿床类型。听起来很简单但真正动手时你会发现三种脏数据混在一起仪器没测到的点填 0低于检出限的记成0.05这样的字符串脉络不清晰的点直接留空。如果不做检查直接喂给 PCA结果会非常难看。import pandas as pd import numpy as np # 原始数据通常是每个测点一行每个元素一列单位 ppm df pd.read_excel(pyrite_trace_elements.xlsx, sheet_nameLA-ICPMS) # 把样品编号和分组信息剥离开只留数值列进模型 meta_cols [Sample, Point, Type] X_raw df.drop(columnsmeta_cols) # 逐一检查每个元素的零值、负值和缺失值这一步不能省 zero_counts (X_raw 0).sum() neg_counts (X_raw 0).sum() print(零值统计:\n, zero_counts[zero_counts 0]) print(负值统计:\n, neg_counts[neg_counts 0]) # 看偏度偏度绝对值超过 2 的元素基本都需要做对数变换 skewness X_raw.skew() print(skewness[skewness.abs() 2])逻辑说明先把元数据与数值剥离开因为 Sample 和 Type 永远不会进模型留着反而可能被 pandas 当数值列处理。零值统计和负值统计是地质数据的“体检报告”——零值意味着仪器没打到信号或者浓度低于检出限负值则说明基线校正出了问题后者往往要回到原始谱图重新处理。参数说明sheet_name不一定是 LA-ICPMS按你自己的 Excel 结构调整如果是 CSV 文件就把read_excel换成read_csv注意中文表头时的编码问题。meta_cols列表也要对应你表里的实际列名这里只是一个通用模板。2.2 log 变换加小常数为什么直接标准化会翻车黄铁矿微量元素的浓度跨度极大Au 可能只有 ppb 级而 Fe 是 wt% 级就算只看微量元素As 上千 ppm 的同时 In 可能只有零点几。如果跳过对数变换直接把原始值扔给 StandardScaler高含量元素会主导整个 PCA 的方差结构低含量但地质意义重要的元素Te、Bi、Tl会被彻底淹没。这不是玄学是数据尺度问题。# 检出限以下的值通常记作 0.05 或 0.00先统一替换为缺失再填 LOD/2 lod_dict {As: 0.05, Sb: 0.03, Te: 0.02, Au: 0.01} # 每个元素自己的检出限 for col in X_raw.columns: lod lod_dict.get(col, 0.01) # 没列出的元素给一个保守默认值 # 小于等于 0 或缺失的替换成 LOD 的一半 mask X_raw[col].isna() | (X_raw[col] 0) X_raw.loc[mask, col] lod / 2 # 对跨数量级的元素做 log10 变换把偏态分布拉回近似正态 X_log np.log10(X_raw 1e-9) # 加极小值兜底防止出现 log(0) # 变换后再看一眼偏度应该大幅下降 print(X_log.skew().abs().max())逻辑说明LOD/2是处理截尾数据的标准做法把低于检出限的值当成“存在但测不准”而不是“不存在”避免特征分布被一堆零拉偏。log10是成分数据最常见的变换方式它把乘法关系变成加法关系正好对应微量元素之间的稀释和富集过程。后面的1e-9只是兜底真正治好零值问题还是要靠前面的 LOD 替换。参数说明lod_dict里的数值来自仪器报告每台 LA-ICP-MS 都不一样不能照抄最稳妥的做法是去原始数据表头里找每个元素的 LOD 行。log 用 10 还是用 e 对 PCA 结果影响很小但论文的方法部分一定要写清楚审稿人最喜欢追这种细节。2.3 离群值处理截断还是保留这是个地质问题微量元素数据里经常出现 Bi 突然冲到 1000 ppm 这种极端点原因可能是一个微小的黄铜矿包裹体被打进去了也可能是热液叠加的晚期阶段确实富集了 Bi。这两个解释的地质意义截然不同所以我不建议直接删点而是先压住它的影响。from sklearn.preprocessing import StandardScaler, RobustScaler # 按元素做 IQR 截断只把上下 3 倍 IQR 的极端值拉回边界 def clip_iqr(df, k3.0): df_clipped df.copy() for col in df.columns: q1, q3 df[col].quantile(0.25), df[col].quantile(0.75) iq q3 - q1 lo, hi q1 - k * iq, q3 k * iq df_clipped[col] df[col].clip(lo, hi) # clip 是截断不是删除 return df_clipped X_clip clip_iqr(X_log, k3.0) # 看截断后偏度是否可控决定用 StandardScaler 还是 RobustScaler if X_clip.skew().abs().max() 1: scaler RobustScaler(quantile_range(10.0, 90.0)) else: scaler StandardScaler() X_scaled scaler.fit_transform(X_clip) X_scaled pd.DataFrame(X_scaled, columnsX_clip.columns)逻辑说明IQR 截断是“压低极端值”而不是“删掉样品”因为热液叠加产生的点往往携带真实的地质过程信息直接删除会丢失那期成矿事件的记录。RobustScaler 用分位数做缩放对残留离群值不敏感适合截断后仍然偏态的数据。参数说明k3.0是经验值样品干净时取 3混入包裹体时我会降到 2.5。注意标准化只对 PCA、PLS-DA 这类基于距离的算法是必需的后面如果直接跑随机森林树模型对尺度不敏感标准化反而无所谓所以要按算法决定是否做这一步。3. PCA 降维主成分载荷与黄铁矿微量元素地球化学指纹3.1 为什么用 PCA 而不是直接画 Au vs As 散点图微量元素的变量之间高度相关——As、Sb、Tl 常常一起升高Co、Ni 也倾向于同步变化。十几个元素两两画散点图会有几十张图每张图都只讲一个局部故事。PCA 把协方差结构压缩成几个相互正交的主成分每个主成分都是一组元素组合相当于把几十张散点图的共识提炼到一张图上。from sklearn.decomposition import PCA import matplotlib.pyplot as plt pca PCA(n_componentsmin(10, X_scaled.shape[1])) pca_scores pca.fit_transform(X_scaled) expl pca.explained_variance_ratio_ cumsum np.cumsum(expl) n_pc np.argmax(cumsum 0.75) 1 print(f达到 75% 方差需要前 {n_pc} 个主成分) # 载荷矩阵每个主成分里各元素的贡献方向和大小 loadings pd.DataFrame( pca.components_.T, indexX_scaled.columns, columns[fPC{i} for i in range(1, pca.n_components_ 1)] ) # 样品投影到 PC1-PC2 平面按类型着色 fig, ax plt.subplots(figsize(8, 6)) for typ in df[Type].unique(): mask df[Type].values typ ax.scatter(pca_scores[mask, 0], pca_scores[mask, 1], labeltyp, alpha0.7) ax.set_xlabel(fPC1 ({expl[0]:.1%})) ax.set_ylabel(fPC2 ({expl[1]:.1%})) ax.legend() plt.show()逻辑说明pca_scores是每个样品在新坐标轴上的位置loadings是原始变量对主成分的贡献。造山型金矿数据的 PC1 载荷里通常是 As、Sb、Tl、Au 同向Co、Ni 反向这种元素组合直接对应流体温度从高到低的变化是后续地质解释的出发点。参数说明n_components10是一个上限样品数少时 sklearn 会自动限制实际成分数。0.75是我的默认阈值造山型金矿的微量元素数据一般 3 到 5 个主成分就能到 75%如果发现需要 8 个以上说明数据质量有问题先回头查预处理。3.2 主成分数的选择不要只看碎石图拐点碎石图看拐点选主成分数是入门做法但地质数据噪声大拐点经常不明显。更稳的做法是结合累计方差贡献率和平行分析parallel analysis——用随机打乱的同等大小矩阵做 PCA取其特征值作为噪声基线只有真实数据的特征值高于基线时才保留该主成分。from sklearn.utils import resample def parallel_analysis(X, n_iter100, alpha0.95): # 记录真实数据的特征值 real_eigvals PCA().fit(X).explained_variance_ fake_eigvals [] for _ in range(n_iter): X_fake resample(X, replaceFalse) # 逐列打乱破坏相关性模拟纯噪声特征值分布 for col in X_fake.T: np.random.shuffle(col) fake_eigvals.append(PCA().fit(X_fake).explained_variance_) fake_mean np.mean(fake_eigvals, axis0) fake_upper np.quantile(fake_eigvals, alpha, axis0) return real_eigvals, fake_mean, fake_upper real, fake_mean, fake_upper parallel_analysis(X_scaled.values) n_keep np.sum(real fake_upper) print(f平行分析建议保留 {n_keep} 个主成分)逻辑说明平行分析的核心是把每一列数据单独打乱破坏变量间的相关性剩下的特征值就是纯噪声水平。真实数据的特征值高于噪声上限才说明这个维度携带了超出随机水平的结构信息。这个方法比只看碎石图拐点可靠审稿人也认可。参数说明n_iter100是模拟次数越多越稳定但耗时更长几百个样品时 100 次已经很够用。alpha0.95是置信水平取 95% 上分位数做阈值对应显著性检验的直觉。3.3 载荷图的正负方向一个最容易读反的细节PCA 的载荷方向是任意的同一个解乘上 -1 还是同一个解。也就是说 PC1 上 As 的载荷是 0.5 还是 -0.5完全取决于算法初始化的方向不代表地质意义上的正相关或负相关。真正要看的是元素之间的相对方向如果 As 和 Sb 的载荷符号相同说明它们在同一主成分上协同变化如果 Co 和 As 符号相反说明它们在此主成分上呈消长关系。# 以 PC1 载荷为例看的是元素之间的相对关系 pc1 loadings[PC1] print(pc1.sort_values(ascendingFalse)) # 如果 PC1 整体反号可以手动翻转让高载荷元素为正方便解释 if pc1.abs().idxmax() 0: pca.components_[0] * -1 pca_scores[:, 0] * -1逻辑说明翻转符号不会改变样品点之间的相对距离和聚类结构只是让载荷图的方向更符合直觉。我在写论文作图时通常会强制让最重要的元素为正这样读者一眼就能看懂元素组合而不是盯着负号怀疑自己读反了。参数说明这里的pc1.abs().idxmax()是取 PC1 载荷绝对值最大的元素名把它和 0 比较是判断整体符号方向。注意翻转要同步作用在components_和scores上只翻一个会出现投影图与载荷图对不上的问题。4. 随机森林与 PLS-DA从特征重要性到矿床类型判别4.1 随机森林做特征重要性排序谁在真正区分不同成因PCA 是无监督的它只看到数据内部的方差结构不管样品标签。而实际研究中我们往往已经知道每个样品的矿床类型造山型、浅成低温热液型、斑岩型这时候随机森林能回答一个更具体的问题哪些微量元素组合最能区分这些类型。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score from sklearn.inspection import permutation_importance y df[Type].values rf RandomForestClassifier( n_estimators500, max_depth6, min_samples_leaf2, random_state42, n_jobs-1 ) # 先看交叉验证均值训练集准确率没有参考价值 cv_scores cross_val_score(rf, X_scaled, y, cv5, scoringaccuracy) print(f5 折 CV 准确率: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}) # Gini importance 与 permutation importance 对照 rf.fit(X_scaled, y) perm permutation_importance(rf, X_scaled, y, n_repeats20, random_state42) imp_df pd.DataFrame({ feature: X_scaled.columns, gini: rf.feature_importances_, perm: perm.importances_mean }).sort_values(perm, ascendingFalse) print(imp_df.head(10))逻辑说明随机森林在这里承担两个职责——对“哪些元素组合能区分不同黄铁矿成因”排序以及对“矿床类型判别”的可行性做预估。Gini importance 计算快但有个已知毛病是会偏向高基数或数值范围大的特征permutation importance 把某一列打乱后重新测准确率下降幅度更贴近真实贡献。参数说明n_estimators500对几百个样品完全够用。max_depth6和min_samples_leaf2是让单棵树别太深地质样品量往往只有几十到几百深度超过 10 基本就是在背样本。random_state42固定是为了复现结果论文里要写清楚。4.2 PLS-DA 做判别有监督降维比 PCA 更能拉开组间差异PCA 不管样品属于哪类投影方向只追求方差最大。PLS-DA 不一样它在降维的同时最大化类别间的分离度相当于把“哪类样品”这个信息直接放进投影方向里。对黄铁矿微量元素这种组间差异可能不大的数据PLS-DA 的判别效果通常比 PCA 后接分类器更直接。from sklearn.cross_decomposition import PLSRegression from sklearn.preprocessing import LabelEncoder, OneHotEncoder from sklearn.model_selection import StratifiedKFold le LabelEncoder() y_enc le.fit_transform(y) y_onehot OneHotEncoder(sparse_outputFalse).fit_transform(y_enc.reshape(-1, 1)) def plsda_cv(X, Y, n_comp, cv5): skf StratifiedKFold(n_splitscv, shuffleTrue, random_state42) accs [] for train_idx, test_idx in skf.split(X, y_enc): pls PLSRegression(n_componentsn_comp) pls.fit(X[train_idx], Y[train_idx]) pred pls.predict(X[test_idx]) pred_class np.argmax(pred, axis1) accs.append(np.mean(pred_class y_enc[test_idx])) return np.mean(accs), np.std(accs) # 从 1 到 6 个成分做网格搜索选出 CV 准确率最高的组合 for nc in range(1, 7): mean_acc, std_acc plsda_cv(X_scaled.values, y_onehot, nc) print(fn_components{nc}: {mean_acc:.3f} ± {std_acc:.3f})逻辑说明PLSRegression 的输出是连续值所以用argmax取最大的那一类作为预测类别。这里每一步都在独立训练集上拟合、在测试集上评估没有信息泄漏。分层 K 折保证每一折里各类样品比例和全量数据一致避免某一类样品恰好全落在测试集里。参数说明n_components从 1 到 6 是经验范围五分类问题一般 4 个成分以内就够了。注意新版本 sklearn 的OneHotEncoder用sparse_outputFalse旧版本是sparseFalse版本差异会直接报错。4.3 VIP 分数从 PLS 模型里提炼元素贡献度交叉验证准确率告诉你“能不能分”但没告诉你“凭什么分”。PLS-DA 的变量投影重要性VIP可以量化每个元素对模型的贡献VIP 大于 1 的元素通常被认为是重要变量这个阈值在文献里很常用。# 用上面选出的最优成分数重新拟合 best_nc 3 # 假设网格搜索得到的最优成分数是 3 pls PLSRegression(n_componentsbest_nc) pls.fit(X_scaled.values, y_onehot) # VIP 公式综合所有成分的权重和解释方差 t pls.x_scores_ w pls.x_weights_ ss np.sum(t**2, axis0) vip np.sqrt(len(X_scaled.columns) * np.sum(ss * (w**2), axis1) / np.sum(ss)) vip_df pd.DataFrame({ feature: X_scaled.columns, VIP: vip }).sort_values(VIP, ascendingFalse) print(vip_df.head(10))逻辑说明VIP 的计算思路是一个元素在某成分里权重高、且该成分解释的方差大那这个元素的重要性就高。它不像随机森林的特征重要性那样依赖打乱数据而是直接从模型参数里推导两组结果可以互相印证。下面这张表是我常用的对照方式排名随机森林 PermutationPLS-DA VIP1AsAs2SbTl3TeSb4CoCo5NiNi参数说明best_nc要替换成上一段网格搜索实际得到的最优值我这里假设是 3。VIP 的值是相对量不同数据集之间不能直接比较但同一数据集内按 1 为阈值筛选是通用的做法。5. 复现避坑指南五个把结果带偏的细节5.1 零值取 log 之后全是 NaNPCA 直接崩溃现象跑np.log10(X)后打印出来一片-inf和NaNPCA 报错说输入包含缺失值。原因原始数据里 0 值没有处理log10(0) 是负无穷pandas 会把它记录为-inf而不是报错后续所有矩阵运算全部失效。这个坑非常隐蔽因为报错信息往往指向 PCA 而不是 log 那一步。解决在 log 之前先统计每个元素的零值数量把零值统一替换成 LOD/2再用np.log10后检查一次np.isfinite。从那以后我做预处理的固定顺序是零值统计 → LOD 替换 → log → 有效性检查四步缺一不可。5.2 检出限以下的值当成缺失值删掉特征分布被扭曲现象某个元素原始数据里 30% 是0.05你把它当成缺失值删掉之后PCA 的载荷图上 Te 和 Bi 的位置完全不符合地质常识。原因LOD不是缺失而是“低于检测能力”意味着元素存在但浓度测不准。全删掉等于系统性地丢掉了低浓度那一端的数据特征分布被人为裁掉了一截方差结构自然失真。解决统一替换为 LOD/2这是地球化学界的常规操作。如果样品量足够大也可以考虑 Tobit 回归或者生存分析这类专门处理截尾数据的方法但对大多数复现场景LOD/2 够用且好写进方法部分。5.3 PCA 载荷符号写反把正相关读成负相关现象第一版图里 As 和 Sb 的载荷符号一正一负你的讨论部分写“As 与 Sb 呈负相关”审稿人质疑与原始数据矛盾。原因PCA 的特征向量方向是任意的同一个解乘上 -1 物理意义完全相同。算法每次运行可能因为数值误差或初始化不同给出相反的符号这不是错误但很容易被解读反。解决解释时永远看元素之间的相对方向不看绝对正负。作图时手动翻转主成分让最重要的元素符号为正并在图注里写清楚“符号已翻转”。这样既能避免误读也能让图更直观。5.4 PLS-DA 不交叉验证就报告 100% 准确率现象训练集上 PLS-DA 判别准确率 100%你激动得差点写进结论换到新数据立刻掉到 55%。原因PLS-DA 是有监督方法它天然会利用类别信息来构造投影方向。如果只用训练集评估模型记住每个样本的标签位置判别率接近 100% 是数学必然不是模型能力。解决所有准确率必须来自分层 K 折交叉验证每一折都重新拟合并预测。我常用 5 折样品少时建议用留一法leave-one-out但要注意留一法方差大结果要报告平均和标准差。5.5 标准化在划分训练集之前做了造成数据泄漏现象交叉验证准确率奇高但换到外部数据集就崩。检查代码发现StandardScaler().fit_transform(X)跑在了train_test_split前面。原因用全量数据的均值和标准差去缩放训练集和测试集测试集的信息已经通过均值和方差渗进了训练过程。这看起来只差一行代码的顺序却会让模型评估结果虚高。解决用 sklearn 的 Pipeline 把标准化和模型串起来让每一折 CV 里只对训练折做 fit再对验证折做 transform。下面的代码是标准写法from sklearn.pipeline import Pipeline pipe Pipeline([ (scale, StandardScaler()), (pls, PLSRegression(n_components3)) ]) cv_scores cross_val_score(pipe, X_scaled, y_onehot, cv5) print(fPipeline CV 准确率: {cv_scores.mean():.3f})逻辑说明Pipeline的核心作用是保证预处理参数只在训练折上估计测试折永远接触不到训练过程中计算出的任何统计量。代码里cross_val_score会把整个 Pipeline 当成一个模型来交叉验证每折自动完成“先标化再 PLS-DA”的完整流程。参数说明n_components还是要回到 4.2 的网格搜索结果来确定Pipeline 只是修数据泄漏不解决超参数选择。6. 完整复现脚本把预处理到模型输出的流程固化下来把前面所有代码串成一个脚本是我拿到任何新数据集都会先搭的骨架。整体流程是读取原始数据 → 零值和 LOD 检查 → log 变换 → IQR 截断 → 标准化 → PCA 投影图 → 随机森林 CV 置换重要性 → PLS-DA CV VIP。这个流程跑完你手里的成果是一张 PC1-PC2 投影图、一份特征重要性表、一份 VIP 表外加一个交叉验证准确率足够支撑一篇短文的核心图件。验证技巧上我习惯把随机森林的置换重要性排名和 PLS-DA 的 VIP 排名放在同一张表里对照两者重合的元素组合就是这篇论文真正要讨论的要素。如果随机森林说 As 最重要而 VIP 说 Te 最重要我倾向于回头检查数据预处理而不是直接采信某一个结果。还有一种常见做法是把 PCA 投影图上明显分群的样品挑出来重新在原始数据里看中位数差异——模型方面的结论最终要能回到原始分析数据上被手动验证这一步别省。我印象最深的一次翻车是一批数据里 Te 在所有分析里高得反常随机森林和 PLS-DA 都把它排在第一图也画得很漂亮。后来核对仪器日志才发现那天标样老化Te 的校正系数偏了整整一个数量级。从那以后我每次动手跑模型之前都强制走一遍元素浓度数量级检查、标样对比和缺失值分布确认模型跑得再快也不如数据本身可靠。希望帮到你。本文还有配套的精品资源点击获取