引言

随着全球经济和社会快速发展, 城市化和工业化强度不断加剧, 环境问题日益凸显, 世界上重要的城市均经历过空气污染事件。近年来, 我国和一些发展中国家由于其快速的经济发展和城市化进程, 空气污染问题也日趋严重, 尤其是颗粒物为主的污染已经成为当前世界发展中国家大城市主要的大气污染物之一(Lelieveld et al, 2015;Mukherjee and Agrawal, 2017), 而我国大城市以可吸入颗粒物PM₁₀与细颗粒物PM_2.5为特征的区域性大气环境问题尤为突出, 严重影响人体健康与社会经济发展(谢元博等, 2014;郭家瑜等, 2017;王继康等, 2017;安林昌等, 2018)。北京奥林匹克体育中心(以下简称奥体中心)作为北京2022年冬奥会主赛区之一, 其空气质量更是备受世界各国人员关注。

城市大气污染中, 颗粒物主要有本地和远距离输送两个来源, 徐祥德(2002)研究亦表明, 北京地区大气污染是本地排放与周边区域输送作用的共同结果。结合气象条件、地形条件和区域污染源分布情况对气象传输条件进行分析是空气质量研究的一个重要方面, 拉格朗日大气扩散轨迹模式结合轨迹统计方法(trajectory statistical methods, TSM)成为研究污染物远距离输送问题的常用手段之一。近年来, 王茜(2013)、刘超等(2017)、周沙等(2017)、葛跃等(2017)和瞿元昊等(2018)大多采用潜在源贡献函数(potential source contribution function, PSCF)和浓度权重轨迹分析(concentration weighted trajectory, CWT)两种轨迹统计方法对不同地区污染过程的颗粒物潜在源区开展了大量分析。花丛等(2016)采用轨迹分析、聚类方法和PSCF方法对北京2013—2014年秋、冬两季大气污染来源进行了研究, 其结果表明:北京地区秋、冬两季的典型污染路径主要为西南和偏东两条, 潜在贡献位置为山东西部、河北中南部及天津等地区；清洁路径为西北一条。王郭臣等(2016)采用PSCF和CWT方法对北京2015年11月26日至12月2日的严重污染过程分析表明, 此次过程内蒙古中西部、山西北部、河北和山东北部等地区对北京PM_2.5质量浓度贡献在150 μg·m^-3以上。从上述研究可以看出, 我国关于颗粒污染物来源TSM主要采用PSCF和CWT两种方法, 而常用的轨迹统计方法有四种(Brereton and Johnson, 2012), 在我国使用另两种方法分析颗粒物潜在源区鲜有研究。

为探讨北京近5年2月颗粒污染物潜在来源与大气污染物质浓度趋势变化情况, 本文以奥体中心为研究区域, 利用2013—2017年连续5年的2月空气质量监测数据与气象观测数据, 分析不同大气污染物的浓度变化情况, 同时结合轨迹模式与四种轨迹统计方法分析北京近5年颗粒物质量浓度在75 μg·m^-3以上时的轨迹路径与传输特征, 为2022年北京冬奥会空气质量预报提供技术参考, 并为冬奥期间北京地区大气污染防治和减排方案提供技术支持和决策依据。

1 数据来源

文中所用污染物(PM_2.5, PM₁₀, CO, NO₂, O₃和SO₂)数据均来自生态环境部(http://106.37.208.233:20035/), 气象数据来源于中国气象局全国气象观测系统及CMISS数据库；气象数据与污染物数据均为逐小时监测数据。

2 研究方法 2.1 HYSPLIT模型

HYSPLIT(Hybrid Single-Particle Lagrangian Integrated Trajectory)是美国海洋大气局(NOAA)空气资源实验室(ARL)开发的混合单粒子拉格朗日综合轨迹模式(Draxler and Hess, 1998), 主要用于大气污染物的快速预测与“溯源”服务, 可以使用多种数值天气模式预报和再分析场数据, 本文驱动HYSPLIT模型的气象数据为欧洲中期天气预报中心(ECMWF)一日四次全球0.75°×0.75°再分析数据(ERA_INTERIM)。

2.2 轨迹统计方法理论

TSM是基于受体模型的反向计算模型, 通过计算从受体位置开始的污染物后向轨迹来确定污染物释放的位置。截至目前, TSM已经广泛应用于污染物识别溯源试验, 这类方法关键优势是计算速度快及实现起来方便灵活。

TSM方法主要包含四种, 分别为:潜在来源贡献函数(PSCF)、浓度加权轨迹(CWT)、停留时间浓度加权(residence time weighted concentration, RTWC)和定量传输偏差分析(quantitative transport bias analysis, QTBA)。PSCF方法(Ashbaugh et al, 1985)是将从受体测点开始后的停留时间与条件概率场联系起来, 用于追踪痕量气体和颗粒物的来源和传输途径(Hsu et al, 2003;Lupu and Maenhaut, 2002;Cheng and Lin, 2001)。依据PSCF理论, Seibert et al(1994)在PSCF方法基础上通过对受体点测量到的污染物浓度对轨迹进行权重划分, 将浓度与轨迹进行线性计算, 从而发展了CWT方法。Stohl(1996)在CWT基础上进一步进行了改进形成了RTWC方法, 该方法可以通过迭代重新调整研究区域内各部分的相对重要性, 以此识别轨迹方向上的重要污染源区。Keeler and Samson(1989)为评估TSM方法中轨迹的不确定性发展了QTBA方法, 假设大气扩散对轨迹的不确定性在轨迹中心线附近呈正态分布, 且轨迹路径的不确定性随着轨迹长度的增加而增加。下面对四种方法分别进行介绍。

2.2.1 潜在源贡献函数方法

PSCF是一种简单的轨迹停留时间概率分析方法, 广泛应用于分析痕量气体和颗粒物的来源和传输路径, 其主要特点是计算简单、速率快, 只需要一次计算, 无需进行迭代, 其计算方法详见王爱平等(2014)。本文将模拟区域定在北半球, 网格水平分辨率为0.5°×0.5°, 设定受体点颗粒物质量浓度大于75 μg·m^-3的轨迹为重污染轨迹, 经过计算, 本研究中网格点内平均轨迹节点数为7.64, 故本文取为8。

2.2.2 浓度权重轨迹方法

CWT基于PSCF方法发展而来。该方法以受体点监测到的污染物浓度对轨迹进行权重划分, 将浓度与轨迹进行线性计算, 具有较高CWT值的格点意味着该格点对受体点污染物浓度具有较高贡献。CWT法虽然与PSCF方法类似, 但不需要对轨迹进行浓度划分, 因此高浓度轨迹与低浓度轨迹相比停留时间更长。假设轨迹总数为M, C_m是轨迹m在受体点测量的浓度值, t_{i, j, m}是轨迹m经过单元格点(i, j)的时间, 则CWT计算公式如下:

$ CW{T_{i, j}} = {\bar C_{i, j}} = \frac{1}{{\sum\limits_{m = 1}^M {{t_{i, j, m}}} }}\sum\limits_{m = 1}^M {{C_m}} {t_{i, j, m}} $

(1)

2.2.3 停留时间浓度加权方法

RTWC在CWT方法基础上发展而来, 对研究区域的CWT进行迭代计算再分配, 其主要思想是一条“干净”轨迹经过的网格单元不可能是主要污染物来源, 通过迭代调整权重系数, 采用收敛条件获得最终RTWC的空间分布。轨迹m被网格单元分为N_m段, X_{n, m}是轨迹m第n段的平均浓度, ${{{\bar X}_m}} $是轨迹m的平均浓度。第n段对轨迹m的浓度贡献(C_{n, m})为:

$ \begin{array}{*{20}{c}} {{C_{n, m}} = {C_m}\frac{{{X_{n, m}}{N_m}}}{{{N_m}}} = {C_m}\frac{{{X_{n, m}}}}{{{{\bar X}_m}}}}\\ {\sum\limits_{j = 1}^N {{X_{j, m}}} } \end{array} $

(2)

基于式(2)和式(1)可得到一个新的浓度场如下:

$ {{\bar C}_{i, j}} = \frac{1}{{\sum\limits_{m = 1}^M {\sum\limits_{n = 1}^{{N_m}} {{t_{i, j, n, m}}} } }}\sum\limits_{m = 1}^M {\sum\limits_{n = 1}^{{N_m}} {{C_{n, m}}} } {t_{i, j, n, m}} $

(3)

通过反复迭代计算式(3)获得更新的浓度场, 直到达到设定的收敛条件(通常为迭代次数限制或者前后两次平均浓度场差别小于0.5%)即可得到格点RTWC分布。

2.2.4 定量传输偏差分析方法

QTBA是为了评估后向轨迹不确定性而提出来的一种方法。一般认为轨迹路径的不确定性随着轨迹长度的增加而增加, 并假设大气扩散对轨迹的不确定性在轨迹中心线附近呈正态分布, 对于轨迹m后向积分T时间内, 其潜在传输系数如式(4)所示:

$ \begin{array}{c} \bar{T}_{m}(v)=\frac{1}{T} \int_{0}^{T} \frac{1}{2 \pi a^{2} t^{2}} \exp \left(-\frac{v}{2 a^{2} t^{2}}\right) \mathrm{d} t= \\ \frac{1}{2 a T \sqrt{2 \pi v}}\left[{erfc}\left(\frac{1}{T a} \sqrt{\frac{v}{2}}\right)\right] \end{array} $

(4)

式中:T是轨迹节点所在时刻到受体监测时刻的时间, a是大气扩散速度(单位:km·h^-1), 常取值为5.4 km·h^-1, v是轨迹节点所在位置距离受体监测点的距离平方,

一般使用大圆距离进行计算。与受体监测浓度结合, 可以得到简化的QTBA计算公式:

$ {{ QTBA}} \left({x, y|{x^\prime }, {y^\prime }} \right) = \frac{{\sum\limits_{m = 1}^M {{{\bar T}_m}} \left({x, y|{x_R}, {y_R}} \right){C_m}}}{{\sum\limits_{m = 1}^M {{{\bar T}_m}} \left({x, y|{x_R}, {y_R}} \right)}} $

(5)

式中, x, y是轨迹节点所在经纬度；x_R, y_R是受体监测点所在经纬度; M为轨迹总条数。

基于以上方法可知, QTBA与CWT方法类似是一种非迭代TSM方法。

3 结果分析 3.1 传输路径

综合HYSPLIT模式和TrajStat方法(Wang et al, 2009), 以奥体中心(39.982°N、116.397°E)为考察对象, 将奥体中心上空500 m高度处设为受体监测点, 选取这个高度作为起点是因为这个高度既可以反映大气边界层流场特性, 又可以减少下垫面对气流的影响(刘超等, 2017)。研究选取2013—2017年的2月轨迹作为计算对象, 每日四个时次(00、06、12和18 UTC), 轨迹后向运动72 h, 时间分辨率设置为1 h。研究对得到的轨迹进行聚类统计分析(图 1), 可以得到2月奥体中心污染物传输路径及其对应的污染物浓度(表 1)。

从图 1可以看出, 2月影响北京奥体中心地区污染物主要轨迹路径分为4类:西北路径、北向路径、偏东路径和偏南路径。结合表 1和环境空气质量标准(GB 3095—2012)可以看出:西北路径(轨迹1)在2月出现概率最高为55.85%, 其输送的PM_2.5和PM₁₀质量浓度分别在72.76和101.31 μg·m^-3, 均达到24 h平均值国家二级标准(PM₁₀为150 μg·m^-3, PM_2.5为75 μg·m^-3；中华人民共和国国家质量监督检验疫总局和中国国家标准化管理委员会, 2016), 对应空气质量良好。该路径主要受西风气流影响, 移动速度较快, 但气团移动路径大部分在蒙古国境内, 大部分为沙质土壤, 可能在合适的环流系统配合下带来沙尘天气(PM_2.5/PM₁₀比值小, 自然粉尘粒径大, 占比高), 故导致颗粒物浓度升高, 可认为该路径为奥体中心2月正常状态。

北向路径对应轨迹4为清洁路径, 占2月气团轨迹的16.49%, 其输送的颗粒物质量浓度最低, 只有12.13 μg·m^-3(PM_2.5)和31.48 μg·m^-3(PM₁₀), 对应自然源路径通道, 主要为西伯利亚冷空气团高压南下, 轨迹移动路径长, 但速度最快, 故该路径气团移动过程中与其路径下垫面接触时间较短, 导致下垫面混合交换作用较弱, 气团不易携带沿途颗粒物南下, 且其上游贝加尔湖是清洁源, 轨迹途经地区排放源少。

对奥体中心而言, 重污染状态主要由轨迹2和轨迹3传输路径所致。轨迹2输送的PM_2.5和PM₁₀平均浓度可分别达到129.43和145.89 μg·m^-3, 对应空气质量轻度或中重度污染。该条路径经过唐山、天津在内的京津冀中东部地区, 工业活动较强, 排放强度较大；与此同时东北东部老工业区的排放颗粒物先也向海上传输后通过大气环流回流到华北平原造成奥体中心颗粒物浓度偏高。

轨迹3带来的PM_2.5和PM₁₀的平均浓度分别为223.8和263.24 μg·m^-3, 对应空气质量状况属于中重度污染。该路径移动速度最慢, 长度最短, 且移动轨迹主要在工业活动频繁的河北省中南部地区, 气团移动慢、有足够长的时间与污染物进行混合, 从而携带大量颗粒物输送至奥体中心, 严重影响奥体中心2月空气质量。

通过以上分析可以看出, 如若确保北京奥体中心2月空气质量达到国家二级标准, 需要对偏东轨迹与偏南轨迹所经过区域潜在污染源进行分析, 切断源头, 确保奥体中心2月空气质量。

3.2 潜在源分析

本文设计计算污染物后向72 h轨迹, 传输距离较远, 故研究的网格设置区域为北半球, PM_2.5重点关注区域范围为30°~50°N、105°~125°E。四种轨迹统计方法对奥体中心2月PM_2.5潜在来源分区结果如图 2所示。

从PSCF方法分析结果可以看出, 对奥体中心高浓度PM_2.5的污染物主要来自河北东部、天津和山东西北部, 在京津冀区域以外, 网格贡献率均小于0.1, 该方法无法确定不同潜在源网格对受体点污染浓度贡献, 这是PSCF方法的缺点之一。CWT与PSCF方法确定的高贡献潜在源区分布基本相似, 这主要是因为这两种方法均基于潜在源区网格点中的轨迹节点数量(或停留时间)来确定, 故一般距离受体测点越近的地区, 潜在源网格点是重污染来源的概率越高。不过由于CWT方法在计算过程使用了受体点的污染物浓度, 故可确定潜在源区对受体测点的污染贡献浓度, 可以看出山西东南部对奥体中心PM_2.5浓度贡献较小, 在50~75 μg·m^-3, 河北东部地区对奥体中心PM_2.5浓度贡献均超过了100 μg·m^-3, 海上某一区域可以达到150 μg·m^-3。整体而言, CWT方法可以在PSCF基础上对其低概率潜在源区进行区分, 但总体上这两种方法在污染物潜在来源追踪上针对性和指示性意义均不强。

针对CWT方法不足, RTWC方法对轨迹进行划分并迭代计算, 得到奥体中心2月大于200 μg·m^-3的潜在源地区主要为陕西东部和黄渤海面上颗粒物回流, 河北中部, 天津、山东半岛中南部等, 这一结果说明除了陆地上人为源传输通道外, 海上颗粒物回流及海盐也对PM_2.5也有较大贡献, 这与王跃思(2015)对北京污染物来源成分分析结果一致。此外, RTWC方法的水平分布可以直接、明确地指出奥体中心2月高贡献潜在污染源区。RTWC方法的计算结果与迭代结束判断标准有关, 本文设置的迭代结束标准是前后两次的差异小于0.5%, 基于以上迭代标准RTWC方法能较好识别出影响北京奥体中心的主要潜在污染源区。

与前面三种方法得到的结果不同, QTBA方法结果特征显著, 北京奥体中心纬度以南的区域对PM_2.5的潜在贡献均比较大, 虽然QTBA方法得到南部地区的贡献均在150 μg·m^-3以上, 但还是可以看出, 河北西南部、山西东南部、陕西北部、山东西部和黄海区域的QTBA值在200 μg·m^-3以上, 可以识别出众多排放源中的重污染源区, QTBA和RTWC方法在识别潜在源区的等级方面具有一致性。QTBA方法会导致潜在污染源范围估计偏大和虚假源的存在, 这主要是由于假设轨迹不确定性呈正态分布, 会产生一个尾端效应。

PM₁₀潜在来源分布结果如图 3所示。PSCF概率分布表明PM₁₀来源向西北方向扩展, 途径新疆、甘肃、宁夏、内蒙古中西部、陕西和山西中部等沙源地区, 上述区域的PSCF值小于0.1, 重点区域为京津冀周边和山东北部, 这些地区的PSCF值可以达到0.4以上。从CWT方法可以看出, 新疆北部和蒙古国浓度在75~100 μg·m^-3, 为奥体中心2月PM₁₀长距离输送源区；京津冀周边为100~150 μg·m^-3, 为2月奥体中心近距离源区；其余地区浓度小于50 μg·m^-3。CWT方法能识别出北京奥体中心PM₁₀的两个主要来源区域, 但对于长距离输送过程中气团裹挟沿途下垫面沙尘这部分污染物来源信息没能较好识别。RTWC方法通过迭代计算出PM₁₀质量浓度大于150 μg·m^-3以上的长距离输送潜在源贡献区为新疆北部及我国与蒙古国边境地区、河西走廊与内蒙古交界处、宁夏南部和陕西中部。短距离输送源区较少, 主要分布在河北中部、山东半岛和黄海部分地区。QTBA方法与前面三种方法的潜在源区的空间分布不同, QTBA方法计算的分布主要分为两个部分:西北路径部分和南向路径部分。其中西北部路径对应的大范围潜在来源区主要位于西北部沙源地和矿区等。南向轨迹路径上潜在污染源主要位于山西东南部、河北中南、山东半岛和黄海地区等地区。

基于以上分析可见, 四种轨迹统计方法都有助于确定污染物来源地区, 但其各有特点。PSCF方法计算简单, 不需要事先确定模拟区域的污染信息, 单纯由大量轨迹通过各点停留时间寻找潜在来源, 但为了保证其计算效率需要设置一个单元网格内轨迹节点阈值并需要确定受体测点的“污染”轨迹质量浓度最低值, 故不能计算所有轨迹；CWT方法则没有设置阈值, 可以计算模拟区域的所有气团轨迹, 但在长距离输送中, 距离受体点越远的位置, 网格点内的轨迹节点越少, 所以CWT方法在确定远距离高贡献源区时显得敏感性不足；针对CWT方法的这个问题RTWC方法可以通过迭代对高浓度网格进行加权来避免, 但美中不足的是迭代计算需要以牺牲计算效率为代价；QBTA方法比较适用于确定一些长期存在的、比较明确的污染来源, 但不确定性假设的引入可能会导致一些虚假源存在(如图 2中的大片红色区域), 不过通过与RTWC方法结合, 可以识别出其中的虚假来源。

4 结论与讨论

本文采用四种轨迹统计方法, 结合污染观测数据和气象观测数据, 通过轨迹模式分析了2013—2017年奥体中心2月污染输送特征和潜在来源分析, 主要结论如下。

(1) 通过聚类和统计分析得到:2月奥体中心高污染浓度轨迹出现概率为27.66%, 为偏东路径和偏南路径；北向路径为清洁路径, 概率为16.49%；西北路径出现概率最高(55.85%), 可认为是2月主导路径, 其对应颗粒物浓度可达到国家二级标准。

(2) 四种轨迹统计方法可确定河北地区、山东半岛和黄渤海区域为2月奥体中心PM_2.5主要来源地；新疆北部、河西走廊、河北地区和黄渤海区域为PM₁₀的来源地区, PM₁₀的远距离输送作用比较明显。2月受东南气流的海上颗粒物回流对奥体中心颗粒物浓度具有较高影响。

(3) 四种轨迹统计方法可分为两类:一类是由PSCF、CWT和QTBA方法组成, 这类的显著特点是计算速度快；另一类为RTWC方法, 其受迭代结束标准影响, 计算速度慢。PSCF和CWT方法在近距离污染物源区识别中具有较大优势, 但区分效果较差。RTWC方法虽然慢, 但可以通过迭代将污染物来源确定在具体网格范围内, 针对性较好, 在污染治理中具有较强的指示意义。QTBA方法适用于大范围污染物来源。RTWC与QTBA方法结合分析, 有助于减少QTBA方法虚假污染物来源。

本文仅以奥体中心单点为研究区域, 采用轨迹统计方法均在二维平面, 与实际的三维气团运动存在一定差异。故今后可考虑多站点、三维轨迹统计模型进行研究, 为今后北京大气污染来源解析和减排提供依据。