CORS站高程非线性速度场及方差波动模型构建方法

张恒璟，崔东东，程鹏飞

1. 辽宁工程技术大学测绘与地理科学学院，辽宁 阜新 123000； 2. 中国测绘科学研究院，北京 100830

CORS站坐标运动特征预测分析是维持2000中国大地坐标系统(China geodetic coordinate system 2000，CGCS2000)框架准确性和现势性的重要基础[1]。通过对国内多个CORS站高程坐标分量序列特征分析发现，CORS站在垂直方向表现出明显的非线性和周期性特征，并且不同的CORS站对应周期性变化模式和运动规律也不尽相同[2]。影响全球定位系统(global positioning system，GPS)高程分量周期性变化的因素很多，主要包括地壳运动、冰雪、土壤水分、大气压负载以及CORS站测量技术的高低等[3-4]。在分析CORS站高程坐标分量运动特征时，只考虑线性运动，直接在线性最小二乘拟合模型中对周期项赋予固定数值[5]，而忽略其非线性运动特征，导致建立的CORS站高程线性速度场模型对于垂向速率估计和误差估计存在很大的偏差[6-8]。同时CORS站周围存在物理、地理因素及测量仪器等产生的噪声，也会影响CORS站高程运动的规律。在对我国CORS站高程时间序列数据噪声模型研究时，文献[9]利用功率谱分析CORS高程序列的噪声性质，并采用极大似然法定量地估计噪声序列中有色噪声的分量。文献[10]利用不同组合的噪声模型对CORS站坐标时间序列的噪声进行分析，并计算了积雪深度、大气压负载、土壤湿度负载等对CORS站位移的影响，得到中国区域CORS站坐标序列的噪声特性主要表现为白噪声+闪烁噪声和白噪声+带通幂律噪声。文献[11]利用AR模型对CORS站高程噪声序列建模，并与CORS站高程非线性拟合模型结合进行高程运动的预报。

研究结果表明，CORS站高程分量建模的残差序列表现为有色噪声，从统计意义上是否为非平稳的随机过程，需要进一步研究。统计学采用异方差法分析时间序列的平稳性特性，如果存在异方差特性，则时间序列是非平稳的随机过程。统计学中对时间序列是否存在异方差特性进行了研究，大多数研究者对截面数据产生的异方差给予了足够多的关注而放松了对时间序列数据也会产生异方差的警惕[12-13]。文献[14]于1982年提出了在时间序列背景下也有可能出现异方差性的观点，并从理论上提出了一种观测时间序列方差变动的方法，自回归条件异方差(autoregressive conditional heteroskedasticity，ARCH)检验。ARCH检验认为时间序列存在的异方差性为 ARCH过程，通过检验该过程是否成立去判断时间序列是否存在异方差[15]。除了ARCH方法检验异方差外，还有残差平方图分析法、帕克检验法、等级相关系数检验法、White检验法等[16]。

国内外CORS站已经积累了20余年的连续观测数据，为研究地壳运动、地震活动等提供了重要的信息。本文采用SOPAC(scripps orbit and permanent array center)全球数据中心6个国内外CORS站20余年的高程时间序列数据，设计非线性周期模型迭代算法，建立CORS站高程坐标分量非线性最小二乘拟合模型，并以均方根误差(root mean square error，RMSE)和拟合优度(goodness of fit，目前没有固定的简写形式，文中记为RNEW)为精度指标衡量模型的优劣。建立CORS站高程非线性拟合模型后，以残差平方序列为研究对象[17]，利用ARCH检验法、残差平方序列图分析法判断残差平方序列的异方差特性，对具有异方差特性的残差平方序列建立广义自回归条件异方差(generalized autoregressive conditional heteroskedasticity，GARCH)波动模型，描述其随着时间变化的非平稳波动情况[18]，验证GARCH模型对非平稳残差平方序列建模的可行性，为CORS站高程残差序列有色噪声建模提供一种途径。

1 CORS站高程时间序列非线性建模方法

y=a+b×t+A1sin(2π×t)+A2sin(4π×t)+A3sin(π×t)

(1)

非线性最小二乘模型将各个周期项作为未知值，加入了每个未知周期项对应的初相和频率未知数，例如式(2)中，M=[abA1f1φ1A2f2

φ2A3f3φ3]T是模型的未知参数。

y=a+b×t+A1sin(2×π×f1×t+φ1)+A2sin(2×

π×f2×t+φ2)+A3sin(2×π×f3×t+φ3)

(2)

线性模型(1)的解作为非线性模型(2)的初值，记为M0。将模型(2)在初值M0处按泰勒级数展开至一阶项，得到线性化后的误差方程式为V=Bm-l。

l=[l1l2…li]

LS求解非线性模型时，给定未知参数初值的精度不同，一次平差结果可能不满足参数的精度要求，采用高斯-牛顿迭代算法，计算流程如图1所示，其中|a|=|an-an-1|。

图1 CORS站高程非线性建模流程Fig.1 The nonlinear modeling process of CORS stations height

为了评价CORS站高程时间序列非线性建模的效果，采用均方根误差RMSE、可决系数(coefficient of determination，记为R2)、拟合优度RNEW 3个精度评价指标。

均方根误差RMSE反映了测量值与真实值之间的偏离程度[19]

(3)

式中，n为测量次数；Δd为测量值与真实值之间的偏差。

可决系数R2是线性回归方程对观测值的拟合程度[20]

(4)

拟合优度在统计学中用于判定非线性回归方程对观测值的拟合程度[21]，在对回归方程拟合程度的解释上，RNEW与R2意义相同，当RNEW越接近于1，非线性模型拟合程度越高

(5)

式中，Q=∑r2；Y=∑y2；r是拟合残差；y是序列观测值。

2 CORS站高程非线性建模试验

从SOPAC全球数据中心获得国内外6个CORS站20余年干净的带有趋势项的高程坐标时间序列数据，对于缺值数据，使用三次样条法进行插值处理[22]。采用线性模型(1)和非线性模型(2)分别建立CORS站高程分量拟合模型。

表1是线性模型(1)的未知参数估计值。作为CORS站高程非线性高斯-牛顿迭代的初值，其中YAR2和AUCK两站只给出了一年和半年周期项的初值，这两个站的两年周期项振幅估计值均小于0.1 mm，两年周期项在高程分量运动模型中所占的比重可以忽略不计。

表1 非线性高斯-牛顿迭代法的初值

表2是非线性模型(2)的未知参数估计值及精度评价指标。为了便于对比线性和非线性模型拟合结果的差异，将线性模型的建模精度指标结果也放在了表2，其中(nl)rmse和(l)rmse分别表示非线性建模与线性建模的均方根误差指标。周期方面：6个CORS站高程非线性建模的周期值(频率f的倒数)均不是整数，其周期性变化并不是经典的CORS站高程时间序列模型所采用的固定整年、整半年或两年；近似年周期的估计误差，除了AUCK站达到12%外，其他5个站均在1%；近似半年周期的估计误差不超过18%；近似两年周期的估计误差不超过6%。振幅方面：近似年周期项对应的振幅值介于5～8 mm，近似半年周期项的振幅介于1～2 mm，近似两年周期项对应的振幅值介于0.1～1.3 mm，说明6个CORS站的高程方向以近似年和半年周期性变化为主，近似两年周期性变化基本忽略不计。精度评价指标方面：6个CORS站的非线性模型RMSE精度在5 mm水平，低于线性模型1～2 mm，线性模型R2明显低于非线性模型RNEW指标0.1左右，非线性建模精度整体优于线性模型。

图2是6个CORS站高程分量线性和非线性建模结果，纵轴表示CORS站高程坐标值，横轴表示高程坐标对应的时间值，纵轴单位是mm，横轴单位是年；绿色散点代表CORS站原始高程值，黑色粗实线代表线性模型拟合曲线，红色粗实线代表非线性模型拟合曲线。从图2可以看出，各个CORS站建立的非线性模型拟合效果优于线性模型。

3 ARCH效应检验与GARCH(p, q)波动模型

为了克服AR模型只能对平稳序列数据建模的不足，引入经济学的序列波动分析方法对非平稳序列数据建模，验证序列波动GARCH分析方法对CORS站高程非平稳残差平方序列建模的可行性，为CORS站高程残差序列建模和非线性速度场重构奠定基础。

3.1 ARCH效应检验方法

异方差性是相对于同方差而言，若一个残差平方序列随时间的变化其值波动幅度较大，则称该残差平方序列具有异方差性。对CORS站高程非线性建模后的残差平方序列检验其异方差性是否存在，即残差平方序列是否存在ARCH效应[23]。图3是CORS站高程残差平方序列ARCH效应检验流程，分为两步：首先检验残差平方序列自相关特性，利用Q检验法及残差平方序列自相关系数(autocorrelation coefficient，ACF)检验如图4所示；接着检验残差平方序列的异方差性，构造t统计量检验流程如图5所示。若CORS站高程非线性建模后的残差平方序列同时存在自相关性和异方差性，则其具有ARCH效应，即CORS站高程非线性建模后的残差平方序列具有非平稳性[24]。

表2 非线性高斯-牛顿迭代法求解出的模型参数值及精度评价值

图2 6个CORS站高程数据建模结果Fig.2 The height data modeling results of 6 CORS stations

图3 CORS站高程残差平方序列非平稳性检验流程Fig.3 The nonstationarity test procedure of CORS stations height residual square sequence

图4 自相关性Q检验流程Fig.4 The autocorrelation of Q test procedure

图4和图5中，滞后数Lags值由残差平方序列自相关系数值对应的最大滞后数确定，试验选择的BJFS和BOGO两个CORS站时间序列长度为21年，其残差平方序列的自相关系数值对应的最大滞后数均为21，即Lags=21。自相关性检验的原假设为残差平方序列不存在自相关性，异方差性检验的原假设为残差平方序列不存在异方差性，显著性水平alpha=0.05，表示残差平方序列相关性和异方差性的检验置信度为95%，h为检验残差平方序列自相关性和异方差性程序中设置的逻辑值，若h=1表示拒绝原假设，h=0表示接受原假设。

自相关性检验的Q统计量可表示为

(6)

式中，ρj是残差平方序列的j阶自相关系数；T是残差平方序列的总长度；p是设定的滞后阶数。

残差平方序列的自相关系数定义为

ACF=r(s,t)/[(DX(t)·DX(s))0.5]

(7)

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

r(s,t)=E[(X(s)-E(X(s)))·(X(t)-

E(X(t)))]

(8)

(9)

残差平方序列的偏自相关系数(partial autocorrelation coefficient，PACF)定义如式(10)，φkk表示滞后数为k的偏自相关系数，由克莱姆法则求得，D为系数行列式，ρi为样本的自相关系数，Dk为将D的第k列换为常数项。

图5 异方差性ARCH检验流程Fig.5 The ARCH test procedure of heteroscedasticity

(10)

图5异方差性t检验时，回归统计矩阵X由残差平方序列值对应的时间构成，反应变量y为残差平方序列值，预测响应yhat为残差平方预测值。t统计量stat定义为

(11)

3.1.1 残差平方序列自相关性检验

自相关性指序列之间不是完全相互独立的，而是存在某种相互关系[25]。BJFS和BOGO两站的残差平方序列如图6所示，纵轴表示高程时间序列数据非线性建模后的残差平方值，横轴表示时间。为了检验CORS站高程残差平方序列的自相关性，绘制两个CORS站的残差平方序列自相关如图7所示，纵坐标ACF表示CORS站残差平方序列的自相关系数，横坐标Lags表示自相关系数对应的滞后阶数值。自相关系数ACF随滞后阶数Lags的变化并不是快速趋近于0，而是缓慢地变化为0，即残差平方序列在不同时刻存在某种相互关系，并不是相互独立的，通过自相关图定性地判断出两站残差平方序列具有自相关性。

图6 两个CORS站高程的残差平方序列Fig.6 The residual squared sequence of two CORS stations height

图7 两个CORS站高程残差平方序列的自相关图Fig.7 The autocorrelation graph of two CORS stations height residual square sequences

为了定量地判断两站残差平方序列的自相关性，采用Q检验法对两站的残差平方序列自相关性进行检验，利用式(7)计算残差平方序列的自相关系数ACF，以最大滞后数(Lags=21)计算式(6)的Q统计量，得到BJFS、BOGO两站的Q统计量分别为4.12和4.23，均大于显著性水平为0.05时的Q临界值3.48，故拒绝原假设，即逻辑量h=1，两站的残差平方序列具有自相关性。

3.1.2 残差平方序列异方差性检验

由残差平方序列图6定性的判断残差平方序列的异方差性[26]：两站的残差平方序列随时间变化，并不是一个常数，并且每隔一年左右都会出现集聚现象(尖峰厚尾)，即残差平方序列随时间波动幅度较大，定性的判断出了两站残差平方序列具有异方差特性。

采用ARCH检验方法定量地判定残差平方序列的异方差特性。ARCH过程如式(12)，i为ARCH过程的阶数；α0为常数项；αi为ARCH过程阶数对应残差平方值的系数；υt为随机误差；α0>0,αi≥0(i=1,2,…,p)

选择2017年4月—2018年4月在本院就诊的膝关节置换术患者330例作为研究对象。根据护理方式的不同，将其分为两组，对照组和观察组。对照组患者中，男82例，女83例，年龄为57～78岁，平均年龄为（72.2±4.6）岁，病程为3～12年，平均病程为（6.3±1.3）；观察组组患者中，男85例，女80例，年龄为59～83岁，平均年龄为（73.2±5.1）岁，病程为2～11年，平均病程为（6.4±1.4）年。两组患者在一般临床资料比较，差异不具有统计学意义（P＞0.05），具有可比性。

(12)

由图5计算两站残差平方序列的t统计量stat：BJFS站stat=2.98、BOGO站stat=2.81，均大于显著性水平为0.05时的t临界值2.10，故拒绝原假设，即逻辑值h=1，残差平方序列存在异方差性，证明了两站高程残差平方序列是非平稳随机过程。

3.2 残差平方序列GARCH(p, q)波动模型

为了定量描述残差平方序列的波动，引入GARCH(p,q)模型对CORS站高程非平稳残差平方序列建模，GARCH模型可以有效地拟合具有长期记忆性的异方差函数，而ARCH模型是GARCH模型的一个特例[27]，即p=0的GARCH(p,q)模型。GARCH (p,q)模型定义如式(13)所示，i=1,2,…,p，j=1,2,…,q，其中p是模型拖尾值，它的阈值由序列自相关系数值基本趋于稳定时对应的滞后数Lags确定；q是模型截尾值，它的阈值由序列偏自相关系数值基本趋于稳定时对应的滞后数Lags确定

(13)

根据残差平方序列自相关图7判断拖尾值p：BJFS站ACF在滞后阶数Lags=6后趋于稳定，模型拖尾值p的阈值为6，即p≤6；BOGO站ACF在滞后阶数Lags=10后趋于稳定，模型拖尾值p的阈值为10，即p≤10。

图8是残差平方序列偏自相关图：纵坐标PACF表示CORS站残差平方序列的偏自相关系数值，横坐标Lags表示偏自相关系数值对应的滞后阶数值。BJFS站PACF在滞后阶数Lags=5后趋于稳定，模型截尾值q的阈值为5，即q≤5。BOGO站PACF在滞后阶数Lags=4后趋于稳定，模型截尾值q的阈值为4，即q≤4。

图8 两个CORS站高程残差平方序列的偏相关图Fig.8 The partial correlogram graph of two CORS stations height residual square sequences

由残差平方序列自相关图和偏自相关图综合判断出了两站GARCH模型的拖尾值p和截尾值q的阈值：BJFS站(p≤6，q≤5)，BOGO站(p≤10，q≤4)。

为了求出GARCH模型最佳的(p，q)值，本文参考文献[28]中引用的赤池信息量准则(Akaike information criterion，AIC)选取最佳模型拖尾和截尾值

AIC=2K-2ln(L)

(14)

式中，K是参数的数量(GARCH(p，q)模型的组合个数)；L是似然函数，即残差平方序列之和除以GARCH(p，q)模型的组合个数。

根据AIC准则求出了最佳(p,q)值，BJFS站(p=1,q=1)，BOGO站(p=1,q=1)，并对两站的残差平方序列建立GARCH(1，1)波动模型

(15)

由极大似然估计准则求解BJFS和BOGO站波动模型GARCH(1，1)的未知参数[29-30]分别见表3和表4，constant代表未知参数M中的常数项α0，GARCH(1)代表β，ARCH(1)代表α，offset代表u，由于4个值为常数或某项的系数，故无单位。

将表3和表4残差平方序列GARCH(1,1)模型参数值分别代入式(15)，解算两站的条件方差序列值，如图9所示，纵轴表示残差平方值与条件方差值，横轴表示时间，黑色实线代表残差平方序列，绿色点线代表条件方差序列。条件方差是在残差平方基础上描述的残差平方序列波动情况，其数值比残差平方值小，为了直观地反映条件方差与残差平方随时间的变化规律，图9将条件方差扩大到100倍，与残差平方值在同一个数量等级。从图9可知，条件方差序列随时间的变化趋势会在某一段时间急剧升高或降低，与残差平方序列随时间的变化趋势一致，利用GARCH(1,1)模型可以对CORS站高程非平稳的残差平方序列波动情况进行描述。

图9 残差平方序列与条件方差序列趋势变化比较Fig.9 Comparison of the trend between the residual squared sequence and the conditional variance sequence

表3 BJFS站残差平方序列GARCH(1, 1)模型参数

Tab.3 GARCH (1, 1) model parameters of BJFS station residual squared sequence

parametervaluestandarderrort statisticconstant0.000000200.000006150.13291100GARCH(1)0.688887000.375217001.83597000ARCH(1)0.167274000.906637000.18450000offset0.000048960.000138280.35401500

表4 BOGO站残差平方序列GARCH(1, 1)模型参数

Tab.4 GARCH (1, 1) model parameters of BOGO station residual squared sequence

parametervaluestandarderrort statisticconstant0.000000200.000001700.11776700GARCH(1)0.521826000.648599000.80454200ARCH(1)0.214017000.0120118017.8172000offset0.000026290.000118600.22171300

4 结 论

(1) 建立了CORS站高程分量非线性速度场周期拟合模型。CORS站高程分量周期性运动以近似年周期项为主，近似半年周期项次之，近似两年变化的比重几乎可以忽略不计。周期项并不是严格的年、半年或两年，固定周期项的线性建模估计结果存在偏差，采用非线性模型对CORS站高程数据建模精度明显优于传统固定周期项的线性模型。

(2) 建立了描述CORS站高程残差平方序列非平稳波动的GARCH(1,1)模型。非线性建模后的CORS站高程残差平方序列自相关性和异方差性检验结果表明，残差平方序列是非平稳的随机过程，即非线性建模后的CORS站高程残差序列具有非平稳性。GARCH(p，q)模型克服了AR模型只能对平稳序列建模的不足，为下一步基于GARCH(p，q)模型对CORS站高程非平稳残差序列建模和非线性速度场重构提供了思路。