基于小波去噪与Catboost的短时天气预报方法与流程

本发明涉及天气预报领域，特别是涉及一种基于小波去噪与Catboost的短时天气预报方法。

背景技术：

气象因素的变化(如风速、温度、湿度、降水等)都深刻地影响着人类的生活。准确预报未来气象要素，可广泛服务于人们日常生活(如穿衣着装)，交通运输(如航班起降)，农林畜牧业(如水产养殖)，致灾天气避险(如台风预警)等领域。随着地球观测卫星数量增长与气候模型日益增强，气象研究者们面临着更大规模的数据。机器学习可以在数据量增长时提升预测性能。一个高分辨率的气候模型的一次运行即可以产生千万亿字节的数据。近年发展迅速的深度学习模型也适用于天气预报中的时空序列预测问题。

目前，数值预报与基于人工智能的预报是天气预报的主要方法。对于数值天气预报方法而言，短期预报需要复杂的物理大气模型仿真。近年来，机器学习与深度学习已开始被应用于天气预报。比如，深层卷积神经网络被应用于检测气候数据集中的极端天气。多层长短期记忆(LSTM)模型也被广泛应用于时间序列问题。机器学习中基于决策树的模型，能有效地解决大数据问题，同时训练时间也较短。但是，现有技术中的利用机器学习和深度学习进行天气预报的方案部分存在模型训练收敛时间较长、影响实际预测效率的问题。

技术实现要素：

发明目的：本发明的目的是提供一种基于小波去噪与Catboost的短时天气预报方法，能够解决现有技术中存在的“模型训练收敛时间长、影响实际预测效率”的技术问题。

技术方案：为达到此目的，本发明采用以下技术方案：

本发明所述的基于小波去噪与Catboost的短时天气预报方法，包括以下步骤：

S1：输入t时刻的历史气候特征数据，包含t时刻模式预测的特征数据M1,…,Mm和t时刻实际观测的特征数据O1,…,On，其中，Ms表示t时刻模式预测的第s个特征数据，1≤s≤m，m表示t时刻模式预测的特征数据的总数，Oi表示t时刻实际观测的第i个特征数据，1≤i≤n，n表示t时刻实际观测的特征数据的总数；对由时刻t、O1-On和M1-Mm组成的输入数据进行数据清洗；

S2：对O1-On和M1-Mm进行排序，按照重要性由高到低依次赋予以下分值：m+n分,m+n-1分,...,1分，然后剔除分值低于Q分的特征数据，Q的值预先设定；

S3：对待预测气候特征序列的P个站点进行one-hot编码，完成空间特征添加；对待预测气候特征序列的时间信息进行时钟投影以得到时间特征；

S4：对待预测气候特征序列中的距地面2米高度处的温度、距地面2米高度处的相对湿度以及距地面10米高度处的风速进行小波去噪；

S5：将模式预测的特征数据M1,…,Mm、待预测气候特征序列、步骤S4得到的小波去噪后的待预测气候特征序列、待预测气候特征序列的真实标签值输入Catboost模型，调整树的深度、树的最大数量与迭代次数，得到训练后的Catboost模型，然后将测试集输入到训练后的Catboost模型中，从而输出距地面2米高度处的温度、距地面2米高度处的相对湿度以及距地面10米高度处的风速的预测结果。

进一步，所述步骤S1中的数据清洗包括缺省值填充和异常值删除这两个步骤。

进一步，所述缺省值填充步骤为：将t时刻实际观测的特征数据用t+1时刻实际观测的特征数据与t-1时刻实际观测的特征数据的均值或者t时刻模式预测的特征数据进行填充，将t时刻模式预测的特征数据用t+1时刻模式预测的特征数据与t-1时刻模式预测的特征数据的均值或者t时刻实际观测的特征数据进行填充。

进一步，所述步骤S3中，时间特征中的月份特征Month_new根据式(1)得到：

式(1)中，Month表示步骤S1中时刻t所对应的月份。

进一步，所述步骤S5中，Catboost模型中的损失函数选择交叉熵损失函数。

进一步，所述步骤S4中，去噪所使用的滤波器包括小波滤波器和尺度滤波器；对待预测气候特征序列中的距地面2米高度处的温度进行小波去噪的过程包括以下步骤：

S41：待预测气候特征序列中的距地面2米高度处的温度所对应的历史序列的第j级小波系数和第j级尺度系数根据式(2)和式(3)得到：

其中，t1表示时间，Lj＝(2^j-1)(L1-1)+1，Lj表示第j级小波滤波器的长度，L1表示第一级小波滤波器的长度，尺度滤波器和小波滤波器的长度相等，hj,l表示第j级小波滤波器的滤波器函数中的第l个函数值，gj,l表示第j级尺度滤波器的滤波器函数中的第l个函数值，表示历史序列中t1-lmodN时刻的元素，N为历史序列中的时刻总数；

S42：对第j级小波系数和第j级尺度系数均进行阈值处理，再将阈值处理后的第j级新的小波系数和第j级新的尺度系数进行逆离散小波变换，从而得到去噪后的距地面2米高度处的温度的历史序列。

进一步，所述步骤S42中，对第j级小波系数进行阈值处理得到第j级新的小波系数的过程如式(4)所示：

式(4)中，λj为第j级小波变换的阈值。

有益效果：本发明公开了一种基于小波去噪与Catboost的短时天气预报方法，相比现有技术，能够提高预测的精确度，减少模型训练的收敛时间，提高预测效率。

附图说明

图1为本发明具体实施方式中步骤S3的示意图；

图2为本发明具体实施方式中步骤S4的示意图；

图3为本发明具体实施方式中方法的流程图；

图4为本发明具体实施方式中实施例1的方法与现有技术中方法的预测结果对比图；

图4(a)为距地面2米高度处的温度的预测结果对比图；

图4(b)为距地面2米高度处的相对湿度的预测结果对比图；

图4(c)为距地面10米高度处的风速的预测结果对比图。

具体实施方式

下面结合具体实施方式和附图对本发明的技术方案作进一步的介绍。

本具体实施方式公开了一种基于小波去噪与Catboost的短时天气预报方法，如图3所示，包括以下步骤：

S1：输入t时刻的历史气候特征数据，包含t时刻模式预测的特征数据M1,…,Mm和t时刻实际观测的特征数据O1,…,On，其中，Ms表示t时刻模式预测的第s个特征数据，1≤s≤m，m表示t时刻模式预测的特征数据的总数，Oi表示t时刻实际观测的第i个特征数据，1≤i≤n，n表示t时刻实际观测的特征数据的总数；对由时刻t、O1-On和M1-Mm组成的输入数据进行数据清洗；

S2：采用递归特性消除、相关特性分析或者基于树模型的特征重要性排序对O1-On和M1-Mm进行排序，按照重要性由高到低依次赋予以下分值：m+n分,m+n-1分,...,1分，然后剔除分值低于Q分的特征数据，Q的值预先设定；

S3：对待预测气候特征序列的P个站点进行one-hot编码，完成空间特征添加；对待预测气候特征序列的时间信息进行时钟投影以得到时间特征；如图1所示；

S4：对待预测气候特征序列中的距地面2米高度处的温度、距地面2米高度处的相对湿度以及距地面10米高度处的风速进行小波去噪；如图2所示；

S5：将模式预测的特征数据M1,…,Mm、待预测气候特征序列、步骤S4得到的小波去噪后的待预测气候特征序列、待预测气候特征序列的真实标签值输入Catboost模型，调整树的深度、树的最大数量与迭代次数，得到训练后的Catboost模型，然后将测试集输入到训练后的Catboost模型中，从而输出距地面2米高度处的温度、距地面2米高度处的相对湿度以及距地面10米高度处的风速的预测结果。

步骤S1中的数据清洗包括缺省值填充和异常值删除这两个步骤。缺省值填充步骤为：将t时刻实际观测的特征数据用t+1时刻实际观测的特征数据与t-1时刻实际观测的特征数据的均值或者t时刻模式预测的特征数据进行填充，将t时刻模式预测的特征数据用t+1时刻模式预测的特征数据与t-1时刻模式预测的特征数据的均值或者t时刻实际观测的特征数据进行填充。

步骤S3中，时间特征中的月份特征Month_new根据式(1)得到：

式(1)中，Month表示步骤S1中时刻t所对应的月份。

步骤S5中，Catboost模型中的损失函数选择交叉熵损失函数。

步骤S4中，去噪所使用的滤波器包括小波滤波器和尺度滤波器；对待预测气候特征序列中的距地面2米高度处的温度进行小波去噪的过程包括以下步骤：

S41：待预测气候特征序列中的距地面2米高度处的温度所对应的历史序列的第j级小波系数和第j级尺度系数根据式(2)和式(3)得到：

其中，t1表示时间，Lj＝(2^j-1)(L1-1)+1，Lj表示第j级小波滤波器的长度，L1表示第一级小波滤波器的长度，尺度滤波器和小波滤波器的长度相等，hj,l表示第j级小波滤波器的滤波器函数中的第l个函数值，gj,l表示第j级尺度滤波器的滤波器函数中的第l个函数值，表示历史序列中t1-lmodN时刻的元素，N为历史序列中的时刻总数；

S42：对第j级小波系数和第j级尺度系数均进行阈值处理，再将阈值处理后的第j级新的小波系数和第j级新的尺度系数进行逆离散小波变换，从而得到去噪后的距地面2米高度处的温度的历史序列。

步骤S42中，对第j级小波系数进行阈值处理得到第j级新的小波系数的过程如式(4)所示：

式(4)中，λj为第j级小波变换的阈值。

下面以一个实施例为例，对本具体实施方式进行进一步的阐述。

实施例1：

本方法验证数据集为2018AI全球挑战赛提供的气候特征数据集。“观测”和“睿图”数据集，均包含北京市10个气象观测站点，约3年多的数据，连续性较好，缺失样本较少。“观测”集逐时记录当前气象观测站点的9个地面气象要素，通过气象仪器实时监测得到；“睿图”集包含地面和特征气压层共计29个气象要素，由数值预报模式在超级计算机上运算产生，其在每天03:00(北京时11:00)启动区域数值模式，预报至第二天15:00(北京时23:00)，共计37个时次(00–36)。

其中训练集的日期为2015年3月1日3时至2018年5月31日3时，验证集的日期为2018年6月1日3时至2018年8月28日3时，测试集为2018年8月29日3时至2018年11月3日3时。预测精度采用均方根误差RMSE和偏差BIAS作为评价指标，评测样本为北京10个观测站整个评测期内每小时产生的数据样本。

其中n为评测样本总数，为第i个样本的实际观测值，为第i个样本的模型预测值，RMSE(M)表示数值天气预报模式数据与真实数据的均方根误差，RMSE(model)表示模型预测数据与真实数据的均方根误差，总得分会先计算三个预测指标的得分后求平均值。上述评价标准中，以RMSE为首选标准，在相同RMSE得分的前提下，进一步参考BIAS评测预报结果的优势。

本方法步骤S1中，输入数据为3年历史气候数据，2015年3月1日3时至2018年5月31日3时，包含29种模式预测的特征数据M1,…,M29，9种实际观测的特征数据O1,…,O9。步骤S2中，对O1-O9和M1-M29进行排序，按照重要性由高到低依次赋予以下分值：38分，37分，……，1分，然后剔除对所需预测的特征影响最小的特征数据。步骤S3中，对待预测气候特征序列的10个站点进行one-hot编码，完成空间特征添加；对待预测气候特征序列的时间信息进行时钟投影以得到时间特征。步骤S5中，树深度设为10，树的最大数量设为1000，迭代次数设为3000次。

图4(a)-图4(c)为本实施例的模型预测结果与其他方法对比图，图中显示时间为UTC世界标准时间，其中Catboost这条曲线表示本实施例方法的预测结果。表1也显示出本实施例方法与现有技术中其他方法预测结果的对比。

表1本实施例预测分数与其他方法对比结果