顾及有色噪声影响的青海地区应变场分析

张幸，魏冠军

（ 1. 兰州交通大学 测绘与地理信息学院, 兰州 730070；2. 地理国情监测技术应用国家地方联合工程研究中心, 兰州 730070；3. 甘肃省地理国情监测工程实验室, 兰州 730070 ）

0 引 言

目前，中国大陆构造环境监测网络(CMONOC)积累了大量连续运行参考站(CORS)的原始观测数据，为研究青海地区的板块运动、地壳形变、速度场模型等提供了十分重要的基础数据[1]. 由于观测数据中存在的有色噪声(CN)降低了测站运动参数的估计精度，且不同地区的噪声特性差异明显[2-7]，因此，选择合适的噪声模型来研究连续站的稳定性及地表规律具有重要意义. ZHANG等[5]在加利福尼亚南部地区时间序列研究中引入噪声分析，发现时间序列噪声具有白噪声+闪烁噪声(WN+FN)的特征. 蒋志浩等[6]的研究表明，CORS站坐标时间序列具有WN、FN和随机漫步噪声(RW)的噪声特征. 管雅慧等[7]对云南地区CMONOC连续站坐标时间序列进行分析，认为时间序列观测数据超过5年时，噪声模型才能逐渐趋于稳定. 王伟等[8]基于CMONOC观测数据获得了青藏高原的水平运动速度场及应变场，认为该地区应变场的整体分布特征与该地区活动构造的空间分布及地震活动十分一致. 综合上述研究表明，应变场分析应当选择较长时间段的观测数据并顾及CN的影响.

针对目前青海地区应变场研究中未顾及CN的影响且未对较长时间段(10年或以上)坐标时间序列噪声进行分析，本文通过分析青海地区14个CMONOC连续站2010—2020年的时间序列数据，分别确定了3个方向的最优噪声模型，得到国际地球参考框架(ITRF)下经最优噪声模型改正后的青海速度场，并基于修正后的速度场建立整体旋转与线性应变模型来分析青海地区的应变特征.

1 时间序列噪声分析方法

目前大部分测绘工作者主要使用时间序列分析软件CATS和HECToR等软件对不同区域的GPS观测数据进行时间序列噪声分析. 其中，CATS软件使用最大似然估计法(MLE)且具有较为精确的算法和模型. 但是，在处理时间跨度较大并且噪声类型较为复杂的时间序列数据时，CATS软件的处理速度十分缓慢. 为更高效处理含有大量数据的时间序列，BOS等[9]研发了HECToR软件，该软件通过算法改进提高了数据处理的速度，很好地克服了CATS软件的弊端. 改进后的MLE即为贝叶斯信息准则(BIC).

1.1 MLE

目前坐标时间序列分析普遍使用MLE，该方法可以得到GPS坐标时间序列中所涉及到的绝大部分参数. 坐标时间序列表示为

式中：ti为历元，以a为单位；a为初始位置；b为线性速率；c、d、e、f分别为周期项的系数；H为海维西特阶梯函数；gj为同震偏移；vi为残差.

1.2 BIC

Hector软件使用BIC来评价拟合结果的优良性，该软件能够对时间序列的线性趋势项、高阶多项式、季节项、其他周期项信号以及多种噪声模型的组合进行估计[10]. 使用BIC分析噪声模型的公式为

式中：N为实际观测数，r为残差矩阵；协方差矩阵C可分解为

其中， σ 为残差序列的标准差，其计算公式为

detcA=cNdetA

由 可得

于是可得

似然函数L越小，BIC值越小，噪声模型相对越好；反之，似然函数L越大，BIC值越大，噪声模型相对越差.

2 时间序列分析

2.1 数据来源

CMONOC是中国地壳运动观测网的延续，在中国境内共设置了260个连续站，青海境内共有14个连续站，选取GAMIT/GLOBK 软件解算得到的青海地区2010—2020年的14个CMONOC连续站为基础数据，测站分布如图1所示. 本文数据来源于中国地震局GNSS数据产品服务平台[11].

图1 青海地区CMONOC站点分布图

此处仅展示QHGE、QHMY和QHTR站的原始坐标时间序列. 如图2所示，原始坐标时间序列的东(E)、北(N)、天顶(U)方向都存在异常值，必须予以剔除. E方向表现为线性递增趋势，斜率较大；N方向表现为线性趋势，斜率较小；U方向表现为周期变化趋势.

图2 测站原始时间序列图

2.2 原始时间序列处理

本文利用最小二乘法估计并剔除坐标时间序列中的线性趋势和周期性趋势来计算残差，并使用四分位数粗差探测方法剔除异常值. 具体步骤包括：1) 采用最小二乘法剔除原始时间序列中的线性趋势和周期性趋势，得到时间序列中的残差；2) 用四分位距探测残差时间序列包含的异常值，然后剔除；3) 拟合剔除异常值之后的残差时间序列并重复步骤2)，直到完全剔除[9]. 图3为处理后的残差时间序列.

图3 处理后测站残差时间序列图

2.3 最优噪声模型的确立

对14个CMONOC连续站E、N、U方向的时间序列数据进行解算，并将WN与FN、PL、GGM及FN+RW进行组合，利用这四种组合模型计算每个测站所对应E、N、U方向的BIC值，对比各个模型所得BIC值即可获得所有站点各方向的最优噪声模型.表1为各方向噪声模型所占百分比，由表1可知，各方向具有不同噪声特征. 在E、N、U方向上，最优噪声模型分别为 WN+PL、WN+GGM和WN+FN.

表1 E、N、U方向噪声模型百分比 %

3 顾及有色噪声的青海地区水平速度场分析

图4为青海境内未经最优噪声模型修正的CMONOC连续站的水平运动速度场. 由图4可知，青海地区的整体运动方向近E方向，东部地区的运动方向近E方向，其余地区的运动方向为E、N方向.未经最优噪声模型修正的青海CMONOC连续站在ITRF14框架下水平运动方向为88°55′10″NEE，平均速度为39.58 mm/a. E方向的平均速度为39.57 mm/a，标准差为5.09 mm；N方向的平均速度为0.75 mm/a，标准差为3.76 mm.

图4 修正前的水平速度场

图5为经最优噪声模型改正后的青海CMONOC连续站的水平运动速度场，修正后的青海CMONOC连续站在ITRF14框架下水平运动方向为88°57′58″NEE，平均速率为39.45 mm/a. 并且，E方向的平均速度为39.44 mm/a，标准差为5.01 mm；N方向的平均速度为0.71 mm/a，标准差为3.76 mm. 将最优噪声模型修正前后的速度场进行对比分析可知，修正后速度场的精度优于修正前速度场的精度，其平均水平运动速率相差0.13 mm/a；水平运动方向相差2′48″NEE.

4 顾及有色噪声的青海地区应变特征分析

4.1 应变分析模型

以往研究表明，地壳板块是弹塑性体或黏弹性体，在周围板块的作用下会产生整体平移、旋转和内部形变[12-13]. 块体内部应变的实质是板块内部质点的运动，并且这些质点运动的应变参数往往是与位置相关的线性函数[14-15]. 顾及板块的整体运动与线性形变可得块体的整体旋转与线性应变模型，其形式为[16]

式中：A0～A2、B0～B2、C0～C2为板块的应变参数； ωx、ωy、ωz为板块的欧拉矢量的3个分量；r为地球半径；λ为质点的大地经度；φ为质点的大地纬度.

4.2 应变场分析

本文利用经最优噪声模型改正后的速度场建立青海地区地壳水平运动的整体旋转与线性应变模型.通过分析计算得到青海地区内部形变的空间分布特征并绘制该区域的主应变图、面膨胀图以及最大剪应变图. 由图6可知，青海地区东北部表现为明显的挤压应变特征，而该地区的西南部主要表现为拉张特征，最大主压应变和最大主张应变分别出现在青海地区的东北角与西南角，其值分别可达-4.31×10-8和4.56×10-8. 由图7可知，青海地区的东北、西南部为最大剪应变高值区，对应构造活动较强烈. 面膨胀率反映挤压与拉张作用的相对大小，正、负值分别表示以张性活动或压性活动为主. 由图8可知，从东北向西南，挤压应变逐渐减小，拉张应变逐渐增大，总体表现为挤压应变.

图6 主应变图

图7 最大剪应变图

图8 面膨胀图

5 结 论

本文基于最优噪声模型，对青海境内14个CMONOC连续站2010—2020年期间的观测数据进行分析. 首先，利用最优噪声模型得到修正后的青海速度场. 其次，使用整体旋转与线性应变模型对速度场进行插值，并求出应变参数. 最后，根据所得结果对青海地区的应变特征进行分析，结论如下：

1)青海地区CMONOC连续站坐标时间序列在E、N、U方向上具有不同的噪声特性，其最优噪声模型分别为组合模型PL+WN、GGM+WN和 FN+WN；

2)考虑CN影响后的速度场精度明显优于未经修正的速度场精度，修正前后的水平运动方向和平均水平运动速率差值可达2′48″NEE和0.13 mm/a，修正后青海地区CMONOC连续站基于ITRF2014框架下水平方向运动的平均速率为39.45 mm/a，运动方向为88°57′58″NEE；

3)青海地区的东北、西南部构造活动较强烈，东北部表现为明显的挤压应变特征，而西南部主要表现为拉张特征. 从东北往西南，挤压应变逐渐减小，拉张应变逐渐增大，总体表现为挤压应变.