川渝地区是我国空气污染最为严重的地区之一, 准确认识区域内空气质量时空格局及形成机制对该地区环境保护政策制定具有重要意义。早期川渝地区空气质量格局研究多集中在成都和重庆两大核心城市^[1-3], 随着空气质量监测网络的逐步完善, 部分学者展开了更大范围的时空格局分析^[4-6]。另外, 纵观已有研究, 学者们普遍认为逆温^[7]、相对湿度^[8]、边界层高度^[9]、地形^[10]、风速^[11-12]、风向^[13]、气温^[14]、太阳辐射^[14]、降水量^[15]等自然因素可对盆地区域的空气质量产生明显影响。

但是上述研究主要基于城市站的监测数据, 受监测站位置的限制, 这种有偏样本难以准确反映川渝地区西部高原与盆地偏远地区的空气质量时空分布格局。而对空气质量影响因素的剖析也多基于重污染事件, 注重气象和地形条件等自然因素的影响^[7-15], 对自然与人为因素的综合性分析较为缺乏, 给准确认识川渝地区空气质量的形成机制带来较大局限。因此, 本文拟基于2018年和2019年覆盖川渝地区全部县域的199个空气质量监测以及自然与社会经济数据, 采用空间分析方法从较为宏观的尺度研究川渝地区全域的空气质量时空格局, 并确定其主要的自然与人为影响因素, 从而为环境管理决策提供依据。

1 数据与方法 1.1 数据来源与数据预处理

(1) 研究区概况

川渝地区行政上包括四川省和重庆市, 整体可分为东部盆地区和西部高原区两大自然地理单元。西部高原区主要位于攀枝花市、阿坝州、甘孜州及凉山州4个地市州, 其余18个地级市多位于东部盆地区。东部盆地区多为平原和丘陵, 四面环山；年平均气温为17.1℃, 夏季最高为26.5℃, 冬季最低为7.9℃；年降水量约850mm, 雨季为7-9月。年均风速1.7m/s, 年均相对湿度78%, 逆温现象时常发生；区内人口密集, 工业发达, 能源消耗和排放量高。西部高原区海拔较高, 全年有7个月为旱季, 日照丰富, 降水量较少；经济活动以畜牧业为主, 只有攀枝花市和凉山州拥有相对发达的工业(图 1)。

(2) 数据来源与预处理

本文采用六项基本空气污染物(PM_2.5、PM₁₀、SO₂、NO₂、CO、O₃)2018年和2019年的逐日监测数据, 共涉及128个城市站和71个县级站。城市站数据来源于中国环境监测总站全国城市空气质量实时发布平台(https://quotsoft.net/air/), 县级站数据来源于国家重点生态功能区县域生态环境质量评价项目。

首先, 根据《环境空气质量标准(GB 3095-2012)》对空气质量监测数据有效性的规定, 对获得的空气质量监测数据去除无效数据, 并剔除离群值。然后依据《环境空气质量指数技术规定(试行)(HJ633-2012)》, 利用各污染物的质量浓度数据计算空气质量指数(Air Quality Index, AQI)以及各项污染物的空气质量分指数(Individual Air Quality Index, IAQI)。

其次, 获取各自然指标栅格数据。气温、降水量、风速、相对湿度等气象站点数据来自于中国气象数据网(http://data.cma.cn/), 包含川渝地区及周边181个气象站点的逐日监测数据。使用ANUSPLIN软件, 利用高程作为协变量对气象数据插值^[16]。高程数据来自资源环境科学与数据中心网站(http://www.resdc.cn/)。年均太阳辐射数据来源于ECMWF气象再分析数据(https://cds.climate.copernicus.eu/)。用于表征植被覆盖状况的归一化植被指数(NDVI)数据来源于MODIS网站(https://modis.gsfc.nasa.gov/data/dataprod/mod13.php)。

然后, 确定各市区县空气质量代表值以及站点代表性位置。由于城市站在部分市区有多个站点分布, 为了减少由于站点较多对区域总体统计量估计造成的偏差, 本文对市区内的多个监测站点数据取平均值以代表某市区的空气质量总体状况, 并以监测站位置中心代表某市区的监测站代表性位置。对于其他区县, 则采用原监测数据及监测站位置作为这些区县空气质量代表值以及监测站代表性位置。最终纳入统计分析的站点为86个, 其中盆地区41个, 高原区45个。

最后, 基于各自然因子指标栅格数据, 以监测站代表性位置周围15km为缓冲区, 统计各自然因子的区域均值作为该位置的指标值。人口城镇化率、人口密度、第二产业占比和人均GDP等人为因素指标为县域尺度统计值, 数据来源于《四川统计年鉴》^[17]、《重庆统计年鉴》^[18]、《中国县域统计年鉴》^[19]和《中国城市统计年鉴》^[20]。

1.2 研究方法 1.2.1 空间格局分析

采用全局莫兰指数(Global Moran′s I)和局部莫兰指数(Local Moran′s I)表征空气质量(分)指数的全局和局部聚集性特征^[21]。全局莫兰指数的计算公式如下:

(1)

式中, n为空间单元的数量, x_i、x_j分别是单元i、j的AQI值或IAQI值, x为所有空间单元的平均值。W_ij为空间单元i、j的空间权重矩阵。采用标准化统计量Z检验空间自相关是否显著, 计算公式为:

(2)

式中, E(I)表示I的均值, Var(I)表示I的方差。

局部莫兰指数的计算公式如下:

(3)

式中各符号含义与公式(1)相同, 显著性检验也由公式(3)计算。若研究指标在单元i上有聚集性, 则可能存在高值-高值聚集、高值-低值聚集、低值-高值聚集和高值-低值聚集四种聚集类型。

1.2.2 自然与人为影响因素的综合分析

由于常用的空气质量影响因子间普遍存在多重相关性, 然而在建模时回归系数估计值对样本数据的微小变化敏感, 给回归系数的统计检验造成了一定困难, 回归系数的一般解释将不再适用。因此, 本文采用偏最小二乘回归(Partial Least Squares Regression, PLSR)分析空气质量的影响因素, 该方法集成了多元线性回归、主成分分析和典型相关分析的优点, 能够在自变量存在严重多重相关性的条件下进行回归建模^[22]。

PLSR模型各因子的回归系数反映了该因子对因变量的影响大小与方向, 变量映射重要性(Variable Importance in Projection, VIP)则反映了自变量对因变量的解释力, 自变量的VIP值越大, 对因变量影响能力越强, 本文将VIP的阈值设置为1, 大于该数值代表自变量对因变量有显著影响^[23], 并采用R²反映回归方程的可靠程度。获得各因子的回归系数后, 计算某个(类)因子回归系数绝对值相对于全部因子的占比, 以确定某个(类)因子的解释能力。限于篇幅本文因变量只选择了三种污染较为严重的污染物(PM_2.5、PM₁₀和O₃)的IAQI值。

2 研究结果 2.1 时空分布格局

川渝地区空气质量整体为良, 西部高原区空气质量明显优于东部盆地区。川渝地区全区AQI为55.2±26.5, 其中盆地区为63.9±31.2, 高原区为47.2±18.1；主要污染物为O₃, 其次为PM₁₀和PM_2.5, NO₂、CO和SO₂对空气质量影响较小。盆地区和高原区主要污染物构成存在较大差异。盆地区主要污染物为PM_2.5, 其次为PM₁₀和O₃, 但三者较为接近, 其均值分别为47.8、47.4和44.3；高原区主要污染物为O₃, 其次为PM₁₀和PM_2.5, O₃明显高于其他两项, 其均值分别为41.1、30.4和21.9。

盆地区与高原区的AQI在季节尺度上存在明显差异。盆地区AQI在冬季最高(81.7), 而春、夏、秋分别为65.6、59.2和49.0。高原区AQI在春季最高(54.5), 而冬、夏、秋分别为48.9、45.8和39.9。两区域内除O₃外各污染物的IAQI都表现出大体一致的季节分异。PM_2.5、PM₁₀、CO以及盆地区的NO₂都表现为冬季最高, 随后依次为春季、秋季和夏季, 只有高原区的NO₂表现为秋冬季略高于春夏季。O₃的IAQI在盆地和高原区都表现为春夏季较高, 秋冬季较低。盆地区O₃的IAQI在夏季略高于春季, 而在高原区春季明显高于夏季。SO₂无论在盆地还是在高原区不同季节间差异均较小。

在月尺度上, PM_2.5、PM₁₀基本为“U”形曲线, 每年1月IAQI较高, 随着气温回升逐渐降低, 到夏季7月份降为最低, 然后逐渐升高, 一直持续到年底(图 2)。O₃则表现为倒“U”型变化, 盆地区O₃峰值在4-8月, 持续时间较长, 高原区峰值在4-5月, 持续时间较短。NO₂和CO都呈现宽浅的“U”形, 只在冬季明显升高, 其他季节变化不大, 而SO₂则全年都没有明显的波动。盆地区AQI变化为颗粒物主导型, 曲线形状与PM₁₀和PM_2.5基本一致, 只是在4-8月份受O₃影响, U形曲线底部更为平缓, 并在4月和8月出现局部的峰值。高原区AQI变化则为O₃主导型, 总体上也呈上凸曲线, 3-10月份曲线形状与O₃一致, 只是在冬季月份由于颗粒物的增加, 曲线两端翘起。

2.1.3 空间分布格局

川渝地区空气污染具有明显的空间聚集性特征(图 3)。AQI以及PM₁₀、PM_2.5、O₃、SO₂、NO₂、CO的IAQI全局莫兰指数分别为0.57、0.53、0.60、0.32、0.39、0.27和0.31, 都呈现出显著的聚集性特征(P < 0.01)。除AQI外, PM₁₀和PM_2.5的全局莫兰指数最高, 这表明在川渝地区内, PM₁₀和PM_2.5较其他污染物更具聚集性。从局部莫兰指数所表现出的聚集性来看, 川渝地区的AQI以及各项污染物的IAQI多表现为高-高聚集和低-低聚集特征, 高-低聚集区和低-高聚集区较少。

AQI表现为高-高聚集和低-低聚集同时共存, 高-高聚集主要分布在盆地西部和中部平原的成都、眉山、乐山、自贡、内江、资阳等地, 低-低聚集区主要分布在甘孜州和阿坝州的中南部。PM₁₀、PM_2.5和NO₂与AQI的空间聚集性类似, 但除上述的高值聚集区外, 在大巴山区东部的部分县市也出现了小规模的聚集区；低值聚集区也比AQI范围更广, 几乎囊括了甘孜州和阿坝州的所有县市。O₃则多表现出高-高聚集性, 低-低聚集区很少出现, 高-高聚集区主要出现在四川盆地的西部、中部和南部。SO₂在攀枝花市和凉山州具有明显的高-高聚集区。CO的空间聚集性相对不明显, 仅在盆地东北边缘有局部高值聚集。

2.2 主要影响因素分析

本文所选择的各影响因素指标对三个主要污染物的IAQI都具有较好的解释能力, 自变量可以解释IAQI变化的43%-75%。其中对PM_2.5变化的解释能力最好, 可解释75%；其次为PM₁₀, 可解释61%；对O₃变化解释能力最小, 为43%, 这说明可能存在其他因素影响着O₃浓度(表 1)。

对PM_2.5而言, 所有自然因素指标的解释力占比为78.2%, 人为因素为21.8%。对PM_2.5能产生重要影响的指标主要包括年均风速、年均相对湿度、年均太阳辐射、年均气温和城镇化率, 各项指标的解释力占比分别为23.2%、15.5%、14.1%、8.5%和10.7%, 风速和太阳辐射越低PM_2.5的IAQI越高；气温、相对湿度与城镇化率越高PM_2.5的IAQI越高(表 1)。

对PM₁₀而言, 全部自然因素指标的解释力占比为90.7%, 人为因素只有9.3%。其中年均风速和年均相对湿度的影响最为突出, 解释力占比分别为36.3%和28.6%, 风速越小PM₁₀的IAQI越高；相对湿度越高PM₁₀的IAQI越高。虽然年均气温和年均太阳辐射对PM₁₀也有重要影响, 但解释能力较低, 分别只有4.4%和2.7%(表 1)。

O₃的IAQI变化全部自然因素指标的解释力占比为38.0%, 人为因素为62.0%。其中, NDVI的影响最为突出, 解释力占比为21.7%, 三种人为因素对O₃的IAQI解释能力分别为20.7%、18.5%和17.4%。NDVI越高O₃的IAQI越低；城镇化率、人口密度和人均GDP越高, O₃的IAQI越高(表 1)。

3 讨论 3.1 基于城市站数据评估结果显著高估了川渝地区空气质量指数

由于城市站多位于地市级城市的市区/郊区, 经济相对发达, 比其他区县的污染物排放量大, 基于城市站数据的分析结果往往显著拉升了区域整体的AQI, 无法全面了解空气污染状况^[24]。本文采用独立样本T检验方法, 检验了只采用城市站和采用全部站点(包括城市站和县级站)计算的AQI和各污染物IAQI的均值差异。在川渝全区、盆地区和高原区采用城市站计算的AQI显著高于采用全部站点的结果, 分别高出31.9%、15.6%和32.4%。PM_2.5、NO₂和O₃的IAQI也表现出同样的规律, 在全区分别被拉升了54.4%、67.0%和25.8%, 在盆地区分别被拉升了14.0%、22.3%和20.3%, 在高原区分别被拉升了68.9%、137.1%和38.0%。而其他污染物的IAQI则表现出不同的规律, SO₂和CO都在盆地区被显著拉低, 在高原则相反。PM₁₀则只在全区和盆地区被显著拉升, 但在高原区在统计上未表现出显著的差异(表 2)。

3.2 川渝地区空气质量影响机制

对PM_2.5来说, 高风速是颗粒物稀释和扩散的重要驱动力^{[8, 25]}, 本文研究结果表明在川渝地区诸多影响污染物时空分布特征的因素中, 风速对PM_2.5的IAQI影响最为突出。强太阳辐射可使近地面层热力对流旺盛, 逆温出现的几率减小, 易于形成有利的扩散条件^[26], 实际上本文研究结果表现的太阳辐射越低其IAQI越高, 本质上也反映了空气扩散条件对PM_2.5的IAQI影响。相对湿度和气温对PM_2.5的IAQI影响主要与这两个因素对PM_2.5的形成能够产生促进作用有关。高相对湿度下气溶胶吸湿增长, 有利于二次颗粒物的生成^[9], 温度升高也可加速PM_2.5的二次产生^[10], 因此川渝地区气温与相对湿度越高其PM_2.5的IAQI越高。人口密度既能表征工业源也能反映生活源的污染物排放能力, 但污染物排放能力对川渝地区PM_2.5污染的影响较自然因素影响要小, 因此确定准确的污染物大气容量对川渝地区空气污染的治理工作也就显得更为重要。

对PM₁₀来说, 风速越小其IAQI越高, 相对湿度越高其IAQI越高, 这可能反应了川渝地区PM₁₀来源的内源特征, PM₁₀污染可能主要与低风速高湿度条件下颗粒物的积累及二次转化过程有关。在风速低、湿度高的气象条件下, 颗粒物多附着、溶解或混合于雾气中, 污染物扩散条件差, 颗粒物容易积累^[27]。另外, 川渝地区PM_2.5和PM₁₀浓度的相关性较强, 在高原区与盆地区内相关系数分别为0.78(R²=0.52, P < 0.001)和0.63(R²=0.83, P < 0.001), PM_2.5在全部PM₁₀中的平均占比分别约为49%和62%。低风速高湿度条件下, PM_2.5容易与其他污染物发生化学反应, 增加二次气溶胶产生的可能性^[27]。因此, 在不利于颗粒物扩散但却利于其积累与二次转化的低风速与高湿度双重自然因子作用下, PM₁₀的IAQI表现为自然因素控制为主的特征, 与人为排放的相关性就表现得不再突出。

对O₃来说, 城镇化率、人口密度和人均GDP越高, 可在一定程度上反映区域内汽车保有量较大^{[28, 29]}, 汽车保有量大意味着前驱气体(VOCs和NO_x)的排放量较多, 从而有利于O₃的产生。而O₃沉积在植被外表面是地面臭氧清除的重要途径, 同时O₃也可从气孔进入植物体内被吸收^[30], 因此植被状况(NDVI)对O₃的产生呈现出显著的负向影响。虽然植物源VOC的主要成分异戊二烯和单萜烯是O₃的重要前体物^[31], 但据Qiao等^[32]对四川盆地西缘森林的研究, 植被源O₃占其总量最多为15%, 因此植被本身排放的萜烯类等VOC对O₃浓度影响并不明显。低风速、低相对湿度、高气温、强太阳辐射有利于O₃进行光化学反应及其前体物的积累^{[5, 33]}, 而本文研究结果显示, 上述因子对O₃的IAQI并没有显著影响, 这可能与研究的时空间尺度不同有关, 在不同时空尺度下, 各影响因素对空气污染的影响机制会有所不同^{[33, 34]}。

3.3 研究结果的不确定性与研究展望

本文研究结果表明, 国家重点功能区作为限制性开发区, 虽然人为活动相对较弱, 但在川渝地区的相关县域O₃污染已经显现, 需要引起足够的重视。相关县域较高的O₃浓度, 除了受人为排放影响外, 与区域植被状况较差也可能有密切关系。因此, 在川渝地区国家重点功能区相关县域O₃污染的防治工作中, 除了控制经济开发规模从而减小排放外, 加强植被保护与生态修复也非常重要。

限于县域数据的可得性, 有些可能对污染物时空分布产生影响的因素如逆温、大气边界层高度等自然因素没有考虑, 污染排放与治理等因素采用了城镇化率、第二产业产值与GDP等宏观社会经济指标来评估, 这都可能给分析结果带来一定的偏差。另外, 由于县级站监测数据时间较短, 目前仅有2018~2019年数据较为完整, 因此, 本文空气质量影响因素分析主要采用空间梯度法基于年尺度进行。然而不同时间维度上的空气质量形成机制可能存在差异。随着县级站空气质量监测数据的积累, 未来研究可进一步深入挖掘川渝空气质量的时空分布特征及形成机制。

4 结论

本文通过对比利用川渝地区全部县域和只采用城市站数据的空气质量评价结果, 发现只基于城市站数据的评估结果显著拉升了川渝地区的空气质量指数(AQI), 在全区、盆地区和高原区分别拉升了31.9%、15.6%和32.4%。盆地与高原区虽然O₃大致都呈倒“U”型变化, 但盆地区峰值期为4-8月, 高原区峰值其只发生在春季的4-5月；SO₂无论在盆地还是在高原区不同月份差异较小。AQI及各污染物均有显著的空间聚集性特征, 空间上多呈高-高聚集或低-低聚集特征。高-高聚集区多分布在盆地西部、中部和南部, 低-低聚集区多分布在高原区。基于偏最小二乘回归方程各自变量的回归系数大小、方向及其对因变量影响显著性等分析表明, PM_2.5受风速、太阳辐射、气温、相对湿度与城镇化率等影响突出, PM₁₀受风速和相对湿度影响明显, O₃与城镇化率、人口密度和人均GDP越高以及植被状况(NDVI)关系密切。盆地区PM_2.5污染可能与高强度经济活动叠加低风速与高温高湿等气候条件有关；PM₁₀污染可能主要与低风速与高湿度条件下颗粒物的积累及二次转化过程有关, 因此PM₁₀浓度在盆地区自然会明显高于高原区。高原区较高的O₃浓度可能主要与植被状况较差有关, 盆地区可能与较高较强的人为活动有关。