空氣質量數值預報優化方法研究

趙俊日 ,肖 昕*,吳 濤,李彥鵬,賈紅霞 (.中國礦業大學環境與測繪學院,江蘇 徐州 6；.徐州環境監測中心站,江蘇 徐州 8；.江西應用技術職業學院,江西 漳州 000；.中國環境新聞工作者協會 北京 00095)

作為當前大氣環境科學研究的熱點與難題,大氣污染物預報預警可通過各類預報方法與手段相結合,對多種大氣污染物在全球尺度下的不同類型污染過程進行模擬預測研究,成為城市.及區域大氣復合污染控制研究的重要手段之一[1-3].大氣污染物排放源清單是模型研究和相關控制策略制定的重要基礎,但排放源估算取決于許多因素,包括社會經濟、能源、土地利用、環境資料等[4].但由于統計數據的滯后性和排放因子及時空分配系數等數據的不確定性,定量估算排放源及對其進行實時更新有較大的困難,有必要通過各類“反演”模型采用自下而上的方法反演污染源提高空氣質量模型的模擬精度[5].徐祥德等[6]首次提出了“nudging”源同化方法,通過迭代計算的方式逐步修正污染源,改進NO2和 SO2濃度預報效果;靳璐濱等[7]利用三維變分同化(3D-var)方法對青奧會期間南京地區的 PM10、PM2.5進行了同時同化實驗,結果表明同化對初始場有顯著改進.

CMAQ模型作為美國環保署推薦的第3代空氣質量模擬系統Models-3的核心模塊[8],對污染物區域分布及其變化趨勢雖具有較強的預報能力,但由于污染源時空特征十分復雜,CMAQ污染濃度預報量與實況相比存在明顯偏低的“系統性”誤差[9],且CMAQ模型的物理、化學機制及參數等需不斷完善,因此很多學者對空氣質量模型采用后處理的方式對其進行訂正來減少模型誤差[10].許建明等[11]提出 CMAQ-MOS動力-統計相結合的空氣質量預報技術,其試驗方案表明可顯著降低CMAQ模型的“系統性”預報誤差;尤佳紅等[12]采用WRF-RTIM的MOS方法進行冬季污染預報和霾的診斷,發現在歷史污染監測資料小樣本條件下,RTIM 的統計預報具有較高的精度.

然而目前國內對中小尺度空氣質量數值預報優化方面的研究較少,本文綜合考慮排放源清單的不確定性及模式在其背景下產生的非客觀性預報偏差,采用源清單修正方法、統計修正方法對其進行優化,分析修正前后預報效果的改善程度,以期為區域空氣質量改善和重污染天氣預報預警提供科學依據.

1 模式設置與數據

1.1 模式設置

本文采用的空氣質量模式為CMAQ(v4.7.1),CMAQ為開源的高度模塊化結構模型,已被廣泛應用于常規空氣質量預報和決策制定當中[13],CMAQ模式所需要的氣象場可由中尺度數值天氣預報模式 WRF(v3.6.1)提供,該模式對多種物理過程均有不同的參數化方案,可以更為真實的模擬大氣空間運動[14-16].源排放處理模式SMOKE (v3.0)可以將排放清單處理成空氣質量模式需要的時空分辨率數據[17],本文中 SMOKE模式主要用于徐州本地排放源清單處理,對其進行時間、空間及化學物種分配,外兩層和徐州版圖區域之外的剩余部分采用清華大學 MEIC[18]清單排放數據.模式采用的主要參數化方案如表1、2所示.

如圖 1所示,本文中所有模式采用三層嵌套方式模擬,WRF模式和CMAQ模式采用蘭伯特投影坐標系,坐標原點為117°E和45°N,兩條真緯線為北緯30°和北緯60°.三層水平網格分辨率和網格數分別為 3km,9km,27km 和 97×97,94×112,76×88.第一層區域覆蓋中國華北,華東及中部地區,部分東亞地區;第二層區域覆蓋蘇北周邊的河南,山東,安徽等省份,第三層區域為目標研究區域覆蓋整個徐州市.

圖1 研究區域Fig.1 Map of Research Regional

表1 WRF模式參數設置Table 1 Parameters setting of WRF model

表2 CMAQ模式參數設置Table 2 Parameters setting of CMAQ model

1.2 使用數據

本文使用的WRF模式輸入數據源于美國國家環境預測中心發布的FNL全球再分析數據,水平分辨率為 1°×1°,地形和下墊面數據分別是USGS 30s全球地形和MODIS土地利用分類產品.徐州市本地排放源清單(基準年2015年)和大氣污染物觀測數據及溫度、濕度、風速、風向、氣壓等氣象觀測數據由徐州環境監測中心提供.

2 研究方案設計

2.1 源清單“Nudging”修正方法

本文采用“Nudging”源同化反演方法,即在空氣質量模式積分方程中構造排放源的“張弛逼近項”[19]以減少模式預報結果與實際觀測濃度的誤差.

假設某時刻的模擬結果誤差完全是由排放源S產生,固定其他大氣物理化學過程,已知實際觀測污染物濃度為C*,經過n次迭代計算,模式預報的污染物濃度逐漸向實際觀測值逼近,即 Cn→C*可以得到如下迭代式子:

但本研究中 CMAQ模式的第三層區域水平網格分辨率為 3km×3km,若不考慮周圍面源影響,只考慮監測站點周圍的點源源強,修正之后的源清單對 CMAQ模式預報結果影響甚微,所以引入反距離加權(IDW) 插值算法[20]來彌補上述缺點.

對迭代式子中設:

假設在CMAQ模式中第三層區域里有k個監測站點,各落在第(xi,yi)網格上,設其網格上的α值為αi,i=1,2,3…,k,IDW 的插值函數可描述為如下:

式中:是網格上的任意點到(xi,yi)的水平距離;p是大于0的加權冪指數

(本文中p取1.5),考慮1,2,3…,m種不同污染物,最后得到以下通式:

2.2 統計修正方法

XGBoost是在 2015年由美國華盛頓大學Chen[21]提出,XGBoost是 Gradient Boosting Machine(梯度提升機器學習算法)的C++的實現,能自動利用 CPU的多線程進行并行計算,且對算法加以改進運行速度和精度得到了很大的提升[22].

假設模型有 k個決策樹,其集成模型可以如下表示:

對其目標函數二階泰勒展開得:

本文基于XGBoost算法對CMAQ模式預報的6種常規大氣污染物進行統計修正.設t時間段為起報時間段,t1為5×24h,取t-t1時間段的大氣污染物觀測值和氣象要素觀測值及CMAQ模式模擬值作為訓練數據,t2為 1×24h,取 t+t2時間段的WRF模式氣象要素預報值及CMAQ模式預報值作為預測數據,建立統計修正回歸模型,可以描述為如下方程:

式中:為t- t1時間段的污染物濃度觀測值矩陣;為t- t1時間段的n×t1維的氣象要素觀測值矩陣和m×t1維的CMAQ模型污染物濃度預報值矩陣;為 t+t2時間段的 n×t2維的WRF氣象要素預報值矩陣和 m×t2維的 CMAQ模型污染物濃度預報值矩陣;本文中n和m分別取溫度、濕度、氣壓、風速、風向等5列和PM2.5、PM10、O3、SO2、NO2、CO 等 6列,為t- t1時間段的(n+m)×t1維矩陣,為常數項矩陣.

3 結果與討論

3.1 源清單“Nudging”修正結果與分析

源同化試驗階段:采用2016年12月2～6日徐州市 13個監測站點的 SO2、NO2、CO實測數據進行迭代計算.

控制試驗階段:采用初始排放源模擬 2016年12月3～31日期間徐州市區SO2、NO2、CO濃度.

預報效果檢驗試驗:采用同化源模擬 2016年12月3～31日期間徐州市區SO2、NO2、CO濃度.

表 3給出了 2016年 12月 SO2、NO2、CO逐時濃度預報、修正值與觀測值的相關系數(R)、平均絕對誤差(MAE)、平均相對偏差(MFB)、平均相對誤差(MFE)、均方根誤差(RMSE). Boylan等[23]建議以平均相對偏差(MFB)和平均相對誤差(MFE)為衡量指標評估模式預報的合理性,假如MFB值在-60%～60%之間且MEF小于75%,則可認為模式模擬結果在合理的可接受的范圍內,若 MFB值在-30%～30%之間且 MEF小于50%,則可認為模式表現優秀,模擬結果在理想水平范圍內.

從表3可以看出,冬季徐州市SO2、NO2逐時預報值與觀測值的相關系數比較低,CO預報值與觀測值的相關系數相對來說比較高,在 0.5左右,3類污染物的平均相對誤差和平均相對偏差都在理想水平范圍內,但平均絕對誤差和均方根誤差均較大.

采用初始源模擬的結果相比,修正之后 3類污染物的相關系數均有所提高,SO2、NO2的修正值與觀測值的平均絕對誤差和均方根誤差明顯減少,平均相對誤差分別為 33.64%和 28.45%,平均相對偏差分別為 3.85%和-11.39%,可認為修正結果均在理想水平范圍內.從修正前后的各項指標的變化幅度來看,3類污染物中 NO2修正效果最好,其次是 SO2,CO最差,其平均絕對誤差和均方根誤差及MFB和MFE取值范圍反而變大,且相關系數提高的幅度也相對較小,優化效果不太理想,可能其模擬效果已經達到最優狀態,預報效果很難進一步得到改善,下一步研究有必要對源排放、湍流擴散、化學反應、干濕沉降等過程進行分析,找出對局地 CO污染貢獻最大的過程進行修正優化.

表3 2016年12月CO、SO2、NO2源清單修正誤差統計表Table 3 Error statistics of simulated CO、SO2、NO2 concentration using initial and inversion emission during December 2016

圖2所示為徐州市冬季2016年12月SO2、NO2、CO初始源和同化源源強空間分布,可以看出 3類污染物排放源的空間分布呈現大市區排放源強較大,其他各市區較小,冬季采暖期間燃煤及工業生產排放的污染物主要集中在市區,且冬季氣象條件不利于污染物的擴散.同化源的源強大小空間分布與初始源基本相符,CO同化源的源強比初始源小,而SO2、NO2同化源的源強相對初始源高出一定比例,說明初始源清單對 CO的排放量存在一定的高估,而對SO2、NO2排放量存在一定的低估.

雖然修正之后預報效果提高了,至于清單修正之后的源強是否接近于真實排放源,與模式內部物理、化學反應參數本地化和排放源清單時間、空間、化學輪廓是否正確及使用的觀測資料站點數等多種因素有關.總體而言,根據上述兩種排放源的模擬結果與源強的空間分布差異一定程度上反映了冬季重污染天氣過程中CO、NO2、SO2排放源強的動態分布特征.

圖2 CO、NO2、SO2的初始源和同化源SMOKE處理結果空間分布(moles/s)Fig.2 CO、NO2、SO2 Spatial distributions of emission source intensity of SMOKE results第一、二行分別為初始源和同化源

3.2 統計修正結果與分析

對于初始排放源模擬的 PM2.5、PM10、O3及經“nudging”方法修正的同化源模擬的SO2、NO2、CO等6種污染物,采用基于XGBoost算法的結合 WRF預報氣象要素的統計修正方法對CMAQ模型進行后處理.為了檢驗其修正效果,對2016年12月3～31日冬季徐州市13個監測站點逐時污染物濃度觀測值、預報值及相應的修正值進行分析對比,本文取2m溫度、2m相對濕度、10m風場、地面氣壓等5種WRF氣象要素作為修正建模的訓練樣本,訓練樣本數選取參考程興宏等[24]的方法,并考慮地區適用差異,本文采用前5d的數據作為訓練樣本進行滾動修正.

從表4可以看出,徐州區域SO2、NO2逐時預報值與觀測值相關系數比較低,PM2.5、PM10的預報值平均絕對誤差和均方根誤差都比較大,而CO、O3的相關系數相對較高.經上述方法進行滾動修正之后 6種污染物預報修正值與觀測值的相關系數得到了顯著的提高,各項修正誤差也減少了許多.對于 CO、O3修正之后其相關系數均提高了 0.15左右,經“nudging”修正之后出現的CO的MFB和MFE兩項指標反而變大的現象也成功消去,平均絕對誤差,均方根誤差也減小了 50%左右.就相關系數而言,NO2改進幅度最大,均方根誤差和平均絕對誤差明顯減少,平均相對誤差和平均相對偏差分別為 17.24%和-2.91%,可認為修正效果非常好.PM2.5和 PM10修正值與觀測值相關系數分別提高到 0.6～0.7左右,平均絕對誤差和均方根誤差分別減少了23.7和 23.2μg/m3,其減少幅度最大,平均相對誤差和平均相對偏差均在理想水平范圍內.對于SO2,修正值與觀測值的相關系數提高了0.25左右,平均絕對誤差和均方根誤差分別減少了12.4和6.8μg/m3,平均相對誤差和平均相對偏差分別為1.60%和24.37%,其修正前后各項統計指標變化幅度較明顯,但相比其他 5個污染物,其預報優化效果不太理想,可能與 CMAQ模式本身對SO2模擬效果不太好有關[25],需要對模式內部物理、化學過程進行進一步的研究,找出對SO2模擬影響最大的因子進行優化.

表4 2016年12月統計修正誤差統計表Table 4 Error statistics of statistical revisions during December 2016

圖3 2016年12月預報優化結果散點Fig.3 The result of the six pollutants optimization scatter plot in December 2016 sim為模擬值,opt為修正

散點圖可以反映模擬值與觀測值的線性關系,可一目了然的看出修正效果好壞程度.圖 3給出了修正前后的預報值與觀測值的散點分布特征,可以看出CO、O3預報值與觀測值的散點多數在對角線上方,大部分預報值有一定的高估,散點分布較集中,相關系數也較高;而SO2、NO2預報值與觀測值的散點多數在對角線下方,擬合線與對角線距離較遠,即大多數預報值明顯低于觀測值,相關系數也偏低;PM2.5、PM10預報值與觀測值的相關系數不低,但是其散點分布最為分散.經上述方法修正后,預報偏低或偏高現象得到很大的改善,除了SO2之外,相關系數均提高到0.6～0.7左右,且模擬值與觀測值的分布較集中在對角線附近,說明修正之后預報值更接近于實際觀測值.

4 結論

4.1 在源清單“Nudging”修正部分,本文結合IDW空間插值算法,以SO2、NO2、CO濃度預報值與觀測值的誤差迭代算子逐步修正污染源分布的方式,對模式預報值進行修正.3類污染物初始排放源的空間分布呈現大市區排放源強較大,其他各市區較小,初始源清單對CO的排放量存在一定的高估,而對SO2、NO2排放量存在一定的低估現象.從修正前后的各項指標的變化幅度來看,3類污染物中 NO2修正效果最好,其次是SO2,CO是最差.

4.2 為了改進模式在具有較大不確定性的排放源清單影響背景下,各種物理、化學過程產生的非客觀性預報偏差,本文采用統計修正的方式對CMAQ模式進行后處理.用集成學習算法建立統計修正模型時,考慮到訓練和預報階段的氣象條件,引入了與污染物濃度有顯著相關的觀測氣象要素和WRF預報氣象要素作為不同時間段的建模自變量.經滾動修正之后,成功消除了經“nudging”修正之后出現的反常現象,除了SO2之外,相關系數均提高到 0.6～0.7左右,優化效果頗為明顯.

[1]唐孝炎,張遠航,邵 敏.大氣環境化學-第 2版 [M]. 北京:高等教育出版社, 2006:447-449.

[2]王占山,李曉倩,王宗爽,等.空氣質量模型CMAQ的國內外研究現狀 [J]. 環境科學與技術, 2013,(S1):386-391.

[3]薛文博,王金南,楊金田,等.國內外空氣質量模型研究進展 [J].環境與可持續發展, 2013,38(3):14-20.

[4]Yumimoto K, Uno I. Adjoint inverse modeling of CO emissions over Eastern Asia using four-dimensional variational data assimilation [J]. Atmospheric Environment, 2006,40(35):6836-6845.

[5]孟 凱,程興宏,徐祥德,等.基于CMAQ源同化反演方法的京津冀局地污染源動態變化特征模擬研究 [J]. 環境科學學報,2017,37(1):52-60.

[6]Xu X, Xie L, Cheng X, et al. Application of an adaptive nudging scheme in air quality forecasting in China [J]. Journal of Applied Meterology and Climatology, 2008,47(8):2105-2114.

[7]靳璐濱,臧增亮,潘曉濱,等.PM2.5和 PM2.5~10資料同化及在南京青奧會期間的應用試驗 [J]. 中國環境科學, 2016,36(2):331-341.

[8]Dennis R L, Byun D W, Novak J H, et al. The next generation of integrated air quality modeling: EPA's models-3 [J]. Atmospheric Environment, 1996,30(12):1925-1938.

[9]謝 敏,鐘流舉,陳煥盛,等.CMAQ 模式及其修正預報在珠三角區域的應用檢驗 [J]. 環境科學與技術, 2012,35(2):96-101.

[10]Djalalova I, Monache L D, Wilczak J. PM2.5analog forecast and Kalman filter post-processing for the Community Multiscale Air Quality (CMAQ) model [J]. Atmospheric Environment, 2015,119:431-442.

[11]許建明,徐祥德,劉 煜,等.CMAQ-MOS區域空氣質量統計修正模型預報途徑研究 [J]. 中國科學, 2005,35(z1):131-144.

[12]尤佳紅,束 炯,陳亦君,等.基于MOS的杭州秋冬季空氣污染預報和霾診斷 [J]. 中國環境科學, 2014,34(7):1660-1666.

[13]薛文博,許艷玲,唐曉龍,等.中國氨排放對PM2.5污染的影響 [J].中國環境科學, 2016,36(12):3531-3539.

[14]胡向軍,陶健紅,鄭 飛,等.WRF模式物理過程參數化方案簡介[J]. 甘肅科技, 2008,24(20):73-75.

[15]Hu X M, Nielsengammon J W, Zhang F. Evaluation of Three Planetary Boundary Layer Schemes in the WRF Model [J].Journal of Applied Meteorology & Climatology, 2010,49(9):1831-1844.

[16]Grell G A, Peckham S E, Schmitz R, et al. Fully coupled “online”chemistry in the WRF model [J]. Atmospheric Environment, 2005,39(37):6957-6975.

[17]王 剛.杭州市二次形成PM2.5的研究 [D]. 浙江大學, 2008.

[18]He K. Multi-resolution Emission Inventory for China (MEIC):model framework and 1990-2010 anthropogenic emissions [C]//AGU Fall Meeting, 2012.

[19]孟智勇,徐祥德,陳聯壽.衛星亮溫資料四維同化方案及其對“7·20”武漢特大暴雨的模擬試驗 [J]. 大氣科學, 2002,26(5):663-676.

[20]余小東,武 瑩,何臘梅.反距離加權網格化插值算法的改進及比較 [J]. 工程地球物理學報, 2013,10(6):900-904.

[21]Chen T, Guestrin C. XGBoost: A Scalable Tree Boosting System[C]. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2016:785-794.

[22]張 鈺,陳 珺,王曉峰,等.Xgboost在滾動軸承故障診斷中的應用 [J]. 噪聲與振動控制, 2017,37(4):166-170.

[23]Boylan J W, Russell A G. PM and light extinction model performance metrics, goals, and criteria for three-dimensional air quality models [J]. Atmospheric Environment, 2006,40(26):4946-4959.

[24]程興宏,刁志剛,胡江凱,等.基于CMAQ模式和自適應偏最小二乘回歸法的中國地區 PM2.5濃度動力-統計預報方法研究 [J].環境科學學報, 2016,36(8):2771-2782.

[25]陸維青,江峰琴,劉麗霞,等.江蘇省空氣質量預報與實測結果比對研究 [J]. 環境監控與預警, 2017,9(1):10-14.