CN103575982A - 一种实时频率稳定度分析方法 - Google Patents

Abstract

本发明提供了一种实时频率稳定度分析方法，依次判别粗大误差并记录状态、数据拟合内插、取样间隔自适应确定和数据迭代，最终得到艾伦方差表征的频率稳定度。本发明能实现频率稳定度实时分析同时，还能实时输出、显示分析结果，简化了用户的工作。本发明还根据频率稳定度分析的特点，设计了数据预处理方法，增强了系统的容错能力，为基于频率计、相位计等常用设备低成本的实现频率标准实时稳定度分析提供解决方案。

Description

一种实时频率稳定度分析方法

技术领域

本发明涉及一种频率稳定度的分析方法。

背景技术

频率标准作为一种电子设备，输出信号不可避免受到内部各种电子噪声的影响，使其频率不是一个固定值，而是在一定范围内随机起伏。频率稳定度就是用来描述频率标准输出频率受噪声影响而产生的随机起伏。所谓时域频率稳定度，就是把受噪声影响的输出频率，用一个时间函数来描述，从时域的角度来分析噪声对输出频率的影响。

输出信号频率是由噪声引起的随机函数，因此只能用统计分析的方法来研究，另外，考虑到信号的瞬时频率是无法直接测量，只能通过测量一段时间的平均频率，这两点决定了时域频率稳定度分析只能采用统计分析方法，分析某取样时间内平均频率的随机起伏情况。

随着对频率标准噪声深入研究发现，频率标准除存在有常见的傅氏频率较高的热噪声和散弹噪声外，还存在着傅氏低频分量很丰富的调频闪变噪声频率随机游动噪声。由于它们的存在，使得在用标准方差分析频率标准的稳定度时，测量次数越多，标准方差越大，理论上当测量次数趋近无穷时，标准方差是发散的。因此，用标准方差来描述频率标准的稳定度是不合适的，因为作为一个有用的统计量，应该是测量次数越多，求出的结果越精确，误差越小。为了解决这一问题，提出了阿伦方差表征方法，其定义式如下：

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)=\frac{1}{2(M-1)}\underset{i=1}{\overset{M-1}{\mathrm{\Σ}}}{[f_{i+1}-f_i]}^2$

其中M表示测量数据组数，τ为取样间隔，f_i表示第i次取样间隔时间的平均频率值。

阿伦方差可以理解为从时域测量噪声强弱的一种方法，其测量原理是将一个不平稳的被测信号通过一个滤波器，使之平稳化，即利用了非平稳过程具有平稳增量的属性，然后从中提取表征噪声强弱的噪声系数。从统计学的角度，阿伦方差是采用时间样本平均的方法，提取所需的噪声信息，由于估计对象是单样本，所以不能采用空间样本平均。采用时间样本平均，必须保证被估计信息在时间上的稳定性，否则便没有空间样本与时间样本的互换性，这就是估计量的平稳遍历性要求，对振荡器噪声，标准方差中样本是每一组测量结果求平均，由于不具有平稳遍历性，每一次测量的统计意义都不一样，所以重复测量没有意义，而阿伦方差是标准方差的平均，由于样本是平稳遍历的，所以可复现测量结果，反映样本的噪声信息。

目前原子频标的频率稳定度分析根据处理时间不同分为离线事后分析和在线实时分析两种，其中离线事后分析较为常见，通过统计大量测试数据，可以实现对任意取样时间的稳定度分析，数据越多，处理耗费时间越长；当前的在线分析为了保证响应速度，只能处理有限长时间内的一段数据，分析其稳定度，取样间隔大小、测量数据组数受计算机的运行速度、存储器的容量限制，只能取某有限长时间内数据，不能实现真正的在线分析。

发明内容

为了克服现有技术的不足，本发明提供一种采用数据递推的方法实现基于频率值、相位差值实现频率稳定度在线分析，并能够实时输出、显示分析结果。

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

步骤1：初始运行时，根据预设经验值判断当前测量值是否为粗大误差；运行稳定后，计算当前测量值之前N个测量值的均值

和方差σ，

其中f_i表示获得的第i个测量值，然后计算当前测量值f_i+1与均值

之差ε_i+1，

检测是否满足ε_i+1＞3σ条件，若满足，则判断f_i+1为粗大误差；若发现粗大误差进入步骤2，否则转入步骤4；

步骤2：将该粗大误差的测量值写入状态日志，同时记录下测量值发生的时刻和通道编号，同时粗大误差计数器数值加1，并计算粗大误差个数所占总测量值比例，得到粗大误差出现概率结果，写入状态日志；

步骤3：采用二次项函数拟合当前测量值之前的N个测量值y_i(t)＝a_it+² _ib+t,c＝1i,2，..得N到常数项c、一次项a、二次项b的值，然后根据二次项函数预测当前理论测量值y_i+1(t)，用y_i+1(t)的值替换步骤1中f_i+1的值；

步骤4：更新M值，M为截止当前所有参与稳定度分析的测量值数据总数量，计算可用的最大取样间隔长度τ_max，取值为等式右边值四舍五入取整数，其中Const为一常数，取大于等于5的数，自适应取小于τ_max取样间隔作为频率稳定度计算基础；

步骤5：计算相邻两个取样间隔τ时间内频率均值之差的平方F_i(τ)＝(f_i+1-f_i)²，然后计算并保持

当新的频率值到来时，仅需计算最新的结果

若有新的取样间隔满足步骤4的取样间隔要求，则需重新计算对应S_i(τ)的值；

步骤6：将步骤4得到的M和τ值，以及步骤5得到的S_i(τ)值代入稳定度分析公式

得到不同取样间隔τ对应的频率稳定度值进一步对该值开平方得到艾伦方差表征的频率稳定度。

本发明的有益效果是：

本发明提出一种高效的频率稳定度实时分析方法，克服现有方法不能同时满足实时稳定度分析、长期稳定度分析的困难，进一步解决频率稳定度实时分析的理论和技术问题，提供了一种采用数据递推的方法，能实现频率稳定度实时分析同时，还能实时输出、显示分析结果，简化了用户的工作。

本发明还根据频率稳定度分析的特点，设计了数据预处理方法，增强了系统的容错能力，为基于频率计、相位计等常用设备低成本的实现频率标准实时稳定度分析提供解决方案。

具体实施方式

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

本发明是一种解决在长期观测时产生海量测量结果数据，通过数据迭代递推方法，降低运算量，能实现稳定度在线分析的方法。

频率稳定度是一项表征频率标准受噪声污染情况的主要技术指标，特别是当频率标准被用来作为系统的基准参考时，频率标准的稳定度关系着整个系统的性能，如卫星导航系统中的原子频率标准为整个系统提供参考时间频率基准信号，若频率标准受噪声影响且不能被及时发现并控制，将导致导航系统的定位、授时出现偏差，严重时能导致系统失效。因此需要对频率标准不同取样间隔的稳定度进行实时分析，对频率标准的短期稳定度实时分析，能帮助及时发现异常干扰，频率长期稳定度实时分析，能准确评估频率标准的频率慢变化性能，为使用者通过附加手段校准频率标准提供依据。本发明主要解决两个方面的关键问题：

1)数据预处理算法，剔除奇异值

考虑到频率稳定度分析对数值波动较为敏感，受个别奇异值影响，可能导致各取样间隔的分析结果偏大，不能真实反映被测信号稳定度；另外奇异值也是被测信号真实状态的反映，如某时刻被环境干扰或是电路的噪声，因此不能简单采用粗大值剔除的方法。

为了解决上述问题，提出了数据预处理算法，预处理分四个步骤：步骤一，获取测量值数据，然后采用3σ法则判断是否为粗大误差，若满足粗大误差判决条件则进入步骤二，将该疑似粗大误差的测量值写入日志，同时记录下发生的时刻、通道编号等该值对应的基本信息，便于用户事后查询；步骤三，计数粗大误差出现次数，计算粗大误差出现概率并提交给用户；步骤四，根据历史测量值以及测量值增量，内插测量值取代粗大误差结果，保证稳定度分析结果的连续性。

2)稳定度实时分析

实时计算稳定度最大的难点在于，当新的测量值到来后，所有参与稳定度计算的测量值需要根据各取样间隔重新计算，当参与计算的数据累积到一定数量后，单次计算所需时间可能远超过获得新测量值的时间，旧值未计算完，新值又堆积到缓冲区中，越积越多直到缓冲区溢出，因此不能实现实时分析结果。

根据频率稳定度分析是通过数据统计实现的特点，采用数据迭代、计算典型取样间隔两种手段，大大降低运算量，使实时频率稳定度分析且不受采样数据长度限制成为可能。

本发明的技术特征在于通过通信端口，实时（实时定义为：以某个时间为周期，定期发生）获取频率计、相位计等专用设备实时测量并输出的某被测频率信号的频率值或相位差值，采用以下步骤处理数据，根据得到的频率值或相位差值，分别实现对被测频率信号的稳定度实时分析，具体步骤如下：

步骤1：粗大误差判别。粗大误差判别分两种情况，一种是初始运行时，由于尚未获得更多测量值，因此根据预设经验值判断当前测量值是否为粗大误差；另一种情形为运行稳定后（至少获得了超过N个测量值），计算当前值之前N个测量值的均值

和方差σ，

其中f_i表示获得的第i个测量值，然后计算最新测量值f_i+1与均值之差ε_i+1，

检测是否满足ε_i+1＞3σ条件，若满足，则判断f_i+1为粗大误差；若发现粗大误差进入步骤2，否则转入步骤4；

步骤2：状态记录。将该粗大误差的测量值写入状态日志，同时记录下测量值发生的时刻、通道编号等该值对应的基本信息，便于用户事后查询，同时粗大误差计数器数值加1，并计算粗大误差个数所占总测量值比例，得到粗大误差出现概率结果，写入状态日志；

步骤3：数据拟合内插。完成步骤2后，启动数据内插计算，以当前测量值之前的N个测量值为基础（当前若不满足N个测量值条件，则运用简单预测模式，以离粗大误差最短时间的测量值替换粗大误差值），采用二次项函数拟合

得到常数项c、一次项a、二次项b的值，然后根据二次项函数预测当前理论测量值y_i+1(t)（认为得到的测量值是粗大误差，且是由非正常因素引起），用y_i+1(t)的值替换步骤1中f_i+1的值；

步骤4：取样间隔自适应确定。为实现实时稳定度分析，最重要的参数是截止当前所有参与稳定度分析的测量值数据总数量M，用来确定可用的取样间隔。因此，每获得一个可用的测量结果，需更新M值，并根据最新的M值计算可用的最大取样间隔长度τ_max，

取值为等式右边值四舍五入取整数，其中Const为一常数，根据经验Const通常取大于等于5的数，以保障稳定度分析结果至少是由超过5组数据统计得到。一旦τ_max确定，则根据频率稳定度取样间隔使用习惯，自适应取小于τ_max取样间隔，如τ＝1,2,4,10,20,100,500,1000...,τ_max，作为频率稳定度计算基础；确定取样间隔后进入步骤5；

步骤5：数据迭代。根据频率稳定度计算公式特点，首先计算相邻两个某取样间隔τ时间内频率均值之差的平方F_i(τ)＝(f_i+1-f_i)²，然后计算

其中M表示获得频率均值的总个数。将步骤4得到的各取样间隔τ对应的S_i(τ)值保存起来，当新的频率值到来时，仅需计算最新的结果更新 $S_{i+1}\left(\mathrm{\τ}\right)=\left(\underset{i=1}{\overset{M-1}{\mathrm{\Σ}}}{F}_i\left(\mathrm{\τ}\right)\right)+F_{i+1}\left(\mathrm{\τ}\right)=S_i\left(\mathrm{\τ}\right)+{(f_{i+2}-f_{i+1})}^2,$ 大大降低了运算量，若有新的取样间隔τ满足步骤4的取样间隔要求，则需重新计算对应S_i(τ)的值；

步骤6：频率稳定度计算，最后将步骤4得到的M和τ值，以及步骤5得到的S_i(τ)值，代入稳定度分析公式，

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)=\frac{1}{2(M-1)}\underset{i=1}{\overset{M-1}{\mathrm{\Σ}}}{[f_{i+1}-f_i]}^2=\frac{1}{2(M-1)}{S}_{i+1}\left(\mathrm{\τ}\right),$ 得到不同取样间隔τ对应的频率稳定度值

进一步对该值开平方得到艾伦方差表征的频率稳定度。

实施例：

当前频率计同时测量8个不同频率标准输出的频率信号，输出结果为实际频率相对于标称频率的相对值，每秒输出1组8个测量值，分别代表8个频率标准的频率值，当输出的频率值超过5组时，自动启动频率稳定度实时分析程序，进入步骤1。

步骤1：粗大误差判别。此时为初始运行，由于尚未获得更多测量值，因此根据经验值判断当前测量值是否为粗大误差，通常设为当前值前一个值的3倍作为粗大误差判别条件，若大于则判别为粗大误差，转入步骤2，否则进入步骤4，此处为了介绍所有情形，假设第6组频率值存在粗大误差，则转入步骤2；

步骤2：状态记录。将该粗大误差的测量值、对应通道编号如“通道2”写入状态日志，同时记录下测量值发生的时刻“20130320185612”，“201303201856”代表2013年3月20日18点56分12秒检测到一个被判断为粗大误差的频率值，同时通道2的粗大误差计数器数值由0变为1，此时通道2总共得到6个频率值，发现1个粗大误差，则粗大误差概率为16.7%，写入状态日志，然后转入步骤3；

步骤3：数据拟合内插。启动数据内插计算，以当前通道2测量值之前的5个测量值为基础，采用二次项函数拟合得到常数项c＝0、一次项a＝0.00182、二次项b＝0.000014，然后根据二次项函数预测当前理论测量值y₆(t)＝8.89219361E-8，用y₆(t)的值替换步骤1中f₆的值，然后转入步骤4；

步骤4：取样间隔自适应确定。当前测量值数据总数量M=6，因此即当前最多能计算取样间隔为1秒的频率稳定度。确定取样间隔后进入步骤5；

步骤5：数据迭代。根据频率稳定度计算公式特点，首先计算相邻两个某取样间隔τ＝1时间内频率均值之差的平方F_i(τ)＝(f_i+1-f_i)²，然后计算

将步骤4得到的各取样间隔τ对应的S_i(τ)值保存起来，当新的频率值到来时，仅需计算最新的结果更新

大大降低了运算量，若有新的取样间隔τ满足步骤4的取样间隔要求，则需重新计算对应S_i(τ)的值；

步骤6：频率稳定度计算，最后将步骤4得到的M=6和τ＝1值，以及步骤5得到的S_i(τ)＝2.62958656E-11值，代入稳定度分析公式，

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)=\frac{1}{2(M-1)}\underset{i=1}{\overset{M-1}{\mathrm{\Σ}}}{[f_{i+1}-f_i]}^2=\frac{1}{2\×5}\×2.62958656E-11=2.62958656E-12,$ 得到τ＝1s时艾伦方差表征的频率稳定度值σ_y(τ)＝1.6216E-6。

各通道同时刻测量值循环执行一次步骤1到步骤5，当第7组测量值到来时，步骤5仅需计算新增的频率值与前一个频率值的差值的平方，将其与上一组得到的S_i(τ)值求和，即S_i+1(τ)＝S_i(τ)+(f_i+2-f_i+1)²，简化了计算

这一过程，特别是当测量值累积较多时，如可用测量结果超过100000个时，频率计仍保持每秒输出1组频率值，普通个人计算机计算1个S_i(τ)耗时可能远超过1秒，若几个原子频标同时测量，则耗费时间成倍的增加，远不能满足对每秒输出的频率值实现实时稳定度分析需求。本方法的提出，使得累积再多的数据也不会大幅度增加运算量，而是一种比较均匀稳定的运算量需求，适合长期的、无人值守的频率标准监测场合，实时可靠的输出频率标准的稳定度分析结果，并能记录或是自动远程报警频率标准发生的异常状态，不需要高性能的服务器，不需要工作人员长时间值守岗位。

Claims (1)

1.一种实时频率稳定度分析方法，其特征在于包括下述步骤：

步骤1：初始运行时，根据预设经验值判断当前测量值是否为粗大误差；运行稳定后，计算当前测量值之前N个测量值的均值

和方差σ，

其中f_i表示获得的第i个测量值，然后计算当前测量值f_i+1与均值之差ε_i+1，

检测是否满足ε_i+1＞3σ条件，若满足，则判断f_i+1为粗大误差；若发现粗大误差进入步骤2，否则转入步骤4；

步骤2：将该粗大误差的测量值写入状态日志，同时记录下测量值发生的时刻和通道编号，同时粗大误差计数器数值加1，并计算粗大误差个数所占总测量值比例，得到粗大误差出现概率结果，写入状态日志；

步骤3：采用二次项函数拟合当前测量值之前的N个测量值y_i(t)＝a_it+² _ib+t,c＝1i,2，..得N到常数项c、一次项a、二次项b的值，然后根据二次项函数预测当前理论测量值y_i+1(t)，用y_i+1(t)的值替换步骤1中f_i+1的值；

步骤4：更新M值，M为截止当前所有参与稳定度分析的测量值数据总数量，计算可用的最大取样间隔长度τ_max，

τ_max取值为等式右边值四舍五入取整数，其中Const为一常数，取大于等于5的数，自适应取小于τ_max取样间隔作为频率稳定度计算基础；

步骤5：计算相邻两个取样间隔τ时间内频率均值之差的平方F_i(τ)＝(f_i+1-f_i)²，然后计算并保持

当新的频率值到来时，仅需计算最新的结果

若有新的取样间隔满足步骤4的取样间隔要求，则需重新计算对应S_i(τ)的值；

步骤6：将步骤4得到的M和τ值，以及步骤5得到的S_i(τ)值代入稳定度分析公式

得到不同取样间隔τ对应的频率稳定度值进一步对该值开平方得到艾伦方差表征的频率稳定度。