BDS卫星钟差短期预报的LSTM算法

吉长东，朱锦帅,，黎 虎，张 萌，吕广涵

（1. 辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000；2. 张家界市测绘院，湖南 张家界 427000）

0 引言

在多导航系统并存的环境下，实时卫星钟差产品的精度直接影响到全球卫星导航系统（global navigation satellite system, GNSS）的服务能力[1-2]。我国北斗卫星导航系统（BeiDou navigation satellite system, BDS）搭载的原子钟均为自主研发，其精度和稳定性与美国全球定位系统（global positioning system, GPS）和欧盟伽利略卫星导航系统（Galileo navigation satellite system, Galileo）的卫星钟不同，因此，BDS钟差产品的精度和稳定性应受到更多的关注[3]。

超快速钟差产品作为国际GNSS服务组织（International GNSS Service, IGS）的核心产品之一，是实现分米级实时定位服务的重要一环。文献[4-5]通过二次多项式加1个周期项拟合，预报IGS超快速产品中的GPS钟差，但其精度与广播星历相比并没有显著改善。文献[6]把小波神经网络（wavelet neural network, WNN）模型应用到钟差短期预报中，将小波分析与神经网络结合，预报精度优于二项式模型和灰色模型。文献[7]以WNN模型为基础，对数据预处理方法进行改进后，以IGS发布的GPS超快速钟差进行实验并得到优于超快速钟差预报部分数据的结果。文献[8]使用径向基函数（radial basis function, RBF）神经网络，对GPS钟差进行超短期5 min和1 h以及短期1 d预报，均可取得高精度预报结果。文献[9]提出的循环神经网络（recurrent neural network, RNN）是1种具有记忆功能的深度学习模型。文献[10]通过修改隐藏层的结构，建立1种基于RNN的长短时记忆（long short-term memory, LSTM）模型，进而削弱了RNN在长序列训练过程中的梯度消失和梯度爆炸问题。

目前，只有德国地学研究中心（Deutsches Geo Forschungs Zentrum, GFZ）和我国的全球连续监测评估系统（international GNSS monitoring and assessment system, iGMAS）提供开放式BDS钟差产品，且现有的钟差预报算法主要集中于稳定的GPS钟差数据以及IGS发布的钟差产品，对于BDS以及iGMAS发布的钟差产品缺少研究。为此，针对BDS超快速钟差，提出1种将修正1次差的钟差预处理方法和LSTM模型相组合的钟差短期预报方法。

1 钟差数据预处理及LSTM预报原理

iGMAS综合钟差产品的时标采用北斗时（BDS time, BDT），用户将iGMAS的BDT产品转换到GPS时（GPS time, GPST）时，不考虑微小量可以时[11]，其计算方法为

式中：涉及物理量的单位均为 s；tGPST为t时刻GPS钟差值；tBDT为t时刻BDS钟差值。iGMAS提供的钟差产品具体数据如表1所示，钟差产品精度用均方根值（root mean square，RMS）表示。

表1 iGMAS钟差产品介绍[11]

1.1 钟差1次差预报方法原理

为了获得高精确度的时间信息值，原始卫星钟差相位数据具有更多的有效位数。文献[6]提出了对钟差相位数据相邻历元间求1次差的方法，这不仅可以降低原始钟差数据序列中趋势项的影响，进而得到1组有效数字位数减少且数值较大的1组数据序列。设L=｛li,i=1,2,…,n｝（li为i时刻钟差值）为任意1组钟差序列，对序列中2个相邻历元的钟差做差，可以得到钟差的1次差分序列ΔL=｛Δli,i=2,3,…,n｝，其中Δli=li-li-1。利用1次差分后的序列ΔL进行学习、建模并预报第n个历元以后m（m＞0）个历元的钟差的1次差分序列ΔLm=｛Δlj,j=n+1,n+2,…,n+m｝，最后根据预报的1次差序列求累加和与第n个历元的钟差值ln相加，即可得第k个历元的钟差相位值[6]为

式中：Δlk为预报的第k个历元的1次差值，k=n+1,n+2,…,n+m。

1.2 基于钟差1次差分序列的粗差处理

针对钟差相位数据相邻历元间求1次差后的数据序列，文献[12]设计了1种修正1次差法数据预处理方法，以改进中位数为基础的抗差估计探测粗差方法：将每个钟差1次差分数据Δli与1次差分序列的中数及绝对中位数（median absolute deviation，MAD）数倍之和进行比较，MAD可表示为

式中：kΔ为1次差序列的中间数；Median为对序列取中位数,即kΔ=Median｛ΔL｝。若钟差1次差分值满足式（4），则将其判定为异常值并剔除该值，然后通过内插或者置零的方法补充该点数据[13]，即

式中：a为根据实际情况需要确定的1个正整数。

本文遵循探测出的异常值个数不超过建模数据的5%的原则，将a设置为3，这样既在一定程度上降低异常值对模型预报性能的影响，又避免了有效信息的剔除。对于剔除数据历元，这里使用线性插值的方法补全历元。

1.3 循环神经网络原理

循环神经网络是1种可以有效提取、利用和处理时间序列的高位非线性动力学模型[14]。在循环神经网络中，神经元不但可以接受其他神经元的信息，也可以接受自身的信息，将神经元最近1次（或几次）的结果导入自己接下来的运算中，形成具有环路的网络结构。设x=(x1,x2,xt)为模型输入，y=(y1,y2,yt)为模型输出，在1个完全连接的RNN模型中，在时刻t时，令xt表示网络的输入，yt表示网络输出，网络模型如图1所示。

图1 RNN模型及其展开示意图

图1中：h为隐藏层状态（即隐藏层神经元的活性值，也称之状态或隐状态），ht-1、ht-1为t-1、t时刻活性值；U、W和V均为网络参数，U为状态—状态的权重矩阵，W为状态—输入的权重矩阵，V为状态—输出的权重矩阵。则RNN模型在时刻t的参数更新公式为：

式中：zt为t时刻隐藏层的净输入；b表示偏置向量；f(zt)表示非线性的激活函数。从图1和式（5）、式（6）可以看出，t时刻隐藏层的状态ht不仅有与时刻t的输入xt相关，也和上一时刻隐藏层的状态ht-1有关。

1.4 LSTM原理

在RNN实际应用中，隐藏层存储的时间是有限的，同时梯度消失的问题也是长程时间序列预测中存在的主要问题。基于此，文献[15]以LSTM模型为基础，在网络中引入1个新的内部状态ct和门控机制，专门进行线性的循环信息传递，同时（非线性的）输出信息给隐藏层的外部状态th。在t时刻，内部状态tc的公式计算为：

式中fsig为sigmoid型激活函数。

2 钟差预报的LSTM模型构造

利用LSTM模型建立钟差预报模型的方式包括钟差数据的训练和预报两部分。在训练阶段，将钟差序列作为神经网络的训练数据，输入的训练数据既作为网络的输入也作为网络的期望输出，通过LSTM神经元的训练，使其掌握数据前后间的映射关系；在预报阶段，将当前时刻的钟差数据作为网络的输入，通过已建立的网络之间前后映射及拓扑关系，依次进行外推得到预报结果。

根据训练输入的钟差序列有限个样本点的数据特征，参考LSTM预报模型的从简设计原则，构建基于LSTM的深度循环神经网络。对于钟差序列预报的整体框架如图2所示。

LSTM模型预报部分由输入层、隐藏层、输出层、网络训练和参数优化5个部分组成。输入层负责对输入的钟差时间序列进行处理，使得输入的钟差时间序列满足LSTM神经元的输入要求；隐藏层使用基于LSTM神经元的深度循环神经网络；网络训练部分使用自适应优化算法中的亚当（Adam）算法结合LSTM神经元对输入数据进行训练；使用网格搜索和交叉验证算法进行超参数的优化；输出层负责预报结果的输出。

图2 LSTM预报模型流程

以预处理后钟差时间序列L0= {l1,l2,…,lr}为例，将其导入网络进行训练、预报计算的具体流程如下：

1）数据标准化。将预处理后的钟差时间序列L0= {l1,l2,…,lr}分 为 训 练集Ltr= {l1,l2,…,ls}和 测试集Lte= {l s+1,ls+2,…,lr}，其中1＜s＜r。之后，对训练集中的每个数据lu(1 ≤u≤s)通过式（14）进行Z-score标准化，得到标注化数据集Lt′r= {l1′,l2′ ,…,ls′}，即

式中：μ为序列Ltr的均值；σ为序列Ltr的标准差。

2）数据分割。根据数据分割长度L，将标准化后的训练数据集Lt′r分割为长度为L的(s-L+1)个数列，使训练集输入变为X= {X1,X2,…,Xs-L+1}，其 中Xp= {l′p,l′p+1,…,l′p+L-1}；则此刻 对应 的 理论 输出Y= {Y1,Y2,…,Ys-L+1}，其中Yp={l′p+1,l′p+1+2,…,l′p+L}，1≤p≤s-L+1。

3）数据训练。当模型的输入为X时，令其经隐 藏 层 计 算 后 的 输 出 为P= {P1,P2,…,Ps-L+1}，则Pp=fLSTM-forward(X p,c p-1,hp-1)，Xp表示当前时 刻 的输入，cp-1表示前一历元LSTM神经元的状态值，hp-1表示前一历元LSTM神经元的输出值，fLSTM-forward表示LSTM神经元的更新公式。网络训练过程中以求得损失函数结果更小为模型优化目标，在这个过程中使用小批量梯度下降法（minibatch gradient descent, M-BGD），每次选择批尺寸大小为B个输入值作为1批输入网络进行训练，通过Adam优化器使得网络不断更新各个参数的权重，最后得到1个对输入序列X拟合效果最好的LSTM模型，记作fLSTM-net。

4）超参数优选。将网格搜索和交叉验证方法结合，对数据分割长度L和批尺寸大小B进行优选。

5）数据预测。预测是采用逐个历元数据迭代的方法，利用训练好的fLSTM-net模型，将训练集的最后1组 数 据输 入fLSTM-net中，得 到 的 输 出 结 果 为Ps-L+1=fLSTM-net(Ys-L+1)={ps-L+2,ps-L+3,…,ps+1}。将s+1时刻的预测值ps+1与YL合 并 为1组新 的 数据将YL+1带入fLSTM-net中，可以得到s+2时刻的预测值ps+2；按照这个方法依次计算，得到预测序列P0={p s+1,ps+2,…,pr}。将预测序列P0进行Z-score反标准化还原，就可以得到预报结果。

3 实验与结果分析

实验数据选取iGMAS超快速（iGMAS ultrarapid, ISU）钟差文件，其中包括超快速观测部分（ISU observed, ISU-O）和超快速预报部分（ISU predicted, ISU-P），数据文件格式为isu*****_00.clk，*****为BDT，前4位是BDT星期，第5位是星期内天）。考虑BDS包含地球静止轨道（geostationary Earth orbit, GEO）、倾斜地球同步轨道（inclined geosynchronous orbits, IGSO）及中圆地球轨道（medium Earth orbit, MEO）三种不同轨道的卫星，不同轨道所处空间环境不同，因此在每个轨道随机选取某一刻卫星进行实验，本文选取的是C04、C07和C12三颗卫星进行实验。最后，以实测部分ISU-O为训练数据和期望输出，采用均方根误差（root mean squared error，RMSE）和极差（Range）2个性能指标来评估预报结果的精度与稳定性[9]。RMSE和R的具体定义为：

式中：R为极差；εi为i时刻期望输出li与预报值lˆi的残差；εmax和εmin分别代表差值序列中的最大值与最小值。

3.1 超参数优选实验分析

结合网格搜索和交叉验证的方法，对LSTM神经网络中使用的超参数窗口分割长度L和批尺寸大小B两种超参数进行优化。根据数据类型设计超参数的取值范围，其中窗口分割长度分别取0.5、1.0、1.5、2.0、3.0、4.0、6.0及12 h内的数据个数，即L∈{2, 4, 6, …, 24, 48}，批尺寸大小B∈{2, 4,… ,48}，以训练集中训练数据的误差函数最小为目标进行训练，运用多层网格搜索对超参数L、B进行优选。

以C12卫星的ISU钟差预报为例，选取2019-01-14的15 min采样率共96个历元的ISU-O数据进行训练，预报96个历元数据。以2019-01-15的ISU-O数据为真值，计算RMSE值，遍历本文所选的L、B两个超参数的取值范围，绘制超参数优选图如图3所示。

图3 参数优选

从图3中可以看出，当批尺寸大小B的取值在20以内时，LSTM网络的预报结果不甚理想，且容易陷入局部最优解，B的取值大于20时，不同的窗口分割长度预报结果的RMSE值变化在0～1.2×10-8的范围内。表2给出了C04、C07和C12三颗卫星使用上述预报方案时，前5组最优参数组合以及预测精度。综合来看，当L=4、B=30时，预报精度更为理想。

表2 C04、C07、C12卫星LSTM模型前5组最优参数组合及预测精度

3.2 不同数据预处理方法的LSTM预报实验

不同数据预处理方法对神经网络的预报精度有着重要影响，使用文中提到的2种1次差数据预处理方法与未使用1次差处理的数据进行LSTM模型预报精度分析。为了便于区分，将未使用1次差处理钟差数据的LSTM预报方案命名为LSTM_1，使用钟差1次差法处理的预报方案命名为LSTM_2，经过改进MAD法粗差处理后的1次差数据预报方案命名为LSTM_3。为了对比这3种方案的预报性能，以日期2019-01-09的ISU-O数据进行预处理，然后训练预报第2天（2019-01-10）的钟差，以第2天的ISU-O数据为真值与预报结果作对比，其部分预报结果的相位数据及1 d的预报残差如图4所示。

将图4中实验结果与对应的真值进行对比后可以看出，LSTM_2方案的预报残差变化范围最大，其前5个历元的预报结果最好，经过1次差处理后，钟差数据中的粗差被放大，导入网络训练的输入数据中粗差对网络训练的影响增大，干扰了模型的训练、拟合过程，造成网络出现过拟合现象，使网络结构极不稳定，因此当预报数据长度超过20个历元后，其预报残差将持续增大（包括向正无穷与负无穷）。对比LSTM_1和LSTM_3这2种实验方案：前5个历元，2种方法的残差大小和变化趋势有明显差距，在前20个历元内，2种方法的预报残差趋于稳定，变化趋势相近；但是LSTM_3的残差波动幅度明显较LSTM_1的残差浮动小。

图4 不同预处理方法后预报结果与残差图

为进一步对比3种预处理方法对预报结果精度的影响，将预报结果与ISU-O真值进行RMSE和R值统计，如表3所示。

从表3可以看出：LSTM_2由于原始数据中粗差被放大造成网络训练过程中不稳定，其R值从6 ns到30 ns，而RMSE值高达16 ns，预报结果的精度和稳定度较差；对比LSTM_1和LSTM_3，这2个模型的RMSE和R值均在纳秒级，其中以LSTM_3方案的结果最好，其RMSE和R值都能保持在2 ns以内。通过对1次差序列进行改进MAD法粗差探测、剔除与补齐，不仅削弱了钟差相位数据中趋势项和粗差数据对网络训练的影响，还能减少输入数据有效数位长度，使输入的数据更能表征数据的变化特点，减少网络的冗余。

表3 不同数据预处理方法预报精度统计表 单位：ns

3.3 ISU钟差预报精度分析

图5 ISU-O短期预报结果与残差

为了进一步验证LSTM在BDS的ISU钟差产品中的预报性能，选用日期2019-01-09的ISU-O数据作为LSTM和WNN模型的训练数据，以ISU钟差产品中的预报部分和WNN预报结果作为对照。选择3.2实验中LSTM_3使用的数据预处理方法对实验数据进行预处理，然后将预处理后的数据，分别通过WNN和LSTM两种神经网络模型进行训练并预报第2天（2019-01-10）的钟差，以第2天的ISU-O数据为真值，对比ISU-P数据以及WNN和LSTM预测结果，其部分预报结果的相位数据及1 d的预报残差如图5所示。

从图5每颗卫星前12个历元的预报钟差值可以发现：LSTM、WNN和ISU-P的结果与ISU-O的数据趋势基本相同，结合预报残差图分析，ISU-P数据的残差偏离程度最大，尤其是C04卫星，其ISU-P数据的残差呈递增趋势，其预报误差达到了30 ns以上；WNN的预报结果具有明显的趋势项，从C04和C07 2颗卫星的预报结果来看，其预报误差呈明显的递增或递减趋势，但是其离散程度与ISU-P的残差相比明显较小；而LSTM的预报结果较前二者明显更好一些，表现了明显的“自我约束能力”。

为了反映LSTM、WNN和ISU-P短期预报的性能，将3种预报结果前1、2、3、4、6及12 h的预报结果精度进行RMSE统计，如表4所示。

表4 LSTM、WNN和ISU-P短期预报RMSE统计表

从表4可以看出：3种预报方案在最初1～2 h的预报精度波动较为明显，随着预报时间增加，预报精度逐步提升。为了整体了解LSTM、WNN和ISU-P 3种结果的预报精度和稳定度，将3颗卫星的3种结果与ISU-O作比较表征其精度与稳定度的RMSE与R值进行统计如表5所示。

表5 LSTM、WNN和ISU-P 1d预报精度统计表 单位：ns

从表5中可以看出，LSTM模型的预报精度和稳定性在3颗卫星中表现得最好，3颗卫星平均精度在1 ns左右，稳定度在2 ns左右，尤其是C12卫星的预报精度可达0.85 ns，稳定度在1 ns左右。而ISU-P的预报精度和稳定性最差，1 d内误差最大可达33 ns，同时精度也达到18 ns，这样的精度在进行单点定位时误差在米级。

对比这3种预报结果，WNN和LSTM的日预报精度均优于ISU-P数据，其中在C04、C07和C12三颗卫星中，LSTM的预报精度较ISU-P数据相比精度分别提高了93%、59%和70%，较WNN分别提高了45%、50%和38%。相较于WNN模型，LSTM更适合作为1种较有效的钟差预报方法应用到钟差短期预报中。

4 结束语

为了提高钟差短期预报性能，提出1种1次差和LSTM模型组合的钟差预报模型。经过超参数优选及2种数据预处理方法对比，搭建适合BDS超快速钟差预报的组合LSTM模型，将LSTM预报结果分别与WNN模型预报结果及ISU-P数据比较，得到如下结论：

1）采用网格搜索和交叉验证结合的方法对L和B共2个超参数进行优选，提高了LSTM模型的预报性能和泛化能力。当L=4、B=30时，LSTM模型预报精度较好。

2）对钟差1次差分数据进行粗差剔除，不仅可以消除钟差数据中的异常点，还减少了有效数据长度，可降低LSTM网络的复杂程度。通过这种预处理方法后，LSTM模型预报结果较直接使用LSTM模型进行预报的精度提高了76%。

3）钟差数据1次差剔除粗差预处理后，再使用LSTM模型预报钟差，其单日预报精度较WNN预报模型提高了45%，与ISU-P数据相比，则提高了74%。该模型至少在日内预报精度优于ISU-P数据，可以作为1种较为有效的钟差预报模型。