基于李雅谱诺夫指数和余弦相似度的Kp指数预报方法与流程

本发明涉及一种基于李雅谱诺夫指数和余弦相似度的kp指数预报方法。

背景技术

kp指数即“行星际小时磁情指数”，通过计算磁纬在48度和63度之间的全球13个地磁台站k指数的加权平均值得到。它与3小时时段内地磁扰动有近似对数关系，是衡量近地空间全球磁扰动强度的重要指标之一。kp指数在空间环境监测中可以为地磁暴的发生提供预警，避免空间和地面技术系统遭到灾害性事件的影响。因此，对于kp指数的准确预报有着重要意义。然而由于行星际向磁层电离层能力输入的复杂性，磁层电离层对扰动响应的不确定性以及磁层内部的动力学过程等原因，使得预报kp指数的难度更大。

一些机构发展了短期预报模型，如gehred、takahashi的kp预报模型。利用太阳风数据和现报kp值作为输入参数，借助人工神经网络方法构建的预报模型主要有costello模型、wing等建立的apl模型，boberg模型和bala模型。刘杨等充分利用ace卫星积累的上游行星际条件的数据，以开磁通生成速率函数dφmp/dt和太阳风磁层粘滞作用项n1/2v2为主要输入参数，应用神经网络方法，构建了三个模型，预报三小时时段的kp值。王庚等提出了一种改进的地磁kp指数现报模式，其可以有效识别地磁规则日变化的逐日变化特性，反映地磁扰动的季节效应和地方时效应，从而提升了kp指数现报的准确性。

然而利用人工神经网络方法预报kp指数的方法具有易震荡、收敛速度慢、易陷入局部最优解，同时包括隐含层数目等在内的网络结构难以确定等缺点。

kp指数是非线性和非平稳的时间序列，在一定周期内还可以表现出混沌特性，而对于kp指数的混沌模型的研究则较为少见。时间序列的混沌预测方法有全域法、局域法、加权零阶局域法、加权一阶局域法等。这些方法均属于混沌相空间的几何预测方法，而加权一阶局域法是较为有效和精确的方法。加权一阶局域法往往采用欧氏距离作为相空间向量相似性的度量。但欧氏距离假定相空间向量各分量无关，而相空间中距离较近的各向量间往往并不是完全无关，因此寻找合适的非欧氏距离作为相似性度量，并应用于混沌一阶局域预测中，可以为混沌预测提供更多的参考相似度量，从而在一定程度上提高混沌预测的能力。

利用加权一阶局域法进行kp指数预测的基本思想是将空间轨迹的最后一点作为中心点，把离中心点最近的若干轨迹点作为相关点，找出并根据“历史上情况最相似的情况”估计轨迹下一点的走向，最后从预测出的轨迹点的坐标中分离出预测值。

然而加权一阶局域法使用邻近点的所有分量进行加权拟合，将同一邻近点的所有分量看成“平等的”，无形中引入了误差。如果仅用延迟矢量的第一个分量进行一阶拟合，而不考虑其他分量，则当时间序列的混沌性较强时，此方法可以保证精度。但当混沌性较弱时，其它延迟分量与预测值的相关性仍然很大，如果舍弃，则会导致预测精度降低。可见，只有让权值体现出邻近点的各分量对预测的贡献，才能选出适当的拟合参考点，从而提高预测精度。

另外，对于混沌时间序列而言，序列值之间的相关性是按照时间间隔呈lyapunov指数衰减的。所以，各维分量对预测贡献大小应当由时间延迟和lyapunov指数共同决定。时间越靠前的分量，其预测贡献应该越大，且各分量所加权值应与序列的lyapunov指数成正比。更为重要的是，在进行混沌分析和预测时，通常采用欧氏距离作为衡量相空间轨迹之间差异的度量，而欧式距离在一些时间序列中难以全面反映和表征空间向量的相关性。

因此，利用混沌预测方法对kp指数进行预测是一种新的方法，而在利用加权一阶局域法这种较为有效和精确方法时，需要解决的问题是：①顾及时间序列的lyapunov指数衰减性，改进邻近点权值的设定方法，以提高各分量对预测的贡献；②欧氏距离是衡量相空间轨迹之间差异的基本度量，但欧式距离在一些时间序列中难以全面反映和表征空间向量的相关性。

技术实现要素：

本发明的目的在于提出一种基于李雅谱诺夫指数和余弦相似度的kp指数预报方法，以提高kp指数预测的简便性、可操作性和精度。

本发明为了实现上述目的，采用如下技术方案：

基于李雅谱诺夫指数和余弦相似度的kp指数预报方法，包括如下步骤：

s1.构建kp指数时间序列

kp指数数据每隔3小时预报一次，并且在1-3小时之内的数据是相同的，因此，将相同的数据滤除，得到每3小时为间隔的kp指数时间序列x＝{xi,i＝1,2,…,n}；

其中，n为kp指数时间序列的数据个数；

s2.kp指数时间序列的混沌特性分析

s2.1.计算时间延迟τ；

其中，时间延迟τ表示相空间中相矢量的各分量在时间序列中时间间隔的大小；

s2.2.计算嵌入维数m；

s2.3.计算最大李雅谱诺夫指数λ；

其中，李雅谱诺夫指数是定量描述混沌运动对初值条件的敏感性的吸引子特征量；

将kp指数时间序列重构相空间x＝{xj}，xj为第j个m维相空间中的相点：

xj＝[xj,xj+τ,…,xj+(m-1)τ],j＝1,2,…,m；

其中，m＝n-(m-1)τ，为m维相空间中嵌入点数目；

最大李雅谱诺夫指数可定义为：

其中，dj＝||xj+1-xj||∞，∞为范数，dj表示相空间xj+1与xj这两个相点之间的范数；

s3.kp指数时间序列的预测

基于李雅谱诺夫指数和余弦相似度的预报方法来预报kp指数：

以相空间轨迹的最后一点作为中心点，把与中心点的修正后余弦相似度最小的若干轨迹点作为相关点，以xj+1＝a+bxj来拟合中心点周围的小邻域；

其中，a、b表示拟合参数；

将分量的权值设定成与时间序列的时间和李雅谱诺夫指数成正比，将分量对预测的贡献通过加权的方法体现在余弦相似度这个度量中，则修正后的余弦相似度公式为：

依据修正后的公式计算出余弦相似度，用余弦相似度来计算每个邻近点对应的权值；

以kp指数时间序列重构的相空间x＝{xj}轨迹的最后一点作为中心点；

把离中心点最近的若干轨迹点作为相关点，然后对这些相关点进行拟合，预测出轨迹点的走向；进行一阶加权局域线性拟合，得到系数a，b，进行预测，具体过程如下：

设中心点周围的邻域包括k个相点，分别为t1,t2,…,tk，则xt+1＝a+bxt表示为：

用最小二乘法求出a和b，再通过xm＝a+bxm-1得到相空间中轨迹的趋势，从而能够从xm中分离出时间序列的预测值。

优选地，所述步骤s2.1中，利用c-c法计算时间延迟τ，具体如下：

m维相空间中的相点的关联函数为：

式中，r为计算中所取的搜索半径；θ为heaviside函数，其表达式为：

y为任意实数；

s(m,n,r,τ)反映了kp指数时间序列x＝{xi,i＝1,2,…,n}的自相关特性；

式中，表示关联函数s(m,n,r,τ)的均值；

式中，δs(m,n,r,τ)＝max{s(m,n,r,τ)}-min{s(m,n,r,τ)}，为s(m,n,r,τ)的最大偏差；

其中，r＝lσ/2，σ为kp指数时间序列的标准差；

根据m和l的取值，寻找第一个零点或的第一个局部极小值点即时间延迟τ。

优选地，所述步骤s2.2中，采用g-p算法计算嵌入维数m，具体如下：

根据公式(4)计算关联函数s(m,n,r,τ)，当嵌入维数m分别取值为m＝2,3,…10时，搜索半径r取值为lσ/2，计算lns(m,n,r,τ)/lnr；

随着嵌入维数m的增加，lns(m,n,r,τ)/lnr在一定误差范围内不变为止，此时得到的嵌入维数m，即为kp指数时间序列x＝{xi,i＝1,2,…,n}的嵌入维数。

本发明具有如下优点：

(1)简便性和可操作性：

神经网络存在易震荡、收敛速度慢、易陷入局部最优解，同时包括隐含层数目等在内的网络结构难以确定等缺点，本发明每步都有相应公式，不需人为设置参数。

(2)中、长期预测精度更高：

对kp指数的2天和3天的预测精度，比神经网络模型和一阶局域预测法精度更高。

附图说明

图1为本发明中基于李雅谱诺夫指数和余弦相似度的kp指数预报方法的流程图；

图2为kp指数时间序列示意图；

图3为本发明中c-c法计算时间延迟示意图；

图4为本发明中时间序列的关联积分示意图；

图5为本发明中时间序列的关联空间维数示意图；

图6为现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预报kp指数1天的结果对比示意图；

图7为现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预报kp指数2天的结果对比示意图；

图8为现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预报kp指数3天的结果对比示意图；

图9为现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预报kp指数的中误差比较示意图；

图10为现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预报kp指数的相对中误差比较示意图。

具体实施方式

本发明主要解决两个问题：

①改进邻近点权值的设定方法：

时间越靠前的分量，各分量的权值与序列的李雅谱诺夫指数成正比；

②选择新的度量：

以全面反映和表征空间向量的相关性，提高寻找离中心点最近的若干相关点的准确性。

本发明解决上述问题的具体思路是：

通过研究发现，地球磁场相关的物理参数的时间序列值之间存在一定的角度相关性。

因此，本发明实施例引入余弦相似度这个度量，用向量空间中两个向量夹角的余弦值作为衡量两个个体间相关性的大小。

原欧氏距离为：dj＝||xj+1-xj||∞∞为范数。

其中，kp指数时间序列x＝{xi,i＝1,2,…,n}重构的相空间x＝{xj}。

xj为第j个m维相空间中的相点：xj＝[xj,xj+τ,…,xj+(m-1)τ],j＝1,2,…,m。

m为嵌入维数，τ为时间延迟，m＝n-(m-1)τ，为m维相空间中嵌入点数目。

原余弦相似度为：

将分量的权值设定成与时间序列的时间和李雅谱诺夫指数成正比，将分量对预测的贡献通过加权的方法体现在余弦相似度这个度量中，则修正后的余弦相似度公式为：

λ为最大李雅谱诺夫指数。

依据修正后的公式计算出余弦相似度，用余弦相似度来计算每个邻近点对应的权值。进行一阶加权局域线性拟合，得到系数a，b，从而进行预测。

下面结合附图以及具体实施方式对本发明作进一步详细说明：

本发明实施例以2018年第1天到第59天的kp指数为例，分别利用2018年第1天到第58天的数据预测第59天的kp指数，以第1天到第57天的数据预测第58天和59天的kp指数，第1天到第56天的数据预测第57、58和第59天的kp指数。

结合图1所示，基于李雅谱诺夫指数和余弦相似度的kp指数预报方法，包括如下步骤：

s1.构建kp指数时间序列

kp指数数据每隔3小时预报一次，并且在1-3小时之内的数据是相同的，因此，将相同的数据滤除，得到每3小时为间隔的kp指数时间序列，如图2所示。

s2.kp指数时间序列的混沌特性分析

s2.1.计算时间延迟τ，时间延迟τ是计算其它参数的基础。

时间延迟τ表示相空间中相矢量的各分量在时间序列中时间间隔的大小；由于时间延迟表征时间序列中数据之间的时间关联性，因此，它是相空间重构的重要参数。

计算时间延迟τ的方法主要有自相关法、平均位移法、复自相关法、去偏自相关法、互信息法和c-c法等。本发明实施例选用c-c法计算时间延迟τ，具体如下：

对于kp指数时间序列x＝{xi,i＝1,2,…,n},n为kp指数时间序列的数据个数。将kp指数时间序列重构相空间x＝{xj}，xj为第j个m维相空间中的相点：

xj＝[xj,xj+τ,…,xj+(m-1)τ],j＝1,2,…,m(1)

m维相空间中的相点的关联函数为：

式中，r为计算中所取的搜索半径；θ为heaviside函数，其表达式为：

y为任意实数。

s(m,n,r,τ)反映了kp指数时间序列x＝{xi,i＝1,2,…,n}的自相关特性。

式中，表示关联函数s(m,n,r,τ)的均值；

δs(m,n,r,τ)＝max{s(m,n,r,τ)}-min{s(m,n,r,τ)}，为s(m,n,r,τ)的最大偏差；

其中，r＝lσ/2，σ为kp指数时间序列的标准差；

根据嵌入维数m和l的取值，寻找第一个零点或的第一个局部极小值点即最优的时间延迟τ，如图3所示。由图3可以看出，本发明实施例中的第一个局部极小值点为8，即为最优时间延迟，由此可知，最优时间延迟的时间为8*3＝24，即1天。

s2.2.计算嵌入维数m。

计算嵌入维数m的方法主要有伪最近临域法、真实矢量法、cao方法和饱和关联维数(g-p)法等。其中g-p法是最常用的计算嵌入维数的方法。

本发明实施例采用g-p算法计算嵌入维数m，其计算步骤有：

根据公式(2)，计算关联函数s(m,n,r,τ)，当嵌入维数m分别取值为m＝2,3,…10时，搜索半径r取值为lσ/2，计算lns(m,n,r,τ)/lnr。

随着m的增加，lns(m,n,r,τ)/lnr在一定误差范围内不变为止，此时得到的嵌入维数m，即为kp指数时间序列x＝{xi,i＝1,2,…,n}的嵌入维数。

时间序列的关联积分和关联空间维数分别如图4和图5所示。

由图4可以看出，当m＝5之后，随着m增加，lns(m,n,r,τ)/lnr不在变化。由图5也可以看出，当m＝5之后，增加嵌入维数，相应的lns(m,n,r,τ)/lnr不再随m的增加，而在一定误差范围内在1.5附近保持不变。由图4和图5可以看出，时间序列的嵌入维数m＝5。

s2.3.计算最大李雅谱诺夫指数λ

李雅谱诺夫指数是定量描述混沌运动对初值条件的敏感性的吸引子特征量。

最大李雅谱诺夫指数可定义为：

其中，dj＝||xj+1-xj||∞，∞为范数，dj表示相空间xj+1与xj这两个相点之间的范数。

若系统的最大李雅谱诺夫指数为正，则系统是混沌的，其值越大，混沌特性越强，对初始值的敏感性也越强，反之则敏感性越弱。最大李雅谱诺夫指数的倒数可作为系统的可预报尺度，即只有在该预报尺度范围内，才能对系统作高精度的预测。

kp时间序列的最大lyapunov指数为λ＝0.0417，因此，预报点的个数为

s3.kp指数时间序列的预测

基于李雅谱诺夫指数和余弦相似度的预报方法来预报kp指数：

以相空间轨迹的最后一点作为中心点，把与中心点的修正后余弦相似度(公式如上)最小的若干轨迹点作为相关点，以xj+1＝a+bxj来拟合中心点周围的小邻域。

其中，a、b表示拟合参数。

设中心点周围的邻域包括k个相点，分别为t1,t2,…,tk，则xt+1＝a+bxt表示为：

用最小二乘法求出a和b，再通过xm＝a+bxm-1得到相空间中轨迹的趋势，从而能够从xm中分离出时间序列的预测值。

根据最大李雅普诺夫指数的倒数，计算的可预报尺度为24个点，每天8个点，24个点即3天的数据。为了比较小于可预报尺度的预报精度，本发明除了比较预报24个点(3天)的精度，还比较了预报8个点(1天)和16个点(2天)的精度，如图6～图8所示。

当以中误差为检核指标时，现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预测精度如图9所示。从图9中看出：

当预测1天时，神经网络模型法预测的中误差最小，约为1.37，加权一阶局域法和李雅谱诺夫指数和余弦相似度法中误差相差很小，都在1.48左右，此时神经网络模型法预测精度最高。由此可以看出，当预测时间较短时，神经网络模型有一定的优势，比加权一阶局域法和李雅谱诺夫指数和余弦相似度法的预测精度略高。当预测2天时，神经网络模型精度最差，约为1.4左右，加权一阶局域法为1.27左右，而李雅谱诺夫指数和余弦相似度法约为1.23左右，加权一阶局域法比神经网络模型的中误差小0.1左右，精度提高约11％，而李雅谱诺夫指数和余弦相似度法又比加权一阶局域法略高。而当预测3天时，加权一阶局域法比神经网络模型的中误差小0.1左右，精度提高约10％，而李雅谱诺夫指数和余弦相似度法的中误差又比加权一阶局域法小0.1左右，精度提高约7％，比神经网络模型精度提高约17％。由此可见，当预测时间为2天和3天时，李雅谱诺夫指数和余弦相似度法的预测精度最高，预测3天时，比加权一阶局域精度提高7％，比神经网络模型精度提高17％。

现有技术中神经网络模型、加权一阶局域法以及本发明中的李雅普诺夫指数和余弦相似度法预测的相对中误差如图10所示。由图10可以看出：

当预测1天时，神经网络模型法预测的相对中误差最小，约为0.5，加权一阶局域法和李雅谱诺夫指数和余弦相似度法相对中误差相差很小，都在0.54左右，此时神经网络模型法预测精度最高。由此可以看出，当预测时间较短时，神经网络模型有一定的优势，比加权一阶局域法和李雅谱诺夫指数和余弦相似度法的预测精度略高。当预测2天时，神经网络模型精度最差，约为0.62左右，加权一阶局域法为0.56左右，而李雅谱诺夫指数和余弦相似度法约为0.55左右，加权一阶局域法比神经网络模型的相对中误差小0.06左右，精度提高约10％，而李雅谱诺夫指数和余弦相似度法又比加权一阶局域法略高。而当预测3天时，加权一阶局域法比神经网络模型的相对中误差小0.08左右，精度提高约12％，而李雅谱诺夫指数和余弦相似度法的相对中误差又比加权一阶局域法小0.05左右，精度提高约8％，比神经网络模型精度提高约20％。由此可见，当预测时间为2天和3天时，李雅谱诺夫指数和余弦相似度法的预测精度最高，其中，预测3天时，比加权一阶局域精度提高8％，比神经网络模型精度提高20％。具体精度指标如表1所示。

表1三种方法不同预测天数的精度比较

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。