GNSS坐标非线性变化的差分长短时记忆网络预测
贾彦锋1，朱新慧1，叶家彬2，纪秀美3
（1.信息工程大学 地理空间信息学院, 郑州 450001；
2.中建七局交通建设有限公司，郑州 450001；
3.31002部队，北京 100094）
摘 要：为了改进GNSS坐标非线性变化的预测方法，获得更高的预测精度，该文提出了一种基于一阶差分和长短时记忆（LSTM）网络的坐标非线性变化预测方法。首先综合利用奇异值分解和最大互信息系数准则对坐标时间序列进行降噪处理得到真实的非线性变化，接着利用谐波模型对其中的周期项进行提取和预测，剩余的未模型化成分经过一阶差分后采用LSTM网络进行预测，然后将两个预测结果进行综合得到坐标非线性变化的高精度预测结果。实验结果显示，该方法在20 d的预测步长内的平均绝对误差达1 mm以内，相比谐波模型、ARIMA模型和未经一阶差分的LSTM模型的预测方法精度至少提升了78%、25%和22%，具有更高的预测精度。同时经过对比也证明了该方法具有更好的适用性。
0 引言
研究表明，全球卫星导航系统（global navigation satellite system, GNSS）基准站坐标时间序列中包含显著的非线性变化，其振幅可达1~2 cm，这使得基于基准站的线性速度来维持参考框架，其精度仅为厘米级，无法满足开展大范围或全球尺度的毫米级地球动态变化的监测及研究的需求[1]。因此，对坐标的非线性变化进行预测对于地球参考框架的高精度维持具有重要的研究价值和现实意义。
近些年来，对基准站坐标非线性变化的分析和建模主要从两个方向展开。①从引起坐标非线性变化的各种物理因素入手，根据影响机制分别建立模型；②不考虑非线性变化产生的各种物理机制，只根据坐标时间序列本身的特点建立纯数据模型[1-2]。
基于物理机制的方法主要通过研究环境负载（大气负载、海洋负载、水文负载等）和热膨胀效应等对基准站非线性运动的影响构建相应的模型。常见的环境负载模型有全球地球物理流体中心(Global Geophysical Fluids Center, GGFC)建立的GGFC模型[3]、QOCA软件的地表负载位移效应计算模型[4](quasi-observation combination analysis based on loading model, QLM[5])以及精密环境负载模型[6](object model diagrams, OMD)等，常见的热膨胀效应模型有半空间热传导模型[7-8]和统一球体热传导模型[9]。基于物理机制的方法能够很好地避免虚假信号的影响，但是由于引起基准站坐标非线性变化的机制比较复杂，目前常见的环境负载模型仅能解释部分周年、半周年变化，并且模型之间的一致性有待进一步提升[1,6,10-11]。坐标时间序列主要包含趋势项、非线性项和噪声[12]。基于坐标时间序列的研究方法主要是对其中的非线性项进行精准建模。非线性项的建模方法主要分为两种:①时间函数模型，常见的有谐波模型[13-14];②递推模型，常见的为差分整合移动平均自回归 (autoregressive integrated moving average, ARIMA) 模型[15]以及神经网络模型[16-17]。谐波模型等时间函数模型方便建模且能够对任一历元的坐标进行预测，但是由于其只能描述坐标时间序列非线性项中能够模型化的部分，对于难以模型化的成分无法预测，因此其精度有限。ARIMA模型在坐标时间序列预测中能够取得一定的效果，但是其本质上是一种线性模型，只能捕捉线性关系，而不能捕捉非线性关系，而且对数据的稳定性要求较高。以长短时记忆（long short term memory , LSTM）为代表的神经网络模型在时间序列预测方面取得了一定的效果[18]，但是坐标时间序列中非线性变化的波动性和噪声会对细小特征的捕获产生影响，进而影响学习的效果，限制了模型的精度。因此，如何利用坐标时间序列建立高精度的非线性变化预测模型，仍需要进一步的研究。
本文针对LSTM网络在对非线性项中难以模型化成分进行预测时精度受限的问题，借鉴ARIMA模型中的差分平稳化方法，提出了一种基于一阶差分和LSTM网络的坐标非线性预测方法，并利用斯克里普斯轨道和永久阵列中心（Scripps Orbit and Permanent Array Center, SOPAC）提供的时间序列数据进行了实验验证。结果表明，该方法能够对坐标非线性变化进行高精度的预测。
1 基于一阶差分和LSTM的坐标非线性变化预测方法
基于一阶差分和LSTM的坐标非线性变化预测方法的流程如图1所示，该算法主要由数据预处理、一阶差分及其还原和LSTM模型的训练和预测3部分组成。
图 1 基于一阶差分和LSTM的坐标非线性变化预测流程
此处，引入差分方法的目的在于削弱坐标时间序列中突变的影响，获得更加平稳的时间序列用于模型的训练和预测。平稳性刻画的是时间序列的统计性质关于时间平移的不变性，其大小决定了在未来一段时间里时间序列能否按照已有的变化规律延续下去[19]。理论上，经过差分得到的坐标时间序列相较于原始时间序列平稳性能够得到提高，进而能够提升训练得到的模型对时间序列未来变化的预测性能。但是随着差分阶次的增大，坐标时间序列中的噪声会变大，因此本文选择了一阶差分的方法。
原始的坐标时间序列在数据预处理部分经过趋势项的提取分离、粗差剔除、缺失值内插和降噪后得到其中的非线性项。接着利用谐波模型分离非线性变化中的模型化成分和未模型化成分。然后将未模型化成分按照一定的比例划分为训练集和测试集，其经过一阶差分后作为LSTM网络的输入，进行模型的训练和测试。最终谐波模型（预测模型化的非线性变化）和LSTM模型（预测未模型化的非线性变化）共同组成坐标非线性变化高精度预测模型。其中趋势项的提取利用的是QOCA软件，降噪采用的是本文提出的一种基于奇异值分解（singular value decomposition, SVD）和最大互信息系数（maximal information coefficient, MIC）准则的降噪方法，该方法能够在有效降低坐标时间序列中噪声的同时，最大限度地保留原序列中的非线性关系。
1.1 基于SVD和MIC的坐标时间序列降噪方法
坐标时间序列中包含的噪声会影响甚至掩盖真实非线性信号的提取[20]，进而影响非线性特征的建模和预测。因此，需要对坐标时间序列进行降噪处理。SVD方法能够通过分解和重构实现时间序列的降噪，其降噪的效果取决于重构阶次的确定[20-21]，若重构阶次较小会损失真实的信号，若重构阶次过大则降噪的效果不够理想，而目前尚未有科学的重构阶次确定方法。针对这个问题，本文引入MIC准则用于重构阶次的确定，进而实现对坐标时间序列的科学降噪。MIC是一种用来衡量两个变量之间关联程度，即线性或非线性关系的方法，MIC值在0~1之间，值越大表示相关性越强[23]。其相比皮尔逊（Pearson）相关系数、阈值相关、最大相关系数、核密度估计等方法具有普适性、公平性、标准化、计算复杂度低以及鲁棒性强等优点。其在样本量足够大时能够捕获各种相关关系，这种相关关系不限于特定的函数类型（如线性函数、指数函数或周期函数等），同时能够在不同类型强度相近的单噪声影响下给出相近相关系数[24]。
1.2 LSTM模型的训练与测试
首先对训练集和测试集进行一阶差分，降低其波动性。然后进行训练样本和测试样本的划分，以获取训练集和测试集。
在获取样本集后，便可以构建LSTM网络进行模型的训练与测试。本文的网络,搭建和训练是基于深度学习框架keras进行的。其中，输入层的神经元个数等于输入时间序列的长度；输出层的神经元个数等于预测步长；隐藏层的层数、各层神经元个数、学习率、激活函数等超参数均通过网格搜索法进行确定，以得到最优的网络构型。另外，网络采用的是Adam优化算法，并通过引入提前停止（early stopping）方法动态选择迭代次数，在保证网络性能的同时降低运算成本。
2 实验与分析
本文选取的实验数据来自于SOPAC提供的全球IGS14框架下的单天解坐标时间序列（http://sopac.ucsd.edu/）。选用的数据集是联合解算的去除粗差后的坐标时间序列，其经过SOPAC和喷气推进实验室（jet propulsion laboratory, JPL）的联合处理，能够消除一定的计算系统误差，具有更高的精度。
2.1降噪效果分析
在对原始的坐标时间序列进行趋势项的提取分离、粗差剔除和缺失值的补全后，综合利用奇异值分解和MIC准则对其进行降噪处理。下面以ZIMM测站U方向为例分析坐标时间序列的降噪过程和效果。
由于坐标时间序列中包含着明显的周年项，因此SVD的窗口长度设置为365[28]，经过SVD分解后各个重构成分的贡献率如图2所示。图3展示了每一个重构成分的MIC值以及重构序列和对应的残差序列的MIC值。
综合图2和图3可以看出：①随着阶次的增加，重构成分的MIC值整体呈下降趋势，说明阶次越高重构成分与时间之间的相关性越弱；②随着阶次k的增大，重构序列的MIC值在36阶以内时保持为1，而后逐渐下降；③随着阶次k的增大，残差序列的MIC值在3阶以后降到0.2以下。根据本文提出的降噪方法即可作为降噪后的时间序列，见图4。经计算，降噪后的时间序列占原序列的82.67%，而残差序列的MIC值为0.068，说明本文提出的降噪方法能够有效提取出坐标时间序列中非线性成分，证明了该方法的有效性。
图 2 各个重构成分的贡献率
图 3 重构成分、重构序列和残差序列MIC值随阶次的变化
图 4 降噪前后的坐标时间序列
2.2 预测结果与分析
依据图1中的预测算法流程，对坐标的非线性变化进行预测，预测结果如图5~图7所示（以ZIMM测站U方向为例）。本文采用平均绝对误差(mean absolute deviation, MAE)来描述单个测站的预测精度。
由图5可知，随着预测步长的增加，预测的平均绝对误差会逐渐变大，但是在20 d的预测步长内仍可取得1 mm以内的预测精度。图6为在17个步长的预测长度下得到的预测结果，其中训练集和测试集的比例为4:1。另外,通过超参数优化得到的超参数设置为：输入时间序列的长度为80；LSTM网络的层数；每层的神经元个数依次为40和80；学习率；激活函数为tan h；损失函数为MAE。从图中可以看出预测结果与真实结果十分相近。
为了进一步验证方法的预测效果，本文从SOPAC数据集中选取了180个测站进一步的分析，结果如图7所示。考虑到神经网络方法需要大量的数据支撑，奇异值分解需要进行缺失值内插[29]，因此,为了充分反映方法自身的性能，选取标准设置为：时间跨度大于10年，数据缺失率小于20%。从图中可以看出，本文提出的预测算法在U方向的预测精度整体在1 mm以内，180个测站MAE的平均值为0.46 mm。另外,N方向和E方向的MAE平均值均为0.07 mm，进一步验证了本文提出方法的有效性。
图 5 预测结果的MAE随预测步长的变化
图 6 预测时间序列（预测步长为17 d）
图7 U方向预测精度统计（预测步长为17 d）
2.3 方法对比
为了进一步验证方法（一阶差分+LSTM）的性能，本文同时利用谐波模型、未经一阶差分处理的LSTM方法和ARIMA模型对非线性变化进行了预测分析。在评价指标上，采用MAE的平均值描述方法的整体精度水平，采用MAE的标准差SDMAE来表示方法对不同测站进行预测时精度的稳定性。为了更加真实地反映方法的优越性，未经一阶差分处理的LSTM方法和ARIMA模型采用与一阶差分+LSTM方法完全相同的输入（经过降噪处理后分离得到的未模型化成分）。首先，以ZIMM测站U方向为例对比4种方法在不同预测步长下的预测精度，结果如图8所示。
图 8 一阶差分+LSTM、LSTM、ARIMA和谐波模型的预测精度对比
其中， LSTM网络的超参数设定与上文相同，同时依据AIC最小信息准则[15]将ARIMA模型的p、 d、q值确定为4、0、4。由于谐波模型是一种时间函数模型，因此其预测精度取决于建模的精度，不随预测步长变化。从图8可以看出，本文提出的方法的预测结果优于谐波模型、未差分处理的LSTM模型和ARIMA模型。
图 9 一阶差分+LSTM、LSTM、ARIMA和谐波模型U方向预测精度（预测步长为17 d）
在此基础上，分别利用4种方法对180个测站的非线性变化进行了预测，U方向上的结果如图9所示，其中本文的方法优于其他3种方法，测站数为166个，占测站总数的92%。表1给出了3个方向，4种模型预测结果的统计数据。首先从平均MAE上来看，3个方向上经过一阶差分处理的LSTM模型均小于其他3种模型，说明该模型预测的绝对精度优于其他3种模型。其中，N方向上相比于谐波模型、ARIMA和LSTM模型分别提高了90%、73%和22%，E方向上分别提高了90%、46%和42%，U方向上分别提高了78%、25%和26%。其次，从预测结果的MAE值的标准差SDMAE来看，本文提出的模型同样是4种模型中最低的，这说明该模型在不同测站上的预测精度更加稳定，适用性更好。
表 1 一阶差分+LSTM、LSTM、ARIMA和谐波模型的预测结果（预测步长为17 d） m
模型
N
E
U
平均MAE
SDMAE
平均MAE
SDMAE
平均MAE
SDMAE
谐波模型
0.69
0.46
0.70
0.74
2.10
0.63
ARIMA
0.26
2.06
0.13
0.23
0.61
0.38
LSTM
0.09
0.11
0.12
0.39
0.62
0.41
一阶差分+LSTM
0.07
0.07
0.07
0.11
0.46
0.33
从上述实验结果和分析可知，一阶差分+LSTM的方法无论是在预测的绝对精度上还是适用性上均优于谐波模型、LSTM模型和ARIMA模型，是一种精度更高，适用性更强的坐标非线性变化预测方法。
3 结束语
本文针对坐标非线性变化的预测问题，提出了一种基于一阶差分和LSTM的坐标非线性变化预测算法，实现了坐标非线性变化的高精度预测。实验结果表明。
1）基于SVD和MIC准则的降噪方法能够在优先保证降噪后时间序列具有较强的非线性关系的同时兼顾更多的重构成分和更少的噪声残留，是一种科学有效的降噪方法。
2）提出的坐标非线性变化预测算法通过对非线性变化中难以利用时间函数模型化的部分进行一阶差分处理和预测，极大地提升了非线性变化的预测精度，相比谐波模型、ARIMA模型和未经一阶差分的LSTM方法其精度分别能够提升78%、25%和22%以上。同时，其对不同测站进行预测时预测精度的标准差相比其他几种方法也有了大幅下降，说明该方法不仅具有更高的预测精度，同时也有更好的稳定性和适用性，能够应用于地球参考框架的高精度维持和形变监测研究。
参考文献
[1] 姜卫平, 马一方, 邓连生,等.毫米级地球参考框架的建立方法与展望[J].测绘地理信息, 2016, 41(4):1-6.(JIANG Weiping, MA Yifang, DENG Liansheng, et al. Establishment of mm-level terrestrial reference frame and its prospect[J]. Journal of Geomatics, 2016, 41(4):1-6.)[2] 程鹏飞,成英燕.我国毫米级框架实现与维持发展现状和趋势[J].测绘学报,2017,46(10)：1327-1335.(CHENG Pengfei, CHENG Yingyan. The current status and tendency of China millimeter coordinate frame implementation and maintenance[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1327-1335.)[3] DAM T V, PLAG H P, FRANCIS O, et al. GGFC special bureau for loading: current status and plans[J]. IERS Technical Note, 2002(30):180-198.[4] DONG D, FANG P, BOCK Y, et al. Anatomy of apparent seasonal variations from GPS-derived site position time series[J]. Journal of Geophysical Research: Solid Earth, 2002, 107(B4): 9-16.
[5] HE X X, MONTILLET J P, HUA X H, et al. Noise analysis for environmental loading effect on GPS position time series[J]. Acta Geodynamica et Geomaterialia, 2017,14(1): 131-142.
[6] JIANG W P，LI Z，DAM T，et al．Comparative analysis of different environmental loading methods and their impacts on the GPS height time series[J]．Journal of Geodesy，2013，87(7)：687-703．
[7] DONG D. Anatomy of apparent seasonal variations from GPS-derived site position time series[J]. Journal of Geophysical Research, 2002, 107(B4):2075.
[8] YAN Haoming，CHEN Wu，ZHU Yaozhong，et al．Contributions of thermal expansion of monuments and nearby bedrock to observed GPS height changes[J]．Geophysical Research Letters，2009，36(13)：L13301．
[9] MING F, DANAN D, HAGER B H. Displacements due to surface temperature variation on a uniform elastic sphere with its centre of mass stationary[J]. Geophysical Journal International, 2014(1):194-203.
[10] 马符讯. 欧洲区域测站非线性运动中环境负载的影响研究[D].武汉: 武汉大学,2017:50-61.(MA Fuxun. The research of the influence of environmental loading on non-linear motion of European GPS stations[D]. Wuhan: Wuhan University, 2017:50-61.)[11] 常航. GPS垂向时序中环境负载效应研究[D].西安: 长安大学, 2019:26-51.(CHANG Hang. Research on environmental loading effect in GPS vertical time series[D]. Xi’an: Chang’an University, 2019:26-51.)[12] 姜卫平,王锴华,李昭,等. GNSS坐标时间序列分析理论与方法及展望[J].武汉大学学报(信息科学版), 2018, 43(12):359-370.(JIANG Weiping, WANG Kaihua, LI Zhao, et al. Prospect and theory of GNSS coordinate time series analysis[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 359-370.)[13] 傅彦博. GPS测站坐标非线性变化的统计改正模型与作用机制研究[D]. 郑州:信息工程大学, 2018:26-39. (FU Yanbo. Study on the statistical correction model and the mechanisms for nonlinear variations of GPS station’s coordinate[D]. Zhengzhou: Information Engineering University, 2018:26-39.)
[14] DAVIS J L，WERNICKE B P，TAMISIEA M E．On seasonal signals in geodetic time series[J]．Journal of Geophysical Research：Solid Earth，2012，117(B1)：B01403．
[15] 田亮．GPS测站坐标非线性变化规律分析与机制研究[D]．郑州：信息工程大学，2011：12-40．(TIAN Liang．Studies on GPS station's coordinate non-linear movement and its causes[D]．Zhengzhou：PLA Information Engineering University，2011:12-40．)
[16] 尹玲,尹京苑,孙宪坤,等.缺失GPS时间序列的神经网络补全[J].测绘科学技术学报, 2018, 35(4):331-336.(YIN Ling, YIN Jingyuan, SUN Xiankun, et al. Reconstruction of gappy GPS coordinate time series based on long short-term memory networks[J]. Journal of Geomatics Science and Technology, 2018, 35(4):331-336.)
[17] 王善辉. 双向循环神经网络在GNSS坐标时间序列插值中的研究[D].太原: 太原理工大学, 2019:27-52.(WANG Shanhui. Reconstruction of gappy GPS coordinate time series using bidirectional LSTM neural network[D]. Taiyuan: Taiyuan University of Technology, 2019:27-52.)
[18] WANG Jian，JIANG Weiping，LI Zhao，et al．A new multi-scale sliding window LSTM framework (MSSW-LSTM)：a case study for GNSS time-series prediction[J]．Remote Sensing，2021，13(16)：3328．
[19] HAMILTON J D. 时间序列分析（上册）[M]. 夏晓华, 译. 北京：中国人民大学出版社，2015:52-53.
(HAMILTON J D. Time series analysis [M].XIA Xiaohua, Translate.Beijing: China Renmin University Press, 2015:52-53.)[20] 马俊.GPS坐标时间序列严密三维噪声模型建立方法[J].测绘学报,2019,48(9):1205.(MA Jun. Establishment of strict three dimensional noise model for GPS coordinate time series[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(9):1205.)
[21] 刘树聃，陈知行．奇异值分解和EEMD的非线性振动信号降噪方法[J]．探测与控制学报，2019，41(3)：37-42． (LIU Shudan，CHEN Zhixing．Nonlinear vibration signal de-noising based on singular value decomposition and EEMD[J]．Journal of Detection & Control，2019，41(3)：37-42．)
[22] 邱小梦，王奉伟，周世健，等．局部均值分解和奇异值分解在GNSS站坐标时间序列信号降噪中的应用[J]．测绘通报，2020(5)：85-89． (QIU Xiaomeng，WANG Fengwei，ZHOU Shijian，et al．Application of local mean decomposition and singular value decomposition in noise reduction of GNSS station coordinate time series signal[J]．Bulletin of Surveying and Mapping，2020(5)：85-89．)
[23] RESHEF D N，RESHEF Y A，FINUCANE H K，et al．Detecting novel associations in large data sets[J]．Science，2011，334(6062)：1518-1524．
[24] TIAN F. CSDN: Maximal information coefficient (MIC) [EB/OL].[ 2021-09-01]. https://blog.csdn.net/FontThrone/article/details/85227239, 2021.
[25] BOCK Y, MOORE A, FANG P, et al. Extended solid earth science ESDR system algorithm theoretical basis document[R]. State of California: Scripps Institution of Oceanography& Jet Propulsion Laboratory, 2019.[26] FERLAND R，PIRASZEWSKI M．The IGS-combined station coordinates，Earth rotation parameters and apparent geocenter[J]．Journal of Geodesy，2009，83(3/4)：385-392．
[27] HE Xiaoxing，HUA Xianghong，YU Kegen，et al．Accuracy enhancement of GPS time series using principal component analysis and block spatial filtering[J]．Advances in Space Research，2015，55(5)：1316-1327．
[28] 周茂盛,郭金运,沈毅,等.基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取[J].地球物理学报,2018,61(11):4383-4395.（ZHOU Maosheng, GUO Jinyun, SHEN Yi, et al. Extraction of common mode errors of GNSS coordinate time series based on multi-channel singular spectrum analysis[J]. Chinese Journal of Geophysics ,61(11): 4383-4395.）
[29] 李慧．奇异值分解在时间序列分析中的应用[D]．北京：北京交通大学，2009：28-36．(LI Hui．Singular-value decomposition in time series analysis[D]．Beijing：Beijing Jiaotong University，2009：28-36．)
作者简介：贾彦锋（1996—），男，河南洛阳人，硕士研究生，主要研究方向为导航时空基准。
E-mail：18603795385@163.com
基金项目：国家自然科学基金项目（41804018）
通信作者：朱新慧 副教授 E-mail：gnss_zzzxh@126.com
转自：“测绘学术资讯”微信公众号
如有侵权，请联系本站删除！
