CN112597622B - 一种用于检测铯原子钟频率异常的方法、系统和介质 - Google Patents
一种用于检测铯原子钟频率异常的方法、系统和介质 Download PDFInfo
- Publication number
- CN112597622B CN112597622B CN202011083475.3A CN202011083475A CN112597622B CN 112597622 B CN112597622 B CN 112597622B CN 202011083475 A CN202011083475 A CN 202011083475A CN 112597622 B CN112597622 B CN 112597622B
- Authority
- CN
- China
- Prior art keywords
- clock
- frequency
- difference
- cesium
- prediction
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 229910052792 caesium Inorganic materials 0.000 title claims abstract description 64
- TVFDJXOCXUVLDH-UHFFFAOYSA-N caesium atom Chemical compound [Cs] TVFDJXOCXUVLDH-UHFFFAOYSA-N 0.000 title claims abstract description 64
- 238000000034 method Methods 0.000 title claims abstract description 49
- 230000002159 abnormal effect Effects 0.000 title description 4
- 238000001514 detection method Methods 0.000 claims abstract description 46
- 230000005856 abnormality Effects 0.000 claims abstract description 31
- 238000012360 testing method Methods 0.000 claims description 31
- 238000009792 diffusion process Methods 0.000 claims description 11
- 230000005653 Brownian motion process Effects 0.000 claims description 10
- 230000000694 effects Effects 0.000 abstract description 9
- 238000009825 accumulation Methods 0.000 abstract description 6
- 235000008694 Humulus lupulus Nutrition 0.000 abstract description 2
- 230000000875 corresponding effect Effects 0.000 description 27
- 238000010586 diagram Methods 0.000 description 13
- 238000004458 analytical method Methods 0.000 description 9
- 238000005315 distribution function Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000010998 test method Methods 0.000 description 4
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 2
- 229910052739 hydrogen Inorganic materials 0.000 description 2
- 239000001257 hydrogen Substances 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000010183 spectrum analysis Methods 0.000 description 2
- 238000003657 Likelihood-ratio test Methods 0.000 description 1
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 101150060820 Pfas gene Proteins 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 229920011301 perfluoro alkoxyl alkane Polymers 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000009255 platelet function activity Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005295 random walk Methods 0.000 description 1
- 230000035484 reaction time Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Manipulation Of Pulses (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种用于检测铯原子钟频率异常的方法、系统和介质。所述方法包括：步骤S1、获取铯钟模型，基于所述铯钟模型确定钟差预测不确定度；步骤S2、根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及步骤S3、利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。所述方法充分利用频率异常在时差上的积累效应来检测铯原子钟频率异常；给出了检测概率和虚警概率、平均频率跳变幅度之间的函数关系；能够分析和预测不同铯钟在不同幅度频率跳变时的检测情况。
Description
技术领域
本发明涉及信号处理领域，尤其是涉及一种用于检测铯原子钟频率异常的方法、系统和介质。
背景技术
守时试验室和全球导航卫星系统(Global Navigation Satellite System，GNSS)需要利用多台原子钟组成钟组，以建立和保持一个时间基准。守时本质上是预测问题，即需要通过原子钟的历史表现来预测其未来表现。国际上已经达成共识：一个好的钟应该是一个可预测的钟；可预测性的一个重要体现是原子钟信号没有异常。建立一个实时的、可预测的纸面时间，需要在异常发生后尽可能短的时间内检测出异常；根据异常的严重等级，合理降低其权重甚至剔除出钟组，尽可能降低异常对于纸面时间的影响。
频率异常检测的性能是和原子钟的噪声水平高度相关的。典型氢钟的闪烁底(flicker floor，即Allan偏差最小值)一般在 1×10-15～1×10-16量级，而典型铯钟的闪烁底只能达到1×10-14量级。这就导致了相比氢钟，铯钟1×10-14量级的频率异常可能淹没在噪声之中，从而更不容易检测，其对纸面时间的影响也更大。
国内外检测频率异常的研究方法包括：动态Allan方差法、时频谱分析法、假设检验法。其中，动态Allan方差法对微小频率异常的反应时间较长；时频谱分析法比较适用于检测时变信号、或是周期性波动等在频谱上出现明显谱线的异常；相比而言，假设检验法比较适用于频率异常的检测。
假设检验法的核心是设计合理的检验统计量和相应的检测门限值。检测门限值是由检验统计量的统计特性以及虚警概率决定的。进一步，可以推导出检测概率和虚警概率、平均频率跳变幅度之间的函数关系。一般认为检验统计量的统计特性是先验已知的，这在大部分情况下都是可行的。因为对于某一台原子钟，其噪声的统计特性一般不会发生很大的改变。当检验统计量的统计特性未知，可以采用广义似然比(Generalized LikelihoodRatio Test， GLRT)的方法，在检验统计量的分布函数中用最大似然估计结果代替未知参数。因此，GLRT法可以看成是假设检验法的特例。
现有技术只分析了不同PFA的试验效果，没有进一步推导不同平均频率跳变幅度下的检测概率；包括采用GLRT方法在内的大部分假设检验法的检验统计量都是频差，例如采用统计量分别为 Kalman滤波器输出的频差的一步和多步预测误差，而没有充分利用频率异常在时差上的积累效应。
发明内容
本发明的目的在于提供一种用于检测铯原子钟频率异常的方案，以解决现有技术中存在的在检测铯原子钟频率异常的过程中未充分利用频率异常在时差上的积累效应的技术问题；以及检测概率和虚警概率、平均频率跳变幅度之间的函数关系未知的技术问题。
本发明第一方面提供了一种用于检测铯原子钟频率异常的方法，所述方法包括：步骤S1、获取铯钟模型，基于所述铯钟模型确定钟差预测不确定度；步骤S2、根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及步骤S3、利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。
根据本发明第一方面提供的方法，在所述步骤S1中：
利用以下公式来表示所述铯钟模型，
其中，所述x0表示时差的初值，所述y0表示频差的初值，所述W1(t)和所述W2(t)分别表示两个独立的维纳过程，所述σ1表示所述W1(t)的扩散系数，所述σ2表示所述W2(t)的扩散系数；
利用以下公式来确定所述钟差预测不确定度u(tp)，
其中，所述T表示观测间隔，所述tp表示未来预测的时间长度，σ代表观测噪声的强度。
根据本发明第一方面提供的方法，在所述步骤S2中，当tp<<T 时，确定所述最优频差观测间隔。
根据本发明第一方面提供的方法，在所述步骤S3中，通过约束虚警概率来确定不同平均频率跳变幅度下的所述检测概率。
本发明第二方面提供了一种用于检测铯原子钟频率异常的系统，所述系统包括：钟差预测不确定度确定单元，被配置为基于铯钟模型确定钟差预测不确定度；最优频差观测间隔确定单元，被配置为根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及检测概率确定单元，被配置为利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。
根据本发明第二方面提供的系统，利用以下公式来表示所述铯钟模型，
其中，所述x0表示时差的初值，所述y0表示频差的初值，所述W1(t)和所述W2(t)分别表示两个独立的维纳过程，所述σ1表示所述W1(t)的扩散系数，所述σ2表示所述W2(t)的扩散系数。
根据本发明第二方面提供的系统，所述钟差预测不确定度确定单元具体被配置为，利用以下公式来确定所述钟差预测不确定度 u(tp)，
其中，所述T表示观测间隔，所述tp表示未来预测的时间长度，σ代表观测噪声的强度。
根据本发明第二方面提供的系统，所述最优频差观测间隔确定单元具体被配置为，当tp<<T时，确定所述最优频差观测间隔。
根据本发明第二方面提供的系统，所述检测概率确定单元具体被配置为，通过约束虚警概率来确定不同平均频率跳变幅度下的所述检测概率。
本发明第三方面提供了一种存储有指令的非暂时性计算机可读介质，当所述指令由处理器执行时，执行根据本发明第一方面的一种用于检测铯原子钟频率异常的方法。
采用上述技术方案，本发明的技术方案具有如下有益效果：其充分利用了频率异常在时差上的积累效应来检测铯原子钟频率异常；并给出了检测概率和虚警概率、平均频率跳变幅度之间的函数关系；以及能够分析和预测不同铯钟在不同幅度频率跳变时的检测情况。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为根据本发明实施例的用于检测铯原子钟频率异常的方法的流程的示意图；
图2为根据本发明实施例的用于钟差预测的基本原理的示意图；
图3为根据本发明实施例的预测不确定度及其表达式中第一项和第二项的平方根的图示；
图4a为根据本发明实施例的Ya＝3u(tp)/tp时铯钟的预测误差的图示；
图4b为根据本发明实施例的Ya＝4u(tp)/tp时铯钟的预测误差的图示；
图5a为根据本发明实施例的铯钟钟差的图示；
图5b为根据本发明实施例的钟差的一次多项式残差的图示；
图6为根据本发明实施例的预测误差曲线以及3σ预测不确定度曲线；以及
图7为根据本发明实施例的用于检测铯原子钟频率异常的系统的结构图。
具体实施方式
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明第一方面提供了一种用于检测铯原子钟频率异常的方法。
图1为根据本发明实施例的用于检测铯原子钟频率异常的方法的流程的示意图。如图1所示，在步骤S1、获取铯钟模型，基于所述铯钟模型确定钟差预测不确定度；在步骤S2、根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及在步骤S3、利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。
步骤S1、获取铯钟模型，基于所述铯钟模型确定钟差预测不确定度。
在一些实施例中，典型铯钟的时差表示为：
其中，等号右边前两项为一次多项式表示的确定性分量，后两项为随机性分量，即原子钟噪声。X1表示原子钟时差；x0和y0分别表示时差和频差的初值；W1(t)和W2(t)分别表示两个独立的维纳过程，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布；σ1和σ2分别是这两个维纳过程的扩散系数，用于表明噪声的强度；等号右边第三项表示调频白噪声(White Frequency Modulation Noise，WFM)，等号右边第四项表示调频随机游走噪声(Walk Random Frequency Modulation Noise，RWFM)，即W1(t)、W2(t)的积分在状态变量X1上分别表现为WFM和RWFM。
WFM和RWFM都是有色噪声，一般采用Allan方差来表征频率稳定度。Allan方差表达式为：
其中，τ为平滑时间。式(2)等号右边第一项为WFM分量，第二项为RWFM分量，表明WFM和RWFM在对数Allan方差图中的斜率分别为-1和1。
观测钟差可以认为是时差叠加上观测噪声，表示为：
Z(tk)＝X1(tk)+σ·ξ(tk) (3)
其中，Z(tk)是在tk时刻的观测钟差；ξ(tk)是在tk时刻的观测噪声，一般认为是调相白噪声(white phase modulation noise，WPM)；σ代表观测噪声的强度。
图2为根据本发明实施例的用于钟差预测的基本原理的示意图。如图2所示，钟差预测算法是通过过去时间段[t0-T,t0]的观测钟差，预测未来时间段[t0,t0+tp]的钟差。在估计得到了预测起点时刻t0的时差和频差/>之后，在t0+tp时刻的预测值表示为：
预测误差由[预测值-真实值]来表示。由式(1)和式(4)，预测误差表示为：
ε(tp)包含两部分，由参数估计不确定度在未来时间段[t0,t0+tp]的传播引起，/>由未来时间段[t0,t0+tp]的原子钟自身噪声引起。
由于W1(t)和W2(t)是随机过程，对于某个具体的tp值，W1(tp) 和W2(tp)是随机变量，显然从式(1)中估计得到的和/>都是随机变量。所以，ε(tp)、ε1(tp)和ε2(tp)都是随机变量。把它们的不确定度分别记为u(tp)、u1(tp)和u2(tp)。显然，ε(tp)～N(0,u2(tp))，且
u2(tp)＝u1 2(tp)+u2 2(tp)+2cov(ε1(tp),ε2(tp)) (6)
时差估计值可以取预测起始时刻的观测钟差Z(t0)，频差估计值/>一般取首尾两点观测钟差的差分，即：
把时差估计值的不确定度记为频差估计值的不确定度记为/>时差与频差估计值的互相关不确定度记为/>
根据维纳过程的特性，推导得到：
根据式(8)、(9)和(10)，得到：
同理，根据维纳过程的特性，推导得到：
2cov(ε1(tp),ε2(tp))＝-σ2 2·(2t0-T)·tp 2 (13)
把式(11)、(12)和(13)代入式(6)，得到：
预测不确定度u(tp)为式(14)的平方根。1σ预测不确定度和2σ预测不确定度从统计意义上分别代表了预测误差的68％和95％的置信区间。
步骤S2、根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔。
本地的钟差测量不确定度可做到优于0.1ns(σ2＝1×10-20s2)，当预测时间tp不是太短时，其影响基本可以忽略。假设观测噪声为零，把式(2)代入式(14)，得到：
式(15)中，u2(tp)包含两项；对于给定的tp，第二项的值不变；当选择观测间隔T等于Allan方差曲线(一般为V字形)取最小值对应的平滑时间值时，第一项取最小值，于是u2(tp)取最小值。对于(σ1 2,σ2 2,σ2)＝(4.8×10-23s,1.9×10-36s-1,1×10-20s2)的典型铯钟，使第1项取最小值的T约等于100天。
图3为根据本发明实施例的预测不确定度及其表达式中第一项和第二项的平方根的图示。取(σ1 2,σ2 2,σ2)＝(4.8×10-23s, 1.9×10-36s-1,1×10-20s2)，观测间隔T＝100天，对于不同预测时间tp，由式(14)计算得到的u(tp)(三角形)，以及前三项的平方根(圆圈)和第4项(即式(15)的第2项)的平方根σy(tp)·tp(方块)。
从图3中可知：预测时间较短即当tp<<T时，u(tp)仅略大于σy(tp)·tp；这时预测性能主要由原子钟噪声决定。
从直观上理解：假设时差和频差的估计误差都一直为零，这相当于估计值和真实值完全相同(实际上是不可能的)，这时按照式(5)，如果忽略观测噪声，预测误差完全就是原子钟噪声，噪声方差随时间变大。
进一步，表1列出了当(σ1 2,σ2 2,σ2)＝(4.8×10-23s,1.9×10-36s-1, 1×10-20s2)时，按照式(14)计算的不同T和tp对应的预测不确定度 u(tp)、式(14)中第四项的平方根σy(tp)·tp。
表1不同T和tp对应的u(tp)(单位：s)
| T | tp | u(tp) | σy(tp)·tp |
| 10天 | 4小时 | 8.44×10-10 | 8.31×10-10 |
| 20天 | 4小时 | 8.41×10-10 | 8.31×10-10 |
| 30天 | 4小时 | 8.40×10-10 | 8.31×10-10 |
| 10天 | 1天 | 2.14×10-9 | 2.04×10-9 |
| 20天 | 1天 | 2.09×10-9 | 2.04×10-9 |
| 30天 | 1天 | 2.08×10-9 | 2.04×10-9 |
从表1可以看出，当tp<<T时，即tp相对而言影响较小的情况下，u(tp)只是略大于σy(tp)·tp。
步骤S3、利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。
基本思路
通过对比分析预测误差和预测不确定度，提供了一种检测原子钟频率异常的方法，即当异常发生在未来时间段[t0,t0+tp]时，在某个t时刻(t>t0)及之后时间段的预测误差一直大于3σ预测不确定度的理论值，认为在t时刻频率发生了跳变。因为正常情况下，预测误差落在±1σ、2σ和3σ预测不确定度范围内的概率分别为 68.26％、95.44％和99.74％；预测误差大于3σ预测不确定度的概率很小。
频率异常检测方法及检测概率的表达式
本方法等价于做如下二元假设检验，H0：未发生频率跳变； H1：发生频率跳变。选取预测误差ε(tp)作为检验统计量。首先需要确定在这两种假设下，检验统计量的数学分布。H0情况下，显然预测误差服从N(0,u2(tp))分布。H1情况下，假设跳变发生在[t0,t0+tp]的中间时刻，记为t0+tp-Tjump时刻，其中0<Tjump<tp；跳变幅度为 Y0，显然，在[t0,t0+tp]时间段内的平均频率跳变幅度为：
所以，在[t0,t0+tp]时间段内积累的时差为：
xjump(T)＝Ya·tp＝Y0·Tjump (17)
这时预测误差服从N(Ya·tp,u2(tp))分布。
综上，在这两种假设下统计量的数学分布为：
把H0为真时判H1成立的概率称为虚警概率，用PFA表示。把 H1为真时判H1成立的概率称为检测概率，用PD表示。PFA的值也被称为显著性水平。
本方法通过约束虚警概率PFA，求解不同的平均频率跳变幅度 Ya对应的检测概率PD。
如前文所述，ε(tp)落在±u(tp)、±2u(tp)和±3u(tp)范围内的概率分别为68.26％、95.44％和99.74％，对应的PFA分别为 1-68.26％＝31.74％、1-95.44％＝4.56％和1-99.74％＝0.26％。
以约束PFA＝0.26％为例。直观上理解，假如|Ya|·tp＝3u(tp)，显然有PD＝50％；那么当|Ya|·tp>3u(tp)时，显然PD>50％；当|Ya|·tp<3u(tp) 时，显然PD<50％。因此，当约束PFA时，PD和|Ya|正相关；|Ya|越大，PD越大。根据式(18)，推导出检测概率的表达式为
由式(19)，可以求解得到任意幅度Ya对应的PD。
以(σ1 2,σ2 2)＝(4.8×10-23s,1.9×10-36s-1)的典型铯钟为例，取σ2＝1×10-20s2，按照T＝20天和tp＝1天进行钟差预测，u(tp)的计算结果对应表1倒数第2行。约束PFA＝0.26％，即式(21)的积分下限 (上式)或上限(下式)分别为3u(tp)和-3u(tp)，然后，根据式(21)，计算不同|Ya|值对应的PD值。
进一步，按照概率论与数理统计的相关理论，服从标准正态分布的随机变量X的分布函数用φ(x)来表示。几个关键数值为φ(-3)＝0.0013，φ(-2)＝0.0228，φ(-1)＝0.1587，φ(0)＝0.5，φ(1)＝ 0.8413，φ(2)＝0.9772，φ(3)＝0.9987。
前文所述的“ε(tp)落在±u(tp)、±2u(tp)和±3u(tp)范围内的概率分别为68.26％、95.44％和99.74％”，对应于φ(1)-φ(-1)＝68.26％、φ(2)-φ(-2)＝95.44％和φ(3)-φ(-3)＝99.74％。
显然，在H1情况下，统计量ε(tp)的分布函数为
所以，对于该实例，在约束PFA＝0.26％时，不同|Ya|对应的PD如表2所示。
表2不同|Ya|对应的PD
由表2看出：当|Ya|>4u(tp)/tp时，可以比较有效地检测出异常；当|Ya|>1.05×10-13时，检测概率在90％以上；当|Ya|较小(例如<3u(tp)/tp)时，检测效果并不显著。
更一般的情况
约束PFA＝0.26％，对应的ε(tp)门限设置为±3u(tp)。在更一般的情况下，可以约束PFA＝5％、PFA＝1％、PFA＝0.1％等。这时，对应的ε(tp)门限值±γ需要查询分布函数φ(x)的函数表来确定，具体方式是使φ(γ)-φ(-γ)＝1-PFA。例如，当约束PFA＝5％时，查阅φ(x)的函数表，得到γ＝1.96u(tp)。然后，需要将式(19)中积分下限(上式)或上限 (下式)分别修改为γ和-γ，得到式(20)；根据式(20)，可以求解得到任意幅度|Ya|对应的PD值。
示例分析
一、仿真试验
仿真生成10000台参数为(σ1 2,σ2 2,σ2)＝(4.8×10-23s,1.9×10-36s-1, 1×10-20s2)的典型铯钟。按照T＝20天和tp＝1天进行钟差预测。选取γ＝3u(tp)，即约束PFA＝0.26％。分别设置式(18)中Y0＝6u(tp)/tp且 Tjump＝tp/2(即在t0+12小时时刻发生频率跳变)，计算得到 Ya＝3u(tp)/tp。然后重复上述步骤，再次仿真生成同样参数的另10000 台典型铯钟。设置Y0＝8u(tp)/tp且Tjump＝tp/2，计算得到Ya＝4u(tp)/tp。按照表2，理论上这两种情况下PD分别为50％和84.13％。
图4a为根据本发明实施例的Ya＝3u(tp)/tp时铯钟的预测误差的图示；图4b为根据本发明实施例的Ya＝4u(tp)/tp时铯钟的预测误差的图示；从图4a和4b中明显看出在t0+tp/2时刻发生的频率跳变。
分别统计在这两种情况下ε(tp)落在±3u(tp)外的铯钟数量，作为本方法检测出频率跳变异常的数量。在这两种情况下，能检测出频率跳变异常的数量分别为5067台和8502台。
表3列出了在这两种情况下，理论PD值、实际检测数量和试验PD值。
表3不同Ya对应的预测不确定度(单位：s)
从表3看出，实际PD值和理论PD值基本符合。
综上，仿真试验验证了理论分析的结论。
二、实测试验
采用某1台5071A铯钟的实测数据进行分析，其参数为(σ1 2,σ2 2, σ2)＝(4×10-23s,6×10-37s-1,1.6×10-21s2)。图5a为根据本发明实施例的铯钟钟差的图示；图5b为根据本发明实施例的钟差的一次多项式残差的图示；长度共46天。
人为在第39.5天处引入幅度为8×10-14的频率跳变。共进行20 次预测，都选取T＝20天，第1次预测从第20天开始，每次预测开始时间往后移1天，最后1次预测开始时间为第39天。每次预测时间选取tp＝3天。共获得20组预测误差曲线。选取γ＝3u(tp)，即约束PFA＝0.26％。
图6为根据本发明实施例的预测误差曲线以及3σ预测不确定度曲线，其中第20组标注圆圈，第19组为标注方块，其余18组为细线无标注，3σ预测不确定度曲线为粗线。从图6看出：对于第 20组，相当于在tp＝0.5天时发生了异常，当tp＝1天时未检测出异常，tp＝2天和tp＝3天时成功检测出了异常；对于第19次预测，相当于在tp＝1.5天时发生了异常，当tp＝2天时未检测出异常，当tp＝3 天时成功检测出了异常。
表4列出了对于第20组和第19组，当tp分别等于1天、2天和3天时，按照式(16)计算得到的γ、按照式(18)计算得到的 Ya、按照式(23)计算得到的理论PD值、实际检测情况。
表4不同Ya对应的预测不确定度(单位：s)
综上，实测数据试验验证了理论分析的结论。
三、试验结果分析
进一步，根据本发明的理论分析和仿真试验，还可以进一步得到以下结论。
(1)PD和PFA、|Ya|之间存在函数关系，约束PFA、|Ya|之中任意一个量的值，可以得到PD和另一个量的函数曲线。
当约束PFA时，可以画出了PD-|Ya|曲线。当约束|Ya|时，还可以画出PD-PFA曲线，这就是著名的接收机工作特性(receiver operating characteristic，ROC)曲线，分析在不同|Ya|情况下不同 PFA值对应的PD值，可以更合理地选取PFA的值。
(2)根据ROC曲线，对于相同的|Ya|，假如PFA值更大，相应的PD也更大。
例如，当T＝20天和tp＝1天时，根据表2，假如把PFA值由0.26％提高至4.56％，相当于把检测门限γ由±3u(tp)降低至±2u(tp)，原来 |Ya|＝3u(tp)/tp＝7.26×10-14时，才有PD＝50％，现在只需要 |Ya|＝2u(tp)/tp＝4.84×10-14时，就有PD＝50％，而当|Ya|依然为3u(tp)/tp＝7.26×10-14时，PD＝84.13％。
(3)当tp<<T时，对于相同的PFA，假设|Ya|保持不变，增加 tp，可以有效增大PD。
例如：当T＝20天，约束PFA＝0.26％时，如果tp＝1天，当Ya＝3u(1 天)/1天＝7.26×10-14时，有PD＝50％；如果tp＝2天，当Ya依然为 7.26×10-14时，根据式(21)计算，有PD＝87.36％；如果tp＝2天，此时只要Ya＝3u(2天)/2天＝5.26×10-14，即有PD＝50％。从直观上理解：预测时间tp越长，由频率跳变引起的积累时差越大，即式(20) 中的|Ya|·tp变大，所以PD增大。
(4)进一步，当tp<<T时，对于同样的Y0，增加tp，使得|Ya| 也相应地变大，|Ya|和tp同时增大使PD增大的效果更显著。
例如，依然采用上述情况，假设在t0+12小时时刻发生频率跳变，幅度为Y0＝6u(1天)/1天＝1.45×10-13，如果tp＝1天，按照(18) 计算得到Ya＝Y0/2＝7.26×10-14，这时根据式(21)计算，有PD＝50％；如果tp＝2天，按照(18)计算得到Ya＝3Y0/4＝1.14×10-13，这时根据式(21)计算，有PD＝99.93％，即基本可以检测出异常。所以，对于同样的Y0，增加tp，相当于使式(20)中的|Ya|和tp都增大，所以|Ya|·tp将变得更大，从而PD更高。
综上，本发明的第一方面提出了一种基于钟差预测的铯原子钟频率异常的检测算法。采用假设检验法时，检测概率(PD)和虚警概率(PFA)、平均频率跳变幅度(|Ya|)之间存在函数关系。为了充分利用频率异常体现在时差上的积累效应，本方法把预测误差作为检验统计量，对检验统计量做二元假设检验。在两种情况下，检验统计量都服从正态分布。通过约束PFA，查询正态分布的分布函数，得到对应的检测门限值。根据在异常情况下的检验统计量分布函数和检测门限值，推导给出了求解|Ya|对应的PD的表达式。仿真试验和实测数据试验都验证了理论分析的结论。
在试验中，由于算法已经约束了PFA，一旦算法给出报警，那么极有可能发生了频率异常(发生异常的概率为1-PFA)，此时需要进行人工干预。|Ya|越大，PD也就越高。根据本发明的分析结论，提高PD的方法有：(1)提高PFA(相当于降低γ，尽管会增加系统报警的频度)；(2)增加tp(当tp<<T时，相当于同时增大了|Ya|和 tp)。
本发明第二方面提供了一种用于检测铯原子钟频率异常的系统，图7为根据本发明实施例的用于检测铯原子钟频率异常的系统的结构图。如图7所示，系统700包括：钟差预测不确定度确定单元701，被配置为基于铯钟模型确定钟差预测不确定度；最优频差观测间隔确定单元702，被配置为根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及检测概率确定单元703，被配置为利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率。
在一些实施例中，利用以下公式来表示所述铯钟模型，
其中，所述x0表示时差的初值，所述y0表示频差的初值，所述W1(t)和所述W2(t)分别表示两个独立的维纳过程，所述σ1表示所述W1(t)的扩散系数，所述σ2表示所述W2(t)的扩散系数。
在一些实施例中，所述钟差预测不确定度确定单元701具体被配置为，利用以下公式来确定所述钟差预测不确定度u(tp)，
/>
其中，所述T表示观测间隔，所述tp表示未来预测的时间长度，σ代表观测噪声的强度。
在一些实施例中，所述最优频差观测间隔确定单元702具体被配置为，当tp<<T时，确定所述最优频差观测间隔。
在一些实施例中，所述检测概率确定单元具体703被配置为，通过约束虚警概率来确定不同平均频率跳变幅度下的所述检测概率。
本发明第三方面提供了一种存储有指令的非暂时性计算机可读介质，当所述指令由处理器执行时，执行根据本发明第一方面的一种用于检测铯原子钟频率异常的方法。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (8)
1.一种用于检测铯原子钟频率异常的方法，其特征在于，所述方法包括：
步骤S1、获取铯钟模型，基于所述铯钟模型确定钟差预测不确定度；
步骤S2、根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及
步骤S3、利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率；
在所述步骤S3中，通过约束虚警概率来确定不同平均频率跳变幅度下的所述检测概率；其中：
进行二元假设检验，H0表示未发生频率跳变；H1表示发生频率跳变，选取预测误差ε(tp)作为检验统计量；
确定在所述二元假设下检验统计量的数学分布：
H0情况下，预测误差服从N(0,u2(tp))分布；
H1情况下，跳变发生在[t0,t0+tp]的中间时刻，记为t0+tp-Tjump时刻，0<Tjump<tp，跳变幅度为Y0，在[t0,t0+tp]时间段内的平均频率跳变幅度为：
则在[t0,t0+tp]时间段内积累的时差为：
xjump(T)＝Ya·tp＝Y0·Tjump
预测误差服从N(Ya·tp,u2(tp))分布；则
在所述二元假设下统计量的数学分布为：
将H0为真时判H1成立的概率作为虚警概率PFA，将H1为真时判H1成立的概率作为检测概率PD。
2.根据权利要求1所述的用于检测铯原子钟频率异常的方法，其特征在于，在所述步骤S1中：
利用以下公式来表示所述铯钟模型，
其中，所述x0表示时差的初值，所述y0表示频差的初值，所述W1(t)和所述W2(t)分别表示两个独立的维纳过程，所述σ1表示所述W1(t)的扩散系数，所述σ2表示所述W2(t)的扩散系数；
利用以下公式来确定所述钟差预测不确定度u(tp)，
其中，所述T表示观测间隔，所述tp表示未来预测的时间长度，σ代表观测噪声的强度。
3.根据权利要求1所述的用于检测铯原子钟频率异常的方法，其特征在于，在所述步骤S2中，当tp<<T时，确定所述最优频差观测间隔。
4.一种用于检测铯原子钟频率异常的系统，其特征在于，所述系统包括：
钟差预测不确定度确定单元，被配置为基于铯钟模型确定钟差预测不确定度；
最优频差观测间隔确定单元，被配置为根据所述钟差预测不确定度，确定钟差预测过程中的最优频差观测间隔；以及
检测概率确定单元，被配置为利用所述最优频差观测间隔对应的钟差预测不确定度，确定用于检测所述铯原子钟频率异常的检测概率；
所述检测概率确定单元具体被配置为，通过约束虚警概率来确定不同平均频率跳变幅度下的所述检测概率；其中：
进行二元假设检验，H0表示未发生频率跳变；H1表示发生频率跳变，选取预测误差ε(tp)作为检验统计量；
确定在所述二元假设下检验统计量的数学分布：
H0情况下，预测误差服从N(0,u2(tp))分布；
H1情况下，跳变发生在[t0,t0+tp]的中间时刻，记为t0+tp-Tjump时刻，0<Tjump<tp，跳变幅度为Y0，在[t0,t0+tp]时间段内的平均频率跳变幅度为：
则在[t0,t0+tp]时间段内积累的时差为：
xjump(T)＝Ya·tp＝Y0·Tjump
预测误差服从N(Ya·tp,u2(tp))分布；则
在所述二元假设下统计量的数学分布为：
将H0为真时判H1成立的概率作为虚警概率PFA，将H1为真时判H1成立的概率作为检测概率PD。
5.根据权利要求4所述的用于检测铯原子钟频率异常的系统，其特征在于，利用以下公式来表示所述铯钟模型，
其中，所述x0表示时差的初值，所述y0表示频差的初值，所述W1(t)和所述W2(t)分别表示两个独立的维纳过程，所述σ1表示所述W1(t)的扩散系数，所述σ2表示所述W2(t)的扩散系数。
6.根据权利要求4所述的用于检测铯原子钟频率异常的系统，其特征在于，所述钟差预测不确定度确定单元具体被配置为，利用以下公式来确定所述钟差预测不确定度u(tp)，
其中，所述T表示观测间隔，所述tp表示未来预测的时间长度，σ代表观测噪声的强度。
7.根据权利要求4所述的用于检测铯原子钟频率异常的系统，其特征在于，所述最优频差观测间隔确定单元具体被配置为，当tp<<T时，确定所述最优频差观测间隔。
8.一种存储有指令的非暂时性计算机可读介质，当所述指令由处理器执行时，执行根据权利要求1-3中任一项所述的用于检测铯原子钟频率异常的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011083475.3A CN112597622B (zh) | 2020-10-12 | 2020-10-12 | 一种用于检测铯原子钟频率异常的方法、系统和介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011083475.3A CN112597622B (zh) | 2020-10-12 | 2020-10-12 | 一种用于检测铯原子钟频率异常的方法、系统和介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112597622A CN112597622A (zh) | 2021-04-02 |
| CN112597622B true CN112597622B (zh) | 2024-01-19 |
Family
ID=75180367
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011083475.3A Active CN112597622B (zh) | 2020-10-12 | 2020-10-12 | 一种用于检测铯原子钟频率异常的方法、系统和介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112597622B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007178363A (ja) * | 2005-12-28 | 2007-07-12 | Seiko Precision Inc | 計時システム、時刻サービス提供装置及び時刻サービス提供システム |
| CN104199058A (zh) * | 2014-09-18 | 2014-12-10 | 中国人民解放军国防科学技术大学 | 基于Kalman滤波器实时预测值调整时间尺度算法 |
| KR20150009125A (ko) * | 2013-07-15 | 2015-01-26 | 한국과학기술원 | 인지 무선에서의 광대역 주파수 검출 방법 및 장치 |
| CN106773610A (zh) * | 2016-11-28 | 2017-05-31 | 北京工业大学 | 一种铯原子钟和氢钟频差预估方法 |
| CN107086901A (zh) * | 2017-05-11 | 2017-08-22 | 中国人民解放军国防科学技术大学 | 一种bdt建立方法及utc(ntsc)建立方法 |
| CN109508510A (zh) * | 2018-12-20 | 2019-03-22 | 国网河南省电力公司焦作供电公司 | 一种基于改进的卡尔曼滤波的铷原子钟参数估计算法 |
| CN109902882A (zh) * | 2019-03-21 | 2019-06-18 | 北京工业大学 | 原子钟钟差预测模型训练方法及装置 |
| CN110687555A (zh) * | 2019-09-23 | 2020-01-14 | 西安空间无线电技术研究所 | 一种导航卫星原子钟弱频率跳变在轨自主快速检测方法 |
| CN111143989A (zh) * | 2019-12-25 | 2020-05-12 | 北京无线电计量测试研究所 | 频率调整量计算方法、模块、系统、存储介质和设备 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9826338B2 (en) * | 2014-11-18 | 2017-11-21 | Prophecy Sensorlytics Llc | IoT-enabled process control and predective maintenance using machine wearables |
| JP6584662B2 (ja) * | 2016-06-07 | 2019-10-02 | 三菱電機株式会社 | 異常診断装置及び異常診断方法 |
-
2020
- 2020-10-12 CN CN202011083475.3A patent/CN112597622B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007178363A (ja) * | 2005-12-28 | 2007-07-12 | Seiko Precision Inc | 計時システム、時刻サービス提供装置及び時刻サービス提供システム |
| KR20150009125A (ko) * | 2013-07-15 | 2015-01-26 | 한국과학기술원 | 인지 무선에서의 광대역 주파수 검출 방법 및 장치 |
| CN104199058A (zh) * | 2014-09-18 | 2014-12-10 | 中国人民解放军国防科学技术大学 | 基于Kalman滤波器实时预测值调整时间尺度算法 |
| CN106773610A (zh) * | 2016-11-28 | 2017-05-31 | 北京工业大学 | 一种铯原子钟和氢钟频差预估方法 |
| CN107086901A (zh) * | 2017-05-11 | 2017-08-22 | 中国人民解放军国防科学技术大学 | 一种bdt建立方法及utc(ntsc)建立方法 |
| CN109508510A (zh) * | 2018-12-20 | 2019-03-22 | 国网河南省电力公司焦作供电公司 | 一种基于改进的卡尔曼滤波的铷原子钟参数估计算法 |
| CN109902882A (zh) * | 2019-03-21 | 2019-06-18 | 北京工业大学 | 原子钟钟差预测模型训练方法及装置 |
| CN110687555A (zh) * | 2019-09-23 | 2020-01-14 | 西安空间无线电技术研究所 | 一种导航卫星原子钟弱频率跳变在轨自主快速检测方法 |
| CN111143989A (zh) * | 2019-12-25 | 2020-05-12 | 北京无线电计量测试研究所 | 频率调整量计算方法、模块、系统、存储介质和设备 |
Non-Patent Citations (5)
| Title |
|---|
| "Optimal Observation Intervals for Clock Prediction Based on the Mathematical Model Method;Yiwei Wu, Xiangwei Zhu, Yangbo Huang, Guangfu Sun, Gang Ou;《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》;第65卷(第1期);第132-143页 * |
| Uncertainty Derivation and Performance Analyses of Clock Prediction Based on Mathematical Model Method;Yiwei Wu, Xiangwei Zhu, Yangbo Huang, Guangfu Sun, and Gang Ou;《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》;第64卷(第10期);全文 * |
| 原子钟模型和频率稳定度分析方法;伍贻威，杨斌，肖胜红，王茂磊;《武汉大学学报· 信息科学版》;第44卷(第8期);全文 * |
| 导航卫星星载原子钟异常监测分析;牛飞;韩春好;张义生;常守峰;;武汉大学学报(信息科学版)(05);全文 * |
| 星载原子钟频率跳变检测方法仿真分析研究;秦晓伟;孙云峰;杜二旺;王国永;;时间频率学报(02);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112597622A (zh) | 2021-04-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP5572877B2 (ja) | 異常な擬似距離測定値から無線ナビゲーション受信機ユーザを保護するための方法 | |
| Betken | Testing for change‐points in long‐range dependent time series by means of a self‐normalized Wilcoxon test | |
| Galleani et al. | Detection of atomic clock frequency jumps with the Kalman filter | |
| Nunzi et al. | Detection of anomalies in the behavior of atomic clocks | |
| CN110687555B (zh) | 一种导航卫星原子钟弱频率跳变在轨自主快速检测方法 | |
| Galleani | The dynamic Allan variance III: Confidence and detection surfaces | |
| Tavella | Statistical and mathematical tools for atomic clocks | |
| Aliev et al. | Algorithms for indicating the beginning of accidents based on the estimate of the density distribution function of the noise of technological parameters | |
| CN112597622B (zh) | 一种用于检测铯原子钟频率异常的方法、系统和介质 | |
| Huang et al. | Detection of weak frequency jumps for GNSS onboard clocks | |
| He et al. | Sequential Bayesian planning for accelerated degradation tests considering sensor degradation | |
| JP2011044592A (ja) | 信頼度判断装置、信頼度判断方法、及び信頼度判断用コンピュータプログラム | |
| Betken et al. | Subsampling for general statistics under long range dependence with application to change point analysis | |
| Dogu | Change point estimation based statistical monitoring with variable time between events (TBE) control charts | |
| Nesterov | Online Prediction of COVID19 Dynamics: Belgian Case Study | |
| Choudhury et al. | Linear or nonlinear? A bicoherence based metric of nonlinearity measure | |
| CN112581727A (zh) | 桥梁的位移漂移预警方法、装置、设备及存储介质 | |
| Betken et al. | Subsampling for General Statistics under Long Range Dependence with application to change point analysis | |
| Das et al. | Data driven approach for performance assessment of linear and nonlinear Kalman filters | |
| Nunzi et al. | The generalized likelihood ratio test for detecting anomalous behaviors of atomic clocks | |
| Abbasi et al. | Monitoring coefficient of variation using progressive mean technique | |
| Chernoyarov | The efficiency of reception of a random pulse signal with unknown parameters under the conditions of detuned duration | |
| Cao et al. | Electromechanical mode estimation validation using recursive residual whiteness testing | |
| Mařík | Thresholding Using Extreme Value Theory Threshold Models | |
| de Carvalho | The establishment of a Brazilian atomic time scale |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |