盧 彬，馬 行，穆春陽(yáng)，張 鄂

(1.北方民族大學(xué)寧夏智能信息與大數(shù)據(jù)處理重點(diǎn)實(shí)驗(yàn)室，寧夏 銀川 750021；2.北方民族大學(xué)電氣信息工程學(xué)院，寧夏 銀川 750021；3.北方民族大學(xué)機(jī)電工程學(xué)院，寧夏 銀川 750021)

近年來(lái)，空氣質(zhì)量問(wèn)題成為了全社會(huì)廣泛關(guān)注的焦點(diǎn)問(wèn)題，空氣質(zhì)量狀況影響著人們的日常生活出行和當(dāng)?shù)厣鐣?huì)生產(chǎn)活動(dòng)。因此，空氣質(zhì)量的預(yù)測(cè)對(duì)當(dāng)?shù)卮髿馕廴颈O(jiān)管和治理有著重要的意義。

我國(guó)學(xué)者楊瑞君等[1]將隨機(jī)森林(BF)模型引入上海市空氣質(zhì)量的評(píng)價(jià)中，結(jié)果表明BF模型在空氣質(zhì)量等級(jí)的評(píng)價(jià)中有較高的準(zhǔn)確率；化虎蝶等[2]將貝葉斯網(wǎng)絡(luò)(BN)引入到大連市的空氣質(zhì)量評(píng)價(jià)與預(yù)測(cè)中，并與模糊綜合評(píng)價(jià)法的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比，結(jié)果表明BN法在預(yù)測(cè)精度上有顯著的提高。這些研究雖對(duì)城市空氣質(zhì)量的預(yù)測(cè)提供了可靠、準(zhǔn)確的方法，但未考慮氣象因素對(duì)城市空氣質(zhì)量的影響。楊正理等[3]將影響城市空氣質(zhì)量的34種氣象因素作為輸入，將互信息與BF算法引入太原市空氣質(zhì)量的預(yù)測(cè)中，結(jié)果表明該方法能夠有效識(shí)別影響城市不同區(qū)域空氣質(zhì)量的氣象因素，并在城市空氣質(zhì)量預(yù)測(cè)上有較高的可靠性和準(zhǔn)確率；郭飛等[4]將溫度、氣壓、濕度、風(fēng)速4種氣象因素作為輸入，并將改進(jìn)的支持向量機(jī)(Support Vector Machine,SVM)應(yīng)用于沈陽(yáng)市空氣質(zhì)量預(yù)測(cè)中，結(jié)果表明改進(jìn)后的SVM有效地降低了預(yù)測(cè)錯(cuò)誤率，提升了預(yù)測(cè)精度。這些研究雖然考慮了氣象等因素對(duì)城市空氣質(zhì)量的影響，但未對(duì)氣象等因素進(jìn)行信息的濃縮，當(dāng)因素較多時(shí)，因素之間可能存在信息的重疊和相關(guān)性，從而增加模型結(jié)構(gòu)和數(shù)據(jù)分析的復(fù)雜度[5]。

為了解決上述研究中所存在的問(wèn)題，本文結(jié)合主成分分析(PCA)方法對(duì)數(shù)據(jù)降維以消除信息冗余的特點(diǎn)以及BN在多因素復(fù)雜系統(tǒng)中具有較強(qiáng)推理能力的優(yōu)勢(shì)[6]，以銀川市2015年1月1日至2019年10月31日的空氣質(zhì)量數(shù)據(jù)和同期的氣象數(shù)據(jù)為基礎(chǔ)，將PCA-BN模型引入到銀川市空氣質(zhì)量的預(yù)測(cè)中，建立了基于PCA-BN的銀川市空氣質(zhì)量預(yù)測(cè)模型，并對(duì)模型預(yù)測(cè)準(zhǔn)確率以及可靠性和可行性進(jìn)行了驗(yàn)證。

1 數(shù)據(jù)來(lái)源與研究方法

1. 1 數(shù)據(jù)來(lái)源

本文通過(guò)網(wǎng)絡(luò)爬蟲(chóng)爬取銀川市2015年1月1日至2019年10月31日的空氣質(zhì)量等級(jí)、空氣質(zhì)量指數(shù)以及空氣中細(xì)顆粒物(PM2.5)、可吸入顆粒物(PM10)、二氧化硫(SO2)、一氧化碳(CO)、二氧化氮(NO2)、臭氧(O3)6種主要污染物濃度數(shù)據(jù)和同期的氣象數(shù)據(jù)日均值。其中，同期的氣象數(shù)據(jù)主要包括平均地表氣溫、日最高地表氣溫、日最低地表氣溫、平均相對(duì)濕度、日照時(shí)數(shù)、平均氣溫、日最高氣溫、日最低氣溫、平均氣壓、平均風(fēng)速10項(xiàng)指標(biāo)?？諝庵?種主要污染物濃度數(shù)據(jù)的時(shí)間序列圖見(jiàn)圖1，10項(xiàng)氣象指標(biāo)數(shù)據(jù)的時(shí)間序列圖見(jiàn)圖2。

圖1 2015年1月1日至2019年10月31日銀川市空氣中6種主要污染物濃度數(shù)據(jù)的時(shí)間序列圖(單位：μg/m3)Fig.1 Time series diagram of the content data of six major pollutants in the air of Yinchuan City from Jan 1st, 2015 to Oct 31st,2019(unit：μg/m3)

1. 2 研究方法

本文研究方法分為三個(gè)步驟進(jìn)行：①對(duì)1 765條樣本數(shù)據(jù)進(jìn)行數(shù)據(jù)清洗，處理樣本數(shù)據(jù)中存在的缺失值；②利用主成分分析方法對(duì)氣象數(shù)據(jù)進(jìn)行降維處理并提取氣象因素的綜合評(píng)價(jià)指標(biāo)；③構(gòu)建基于PCA-BN的銀川市空氣質(zhì)量預(yù)測(cè)模型。

1.2.1 數(shù)據(jù)清洗

首先對(duì)1 765條樣本數(shù)據(jù)進(jìn)行數(shù)據(jù)清洗，處理樣本數(shù)據(jù)中存在的缺失值。由于樣本數(shù)據(jù)中存在缺失值的數(shù)據(jù)量占樣本數(shù)據(jù)總量的比例較小，因此對(duì)樣本數(shù)據(jù)中連續(xù)缺失的數(shù)據(jù)做剔除處理，對(duì)于一個(gè)有數(shù)據(jù)缺失的缺失值， 采用相鄰日期前后兩項(xiàng)數(shù)據(jù)的平均值進(jìn)行插補(bǔ)。將經(jīng)過(guò)數(shù)據(jù)清洗后得到的1 677條數(shù)據(jù)作為新的樣本數(shù)據(jù)，以樣本數(shù)據(jù)中的70%作為訓(xùn)練數(shù)據(jù)集，30%作為驗(yàn)證數(shù)據(jù)集，構(gòu)建PCA-BN模型。

圖2 2015年1月1日至2019年10月31日銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)的時(shí)間序列圖Fig.2 Time series chart of the content data of ten meteorological indicators of Yinchuan City from Jan 1st, 2015 to Oct 31st,2019

1.2.2 主成分分析

主成分分析(Principal Component Analysis,PCA)方法常用于處理變量之間具有高度相關(guān)性的數(shù)據(jù)[7]。該方法是運(yùn)用降維的思想去除冗余的信息，將多個(gè)變量轉(zhuǎn)化為由多個(gè)變量線(xiàn)性表示的少數(shù)主成分的過(guò)程[8]。根據(jù)PCA原理，可分為如下幾個(gè)步驟進(jìn)行計(jì)算[9]：

(1) 對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理：

(1)

(2) 計(jì)算相關(guān)系數(shù)矩陣R=(rjk)m×m：

(2)

(3) 求解相關(guān)系數(shù)矩陣的特征值與特征向量。

(4) 提取主成分：主要根據(jù)特征根和累計(jì)貢獻(xiàn)率的值來(lái)提取主成分，一般選取特征根大于1的主成分[10]。主成分zt的方差貢獻(xiàn)率wt為

(3)

wt值越大，表明主成分Zt對(duì)原始變量信息的解釋能力越強(qiáng)。前q個(gè)主成分的累計(jì)方差貢獻(xiàn)率Wq為

(4)

Wq表示所提取的q個(gè)主成分能夠解釋原始變量信息量的多少。

1.2.3 貝葉斯網(wǎng)絡(luò)

貝葉斯網(wǎng)絡(luò)(Bayesian Networks,BN)，又稱(chēng)信度網(wǎng)絡(luò)，是一種結(jié)合了概率論和圖論的不確定性知識(shí)表達(dá)與推理的方法。BN作為一種圖形化模型，能夠表達(dá)出因果關(guān)系和相關(guān)的不確定性知識(shí)，該模型由節(jié)點(diǎn)、有向弧和條件概率表(Conditional Probability Tables,CPT)組成[11-12]。由于BN要求各節(jié)點(diǎn)之間不能形成閉環(huán)，因此BN符合有向無(wú)環(huán)圖(Directed Acyclic Graph,DAG)結(jié)構(gòu)[13]。BN因其具有的簡(jiǎn)潔、直觀、準(zhǔn)確、方便的特點(diǎn)，在實(shí)際工程中得到了廣泛的應(yīng)用[14]。在BN中，貝葉斯公式為概率計(jì)算的基礎(chǔ)，貝葉斯公式[15]如下：

(5)

其中P(Xi)>0，P(E|Xi)>0，P(Xi|E)>0。

式中：X1，X2，…，Xn為樣本空間Ω中兩兩互斥的事件;P(Xi)為事件Xi的先驗(yàn)概率[16]；P(Xi|E)為后驗(yàn)概率，表示在已知事件E發(fā)生的概率下，事件Xi發(fā)生的概率；P(E|Xi)為條件概率，表示在已知事件Xi發(fā)生的概率下，事件E發(fā)生的概率。

BN的學(xué)習(xí)分為BN的結(jié)構(gòu)學(xué)習(xí)和BN的參數(shù)學(xué)習(xí)。其中，BN結(jié)構(gòu)學(xué)習(xí)是根據(jù)數(shù)據(jù)經(jīng)驗(yàn)和先驗(yàn)知識(shí)確定BN結(jié)構(gòu)，在得到確定的BN結(jié)構(gòu)的基礎(chǔ)上，采用合理的參數(shù)學(xué)習(xí)方法確定BN拓?fù)浣Y(jié)構(gòu)中各節(jié)點(diǎn)之間的條件概率[17]；BN參數(shù)學(xué)習(xí)是以BN結(jié)構(gòu)為基礎(chǔ)，如圖3所示為一個(gè)簡(jiǎn)單結(jié)構(gòu)的BN結(jié)構(gòu)圖，根據(jù)先驗(yàn)概率、后驗(yàn)概率以及條件概率等知識(shí)對(duì)構(gòu)成BN結(jié)構(gòu)的各節(jié)點(diǎn)進(jìn)行定量分析，并得到各節(jié)點(diǎn)間相關(guān)關(guān)系的強(qiáng)弱[18]。

圖3 簡(jiǎn)單結(jié)構(gòu)的貝葉斯網(wǎng)絡(luò)(BN)結(jié)構(gòu)圖Fig.3 Diagram of Bayesian Network (BN) with simple structure

如圖3所示，節(jié)點(diǎn)X1、X2均通過(guò)有向線(xiàn)段指向節(jié)點(diǎn)Y，有向線(xiàn)的起點(diǎn)為子節(jié)點(diǎn)，終點(diǎn)指向父節(jié)點(diǎn)，利用參數(shù)學(xué)習(xí)方法可分別求得父節(jié)點(diǎn)Y與各子節(jié)點(diǎn)X1、X2的條件概率。

為了提高BN局部因果結(jié)構(gòu)的效率并有效地對(duì)原始變量的特征空間進(jìn)行降維，本文利用馬爾科夫毯(Markov Blanket,MB)算法進(jìn)行BN結(jié)構(gòu)學(xué)習(xí)。馬爾科夫毯的定義為：集合U為隨機(jī)變量的全集，對(duì)于一個(gè)給定的變量X(X∈U)和變量集MB(MB∈U且X?MB)，若滿(mǎn)足[19]：

X⊥{U-MB-{X}}|MB

則稱(chēng)最小變量集MB為變量X的馬爾科夫毯。

本文所采集的數(shù)據(jù)為各項(xiàng)指標(biāo)的日均值且數(shù)據(jù)量充足、較為完整，故采用最大似然法作為BN結(jié)構(gòu)的參數(shù)學(xué)習(xí)方法。假定數(shù)據(jù)集C={C1，C2，…，Cn}，節(jié)點(diǎn)Xi共有hi個(gè)取值，且節(jié)點(diǎn)Xi的參數(shù)為θi，節(jié)點(diǎn)Xi的父節(jié)點(diǎn)α(Xi)的取值共有qi個(gè)組合，則節(jié)點(diǎn)Xi的參數(shù)θi的對(duì)數(shù)似然函數(shù)表達(dá)式如下[20]：

(6)

(7)

式中：mijk為數(shù)據(jù)集C中滿(mǎn)足Xi=k且α(Xi)=j的樣本數(shù)量。

2 研究結(jié)果與分析

2.1 PCA結(jié)果分析

本文利用PCA方法對(duì)銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)進(jìn)行降維處理，得到KMO和Bartlett的球形度檢驗(yàn)結(jié)果見(jiàn)表1。

表1 KMO和Bartlett的球形度檢驗(yàn)結(jié)果

由表1可知，KMO的檢驗(yàn)結(jié)果為0.790，大于0.500；Bartlett的球形度檢驗(yàn)拒絕單位相關(guān)陣的原假設(shè)(Sig.<0.001),說(shuō)明本文選取的10項(xiàng)氣象指標(biāo)數(shù)據(jù)適合進(jìn)行主成分分析。

銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)經(jīng)過(guò)PCA法降維處理后提取的3個(gè)主成分的方差百分比和累計(jì)方差貢獻(xiàn)率，見(jiàn)表2。

表2 銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)經(jīng)過(guò)PCA法降維處理后提取的3個(gè)主成分的方差百分比和累計(jì)方差貢獻(xiàn)率

由表2可知，所提取的3個(gè)主成分解釋了原10項(xiàng)氣象指標(biāo)數(shù)據(jù)信息量的91.563%，大于85%，說(shuō)明提取的3個(gè)主成分是合適的。

銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)經(jīng)過(guò)主成分分析后的碎石圖，見(jiàn)圖4。

圖4 銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)經(jīng)過(guò)主成分分析后的 碎石圖Fig.4 Scree plot in PCA analysis of ten meteorological indicators of Yinchuan City

由圖4可見(jiàn)，第三個(gè)主成分之后，曲線(xiàn)變得平滑，沒(méi)有出現(xiàn)明顯的變化趨勢(shì)，這從側(cè)面說(shuō)明PCA法將銀川市10項(xiàng)氣象指標(biāo)數(shù)據(jù)降維處理后提取3個(gè)主成分是合適的。

由主成分因子負(fù)荷矩陣(見(jiàn)表3)可知，第一主成分主要包含原指標(biāo)中各氣溫和平均氣壓指標(biāo)，因此第一主成分可作為總體氣溫狀況的綜合描述指標(biāo)；第二主成分主要包含原指標(biāo)中平均相對(duì)濕度和日照時(shí)數(shù)指標(biāo)，且平均相對(duì)濕度指標(biāo)相較于日照時(shí)數(shù)指標(biāo)在第二個(gè)主成分上具有更大的載荷，說(shuō)明在第二主成分中平均相對(duì)濕度相比于日照時(shí)數(shù)具有更高的重要性，因此第二主成分可作為濕度狀況的綜合描述指標(biāo)；第三主成分主要包含原指標(biāo)中平均風(fēng)速指標(biāo)，因此第三主成分可作為風(fēng)速狀況的綜合描述指標(biāo)。

由表2和表3可知，第一主成分方差百分比較高，方差貢獻(xiàn)率達(dá)到65.014%，說(shuō)明氣象因素中氣溫對(duì)空氣質(zhì)量有較大的影響。銀川市10項(xiàng)氣象指標(biāo)的3個(gè)主成分序列圖，見(jiàn)圖5。

表3 主成分因子負(fù)荷矩陣

圖5 銀川市10項(xiàng)氣象指標(biāo)的3個(gè)主成分的序列圖Fig.5 Sequence diagram of three principal components extracted from the ten meteorological indicators of Yinchuan City

2.2 構(gòu)建PCA-BN模型

銀川市原10項(xiàng)氣象指標(biāo)經(jīng)過(guò)主成分分析后提取得到3項(xiàng)氣象綜合描述指標(biāo)，不僅降低了氣象指標(biāo)的維數(shù)，并且反映了原氣象指標(biāo)的絕大部分信息量。以3項(xiàng)氣象綜合描述指標(biāo)和大氣中6項(xiàng)主要污染物濃度指標(biāo)共9項(xiàng)指標(biāo)數(shù)據(jù)作為輸入，并對(duì)指標(biāo)數(shù)據(jù)進(jìn)行離散化處理,構(gòu)建基于PCA-BN的銀川市空氣質(zhì)量預(yù)測(cè)模型，即PCA-BN模型。

構(gòu)建PCA-BN模型時(shí)，首先利用MB算法進(jìn)行BN的結(jié)構(gòu)學(xué)習(xí)。以質(zhì)量等級(jí)為目標(biāo)變量，先得到目標(biāo)變量的MB;再根據(jù)目標(biāo)變量的MB進(jìn)一步確定所研究的預(yù)測(cè)變量的條件概率分布，由于目標(biāo)變量的MB包含了目標(biāo)變量的所有信息，則MB以外的變量可被視為冗余變量并排除了其余變量的影響[21]；最后基于條件獨(dú)立性測(cè)試判斷兩變量間是否相互獨(dú)立，若相互獨(dú)立，則刪去連接兩變量節(jié)點(diǎn)的線(xiàn)。

然后，利用最大似然估計(jì)進(jìn)行BN的參數(shù)學(xué)習(xí)，確定每個(gè)變量的條件概率分布。PCA-BN模型中輸入變量的重要性見(jiàn)表4，PCA-BN模型的網(wǎng)絡(luò)結(jié)構(gòu)圖見(jiàn)圖6，其中輸入變量節(jié)點(diǎn)的顏色深淺表示預(yù)測(cè)變量重要性的高低。

表4 PCA-BN模型中輸入變量的重要性

圖6 PCA-BN模型的網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.6 Network structure diagram of the PCA-BN model

由表4和圖6可知，PCA-BN模型中所輸入的9項(xiàng)指標(biāo)均會(huì)對(duì)銀川市空氣質(zhì)量等級(jí)產(chǎn)生影響，且各項(xiàng)指標(biāo)間也會(huì)有相互的影響關(guān)系。其中，SO2對(duì)于銀川市空氣質(zhì)量的影響較大，故必須注意對(duì)于SO2污染源頭的控制。

PCA-BN模型預(yù)測(cè)結(jié)果的準(zhǔn)確率見(jiàn)表5。

表5 PCA-BN模型預(yù)測(cè)結(jié)果的準(zhǔn)確率

由表5可知，PCA-BN模型在訓(xùn)練數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集中的正確率接近，且預(yù)測(cè)結(jié)果的準(zhǔn)確率達(dá)到97.88%。

PCA-BN模型預(yù)測(cè)結(jié)果的混淆矩陣見(jiàn)表6，其中行為空氣質(zhì)量實(shí)際值，列為空氣質(zhì)量預(yù)測(cè)值。

由表6可知，在驗(yàn)證數(shù)據(jù)集中，空氣質(zhì)量實(shí)際值為良的情況下，PCA-BN模型的預(yù)測(cè)值為輕度污染的個(gè)數(shù)為1個(gè)；空氣質(zhì)量實(shí)際值為輕度污染的情況下，PCA-BN模型的預(yù)測(cè)值為良的個(gè)數(shù)為2個(gè)，為中度污染的個(gè)數(shù)為3個(gè)；空氣質(zhì)量實(shí)際值為中度污染的情況下，PCA-BN模型的預(yù)測(cè)值為輕度污染的個(gè)數(shù)為4個(gè)，為重度污染的個(gè)數(shù)為1個(gè)；其余情況下均未出現(xiàn)預(yù)測(cè)錯(cuò)誤。

表6 PCA-BN模型預(yù)測(cè)結(jié)果的混淆矩陣

根據(jù)上述預(yù)測(cè)結(jié)果的準(zhǔn)確率和混淆矩陣，說(shuō)明PCA-BN模型在銀川市空氣質(zhì)量的預(yù)測(cè)中有著較高的預(yù)測(cè)精度，具有一定的應(yīng)用價(jià)值。

2.3 PCA-BN模型與隨機(jī)森林模型預(yù)測(cè)結(jié)果的比較

為了進(jìn)一步說(shuō)明PCA-BN模型在銀川市空氣質(zhì)量預(yù)測(cè)中的可靠性和可行性。本文利用隨機(jī)森林(Random Forests,RF)模型并將空氣中6種大氣主要污染物濃度數(shù)據(jù)和同期的氣象數(shù)據(jù)共16項(xiàng)指標(biāo)作為輸入，得到基于RF的銀川市空氣質(zhì)量預(yù)測(cè)結(jié)果，并與基于PCA-BN模型的銀川市空氣質(zhì)量預(yù)測(cè)結(jié)果進(jìn)行了比較，其結(jié)果見(jiàn)表7。

表7 PCA-BN模型與RF模型預(yù)測(cè)結(jié)果的對(duì)比

由表7可知，RF模型和PCA-BN模型預(yù)測(cè)結(jié)果的準(zhǔn)確率分別為88.03%和97.88%，PCA-BN模型對(duì)銀川市空氣質(zhì)量預(yù)測(cè)結(jié)果的準(zhǔn)確率有明顯的提高，說(shuō)明該模型在空氣質(zhì)量的預(yù)測(cè)中有一定的可行性和可靠性。

3 結(jié) 語(yǔ)

本文將PCA-BN模型引入銀川市空氣質(zhì)量的預(yù)測(cè)中，結(jié)果表明：該方法預(yù)測(cè)結(jié)果的準(zhǔn)確率達(dá)到了97.88%，并且有效降低了數(shù)據(jù)分析的復(fù)雜性。今后針對(duì)城市空氣質(zhì)量的預(yù)測(cè)還需要考慮更多的氣象條件、區(qū)域環(huán)境、污染源排放、城市結(jié)構(gòu)、城市發(fā)展水平等有可能影響到城市空氣質(zhì)量的因素，并提出合理的建議，以達(dá)到對(duì)空氣質(zhì)量狀況更加準(zhǔn)確的評(píng)價(jià)，促進(jìn)城市空氣質(zhì)量的提升。