自2013年国务院颁布《大气污染防治行动计划》以来, 我国大气污染防治工作取得了巨大进展^[1].近年来, 全国颗粒物浓度呈逐渐下降的趋势, 优良天数和比例增加, 重污染天数显著减少^[2].然而, 我国臭氧(O₃)浓度逐渐上升, 长时间、大范围的O₃污染事件屡屡发生, 影响了数亿人口的生产和生活^{[3, 4]}. 2019年, 全国337个城市年ρ(O₃)平均值为138 μg·m^-3[5], 其中, 超过一半的城市O₃浓度低于国家二级标准限值.O₃已成为部分城市空气质量不达标的首要因素, 是我国空气质量持续改善的重要问题^[6~8].

O₃是一种二次污染物, 其化学形成过程具有很强的非线性^{[9, 10]}.O₃浓度的变化高度依赖于气象因素, 具有局地、区域和季节尺度的特征.局地尺度看, 短期O₃浓度的变化和当地气象因素有关, 晴空、高温和低湿条件下, O₃浓度较高^[11~13].全国大部分地区臭氧浓度与温度和边界层高度具有显著的正相关关系, 与相对湿度和云量具有显著的负相关相关^{[14, 15]}.区域尺度, O₃浓度的形成受盛行风向、局地环流和大气边界层结构的共同影响^[16].风向和风速的变化对臭氧的传输和消散也会产生影响^{[17, 18]}.季节上看, 春季和秋季, 温度是我国O₃浓度变化的主导因子, 相对湿度对南方城市的O₃浓度有影响; 夏季, 部分区域O₃浓度变化受温度影响, 少部分区域受相对湿度影响; 冬季, 温度和湿度为主导的城市数量较少^[19].臭氧浓度还和前体物, 即氮氧化物(NO_x)和可挥发性有机化合物(VOCs)有很强的关系^[20~22].臭氧与前体物敏感性研究已受到国内外学者的广泛关注.有研究发现, 我国典型城市群区域VOCs控制区和NO_x控制区分界线对应的VOCs反应活性与NO_x反应活性的比值约为3∶1~4∶1^[23].不同时间段, 臭氧对前体物的敏感性差异较大.在上午, 京津冀、长三角和珠三角等地区臭氧生成的控制区为VOCs控制区, 在下午, 这些区域臭氧生成的控制区为NO_x控制区^[24].不同区域臭氧对前体物的敏感性也不尽相同.珠三角地区和台湾地区臭氧生成的控制区主要为VOCs控制区^[25~27].一般而言, 城市属于VOCs控制区, 远郊区属于NO_x控制区, 郊区则处于二者的过渡区^[28].此外, 天气系统和气候变化对臭氧浓度的变化也有一些影响.全球气候变暖影响局地对流层臭氧浓度, 对流层臭氧污染也反作用于气候变化^[29].热带气旋环流影响我国东南沿海臭氧浓度, 其外围大范围的下沉气流易引起高浓度臭氧污染事件^[30]; 东亚季风变化能够引起我国南部地区10%左右的臭氧年际差异^[31].

上述因素使得O₃的时间动态分析成为一项非常复杂的任务, 也导致了在研究复杂的和高度变化的时间序列时, 基于传统的统计分析方法描述O₃变化存在局限性.原因是这些方法通过拟合函数来获得近似和平滑的信号, 拟合过程中存在信息损失^[28].而复杂网络方法可以避免这一缺陷, 因为它们直接处理原始数据, 能够测度多元数据集内相互依赖结构的全部复杂性.

近20年来, 复杂系统科学研究领域最重要的发展就是在数学图论的基础上逐步完善了复杂网络的思想和方法.在复杂网络中, 网络节点是空间某一位置的时间序列, 连边是任意两个空间位置上时间序列的统计关联, 边的强度由节点之间关联性的强度量化^[32~35].复杂网络方法已经被国内外大气环境相关领域专家和学者所关注, 并在大气环境研究中得到良好应用. Zhang等^[36]的研究通过复杂网络方法建立了PM_2.5浓度和位势高度耦合网络, 揭示Rossby波对PM_2.5浓度变化的影响.Li等^[37]的研究构建京津冀地区PM_2.5浓度网络, 分析了京津冀地区PM_2.5浓度的空间的传输特征.肖琴等^[38]的研究应用复杂网络方法对中国城市PM_2.5浓度网络进行划分. Carmona-Cabezas等^[39]的研究应用复杂网络方法, 对安达卢西亚城市和农村O₃浓度进行区分.

复杂网络可以看作是研究大气环境相互关系的一种补充方法, 它继承了时间序列的特征, 最终通过度/加权度指标提供有用的信息.本研究将采用复杂网络的方法, 以城市为节点, 构建了包含全国335个城市的O₃浓度加权网络, 并分析近地面O₃浓度拓扑特征, 来深化对O₃浓度变化的认识.

1 材料与方法 1.1 研究数据

本研究使用的数据为2015年1月1日至2019年12月31日中国O₃的日最大8 h滑动平均浓度(O₃-8h, 下文O₃浓度均指O₃-8h浓度), 数据来源于中国环境监测总站全国城市空气质量实时发布平台(中国香港、澳门和台湾数据暂缺).经过对数据的处理, 共获得335个城市, 1 800 d时间序列.

1.2 研究方法

本文中, 以各个城市为节点, 城市间O₃浓度关系为连边.通过设置合适的阈值, 构建中国O₃浓度网络模型.

1.2.1 数据预处理

首先去除O₃浓度数据中包含的季节性变化, 处理后的数据记为T_i(l), 其中, i代表节点, i=1, 2, 3, …, 335, l代表总天数(d), l =1, 2, 3, …, 1 800.

1.2.2 计算节点间关联

网格上每对节点i和j之间的相关性连接可通过计算相关系数X_{i, j} (τ)来得到.计算方法参考已有研究^{[40, 41]}, 具体公式见式(1)和式(2).

(1)

(2)

式中, τ为时滞(以天数计). τ=-τ_max; -τ_max +1, …, τ_max -1; τ_max ×τ_max =30 d.为了避免边沿效应, 此处t的范围是1~l-τ_max.

计算正连接权重和负连接权重, 公式定义如式(3)和式(4)：

(3)

(4)

式中, max(X_{i, j})、min(X_{i, j})、(X_{i, j})和σ(X_{i, j})分别为τ范围内相关系数的最大值、最小值、平均值和标准差.除以σ(X_{i, j})是为了去除时间序列中由于自相关效应而导致的偏差. X_{i, j} (τ)的最大值定义为: P_{i, j}⁺ [P_{i, j}⁺ =max(X_{i, j})]. X_{i, j} (τ)的最小值定义为: P_{i, j}^－ [P_{i, j}^－ =min(X_{i, j})]. τ_{i, j}⁺和τ_{i, j}^－分别为X_{i, j}取得最大值和最小值所对应的时滞, 它们代表了正向连接或负向连接的方向.

1.2.3 显著性连边识别

为了识别显著性连边, 本文将原始数据进行随机重排(shuffling), 构建方法参考以往研究^{[41, 42]}.为了便于区分, 将基于原始数据构建的网络称为原始网络, 后者称为随机网络.

通过分析原始网络和随机网络的性质, 设置合适的阈值Q来分析原始网络中具有统计意义的连边.具体计算根据式(5)和式(6)：

(5)

(6)

当两个节点连接权重大于Q时, 将两个节点进行连接, 否则, 不建立连接.根据网络统计指标, 分析中国O₃浓度网络的拓扑特征.

1.2.4 加权度

度是复杂网络中最常见的统计指标.入度是指向节点的连接, 出度是指向远离节点的连接.节点i的加权入(或出)度被视为其传入(或传出)连接权重的数目.基于正向连接权重得到的加权入度表示为In⁺(w)_i, 加权出度表示为Out⁺(w)_i, 具体计算根据式(7)和式(8)：

(7)

(8)

式中, In和Out分别描述节点对其周围环境的依赖程度和其对周围环境的影响程度.基于负向连接权值得到的加权出入度类似于正向连接权重.

2 结果与讨论 2.1 中国O₃浓度网络性质分析

图 1描述了正/负向连接权重和最大/小相关系数和距离D_{i, j}的关系.从W_{i, j}⁺和D_{i, j}的关系可以看出, 原始网络的W_{i, j}⁺在近距离(< 2 000 km)较强. W_{i, j}⁺在随机网络中显著地小于其在原始网络, 且不随距离的改变而发生变化, 且99%的数据点主要分布在绝对值4以下. P_{i, j}⁺和D_{i, j}的分布图展示了随着D_{i, j}增加, P_{i, j}⁺呈现逐渐减少的模式; 并且在近距离内, 原始网络的P_{i, j}⁺显著高于随机网络.从W_{i, j}^－、P_{i, j}^－和D_{i, j}的分布可以看出, 原始网络和随机网络没有什么差别, 且不随距离的改变而发生变化, 这说明O₃浓度存在随机行为.

图 2(a)是W_{i, j}⁺的概率密度函数分布, 可以看出, 在[0, 4]范围内, 原始数据和打散重排数据都有一个峰值, 两者差别很小.这说明原始网络中的确存在一些随机连接.但在较高的值[4, 6]范围内, 原始网络的W_{i, j}⁺存在一个长的拖尾, 随机网络中不存在这种情况.说明这些较强的联系是由于某些确定物理相互作用而产生的, 而不是由于随机的因素导致的.真实网络和随机网络分布的差异性表明, 原始网络中存在许多显著的正相关关系, 在分析网络特征时, 需要去除随机连接, 以识别有重要意义的连接.从W_{i, j}^－的概率分布函数分布来看[图 2(b)], 原始网络和随机网络类似, W_{i, j}^－是一个随机行为.进一步分析τ_{i, j}⁺和τ_{i, j}^－的概率密度函数的关系可知[图 2(c)和2(d)], O₃浓度传播天数为3~4 d.

2.2 中国O₃浓度网络结构分析

通过上节分析, 随机网络中99%的W_{i, j}⁺主要分布在绝对值4以下, 因此本研究将4设为阈值.接下来, 主要考虑W_{i, j}⁺ ≥4的连边.图 3描述了加权度和连边长度的概率分布函数.可以看出, 加权度和连边长度的概率分布函数符合幂律分布, 拟合函数的相关性R>0.5.

为了研究每个节点对其他节点的依赖和影响关系, 本文分析了中国O₃浓度网络的加权入度、出度和其差异的空间分布.可以看出, 加权入度和加权出度的高值区均位于中国京津冀和周边区域, 和O₃浓度高值区分布类似.这说明浓度高的区域和周边城市关系紧密, 具有很强的传输能力.图 4(c)给出了加权入度和加权出度差异的空间分布, 可以看出, 内蒙古自治区和甘肃省O₃浓度受到周边区域的影响较大, 而北京、天津和河北东部(渤海湾沿线)城市可以将其O₃浓度传输到其他城市.具体地, 加权出度排名前10的城市为: 天津市、北京市、唐山市、沧州市、廊坊市、秦皇岛市、承德市、湖州市、常州市和南通市, 其中7个城市位于京津冀地区, 3个城市位于长三角区域.加权入度排名前10的城市为天津市、北京市、唐山市、湖州市、沧州市、常州市、南通市、廊坊市、无锡市、秦皇岛市和苏州市, 这些城市的出度值显著高于入度值.说明这些城市的O₃污染会影响到其他城市的O₃浓度.亟需在这些城市采取措施控制NO_x和VOCs的排放.

为了理解不同季节城市间O₃浓度传输和影响, 分析了不同季节下O₃浓度的加权出度, 具体如图 5所示.冬季, 加权出度值最高, 传输最强的城市主要集中在东部城市区域.尽管夏季O₃浓度值较高, 夏季城市间O₃浓度传输作用小于冬季, 春季和秋季加权度的强度较小.冬季传输能力强可能和较强的冬季风有关.当冬季风速偏大时, 有利于O₃污染严重地区污染物向周围城市扩散, 这加强了城市间O₃浓度的相关性, 进而增加了O₃浓度加权出度值.

全年来看, O₃浓度传输路径主要为西北-东南向传输[图 6(a)], 这反映了O₃浓度变化受到冬季风的影响.冬季, O₃浓度传输有两条典型的传输路径[图 6(b)]：①从内蒙古自治区向西南传输到京津冀地区; ②从京津冀地区向南传输到长三角区域.其中, 北京和天津O₃污染和河北省O₃传输有关系, 这和以往北京地区O₃重污染时段研究结果类似^{[43, 44]}.安徽省、山东省、河南省和河北省对长三角区域O₃污染有一定影响, 这和以往O₃污染过程中区域传输分析结果一致^{[45, 46]}.总体上, 全年和冬季传输路径和冬季风有关, 反映了风向对O₃浓度传输的影响.

3 结论

(1) 基于2015~2019年中国O₃浓度构建的原始网络中正向连接权重、最大相关系数在近距离显著高于随机网络中的正向连接权重.而原始网络中负向连接权重、最小相关系数和随机网络的差别不大.这说明O₃浓度网络既包含真实情况, 也包含随机的特征.

(2) 显著统计意义的连边的分析发现, 加权度和边长的概率分布函数符合幂律分布, 说明O₃浓度变化有一定规律性.应用加权度指标分析了O₃浓度的加权出度和加权入度的空间分布, 结果表明, 加权度的高值分布和O₃浓度高值区分布一致, 这说明浓度高的区域和周边城市关系紧密, 具有很强的传输能力.

(3) 北京、天津和河北东部(渤海湾沿线)城市向其他城市的传输能力较强, 并且包含很强的季节特征.冬季城市间O₃传输能力最强, 夏季次之, 春季和秋季传输较小.

(4) 全年和冬季传输路径分析发现, O₃浓度变化和冬季风有很强的关系.