兩種空氣質量數值模式的應用評估與集合改進研究

孫弦，孫磊，聶會文，梁秀姬，蘇燁康，王靜，夏冬

(珠海市公共氣象服務中心，廣東珠海 519000)

1 引 言

隨著城市化和工業化的不斷發展，近些年來我國區域大氣污染事件頻發[1]，其中以臭氧為代表的光化學污染事件[2]和以細顆粒物為代表的灰霾污染[3]最為突出，給人體健康、生態環境和氣候等方面造成很大威脅[4-6]。作為我國城市化進程最高的城市集群之一，珠三角地區的空氣污染呈現出明顯的區域性和復合型特點[7-8]。其中，廣東珠海作為粵港澳大灣區核心城市，社會經濟的飛速發展和人口的快速增長給空氣質量帶來嚴峻考驗。近幾年來一次污染防治雖有所成效，但以臭氧為代表的二次污染仍有加重趨勢[9]。因此，空氣質量的準確預報作為聯防聯控工作的基礎，對切實有效改善空氣質量具有重要意義。

目前，空氣質量預報方法主要分為人工研判、統計模型和數值模式三種[10]。其中人工研判依賴觀測數據和主觀判斷，在時間精度和空間尺度存在局限性，且不具備繼承性。統計模型雖運行操作簡單，整體預測精度較高[11-14]，但對一些極端污染事件的捕捉能力較差[15-16]，并缺乏可解釋性。另一方面，空氣質量數值模式利用數學方法定量描述污染物從排放、平流輸送、湍流擴散、化學反應到清除的完備過程，精細化模擬區域內污染物的時空分布特征與未來演變趨勢，已成為當前預報最主流的方法[17-20]。CMAQ、CAMx、WRF-CHEM和NAQPMS 等第三代空氣質量數值預報模式，自身基于“一個大氣”理念，將各類大氣問題、物理化學機制和相互作用統一考慮，已成為大氣污染科學研究的主要工具，并得到廣泛業務應用[21-24]。目前，分別基于CMAQ 和CAMx建立的華南區域大氣成分數值預報系統，通過使用我國自主研發的區域天氣模式作為氣象驅動場，并充分融合多套排放源，已經順利業務運行多年[21,25]。李婷苑等[26]評估了CMAQ 業務模式在廣東的模擬能力，但只重點關注PM2.5、O3及其前體物NO2（其他三種主要污染物PM10、SO2和CO 未評估），且模式在珠海的局地表現尚不清楚。此外，CAMx 作為平行運作的另一套模式，未有公開研究對其進行詳細評估。

不同數值模式對于不同空氣污染物在不同地區的預報效果存在明顯的差異[27-28]，因此，開展本地預報效果系統評估是業務應用的前提。空氣質量模式構成復雜，外部和內部皆具有較大不確定性[23,29]，外部主要包括模型運行所需的氣象初始場（包括土壤、地面和高空）、側邊界強迫和不同污染源排放清單等[30-31]，內部主要源自物理和化學過程參數化方案的不確定性[32-33]，使得預報結果必然存在一定程度的偏差。因此，利用數學統計方法對多個預報結果進行集合預報，對于衡量模式不確定性和提高預報能力具有關鍵作用[34-36]。

集合方法通常分為線性與非線性兩大類。其中，多元線性回歸因其構造簡單且考慮不同模式的權重，在研究應用中取得明顯改進效果。例如，潘錦秀等[37]利用多元線性回歸方法將CMAQ、CAMx 和NAQPMS 等三個模式進行集成，消除了單個模式系統性偏差，顯著提高了北京市2016 年PM2.5日均預報準確率。另一方面，以機器學習（例如BP 神經網絡、隨機森林和支持向量機等）為代表的非線性算法日益得到關注，但在空氣質量集合預報方面的應用卻不足。楊關盈等[38]綜合評估了多種集合方法對安徽地區PM2.5預報的改進，發現BP 神經網絡雖有一定訂正效果，但其效果卻不如多元回歸。但最近李娟等[39]卻揭示出相較于線性回歸，隨機森林和支持向量機方法對西安市O3和PM2.5預報的改進。湯靜等[40]采用主成分分析結合機器學習算法K 近鄰方法，有效地改進了CMAQ 模式對于廣州市PM2.5的預報水平。但需要說明的是，以上研究是針對單一預報模式，通過引入驅動氣象場進行回歸改進或者直接對模式預報進行后訂正，與多模式集合優化的思路有所區別。因此，評估檢驗以隨機森林為代表的機器學習方法在多模式集合方面的應用存在較高必要性。此外，以往研究多關注1～2 種污染物，并未實現對6項主要污染物的全面覆蓋。

綜上所述，本研究選取珠海市為研究對象，基于CAMx 和CMAQ 模型兩套獨立運行的空氣質量業務預報系統和國控點觀測數據，首先檢驗評估兩者對六項主要污染物的時空分布和演變特征的模擬能力，然后分別利用線性和非線性方法（即多元線性回歸和隨機森林方法）進行多模式集合，探究不同方法的改進能力，以期提高珠海市空氣質量預報水平，并為今后空氣質量多模式集合的研究與業務應用提供重要參考。

2 資料與方法

2.1 觀測與模式

本文選取珠海市四個國控站（環境空氣質量國控自動監測站，唐家、吉大、前山和斗門，圖1）為研究站點，收集2018—2019 年CO、PM2.5、PM10、O3、SO2和NO2等六種主要空氣污染物逐小時濃度觀測數據（缺測率約為20%），并在此基礎上計算不同時間尺度（日-月）的均值。需要說明的是，依據《環境空氣質量評價技術規范（試行）》（HJ663-2013）要求，O3日均值為當天8 小時滑動平均最大值（記為O3_8 h）。此外，珠海市平均污染物濃度近似認為是四個站點的平均。最后，2018—2019年珠海市空氣質量持續下行（年AQI 達標率均低于90%），所以被選為具體研究時段。

圖1 珠海市地形高度空間分布 其中紅星代表四個國控站（分別為唐家、吉大、前山和斗門）所在位置。

目前，中國氣象局廣州熱帶海洋氣象研究所和廣東省生態氣象中心分別基于CMAQ 和CAMx空氣質量模型，在華南區域構建了兩套大氣成分業務數值預報系統[21]，于每日08 時和20 時開始起報，預報未來72 小時逐小時空氣質量產品。兩者區域設置保持一致，水平為三重（27-9-3 km）單向嵌套，垂直分層數為25，并都使用國產自主高精度區域氣象模式CMA-GD（自身已同化多種實時氣象觀測）預報產品作為氣象輸入。對于排放清單，CMAQ 充分應用了清華大學的源清單、廣東EPA的珠三角排放源清單與中山大學的廣東交通排放源清單，并使用大氣成分衛星遙感資料和本地區地面站點觀測資料，對排放源分布和量級進行優化[21]。CAMx 的源清單也來自于多套源清單的融合，但并未進行觀測同化與人工訂正。兩個模式使用的物理化學方案也存在異同，主要設置詳見表1。

表1 CMAQ與CAMx模式設置

本文選取研究時段內（2018—2019年）兩套模式每日20 時起報的未來24 小時逐小時最內層（3 km）污染物濃度預報數據，并使用最臨近插值方法將模式格點數據插值到四個國控站點（圖1）以方便比較。

2.2 集合方法與實驗設計

（1）多元線性回歸。

多元線性回歸（multiple linear regression，MLR）方法通過將因變量Y（即集合預報）與多個自變量X1，X2，...，Xn（即多個模式預報）聯系起來，構建如下線性數學關系：

其中ai和b分別為回歸系數（可認為是第i個模式的權重系數）和回歸常數，可通過使用最小二乘估計進行求解。

（2）隨機森林。

隨機森林（random forest，RF）是一種監督學習算法[41]，由多個決策樹{y(x,θn),n= 1,2,……,N}組成的統計模型，其中θ為隨機變量（服從獨立分布），x為自變量，N為決策樹的數量。每一棵決策樹包含根節點、中間節點和葉節點，構建時首先在根節點進行分裂成各個分支，分裂過程需經過多個中間節點，最終達到樹的末端（即葉節點）為止。隨機森林里的每棵樹都利用訓練數據的子集（隨機選取樣本和特征）開展訓練，對于某一輸出規則，其輸出值是唯一的，最終輸出結果由各決策樹共同確定，因此具有不易過擬合、對異常值不敏感、解釋性強（可追溯），結果較為穩健等優點，因此廣泛應用于分類與回歸問題分析。對于在模式集合方面的應用，預測結果由各決策樹輸出值均值所確定，即：

其中，表示集合預報結果，y表示某一決策樹基于x和θ的輸出。此外，隨機森林是一種非參數算法，可以對每個輸入特征（即模型結果）相對于預測結果（即集合結果）的重要性（PIM，也稱為貢獻度）進行計算和排序。重要性基于袋外數據（outof-bag，OOB）計算，對于某一輸入特征，通過隨機置換（permute）輸入特征來計算該變化引起的平均準確度的下降（變化越大則該特征越重要），具體表達如下：

其中，i表示某一輸入特征，N為構造決策樹的數量，p表示置換后特征，MSE（mean square error）為均方誤差。

（3）實驗設計。

本文選取研究時段內模式預報與觀測分別作為兩種集合方法的輸入和輸出。為更好驗證集合方法的可靠性和泛化能力，本文采用5折交叉驗證法（5-fold cross validation）去開展模型訓練與測試。首先將2 年樣本劃分成5 個長度相等的樣本子集，然后依次遍歷5 個子集，每次選取其余所有樣本進行模型訓練，當前子集則作為測試集進行輸入驗證，最后合并5組驗證結果進行后續分析評估。集合模型基于不同污染物而獨立構建，并默認使用全部站點作為樣本數據。

2.3 評價指標

為定量評估兩個空氣質量模式及其集合方法的預報結果，本研究選取均方根誤差（Root Mean Square Error，RMSE）、相 關 系 數（Correlation Coefficient，R）和標準化平均偏差（Normalized Mean Bias，NMB）這三個統計檢驗，計算公式如下：

上式中，O代表觀測值，P代表預報值，N為樣本總數，為觀測值樣本平均，為預報值樣本平均。具體利用RMSE 來衡量預報準確程度，利用R來表明預報與觀測之間線性相關程度，以及利用NMB來反映預報系統偏差情況。

3 結果與分析

3.1 季節變化

首先，各污染物（除O3外）均呈現出明顯的冬高夏低特征（圖2），這與冬季化石燃料的加劇燃燒有關，而O3的產生主要依賴于光化學反應，因此高值出現在8—10 月。總體而言，CMAQ 模式較為合理地還原了各污染物季節變化，相關系數R介于0.72～0.84 之間，但存在明顯系統偏差，CO、PM2.5、PM10、SO2、O3和NO2的NMB 分 別 達 到-0.58、-0.18、-0.30、1.52，-0.16 和-0.20。CAMx 模式整體表現為顯著降低，各污染物相關系數均低于CMAQ（SO2甚至未通過0.05顯著性檢驗），低估了CO、PM10和NO2濃度（NMB 分別為-0.49、-0.53和-0.87），而對SO2則明顯高估（NMB為1.99）。需要注意的是，模式RMSE 和NMB 數值差異較大（特別是臭氧），這主要是NMB 在計算時進行了標準化(公式(6))，但正負偏差的相互抵消也對其NMB 的表現有所提升。例如，CAMx 整體低估了臭氧的平均濃度，但在2018 年11 月—2019 年2 月期間卻存在高估。

圖2 兩種數值模式（CAMx和CMAQ）及其集合方法（MLR和RF）2018—2019年珠海市六種空氣污染物濃度月均值變化與觀測（OBS）對比

通過多元線性回歸進行集合優化，CO、PM2.5、PM10、SO2、和NO2等污染要素的系統偏差得到有效糾正，NMB 降低到0.01～0.04，RMSE 分別降低到0.08 mg/m3、6.42、10.86、1.75和9.93 μg/m3，但在CMAQ 較好還原季節變化的基礎上，相關系數R無明顯改進，其中SO2相關性下降到0.71。更為重要的是，O3作為近幾年影響珠三角乃至全國最主要的污染物[8-9]，該方法對其季節變化的預報能力并未產生改進，RMSE 相較于CMAQ 模型，反而有所增加，這體現出線性方法的局限性。另一方面，非線性方法隨機森林表現明顯更為出色（表2），在其基礎上將各污染物（包括O3）的預報誤差RMSE進一步縮小到0.08 mg/m3、5.17、8.68、1.57、22.44和9.37 μg/m3，相關系數R提高到0.81、0.93、0.90、0.78、0.76 和0.78，這歸功于該方法基于集合算法（即基于多個獨立決策樹平均結果），準確性較單一算法（如多元線性回歸）有所提高[42]。另外，其在樣本和特征選擇時的雙隨機性，降低了模型產生過擬合的風險，使得研究時間段內表現均較為穩定。但是，包括隨機森林在內的兩種集合方法仍有缺陷，比如對O3和PM2.5高值月份的還原存在低估，這主要是因為樣本數量有限，未根據不同季節（或不同月份）對模型進行訓練所導致的，隨著模式和觀測數據的不斷積累，可在后續應用中得到優化。

表2 珠海市六種污染物季節變化統計參數

3.2 逐日變化

總體而言，CMAQ 對多數污染物日變化的預報能力都明顯優于CAMx（圖3）。對于CO，兩者表現接近，均可較好還原CO 的逐日變化趨勢（R為0.7 左右），但卻存在明顯系統性低估（NMB 分別 為-0.51 和-0.53）。CMAQ 不 但 有 效 減 輕 了CAMx對顆粒物的低估，PM2.5和PM10的NMB分別降低至-0.06 和-0.21，而且提高了年初污染天氣（即PM2.5日均值＞75或PM10日均值＞150）的捕捉能力，從而降低了預報誤差（RMSE分別降低了12.19和6.08 μg/m3），相關系數也得到提升。對于SO2，兩者表現均不理想，存在上述指出的嚴重正偏差，CMAQ 表現稍好，體現在演變趨勢的合理還原（R為0.55）。對于NO2，CMAQ 大幅糾正了CAMx 預報負偏差，NMB 從-0.88 提升至-0.12，但預報偏差仍較為明顯，RMSE 高達16.84 μg/m3。此外，NO2作為O3生成的前體物，CMAQ 對其模擬能力的改進，間接提高了O3的預報能力，O3相關性提高至0.56，預報偏差也降低了4.93 μg/m3，但對夏秋季易發的O3污染事件（即O3_8 h＞160 μg/m3）的捕捉能力仍有待加強[26]。

圖3 兩種數值模式（CAMx和CMAQ）及其集合方法（MLR和RF）2018年珠海市六種空氣污染物濃度日均值變化與觀測（OBS）對比

對于存在明顯系統偏差的污染物（即CO、SO2和NO2），多元線性回歸大幅糾正偏差，NMB 分別緩解至0.05、-0.09 和0.01，但SO2的相關性出現小幅降低。此外，該方法雖有效地提高了顆粒物統計評分，但對極端污染情況的還原能力卻不如CMAQ，這是由于顆粒物濃度在冬季明顯偏高，而模型基于所有時間段進行訓練，因此在該種情況下的表現受到了限制。最后，該方法對O3日變化的模擬未有改進，表現與CMAQ 基本相當。相較于線性回歸，隨機森林方法進一步提高了各污染物模擬的整體表現，各污染物的多項統計指標幾乎均為最優。另外，隨機森林同樣對冬季顆粒物污染事件還原能力有限，進一步驗證了利用所有季節樣本進行訓練的局限性。需要注意的是，臭氧作為近些年來珠三角空氣污染的首要威脅，隨機森林一定程度上彌補了線性方法的缺陷，不僅提高了其各項預報指標，而且加強了對極端污染事件的捕捉能力。

圖4 進一步給出了各要素逐日觀測與不同模式和集合方法的散點分布。CMAQ 雖明顯優于CAMx，但同樣對包括SO2、NO2在內的一些污染物存在明顯偏差，因此擬合斜率k距完美值1 差距較大。兩種集合方法明顯提高了各要素預報能力，尤其是隨機森林方法，各要素的擬合斜率k和決定系數R2都與完美值1 最為接近，展示出該模式優秀的集合預報能力。

圖4 兩種數值模式（CAMx和CMAQ）及其集合方法（MLR和RF）2018—2019年珠海市六種空氣污染物（a～f）濃度日均值（x軸）與對應觀測（y軸）對比散點圖（不同顏色代表不同模式或方法） 其中k和R2分別為擬合線的斜率和決定系數（兩者越接近于1，模擬效果越好，顏色與點相對應）。

3.3 晝夜變化

人為活動作為主要排放源，排放強度和類型具有明顯晝夜變化特征。并且，污染物的擴散活動主要受到大氣邊界層湍流活動的支配，而大氣邊界層高度也存在明顯晝夜變化[43]。因此，各污染要素也存在明顯的晝夜變化[44]。圖5 給出了模式和不同集合方法預報的各要素濃度晝夜變化（已減去自身均值）對比。據觀測，NO2晝夜變化為雙峰型外，其他污染物的日變化均為單峰型。總體而言，CAMx 模式幾乎無法還原各污染物的晝夜變化，出現明顯偏差，其中顆粒物和NO2的相關系數甚至為負，且CO、PM 和SO2均表現出類似的晝夜變化，揭示出排放清單的明顯缺陷。CMAQ 能較為準確還原O3晝夜變化（相關系數達到0.96），并大致表現出NO2的雙峰型特征，但對其他污染物的表現也不太理想，例如顯著高估了PM10和SO2的晝夜變化幅度，誤差分別達到13.18和9.5 μg/m3。另一方面，兩種集合模型對多數污染物（除CO 和SO2）晝夜變化并無明顯改進。這主要是由于集合方法均以減小誤差（如最小二乘法）為單一訓練目標，雖能有效減小模式的系統偏差，但未能對晝夜變化的還原產生附加價值。因此，污染物晝夜變化預報能力的改進主要在于模型自身的提高，并可嘗試在非線形算法中引入多目標函數進行多模式集合優化。

圖5 觀測（OBS）、兩種數值模式（CAMx和CMAQ）及其集合方法（MLR和RF）給出的珠海市2018—2019年六種空氣污染物濃度（已減去自身均值）晝夜變化對比

3.4 空間變化

排放源與氣象要素的空間差異，在擴散條件進一步作用下，各污染物要素呈現明顯的空間變化（圖6）。對于多數站點，PM 異常的符號與O3相反，這體現出兩者之間的“蹺蹺板”效應，即較高的PM 濃度削弱了太陽輻射，從而抑制了臭氧生成依賴的光化學反應。但PM10和O3在唐家站同為正異常，揭示了珠三角頻發的復合型污染[45]。總體而言，兩個模式合理還原珠海O3“東多西少”的空間特征，但對PM 和NO2空間差異的模擬卻存在明顯缺陷，這主要是由于氣象驅動模型GRAPES 能真實地模擬氣象條件（尤其是太陽輻射）的空間差異，為O3的生成與擴散提供了良好基礎，但排放清單由于空間分辨率和較大不確定性的限制，嚴重制約了PM 和NO2空間變化的模擬能力。同樣，基于所有站點樣本進行訓練的集合模型未能對空間差異的模擬產生效果。但以隨機森林方法為例，當基于不同站點構建模型，大幅改進了各污染物空間變化的預報水平。但是，空間技巧的提升也部分抑制了多尺度時間變化的還原能力（圖未展示），這同樣是由于訓練樣本長度不夠充分，因此無法支持模式基于不同維度（如不同季節和站點）開展優化。

圖6 兩種數值模式（CAMx和CMAQ）及其集合方法（MLR、RF和RF-sta）預報的吉大站（第1列）、斗門站（第2列）、前山站（第3列）、唐家站（第4列）2018—2019年四種主要空氣污染物年日濃度均值（減去站點平均，柱狀線，對應左側縱坐標，單位為μg/m3）、標準差（除以站點平均，三角形，對應右側縱坐標，單位為μg/m3）與實測對比

3.5 模式重要性

圖7進一步利用隨機森林模型的算法特點，展現了兩個模型對于不同要素重要性。以上分析表明CMAQ對于多數污染物的預報水平雖明顯優于CAMx，但兩者對于多數污染物的重要性未存在明顯差異，CMAQ 僅在O3方面展現出60%左右的較大優勢，而CAMx卻在CO預報方面占據明顯優勢，重要性達到64.6%。該結果揭示出模型自身的線性偏差對于隨機森林算法的結果并不產生影響[46]，另外的測試首先利用線性回歸對兩個模型進行誤差訂正，然后通過隨機森林進行訓練，其預報結果與未訂正相比也幾乎沒有差異。本研究僅使用兩個數值模型進行集成，因此，進一步提高集合預報結果的關鍵在于代表性集合成員的增加，而隨機森林多個獨立決策樹對特征的隨機選取，也極大程度上避免了過擬合發生，從而無需考慮集合成員過多對模擬結果產生負面影響。

圖7 隨機森林集合方法中CAMx和CMAQ模型對于各空氣污染物的重要性

4 結論與討論

本研究利用2018—2019 年國控站觀測資料，評估CAMx 和CMAQ 模式對珠海主要污染物時空分布與演變特征的預報能力，并引入多元線性回歸和隨機森林方法對預報結果進行集成，探究不同集合方法的改進能力。得出如下結論。

CMAQ 表現明顯優于CAMx，合理地還原了CO、PM2.5、PM10、SO2、O3和NO2的季節變化，相關系數介于0.72～0.84，但存在明顯系統偏差，NMB分別達到-0.58、-0.18、-0.30、1.52，-0.16 和-0.20，RMSE 分 別達 到0.40 mg/m3、6.86、16.02、10.71、25.05 和10.21 μg/m3。對于日變化，兩者對CO 和SO2技巧相當，但CMAQ 大幅修正了CAMx 模擬PM 和NO2的負偏差，提高了對冬季PM 污染事件的捕捉能力。由于對NO2預報的改進，CAMQ 提高了O3日變化的預報能力，相關性提升至0.56，預報偏差降低了4.93 μg/m3，但對夏秋季O3污染事件的預報能力存在不足。對于晝夜變化，CAMx 模式幾乎無法再現，而CMAQ 較為合理地還原了O3的晝夜變化（相關系數達到0.96），同時再現了NO2的雙峰型特征，但對其余污染要素存在明顯不足。并且，兩者對多數污染物（除O3之外）的晝夜和空間變化的模擬能力仍存在明顯缺陷，這主要來自于排放清單和氣象條件兩者的不確定性[23,29]。關于模式表現的差異，可以部分歸因于兩者基本架構和所使用參數化方案（如干沉降、氣象化學機理）[47]。此外，空氣質量模式的準確性依賴于合理精確的排放源清單數據[26]。CMAQ 所使用的排放清單在融合多種源清單的基礎上，進一步結合衛星遙感和觀測進行優化[21]，而CAMx 使用的的源清單則未經觀測同化和人工訂正，因此可以合理解釋CMAQ較優的預報能力。

基于不同污染物構建的兩種集合方法，均有效提高了季節-日尺度上的預報水平，其中隨機森林表現更優，對于各污染物的多項技巧評分幾乎均為最佳，但均對模式缺陷無明顯改進。這主要是由于線性模型為單個（或多個）輸入自變量和輸出因變量創建線性關系，但不同模型的結果通常是復雜的且具有高度非線性的關系。另一方面，隨機森林在解析非線性問題方面的優勢，配合在樣本和特征選擇時的雙隨機性，降低了模型產生過擬合的風險，因此展現出更為優秀的預報能力。但是，集合方法對污染物的晝夜與空間變化并無明顯改進，這表明集合預報雖具備優秀的附加價值，但預報水平受到集合成員預報能力制約。進一步基于不同地點對模型進行訓練，顯著提升了各污染物空間差異的還原能力，但其他方面表現受限于樣本長度而有所下降，這體現出集合方法對數據量的依賴性。隨著預報數據和觀測的積累，集合方法的實際應用中基于多維度（如季節和地點）展開較為必要。此外，隨機森林算法中CMAQ 與CAMx 的重要性基本相當，表明集合方法的預報能力與集合成員的線性偏差無關，主要取決于不同成員的代表性。

綜上所述，本研究揭示以隨機森林為代表的集合方法雖有效改進了污染物的預報能力，但提高數值模式自身能力和增加具有代表性的集合成員對預報水平的進一步提升非常關鍵。后續研究可以綜合利用多種機器學習算法（如卷積神經網絡），構建以多氣象要素為主要自變量的空氣質量統計預報模型，在評估其預報能力的基礎上，將其作為成員進行集合預報，以期進一步提高珠海市（乃至大灣區）污染物預報能力。