随着社会和经济的快速发展, 我国的大气环境污染日益严重, 并带来了显著的环境、生态和健康问题, 因而受到全社会的广泛关注, 大气污染研究也因此成为非常重要而必要的研究课题(上官铁梁等, 2000).局地大气污染或污染物排放在大气环流的作用下会产生输送和扩散, 使得大气环境污染表现出明显的区域性特征, 因此, 在关注某个地区的污染问题时, 不仅需要关注当地的排放, 还要考虑周边及更远地区的不同程度的贡献.

在过去的研究中, 有2种方法被经常用于分析大气污染的源区分布：① 采用气象和大气化学模式相结合的方法定量计算；② 使用轨迹模式在一定气象条件下定性估计.前者如李锋等(2015)利用WRF-CMAQ模型探究了长江三角洲在灰霾天气下大气细颗粒物(PM_2.5)的时空分布和输送, 并定量研究了外部源区及本地源对长江三角洲PM_2.5的贡献；王自发等(2008)利用嵌套网格空气质量预报模式系统(NAQPMS)模拟了北京及周边地区臭氧污染情况, 并利用质量追踪法计算了周边地区对北京的贡献率；Fountoukis等(2011)利用三维化学传输模式PMCAMx研究了2008年5月欧洲的颗粒物的组成和质量浓度.后者如王爱平等(2014)利用潜在源贡献因子分析法PSCF定性分析了不同气团背景下黄山顶积聚模态颗粒物数浓度的潜在源区；王茜等(2013)为分析不同季节PM₁₀和NO₂的潜在源区分布概率特征, 利用HYSPLIT4模式结合全球气象数据确定了影响上海气流的后向轨迹；Karaca等(2009)利用HYSPLIT的后向轨迹来估量可能对土耳其伊斯坦布尔PM₁₀浓度有远距离输送影响的区域.

FLEXPART是一个拉格朗日粒子扩散模式, 自1996年被用于计算长距离中尺度的点源空气扩散问题以来, 越来越多地被应用于大气污染研究.Stohl等(2003)运用该模式验证了城市区域间的污染输送；Halse等(2013)利用FLEXPART对持久性有机污染物的长距离运输进行了预测, 并根据预报结果收集样本, 结果验证了预测结果的正确性；DeCarlo等(2010)在MILAGRO综合外场观测期间研究墨西哥高原有机气溶胶的来源和过程中运用该模式验证了在盆地地区, 火灾影响因素与生物质燃烧有机气溶胶质量浓度具有很好的相关性.这些工作证实了FLEXPART模式对于污染物来源的模拟具有较高的可靠性.同时, Foy等(2009)比较了使用WRF和MM5两个区域气象模式结果驱动FLEXPART的效果, 结果显示, WRF气象场的效果更好, 因此, 目前多用WRF模式模拟结果作为驱动FLEXPART的气象场.

为解决气象与大气化学模式定量计算大气污染物来源的方法存在耗费大量计算资源及单纯气象条件分析又只能得到定性结果的不足, 本文利用FLEXPART模式、WRF模式和清华大学MEIC人为排放源清单, 建立一种快速定量估计某个目标地区大气污染物源区时空分布的方法, 并通过WRF-Chem模式模拟证实该方法的有效性.

FLEXPART模式是由挪威大气研究所(Norwegian Institute for Air Research, NILU)开发的一种拉格朗日(Lagrangian)粒子扩散模式, 它通过计算排放源释放的粒子的轨迹来描述示踪物在大气中的传输、扩散、干湿沉降和衰减等过程, 其中, 粒子排放源可以是点、线、面或体积源.该模式支持前向轨迹和后向轨迹计算, 前向轨迹可以模拟示踪物释放后随着时间的变化, 刻画出自源区释放后随气象条件的输送和扩散情况, 常用于判断特殊状况发生(如大气污染物泄露)后污染物演变情况；应用后向运算则可以确定对目标点或区域有影响的潜在源区的分布.FLEXPART的一个特点是描述一定气象条件下大气污染物的源与汇的关系, 目标区域为“汇”, 类似于受体, 大气污染物排放区域则为“源”(李岩等, 2010).事实上, FLEXPART模拟得到的是潜在源区, 即印痕(蔡旭晖, 2008), 而不是真正的大气污染物的源区.

WRF-Chem模式是一个当前广泛应用的气象和大气化学在线耦合模式(Skamarock et al., 2005；Grell et al., 2005), 由美国多所研究部门的学者共同参与研发而成.该模式包含了详细的大气物理和化学过程, 其化学部分(Chem)和物理部分(WRF)在空间和时间上完全在线耦合, 使用相同的模式网格、平流、对流和扩散方案, 计算化学时避免了对气象场插值而带来的误差, 因而得到广泛运用(Gao et al., 2016；周广强等, 2015).

清华大学MEIC人为排放源清单(http://www.meicmodel.org/)基于700多种人为排放源资料研制而成, 涵盖了10种主要大气污染物和温室气体, 分辨率为0.25°, 它按照电力、交通、工业、民用和农业5个行业进行逐月编制, 可以适应不同化学机制的需要, 应用广泛.本文使用基于2010年排放情景编制的版本, 常炉予等(2016)利用该清单研究了上海持续性重污染天气过程, 具有较好的模拟效果.

利用WRF模式模拟得到气象场, 利用该结果驱动FLEXPART模式后向运行, 得到针对目标区域和目标时间段内的大气污染物的潜在源区时空分布.在此基础上, 根据目标污染物的特性, 确定示踪物, 使其能够较好地反映目标污染物在大气中的生成和传输过程.以该前体物的排放强度(排放率)作为权重, 进而计算得到任意时刻和地点的“源”(所求源)对目标区域在目标时段内的选定大气污染物(目标物)的贡献率, 具体算法见式(1), 并进而获得某个时段和范围内的排放源对目标区域和目标时段的累计贡献率(式(2)).

(1)

(2)

式中, (i, j)和t分别代表某个空间位置和时刻；P为所求源在该时刻和位置上的敏感性系数；E为对应的代表示踪物排放强度；t1、t2为所求时段起止时刻；-∞代表研究目标的最早时间；(N1, S1) 和(N2, S2) 分别表示累计贡献的起始和终止位置, 当为(0, 0) 和(N, S)且t∈(-∞, 0]时, R=1；(N, S)为模式区域最大位置.E为1时, 其结果为潜在贡献率.

本文以上海地区2015年12月22—23日一次PM_2.5污染过程为例, 验证了前文建立的快速定量估计算法.此次污染时段涵盖了从较清洁到较重污染的不同程度PM_2.5状况, 污染过程较完整, 有利于凸显算法.徐家汇、闵行、宝山、崇明、浦东、金山、佘山和临港8个观测站平均结果(图 1)显示, 上海地区在此过程中, 22日4：00(本文未特别指明均为北京时间)PM_2.5浓度较低, 为40.7 μg·m^-3；之后迅速增大, 至23日2：00达到最大浓度213.1 μg·m^-3；之后快速下降, 至23日22：00降至29.3 μg·m^-3, 污染过程结束.22日20：00—23日12：00期间为PM_2.5高浓度时段, 平均浓度达到179.9 μg·m^-3.

图 1(Fig. 1)

图 1 20日16：00—24日16:00期间上海PM2.5浓度模拟值与观测值对比 Fig. 1 Comparison of simulated and observed PM2.5 concentrations in Shanghai from 16:00, 20 December to 16:00, 24 December 2015

模拟采用WRF模式3.2版, 其初始和边界条件使用1°×1°分辨率NCEP FNL全球再分析气象场资料.模拟区域以(31.5°N, 118°E)为中心, 水平网格数为360(东西)×400(南北), 分辨率为6 km, 垂直层数28层；时间步长为30 s.其主要的物理过程为YSU边界层物理方案、Monin_Obukhov近地面层方案、Unified Noah陆面过程方案、RRTM长波和Dudhia短波辐射方案及WSM6微物理方案(周广强等, 2015).模拟积分时段为2015年12月19日8：00—24日20：00, 模式结果逐小时输出.

与观测资料的对比表明, 模拟结果具有较好的可靠性.图 2显示了贡献率最大时刻22日20：00近地面(10 m)、925 hPa及850 hPa 3个层次的风场对比, 结果显示, 模拟风场和观测接近, 风速和风向整体上有较好的一致性.其中, 近地面风由于受下垫面的影响较大, 而且风向和风速是瞬时值, 变化剧烈, 模拟结果偏差较大；925和850 hPa 2个层次的一致性更好, 特别是上海及其临近区域.因此, 模拟的气象场可以用于污染物输送的模拟和分析.

图 2(Fig. 2)

图 2 22日20：00近地面10 m(a)、925 hPa(b)和850 hPa(c)模拟与观测风速风向对比 Fig. 2 Comparison of simulated and observed wind at 10 m (a), 925 hPa (b), and 850 hPa (c), at 20:00, 22 December 2015

FLEXPART采用后向模拟, 模拟时段为20日0：00—23日12：00, 垂直方向分5、10、100、500和1000 m 5层；取22日20：00—23日12：00作为粒子释放时间段, 上海部分区域(31.0°~31.4°N, 121.2°~121.7°E, 以下简称本地)、离地高度0~100 m作为粒子的释放区域, 连续释放粒子10000个, 离子在释放区域和时段内均匀分布.Zhou等(2014)研究发现, 上海在PM_2.5污染时硫酸盐、硝酸盐和铵盐占质量浓度的41.6%, 在重污染时铵盐浓度较之前不污染时增大10倍；同时, 一般区域大气化学模式不考虑金属离子, 铵根离子是其中唯一的阳离子, 因此, 本文选择铵盐的前体物NH₃作为PM_2.5示踪物, 其排放强度分布如图 3所示.考虑到人为排放源的高度, 选择1~3层(100 m)高度内的敏感性系数参与定量估计计算.

图 3(Fig. 3)

图 3 12月NH3排放强度分布(20：00) Fig. 3 Distribution of NH3 emission rate at 20:00 in December

首先分析了贡献率的时间变化.如图 4所示, 22日20：00—23日12：00, 本地和外地(本地之外所有区域)对上海PM_2.5污染的贡献总体相当；目标时段内, 本地的贡献略高于外地, 而由于输送作用的存在, 目标时段之前则是外地的贡献略高.针对本次过程, 22日10：00之前外地的贡献基本为0(图 4).22日之前上海主要受东北风的影响, 气团来源于海上, 因此, 22日之前的外地贡献几乎没有.

图 4(Fig. 4)

图 4 目标时段内的本地、外地及总小时贡献率随时间变化 Fig. 4 Temporal variation of local, non-local, and total contribution rates to the targeted period

空间分布上, 图 5a、c、e分别显示了23日12：00起后向依次前推12 h的贡献率分布, 反映了气象输送的踪影, 发现上海以东海洋区域也存在明显的贡献, 且随着时间的延长, 贡献区域范围增大, 最远延伸至朝鲜半岛；高贡献率区域主要集中在本地及临近地区, 远离本地的贡献率很低；主要气流方向上的贡献率明显高于其侧向上.加入示踪物NH₃排放强度作为权重后, 贡献率的分布发生了显著的变化(图 5b、d、f), 由于海上没有人为NH₃排放, 贡献率全部集中在本地和临近地区, 远距离输送的贡献极少.综合贡献率的时间变化和空间分布, 造成本次目标区域和目标时段内PM_2.5污染的主要来源为上海及临近的江苏地区.

图 5(Fig. 5)

图 5 23日0：00—12：00 (a)、22日12：00—23日12：00 (c)、22日0：00—23日12：00 (e)单位网格(6 km×6 km)潜在贡献率与贡献率(b、d、f)分布 Fig. 5 Distribution of potential contribution rate during 0:00—12:00, 23 December(a), 12:00, 22 December—12:00, 23 December (c), and 0:00, 22 December—12:00, 23 December (e) and the corresponding contribution rate (b, d, and f) for the same periods

为验证定量估计方法的可靠性, 使用WRF-Chem模式进行了包含完整大气化学过程的模拟试验.WRF-Chem试验采用与前文WRF试验相同的框架和物理方案.大气化学部分, 时间步长为60 s, 气溶胶化学为SORGAM有机气溶胶机制和ISORROPIA动力平衡无机气溶胶机制, 气相化学为RADM2机制, 生物排放源由MEGAN2模式在线计算, 化学侧边界条件采用包含了日变化的MOZART全球模拟月平均值.

首先开展了基本试验, 使用完整的人为排放源数据, 其在上海8个站点的PM_2.5质量浓度模拟结果如图 1所示.模拟结果与观测值相比有较明显的低估, 其中, 22日20：00—23日12：00平均PM_2.5浓度为126.1 μg·m^-3, 比观测值(179.9 μg·m^-3)偏低29.9%；21日16：00—23日12：00的平均PM_2.5浓度为82.0 μg·m^-3, 比观测值(119.4 μg·m^-3)偏低31.3%, 与周广强等(2016)得出的WRF-Chem模式倾向于低估污染时PM_2.5浓度的结论一致.但模拟的PM_2.5浓度的变化趋势与观测值具有非常好的一致性, 21日0：00—23日23：00的相关系数达到0.95, 合理地模拟出了该时段内PM_2.5浓度的变化规律, 说明其可以用于本文相关定量数值试验.

为验证WRF-FLEXPART定量估计结果的可靠性, 设计了一组敏感性试验, 对选定源区和时段内的人为排放源进行“清零”(常炉予等, 2016).根据贡献率的相对大小, 确定了A区(31.0°~31.4°N, 121.2°~121.7°E)、B区(32.0°~33.0°N, 120.0°~122.0°E)和C区(29.5°~30.5°N, 120.0°~122.0°E)区域(图 5), 分别代表本地、有一定贡献区域和无贡献区, B和C区为随机选择；分别将22日8：00—23日7：00和21日8：00—23日7：00内(表 1)的所有人为排放源调整为0.图 6显示了基本试验中上海8个站点平均PM_2.5浓度及22日8：00—23日7：00清零试验前后PM_2.5浓度的变化.A区清零后, 22日8：00上海PM_2.5浓度开始变化, 于22：00达到最大, 较基本试验偏低5.3%；23日11：00后PM_2.5浓度变化趋于零.B区清零后, 22日8：00上海PM_2.5浓度开始变化, 于22：00达到最大, 较基本试验偏低5.0%；23日11：00后PM_2.5浓度变化趋于零.C区域清零后, 上海PM_2.5浓度与基本试验相同, 无变化.

图 6(Fig. 6)

图 6 清零试验前后PM2.5浓度对比 Fig. 6 Simulated PM2.5 concentration and the relative changes of "zero-out" experiments

表 1(Table 1)

表 1 WRF-FLEXPART与WRF-Chem对目标区域PM2.5来源的定量计算结果对比 Table 1 Calculated contribution rates to targeted PM2.5 concentration by WRF-FLEXPART and WRF-Chem

表 1 WRF-FLEXPART与WRF-Chem对目标区域PM2.5来源的定量计算结果对比 Table 1 Calculated contribution rates to targeted PM2.5 concentration by WRF-FLEXPART and WRF-Chem

WRF-FLEXPART和WRF-Chem试验结果(表 1)对比分析表明, 2种方法得到的不同区域和时段对上海PM_2.5的贡献率有良好的一致性, 高、中、低贡献率的分布趋势相同；与WRF-Chem结果相比, WRF-FLEXPART方法低估了本地的贡献而高估了外来输送的作用, 而且差异显著.由于模拟期间发生了降水, 湿沉降对PM_2.5的清除作用也是一个重要因素, 因此, 对WRF-FLEXPART进行了湿沉降的敏感性试验, 采用铵盐的湿沉降特性参数进行计算.增加湿沉降作用后, WRF-FLEXPART结果更接近WRF-Chem结果, 偏差显著降低：本地的贡献明显上升, 外来输送的贡献下降, 无贡献的C区域保持不变.湿沉降作用的存在使得外地污染物在输送过程中通过降水而被部分清除, 输送作用对上海PM_2.5的贡献降低；与此同时, 本地排放的污染物停留时间相对较短, 湿沉降的作用较小.2个因素综合作用下, 外地贡献率明显下降, 本地的贡献明显上升.虽然WRF-FLEXPART的模拟结果仍与WRF-Chem模拟结果存在一定偏差, 但已显著改善, 具有可比性.

1) 本文利用FLEXPART拉格朗日粒子扩散模式、WRF模式和清华大学MEIC人为排放源清单, 建立了一个快速定量估计大气污染物来源时空分布的方法.此方法弥补了气象条件分析只能定性估计而完整大气化学模式方法计算量过大的不足, 可以快速获得不同区域、不同时段的污染物排放对目标区域、目标时段大气污染物浓度的贡献比率.

2) 本研究建立的方法的计算结果具有较好的可用性, 与利用WRF-Chem模式进行的人为排放源“清零”试验结果相比, 该方法得到的不同源区和不同时段的贡献率在时空分布上具有良好的一致性, 贡献率在数值上存在一定偏差但具有可比性.

3) 本研究所建立方法的定量估计结果存在一定不确定性, 应用FLEXPART模式时, 示踪物的选择、传输过程参数的确定(如湿沉降特性参数)等都可能带来偏差, 需要根据研究对象进行合理的设置.