顾及噪声的CORS 站高程非线性建模方法

袁壮志，杨贵军，刘明星

（1.辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000；2.北京农业信息技术研究中心，北京 100097）

0 引言

连续运行参考站(continuously operating reference stations, CORS）的水平与高程坐标分量运动特性预测分析，是建立和维持2000 国家大地坐标系（China geodetic coordinate system 2000,CGCS2000）准确系统框架和现势性的基础[1]。现有的研究结果表明：CORS 站高程序列分量的变化同时具有线性和非线性特征[2]；不仅包含白噪声，也包含有色噪声[3]。因此，建立1 个完整、准确的CORS 站高程分量变化规律模型，需要包含噪声项。

现有文献研究主要集中在对CORS 站高程分量数据建模预测、噪声模型分析等方面。文献[4]引入了奇异谱分析的技术进行高程数据非线性建模，并指出该技术相比地球物理效应分析建模方法具有很明显的优势；文献[5]通过功率谱分析的方法，探测得到CORS 站垂直方向主要存在半年周期运动和年周期运动，在此基础上建立了CORS 站高程坐标分量序列的非线性速度场模型；文献[6]首先利用给定固定周期项的方法，建立线性最小二乘模型，然后将其解作为非线性最小二乘模型的迭代初值，实现了CORS 站高程时间序列数据的非线性周期拟合建模；文献[7]采用功率谱分析高程残差序列的噪声性质，并利用极大似然法定量地估计残差序列中有色噪声的分量；文献[8]首先利用非线性拟合的方法对高程数据建模，然后对未来时间的高程数据进行预测，该预测方法存在的问题是，对高程数据建模时只考虑了周期项和趋势项，忽视了噪声项的影响，使得预测的模型仅适用于短时间内的预测，而在较长时间的预测数据，往往与实际数据会存在较大的误差；文献[9]利用不同组合的噪声模型，分析了CORS 站高程坐标时间序列的噪声，并且考虑了积雪深度、土壤湿度、大气压负载等外界因素对CORS 站位移的影响，其结果证明了中国区域内CORS 站坐标序列的噪声特性主要有2 类，即白噪声与闪烁噪声，白噪声与带通幂律噪声。

截至目前，国内外众多CORS 站已经记录了长达20 余年的连续观测数据，该数据为全球研究地壳板块运动与地质灾害等的学者提供了重要的信息[10]，本论文从斯克里普斯轨道和常驻阵列中心(Scripps Orbit and Permanent Array Center,SOPAC）获取CORS 站高程数据，分别对CORS 站高程分量数据建立非线性最小二乘周期拟合模型（nonlinear least square periodic fitting model, NL）、噪声项自回归（auto-regressive, AR）模型。然后将NL 模型和AR 模型进行组合，对CORS 站的高程坐标序列进行预测[11]。

1 CORS 站高程序列建模预测方法

目前CORS 站高程分量建模方法，包括固定周期项的线性最小二乘建模和将周期值作为未知项的非线性最小二乘建模，这2 种方法主要对高程数据的趋势项和周期项进行建模，未考虑噪声项的影响。本文首先利用非线性最小二乘周期拟合方法对高程数据建模，其次对得到的残差序列（噪声项）建立AR 模型，最后将NL 模型和AR 模型进行组合，对CORS 站的高程坐标序列进行预测，后文中将组合模型记为NL/AR 模型。

1.1 非线性最小二乘周期拟合建模方法

非线性最小二乘周期拟合模型是将各个周期项作为未知值，加入了每个未知周期项对应的初相和频率未知数，得到非线性最小二乘拟合模型为

式中：a 为常数项；b 为线性速度值；t 为CORS 站高 程 值 对 应 的 单 日 历 元 时 间； A1、 A2、 A3分 别 对 应了各自周期项的振幅值； f1、 f2、 f3分别对应了各自周期项的频率值； φ1、 φ2、 φ3分别对应了各自周期项初相值。

将线性最小二乘模型的解作为非线性最小二乘周期拟合模型式（1）的初始值，记为H0，并将非线性最小二乘模型式（1）在初始值H0处，按照泰勒级数展开至1 阶项，即得到线性化之后的误差方程式为

然后根据最小二乘原理TV PV =min，求解出未知参数。

1.2 噪声项AR 建模方法

对利用非线性最小二乘周期拟合模型计算得到得残差序列（噪声项），建立AR 模型，具体表达式为

式中： Wt为噪声序列的实际值； εt为误差；α 、β 为对应的系数项；p、q 为利用自相关函数和偏自相关函数确定，残差序列的自相关函数定义为

式中：ACF 表示残差平方的自相关系数值； r ( s, t)为序列自协方差； D( X (t )) 和D( X ( s)) 表示不同时刻方差。其中序列自协方差、不同时刻方差定义为：

偏自相关函数的数学表达式定义为

式中：kkφ 表示滞后数为k 的偏自相关函数值；D 为系数行列式；kD 为将D 的第k 列换为常数项，ρi为样本的自相关系数。

本文采用均方根误差（root mean squared error，RMSE）作为指标，来评价NL/AR 模型的高程预测效果。

RMSE 反映了预测值与真实值之间的偏差程度[12]，其计算方法为

式中：n 为预测数据个数； dΔ 为预测值与真实值之间的偏差程度。

2 数据预处理

本文从SOPAC 网站上下载了近20 余年的高程分量数据。发现每个CORS 站高程数据均有缺值现象，利用3 次样条法[13]对缺省的数据进行插值处理。插值情况如下：BJFS 站插值110 个，WUHN 站插值210 个，BARH 站插值150 个，IENG 站插值189 个。对插值后的数据进行粗差剔除，利用中误差公式剔除了BJFS 站的12 个数据、WUHN 站的25 个数据、BARH 站的18 个数据、IENG 站的14 个数据。剔除后的数据作为后续建模预测的研究对象，图1 给出了剔除数据的结果，其中黑色小圆圈代表超出3 倍中误差的高程数据（粗差）。

图1 粗差剔除结果

3 实验结果与分析

本文以进行插值和粗差剔除后的CORS 站高程数据为研究对象，分别做了BJFS 站、WUHN 站、BARH 站、IENG 站等4 个CORS 站的高程数据非线性最小二乘周期拟合建模、噪声项AR 建模、NL/AR 模型的高程数据预测。最后对比了NL/AR模型与未考虑噪声项的非线性最小二乘周期拟合模型预测高程数据的效果。

3.1 非线性最小二乘周期拟合建模分析

BJFS 站建模区间（单日历元时间）为1999.804～2018.741（单位为 d），WUHN 站为 1999.067～2013.999，BRAH 站为1998.752～2018.673，IENG 站为2004.001～2018.640。表1 为4 个CORS 站高程数据非性最小二乘周期拟合建模得到的均方根误差和未知参数值。

表1 非线性建模结果值

图2 为建模结果，其中纵轴代表CORS 站高程坐标，单位为mm；横轴代表与各个CORS 站高程坐标所相应的时间；灰色散点代表CORS 站原始高程值，黑色粗线代表模型拟合得到的高程值对应曲线，4 个CORS 站的高程数据建模精度在6 mm 左右。

图2 非线性建模结果值

3.2 噪声项AR 建模预测

3.1 节建立的非线性最小二乘周期拟合模型实际上是对CORS 站高程数据中的周期项、趋势项进行的建模分析，未对噪声项建模，故模型的残差项（真实高程值与拟合高程值之差）即为噪声项。利用AR 对4 个CORS 站高程数据的噪声序列进行建模预测：表2 给出了建模预测精度值；图3 给出了预测时间分别为0.5、1 a 的残差序列预测结果，纵轴为残差真值和预测值，单位是 m，横轴为预测时间，其中残差真值为非线性最小二乘周期拟合建模得到的残差值，残差预测值表示利用AR模型对残差序列建模得到的预测值。

表2 AR 建模精度值

图3 残差序列AR 建模结果值

3.3 高程数据预测

本节对CORS 站高程数据建立了NL/AR 预测模型，并与非线性最小二乘周期拟合模型预测进行对比，预测时长为0.5 和1 a。表3 为0.5 a 和1 a 预测精度结果精度值，评价指标为RMSE，其中NL表示非线性最小二乘周期拟合预测，NL/AR 表示组合模型预测。图4 给出了非线性最小二乘周期拟合和NL/AR 预测模型的结果，纵轴表示高程真实值和2 种不同方法的预测高程值，单位是mm，横轴表示预测时间。由表3 的预测结果精度值统计可得：对于2 种预测方法，预测区间越长，相应的预测结果精度越低，反之，预测时间短，则预测精度高；采用NL/AR 模型的预测方法相对于NL模型的预测提高了精度（RMSE 越小，精度越高），说明文章提出的顾及噪声的CORS 站高程非线性建模方法是可行的，由于加入了噪声项的AR 模型，使得高程预测模型更加完整（包含了周期项、趋势项、噪声项）。

图4 4 个CORS 站高程数据预测结果

表3 预测结果精度值（RMSE）

4 结束语

本文以国内外4 个CORS 站近20 余年的高程时间数据为研究对象，提出了1 种顾及噪声的CORS站高程非线性建模方法，分别建立了非线性最小二乘周期拟合模型和AR 模型，最后将2 种模型组合进行CORS 站高程数据的预测，得出了如下结论：

1）未考虑噪声项的非线性最小二乘周期拟合模型，对CORS 站高程数据建模是不完善的，因为高程数据中一定会含有噪声，建立包含噪声项、周期项、趋势项数学模型才能进行高程数据的准确预测。

2）从4 个CORS 站预测结果的精度值统计可知，提出的NL/AR 模型对CORS 站高程数据的预测较NL 模型预测精度高。