随着当今社会城市化进程的不断推进, 人类赖以生存的大气环境不断恶化, 大气污染问题已成为影响城市环境的重要因素^[1~3].大气污染具有区域性, 污染物浓度不仅受本地排放源的影响, 还在一定程度上受区域传输的影响^[4], 大气污染物以边界层为载体, 在大气多尺度时空系统作用下, 通过化学转化及光化学作用, 形成污染时空多尺度分布^[5].经济发达区域的多个城市可构成“城市群”, 其污染物扩散混合构成城市之间的羽流影响效应及其区域性大范围污染扩散现象^[6].研究显示, 城市间存在着大气污染的相互影响和输送污染物的传输路径^[4], 苏福庆等^{[7, 8]}利用统计方法研究大气污染物的输送特征, 分析了静稳天气下北京市外来污染物各类尺度的输入通道流场及典型边界层输送汇以及跨区输送对北京市空气质量影响, 陈朝晖等^[9]分析研究了华北区域大气污染过程中天气型及污染物输送路径.

目前运用后向轨迹模型HYSPLIT4(hybrid single-particle lagrangian integrated trajectory)^[10~12]通过模拟气流来源来判断上风向源区对目标地区是否存在影响是一种非常有效的方式.运用该模型, 国内外学者开展了大量研究工作, Anil等^[13]对伊斯坦布尔的可吸入颗粒物的远距离大气迁移影响进行了研究; Park等^[14]对亚洲大陆到东北亚和西北太平洋的SO₂远距离输送进行了研究; Fishman等^[15]对美国东部的O₃污染事件进行了研究; 石春娥等^[16]结合聚类分析的方法研究了合肥地区不同输送轨迹与该地区PM₁₀浓度之间的关系; 赵恒等^[17]研究了TRACE-P期间中国香港大气污染物的来源; 王郭臣等^[18]分析了地面和高空反向气流轨迹对北京的影响; 占明锦等^[19]分析了不同气团来源对于瓦里关地区颗粒物的数浓度及其谱分布的影响; 高晋徽等^[20]分析了长/近距离输送对南京地区O₃、NO₂和SO₂浓度变化的影响; 陈霖等^[21]研究了半干旱地区黑碳气溶胶和含碳气体来源; 杨阳等^[22]对北京沙尘前后气溶胶污染物的来源进行了分析; 此外, 相关学者利用后向轨迹并结合聚类分析、潜在源贡献因子分析法以及浓度权重轨迹分析法分析了西安^[23]、兰州^[24]、重庆^[25]、南昌^[26]、上海^[27]和郑州^[28]等地区的大气污染输送来源以及污染贡献源区.

上述研究多关注经济发达城市以及京津冀、长三角和珠三角地区, 且多为针对某次重大污染事件或者季节性某个要素的浓度值进行后向轨迹的分析, 对于较长时间尺度的后向轨迹综合聚类分析和污染源区的探究较少.太原位于汾渭平原北端, 是京津冀大气污染传输通道城市之一, 本研究选取太原2014~2017年污染较严重的秋冬季节进行后向轨迹的聚类分析, 并应用潜在源区贡献因子(PSCF)和浓度权重轨迹(CWT)进行综合分析, 研究太原市秋冬季大气污染物的传输来向和潜在源区, 从而有效地进行工业企业的合理布局, 控制入侵气团的污染浓度, 减少外来污染物对太原市城市空气质量的影响.

用于后向轨迹模式计算的气象场资料(2014年1月至2017年12月)为美国国家环境预报中心(NCEP)提供的全球资料同化系统(global data assimilation system, GDAS)数据, 每日4个时次, 即00:00、06:00、12:00和18:00(UTC), 水平分辨率为1°×1°.

用于监测资料分析的数据来源于原中华人民共和国环境保护部数据中心实时发布的AQI(air quality index, 空气质量指数；2014年1月至2018年12月每日AQI以及2014年1月至2017年12月每小时AQI).

本研究采用美国NOAA空气资源实验室研制的HYSPLIT4模式来计算和分析大气污染物的来源、输送、扩散轨迹, 识别大气污染物随气团的移动方向.该模式模拟精度高, 时间连续变化, 被广泛应用于大气输送研究及污染过程分析中.由于模拟高度和轨迹运行时间的选取对HYSPLIT模式的模拟精度有较大影响, 本研究的后向轨迹计算起始高度选取为300 m, 能反映气流的区域性流动特征, 又可减小下垫面摩擦力的影响, 能够更好地代表边界层底层的情况^[5].轨迹运行时间选择为72 h, 以便包含二次污染物的生命周期^[29], 时间间隔取为1 h.

聚类分析是一种多元统计技术, 可对大量的数据进行分组, 合理有效地简化数据, 提高信息的利用率, 广泛地应用于大气污染研究中.根据气团移动速度和方向对大量轨迹进行分组, 得出不同的输送轨迹组, 可以比较直观地了解气团的来源方向和传输距离, 从而估计大气污染物的潜在源区.本研究采用二分k均值算法进行聚类分析, 该算法首先将所有的轨迹作为一组, 然后将该组一分为二, 之后选择其中的一个组继续划分, 选择哪一个组进行划分取决于对其划分是否可以最大程度地降低误差平方和.重复上述基于距离平方和的划分过程, 直到组数目和指定的要聚类的数量相等.

每个轨迹组的空间方差(total spatial variance, TSV)为组内每条轨迹与组平均轨迹对应点的距离平方和.轨迹间距离采用欧式距离进行计算, 计算公式如下：

(1)

式中, X₁(i)和Y₁(i)表示轨迹一的坐标, X₂(i)和Y₂(i)表示轨迹二的坐标.

潜在源贡献因子分析法(potential source contribution function, PSCF)是基于条件概率函数发展而来的一种方法, 它基于气流轨迹分析来识别可能源区^[30], 通过结合气团轨迹和某特征值(如AQI、污染物浓度值等)来给出可能的污染排放源方位.

PSCF函数定义为经过研究区域的气团到达观测点时对应的某要素值超过设定的阈值的条件概率.将研究区域划分为一定分辨率的若干个网格, 每个网格记为网格(i, j), 对研究的要素设定一个阈值, 当轨迹所对应的要素值高于这个阈值时, 认为该轨迹是污染轨迹, 经过网格(i, j)的污染轨迹端点数为m_ij, 而落在网格(i, j)内的所有轨迹端点数为n_ij, 则PSCF可以定义为公式：

(2)

对于某些偏远网格, 由于轨迹总体停留时间较少及分母n_ij较小, 可能会引起PSCF计算结果的较大不确定性.为了降低这些特殊网格对于计算结果的影响, 引入权重系数W_ij, 当某一网格内的n_ij小于研究区域内每个网格的平均轨迹端点数的一定位数时, 就要使用W_ij来减小PSCF的不确定性.

(3)

(4)

浓度权重轨迹分析法(concentration weighted trajectory method, CWT)^{[31, 32]}是一种计算潜在源区气流轨迹权重浓度, 反映不同轨迹的污染程度的方法, 可以定量给出每个网格的平均权重浓度, 计算公式为：

(5)

(6)

式中, C_ij是网格(i, j)上的平均权重浓度; k是轨迹; M是轨迹总数; C_k是轨迹k经过网格(i, j)时对应的某要素值(如AQI、污染物质量浓度等); τ_ijk是轨迹k在网格(i, j)停留的时间.在PSCF分析法中所用的权重函数W_ij也适用于CWT分析法.

依据历年中国生态环境状况公报, 新标准第一阶段监测实施城市(共74个城市)中按照环境空气质量综合指数评价, 太原市的环境空气质量指数2015~2017年分别位于倒数的第6、10和13位, 环境空气污染较为严重^[33]. 图 1给出了太原市2014~2018年逐日AQI变化, 年初1月前后出现最高值, 在8月前后出现最低值, 2014~2018年逐日AQI达标天数分别为202、233、236、176和182 d, 近两年达标天数显著下降, 表明太原地区大气环境质量改善的基础仍十分脆弱, 环境保护面临的形势依然严峻.

图 1

Fig. 1

图 1 太原市2014~2018年环境空气质量指数逐日分布 Fig. 1 Daily distribution of air quality index in Taiyuan from 2014 to 2018

基于逐日AQI数据, 分析了2014~2018年AQI月变化及不同AQI等级占比情况(图 2).从图 2(a)可以看出, 太原市AQI月均值呈现出“凹”型分布, 其中8和9月为全年AQI最低的月份, 同时也是优良天数占比最高的月份[图 2(b)], 9月空气质量为优的天数为全年最多; 图 2中也显示, 11月虽不是AQI最高的月份, 但优良天数最少; 太原以及周边地区从10月中下旬进入采暖期, 因此, 10月至翌年2月AQI一直保持在较高的状态, 同时也是中度污染以上天数占比较高的月份, 本研究将重点关注此时段污染物传输特征, 以下秋冬季均指这5个月.

图 2

Fig. 2

图 2 2014~2018年AQI逐月变化以及不同AQI等级逐月占比情况 Fig. 2 Monthly variation of AQI and different AQI grades in Taiyuan from 2014 to 2018

图 3为不同等级AQI全年及秋冬季占比变化情况.由图 3(a)可以看出, 2015~2016年优良天数比例较高超过60%, 2017年又下降到不足50%;伴随优良天数的减少, 中度、重度污染以上出现天数比例明显上升.从秋冬季来看[图 3(b)], 2014~2018年秋冬季优良天数占比逐步上升, 2018年超过50%, 空气质量好转趋势明显, 可能得益于政府在秋冬季实行的一系列污染管控措施, 2017和2018年秋冬季节多次发布重污染天气预警, 重污染天气期间对重点污染企业进行关停和限产措施, 并多次实施了车辆限行.结合全年及秋冬季优良天数年变化特征分析表明, 近年来全年优良天数下降可能来自于春季沙尘天气影响和夏季O₃污染.

图 3

Fig. 3

图 3 2014~2018年全年及秋冬季各AQI等级逐年占比情况 Fig. 3 Annual, autumn, and winter variation in different AQI grades in Taiyuan from 2014 to 2018

图 4给出了不同年份及逐年秋冬季首要污染物占比情况.从全年来看, 颗粒性污染物为首要污染物的天数占主导地位, 但呈逐年下降趋势.其中, PM_2.5为首要污染物的天数占比下降趋势非常显著, PM₁₀占比则有上升趋势; 此外, O₃为首要污染物的天数占比逐年上升趋势明显, SO₂则有明显地下降趋势.从秋冬季来看, 总体变化趋势同全年类似, PM_2.5为首要污染物的天数占比下降趋势明显, 而PM₁₀则有明显的上升趋势; 不同的是秋冬季PM_2.5为首要污染物的天数比例明显要高, 接近全年均值的2倍, 而PM₁₀为首要污染物的天数占比在近两年增幅明显要高于全年, 表明太原地区近年来污染类型可能发生变化, 因此在污染治理过程中不仅要继续加大秋冬季颗粒物治理力度, 也要进一步关注其他季节不同类型的污染情况.

图 4

Fig. 4

图 4 2014~2018年全年及秋冬季各首要污染物占比情况 Fig. 4 Annual, autumn, and winter variation in the proportion of primary pollutants in Taiyuan from 2014 to 2018

采用后向轨迹模式计算了到达太原(37.73°N, 112.55°E)的2014年1月1日至2017年12月31日这4a逐时后向轨迹.采用二分k值法对其进行聚类, 得到太原市秋冬季典型的8类输送轨迹, 水平方向轨迹分布见图 5.轨迹的路线和方向表示气团在到达计算点(太原市)所经过的地区, 根据其长短可以判断出气团移动的速度, 长轨迹对应移动快速的气团, 短轨迹对应移动缓慢的气团, 垂直方向的变化表示气团在垂直方向混合和垂直运动的信息. 表 1给出了8类轨迹的平均速度, 其中第1、3和7类轨迹平均移动速度均超过10m·s^-1, 为远距离输送轨迹; 第2和5类轨迹平均移动速度在8 m·s^-1左右, 为中距离输送; 第4、6和8类轨迹移动速度低于5m·s^-1, 为近距离输送轨迹.整体来看(图 5), 3类西北轨迹(第1、2和3类)占21.1%, 3类偏西轨迹占53%(第5、6和7类), 西南轨迹占12%(第8类), 偏东方向轨迹占13.9%(第4类).轨迹结果与太原所处的地理位置和季风气候是相适应的, 太原市地处中纬度地带大陆内部, 其西北部为广阔的欧亚大陆腹地, 西北风盛行, 距东部海岸线较远, 且有太行山脉屏障, 偏东气流对本区域影响较小.

图 5

Fig. 5

图 5 太原市2014~2017年秋冬季8类聚类轨迹 Fig. 5 Eight trajectory clusters to Taiyuan in autumn and winter from 2014 to 2017

表 1 (Table 1)

表 1 太原市秋冬季8种聚类轨迹气团平均移动速度 Table 1 Air mass movement velocity of eight trajectory clusters to Taiyuan in autumn and winter 轨迹分类 1 2 3 4 5 6 7 8 平均风速/m·s-1 17.6 7.1 11.6 1.6 8.9 4.7 14.1 2.6

表 1 太原市秋冬季8种聚类轨迹气团平均移动速度 Table 1 Air mass movement velocity of eight trajectory clusters to Taiyuan in autumn and winter

图 6给出了太原市秋冬季8种聚类轨迹垂直分布.结合轨迹水平和垂直分布情况可以看出, 第1、3和7类长距离轨迹在垂直方向上也有较大跨度, 传输过程有较为明显的先上升后下沉的过程, 说明这3类气团先是从源地的抬升, 然后以较快的速度向下游输送, 并在垂直混合和下沉运动的作用下达到目标城市; 第5、6和7类轨迹为方向一致但水平长度显著不同的轨迹, 其中第6类轨迹的气团速度和整体高度(垂直速度)均显著低于第5和7类轨迹, 轨迹越短、垂直高度越低, 更易将沿途的污染裹挟向下游输送, 因此第6类轨迹对目标地区的影响程度可能大于其他两类; 中距离轨迹第2和5类在垂直高度较为接近, 近距离轨迹第4、6和8类在垂直高度也较为接近, 但水平方向上经过区域不同, 对目标区域的影响也可能产生差异. 图 6中还可以看出, 第8类轨迹在垂直方向有明显的先下沉后上升的现象, 可能更易于传输气流中污染物的积累.

图 6

Fig. 6

图 6 太原市秋冬季8种聚类轨迹垂直分布 Fig. 6 Profile of eight trajectory clusters to Taiyuan in autumn and winter

图 7给出了太原市秋冬季8种聚类轨迹下AQI分布箱式图, 表 2给出了不同聚类轨迹下优良天气统计特征.结果显示, 不同轨迹下对应的太原市AQI均值有显著差异.第1、2和3类轨迹下太原市AQI的优良比例均在70%以上, 其中偏北方向的第2类轨迹下太原AQI优良比例最高, 为最清洁轨迹.偏西方向的第5、6和7类轨迹下AQI分布特征较为相似, 污染天气比例均在50%以上, 这3条轨迹均可能造成太原市AQI偏高, 其中第6类轨迹下太原污染天气比例均略高于第5和7类, 说明气团高度较低速度较慢时, 更容易将气团移动路径上的污染传输至下游.第4类轨迹是唯一的偏东方向轨迹, 4类轨迹下太原污染天气比例同第6类轨迹接近, 但优良比例低于第6类轨迹, 因此该轨迹对目标城市的影响较第6类轨迹影响偏小.来自西南方向的第8类轨迹下太原市污染比例超过84%, 该轨迹下太原市极易发生污染, 且污染程度也最重, 表明该轨迹为典型的污染轨迹, 其在垂直运动过程中的下沉后再抬升现象更容易将沿途路径的低层污染裹挟向下游输送, 造成目标城市的严重污染.

图 7

Fig. 7

图 7 太原市秋冬季8种聚类轨迹AQI分布情况 Fig. 7 AQI of eight trajectory clusters to Taiyuan in autumn and winter

表 2 (Table 2)

表 2 不同聚类轨迹统计分析结果 Table 2 Statistical analysis of eight trajectory clusters to Taiyuan in autumn and winter 轨迹类型 优良比例/% 污染比例/% 1 71.32 28.67 2 82.69 17.31 3 81.16 18.84 4 46.52 53.48 5 47.53 52.47 6 46.70 53.31 7 48.41 51.59 8 15.95 84.05

表 2 不同聚类轨迹统计分析结果 Table 2 Statistical analysis of eight trajectory clusters to Taiyuan in autumn and winter

表 3为太原市在不同AQI等级下各聚类轨迹分布情况.发生轻度污染时, 第6和5类轨迹出现几率较高均超过20%, 第4类轨迹出现的几率也接近20%;发生中度污染时, 以第6类轨迹出现的几率最高达到30%, 第4、5和8类轨迹出现的几率约20%, 较为接近; 当发生重度污染时, 第6和8类轨迹出现的比重相当均为28%, 其次是第5类轨迹比重为22%;发生严重污染时, 第8类轨迹出现的比重达到34%, 其次是第6类轨迹28%和第5类轨迹19%.

表 3 (Table 3)

表 3 太原市不同AQI等级发生时各聚类轨迹所占的比重/% Table 3 Proportion of trajectory clusters under different AQI grades/% AQI等级 轨迹分类 合计 1 2 3 4 5 6 7 8 优 3.01 19.49 14.45 18.64 14.58 21.78 4.32 3.73 100 良 4.65 14.64 13.59 11.59 19.86 25.70 6.12 3.85 100 轻度污染 2.46 5.05 4.62 17.89 21.78 28.40 6.55 13.24 100 中度污染 1.31 2.61 2.75 17.52 19.15 30.26 6.21 20.20 100 重度污染 1.33 1.67 2.25 9.92 21.74 28.32 5.94 28.84 100 严重污染 0.37 1.12 0.19 12.34 18.69 28.04 4.86 34.39 100

表 3 太原市不同AQI等级发生时各聚类轨迹所占的比重/% Table 3 Proportion of trajectory clusters under different AQI grades/%

综合以上分析, 第8类轨迹是外来污染物输送入太原的主要轨迹, 其次为第6类轨迹和第5类轨迹, 第4类轨迹主要造成目标城市轻度和中度污染.

对太原市气流轨迹的聚类分析能够确定出影响太原市的气流轨迹的大致方向和贡献, 但难以明确地反映潜在源区的位置和贡献大小, 而潜在源贡献因子分析法(PSCF)和浓度权重轨迹分析法(CWT)可以分析潜在污染来源和贡献大小.潜在源区一般结合某种污染物浓度来分析, 例如对某地区PM_2.5^[24]或PM₁₀^[31]的潜在源区进行分析, AQI作为一个综合指标可以更宏观地表征一个地区的总体污染程度, 本研究尝试结合太原地区AQI指数来对2014~2017年污染物潜在源区进行分析, 从而获得引起太原市AQI增高的潜在源区, 以期对引起太原市秋冬季污染的潜在源区有更深入地了解.

PSCF基于条件概率函数来判断污染源可能方位, 因此从一次污染过程或者短期污染过程来进行潜在源区贡献的计算并不能准确反映影响某地的潜在源区贡献, 或者说只能局限于某次过程的潜在源区贡献, 因此本研究基于2014~2017年秋冬季的后向轨迹来进行PSCF计算, 以更全面反映潜在源区对太原市的长期影响和贡献.太原市秋冬季污染的潜在源区贡献计算结果见图 8, 颜色越深, 表示PSCF的计算值越大, 表明该网格对计算点太原市AQI的影响越大.由图可见, 对太原市AQI影响较大即PSCF高值区(潜在贡献0.8以上)主要位于甘肃庆阳以南、陕西汉中以北包括西安、铜川、陕西延安以南和山西的吕梁、临汾以及河南洛阳以北、山西晋城以西等地, 这些潜在源区的气团主要沿着第8类轨迹通过近距离输送到达太原; 通过第6和5类轨迹进行中远距离输送的张掖以北、银川以北等区域, 潜在贡献0.6左右; 通过第4类轨迹进行近距离输送的河北石家庄、邢台等区域潜在贡献为0.5左右.

图 8

Fig. 8

图 8 太原市潜在源区贡献分析结果 Fig. 8 PSCF analysis in autumn and winter from 2014 to 2017 in Taiyuan

运用CWT分析可以进一步确定潜在源区的权重浓度(本文基于AQI计算, 即AQI的大小), 见图 9.结合图 8和图 9可以看出, 两种方法结论基本一致. CWT计算结果表明, 对太原市AQI值贡献超过200以上的区域主要是甘肃庆阳以南、陕西汉中以北包括西安和铜川等地; 对太原市AQI值贡献超过150以上的地区是内蒙阿拉善盟以南、陕西延安以南、山西的吕梁以南、临汾以及河南洛阳、郑州以北、山西晋城等地; 对太原市AQI值超过100以上贡献的区域主要是内蒙中西部、陕西榆林、山西晋中、长治以及河北西南、河南郑州、开封以北地区.以上这些区域对太原市秋冬季AQI的贡献在100以上, 是重要的潜在源区.

图 9

Fig. 9

图 9 太原市浓度权重轨迹分析结果 Fig. 9 CWT analysis in autumn and winter from 2014 to 2017 in Taiyuan

通过对太原2014~2017年输送路径和潜在源区分析表明, 偏西南方向的第8类轨迹是造成太原秋冬季大气严重污染的主要输送途径, 太原市大气严重污染主要的潜在源区就是汾渭平原的核心区域, 包括陕西省西安、宝鸡、汉中等区域以及山西省汾河谷地的吕梁、晋中、临汾和运城等地区.杨晓春等^[23]和王珊等^[34]对西安地区一次重污染过程的后向轨迹进行了模拟, 段时光^[28]等对郑州冬季PM_2.5的输送路径进行了分析, 结果显示, 山西汾河谷地也是西安、郑州等地区重污染过程中输入污染物的重要来源之一, 因此, 太原和周边地区的影响是显著且相互的, 建立起有效的区域联防联控机制, 才能提高空气污染治理效果.

(1) 太原市2014~2018年优良天数波动较大, 2014年优良天占55%, 2015、2016年上升到64%, 2017、2018年下降到不足50%;秋冬季变化则不同, 2014~2018年秋冬季优良天数占比逐步上升, 2018年超过50%, 空气质量好转趋势明显.太原8、9月空气质量最佳, 从10月开始转差一直持续到翌年2月.污染类型发生变化, 全年PM_2.5为首要污染物的天数显著下降, PM₁₀和O₃为首要污染物的天数明显上升, 秋冬季亦是PM_2.5为首要污染物的天数下降, PM₁₀上升.

(2) 利用HYSPLIT4模型和气流轨迹聚类分析2014~2017年秋冬季太原的后向轨迹, 得到8类轨迹, 21%来自西北方向, 53%来自偏西方向, 12%来自西南方向, 14%来自偏东方向.远距离输送对太原影响较小, 近距离输送对太原影响较大, 西南方向轨迹下太原有84%的比例出现污染天气, 且多重污染, 是外来污染物输送进入太原的主要轨迹, 偏西方向近距离轨迹和偏东方向轨迹下太原发生污染天气的比例也高于50%.

(3) 潜在源贡献因子分析和浓度权重轨迹分析结果具有一致性, 甘肃庆阳以南, 陕西的汉中、西安、铜川、延安和山西的吕梁、临汾、晋城以西以及河南洛阳以北等地是对太原污染贡献最大的区域, 表明建立汾渭平原及周边区域联防联控机制对控制区域污染有着重要意义.