附有周期项的二次多项式LASSO钟差预报模型

李坤，王潜心，闵扬海，龚佑兴，苗伟，程彤

(1.中国矿业大学 环境与测绘学院，江苏 徐州 221116；2.苏州市房地产市场与交易管理中心，江苏 苏州 215002；3.国防科学技术大学 指挥军官基础教育学院，湖南 长沙 410072)

0 引 言

北斗卫星导航定位系统星载原子钟因自身时频特性的复杂性和外部因素影响，导致其系统时间与地面控制系统的原子钟钟面值存在偏差，该偏差即为卫星钟差[1-2]。目前，国际GNSS监测评估系统(international GNSS monitoring &assessment system,IGMAS)最终解算的钟差产品虽精度较高[3]，但解算具有一定滞后性,无法满足实时或近实时用户，因此，开发预报高精度钟差的产品非常重要。

目前，关于钟差预报的模型主要有二次多项式模型(quadratic polynomial,QP)、灰色模型(gray model,GM)、谱分析模型(spectrum analysis,SA)、时间序列模型(autoregressive integrated moving average,ARIMA)、卡尔曼滤波模型(kalman filtering,KF)[1,4]、神经网络模型(back propagation network,BP)等。基于卫星钟的物理特性，国内外学者做了大量关于二次多项式预报模型的研究，并取得了丰硕成果[5-6]。郑作亚等[7]、杨定江等[8]在二次多项式模型基础上，增加了周期项因素，构造附有周期项的二次多项式预报模型；王甫红等[9]在附有周期项的二次多项式模型基础上，提出了一种基于钟差变化率拟合建模的卫星钟差预报方法；MAO Y等[10]在二次多项式模型基础上，采用一种顾及卫星间相关性的钟差预报模型，并对模型中涉及的主要算法进行了推导；HANG G等[11]、艾青松等[12]针对各天体间钟差基准偏差，在二次多项式模型基础上提出了起点偏差修正法；于烨等[13]、王宇谱等[14]考虑了二次多项式预报模型随机性变化的特点，对预报模型进行了进一步优化。常规的二次多项式模型在求取模型系数时，常采用最小二乘估计(least squares estimation,LSQ)算法解算模型系数，但往往导致预报模型过拟合，因此，本文采用LASSO算法进行模型系数的解算。先分别使用LSQ和LASSO算法求得模型参数，比较两种模型拟合残差，分析模型拟合精度，求取模型周期项，然后对附有周期项的各模型进行预报，比较6，12，18，24 h预报精度。

1 二次多项式模型参数求解方法

1.1 LSQ算法

不考虑残差时，可表示为

y=a0+a1t+a2t2，

(1)

式中：a0为钟差；a1为钟速或频偏；a2为频漂；t为时间。

二次多项式预报模型的误差方程可表示为

v=Bx-l，

(2)

式中：v为误差项；B为钟差历元系数矩阵；l为已知钟差；x为待求模型参数。

针对上述预报模型和误差方程，根据最小二乘原理求得模型系数为[15]

x=(BTPB)-1BTPl，

(3)

其中，P为权阵(本文所取得的权阵为单位阵)。

1.2 LASSO算法

LASSO算法是T.Robert[16]提出的一种处理复共线性数据的有偏估计。LASSO算法的代价函数可表示为

J(x)=(Y-Ax)2+λ‖x‖1,

(4)

式中：A为钟差历元系数矩阵；Y为已知钟差；λ为惩罚系数；‖x‖1为x的1范数。

对式(4)右边第一项求偏导,得

-2mj+2xjnj，

(5)

因式(4)第二项不可导，对其使用次导数，

令式(5)，式(6)相加等于0，则

(7)

由式(7)可得

(8)

综合考虑惩罚系数的有效性和算法的运行速率，选取经验系数0.01,0.1,1,2,3,5,10,100,200,500为惩罚系数。为了达到理论最优，惩罚系数应越大，相应附加惩罚系数的参数越小，甚至为0，从而达到减少变量和降低方程维度的目的。相关学者研究了穷举法选取最优惩罚系数，但耗时较长，本文选取常用经验系数，通过模型迭代和筛选，获得最佳惩罚系数，并将惩罚系数回代到模型中进行计算。

2 二次多项式钟差预报策略

2.1 数据预处理

GNSS原始钟差数据存在数据异常问题，同时由于外部因素影响，如设备老化、外部环境变化等导致原始钟差数据中存在粗差和钟跳等情况，对后续钟差预报和分析带来影响[18]。针对上述情况，需要对钟差进行预处理，得到有效数据，为后续数据的正确预报和分析提供可靠基础。

钟差数据预处理方法包括中位数法(median absolute deviation,MAD)和Baarda粗差探测法。中位数法简单、方便，但对粗差的大小不敏感，不能探测出小粗差。Baarda粗差探测法通过对验后残差进行处理，能有效探测出小粗差，但在进行钟差系数拟合求取钟差残差时，常因钟跳等异常使得钟差值无法准确获得拟合系数[19]。

针对Baarda粗差探测法和中位数法的优缺点，对Baarda算法进行优化和改进，先将钟差相位数据转化为频率数据，具体数学模型为

(9)

式中：x(t)，y(t)分别为相位和频率数据；εt和ωt分别为相位和频率数据的时间序列拟合残差。

对频率数据进行最小二乘平差估计得到模型的拟合残差，计算单位权方差，判断最大残差是否超过阈值。实际运算中有部分残差超过阈值，但最小二乘算法具有残差平均功能，不能剔除超限的其他残差，而是重新组建方程，然后再次进行残差的剔除。

最后对剔除残差的钟差频率数据进行二次多项式拟合，补全缺失的数据，同时将频率数据恢复为相位数据，便于后续钟差预报等工作。图1为钟差数据预处理流程。

图1 钟差数据预处理流程Fig.1 clock offset data preprocessing

为检验改进Baarda算法的有效性和可行性，本文采用德国地学中心(GFZ)发布的GBM钟差数据，选取2017年10月1日，即gbm19690.clk钟差数据，采样间隔为30 s，进行钟差数据预处理。如图2所示，为了准确分析和直观显示改进的Baarda算法对数据处理的效果，对GBM钟差数据中的C02号卫星人工添加粗差数据，利用Barrda算法进行粗差探测，并修复钟差数据中的粗差和钟跳数，对于频漂大的钟差数据，很难通过常规方法进行预处理，而通过图2可以看出，改进后的Baarda算法能够很好地识别和处理粗差剔除数据。

图2 C02卫星频率数据预处理效果Fig.2 Frequency data preprocessing result of C02 satellite

2.2 拟合残差分析

自郑作亚等[7]分析了附有周期项的二次多项式预报模型以来，周期项误差已成为钟差基本属性之一，钟差主要分为趋势项和周期项两部分，鉴于此，本文采用附有周期项的二次多项式预报模型，具体模型为

(10)

式中：P为周期项个数；Ai，Bi，t′分别为正、余弦函数的振幅和周期项；εi为附加周期项后的拟合残差。

采用IGMAS 2018年6月18日快速钟差产品进行模型参数求解。首先对钟差数据预处理，获得有效的钟差数据，然后分别采用LSQ和LASSO算法求得附有周期项的二次多项式预报模型系数，用求解的模型参数分别拟合当天钟差，并与当天实际钟差作差，求解模型拟合精度。图3～4为各取两颗GEO卫星(C01，C04)，IGSO卫星(C07，C09)，MEO卫星(C12,C14)进行模型精度分析。

图3 模型拟合残差时间序列Fig.3 Fitting residual time series diagram

由图3可以看出，LSQ算法求得的自拟合残差分布在[-2,2]，且相对集中，而LASSO算法求得的自拟合残差则明显分散，且大于LSQ算法的。可以得出，LSQ算法求取的模型系数自拟合程度明显优于LASSO算法的。

为了定性分析LSQ和LASSO算法的自拟合精度，分别统计图3卫星精度情况，并绘制模型拟合精度直方图(图4)。由图4可以得出，LSQ算法拟合精度整体优于LASSO算法的，除C07卫星外，LSQ算法求解的模型拟合精度较高，小于0.5 ns，可以看出，使用LSQ算法进行模型拟合，求取的模型具有很高的拟合精度。但是否存在过拟合导致模型错误，从而使得模型在预报时偏离太大，下文将对该现象进行具体分析。

图4 模型拟合精度Fig.4 Model fitting accuracy

3 实验结果分析

为验证LASSO算法求得预报模型的可行性，采用IGMAS 2018年6月18日的快速钟差产品进行模型系数求解，以2018年6月19日的快速钟差产品为真值，进行预报精度的分析。

首先对2018年6月18日的钟差数据进行预处理，获取有效频率数据，然后对钟差数据进行附有周期项的二次多项式建模，用LSQ和LASSO算法分别求解模型参数，用获得的模型求解钟差预报产品，对求取的钟差预报产品与6月19日的钟差产品进行对比分析，获取钟差预报精度数据。最后联合拟合残差数据和钟差预报产品精度数据进行过拟合分析，分析LSQ算法是否存在过拟合导致模型错误。实验分析流程见图5。

图5 实验分析流程图Fig.5 Flow chart of experimental analysis

统计分析上述附有周期项的LSQ和LASSO算法的钟差预报模型中预报的钟差，用于后续过拟合现象的具体分析。

根据实际情况，分别统计各卫星6，12，18，24 h钟差预报精度。图6为两种方法求解的模型预报残差时间序列，为更直观地分析两种算法的预报结果，分别统计两种方法的预报均方根误差(root mean square，RMS)，并绘制直方图(图7)。两种算法不同时间段的钟差预报精度数据如表1所示。

表1 两种方法预报精度统计Tab.1 Prediction accuracy statistics of the two methods

由图6可知，LSQ和LASSO算法求解的钟差预报残差为[-40，40]ns，LASSO算法求解的模型预报发散程度明显小于LSQ算法的。由图7可知，LASSO算法求解的预报模型，各卫星的预报精度都优于LSQ算法的。根据图3～4和图6～7对比分析可得，LSQ算法求解的模型拟合精度高，但预报精度差。

图6 钟差预报残差时间序列图Fig.6 Time series diagram of residuals of clock offset prediction

图7 钟差预报精度均方根误差直方图Fig.7 RMS histograms of clock offset prediction accuracy

由于本文采用附有周期项的二次多项式预报模型，周期项数据未做筛选，导致模型变量非常多。常规无偏LSQ算法在处理多变量数据时易产生过拟合现象，即模型本身拟合精度非常高，但将模型外推时其精度很差，因而不能很好地表现模型的一般性。LASSO算法的最大优势为添加了惩罚项，可以降低模型复杂度，即并非所有变量都放入模型中，而是将部分参数的系数变为0，有选择地将一些性能更好的参数加入模型中，从而降低模型维度，避免过拟合现象。

比较图3、图6可以看出，考虑所有周期项变量的LSQ算法自拟合精度非常高，是因为LSQ算法本身对所有周期进行建模，相比LASSO算法剔除一些变量而言，LSQ算法自拟合精度显著高于LASSO算法的，但数据变量庞大和观测数据有限之间的矛盾导致LSQ算法求解预报模型时存在过拟合现象，造成求解的模型预报精度不高，而LASSO算法能很好地避免这一现象发生。

由表1可得，LASSO算法求解的模型预报精度提高较显著，除预报时长6 h的C09，C12卫星和预报时长12 h的C12卫星外，其余的提升效率均优于10%，且随着预报时间增加，LASSO算法求解的模型预报精度提升效果越显著。

4 结 论

(1)LSQ算法求解的模型系数自拟合精度显著优于LASSO算法的，但模型预报精度在预报时长6，12，18，24 h时均低于LASSO算法的，因此，使用LSQ算法进行模型参数求解，存在过拟合现象，从而导致预报精度不高。

(2)LASSO算法求解的模型预报精度较LSQ算法提升明显，除预报时长6 h的C09，C12卫星和预报时长12 h的C12卫星外，其余均提升超过10%。从总体趋势上看，随着预报时间增加，LASSO算法求解的模型优势更为明显，其预报发散程度显著小于LSQ算法解算的模型。因此，LASSO算法可以抑制过拟合现象，大幅提升预报模型的准确性，对钟差预报的改进具有显著意义。