数据科学在文物成分分析中的应用：从数据预处理到分类建模

1. 项目背景与核心挑战：当现代分析化学遇见千年文物

去年，我参与了一个非常有意思的交叉学科项目，核心任务是对一批出土的古代玻璃制品进行成分分析，并在此基础上完成鉴别与分类。这听起来像是考古学家的工作，但实际上，它是一次典型的“用数据科学解决人文历史问题”的实战。项目源于一个真实的竞赛题目，其核心挑战在于：我们面对的不是标准化的工业品，而是历经千年埋藏、成分复杂且保存状况各异的文物。如何从一堆看似杂乱的光谱或化学分析数据中，提炼出能够区分玻璃类型、产地甚至年代的有效信息，是贯穿始终的难题。

玻璃，作为一种非晶态硅酸盐材料，其成分直接反映了古代的原料来源、工艺水平和文化交流。例如，高钾玻璃和铅钡玻璃就代表了截然不同的技术体系。我们的工作，就是充当这些“沉默文物”的翻译官，用定量的数据还原定性的历史。整个过程涉及化学分析、数据清洗、特征工程、统计建模和考古学解释等多个环节，任何一个环节的疏漏都可能导致结论的偏差。这不仅考验数据分析的基本功，更考验对领域知识的理解和跨学科思维的融合能力。

2. 数据获取与预处理：从“脏数据”到“干净特征”

拿到手的原始数据，通常来自能量色散X射线荧光光谱（ED-XRF）、电感耦合等离子体质谱（ICP-MS）等无损或微损分析技术。数据表可能包含几十个样本，每个样本对应几十种化学元素（如SiO₂, Na₂O, K₂O, PbO, BaO等）的氧化物含量百分比。第一步的预处理，就充满了“坑”。

2.1 缺失值处理：不是简单的删除或填充

文物数据的缺失值非常普遍。一个元素含量显示为“ND”（未检出），可能意味着三种情况：1）该元素确实不存在；2）含量低于仪器检测限；3）样品该区域风化严重导致元素流失。粗暴地用0或均值填充，会严重扭曲数据分布。

我的处理策略是分层处理：

1. 对于主要成分（如SiO₂, Na₂O, K₂O, CaO等）：若缺失，需结合同类型其他样本和考古背景判断。如果是关键判别元素缺失，该样本在后续建模中可能需要谨慎使用或单独标记。
2. 对于微量元素：若大部分样本都未检出（例如，超过80%的样本某元素含量为ND或0），可以考虑直接删除该特征（元素），因为它可能不包含有效判别信息。
3. 对于部分缺失的微量元素：我常用“检测限的一半”（LOD/2）进行填充。这是一种在环境化学和地球化学中常用的保守估计方法，比用0更合理，因为它承认了该元素可能存在，只是浓度很低。你需要从分析报告中找到或估算出每种元素的检测限。

注意：预处理决策必须有记录。在报告中明确写出：“针对XX元素的缺失，采用LOD/2（具体数值）填充，原因是……”。这体现了分析过程的严谨性。

2.2 成分数据的特殊性：定和约束

玻璃成分数据是典型的“成分数据”。所有氧化物的百分比之和应为100%（或接近100%，因分析误差）。这导致数据处于一个“单纯形”空间，变量之间存在固有的相关性（一个成分增加，其他成分之和必然减少）。直接使用欧氏距离进行聚类或分类可能会得出误导性结果。

解决方案是进行数据转换：

• 对数比转换：这是处理成分数据的标准方法。常用的是中心对数比转换（CLR）或等距对数比转换（ILR）。CLR计算相对简单，即对每个成分取对数后，减去所有成分对数的均值。clr(x) = ln(x_i / g(x))，其中g(x)是所有成分的几何平均数。经过CLR转换后的数据，消除了定和约束，可以在欧氏空间中进行更合理的统计分析。
• 以SiO₂为参比：在考古玻璃分析中，由于SiO₂是玻璃形成体，含量高且相对稳定，有时也采用以SiO₂为分母，计算其他成分与SiO₂的比值，这也能在一定程度上减弱定和约束的影响。

在实际操作中，我通常会准备两套数据：原始百分比数据（用于描述性统计和可视化）和经过CLR转换后的数据（用于后续的多元统计分析建模），并在报告中说明选择依据。

2.3 异常值甄别：是分析误差还是重要线索？

一个样本的PbO含量异常高，是分析错误，还是发现了一件罕见的铅玻璃珍品？异常值处理需要格外谨慎。

我采用的流程是：

1. 可视化初筛：绘制每个元素的箱线图或散点图，直观查看远离主体的样本点。
2. 统计方法确认：使用稳健的统计量，如基于中位数和绝对中位差（MAD）的方法来定义异常值阈值，这比基于均值和标准差的方法对异常值更不敏感。
3. 领域知识复核：这是最关键的一步。将疑似异常样本的所有成分列出，结合考古类型学知识判断。例如，如果一个被初步判定为“高钾玻璃”的样本含有显著高于同类的铅，我们不应轻易将其删除，而应思考：它是否是钾铅玻璃？还是代表了某种过渡类型或外来影响？这个“异常值”本身可能就是最重要的发现。

3. 探索性数据分析与可视化：看见数据背后的故事

在正式建模前，花时间进行探索性数据分析（EDA）至关重要，它能帮助我们形成初步假设，并指导后续的模型选择。

3.1 描述性统计与成分模式识别

首先，按考古学上预设的类别（如“高钾玻璃”、“铅钡玻璃”、“混合类型”或“未知”）分组计算主要成分的均值、标准差、中位数等。制作一个汇总表格：

从这个表格中，我们可以直观看到：K₂O是高钾玻璃的标签性元素，而PbO和BaO则是铅钡玻璃的标签性元素。未知样本的成分值介于两者之间，提示它们可能是混合类型、风化产物或不同的亚类。

3.2 降维可视化：PCA与t-SNE的实战选择

面对十几种元素，我们需要降维到二维或三维来观察样本的整体分布。这里有两个主力工具：

• 主成分分析（PCA）：我最先使用的往往是PCA。它对经过CLR转换后的数据效果很好。PCA能告诉我们哪几个元素（主成分）对数据变异的贡献最大。例如，PC1可能主要由PbO和BaO驱动（代表铅钡玻璃特征），PC2可能主要由K₂O驱动（代表高钾玻璃特征）。在PCA得分图上，如果高钾玻璃和铅钡玻璃能清晰地分居两侧，而未知样本落在中间或形成独立簇，那么我们的分类任务就成功了一大半。

• t-分布随机邻域嵌入（t-SNE）：当PCA图上的分类效果不明显，样本点混杂在一起时，我会尝试t-SNE。t-SNE擅长捕捉复杂的非线性结构，能将局部相似性高的样本聚集得更紧密。但是，t-SNE有一个重大陷阱：它的结果具有随机性，且不能像PCA那样解释每个维度（轴）的物理意义。我的经验是：用PCA做可解释的初步观察，用t-SNE作为辅助验证和寻找复杂簇结构的工具。运行t-SNE时，务必固定随机种子（如random_state=42），以保证结果可复现，并多次调整“困惑度”参数以观察结构的稳定性。

3.3 相关性分析与风化作用考量

绘制所有元素之间的相关性热图。你可能会发现Na₂O和K₂O呈负相关（因为它们是不同的助熔剂），CaO和MgO正相关（可能来自同一类原料，如草木灰）。更重要的是，要关注风化指示元素。

古代玻璃埋藏环境中，碱金属离子（K⁺, Na⁺）容易析出，而土壤中的钙、镁等离子会渗入，发生“离子交换”。这会导致数据中K₂O、Na₂O含量偏低，而CaO、MgO含量偏高。在EDA阶段，如果发现“未知”或“难以分类”的样本普遍具有低钾低钠、高钙高镁的特征，就需要高度怀疑其成分受到了严重风化干扰，而非原始配方。此时，可能需要引入风化校正模型，或使用抗风化元素（如Al₂O₃, Zr等）的比值来进行分析。

4. 分类与鉴别模型构建：从逻辑回归到集成学习

我们的终极目标是建立一个模型，能够根据化学成分准确地将未知玻璃样本归类。这是一个有监督的分类问题。

4.1 特征工程：创造更有判别力的“尺子”

原始元素含量是基础特征，但我们可以创造更有力的衍生特征：

1. 元素比值：如K₂O/PbO， K₂O/BaO， PbO/BaO。这些比值能放大不同类型之间的差异，且对绝对含量的测量误差有一定鲁棒性。
2. 类型指数：定义一个“铅钡玻璃指数” = (PbO + BaO) / (K₂O + 1)，或“高钾玻璃指数” = K₂O / (PbO + BaO + 1)。加1是为了防止分母为零。这些指数能直观地将样本映射到一个标尺上。
3. 风化校正特征：如果确定风化是主要干扰，可以尝试计算“碱金属总量”（R₂O = K₂O + Na₂O）与“碱土金属总量”（RO = CaO + MgO）的比值，或使用Al₂O₃作为内标进行标准化。

4.2 模型选型与实战步骤

我通常会构建一个从简到繁的模型流水线进行对比：

第一步：逻辑回归（LR）或线性判别分析（LDA）

• 为什么先用它们？它们简单、可解释性强。如果数据线性可分，它们的性能不会差，而且LDA本身就是为分类设计的降维方法，我们可以直接观察其判别函数（类似于PCA的主成分）由哪些元素权重构成。
• 实操：将数据（建议用CLR转换后的数据或关键比值特征）分为训练集和测试集（如7:3）。用训练集拟合LR/LDA模型，在测试集上评估准确率、精确率、召回率和F1分数。绘制混淆矩阵，看清模型具体在哪些类别上容易犯错。

第二步：支持向量机（SVM）与核技巧

• 何时使用？当LR/LDA效果不佳，且PCA/t-SNE图显示分类边界可能是非线性时。
• 关键参数：重点是核函数的选择。线性核（kernel='linear'）首先尝试；如果不行，改用径向基函数核（kernel='rbf'），并小心调整惩罚参数C和核系数gamma。过大的gamma会导致过拟合，模型会去拟合每一个训练样本，在未知数据上表现很差。我习惯用网格搜索（GridSearchCV）配合交叉验证来寻找最优参数。

第三步：随机森林（RF）或梯度提升树（如XGBoost）

• 为什么需要树模型？树模型能自动捕捉特征间的交互作用和非线性关系，且对特征的量纲不敏感，无需CLR转换也可以直接使用原始百分比数据（但缺失值需处理）。它们还能给出特征重要性排序，这是一个非常宝贵的副产品，可以告诉我们哪些元素对分类的贡献最大，这与考古学家的经验判断可以相互印证。
• 实操心得：
  • 随机森林的n_estimators（树的数量）通常设置大一些（如200、500），效果会随着树增多而提升并趋于稳定。
  • 控制单棵树的深度（max_depth）以防止过拟合，可以通过交叉验证来调整。
  • 一定要看特征重要性！如果发现某个我们认为是关键判别元素的特征（如PbO）重要性排名很低，就需要回头检查数据或思考我们的领域假设是否正确。

4.3 模型评估与考古学解释的融合

模型准确率达到95%固然好，但更重要的是模型犯错的地方。仔细分析被错误分类的样本：

• 是被分到了另一个主要类别，还是被分到了“未知/其他”类别？
• 这些错误样本的化学成分有什么共同点？是否都表现出强烈的风化特征？
• 从考古学背景看，这些样本是否来自特殊的墓葬、地域或年代？

我曾遇到一个案例：一个本应属于“铅钡玻璃”的样本被模型反复归为“高钾玻璃”。检查数据发现其K₂O含量确实异常高，而PbO含量处于该类别的下限。查阅原始考古记录后发现，该样本出土于一个不同文化交汇地带的遗址。最终，我们并未简单地将此视为模型错误，而是提出了“该样本可能采用了混合配方或受到了外来技术影响”的考古学假设。数据分析的终点不是模型指标，而是产生一个能够被考古学逻辑所检验和解释的合理推论。

5. 聚类分析探索未知：当没有标签时如何发现新类

对于那部分没有任何先验类型信息的“未知”样本，或者当我们想验证已有的分类是否合理时，无监督的聚类分析就派上了用场。

5.1 K-Means与系统聚类法的配合使用

• K-Means：需要指定聚类数K。我们可以借助“肘部法则”绘制不同K值下的误差平方和（SSE）曲线，寻找拐点。但更有效的方法是结合领域知识：如果我们从历史上知道该批文物可能主要来自2-3个技术传统，那么K可以从2、3、4开始尝试。
• 系统聚类（层次聚类）：它的优势是不需要预先指定K，并且可以通过树状图（Dendrogram）直观展示样本逐步聚合的过程。树状图的“纵轴”距离可以帮助我们决定在何处切割以形成类簇。我通常的做法是：先用系统聚类得到一个初步的、可视化的分类趋势，再用这个趋势来指导K-Means中K的选择。

一个实用的流程是：

1. 对“未知”样本数据（或全部样本数据）进行标准化（如Z-score标准化或使用CLR数据）。
2. 绘制系统聚类树状图，观察在某个距离阈值下，样本大致分成了几组。
3. 将这个组数作为K-Means的输入，进行聚类。
4. 比较K-Means的聚类结果与系统聚类在指定K下的结果，看是否一致。
5. 最关键的一步：提取每个簇的化学成分均值剖面，与已知的“高钾玻璃”、“铅钡玻璃”剖面进行对比。看看新发现的簇是更接近其中某一类，还是具有独特的成分特征（例如，中等钾含量、同时含有一定量的铅和钡），从而可能定义出一个新的亚型。

5.2 聚类结果的验证与解读

聚类结果没有千篇一律的正确。需要从多个角度进行验证：

• 内部指标：轮廓系数（Silhouette Score）。值越接近1，说明聚类效果越好。可以计算不同K值下的轮廓系数，辅助选择K。
• 外部指标（如果有部分真实标签）：调整兰德指数（Adjusted Rand Index, ARI）或归一化互信息（NMI）。即使只有部分样本有标签，也可以评估聚类结果与已知标签的一致性。
• 考古学合理性：这是最终的试金石。聚类出的不同群体，是否与出土单位（墓葬、灰坑）、器型、颜色等考古学信息存在关联？如果聚类结果完全无法与任何考古背景信息关联，我们就需要反思特征选择、数据预处理或聚类算法本身是否合适。

6. 综合分析流程与报告撰写：将数据结果转化为考古结论

完成所有分析后，需要将零散的结果整合成一个逻辑连贯、证据链完整的报告。

我的报告结构通常如下：

1. 摘要与问题重述：简明扼要地说明分析了什么、用了什么方法、主要结论是什么。
2. 数据与方法：
   • 数据来源与描述（样本数、元素列表）。
   • 详细的预处理步骤（缺失值、定和约束、异常值处理方案及理由）。
   • 采用的分析方法清单（描述性统计、PCA、t-SNE、逻辑回归、SVM、随机森林、K-Means等）及选用理由。
3. 结果与分析（核心部分）：
   • 整体成分模式：通过描述性统计表和箱线图，展示不同类型玻璃的成分特征。
   • 分类模型结果：以混淆矩阵和性能指标表呈现最佳模型的分类效果。用图表展示模型如何将“未知”样本归类，并列出分类置信度较低的样本以供讨论。
   • 聚类发现：展示树状图和K-Means聚类结果，描述新发现簇的成分特征。
   • 关键判别元素：通过模型的特征重要性排序或LDA的判别函数系数，明确指出哪些元素起到了最关键的分类作用。
4. 讨论：
   • 综合解释分类与聚类结果。例如：“模型成功地将XX号等10个未知样本归入高钾玻璃类，其高K₂O、低PbO/BaO的特征与训练集一致。同时，聚类分析发现第3组样本具有独特的中等铅钡含量，可能代表了本地对铅钡玻璃工艺的一种改良。”
   • 分析错误分类和特殊样本。结合风化迹象、考古背景进行探讨，提出合理的假设。
   • 指出本分析的局限性（如样本量、风化影响、元素分析范围等）。
5. 结论与建议：
   • 给出明确的鉴别结论（每个未知样本的推荐类型）。
   • 提出对古代玻璃工艺源流、文化交流的见解。
   • 建议后续可开展的工作（如扩大样本量、增加同位素分析等）。

在整个项目中，我最大的体会是：技术是手段，人文是灵魂。PCA散点图上的一个点，对应的是古人手中一件真实的器物。我们的每一个参数选择、每一次模型调用，最终都是为了更清晰地听见文物所诉说的历史。保持对数据的严谨，对领域知识的敬畏，在“算力”与“人力”之间找到平衡点，是完成这类交叉学科课题的不二法门。最后一个小建议：所有分析代码（Python/R）必须做好注释，并使用Jupyter Notebook或R Markdown等工具将代码、结果和文字叙述整合在一起，确保整个分析流程完全可复现，这既是专业性的体现，也方便日后回顾或与他人合作。