基于Delauany三角网与反距离内插实现GPS区域网中噪声剔除

纪海源 张伟佳 梁 磊 何远梅 郝 明

(1. 陕西工业职业技术学院 土木工程学院, 陕西 咸阳 712000; 2. 陕西新旅程测绘科技有限公司,陕西 西安 712000; 3. 中国地震局第二监测中心, 陕西 西安 710054)

0 引言

中国陆态网络工程的区域站已经有两千多个了,这些区域站的观测方式不是连续的,而是采取流动观测的方式,其优点是可以节约人力、财力、物力,观测效率高,具有很强的灵活性;缺点是得到的坐标时间序列不连续,而在全球定位系统(Global Positioning System，GPS)观测结果中存在非构造噪声是有目共睹的,时间序列的不连续无法使用QOCA(Quasi-Observation Combination Analysis)地壳形变分析软件进行非构造噪声的剔除。因此区域站速度所反映的地壳形变中包含有大量的误差,剔除区域站时间序列中的非构造噪声成为国内外学者研究的热点[1]。田云峰等认为GPS坐标中非构造噪声存在区域相关性[2-3],于是设想到区域站分布和连续站分布的关系,是否能寻找某种方法将连续站和区域站关联起来。 Delaunay三角形能把三角形顶点和内部的点联系起来,用连续站构建Delaunay三角形成本文新的探索方法;其次采用反距离权内插将区域站和Delaunay三角形的顶点联系起来,用连续站的周年项和半周年项参数(QOCA软件结果)去拟合区域站的周年项和半周年项参数,最终剔除区域站中的部分非构造噪声,提高GPS流动观测区域站三维运动速度场的精度,为深入认识地块动力学机理提供重要基础资料。

1 连续站时间序列参数的提取

本文研究区域为中国地震监测陆态网在青藏高原东北缘地区(101°～107°E,32°～40°N)分布的GPS连续站和流动观测区域站,分布如图1所示;采用GAMIT/GLOBK软件处理13个连续站和153个流动观测区域站,得到连续站高精度位移序列,再使用QOCA地壳分析软件[2]的analyze_tseri和pca模块剔除GPS连续站中的非构造噪声,在QOCA软件的结果文件中,将三个方向分析拟合出来的周年项、半周年项参数等提取出来(表1～3)。s1、s2和s3、s4分别为年、半年周期项系数。

图1 汾渭地区GPS连续站和区域站分布图注：五角形为GPS连续站；圆形为区域站。

表1 北向参数提取结果

表2 东向参数提取结果

表3 垂向参数提取结果

2 Delaunay三角网的构建

2.1 VB实现Delaunay三角网

狄洛尼(Delaunay)三角网是Voronoi图的对偶图,由对应的Voronoi多边形共边的点连接而成。狄洛尼三角形由三个相邻点连接而成,这三个相邻点对应的Voronoi多边形有一个公共顶点,此顶点同时也是狄洛尼三角形外接圆的圆心[4-5]。本文采用VB(Visual Basic)编程实现了13个连续站狄洛尼(Delaunay)三角网的生成(图2)。

图2 Delaunay三角网

2.2 区域站在Delaunay三角网中的识别

在生成三角网之后,需要判断每个区域站隶属于哪个三角形。判断点在三角形内的算法有好多种,在这里选择面积法:如果区域站与三角形三个顶点的连线构成三个小三角形的面积之和与大三角形的面积相等,则确定点落在这个三角形里面(图3)。

图3 部分区域站在Delaunay三角网中的分布

3 根据连续站拟合区域站参数

3.1 反距离权算法简介

“反距离加权算法”又称为倒数距离加权插值或“Shepard方法”[6]。设有n个点,平面坐标为(xi,yi),垂直高度为zi,i=1,2,…,n,倒数距离加权插值的插值函数为

f(x,y)=

3.2 VB实现区域站参数拟合

将提取的连续站周年项、半周年项参数准备好,VB编写程序实现反距离权内插的程序,用三个顶点的参数拟合出区域站的周年项、半周年项参数。表4～6是提取的部分三角形三个方向拟合的区域站参数结果。

表4 北向参数拟合结果

表5 东向参数拟合结果

表6 垂向参数提取结果

4 区域站非构造噪声的剔除

在上一节已经将区域站的周年项、半周年项参数拟合出来了,可以应用数学模型计算其中包含的非构造噪声,然后从位移序列中剔除;对位移序列分量每日解观测序列建立参数模型:

y(ti)=a+bti+s1×sin(2πti)+s2×

cos(2πti)+s3×sin(4πti)+s4×cos(4πti)(2)

式中,ti(i=1…N)为以年为单位的时间,在这里我们以每日的单日解构成时间序列;a为初始位置;b为速率;s1、s2和s3、s4分别为年、半年周期项系数[7]。

将表5中拟合的参数按照测站代入位移序列的数学模型,把计算出来的非构造噪声从位移序列y(ti)中剔除,这样可以得到相对精确的地壳运动形变,下面举例G120站的剔除结果。

从表7的噪声列可以看到北南方向误差最大可以剔除掉1.2 mm,最小为0.5 mm,噪声的均方根为0.89 mm;东西方向误差最大可以剔除掉1.0 mm,最小为0.2 mm,噪声的均方根为0.66 mm;垂直向最大可以剔除掉3.5 mm,最小1.2 mm,噪声的均方根为2.44 mm。噪声在水平上的均方根为1.11 mm。垂直方向的噪声均方根是水平方向的两倍多,说明这种方法在垂直方向上影响比水平向影响更为明显;这也能说明GPS在垂直方向监测中包含的误差要比水平方向包含的误差大。

表7 G120站剔除非构造噪声结果 单位：mm

5 剔除前后结果对比分析

将剔除非构造噪声前后的153个区域站位移序列分别准备在两个文件里面,用Matlab编写程序实现批处理拟合区域站的线性速度[8-9],结果见表8。

表8 部分区域站剔除非构造噪声前后速度对比表 单位：mm/年

对剔除后的时间序列用加权最小二乘拟合速度[10],经过对比153个站发现(表10只列出部分):在南北方向上速度差最大为0.8 mm最小为0 mm,速度差均方根0.33 mm;在东西方向上速度差最大为1.0 mm,最小为0.1 mm,速度差均方根为0.55 mm;在垂直方向上速度差最大为-3.0 mm,最小为 0 mm,速度差均方根为1.71 mm;水平方向上的速度差均方根为0.58 mm;垂直方向上的速度差均方根是水平方向上的2.13倍;这个结论也说明在GPS在高程方向上包含的非构造噪声要大于水平方向上的。

通过Delaunay三角网和反距离内插算法相结合的处理,发现此方法能够从区域站中剔除部分非构造误差;同时,与QOCA剔除连续站的结果对比发现:前者在水平速度差均方根上为后者的2.0倍,在垂直速度差均方根上为后者的2.7倍,说明Delaunay三角网和反距离内插方法剔除区域站非构造噪声明显比QOCA软件剔除连续站的非构造噪声大两倍多。至于这两种方法存在差距的原因和这种新方法在区域站非构造噪声剔除中应用的可靠性还有待于深入研究,笔者将继续改变网形以及内插的方法,对比分析区域站和连续站位移时间序列中存在的关系。

6 结束语

使用QOCA软件能够提取GPS连续站的时间序列,借助连续站时间序列参数构建Delaunay三角网,结合反距离内插算法能够实现GPS流动观测区域站中非构造噪声的剔除;GPS区域站中非构造噪声的剔除,提高了地壳三维运动速度的精度,为地壳运动研究以及地震预报提供了基础数据。