大气污染论文数学建模

大气污染评价与预报模型

摘要

本文对空气质量的评价及污染预报问题进行了分析，运用层次分析法依据处理后的数据对六个城市的空气质量进行了具体细致的排序；对2010年9月15日至9月21日的各项污染物浓度、各气象参数运用一元多项式回归模型进行了预测；就气象参数所属城市问题及污染物浓度与其的关系建立了相关性分析模型和多元线性回归模型；最后，根据建模过程和结果，我们对相关部门提出了几个具体的建议。

通过将数据附件所给有效数据，即日污染物浓度，转化为对应的月污染物浓度的均值，根据各城市月均污染浓度做出其随时间的走势折线图，分析了各个城市2SO 、2NO 、PM10之间的特点。我们拟根据API 指数值，以二级达标次数为准，对各城市之间的空气质量进行排名，但由于依据API 的区分空气质量等级时灵敏度较低，故采用了层次分析法对空气质量进行排名。由于我们采用了全部数据进行排名，而E 、F 数据较少，故只对ABCD 进行了排名。依据层次分析法得出的排名为：A 、B 、D 、C 。

为了精确预测各城市短期内的数据，本文选用一元多项式回归模型。对2010年的数据进行分析整理，依据回归模型得出其与时间的关系，得出预测值，并得出其置信度为95%的置信区间，结果显示模型的预测效果尚能接受，能够对所要预测数据进行预测。但由于F 城市数据缺失，根据假设做了合理的定性分析，并未对其进行定量预测。

分析空气质量与气象参数之间的关系时，首先根据数据完整性，气象参数应只属于其中一个城市，排除了D 、E 、F 的可能性，再根据相关性分析的方法，确定了气象参数属于A 城市。根据污染物与气象参数之间的因果关系，建立了多元线性回归模型，由于季节对污染物的浓度存在影响，分季节得出各污染物与各气象参数之间的相关系数，定性分析该相关系数，得出污染物与气象参数之间的关系。最后对该系数的理论与实际意义做了检验。

根据以上分析及结果，确定部分与空气质量控制相关的部门，针对其职能提出了诚恳建议。

关键词：API 评价模型 层次分析 一元多项式回归模型 相关性分析 多元回归

一、 问题重述

大气是指包围在地球外围的空气层，是地球自然环境的重要组成部分之一。人类生活在大气里，洁净大气是人类赖于生存的必要条件。随着地球上人口的急剧增加，人类经济增长的急速增大，地球上的大气污染日趋严重，其影响也日趋深刻。因此，加强大气质量的监测和预报是非常必要。

目前对大气质量的监测主要是监测大气中2SO 、2NO 、悬浮颗粒物（主要为PM10）等的浓度，研究表明，城市空气质量好坏与季节及气象条件的关系十分密切。

现有城市A 、B 、C 、D 、E 、F 从2003年3月1日至2010年9月14日测量的污染物含量及气象参数的数据。

1. 找出各个城市2SO 、2NO 、PM10之间的特点，并将几个城市的空气质量

进行排序。???

2. 对未来一周即2010年9月15日至9月21日各个城市的2SO 、2NO 、PM10

以及各气象参数作出预测。

3. 分析空气质量与季节、气象参数之间的关系。

4. 就空气质量的控制对相关部门提出你的建议。

二、 问题分析

问题一寻找各城市SO2、NO2、PM10的特点，最直观的方法就是作图，把各城市的三种污染物浓度做到一张图中进行比较分析可较容易的得到其特点， 而排序题目中给出的是三种污染物浓度，必须先用一个指标将它们统一起来综合的对城市的空气质量进行评价，用同一个指标进行排序。

问题二是依据所给的2010年1月1日至9月14日的数据，预测2010年9月15日至9月21日各个城市的SO2、NO2、PM10以及各气象参数，预测的时期较短，数目多，选择时间序列进行预测。同时将数据序列选取为2010年所有测量日的数据。根据给定的数据，利用一元多项式回归，求得回归模型，从而预测出需要的7项数据，并进行预测误差估计。

问题三是研究2SO 、2NO 、PM10的浓度与气象参数之间的关系。首先应对六个城市的2SO 、2NO 、PM10的浓度同气象参数进行相关性检验，以找出气象参数与2SO 、2NO 、PM10的浓度的对应关系。又由于城市空气质量好坏与季节及气象

条件的关系十分密切。故分季节对2SO 、2NO 、PM10的浓度与气象参数进行回归分析，并进行检验。

问题四则是通过对气象参数，季节等因素的考量，分析得出提高空气质量减小2SO 、2NO 、PM10的浓度的方法，从而给有关部门提出合理性的建议。

三、 模型假设

1. 各城市的检测时间具有随机性；

2. F 城市的发展是平衡发展，政府对环境治理干预较小，即F 城市的环境

不会出现强烈波动；

3. 数据附件所给的各个城市的污染物浓度及气象参数不随测量地点的改

变而改变，且有效数据都准确可靠。同时不考虑人为因素，检测仪器精确度不同的影响具有统计、预测意义；

4. 在对预测期内即2010年9月15日至9月21日时间段内，各个城市不

会出现重大空气污染事故，或环境不会变好；

5. 城市空气质量好坏与气象条件的关系十分密切，与其他因素关系不大。

四、 符号说明

五、 模型的建立和求解

(一)

找出各个城市2SO 、2NO 、PM10之间的特点，并将几个城市的空气质量

进行排序

数据处理：剔除1.第86项2005年的，与前后年份不符很明显是错误数据；2.第986项气压值的明显错误。

大气污染程度与空气中有害气体的含量有关，根据题目要求我们只考虑SO2、NO2、PM10的影响，由于数据量大我们按月为周期计算出了每个城市每种污染物浓度的月平均值

，据此画出折线图从而观察分析其特点。

对于排序问题，我们参考相关资料采用了判断大气污染的空气污染指数（API ），API 的计算依据为

其中，I 为某污染物的污染指数，C 为该污染物的污染浓度。C 大与C 小分别为上表中最贴近C 值的两个限值，C 大为大于C 的限值，C 小为小于C 的限值，同样，

I 大与I 小也是最靠近C 值的两个限值。

得出每种污染物的月平均污染指数。取三种污染物中API 中的最大值作为该市的月平均污染指数。 等级判断标准：

我们采用层次分析法对来处理，最终得到方案层对目标层权重，进行排序。

①建立以各个城市为方案层，空气质量为准侧层，空气质量排名为目标层的层次分析模型,由于E 、F 市数据过少，这里我们只考虑A 、B 、C 、D 市的排名，层次图如下

②根据两两比较法建立准则层对目标层的判断矩阵A

根据等级的1—9比较尺度定性的两两比较得到判断矩阵A 计算出权重向量E

②建立方案层对准则层的判断矩阵B

根据我们求得的各城市每月空气污染指数统计得

从而再次利用两两比较法分别得到各城市对不同空气等级的判断矩阵Bi ，根据特征根法确定权重向量F=（w1,w2,w3,w4,w5）通过一致性检验后，由

W=T

E F

得到最终方案层对目标层的权重，它表示的是受污染程度的权重，所以权重越大，受污染程度越高，空气质量越差。

求解

用Matlab 作图得到各城市直线图，下图中蓝色的实线为SO2浓度折线图，红色的虚线为PM10的浓度折线图，黑色的点线为NO2浓度折线图。

空气质量排名

Ⅰ

Ⅱ

Ⅴ

Ⅲ Ⅳ 城市D

城市A 城市B

城市C

从图中可看出A城市中SO2指数除前几个月有明显上升外，有明显的下降趋势，大部分低于PM10的指数，且变化与PM10有一定的一致性。PM10指数有较明显的振荡，但总体变化不是很大。NO2指数普遍最低，较平稳。

B城市的SO2、PM10变化很是相似，都以较大的幅度振荡且总趋势是下降的。SO2的振荡更明显，NO2指数均低于前两种污染物指数，较平稳。

与A、B一样C市中SO2、PM10的振荡具有一致性，但C市的PM10普遍比SO2高，且在20~25月左右有十分明显的上升，SO2振荡幅度不是很大，较平稳。

D市中SO2和PM10又出现大幅度的振荡总体趋于降低的趋势，NO2指数还是最低最平稳。

E 市中PM10指数高于另两种指数，在5月和22月左右各出现一次大幅度振荡，且下降趋势不明显，SO2指数在PM10下振荡下降，NO2出项小幅的振荡但总体平稳。

F 市数据太少，仅有的三个月数据看较平稳，但数据太少不足信。

1根据两两比较的到准则层对目标层的判断矩阵A A=11/21/31/41/5212/32/42/533/213/43/544/24/314/555/25/35/41????????????????

计算的max)（λ=，CI=0，CR1=0,一致性通过，所以得到权向量 E= (0.1238,0.1599,0.3857,0.3317)T

又根据表3得到方案层对准则层的判断矩阵B

B1=1

5/35/1153/513/11311/511/31111/51/31/111?????

???????

B2= 1

40/3240/234/331/40131/2331/3023/4023/31123/303/430/3130/231????????????

B3=115/2515/2115/1825/15125/2125/1821/1521/25121/1818/1518/2518/211?????

??????? B4=11/21/41/4211/21/2421142

11??

????????

?? 由于第五等级的数据较少，可以定性的判断其权重

运用Matlab 编程算出其权向量i w ，最大特征根i λ和一致性指标i CI ,结果列入下表

由表中i CI 值可知一致性检验全部通过。最后根据公式W=T

E F ，算得组合权向量W=(0.1238,0.1599,0.3857,0.3317)T

组合一致性检验

CR= 0*0.12380*0.15990.066*0.38570*0.3317+-+=<< 所以依权重，空气质量从优到劣的顺序为：A 、B 、D 、C

(二) 一元多项式回归分析预测2010年9月15日至9月21日的气象参数和污

染物浓度

环境中污染物的浓度、气象因子与时间之间存在着某种振荡关系。例如A 城市SO2的浓度随着时间呈现先减少后增加再减少的趋势。为了更准确的刻画这种污染物（气象因子）与时间的关系，我们建立了一元多项式回归模型。进而预测2010年9月15日至9月21日的污染物浓度和气象因子，最后进行预测误差估计，分析估计数据的置信度。

一般的，一元多项式回归分析的模型为：

012

0,n

n y x x E D βββε

εεσ

?=++++??==??L （**）

式中

2

01,,,,n βββσL 都是与

x

无关的未知参数，其中

201,,,,n βββσL 称为回归系数。本题中，由于观测数据为污染物浓度，

故（**）式中所示的

y 应是相应的污染物浓度。自变量元素为时间（日期

转化为对应的1，2，…），则具体模型为：

5

0152

0,y x x E D βββε

εεσ?=++++??==??

L （***） 用试验值（污染物浓度，气象参数），对201,,,,n βββσL 作点估计，

并对回归系数01,,,n βββL 作假设检验，最后在0x x =处对y 作预测，并

对

y 作区间估计。

利用

[],(,,,)Y DELTA polyconf p x s alpha =求polyfit 所得的

回归多项式在x 处的预测值y 及预测值的显着性为1-alpha 的置信区间

Y DELTA ±，从而能够分析出预测数据的准确性。

因为数据组数较多，仅以A 城市的22,,10SO NO PM 的浓度为例，分析

22,,10SO NO PM 浓度与时间的关系，估计出将来时间A 城市的22,,10SO NO PM 浓度，其他的数据和图像将在附录中给出。

首先读入A 城市2SO 数据，对读入数据进行拟合求解并分析。 以天为单位将日期转化成1,2,L ，运用一元多项式回归分析得到A 城市的2SO 随时间的变化曲线。A 城市2SO 浓度随时间变化曲线回归方程为： y =*x^5 + *x^4 +*x^3 + *x^2 + *x +

计算可得2010年9月15日至9月21日A城市

2

SO浓度为：

日期9/15 9/16 9/17 9/18 9/19 9/20 9/21 SO2浓度

然后进行预测的误差估计可得当alpha=时，求得的

2

SO的置信区间为：

日期9/15 9/16 9/17 9/18

SO2

浓度置信区间0.02190.0044

±0.02150.0044

±0.02110.0044

±0.02060.0045

±

日期9/19 9/20 9/21 SO2

浓度置信区间0.02010.0046

±0.01960.0046

±0.01910.0047

±

从以上表格可知

2

SO的置信区间较小，其估计值可信度比较高。

同时对拟合所得的曲线进行趋势分析，统计数据呈较大的振荡趋势，而通过

回归分析所得的A城市的

2

SO与时间的图线可得知，曲线虽然仍有小幅的振荡，但是该曲线仍旧呈整体下降趋势，这与预测所得的七天数据的趋势是相同的。

A市的

2,10

NO PM与时间的关系为：

NO浓度为：

计算可得2010年9月15日至9月21日A城市

2

日期9/15 9/16 9/17 9/18 9/19 9/20 9/21 NO2浓度

数据2010

年9月15日至9月21日的各个污染物浓度进行预测显然不合理，因此对影响F 城市空气质量的各污染物采取定性分析。由假设2， F 城市的环境不会出现强烈了变化，因此可认为F 城市空气中各个因素的浓度只在很小范围内波动，并且由第一个问题的综合评价模型得出，F 的排名在A 排名之前，可以对A 、F 进行对比作出定性的分析，分析如下：

在F 城市的环境不会出现强烈的变化的情况下，分析城市A 中各项污染物的API 的值，发现均有所下降，可以认为城市A 的空气质量有所改善，即城市A 的空气质量综合评价有所提高，缩小了与F 城市的差距。由此可认为F 城市各个污染物浓度变化走势与A 城市接近。

(三) 多元函数回归分析空气质量与气象参数之间的关系。

题目中所给定了六个城市的22,,10SO NO PM 浓度数据，但是只给定了一组气象参数，故在进行多元函数拟合前，必须进行相关性检验。

由概率论知识可知，两个随机变量ξ和η的相关性可以用相关系数(),ρξη来衡量

(),,Cov ξηρξη=

(),ρξη越大，

说明ξ和η之间的线性关系越显着，一般认为(),0.05ρξη≤可认为两者关系不大。经过求解，得到六个城市的22,,10SO NO PM 浓度的相关系数

可见气象参数与A 市的22,,10SO NO PM 浓度相关，故利用A 组数据与气象参数进行回归分析，得出22,,10SO NO PM 浓度与气象参数的关系。

(四) 就空气质量的控制对相关部门提出的建议

一方面，环保部门应对空气中污染物浓度数据的真实，有效，准确的检测并记录这对评估空气质量有着重要作用；另一方面，由第三个问题中的分析与求解可以得出，各个气象因素对每种污染物都有一定程度的影响，因此对气象因素的监测有利于环保部门做出对当前时期的主要污染物以及其未来一定时期的变动做出有效地评估分析；最后由第三问的求解结果分析，发现气体污染物的污染程度与风速，大气压等密切相关，因此，在城市规划中可以考虑采取人为措施，对环境进行改变。据此，我们对以下有关部门提出建议：

环保部门：应有效地做好对空气质量评价各个指标的监测，切实提高数据的真实性与可靠性，并根据各项污染物在一段时间内的变化趋势采取相应的治理措施。同时环保部门也应该做好监督的责任，按时查访各个工厂，检验废气的排放是否达标，工厂烟囱设计是否合理。对不达标，不合理的工厂要依法处理。

气象部门：根据当前气象参数对空气的污染程度，及时做出可行的人为调整，以提高空气质量。例如：当前和未来一定时期内为强烈干燥气候时，空气的相对湿度较低，部分气体污染物的污染相对增强，应采取人工降雨的措施增加空气相对湿度，以此降低空气中气体污染物，如SO2、NO2、PM10的浓度，以提高空气质量。

车辆监管部门及城市道路规划部门：在城市原有车辆的基础上，采取措施控制车辆的增长，以减少汽车尾气中硫化物、氮氧化物的排放，并对城市道路规划做出改善与完善，以避免车辆引起的二次扬尘增加空气中可吸入颗粒的浓度。

六、 模型的检验与评价

对于问题一一开始拟用现在通用的空气质量标准API ，后发现API 值不够精确，而采用了层次分析法，层次分析法是一种定性和定量相结合的、系统化、层次化的分析方法。可较全面的得到其结果，但由于数据问题，只对ABCD 四个城市进行了排序，若舍掉前四个城市钱30个月的数据就可以将E 城市加入排序。而如

果对模型进一步完善，将F 市也加入评价模型就更完美了。 对于问题二，利用了一元多项式回归分析得出六个城市了2SO 、2NO 、PM10的

浓度与事件的关系，预测出以后一周的2SO 、2NO 、PM10的浓度。虽然预测所得数据较符合变化规律，但是对数据进行了误差估计，置信度为95%的置信区间范围较大，数据不能令人满意。但是如果利用BP 网络进行预测，虽然函数关系较符合给定的数据，但是其存在较大的随机性，是预测出来的数值不准确，也不能称之为一个好方法。 对于问题三，首先进行相关性检验，分析解决了气象参数归属哪一个城市的问题，从总体上把握了气象参数与各城市三种污染物指标之间的相关关系，之后建立多元线性回归模型，定量的解释了各个气象参数与三种污染物浓度之间的关系，并能从实际意义与理论意义做出检验。由于采用多元线性分析，确定的是四个气象参数对每个污染物浓度的影响，并没有做剔除某一参数或控制某一参数对污染物浓度影响的分析，因此为了更好的说明气象参数与空气质量的关系，应该对此做出更为细致的分析。

参考文献：

[1] 赵静 但琦，数学建模与数学实验（第3版），北京：高等教育出版社，2003。 [2] 姜启源 谢金星 叶俊，数学模型（第三版），北京：高等教育出版，2003。 [3] 赵东方，数学模型与计算，北京：科学出版社，2007。

[4] 汤海波 肖培平 杨文增 郭新科,菏泽市气象因子与空气质量相关性研究与应用,中国环境监测，（5）.75-78。

[5] 王庆梅，大气污染预报技术及有关防治对策的研究，中国环境监测，(2).56-58。

[6] 匿名，空气污染指数，。

一、第一问的程序

求A市的API值,其他城市的类似

aapi=zeros(61,1);

apz=zeros(61,1);

for i=1:61

apz(i)=max(iso2(a(i,1)),ino2(a(i,2)));

aapi(i)=max(apz(i),ipm10(a(i,3)));

end

function I=ino2(t)

if 0<=t<

I=(50-0)/*(t-0);

else if t<

I =(100-50)/ else if t<

I =(200-100)/ else if t<

I=(300-200)/ else if t<

I=(400-300)/ else if t<

I=(500-400)/ end

end

end

end

end

end

end

function I=ipm10(t)

if t<

I=(50-0)/*(t-0);

else if t<

I =(100-50)/ else if t<

I =(200-100)/ else if t<

I=(300-200)/ else if t<

I=(400-300)/ else if t<=

I=(500-400)/ else if t> I=600;

end

end

end

end

end

end

end

function I=iso2(t)

if 0<=t<

I=(50-0)/*(t-0);

else if t<

I =(100-50)/ else if t<

I =(200-100)/ else if t<

I=(300-200)/ else if t<

I=(400-300)/ else if t< I=(500-400)/ end

end

end

end

end

end

end

统计A市的各个等级的数量，其他城市类似

agrade1=0;

agrade2=0;

agrade3=0;

agrade4=0;

agrade5=0;

agrade6=0;

for i=1:61

w=aapi(i);

if w<50

agrade1=agrade1+1;

else if w<100

agrade2=agrade2+1;

else if w<150

agrade3=agrade3+1;

else if w<200

agrade4=agrade4+1;

else if w<250

agrade5=agrade5+1;

else if w<300

agrade6=agrade6+1;

end

end

end

end

end

end

计算判断矩阵A，其他判断矩阵类似

A=[1 1/2 1/3 1/4 1/5;2 1 2/3 2/4 2/5;3 3/2 1 3/4 3/5;4 4/2 4/3 1 4/5;5 5/2 5/3 5/4 1];

[n,n]=size(A);

x=ones(n,100);

y=ones(n,100);

m=zeros(1,100);

m(1)=max(x(:,1));

y(:,1)=x(:,1);

x(:,2)=A*y(:,1);

m(2)=max(x(:,2));

y(:,2)=x(:,2)/m(2);

p=;i=2;k=abs(m(2)-m(1));

while k>p

i=i+1;

x(:,i)=A*y(:,i-1);

m(i)=max(x(:,i));

y(:,i)=x(:,i)/m(i);

k=abs(m(i)-m(i-1));

e nd

a=sum(y(:,i));

w=y(:,i)/a;

t=m(i);

disp('è¨?òá?:');disp(w);

disp('×?′óì??÷?μ:');disp(t);

CI=(λ-n)/(n-1);

RI=[0 0 ];

CR=CI/RI(n);

if CR<

disp('′????óμ?ò???D??éò??óêü!');

disp('CI=');disp(CI);

disp('CR=');disp(CR);

else

disp('′????óμ?ò???D?2??éò??óêü!');

end

二、第二问的求解数据与程序