删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

原子钟噪声变化时改进的Kalman滤波时间尺度算法

本站小编 Free考研考试/2021-12-29
摘要:Kalman滤波时间尺度算法是一种实时的原子钟状态估计方法, 在时间保持工作中具有重要的实用价值. 本文引入自适应因子改进Kalman滤波时间尺度算法, 对Kalman滤波时间尺度算法中状态预测协方差矩阵引入自适应因子, 利用统计量实时计算自适应因子的量值, 控制状态预测协方差矩阵的增长, 降低原子钟状态估计的扰动, 从而提高时间尺度的准确度和稳定度. 模拟数据和实测数据表明, 原子钟噪声强度变化时或噪声强度估计存在误差时改进的Kalman滤波时间尺度算法有效地提高了时间尺度的准确度和长期稳定度.
关键词: 原子钟噪声/
Kalman滤波/
自适应因子/
时间尺度

English Abstract


全文HTML

--> --> -->
1.引　言
Kalman滤波作为一种经典的估计技术, 广泛应用于控制理论、信号分析和卫星导航星载原子钟的异常检测等领域[1,2]. 作为最小二乘方法的动态形式, Kalman滤波是一种实时算法, 通过一组观测序列及系统的动力学模型来求解状态向量估值. 这个估计值是通过混合两种信息得到的, 分别为基于动力学模型的预测信息和基于测量方程的新息向量. 最后通过对这两部分信息进行加权平均得到Kalman滤波估计. 当动力学模型噪声和测量方程噪声服从高斯分布时, 那么在最小二乘意义下, 估计是线性最优的, 在其他情况下, Kalman滤波算法是线性次优的. 除此之外, 为了保证估计的无偏性, 每一种噪声均值为0. 在时间尺度计算中, 要对原子钟的噪声强度进行估计. 文献[3-5]研究了Kalman滤波算法中噪声模型及协方差矩阵的确定方法. 文献[6-11]研究了基于Kalman滤波的时间尺度算法, 其中包含原子钟噪声模型的建立. 文献[6]基于原子钟模型研究了三类Kalman滤波时间尺度, 并分析了降低协方差算法的最优性. 其中Kalman加权时间尺度算法需要考虑原子钟的基本时间尺度方程(BTSE), 考虑了原子钟的频率状态估计和频漂状态估计, 权重算法考虑了时间尺度的短期稳定度. 本文通过Kalman滤波估计原子钟相位状态, 并根据“可预测性”权重算法计算时间尺度, 保证了时间尺度长期稳定度前提下, 改善其短期稳定度. 文献[7]基于氢原子钟研究了噪声类型主要为调频闪变噪声的一种Kalman滤波时间尺度算法. Kalman滤波算法中调频闪变噪声模拟为马氏噪声过程的线性组合. 本文考虑原子钟噪声类型主要为调频白噪声、调频随机游走噪声和频漂随机游走噪声. 文献[7]通过对状态估计误差协方差矩阵进行改正, 提高了Kalman滤波时间尺度的短期稳定度. 本文是在保证时间尺度长期稳定度的基础上, 改进时间尺度的短期稳定度. 文献[8]基于Kalman滤波算法研究了远距离不同实验室比对噪声情况下的时间尺度算法, 是对同一实验室无比对噪声情况的一种扩展. 文献[9]研究了Kalman滤波时间尺度算法中组合Brown和Greenhall 的协方差降低运算, 进一步优化了时间尺度计算结果, 提高了时间尺度的短期稳定度. 文献[10]探讨了将Kalman滤波算法用于计算协调世界时(coordinated universal time, UTC), 并与目前原子时尺度算法做了比较, 表示Kalman时间尺度具有较高的短期稳定度. 文献[11]研究了使用Kalman滤波器调整预测值的时间尺度算法. 文中建立原子钟的两状态方程, 通过Kalman滤波器估计原子钟的状态, 根据状态预测值计算时间尺度. 本文基于原子钟的三状态模型, 通过Kalman滤波器估计原子钟状态, 并根据原子钟的“可预测性”权重算法计算时间尺度, 保证时间尺度长期稳定度, 同时根据Kalman滤波算法改进时间尺度的短期稳定度. 以上文献研究的前提是原子钟噪声强度不变并且噪声强度能够准确估计, 但是外界环境影响与原子钟长期运行导致其性能下降, 原子钟噪声可能发生变化或噪声强度不能够准确估计, 这时根据Kalman滤波时间尺度算法产生的时间尺度性能降低. 文献[12,13]分别是基于渐消因子的改进Kalman滤波时间尺度估计算法研究和Sage窗的自适应Kalman滤波用于钟差预报的研究. 文献[12]对状态预测协方差矩阵引入渐消因子, 降低原子钟状态估计的扰动, 通过原子钟的时差状态和时差状态估计之差定义时间尺度, 提高了时间尺度的短期稳定度. 而本文研究自适应因子改进Kalman滤波状态估计, 原子钟的状态估计作为原子钟的预测部分进行扣除, 基于“可预测性”取权方法计算时间尺度. 这样既保证了时间尺度的长期稳定度, 同时又改善了时间尺度的短期稳定度. 文献[13]研究了基于Sage窗的自适应Kalman滤波用于钟差预报, 基于Sage窗构造自适应因子, 取决于窗口内各历元状态参数与当前状态的适应程度, 若所开的窗口内的状态参数改正量正好反映载体当前状态扰动水平, 则求得的状态方程协方差矩阵能调节状态方程的噪声水平; 反之由于对历史信息的平滑, 反而使状态参数偏离. 本文时间尺度算法采用预测残差构造自适应因子, 能够较好地反映模型变化, 但对于稳定状态, 基于Sage窗的自适应Kalman滤波具有更好的效果. 鉴于此, 本文研究了噪声强度变化时改进的Kalman滤波时间尺度算法, 通过模拟数据和实测数据验证并分析了该算法的性能.
Kalman 滤波时间尺度算法需要确定原子钟状态模型噪声, 通常情况下模型噪声由调频白噪声、调频随机游走噪声和频漂随机游走噪声组成, 并且假定原子钟噪声强度恒定, 但是受环境因素或原子钟老化影响, 原子钟噪声强度会发生改变. Kalman滤波时间尺度算法需要估计噪声强度, 估计方法的准确与否会影响时间尺度的性能, 通常根据原子钟测得的Allan 方差和Hadamard 方差来估计原子钟噪声强度, 但是受测量条件的影响, 原子钟各类方差估计存在误差, 因此原子钟噪声强度估计存在偏差, 最终影响时间尺度的性能. 文中改进Kalman滤波时间尺度算法的主要思想是原子钟噪声强度变化或噪声强度估计存在误差时, 计算自适应因子控制Kalman滤波预测协方差矩阵的增长, 降低状态模型信息不准确对状态估计的影响. 对于原子钟噪声强度恒定与噪声强度准确估计时, 文中改进的Kalman滤波时间尺度算法变为经典Kalman滤波时间尺度算法.
文章组织如下: 第2节建立基于Kalman滤波算法的原子钟组模型; 第3节推导改进的Kalman滤波状态估计; 第4节介绍“可预测性”取权方法; 第5节为时间尺度的计算及分析; 第6节为本文的结论.
2.钟组模型
在具有$N$台原子钟的钟组中, 状态向量包含每一台原子钟的相位偏差, 频率偏差和频率漂移偏差信息. 定义${t_k}$时刻状态向量的顺序依次为相位偏差、频率偏差和频率漂移偏差
$\begin{split}& X(t_k) = \\& [x_1(t_k)~y_1(t_k)~z_1(t_k), \cdots, x_N (t_k) y_N(t_k) z_N(t_k)], \end{split} $
式中, ${x_i}\left( t \right)$, ${y_i}\left( t \right)$和${z_i}\left( t \right)$分别表示钟$i$在$t$时刻的相位偏差、频率偏差和频率漂移偏差. 第$i$台钟从时刻${t_{k - 1}}$到时刻${t_k}$的状态由下面方程描述[6]:
$\begin{split}x_i(t_k)={}& x_i(t_{k -1}) + y_i(t_{k - 1}) \cdot \tau \\&+ \frac{1}{2} \cdot {\tau ^2} \cdot {z_i} (t_{k - 1}) + w_{xi}(t_k), \end{split} $
${y_i}\left( {{t_k}} \right) = {y_i}\left( {{t_{k - 1}}} \right) + {z_i}\left( {{t_{k - 1}}} \right) \cdot \tau + {w_{yi}}\left( {{t_k}} \right), $
${z_i}\left( {{t_k}} \right) = {z_i}\left( {{t_{k - 1}}} \right) + {w_{zi}}\left( {{t_k}} \right).$
(2)式和(3)式中$\tau = {t_k} - {t_{k - 1}}$, 表示相邻时刻的时间间隔; ${w_{xi}}\left( {{t_k}} \right)$, ${w_{yi}}\left( {{t_k}} \right)$和${w_{zi}}\left( {{t_k}} \right)$分别表示${t_k}$时刻的调频白噪声、调频随机游走噪声和频漂随机游走噪声; ${{{W}}_i}({t_k}) = {\left[ {{w_{xi}}({t_k}), {w_{yi}}({t_k}), {w_{zi}}({t_k})} \right]^{\rm T}}$表示原子钟状态估计的过程噪声向量, ${{{W}}_i} \left( {{t_k}} \right)\sim N\left( {0, {{{Q}}_i}} \right)$, ${{{Q}}_i}$表示第$i$台钟的过程噪声协方差矩阵,
$ \begin{split}& {{{Q}}_{i}} (\tau ) = \\& \left[ {\begin{array}{*{20}{c}} {{q_{xi}}\tau + {q_{yi}}\dfrac{{{\tau ^3}}}{3} + {q_{zi}}\dfrac{{{\tau ^5}}}{{20}}}&{{q_{yi}}\dfrac{{{\tau ^2}}}{2} + {q_{zi}}\dfrac{{{\tau ^4}}}{8}}&{{q_{zi}}\dfrac{{{\tau ^3}}}{6}} \\ {{q_{yi}}\dfrac{{{\tau ^2}}}{2} + {q_{zi}}\dfrac{{{\tau ^4}}}{8}}&{{q_{yi}}\tau + {q_{zi}}\dfrac{{{\tau ^3}}}{3}}&{{q_{zi}}\dfrac{{{\tau ^2}}}{2}} \\ {{q_{zi}}\dfrac{{{\tau ^3}}}{6}}&{{q_{zi}}\dfrac{{{\tau ^2}}}{2}}&{{q_{zi}}\tau }\end{array}} \right] .\end{split}$
(5)式中${q_{xi}}$, ${q_{yi}}$和${q_{zi}}$分别表示钟$i$的调频白噪声、调频随机游走噪声和频漂随机游走噪声的强度, 与Hadamard方差的关系为[14]
$H\sigma _{yi}^2\left( \tau \right) = \frac{{{q_{xi}}}}{\tau } + \frac{{{q_{yi}}\tau }}{6} + \frac{{11{q_{zi}}{\tau ^3}}}{{120}}.$
假设以钟1为参考钟, ${t_k}$时刻钟$i$的钟差测量表示为
${z_{i1}}\left( {{t_k}} \right) = {x_i}\left( {{t_k}} \right) - {x_1}\left( {{t_k}} \right) + {v_{1i}}\left( {{t_k}} \right).$
为了利用Kalman滤波算法, 需要建立原子钟组的矩阵向量方程, 其中状态向量方程表达为[15]
${{X}}\left( {{t_k}} \right) = {{\varPhi}} \left( \tau \right){{X}}\left( {{t_{k - 1}}} \right) + {{W}}\left( {{t_k}} \right).$
文献[15]说明GPS是一个由原子钟组成的网络, 早期的GPS设计将其中的一台原子钟视为“主钟”, 当“主钟”出现故障时会影响系统的性能, 因此这种设计过多依赖于主钟. 提出了一种称为GPS组合时钟的解决方案, 本质上它是一个Kalman滤波应用到原始的轨道/时钟问题, 包括所有时钟, 没有做主钟的选择, 提高了系统的可靠性. (8)式中${{\varPhi}} \left( \tau \right)$为状态转移矩阵, 有对角块形式
${{\varPhi}} \left( \tau \right) = \left[ {\begin{array}{*{20}{c}} 1&\tau &{{{{\tau ^2}}}/{2}} \\ 0&1&\tau \\ 0&0&1 \end{array}} \right].$
${t_k}$时刻过程噪声向量${{W}}({{t_k}}) \!=\! [ {{{W}}_1}( {{t_k}}), \cdots\!, {{{W}}_N} (t_k) ]^{\rm{T}}$,其中${\rm{E}}{{W}}\left( {{{{t}}_k}} \right) \!=\! 0$, ${\rm{E}}{{W}}({{t_k}}){{W}}{({{t_j}})^{\rm{T}}} \!=\! {{Q}}( \tau ) {\delta _{kj}}$, 协方差矩阵${{Q}}\left( \tau \right)$是由${{{Q}}_{i}}\left( \tau \right)$组成的对角块矩阵.
测量方程表达式为
${{Z}}\left( {{t_k}} \right) = {{H}}\left( {{t_k}} \right){{X}}\left( {{t_k}} \right) + {{V}}\left( {{t_k}} \right).$
(9)式中${{Z}}\left( {{t_k}} \right) = {\left[ {{z_{21}}\left( {{t_k}} \right), L, {z_{N1}}\left( {{t_k}} \right)} \right]^{\rm T}}$, ${{H}}\left( {{t_k}} \right)$是$\left( {N - 1} \right) \times 3 N$矩阵. 对于三台钟的情况:
${{H}}\left( {{t_k}} \right) = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&1&0&0&0&0&0 \\ { - 1}&0&0&0&0&0&1&0&0 \end{array}} \right].$
测量噪声向量${{V}}\left( {{t_k}} \right) = {\left[ {{v_{12}}\left( {{t_k}} \right), \cdots, {v_{1 N}}\left( {{t_k}} \right)} \right]^{\rm T}}$, 其中${\rm{E}}{{V}}\left( {{t_k}} \right) = 0$, ${\rm{E}}{{V}}\left( {{t_k}} \right){{V}}{\left( {{t_j}} \right)^{\rm{T}}} = {\varSigma _k}{\delta _{kj}}$, ${\varSigma _k}$是测量噪声方差矩阵, 并且${\rm{E}}{{W}}\left( {{t_k}} \right){{V}}{\left( {{t_j}} \right)^{\rm{T}}} = 0$.
3.改进Kalman滤波状态估计
改进Kalman滤波时间尺度算法与传统Kalman滤波时间尺度算法步骤基本相同, 分为两部分, 包括利用原子钟状态方程得到的状态预测估计和利用测量数据对状态预测估计的修正, 即状态验后估计, 对于时刻${t_k}$, 具体计算步骤如下.
2
3.1.算　法
-->

3.1.算　法

1)预测. 利用状态转移矩阵${{\varPhi }}\left( \tau \right)$和原子钟模型噪声${{W}}\left( {{t_k}} \right)$计算得到${t_k}$时刻原子钟状态的线性最小方差预测, 表示为
$\bar {{X}}\left( {{t_k}} \right) = {{\varPhi}} \left( \tau \right) \cdot \hat {{X}}\left( {{t_{k - 1}}} \right), $
式中, $\bar {{X}}\left( {{t_k}} \right)$表示${t_k}$时刻原子钟的预测状态估计, $\hat {{X}}\left( {{t_{k - 1}}} \right)$表示${t_{k - 1}}$时刻的原子钟的状态估计, 预测协方差矩阵表示为
${\varSigma _{{{{\bar X}}_k}}} = {{\varPhi}} \left( \tau \right) \cdot {\varSigma _{{{{\hat X}}_{k - 1}}}} \cdot {{{\varPhi}} ^{\rm T}}\left( \tau \right) +{{ Q}}\left( \tau \right), $
式中, ${\varSigma _{{{{{\bar X}}}_{{k}}}}}$表示${t_k}$时刻原子钟的预测状态协方差矩阵, ${\varSigma _{{{{{\hat X}}}_{k - 1}}}}$表示${t_{k - 1}}$时刻原子钟的状态估计协方差矩阵.
2)测量. 通过${t_k}$时刻原子钟测量数据${{Z}}\left( {{{{t}}_k}} \right)$估计新息向量, 新息向量表示为测量值与预测值之差, 表达式为
${{i}}\left( {{t_k}} \right) = {{Z}}\left( {{t_k}} \right) - {{H}}\left( {{t_k}} \right) \cdot \bar {{X}}\left( {{t_k}} \right).$
新息向量的协方差矩阵表达式为
${{C}}\left( {{t_k}} \right) = {{H}}\left( {{t_k}} \right) \cdot {\varSigma _{{{{{\bar X}}}_{{k}}}}} \cdot {{H}}{\left( {{t_k}} \right)^{\rm{T}}} + {\varSigma _k}, $
式中, ${\varSigma _k}$表示${t_k}$时刻测量值${{Z}}\left( {{t_k}} \right)$的误差协方差矩阵. 改进Kalman滤波时间尺度算法的增益矩阵表示为[16]
$\begin{split}{\bar {{K}}_k} =& \dfrac{1}{{{\alpha _k}}} \cdot {\varSigma _{{{\bar X}_k}}} \cdot {{{H}}^{\rm T}}\left( {{t_k}} \right) \! \Big(\frac{1}{{{\alpha _k}}} \cdot {{H}}({{t_k}}) \\& \times {\varSigma _{{{{{\bar X}}}_k}}} \cdot {{{H}}^{\rm T}}({{t_k}}) \!+\! {\varSigma _k} \!\Big)^{-1}, \end{split}$
式中, ${\alpha _k}$为自适应因子, 用于原子钟状态预测信息的自适应估计. 文献[16]表明运动定位与导航一般应用Kalman滤波算法, 可靠的Kalman滤波算法要求有可靠的函数模型、随机模型以及合理的估计方法. 然而精确的函数模型的构造十分困难, 随机模型先验信息的获取一般都是基于验前统计信息, 而任何统计信息都有可能失真. 本文围绕如何利用当前观测信息和状态估计更新先验信息, 提高状态估计的准确性展开研究. 在Kalman滤波中, 这种由估计过程自适应地调整、更新先验信息的算法称为自适应Kalman滤波.
根据以上信息得到${t_k}$时刻原子钟的状态估计, 表达式为
$\hat {{X}}\!\left( {{t_k}} \right) \!=\! \bar {{X}}\!\left( {{t_k}} \right)\! +\! {\bar K_k} \cdot\! \left[ {{{Z}}\!\left( {{t_k}} \right) \!-\! {{H}}\left( {{t_k}} \right) \cdot \bar {{X}}\!\left( {{t_k}} \right)} \right].$
${t_k}$时刻原子钟状态估计的协方差矩阵表达式为
${\varSigma _{{{\hat X}}\left( {{t_k}} \right)}} = \left[ {{{I}} - {{\bar {{K}}}_k} \cdot {{H}}\left( {{t_k}} \right)} \right] \cdot {\varSigma _{{{{{\bar X}}}_{{k}}}}}, $
式中, I表示单位矩阵, (14)式对预测状态协方差矩阵添加自适应因子${\alpha _k}$, 自适应因子用于控制模型误差对原子钟状态估计的影响, 通过构造模型误差判别统计量, 判断统计量与给定阈值的大小建立自适应因子函数, 当模型误差判别统计量小于给定阈值时, 自适应因子${\alpha _k}=1$, 改进的Kalman滤波算法变为传统Kalman滤波算法. 当模型误差判别统计量大于给定阈值时, 自适应因子增大, 模型预测状态协方差矩阵降低, 从而减小模型误差对原子钟状态估计的影响.
2
3.2.自适应因子的求解
-->

3.2.自适应因子的求解

改进Kalman滤波算法中自适应因子的确定, 对于原子钟状态估计, 由于测量信息不足以解算当前历元的状态参数, 故不能采用状态参数不符值构造统计量, 由于预测残差向量能够较好地反应模型变化, 考虑采用预测残差构造判别统计量, 利用三段函数构造自适应因子, 自适应因子表达式为[17-19]
${\alpha _k}\! = \!\left\{ {\begin{array}{*{20}{l}} {1,} \\ {\dfrac{{{c_0}}}{{| {{{\tilde V}_k}}|}} \cdot {{\left( {\dfrac{{{c_1} - | {{{\tilde V}_k}}|}}{{{c_1} - {c_0}}}} \right)}^2}\!,} \\ {0,} \end{array}} \right.\;\; \begin{array}{*{20}{l}} {| {{{\tilde V}_k}}| \leqslant {c_0},} \\ {{c_0} < | {{{\tilde V}_k}}| \leqslant {c_1},} \\ {| {{{\tilde V}_k}} | > {c_1},} \end{array}$
式中, ${c_0}$可取为1.0—1.5; ${c_1}$可取为3.0—8.5.
${\tilde {{V}}_k} = {{\left\| {{{\bar {{V}}}_k}} \right\|}\Big/ {\sqrt {{\rm{tr}}\left( {{\varSigma _{{{{{\bar V}}}_k}}}} \right)} }}, $
式中, ${\bar {{V}}_k}$表示${t_k}$时刻预测残差向量, ${\rm{tr}}\left( {{\varSigma _{{{{{\bar V}}}_k}}}} \right)$表示预测残差向量协方差矩阵的迹.
4.取权方法
原子钟权重算法的原则为一台好的钟是一台“可预测”的钟[20,21], 主要思想是原子钟的预测频率与实际频率的差异在每个月的计算区间内进行估计, 使用这些值来定义权重, 并适当设定权重的上限. 算法进行4次迭代, 每次迭代运行如下.
1) 对于值$\left[ {{{\bar h}_0} - {h_i}} \right]$是使用一组相对权重计算得到的, 式中${\bar h_0}$表示时间尺度, ${h_i}$表示第$i$台钟的读数. 第一次迭代的权重使用最近区间归一化后的权重值. 在接下来的迭代中, 权重根据上一次迭代计算得到.
2) 频率值${y_{i, {I_k}}}$与预测频率${\hat y_{i, {I_k}}}$之差的绝对值表示为
${\varepsilon _{i,{I_k}}} = \left| {{y_{i,{I_k}}} - {{\hat y}_{i,{I_k}}}} \right|, $
式中, $i$表示时钟的索引, ${I_k}$表示时间间隔, ${y_{i, {I_k}}}$表示频率实测值, ${\hat y_{i, {I_k}}}$表示频率预测值.
3) 考虑到新的测量数据具有更可靠的统计性质, 通过滤波器使得较新的测量数据具有更重要的作用, 计算
$\sigma _i^2 = \frac{{\displaystyle\sum\limits_{j = 1}^M {\left( {\frac{{M + 1 - j}}{M}} \right)\varepsilon _{i,j}^2} }}{{\displaystyle\sum\limits_{j = 1}^M {\left( {\frac{{M + 1 - j}}{M}} \right)} }}, $
其中$i$表示钟; $j$表示计算区间; $M$为计算区间的数目.
4) 计算时钟${H_i}$的权重为
${\omega _{i,{\rm{temp}}}} = \frac{{{1/ {\sigma _i^2}}}}{{\displaystyle\sum\nolimits_{i = 1}^N {{1/ {\sigma _i^2}}} }}.$
除了时钟${H_i}$的权重大于设置的最大权和时钟${H_i}$在计算区间出现异常外, 时钟${H_i}$的相对权重${\omega _i}$$={\omega _{i, {\rm{temp}}}}$.
5.时间尺度计算与分析
为了检验噪声强度变化时改进Kalman滤波时间尺度算法的性能, 首先模拟原子钟的数据, 通过改变原子钟运行过程中的噪声强度, 模拟原子钟噪声强度变化的情况, 研究改进Kalman滤波时间尺度算法的性能, 最后根据实测数据验证噪声强度估计存在误差时算法的性能并分析.
2
5.1.模拟数据
-->

5.1.模拟数据

根据原子钟的数学模型模拟原子钟的数据[22,23], 原子钟系统随机微分方程的解能够表达为迭代形式:
${{X}}\left( {{t_{k + 1}}} \right) = {{\varPhi}} \left( \tau \right){{X}}\left( {{t_k}} \right) + {{W}}\left( {{t_{k + 1}}} \right).$
对于单台原子钟, (22)式中${{X}}(t) = \left[ {x(t), y(t), z(t)} \right]$, $x(t)$表示模拟原子钟相对于参考原子钟的相位偏差, $y(t)$表示频率偏差的一部分, $z(t)$表示频率漂移偏差的一部分. ${{\varPhi}} \left( \tau \right)$表达式与(8)式相同, ${{W}}\left( {{t_k}} \right)$是原子钟的噪声向量; ${{Q}}\left( \tau \right)$为协方差矩阵, 与(8)式中的定义相同, ${q_{xi}}$, ${q_{yi}}$和${q_{zi}}$与(5)式中表达的意义相同.
产生钟的数据后, 计算钟差测量值, 钟差测量值为钟组中原子钟与某台参考钟的差值, 文中以基准钟为参考, 计算各台原子钟与参考钟的钟差, 通过测量协方差矩阵${\Sigma _k}$产生测量噪声加入到钟差, 计算方法如(9)式. 模拟中假设原子钟测量方差均为$1 \times {10^{ - 20}}$ s2, 且测量值互不相关. 根据原子钟系统模型, 钟的噪声强度和测量噪声对原子钟进行模拟, 文中模拟5台原子钟, 包括2台氢原子钟, 2台铯原子钟和1台基准钟[24].
为了能够清楚的反映原子钟噪声强度变化的影响, 假设原子钟比对测量间隔为300 s, 即每间隔300 s测量一次, 测量周期为3000次. 当原子钟运行至1001次时改变原子钟的噪声强度, 取两台模拟氢原子钟为例, 改变模拟氢原子钟编号为AHM1和编号为AHM2的噪声强度, AHM1和AHM2的噪声强度分别为: ${q_{x1}}$= 5.8×10–26 (s2/s), ${q_{y2}}$= 5.1×10–34 (s2/s3), ${q_{z1}}$= 9.4×10–40 (s2/s5)和${q_{x2}}$= 5.9×10–26 (s2/s), ${q_{y2}}$= 6.2×10–34 (s2/s3), ${q_{z2}}$= 9.5×10–51 (s2/s5). 当原子钟运行至2001次时再次改变模拟氢原子钟的噪声强度, AHM1和AHM2的噪声强度分别为: ${q_{x1}}$= 5.8×10–26 (s2/s), ${q_{y1}}$= 9.1×10–34 (s2/s3), ${q_{z1}}$= 9.4×10–51 (s2/s5)和${q_{x2}}$= 5.9×10–26 (s2/s), ${q_{y2}}$= 9.2×10–34 (s2/s3), ${q_{z2}}$= 9.5×10–40 (s2/s5). 分别基于经典Kalman滤波时间尺度算法和改进Kalman滤波时间尺度算法估计原子钟的状态, 并计算各台原子钟的时间偏差改正值(corrected time deviations), 时间偏差改正值表达式为
${\tilde x_i}({t_k}) = {x_i}\left( {{t_k}} \right) - {\hat x_i}\left( {{t_k}} \right), $
式中, ${\tilde x_i}({t_k})$表示${t_k}$时刻原子钟$i$的时间偏差改正值, ${\hat x_i}\left( {{t_k}} \right)$表示${t_k}$时刻原子钟$i$的时间偏差估计值. AHM1和AHM2的时间偏差改正值分别如图1和图2所示. 计算表明, 基于经典Kalman滤波时间尺度算法, AHM1与AHM2的时间偏差改正值分别约有4 ns和3 ns的扰动, 而基于改进Kalman滤波时间尺度算法得到的时间偏差改正值不存在扰动现象, 说明基于经典Kalman滤波时间尺度算法得到的原子钟状态估计受原子钟噪声强度变化的影响, 而改进Kalman滤波时间尺度算法有效抵制了原子钟噪声强度变化对状态估计的影响, 提高了状态估计的准确度. 基于原子钟时间偏差改正值计算时间尺度, 计算公式为
图 1 基于两种算法模拟的AHM1时间偏差改正值　(a) 基于Kalman滤波算法模拟的AHM1的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM1时间偏差改正值
Figure1. The corrected time deviations of modeling AHM1 based on two algorithms:　(a) The corrected time deviations of modelling AHM1 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM1 based on modified Kalman filter algorithm.

图 2 基于两种算法的模拟AHM2的时间偏差改正值　(a) 基于Kalman滤波算法模拟的AHM2的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM2的时间偏差改正值
Figure2. The corrected time deviations of modelling AHM2 based on two algorithms:　(a) The corrected time deviations of modelling AHM2 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM2 based on modified Kalman filter algorithm.

${\bar x_0}\left( {{t_k}} \right) = \sum\limits_{i = 1}^N {{w_i}} {\tilde x_i}\left( {{t_k}} \right).$
不妨假设${w_i}$取等权, 两种Kalman滤波算法的时间尺度如图3所示, 图3中经典Kalman滤波时间尺度用KF Scale表示, 改进Kalman滤波时间尺度用R-KF Scale表示, 图4和表1中的表列与图3相同. 基于经典Kalman滤波算法的时间尺度受噪声强度变化的影响, 出现明显的扰动, 降低了时间尺度的准确度. 基于改进Kalman滤波算法的时间尺度有效抵制了噪声强度变化对时间尺度的影响, 提高了时间尺度的准确度.
取样时间/102 sKF scale/10–14R-KF scale/10–14
39.8513.80
66.738.75
124.385.20
302.502.70
601.761.85
1201.221.24
3000.850.73
6000.880.57
12000.890.63


表1基于两种Kalman滤波算法的时间尺度的重叠Allan偏差
Table1.The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

图 3 基于两种Kalman滤波算法的时间尺度
Figure3. Time scale based on two Kalman filter algorithms.

图 4 基于两种Kalman滤波算法的时间尺度的重叠Allan偏差
Figure4. The overlapping Allan deviation of time scale based on two Kalman filter algorithms.

基于重叠阿伦偏差(overlapping Allan deviation)分析频率稳定度, 两种Kalman滤波算法的时间尺度的重叠阿伦偏差如图4与表1所列, 分析得出相比于经典Kalman滤波算法, 基于改进Kalman滤波算法得到的时间尺度长期稳定度提高, 短期稳定度有所降低, 分析原因, 对Kalman滤波预测协方差矩阵引入自适应因子, 控制了噪声强度变化对时间尺度的影响, 降低了扰动, 提高了长期稳定度, 但同时控制作用影响了相位偏差改正值的短期稳定度, 进而影响了时间尺度的短期稳定度.
2
5.2.实测数据
-->

5.2.实测数据

不同于模拟钟情况, 对于原子钟实际测量数据定义的Kalman滤波时间尺度是非观测的. 比如, 不能得到单台原子钟的读数, 但是可以得到原子钟之间的钟差数据, 原子钟相对于参考钟或时间尺度的钟差数据. 因此, 只能通过计算得到Kalman滤波时间尺度与钟组中单台原子钟的钟差值. 原子钟实测数据的Kalman滤波时间尺度的计算式为
$\begin{split} {\bar h_0}\left( t \right) - {h_j}\left( t \right)& =\sum\limits_{i = 1}^N {{w_i}} \left( {{{\tilde h}_i}\left( t \right) - {h_j}\left( t \right)} \right)\\ &= \sum\limits_{i = 1}^N {{w_i}\left[ {{x_{ij}}(t) + {{\hat x}_i}(t)} \right]}, \end{split}$
式中, ${\bar h_0}\left( t \right)$表示$t$时刻Kalman滤波时间尺度; ${h_j}\left( t \right)$表示$t$时刻原子钟$j$的读数; ${\tilde h_i}\left( t \right)$表示改正钟, 有关系式${\tilde h_i}\left( t \right) = {h_0}(t) - {\tilde x_i}\left( t \right)$, ${h_0}(t)$表示$t$时刻理想钟的读数; ${w_i}$表示第$i$台原子钟的权重; ${x_{ij}}(t)$表示$t$时刻钟$i$与钟$j$的钟差值.
为了检验噪声强度估计存在误差时改进Kalman滤波时间尺度算法的性能, 随机选取中国科学院国家授时中心(national time service center, NTSC)时频基准实验室的5台原子钟, 在国际原子时(international atomic time, TAI)计算中的编号分别为HM4926, HM4967, HM0296, Cs3436和Cs2098. 钟差UTC(NTSC)-Clock(i)为用于时间尺度计算的数据形式, 数据长度为MJD58484 —MJD58726.96 (2019年 1月1日0时—8月31日23时), 数据采样间隔为1 h, 时间尺度根据第4节“可预测性”方法计算相对权重, 并且最大权设为0.3, 计算结果如表2所列, 得到HM4926与HM4967取得最大权重为0.3, 说明这2台氢原子钟的可预测性能优于其余3台原子钟, Cs2098取得最小权重0.07, 表明相对于钟组中其余原子钟, 其可预测性能低.
原子钟编号HM4926HM4967HM0296Cs3436Cs2098
相对权重0.300.300.140.190.07


表2原子钟的相对权重
Table2.Relative weights of atomic clocks.

研究原子钟噪声强度估计存在误差时改进Kalman滤波时间尺度算法的性能. 随机改变原子钟的噪声强度, 如表3所列. 改变氢原子钟的噪声强度, 2台铯原子钟的调频随机游走噪声与频漂随机游走噪声强度为0. 分别基于经典Kalman滤波时间尺度算法和改进Kalman滤波算法计算时间尺度, 计算结果与UTC (NTSC) 进行比较, 如图5所示. 基于经典Kalman滤波算法的时间尺度明显偏离UTC (NTSC), 而基于改进Kalman滤波算法的时间尺度保持与UTC (NTSC)一致. 表4给出了通过两种Kalman滤波算法得到的时间尺度与UTC (NTSC)的最大偏差值, 最小偏差值与平均偏差值. 相比于经典Kalman滤波算法, 基于改进Kalman滤波算法得到的时间尺度与UTC (NTSC) 相比, 最大偏差值降低了3.25 ns, 最小偏差值改进了0.48 ns, 平均绝对偏差值改进了0.51 ns, 提高了时间尺度的准确度.
原子钟
编号		q1WFM/10–36 [s2/s]q2RWFM/10–50 [s2/s3]q3RWFM/[s2/s5]
HM49264.644.994.27 × 10–60
HM49675.884.603.87 × 10–60
HM02962.014.106.80 × 10–49
Cs343626.9000
Cs20985.4300


表3原子钟的噪声强度
Table3.Noise intensity of atomic clocks.

计算方法最大偏
差值/ns		最小值偏
差值/ns		平均偏
差值/ns
KF scale4.36–3.130.98
R-KF scale1.11–2.65–0.47


表4基于两种Kalman滤波算法的统计值
Table4.Statistical values based on two Kalman filter algorithms.

图 5 基于两种Kalman滤波算法的时间尺度
Figure5. Time scales based on two Kalman filter algorithms.

基于重叠阿伦偏差估计时间尺度的频率稳定度, 分别计算基于两种Kalman滤波算法得到的时间尺度的重叠阿伦偏差. 因为系统中原子钟的采样时间为1 h, 因此重叠阿伦偏差的最小平滑时间也为1 h, 计算结果如图6与表5所列. 相比于经典Kalman滤波算法, 基于改进Kalman滤波算法的时间尺度在相应平滑时间的重叠阿伦偏差有所降低, 说明稳定度有所提高. 分析其原因是受原子钟噪声强度估计误差的影响, 基于经典Kalman滤波的时间尺度的频率稳定度降低.
图 6 基于两种Kalman滤波算法的时间尺度的稳定度
Figure6. Time scale stabilities based on two Kalman filter algorithms.

Sample time/103 sKF scale/10–14R-KF scale/10–14
3.601.941.64
7.201.341.17
14.401.030.92
36.000.920.85
72.000.740.69
144.000.530.51
360.000.270.25


表5基于两种Kalman滤波算法的时间尺度的重叠阿伦偏差值
Table5.The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

原子钟噪声强度估计存在误差时, 基于经典Kalman滤波算法的时间尺度准确度降低, 噪声强度估计误差影响了时间尺度的准确度. 改进Kalman滤波算法对状态预测协方差矩阵引入自适应因子, 控制噪声强度估计误差引起的状态预测协方差矩阵的增长, 降低原子钟状态估计的扰动, 提高了时间尺度的准确度.
6.结　论
针对原子钟运行过程中的噪声强度变化和噪声强度估计存在误差的情况, 本文研究了一种改进的Kalman滤波时间尺度算法, 通过原子钟模拟数据和实测数据验证了算法的性能并进行了分析. 该算法使用预测残差构造判别统计量, 并构造自适应因子, 通过对经典Kalman滤波预测协方差矩阵引入自适应因子, 实时控制原子钟噪声强度变化与噪声强度估计误差引起的状态估计扰动, 提高了时间尺度的准确度, 保证了时间尺度的长期稳定度. 同时结合模拟数据分析, 实时控制噪声强度变化引起的原子钟状态估计扰动过程也影响了时间尺度的短期稳定度. 进一步研究自适应因子的构造方法, 降低由于控制作用影响时间尺度的短期稳定性能是下一步的重点工作.
模拟数据和实测数据的结果都计算了重叠Allan偏差, 如图4和图6所示. 图4和图6中的结果特征有一些差别, 分析原因, 与原子钟的取权算法有关, 文中实测数据采用“可预测性”取权算法是为了保证时间尺度的长期稳定度.
根据时间尺度的计算过程可以得出, 基于经典Kalman滤波时间尺度的性能也与发生噪声强度变化原子钟的权重有关, 噪声强度变化的原子钟权重越大, 对时间尺度性能的影响越大.
相关话题/计算 数据 测量 文献 信息
  • 领限时大额优惠券,享本站正版考研考试资料!
    大额优惠券
    优惠券领取后72小时内有效，10万种最新考研考试考证类电子打印资料任你选。涵盖全国500余所院校考研专业课、200多种职业资格考试、1100多种经典教材，产品类型包含电子书、题库、全套资料以及视频，无论您是考研复习、考证刷题，还是考前冲刺等，不同类型的产品可满足您学习上的不同需求。 ...
    考试优惠券 本站小编 Free壹佰分学习网 2022-09-19
  • 光场高阶光子关联的分析与测量
    摘要:利用双HanburyBrown-Twiss探测系统理论分析并在实验上精确测量了不同光场的高阶光子关联g(n)(n>2).探测系统通过四个单光子计数模块,探测分析光子时间关联的联合分布概率.在理论上,考虑实际探测系统背景噪声和系统效率的影响,分析研究了热态、相干态、压缩真空态和Fock态的三阶及 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 质子成像法测量电容线圈靶磁场
    摘要:质子背光成像技术是一种诊断等离子体电磁场的重要诊断手段.当质子穿过等离子体的电磁场,质子受洛伦兹力影响在成像板上重新分布.如何从质子成像结果中重构电磁场是一个非常重要的研究课题.本文以激光驱动电容线圈靶为例介绍和对比了粒子追踪法和流量分析法这两种通过质子成像结果重构磁场的方法.激光驱动电容线圈 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 荧光寿命数据的相量分析及其应用
    摘要:荧光寿命显微成像技术(fluorescencelifetimeimagingmicroscopy,FLIM)具有特异性强、灵敏度高、可定量测量等优点,被广泛应用于生物医学、材料学等领域的研究.为使FLIM技术更好地适用于高通量数据的快速分析,近年来涌现出多种荧光寿命分析的新算法.其中,相量分析 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 空域模拟光学计算器件的研究进展
    摘要:空域模拟光学计算器件具备高通量、实时性和低能耗的信息处理能力.光学超构材料结构紧凑、对光波具有强大调控能力,可被用于构建小型化、集成化的空域模拟光学计算器件.目前空域模拟光学计算器件的研究根据设计方法主要分为4F系统法和格林函数法两类.4F系统法需要两个傅里叶变换透镜和一个空间频率滤波器,实际 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 信息超材料研究进展
    摘要:超材料是物理和信息领域的研究热点之一,本文主要介绍信息超材料的研究进展.不同于传统超材料的等效媒质参数表征,信息超材料由物理单元的数字编码来描述,通过控制不同的编码序列来实时地调控电磁波,进而实现超材料的现场可编程功能.由于在超材料的物理空间上构筑起数字空间,因此可在超材料的物理平台上直接处理 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于高分辨率激光外差光谱反演大气CO<sub>2</sub>柱浓度及系统测量误差评估方法
    摘要:利用实验室研制的近红外激光外差光谱仪,开展了基于最优估计算法的温室气体柱浓度反演和系统测量误差的近似评估等相关工作.首先,通过光谱数据库、参考正向模型计算结果与傅里叶变换红外光谱技术探测结果筛选出了探测窗口,并以此为依据选择了相应的激光器和探测器;其次,建立了基于参考正向模型最优估计浓度反演算 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 绝热跃迁方法测量铯喷泉钟冷原子碰撞频移
    摘要:冷原子碰撞频移是限制铯原子喷泉钟频率不确定度性能的主要因素之一.在使用外推法测量冷原子碰撞频移时,制备密度均匀成比例的原子团是减小系统误差的关键.绝热跃迁方法可以用来实现均匀跃迁比例,均匀度可达10–3.通过理论分析Bloch矢量的演化,导出了误差满足的方程,实验测量了不同参数对跃迁几率的影响 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于增强型视觉密码的光学信息隐藏系统
    摘要:提出了一种基于增强型视觉密码的光学信息隐藏系统.该系统可将秘密图像分解为多幅有实际意义的分享图像,然后将这些分享图像隐藏在相位密钥中,相位密钥可以制成衍射光学元件,以实体的形式保存和传输,扩展了视觉密码的应用范围.在提取过程中,只需要使用激光照射衍射光学元件,再现分享图像,然后只需要将一定数量 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 脉冲硬X射线能注量测量技术
    摘要:阐述了全吸收法测量脉冲硬X射线能注量的基本原理,选择了光电管配合硅酸镥(LSO)闪烁体作为探测系统的核心部件,研制了脉冲硬X射线能注量测量系统.利用该系统测量了“闪光二号”加速器产生的脉冲硬X射线强度,结合灵敏度的实验标定结果,计算得到了脉冲硬X射线的能注量.在连续5发次的实验中,能注量平均测 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 蛋白质“液-液相分离”的理论和计算方法进展
    摘要:蛋白质分子的“液-液相分离”是近几年生物物理学领域迅速发展起来的研究热点.蛋白质相分离在一系列生物学过程中发挥着重要的作用.蛋白质分子序列和构象的多样性和复杂性,给蛋白质分子的理论研究、计算机模拟和实验研究都带来了巨大的挑战.当前,多尺度理论模型和多分辨率计算方法被广泛地用于蛋白质分子的“液- ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29

Free考研考试FreeKaoYan.Com
欢迎来到Free考研考试，"为实现人生的Free而奋斗"

© 2020 FreeKaoYan! . All rights reserved.