基于广义高斯过程回归的北京市近五十年能见度分析

李海蓉，朱晓欣，曹春正

（南京信息工程大学 数学与统计学院，南京 210044）

0 引言

大气能见度的降低逐渐成为气候和空气污染研究领域的关键问题，引起了政府部门和公众的高度重视。目前，大气能见度的研究主要集中在影响因素分析、趋势研究和预测推断上。需要指出的是，虽然大气能见度明确地受到气象条件和污染物状况的影响，但这种复杂的影响关系很难通过机理分析或常规的统计分析方法获取。Lin等[1]研发了一种研究能见度、空气质量及气象条件间因果关系的最优经验回归模型；Yang等[2]通过PLAM/h指标建立了空间大尺度上低能见度的对数线性模型。此类方法的明显缺陷是忽视了数据的序列相关性以及变量之间的非线性相关性。各种时间序列模型可以一定程度上刻画这些特征。例如，Pedrycz等[3]运用模糊时间序列模型对上海的大气能见度进行了预测；Yang等[4]研究了雾霾变化的影响因子及作用机制，并运用时间序列分析提出一种对北京雾霾的长期预测模型。另外一种处理时空数据的常用方法是通过预报方程进行数值模拟，比如CALPUFF模式[5,6]。但是数值预报方法存在几个问题：首先，预报方程的建立需要机理支撑，这对于很多情况难以实现；其次，由于预报方程的动态性和复杂性，当数据信息不够充分时，预报的准确性不足；另外，模式方法往往还会在数据空间尺度范围上受到限制。支持向量机也是数值预报的常用方法[7]。但是单纯使用支持向量机方法进行预报容易出现空报、漏报情况。也有文献将支持向量机和粒子群算法结合对能见度进行预测[8]，但只适用于低能见度情形，且预测结果不稳健。

鉴于能见度的数据格式以及时空结构关系，本文基于广义高斯过程回归（GGPFR）[9]，从函数型数据分析[10]的角度对能见度进行研究。GGPFR模型能够针对非高斯型函数数据，通过联系函数和高斯过程，刻画响应变量与时间以及其他多维函数型协变量之间的非线性相关性，本质上属于一种非参数回归方法。另外，相比单一站点的预报模式，该模型能够对多个站点数据进行整体建模，充分挖掘数据信息之间的关联性。为了获得更好的拟合和预测效果，本文首先在模型上做出改进，建立带有混合效应的GGPFR（M-GGPFR）模型，充分利用参数结构的可释性和非参数结构的灵活性。本文首先阐述模型结构和估计、预测方法，并通过数值模拟研究展示推断效果，然后对我国20个城市进行建模，并重点分析北京市近五十年的能见度变化规律和趋势。

1 M-GGPFR模型

1.1 模型介绍

令{zm（t），t∈T}，m=1，2，…，M为第m批数据的函数型响应变量，服从指数族分布：

其中，αm（t） 和ϕm（t）都是函数型变量，αm（t）是点则参数，ϕm（t）为散度参数。

令xm（t）为一组函数型协变量，t表示时间或其他一维协变量。由此，函数型响应变量zm（t）与函数型协变量之间的关系刻画如下：

其中，βm（t）采用B样条近似估计，vm（t）、wm（t）和xm（t）为函数型协变量的子集，维数分别为r、k、Q。τm（t）为潜在随机过程，通过高斯过程回归（GPR）训练学习，协方差核函数采用平方指数形式，即

其中，θ=（vc，wc1，…，wcq）为一组超参数，可以通过训练数据集进行估计。

1.2 参数估计

设训练集共有M批数据，其中第m批数据观测的时间节点为为方便起见，记zm（tmi）、τm（tmi） 和xm（tmi） 为zmi、τmi和xmi，记Zm= （zm1，… ，zmNm）T，由此给出离散数据下的对应模型：

其中，EF（，）指代（1）式的指数族分布 ，这里为协方差矩阵。结构上，上述模型比GGPFR模型多了混合效应成分，其参数估计的具体实施仍可采用极大似然估计结合Laplace近似方法，详见文献[10]。

1.3 预测

以下考虑两种不同类型的预测问题。

（1）若预测对象来自训练样本，即已获得预测对象的部分数据，目的是预测未来时刻或缺失时刻数据。记在时刻{tk1，…，tkN}的观测值为Zk={zki，i=1，…，N} ，对应的输入向量为xk={xk1，…，xkN} 、vk（t）={vk1（t），…，vkN（t）}和wk（t）={wk1（t），…，wkN（t） } ，希望预测新时刻t*下的z*。令τ*=τk（t*）为t*时的潜在过程值，则给定τ*下z*的条件期望为：

于是有：

而Var（z*|D）可由下式计算

（2）若预测对象为全新样本，即没有该样本的响应变量信息，需要给出t*时刻下预测z*，上述方法不再适用。此时可以考虑一种加权预测的思想：先将该样本视为训练集中的已知样本，由此获得预测值，再通过z*属于第m样本的概率ωm进行如下加权预测：

如果没有信息推断ωm时，可以取ωm=1/M，m=1，…，M。

2 模拟研究

通过数值模拟检测模型估计和预测效果。数据产生自：

其中，对于每个m，xmi（=tmi）为将区间[- 4，4] 等分所得，vmi为取值为0或1的列向量，γ=1.0，η=100分别为vmi和wmi的系数。τmi是均值为0，协方差为（3）式所示的高斯过程，其中vc=0.5，wcq=0.1。响应变量zmi由（12）式生成：

由此生成40条曲线的样本，每条曲线包含41个样本点。为方便运算且避免协方差矩阵出现奇异，在模型最后加入一项噪声扰动。一般来讲，对于n类有序数组，若定义模型：

其中，b0=-∞，ba=∞，j=1，…，n-1为需要估计的阈值的密度函数可由下式得出：

进而运用经验贝叶斯方法估计出B样条系数、函数型协变量系数、超参数和阈值。

模拟部分n=4且阈值b12和b3均为未知参数。图1（a）展示了M-GGPFR模型对均值曲线的估计（选取第40条曲线进行展示），实线为真实曲线，虚线为估计曲线，表明估计有效。图1（b）比较了GGPFR模型（点线）与M-GGPFR模型（虚线）对β（t）（实线）的估计结果，结果表明混合效应模型估计更为准确。

图1 模型估计效果：（a）均值的拟合,（b）两种模型对 β（t）的拟合

进一步考虑预测问题，首先生成一条由40个点组成的全新曲线，随机选取一部分观测值估计模型，剩下的部分作为测试集进行预测（内含预测），同时比较两种模型的预测情况，如图2所示。不难发现，M-GGPFR预测效果明显优于GGPFR模型。

图2内含预测效果比较

图2 中实线（星型）为真实值，虚线（空圈）为预测值。其中：（a）GGPFR模型对ym（t） 的预测；（b）GGPFR模型对zm（t）的预测及95%的置信区间；（c）M-GGPFR模型对ym（t） 的预测；（d）M-GGPFR模型对zm（t）的预测及95%的置信区间。

接下来，考虑外延预测问题，选取[-4,0]时刻对应数据作为观测值，预测（0,4]时刻部分。图3展示两种模型的预测效果。相比而言，混合GGPFR模型优于GGPFR模型。

图3 外延预测效果比较（记号同前）

3 实例分析

本文分析中国20个城市1960—2012年的日观测能见度数据，包括能见度、降雨、风速和相对湿度。由于测量标准不同，能见度在1980年之前按0～9等级[11]记录，而1980年后记录具体的能见度值，2010年开始采取0～6级标准（见表1）。

表1 能见度等级和划分标准

表2给出样本下的能见度、风速、相对湿度及降雨之间相关系数，表明能见度与风速呈正相关，与相对湿度呈负相关。由此，将能见度作为响应变量，风速作为固定效应，降雨为随机效应，而时间和相对湿度用于训练（3）式中的协方差核函数。

表2 能见度、风速、相对湿度及降雨之间的相关系数

本文运用1981—2002年20个城市1月份的月平均数据数据作为训练集，以此预测后十年的能见度。下页图4是不同方法下北京市能见度趋势预测效果，包括普通线性回归、ARMA模型和M-GGPFR模型预测。

图4 北京市能见度趋势预测比较（单位：0.1km）

实线和星型表示真实测量值，虚线为预测值；实心圆表示M-GGPFR模型预测，三角形为时间序列预测，正方形为线性回归预测。

观察预测曲线并比较三者的均方根误差（RMSE）发现，M-GGPFR模型有最好的预测效果，其均方根误差为0.12，远小于时间序列模型的61.57和线性回归的72.02。

接下来，选取北京市1969—1989年冬季的数据展开研究。参照Zhang等[12]提出的能见度评价体系，日能见度20km以上视为“极好”天气，10～20km为“好”天气，2～10km为“坏”天气，低于2km为“特坏”天气。将北京市1969—1989年冬季的能见度进行频率划分，见下页图5。结果表明，北京市冬季约30%的天气为“好”天气，23%为“坏”天气。从1969—1989年的20年时间里，“特坏”天气从无到有，比重保持在5%以下，“坏”天气有所增加，而“好”天气在逐步减少，“极好”天气的频率基本保持稳定。可见，北京市的环境状况从80年代起已经开始受到一定的污染。

图5 北京市历年冬季不同可视范围下能见度频率柱状图

4 结论

本文基于M-GGPFR模型构建能见度与气象因子的关联并进行估计和预测分析。模型可以刻画气象因子与能见度的非线性相关性，而且可以借助不同站点数据信息，对目标站点进行多类型预测。数值分析说明模型很大程度改进了GGPFR模型的估计和预测效果，并且优于普通的回归模型和时间序列模型。在应用上，选取了部分城市能见度历史数据建模，并重点分析了北京市的能见度年际演变规律。由于能见度数据前后记录格式不同，常规方法很难获得年际演变规律。本文模型可以将等级数据背后的潜在数值进行恢复，进而在同一尺度下对数据进行比较分析。从结构上看，该模型可以视为广义线性模型与函数型数据模型的结合。除等级数据外，本文模型也适用于其他非高斯型函数数据，并且可以利用其进行聚类和判别分析。