CN114545764B - 一种双重滤波加权平均的星站钟差数据融合方法 - Google Patents

Abstract

本发明公开了一种双重滤波加权平均的星站钟差数据融合方法，采用前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值；用动态加权平均算法实现钟差数据融合；采用α‑β滤波器进行星站钟差跟踪滤波，生成钟差融合终值，并为加权平均时的动态权值分配提供依据，同时还为前置滤波器的阈值调优提供参考。该方法适用于存在奇异值和空缺值的卫星钟差观测数据处理，能够获得标准差较低的星站钟差融合结果，即钟差波动程度降低，使融合结果更接近钟差真实值。

Description

一种双重滤波加权平均的星站钟差数据融合方法

技术领域

本发明属于卫星共视授时技术领域，具体涉及一种双重滤波加权平均的星站钟差数据融合方法。

背景技术

时间为一切时序过程的测量和定量研究提供必不可少的时基坐标。不仅在天体物理、量子信息等基础研究领域广泛应用，而且在信息通信、电力传输、公共交通、航空航天等工程技术领域同样发挥重要作用。各领域对高精度时间标准的需求不断对授时技术提出新的要求。人造卫星的出现促使卫星授时技术形成并快速发展。卫星授时的信号覆盖范围大、传递精度高、传播衰减小，是目前精度最高的授时方法之一。在新一代基于GNSS(GlobalNavigation Satellite System)的卫星授时系统中，单个基准站通常配置多个基准终端。这些基准终端同时观测GNSS卫星群，通过测算伪距得到本地参考时间与卫星钟之间的偏差，即星站钟差。由于观测位置和角度不同等原因，不同的基准终端在同一时刻观测到的卫星及其数量不尽相同。对同一基准站中不同基准终端测得的星站钟差进行数据融合处理，是保证授时精度的重要前提。而对这些基准终端的观测数据进行融合处理时，数据源具有以下特点：

(1)每个基准终端的观测数据通常包含多颗卫星的钟差数据。

(2)每一颗卫星的钟差观测数据中不可避免会含有奇异值和空缺值。星站钟差中的奇异值会影响加权平均数据融合的精度，而星站钟差中空缺值会对跟踪滤波环节造成影响，必须在此之前进行填补。

(3)多个基准终端观测数据中的同一颗卫星的钟差才可以进行数据融合处理。

发明内容

综合考虑共视授时基准站的卫星观测数据特点和实际需求，本发明的目的在于，提供一种双重滤波加权平均的星站钟差数据融合方法。

为了实现上述任务，本发明采取如下的技术解决方案：

一种双重滤波加权平均的星站钟差数据融合方法，其特征在于，该方法采用前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值；用动态加权平均算法实现钟差数据融合；采用α-β滤波器进行星站钟差跟踪滤波，生成钟差融合终值，并为加权平均时的动态权值分配提供依据，同时还为前置滤波器的阈值调优提供参考。

根据本发明，所述前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值的具体实现方法是：

对输入星站钟差序列D_k，前置滤波器的输出响应如下：

式中，M_k是当前星站钟差样本及其前K-1个样本组成的滑动窗口的中值，定义为：

M_k＝median{D_k-(K-1),…D_k-2,D_k-1,D_k} (2)

其中，median是求给定数值的中值的函数；K是一个正整数，称为窗口宽度；S_k是中位数绝对偏差尺度估计，定义为：

S_k＝1.4826×median_j∈[1,K]{|D_k-j-M_k|} (3)

其中，1.4826是一个工程经验值，它使中位数绝对偏差尺度估计成为高斯数据标准偏差的无偏估计；T是一个动态阈值调优参数，其大小由跟踪滤波器输出最新的钟差融合终值的方差大小决定；当前星站钟差样本D_k与中值M_k相差超过TS_k，则视该样本为奇异值，并用中值M_k替换该样本；

同时，前置滤波器的滑动窗口可将钟差数据源中的空缺值与邻近的非空数值联系起来，只要滑动窗口内不全是空缺值，即只需滑动窗口内有一个及以上非空数值，则当前正被处理的空缺值可被填补为窗口内数据的中值；而对于持续时间大于窗口长度的空缺值，取同一时刻其它基准终端观测同一颗卫星的非空钟差数据的平均值填补。

进一步地，所述用动态加权平均算法实现钟差数据融合的方法是：

将不同基准站观测到的同一颗卫星钟差数据进行分类融合，进而得到地面观测站与每颗卫星之间更精确的星站钟差；

来自N个基准终端的N组原始数据的融合过程可用公式(4)来描述：

其中，

表示k时刻的融合值，它是一个估计值；^代表估计，Y_ki表示各基准终端k时刻通过前置滤波之后的星站钟差数据，a_ki表示k时刻各基准终端的加权系数，并满足下式(5)的关系：

如果该基准终端的观测数据与跟踪滤波器输出的滤波值相差大，则该基准终端观测数据的均方根误差大，可靠度小，分配的加权系数就小；否则认为可靠度高，分配到的加权系数就大；

式(6)用以计算基准终端观测数据的均方根误差：

其中，RMSE代表均方根误差，Y_k代表k时刻的观测数据，X_k-1代表k-1时刻的跟踪滤波值，n代表测量数据的个数；

各个基准终端的观测数据可看成是由真实值叠加了各种测量噪声得到的，即：

其中，

表示k时刻第i个基准终端的观测值，

表示k时刻第i个基准终端的观测噪声，Y_k视为待估计的真值；使估计值与真值的偏差平方和的均值达到最小值，即满足下述指标：

将式(4)代入到式(8)中，得到：

由于加权系数还满足式(5)，那么联合式(9)，应用拉格朗日乘数法求解加权系数，引入拉格朗日函数函数L：

对函数中的a₁，a₂，…a_N分别求偏导，并令其等于零，最终可求得各基准终端加权系数的计算公式：

其中，σ_i取各基准终端观测数据的均方根误差RMSE。

优选地，所述α-β滤波器包含预测和更新两个步骤，用下面的两个公式描述：

式中，-代表先验，^代表估计；X_k(X_k∈Rⁿ)：系统状态变量；Z_k(Z_k∈R^m)：观测变量；F：状态转移矩阵；H：观测矩阵；K：常增益矩阵；

系数α和β是无量纲的量，分别为状态向量不同分量的常滤波增益，α和β的值是一个很小的正值，满足以下关系：

0.3＜α＜0.5 (14)

将α代入式(15)可得β的范围是：0.05＜β＜0.18；α和β的值确定之后，α-β滤波器增益就可以按照下面公式确定：

式中，T为滤波周期。

本发明的双重滤波加权平均的星站钟差数据融合方法，适用于存在奇异值和空缺值的卫星钟差观测数据处理，能够获得标准差较低的星站钟差融合结果，即钟差波动程度降低，使融合结果更接近钟差真实值。

附图说明

图1是本发明的双重滤波加权平均的星站钟差数据融合方法流程图。

图2是授时基准站中基准终端1针对同一颗卫星的钟差观测数据；

图3是授时基准站中基准终端2针对同一颗卫星的钟差观测数据；

图4是授时基准站中基准终端3针对同一颗卫星的钟差观测数据；

图5是授时基准站中基准终端4针对同一颗卫星的钟差观测数据；

图6是授时基准站中基准终端5针对同一颗卫星的钟差观测数据；

图7是对五个基准终端观测同一颗卫星的钟差数据的融合处理结果。

下面结合附图和实施例对本发明进行进一步地详细说明。

具体实施方式

本实施例给出一种双重滤波加权平均的星站钟差数据融合方法，流程如图1所示，主要包括三个核心部分：星站钟差前置滤波、星站钟差加权平均和星站钟差跟踪滤波。具体地，设计前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值；设计动态加权平均算法实现钟差数据融合；采用α-β滤波器进行星站钟差跟踪滤波，生成钟差融合终值，并为加权平均时的动态权值分配提供依据，同时还为前置滤波器的阈值调优提供参考。

具体实现的步骤如下：

(1)前置滤波器进行星站钟差的奇异值剔除和空缺值填补。具体来说，对输入星站钟差序列D_k，前置滤波器的输出响应如下：

式中，M_k是当前星站钟差样本及其前K-1个样本组成的滑动窗口的中值，定义为：

M_k＝median{D_k-(K-1),…D_k-2,D_k-1,D_k} (2)

其中，median是求给定数值的中值的函数；K是一个正整数，称为窗口宽度；S_k是中位数绝对偏差尺度估计，定义为：

S_k＝1.4826×median_j∈[1,K]{|D_k-j-M_k|} (3)

其中，1.4826是一个工程经验值，它使中位数绝对偏差尺度估计成为高斯数据标准偏差的无偏估计；T是一个动态阈值调优参数，其大小由跟踪滤波器输出最新的钟差融合终值的方差大小决定。当前星站钟差样本D_k与中值M_k相差超过TS_k，则视该样本为奇异值，并用中值M_k替换该样本。

同时，前置滤波的滑动窗口可将钟差数据源中的空缺值与邻近的非空数值联系起来。只要滑动窗口内不全是空缺值，即只需滑动窗口内有一个及以上非空数值，则当前正被处理的空缺值可被填补为窗口内数据的中值。而对于持续时间大于窗口长度的空缺值，取同一时刻其它基准终端观测同一颗卫星的非空钟差数据的平均值填补。

(2)星站钟差动态加权平均是将不同基准站观测到的同一颗卫星钟差数据进行分类融合，进而得到地面观测站与每颗卫星之间更精确的星站钟差。来自N个基准终端的N组原始数据的融合过程可用公式(4)来描述：

其中，

表示k时刻的融合值，它是一个估计值(^代表估计)，Y_ki表示各基准终端k时刻通过前置滤波之后的星站钟差数据，a_ki表示k时刻各基准终端的加权系数，并满足式(5)的关系：

如果该基准终端的观测数据与跟踪滤波器输出的滤波值相差大，则该基准终端观测数据的均方根误差大，可靠度小，分配的加权系数就小；否则认为可靠度高，分配到的加权系数就大。式(6)用以计算基准终端观测数据的均方根误差：

其中，RMSE代表均方根误差，Y_k代表k时刻的观测数据，X_k-1代表k-1时刻的跟踪滤波值，n代表测量数据的个数。各个基准终端的观测数据可以看成是由真实值叠加了各种测量噪声得到的，即：

其中，

表示k时刻第i个基准终端的观测值，

表示k时刻第i个基准终端的观测噪声，Y_k视为待估计的真值。使估计值与真值的偏差平方和的均值达到最小值(min)，即满足下述指标：

将式(4)代入到式(8)中，得到：

由于加权系数还满足式(5)，那么联合式(9)，应用拉格朗日乘数法求解加权系数，引入拉格朗日函数函数L：

对函数中的a₁，a₂，…a_N分别求偏导，并令其等于零，最终可求得各基准终端加权系数的计算公式：

其中，σ_i取各基准终端观测数据的均方根误差RMSE。

(3)采用α-β滤波器进行星站钟差跟踪滤波，生成钟差融合终值，并为加权平均时的动态权值分配提供依据，同时还为前置滤波器的阈值调优提供参考。

所述α-β滤波器包含预测和更新两个步骤，可以用下面的两个公式描述：

式中，-代表先验，^代表估计；X_k(X_k∈Rⁿ)：系统状态变量；Z_k(Z_k∈R^m)：观测变量；F：状态转移矩阵；H：观测矩阵；K：常增益矩阵。

系数α和β是无量纲的量，分别为状态向量不同分量的常滤波增益，α和β的值对滤波效果的影响比较大，且需要根据实际的应用背景和数据特征来确定。这两个系数一旦确定，增益就是个确定的量。为了保证系统的稳定性，α和β的值是一个很小的正值，满足以下关系：

0.3＜α＜0.5 (14)

将α代入式(15)可得β的范围是：0.05＜β＜0.18。α和β的值确定之后，滤波器增益就可以按照下面公式确定：

公式(16)中，T为滤波周期。

实验例：

某授时基准站有五个基准终端，共同观测同一颗卫星、连续观测500次获得的原始星站钟差数据分别如图2、图3、图4、图5和图6所示。

图7是采用本实施例所提出的双重滤波加权平均的星站钟差数据融合方法，对五组基准终端的星站钟差数据融合处理结果。与原始星站钟差观测数据相比，在存在钟差奇异值和空缺值的前提下，该方法获得了标准差较低的星站钟差融合结果，即钟差波动程度降低，使融合结果更接近钟差真实值。

Claims (4)

1.一种双重滤波加权平均的星站钟差数据融合方法，其特征在于，该方法采用前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值；用动态加权平均算法实现钟差数据融合；采用α-β滤波器进行星站钟差跟踪滤波，生成钟差融合终值，并为加权平均时的动态权值分配提供依据，同时还为前置滤波器的阈值调优提供参考。

2.如权利要求1所述的方法，其特征在于，所述前置滤波器剔除各个基准终端观测钟差数据中的奇异值，并填补钟差数据空缺值的具体实现方法是：

对输入星站钟差序列D_k，前置滤波器的输出响应如下：

式中，M_k是当前星站钟差样本及其前K-1个样本组成的滑动窗口的中值，定义为：

M_k＝median{D_k-(K-1)，…D_k-2，D_k-1，D_k} (2)

其中，median是求给定数值的中值的函数；K是一个正整数，称为窗口宽度；S_k是中位数绝对偏差尺度估计，定义为：

S_k＝1.4826×median_j∈[1，K]{|D_k-h-M_k|} (3)

其中，1.4826是一个工程经验值，它使中位数绝对偏差尺度估计成为高斯数据标准偏差的无偏估计；T是一个动态阈值调优参数，其大小由跟踪滤波器输出最新的钟差融合终值的方差大小决定；当前星站钟差样本D_k与中值M_k相差超过TS_k，则视该样本为奇异值，并用中值M_k替换该样本；

同时，前置滤波器的滑动窗口可将钟差数据源中的空缺值与邻近的非空数值联系起来，只要滑动窗口内不全是空缺值，即只需滑动窗口内有一个及以上非空数值，则当前正被处理的空缺值可被填补为窗口内数据的中值；而对于持续时间大于窗口长度的空缺值，取同一时刻其它基准终端观测同一颗卫星的非空钟差数据的平均值填补。

3.如权利要求1所述的方法，其特征在于，所述用动态加权平均算法实现钟差数据融合的方法是：

将不同基准站观测到的同一颗卫星钟差数据进行分类融合，进而得到地面观测站与每颗卫星之间更精确的星站钟差；

来自N个基准终端的N组原始数据的融合过程可用下式(4)来描述：

其中，

表示k时刻的融合值，它是一个估计值；^代表估计，Y_ki表示各基准终端k时刻通过前置滤波之后的星站钟差数据，a_ki表示k时刻各基准终端的加权系数，并满足下式(5)的关系：

如果该基准终端的观测数据与跟踪滤波器输出的滤波值相差大，则该基准终端观测数据的均方根误差大，可靠度小，分配的加权系数就小；否则认为可靠度高，分配到的加权系数就大；

式(6)用以计算基准终端观测数据的均方根误差：

其中，RMSE代表均方根误差，Y_k代表k时刻的观测数据，X_k-1代表k-1时刻的跟踪滤波值，n代表测量数据的个数；

各个基准终端的观测数据可看成是由真实值叠加了各种测量噪声得到的，即：

其中，

表示k时刻第i个基准终端的观测值，

表示k时刻第i个基准终端的观测噪声，Y_k视为待估计的真值；使估计值与真值的偏差平方和的均值达到最小值，即满足下述指标：

将式(4)代入到式(8)中，得到：

由于加权系数还满足式(5)，那么联合式(9)，应用拉格朗日乘数法求解加权系数，引入拉格朗日函数函数L：

对函数中的a₁，a₂，...a_N分别求偏导，并令其等于零，最终可求得各基准终端加权系数的计算公式：

其中，σ_i取各基准终端观测数据的均方根误差RMSE。

4.如权利要求1所述的方法，其特征在于，所述α-β滤波器包含预测和更新两个步骤，用下面的两个公式描述：

式中，-代表先验，^代表估计；X_k(X_k∈Rⁿ)：系统状态变量；Z_k(Z_k∈R^m)：观测变量；F：状态转移矩阵；H：观测矩阵；K：常增益矩阵；

系数α和β是无量纲的量，分别为状态向量不同分量的常滤波增益，α和β的值是一个很小的正值，满足以下关系：

0.3＜α＜0.5 (14)

将α代入式(15)可得β的范围是：0.05＜β＜0.18；α和β的值确定之后，α-β滤波器增益按照下面公式确定：

式中，T为滤波周期。

Images