基于多變量隨機(jī)森林模型的空氣污染過(guò)程分析及應(yīng)用

張玲玲，張亞一，盛 夏，章許云，吳 劍

（江蘇環(huán)保產(chǎn)業(yè)技術(shù)研究院股份公司，江蘇 南京 210019）

0 引言

空氣污染對(duì)人類(lèi)健康和社會(huì)發(fā)展均具有顯著的負(fù)面影響[1]。根據(jù)2020年環(huán)境公報(bào)，全國(guó)337個(gè)城市環(huán)境空氣質(zhì)量超標(biāo)率達(dá)43.3%，累計(jì)發(fā)生嚴(yán)重污染345 d，重度污染1 152 d，其中，以PM2．5，PM10，O3為首要污染物的天數(shù)分別占重度及以上污染天數(shù)的77.7%，22.0%和1.5%，空氣污染形勢(shì)嚴(yán)峻。

目前，空氣質(zhì)量觀(guān)測(cè)主要以地面監(jiān)測(cè)形式為主[2-3]，我國(guó)雖已建成覆蓋全國(guó)主要城市的監(jiān)測(cè)網(wǎng)，但受監(jiān)測(cè)站地理位置和數(shù)量的限制，無(wú)法對(duì)整體區(qū)域的污染程度進(jìn)行評(píng)估。相對(duì)而言，衛(wèi)星遙感技術(shù)不受地面監(jiān)測(cè)站點(diǎn)選址的限制，空間覆蓋范圍廣、時(shí)空分辨率高，但衛(wèi)星傳感器受云霧等不利天氣、設(shè)備故障等影響大，數(shù)據(jù)易缺失和出錯(cuò)，所獲得的數(shù)據(jù)與近地面實(shí)際差異較大，難以真實(shí)反映人類(lèi)活動(dòng)對(duì)環(huán)境的影響[4-5]。線(xiàn)性統(tǒng)計(jì)模型廣泛應(yīng)用于近地面NO2，O3和PM2.5等濃度的模擬，但該模型納入的參數(shù)較少，精度低，無(wú)法捕獲污染物濃度與更多參數(shù)之間的復(fù)雜關(guān)系。

隨著機(jī)器學(xué)習(xí)的快速發(fā)展，非線(xiàn)性機(jī)器學(xué)習(xí)模型被廣泛應(yīng)用，其中，隨機(jī)森林作為一種新型的非線(xiàn)性機(jī)器學(xué)習(xí)模型，可高效處理多變量和大數(shù)據(jù)量的問(wèn)題，具有運(yùn)算效率高、數(shù)據(jù)挖掘能力強(qiáng)及預(yù)測(cè)精度高等特點(diǎn)，被成功用于反演區(qū)域空氣污染物時(shí)空分布。游介文等[6]基于隨機(jī)森林算法、融合多源地理要素開(kāi)展了近地面NO2濃度空間分布模擬研究，模型精度高,月均模型整體擬合度R=0.85。MARM等[7]構(gòu)建了高性能隨機(jī)森林模型，以1 km的空間分辨率估算2005年～2017年中國(guó)PM2.5日平均濃度和O3日最大值8 h平均濃度，PM2.5的每日、每月和每年平均模型擬合R2值分別為0.85，0.88和0.90，O3～8h擬合R2值分別為0.77，0.77和0.69。目前，非線(xiàn)性機(jī)器學(xué)習(xí)模擬近地面污染物濃度的方式多針對(duì)大尺度區(qū)域，極少用于城市街道等小尺度的污染過(guò)程分析，本研究利用隨機(jī)森林模型，耦合氣象變量、地理變量、社會(huì)經(jīng)濟(jì)變量等，選擇南京市北部區(qū)域進(jìn)行空氣污染過(guò)程研究和實(shí)踐，為機(jī)器學(xué)習(xí)模型在城市街道污染防治工作中的應(yīng)用提供參考。

1 觀(guān)測(cè)區(qū)概況及研究方法

1.1 觀(guān)測(cè)區(qū)概況

南京市為重要的工業(yè)企業(yè)集中區(qū)，尤其長(zhǎng)江沿岸，環(huán)境污染指標(biāo)濃度較高，以南京北部六合區(qū)某區(qū)域?yàn)橛^(guān)測(cè)區(qū)，該區(qū)域位于南京、揚(yáng)州主城區(qū)之間，且鄰近多個(gè)工業(yè)集中區(qū)。監(jiān)測(cè)點(diǎn)位見(jiàn)圖1。

圖1 觀(guān)測(cè)區(qū)位置及對(duì)比區(qū)域空氣監(jiān)測(cè)站點(diǎn)分布

1.2 數(shù)據(jù)來(lái)源

本研究中空氣質(zhì)量數(shù)據(jù)包含國(guó)控站、省控站、區(qū)（縣）站點(diǎn)數(shù)據(jù)。其中，國(guó)控站數(shù)據(jù)來(lái)源于中國(guó)環(huán)境監(jiān)測(cè)總站的全國(guó)城市空氣質(zhì)量實(shí)時(shí)發(fā)布平臺(tái)（https://air.cnemc.cn:18007/），省控站數(shù)據(jù)來(lái)源于江蘇省環(huán)境監(jiān)測(cè)中心，區(qū)（縣）數(shù)據(jù)來(lái)源于南京大氣自動(dòng)監(jiān)測(cè)數(shù)據(jù)業(yè)務(wù)平臺(tái)。

氣象數(shù)據(jù)來(lái)自歐洲中期天氣預(yù)報(bào)中心（ECMWF）ERA5的0.25°×0.25°逐小時(shí)氣象數(shù)據(jù)，并綜合考慮污染物的產(chǎn)生、擴(kuò)散、消散機(jī)制及以往模擬研究的經(jīng)驗(yàn)，在數(shù)據(jù)可得的前提下，納入2 m露點(diǎn)溫度、2 m溫度、下行短波輻射、低云量、中云量、高云量、地表溫度、地表輻射、降水量、地面氣壓、100 m U風(fēng)分量、100 m V風(fēng)分量、10 m U風(fēng)分量、10 m V風(fēng)分量。同時(shí)，考慮到氣象對(duì)污染物濃度影響的滯后性，研究還納入了滯后1 h的氣象變量。

納入模型的地理類(lèi)參數(shù)包括高程、人口、GDP、夜間燈光數(shù)據(jù)，均來(lái)自于中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心（http://www.resdc.cn），使用2013年DMSP/OLS夜間燈光、2019年人口和GDP數(shù)據(jù)代表2021年情況。

1.3 模型建立

隨機(jī)森林算法是由BREIMAN L[8]在2001年提出的基于決策樹(shù)的高階機(jī)器學(xué)習(xí)算法，基本思想為構(gòu)造一定數(shù)量的決策樹(shù)并按照一定準(zhǔn)則對(duì)決策樹(shù)進(jìn)行組合生成隨機(jī)森林。由于多層隨機(jī)過(guò)程的存在，使得隨機(jī)森林可隨機(jī)生成幾百甚至上千棵決策樹(shù)，并保證每次構(gòu)建的決策樹(shù)均可能因隨機(jī)性而不同，可用于模擬多重非線(xiàn)性關(guān)系，構(gòu)成復(fù)雜的隨機(jī)森林模型[9]。

本研究對(duì)所有參數(shù)數(shù)據(jù)進(jìn)行預(yù)處理后，獲得2021年江蘇省5 km×5 km分辨率下逐小時(shí)參數(shù)數(shù)據(jù)集，隨機(jī)選擇其中90%的數(shù)據(jù)作為訓(xùn)練集，使用Python 2.7的scikit-learn庫(kù)基于十折交叉驗(yàn)證算法建立隨機(jī)森林模型，通過(guò)網(wǎng)格參數(shù)尋優(yōu)實(shí)現(xiàn)最優(yōu)參數(shù)的選擇。最終模型參數(shù)設(shè)置如下：最大深度（max_depth）設(shè)置為55，決策樹(shù)個(gè)數(shù)（n_estimators）設(shè)置為250，時(shí)間分辨率為每小時(shí)，空間分辨率為5 km×5 km。具體流程見(jiàn)圖2。

圖2 污染物時(shí)空分布模擬流程

1.4 模型驗(yàn)證方法

根據(jù)模型中的變量重要性排序，挑選重要性排序在前的變量建立最終模型，對(duì)未進(jìn)入模型建立過(guò)程10%的數(shù)據(jù)進(jìn)行模擬，通過(guò)線(xiàn)性擬合模型計(jì)算決定系數(shù)（R2）驗(yàn)證模型性能，并進(jìn)一步計(jì)算均方根誤差（RMSE）和平均絕對(duì)誤差（MAE）評(píng)估模型性能。

2 結(jié)果與討論

2.1 模型擬合和驗(yàn)證

隨機(jī)森林最終模型對(duì)驗(yàn)證集數(shù)據(jù)進(jìn)行模擬結(jié)果見(jiàn)圖3。由圖3可知，隨機(jī)森林模型估算的PM2.5，PM10，NO2，O3時(shí)均濃度與地面實(shí)測(cè)值吻合度整體較高，R2值分別達(dá)到0.82，0.85，0.77和0.89，其中，O3濃度模擬R2最優(yōu)，PM10濃度次之。

圖3 污染物時(shí)均濃度預(yù)測(cè)的隨機(jī)森林模型驗(yàn)證結(jié)果散點(diǎn)

隨機(jī)森林模型性能參數(shù)RMSE和MAE情況見(jiàn)表1。由表1可知，PM2.5，PM10，NO2，O3的RMSE值分別為10.68，23.01，9.38，14.29μg/m3，MAE值分別為6.25，11.50，6.27，9.95μg/m3。NO2和PM2.5的R2值雖然較PM10和O3低，但RMSE，MAE值偏小，模型的穩(wěn)定性更高。

表1 不同污染物預(yù)測(cè)的隨機(jī)森林主模型模擬精度

2.2 顆粒物污染過(guò)程分析

2021年3月29日～31日，南京市發(fā)生了連續(xù)重污染過(guò)程，期間，觀(guān)測(cè)區(qū)的空氣質(zhì)量指數(shù)（AQI）均值為172，共出現(xiàn)2次明顯的污染過(guò)程，首要污染物均為PM10，濃度變化情況見(jiàn)圖4。由圖4可知，3月29日6：00起，觀(guān)測(cè)區(qū)PM10質(zhì)量濃度開(kāi)始出現(xiàn)波動(dòng)，濃度數(shù)值呈上升趨勢(shì)，當(dāng)日21：00時(shí)達(dá)到395μg/m3后回落。3月30日3：00，PM10質(zhì)量濃度降至90μg/m3后隨即上升，期間質(zhì)量濃度最高達(dá)576μg/m3，污染過(guò)程持續(xù)至31日中午逐漸消退。PM10濃度高值期間，PM2.5濃度雖同步波動(dòng)，但整體穩(wěn)定，污染特征以沙塵污染為主。對(duì)比各區(qū)域的空氣環(huán)境監(jiān)測(cè)站點(diǎn)，觀(guān)測(cè)區(qū)與其它區(qū)域站點(diǎn)的PM10濃度變化趨勢(shì)基本一致。

圖4 觀(guān)測(cè)區(qū)與其它區(qū)域PM10濃度變化

基于隨機(jī)森林算法模擬南京及周邊區(qū)域的PM10濃度時(shí)空分布等值線(xiàn)結(jié)果見(jiàn)圖5。由圖5可知，2021年3月29日～31日污染整體表現(xiàn)為區(qū)域型。污染期間，區(qū)域的主導(dǎo)風(fēng)向?yàn)闁|北風(fēng)，位于觀(guān)測(cè)區(qū)上風(fēng)向的六合冶山站率先受到污染氣團(tuán)的傳輸影響，且PM10濃度在所有站點(diǎn)中最高，位于觀(guān)測(cè)區(qū)下風(fēng)向的江寧彩虹橋站受到的區(qū)域傳輸影響最小，該站點(diǎn)PM10濃度波動(dòng)相對(duì)滯后。潔凈空氣的輸入有利于降低站點(diǎn)的污染指標(biāo)濃度，30日3：00左右，風(fēng)向由東北風(fēng)轉(zhuǎn)為東風(fēng)，東側(cè)空氣相對(duì)潔凈，東風(fēng)環(huán)境下南京市各站點(diǎn)PM10濃度均大幅下降。但隨著污染擴(kuò)散和風(fēng)向的進(jìn)一步轉(zhuǎn)變，各站點(diǎn)又再次受污染氣團(tuán)影響，PM10濃度再次上升，且污染過(guò)程持續(xù)至31日中午，南京東北部站點(diǎn)空氣質(zhì)量率先轉(zhuǎn)好。

圖5 2021年3月29日～31日PM10濃度時(shí)空分布等值線(xiàn)

2.3 NO2污染過(guò)程分析

2021年1月12日～15日，觀(guān)測(cè)區(qū)的NO2濃度波動(dòng)明顯，變化情況見(jiàn)圖6。由圖6可知，共出現(xiàn)3次持續(xù)的濃度高值（1月12日17：00～13日9：00，1月13日19：00～14日11：00，1月14日18：00～15日9：00），且多集中在夜間。

圖6 2021年1月12日～15日NO2濃度變化

對(duì)比觀(guān)測(cè)區(qū)與其它區(qū)域站點(diǎn)NO2濃度變化，12日13：00起，觀(guān)測(cè)區(qū)NO2濃度不斷升高，而仙林大學(xué)城站同時(shí)段的NO2濃度數(shù)值波動(dòng)較弱，污染發(fā)生期間的污染物時(shí)空分布等值線(xiàn)見(jiàn)圖7。由圖7可知，觀(guān)測(cè)區(qū)NO2濃度主要受上風(fēng)向的主城區(qū)影響，而非大范圍的區(qū)域污染傳輸。13日15：00起，南京市多個(gè)站點(diǎn)NO2濃度變化趨勢(shì)較一致，均出現(xiàn)大幅上升，整體表現(xiàn)為區(qū)域型污染。14日12：00，隨著大氣擴(kuò)散條件好轉(zhuǎn)，區(qū)域NO2濃度下降，18：00左右再次迅速上升。日間，部分NO2參加光化學(xué)反應(yīng)被消耗，而夜間難以被轉(zhuǎn)化，且濕度增大、氣壓變低、風(fēng)速減弱等導(dǎo)致空氣擴(kuò)散條件變差，使得NO2積聚，濃度升高。15日12：00后，隨著污染消散，各站點(diǎn)濃度降至正常水平。

圖7 2021年1月12日～15日NO2濃度時(shí)空分布等值線(xiàn)

2021年2月6日～7日，受主城區(qū)及附近工業(yè)區(qū)共同影響，觀(guān)測(cè)區(qū)NO2夜間濃度較高，監(jiān)測(cè)結(jié)果見(jiàn)圖8。由圖8可知，6日18：00起，觀(guān)測(cè)區(qū)NO2濃度上升，濃度高值持續(xù)至7日12：00。對(duì)比觀(guān)測(cè)區(qū)與南京市、揚(yáng)州市部分站點(diǎn)，除北側(cè)的馬鞍街道外，草場(chǎng)門(mén)、奧體中心、揚(yáng)州儀征紅旗閘的NO2濃度數(shù)值變化趨勢(shì)較一致。

圖8 2021年2月6日～7日NO2濃度變化

2021年2月6日～7日NO2時(shí)空分布等值線(xiàn)見(jiàn)圖9。由圖9可知，6日夜間風(fēng)速較低，積聚的NO2消散慢，7日上午進(jìn)一步受早高峰、工業(yè)企業(yè)生產(chǎn)等影響，NO2濃度凌晨回落后有所上升，揚(yáng)州儀征紅旗閘站NO2質(zhì)量濃度升高至120μg/m3，遠(yuǎn)高于周邊站點(diǎn)，此時(shí)風(fēng)向?yàn)闁|北風(fēng)，污染可能傳輸至下風(fēng)向的觀(guān)測(cè)區(qū)及南京部分其它地區(qū)。

圖9 2021年2月6日～7日NO2濃度時(shí)空分布等值線(xiàn)

SO2濃度變化可側(cè)面表征工業(yè)生產(chǎn)現(xiàn)狀，揚(yáng)州儀征紅旗閘站鄰近工業(yè)集中區(qū)，不同區(qū)域的NO2和SO2濃度關(guān)系對(duì)比見(jiàn)圖10。由圖10可知，對(duì)比污染時(shí)段觀(guān)測(cè)區(qū)和上風(fēng)向揚(yáng)州儀征紅旗閘站SO2濃度變化，2個(gè)站點(diǎn)的SO2濃度變化趨勢(shì)較一致，且揚(yáng)州儀征紅旗閘站SO2濃度明顯高于觀(guān)測(cè)區(qū)，可進(jìn)一步明確7日上午觀(guān)測(cè)區(qū)NO2濃度升高為受工業(yè)排放影響。

圖10 不同區(qū)域的NO2和SO2濃度變化

3 結(jié)論

（1）基于多變量隨機(jī)森林模型的建立可較好地模擬近地面PM2.5，PM10，NO2，O3時(shí)均濃度，模型R2值分別達(dá)到0.82，0.85，0.77和0.89。NO2和PM2.5的R2值雖然較PM10和O3低，但RMSE，MAE值偏小，模型的穩(wěn)定性更高。

（2）2021年3月29日～31日，南京市發(fā)生了連續(xù)的重污染過(guò)程，首要污染物為PM10，通過(guò)隨機(jī)森林模型模擬江蘇省范圍PM10時(shí)均濃度時(shí)空分布更直觀(guān)地展現(xiàn)污染變化過(guò)程。污染期間主導(dǎo)風(fēng)向?yàn)闁|北風(fēng)，污染整體表現(xiàn)為區(qū)域型，觀(guān)測(cè)區(qū)站點(diǎn)PM10濃度隨污染擴(kuò)散和風(fēng)向轉(zhuǎn)變的影響顯著，污染過(guò)程持續(xù)至31日中午，南京東北部站點(diǎn)空氣質(zhì)量率先轉(zhuǎn)好。

（3）2021年1月12日～15日，觀(guān)測(cè)區(qū)的NO2濃度波動(dòng)明顯，共出現(xiàn)3次持續(xù)的濃度高值，且多集中在夜間，結(jié)合污染物時(shí)空分布模擬圖，觀(guān)測(cè)區(qū)NO2濃度除受夜間擴(kuò)散條件減弱等因素影響外，主要受上風(fēng)向主城區(qū)的污染傳輸。

（4）2021年2月6日～7日，污染物時(shí)空分布模擬圖顯示觀(guān)測(cè)區(qū)NO2濃度高值受主城區(qū)及附近工業(yè)區(qū)共同影響，尤其是7日上午，觀(guān)測(cè)區(qū)上風(fēng)向的揚(yáng)州儀征紅旗閘站NO2質(zhì)量濃度升高至120μg/m3，遠(yuǎn)高于周邊站點(diǎn)。當(dāng)時(shí)風(fēng)向?yàn)闁|北風(fēng)，污染傳輸至下風(fēng)向的觀(guān)測(cè)區(qū)，對(duì)比觀(guān)測(cè)區(qū)與揚(yáng)州儀征紅旗閘站SO2濃度波動(dòng)，進(jìn)一步確定污染源為工業(yè)排放。