主要城市廢氣中主要污染物排放情況(2013年)

1、我國主要城市廢氣中主要污染物排放情況摘要近幾年來環(huán)境問題成為全社會極為關(guān)注的熱點, 空氣污染是其中最熱門的話題，同時也是最重要的民生問題。本文針對這個現(xiàn)狀，搜集了全國有代表性的31個城市的主要大氣污染物的排放情況，先利用主成分分析評價了31個城市的綜合空氣質(zhì)量，然后又分別用最短距離法和離差平方和法進(jìn)行聚類分析，最終結(jié)果為北京、天津、石家莊等城市的空氣質(zhì)量較差；而?？凇⒗_、南寧等城市的空氣較好。特別需要說明的是北京的空氣污染與其它城市相比有很大的不同，在最短距離法中被單獨聚為一類且與其它類相距較遠(yuǎn)，這與北京目前空氣現(xiàn)狀是相吻合的。在本文的最后還根據(jù)實際情況對模型的優(yōu)缺點做了評價，并指出了需要改

2、進(jìn)的地方。關(guān)鍵詞：大氣污染；主成分分析；聚類分析1、數(shù)據(jù)資料本文的原始數(shù)據(jù)取自中國統(tǒng)計年鑒，2014，表1 我國主要城市廢氣中主要污染物排放情況地區(qū)工業(yè)二氧化硫工業(yè)氮氧化物工業(yè)煙(粉)塵生活二氧化硫生活氮氧化物生活煙塵 北 京520417592727182349671363828258 天 津207793250646627668959522118400 石 家莊17646920030199806956428026635 太 原88880960183700333396673826727呼和浩特961901316654882242576653763 沈 陽130672833486042514389

3、515415276 長 春572469519072970734415457919 哈 爾濱659878551582323500122298580792 上 海1728672623466717442947234746451 南 京1106651096936525617504001000 杭 州820216728340243633335135 合 肥41483703114238727101303188 福 州7604372284434831279169547 南 昌40756185971141364158254 濟(jì) 南8111872969471172608736298355 鄭 州10612313

4、4120338281197517809150 武 漢962229561220020572014161001 長 沙21173159511954523661532946 廣 州655895716416660663276214 南 寧330453479720950874810684631 海 口179886114911175 重 慶4944152479051798425326144874401 成 都52040444112145248912109661 貴 陽7060330450242333549317535530 昆 明10284268213573665263970328 拉 薩93020165

5、3867840199 西 安691033491715893238311095114012 蘭 州721487991540109741319501088 西 寧718395328052765712914194793 銀 川923698432127170569712373016烏魯木齊7421611380352441669114254920用表示工業(yè)二氧化硫排放量，表示工業(yè)二氧化硫排放量，表示工業(yè)煙(粉)塵排放量，表示生活二氧化硫排放量，表示生活氮氧化物排放量，表示生活煙塵排放量。2、主成分分析2.1主成分分析的步驟（1）計算相關(guān)系數(shù)矩陣有（2）計算特征值和特征向量。計算相關(guān)系數(shù)矩陣R的特征值，以

6、及對應(yīng)的特征向量由特征值組成m個新的指標(biāo)變量：其中：是第一主成分，是第二主成分，是第m主成分。（3）計算特征值的信息貢獻(xiàn)率和累積貢獻(xiàn)率。為主成分的信息貢獻(xiàn)率，同時有為主成分的累積貢獻(xiàn)率。（4）根據(jù)累積貢獻(xiàn)率選取幾個主成分作為新的評價指標(biāo)。2.2 主成分分析構(gòu)建評價指標(biāo)定性地考慮反應(yīng)各個城市空氣質(zhì)量的6個評價指標(biāo)， 不難看出某些指標(biāo)可能存在較強(qiáng)的相關(guān)性，比如汽車的尾氣中既含有二氧化硫也含有氮氧化物， 這兩個指標(biāo)之間可能存在相關(guān)性。為了驗證這個想法用MATLAB計算指標(biāo)之間的相關(guān)系數(shù)矩陣的特征值以及貢獻(xiàn)率，如下表所示：表2 主成分分析結(jié)果序號特征值貢獻(xiàn)率/%累積貢獻(xiàn)率/%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 從結(jié)果中我們可以看出某些指標(biāo)之間確實存在很強(qiáng)的相關(guān)性，比如生活煙塵對空氣質(zhì)量的貢獻(xiàn)率接近0，這說明空氣煙塵的值幾乎可以由前面5個變量的值完全確定，這也是意料之中的結(jié)果。如果直接用這些指標(biāo)進(jìn)行綜合評價，必然造成信息重疊，影響評價的客觀性，因此可以考慮用主成分分析的方法來進(jìn)行分析。從主成分分析結(jié)果中可以看出前3個主成分的累積貢獻(xiàn)率就達(dá)到了88.42%，因此選取前3個主成分進(jìn)行綜合評價。前3個特征值對應(yīng)的特征向量為

8、見表3。表3 前3個主成分對應(yīng)的特征向量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個主成分的貢獻(xiàn)率為權(quán)重，構(gòu)建主成分綜合評價模型，即把各地區(qū)的3個主成分值帶入上式，可以得到各地空氣質(zhì)量的綜合評價，見表4。表4 排名和綜合評價結(jié)果2.3 結(jié)論從以上的綜合評價結(jié)果可以看出。2013年全國各地區(qū)廢氣中污

9、染物情況存在較大的差異，重慶作為我國的國家中心城市是經(jīng)濟(jì)中心、金融中心和創(chuàng)新中心，同時也是全球著名的6大霧都城市年平均霧日是104天，在這次主成分的得分中位居第一，上海作為國際著名的經(jīng)濟(jì)中心，2014年GDP總量居中國城市第一、亞洲第二，工商業(yè)活動頻繁，空氣質(zhì)量較差。北京作為我國的首都城市空氣質(zhì)量不容樂觀，在31個城市中排名23；海口是海南省省會，空氣質(zhì)量排名第一，也在2011年世界衛(wèi)生組織(WHO)發(fā)布首份全球城市空氣污染調(diào)查報告，獲中國空氣最清潔的城市，其次拉薩、長沙、南寧等城市的空氣質(zhì)量也一直位居全國前列。3、最短距離法聚類分析3.1 最短距離法的步驟選取主成分分析中的前3個主成分得分作

10、為聚類指標(biāo)，定義類與類之間的距離為兩類最近樣品間的距離，即（1） 計算n個樣品的距離矩陣。（2） 選擇中的最小元素，設(shè)為將合成一個新類記為。（3） 計算新類與任一類之間的距離的遞推公式為在中所在的行和列合并成一個新行新列，對應(yīng)，該行列上的新距離由上述遞推公式求得，其余行列上的距離值不變，得到新的距離矩陣 記作。（4）對重復(fù)上述的2步得到，如此下去直到所有的元素合并成一類為止。3.2 最短距離法聚類模型(1) MATLAB的算法流程圖如下：圖1 程序流程圖開始讀取excel文件計算（31-）個樣本的距離D（-1）找出D（-1）的最小值并返回其行row和列col將樣本row和col合成一類=+1是

11、否聚類結(jié)束畫聚類圖(2)計算結(jié)果第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

12、=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(

13、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.22

14、048 合并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） 結(jié)論從圖中可以看出，如果根據(jù)空氣質(zhì)量把31個地區(qū)分為3類結(jié)果為：第一類北京第二類天津石家莊呼和浩特第三類其他城市如果分為4類結(jié)果為：第一類北京第二類呼和浩特第三類天津石家莊第四類其他城市從聚類結(jié)果可以看出，北京的空氣質(zhì)量

15、與其他城市相比有很大的不同，可能的原因是生活煙塵的排放量過高，這與北京近年來的空氣質(zhì)量狀況是相符的，其次石家莊、天津、太原等城市的空氣狀況和北京類似。4、離差平方和法聚類分析直接調(diào)用MATLAB中的函數(shù)，得到聚類圖如下圖所示。從圖中可以看出離差平方和法得到的結(jié)果仍把北京、天津、石家莊、呼和浩特歸為一類，與最短距離法得到的結(jié)果相同。不同點在于離差平方和法使得兩個大的類不容易合并，兩個小的類容易合并，因而能達(dá)到分離的開的聚類結(jié)果，這符合我們對聚類的實際要求，因而更具有優(yōu)越性。5、模型分析本模型的優(yōu)點在于同時利用了主成分分析，最短距離法和離差平方和法對城市的空氣質(zhì)量進(jìn)行綜合評價， 三種方法得到的結(jié)果

16、相類似，都認(rèn)為北京、天津等城市的空氣較差， 也與這些城市的實際空氣情況相符，從而說明了模型的可靠性。本模型的缺點是只考慮了城市的污染物排放量，并沒有考慮每個城市的面積，而且實際情況中相鄰近的城市還可能存在污染物的相互傳播問題，進(jìn)一步的改進(jìn)模型中可以用污染物排放量除以城市的面積作為指標(biāo)變量進(jìn)行分析。6、參考文獻(xiàn)1 王學(xué)民應(yīng)用多元分析M上海財經(jīng)大學(xué)出版社，2014：153-197.2 司守奎王數(shù)學(xué)建模算法與應(yīng)用M國防工業(yè)出版社，2014：193-207.7、附件(1) 主要城市廢氣中主要污染物排放情況 (2013年).xls(2) principal_component_analysis.m(3)

17、 Single_linkage_method.m(4) Ward_method.m8、附錄1、 主成分分析程序。clc,clearx=xlsread(主要城市廢氣中主要污染物排放情況 (2013年),2);x0=x;r=corrcoef(x0);%求相關(guān)系數(shù)矩陣vec1,lamda,rate=pcacov(r)%求相關(guān)系數(shù)矩陣的特征值以及特征向量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);%把得分按高

18、到低排序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:M for j=1:N if x0(i,j)=m(i) m(i)=x0(i,j); end if x0(i,j)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),),記