CN110350996B - 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统 - Google Patents

基于交互式多模型滤波器的时钟漂移率跟踪方法及系统 Download PDF

Info

Publication number
CN110350996B
CN110350996B CN201910574410.XA CN201910574410A CN110350996B CN 110350996 B CN110350996 B CN 110350996B CN 201910574410 A CN201910574410 A CN 201910574410A CN 110350996 B CN110350996 B CN 110350996B
Authority
CN
China
Prior art keywords
time
clock
kalman filter
drift rate
model
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
Application number
CN201910574410.XA
Other languages
English (en)
Other versions
CN110350996A (zh
Inventor
朱敏
刘迪
武岩波
房小芳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Acoustics CAS
Original Assignee
Institute of Acoustics CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Acoustics CAS filed Critical Institute of Acoustics CAS
Priority to CN201910574410.XA priority Critical patent/CN110350996B/zh
Publication of CN110350996A publication Critical patent/CN110350996A/zh
Application granted granted Critical
Publication of CN110350996B publication Critical patent/CN110350996B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04J—MULTIPLEX COMMUNICATION
    • H04J3/00—Time-division multiplex systems
    • H04J3/02—Details
    • H04J3/06—Synchronising arrangements
    • H04J3/0635—Clock or time synchronisation in a network
    • H04J3/0682—Clock or time synchronisation in a network by delay compensation, e.g. by compensation of propagation delay or variations thereof, by ranging

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Synchronisation In Digital Transmission Systems (AREA)

Abstract

本发明提出了一种基于交互式多模型滤波器的时钟漂移率跟踪方法及系统，所述方法包括：根据时钟偏差和漂移特性，计算初始状态参数；将初始状态参数进行输入交互，得到混合初始参数；将混合初始参数输入多模型Kalman滤波器进行滤波；输出变时钟漂移率的估计结果；根据变时钟漂移率的估计结果修正系统的时钟偏差。本发明的方法使用交互式多模型Kalman滤波器对变时钟状态向量进行跟踪和估计，并在滤波过程中采用Sage‑Husa自适应方法动态调整滤波参数，提高了水声传感器网络授时精度。基于上述方法，本发明还提出提出了一种基于交互式多模型滤波器的时钟漂移率跟踪的系统，一种计算机设备和一种计算机可读存储介质。

Description

基于交互式多模型滤波器的时钟漂移率跟踪方法及系统
技术领域
本发明涉及一种水声传感器网络的时钟漂移率跟踪方法，特别涉及一种基于交互式多模型滤波器的时钟漂移率跟踪方法及系统。
背景技术
基于水声通信技术的无线传感器网络作为海底观测网络的重要组成部分，在维护国家海洋权益、开发和保护海洋资源方面有重要意义，在民用和军用海洋信息开发领域具有巨大的发展潜力，吸引了大量学者进行研究。时钟同步是水声无线网各个节点协同观测的关键技术，对某些声学或地震信号传感器，在长期工作过程中需要保持较高的同步精度才能检测和分析地震、海啸等突发事件。另外，为节省传感器节点能量，延长水下无线传感网络的生命周期，大多数节点工作采用睡眠-唤醒机制，因此水下授时显得更为重要。
传感器节点使用晶振产生本地时钟。在实际使用中，出于对降低成本的考虑，普遍使用频率准确度低和稳定性差的廉价晶体振荡器。因制造工艺、温度等因素的影响，不同传感器节点晶振的输出频率会不同。节点时钟输出频率的差异称为时钟漂移率。若两时钟之间时钟漂移率为α意味着记录相同一段时间T，第一个时钟相对于第二个时钟记录的偏差为αT。时钟漂移率的存在会导致本地时钟之间产生累积误差，称为时钟偏差。这就要求授时协议要同时获得时钟漂移率和时钟偏差并进行补偿。
由于水声通信带宽窄，通信速率低，传播时延大，现有的许多针对陆地无线传感网络的授时算法，如ERBS(Energy-efficient Reference Broadcast Synchronization)、WPTP(Wireless Precision Time Protocol)、TPSN(Timing-sync Protocol for SensorNetworks)等算法，不能直接应用于水下。针对水声环境特性提出的多种授时协议，如TSHL(Time Synchronization for High Latency Acoustic Networks)授时协议、Tri-message授时协议等，考虑了时钟漂移率高时延、信道多径等因素对授时精度的影响，但没有相关研究考虑长时时钟漂移率变化对授时误差影响，导致水声传感器网络授时精度有很大的时钟偏差。
发明内容
本发明的目的在于解决现有的许多针对陆地无线传感网络的授时算法不能直接应用于水下以及长时时钟漂移率变化影响水声传感器网络授时精度的问题。
为实现上述目的，本发明提出一种基于交互式多模型滤波器的时钟漂移率跟踪方法，包括：
根据时钟偏差和漂移特性，计算初始状态参数；
将初始状态参数与当前时刻的转移概率进行交互运算，得到混合初始参数；
将混合初始参数输入多模型Kalman滤波器进行估计计算；输出变时钟漂移率的估计结果；
根据变时钟漂移率的估计结果修正系统的时钟偏差。
作为所述方法的一种改进，所述根据时钟偏差和漂移特性，计算初始状态参数，具体包括：
步骤1-1)设在k时刻的时钟离散状态向量为c(k)：
c(k)＝[θ(k)，α(k)，ρ(k)] (3)
其中，θ(k)表示k时刻的时钟偏差，α(k)表示k时刻的时钟漂移率，ρ(k)表示k 时刻的时钟漂移率变化率；
步骤1-2)k时刻时钟漂移率α(k)为：
α(k)＝α(k-1)+τ(k)×ρ(k)+ωα(k) (4)
其中，α(k-1)为k-1时刻的时钟漂移率，τ(k)表示k时刻采样间隔，ωα(k)表示 k时刻时钟漂移率为α(k)时的噪声向量；
步骤1-3)将时钟离散状态向量c(k)转化为k时刻的矩阵向量表达式：
c(k)＝Ac(k-1)+ω(k-1) (5)
其中，A为状态转移矩阵，ω(k-1)为k-1时刻的噪声向量；
步骤1-4)k时刻的观测向量z(k)为：
z(k)＝Hc(k)+v(k) (6)
其中H为观测矩阵，v(k)为观测噪声向量；
步骤1-5)对于输入的n个Kalman滤波器，对应n个估计模型mn，在k时刻的第j个Kalman滤波器对应的估计模型为mj(k)，j＝1，2，...，n；mj(k)的下一个时刻时钟离散状态向量cj(k+1)和当前时刻的观测向量zj(k)分别为：
cj(k+1)＝Ajcj(k+1)+ωj(k) (7)
zj(k)＝Hjcj(k)+vj(k) (8)
n为自然数，ωj(k)和vj(k)分别表示k时刻第j个估计模型mj(k)的噪声向量和观测噪声向量，其协方差矩阵分别为Q和R，将下一时刻时钟离散状态向量cj(k+1) 和当前的观测向量zj(k)作为初始状态参数。
作为所述方法的一种改进，将初始状态参数与当前时刻的转移概率进行交互运算，得到混合初始参数，具体包括：
步骤2-1)以μj(k)表示第j个Kalman滤波器对应的估计模型mj(k)在k时刻的转移概率，Pij表示第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)的状态转移矩阵，并设Pij为马尔科夫矩阵；则k-1时刻的模型预测概率
Figure BDA0002111708970000031
为：
Figure BDA0002111708970000032
其中，i＝1，2，...，n；第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)的转移概率μi/j(k-1|k-1)为：
Figure BDA0002111708970000033
步骤2-2)将初始状态参数在第i个Kalman滤波器对应的估计模型mi(k)到第j 个Kalman滤波器对应的估计模型mj(k)之间进行交互计算后得到初始混合状态参数包括：
Figure BDA0002111708970000034
Figure BDA0002111708970000041
其中，
Figure BDA0002111708970000042
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)的时钟离散状态向量估计，
Figure BDA0002111708970000043
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)输出误差的协方差矩阵估计，
Figure BDA0002111708970000044
为k-1时刻交互计算后的输入到滤波器j的时钟状态向量估计，
Figure BDA0002111708970000045
为k-1时刻交互计算后的输入到滤波器j的初始状态协方差矩阵估计。
作为所述方法的一种改进，所述将混合初始参数输入多模型Kalman滤波器进行滤波；输出变时钟漂移率的估计结果，具体包括：
步骤3-1)将混合初始参数
Figure BDA0002111708970000046
及z(k)输入n个Kalman滤波器，根据协方差矩阵Q和R为滤波参数，分别实时更新滤波参数向量
Figure BDA0002111708970000047
和
Figure BDA0002111708970000048
Figure BDA0002111708970000049
Figure BDA00021117089700000410
Figure BDA00021117089700000411
其中，dn为中间系数，b为遗忘因子，ε(k)为新息矩阵，上标T表示转置，K0j(k)为 k-1时刻的第j个Kalman滤波器的滤波增益矩阵，大小为n*n，下标j表示第j个 Kalman滤波器；
步骤3-2)对n个Kalman滤波器，首先分别计算k时刻第j个Kalman滤波器对应的估计模型mj(k)的更新概率μj(k)：
Figure BDA00021117089700000412
其中，Sj(k)为基于似然函数计算的k时刻第j个Kalman滤波器对应的估计模型mj(k)的动态权重向量，r为：
Figure BDA0002111708970000051
步骤3-3)输出变时钟漂移率的状态估计结果
Figure BDA0002111708970000052
Figure BDA0002111708970000053
本发明还提出一种基于交互式多模型Kalman滤波器的时钟漂移率跟踪系统，所述系统包括：输入的初始状态参数模块、输入交互模块和模型概率更新模块和输出组合模块；
所述输入的初始状态参数模块，用于根据k时刻发送与接收的时钟偏差和漂移特性，计算初始状态参数；
所述输入交互模块，用于将初始状态参数进行输入交互，得到混合初始参数；
所述多模型Kalman滤波器滤波模块，用于将混合初始参数输入多模型Kalman 滤波器进行滤波；输出变时钟漂移率的估计结果；
所述修正偏差模块，用于根据变时钟漂移率的估计结果修正系统的时钟偏差。
本发明还提出一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序所述处理器执行所述计算机程序时实现上述任一项所述的方法。
本发明还提出一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行上述任一项所述的方法。
本发明的优势在于：
本发明的方法使用交互式多模型Kalman滤波器对变时钟状态向量进行跟踪和估计，并在滤波过程中采用Sage-Husa自适应方法动态调整滤波参数，提高了水声传感器网络授时精度。
附图说明
图1为现有技术的混合使用TSHL和Tri-message授时协议计算时钟漂移率示意图；
图2为本发明的基于交互式多模型滤波器的时钟漂移率跟踪方法的交互式多模型Kalman滤波结构图；
图3为本发明的基于交互式多模型滤波器的时钟漂移率跟踪方法、混合使用 TSHL和Tri-message授时协议计算时钟漂移率方法、使用不调整滤波参数的Kalman 滤波时钟漂移率计算方法三种不同授时方法估计结果对比；
图4为本发明的基于交互式多模型滤波器的时钟漂移率跟踪方法、混合使用 TSHL和Tri-message授时协议计算时钟漂移率方法、使用不调整滤波参数的Kalman 滤波时钟漂移率计算方法三种不同授时方法累积授时误差对比；
图5为本发明的基于交互式多模型滤波器的时钟漂移率跟踪方法、混合使用 TSHL和Tri-message授时协议计算时钟漂移率方法、使用不调整滤波参数的Kalman 滤波时钟漂移率计算方法三种不同授时方法MSE对比。
具体实施方式
本发明的一种基于交互式多模型滤波器的时钟漂移率跟踪方法能够提高水声传感器网络授时精度。该方法使用现有水声授时协议TSHL和Tri-message计算获得时钟漂移率，作为交互式多模型Kalman滤波器的输入；使用交互式多模型Kalman滤波器对变时钟状态向量进行跟踪和估计，并在滤波过程中采用Sage-Husa自适应方法动态调整滤波参数，得到变时钟漂移率的估计结果，用估计结果对时钟漂移率进行修正和补偿，修正时钟偏差。
如图1所示，在现有技术中混合采用TSHL协议计算时钟漂移率和Tri-message 授时协议计算时钟漂移率；
所述采用TSHL计算时钟漂移率的方法为：
标准节点A向待同步节点B连续发送M个标准时间戳
Figure BDA0002111708970000061
其中，l表示第l次TSHL授时过程，l为自然数；
B节点根据接收到的M个消息时的本地时间戳
Figure BDA0002111708970000062
和相应的标准时间戳
Figure BDA0002111708970000063
在以下数据点进行线性拟合：
Figure BDA0002111708970000064
从而得到本地时钟B的时钟漂移率。
所述采用Tri-message计算时钟漂移率的方法为：
每一次Tri-message授时过程采用三次消息交换，B节点获得6个时间戳
Figure BDA0002111708970000071
后，用下式估计本地时钟漂移率αTri：
Figure BDA0002111708970000072
其中，s表示第s次Tri-message授时过程，s为自然数。
混合使用TSHL和Tri-message授时协议的方法为：
假设在混合使用TSHL和Tri-message授时协议的一次授时过程中，使用TSHL 方式估计时钟漂移率的时间间隔为4小时，标准节点A向待同步节点B连续发送 M＝25个标准时间戳，使用Tri-message方式估计时钟漂移率的时间间隔为0.5小时， TSHL与Tri-message两种方式交替使用，具体通信过程如图1所示。
本发明的交互式多模型Kalman滤波算法中的多个滤波器对应于不同的模型，各模型描述不同的时钟漂移率的变化特性，计算输入的初始状态参数：
设在k时刻时钟离散状态向量为c(k)：
c(k)＝[θ(k)，α(k)，ρ(k)]T (3)
其中，θ(k)表示k时刻时钟偏差，α(k)表示k时刻时钟漂移率，ρ(k)表示k时刻时钟漂移率变化率；使用恒定加速模型表征时钟漂移率的变化情况为：
α(k)＝α(k-1)+τ(k)×ρ(k)+ωα(k) (4)
式中τ(k)表示采样间隔，ωα(k)表示噪声量；
时钟离散状态向量表达式为：
c(k)＝Ac(k-1)+ω(k-1) (5)
其中，式中A为状态转移矩阵，ω(k-1)为k-1时刻的系统噪声向量。
k时刻观测向量z(k)为：
z(k)＝Hc(k)+v(k) (6)
式中H为观测矩阵，v(k)为观测噪声向量。
对具有n个模型(mj|j＝1，2，...，n)的系统，在k时刻的第j个模型mj(k)的下一个时钟离散状态向量cj(k+1)和当前的观测向量zj(k)为：
cj(k+1)＝Ajcj(k+1)+ωj(k) (7)
zj(k)＝Hjcj(k)+vj(k) (8)
n为自然数，ωj(k)和vj(k)分别表示k时刻第j个模型mj(k)的系统噪声向量和观测噪声向量，其协方差矩阵分别为Q和R。
如图2所示，使用交互式多模型估计时钟漂移率的过程，主要分成三步：输入交互，Kalman滤波和输出组合。
第一步：进行输入交互，计算k时刻匹配模型mj(k)的混合初始条件，μj(k)表示第j模型mj(k)在k时刻的概率，Pij表示第i模型mi(k)到第j模型mj(k)的状态转移矩阵，并设Pij为马尔科夫矩阵。k时刻的mj(k)模型预测概率
Figure BDA0002111708970000081
和模型转移概率μi/j(k-1|k-1)分别表示为：
Figure BDA0002111708970000082
Figure BDA0002111708970000083
输入交互后的混合输入表示为：
Figure BDA0002111708970000084
Figure BDA0002111708970000085
其中，
Figure BDA0002111708970000086
为k-1时刻输入交互后的输入到滤波器j的时钟状态向量估计，
Figure BDA0002111708970000087
为k-1时刻输入交互后的输入到滤波器j的初始状态协方差矩阵。
第二步：Kalman滤波：
将
Figure BDA0002111708970000088
及z(k)作为输入进行Kalman滤波。采用Sage-Husa方法，根据新获得的观测数据，在线估计出系统噪声特性，分别实时更新第j 个滤波器的系统噪声协方差矩阵
Figure BDA0002111708970000091
和第j个滤波器的测量噪声协方差矩阵
Figure BDA0002111708970000092
保证滤波稳定性，以提高滤波收敛速率和效果。更新参数过程如下：
Figure BDA0002111708970000093
Figure BDA0002111708970000094
Figure BDA0002111708970000095
其中，dn为中间系数，b为遗忘因子，ε(k)为新息矩阵，上标T表示转置，K0j(k)为 k时刻的第j个Kalman滤波器的滤波增益矩阵；
第三步：输出组合：基于似然函数计算k时刻第j模型mj(k)的动态权重向量 Sj(k)，并由此获得多模型滤波器的输出。对n个模型，首先k时刻分别计算第j模型 mj(k)的更新概率μj(k)：
Figure BDA0002111708970000096
其中，r通过以下公式得到：
Figure BDA0002111708970000097
则最终输出的变时钟漂移率状态估计
Figure BDA0002111708970000098
为
Figure BDA0002111708970000099
采用该滤波算法可以得到变时钟漂移率的估计结果，用估计结果对时钟漂移率进行修正和补偿，之后可进行双程信息交换修正时钟偏差。
比较方法1、方法2及本文授时方法的时钟漂移率估计结果。其中，方法1为混合使用TSHL和Tri-message授时协议计算时钟漂移率方法，方法2为在混合使用 TSHL和Tri-message授时协议计算时钟漂移率的基础上，使用不调整滤波参数的 Kalman滤波时钟漂移率估计方法。
设待同步节点初始时钟偏差10μs，时间戳接收抖动分布由0均值、标准差为 15μs的高斯分布引入，时钟的仿真粒度为1μs，标准节点发送时间戳的时间间隔为 1s。
如图3所示，一天内采用本发明授时算法、方法1和方法2对变时钟漂移率的估计结果对比。由图可以看出，在0～4h、8～12h、16～20h等时钟漂移率出现较大程度变化的时间段，本发明的时钟漂移率计算方法可以进行较好的估计和跟踪，估计结果优于方法1和方法2。
如图4所示，为三种不同授时方法累积授时误差对比。由图可以看出，三种方法的授时误差随着时间的推移都会有不同程度的增加，采用本发明的时钟漂移率估计方案对其进行估计和修正，授时精度最高，一天内现有授时协议累积授时误差大约为本文提出的估计方法的1.5倍。
如图5所示，为三种不同授时方法MSE对比。本发明的时钟漂移率估计方法与其他两种方法相比，误差增加速度最小，一天内的授时均方误差可由3.0E-9降低至 5E-10。
本发明还提出一种基于交互式多模型Kalman滤波器的时钟漂移率跟踪系统，所述系统包括：输入的初始状态参数模块、输入交互模块和模型概率更新模块和输出组合模块；
所述输入的初始状态参数模块，用于计算输入的初始状态参数；
所述输入交互模块，用于将初始状态参数进行输入交互，得到混合初始参数；
所述模型概率更新模块模块，用于将混合初始参数输入多模型Kalman滤波器进行滤波；
所述输出组合模块，用于输出变时钟漂移率的估计结果。
本发明还提出一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的方法。
本发明一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行上述的方法。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (4)

1.一种基于交互式多模型滤波器的时钟漂移率跟踪方法，包括：
根据时钟偏差和漂移特性，计算初始状态参数；
将初始状态参数与当前时刻的转移概率进行交互运算，得到混合初始参数；
将混合初始参数输入多模型Kalman滤波器进行估计计算；输出变时钟漂移率的估计结果；
根据变时钟漂移率的估计结果修正系统的时钟偏差；
所述根据时钟偏差和漂移特性，计算初始状态参数，具体包括：
步骤1-1)设在k时刻的时钟离散状态向量c(k)为：
c(k)＝[θ(k),α(k),ρ(k)]T (3)
其中，θ(k)表示k时刻的时钟偏差，α(k)表示k时刻的时钟漂移率，ρ(k)表示k时刻的时钟漂移率变化率；
步骤1-2)k时刻时钟漂移率α(k)为：
α(k)＝α(k-1)+τ(k)×ρ(k)+ωα(k) (4)
其中，α(k-1)为k-1时刻的时钟漂移率，τ(k)表示k时刻采样间隔，ωα(k)表示k时刻时钟漂移率为α(k)时的噪声向量；
步骤1-3)将时钟离散状态向量c(k)转化为k时刻的矩阵向量表达式：
c(k)＝Ac(k-1)+ω(k-1) (5)
其中，A为状态转移矩阵，ω(k-1)为k-1时刻的噪声向量；
步骤1-4)k时刻的观测向量z(k)为：
z(k)＝Hc(k)+v(k) (6)
其中H为观测矩阵，v(k)为观测噪声向量；
步骤1-5)对于输入的n个Kalman滤波器，对应n个估计模型mn，在k时刻的第j个Kalman滤波器对应的估计模型为mj(k)，j＝1,2,...,n；mj(k)的下一个时刻时钟离散状态向量cj(k+1)和当前时刻的观测向量分别为：
cj(k+1)＝Ajcj(k+1)+ωj(k) (7)
zj(k)＝Hjcj(k)+vj(k) (8)
n为自然数，ωj(k)和vj(k)分别表示k时刻第j个估计模型mj(k)的噪声向量和观测噪声向量，其协方差矩阵分别为Q和R，将下一时刻时钟离散状态向量cj(k+1)和当前的观测向量zj(k)作为初始状态参数；
所述将初始状态参数与当前时刻的转移概率进行交互运算，得到混合初始参数，具体包括：
步骤2-1)以μj(k)表示第j个Kalman滤波器对应的估计模型mj(k)在k时刻的转移概率，Pij表示第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)的状态转移矩阵，并设Pij为马尔科夫矩阵；则k-1时刻的模型预测概率
Figure FDA0002552915020000021
为：
Figure FDA0002552915020000022
其中，i＝1,2,...,n；第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)在k-1时刻的转移概率μi/j(k-1|k-1)为：
Figure FDA0002552915020000023
步骤2-2)将初始状态参数在第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)之间进行交互计算后得到初始混合状态参数包括：
Figure FDA0002552915020000024
Figure FDA0002552915020000025
其中，
Figure FDA0002552915020000031
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)的时钟离散状态向量估计，
Figure FDA0002552915020000032
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)输出误差的协方差矩阵估计，
Figure FDA0002552915020000033
为k-1时刻交互计算后的输入到滤波器j的时钟状态向量估计，
Figure FDA0002552915020000034
为k-1时刻交互计算后的输入到滤波器j的初始状态协方差矩阵估计；
所述将混合初始参数输入多模型Kalman滤波器进行滤波；输出变时钟漂移率的估计结果，具体包括：
步骤3-1)将混合初始参数
Figure FDA0002552915020000035
及z(k)输入n个Kalman滤波器，根据协方差矩阵Q和R为滤波参数，分别实时更新滤波参数向量
Figure FDA0002552915020000036
和
Figure FDA0002552915020000037
Figure FDA0002552915020000038
Figure FDA0002552915020000039
Figure FDA00025529150200000310
其中，dn为中间系数，b为遗忘因子，ε(k)为新息矩阵，上标T表示转置，K0j(k)为k时刻的第j个Kalman滤波器的滤波增益矩阵，大小为n*n，下标j表示第j个Kalman滤波器；
步骤3-2)对n个Kalman滤波器，首先分别计算k时刻第j个Kalman滤波器对应的估计模型mj(k)的更新概率μj(k)：
Figure FDA00025529150200000311
其中，Sj(k)为基于似然函数计算的k时刻第j个Kalman滤波器对应的估计模型mj(k)的动态权重向量，r为：
Figure FDA0002552915020000041
步骤3-3)输出变时钟漂移率的状态估计结果
Figure FDA0002552915020000042
Figure FDA0002552915020000043
2.一种基于交互式多模型Kalman滤波器的时钟漂移率跟踪系统，其特征在于，所述系统包括：输入的初始状态参数模块、输入交互模块、多模型Kalman滤波器滤波模块和修正偏差模块；
所述输入的初始状态参数模块，用于根据k时刻发送与接收的时钟偏差和漂移特性，计算初始状态参数；具体包括：
步骤1-1)设在k时刻的时钟离散状态向量c(k)为：
c(k)＝[θ(k),α(k),ρ(k)]T (3)
其中，θ(k)表示k时刻的时钟偏差，α(k)表示k时刻的时钟漂移率，ρ(k)表示k时刻的时钟漂移率变化率；
步骤1-2)k时刻时钟漂移率α(k)为：
α(k)＝α(k-1)+τ(k)×ρ(k)+ωα(k) (4)
其中，α(k-1)为k-1时刻的时钟漂移率，τ(k)表示k时刻采样间隔，ωα(k)表示k时刻时钟漂移率为α(k)时的噪声向量；
步骤1-3)将时钟离散状态向量c(k)转化为k时刻的矩阵向量表达式：
c(k)＝Ac(k-1)+ω(k-1) (5)
其中，A为状态转移矩阵，ω(k-1)为k-1时刻的噪声向量；
步骤1-4)k时刻的观测向量z(k)为：
z(k)＝Hc(k)+v(k) (6)
其中H为观测矩阵，v(k)为观测噪声向量；
步骤1-5)对于输入的n个Kalman滤波器，对应n个估计模型mn，在k时刻的第j个Kalman滤波器对应的估计模型为mj(k)，j＝1,2,...,n；mj(k)的下一个时刻时钟离散状态向量cj(k+1)和当前时刻的观测向量分别为：
cj(k+1)＝Ajcj(k+1)+ωj(k) (7)
zj(k)＝Hjcj(k)+vj(k) (8)
n为自然数，ωj(k)和vj(k)分别表示k时刻第j个估计模型mj(k)的噪声向量和观测噪声向量，其协方差矩阵分别为Q和R，将下一时刻时钟离散状态向量cj(k+1)和当前的观测向量zj(k)作为初始状态参数；
所述输入交互模块，用于将初始状态参数进行输入交互，得到混合初始参数；
具体包括：
步骤2-1)以μj(k)表示第j个Kalman滤波器对应的估计模型mj(k)在k时刻的转移概率，Pij表示第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)的状态转移矩阵，并设Pij为马尔科夫矩阵；则k-1时刻的模型预测概率
Figure FDA0002552915020000051
为：
Figure FDA0002552915020000052
其中，i＝1,2,...,n；第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)在k-1时刻的转移概率μi/j(k-1|k-1)为：
Figure FDA0002552915020000053
步骤2-2)将初始状态参数在第i个Kalman滤波器对应的估计模型mi(k)到第j个Kalman滤波器对应的估计模型mj(k)之间进行交互计算后得到初始混合状态参数包括：
Figure FDA0002552915020000054
Figure FDA0002552915020000061
其中，
Figure FDA0002552915020000062
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)的时钟离散状态向量估计，
Figure FDA0002552915020000063
为k-1时刻第i个Kalman滤波器对应的估计模型mi(k)输出误差的协方差矩阵估计，
Figure FDA0002552915020000064
为k-1时刻交互计算后的输入到滤波器j的时钟状态向量估计，
Figure FDA0002552915020000065
为k-1时刻交互计算后的输入到滤波器j的初始状态协方差矩阵估计；
所述多模型Kalman滤波器滤波模块，用于将混合初始参数输入多模型Kalman滤波器进行滤波；输出变时钟漂移率的估计结果；
具体包括：
步骤3-1)将混合初始参数
Figure FDA0002552915020000066
及z(k)输入n个Kalman滤波器，根据协方差矩阵Q和R为滤波参数，分别实时更新滤波参数向量
Figure FDA0002552915020000067
和
Figure FDA0002552915020000068
Figure FDA0002552915020000069
Figure FDA00025529150200000610
Figure FDA00025529150200000611
其中，dn为中间系数，b为遗忘因子，ε(k)为新息矩阵，上标T表示转置，K0j(k)为k时刻的第j个Kalman滤波器的滤波增益矩阵，大小为n*n，下标j表示第j个Kalman滤波器；
步骤3-2)对n个Kalman滤波器，首先分别计算k时刻第j个Kalman滤波器对应的估计模型mj(k)的更新概率μj(k)：
Figure FDA0002552915020000071
其中，Sj(k)为基于似然函数计算的k时刻第j个Kalman滤波器对应的估计模型mj(k)的动态权重向量，r为：
Figure FDA0002552915020000072
步骤3-3)输出变时钟漂移率的状态估计结果
Figure FDA0002552915020000073
Figure FDA0002552915020000074
所述修正偏差模块，用于根据变时钟漂移率的估计结果修正系统的时钟偏差。
3.一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1所述的方法。
4.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行权利要求1所述的方法。
CN201910574410.XA 2019-06-28 2019-06-28 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统 Active CN110350996B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910574410.XA CN110350996B (zh) 2019-06-28 2019-06-28 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910574410.XA CN110350996B (zh) 2019-06-28 2019-06-28 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统

Publications (2)

Publication Number Publication Date
CN110350996A CN110350996A (zh) 2019-10-18
CN110350996B true CN110350996B (zh) 2020-10-23

Family

ID=68177105

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910574410.XA Active CN110350996B (zh) 2019-06-28 2019-06-28 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统

Country Status (1)

Country Link
CN (1) CN110350996B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110955966A (zh) * 2019-11-25 2020-04-03 北京无线电计量测试研究所 综合原子时数据模拟方法及系统
CN113472318B (zh) * 2021-07-14 2024-02-06 青岛杰瑞自动化有限公司 一种顾及观测模型误差的分级自适应滤波方法及系统
CN114415157A (zh) * 2021-12-30 2022-04-29 西北工业大学 一种基于水声传感器网络的水下目标多模型跟踪方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7705780B1 (en) * 2007-12-20 2010-04-27 The United States Of America As Represented By The Secretary Of The Navy Electronic support measures (ESM) tracking system and method
CN103776453A (zh) * 2014-01-22 2014-05-07 东南大学 一种多模型水下航行器组合导航滤波方法
CN104020480A (zh) * 2014-06-17 2014-09-03 北京理工大学 一种带自适应因子的交互式多模型ukf的卫星导航方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105513091A (zh) * 2015-11-26 2016-04-20 哈尔滨工程大学 一种基于双贝叶斯估计的水下运动体运动状态估计方法
CN107966145B (zh) * 2017-12-21 2020-12-15 中国船舶重工集团公司第七0七研究所 一种基于稀疏长基线紧组合的auv水下导航方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7705780B1 (en) * 2007-12-20 2010-04-27 The United States Of America As Represented By The Secretary Of The Navy Electronic support measures (ESM) tracking system and method
CN103776453A (zh) * 2014-01-22 2014-05-07 东南大学 一种多模型水下航行器组合导航滤波方法
CN104020480A (zh) * 2014-06-17 2014-09-03 北京理工大学 一种带自适应因子的交互式多模型ukf的卫星导航方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Adaptive Clock Skew Estimation with Interactive Multi-Model Kalman Filters for Sensor Networks;Zhe Yang et al.;《2010 IEEE International Conference on Communications》;20100701;1-5 *
Zhe Yang et al..Adaptive Clock Skew Estimation with Interactive Multi-Model Kalman Filters for Sensor Networks.《2010 IEEE International Conference on Communications》.2010, *
组合导航自适应交互多模型算法研究;徐田来等;《系统工程与电子技术》;20081130;第30卷(第11期);2070-2074 *
自适应Kalman滤波的战术数据链自主时问同步算法;顾仁财等;《火力与指挥控制》;20190228;第44卷(第2期);76-79 *

Also Published As

Publication number Publication date
CN110350996A (zh) 2019-10-18

Similar Documents

Publication Publication Date Title
CN110350996B (zh) 基于交互式多模型滤波器的时钟漂移率跟踪方法及系统
Yin et al. A GNSS/5G integrated positioning methodology in D2D communication networks
CN103281772B (zh) 一种无线传感器网络的时间同步方法及系统
Gholami et al. TDOA based positioning in the presence of unknown clock skew
CN102968552B (zh) 一种卫星轨道数据预估与修正方法
US20130337825A1 (en) Space Time Calibration for Networks Using State Model of Node Clock Parameters
CN108668356A (zh) 一种水下传感器时间同步方法
CN110572232B (zh) 一种基于动态响应的免时间戳同步频率偏移跟踪方法
Huang et al. An accurate on-demand time synchronization protocol for wireless sensor networks
CN110940332B (zh) 考虑航天器轨道动态效应的脉冲星信号相位延迟估计方法
CN109557371B (zh) 一种用于配电网相量测量的同步授时守时方法
Wang et al. Timestamp-free clock parameters tracking using extended Kalman filtering in wireless sensor networks
CN114362140B (zh) 一种适用于配电网量测装置的高精度守时方法及装置
CN114297833A (zh) 一种面向分布式仿真的时间同步方法
Koo et al. Novel control theoretic consensus-based time synchronization algorithm for WSN in industrial applications: Convergence analysis and performance characterization
CN111812970B (zh) 一种基于ieee1588协议的双补偿时钟同步方法
Yang et al. A proportional integral estimator-based clock synchronization protocol for wireless sensor networks
Xie et al. A practical clock synchronization algorithm for UWB positioning systems
CN113055114A (zh) 一种基于动态补偿与分级传递式的输油管网时间同步方法
CN116506802A (zh) 面向分布式海洋流场声学测量的水下节点时间同步方法
Ting et al. Clock synchronization in wireless sensor networks: a new model and analysis approach based on networked control perspective
CN114598611B (zh) 面向二集值fir系统事件驱动辨识的输入设计方法及系统
CN109525349A (zh) 一种基于噪声估计与信任加权的分布式时间同步方法
Qi et al. A clock synchronization method for ad hoc networks
CN106211309A (zh) 一种高效的水声传感器网络时间同步方法

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