岭回归算法在臭氧浓度集合预报中的应用

季 冕，程 龙

(安徽省环境监测中心站，安徽 合肥 230071)

臭氧(O3)是大气中常见的微量气体，主要分布在10 km至50 km的平流大气层中.高浓度的臭氧会刺激人体组织黏膜，对人体造成伤害.随着工业化发展进程的加快，臭氧污染已经成为一个严峻的社会问题，甚至在某些城市取代PM2.5成为空气污染的“罪魁祸首”[1].近地大气中臭氧的浓度变化有很显著的季节规律，很大程度上也有人为因素的影响，这使每日臭氧浓度变化具有很大的不确定性.

空气质量模型的研究发展经历了3个阶段[2]，第一代的箱式模型、高斯烟团等，第二代考虑复杂扩散过程的欧拉网格模型，到多区域多尺度的综合空气质量模型.目前，对臭氧等污染物的浓度预报的第三代模型，主要有两种[3]：一种是以大气动力学为基础，建立关于大气污染物浓度稀释扩散的数值模型，通过计算来模拟和预测大气污染物的动态分布，亦即数值预报；另一种是基于统计学方法，建立污染物浓度与气象参数间的统计预报模型，用来预测大气污染浓度，亦即统计预报.

目前，主流的空气预报模式主要有：美国环境技术公司(ENVIRON)提出的CAMx[4]、美国国家环境保护局(USEPA)[3]开发的第三代空气质量预报和评估系统(Models-3)中的CMAQ(community multiscale air quality model)、中科院大气物理研究所开发的嵌套网格空气质量预报模式NAQP(nested air quality prediction)[5]、基于美国环境预测中心(NCEP)和美国国家大气研究中心(NCAR)等科研机构开发的WRF(weather research forecast)的WRFC模式[3]等.

1 研究方法和数据来源

使用某种单一的模式独立进行气象数据预测，通常会带来很大的不确定性，预测偏差和方差往往比较大[6].因为每种方法有其特殊的适应性，为方便结合各模式的优点，很多方法将这些模式组合起来，构成一个新的预测模式，即集合预报模式.大气科学界比较著名的最优化集合预测方法(operational consensus forecasts, 简称OCF)是Woodcock和Engel提出的一种自动化的集合预测系统[7].通常选取一定时间作为预测的滑动窗口，对于一些基础模式的预测数据进行整合升级.

笔者以多个气象模式预测数据和真实观测数据为样本，利用已有预报模型的某种线性组合或凸组合来构成一种全新的臭氧预报模型，并对其性能进行量化评测.总体来讲，主要工作如下：首先实现基于Boosting的集合预测方法岭回归(ridge regression，简称RR)算法，针对算法中参数的特点进行分析，并给出参数的调整方案.其次进行算法性能的测试和评价.根据我国站点观测值和模式预报值，对OCF算法和RR算法分别给出相应的集合预报值，并通过均方根误差(RMSE)和时间序列进行性能评估.

原始数据来自中国环境监测总站的环境空气质量数值预报模式系统，原始数据集包含2015年9月1日至2017年2月19日的各个模式的预报结果和相应的观测值资料.

2 Ridge Regression算法介绍及实现策略

机器学习方法是人工智能领域研究的重要内容，集成机器学习的相关算法往往可以根据有限的观测数据，绕过难以通过数值方法直接得出的公式描述，得到所需的结果数据[8-9].笔者介绍的机器学习集合预测方法，可以很好地集成这些模式的优势所在，从而给出一个更加有效的预测结果.基于PAC可学习性原理中的Boosting算法[10-11]，该文给出了一个较好的思路.

机器学习旨在实现自动决策或预测.通常，需要基于初步估计步骤做出良好决策.然而，并不是所有算法都依赖于估计统计参数，笔者应用的序列集成技术即不采用估计的算法.

2.1 概 述

2.1.1 符号定义

2.1.2 序列集成技术

集成学习方法中，组合中每个模式获得的权重与模式过去预测值与观测值的差距有关系，也就是说理论上预测效果越好的模式将获得更高的权重.序列集成技术类似一个黑盒子，根据PAC可学习性理论，集成后的预测结果会优于单个预测模式(基学习器).单个的预测模式可以来自统计建模、数值模式或者决策树、神经网络等复杂的方法.所以，该方法实际是工作在元模型级别的.

(1)

2.2 算法详细说明

RR算法即岭回归算法[8]，通常设置1个参数λ≥0，并初始化向量u1=(0,…,0)，对于t≥2有

(2)

RR算法流程如下：

fort=1,2,…,T:

输出：uT=(uT[1],…,uT[N])

2.3 算法实现策略

RR算法实现的关键是通过迭代来不断更新模型权重，算法的核心代码如图1所示.在确定参数λ的值时，通常可以初始化一个遍历范围，然后以一定的步长进行遍历，找到给定精确度下最优的参数取值范围.

算法核心模块的输入可以是4个模式矩阵和1个观测矩阵，其格式可以为s×t(站点数、日期数)、罚分因子λ.算法输出为计算出的集成模式数值，其规模也为s×t，或者进一步根据实际观测矩阵(向量)给出该集成模式的均方根误差值(RMSE).

算法流程如下：

输入：P1、P2、P3、P4：4个模式矩阵，OB:观测矩阵，s：矩阵站点数，t：日期数，u：可初始化为0，保存加权向量，R：保存计算结果，A：保存每一步的中间变量Ai，AP:保存每一步Ai的广义逆矩阵，λ：罚分因子

fori=1,…,t:

v=u[i];#保存当天加权向量

forj=1,…,s

p=[P1[j,i];P2[j,i];P3[j,i];P4[j,i]]; #保存模型在特定t和s下的预测值

R[j,i]=v*p;

Re=R[j,i]-OB[j,i];

e_sum=e_sum+Re*p;

pd=p*p的逆;

pd_sum=pd_sum+pd;

end for

A[i+1]=A[i]+pd_sum; #用每天累加后的pd_sum更新A[i+1]

AP[i+1]=A[i+1]的广义逆矩阵；

u[i+1]=u[i]-(AP[i+1]*e_sum)的逆阵；

end for

OBA=t-30至t的OB值；

RA=t-30至t计算出的RR值；

输出：该参数λ下的RR矩阵R和其RMSE值res.

外层可以通过一个脚本方法，通过步长和区间的设置，寻找到合适的罚分因子λ.如果最后输出的是RMSE值，可以将每次算出的RMSE值存储到1个数组中，然后求得其中最小的RMSE罚分因子的λ值.RR算法的性能具有理论上的保证，算法性能的评测标准在下一部分详细说明.

3 性能评测与分析

3.1 评测指标与度量标准

均方根误差(root mean square error，简称RMSE)亦称为标准误差，其定义为：观测值与真实值差值的平方除以观测次数的平方根.它对一组测量中特大值和特小值比较敏感，故常用来反映精度.其表达式如下

(3)

其中：n为观测次数，xi为第i次的模式值，yi为第i次的真实值.RMSE值越小，说明模式的精度越高，与真实值贴近的越好.

(4)

其中：t0是评估开始时的第1个时间，T是评估结束的时间.

针对某个观测数据集和模式预测集，有以下几个度量标准：Bp，B，BX，BM[12].其中，Bp是所能获得的最小点误差，B是权重向量取值范围为RN时的误差，BX是权重向量取值为凸空间时的误差，BM是表现最好的模型M的误差.以上4种度量标准之间在样本不是太少的情况下，恒存在关系Bp≤B≤BX≤BM.

3.2 算法性能评测与分析

进行算法性能评测时，因为RR算法既可以应用在多个站点，也可以应用在单个站点.所以，分别对多站点的集合预测和单站点的集合预测进行了比对分析.最后，针对OCF算法、RR算法的预报性能，选取了几个站点进行时间序列的分析比对.

3.2.1 程序运行环境与整体架构

程序主要分为数据清洗模块和算法应用两大部分，所用的编程工具是Python 2.7和MATLAB 2016b，其中数据处理时还需要以下Python库pandas-0.19.2、numpy-1.12.0+mkl、nltk-3.2.2、scipy-0.18.1、XlsxWriter-0.9.5.Python的pandas库中包含的数据框格式(DataFrame)非常适于进行二维的数据处理.程序数据的处理流程如图1所示.

图1 程序数据处理流程图

3.2.2 多站点集合评测分析

多站点集合评测分析即是将多个站点的数据放到同一个数据集中，然后对这个数据集进行计算和处理.理想情况下，多个站点的选取最好是在同一区域中，这样可以对某个地区进行集合预测，理论上如果监测站点比较多，其他相邻地区可以通过插值的方法得到预测值，对区域的预测会更有针对性.

数据清洗采用两种力度：一种不除去重复模式值，一种除去重复模式值.以s代表参与计算的站点数，T代表日期数，BOCF和BRR分别代表OCF算法和RR算法对该数据集给出预测的预测矩阵的RMSE值.实验中OCF算法的窗口期为7 d，计算OCF算法和RR算法RMSE值的评测时段均为第8天至第T天(第1天权重都为0，无法给出预测值).对此，进行了如下几组测试.

(1) 数据初始规模s=30，T=90；不去重清洗后规模为s=11，T=89.此时得到表1中所列数据.

去重后数据规模为s=7，T=89.此时得到表2中所列数据.

(2) 数据初始规模s=90，T=90；不去重清洗后规模为s=40，T=89.此时得到表3中所列数据.

表3 s=90，T=90数据RMSE值计算(不去重)

去重后数据规模为s=22，T=89.此时得到表4中所列数据.

表4 s=90，T=90数据RMSE值计算(去重)

(3) 数据初始规模s=120，T=120；不去重清洗后规模为s=50，T=119.此时得到表5中所列数据.

表5 s=120，T=120数据RMSE值计算(不去重)

去重后数据规模为s=22，T=119.此时得到表6中所列数据.

表6 s=120，T=120数据RMSE值计算(去重)

通过以上实验结果分析可以得到：对于较长时间维度(90 d以上)，无论数据集是否去重，RR算法性能比OCF算法性能都要好很多，相比最优模式值BM,其优势也很明显，甚至比最优常数线性组合B的性能还要好.实际上，如果在较短时间维度上(如30 d)，RR算法的表现或许会比OCF算法差.根据集成学习的理论，这种情况往往是由训练数据过少所导致.

3.2.3 单站点效果评测分析

RR算法和OCF算法都可以应用于单站点.对于某些站点，将其作为一个集合得到的结果更优还是单个站点独立应用该模式得到的方法更优呢？为了探讨这个问题，进行如下测试实验.

随机选取10个站点的2015年9月1日起连续120 d的数据，进行清洗后，数据规模为s=10，T=119.评测RMSE值的时段均为第16天至第T天(对于RR算法，第1天权重都为0，无法给出预测值；而OCF算法需要相应的窗口期).

(1) 将多个站点作为1个集合应用RR算法和OCF算法(窗口期W分别取7和15)，可以得到表7中所列数据.

表7 多站点集合应用RR算法和OCF算法的RMSE值

(2) 对这些站点分别应用RR算法和OCF算法，则可以得到表8中所列数据.

表8 单站点分别应用RR算法和OCF算法的RMSE值

经过分析可以发现：

(1) OCF算法(窗口期为7 d)的性能比OCF算法(窗口期为15 d)的好.

(2) 从整体上看，RR算法性能比OCF算法要好.多站点集合应用时RR算法比OCF算法性能优势明显，单站点应用时，除第9号单站外，其他单站RR算法的应用效果相比OCF算法均有较大优势.

(3) RR算法可以在单站点上使用，多站点集合应用时计算出的RMSE值处于最优单站和最差单站之间(最差单站RMSE值比集合应用时要差).就所选样本而言，大多数站点单独应用RR算法时，性能都较好.

(4) 对于RR算法中参数λ的取值，其必然远离0，然而对于不同的站点、输入数据或不同的评测时长，λ值通常不相同.所以，应该根据实际情况确定该值的大小.

3.2.4 时间序列比对分析

为了对OCF算法和RR算法的性能有一个更加直观的表示，选取了几个站点，对模式数据和观测数据进行时间序列的效果比对和分析，结果如下.

(1) 随机抽取站点A，从2015年9月1日起，90 d时间维度内(实际清洗后为89 d)针对OCF算法、RR算法、实际观测值进行时间序列比对分析，结果如图2所示.

图2 站点A 90 d预测值与观测值对照图

(2) 随机抽取站点B，从2015年9月1日起，90 d时间维度内(实际清洗后为89 d)针对OCF算法、RR算法、实际观测值进行时间序列比对分析，结果如图3所示.

图3 站点B 90 d预测值与观测值对照图

(3) 随机抽取站点C，从2015年9月1日起，120 d时间维度内(实际清洗后为119 d)针对OCF算法、RR算法、实际观测值进行时间序列比对分析，结果如图4所示.

图4 站点C 120 d预测值与观测值对照图

(4) 随机抽取站点D，从2015年9月1日起，120 d时间维度内(实际清洗后为119 d)针对OCF算法、RR算法、实际观测值进行时间序列比对分析，结果如图5所示.

图5 站点D 120 d预测值与观测值对照图

(5) 随机抽取站点E，在2016年1月1日至2016年12月31日(实际清洗后为351 d)针对OCF算法、RR算法、实际观测值进行时间序列比对分析，结果如图6所示.

图6 站点E 2016年预测值与观测值对照图

(6) 随机抽取站点E，对2016年1月1日至2016年12月31日(实际清洗后为351 d)针对RR算法、4个模式值、观测值进行时间序列比对分析，结果如图7所示.

图7 站点E 2016年RR算法值、模式值与观测值对照图

通过以上时间序列的比对，可以看出，RR算法比OCF算法更加贴近于真实值，对于臭氧浓度趋势的变化适应性更强.同时，也印证了RR算法的RMSE值为何相比OCF小很多.然而，RR算法预测的精确度依赖于4个子模式，即对于某种极端值，若子模式都没有给出合适范围内的取值，则RR算法也很难直接给出这种极端值的合理预测.

3.2 算法的性能和优势讨论

在气象上，通常认为某个时次(某天)的误差是由一个较为稳定的系统偏差和一个确定性的扰动误差组成.OCF算法的偏差校正其实是对系统偏差的校正，系统偏差可以通过计算一段时间的误差平均获得，为了得到一个鲁棒性较好的平均误差(不易受到极端值的影响)，OCF算法往往选用分位数来获得系统偏差.但是，实际应用中，OCF算法并不能很好地满足假设条件，因此，只能从某种程度上减少极值带来的影响.

RR算法在算法性能上，相比OCF算法和传统模式(在这里作为集成模式的子模式)都有很大的提升(见图8所示)，而且其算法仅需要一个参数，避免了很多复杂的参数条件问题.尽管该参数在指定过程中可能存在一定的不确定性，但是通过遍历等方式并不难找到.而且RR算法的性能可以有理论上的保证，因此可以应用的范围也比较广.

RR算法通常可以把其他预测模式集成起来，提高整体性能.其优势可以归结为：对于回归问题，不需要构造拟合精度和预测能力都很好的回归算法，只要基预测器比随机猜测略好即可；它具有很好的通用性和鲁棒性，通常可以应用于任何基础回归算法，神经网络、决策树、支持向量机、线性回归、数值模式等方法都可以作为基学习器进行集成，而且理论上可以保证集成效果比单个学习器更好，不容易出现过拟合问题.当然，参与集成的个体数并不是越多越好，更多的个体需要更大的计算和存储开销，个体间的差异更难以获得.

图8 数据集RMSE对比图(不去重清洗/去重清洗)

4 结束语

论文实现了基于Boosting的集合预测方法岭回归Ridge Regression (RR)算法，针对算法中参数的特点进行了分析，并给出了调整参数的方案.进行了算法性能的测试和评价.根据我国站点观测值和模式预报值，对OCF算法和RR算法分别给出相应的集合预报值，并通过均方根误差(RMSE)和时间序列进行性能评估.发现通常情况下，与已有的预报模式相比，RR算法具有更高的准确度和更好的稳定性.

基于机器学习技术的Ridge Regression算法是一个元模式级别的算法，它可以将多种模式集成在一起，起到“博采众长”的作用，它的集成性能和目标预测的空气污染物没有相关性，即也可以进行其他气象指标的预测应用.同时，它的预测准确性与这些基学习器紧密相关.现阶段，传统气象预报模式大都基于数值模式和物理过程的模拟，如果能够通过深度学习等机器学习技术，探讨出一个更强的分类器，无疑将会对集合预测的精度带来很大提升.

进一步的研究如果能够从算法层面更好地解决权重更新和参数选取问题，将对更大规模的数据设计高效的并行算法，为臭氧等大气污染物的预测带来新的机遇与挑战.当预测精度达到一定程度后，可以用来反推污染源位置等信息，跟踪生产方式带来的影响等，这将为各种空气污染物的治理提供解决思路.