CN105974495B - 利用分类拟合法预判目标区域未来平均云量的方法 - Google Patents

Abstract

本发明提供了一种利用分类拟合法预判目标区域未来旬平均云量的方法。该方法中，对FY‑2E数据进行小波分解，分别提取低频趋势时间序列与高频随机时间序列，分别将ARMA系列时间序列模型应用于低频趋势时间序列的预测，将神经网络模型应用于高频随机时间序列的预测，通过整合得到预测的未来旬平均云量，结果可以供遥感卫星订购用户做参考，预判目标区域有效数据获取的可行性，从而提高了预测的准确度。

Description

利用分类拟合法预判目标区域未来平均云量的方法

技术领域

本发明涉及卫星遥感技术领域，尤其涉及一种利用分类拟合法预判目标区域未来平均云量的方法。

背景技术

云量预测对于光学遥感增加获取无云或者少云卫星数据的可能性，获取感兴趣区域无云或者少云时间段的卫星数据，提高卫星资源利用率以及卫星数据获取率具有重要的指导意义。在传统的卫星数据成像规划和数据获取能力评估等情况应用中，国际卫星云气候计划(ISCCP)提供的云量数据一直是国内外研究者广泛研究的资料集。近年来，随着遥感技术的发展，越来越多的遥感应用需要高质量、高分辨率的卫星遥感数据为其提供空间信息支持。面对新的数据需求，ISCCP的低空间分辨率、低时间分辨率的特性大大限制了云量预测的应用，与此同时，高空间分辨率、高时间分辨率的FY-2E数据逐渐成为云量预测的研究热点。然而，FY-2E数据集与传统ISCCP数据集的差异性使得基于ISCCP的时间序列云量预测模型不再适用于FY-2E，这制约了新数据在云量预测上的广泛应用。因此，基于FY-2E数据集的云量预测模型的研究对于高分辨率卫星数据获取和应用研究具有重要的意义。

FY-2即风云二号气象卫星是我国自行研制的第一代地球同步轨道气象卫星，与极地轨道气象卫星相辅相成，构成我国气象卫星应用体系。这些卫星作用是获取白天可见光云图、昼夜红外云图和水气分布图，进行天气图传真广播，收集气象、水文和海洋等数据收集平台的气象监测数据，供国内外气象资料利用站接收利用，监测太阳活动和卫星所处轨道的空间环境，为卫星工程和空间环境科学研究提供监测数据。FY-2E云量数据集以经纬度投影以及等面积投影两种投影方式，将全球按照不同的空间分辨率划分为网格阵列，记录了日尺度数量级、旬尺度数量级以及月尺度数量级三种时间分辨率下的云量数据。

以FY-2E云量数据集为例，选取等面积投影方式，对全球范围按照标准网格划分，例如140km*140km分辨率下，全球共划分为26410个网格，其中，覆盖中国的有效网格为626个。根据FY系列卫星多年来收集的数据，国内外学者对云的空间和时间分布及其变化规律进行了大量的研究，但基本不涉及未来云量的预判问题。

现有技术一般采用基于时间序列特征分类的ARMA系列时间序列模型拟合方法计算未来一年月平均云量的方法。如图1所示，以FY-2E的140km*140km分辨率数据集为例，该方法包括：

步骤S101，从数据集中获取某一网格的n年某年的旬平均云量序列C_i，j。

其中，i为网格编号，取值范围[1，626]；j为旬序列，取值范围[1，36]，C_i，j为36×n个数组成的序列；

步骤S102，采用ADF单位根检验法检验所述旬平均云量序列C_i，j的平稳性，如果所述旬平均云量序列C_i，j的数据平稳，则为平稳时间序列，执行子步骤S103；如果所述旬平均云量序列C_i，j的数据不平稳，则执行步骤S105

步骤S104，利用ARMA自回归滑动平均模型对旬平均云量进行处理，得到未来旬平均云量。

步骤S105，采用36步差分运算对判定为非平稳序列的旬平均云量序列C_i，j进行差分计算，检验序列的周期性。如果序列具有周期性特点，并且在滞后期36的整数倍出现峰值，则存在季节特征，执行步骤S106，否则执行步骤S107

步骤S106，利用SARIMA季节性差分自回归滑动平均模型对旬平均云量进行处理，得到未来旬平均云量。

步骤S107，利用ARIMA自回归积分滑动平均模型对旬平均云量进行处理，得到未来旬平均云量。

然而，在实现本发明的过程中，申请人发现上述的预判未来月云量的方法存在如下缺陷：上述方法仅适用于时间分辨率与空间分辨率都较低的ISCCP数据，其数据空间分辨率为280km*280km，时间分辨率为月，对于时间分辨率与空间分辨率大幅度提高，空间分辨率最高达到35km*35km，时间分辨率为旬，数据集体量明显丰满的FY-2E数据集，预测准确度较低。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种利用分类拟合法预判目标区域未来平均云量的方法。

(二)技术方案

本发明利用分类拟合法预判目标区域未来平均云量的方法包括：

步骤A：读取云气候学数据集中目标区域所在网格的往年的平均云量序列C；

步骤B：采用小波分解的分析方法，将平均云量序列C进行分解与重构，得到低频趋势时间序列与高频随机时间序列σ^*(s)，其中，s＝1，2，3...，S，S为平均云量序列C中的云量数据的个数；

步骤C：对于低频趋势时间序列和高频趋势时间序列σ^*(s)，分别建立模型进行拟合预测，得到相应的低频趋势时间序列的预测结果和高频趋势时间序列的预测结果以及

步骤D：对低频趋势时间序列的预测结果和高频趋势时间序列的预测结果进行整合，得到最终云量预测的时间序列实现目标区域未来平均云量的预测。

(三)有益效果

从上述技术方案可以看出，本发明利用分类拟合法预判目标区域未来旬平均云量的方法具有以下有益效果：

(1)对FY-2E数据进行小波分解，分别提取低频趋势时间序列与高频随机时间序列，能够规避密集时间节点带来的大量云量数据对未来云量变化趋势的掩盖，也能够充分利用数据量增加带来的提高预测准确度的优势；

(2)自回归滑动平均模型(ARMA)、自回归积分滑动平均模型(ARIMA)和季节性求和自回归滑动平均模型(SARIMA)能够充分利用时间序列自身的规律性，得出较为准确的预判值；神经网络对随机无序时间序列的高拟合性能保证对高频随机时间序列进行拟合并作出较为准确的预判。本发明结合历史旬平均云量序列特征，分别将ARMA系列时间序列模型应用于低频趋势时间序列的预测，将神经网络模型应用于高频随机时间序列的预测，通过整合得到预测的未来旬平均云量，结果可以供遥感卫星订购用户做参考，预判目标区域有效数据获取的可行性。

附图说明

图1为本发明实施例利用FY-2E云气候学数据预测未来旬云量的方法的流程图；

图2为对于北京地区利用2009年到2012年4年数据预判2013年36个旬的预判图；

图3为北京地区2013年采用传统方法与本发明WTNC方法进行云量预测结果比较图；

图4为北京地区2013年本发明WTNC方法云量变化斜率95％置信区间评判图。

具体实施方式

本发明提供了一种适用于高空间分辨率云气候学数据利用分类拟合法预测未来平均云量的方法，提高原有方法预判结果的准确性。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

在本发明的一个示例性实施例中，提供了一种利用FY-2E云气候学数据通过分类拟合法预测未来旬云量的方法。请参照图1，该方法包括四个步骤：

步骤A：读取FY-2E等面积投影数据集中140km*140km分辨率数据目标区域所在网格的预设年份段【2009-2012】的旬平均云量序列C：

其中，C为4×36行1列的向量，i为年份，j代表旬编号。

读取目标区域所在网格的云气候学数据中的旬平均云量数据可以按照FY-2E数据集的说明进行。

步骤B：采用小波分解的分析方法，将旬平均云量序列C进行分解与重构，选用db4小波，得到两类时间序列，即低频趋势时间序列与高频随机时间序列σ^*(t)；

该将云气候数据进行小波分解的步骤B进一步包括：

子步骤B1：指定分解尺度n，通过小波分解，利用Mallat算法将旬平均云量序列C逐级极限逼近，分解得到低频序列的小波系数数据与高频序列的小波系数数据σ；

其中，L为系数序列的项数，具体取值与旬平均云量序列C计算有关。

子步骤B2：对两类小波系数数据和σ分别使用Mallat算法逆运算，重构到原尺度上，分别得到低频趋势时间序列和高频随机时间序列σ^*(s)。

此步骤实现了旬平均云量序列C的分解拆分，其中，低频趋势时间序列与高频随机时间序列σ^*(s)的叠加可以得到完整的旬平均云量序列C，即：

步骤C：对于低频趋势时间序列和高频趋势时间序列σ^*(s)，分别建立模型进行拟合预测，预测未来t时间节点的平均云量，得到相应的预测结果和其中，t＝1，2，3...，T，T为预测未来云量数据的个数；

该步骤C包括：

子步骤C1：将低频趋势时间序列的数据作为模型训练样本，利用ARMA系列时间序列模型拟合预测，预测未来t时间节点的平均云量，得到低频趋势时间序列的预测结果

该子步骤C1进一步包括：

子分步骤C1a：建立检验的原假设H₀：低频趋势时间序列是非平稳序列；相应的备择假设H₁：低频趋势时间序列是平稳序列；

子分步骤C1b，查找ADF检验法统计表，对比检验低频趋势时间序列的统计量t_τ与显著性水平α相应的分位点，如果小于分位点，拒绝原假设，低频趋势时间序列为平稳序列；如果大于该分位点，则接受该原假设，低频趋势时间序列为非平稳序列；

其中，统计量t_τ是ADF检验得出的，显著性水平α是查找ADF检验法统计表得出的，本领域技术人员非常清楚其计算过程，此处不再赘述。

子分步骤C1c，当低频趋势时间序列为非平稳序列时，对低频趋势时间序列进行s＝36的差分计算，得到新的时间序列对重复子分步骤C1a与C1b，如果为平稳序列，则低频趋势时间序列为季节性非平稳序列；如果为非平稳序列，则低频趋势时间序列为普通非平稳序列；

子分步骤C1d，对平稳序列的选择自回归滑动模型，即ARMA模型，进行以T为预测步长的拟合；对于普通非平稳序列选择自回归积分滑动平均模型，即ARIMA模型，进行以T为预测步长的拟合；或对于季节性非平稳序列选择季节性差分自回归滑动平均模型，即SARIMA模型，进行以T为预测步长的拟合，从而得到低频趋势时间序列的预测结果

子步骤C2：对高频随机时间序列σ^*(s)进行神经网络模型拟合预测，预测未来t时间节点的平均云量，得到高频趋势时间序列的预测结果

该子步骤C2进一步包括：

子分步骤C2a，提取高频随机时间序列σ^*(s)中的144个时间节点的云量数据；

子分步骤C2b，以该144个时间节点的云量数据为初始输入，选用BP算法进行权值修正，设第k步系统的实际输出为y_d(k)，计算网络的学习指标函数即误差平方和函数，得到误差平方和E(k)；

误差平方和E(k)的计算过程参见下述方程：

其中，y(k)为第k步神经网络系统输入的数值，y_d(k)为第k步神经网络系统输出的数值，n为总步数。

子分步骤C2c，根据梯度下降法，分别计算误差平方和E(k)对权值的偏导数，使误差沿梯度方向下降，即令偏导数为0，得到Elman网络的学习算法；

式中，w¹，w²和w³分别表示神经网络中间层到输出层、输入层到中间层，承接层到中间层的连接权值，η₁，η₂，η₃分别是Elman网络中w¹，w²和w³的学习步长，x_j(k)为第j次迭代在第k步的初始值；代表i，j次迭代中承接层到中间层的权值差；u_q(k-1)为第(k-1)步(即上一步)q层中间值；为i输出层的校验系数；为中间层的校验系数；y_d，i(k)为k步神经网络输出数值；g′_i(·)为输出神经元传递函数；f′_j(·)代表中间层神经元的传递函数。

子分步骤C2d：使用以下参数为Elman网络的学习算法赋阈值：输入单元层m＝144，输出单元层n＝3，隐含层单元数r＝10，最大迭代次数mix＝2000，训练误差指标θ＝0.15；

子分步骤C2e：以2009至2012年4年的144个云量数据作为144个初始值输入数据，让神经网络模型对时间序列进行训练学习，当达到设定的阈值标准时，即存在确定的训练模型权值，由此权值得到36维结果输出，即获得高频随机时间序列预测结果

步骤D：对低频趋势时间序列的预测结果和高频趋势时间序列的预测结果进行整合，得到最终云量旬平均预测的时间序列实现目标区域未来平均云量的预测。

该步骤D包括：

子步骤D1，整理低频趋势时间序列预测结果与高频随机时间序列预测结果

此子步骤属于数据准备，本领域技术人员很清楚数据准备的方法，此处不再重述，另外，在某些情况下，该数据准备步骤也可以省略。

子步骤D2，通过叠加计算，将低频趋势时间序列的预测结果与高频趋势时间序列的预测结果整合在一起，得到最终云量旬平均预测的时间序列

至此，得到最终云量旬平均预测的时间序列预测过程结束。

依据以上描述，本领域技术人员应当对本发明利用FY-2E云气候学数据预测未来旬平均云量的方法有了清楚的认识。以下给出应用本实施例方法进行北京地区的云量预判具体结果。

实际应用实例：

以北京所在网格为例，在时间上，选取2013年份进行预判和验证。用2009年到2012年5年数据预判2013年36个旬云量。

将传统方法(即单纯使用ARMA系列时间序列模型进行拟合预测的方法)的预测结果、本发明方法(a Wavelet Transform，Time series model，Neural network combinedmethod，简称WTNC方法)的预测结果和实际云量值进行比较，结果如图3所示。

请参照图3，横坐标轴表示时间的旬节点，纵坐标轴表示云量所占百分比，实线线段为实际云量值，虚线线段为传统方法预测值，点线段为本发明WTNC方法的预测值。从图中可以看出，整体上，传统方法并没有反映出云量的变化趋势，尤其在第20旬到30旬部分，预测效果偏差较大。另一方面，在细节上，传统方法的预测结果并没有体现细节信息，多以平滑的曲线描绘掩盖了波动细节，同时出现了大量错误的波动预测，而本发明WTNC方法不仅如实反映了变化趋势，而且在大部分点上准确地刻画了波动变化信息，相较传统方法，具有明显的优势。

通过云量趋势符合度K值对预测云量趋势上升或下降变化正确的节点数进行定量的判断。由定义可知，云量趋势符合度K值反映了预测结果符合真实波动情况的百分比，K值越高，云量预测变化的符合度也就越高，反之，K值越低，云量预测变化的结果就越不能反映真实云量的变化趋势。以下是北京458网格传统方法与本发明WTNC方法预测结果的云量趋势符合度K值的统计表。

表1北京地区传统方法与本发明WTNC方法预测结果参数比较

由上述计算结果分析可知，与传统方法相比，本发明WTNC方法在平均绝对误差，平均相对误差，以及均方根误差上都有所降低，其中，平均绝对误差降低了1.36，平均相对误差降低了0.03，均方根误差降低了2.66。平均绝对误差、平均相对误差以及均方根误差都是预测得到的一整年结果中每一个时间节点，与当年实际云量中对应时间节点云量的误差计算累计结果，因此，这三个数值体现了整体预测结果的误差情况。对比结果表明，从以上三个方面的全方位考量，本发明WTNC方法整体预测结果的符合度较传统方法有一定的提升。对于云量趋势符合度K值，本发明WTNC方法高出近23％，趋势符合度显著提高，证明本发明WTNC方法在时间序列节点方向变化的预测效果上有明显的优势。

由符合度K值得到了预测时间节点变化趋势正确率的统计，再对预测时间节点上升与下降的程度进行进一步定量的评估。计算预测时间节点变化的斜率，并与实际真实变化的斜率进行比较，结合95％置信区间进行进一步验证。结果如图4所示。

图4中，横坐标轴表示时间的旬节点，纵坐标轴表示云量所占百分比，实线线段为实际云量值，虚线线段为实际云量数值95％置信区间，点为本发明WTNC方法的预测结果各时间节点斜率值。从图中可知，实际云量变化斜率值绝大部分落在预测云量斜率值的95％置信区间内，置信区间内落点达到91.27％，仅有3个点落在置信区间之外，也就是说，有95％的把握认为，预测的斜率值是在可接受范围之内。预测结果的斜率值十分理想。

由以上试验可知，对于北京地区，本发明WTNC方法的预测结果优于传统方法，并且与真实值的变化趋势相一致，取得了理想的效果。

此外，上述方法并不仅限于实施方式中提到的各种具体算法，本领域的普通技术人员可对其进行简单地熟知地替换，例如：

(1)步骤B中db4小波还可以使用其他小波代替；

(2)步骤C2中参数估计环节提到的最速梯度法还可以使用高斯牛顿法或阻尼最小二乘法代替。

(3)步骤C2中的阈值设定环节可以设置其他阈值修正模型权值；

(4)除FY-2E数据集外，本发明方法还能适用于其他的同类数据，空间分辨率小于280km*280km，时间分辨率小于月，以时间序列形式呈现的云量数据。

综上所述，本发明基于小波分解的分析理论，将云量时间序列分解为低频趋势时间序列与高频随机时间序列，选用不同的模型来分别预测未来旬平均云量，具有适用于高分辨率数据、可靠性高、应用广泛的优势，可广泛应用于经济、气象、农业、生物科学等领域，具有较好的应用前景。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种利用分类拟合法预判目标区域未来平均云量的方法，其特征在于，包括：

步骤A：读取云气候学数据集中目标区域所在网格的往年的平均云量序列C；

步骤B：采用小波分解的分析方法，将平均云量序列C进行分解与重构，得到低频趋势时间序列与高频随机时间序列σ^*(s)，其中，s＝1，2，3...，S，S为平均云量序列C中的云量数据的个数；

步骤C：对于低频趋势时间序列和高频趋势时间序列σ^*(s)，分别建立模型进行拟合预测，得到相应的低频趋势时间序列的预测结果和高频趋势时间序列的预测结果以及

步骤D：对低频趋势时间序列的预测结果和高频趋势时间序列的预测结果进行整合，得到最终云量预测的时间序列实现目标区域未来平均云量的预测；

其中，所述步骤B的采用小波分解的分析方法中，选用db4小波，所述步骤B包括：子步骤B1：指定分解尺度n，通过小波分解，利用Mallat算法将平均云量序列C逐级极限逼近，分解得到低频序列的小波系数数据与高频序列的小波系数数据σ；以及子步骤B2：对两类小波系数数据和σ分别使用Mallat算法逆运算，重构到原尺度上，分别得到低频趋势时间序列和高频随机时间序列σ^*(t)；

所述步骤C包括：子步骤C1：将低频趋势时间序列的数据作为模型训练样本，利用ARMA系列时间序列模型拟合预测，预测未来t时间节点的平均云量，得到低频趋势时间序列的预测结果以及子步骤C2：对高频随机时间序列σ^*(s)进行神经网络模型拟合预测，预测未来t时间节点的平均云量，得到高频趋势时间序列的预测结果其中，t＝1，2，3...，T，T为预测未来云量数据的个数；

其中，所述平均云量序列C和最终云量预测的时间序列中时间分辨率均为旬尺度数量级；

所述步骤A中的云气候学数据集为风云二号气象卫星记录的云量数据集，其云量数据的空间分辨率小于280km*280km。

2.根据权利要求1所述的方法，其特征在于，所述子步骤C1包括：

子分步骤C1a：建立检验的原假设H₀：低频趋势时间序列是非平稳序列；相应的备择假设H₁：低频趋势时间序列是平稳序列；

子分步骤C1b，查找ADF检验法统计表，对比检验低频趋势时间序列的统计量t_τ与显著性水平α相应的分位点，如果小于分位点，拒绝原假设，低频趋势时间序列为平稳序列；如果大于该分位点，则接受该原假设，低频趋势时间序列为非平稳序列；

子分步骤C1c，当低频趋势时间序列为非平稳序列时，对低频趋势时间序列进行s＝36的差分计算，得到新的时间序列对重复子分步骤C1a与C1b，如果为平稳序列，则低频趋势时间序列为季节性非平稳序列；如果为非平稳序列，则低频趋势时间序列为普通非平稳序列；

子分步骤C1d，对平稳序列的选择自回归滑动模型进行以T为预测步长的拟合；对于普通非平稳序列选择自回归积分滑动平均模型进行以T为预测步长的拟合；或对于季节性非平稳序列选择季节性差分自回归滑动平均模型进行以T为预测步长的拟合，从而得到低频趋势时间序列的预测结果

3.根据权利要求1所述的方法，其特征在于，所述子步骤C2包括：

子分步骤C2a，提取高频随机时间序列σ^*(s)中的S个时间节点的云量数据；

子分步骤C2b，以该S个时间节点的云量数据为初始输入，选用BP算法进行权值修正，设第k步系统的实际输出为y_d(k)，计算网络的学习指标函数即误差平方和函数，得到误差平方和E(k)；

子分步骤C2c，根据梯度下降法，分别计算误差平方和E(k)对权值的偏导数，使误差沿梯度方向下降，即令偏导数为0，得到Elman网络的学习算法；

子分步骤C2d：使用以下参数为Elman网络的学习算法赋阈值；以及

子分步骤C2e：以平均云量序列C中的云量数据作为初始值输入数据，让神经网络模型对时间序列进行训练学习，当达到设定的阈值标准时，即存在确定的训练模型权值，由此权值得到T维结果输出，即获得高频随机时间序列预测结果

4.根据权利要求3所述的方法，其特征在于，所述子分步骤C2b中，误差平方和E(k)的计算过程参见下述方程：

其中，y(k)为第k步神经网络系统输入的数值，y_d(k)为第k步神经网络系统输出的数值，n为总步数。

5.根据权利要求3所述的方法，其特征在于，所述子分步骤C2c中，得到的Elman网络的学习算法如下：

式中，w¹，w²和w³分别表示中间层到输出层、输入层到中间层，承接层到中间层的连接权值，η₁，η₂，η₃分别是Elman网络中w¹，w²和w³的学习步长，x_j(k)为第j次迭代输入初始值；代表i，j次在承接层到中间层迭代的权值差；u_q(k-1)为上一步q层中间值；为输出层校验系数；为中间层校验系数；y_a，i(k)为k步神经网络输出数值；g′_i(.)为输出神经元传递函数；f′_j(.)代表中间层神经元的传递函数。

6.根据权利要求3所述的方法，其特征在于，对于试验中FY-2E以旬为时间分辨率的云量数据，预测未来一年的云量预测，所述子分步骤C2d中，Elman网络的学习算法的阈值为：输入单元层m＝144，输出单元层n＝3，隐含层单元数r＝10，最大迭代次数mix2000，训练误差指标θ＝0.15。

7.根据权利要求1所述的方法，其特征在于，所述步骤D包括：

通过叠加计算，将低频趋势时间序列的预测结果与高频趋势时间序列的预测结果整合在一起，得到最终云量预测的时间序列