东北陆态网坐标时序有色噪声及环境负载分析

黄裕，隋哲民，金泽林

(1.辽宁工程技术大学 测绘与地理科学学院,辽宁 阜新 123000；2.兰州交通大学 测绘与地理信息学院,兰州 730070；3.中国能源建设集团辽宁电力勘测设计院有限公司,沈阳110000)

0 引言

中国大陆环境监测网络(CMONOC)在全国范围内已陆续建成260 个连续运行参考站(CORS) 目前已经积累了数量可观且连续的观测数据[1].通过对观测数据进行处理，可以分析相应地区的速度场，进而为地壳运动形变规律研究以及地震分析预报等提供数据基础[1].谢方等[2]通过最小二乘配置移动拟合推估法以及克里金插值法于2015 年初步建立东北地块的速度场.于吉鹏等[3]通过对东北地区2011—2013 年的GAMIT 结果的解算初步获得了该地区在全球框架下的速度场.朱李忠等[4]选取2012—2017 年中国东北及俄罗斯远东地区东南部的GPS 观测数据初步得到了该地区基于ITRF2014 下的速度场.但是，由于东北地区陆态网络建立时间较晚，上述研究并未使用较长时间段即5 a 以上的时间序列，或者未对其有色噪声(CN)及环境负载进行分析研究.倘若忽视坐标时间序列数据中存在的CN，将无法有效分离CN 与观测的时间序列数据，进而会对速度场的解算产生很大的影响.管雅慧等[5]的研究结果也证实了在分析速度场时，需要积累超过5 a 的时间序列观测数据才能获得较为可靠且稳定的时间序列最优噪声模型.同时，姜卫平等[6]也指出在研究测站时间序列及其速度场时环境负载不可忽视，尤其是在沿海地区需要重视非潮汐海洋负载.其研究成果也证明，采用环境负载数据对坐标时间序列数据进行改正后，得到的非潮汐海洋负载位移对沿海区域测站的 GPS 会进一步优化[6].

为此，本文选取东北地区16 个陆态网络连续站2010—2021 年间 10 a 的时间序列观测数据及同一时段的环境负载数据进行研究，确定了三个方向分量的最优噪声模型类型，并剔除环境负载的影响，进而得到了东北地区陆态网络连续观测站经改正后基于国际地球参考框架(ITRF)下的速度场.为东北地区高精度坐标框架的研究提供了数据基础，为我国东北地区板块运动、地壳形变监测、矿区生产建设等提供了重要的参考资料.

1 时间序列获取

1.1 坐标时间序列的构建

坐标时间序列函数模型可表示为[7-8]

式中：时间ti的单位为年，y(ti) 为测站ti时刻下的位移；参数a为截距，b为时间序列中线性变化的速率；第j个谐波对应的振幅、角频率以及相位分别表示为Aj、ωj和 φj；n为谐波的数量；gk为阶跃项；Tk为产生阶跃的时刻；H为阶跃函数；参数gi为同震偏移(offset)；ε(ti)为随机过程[9].

本文选取GAMIT/GLOBK 软件[10]解算得到的基于ITRF14 框架下的东北地区16 个陆态网络连续站近10 a 的原始坐标时间序列数据，具体的数据解算流程可参考文献[11-12]，陆态网络连续观测站分布如图1 所示.

图1 东北地区CMONOC 测站分布图

鉴于相同测站在不同时段所求得最优噪声模型特性不同，进行白噪声(WN)分析时需选取10 a 以上的时间序列数据才能获得较为稳定的噪声模型[13].本文选取的连续站时间序列数据均来源于中国地震局GNSS 数据产品服务平台(http://www.cgps.ac.cn)，所选测站对应的时间段如表1 所示.

表1 东北地区CMONOC 站点概况

在评价GNSS 时间序列数据精度时，通常使用加权均方根误差(WRMS)来评价，WRMS 的值越小，表明对应时间序列数据的重复性更佳[14].对应站点北(N)、东(E)、天顶(U)方向的WRMS 值越小说明时间序列离散程度越小，重复性越好，精度越高.计算公式为

式中：nexp为预测值；ncale为实际观测值；σexp为预测值中误差.

其中，水平方向上的理想WRMS 值约为1～2 mm，垂直方向上的理想WRMS 值约为3～10 mm[14].东北地区陆态网络连续站时间序列在N、E、U 方向上的WRMS 值统计结果如表2 所示.由表2 可知，东北地区大部分站点原始时间序列均满足理想情况，但HLFY、HLMH 以及SUIY 站在N 方向上WRMS 值偏大，CHUN、HLMH、HRBN、JLCB、JLYJ 以及SUIY站在E 方向上WRMS 值偏大.

表2 东北连续站时间序列各分量WRMS 值 mm

1.2 残差时间序列获取

以LNDD 站为例，图2 为该站的残差时间序列图像.可以看出，N、E 方向分量观测数据均有明显的线性变化趋势，N 方向随时间变化呈现向南运动趋势且斜率较小，E 方向随时间变化呈现向东运动趋势且斜率较大；U 方向的变化趋势呈现出明显的周年项波动特性；原始时间序列观测数据中存在明显的奇异值(或称为外野值)，必须在数据处理前对其进行剔除.

图2 LNDD 站残差时间序列图像

Hector 软件使用最小二乘法估计原始坐标时间序列的线性趋势与周期性趋势来剔除时间序列奇异值，进而通过去除线性趋势和周期性趋势计算残差.在去趋势项后利用四分位数粗差探测方法剔除残差.其具体步骤为：1)基于最小二乘法，利用式(1)获取残差时间序列；2)基于残差时间序列，利用四分位距探测原始数据所包含的粗差，并剔除其中的奇异值；3)重新拟合新获取的残差时间序列并再次进行步骤2，直至完全剔除粗差[11].

2 噪声分析

目前，速度场研究大多采用CATS 软件来分析坐标时间序列噪声，但是在处理较大量的数据时，解算速度将十分缓慢，无法为数据解算提供便利[6].为了克服此难题，BOS 等[15]研究人员于2012 年开发了Hector 时间序列分析软件，该软件通过优化算法大幅提高了数据处理速率.这种经过算法优化后的极大似然估计法(MLE)即为贝叶斯信息准则(BIC)数值分析法.

2.1 BIC 数值分析法

Hector 软件可以用来估计线性趋势项、周期项信号以及不同组合噪声模型.在解算最优的噪声模型时，Hector 软件采用BIC 信息判别准则，该准则是评价统计模型拟合结果优良性的一种标准，采用对数似然函数来解算数据，并在解算过程中添加参数k来避免过度拟合[15].此对数函数的定义为

式中：N为实际观测数据的数量(不包括缺失的数据)；r为残差矩阵；协方差阵C可分解为

式中：C为各类噪声的总和；σ 为残差序列的标准差

由 d etCA=CNdetA，式(3)～(5)可得

通过上式可以粗略知晓所解算模型的复杂程度以及该模型拟合数据的准确程度[15]，具体可定义为

式中，k为模型参数的个数，似然函数L越大对应的模型更优，即具有(相对)最小BIC 值的模型为最优模型.在对不同组合模型进行解算后，比对这些模型的BIC 值，进而可以得到最优的组合噪声模型.

2.2 噪声模型分析

通过Hector 软件对16 个连续站三个方向分量的数据进行解算，分别使用这4 种组合模型计算出每个站所对应的N、E、U 三个方向的BIC 值，然后在同一分量中找到每个测站所对应的最小BIC 值.全部测站E 方向上的噪声模型解算所得BIC 差值，如图3所示.

图3 各测站E 方向BIC 差值

通过对不同组合噪声模型对应的BIC 值比较可知三个方向时间序列的不同组合噪声模型间的关系，且N、E、U 分量具有不同的噪声特性：在N 方向上，最优噪声模型为WN+FN，占比为56.25%；在E 方向，最优噪声模型为WN+PL，占比为56.25%；U 方向上，最优噪声模型为WN+FN，占比为87.5%.

3 顾及CN 的东北地区速度场分析

绘制经最优噪声模型改正后的ITRF14 框架下东北地区CMONOC 基准站的水平运动速度场，如图4 所示.东北地区整体呈现近SE 向运动，只考虑CN 修正后东北陆态网络基准站在ITRF14 框架下N 方向上运动的平均速率为-12.902 mm/a，E 方向上运动的平均速率为26.852 mm/a，U 方向上运动的平均速率为0.312 mm/a.

图4 修正后的水平速度场

绘制经最优噪声模型改正后的ITRF14 框架下东北地区水平运动速度场，如图4 所示.

修正后东北陆态网络基准站的水平速度精度较高，E 方向速度的标准差为2.462 mm，N 方向标准差为1.369 mm，且整体精度优于未经修正的原始速度场模型(原始E 方向速度标准差为2.664 mm，N 方向速度标准差为1.511 mm).

绘制经最优噪声模型改正后的ITRF14 框架下东北地区垂向运动速度场，如图5 所示.

图5 修正后的垂向速度场

本文获取的东北CMONOC 基准站的垂向速度精度较高，U 方向速度标准差为0.892 mm，且整体精度优于未经修正的原始速度场模型(原始U 方向速度标准差为0.981 mm).

对比修正前后的速度场模型可知：经最优噪声模型修正后的速度场精度明显优于未修正的速度场且二者在水平方向上的运动速率最大差值为0.564 mm/a；垂直方向上的运动速率最大差值为1.285 mm/a.因而，在对时间序列数据进行处理以及解算速度场之前需要考虑CN 的最优模型.

4 顾及CN 及环境负载的东北地区速度场估计

在计算环境负载时采用GFZ (德国地学中心)的数据：非潮汐大气负载时间分辨率为3 h，空间分辨率为0.5°×0.5°；非潮汐海洋负载时间分辨率为3 h，空间分辨率为0.5°×0.5°；水文负载时间分辨率为24 h，空间分辨率为0.5°×0.5°[16].在解算时选取与各站点一致的观测时段，这样可以避免环境负载中非线性趋势对原始时间序列数据的影响，进而避免在测站速度估计时产生偏差.图6 为各测站不同方向上时间序列环境负载的平均值.

图6 各测站环境负载时间序列图像

结合GFZ 的负载数据分析，由图6 可知U 方向上环境负载对东北地区时间序列引起的负荷形变尤为严重.

由表3 可知，在考虑环境负载后，部分测站的噪声会发生变化：在N 方向上，最优噪声模型WN+FN的占比提高至75%；在E 方向，最优噪声模型WN+PL的占比提高至68.75%；U 方向上最优噪声模型WN+FN的占比降低至81.25%.结合环境负载影响对东北地区16 个陆态网络连续观测站的N、E、U 方向速度场进行估计，结果如表4 所示.

表3 环境负载改正后各方向最优噪声模型变化

表3 (续)

表4 东北陆态网络连续站区域速度场估计 mm/a

东北地区陆态网络连续站总体上向E 方向沿顺时针运动这与文献[3]的研究结果一致.在垂直方向上，可以看出北部地区垂直运动速率相对较低且垂直向运动速率差值较大，进而可以推断西南部地区沉降速率大于东部地区[4]，进一步证明了文献[4]研究结果的可靠性.

5 结束语

本文采用东北地区16 个陆态网连续站10 a的坐标时间序列观测数据，利用BIC 信息准则确定了东北地区各坐标分量对应的最优噪声模型，并在考虑环境负载的基础上确定了修正后的东北地区陆态网络速度场.最后得到如下结论：

1)通过10 a 的观测数据得到了较为稳定且可靠的三个方向分量时间序列与最优噪声模型之间的关系，且N 分量和E、U 分量具有明显不同的噪声特性：在N 方向和U 方向上，最优噪声模型均为WN+FN；在E 方向上，最优噪声模型为WN+PL.

2)考虑CN 及环境负载后，东北地区陆态网络基于ITRF14 框架下在N 方向上运动的平均速率为-13.003 mm/a，E 方向上运动的平均速率为27.020 mm/a，U 方向上运动的平均速率为0.528 mm/a.

3)考虑CN 及环境负载的速度场精度明显优于未修正的速度场，二者在水平方向上的运动速率最大差值为0.802 mm/a；垂直方向上的运动速率最大差值为1.377 6 mm/a.因此，在处理东北地区时间序列数据以及解算速度场时考虑CN 及环境负载影响是十分必要的.