两种空气质量数值模式的应用评估与集合改进研究

孙弦，孙磊，聂会文，梁秀姬，苏烨康，王静，夏冬

(珠海市公共气象服务中心，广东珠海 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 的重要性基本相当，表明集合方法的预报能力与集合成员的线性偏差无关，主要取决于不同成员的代表性。

综上所述，本研究揭示以随机森林为代表的集合方法虽有效改进了污染物的预报能力，但提高数值模式自身能力和增加具有代表性的集合成员对预报水平的进一步提升非常关键。后续研究可以综合利用多种机器学习算法（如卷积神经网络），构建以多气象要素为主要自变量的空气质量统计预报模型，在评估其预报能力的基础上，将其作为成员进行集合预报，以期进一步提高珠海市（乃至大湾区）污染物预报能力。