2022年数学建模国赛C题(古文物风化)助攻论文&支撑材料.doc_2022年国赛c题题目下载

本博客下载链接包含修改的word版本, 可免费下载阅览学习, 也可作为数学建模相关课程作业修改上交:

链接：https://pan.baidu.com/s/1HxzDk3q0p6y2xpuJyxPgvw?pwd=qtnc

提取码：qtnc

       古代玻璃在风化过程中，由于内部元素与环境元素进行大量交换，导致其成分比例发生变化，从而影响对其类别的正确判断，本文通过对已知数据进行统计性分析进行相关的分类研判并采用无监督学习的K-means++算法对玻璃制品成分进行分析与鉴别

       首先进行数据预处理，根据题目要求，删除编号为15和17的两条错误数据，将颜色为空值的部分根据风化程度规律设定为黑色处理，将其他空值进行填“0”处理，进行接下来的计算

       针对问题一，首先针对四个定类变量使用卡方检验进行分析，根据显著性P值是否小于0.05判断出玻璃类型和颜色与表面风化存在显著性差异，与纹饰不存在显著性差异，在此基础上进行效应量化分析，其中包含phi、Crammer's V、列联系数与lambda。接着针对不同玻璃类型风化前后的统计规律进行对比分析，首先对风化前后的化学含量指标进行描述性统计分析，接着绘制风化前后频率分布直方图，观察数据的波动情况，最终得出高钾玻璃在风化后主要化学成分含量呈下降趋势；铅钡类玻璃在风化后主要化学成分含量呈上升趋势，根据风化点检测数据的特殊性，建立加权平均预测占比模型，通过采用标准正态分布函数进行赋权，使数据分布较为均匀，计算出风化前后不同种类玻璃的化学成分含量所占的比例的线性映射关系，最终预测出风化前的化学成分含量

       针对问题二，首先将玻璃类型仅分为风化前和风化后，然后使用系统聚类方法进行聚类分析，分析当聚类数为2的时候两种聚类类别属于高钾玻璃和铅钡玻璃的情况，最终计算出在高钾玻璃风化前后对比中，仅有11号和48号判断错误，准确率为93.75%；铅钡玻璃风化前后对比中，准确率为100%。在此基础上建立亚类划分分类研判模型，使用使用K-Means算法对不同类型玻璃是否风化4种类别进行化学指标的选择，分别讨论聚类数为2和3时，每一次聚类结果中的差异性，根据显著性P值得出相对合适的指标进行亚类划分，显著性越强说明模型的合理性好，通过欧式聚类进行判断样本的类别。最后进行灵敏度分析，对数据进行扰动，重新判断，最终得出当扰动范围为30%时，模型的准确率为91%。

       针对问题三，使用问题二中的模型进行分类，给出分类结果，在此基础上对数据进行扰动处理，将变量值在范围（-130%，+130%）内进行随机缩放，观察结果与未进行扰动时结果的对比情况，最终得出扰动后，模型的准确率为100%，说明了模型具有良好的鲁棒性以及普适性。

       针对问题四，建立玻璃种类分类判研判模型，采用无监督学习的kmeans++算法进行求解，聚类范围选择2-5，使用启发式方法找到K均值聚类的质心种子，采用C-H值选择最佳聚类个数，判断出不同类型玻璃风化前后的聚类情况以及各个化学成分所属类别的变化情况，选择CHI、DBI、轮廓系数进行比较分析高钾玻璃与铅钡玻璃之间的差异性，最终计算出高钾类玻璃CHI指数风化后明显增加，铅钡类玻璃风化后相对减少；针对DBI指数，高钾类玻璃DBI指数增加，铅钡玻璃风化后减少；针对轮廓系数，分析得出铅钡玻璃风化后化学成分关联关系较强，而高钾玻璃风化后化学成分关联关系相对减弱。

关键词：玻璃成分鉴别 ；卡方检验 ；加权平均预测 ；系统聚类 ；K-Means++

正文思路展示(部分)

对本文提出的四个相关问题，我们逐一做如下分析：

问题一的分析：首先需要对玻璃表面风化情况与玻璃类型，纹饰和颜色的差异性进行分析，并结合玻璃的类型分析化学成分含量的变化规律以及预测风化前的化学成分含量。

第一步分析：针对定类变量进行卡方检验分析确定自变量与因变量之间的关系，带入SPSS软件中进行求解。分析显著性p值是否小于0.05，进而分析差异性关系

第二步分析：分别讨论铅钡玻璃与高钾玻璃风化前后的变化差异进行描述性统计分析、频率直方图统计分析、正态分布检验等，总结变化情况。

第三步分析：根据风化前后的数据规律，总结出各个化学成分的变化情况，使用加权平均法找到其映射关系并预测风化前的含量。

问题二的分析：需要我们针对高架玻璃和铅钡玻璃进行分类以及亚类划分，并分析模型的合理性和敏感性。

第一步分析：针对高钾玻璃和铅钡玻璃不同化学成分的数值进行统计，找到其具有代表性的化学指标的变化情况作为分类的依据。

第二步分析：在此基础上进行亚类划分，观察化学成分在风化前后的变化情况，颜色变化，纹理变化等，并给出相应的分类依据

第三步分析：在此基础上对数据进行扰动处理（灵敏性检验），并给出相关的合理性依据。

问题三的分析：需要我们对表单3中未知玻璃文物的化学成分进行分析，并预测其所属的类型，并进行敏感性分析。

第一步分析：将表单三中的数据中有无风化的情况进行分类讨论，结合问题2中模型的结论，对表单三中不同类型的玻璃进行分类研判，分析模型的鲁棒性。

第二步分析：将某一类化学元素含量增加一个扰动（-5%，-10%，10%，20%）带入问题2的模型中，观察分类情况是否会变化，并给出模型的稳定性结论。

问题四的分析：需要我们针对不同类别的玻璃样品分析化学成分之间的关联关系。

第一步分析：需要我们针对高架玻璃和铅钡玻璃进行分类，使用无监督学习的K-Means++聚类算法进行聚类分析，根据聚类结果分析化学成分之间的关系。 

第二步分析：采用不同评价指标（CHI、DBI、轮轮廓系数等）比较不同类别之间的化学成分关联关系的差异性

 

附录一：频率直方图绘制Matlab程序 

clc
clear
close all
fName = ["GKN.xlsx","GKF.xlsx","KBN.xlsx","KBF.xlsx"];
yLabels = ["高钾类未风化","高钾类风化","铅钡类未风化","铅钡类风化"];
xLabels = ["二氧化硅(SiO2)","氧化钾(K2O)","氧化钙(CaO)","氧化铝(Al2O3)","氧化铁(Fe2O3)","氧化铜(CuO)","五氧化二磷(P2O5)"];
for t=1:4
    A = xlsread(fName(t));
    [n,m] = size(A);
    
    idx = [2,4,5,7,8,9,12]
    for k=1:length(idx)
        for i = 1:n
            x(i)=A(i,idx(k));
        end
        
        picPos = k + (t-1)*7;
        subplot(4,7,picPos);
       
        
        h = histogram(x,5)
        if t ==2 || t == 4
            h.FaceColor = [0.1 0.1 0.5];
            h.EdgeColor = 'r';
        end
        hold on
        
        [counts,centers] = hist(x,5);
        x2 = centers(1)*0.5:((centers(end)-centers(1)))/1000:centers(end)*1.5;
        
        [mu,sigma]=normfit(x);%用正态分布拟合出平均值和标准差
        delta = centers(2)-centers(1);
        y2 = pdf('Normal', x2, mu,sigma)*10;%求在x2处的pdf值
        
        axis([centers(1)-delta,centers(end)+delta,0,max(counts)+0.5]);%限定x坐标范围
        
        hh = plot(x2,y2)
        set(hh,'LineWidth',5);
        
        if picPos==4
            title("不同玻璃类型风化前后不同化学成分频率分布直方图")
        end
        if t==4
           xlabel(xLabels(k)); 
        end
        if k==1
           ylabel(yLabels(t));
        end
    end
end

附录二： 加权平均预测Matlab程序

clc
clear
close all
fName = ["GKN.xlsx","GKF.xlsx","KBN.xlsx","KBF.xlsx"];
average = zeros(4,16);
for t=1:4
    A = xlsread(fName(t));
    [n,m] = size(A);
    for i = 1:m
       temp = sort(A(:,i));
       ave = 0;
       sum = 0;
       for j=1:n
           alpha = normpdf((j/n-0.5)*3);
           ave = ave + temp(j)*alpha;
           sum = sum + alpha;
       end
       ave = ave / sum;
       average(t,i)=ave;
    end
    
end
average

附录三：问题二灵敏度分析数据扰动Matlab程序

fName = ["高钾类风化","高钾类无风化","铅钡类风化","铅钡类无风化"];
val = [1,2,5,10,20,30,50,100];
for i=1:length(val)
    percent = val(i);
    dir = "rand" + percent;
    mkdir(dir);
    for t=1:4
        fname = fName(t) + ".xlsx";
        %target = fullfile
        copyfile(fname,dir);
        fname = dir + "/" + fname;
        A = xlsread(fname);
        [n,m] = size(A);
        R = rand(n,m-3).*0.01*percent + 1-percent*0.005;
        A(:,3:m-1) = A(:,3:m-1).*R;
        A(:,m) = sum(A(:,3:m-1),2);
        xlswrite(fname,A,"E2:U" + (n+1));
    end
end

附录四： 问题二亚类划分Matlab程序

clc
clear
close all
fName = ["高钾类风化.xlsx","高钾类无风化.xlsx","铅钡类风化.xlsx","铅钡类无风化.xlsx"];
dicName = [1,2,5,10,20,30];
T = xlsread("所有中心点.xlsx");
P = [[4,5],[2,3,10],[7,8,9,13],[11]];
B = zeros(4,2,13);
B(1,:,:) = T(1:1:2,2:1:14);%高钾类风化
B(2,:,:) = T(3:1:4,2:1:14);%高钾类无风化
B(3,:,:) = T(5:1:6,2:1:14);%铅钡类风化
B(4,:,:) = T(10:1:11,2:1:14);%铅钡类无风化
error = 0;
total = 0;
for t=1:4
    p = P(t);
    A = xlsread(fName(t));
    [n,m] = size(A);
    Data = A(:,4:m-1);
    for i=1:n      
       t1 = GetNorm(p,Data(i,:),B(t,1,:));
       t2 = GetNorm(p,Data(i,:),B(t,2,:));
       if t1>t2
           Result(t,i) = 1;
       else
           Result(t,i) = 2;
       end
    end
    for j = 1:length(dicName)
        D =  xlsread( "扰动后数据\rand" + dicName(j) + "\" + fName(t));
        [n,m] = size(D);
        Data = D(:,4:m-1);
        for ii=1:n
            t1 = GetNorm(p,Data(ii,:),B(t,1,:));
            t2 = GetNorm(p,Data(ii,:),B(t,2,:));
            if t1>t2
                res = 1;
            else
                res = 2;
            end
            if res ~= Result(t,ii)
                error = error + 1;
            end
            total = total + 1;
        end
        
        er(t,j) =(1- error/total)*100;
        total = 0;
        error = 0;
    end
end
er
 
function val = GetNorm(p,a,b)
    val = 0;
    for i=1:length(p)
        val = val + (a(1,p(i))-b(1,1,p(i)))*(a(1,p(i))-b(1,1,p(i)));
    end
%     for i=1:13
%         val = val + (a(1,i)-b(1,1,i))*(a(1,i)-b(1,1,i));
%     end
end