基于機器學習方法的西安市數值模擬優化研究

李 娟， 尉 鵬， 戴學之， 趙 森， 張博雅， 呂玲玲， 胡京南

1.西安科技大學測繪科學與技術學院， 陜西 西安 710054 2.中國環境科學研究院大氣環境研究所, 北京 100021 3.合肥市氣象局， 安徽 合肥 230041

近年來，我國大氣污染問題突出，引起社會高度關注，大氣污染物對人體健康、大氣能見度以及氣候變化等都有重要影響[1-2]，因此對大氣污染物進行預報、調控及污染機理分析成為當下亟待解決的科學問題. 大氣化學模式是研究大氣污染的重要工具，它以大氣動力學為基礎，考慮多種物理和化學過程，定量描述污染物的擴散和輸送規律，但由于物理化學機理復雜、排放源不確定性等原因，使得數值預報結果存在不確定性[3-4].

為提高模式預報能力，結合數值模式與統計方法對模式預報結果進行統計修正的應用也較為普遍. 研究表明，利用機器學習方法可以有效提高模式預報準確率[5-6]，如利用多元線性回歸(MR)、隨機森林、BP神經網絡、深度學習以及集合預報等方法對模式輸出結果進行優化[7-9]. 王茜等[10]采用學習型線性回歸方法對上海市ρ(PM2.5)數值預報結果進行修正，統計檢驗結果明顯提高了ρ(PM2.5)預報效果；張岳軍等[11]在BREMPS預報結果的基礎上，結合MR、BP及多層遞階建立10 d滾動修正模型，PM2.5預報準確率升高；潘錦秀等[12]利用多元線性回歸方法對北京市多模式系統(CMAQ、CAMx、NAQPMS)預報值和觀測值進行集合預報，改進了模式的高估或低估現象. 但線性回歸模型無法解決非線性問題，隨著機器學習的發展，許多學者利用機器學習方法對污染物預報值進行訂正優化，如張恒德等[13]基于CUACE、BREMPS、WRF-Chem模式預報產品，利用BP神經網絡方法建立多模式集成預報模型，集成預報結果的歸一化平均偏差及均方根誤差均下降；戴李杰等[14]應用支持向量機(SVM)和粒子群優化算法建立PM2.5滾動預報模型，使得預報準確度提高；馬井會等[15]結合WRF氣象預報數據、地面及高空氣象觀測數據、ρ(PM2.5)觀測數據，基于人工智能深度學習序列到序列的算法建立了上海市PM2.5統計修正預報模型，能有效提升ρ(PM2.5)預報精度. 機器學習可以很好地解決非線性問題，但目前大多數研究集中在單個機器學習模型對預報值進行優化，而對多種機器學習模型優化結果的對比以及模型對不同污染物優化效果評估的研究相對較少.

為提高CAMx污染物模擬值精確度，該研究以西安市為研究對象，基于CAMx模式污染物濃度預報數據、WRF模式氣象要素預報數據以及環境空氣質量監測國控點污染物濃度觀測數據，利用6種機器學習優化模型(多元線性回歸、嶺回歸、lasso回歸、決策樹、隨機森林以及支持向量機)，對CAMx污染物濃度預報數據進行優化，旨在修正PM2.5、O3模擬值預報誤差，以期為西安市大氣空氣質量預報預警方法做出改進.

1 數據與方法

1.1 數據來源

該研究所用氣象數據來自中尺度氣象模式WRF(V4.1)，其中WRF模式每日對初始場進行初始化，每次模擬時長為31 h，Spin-up時間設置為6 h. WRF模式中心點坐標為35 °N、110 °E，采用三層區域嵌套網格，水平分辨率分別為27、9、3 km，依次覆蓋整個中國地區、陜西省和關中地區，與CAMx模式對應. WRF三層嵌套網格數分別為238×168、94×121、169×121，覆蓋垂直層為35層. 氣象背景場和邊界條件資料來自NCEP(National Centers for Environmental Prediction)的再分析逐日資料FNL，分辨率為1°×1°，時間分辨率為6 h. 地形和下墊面輸入資料分別來自USGS30全球地形和MODIS下墊面分類資料. 積云參數化方案采用Kain-Fritsch(new Eta)方案，邊界層參數化方案采用Mellor-Yamada-Janjic(Eta)湍流動能方案，大氣輻射方案為RRTM長波和云(Dudhia)短波輻射方案.

ρ(PM2.5)及ρ(O3)小時模擬值由綜合空氣質量模式CAMx提供，模擬時間為2019年1月1日—12月31日. CAMx模式采用Lambert投影，氣相及液相化學機理分別為CB05、RADM-AQ；氣溶膠熱力學平衡模式為ISORROPIA，干沉降參數化方案采用WESELY89；水平平流及垂直擴散方案分別采用PPM、隱式歐拉方案. 模式第1次運行時模擬預測5 d的初始場，以消除初始條件的影響，后續初始場采用前一次的模擬結果，以消除排放源及初始條件累積誤差.

人為源采用2016年MEIC (中國多尺度排放清單模型，Multi Resolution Emission Inventory for China)排放清單，分辨率為0.25°×0.25°；天然源運用陸地生態系統估算模型MEGAN[16]計算，整合、處理后的天然源、人為源排放清單共同輸入到排放源處理模型SMOKE，轉化為網格化的模式源排放文件.

西安市ρ(PM2.5)及ρ(O3)小時觀測數據來自環境專業知識服務系統(http://envi.ckcest.cn/environment)，空氣質量監測站點選擇西安市長安氣象監測站(站點號:57039).

1.2 模型建立

1.2.1特征提取及特征縮放

研究表明，氣象因素可以影響大氣污染物的稀釋擴散、積聚清除，是影響污染物濃度的重要因素之一[1]，其中，較低的風速不利于污染物擴散，主導風向影響著污染物的區域輸送，較高的相對濕度有利于顆粒物的吸濕增長，而高溫則加劇了臭氧的光化學反應，同時天氣形勢的演變也是污染現象的一個重要原因[17-20]. 可見，風速、風向、相對濕度、溫度、氣壓等氣象因子對大氣污染物濃度分布有較大影響. 結合西安市氣象特點及相關研究[21-25]，該研究從CAMx模式模擬結果中分別選取PM2.5、O3小時質量濃度模擬值，從WRF中提取對應時間的溫度、相對濕度、海平面氣壓、風速及風向小時值5個氣象因子，共計6個因子作為建模訓練的輸入特征，其中回歸目標值分別為相應時刻的ρ(PM2.5)、ρ(O3)小時觀測值，并在數據集中隨機選擇80%作為訓練集進行優化模型訓練，20%作為測試集進行優化結果驗證.

特征縮放可以消除各維數據之間的量級差別，提升分類算法和優化算法的性能，訓練之前需要對數據進行標準化預處理[26-27]. 特征縮放的主要方法有歸一化和標準化，其中，標準化方法使得特征值呈正態分布，利于訓練階段的權重更新；同時，標準化方法還可以保持異常值的信息，進一步減少異常值對算法的影響[27]. 標準化函數為

z=(x-u)/s

(1)

式中，x為訓練樣本，u為訓練樣本平均值，s為訓練樣本標準偏差.

1.2.2模型選取

在大氣化學模式預報結果的基礎上，利用機器學習模型對模擬結果進行優化，可以有效提高預報結果精確度[11]. 目前，機器學習模型中常用模型有多元線性回歸[12]、嶺回歸[28]、lasso回歸、決策樹、隨機森林[29]以及支持向量機[14].

多元線性回歸(普通最小二乘法，ordinary least squares，OLS)通過尋找回歸參數、回歸常數及回歸殘差，使訓練集預測值與真實回歸目標值y之間的均方誤差最小，利用訓練集計算出特征值與觀測值之間的函數關系，并將得出的方程應用于測試集，其回歸方程表示為

Y=β0+β1X1+β2X2+…+βkXk+ε

(2)

式中：Y為模型預測結果;β0為回歸常數;β1～βk均為回歸系數;X1～Xk均為特征值，該研究中有6個特征值;ε為回歸殘差.

嶺回歸的預測公式與多元線性回歸相同，通對預測模型進行顯式約束，避免系數過擬合. Lasso回歸通過正則化產生稀疏權值矩陣，使某些回歸常數βk剛好為0，對特征進行自動化選擇. 決策樹通過反復切分數據進行回歸或分類，與支持向量機、神經網絡相比計算量大大降低，而且算法不受數據縮放影響. 隨機森林通過隨機的方式建立許多決策樹并組成森林，是決策樹的集成模型，隨機森林相較于其他機器學習模型，可以量化的方式體現出模型參數對模型預報效果的影響，且預報效果較好、預報結果穩定. 支持向量機(Support Vector Machine，SVM)基于VC維理論和結構分析風險最小原則的理論，在解決小樣本、非線性及高維模式識別中表現出許多特有的優勢[30-32].

1.3 評估方法

對CAMx模式模擬結果與模型優化結果進行時間序列分析及統計檢驗. 統計檢驗指標選取均值偏差(MB)、標準化均值偏差(NMB)、平均相對偏差(MFB)、標準平均誤差(NME)、平均相對誤差(MFE)、均方根誤差(RMSE)、平均絕對誤差(MAE)及相關性系數(R)，計算公式見式(3)～(8). MB與RMSE反映了模擬值和觀測值之間的偏差和誤差大小，其絕對值越小，表明數值模擬結果與觀測結果越接近. NMB、NME反映模擬值與觀測值之間相對偏差大小，一般情況下，如果二者均小于50%，則認為模型模擬效果較好[33]. 當MFB≤±60%、MFE≤75%時，符合文獻[34]中空氣污染物模擬標準，則認為模式結果可靠；當MFB值在-30%～30%之間且MEF小于50%時，則可認為模式表現優秀，模擬結果在理想水平范圍內.

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

2 結果與討論

2.1 數值模擬結果

2.1.1WRF模擬結果

氣象小時觀測數據來自西安市長安氣象監測站(站點號:57039)，利用觀測數據對WRF模擬數據進行驗證，檢驗參數包括均方根誤差(RMSE)及相關性系數(R). 其中，溫度及氣壓模擬值均與觀測值的相關性較好，R分別為0.93、0.97，誤差較小，RMSE值分別為2.96 ℃、2.88 hPa；而風速、相對濕度模擬結果均與其觀測值相關性稍低，分別為0.70、0.67，RMSE值分別為1.38 m/s、23.1%. WRF模擬值誤差分析可知，WRF整體模擬誤差較小、精度較高，模擬結果可靠.

2.1.2CAMx模擬結果

為了解CAMx模式的預報特點，首先評估測試集中ρ(PM2.5)、ρ(O3)模擬值與觀測值的相關性及時間序列分布. 由圖1可見：ρ(PM2.5)模擬值與觀測值的散點分布較為散亂，R為0.63；ρ(O3)模擬值與觀測值分布較ρ(PM2.5)集中，R為0.78，模擬效果較好. 由圖2可見：ρ(O3)模擬值與觀測值時間變化序列一致，ρ(O3)模擬值與觀測值較為接近，模擬效果較為理想；ρ(PM2.5)模擬值存在高估現象，模擬偏差主要出現在峰值階段，而排放源的不確定性以及模式本身在多氧化過程和濕清除過程的不確定性對模式模擬結果有較大影響[34-35]. 分析其原因，一方面CAMx模式所用氣象場來自WRF，WRF模擬氣象場具有不確定性，而氣象與污染物的相互反饋作用對污染物的分布有重要影響，從而造成CAMx模擬值偏差；另一方面，人為源排放量及排放類型變化速度快，而排放清單編制時效滯后也會帶來預報的不確定性，同時MEIC清單本身存在誤差，且排放清單分辨率低于模式網格分辨率，從而導致造成模式值的污染物空間分布具有差異性[36-39].

圖1 ρ(PM2.5)、ρ(O3)模擬值與觀測值的散點分布Fig.1 Scatter distribution of simulated and observed values of ρ(PM2.5) and ρ(O3)

表1為測試集中ρ(PM2.5)及ρ(O3)模擬值與觀測值統計評估參數. 由表1可見：ρ(PM2.5)的MB值為74.07 μg/m3，表明模式值對觀測值高估；NMB、NME值均大于1，表明模式值高估程度較高；MFB值較小而MFE值偏大，說明存在部分ρ(PM2.5)模擬值被低估現象；RMSE、MAE值分別為174.00、99.85 μg/m3，說明模式模擬結果誤差及偏離程度較大.ρ(O3)整體模擬效果較好，其MB值為-3.99 μg/m3，說明模式值對ρ(O3)觀測值整體略微低估；NMB值為-8.02%，表明模式值低估現象不明顯，模式值與觀測值較為接近；MFB、MFE值分別為-19.12%、52.86%，說明模擬結果可靠；RMSE、MAE值分別為37.11、26.81 μg/m3，表明模擬結果誤差較小，模擬結果精度較高；R值為0.78，說明模擬值與觀測值之間相關性較好. 總體來看，CAMx模式模擬結果能夠較好地反映西安市ρ(O3)狀況及變化趨勢.

2.2 機器學習優化結果檢驗

2.2.1ρ(PM2.5)優化效果檢驗

該研究采用機器學習模型對PM2.5小時模擬結果進行優化，結果如圖3所示. 由圖3可見，對比ρ(PM2.5)觀測值及相對應時間的CAMx模擬值、CAMx模型優化值發現，6種機器學習模型對CAMx的ρ(PM2.5)修正效果明顯，相關性系數可提至0.70～0.78，散點由分散分布變為線性集中分布，其中，多元線性回歸、嶺回歸及lasso回歸優化后相關性系數提至0.70，決策樹及支持向量機優化后R值分別提至0.72、0.74，隨機森林優化后R值提至0.78.

圖2 ρ(PM2.5)、ρ(O3)模擬值與觀測值的時間序列分布Fig.2 Time series distribution of simulated and observed values of ρ(PM2.5) and ρ(O3)

表1 ρ(PM2.5)、ρ(O3)模式值檢驗評估參數統計

圖3 ρ(PM2.5)優化結果散點分布Fig.3 Scatter distribution of ρ(PM2.5) optimization results

由表2可見，經機器學習模型優化后，各統計參數明顯降低，ρ(PM2.5)訂正結果顯著. 多元線性回歸優化后，MB、NMB、NME值分別從74.07 μg/m3、132.6%、178.75%降至1.37 μg/m3、2.5%、43.94%，模式值高估的現象得到明顯改善，RMSE、MAE分別下降了137.12、75.31 μg/m3，離散程度大大降低；嶺回歸、lasso回歸與多元線性回歸的優化結果差別較小，模式優化結果達理想水平. 經決策樹優化后，MAE、MFB、MFE值分別降至20.91、5.54、26.21 μg/m3，優化結果合理性提高；但RMSE為39.37 μg/m3，在6種模型中最大，說明優化值與觀測值之間的偏差較大. 隨機森林優化后MB、NMB值分別為-5.70 μg/m3、-10.21%，支持向量機的優化后MB、NMB值分別為-7.87 μg/m3、-14.09%，表明隨機森林與支持向量機的優化結果均出現低估現象，其中支持向量機低估現象較為明顯. 隨機森林優化后RMSE、MAE值分別降至34.36、16.24 μg/m3，優化結果誤差最小，檢驗指標達到要求標準，且該模型得到優化結果的相關性系數為0.78，R值在6個模型中最大. 綜上，各優化模型有不同優勢，綜合分析可得，隨機森林對ρ(PM2.5)的優化結果最優.

表2 ρ(PM2.5)優化模型評估參數統計Table 2 Evaluation parameters of ρ(PM2.5) optimization results

圖4 ρ(O3)優化結果散點分布Fig.4 Scatter distribution of ρ(O3) optimization results

2.2.2ρ(O3)優化效果檢驗

由圖4可見，6種機器學習模型對ρ(O3)優化效果表現優秀，相關性系數提高了6.4%～12.8%，優化結果分布更為合理，擬合度較高，其中，決策樹、隨機森林優化后R值分別提至0.83、0.85，多元線性回歸、嶺回歸及lasso回歸優化后R值提至0.86，支持向量機R值提至0.88.

結合表3的ρ(O3)優化結果發現：6種機器學習模型對ρ(O3)的優化效果理想，ρ(O3)模擬結果誤差減小. 多元線性回歸及嶺回歸對ρ(O3)優化效果一致，NB、NMB值均較小而NME值均較大；同時，MFB、MFE值分別為-17.27%、68.25%，對ρ(O3)部分數據有低估現象. Lasso回歸MFE值降至47.54%，優化結果合理性提高. 隨機森林的MB及NMB值均最大，低估現象較明顯. 決策樹RMSE值為28.82 μg/m3，誤差最大，R值為0.83，小于其他5種優化模型，優化性能相對較差. 與其他5種機器學習模型相比，支持向量機優化后的RMSE、MAE值分別為23.76、17.40 μg/m3，預測值與觀測值之間偏離程度低、誤差小，MB、NMB、NME值較為合理，MFB、MFE值差異小，且擬合程度最好. 綜上，支持向量機對ρ(O3)優化結果更為合理.

表3 ρ(O3)優化模型評估參數統計Table 3 Evaluation parameters of ρ(O3) optimization results

2.3 機器學習優化效果評估

2.3.1優化原理分析

CAMx模型的計算過程實質上是求解每個網格中每種污染物的物理化學變化連續方程，連續方程通過分子分裂法積分，進而計算每個主要過程(包括排放、平流、輸送、擴散、去除和化學反應)對每個網格單元格內污染物濃度改變的獨立貢獻，模型核心表達式[40]：

(11)

2.3.2優化效果評估

圖5為ρ(PM2.5)及ρ(O3)觀測值、CAMx模擬值及優化模型優化結果分布對比. 由圖5(a)可見：ρ(PM2.5)的CAMx模擬數據分布離散，與觀測值差異較大，經機器學習模型優化后，CAMx模擬數據高估現象明顯改善，數據分布形態與觀測值更為接近. 對比ρ(PM2.5)觀測值，多元線性回歸、嶺回歸、lasso回歸及支持向量機的優化結果對高濃度觀測值出現低估現象，決策樹優化后對高濃度觀測值略有高估，而隨機森林優化結果與觀測值形態分布最為接近. 由圖5(b)可見：與ρ(O3)觀測值相比，CAMx模擬值出現峰值高估、整體低估現象，經機器學習模型優化后CAMx模擬數據整體更加符合觀測值分布. 多元線性回歸、嶺回歸、lasso回歸優化后，ρ(O3)觀測峰值被低估、ρ(O3)平均值提高，與觀測值數據分布存有差異；隨機森林優化后改善了觀測峰值濃度被高估或低估現象，但整體ρ(O3)值偏低；決策樹優化結果與ρ(O3)觀測值整體較為接近，但其誤差較大；支持向量機對ρ(O3)峰值略有低估，但優化結果與觀測值最為接近.

圖5 ρ(PM2.5)、ρ(O3)模擬值與優化值統計分布對比Fig.5 Statistical distribution comparison of simulated and optimized values of ρ(PM2.5) and ρ(O3)

各優化模型對ρ(PM2.5)、ρ(O3)優化結果精確度提高率如表4所示. 優化前，ρ(PM2.5)的CAMx模擬值與觀測值之間的RMSE為174.00 μg/m3(見表1)，CAMx模擬結果優化后，RMSE大幅降低，預報精度整體提高了77%～80%，其中經隨機森林優化后的ρ(PM2.5)精度提高了80%，對ρ(PM2.5)的CAMx模擬值優化效果最好. 優化前，ρ(O3)的CAMx模擬值與觀測值之間的RMSE為37.11 μg/m3(見表1)，CAMx模擬效果較ρ(PM2.5)模擬效果要好，機器學習模型優化后的預報精度提高了22%～36%，其中支持向量機優化效果最好. 對比各機器學習模型對ρ(PM2.5)與ρ(O3)的優化結果發現，決策樹對ρ(PM2.5)和ρ(O3)的優化效果均較差；隨機森林對ρ(PM2.5)的優化效果最好，但對ρ(O3)優化效果較差；多元線性回歸、嶺回歸及lasso回歸對ρ(PM2.5)的優化效果較差，但對ρ(O3)的優化效果較好.

表4 ρ(PM2.5)、ρ(O3)優化模型結果精度提高率Table 4 The accuracy improvement rate of ρ(PM2.5) and ρ(O3) optimization results

3 結論

a) 2019年西安市污染物模擬濃度評估結果表明，CAMx模式在一定程度上能夠反映ρ(PM2.5)、ρ(O3)的變化趨勢，由于排放源滯后、模式分辨率較低等原因，使預報結果存在系統性偏差，ρ(PM2.5)、ρ(O3)模擬均方根誤差分別為174.00、37.11 μg/m3；其中ρ(PM2.5)的CAMx模擬值對觀測值出現高估現象，高估現象主要出現在ρ(PM2.5)峰值階段，ρ(O3)的CAMx模擬值對觀測值則略有低估.

b) 利用機器學習模型對CAMx模擬結果進行優化后，ρ(PM2.5)模擬值與觀測值的相關性系數提至0.70～0.78，ρ(O3)模擬值與觀測值的相關性系數提至0.83～0.88；不同優化模型修正程度不同，其中隨機森林對ρ(PM2.5)優化后，其結果的穩定性和修正趨勢整體優于其他方法，優化后精度提了80%，但對ρ(O3)模擬值優化效果相對較差；而支持向量機對ρ(O3)模擬值優化結果最接近觀測值，優化后精度提高了36%.