1. 项目背景与核心挑战当数学建模遇上文物科学去年国赛C题直接把一堆古代玻璃的数据甩到了我们面前。说实话刚拿到题的时候我们队三个人面面相觑这玩意儿跟数学建模有啥关系但仔细一琢磨这恰恰是国赛的魅力所在——它不考你死记硬背的算法而是考你如何用数学工具去解决一个真实的、跨学科的复杂问题。题目要求我们根据古代玻璃文物表面的化学成分数据去分析它们的风化情况、鉴别它们的类型甚至追溯它们的亚类划分。这背后其实是文物保护、考古研究中的一个经典难题如何通过有限的、可能已经发生变化的现代检测数据去反推文物在千百年前的状态和归属。玻璃尤其是中国古代的玻璃成分非常复杂。它不是现代那种纯净的二氧化硅古人会加入草木灰富含钾、硝石引入钠、铅丹、氧化钡等各种助熔剂和着色剂。题目给出的数据就是通过现代仪器比如X射线荧光光谱仪测得的各种氧化物如SiO₂, Na₂O, K₂O, CaO等的百分比含量。但问题来了这些玻璃埋藏地下几百年甚至上千年表面的化学成分早就和内部不一样了这就是“风化”。风化会使得某些成分如K₂O, Na₂O流失而另一些成分如SiO₂相对富集。我们拿到的数据是“风化点”和“未风化点”的混合第一步就得先把它们区分开这本身就是一个基于成分模式的分类问题。更核心的挑战在于题目要求我们对“高钾玻璃”和“铅钡玻璃”这两大类进行亚类划分。这就像给你一堆人脸照片让你不仅分出男女大类还要把同一性别里的人按照某种特征比如脸型、发型再分成几个小群。数据里没有标签全凭数据本身说话。这就需要我们深入理解玻璃的工艺古代工匠的配方不是标准化的不同产地、不同时期的工匠可能用不同比例的铅和钡或者加入不同的铜、铁等着色剂从而形成化学成分特征有微妙差异的亚类。我们的任务就是用数学方法把这些“配方特征”从数据中挖掘出来。所以这个题目的本质是一个无监督学习和有监督学习相结合的数据挖掘问题。前半部分风化判别、类型鉴别可以看作有监督分类如果利用已知部分样本的话或模式识别后半部分亚类划分则是纯粹的无监督聚类分析。而Matlab以其强大的矩阵运算能力、丰富的统计与机器学习工具箱以及出色的可视化功能成为了解决这类问题的绝佳平台。它能让我们的分析过程从数据清洗、特征工程到模型构建、结果可视化形成一个流畅的闭环。2. 数据预处理清洗、变换与特征工程的“炼丹”过程拿到数据通常是附件表单1.xlsx和附件表单2.xlsx后千万别急着跑模型。原始数据里埋着不少“坑”预处理这一步没做好后面所有高级算法得出的结论都可能是在垃圾上建高楼。我们的预处理流程可以概括为“查、补、删、变”四步。2.1 数据探查与缺失值处理首先用Matlab的readtable函数读入数据然后立刻用summary或ismissing函数查看数据概况。玻璃成分数据常见的第一个问题就是缺失值在数据表中可能显示为“NaN”或空值。这些缺失值产生的原因很多比如检测时该成分含量低于仪器检出限或者该成分根本不存在。对于缺失值不能简单地删除整行或整列因为每个样本都极其珍贵。我们采用的策略是分情况填补对于主要成分如SiO₂的缺失这几乎是不可接受的因为它是玻璃的骨架。如果出现需要结合文物编号回溯原始记录或视为无效样本。在本题数据中通常不会出现。对于微量成分如SrO, SnO₂的缺失这些元素含量本身很低缺失很可能意味着“未检出”即含量极低接近于0。因此一个合理且常用的方法是将其填补为0。在Matlab中可以简单地使用data_matrix(isnan(data_matrix)) 0;对于助熔剂成分如K₂O, Na₂O, CaO的缺失需要谨慎。如果同一个大类如所有高钾玻璃中其他样本的该成分都有值唯独某个样本缺失可以考虑用该类别的中位数进行填补这比均值更抗干扰。% 假设 data 是表格K2O是列名Type是玻璃类型列 k_glass_idx strcmp(data.Type, ‘高钾’); k2o_values data.K2O(k_glass_idx); median_k2o median(k2o_values, ‘omitnan’); data.K2O(isnan(data.K2O) k_glass_idx) median_k2o;2.2 成分数据的特殊性定和约束玻璃化学成分数据有一个根本特性——定和约束。即所有氧化物的百分比含量加起来应该等于100%或接近100%因检测误差。这个约束带来了两个关键问题冗余性知道了其中n-1种成分的含量第n种就能推算出来。这意味着我们的变量特征之间存在完美的线性关系会直接导致协方差矩阵奇异让许多基于距离或协方差的模型如PCA、LDA出问题。闭合效应成分数据位于一个单纯形空间Simplex其几何结构与欧氏空间不同。任意两种成分的变化都不是独立的一个成分的增加必然导致其他一个或多个成分的减少。为了解决这个问题我们必须进行数据变换。最经典、最适合后续统计建模的方法是中心对数比变换。原理CLR变换能消除定和约束将数据映射到欧氏空间。其公式为对于样本i的某个成分x_i变换后为ln(x_i / g(x))其中g(x)是该样本所有成分的几何平均数。Matlab实现function data_clr clr_transform(data) % data: n_samples x n_features 矩阵 已处理缺失值无零值 geo_mean exp(mean(log(data), 2)); % 按行求几何平均 data_clr log(data ./ geo_mean); % 中心对数比变换 end注意CLR变换要求数据中不能有零值因为要取对数。这就是为什么之前我们把微量成分的缺失值填为0会带来麻烦。更严谨的做法是对于这些“未检出”的零值在变换前用一个极小的值如检测限的一半进行替换这个技巧在成分数据分析中称为“零值处理”。2.3 特征选择与构建原始特征就是各种氧化物。但我们可以根据化学知识构建更有意义的衍生特征这往往能显著提升模型效果。比值特征PbO/BaO铅钡比是区分铅钡玻璃亚类的关键指标。K2O/Na2O钾钠比可能反映高钾玻璃的不同原料来源。SiO2/(K2ONa2OCaOPbOBaO)玻璃形成体与助熔剂之比可能和玻璃的稳定性、耐风化性相关。类别特征编码如果后续模型需要将“类型”高钾/铅钡、“风化情况”风化/未风化等文本标签转化为数值标签如1, -1或0, 1。经过这一系列预处理我们得到的data_clr矩阵才是真正适合进行后续统计分析和机器学习建模的“干净”数据。这个过程就像炼丹前的药材处理虽繁琐但决定了丹药的成败。3. 风化与类型鉴别从有监督模式识别到模型对比预处理后的数据我们首先解决两个有相对明确目标的问题判断一个样本点是否风化以及鉴别玻璃的类型是高钾还是铅钡。题目中部分样本给出了这些标签这为我们提供了训练有监督模型的可能。3.1 问题拆解与思路选择这里实际上包含两个子任务风化情况判别对于表单1有类型和风化标签这是一个二分类问题。我们可以用这部分数据训练模型然后去预测表单2无标签数据的风化情况。但要注意表单1的数据是“表面成分”其风化特征可能已经包含在成分模式中。玻璃类型鉴别对于表单2中未风化点的数据我们需要判断其类型。一个直接的思路是用表单1中未风化样本的数据已知类型作为训练集构建分类器对表单2的未风化点进行分类。3.2 特征空间的可视化初探在动手建模前一定要先“看”数据。主成分分析PCA是降维可视化的利器。我们对CLR变换后的数据做PCA并取前两个或三个主成分画散点图。[coeff, score, latent, ~, explained] pca(data_clr); figure; gscatter(score(:,1), score(:,2), glass_type); % glass_type是类型标签 xlabel([‘PC1 (‘, num2str(explained(1)), ‘%)’]); ylabel([‘PC2 (‘, num2str(explained(2)), ‘%)’]); legend(‘Location’, ‘best’); title(‘PCA Plot of Glass Samples (by Type)’);通过这个图我们可以直观地看到“高钾”和“铅钡”两类样本在成分空间上是否已经自然分开。如果分开得比较好说明成分本身对类型有很强的区分度后续分类任务会比较容易。我们还可以用不同标记形状表示“风化”与“未风化”观察风化作用是否会使样本点在PCA空间中发生系统性偏移。3.3 分类模型的选择、实现与对比仅仅看图不够我们需要定量的分类模型。这里的关键不是追求最复杂的模型而是理解不同模型的特性并选择适合数据特性的那一个。线性判别分析LDA为什么用LDA假设不同类别数据服从同协方差的正态分布其目标是最大化类间散度与类内散度之比。对于成分数据经过CLR变换后近似于正态分布且两类玻璃的化学成分差异可能体现在线性组合上因此LDA是一个强有力的基准模型。Matlab实现% 假设 train_data, train_label 是训练集 MdlLinear fitcdiscr(train_data, train_label, ‘DiscrimType’, ‘linear’); pred_label predict(MdlLinear, test_data); accuracy sum(pred_label test_label) / numel(test_label);注意事项LDA对特征间的多重共线性敏感。虽然CLR变换缓解了定和约束但某些氧化物之间仍可能存在较强相关性。可以先用PCA降维再用LDA即PCA-LDA或者使用正则化的判别分析‘DiscrimType’, ‘diagLinear’。支持向量机SVM为什么用如果PCA图显示两类样本的边界是非线性的SVM特别是带高斯核的SVM就能大显身手。它能通过核函数将数据映射到高维空间找到一个最优的非线性分割超平面。Matlab实现与调参% 使用高斯核SVM并优化关键参数 rng(1); % 重现性 Mdl fitcsvm(train_data, train_label, ‘KernelFunction’, ‘rbf’, … ‘Standardize’, true, … ‘OptimizeHyperparameters’, ‘auto’); % 或者手动尝试不同参数组合 box_constraint [0.1, 1, 10]; kernel_scale [0.1, 1, 10]; for b box_constraint for k kernel_scale Mdl fitcsvm(train_data, train_label, ‘KernelFunction’, ‘rbf’, … ‘BoxConstraint’, b, ‘KernelScale’, k); % … 交叉验证计算准确率 … end end核心技巧SVM一定要进行标准化‘Standardize’, true因为各氧化物含量量纲一致但数值范围不同。核函数的选择和参数如BoxConstraint和KernelScale对结果影响巨大必须通过交叉验证网格搜索来确定。集成学习模型如随机森林为什么用随机森林能自动处理非线性关系对异常值和过拟合相对鲁棒还能给出特征重要性排序。这对于我们理解“究竟是哪种化学成分对区分类型/风化起决定性作用”非常有帮助。Matlab实现与特征重要性分析Mdl fitcensemble(train_data, train_label, ‘Method’, ‘Bag’, … ‘NumLearningCycles’, 200, ‘Learners’, ‘tree’); % 计算特征重要性基于OOB误差 imp oobPermutedPredictorImportance(Mdl); % 可视化 figure; bar(imp); xlabel(‘Predictor’); ylabel(‘Importance’); title(‘Feature Importance by Random Forest’); set(gca, ‘XTickLabel’, feature_names); % feature_names是特征名称3.4 模型评估与结果融合我们不能只用一个模型。标准的做法是采用k折交叉验证例如5折或10折来评估每个模型的泛化性能避免因数据划分偶然性导致的过拟合评价。cv cvpartition(train_label, ‘KFold’, 5); accuracy_lda crossval(‘mcr’, train_data, train_label, ‘Predfun’, lda_predict, ‘partition’, cv); % 类似地计算SVM和RF的交叉验证准确率比较LDA、SVM、随机森林在交叉验证下的平均准确率。通常我们会选择性能最优的模型作为最终预测模型。但还有一种更稳健的策略是模型集成例如对几个表现相近的模型的预测结果进行投票硬投票或平均概率软投票这往往能获得比单一模型更稳定、更准确的结果。实操心得在这个问题上我们当时发现一个有趣的现象。对于“类型鉴别”LDA和线性SVM的效果就已经非常好交叉验证准确率95%这说明“高钾”和“铅钡”在化学成分上存在近乎线性的可分边界。而“风化判别”则稍难一些非线性模型如RBF-SVM表现更优这可能意味着风化过程对成分的影响模式更为复杂。特征重要性分析显示对于区分类型PbO、BaO、K2O、SiO2的权重最高这与化学常识完全吻合。4. 亚类划分无监督聚类与古代“配方”的揭秘这是本题最精彩也最具挑战性的部分。我们需要在没有先验标签的情况下对“高钾玻璃”和“铅钡玻璃”分别进行亚类划分。这完全依赖于数据内在的结构属于典型的聚类分析。4.1 聚类前的关键准备特征子空间选择不能直接把所有CLR变换后的成分扔给聚类算法。因为噪声干扰像P2O5、SrO这些含量极低且波动大的微量元素可能对聚类产生噪声。相关性高度相关的特征如PbO和BaO在铅钡玻璃中可能此消彼长会给基于距离的算法带来误导。我们的策略是基于化学知识筛选聚焦主要成分和关键比值。对于铅钡玻璃核心特征是PbO、BaO及其比值PbO/BaO可能还有SiO2、CuO着色剂等。对于高钾玻璃核心是K2O、SiO2、CaO、MgO等。基于有监督分析的结果利用上一节随机森林得出的特征重要性排名选取对区分大类贡献大的特征用于该大类的亚类划分。降维对筛选后的特征子集进行PCA并选择累积贡献率超过85%的前几个主成分作为聚类输入。这既能去除噪声和相关性又能降低维度使聚类更稳定。4.2 聚类算法的选择与实战我们对比了两种最主流的聚类算法K-means聚类原理与适用性试图将样本划分成K个簇使得每个样本到其所属簇中心的平方欧氏距离之和最小。它适用于球形簇、簇大小相近的情况。对于古代玻璃配方如果不同亚类的配方在成分空间形成相对紧凑、分离的“球团”K-means会是一个好选择。Matlab实现与确定K值% 假设 data_for_cluster 是经过筛选和PCA降维后的高钾玻璃数据 max_k 6; % 假设我们最多尝试划分6个亚类 eval evalclusters(data_for_cluster, ‘kmeans’, ‘CalinskiHarabasz’, ‘KList’, 2:max_k); figure; plot(eval); xlabel(‘Number of Clusters’); ylabel(‘Criterion Value (CalinskiHarabasz)’); title(‘Optimal K Evaluation’); optimal_k eval.OptimalK; [idx_kmeans, C] kmeans(data_for_cluster, optimal_k, ‘Replicates’, 10, ‘Display’, ‘final’);关键参数‘Replicates’重复次数非常重要K-means对初始簇中心敏感多次运行取最优结果可以避免局部最优。‘CalinskiHarabasz’方差比准则是常用的内部评估指标值越大通常意味着聚类结构越好。层次聚类Hierarchical Clustering原理与适用性通过计算样本间的距离逐步合并聚合式或分割分裂式来构建一个树状图系统树图。它不需要预先指定K值并且可以通过树状图直观地看到不同层次的分组情况特别适合探索性分析以及当亚类结构可能存在嵌套关系时比如一个大亚类下还可细分。Matlab实现与树状图解读% 计算距离矩阵例如欧氏距离 Y pdist(data_for_cluster); % 创建系统树图使用沃德法连接其对球形簇效果较好 Z linkage(Y, ‘ward’); figure; dendrogram(Z); title(‘Dendrogram for Hierarchical Clustering’); xlabel(‘Sample Index’); ylabel(‘Distance’); % 根据树状图在某个距离阈值上切割得到聚类标签 T cluster(Z, ‘MaxClust’, optimal_k); % 或使用 ‘Cutoff’ 参数如何看树状图纵轴表示合并时的距离。寻找纵轴上距离较长的水平线段然后画一条横线穿过它这条线与垂直线的交点数量就大致是合理的聚类数。水平线段越长说明合并的两个簇差异越大划分越清晰。4.3 聚类结果的验证、解读与化学溯源得到聚类标签后工作只完成了一半。更重要的是解读聚类结果并赋予其化学或考古学意义。可视化验证使用t-SNE一种更擅长保持局部结构的非线性降维方法将数据降至2维或3维并用聚类标签着色直观检查聚类是否分离良好有无明显重叠。Y_tsne tsne(data_for_cluster, ‘NumDimensions’, 2, ‘Perplexity’, 30); gscatter(Y_tsne(:,1), Y_tsne(:,2), idx_kmeans); title(‘t-SNE visualization of Clustering Result’);剖面分析计算每个簇亚类在所有原始化学成分上的均值或中位数绘制雷达图或平行坐标图。这是最关键的一步。% 假设 original_data 是原始成分数据 idx 是聚类标签 clusters unique(idx); figure; for i 1:length(clusters) cluster_data original_data(idx clusters(i), :); mean_composition mean(cluster_data, 1); % 绘制雷达图 (需要自定义或使用其他工具箱函数如polarplot) % 这里以条形图示例对比 subplot(2, ceil(length(clusters)/2), i); bar(mean_composition); set(gca, ‘XTickLabel’, chemical_names, ‘XTickLabelRotation’, 45); title([‘Cluster ‘, num2str(clusters(i)), ‘ Mean Composition’]); ylabel(‘Content (%)’); end通过对比这些“成分剖面图”我们可以描述每个亚类的特征铅钡玻璃亚类1可能是“高铅低钡”型PbO/BaO比值很高可能SiO2含量较低颜色可能因高铅而呈黄色或褐色。铅钡玻璃亚类2可能是“铅钡均衡”型PbO和BaO含量相当可能含有特定比例的CuO绿色或Fe2O3褐色。高钾玻璃亚类1可能是“高钾钙”型K2O和CaO含量都高稳定性较好。高钾玻璃亚类2可能是“高钾镁”型K2O和MgO含量高可能使用了不同的植物灰原料。统计检验使用ANOVA方差分析检验各个化学成分在不同亚类间的均值是否存在显著差异从统计上确认划分的有效性。踩坑实录我们最初对所有特征进行K-means聚类结果非常混乱t-SNE图上一片模糊。后来意识到是SrO、Cl等微量元素噪声太大。当我们只聚焦PbO、BaO、SiO2、CuO、K2O等主要和特征性成分后聚类结果立刻变得清晰树状图上出现了明显的长距离分支。雷达图分析显示我们为铅钡玻璃划分的三个亚类其PbO/BaO比值形成了低、中、高三个清晰的梯度这很可能对应了三种不同的配方传统或工艺阶段。5. 关联分析与风化规律探索数据背后的化学故事完成分类和聚类后我们还可以深入挖掘数据中更深层的关系回答诸如“风化究竟改变了什么”、“不同类型玻璃的风化行为有何不同”等问题。这需要用到相关性分析、统计检验和回归建模。5.1 风化前后的成分对比分析对于表单1中同时有“风化”和“未风化”检测点的样本我们可以进行配对样本分析。这是揭示风化规律最直接的证据。% 假设 weathered_data 和 unweathered_data 是配对样本的风化与未风化点数据 % 计算每个成分的风化前后差值 diff_data weathered_data - unweathered_data; % 计算各成分差值的均值及置信区间 mean_diff mean(diff_data, 1, ‘omitnan’); [~, ~, CI] ttest(diff_data); % 对每列每个成分进行配对t检验 % 可视化绘制带有置信区间的均值差值条形图 figure; bar(1:size(diff_data,2), mean_diff); hold on; errorbar(1:size(diff_data,2), mean_diff, mean_diff - CI(1,:), CI(2,:) - mean_diff, ‘k.’, ‘LineWidth’, 1.5); set(gca, ‘XTick’, 1:size(diff_data,2), ‘XTickLabel’, chemical_names, ‘XTickLabelRotation’, 45); ylabel(‘Mean Change after Weathering (%)’); title(‘Compositional Changes Due to Weathering’); hline refline(0,0); hline.Color ‘r’; hline.LineStyle ‘—’;通过这个分析我们可以清晰地看到哪些成分显著流失了如Na2O,K2O它们的差值均值显著为负且置信区间不包含0。这些是易溶的碱金属氧化物。哪些成分相对富集了如SiO2,Al2O3差值均值为正。因为它们不易溶在碱金属流失后其相对百分比升高。哪些成分基本不变如PbO,BaO置信区间包含0。这表明铅钡玻璃中的铅、钡元素相对稳定。5.2 风化程度与成分的关联模型我们还可以尝试量化“风化程度”。一个简单的代理指标是风化指数例如定义为(SiO2 Al2O3) / (Na2O K2O CaO)在风化点与未风化点的比值。比值越大风化越严重。 接下来可以探索哪些初始成分未风化点的成分能够预测该文物的风化指数。这可以通过多元线性回归或回归树来实现。% 假设 weathering_index 是风化指数 initial_composition 是未风化点成分 % 使用逐步回归自动选择重要预测因子 mdl stepwiselm(initial_composition, weathering_index, ‘Upper’, ‘linear’, ‘Criterion’, ‘aic’); disp(mdl); % 查看回归方程和显著性 % 或者使用回归树可视化决策路径 tree fitrtree(initial_composition, weathering_index, ‘PredictorNames’, chemical_names); view(tree, ‘Mode’, ‘graph’);如果回归模型显示高含量的K2O或特定的PbO/BaO比值与更高的风化指数相关那么我们就能得出“某种配方可能更易风化”的推论这对文物保护具有指导意义。5.3 亚类与风化、纹饰的关联性分析最后我们可以做一个综合性的探索将我们聚类得到的亚类标签与“是否风化”、“纹饰类型”等属性进行列联表分析或卡方检验看看是否存在统计上的关联。% 创建亚类与风化情况的列联表 contingency_tbl crosstab(cluster_labels, weathering_labels); % 进行卡方检验 [h, p, stats] chi2gof(…); % 需要适当的数据准备如果发现某个亚类的玻璃风化比例显著高于其他亚类或者某个亚类与特定的纹饰如“弦纹”、“素面”强相关这就能将化学成分分析、工艺分类与考古学的类型学、艺术史研究联系起来形成更完整的证据链。个人体会做到这一步数学建模才真正与文物研究产生了深度共鸣。我们不再仅仅是报告“聚类为3类”而是可以提出“根据成分分析这批铅钡玻璃至少存在三种配方传统。其中高铅型的配方亚类A在埋藏环境中表现出更强的抗风化特性而含有特定比例铜元素的亚类B则与‘蜻蜓眼’纹饰有较高的关联性这可能指示了其特定的文化来源或用途。” 这样的结论才是赛题所期望的、有洞察力的分析。整个过程中Matlab从数据清洗、统计分析到可视化提供了一站式的解决方案其交互式环境和丰富的图形功能让我们在迭代分析和呈现结果时效率倍增。