利用 GPS阵列观测数据提取地震波信息的波束聚焦法＊

陆亚英 伍吉仓 孟国杰 胡丛玮

(1)同济大学测量与国土信息工程系,上海 200092 2)中国地震局地震预测研究所,北京 100036)

利用 GPS阵列观测数据提取地震波信息的波束聚焦法＊

陆亚英1)伍吉仓1)孟国杰2)胡丛玮1)

(1)同济大学测量与国土信息工程系,上海 200092 2)中国地震局地震预测研究所,北京 100036)

介绍一种利用 GPS阵列高频观测数据提取地震面波信息的波束聚焦法。通过在 GPS坐标序列的垂直分量中加入模拟的瑞利波信号,得到 GPS阵列观测站模拟地震波记录,分析了波束聚焦法在不同信噪比下提取地震波信号的能力。结果表明,当信噪比为 1.5时,不能获得正确的慢度估值;当信噪比为 3时,提取的慢度估值接近理论值。

波束聚焦;GPS;高频数据;地震波;慢度

1 概述

波束聚焦法是 20世纪 60年代发展起来的地震波成像方法[1,2]。波束聚焦最初利用地震台站阵列相对于震源的距离和方向不同,把每个台站记录到的地震波形在时间轴上做相应的平移后再叠加,从而使相干波源得到增强,便于提取。近年来,随着GPS记录的频率不断提高,GPS观测得到的坐标序列同附近地震台站数字地震仪记录得到的波形非常接近,进而有人提出了 GPS地震仪[3-5]的方法。GPS地震仪相对于传统的数字地震仪有不受振幅限制,直接测量位移等优势。但采用 GPS观测时,由于多路径影响、观测噪音以及参考站稳定性等都有可能使得表征位移的坐标时间序列产生误差,不能真实地记录地震波形。考虑到来自震源的面波的相干性,它们到达各个 GPS站的时间因路径长度不同而存在差异,这就为采用波束聚焦法提取地震波形提供了可能性。因此,将通过在 GPS实际观测数据中加入模拟地震波信号,研究波束聚焦法提取地震波信息的能力。

2 基于 GPS阵列提取地震波慢度信息

2.1 基于 GPS阵列的波束聚焦法

假定从震源传来的地震波是一个平面波,波速为 v,其波束相对于 GPS测站阵列所在水平面的入射角为 i,波束在水平面上的投影的后视方位角为θ (图 1)。图 1中 r为 GPS站的位置矢量,u为平面波的慢度矢量,(ux,uy,uz)为其在 N、E、U 3个方向的分量,可以由公式 (1)计算得到,u//为慢度矢量在水平面内的投影。

图1 波束慢度矢量和 GPS站点位置矢量Fig.1 Beam slowness vector and position vector of GPSstation

一般我们选定 GPS网中心附近的站为参考站,假定到达参考站的波为：

式中,f(t)为待求的地震波信号,n0(t)为噪声信号,包含多路径、对流层影响以及观测噪声等。考虑到其他 GPS站点记录的波形相同只是相位不同,起因于波束到达该观测站相对于到达参考站的时间差,假定坐标原点位于参考站,则位置向量为 ri的 GPS站点记录的波形应为：

其中,时间差为:

如果不考虑参考站,并将每个 GPS站点记录得到的波形平移一个时间,得到与参考站同相位的波形,记为：

对式(5)相加取平均,这样慢度为 u//的波束得到增强,从而抑制了噪声和其他慢度的波束。波束聚焦后得到的波形为：

式中,N表示 GPS站点的个数,˜x0(t)为参考站记录的波形数据,即为 x0(t)。则慢度 u对应的聚焦波束能量为:

其中,tj为 GPS历元 j所对应的时刻,M为历元个数。实际利用波束聚焦法提取地震波信号时,在慢度分量 ux、uy的取值区间内逐个计算其所对应的能量值,按式(4)～(7)依次计算时间平移量和聚焦波束能量,把能量最大对应的慢度信号提取出来,这个信号就是地震波的面波信号。

2.2 利用波束聚焦法提取慢度信息

理想情况下,正确的地震波慢度对应聚焦后能量最大,但由于 GPS数据中存在误差,使从能量分布图中得到的慢度与理论值之间存在偏差,聚焦后能量最大点所对应的值不一定就是地震波的理论慢度,最大值及其周围能量比较大的若干点,都有可能是所要求的慢度。为此我们以能量最大点为中心取一个窗口,将窗口内所有能量大于某一阈值的点提取出来,利用这些点计算地震波的慢度。则地震波的慢度为:

其中,K为所提取的点的个数,Pi为第 i个点的权,记为:

则慢度的精度可以表示为:

3 GPS高频数据处理方法

GPS高频数据处理一般选择一个测站作为参考站,将其位置约束到某一个参考框架上,然后解算其他测站的坐标[3]。本文采用天津市某 GPS连续观测阵列在 2005年,年积日为 357的数据,取其中 8个站进行分析,其点位的分布如图 2所示,数据的采样频率为 1 Hz。选取 YC01为参考点,用 GAM IT软件中的 TRACK运动学模块[6-8],计算得到 GPS阵列相对参考站的坐标序列。

图2 GPS观测站的分布Fig.2 Distribution of GPS stations

3.1 粗差剔除与补齐

由于外界条件、硬件因素或软件原因,所得到的坐标序列中存在粗差,应予以剔除,我们采用 3σ准则来探测并剔除粗差:

其中,xi为某一方向的坐标,为该方向坐标序列的平均值,σ为其标准差。此时得到的坐标序列是不连续的,因此还需要将其补齐,文中采用线性插值方法进行处理。

3.2 公共误差的处理

GPS站点相对于参考站的坐标序列中都含有由参考站所引起的公共误差[9,10],各条基线的误差可表示为：

其中,Δi是第 i个 GPS观测站上的随机误差,Δ0是由参考站所引起的误差。把所有基线对应方向的坐标序列相加,得到其对应方向的误差总和为:

表 1为削弱公共误差前后均方根的对比,从表1中可以看出,削弱公共误差后,3个方向的精度都有明显的提高。

表 1 削弱公共误差前后均方根的比较(单位:mm)Tab.1 Comparison of RM S between before and after weaken i ng common errors(un it:mm)

4 模拟实验分析

4.1 模拟地震波

瑞利波是地震面波的一种,在地球内部的传播速度为 2.0～4.2 km/s,质点的运动垂直于传播方向。1998年 9月 3日,智利中海岸发生了 6.5级地震,本文利用地震台站NNA观测到的地震波数据进行模拟实验,瑞利波在该站持续的时间大约为160 s。模拟计算时把 YC01站视作NNA站,将瑞利波数据加到 GPS解算得到的垂直分量坐标序列中,模拟得到包含地震波信息的坐标序列。

具体模拟过程如下：首先取表 1中列出的 7条基线U分量 2小时的坐标序列,经过粗差剔除、补齐和削弱公共误差后,去掉 U分量的均值,作为地震波信号记录的背景噪声。假定地震波从序列开始后 15分钟到达参考站,记为 t=900 s,瑞利波的传播速度 v=3 km/s,入射角 i=30°,后视方位角θ= 30°。通过式 (1)计算得到的慢度分量理论值为 ux=0.14 s/km,uy=0.08 s/km。考虑到每个 GPS站点相对于参考站的位置不同,地震波到达的时间不同,在站的U分量坐标序列中加入地震波的时候要考虑时间差,时间差Δti可由公式 (4)计算得到。另外取瑞利波振幅为 7条基线 U分量坐标序列中误差平均值的3倍,加入 GPS坐标序列中,得到的 GPS模拟地震波记录的时间序列如图 3所示。图 4为包含地震波信号的局部放大图,从图 4可以看出,由于各 GPS站点相对于参考站的位置不同,地震波到达GPS测站的时间是不同的。

4.2 慢度信息的提取

为了研究该方法在不同信噪比的情况下提取地震波信号的能力,我们分别试验了以 3种不同大小振幅的瑞利波加入到 GPS坐标序列的情况。3种情况的信噪比分别为 1.5、2和 3(这里信噪比是指地震波的振幅与所有 GPS的U分量坐标序列中误差平均值的比值)。

在准备好模拟 GPS地震波观测时间序列后,首先要给水平方向的两个慢度分量一个搜索范围。本文以已知的慢度值为中心,慢度变化区间取(-0.3,0.3),循环间隔为 0.05,利用公式 (4)依次计算每次循环中每个站点相对参考站的平移时间Δt,然后把每个序列平移相应的时间,再由式 (7)得到每次波束聚焦后的波形,然后计算每次循环相应慢度对应的能量 d。

图 5为不同信噪比情况下每次实验得到的能量图,其中图 5(a)是信噪比为 0时,即没有地震波信号时的能量图,图 5(b、c、d)分别是信噪比为 1.5、2、3 3种情况的能量图。从图 5可以看出,当信噪比达到一定的数值时,我们在能量图中可以很明显地发现有一块区域的能量强于其他地区,在这个区域中有地震波的存在。随着信噪比的增强,这个地震波信息会反应得更加明显。

表 2列出了 3种信噪比下提取出的地震波慢度值及其中误差,其中慢度理论值 ux=0.14 s/km,uy=0.08 s/km。从表 2可以看出,随着信噪比的增大,慢度的精度逐渐提高。当信噪比为 1.5时,不能获得有效的慢度信息,当信噪比为 3时提取的慢度信息接近理论值。

图3 7条基线U分量的坐标序列Fig.3 Coordinate series of seven baselines in U direction

图4 坐标序列局部Fig.4 Partial graph of coordinate series

图5 不同信噪比的能量图Fig.5 Energy diagrams of different SNR

表 2 不同信噪比下提取的慢度信息Tab.2 Slowness information extracted from energy di agram s of different SNR

5 结论

研究了利用 GPS高频数据提取地震波信息的波束聚焦法。用波束聚焦法处理天津市某连续观测GPS阵列 8个站点 1 Hz的 GPS模拟地震波观测数据,提取了地震波的慢度信息。研究表明,随着信噪比的增大,地震波信号在能量图中表现越明显,提取的慢度信息越准确。当信噪比足够大时,可以精确的提取出慢度信息,但当信噪比较小时,该方法不能获得准确的慢度信息。

1 Rost S and Thomas C.Array seismology：Methods and applications[J].Rev.Geophys.,2002,40(3),1008,doi：10.1029/2000RG000100.

2 Davis J P and Smalley R.Love wave dispersion in central North America deter mined using absolute displacement seismograms from high-rate GPS[J]. J. Geophys. Res., 2009,114,B11303,doi：10.1029/2009JB006288.

3 孟国杰,等.GPS高频数据处理方法及其在地震学中的应用研究进展[J].国际地震动态,2007,343(7)：26-31. (Men Guojie,et al.Data processing methods of high rate GPS and its application to seismology[J].Recent development inWord Seis mology,2007,343(7)：26-31)

4 殷海涛,等.利用高频 GPS技术进行强震地面运动监测的研究进展[J].地球物理学进展,2009,24(6)：2 012 -2 019.(Yin Haotao,et al.Progress on monitoring strong earthquake ground motions using high-rate GPS[J].Progress in Geophys,2009,24(6)：2 012-2 019)

5 Nikolaidis R,et al.Seis mic wave observationwith the Global Positioning System[J].J.Geophys.Res.,2001,106：218 097-21 916.

6 Gang Chen.GPS kinematic positioning for the airborne laser alti metry at Long Valley,California[D].Mass. Inst. of Tech.,USA,1998.

7 King R W and Bock Y.Documentation for the GAM IT GPS analysis software[M].Mass.Inst.of Tech.,Scripps Inst. Occeangr.,2003.

8 苏小宁,等.基于 TRACK进行 GPS单历元定位[J].大地测量与地球动力学,2009,(3)：100-103.(Su Xiaoning, et al.Single epoch GPS positioning based on track module [J].Journal of Geodesy and Geodynamics,2009(3)：100-103)

9 伍吉仓,孙亚峰,刘朝功.连续 GPS站坐标序列共性误差的提取与形变分析[J].大地测量与地球动力学,2008, (4)：97-101.(Wu Jichang,Sun Yafeng and Liu Chaogong.Extraction of common mode errors for continuous GPS net works and deformation analysis[J].Journal of Geodesy and Geodynamics,2008,(4)：97-101)

10 胡守超,伍吉仓,孙亚峰.区域 GPS网 3种时空滤波方法的比较[J].大地测量与地球动力学,2009,(3)：95-99.(Hu Shouchao,Wu Jicang and Sun Yafeng.Comparison among three spatiotemporal filtering methods for regional GPS networks analysis[J].Journal of Geodesy and Geodynamics,2009,(3)：95-99)

M ETHOD OF BEAM FORM ING FOR EXTRACTING SEISM IC WAVE INFORMATION BY USING GPS ARRAY DATA

Lu Yaying1),Wu Jicang1),Meng Guojie2)and Hu Congwei1)

(1)Depart m ent of Surveying and Geo-infor m atics,Tongji University,Shanghai 200092 2)Institute of Earthquake Science,CEA,B eijing 100036)

A method of beam for ming forprocessing high-rate GPS data and extracting seis mic wave information is introduced.Simulated seismicwavesof each GPS station are synthesized by adding simulated Rayleighwaves into vertical component of GPS coordinates series.The capacity of beam for ming to extract seis mic wave infor mation of different SNR(signal to noise ratio)is analyzed.The results indicate that,the correct value of slowness cannot be obtained when the SNR is 1.5,while it is close to the theoretical value when the SNR ratio is 3.

beam for ming;GPS;high-rate data;seis mic wave;slowness

1671-5942(2010)05-0068-05

2010-03-10

公益性行业科研专项(200708030);国家自然科学基金(40674004,40671155)

陆亚英,女,1986年生,硕士研究生,研究方向为 GPS数据处理.E-mail：yayinglu1986@163.com

P315.61

A