2. 南京智绘星图信息科技有限公司,南京市玄武大道699号,210023
基于长期观测获得的GPS坐标时间序列能够确定GPS速度场以及地球物理过程导致的非线性变形,因此GPS观测是现今地壳运动监测和研究的主要手段之一。垂向分量通常由地球表面大气、流体等负载变化引起,呈现显著的季节性变化[1]。大量研究表明,中国大陆垂向线性运动主要表现为山区上升、盆地和平原下降[2-3]。符养[4]分析全球178个基准站垂向时间序列的线性和季节性变化发现,地球垂向上存在周期性震荡运动;王敏[2]发现1998~2007年除东北和华北地区外,中国大陆在垂向上均存在上升趋势,青藏高原在整体隆升。随着陆态网络二期工程的实施,中国大陆连续站增加至260多个。畅柳等[5]选择263个站点数据分析中国连续基准站坐标时间序列3个分量不同周期的振幅在空间上存在的规律性差异,研究表明,垂向周期振幅与相位的空间分布具有一定的相关性与波动性;占伟[6]选取235个连续站2010~2016年数据对中国大陆垂向坐标时间序列进行研究,得到中国大陆年周期振幅的空间特征为“两边小中间大”、半年周期运动性较弱等结论;傅彦博等[7]基于461个GPS站分析E、N、U三个方向上非线性变化的振幅和初相全球分布,发现周年项振幅随纬度的分布规律,并建立了GPS测站三维周年变化统计的改正模型。
目前GPS时间序列垂向线性速率、周期项振幅和相位空间分布特征的相关研究较多,但没有建立基于中国大陆的GPS测站周年变化改正模型。因此,本文基于中国大陆构造环境监测网络数据,分析趋势项及季节项振幅的分布特征。
1 数据预处理 1.1 GPS数据处理本文选用中国大陆构造环境监测网络260个基准站数据(图 1),利用美国航空航天局NASA的喷气推进实验室(Jet Propulsion Laboratory, JPL)开发研制的GIPSY-OASIS软件包,采用精密单点定位PPP方法计算获得高精度坐标时间序列。数据处理策略如下:首先采用GIPSY-OASIS和PPP模式进行每日24 h时段的处理,获得单日坐标松弛约束解;然后采用st_filter对所有站点的单日松弛约束解进行联合求解,通过Helmert转换将其整体旋转、平移至统一的ITRF2014参考框架中;最后获得各站点的坐标时间序列。
由于天气变化、树木遮挡等原因,各站坐标时间序列中存在一些明显偏离线性趋势或周期变化的数据。为得到准确的趋势项和周期项参数,采用3σ准则进行剔除。图 2为260个站点的历元数与野值剔除率,剔除野值后,所有站点的数据完整率大于99%,不影响后续的拟合分析。
将时间序列分解为相互独立的趋势项、周期项和残差项:
$ Y=Y_k+Y_t+v_t $ | (1) |
式中,Y为历元坐标,Yk为趋势项,Yt为季节项,vt为不规则残差。
2.1 线性拟合获取趋势项由于对时间序列周期项进行分析时要求去除趋势项和阶跃项,因此需要对去除野值后的时间序列进行最小二乘法拟合,其拟合模型为:
$ y_t=a t+b+x_t+\sum\limits_{j=1}^{n_j} g_i H\left(t_i-T_{g_j}\right) $ | (2) |
式中,yt为坐标时间序列,t为时间,a为线性变化率,b为初始位置,xt为时间序列季节项,gi为各种原因引起的坐标突变,Tgj为坐标发生突变的时刻。阶跃发生前H取值为0,阶跃发生后H取值为1。拟合后在时间序列中减去拟合值,得到用于小波分析的新时间序列。以YNYL站为例,图 3为YNYL站的垂向时间序列。
小波分析能够在时域和频域上同时进行,并且同时具有时间和频率特性。小波分析的基础是小波变换,小波变换就是通过小波函数关系去表述或逼近一个信号,该信号能够由基小波函数通过平移或伸缩构成[8]。设φ(t)∈L2(R)(L2(R)为平方可积的实数空间,即能量有限的信号空间),其傅里叶变换为φ(t),当φ(t)满足式(3)的条件时,则称φ(t)为一个基本小波:
$ C_{\varphi}=\int_R \frac{|\varphi(t)|^2}{|\omega|} \mathrm{d} \omega <\infty $ | (3) |
对基本小波进行平移或伸缩,得到一个小波序列:
$ \varphi_{a b}(t)=\frac{1}{\sqrt{|a|}} \varphi\left(\frac{t-b}{a}\right), a, b \in R, a \neq 0 $ | (4) |
式中,a为伸缩因子,b为平移因子。本文选取coif 5小波对时间序列数据进行小波分解。以YNYL站为例,将去除趋势项的时间序列解算为8层,结果如图 4所示,其中d1~d6为分解后的低频项。
由图 4可见,通过小波分析分解的d7、d8信号中可以得到较为稳定的周期信号,分别对应着半年周期信号和年周期信号。YNYL站小波分析获得的年周期信号振幅为8~12 mm,半年信号振幅为2~5 mm,YNYL站年振幅最大值出现在5月份,最小值出现在6月,这可能与大气压以及其他季节性因素的变化周期有关。对文中所有站点数据进行相同的处理,可获得每个站点的线性速率及周期特征。
2.3 周期项时间序列振幅拟合GPS时间序列通常将周期信号视为具有年周期和半年周期的简谐波信号,拟合模型为:
$ \begin{gathered} U_t=c \cdot \sin (2 \pi t)+d \cdot \cos (2 \pi t)+ \\ e \cdot \sin (4 \pi t)+f \cdot \cos (4 \pi t)+v_t \end{gathered} $ | (5) |
式中,c、d为周年项系数,e、f为半周年项系数,vt为残差项。时间序列中的周年项、半周年项的振幅为:
$ A_{\mathrm{ann}}=\sqrt{c^2+d^2}, A_{\mathrm{semi}-\mathrm{ann}}=\sqrt{e^2+f^2} $ | (6) |
采用RMSE来衡量拟合效果,RMSE表示为:
$ \mathrm{RMSE}=\sqrt{\frac{1}{N} \sum\limits_{i=1}^N\left(y_i-\bar{y}_i\right)} $ | (7) |
RMSE值越小表示其拟合效果越好。以YNYL站为例,其周期项拟合结果如图 5所示。为获取周期信号振幅的全国分布特征,对小波分析提取的主周期信号进行常振幅谐波拟合获取振幅值。
图 6为各站对主周期项谐波拟合的RMSE。在所有站点垂向周期信号拟合中,RMSE最大为6.842 mm(GDST),最小为0.463 mm(LHAS),RMSE均值为3.65 mm,超过80%站点的拟合RMSE小于3 mm,上述模型能够较好地反映GPS垂向时间序列季节变化情况。
获得我国各基准站垂向坐标时间序列的趋势项及周期项振幅后,利用克里金插值法对各站数据进行空间插值获得各参数的空间分布图。
3.1 垂向速率分布特征由图 7可见,在260个站点中,90%站点的垂向年变化速率范围为-2~2 mm/a,31%站点的年速率为负,69%站点的年速率为正。塔里木盆地中部有轻微下沉现象,两侧有一定程度的隆升,青藏块体南部有大幅度隆升,可能与欧亚板块和印度洋板块的挤压有关。XZZF(西藏珠峰)、XZCY(西藏察隅) 的线性速率分别为1.04 mm/a和1.60 mm/a,青藏块体的东北部有隆升现象,与占伟[6]结论的一致性较高(1.89 mm/a、1.19 mm/a)。西域块体中天山地区在逐年上升,天山范围内的站点XJZS(新疆昭苏)、XJBY(新疆巴音布鲁克)、XJXY(新疆新源)、XJWL(新疆乌拉斯台)的垂向线性速率分别为0.80 mm/a、0.73 mm/a、0.34 mm/a、1.02 mm/a,略微大于占伟[6]的结论(0.59 mm/a、0.55 mm/a、-0.09 mm/a、-0.10 mm/a),略微小于王伟[3]得到的平均速率(1.8 mm/a)。华北地块中,天津沿海区域在垂向上具有较大的下沉速率,鄂尔多斯块体具有较大的上升速率。东北块体南部鄂尔多斯块体北方有微量上升,但XIAA(西安)站出现下降,垂向下降速率为4.24 mm/a,该站点异常可能与市政府建设有关[2]。本文得到的3个垂向沉降速率最大站点与占伟[4]结果一致,分别为HECX(河北沧州)、TJBH(天津滨海)、TJWQ(天津武清),其沉降速率分别为-28.56 mm/a、-14.25 mm/a、-38.38 mm/a。华南沿海地区多地呈轻微下沉趋势,青藏地块和华南地块交界处有部分区域微量隆升。全国垂向线性变化呈区域性分布特征,其中华北地块特征明显,西部轻微上升、东部轻微下沉,垂向变化速率范围波动在0.5 mm/a内。
周期项拟合后获得260个站点的周期项振幅,利用克里金插值法得到全国周期项振幅的分布特征。
3.2.1 年周期振幅分布特征由图 8可见,约90%站点的周年振幅大于2 mm,约96%站点的周年振幅小于10 mm,60% 站点的周年振幅小于5 mm,大多数站点的振幅为2~8 mm。在260个站点中,最大振幅为12.40 mm(YNLA),最小振幅为0.89 mm(FJLY),平均振幅为4.84 mm。中国GPS站年周期振幅的空间分布特征为东西较小、中部较大。川滇地区的垂向周期性运动最显著,振幅为6~14 mm,研究表明,降水、地下水负荷变化是影响该地区垂向周期运动变化的主要因素[9];华南地块年周期振幅由沿海向内陆逐渐递增,振幅为2~6 mm,东部沿海区域海潮负荷的影响大于内陆区域,这也是造成沿海与内陆年周期运动差异的因素之一;华北地块年周期振幅由中心向四周逐渐递减,振幅为3~6 mm;东北地块年周期振幅随纬度的增加而增大,振幅为0.89~5.40 mm;西域地块天山地区年周期振幅为4.68~7.13 mm,除天山区域外西域块体其他区域年周期季节特征不明显,年周期振幅为2.74~3.89 mm;青藏地块年周期振幅值随纬度的增加而减小,年周期振幅跨度较大,振幅为2.03~10 mm。
由图 9可见,95%站点的年周期项最大值出现在4~7月份,60%站点的年周期项最大值出现在5~6月份。随着纬度的增加,周期项最大值出现的时间逐渐延后,年周期最大值多出现在上半年,只有西域地块最北区域年周期项最大值出现在8~9月。全国南部地区年周期项最大值出现在2~3月,东北区域年周期项最大值出现在4~6月。沿海区域和内陆区域的年周期振幅峰值出现时间具有不同的分布特征,同一纬度上沿海区域的年周期项最大值出现时间会比内陆区域稍晚一些。沿海区域温度、降水量可能是影响沿海区域和内陆区域年周期振幅最大值出现时间差异的因素之一。
由图 10可见,超过80%站点的半年周期振幅小于2 mm,约33%站点的半年周期振幅小于1 mm,12%站点的半年周期振幅小于0.5 mm。在半年周期振荡中,最大振幅为5.056 mm(YNGM),最小振幅为0.069 mm(HLWD),平均振幅为1.380 mm。中国GPS站半年周期项振幅的空间分布特征为南北较小、中部较大。川滇地区半年周期项特征明显,振幅为1.67~5.06 mm,其他区域半年周期项振幅均小于3 mm;华南地区半年周期项振幅值由沿海向内陆递增,与该区域年周期项振幅分布特征相似,其半年周期项振幅为0.5~1.5 mm;华北地块半年周期项振幅为1.4~2.1 mm,中部振幅值较大、四周较小;东南沿海区域半年周期项振幅较小,振幅为0.51~1.16 mm。
由图 11可见,30%站点的半年周期项最大值出现在1~2月,58%站点的半年周期项最大值出现在5~6月。从全国范围来看,青藏和川滇地块半年周期项最大值多出现在5~6月;随着纬度的增加,西域块体半年周期项最大值出现的时间逐渐提前,并且出现时间存在梯度变化;中部地区半年周期项最大值的出现时间由中心向四周递增,大致在1~3月;东南沿海区域半年周期项最大值出现时间较早,大致在1~2月;东北区域半年周期项最大值出现较晚,大致在5~6月。综合分析可知,年周期项最大值与半年周期项最大值出现的月份空间分布特征截然不同。
由表 1可见,半年周期振幅平均值约为年周期振幅平均值的28.51%。垂向年周期振幅为4~10 mm,半年周期振幅为0.5~3.5 mm,说明年周期项为主要季节特征。结果表明,中国大陆连续GPS观测时间序列垂向周期运动具有明显的区域差异性特征。
通过前文获得各站的周期项振幅信号,年周期振幅和半年周期振幅在100°~110°E、30°~40°N范围内有聚集特征(图 12)。傅彦博等[7]研究全球季节项发现,北半球周年项振幅和纬度之间存在二项式拟合关系。基于这一结论,结合经纬度建立季节项振幅模型:
$ \begin{gathered} Y=a x+b y+c x^2+ \\ d x y+e y^2+f \end{gathered} $ | (12) |
式中,x为纬度,y为经度,a、b、c、d、e、f为待求参数。通过最小二乘法拟合得到拟合系数如表 2、3所示。
采用误差平方和(SSE)和均方根误差(RMSE) 来衡量模型精度,SSE的表达式为:
$ \mathrm{SSE}=\sum\limits_{i=1}^m\left(y_i-\tilde{y}_i\right)^2 $ | (9) |
式中,yi为模型值,
根据分块对各区域进行模型拟合,其拟合系数如表 2、3所示。在分区域模型中年周期模型、半年周期模型的SSE区域平均值分别为68.87 mm2、14.17 mm2,RMSE区域平均值分别为1.332 3 mm、0.532 1 mm,相比于全国模型,其拟合效果更好。对比全国模型与分区域模型可见(图 14),年周期模型、半年周期模型的分区域模型的振幅差值更接近于0,分布更加聚集,相较于全国模型有一定程度的提高。
表 4(单位%)中分区域模型的RMSE值大多低于全国模型,说明分区域模型拟合效果优于全国模型。其中,华北块体年周期模型振幅拟合效果提高最明显,相较于全国模型提高了48.76百分点,华南块体和川滇块体的半年周期模型拟合效果低于全国模型。对各区域拟合效果取平均值,年周期项区域模型较全国模型提高31.59百分点,半年周期项区域模型较全国模型提高9.41百分点。
1) 时间序列分析中得到的线性趋势项能够表示中国大陆垂向运动的大致趋势,西域块体天山地区、青藏高原南部、鄂尔多斯块体等地有明显上升趋势,天津沿海地区、华南沿海地区等地有下沉趋势。
2) 通过对选取的260个基准站的周期项振幅进行对比分析可知,年周期振幅大于半年周期振幅,前者通常是后者的3倍以上,年周期运动和半年周期运动具有不同的空间分布特征。季节项中川滇地区的年周期运动和半年周期运动最为明显。从整体上来看,年周期项最大值出现月份在沿海地区和内陆地区具有一定的差异性,相同纬度条件下沿海地区的年周期项最大值出现月份晚于内陆。
3) 年周期模型、半年周期模型相较于全国模型分别提高31.59百分点、9.41百分点,但二项式拟合方法容易受到异常值的影响,今后可对模型进行改进。通过小波分析可知,GPS时间序列的季节信号具有调制特征,没有被准确拟合的季节信号会进入到残差序列中。
[1] |
张诗玉, 钟敏, 唐诗华. 我国GPS基准站地壳垂直形变的大气负荷效应[J]. 武汉大学学报: 信息科学版, 2006, 31(12): 1090-1093 (Zhang Shiyu, Zhong Min, Tang Shihua. Vertical Crustal Displacements due to Atmospheric Loading Effects at GPS Fiducial Stations in China[J]. Geomatics and Information Science of Wuhan University, 2006, 31(12): 1090-1093)
(0) |
[2] |
王敏. GPS观测结果的精化分析与中国大陆现今地壳形变场研究[D]. 北京: 中国地震局地质研究所, 2009 (Wang Min. Analysis of GPS Data with High Precision and Study on Present-Day Crustal Deformation in China[D]. Beijing: Institute of Geology, CEA, 2009)
(0) |
[3] |
王伟. 中国大陆现今地壳运动的GPS分析与构造变形模拟[D]. 北京: 中国地震局地球物理研究所, 2013 (Wang Wei. Study on Present-Day Crustal Movement with GPS Data and Modeling the Active Deformation in the Chinese Mainland[D]. Beijing: Institute of Geophysics, CEA, 2013)
(0) |
[4] |
符养. 中国大陆现今地壳形变与GPS坐标时间序列分析[D]. 上海: 中国科学院上海天文台, 2002 (Fu Yang. Present-Day Crustal Deformation in China and GPS-Derived Coordinate Time Series Analysis[D]. Shanghai: Shanghai Astronomical Observatory, CAS, 2002)
(0) |
[5] |
畅柳, 杨国华, 朱爽. 中国大陆导航连续站周期信号空间规律研究[J]. 测绘科学, 2016, 41(1): 50-54 (Chang Liu, Yang Guohua, Zhu Shuang. Study on Spatial Distribution of GNSS Seasonal Signals for Continuous GNSS Stations in China[J]. Science of Surveying and Mapping, 2016, 41(1): 50-54)
(0) |
[6] |
占伟. 基于GPS连续观测的中国大陆典型区域地壳垂直运动研究[D]. 武汉: 武汉大学, 2017 (Zhan Wei. Study on Vertical Crustal Motion in Chinese Mainland and Typical Areas Based on Continuous GPS[D]. Wuhan: Wuhan University, 2017)
(0) |
[7] |
傅彦博, 孙付平, 朱新慧, 等. 全球GPS测站垂向周年变化统计改正模型的建立[J]. 测绘学报, 2018, 47(10): 1337-1345 (Fu Yanbo, Sun Fuping, Zhu Xinhui, et al. Establishment of Statistical Correction Model for Vertical Annual Variations of Global GPS Stations[J]. Acta Geodaetica et Cartographica Sinica, 2018, 47(10): 1337-1345)
(0) |
[8] |
范玉磊, 王贺, 黄声享, 等. 利用小波分析IGS跟踪站的非线性运动特征[J]. 测绘工程, 2014, 23(9): 5-8 (Fan Yulei, Wang He, Huang Shengxiang, et al. Wavelet Approach to Studying the Nonlinear Motional Characteristics of IGS Stations[J]. Engineering of Surveying and Mapping, 2014, 23(9): 5-8)
(0) |
[9] |
盛传贞, 甘卫军, 梁诗明, 等. 滇西地区GPS时间序列中陆地水载荷形变干扰的GRACE分辨与剔除[J]. 地球物理学报, 2014, 57(1): 42-52 (Sheng Chuanzhen, Gan Weijun, Liang Shiming, et al. Identification and Elimination of Non-Tectonic Crustal Deformation Caused by Land Water from GPS Time Series in the Western Yunnan Province Based on GRACE Observations[J]. Chinese Journal of Geophysics, 2014, 57(1): 42-52)
(0) |
[10] |
张培震, 邓起东, 张国民, 等. 中国大陆的强震活动与活动地块[J]. 中国科学D辑: 地球科学, 2003, 33(增1): 12-20 (Zhang Peizhen, Deng Qidong, Ma Jin, et al. Strong Earthquake Activity and Active Blocks in Chinese Mainland[J]. Science in China Ser D Earth Sciences, 2003, 33(S1): 12-20)
(0) |
2. Nanjing Zhixing Map Information Technology Co Ltd, 699 Xuanwu Road, Nanjing 210023, China