基于主成分分析法的远处燃烧体物质元素含量的测定模型

基于主成分分析法的远处燃烧体物质元素含量的测定模型---------------------------------------------------------------------------------------------------
基于主成分分析法的远处燃烧体物质元素含量的测定模型
2015大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料(包括网上查到的资料)，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示(包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等)。

我们参赛选择的题号是(从A/B中选择一项填写): B
我们的参赛报名号为(报名编号):所属学院(请填写完整的全名):
参赛队员 (打印并签名) :1.
日期: 2015 年 5 月 18 日
评阅编号(由组委会评阅前进行编号):
摘要:利用主成分分析法对远处燃烧体物质元素含量进行分析，———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
将得到的综合主成分指标代替原先较多的评价指标建立数学模型，从而对各物质元素含量测定进行评价分析。

结果表明，该方法能用降维的思想将问题简化，来得到可信客观的结论。

关键词:主成分分析法;碳元素含量;综合评价;指标
1
0 引言
对于燃烧体而言，若其碳元素含量不同，则各个光波长对应的光强度也不同。

因此，建立碳元素的含量与各个光波长对应的光强度之间的关系模型，对研究不同光强时的元素含量有着重要的意义。

然而在研究此问题的时候，总会借助于有一定相关性的多个指标，这样便增加了分析此问题的复杂性，而我们可以利用主成分分析法将这些多个相关的指标转换为少数互相无关的综合指标进行更方便的研究。

1 主成分分析法的原理
主成分分析法是利用一种降维的方法，用较少的新变量指标来代替众多的原变量，使这些新变量能尽可能的代表原变量所蕴含的信息，且彼此之间互不相关。

这种将多个变量用较少的互不相关的变量来表示的统计分析方法叫做主成分分析法。

定义:记X1,X2...XP为原变量指标，F1,F2...Fm为新变量指标，
?F1?a11X1?a12X2?...?a1pXp?F?aX?aX?...?aX?22112222pp ?......???
Fm?am1X1?am2X2?...?ampXp
(1) Fi与Fj互不相关;
(2) F1是X1,X2...XP的一切线性组合中方差最大者，F2是与F1不———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
相关X1,X2...XP的所有线性组合中方差最大者……。

所以新变量
F1,F2...Fm分别称为原变量指标X1,X2...XP的第1，第2，……第m主
成分。

2 主成分分析法的计算步骤
(1) 计算相关系数矩阵
?r11r12?r1p? ?r?r?r21222p ?R??????? ?? rr?r??p1p2pp??
rij(i，j=1，2，…，p)为原变量xi与xj的相关系数， rij?rji，
其计算公式为
rij??(xk?1nki?i)(xkj?j)2?(xk?1nki?i)?(xk?1nkj?j)22
(2) 计算特征值和特征向量
特征方程 ?I?R?0
用雅可比法求出特征值，并排序?1 ??2????p?0
而特征向量eij可用公式 e2?1?ij
(3) 计算主成分贡献率和累计贡献率
贡献率: ?
pk?1ikpj?1(i?1,2,?,p)?? 累积贡献率: ????k?1k?1pik(i?1,2,?,p)k
(4) 计算主成分载荷
lij?p(zi,xj)?ieij(i,j?1,2,?,p)
(5) 各主成分得分
?z11z12?z1m? ?z?z?z21222m Z???????? ??zz?z nm??n1n2
3 远处燃烧体物质含量的测定
问题(1):针对data0.csv包含的数据建立变量data由变量
———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
x1-x3648计算或者估计的数据模型:
A(对数据进行主成分分析:
本题所给的data0.csv数据中为第一行为格式数据，其中no，data，x1，
x2，…，x3648分别为数据的序号，碳元素的含量，245nm光强度，…，1044nm光强度。

此处光强度值是某一种固定参考强度的倍数，且x1-x3648变量的光波长差是一个常数。

在原数据中，每一次的样本包含3648个变量，要从19次的样本中直接得出碳含量与3648个变量的关系是不合实际的，故需要采用主成分分析法对3048个变量进行降维处理。

在实际的操作中，我们先使用matlab读入csv文件中的数据，再使用zscore 函数对数据矩阵进行标准化，然后使用matlab自带的主成分分析法的函数princomp得出从原始数据到各主成分的变换矩阵、降序排列的各成分的贡献率及变换后得出的各成分所构成的19*3548矩阵。

在对此数据的分析中我们得出第一主成分的贡献率达到了85.82%，故我们在之后的数据拟合中只采用19个样本的第一主成分作为自变量。

3
B(拟合
得到了十九个样本对应的第一主成分与碳含量之后，我们对这十九个点进行了
拟合。

有两种方案，
———————————————————————————————————————————————
--------------------------------------------------------------------
-------------------------------
a(多项式拟合
f(x) = p1*x^8 + p2*x^7 + p3*x^6 + p4*x^5 + p5*x^4 + p6*x^3 + p7*x^2
+ p8*x + p9
其中:
p1 = 6.55e-14-1.115e-13,
p2 = -3.935e-12
p3 = -1.212e-09
p4 = 1.507e-08
p5 = 4.738e-06
p6 = -9.252e-06
p7 = -0.004504
p8 = -0.001392
p9 = 21.69
b(fourier拟合
f(x) = a0 + a1*cos(x*w) + b1*sin(x*w) + a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w) + a4*cos(4*x*w) + b4*sin(4*x*w) +
a5*cos(5*x*w) + b5*sin(5*x*w) + a6*cos(6*x*w) + b6*sin(6*x*w) 其中:
a0 = 697.6 ，a1 = -1230 ，b1 = 292.1 ，a2 = 922.8 ，b2 =
-460.4 ，a3 = -564.6 ，b3 = 463.9 ，a4 = 270.5 ，b4 = -339.7 ，a5 = -88.85 ，b5 = 176.9 ，a6 = 12.92 ，b6 = -53.95 ，w = 0.02808。

问题(2):预测data1.csv包含的碳元素的含量data
———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
A(采用(1)中多项式拟合得到的模型，记x为第一主成分，则
4
data = p1*x^8 + p2*x^7 + p3*x^6 + p4*x^5 + p5*x^4 + p6*x^3 +
p7*x^2 + p8*x + p9
其中:
p1 = 6.55e-14-1.115e-13,，p2 = -3.935e-12 ，p3 = -1.212e-09 ，p4 =
1.507e-08 ，p5 = 4.738e-06 ，p6 = -9.252e-06 ，p7 = -0.004504 ，p8 = -0.001392 ，p9 = 21.69。

我们对data1.csv中的数据进行主成分分析时的变换得到另一组所需坐标下的数据，将第一列的数据代入上式进行验证:
Data =
21.3193411235572
21.0564220326128
20.9156753178776
21.0991331400088
19.4704661941218
21.2224487841091
20.5879370003697
20.8606104514896
21.0907562532152
20.4995893413615
———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
21.0348451838178
21.3032179740113
-181.859303329338
13.3428433132352
20.6158443047972
21.1401220072465
20.3668406319181
-9157.46356359878
-16350.4653108800
而data1.csv中
Data =
35.3900000000000
20.1000000000000
12.4800000000000
14.2700000000000
17.7200000000000
34.9000000000000
12.8900000000000
22.5100000000000
12.5800000000000
21.1400000000000
29.6200000000000 ———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
23.9100000000000
22.6000000000000
23.6800000000000
5
62.5300000000000
35.9500000000000
30.2200000000000
12.2500000000000
32.4500000000000
B(建立评估标准
记从模型中得出的data为d，所给的实际data为D。

则我们记S = sum(dij-Dij)2为评估模型准确性的参数。

对于通过主成分分析与多项式拟合得到的此模型，其在data1.csv中数据参考下的S为3.5253e+08;不难看出此模型在data1.csv中数据变换后的第一主成分下
三个点处存在较大的瑕疵，对S有着很大影响，需在后面进一步引入对信息的过滤和对模型的优化。

4 模型的修正
由于获得光谱数据时会受到许多不确定因素的影响，如样品不同成分之间相互干扰导致的光谱谱线的重叠、信号噪声等问题，因此可以对data0.csv所包含的数据进行消除噪声、波长优化选择等处理后可以更好的提高模型的性能。

1)噪声消除
———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
在实际问题中，通常使用小波变换进行降噪消噪。

小波变换通过伸缩平移运算对信号进行多尺度的变化，达到高频处时间细分，低频处频率细分，并能自动适应时频信号分析的要求，最终聚焦到信号的任意细节。

小波变换具有时频分析的优势，因此其在降噪消噪的同时也能很好地保留信号的突变部分;它也可以被看作低通滤波，很好地保留原始信号的特征;它还具有选基灵活性，针对不同的对象选取不同的小波基。

目前常用的小波变换降噪消噪的方法主要有空域相关消噪法、模极大值消噪法以及小波阈值消噪法。

其中小波阈值消噪法将绝对值较小的系数置零，绝对值较大的保留，最后重新构造出更理想的信号。

该方法简单、计算量小，因此在实际中得到了很大的应用。

在当前，小波变换也经常被用于图像去噪、地震数据处理等。

2)波长优化选择
波长优化选择是指对所获得的波长数据进行一定的筛选之后得到较为理想的数据，然后再参与数学建模之中。

与未经过波长优化选择而直接进行建模相比，其更
简单，且具有更高的稳定性。

目前常用到的波长优化选择的方法主要有遗传算法、逐步回归分析法以及区间偏最小二乘法等。

对于给出的数据，遗传算法将波长分布平均分为多个子空间，然后子空间进行不同的组合，最后找到一个能够反映全部数据的最佳组合来参与我们的数学建模。

经查阅资料以及相关的基本计算可知，遗传算法可以使波长点数减少，缩短建模所需时间，并且去除噪声比较———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
大的区域，从而提高模型的性能。

同时也可以在遗传算法中加入一种可以进行自我调节和交叉因子的值的算法，经过交叉因子形成正反馈，提升了遗传算法的运算速度。

这种方法在当前有很大的应用，如用于全光网中静态路由与波长分配的优化、选择近红外光谱谱区等。

5 燃烧物体中铁元素含量的研究
6
首先可以在网上查到铁元素的特征谱线，然后和燃烧物体的各元素的特征谱线进行比较，若燃烧物体的特征谱线不包含铁元素的特征谱线，则燃烧物体中铁元素含量为零;否则，燃烧物体含有铁元素。

接着可通过在网上查阅资料得到当燃烧物体含有不同含量的铁元素时，其各个波长对应的光强度。

参考(1)所建立的主成分分析模型，可建立基于主成分分析法的远处燃烧体物质铁元素含量的测定模型。

然后便可根据所建立的模型预测data0.csv、data1.csv包含的碳元素的含量data。

参考文献:
[1] 王凤花,朱海龙,戈振杨，近红外光谱数据建模方法的研究进展，56-60，2009。

[2] 田兵，主成分分析的数学模型及实际应用，18-21，2010。

[3] 李蔚,何军,刘德明,黄德修，改进遗传算法应用于全光网中静态路由与波长
分配的优化，133-136，2004。

[4]胡良剑,孙晓君，MATLAB数学实验，高等教育出版社，2008。

———————————————————————————————————————————————
---------------------------------------------------------------------------------------------------
附录:
mat文件中各变量的说明:
M: data0.csv中包含碳元素含量的数据;
X: data0.csv中不包含碳元素含量的数据;
X1:data0.csv中不包含碳元素含量的数据的标准化矩阵;
a:data0.csv中不包含碳元素含量的数据的标准化矩阵到另一坐标下的变换矩阵;
b:data0.csv中不包含碳元素含量的数据的标准化矩阵在另一坐标下的矩阵; c:b中各成分的贡献率;
data:根据已建立的模型从data1.csv中数据得出的19个样本的碳含量;
data1zhuchengfen:从data1.csv中数据得到的第一主成分的19个值; p1~9:模型中参数;
t:临时变量;
x:data0.csv中的第一主成分的19个值;
x1:data1.csv中的第一主成分的19个值;
y:data0.csv中的第一主成分的19个值对应的碳含量;
y1:data1.csv中的第一主成分的19个值对应的碳含量;
yz:data1.csv中不包含碳元素含量的数据;
yyz:data1.csv中不包含碳元素含量的数据的标准化矩阵;
yyyz:data1.csv中不包含碳元素含量的数据的标准化矩阵在另一坐标下的矩阵;
程序代码:
———————————————————————————————————————————————
--------------------------------------------------------------------
-------------------------------
M = csvread('data0.csv',1,1)
X = csvread('data0.csv',1,2);
X1 = zscore(X);
[a,b,c]=princomp(X1);
c=100.*c/sum(c)
x=b(:,1)
7
t=csvread('data0.csv',1,1)
t=t(:,1)
y=t
yz=csvread('data1.csv',1,2)
yyz=zscore(yz)
yyyz=yyz*a
x1=yyyz(:,1)
y1=csvread('data1',1,1)
y1=csvread('data1.csv',1,1)
y1=y1(:,1)
data1zhuchengfen=yyyz(:,1)
p1 = 6.55e-14-1.115e-13,
p2 = -3.935e-12
p3 = -1.212e-09
p4 = 1.507e-08
———————————————————————————————————————————————
--------------------------------------------------------------------
-------------------------------
p5 = 4.738e-06
p6 = -9.252e-06
p7 = -0.004504
p8 = -0.001392
p9 = 21.69
data = p1*data1zhuchengfen^8 + p2*data1zhuchengfen^7 +
p3*data1zhuchengfen^6 + p4*data1zhuchengfen^5 + p5*data1zhuchengfen^4 +
p6*data1zhuchengfen^3 + p7*data1zhuchengfen^2 + p8*data1zhuchengfen + p9 data = p1*data1zhuchengfen.^8 + p2*data1zhuchengfen.^7 +
p3*data1zhuchengfen.^6 + p4*data1zhuchengfen.^5 + p5*data1zhuchengfen.^4 + p6*data1zhuchengfen.^3 + p7*data1zhuchengfen.^2 + p8*data1zhuchengfen + p9
a=sum((y1-data).^2)
[a,b,c]=princomp(X1);
8
———————————————————————————————————————————————。