数学建模—大气污染预报问题

1、学生数学建模竞赛第一次预选赛(必做题)(1)油罐的体积(本题10分)一平放的椭圆柱体形状的油罐，长度为 L,椭圆的长半轴为a,短半轴为b,油的密度 为P，冋当油罐中油的咼度为h时油量是多少？解：由题意可话画出画出几何图形如图 1所示图1."x = acost椭圆方程为Jrbsint如图2,设阴影部分面积为S/2，则油桶的底面积为S。aX图下面将会利用mathematics5.0软件进行求解,求解的程序如下:In tegrate2*a*b*CostF2,t,ArcSi n1-h/b,Pi/2 解得结果为：S =a( J(2b ；h)h (b +h) +barccos(1 h)Vb2b当

2、h b时，由椭圆对称性，A中的h用2b-h代替得到： 所以油液质量M为：(2)光的反射定律(本题10分)费马原理：光总是沿用时最短的光程传播。试根据这一原理利用极值的有关知识证明光的反射定律：入射角等于反射角。解：由于光在同一介质中的速度为常数，所以在同一介质中光总是沿直线传播。如图3,现假设有两种介质1、2相接，光线在介质1中的传播速度为v,取两介质的分界线上的一条直线为X轴，设有一束光线从介质1中的A(O,a)点经X轴上的P(x,O)点反射，并沿直线方向行进到B(d,b)点。设直线AP与X轴法线的夹角为“，PB直线与X轴法线的夹角为匕，下面，根据最短时间效应来推导出光学中的反射定理P图三光

3、线由A点传到P点所需的时间为：光线由P点传到B点所需的时间为：故光线由A传到B所需的总时间为：根据费马定理，最短时间效应对应的优化问题为:dt 1xdx v a2 x2d x.b2 (d - x)2于是可以得到:xd _xa2x2. b2 (d -x)2又由于x2d x.b2 (d -x)2二 sinr2所以有：这就是光学中的反射定理。证毕大气污染预报问题摘要本文通过对四个城市的空气质量的排名以及城市 A的空气质量，利用C语言、Excel、 Mathematics和MATLAB!工具，分别建立了层次模型、多元线性回归预测模型 进行了合理 地分析。最后，我得到了以下一些比较满意的结果。问题（1）

4、:通过对问题（1）问题的分析，得出了这是一个比较典型的层次模型，目标层是空气质量的排名，因素是三种污染物的浓度情况，对象是题目给出的4个城市。查 找资料后，我找到了一个非常关键的东西一一 空气污染指数的计算方法，于是利用C语言 的编程知识我很快求出了这些城市的污染情况，结合层次模型的相关知识，建立层次模型 后很快得到了我们所期望的答案：总权重：A-0.1374,B-0.1301 ，C-0.1028，D-0.6298.于是城市环境排名也就解决了，由优到劣的排名情况：C B A B而且层次的模型的一次 性检验也顺利通过。问题（2）:问题（2）要求我们找出空气质量与气象因素之间的关系，一开始查阅了

5、很多资料，本想借助灰色预测模型进行求解，可是灰色预测模型的使用条件和咱们这个题 目的要求似乎没有什么关联，后来在网上浏览一片文章的时候，我找到了问题的突破口， 便是利用多元回归预测模型进行求解。然后根据这个模型的所要的处理数据，利用MATLAB Excel等工具，求解到本题的回归系数。得到了三种污染物与气象因子之间的关系：SO2的浓度与气象因素的关系：q =1.0173-0.0015p 0.0019 - 0.0004f - 0.0044/NO2的浓度与气象因素的关系:c2 二-0.2022 0.0003p 0.0017 0.0001f -0.0077/PM10的浓度与气象因素的关系：c3 =0

6、.9329-0.0014p + 0.0033-0.0022f +0.0019/接着，我又利用了 F检验和复相关系数R用来判别回归方程在统计上是否合理。结果还是 很让人满意的，回归模型的拟合度还是很高的。关键词：层次模型多元线性回归预测模型空气污染指数问题提出大气污染预报问题大气是指包围在地球外围的空气层，是地球自然环境的重要组成部分之一。人类生活在大气里，洁净大气是人类赖于生存的必要条件。一个人在五个星期内不吃饭或5天内不喝水，尚能维持生命，但超过5分钟不呼吸空气，便会死亡。随着地球上人口的急剧增加， 人类经济增长的急速增大，地球上的大气污染日趋严重，其影响也日趋深刻，如由于一些 有害气体的大

7、量排放，不仅造成局部地区大气的污染，而且影响到全球性的气候变化。因 此，加强大气质量的监测和预报是非常必要。目前对大气质量的监测主要是监测大气中SO2、N02、悬浮颗粒物（主要为PM10等的浓度。附件一给出了城市 A、B、C、D从2009年6月1日至2009年7月25日测量的污染物 含量及城市A的气象参数的数据；附件二给出了城市 A从2009年7月26日至2009年7 月30日测量的污染物含量及气象参数的数据。请解决下面两个问题：（1）建立由污染物浓度评价空气质量 的数学模型，然后利用附件一中的数据对四个 城市的空气质量进行排序。（2）分析城市A的空气质量（指SO2、NO2、PM10的浓度）与

8、气象参数之间的关系， 并利用附件二中的数据进行检验。二、基本假设1、题目所给的四个城市的污染物含量及城市 A的气象参数等数据都准确可靠。2、根据内地空气污染指数（API）来划分为5个等级：API值小于等于50,空气质量为优, 相当于国家空气质量一级标准；API值大于50且小于等于100,表明空气质量良好，相当 于达到国家质量二级标准；API值大于100且小于等于200,表明空气质量为轻度污染， 相当于国家空气质量三级标准；API值大于200表明空气质量差，称之为中度污染，为国 家空气质量四级标准；API大于300表明空气质量极差，已严重污染。三、符号说明符号意义备注入(max)对角矩阵的最大特

9、征值n矩阵的阶数CI层次模型的一致性指标CI=(入-n)/(n-1)RI随即一致性指标CR一次性比率P大气压强的数值单位：mmhgt温度的数值Cf空气湿度的数值v风速的数值m/sFF检验统计量R2预测模型的复相关系数四、问题分析1问题(1)的分析:要对A、B、C、D四个城市的空气质量进行排序，可从题目的要求中获知利用污染物 浓度来进行四个城市的排名。经过分析和查阅相关资料，这个问题应该属于典型的层次模 型的运用。下面将层次模型的相关内容说明如下：(1) 最大特征值入(max)的MATLA计算方法：V,D=eig(A),其中A为待计算特征值的矩阵, D为对角矩阵，其对角元素为A的特征值，最大的即

10、为入(max)。(2) 一致性指标CI计算方法：CI=(入-n)/(n-1);(其中入为矩阵A的最大特征值，n为矩阵的阶数)随即一致性指标RI的计算方法：RI与n有如下关系，如表n12345678910RI000.580.901.121.241.321.411.451.49(4)权重计算方法计算矩阵A的特征根及特征向量，将所求的特征向量单位化后得到的就是权重值。2、问题(2)的分析：通过仔细分析题目的要求，得知题目要求我们找出空气质量与气象因素的之间的的关 系。于是我首先想到了用预测模型去处理，但是由于变量太多而且，处理起来的拟合度太低了，达不到我们满意的要求。后来，通过仔细阅读相关资料找到了

11、以个比较好的模型一 多元线性回归预测模型，去处理五、模型的建立与求解5.1问题一模型建立与求解5.1.1问题一的分析要得到城市的空气质量的排名，首先我们得找一个量去权衡它们的关系才能得出比较 理想的结果，在前面的假设中我们便得到了，一个空气污染指数，我们可以以这个关键因 素作为突破口求解。5.1.2问题一模型的建立1、将研究目标（Z）、因素（P）、对象（C）按相关关系分成目标层Z、准则层P、对象层 Co层次结构图如图所示：城市城市A城市C城市B空气质量的排名2、给出空气质量一级，二级，三级两两成对比较的判断矩阵P污染级别一级二级三级123根据上图得出如下例两两成对比较的判断矩阵P一级二级三级权

12、重一级11/21/30.500二级212/30.333三级33/210.167由表中数据,计算可知：入（max）=3.00 , CI=0.00 , RI=0.58 , CR=0.00<0.1。因为 CR=0.00<0.1, 所以此排序有满意的一致性。3、给出对象层对准则层的各个因素的判断矩阵并进行分析。由于各个城市只存在污染程度的不同，所以它们两者之间各因素之间的关系。在这里我们利用了 C语言的相关知识求解出了这55天中各个城市的空气污染指数：四个城市的空气污染指数统计（单位：天）优良轻微污染轻度污染中度污染重污染A11430000B21340000C4870000D1637200

13、0根据表中数据，类比（2）中方法，计算出各种不同污染等级对不同城市的权重ABCD权重A111/2111/4811/160.115B21/11121/4821/160.219C48/1148/21148/160.500D16/1116/3916/4810.167由表中数据，计算可知：入（max）=4.00,CI=0.00,RI=0.90, CR=0.00<0.1所以此排序有满意的一致性。同理，可以计算出其余空气质量等级 4个城市的不同权重。计算方法类似，用 MATLAB 软件的计算过程详见本文附表。空气质量“良”级对4个城市的不同权重表城市ABCD权重0.3550.2810.0580.30

14、6由表中数据，计算可知：入（max）=4.00 , CI=0.00 , RI=0.90 , CR=0.00<0.1空气质量“轻微污染”级对3个城市的不同权重表城市ABCD权重0.0000.0000.0001.000由表中数据，计算可知：入（max）=1.00 , CI=0.00 , RI=0.00 , CR=0.00<0.1由于其他的污染指数均为零，在这里不再考虑了。5.1.3问题一模型的求解进行层次总排序，方法：将上面3个空气质量等级对4个城市的不同权重表单位化后作为列向量构成4X 3矩阵，和空气质量一级，二级，三级两两成对比较的判断矩阵P相乘，结果便得到4个城市的权重值。根据上

15、述问题的分析中的假设可知，权重值越大，表明空气污染情况越严重。 因此，将4个城市的权重值，按照从小到大依次排序，得出的结果便是4个城市的空气污染严重程度的排名。最终结果如下表所示：优（0.167）良（ 0.333）轻微污染（0.500）总权重A0.1150.3500.0000.1374B0.2190.2760.0000.1301C0.5000.0570.0000.1028D0.1670.3171.0000.6298根据4个城市的总权重值进行从小到大依次排序，空气污染严重程度的排名如下：C、B、A D5.1.4问题二结果的分析及验证总的一致性检验：CR=0.167X 0+0.333 X 0+0.

16、500 X 0=0<<0.1。此结果说明排序结有非常满意的一致性。结论显示城市C的空气质量状况最好，而D的空气质量状况最差。而且我们把 A和D做个 比较可以看出：虽然D的“优级”天数比A还多，但是由于D出现了两个轻度污染而造成D的 污染权重显着增加。5.2问题二模型建立与求解5.2.1问题二的分析我们首先利用Excel对SO2与各个气象因子之间的关系图如下面四个图所示：由于S02的浓度高低并不是由单一因素决定的，而是由于大气压强、地面风速、温度 以及湿度等气象因素共同影响的结果。因此，可以建立多元回归预测模型，对其浓度变化 规律进行分析和预测，从而实现对可吸入颗粒物(PMio)浓度

17、的最优控制522问题二模型的建立当前，对于大气污染物浓度预测所采取的方法主要是从污染物排放量高低为基础进行 预测的，典型的预测模型有：箱式模型、高斯扩散模式、多源扩散模式、线源扩散模式、 面源扩散模式和总悬浮微粒扩散模式。随着灰色系统、模糊数学和人工神经网络的发展， 预测方法又出现了以污染物排放相关因素为基础的模型，如：灰色预测模型(GM)、多元统计分析理论、模糊识别方法和人工神经网络预测方法。本题研究主要采用多元统计的方法 进行分析。在许多实际问题中，影响结果y的因素往往不止一个，而是多个变量 xi, X2 , xp与y之间存在着如下线性关系：y 0 ：必1 ：2%2 n “ HFpXip

18、；2 P P (1)I E ; = 0, D ( ;) = ；其中：® w '2, ，'p是回归系数；xi, x2, ，xp是p个可以精确测量或控制的变量，及回归因子；是不可观测的随机误差，满足E ； =0, D(；)=宀：住2未知)(2)一般地，我们称由公(1)和(2)确定的模型为多元线性回归模型，记为：jy =+久冷+內心+ %Xip +呂i2(3)E ；： =0, D(；)二；具体方法为：(1) 计算各变量的平均值：1 n - 1 n人&, yyt (4)n t壬n tv(2) 根据公式(5)计算出矩阵Lj和矩阵L：n Lij (xti - Xi )(

19、Xtj - X j )y(i,j =1,2j|l, p)(5)n_57 (禺-Xi)( yt y)t壬(3) 根据公式求出回归系数的估计值：即可求出回归模型：根据本题的特点，可以得到这样一个模型:523问题二模型的求解根据多元线性回归法的基本理论，分别考虑大气压强、温度、湿度和地面平均风速 个自变量，自变量分别以p、t、f、v表示，变量用ci表示，即武汉城区吸入颗粒物（PMio） 浓度，mg/m3。贝可设数学模型为：以环境空气质量自动监测子站监测的城区可吸入颗粒物 （PMio ）浓度数据，和相应的地 面平均风速、气温、相对湿度 3个气象因子为原始数据，先根据公式（4）利用Excel计算出 各变

20、量的平均值： 再按公式（5）利用MATLAB计算出Lj和Liy：23.15”-0.4879-6.211.2595L = §-1740， iy - 6.344710.78-0.0219890048694-195086.9446677-11249Lj = J -1950 -112491900123.15-6.21-仃 40最后根据公式（6）计算出回归系数的估计值：故根据多元线性回归方法，建立的城市 A的SO2的浓度拟合模型为:其中：&为SO2的预测浓度，mg/m； p为大气压强，mmhg t为地面温度，°C； f为 近地面空气中的湿度，% v为地面平均风速，m/s。利用

21、上面类似的方法可以求到：城市A的NO2的浓度的拟合模型为：其中：c2为NO2浓度，mg/m； p为大气压强，mmhg t为地面温度，C； f为近地面空气 中的湿度，% v为地面平均风速，m/s。城市A的PM10勺浓度的拟合模型为其中：C3为PM10测浓度，mg/m ； p为大气压强，mmhg t为地面温度，C； f为近地面空 气中的湿度，% v为地面平均风速，m/s。5.2.4问题二结果的分析及验证1、首先利用ExceI做出预测值与实际值之间的折线图：2、利用附表二中的数据结合 Excel表格进行检验:检验结果如下表所示:SO2NO2PM10预测值实际值预测值实际值预测值实际值0.031420

22、.0310.0319410.0370.0623430.0470.03110.0210.028690.0220.0653790.0300.017120.0250.0300370.0340.0468690.0340.0222360.0240.0328470.0350.0536870.0350.0125370.0260.0331380.0330.050110.081上面的见表格中：我们可以看出预测值和实际值之间还是比较吻合的。特别是S02与N02的预测值和实际值之间还是很吻合的，只有PM10的值稍差了一点。3、下面分别用F检验和复相关系数R用来判别回归方程在统计上是否合理。F检验统计量F的计算公式见

23、式（7）:UF = 一m（7）Qn - m -1其中，m为回归变量的自由度，n为观察值的组数，回归平方和 U和残差平方和Q的 计算公式见公式（8）:nU 八（？-y）27（8）Q 八（ - ?）2i 1复相关系数R的计算公式见式（9）:其中，回归平方和U和残差平方和Q的计算公式见公式（8）。1） S02函数关系的检验：选择所建预测模型的显着性水平为 0.05，而F检验的统计了 F=17.18>Fo.o5，预测模型 在统计意义上是显着成立的。预测模型的复相关系数 R2为0.9357，表明SO2浓度与气象因子（p、t、f、v）之间的 关系为高度正相关。预测模型的标准误差由相关表达式计算得0.

24、0139，因此，表明预测模型的拟合程度很2）各个污染物与气象参数之间关系式的检测情况表:三个个污染物的检验情况表相关指标污染物F0.05FR2标准误差SO20.0517.180.93570.0139NO20.0529.130.74260.0051PM100.0517.070.60150.9511由上面的表可以看出，我们建立的各个污染物与天气参数之间的关系式都是合理的六、模型的评价与推广6.1 模型的评价本文通过对大气污染预报问题的研究，建立了层次模型和多元线性回归预测模型，使 得问题得到了比较满意的解决，而且还得出三个污染物的预测方程，拟合度也满足要求。 但是通过后面五天的检验，我发现S02的

25、拟合度和让人满意，很多预测值和实际值差距很 小，但是PM10的拟合就不太好；所以模型仍然需要进一步的改进。6.2 模型的推广层次模型可以运用来解决我们日常生活中很多决策方面的问题，而且比较简单处理， 特别适合运用到政府部门对人口、交通、经济、环境等领域的发展规划做出决策。多元线性回归预测模型适合于类似与这种浓度预测中出现多因素的问题，可以使这类 问题得到很好的解决。七、参考文献1 姜启源等 , 数学模型（第三版） ，高等教育出版社， 2003年 8 月2 内空气污染指数计算方法：3 多元回归在武汉市城区可吸入颗粒物（PM10）浓度预测中的应用：4 数学建模： _城市空气质量评估及预测 （省级优

26、秀奖 ）八、附录8.1 附录清单 求解问题一的 求解问题一的 求解问题二的8.2 附录正文 附录 1： 求解问题一的C语言程序MATLA程序MATLA程序C 语言程序：#include<stdio.h> main()doubleS55,N55,M55,IS55,IN55,IM55,I55;doubleCS6=0.05,0.15,0.8,1.6,2.10,2.62,CN6=0.08,0.12,0.28,0.565,0.750,0.940, CM6=0.05,0.15,0.35,0.42,0.50,0.60;intA6=50,100,200,300,400,500,i;prin tf(

27、"输入 S02勺浓度 n");for(i=0;i<=54;i+)scanf("%lf",&Si);printf("输入 N02勺浓度 n");for(i=0;i<=54;i+)scanf("%lf",&Ni);printf(" 输入 PM1的浓度 n");for(i=0;i<=54;i+)scanf("%lf",&Mi);for(i=0;i<=54;i+)if(Si<=CS0)ISi=A0;elseif(Si<=CS

28、1)ISi=(A1-A0)/(CS1-CS0)*(Si-CS0)+A0;elseif(Si<=CS2)ISi=(A2-A1)/(CS2-CS1)*(Si-CS1)+A1;elseif(Si<=CS2)ISi=(A3-A2)/(CS3-CS2)*(Si-CS2)+A2;elseif(Si<=CS3)ISi=(A4-A3)/(CS4-CS3)*(Si-CS3)+A3;elseif(Si<=CS4)ISi=(A5-A4)/(CS5-CS4)*(Si-CS4)+A4;elseif(Si<=CS5)ISi=CS5+1;for(i=0;i<=54;i+)if(Ni<

29、;=CN0)INi=A0;elseif(Ni<=CN1)INi=(A1-A0)/(CN1-CN0)*(Ni-CN0)+A0;elseif(Ni<=CN2)INi=(A2-A1)/(CN2-CN1)*(Ni-CN1)+A1;elseif(Ni<=CN2)INi=(A3-A2)/(CN3-CN2)*(Ni-CN2)+A2;elseif(Ni<=CS3)INi=(A4-A3)/(CS4-CS3)*(Ni-CS3)+A3; elseif(Ni<=CS4)INi=(A5-A4)/(CS5-CS4)*(Ni-CS4)+A4;elseif(Ni>CN5)INi=CN5+1

30、;for(i=0;i<=54;i+)if(Mi<=CM0)IMi=A0;elseif(Mi<=CM1)IMi=(A1-A0)/(CM1-CM0)*(Mi-CM0)+A0;elseif(Mi<=CM2)IMi=(A2-A1)/(CM2-CM1)*(Mi-CM1)+A1;elseif(Mi<=CM2)IMi=(A3-A2)/(CM3-CM2)*(Mi-CM2)+A2;elseif(Mi<=CS3)IMi=(A4-A3)/(CS4-CS3)*(Mi-CS3)+A3;elseif(Mi<=CS4)IMi=(A5-A4)/(CS5-CS4)*(Mi-CS4)+A

31、4;elseif(Mi>CM3)IMi=CM5+1;printf("该城市的API分别为：n");for(i=0;i<=54;i+)Ii=ISi;if(Ii<INi)Ii=INi;if(Ii<IMi)Ii=IMi;printf("%3gt",Ii);printf("n");求解问题一的MATLAB?序A=111/2111/4911/14;21/1113/73/2;49/117/317/2;14/112/32/71A=1.00000.52380.22450.78571.90911.00000.42861.5000

32、4.45452.33331.00003.50001.27270.66670.28571.0000>>V,D=eig(A)V=-0.5137-0.1957-0.0040-0.07760.3269-0.37360.3620-0.57090.7627-0.8718-0.9317-0.52250.2179-0.24910.02990.6286D=-0.000000004.000000000.000000000.0000A=143/3443/643/40;34/43134/634/40;6/436/3416/40;40/4340/3440/61A=1.00001.26477.16671.07500.79071.00005.66670.85000.13950.17651.00000.15000.93021.17656.66671.0000>>V,D=eig(A)V=-0.92540.6312-0.3714-0.20190.24390.49910.7910-0.18360.04300.0881-0.0149-0.08320.28690.5872-0.48590.9584D=-0.