京津冀及其周邊地區(qū)“2+26”城市PM2.5的蒙特卡羅集合預報試驗

曹凱 唐曉 孔磊 王威 吳倩 黃樹元 ,4 張佩文 ,4韓麗娜 ,4 吳其重 王自發(fā)

1 中國科學院大氣物理研究所大氣邊界層物理與大氣化學國家重點實驗室，北京 100029

2 中國科學院大學，北京 100049

3 中國環(huán)境監(jiān)測總站，北京 100012

4 成都信息工程大學，成都 610225

5 北京師范大學全球變化與地球系統(tǒng)科學學院，北京 100875

1 引言

京津冀及其周邊地區(qū)“2+26”城市是當前中國大氣污染最嚴重的地區(qū)之一（王恰和鄭世林,2019），近年來空氣質量雖有所改善，但大氣污染整體形勢仍十分嚴峻（程鈺等,2019），給社會經濟造成巨大損失的同時也對居民健康產生巨大影響（謝楊等,2016）。對空氣質量變化進行高精度預報和預測是防控大氣污染的關鍵環(huán)節(jié)，以多模式集合預報系統(tǒng)為核心的空氣質量預報方法已成為我國空氣質量業(yè)務預報的重要工具之一（王自發(fā)等,2009; 王茜等,2010; 陳煥盛等,2013）。

多模式集合預報是基于不同團隊研發(fā)的模型構建（王自發(fā)等,2009; Marécal et al.,2015），通過采用統(tǒng)一的輸入數據得到預報集合。在多模式集合預報中，采用合適的集成方法將集合預報系統(tǒng)中不同模式成員集成起來能顯著提升空氣質量預報效果（Pagowski et al.,2005; 王自發(fā)等,2008）。Monteiro et al.（2013）利用多模式集合和中位數法、靜態(tài)線性回歸法、動態(tài)線性回歸法、貝葉斯模型平均法等集成方法，顯著提升了PM10和臭氧的預報效果。黃思等（2015）及潘錦秀等（2019）將多模式集合預報與線性回歸集成方法結合起來提升了城市空氣質量預報效果。吳劍斌等（2017）評估了多模式空氣質量數值預報業(yè)務系統(tǒng)對城市臭氧的預報效果，發(fā)現(xiàn)采用最優(yōu)化集成方法的預報效果明顯優(yōu)于單個模式。張?zhí)旌降龋?019）采用了均值集成、權重集成、多元線性回歸集成、BP神經網絡集成和最優(yōu)集成等方法提升了多模式PM2.5濃度的預報準確率。謝磊等（2019）基于多個空氣質量預測模式，利用最優(yōu)定權組合法顯著提升了SO2預測精度。

多模式集合預報的優(yōu)點在于考慮各模式物理參數化方案的不確定性，但集合樣本數量相對較少，且通常無法充分考慮模式輸入場的不確定性，而模式輸入場不確定性往往是模式不確定的重要來源（唐曉等,2010）。因此 Hanna et al.（1998）利用蒙特卡羅方法（Mullen and Baumhefner,1994; Du et al.,1997）分析了輸入參數不確定性對UAM-IV模型預測結果不確定性的影響。此外，唐曉等（2010）針對模式輸入數據不確定性概率分布特征，利用蒙特卡羅方法對模式輸入場進行集合擾動，構建了包含50個集合樣本的蒙特卡羅集合預報系統(tǒng)，并開展臭氧概率預報的嘗試。

相較于單個空氣污染數值模式，蒙特卡羅集合預報系統(tǒng)的集合樣本數量較多，且能考慮模式輸入場的不確定性，為結果提供更為豐富和全面的預報信息。因此本文首先利用蒙特卡羅模擬方法搭建了空氣質量的多擾動集合預報系統(tǒng)，并評估此系統(tǒng)對“2+26”城市PM2.5質量濃度的預報效果，然后根據統(tǒng)計參數對集合樣本進行篩選，最后利用“集合樣本優(yōu)選”均值集成法對“2+26”城市PM2.5質量濃度進行確定性預報并評估其預報效果。

2 集合預報系統(tǒng)及數據來源簡介

2.1 集合預報系統(tǒng)介紹與模擬設置

2.1.1 集合預報系統(tǒng)介紹

嵌套網格空氣質量預報模式系統(tǒng)NAQPMS（王自發(fā)等,2006）是中國科學院大氣物理研究所自主研發(fā)的三維歐拉化學傳輸模式，能模擬大氣中污染物的平流、擴散、氣相化學、液相化學、干濕沉降等物理化學過程，在大氣污染的科研和業(yè)務預報上得到了廣泛應用。NAQPMS氣相化學機制采用 CBM-Z（Carbon Bond Mechanism Z）碳鍵反應機制，考慮了71種化學物質以及133個核心化學反應；濕沉降和液相化學過程基于RADM中的機制，包含了22種氣體和氣溶膠；干沉降采用Wesely 方案（Wesely,1989）；氣溶膠熱力學模塊基于 ISORROPIA（Li et al.,2012）。本研究基于NAQPMS搭建了PM2.5的蒙特卡羅集合預報系統(tǒng)。蒙特卡羅集合預報系統(tǒng)框架如圖1所示。由于受排放因子、水平活動數據的不確定性以及排放清單更新緩慢的影響，排放清單存在較大不確定性（Cao et al.,2011），因此本次試驗僅將排放源中的污染物排放按照其不確定性特征進行隨機擾動，產生一組代表排放源不同可能狀態(tài)的隨機集合樣本，然后將集合擾動的排放源樣本輸入至NAQPMS中進行模式積分，最后得到不同輸入樣本狀態(tài)下的PM2.5質量濃度集合預報樣本。同時考慮到計算資源有限，本次試驗僅隨機抽取50組擾動樣本，從而得到50組擾動狀態(tài)下PM2.5濃度預報樣本。對排放源的擾動方法如下所示：

其中，X0表示擾動前不同物種的排放初值，Xi表示排放物種i擾動后的排放大小，Pi表示排放物種i的擾動系數，其遵循均值為1、標準差為 σ的對數正態(tài)分布， σ則通過表1中不同排放物種不確定度大小確定，具體數值來自 Zhang et al.（2009）。

表1 不同排放物種不確定度Table 1 Uncertainties of emissions for different species

2.1.2 氣象場模擬設置

本次試驗利用中尺度氣象模式WRF模擬得出的氣象場作為NAQPMS的動力驅動。在設置WRF參數化方案時，分別選用RRTM和Dudhia方案作為長波輻射方案和短波輻射方案，陸面過程選擇Noah方案，邊界層方案采用YSU方案，微物理方案采用WSM3方案。運行WRF所需的初始和邊界條件來自美國國家環(huán)境預報中心的FNL再分析數據集，空間分辨率為1°（緯度）×1°（經度），時間分辨率為6 h。另外在每天的氣象場模擬中，WRF自由運行36 h，其中前12 h作為模式啟動時間，后24 h為NAQPMS提供氣象場數據。

2.1.3 NAQPMS 模擬設置

NAQPMS模擬區(qū)域及“2+26”城市地理位置如圖2所示。采用一層嵌套，以（34°N，105°E）為中心覆蓋東亞大部分地區(qū)，網格分辨率為15 km，網格個數為432（經向） ×339（緯向）。垂直方向上采用α?z地形追蹤坐標系，不等距分為20層，其中2 km以下設置9層，模式層頂海拔高度為20 km。由全球大氣化學模式MOZART提供初始場和邊界條件。排放源數據中，人為源選用HTAP_v2清單（Janssens-Maenhout et al.,2015）、生物質燃燒源來自 GFED_v4 清單（Randerson et al.,2017）、生 物 源 來 自 MEGAN-MACC（ Sindelarova et al.,2014）、海洋 VOCs來自 POET（http://www.aero.jussieu.fr/projet/ACCENT/POET.php[2019-06-11]）、土壤及閃電 NOx數據分別采用 Yan et al.（2005）和 Price et al.（1997）。試驗時段為 2017 年 9～12月，在對污染物進行逐日循環(huán)模擬時，集合預報系統(tǒng)都會在初時時刻生成一組排放源的擾動集合，隨后將排放源的擾動集合輸入至NAQPMS中進行積分，從而得到PM2.5的預報集合。

圖1 蒙特卡羅集合預報系統(tǒng)框架Fig.1 Framework of the Monte Carlo ensemble forecast system

圖2 模擬區(qū)域范圍及“2+26”城市地理位置Fig.2 Simulated domain and locations of "2+26" cities

2.2 數據來源簡介

PM2.5質量濃度觀測數據來自中國環(huán)境監(jiān)測總站，其中包括“2+26”城市共計161個監(jiān)測站點2017年9～12月逐時觀測濃度，并以該市PM2.5濃度觀測值的站點平均代表該市PM2.5質量濃度觀測值；PM2.5質量濃度預報值數據來自蒙特卡羅集合預報系統(tǒng)50個集合樣本的預報結果。為方便與觀測值比較，根據161個監(jiān)測站點的地理位置，從集合系統(tǒng)預報結果中提取對應網格的PM2.5質量濃度預報值，同樣以該市PM2.5質量濃度預報值的站點平均代表該市PM2.5質量濃度預報值。

3 統(tǒng)計集成方法及評估指標簡介

3.1 統(tǒng)計集成方法簡介

3.1.1 均值集成法

均值集成法是一種簡單且常用的集成方法，其預報值可表示為各集合樣本預報值的算數平均，公式如下：

其中，F(xiàn)mean、m、Fk分別表示均值集成法的預報值、集合樣本個數以及第k個集合樣本的預報值。

3.1.2 “集合樣本優(yōu)選”均值集成法

Boylan and Russell（2006）曾指出平均分數偏差（Mean Fractional Bias，MFB）和平均分數誤差（Mean Fractional Error，MFE）可作為模式對顆粒物質量濃度模擬準確性的衡量指標，具體計算公式如下：

其中，MFB和MFE分別表示MFB和MFE，Mi、Oi、n分別表示第i個時刻模擬值、第i個時刻觀測值、有效樣本對數，并認為最優(yōu)模式能夠達到的準確性標準為?30%≤MFB≤30%、MFE≤50%。當存在多個集合樣本時，利用MFB和MFE可剔除不確定性較大的樣本以提高集合樣本的整體可靠性。在此基礎上，對優(yōu)選出的集合樣本采用均值集成法作為集合預報系統(tǒng)的確定性預報結果。為方便描述，對集合樣本優(yōu)選后采用均值集成的統(tǒng)計方法稱為“集合樣本優(yōu)選”均值集成法。

3.2 評估指標簡介

均方根誤差（RMSE）和相關系數（r）分別表示預報值與觀測值之間平均偏離程度和變化趨勢相似程度的統(tǒng)計量，具體計算公式如下：

其中，RMSE表示RMSE，Pi、Oi、n、、分別表示第i時刻的預報值、第i時刻的觀測值、有效樣本對數、預報均值以及觀測均值。

4 試驗結果

4.1 集合均值預報

為評估蒙特卡羅集合預報系統(tǒng)的預報能力和特點，結合觀測資料首先分析各集合樣本的PM2.5預報時間序列特征，然后利用RMSE和r等統(tǒng)計參數定量評估均值集成法的預報效果。

圖3是保定、滄州、邢臺、邯鄲、鄭州、太原、陽泉、德州等城市利用蒙特卡羅集合預報系統(tǒng)得出的各集合樣本PM2.5質量濃度預報值與觀測值時間序列對比。由于蒙特卡羅集合預報系統(tǒng)是對排放源進行擾動，并將擾動后的結果入至同一模式中，因此各集合樣本間PM2.5質量濃度變化趨勢比較一致。此外，各集合樣本所構成的預報值集合雖在絕大部分時間段能包含觀測值，但大部分樣本存在較大的先驗模擬偏差，這使得均值集成法的預報值偏高，其余城市也有類似現(xiàn)象。

圖3 2017年9～12月“2+26”城市中部分城市蒙特卡羅集合預報系統(tǒng)PM2.5質量濃度預報值與觀測值時間序列對比：（a）保定；（b）滄州；（c）邢臺；（d）邯鄲；（e）鄭州；（f）太原；（g）陽泉；（h）德州Fig.3 Comparison of time series of PM2.5 concentration observations and forecast values of the Monte Carlo ensemble forecast system from September to December 2017 in some cities of "2+26" cities: (a) Baoding; (b) Cangzhou; (c) Xingtai; (d) Handan; (e) Zhengzhou; (f) Taiyuan; (g)Yangquan; (h) Dezhou

圖4 2017年9～12月“2+26”城市均值集成法PM2.5濃度均方根誤差（RMSE）、預報值與觀測值之間相關系數（r）Fig.4 The RMSE (Root Mean Squared Error) of PM2.5 concentration and correlation between observations and forecast values of PM2.5 concentration with the ensemble mean method in "2+26" cities from September to December 2017

RMSE與r能較好表征出預報值與觀測值之間平均偏離程度和變化趨勢相似程度，因此利用其定量分析均值集成法對PM2.5質量濃度的預報效果，結果如圖4所示。對所有集合樣本采用均值集成后，各城市PM2.5濃度預報值與觀測值之間相關系數r均大于0.5，r城市均值接近0.7。從平均偏離程度來看，各城市 RMSE 主要集中在 45～65 μgm-3，均值集成法的 RMSE 城市均值接近 55 μgm-3。可以看出，均值集成法雖在PM2.5質量濃度變化趨勢上有著較好的把控能力，但由于大部分集合樣本存在較明顯的先驗模擬偏差，使得均值集成法的預報結果高估。

4.2 “集合樣本優(yōu)選”均值預報

考慮到排放清單不確定性對模式結果影響較大（Kong et al.,2019），使得均值集成法預報存在較顯著的先驗模擬偏差，我們嘗試了另外一種方法來集成集合樣本。該方法采用滾動預報的方式，首先結合MFB、MFE等統(tǒng)計參數從前30天集合預報樣本中優(yōu)選出符合最優(yōu)模式準確性標準的樣本，然后針對優(yōu)選出的集合樣本，計算其第31天的集合平均值，作為未來1天PM2.5的確定性預報結果。

為了評估兩種方法的差異，首先對比分析均值集成法和“集合樣本優(yōu)選”均值集成法的總體預報效果，預報時段均為2017年10～12月，再分析其對各污染等級的預報能力。圖5給出了北京、天津、唐山、鄭州、鶴壁、新鄉(xiāng)、濟南、濟寧等城市采用兩種統(tǒng)計集成方法后的PM2.5質量濃度預報值與觀測值的時間序列對比?？梢钥闯觯c所有集合樣本均值預報相比，“集合樣本優(yōu)選”均值集成法能有效降低預報偏差，使得預報值更加貼近實際觀測值。

為進一步分析兩種方法的預報技巧，利用RMSE、r和觀測—模擬兩倍因子百分比（FAC2）等統(tǒng)計參數來評估兩種集成統(tǒng)計方法的預報效果。由圖6a可見，采用“集合樣本優(yōu)選”均值集成法后，各城市PM2.5預報的RMSE均顯著減小，RMSE 城市均值由 58.0 μg m?3降低至 34.7 μg m?3，r城市均值由0.69提升至0.70。FAC2表示預報值落于0.5～2倍觀測值范圍內的比例（Chang and Hanna,2004），若FAC2越大則表明預報精度越高。從圖6b、6c看出，均值集成法的FAC2為67%，有相當一部分預報值大于2倍觀測值，高估較為嚴重。而采用“集合樣本優(yōu)選”均值集成法后，高估現(xiàn)象明顯改善，將PM2.5預報的FAC2指標提升至87%。

相對作用特征（ROC）是基于雙態(tài)分類聯(lián)列表對雙態(tài)事件進行檢驗。對于一次事件的發(fā)生與否，預報可分為預報正確、漏報、空報和正確否定四種情況。命中率（hit rate）=預報正確數/（預報正確數+漏報數），假警報率（false alarm rate）=空報數/（空報數+正確否定數）。以假警報率為橫坐標，命中率為縱坐標即可構成ROC散點圖。若散點落在隨機猜測線（即對角線）上或右側，則表明不具備任何預報技巧；若落在隨機猜測線左側且離隨機猜測線垂直距離越遠，則表明對該事件預報技巧越高。

圖5 2017年10～12月“2+26”城市中部分城市不同集成統(tǒng)計方法預報值與觀測值時間序列對比：（a）北京；（b）天津；（c）唐山；（d）鄭州；（e）鶴壁；（f）新鄉(xiāng)；（g）濟南；（h）濟寧Fig.5 Comparison of time series of forecast values and observations by different ensemble statistical methods from October to December 2017 in some cities of "2+26" cities: (a) Beijing; (b) Tianjin; (c) Tangshan; (d) Zhengzhou; (e) Hebi; (f) Xinxiang; (g) Jinan; (h) Jining

圖6 2017年10～12月“2+26”城市不同集成方法的（a）PM2.5濃度均方根誤差（RMSE）、PM2.5濃度預報值與觀測值之間相關系數（r）；2017年10～12月“2+26”城市（b）均值集成法、（c）“集體樣板優(yōu)選”均值集成法PM2.5濃度預報值與觀測值對比Fig.6 (a) RMSE of PM2.5 concentration of different ensemble methods from October to December 2017 in "2+26" cities.Comparison of observations and forecast values of PM2.5 concentration with (b) the ensemble mean method and (c) "collective sample selection" ensemble mean method from October to December 2017 in "2+26" cities

圖7 不同集成方法各污染等級的ROC單點分布（DROC表示散點到對角線的垂直距離，負值表示散點位于對角線右側）Fig.7 Single point distribution of ROC (Receiver Operating Characteristic) for each pollution level with different ensemble methods(DROC is the vertical distance from scatter to the diagonal,and a negative value indicates that the scatter is located to the right of the diagonal)

結合《環(huán)境空氣質量指數（AQI）技術規(guī)定（試行）》（中華人民共和國生態(tài)環(huán)境部，2016），利用ROC散點圖對比分析兩種集成統(tǒng)計方法對空氣污染等級分別“優(yōu)”、“良”、“輕度污染”、“中度污染”、“重度污染”、“嚴重污染”等事件的預報技巧，結果如圖7所示。采用所有集合樣本的均值集成法時，空氣污染等級為“良”、“輕度污染”和“中度污染”的DROC值分別為0.01、－0.02和0.06，說明均值集成法對上述事件的預報技巧非常低。而采用“集合樣本優(yōu)選”均值集成法后，“優(yōu)”、“良”、“輕度污染”、“中度污染”、“重度污染”等事件的DROC值明顯增大，預報技巧顯著提升，說明“集合樣本優(yōu)選”均值集成法對各污染等級的整體預報技巧要高于均值集成法。值得注意的是，“集合樣本優(yōu)選”均值集成法對“嚴重污染”事件的預報技巧略低于均值集成法，這種現(xiàn)象可能與所有集合樣本的均值集成法存在較大模擬偏差有關，這種顯著高估導致其對“嚴重污染”事件預報正確的概率較高。

5 結論

本文在嵌套網格空氣質量預報模式系統(tǒng)NAQPMS的基礎上，結合蒙特卡羅模擬方法搭建多擾動空氣質量集合預報系統(tǒng)。從該系統(tǒng)對京津冀及其周邊地區(qū)“2+26”城市PM2.5質量濃度預報結果可以看出，在所有集合樣本均值集成法的預報值存在較大先驗模擬偏差的情況下，采用“集合樣本優(yōu)選”均值集成法能顯著改善各城市PM2.5預報效果，減小PM2.5預報偏差。

此外，多擾動空氣質量集合預報系統(tǒng)具有一定局限性，PM2.5預報效果存在較大提升空間。首先，本次試驗主要考慮排放源的不確定性，對模式自身物理、化學過程不確定性以及氣象場不確定性考慮較少，根據 Li et al.（2019）研究表明，氣象場不確定性對京津冀地區(qū)重霾期間的PM2.5模擬效果產生較大影響，因此未來可在多模式集合預報系統(tǒng)基礎上，結合蒙特卡羅模擬方法對排放源和氣象場進行擾動，構建多模式多擾動的超級集合預報系統(tǒng)，從而更好考慮輸入場和模型動力框架不確定性對PM2.5預報效果的影響；其次，對集合樣本采用合適方法進行集成之前可利用偏差訂正等技術降低集合系統(tǒng)所存在的先驗模擬偏差；最后，本研究僅針對符合最優(yōu)模式準確性標準的集合樣本采用均值集成法來提升PM2.5質量濃度預報效果，即賦予每個優(yōu)選集合成員相同的權重系數，未來可考慮進一步利用多元線性回歸、神經網絡模型等方法對優(yōu)選出的集合樣本進行集成以此來提升PM2.5預報精度。