郭映映,齊賀香,李素文,牟福生
(淮北師范大學(xué)物理與電子信息學(xué)院,安徽 淮北 235000)
隨著我國(guó)從制造業(yè)大國(guó)向制造業(yè)強(qiáng)國(guó)轉(zhuǎn)變,能源消耗與污染問(wèn)題變得愈發(fā)嚴(yán)重,空氣污染問(wèn)題受到越來(lái)越多的關(guān)注。其中,NO2在酸雨、臭氧和光化學(xué)煙霧的形成過(guò)程中都扮演著重要角色[1],不僅影響著生態(tài)平衡與綠色環(huán)境,更影響著社會(huì)發(fā)展與人類(lèi)安全,因此NO2空氣污染問(wèn)題是亟需解決的重要難題之一。
近年來(lái),關(guān)于PM2.5、O3、CO2等污染物建模預(yù)測(cè)的研究較多。國(guó)內(nèi)彭巖等[2]構(gòu)建了集成樹(shù)-梯度決策樹(shù)的模型對(duì)PM2.5濃度進(jìn)行了預(yù)測(cè);張怡文等[3]利用PCA-BP神經(jīng)網(wǎng)絡(luò)也進(jìn)行了PM2.5濃度的預(yù)測(cè);黃鴻等[4]利用深度信念網(wǎng)絡(luò)和極限學(xué)習(xí)機(jī)對(duì)SO2濃度進(jìn)行檢測(cè);李素文等[5]基于改進(jìn)Elman網(wǎng)絡(luò)對(duì)污染物濃度進(jìn)行實(shí)時(shí)預(yù)測(cè);Zeng等[6]預(yù)測(cè)車(chē)輛排放的CO2及其在生態(tài)路線導(dǎo)航中的應(yīng)用。國(guó)外Luna等[7]人利用神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)預(yù)測(cè)出巴西里約熱內(nèi)盧對(duì)流層O3濃度。但關(guān)于NO2的研究較少,因此本文以合肥地區(qū)NO2氣體濃度為研究對(duì)象,利用反向傳播(BP)神經(jīng)網(wǎng)絡(luò)模型來(lái)預(yù)測(cè)NO2濃度。粒子群算法(PSO)是一種通過(guò)種群行為找出目標(biāo)值最優(yōu)解的方法,具有全局搜索、收斂速度快等優(yōu)點(diǎn)[8]。BP神經(jīng)網(wǎng)絡(luò)局部搜索能力較強(qiáng),但易受初始權(quán)值的影響,產(chǎn)生局部最優(yōu)[9]。因此將二者結(jié)合起來(lái)構(gòu)建粒子群優(yōu)化的反向傳播(PSO-BP)網(wǎng)絡(luò)預(yù)測(cè)模型,先用PSO算法對(duì)BP網(wǎng)絡(luò)的權(quán)重和閾值進(jìn)行初始化,范圍縮小后,再用BP網(wǎng)絡(luò)精確求解,從而實(shí)現(xiàn)NO2氣體濃度的準(zhǔn)確預(yù)測(cè)。
BP神經(jīng)網(wǎng)絡(luò)包括輸入層、隱含層和輸出層三部分(圖1),輸入層用來(lái)接收外界信息,隱含層用來(lái)對(duì)信息進(jìn)行表示和保存,輸出層用來(lái)對(duì)輸入層的信息進(jìn)行識(shí)別和判斷[10]。各層神經(jīng)元之間不存在連接,但與相鄰層的神經(jīng)元之間存在連接。學(xué)習(xí)的過(guò)程分為信號(hào)的正向傳輸和誤差信號(hào)的反向傳輸兩部分。信號(hào)正向傳輸?shù)倪^(guò)程中,信號(hào)從輸入層經(jīng)隱含層傳輸并處理后,最終到達(dá)輸出層。當(dāng)輸出的預(yù)測(cè)值與真實(shí)值之間的誤差不滿足要求時(shí),則誤差信號(hào)開(kāi)始反向傳輸,誤差信號(hào)由輸出層向隱含層傳輸后重新回到輸入層,各層神經(jīng)元得到的誤差信號(hào),被用來(lái)調(diào)整該層神經(jīng)元的權(quán)值和閾值。信號(hào)的正向傳輸和誤差信號(hào)的反向傳輸不斷循環(huán),因此權(quán)值和閾值被不斷調(diào)整,直至網(wǎng)絡(luò)輸出的誤差滿足要求或?qū)W習(xí)次數(shù)達(dá)到設(shè)定時(shí),結(jié)束學(xué)習(xí)[11]。
圖1 神經(jīng)網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)圖Fig.1 Topological structure of BP neural network
BP神經(jīng)網(wǎng)絡(luò)可以看成自變量是輸入值、因變量是網(wǎng)絡(luò)模型預(yù)測(cè)值的非線性函數(shù)。若網(wǎng)絡(luò)的輸入節(jié)點(diǎn)有n個(gè),輸出節(jié)點(diǎn)有m個(gè),則BP神經(jīng)網(wǎng)絡(luò)就可以看成是一個(gè)包含n個(gè)自變量和m個(gè)因變量的非線性函數(shù)。網(wǎng)絡(luò)的學(xué)習(xí)過(guò)程如下:
1)初始化網(wǎng)絡(luò)的訓(xùn)練次數(shù)、誤差下限、學(xué)習(xí)速率及連接權(quán)值和閾值。
2)通過(guò)連接權(quán)值ωi,j,輸入信號(hào)xi及閾值aj,得到隱含層的輸出結(jié)果:
式中 f為隱含層的激勵(lì)函數(shù),l為隱含層的節(jié)點(diǎn)個(gè)數(shù)。
3)利用輸出層的連接權(quán)值ωj,k和閾值bk及2)中隱含層的輸出結(jié)果,得到輸出層的輸出結(jié)果:
4)利用網(wǎng)絡(luò)模型的輸出結(jié)果Yk和期望輸出結(jié)果Vk,得到網(wǎng)絡(luò)模型的誤差:
5)通過(guò)網(wǎng)絡(luò)模型的誤差ek對(duì)網(wǎng)絡(luò)的連接權(quán)值ωi,j,ωj,k和閾值aj、bk進(jìn)行更新,即
6)判斷是否滿足終止條件,若不滿足條件,則返回2),若滿足條件,則結(jié)束學(xué)習(xí)。
粒子群算法是一種新型的進(jìn)化算法,它來(lái)源于對(duì)鳥(niǎo)類(lèi)捕食行為的研究[12]。其核心思想是利用搜索空間中的每一個(gè)粒子來(lái)代表問(wèn)題的每一個(gè)候選解,由目標(biāo)函數(shù)來(lái)確定每個(gè)粒子的適應(yīng)度值,粒子下一步運(yùn)動(dòng)的位置和速度由該粒子本身和已經(jīng)確定好適應(yīng)度值的粒子決定,最后找出個(gè)體粒子的最優(yōu)解和種群的最優(yōu)解。在粒子中不斷進(jìn)行尋優(yōu),直到滿足迭代次數(shù),結(jié)束尋優(yōu)。假設(shè)搜索空間為D,n個(gè)不同的粒子構(gòu)成種群X=(X1,X2,...,Xn),用Xi=(xi,1,xi,2,...,xi,D)T表示第i個(gè)粒子在這個(gè)D維搜索空間中的位置,即每個(gè)問(wèn)題的候選解,將Xi代入到適應(yīng)函數(shù) f(xi)中,求適應(yīng)值。Vi=(vi,1,vi,2,...,vi,D)T表示第i個(gè)粒子的速度,Pbesti,d=(pi,1,pi,2,...,pi,D)T表示個(gè)體的極值,Gbestd=(g1,g2,...,gD)T表示群體的極值。在粒子每一次的搜索過(guò)程中,根據(jù)種群的群體極值和個(gè)體極值來(lái)更新粒子的位置和速度,更新公式為
式中Vi,d為粒子速度,k為迭代次數(shù),i=1,2,...,n,d=1,2,...,D,c1和c2為學(xué)習(xí)因子,r1和r1為[0,1]之間的隨機(jī)數(shù),xi,d為粒子位置。
在粒子群算法中存在著粒子向鄰域或群體歷史最優(yōu)位置和粒子自身最優(yōu)位置靠近的現(xiàn)象,從而導(dǎo)致粒子的快速趨同效應(yīng),引起種群產(chǎn)生過(guò)早收斂、局部極值等不良現(xiàn)象。根據(jù)遺傳算法中的變異思想,在粒子群算法中加入自適應(yīng)變異算子,以一定的可能性選取變量并重新對(duì)其進(jìn)行初始化,來(lái)避免上述不良現(xiàn)象的產(chǎn)生。變異的操作使粒子在迭代過(guò)程中通過(guò)種群搜索空間的不斷減小而跳出當(dāng)前的局部搜索空間,在更大的空間中進(jìn)行搜索,這不僅保證了粒子的多樣性不會(huì)減少,還確保粒子不會(huì)產(chǎn)生局部最優(yōu)的現(xiàn)象。因此,BP神經(jīng)網(wǎng)絡(luò)連接權(quán)值和閾值經(jīng)粒子群算法對(duì)其進(jìn)行初始化后,再利用BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型進(jìn)行訓(xùn)練,最后輸出的NO2濃度就是預(yù)測(cè)的最優(yōu)解。粒子群算法改進(jìn)BP神經(jīng)網(wǎng)絡(luò)模型的流程圖如圖2。
圖2 粒子群算法改進(jìn)BP神經(jīng)網(wǎng)絡(luò)算法的流程圖Fig.2 Flow diagram of the improved BP neural network algorithm by particle swarm optimization
參與空氣質(zhì)量評(píng)價(jià)的主要污染物為NO2、PM2.5、PM10、SO2、CO和O3,因此下面以合肥市2017年1月1日至2019年12月31日的污染物PM2.5、PM10、SO2、CO和O3濃度為影響因子,對(duì)NO2濃度進(jìn)行預(yù)測(cè),以制定有效的環(huán)境改善措施。實(shí)驗(yàn)數(shù)據(jù)來(lái)源共分兩類(lèi),一類(lèi)是來(lái)源于中國(guó)空氣質(zhì)量在線監(jiān)測(cè)分析平臺(tái)的大氣污染數(shù)據(jù),該數(shù)據(jù)收集了研究期間每日的NO2、PM2.5、PM10、SO2、CO平均濃度和O3_8h(臭氧的日最大值8 h滑動(dòng)平均濃度)。另一類(lèi)是氣象數(shù)據(jù),來(lái)源于中國(guó)氣象科學(xué)數(shù)據(jù)中心,選取同期與NO2濃度相關(guān)性較大的氣象數(shù)據(jù),即每日平均溫度、平均相對(duì)濕度、平均風(fēng)速、平均氣壓和總降水量。合肥市2017–2019年樣本個(gè)案數(shù)、缺失計(jì)數(shù)及缺失百分比如表1所示,空氣質(zhì)量與氣象要素描述性統(tǒng)計(jì)結(jié)果如表2所示。
表1 合肥市2017–2019年樣本個(gè)案數(shù)、缺失計(jì)數(shù)及缺失百分比統(tǒng)計(jì)Table 1 Statistics of sample cases,missing count and missing percentage in Hefei City from 2017 to 2019
表2 合肥市2017–2019年空氣質(zhì)量與氣象要素描述性統(tǒng)計(jì)Table 2 Descriptive statistics of air quality and meteorological elements in Hefei City from 2017 to 2019
從表1中可以看出,本次實(shí)驗(yàn)選取觀測(cè)數(shù)據(jù)共有1095組,其中有27組存在缺失值,單個(gè)變量中缺失個(gè)數(shù)最多的是O3_8h,僅占其樣本數(shù)的1.1%。利用插值法中的分段三次樣條插值[13]對(duì)缺值進(jìn)行補(bǔ)充,補(bǔ)充后的各要素的樣本個(gè)案數(shù)均為1095,插值后的完整數(shù)據(jù)描述性統(tǒng)計(jì)特征見(jiàn)表3。
表3 插值后的完整數(shù)據(jù)描述性統(tǒng)計(jì)Table 3 Descriptive statistics of complete data after interpolation
圖3(a)展示了合肥地區(qū)2017–2019年NO2濃度的日均值,可以看出,合肥地區(qū)NO2日均濃度在一年內(nèi)呈現(xiàn)“U”型變化,即冬季濃度最高,夏季濃度最低的特征。三年來(lái)最高NO2濃度出現(xiàn)在2017年3月10日,為137μg·m-3,最低NO2濃度出現(xiàn)在2018年8月17日,僅為9μg·m-3。圖3(b)展示了合肥地區(qū)NO2濃度2017–2019年年度每周工作日和非工作日的日變化平均值及其在三年期間總平均值,可以看出,NO2濃度穩(wěn)定在總體年均值附近,與具體日期關(guān)系不大,沒(méi)有明顯的“周末效應(yīng)”[14]。
圖3 2017–2019年間合肥地區(qū)NO2濃度日均值(a)和年度每周工作日和非工作日的日變化平均值及其在三年期間總平均值(b)Fig.3 Daily mean of NO2concentration in Hefei from 2017 to 2019(a)and the daily variation mean value of working and non working days per week in each year and its total mean value over the three-year period(b)
圖4給出了三年來(lái)合肥地區(qū)NO2日均濃度的月均值變化圖,發(fā)現(xiàn)這三年變化趨勢(shì)一致,均呈現(xiàn)“U”型曲線分布特征。其中,11月份平均濃度最高,為60.47μg·m-3,隨后再逐漸降低,到7月份平均濃度達(dá)到最低,為27.19μg·m-3,隨后再逐漸升高。觀察2019年NO2月均濃度變化,發(fā)現(xiàn)2月份下降迅速,這是由于2月份是疫情管控期間,交通和工廠管制使得NO2排放量減少。3月份開(kāi)始,疫情得到初步控制,工廠復(fù)工、交通恢復(fù)導(dǎo)致NO2濃度又開(kāi)始增加。
圖4 2017–2019年合肥地區(qū)NO2濃度的月均值Fig.4 Monthly mean NO2concentration in Hefei from 2017 to 2019
為分析各影響因子的變動(dòng)對(duì)NO2濃度變化的影響程度,選取2017年1月1日至2019年12月31日合肥地區(qū)6項(xiàng)大氣污染物和5項(xiàng)氣象因素的觀測(cè)數(shù)據(jù),對(duì)NO2日均濃度構(gòu)建逐步回歸模型[15]。進(jìn)行變量篩選時(shí),被模型選入或剔除的顯著性水平分別為0.05和0.1,模型篩選過(guò)程的匯總見(jiàn)表4,因變量為NO2日濃度。
表4 逐步回歸模型篩選過(guò)程匯總Table 4 Summary of the screening process of stepwise regression model
模型最先引進(jìn)自變量PM2.5、PM10,緊接著依次引入O3_8h、SO2、CO、平均氣壓、平均2 min風(fēng)速和平均氣溫。模型3剔除自變量PM2.5,模型4重新引入PM2.5,模型7剔除平均氣壓,至此逐步回歸分析建模結(jié)束。模型7為逐步回歸分析得到的最終模型,其對(duì)應(yīng)的回歸系數(shù)和統(tǒng)計(jì)特征見(jiàn)表5。該模型的表達(dá)式為
表5 模型7對(duì)應(yīng)的回歸系數(shù)和t統(tǒng)計(jì)量及其概率Table 5 Regression coefficients,t statistics and its probabilities corresponding to Model 7
式中Y表示NO2日均濃度,X1、X2、X3、X4分別表示PM2.5、PM10、SO2、CO日均濃度,X5、X6分別表示平均2 min風(fēng)速和平均氣溫。該逐步回歸模型可以解釋為PM2.5等自變量的變動(dòng)會(huì)引起因變量NO2發(fā)生相應(yīng)的變動(dòng)。例如,當(dāng)其他條件不變時(shí),PM2.5的濃度每增加1單位,NO2的濃度就會(huì)減少0.2218單位。由表4可知,該逐步回歸模型F-stat為414.555,其對(duì)應(yīng)概率為0.0,表明在0.05顯著水平下,該模型通過(guò)了顯著性檢驗(yàn),并且模型調(diào)整后的R2為0.694,說(shuō)明該模型具有較好的擬合效果。此外模型t統(tǒng)計(jì)量的概率均小于0.05顯著性水平,表明這些變量均與NO2日均值濃度的相關(guān)性較大。
首先對(duì)采集到的數(shù)據(jù)進(jìn)行預(yù)處理操作,通過(guò)逐步回歸方法對(duì)預(yù)處理后的數(shù)據(jù)進(jìn)行特征篩選,篩選出與NO2濃度值關(guān)聯(lián)度較高的因子,將這些因子構(gòu)成的數(shù)據(jù)集作為BP神經(jīng)網(wǎng)絡(luò)的輸入樣本,分別利用遺傳算法和粒子群算法對(duì)BP神經(jīng)網(wǎng)絡(luò)的權(quán)值、閾值進(jìn)行調(diào)整,再將這三種預(yù)測(cè)模型的結(jié)果進(jìn)行對(duì)比。利用合肥地區(qū)2017年1月1日至2019年12月31日的數(shù)據(jù)構(gòu)建預(yù)測(cè)模型。其中2017年1月1日至2019年9月22日共995天的NO2日均濃度數(shù)據(jù)做為訓(xùn)練集樣本,2019年9月23日至12月31日共100天的數(shù)據(jù)做為測(cè)試集樣本。由于預(yù)測(cè)是對(duì)未知數(shù)據(jù)進(jìn)行估計(jì),所以預(yù)測(cè)結(jié)果與真實(shí)數(shù)據(jù)之間必定存在誤差。因此引入平均絕對(duì)誤差EMA、平均絕對(duì)百分比誤差EMAP、均方誤差EMS和均方根誤差ERMS,對(duì)模型的預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)。各模型的參數(shù)設(shè)置為:
式中yj為模型的預(yù)測(cè)值,vj為模型的真實(shí)值。
1)BP神經(jīng)網(wǎng)絡(luò)模型參數(shù)設(shè)置:采用5個(gè)隱層節(jié)點(diǎn),傳輸函數(shù)隱層采用默認(rèn)的tansig函數(shù),輸出層采用purelin函數(shù),迭代次數(shù)100次,最小誤差限0.001,學(xué)習(xí)率0.01,歸一化算法為y=(x-xmin)/(xmax-xmin)。
2)遺傳算法改進(jìn)的反向傳播(GA-BP)神經(jīng)網(wǎng)絡(luò)模型參數(shù)設(shè)置:遺傳算法進(jìn)化30代,種群規(guī)模20,交叉選擇概率0.3,變異選擇概率0.1,其他參數(shù)設(shè)置與1)相同。
3)PSO-BP神經(jīng)網(wǎng)絡(luò)模型參數(shù)設(shè)置:粒子群進(jìn)化30代,種群規(guī)模20,學(xué)習(xí)因子C1=C2=1.5,其他參數(shù)設(shè)置與1)相同。
表6給出在三種不同輸入樣本的情況下,利用三種不同的模型對(duì)NO2濃度進(jìn)行預(yù)測(cè)的誤差。三種不同輸入樣本的情況分別為NO2濃度、NO2濃度及全部影響因子和NO2濃度及相關(guān)性較高的影響因子。
表6 三種預(yù)測(cè)模型的預(yù)測(cè)誤差對(duì)比Table 6 Comparison of prediction errors of the three prediction models
對(duì)比表6中三種不同輸入樣本的結(jié)果,發(fā)現(xiàn)NO2濃度及相關(guān)性較高的影響因子作為三種預(yù)測(cè)模型的輸入樣本時(shí),EMA、EMAP、EMS和ERMS的值均最小。再對(duì)比不同模型的結(jié)果,發(fā)現(xiàn)PSO-BP模型輸出結(jié)果的EMA、EMAP、EMS和ERMS的值最小,GA-BP模型輸出結(jié)果次之,BP模型輸出結(jié)果最大。表明NO2及相關(guān)性較高的影響因子作為PSO-BP模型的輸入樣本時(shí)的誤差最小,EMA、EMAP、EMS和ERMS分別為0.039、0.072、0.002、0.046。
圖5給出了NO2濃度及相關(guān)性較高的影響因子作為輸入樣本時(shí)不同預(yù)測(cè)模型的預(yù)測(cè)值與真實(shí)值的對(duì)比結(jié)果。由圖可知,三種預(yù)測(cè)模型均能預(yù)測(cè)出NO2日均濃度的變化趨勢(shì),BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型的預(yù)測(cè)結(jié)果普遍高于真實(shí)值,GA-BP模型的預(yù)測(cè)結(jié)果略低于真實(shí)值,只有PSO-BP模型的預(yù)測(cè)結(jié)果與真實(shí)值的擬合程度、擬合曲線的吻合度都比較高。對(duì)比三種不同的預(yù)測(cè)模型,發(fā)現(xiàn)PSO-BP模型預(yù)測(cè)的結(jié)果明顯優(yōu)于BP模型和GA-BP模型。表明PSO-BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)精度高、誤差小,能夠廣泛應(yīng)用于大氣NO2濃度的實(shí)際預(yù)測(cè)中。
圖5 BP(a)、GA-BP(b)和PSO-BP(c)三種預(yù)測(cè)模型的預(yù)測(cè)值與真實(shí)值的對(duì)比圖Fig.5 Comparison between the real value and the predicted value of BP(a),GA-BP(b)and PSO-BP(c)respectively
以合肥地區(qū)2017年1月1日至2019年12月31日的大氣污染數(shù)據(jù)和氣象數(shù)據(jù)為基礎(chǔ),利用逐步回歸的方法篩選出6個(gè)與NO2濃度相關(guān)性較大的影響因子作為輸入樣本。構(gòu)建了PSO-BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型,與BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型和GA-BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型相比,PSO-BP模型具有全局搜索、收斂速度快、預(yù)測(cè)精度高、算法簡(jiǎn)單等特點(diǎn)。將PSO-BP模型的預(yù)測(cè)結(jié)果與其他兩種模型的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比,發(fā)現(xiàn)PSO-BP模型的擬合效果最好,預(yù)測(cè)值與真實(shí)值之間的均方根誤差最小,僅為0.046。表明PSO-BP模型能夠較準(zhǔn)確地預(yù)測(cè)出NO2濃度的動(dòng)態(tài)變化規(guī)律,可為治理空氣污染問(wèn)題提供基本依據(jù),具有一定的實(shí)際應(yīng)用價(jià)值。