基于隨機森林算法的山東省區域臭氧污染事件中氣象條件和排放貢獻影響研究

龔安保，解 歡，于陽春，沈楠馳，杜雨飛，張 坤，陳姝芮，吳 彤，管 旭*，李 莉

1. 山東省生態環境規劃研究院，國家環境保護陸海統籌生態治理與系統調控重點實驗室，山東 濟南 250101

2. 上海大學環境與化學工程學院，上海 200444

近年來，隨著經濟的快速發展以及工業化和城市化的廣泛推進，能源消費持續增加，大氣污染問題愈發凸顯. 自《大氣污染防治行動計劃》[1]實施以來，以細顆粒物(PM2.5)為首要污染物的大氣污染得到一定控制，PM2.5和二氧化硫(SO2)質量濃度明顯下降，但O3污染日益嚴重[2-4]. O3對于人體健康、生態系統和氣候變化都有重要影響[5-9]. 近地面O3主要由其前體物(如NOx和VOCs) 在光照條件下通過一系列光化學反應生成[10]，此外，氣象條件對于O3污染的擴散與傳輸也起到至關重要的作用. 因此，O3觀測值往往受到排放及氣象條件的共同影響[11-14]. 如何識別排放及氣象條件對O3污染的影響，對于評估已有管控政策的實施效果、制定新的政策具有重要意義.

當前已有多種方法用來評估排放和氣象條件對O3的貢獻，如化學傳輸模型、數理統計和機器學習等. 化學傳輸模型通過調整排放清單和氣象設置，設計多情景分析以表征排放和氣象條件對空氣質量的影響[15-16]. 但化學傳輸模型作為一種數值模擬，其對計算機資源需求相對較高，此外，排放清單和氣象場數據的不確定性均在一定程度上弱化了解析結果的可靠性. 數理統計方法主要利用多元線性回歸、經驗正交函數或KZ 濾波等方法構建氣象或排放數據和空氣質量觀測濃度之間的關系[17-20]. 這些方法通常旨在通過數學方程擬合或剝離單個影響因素，但由于O3與其前體物之間存在強烈的非線性關系，且受到多個氣象因素的影響，因此往往難以用統計數據準確捕獲它們之間的關系. 近年來隨著大數據分析技術的快速崛起，機器學習方法開始被廣泛用于剝離氣象條件和排放對空氣質量的影響[21-25]. 例如，Grange 等[21]利用隨機森林模型研究了1997-2016 年瑞士PM10的排放趨勢；Vu 等[22]利用氣象標準化的方法評估了階段性的“清潔空氣行動計劃”的效果和影響；Dai 等[24]運用機器學習的方法量化了特定事件中氣象條件和排放的貢獻. 這些研究已充分證明了機器學習方法可以有效剝離氣象條件和排放對空氣質量的影響. 相較于其他方法，機器學習方法低廉的使用成本無疑是其另一優勢，這也使得機器學習方法有著更廣泛的應用潛力.

山東省地處我國東部沿海地區，東部毗鄰日本和韓國，南北承接長三角、京津冀兩個國家級城市群，是中國重要的城鎮密集地區[26]. 山東省經濟發展迅速，隨之而來的大氣O3污染問題日益嚴峻，這也嚴重制約了山東的空氣質量持續改善和經濟可持續發展[27-28].目前山東省O3污染的相關研究多集中于時空變化特征和源解析等[28-30]，較缺乏長時間周期、大空間尺度的評估排放和氣象條件對區域O3污染的貢獻. 基于此，本研究利用2020-2022 年山東省地面監測站點的污染和氣象數據，運用隨機森林的方法剝離了排放和氣象條件對于O3污染的貢獻，進一步運用沙普利加和的方法分析了不同氣象因子對O3的影響程度，并定義了區域O3污染事件，對區域O3污染事件發生時段的排放和氣象條件貢獻進行了剖析，以期為區域O3污染成因剖析和防控治理提供參考.

1 數據與方法

1.1 研究區概況

山東省位于中國東部沿海地區，瀕臨渤海和黃海，地處黃河下游，地跨34°22.90′N～38°24.01′N、114°47.50′E～122°42.30′E 之間，自北至南與河北、河南、安徽、江蘇四省接壤. 山東省地形以平原丘陵為主，中南部山地突出，西南、西北低洼平坦，東部是緩丘起伏的山東半島.

1.2 數據來源與數據處理

山東省近地面國控站點2020-2022 年O3逐小時濃度數據來自中國環境監測總站(http://www.cnemc.cn).地面氣象及污染物觀測站點分布見圖1. 對于部分異常數據，根據《環境空氣質量標準》(GB 3095-2012)[31]中針對數據有效性的規定進行剔除；對于部分缺失數據，則采用線性內插法[32]補全被剔除數據和缺失數據，即

圖1 山東省氣象監測站和空氣質量監測站分布Fig.1 The distribution of meteorological and air quality monitoring stations in Shandong Province

式中，Xn為缺失數據，XP和Xq分別代表缺失前、后的值.氣象數據來源于國家氣象數據中心(http://data.cma.cn)，觀測因子包括溫度(T)、相對濕度(RH)、風速(WS)、風向(WD)、降水量(RP)和大氣壓強(Press). 每個城市選取所有空氣質量監測站的平均值代表該城市的污染水平，選取以城市命名的氣象站點數據代表該城市的氣象條件.

1.3 研究方法

1.3.1 隨機森林算法和沙普利加和解釋

隨機森林是一種機器學習算法，常被用于預測PM2.5和O3等污染物濃度，具有較高的準確性和泛用性[16]. 該算法主要通過在訓練數據集上訓練多個決策樹來實現預測. 按照時間先后順序對歷史氣象數據進行隨機采樣，將采樣結果輸入數學模型中，預測O3濃度. 按照“隨機采樣-模型預測”的步驟重復完成1 000 次，將1 000 次的預測結果進行平均，得到氣象標準化后的O3濃度(O3_rmw). 氣象標準化的O3濃度(O3_rmw)和原始的O3濃度(O3_obs) 的差值就可以認為是氣象所引起的O3濃度(O3_met). 但隨機森林模型作為一種集成的機器學習方法，難以對模型中各輸入參數的影響進行解釋. 因此，本研究采用沙普利加和解釋方法以彌補上述缺陷. 沙普利加和解釋方法作為一種基于合作博弈論的機器學習解釋方法，可將預測結果合理分配到每個特征上，從而克服了機器學習可解釋性差的問題[33]，其具體原理如下：

式中：f為原始隨機森林模型； Φ0(f)為模型預測結果的平均值，是模型的基線值；M為模型特征的數量；Φi(f,x)為特征i的SHAP 值，是特征i對預測結果f(x) 的貢獻. 當 Φi(f,x) ＞0〔或 Φi(f,x)＜0〕時，意味著特征i會增加(或減少)f(x)， |Φi(f,x)|越高，表示特征i越重要；R是所有特征排序的集合，PiR是在排序R中排在特征i之前的所有特征的集合.

1.3.2 區域O3污染事件的定義

根據《環境空氣質量評價技術規范(試行)》(HJ 663-2013)[34]對山東省各監測站點O3小時濃度數據進行審核統計，基于審核后的有效數據計算O38 h滑動平均濃度，將O3日最大8 h 平均值 (MDA8 O3)超過160 μg/m3(GB 3095-2012)[31]視為O3超標日.本研究定義單個城市O3超標日持續3 d 及以上，且同時涉及城市超過8 個(山東半數以上城市)為一次污染事件.

2 結果與討論

2.1 氣象條件和排放對于O3 污染貢獻分析

2.1.1 氣象標準化

鑒于山東省行政區域面積廣闊及其半島特征，地理位置及經濟發展水平的差異可導致不同城市O3污染形成機制存在較大差異. 因此在進行氣象條件和排放貢獻剝離分析前，首先利用聚類分析將山東省各城市聚類. 聚類計算基于2020-2022 山東省大氣氧化劑Ox(即O3+NO2)觀測數據，應用K-means 聚類分析方法，將山東省16 市分為三個亞區域(見圖2)，分別為魯西(包括濟南市、德州市、聊城市、菏澤市、濟寧市、泰安市、淄博市和棗莊市)、魯中(包括濱州市、東營市、濰坊市、日照市和臨沂市)和魯東(青島市、煙臺市和威海市). 這一區域劃分也與Wang 等[35]研究的山東省O3污染空間分異結果相一致.

圖2 2020-2022 年山東省Ox 聚類結果空間分布Fig.2 Spatial distribution of the cluster of Ox in Shandong Province

按照緯度近似原則，在山東省魯西、魯中及魯東三個亞區域各選取一個典型城市-濟南市、濰坊市及煙臺市作為代表，對其2020-2022 年O3觀測值和氣象標準化值進行統計(見圖3). 從O3氣象標準化值看，三個亞區域O3前體物排放整體變化上呈逐年上升趨勢，年變化速率分別為3.3%(魯西)、7.8%(魯中)和4.4%(魯東). 2020 年氣象標準化值普遍偏低，也與當年前半年的新型冠狀病毒感染疫情影響有關；而后2021 年各地普遍復工復產，O3前體物排放明顯抬升. 從O3觀測值和氣象標準化值看，2020 年，魯西(濟南市)、魯中(濰坊市)、魯東(煙臺市)O3氣象標準化結果均分別高于觀測值(分別為8.26、6.14、6.60 μg/m3)，說明氣象條件有利于高濃度O3的消減，屬于“天幫忙”；2021 年三市氣象條件對3 年平均濃度的貢獻分別為-4.39、-4.74、-2.27 μg/m3，均略低于2020 年同期，說明氣象條件的有利效果有所減弱；而2022 年三地不利氣象條件影響分別為8.73、10.69、5.16 μg/m3，氣象條件有利于高濃度O3的生成.

圖3 2020-2022 年濟南市、濰坊市、煙臺市O3 氣象標準化結果與觀測值Fig.3 Meteorological standardization results and observed values of O3 in Jinan, Weifang and Yantai from 2020 to 2022

為進一步厘清不同氣象因素對于O3濃度的貢獻，選取山東省O3污染高發的4-9 月(即夏半年)，基于上述三個典型城市，分析2021-2022 年4-9 月山東省O3觀測值、氣象標準化值和氣象條件貢獻的平均濃度，結果如圖4 所示. 相較于2020 年上半年受新型冠狀病毒感染疫情較大的影響，2021 年和2022 年的月際變化更加客觀，因而選擇2021-2022 年來觀察各指標的月際變化. 從觀測值來看，各市各月O3觀測值多表現為2022 年高于2021 年，少數月份中兩年O3觀測值濃度相仿；8 月較為特殊，2022 年的O3觀測值月均濃度明顯下降，這一變化或與當年天氣形勢變化(諸如臺風等)有關[36]. 氣象標準化結果中，各市各月O3平均濃度差異不大，從年際變化看，2022年略高于2021 年. 從氣象條件貢獻看，2021-2022年各市氣象條件貢獻均為正值，說明在山東地區2021-2022 年夏半年氣象條件有利于O3的生成. 比較各月間氣象貢獻，發現濟南市(魯西)、濰坊市(魯中)氣象條件對O3濃度的貢獻最強，而在煙臺市(魯東)氣象條件對于O3濃度貢獻的最高值出現在5 月，這可能是因為在山東東部地區城市因沿海而有著更高的濕度水平，故而在更早的月份達到了氣象條件貢獻的峰值水平.

圖4 2021-2022 年濟南市、濰坊市、煙臺市的O3 觀測值、氣象標準化結果、氣象條件貢獻的月際變化Fig.4 Inter-monthly changes of the observed, meteorological standardized O3 levels and meteorological condition contribution to O3 in Jinan, Weifang and Yantai from 2021 to 2022

2.1.2 氣象條件影響歸因

不同氣象因子對于當地O3貢獻的SHAP 值如圖5 所示，可以看出在典型城市的所有年份中相對濕度(RH)和溫度(T) 都是影響O3小時濃度的主控因素，這與劉杰等[37]的研究結果一致. 溫度通過控制大氣中的化學反應過程可直接影響O3濃度[38]，也可通過調節氣孔的吸收強度影響O3干沉降過程[39]. 作為沿海城市的煙臺市，其風向(WD)對于O3小時濃度的貢獻較為顯著. 由各亞區域間對比可以發現，魯西、魯中地區氣象條件對于O3的影響水平相仿，且明顯大于魯東地區，這可能由于魯西及魯中地區前體物排放強度更高，在不利氣象條件下更容易導致O3的生成. 對比各亞區域3 年的結果可以發現，氣象條件對當地O3的整體貢獻逐年趨于正向影響. 2020-2021年00：00-09：00 氣象因子對當地O3小時濃度呈負貢獻，隨后在11：00-19：00 呈正向影響，20：00 后轉為負向；2022 年00：00-08：00 呈負貢獻，09：00 開始氣象因子對當地O3貢獻呈正向，且正貢獻持續時間增加. 為進一步對比不同城市間氣象條件影響的差異，統計了山東省16 市日間(07：00-17：00) 各氣象因子對O3小時濃度貢獻的平均值(見圖6). 整體來看，溫度(T) 和相對濕度(RH) 對O3小時濃度的貢獻最為重要；而部分沿海城市風向(WD)的貢獻比重明顯上升，這與沿海地區頻繁的海陸風影響相吻合[40].

圖5 2020-2022 年濟南市、濰坊市、煙臺市氣象條件對O3 小時濃度的影響Fig.5 Effect of meteorological conditions on O3 hourly concentration in Jinan, Weifang and Yantai during 2020-2022

圖6 2020-2022 年山東省16 市日間各氣象因子對于O3 濃度的貢獻Fig.6 Contribution of each meteorological factors to daytime O3 concentrations in 16 cities in Shandong Province

2.2 區域O3 污染事件中氣象條件和排放影響剖析

根據區域O3污染事件發生的定義進行統計分析，山東省2020-2022 年共經歷24 次區域O3污染事件，2020 年、2021 年、2022 年分別為9、7、8 次，全部出現在4-9 月. 區域污染期間，全省16 個城市MDA8 O3范圍為70～319 μg/m3，平均值為179.4 μg/m3. 對O3污染事件時段利用隨機森林的機器學習方法進行氣象標準化分析，并對各年度O3污染事件時段O3平均濃度的觀測值和氣象標準化值進行統計，結果如圖7 所示. 總體來看，在O3污染事件發生時，O3_obs明顯高于O3_rmw，說明不利氣象條件是區域O3污染發生的主要貢獻者. 從年度變化看，氣象條件影響整體呈現逐年上升的趨勢，特別是在2022 年，濟南市、濰坊市和煙臺市不利氣象條件貢獻分別達到59.19、51.93 和41.45 μg/m3，遠高于全年平均水平. 為進一步分析不同氣象因子的貢獻，對2020-2022 年O3污染事件時段不同氣象因子對于O3貢獻的SHAP 值進行統計(見圖8)，可以看出溫度(T)和相對濕度(RH)是影響O3小時濃度的主控因素，其中溫度(T)的影響最為顯著，濟南市、濰坊市和煙臺市三市日間峰值影響分別為89.61、66.62 和33.36 μg/m3. 這與王治非等[41]的研究結果一致，溫度在影響O3濃度的氣象因子中占據主導地位，且魯西大于魯東. 氣溫升高時O3前體物的排放增加，O3濃度也隨之增加[42]. 此外，對比3 年間結果可以發現，夜間氣象因子的負貢獻逐漸減弱，而日間的正貢獻逐漸增強，整體上對于當地O3貢獻逐年趨于正向影響. 比較不同區域，則發現不利氣象條件的影響在三個亞區域呈現魯西＞魯中＞魯東.

圖7 2020-2022 年濟南市、濰坊市、煙臺市在區域O3 污染事件期間氣象標準化結果與觀測值比較Fig.7 Comparison of meteorological standardization results and observed values during regional O3 pollution incidents in Jinan, Weifang and Yantai from 2020 to 2022

圖8 2020-2022 年濟南市、濰坊市、煙臺市區域O3 污染事件期間氣象條件對O3 濃度的影響Fig.8 Influence of meteorological conditions on O3 concentrations during regional O3 pollution events in Jinan, Weifang and Yantai from 2020 to 2022

3 結論

a) 隨機森林模型分析結果顯示，2020-2022 年山東省不利氣象條件是導致其O3污染的重要誘因.從氣象條件影響看，2020 年及2021 年氣象條件有利于高濃度O3的消減，屬于“天幫忙”；而2022 年氣象條件有利于高濃度O3的生成.

b) 相對濕度(RH) 和溫度(T) 都是影響O3小時濃度的主控因素，而作為沿海地區的魯東城市煙臺市，其風向(WD)對于O3小時濃度的影響較為突出.

c) 區域O3污染事件期間，排放的貢獻與非污染期間接近，而不利氣象條件是O3污染形成的主要誘因. 從單氣象因子看，區域O3污染事件期間溫度和相對濕度是影響O3濃度的主控氣象因素，其中溫度的影響最為顯著，對O3小時濃度的貢獻為33.36～89.61 μg/m3. 后續應更加重視在高溫、低濕條件下的O3污染防控工作.