董紅召,廖世凱,楊 強,應 方 (.浙江工業(yè)大學,智能交通系統(tǒng)研究所,浙江 杭州 004;.杭州環(huán)研科技有限公司,浙江 杭州 ;.浙江省杭州生態(tài)環(huán)境監(jiān)測中心,浙江 杭州 0004)
研究表明,工業(yè)生產和廢氣排放與區(qū)域PM2.5濃度高度正相關[1].從大氣環(huán)境容量(AEC)的角度看,大氣污染程度取決于污染物排放量與大氣環(huán)境容量間的平衡,在無源的情景下,大氣中污染物濃度衰減的快慢取決于大氣自凈能力的大小[2-4].通過大氣自凈能力可得到大氣污染物排放量與污染物濃度間的影響規(guī)律,繼而在大氣自凈能力不變時,由污染物濃度反演得到理想排放限值,實現(xiàn)對工業(yè)園區(qū)排放量的精準管控.
目前 PM2.5濃度預測模型主要包括機理模型和統(tǒng)計學模型[5].機理模型通過物理化學原理對大氣污染物的理化反應進行模擬分析,需要有各類氣象數據和污染源排放清單.主要包括通用多尺度空氣質量模型[6](CAMQ模式)、氣象-化學耦合模型[7](WRF-Chem模式)和嵌套空氣質量預報模式系統(tǒng)[8](NAQPMS模式)等.在缺少各類氣象數據和污染源排放清單時,機理模型難以精準預測 PM2.5濃度.統(tǒng)計學模型則利用機器學習、深度學習和統(tǒng)計學分析其歷史監(jiān)測數據發(fā)現(xiàn)變化規(guī)律并給出合理預測.模型預測精度高,在缺少探空氣象數據和污染源排放清單時也能保證較好的預測效果.目前統(tǒng)計模型主要有:支持向量回歸模型[9-11]、隨機森林[12-14]、BP神經網絡[15]以及 LSTM[16-19]、BILSTM[20-21]神經網絡等單機器學習模型.為進一步提高預測精度,多機器學習模型應運而生.Liu等[22]構建了DBN、LSTM和多層神經網絡(MLP)的組合模型.康俊鋒等[23]建立XGBoost、LSTM組合模型進行PM2.5濃度短期預測.梁澤等[24]耦合遺傳算法和RBF神經網絡構建PM2.5預測模型,具有依賴變量少、預測精度高和運算效率高等特征.AEC源于機理模型,其變化反映出污染物在大氣中的物理化學反應,但現(xiàn)有研究未實現(xiàn)機理模型和非機理模型的有機結合,未考慮 AEC變化對污染物濃度的影響,且往往籠統(tǒng)的將所有站點看作獨立個體,未考慮大氣污染物濃度的時空特征,即站點的空間相關性和污染物濃度變化的周期性[25-26].
本研究提出一種集成 AEC和時空特征的工業(yè)園區(qū)PM2.5預測模型,將工業(yè)園區(qū)氣象數據、空氣污染物數據、PM2.5歷史濃度數據及排放數據結合,在大氣自凈能力-排放源-PM2.5濃度的大氣環(huán)境容量理念下,構建 CNN-BILSTM 網絡來捕捉工業(yè)園PM2.5濃度的變化規(guī)律,對未來24h平均PM2.5濃度進行預測,為工業(yè)園區(qū)大氣重污染預警和企業(yè)排放控制提供理論支持.
濮陽市工業(yè)園區(qū)坐落于河南省濮陽市,濮陽市為石油化工型城市,污染物排放量較大,重污染天氣頻發(fā),化工、新能源新材料和裝備制造為其三大主導產業(yè).濮陽市工業(yè)園區(qū)二、三類工業(yè)企業(yè)數量較多,對城區(qū)環(huán)境質量影響較大[27].工業(yè)園區(qū)企業(yè)粉塵、廢氣無組織排放的監(jiān)管督辦力度自2018年起逐步加大,違者停業(yè)整頓.但污染排放量基數龐大,空氣質量狀況不容樂觀[28-29].根據濮陽市工業(yè)園區(qū)2018年大氣源清單,各類源 PM2.5主要前體物排放量如表1所示.
表1 工業(yè)園區(qū)2018年PM2.5主要前體物排放量統(tǒng)計(t)Table 1 Emission statistics of main precursor of PM2.5 in industrial parks in 2018 (t)
由表1可知,工業(yè)園區(qū)污染排放主要來自工業(yè)源, PM2.5主要前體物,如煙塵、二氧化碳、氮氧化物的排放量均為最高,平均占比達到 90%.工業(yè)源排放主要來自于各化工或其它制造企業(yè)的生產活動,相比于移動源、揚塵源等,可控性更強,將其作為輸入變量構建模型更利于實現(xiàn)工業(yè)園區(qū)排放管控.
使用的大氣污染物濃度數據為2018年1月1日~2020年12月31日濮陽市工業(yè)園區(qū)和周圍共計12個空氣自動監(jiān)測站點(圖1)逐小時監(jiān)測數據,包括 PM2.5、PM10、CO、SO2及 O3質量濃度.其中氣象數據來源于濮陽市地面氣象站,主要包括氣溫、相對濕度、氣壓、降水、風速風向及露點溫度等.排放數據來源于河南省重點污染源自動監(jiān)控系統(tǒng),監(jiān)測排口主要為園區(qū)內大型石油化工、火力發(fā)電企業(yè)等重點污染源的窯爐排放口、煙囪排放口和脫硫出口等,主要排放物包括二氧化硫、氮氧化物及煙塵等.
圖1 空氣自動監(jiān)測站點位置分布Fig.1 Position of air automatic monitoring points
1.2.1 大氣自凈能力指數及其計算方法 大氣自凈指數用來表征大氣自身運動中對大氣中污染物的擴散、稀釋和清除能力.根據平均濃度預報方程和箱模型,采用平均體積法可推導出單位時間、單位面積上大氣平流擴散和降水所能清除的最大污染物總量大氣自凈能力指數ASI[3],即
式中:Q為大氣污染物在一定空氣體積內的排放量,kg;S為底面積,km2;VE為大氣通風量,m2/s;Wt為雨洗常數,無量綱;R為單位時間內的降水量,mm/h; CS為污染物達標濃度,m g/m3.
大氣自凈能力指數越大,代表大氣對污染物的清除能力越強,反之,代表其越弱.
1.2.2 CNN-BILSTM 預測模型 氣象、大氣污染物濃度和污染源強之間的變化呈高度復雜性、非線性,采用傳統(tǒng)模擬方法難以快速準確的達到模擬效果.因此采用融合歷史氣象、排放和空氣質量數據的深度學習模型來擬合三者變化規(guī)律.大氣污染物濃度的影響因子變化復雜且數量繁多,包括氣象因子、排放因子和及其他大氣污染物,故利用卷積神經網絡(CNN)進行關鍵特征提取,減少冗余信息;污染物濃度變化具有強自相關性和周期性,污染物濃度變化不僅受過去的信息特征影響,而且與未來的特征變化亦息息相關.而雙向長短期記憶網絡(BILSTM)兼顧雙向時間序列中所有值的關鍵歷史長期和短期依賴關系.因此采用CNN-BILSTM組合模型擬合工業(yè)園區(qū)內PM2.5濃度變化規(guī)律.
如圖2所示,數據建模包括時間序列建模和歸一化處理,時間序列建模將原始數據劃分為不同時間步樣例,與PM2.5濃度相關的氣象、排放和其他大氣污染物濃度數據作為特征集,數據歸一化采用Min-Max標準化方法.CNN包括卷積層、池化層和全連接層;BILSTM包括前向LSTM層和后向LSTM層.模型輸入為氣象、排放和大氣污染物濃度時序數據,輸出為未來1hPM2.5濃度.預測建模包括CNN卷積層將時序數據映射為n×s矩陣,其中,n代表時序數據樣例數,s代表數據維度,經過濾波器處理后進入池化層,提取局部特征,減少冗余信息,將提取后的多組特征向量輸入BILSTM層、Dropout層和全連接層中,得到未來 PM2.5濃度輸出結果.此外,為了增強 CNN-BILSTM 網絡模型的泛化性和魯棒性,選擇Adam作為優(yōu)化算法.
圖2 CNN-BILSTM框架Fig.2 The framework of CNN-BILSTM
1.2.3 組合預測模型流程 預測模型構建流程如圖3所示,包括預處理、模型構建和模型預測.
圖3 CNN-BILSTM模型預測流程Fig.3 CNN-BILSTM model prediction process
(1)預處理:預處理包括數據預處理、特征獲取和構建訓練集.數據預處理主要包括對氣象數據、大氣污染物數據和排放數據進行異常值處理、缺失值填補和歸一化處理,異常值處理方法為3σ準則,缺失值填補采用均值填充和 MICE算法,歸一化處理采用Min-Max標準化方法;特征獲取包括通過有限體積法得到大氣自凈能力指數,結合排放數據組成大氣自凈能力-排放變化特征,通過小波分析得到PM2.5濃度變化的主周期,作為時間特征,相關性分析得到其濃度變化相似的站點數據和輔助氣象及大氣污染物,作為空間特征和輔助特征;最后利用獲取的特征構建訓練集.
(2)模型構建:秋冬季重污染天氣頻發(fā),為了保證模型污染預警能力,測試集選擇在秋冬季,數據集按照訓練集、驗證集和測試集為6:2:2的比例劃分,劃分方法為留出法,算法優(yōu)化器為Adam,以實測值和預測值的相關系數作為標準對結果進行評估,最終獲取模型的最佳參數,保存最優(yōu)模型.
(3)模型預測:將測試集數據輸入模型中,對預測值和實測值進行數學分析,統(tǒng)計其對應污染等級情況并分析預測結果.
1.2.4 評價指標 組合預測模型精度通過均方根誤差(RMSE)、平均絕對誤差(MAE)、相關系數(R2)來對預測結果進行評價.
2.1.1 AEC特征 如表2所示, ASI與PM2.5濃度呈負相關,當大氣自凈能力增強時,大氣擴散清除能力增強, PM2.5濃度呈下降趨勢; 3種主要污染排放與PM2.5呈顯著正相關, 表明當園區(qū)內排放量增加會導致PM2.5濃度上升.由圖4可知,污染物排放量與 PM2.5濃度之間呈明顯的正比關系,當污染物排放量增大時,PM2.5濃度呈上升趨勢;當排放量在一定范圍內時,隨著污染物濃度增大,ASI呈減小趨勢;即當 PM2.5濃度呈上升趨勢時,排放量呈增大趨勢,而 ASI逐漸減小.此外,大氣自凈能力與排放都呈現(xiàn)明顯的季節(jié)變化,濮陽市工業(yè)園區(qū)屬于暖溫帶半濕潤季風型大陸性氣候,不同季節(jié)氣候差異明顯,園區(qū)企業(yè)不同季節(jié)生產計劃不同,污染排放量也不盡相同.
圖4 PM2.5、ASI和排放量的日均數據分布Fig.4 Daily data distribution of PM2.5, ASI and emissions
表2 PM2.5日均濃度與ASI和排放量相關性統(tǒng)計Table 2 Correlation statistics of PM2.5 daily average concentration with ASI and emissions
季節(jié)劃分:3~5月為春季,6~8月為夏季,9~11月為秋季,12月和次年的1、2月為冬季[30].如圖5所示,隨著春夏秋冬時間線的推進,PM2.5濃度逐漸上升,污染物排放量逐漸增大.將ASI、大氣污染物排放和季節(jié)作為AEC特征輸入模型.
圖5 PM2.5和大氣污染物排放的季節(jié)變化Fig.5 Seasonal variations in PM2.5 and air pollutant emissions
2.1.2 時間特征 如圖6所示,2018~2020年濮陽市工業(yè)園區(qū) PM2.5日均濃度時間序列呈現(xiàn)明顯的周期性特征.對 PM2.5時間序列進行小波分析,繪制小波方差圖(圖7),小波方差圖中波峰對應的尺度為PM2.5的主周期,由圖7可知,在70d以內,PM2.5時間序列存在 2個振蕩周期,分別為10~20d和 30~50d,第 1主周期為37d,第 2主周期為14d.結合圖8可知,PM2.5時間序列具有較強的自相關性,但是在滯后天數大于 52d時,相關系數低于 0.2,在滯后天數為14d時相關性最高, 因此選擇日前14d的PM2.5濃度作為時間特征輸入模型.
圖6 2018~2020年工業(yè)園區(qū)PM2.5日均濃度時間序列Fig.6 Time series of PM2.5 daily average concentration in industrial parks from 2018 to 2020
圖7 PM2.5時間序列小波方差Fig.7 Wavelet variance diagram of PM2.5 time series
圖8 PM2.5時間序列時間滯后自相關系數Fig.8 Lagging self-correlation coefficient of PM2.5 time series
2.1.3 空間特征 空間特征分析主要是工業(yè)園區(qū)內監(jiān)測站點與周圍監(jiān)測站點PM2.5濃度序列的相關性分析.由于工業(yè)園區(qū)緊鄰的華龍區(qū)內有其小工業(yè)園,區(qū)域污染可能有其他貢獻源,故華龍區(qū)3個站點數據均用于分析;工業(yè)園區(qū)北面、南面及東面的清豐縣、范縣及濮陽縣均有 3個空氣質量監(jiān)測站點,由于監(jiān)測站點位置集中及與工業(yè)園區(qū)站點較遠,因此取各縣區(qū)域內 3個站點的平均值用于分析.如表3所示,華龍區(qū)的3個站點PM2.5濃度變化與工業(yè)園區(qū) PM2.5濃度變化有極強的相關性,經過實驗分析,相關性過高會產生冗余,導致精度降低[24].故只將清豐縣、范縣和濮陽縣3個縣區(qū)的站點平均值作為空間特征輸入.
表3 工業(yè)園區(qū)站點與相鄰站點PM2.5濃度相關性統(tǒng)計Table 3 Correlation statistics of PM2.5 concentration between industrial park sites and adjacent sites
2.1.4 輔助特征 PM2.5濃度變化的輔助因素包括氣象因素和其他大氣污染物濃度.如表4所示,PM2.5與降水、相對濕度和風向相關性不高,濮陽市位于中緯地帶,常年受東南季風環(huán)流的控制和影響,屬暖溫帶半濕潤季風型大陸性氣候.年平均日照時數2300~2500h,是河南省日照高值區(qū).由于氣溫變化大,降水時空分布不均,氣象災害頻繁,干旱占突出地位[31];年平均風速 2.7m/s,常年主導風向是南風、東南風,因此降水、相對濕度和風向的變化與 PM2.5相關性不高.PM2.5與PM10、NO2和CO相關性很強,與SO2的相關性較低,這與前人的研究結果相同[32-33].PM2.5與O3呈負相關,這是因為高濃度的PM2.5和 PM10會導致氣溶膠光學厚度增大而降低光化學速率,從而減少O3生成[26].
表4 PM2.5濃度與氣象因素和輔助大氣污染物相關性統(tǒng)計Table 4 Correlation statistics of PM2.5 concentration with meteorological factors and air pollutants
綜上所述,將氣溫、氣壓、風速和露點溫度4個氣象因素及PM10、NO2、CO和O34個污染物變量作為輔助特征輸入.
2.2.1 PM2.5預測結果 為了實現(xiàn)大氣污染預警且保證工業(yè)園區(qū)管控策略有效實施,選擇未來1d PM2.5的日均濃度作為預測變量.考慮到計算時長及BILSTM 的雙向搜索,設置當前模型輸入步長為1,BILSTM 中前向和后向 LSTM 隱層狀態(tài)的維度l=l'=32,64,128,測試訓練后的編碼器-解碼器結構,結果如表5所示.當 BILSTM 隱層狀態(tài)維度最大為128時,預測精度最高.隨著維度增加,結構參數和訓練難度均在增加,不過,當維度從32升高到64時,精度略有下降,因此,簡單提升網絡規(guī)模不能夠顯著增加預測準確度.
表5 不同隱層狀態(tài)維度的預測結果評價Table 5 Prediction results evaluation of different hidden layer state dimensions
CNN模型和LSTM模型在短時序列預測中表現(xiàn)優(yōu)異[34-37],為檢驗集成 AEC和時空特征的CNN-BILSTM預測模型能否提高PM2.5濃度預測精度,分別建立 CNN傳統(tǒng)模型、LSTM 傳統(tǒng)模型、CNN-BILSTM 傳統(tǒng)模型、含時空特征的 CNNBILSTM 預測模型、集成 AEC特征的 CNNBILSTM 預測模型、集成 AEC和時空特征的CNN-BILSTM預測模型進行對比驗證.如圖9所示,其中CNN模型、LSTM模型和集成AEC和時空特征的CNN-BILSTM預測模型作縱向對比,用于驗證組合模型優(yōu)越性;CNN-BILSTM傳統(tǒng)模型、僅含時空特征或AEC特征的CNN-BILSTM預測模型及集成AEC和時空特征的CNN-BILSTM預測模型作橫向對比,驗證AEC和時空特征對模型的精度增益.
圖9 預測和實測結果對比Fig.9 Comparison of predicted and true values
由圖9可見,CNN和LSTM模型在濃度預測精度上遜色于 CNN-BILSTM 模型; 在加入時空特征后,預測值的波動減小,加入 AEC和時空特征后,濃度預測精度在高值區(qū)和低值區(qū)表現(xiàn)均優(yōu)異,集成AEC特征的CNN-BILSTM模型表現(xiàn)較好,集成AEC特征而未集成時空特征的CNN-BILSTM模型發(fā)生預測值大于實測值的情況較多,而集成時空特征未集成AEC特征的 CNN-BILSTM 模型捕捉濃度高值的能力較差,集成AEC和時空特征的CNN- BILSTM模型既能在實測濃度范圍內預測,又可以捕捉污染突變,綜合表現(xiàn)最好.同時,實驗表明增加了工業(yè)源排放數據的深度學習模型預測 PM2.5濃度的精度大幅提升,顯然,只有深度學習的數據集完備且準確的情況下,預測精度才能明顯提升,這也側面印證了研究的園區(qū)內工業(yè)源是造成PM2.5污染峰值的主要來源.
如圖10所示,集成 AEC和時空特征的 CNNBILSTM 模型擬合效果最好,各項評價指標最優(yōu).從擬合線和輔助線的偏離情況來看,集成AEC和時空特征的CNN-BILSTM模型在高值區(qū)的擬合效果最好,其 RMSE值為12.912,相比傳統(tǒng)的 LSTM 模型,精度提高10%.證明集成AEC和時空特征的CNNBILSTM模型有更高的預測精度.
圖10 不同模型預測性能對比Fig.10 Comparison of prediction performance of different models
2.2.2 PM2.5污染等級預測結果 根據國家環(huán)境保護標準空氣質量指數技術規(guī)定[38],分別對PM2.5濃度的真實值和預測值進行污染程度劃分,統(tǒng)計污染天數和中度污染(4級)以上的預測準確率,為了量化估計模型的不確定性,利用置信區(qū)間來獲得模型的污染等級分類的真實誤差,置信半徑的計算公式如下所示:
式中,Rint是置信空間半徑;vacc是分類準確率;n是樣本大小;z是高斯分布臨界值,當顯著性水平為95%時取值z=1.96.
如表6所示,CNN-BILSTM模型對PM2.5的濃度高值捕捉能力明顯提升,在污染天數和重度污染及嚴重污染的預測準確率均為最高,分別為93%和83%.根據其置信區(qū)間可知其真實誤差為±3.4%和±0.3%.其中集成 AEC和時空特征后的 CNNBILSTM預測模型比傳統(tǒng)模型在預測污染天數上的準確率平均提高 13.5%;集成 AEC和時空特征的CNN-BILSTM 模型比傳統(tǒng)模型在重度污染和嚴重污染以上的預測準確率平均提升50%.
表6 不同模型污染預測準確率對比Table 6 Comparison of pollution prediction accuracy of different models
2.2.3 不同季節(jié)預測 由圖11可知,模型在秋季和冬季的預測結果較好,其中冬季的預測精度最高,RMSE、MAE和R2各指標值分別為13.846、11.482和 0.889;而在春夏季的預測結果較差,其中夏季的預測結果最差,各指標值為12.569、10.033、0.253.結合濮陽市氣候特點即工業(yè)園區(qū)排放特點,造成夏季預測結果較差原因是由于夏季雨量充沛,大氣擴散清除能力好;而在秋冬季預測精度高是由于大氣擴散條件差,園區(qū)排放量大,PM2.5濃度與影響因子的相關性更好.
圖11 組合模型不同季節(jié)預測結果Fig.11 Combination model forecast results in different seasons
3.1 工業(yè)園區(qū)日均PM2.5濃度存在顯著的周期性變化且時間變化主周期為14d;空間上PM2.5日均濃度變化存在差異性,相鄰站點的PM2.5日均濃度變化具有高度相似性.
3.2 CNN-BILSTM組合模型具有更高的預測精度和穩(wěn)定性.與CNN和LSTM模型相比, 在污染天數上的預測準確率平均提高 13.5%,在重度污染和嚴重污染天數的預測準確率平均提升 50%.可用于工業(yè)園區(qū)重污染事件預警預報.
3.3 AEC特征和時空特征不僅考慮了氣象、工業(yè)排放和 PM2.5濃度間的理化反應,又兼顧 PM2.5濃度變化的周期性、自相關性和空間相關性,因此集成AEC和時空特征的工業(yè)園區(qū)PM2.5預測模型相對于未集成 AEC和時空特征的模型具有更高的預測精度,且由于模型建立了工業(yè)排放和PM2.5濃度的映射關系,既提高了模型對 PM2.5峰值的捕捉能力,亦有助于實現(xiàn)工業(yè)園企業(yè)污染排放精準管控.
3.4 由于季節(jié)特征差異,模型在不同季節(jié)的預測能力有所差異,秋冬季預測結果較好,而在春夏季預測能力較差.