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

1、中国各城市空气质量聚类和判别分析摘要 中国经济的快速增长导致环境污染不断加重，其中空气污染与人类的健康密切相关。结合全国74个城市的空气污染物浓度数据对各城市进行聚类分析，分类方法包括等价关系法和Kmeans分析两种方法。结果表明，海口是全国空气质量最好的城市，石家庄和邢台是全国空气质量最差的城市，而武汉、成都和乌鲁木齐的空气质量与北京最为接近。关键词 聚类分析 空气质量 等价关系 Kmeans1. 介绍随着中国经济的高速发展和工业化、城市化进程的加快，能源的消耗速度也不断提高。中国的工业发展大量依赖煤炭、石油等化石燃料，其燃烧产生的废气严重污染空气，导致中国各地区空气质量不断下降。90年代初

2、期，中国的500个城市当中，达到国家空气质量I级标准的仅占1%；此外，近年来的数据显示，暴露于未达标空气中的城市人口占统计城市人口的三分之二1。城市的空气污染对人体健康构成极大威胁，研究表明，即便暴露于污染物密度较低的空气中也会提高慢性呼吸系统的发病率以及多种癌症的患病概率2。因此，有必要对全国各大城市的空气质量进行数据收集和分析，确定不同城市的污染程度及相互之间的关系，为相关部门制定政策提供有力的数据支撑。城市的空气污染程度主要受经济发展水平影响，但二者不是呈简单的倒U型曲线关系，不同的污染物与经济水平之间有不同的关系3，因此需要对各种污染物进行综合分析和评价。而目前对环境进行综合评价的方法

3、包括模糊数学法、距离判别法和物元分析法4，本文即采取其中的模糊数学法对全国74个主要城市的空气污染数据进行聚类和判别分析，以研究目前中国各大城市的空气污染水平和特点。2. 原始数据及聚类分析方法本文所用到的城市空气污染数据来自环保部相关统计数据5，参见表1。为了便于分析，选取空气污染指标中量纲相同的三个指标进行考察，分别是SO2浓度、NO2浓度和PM2.5浓度。采用两种聚类分析方法对这74个城市进行分类，分别是等价关系法和Kmeans分类方法。表1 全国主要城市空气污染物浓度城市SO2浓度NO2浓度PM2.5浓度北 京265689天 津595496石 家 庄10568154唐 山1146911

4、5秦 皇 岛604765邯 郸9558139邢 台11369160保 定6956135张 家 口513240承 德373549沧 州5431102廊 坊4648110衡 水6846122太 原804381呼和浩特564057沈 阳904378大 连343152长 春444473哈 尔 滨445681上 海244862南 京375578无 锡404775徐 州524777常 州414872苏 州315370南 通283672连 云 港343567淮 安343279盐 城282965扬 州354370镇 江304272泰 州262677宿 迁333574杭 州285370宁 波224454温 州2

5、35158嘉 兴304768湖 州295274绍 兴384971金 华344170衢 州363768舟 山102233台 州173453丽 水193249合 肥223988福 州114336厦 门204436南 昌404069济 南9561110青 岛584367郑 州5952108武 汉336094长 沙334683广 州205253深 圳114040珠 海133738佛 山325352江 门273351肇 庆283854惠 州162938东 莞234547中 山194248南 宁193857海 口71727重 庆323870成 都316396贵 阳313353昆 明284042拉 萨922

6、26西 安4657105兰 州333567西 宁484170银 川774351乌鲁木齐2961882.1 等价关系法将每个城市的三种污染物浓度值作为其数据表征，根据最小最大法6建立74个城市之间的相似关系矩阵，得74×74矩阵，然后利用平方法求得相似矩阵的传递闭包，最后取不同的值将矩阵中的元素分为不同的类。相关程序代码见附录I，部分分类结果如表2所示。表2 等价关系法部分分类结果=0.0423所有城市为一类=0.4362海口北京，··· ，乌鲁木齐=0.87石家庄 唐山邢台=0.89石家庄 邢台=0.87北京武汉成都乌鲁木齐=0.89北京武汉成都乌鲁木齐

7、=0.9北京武汉成都乌鲁木齐=0.91北京武汉成都乌鲁木齐上表中，位于同一单元格内的城市为一类，位于同一行不同单元格内的城市属不同的类。从表中可看出，随着值从0开始增加，海口为第一个脱离大类的城市，这意味着海口的污染程度与其它城市的平均污染程度相差最远。由于污染物浓度值越大意味着污染程度越严重，这里定义城市的污染指数为3种污染物浓度值的代数平均，则得海口市的污染指数为17，而其它城市的平均污染指数为52.5，因此可认为海口是全国空气质量最好的城市。另一方面，当值增大至0.9附近时，石家庄和邢台依然分在同一类中，说明石家庄和邢台的空气污染程度极为接近，而邢台的污染指数为114，位于所有城市之首，

8、因此石家庄和邢台是全国空气污染最严重的两个城市。当在0.870.91之间时，北京、武汉、成都、乌鲁木齐四个城市始终属于同一类，说明武汉、成都、乌鲁木齐的空气污染状况与北京比较接近。图1给出了几个典型城市的污染物浓度值与全国平均水平的比较情况，从图中依然可以看出，石家庄、唐山、邢台和邯郸的污染程度明显高于全国城市的平均水平，而北京的PM2.5指标值也高于全国平均水平，北京近年来的雾霾情况反映了这一点。图1 几个典型城市的污染物浓度值2.2. Kmeans分类方法Kmeans算法7是一种基于质心簇的分类方法，首先在元素当中任意挑选k个元素作为初始质心，根据欧氏距离确定与每个元素相距最近的质心，并将

9、每个元素与最近的质心归为一类，接下来再次计算每一簇的质心，再次重新归类，以此循环直至质心不再变化，算法可表示如下：1) 初始化-确定聚类数和聚类中心2) 计算所有样本与聚类中心的距离3) 按照最小距离原则对各个样本进行聚类4) 重新计算聚类中心5) 重复1-4，直到聚类中心不再改变由于初始质心的数量由用户自定义，所以最终分成的组数也是随意的，表3给出了划分为9组和5组的部分结果，相关程序代码参见附录。表3 Kmeans法部分分类结果承德大连上海宁波温州台州丽水舟山福州厦门深圳珠海惠州海口石家庄唐山邯郸邢台天津沧州廊坊郑州西安北京哈尔滨南京无锡武汉成都乌鲁木齐秦皇岛张家口呼和浩特长春徐州青岛西宁

10、保定衡水济南南通连云港淮安盐城扬州镇江泰州太原沈阳银川北京天津沧州廊坊哈尔滨南京成都武汉乌鲁木齐石家庄唐山邯郸邢台保定衡水济南秦皇岛张家口太原呼和浩特沈阳青岛银川舟山台州丽水福州厦门深圳珠海江门惠州承德大连长春上海无锡常州苏州南通连云港从表中可看出，Kmeans方法分类结果中，北京、武汉、成都、乌鲁木齐同样属于同一类，类似地，石家庄和邢台也属同类，这一结果与等价关系法形成互相印证。3. 结论1) 等价关系法和kmeans算法能够有效地将样本分类2) 聚类分析显示海口是空气最好的城市3) 石家庄、邯郸、邢台和唐山空气质量较差的城市4) 北京的空气质量同武汉、成都、乌鲁木齐相近4. 人员分工数据的

11、收集和整理张文平、于鹏程序的编制任文希，彭炽，于鹏结果分析彭可文、张文平PPT制作彭炽，李冉，任文希报告编写李冉、彭可文附录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 endend% 聚类m1,n1=size(r_tru); f

12、or 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 endendvar5=0;for i5=1:m2 if (sum(r_find(i5,:)=0 var5=var

13、5+1; r_sol(var5,:)=r_find(i5,:); endendendm 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); endend% 建立相似关系最大最小法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,inde

14、x); end r(i1,j1)=var1/var2; endend% 改造相似关系为等价关系平方法求传递闭包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 endend% 方法1利用等价关系聚类% 求lamda% r_equ=1 0.8 0.7;0.8 1 0.3;0.7 0.3 1m_1,n_1=size(r_e

15、qu);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 endendr_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中

16、0的个数 endendr_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 0m_2 n_2=size(Res

17、ult);for i_4=1:m_2 for j_4=1:n_2 if Result(i_4,j_4)=0 Result1i_4,j_4=cnResult(i_4,j_4); end endend附录 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 endend% 以城市名称表示结果% Result=1 3 0;2 0 0;7 0 0m_2 n_2=size(Result);for i_4=1:m_2 for j_4=1:n_2 if Result(i_4,j_4)=0 Result1i_4,j_4=cnResult(i_4,j_4); end