上海CORS站高程数据非线性建模方法研究

2021-06-21姜涛

现代测绘 2021年2期

姜涛

(广东省核工业地质局二九二大队，广东河源 517001)

0 引言

国家连续运行参考站(continuously operating referencestations，CORS)目前已在国内大部分省市建立，且积累了大量的连续观测数据[1-2]。这些数据为各省市的测量、导航、地震研究等提供了重要的数据信息，其中CORS站高程坐标数据的信息最为重要[3-5]。为此，研究CORS站高程坐标方向的运动规律具有一定的实际价值和意义。

从国内不同的CORS站高程坐标时序图可以看出，CORS站高程坐标运动规律存在周期性变化，且不同的CORS站其高程坐标运动规律是不同的[6]。造成不同CORS站高程坐标运动周期性变化规律的因素有很多种，其中包括地球潮汐、地壳运动、大气层变化、GPS测量技术等[7]。先前对CORS站高程数据建立数学模型时，常采用线性最小二乘的方法，该方法事先将固定的整年周期值赋给线性最小二乘模型周期项，导致建立的模型周期值与实际周期值存在一定的偏差[8]。这种建立CORS站高程数据数学模型的方法明显地认为CORS站高程坐标运动只存在线性运动，而忽视了其非线性运动特征[9]。

目前GPS测量技术已经应用于交通、建筑、导航、定位等领域，起到了越来越重要的作用[10]。为了对CORS站高程坐标建立符合其实际运动规律的数学模型，本文选取了国内上海CORS站10余年高程数据作为研究对象，对其建立非线性速度场模型[11]。建立非线性速度场模型的前提是先给CORS站高程数据建立线性最小二乘模型，其次求出线性最小二乘模型的解，最后将其作为非线性最小二乘模型的迭代初值进行未知参数的求解，得到反映CORS站高程坐标运动规律的非线性速度场模型。

1 CORS站高程数据建模方法

1.1 线性最小二乘建模原理

通过对CORS站高程坐标时间序列的时序图分析发现，CORS站高程方向运动规律存在明显的年周期运动。由于不同因素的影响，CORS站高程坐标运动周期性变化中可能存在的月、半年、两年变化在时序图中不能被准确发现，并且不同的CORS站其周期性变化规律也不尽相同[12]。为了得到符合各个CORS站高程运动规律的数学模型，需进行多次的数据拟合。包含年、半年周期项的线性最小二乘拟合模型如式(1)所示。

Yi=A+B×ti+C1×sin(2π×ti)+C2×sin(4π×ti)

(1)

式(1)中，Yi为模型拟合值，A为常数项，B为线性速度值，ti代表时间，C1表示一年周期项对应的振幅值，C2表示半年周期项对应的振幅值，周期值分别为T1=1,T2=0.5，该周期值表示线性模型给定的固定周期值为一年和半年。

A,B,C1,C2为4个未知数，设CORS站高程坐标序列观测值为Z。利用间接平差原理得到的误差方程如式(2)所示。

V=AN-Z

(2)

式(2)中，

(3)

(4)

按照最小二乘原理VTPV=min，P为单位矩阵，求解未知参数N的值。

1.2 非线性最小二乘建模原理

非线性最小二乘拟合模型是在线性模型(1)的基础上加入相应的初相值，且周期项不是固定的数值[13]。具体模型表达式如式(5)所示。

Yi=A+B×ti+C1×sin(2π×f1×ti+φ1)+C2×sin(4π×f2×ti+φ2)

(5)

式(5)中A,B,C1,f1,φ1,C2,f2,φ2为未知参数，记为S。

将非线性模型(5)在近似值N0处线性化，线性化后的误差方程如式(6)所示。

V=BS-Z

(6)

式(6)中，S=[ABC1f1φ1C2f2φ2]T，

qi=sin(2×π×ti+φ1)，wi=2×π×ti×cos(2×π×ti+φ1)，ei=cos(2×π×ti+φ1)，

ri=sin(4×π×ti+φ2)，di=4×π×ti×cos(4×π×ti+φ2)，fi=cos(4×π×ti+φ2)，

最后按照最小二乘原理VTPV=min，P为单位矩阵，求解未知参数S的值。

2 CORS站高程数据建模实验分析

本文从SOPAC网站下载上海CORS站2005～2017年高程坐标数据，并绘制上海CORS站高程坐标时序如图1所示。其中，纵轴表示的是上海站高程数据值，单位为米，该高程数据由原始高程数据减去平均值得到。横轴表示的是高程数据对应的观测时间为2005～2017年。从图1上海站高程坐标时序图可以看出：高程坐标数据随着时间的变化具有周期性，且不是线性变化，即不存在直线上升或下降的趋势[14]。这就为高程数据建立数据模型提供了相应的思路。

图1 上海站高程坐标时序图

2.1 线性最小二乘建模

通过对上海站高程坐标时序图分析发现，高程坐标数据随着时间的变化其周期性可能存在一年、半年变化。为此，本文在建立线性最小二乘模型时，首先假定模型周期值为固定一年、半年数值，即T1=1,T2=0.5。然后按照线性最小二乘模型(1)进行高程数据的拟合，得到的上海站高程数据线性最小二乘拟合结果如图2所示，其中纵轴表示高程坐标数据(包括观测高程数据和线性最小二乘拟合建模所求解的高程拟合值)黑色的点表示的是原始高程数据，红色曲线表示的是线性最小二乘拟合曲线图。从图2可以看出，上海站高程坐标运动存在周期性的变化规律，且周期性变化主要是一年变化为主、半年变化占比较少。表1给出了线性最小二乘建模初值、线性最小二乘解及建模精度评价指标均方根误差的值[15]。

图2 上海站高程坐标线性最小二乘拟合图

表1 线性最小二乘建模参数值

由表1线性最小二乘建模参数值可以看出：线性最小二乘拟合建模所需的周期项为给定的周期值，分别是一年、半年数值，由计算出的振幅值C1和C2数值可得，上海CORS站高程数据周期性变化一年周期项对应的振幅值明显大于半年周期项对应的振幅值，这就说明上海CORS站高程数据周期性变化主要为一年变化，半年变化占比非常少。

2.2 非线性最小二乘建模预测

2.2.1 非线性最小二乘建模

前文中已经提到线性最小二乘建模方法中给定固定周期值的方法与实际CORS站高程数据周期性变化值存在一定的偏差，为此，本文建立了非线性最小二乘拟合模型(5)。该模型以线性最小二乘模型的解作为迭代初值，其中非线性最小二乘模型迭代初值如表2所示。

表2 非线性最小二乘建模迭代初值

将表2中给出的非线性最小二乘模型迭代初值代入非线性最小二乘模型(5)，进行不断地迭代解算得到上海站高程数据非线性最小二乘拟合建模结果如图3所示：纵轴代表高程数据值，横轴表示的是高程值对应的时间。

图3 上海CORS站高程坐标非线性速度场拟合模型图

从图3上海CORS站高程坐标非线性速度场拟合模型图可以看出：拟合曲线能够较好地描述CORS站高程数据的变化情况，即拟合曲线能够反映出CORS站高程数据的周期性变化的规律。并且从图3中也可以看出，拟合模型计算出来的值与实际观测的高程值误差在2 cm以内。为了说明非线性速度场模型可以反映CORS站高程数据周期性变化的实际情况，本文计算出非线性最小二乘拟合模型中的未知参数值，如表3所示。从表3可以看出，非线性最小二乘模型中的未知参数f1，f2的数值并不是整数值，即周期性变化的数值并不是固定的整年或整半年数值，这就说明了非线性最小二乘模型将周期项作为未知数进行迭代是符合CORS站高程数据周期性变化规律的。另外，通过对比非线性最小二乘拟合模型残差均方根误差与线性最小二乘拟合模型误差发现，非线性最小二乘拟合模型相较于线性最小二乘拟合模型建模精度提高了32%。

表3 非线性最小二乘建模解算未知参数值

2.2.2 非线性最小二乘预测

为了验证本文提出的非线性最小二乘速度场建模方法的可行性，作者在本小节以上海CORS站2005～2014年的高程数据利用非线性速度场模型进行拟合建模，然后求解出非线性模型的未知参数，得到非线性速度场模型的数学表达式。将上海CORS站2015～2017年的高程数据值作为实际观测值，并与非线性预测模型求解出的预测值进行比较，得到的预测结果如图4所示，前一段绿色的点代表的为2005～2014年的CORS站高程数据值，黑色的曲线代表的为利用非线性最小二乘方法对该段高程数据的拟合。后一段黑色的点代表的是待预测的2015～2017年的CORS站高程数据观测值，红色的曲线为非线性速度场模型预测值。相应的非线性模型参数值和预测结果精度值见表4，其中RMSE/m为预测值于实际值误差的均方根误差统计值。从图4的预测曲线可以看出，预测曲线数据值的变化与CORS站高程数据实际值的变化基本一致，且预测时间越短，预测的精度越高，但从表4可看出总体预测误差在1.7 cm以内，说明本文建立的非线性速度场模型具有一定可行性与正确性。