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

1、 学生数学建模竞赛第一次预选赛一、（必做题）（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分）费马原理：光

2、总是沿用时最短的光程传播。试根据这一原理利用极值的有关知识证明光的反射定律：入射角等于反射角。解：由于光在同一介质中的速度为常数，所以在同一介质中光总是沿直线传播。如图3，现假设有两种介质1、2相接，光线在介质1中的传播速度为v,取两介质的分界线上的一条直线为X轴，设有一束光线从介质1中的点经X轴上的点反射，并沿直线方向行进到点。设直线AP与X轴法线的夹角为，PB直线与X轴法线的夹角为，下面，根据最短时间效应来推导出光学中的反射定理。PYX0APB（d,b）图三光线由A点传到P点所需的时间为:光线由P点传到B点所需的时间为:故光线由A传到B所需的总时间为：根据费马定理，最短时间效应对应的优化问

3、题为： 令于是可以得到：又由于，所以有：这就是光学中的反射定理。 证毕大气污染预报问题摘要 本文通过对四个城市的空气质量的排名以及城市A的空气质量，利用C语言、Excel、Mathematics和MATLAB等工具，分别建立了层次模型、多元线性回归预测模型进行了合理地分析。最后，我得到了以下一些比较满意的结果。问题（1）：通过对问题（1）问题的分析，得出了这是一个比较典型的层次模型，目标层是空气质量的排名，因素是三种污染物的浓度情况，对象是题目给出的4个城市。查找资料后，我找到了一个非常关键的东西空气污染指数的计算方法，于是利用C语言的编程知识我很快求出了这些城市的污染情况，结合层次模型的相关

4、知识，建立层次模型后很快得到了我们所期望的答案：总权重城市环境排名也就解决了，由优到劣的排名情况：C、B、A、D。而且层次的模型的一次性检验也顺利通过。问题（2）：问题（2）要求我们找出空气质量与气象因素之间的关系，一开始查阅了很多资料，本想借助灰色预测模型进行求解，可是灰色预测模型的使用条件和咱们这个题目的要求似乎没有什么关联，后来在网上浏览一片文章的时候，我找到了问题的突破口，便是利用多元回归预测模型进行求解。然后根据这个模型的所要的处理数据，利用MATLAB、Excel等工具，求解到本题的回归系数。得到了三种污染物与气象因子之间的关系：SO2的浓度与气象因素的关系：NO2的浓度与气象因素

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

6、响到全球性的气候变化。因此，加强大气质量的监测和预报是非常必要。目前对大气质量的监测主要是监测大气中、悬浮颗粒物（主要为PM10）等的浓度。附件一给出了城市A、B、C、D从2009年6月1日至2009年7月25日测量的污染物含量及城市A的气象参数的数据；附件二给出了城市A从2009年7月26日至2009年7月30日测量的污染物含量及气象参数的数据。请解决下面两个问题：（1）建立由污染物浓度评价空气质量的数学模型，然后利用附件一中的数据对四个城市的空气质量进行排序。（2）分析城市A的空气质量(指、PM10的浓度)与气象参数之间的关系，并利用附件二中的数据进行检验。二、基本假设1、题目所给的四个城

7、市的污染物含量及城市A的气象参数等数据都准确可靠。2、根据内地空气污染指数（API）来划分为个等级：API值小于等于50，空气质量为优，相当于国家空气质量一级标准；API值大于50且小于等于100，表明空气质量良好，相当于达到国家质量二级标准；API值大于100且小于等于200，表明空气质量为轻度污染，相当于国家空气质量三级标准；API值大于200表明空气质量差，称之为中度污染，为国家空气质量四级标准；API大于300表明空气质量极差，已严重污染。三、符号说明符号意义备注(max)对角矩阵的最大特征值n矩阵的阶数CI层次模型的一致性指标CI=(-n)/(n-1)RI随即一致性指标CR一次性比率

8、CR=CIRIp大气压强的数值单位：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)；（其中为

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

10、题一模型建立与求解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二级

11、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根据表中数据，类比

12、（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

13、由表中数据，计算可知：(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个城市的不同权重表单位化后作为列向量构成4×3矩阵，和空气质量一级，二级，三级两两成对比较的判断矩阵P相乘，结果便得到

14、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&

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

16、度变化规律进行分析和预测，从而实现对可吸入颗粒物(PM10)浓度的最优控制5.2.2 问题二模型的建立当前，对于大气污染物浓度预测所采取的方法主要是从污染物排放量高低为基础进行预测的，典型的预测模型有：箱式模型、高斯扩散模式、多源扩散模式、线源扩散模式、面源扩散模式和总悬浮微粒扩散模式。随着灰色系统、模糊数学和人工神经网络的发展，预测方法又出现了以污染物排放相关因素为基础的模型，如：灰色预测模型(GM)、多元统计分析理论、模糊识别方法和人工神经网络预测方法。本题研究主要采用多元统计的方法进行分析。在许多实际问题中，影响结果y的因素往往不止一个，而是多个变量x1，x2，··&

17、#183;，xp与y之间存在着如下线性关系： (1)其中：···，是回归系数；x1，x2，···，xp是p个可以精确测量或控制的变量，及回归因子；是不可观测的随机误差，满足 (2)一般地，我们称由公(1)和(2)确定的模型为多元线性回归模型，记为： (3)具体方法为：(1)计算各变量的平均值： (4)(2)根据公式(5)计算出矩阵Lij和矩阵Li： () (5)(3)根据公式(6)求出回归系数的估计值： (6)即可求出回归模型：根据本题的特点，可以得到这样一个模型：5.2.3 问题二模型的求解根据多元线性回归法的基本理论，分别考虑大气

18、压强、温度、湿度和地面平均风速4个自变量，自变量分别以p、t、f、v表示，变量用表示，即武汉城区吸入颗粒物(PM10)浓度，mg/m3。则，可设数学模型为：以环境空气质量自动监测子站监测的城区可吸入颗粒物(PM10)浓度数据，和相应的地面平均风速、气温、相对湿度3个气象因子为原始数据，先根据公式(4)利用Excel计算出各变量的平均值：再按公式(5)利用MATLAB计算出Lij和Liy： ，最后根据公式(6)计算出回归系数的估计值：故根据多元线性回归方法，建立的城市A的SO2的浓度拟合模型为：其中：为SO2的预测浓度，mg/m3；p为大气压强，mmhg，t为地面温度，；f为近地面空气中的湿度，

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

20、031420.031 0.0319410.037 0.0623430.047 0.03110.021 0.028690.022 0.0653790.030 0.017120.025 0.0300370.034 0.0468690.034 0.0222360.024 0.0328470.035 0.0536870.035 0.0125370.026 0.0331380.033 0.050110.081 上面的见表格中：我们可以看出预测值和实际值之间还是比较吻合的。特别是SO2与NO2的预测值和实际值之间还是很吻合的，只有PM10的值稍差了一点。3、下面分别用F检验和复相关系数R用来判别回归方程在

21、统计上是否合理。F检验统计量F的计算公式见式(7): （） 其中，m为回归变量的自由度，n为观察值的组数，回归平方和U和残差平方和Q的计算公式见公式(8)： （8） 复相关系数R的计算公式见式(9)： （9）其中，回归平方和U和残差平方和Q的计算公式见公式(8)。1) SO2函数关系的检验：选择所建预测模型的显著性水平为0.05，而F检验的统计了F=17.18>F0.05，预测模型在统计意义上是显著成立的。预测模型的复相关系数R2为0.9357，表明SO2浓度与气象因子(p、t、f、v)之间的关系为高度正相关。预测模型的标准误差由相关表达式计算得0.0139，因此，表明预测模型的拟合程度

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

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

24、求解问题一的MATLAB程序求解问题二的MATLAB程序:8.2 附录正文附录1：求解问题一的C语言程序：#include<stdio.h>main()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;

25、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;

26、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<=C

27、N1)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(Ni<=CS4) INi=(A5-A4)/(CS5-CS4)*(Ni-CS4)+A4;else if(Ni>CN5)INi=CN5+1;for(i=0;i<=54;i+)i

28、f(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(Mi<=CS4) IMi=(A5-A4)/(CS5-CS4)*(Mi-CS4)+A4;else if(Mi>

29、;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=1 11/21 11/49 11/14;21/11 1 3/7 3/2;49/11 7/3 1 7/2;14/11 2/3 2/7 1A = 1.0000 0.5238 0.2245 0.7857 1.9091 1.0000 0.4286

30、 1.5000 4.4545 2.3333 1.0000 3.50001.2727 0.6667 0.2857 1.0000>> 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/4

31、0;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.0000>>B=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.2