基于BDS-2/BDS-3联合处理的北斗超快速钟差预报优化策略

2020-09-23马昌忠王潜心闵扬海王泽杰

导航定位与授时 2020年5期

马昌忠，王潜心，胡超，闵扬海，王泽杰

(1. 中国矿业大学自然资源部国土环境与灾害监测重点实验室，徐州 221116;2. 中国矿业大学环境与测绘学院，徐州 221116)

0 引言

目前，国际GNSS监测与评估系统(interna-tional GNSS Monitoring and Assessment System, iGMAS)发布的超快速钟差产品中，中国北斗卫星导航系统(BeiDou Navigation Satellite System, BDS)观测和预报部分精度分别约为0.60ns和6.00ns，远远不能满足实时cm级应用需求[1]。因此，在提高全球导航卫星系统产品可用性、稳定性和可靠性的前提下，有必要进一步完善BDS超快速卫星钟差产品质量，尤其是预报钟差精度。

当前，针对卫星钟差预报的研究主要集中在3个方面：1)卫星钟差序列预处理策略[2-4]；2)预报模型精化处理[5-7]；3)环境因素对钟差建模的影响分析[8-9]。此外，在钟差预报算法和策略研究中，学者主要关注钟差序列长期预报模型研究，如扩展状态模型[10]和人工神经网络预报[11]；多个短周期项叠加预报，如改进的迭代法以及基于恒星日滤波的单天内钟差变化量预报等[5]。通常卫星钟差预报模型采用趋势项(多项式)与周期项表示[5,12]。同时，毛亚等提出了一种改进的预报策略[13]，即考虑星间相关性对模型参数估计的影响，一定程度上提升了BDS-2钟差预报精度。

然而，BDS超快速钟差建模中广泛采用组合趋势项与周期项的函数模型[2,14]，主要以全球定位系统(Global Positioning System，GPS)相关研究为参考。此外，黄观文等提出了一种改进的BDS超快速钟差预报策略，其中非线性项通过BP神经网络(Back Propagation Neural Network，BPNN)进行建模处理[1]。陈金平等评估了基于星间链路观测数据的BDS-3试验星钟差预报特性[15]。此外，研究表明BDS-3配备的星载原子钟频率稳定性较BDS-2提升了1个量级以上。胡超等通过引入星间相关性[16]，联合预报了BDS-2和BDS-3超快速产品，并评估了星间相关性对超快速钟差预报产品精度的影响[17]；基于1个月的实验数据表明，引入星间相关性可实现BDS-2和BDS-3卫星18h预报钟差精度分别提高30.7%～47.3%和49.9%～59.3%。因此，为进一步提高北斗超快速预报钟差产品质量，后续研究中有必要对BDS-2和BDS-3钟差联合处理中的相关性系数进行深入挖掘分析与处理。

与当前广泛使用的钟差预报模型[1,13]相比，北斗卫星钟差建模中，模型参数估计存在2个主要问题需要优化：1)建模过程中，通常将趋势项与周期项分为两部分进行独立处理，不可避免地忽略了未知参数之间的相关性[16]；2)在建模中，只考虑显著周期项(如BDS-2 2个周期、BDS-3 3个周期)[17]，将导致模型残差无法建模。因此，为了获得更准确且可靠的北斗钟差预报模型，需要对BDS钟差建模策略进行优化处理，如精化钟差模型项以及模型系数等。本文将利用一步求解策略，实现北斗各卫星钟差模型的统一估计；同时在一步模型参数估计中，对预报模型的精度、参数数目和参数类型(周期项和趋势项个数)进行权衡；为获得最优钟差模型，模型参数估计中借助机器学习中稀疏建模的思想进行优化处理[18]。

本文主要针对传统钟差预报模型建模策略的缺点，采用BDS-2/BDS-3钟差序列联合处理的方法，优化北斗超快速卫星钟差产品。首先，基于稀疏建模方法提出了一种模型选择方法，对传统的钟差预报两步策略进行改进；其次，设计了一种BDS-2/BDS-3联合处理策略，以充分利用BDS-2与BDS-3星间相关性实现参数解算增强；然后，考虑预报模型残差序列的时空相关性，利用半变异函数构建经验模型，精化北斗钟差模型参数；最后，通过大量实验对本文的优化策略进行验证。

1 北斗超快速钟差预报优化策略

1.1 BDS-2/BDS-3联合处理模型

在对北斗原始钟差序列进行严格预处理的基础上，如粗差探测与修复和降噪等[16]，可将超快速解算钟差序列设为L，则钟差模型可表示为[1,13,16]

(1)

式中，ti和k分别表示第ti个历元时间和第k颗卫星；a0、a1和a2为趋势项(多项式)系数；n为周期项总数，j为第j个周期；Aj、Bj和Tj表示周期项振幅及其相应周期大小；ε(ti)表示模型残差项。

式(1)中，周期项通常采用快速傅立叶变换算法计算得到。因此，根据得到的周期项，相应的钟差模型误差方程可表示为

(2)

式中，Lk为钟差减去周期项结果。需要注意的是，传统钟差序列建模策略中[16-17]，解算钟差模型参数(趋势项与周期项)分为2个独立处理过程分别进行平差估计；理论上单个模型的所有参数需经统一解算得到。同时，由于星载原子钟的特性存在差异，所有钟差模型参数的个数和类型不可能表示为统一形式；因此，为获得更精确的钟差模型，有必要细化北斗钟差模型参数。

首先，一步估计式(2)中所有模型参数。为提高计算效率，同时解算所有卫星多项式参数，假设有b颗卫星、s个历元，式(1)矩阵形式为

(3)

(4)

式(4)中，权矩阵P(t1)可表示为

(5)

rkb=

(6)

在钟差模型参数估计中，有必要对权阵进行调整，以提高预报模型的精度，本文在权阵中引入时间变量，即

(7)

式中，Δt表示相邻历元的时间间隔。式(2)～式(7)，即为基于BDS-2/BDS-3星间相关性的北斗钟差预报模型一步建模策略。

1.2 基于稀疏建模的钟差模型估计

在1.1节的基础上，钟差模型参数估计考虑了卫星间相关性，间接地优化了模型参数求解过程。然而，参数估计中仍存在2个关键问题会导致待估参数无法获得最优解：1)由于建模中观测值(钟差序列)有限和模型过度拟合等因素影响，模型参数的稳定性和精度不可避免地受到限制；2)不同类型星载原子钟呈现出不同的特性(钟差模型)[19]，无法基于统一的模型进行描述(相同的待估参数)。为提高钟差建模精度与稳定性，有必要进一步优化北斗卫星钟差建模策略。

本节利用机器学习中稀疏建模的思想进行钟差模型精确构建。理论上，引入所有的自变量进行建模处理可有效提高建模精度，但同时也增加了模型复杂性并导致过度拟合。在实际应用中，为了提高模型的可靠性与预报精度，必须考虑模型自变量的合理选取，即模型选择(或变量选择)[20-21]。通常，模型选择只能在有限的备选模型范围内进行，而稀疏建模可以有效地选择建模数据与自变量(模型项)，从而实现模型预报精度和计算效率的提升。为对北斗钟差序列进行稀疏建模处理，本文采用Tibshirani提出的最小绝对值收敛和选择算子(Least absolute shrinkage and selection operator，Lasso)算法[22]。在Lasso算法中，通过最小化多变量线性回归模型的损失函数，对回归系数绝对值施加约束，实现了参数估计和模型选择同步进行。根据Lasso算法原理，可实现对部分变量进行压缩和消除；Lasso算法的具体实现是通过L1范数的最小二乘法，即

(8)

(9)

式中，λ为正则化系数。等式右边相当于模型残差平方和与正则化函数的组合；区别于Tikhonov正则化的L2范数，L1正则化范数通过自动模型选择达到稀疏性的同时，不失函数模型的凸性(产生少量特征)。为求解式(9)，可通过迭代加权最小二乘法求解Lasso问题的数值解[23]，即式(9)可表示为

(10)

(11)

式中，U(m)是权阵(对角阵)，其元素根据第m-1次迭代解确定。因此

(12)

通过矩阵求逆公式，式(12)的解为

(13)

1.3 基于半变异函数的模型残差处理

通过上述基于稀疏建模思想的一步求解钟差模型参数，可较好地实现北斗钟差预报模型构建。但是受模型残差以及随机噪声影响，建模过程中需对模型残差序列进行深入分析。前期研究中，针对预报模型残差序列提出了几种改进的处理策略[1,13]；此外，学者组合了偏最小二乘法(Partial Least Square，PLS)和BPNN算法，并将其应用于超快速钟差预报中[14,17]；但缺少卫星钟差精度的相关信息(方差-协方差阵)，残差序列建模及其预报难以达到预期精度[17]。本文从钟差序列特点的角度考虑，即所有卫星建模残差序列不可避免地呈现时空相关性；因此，为充分挖掘BDS-2与BDS-3钟差序列各卫星潜在的时空特性，本节提出了利用半变异函数进行残差序列优化处理。有关半变异函数的更多细节已在相关研究中进行了详细讨论与总结[17,24-25]。

在实际建模处理中，通常使用改进的残差半变异函数模型[19,24]

(14)

2γ(ti,tl)=2(C(0)-C(ti-tl))

(15)

式中，C表示参数协变差运算。C(tl-ti)可以通过式(16)计算得到

C(ti-tl)=Cov(ε(ti),ε(tl))

(16)

其中，Cov表示协方差运算；模型残差的时间相关性将由半变异函数进行数值量化。

在得到每颗卫星量化的残差序列时间相关性之后，可拟合出经验半变异函数模型，如球形和指数型等。基于得到的经验半变异函数模型(如式(17)与式(18))，当h→∞时，2γ(h)→2C(0)；由式(15)可知，协变差C(h)可由经验函数模型确定。因此，在式(16)的基础上可进一步构造参数估计的权阵，并将其代入式(3)中，从而实现模型系数估计的精化处理。

综上，本文从星间相关性、稀疏建模以及半变异函数等3个方面对北斗超快速预报钟差模型进行了精化处理。首先，基于稀疏建模对钟差模型参数进行一步估计，实现了模型的自动选择，并改进了传统钟差分步建模策略；其次，考虑BDS-3星载原子钟具有更强性能，通过提取星间相关性实现了模型参数的解算增强；最后，针对模型残差序列，利用半变异函数对参数估计随机模型进行了精化处理。

2 BDS-2/BDS-3卫星钟差预报实验分析

本节将对改进的BDS-2/BDS-3卫星超快速钟差联合预报策略进行实验分析。实验前，首先基于连续1个月(2019年第41～70天)的BDS-2和BDS-3观测数据进行超快速联合定轨实验，得到了解算的超快速卫星钟差产品。钟差预报实验中，根据本文提出的改进策略，将实验分为3组，分别从星间相关性、稀疏建模和半变异函数模型的角度对钟差预报模型进行验证分析。考虑到北斗地球静止轨道(Geostationary Earth Orbits，GEO)卫星钟差精度较低的特点，在实验中剔除了相应部分。

第1组实验中包括4个对比方案，验证了基于稀疏建模的钟差序列一步建模策略：

方案1：基于趋势项和周期项，建立钟差预报模型，该方案中，趋势项系数通过拟合钟差序列得到，而周期项则通过残差的快速傅里叶变换确定；

方案2：与方案1类似，其中BDS-2卫星周期项参考黄观文等[1]以及毛亚等[13]的研究结论，BDS-3则利用3个主要周期项建立钟差预报模型；

方案3：与方案1和方案2相比，将快速傅里叶变换获得的所有周期项加入钟差序列建模中；

方案4：在方案3的基础上，采用Lasso算法对模型参数进行一步估计，以验证稀疏建模的有效性。

根据上述四种方案分别进行了1天弧长的BDS钟差预报实验，并以武汉大学iGMAS分析中心发布的BDS快速钟差为参考进行精度分析。各方案中BDS-2与BDS-3卫星的钟差建模残差以及相应预报精度如图1和图2所示。

通过上述4套对比实验方案可知：1)BDS-3卫星钟差建模残差明显小于BDS-2卫星，间接说明了BDS-3星载原子钟性能优于BDS-2；2)方案1与方案2的拟合残差近似相等，而方案3中利用所有周期项建模得到的模型残差最小；3)基于图2中钟差预报精度的统计结果，对BDS-2而言，方案3较传统方法(方案1与方案2)预报精度有所提高，且方案4(一步建模)略优于方案3(除C13外)；对于BDS-3卫星，方案3一步估计所得到的预报精度有所降低，但在稀疏建模之后(方案4)，在方案1与方案2的基础上，一定程度实现了预报精度的提高。

图1 年积日41天(2019)不同北斗卫星钟差拟合残差Fig.1 Fitting residuals of BDS clock offsets for different satellites on DOY 41

图2 四种方案连续30天北斗钟差预报精度(18h)平均值Fig.2 Average daily RMS values of BDS predicted clock offsets for different satellites in 30-day experiments with four schemes (18 hours)

第2组实验包括2个对比方案，验证了引入星间相关性求解模型参数的可行性：

方案5：类似方案1，但在模型参数求解中加入了星间相关性系数，该方案具体实验细节可参考作者前期研究成果[16-17]；

方案6：同样地，在方案4的基础上，将获取的星间相关性系数加入卫星钟差模型参数估计中。

图3所示为连续10天C14与C24卫星预报钟差的精度；表1则给出了1个月钟差精度平均值(18h)及其精度提升率。实验结果表明，星间相关性可显著改善北斗卫星预报钟差精度；且与传统建模方法(方案1)相比，BDS-2与BDS-3卫星钟差预报精度分别提升了32.9%和16.9%；而方案4考虑了星间相关性后，BDS-2与BDS-3卫星钟差18h预报精度分别提升了27.2%与28.6%。

图3 不同方案连续10天钟差预报精度Fig.3 Prediction accuracy of clock offsets in 10 consecutive days for different schemes

表1 不同方案北斗卫星钟差预报精度(ns)及其提升率(18h)

通过上述引入星间相关性与稀疏建模处理后，模型残差中仍包含了钟差序列的有效成分，其一定程度上降低了北斗钟差预报的精度。第3组实验主要分析了本文提出的超快速钟差预报策略的正确性，并验证了利用半变异函数精化模型参数估计的随机模型的可行性。在钟差预报实验前，首先基于方案4中的模型残差序列进行时空相关性估计，结果如图4所示，其中虚线表示基于1个月残差序列计算的实验半变异函数值；其次，构建了一个经验模型，即半变异函数模型，实验中选择球面模型进行残差半变异函数模型构建，即如式(17)与(18)所示

γ(h)BDS-2=

(17)

γ(h)BDS-3=

(18)

图4 BDS-2和BDS-3方案4钟差残差序列半变异函数Fig.4 Variogram from the satellite clock offsets residuals of scheme 4 for BDS-2 and BDS-3 satellites

方案7：基于BDS-2与BDS-3联合超快速钟差序列，考虑星间相关性，通过方案1进行钟差预报；

方案8：基于BDS-2与BDS-3联合超快速钟差序列，考虑星间相关性，通过方案2进行钟差预报；

方案9：基于BDS-2与BDS-3联合超快速钟差产品，通过方案6进行钟差预报，每个预报弧段长度为18h，并利用BPNN算法对模型残差进行处理与预报，最后将两部分预报钟差序列合并；

方案10：类似于方案9，针对模型残差部分采用灰色模型进行预报处理；

方案11：基于方案6，采用PLS+BPNN策略对模型残差进行建模与预报处理[24]；

方案12：基于方案6，采用经验半变异函数模型计算的模型残差协方差阵对式(4)中的权阵非对角线元素进行更新处理，同时采用PLS+BPNN策略对模型残差进行建模与预报处理。

为对本文提出的改进的钟差预报模型进行全面分析，将实验中12个预报方案结果与WHU的快速钟差产品进行对比。为讨论不同方案的钟差预报精度，图5中选取了6套方案中C14与C24卫星对比分析建模残差的差异，通过模型残差可明显发现方案11的建模效果最优。

图5 基于2019年第41天的C14和C24不同方案模型残差Fig.5 Model residuals of different schemes based on DOY 41, 2019 for C14 and C24 satellites

同时，为具体说明方案7～12中不同实验的钟差预报精度，表2给出了不同方案下钟差预报精度平均值。由于超快速钟差产品以3h时延和6h间隔进行更新发布[16]，因此，钟差预报结果中仅分析了1天内18h的钟差预报精度。

表2 基于传统方法、星间相关性和Lasso算法的18h超快钟差预报精度(ns)及不同方案提升率

基于图5与表2中各方案的统计结果，可以得出：1)与方案2相比，考虑BDS-2/BDS-3卫星钟差相关性对模型残差的影响不明显；2)基于广泛使用的周期项选取方法进行钟差建模，将导致BDS-2与BDS-3预报钟差精度分别下降4.3%和21.1%；3)通过灰色模型对模型残差进行处理，BDS-2与BDS-3钟差预报精度分别提高了23.3%与16.9%。需要说明的是，方案8是在传统方法的基础上进行建模优化处理，需对其模型残差进行补偿。因此，实验中设置了基于BPNN和PLS+BPNN算法的钟差残差序列预报方案。实验结果表明，PLS+BPNN策略可实现模型残差较传统BPNN算法建模降低0.5%和6.2%。同时，由于模型残差序列中存在时空相关性，在系数估计中引入半变异函数对权阵进行精化处理，方案12较方案11可进一步实现BDS-2与BDS-3预报钟差精度分别提升8.0%与11.1%。