陸亞英 伍吉倉 孟國杰 胡叢瑋

(1)同濟(jì)大學(xué)測量與國土信息工程系,上海 200092 2)中國地震局地震預(yù)測研究所,北京 100036)

利用 GPS陣列觀測數(shù)據(jù)提取地震波信息的波束聚焦法＊

陸亞英1)伍吉倉1)孟國杰2)胡叢瑋1)

(1)同濟(jì)大學(xué)測量與國土信息工程系,上海 200092 2)中國地震局地震預(yù)測研究所,北京 100036)

介紹一種利用 GPS陣列高頻觀測數(shù)據(jù)提取地震面波信息的波束聚焦法。通過在 GPS坐標(biāo)序列的垂直分量中加入模擬的瑞利波信號(hào),得到 GPS陣列觀測站模擬地震波記錄,分析了波束聚焦法在不同信噪比下提取地震波信號(hào)的能力。結(jié)果表明,當(dāng)信噪比為 1.5時(shí),不能獲得正確的慢度估值;當(dāng)信噪比為 3時(shí),提取的慢度估值接近理論值。

波束聚焦;GPS;高頻數(shù)據(jù);地震波;慢度

1 概述

波束聚焦法是 20世紀(jì) 60年代發(fā)展起來的地震波成像方法[1,2]。波束聚焦最初利用地震臺(tái)站陣列相對(duì)于震源的距離和方向不同,把每個(gè)臺(tái)站記錄到的地震波形在時(shí)間軸上做相應(yīng)的平移后再疊加,從而使相干波源得到增強(qiáng),便于提取。近年來,隨著GPS記錄的頻率不斷提高,GPS觀測得到的坐標(biāo)序列同附近地震臺(tái)站數(shù)字地震儀記錄得到的波形非常接近,進(jìn)而有人提出了 GPS地震儀[3-5]的方法。GPS地震儀相對(duì)于傳統(tǒng)的數(shù)字地震儀有不受振幅限制,直接測量位移等優(yōu)勢(shì)。但采用 GPS觀測時(shí),由于多路徑影響、觀測噪音以及參考站穩(wěn)定性等都有可能使得表征位移的坐標(biāo)時(shí)間序列產(chǎn)生誤差,不能真實(shí)地記錄地震波形?？紤]到來自震源的面波的相干性,它們到達(dá)各個(gè) GPS站的時(shí)間因路徑長度不同而存在差異,這就為采用波束聚焦法提取地震波形提供了可能性。因此,將通過在 GPS實(shí)際觀測數(shù)據(jù)中加入模擬地震波信號(hào),研究波束聚焦法提取地震波信息的能力。

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

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

假定從震源傳來的地震波是一個(gè)平面波,波速為 v,其波束相對(duì)于 GPS測站陣列所在水平面的入射角為 i,波束在水平面上的投影的后視方位角為θ (圖 1)。圖 1中 r為 GPS站的位置矢量,u為平面波的慢度矢量,(ux,uy,uz)為其在 N、E、U 3個(gè)方向的分量,可以由公式 (1)計(jì)算得到,u//為慢度矢量在水平面內(nèi)的投影。

圖1 波束慢度矢量和 GPS站點(diǎn)位置矢量Fig.1 Beam slowness vector and position vector of GPSstation

一般我們選定 GPS網(wǎng)中心附近的站為參考站,假定到達(dá)參考站的波為：

式中,f(t)為待求的地震波信號(hào),n0(t)為噪聲信號(hào),包含多路徑、對(duì)流層影響以及觀測噪聲等?？紤]到其他 GPS站點(diǎn)記錄的波形相同只是相位不同,起因于波束到達(dá)該觀測站相對(duì)于到達(dá)參考站的時(shí)間差,假定坐標(biāo)原點(diǎn)位于參考站,則位置向量為 ri的 GPS站點(diǎn)記錄的波形應(yīng)為：

其中,時(shí)間差為:

如果不考慮參考站,并將每個(gè) GPS站點(diǎn)記錄得到的波形平移一個(gè)時(shí)間,得到與參考站同相位的波形,記為：

對(duì)式(5)相加取平均,這樣慢度為 u//的波束得到增強(qiáng),從而抑制了噪聲和其他慢度的波束。波束聚焦后得到的波形為：

式中,N表示 GPS站點(diǎn)的個(gè)數(shù),?x0(t)為參考站記錄的波形數(shù)據(jù),即為 x0(t)。則慢度 u對(duì)應(yīng)的聚焦波束能量為:

其中,tj為 GPS歷元 j所對(duì)應(yīng)的時(shí)刻,M為歷元個(gè)數(shù)。實(shí)際利用波束聚焦法提取地震波信號(hào)時(shí),在慢度分量 ux、uy的取值區(qū)間內(nèi)逐個(gè)計(jì)算其所對(duì)應(yīng)的能量值,按式(4)～(7)依次計(jì)算時(shí)間平移量和聚焦波束能量,把能量最大對(duì)應(yīng)的慢度信號(hào)提取出來,這個(gè)信號(hào)就是地震波的面波信號(hào)。

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

理想情況下,正確的地震波慢度對(duì)應(yīng)聚焦后能量最大,但由于 GPS數(shù)據(jù)中存在誤差,使從能量分布圖中得到的慢度與理論值之間存在偏差,聚焦后能量最大點(diǎn)所對(duì)應(yīng)的值不一定就是地震波的理論慢度,最大值及其周圍能量比較大的若干點(diǎn),都有可能是所要求的慢度。為此我們以能量最大點(diǎn)為中心取一個(gè)窗口,將窗口內(nèi)所有能量大于某一閾值的點(diǎn)提取出來,利用這些點(diǎn)計(jì)算地震波的慢度。則地震波的慢度為:

其中,K為所提取的點(diǎn)的個(gè)數(shù),Pi為第 i個(gè)點(diǎn)的權(quán),記為:

則慢度的精度可以表示為:

3 GPS高頻數(shù)據(jù)處理方法

GPS高頻數(shù)據(jù)處理一般選擇一個(gè)測站作為參考站,將其位置約束到某一個(gè)參考框架上,然后解算其他測站的坐標(biāo)[3]。本文采用天津市某 GPS連續(xù)觀測陣列在 2005年,年積日為 357的數(shù)據(jù),取其中 8個(gè)站進(jìn)行分析,其點(diǎn)位的分布如圖 2所示,數(shù)據(jù)的采樣頻率為 1 Hz。選取 YC01為參考點(diǎn),用 GAM IT軟件中的 TRACK運(yùn)動(dòng)學(xué)模塊[6-8],計(jì)算得到 GPS陣列相對(duì)參考站的坐標(biāo)序列。

圖2 GPS觀測站的分布Fig.2 Distribution of GPS stations

3.1 粗差剔除與補(bǔ)齊

由于外界條件、硬件因素或軟件原因,所得到的坐標(biāo)序列中存在粗差,應(yīng)予以剔除,我們采用 3σ準(zhǔn)則來探測并剔除粗差:

其中,xi為某一方向的坐標(biāo),為該方向坐標(biāo)序列的平均值,σ為其標(biāo)準(zhǔn)差。此時(shí)得到的坐標(biāo)序列是不連續(xù)的,因此還需要將其補(bǔ)齊,文中采用線性插值方法進(jìn)行處理。

3.2 公共誤差的處理

GPS站點(diǎn)相對(duì)于參考站的坐標(biāo)序列中都含有由參考站所引起的公共誤差[9,10],各條基線的誤差可表示為：

其中,Δi是第 i個(gè) GPS觀測站上的隨機(jī)誤差,Δ0是由參考站所引起的誤差。把所有基線對(duì)應(yīng)方向的坐標(biāo)序列相加,得到其對(duì)應(yīng)方向的誤差總和為:

表 1為削弱公共誤差前后均方根的對(duì)比,從表1中可以看出,削弱公共誤差后,3個(gè)方向的精度都有明顯的提高。

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

4 模擬實(shí)驗(yàn)分析

4.1 模擬地震波

瑞利波是地震面波的一種,在地球內(nèi)部的傳播速度為 2.0～4.2 km/s,質(zhì)點(diǎn)的運(yùn)動(dòng)垂直于傳播方向。1998年 9月 3日,智利中海岸發(fā)生了 6.5級(jí)地震,本文利用地震臺(tái)站NNA觀測到的地震波數(shù)據(jù)進(jìn)行模擬實(shí)驗(yàn),瑞利波在該站持續(xù)的時(shí)間大約為160 s。模擬計(jì)算時(shí)把 YC01站視作NNA站,將瑞利波數(shù)據(jù)加到 GPS解算得到的垂直分量坐標(biāo)序列中,模擬得到包含地震波信息的坐標(biāo)序列。

具體模擬過程如下：首先取表 1中列出的 7條基線U分量 2小時(shí)的坐標(biāo)序列,經(jīng)過粗差剔除、補(bǔ)齊和削弱公共誤差后,去掉 U分量的均值,作為地震波信號(hào)記錄的背景噪聲。假定地震波從序列開始后 15分鐘到達(dá)參考站,記為 t=900 s,瑞利波的傳播速度 v=3 km/s,入射角 i=30°,后視方位角θ= 30°。通過式 (1)計(jì)算得到的慢度分量理論值為 ux=0.14 s/km,uy=0.08 s/km。考慮到每個(gè) GPS站點(diǎn)相對(duì)于參考站的位置不同,地震波到達(dá)的時(shí)間不同,在站的U分量坐標(biāo)序列中加入地震波的時(shí)候要考慮時(shí)間差,時(shí)間差Δti可由公式 (4)計(jì)算得到。另外取瑞利波振幅為 7條基線 U分量坐標(biāo)序列中誤差平均值的3倍,加入 GPS坐標(biāo)序列中,得到的 GPS模擬地震波記錄的時(shí)間序列如圖 3所示。圖 4為包含地震波信號(hào)的局部放大圖,從圖 4可以看出,由于各 GPS站點(diǎn)相對(duì)于參考站的位置不同,地震波到達(dá)GPS測站的時(shí)間是不同的。

4.2 慢度信息的提取

為了研究該方法在不同信噪比的情況下提取地震波信號(hào)的能力,我們分別試驗(yàn)了以 3種不同大小振幅的瑞利波加入到 GPS坐標(biāo)序列的情況。3種情況的信噪比分別為 1.5、2和 3(這里信噪比是指地震波的振幅與所有 GPS的U分量坐標(biāo)序列中誤差平均值的比值)。

在準(zhǔn)備好模擬 GPS地震波觀測時(shí)間序列后,首先要給水平方向的兩個(gè)慢度分量一個(gè)搜索范圍。本文以已知的慢度值為中心,慢度變化區(qū)間取(-0.3,0.3),循環(huán)間隔為 0.05,利用公式 (4)依次計(jì)算每次循環(huán)中每個(gè)站點(diǎn)相對(duì)參考站的平移時(shí)間Δt,然后把每個(gè)序列平移相應(yīng)的時(shí)間,再由式 (7)得到每次波束聚焦后的波形,然后計(jì)算每次循環(huán)相應(yīng)慢度對(duì)應(yīng)的能量 d。

圖 5為不同信噪比情況下每次實(shí)驗(yàn)得到的能量圖,其中圖 5(a)是信噪比為 0時(shí),即沒有地震波信號(hào)時(shí)的能量圖,圖 5(b、c、d)分別是信噪比為 1.5、2、3 3種情況的能量圖。從圖 5可以看出,當(dāng)信噪比達(dá)到一定的數(shù)值時(shí),我們?cè)谀芰繄D中可以很明顯地發(fā)現(xiàn)有一塊區(qū)域的能量強(qiáng)于其他地區(qū),在這個(gè)區(qū)域中有地震波的存在。隨著信噪比的增強(qiáng),這個(gè)地震波信息會(huì)反應(yīng)得更加明顯。

表 2列出了 3種信噪比下提取出的地震波慢度值及其中誤差,其中慢度理論值 ux=0.14 s/km,uy=0.08 s/km。從表 2可以看出,隨著信噪比的增大,慢度的精度逐漸提高。當(dāng)信噪比為 1.5時(shí),不能獲得有效的慢度信息,當(dāng)信噪比為 3時(shí)提取的慢度信息接近理論值。

圖3 7條基線U分量的坐標(biāo)序列Fig.3 Coordinate series of seven baselines in U direction

圖4 坐標(biāo)序列局部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 結(jié)論

研究了利用 GPS高頻數(shù)據(jù)提取地震波信息的波束聚焦法。用波束聚焦法處理天津市某連續(xù)觀測GPS陣列 8個(gè)站點(diǎn) 1 Hz的 GPS模擬地震波觀測數(shù)據(jù),提取了地震波的慢度信息。研究表明,隨著信噪比的增大,地震波信號(hào)在能量圖中表現(xiàn)越明顯,提取的慢度信息越準(zhǔn)確。當(dāng)信噪比足夠大時(shí),可以精確的提取出慢度信息,但當(dāng)信噪比較小時(shí),該方法不能獲得準(zhǔn)確的慢度信息。

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高頻數(shù)據(jù)處理方法及其在地震學(xué)中的應(yīng)用研究進(jìn)展[J].國際地震動(dòng)態(tài),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技術(shù)進(jìn)行強(qiáng)震地面運(yùn)動(dòng)監(jiān)測的研究進(jìn)展[J].地球物理學(xué)進(jìn)展,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進(jìn)行 GPS單歷元定位[J].大地測量與地球動(dòng)力學(xué),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 伍吉倉,孫亞峰,劉朝功.連續(xù) GPS站坐標(biāo)序列共性誤差的提取與形變分析[J].大地測量與地球動(dòng)力學(xué),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 胡守超,伍吉倉,孫亞峰.區(qū)域 GPS網(wǎng) 3種時(shí)空濾波方法的比較[J].大地測量與地球動(dòng)力學(xué),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

公益性行業(yè)科研專項(xiàng)(200708030);國家自然科學(xué)基金(40674004,40671155)

陸亞英,女,1986年生,碩士研究生,研究方向?yàn)?GPS數(shù)據(jù)處理.E-mail：yayinglu1986@163.com

P315.61

A