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

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

摘要
近几年来环境问题成为全社会极为关注的热点, 空气污染是其中最热门的话题，同时也是最
重要的民生问题。本文针对这个现状，搜集了全国有代表性的31个城市的主要大气污染物的排放
情况， 先利用主成分分析评价了31个城市的综合空气质量，然后又分别用最短距离法和离差平方
和法进行聚类 分析，最终结果为北京、天津、石家庄等城市的空气质量较差；而海口、拉萨、南宁
等城市的空气较好。 特别需要说明的是北京的空气污染与其它城市相比有很大的不同，在最短距离
法中被单独聚为一类且与其 它类相距较远，这与北京目前空气现状是相吻合的。
在本文的最后还根据实际情况对模型的优缺点做了评价，并指出了需要改进的地方。
关键词：大气污染；主成分分析；聚类分析
1、数据资料
本文的原始数据取自《中国统计年鉴，2014》，
表1 我国主要城市废气中主要污染物排放情况
地区
北 京
天 津
石 家庄
太 原
呼和浩特
沈 阳
长 春
哈 尔滨
上 海
南 京
杭 州
合 肥
福 州
南 昌
济 南
郑 州
武 汉
长 沙
广 州
南 宁
海 口
重 庆
成 都
工业二氧
化硫
52041
207793
176469
88880
96190
130672
57246
65987
172867
110665
82021
41483
76043
40756
81118
106123
96222
21173
65589
33045
1798
494415
52040
工业氮氧
化物
75927
250646
200301
96018
131665
83348
95190
85515
262346
109693
67283
70311
72284
18597
72969
134120
95612
15951
57164
34797
86
247905
44411
工业烟
(粉)尘
27182
62766
99806
37003
48822
60425
72970
82323
67174
65256
40243
42387
43483
11413
47117
33828
20020
19545
16660
20950
1149
179842
21452
生活二氧
化硫
34967
8959
9564
33396
4257
14389
7344
50012
42947
1750
633
2710
1279
641
26087
11975
5720
2366
663
8748
11
53261
4891
生活氮氧
化物
13638
5221
2802
6738
665
5154
1545
22985
23474
400
335
130
169
58
3629
1780
1416
153
276
1068
17
4487
2109
生活烟尘
28258
18400
6635
26727
3763
15276
7919
80792
6451
1000
135
3188
547
254
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
15893
40109
52765
27170
52441
35493
5263
678
23831
7413
7129
5697
6691
1753
970
40
10951
1950
1419
1237
1425
5530
328
199
14012
1088
4793
3016
4920
用
x
1
表示工业 二氧化硫排放量，
x
2
表示工业二氧化硫排放量，
x
3
表示 工业烟(粉)尘排放量，
x
4
表示生活二氧化硫排放量，
x
5
表示生活氮氧化物排放量，
x
6
表示生活烟尘排放量。

2、主成分分析
2.1主成分分析的步骤
（1）计算相关系数矩阵
R?(r
ij
)
m?m
有
（2）计算特征值和特征向量。计算相关系数矩阵R的特征值
?
1
?
?2
????
?
m
?0
，以及对应的特
征向量
u
1
,u
2
,???u
m
由特征值组成m个新的指标变量：
其中：
y
1
是第一主成分，
y
2
是第二主成分，< br>（3）计算特征值的信息贡献率和累积贡献率。
为主成分
y
j
的信息贡献率，同时有
为主成分
y
1
,y
2
,,y
p
的累积贡献率。
，
y
m
是第m主成分。
（4）根据累积贡献率选取几个主成分作为新的评价指标。
2.2 主成分分析构建评价指标
定性地考虑反应各个城市空气质量的6个评价指标， 不难看出某些指标可能存在较强的相关
性，比如汽车的尾气中既含有二氧化硫也含有氮氧化物， 这两个 指标之间可能存在相关性。为了
验证这个想法用MATLAB计算指标之间的相关系数矩阵的特征值以及 贡献率，如下表所示：
表2 主成分分析结果
序号
1
2
3
4
5
6
特征值 贡献率% 累积贡献率%
2.30 38.27 38.27
1.57 26.20 64.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个主成分对应的特征向量

x
1

1
0.425989
2
-0.43086
3
-0.17122
x
2

0.429225
-0.31801
0.634299
x
3

0.447549
-0.32223
-0.42192
x
4

0.45322
0.269066
-0.22797
x
5

0.378817
0.471933
0.488375
x
6

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 )
中的最小元素，设为
D
KL
将
G
K
,G
L
合成一个新类记为
G
M
。
（3） 计算新类
G
M
与任一类之间的距离的递推公式为

在
D
(0)
中
G
K
,G
L
所在的行和列合并成一个 新行新列，对应
G
M
，该行列上的新距离由上述递
推公式求得，其余行列上的 距离值不变，得到新的距离矩阵 记作
D
(1)
。
（4）对
D< br>(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.021644>> 合并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(54)
第24步：d=0.15638>> 合并G(36)和G(54)记为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、模型分析

本模型的优点在于同时利用了主成分分析，最短距离法和离差平方和法对城市的空气质量进行< br>综合评价， 三种方法得到的结果相类似，都认为北京、天津等城市的空气较差， 也与这些城市的
实际空气情况相符，从而说明了模型的可靠性。
本模型的缺点是只考虑了城市 的污染物排放量，并没有考虑每个城市的面积，而且实际情况中
相邻近的城市还可能存在污染物的相互传 播问题，进一步的改进模型中可以用污染物排放量除以城
市的面积作为指标变量进行分析。
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=ze ros(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(sigma0N);
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(',num2 str(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)