基于回归和聚类的青岛市空气质量分析

Analysis of Air Quality in Qingdao Based on Regression and Cluster

摘要: 本文通过研究山东省青岛市空气质量指数和空气质量相关指标，分析得到青岛市主要污染物并对该污染物于2020年在青岛市的空间分布情况进行分析。首先通过描述性统计对2016~2020年青岛市空气质量情况进行分析，随后利用回归分析进行进一步的判断，综合分析后得出五年间青岛市空气主要污染物类别。而后利用聚类分析以及方差分析的方法，对主要污染物空间分布进行分析。本文最终得出的结论：近五年来，2018年青岛市空气质量最优；主要污染物为PM_2.5。全年空气污染较严重时期为冬季；青岛市空气污染较为严重的区域是黄岛区、市北区以及城阳区。

Abstract: This paper analyzes the main pollutants in Qingdao and the spatial distribution of the pollutants in Qingdao in 2020 by studying the air quality index and air quality related indicators. Firstly, the air quality of Qingdao from 2016 to 2020 was analyzed through descriptive statistics, and then the regression analysis was used to make further judgment. After comprehensive analysis, the main air pollutant categories in Qingdao in five years were obtained. Then, cluster analysis and variance analysis are used to analyze the spatial distribution of main pollutants. The final conclusion of this paper is that the air quality of Qingdao in 2018 is the best in the past five years. The main pollutant is PM_2.5 and the serious period of air pollution throughout the year is winter. The areas with heavy air pollution in Qingdao are Huangdao District, Shibei District and Chengyang District.

1. 引言

随着人类社会的发展，近年来人类对环境的过度干涉急剧增长，已经使得全球大气环境发生了较为严重的变化，对人类的生存和发展带来了现实性的危害和潜在的威胁。空气质量相关问题受到世界各国专家学者的关注。

刘丽等 [1] 认为春节期间烟花爆竹燃放等人为活动使空气质量变差；郭保林等 [2] 利用Mann-Kendall趋势检验法分析时间序列数据的可信度和变化趋势，认为2015~2019年五年期间东北地区各区域的空气质量指数以及空气中颗粒物(PM_2.5、PM₁₀)、SO₂、CO、NO₂浓度呈现出下降趋势；陈美如等 [3] 对秦皇岛市环境空气中主要污染物通过绘制折线图等描述性分析，认为近年以来秦皇岛市空气环境污染物以PM_2.5和O₃最为突出。周作明等 [4] 利用空气污染物负荷系数作为评价指标，认为湘潭市大气污染属于煤烟型污染，并且空气污染物的浓度呈现出逐年上升的趋势。Grange和David [5] 指出时间序列的明显特征可以被气象正态化程序所阐明，使得该技术具有灵活性，适合于空气质量干预勘探的一系列潜在应用。

多元回归分析，是基于空气质量监测站点数据与其它影响因素变量，借助最小二乘法建立的多变量回归建模方法 [6] 。尹杰 [7] 利用线性回归模型模拟PM_2.5的空间分布特征。方差分析是英国统计学家兼遗传学家Gitbert提出的一种分析方法，在科学试验、医学研究等众多领域有极其广泛的应用 [8] 。聚类分析同样广泛应用于多种领域，聚类可以分为层次聚类方法及系统聚类方法 [9] ，它是数据挖掘中的一种经典算法，将数据分别划分为多个类，具有较高相似度的对象被划分为一类，不同类之间具有较大的差别。其中，系统聚类法是较为常用的聚类分析方法。系统聚类的中心思想是确定点和点之间的距离，根据距离确定具体分类，距离测算可采用欧式距离、平方欧式距离等；类间距离还可使用最短距离法、最长距离法等 [10] 。

基于上述原因，深入了解、认识空气污染情况，探索青岛市主要空气污染物，准确、可靠的预测空气污染物情况是具有十分重要的意义的。本文最终所得到的结论可为青岛市大气污染防治提供科学依据，对进一步预测、控制及治理大气污染提供参考。

2. 方法

2.1. 线性回归

多元线性回归模型 $Y={\beta }_0+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p+ϵ,ϵ~N\left(0,{\sigma }^2\right)$ ，定义模型的总离差平方和、回归平方和与残差平方和为：

$S{S}_T={\displaystyle \underset{i}{\sum }{\left(y_i-\bar{y}\right)}^2},S{S}_R={\displaystyle \underset{i}{\sum }{\left({\hat{y}}_i-\bar{y}\right)}^2},S{S}_E={\displaystyle \underset{i}{\sum }{\left(y_i-{\hat{y}}_i\right)}^2}.$

$S{S}_T=S{S}_R+S{S}_E$ (1)

其中 ${\hat{y}}_i$ 为样本拟合值， $\bar{y}$ 为样本均值 [9] 。易证，即回归模型的方差解。

2.2. 单因素方差分析

单因素方差分析的模型为其中 $u_i$ 代表Y在因素A (共r个)

$\{\begin{array}{l}{y}_{ij}=u_{ij}+{\epsilon }_{ij},i=1,\cdots ,r;j=1,\cdots ,n_j,\\ {\epsilon }_{ij}~N\left(0,{\sigma }^2\right)且相互独立\end{array}$

水平第i个水平下的总体均值．引入r − 1个虚拟变量得到方差分析的虚拟变量回归模型 [10] ：

$A_i=\{\begin{array}{l}1,水平i\\ 0,其他\end{array},i=1,\cdots ,r-1,$

$Y={\beta }_0+{\beta }_1{A}_1+\cdots +{\beta }_{r-1}{A}_{r-1}+\epsilon ,\epsilon ~N\left(0,{\delta }^2\right).$ (2)

其中 ${\beta }_0={\mu }_r,{\beta }_i={\mu }_i-{\mu }_r\left(i=1,\cdots ,r-1\right)$ 。因素显著性可以用原假设 $H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_r$ 或 $H_0:{\beta }_1=\cdots ={\beta }_{r-1}=0$ 进行检验。由于因素A将数据Y分组，定义Y的总离差平方和、组间离差平方和、组内离差平方和分别为：

$S{S}_T\left(Y\right)={\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left(y_{ij}-\bar{y}\right)}^2}},S{S}_A\left(Y\right)={\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left(y_{i\cdot }-\bar{y}\right)}^2}},S{S}_E\left(Y\right)={\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left(y_{ij}-{\bar{y}}_{i\cdot }\right)}^2}},$

其中 ${\bar{y}}_{i\cdot }$ 为水平i下的样本均值．易证在方差分析的虚拟变量回归模型中 ${\hat{y}}_{ij}={\bar{y}}_{i\cdot }$ ，故模型的回归平方和与残差平方和为：

$S{S}_R=S{S}_A\left(Y\right),S{S}_E=S{S}_E\left(Y\right).$ (3)

2.3. 聚类分析

聚类分析(cluster analysis)是数据挖掘中的一种经典算法，将数据分别划分为多个类。具有较高相似度的对象被划分为同一类，不同类之间具有较大的差别。聚类分析广泛应用于多种领域，聚类可以分为层次聚类方法及系统聚类方法 [11] 。

2.4. Spearman秩相关

Spearman秩相关系数计算公式如式(4)所示。

$r_s=1-\frac{6}{n\left(n^2-1\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_j-Y_j\right)}^2}.$ (4)

该式中：n——时间周期的数量； $X_j$ ——周期j按时间序列排序的序号； $Y_j$ ——周期j内污染物浓度值按升序排序的序号。在Spearman秩相关系数分析中，如果得到 $r_s$ 的绝对值小于条件n下的临界值，则认为所得变化趋势不具有统计意义，反之则具有统计意义。若 $r_s$ 大于0则认为所研究数据具有上升趋势，反之则具有下降趋势。

3. 实证分析

3.1. 数据说明

根据我国相关政府文件，一个城市各点位24小时平均浓度的算数平均值作为该城市主要污染物当日的平均浓度值；各监测点位臭氧日最大8小时平均浓度值的算术平均值为该城市当日最大八小时平均臭氧浓度值。

本文数据来源为青岛市空气质量实时发布系统，数值单位为μg/m³ (CO的数值单位为mg/m³)。其中每日AQI数据和主要污染物浓度数据为当日各观测点24小时全部数据计算平均值的结果，符合国家相关文件统计方法要求并且保证数据准确性。本文研究所用数据在数据预处理过程中发现少量缺失值，因缺失值极少故采取剔除缺失值的处理方式，取均值后最终得到1806条有效数据。

3.2. 结果分析

3.2.1. 空气污染物变化趋势分析

为大致了解2016年~2020年青岛市空气质量变化情况，求得青岛市2016~2020年各年AQI值及主要污染物浓度值均值，绘制折线图如图1。可以看出，5年来青岛市AQI值及各个主要污染物浓度的年均值均在2018年达到最低，一定程度上可以认为五年中2018年青岛市的空气质量最优；2020年六类污染物浓度及AQI值均低于2016年，可一定程度认为2020年青岛市空气质量较5年前有所改善；五年间空气中SO₂浓度呈现逐年降低的较好趋势，2018年后下降速度变缓；AQI值、NO₂、O₃-8h、PM₁₀、PM_2.5的浓度总体呈现稳定，虽然有一定的下降趋势，但是其下降趋势不如空气中SO₂的浓度下降的趋势明显。

Figure 1. Change trend of AQI and major air pollutant concentrations in Qingdao from 2016 to 2020

图1. 2016~2020年青岛市AQI及主要空气污染物浓度变化趋势图

进一步分析2016~2020年青岛市空气质量的变化趋势，当时间周期为5时，统计量的临界值 $\gamma =0.\text{9}$ 。根据公式计算所得近五年青岛市6类环境评价监测指标的Spearman秩相关系数如下表1所示：

Table 1. Calculated Spearman Rank correlation coefficients for six monitoring indicators

表1. 计算所得6类监测指标的Spearman秩相关系数表

计算所得青岛市近5年来6类指标的Spearman秩相关系数如表1所示，可以看出，SO₂的 $r_s=-1$ ， $\left|r_s\right|>\gamma$ (0.9)；AQI、NO₂、O₃-8h的 $r_s=-0.\text{3}$ ， $\left|r_s\right|<\gamma$ (0.9)；CO、PM₁₀、PM_2.5的 $r_s=-0.7$ ， $\left|r_s\right|<\gamma$ (0.9)。故得到以下结论：2016~2020年青岛市空气中SO₂浓度呈显著下降趋势；空气中NO₂、O₃-8h的浓度有下降趋势，但并不显著；CO、PM₁₀、PM_2.5的浓度有下降趋势同样不显著，但其下降趋势相较于NO₂、O₃-8h而言，是显著的；另外，从AQI角度来看，5年来青岛市空气质量有提升趋势但并不显著。

综上所述，通过计算Spearman秩相关系数的方法对青岛市2016~2020年5年来空气质量变化趋势进行量化分析得到的结果与上节描述性分析所得的结论大致相同，可以得到以下结论：2016~2020年间青岛市空气质量得到较大幅度的改善；空气中SO₂的浓度呈现出显著的下降趋势；除SO₂以外的五种主要监测空气污染物的浓度均有所下降，但该趋势并不显著。

3.2.2. 空气主要污染物分析

因为空气质量状况与六类污染物的浓度状况全部相关，故本节利用Eviews软件构建多元线性回归模型，目的在于判断哪一种空气污染物的浓度变化更能影响当日的空气质量，影响越大，则认为该污染物为主要污染物。此处利用AQI值刻画每日空气质量状况。另外，由于研究过程中CO的浓度与其余五类污染物的单位不同，且AQI是无单位的，故本节构建多元线性回归模型的过程中将全部数据进行标准化处理，变成无量纲形式进行分析。将六类污染物全部引入作为自变量，AQI值作为因变量，依据最小二乘原理构建多元线性回归模型，输出结果如下表2：

Table 2. Multiple linear regression results under least squares

表2. 最小二乘下多元线性回归结果

表2中可以看出F统计量的值为2384.357，P值远小于0.05，所以认为模型通过检验，显著有效；观察每个变量的t检验结果发现，NO₂的t统计量值为−1.218，P值大于0.05，故该自变量不显著。调整后的R²约为0.888，大于0.7，可以认为模型拟合效果比较好。接下来检验模型是否符合模型设立的各个假设。首先判断是否符合各自变量独立的假设，检验多重共线性。若存在多重共线性，则有可能导致自变量的t检验失去意义、变量的系数值的含义不合理等情况，应考虑消除多重共线性。

多重共线性检验结果如下，此处根据方差膨胀系数VIF作为判断标准。由表3可以看出，自变量PM_2.5的VIF值大于10，所以是存在多重共线性的。所以该模型需要进行进一步的完善处理消除多重共线性。

Table 3. Results of multicollinearity test

表3. 多重共线性检验结果

首先利用向前逐步回归的方法再次构建模型，试图消除多重共线性，然而结果并不理想，仍然存在VIF大于10的自变量，所以又考虑到删除变量消除共线性。删除哪一个变量是当下面对的问题，经过上面的表可以看出PM_2.5的VIF较大，可以删除PM_2.5。然而，根据PM_2.5和PM₁₀的定义我们可知，PM₁₀包含PM_2.5且PM₁₀中接近70%为PM_2.5，故此处删除变量PM₁₀，再次拟合回归模型。

Table 4. Regression model after reducing variables

表4. 减少变量后的回归模型

根据表4可以看出模型通过了F检验，模型显著有效；五个自变量均通过了t检验，对空气质量水平有显著影响；调整后的R²为0.86，大于0.7，故模型拟合程度也比较好。根据表5判断是否存在多重共线性，发现所有变量的VIF均小于10，故多重共线性得以消除。

接下来对误差项进行检验，判断误差项是否符合正态分布，即是否符合误差正态分布的假定。根据图2残差图中下方的折线图可以看出，残差并无较明显的变化趋势，残差服从正态分布，符合模型假定。

Table 5. Multiple collinearity test results after removing variables

表5. 去除变量后多重共线性检验结果

Figure 2. Error chart

图2. 误差图

最后，对模型进行基于怀特检验异方差性检验。结果如表6，可以看出P值均小于0.05，所以拒绝原假设，认为存在异方差性，不符合异方差假设，所以对模型进行进一步修正。修正后模型表7。

Table 6. Heteroscedasticity test results

表6. 异方差检验结果

根据模型修正后的结果可以看出，五个变量均可以通过t检验，且模型通过F检验，模型拟合程度较好。表6中发现自变量PM_2.5的系数接近于0.818，大于0且是五个变量中系数最大的变量，故可以认为五个自变量中PM_2.5对于AQI值变化影响最大，并且其浓度越高，AQI值越大，空气污染越严重。一定程度可以认为PM_2.5是主要污染物。

Table 7. Model heteroscedasticity correction results

表7. 模型异方差修正结果

3.2.3. 主要空气污染物空间分布分析

经过上一节的研究可以知道，近些年来青岛市主要空气污染物为颗PM_2.5，本届主要研究PM_2.5在青岛市的时空分布，判断青岛市哪一个区域污染更加严重。本节所研究的数据为2020年全年青岛市九个观测点逐时观测数据取平均后得到的每日PM_2.5浓度以及AQI值。本文选七个观测站点进行分析，分别是：城阳区子站、李沧区子站、四方区子站、崂山区子站、市南区子站、市北区子站和黄岛区子站。分析过程中，根据国控监测点名称以及位置作为一个行政区的空气污染浓度，且由于老四方区是老工业集中区，具有较强的研究意义，故此处研究时暂不将其归于市北区的空气污染浓度。

利用青岛市七个观测站点2020年每日平均PM_2.5浓度的数据做聚类分析，目的是将含有相关关系的站点进行合并，达到数据降维的目的。分类情况整理如表8所示。

Table 8. System clustering results

表8. 系统聚类结果

为判断上小节分成的四类观测点的数值之间是否存在显著差异，本小节拟采用方差分析的方法进行进一步的判断。为判断数据是否能够进行方差分析，要判断其是否满足方差分析的假定，故首先对数据进行齐次性检验。由于数据原始分布未知，故此处采用的是列文检验的方法。根据上表可以看出，Levene统计量值为1.508，显著性为0.211，显著性大于0.05，故拒绝原假设，认为数据方差齐次，进行方差分析。

Table 9. Information table of analysis of variance results

表9. 方差分析结果信息表

表9是方差分析所得到的结果，可以看出方差分析的P值为0.015，小于0.05，故在显著性为0.05的水平下，认为所研究四类别之间存在显著差异，即各个观测站点分组后的新类别之间至少两组之间存在显著差异，但具体差异的组别并不知晓，故进行多重比较，此处采用最小显著差异法进行比较，结果如表10。

Table 10. Multiple comparison results based on LSD method

表10. 基于LSD方法下的多重比较结果

^*表示该值通过95%显著性检验。

可以清楚地看出每一类别与其余三个类别之间进行两两比较的结果，根据表中标示可以得到结论：在0.05显著性水平下第二类别和第一类别、第二类别和第四类别之间存在显著的差异，第三类别与其他类别均无显著差异。

根据表9可以看到，第二组是市北区站点、四方区站点、黄岛区站点；第一组是崂山区站点和市南区站点；根据表10可以看出第二组三个站点污染程度高于轻度污染的比例相对较高，一定程度认为该三个区域PM_2.5浓度相对较高，污染较为严重，第一组和第四组均与第二组存在显著差异，故认为市南区站点、崂山区站点以及李沧区站点所监测范围污染程度较低；而又因为第三组与第二组无差异，故认为李沧区站点所监测范围污染同样较为严重。

4. 结论

2016~2020年青岛市总体空气质量较好，空气质量水平高于山东省平均水平。六类监测污染物中，空气中SO₂的浓度连续5年下降，下降趋势最为显著；其余五类空气污染物浓度总体呈现下降趋势，但趋势并不明显。对该现象进行分析认为空气中SO₂浓度的持续走低与青岛市政府近年来大力推动清洁能源采暖实施密不可分。5年中2018年空气质量最优，对原因进行分析认为除了长期的防治举措以外，2018年大型会议举办期间，部分地区停工停产措施对当年空气质量改善起到十分重要的作用。

通过2016~2020年青岛市六类监测空气污染物的日均浓度以及每日AQI值进行可视化分析，得到如下结论：主要空气污染物和AQI值都呈现出来了明显的季节特征，SO₂、NO₂、CO、颗粒物(PM₁₀、PM_2.5)日平均浓度的变化趋势与AQI值变化趋势相同，在冬季较高，夏季较低；O₃的日平均浓度的变化特征与上述五种污染物相反，呈现出夏季和秋季明显高于冬季和春季的特点。

本文通过描述性分析、回归分析的分析方法最终得到2016~2020年青岛市最主要污染物为颗粒物(PM₁₀、PM_2.5)。对青岛市主要污染物即颗粒物进行空间分析，最终认为四方区、市北区、城阳区、黄岛区站点所监测范围的污染程度较严重，市南区、崂山区站点以及李沧区站点所监测范围污染相对不严重。一定程度反映各个行政区域受颗粒物污染程度。