基于经验模态分解的大气能见度预测

Atmospheric Visibility Prediction Based on Empirical Mode Decomposition

摘要: 经验模态分解(EMD)可以处理非线性、非平稳的时间序列数据，使之成为一组便于提取内部特征的分量。本文提出了一种结合经验模态分解(EMD)和AR模型的大气能见度预测方法。利用EMD将大气能见度数据分解为一组便于提取内部特征的本征模态函数(IMFs)。分别计算IMFs与原数据的相关系数，剔除相关性弱的本征模态函数以达到去噪目的。使用AR模型分别预测各个IMF，将各个预测值相加，得到最终的大气能见度预测值。为了验证所使用模型的预测性能，将只用AR模型的预测方法作为对比模型，并对比了两种模型的均方误差(RMSE)等指标，结果表明，结合了EMD的AR模型的预测性能优于只用AR模型的预测方法。

Abstract: Empirical Modal Decomposition (EMD) can deal with nonlinear, nonsmooth time series data into a set of components that facilitate the extraction of internal features. In this paper, we propose an atmospheric visibility prediction method that combines empirical modal decomposition (EMD) and AR modeling. The atmospheric visibility data are decomposed into a set of intrinsic modal functions (IMFs) that facilitate the extraction of internal features using EMD. The correlation coefficients between the IMFs and the original data were calculated separately, and the weakly correlated intrinsic modal functions were eliminated for denoising purposes. Each IMF was predicted separately using the AR model, and the predicted values were summed to obtain the final atmospheric visibility prediction. In order to verify the prediction performance of the model used, the prediction method using only the AR model was used as a comparison model, and the mean square error (RMSE) and other indexes of the two models were compared, and the results showed that the prediction performance of the AR model combined with the EMD was better than that of the prediction method using only the AR model.

1. 引言

大气能见度的预测对于保障交通运输安全、环境监测与保护、城市规划与建设、公共健康以及气象预报与服务等方面都具有重要的意义。低能见度条件下，如雾、霾等天气，会严重降低驾驶员的视线范围，增加交通事故的风险。因此，提前预测大气能见度可以帮助交通管理部门及时采取应对措施，如限速、封路等，从而保障交通安全。

目前，大气能见度的预测方法可以分为两类：一类是基于天气研究和预报(Weather Research and Forecasting Model, WRF)模型的大气能见度预测。Huang等人[1]使用WRF模型预测了乌鲁木齐市气象台的能见度。Dimitrova等人[2]为了探究大气在陆地和海面分布的差别，使用WRF模型分别预测了加拿大大西洋沿岸海面上和陆地上的能见度，发现WRF模型预测海面上能见度的性能优于陆地上的性能，从而得出海洋表面的大气覆盖得更加均匀的结论。Zhang等人[3]采用改进后的WRF模型，即高时空分辨率的WRF-Chem (WRF模型与化学相结合)模型预测了中国中东部地区(具体位置为东经112~122度，北纬34~42度)的能见度，预测效果较好。尽管WRF模型广泛应用于大气能见度预测中，并取得了较好的预测效果，但文献[4]认为WRF模型通常基于已知的天气规律进行建模，对不同地区季节的参数模拟预测效果差异较大。

另一类是基于多种机器学习方法的大气能见度预测，通过训练算法和模型优化，机器学习方法可以从数据中学习到更准确的模型和规律。Wang等人[5]使用了主成分分析和深度信念网络(Deep Belief Network, DBN)来预测短期和长期序列中的大气能见度，预测结果优于BP神经网络和卷积神经网络。Zang等人[6]使用了基于跳帧传输门(Frame-Hopping Transmission Gate, FHTG)、特征融合模块(Feature Fusion Module, FFM)和反向调度采样(Reverse Scheduled Sampling, RSS)的递归神经网络(Recursive Neural Network, RNN)预测模型SwiftRNN，对我国中东部地区的大气能见度进行了预测，预测效果优于ConvLSTM和PredRNN模型。Chen等人[7]提出了用于短期能见度预测的图卷积网络和门控循环单元GCN-GRU (Graph Convolutional Network and Gated Recurrent Unit)模型，预测了2017~2018年江苏省气象站和环境站的监测数据，预测效果优于独立的GRU模型。Yang等人[8]使用了基于BP神经网络[9]和光梯度增强机(Light Gradient Boosting Machine, LGBM)分类器的能见度预测模型——BP-LGBM组合方法(BP-LGBM Combination Method, BLCM)，预测了美国气象站2010年至2013年的大气能见度数据，预测效果优于BP神经网络。Baran等人[10]使用了比例优势逻辑回归(Proportional Odds Logistic Regression, POLR)模型，预测了欧洲中期预测中心(European Centre for Medium-Range Forecasts, ECMWF)的能见度数据，预测效果优于多层感知器(Multilayer Perceptron, MLP)模型。

2023年，Li等人[11]利用经验模态分解(Empirical Mode Decomposition, EMD)处理了大气相干长度数据，并结合了基于LSTM网络的编码器–解码器模型对处理后的大气相干长度数据进行预测，得到了较好的大气相干长度预测效果。EMD的主要思想是使用希尔伯特–黄变换(Hilbert-Huang Transform, HHT)来分解非线性和非平稳数据集，直到最终数据集平稳。EMD方法也适用于大气能见度数据。可以将其分解为平稳的分量数据，便于分别学习各分量的特征来进行预测工作。

综上，本文基于EMD，建立了一种新的大气能见度预测模型。先利用EMD将数据分解为平稳线性的本征模态函数(Intrinsic Mode Functions, IMFs)，针对数据含有噪声的问题，计算各个IMF与原数据集的相关系数，通过去掉相关性弱的IMF以达到去噪效果。再使用AR自回归模型分别对各个IMF进行预测，最后将各预测结果进行相加得出最终预测结果。选用深圳市气象局网站的数据进行实证分析，对深圳市大浪街道的大气能见度进行预测，通过均方误差等指标得出，该模型具有较好的预测效果。

2. 大气能见度的预测模型

记t时刻的大气能见度数据为$y\left(t\right)$ 。首先对能见度观测数据使用经验模态分解进行处理，求出分解后的各本征模态函数与原数据的相关系数，将相关系数低的本征模态函数视为噪声并剔除，最后用AR模型预测各个本征模态函数，将得到的预测结果相加，得到最终的预测结果。

2.1. 大气能见度的经验模态分解

通过EMD将$y\left(t\right)$ 分解为本征模态函数IMFs与残差的和：

$y\left(t\right)={\displaystyle \sum _{i=1}^{n}IM{F}_i\left(t\right)}+r\left(t\right)$ (1)

其中$IM{F}_i\left(t\right)$ 表示分解后的本征模态函数，n表示IMF的数量，$r\left(t\right)$ 表示残差。

大气能见度数据的经验模态分解处理过程如下：

在原始大气能见度数据中，找到所有的极大值点和极小值点。根据找到的极值点，构建上包络线和下包络线。上包络线连接所有的极大值点，下包络线连接所有的极小值点。将上包络线和下包络线的平均值计算出来。将原始数据减去均值，得到去趋势后的数据。如果去趋势后的数据满足整个数据序列中极值的数目(极大值和极小值的数目总和)和零点的数目相等或相差不超过1，而且在任意位置，局部最大包络线与局部最小包络线的平均值为0，就将去趋势后的数据作为第一个IMF。如果不满足条件，将去趋势后的数据作为新的原始数据，重复以上步骤，直到满足条件为止。每一次迭代都会产生一个新的IMF。将所有的IMF相加，得到一个近似原始数据的重构序列。将原始数据与重构序列之间的差值视为残差。如果残差足够小，则认为分解结果已经足够精确。参照文献[12]中Huang等人对EMD的残差研究，如果残差仍然较大，则将残差视为新的原始数据，重复以上步骤，直到满足

$\text{0}\text{.2}\le {\displaystyle \sum _{t=1}^N\frac{{\left[IM{F}_{i-1}\left(t\right)-IM{F}_i\left(t\right)\right]}^2}{{\left[IM{F}_{i-1}\left(t\right)\right]}^2}}\le 0.3$ (2)

为止，其中N表示序列长度。

2.2. 大气能见度的去噪

在观测到的大气能见度数据集中存在噪声，因此，在预测之前，要先对大气能见度数据集进行去噪。将分解出的高频分量视为噪声并直接剔除是常见的做法。这种方法常常会丢失有用信息。为了减少信息损失，计算每个IMF与观测到的大气能见度数据集之间的相关系数

$\rho =\frac{{\displaystyle \sum _{t=1}^n\left[y\left(t\right)-\overline{y}\right]\left[IM{F}_i\left(t\right)-\overline{IM{F}_i}\right]}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left[IM{F}_i{\left(t\right)}^2-\overline{IM{F}_i}\right]}^2}}\sqrt{{\displaystyle \sum _{i=1}^n{\left[y{\left(t\right)}^2-\overline{y}\right]}^2}}}$ (3)

去除相关系数小于0.2的弱相关分量[13]以达到去噪效果。

2.3. AR模型预测大气能见度

假设能见度过去值$IM{F}_i\left(t-1\right),IM{F}_i\left(t-2\right),\cdots ,IM{F}_i\left(t-p\right)$ 与t时刻IMF值$IM{F}_i\left(t\right)$ 为一线性关系。基于AR自回归模型的t时刻IMF值$IM{F}_i\left(t\right)$ ($i=1,\cdots ,n$ ，n为IMF的数量)的表达式如下：

$IM{F}_i\left(t\right)=a_{i1}IM{F}_i\left(t-1\right)+a_{i2}IM{F}_i\left(t-2\right)+\cdots +a_{i{p}_i}IM{F}_i\left(t-p_i\right)+{\epsilon }_{it},t\in \mathbb{Z}$ (4)

其中$p_i$ 表示AR模型的阶数，$a_i={\left(a_{i1},a_{i2},\cdots ,a_{i{p}_i}\right)}^{\text{T}}$ 表示模型的自回归系数，$a_{i1},a_{i2},\cdots ,a_{i{p}_i}\left(a_{i{p}_i}\ne 0\right)$ 为实数，$\left\{{\epsilon }_{it}\right\}$ 表示白噪声序列$WN\left(0,{\sigma }_i^2\right)$ 。

由最大似然估计，得自回归系数$a_i$ 的估计量为

${\widehat{a}}_i=\left({\widehat{a}}_{i1},{\widehat{a}}_{i2},\cdots ,{\widehat{a}}_{i{p}_i}\right)={\left(X_i^{\text{T}}{X}_i\right)}^{-1}{X}_i^{\text{T}}{Y}_i$ (5)

方差${\sigma }_i^2$ 的最大似然估计为

${\widehat{\sigma }}_i^2=\frac{1}{N-p_i}{\left(Y_i-X_i{\widehat{a}}_i\right)}^{\text{T}}\left(Y_i-X_i{\widehat{a}}_i\right)$ (6)

其中$Y_i=\left[\begin{array}{c}IM{F}_i\left(p_i+1\right)\\ IM{F}_i\left(p_i+2\right)\\ ⋮\\ IM{F}_i\left(N\right)\end{array}\right]$ ，$X_i=\left[\begin{array}{cccc}IM{F}_i\left(p_i\right)& IM{F}_i\left(p_i-1\right)& \cdots & IM{F}_i\left(1\right)\\ IM{F}_i\left(p_i+1\right)& IM{F}_i\left(p_i\right)& \cdots & IM{F}_i\left(2\right)\\ ⋮& ⋮& & ⋮\\ IM{F}_i\left(N-1\right)& IM{F}_i\left(N-2\right)& \cdots & IM{F}_i\left(N-p_i\right)\end{array}\right]$ 。

由于每个IMF具有不同的数据特征，因此需要为每个IMF独立选择合适的AR模型阶数。这里使用贝叶斯信息准则(Bayesian Information Criterion, BIC) [14]，其定义如下：

$BIC\left(p_i\right)=\ln {\widehat{\sigma }}_i^2+\frac{\ln N\left(p_i+1\right)}{N}$ (7)

BIC准则确定的模型阶数是其真值的一致估计。最终得到的$AR\left(p_i\right)$ 模型的预测值为

$\widehat{IM{F}_i}\left(t\right)={\widehat{a}}_{i1}IM{F}_i\left(t-1\right)+{\widehat{a}}_{i2}IM{F}_i\left(t-2\right)+\cdots +{\widehat{a}}_{i{p}_i}IM{F}_i\left(t-p_i\right)+{\epsilon }_{it},t\in \mathbb{Z}$ (8)

最后，将各$\widehat{IM{F}_i}\left(t\right)$ 连同残差$r\left(t\right)$ 相加，得到如下式所示的大气能见度数据$y\left(t\right)$ 的预测结果$\widehat{y}\left(t\right)$ ：

$\widehat{y}\left(t\right)={\displaystyle \sum _{i=1}^n\widehat{IM{F}_i}\left(t\right)}+r\left(t\right)$

3. 实证分析

选取深圳市气象台官网上2024年1月1日至2024年3月15日深圳市大浪街道的大气能见度数据进行实证分析，该数据集共1798个。将本文所提出的基于经验模态分解的AR模型预测方法(简记为EMD_AR模型)与只用AR模型的预测方法(简记为AR模型)进行比较，对比平均绝对误差(MAE)、平均绝对百分比误差(MAPE)、均方根误差(RMSE)和决定系数(R²)等指标以比较预测值与真实值的差异以评估EMD_AR模型的预测性能。

用EMD将大气能见度数据分解成IMFs，计算各IMF与原数据的相关系数，将弱相关的IMF剔除，接着将剩余各IMF分为两部分：用于训练的训练集(70%，1259个数据)和用于验证模型预测效果的测试集(30%，539个数据)，分别使用AR自回归模型对各IMF进行预测，最后将各IMF的预测结果连同残差$r\left(t\right)$ 相加得出最终的预测结果。

MAE、MAPE和RMSE是时间序列预测中比较标准的评价指标，它们用于评估预测结果与原数据的差距，越小越好。R²用于评估模型的拟合效果，值越接近1，模型的拟合效果越好。通过计算这些指标，可以验证所用预测方法的预测效果。

使用EMD算法分解观测到的大气能见度数据集，当本征模态函数IMFs满足(2)式条件时，将观测到的大气能见度数据集分解为如图1所示的一系列具有信息成分和噪声成分的8个本征模态函数IMFs。

Figure 1. Empirical modal decomposition results

图1. 经验模态分解结果

通过求各个IMF与原数据的相关系数，将那些相关系数低于0.2的IMFs视为噪声成分并予以剔除，以达到数据去噪的目的。每个IMF和原数据的相关系数如表1所示。

Table 1. Correlation coefficients for each IMF and original data

表1. 各IMF和原数据的相关系数

IMFs	相关系数
IMF1	0.2828
IMF2	0.4286
IMF3	0.4493
IMF4	0.4705
IMF5	0.3346
IMF6	0.3321
IMF7	0.2489
IMF8	0.2549

根据得出的相关系数结果，保留所有本征模态函数IMF1~IMF8。

在将数据分解和去噪后，对每个剩下的IMF使用AR自回归模型进行预测，最后对每个预测的IMF进行求和，得到最终预测结果。AR模型的阶数通过BIC信息准则进行确定。

以IMF1的预测过程为例，可通过绘制IMF1的BIC热度图来选择自回归模型的阶数$p_{\text{1}}$ 。

Figure 2. BIC heat map for IMF1

图2. IMF1的BIC热度图

BIC值越小，表示模型在拟合数据的同时，复杂度也较低。由图2，可知BIC最小值为3.713，对应AR模型阶数为4，于是，选用AR(4)模型对IMF1进行预测。最终，得出预测结果如图3。

Figure 3. AR(4) model predicts IMF1 results

图3. AR(4)模型预测IMF1结果

同理，通过BIC准则，分别选用AR(8)模型、AR(7)模型、AR(8)模型、AR(2)模型、AR(1)模型、AR(1)模型、AR(1)模型得出IMF2~IMF8的预测结果，根据公式(8)，将各IMF预测结果连同残差$r\left(t\right)$ 相加，得出最终预测结果。对照实验选用不经过EMD数据处理的AR自回归模型，预测结果对比如图4所示。

Figure 4. Prediction results of both EMD_AR and AR models

图4. EMD_AR和AR两种模型预测结果

Table 2. Correlation coefficients for each IMF and original data

表2. 各IMF和原数据的相关系数

模型	R^2	RMSE	MAE	MAPE
EMD_AR	0.93104	1.12030	0.76942	0.11808
AR	0.82207	1.78790	1.21640	0.20394

分别计算EMD_AR和AR两种模型预测值和测试集数据的决定系数(R²)、均方根误差(RMSE)、平均绝对误差(MAE)和平均绝对百分比误差(MAPE) 4个指标，得出了如表2所示的结果。根据表2中不同指标的对比可以得出：EMD_AR模型的R²值更大且更接近于1，EMD_AR模型的其它指标均小于AR模型对应的其它指标。例如：EMD_AR模型和AR模型的RMSE分别为1.1203和1.7879。由此可见，EMD_AR模型预测大气能见度数据的效果比AR模型更好。

4. 结论

本文提出了一种结合经验模态分解(EMD)和AR自回归模型的大气能见度预测方法——EMD_AR模型，用深圳市气象台的大气能见度数据进行了验证。利用EMD将大气湍流能见度数据集分解为一组称为本征模态函数(IMFs)的平稳分量，计算各个IMF与原数据集的相关系数，通过去掉相关性弱的IMF以达到去噪效果。AR模型的阶数可以通过BIC准则进行选择，BIC准则有助于在大气能见度样本量较大时选择拟合更加精确的AR模型，从而避免过拟合。为了研究所使用模型的性能，对比了只用AR模型对大气能见度数据的预测实验，所使用的EMD_AR模型在均方误差(RMSE)等评估标准上均明显优于只用AR模型的预测方法。由于AR模型对复杂关系的处理能力有限，当时间序列中存在复杂的非线性关系时，AR模型可能无法充分描述这些关系，从而导致预测效果下降，其预测精度不如ARIMA、支持向量机等模型。探索EMD和ARIMA、支持向量机等其他模型结合对大气能见度的预测效果是下一阶段的研究目标。

基金项目

国家自然科学基金面上项目(12271271)。

NOTES

^*通讯作者。