基于PCA＿HCM的UAQ影響因素及區(qū)域性研究

劉 敏,羅小玲,潘 新,張立倩

(內(nèi)蒙古農(nóng)業(yè)大學(xué)計算機與信息工程學(xué)院,內(nèi)蒙古 呼和浩特 010018)

1 引言

當(dāng)前,我國經(jīng)濟高速發(fā)展,但隨之而來的是環(huán)境空氣質(zhì)量的下降[1-3],嚴(yán)重威脅人類健康,影響植物生長。《環(huán)境空氣質(zhì)量標(biāo)準(zhǔn)》(GB3095-2012)中規(guī)定參與環(huán)境空氣質(zhì)量優(yōu)劣評價的主要污染物指標(biāo)為SO2(二氧化硫)、N02(二氧化氮)、PM10(可吸入顆粒物)、C0(一氧化碳)、O3(臭氧)和PM2.5(細(xì)顆粒物)6項[4,5]。《環(huán)境空氣質(zhì)量指數(shù)(AQI)技術(shù)規(guī)定(試行)》中指出空氣質(zhì)量分為6個等級,6級屬于嚴(yán)重污染,1級屬于優(yōu),級別越高說明污染的情況越嚴(yán)重[6]。

近年來,國內(nèi)外學(xué)者圍繞環(huán)境空氣質(zhì)量的污染物指標(biāo)[7-9]、預(yù)測模型[10,11]和地區(qū)差異[12,13]開展了廣泛的研究,如文獻(xiàn)[8]分析了火災(zāi)排放對地表細(xì)顆粒物濃度和空氣質(zhì)量的影響;文獻(xiàn)[9]根據(jù)可吸入顆粒物、二氧化氮和臭氧對歐洲6個城市的空氣污染情況進行評估;文獻(xiàn)[10,11]基于貝葉斯模型預(yù)測和診斷城市空氣質(zhì)量;文獻(xiàn)[12]研究COVID-19大流行病前后印度3個城市空氣質(zhì)量的差異;文獻(xiàn)[13]采用K均值聚類法對我國113個城市的空氣質(zhì)量進行了區(qū)域性研究。上述文獻(xiàn)沒有結(jié)合污染物綜合指標(biāo)研究區(qū)域性污染,提取環(huán)境空氣質(zhì)量評價中的污染物綜合指標(biāo)既可以降低問題的復(fù)雜度,又便于環(huán)境保護部門快速了解區(qū)域性污染來源,科學(xué)制定區(qū)域性大氣污染防治措施,為此,本文擬利用中國統(tǒng)計年鑒提供的2015-2019年全國31個主要城市環(huán)境空氣質(zhì)量情況的155條樣本數(shù)據(jù),采用主成分分析法(PCA)[14]研究表征環(huán)境空氣質(zhì)量狀況的污染物綜合指標(biāo),并基于綜合指標(biāo)結(jié)合譜系聚類法(HCM)對31個城市的環(huán)境空氣質(zhì)量進行分類。

2 材料與方法

2.1 數(shù)據(jù)來源

本文數(shù)據(jù)選自于《中國統(tǒng)計年鑒》,是關(guān)于2015年至2019年全國31個主要城市環(huán)境空氣質(zhì)量情況的數(shù)據(jù),共155個樣本,每個樣本有9項指標(biāo),依次為: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:空氣質(zhì)量達(dá)到及好于二級的天數(shù)(天)。使用的軟件為SAS9.0,部分樣本數(shù)據(jù)見表1。

表1 部分樣本數(shù)據(jù)

2.2 主成分分析法

主成分分析(Principal Component Analysis)也稱PCA法,是處理多個具有相關(guān)性指標(biāo)的一種統(tǒng)計方法。該方法運用降維的思想,通過正交變換對原始指標(biāo)作線性組合,獲得盡可能少的互不相關(guān)的綜合指標(biāo)即主成分去盡可能多地反映原始指標(biāo)信息[15],設(shè)有n個樣品,每個樣品測p項指標(biāo),原始指標(biāo)觀測數(shù)據(jù)陣記為X,每個觀測值記為xij,i=1,2,…,n,j=1,2,…,p,計算步驟如下:

1)計算X的相關(guān)陣R

i,j=1,2,…,p

(1)

當(dāng)p個原始指標(biāo)取值范圍彼此相差很大時,需要標(biāo)準(zhǔn)化X,記為X*,從X的相關(guān)陣出發(fā)計算主成分等價于標(biāo)準(zhǔn)化處理。

2)計算相關(guān)陣R的特征值及單位正交化特征向量

R的特征值按降序排列為:λ1≥λ2≥…λp>0,第i個主成分為Fi,λi是Fi的方差,特征值相應(yīng)的正交化單位特征向量記為

X的第i個主成分為

(2)

3)選取主成分

4)解釋主成分

5)計算主成分得分

計算n個樣品在m個主成分上的得分

j=1,2,…,m

(3)

2.3 譜系聚類法

譜系聚類法(Hierarchical Clustering Method)也稱HCM法,用來研究樣品分類的一種統(tǒng)計方法,它的思想是用距離尺度衡量樣品之間的親疏程度并以此來實現(xiàn)分類[16]。設(shè)有n個樣品觀測值,每個觀測值測p項指標(biāo)(變量),得到觀測數(shù)據(jù)xij,i=1,2,…,n,j=1,2,…,p,Xj=(X1j,X2j,…,Xnj)T表示第j項指標(biāo),X(i)=(Xi1,Xi2,…,Xip)表示第i個樣品,基本步驟如下:

1)標(biāo)準(zhǔn)化數(shù)據(jù)

(4)

(5)

(6)

2)計算n個樣品兩兩間的距離

可以使用明氏、蘭氏或者馬氏距離公式計算n個樣品兩兩間的距離,得樣品間的距離矩陣D(0)。定義樣品X(i)到樣品X(j)的歐氏距離為

(7)

開始每個樣品自成一類,此時Dij=dij。

3)合并類間距離最小的兩類為一新類

找出D(0)的非對角線最小元素,設(shè)為Dpq,則將Gp和Gq合并成一個新類Gr={Gp,Gq}。

4)計算新類Gr與其它類Gk的距離

可以使用最短距離法、類平均法、離差平方和法等方法計算新類Gr與其它類Gk的距離,以類平均法為例定義

(8)

將D(0)中第p、q行及p、q列合并成新行新列,新行新列對應(yīng)Gr,此時距離陣記為D(1)。

5)對D(1)重復(fù)2)、3)兩步得D(2),如此下去,直到所有的元素并成一類為止。

6)繪制譜系聚類圖

以每一步合并類的最小類間距離為橫軸,樣品序號為縱軸,繪制橫向聚類圖,從聚類圖上可以清晰地描述各個類的樣本點。

7)決定分類個數(shù)及各類樣本點

偽F統(tǒng)計量用于評價分為k個類的效果,設(shè)已將n個樣品分為k類。

(9)

其中Pk為分類數(shù)為k個類時的總類內(nèi)離差平方和,T為所有樣品或變量的總離差平方和,取偽F統(tǒng)計量較大而類數(shù)較小的聚類水平。

3 結(jié)果與分析

3.1 不同時期環(huán)境空氣質(zhì)量差異的方差分析

雙因素方差分析模型中選取y(空氣質(zhì)量達(dá)到及好于二級的天數(shù))為觀測值,year(年份)和city(城市)為因素A和因素B,這里主要分析因素A(year)對y的作用是否顯著,對于因素A(year)的顯著性F檢驗結(jié)果中,FA=10.85,p<0.0001,在0.05的顯著性水平下,因素A的作用顯著,說明2015年到2019年間31個城市環(huán)境空氣質(zhì)量有顯著差異,接下來使用鄧肯法進一步做組間多重比較,比較結(jié)果見表2。

表2 鄧肯法組間多重比較結(jié)果

由表2可知,2019年31個城市空氣質(zhì)量達(dá)到及好于二級的平均天數(shù)為284天,環(huán)境空氣質(zhì)量與2015年至2018年顯著不同,是這5年中空氣質(zhì)量最優(yōu)的,2018年次之,2015年最差,這說明全國環(huán)境空氣質(zhì)量逐漸改善,呈持續(xù)向好局面,這得益于全國各地各部門積極優(yōu)化產(chǎn)業(yè)布局,節(jié)能減排,推進煤炭清潔化利用等一系列對大氣污染聯(lián)防聯(lián)控的措施。

3.2 使用PCA法分析環(huán)境空氣質(zhì)量污染物綜合指標(biāo)

基于方差分析的結(jié)論,PCA法選取了2019年全國31個主要城市環(huán)境空氣質(zhì)量的數(shù)據(jù),樣本數(shù)為31個,指標(biāo)為6項,分別是x1-x6(即SO2、N02、PM10、C0、O3、PM2.5)。

1)標(biāo)準(zhǔn)化數(shù)據(jù)

由表1可知,樣本數(shù)據(jù)6項指標(biāo)取值范圍彼此相差大,所以首先利用式(6)對樣本數(shù)據(jù)進行標(biāo)準(zhǔn)化,標(biāo)準(zhǔn)化后的部分?jǐn)?shù)據(jù)見表3。

2)計算相關(guān)陣R

為了驗證6項污染物指標(biāo)的相關(guān)性,需要利用式(1)計算原始指標(biāo)pearson相關(guān)陣R,計算結(jié)果如下

根據(jù)相關(guān)陣R可知SO2(x1)與C0(x4)相關(guān)性最強,相關(guān)系數(shù)為0.66411;N02(x2)與PM10(x3)、PM2.5(x6)相關(guān)性最強,相關(guān)系數(shù)為0.76444和0.76980;PM10(x3)與PM2.5(x6)相關(guān)性最強,相關(guān)系數(shù)為0.91405,指標(biāo)之間存在一定的相關(guān)性。

3)計算相關(guān)陣R特征值和特征向量

從相關(guān)陣R出發(fā),計算特征值和特征向量,并進一步找到主成分。相關(guān)陣R的特征值和主成分貢獻(xiàn)率見表4。

表4 相關(guān)陣的特征值和主成分貢獻(xiàn)率

由表4可知,第一、第二主成分的貢獻(xiàn)率分別為63.59%和21.25%,累積貢獻(xiàn)率為84.83%。

4)選取并解釋主成分

按照累積貢獻(xiàn)率達(dá)到85%的原則,本文選取前2個主成分,這2個主成分可以反映原始指標(biāo)84.84%的信息量,約等于85%。前2個主成分的特征向量見表5。

表5 前2個主成分的特征向量

由表5及式(2)可以寫出主成分的表達(dá)式如下(保留3位小數(shù)):

在第一主成分F1中,x2、x3和x6的系數(shù)絕對值是最大的前三項,分別是0.446、0.491和0.466,因此F1主要綜合了N02、PM10和PM2.53項污染物指標(biāo),PM10主要來自燃煤排放的煙塵、建筑工地和地面揚起的灰塵等一次污染物,PM2.5主要來自二次顆粒物[17],可以把F1稱為顆粒物污染綜合指標(biāo),它能反映原始指標(biāo)63.59%的信息量;在第二主成分F2中,x1、x4和x5的系數(shù)絕對值分別是0.631、0.464和0.569,因此F2主要綜合了SO2、C0和O33項污染物指標(biāo),二氧化硫主要來自燃燒廢氣,氮氧化物主要來自汽車尾氣[18],可以把F2稱為廢氣污染綜合指標(biāo),它能反映原始指標(biāo)21.25%的信息量。

5)計算主成分得分并繪制主成分散點圖

將31個城市的6項污染物指標(biāo)的觀測數(shù)據(jù)標(biāo)準(zhǔn)化后分別代入兩個主成分表達(dá)式,利用式(3)計算每個城市的主成分得分并按降序輸出,輸出結(jié)果見表6。

表6 部分主成分得分降序排列結(jié)果

由表6可知,石家莊、太原、濟南、鄭州和西安這5個城市主成分得分較高,排在前5名,說明這5個城市顆粒物和廢氣污染較嚴(yán)重;昆明、貴陽、福州、海口和拉薩這5個城市主成分得分較低,空氣質(zhì)量好。以第一主成分為縱軸,第二主成分為橫軸,繪制31個城市的主成分得分散點圖,如圖1所示。

圖1 31個城市主成分得分散點圖

從圖1可以看出來,散點圖越靠左上角的地區(qū),顆粒物污染越嚴(yán)重(以下結(jié)論給出的城市名稱均按污染程度遞減排序),如顆粒物污染最嚴(yán)重的城市有6個,分別是石家莊、太原、濟南、鄭州、西安和天津;較嚴(yán)重的城市有12個,分別是武漢、南京、北京、合肥、杭州、成都、長沙、廣州、重慶、南昌、長春和上海;較輕的城市6個,分別是南寧、昆明、貴陽、福州、海口和拉薩。越靠右側(cè)的地區(qū),廢氣污染越嚴(yán)重,廢氣污染較嚴(yán)重的7個城市有西寧、蘭州、沈陽、哈爾濱、呼和浩特、銀川和烏魯木齊。越靠右上角的地區(qū),顆粒物廢氣污染越嚴(yán)重,相對來說,石家莊和太原這兩個城市環(huán)境空氣質(zhì)量較差。越靠左下角的地區(qū),顆粒物廢氣污染越少,環(huán)境空氣質(zhì)量越好,相較于其它城市,福州和海口的環(huán)境更宜人。

3.3 使用HCM法基于污染物綜合指標(biāo)進行聚類分析

根據(jù)前2個主成分對2019年31個城市的環(huán)境空氣質(zhì)量數(shù)據(jù)進行主成分聚類分析,聚類歷史的輸出結(jié)果見表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類能較準(zhǔn)確地體現(xiàn)城市環(huán)境空氣質(zhì)量的區(qū)域特性,分類結(jié)果如圖2所示。

圖2 31個城市類平均法橫向聚類圖

表7 聚類歷史輸出結(jié)果

在圖2的聚類圖上進行標(biāo)識,可以看出,第一類城市群有{濟南、天津、西安、鄭州、石家莊、太原};第二類城市群有{哈爾濱、銀川、西寧、呼和浩特、沈陽、烏魯木齊、蘭州};第三類城市群有{北京、成都、杭州、合肥、廣州、長沙、重慶、南京、武漢、南昌、長春、上海};第四類城市群有{福州、海口、拉薩、貴陽、昆明、南寧},與主成分得分的散點圖分析結(jié)論一致。綜合比較而言,第一類城市群是大氣污染的重災(zāi)區(qū)[19],這些地區(qū)在京津冀周邊,屬于我國內(nèi)陸城市,沙塵天氣多,城市大風(fēng)日數(shù)較少,不利于污染物擴散,冬季通過燃煤取暖,產(chǎn)業(yè)結(jié)構(gòu)以重工業(yè)為主,復(fù)合型大氣污染比較突出,尤其是石家莊和太原兩個城市,顆粒物和廢棄污染問題嚴(yán)峻。第二類城市群大部分位于我國的西北部,氣候干燥,春秋風(fēng)沙大,雖然能源結(jié)構(gòu)也是以煤炭為主,但是相較于前兩類城市群,經(jīng)濟欠發(fā)達(dá),地廣人稀,汽車保有量逐年增加,目前亟需解決的是廢氣污染。第三類城市群中大部分城市屬于長三角區(qū)域,處于我國南方,常年雨水多,對空氣污染能起到一定減少的作用,但是這些地區(qū)土地面積狹小,資源消耗大,人類活動強度高,目前主要面臨顆粒物污染[20]。第四類城市群環(huán)境宜人,有的城市依江面海,自然植被密集,有的城市海拔高,全年日照時間長,人口密度低,這些都有利于形成優(yōu)良的城市環(huán)境空氣質(zhì)量。

4 結(jié)論

本文采用雙因素?zé)o交互作用的方差分析法、主成分分析法(PCA)和譜系聚類法(HCM)探討了2015-2019這五年不同時期全國環(huán)境空氣質(zhì)量的總體差異、影響環(huán)境空氣質(zhì)量的污染物綜合指標(biāo)以及城市環(huán)境空氣質(zhì)量的區(qū)域性特征,結(jié)果表明:

1)2019年全國空氣質(zhì)量達(dá)到及好于二級的平均天數(shù)為284天,與2015年-2018年顯著不同,環(huán)境空氣質(zhì)量最優(yōu);2018年與2016年、2015年與2017年的環(huán)境空氣質(zhì)量無差異;2015年環(huán)境空氣質(zhì)量最差。

2)影響環(huán)境空氣質(zhì)量污染物的6項原始指標(biāo)之間具有不同程度的相關(guān)性,經(jīng)過PCA法獲得了顆粒物和廢氣2個主成分,分別提取原始指標(biāo)63.59%和21.25%的信息量,累積提取原始指標(biāo)約85%的信息量;

3)基于顆粒物和廢氣2個主成分,使用HCM法對31個城市環(huán)境空氣質(zhì)量由好至差分為了4類,第一類城市群有{濟南、天津、西安、鄭州、石家莊、太原};第二類城市群有{哈爾濱、銀川、西寧、呼和浩特、沈陽、烏魯木齊、蘭州};第三類城市群有{北京、成都、杭州、合肥、廣州、長沙、重慶、南京、武漢、南昌、長春、上海};第四類城市群有{福州、海口、拉薩、貴陽、昆明、南寧}。

5 討論

1)本文的研究方法將環(huán)境空氣質(zhì)量影響因素的問題研究空間從6維降到了2維,雖然損失了15%的信息量,但是卻保留了85%的信息量,抓住了主要矛盾,

而且明顯降低了問題的復(fù)雜性。

2)基于環(huán)境空氣質(zhì)量評價中的污染物綜合指標(biāo)可以客觀準(zhǔn)確地反映區(qū)域性空氣污染現(xiàn)狀,為改善和控制區(qū)域性空氣質(zhì)量提供理論依據(jù)。

3)可以繼續(xù)研究基于主成分聚類的判別分析,這將對城市環(huán)境空氣質(zhì)量的預(yù)判和科學(xué)精準(zhǔn)地防治環(huán)境空氣污染有一定的借鑒意義。