在研究地壳形变问题时,常利用流动GNSS复测资料计算指定区域的位移场、速度场和应变场,如利用GNSS流动资料研究中国大陆应变场及块体运动,取得了丰硕的成果(江在森等,2003; 王敏等,2003; 杨国华等,2002; 张培震等,2003,刘志广等,2013)。但是流动GNSS资料复测周期长,且其观测资料的质量、点位的稳定性等各方面都与连续观测资料存在较大的差距。随着《中国地壳运动观测网络》二期工程的实施,GNSS连续观测站由33个增至253个,由地方或其他部门建设的GNSS连续站也已超过千个,这使得利用GNSS连续观测资料进行形变场的分析成为可能。若发挥GNSS连续观测资料时间上的优势,计算并分析月、季度、半年和年各个时段尺度下的形变场的动态演变过程,可使GNSS连续观测资料更有效地服务于震情会商和形变分析,因此利用GNSS连续观测资料及时计算获取各个时段尺度下稳定可靠的速度场是非常必要的。GNSS连续观测时间序列中除了线性成分外,还包含多种周期运动(江在森等,2006; 王敏等,2005a; 朱文耀等,2003)。白噪声、闪烁噪声和随机漫步噪声的噪声性质是GNSS连续观测时间序列的基本特征(黄立人,2006; 黄立人,符养,2007; 蒋志浩等,2010,夏峰等,2014)。如何用GNSS时间序列计算不同时段下的高精度速度场,获取稳定速度场与GNSS时间序列时段长度的关系是怎样的?周期性非构造运动、随机漫步噪声和闪烁噪声等有色噪声对不同时段长度的GNSS时间序列速度场的影响有多大?本文以上述问题为研究背景,采用小波分析方法和加权最小二乘法评定顾及周期性非构造运动和不同噪声影响的中国GNSS连续观测网的月时段、半年时段、年时段的速度估值和精度,计算不同时间尺度下稳定的GNSS连续观测速度场。

本文使用的GNSS数据是中国地壳运动观测网络253个基准站运行以来记录的连续观测资料,其中一期网络基准站33个,观测时间为1999年初至2013年4月的数据,二期网络基准站220个,观测时间为2010年至2013年4月的数据,参与解算的还有90个IGS基准站在相应时间段内的观测数据。GNSS连续站观测值的数据处理采用GAMIT/QOCA软件完成(Herring et al.,2010)。数据处理的基本流程(王敏等,2005b; 王敏,2007)首先利用GAMIT获得网络工程253个连续站及90个IGS测站的单日松弛解,同步观测的测站较多时,采用分区处理。完成GAMIT计算之后,利用QOCA软件将计算所得的各区单日松弛解进行综合平差,在此基础上通过IGS核心站求解相对于全球参考框架ITRF2005的相似变换7参数,从而获得ITRF2005下的单日解,即GLOBK的NEU坐标值,不同于通常的站心坐标,其具体含义可参阅黄立人等(2006)。对计算得到的各测站坐标时间序列,在剔除粗差的基础上,我们进行了2001年昆仑山口西8.1级地震、2004年苏门答腊9级地震、2008年汶川8.0级地震和2013年庐山7.0级地震同震位移改正的处理。

通过剔除GNSS站坐标时间序列的粗差,仅可以削弱观测值中的偶然成分,还需要对信息进一步分离。GNSS连续观测时间序列中包含着多种有规律的、周期性成分,研究表明,在GNSS连续观测时间序列中比较显著且共性较好的周期成分是年周期和半年周期(江在森等,2006)。可采用谐波拟合法(王敏等,2005b; 韩鹏等,2009; 谢凡等,2011)剔除GNSS连续观测时间序列中的年周期和半年周期成分。谐波拟合方法是在傅里叶级数及傅里叶变换的基础上,把满足一定条件的周期函数分解为组成它的各次谐波成分。对GNSS连续观测时间序列去线性趋势后得到的残差时间序列的拟合函数f(t)通过傅里叶级数表示成如下形式:

f(t)=∑^N_m=1(A_mcos2^mπt+B_msin2^mπt).(1)

式中,f(t)可以看成是这些正弦函数和余弦函数的“加权和”,其中m为谐波的次数,系数A_m、B_m是权重因子,表示各次谐波对总序列的贡献,可以通过最小二乘法计算得到。对GNSS连续观测时间序列而言,m取2,即用2次谐波拟合就可以较好的改正GNSS时间序列中的年周期和半年周期成分。

对GNSS连续观测时间序列中的年周期和半年周期成分用谐波拟合的方法进行改正后,笔者采用加权最小二乘法计算并评定未剔除和剔除周期性运动后的测站年速度及精度。最小二乘准则为(武汉大学测绘学院测量平差学科组,2009)

V^TPV=(v₁,v₂,…,v_n)(p₁ 0 … 0

0 p₂ … 0

   

0 0 … p_n)(v₁

v₂



v_n)(2)

=p₁v²₁+p₂v²₂+…+p_nv²_n=[pvv]=min.

式中,p_i为观测值的权,p₁=1/(σ²_i),σ_i为站坐标的精度,v_i为观测值y_i的改正数,v_i=a+bt_i-y_i,(i=1,2,…,n)。由式(2)可求出参数a、b,b即为测站速度。表1给出了未剔除和剔除年、半年周期后计算得到的测站年速度精度对比,由表中可以发现,通过剔除GNSS连续观测站的年周期和半年周期运动,可得到更稳定且精度更高的测站年速度。在剔除年周期和半年周期后,计算得到的测站东向和北向年速度的平均精度比未剔除年、半年周期的结果分别提高约26%和20%。

剔除年、半年周期后计算得到半年速度的精度和年速度一样也有较大幅的提高。表2给出了2003～2012年(共10年)哈尔滨站和泰安站以每半年为一个时段计算所得的东向速度对比。由表中可见,2003～2012年哈尔滨站每年上半年未剔除周期性运动的速度明显小于下半年速度,而剔除周期性运动后各年份半年速度间的差异明显减小。哈尔滨站和泰安站剔除年周期和半年周期后得到的半年时段速度的离散度相对未剔除的结果显著减小,说明剔除年周期和半年周期后的半年时段速度更稳定,可靠性更高。

剔除GNSS连续观测时间序列中的年、半年周期,并计算相对于区域无旋转参考基准(杨国华等,2005; 胡新康等,2007)的速度,笔者计算得到了2010～2013年中国大陆地壳运动观测网络工程GNSS连续观测站的背景速度场(图1中蓝色箭头),其中一期网络33个、二期网络220个GNSS连续站计算所采用时间序列的长度分别为1999～2013年和2010～2013年。为了反映GNSS连续观测站速度场的年度差异运动,笔者利用2011年4月至2012年4月的GNSS连续观测资料,计算了2011～2012年中国大陆GNSS连续观测站速度场(图1中黄色箭头)。由图1可以发现,计算得到的GNSS连续观测站背景速度场与我国地质

构造格局较为相关,基本反映了中国大陆的形变背景。2011年日本9.0级地震后,2011～2012年速度场相对背景速度场的差异变化在一定程度上说明了GNSS连续观测时间序列年速度场能够显示中国大陆地壳水平差异运动的空间分布信息。

由上述分析,剔除GNSS时间序列中的年周期

表1 未剔除和剔除周期后计算得到的GNSS测站年速度精度对比
Tab.1 Precision comparison of the calculated annual velocity of GNSS stations after rejecting the cycle motion or not

和半年周期成分,可以得到稳定且高精度的年时段和半年时段速度场。但是,用上述方法计算了2011年4月初至2013年3月底(共24个月)GNSS连续观测站各月时段的速度,却发现各个月时段速度离散度大,很不稳定,不同年份相对应月份的速度差别也很大,因此要得到稳定的月时段速度,还需对GNSS连续观测时间序列进行更深入的信息分离。

表2 哈尔滨和泰安站半年时段东向速度离散度比较
Tab.2 Comparison of the dispersion of velocity in east recorded by Harbin and Taian Stations in the period of half a year

图1 2011～2012年中国大陆GNSS连续观测站速度场相对背景场的变化
Fig.1 Variation of velocity field compared with background field of GNSS continuous observation stations in China Mainland from 2011 to 2012

为进一步分离GNSS站坐标时间序列所包含的有色噪声,我们采用小波变换的方法。小波分析可以提取时间序列的时频特征,具有分析信号的局部特征和多尺度分辨率的特点,借助小波方法可以将不同特性的噪声进行分离与估计(飞思科技产品研发中心,2005)。

设f(t)是观测点变形序列,采用Mallat算法,将信号分解成不同频率成分:

A_j-1f(x)=A_jf(x)+D_jf(x).(3)

其中,信号在空间V_j上的投影为

A_jf(t)=∑_k∈ZC_j,kφ_j,k(x).(4)

在空间W_j上的投影为

D_j(x)=∑_k∈ZD_j,kψ_j,k(x).(5)

其中,A_j-1 f(x)和A_j f(x)分别是信号f(x)的频率不超过2^-j+1和2^-j的成分,而D_j f(x)是频率介于2^-j与2^-j+1之间的成分。

上述小波分解式实际可写成如下的矩阵形式:

C_j+1=HC_j

D_j+1=GC_j,(j=1,2,…,J).(6)

式中:H是尺度函数对应的低通滤波器,H=(h_k-2n); G是小波函数对应的带通滤波器,G=(g_k-2n); C_j是2^j在分辨率下的离散逼近; D_j是2^j在分辨率下的离散细节; {h_k}_k∈Z和{g_k}_k∈Z是一对离散正交镜像滤波器(低通和高通)。式(6)即为Mallat塔式分解算法。

小波变换要求数据是等间隔采样的,而GNSS时间序列的日观测值中有一定的数据缺失,我们采用最小二乘配置方法对GNSS时间序列进行插值。由于最小二乘配置方法不像谱分析或是小波变换要求数据等间隔,其数据输入可以较为松散,而其输出可以根据需求决定是否等间隔(武艳强等,2007)。这就决定了该方法可以作为谱分析、小波变换等方法的预处理过程。即先由最小二乘配置输出滤波、插值后的等间隔值,然后再由其它方法进行后继处理。

武汉大学测绘学院测量平差学科组(2009)详细论述了最小二乘配置,其误差方程为

{V=B_ZZ^{^}+GY^{^} -L,

V_Z=Z^{^}-L_Z.(7)

式中,B_Z为信号的系数阵(通常观测点部分为单位阵,推估点部分为零阵),Z^{^}为随机信号估值(包括观测点和推估点),G为经典平差问题的系数阵,Y^{^}为经典平差问题需要求解的待定参数,L为观测值,L_Z为信号的虚拟观测值(一般为信号的数学期望)。

根据间接平差理论和矩阵反演理论对式(7)进行求解得到

{Y^{^}={G^T(D_X+D_Δ)^-1}^-1G^T(D_X+D_Δ)^-1L,

Z^{^}=D_Z(D_X+D_Δ)^-1(L-GY^{^}). (8)

式中,D_X为已测点信号的协方差阵,D_Δ为观测值的协方差阵。

应用最小二乘配置方法的前提是获取信号的协方差阵,推荐使用如式(9)的高斯经验协方差函数:

f(t)=Ae^-k2t2.(9)

式中,t为时间间隔,A和k为待定参数。采用赫尔默特方差估计方法获得参数A,而对于参数k,可以根据关注的信息频段人为设定。如果k值较大,得到的信息就接近于原始序列; 如果k值较小,得到的信息就是相对低频的数据。

对最小二乘配置插值后的数据采用具有对称性和正则性的sym8小波进行小波分解(张风霜,杨国华,2011)。以2011～2013年哈尔滨站东向站坐标时间序列小波分解的结果为例(图2),随着阶数的增加,信号的细节部分从最高频向低频分解,当阶数为0时,信号为采样频率; 当阶数为1时,频率二等分,对于GNSS连续观测时间序列,1阶信号的周期为1～2 d,同样,2阶信号的周期为2～4 d,以此类推,当阶数为8时,8阶信号的周期为128～256 d。研究表明,GNSS时间序列小波分解的1～3阶的高频信号主要是白噪声,4、5阶的高频信号主要是闪烁噪声,6、7阶的高频信号主要是随机漫步噪声(杨国华等,2007; 杨博等,2010)。

对哈尔滨站和泰安站已剔除年周期和半年周期的东向站坐标时间序列,分别剔除白噪声、闪烁噪声和随机漫步噪声后再计算2011年4月初至2013年3月底(共24个月)各月时段速度,结果见表3、图4和图3。由图3可见,随着白噪声、闪烁噪声和随机漫步噪声的逐步剔除,24个月各月的速度由离散逐步趋于稳定。由表3、4可知,只剔除白噪声,计算得到的哈尔滨站和泰安站24个月的各月速度离散度较大,分别为0.795和0.464,说明各月速度并不稳定,进一步剔除闪烁噪声后,哈尔滨站和泰安站的月速度离散度分别减小至0.307和0.247,再进一步剔除随机漫步噪声后,哈尔滨站和泰安站离散度分别进一步减小至0.165和0.180,说明剔除白噪声、闪烁噪声、随机漫步噪声后得到的月时段速度较稳定。

图2 HRBN站东向站坐标时间序列小波分析结果
(a)原始信号;(b)低频A8;(c)低频A7;(d)低频A6;(e)高频D8;(f)高频D7;
(g)高频D6;(h)高频D5;(i)高频D4;(j)高频D3;(k)高频D2;(l)高频D1
Fig.2 Wavelet analysis result of time series of the eastward coordinate recorded by Harbin Station(a)original signal;(b)low-frequency A8;(c)low-frequency A7;(d)low-frequency A6;(e)high-frequency D8; (f)high-frequency D7;(g)high-frequency D6;(h)high-frequency D5;(i)high-frequency D4; (j)high-frequency D3;(k)high-frequency D2;(l)high-frequency D1

表3 哈尔滨站东向剔除不同频段噪声后月时段速度对比
Tab.3 Comparison of monthly velocity in east after rejecting different noises recorded by Harbin Station

表4 泰安站东向剔除不同频段噪声后月时段速度对比
Tab.4 Comparison of monthly period velocity of rejecting different noises of the east recorded by Taian Station

图3 剔除不同噪声的哈尔滨东向月时段
速度统计直方图对比
(a)剔除白噪声;(b)剔除白噪声和闪烁噪声;
(c)剔除白噪声、闪烁噪声和随机漫步噪声
Fig 3 Statistical histogram comparison of monthly velocity after rejecting the different noises recorded by Harbin Station (a)rejecting white noise;(b)rejecting white and

通过剔除GNSS连续观测时间序列的年周期和半年周期成分,可得到稳定、高精度的年时段和半年时段速度场,但由于GNSS连续观测时间序列的噪声性质,短时段稳定速度场的获取需要对GNSS连续观测时间序列进行滤波剔除有色噪声。本文采用小波分析方法分别剔除了GNSS连续观测时间序列中的白噪声、闪烁噪声和随机漫步噪声,并比较了逐步剔除不同频段噪声后计算得到的各月时段速度的精度和离散度,结果表明,剔除白噪声、闪烁噪声、随机漫步噪声后得到的月时段速度更稳定可靠。利用本文的方法,实时计算并分析月、半年和年各个时段尺度下的形变场的动态演变过程,有利于发挥GNSS连续观测资料时间上的优势,从而使GNSS连续观测资料更有效地服务于震情会商和形变分析。