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

学生数学建模竞赛第一次预选赛一、（必做题）（1）油罐的体积（本题10分）一平放的椭圆柱体形状的油罐，长度为L，椭圆的长半轴为a，短半轴为b，油的密度为，问当油罐中油的高度为h时油量是多少？解：由题意可话画出画出几何图形如图1所示2aL2bh图 1.1椭圆方程为如图2，设阴影部分面积为S/2，则油桶的底面积为S。bXYab-h图 2下面将会利用mathematics 5.0软件进行求解,求解的程序如下：Integrate2*a*b*Cost2,t,ArcSin1-h/b,Pi/2解得结果为： 当时，由椭圆对称性，A中的h用代替得到：所以油液质量M为：（2）光的反射定律（本题10分）费马原理：光总是沿用时最短的光程传播。试根据这一原理利用极值的有关知识证明光的反射定律：入射角等于反射角。解：由于光在同一介质中的速度为常数，所以在同一介质中光总是沿直线传播。如图3，现假设有两种介质1、2相接，光线在介质1中的传播速度为v,取两介质的分界线上的一条直线为X轴，设有一束光线从介质1中的点经X轴上的点反射，并沿直线方向行进到点。设直线AP与X轴法线的夹角为，PB直线与X轴法线的夹角为，下面，根据最短时间效应来推导出光学中的反射定理。PYX0APB（d,b）图三光线由A点传到P点所需的时间为:光线由P点传到B点所需的时间为:故光线由A传到B所需的总时间为：根据费马定理，最短时间效应对应的优化问题为： 令于是可以得到：又由于，所以有：这就是光学中的反射定理。 证毕大气污染预报问题摘要 本文通过对四个城市的空气质量的排名以及城市A的空气质量，利用C语言、Excel、Mathematics和MATLAB等工具，分别建立了层次模型、多元线性回归预测模型进行了合理地分析。最后，我得到了以下一些比较满意的结果。问题（1）：通过对问题（1）问题的分析，得出了这是一个比较典型的层次模型，目标层是空气质量的排名，因素是三种污染物的浓度情况，对象是题目给出的4个城市。查找资料后，我找到了一个非常关键的东西空气污染指数的计算方法，于是利用C语言的编程知识我很快求出了这些城市的污染情况，结合层次模型的相关知识，建立层次模型后很快得到了我们所期望的答案：总权重：A-0.1374,B-0.1301，C-0.1028，D-0.6298.于是城市环境排名也就解决了，由优到劣的排名情况：C、B、A、D。而且层次的模型的一次性检验也顺利通过。问题（2）：问题（2）要求我们找出空气质量与气象因素之间的关系，一开始查阅了很多资料，本想借助灰色预测模型进行求解，可是灰色预测模型的使用条件和咱们这个题目的要求似乎没有什么关联，后来在网上浏览一片文章的时候，我找到了问题的突破口，便是利用多元回归预测模型进行求解。然后根据这个模型的所要的处理数据，利用MATLAB、Excel等工具，求解到本题的回归系数。得到了三种污染物与气象因子之间的关系：SO2的浓度与气象因素的关系：NO2的浓度与气象因素的关系：PM10的浓度与气象因素的关系：接着，我又利用了F检验和复相关系数R用来判别回归方程在统计上是否合理。结果还是很让人满意的，回归模型的拟合度还是很高的。关键词：层次模型 多元线性回归预测模型 空气污染指数问题提出大气污染预报问题大气是指包围在地球外围的空气层，是地球自然环境的重要组成部分之一。人类生活在大气里，洁净大气是人类赖于生存的必要条件。一个人在五个星期内不吃饭或5天内不喝水，尚能维持生命，但超过5分钟不呼吸空气，便会死亡。随着地球上人口的急剧增加，人类经济增长的急速增大，地球上的大气污染日趋严重，其影响也日趋深刻，如由于一些有害气体的大量排放，不仅造成局部地区大气的污染，而且影响到全球性的气候变化。因此，加强大气质量的监测和预报是非常必要。目前对大气质量的监测主要是监测大气中、悬浮颗粒物（主要为PM10）等的浓度。附件一给出了城市A、B、C、D从2009年6月1日至2009年7月25日测量的污染物含量及城市A的气象参数的数据；附件二给出了城市A从2009年7月26日至2009年7月30日测量的污染物含量及气象参数的数据。请解决下面两个问题：（1）建立由污染物浓度评价空气质量的数学模型，然后利用附件一中的数据对四个城市的空气质量进行排序。（2）分析城市A的空气质量(指、PM10的浓度)与气象参数之间的关系，并利用附件二中的数据进行检验。二、基本假设1、题目所给的四个城市的污染物含量及城市A的气象参数等数据都准确可靠。2、根据内地空气污染指数（API）来划分为个等级：API值小于等于50，空气质量为优，相当于国家空气质量一级标准；API值大于50且小于等于100，表明空气质量良好，相当于达到国家质量二级标准；API值大于100且小于等于200，表明空气质量为轻度污染，相当于国家空气质量三级标准；API值大于200表明空气质量差，称之为中度污染，为国家空气质量四级标准；API大于300表明空气质量极差，已严重污染。三、符号说明符号意义备注(max)对角矩阵的最大特征值n矩阵的阶数CI层次模型的一致性指标CI=(-n)/(n-1)RI随即一致性指标CR一次性比率p大气压强的数值单位：mmhgt温度的数值f空气湿度的数值v风速的数值m/sFF检验统计量R2预测模型的复相关系数四、问题分析1、问题（1）的分析：要对A、B、C、D四个城市的空气质量进行排序，可从题目的要求中获知利用污染物浓度来进行四个城市的排名。经过分析和查阅相关资料，这个问题应该属于典型的层次模型的运用。下面将层次模型的相关内容说明如下： (1) 最大特征值(max) 的MATLAB计算方法：V,D=eig(A),其中A为待计算特征值的矩阵，D 为对角矩阵，其对角元素为A的特征值，最大的即为(max) 。(2)一致性指标CI 计算方法： CI=(-n)/(n-1)；（其中为矩阵A的最大特征值，n为矩阵的阶数）(3)随即一致性指标RI 的计算方法：RI 与n 有如下关系，如表n12345678910RI000.580.901.121.241.321.411.451.49(4)权重计算方法计算矩阵A 的特征根及特征向量，将所求的特征向量单位化后得到的就是权重值。2、问题（2）的分析：通过仔细分析题目的要求，得知题目要求我们找出空气质量与气象因素的之间的的关系。于是我首先想到了用预测模型去处理，但是由于变量太多而且，处理起来的拟合度太低了，达不到我们满意的要求。后来，通过仔细阅读相关资料找到了以个比较好的模型多元线性回归预测模型，去处理。五、模型的建立与求解5.1 问题一模型建立与求解5.1.1 问题一的分析要得到城市的空气质量的排名，首先我们得找一个量去权衡它们的关系才能得出比较理想的结果，在前面的假设中我们便得到了，一个空气污染指数，我们可以以这个关键因素作为突破口求解。5.1.2 问题一模型的建立1、 将研究目标（Z）、因素（P）、对象（C）按相关关系分成目标层Z、准则层P、对象层C。层次结构图如图所示：城市A城市B城市C良轻微污染轻度污染空气质量的排名城市D 优 中度污染 重 污染2、给出空气质量一级，二级，三级两两成对比较的判断矩阵P污染级别一级二级三级123根据上图得出如下例两两成对比较的判断矩阵P一级二级三级权重一级1/21/30.500二级22/30.333三级33/20.167由表中数据, 计算可知：(max) = 3.00 ，CI = 0.00 ，RI = 0.58 ，CR = 0.00 0.1。因为CR = 0.00 0.1，所以此排序有满意的一致性。3、给出对象层对准则层的各个因素的判断矩阵并进行分析。由于各个城市只存在污染程度的不同，所以它们两者之间各因素之间的关系。在这里我们利用了C语言的相关知识求解出了这55天中各个城市的空气污染指数：四个城市的空气污染指数统计（单位：天）优良轻微污染轻度污染中度污染重污染A11430000B21340000C4870000D 16372000根据表中数据，类比（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.306由表中数据，计算可知：(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 由于其他的污染指数均为零，在这里不再考虑了。 问题一模型的求解进行层次总排序，方法：将上面3个空气质量等级对4个城市的不同权重表单位化后作为列向量构成43矩阵，和空气质量一级，二级，三级两两成对比较的判断矩阵P相乘，结果便得到4个城市的权重值。根据上述问题的分析中的假设可知，权重值越大，表明空气污染情况越严重。因此，将4个城市的权重值,按照从小到大依次排序，得出的结果便是4个城市的空气污染严重程度的排名。最终结果如下表所示:优（0.167）良（0.333）轻微污染（0.500）总权重A0.115 0.3500.0000.1374 B0.2190.2760.0000.1301C0.5000.0570.0000.1028D0.1670.3171.0000.6298根据4个城市的总权重值进行从小到大依次排序，空气污染严重程度的排名如下：C、B、A、D5.1.4 问题二结果的分析及验证总的一致性检验：CR = 0.167 0 + 0.333 0 +0.500 0 = 0 F0.05，预测模型在统计意义上是显著成立的。预测模型的复相关系数R2为0.9357，表明SO2浓度与气象因子(p、t、f、v)之间的关系为高度正相关。预测模型的标准误差由相关表达式计算得0.0139，因此，表明预测模型的拟合程度很高。2) 各个污染物与气象参数之间关系式的检测情况表：三个个污染物的检验情况表 相关指标污染物F0.05FR2标准误差SO20.0517.180.93570.0139NO20.0529.130.74260.0051PM100.0517.070.60150.9511由上面的表可以看出，我们建立的各个污染物与天气参数之间的关系式都是合理的。六、模型的评价与推广6.1 模型的评价本文通过对大气污染预报问题的研究，建立了层次模型和多元线性回归预测模型，使得问题得到了比较满意的解决，而且还得出三个污染物的预测方程，拟合度也满足要求。但是通过后面五天的检验，我发现SO的拟合度和让人满意，很多预测值和实际值差距很小，但是PM的拟合就不太好；所以模型仍然需要进一步的改进。6.2 模型的推广层次模型可以运用来解决我们日常生活中很多决策方面的问题，而且比较简单处理，特别适合运用到政府部门对人口、交通、经济、环境等领域的发展规划做出决策。多元线性回归预测模型适合于类似与这种浓度预测中出现多因素的问题，可以使这类问题得到很好的解决。七、参考文献1 姜启源等, 数学模型（第三版），高等教育出版社，2003年8月2 内空气污染指数计算方法： 3 多元回归在武汉市城区可吸入颗粒物(PM10)浓度预测中的应用： 4 数学建模：_城市空气质量评估及预测(省级优秀奖) 八、附录8.1 附录清单求解问题一的C语言程序求解问题一的MATLAB程序求解问题二的MATLAB程序:8.2 附录正文附录1：求解问题一的C语言程序：#includemain()double S55,N55,M55,IS55,IN55,IM55,I55;double CS6=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;int A6=50,100,200,300,400,500,i;printf(输入SO2的浓度n);for(i=0;i=54;i+)scanf(%lf,&Si);printf(输入NO2的浓度n);for(i=0;i=54;i+)scanf(%lf,&Ni);printf(输入PM10的浓度n);for(i=0;i=54;i+)scanf(%lf,&Mi);for(i=0;i=54;i+)if(Si=CS0)ISi=A0;else if(Si=CS1)ISi=(A1-A0)/(CS1-CS0)*(Si-CS0)+A0;else if(Si=CS2)ISi=(A2-A1)/(CS2-CS1)*(Si-CS1)+A1;else if(Si=CS2)ISi=(A3-A2)/(CS3-CS2)*(Si-CS2)+A2;else if(Si=CS3) ISi=(A4-A3)/(CS4-CS3)*(Si-CS3)+A3;else if(Si=CS4) ISi=(A5-A4)/(CS5-CS4)*(Si-CS4)+A4;else if(Si=CS5)ISi=CS5+1;for(i=0;i=54;i+)if(Ni=CN0)INi=A0;else if(Ni=CN1)INi=(A1-A0)/(CN1-CN0)*(Ni-CN0)+A0;else if(Ni=CN2)INi=(A2-A1)/(CN2-CN1)*(Ni-CN1)+A1;else if(Ni=CN2)INi=(A3-A2)/(CN3-CN2)*(Ni-CN2)+A2;else if(Ni=CS3) INi=(A4-A3)/(CS4-CS3)*(Ni-CS3)+A3;else if(NiCN5)INi=CN5+1;for(i=0;i=54;i+)if(Mi=CM0)IMi=A0;else if(Mi=CM1)IMi=(A1-A0)/(CM1-CM0)*(Mi-CM0)+A0;else if(Mi=CM2)IMi=(A2-A1)/(CM2-CM1)*(Mi-CM1)+A1;else if(Mi=CM2)IMi=(A3-A2)/(CM3-CM2)*(Mi-CM2)+A2;else if(Mi=CS3) IMi=(A4-A3)/(CS4-CS3)*(Mi-CS3)+A3;else if(MiCM3)IMi=CM5+1;printf(该城市的API分别为：n)；for(i=0;i=54;i+)Ii=ISi;if(IiINi)Ii=INi;if(Ii V,D=eig(A)V = -0.5137 -0.1957 -0.0040 -0.0776 0.3269 -0.3736 0.3620 -0.5709 0.7627 -0.8718 -0.9317 -0.5225 0.2179 -0.2491 0.0299 0.6286D = -0.0000 0 0 0 0 4.0000 0 0 0 0 0.0000 0 0 0 0 0.0000A=1 43/34 43/6 43/40;34/43 1 34/6 34/40;6/43 6/34 1 6/40;40/43 40/34 40/6 1A = 1.0000 1.2647 7.1667 1.0750 0.7907 1.0000 5.6667 0.8500 0.1395 0.1765 1.0000 0.1500 0.9302 1.1765 6.6667 1.0000 V,D=eig(A)V = -0.9254 0.6312 -0.3714 -0.2019 0.2439 0.4991 0.7910 -0.1836 0.0430 0.0881 -0.0149 -0.0832 0.2869 0.5872 -0.4859 0.9584D = -0.0000 0 0 0 0 4.0000 0 0 0 0 0.0000 0 0 0 0 -0.0000B=1 1/2 1/3;2 1 2/3;3 3/2 1B = 1.0000 0.5000 0.3333 2.0000 1.0000 0.6667 3.0000 1.5000 1.0000 A=0.116 0.352 0;0.221 0.279 0;0.516 0.049 0;0.147 0.320 1A = 0.1160 0.3520 0 0.2210 0.2790 0 0.5160 0.0490 0 0.1470 0.3200 1.0000 A*Bans = 0.8200 0.4100 0.2733 0.7790 0.3895 0.2597 0.6140 0.3070 0.2047 3.7870 1.8935 1.2623 0.8200 0.4100 0.2733 0.7790 0.3895 0.2597 0.6140 0.3070 0.20473.7870 1.8935 1.2623 B=0.167;0.333;0.500B = 0.1670 0.3330 0.5000 A=0.115 0.355 0;0.219 0.281 0;0.500 0.058 0;0.167 0.306 1A = 0.1150 0.3550 0 0.2190 0.2810 0 0.5000 0.0580 0 0.1670 0.3060 1.0000 A*Bans = 0.1374 0.1301 0.1028 0.6298附录2：求解问题二的MATLAB程序:A=890.042386.9409-1.95E+0323.146286.9409466.7687-1124.9-6.2106-1.95E+03-1124.919001-173.994823.1462-6.2106-173.994810.7827 %矩阵LijA = 1.0e+004 * 0.0890 0.0087 -0.1950 0.0023 0.0087 0.0467 -0.1125 -0.0006 -0.1950 -0.1125 1.9001 -0.0174 0.0023 -0.0006 -0.0174 0.0011B=inv(A) %求A的逆矩阵B = 0.0015 0.0001 0.0001 -0.0007 0.0001 0.0027 0.0002 0.0048 0.0001 0.0002 0.0001 0.0013 -0.0007 0.0048 0.0013 0.1179B*C %求回归系数ans = -0.0015 0.0003 -0.0014 0.0019 0.0017 0.0033 -0.0004 0.0001 -0.0002 -0.0044 -0.0077 0.0019