COVID- 19疫情期间城郊森林公园O3变化对NO2减排响应的多重分形模式

张 娇, 刘春琼, 吴 波,杜 娟, 史 凯,*

1 吉首大学数学与统计学院, 吉首 416000 2 吉首大学生物资源与环境科学学院, 吉首 416000

2020年始,新型冠状病毒肺炎(COVID- 19)疫情在全世界蔓延[1- 2],中国各地区广泛实施了包括封闭城市在内的严格防控措施,以遏制疫情的蔓延[3- 4]。人类生产生活的急剧减少导致大气污染物人为排放大大降低[5- 6]。因此对生态环境来说,疫情的发生相当于进行了一场代价极为昂贵的污染减排实验,这也为探究自然生态系统中大气污染的本征演化规律提供了极为难得的机会[7]。

张家界国家森林公园拥有独特的原始自然生态系统和地貌景观,是世界自然遗产地和国家5A级风景名胜区,成为中外著名的旅游胜地,旅游产业已成为张家界主要的经济支柱产业[8- 9]。旅游旺季期间,随着旅游人数猛增,机动车使用量和餐饮油烟排放都将大大增加,导致二氧化氮(NO2)污染排放急剧增加。而淡季期间,人类旅游活动的减少使得NO2的排放大大降低。因此张家界大气NO2浓度波动的多尺度变异性特征,尤其是局部极端波动过程,蕴含着人类旅游活动的印迹和信息。同时,NO2是臭氧(O3)的重要前体物之一[10- 11]。直接排入大气中的氮氧化物NOx(NO和NO2)和挥发性有机物(VOCs)等在光照及高温条件下经过一系列复杂光化学反应生成二次污染物O3[12- 15]。2020年,由于疫情原因导致春季旅游旺季期间张家界游客锐减,城区与森林公园内NO2浓度也随之降低。因而疫情期间NO2锐减有可能直接影响着城区和森林公园大气O3的污染演变特征。对比分析疫情期间与过去同期大气观测数据,研究张家界不同生态功能区O3浓度变化对疫情期间NO2污染急剧减排的响应特征及演变动力机制,有助于进一步理解张家界不同生态功能区的O3生成机制差异,这将为森林生态系统功能维护及张家界景区污染排放控制策略的制定提供科学依据。

然而,O3生成是一个复杂的非线性光化学过程[16-17]。不仅与其前体物NOx和VOCs的排放有关,还受到气象条件、地形等下垫面复杂特征的影响,各种影响要素之间相互作用和关联[18- 19],使得O3浓度变化与NO2排放在不同时间尺度下存在高度非线性关系[20-23]。传统的统计学方法难以精准刻画O3及其前体物之间相互作用的复杂非线性特征。分形方法能从宏观、整体的角度表征和分析大气污染物随时间演化的内在标度规律及其复杂非线性特征[24- 25]。多重分形去趋势互相关分析(MFDCCA)[26]是多重分形中重要的分析方法之一,可为研究不同浓度波动下O3与其前体物之间互相关性的复杂非线性演化特征提供新的研究手段。目前,该方法已被应用到大气环境领域中。例如,Zhang等[27]基于MFDCCA方法研究发现北京和香港细颗粒物浓度与4个气象因子的互相关关系具有多重分形和反长期持续性特征。Shen等[28]基于MFDCCA方法研究发现南京市空气污染指数与不同气象因素之间的互相关关系具有不同的长期持续性特征。目前尚未将MFDCCA方法应用于研究森林生态系统内NO2与O3之间复杂非线性相关关系,应用该方法能有效的提取二者互相关的长期持续性等非线性特征。Bak等[29]创建自组织临界态(SOC)理论为阐释大气系统中污染物浓度演化的长期持续性主导机制提供可靠的工具。因此引入MFDCCA和SOC理论研究张家界不同生态功能区大气O3浓度演变对NO2减排响应特征具有重要意义。

为揭示张家界不同生态功能区O3浓度变化对疫情期间NO2污染急剧减排的响应机制,本文以疫情发生期间与过去同期张家界各监测点NO2、O3小时平均浓度数据为基础资料,应用MFDCCA方法对比分析不同时期O3和NO2之间相互作用的多重分形特征,并基于SOC理论阐明不同生态功能区NO2污染急剧减排对O3生成的非线性动力学影响。以期揭示自然生态系统大气环境质量随人类旅游活动变化响应机制,并科学评估人类旅游活动对自然景区的影响。

1 材料与方法

1.1 研究区域与对象

COVID- 19爆发之后,张家界采取了相对严格的交通管制措施,使得2020年3—5月旅游旺季期间旅游人数锐减。本文选取2020年3月1日—5月31日疫情期间与2019年3月1日—5月31日非疫情期间张家界不同生态功能区NO2、O3小时平均浓度数据为研究资料。数据来源于全国城市空气质量实时发布平台。由于停电、设备故障等非人为原因,研究时段内NO2、O3浓度数据有所缺失,故选取数据缺失较少的未央路、永定新区和袁家界站点进行研究,其缺失率均低于2.25%。其中,未央路监测点位于森林公园景区入口处,车流量较大,私家车和旅游客运汽车尾气排放是主要大气污染源；永定新区位于张家界市区内,交通排放和生活餐饮油烟排放是其主要大气污染源；袁家界监测点位于张家界国家森林公园内,森林覆盖率高达98%以上,天然植被类型多样化[9],主要反映了张家界自然植被生态系统大气成分背景特征。此外对于缺失数据,采用缺失数据前一时刻和后一时刻监测浓度值的算术平均值进行填补,如果前后时刻的数据也缺失,使用缺失数据前一天和后一天同时刻的数据取算术平均值代替。3个站点各年NO2和O3研究数据均为2208个。

1.2 NO2和O3时间序列基本统计量

NO2和O3浓度序列基本统计量如表1所示。受疫情影响,2020年未央路、永定新区和袁家界站点NO2平均浓度较2019年分别下降47.70%、2.21%、25.10%。这显示出疫情的发生造成了NO2的急剧减排。未央路站点NO2的平均浓度减幅远大于其他两个站点,这主要是未央路站点位于森林公园景区入口处,疫情期间旅游人数锐减,使得私家车和旅游客运汽车等尾气排放显著减少。此外,2020年未央路、永定新区站点O3平均浓度较2019年分别降低2.79%和3.40%,而袁家界站点O3浓度却升高了9.70%。COVID- 19期间前所未有的氮氧化物减排反而促进了自然森林生态系统中O3的生成。从峰度和偏度值来看,NO2和O3浓度时间演化并不服从正态分布,这暗示着NO2和O3浓度波动变化可能源于非线性动力系统。

表1 各监测站点2019—2020年3至5月NO2、O3浓度的基本统计量

各站点NO2、O3小时平均浓度的日变化规律如图1所示,2019年所有站点NO2日变化均呈现出明显的双峰状“M”型,这是典型人类旅游活动日变化规律造成的。永定新区NO2浓度最高,因为城区内机动车流量大使得交通尾气排放量高于其他站点。袁家界站点NO2浓度总体最低,这是由于张家界国家森林公园内主要采用电动车辆出行,受交通污染影响最小,其NO2主要可能来源于景区外排放源的污染输送。2020年袁家界站点NO2日变化规律几乎表现为一条直线,基本上不展示出人类旅游活动的影响。

各年各站点O3的日变化均呈现类似的单峰状特征。其中,未央路和永定新区站点O3日变化具有一致性,而袁家界站点O3日变化较为平缓,且浓度值明显高于其他站点。这可能与森林公园内NO2浓度日变化规律较为平缓以及森林系统中天然VOCs排放量较大有关。2020年袁家界站点O3浓度值较高,这是由于疫情期间氮氧化物的减排使得NO排放大幅减少,降低了对O3的滴定消耗作用。

由于疫情防控措施的实施,大气污染浓度发生了一定的变化,这说明污染物的排放和空气质量明显受到人类生产生活活动的变化影响。春季旅游旺季期间,张家界平均气温为11.8—12.2℃,平均降水量在18.1—22.1 mm之间,平均风力3级左右。在这两个研究时段内,张家界气象条件变化不大。因此通过对比2019年和2020年3—5月期间张家界3个站点O3和NO2的多尺度非线性相关作用的变化,可以突出反映疫情导致的人为旅游活动减少对不同生态系统功能区造成的影响。

图1 三个监测点2019—2020年春季期间NO2、O3浓度日变化规律Fig.1 Diurnal variation of NO2 and O3 concentrations in three monitoring stations in spring from 2019 to 2020

1.3 研究方法

1.3.1多重分形去趋势互相关分析(MFDCCA)方法

多重分形去趋势互相关分析(MFDCCA)由Zhou等在DCCA方法基础上建立的一种新的统计方法,用以识别两组非平稳时间序列之间在不同时间尺度上互相关性的多重分形特征[27,30]。

首先,通过两组非平稳序列x(i),y(i)重构序列X(i)和Y(i)：

(1)

其次,将新序列X(i)和Y(i)划分为互不重叠的等长区间Ns

(2)

其中s是时间尺度。

(3)

此时v=1,2,…,Ns,和

(4)

此时v=Ns+1,Ns+2,…,2Ns。

然后计算序列q(q≠0)阶互相关波动函数Fq(s),即

(5)

(6)

若x(i)=y(i),则MFDCCA等同于多重分形去趋势波动分析(MFDFA)。当q=2时,该算法等效于去趋势互相关分析(DCCA)方法。最后,若x(i)和y(i)具有长期幂律互相关性,那么s、Fq(s)及广义互相关指数h(q)应满足以下幂律关系：

Fq(s)∝sh(q)

(7)

当q=2时,DCCA标度指数h(2)与Hurst指数具有相似的性质。h(2)=0.5表示两序列之间没有互相关关系。h(2)>0.5代表两序列之间呈现出长期幂律互相关性,即一个变量极大值可能会在未来一定时间尺度内导致另一个变量极大值的产生,两变量之间的相关性具有长期持续性特征。h(2)<0.5则相反。

如果h(q)随着q的增加而单调递减,则表明两序列之间互相关性具有多重分形特征。多重分形的强度可通过h(q)的范围来计算,

Δh=Maxh(q)-Minh(q)

(8)

当Δh越大,多重分形性越强。即不同波动程度下,序列之间互相关性的长期持续性特征具有更强的多时间尺度变异性。与单一分形的长期持续性特征相比,多重分形结构赋予了多时间尺度上两个序列之间长期互相关作用子集的复杂变异性特征。不同时间尺度上,两个变量长期互相关作用子集相互嵌套和集聚,反映着复杂相互作用模式在多时间尺度演化的高度非线性和自相似特征。因此,h(2)和Δh参量为更准确和精细的刻画两个变量相互作用的非线性关系提供了一种异于传统随机过程中正态分布统计量的有效途径。

1.3.2自组织临界理论(SOC)

Bak等[29]于1987年创建自组织临界(SOC)理论,用于解释复杂系统中长期持续性特征的产生动力学机制。SOC系统包含许多发生非线性和短程近邻相互作用的组元,在外界物质及能量的输入驱动下,系统将自发地向临界状态演化。当一个系统达到临界状态时,系统内各组元相对稳定,整个系统表现出满足幂律分布的时空关联。外界发生一个微小的局部扰动可能会引发连锁反应,使该复杂系统发生各种规模的“雪崩”,从而导致大规模事件的发生,影响整个系统。这种临界状态下各组元间的非线性动态行为,服从具有长期持续性及标度不变性的幂律统计分布[31-32]。因此,长期持续性特征和幂律统计分布可以作为系统具有 SOC 特性的表征。如果O3浓度波动演化具有SOC特性,则意味着O3浓度波动满足幂律统计规律,

P(Δc)=A×Δc-λ, 即ln(P(Δc))∝-λln(Δc)

(9)

其中Δc(Δc=cn+1-cn)是O3浓度波动值,cn为n时刻O3浓度,A为无纲量参量,P(Δc)为Δc0所出现的概率,λ是衡量标度不变性结构的幂律指数。此时,O3浓度波动大小的发生频率随波动大小呈现幂指数下降,这种情况下O3浓度波动没有典型的特征浓度,具有标度不变性特征。这种特征与随机变量的正态分布结构完全不同,在具有正态分布的特征变量中,都存在典型特征值(即平均值)。

如果O3浓度波动还未达到自组织临界状态,而处于亚临界状态时,则不满足幂律统计分布规律。有研究认为,拉伸指数分布是系统处于亚临界状态的特征,O3浓度波动分布满足

P(Δc)=D×exp(-BΔcσ), 即ln(ln(P(Δc))∝σln(Δc)

(10)

其中B、D为无纲量参数,σ为拉伸指数。O3浓度波动的拉伸指数概率分布表明O3浓度演化尚未具有标度不变性结构,类似于正态分布,其某些典型特征浓度值出现的概率相对较大。

2 结果和讨论

2.1 NO2和O3互相关的长期持续性特征

首先应用MFDCCA方法对2019和2020年3个站点O3和NO2浓度互相关的长期持续性特征进行研究,计算结果如图2所示。q=2时,各站点lnFq(s)-ln(s)均满足线性关系。其斜率拟合结果见表2所示。各年各站点NO2和O3互相关性的h(2)值均大于0.5。这意味着NO2和O3互相关性存在较强的长期持续性特征。该特征表现为,在一定时间尺度内NO2和O3之间的互相关函数随时间的衰减不遵循经典的马尔可夫函数(Markov Function),即随时间变化的互相关函数不呈现指数形式的快速衰减模式,而是以一种缓慢衰减的幂律形式进行。这也意味着过去某一时段内NO2浓度的变化模式将在未来一定尺度内持续性的影响着O3的变化趋势。

图2 各站点2019—2020年NO2和O3波动MFDCCA分析的双对数波动曲线图Fig.2 The lnFq(s)- ln(s) between O3 and NO2 at different stations in 2019 and 2020, Zhangjiajie

从表2还可发现,2020年3个站点的h(2)值显著高于2019年。3个站点增高幅度分别为33.3%,17.1%和7.7%,这一现象可以结合疫情期间人类活动排放剧减进行理解。

疫情期间交通尾气排放的NO也大大减低,造成O3的滴定消耗随之大大减弱。同时,O3生成前体物NOx的急剧减少势必减弱了大气光化学反应的程度,这样O3及NO2在大气中将更加持久性的累积存在,从而造成过去某一时刻NO2将更加持续性的影响未来O3浓度变化趋势。这也是2020年NO2和O3之间互相关性的长期持续性增强的主要原因。同时,袁家界站点h(2)增值明显小于其他两个站点。这是因为和城区相比,森林公园内丰富的植被会排放大量的高活性生物质VOCs。即使在疫情期间人为源大幅度削减的情况下也可提供丰富的O3生成前体物质以继续维持相对活跃的大气光化学反应,削弱了O3及NO2在大气中的累积效应,导致袁家界站点h(2)增值不大。

表2 各站点2019年—2020年MFDCCA计算h(2)和Δh值

2.2 NO2和O3之间互相关性的多重分形特征

根据图2,各年各站点lnFq(s)-ln(s)拟合直线在不同的q值下斜率不同,这种不同波动程度下的长期持续性变异性需要进一步利用多重分形方法进行定量刻画。计算结果见图3。当q从-20到20变化时,h(q)随着q的增大呈单调递减趋势,这表明NO2和O3之间存在不同的非线性依赖关系,具有明显的多重分形性特征。

各站点NO2和O3互相关性的Δh计算结果见表2。各年Δh值均大于0,表明NO2和O3互相关性演化特征具有多重分形本质。从表2中可以看出两个显著性的规律。首先,从整体来看,2019和2020年袁家界站点Δh值比其他两个站点高。主要原因在于：袁家界站点植被覆盖率远高于未央路和永定新区站点。丰富的植物排放的植物源挥发性有机化合物(BVOCs)远超过人为源挥发性有机化合物 (AVOCs),占总挥发物的90%以上[33-34]。因此,袁家界站点BVOCs总量往往处于饱和状态,使得其NO2和O3之间的相互作用程度增强。另一方面,由于树种、树龄、气温等因素影响,不同树种BVOCs的化学成分、排放速率和排放量也有所不同[35-37],这将导致NO2和O3互相关性的变异性增强,从而袁家界站点Δh值比其他站点高。

其次,2020年3个站点的Δh值均低于2019年,下降幅度分别为2.6%,2.5%和12.1%。主要原因在于：疫情期间,NO2和VOCs等前体物浓度的大大减少,减弱了大气光化学反应的程度。这使得NO2和O3浓度波动在时间尺度上相对均匀,不同波动程度下二者相互作用的时间变异性减弱,从而2020年各站点Δh值相对于2019年有所下降。

为说明NO2和O3互相关的长期持续性在多重分形结构产生过程中的作用,本文进一步对原始NO2和O3序列进行随机重排处理,破坏序列内所有互相关性性质,保留其尖峰胖尾概率分布特征,产生NO2和O3波动的随机重排序列,并用MFDCCA方法对其进行分析,计算结果如图3所示。当q=2时,2019和2020年各站点随机序列的h(2)值接近0.5,表现出完全随机的特征,序列之间不存在内在的相关性。同时,随机和原始序列的多重分形结构差异巨大,这说明NO2与O3之间互相关的长期持续性是导致其多重分形结构形成的重要因素。

图3 各站点2019—2020年NO2和O3波动q阶广义Hurst指数图Fig.3 The generalized Hurst index between O3 and NO2 at different stations in Zhangjiajie

2.3 O3生成的自组织临界机制

上述讨论发现NO2与O3互相关的多时间尺度长期持续性特征对O3的生成发挥重要作用。然而究竟是什么宏观动力机制控制着长期持续性的形成？为了说明该宏观动力机制,本文首先对2019和2020年各站点O3浓度波动进行了频度统计分析。图4给出了各年各站点O3小时平均浓度波动的累积频率统计分布,具体分布函数关系见表3。袁家界站点O3的累积频率统计分布服从类似方程(9)的负幂律分布,O3演化处于自组织临界状态。此时高浓度O3波动和低浓度O3波动的出现概率相等。而未央路和永定新区站点的O3累积频率统计分布服从类似公式(10)的拉伸指数分布,O3浓度波动还处于亚临界状态,某些典型特征浓度值出现的概率相对较大。

图4 2019和2020年各站点O3小时平均浓度波动的累积统计分布Fig.4 Cumulative distribution of O3 pollution at different stations in 2019 and 2020

表3 张家界未央路、永定新区、袁家界站点O3波动累积统计分布拟合函数关系

袁家界站点O3浓度波动是开放、耗散大气巨系统,在前体物质作用下的复杂动力现象,其幂律统计分布的形成与演化既受到微观的大气光化学机制作用,同时也是多尺度宏观系统动力学长期相关作用的结果。

为进一步说明张家界国家森林公园森林生态系统中O3浓度波动涌现出负幂律统计分布的SOC 机制,将O3演化与自组织临界的沙堆系统进行类比分析。真实物理沙堆系统是SOC的典型范例[38]。在圆盘上通过逐粒加沙构造沙堆,当沙堆的倾角达到临界状态时,系统中加入的沙粒数量与沙堆崩塌掉落在圆盘之外的沙粒数量总体平衡。新添加的沙粒(系统的输入)可能停留在沙堆上,也有可能引起沙堆表面的小范围滑动,甚至造成更大规模的崩塌。宏观上,沙粒崩塌规模与其出现频率呈现负幂律统计分布[39- 40]。对于森林生态系统中的O3浓度演化来说,其前体物VOCs主要来自于植物排放的BVOCs, NO2主要来自于景区外人类活动中机动车排放。前体物之间的复杂光化学作用造成O3生成,这个过程类似于一个平板,持续投入沙粒,这些输入的物质或能量是森林生态系统中O3演化的直接驱动力。同时,森林大气系统中局域光化学作用导致的O3浓度上升将增强大气氧化能力,加剧光化学反应生成其他光化学产物(如二次气溶胶),该过程使得O3迅速消耗。二次气溶胶也可通过碰并凝结、降水洗刷等作用脱离大气系统。这正如沙堆达到临界倾斜角后,以崩塌方式清除多余沙粒,维持倾斜角的稳定临界状态。这过程中,O3与前体物相互作用的长期持续性可视为沙堆系统中局域沙粒之间相互挤压的短程相互作用扩散影响到整个系统的作用机制。最终,森林生态系统中,前体物之间非线性光化学作用持续累积生成的O3在一定时间尺度上并不产生持续稳定的O3浓度波动。相反,却产生类似于沙堆系统的规模大小遵循负幂律分布的非线性浓度波动演化。这样张家界森林生态系统中O3浓度波动演化过程与真实SOC沙堆系统在物理原型上非常相似。

对于景区以外的未央路和永定新区站点来说,由于城区森林植被覆盖远低于森林生态系统,再加上树种、树龄、气温等因素的影响,大气系统中VOCs的化学成分、排放速率和排放量将远低于森林生态系统,使得城区大气系统中O3前体物难以通过包含更多化学组分的复杂非线性光化学反应让O3浓度波动快速达到临界状态。同时,相对于森林生态系统来说,城区大气环境更容易受到人类活动的干扰。由于SOC机制是系统内禀机制,当外界力的干扰达到一定程度时,可能破坏系统内在的SOC特征。上述因素导致城市生态系统内的未央路和永定新区站点O3演化尚未达到自组织临界状态,而仅处于亚临界状态。

从复杂性理论来看,森林生态系统中O3浓度波动的长期动态可表征为多因素局部相互作用、相互影响所导致的宏观动力学效应。一方面,O3浓度波动表现为城市特定污染排放(包括人类旅游活动过程中交通、餐饮等污染排放)和自然环境要素(包括地形、气象、太阳辐射等)的综合影响,这些影响要素具有一定程度的确定性。另一方面,在开放的、耗散的大气光化学系统中,影响O3生成的因素之间相互作用非常复杂,在不同时间尺度上具有长期持续的非线性相关影响,因而导致O3生成在一定时间尺度上表现为不规则、非线性的变化。这些都是复杂系统的基本特性。长时间尺度宏观来看,各种确定性和不确定性的影响因素对森林生态系统中O3的生成产生了宏观“有序”的统计结构,即前述O3演化的幂律指数λ及其与前体物相互作用的DCCA标度指数h(2)和多重分形强度Δh。这些有序性的统计结构越稳定说明自组织演化的动力作用越强。根据2.1节的讨论,2020年袁家界站点NO2和O3之间互相关性的长期持续性较2019年显著增强,这就使得2020年森林生态系统中O3演化的自组织临界动力机制增强,从而导致高浓度O3的涌现,这是造成2020年袁家界站点O3浓度上升的主要动力原因。而其余站点由于仅达到亚临界状态,就无法通过SOC机制涌现形成高浓度O3。因此利用SOC动力机制可以更清晰的阐明2020年疫情期间张家界国家森林公园内O3浓度增高而城区站点O3浓度减低的动力学原因。

高浓度的O3具有很强的氧化性,会对森林生态系统中各种野生动植物造成生理伤害[41- 43]。众多学者已在野外观测到森林生态系统受到高浓度O3危害的症状[44- 46]。疫情期间张家界城郊型森林公园中O3浓度对NO2集中减排响应的SOC机制的准确识别,这对于高浓度O3的风险评估,有助于科学评估人类旅游活动对森林生态系统造成的影响。

3 结论

本研究基于2020年3月1日—5月31日疫情期间与2019年3月1日—5月31日非疫情期间张家界不同生态功能区NO2、O3小时平均浓度观测数据,探究了O3浓度变化对疫情期间NO2污染急剧减排的响应机制,结果表明：

(1)张家界森林及城区O3和NO2浓度波动的多时间尺度互相关性均存在较强的长期持续性特征,相对于2019年,2020年各站点NO2和O3之间互相关性的长期持续性分别增强33.3%,17.1%和7.7%,说明NO2在大气中将更加持久性的影响未来O3的浓度演化。

(2)张家界各站点NO2和O3之间互相关性存在显著的多重分形特征,其多重分形性的产生根源主要来自二者之间互相关性演化的长期持续性动力机制。2020年各站点NO2和O3之间互相关性的多重分形强度分别下降了2.6%,2.5%和12.1%,说明疫情期间NO2和O3之间互相关性的时间变异性减弱。

(3)袁家界站点O3浓度波动具有负幂律统计分布结构,未央路和永定新区站点O3频率统计分布服从拉伸指数分布。这表明森林生态系统中O3演化处于自组织临界状态,而城区生态系统内O3演化还处于亚临界状态。在NO2急剧减排的情况下,该内禀机制使得张家界国家森林公园内O3浓度增高,而城区站点O3浓度减低。准确识别NO2减排下O3演化的自组织临界特征,将有助于科学评估未来高浓度O3的发生风险。