近年来我国区域大气污染形势日益严峻，以细颗粒物(PM2.5)污染为关键污染物的区域灰霾污染形势尤其严峻(吴兑，2012；贺泓等，2013)，对人体健康以及人民日常生活带来很大不利影响。为了对细颗粒物污染进行动态监测，我国从2012年正式开始全国性PM2.5浓度业务化监测工作。在此之前，基于卫星遥感反演PM2.5浓度是了解我国PM2.5空间分布特征的主要手段(马国欣等，2008；van Donkelaar et al., 2010)，其主要缺点是时间分辨率较低，难以提供小时分辨率的连续观测信息。此外，反演过程中引入的偏差也使其存在较大局限性(Kahn et al., 2005；Lin et al., 2014)。因此，目前对重点灰霾污染区PM2.5污染的时空分布及其生消演变特征仍不清晰，给区域灰霾污染控制及其预报带来很大障碍。

全国PM2.5监测网络建立为了解我国PM2.5污染状态提供了良好数据基础，如何利用目前在空间分布上仍较为稀疏的观测站点数据评估PM2.5污染的时空演变特征，其潜力亟待挖掘。现代先进大气污染数值模型(如NAQPMS，CMAQ等)包含了较为完整的大气污染物理化学演变过程，能提供主要污染物在三维空间的分布数据，为评估污染演变特征提供了强有力的工具(Binkowski and Roselle, 2003；王自发等，2006)。大气化学资料同化技术则为观测系统和数值模式的融合提供了可能，这种技术基于观测和模拟数据的误差估计，对污染物浓度场进行最优状态估计，可以有效减小数值模式初始场输入数据误差，提高大气污染预报准确率。在国际上，一些国家已经开展了大气化学资料同化研究(Carmichael et al., 2008；Sandu and Chai, 2011)，采用同化算法主要包括最优插值(如Tombette et al., 2009)、三维变分(如Pagowski et al., 2014)、四维变分方法(Errera et al., 2008)、集合卡尔曼滤波(如Gaubert et al., 2014)等。在国内，崔应杰等(2006)利用最优插值方法将PM10(空气动力学当量直径小于等于10 μm的颗粒物，即可吸入颗粒物)、SO2和NO2观测资料进行了同化，有效提高城市尺度大气污染预报精度。Lin et al.(2008)和Niu et al.(2008)利用资料同化方法对沙尘暴预报进行了改进。白晓平等(2008)利用最优插值法和卡尔曼滤波法改进了南京地区NOx、SO2的数值预报。Tang et al.(2011)基于集合卡尔曼滤波算法开展了京津冀臭氧资料同化试验，并建立了臭氧资料集合同化系统。然而，利用同化数据对大气污染状态特别是我国PM2.5区域污染时空演变特征进行评估分析的研究还非常少，其与观测和模式数据的作用差异还不清晰。

因此，本研究将尝试利用地面观测、数值模拟以及资料同化得到的PM2.5数据，分别评估京津冀地区在一次灰霾污染过程中PM2.5污染的时空演变特征，系统比较分析三种数据在PM2.5区域污染回顾分析上的优缺点，为三种数据在PM2.5区域污染时空演变特征分析上的应用提供研究基础。

本次研究采用中国科学院大气物理研究所自主研发的嵌套网格空气质量预报模拟系统(NAQPMS)进行数值模拟，该系统已用于区域和城市尺度的污染物形成和输送机制的研究(Li et al., 2012)，也在业务化空气质量预报中得到了广泛的运用，并在北京奥运会，上海世博会，广州亚运会和南京青奥会等重大活动空气质量保障上发挥了重大作用(吴其重等，2010)。

NAQPMS模式系统包含污染物的动态污染源排放、平流、扩散、干湿沉降和气相、液相和非均相反应等物理化学过程，可同时模拟SO2、NOx、O3、CO、NH3、硫酸盐、硝酸盐、沙尘等多种污染物，同时还包含污染来源解析、敏感性分析和资料同化等技术(Wang et al., 2001)。模式为三维欧拉输送模型，垂直方向采用地形追随坐标。水平方向采用双向嵌套网格实现从区域到城市污染模拟的研究，并可以反映不同尺度污染物的相互作用。NAQPMS对物理、化学过程的处理均采用目前国际上较为成熟的方案：气相化学机制为CBM-Z碳键反应机制，包含了核心化学反应机制134个，可模拟71种污染物，液相化学采用RADM模式的反应机制，包括气体和气溶胶物种约22个，非均相化学过程包括不同类型的气溶胶表面的14个物种及28个反应(Li et al., 2011, 2012)。模式还包括起沙机制，可计算亚洲沙尘的传输。

模式区域设置为三层网格嵌套，如图 1a。第一区域覆盖东亚地区，水平格距45 km×45 km，网格数182×172；第二区域包含中国东部、韩国及日本，水平格距15 km×15 km，网格数210×198；第三区域为京津冀及周边地区，水平格距5 km×5 km，网格数219×213；垂直方向采用地形跟随坐标，不等距分为20层，模式顶层高度20 km。

模式模拟采用中尺度气象模式WRFv3.5为NAQPMS污染模拟提供逐时的气象场。区域排放源为清华大学提供的MEIC多尺度排放清单模型(http://www.meicmodel.org/[2013-07-05])，在此基础上利用河北环保局提供的污染源普查和排污申报数据对河北局地点源做了更新。模式的初始和边界条件来源于全球大气化学输送模式MOZART的模拟结果。

模拟分析时段为2014年冬季京津冀地区一次细颗粒物重污染事件(2014年2月18~28日)，分析区域为北京、天津和河北。模拟分析数据采用第三区域数据，分辨率为5 km。模拟的PM2.5包括硫酸盐、硝酸盐、铵盐、黑炭、有机物气溶胶、一次排放PM2.5等。

观测数据采用中国环境监测总站的全国城市环境空气自动监控系统数据，该系统覆盖全国113个城市956个站点，本次研究区域将用到京津冀区域13个城市80个监测站点的PM2.5浓度数据。各站点位置见图 1b。

本次研究将用到三种数据，观测空间插值数据，模式模拟数据和基于最优插值的同化数据。模式模拟数据由2.1介绍的NAQPMS模式直接积分算出。本节将简要介绍观测空间插值方法和最优插值方法。

为了了解PM2.5的空间分布特征，本文利用Cressman插值算法对站点观测的PM2.5浓度数据进行了空间插值。客观分析方法中的Cressman插值算法是将离散点内插到规则格点引起误差较小的一种逐步订正的内插方法，在气象领域中有广泛运用(张德二等，2003)。

Cressman插值算法采用逐步订正方法进行插值，用实际资料与初始猜测场之差来订正预备场，再用实际资料与新场之差来订正新场，直到订正的新场逼近实际资料为止。对于每一步订正有如下定量描述公式：

其中，V为任一气象要素场，V0为准备场，ΔV为观测资料与准备场的插值，V'为订正后的新场。Δvij为格点(i, j)上观测资料和猜测值的差：

其中，Δvk表示观测站点k的监测值与初始准备场的差值，wijk是权重因子，在0~1之间变化，K是影响半径R内的站点的数量。权重wijk的确定由下面的公式确定：

其中，影响半径R的一般是人为选取的一个常数，常用1、2、4、7和10，dijk是格点(i, j)到观测站点k的距离。

本文利用最优插值对观测数据和模拟数据进行融合，获得PM2.5浓度的二维空间场数据。最优插值方法(OI)是基于统计理论的一种分析方法，最优插值的分析方程为

分析场是在方差最小的意义下背景场和观测场的一种线性组合。其中，P为背景误差协方差矩阵，认为是静态的；R为观测误差协方差矩阵。利用所求增益矩阵K，以及背景场xb和观测场y直接求解分析场xa。一般还假定在决定分析增量的过程中只有少数的观测数据是重要的，也就是只同化模式点周围一定区域内的观测信息。观测误差协方差矩阵R的选取，观测的误差方差为δ2=(0.1Ytobs)2，即各点的观测误差均方差为该点观测值的15%，并且不考虑空间相关性，即R为对角矩阵。最优插值的背景场误差协方差矩阵是静态的，本文采用各向同性的背景场误差协方差矩阵，即不同模式网格间状态变量的误差相关性只依赖于网格间的空间距离，其水平相关性函数采用高斯分布函数：

其中，L是特征空间相关性尺度，根据模式分辨率和站点空间分布选定为125 km，即第三模拟区域的25个模式网格距；d是两个网格点间的距离；C是两个网格点间的误差协方差；V是模式网格数据模拟误差的方差，根据长时间观测模拟对比结果，将其设定为基准模拟浓度的80%。

通过观测数据Cressman空间插值、NAQPMS数值模拟以及模拟观测数据最优插值生成产生了三种类型的数据，分别是观测空间插值数据、模拟数据和最优插值分析数据。三种数据由于产生的方法不同，对污染事件表征能力的特点亦不同。下面将从空间和时间两方面，给出三类数据的分辨率和准确度的评估。

图 2给出了根据不同种类的数据得出的空间分布和在北京、天津、石家庄城市站点的时间序列的特征。图 2a为观测数据的插值分布，通过插值可以观察出污染物的空间分布特征，然而观测数据虽然精度较高，但基于观测数据的空间插值并没有考虑大气污染自身的各种物理化学特性约束，在观测没有覆盖地区，该插值结果存在很大不确定性。例如在京津冀西北部地区，由于观测站点稀疏，导致了虚假的低值区域(图 2a白色空洞)。图 2b为NAQPMS模式模拟数据，由于模式所采用的源数据、模式参数等仍存在较大不确定性，并且单个模式模拟不确定性较大，因此其模拟PM2.5浓度与观测值也存在较大差异，具有较大误差，对比图 2a和图 2b，可看出模式模拟的京津冀中部和西部区域的污染物相比观测插值明显偏高50~100 μg/m3。图 2c为同化数据的空间分布特征，由于同化数据能将观测在局地高精度信息以及模式的各种物理化学约束考虑进来，并利用局地观测信息对模式模拟结果进行较大范围的订正，使得同化结果在观测未覆盖地区精度也较高。对比图 2b和图 2c，同化后的PM2.5浓度场分布特征与观测插值场更加接近，而模拟场的浓度相对观测场有明显偏高。表明同化能很好修正模式模拟的偏差，能较好再现污染的整体分布特征。

图 2d-2f给出了三类数据在北京、天津、石家庄三个城市代表站点的事件序列，观测数据的时间序列观测值是代表实际的污染物时间变化，模拟数据在天津和石家庄的污染有明显的偏高，同化数据修正了模拟数据的偏高现象。由表 1模拟和同化数据与观测值的均方根误差和相关系数对比可以看出，除了廊坊站点外，京津冀区域各个站点模拟结果的数据精度在同化后得到了显著提升，大部分站点数据的均方根误差降至100 μg/m3以下，与观测值相关系数在同化前在0.6~0.8范围内，同化后大部分站点的相关系数提高到0.85以上。对于整个京津冀地区平均水平而言，基于最优插值的同化技术使均方根误差降低约47%，相关系数提高了0.22，说明同化资料显著提高了模拟数据与观测数据的吻合度，更好的反映了PM2.5浓度的实际的变化趋势和特征。需要指出的是，廊坊站点在同化后数据均方根误差大于同化前，主要是在个别时次同化对数据订正不恰当导致的，这可能与最优插值算法采用静态、各向同性的背景场误差协方差有关。表 2比较了观测数据，模拟数据和同化插值数据的各个特征，综合各个指标来看，同化数据兼顾精确度较高、数据连续性好且获取较为方便，相比观测数据和模拟数据有较大优势。

上述同化试验将所有站点观测数据纳入同化系统，为了进一步论证同化资料对于观测缺失站点污染数据重构的意义，我们专门设计了单独的交叉检验(cross-validation)试验，评估同化方法对观测没有覆盖地区数据精度的影响，其主要思路是利用没有在同化过程中使用的观测数据来评估同化数据的精度。

交叉验证(Cross Validation，简称CV)，有时候也称之为循环估计，是用来作分类器性能评估比较的一种统计分析方法，其基本思想是把某种意义下原始数据进行分组，用其中一组作为验证集，其他数据作训练集，首先用训练集对分类器进行训练，在用验证集对得到的模型做统计评价，得到分类器的性能指标。交叉验证应尽量满足均匀抽样，其一般的做法为随机抽样，当样本数量足够时，便可达到均匀抽样的效果。如样本数为N，分组数为G，按照分组的个数，交叉验证可分为2折验证(Hold-Out Method，G=2)，K折验证(K-fold cross-validation，G=K)，和留一验证(LOOCV，G=N)。为了减小计算量，本次试验采用2折验证，将京津冀80个站点数据随机分成2组(CV01、CV02)，每组40个站点，分别开展2次资料同化试验，每次将其中1组观测数据作为训练集进行同化，剩下1组观测数据作为验证集对同化效果进行评估。这种检验好处是所有观测数据都在独立验证中得到了应用，能更为全面反映再分析数据的精度。

表 3算出了每组交叉验证试验模拟数据和同化数据的均方根误差和相对系统偏差(relative error)。与4.1中的结论类似，相对于模拟数据，同化数据的均方根误差低近50%，为64~69 μg/m3。两组试验中，模拟数据与观测数据的平均相关系数都为0.65，同化数据的平均相关系数为0.9左右，后者比前者提高了36%。通过统计参数的比较可以认为，在缺少观测站点的区域，同化数据可明显提高对PM2.5的状态估计，其均方根误差和相关系数都明显优于简单的模式模拟结果。

前面的分析表明同化数据在空间连续性上能比观测提供更为丰富的信息，在数据精度上比模拟数据精度更高。针对本文研究的2014年2月19~28日京津冀地区重灰霾污染过程，张恒德等(2016)对这次污染过程的天气成因进行了详细分析，指出这次污染过程中京津冀及周边地区受弱高压脊影响，处于地面高压后部，地面受弱气压场控制，地面风速小，混合层高度低，存在逆温现象，大气扩散条件差，是这次京津冀重污染过程形成的主要气象成因。本文将不再详细分析这次污染过程的气象条件，而侧重利用最优插值同化数据来分析其时空演变特征，进一步了解这种数据在污染事件回顾分析上的潜力。

将京津冀地区80个站点PM2.5浓度数据和NAQPMS模拟的PM2.5数据进行最优插值的同化融合，利用同化数据分析污染事件期间地面PM2.5日均浓度的区域分布特征和演变。图 1是重污染期间PM2.5日均浓度空间分布变化，表 4是根据中国环境监测总站污染浓度等级标准计算出的各污染等级面积所占第三区域(包括京津冀和周边地区)面积的百分比。可以看到PM2.5在19日开始累积生成，在石家庄、邢台和邯郸地区有轻度污染，污染面积仅占6.3%。之后由于偏北气流的影响，污染气团移到京津冀南侧山东西部，20日早上开始由于气流转为偏南气流。污染气团由再山东输送到京津冀南部(石家庄、邢台、衡水、邯郸等)并在21日扩散到京津冀中部(北京、保定、廊坊等)，污染程度初步增强，重污染面积扩散至总区域的10%，石家庄和北京PM2.5浓度大于200 μg/m3(严重污染)。随后由于京津冀地区风速较小，天气静稳，污染物呈现夜间污染物累积加重，白天污染物清除减弱的日变化规律，轻度污染和重度污染区域明显扩大。23日夜间污染范围扩大至整个京津冀地区，全境达到重区域高达15%。24~26日是本次污染事件峰值，河北中部、南部，北京和天津PM2.5浓度均超过250 μg/m3(严重污染)，石家庄和北京污染浓度超过350 μg/m3。其中25日污染范围区域最大，发生污染的区域占总区域的60.5%；26日总污染区域稍有缩小，但污染严重的区域却进一步扩大，严重污染区域达6.6%。27日受强西北部冷空气南下，污染气团迅速南退衰减，27日凌晨北京PM2.5浓度降低到35 μg/m3以下，空气质量达优，27日傍晚污染气团被彻底清除，京津冀空气质量达优。至28日污染区域基本得到完全清除，清洁区域高达98.9%。整体看来本次污染的气象条件特征为：20日之前污染主要在山东地区累积，并通过南风传输到京津冀南部地区，20~26日由于京津冀静稳的天气条件造成污染物在京津冀进一步累计，偶有的较小的偏南风也加剧了污染物的扩散。27日之后的清除则主要是靠较大的偏北风造成。

图 3中可以看出污染有明显的由南向北的传输路径。图 4给出了京津冀地区由南向北8个城市的PM2.5同化浓度的时间序列。可以看出，这次污染过程可大致分为4个阶段，18日到20日为污染累计阶段，20日00:00(北京时间)，南部的邯郸、石家庄、邢台和保定最先达到中度污染；南部20~22日为污染局部扩散阶段，南部污染城市污染物持续累计升高，中部的保定、廊坊、北京的污染物也达到较高水平。22~27日是污染最严重的阶段，污染物全面扩散，除张家口以外所有的城市在该阶段的污染物平均浓度都在200μg/m3以上，并经历了相似的平滑累积到污染峰值爆发的过程，25~26日是大部分城市污染浓度最高的时段，由于污染气团向北输送，最南部的邯郸市在26日略有降低；随后在27日由于强北风入侵，北部城市(承德、北京、天津、廊坊)的污染物得到迅速清除，浓度迅速下降至60μg/m3以下，污染气团在南下的途中经过依次经过保定、邢台、石家庄、邯郸，使得这些城市的污染浓度依次经历短暂的升高然后被清除的过程。本次污染过程的传输特征为：污染物由南向北输送，由北向南清除。造成南部城市先污染后清除，北部城市后污染先清除。且有慢累积、快清除的特征。

图 5给出了基于同化数据估计的污染期间京津冀13个城市污染物浓度的平均值，它反映了整个污染事件期间，各个城市的污染程度。按照污染等级划分标准可将13各城市划分为4类。石家庄、邢台、衡水为重度污染；邯郸、保定、唐山、北京、廊坊为中度污染；沧州、天津、秦皇岛为轻度污染；承德、张家口为良。各城市污染等级的排序大致与城市由南到北的位置相对应。说明京津冀地区由北向南污染程度递增。

本次污染过程输送特征较为典型，其污染程度重，影响区域广。污染气团从18号开始累积，持续到27号凌晨才得到清除，影响了整个京津冀所有13个城市。其中张家口、邯郸、邢台、廊坊、保定等城市出现了持续多天的PM2.5浓度超过400 μg/m3的污染。

针对PM2.5污染事件的回顾分析和评估，本文通过不同方法得到了三类数据，分别是观测插值数据、模拟数据和同化数据，并利用时空分布图和交叉验证试验比较了各类数据的时空特征。结果认为：(1)观测数据具有较好时间连续性，但观测站点在空间上离散分布，采用Cressman空间插值构建填补没有观测数据覆盖的区域，插值数据在有观测站点附近精度较高，但在没有观测覆盖的区域误差较大甚至出现虚假的高低值中心。(2)模拟数据的时空分辨率高，但单模式预报数据仍存在较大不确定性，模拟误差较大。在本次模拟中，模拟数据在京津冀中部和西部区域的污染物相比观测存在明显偏高，各城市污染浓度与实际观测浓度的相关系数0.49~0.78，均方根误差大于100 μg/m3。(3)同化数据的时空分辨率与模拟数据相同，可以提供PM2.5浓度时空分布的动态演变数据。由于同化数据纳入了观测信息，并考虑到了大气中的物理化学过程，其数据精度比模拟数据更高。同化数据空间分布的量级和高低值中心与实际观测更加匹配，验证站点浓度随时间的变化趋势与观测数据相关系数最高可到0.9以上，大部分站点的均方根误差相比模拟数据有大幅降低。

同化数据较高的时空分辨率和数据精度，可以很好地重构观测缺失区域数据。本文利用同化数据描述并分析了2014年2月19~28日污染事件。并将这次污染过程大致分为累积、局地扩散、全面扩散、清除4个阶段。污染气象条件成因归结为偏南气流输送和本地静稳天气的共同作用，污染气团由南向北较慢地积累和传输，由北向南快速清除。该次污染过程演变过程典型，污染程度高，影响范围广。京津冀13个城市均受到了不同程度的影响，京津冀地区由北向南各城市污染程度递增，其中石家庄、邢台、衡水达到了重度污染。

值得指出的是，与同化数据在天气预报领域的广泛应用相比，大气污染同化资料在大气污染研究中的实际应用仍然非常有限。本文研究只是针对一次重污染过程的个例研究，通过比较分析了三种不同方法获得数据在大气污染事件评估上的优缺点。虽然这种比较揭示了同化资料在大气污染评估上的巨大应用潜力，但是目前研究仍有一定局限性。首先，采用同化方法为最优插值算法，这种算法具有计算效率高的优点，但其假设背景场误差协方差是静态的、各向同性的，未能有效考虑模式背景场误差的动态变化特征。因此，亟待发展和应用性能更高的同化方法，如能实现背景场误差协方差的流依赖动态估计的集合卡尔曼滤波和四维变分算法。其次，随着大气污染监测数据从少量稀缺开始向海量多源转变，地面、卫星、雷达等大气污染监测在空间上向立体多维发展，具有互补性的多源观测数据与数值模拟数据进行同化融合，有望为大气污染研究提供更高精度的三维大气污染同化数据，支撑污染事件分析、健康效果评估等关键研究。