基于两阶段自适应Lasso方法的乌鲁木齐市PM_2.5浓度建模研究

The Modeling Study of PM_2.5 Concentration in Urumqi City Based on the Two-Stage Adaptive Lasso Method

摘要: 细颗粒物浓度是评估空气质量的重要指标之一，精确预测PM_2.5浓度仍是国内外研究的挑战，也是政府和学术界关注的重点。文章以2014年1月1日至2024年9月30日的乌鲁木齐市PM_2.5日平均浓度为研究对象，对其线性回归模型应用两阶段自适应Lasso多变点检测与估计方法检测并估计回归系数的多变点，建立了更精准的分段线性回归模型来刻画乌鲁木齐市PM_2.5日平均浓度。

Abstract: The concentration of fine particulate matter (PM_2.5) is one of the key indicators for assessing air quality. Accurate prediction of PM_2.5 concentration remains a challenge in both domestic and international research and is a major focus of attention for governments and scholars. This study takes the daily average PM_2.5 concentrations in Urumqi City from January 1, 2014 to September 30, 2024 as the research object. A two-stage adaptive Lasso multiple change-point detection and estimation method is applied to detect and estimate multiple change points in the regression coefficients of the linear regression City daily average PM_2.5 concentrations in Urumqi City.

1. 引言

随着全球经济和工业化的迅速发展，空气污染问题日益严重，已对居民的日常生活和幸福感产生了显著影响。2021年是“十四五”规划的开局之年，各省和地区纷纷出台生态环境保护规划，明确以改善生态环境为核心，推动绿色发展。而空气质量问题成为其中亟需解决的关键问题。细颗粒物(Particulate Matter 2.5，PM_2.5)浓度是评估空气质量的重要指标之一，也是城市空气质量的主要污染物之一。其来源可分为自然源和人为源，自然源包括沙尘暴、细菌、动植物屑片等；人为源则主要来自交通、灰尘、燃煤和工业排放。PM_2.5属于细颗粒物，指直径小于或等于2.5微米(μm)的固体颗粒。由于PM_2.5颗粒极其细小，又在空气中悬浮时间较长，可渗入肺部甚至血液，流行病学研究表明，长期生活在PM_2.5浓度较高的环境中会增加心血管疾病的风险甚至引起死亡；此外，PM_2.5浓度过高还会引发雾霾，影响交通安全，阻碍城市发展。为减轻空气污染对身体健康和社会造成的危害，保障公众健康和日常出行，精确预测PM_2.5浓度并采取合理措施是至关重要的。PM_2.5浓度的准确预测是国内外研究面临的重大挑战，也是政府和学术界关注的重点，近年来，学者们在PM_2.5预测的理论研究和实证分析方面已开展了大量研究。2012年，朱少钧等人利用美国空气资源实验室的混合单粒拉格朗日积分轨迹(Hybrid Single-Particle Lagrangian Integrated Trajectory, HYSPLIT)模型进行颗粒物的溯源和追踪分析，并做了相关预测[1]；2014年，蒋雪峰等人通过建立多元线性回归模型和shepard二维插值模型分析了PM_2.5浓度的演变规律并进行了预测[2]；2015年，胡玉筱与段显明结合气象因素和其他空气污染物构建多元线性回归模型预测西安市PM_2.5浓度的演变[3]；2015年，刘杰等人应用多元回归模型和机器学习模型对PM_2.5浓度进行了预测[4]；2018年，盛永财等人通过相关分析和主成分分析方法探讨了乌鲁木齐市气象因素对PM_2.5浓度的影响[5]；2020年，张纯曦等人运用SPSS数据分析和GIS插值方法探讨了乌鲁木齐市PM_2.5浓度的时空分布规律，建立了PM_2.5浓度预测模型[6]；2023年，张一准和颜七笙提出了一种基于最大相关最小冗余算法(Minimum Redundancy Maximum Relevance, MRMR)和麻雀搜索算法(Sparrow Search Algorithm, SSA)优化的BP神经网络模型，并利用该模型预测济南市2019年的PM_2.5浓度[7]；2024年，张晗建立广义可加模型研究了PM_2.5浓度的影响因素并进行了预测[8]。

线性回归模型因其解释性强、应用性广，常被用作基准模型。然而当模型结构发生变化时，其参数会随之发生变化，如果忽视了这些参数的变化继续进行统计分析，可能会导致模型构建不当，从而引发预测误差，使模型失去可靠性，造成不必要的损失，因此进行统计分析与建模之前对模型的变点进行检测和估计尤为重要。分段线性回归模型的关键要素是分割点(即变点)，变点的存在反映了数据中潜在的结构性变化，这些变点是模型中各段之间的“切换点”，分段线性回归模型的最大优势之一是其具有较高的解释性。在传统的线性回归模型中，因变量与自变量之间的关系为单一的线性关系，而在分段线性回归模型中，通过引入多个区间和相应的线性回归模型，可以揭示数据中存在不同的关系，能够更清晰地展示因变量在不同自变量取值范围内的变化趋势。

变点检测与估计是统计学和数据分析中的重要问题，旨在识别数据序列中发生结构变化的时刻，这些变化可能由内部机制的改变或外部冲击引起。近年来，随着数据量的激增和计算能力的提升，变点检测方法得到了广泛研究。但随着高维数据的出现，传统的多变点检测与估计方法表现不佳(少估、错估或过度估计)。因此，研究者开始探索更灵活的多变点检测与估计方法，基于变量选择技术的多变点检测与估计方法是目前十分有效且流行的方法之一。2006年，Davis等人针对相依数据，通过不同自回归过程分割非平稳时间序列，并利用最小描述长度原理(Minimum Description Length, MDL)估计结构变点位置[9]。2013年，Jin等将信号分割后，应用带平滑削边绝对偏离(Smoothly Clipped Absolute Deviation, SCAD)和极大极小凹惩罚(Minimax Concave Penalty, MCP)方法将分段平稳自回归(Piecewise Stationary Autoregressive, PSAR)模型转化为高维线性回归模型后进行变量选择，并检测出了多变点的位置[10]。在后续研究中，2016年，Jin等人又提出了基于自适应Lasso (adaptive-Lasso)、MCP、SCAD等变量选择技术的两阶段多变点检测与估计方法，并对线性回归模型中多变点进行了检测与估计[11]。2020年，周佳琪与金百锁将两阶段多变点检测与估计方法推广至空间网络自回归模型，并分析了合肥市房地产价格[12]。2020年，Sun和Wu研究了广义线性回归模型中的两阶段多变点估计问题[13]。2021年，吕丽和金百锁应用随机加权自助法和高斯混合模型，得到了线性模型中多变点的置信区间[14]。2023年，李美琪等人结合两阶段程序与高维信息准则提出了组正交贪婪算法，并得出具有相依数据结构线性回归模型多结构变点的相合估计[15]。2023年，Ding等对具有分组结构的大型数据集提出组方差膨胀因子(Variance Inflation Factor, VIF)回归模型，结合两阶段程序获得线性回归模型的多变点估计[16]。2024年，廖红怡与金百锁提出多阈值负二项回归模型，通过两阶段迭代程序估计多阈值负二项回归模型的阈值位置和模型参数[17]。

根据《2019年中国生态环境状况公报》，全国超过一半的城市存在大气污染物超标问题，严重或重度污染的天数占全年总天数的78.8%，PM_2.5是空气质量和居民健康的主要影响因素，研究其影响因素并构建更精确的PM_2.5浓度预测模型，对于有效治理大气污染和科学应对雾霾问题具有重要意义。本文以2014年1月1日至2024年9月30日的乌鲁木齐市PM_2.5日平均浓度为研究对象，应用两阶段多变点检测与估计方法对乌鲁木齐市PM_2.5日平均浓度的线性回归模型系数进行多变点检测与估计，建立了更精准的分段线性回归模型。

2. 两阶段多变点检测与估计方法的介绍

检测和估计线性回归模型的多变点时，传统的变点检测方法常常面临着一些挑战，如计算复杂度高、估计精度低(少估、错估或过度估计)等问题。针对线性回归模型的多变点检测问题，为了提高多变点检测与估计的运算速度及估计的精度，2016年，Jin等人提出了一种用于线性回归模型的两阶段多变点检测与估计方法，其中第一阶段切割阶段(the cutting stage)将数据序列切割为$p_n+1$ 个子段，并假设各子段最多包含一个变点(At Most One Change-point, AMOC)，利用特殊的下三角矩阵将线性回归模型的多变点估计问题转为高维线性回归中的变量选择问题，并应用变量选择技术(Lasso, adaptive Lasso, SCAD, MCP)从$p_n+1$ 个子段中选择最有可能包含变点的子段，通过分割数据及应用变量选择技术提高了多变点检测的运算速度；第二阶段精炼阶段(the refining stage)，对切割阶段中筛选出的可能存在变点的子段应用拟似然比检验来确定该子段中是否存在变点，从而提高了多变点估计的精度。该方法能同时准确并有效的估计出线性回归模型中的多变点，且保证了估计结果的相合性。两阶段多变点检测与估计方法的具体步骤如下：

考虑具有s个变点的线性回归模型：

$\begin{array}{c}{y}_i={\displaystyle \sum }_{j=1}^q{x}_{i,j}{\beta }_0+{\displaystyle \sum }_{l=1}^s{\displaystyle \sum }_{j=1}^q{x}_{i,j}{\delta }_{j,l}I\left(a_l<i\le n\right)+{\epsilon }_i,\\ =x_i^{\text{T}}\left[{\beta }_0+{\displaystyle \sum }_{l=1}^s{\delta }_{l}I\left(a_l<i\le n\right)\right]+{\epsilon }_i i=1,2,\cdots ,n,\\ =\{\begin{array}{l}{x}_i^{\text{T}}{\beta }_0+{\epsilon }_{i}1\le i\le a_1,\hfill \\ x_i^{\text{T}}\left({\beta }_0+{\delta }_1\right)+{\epsilon }_i a_1<i\le a_2,\hfill \\ \dots \hfill \\ x_i^{\text{T}}\left({\beta }_0+{\displaystyle \sum }_{l=1}^s{\delta }_l\right)+{\epsilon }_i a_s<i\le n.\hfill \end{array}\end{array}$ (1)

其中，$y_i$ ($i=1,2,\cdots ,n$ )是观测数据序列，$x_i={\left(x_{i,1},\cdots ,x_{i,q}\right)}^{\text{T}}$ 是$q$ 维的解释变量，${\beta }_0={\left({\beta }_{1,0},\cdots ,{\beta }_{q,0}\right)}^{\text{T}}\ne 0$ 是q维回归系数向量，s是变点的个数，$1<a_1<\cdots <a_s<n$ 是s个变点的位置，${\delta }_l={\left({\delta }_{1,l},\cdots ,{\delta }_{q,l}\right)}^{\text{T}}$ ，$l=1,2,\cdots$ ，$s$ 是回归系数在第l个变点处的增量，${\epsilon }_1,\cdots ,{\epsilon }_n$ 表示随机误差，其中${\beta }_0,s,a_1,a_2,\cdots ,a_s,{\delta }_l$ 都是未知的。

2.1. 切割阶段(the Cutting Stage)

2.1.1. 数据切割阶段

将模型(1)的观测数据序列切割成$p_n+1$ 个子段，其中第一个子段的长度为$n-p_{n}m$ ，其余$p_n$ 个子段的长度均为m，且当$n\to \infty$ 时$p_n\to \infty ,m\to \infty$ 。记${\mathcal{S}}_j$ 为第j个子段观测值的指标集的集合，$j=1,2,\cdots ,p_n+1$ 。即第一个子段观测值的指标集为${\mathcal{S}}_1=\left\{1,2,\cdots ,n-p_{n}m\right\}$ ，第j个子段观测值的指标集为${\mathcal{S}}_j=\left\{n-\left(p_n-j+2\right)m+1,\cdots ,n-\left(p_n-j+1\right)m\right\},j=2,\cdots ,p_n+1$ ，其中$m=n/p_n+1$ ($\lfloor ·\rfloor$ 表示向下取整)。

令

$\begin{array}{c}X={\left(X^{\left(1\right)},X^{\left(2\right)},\cdots ,X^{\left(p_n+1\right)}\right)}_{n\times \left(p_n+1\right)q}\\ ={\left(\begin{array}{ccccc}{X}_{\left(1\right)}& 0& 0& \cdots & 0\\ X_{\left(2\right)}& X_{\left(2\right)}& 0& \cdots & 0\\ X_{\left(3\right)}& X_{\left(3\right)}& X_{\left(3\right)}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ X_{\left(p_n+1\right)}& X_{\left(p_n+1\right)}& X_{\left(p_n+1\right)}& \cdots & X_{\left(p_n+1\right)}\end{array}\right)}_{n\times \left(p_n+1\right)q}\end{array}$ (2)

其中，$X^{\left(1\right)}={\left(X_{\left(1\right)}^{\text{T}},X_{\left(2\right)}^{\text{T}},\cdots ,X_{\left(p_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$X^{\left(j\right)}={\left(0_{q\times \left(n-\left(p_n-j+2\right)m\right)},X_{\left(j\right)}^{\text{T}},\cdots ,X_{\left(p_n+1\right)}^{\text{T}}\right)}^{\text{T}}$ ，$X_{\left(1\right)}={\left(x_1,\cdots ,x_{n-p_{n}m}\right)}^{\text{T}}$ ，是第一子段的解释变量，$X_{\left(j\right)}={\left(x_{n-\left(p_n-j+2\right)m+1},\cdots ,x_{n-\left(p_n-j+1\right)m}\right)}^{\text{T}}$ ，$j=2,3,\cdots ,p_n+1$ 是第j子段的解释变量，且$x_i={\left(x_{i,1},\cdots ,x_{i,q}\right)}^{\text{T}},i=1,2,\cdots ,n$ 。这样原来$n\times q$ 维的自变量矩阵被拓展为$n\times \left(p_n+1\right)q$ 维的特殊下三角矩阵。假设每个子段至多有一个变点(AMOC)，若变点$a_j$ 位于第$k_j$ 个子段中，即$a_j\in {\mathcal{S}}_{k_j}=\left\{n-\left(p_n-k_j+2\right)m+1,\cdots ,n-\left(p_n-k_j+1\right)m\right\},j=1,\cdots ,s$ ，记$A=\left\{1,k_1,k_1+1,\cdots ,k_s,k_s+1\right\}$ ，即A表示由第一个子段和变点$a_j$ 所在的子段${\mathcal{S}}_{k_j}$ 及其后一个子段${\mathcal{S}}_{k_j}+1$ 的指标集构成的集合。

通过对数据的切割及特殊的下三角矩阵(2)的构造，模型(1)的矩阵形式为：

$y=X{\theta }^{*}+X_{\omega }\overrightarrow{\omega }+\epsilon ,$ (3)

其中，$y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$ ，$X$ 是由(2)式定义的下三角矩阵，${\theta }^{\text{*}}={\left({\theta }_1^{\text{T}},{\theta }_2^{\text{T}},\cdots ,{\theta }_{p_n+1}^{\text{T}}\right)}^{\text{T}}$ 。(3)式中${\theta }^{*},X_{\omega },\overrightarrow{\omega }$ 的定义以及模型(3)的推导过程详见文献[11]。

2.1.2. 回归系数的自适应Lasso估计

1996年，Tibshiran对于线性回归模型的回归系数提出了Lasso估计[18]。在此基础上，2006年，Zou提出了自适应Lasso估计并讨论了估计量的Oracle性质[19]。由文献[19]知模型(3)中回归系数${\theta }^{*}={\left({\theta }_1^{\text{T}},{\theta }_2^{\text{T}},\cdots ,{\theta }_{p_n+1}^{\text{T}}\right)}^{\text{T}}$ 的自适应Lasso估计为：

${\widehat{\theta }}^{*}=\arg \underset{{\theta }^{*}}{\min }\left\{\frac{1}{n}{\Vert y-X{\theta }^{*}\Vert }^2+\frac{\lambda }{n}{\displaystyle \sum }_{j=1}^{p_n+1}\frac{1}{{\Vert {\tilde{\theta }}_j\Vert }^{\gamma }}{\displaystyle \sum }_{i=1}^q\left|{\theta }_{ji}\right|\right\},$ (4)

其中，${\widehat{\theta }}^{*}=\left({\widehat{\theta }}_1^{\text{T}},{\widehat{\theta }}_2^{\text{T}},\cdots ,{\widehat{\theta }}_{p_n+1}^{\text{T}}\right)$，${\widehat{\theta }}_j={\left({\widehat{\theta }}_{j1},{\widehat{\theta }}_{j2},\cdots ,{\widehat{\theta }}_{jq}\right)}^{\text{T}},{\theta }^{*}={\left({\theta }_1^{\text{T}},{\theta }_2^{\text{T}},\cdots ,{\theta }_{p_n+1}^{\text{T}}\right)}^{\text{T}}$，${\theta }_j={\left({\theta }_{j1},{\theta }_{j2},\cdots ,{\theta }_{jq}\right)}^{\text{T}}$ ，$\lambda$ 是阈值参数，通过BIC信息准则获取，$\gamma >0$ 的常数，${\tilde{\theta }}_j$ 是${\theta }_j$ 的初始值，可以通过最小二乘法得到，即${\tilde{\theta }}_1={\left(X_{\left(1\right)}^{\text{T}}{X}_{\left(1\right)}\right)}^{-1}{X}_{\left(1\right)}^{\text{T}}{y}^{\left(1\right)}$，${\tilde{\theta }}_j={\left(X_{\left(j\right)}^{\text{T}}{X}_{\left(j\right)}\right)}^{-1}{X}_{\left(j\right)}^{\text{T}}{y}^{\left(j\right)}-{\left(X_{\left(j-1\right)}^{\text{T}}{X}_{\left(j-1\right)}\right)}^{-1}{X}_{\left(j-1\right)}^{\text{T}}{y}^{\left(j-1\right)}$，这里$y^{\left(1\right)}={\left(y_1,\cdots ,y_{n-p_{n}m}\right)}^{\text{T}}$ 表示第一个子段的观测值，$y^{\left(j\right)}={\left(y_{n-\left(p_n-j+2\right)m+1},\cdots ,y_{n-\left(p_n-j+1\right)m}\right)}^{\text{T}}$ 表示第$j$ 个子段的观测值，$j=2,\cdots ,p_n+1$ 。

2.2. 精炼阶段(the Refining Stage)

记$\widehat{A}=\left\{j:{\widehat{\theta }}_j\ne 0,j=1,\cdots ,p_n+1\right\}$ 为模型(3)回归系数不为零的子段的指标集的集合。${\widehat{A}}^{*}=\left\{j:j\in \widehat{A},j-1\notin \widehat{A},j=2,\cdots ,p_n+1\right\}=\left\{{\widehat{k}}_1,\cdots ,{\widehat{k}}_{\widehat{s}}\right\}$ 为可能存在变点的子段指标集的集合。为了准确估计变点$a_l$ 的位置${\widehat{a}}_l$ ，应用拟似然比检验方法在可能存在变点的子段及其相邻的上一个子段，即在

${\widehat{\mathcal{S}}}_{\left(l\right)}=\left\{n-\left(p_n-{\widehat{k}}_l+3\right)m+1,\cdots ,n-\left(p_n-{\widehat{k}}_l+1\right)m\right\},l=1,2,\cdots ,\widehat{s}$

中检验是否存在变点，由文献[20]知该变点检验问题等价于以下假设检验问题

$y_i=x_i^{\text{T}}{\beta }_{l}I\left(n_l^{\left(u\right)}\le i\le {\xi }_l\right)+x_i^{\text{T}}{\beta }_{l+1}I\left({\xi }_l<i\le n_l^{\left(v\right)}\right)+{\epsilon }_i$ (5)

$H_0:{\beta }_l={\beta }_{l+1}\leftrightarrow H_1:{\beta }_l\ne {\beta }_{l+1}.$ (6)

且检验问题(6)的拟似然比检验统计量为：

$T_l=N_l\left({\widehat{\sigma }}_l^2-\underset{\beta }{\min }{\displaystyle {\sum }_{i=n_l^{\left(u\right)}}^{{\widehat{\xi }}_l}{\left(y_i-x_i^{\text{T}}\beta \right)}^2}-\underset{\beta }{\min }{\displaystyle {\sum }_{i={\widehat{\xi }}_l+1}^{n_l^{\left(v\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2}\right)/{\widehat{\sigma }}_l^2.$ (7)

其中，$n_l^{\left(u\right)}= n-\left(p_n-{\widehat{k}}_l+3\right)m+1$ ，$n_l^{\left(v\right)}= n-\left(p_n-{\widehat{k}}_l+1\right)m$ 。${\beta }_l$ 与${\beta }_{l+1}$ 是未知的$q$ 维回归系数向量，$i=n_l+1,\cdots ,n_{l+1}$ ，$l=1,2,\cdots ,\widehat{s}$ 。$N_l=n_l^{\left(v\right)}-n_l^{\left(u\right)}$ 为${\widehat{\mathcal{S}}}_{\left(l\right)}$ 的段长，${\widehat{\sigma }}_l^2={\min }_{\beta }{\displaystyle {\sum }_{i=n_l^{\left(u\right)}}^{n_l^{\left(v\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2}$是模型(5)的残差平方和的估计值。进一步，由文献[20]的引理3.1.9知该检验统计量的临界值为$b_l={\left(2\log \log N_l+q\left(\log \log \log N_l\right)/2-\log \Gamma \left(q/2\right)\right)}^2/2\log \log N_l$ ，且当$T_l>b_l+2c_l\text{log}\left(-2/\text{log}\left(1-\alpha \right)\right)$ 时拒绝原假设H₀，认为模型(5)的第$l$ 个子段回归系数存在一个变点${\xi }_l\in \left\{n_l^{\left(u\right)},\cdots ,n_l^{\left(v\right)}\right\}$ ，使得在模型(5)中${\beta }_l\ne {\beta }_{l+1}$ ，其中$c_l={\left(b_l/\left(2\log \log l\right)\right)}^{1/2}$ ，$\alpha$ 为显著性水平。变点${\xi }_l$ 的相合估计量${\widehat{\xi }}_l$ 为：

${\widehat{\xi }}_l=\underset{n_l^{\left(u\right)}+q<k<n_l^{\left(v\right)}-q}{\arg \min }\left[\underset{\beta }{\min }{\displaystyle \sum }_{i=n_l^{\left(u\right)}}^k{\left(y_i-x_i^{\text{T}}\beta \right)}^2+\underset{\beta }{\min }{\displaystyle \sum }_{i=k+1}^{n_l^{\left(v\right)}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2\right].$ (8)

如果$T_l<b_l+2c_l\text{log}\left(-2/\text{log}\left(1-\alpha \right)\right)$ ，接受原假设$H_0$ ，认为第$l$ 个子段不存在变点，此时删除该子段。记由(8)式得到的所有变点的估计量为${\tilde{\xi }}_1<\cdots <{\tilde{\xi }}_{\tilde{s}}$ 。为了提高变点估计的准确性，令${\tilde{n}}_0=0$ ，${\tilde{n}}_l={\tilde{\xi }}_l,l=1,2,\cdots ,\tilde{s},{\tilde{n}}_{\tilde{s}+1}=n$ ，对子段$\left\{{\tilde{n}}_{l-1}+1,\cdots ,{\tilde{n}}_{l+1}\right\}$ 再次应用(7)式进行拟似然比检验，如果存在变点，则最终的变点估计量为：

${\widehat{a}}_l=\underset{{\tilde{n}}_{l-1}+q<k<{\tilde{n}}_{l+1}-q}{\arg \min }\left[\underset{\beta }{\min }{\displaystyle \sum }_{i={\tilde{n}}_{l-1}+1}^k{\left(y_i-x_i^{\text{T}}\beta \right)}^2+\underset{\beta }{\min }{\displaystyle \sum }_{i=k+1}^{{\tilde{n}}_{l+1}}{\left(y_i-x_i^{\text{T}}\beta \right)}^2\right].$ (9)

变点估计的相合性及相关引理、定理和其更详细内容请参见文献[11]及其补充材料。

3. 实证分析

3.1. 研究区域概括

本文研究区为乌鲁木齐市，位于新疆维吾尔自治区，是该省首府。地处我国西北内陆干旱区的新疆中北部，天山北麓和准格尔盆地南部。市区南北宽度约153公里，东西最长约190公里，海拔高度在580米至920米之间。乌鲁木齐市的山区面积占总面积的50%以上，南部和东北部地势较高，中部和北部则较低，气候类型为中温带大陆性干旱气候，属于典型的“河谷型城市”，年降水量稀少。由于地理位置及扬尘、雾霾等因素的影响，大气污染物难以扩散容易导致严重的空气污染。

3.2. 变量选取与数据来源

根据以往的研究表明影响PM_2.5浓度的因素大致分为两类：一类是空气污染物，主要包括可吸入颗粒物(PM₁₀)、一氧化碳(CO)、二氧化硫(SO₂)、二氧化氮(NO₂)和臭氧8小时滑动平均浓度(O₃)等微观变量；另一类是气象因素，包括气温、气压、湿度、风速和降水量等宏观变量；此外，还有其他宏观变量，如人口密度、机动车尾气排放量、GDP增长率、植被覆盖率等。因此，本文使用微观变量与宏观变量对乌鲁木齐市PM_2.5日平均浓度进行建模。考虑到上述变量数据的可获得性与公开性，本文选取了影响乌鲁木齐市PM_2.5日均浓度的因素为空气污染物¹和气象因素²，如表1所示。以2014年1月1日至2024年9月30日的乌鲁木齐市PM_2.5日平均浓度为研究对象(共3924个数据)，其中2014年1月1日至2024年7月1日的数据(共3833个数据)作为训练集，2024年7月2日至2024年9月30日的数据(共91个数据)作为测试集，用作评估模型的短期预测能力。

Table 1. PM_2.5 and the pollutants and meteorological factors affecting PM_2.5 concentration and their units

表1. PM_2.5及影响PM_2.5浓度的污染物和气象因素与其单位

空气污染物(微观变量)			气象因素(宏观变量)
变量名称	代表符号	单位	变量名称	代表符号	单位	备注说明
细微颗粒	PM2.5	μg/m^3	平均气温	T	˚C	地面以上两米处大气温度
可吸入颗粒	PM10	μg/m^3	平均气压	P	mmHg	平均海平面的大气压
一氧化碳	CO	mg/m^3	平均相 对湿度	U	%	地面以上两米处相对湿度
二氧化硫	SO2	μg/m^3
二氧化氮	NO2	μg/m^3	平均风速	Ff	m/s	地面以上10米处平均风速
臭氧	O3	μg/m^3	总降水量	Rain	mm

3.3. 数据预处理

由于数据来源不统一、存在缺失值、量纲不一致、数值范围差异大等现象，在分析建模之前对数据进行预处理。本文采用以下方法：

1) 逐时数据转化为逐日数据：由于获取的气象因素为逐时数据(每日02:00至23:00每3小时采集一次的逐时数据)，空气污染物为逐日数据，因此需要把气象逐时数据转换为日数据，使得气象逐时数据与污染物逐日数据相匹配。本文将一天中在不同时段的平均气温、相对湿度、风速和气压等变量的平均值作为相应变量的日平均数据。

2) 缺失值处理：由于监测站设备故障或维护，部分污染物的浓度数据可能会缺失。本文采用均值法填补缺失值，即应用缺失值前后数据的算术平均值进行填补，以保证数据的完整性。

3) 数据标准化：由于自变量的量纲不同且变量之间数值范围差异较大，这可能导致模型训练时出现偏差，因此对数据进行标准化处理，以消除量纲差异，有助于提高模型的收敛速度，提升模型的性能和预测准确性。

3.4. 自变量的筛选

为了模型的稳定性及预测的准确性，在建立模型之前对自变量进行筛选。以PM_2.5的日平均浓度为因变量(y)，以PM₁₀ (x₁)、CO (x₂)、SO₂ (x₃)、NO₂ (x₄)、O₃ (x₅)的日平均浓度及日平均气温(x₆)、气压(x₇)、相对湿度(x₈)、风速(x₉)、总降水量(x₁₀)为自变量，应用R软件通过逐步回归对自变量进行筛选，结果如表2所示。

Table 2. Selection of significant variables in the training set using stepwise regression method

表2. 通过逐步回归法选取训练集中的显著变量

变量	系数(Estimate)	标准误差(Std. Error)	t值(t value)	p值(Pr (>|t|))	显著性
截距(Intercept)	2.000e−15	4.623e−03	0	1
x1 (PM10)	4.181e−01	7.826e−03	53.43	<2e−16	***
x2 (CO)	5.768e−01	1.228e−02	46.963	<2e−16	***
x3 (SO2)	−1.31e−01	6.441e−03	−20.328	<2e−16	***
x4 (NO2)	6.746e−02	1.014e−02	6.655	3.23e−11	***
x5 (O3)	1.965e−01	8.899e−03	22.086	<2e−16	***
x6 (平均气温)	−1.78e−01	1.033e−02	−17.21	<2e−16	***
x7 (平均气压)	1.112e−04	5.468e−03	2.033e−02	9.838e−01
x8 (平均相对湿度)	9.633e−02	9.125e−03	10.556	<2e−16	***
x9 (平均速度)	6.01e−03	4.962e−03	1.2131	2.251e−01
x10 (总降水量)	−4.470e−03	5.149e−03	−8.680e−01	3.854e−01

注：***表示显著性水平等于0.001的情况下p < 0.001，说明该变量对模型具有显著影响。

从逐步回归的结果可知，主要影响乌鲁木齐市PM_2.5日平均浓度(y)的污染物有PM₁₀ (x₁)、CO (x₂)、SO₂ (x₃)、NO₂ (x₄)、O₃ (x₅)；气象因素有日平均气温(x₆)和日平均相对湿度(x₈)，采用选取的显著变量建立的全局线性回归模型如下：

$\begin{array}{l}y=0.4181x_{1t}+0.5768x_{2t}-0.1309x_{3t}+0.06746x_{4t}+0.1965x_{5t}\\ -0.1778x_{6t}+0.09633x_{8t}1\le t\le 3832\end{array}$ (10)

接下来应用两阶段自适应Lasso多变点检测与估计方法对模型(10)的回归系数进行多变点检测和估计，得出的变点位置分别为1065、1588、2142、2944，对应时间分别为2016年12月2日、2018年5月09日、2020年1月30日、2022年1月24日，变点位置如下图1红色虚线所示。

Figure 1. Location of multiple change-point in the coefficients of the linear regression model for daily average PM_2.5 concentration in Urumqi City

图1. 乌鲁木齐市PM_2.5日平均浓度线性回归模型系数的多变点位置

从气象角度来看出现上述变点的原因是2016年冬季受大气环流异常影响，冷空气势力弱，暖气团控频繁且时间较长导致逆温现象、强度比往年强，且底层水汽充沛，使得雾日明显偏多，污染物更易在低空积聚不利于扩散，增加了大气中PM_2.5浓度，此外当时乌鲁木齐市能源结构中煤炭占比较大，冬季采暖与工业燃煤排放大量的烟尘、PM₁₀、CO、SO₂、NO₂、O₃等污染物加剧了PM_2.5浓度的上升；2017年起乌鲁木齐市逐步开始改造小燃煤供热设施，煤炭燃烧量减少，致使PM₁₀、CO、SO₂、NO₂等污染物排放量显著降低，直接削减了PM_2.5的生成源，直至2018初显成效，又由于冬季结束，停止供暖，燃煤排放量进一步减少，气温回升使大气垂直扩散条件改善，有利于污染物向上输送稀释，打破了冬季污染物积聚态势，相对湿度也处于较适宜的范围，降低了PM_2.5浓度上升的风险，另外O₃浓度因光照增强，不再助推PM_2.5浓度增加，综合以上因素2018年5月PM_2.5浓度相比之前有所下降；疫情突发后，在2020年1月26日起多数工业企业停工停产导致PM₁₀、CO、SO₂、NO₂等污染物的浓度下降，扬尘等粗颗粒物排放归零，切断了PM_2.5生成的工业渠道，此外民众居家、机动车出行锐减，尾气排放几乎停止，减少了NO₂、CO及碳氢化合物等尾气成分，抑制了生成PM_2.5的路径。气温与相对湿度虽处于冬季常态范围，但因污染源骤减，PM_2.5浓度随之显著下降；2021年7月1日，乌鲁木齐市实施由乌鲁木齐市人民代表大会常务委员会公布修改的《乌鲁木齐市大气污染防治条例》，条例明确推进能源结构优化、扩大清洁能源占比、削减工业燃煤、大力淘汰老旧柴油货车、管控机动车尾气等一系列政策，从多方面对传统污染源进行有效压制，从源头降低了PM₁₀、CO、SO₂的排放量，与此同时乌鲁木齐市入选冬季清洁取暖项目，天然气和电力等清洁能源逐步替代传统煤炭供暖，通过以上措施最终促成2022年初乌鲁木齐市PM_2.5浓度明显下降，空气质量显著改善。

上述4个变点将整个数据分成五个子段(第一子段[1~1064]、第二子段[1065~1587]、第三子段[1588~2141]、第四子段[2142~2943]、第五子段[2944~3832])，对每个子段的数据应用R软件得到乌鲁木齐市PM_2.5日平均浓度(y)的分段线性回归模型，如下：

$y=\{\begin{array}{l}0.4433x_{1t}+0.4492x_{2t}-0.03339x_{3t}+0.0506x_{4t}+0.0555x_{5t}\\ +0.0123x_{6t}+0.1494x_{8t}-0.16651\le t\le 1064，\\ 0.5193x_{1t}+0.5105x_{2t}+0.2272x_{3t}+0.00634x_{4t}+0.1860x_{5t}\\ -0.0676x_{6t}+0.2604x_{8t}-0.07751065\le t<1588，\\ 0.4245x_{1t}+0.6953x_{2t}+0.1776x_{3t}-0.0860x_{4t}+0.1604x_{5t}\\ -0.0995x_{6t}+0.1113x_{8t}-0.01631588\le t<2141，\\ 0.2261x_{1t}+1.0800x_{2t}-0.2657x_{3t}+0.0304x_{4t}+0.1803x_{5t}\\ -0.0122x_{6t}+0.0856x_{8t}-0.04882142\le t<2943，\\ 0.3206x_{1t}+1.2016x_{2t}-{0.1724}_{3t}-0.1216x_{4t}+0.0672x_{5t}\\ -0.0339x_{6t}+0.0180x_{8t}+0.29002944\le t\le 3832，\end{array}$ (11)

表3列出了模型(10)与模型(11)的均方误差(MSE)、均方根误差(RMSE)，以及修正后的拟合优度R²。

Table 3. Stability statistical indicators of model (10) and model (11)

表3. 模型(10)与模型(11)的稳定性统计指标

模型	均方误差(MSE)	均方根误差(RMSE)	Adjusted R-squared
全局线性回归模型(10)	0.0817	0.2858	0.9181
分段线性回归模型(11)的第一子段	0.0797	0.2822	0.9311
分段线性回归模型(11)的第二子段	0.0639	0.2527	0.9606
分段线性回归模型(11)的第三子段	0.0384	0.1960	0.9383
分段线性回归模型(11)的第四子段	0.0511	0.2261	0.9317
分段线性回归模型(11)的第五子段	0.0320	0.1790	0.9484

图2给出了由模型(10)与模型(11)得到的乌鲁木齐市PM_2.5日平均浓度的拟合值与真实值的对比图，其中蓝色实线表示真实值，绿色虚线和红色虚线分别表示由模型(10)与模型(11)得到的拟合值。

Figure 2. Comparison chart of fitted values and actual values obtained from model (10) and model (11)

图2. 由模型(10)与模型(11)得到的拟合值与真实值的对比图

从表1和图2可以看出，模型(11)在均方误差、均方根误差和拟合优度等方面都均优于模型(10)，拟合值更贴近真实值，能够更准确地反映乌鲁木齐市PM_2.5日平均浓度的实际情况。

现将模型(10)和模型(11)的第五子段应用于测试集(2024年7月2日至2024年9月30日，共91天的乌鲁木齐市PM_2.5日平均浓度)上，进行短期预测，以评估模型的短期预测能力。图3展示了真实值与短期预测结果的对比图，其中蓝色实心圆线表示真实值，绿色实心三角形线和红色空心圆线分别表示由模型(10)和模型(11)的第五子段得到的短期预测结果。

Figure 3. Comparison chart of predicted results and actual values obtained from model (10) and the fifth subsection of model (11)

图3. 由模型(10)与模型(11)的第五子段得到的预测结果与真实值的对比图

表4列出了由模型(10)与模型(11)的第五子段得到的短期预测结果与真实值的残差均值、残差标准差、平均绝对误差(MAE)及平均绝对百分比误差(MAPE)。

Table 4. Statistical indicators of short-term prediction results obtained from model (10) and the fifth subsection of model (11)

表4. 由模型(10)与模型(11)的第五子段得出的短期预测结果的统计指标

统计量	全局线性回归模型(10)	分段线性回归模型(11)的第五子段
残差均值	6.4468	0.7862
残差标准差	4.8369	4.7227
平均绝对误差(MAE)	6.8862	3.8177
平均绝对百分比误差(MAPE)	62.8519	31.7740

图4给出了模型(10)与模型(11)的第五子段得到的短期预测结果与真实值的残差图。

Figure 4. Comparison chart of residuals for short-term prediction values from model (10) and the fifth subsection of model (11)

图4. 由模型(10)与模型(11)的第五子段得到的短期预测值的残差对比图

表5给出了基于全局线性回归模型(10)与分段线性回归模型(11)中第五子段得到的短期预测残差的方差和均值显著性检验结果。

Table 5. Variance and mean significance test of prediction residuals obtained from model (10) and the fifth subsection of model (11)

表5. 由模型(10)与模型(11)的第五子段得出的预测残差的方差和均值显著性检验

模型	F检验		T检验
	F检验统计量	P值	T检验统计量	P值
全局线性回归模型(10)	1.0489	0.4106	12.715	2.2e−16
分段线性回归模型(11)的第五子段			1.588	0.1158

从表5中的F检验可知，在显著性水平$\alpha =0.05$ 的情况下，全局线性回归模型(10)与分段线性回归模型(11)的第五子段在预测残差的方差上没有显著性差异；随后对预测残差的均值是否为零进行显著性检验，由T检验的结果可知，在显著性水平$\alpha =0.05$ 的情况下，对于全局线性回归模型(10)预测残差均值的T检验，其p值为2.2e−16 < 0.05，表明模型(10)的预测残差均值不为零，说明模型(10)存在系统性偏差；而分段线性回归模型(11)的第五子段的T检验p值为0.1158 > 0.05，表明模型(11)第五子段预测残差的均值为零，说明模型(11)的第五子段没有系统性偏差，模型具有稳定性。从表4和图4可以看出，由分段线性回归模型(11)的第五子段得到的预测值与真实值比较接近，且残差较小，说明能够更准确地预测乌鲁木齐市PM_2.5的日平均浓度。这表明根据变点位置建立的分段线性回归模型比全局线性回归模型在模型稳定性及预测精度上表现更优，说明应用线性回归模型研究PM_2.5浓度的变化及走势前检测与估计回归系数的变点是至关重要的。

4. 结论

2024年3月18日，生态环境部公布2月全国168个重点城市空气质量排名，乌鲁木齐市处于后20位，这意味着乌鲁木齐城市发展速度加快的同时，其空气质量也随之恶化，为了客观分析乌鲁木齐市的环境污染问题，并采取有效的治理措施，对乌鲁木齐市大气环境的相关研究是非常有必要的。本文选取2014年1月1日至2024年9月30日的乌鲁木齐市PM_2.5日平均浓度数据作为研究对象，以提高PM_2.5浓度线性回归模型的拟合优度及预测能力为出发点对其进行研究。采用两阶段自适应Lasso多变点检测与估计方法对乌鲁木齐市PM_2.5日平均浓度线性回归模型的系数进行多变点检测与估计，从而可以根据变点位置建立更精准的分段线性回归模型来刻画乌鲁木齐市PM_2.5日平均浓度的情况，进一步说明进行统计分析与建模之前进行变点检测与估计是至关重要的，对提升模型稳定性和预测精度具有显著作用。

基金项目

新疆维吾尔自治区自然科学基金项目(2023D01A37)。

NOTES

^*通讯作者。

¹空气污染物数据来源于中国空气质量在线监测分析平台：https://www.aqistudy.cn/。

²气象数据来源于“Reliable Prognosis”(RP5)网站：https://rp5.ru/。气象数据为每日02:00至23:00每3小时采集一次的逐时数据。