基于改进指数平滑法的钟差预报方法与流程

本发明属于原子钟数据处理
技术领域：
，尤其涉及一种卫星星载原子钟的钟差预报方法。
背景技术：
：高精度的时间同步技术直接影响着卫星的导航、定位和授时精度，不仅对卫星的正常运转意义重大，而且对如分布式组网雷达系统、多基地雷达系统等分布式武器系统的预警、制导和定位等性能也有着重要的影响。时间同步是卫星、分布式武器系统正常运行的基础，钟差预报是卫星、分布式武器系统在无法进行时间比对的情况下保持时间同步的一种有效手段。为了在卫星无法与地面站点进行比对的时段中保持星载钟源与系统钟源的时间同步，一般采用钟差预报的方法对二者钟差进行预报。此外，分布式武器系统在战时条件下，部分站点可能会因干扰、故障、部分损伤等原因导致该站点无法与其他站点进行比对来获得站间钟差，在这种情况下，也采用钟差预报的方法对该站点原子钟时间与系统时间的钟差进行预报，来获得该站点与系统的时间同步。卫星无法与地面站点进行比对的时长一般为几小时，虽然也存在因钟源故障导致长达几天无法比对的情况，但总体而言，钟差预报一般为中短期预报，钟差预报的时长一般为若干小时。准确的钟差预报是后续钟差校准的前提，也是在无法进行时间比对的情况下保证站点与系统时间同步的基础，如何在分布式系统的部分站点无法进行比对的情况下保持站点与系统的高精度时间同步一直是本领域的研究热点。技术实现要素：本发明的目的在于提供一种基于改进指数平滑法的钟差预报方法。为了实现上述目的，本发明采取如下的技术解决方案：基于改进指数平滑法的钟差预报方法，包括以下步骤：数据预处理步骤、获取某一时间段的原始钟差序列，对原始钟差序列进行差分处理，得到钟差差分序列；参数初始化步骤、将钟差差分序列分为若干个子序列，每一个子序列均分为学习段序列和预报段序列，对指数平滑法的初始值进行初始化；平滑系数优化步骤、对钟差差分序列各子序列的平滑系数进行优化，确定各子序列的最优平滑系数：在(0，1)范围内均分m份，分别得到m-1个平滑系数的分值，平滑系数的第m个分值＝＝m/m，m＝1，2，...，m-1，将平滑系数的每一个分值作为指数平滑法的平滑系数，采用指数平滑法对各个子序列进行钟差预报，对于某一子序列，计算与该子序列的平滑系数的每一个分值对应的预报误差，从m-1个预报误差中选出最小值，将与该最小值对应的平滑系数的分值作为该子序列的最优平滑系数；预报序列获取步骤、基于各子序列的最优平滑系数，采用指数平滑法获取各子序列的预报序列；误差预报步骤、将预报序列与钟差差分序列作差，得到误差序列，利用灰色预报模型计算误差预报值；误差预报值式中的为灰色预报模型的拟合参数，k为子序列的学习段长度，y(0)(1)为误差序列的第1个元素；预报序列融合步骤、将预报序列获取步骤得到的结果与误差预报步骤得到的结果相加获得融合预报序列；数据后处理步骤、将融合预报序列进行差分逆运算，得到最终的钟差预报序列。进一步的，所述指数平滑法为二阶指数平滑法或三阶指数平滑法。进一步的，在获得融合预报序列之后，各子序列利用其前一子序列的融合预报序列对自身的学习序列进行更新，然后根据滑动更新后的序列，重复进行分段钟差预报和滑动更新，直到到达预报时长，将得到的各段融合预报序列拼接得到完整的钟差预报序列，然后进行数据后处理。由以上技术方案可知，本发明针对传统指数平滑法进行钟差预报时，平滑系数固定带来的预报误差累积较大的缺点，引入基于滑动窗的动态平滑系数进行预报，同时利用灰色预报模型进行误差学习，从而提高钟差预报精度。附图说明图1为实施例1的流程图；图2为实施例1钟差差分序列的分段示意图；图3为实施例2的流程图；图4为本发明实施例2钟差差分序列的分段示意图；图5a至图5c分别为采用常规的des法、tes法和gm法进行钟差预报的性能曲线图；图6a和图6b分别为采用本发明实施例1的方法进行钟差预报的性能曲线图；图7a和图7b分别为采用常规的des法和tes法与本发明实施例2的方法以及常规的gm法来进行钟差预报的性能对比图。下面结合附图和各实施例对本发明进一步详细说明。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下获得的所有其他实施例，都属于本发明保护的范围。实施例1图1为本实施例的流程图，如图1所示，本实施例钟差预报方法包括以下步骤：数据预处理步骤、获取某一时间段的原始钟差序列{z1，z2，...，zn，zn+1}，其中，z1、z2...、zn、zn+1分别表示该时间段内不同时刻的原始钟差值，n为原始钟差序列的时长，对原始钟差序列进行差分处理，得到钟差差分序列{x1，...，xi，...，xn}，xi＝zi+1-zi，i＝1，...，n；由于原始钟差序列可能存在部分钟差数据绝对值很大，为了减少指数运算的复杂度，本发明对原始钟差序列采用差分的方法进行预处理，再基于差分序列进行后续的步骤；参数初始化步骤、将钟差差分序列分为若干个子序列，每一个子序列均分为学习段序列和预报段序列(图2)，每一子序列的学习段序列构成一个滑动窗(slidingwindow，简称sw)，对指数平滑法的初始值(s1′、s1″、s1″′)进行初始化，令指数平滑法的初始值等于钟差差分序列中第1个元素x1，即令s1′＝s1″＝s1″′＝x1；平滑系数优化步骤、为了减小误差累积，本发明要对各子序列的平滑系数(α1，α2，...，αk)进行优化，以确定各子序列的最优平滑系数下面对子序列的平滑系数的优化过程作进一步的说明：平滑系数α为经验值，其取值范围为0＜α＜1，在(0，1)范围内均分m份，分别得到m-1个平滑系数的分值，平滑系数的第m个分值＝m/m，m＝1，2，...，m-1，m为整数，基于平滑系数的每一个分值，分别采用指数平滑法对各个子序列进行钟差预报，对于某一子序列，计算与该子序列的平滑系数的每一个分值对应的预报误差，从m-1个预报误差中选出最小值，将与该最小值对应的分值作为该子序列的最优平滑系数；以本实施例中第1个子序列part1的平滑系数为例，将(0，1)区间均分为m＝100份，则第1个子序列的平滑系数α1的第m个分值α1(m)＝m/100，m＝1，2，...，99，即第1个子序列的平滑系数α1的初值α1(1)(第1个分值)为0.01，终值α1(99)(第99个分值)为0.99；以平滑系数的每一个分值α1(m)作为指数平滑法中的平滑系数分别进行第1个子序列的钟差预报，计算得到与各分值对应的预报误差rmseα，从99个预报误差中找到最小值，与该预报误差最小值对应的分值即为第1个子序列的最优平滑系数预报误差式中的zjj为原始钟差序列中的第jj个钟差值，为第jj个预报钟差值，n为子序列的序列长度；常用的指数平滑法(exponentialsmoothing，简称es)包括一阶指数平滑法(ses)、二阶指数平滑法(des)和三阶指数平滑法(tes)，下面分别对这三种方法进行说明：对于钟差差分序列{x1，...，xi，...，xn}的各子序列，当采用一阶指数平滑法对各子序列进行预报时，预报值式中的α为平滑系数，本步骤中的平滑系数为平滑系数的各分值，j表示求和序列变量，t为预报时刻，xt-j为第t-j个钟差差分值，s1′为一阶指数平滑法的初始值；一阶指数平滑法的预报值是二阶指数平滑法和三阶指数平滑法预报的基础；指数平滑法预报从t+1时刻开始，基于前面t个钟差差分值平滑外推下一时刻的预报值当采用二阶指数平滑法进行预报时，预报值ft+ξ＝at+btξ，式中的st′为t时刻一阶指数平滑法的钟差预报值，st″为t时刻二阶指数平滑法的钟差预报值，ξ为预报长度；二阶指数平滑法的初始值s1″已知，预报值st″可根据st″＝αst′+(1-α)st-1″依次递推得到；当采用三阶指数平滑法进行预报时，预报值式中的st′为t时刻一阶指数平滑法的钟差预报值，st″为t时刻二阶指数平滑法的钟差预报值，st″′为t时刻三阶指数平滑法的钟差预报值，η为预报长度；预报序列获取步骤、基于各子序列的最优平滑系数采用指数平滑法获取各子序列的预报序列，在本步骤中采用指数平滑法获取预报序列时，平滑系数为平滑系数优化步骤所确定的对应子序列的最优平滑系数，采用指数平滑法得到的预报值组成预报序列；误差预报步骤、将预报序列与钟差差分序列作差，得到误差序列，然后利用灰色预报模型(greymodel，gm)计算误差预报值；采用灰色预报模型对误差序列进行学习预报：式中的a为灰色预报模型的灰色系数，反映序列之间的发展态势，b为灰色预报模型的灰色作用量，反映序列的变化关系，e为自然常数，y(1)为误差序列的一阶累加序列；根据灰色预报模型计算p时刻的误差预报值，式中的为灰色预报模型的拟合参数，k为子序列的学习段长度，即用k个误差序列进行学习，得到拟合参数，预报后续p个误差序列，y(0)(1)为误差序列的第1个元素；灰色预报模型的拟合参数可采用最小二乘估计法获得：其中，y(1)(1)、y(1)(2)...分别表示误差序列的一阶累加序列的第1个、第2元素，y(0)(2)、y(0)(3)...分别表示误差序列的第2、第3个元素；采用灰色预报模型进行预报前，先检查待误差序列中各元素的正负号是否一致，若一致则直接进行预报，若不一致，则给每个元素加一个常数，该常数与误差序列中绝对值最大的元素的符号一致，且该常数的绝对值大于误差序列中最大元素绝对值，通过灰色预报模型获得预报值后减去该常数获得最终灰色预报值；预报序列融合步骤、将预报序列获取步骤得到的结果与误差预报步骤得到的结果相加获得融合预报序列；数据后处理步骤、将上一步得到的融合预报序列进行差分逆运算，得到最终的钟差预报序列。实施例2实施例1的方法用于短期钟差预报，当进行中期钟差预报时，为了改善预报精度，本实施例在实施例1的基础上，在获得预报序列后，利用预报序列对各子序列的学习序列进行更新，然后再进行误差预报，再将预报误差的结果和预报序列进行融合，得到最后的钟差预报结果。图3为实施例2的流程图，如图3所示，本实施例的方法包括以下步骤：数据预处理步骤、获取某一时间段的原始钟差序列{z1，z2，...，zn，zn+1}，对原始钟差序列进行差分处理，得到钟差差分序列{x1，...，xi，...，xn}，xi＝zi+1-zi，i＝1，...，n；参数初始化步骤、将钟差差分序列分为若干个子序列，如图4所示，钟差差分序列被分为k个子序列(part1，part2，...，partk)，每一个子序列均包括学习段序列和预报段序列，图4中的l1、l2...分别表示各子序列的学习段序列的时长，对指数平滑法的初始值及滑动窗长度等参数(s1″′、s1″、s1′、lk)进行初始化，令指数平滑法的初始值等于钟差差分序列中第1个元素x1，即令s1′＝s1″＝s1″′＝x1，lk根据预报时长确定；平滑系数优化步骤、对钟差差分序列的各子序列的平滑系数(α1，α2，...，αk)进行优化，确定各子序列的最优平滑系数在(0，1)范围内均分m份，分别得到m-1个平滑系数的分值，平滑系数的第m个分值＝m/m，m＝1，2，...，m-1，m为整数，以每一个分值作为平滑系数，分别采用指数平滑法对各个子序列进行钟差预报，计算与各子序列的平滑系数的每一个分值对应的预报误差，从m-1个预报误差中选出最小值，将与该最小值对应的分值作为该子序列的最优平滑系数；分段钟差预报步骤、基于各子序列的最优平滑系数采用指数平滑法获取各子序列的预报序列，将预报序列与钟差差分序列作差，得到误差序列，然后利用灰色预报模型计算误差预报值，得到误差预报序列，将指数平滑法获得的各子序列的预报序列和灰色预报模型计算得到误差预报序列相加，获得融合预报序列；滑动更新步骤、各子序列利用其前一子序列的融合预报序列对自身的学习序列进行更新；如图4所示，各子序列利用其前一子序列的融合预报序列中的靠近该子序列的后段数据对自身的学习序列进行更新；完整钟差预报序列获取步骤、根据滑动更新后的序列，重复执行分段钟差预报步骤和滑动更新步骤，直到到达预报时长，将得到的各段融合预报序列拼接得到完整的钟差预报序列；数据后处理步骤、将上一步得到的钟差预报序列进行差分逆运算，得到最终的钟差预报序列。为了验证本发明方法的预报效果，随机选用第1980个和第1981个gps周(2017.12.17-2017.12.30)的钟差数据进行预报，同时采用现有技术的基于灰色预报模型的预报方法进行对比。所选择的星钟情况如表1所示。表1星钟类型星钟编号iirrbpg13(no.1)pg23(no.2)iir-mrbpg5(no.3)pg17(no.4)iifcspg8(no.5)pg24(no.6)方案1分别采用了常规的des法、tes法和gm法进行短期钟差预报，预报时长为一天。采用一天的钟差序列来进行学习，之后获得预报序列，选择rmse作为评价指标，同时列出了预报误差的最大值和最小值，最后统计分析了本方案的均值和标准差，统计结果分别如图5a至图5c和表2所示。从图5a至图5c以及表2的结果可知，des法和tes法短期预报性能大体上在同一精度，des法略优于tes法，gm法的短期预报性能差于本发明的tes法和des法，同时也可得出铷钟的短期预报性能优于铯钟。方案2为采用实施例1的方法进行短期钟差预报，与方案1的常规方法进行对比分析，同样采用一天的钟差序列来进行学习，之后获得预报序列，统计结果分别如图6a、图6b和表2所示。对比图5a、图5b、图5c和图6a、图6b，及表2，可以看出在常规的二阶、三阶指数平滑法中引入灰色预报模型进行误差学习能够一定程度上改善短期指数平滑法的短期预报性能，但由于指数平滑法的短期预报精度已经很高，灰色预报模型学习改善效果不明显。方案3分别采用常规的des法、常规的tes法、实施例2的方法以及常规的gm来进行中期钟差预报，预报时长为一星期。将所选时段中的第一天钟差数据作为des/tes/gm的学习样本，目标预报时长为7天。在滑动窗方面，将预报时长均分为2个子序列，part1的预报时长为3.5天，选择part1预报序列中的最后一天时长数据作为part2的学习样本进行更新，之后进行part2的预报，预报统计性能如图7a、图7b及表2所示。从图7a、图7b和表2分析可知，des的中期预报性能大体上优于tes，加入滑动更新后能够减小预报误差，gm的误差学习能够改善预报性能。与gm算法相比，des与gm预报精度大致在同一水平上，而tes差于gm，在des下，铷钟预报性能优于铯钟，而在tes下则相反。从表2统计分析可知，和单纯des(tes)相比，本发明实施例2的方法能够将预报误差减小35.37％(66.34％)。表2预报性能统计表(单位：ns)以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。当前第1页12