基于奇异谱分析的BDS 卫星钟差周期项提取

雷 雨，赵丹宁

（1.西安邮电大学 计算机学院，西安 710121；2.宝鸡文理学院 电子电气工程学院，宝鸡 721016）

卫星导航系统的定位、导航和授时服务依赖于准确的时间基准，目前北斗卫星导航系统（Beidou Navigation Satellite System,BDS）时间基准的建立和保持由地面监测站原子钟通过综合原子时算法计算产生，这种模式过度依赖地面站，一旦地面监测站出现异常，系统时间基准就会不连续。为保证系统时间基准的连续性和稳定性，必须降低对地面监测站的依赖。随着BDS 卫星数量与星载钟性能的不断提高，利用星载钟建立和保持星座自主时间基准，对提高系统的生存能力具有重要的现实意义[1,2]。

受星载钟自身因素和太空环境的影响，卫星钟差序列一般含有趋势项、周期项和随机项[3]，其中，星载钟受太空环境影响的主要表现为钟差的周期波动，导致星载钟在轨性能逊于地面钟。卫星钟差的周期波动不仅和其轨道周期相关，而且和日食周期、测量性能等有关[4,5]。在未扣除卫星钟差的周期波动之前，不应将星载钟纳入星座自主时间基准的建立和保持，以避免将周期波动引入系统时间。现有的卫星钟差特征分量提取方法有多项式拟合、傅里叶变换、小波变换和经验模态分解等，Huang 等构建多项式模型拟合提取钟差周期项，提升了BDS 卫星钟差预报精度[3]；李骁逸等基于傅里叶变换的频谱分析方法校正BDS 卫星钟的周期波动，校正后BDS 各类星载钟的万秒稳均得到提高[6]；雷雨等利用小波变换将钟差序列分解成具有不同频率特征的信号分量，并根据各分量的特征构建预报模型，改进了钟差预报精度[7]；梁益丰等提出以完备集合经验模态分解为基础的钟差信噪分离方法，并将其应用于BDS 卫星钟差周期项识别，扣除周期项后卫星钟万秒稳获得提升[8]。上述方法虽然在钟差特征分量分离和提取方面取得了较好的效果，但均存在一定的局限性。传统的多项式拟合方法无法提取钟差序列中的时变信号，提取的钟差趋势项、周期项和随机项不够准确；傅立叶变换将信号从时间域转化到频率域，以分离信号中的周期分量，但该方法不适用于非线性、非平稳信号；经验模态分解能将非线性、非平稳信号分解为一系列固有模态分量，并根据不同模态分量的特征识别周期成分，但这类方法存在端部效应和模态混叠的问题。

奇异谱分析（Singular Spectrum Analysis,SSA）是20 世纪90 年代兴起的一种研究非线性、非平稳信号的有效方法。根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解和重构，从而提取出代表原时间序列不同成分的信号，如长期趋势信号、周期信号和噪声信号，从而对时间序列的波动特征进行分析[9,10]。SSA 方法提取信号不需要先验信息，不受正弦波假定的约束，能够较好地从含噪声的时间序列中提取时变信号，目前已在GPS 坐标时间序列分析、多路径效应消除、钟差分析和预报[11,12]等方面取得了成功应用。其中，肖胜红等[11]提出SSA 和傅立叶带通滤波器相结合的卫星钟差周期项提取方法，有效抑制了重建周期项时频率混叠和边界效应；Xue 等[12]利用SSA 分离卫星钟差的趋势项和随机项，提高了钟差短期预报精度。上述研究表明，SSA 可以较好地从复杂卫星钟差序列中提取趋势项和周期项信息。

本文在分析卫星钟差特性的基础上，提出一种基于SSA 的“分解—辨识—重建”钟差周期项提取方法。首先应用SSA 对钟差序列进行分解，获得时间序列的多个重构成分，然后利用重标极差分析（Rescaled Range Analysis,R/S）方法计算各重构成分的Hurst 指数，最后根据Hurst 指数对各重构成分的波动特征进行辨识，进而完成钟差趋势项、周期项和随机项的分离和提取。利用BDS 卫星钟差数据分析和检验了本文方法的有效性与实用性。

1 星载原子钟物理特性

1.1 星载原子钟物理模型

受星载原子钟自身因素及空间外界环境的影响，星载原子钟信号不仅存在趋势性变化，而且表现出复杂的周期性波动。星载原子钟模型通常可表示为[13]：

其中，x(t) 为t时刻原子钟相位；x0为原子钟初始相位；y0为原子钟钟速；d为原子钟频漂；A、f0和φ分别为原子钟相位周期分量的振幅、频率和初始相位；σ1W1(t)和表示两种起主导作用的原子钟噪声，W1(t) 和W2(s) 分别为两个独立的维纳过程，σ1和σ2分别为两个维纳过程的扩散系数，用来表示两种原子钟噪声强度；σ3ε(t) 为原子钟测量噪声，σ3表示测量噪声强度。与趋势分量相比，周期分量和随机分量的数量级较小，但对频率稳定度的影响会逐渐累积。

1.2 各分量对频率稳定度的影响

阿伦方差与原子钟频漂、相位白噪声、维纳过程的扩散系数、周期分量之间的关系可表示为[13]：

2 基于SSA 的卫星钟差周期项分析方法

基于SSA 的卫星钟差周期项分析方法主要包括钟差数据预处理、钟差分解与重构以及钟差重构成分辨识等过程，图1 给出卫星钟差周期项分析流程。

已知卫星钟差序列为{x1,x2,…,xN}，钟差周期项分析的具体步骤为：

1）钟差数据预处理

运用中位数法检测和剔除钟差序列中的粗差，若钟差的一次差分数据Δxi(i=1,2…N-1)满足式(3)即认为Δxi为粗差[14]：

其中，m表示钟差一次差分序列的中位数；n表示检测阈值，为保证粗差检测的准确性和可靠性，本文取n=6。对于粗差点，利用线性内插法对其进行插补。

2）钟差序列分解

利用SSA 对数据预处理后的钟差序列进行分解，主要计算过程为[9,10]：

②对角平均。计算时滞矩阵X在Uj上的投影：

其中，ai,j称为时间主成分。根据时间正交函数与时间主成分进行对角平均，得到钟差序列的重构成分（Reconstructed Component,RC）：

其中，M*=min(M,N-M+1)，K*=max(M,N-M+1)。

③分组。卫星钟差序列中通常包含多种周期性变化、趋势性变化和随机性变化，利用w-correlation 分析RC 成分之间的相关性，将信号特征相似的RC 分组。将不同RC 分别用Z(i)、Z(j)表示，则重构分量之间的w-correlation 可表示为[15]：

其中，wk为权重系数，其定义为wk=min(k,M,N-k)；ρi,j的绝对值越接近于1，说明Z(i)、Z(j)两者之间的相关性越强。根据经验，ρi,j＞0.6即认为两者之间存在相关性，因此，将RC成分之间相关性大于0.6 的两者视为是同一组信号。

3）RC 成分辨识

一个时间序列在一段时间内的波动如何随时间跨度大小而变化往往可以揭示该时间序列的特性，对于卫星钟差的趋势项和周期项，其时间序列当前或过去的取值以远超随机扰动所能达到的程度影响该序列在未来的取值，统计学上称为时间序列存在长期记忆性（Long-term Memory）和长时间相关效应，而随机项的未来值和当前或过去的取值相关性不强，长期记忆性较弱，属于均值回复过程。由英国水利学家赫尔斯特提出Hurst 指数是用来衡量时间序列是否有长期记忆性的一个指标[16,17]，体现了时间序列的自相关性，尤其能够反映时间序列中隐藏的长期趋势。本文采用Hurst 指数H来辨识RC 成分序列特征，以定性分析钟差序列中的信号主导分量和噪声主导分量。若RC成分的Hurst 指数0.5＜H＜1，则RC 成分属于趋势分量或周期信号主导分量，且H越大，规律性越强；若0＜H≤0 .5，则判定RC 成分为随机噪声主导分量，且H越小，随机性越强。

基于R/S 分析法计算Hurst 指数的原理为[16,17]：对于时间序列{x1,x2…xN}，将其分为C个长度为L的等长子区间，L=N/C，则累积离差为：

其中，ec为第c个子区间序列的平均值，c=1,2…C。令第c个子区间序列的标准差为Sc，则极差Rc为：

定义重标极差γL为：

其中，H为Hurst 指数。对L和γL进行双对数线性回归拟合，回归方程的截距就是式(11)中的常数O，而斜率就是H。

3 算例分析

采用德国地球科学研究中心（Deutsches Geoforschungs Zentrum,GFZ）提供的BDS 精密卫星钟差数据（ftp://ftp.gfz-potsdam.de/GNSS/products/mgex），选取钟差序列连续性较好的C05、C09、C11、C43 星载钟，包括GEO（C05）、IGSO（C09）与MEO（C11与C43）轨道类型，铷原子钟（C05、C09、C11）与氢原子钟（C43）类型，选择2022 年7 月29 日至8月11 日共14 天4032 个数据点，数据采样间隔为5 min。利用SSA 对时间序列进行分析，重点是确定窗口长度M，一般1＜M＜N/2且为周期的最小公倍数。BDS 卫星钟差序列中通常存在6 h、12 h 和24 h的周期[18-20]，因此，选择窗口长度M为卫星钟差周期的最小公倍数，即M=288。

3.1 钟差序列分解和重建试验

由于篇幅限制，以C05 卫星钟差为例展示钟差序列分解和重建的过程及效果。卫星钟差数据经SSA 分解获得288 个RC 序列，频率按照RC 阶次依次从低到高排列。原始钟差序列与所有RC 之和的残差均方根为1.17×10-19s，表明原始钟差序列与重建序列之差接近0，证明SSA 用于钟差序列分解和重建是可行的。C05 卫星钟差前20 阶RC 成分之间的w-correlation 分析结果如图2 所示。根据SSA 分组原理，当RC 之间的相关系数大于一定的阈值，则判断属于同一周期信号。若阈值选择过大，会导致一些周期信号被过滤掉，反之则不能较好地排除一些周期性不明显的RC 成分。经实验，本文选择相关性阈值为0.6。从图2 可以看出，当阶次大于14 时，各RC 之间就不能很好地相互分离，说明非周期变化占较大部分。

图2 前20 阶RC 成分的w-correlationFig.2 W-correlation of the first 20 RCs

根据w-correlation 分析结果，将ρi,j＞0.6的RC成分合并为同一周期信号，C05 卫星钟差的分析结果如图3 所示。

图3 C05 卫星钟差序列及其RC1～RC22 序列Fig.3 Time-series of the C05 satellite clock offset and the first twenty-two RCs

从图3 可以看出，第1 阶RC1序列与原始钟差序列变化趋势非常相似，代表卫星钟差趋势项；RC2+RC3、RC4+RC5+RC6、RC7+RC8+RC9+RC10、RC11+RC12+RC13合成信号表现为显著的周期性变化，代表卫星钟差周期项；第14～22 阶RC14～RC22序列表现为不规则的拟周期性变化，RC22以后的高阶RC序列则表现为随机性变化，限于篇幅未在图3 中给出。

为分析RC 序列的波动特性，利用R/S 分析法计算C05 卫星钟差的各阶RC 序列的Hurst 指数，结果如图4 所示。

图4 C05 卫星钟差各RC 序列的Hurst 指数Fig.4 Hurst exponent of each RC sequence for C05 satellite clock offset

从图4 可以看到，第1 阶RC1序列的Hurst 指数为0.99，与原始钟差的Hurst 指数0.98 几乎相等，说明RC1序列几乎不含波动分量，可直接作为趋势项；第2～13阶RC2～RC13序列的Hurst指数在0.51至0.94之间，且RC4和RC5序列的Hurst 指数大于RC3，这是由于RC4和RC5的周期性特征比RC3更显著，RC2～RC13属于钟差低频信号主导分量；其余RC 序列的Hurst 指数小于0.5，属于钟差高频不规则信号分量。图4 的Hurst 指数计算结果与图2 的RC 序列变化特征具有很好的对应性，说明本文钟差分量辨识方法具有合理性。

为进一步说明本文方法的合理性，利用快速傅立叶变换（Fast Fourier Transform,FFT）对C05 卫星钟差的RC 成分进行频谱分析，分析合并信号的周期性特征，结果如图5 所示。从图5 可以发现，RC2+RC3、RC4+RC5+RC6、RC7+RC8+RC9+RC10、RC11+RC12+RC13合成信号的频谱具有周期性特征，而RC14+RC15+…+RC288合成信号的频谱是非周期性，说明RC14以后的高阶RC 序列属于高频不规则分量，进一步验证了本文钟差分量辨识方法的有效性。

图5 C05 卫星钟差RC 序列的频谱Fig.5 Frequency spectrum of the RC sequence of C05 satellite clock offset

将C05 卫星钟差的RC2～RC13序列相加获得钟差低频周期分量，将第14 阶RC14以后剩余的RC 序列相加获得钟差高频随机分量，其他三颗卫星钟差也作类似处理，结果如图6 所示。从图6 可以看出，C05、C09、C11 和C43 四颗卫星钟差的变化趋势均为递增，表明SSA 方法能准确地分解和重建卫星钟差序列中的趋势项、周期项和随机项，即使近似平稳的微小波动也能检测；卫星钟差呈复杂的多周期性变化，这是由于卫星钟周期性变化不仅与卫星轨道周期相关，还与日食周期、测量性能等多种因素有关，卫星钟差周期项的数值普遍大于随机项，其中，C05 和C09 卫星钟差周期项的数值在ns 量级，C11 和C43 卫星钟差周期项的数值在亚ns 量级，说明GEO 和IGSO 卫星钟受周期项的影响更大；随机项呈不规则性变化，其中，C43 卫星钟随机项的数值小于其他三颗卫星，反映了BDS-3 星载氢钟的优良性能。

图6 四颗卫星钟差趋势项、周期项和随机项的重建结果Fig.6 Reconstruction results of the trend term,period term and random term of the four satellite clock offset

3.2 周期分量特性分析

目前常用多项式拟合残差频谱分析方法对BDS卫星钟的周期特性进行分析，研究表明卫星钟的主周期与轨道周期耦合[18-20]（通常近似为其卫星轨道周期的1 倍或1/2 倍）。为更好地分析BDS 卫星钟的周期特性，利用SSA+FFT 方法对BDS 卫星钟差序列中存在的周期分量进行检测，四颗卫星钟差SSA 重建周期项的频谱分析结果如图7 所示。为进行对比分析，二次多项式拟合残差频谱分析结果也在图7 中给出。

图7 四颗卫星钟差的频谱Fig.7 Frequency spectrum of the four satellite clock offset

通过对比分析卫星钟差重建周期项和拟合残差的频谱分析结果，可以得到如下结论：

1）多项式拟合方法无法完整地提取卫星钟差趋势项，导致频谱图中靠近右侧纵轴处出现残余的趋势信号。拟合残差中包含的剩余随机分量导致靠近左侧纵轴处出现大量不规则波动，而钟差SSA 重建方法提取的周期信号低频处没有出现不规则峰值，高频部分没有出现异常波动，周期分量的频谱特征明显，证明SSA重建方法提取的周期信号不包含趋势分量和随机分量，周期信号提取能力优于多项式拟合方法。

2）C43 卫星钟差拟合残差和重建周期项频谱分析检测的周期存在区别，钟差拟合残差频谱分析检测的周期按幅值排序约为12 h、17 h、21 h、8 h、6 h、4.3 h和3.4 h 等，重建周期项频谱分析检测的周期排序约为12 h、17 h、8 h 和21 h。卫星钟差的周期性变化主要由卫星运行所致，拟合残差的高频周期约为主周期12 h 的公约数，这些高频周期分量可能是在随机分量的干扰下耦合产生的，SSA 信号重建方法能有效地分离和提取钟差周期分量和随机分量，避免了这种耦合现象。

3）卫星钟周期与其卫星轨道周期密切相关，C05卫星钟周期为24 h 和12 h，约为GEO 卫星轨道周期24 h 的1 倍和1/2 倍；C09 卫星钟周期为24 h、16 h和12 h，约为IGSO 卫星轨道周期24h 的1 倍、2/3 倍和1/2 倍；C11 和C43 卫星钟周期为24 h、17 h、12 h和8 h，近似为MEO 卫星轨道周期12 h 的2 倍、4/3倍、1 倍和2/3 倍。其中，C11 卫星钟差拟合残差频谱分析未检测到8 h 周期信号，这是由于钟差拟合残差中的高频不规则信号干扰了频谱分析结果，周期分量被干扰信号淹没，重建周期项频谱中部分周期信号的波峰附近也存在多余信号，这可能由多周期频谱旁瓣引起，但相比钟差拟合残差频谱有明显改善；MEO 卫星钟主周期分量的幅值比其他轨道类型的卫星钟低一个数量级，说明MEO 卫星钟受周期项的影响较小。

4）本文方法与多项式拟合方法对于C05 和C09卫星钟差的周期项分析结果相同，但对于C11 和C43卫星钟差的分析结果却存在差异。具体而言，C11 和C43 卫星钟差的多项式拟合残差频谱不仅受高频信号干扰，而且存在残余趋势信号，这是由于多项式拟合方法仅能提取钟差序列中的“平均”趋势性效应，无法反映钟差信号的时变特性，对于变化复杂的C11 和C43 卫星钟差（图6）多项式拟合残余的趋势分量和高频分量势必影响钟差的频谱分析结果。

3.3 频率稳定度分析

为进一步检验本文方法分离和提取卫星钟差周期项的准确性，分别计算四颗卫星钟差扣除周期项前后的频率稳定度，定量分析卫星钟周期项对频率稳定度的影响。为消除铷钟和氢钟频漂对频率稳定度分析的影响，同时提高方差估计值的置信度，选择能较好消除频漂影响的重叠哈达玛方差计算卫星钟的频率稳定度[14]，同时对比本文方法和多项式拟合方法的周期性波动校正效果，结果如图8 所示。从图8 可以看出，由四颗卫星的原始钟差计算的重叠哈达玛方差曲线均存在不同程度的隆起“鼓包”现象，其中C05 和C09卫星钟差的重叠哈达玛方差曲线隆起现象更为明显。根据式(2)易知，这种隆起异常程度与周期项大小密切相关，即周期项越大，隆起越明显。对比图7 中的四颗卫星钟差的频谱不难发现，C05 和C09 卫星钟差的周期项幅值大于C11 和C43 的周期项幅值；重叠哈达玛方差曲线“鼓包”与重建周期项的重叠哈达玛方差曲线的凸起比较吻合，表明SSA 重建钟差周期项能表征周期项对频率稳定度的影响；在扣除钟差周期项后，各卫星钟在不同取样时间内频率稳定度均有一定程度的提高。从卫星轨道类型而言，C05 GEO 和C09 IGSO卫星钟频率稳定度受周期项的影响最为明显，峰值约为7.5×10-14，C11 MEO 和C43 MEO 卫星钟频率稳定度受周期项的影响较小，峰值分别约为6×10-14和3.5×10-14；从卫星钟类型而言，BDS-3 星载氢钟频率稳定度受周期项的影响最小，这与卫星钟周期项幅值大小密切相关（图7）。多项式拟合校正方法对C05和C09 卫星钟差的频率稳定度提高效果较好，但对C11 和C43 卫星钟差的提高效果有限，结合图7 可以发现，这是由于C05 和C09 卫星钟差的多项式拟合残差的频谱受残余趋势信号的影响较小，而C11 和C41卫星钟差的多项式拟合残差中则包含较强的残余趋势信号；相对于多项式拟合校正方法，无论何种轨道类型和何种星载钟类型，本文方法扣除周期项后的重叠哈达玛方差曲线均有效地削去了原始钟差稳定度曲线的隆起“鼓包”，频率稳定度均得到提高。

图8 四颗卫星钟差的重叠哈达玛方差Fig.8 Overlapping Hadamard variance of the four satellite clock offset

四颗卫星钟差扣除周期项前后的频率稳定度数值结果如表1 所示。由表1 可以看到，利用多项式拟合方法扣除周期项后，C05、C09、C11 和C43 卫星钟的万秒稳分别提高6.5%、6.4%、0.7%和2.8%，C11 和C43 卫星钟日稳分别提高7.7%和1.9%，扣除周期项后C05 和C09 卫星钟的日稳没有得到改善，平均而言，万秒稳和日稳平均分别提升4.1%和2.4%；而利用本文方法扣除周期项后，四颗卫星钟的万秒稳和日稳都得到一定的改善，且相对于多项式拟合方法对卫星钟的万秒稳和日稳的提高效果更为明显，C05、C09、C11和C43 卫星钟的万秒稳分别提高21.0%、23.1%、17.7%和21.6%，日稳分别提高48.8%、54.7%、20.0%和13.1%，万秒稳和日稳平均分别提升20.9%和34.1%，其中，C09 IGSO 卫星钟的提高幅度最大，万秒稳和日稳分别由4.98×10-14和2.87×10-14提高至3.83×10-14和1.30×10-14，MEO 卫星钟的日稳提高幅度较小，这与MEO 卫星钟的主周期（12 h）有关，但由图8 可知，MEO 卫星钟频率稳定度提高最明显的取样时间在3～12 h 之间，对于以MEO 卫星为主的BDS-3 系统，这种稳定度提高非常有利于星座时间基准的建立和保持。

表1 四颗卫星钟的万秒稳和日稳Tab.1 Ten thousand second stability and daily stability of four satellite clocks

4 结论

本文分析了卫星钟周期项对频率稳定度的影响，阐述了周期项分离和提取的必要性，提出了融合SSA分解和R/S 分析方法的卫星钟周期项分离和提取方法。利用不同类型的BDS 卫星钟差数据进行分析和研究，结果表明，所提方法能有效地分离卫星钟差序列钟的趋势项、周期项和随机项，所提取的周期项频谱图比钟差拟合残差频谱图更加清晰，周期性特征更为明显。通过对比分析扣除周期项前后卫星钟频率稳定度的差异发现，原始钟差的频率稳定度曲线存在隆起“鼓包”现象，且与周期项稳定度曲线的异常隆起吻合度很高，印证了所提方法的有效性；利用所提方法扣除周期项后，卫星钟的频率稳定度提高显著，万秒稳和日稳平均分别提升20.9%和34.1%，这非常有利于高稳定度星座时间基准的建立和保持。本文方法不仅可以应用于卫星钟周期项提取，而且在原子钟信号降噪、钟差预报和时间尺度算法等方面具有潜在的应用价值。