基于X-11-ARIMA 模型在GNSS 定位数据后处理的应用

匡宇龙，雷孟飞

( 湖南联智科技股份有限公司, 长沙 410200 )

0 引 言

导航定位在测绘和监测领域受到广泛应用，其全天候提供位置服务的特性让测绘和监测工作能够实现24 h 实时监控位置信息，获取长时间的位置变化信息或者通过长时间静态定位获取高精度的测绘点位坐标. 在理想情况下，所选取的静态监测点位不存在位移或者位移可以忽略. 但是在地灾监测的点位，比如边坡及桥梁等区域，因其通行车辆的特点使得震动不可避免，如何区分正常波动和灾害发生导致的位移将是导航定位在地灾监测领域中的关键部分[1].

在导航定位后处理算法领域中，可行的数据处理步骤是：首先通过差分定位获得基线向量改正值；然后通过卡尔曼滤波或其改进算法对坐标参数进行滤波平滑处理，能够实现高精度实时数据输出，在监测中能够满足探测位移的发生，并且能够剔除粗差等干扰[2]. 但是以年为单位的数据记录表明，单纯通过一次滤波，或单纯使用基于原始数据的滤波处理虽然能够剔除粗差和异常位移的干扰，但切实发生的长期性波动如车辆通行不同带来的震动幅度变化、潮汐变化、电离层变化、对流层变化、卫星星座分布变化及可能存在的地壳变化带来的干扰则不属于“误差”，而是切实发生的位移，这些变化难以通过滤波算法进行剔除. 而对该数据进行时间序列分析则能够针对这些长期性变化的波动干扰进行分离和处理，并获得在长时间尺度上提高定位精度的效果.

时间序列分析被广泛应用于数据流的分析和预测[3]. 在各种时间序列分析方法中，差分自回归移动平均模型(ARIMA)模型则是由Box 与Jenkins 提出的一种著名时间序列预测模型. 游贤菲等[4]. 利用X-11-ARIMA 模型对手足口疫病发病数预测和张立欣等[5]利用X-11-ARIMA 模型铁路货运量进行分析的应用. 而在GNSS 领域的运用，文献[6-8]利用ARMA 模型改善最小二乘解算观测量数据，文献[8-9]利用ARMA 模型分析预测大气电离层电子总量(TEC). 关于ARIMA 和季节性拆分方法的运用，有HILLMER S C[9]和JUNAIDI[10]以及PENG Y 等[11]运用ARIMA 及数据拆分进行时间序列分析. Xin J Z 等[12]针对桥梁结构监测中卫星定位的应用，采用了卡尔曼滤波结合ARIMA 模型的方法预测桥梁形变，Li Q S 等[13]利用ARIMA-GRACH 模型探测周跳，JANPAULE I 等[14]分析了十年跨度以日为间隔的定位解算数据，指出地磁暴对全球卫星导航系统(GNSS)数据影响的时间跨度为一天.

文中综合分析了国内外研究情况，结合卫星定位监测数据的特点，采用时间序列分析的方式拆分监测数据. 通过ARIMA 算法预测补偿平滑算法数据缺失，得到原始的完整趋势项数据. 通过实验对比分离周期性波动和误差后的数据改善情况.

1 数据处理方案

1.1 基于季节结构的分解方法

在分离出趋势因素Tt后，将乘法模型转换可得Xt/Tt=St·It，即将原序列除去趋势项后便得到季节因素和不规则因素项. 显然继续对St·It进行移动平均，便能够剔除不规则因素项而分离出季节因素项.不同的平均项数的选择将会对剔除结果产生显著的影响，当项数越大则剔除不规则因素的效果越好，但是越大的项数意味着损失更多的信息，尤其以尾项的数据损失最为严重. 所以移动项数既不能过大，也不宜过小. 在X-11 方案中，一般选择(3×3)、(3×5)或者(3×9)项移动平均法，视原始序列的不规则因素大小和分析目的而定. 由于在剔除不规则因素之前，已经通过十二项移动平均法分离趋势项，不规则因素经过两次移动平均处理几乎趋近于零，所以一般不用担心存在不规则因素消除得不够彻底的问题. 对于一组时间序列的原始数据进行分析需要循环多次进行分解运算，即重复上述的分离趋势和剔除不规则项运算，在多次迭代后得到被平滑和拆分的数据，以进行后续的分析预测工作.

1.2 季节性ARIMA 预测模型

时间序列分析模型ARIMA 在GNSS 领域的应用按照其参与的数据运算环节，参考GNSS 于INS融合的类型，将其分为松耦合与紧耦合. 文中提出的利用X-11 方案的ARIMA 模型进行季节性分析后处理平滑，因其不参与原始数据定位解算而定义为松耦合. 在文献[6-7]的研究中， 利用ARMA 模型进行随机模型参数估计，实时修正GNSS 观测数据. ANSARI K 等[8]利用ARMA 模型分析TEC，研究电离层延迟与卫星角度关系，这类运用将ARMA 模型用于协助原始数据解算获取坐标参数，可认为是紧耦合的模式. 从国内外研究可知，GNSS 定位的确存在一定的周期性波动，无论该波动是来自于电离层周期性变化还是其他因素. 于是可以推论基于一种时间序列分析的数据平滑后处理在原理上是可行有效的.

ARIMA 模型建模需要确定两组参数，一组是自回归移动参数 (p,d,q) ，另一组是季节性自回归移动参数 (P,D,Q)S.

通常季节性ARIMA 模型被写作SARIMA(p,d,q)(P,D,Q)S，S表示Seasonal，即季节项数据序列，其数学算式可写作如式(10)形式：

同位素示踪法则是一种以同位素作为示踪物质，对研究对象的特征和行为进行示踪观察的信息获取方法，分为两类： 放射性同位素示踪法和稳定同位素示踪法。

式中：p为自回归阶数；d为差分阶数；q为移动平均阶数；P为季节性自回归阶数；D为季节性差分阶数；Q为季节性移动平均阶数.

式(10)的符号含义如式(11)～(16)所示：

1.3 数据处理

对北斗定位数据进行时间序列分析，将定位数据中的趋势项和周期项以及不规则项通过X-11 方法进行拆分. 其中趋势项数据可以视为剥离干扰和周期性波动的稳定平滑数据. 但是由于X-11 固有特点导致首尾各缺失6 历元数据，无法满足监测的实时性需求.

本算法包括三个部分，首先对原始数据进行拆分，得到原始数据的趋势项数据、季节项数据和不规则项数据；接下来通过ARIMA 算法对趋势项和不规则项数据进行建模，并给出12 个历元的预测值；最后将预测得到的趋势项和不规则项与季节项数据进行结合，季节项相应的12 个历元数值通过季节项自身周期性推导. 趋势项、不规则项以及季节项重组的数据再次进行拆分运算，得到缺失的6 个历元趋势项.对首部缺失的6 个历元同样可以通过该方法进行补充. 由此可以得到原始数据完整历元的趋势项数据，达到剥离干扰并获得贴近真实位移的完整数据. 具体流程如下.

1.3.1 数据拆分

针对北斗基线向量残差观测数据进行X-11 拆分，由于数据过零点且有负值，不适合使用乘法模型，故而采用Xt=Tt+St+It加法模型.

X-11 拆分算法采用(12×2)的模型. 进行两轮运算，首轮运算窗口宽度为12 个历元，从头到尾遍历原始数据. 按式(3)获得新时间序列，排列可记为6.5、7.5、···、[(n–6)+0.5].

接着对新数列进行第二轮运算，按式(5)得到序列T：7、8、···、n–6，即一次拆分的趋势项数据序列，可以看到首尾各有6 个历元的数据缺失.

1.3.2 数据预测

在前一步骤中通过拆分运算依次获得的T趋势项数据序列、S季节项数据序列以及I不规则项数据序列中，均存在由于拆分算法特点而导致的数据缺失问题. 由于北斗观测数据并非真正意义上的周期性时间序列，所以传统时间序列分析中选择滞后历元替代缺失历元的方法并不适用. 为了补充该部分数据，同时保证数据准确性，采用数据预测的方法对缺失数据进行填充.

以趋势项T为例，其数据序列按对应历元记为：7、8、···、n–6. 对序列T进行ACF 和PACF 分析，得到ARMA(p,d,q)模型的对应参数(3,1,3). 结合数据T，通过模型ARMA(3,1,3)对前12 历元和后12 历元进行预测，补足的T数据序列可记为：–5、···、0、1、···、6、7、8、···、n–6、n–5、···、n、n+1、···、n+6. 其中–5、···、0、1、···、6 为前向补足数据，n–5、···、n、n+1、···、n+6 为后向补足数据.

选取补足数据T中序列为：1～n的数据序列与原始序列X进行差分运算，便能获得序列为1～n的混合序列Y. 随后Y序列拆分得到季节项序列S同样首尾缺失6 历元数据，由于季节项序列S具有周期性，便能推导出相应历元的数值，从而得到完整序列，进而与序列Y差分得到完整长度的不规则项序列I. 由于I仅通过差分剥离季节项序列S得到，所以没有数据缺失.

1.3.3 数据重组与再拆分

经过上述步骤得到了与原始数据序列长度相同的趋势项序列T、季节项序列S和不规则序列I，但由于趋势项T仅经过一次预测得到缺失的6 历元数据，其可靠性相对而言不够充分. 所以需要对三组序列进一步延伸，得到未来6 历元的数据. 由于在上一步骤中序列T已经获得了后向12 历元的数据，分为6 历元补足数据，6 历元预测数据. 季节项序列S按其周期性推导6 历元预测数据，而不规则序列I则通过ARMA 预测得到相应的6 历元预测数据.

将得到的三组数据进行叠加还原为原始数据类型，随后对“还原的”原始数据进行拆分运算，得到缺失的6 个预测历元的趋势项序列，易见该数据长度贴合真实数据长度，从而达到了获取当期历元趋势项序列数据的目的. 由于趋势项数据剥离了其他波动，可认为该数据为贴近真实位移的平滑数据.

经实验分析，拆分后预测能够将数据补充到满足拆分后处理得到当下历元数据. 与直接ARIMA 预测对比，精度有一定程度提升.

2 实验分析

2.1 数据获取

本文实验通过基于ublox-f9p 芯片研制的接收机所接收的GNSS 数据进行测试，接收机接收北斗B1、B2、B3；GPS L1、L2 信号，采用单基站实时动态(RTK)工作模式，数据经过自适应卡尔曼滤波处理. 选取3 个项目上的历史数据进行对比测试，如表1 所示，均是经过自适应卡尔曼滤波[14]后的长时间跨度数据样本. 实验分为两大部分，第一部分为X-11 季节性拆分，其原理是加权滑动平均滤波. 展示拆分季节项和不规则项后的趋势曲线相较于仅经过卡尔曼滤波的数据. 第二部分为ARIMA 预测数据，包括直接对原始数据进行ARIMA 预测，以及将原始数据进行季节性拆分后，分别对趋势项和不规则项进行ARIMA 预测，并将预测结果结合季节项还原为原类型数据.

表1 数据样本分类

2.2 X-11 拆分

如图1 所示，Y轴为位移刻度，反映的是基线向量在X方向的坐标残差，单位为mm，X轴为历元，历元间隔为1 h. 图中蓝色曲线为只经过卡尔曼滤波处理的卫星监测数据，红色曲线为在卡尔曼滤波基础之上利用X-11 拆分后得到的趋势项. 由于长时间跨度，图1 中蓝色曲线所代表的卡尔曼滤波本身每个点间隔为1 h. 而蓝色曲线的“粗差”实际上跨度为1～2 个点，即1～2 h，这种情况可以视为动态定位数据本身反映的“位移”而非误差所造成的. 无论引起该误差的原因何在，如此长时间跨度的位置误差无法用短窗口的实时平滑算法处理，如果并非真实位移所引起，则数值会在该类误差源影响消失后回到正常水平. 利用这点便能统计出误差源影响周期，并能够剔除该类长跨度周期性的误差干扰.

图1 数据“波动”的时间跨度

经过卡尔曼滤波数据的方差为4.733，X-11 滤去卡尔曼滤波的季节项和不规则项后数据的方差为2.718. 整体精度提升了42.6%，应对存在的“粗差”也能进行较为明显的消除. 但是由于自回归及滑动平均算法的固有特点，即生成一个平滑数据点需要该点之后6 个数据点的支持(以X-11(12×2)方法为例)，如图2 所示(假设平滑1 000 历元)，反映到实际情况就是6 h 的延迟，这在地灾监控中是无法使用的. 而为了尝试解决该算法固有的滞后性问题，运用ARIMA算法的预测能力，基于现有数据给出需要平滑的历元之后6 历元的数据，以供数据平滑算法保证实时性.

图2 算法原理

2.3 ARIMA 预测分析

按数据不同构建不同的拟合模型，原始数据均是SARIMA 模型，即季节性ARIMA 模型，季节周期设为13，其根据是对原始数据进行自相关函数ACF分析，由图3 中所示，当滞后历元设置为13 时自相关到达一个高峰，随后自相关程度下降，即该时间序列约存在13 历元的周期，判断周期项数值设为13 较为合适.

图3 坐标残差数据ACF 分析

对趋势项和不规则项的建模均为ARIMA 模型，由于原始数据存在零点与负值数据，所以分解采用的方法为相加模型(Additive Decompose). 显然，由于剔除了季节项，所以对趋势项和不规则项的建模不需要引入季节性因素，只需要一般ARIMA 算法进行模型的构建. 构建的参数如表2 所示. 由于季节项数据属于周期性数据，所以不进行预测，直接按历元所处周期进行数据获取.

表2 预测模型

趋势项以及不规则项的预测效果如图4(a)所示.蓝色曲线为原始数据拆分而来的各项数据，包括Trend 趋势项和Resid 不规则项，红色曲线为模型拟合后进行预测给出的数据.

将上述两项数据结合季节项数据后，便能还原为原始数据. 如图4(c)所示，蓝色曲线为原始数据，红色曲线表示对原始数据直接进行SARIMA 建模所输出的预测数据，绿色曲线则表示拆分项分别预测的数据加以整合后还原的数据.

图4 数据预测结果

图4(c)中两种预测方式的性能经量化后，得到的直观数据如表3 所示. 三组数据样本的预测结果分别进行四项误差分析，并对比二者性能差距.

表3 直接预测与拆分后预测效果对比

直接预测是直接对原始数据进行ARIMA 建模并预测，拆分后预测是通过X-11 拆分原始数据为趋势项、季节项和不规则项，然后分别进行ARIMA 建模并预测，将各自预测结果按拆分原理进行还原，得到经拆分后预测的结果.

数据对比方法采用四类方法，分别是平均绝对误差(MAE)、均方误差(MSE)、平均绝对百分比误差(MAPE)和均方根误差(RMSE). 对比直接预测和拆分预测的改进程度采用如下公式计算：

当 Δ 数值为正表明拆分预测较于直接预测没有改善，当 Δ 数值为负时则表示拆分预测对比直接预测误差程度得到了改进. 可以看出三组数据在拟合上均是拆分更优，而拆分预测相较于直接预测都有着一定程度的改善. 在串丝斜坡数据上拆分预测提升显著，整体提升约在50%. 可见在数据存量在300 历元左右时，拆分预测较直接预测有着显著的改善，而2 600历元以上的数据中，拆分预测改进程度较小. 据此可以大致推断用作平滑的数据范围可以选定当前历元往前300 个历元的数据为最佳.

3 结 论

时间序列分析法在其他领域有着成熟的运用[15]，基于该算法的特性，结合GNSS 行业中通过卡尔曼滤波等算法进行数据处理[16]的运用经验和总结出的一些问题，尝试使用本文对时间序列分析法的运用方式解决实时性不足的问题. 理论上经过预测获得所需数据后，能够实现无数据缺失的完整平滑处理. 但是在地灾监测中，重要的是真实和迅速，所以不能单纯依赖预测，本文尚未解决预测数据与实时数据结合并互相修正的问题. 不过根据GNSS 行业其他方向的经验，可参考惯性导航系统(INS)/GNSS耦合方式，将ARIMA 和GNSS 经卡尔曼滤波滤波后的数据进行解算层面的耦合，构建一套以预测数据辅助后处理的运行系统.