顾及有色噪声的GNSS时间序列时域信号提取
时间：2023-04-07 19:10:05 来源：千叶帆 本文已影响人
任安康，徐克科
河南理工大学测绘与国土信息工程学院，河南焦作 454000
GNSS时间序列可以为动态地球坐标框架的建立与维持以及板块运动、冰后回弹等地球动力学现象提供丰富的数据支撑，因此被广泛地应用于大地测量和地球动力学研究(姜卫平等，2018；
姚宜斌等，2019；
李斐等，2019；
马俊等，2021).受到大型地震等因素影响，GNSS时间序列中通常存在震后瞬态，利用标准线性函数模型分析GNSS时间序列无疑会带来函数模型误差，进而导致参数估计有偏.因此，构建合适的GNSS时间序列函数模型是GNSS时间序列分析的先决条件.
对于存在震后瞬态的GNSS时间序列，合适的GNSS时间序列函数模型包括线性趋势、周期信号、偏移和震后瞬态.其他未建模效应通常使用WN+FN来描述.然而，震后瞬态模型中的特征时间尺度与观测量非线性相关，该值无法使用线性回归方法进行估计.为了构建震后瞬态模型，许多学者(王东振等，2017；
苏利娜和张勇，2018；
姚未正等，2021)使用对数衰减函数建模震后瞬态，并采取试错法确定特征时间尺度.试错法通过遍历一定范围内的时间区间，在所有遍历区间中寻求拟合残差最小的时间常数作为最佳特征时间尺度.尽管该方法计算较为简便，但最终计算结果取决于用户的假定区间，因此该方法存在随机性.Sobrero等(2020)探究了最佳双瞬态模型并使用网格搜索法确定最佳特征时间尺度.网格搜索法计算精度很高，但仍受限于假定区间，其本质与试错法类似.Tobita(2016)提出了利用NLS算法估计特征时间尺度的方法，但该算法基于序列误差项中仅存在WN的假设，对于GNSS时间序列，这种假设显然不成立.因此有必要顾及有色噪声影响，开发出准确计算特征时间尺度的方法来构建GNSS时间序列函数模型，对于获取GNSS时间序列未知参数最优线性无偏估计具有重要意义.
MLE可以同时估计GNSS时间序列函数模型和随机模型中各参数及其不确定度，被认为是目前最准确的GNSS噪声分析方法(Williams 2003).GNSS时间序列中的最佳噪声模型可以用WN+FN模型来描述，WN和FN的比率因站而异，MLE能够确定不同噪声之间的比率.基于此，本文提出顾及有色噪声的GNSS时间序列时域信号提取法.首先基于WN+FN模型，使用MLE对震前GNSS时间序列进行参数估计，并根据估计的参数估值来去除震后GNSS时间序列中的初始位置、线性趋势、周期信号和偏移等信号，以此获取残差序列；
然后将残差序列作为求解特征时间尺度的观测量，WN+FN模型作为观测量的随机模型，并采取NLS算法来估计特征时间尺度；
最后，利用估计的特征时间尺度来构建GNSS时间序列函数模型，并采用MLE估计其未知参数，进而实现GNSS时间序列中时域信号的提取.
1.1 GNSS时间序列函数模型构建
对于存在震后瞬态的GNSS时间序列，其函数模型可以描述为(Bevis et al., 2020)：
+d(t)+ε(t),
(1)
式中，y(t)为t时刻下的基准站位移；
y0为参考历元t0时刻下的基准站初始位置；
v为线性速度；
sk和ck是角频率为ωk的谐波对应的傅里叶系数，nK为谐波个数；
bj为发生在tj时刻的偏移量，nJ为偏移个数；
H代表Heaviside阶梯函数，阶跃前为0，阶跃后为1；
d(t)为震后位移；
ε(t)为随机过程.根据Tobita(2016)的研究，顾及震后余滑和黏弹性效应双重影响，震后位移d(t)可以描述为
(2)
式中，Al和Ae分别为对数、指数衰减函数的振幅；
Tl和Te分别为对数、指数衰减函数的特征时间尺度；
Δt为地震发生以来的相对时间，要求满足Δt≥0.特征时间尺度与观测值非线性相关，其值无法使用线性回归估计器进行估计.为了估计特征时间尺度，对式(2)进行泰勒级数展开，其一阶展开式为
式中，dTl和dTe为待估参数.根据式(3)中参数和测量值之间的关系，构建如下观测方程：
(4)
(5)
式中，N=ATC-1A；
L=ATC-1V；
C为观测值协方差阵.Levenberg-Marquardt算法用于解决非线性参数估计，使用以下规则定义的新矩阵N′=[n′ij]替换矩阵N(Mao et al., 1999)：
(6)
式中，λ为阻尼因子.当阻尼因子λ非常大，矩阵N′被迫接近对角，此时由于正规矩阵的对角元素被放大，因此逼近解的步长被缩放.相反，当λ接近零时，N′将接近矩阵N.
利用上述方程，以下给出估计特征时间尺度的伪代码
①选择噪声模型(本文取WN+FN)，使用MLE估计地震前信号.地震前信号包括初始位置、线性趋势、周期信号和偏移.
②去除震后GNSS时间序列中的地震前信号，获取震后残差序列.
③将震后残差序列作为求解特征时间尺度的观测值，并利用震前噪声参数构建观测值的随机模型.
④选择一个适当的值λ，比如λ=0.01.
⑤给定特征时间尺度初值向量ξ，使用线性最小二乘法(随机模型为WN+FN)计算Al、Ae.
1.2 GNSS时间序列随机模型构建
信号建模的最大问题之一是测量噪声性质(Sobrero et al., 2020).众所周知，GNSS时间序列中的随机过程可以表述为WN+FN.MLE可以同时估计时间序列中的函数模型和随机模型中各项参数，因此用来构建GNSS时间序列随机模型.
假设GNSS时间序列中的随机过程ε(t)由振幅分别为σw和σf的白噪声α及闪烁噪声β组成(Zhang et al.1997; Amiri-Simkooei et al., 2019)：
ε(t)=σw·α(t)+σf·β(t)，
(7)
其观测值协方差阵为
(8)
式中：I为单位阵；
Jf为闪烁噪声的协方差阵，Jf的构造方式可以参考相关文献(Zhang et al.1997).
1.3 GNSS时间序列模型参数估计
GNSS时间序列函数模型未知参数可以使用加权最小二乘法进行估计，其最小二乘解为(Bos et al., 2013)
(9)
(10)
(11)
(12)
为了从线性函数中估计噪声分量和参数，对于给定观测量x，必须使这些值的可能性l最大化.假设为高斯分布，可能性为(Williams，2003)
(13)
式中：det为矩阵行列式；
N为历元数.求解最大似然问题的算法可以选择下山单纯形法.
为了研究所提出方法的计算性能，本节使用了合成GNSS时间序列进行分析.本节强调了两个问题.(1)研究随机模型对特征时间尺度的影响.(2)研究随机模型对GNSS时间序列函数模型未知参数的影响.研究均表明，GNSS时间序列的最佳噪声模型可以用WN+FN(Williams et al., 2004; 李昭等, 2012).因此本文采用WN+FN模型来模拟噪声项.表1中，噪声参数来源于文献(Williams et al., 2004)，该值是对区域GNSS时间序列噪声分析得到的经验值.信号参数来源于日本全球定位系统永久性跟踪站网(GPS Earth Observation Network System, GEONET)实测数据，该值是利用MLE对GEONET GNSS时间序列进行参数估计获取.其中，构建GNSS时间序列函数模型的特征时间尺度参数来源于文献(Sobrero et al., 2020)，该值代表了东日本大地震震后位移特征时间尺度的最小值、中值和最大值.根据表1中参数设置，本文创建了三个GNSS永久站水平分量合成时间序列，这些合成时间序列包括线性趋势、年和半年周期信号、偏移、震后位移和噪声.对于每个时间序列，其噪声项基于表1中规定的噪声振幅利用Hector软件生成了50个随机误差向量.然后，将模拟的随机误差向量分别添加到上述确定性模型中构成模拟数据集.图1显示了合成时间序列的一个典型示例以及模拟同震位移发生时刻.
图1 3个永久GNSS站合成时间序列典型示例
表1 模拟时间序列参数
第一个目标是研究随机模型对特征时间尺度估计的影响.与真实数据相比，模拟数据的优势在于，特征时间尺度的大小是完全已知的.我们将对前一节中描述的NLS算法估计特征时间尺度的计算性能进行统计分析.表2中设置了2种解算方案.在案例1中，使用NLS算法估计特征时间尺度.在这种情况下，时间序列的协方差阵仅由WN组成.在案例2中同样使用NLS算法计算特征时间尺度，但其协方差阵同时考虑了WN和FN.最后，根据表3中参数初值设置利用NLS算法计算了每个合成时间序列的特征时间尺度.第二个目标是研究随机模型对GNSS时间序列函数模型未知参数估计的影响.类似地，根据表2中的解算方案，利用MLE法计算了每个合成时间序列的未知参数.以下重点介绍模拟实验结果：
表2 GNSS时间序列分析解决方案
表3 特征时间尺度初值
图2 NLS算法估计特征时间尺度收敛性百分比直方图
图3 使用NLS算法执行50次独立运行计算获取的特征时间尺度散点图及其设计值
最后，为考虑解决方案的整体性能，根据表2中的解算方案对所有合成时间序列的未知参数进行了估计.图4显示了估计值的误差分布.经统计，案例2中99.89%的线性速率误差为-0.5～0.5 mm，而案例1中89.53%的线性速率误差为-1～1 mm.案例2中线性速率的STD为0.14 mm，而在案例1中其值为0.48 mm.案例1中线性速率的STD比案例2大2.43倍.对于周期项振幅而言，两种案例的结果无显著差异.然而，对于同震位移和震后瞬态振幅，两种案例的结果具有显著差异.在案例2中99.11%的同震位移相对误差为-0.5%～0.5%，而案例1中86.70%的同震位移相对误差为-20%～1%.在案例2中同震位移相对误差的STD为0.18%，而在案例1中其值高达7.83%，这表明忽略有色噪声可能会带来严重的模型偏差.对于对数项振幅，案例2中96.89%的对数项振幅误差为-20～20 mm，而案例1中86.63%的对数项振幅误差为-40～40 mm.对于指数项振幅，案例2中84.56%的指数项振幅误差为-20～20 mm，而案例1中仅74%的指数项振幅误差为-40～40 mm.案例1中对数项振幅和指数项振幅的STD分别比案例2大2.55倍和1.29倍.综上分析表明有色噪声对GNSS时间序列信号提取至关重要.
图4 未知参数估值误差分布(设计值为基准)
3.1 震后瞬态信号检验与数据处理
为验证顾及有色噪声的GNSS时间序列时域信号提取算法的实际应用效果，利用美国内华达大学大地测量测量实验(http:∥geodesy.unr.edu/NGLStationPages/)提供的日本GEONET中83个水平分量坐标时间序列(截止观测时间为2020.15年).这些台站所在区域于2011年3月11日发生MS9.0级大地震，致使震后时间序列存在明显的同震位移和震后形变，其位置分布见图5.
图5 实测数据站点分布(截止观测时间为2020.15年)
为进一步定量分析测站水平分量震后余滑及黏弹性效应的显著性，本文采用DIA估计器对其检验.DIA的检验方法为(Amiri-Simkooei et al., 2019)
(14)
式中，零假设H0为不存在显著的震后余滑及黏弹性效应,备选假设H1为存在显著的震后余滑及黏弹性效应.E为期望算子；
A为m×n设计矩阵；
y为m维时间序列观测向量；
x为n维未知参数向量；
D(y)为描述观测量的随机模型；
在备选假设H1中，xk=[AlAe]T；
Ak为对数和指数衰减函数对应的m×2维矩阵，其定义为
(15)
为了确定式(14)中的两个备选假设，构造如下检验统计量(Amiri Simkooei 2013)：
(16)
图6 实测数据DIA检验结果
分析图6可以发现各测站的T2值均超过卡方分布的临界值，表明在99.9%的置信区间内，各测站均存在显著的震后余滑及黏弹性效应.最后，根据表2中的解算方案对实测数据进行信号提取以验证该方法的实际应用效果.图7总结了两种方案的残差RMS统计数据.结果表明案例2的拟合性能显著优于案例1.当使用案例2时，其RMS对应的四分位数明显小于案例1，这表明忽略有色噪声会降低解决方案拟合性能.为进一步分析两种方案的拟合性能，图7显示了RMS的经验累积分布函数.经统计，案例2中约70%的RMS比案例1大1 mm，因此精确的GNSS时序拟合需要顾及有色噪声.对于提取的各类信号，下文将重点进行分析.
图7 残差序列RMS统计结果
3.2 GNSS测站速度和速度不确定度分析
根据表2中解算方案，本文估计了83个GNSS站速度.图8给出了两种方案计算的日本区域水平速度场.从整体上分析，两种速度解呈现向东南方向运动趋势，但部分站点的结果差异较大.其中，J193测站(见图8黑色三角标记)的运动趋势均发生了显著的反向变化.在案例2中J193站估计的水平运动速度为13.90 mm·a-1，方向为EN2°，而案例1中其水平运动速度为23.01 mm·a-1，方向为SW3°.为了进一步分析该测站解算结果，我们计算了该测站震前站速度(见图8绿色矢量，已等比例增大2倍).与案例1相比，案例2估计的该测站结果明显更接近震前速度，这说明考虑有色噪声更有利从震后时间序列中获取稳态速度.除J193测站外，其他测站同样存在不同程度差异.为了对比两种方案估计的速度结果，本文统计了案例1与案例2获得的GNSS站速度差值，其结果见图9.对于东分量时间序列，92%的速度差值为-10～10 mm·a-1，最大的速率差值为36.86 mm·a-1；
对于北分量时间序列，96%的速度差值为-4～4 mm·a-1，最大的速率差值为8.07 mm·a-1.综合上述分析表明：(1)对于日本区域的GNSS速度场，必须顾及噪声模型带来的差异.(2)有色噪声对日本区域水平分量GNSS站速度的影响并不均衡，其中有色噪声对东分量GNSS站速度的影响明显大于北分量.
图8 日本区域GNSS速度场；
蓝色矢量为表2中案例1；
红色矢量为表2中案例2
图9 案例1与案例2解算的GNSS速度差值统计结果
进一步从速度误差来分析，案例1计算的速度不确定度均小于案例2的结果.经统计，顾及有色噪声所计算的速度不确定度要比仅考虑WN计算的速度不确定度大1～8倍，该值大于中国香港地区(袁林果等，2008)及国家连续运行参考站系统(蒋志浩等，2010)结果，主要原因可能是日本区域震后时间序列中存在大量非线性信号导致有色噪声振幅偏大所致.
3.3 季节信号分析
现有研究表明GNSS时间序列中除趋势项外，还存在一定的季节信号.提取GNSS时间序列中的季节信号具有双重意义：其一，对GNSS时间序列中季节信号进行建模可以改善GNSS时间序列拟合性能，提高趋势项不确定度.其二，联合提取的季节信号和地球物理模型能够用于探究地表季节性形变的已知物理成因.本文利用表2中的解算方案，提取日本区域83个基准站季节信号，其统计结果见表4—5.
表4 基于WN模型估计的周期信号振幅统计(mm)
表5 基于WN+FN模型估计的周期信号振幅统计(mm)
从统计结果来分析，在两种方案中东分量的周年项振幅普遍大于半周年项振幅，而北分量的周年项和半周年项振幅差异不大.在案例1中东分量的周年项振幅中值分别为0.8 mm，而半周年项振幅中值为0.35 mm；
在案例2中东分量的周年项振幅中值为0.78 mm，而半周年项振幅中值分别为0.37 mm.因此，平均而言东分量的周年项振幅比半周年项振幅大.
进一步从噪声对季节信号的影响来分析，本文统计了案例1与案例2获得的周期项振幅差值，其结果见图10.从统计结果来看，98%的周年项振幅差值为-0.40～0.40 mm，最大的差值为0.80 mm(见图10a,b,c)；
相对周年项振幅，半周年项振幅差值普遍较小，约93%的振幅差值-0.2～0.2 mm，最大的差值为0.35 mm，这表明有色噪声对周年项振幅的影响大于半周年项振幅.此外，观察图10d,e,f可以发现北分量周期项振幅差值分布呈现明显的不均衡迹象.其中，在北分量中71.08%的周年项振幅差值为负值(东分量为48.19%)，89.16%的半周年项振幅差值为负值(东分量为54.22%)，这表明忽略时间序列中的有色噪声会引起日本区域北分量周期项振幅低估.因此，在对日本区域进行地表季节性形变分析时，需顾及有色噪声影响以避免北分量上季节性形变信号的掩盖.
图10 周期信号振幅统计结果
3.4 同震和震后位移分析
利用估计的参数，本文分析了各测站的同震位移.图11显示了两种方案的估计结果.从整体上分析，两种方案的同震位移均指向震中方向，且同震位移的量级随着距震中距离的增加而急剧减小，其中高值区约为2.5～4 m，中值区约为1.5～2.5 m，低值区为0.5～1.5 m，平均方向为ES30°.为了进一步对比两种方案的差异，本文还计算了两种方案的同震位移差值(见图11灰色矢量).结果表明，同震位移的差值具有明显的空间模式.在靠近震中区域，J179和J913测站的同震位移差值明显大于相邻测站(见图11黑色三角标记).其中，J179和J913站的同震位移差值分别为32.18 mm和34.69 mm，而相邻测站的同震位移差值普遍较小，说明两个站可能存在强烈局部效应，其原因还需进一步分析.总体而言，靠近震中区域的同震位移差值普遍较小.相反，在同震形变最小的北部区域，但两种方案估计的同震位移差值普遍较大.同震位移差值约为-30～-10 mm，最大的同震位移差值为-58.54 mm，这表明忽略有色噪声会使微量级同震信号被低估.因此，对于高精度地壳形变分析有必要考虑有色噪声的影响.
图11 同震位移场
根据上述分析结果，日本区域GNSS时间序列中除同震位移外，还存在显著的瞬态信号.利用估计的参数，本文计算了该区域震后9年的水平地表位移，其结果见图12.从整体上分析，震后位移与同震位移趋势一致即整体指向震中方向，但其变化梯度相对较小，整体量级约为0.5～1.7 m，平均方向为ES26°.为了进一步对比两种方案的差异，本文还计算了两种方案的同震位移差值(见图12灰矢量).从震后位移差值的空间分布来分析，最大的震后位移差值分布于震中区域，分别为测站G145、G159和J193(见图12黑色三角标记)，其量级分别为365.72 mm、287.32 mm和527.93 mm.除上述3个测站外，其余震后位移差值的大小和方向呈现随机性，其量级约为-100～100 mm.因此，在对GNSS时间序列进行瞬态信号提取时，其噪声特性不容忽视.
图12 震后位移场(震后9年)
GNSS时间序列中可能存在震后瞬态，忽略震后瞬态无疑会引起参数估计有偏.为此，本文结合GNSS时间序列噪声特性，提出了顾及有色噪声的GNSS时间序列时域信号提取法.该方法首先基于WN+FN模型，使用MLE对震前GNSS时间序列进行参数估计，并根据参数估值来去除震后GNSS时间序列中震前信号，以此获取残差序列；
然后将残差序列作为求解特征时间尺度的观测量，WN+FN模型作为观测量的随机模型，并采取NLS算法来估计特征时间尺度；
最后利用估计的特征时间尺度构建GNSS时间序列函数模型，并采用MLE估计其未知参数，进而实现时域信号的提取.经模拟数据和实测数据将该方法与传统方法进行了比较分析，其主要结论如下：
(1)对日本区域扣除震前信号后的水平残差时间序列进行DIA检验表明，在99.9%的置信区间内日本区域83个GNSS站均存在显著的震后余滑及黏弹性效应.
(2)有色噪声对日本区域GNSS站速度具有显著影响.相对WN+FN模型，在东分量仅考虑WN模型引起的速度差值约为-10～10 mm·a-1，最大差值高达36.86 mm·a-1；
在北分量仅考虑WN模型引起的速度差值约为-4～4 mm·a-1，最大差值为8.07 mm·a-1.
(3)有色噪声对日本区域周年项振幅的影响大于半周年项振幅.相对WN+FN模型，仅考虑WN模型引起的周年项振幅差值约为-0.40～0.40 mm，最大差值为0.80 mm；
相对周年项振幅，假设WN模型引起的半周年项振幅差值普遍较小，其值约为-0.2～0.2 mm，最大差值为0.35 mm.
(4)忽略有色噪声会引起日本区域微量级同震信号的低估，其量级约为10～30 mm，最大量级为58.54 mm；
相反，有色噪声对震后位移的影响则呈现为随机性，其量级约为-100～100 mm.
