中国各城市空气质量的聚类分析

中国各城市空气质量聚类和判别分析

摘要 中国经济的快速增长导致环境污染不断加重，其中空气污染与人类的健康密切相关。结合全国74个城市的空气污染物浓度数据对各城市进行聚类分析，分类方法包括等价关系法和Kmeans分析两种方法。结果表明，海口是全国空气质量最好的城市，石家庄和邢台是全国空气质量最差的城市，而武汉、成都和乌鲁木齐的空气质量与北京最为接近。

关键词 聚类分析 空气质量 等价关系 Kmeans

1. 介绍

随着中国经济的高速发展和工业化、城市化进程的加快，能源的消耗速度也不断提高。中国的工业发展大量依赖煤炭、石油等化石燃料，其燃烧产生的废气严重污染空气，导致中国各地区空气质量不断下降。90年代初期，中国的500个城市当中，达到国家空气质量I级标准的仅占1%；此外，近年来的数据显示，暴露于未达标空气中的城市人口占统计城市人口的三分之二[1]。城市的空气污染对人体健康构成极大威胁，研究表明，即便暴露于污染物密度较低的空气中也会提高慢性呼吸系统的发病率以及多种癌症的患病概率[2]。因此，有必要对全国各大城市的空气质量进行数据收集和分析，确定不同城市的污染程度及相互之间的关系，为相关部门制定政策提供有力的数据支撑。城市的空气污染程度主要受经济发展水平影响，但二者不是呈简单的倒U型曲线关系，不同的污染物与经济水平之间有不同的关系[3]，因此需要对各种污染物进行综合分析和评价。而目前对环境进行综合评价的方法包括模糊数学法、距离判别法和物元分析法[4]，本文即采取其中的模糊数学法对全国74个主要城市的空气污染数据进行聚类和判别分析，以研究目前中国各大城市的空气污染水平和特点。

2. 原始数据及聚类分析方法

本文所用到的城市空气污染数据来自环保部相关统计数据[5]，参见表1。为了便于分析，选取空气污染指标中量纲相同的三个指标进行考察，分别是SO2浓度、NO2浓度和PM2.5浓度。采用两种聚类分析方法对这74个城市进行分类，分别是等价关系法和Kmeans分类方法。

表1 全国主要城市空气污染物浓度

2.1 等价关系法

将每个城市的三种污染物浓度值作为其数据表征，根据最小最大法[6]建立74个城市之间的相似关系矩阵，得74×74矩阵，然后利用平方法求得相似矩阵的传递闭包，最后取不同的λ值将矩阵中的元素分为不同的类。相关程序代码见附录I，部分分类结果如表2所示。

表2 等价关系法部分分类结果

上表中，位于同一单元格内的城市为一类，位于同一行不同单元格内的城市属不同的类。从表中可看出，随着λ值从0开始增加，海口为第一个脱离大类的城市，这意味着海口的污染程度与其它城市的平均污染程度相差最远。由于污染物浓度值越大意味着污染程度越严重，这里定义城市的污染指数为3种污染物浓度值的代数平均，则得海口市的污染指数为17，而其它城市的平均污染指数为52.5，因此可认为海口是全国空气质量最好的城市。另一方面，当λ值增大至0.9附近时，石家庄和邢台依然分在同一类中，说明石家庄和邢台的空气污染程度极为接近，而邢台的污染指数为114，位于所有城市之首，因此石家庄和邢台是全国空气污染最严重的两个城市。当λ在0.87~0.91之间时，北京、武汉、成都、乌鲁木齐四个城市始终属于同一类，说明武汉、成都、乌鲁木齐的空气污染状况与北京比较接近。图1给出了几个典型城市的污染物浓度值与全国平均水平的比较情况，从图中依然可以看出，石家庄、唐山、邢台和邯郸的污染程度明显高于全国城市的平均水平，而北京的PM2.5指标值也高于全国平均水平，北京近年来的雾霾情况反映了这一点。

图1 几个典型城市的污染物浓度值

2.2. Kmeans分类方法

Kmeans算法[7]是一种基于质心簇的分类方法，首先在元素当中任意挑选k个元素作为初始质心，根据欧氏距离确定与每个元素相距最近的质心，并将每个元素与最近的质心归为一类，接下来再次计算每一簇的质心，再次重新归类，以此循环直至质心不再变化，算法可表示如下：

1) 初始化-确定聚类数和聚类中心

2) 计算所有样本与聚类中心的距离

3) 按照最小距离原则对各个样本进行聚类

4) 重新计算聚类中心

5) 重复1-4，直到聚类中心不再改变

由于初始质心的数量由用户自定义，所以最终分成的组数也是随意的，表3给出了划分为9组和5组的部分结果，相关程序代码参见附录Ⅱ。

表3 Kmeans法部分分类结果

从表中可看出，Kmeans方法分类结果中，北京、武汉、成都、乌鲁木齐同样属于同一类，类似地，石家庄和邢台也属同类，这一结果与等价关系法形成互相印证。

3. 结论

1) 等价关系法和kmeans算法能够有效地将样本分类

2) 聚类分析显示海口是空气最好的城市

3) 石家庄、邯郸、邢台和唐山空气质量较差的城市

4) 北京的空气质量同武汉、成都、乌鲁木齐相近

4. 人员分工

附录I 等价关系法程序代码

function r_sol= Cluster(lamda,r_equ)

% Written by Wensi Jen, 2016.03.22

% 求截矩阵

[m n]=size(r_equ);

r_tru=0*ones(m);

for i3=1:m

for j3=1:m

if r_equ(i3,j3)>=lamda

r_tru(i3,j3)=1;

else

r_tru(i3,j3)=0;

end

end

end

% 聚类

[m1,n1]=size(r_tru);

for i3=1:m1

var4=0;

for j3=1:n1

if(r_tru(i3,j3)==1)%找到每一行的元素1

var4=var4+1;

r_find(i3,var4)=j3;%把每一行找到的元素1的列下标顺次储存在矩阵的每一行中

end

end

end

[m2,n2]=size(r_find);

for i4=1:m2

for j4=2:n2

if (r_find(i4,j4)~=0)

r_find(r_find(i4,j4),:)=0;%将相同行只保留一个其他变为零行

end

end

end

var5=0;

for i5=1:m2

if (sum(r_find(i5,:)))~=0

var5=var5+1;

r_sol(var5,:)=r_find(i5,:);

end

end

end

[m n]=size(u);

% 数据标准化—极差变换

for i=1:m

for j=1:n

u_nor(i,j)=(u(i,j)-min(u(:,j)))/(max(u(:,j))-min(u(:,j)));

end

end

% 建立相似关系—最大最小法

for i1=1:m

for j1=1:m

var1=0;

var2=0;

for index=1:3

var1=var1+min(u_nor(i1,index),u_nor(j1,index));

var2=var2+max(u_nor(i1,index),u_nor(j1,index));

end

r(i1,j1)=var1/var2;

end

end

% 改造相似关系为等价关系—平方法求传递闭包

num=floor(log2(m))+1;

r_equ=r;

for k2=1:num

r=r_equ;

for i2=1:m

for j2=1:m

for index2=1:m

var3(index2,1)=min(r(i2,index2),r(index2,j2));

end

r_equ(i2,j2)=max(var3);

end

end

end

% 方法1—利用等价关系聚类

% 求lamda

% r_equ=[1 0.8 0.7;0.8 1 0.3;0.7 0.3 1]

[m_1,n_1]=size(r_equ);

tot=m_1*n_1;%矩阵r_equ中总的元素个数

r_row=r_equ(1:tot);%将矩阵r_equ转化为一个行向量

r_rows=sort(r_row);%将向量r_row中的元素按照从小到大排序

for i_1=1:tot

for j_1=i_1+1:tot

if (r_rows(i_1)==r_rows(j_1))

r_rows(j_1)=0;

end

end

end

r_rows1=sort(r_rows);%将r_rows中的元素按照从小到大排序

var_1=0;

for i_2=1:tot

if (r_rows1(i_2)==0)

var_1=var_1+1;%统计r_rows1中0的个数

end

end

r_rows2=r_rows1(tot:-1:var_1+1);%从r_rows2中挑出非零元素并按从大到小的顺序排列

s=length(r_rows2);

for i_3=1:s

lamda_row(i_3,1)=r_rows2(i_3); %lamda的列向量

end

% 直接聚类法

% 根据lamda_row中的lamda值和等价矩阵r_equ，调用函数Cluster

% 调用格式 Cluster(lamda,r_equ)

Result=Cluster(0.7,r_equ);

% 以城市名称表示结果

% Result=[1 3 0;2 0 0;7 0 0]

[m_2 n_2]=size(Result);

for i_4=1:m_2

for j_4=1:n_2

if Result(i_4,j_4)~=0

Result1{i_4,j_4}=cn{Result(i_4,j_4)};

end

end

end

附录Ⅱ Kmeans算法程序代码

K=9;%指定分为多少类

[Idx,C]=kmeans(u,K);%Idx N*1的向量，存储的是每个点的聚类标号

[m n]=size(u);

Result=0*ones(K,m);

[totr torc]=size(Idx);

for i=1:m

for j=1:m

if Result(Idx(i),j)==0

Result(Idx(i),j)=i;

break

end

end

end

% 以城市名称表示结果

% Result=[1 3 0;2 0 0;7 0 0]

[m_2 n_2]=size(Result);

for i_4=1:m_2

for j_4=1:n_2

if Result(i_4,j_4)~=0

Result1{i_4,j_4}=cn{Result(i_4,j_4)};

end

end

end

参考文献

[1] 张继娟，魏世强,我国城市大气污染现状与特点,四川环境,2006(03):104-108.

[2] 张莉，王五一，廖永丰,城市空气质量健康风险评估研究进展,地理科学进展,2006(03):39-47.

[3] 李茜，宋金平，张建辉 等,中国城市化对环境空气质量影响的演化规律研究,环境科学学报,2013(09):2402-2411.

[4] 鲁然英.城市环境空气质量及其评价方法研究:(硕士).兰州大学,2006.

[5] 2013中国主要城市的空气质量, .

[6] 杨纶标，高英仪,模糊数学原理及应用,华南理工大学出版社,2005.

[7] 基本Kmeans算法介绍及其实现, .

本文来源：https://www.2haoxitong.net/k/doc/6d191ebb49d7c1c708a1284ac850ad02df80070b.html