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

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

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

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

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

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

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

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

1 概述

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

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

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

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

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

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

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

其中,時間差為:

如果不考慮參考站,并將每個 GPS站點記錄得到的波形平移一個時間,得到與參考站同相位的波形,記為：

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

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

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

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

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

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

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

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

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

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

3.1 粗差剔除與補齊

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

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

3.2 公共誤差的處理

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

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

表 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 模擬地震波

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

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

4.2 慢度信息的提取

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

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

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

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

圖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個站點 1 Hz的 GPS模擬地震波觀測數(shù)據(jù),提取了地震波的慢度信息。研究表明,隨著信噪比的增大,地震波信號在能量圖中表現(xiàn)越明顯,提取的慢度信息越準(zhǔn)確。當(dāng)信噪比足夠大時,可以精確的提取出慢度信息,但當(dāng)信噪比較小時,該方法不能獲得準(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].國際地震動態(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)行強震地面運動監(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].大地測量與地球動力學(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].大地測量與地球動力學(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種時空濾波方法的比較[J].大地測量與地球動力學(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è)科研專項(200708030);國家自然科學(xué)基金(40674004,40671155)

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

P315.61

A