一种基于经验模态分解的短期风速预测方法与流程

本发明涉及风速预测
技术领域：
，具体涉及一种基于经验模态分解的短期风速预测方法。
背景技术：
：人类文明的快速发展离不开能源的开发，现有的化石能源不仅有限而且会对我们的生态环境造成很大的破坏，为了社会的可持续发展，发展风能尤其重要。2016年中国风电新增装机容量占到了全球市场份额的42.7％，但是短期风力的随机性、非平稳性制约着我们对风能的利用水平。因此，对短期风速预测的精度对于我们发展风能来说是非常重要的。风速是一个非平稳的时间序列，针对于短期风速预测，传统的时间序列分析方法是通过直接差分的方法对风速进行平稳化处理，差分采用的阶数随着风速信号的复杂性而上升，虽然最终能够达到平稳效果，但是每次的差分都会损失风速序列的信息，从而影响风速预测的精度。技术实现要素：针对上述现有短期风速预测中信息利用率较低且精度不足的缺点，本发明提出一种基于经验模态分解的短期风速预测方法。为了实现上述任务，本发明采用以下技术方案：一种基于经验模态分解的短期风速预测方法，包括以下步骤：步骤1，对原始的风速序列进行预处理对原始的风速序列中的缺失值、异常值进行预处理，采用的策略为：对于异常值和缺失值，用该数值前后的风速平均值代替；步骤2，分解预处理后的风速序列，得到风速子序列步骤2.1，对于预处理后的风速序列e(t)，找出该风速序列上的所有极大值点，通过三次样条插值函数形成极大值点的包络线emax(t)，采用同样的方法找出该风速序列上所有的极小值点形成极小值点包络线emin(t)；极小值点包络线和极大值点包络线的均值记为m1，用风速序列e(t)减去m1，得到h1(t)：h1(t)＝e(t)-m1式2步骤2.2，将h1(t)看成是一个新的信号序列，计算其系数dk；其中，第k次计算h1(t)计算系数dk的公式如下：其中，t为风速数据的序号，t表示信号序列h1(t)中数据个数，h1k-1(t)、h1k(t)为第k-1次、第k次计算dk系数时的h1(t)；判断系数dk的值是否在0.1～0.2之间，如果不在则将h1(t)作为e(t)并重复步骤2.1，之后将k的值加1然后执行步骤2.2；步骤2.3，如经过k次迭代后满足dk在0.1～0.2之间，则此时得到一个经验模态函数：上式中，为第k次迭代时的h1(t)，c1(t)表示第1个经验模态函数；计算经验模态函数的余项：r1(t)＝e(t)-c1(t)式5步骤2.4，判断r1(t)函数是否为单调函数或者为一个常数，若是，则说明不能再分解，分解结束；否则，令e(t)＝r1(t)，重复步骤2.1至步骤2.3，继续得到新的经验模态函数；，则所有经验模态函数以及最后一次分解所得到的经验模态函数的余项即为风速子序列；则风速序列e(t)表示为：上式中，ci(t)表示第i个经验模态函数，rn(t)表示最后一次分解得到的经验模态函数的余项，n表示分解次数，则风速子序列有n+1个；步骤3，确定风速子序列进行arma建模所需的阶数p0、q0；步骤4，确定最佳阶数并进行arma建模取迭代区间为p∈[p0-5,p0+5]，q∈[q0-5,q0+5]，将阶数p、q组合代入matlab中的aic(p,q)信息系数计算函数，计算p、q组合对应的信息系数aic；取最小的信息系数aic对应的阶数p、q作为最佳阶数，记为pbest,qbest，利用该参数进行arma建模；步骤5，利用步骤3和步骤4相同的方法，得到所有风速子序列的arma模型；每一个arma模型的预测结果分别记为a1,a2,...,an+1，将结果叠加即为风速序列最终的预测结果，表示为：进一步地，步骤3所述的确定风速子序列进行arma建模所需的阶数，包括：做出风速子序列的采样自相关系数图和采样偏相关系数图，找到采样自相关系数图、采样偏相关系数图在收敛至置信区间时对应的滞后数作为所需阶数。进一步地，步骤4中，若p0-5或q0-5的值小于0，则将其置为0。本发明与现有技术相比，具有以下技术特点：为了克服差分造成的信息损失对预测结果的影响，本发明采用了经验模态分解的方法将非平稳复杂的原始的风速先进行经验模态分解，将复杂的风速分解为多个相对平稳的子风速序列，选择平稳的风速信号进行时间序列建模，这样一来既能够减少由于差分阶数过大带来的信息损失，又能达到预测的需求，从而达到一个较高的预测水平，有效地解决了由于高阶差分平稳化处理造成的预测精度不足的问题。附图说明图1为本发明方法的流程示意图；图2的(a)为风速序列中存在异常值的示意图，(b)为风速序列中存在缺失值的示意图，单位：m/s；图3为预处理后的风速序列的系数图；图4的(a)为采样自相关系数图，(b)为采样偏相关系数图；图5为本发明方法的风速预测效果图。具体实施方式下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。本实施例公开了一种基于经验模态分解的短期风速预测方法，如图1所示，具体包括以下步骤：步骤1，对原始的风速序列进行预处理对原始的风速序列中的缺失值、异常值进行预处理，采用的策略为：对于异常值和缺失值，用该数值前后的风速平均值代替，以保证数据的完整性。其中，异常值是指与其他数值相比明显偏大或偏小的数据，如将与平均数据偏差n倍以上的数据作为异常值。在本实施例中，数据来源于美国林肯首都机场2017年12月1日到2017年12月17日逐小时的共400小时的风速数据作为对象，预测利用前300小时测量值，预测逐小时的风速数据。如图2所示，在原始的风速序列中，风速值为999.99为异常值，风速值为nan的数据表示缺失值，对其前后两个小时数据取平均值替代之即：(3.1+5.1)/2＝4.5替代999.99，nan用(5.7+5.7)/2＝5.7替代，处理结果如图3所示。最终处理的结果是完整的、可信的。步骤2，分解预处理后的风速序列，得到风速子序列步骤2.1，对于预处理后的风速序列e(t)，找出该风速序列上的所有极大值点，通过三次样条插值函数形成极大值点的包络线emax(t)，采用同样的方法找出该风速序列上所有的极小值点形成极小值点包络线emin(t)；极小值点包络线和极大值点包络线的均值记为m1，用风速序列e(t)减去m1，得到h1(t)：h1(t)＝e(t)-m1式2步骤2.2，将h1(t)看成是一个新的信号序列，计算其系数dk；其中，第k次计算h1(t)计算系数dk的公式如下：其中，t为风速数据的序号，t表示信号序列h1(t)中数据个数，h1k-1(t)、h1k(t)为第k-1次、第k次计算dk系数时的h1(t)。判断系数dk的值是否在0.1～0.2之间，如果不在则将h1(t)作为e(t)并重复步骤2.1，之后将k的值加1然后执行步骤2.2。步骤2.3，如经过k次迭代后满足dk在0.1～0.2之间，则此时得到一个经验模态函数：上式中，为第k次迭代时的h1(t)，c1(t)表示第1个经验模态函数。计算经验模态函数的余项：r1(t)＝e(t)-c1(t)式5步骤2.4，判断r1(t)函数是否为单调函数或者为一个常数，若是，则说明不能再分解，分解结束，否则，令e(t)＝r1(t)，重复步骤2.1至步骤2.3，继续得到新的经验模态函数；则所有经验模态函数以及最后一次分解所得到的经验模态函数的余项即为风速分量，也即风速子序列；则风速序列e(t)可表示为：上式中，ci(t)表示第i个经验模态函数，rn(t)表示最后一次分解得到的经验模态函数的余项，n表示分解次数，则风速子序列有n+1个。本实施例中，风速序列最终分解得到8个风速子序列。步骤3，确定风速子序列进行arma建模所需的阶数这里以本实施例得到的第2个风速子序列为例，做出风速子序列的采样自相关系数图(sampleautocorrection)和采样偏相关系数图(samplepartialautocorrelation)，如图4所示，采样自相关系数图和采样偏相关系数图是收敛的，采样自相关系数图在滞后数lag＝10时收敛至置信区间，而采样偏相关系数图在滞后数lag＝2时收敛至置信区间，因此该风速子序列进行经验arma建模所需的阶数p、q的初步取值分别为10、2，记为p0＝10，q0＝2。该步骤采用经验法进行预估，最终预估得到的最佳阶数p、q的取值将作为步骤4的重要参考。步骤4，确定最佳阶数并进行arma建模通常来讲，步骤3基于经验法得到的风速子序列进行经验arma建模所需的阶数p、q不一定是最佳的，需要进行进一步处理。取迭代区间为p∈[p0-5,p0+5]，q∈[q0-5,q0+5]，若p0-5或q0-5的值小于0，则将其置为0；由此产生最多10*10＝100个阶数p、q组合，代入matlab中的aic(p,q)信息系数计算函数，计算p、q组合对应的信息系数aic；取最小的信息系数aic对应的阶数p、q作为最佳阶数，记为pbest,qbest，则该风速子序列进行arma建模所需的最佳阶数为pbest,qbest，即arma(pbest,qbest)，利用该参数进行arma建模。本实施例中，第2个风速子序列的迭代区间取p∈[5,15],q∈[0,7]，共产生10*7＝70个阶数组合代入matlab中的aic(p,q)计算信息系数，最小的信息系数对应的p、q为8、2，则pbest＝8,qbest＝2，因此这个风速子序列模型可以取arma(8,2)。步骤5，利用步骤3和步骤4相同的方法，得到所有风速子序列的arma模型；每一个arma模型的预测结果分别记为a1,a2,...,an+1，将结果叠加即为风速序列最终的预测结果，表示为：本实施例最终的预测结果如图5所示。利用本发明方法和传统的arma建模的参数比较结果如表1所示，可以看出本发明方法有效地提提升了预测精度。表1本发明方法和传统的arma建模参数比较模型\指标msemape传统arma建模0.700114.8204本发明方法0.49399.8496上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。当前第1页12