基于Kolmogorov-Zurbenko滤波法分析2015—2021年京津冀地区大气颗粒物变化趋势

占 青,张运江,陈 红,张可馨,盖鑫磊

(南京信息工程大学环境科学与工程学院/ 江苏省大气环境监测与污染控制高技术研究重点实验室/ 大气环境与装备技术协同创新中心,江苏 南京 210044)

大气颗粒物是影响人体健康和环境质量的主要污染物,按照颗粒物粒径进行划分,可分为细颗粒物(PM2.5)、可吸入颗粒物(PM10)和粗颗粒物(PM2.5-10)。我国城市规模扩大、工业化迅速扩展以及长期积累的能源结构矛盾,使得区域大气状况日益恶化[1]。为了统筹推进区域空气环境治理工作,2013年国家开始实施《大气污染防治行动计划》《打赢蓝天保卫战三年行动计划》[2]。北京、天津、河北三地在积极响应国家号召的同时也在积极探讨区域环境污染的协调机制,以降低环境中有机物排放量和改善三地空气环境质量;2018和2019年连续发布《京津冀及周边地区秋冬季大气污染综合治理攻坚行动方案》;2020年出台《天津市大气污染防治条例》[3]。京津冀作为国家经济发展中心,同时也是连接中原和华北平原的主要通道,通过分析京津冀大气环境质量和污染物空间分布特征,对实现区域空气质量联合控制起到至关重要的作用。

自2013年以来,我国各大城市空气环境监测中心陆续公布各大城市大气污染指数,并发表了一批关于京津冀、长三角、珠三角、成渝和汾渭平原等区域的研究成果[4-5]。基于大气化学传输模型和污染源排放清单数据,近年我国学者研究发现减排行动是主导2013—2017年京津冀地区大气颗粒物污染改善的重要因素,贡献达80%[6-9]。采用大气化学传输模式评估排放对大气颗粒物污染贡献的优点之一是可以定量掌握格点化的大气污染物排放总量与空气质量变化的直接响应关系。然而,基于大气化学传输模式的方法需要最新的排放清单数据作为模型输入项。因此,近年来很多学者选用基于观测数据的统计学方法来评估气象条件和排放对大气污染的贡献[10-12]。基于观测数据的统计学方法可以直接采用观测数据进行评估,而不依赖于排放数据,有助于快速掌握排放、气象条件对污染物变化趋势的贡献,目前该方法已被国内外学者广泛使用。例如,ZHAI等[13]基于逐步多元线性回归模型(MLR)量化气象条件对中国各地PM2.5的贡献趋势。在国外,该方法还在美国亚利桑那州图森地区[14]和韩国[15]得到应用,研究这些地区O3和颗粒物的污染趋势、滤波分量与气象参数之间的相关性,以及耦合其他方法开展O3污染预测研究[16]。但这些研究基本集中在2019年及之前,而近年来减排对大气颗粒物的影响尚不清楚,对PM10尤其是PM2.5-10变化趋势的了解也较少。虽然目前京津冀地区大气环境质量得到一定改善,但区域颗粒物污染状况依然不容乐观[17-18]。因此,以2015—2021年京津冀地区典型城市大气颗粒物地面观测数据和气象要素数据为基础,基于观测数据的统计学方法,即Kolmogorov-Zurbenko(KZ)滤波法,从时间、空间角度剖析区域大气颗粒物(包括PM2.5、PM10和PM2.5-10)污染特征及变化趋势,同时定量评估气象条件和排放对变化趋势的贡献,为制定京津冀地区大气颗粒物减排对策提供科学参考。

1 材料与方法

1.1 研究区域及数据来源

选取京津冀地区13个城市作为研究对象,即保定、北京、沧州、承德、邯郸、衡水、廊坊、秦皇岛、石家庄、唐山、天津、邢台和张家口。采用KZ滤波方法,选取2015年1月1日至2021年12月31日13个城市PM2.5、PM10和PM2.5-10数据进行分析。空气质量数据来自于全国大气环境质量监控网格,包括PM2.5、PM10、SO2、NO2、O3和CO数据。相应气象数据来自于欧洲中期天气预报中心(ECMWF)ERA5全球再分析数据,使用的气象要素包括:经纬风(U10和V10)、边界层高度(BLH)、相对湿度(RH)、总降水量(TP)、太阳辐射(SR)、距地面2 m温度(T2m)、气压(SP)、云覆盖(TCC)和850百帕经纬风和垂直输送变量(u850、v850和w850)。

1.2 研究方法

1.2.1KZ滤波法

KZ滤波方法计算过程简便,且无需对缺失值进行特殊处理,已广泛用于北京、天津、河北等地的大气污染指数、O3和颗粒物的相关分析[19-23]。KZ滤波器通常被应用于消除高频噪声、将低频成分与原始信号分开以及消除季节效应3个典型场合[24-26]。

在此基础上,空气质量数据的时间序列可以表示为

A(t)=e(t)+S(t)+W(t)。

(1)

式(1)中,A(t)为原始序列;e(t)为长期分量;S(t)为季节分量;W(t)为短期分量。长期趋势是由总排放、污染物传输、气候、政策或经济方面的改变引起,季节分量是由太阳角度的变化引起,短期分量可以解释为天气和污染物排放的短期波动。

KZ滤波器是一种通过移动平均的重复迭代产生的低通滤波器。将KZ(m,p)滤波器(表示窗口长度为m和迭代次数为p的滤波器)的移动平均定义为

(2)

式(2)中,Yi为经KZ滤波后的时间序列;k为Ai在滤波时其两端的滑动窗口长度;m为滑动窗口长度,m=2k+1;i为序列的时间间隔;j为滑动窗口变量;A为输入时间序列。迭代时,上一次传递的输出将成为下一次传递的输入。通过调整窗口长度和迭代次数,可以控制不同运动尺度的滤波。对于小于N天的周期,采用以下准则来确定滤波器有效宽度:

m×p1/2≤N。

(3)

因此,KZ(15,5)滤波器(表示窗口长度为15和迭代次数为5的滤波器)将去除<33 d(15×51/2≤33)[24]的周期。Eq 3方法可用来确定滤波器有效宽度,即滤波器权值约为0.5时的窗口大小,而不是截止周期,这是滤波器传递函数的一个性质[14]。在以往KZ滤波器研究中,有效宽度和截止时间都被使用[24,26]。采用KZ(15,5)过滤器过滤后的气象和空气质量时间序列(分别为MBL和ABL)被称为基线分量。将基线空气质量定义为长期和季节分量的总和,计算公式为

KZ(15,5)=e(t)+S(t)。

(4)

每年的周期和更小的时间尺度可以通过选择更大的窗口来移除。KZ(365,3)滤波器可以通过创建包含周期为632 d(365×31/2≤632)或1.7 a的周期的时间序列来从数据中提取长期分量,其结果是原始空气质量时间序列的长期成分:

e(t)=KZ(365,3)。

(5)

根据长期分量,可以计算出季节和短期分量。数据的季节分量是由空气质量和气象资料的基线序列减去长期成分而得到:

S(t)=KZ(15,5)-KZ(365,3)。

(6)

短期时间序列是由原始时间序列减去基线而得到:

W(t)=A(t)-KZ(15,5)。

(7)

1.2.2逐步多元回归模型

通过对浓度-时间序列的短期成分〔M(t)〕和基线成分〔XB(t)〕的多元线性回归模型进行研究[27],结果表明:

M(t)=a0+∑aiSi(t)+εM(t),

(8)

XB(t)=b0+∑bjBj(t)+εB(t)。

(9)

式(8)～(9)中,εM(t)和εB(t)分别为短期分量回归残差和基线分量回归残差;Si(t)为短期分量回归的第i个气象因子;Bj(t)为基线分量回归的第j个气象因子;a0、ai、b0和bj为回归系数。

污染物长期分量在均值附近的波动主要有2个方面原因:一是污染物排放的变动,二是天气状况的改变。在进行长期成分研究时,应将其与大气和污染物的环境因子区分开。所以,采用逐步多元线性回归模型是一个较好的选择。采用SPSS 26.0建立各个时期大气污染物与气象条件之间的对应关系模型。

浓度时间序列回归的总残差〔ε(t)〕表达式为

ε(t)=εM(t)+εB(t)={M(t)-[a0+∑aiSi(t)]}+{B(t)-[b0+∑bjBj(t)]}。

(10)

(11)

评价指标(K)为协方差总和与总方差之比,K值越小,分解效果就越好,其计算公式为

(12)

式(12)中,D(X)为ρ(PM2.5)、ρ(PM10)和ρ(PM2.5-10)时间序列总方差;D(W+S+e)为污染物短期分量、季节分量和长期分量的总方差;cov(e,S)、cov(e,W)和cov(S,W)分别为长期分量与季节分量、长期分量与短期分量、季节分量与短期分量的协方差。

气象条件贡献率计算公式为

(13)

式(13)中,Pcontrib为气象条件贡献率;Korigin和Kfilter分别为原始数据趋势斜率和滤波后趋势斜率。

2 结果与分析

2.1 总体污染概况

京津冀地区西侧为太行山脉,北侧为燕山山脉,东临渤海湾,地形呈现西北高、东南低的特征,因复杂地形产生的山谷风是造成该地区颗粒物污染的重要原因[28]。已有调查结果显示,京津冀地区主要受到北京西侧太行山东侧的西南向传输路径的影响,而西南部和东南部的传输路径为北京地区典型的污染物传输路径,PM2.5在西南部的传播通量最大,有时甚至超过60%[29]。同时,京津冀地区极易出现严重的沙尘暴[30]。京津冀及其周围区域位于华北大草原,受到西北风沙的影响较大,而风沙不但增加PM10浓度,还增加PM2.5浓度[31]。图1为京津冀地区2015—2021年PM2.5、PM10和PM2.5-10年际变化趋势及平均浓度分布,京津冀13个城市PM2.5浓度均呈下降趋势,其中PM2.5浓度高值区以河北省南部地区(邯郸、邢台、石家庄和保定)为主,该区域PM2.5浓度下降趋势也尤为明显,下降幅度介于32.3%～53.5%之间;河北省南部地区(张家口、承德和秦皇岛)PM2.5浓度最低,该区域PM2.5浓度下降幅度约为30%。13个城市PM10浓度下降幅度介于22.7%～54.7%之间,其分布特征与PM2.5浓度分布特征一致,PM10浓度高值区也以河北省南部地区(邯郸、邢台、石家庄和保定)为主,该区域PM10浓度下降幅度介于33.6%～50.9%之间。13个城市PM2.5-10浓度下降幅度介于9.8%～56.1%之间,其分布特征与PM2.5和PM10浓度分布特征一致,高值区以河北省南部地区(邯郸、邢台和石家庄)为主,其下降趋势介于17.1%～43.4%之间。

图1 2015—2021年京津冀地区PM2.5、PM10和PM2.5-10年际变化趋势及平均浓度分布特征

图2为京津冀地区PM2.5、PM10年均浓度和不同程度污染时段(优、良、轻度污染、中度污染和重度污染)占比。如图2所示,京津冀地区PM2.5年均浓度从2015年的75.32 μg·m-3迅速下降至2021年的39.82 μg·m-3,年均下降约5.07 μg·m-3;PM10年均浓度从2015年的136.52 μg·m-3迅速下降至2021年的79.82 μg·m-3,年均下降约8.10 μg·m-3。

图2 2015—2021年京津冀地区PM2.5和PM10年平均浓度和各污染时段占比

根据GB 3095—2012《环境空气质量标准》和HJ 633—2012《环境空气质量指数(AQI)技术规定(试行)》评价标准,进一步分析颗粒物浓度发现,2015年PM2.5和PM10污染(即24 h PM2.5平均浓度>75 μg·m-3、24 h PM10平均浓度>150 μg·m-3)天数分别为99(占全年有效数据天数的43%)和70 d(占全年的31%),其中,PM2.5和PM10重度污染(即24 h PM2.5平均浓度>150 μg·m-3、24 h PM10平均浓度>350 μg·m-3)天数分别为15(占全年的7%)和1 d(占全年的不到1%)。2021年,PM2.5和PM10污染天数分别下降为43(占全年的12%)和30 d(占全年的8%),其中,PM2.5和PM10重度污染天数分别为1和2 d。有关研究结果表明,京津冀城市群颗粒物浓度分布不是完全随机的,而是表现出空间聚集性,城市群相似的高能耗工业结构布局是影响区域颗粒物污染的重要因素,因此,调整产业结构、统筹工业布局是降低京津冀区域颗粒物污染的主要措施之一[32]。

2.2 典型城市PM2.5、PM10和PM2.5-10浓度时间序列分解

图3为采用KZ滤波器得到的PM2.5、PM10和PM2.5-10时间序列各分量。在KZ滤波过程中,由于移动平均的重复迭代,时间序列开始和结束部分数据被截断。更长的窗口长度会导致额外数据丢失。因此,长期趋势分量比短期时间序列有更多数据被截断。需要注意的是,如果不应用时间序列分离方法,长期分量的相对较小的变化范围将很难从原始时间序列中分离出来[25]。

如图3所示,PM2.5、PM10和PM2.5-10长期分量均略有波动,但总体呈降低趋势,其中PM2.5、PM10和PM2.5-10浓度分别降低36.30、58.57和22.77 μg·m-3。短期分量变化幅度明显,尤其冬天变化幅度最大。在短时间内,排放源变化不会太大,所以短期分量的剧烈变动可看作是气候的短期变动。季节分量表现出明显周期性变化,主要表现为6—12月和12月—次年6月的波动。冬季大气颗粒物浓度高,可能是受到冬季取暖的影响[33]。冬季供暖能源大多为一次能源,燃烧会产生大量颗粒物污染。除受供暖影响外,还受到冬季气象条件影响[34],例如,降雪对颗粒物具有高达79%～100%的去除率[35]。另外,2019年11月15日至2020年3月15日为停工期(暂停项目包括水利工程、土石工程、各种道路施工和房屋拆除),同时2020年初新冠疫情爆发,综合各种情况,京津冀地区颗粒物浓度在冬季虽有所增加,但在一定程度上得到有效控制,未出现长期严重污染现象。在污染源一定的条件下,由于各季节气象因子不同,导致各地区颗粒物浓度存在明显季节性变化[36],如,秋季较弱的冷空气活动和低地面风速导致大气边界层降低,从而使得大气低层水汽积累,形成重污染天气[37]。在春季,当风速大于2级时,强风使空气中颗粒物浓度急剧上升[38],同时由于重污染天气,颗粒物浓度也随之升高。夏季降水对颗粒物的去除效果显著[39],但夏季高温高湿天气也会使空气中污染物扩散条件出现不同程度恶化,从而使空气中颗粒物浓度增加,因此总体上夏季空气中颗粒物浓度变化不大[40]。

通过分析短期分量、季节分量和长期分量对原始序列的影响,可以得到各分量对原始序列的贡献率。在理想条件下,3个成分是彼此独立的,所以原始序列变化应等于3个成分的变化之和。图4为京津冀地区PM2.5、PM10和PM2.5-10各分量方差对原始序列方差的贡献,其中未知项为原始序列方差与3个分量方差和的差值。如图4所示,京津冀地区各城市PM2.5、PM10和PM2.5-10的短期、长期和季节3个分量的总和均超过90%,且K值在0.04～0.10之间(表1),说明3个分量基本上都符合独立性,证明KZ滤波法适用于京津冀城市群13个城市PM2.5、PM10和PM2.5-10浓度的序列分解。3个分量中,短期分量贡献最大,均在50%以上,季节分量次之,长期分量最少。由此可以看出,原始时间序列的波动主要与短期和季节影响有关,即与污染源排放和气象条件的短期和季节变化有关。因此,为了更好地研究长期变化趋势,必须将长期分量与原始序列分离。

表1 2015—2021年京津冀13个城市污染物评价指标(K)

图4 京津冀13个城市PM2.5、PM10和PM2.5-10分量贡献

2.3 减排措施与气象条件对颗粒物长期分量变化趋势的影响

图5为2015—2021年京津冀13个城市经气象调整后PM2.5、PM10和PM2.5-10浓度长期分量序列以及气象条件对PM2.5、PM10和PM2.5-10浓度长期分量的影响。如图5所示,京津冀13个城市PM2.5、PM10和PM2.5-10长期分量变化趋势在调整后与调整前基本一致,只在细节上有一些差异,长期分量与经气象调整后长期分量都呈现波动下降趋势,说明长期分量变化趋势主要与污染物排放变化有关。为了得到气象条件对长期分量的影响,且能更直观表现气象条件对PM2.5浓度的影响,将长期分量与经气象调整后长期分量做差值。当差值为正时,表明污染源导致的长期分量小于污染源和气象因素影响的长期分量,气象条件对污染状况改善是有利的;当差值为负时,表明由污染源导致的长期分量大于污染源和气象因素影响的长期分量,表明气象条件对环境改善不利[20]。如图5所示,京津冀13个城市气象条件对PM2.5、PM10和PM2.5-10长期分量的影响均在“0”上下波动,这表明气象条件在某些时段有利,在有些时段则为不利。2015—2019年京津冀13个城市气象条件对PM2.5、PM10和PM2.5-10长期分量的影响呈上升趋势,说明该时段气象条件有利于污染状况改善,但是至2021年底,气象条件的影响呈现下降趋势,且逐渐趋于“0”以下,说明这时气象条件从有利于转为不利于大气污染状况改善。

图5 2015—2021年京津冀13个城市经气象调整后PM2.5、PM10和PM2.5-10浓度长期分量序列以及气象条件对PM2.5、PM10和PM2.5-10浓度长期分量的影响

表2为京津冀城市群减排和气象条件分别对PM2.5、PM10和PM2.5-10浓度变化趋势的具体贡献。由于京津冀各城市气象条件和大气环境污染防治措施实施效果的不同,京津冀各城市气象条件和减排措施对PM2.5、PM10和PM2.5-10长期分量改变的贡献也不同。就PM2.5而言,气象条件的贡献介于3.23%～14.86%之间,减排措施的贡献介于85.14%～96.77%之间,其中,秦皇岛市气象条件对PM2.5的影响最大,占比为14.86%,衡水市减排措施影响最大,占比为96.77%。

表2 2015—2021年京津冀13个城市减排措施和气象条件对颗粒物浓度的贡献

就PM10而言,气象条件的贡献介于6.04%～13.69%之间,减排措施的贡献介于86.31%～93.96%之间,其中,秦皇岛市气象条件影响最大,占比为13.69%,廊坊市减排措施影响最大,占比为93.96%。就PM2.5-10而言,气象条件的贡献介于8.20%～26.23%之间,减排措施的贡献介于73.77%～91.80%之间,其中,沧州市气象条件影响最大,占比为26.23%,廊坊市减排措施影响最大,占比为91.80%。为了验证KZ滤波模型的分离成果,CHEN等[9]将基于WRF-CMAQ模型分离出的结果与基于KZ滤波模型分离出的结果进行对比,发现2种方法的结果相似:基于KZ滤波模型的结果显示,气象条件和减排措施对北京PM2.5浓度减少的贡献分别为19%和81%;基于WRF-CMAQ模型结果显示,气象条件和减排措施对PM2.5浓度变化的贡献分别为21%和79%。这2个模型都表明减排措施对大气颗粒物污染的改善至关重要。此外,张小曳等[4]和ZHAI等[13]研究结果也与之类似。这表明笔者研究结果具有一定可信度。

3 结论

基于KZ滤波法评估京津冀13个城市大气PM2.5、PM10和PM2.5-10变化趋势特征,以及气象条件和减排对变化趋势的贡献,主要结论如下:

(1)京津冀区域PM2.5年均浓度从2015年的75.32 μg·m-3迅速下降至2021年的39.82 μg·m-3,年均下降5.07 μg·m-3;PM10年均浓度从2015年的136.52 μg·m-3迅速下降至2021年的79.82 μg·m-3,年均下降8.10 μg·m-3。PM2.5和PM10污染天数比例从2015年的43%和31%分别下降至2021年的12%和8%。PM2.5、PM10和PM2.5-10浓度高值区主要分布在河北省南部邯郸、邢台、石家庄和保定等地,其变化趋势也最为显著。

(2)采用KZ滤波法对京津冀城市群13个城市大气PM2.5、PM10和PM2.5-10浓度时间序列进行分解,发现短期分量贡献最大,表明原始序列的波动主要是受污染物排放量、天气和季节变化的影响。京津冀城市群13个城市大气PM2.5、PM10和PM2.5-10浓度长期分量呈下降趋势。基于KZ滤波分析结果建立多元气象因子回归模型,得到气象条件和减排措施对颗粒物长期分量变化的影响规律。结果表明,京津冀地区长期分量变化趋势主要与污染物排放变化有关,说明2015至2021年京津冀地区实施的大气污染治理措施取得良好成效,同时,气象条件对大气环境改善也存在一定影响。

(3)2015—2021年减排措施对京津冀地区大气颗粒物污染状况改善起着至关重要作用。就PM2.5而言,气象条件和减排措施的贡献分别介于3.23%～14.86%和85.14%～96.77%之间,其中,秦皇岛市气象条件对PM2.5的影响最大,为14.86%,衡水市减排措施影响最大,为96.77%。就PM10而言,气象条件和减排措施的贡献分别介于6.04%～13.69%和86.31%～93.96%之间,其中,秦皇岛市气象条件影响最大,占比13.69%,廊坊市减排措施影响最大,占比93.96%。就PM2.5-10而言,气象条件和减排措施的贡献分别介于8.20%～26.23%和73.77%～91.80%之间,其中,沧州市气象条件影响最大,为26.23%,廊坊市减排措施影响最大,为91.80%。