简介面向参加全国大学生数学建模竞赛的学生与指导教师聚焦2022年C题古代玻璃文物成分分析与鉴别可作为赛题复盘、建模方法学习和论文写作参考。包体共1个PDF文件约3.55MB内容围绕完整赛题求解展开涵盖数据预处理与探索、风化系数定义、主成分分析降维、因子分析预测未风化成分、亚类划分、敏感性分析及灰色关联度分析等关键环节并给出问题重述、模型构建、结果解释与结论表述。内容预览显示对高钾玻璃与铅钡玻璃的风化差异、分类阈值、未知类别鉴别和成分关联度排序均有较细说明。已有928人学习适合需要梳理C题思路、理解统计建模流程或准备同类考古与文物保护数据题目的读者参考也可借鉴论文结构、公式表达与结果分析写法。1. 风化系数一块玻璃「风化没风化」为什么不能用眼睛判断拿到这道数学建模赛题的第一眼很多人会把「风化」当成一个现成的分类标签——表单里不是明明白白标了风化、无风化吗真正动手算才发现标签和成分对不上有的样品标着无风化氧化铅占比却高得离谱。风化其实是成分比例的定向漂移并非外观特征。这道题要建的模型是把「漂移了多少」压成一个可比较的数再反过来把漂移推回去。搞懂这套流程你拿到的就不只是一份答案而是一套处理「观测被环境改写过」的通用思路——从考古成分数据到设备老化后的检测值内核是同一个。2. 数据探索纹饰、类型、颜色三张牌怎么把风化规律锁死在建模之前这道题最值钱的部分其实发生在 Excel 和几张交叉表里。先把「谁更容易风化」这件事从数据里读出来后面的因子分析、区间划分才有落脚点。跳过这一步直接上模型等于闭着眼睛调参。2.1 交叉表先于模型三个观察事实怎么被读出来表单 1 只有类型、纹饰、颜色、是否风化四个字段看起来信息量很小但按纹饰分组切开之后规律强得惊人。纹饰 B 的样品全部是高钾玻璃、全部风化、颜色全部是蓝绿色纹饰 A 的高钾玻璃全部无风化铅钡玻璃则是风化比无风化约 11:5纹饰 C 的高钾玻璃同样全部无风化铅钡玻璃约 17:7。这三条事实指向两个结论同等条件下铅钡玻璃比高钾玻璃更容易被风化纹饰 B 所对应的埋藏环境风化效应更强强到连相对耐风化、以二氧化硅为主的高钾玻璃都全军覆没。值得单独说一句的是卡方检验的结果。纹饰与表面风化的 χ² 约 16.4、p 值在千分之一量级颜色与风化的 χ² 约 26.0、p 值约 0.001效应量 Phi 分别到 0.561 和 0.707属于强关联。而「类型」与「表面风化」的相关系数只有 0.272p 值 0.051刚好卡在 5% 的边缘。这个「不显著」恰恰是重要信息玻璃类型本身并不决定风化与否真正决定风化强度的是材质与环境两个因素的叠加。后面所有分组建模本质上都是在给这个判断补细节。2.2 风化系数把「漂移方向相反」写成一个比值数据分析里最容易被忽略的一步是把化学常识翻译成可计算的指标。两类玻璃的风化方向完全相反铅钡玻璃风化后二氧化硅大幅流失氧化铅、氧化钡相对富集高钾玻璃风化后二氧化硅占比显著上升氧化钾、氧化钙、氧化镁、氧化铝等助熔剂成分显著下降。如果只用某一种成分的变化量做判别两类玻璃就会打架。于是有了风化系数这个构造分子取「风化后富集的那一类成分之和」分母取「风化后流失的那一类成分之和」。铅钡玻璃的分子是 PbO、BaO分母是 SiO₂高钾玻璃反过来分子是 SiO₂分母是 K₂O、CaO、MgO、Al₂O₃、Fe₂O₃、CuO。用比值而非差值是这个指标最核心的设计动机——成分数据是百分比各项之和恒等于 100一个成分涨必然有别的成分跌差值会互相抵消比值才能把「相对富集」表达清楚而且天然无量纲跨样品可比。import numpy as np # 风化系数分子取风化后富集的成分分母取风化后流失的成分 # 两类玻璃的风化方向相反所以分子分母要对调不能共用一个公式 LEAD_BA_RICH [PbO, BaO] # 铅钡玻璃风化后富集项 POTASH_RICH [SiO2] # 高钾玻璃风化后富集项 LEAD_BA_LOSS [SiO2] # 铅钡玻璃风化后流失项 POTASH_LOSS [K2O, CaO, MgO, Al2O3, Fe2O3, CuO] # 高钾玻璃风化后流失项 def weathering_coeff(row, kind, eps1e-6): row: dict键为成分名值为百分比kind: 铅钡 或 高钾 if kind 铅钡: num sum(row.get(c, 0.0) for c in LEAD_BA_RICH) den sum(row.get(c, 0.0) for c in LEAD_BA_LOSS) eps # 防分母为 0 else: num sum(row.get(c, 0.0) for c in POTASH_RICH) den sum(row.get(c, 0.0) for c in POTASH_LOSS) eps return num / den这里要提醒一句原始文档里给出的系数在组内还做过一次尺度缩放所以你不能拿表格里的百分比手算复现出 0.818、3.2369 这种数。能复现的是同一组样品内部的排序以及下面要讲的跳跃带位置——那才是这个指标真正承载信息的地方。2.3 阈值不是拍的从 0.879 到 1.3 的跳跃带系数算出来之后紧接着的问题是怎么切分风化与未风化。常见的坑是直接拿均值加减若干倍标准差或者上 KMeans 聚两类。前者会被极端风化样本拖偏后者强行把连续过程切成两半切点位置全看初始化。这道题用的是最朴素也最稳的办法把同一组样品的系数排序看相邻两值的间隔在哪里突然放大。铅钡玻璃 A 纹的系数在 0.879 与 1.3 之间出现明显断层C 纹在 0.8 与 1.2 之间同样断层于是铅钡组的阈值统一放在 1.0高钾玻璃在 7.23 与 12.41 之间跳跃阈值放在 10.0。断层之所以存在是因为风化本身是个阈值型过程——所有玻璃都在缓慢风化只有超过某个程度才被称为风化玻璃过渡区样本天然稀少。def find_threshold(coefs, round_toNone): 在排序后的系数序列里找最大跳跃带返回带中点作为候选阈值 c np.sort(np.asarray(coefs, dtypefloat)) gaps np.diff(c) i int(np.argmax(gaps)) # 断层位置 mid 0.5 * (c[i] c[i 1]) # 断层带中点 return (round(mid) if round_to else mid), gaps[i] # 铅钡组断层中点约 1.09取整到 1.0 # 高钾组断层中点约 9.82取整到 10.0取整到 1.0 和 10.0 是刻意的阈值落在断层带的任何位置判别结果都一样取整只是让规则好记、好写进文档。还有一个细节值得学——原始文档里有三个样本编号 20、48、30 部位 1的人工标签与系数判别结果相反作者选择以系数为准手工修正标签。这个决定的依据是标签来自现场目视或记录系数来自成分实测当两者冲突时物理测量优先。3. 因子分析对负荷矩阵求逆为什么能把风化后的成分推回原配方问题一真正的难点在最后一问给定一个已经风化的样品预测它风化之前的成分。这是典型的反问题——你观测到的是被环境改造过的结果要倒推原始输入。因子分析在这里承担的角色是把「成分向量」拆成「公共因子 特殊误差」再对这个映射求逆。3.1 PCA 在这条路线里只做一件事压掉共线性十三四个成分指标之间高度相关二氧化硅涨、助熔剂就跌直接建模会撞上多重共线性矩阵接近奇异求逆时数值误差被放大到不可用。PCA 的作用是把这十几个相关变量旋转成若干个互不相关的公共因子同时把维度从十几维降到八维左右累计贡献率还能保住 97% 以上。顺序很重要先 PCA 拿主成分再用主成分构造因子负荷矩阵最后才谈反推。跳过降维直接对原始成分矩阵做因子分析特征值谱尾部全是噪声因子数的选择会变得非常主观。这里的选择标准很直接——保留累计方差贡献率 85% 以上问题一对精度要求更高取到 97%对应的特征值个数就是因子数 k。降维带来的另一个好处是抑制干扰末尾那些只解释 1% 方差的主成分装的几乎全是采样误差和测量噪声把它们留在反推链路里只会让预测值抖得没法看。3.2 从 x μ A f ε 到 f 的线性表达因子模型写成矩阵形式是 x μ A f ε其中 x 是成分特征向量μ 是均值向量f 是公共因子这里对应未风化时的成分结构ε 是特殊因子即随机误差A 是负荷矩阵元素 a_ij 表示第 i 个成分在第 j 个公共因子上的投影。关键在于 A 怎么来。对相关系数矩阵 R 做特征分解得到降序排列的特征值 λ₁ ≥ λ₂ ≥ … 和对应的标准正交特征向量 u₁, u₂, …负荷矩阵取 A [√λ₁u₁, √λ₂u₂, …]。这一步有个容易被抄错的细节原始文档里写的是 A [λ₁u₁, λ₂u₂, …]但落到代码里必须开根号否则重构出来的协方差矩阵 AAᵀ Ψ 与 R 对不上后面的反推系数全部跑偏。这是因子分析实现里最常见的翻车点之一。import numpy as np def factor_model(X, var_ratio0.97): X: (n_samples, n_features)返回均值、负荷矩阵、因子得分系数矩阵 X np.asarray(X, dtypefloat) mu X.mean(axis0) R np.corrcoef(X, rowvarFalse) # 相关系数矩阵先标准化再算 vals, vecs np.linalg.eigh(R) # 对称矩阵用 eigh快且稳 order np.argsort(vals)[::-1] vals, vecs vals[order], vecs[:, order] k int(np.searchsorted(np.cumsum(vals) / vals.sum(), var_ratio) 1) A vecs[:, :k] * np.sqrt(vals[:k]) # 开根号别写成 vals[:k] B np.linalg.pinv(R) A # 因子得分系数矩阵形状 (m, k) return mu, A, B, k def topc_from_observed(x, mu, B): 由风化后的成分反推公共因子即未风化成分的低维表示 return B.T (x - mu)B 就是原始文档里说的「转换系数矩阵」。注意它并不是把 A 直接求逆得到的——A 是 m×k 的长方阵m 个成分、k 个因子当 k m 时根本不存在逆矩阵np.linalg.inv 会直接抛 Last 2 dimensions of the array must be square。合法写法有两条路一是回归法 B R⁻¹A即上面代码用的伪逆形式二是主成分法 B A 本身当特殊因子方差 Ψ 取 0 时二者等价。原始文档说的「对 A 求逆」实际落地时走的是前者这个差别不搞清楚代码就跑不起来。3.3 预测表怎么读、误差从哪来拿到 B 之后铅钡玻璃 C 纹那一组用 8 个主成分累计贡献率 97% 以上代入得到一张从风化成分到未风化成分的转换系数表逐行乘加就得到每个采样点的未风化成分预测值。表里能看到几个有意思的现象严重风化点的二氧化硅预测值从观测的二十几回补到五六十氧化铅从四十几回补到二十上下方向与化学机理一致同一采样点的「部位 1」「部位 2」预测值高度接近说明反推稳定。但误差必须交代清楚它来自四个地方。一是假设本身模型假定风化过程中成分变化基本线性、同类型玻璃的风化反应相同一旦某个样品风化到非线性区间反推必然偏。二是因子数截断扔掉的那部分主成分虽然方差小但反推时它们的方向可能被 B 放大。三是逆映射的条件数R 接近奇异时微小的观测误差经过 R⁻¹ 会被指数级放大这也是必须先做 PCA 的核心理由。四是采样点本身的空间异质性——同一件文物的不同部位风化程度不同把它们当同分布样本处理本身就是近似。实际使用时最稳妥的做法是给出预测区间而非点估计并在报告里明确标注哪些样本落在因子空间的边缘。4. 区间分类没有标签的亚类怎么用「掐头去尾」划出来问题二要解决的是无监督划分类别。给了一堆高钾、铅钡玻璃的成分数据要求在每个大类内部再分亚类还得说明划分方法合理、对扰动不敏感。既有聚类算法可以选也可以像这道题一样走「主成分筛特征 极值区间 一维聚类切分」的路线。4.1 贡献率筛选88.89% 和 95.41% 分别意味着什么先把数据按「高钾风化、高钾未风化、铅钡风化、铅钡未风化」四组切开每组单独做主成分分析取贡献率最高的几个成分作为判别特征。铅钡未风化组选出 SiO₂、Na₂O、K₂O、CaO、MgO、Al₂O₃ 六个成分累计贡献率 88.89%铅钡风化组选五个成分87.36%高钾风化组只用 SiO₂、Na₂O、K₂O 三个成分就达到 95.41%。高钾风化组维度这么低是因为严重风化之后其他成分被稀释到接近零方差几乎全集中在二氧化硅上。原始文档里高钾未风化组的贡献率也写成 88.89%与前文铅钡组完全一致这处大概率是复用了上一段的数字。复现时以实际算出的累计贡献率为准不要照抄。这个细节看着小但它提醒你一件事文档里的数字要交叉验算尤其是同一段落里反复出现的同一个百分比。4.2 去掉最大值最小值之后的区间为什么能当判别器特征选定后对每个特征做一件很土但很有效的事排序去掉一个最大值和一个最小值用剩下的最小值和最大值构成区间。四个组、几个关键成分最后得到一张区间表比如铅钡未风化组的 SiO₂ 区间是 [50.61, 69.71]高钾未风化组是 [59.81, 69.33]铅钡风化组是 [4.61, 37.36]高钾风化组是 [92.35, 95.02]。判别规则就一句话待测样品的成分落在哪个组的区间里或足够靠近就归到哪一组。掐头去尾的动机值得展开。区间法的分辨力完全取决于区间宽度而极值点往往是两类样本——极端风化的样品或者采样、测量环节出了偏差的样品。把它们留在里面区间会被撑宽一大截本来能分开的两组就糊在一起了。去掉两端各一个点是一种极低成本的稳健化处理比用什么 3σ 准则更保守因为它完全不假设分布形态。这也解释了为什么区间法的适用前提是「组内样本不能太少」样本只有七八个的时候去掉两个点剩下的范围已经很窄判别会过于激进。4.3 最短距离法聚点把一维成分切成两段亚类划分选的是数据差异最大的成分——二氧化硅。四个组各自用最短距离法聚成两类得到两个聚点铅钡未风化 64.82 和 52.49高钾未风化 71.22 和 62.77铅钡风化 33.96 和 15.63高钾风化 94.66 和 92.57。切分阈值取两聚点的算术平均于是得到 58.70、67.00、24.80、93.62 这四个分界值——把它们和原始文档里的子区间端点对照完全吻合。这条对应关系是验证整条链路有没有理解对的最好证据。import numpy as np def single_linkage_1d(x): 一维最短距离法聚成 2 类等价于在最大间隔处切开 x np.sort(np.asarray(x, dtypefloat)) gaps np.diff(x) i int(np.argmax(gaps)) return x[:i 1], x[i 1:] def centroids_and_split(x): left, right single_linkage_1d(x) c1, c2 left.mean(), right.mean() # 两个簇的聚点 return sorted([c1, c2]), 0.5 * (c1 c2) # 切分阈值取聚点平均 # 铅钡未风化[52.49, 64.82] - 阈值 58.66取两位小数 58.70 # 高钾未风化[62.77, 71.22] - 阈值 67.00 # 铅钡风化 [15.63, 33.96] - 阈值 24.80 # 高钾风化 [92.57, 94.66] - 阈值 93.62用聚点平均而不是直接用断层位置当阈值是有意为之断层位置对单侧极端值很敏感一个离群点就能把切点拽走聚点平均是两个簇中心的中点受簇内分布影响抗噪能力更强。代价是当两个簇大小严重不平衡时中点会偏向大簇这一点在样本量差异大的组里要额外检查。4.4 3% 扰动问题二敏感、问题三不敏感矛盾吗敏感性分析的做法很统一给原始数据乘上 (1 ± 3%) 的随机扰动重新跑一遍分类看结果变了多少。问题二的结果是四类中有三类的分布比例误差超过 20%判定为敏感性强问题三把表单 3 的数据同样扰动 3% 再判别类别完全不变判定为敏感性不强。乍看矛盾其实两件事衡量的根本不是同一个东西。问题二统计的是「样本落在第一子区间还是第二子区间」的比例变化而切点附近本来就挤着一批成分接近中点的样本3% 的扰动足以让它们跨过阈值比例误差自然大——这是区间法固有的边界效应不代表分类逻辑错了。问题三判别的是四个大类的主区间组与组之间的间隔极宽比如高钾风化组 SiO₂ 在 92 以上铅钡风化组在 4 到 37 之间3% 的扰动根本跨不过去。这个对比给出的实践启示是敏感性分析必须说清楚「扰动作用在哪个判别层级上」。粗粒度判别稳、细粒度切分抖是正常现象。如果反过来——粗判别也大面积翻转那才是模型真出了问题。5. 避坑与排查这类成分数据建模最容易翻车的五处下面五条都是这套流程里真实会卡住人的地方按「现象 → 原因 → 解决」写遇到时可以直接照着排查。5.1 数据关联与阈值划分阶段的坑其一两张表按编号 merge 之后行数变少、出现大量 NaN。现象是合并后的表比左表短或者成分列整片缺失。原因是表单 1 的采样点编号带「部位 1」「部位 2」「未风化点」这类后缀而表单 2 的编号是另一套写法空格、全半角、括号都可能不一致字符串直接匹配必然漏。解决办法是把编号拆成「纯数字主体 部位标识」两列做复合键merge 时带上 indicatorTrue合并后先断言左表行数没变再检查 _merge 列里有没有 left_only。其二阈值硬编码成 1.0 用在所有组上。现象是铅钡 A 纹判得挺准换到 C 纹或高钾组就大面积误判。原因是各组系数的尺度不同同一套组内缩放在不同组之间不可比铅钡组阈值 1.0、高钾组阈值 10.0 差了整整一个数量级。解决办法是按「类型 纹饰」分组每组各自排序找断层带把阈值写进配置表而不是写死在函数里。其三把人工标签当成金标准。现象是模型判别结果与标签冲突时第一反应是回去调阈值迁就标签。原因是标签来自目视或现场记录本身就有误差而系数来自成分实测。正确做法是冲突样本逐个复核以成分事实为准修正标签并且在报告里列出所有被修正的样本编号和修正理由——这既是可复现性要求也是审稿人最关心的地方。5.2 建模与验证阶段的坑其四直接对负荷矩阵调 np.linalg.inv。现象是抛出 Last 2 dimensions of the array must be square。原因是 A 的形状是 m×k成分数大于因子数压根不是方阵。解决办法是改用回归法 B pinv(R) A或者用主成分法令 B A同时在求逆前后打印 R 的条件数条件数超过 1e4 就说明共线性没压干净得回头加大降维力度。其五扰动实验不固定随机种子也不做重复。现象是同一份代码跑两次误差率从 18% 跳到 26%结论没法写进文档。原因是一次性随机扰动的方差本身就大。解决办法是固定 np.random.default_rng(seed)扰动用乘性因子 (1 U(-0.03, 0.03)) 而不是加性并且重复跑 100 次报告误差率的均值与 95% 分位区间用区间说话而不是用单次结果下结论。再补一条容易被忽略的灰色关联度分析前不做无量纲化。现象是所有成分的关联度都挤在 0.9 以上排序几乎等于含量排序。原因是 SiO₂ 的数值在 60 上下Na₂O 只有 1 左右绝对差值被大数量级成分完全主导。解决办法是先做均值化 x_i(k) / mean(x_i)让各序列处在同一量级再计算关联系数——这一步不做后面的排序结论基本没有意义。6. 灰色关联度与设计取舍ρ 0.5 到底在权衡什么问题四要比较不同类别内部化学成分之间的关联强度样本量最小的组只有 7 条记录。这个规模下皮尔逊相关系数的置信区间宽到没法用而且成分数据是否近似正态本身就存疑。灰色关联度分析绕开了分布假设把某一成分序列当参考数列其余成分当比较数列做无量纲化后逐点计算绝对差用 ξ (Δmin ρΔmax) / (Δ ρΔmax) 得到关联系数再取均值作为关联度最后排序。分辨系数 ρ 的取值是整个方法里唯一需要拍板的地方它控制的是「最大差在分母里的权重」。ρ 越小Δmax 项占比越低关联系数之间的差距被拉大分辨力上升但噪声也被同步放大ρ 越大Δmax 主导分母所有关联系数被压向 1排序趋于平坦。工程上通常取 0.5此时分辨力与稳定性折中有经验说法认为 ρ ≤ 0.5463 时分辨力较好这条界线可以当作调参的起点但不必当成硬约束。对比项灰色关联度皮尔逊相关系数最小样本量4 条左右即可计算小样本下置信区间极宽分布假设不要求近似正态、线性关系对异常值中等敏感受 Δmax 牵制高度敏感一个点就能翻转输出含义序列形状的相似度排序方向 强度主要代价对无量纲化方式极敏感对非线性关系无能为力实测的排序结果里四个类别都是二氧化硅关联度最强这符合化学直觉——二氧化硅是骨架成分占比六成以上其余成分的含量变化本质上是在它的稀释效应下发生的。有意思的是风化前后的差异铅钡未风化的排序是 SiO₂ CaO K₂O Al₂O₃ MgO Na₂O高钾风化变成 SiO₂ Na₂O K₂O CaO MgO钠、钙的位置发生了对调。这说明风化不只是等比例稀释它改变了成分之间的耦合关系这一点用简单的相关系数矩阵是看不出来的。要提醒的是灰色关联度的排序结论不能跨数据集比较。换一种无量纲化方式初值化换成均值化、换一个 ρ排序尾部几个关联度接近的成分就可能互换位置。所以写报告时应当同时给出关联度数值和排序而不是只给排序对数值差距小于 0.02 的相邻项直接标注「差异不显著」比强行排名更诚实。这套流程走下来我最深的一条教训是凡是处理「观测值被环境改写过」的数据都别急着上模型先花时间把「改写方向」用一个比值量化出来再用分布里的断层确定切点。阈值定错了后面因子分析做得再漂亮反推出来的原始配方也是错的。从那以后我每次做类似的成分或指标反推题都强制先跑一遍断层搜索把切点位置和断层宽度打印出来看一眼再动模型。希望帮到你。本文还有配套的精品资源点击获取