一种利用GPS信号驯服晶振频率的方法、系统及介质与流程

本发明涉及gps驯服晶振信号的频率源误差滤波技术，具体涉及一种利用gps信号驯服晶振频率的方法、系统及介质，用于实现gps信号对晶振频率的高精度驯服。

背景技术：

随着网络技术的快速发展，系统对时间基准的要求也越来越高。电力系统中时间基准对事故分析、提高电能质量、电网调度操作以及电能分时计费有着重要的意义。

目前时间(频率)基准按照应用领域和性能指标可划分为一级频率标准(铯原子频标)、二级频率标准(包括铷原子频标和高稳晶体振荡器)和其它频率标准(高稳晶体振荡器意外的其他晶体振荡器)。虽然一级频标长期稳定度和准确度较高，但是其价格十分昂贵，难以应用在对成本有严格要求的民用场合；二级频标的价格相对较低，其长期稳定度和准确度却大不如前者，故难以应用于高精度场合。

现有技术利用gps锁定二级频标(高稳晶振)原理如附图1所示：首先测量gps接收机输出的秒脉冲信号pps与晶振分频后输出的脉冲信号之间的相对频差，然后对频差信号进行滤波处理，最后生成控制信号，不断调节晶振频率从而输出高准确度和高稳定度的频率输出。该技术在一定程度上解决了上述问题，即利用较低的成本实现高长期稳定度和平均准确度的频率标准。但是在实际应用中由于gps的1pps信号存在随机误差和跳变野值、晶振因老化和温度等特性产生的频率漂移为频率的校准造成了一定的困难。

为了解决以上问题，有的学者提出采用多次采样取平均的方法，对测量的频差数据进行处理，计算简单但是校频时间较长。有学者通过建立晶振信号和gps信号之间的频差数学模型，利用最小二乘法对其误差进行计算和补偿，但该方法计算较为复杂且没有考虑gps信号的野值对校频精度的影响。还有学者提出利用卡尔曼滤波算法来消除频差信号中的噪声，但是现有利用卡尔曼滤波算法来消除频差信号中的噪声的方式没有考虑晶振频率漂移的影响，导致校频精度有限。

技术实现要素：

本发明要解决的技术问题：针对现有技术的上述问题，提供一种利用gps信号驯服晶振频率的方法、系统及介质，本发明基于差分进化算法的新息序列加权扩展卡尔曼滤波算法，能够对gps信号中的跳变野值和随机误差、晶振信号中累积误差和非线性频率漂移进行滤波处理，能够提高gps锁定二级频标的最终校频精度。

为了解决上述技术问题，本发明采用的技术方案为：

一种利用gps信号驯服晶振频率的方法，该方法的步骤包括：

1)根据gps秒信号、晶振信号的误差特性建立两者的频差数学模型；

2)根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差、得到受跳变野值和系统噪声影响较大的滤波模型；

3)引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整以减少其对于频差信号滤波结果的影响；

4)整合新息序列加权到完成过程噪声动态调整后的扩展卡尔曼滤波处理中，以减少跳变野值对于频差信号滤波结果的影响。

可选地，步骤1)的详细步骤包括：

1.1)获取gps秒信号序列y，gps秒信号序列y中的任意k时刻的gps秒信号的世界统一时间utc的函数表达式如式(1)所示；

yk＝k-εk(1)

上式中，yk表示gps秒信号序列y中的k时刻的gps秒信号的世界统一时间utc，εk为k时刻的gps秒信号的误差，且满足ε～n(0,σ²)；

获取晶振信号序列y′，晶振信号序列y′中的任意k时刻的晶振信号的世界统一时间utc的函数表达式如式(2)所示；

y′k＝k+a+bk+ck²(2)

上式中，y′k表示晶振信号序列y′中的k时刻的晶振信号的世界统一时间utc，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；

1.2)根据gps秒信号序列y、晶振信号序列y′做差得到频差序列，且频差序列中任意k时刻的频差信号的函数表达式如下式所示；

xk＝y′k-yk＝a+bk+ck²+εk(3)

上式中，xk表示频差序列中k时刻的频差信号。

可选地，步骤2)中根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理的步骤包括：

2.1)建立频差的状态方程和观测方程的函数表达式如式(4)所示；

上式中，x(k+1)表示k+1时刻的频差状态向量，x(k)表示k时刻的频差状态向量，f和h分别表示状态函数和量测函数，z(k+1)表示k+1时刻的量测值；ω(k)表示k时刻的过程噪声向量，其协方差矩阵为q；γ(k)表示k时刻的量测噪声向量，其协方差矩阵为r；

2.2)令k时刻到k+1时刻的状态转移矩阵f(k+1,k)、k时刻的观测矩阵h(k)的函数表达式如式(5)所示；

上式中，f表示状态函数，表示k时刻频差状态向量的估计，h表示量测函数，x(k)表示k时刻的频差状态向量；

2.3)建立常规扩展卡尔曼滤波模型的预测方程如式(6)～(8)所示、更新方程如式(9)～(10)所示，其中过程噪声如式(11)所示，过程噪声方差如式(12)所示；

式(6)中，表示基于k时刻预测值的k+1时刻频差状态向量预测值，f表示状态函数，表示k时刻频差状态向量的最优估计值；

p(k+1,k)＝f(k+1,k)p(k,k)f^t(k+1,k)+q(k)(7)

式(7)中，p(k+1,k)表示k+1时刻的协方差向量预测值，f(k+1,k)表示k时刻到k+1时刻的状态转移矩阵，p(k,k)表示k时刻的协方差向量预测值，f^t(k+1,k)表示k时刻到k+1时刻的状态转移矩阵的转置，q(k)表示k时刻的过程噪声方差；

k(k+1)＝p(k+1,k)h^t(k+1)[h(k+1)p(k+1,k)h^t(k+1)+r(k+1)]^-1(8)

式(8)中，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，p(k+1,k)表示k+1时刻的协方差向量预测值，h^t(k+1)表示k+1时刻的观测矩阵h(k+1)的转置，p(k+1,k)表示k+1时刻的协方差向量预测值，r(k+1)表示k+1时刻的观测噪声方差向量；

式(9)中，表示k+1时刻频差状态向量的最优估计值，表示k+1时刻频差状态向量的预测值，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，z(k+1)表示k+1时刻的量测值，表示由k+1时刻频差状态向量预测值得到的k+1时刻量测向量预测值；

p(k+1)＝[i-k(k+1)h(k+1)]p(k+1,k)(10)

式(10)中，p(k+1)表示k+1时刻更新的协方差向量估计值，i表示维数与k(k+1)h(k+1)相同的单位矩阵，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，h(k+1)表示k+1时刻的观测矩阵，p(k+1,k)表示k+1时刻的协方差向量预测值；

式(11)中，ω(k+1)表示k+1时刻的过程噪声，表示k+1时刻频差状态向量的最优估计值，z(k+1)表示k+1时刻的量测值；

qk+1＝var(ω(k+1))(12)

式(12)中，qk+1表示k+1时刻的过程噪声方差，ω(k+1)表示k+1时刻的过程噪声，var表示求过程噪声方差的操作。

可选地，步骤3)中引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整时，通过差分进化算法生成任意k时刻的过程噪声方差的步骤包括：

3.1)进行种群初始化，且函数表达式如下式(13)所示：

xi,k(0)＝lk+rand(0,1)*(uk-lk)(13)

式(13)中，xi,k(0)表示种群中第0代的第i条“染色体”的第k个基因；rand(0,1)为0～1之间均匀分布的随机数；uk和lk为搜索的上界和下界，且uk和lk分别取过程噪声方差的最大值qmax和最小值qmin；

3.2)针对种群进行变异操作，且变异操作的函数表达式如下式(14)所示：

xi(g+1)＝xr1(g)+f*[xr2(g)-xr3(g)](14)

式(14)中，xi(g+1)为得到的g+1代变异个体；f为压缩比例因子，取值范围0～1；xr1(g)、xr2(g)和xr3(g)为3个父代；

3.3)针对种群进行交叉操作以保留较优良的变量，且交叉操作采用二项交叉方式，其函数表达式如下式(15)所示：

式(15)中，yi,j(g+1)表示当前个体对应的交叉操作结果，xi,j(g+1)表示变异个体对应的分量，xi,j(g)表示当前个体对应的分量，r为每个变量生成的一个0～1之间均匀分布的随机数；cr为变量的交叉概率；rnd为1～d之间均匀分布的整数，如果r＜cr则接受变异个体对应的分量xi,j(g+1)，否则保留当前个体对应的分量xi,j(g)；

3.4)针对交叉操作得到的结果，采用如下式(16)所示贪婪选择的方式选择最优解并作为k时刻的过程噪声方差；

式(16)中，xi(g+1)为选择结果，yi(g+1)表示交叉操作结果对应的个体，xi(g)表示当前个体，f[yi(g+1)]表示个体yi(g+1)对应的种群适应度值，f[xi(g)]表示个体xi(g)对应的种群适应度值。

可选地，步骤4)整合新息序列加权到完成过程噪声进行动态调整后的扩展卡尔曼滤波处理中后，得到的扩展卡尔曼滤波处理的更新方程的函数表达式如下式(17)～(21)所示；

dk+1＝p(k+1,k)+rk+1(21)

式(17)～(21)中，表示调整后的k+1时刻频差状态向量的最优估计值，表示k+1时刻频差状态向量预测值，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，φk+1表示修正函数，τk+1表示衡量k+1时刻量测误差大小的变量，ek+1表示k+1时刻的量测误差，z(k+1)表示k+1时刻频差向量的量测值，表示k+1时刻量测向量的预测值，τk+1服从χ²分布，当出现野值时τk+1增大、φk+1(·)减小；ck+1为选取的门限值或常数序列，λk+1为矩阵的最大特征值，kk+1表示k+1时刻的卡尔曼滤波增益，dk+1表示k+1时刻衡量量测误差大小的权矩阵，p(k+1,k)表示k+1时刻的协方差预测值，r(k+1)表示k+1时刻的量测噪声方差。

此外，本发明还提供一种利用gps信号驯服晶振频率的系统，包括：

频差模型生成程序单元，用于根据gps秒信号、晶振信号的误差特性建立两者的频差数学模型；

常规滤波程序单元，用于根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差、得到受跳变野值和系统噪声影响较大的滤波模型；

噪声调整程序单元，用于引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整以减少其对于频差信号滤波结果的影响；

滤波整合程序单元，用于整合新息序列加权到完成过程噪声进行动态调整后的扩展卡尔曼滤波处理中。

此外，本发明还提供一种利用gps信号驯服晶振频率的系统，包括计算机设备，该计算机设备被编程或配置以执行所述利用gps信号驯服晶振频率的方法的步骤。

此外，本发明还提供一种利用gps信号驯服晶振频率的系统，包括计算机设备，该计算机设备的存储器上存储有被编程或配置以执行所述利用gps信号驯服晶振频率的方法的计算机程序。

此外，本发明还提供一种带有gps的智能终端，包括gps模块、微处理器和存储器，该微处理器被编程或配置以执行所述利用gps信号驯服晶振频率的方法的步骤，或该存储器上存储有被编程或配置以执行所述利用gps信号驯服晶振频率的方法的计算机程序。

此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行所述利用gps信号驯服晶振频率的方法的计算机程序。

和现有技术相比，本发明具有下述优点：本发明利用gps信号驯服晶振频率的方法主要包括四个部分：一是根据gps的1pps信号和晶振信号的误差特性建立两者的频差数学模型。通过对频差信号的数学建模为下一步建立频差信号的状态方程和量测方程提供基础；二是根据频差信号数学模型，对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差，得到受跳变野值和系统噪声影响较大的滤波模型；三是引入差分进化算法对常规扩展卡尔曼滤波中过程噪声进行动态调整，从而减少其对于频差信号滤波结果的影响；四是整合新息序列加权技术到上述扩展卡尔曼滤波算法中，通过在频差向量的最优估计方程中引入修正函数，来减少gps输出pps信号中包含的跳变野值对于滤波结果的影响。因此，本发明基于差分进化算法的新息序列加权扩展卡尔曼滤波算法，能够对gps信号中的跳变野值和随机误差、晶振信号中累积误差和非线性频率漂移进行滤波处理，能够提高gps锁定二级频标的最终校频精度。

附图说明

图1为现有技术利用gps锁定二级频标(高稳晶振)原理。

图2为本发明实施例方法的基本流程示意图。

图3为本发明实施例中引入差分进化算法的基本流程示意图。

具体实施方式

下面通过借助实施例更加详细地说明本发明，但以下实施例仅是说明性的，本发明的保护范围并不受这些实施例的限制。

如图2所示，本实施例利用gps信号驯服晶振频率的方法的步骤包括：

1)根据gps秒信号(1pps信号)、晶振信号的误差特性建立两者的频差数学模型；

2)根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差、得到受跳变野值和系统噪声影响较大的滤波模型；

3)引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整以减少其对于频差信号滤波结果的影响；

4)整合新息序列加权到完成过程噪声动态调整后的扩展卡尔曼滤波处理中，以减少跳变野值对频差信号滤波结果的影响。

本实施例利用gps信号驯服晶振频率的方法以滤除gps锁定晶振频率技术中频差信号的噪声为为目标，提出了一种基于差分进化算法的新息序列加权扩展卡尔曼滤波方法，可以对gps信号中的跳变野值和随机误差以及晶振信号中随机误差和非线性频率漂移进行处理，从而提高最终的校频精度。该方法主要包括四个部分：一是根据gps的1pps信号和晶振信号的误差特性建立两者的频差数学模型。通过对频差信号的数学建模为下一步建立频差信号的状态方程和量测方程提供基础；二是根据频差信号数学模型，对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差，得到受跳变野值和系统噪声影响较大的滤波模型；三是引入差分进化算法对常规扩展卡尔曼滤波中过程噪声进行动态调整，从而减少其对于频差信号滤波结果的影响；四是整合新息序列加权技术到上述扩展卡尔曼滤波算法中。通过在频差向量的最优估计方程中引入修正函数，来减少gps输出pps信号中包含的跳变野值对于滤波结果的影响。

本实施例中，设定频差信号及滤波相关参数时，设置过程噪声向量ω(k)，且其协方差矩阵为q＝0.0005。设置量测噪声向量γ(k)，设其协方差矩阵为r＝0.01；此外，频差初值设为x0＝5×10^-9s。

获取gps秒信号序列y时，gps秒信号序列y(大小为n)中的gps秒信号的世界统一时间utc可分别表示为：

1-ε1，2-ε2，…，k-εk，…，n-εn

获取晶振信号序列y′时，晶振信号序列y′(大小为n)中的晶振信号的世界统一时间utc可分别表示为：

1+a+b+c，1+a+2b+4c，…，k+a+bk+ck²，…，n+a+bn+cn²

其中，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；因此，本实施例中，步骤1)的详细步骤包括：

1.1)获取gps秒信号序列y，gps秒信号序列y中的任意k时刻的gps秒信号的世界统一时间utc的函数表达式如式(1)所示；

yk＝k-εk(1)

上式中，yk表示gps秒信号序列y中的k时刻的gps秒信号的世界统一时间utc，εk为k时刻的gps秒信号的误差，且满足ε～n(0,σ²)；

获取晶振信号序列y′，晶振信号序列y′中的任意k时刻的晶振信号的世界统一时间utc的函数表达式如式(2)所示；

y′k＝k+a+bk+ck²(2)

上式中，y′k表示晶振信号序列y′中的k时刻的晶振信号的世界统一时间utc，a为初始相位误差，b为考虑频率偏差的误差系数，c为考虑频率线性漂移的误差系数；

1.2)根据gps秒信号序列y、晶振信号序列y′做差得到频差序列，且频差序列中任意k时刻的频差信号的函数表达式如下式所示；

xk＝y′k-yk＝a+bk+ck²+εk(3)

上式中，xk表示频差序列中k时刻的频差信号，频差序列为x1，x2，…，xk，…，xn。

本实施例中，步骤2)用于根据频差信号数学模型，对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差，得到受跳变野值和系统噪声影响较大的滤波模型。步骤2)中根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理的步骤包括：

2.1)建立频差的状态方程和观测方程的函数表达式如式(4)所示；

上式中，x(k+1)表示k+1时刻的频差状态向量，x(k)表示k时刻的频差状态向量，f和h分别表示状态函数和量测函数，z(k+1)表示k+1时刻的量测值；ω(k)表示k时刻的过程噪声向量，其协方差矩阵为q；γ(k)表示k时刻的量测噪声向量，其协方差矩阵为r；

2.2)令k时刻到k+1时刻的状态转移矩阵f(k+1,k)、k时刻的观测矩阵h(k)的函数表达式如式(5)所示；

上式中，f表示状态函数，表示k时刻频差状态向量的估计，h(k)表示k时刻的观测矩阵，h表示量测函数，x(k)表示k时刻的频差状态向量；

2.3)建立常规扩展卡尔曼滤波模型的预测方程如式(6)～(8)所示、更新方程如式(9)～(10)所示，其中过程噪声如式(11)所示，过程噪声方差如式(12)所示；

式(6)中，表示基于k时刻预测值的k+1时刻频差状态向量预测值，f表示状态函数，表示k时刻频差状态向量的最优估计值；

p(k+1,k)＝f(k+1,k)p(k,k)f^t(k+1,k)+q(k)(7)

式(7)中，p(k+1,k)表示k+1时刻的协方差向量预测值，f(k+1,k)表示k时刻到k+1时刻的状态转移矩阵，p(k,k)表示k时刻的协方差向量预测值，f^t(k+1,k)表示k时刻到k+1时刻的状态转移矩阵的转置，q(k)表示k时刻的过程噪声方差；

k(k+1)＝p(k+1,k)h^t(k+1)[h(k+1)p(k+1,k)h^t(k+1)+r(k+1)]^-1(8)

式(8)中，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，p(k+1,k)表示k+1时刻的协方差向量预测值，h^t(k+1)表示k+1时刻的观测矩阵h(k+1)的转置，p(k+1,k)表示k+1时刻的协方差向量预测值，r(k+1)表示k+1时刻的观测噪声方差向量；

式(9)中，表示k+1时刻频差状态向量的最优估计值，表示k+1时刻频差状态向量的预测值，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，z(k+1)表示k+1时刻的量测值，表示由k+1时刻频差状态向量预测值得到的k+1时刻量测向量预测值；

p(k+1)＝[i-k(k+1)h(k+1)]p(k+1,k)(10)

式(10)中，p(k+1)表示k+1时刻更新的协方差向量估计值，i表示维数与k(k+1)h(k+1)相同的单位矩阵，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，h(k+1)表示k+1时刻的观测矩阵，p(k+1,k)表示k+1时刻的协方差向量预测值；

式(11)中，ω(k+1)表示k+1时刻的过程噪声，表示k+1时刻频差状态向量的最优估计值，z(k+1)表示k+1时刻的量测值；

qk+1＝var(ω(k+1))(12)

式(12)中，qk+1表示k+1时刻的过程噪声方差，ω(k+1)表示k+1时刻的过程噪声，var表示求过程噪声方差的操作。

本实施例中，步骤3)引入差分进化算法用于对常规扩展卡尔曼滤波算法中过程噪声进行动态调整，从而减少其对于频差信号滤波结果的影响。常规扩展卡尔曼在滤波过程中一般假设过程噪声均值为零，方差为q。由于系统模型存在不确定因素及异常扰动，因此过程噪声也是变化的；固定的过程噪声将会增大误差使滤波精度降低。为此引入差分进化算法对过程噪声的最优方差进行求解，然后将其应用在扩展卡尔曼滤波算法之中提高滤波精度，最终提高校频精度。如图3所示，步骤3)中引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整时，通过差分进化算法生成任意k时刻的过程噪声方差的步骤包括：

3.1)进行种群初始化，且函数表达式如下式(13)所示：

xi,k(0)＝lk+rand(0,1)*(uk-lk)(13)

式(13)中，xi,k(0)表示种群中第0代的第i条“染色体”的第k个基因；rand(0,1)为0～1之间均匀分布的随机数；uk和lk为搜索的上界和下界，且uk和lk分别取过程噪声方差的最大值qmax和最小值qmin；本实施例中设定差分进化算法初值时，变异率f＝0.5，交叉概率cr＝0.9，种群规模np＝10，迭代次数g＝100，搜索边界uk＝0.005,lk＝-0.005；

3.2)针对种群进行变异操作，且变异操作的函数表达式如下式(14)所示：

xi(g+1)＝xr1(g)+f*[xr2(g)-xr3(g)](14)

式(14)中，xi(g+1)为得到的g+1代变异个体；f为压缩比例因子，取值范围0～1；xr1(g)、xr2(g)和xr3(g)为3个父代；

3.3)针对种群进行交叉操作以保留较优良的变量，且交叉操作采用二项交叉方式，其函数表达式如下式(15)所示：

式(15)中，yi,j(g+1)表示当前个体对应的交叉操作结果，xi,j(g+1)表示变异个体对应的分量，xi,j(g)表示当前个体对应的分量，r为每个变量生成的一个0～1之间均匀分布的随机数；cr为变量的交叉概率；rnd为1～d之间均匀分布的整数，如果r＜cr则接受变异个体对应的分量xi,j(g+1)，否则保留当前个体对应的分量xi,j(g)；

3.4)针对交叉操作得到的结果，采用如下式(16)所示贪婪选择的方式选择最优解并作为k时刻的过程噪声方差；

式(16)中，xi(g+1)为选择结果，yi(g+1)表示交叉操作结果对应的个体，xi(g)表示当前个体，f[yi(g+1)]表示个体yi(g+1)对应的种群适应度值，f[xi(g)]表示个体xi(g)对应的种群适应度值。将最优解赋值给qi,代入协方差预测方程，即可提高扩展卡尔曼的滤波精度。

本实施例中，步骤4)整合新息序列加权技术到上述扩展卡尔曼滤波算法中，以减少gps输出pps信号中包含的跳变野值对于滤波结果的影响。在实际应用中，由于gps信号受噪声干扰还会产生较大的跳变野值，带有野值的频差信号将会使滤波结果发生偏移甚至发散。为了解决上述扩展卡尔曼算法的抗野值问题，有学者将新息序列加权的方法引入卡尔曼滤波来消除野值的影响。当观测数据不包含野值时，滤波算法能够正常运行，对观测噪声进行滤除；当观测数据包含野值时，能够克服或减少野值对于滤波结果的影响，使滤波结果更加准确。本实施例中，步骤4)整合新息序列加权到完成过程噪声进行动态调整后的扩展卡尔曼滤波处理中后，得到的扩展卡尔曼滤波处理的更新方程的函数表达式如下式(17)～(21)所示；

dk+1＝p(k+1,k)+rk+1(21)

式(17)～(21)中，表示调整后的k+1时刻频差状态向量的最优估计值，表示k+1时刻频差状态向量预测值，k(k+1)表示k+1时刻卡尔曼滤波增益矩阵，φk+1表示修正函数，τk+1表示衡量k+1时刻量测误差大小的变量，ek+1表示k+1时刻的量测误差，z(k+1)表示k+1时刻频差向量的量测值，表示k+1时刻量测向量的预测值，τk+1服从χ²分布，当出现野值时τk+1增大、φk+1(·)减小；ck+1为选取的门限值或常数序列，λk+1为矩阵的最大特征值，kk+1表示k+1时刻的卡尔曼滤波增益，dk+1表示k+1时刻衡量量测误差大小的权矩阵，p(k+1,k)表示k+1时刻的协方差预测值，r(k+1)表示k+1时刻的量测噪声方差。

综上所述，传统应用于低成本高精度频率标准的主要技术为gps锁定二级频标(晶振)技术，为减少gps锁定二级频标技术中gps信号的随机误差和跳变野值以及晶振信号的非线性频率偏移对于频差信号的影响。本实施例提供了一种对gps锁定晶振频率技术中的频差信号进行处理的基于差分进化算法的新息序列加权扩展卡尔曼滤波方法，本实施例方法主要包括四个部分：一是根据gps的1pps信号和晶振信号的误差特性建立两者的频差数学模型。通过对频差信号的数学建模为下一步建立频差信号的状态方程和量测方程提供基础；二是根据频差信号数学模型，对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差，得到受跳变野值和系统噪声影响较大的滤波模型；三是引入差分进化算法对常规扩展卡尔曼滤波中过程噪声进行动态调整，从而减少其对于频差信号滤波结果的影响；四是整合新息序列加权技术到上述扩展卡尔曼滤波算法中。通过在频差向量的最优估计方程中引入修正函数，来减少gps输出pps信号中包含的跳变野值对于滤波结果的影响。而且，本实施例方法采用扩展卡尔曼滤波算法对gps锁定晶振频率过程中的频差信号进行处理，在对gps信号的随机误差进行处理的同时也考虑了晶振由于老化和温度影响而产生的非线性误差。本实施例方法将新息序列加权和差分进化算法引入扩展卡尔曼滤波过程，可以减少gps信号中较大跳变野值以及系统噪声不断变化对于滤波结果的影响。

此外，本实施例还提供一种利用gps信号驯服晶振频率的系统，包括：

频差模型生成程序单元，用于根据gps秒信号、晶振信号的误差特性建立两者的频差数学模型；

常规滤波程序单元，用于根据频差数学模型对频差信号进行常规的扩展卡尔曼滤波处理，从而消除随机误差和晶振频率非线性漂移引起的误差、得到受跳变野值和系统噪声影响较大的滤波模型；

噪声调整程序单元，用于引入差分进化算法对常规的扩展卡尔曼滤波处理中的过程噪声进行动态调整以减少其对于频差信号滤波结果的影响；

滤波整合程序单元，用于整合新息序列加权到完成过程噪声进行动态调整后的扩展卡尔曼滤波处理中。

此外，本实施例还提供一种利用gps信号驯服晶振频率的系统，包括计算机设备，该计算机设备被编程或配置以执行前述利用gps信号驯服晶振频率的方法的步骤。

此外，本实施例还提供一种利用gps信号驯服晶振频率的系统，包括计算机设备，该计算机设备的存储器上存储有被编程或配置以执行前述利用gps信号驯服晶振频率的方法的计算机程序。

此外，本实施例还提供一种带有gps的智能终端，包括gps模块、微处理器和存储器，该微处理器被编程或配置以执行前述利用gps信号驯服晶振频率的方法的步骤，或该存储器上存储有被编程或配置以执行前述利用gps信号驯服晶振频率的方法的计算机程序。该智能终端可以为智能手机、平板电脑或者其他类型的终端设备。

此外，本实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行前述利用gps信号驯服晶振频率的方法的计算机程序。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。