CN104899467A - 一种基于频率驾驭的频率偏移校正方法及系统 - Google Patents

Abstract

本发明公开一种基于频率驾驭的频率偏移校正方法及系统，该方法包括如下步骤：S1、比对得出被控频率信号与参考信号之间的历史时差序列，并基于支持向量机原理建立历史时差序列模型；S2、利用历史时差序列模型，预测下一阶段被控频率信号与参考信号之间的时差序列；S3、基于下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量；S4、基于被控频率信号的调整量校正被控频率信号。本发明所述技术方案实施简单、推广性能好、学习速度快、优化求解时具有唯一的极小点，而且所构造的模型有很好的预测性能。

Description

一种基于频率驾驭的频率偏移校正方法及系统

技术领域

本发明涉及频率源控制领域。更具体地，涉及一种基于频率驾驭的频率偏移校正方法及系统。

背景技术

维持一个稳定可靠的时间基准系统，需要一个精准的频率源。而自由运行的频率源在其自身速率的影响下，其产生的频率信号或多或少存在着频率偏移现象。频率驾驭是指将一个有着更高准确度和稳定度的信号作为参考信号，通过一定的算法获得被控频率信号相对于参考信号的偏移补偿量，以此干预被控频率信号的自由走势，在保证其短期稳定度的前提下，提高其准确度和长期稳定度。

频率调整量的计算，其关键在于利用当前阶段获取的历史时差数据合理建模，来预测下一阶段的时差。传统的频率驾驭方法中，通常采用最小二乘拟合模型和卡尔曼滤波模型对频率信号的走势进行预测。最小二乘拟合模型需要采集大量样本，事先确定样本的多项式次数，再利用历史数据拟合确定各次项的系数，进而作外推预报估计下一时刻的时差。卡尔曼滤波模型则需要事先确定模型噪声协方差矩阵和观测噪声协方差矩阵，这两个矩阵的获取并没有固定的方法，只能采用经验估计。这两种预测方法除了需要很强的先验知识外，还需要预先知道模型的结构形式。但在实际预测问题中，样本的数量都是有限的，根据有限样本事先确定的模型结构并不能很好的总结归纳实际情况；因此以上两种预测模型存在共同的弊端，固定结构的预测模型缺乏对实际情况的自适应能力，这将导致预报误差随时间的增加而不断增大。

因此，需要提供一种基于频率驾驭的频率偏移校正方法及系统。

发明内容

本发明的目的在于提供一种基于频率驾驭的频率偏移校正方法及系统，采用支持向量机原理建立预测模型，解决传统的频率偏移校正中需事先确定预测模型结构，自适应性差的问题。

为达到上述目的，本发明采用下述技术方案：

一种基于频率驾驭的频率偏移校正方法，该方法包括如下步骤：

S1、比对得出被控频率信号与参考信号之间的历史时差序列，并基于支持向量机原理建立历史时差序列模型，模型公式为：

$y_i=f\left(x\right)=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i-{\mathrm{\α}}_i^{*})K(x,x_i)+b;$

公式中，l为历史时差序列的个数，y_i为第i个时差值，x_i为第i个时差值对应的时标，α_i为拉格朗日乘子；K(x,x_i)为高斯径向基核函数；b为常数项；

S2、利用历史时差序列模型，预测下一阶段被控频率信号与参考信号之间的时差序列；

S3、基于下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量；

S4、基于被控频率信号的调整量校正被控频率信号。

优选地，步骤S1中的参考信号包括内部参考信号和外部参考信号，内部参考信号为综合时间尺度源生成的综合时间尺度；外部参考信号为国际时钟源生成的国际时钟信号。

优选地，步骤S2中下一阶段被控频率信号与参考信号之间的时差序列为：

{(x_l+1,y_l+1),…,(x_l+i,y_l+i)}∈(X,Y)^N-l；

其中，第l+i个时差值对应的时标x_l+i＝x_l+i×τ，i＝1,…,N-l，N为历史时差序列与下一阶段时差序列的总个数，τ为固定的采样间隔；第l+i个时差值y_l+i，i＝1,…,N-l通过将第l+i个时差值对应的时标x_l+i代入历史时差序列模型计算得到。

优选地，步骤S3中计算被控频率信号的调整量的方法为：

对下一阶段被控频率信号与参考信号之间的时差序列做线性拟合，拟合曲线的一次项系数a₀即为被控频率信号相对于参考信号的预测速率：

f`(x)＝a₀*x+a₁

公式中，a₀为拟合曲线的一次项系数，a₁为拟合曲线的常数项；

预测速率a₀的相反数-a₀即为被控频率的频率调整量。

一种实施上述方法的基于频率驾驭的频率偏移校正系统，该系统包括：

被控频率源、综合时间尺度源、国际时钟源、内部时间比对模块、外部时间比对模块、频率调整量计算模块、频率调整控制模块；

被控频率源生成被控频率信号；

综合时间尺度源生成综合时间尺度，作为内部参考信号；国际时钟源生成国际时钟信号，作为外部参考信号；

内部时间比对模块获取被控频率信号与内部参考信号之间比对的历史时差；

外部时间比对模块获取被控频率信号与外部参考信号之间比对的历史时差；

频率调整量计算模块根据内部时间比对模块和外部时间比对模块获取的历史时差数据，预测下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量；

频率调整控制模块根据调整量对被控频率信号进行调整控制。

本发明的有益效果如下：

本发明所述技术方案中建立预测模型采用的支持向量机是建立在结构风险最小化原理基础上的结构化学习方法，它能利用核函数很好地解决有限数量样本的高维建模问题，实施简单、推广性能好、学习速度快、优化求解时具有唯一的极小点，而且所构造的模型有很好的预测性能。此外，本发明所述技术方案中的参考信号是相较于被控频率信号有着更高准确度和更好稳定性的信号，该参考信号包括多台原子钟组合计算得出的综合时间尺度TA和协调世界时UTC。

附图说明

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

图1示出基于频率驾驭的频率偏移校正方法流程图。

图2示出基于频率驾驭的频率偏移校正系统示意图。

具体实施方式

为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。

本实施例提供的基于频率驾驭的频率偏移校正方法，包括如下步骤：

Step1、比对得出被控频率信号与参考信号之间的历史时差序列，并基于支持向量机原理建立历史时差序列模型：

被控频率信号由被控频率源生成，

参考信号包括内部参考信号和外部参考信号两种，内部参考信号为综合时间尺度源生成的综合时间尺度，综合时间尺度的生成过程为利用多台原子钟生成的信号之间的比对数据进行基于AT1算法的计算得到；外部参考信号为国际时钟源生成的国际时钟信号，

进而历史时差序列包括：被控频率信号与综合时间尺度之间的历史时差序列、被控频率信号与国际时钟信号之间的历史时差序列。

为保证模型对实际情况的适应性，至少需要积累一个月的数据。这些历史数据称为训练集T＝{(x₁,y₁),…,(x_i,y_i)}∈(X,Y)^l，其中，l为历史时差序列的个数，y_i∈R表示第i个时差值，单位为ns，x_i∈Rⁿ表示第i个时差值对应的时标，用约化儒略日表示。历史时差序列模型为利用实值函数f(x)来表达时差值y对时差值对应的时标x的依赖关系。

在一维空间中，训练集T中的(x_i,y_i)是一种非线性关系，需要使用一个非线性映射φ(·)将训练集的x_i由输入空间X映射到高维特征空间H，在特征空间H中构造线性支持向量回归模型，该模型对应一维输入空间中的非线性回归。这个映射由满足Mercer条件下构成的核函数K(x_i,x_j)＝(φ(x_i)·φ(x_j))来实现。考虑到可能存在误差，引入两个松弛变量：求解最优实值函数f(x)＝<ω·φ(x)>+b，即最小化

$\frac{1}{2}\mathrm{\&omega;}\&CenterDot;\mathrm{\&omega;}+C\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})---\left(1\right)$

约束条件为

$\left\{\begin{array}{c}f\left(x_i\right)-y_i\&le;{\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}\\ f\left(x_i\right)-y_i\&le;{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;},\\ {\mathrm{\&xi;}}_i^{*},{\mathrm{\&xi;}}_i\&GreaterEqual;0\end{array}\right.,i=1,...,l---\left(2\right)$

公式(1)中，第一项使函数更为平坦，提高泛化能力；第二项为减小误差，常数C为惩罚参数，对两者做折中。公式(2)中ε为一正常数，f(x_i)与y_i的差别小于ε时不计入误差，大于ε时误差计为|f(x_i)-y_i|-ε。采用拉格朗日乘子解决这个二次规划问题，即

$\begin{array}{c}L(\mathrm{\&omega;},b,\mathrm{\&alpha;},{\mathrm{\&alpha;}}^{*})=\frac{1}{2}{||\mathrm{\&omega;}||}^2+C\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})-\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i({\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+\&lang;\mathrm{\&omega;}\&CenterDot;\mathrm{\&phi;}\left(x_i\right)\&rang;+b)\\ -\underset{i=1}{\overset{l}{\&Sigma;}}{\mathrm{\&alpha;}}_i^{*}({\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}+y_i-\&lang;\mathrm{\&omega;}\&CenterDot;\mathrm{\&phi;}\left(x_i\right)\&rang;-b)-\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i{\mathrm{\&gamma;}}_i+{\mathrm{\&xi;}}_i^{*}{\mathrm{\&gamma;}}_i^{*})\end{array}---\left(3\right)$

函数L的极值应满足条件：

$\left\{\begin{array}{c}\frac{\&part;L}{\&part;\mathrm{\&omega;}}=0\&DoubleRightArrow;\mathrm{\&omega;}=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)\mathrm{\&phi;}\left(x_i\right)\\ \frac{\&part;L}{\&part;b}=0\&DoubleRightArrow;\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)=0\\ \frac{\&part;L}{\&part;{\mathrm{\&xi;}}_i}=0\&DoubleRightArrow;C-{\mathrm{\&alpha;}}_i-{\mathrm{\&gamma;}}_i=0\\ \frac{\&part;L}{\&part;{\mathrm{\&xi;}}_j^{*}}=0\&DoubleRightArrow;C-{\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&gamma;}}_i^{*}=0\end{array}\right.---\left(4\right)$

于是得到优化问题的对偶形式，最大化函数：

$W(\mathrm{\&alpha;},{\mathrm{\&alpha;}}^{*})=-\frac{1}{2}\underset{i,j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})({\mathrm{\&alpha;}}_j-{\mathrm{\&alpha;}}_j^{*})K(x_i,x_j)+\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})y_i-\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})\mathrm{\&epsiv;}---\left(5\right)$

约束条件为

$\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})=0,0\&le;{\mathrm{\&alpha;}}_i,{\mathrm{\&alpha;}}_i^{*}\&le;C,i=1,...,l---\left(6\right)$

从而构造非线性回归函数，即历史时差数据模型，公式如下：

$y_i=f\left(x\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x,x_i)+b---\left(7\right)$

公式中，α_i为拉格朗日乘子；K(x,x_i)为高斯径向基核函数：σ`为函数的宽度参数；b为模型的常数项；

Step2、利用历史时差序列模型，预测下一阶段被控频率信号与参考信号之间的时差序列：

内部时间比对模块和外部时间比对模块存储的时差序列T＝{(x₁,y₁),…,(x_i,y_i)}∈(X,Y)^l，时标x_i有唯一固定的采样间隔τ。步骤Step1得到了该序列的支持向量机回归模型，即时标x_i与时差值y_i之间的依赖关系。根据训练集T的最后一个采样时刻x_l和采样间隔τ，可以得到下一阶段采样点{(x_l+1,y_l+1),…,(x_l+i,y_l+i)}∈(X,Y)^N-l的时标：

x_l+i＝x_l+i×τ，i＝1,…,N-l     (8)

公式中，N为历史时差序列与下一阶段时差序列的总个数；

将公式(8)代入公式(7)，得到下一阶段各采样时刻对应的时差预测值：

y_l+i，i＝1,…,N-l；

从而预测出下一阶段被控频率信号与参考信号之间的时差序列{(x_l+1,y_l+1),…,(x_l+i,y_l+i)}∈(X,Y)^N-l；

Step3、基于下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量：

利用步骤Step2得到下一阶段被控频率信号与参考信号之间的时差预测序列T＝{(x_l+1,y_l+1),…,(x_l+i,y_l+i)}∈(X,Y)^N-l。对预测序列做线性拟合，拟合曲线(公式(9))的一次项系数a₀即被控频率信号相对于参考信号的预测速率：

f`(x)＝a₀*x+a₁     (9)

公式中，a₀为拟合曲线的一次项系数，a₁为拟合曲线的常数项；

则预测速率a₀的相反数-a₀即被控频率的频率调整量。

Step4、基于被控频率信号的调整量校正被控频率信号。

本实施例提供的实施上述方法的基于频率驾驭的频率偏移校正系统，包括：

自由运行的被控频率源、综合时间尺度源、国际时钟源、内部时间比对模块、外部时间比对模块、频率调整量计算模块、频率调整控制模块；

自由运行的被控频率源生成被控频率信号；

综合时间尺度源生成的综合时间尺度由多台原子钟组合计算得到，综合时间尺度作为内部参考信号；

国际时钟源生成国际时钟信号，作为外部参考信号；

内部时间比对模块获取被控频率信号与内部参考信号之间比对的历史时差，外部时间比对模块获取被控频率信号与外部参考信号之间比对的历史时差；

频率调整控制模块用于对被控频率信号进行调整控制；

内部时间比对模块和外部时间比对模块存储历史时差序列，以备计算频率调整量时所用；

频率调整量计算模块提取内部时间比对模块和外部时间比对模块的历史时差数据，经过数据预处理后，建立数学模型，解算频率调整量，并实时将该调整量送入频率调整控制模块，校正被控频率信号的频率偏移；

自由运行的频率源输出端与频率调整控制设备的输入端电缆线连接，频率调整控制设备的输出端分别与内部时间比对模块和外部时间比对模块导线连接，内部时间比对模块数据和外部时间比对模块数据分别与频率调整量计算模块数据导线连接，频率调整量计算模块数据与频率调整控制设备输入端串口线连接，频率驾驭控制系统形成反馈回路。

下面将本实施例提供的基于频率驾驭的频率偏移校正方法代入具体的计算环境进行进一步说明：

由于测量环境影响，内部时间比对模块和外部时间比对模块获取的时差序列可能会存在异常值。频率调整量计算方法第一步所述的数据预处理，即对时差序列中异常值的剔除和补足。剔除异常值采用3σ法则，对剔除的异常值补足采用平均值插值法。对于训练集的时差序列T＝{(x₁,y₁),…,(x_i,y_i)}∈(X,Y)^l，对异常值的处理步骤如下：

首先，计算训练集T中所有时差{y₁,…,y_i}∈Y^l的均值

$\stackrel{\&OverBar;}{y}=\frac{1}{l}\underset{i=1}{\overset{l}{\&Sigma;}}{y}_i$

其次，计算训练集T中所有时差{y₁,…,y_i}∈Y^l的方差σ²；

${\mathrm{\&sigma;}}^2=\frac{1}{l}\underset{i=1}{\overset{1}{\&Sigma;}}{(y_i-\stackrel{\&OverBar;}{y})}^2$

再次，判断训练集T中的时差{y₁,…,y_i}∈Y^l是否在3σ控制线以内。若则y_i为异常值，将其剔除。

最后，用异常值前后两个正常值的均值补足剔除的异常值。

频率调整量计算方法第二步和第三步所述的利用支持向量机原理对历史时差序列建模并预测是其中的关键技术。建模及预测过程利用支持向量机加强版工具箱libsvm-mat-2.89-3在matlab环境下编程实现，步骤如下：

首先，确定模型的输入量、输出量。输入量为经过预处理后的训练集T；输出量为下一阶段的时差预报值，称其为测试集。

其次，将训练集中的时差序列转换为libsvm-mat-2.89-3工具箱所需要的格式。libsvm-mat-2.89-3工具箱要求的数据格式为：<label><index1>:<value1><index2>:<value2>…。其中，<label>是目标值，在训练集中，<label>有确切实数值，在测试集中，由于目标值未知，<label>可设为任意实数，也可空着不填；<index>是以1开始的整数，可以是不连续的；<value>为实数。对应到训练集和测试集的数据对(x,y)，<label>即(x,y)中的因变量y，<value>即(x,y)中的自变量x，由于训练集T中自变量只有x一个，所以<index>为1。例如，训练集T中某个时差数据对(57069.25,42.3)转换为libsvm-mat-2.89-3工具箱要求的格式为：(42.31:57069.25)；测试集中的未知时差数据对(57069.5,)可转换为(1:57069.5)或(21:57069.5)，目标值为空或任意实数，如“2”。

再次，选择高斯径向核函数，利用libsvm-mat-2.89-3工具箱建立支持向量机模型，利用网格法进行参数寻优，并根据训练集采用交叉验证法(CrossValidation)对模型进行训练。该过程结束后，得到训练集T＝{(x₁,y₁),…,(x_i,y_i)}∈(X,Y)^l的模型model。

最后，利用预测函数svmpredict，根据模型model对测试集中的目标值<label>进行预测。

对于频率调整量计算方法第四步所述的线性拟合，即对测试集中的预测序列T'＝{(x_l+1,y′_l+1),…,(x_l+i,y′_l+i)}∈(X,Y')^N-l，找到一个线性函数y'＝S^*(x)，使得误差平方和最小，

其中

S(x)＝a₀x+a₁

实现过程如下：

首先，将考虑成加权平方和；

然后，将求解最小的问题转换为多元函数求极小值；

$\left\{\begin{array}{c}I(a_0,a_1)=\underset{i=l+1}{\overset{N}{\&Sigma;}}\mathrm{\&omega;}\left(x_i\right){\&lsqb;\underset{j=0}{\overset{1}{\&Sigma;}}{a}_j{x}_i-y_i^{\&prime;}\&rsqb;}^2\\ \frac{\&part;I}{\&part;a_k}=2\underset{i=l+1}{\overset{N}{\&Sigma;}}\mathrm{\&omega;}\left(x_i\right)\&lsqb;\underset{j=0}{\overset{1}{\&Sigma;}}{a}_j{x}_i-y_i^{\&prime;}\&rsqb;x_i=0,(k=0,1)\end{array}\right.$

从上式中解出唯一解，即得测试集的线性拟合曲线。其一次项系数a₀就是时差序列的预测速率，被控频率的调整量为-a₀。

如频率调整量计算方法第一步所述将调整量-a₀送入频率调整控制设备，修正频率信号的偏移，数据的传送过程通过RS232串口实现。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (5)

1.一种基于频率驾驭的频率偏移校正方法，其特征在于，该方法包括如下步骤：

S1、比对得出被控频率信号与参考信号之间的历史时差序列，并基于支持向量机原理建立历史时差序列模型，模型公式为：

$y_i=f\left(x\right)=\underset{i=1}{\overset{1}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x,x_i)+b;$

公式中，l为历史时差序列的个数，y_i为第i个时差值，x_i为第i个时差值对应的时标，α_i为拉格朗日乘子；K(x,x_i)为高斯径向基核函数；b为常数项；

S2、利用历史时差序列模型，预测下一阶段被控频率信号与参考信号之间的时差序列；

S3、基于下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量；

S4、基于被控频率信号的调整量校正被控频率信号。

2.根据权利要求1所述的基于频率驾驭的频率偏移校正方法，其特征在于，步骤S1中的参考信号包括内部参考信号和外部参考信号，内部参考信号为综合时间尺度源生成的综合时间尺度；外部参考信号为国际时钟源生成的国际时钟信号。

3.根据权利要求1所述的基于频率驾驭的频率偏移校正方法，其特征在于，步骤S2中下一阶段被控频率信号与参考信号之间的时差序列为：

{(x_l+1,y_l+1),…,(x_l+i,y_l+i)}∈(X,Y)^N-l；

其中，第l+i个时差值对应的时标x_l+i＝x_l+i×τ，i＝1,…,N-l，N为历史时差序列与下一阶段时差序列的总个数，τ为固定的采样间隔；第l+i个时差值y_l+i，i＝1,…,N-l通过将第l+i个时差值对应的时标x_l+i代入历史时差序列模型计算得到。

4.根据权利要求1所述的基于频率驾驭的频率偏移校正方法，其特征在于，步骤S3中计算被控频率信号的调整量的方法为：

对下一阶段被控频率信号与参考信号之间的时差序列做线性拟合，拟合曲线的一次项系数a₀即为被控频率信号相对于参考信号的预测速率：

f`(x)＝a₀*x+a₁

公式中，a₀为拟合曲线的一次项系数，a₁为拟合曲线的常数项；

预测速率a₀的相反数-a₀即为被控频率的频率调整量。

5.一种实施权利要求1所述方法的基于频率驾驭的频率偏移校正系统，其特征在于，该系统包括：

被控频率源、综合时间尺度源、国际时钟源、内部时间比对模块、外部时间比对模块、频率调整量计算模块、频率调整控制模块；

被控频率源生成被控频率信号；

综合时间尺度源生成综合时间尺度，作为内部参考信号；国际时钟源生成国际时钟信号，作为外部参考信号；

内部时间比对模块获取被控频率信号与内部参考信号之间比对的历史时差；

外部时间比对模块获取被控频率信号与外部参考信号之间比对的历史时差；

频率调整量计算模块根据内部时间比对模块和外部时间比对模块获取的历史时差数据，预测下一阶段被控频率信号与参考信号之间的时差序列，计算被控频率信号的调整量；

频率调整控制模块根据调整量对被控频率信号进行调整控制。