采用方差-协方差分量估计GPS时间序列噪声特性
张旭飞
【摘 要】本文采用方差-协方差分量估计分析GPS残差时间序列噪声特性.介绍了该方法如何运用于GPS时间序列分析,详细的推导了函数模型,建立了数据处理流程.对比传统的极大似然估计,该方法可以定量计算各噪声分量的大小,并且具有计算速度快,数学模型严谨等优点.
【期刊名称】《北京测绘》
【年(卷),期】2017(000)004
【总页数】5页(P33-37)
【关键词】方差分量估计;GPS时间序列噪声分析;极大似然估计
【作 者】张旭飞
【作者单位】河北省地矿局石家庄综合地质大队,河北石家庄050085
【正文语种】中 文
【中图分类】P228.4
运用GPS对地壳进行监测,是以一系列固连在稳定基岩上的地面观测墩进行观测,其主要产品是测站的精确位置信息,即时间序列和速率场[1,2];这些产品反映的是地表与地下构造的运动情形,可以提供有关更多的地壳形变信息。连续GPS观测可以提供许多震前、同震及震后变形资料,让我们更加了解地震的震源特性、地壳的应变累积与能量释放过程。
在多数GPS时间序列的研究中,仅将观测资料的误差视为与观测时间无关的白噪声(White Noise,WN),因为白噪声的数值模型及计算较为容易。20世纪90年代末开始对GPS时间序列的噪声特性分析进行大量的研究[3-6],这些研究主要使用两种计算方法:频谱分析法和最大似然估计(Maximum Likelihood Estimation,MLE)法。结果表明各站的GPS残差序列在空间和时间上不完全独立,除了白噪声还有明显的幂律噪声(power law noise)成分;GPS残差序列的幂律噪声成分主要是闪烁噪声(Flicker Noise,FN)或随机游走噪声(Random walk Noise,RN)。根据前人研究指出,若不考虑与时间相关的有噪声(Color Noise)的存在将会低估地壳形变速度的误差。国内的黄立人[7-8],田云峰[9]等利用国内大量的GPS资料,验证了GPS时间序列的最
佳噪声模型为“WN+RN”形式。
本文试图采用一种新的方法来探求GPS时间序列中的噪声特性。对比MLE方法,方差-协方差分量估计法具有如下优势[10-11]:1)是一种无偏估计且有最小方差估值;2)计算速度快;3)任何噪声数据都可以建立随机模型。
GPS坐标时间序列单分量yti可以用如下公式表达[12]:
yti=a+bti+csin2πti+dcos2πti+esin4πti
+fcos4πti+ j=1nggjHti-Teq + j
=1nhhjHti-Tpostti+j=1nkkjexp-ti
-TpostτjHti-Tpost + vi
上式中,ti(i=1…N)表示观测时间,单位为年;系数a、b分别表示测站时间序列起始位置和速度;系数c、d表示测站周年项运动,e、f表示本周年运动;H*表示阶跃函数;系数gj表示地震发生时刻Teq发生的突变值,ng表示有n个地震发生(nh、nk是一样的含义);hj表示地震发生之后,对测
站速度的变化率;kj表示的是在震后Tpost时刻松弛或者滑移位移大小;τj表示测站震后松弛时间常数；vi表示测量误差,其一般假定为与时间无关,也就是期望Evi=0。
2.1 函数模型
对于某一个特定的地震(比如,尼泊尔地震),发生时刻和震后时间是已知的,也就是上式中Teq、Tpost已知,ng=nh= nk=1 。测站震后松弛时间常数τj的最佳值可以利用模型不拟合度最小(Model Misfit)来确定;如此,上式可以利用最小二乘估计求解各个参数的最佳估值。将上式写出矩阵形式,如下:
其中,
y=[y1,y2,…,yn]T
B = 1t1sin2πt11t2sin2πti⋮⋮⋮cos2
πt1sin4πt1cos4πt1cos2πt2sin4πt2cos4πt2⋮⋮⋮
Ht1-TeqHt1-Tpostt1expf()(-(t1-
Tpost)/τ)Ht1-TpostHt2-TeqHt2-
Tpostt2expf()(-(t2-Tpost)/τ)
Ht2-Tpost⋮⋮⋮1tnsin2πtn
cos2πtnsin4πtncos4πtnHtn-TeqHtn-
Tposttnexpf()(-(tn-Tpost)/τ)Htn-Tpost
x= a,b,c,d,e,f,g,h,kT
∈=v1, v2,…,vnT
构建平差误差方程如下:
V=Bx-y
按最小二乘原理,上式的x必须满足VTPV=min的要求,因为认为上面7个待估参数为独立量,故可以按数学上求函数自由极值的方法,得
∂VTPV∂x=2VTP∂V∂x= VTPB=0
转置后可得
BTPV=0
将误差方程带入上式,可得
BTPBx-y= BTPBx-BTPy=0
x=BTPB-1BTPy
2.2 随机模型
一般认为GPS信号中的噪声时间序列可以利用功率谱(power-spectrum)过程来描述[3,5-6]:
Pxf= P0ff0k
其中,f是时间频率(temporal frequency);P0和f0是正则化常数;k是谱指数(spectral index),通常介于-3到1之间,整数k代表一些特殊的噪声类型:k = 0时是标准的白噪声、k =-1时是标准的闪
烁噪声、而k =-2时则是标准的随机漫步噪声。另外,除了标准的白噪声以外,其余的部分统称为有噪声。
因此在构建平差模型时,随机模型表示如下
D= σw2Qw+ σf2Qf+ σrw2Qrw
上式中,σw2、σf2和σrw2分别表示白噪声、闪烁噪声和随机游走噪声的方差;Qw、Qf和Qrw分别表示白噪声、闪烁噪声和随机游走噪声协因数矩阵(cofactor matrix);根据文献[13]各类噪声的协方差阵可表示成转换矩阵Tk与其转置的乘积,转换矩阵通常可以表示为:
Tk= ΔT-k/4φ000…0φ1φ00…0φ2φ1φ0…0
转换矩阵中ΔT为采样间隔(GPS时间序列日解,时间单位为1天)
1)k = 0为白噪声,白噪声与时间无关,对应的转换矩阵
Tk=0=1⋱1 , Qw= Tk=0Tk=0T= In
2) k =-1为闪烁噪声,按照转换的形式,闪烁噪声协因数阵形式复杂,因为绝大多数的空间大地测量时间序列维数小于222,即i-j << 222,所以可以用其简单形式
qijf= 342×2
i=j342×(2-logi-jlog2+212) i≠j
例如,包含1000个点的时间序列的闪烁噪声协因数阵如下(对称阵)
Qf=1.1251.0310.9840.957…
0.5641.1251.0310.984…0.5641.1251.031…
0.564⋱⋱⋮1.1251.0311.125
3)k =-2为随机游走噪声,于是φn = 1,则转换矩阵如下
Tk=-2= ΔT1/2100…0110…0111…
Qrw= Tk=-2Tk=-2T= ΔT111…
1122…2123…3⋮⋮⋮⋱⋮123…n
如上所述,确定了各噪声分量的协因数阵以后,可以利用方差-协方差分量法求σw2、σf2和σrw2最佳估值。
下面介绍具体求解σw2、σf2和σrw2的步骤[14,15],为了叙述方便,表示如下:
θ= σw2σf2σrw2T=θ1θ2θ3T
Q= QwQfQrwT= Q1Q2Q3T
D= σw2Qw+ σf2Qf+
σrw2Qrw= i=13Qiθi
为平差时定权的需要,应取θ的初值,一般取
θ0= 111T
带入随机模型可得D0,由此可得权阵P0=D0-1,并且协因数阵与权阵互为逆矩阵。
D0= Qw+ Qf+ Qrw= i=13Qi正则化协方差
P0=D0-1·D0D0-1=
D0-1Qw+ Qf+ QrwD0-1
P0=D0-1QwD0-1+D0-1QfD0-1+D0-
1QrwD0-1= P1+ P2 + P3 = i=13Pi
参照式(10),x= N-1BTP0y
N=BTP0B
则带入误差方程可得
V=Bx-y= BN-1BTP0y-y=
BN-1BTP0-Ey=R0y
式中,R0=BN-1BTP0-E;其为幂等矩阵,具有如下性质
R0R0=R0 ; R0B=0 ;

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。