山东GNSS地壳运动监测网自2007年建立至今，已积累了大量的观测数据，成为获取地壳形变与地球动力学信息的最佳数据源之一。GNSS基准站由于板块运动、海洋潮汐、温度变化、地下水升降等因素的影响，会在水平和高程方向产生位移^[1]。通过时间序列分析，可获取测站精确的运动趋势，进而反映地壳形变特征^[2]。

国内外学者对GNSS基准站的垂直运动特征进行了大量研究。张飞鹏等^[3]分析GNSS基准网1999~2001年的高程时间序列，探测到3~10 mm的季节性变化，然而由于观测跨度较短，未给出垂直速率；符养^[4]利用200多个基准站的非线性结果，发现整体地球高程方向存在周期性振荡运动；顾国华^[5]分析基准站垂直形变的年周期变化特点，讨论获得中国大陆垂直形变长趋势速率的条件，表明由基准站多期GPS观测得到的长趋势垂直形变速率较为可靠；杨强等^[6]利用ARMA模型对部分IGS跟踪站GPS高程时间序列建模，改进GPS高程时间序列；梁诗明^[7]采用地球物理模型与数学拟合模型相结合的方法，有效提高非连续GPS观测站的垂直速率估计精度。

上述研究均利用不同时段的GPS观测资料，主要研究基准站垂直形变的周期性情况，但还存在一些局限性，主要体现在2个方面：一是所采用的数据观测时间较短，不能体现区域的长时间垂直形变，垂直精度也存在一定的局限性；二是未利用大观测跨度的基准站时间序列数据拟合出垂直形变场，无法对区域垂直形变进行综合分析。本文采用山东地区38个GNSS基准站2011~2016年的连续观测数据进行时间序列分析，给出垂直形变速率和垂直形变场，为区域垂直形变特征研究及地球动力学分析提供依据。

1 研究区资料概况

山东地区位于中国东部沿海，黄河下游，有多种地质构造。山东东部只有古老的结晶岩基底，属单层结构型；西部地区在古老的结晶岩基底之上，还存在着较新的沉积岩层覆盖，属双层结构上的地质构造主要有褶皱构造和断裂构造。山东地区的山地地势的形成，明显受地质基础的控制。平原地区主要分布在山东省的西部和北部，属于断块构造差异沉降区；山地丘陵地区分布于山东的中南和东部地区，属于断块构造差异隆起区^[8]。

为监测山东地区地壳运动变化，研究其主要板块运动及断裂带的活动状况，山东省自2007年开始建立GNSS基准站观测网络^[9]，现有119个GNSS连续观测基准站。为获取山东地区垂直形变场的分布特征，筛选出数据完整度超过90%的站点，选取2011-01-01~2016-12-31共6 a的时间序列数据资料，最终采用38个均匀分布且时间序列数据资料较为可靠的GNSS基准站(图 1)。对比分析现有资料，给出山东地区近期垂直形变分布特征。

2 数据处理方法

本文通过以下方法对时间序列数据资料进行预处理：1)以2011-01-01为时间起点将时间序列数据统一到同一时间基准，计算分析6 a内站点的垂直形变；2)以3倍中误差为阈值，剔除由于外界条件突然的变化(更换基准站天线或地震)产生的阶跃项数据，减小粗差对实验结果造成的偏差；3)GNSS基准站由于受停电或其他因素的影响导致坐标时间序列部分数据缺失、采样不均，选取三次多项式插值法进行处理，以获得连续、均匀采样的时间序列数据；4)采用趋势面插值法对山东地区沉降量进行插值，获得山东地区垂直形变场。

2.1 年平均速率计算

利用预处理得到的时间序列数据计算GNSS基准站的年垂直位移和年均垂直矢量位移速率。GNSS基准站的年垂直位移为：

$ \begin{gathered} {H_i}{\text{ = }}\frac{1}{m}\sum\limits_{j = 1}^m {\left( {{H_{ij}} - {H_{ij - 1}}} \right)} \hfill \\ i = 1, 2 \cdots n;j = 1, 2 \cdots , m \hfill \\ \end{gathered} $

(1)

式中，H_ij为第i年第j个观测值，m为每年的观测值总数，n为年时间序列总数，H_i为每个基准站的年垂直形变均值。

GNSS基准站垂直位移为：

$ \Delta H = \sum\limits_{i = 1}^n {{H_i}} $

(2)

式中，ΔH指每个基准站的垂直位移。

GNSS基准站年均垂直矢量位移速率为：

$ \Delta V = \frac{1}{n}\Delta H $

(3)

式中，ΔV指每一基准站的年矢量位移速率。ΔV为正，表示基准站上升；ΔV为负，表示基准站沉降。

2.2 趋势面插值法

趋势面插值用多项式表示的面按最小二乘法原理对数据点进行拟合，并用于估算其他值的点。假设二维空间中有n个观测点(x_l, y_l)(l=1, 2, …, n)，观测值为z_l(l=1, 2, …, n)，则空间分布z的趋势面可表示为N次多项式：

$ \hat z = \sum\limits_{i, j = 0}^N {{a_{ij}}{x^i}{y^j}} , i + j \leqslant N $

(4)

式中，a_ij为常数项。

2.3 趋势面拟合适度检验

利用趋势面拟合适度R²和显著性F分布对趋势面插值结果进行检验。R²检验公式为：

$ {R^2} = \frac{{{\text{S}}{{\text{S}}_R}}}{{{\text{S}}{{\text{S}}_T}}} = 1 - \frac{{{\text{S}}{{\text{S}}_D}}}{{{\text{S}}{{\text{S}}_T}}} $

(5)

式中，SS_T=SS_D+SS_R，其中$ {\text{S}}{{\text{S}}_D} = \sum\limits_{i = 1}^n {{{({z_i} - {{\hat z}_i})}^2}} $为剩余平方和，表示随机因素对z的离差；$ {\text{S}}{{\text{S}}_R} = \sum\limits_{i = 1}^n {{{({{\hat z}_i} - {{\bar z}D_i})}^2}} $为回归平方和，它表示p个自变量对因变量z的离差的总影响。R²越接近1，趋势面的拟合度就越高。

趋势面拟合适度的显著性F检验公式为：

$ F = \frac{{{\text{S}}{{\text{S}}_R}/p}}{{{\text{S}}{{\text{S}}_D}/(n - p - 1)}} $

(6)

式中，p为多项式项数(不包括常数项)。在显著性水平下，若计算值F大于临界值F_α，则认为趋势面方程显著，否则不显著。

3 数据分析

以2011-01-01数据为基准，绘制2011~2016年的垂直位移时间序列(图 2)。为获得站点的长期趋势项，利用线性直线拟合法对时间序列数据进行最佳直线拟合，剔除基准站垂直形变周期项的影响。由图 2看出，GNSS基准站有明显的周期性变化，表明基准站周边应变状态较为平稳。拟合直线与时间序列数据关联度较大，说明可以利用拟合直线表示基准站的垂直变化情况。从时间和空间尺度上看，基准站垂直位移随着时间平缓地上升或下降。GNSS基准站的垂直位移也可以反映出山东地区地壳的垂直形变，进而了解山东地区地壳的活动情况。

表 1给出了6 a内山东地区GNSS基准站的年均垂直形变速率。可知，山东地区大部分GNSS基准站呈现沉降的趋势；山东半岛地区及南部地区基准站沉降速率平缓，属于正常沉降，速率为1.0~8.0 mm/a；山东地区北部及西部地区沉降量较大，可能是受过度开采地下水的影响。

图 3给出了GNSS基准站的垂直形变速率场。其中，山东地区沉降速率最大的地区为菏泽(HEZE)，速率为-24.2 mm/a；其次是聊城(LICH)，速率为-22.6 mm/a；沉降最小的地区为石岛(SHID)，仅为-0.6 mm/a。在山东地区沿海地区以及山东地区南部地区的GNSS基准站沉降速率较慢，为1~7 mm/a；而位于山东地区中部地区的GNSS基准站有小幅度的隆升，上升速率约为0.5~3.0 mm/a。通过以上分析可知，GNSS基准站的垂直形变速率场能反映山东地区整体上的垂直变化规律。

对38个GNSS基准站时间序列数据处理后，利用趋势面插值方法对时间序列数据进行加密处理，获取山东地区垂直形变场(图 4)。可见，山东地区整体上呈现差异性的下降趋势，以沂沭断裂带(郯庐断裂带山东段)为界线，山东地区的垂直形变构造分区也较为显著。山东东部地区垂直形变较为平缓，总体表现为下降趋势，但其中也存在着明显的差异运动；在沂沭断裂带西部地区垂直形变更为明显，在山东西部地区呈现沉降趋势，而在中部山区出现上升现象，这与地质构造密切相关。山东地区的垂直形变场与地形地貌大体一致，一定程度上表明垂直形变场能反映出地壳垂直差异性运动趋势。

自滨州(SDBZ)向西南方向延伸，经过聊城(LICH)和菏泽(HEZE)存在一个沉降带，近6 a沉降量由15.0 mm逐渐增大至140.0 mm。该沉降带内存在2个较大的沉降区，其中西北部地区主要分布在滨州(SDBZ)、阳信(YAXI)以及垦利(SDKL)周围，该沉降区自东向西呈带状分布，年沉降速率为10.0~15.0 mm/a，平均沉降速率为12.1 mm/a；西部地区的沉降区分布在聊城(LICH)和菏泽(HEZE)周围，沉降区由北向南呈带状分布，年沉降速率为10.0~24.0 mm/a，平均沉降速率为16.7 mm/a。这两个沉降区走向与此地区构造断裂带走向一致，其沉降原因与过度开采地下水有关，也可能受断裂带运动的影响，另文再详细描述。在山东东部沿海以及南部地区沉降量较小，年沉降速率为0~6.0 mm/a，平均沉降速率为3.2 mm/a。在沂沭断裂带及周围较为稳定，没有明显的异常沉降。在山东中部泰山山区由临沂(SDLY)沿东北方向至泰安(TAIN)有明显的隆升趋势，年隆升速率为0.5~3.0 mm/a，平均上升速率为1.3 mm/a。

对趋势面拟合适度进行R²检验结果表明，R²检验为0.989，趋势回归模型的显著性较高，拟合程度较好。对趋势面适度的显著性F检验，在置信水平α=0.05下查F分布表得F_α=F_0.05(9, 28)=2.24，显然，F=297.4 > F_α，故趋势面的回归方程显著。因此，利用趋势面插值所得的垂直形变量具有一定的可信度。

4 结语

选取山东地区GNSS基准站2011-01-01~2016-12-31共6 a、38个站的时间序列数据，获取了山东地区高精度、高时空分辨率的地壳垂直形变速率。分析认为：

1) 山东大部分地区呈现差异性下降状态，只有中部泰山山脉地区出现上升现象，与平原洼地下降、山地丘陵地区上升的继承性构造活动相一致。

2) GNSS基准站时间序列有明显的周期性变化，表明基准站较为稳定。山东半岛区域的基准站呈现正常的沉降趋势，而山东省的西北部地区和西部地区由于过度开采地下水，沉降较大。

3) 山东地区以沂沭断裂带为界，构造分区较为显著，沂沭带东部地区垂直形变平缓，年均沉降速率为-3.2 mm/a；西部地区垂直形变特征较为明显，存在两个呈带状分布的沉降区，形成原因与过度开采地下水有关。其中，分布在滨州(SDBZ)、阳信(YAXI)以及垦利(SDKL)周围的沉降区范围较小，年均沉降速率为-12.1 mm/a；分布在聊城(LICH)以及菏泽(HEZE)周围的沉降区范围较大，年均沉降速率为-16.7 mm/a。而在山东中部地区(泰山山脉)出现抬升现象，年均上升速率为1.3 mm/a。