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

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

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

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

1 模式設(shè)置與數(shù)據(jù)

1.1 模式設(shè)置

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

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

圖1 研究區(qū)域Fig.1 Map of Research Regional

表1 WRF模式參數(shù)設(shè)置Table 1 Parameters setting of WRF model

表2 CMAQ模式參數(shù)設(shè)置Table 2 Parameters setting of CMAQ model

1.2 使用數(shù)據(jù)

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

2 研究方案設(shè)計

2.1 源清單“Nudging”修正方法

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

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

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

對迭代式子中設(shè):

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

式中:是網(wǎng)格上的任意點到(xi,yi)的水平距離;p是大于0的加權(quán)冪指數(shù)

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

2.2 統(tǒng)計修正方法

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

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

對其目標函數(shù)二階泰勒展開得:

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

式中:為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維矩陣,為常數(shù)項矩陣.

3 結(jié)果與討論

3.1 源清單“Nudging”修正結(jié)果與分析

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

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

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

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

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

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

表3 2016年12月CO、SO2、NO2源清單修正誤差統(tǒng)計表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類污染物排放源的空間分布呈現(xiàn)大市區(qū)排放源強較大,其他各市區(qū)較小,冬季采暖期間燃煤及工業(yè)生產(chǎn)排放的污染物主要集中在市區(qū),且冬季氣象條件不利于污染物的擴散.同化源的源強大小空間分布與初始源基本相符,CO同化源的源強比初始源小,而SO2、NO2同化源的源強相對初始源高出一定比例,說明初始源清單對 CO的排放量存在一定的高估,而對SO2、NO2排放量存在一定的低估.

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

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

3.2 統(tǒng)計修正結(jié)果與分析

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

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

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

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

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

4 結(jié)論

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

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

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

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

[3]薛文博,王金南,楊金田,等.國內(nèi)外空氣質(zhì)量模型研究進展 [J].環(huán)境與可持續(xù)發(fā)展, 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源同化反演方法的京津冀局地污染源動態(tài)變化特征模擬研究 [J]. 環(huán)境科學學報,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]. 中國環(huán)境科學, 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 模式及其修正預報在珠三角區(qū)域的應用檢驗 [J]. 環(huán)境科學與技術(shù), 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區(qū)域空氣質(zhì)量統(tǒng)計修正模型預報途徑研究 [J]. 中國科學, 2005,35(z1):131-144.

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

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

[14]胡向軍,陶健紅,鄭 飛,等.WRF模式物理過程參數(shù)化方案簡介[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]孟智勇,徐祥德,陳聯(lián)壽.衛(wèi)星亮溫資料四維同化方案及其對“7·20”武漢特大暴雨的模擬試驗 [J]. 大氣科學, 2002,26(5):663-676.

[20]余小東,武 瑩,何臘梅.反距離加權(quán)網(wǎng)格化插值算法的改進及比較 [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模式和自適應偏最小二乘回歸法的中國地區(qū) PM2.5濃度動力-統(tǒng)計預報方法研究 [J].環(huán)境科學學報, 2016,36(8):2771-2782.

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