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

我国主要城市废气中主要污染物排放情况摘要近几年来环境问题成为全社会极为关注的热点,空气污染是其中最热门的话题，同时也是最重要的民生问题。

本文针对这个现状，搜集了全国有代表性的31个城市的主要大气污染物的排放情况，先利用主成分分析评价了31个城市的综合空气质量，然后又分别用最短距离法和离差平方和法进行聚类分析，最终结果为北京、天津、石家庄等城市的空气质量较差；而海口、拉萨、南宁等城市的空气较好。

特别需要说明的是北京的空气污染与其它城市相比有很大的不同，在最短距离法中被单独聚为一类且与其它类相距较远，这与北京目前空气现状是相吻合的。

在本文的最后还根据实际情况对模型的优缺点做了评价，并指出了需要改进的地方。

关键词：大气污染；主成分分析；聚类分析1、数据资料本文的原始数据取自《中国统计年鉴，2014》，表1我国主要城市废气中主要污染物排放情况用1x 表示工业二氧化硫排放量，2x 表示工业二氧化硫排放量，3x 表示工业烟(粉)尘排放量，4x 表示生活二氧化硫排放量，5x 表示生活氮氧化物排放量，6x 表示生活烟尘排放量。

2、主成分分析2.1主成分分析的步骤（1）计算相关系数矩阵()ij m m R r ⨯=有~~11nkikjk ij aa r n =⋅=-∑（2）计算特征值和特征向量。

计算相关系数矩阵R 的特征值120m λλλ≥≥⋅⋅⋅≥，以及对应的特征向量12,,m u u u ⋅⋅⋅由特征值组成m 个新的指标变量：~~~12111211~~~12212222~~~1212,,,m m m m m m m m mm y u x u x u x y u x u x u x y u x u x u x =++⋅⋅⋅=++⋅⋅⋅=++⋅⋅⋅其中：1y 是第一主成分，2y 是第二主成分， ，m y 是第m 主成分。

（3）计算特征值的信息贡献率和累积贡献率。

1jj mkb λλ=∑为主成分j y 的信息贡献率，同时有11pkk p mkk λαλ===∑∑为主成分12,,,p y y y 的累积贡献率。

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

2.2 主成分分析构建评价指标定性地考虑反应各个城市空气质量的6个评价指标，不难看出某些指标可能存在较强的相关性，比如汽车的尾气中既含有二氧化硫也含有氮氧化物， 这两个指标之间可能存在相关性。

为了验证这个想法用MATLAB 计算指标之间的相关系数矩阵的特征值以及贡献率，如下表所示：表2主成分分析结果从结果中我们可以看出某些指标之间确实存在很强的相关性，比如生活烟尘对空气质量的贡献率接近0，这说明空气烟尘的值几乎可以由前面5个变量的值完全确定，这也是意料之中的结果。

如果直接用这些指标进行综合评价，必然造成信息重叠，影响评价的客观性，因此可以考虑用主成分分析的方法来进行分析。

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

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

表3 前3个主成分对应的特征向量于是得到3个主成分分别为1126212631260.42600.42920.29170.43090.31800.56060.17120.63430.316y x x x y x x x y x x x =+++=--++=-++- 分别以3个主成分的贡献率为权重，构建主成分综合评价模型，即1230.38270.26200.2395Z y y y =++把各地区的3个主成分值带入上式，可以得到各地空气质量的综合评价，见表4。

表4 排名和综合评价结果2.3结论从以上的综合评价结果可以看出。

2013年全国各地区废气中污染物情况存在较大的差异，重庆作为我国的国家中心城市是经济中心、金融中心和创新中心，同时也是全球著名的6大雾都城市年平均雾日是104天，在这次主成分的得分中位居第一，上海作为国际著名的经济中心，2014年GDP 总量居中国城市第一、亚洲第二，工商业活动频繁，空气质量较差。

北京作为我国的首都城市空气质量不容乐观，在31个城市中排名23；海口是海南省省会，空气质量排名第一，也在2011年世界卫生组织(WHO)发布首份全球城市空气污染调查报告，获中国空气最清洁的城市，其次拉萨、长沙、南宁等城市的空气质量也一直位居全国前列。

3、最短距离法聚类分析3.1 最短距离法的步骤选取主成分分析中的前3个主成分得分作为聚类指标，定义类与类之间的距离为两类最近样品间的距离，即,min K LKL ij i G j G D d ∈∈=（1） 计算n 个样品的距离矩阵(0)D 。

（2） 选择(0)D 中的最小元素，设为KL D 将,K L G G 合成一个新类记为M G 。

（3） 计算新类M G 与任一类之间的距离的递推公式为{},,,min minmin,min min{,}M JK JL JMJ ij ij iji G j G i G j G i G j G KJ LJ D d d d D D ∈∈∈∈∈∈===在(0)D 中,K L G G 所在的行和列合并成一个新行新列，对应M G ，该行列上的新距离由上述递推公式求得，其余行列上的距离值不变，得到新的距离矩阵 记作(1)D 。

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

3.2 最短距离法聚类模型(1) MATLAB 的算法流程图如下：图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、模型分析本模型的优点在于同时利用了主成分分析，最短距离法和离差平方和法对城市的空气质量进行综合评价，三种方法得到的结果相类似，都认为北京、天津等城市的空气较差，也与这些城市的实际空气情况相符，从而说明了模型的可靠性。

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

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,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、最小距离法程序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:Mfor j=1:Nif x0(i,j)>=m(i)m(i)=x0(i,j);endif x0(i,j)<=n(i)n(i)=x0(i,j);ends(i)=s(i)+x0(i,j);endeq(i)=s(i)/N;end%计算sigma，它是标准差的意思sigma0=zeros(M);for i=1:Mfor j=1:Nsigma0(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:Mfor j=1:Nx0_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);for i=1:Mfor j=1:Mfor k=1:Nd_abs(i,j)=d_abs(i,j)+abs(test(i,k)-test(j,k));endendendtest=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;endendendt=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;endg=zeros(1,M-2); %生成新的距离矩阵ii=0;for i=1:Mif(i==x||i==y);elseii=ii+1;g(ii)=(test(x,i)<test(y,i))*test(x,i)+(test(x,i)>=test(y,i))*test(y,i );end %...两列变一列：两两比较取其小者endtest([x,y],:)=[];test(:,[x,y])=[];test=cat(1,test,g);test=cat(2,test,[g';0]);M=length(test(1,:));enddendrogram(Z,0); %画树状聚类图title('最短距离聚类');ylabel('lambda');clear M MM i ii j x y g a t min str Z3、离差平方和法聚类程序clc,clearx=xlsread('主要城市废气中主要污染物排放情况 (2013年)',4,'B2:D32'); x0=x;y=pdist(x0,'cityblock'); yc=squareform(y)Z=linkage(y,'ward') [h,t]=dendrogram(Z)。