1. 项目概述一次从混沌到清晰的解题之旅2022年的全国大学生数学建模竞赛C题题目是“古代玻璃制品的成分分析与鉴别”。当时拿到这个题目很多队伍的第一反应可能是懵的——这看起来既不像传统的优化问题也不像典型的预测模型它横跨了化学、材料学、考古学和数据分析多个领域。我记得当时我们团队内部也有过短暂的迷茫但很快我们就意识到这道题的精髓恰恰在于其“跨界”特性它考察的不是某个单一的数学模型而是一套完整的数据驱动型问题解决框架从数据预处理、特征工程、模型构建到结果解释环环相扣。最终我们凭借一套清晰的解题逻辑拿到了不错的成绩。今天我就把这套从实战中打磨出来的思路、踩过的坑以及核心代码实现毫无保留地分享出来。无论你是正在备赛的新手还是对数据科学感兴趣的朋友这篇文章都能为你提供一个从原始数据到最终论文的完整复现路径。简单来说这道题给了我们两组数据一批已知分类的玻璃文物化学成分数据训练集和一批待分类的未知玻璃文物数据测试集。我们的核心任务就是根据化学成分的差异对未知文物进行科学分类高钾玻璃或铅钡玻璃并进一步分析其亚类以及风化前后的成分规律。这本质上是一个有监督的分类问题但夹杂着大量的数据清洗、统计分析和机理探索。它适合所有对数学建模、机器学习、数据分析以及跨学科应用感兴趣的同学。接下来我将彻底拆解我们的解题全流程。2. 解题整体设计与核心思路拆解面对这样一个多任务、数据质量不高的题目盲目套模型是死路一条。我们的整体设计遵循了“先治理后分析先探索后建模先验证后深挖”的原则。2.1 问题界定与任务分解首先我们必须把赛题模糊的描述转化为清晰、可执行的数据科学任务。题目要求可以分解为三个核心子问题主成分分析与相关性探索分析化学成分之间的关联以及化学成分与玻璃类型、风化状态之间的关系。这属于描述性统计和探索性数据分析EDA的范畴。玻璃类型判别根据化学成分对未知样品进行高钾/铅钡玻璃的分类。这是典型的二分类问题。亚类划分与风化分析对分类后的玻璃进行合理的亚类划分并分析风化对化学成分的影响规律。这涉及到无监督学习聚类和差异性检验。我们的思路是任务1为任务2和3提供特征理解和筛选依据任务2是核心建模环节任务3是在任务2结果上的深化分析。三者呈递进关系而非并列。2.2 技术路线选型与理由为什么选择这样的技术路线这是基于数据特性和问题目标深思熟虑的结果。数据预处理优先原始数据存在大量缺失值记为“NaN”或“-”且各成分量纲差异巨大SiO2含量可达70%而某些微量元素含量不足1%。不处理缺失值大多数模型无法运行不进行标准化模型会被大数值特征主导。因此数据清洗和标准化是无可争议的第一步。EDA引导特征工程在建模前我们通过相关性分析、箱线图、主成分分析PCA等方法直观地看数据。例如通过PCA我们发现前两个主成分就能解释大部分方差并且能在二维图上较好地区分两类玻璃这增强了我们使用成分数据进行分类的信心。同时EDA能帮助我们发现异常点为后续模型稳定性提供保障。分类模型选择逻辑回归与支持向量机SVM为主逻辑回归模型简单可解释性强。我们可以直接得到每个化学成分对“是铅钡玻璃”这一事件的贡献度系数这对于回答“哪些成分是主要鉴别因素”至关重要。这是论文中的一个亮点。支持向量机SVM特别是线性SVM在特征维度不高、样本量不大的情况下寻找最大间隔分类超平面通常效果稳定且泛化能力好。我们将其作为与逻辑回归对比的模型。为什么不直接用复杂的树模型或神经网络样本量有限训练集仅数十个复杂模型极易过拟合。且赛题强调“分析”需要模型具有一定的可解释性黑箱模型在此处不占优势。亚类划分使用聚类分析K-Means在完成大类分类后对每一类玻璃内部进行聚类寻找自然分组。选择K-Means是因为其经典、高效且对于这种连续型特征的数据表现良好。关键点在于如何确定最佳的聚类数目K我们会使用手肘法和轮廓系数共同判定。风化分析采用统计检验比较风化与未风化玻璃在各项成分上的差异使用Mann-Whitney U检验非参数检验而非T检验。因为部分成分数据可能不服从正态分布且样本量小非参数检验更加稳健。注意整个解题过程不是单向流水线而是一个循环迭代的过程。例如在建模后发现分类效果不佳可能需要返回去重新审视特征工程或数据预处理步骤。3. 核心环节一数据预处理与探索性分析这是所有工作的基石也是最容易出错、最考验耐心的地方。我们花了近三分之一的时间在这里。3.1 数据清洗实战原始数据通常以Excel或CSV格式提供。我们使用Python的Pandas库进行操作。import pandas as pd import numpy as np # 读取数据 train_data pd.read_excel(附件表单1.xlsx) # 已知类别的训练集 test_data pd.read_excel(附件表单2.xlsx) # 待分类的测试集 # 1. 缺失值处理 # 首先查看缺失情况 print(train_data.isnull().sum()) # 发现缺失值多为微量元素且可能是“未检测到”而非“数据丢失”。 # 策略对于化学成分用0或一个极小的值如检测下限的一半填充表示“含量极低可忽略”。 # 这里采用0填充因为后续标准化会处理量纲且0具有明确的物理意义“无”。 train_data_filled train_data.fillna(0) test_data_filled test_data.fillna(0) # 2. 数据格式化 # 确保所有成分列为数值型 component_columns [SiO2, Na2O, K2O, ...] # 所有化学成分列名 for col in component_columns: train_data_filled[col] pd.to_numeric(train_data_filled[col], errorscoerce) test_data_filled[col] pd.to_numeric(test_data_filled[col], errorscoerce) # 再次用0填充因转换产生的NaN train_data_filled[component_columns] train_data_filled[component_columns].fillna(0) test_data_filled[component_columns] test_data_filled[component_columns].fillna(0) # 3. 特征与标签分离 X_train_raw train_data_filled[component_columns] y_train train_data_filled[类型] # ‘高钾’ or ‘铅钡’ X_test_raw test_data_filled[component_columns]实操心得关于缺失值填充我们尝试过均值填充、中位数填充和0填充。最终选择0填充是因为在化学成分分析中未检出通常意味着含量低于仪器检测限近似为0。这在物理意义上最合理并且简化了模型。一个关键的检查步骤填充后务必计算每个样品的成分总和。玻璃的主要成分总和应接近100%允许有少量误差。如果某个样品总和偏差巨大如远大于110%或低于80%说明数据可能存在严重错误或异常需要单独审查。3.2 特征工程与标准化化学成分数据是典型的“成分数据”其特征是所有特征之和为常数约100%。这会导致“定和约束”使特征间存在虚假的负相关性。一种常见的处理方法是进行中心对数比变换但为了简化并优先保证模型可解释性我们采用了更通用的方法。from sklearn.preprocessing import StandardScaler # 方法一直接标准化Z-Score # 优点消除量纲使所有特征均值为0方差为1。 scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train_raw) X_test_scaled scaler.transform(X_test_raw) # 注意使用训练集的均值和方差来转换测试集 # 方法二备选归一化Min-Max Scaling # from sklearn.preprocessing import MinMaxScaler # scaler MinMaxScaler() # X_train_scaled scaler.fit_transform(X_train_raw) # X_test_scaled scaler.transform(X_test_raw)为什么测试集要使用训练集的scaler这是机器学习中的基本原则——信息泄露。测试集模拟的是未知的、未来的数据。我们必须假设在“预测时”是不知道测试集整体分布的因此标准化参数必须完全从训练集中学习。用fit_transform处理训练集用transform处理测试集。3.3 探索性数据分析可视化我们用图表来“感受”数据这部分结果可以直接放入论文。import matplotlib.pyplot as plt import seaborn as sns from sklearn.decomposition import PCA # 1. 相关性热力图 plt.figure(figsize(12, 10)) corr_matrix pd.DataFrame(X_train_scaled, columnscomponent_columns).corr() sns.heatmap(corr_matrix, annotTrue, fmt.2f, cmapcoolwarm, center0) plt.title(化学成分相关性热力图) plt.tight_layout() plt.show() # 从热力图可以发现K2O和PbO可能呈现强负相关这符合高钾玻璃和铅钡玻璃的化学定义。 # 2. PCA降维可视化 pca PCA(n_components2) X_train_pca pca.fit_transform(X_train_scaled) plt.figure(figsize(8, 6)) colors {高钾: red, 铅钡: blue} for glass_type in colors.keys(): idx y_train glass_type plt.scatter(X_train_pca[idx, 0], X_train_pca[idx, 1], labelglass_type, colorcolors[glass_type], alpha0.7) plt.xlabel(fPC1 (方差解释度: {pca.explained_variance_ratio_[0]:.2%})) plt.ylabel(fPC2 (方差解释度: {pca.explained_variance_ratio_[1]:.2%})) plt.title(训练集PCA二维投影) plt.legend() plt.grid(True) plt.show() # 如果两类点在PC1-PC2平面上能较好分离说明成分数据确实携带了强烈的分类信息。踩坑记录最初我们曾尝试对原始百分比数据直接做PCA结果发现前几个主成分完全被SiO2这种高含量成分主导其他重要鉴别成分如PbO, K2O的贡献被淹没。这再次证明了标准化的重要性。标准化后每个特征都被赋予了同等的重要性单位方差PCA才能公平地评估所有成分的贡献。4. 核心环节二玻璃类型判别模型构建与优化这是解题的核心我们构建了多个模型进行对比并深入进行了模型解释。4.1 逻辑回归模型逻辑回归不仅能分类还能给出概率和特征重要性。from sklearn.linear_model import LogisticRegression from sklearn.model_selection import cross_val_score, StratifiedKFold from sklearn.metrics import classification_report, confusion_matrix, accuracy_score # 将标签转化为数值 y_train_numeric y_train.map({高钾: 0, 铅钡: 1}) # 创建并训练逻辑回归模型 # 参数C是正则化强度的倒数较小的C值意味着更强的正则化可以防止过拟合。 lr_model LogisticRegression(C1.0, penaltyl2, solverliblinear, random_state42, max_iter1000) lr_model.fit(X_train_scaled, y_train_numeric) # 交叉验证评估模型稳定性 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) cv_scores cross_val_score(lr_model, X_train_scaled, y_train_numeric, cvcv, scoringaccuracy) print(f逻辑回归5折交叉验证准确率: {cv_scores.mean():.4f} (/- {cv_scores.std()*2:.4f})) # 查看模型在训练集上的表现仅供参考交叉验证更可靠 y_train_pred lr_model.predict(X_train_scaled) print(训练集分类报告:) print(classification_report(y_train_numeric, y_train_pred, target_names[高钾, 铅钡]))模型解释与特征重要性逻辑回归最大的优势在此体现。# 获取特征系数 coefficients lr_model.coef_[0] feature_importance pd.DataFrame({ 特征: component_columns, 系数: coefficients, 绝对值系数: np.abs(coefficients) }).sort_values(by绝对值系数, ascendingFalse) print(逻辑回归特征系数按绝对值排序:) print(feature_importance) # 系数为正意味着该成分含量越高越倾向于预测为‘铅钡玻璃’标签1。 # 系数为负则意味着该成分含量越高越倾向于预测为‘高钾玻璃’标签0。 # 例如PbO的系数很可能为正且绝对值很大K2O的系数很可能为负且绝对值很大这与化学常识完全吻合。4.2 支持向量机模型我们尝试了线性核SVM并与逻辑回归对比。from sklearn.svm import SVC svm_model SVC(kernellinear, C1.0, random_state42) # 线性核 svm_model.fit(X_train_scaled, y_train_numeric) svm_cv_scores cross_val_score(svm_model, X_train_scaled, y_train_numeric, cvcv, scoringaccuracy) print(fSVM(线性核)5折交叉验证准确率: {svm_cv_scores.mean():.4f} (/- {svm_cv_scores.std()*2:.4f}))模型对比与选择通过交叉验证准确率、稳定性和可解释性综合判断。通常逻辑回归和线性SVM在这个规模的数据集上表现接近。如果逻辑回归的交叉验证准确率与SVM相当甚至略高我们倾向于选择逻辑回归因为其系数解释更直接。如果SVM明显更优则选择SVM。在我们的实际解题中两者相差无几都在94%-96%左右因此我们最终选择逻辑回归作为主模型以便在论文中详细分析成分贡献。4.3 对测试集进行预测使用训练好的最佳模型对未知样品进行分类。# 假设我们最终选定逻辑回归模型 best_model lr_model # 预测测试集类别 y_test_pred_numeric best_model.predict(X_test_scaled) y_test_pred pd.Series(y_test_pred_numeric).map({0: 高钾, 1: 铅钡}) # 预测测试集属于各类别的概率对于逻辑回归 y_test_prob best_model.predict_proba(X_test_scaled) # y_test_prob 是一个二维数组第一列是类别0高钾的概率第二列是类别1铅钡的概率 # 将预测结果保存 test_data_filled[预测类型] y_test_pred test_data_filled[预测为铅钡玻璃的概率] y_test_prob[:, 1] test_data_filled[[文物编号, 预测类型, 预测为铅钡玻璃的概率]].to_excel(测试集预测结果.xlsx, indexFalse)注意事项predict_proba方法只有部分分类器具备如逻辑回归。SVC默认不支持需要设置probabilityTrue参数但这会显著增加计算时间。如果只需要分类标签用predict即可。5. 核心环节三亚类划分与风化规律分析完成大类判别后我们需要对结果进行更精细的解读。5.1 基于聚类分析的亚类划分我们将预测后的测试集按“高钾”和“铅钡”两类分开分别进行聚类分析。from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score # 分离数据 high_k_test X_test_scaled[y_test_pred 高钾] lead_barium_test X_test_scaled[y_test_pred 铅钡] def find_optimal_k(data, max_k8): 寻找最优聚类数K inertias [] silhouette_scores [] K_range range(2, max_k1) for k in K_range: kmeans KMeans(n_clustersk, random_state42, n_init10) kmeans.fit(data) inertias.append(kmeans.inertia_) # 肘部法则簇内误差平方和 if len(data) k: # 轮廓系数要求聚类数小于样本数 silhouette_scores.append(silhouette_score(data, kmeans.labels_)) else: silhouette_scores.append(None) # 绘制肘部法则图 plt.figure(figsize(12,4)) plt.subplot(1,2,1) plt.plot(K_range, inertias, bo-) plt.xlabel(聚类数 K) plt.ylabel(簇内误差平方和 (Inertia)) plt.title(肘部法则图) plt.grid(True) # 绘制轮廓系数图 plt.subplot(1,2,2) valid_scores [s for s in silhouette_scores if s is not None] valid_K K_range[:len(valid_scores)] plt.plot(valid_K, valid_scores, ro-) plt.xlabel(聚类数 K) plt.ylabel(轮廓系数 (Silhouette Score)) plt.title(轮廓系数图) plt.grid(True) plt.tight_layout() plt.show() # 综合建议选择轮廓系数最高且inertia下降趋势变缓的K print(f轮廓系数: {dict(zip(valid_K, valid_scores))}) return valid_K, valid_scores, inertias # 对铅钡玻璃测试集寻找最优K print(--- 铅钡玻璃亚类划分 ---) k_range, scores, inertias find_optimal_k(lead_barium_test, max_k6) # 假设我们根据图表选择 K3 k_optimal_lead 3 kmeans_lead KMeans(n_clustersk_optimal_lead, random_state42, n_init10).fit(lead_barium_test) # 将聚类标签赋回原数据框 lead_indices test_data_filled[y_test_pred 铅钡].index test_data_filled.loc[lead_indices, 亚类] kmeans_lead.labels_ # 对高钾玻璃进行类似操作...实操心得聚类分析没有绝对正确的答案。肘部法则看的是inertia下降的拐点轮廓系数越接近1越好。需要结合两者并考虑实际的化学或考古学意义。例如如果K3和K4的轮廓系数相近但K3时每个簇的样本量更均衡且能对应到不同的化学成分模式如高铅、中铅低钡、高钡等那么K3可能是更合理、更可解释的选择。一定要在论文中阐述你选择K值的理由。5.2 风化影响的统计分析我们利用训练集中有风化标识的样品分析风化前后成分变化。from scipy.stats import mannwhitneyu # 假设训练集数据框中有‘风化’列取值为‘风化’或‘无风化’ weathered_data train_data_filled[train_data_filled[风化] 风化][component_columns] unweathered_data train_data_filled[train_data_filled[风化] 无风化][component_columns] significant_changes [] for component in component_columns: # 进行Mann-Whitney U检验 stat, p_value mannwhitneyu(weathered_data[component].dropna(), unweathered_data[component].dropna(), alternativetwo-sided) if p_value 0.05: # 显著性水平设为0.05 median_weathered weathered_data[component].median() median_unweathered unweathered_data[component].median() change_rate (median_weathered - median_unweathered) / median_unweathered if median_unweathered ! 0 else np.inf significant_changes.append({ 成分: component, p值: p_value, 风化中位数: median_weathered, 未风化中位数: median_unweathered, 变化率: change_rate }) # 将显著变化的成分整理成表格 significant_df pd.DataFrame(significant_changes).sort_values(byp值) print(风化前后成分显著变化分析:) print(significant_df)结果解读通过这个分析我们可以得出诸如“风化后K2O含量显著降低而SiO2相对含量升高”这样的结论。在论文中可以用表格清晰展示这些显著变化的成分及其变化方向并结合化学知识如钾离子易淋溶进行机理解释这是论文的另一个加分点。6. 常见问题、避坑技巧与实战心得回顾整个解题过程我们遇到了不少典型问题以下是总结出的“避坑指南”。6.1 数据预处理中的陷阱陷阱一忽视成分数据的“定和约束”。直接对原始百分比做相关性分析或PCA结果可能是扭曲的。我们的处理方式是优先采用标准化后的数据进行分析和建模这在一定程度上缓解了问题。更严谨的做法是使用中心对数比变换但需注意变换后数据中存在负无穷值当成分为0时的处理。陷阱二测试集标准化参数错误。务必记住StandardScaler的fit只用在训练集上然后用这个scaler去transform训练集和测试集。用训练集fit_transform测试集只transform。这是一个会导致模型评估完全失真的低级错误。陷阱三异常点处理不当。在EDA阶段通过箱线图发现某个样品的PbO含量奇高无比。我们没有直接删除而是先检查了文物编号和背景知识。确认可能是特殊工艺制品后我们选择保留该样本但在后续聚类分析时注意到它自成一类并对此进行了说明。不要轻易删除异常点要探究其产生原因。6.2 模型选择与评估的误区误区一盲目追求高精度模型。在训练集上复杂的模型如带RBF核的SVM或随机森林可能达到100%准确率但这极有可能是过拟合。通过5折或10折交叉验证你会发现它们的泛化性能可能还不如简单的逻辑回归。小样本下模型简单就是美。误区二忽略模型的可解释性。国赛评阅看重解题思路和过程。逻辑回归的系数、SVM的支持向量、决策树的特征重要性都是你可以用来“讲故事”的素材。在论文中清晰地展示“是PbO和K2O的含量对比起到了决定性作用”比单纯给出一个高精度黑箱模型更有说服力。误区三未考虑类别不平衡。虽然本题两类样本大致平衡但养成检查的习惯很重要。如果发现类别不平衡如铅钡玻璃样本远多于高钾玻璃在逻辑回归中可以使用class_weightbalanced参数在SVM中也可以设置class_weight让模型更关注少数类。6.3 论文写作与结果呈现要点要点一图表清晰专业。热力图、PCA散点图、聚类结果图、成分变化条形图……一图胜千言。确保图表有清晰的标题、坐标轴标签、图例。使用专业的配色如Set2, Set3, tab20c等色盲友好配色。要点二分析层层递进。论文结构应反映你的解题思路数据预处理 - 探索性分析 - 模型建立与对比 - 亚类划分 - 风化分析。每一部分的分析结论都应为下一部分提供依据。要点三结论明确回答赛题。最终结论要直接回应题目问的所有问题1. 哪些成分关联性强2. 未知玻璃如何分类给出具体编号和类别3. 亚类如何划分给出划分标准和各类特征4. 风化规律是什么列出显著变化的成分及变化方向。避免使用“可能”、“大概”等模糊词汇基于数据和模型给出肯定结论。要点四代码与论文分离。论文正文中只展示关键的流程图、结果图和核心表格。将完整的代码、数据清洗过程、模型参数调优细节等以附录或附件形式提交。保持正文的简洁和可读性。最后我想分享一点最深的体会数学建模竞赛尤其是国赛比拼的从来不是谁用的模型最高深、最时髦而是谁对问题的理解更透彻谁的解题逻辑更严谨谁的故事讲得更完整。从一团乱麻的数据中通过一步步扎实的分析抽丝剥茧得出令人信服的结论这个过程本身带来的成就感远比奖项更重要。希望这份详细的复盘能帮助你少走弯路更自信地面对未来的挑战。