孫鳳娟 ，田勇，孫開爭，付華軒，張文娟，李敏，呂晨

(山東省濟(jì)南生態(tài)環(huán)境監(jiān)測中心，山東 濟(jì)南 250101)

全球臭氧(O3)污染已經(jīng)存在較長的歷史，近年來隨著經(jīng)濟(jì)的快速發(fā)展，NOx、VOCs排放量持續(xù)增加[1]，導(dǎo)致二次污染物O3濃度逐年升高[2-3]，臭氧污染不僅會造成巨額經(jīng)濟(jì)損失，對人體健康也具有一定危害[4-6]。2013—2020年，全國及重點(diǎn)區(qū)域O3污染問題呈加劇態(tài)勢，且具有“時間長、范圍大”的污染擴(kuò)散特點(diǎn)，在PM2.5濃度仍未實(shí)現(xiàn)全面達(dá)標(biāo)的背景下，我國大氣污染防治轉(zhuǎn)向了以PM2.5和O3為主的“二次污染”治理，因此亟需及時、準(zhǔn)確的O3預(yù)報模型，以應(yīng)對城市O3預(yù)報和污染防治策略的新要求[7-8]，減輕O3污染對人民群眾健康的危害。

空氣質(zhì)量預(yù)報方法分為統(tǒng)計模型和數(shù)值模式兩類[9]，目前數(shù)值模式在京津冀及周邊地區(qū)、長三角、珠三角等地得到了較好地應(yīng)用[10]，但數(shù)值模式方法預(yù)報結(jié)果的優(yōu)劣很大程度上依賴于模型輸入的排放源清單的精度[11]。由于建立高精度、高分辨率排放源清單難度大、模式運(yùn)行復(fù)雜且計算量龐大，數(shù)值預(yù)報方法應(yīng)用難度較大[12]。統(tǒng)計預(yù)報模型以污染物濃度數(shù)據(jù)及氣象觀測資料為基礎(chǔ)，通過多元線性回歸法、神經(jīng)網(wǎng)絡(luò)法、決策樹法等統(tǒng)計方法建立污染物和氣象因子之間的數(shù)值關(guān)系模型，具有運(yùn)算量小、使用方便、不需要污染源排放清單、經(jīng)濟(jì)高效等優(yōu)點(diǎn)[13-14]，已成為許多研究者普遍使用的方法[15-17]，在區(qū)域、地級市及區(qū)縣空氣質(zhì)量預(yù)報中得到廣泛應(yīng)用[18-19]。但傳統(tǒng)的數(shù)值預(yù)報和統(tǒng)計預(yù)報均是“單一”的決定論的預(yù)報，而由于觀測的不準(zhǔn)確以及資料分析、同化中引入的誤差，我們所得到的氣象場資料總是含有不確定性，即實(shí)際大氣的真正狀態(tài)永遠(yuǎn)也不可能被完全精確地描述出來，模擬所用的氣象場也僅僅是實(shí)際大氣可能出現(xiàn)情況的一個可能的解[20]，同時由于排放因子、水平活動等數(shù)據(jù)存在較大不確定性以及排放清單更新滯后的影響，排放清單存在較大不確定性[21]。集合預(yù)報是一種為了體現(xiàn)初值不確定性并將擾動的隨機(jī)性系列化的方法，相對于單一預(yù)報，集合預(yù)報可以得到“一群”預(yù)報值，從而可以提供一定事件的發(fā)生概率，便于預(yù)報員進(jìn)行決策分析，提高預(yù)報正確率。

目前，采用隨機(jī)擾動方法進(jìn)行的集合預(yù)報多用于溫度、風(fēng)、降水等[22-25]氣象模擬預(yù)報中，針對空氣質(zhì)量模擬預(yù)報的研究較少。針對這一情況，本研究收集了2013年—2020年5月—9月共40個月濟(jì)南市O3日最大8 h滑動平均值與氣象條件，建立了人工神經(jīng)網(wǎng)絡(luò)O3預(yù)報模型，并將神經(jīng)網(wǎng)絡(luò)模型與集合預(yù)報方法相結(jié)合，設(shè)計并實(shí)現(xiàn)了基于集合預(yù)報方法的神經(jīng)網(wǎng)絡(luò)預(yù)報模型，改進(jìn)了單一人工神經(jīng)網(wǎng)絡(luò)模型的固定性，降低了預(yù)報結(jié)果的隨機(jī)性，為O3預(yù)報方法研究開拓了新思路，為有效開展O3污染防治工作提供了決策支持。

1 數(shù)據(jù)與方法

1.1 資料數(shù)據(jù)

環(huán)境空氣質(zhì)量監(jiān)測資料來源于山東省濟(jì)南生態(tài)環(huán)境監(jiān)測中心發(fā)布的2013年—2020年5月—9月11個國控點(diǎn)全市平均O3濃度數(shù)據(jù)，為保持?jǐn)?shù)據(jù)一致性，均采用標(biāo)況數(shù)據(jù)；氣象資料來源于美國國家環(huán)境預(yù)報中心(national centers for environmental prediction，NCEP)提供的FNL(final operational global analysis)再分析資料，時間分辨率4次/d，水平分辨率為1°×1°，垂直分辨率為26層[26]，提取北京時間8：00一個時次的氣象數(shù)據(jù)。

1.2 研究區(qū)域

濟(jì)南市作為省會城市，位于山東省中西部，南依泰山，北跨黃河，背山面水，地處魯中南低山丘陵與魯西北沖積平原的交接帶上，地勢南高北低，呈淺碟狀，特殊的地形條件導(dǎo)致其空氣質(zhì)量狀況受氣象條件變化影響顯著。濟(jì)南是全國74個重點(diǎn)城市中空氣污染最為嚴(yán)重的城市之一，由于經(jīng)濟(jì)高速發(fā)展，能源、交通(特別是城市汽車尾氣)規(guī)模的擴(kuò)大、城市人口的膨脹、大型工業(yè)開發(fā)區(qū)的發(fā)展等人為活動使污染物排放量不斷增加，污染范圍不斷擴(kuò)大，城市環(huán)境空氣呈現(xiàn)出典型區(qū)域性大氣復(fù)合污染特征，O3對空氣質(zhì)量的影響程度逐年增加，已成為夏秋季影響大氣優(yōu)良率的重要因素。

濟(jì)南市2013—2015年O3污染略有減輕(圖1(a))，2016—2019年O3污染呈逐年加重趨勢，2016年O3輕度污染及以上天數(shù)為59 d，2019年增至112 d，2020年O3污染天數(shù)略有減少，為104 d，但O3重度污染天數(shù)近年來明顯增多，2017—2020年O3重度污染分別為2、3、5和7 d。伴隨全球變暖及城市化進(jìn)程導(dǎo)致的城市熱島效應(yīng)加劇，濟(jì)南市O3污染出現(xiàn)的月份整體出現(xiàn)提前現(xiàn)象(圖1(b))，2018年3月出現(xiàn)O3污染5 d，2019年3月出現(xiàn)O3污染2 d；O3污染過程(持續(xù)污染超過2 d)和最長持續(xù)時間整體上也呈上升趨勢，O3污染最長持續(xù)時間18 d，出現(xiàn)在2019年6月，2020年7月初和8月底出現(xiàn)連續(xù)2 d的O3重度污染。

圖1 濟(jì)南市2013—2020年O3污染狀況Fig.1 O3 pollution status in Jinan during 2013 to 2020

1.3 數(shù)據(jù)處理

對于神經(jīng)網(wǎng)絡(luò)模型，原始數(shù)據(jù)集的好壞直接決定了預(yù)測結(jié)果的準(zhǔn)確性，數(shù)據(jù)傳輸誤差、設(shè)備故障等都會影響原始數(shù)據(jù)集，必須對原始O3濃度數(shù)據(jù)及氣象數(shù)據(jù)集進(jìn)行刪除空值、異常值等數(shù)據(jù)清洗工作，才能有效地提升神經(jīng)網(wǎng)絡(luò)模型的預(yù)測精度。

1.4 評估方法

將模擬輸出的結(jié)果與觀測值采用兩分類法進(jìn)行統(tǒng)計分析，當(dāng)事件發(fā)生時為“有”，事件未發(fā)生時為“無”[27]。模擬結(jié)果檢驗方法包括準(zhǔn)確率(PC)、TS評分(TS)、空報率(F)、漏報率(PO)、預(yù)報偏差(B)及命中率(POD)，如式(1)～(6)所示。

(1)

(2)

(3)

(4)

(5)

(6)

式中:NA為擊中天數(shù)，表示預(yù)報發(fā)生O3污染，實(shí)況也發(fā)生的O3污染的天數(shù)；NB為漏報天數(shù)，表示預(yù)報不發(fā)生O3污染，但實(shí)況發(fā)生O3污染的天數(shù)；NC為空報天數(shù)，表示預(yù)報發(fā)生O3污染，但實(shí)況未發(fā)生O3污染的天數(shù)；ND為反擊中天數(shù)，表示預(yù)報不發(fā)生O3污染，實(shí)況也未發(fā)生O3污染的天數(shù)?？梢钥闯?，PC數(shù)值范圍為0～1，理想值為1，用以表示所有預(yù)報正確的比率；TS評分?jǐn)?shù)值范圍為0%～100%，理想值為100%，用以定量評價對O3污染的預(yù)報效果；F數(shù)值范圍為0%～100%，理想值為0%，用以表征O3污染空報比率；PO數(shù)值范圍為0%～100%，理想值為0%，用以表征O3污染漏報比率；B數(shù)值范圍為0～+∞，理想值為1，用以反映預(yù)報NA事件過多或過少的傾向；POD數(shù)值范圍為0～1，理想值為1，用以表征O3污染事件中被正確預(yù)報的概率[28]。

2 實(shí)驗設(shè)計

2.1 模擬因子篩選

2.1.1 氣象因子篩選

氣象條件是決定一次空氣污染能否形成的重要因素[29]，本研究提取FNL數(shù)據(jù)，并計算得到共103組氣象因子數(shù)據(jù)，分別計算其與O3日最大8 h滑動平均值(簡稱:O3-8 h濃度)相關(guān)性，同類氣象因子選取相關(guān)性最高的一組，最終選取相關(guān)性最高的整層相對濕度(RHclm)、地面溫度(TMPsfc)、 整層可降水(PWATclm)、整層云水(CWATclm)、1 000 hPa垂直速度(VVELprs-1 000)、500 hPa緯向風(fēng)(UGRDprs-500)、10 m經(jīng)向風(fēng)(VGRD10 m)、風(fēng)向(WD)等8組氣象因子。

圖2 O3-8 h濃度與不同氣象因子關(guān)系Fig. 2 Relation between the O3-8 h concentration and different meteorological factors

圖2(續(xù))

分析這8組氣象因子和O3-8 h濃度關(guān)系可以看出(圖2)，O3-8 h濃度與地面溫度相關(guān)性最高，呈明顯正相關(guān)關(guān)系，其次為整層相對濕度，呈明顯負(fù)相關(guān)關(guān)系，其相關(guān)性分別為0.504和-0.463，與風(fēng)向的相關(guān)性最差，為-0.181，且偏北風(fēng)時O3-8 h濃度水平相對較高，可能是由于夏季偏北風(fēng)時天況多較好，太陽輻射增強(qiáng)，加之濟(jì)南特殊的南高北低的地形，增強(qiáng)了O3的生成及本地累積。

2.1.2 污染物因子篩選

由于前期空氣污染物的累積程度直接影響后期的空氣污染程度[30]，同時考慮實(shí)際預(yù)報制作時間限制，選取前一日14時6種污染物質(zhì)量濃度。同時，光化學(xué)氧化劑Ox(NO2+ O3)是大氣氧化能力的主要指標(biāo)之一[31]，代表著大氣氧化能力的強(qiáng)弱。因此，統(tǒng)計分析前一日14時PM10、PM2.5、SO2、NO2、CO、O3及Ox濃度水平與當(dāng)日O3-8 h濃度相關(guān)性，發(fā)現(xiàn)O3-8 h濃度與14時O3及Ox濃度均存在較好的相關(guān)性(圖3)，相關(guān)系數(shù)分別達(dá)到0.501、0.465。

圖3 O3-8 h濃度與前一日14時 O3、Ox濃度散點(diǎn)圖Fig.3 Scatter diagram of the O3-8 h concentration and O3 and Oxconcentration at 14 o′clock on the previous day

2.2 人工神經(jīng)網(wǎng)絡(luò)模型搭建

多層感知器(multi-layer perceptron，MLP)也叫人工神經(jīng)網(wǎng)絡(luò)(artificial neural network，ANN)，是一種運(yùn)算模型，旨在模擬神經(jīng)系統(tǒng)構(gòu)造與功能進(jìn)行數(shù)據(jù)處理，從信息處理角度對人腦神經(jīng)元網(wǎng)絡(luò)進(jìn)行抽象模擬，由大量的節(jié)點(diǎn)(或稱神經(jīng)元)相互聯(lián)接構(gòu)成，通過不斷調(diào)整模擬神經(jīng)元之間鏈條的權(quán)值，使得整個網(wǎng)絡(luò)可以較好地擬合訓(xùn)練數(shù)據(jù)的關(guān)系。

本文構(gòu)建的人工神經(jīng)網(wǎng)絡(luò)包含2個隱含層，第一個隱含層均包含1 000個神經(jīng)元，第二個隱藏層包含500個神經(jīng)元，最優(yōu)化方案采用擬牛頓法，懲罰系數(shù)0.000 1，最大迭代次數(shù)1 000次。同時，根據(jù)上文分析結(jié)果，分別設(shè)計三套模擬方案，使用2013—2019年共1 071組樣本進(jìn)行模擬訓(xùn)練，對比不同方案的訓(xùn)練效果，選取相關(guān)性最高的一套方案(表1)。方案一選取RHclm、TMPsfc、PWATclm、CWATclm、VVELprs-1 000、UGRDprs-500、VGRD10 m、WD共8類氣象因子，訓(xùn)練結(jié)果與監(jiān)測結(jié)果相關(guān)性為0.761，為顯著相關(guān)；方案二在方案一基礎(chǔ)上加入前一日14時的光化學(xué)氧化劑Ox濃度數(shù)據(jù)，相關(guān)性明顯提高，為0.801；方案三在方案二基礎(chǔ)上加入前一日14時O3小時濃度數(shù)據(jù)，相關(guān)性為0.826，最終選取方案三作為神經(jīng)網(wǎng)絡(luò)訓(xùn)練最終方案，模型訓(xùn)練結(jié)果較好(圖4)。

表1 三種方案與O3-8 h監(jiān)測值的相關(guān)性Table 1 Correlations between monitoring and simulated values of three schemes

圖4 2013年—2019年5月—9月O3-8 h濃度監(jiān)測值(實(shí)線)及模擬值(虛線)Fig.4 Monitoring values (solid line) and neural network simulation values (dotted line) of O3-8 h during May to September of each year,starting from 2013 and ending with 2019

2.3 人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報模型搭建

與單一的確定性預(yù)報相比，通過提供定量的概率預(yù)報，集合預(yù)報可以給出多種可能性及其發(fā)生的概率，能為預(yù)報員提供包括不確定性在內(nèi)的更多預(yù)報信息[32]。通過集合成員預(yù)報結(jié)果的分布可以計算出O3-8 h濃度各等級出現(xiàn)的概率，包含了該集合預(yù)報系統(tǒng)所能提供的所有信息。

本文集合預(yù)報方法采用隨機(jī)擾動法，即經(jīng)典的蒙特卡羅法，對2020年5月—9月153組因子進(jìn)行隨機(jī)擾動，產(chǎn)生15組相互獨(dú)立的隨機(jī)擾動氣象場，連同原始?xì)庀髨?，?6組氣象場，人工神經(jīng)網(wǎng)絡(luò)模型方案采用上文選取的方案三，最終得到16組O3-8 h濃度模擬結(jié)果。根據(jù)《HJ633-2012環(huán)境空氣質(zhì)量指數(shù)(AQI)技術(shù)規(guī)定》[33]，將O3-8 h濃度模擬結(jié)果換算為對應(yīng)的空氣質(zhì)量分指數(shù)(IAQIP)及空氣質(zhì)量指數(shù)類別，得到優(yōu)、良、輕度污染、中度污染、重度污染、嚴(yán)重污染6個類別分別對應(yīng)的概率，選取預(yù)測比例最高的兩個類別作為預(yù)測類別，與實(shí)測O3-8 h分指數(shù)類別進(jìn)行對比，并做準(zhǔn)確率(PC)、TS評分(TS)、空報率(F)、漏報率(PO)、預(yù)報偏差(B)及命中率(POD)檢驗，對人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報(ENF)的預(yù)測效果和精度進(jìn)行評價。

3 實(shí)例驗證與結(jié)果分析

3.1 實(shí)例驗證

3.1.1 人工神經(jīng)網(wǎng)絡(luò)模型驗證

將2020年5月—9月共153組測試集數(shù)據(jù)輸入模型進(jìn)行計算，預(yù)測O3-8 h濃度結(jié)果與實(shí)際監(jiān)測值之間相關(guān)系數(shù)達(dá)到0.665，為顯著相關(guān)。監(jiān)測值與預(yù)測值對比情況如圖5所示，O3-8 h濃度預(yù)測值與監(jiān)測值變化趨勢呈現(xiàn)較好的一致性，但污染較重時段尤其是重污染時預(yù)測值較監(jiān)測值偏低。2020年共出現(xiàn)O3重度污染7 d，分別為6月4日、6月6日、6月10日、7月7日—8日、8月29日—30日，該模型對6月4日—10日、7月7日—8日及8月29日—30日污染過程預(yù)測效果較好，能夠明顯地反映出O3污染累積、持續(xù)及快速消散的不同階段，但預(yù)測濃度均偏低。

圖5 2020年5月—9月O3-8 h濃度監(jiān)測值及模擬值Fig.5 Monitoring values and neural network simulation values of O3-8 h during May to September 2020

3.1.2 人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報模型驗證

圖6給出了16組O3-8 h濃度模擬值箱線圖?？梢钥吹剑蓡T間的離散度很高，集合中位數(shù)變化趨勢與監(jiān)測值趨勢基本一致，且監(jiān)測值數(shù)據(jù)56%落入10百分位值至90百分位值區(qū)間內(nèi)。

根據(jù)16組O3-8 h濃度模擬結(jié)果，對比人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報及單一人工神經(jīng)網(wǎng)絡(luò)預(yù)報的PC、TS、F、PO、B及POD檢驗結(jié)果可以看出(表2)，ENF整體準(zhǔn)確率較NF明顯提高，O3污染漏報率明顯減少，O3污染命中率明顯提高，O3污染空報率較NF略有增加，B為1.26，大于1，表明ENF對O3污染預(yù)報有過多傾向，而NF則有過少傾向。

表2 人工神經(jīng)網(wǎng)絡(luò)集合及單一人工神經(jīng)網(wǎng)絡(luò)預(yù)報結(jié)果檢驗Table 2 Test results of the ensemble neural network and neural network prediction

3.2 典型案例分析

2020年7月3日—9日濟(jì)南市經(jīng)歷了一次O3典型污染過程，3日，降水過程，O3濃度最低，之后天況轉(zhuǎn)好，氣溫升高，O3生成加快，O3污染逐步加重，7日—8日，氣溫達(dá)到最高，O3達(dá)到重度污染水平，9日，弱降水過程，O3緩解，本次O3污染過程結(jié)束。

以2020年7月3日—9日O3重污染過程為例(圖7)，可以看出ENF及NF對本次污染過程的預(yù)報趨勢均較好，但污染累積階段，相較于NF，ENF能夠更好地反映出污染的迅速累積上升，對于O3污染最重的7日、8日，ENF 預(yù)測的O3-8 h濃度中位數(shù)更高。7月3日—7日污染累積階段，監(jiān)測值均落入10百分位值至90百分位值區(qū)間內(nèi)，8日污染持續(xù)階段，ENF 預(yù)測結(jié)果略低，9日污染消散階段ENF預(yù)測結(jié)果下降趨勢落后于監(jiān)測值。主要是由于參與模型訓(xùn)練的2013—2019年共1 071組數(shù)據(jù)中，O3重度污染天數(shù)僅有11 d，訓(xùn)練樣本數(shù)量的不足導(dǎo)致對O3重度污染過程的預(yù)報偏差較大，后續(xù)伴隨參與訓(xùn)練的樣本數(shù)增多，這種偏差將得以彌補(bǔ)。

注：圖中數(shù)值從上到下分別標(biāo)注10百分位值、中位數(shù)、90百分位值，星號表示離散點(diǎn)。圖7 2020 年7月3日—9日各集合成員O3-8 h預(yù)測箱線圖分布特征Fig.7 Box plot of O3-8 h simulation values for each member forecast during July 3 to 9, 2020

3 結(jié)論

(1)人工神經(jīng)網(wǎng)絡(luò)方法預(yù)測O3-8 h濃度結(jié)果與實(shí)際監(jiān)測值之間相關(guān)系數(shù)達(dá)到0.665，為顯著相關(guān)，兩者變化趨勢呈現(xiàn)較好的一致性，但污染較重時段尤其是重污染時預(yù)測值較監(jiān)測值偏低。

(2)與單一人工神經(jīng)網(wǎng)絡(luò)預(yù)報相比，集合預(yù)報準(zhǔn)確率明顯提高，O3污染漏報率明顯減少，O3污染命中率明顯提高，O3污染空報率較單一預(yù)報略有增加，人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報對O3污染預(yù)報有過多傾向，而單一人工神經(jīng)網(wǎng)絡(luò)預(yù)報則有過少傾向。通過提供定量的概率預(yù)報，人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報可以給出多種可能性及其發(fā)生的概率，能為預(yù)報員提供包括不確定性在內(nèi)的更多預(yù)報信息，對于提高O3預(yù)報準(zhǔn)確率具有較高的參考價值，該模型具有一定的實(shí)際應(yīng)用及推廣價值。

(3)以2020年7月3日—9日一次O3重污染過程為例，與單一人工神經(jīng)網(wǎng)絡(luò)的確定性預(yù)報相比，人工神經(jīng)網(wǎng)絡(luò)集合預(yù)報能夠更好地反映出污染的迅速累積上升及持續(xù)過程，但由于訓(xùn)練樣本數(shù)量的不足，導(dǎo)致對O3重度污染過程的預(yù)報偏差仍然較大，后續(xù)伴隨參與訓(xùn)練的樣本數(shù)增多，這種偏差將得以彌補(bǔ)。