|
《测绘学报》 构建与学术的桥梁 拉近与权威的距离 CORS站高程非线性速度场及方差波动模型构建方法 张恒璟1,2 1. 辽宁工程技术大学测绘与地理科学学院, 辽宁 阜新 123000; 收稿日期：2019-01-14；修回日期：2019-05-30 基金项目：高分辨率对地观测系统重大专项(42-Y20A09-9001-17/18);辽宁省教育厅青年项目(LJ2017QL008);中国测绘科学研究院基本科研业务费(AR1918);武汉大学地球空间环境与大地测量教育部重点实验室开放基金(18-01-05) 第一作者简介：张恒璟(1982-), 男, 博士, 副教授, 研究方向为空间大地测量数据处理与GNSS高程非线性运动。E-mail:sun_winter2009@163.com 摘要：本文描述了CORS站高程时间序列采用固定周期项建立线性速度场的基本概念，指出给定周期与实际周期存在偏差的问题。为此，提出了CORS站高程时间序列非线性速度场建模方法。该方法以线性最小二乘解作为迭代初值，利用高斯-牛顿迭代算法求解非线性速度场模型的未知参数，实现了CORS站高程时间序列数据的非线性建模; 给出了拟合模型残差平方序列异方差特性的检验方法，阐述了建立GARCH(p，q)模型反映CORS站高程非平稳序列波动情况的基本准则。本文以国内外6个CORS站20余年的高程时间序列为研究对象，建立非线性速度场运动模型。结果表明:CORS站高程运动并不存在严格的整年或半年周期项，近似年周期运动最为明显，近似两年周期运动占比最小，年周期估计偏差达到12%，半年周期偏差18%，两年周期偏差6%，非线性建模精度和效果整体优于线性模型。利用ARCH检验法得到CORS站高程非线性模型的残差平方序列存在异方差性，即残差平方序列具有非平稳特性，引入GARCH(p，q)模型对CORS站高程分量非平稳残差平方序列建模，反映了残差平方序列的非平稳波动。验证了GARCH(p，q)模型对CORS站高程非平稳残差平方序列建模的可行性，为今后利用GARCH(p，q)模型对CORS站高程非平稳噪声序列建模和非线性速度场重构提供了思路。 关键词：CORS站高程时间序列 非线性速度场 异方差 GARCH模型 Height nonlinear velocity field and variance fluctuation model construction method for CORS stations ZHANG Hengjing1,2 1. School of Geomatics, Liaoning Technical University, Fuxin 123000, China; Foundation support: The Major Project of High Resolution Earth Observation System (No. 42-Y20A09-9001-17/18); The Youth Project of Liaoning Education Department (No. LJ2017QL008); The Basic Research Fund of Chinese Academy of Surveying and Mapping(No. AR1918); The Open Fund of the Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University(No. 18-01-05) First author: ZHANG Hengjing(1982—), male, PhD, associate professor, majors in spatial geodetic data processing and nonlinear motion of GNSS height time series. E-mail:sun_winter2009@163.com. Abstract: The basic concept of establishing linear velocity field with fixed period term in CORS stations height time series is described, and the problem of deviation between given period and actual period is pointed out. Therefore, this paper proposes a nonlinear velocity field modeling method for CORS stations height time series. This method takes the linear least squares solution as the initial value of iteration, and uses the gauss-newtoni-teration algorithm to solve the unknown parameters of the nonlinear velocity field model, realizing the nonlinear fitting of height time series data of CORS stations. The test method for the heteroscedasticity of the residual square sequence of the fitting model is given, and the basic criteria for the establishment of GARCH (p, q) model to reflect the fluctuations of non-stationary sequence are expounded. The six CORS stations at home and abroad more than 20 years as the research object, the height of time series nonlinear velocity motion model is set up, the results show that the CORS stations height movement does not exist strict year or half year cycle, the approximate periodic motion is most obvious, the approximate period of two years than the minimum, cycle in deviation of 12%, half year cycle is 18%, two year period is 6%, nonlinear modeling accuracy and effect is better than the linear model as a whole. The ARCH test method is used to obtain the heteroscedasticity of the residual square sequence of the nonlinear model of CORS stations height, that is, the residual square sequence has non-stationary characteristics. GARCH (p, q) model is introduced to model the non-stationary residual square sequence of the height component of CORS stations, which reflects the non-stationary fluctuation of the residual square sequence. The feasibility of GARCH (p, q) model in modeling non-stationary residual square series of CORS stations height is verified, which provides an idea for future modeling of non-stationary noise series of CORS stations height and reconstruction of nonlinear velocity field with GARCH (p, q) model. Key words: height time series of CORS stations nonlinear velocity field heteroskedasticity GARCH model 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站高程时间序列非线性建模方法 经典的CORS站高程时间序列线性周期拟合模型采用固定周期项。例如式(1)包含了年、半年和两年周期项，对应周期值分别为T1=1、T2=0.5、T3=2，利用线性最小二乘法求解模型未知参数M =[a b A1 A2 A3]T, 其中a为常数项，b为线性速度，A1、A2、A3为对应不同周期项的振幅值。
非线性最小二乘模型将各个周期项作为未知值，加入了每个未知周期项对应的初相和频率未知数，例如式(2)中，M =[a b A1f1 φ1 A2 f2 φ2 A3 f3 φ3]T是模型的未知参数。
线性模型(1)的解作为非线性模型(2)的初值，记为M0。将模型(2)在初值M0处按泰勒级数展开至一阶项，得到线性化后的误差方程式为 LS求解非线性模型时，给定未知参数初值的精度不同，一次平差结果可能不满足参数的精度要求，采用高斯-牛顿迭代算法，计算流程如图 1所示，其中
为了评价CORS站高程时间序列非线性建模的效果，采用均方根误差RMSE、可决系数(coefficient of determination，记为R2)、拟合优度RNEW 3个精度评价指标。 均方根误差RMSE反映了测量值与真实值之间的偏离程度[19]
式中，n为测量次数；Δd为测量值与真实值之间的偏差。 可决系数R2是线性回归方程对观测值的拟合程度[20]
式中， 拟合优度在统计学中用于判定非线性回归方程对观测值的拟合程度[21]，在对回归方程拟合程度的解释上，RNEW与R2意义相同，当RNEW越接近于1，非线性模型拟合程度越高
式中，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 非线性高斯-牛顿迭代法的初值Tab. 1 The initial value of nonlinear Gauss-Newton iterative method
表选项 表 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 非线性高斯-牛顿迭代法求解出的模型参数值及精度评价值Tab. 2 The model parameter value and precision evaluation value by nonlinear Gauss-Newton iterative method
表选项 图 2是6个CORS站高程分量线性和非线性建模结果，纵轴表示CORS站高程坐标值，横轴表示高程坐标对应的时间值，纵轴单位是mm，横轴单位是年；绿色散点代表CORS站原始高程值，黑色粗实线代表线性模型拟合曲线，红色粗实线代表非线性模型拟合曲线。从图 2可以看出，各个CORS站建立的非线性模型拟合效果优于线性模型。
为了克服AR模型只能对平稳序列数据建模的不足，引入经济学的序列波动分析方法对非平稳序列数据建模，验证序列波动GARCH分析方法对CORS站高程非平稳残差平方序列建模的可行性，为CORS站高程残差序列建模和非线性速度场重构奠定基础。 3.1 ARCH效应检验方法 异方差性是相对于同方差而言，若一个残差平方序列随时间的变化其值波动幅度较大，则称该残差平方序列具有异方差性。对CORS站高程非线性建模后的残差平方序列检验其异方差性是否存在，即残差平方序列是否存在ARCH效应[23]。图 3是CORS站高程残差平方序列ARCH效应检验流程，分为两步:首先检验残差平方序列自相关特性，利用Q检验法及残差平方序列自相关系数(autocorrelation coefficient，ACF)检验如图 4所示；接着检验残差平方序列的异方差性，构造t统计量检验流程如图 5所示。若CORS站高程非线性建模后的残差平方序列同时存在自相关性和异方差性，则其具有ARCH效应，即CORS站高程非线性建模后的残差平方序列具有非平稳性[24]。
图 4和图 5中，滞后数Lags值由残差平方序列自相关系数值对应的最大滞后数确定，试验选择的BJFS和BOGO两个CORS站时间序列长度为21年，其残差平方序列的自相关系数值对应的最大滞后数均为21，即Lags=21。自相关性检验的原假设为残差平方序列不存在自相关性，异方差性检验的原假设为残差平方序列不存在异方差性，显著性水平alpha=0.05，表示残差平方序列相关性和异方差性的检验置信度为95%，h为检验残差平方序列自相关性和异方差性程序中设置的逻辑值，若h=1表示拒绝原假设，h=0表示接受原假设。 自相关性检验的Q统计量可表示为
式中，ρj是残差平方序列的j阶自相关系数；T是残差平方序列的总长度；p是设定的滞后阶数。 残差平方序列的自相关系数定义为
式中，ACF表示残差平方的自相关系数值；r(s, t)为序列自协方差；DX(t)和DX(s)表示不同时刻方差。其中，序列自协方差、不同时刻方差定义为
残差平方序列的偏自相关系数(partial autocorrelation coefficient，PACF)定义如式(10)，φkk表示滞后数为k的偏自相关系数，由克莱姆法则求得，D为系数行列式，ρi为样本的自相关系数，Dk为将D的第k列换为常数项。
图 5异方差性t检验时，回归统计矩阵X由残差平方序列值对应的时间构成，反应变量y为残差平方序列值，预测响应yhat为残差平方预测值。t统计量stat定义为
式中，X为样本平均数；μ为总体残差平方序列值的平均数；σx为样本标准差；n为样本容量。本文研究中样本容量选取为总体的0.2倍。 3.1.1 残差平方序列自相关性检验 自相关性指序列之间不是完全相互独立的，而是存在某种相互关系[25]。BJFS和BOGO两站的残差平方序列如图 6所示，纵轴表示高程时间序列数据非线性建模后的残差平方值，横轴表示时间。为了检验CORS站高程残差平方序列的自相关性，绘制两个CORS站的残差平方序列自相关如图 7所示，纵坐标ACF表示CORS站残差平方序列的自相关系数，横坐标Lags表示自相关系数对应的滞后阶数值。自相关系数ACF随滞后阶数Lags的变化并不是快速趋近于0，而是缓慢地变化为0，即残差平方序列在不同时刻存在某种相互关系，并不是相互独立的，通过自相关图定性地判断出两站残差平方序列具有自相关性。
为了定量地判断两站残差平方序列的自相关性，采用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)
由图 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确定
根据残差平方序列自相关图 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。
由残差平方序列自相关图和偏自相关图综合判断出了两站GARCH模型的拖尾值p和截尾值q的阈值:BJFS站(p≤6，q≤5)，BOGO站(p≤10，q≤4)。 为了求出GARCH模型最佳的(p，q)值，本文参考文献[28]中引用的赤池信息量准则(Akaike information criterion，AIC)选取最佳模型拖尾和截尾值
式中，K是参数的数量(GARCH(p，q)模型的组合个数)；L是似然函数，即残差平方序列之和除以GARCH(p，q)模型的组合个数。 根据AIC准则求出了最佳(p, q)值，BJFS站(p=1, q=1)，BOGO站(p=1, q=1)，并对两站的残差平方序列建立GARCH(1，1)波动模型
式中，yt=u+ut为均值方程；ut=εt－θεt－1方差方程；zt为服从标准正态分布的随机变量；σt2为εt的条件方差；εt－12为ARCH项；σt－12为GARCH项；α0为常数项；u是残差平方序列均值；α为ARCH项系数；β为GARCH项系数；模型待求的未知参数是M =(α0, α, β, u)。 由极大似然估计准则求解BJFS和BOGO站波动模型GARCH(1，1)的未知参数[29-30]分别见表 3和表 4，constant代表未知参数M中的常数项α0，GARCH(1)代表β，ARCH(1)代表α，offset代表u，由于4个值为常数或某项的系数，故无单位。 表 3 BJFS站残差平方序列GARCH(1, 1)模型参数Tab. 3 GARCH (1, 1) model parameters of BJFS station residual squared sequence
表选项 表 4 BOGO站残差平方序列GARCH(1, 1)模型参数Tab. 4 GARCH (1, 1) model parameters of BOGO station residual squared sequence
表选项 将表 3和表 4残差平方序列GARCH(1, 1)模型参数值分别代入式(15)，解算两站的条件方差序列值，如图 9所示，纵轴表示残差平方值与条件方差值，横轴表示时间，黑色实线代表残差平方序列，绿色点线代表条件方差序列。条件方差是在残差平方基础上描述的残差平方序列波动情况，其数值比残差平方值小，为了直观地反映条件方差与残差平方随时间的变化规律，图 9将条件方差扩大到100倍，与残差平方值在同一个数量等级。从图 9可知，条件方差序列随时间的变化趋势会在某一段时间急剧升高或降低，与残差平方序列随时间的变化趋势一致，利用GARCH(1, 1)模型可以对CORS站高程非平稳的残差平方序列波动情况进行描述。
(1) 建立了CORS站高程分量非线性速度场周期拟合模型。CORS站高程分量周期性运动以近似年周期项为主，近似半年周期项次之，近似两年变化的比重几乎可以忽略不计。周期项并不是严格的年、半年或两年，固定周期项的线性建模估计结果存在偏差，采用非线性模型对CORS站高程数据建模精度明显优于传统固定周期项的线性模型。 (2) 建立了描述CORS站高程残差平方序列非平稳波动的GARCH(1, 1)模型。非线性建模后的CORS站高程残差平方序列自相关性和异方差性检验结果表明，残差平方序列是非平稳的随机过程，即非线性建模后的CORS站高程残差序列具有非平稳性。GARCH(p，q)模型克服了AR模型只能对平稳序列建模的不足，为下一步基于GARCH(p，q)模型对CORS站高程非平稳残差序列建模和非线性速度场重构提供了思路。 【引文格式】张恒璟, 崔东东, 程鹏飞. CORS站高程非线性速度场及方差波动模型构建方法. 测绘学报，2019，48(9)：1096-1106. DOI: 10.11947/j.AGCS.2019.20190017 |
|
|
, 崔东东1,2, 程鹏飞2 
(1)
(2)
。
。
(3)
(4)
表示拟合值；yi表示CORS站高程原始数据；y表示CORS站高程原始数据平均值。R2越接近于1，表示线性模型对原始数据拟合程度越高。
