基于Vondark-Cepek滤波的氢铯时间尺度融合方法与流程

本发明属于时间频率、信号处理领域，涉及一种时间尺度融合方法。

背景技术：

目前，高精度时间已成为一个国家科技、经济、军事和社会生活中至关重要的参量，广泛应用于导航、电力、通信、航空、国防等领域，作为最基本的物理量，对提升国家科研水平也有基础性作用。为此，世界各时间频率实验室都拥有自己的原子钟组，每台原子钟都可以产生一个时间尺度。但守时不能仅依赖单台原子钟，因为每个物理装置都有可能出现异常。时间尺度算法则是将守时钟组各台钟的时间根据需求，计算出一个标准时间，目的在于寻找一种方式使得算法的不确定度最小，稳定度最高。

氢铯原子钟是两种不同类型的原子频标，二者分别有着很好的短期稳定度和长期稳定度，它们是当前全球各时间实验室最常见的两种守时型原子钟。氢铯时间尺度的融合方法主要是研究两类不同原子钟有效组合而产生更稳定更准确的时间尺度。一个时间实验室所强调的是时间尺度的稳定度。

目前很多守时实验室采用氢铯互为参考的方式进行两种不同类型原子钟时间尺度的融合：分别以氢钟、铯钟或者铯钟组时间作为参考，对另一类原子钟进行钟差测量形成时间序列后，对守时钟组中的铯钟短期波动进行修正或氢钟的速率和频率漂移进行扣除，然后通过加权平均来获得最终的时间尺度。在上述的研究方法中，主要是以氢或铯为基准，对另一种原子钟的单台钟进行评估和修正，然后将守时钟组中的所有钟一起加权平均来归算最终的时间尺度。该方法的不足在于：当氢原子钟作为参考对铯原子钟进行测量的主要噪声是相位白噪声，经数学方法滤波后，其时间尺度短期稳定度仍会受铯原子钟噪声的影响。而以铯钟组时间尺度来修正氢钟，当氢钟出现故障时，产生时间尺度的可靠性不能保证。

技术实现要素：

为了克服现有技术的不足，本发明提供一种基于vondark-cepek滤波的氢铯时间尺度融合方法，有效抑制噪声的同时能够兼顾氢铯两种类型原子钟特性，在钟组尺度层面进行融合，产生长短期稳定度都良好的时间尺度。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1，将铯原子钟和氢原子钟分组，分别经钟差预测、频率估计及权重估计步骤产生指数滤波的铯钟组时间尺度tacs和氢钟组时间尺度tah；

步骤2：将氢钟组时间尺度序列进行一次差分，获得对应频率数据；

步骤3，将铯钟组时间尺度序列经傅里叶变换转化到频域，经频谱分析确定需抑制信号f或周期p；根据最小二乘理论计算铯钟组时间尺度平滑因子和氢钟组时间尺度平滑因子其中t和t’分别是tacs和tah对应的频率响应；

步骤4，将铯钟组时间尺度和氢钟组频率序列作为输入，将平滑因子用于vondark-cepek组合滤波方法，产生时间尺度融合结果，实现铯钟组和氢钟组时间尺度的融合。

所述的步骤1之前对实验室测得的原子钟钟差数据进行预处理，分别选择铯原子钟组和氢原子钟组中长短期稳定度最佳的原子钟作为参考钟，利用时间间隔计数器测量同一钟组内其他原子钟与该参考钟的钟差数据。

所述的预处理得到的钟差数据采用3σ法则去除异常点，同时检测钟差数据的连续性，对于连续缺失数据小于5个测量周期的钟差序列采用插值法进行修补；对于钟差数据连续缺失大于5个测量周期的，则将该台钟在守时钟组中予以剔除。

所述的步骤1中，若每台原子钟的权重大于设定最大权重上限，则该台原子钟的权重设置为最大权；设h’i(t)为t时刻原子钟hi时间改正量，则将权重、原子钟数据、h’i(t)通过公式计算，其中，xij(t)为原子钟hj和原子钟hi在t时刻的钟差值，n为参与计算原子钟数量，xi(t)为实验室单个原子钟与实验室原子钟组平均时间尺度之差，wi为每台原子钟的权重。

所述的步骤1中，最大权重设为2.5/n。

所述的步骤4中，定义铯钟组时间尺度tacs数据为m(ti)，是时间ti的函数，i＝1,2,3,…，其在ti时刻的导数m′(t)＝[m(ti+1)-m(ti)]/(ti+1-ti)；定义m(tj)是氢钟组时间尺度(tah)的数据，其导数m′(t)＝[m(tj+1)-m(tj)]/(tj+1-tj)，其中，j＝1,2,3,…；当ti＝tj时，m′(t)＝m′(t)，m(mjd2)＝m(mjd1)+m′(tj)δtj，mjd表示简化儒略日；假设输入观测序列1表示为y′j，对应时标为xj，权重为pj；观测序列2一阶导数表示为对应时标为权重为输出序列则是要得到的平滑曲线yi，对应时标为xi，定义曲线的平滑度式中，的表达式是未知的，需根据平滑的数据对其三阶导数进行估计；在两点[xi+1,yi+1]和[xi+2,yi+2]之间的平滑曲线被定义为由邻近的四个点i、i+1、i+2、i+3得到的一个三阶的拉格朗日多项式其三阶导数定义平滑曲线对观测值的逼真度定义平滑曲线对观测值一阶导数的逼真度

本发明的有益效果是：实现了氢铯时间尺度层面的融合，讨论了铯钟组时间尺度的绝对平滑和绝对逼真以及氢钟组时间尺度一阶导数的绝对拟合这三种情况下的最佳取权。在指数滤波算法的基础上，进一步实现对钟组时间尺度的滤波，有效地改善了铯钟短期波动对融合结果的影响。同时，将两种时间尺度的融合，进一步提升了仅一次加权平均时间尺度的可靠性。该方法既能够有效地抑制噪声，又能充分利用氢铯原子钟的特性，使融合时间尺度长短期稳定度都好。

附图说明

图1为本发明产生融合时间尺度的流程示意图(以三台铯钟和三台氢钟为例)。

图2为时间尺度示意图，其中(a)为铯钟组时间尺度，(b)为氢钟组时间尺度。

图3为本发明时间尺度融合结果示意图，，其中(a)为融合结果，(b)为allan偏差。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明提供一种基于vondark-cepek滤波的氢铯时间尺度融合方法，首先利用指数滤波法分别对氢钟组和铯钟组各产生一个钟组时间尺度，然后根据最小二乘原则对vondark-cepek组合滤波关键参数进行选取，进而通过氢钟组时间尺度时间序列的差分信息对铯钟组时间尺度进行性能增强，从而获得氢铯融合时间尺度。

本发明具体包括以下步骤：

步骤1：对实验室测得的原子钟钟差数据进行预处理。首先选择钟组中长短期稳定度最佳的原子钟作为参考钟，利用时间间隔计数器测量其他原子钟与该参考钟的钟差数据。钟差测量每小时一次。采用3σ法则去除异常点，同时检测钟差数据的连续性。对于连续缺失数据小于5小时的钟差序列采用插值法进行修补；对于钟差数据连续缺失5小时以上的情况，则将该台钟在守时钟组中予以剔除。

步骤2：将铯原子钟和氢原子钟分组，分别经钟差预测、频率估计及权重估计步骤(每小时计算一次)产生指数滤波的铯钟组时间尺度tacs和氢钟组时间尺度tah。(以下以计算3个月的时间尺度为例)

若每台原子钟的权重大于设定最大权重上限，则该台原子钟的权重设置为最大权。根据本方法的优选，我们最大权重设为2.5/n。其中n为参与计算原子钟数量。设h’i(t)为t时刻原子钟hi时间改正量，则将权重、原子钟数据、h’i(t)通过如下公式计算：

其中，xij(t)为原子钟hj和原子钟hi在t时刻的钟差值；xi(t)为实验室单个原子钟与实验室原子钟组平均时间尺度之差，wi为每台原子钟的权重。

步骤3：将步骤2计算得到的氢钟组时间尺度序列进行一次差分，获得对应频率数据。

步骤4：将铯钟组时间尺度序列经傅里叶变换转化到频域，经频谱分析确定需抑制信号f或周期p。根据该频率f或周期p、经验及实际算法确定tacs和tah对应的频率响应t和t’，根据最小二乘理论计算铯钟组时间尺度平滑因子ε和氢钟组时间尺度平滑因子

步骤5：将步骤2计算产生的铯钟组时间尺度和步骤3产生的氢钟组频率序列作为输入，将步骤4获得的平滑因子用于vondark-cepek组合滤波方法，产生时间尺度融合结果，实现铯钟组和氢钟组时间尺度的融合。

建立氢铯时间尺度融合模型:定义铯钟组时间尺度(tacs)数据为m(ti)，是时间ti的函数，i＝1,2,3,…，其在ti时刻的导数定义为m′(t)＝[m(ti+1)-m(ti)]/(ti+1-ti)。m(tj)是氢钟组时间尺度(tah)的数据，它的导数定义为m′(t)＝[m(tj+1)-m(tj)]/(tj+1-tj)。其中，j＝1,2,3,…。m′(t)和m′(t)的物理意义是一致的，都是时间尺度序列的速率。因此，当ti＝tj时，m′(t)＝m′(t)。

m(mjd2)＝m(mjd1)+m′(tj)δtj(3)

其中mjd为简化儒略日。

因所选取铯钟组和氢钟组时间尺度的时间间隔相同，故有公式(3)铯钟组下一时刻的钟差估计可用本时刻铯钟组钟差与氢钟组时间尺度速率和时间间隔乘积的和表示。

假设输入观测序列1表示为y′j，对应时标为xj，权重为pj；观测序列2一阶导数表示为对应时标为权重为输出序列则是要得到的平滑曲线表示为yi，对应时标为xi。定义三个量：

1)曲线的平滑度：

式中，的表达式是未知的，需根据平滑的数据对其三阶导数进行估计。在两点[xi+1,yi+1]和[xi+2,yi+2]之间的平滑曲线被定义为由邻近的四个点i，i+1，i+2，i+3得到的一个三阶的拉格朗日多项式li(x)：

其三阶导数为：

每对数据点之间的三阶导数设定为常量，式(1)可以表示为：

其中

这个平滑的定义表示理想的平滑函数是一个时间的二次函数(其一阶导数是一个线性函数)。

2)平滑曲线对观测值的逼真度：

3)平滑曲线对观测值一阶导数的逼真度：

一阶导数的平滑值可以根据平滑函数值yi表示。

使用用来定义s的拉格朗日多项式li(x)的一阶导数，表示为l′i(x)：

l′i(x)＝ai(x)yi+bi(x)yi+1+ci(x)yi+2+di(x)yi+3(10)

其中，

为用平滑函数值yi表示一阶导数的平滑值对于每个时刻xi的周围四个点的选择有着一定的自由。以下公式的约束保证了的值处于平滑曲线上。

1)为了计算使用时刻x1，x2，x3，x4的值，也就是说可以得到公式(11)

其中，

2)对于输入数据的前半部分的的计算，i＝2，3，...n/2，使用时刻xi-1，xi，xi+1，xi+2，即

其中，

3)对于输入数据的后半部分的的计算，i＝n/2+1，n/2+2，…n-1,使用时刻xi-2，xi-1，xi，xi+1，即

其中，

4)对于的计算，使用时刻xn-3，xn-2，xn-1，xn，即

其中，

则可以表示为：

v-c组合滤波方法旨在确定平滑值yi在三种不同条件下的均衡，每种不同的设置会导致不同的结果：

(a)曲线需要被平滑(最小化s)

(b)平滑值需要接近函数的观测值(最小化f)

(c)平滑曲线的一阶导数需要接近于一阶导数的观测值(最小化)

采用最小二乘法最小化以上约束条件对各参数进行调节，表示为：

其中ε≥0，三种条件的均衡程度可以通过选择这两个参数的值来确定：平滑系数ε和平滑系数值越大，逼真度相对于观测函数值或者其一阶导数的权值越大，平滑值与观测值越接近。

本发明实施例提供一种利用vondark-cepek组合滤波实氢铯时间尺度融合的方法。首次提出对于氢铯钟组时间尺度层面的融合，避免了氢铯互为参考产生的白噪声影响时间尺度结果的问题。氢钟组时间尺度和铯钟组时间尺度分别具有良好的短期和长期稳定性，因此在不影响铯钟组时间尺度长期稳定性情况下，利用氢钟组时间尺度的优良短期稳定性对其进行性能增强，可产生长短稳都良好的融合时间尺度。

如图1所示，产生融合时间尺度的方法包括以下步骤：

步骤1、对实验室测得的原子钟钟差数据进行预处理，采用3σ法则去除异常点。

所述粗差剔除是利用3σ准则，检测的原子钟频差数据为{x1,x2,x3,…,xn-1,xn}，样本均值为标准偏差为3σ准则为若测量值则该值为奇异点，应予以剔除。

步骤2、将铯原子钟和氢原子钟分组，并采用步骤1处理后的钟差数据采用指数滤波方法进行各自钟组时间尺度的计算，得到tacs和tah，如图2所示。

钟差估计：每台守时钟与参考钟的钟差数据xij(t)，测量间隔为τ。在(t+τ)时刻的钟差估计值根据时刻t的钟差和钟速由公式(1)计算。参考钟相对于时间尺度钟差的最优估计则依据加权平均原则，由公式(2)计算。

其中，xi(t),yi(t)分别是钟i在时刻t相对于某个参考时间尺度，对于时差和频差的估计；是钟i在时刻t的钟差预报。

频率估计：钟i在(t+τ)时刻的平均频差估计是采用钟差数据的一次差分，再经指数滤波器计算得到，如公式(3)、(4)。

其中是钟i在时刻t到t+τ的时间间隔上的频差预报；mi是指数频率平均时间常数，是最小化时差估计公式(1)中引入的误差，根据式(5)来估算。

其中τmini为钟i最稳定的时间间隔，取阿伦方差曲线上最小值时τ的值。

权重估计：利用每台钟相对于纸面钟的钟差与当前时刻计算值之差来确定权重。

如公式(6)：

其中ki为误差估计的偏差，根据公式(10)计算。每台钟的均方时间误差是由式(7)的指数滤波器决定。

其中εi(τ)是在时间间隔τ上钟i时差估计的累积误差；<ε²x(τ)>是在时刻t在时间间隔为τ的纸面时均方差，初始值一般取nτ是指数滤波器的时间常数，用来估计当前均方差，nτ一般取20天。

n为原子钟数量。权重在公式(2)中应用，其计算如式(9)所示。

首先读取钟差文件，对xi(t)、yi(t)取初值。对于初值的选取，为了简便选择原子钟在起始时刻相对于utc的钟差值。bipm的t公报发布的是utc与utc(ntsc)的钟差，已知该时刻三台氢钟相对于utc(ntsc)钟差，两差相减，便得到三台氢钟相对于utc的钟差。钟速yi(t)和漂移参数d的初值选择是根据前十天的钟差数据二次拟合得到。其次，根据公式(1)计算下一时刻预测钟差然后计算的初始值，一般取那么首先应计算三台钟的阿伦方差。由于我们测得的是每台种与主钟的钟差，要使用三角帽方法解出每台钟各自的阿伦方差，故需将已知的钟差数据两两相减，去掉主钟，得到三台氢原子钟相互之间的钟差。再计算原子钟两两间阿伦方差，根据三角帽计算公式(11)解出每台钟的阿伦方差，按照公式计算初值。再根据公式(8-9)计算权重初值。然后，根据公式(2)及每台原子钟钟数据计算下一时刻主钟相对于纸面钟的钟差，并进行加权平均，得到最终的钟差估计xj(t+τ)即为计算的平均时间在时刻(t+τ)与参考钟的差，即tacs或tah。接着根据已估计的钟差，进行下一时刻(t+τ)的钟速预报。根据xi＝xj+xij计算每台钟与纸面钟的钟差，代入公式(3)计算频差估计。再根据公式(5)计算出指数频率常数，代入公式(4)计算出yi(t+τ)，以供接下来的时刻(t+2τ)钟差预报的计算。最后，进行接下来时刻(t+2τ)的钟差，权重，频率估计。此时，时差估计与之前计算方法一致，权重计算则根据公式(6)-(10)来计算，先根据公式(10)确定参数ki，再根据公式(6)-(9)计算本次取权。依次循环，再进行下一时刻的预测计算。

步骤3、将步骤2计算得到的氢钟组时间尺度序列进行一次差分，获得对应频率数据。

步骤4、根据所需减弱噪声的频率f或周期p，确定频率响应t/t’，并计算平滑因子ε和

将铯钟组时间尺度序列tacs经傅里叶变换转变到频域，经谱分析确定噪声周期p或频率f，同时确定频率响应t/t’，最终按照公式计算平滑因子。

步骤5、将步骤2计算产生的铯钟组时间尺度和步骤3产生的氢钟组频率序列作为输入，将步骤4获得的平滑因子用于vondark-cepek组合滤波方法，产生时间尺度融合结果。

首先，系统方程以矩阵方式可以表示为

其中，

其中，

由于a是正定矩阵，因此，cholesky分解可用来解系统方程(14)。矩阵a经cholesky分解可得：

a＝u^tdu(15)

令u^tw＝by'，可求得向量w。再根据uy＝d^-1w，可求得y。即：

y＝a^-1by'(16)

因此可求得氢铯融合时间尺度结果即序列y，结果见图3，(a)为时间尺度融合结果，(b)为allan偏差。

以上实施仅为本发明的示例型实施例，不用于限制本发明。本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或同等替换，这种修改或同等替换也应视为落在本发明的保护范围内。