卫星双向和gps的不同时间频率传递链路的融合方法
技术领域
1.本发明涉及时间频率数据的处理技术领域，具体涉及到卫星双向和gps的不同时间频率传递链路的融合方法。
背景技术：
2.国际权度局(bureau international des poidset measures,bipm)利用不同国家的80多个守时实验室的420多台原子钟通过加权算法计算原子时。这需要一条准确的时间传递链路连接这些原子钟。卫星双向时间频率传递(twstft)和gps精密单点定位(gps ppp)时间传递是当前用于协调世界时(coordinated universal time)产生的两种主要技术。主要守时实验室同时运行着两种技术。twstft是一种双向方式，从1999年以来用于utc产生，是两个地面站通过一颗地球同步卫星传递时间信号。twstft用于时间传递具有准确度高和长期稳定度高的优势。但是分辨率低并且受周日效应扰动的影响。周日效应是twstft链路的主要误差源。gps ppp是一种单项方法，利用国际全球导航卫星系统(gnss)产生的精密卫星轨道和钟差值，得到当地时钟与国际gnss服务系统时间(igs)的差值。gps ppp具有高的短期稳定度和高分辨率的优势。因此twstft和gps ppp融合是有效提高utc时间链路质量理想j9九游会真人的解决方案。
3.随着多种gnss系统的建设，utc时间传递网络变得更加冗余。同时随着高精度时间频率标准的发展，比如光钟和喷泉钟，对远程时间链路的精度提出了越来越高的要求。怎样充分利用冗余链路提高时间传递的准确度，稳定度一直是研究的重点。目前主要的研究工作有对twstft和gps ppp数据加权平均；基于vondrak-cepeck融合平滑算法将twstft数据和gps ppp数据进行融合；将twstft数据作为gps共视最小二乘的附加观测方程。这些研究成果主要是从单尺度方面对数据进行融合。然而实际中，不同时间传递链路系统的比对数据具有不同的分辨率，并且不同时间传递链路的噪声特性不同，通常不同频率段表现出不同类型的噪声。需要同时降低多种噪声的影响提高融合比对结果的准确性和稳定性。因而需要解决多分辨率的融合技术，同时降低不同频率段的不同类型噪声的影响，以便更好地利用不同分辨率数据地互补信息，达到更佳的融合效果。因此，需要提出一种针对twstft和gps ppp的不同时间频率传递链路的融合方法。
技术实现要素：
4.本发明所要解决的技术问题在于克服现有技术的缺点，提供一种使链路数据具有时间频率传递性能好、准确度高、稳定度高的卫星双向和gps的不同时间频率传递链路的融合方法。
5.解决上述技术问题所采用的技术方案是：一种卫星双向和gps的不同时间频率传递链路的融合方法，包括以下步骤：
6.步骤1.数据预处理
7.每间隔30分钟采样卫星双向时间频率传递链路数据twstft，得到链路数据集w，同
时每间隔5分钟采样gps精密单点定位链路数据gps ppp，得到链路数据集g，分别将链路数据集w和链路数据集g中的数据换算成以纳秒为单位，然后对数据进行完整性、异常值检验；
8.步骤2.分别建立链路数据集w与链路数据集g的状态模型和测量模型
9.所述链路数据集w的状态模型为：
[0010][0011][0012][0013]
式中，为链路数据集w中k时刻数据的相位偏差，k为数据计算的时刻，a为常数，为链路数据集w中k时刻数据的频率偏差，为链路数据集w中状态模型k时刻的噪声向量，t为转置符号；
[0014]
所述链路数据集w的测量模型为：
[0015][0016][0017]
式中，为链路数据集w的测量模型k时刻的观测值，为链路数据集w的测量模型k时刻的观测误差；
[0018]
所述链路数据集g的状态模型为：
[0019][0020][0021][0022]
式中，为链路数据集g中k时刻数据的相位偏差，b为常数，为链路数据集g中k时刻数据的频率偏差，为链路数据集g中状态模型k时刻的噪声向量，t为转置符号；
[0023]
所述链路数据集g的测量模型为：
[0024][0025][0026]
式中，为链路数据集g的测量模型k时刻的观测值，为链路数据集g的测量模型k时刻的观测误差；
[0027]
步骤3.利用小波分解法对链路数据集w的状态模型和测量模型与链路数据集g的状态模型和测量模型进行尺度分解，将链路数据集w与链路数据集g分解到相同的分辨率，得到低分辨率上的链路数据集w的状态模型和测量模型及低分辨率上的链路数据集g的状态模型和测量模型；
[0028]
步骤4.在相同分辨率上对链路数据集w与链路数据集g按照下式进行卡尔曼滤波融合
[0029][0030][0031][0032][0033][0034][0035]
式中，为k时刻的一步预测结果，φ
k,k-1
为k时刻状态转移矩阵，为k-1时刻融合结果，为k时刻融合结果，p
k/k-1
为k时刻的一步预测协方差矩阵，p
k-1
为k-1时刻的误差协方差矩阵，q
k-1
为k-1时刻模型噪声协方差矩阵，为链路数据集w的卡尔曼滤波增益矩阵，为链路数据集w的测量模型的观测值，为链路数据集g的卡尔曼滤波增益矩阵，为数据集g的测量模型的观测值，pk为k时刻的误差协方差矩阵，为链路数据集w的协方差矩阵，为链路数据集g的协方差矩阵，i为单位矩阵；
[0036]
步骤5.将步骤4得到的卡尔曼滤波融合结果利用小波分解系数重构方法得到原始数据融合后的结果即原始twstft链路和gps ppp链路数据的融合结果；
[0037]
所述小波分解系数重构方法为：
[0038][0039]
式中，xj(m)为第m个尺度系数，j为尺度，m为尺度系数的排序，r为正整数，t为小波系数的排序，h()为第一滤波器，g()为第二滤波器，y
j-1
(t)为第t个小波系数；
[0040]
步骤6.按照下式得到twstft链路和gps ppp链路融合后的增益，用于判断融合后链路噪声的降低情况；
[0041]
q＝(σ
a-σb)/σa[0042]
式中，σa是原始twstft链路和gps ppp链路数据的标准偏差，σb是融合数据的标准偏差，a为常数，b为常数；
[0043]
步骤7.按照下式得到twstft链路和gps ppp链路融合后的allan方差，用于判断融合后链路的稳定度；
[0044][0045]
p2＝m 1-2n
[0046]
式中，为allan方差，y为常数，x
k 2n
为k 2n时刻钟差，x
k n
为k n时刻钟差，xk为k时刻钟差，m为测量值的个数，n为取样周期，τ0为基本采样间隔。
[0047]
作为一种优选的技术方案，所述步骤3中小波分解法为daubechies 2小波分解法。
[0048]
作为一种优选的技术方案，所述数据完整性检验方法为：依据数据比对日期后项
减去前项，判断有无跳变，有跳变则是数据有缺失，反之，没有缺失，对于缺失的数据运用最近邻差值法或线性插值法进行补充；
[0049]
所述异常值检验方法为：采用莱特准则进行检验并依据最近邻差值法插值替换。
[0050]
本发明的有益效果如下：
[0051]
本发明在考虑不同链路数据集的特点的基础上，结合小波变换的分解算法将不同采样率的链路数据分解到相同的尺度进行融合处理，相较于仅用单一链路的数据提升时间频率传递的可靠性。
[0052]
本发明在不同链路数据集进行融合处理时，利用卡尔曼滤波对小波分解尺度的数据进行融合，依据本发明的算法同时降低了链路数据中多类噪声的影响，提高了链路数据的短期稳定性和长期稳定性。
[0053]
最后通过融合结果相对于链路数据的增益得出融合链路降低了链路数据集w的周日效应和链路数据集g的误差。相较于单一链路数据，提高时间频率传递的准确性。
附图说明
[0054]
图1是本发明的流程示意图。
[0055]
图2是本发明关于链路数据集w和链路数据集g的小波三尺度分解的组织形式图。
[0056]
图3是链路数据集w和链路数据集g基于卡尔曼滤波的数据融合流程图。
[0057]
图4是是融合结果与twstft和gps ppp的对比图。
[0058]
图5是是融合结果的allan偏差与twstft的allan偏差和gps ppp的allan偏差的对比图。
具体实施方式
[0059]
下面结合附图和实施例对本发明进一步详细说明，但本发明不限于下述的实施方式。
[0060]
在图1、2、3中，本实施例的卫星双向和gps的不同时间频率传递链路的融合方法，其特征在于，包括以下步骤：
[0061]
步骤1.数据预处理
[0062]
每间隔30分钟采样卫星双向时间频率传递链路数据twstft，得到链路数据集w，同时每间隔5分钟采样gps精密单点定位链路数据gps ppp，得到链路数据集g，分别将链路数据集w和链路数据集g中的数据换算成以纳秒为单位，然后对数据进行完整性、异常值检验；
[0063]
所述数据完整性检验方法为：依据数据比对日期后项减去前项，判断有无跳变，有跳变则是数据有缺失，反之，没有缺失，对于缺失的数据运用最近邻差值法或线性插值法进行补充；
[0064]
所述异常值检验方法为采用莱特准则进行检验并依据最近邻差值法插值替换，具体过程为先对链路数据进行拟合得到拟合值，以实测值与拟合值的残差值为基础值，计算方差值，并得到标准值，如果在tj时刻，满足式|xi(tj)-si(tj)|＞3σ，xi(tj)为tj时刻的实测值，si(tj)为tj时刻的拟合值，σ为标准值，则xi(tj)为检测出的异常值，随后对拟合值进行插值替代异常值。
[0065]
步骤2.分别建立链路数据集w与链路数据集g的状态模型和测量模型
[0066]
所述链路数据集w的状态模型为：
[0067][0068][0069][0070]
式中，为链路数据集w中k时刻数据的相位偏差，k为数据计算的时刻，a为常数，为链路数据集w中k时刻数据的频率偏差，为链路数据集w中状态模型k时刻的噪声向量，t为转置符号；
[0071]
所述链路数据集w的测量模型为：
[0072][0073][0074]
式中，为链路数据集w的测量模型k时刻的观测值，为链路数据集w的测量模型k时刻的观测误差；
[0075]
所述链路数据集g的状态模型为：
[0076][0077][0078][0079]
式中，为链路数据集g中k时刻数据的相位偏差，b为常数，为链路数据集g中k时刻数据的频率偏差，为链路数据集g中状态模型k时刻的噪声向量，t为转置符号；
[0080]
所述链路数据集g的测量模型为：
[0081][0082][0083]
式中，为链路数据集g的测量模型k时刻的观测值，为链路数据集g的测量模型k时刻的观测误差；
[0084]
步骤3.利用daubechies 2小波分解法对链路数据集w的状态模型和测量模型与链路数据集g的状态模型和测量模型进行尺度分解，将链路数据集w与链路数据集g分解到相同的分辨率，得到低分辨率上的链路数据集w的状态模型和测量模型及低分辨率上的链路数据集g的状态模型和测量模型；
[0085]
步骤4.在相同分辨率上对链路数据集w与链路数据集g按照下式进行卡尔曼滤波融合
[0086][0087]
[0088][0089][0090][0091][0092]
式中，为k时刻的一步预测结果，φ
k,k-1
为k时刻状态转移矩阵，为k-1时刻融合结果，为k时刻融合结果，p
k/k-1
为k时刻的一步预测协方差矩阵，p
k-1
为k-1时刻的误差协方差矩阵，q
k-1
为k-1时刻模型噪声协方差矩阵，为链路数据集w的卡尔曼滤波增益矩阵，为链路数据集w的测量模型的观测值，为链路数据集g的卡尔曼滤波增益矩阵，为数据集g的测量模型的观测值，pk为k时刻的误差协方差矩阵，为链路数据集w的协方差矩阵，为链路数据集g的协方差矩阵，i为单位矩阵；
[0093]
步骤5.将步骤4得到的卡尔曼滤波融合结果利用小波分解系数重构方法得到原始数据融合后的结果即原始twstft链路和gps ppp链路数据的融合结果；
[0094]
所述小波分解系数重构方法为：
[0095][0096]
式中，xj(m)为第m个尺度系数，j为尺度，m为尺度系数的排序，r为正整数，t为小波系数的排序，h()为第一滤波器，g()为第二滤波器，y
j-1
(t)为第t个小波系数；
[0097]
步骤6.按照下式得到twstft链路和gps ppp链路融合后的增益，用于判断融合后链路噪声的降低情况；
[0098]
q＝(σ
a-σb)/σa[0099]
式中，σa是原始twstft链路和gps ppp链路数据的标准偏差，σb是融合数据的标准偏差，a为常数，b为常数；
[0100]
步骤7.按照下式得到twstft链路和gps ppp链路融合后的allan方差，用于判断融合后链路的稳定度，其结果如图5；
[0101][0102]
p2＝m 1-2n
[0103]
式中，为allan方差，y为常数，x
k 2n
为k 2n时刻钟差，x
k n
为k n时刻钟差，xk为k时刻钟差，m为测量值的个数，n为取样周期，τ0为基本采样间隔。
[0104]
本发明在不同链路数据集进行融合处理时，利用卡尔曼滤波对小波分解尺度的数据进行融合，本发明同时降低了链路数据中多类噪声的影响，提高了链路数据的短期稳定性和长期稳定性。最后通过融合结果相对于链路数据的增益得出融合链路降低了链路数据集w的周日效应和链路数据集g的误差，如图4和图5，相较于单一链路数据，提高时间频率传递的准确性。