主要城市废气中主要污染物排放情况-(2013年)

摘要

近几年来环境问题成为全社会极为关注的热点, 空气污染是其中最热门的话题，同时也是最重要的民生问题。本文针对这个现状，搜集了全国有代表性的31个城市的主要大气污染物的排放情况，先利用主成分分析评价了31个城市的综合空气质量，然后又分别用最短距离法和离差平方和法进行聚类分析，最终结果为北京、天津、石家庄等城市的空气质量较差；而海口、拉萨、南宁等城市的空气较好。特别需要说明的是北京的空气污染与其它城市相比有很大的不同，在最短距离法中被单独聚为一类且与其它类相距较远，这与北京目前空气现状是相吻合的。

在本文的最后还根据实际情况对模型的优缺点做了评价，并指出了需要改进的地方。 关键词：大气污染；主成分分析；聚类分析

1、数据资料

本文的原始数据取自《中国统计年鉴，2014》，

表1 我国主要城市废气中主要污染物排放情况

地区 北 京 天 津 石 家庄 太 原 呼和浩特 沈 阳 长 春 哈 尔滨 上 海 南 京 杭 州 合 肥 福 州 南 昌 济 南 郑 州 武 汉 长 沙 广 州 南 宁 海 口 重 庆 成 都 工业二氧化硫 52041 207793 1769 88880 96190 130672 57246 65987 172867 110665 82021 41483 76043 40756 81118 106123 96222 21173 655 33045 1798 494415 52040 工业氮氧化物 75927 2506 200301 96018 131665 83348 95190 85515 262346 109693 67283 70311 72284 18597 72969 134120 95612 15951 571 34797 86 247905 44411 工业烟(粉)尘 27182 62766 99806 37003 48822 60425 72970 82323 67174 65256 40243 42387 43483 11413 47117 33828 20020 195 16660 20950 1149 179842 21452 生活二氧化硫 34967 59 95 33396 4257 143 7344 50012 42947 1750 633 2710 1279 1 26087 11975 5720 2366 663 8748 11 53261 41 生活氮氧化物 13638 5221 2802 6738 665 51 15 22985 23474 400 335 130 169 58 3629 1780 1416 153 276 1068 17 4487 2109 生活烟尘 28258 18400 6635 26727 3763 15276 7919 80792 51 1000 135 3188 7 2 8355 9150 1001 2946 214 4631 5 4401 661 贵 阳 昆 明 拉 萨 西 安 兰 州 西 宁 银 川 乌鲁木齐 70603 102842 930 69103 72148 71839 92369 74216 30450 68213 2016 34917 79915 53280 84321 113803 24233 57366 538 153 40109 52765 27170 52441 393 5263 678 23831 7413 7129 5697 6691 1753 970 40 10951 1950 1419 1237 1425 5530 328 199 14012 1088 4793 3016 4920 用x1表示工业二氧化硫排放量，x2表示工业二氧化硫排放量，x3表示工业烟(粉)尘排放量，

x4表示生活二氧化硫排放量，x5表示生活氮氧化物排放量，x6表示生活烟尘排放量。

2、主成分分析

2.1主成分分析的步骤

（1）计算相关系数矩阵R(rij)mm有

（2）计算特征值和特征向量。计算相关系数矩阵R的特征值12m0，以及对应的特征向量u1,u2,um由特征值组成m个新的指标变量：

其中：y1是第一主成分，y2是第二主成分，（3）计算特征值的信息贡献率和累积贡献率。 为主成分yj的信息贡献率，同时有 为主成分y1,y2,,yp的累积贡献率。

，ym是第m主成分。

（4）根据累积贡献率选取几个主成分作为新的评价指标。

2.2 主成分分析构建评价指标

定性地考虑反应各个城市空气质量的6个评价指标， 不难看出某些指标可能存在较强的相关性，比如汽车的尾气中既含有二氧化硫也含有氮氧化物， 这两个指标之间可能存在相关性。为了验证这个想法用MATLAB计算指标之间的相关系数矩阵的特征值以及贡献率，如下表所示：

表2 主成分分析结果

序号 1 2 3 4 5 6 特征值 贡献率/% 累积贡献率/% 2.30 38.27 38.27 1.57 26.20 .47 1.44 23.95 88.42 0.48 7.99 96.41 0.22 3.59 100.00 0.00 0.00 100.00 从结果中我们可以看出某些指标之间确实存在很强的相关性，比如生活烟尘对空气质量的贡献率接近0，这说明空气烟尘的值几乎可以由前面5个变量的值完全确定，这也是意料之中的结果。如果直接用这些指标进行综合评价，必然造成信息重叠，影响评价的客观性，因此可以考虑用主成分分析的方法来进行分析。

从主成分分析结果中可以看出前3个主成分的累积贡献率就达到了88.42%，因此选取前3个主成分进行综合评价。

前3个特征值对应的特征向量为见表3。

表3 前3个主成分对应的特征向量

x1 1 0.4259 2 -0.43086 3 -0.17122 x2 0.429225 -0.31801 0.634299 x3 0.4479 -0.32223 -0.42192 x4 0.45322 0.269066 -0.22797 x5 0.378817 0.471933 0.488375 x6 0.291699 0.560609 -0.316 于是得到3个主成分分别为

分别以3个主成分的贡献率为权重，构建主成分综合评价模型，即

把各地区的3个主成分值带入上式，可以得到各地空气质量的综合评价，见表4。

表4 排名和综合评价结果

2.3 结论

从以上的综合评价结果可以看出。2013年全国各地区废气中污染物情况存在较大的差异，重庆作为我国的国家中心城市是经济中心、金融中心和创新中心，同时也是全球著名的6大雾都城市年平均雾日是104天，在这次主成分的得分中位居第一，上海作为国际著名的经济中心，2014年GDP总量居中国城市第一、亚洲第二，工商业活动频繁，空气质量较差。北京作为我国的首都城市空气质量不容乐观，在31个城市中排名23；海口是海南省省会，空气质量排名第一，也在2011年世界卫生组织(WHO)发布首份全球城市空气污染调查报告，获中国空气最清洁的城市，其次拉萨、长沙、南宁等城市的空气质量也一直位居全国前列。

3、最短距离法聚类分析

3.1 最短距离法的步骤

选取主成分分析中的前3个主成分得分作为聚类指标，定义类与类之间的距离为两类最近样品间的距离，即

（1） 计算n个样品的距离矩阵D(0)。

（2） 选择D(0)中的最小元素，设为DKL将GK,GL合成一个新类记为GM。 （3） 计算新类GM与任一类之间的距离的递推公式为

在D(0)中GK,GL所在的行和列合并成一个新行新列，对应GM，该行列上的新距离由上述递

推公式求得，其余行列上的距离值不变，得到新的距离矩阵 记作D(1)。

（4）对D(1)重复上述D(0)的2步得到D(2)，如此下去直到所有的元素合并成一类为止。

3.2 最短距离法聚类模型

(1) MATLAB的算法流程图如下：

图1 程序流程图

开始 读取excel文件 计算（31-i）个样本的距离D（i-1） 找出D（i-1）的最小值 将样本row和col合成一类 是 i=i+1 否 聚类结束 (2)计算结果 画聚类图 第1步：d=0.013451>> 合并G(30)和G(31)记为G(32)

第2步：d=0.0214>> 合并G(24)和G(25)记为G(33) 第3步：d=0.023779>> 合并G(22)和G(23)记为G(34) 第4步：d=0.048656>> 合并G(13)和G(15)记为G(35) 第5步：d=0.051912>> 合并G(8)和G(10)记为G(36) 第6步：d=0.059749>> 合并G(28)和G(29)记为G(37) 第7步：d=0.060741>> 合并G(17)和G(19)记为G(38) 第8步：d=0.067867>> 合并G(34)和G(38)记为G(39) 第9步：d=0.070466>> 合并G(27)和G(37)记为G(40) 第10步：d=0.073204>> 合并G(14)和G(39)记为G(41) 第11步：d=0.095139>> 合并G(7)和G(16)记为G(42) 第12步：d=0.10012>> 合并G(35)和G(41)记为G(43) 第13步：d=0.10023>> 合并G(18)和G(40)记为G(44) 第14步：d=0.10111>> 合并G(6)和G(43)记为G(45) 第15步：d=0.10566>> 合并G(33)和G(45)记为G(46) 第16步：d=0.10831>> 合并G(26)和G(44)记为G(47) 第17步：d=0.11077>> 合并G(12)和G(42)记为G(48) 第18步：d=0.11078>> 合并G(46)和G(47)记为G(49) 第19步：d=0.11621>> 合并G(20)和G(49)记为G(50)

得到聚类图如下图所示第20步：d=0.11995>> 合并G(21)和G(48)记为G(51) 第21步：d=0.12293>> 合并G(50)和G(51)记为G(52) 第22步：d=0.13024>> 合并G(32)和G(52)记为G(53) 第23步：d=0.15135>> 合并G(9)和G(53)记为G() 第24步：d=0.15638>> 合并G(36)和G()记为G(55) 第25步：d=0.18692>> 合并G(11)和G(55)记为G(56) 第26步：d=0.22048>> 合并G(2)和G(3)记为G(57) 第27步：d=0.30904>> 合并G(4)和G(56)记为G(58) 第28步：d=0.40013>> 合并G(5)和G(57)记为G(59) 第29步：d=0.45985>> 合并G(58)和G(59)记为G(60) 第30步：d=1.5322>> 合并G(1)和G(60)记为G(61)

图2 最短距离法聚类图

表5 城市编号

（4） 结论

从图中可以看出，如果根据空气质量把31个地区分为3类结果为： 第一类 北京 第二类 天津 石家庄 呼和浩特 第三类 其他城市 如果分为4类结果为：

第一类 第二类 第三类 第四类 北京 呼和浩特 天津 石家庄 其他城市 从聚类结果可以看出，北京的空气质量与其他城市相比有很大的不同，可能的原因是生活烟尘的排放量过高，这与北京近年来的空气质量状况是相符的，其次石家庄、天津、太原等城市的空气状况和北京类似。

4、离差平方和法聚类分析

直接调用MATLAB中的函数，得到聚类图如下图所示。

从图中可以看出离差平方和法得到的结果仍把北京、天津、石家庄、呼和浩特归为一类，与最短距离法得到的结果相同。不同点在于离差平方和法使得两个大的类不容易合并，两个小的类容易合并，因而能达到分离的开的聚类结果，这符合我们对聚类的实际要求，因而更具有优越性。

5、模型分析

本模型的优点在于同时利用了主成分分析，最短距离法和离差平方和法对城市的空气质量进行综合评价， 三种方法得到的结果相类似，都认为北京、天津等城市的空气较差， 也与这些城市的实际空气情况相符，从而说明了模型的可靠性。

本模型的缺点是只考虑了城市的污染物排放量，并没有考虑每个城市的面积，而且实际情况中相邻近的城市还可能存在污染物的相互传播问题，进一步的改进模型中可以用污染物排放量除以城市的面积作为指标变量进行分析。

6、参考文献

[1] 王学民．应用多元分析[M]．上海财经大学出版社，2014：153-197. [2] 司守奎王．数学建模算法与应用[M]．国防工业出版社，2014：193-207.

7、附件

(1) 主要城市废气中主要污染物排放情况 (2013年).xls (2) principal_component_analysis.m (3) Single_linkage_method.m (4) Ward_method.m

8、附录

1、主成分分析程序。 clc,clear

x=xlsread('主要城市废气中主要污染物排放情况 (2013年)',2); x0=x';

r=corrcoef(x0');%求相关系数矩阵

[vec1,lamda,rate]=pcacov(r)%求相关系数矩阵的特征值以及特征向量 f=repmat(sign(sum(vec1)),size(vec1,1),1); num=3;

df=x0'*vec1(:,1:num);

tf=df*rate(1:num)/100;%计算综合得分

[stf,ind]=sort(tf,'descend');%把得分按高到低排序 stf=stf',ind=ind 2、最小距离法程序 clc,clear

x=xlsread('主要城市废气中主要污染物排放情况 (2013年)',4,'B2:D32'); x0=x';

[M,N]=size(x0);m=zeros(1,M);n=9999*ones(1,M);s=zeros(1,M);eq=zeros(1,M); for i=1:M for j=1:N

if x0(i,j)>=m(i) m(i)=x0(i,j); end

if x0(i,j)<=n(i) n(i)=x0(i,j);

end

s(i)=s(i)+x0(i,j); end

eq(i)=s(i)/N; end

%计算sigma，它是标准差的意思 sigma0=zeros(M); for i=1:M for j=1:N

sigma0(i)=sigma0(i)+(x0(i,j)-eq(i))^2; end end

sigma=sqrt(sigma0/N); jicha=m-n; he=sum(x0,2);

x0_jc0=zeros(M,N); for i=1:M for j=1:N

x0_jc0(i,j)=x0(i,j)/jicha(i); end end

test=x0_jc0'; %test为标准化后矩阵 [M,N]=size(test); a='?';

d_abs=zeros(M,M);d_ou0=zeros(M,M); for i=1:M for j=1:M for k=1:N

d_abs(i,j)=d_abs(i,j)+abs(test(i,k)-test(j,k)); end end end

test=d_abs;t=0;

M=length(test(1,:));MM=M;a=1:MM; Z=zeros(MM-1,3);

disp('最短距离聚类分析结果：') while(sum(sum(test))) min=9999;

for i=1:M %在test中找出最大的相关系数及其下标... for j=1:M %...并提示：合并下标对应的两组数据 if(min>test(i,j)&&test(i,j)~=0) min=test(i,j); x=i;y=j; end end end t=t+1;

str=['第',num2str(t),'步：d=',num2str(min),'>> 合并G(',...

num2str(a(x)),')和G(',num2str(a(y)),')','记为G(',num2str(t+MM),')']; disp(str) %提示操作步骤

Z(t,:)=[a(x),a(y),min]; %收集dendrogram()画聚类图时所需的数据

a([x,y])=[]; %每执行一步，在a中删除被合并的数据号，在末尾顺次新增一个 if(~isempty(a)) a(end+1)=t+MM; end

g=zeros(1,M-2); %生成新的距离矩阵 ii=0; for i=1:M

if(i==x||i==y); else

ii=ii+1;

g(ii)=(test(x,i)=test(y,i))*test(y,i); end %...两列变一列：两两比较取其小者 end

test([x,y],:)=[]; test(:,[x,y])=[]; test=cat(1,test,g);

test=cat(2,test,[g';0]); M=length(test(1,:)); end

dendrogram(Z,0); %画树状聚类图 title('最短距离聚类'); ylabel('lambda');

clear M MM i ii j x y g a t min str Z 3、离差平方和法聚类程序 clc,clear

x=xlsread('主要城市废气中主要污染物排放情况 (2013年)',4,'B2:D32'); x0=x;

y=pdist(x0,'cityblock'); yc=squareform(y) Z=linkage(y,'ward') [h,t]=dendrogram(Z)