劉 敏,羅小玲,潘 新,張立倩
(內蒙古農業(yè)大學計算機與信息工程學院,內蒙古 呼和浩特 010018)
當前,我國經濟高速發(fā)展,但隨之而來的是環(huán)境空氣質量的下降[1-3],嚴重威脅人類健康,影響植物生長?!董h(huán)境空氣質量標準》(GB3095-2012)中規(guī)定參與環(huán)境空氣質量優(yōu)劣評價的主要污染物指標為SO2(二氧化硫)、N02(二氧化氮)、PM10(可吸入顆粒物)、C0(一氧化碳)、O3(臭氧)和PM2.5(細顆粒物)6項[4,5]。《環(huán)境空氣質量指數(AQI)技術規(guī)定(試行)》中指出空氣質量分為6個等級,6級屬于嚴重污染,1級屬于優(yōu),級別越高說明污染的情況越嚴重[6]。
近年來,國內外學者圍繞環(huán)境空氣質量的污染物指標[7-9]、預測模型[10,11]和地區(qū)差異[12,13]開展了廣泛的研究,如文獻[8]分析了火災排放對地表細顆粒物濃度和空氣質量的影響;文獻[9]根據可吸入顆粒物、二氧化氮和臭氧對歐洲6個城市的空氣污染情況進行評估;文獻[10,11]基于貝葉斯模型預測和診斷城市空氣質量;文獻[12]研究COVID-19大流行病前后印度3個城市空氣質量的差異;文獻[13]采用K均值聚類法對我國113個城市的空氣質量進行了區(qū)域性研究。上述文獻沒有結合污染物綜合指標研究區(qū)域性污染,提取環(huán)境空氣質量評價中的污染物綜合指標既可以降低問題的復雜度,又便于環(huán)境保護部門快速了解區(qū)域性污染來源,科學制定區(qū)域性大氣污染防治措施,為此,本文擬利用中國統(tǒng)計年鑒提供的2015-2019年全國31個主要城市環(huán)境空氣質量情況的155條樣本數據,采用主成分分析法(PCA)[14]研究表征環(huán)境空氣質量狀況的污染物綜合指標,并基于綜合指標結合譜系聚類法(HCM)對31個城市的環(huán)境空氣質量進行分類。
本文數據選自于《中國統(tǒng)計年鑒》,是關于2015年至2019年全國31個主要城市環(huán)境空氣質量情況的數據,共155個樣本,每個樣本有9項指標,依次為:city:城市;year:年份;x1:SO2(μg/m3);x2:N02(μg/m3);x3:PM10(μg/m3);x4:C0(mg/m3);x5:O3(μg/m3);x6:PM2.5(μg/m3);y:空氣質量達到及好于二級的天數(天)。使用的軟件為SAS9.0,部分樣本數據見表1。
表1 部分樣本數據
主成分分析(Principal Component Analysis)也稱PCA法,是處理多個具有相關性指標的一種統(tǒng)計方法。該方法運用降維的思想,通過正交變換對原始指標作線性組合,獲得盡可能少的互不相關的綜合指標即主成分去盡可能多地反映原始指標信息[15],設有n個樣品,每個樣品測p項指標,原始指標觀測數據陣記為X,每個觀測值記為xij,i=1,2,…,n,j=1,2,…,p,計算步驟如下:
1)計算X的相關陣R
i,j=1,2,…,p
(1)
當p個原始指標取值范圍彼此相差很大時,需要標準化X,記為X*,從X的相關陣出發(fā)計算主成分等價于標準化處理。
2)計算相關陣R的特征值及單位正交化特征向量
R的特征值按降序排列為:λ1≥λ2≥…λp>0,第i個主成分為Fi,λi是Fi的方差,特征值相應的正交化單位特征向量記為
X的第i個主成分為
(2)
3)選取主成分
4)解釋主成分
5)計算主成分得分
計算n個樣品在m個主成分上的得分
j=1,2,…,m
(3)
譜系聚類法(Hierarchical Clustering Method)也稱HCM法,用來研究樣品分類的一種統(tǒng)計方法,它的思想是用距離尺度衡量樣品之間的親疏程度并以此來實現(xiàn)分類[16]。設有n個樣品觀測值,每個觀測值測p項指標(變量),得到觀測數據xij,i=1,2,…,n,j=1,2,…,p,Xj=(X1j,X2j,…,Xnj)T表示第j項指標,X(i)=(Xi1,Xi2,…,Xip)表示第i個樣品,基本步驟如下:
1)標準化數據
(4)
(5)
(6)
2)計算n個樣品兩兩間的距離
可以使用明氏、蘭氏或者馬氏距離公式計算n個樣品兩兩間的距離,得樣品間的距離矩陣D(0)。定義樣品X(i)到樣品X(j)的歐氏距離為
(7)
開始每個樣品自成一類,此時Dij=dij。
3)合并類間距離最小的兩類為一新類
找出D(0)的非對角線最小元素,設為Dpq,則將Gp和Gq合并成一個新類Gr={Gp,Gq}。
4)計算新類Gr與其它類Gk的距離
可以使用最短距離法、類平均法、離差平方和法等方法計算新類Gr與其它類Gk的距離,以類平均法為例定義
(8)
將D(0)中第p、q行及p、q列合并成新行新列,新行新列對應Gr,此時距離陣記為D(1)。
5)對D(1)重復2)、3)兩步得D(2),如此下去,直到所有的元素并成一類為止。
6)繪制譜系聚類圖
以每一步合并類的最小類間距離為橫軸,樣品序號為縱軸,繪制橫向聚類圖,從聚類圖上可以清晰地描述各個類的樣本點。
7)決定分類個數及各類樣本點
偽F統(tǒng)計量用于評價分為k個類的效果,設已將n個樣品分為k類。
(9)
其中Pk為分類數為k個類時的總類內離差平方和,T為所有樣品或變量的總離差平方和,取偽F統(tǒng)計量較大而類數較小的聚類水平。
雙因素方差分析模型中選取y(空氣質量達到及好于二級的天數)為觀測值,year(年份)和city(城市)為因素A和因素B,這里主要分析因素A(year)對y的作用是否顯著,對于因素A(year)的顯著性F檢驗結果中,FA=10.85,p<0.0001,在0.05的顯著性水平下,因素A的作用顯著,說明2015年到2019年間31個城市環(huán)境空氣質量有顯著差異,接下來使用鄧肯法進一步做組間多重比較,比較結果見表2。
表2 鄧肯法組間多重比較結果
由表2可知,2019年31個城市空氣質量達到及好于二級的平均天數為284天,環(huán)境空氣質量與2015年至2018年顯著不同,是這5年中空氣質量最優(yōu)的,2018年次之,2015年最差,這說明全國環(huán)境空氣質量逐漸改善,呈持續(xù)向好局面,這得益于全國各地各部門積極優(yōu)化產業(yè)布局,節(jié)能減排,推進煤炭清潔化利用等一系列對大氣污染聯(lián)防聯(lián)控的措施。
基于方差分析的結論,PCA法選取了2019年全國31個主要城市環(huán)境空氣質量的數據,樣本數為31個,指標為6項,分別是x1-x6(即SO2、N02、PM10、C0、O3、PM2.5)。
1)標準化數據
由表1可知,樣本數據6項指標取值范圍彼此相差大,所以首先利用式(6)對樣本數據進行標準化,標準化后的部分數據見表3。
2)計算相關陣R
為了驗證6項污染物指標的相關性,需要利用式(1)計算原始指標pearson相關陣R,計算結果如下
根據相關陣R可知SO2(x1)與C0(x4)相關性最強,相關系數為0.66411;N02(x2)與PM10(x3)、PM2.5(x6)相關性最強,相關系數為0.76444和0.76980;PM10(x3)與PM2.5(x6)相關性最強,相關系數為0.91405,指標之間存在一定的相關性。
3)計算相關陣R特征值和特征向量
從相關陣R出發(fā),計算特征值和特征向量,并進一步找到主成分。相關陣R的特征值和主成分貢獻率見表4。
表4 相關陣的特征值和主成分貢獻率
由表4可知,第一、第二主成分的貢獻率分別為63.59%和21.25%,累積貢獻率為84.83%。
4)選取并解釋主成分
按照累積貢獻率達到85%的原則,本文選取前2個主成分,這2個主成分可以反映原始指標84.84%的信息量,約等于85%。前2個主成分的特征向量見表5。
表5 前2個主成分的特征向量
由表5及式(2)可以寫出主成分的表達式如下(保留3位小數):
在第一主成分F1中,x2、x3和x6的系數絕對值是最大的前三項,分別是0.446、0.491和0.466,因此F1主要綜合了N02、PM10和PM2.53項污染物指標,PM10主要來自燃煤排放的煙塵、建筑工地和地面揚起的灰塵等一次污染物,PM2.5主要來自二次顆粒物[17],可以把F1稱為顆粒物污染綜合指標,它能反映原始指標63.59%的信息量;在第二主成分F2中,x1、x4和x5的系數絕對值分別是0.631、0.464和0.569,因此F2主要綜合了SO2、C0和O33項污染物指標,二氧化硫主要來自燃燒廢氣,氮氧化物主要來自汽車尾氣[18],可以把F2稱為廢氣污染綜合指標,它能反映原始指標21.25%的信息量。
5)計算主成分得分并繪制主成分散點圖
將31個城市的6項污染物指標的觀測數據標準化后分別代入兩個主成分表達式,利用式(3)計算每個城市的主成分得分并按降序輸出,輸出結果見表6。
表6 部分主成分得分降序排列結果
由表6可知,石家莊、太原、濟南、鄭州和西安這5個城市主成分得分較高,排在前5名,說明這5個城市顆粒物和廢氣污染較嚴重;昆明、貴陽、福州、??诤屠_這5個城市主成分得分較低,空氣質量好。以第一主成分為縱軸,第二主成分為橫軸,繪制31個城市的主成分得分散點圖,如圖1所示。
圖1 31個城市主成分得分散點圖
從圖1可以看出來,散點圖越靠左上角的地區(qū),顆粒物污染越嚴重(以下結論給出的城市名稱均按污染程度遞減排序),如顆粒物污染最嚴重的城市有6個,分別是石家莊、太原、濟南、鄭州、西安和天津;較嚴重的城市有12個,分別是武漢、南京、北京、合肥、杭州、成都、長沙、廣州、重慶、南昌、長春和上海;較輕的城市6個,分別是南寧、昆明、貴陽、福州、??诤屠_。越靠右側的地區(qū),廢氣污染越嚴重,廢氣污染較嚴重的7個城市有西寧、蘭州、沈陽、哈爾濱、呼和浩特、銀川和烏魯木齊。越靠右上角的地區(qū),顆粒物廢氣污染越嚴重,相對來說,石家莊和太原這兩個城市環(huán)境空氣質量較差。越靠左下角的地區(qū),顆粒物廢氣污染越少,環(huán)境空氣質量越好,相較于其它城市,福州和海口的環(huán)境更宜人。
根據前2個主成分對2019年31個城市的環(huán)境空氣質量數據進行主成分聚類分析,聚類歷史的輸出結果見表7,在類別控制在4類以下的前提下,利用式(9)計算的偽F統(tǒng)計量最大和次大依次為49.4和24.5,建議分為4類或3類是較合適的;偽T2最大和次大依次為47和21,建議分為4類或2類;半偏R2最大和次大依次為0.4199和0.2162,建議分為2類或3類;R2最大和次大依次為0.846和0.636,建立分為4類或3類,綜合以上統(tǒng)計量及主成分得分的信息,最終決定分為4類能較準確地體現(xiàn)城市環(huán)境空氣質量的區(qū)域特性,分類結果如圖2所示。
圖2 31個城市類平均法橫向聚類圖
表7 聚類歷史輸出結果
在圖2的聚類圖上進行標識,可以看出,第一類城市群有{濟南、天津、西安、鄭州、石家莊、太原};第二類城市群有{哈爾濱、銀川、西寧、呼和浩特、沈陽、烏魯木齊、蘭州};第三類城市群有{北京、成都、杭州、合肥、廣州、長沙、重慶、南京、武漢、南昌、長春、上海};第四類城市群有{福州、???、拉薩、貴陽、昆明、南寧},與主成分得分的散點圖分析結論一致。綜合比較而言,第一類城市群是大氣污染的重災區(qū)[19],這些地區(qū)在京津冀周邊,屬于我國內陸城市,沙塵天氣多,城市大風日數較少,不利于污染物擴散,冬季通過燃煤取暖,產業(yè)結構以重工業(yè)為主,復合型大氣污染比較突出,尤其是石家莊和太原兩個城市,顆粒物和廢棄污染問題嚴峻。第二類城市群大部分位于我國的西北部,氣候干燥,春秋風沙大,雖然能源結構也是以煤炭為主,但是相較于前兩類城市群,經濟欠發(fā)達,地廣人稀,汽車保有量逐年增加,目前亟需解決的是廢氣污染。第三類城市群中大部分城市屬于長三角區(qū)域,處于我國南方,常年雨水多,對空氣污染能起到一定減少的作用,但是這些地區(qū)土地面積狹小,資源消耗大,人類活動強度高,目前主要面臨顆粒物污染[20]。第四類城市群環(huán)境宜人,有的城市依江面海,自然植被密集,有的城市海拔高,全年日照時間長,人口密度低,這些都有利于形成優(yōu)良的城市環(huán)境空氣質量。
本文采用雙因素無交互作用的方差分析法、主成分分析法(PCA)和譜系聚類法(HCM)探討了2015-2019這五年不同時期全國環(huán)境空氣質量的總體差異、影響環(huán)境空氣質量的污染物綜合指標以及城市環(huán)境空氣質量的區(qū)域性特征,結果表明:
1)2019年全國空氣質量達到及好于二級的平均天數為284天,與2015年-2018年顯著不同,環(huán)境空氣質量最優(yōu);2018年與2016年、2015年與2017年的環(huán)境空氣質量無差異;2015年環(huán)境空氣質量最差。
2)影響環(huán)境空氣質量污染物的6項原始指標之間具有不同程度的相關性,經過PCA法獲得了顆粒物和廢氣2個主成分,分別提取原始指標63.59%和21.25%的信息量,累積提取原始指標約85%的信息量;
3)基于顆粒物和廢氣2個主成分,使用HCM法對31個城市環(huán)境空氣質量由好至差分為了4類,第一類城市群有{濟南、天津、西安、鄭州、石家莊、太原};第二類城市群有{哈爾濱、銀川、西寧、呼和浩特、沈陽、烏魯木齊、蘭州};第三類城市群有{北京、成都、杭州、合肥、廣州、長沙、重慶、南京、武漢、南昌、長春、上海};第四類城市群有{福州、???、拉薩、貴陽、昆明、南寧}。
1)本文的研究方法將環(huán)境空氣質量影響因素的問題研究空間從6維降到了2維,雖然損失了15%的信息量,但是卻保留了85%的信息量,抓住了主要矛盾,
而且明顯降低了問題的復雜性。
2)基于環(huán)境空氣質量評價中的污染物綜合指標可以客觀準確地反映區(qū)域性空氣污染現(xiàn)狀,為改善和控制區(qū)域性空氣質量提供理論依據。
3)可以繼續(xù)研究基于主成分聚類的判別分析,這將對城市環(huán)境空氣質量的預判和科學精準地防治環(huán)境空氣污染有一定的借鑒意義。