魏颖 赵秀娟 张自银 徐敬 刘志权 孙溦 陈丹

1 北京城市气象研究院，北京 100089

2 美国大气研究中心，科罗拉多州博尔德 80301

3 中国气象科学研究院灾害天气国家重点实验室，北京 100081

1 引言

近年来，中国的大城市面临以细颗粒物（PM2.5）和高浓度臭氧（O3）为典型特征的大气复合污染问题（朱彤等, 2010; Zhang and Cao, 2015;Sun and Chen, 2017; Zhao et al., 2020）。大气颗粒物污染影响大气能见度，并通过气溶胶—辐射—云相互作用产生显著的气候效应和生态环境效应（Zhang et al., 2010; Boucher et al., 2014），同时影响人体健康（阚海东和陈秉衡, 2002; Liang et al.,2018）。当前国际主流大气化学模式对我国大气污染模拟预报存在明显偏差，秋冬季重污染期间PM2.5 浓度的预报偏差高达30%～50%（Zheng et al., 2015; Gao et al., 2016; Wu et al., 2020）。大气化学模式的模拟结果与观测存在较大偏差，其不确定性有多方面来源，从模式对污染形成和传输的表达角度，主要来自于以下3 个方面：1）排放清单的不确定性（Zheng et al., 2009; 曹国良等, 2011; 魏巍等, 2011）；2）输入气象场的不确定性，尤其风场、湿度、降水、温度等气象因子对大气化学过程模拟会产生显著影响（Zhang, 2008; Li et al., 2020）；3）模式理化过程方案的不确定性（王颖等, 2010; 张颖等, 2016）。除了以上与污染形成和传输等过程密切相关的3 个方面之外，数值模式技术本身的局限，如空间分辨率和理化方案适用尺度的匹配等，也会影响模拟结果的可靠性。减小空气质量预报不确定性对我国当前大气污染防控具有非常重要的现实意义和科学价值，在多方面不确定性来源混杂而无法快速诊断并提升模式性能的情况下，基于观测约束的资料同化是提高预报效果的核心技术。

早期限制大气污染资料同化的主要因素是观测资料的有限性（朱江等, 2018），近年随着观测数据增加，世界各国已开展了基于不同方法的大气污染资料同化方法研究：1）较为简便易行的最优插值法，如应用于中国上海地区NO2、SO2及PM10的三维分析场预报（崔应杰等, 2006）、欧洲地区PM10 的总浓度预报（Tombette et al., 2009）等；2）基于集合卡尔曼滤波系统的初始场、排放源或二者的同时同化，如基于嵌套网格空气质量预报模式系统（Nested Air Quality Prediction Modeling System,NAQPMS）开展的地面观测对O3及其前体物的初始场和排放源同时同化（Tang et al., 2010），Miyazaki et al.（2012）同化NO2柱浓度资料得到全球NOx排放源，基于天气大气化学模式（Weather Research and Forecasting model coupled with Chemistry, WRF-Chem）开展对PM2.5 初始场和排放源的同时同化预报（Peng et al., 2017），基于同化卫星资料的气溶胶同化（Dai et al., 2019）等；3）基于全局优化的三维变分同化技术，如中国气象科学研究院基于区域环境气象数值模式（Chinese Unified Atmospheric Chemistry Environment,CUACE）开展的气溶胶同化研究（Niu et al., 2008），Pagowski et al.（2010）对北美地区PM2.5 和O3浓度的同化预报试验，Li et al.（2013）同化洛杉矶地面PM2.5 及其组分质量浓度并进行验证以及美国大气研究中心开展的一系列基于格点统计插值（Gridpoint Statistical Interpolation, GSI）同化系统针对北美和中国地区的初始场同化工作（Liu et al., 2011; Jiang et al., 2013; Chen et al., 2019），汪代春等（2020）基于WRF-Chem 开展的中国气溶胶和气态前体物协同同化等；4）将空间全局最优扩展到空间和时间全局最优的四维变分技术，如Yumimoto et al.（2008）通过同化雷达消光系数反演起沙通量，Guerrette and Henze（2015）构建伴随模式实现对黑碳气溶胶的同化研究，Zhang et al.（2015）利用伴随模式进行了颗粒物污染源解析等。随着研究的开展，同化技术也开始应用于业务系统中，例如，中国环境监测总站的NO2、SO2及PM2.5 空气质量预报采用最优插值法同化，24 h 预报均方根误差能够减小23%（Zheng et al., 2018）。但总体来说在国内关于大气污染业务预报系统中的资料同化研究还较为缺乏。

在资料同化方法中，三维变分将最小方差估计问题转化为求代价函数最小值的问题（Sasaki,1970），可以考虑较为复杂和弱非线性的观测算子，且计算代价小，往往作为常规业务资料同化的首选方法（Niu et al., 2008; Zhang et al., 2011）。由于气溶胶多粒径段、多化学成分组成的特点，观测资料形式多样（如不同观测手段、不同粒径段和物种的气溶胶观测数据以及间接观测产品气溶胶光学厚度AOD），而不同模拟方案中的气溶胶的粒径/组分分类不同，因此面向不同的观测系统和不同的大气化学模式变量，必须建立不同的气溶胶资料同化系统（Benedetti et al., 2009; Pagowski et al., 2010; Li et al., 2013; 靳璐滨等, 2016; Zang et al., 2016），增加了同化系统开发和移植的技术难度。

为了应对多组分多粒径气溶胶表达在模式中的复杂性，近年国内外学者基于“一步法”（直接对所有的气溶胶模式变量实施同时同化），在GSI同化平台做了一系列工作，如Liu et al.（2011）针对GOCART 气溶胶方案，在GSI 同化平台上发展了对MODIS AOD 产品的三维变分同化算法，实现了对模式14 个化学变量的直接同化；Jiang et al.（2013）将WRF/Chem 中的GOCART 气溶胶方案引入GSI 同化平台进行地面PM10 观测资料同化；Schwartz et al.（2012）在同化MODIS 的AOD 基础上新增了地面PM2.5 总量的三维变分同化；Chen et al.（2019）在GSI 同化平台引入针对WRF/Chem中MOSAIC 气溶胶机制的同化方案，实现对PM2.5和PM10 包含的32 个分粒径段气溶胶变量的直接同化，并将其应用于中国冬季重污染成因分析。

为了提升大气污染预报和能见度预报的业务服务能力，并推动大气污染业务预报系统（RMAPSChem，WRF-Chem 内核）与气象业务预报系统（RMAPS-ST，WRF 内核，气象同化核心为WRFDA）的耦合，并减轻业务系统维护压力，北京城市气象研究院与美国大气研究中心合作，在RMAPS-ST所采用的气象同化平台WRFDA 基础上开发了与GSI 化学同化模块类似的化学同化平台（WRFDAChem），设计实现能灵活接入常用气溶胶方案，提高气溶胶初始场模拟准确性，通过同化对比试验指明不同同化方案对气溶胶模拟的影响，同时利用同化试验诊断发现模式对粗颗粒的显著低估与排放信息表达缺失有关以及模式常用的非均相反应摄取系数在新增同化前体物的情形下应予以新的调整，以促进模拟同化系统的协同发展（Sun et al., 2020）。本文基于RMAPS-Chem + WRFDA-Chem 同化预报系统，以2016 年秋冬季为例进行了为期1 个月的同化对照试验，研究WRFDA-Chem 同化地面观测资料对中国不同区域不同污染特点的PM2.5 和PM10预报的影响，同时对比分析了不同同化频次对同化效果的影响；在此基础上，结合模式本身的理化方案和偏差来源进一步探讨了在业务系统中开展同化预报的制约因素和发展方向，为大气污染业务预报和同化的相关研究提供参考。

2 模式与同化系统

2.1 RMAPS-Chem 系统介绍

睿图—化学预报系统（RMAPS-Chem）是基于北京区域环境气象数值预报系统（BREMPS）改进发展而来的业务模式系统，其采用全球通用的大气化学模式WRF/Chem 作为主积分模式，即中尺度气象模式（WRF）和大气化学模式（Chem）在线耦合同步积分计算，并由采用快速更新同化技术的睿图—短期预报系统（RMAPS-ST）提供气象背景场。RMAPS-Chem 采用CBMZ 气相化学机制和MOSAIC气溶胶机制（未开启高湿度下SO2-NO2-NO3通过快速非均相过程生成硫酸盐—硝酸盐—铵盐的过程表达）（Chen et al., 2016），经过多次升级，基于京津冀地区观测数据进行了能见度参数化方案改进和SO2排放清单优化（徐敬等, 2019; 赵秀娟等, 2019），近年来为北京重大活动和华北地区重霾过程和污染气象条件预报提供技术支撑。本研究所采用的RMAPSChem v2.0 是在RMAPS-Chem v1.0 的基础上进行配置更新和升级后的新版本，主要围绕扩大预报区域范围、保证预报时效的目的，在水平网格范围、垂直分层和模式版本等配置上进行升级，同时对排放清单和地表静态数据集做了更新。RMAPSChem v2.0 目前使用WRF v4.0.2 版本，采用兰勃特投影水平分辨率9/3 km 的嵌套网格，外层覆盖我国东部大部分地区（即研究区域，见图1），内层覆盖京津冀与汾渭平原地区。模式垂直分层为30 层，其中1500 m 以下为12 层，保障边界层内污染物排放、输送扩散和化学转化等过程的精细化模拟和污染物的可靠预报。静态排放源则引入了清华大学最新发布的0.1°分辨率的MEIC-2016 清单替代原来的业务清单（Zhang et al., 2019）。地表分类采用了基于2000～2010 年MODIS 遥感数据建立的Landuse 数据集，能更好地表征中国城市地表特征和陆地利用情况。在业务系统中，每天从00:00（协调世界时，下同）循环热启以提供96 h预报，初始场来自于前一天24 h 预报结果。

图1 研究区域范围及同化试验中采用观测站点（圆点）位置。方框为3 个分析区域：华北平原（NCP）、西北地区（NWC）、长三角地区（YRD）Fig. 1 Model domains and locations of the observational sites (dots).The regions defined in rectangles represent the three regions used in the analysis: North China Plain (NCP), Northwest China (NWC), and Yangtze River Delta (YRD)

2.2 WRFDA-Chem 同化系统

北京城市气象研究院联合美国大气研究中心（NCAR）在WRFDA 气象同化平台基础上开发了全新的多物种多粒径段气溶胶3D-VAR 同化系统WRFDA-Chem（Sun et al., 2020），采用和现有气象业务同化系统RMAPS-ST 完全相同的同化内核（WRFDA），简化了未来耦合系统研发和业务运行成本。在采用3D-VAR 同化初始场中，通过构造代价函数J(x)将最小方差估计问题转化为求代价函数最小值的问题，由此避免了直接求解权重矩阵，代价函数J(x)的表达式如下：

其中，x为分析变量，xb为背景场矢量；B为背景误差协方差矩阵，分析中决定背景场误差和观测误差权重分配；y为观测场矢量；R是观测误差协方差矩阵，H为观测算子。其中，对背景场误差协方差矩阵B的估计问题是资料同化的核心，它决定了同化时观测信息对分析变量（即各个组分/粒径变量）的订正范围和程度。WRFDA-Chem 同化系统通过水平、垂直和物理三方面的控制变量变换实现简化处理，并通过标准方差和水平/垂直相关尺度参数反映各变量之间的具体相关和结构；为B矩阵服务的气溶胶和气体物种变量的相关统计参数（Descombes et al., 2015），通过GEN_BE v2.0 模块计算，采用NMC 方法（Parrish and Derber, 1992）生成。本文基于研究时段2016 年11 月为期1 个月的WRF/Chem 预报，计算B统计量。

针对颗粒物的多种观测资料，WRFDA-Chem系统设计了多种同化方案以表达不同的观测算子，既可单独同化PM2.5 和PM10，也可同时同化二者；在同时同化的过程中，为了更好地表达粒径间关系，还设计了针对PM2.5 和PMcoarse（粗颗粒，PM10-2.5）观测算子的同时同化方案。在颗粒物观测资料外，WRFDA-Chem 中还加入了针对4 种常规气态污染物（SO2、NO2、CO、O3）观测资料的观测算子，能实现对颗粒物观测和气态污染物观测的同时同化，以更好地表达气体污染物对二次生成颗粒物的过程。为了适应不同气溶胶方案，WRFDAChem 设计采用注册表输入方式灵活切换分析变量；目前已接入了GOCART 和MOSAIC 两种气溶胶机制。

与RMAPS-Chem 匹配，本研究中同化过程所采用的气溶胶分析变量是适用于中国地区尤其是人为源大气污染状况的MOSAIC 机制，气溶胶成分包含硫酸盐（SO）、硝酸盐（N）、铵盐（N）、氯盐（NA）、钠盐（Cl）、黑碳（BC）、有机碳（OC）和其他无机成分（OIN）等8 个物种，每个物种采用0.039～0.1、0.1～1.0、1.0～2.5 和2.5～10 µm 4 个分段粒径，共计32 个变量。其中颗粒物同化过程中，PM2.5 总质量浓度同化影响前3 个粒径段8 种气溶胶物种，共计24 个变量（见公式2），PM10总质量浓度观测的影响则是通过PMcoarse 的观测算子，影响第4 个粒径段共计8 个变量（见公式3和4）。同化时对观测PM2.5质量浓度的分配采用直接同化的“一步法”，在考虑背景误差和观测误差的情况下，通过对背景误差协方差的统计，考虑多物种多粒径段组分为变量的整体代价函数最小化计算得来：

其中ρd为干空气密度。

在计算代价函数最小值时，对观测误差协方差的估计考虑仪器误差和观测代表性误差，其中观测代表性误差估计采用下式计算：

其中，γ 为比例因子，设定为0.5；ε0为仪器误差，为1.0+0.0075Mi，Mi为污染物浓度（单位：µg m-3）；Δx为网格分辨率，本文为9 km；L为观测影响半径，设定为2 km。上述参数根据Chen et al.（2019）的同化试验敏感性测试结果设定。

在利用地面观测数据进行大气污染同化时，单个观测点与模式网格面不匹配是在同化结果中引入误差的重要因素。当前对观测代表性误差的估计一般考虑网格分辨率、观测站点类别和物种等要素的影响（Elbern et al., 2007）。有研究利用京津冀地面观测网络对大气污染物地面观测代表误差进行了估计，结果表明网格内PM2.5 和PM10 的代表性误差分别为11.1 µg m-3和22.2 µg m-3（不确定性11%）（李飞等, 2019）。本文暂未针对观测代表性误差对同化结果的影响进行单独评估，将在未来的工作中考虑进行相关敏感性试验评估其影响。

3 背景误差协方差统计分析

气溶胶同化的难点就是利用有限散点的PM2.5地面总质量浓度资料对模式中多物种多粒径段气溶胶组分进行约束。如前文所述，是通过控制变量变换简化处理，通过水平、垂直和物理三方面的转化实现简化，通过标准方差和水平/垂直相关尺度参数反映各变量之间的具体相关和结构。我们采用了NMC 方法生成了背景误差协方差矩阵中简化的各统计参数。图2 为全区域平均的各变量的水平相关尺度参数，可以看到4 个粒径段内的气溶胶各组分，除量级较低的钠盐（Na+）和氯盐（Cl-）外，其他组分（硫酸盐SO、硝酸盐N、铵盐N、黑碳BC、有机碳OC、其他无机成分OIN）的水平相关尺度在地面附近大多在15～25 km 左右，结合模式网格9 km 的分辨率，以上水平相关尺度较为符合预期。在对流层内，各变量的水平相关尺度随着高度增加逐步减小。图3 给出了全区域平均的各变量的模式背景误差标准偏差，可以看到不同粒径段内不同组分的背景误差的标准偏差存在很大的差异。在构成PM2.5 的小粒径段内，硝酸盐具有最大的标准偏差，其次是铵盐，地面附近OIN 也具有较大的标准偏差，这反映了模式在没考虑新增非均相化学过程时，模拟区域中PM2.5 浓度的表达以硝酸盐成分为主的情况，而在大粒径段内OIN 组分的标准偏差量级显著大于其他组分，这是由于在大粒径段气溶胶成分以沙尘等自然源为主，二次无机盐组分占比较小。由于大部分排放源分布在地面层，气溶胶质量浓度在边界层以上逐渐降低，因而各粒径段各物种变量的背景误差标准差随着垂直高度增加而降低，在对流层顶部下降至接近0。上述气溶胶组分在不同粒径段的标准差的分布变化体现了中国中东部地区颗粒物由多组分多粒径段气溶胶构成的特点。

图2 2016 年11 月模式区域颗粒物粒径段（a）bin1、（b）bin2、（c）bin3 和（d）bin4 控制变量的水平相关尺度（单位：km）Fig. 2 Length scale of the aerosol species (units: km) of the (a) first size bin, (b) second size bin, (c) third size bin, and (d) fourth size bin over the domain in Nov 2016

图3 同图2，但为PM2.5 控制变量的背景误差标准差垂直廓线（单位：µg m-3）Fig. 3 Background error standard deviations (µg m-3) of the aerosol species of the (a) first size bin, (b) second size bin, (c) third size bin, and (d)fourth size bin over the domain in Nov 2016

4 资料与方法

4.1 模式设置与试验设计

本研究针对RMAPS-Chem 外层9 km×9 km区域进行同化试验，范围为（24°N～46°N，97°E ～122°E），覆盖中国大陆大部分区域，网格数为304×304（研究区域范围详见图1）。研究时间段为2016 年11 月，其中初始场模拟的预热时间为2016 年10 月24～31 日共计8 d。气象场初始条件和边界条件采用美国环境预报中心（NCEP）的1°（纬度）×1°（经度）FNL 全球再分析数据，化学场初始条件和边界条件采用WRF-Chem 提供的默认方案。模式采用的主要物理方案信息见表1，人为排放源采用清华大学0.1°（纬度）×0.1°（经度）的MEIC-2016 清单（Zhang et al., 2019）经插值匹配到RMAPS-Chem 的分辨率，包含甲烷、一氧化碳、二氧化硫（SO2）、氮氧化物（NOx）、挥发性有机物、氨、一次颗粒物（含元素碳、有机碳）等30 个物种。由于人为源相关的无组织扬尘排放不仅受自然条件影响，同时也受人为活动和减排管控措施影响较大，准确的排放总量统计和时空分配仍然存在较大困难，为了保证业务运行的一致性和稳定性，在业务模式中尚没有考虑人为源无组织扬尘排放。另外，由于现有模式的自然源沙尘模块中，采用了半经验的起沙参数化方案，起沙过程并未对不同起沙机制（空气动力直接拖曳作用、沙粒的跃移撞击、土壤聚合体的分裂）进行区分（张宏升和李晓岚, 2014），也未考虑不同粒径的沙尘粒子微观运动方式（蠕移、跃移、悬浮运动）的差异，对于起沙参数（如临界起沙摩擦速度、土壤湿度等）并没有开展很好的本地化验证，对垂直方向的抬升过程和分布表达也较为简易，容易造成虚假起沙和局地近地面层沙尘通量过大，从而带来业务模拟效果的不稳定，所以在业务系统中也未开启WRF-Chem 模式的沙尘自然源排放过程。

为了分析同化地面观测资料对颗粒物预报的影响，本研究设置了同化和不同化2 组对比试验，同化试验采用逐6 小时循环同化（DA_6h），另一组则采用不同化（NODA）作为对照试验。同时为了探究不同频次的同化对预报的影响，另外设计了逐3 小时循环同化（DA_3h）和逐1 小时循环同化（DA_1h）2 组试验。研究采用同时同化颗粒物PM2.5 和PMcoarse 的方案，并用同化得到的00:00的结果作为初始场进行24 h 预报试验，探讨同化初始场对于污染预报的作用，NODA 则采用原来预报循环热启的结果作为初始场，同化方案和试验设置详见表1。

表1 模拟采用的参数化方案和试验设计Table 1 Information regarding the parameterization and the detailed settings of the four experiments

4.2 观测资料与统计

本文观测数据来自中国环境监测总站发布的国控站常规6 项污染物浓度观测结果，覆盖研究区域共计836 个站点（分布见图3），观测时间段为2016 年11 月。为了保证观测资料的质量，同化前首先对观测数据进行了质量控制，包括极值控制和异常值剔除：对PM2.5 大于500 µg m-3，PM10 大于800 µg m-3的数据进行剔除；同时剔除结果小于0 以及PM10 观测小于PM2.5 的异常值。为保证同化的稳定运行，在同化过程中对与背景值（模式模拟值投射到观测空间）相比差异大于120 µg m-3的观测值也进行了剔除。

在定量化评估WRFDA-Chem 同化系统对预报的效果改善情况时，引入平均偏差（MB），归一化平均偏差（NMB），均方根误差（RMSE）和相关系数（R）等统计指标进行检验分析，且仅在观测有效数据大于总时次的2/3 时认为该站点为有效站点。统计时根据观测站点经纬度插值提取相近格点的模拟值，将区域内所有有效的单个站点的观测值和模拟值按时间序列依次接续作为样本计算相关的统计参数，作为区域平均值。

5 结果分析

5.1 分析场效果

对大气化学数值模式而言，初始场的准确性直接影响模式预报的效果，下面对同化前/后00:00 的初始场/分析场进行对比分析。图4a 和4b 分别给出了模式同化前的00:00 初始场（来自于前一天00:00 热启的24 h 预报结果）PM2.5 和PM10 月平均值以及与观测值的对比。如图所示，模式对我国北方大部分区域尤其是华北的PM2.5 浓度整体预报能力较好，京津冀—山东—河南污染带和汾渭平原地区的PM2.5 浓度在100～200 µg m-3之间，河西地区和京津冀北部城市PM2.5 浓度在30～100 µg m-3之间。然而模式对我国中南部地区的PM2.5 浓度存在明显的高估，例如对湖南、湖北地区高估超过50 µg m-3，对四川盆地的高估则超过100 µg m-3，可能反映了排放清单对这些区域的冬季排放量存在高估，或者模式本身对四川盆地的复杂地形和冬季污染表达仍然有较大不确定性。PM10方面，由于我国北方下垫面（地表粗糙度、可风蚀程度、土壤湿度等）和空气动力学（风速）条件比南方更利于起沙，而且降水相对较少，除了春季的长距离传输的自然沙尘天气过程，也易于产生小范围起沙和局地无组织扬尘，即使在秋冬季节我国北方观测的PM10 浓度也达到了150 µg m-3以上。而这部分小范围起沙和局地无组织扬尘，不仅受自然条件影响，同时也受人为活动和减排管控措施影响，其排放总量和时空分配都存在较大不确定性；目前在模式的静态排放源和起沙模块中都没有表达，因此模式对PM10 浓度的北方特异性很难准确把握。如图4d 所示，大部分北方地区模拟低估0～100%（约100～200 µg m-3），而中南部地区则高估20%～100%（约100～200 µg m-3）。由于人为源排放、自然源过程、以及模式本身理化过程的多方面不确定性，目前模式对于南北方的PM10 偏差存在方向性差异。

图4 2016 年11 月模式区域分析时刻（00:00）未同化的（a）PM2.5 和（b）PM10 月平均浓度（阴影）与地面站点（圆点）观测值对比、同化前00:00（c）PM2.5 和（d）PM10 月平均浓度（阴影）归一化平均偏差、同化前后（DA_6h-NODA）（e）PM2.5 和（f）PM10 月平均浓度（阴影）差异Fig. 4 Surface distributions of monthly mean (a) PM2.5 and (b) PM10 concentrations (shandings) at 0000UTC for NODA against the observations(circles), the normalized mean bias of (c) PM2.5 and (d) PM10 concentrations (shandings) at 0000UTC for NODA, and the difference of (e) PM2.5 and(f) PM10 monthly mean concentrations (shandings) between DA_6h and NODA over the domain in Nov 2016

经过DA_6h 同化得到的PM2.5 浓度00:00 分析场如图4e 所示，同化后的PM2.5 浓度在空间分布上得到调整，西北部分城市PM2.5 浓度提升20～50 µg m-3，中国南方大部分地区降低20～100 µg m-3，结合图4a 可见同化后的PM2.5 浓度更接近观测结果。图4f 表明，同化系统对PM10模拟的改善作用也同样明显，同化地面观测信息对我国北方地区PM10 浓度的提升可超过100 µg m-3，对南方PM10 浓度高估的问题也能降低50 µg m-3左右。图5 和表2 分别给出了同化前后分析时刻的观测—模拟散点对比和统计参数。如图5a 所示，同化后模式对不同站点的高估和低估现象都得到明显改善，DA_6h 将初始场的模拟偏差减少超过50%，NMB 和RMSE 也均减小一半左右。与此同时，分析时刻的PM2.5 浓度相关性也增加了0.27，说明WRFDA-Chem 同化系统对污染物在分析时刻的时空相关性能起到明显的调整作用。从图5b 的散点图来看，PM10 浓度的空间差异更大，同化后能够减小这种差异，尤其是对PM10 浓度低估的问题有明显改善。虽然由于区域内东西部区域的偏差方向不同导致同化前后区域平均偏差变化较小，但是同化过程能显著降低PM10 浓度均方根误差30%～40%左右，且相关性显著提升0.37。

表2 同化前后各试验分析时刻PM2.5 和PM10 浓度的模拟统计参数Table 2 Statistics between model results and observations for PM2.5 and PM10 for different experiments at 0000UTC

图5 2016 年11 月模式区域同化前后分析时刻（00:00）的模拟与观测（a）PM2.5 和（b）PM10 月平均浓度散点图Fig. 5 Scatter plots between the model results and the observations for monthly mean concentration of (a) PM2.5 and (b) PM10 at 0000UTC for the experiments of NODA and DA_6h over the domain in Nov 2016

需要指出的是，考虑到业务同化运行的稳定性，本研究中采用同化系统内外两重质控来筛选观测数据，质控过程的不同可能影响同化效果。第一重质控是在进同化系统前，对PM2.5 和PM10 观测分别大于500 和800 µg m-3的观测进行剔除；第二重质控是在同化过程中，对偏离背景值（模式模拟值投射到观测空间）超过120 µg m-3的观测进行了剔除。由于某些地区（如京津冀地区、山西省、汾渭平原）的PM2.5 模拟结果在污染过程严重偏低，带来部分观测数据在第二重质控过程中被剔除而无法用于同化过程，导致PM2.5 同化效果在这些地区不明显。而由于PM10 观测数据的同化是通过PMcoarse 影响第四个粒径段的变量来开展；对PMcoarse 而言，120 µg m-3的剔除标准已较为宽泛，剔除数据量较少，因而PM10 在京津冀地区的同化效果优于PM2.5。后期需要对双重质控，特别是第二重质控中的剔除标准开展进一步测试，以进一步提高PM2.5 的同化效果。

5.2 24 小时预报效果

模式预报的准确性不仅受到初始场的影响，同时也受排放源和气象场的不确定性影响，除此之外，模式的系统性误差也是影响预报结果的重要因素，因此随着预报时效的增加，同化初始场对于改善预报的正向作用会逐渐减弱。为了检验WRFDAChem 同化初始场对预报的改进效果和持续影响时间，采用6 h 循环同化前后获得的00:00 初始场进行了为期1 个月的逐日预报，并对0～24 h 预报的平均结果进行了对比。如图6 所示，24 h 以内WRFDAChem 同化系统对PM2.5 浓度预报偏差的调整能降低10～25 µg m-3左右，RMSE 能降低20～40 µg m-3，相关性提高～0.3。表3 给出了24 h 内平均统计参数，同化后预报偏差降低50%以上，RMSE 降低25%，相关性提升14%。相较PM2.5 浓度而言，同化对PM10 浓度预报偏差的改善相对弱一些。由于PM10 浓度模拟值在我国北方和南方的偏差方向相反，同化后整个研究区域的平均预报偏差反而有所增加，但同化后24 h 内平均RMSE 降低10%以上，且相关性明显提升，超过25%。总体而言，同化对PM2.5 和PM10 浓度预报效果的改善能持续24 h 以上，前12 小时改善效果更为明显，其中PM2.5 的RMSE 能降低30%～50%。

图6 同化前、后的整个模拟区域（a1-a3）PM2.5 浓度和（b1-b3）PM10 浓度0～24 h 预报统计参数比较Fig. 6 Averaged bias, RMSE, and correlation for the simulated concentrations of (a1-a3) PM2.5 and (b1-b3) PM10 against the surface observations in different experiments over the 0-24 forecast hours

表3 同化前后各试验PM2.5 和PM10 在0～24 h 预报中的统计参数Table 3 Statistics over the 0-24 forecast hours between model results and observations for PM2.5 and PM10 for different experiments

由于颗粒物浓度的空间分布差异较大，且观测站点东密西疏、空间分布不均，在同化系统迭代求解代价函数最小值进行全局最优化时，对不同地区的同化效果也可能存在较大的差异。如图7 所示，同化后华北平原、华中和四川盆地的PM2.5 浓度均方差显著降低约30～60 µg m-3，华中地区和预报偏差较大的西北地区PM10 浓度均方差降低约20～40 µg m-3。为进一步分析这一影响，本文对研究区域内较受关注的华北平原、长三角地区和西北地区的24 h 预报进行了统计分析（图8）。未同化时，华北平原PM2.5 浓度模拟偏差随着预报时效（0～24 h）增加经历了高估—低估—高估的变化，这反映了在00:00 至24:00 的日内变化过程中，原模式模拟在排放、气象场、化学过程等多方面产生偏差的过程中，可能出现了正/负两个不同方向的影响在不同时次存在协同和制约两种相互作用的现象。由于初始场偏差和中间时次预报场偏差相反，导致同化后的预报偏差在某些时次反而有所增大。而长三角和西北地区的PM2.5 浓度偏差方向单一，能看到图8b 和8c 中初始时刻的同化对减小0～24 h预报整体偏差能起到不同程度的正向作用，其中长三角地区的偏差改善更为明显，24 h 预报时效内偏差减小20%～50%，RMSE 减小12%～80%。对PM10而言，西北地区是重要的粗颗粒物源地，同化对于调整该地的PM10 浓度低估具有更为明显的效果。相较长三角地区来说，同化系统对我国北方的PM10浓度预报改善效果更好，偏差减小更明显，相关性提升的优势也能保持超过12 h。

图7 同化前（左）、DA_6h 同化后（右）24 h 预报与地面站点观测（a、b）PM2.5 和（c、d）PM10 浓度的均方根误差对比Fig. 7 RMSE maps at the surface validation sites for the concentration of (a, b) PM2.5 and (c, d) PM10 for 24-h forecasts of NODA (left column) and DA_6h (right column)

图8 同化前、后华北平原（左列）、长三角地区（中列）、西北地区（右列）（a1-a3、b1-b3、c1-c3）PM2.5 浓度和（d1-d3、e1-e3、f1-f3）PM10 浓度0～24 h 预报检验统计参数Fig. 8 Averaged bias, RMSEs, and correlation coefficients for the simulated concentrations of (a1-a3, b1-b3, c1-c3) PM2.5 and (d1-d3, e1-e3,f1-f3) PM10 in the North China Plain (left colunm), the Yangtze River Delta (middle panel), and the Northwest China (right panel) against the surface observations in different experiments over the 0-24 forecast hours

需要指出的是，在利用观测同化来更新初始场的过程中，其本身的物理意义是将模式在一段时间内来自多方面不确定性产生的累积偏差，在这一时刻予以强迫修正；而由于模式多方面不确定性（排放、气象场、化学过程）本身在预报过程中（这里体现为日内变化上）可能存在正/负两个不同方向的协同和制约作用，因而在某一时刻对浓度场的强迫修正，可能带来整体累积偏差在后续预报中反而增大的现象；在这种情况下，选择哪个时次对浓度场进行修正可能影响后续预报效果。整体而言，同化预报效果的提升其实与模式本身性能息息相关，需要在深入理解模式不确定性和偏差来源的情况下，开展模式和同化系统的协同发展。由于本研究主要目的是服务于业务系统，因而在试验设计上仅考虑了采用00:00 的同化分析场作为初始场来开展后续24 h 预报，而未考虑其他时次，未来需要开展更多的试验以优化同化预报流程。

5.3 同化频次对同化效果的影响

对3D-VAR 而言，同化频次也是影响同化效果的重要因素。为了探究RMAPS-Chem + WRFDAChem 同化预报的最佳方案，将3 h 和1 h 循环同化（DA_3h、DA_1h）与6 h 循环同化（DA_6h）的00:00 分析场及24 h 预报场进行了对比分析（见图6）。3 套循环同化试验在00:00 同化分析时刻的统计对比见表2，基于00:00 的初始场/同化分析场开展的24 h 预报结果见图6 和表3。对比表明，更高的同化频次能够通过逐步减小背景场与观测的差异从而进一步降低00:00 时刻PM2.5 分析场的偏差，但是PM10 分析场偏差变化不大，DA_1h 中PM2.5 和PM10 的相关性都得到进一步提升。随着同化频次提高对分析场的改进，PM2.5 的24 h 预报改进也较为明显，DA_1h 的24 h 预报偏差（MB、NMB）比DA_6h 减小了将近50%，相关性略提升0.01。PM10 方面，基于更高频次的同化初始场开展的24 h 预报相关性有所提高，DA_1h 相关性较DA_6h 提高0.02，但偏差变化不大，这是由于模式对PM10 模拟在我国南北方偏差方向明显相反，对整个区域模拟偏差的统计上可能存在抵消。而如图6 所示，相较于DA_6h，DA_1h 能够进一步降低PM2.5 和PM10 的均方差，对12 h 后的预报相关性也能进一步提高。因此在业务同化的应用中，采用逐小时循环同化高频次地引入观测信息，可以获得更为准确的初始场，从而有利于更长时效的提高预报准确性。

6 总结与讨论

为了服务于未来的气象—化学一体化业务系统，北京城市气象研究院与NCAR 合作，在WRFDA气象同化平台基础上新增了3D-VAR 大气化学三维变分同化系统WRFDA-Chem，实现了对多物种多粒径段气溶胶观测的同化功能。本研究将WRFDAChem 引入RMAPS-Chem 环境气象数值预报系统，利用地面颗粒物浓度观测资料开展了多组同化试验的对比分析研究，结论如下：

（1）利用WRFDA-Chem 同化地面PM2.5 和PM10 总质量浓度观测，对污染物分析场的时空相关性能起到明显的调整作用，DA_6h 同化后PM2.5和PM10 浓度相关性分别提升0.27 和0.37，同时有效减小污染物的空间分布偏差，同化后PM2.5浓度的平均偏差（MB）和归一化平均偏差（NMB）降低50%以上，RMSE（均方根误差）减小40%，对PM10 浓度的RMSE 也减少约40%。

（2）WRFDA-Chem 通过同化初始浓度场对颗粒物预报效果的改进能持续24 h 以上，PM2.5预报偏差降低50%以上，RMSE 降低25%，相关性提升14%，PM10 同化后24 h 内平均RMSE 降低10%以上，相关性提升超过25%。同化在前12小时改善效果更为明显，PM2.5 的RMSE 能降低30%～50%。

（3）在具有不同污染特征的区域，同化效果存在较为明显的差异，同化地面PM2.5 观测对长三角地区的PM2.5 偏差改善更为明显，24 h 预报时效内偏差减小20%～50%，但同化地面PM10 观测则是北方的预报效果提升更显著；另外，对于京津冀地区而言，原模拟偏差随预报时效在日变化上呈现正、负两个方向的差异，显示模式本身的排放、气象和理化过程等多方面不确定性作用存在协同和制约作用，因而在某个时刻引入同化后初始场可能会导致后续某些时刻的偏差反而增大。

（4）DA_1h 能够在DA_6h 基础上进一步降低PM2.5 和PM10 的均方差并提高相关性，根据Sun et al.（2020）和本研究分别采用45 km 和9 km分辨率在中国地区的WRFDA-Chem 同化试验，逐1 小时循环同化获得的分析场能够在不同程度上进一步优化预报效果并持续更长时间，对科学试验和构建再分析数据集而言，采用逐1 小时同化是更好的选择。但对于业务预报服务而言，需要同时考虑预报时效性和节约计算存储资源，采用6 h 循环同化获得的初始场启动预报已经可以获得较好的改善。

由于云—气溶胶—辐射相互作用，提高PM2.5/PM10 预报准确性能够有效减小气象要素的预报偏差，而对气象要素预报的改进也能反过来改进PM2.5/PM10 预报。Yang et al.（2020）将RMAPSChem 实时模拟的气溶胶浓度离线耦合到气象预报模式RMAPS-ST，对中国华北地区进行气象预报试验发现，采用实时模拟的气溶胶信息比采用模式默认AOD 固定值能够更好地预报2 m 气温和10 m风速，与不考虑气溶胶反馈相比，实时更新气溶胶信息能够比采用AOD 固定值将2 m 气温和10 m风速预报均方差分别进一步较低2.8%和0.3%。王哲等（2014）利用WRF-NAQPMS 双向耦合模式对京津冀秋季重霾的模拟表明，考虑气溶胶—化学双向反馈效应能够同时改善2 m 温度、2 m 相对湿度、10 m 风速和PM2.5 浓度模拟相关性和均方根误差。本文聚焦同化在颗粒物污染预报中的作用，没有考虑气溶胶反馈的影响，但在后期研究中，将利用RMAPS-Chem 气象—化学双向反馈的模式框架与WRFDA-Chem 同化技术的应用探讨同化气溶胶观测资料对气象预报和污染预报的影响。而在针对同化技术发展的未来工作中，一方面，需要对同化过程中观测数据质量控制环节进行更多测试，以优化同化效果；另一方面，同化预报效果的提升其实与模式本身性能息息相关，需要在深入理解模式不确定性和偏差来源的情况下，进一步开展模式和同化系统的协同发展。

致谢感谢审稿专家和编辑为论文改进提出的有益意见和建议，感谢中国环境监测总站CNEMC 观测网络提供的数据。