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

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

2、进的地方。关键词：大气污染；主成分分析；聚类分析1、数据资料本文的原始数据取自中国统计年鉴，2014，表1 我国主要城市废气中主要污染物排放情况地区工业二氧化硫工业氮氧化物工业烟(粉)尘生活二氧化硫生活氮氧化物生活烟尘 北 京520417592727182349671363828258 天 津207793250646627668959522118400 石 家庄17646920030199806956428026635 太 原88880960183700333396673826727呼和浩特961901316654882242576653763 沈 阳130672833486042514389

3、515415276 长 春572469519072970734415457919 哈 尔滨659878551582323500122298580792 上 海1728672623466717442947234746451 南 京1106651096936525617504001000 杭 州820216728340243633335135 合 肥41483703114238727101303188 福 州7604372284434831279169547 南 昌40756185971141364158254 济 南8111872969471172608736298355 郑 州10612313

4、4120338281197517809150 武 汉962229561220020572014161001 长 沙21173159511954523661532946 广 州655895716416660663276214 南 宁330453479720950874810684631 海 口179886114911175 重 庆4944152479051798425326144874401 成 都52040444112145248912109661 贵 阳7060330450242333549317535530 昆 明10284268213573665263970328 拉 萨93020165

5、3867840199 西 安691033491715893238311095114012 兰 州721487991540109741319501088 西 宁718395328052765712914194793 银 川923698432127170569712373016乌鲁木齐7421611380352441669114254920用表示工业二氧化硫排放量，表示工业二氧化硫排放量，表示工业烟(粉)尘排放量，表示生活二氧化硫排放量，表示生活氮氧化物排放量，表示生活烟尘排放量。2、主成分分析2.1主成分分析的步骤（1）计算相关系数矩阵有（2）计算特征值和特征向量。计算相关系数矩阵R的特征值，以

6、及对应的特征向量由特征值组成m个新的指标变量：其中：是第一主成分，是第二主成分，是第m主成分。（3）计算特征值的信息贡献率和累积贡献率。为主成分的信息贡献率，同时有为主成分的累积贡献率。（4）根据累积贡献率选取几个主成分作为新的评价指标。2.2 主成分分析构建评价指标定性地考虑反应各个城市空气质量的6个评价指标， 不难看出某些指标可能存在较强的相关性，比如汽车的尾气中既含有二氧化硫也含有氮氧化物， 这两个指标之间可能存在相关性。为了验证这个想法用MATLAB计算指标之间的相关系数矩阵的特征值以及贡献率，如下表所示：表2 主成分分析结果序号特征值贡献率/%累积贡献率/%12.30 38.27 3

7、8.27 21.57 26.20 64.47 31.44 23.95 88.42 40.48 7.99 96.41 50.22 3.59 100.00 60.00 0.00 100.00 从结果中我们可以看出某些指标之间确实存在很强的相关性，比如生活烟尘对空气质量的贡献率接近0，这说明空气烟尘的值几乎可以由前面5个变量的值完全确定，这也是意料之中的结果。如果直接用这些指标进行综合评价，必然造成信息重叠，影响评价的客观性，因此可以考虑用主成分分析的方法来进行分析。从主成分分析结果中可以看出前3个主成分的累积贡献率就达到了88.42%，因此选取前3个主成分进行综合评价。前3个特征值对应的特征向量为

8、见表3。表3 前3个主成分对应的特征向量x1x2x3x4x5x610.4259890.4292250.4475490.453220.3788170.2916992-0.43086-0.31801-0.322230.2690660.4719330.5606093-0.171220.634299-0.42192-0.227970.488375-0.316于是得到3个主成分分别为分别以3个主成分的贡献率为权重，构建主成分综合评价模型，即把各地区的3个主成分值带入上式，可以得到各地空气质量的综合评价，见表4。表4 排名和综合评价结果2.3 结论从以上的综合评价结果可以看出。2013年全国各地区废气中污

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

10、为聚类指标，定义类与类之间的距离为两类最近样品间的距离，即（1） 计算n个样品的距离矩阵。（2） 选择中的最小元素，设为将合成一个新类记为。（3） 计算新类与任一类之间的距离的递推公式为在中所在的行和列合并成一个新行新列，对应，该行列上的新距离由上述递推公式求得，其余行列上的距离值不变，得到新的距离矩阵 记作。（4）对重复上述的2步得到，如此下去直到所有的元素合并成一类为止。3.2 最短距离法聚类模型(1) MATLAB的算法流程图如下：图1 程序流程图开始读取excel文件计算（31-）个样本的距离D（-1）找出D（-1）的最小值并返回其行row和列col将样本row和col合成一类=+1是

11、否聚类结束画聚类图(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(

12、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>> 合并

13、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>> 合并

14、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

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

16、原等城市的空气状况和北京类似。4、离差平方和法聚类分析直接调用MATLAB中的函数，得到聚类图如下图所示。从图中可以看出离差平方和法得到的结果仍把北京、天津、石家庄、呼和浩特归为一类，与最短距离法得到的结果相同。不同点在于离差平方和法使得两个大的类不容易合并，两个小的类容易合并，因而能达到分离的开的聚类结果，这符合我们对聚类的实际要求，因而更具有优越性。5、模型分析本模型的优点在于同时利用了主成分分析，最短距离法和离差平方和法对城市的空气质量进行综合评价， 三种方法得到的结果相类似，都认为北京、天津等城市的空气较差， 也与这些城市的实际空气情况相符，从而说明了模型的可靠性。本模型的缺点是只考虑

17、了城市的污染物排放量，并没有考虑每个城市的面积，而且实际情况中相邻近的城市还可能存在污染物的相互传播问题，进一步的改进模型中可以用污染物排放量除以城市的面积作为指标变量进行分析。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.m8、附录1、 主成分分析程序。clc,

18、clearx=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=ind2、

19、 最小距离法程序clc,clearx=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%计算s

20、igma，它是标准差的意思sigma0=zeros(M);for i=1:M for j=1:N sigma0(i)=sigma0(i)+(x0(i,j)-eq(i)2; endendsigma=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); endendtest=x0_jc0' %test为标准化后矩阵M,N=size(test);a='?'d_abs=zeros(M,M);d_ou0=zeros(M,M)

21、;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 endendtest=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 endendt=t+1;str='第',num2str(t),'步：d=',num2str(min),'>> 合并G(',. num2str(a(x),')和G(',num2str(a(y),')','记为