基于台州及周边地区CORS站网地表形变时空特征分析

徐创富,王恒,吉渊明,刘宇

(1.台州市空间地理信息,浙江 台州 318000;2.丽水市测绘中心,浙江 丽水 323000;3.浙江省测绘科学技术研究院, 浙江 杭州 310000;4.山东科技大学测绘科学与工程学院, 山东 青岛 266000;5.中国测绘科学研究院, 北京 100830)

0 引 言

随着连续运行参考站(CORS)站点观测精度的不断提高,国内外大量学者对CORS站时间序列的分析已经进行了深入的研究,其研究内容主要集中于对时间序列进行拟合确定时间序列模型[1-5],通过滤波等方法研究其共模误差分析其噪声模型,通过功率谱、小波变换等方法研究其周期特征[6-8]．目前大多数研究都是对国际GNSS服务(IGS)站及一些小区域的CORS站进行分析[9-11],对于浙江区域CORS站的分析还非常少．

浙江省卫星导航定位基准服务系统(ZJCORS)于2008年向全省用户提供服务,这些覆盖全省的CORS站已持续运行10年,一直在全天候连续监测地面垂直形变．CORS站点能以毫米级精度全天候、连续监测地表点位的三维位置变化．所以本文基于台州及周边地区15个CORS站2017-2019年的数据,采用线性拟合、功率谱、小波变换等方法对时间序列进行分析．

1 数据来源与预处理

1.1 数据来源

采用GAMIT/GLOBK软件对数据进行处理,首先使用GAMIT软件进行基线解算,解算过程中主要参数历元间隔30 s,固体潮模型采用IERS2010,海潮模型采用FES2004,IGS站先验坐标为SOPAC公布的ITRF2014下坐标．然后利用GLOBK平差软件进行整网平差,得到ITRF2014框架下的站点坐标变化时间序列及运动速度值,在进行解算时,选择了11个IGS 站点参与计算,站点包括:AIRA、ARTU、BJFS、CHAN、DAEJ、LHAZ、PIMO、POL2、SHAO、URUM、YSSK,在计算过程中对于IGS站点E、N、U三个方向的约束量分别为5 mm、5 mm、10 mm,对于CORS站点三个方向的约束均为10 m,图1为CORS站点分布图．

图1 CORS站分布

经GAMIT/GLOBK软件处理之后获得15个CORS站2017年1月1日-2019年8月7日的时间序列,其大部分测站时间序列完整度在95%以上,仅有少部分站时间长度在80%以下．

1.2 数据预处理

通过观察CORS站点时间序列可以发现时间序列存在少量的间断点以及突变点,在对时间序列进行分析时应该首先将突变数据剔除,剔除标准是将大于三倍标准差的值删掉,然后通过插值对数据进行修复,以获得均匀采样的时间序列,以此才能正确地提取时间序列中的周期项并对其分析．

常用的插值方法主要有线性插值和三次样条插值,线性插值及三次样条插值法适用于缺失值较少的时间序列[12],且缺失值越少插值的效果越好,目前有学者新提出正交多项式拟合[12]方法,对于数据连续缺失值较多的问题提供了更好的结局方案,本文根据15个CORS站数据质量情况,分别采用三次样条插值以及正交多项式拟合的方法对数据进行插值处理,

图2为测站数据预处理前后的对比图．图2(a)为预处理前时间序列,图2(b)为预处理后时间序列,因篇幅所限只展示DONT站U方向处理结果,可以发现时间序列中突变点有明显减少,数据间断也得到处理．

(a)DONT站U方向预处理前时间序列

(b)DONT站U方向预处理后时间序列图2 预处理前后时间序列对比

1.3 线性项提取

使用谱分析进行周期探测时,要求时间序列的均值为零,即没有趋势项,否则会影响功率谱分析的结果,根据式(1)对时间序列进行趋势项拟合．

y=a+bt+v，

(1)

式中：y为原始时间序列；a、b为趋势项参数；b为线性速率；t为时间；v为时间序列残差．图3为CORS站时间序列插值与去线性项之后的结果图,红线为插值之后的时间序列,蓝线为线性项,绿线为残差时间序列,由于篇幅所限,只展示DONT站的处理结果,由图3可以看出台州CORS站的线性运动主要存在于水平方向,垂直方向线性运动相较于水平方向较小．去线性项之后DONT站水平东方向时间序列残差范围在-10～10 mm,水平北方向时间序列残差范围在-5～5 mm,垂直方向时间序列残差范围在-20～20 mm．

(a)DONT站E方向时间序列

(c)DONT站U方向时间序列图3 去趋势项时间序列图

2 三维速度场

利用15座CORS站线性项,以时序长度为权,变差函数选择高斯模型,按加权克里金进行格网化,再对格网采用滑动平均滤波,生成的大地高年变率格网模型．如图4所示，为了更好地显示水平速度场特征，绘图过程中，移去了台州市整体水平速度，由图可以看出台州地区地面沉降呈明显的空间分布不均匀特点,从总体上来看,台州市西北部天台临海地区沉降明显,台州东南部地区隆升明显．

图4 大地高年变率格网模型

表1为台州地区附近部分CORS站运动速度统计表,根据统计2017-2019年三年时间中台州及其附近地区各基准站平均运动速度值为北方向-12.9 mm/a,东方向31 mm/a,总体运动趋势为东南方向,这与前人所进行的研究也是相一致的．CORS站高程方向呈不同的运动趋势,LINH站沉降最为明显其沉降速率超过1 mm/a,台州站抬升最为明显,隆升速率超过3 mm/a．

表1 台州CORS站运动速度统计

3 时序周期特征分析

3.1 功率谱分析

功率谱是功率谱密度函数的简称,它定义为单位频带内的信号功率,功率谱表示了信号功率随着频率的变化关系[13-14]．本文使用周期图法对时间序列进行功率谱分析,对于随机信号s(ω)其相应的功率谱表达式如式(2)所示:

(2)

由公式(2)可知,功率谱分析只是在频域对信号进行分析,并不包含任何时频的信息,因此通过功率谱分析仅能获得CORS站时间序列的周期信息,而不能得到其各个周期具体的持续时间以及振幅等信息．图5为DNOT站E、N、U三个方向的功率谱密度图,横坐标为周期,纵坐标为功率谱能量密度．

根据计算,台州地区CORS站E、N、U三个方向均具有一定的周期性,但是U方向上的功率谱能量明显高于E、N方向,约为水平方向的十倍,这说明高程方向上的周期性变化明显大于水平方向．根据图中两个CORS站的功率谱密度可以看出U方向主要的年周期大概在320天左右,半年周期在200天左右,季节周期在90～120天．

图5 功率谱密度

3.2 小波变换

对CORS站进行功率谱分析只能得到其周期信号的频率范围,如果想要得到其精确的周期、振幅及持续时间,我们应该使用小波分析对基准站的数据进行处理,小波分析克服了功率谱分析的缺点,具有多分辨率分析的特点,在时域和频域都有表征信号局部信息的能力,其可以根据信号的频率改变时间分辨率,从而探测正常信号中的瞬态,并展示其频率成分．小波变换的含义是把某一被称为基本小波的函数作位移b后,再在不同尺度a下,与待分析信号f(x)作内积,其原理如式(3)所示:

f(x)∈L2(R).

(3)

上述积分变换为f(x)以ψ(x)为基的积分连续小波变换(CWT),a为尺度因子,表示与频率相关的伸缩,b为时间平移因子．

在实际应用中,连续小波变换需要离散化,相应的离散化小波变换系数为.

(4)

式中,ψj,k为离散小波函数．

根据之前功率谱密度得到的结果,本文将主要的周期信号分为年周期、半年周期以及季节周期,对于其余的短周期以及不规则周期信号不进行处理,图6为通过小波变换所获得的DONT站三个方向各周期振幅及时间信息结果图．

通过小波分析结果图6可以看出,台州地区及其附近CORS站U方向年周期振幅在5～10 mm,波峰出现在每年8月份左右,波谷出现在每年3月份左右,半年周期振幅在0～5 mm,波峰出现在每年6月份和12月份,波谷出现在每年3月份和9月份,季节周期振幅在5～10 mm,波峰出现在每年的3、6、9、12月份,波谷出现在每年的2、5、8、11月份,CORS站水平方向的年周期振幅在0～2 mm,半年周期以及季节周期振幅均在0～1 mm,由此可见,CORS站水平方向上的周期性明显弱于垂直方向,垂直方向相较于水平方向受到地球物理因素的影响更加明显．

根据研究表明,引起CORS站周期性变化的原因非常复杂,通常来说主要分为以下三类,第一类是由日月引力摄动引起的,包括极移、固体潮、海潮以及大气负荷等,第二类是由流体引起的,包括大气压的季节性变化、非潮汐负荷的波动．第三类是由各种误差引起的,包括卫星轨道模型、电离层模型,相位中心的变化等[15-16]．

台州地处浙江东南部沿海地区,属于亚热带季风气候,四季分明,夏季气温炎热,降雨量丰富,冬季气温低,降雨量少,因此,台州CORS站受到的非构造负荷的影响十分明显,CORS站周期性变化更加突出．

图6 小波分析结果图

4 结 论

本文利用线性拟合法、功率谱分析和小波变换技术,对台州及其附近地区15个CORS站2017年1月-2019年8月的连续观测数据进行时间序列分析,结果表明,该地区CORS站的水平运动主要以东南方向的线性运动为主,平均运动速度值为东方向31 mm/a、南方向12.9 mm/a,总体运动速度为33.58 mm/a．在高程方向,台州地区沉降具有明显的区域性特征,具体表现为台州西北部CORS站沉降明显,最大沉降速率超过1 mm/a,台州东南部CORS站抬升明显,最大隆升速率超过3 mm/a,位于台州市区内．