Analysis of strain field considering the influence of colored noise in Qinghai area
-
摘要: 基于青海地区14个中国大陆构造环境监测网络(CMONOC)连续站10年的坐标时间序列数据，利用贝叶斯信息准则(BIC)确定各连续站的最优噪声模型，进而得出修正后的水平速度场，并在此基础上建立整体旋转与线性应变模型来分析青海地区的应变特征. 结果表明：青海地区CMONOC连续站坐标时间序列各方向的噪声特性存在较大差异，东(E)、北(N)、天顶(U)方向的最优噪声模型分别为 白噪声+幂律噪声(WN+PL)、白噪声+高斯-马尔科夫噪声(WN+GGM)和白噪声+闪烁噪声(WN+FN). 考虑有色噪声(CN)影响青海地区CMONOC连续站基于ITRF2014框架下的平均水平运动速率为39.45 mm/a，运动方向为88°57′58″NEE. 青海地区构造活动相对强烈的东北与西南部分别表现为挤压应变特征和拉张应变特征；从东北向西南，挤压应变逐渐减小，拉张应变逐渐增大，总体表现为挤压应变.
-
关键词:
- 全球卫星导航系统(GNSS) /
- 有色噪声(CN) /
- 时间序列 /
- 贝叶斯信息准则(BIC) /
- 速度场 /
- 应变场
Abstract: Based on the 10-years coordinate time series data of crustal movement observation network of 14 china (CMONOC) continuous stations in Qinghai, the optimal noise model of each continuous station was determined by using Bayesian information criterion (BIC), and then the modified horizontal velocity field was obtained. On this basis, the global rotation and linear strain model was established to analyze the strain characteristics in Qinghai. The results show that the noise characteristics of CMONOC continuous station in Qinghai are different in different directions. The optimal noise models in east (E), north (N) and up (U) directions are “white noise + power law noise (WN+PL)”, “white noise + Gaussian Markov noise (WN+GGM)” and “white noise + flicker noise (WN+FN)” respectively. After considering the influence of colored noise (CN), the average horizontal motion velocity of CMONOC continuous station in Qinghai based on the framework of ITRF2014 is 39.61 mm/a and the motion direction is 88°57'58"NEE. The northeastern and southwestern parts of the Qinghai region with relatively strong tectonic activities are characterized by compressive strain and tensile strain, respectively. From northeast to southwest, the compression strain gradually decreases and the tensile strain gradually increases, which is generally manifested as the compression strain. -
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坐标时间序列中所涉及到的绝大部分参数. 坐标时间序列表示为
$$ \begin{split}y\left({t}_{i}\right)=&a+b{t}_{i}+c\;\mathrm{sin}\left({2\text{π}}{t}_{i}\right)+d\;\mathrm{cos}\left({2\text{π}}{t}_{i}\right)+e\;\mathrm{sin}\left({4\text{π}}{t}_{i}\right)+ \\& f\;\mathrm{cos}\left({4\text{π}}{t}_{i}\right)+{\displaystyle \sum _{j=1}^{{n}_{j}}{g}_{j}H\left({t}_{i}-{T}_{j}\right)}+{v}_{i}.\end{split} $$ (1) 式中：ti为历元，以a为单位；a为初始位置；b为线性速率；c、d、e、f分别为周期项的系数；H为海维西特阶梯函数；gj为同震偏移；vi为残差.
1.2 BIC
Hector软件使用BIC来评价拟合结果的优良性，该软件能够对时间序列的线性趋势项、高阶多项式、季节项、其他周期项信号以及多种噪声模型的组合进行估计[10]. 使用BIC分析噪声模型的公式为
$$ \ln \left( L \right) = - \frac{1}{2}\left[ {N\ln \left( {2{\text{π }}} \right) + \ln \det \left( {\boldsymbol{C}} \right) + {{\boldsymbol{r}}^{\text{T}}}{\bar{\boldsymbol{C}}^{ - 1}}{\boldsymbol{r}}} \right] . $$ (2) 式中：N为实际观测数，r为残差矩阵；协方差矩阵C可分解为
$$ {\boldsymbol{C}} = {\sigma ^2}\bar {\boldsymbol{C}} . $$ (3) 其中，
$ \sigma $ 为残差序列的标准差，其计算公式为$$ \sigma = \sqrt {\frac{{{{\boldsymbol{r}}^{\text{T}}}{{\bar {\boldsymbol{C}} }^{ - 1}}{\boldsymbol{r}}}}{N}} . $$ (4) 由
$\det cA = {c^N}\det A$ 可得$$ {\text{ln}}\left( L \right) = - \frac{1}{2}\left[ {N\ln \left( {2{{\text{π} }}} \right) + \ln \det \left( {\bar {\boldsymbol{C}} } \right) + 2N\ln \left( \sigma \right) + N} \right] . $$ (5) 于是可得
$$ {\text{BIC}} = k\ln \left( N \right) + 2\ln \left( L \right). $$ (6) 似然函数
$ L $ 越小，BIC值越小，噪声模型相对越好；反之，似然函数$ L $ 越大，BIC值越大，噪声模型相对越差.2. 时间序列分析
2.1 数据来源
CMONOC是中国地壳运动观测网的延续，在中国境内共设置了260个连续站，青海境内共有14个连续站，选取GAMIT/GLOBK 软件解算得到的青海地区2010—2020年的14个CMONOC连续站为基础数据，测站分布如图1所示. 本文数据来源于中国地震局GNSS数据产品服务平台[11].
此处仅展示QHGE、QHMY和QHTR站的原始坐标时间序列. 如图2所示，原始坐标时间序列的东（E）、北（N）、天顶（U）方向都存在异常值，必须予以剔除. E方向表现为线性递增趋势，斜率较大；N方向表现为线性趋势，斜率较小；U方向表现为周期变化趋势.
2.2 原始时间序列处理
本文利用最小二乘法估计并剔除坐标时间序列中的线性趋势和周期性趋势来计算残差，并使用四分位数粗差探测方法剔除异常值. 具体步骤包括：1) 采用最小二乘法剔除原始时间序列中的线性趋势和周期性趋势，得到时间序列中的残差；2) 用四分位距探测残差时间序列包含的异常值，然后剔除；3) 拟合剔除异常值之后的残差时间序列并重复步骤2)，直到完全剔除[9]. 图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方向噪声模型百分比% 方向 WN+GGM WN+PL WN+FN WN+RW+FN E 21.43 42.86 28.57 7.14 N 35.72 28.57 28.57 7.14 U 21.43 21.43 50.00 7.14 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.
图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]
$$ \begin{array}{l}\left[\begin{array}{c}{V}_{e}\\ {V}_{n}\end{array}\right]=\left[\begin{array}{ccc}\begin{array}{c}-r\;\mathrm{cos}\;\lambda\;\mathrm{sin}\;\varphi \\ r\;\mathrm{sin}\;\lambda \end{array}& \begin{array}{c}-r\;\mathrm{sin}\;\lambda \;\mathrm{sin}\;\varphi \\ -r\;\mathrm{cos}\;\lambda \end{array}& \begin{array}{c}r\;\mathrm{cos}\;\varphi \\ 0\end{array}\end{array}\right]\left[\begin{array}{c}{\omega }_{x}\\ {\omega }_{y}\\ {\omega }_{z}\end{array}\right]+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left[\begin{array}{cc}{A}_{0}& {B}_{0}\\ {B}_{0}& {C}_{0}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\displaystyle\frac{1}{2}\left[\begin{array}{cc}{A}_{1}& {B}_{2}\\ {B}_{1}& {C}_{2}\end{array}\right]\left[\begin{array}{c}{x}^{2}\\ {y}^{2}\end{array}\right]+\left[\begin{array}{cc}{A}_{2}& {B}_{1}\\ {B}_{2}& {C}_{1}\end{array}\right]xy.\end{array} $$ (7) 式中：
${A_0} \text{～} {A_2}$ 、${B_0} \text{～} {B_2}$ 、${C_0} \text{～} {C_2}$ 为板块的应变参数；$ {\omega _x} $ 、$ {\omega _y} $ 、$ {\omega _z} $ 为板块的欧拉矢量的3个分量；$ r $ 为地球半径；$ \lambda $ 为质点的大地经度；$ \varphi $ 为质点的大地纬度.4.2 应变场分析
本文利用经最优噪声模型改正后的速度场建立青海地区地壳水平运动的整体旋转与线性应变模型. 通过分析计算得到青海地区内部形变的空间分布特征并绘制该区域的主应变图、面膨胀图以及最大剪应变图. 由图6可知，青海地区东北部表现为明显的挤压应变特征，而该地区的西南部主要表现为拉张特征，最大主压应变和最大主张应变分别出现在青海地区的东北角与西南角，其值分别可达−4.31×10−8和4.56×10−8. 由图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）青海地区的东北、西南部构造活动较强烈，东北部表现为明显的挤压应变特征，而西南部主要表现为拉张特征. 从东北往西南，挤压应变逐渐减小，拉张应变逐渐增大，总体表现为挤压应变.
-
表 1 E、N、U方向噪声模型百分比
% 方向 WN+GGM WN+PL WN+FN WN+RW+FN E 21.43 42.86 28.57 7.14 N 35.72 28.57 28.57 7.14 U 21.43 21.43 50.00 7.14 
下载: 导出CSV
-
[1] 瞿伟, 高源, 陈海禄, 等. 利用GPS高精度监测数据开展青藏高原现今地壳运动与形变特征研究进展[J]. 地球科学与环境学报, 2021, 43(1): 182-204. [2] 黄立人, 符养. GPS连续观测站的噪声分析[J]. 地震学报, 2007, 29(2): 197-202. DOI: 10.3321/j.issn:0253-3782.2007.02.009 [3] MAO A, HARRISON C G A, DIXON T H. Noise in GPS coordinate time series[J]. Journal of geophysical research atmospheres, 1999, 104(B2): 2797-2816. DOI: 10.1029/1998JB900033
[4] WILLIAMS S D P, BOCK Y, FANG P, et al. Error analysis of continuous GPS position time series[J]. Journal of geophysical research solid earth, 2004, 109(B3): B03412. DOI: 10.1029/2003JB002741
[5] ZHANG J, BOCK Y, JOHNSON H, et al. Southern california permanent GPS geodetic array: error analysis of daily position estimates and site velocities[J]. Journal of geophysical research solid earth, 1997, 102(B8): 18035-18055. DOI: 10.1029/97JB01380
[6] 蒋志浩, 张鹏, 秘金钟, 等. 顾及有色噪声影响的CGCS2000下我国CORS站速度估计[J]. 测绘学报, 2010, 39(4): 355-363. [7] 管雅慧, 王坦, 胡顺强, 等. 云南地区GPS基准站噪声模型分析[J]. 全球定位系统, 2020, 45(2): 68-73,79. [8] 王伟, 王迪晋, 陈正松, 等. 用GPS资料分析青藏高原现今应变率场[J]. 大地测量与地球动力学, 2017, 37(9): 881-883,897. [9] BOS M S, FERNANDES R M S, WILLIAMS S D P, et al. Fast error analysis of continuous GNSS observations with missing data[J]. Journal of geodesy, 2013, 87(4): 351-360. DOI: 10.1007/s00190-012-0605-0
[10] 何月帆. 高精度GPS数据处理及时间序列分析[D]. 西安: 长安大学, 2019. [11] 中国地震局GNSS数据产品服务平台 [DB/OL]. (2017-09-06)[2021-05-01]. http://www.eqhb.gov.cn/info/1435/12636.htm [12] 吴啸龙, 向洋, 汤伏全. 基于GPS应变与震源机制解应力反演喜马拉雅构造带现今地壳形变特征[J]. 地球物理学报, 2020, 63(8): 2924-2939. DOI: 10.6038/cjg2020N0362 [13] 瞿伟, 张勤, 张冬菊, 等. 基于GPS速度场采用RELSM模型分析青藏块体东北缘的形变—应变特征[J]. 地球物理学进展, 2009, 24(1): 67-74. [14] 雷前坤. 环渤海区域地壳水平运动形变模型构建及特征分析[D]. 贵阳: 贵州大学, 2020. [15] 张俊, 李屹旭, 张显云. 利用半参数模型精化REHSM和RELSM两种弹性地壳运动模型[J]. 大地测量与地球动力学, 2018, 38(4): 362-365. [16] 李延兴, 张静华, 何建坤, 等. 菲律宾海板块的整体旋转线性应变模型与板内形变-应变场[J]. 地球物理学报, 2006, 49(5): 1339-1346. DOI: 10.3321/j.issn:0001-5733.2006.05.013
计量
- 文章访问数: 382
- HTML全文浏览量: 86
- PDF下载量: 31
