数学建模C题玻璃文物风化反推：风化系数与PCA因子分析

简介面向参加全国大学生数学建模竞赛的学生与指导教师这份资料围绕2022年C题“古代玻璃文物的成分分析与鉴别”提供赛题解析、建模思路与求解参考适合赛前训练、复盘建模流程及查漏补缺。包内为1个PDF文件约3.55MB以文档形式集中呈现题目分析、数据处理与模型构建内容。已有928人学习说明其在建模备赛群体中具有一定参考价值。资料覆盖数据预处理与探索、风化系数定义、主成分分析降维、因子分析预测风化前成分、亚类划分、敏感性分析、灰色关联度分析以及模型评估等关键环节并给出问题一至问题四的求解脉络。读者可据此理解高钾玻璃与铅钡玻璃的风化差异、未知类别鉴别方法和分类结果稳健性检验也可借鉴其中的建模组织方式与结果解释思路。1. 这道C题真正拉开差距的是“反推”而不是“预测”数学建模竞赛C题给出一批古代玻璃文物的化学成分检测数据四个问题表面上是四道独立小题底下却共用一条主线先把“风化”这件事量化成一个可排序的指标再用降维和因子分析把风化后的成分反推回风化前最后做分类、做关联度。多数人栽在第一步——看到“预测风化前成分”就直接上回归翻遍数据才发现没有配对样本方向从一开始就错了。这场考试考的不是化学知识而是能不能把一个模糊的现实问题拆成可计算的步骤并在论文里把每一步的理由写到评委没法追问。四问里问题一、问题二占分最重问题三、问题四是拉开档次的区间。正在准备竞赛、准备数据分析岗面试或者要给统计建模类作业答辩的人都适合按这条线练一遍。2. 数据预处理与探索风化系数阈值是第一个硬得分点2.1 三张表单怎么对齐字段关联决定后面模型能不能跑原始数据分成三张表单第一张是文物级信息一件文物一行记着类型、纹饰、颜色、表面有无风化第二张是采样点级成分一件文物可能取好几个点每个点十几种氧化物含量第三张是待鉴别的未知样品。三张表能不能拼到一起全看采样点编号——编号里前几位是文物编号后面才是部位号需要先切出来再连接。这一步做错后面所有模型都是在错误的数据上跑写得再漂亮也拿不到分。成分合计列是自查工具。把所有氧化物百分比加起来接近100说明数据基本完整明显低于100说明有成分未检出或者录入缺失。未检出按0处理整行关键成分缺失才剔除并在论文里写一句处理规则。别小看这句话评委翻论文时最先看的就是数据清洗有没有交代清楚血泪经验是清洗步骤写不清楚的论文后面模型再花哨也会被怀疑结果来源。import pandas as pd import numpy as np # 表单1文物级信息表单2采样点级成分 df1 pd.read_excel(表单1.xlsx) df2 pd.read_excel(表单2.xlsx) # 采样点编号形如“01部位1”抽出前导数字作为文物编号并补齐两位 df2[文物编号] df2[文物采样点].astype(str).str.extract(r(\d))[0].str.zfill(2) # 按实际列名替换所有氧化物成分列 comp_cols [c for c in df2.columns if c not in [文物采样点, 文物编号]] df2[成分合计] df2[comp_cols].sum(axis1) # 左连接文物信息并到成分表右侧保证成分表的行数不变 df df2.merge(df1[[文物编号, 类型, 纹饰, 颜色, 是否风化]], on文物编号, howleft) # 未检出成分记0类型或风化标签缺失的行无法用于建模直接剔除 df[comp_cols] df[comp_cols].fillna(0) df df.dropna(subset[类型, 是否风化]) print(df[成分合计].describe()) # 合计明显偏离100的样品要单独核查代码里三个动作值得说清楚。str.extract(r(\d))抓的是编号开头的连续数字如果编号格式不统一这条正则就会抓错失败时先打印十几行原始编号肉眼核对。howleft保证以成分为主表不会因为某件文物信息缺失就把它的检测数据丢掉。fillna(0)的前提是“未检出即含量为0”如果原始表里用的是空白表示未检测就不能一概填0要回去看数据字典。这些细节在论文的“数据预处理”一节各写一两句比堆公式管用。2.2 风化系数阈值要找“断崖”不是取平均风化程度需要量化否则“有无风化”只是一列标签没法参与计算。思路是构造一个比值铅钡玻璃风化后氧化铅、氧化钡占比明显升高、二氧化硅下降就把铅钡氧化物之和放在分子、二氧化硅放在分母高钾玻璃正好相反风化后二氧化硅占比冲到九成以上、助熔剂氧化物大幅下降就把二氧化硅放分子、助熔剂氧化物之和放分母。两类玻璃的化学机理相反系数定义也必须分开混用一套公式是这道题最常见的翻车点之一。算完系数不要急着取0.5或者平均值当阈值。正确做法是按“类型纹饰”分组把系数从小到大排序看相邻两个值之间有没有明显的断崖。铅钡玻璃A纹的数据在0.879和1.3之间突然跳开阈值定1.0高钾玻璃在7.23和12.41之间跳开阈值定10.0。这样定出来的阈值有数据支撑论文里可以直接写“依据排序后的断崖区间确定”比拍脑袋有说服力。# 铅钡玻璃铅钡氧化物之和 / 二氧化硅 mask_pb df[类型] 铅钡 df.loc[mask_pb, 风化系数] ( df.loc[mask_pb, PbO] df.loc[mask_pb, BaO] ) / df.loc[mask_pb, SiO2].replace(0, np.nan) # 高钾玻璃二氧化硅 / 助熔剂氧化物之和列名按实际数据调整 mask_k df[类型] 高钾 df.loc[mask_k, 风化系数] df.loc[mask_k, SiO2] / ( df.loc[mask_k, K2O] df.loc[mask_k, CaO] df.loc[mask_k, Al2O3] ).replace(0, np.nan) # 分组排序找相邻差值的最大跳变作为候选阈值 for (t, w), g in df.groupby([类型, 纹饰]): s g[风化系数].dropna().sort_values().reset_index(dropTrue) if len(s) 3: gap s.diff() i gap.idxmax() print(t, w, 候选阈值区间:, round(s[i-1], 2), ~, round(s[i], 2))replace(0, np.nan)是防止分母为零直接报错或产生无穷大这类异常值一旦混进排序断崖位置会被带偏。分组样本少于五个时跳变没有统计意义不要硬套改回人工核对原始标签。还有一种情况必须提前想好算出来的系数与原始“是否风化”标签冲突怎么办。数据里确实有几个样品对不上处理原则是以定量指标为准、原标签作参考并在论文里写明调整了哪几个、为什么调整。藏着不写答辩时被问一句就露馅。2.3 卡方检验和相关性分析别把p值当结论问题一前半段要回答“表面风化与类型、纹饰、颜色有没有关系”标准做法是相关系数矩阵加卡方检验。很多人算出p小于0.05就写“显著相关说明类型决定风化”这句话会被直接扣分因为p值只回答“有没有关联”不回答“关联有多强”。还要看效应量Phi、Cramers V、列联系数都在0.3上下属于小到中等强度说明风化还受埋藏环境、纹饰对应的年代地域等因素影响。结论句可以照这个模板写表面风化与类型存在边缘显著的关联p约0.049效应量约0.27强度偏弱说明单一因素解释力有限因此后续按玻璃类型分组分别建模是必要的。这一句同时完成了三件事——报了统计量、承认了关联强度、把结论接到了下一步建模上。纹饰和颜色的p值更小、效应量更高可以补一句“纹饰所代表的埋藏环境对风化的影响强于类型本身”正好呼应数据探索里“某种纹饰的样品全部风化”的观察。评委想听的不是“显著”两个字而是你能不能从统计结果里读出建模策略。3. 主成分分析因子分析把风化后的成分反推回风化前3.1 降维不是走过场保留几个主成分要有依据十几种化学成分彼此高度相关二氧化硅和助熔剂氧化物此消彼长直接扔进模型会共线系数估计不稳定。主成分分析在这里的作用就两个压掉共线性把维度降到样本量撑得住的水平。做之前必须标准化否则占比几十的二氧化硅会独吞方差其他成分的贡献被淹没主成分就变成“二氧化硅及其影子”。保留几个主成分要同时给三条依据累计贡献率达到85%以上、特征值大于1、每个主成分能解释出一种业务含义。铅钡玻璃某类纹饰取到8个主成分时累计贡献率超过97%可以写“前8个主成分已覆盖绝大部分变异且第1主成分与二氧化硅高度正相关第2主成分与铅钡氧化物正相关具备可解释性”。只写“取累计贡献率85%”会被追问为什么不取90%三条一起给追问就落不到地上。from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA X df_group[comp_cols].values # 某一类如铅钡-未风化的成分矩阵 X_std StandardScaler().fit_transform(X) # 标准化均值0方差1 pca PCA(n_components0.85) # 先让累计贡献率决定下限 pca.fit(X_std) print(保留主成分数:, pca.n_components_) print(累计贡献率:, pca.explained_variance_ratio_.cumsum()) # 若需要固定个数例如与因子分析对齐再显式指定 pca_k PCA(n_components8).fit(X_std) scores pca_k.transform(X_std) # 主成分得分供聚类与关联度使用n_components0.85是按累计贡献率自动截断输出的个数就是下限真正写进论文的个数还要结合特征值和解释性往上调。scores是降维后的样本坐标后面做亚类聚类、灰色关联度都基于它不要再拿原始成分重复算一遍否则前后结果对不上。失败排查看两个地方标准化后如果某个成分方差接近0说明这类样品里该成分几乎不变应从特征里剔除主成分载荷正负号颠倒不影响结论但论文里描述方向时要按载荷符号统一说。3.2 因子分析求逆拿到从风化后回到风化前的转换矩阵因子分析在这道题里承担的是“后悔药”角色。把风化后的成分向量看成原始成分经过风化演进的结果写成“成分 均值 载荷矩阵 × 公共因子 误差”公共因子就是未风化时的成分结构。载荷矩阵由相关系数矩阵的特征值和特征向量构造每个特征向量乘上对应特征值的平方根拼起来就是载荷矩阵。对载荷矩阵求逆就得到从观测成分反推公共因子的系数矩阵代进去就能算出未风化成分的估计值。这一步有两个坑。第一载荷矩阵不一定是方阵主成分个数少于成分个数时要用伪逆直接求逆会报错。第二反推出来的成分可能出现负值化学上不可能处理办法是负值截断为0再按行归一化到100%并在论文里说明这是模型边界处理。预测结果还要给一组验证拿未风化样品回代算预测值与实测值的相对误差误差在可接受范围内才说明转换矩阵靠谱。import numpy as np X df_group[comp_cols].values mu X.mean(axis0) # 均值向量 Xc X - mu R np.corrcoef(Xc, rowvarFalse) # 相关系数矩阵 eig_val, eig_vec np.linalg.eigh(R) idx np.argsort(eig_val)[::-1] # 特征值从大到小 eig_val, eig_vec eig_val[idx], eig_vec[:, idx] k 8 # 与主成分个数对齐 A eig_vec[:, :k] * np.sqrt(eig_val[:k]) # 载荷矩阵 A_inv np.linalg.pinv(A) # 非方阵用伪逆避免直接求逆报错 factor (Xc) A_inv.T # 反推公共因子未风化成分结构 # 还原到原始量纲并做非负、归一化处理 recon factor A.T mu recon np.clip(recon, 0, None) recon recon / recon.sum(axis1, keepdimsTrue) * 100np.linalg.eigh专门用于对称矩阵相关系数矩阵一定对称用它比通用求解器稳。pinv是伪逆k小于成分数时这是唯一能跑通的写法。np.clip那行是负值处理别省略否则归一化后会出现更离谱的数。最后的归一化按行做保证每个样品所有成分加起来还是100%这一步不做后面和实测值比较时量纲对不上。3.3 答辩和面试最常追问的三个问题第一个问题几乎必问为什么不用回归直接拟合风化前后成分。标准答法是——回归需要配对样本也就是同一件文物风化前和风化后的两组观测实际数据里只有风化后的检测值和少量未风化样品配对关系不存在因子分析利用的是变量之间的相关结构不需要配对标签适合这种“只有结果、没有起点”的场景。答完补一句“如果后续能获得同一墓葬出土的风化与未风化配对样品可以用回归做交叉验证”显得你知道方法的边界。第二个问题主成分个数凭什么定8个。答法是把三条依据报一遍再加一句敏感性——取6个和取10个分别重跑看结论方向是否改变如果关联度排序基本稳定说明个数选择不影响主要结论。第三个问题反推结果出现负值或者超过100怎么解释。答法要坦诚这是线性模型外推的边界效应处理方式是截断加归一化同时在论文里给出预测区间提醒使用者不要把它当成精确测量值。面试官最想听的不是“我的模型很完美”而是你知道它在哪里会失效。4. 亚类划分、未知样品鉴别与敏感性分析问题二到问题三的一条线4.1 区间法去极值、取区间、给归属规则问题二要把高钾、铅钡两大类别再往下切亚类。做法是先按“类型是否风化”分成四组每组单独做主成分分析挑贡献率高的几个成分去掉最大值和最小值之后用剩余数据的上下界构成区间。新样品的成分落进区间、或者离区间足够近就归入该类。这套方法不依赖分布假设样本量小也能用缺点是区间边界对极端值敏感所以去极值这一步不能省。四组区间可以整理成一张表放进论文评委一眼就能看到划分依据。表里只放贡献率最高的几个成分即可全放反而看不清重点。归属规则里的“靠近”必须量化比如规定“超出区间不超过区间宽度的10%视为靠近”否则评委会追问到底多近算近。这一步的表述越具体问题三的可复现性越强。类别二氧化硅区间氧化钾区间氧化钙区间其他主要成分区间铅钡-未风化50.61 ~ 69.710.15 ~ 0.320.38 ~ 2.82氧化铝 1.44 ~ 13.65高钾-未风化59.81 ~ 69.337.37 ~ 12.534.71 ~ 8.23氧化镁 0.52 ~ 1.77铅钡-风化4.61 ~ 37.360.14 ~ 1.050.37 ~ 4.49氧化镁 0.47 ~ 1.18高钾-风化92.35 ~ 95.020.59 ~ 1.01—样本少仅取主成分注意区间上界与聚类得到的聚点如果不一致先回头检查去极值是不是做在了错误的组里别急着改区间。4.2 最短距离法聚类一维数据也能切出亚类区间法给出的是“属于哪一大类”亚类还需要再切一刀。数据里差异最大的是二氧化硅把它单独拿出来做一维聚类用最短距离法把样品并成两类得到两个聚点取两聚点的中点作为分割点。铅钡风化类的两个聚点大约在33.96和15.63中点24.80正好把区间分成两段和区间法的结果能对上这种互相印证要写进论文是合理性的证据。from scipy.cluster.hierarchy import linkage, fcluster x df_group[SiO2].values.reshape(-1, 1) # 最短距离法类间距离取两类最近样本的距离 Z linkage(x, methodsingle, metriceuclidean) labels fcluster(Z, t2, criterionmaxclust) # 切成两类 # 两个聚点取各类均值分割点取中点 c1, c2 x[labels 1].mean(), x[labels 2].mean() split (c1 c2) / 2 print(聚点:, round(c1, 2), round(c2, 2), 分割点:, round(split, 2))methodsingle就是最短距离法maxclust控制切几类。一维数据用层次聚类比kmeans稳因为不需要预设初始中心结果可复现。切完之后一定要画分布图看两类比例如果一类占九成、另一类只有一两个样品说明这个亚类划分没有实际意义要换成分或者换聚类数。合理性论证就靠这张比例图加区间互相对照两句话说完即可不用长篇大论。4.3 3%扰动敏感性分析两问结论方向不同别写反敏感性分析的做法是给原始数据加3%的随机扰动重新跑一遍分类看标签变了多少。这里有个极容易写反的地方问题二的亚类划分扰动后有几组的标签变化比例超过20%结论应该是“亚类划分对数据扰动较敏感”问题三的未知样品鉴别扰动后归属类型基本不变结论是“类型鉴别稳健”。一个敏感一个稳健方向相反因为亚类边界本来就是在小样本上人为切出来的而大类之间的成分差异足够大。论文里要给出判定标准不能只说“变化不大”。可以约定标签变化率低于5%为稳健5%到20%为中等超过20%为敏感。再补一句原因分析——亚类敏感不代表模型错而是说明亚类边界不宜作为硬性鉴定依据实际应用中以大类鉴别为准。这句话既是解释也是自我保护评委追问时直接引用。做扰动实验要固定随机种子重复多次取平均变化率只跑一次就下结论结果不可复现。4.4 灰色关联度五步问题四照这个框架写就不会跑偏问题四要分析不同类别内部化学成分之间的关联关系再比较类别之间的差异。灰色关联度分析是标准工具五步固定确定参考数列和比较数列、无量纲化、算关联系数、算关联度、排序。分歧点在第一歩——参考数列选什么。一种写法是把样品序号作为参考、各成分作为比较看哪种成分跟整体水平更同步另一种是把某个主成分或综合得分作为参考、各成分作为比较看成分对整体结构的贡献。两种都能做但必须在论文里说清楚选了哪种、为什么选否则结果没法解释。def grey_relational(ref, comp, rho0.5): ref: 参考数列 (n,)comp: 比较数列矩阵 (m, n) ref np.asarray(ref, dtypefloat) comp np.asarray(comp, dtypefloat) # 初值化无量纲处理消除百分比量纲差异 ref_n ref / ref[0] comp_n comp / comp[:, [0]] diff np.abs(ref_n[None, :] - comp_n) dmin, dmax diff.min(), diff.max() # 关联系数rho取0.5是常用分辨系数 xi (dmin rho * dmax) / (diff rho * dmax) return xi.mean(axis1) # 按行求均值得到关联度 # 示例以第一主成分得分为参考各主要成分为比较数列 scores pca_k.transform(StandardScaler().fit_transform(X))[:, 0] degree grey_relational(scores, X[:, :5].T) print(dict(zip(comp_cols[:5], degree.round(3))))rho取0.5是惯例取值越小分辨力越强论文里写一句“取0.5兼顾分辨力与稳定性”即可。初值化用第一个样品作基准如果第一个样品某成分接近0比值会爆炸这时改用均值化更稳。关联度排序出来后结论要克制四类里二氧化硅关联性最强、其他成分差别不大就如实写“类别间差异有限”不要硬编出“某某成分显著主导”。灰色关联度本身是灰色系统里的近似方法过度解读反而暴露理解不到位。5. 易错点五个最容易在阅卷时被扣分的操作5.1 把“是否风化”当因变量做监督分类错在哪拿到数据就直接用“是否风化”当标签训练分类器然后宣称完成了问题一。为什么错问题一要回答的是风化与类型、纹饰、颜色的统计关系以及反推未风化成分分类器只能告诉你“能不能分开”回答不了“风化前后成分怎么变”而且样本量只有几十条交叉验证都做不稳。怎么答对先用卡方和相关系数回答关系问题再用因子分析做反推分类模型最多作为补充验证放在讨论里。顺序反了整道题的方向就偏了。5.2 两类玻璃共用一套风化系数公式错在哪图省事写一个统一公式把铅钡和高钾放在一起算系数、定阈值。为什么错两种玻璃的风化机理相反高钾玻璃风化后二氧化硅占比升高铅钡玻璃风化后二氧化硅占比下降同一个比值在两类样品里含义完全相反阈值更不可能通用。怎么答对分开定义、分开排序、分开定阈值论文里把两个公式并列写出来各配一句化学解释。这一步多花十分钟能省掉答辩时十分钟的追问。5.3 主成分分析前不标准化或者四类混在一起做错在哪直接把原始百分比矩阵丢进主成分分析或者把四类样品合并成一个大矩阵降维。为什么错不标准化会让占比最大的二氧化硅主导所有主成分其他成分的信息被压掉四类样品成分分布差异极大混在一起降维主成分解释的是“类别差异”而不是“类内结构”后面的区间划分全乱。怎么答对先标准化再按“类型是否风化”四组分别做每组单独报告贡献率和载荷。分组做还有一个好处问题三的未知样品可以直接套用对应组的区间。5.4 敏感性分析的结论写反或者没给判定标准错在哪问题二扰动后标签变化很大却写“模型稳健”或者只写“变化不大”不给出变化率数字。为什么错结论和数字对不上评委一核对就发现没有判定标准“不大”就是主观判断无法复现。怎么答对先定阈值5%、20%两档再报每组的标签变化率最后按阈值给结论。问题二写敏感、问题三写稳健并解释原因——亚类边界是小样本人为划分大类差异是化学本质差异。5.5 灰色关联度里参考数列和比较数列颠倒错在哪把成分当参考、样品当比较算出来的关联度排序和论文结论对不上或者换个数据集结果完全变样。为什么错参考数列代表“系统行为特征”比较数列代表“影响因素”位置颠倒后关联系数的物理含义就变了排序自然不可比。怎么答对先明确这一问要解释什么——如果是成分对整体结构的贡献就用主成分得分或综合得分作参考、成分作比较如果是个体差异就反过来。选定之后在论文里写清楚四类样品用同一套规则保证横向可比。6. 答题节奏四问的时间分配与答辩三句话6.1 时间怎么切把力气压在能拿分的地方建模时间通常三天四问的投入不能平均。问题一的数据清洗和风化系数是地基值得花掉三成时间因为后面三问全部依赖这套清洗结果和分组问题二的主成分加区间划分再占三成这一问的图表最多、最容易做出视觉效果问题三本质是套用问题二的区间做判别两成时间足够问题四的灰色关联度有固定五步一成时间就能跑完剩下的时间留给摘要和图表美化。摘要里必须出现风化系数、主成分、区间划分、敏感性、灰色关联度这几个词评委翻摘要就在找它们。环节建议占比必须交付的内容常见扣分数据清洗与探索30%表单关联说明、风化系数阈值依据阈值来源没写主成分与因子分析30%贡献率表、载荷矩阵、反推结果未标准化、未分组亚类划分与鉴别20%四组区间表、聚类分割点归属规则模糊灰色关联度10%五步流程、关联度排序参考数列选错摘要与图表10%每问一句结论加一个关键数字摘要无数字6.2 答辩三句话把追问挡在结论外面答辩或面试被追问时准备三句话轮流用。第一句讲方法边界“这个方法在样本量小、没有配对标签的条件下是合适的选择如果换成配对数据我们会改用回归做交叉验证。”第二句讲结论强度“统计上显著但效应量偏小所以我们的结论是分组建模有必要而不是类型单独决定风化。”第三句讲失败处理“反推出现负值时我们做了截断和归一化并在论文里标注了这是模型外推的边界。”三句话的共同点是不吹模型、先认边界评委听完通常会转去问下一个问题而不是继续深挖。从那以后我每次带人做这类题都强制先走一遍“清洗—分组—定阈值—再建模”的顺序谁想跳过阈值直接上模型就让他把排序后的系数表打印出来自己看断崖在哪。这套流程不玄学就是多花半天把地基打牢后面每一步都能少返工一次。希望帮到你。本文还有配套的精品资源点击获取