2022年用多道反褶积方法测定台站接收函数 .pdf

《2022年用多道反褶积方法测定台站接收函数 .pdf》由会员分享，可在线阅读，更多相关《2022年用多道反褶积方法测定台站接收函数 .pdf（6页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/吴庆举,李永华,张瑞青等.用多道反褶积方法测定台站接收函数.地球物理学报,2007,50(3):791796Wu Q J,Li Y H,Zhang R Q,et al.Receiver function estimated by multi2channel deconv olution.ChineseJ.Geophys.(in Chinese),2007,50(3):791796用多道反褶积方法测定台站接收函数吴庆举,李

2、永华,张瑞青,张乃铃中国地震局地球物理研究所,北京 100081摘 要 本文提出了一种在时间域用多道反褶积测定台站接收函数的方法,以提高台站接收函数的测量精度与分辨率.在单道反褶积的基础上,选取若干个质量较好的远震P波波形事件,构成多道信号,以垂直分量为输入,径向和切向分量作为期望输出,依据最小二乘,设计多道滤波器,提取接收函数.合成地震图与观测地震图的检验表明,多道滤波方法能有效地测定台站接收函数,特别是多道反褶积能够有效地恢复地壳上地幔间断面所产生的弱转换波震相.关键词 接收函数,反褶积,多道滤波器文章编号 0001-5733(2007)03-0791-06中图分类号 P315收稿日期 2

3、005-11-28,2006-12-09收修定稿基金项目 国家自然科学基金(40274029)资助.中国地震局地球物理研究所论著号:06AC1012.作者简介 吴庆举,男,1966年生,固体地球物理学博士,研究员,主要从事地震学、地球内部结构和地球动力学研究.E2mail:wuqj Receiver function estimated by multi2channel deconvolutionWU Qing2Ju,LIYong2Hua,ZHANG Rui2Qing,ZHANG Nai2LingInstitute of Geophysics,China Earthquake Administ

4、ration,Beijing100081,ChinaAbstractTime2domain multi2channel deconv olution is put forward to estimate receiver function,to improveresolution.Based on the deconv olution results of individual teleseismic P waveforms,a number of eventswithgood quality are selected to form multi2channel signals,the ver

5、tical components are regarded as inputs,andradial and transversecomponents are taken as expectedoutputs.Least2square error is usedto designthe multi2channel filter,to get the comm on filterfactors,i.e.,receiver function.Bothsynthetic and real dataexperiments show that multi2channel deconv olution is

6、 an effective approach to measure receiver function,especially,it can recover the weak converted phasesfrom the upper mantle.KeywordsReceiverfunction,Deconv olution,Multi2channel filter1 引 言接收函数16是用远震P波波形的垂直分量对径向分量和切向分量作反褶积处理后得到的时间序列,被认为基本消除了震源时间函数和深地幔传播路径效应的影响,可以近似为台站下方地壳上地幔速度结构对近垂直入射平面波的脉冲响应.为压制噪音

7、,提高接收函数的信噪比,通常要对同一方向、震中距变化较小的多个远震事件的接收函数进行叠加,得到该方向震中距对的一个接收函数的平均估计.考察台站接收函数的方位性变化特征,可以研究台站下方地壳上地幔速度结构的横向变化以及各向异性.接收函数实质上是近垂直入射的脉冲平面波在台站下方速度界面产生的PS转换波及地表与速度界面之间的多次反射 转换波.与其他地震学研究方法相比,接收函数具有较小的横向采样范围,受介质第50卷 第3期2007年5月地 球 物 理 学 报CHINESEJOURNALOFGEOPHYSICSVol.50,No.3May,2007名师资料总结-精品资料欢迎下载-名师精心整理-第 1 页

8、，共 6 页 -?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/横向非均匀性的影响相对较小,具有较高的横向分辨率,特别是接收函数对S波速度结构非常敏感,而S波速度结构是其他地震学方法难以很好地加以约束的重要参数.因此,接收函数成为观测台站下方地壳上地幔 S波速度结构最为理想的研究方法之一,被广泛地应用于地壳上地幔间断面的探测与研究16.目前,接收函数大致沿着两个相互平行、相对独立的研究方向发展:一是分层结构研究,旨在揭示地壳及上地幔顶部S波速度的垂向变化7、倾斜界面的

9、产状8、介质的各向异性特征9、地壳泊松比及地壳厚度等10;二是偏移成像研究,试图揭示地壳上地幔间断面的横向变化1116.测定精确的台站接收函数是接收函数研究的首要步骤.传统的接收函数测定是在频率域用频谱相除并反变换回时间域得到.由于实际的远震资料是有限带宽的,且含有各种噪音,垂直分量一般都存在近零的频谱成分,导致频率域反褶积不够稳定,故要对垂直分量作预白化处理,压制其近零值的频谱成分.但与此同时,又不可避免地压制了能量较弱的转换波震相,特别是上地幔间断面的微弱震相.为此,人们提出了多种反褶积方法,试图提高台站接收函数的测定精度和分辨率.如最大或然性估计2、Wiener滤波17,18、最大熵谱反

10、褶积19、多重镶窗20、迭代求解21 等.本文提出了一种用多道滤波方法测定台站接收函数的方法,试图以多个远震事件的垂直分量作为输入,水平分量作为输出,设计多道滤波因子,来测定台站接收函数.2 接收函数和反褶积接收函数是用三分量远震P波波形的垂直分量对径向和切向分量作如下的数学运算后得到的时间序列18,19,R(t)=12DH()DV()eitd,(1)其中,R(t)为接收函数,DH()和DV()分别是水平分量和垂直分量的频谱,是圆频率,t是时间.上述数学运算过程称之为反褶积,所得到的接收函数可以看成是台站下方接收介质的近垂直入射平面波的脉冲响应,而与震源时间函数及路径传播效应基本无关.从三分量

11、远震P波波形中分离出台站接收函数是接收函数研究的关键.人们提出了多种方法来提高接收函数的测量精度和分辨率2,1721,可以概括为频率域和时间域两类方法.频率域算法通过改善频谱估计的精度和分辨率来测定接收函数,时间域算法则依据一定的准则,设计相应的滤波器来计算接收函数.我们则提出了一种多道滤波的方法,通过设计多道滤波器,在时间域提取接收函数.设计多道滤波器的关键在于滤波器法则的选取.不同的法则决定了不同的滤波器,如最小二乘滤波器、最大似然滤波器、最大熵滤波器等.我们将以最小二乘滤波器为例,也就是通常所说的 Wiener滤波,介绍多道滤波的基本原理,并可以非常自然地推广应用于其他形式的多道滤波器.

12、3多道滤波台站接收函数一般是同一方向、震中距变化不大的多个远震事件接收函数的叠加平均,构成某一方向 震中距对的接收函数平均估计,以压制噪音,提高接收函数的信噪比.既然多个事件接收函数的叠加平均可以看成是某一方向 震中距对的接收函数,那么,我们可以认为那些方位角震中距变化不大的远震事件具有非常相似的接收函数,可以用一个接收函数代替.再者,无论台站下方地壳上地幔速度结构多么复杂,对于方位角震中距变化不大的若干远震事件而言,理论上,它们应具有非常相似的接收响应.因此,在提取台站接收函数时,我们可以直接从多个远震事件中提取一个共同的接收响应,也就是接收函数,而不必先一一求取单个远震事件的接收函数,再通

13、过叠加平均的方法,获取一个接收函数的平均估计.我们不妨将前一种方法称之为接收函数同步测定法,而将后一种方法称之为接收函数异步测定法.较之于先提取后叠加的接收函数异步测定法,接收函数的同步测定法可以在接收函数的计算过程中:从若干个看似不相干的远震波形数据中,根据波形的相干性,消除震源时间函数及路径效应,并压制噪音,保留并增强相干信号,得到待求的接收函数.我们以最小二乘为例,介绍一下多道滤波的基本思想,这一思路可以自然地推广到其他形式的多道滤波器.基于最小二乘法则的滤波实质上就是著名的 Wiener滤波.我们首先简单介绍一下单道Wiener滤波.用向297地 球 物 理 学 报(ChineseJ.

14、Geophys.)50卷 名师资料总结-精品资料欢迎下载-名师精心整理-第 2 页，共 6 页 -?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/量xit表示第i个远震 P波波形垂直分量的时间序列,yit表示第i个远震 P波波形水平分量的时间序列,fit表示第i个接收函数的时间序列,在频率域分别用Xi(),Yi()和Fi()表示.根据接收函数定义可知,远震 P波波形的水平分量和垂直分量以及接收函数这三个时间序列在频率域可表示为Yi()=Xi()Fi();(2)在时间域

15、可表示成如下的褶积形式:yit=fit3xit.(3)显然可以将接收函数时间序列视为滤波器的单位脉冲响应,而将远震P波波形的垂直分量和水平分量分别当成滤波器的输入和输出,但事实上不可能得到这样的滤波器,即在远震P波波形的垂直分量的激励下,其输出与远震P 波波形的水平分量完全一致,而只能得到它的近似估计,即yit=fit3xit.(4)我们期望在xit的激励下滤波器fit的实际输出yit与期望输出yit尽可能地相似.用eit代表实际输出与期望输出之间的误差向量:eit=yit-yit=yit-fit3xit.(5)以误差向量的l2模取极小作为滤波器的原则,就可得到所谓的Wiener滤波器,误差向

16、量的l2范数可表示成Ei=t(eit)2.(6)对(6)式求极小,可得到如下的矩阵方程:a0a1a2ana1a0a1an-1anan-1an-2a0fi0fi1fin=c0c1cn,(7)其中fi为待定的第i个接收函数时间序列,a为第i个远震 P波波形垂直分量的自相关系数,c为第i个远震 P 波波形垂直分量与水平分量的互相关系数,即a=nt=0 xitxit+,c=nt=0 xityit+.(8)自相关系数矩阵是Toeplitz 矩阵,由于其对称性,可采用 Levinson 递推算法求解18.(7)式可以简化为如下的表达式:Aifi=Ci.(9)对于多个远震事件,接收函数自身独立于地震事件,换

17、言之,所有远震事件均具有非常相似的接收函数,而与地震事件无关.有两种设计多道滤波器的方法,分别加以介绍.(1)将(7)式或(9)式自然扩展成(MN)N维的超定矩阵方程,其中,M为参与运算的远震事件数,N为待求接收函数的数据长度,表示如下:Af=C,(10)其中,矩阵A为形如(7)式或(9)式的M个远震事件垂直分量的自相关系数所组成的(MN)N维的矩阵,C则为M个远震事件垂直分量与水平分量的互相关系数所组成的长度为MN的向量,f为多道滤波所得到的N维向量.超定矩阵方程(10)式的求解可以用奇异值分解(SVD)方法.(2)直接从误差表达式(6)式出发,计算M个远震事件的预测误差平方和,E=Mi=1

18、Ei=Mi=1t(eit)2.(11)对(11)式求极小,即可得到待求的接收函数.本文所采用的是第二种方法.4资料检验图1合成地震图及其接收函数(a)径向分量;(b)垂直分量;(c)接收函数.Fig.1The synthetic seism ogram and receiver function(a)Radial com ponent;(b)Vertical com ponent;(c)Receiverfunction.我们分别用合成地震图和实际观测到的远震资料对本文所提出的多道滤波反褶积方法进行了数值检验,以考察多道滤波反褶积的实际应用效果.首先,我们用合成地震图检验多道滤波反褶积方法.我们

19、采用的是单层半无限空间介质,但在地壳3973期吴庆举等:用多道反褶积方法测定台站接收函数名师资料总结-精品资料欢迎下载-名师精心整理-第 3 页，共 6 页 -?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/和上地幔中分别加了一个低速薄层和高速薄层,薄层厚度均为2 km,速度变化量均为012 km s,以模拟小速度差界面产生的弱转换波震相.图 1 是脉冲响应的合成地震图及其理论接收函数.图3合成地震图的比较(a)实线:合成接收函数,虚线:单道反褶积接收函数;(b)实线

20、:合成接收函数,虚线:多道反褶积接收函数.Fig.3C omparison for synthetic receiver function(a)Solid:Synthetic receiverfunction,Dot:the stackedreceiverfunction from singleevent deconv olution;(b)S olid:Synthetic receiverfunction,Dot:the receiverfunction from multi2channeldeconv olution.图2远震波形数据的垂直分量Fig.2The vertical com p

21、onent of teleseimsic waveform为了在合成地震图中引入震源时间函数和随机噪音,更为逼真地模拟实际的观测地震图,我们选取了 10 个不同远震事件的P 波波形的垂直分量作为等效震源时间函数(图 2),随机噪音自然也包括其中.分别用图2 所示的等效震源时间函数与图1 所示的合成地震图的径向和垂直分量做褶积运算,得到了 10 个虚拟的远震波形记录.我们首先用单道Wiener 滤波反褶积方法计算出相互独立的10 个接收函数,然后叠加平均,得到10 个虚拟远震事件接收函数的均值.我们再用本文所论述的方法对上述10 个虚拟远震事件进行多道反褶积处理,得到一个统一的接收函数.两种方法

22、所得结果示于图3.对比图 3a 和图 3b 不难看出,多道滤波反褶积497地 球 物 理 学 报(ChineseJ.Geophys.)50卷 名师资料总结-精品资料欢迎下载-名师精心整理-第 4 页，共 6 页 -?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/能得到与理论值近乎完全一致的接收函数,其误差要远远小于单道反褶积再叠加的结果,表明多道滤波反褶积的有效性.需要强调的是,单道反褶积再叠加法造成了地壳和上地幔弱转换波震相(分别位于Moho界面强转换波前后,介于

23、0715 s 之间)的明显畸变,而多道滤波反褶积则能非常精确地恢复弱转换波震相.我们又进一步用实际地震数据对单道反褶积方法和多道反褶积方法作了对比.我们选取了同一台站记录到的来自同一方向、震中距变化不大的10 个远震事件的三分量P波波形数据,分别用单道反褶积或多道反褶积提取了接收函数,结果示于图4.从图 4 可以看出,与单道反褶积相比较,多道滤波反褶积所提取的接收函数的波形更加尖锐,振幅更强.值得指出的是,多道滤波反褶积所提取的接收函数在35s以后明显强于单道反褶积,而上地幔过渡带之间的弱转换波震相一般出现在3580 s这一时窗范围内.图4实测地震图的对比虚线:单道反褶积接收函数,实线:多道反

24、褶积接收函数.Fig.4Comparison of real dataDot:the stackingreceiverfunction from single event deconv olution;S olid:the receiverfunction from multi2channel deconv olution.5 结 论多道滤波反褶积方法依据一定的准则,以来自同一方向、震中距变化不大的多个远震事件的垂直分量作为输入,水平分量作为输出,设计多道滤波因子,从而从多个远震事件中测定出统一的台站接收函数.在接收函数的计算过程中,若干个看似不相干的远震波形数据,在多道滤波算子的作用下,根据

25、波形的相干性,消除震源时间函数及路径效应,压制噪音,保留并增强相干信号.合成资料与实际资料的检验表明,多道滤波反褶积不仅能够有效地恢复接收函数,而且具有一定的、从背景噪音中提取弱间断面转换波的本领,提高接收函数的分辨率,在提取台站接收函数方面具有潜在的优势,值得加以应用和推广.参考文献(References)1 OwensT J,Zandt G,TaylorS R.Seism ic evidencefor an ancientriftbeneaththe CumberlandPlateau,Tennessee:A detailed analysisofbroadbandteleseism ic

26、 P waveforms.J.Geophys.Res.,1987,89(B9):77837795 2 刘启元,Rainer K,李顺成.接收函数复谱比的最大或然性估计及非线性反演.地球物理学报,1996,39(4):502511Liu Q Y,Rainer K,LiS C.Maximal likelihood estimation andnonlinearinversionof the com plex receiverfunction spectrumratio.ChineseJ.Geophys(in Chinese),1996,39(4):502511 3 刘启元,Rainer K,李顺成

27、.中国数字地震台网的接收函数及5973期吴庆举等:用多道反褶积方法测定台站接收函数名师资料总结-精品资料欢迎下载-名师精心整理-第 5 页，共 6 页 -?1994-2009 China Academic Journal Electronic Publishing House.All rights reserved.http:/其非线性反演.地球物理学报,1997,40(3):356368Liu Q Y,Rainer K,Li S C.The receiverfunctionsat the stationsofthe Chinese digital seism icnetw ork(CDS N

28、)and their nonlinearinversion.ChineseJ.Geophys.(in Chinese),1997,40(3):356368 4 刘启元,李顺成,沈 扬等.延怀盆地及其邻区地壳上地幔速度结构的宽频带地震台阵研究.地球物理学报,1997,40(6):763773Liu Q Y,Li S C,Shen Y,et al.Broadbandseism ic array studyofthe crust and upper mantle velocity structure beneath Yanhuai basinand its neighbouring region.C

29、hinese J.Geophys.(in Chinese),1997,40(6):763773 5 刘启元,陈九辉,李顺成等.新疆伽师强震群区三维地壳上地幔S波速度结构及其地震成因的探讨.地球物理学报,2000,43(5):356365Liu Q Y,ChenJ H,Li S C,et al.Passiveseism ic experimentinXingjiang2Jiashistrong earthquake region and discussion on itsseism ic genesis.Chinese J.Geophys.(in Chinese),2000,43(5):3563

30、65 6 刘启元,Rainer Kind,陈九辉等.大别造山带壳幔界面的断错结构和壳内低速体.中国科学(D),2005,35(4):304313Liu Q Y,Rainer K,ChenJ H,et al.Dislocation structureof thecrust2mantle boundaryand low2velocity bodywithin the crust beneaththe Dabie Shancollisionorogen.Science in China(SeriesD),2005,48(7):875885 7 吴庆举,曾融生.用宽频带远震接收函数研究青藏高原的地壳结

31、构.地球物理学报,1998,41(5):669679Wu Q J,Zeng R S.Thecrustal structureof Qinghai2Xizang plateauinferred from broadbandteleseism ic waveform.ChineseJ.Geophys.(in Chinese),1998,41(5):669679 8 ZhuL P,Owens T J,Randall G E.Lateral variation in crustalstructureof the Northern Tibetan Plateau inferred from telese

32、ism icreceiver functions.Bull.Seismol.Soc.Am.,1995,85(6):15311540 9 Levin V,ParkJ.Crustalanisotropy beneaththe Ural mtns foredeepfrom teleseism ic receiver functions.Geophys.Res.Lett.,1997,24:1283128610 ZhuL P,Kanam ori H.Moho depth variation in s outhern Californiafrom teleseism ic receiverfunction

33、s.J.Geophys.Res.,2000,105:2969298011 ZhuL P.Crustal structure acrossthe San Andreas Fault,S outhernCalifornia from teleseism ic converted waves.EarthandPlanetarySci.Lett.,2000,179:18319012 吴庆举,曾融生,赵文津.喜马拉雅 青藏高原的上地幔倾斜构造与 陆-陆碰 撞过 程.中 国科 学(D),2004,34(10):919925Wu Q J,Zeng R S,Zhao W J.The upper mantle

34、structure of theTibetan Plateau and itsimplication forthe continent2continentcollision.Science in China(SeriesD),2005,48(8):1158116413 Ai Yinshuang,ZhengTianyu,Xu W W,et al.A com plex 660 kmdiscontinuity beneath northeast China.EarthandPlanet.Sci.Lett.,2003,212:637114 Ai Yinshuang,Zheng Tianyu.Theup

35、permantle discontinuitystructure beneath easternChina.Geophys.Res.Lett.,2003,30(21):2089,doi:10110292003G L01767815 Ai Y S,Zhao D P,Gao X,et al.The crust and upper mantlediscontinuitystructurebeneath Alaskainferredfromreceiverfunctions.Phys.Earth Planet.Inter.,2005,150:33935016 Tian X B,Wu Q J,Zhang

36、Z J,et al.Joint imaging by teleseism icconvertedand multiple wavesandits application in the INDEPTH2passiveseism icarray.Geophys.Res.Lett.,2005,32:L21315101102917 G urrola H,Baker G E,Minster J B.S imultaneous time2domaindeconv olution withapplication tothecom putationofreceiverfunctions.Geophys.J.I

37、nt.,1995,120:53754318 吴庆举,田小波,张乃铃等.用Wiener滤波方法提取台站接收函数.中国地震,2003,19(1):4147Wu Q J,Tian X B,ZhangN L,et al.Receiverfunction estimatedby Wiener filtering.Earthquake Research inChina(in Chinese),2003,19(1):414719 吴庆举,田小波,张乃铃等.计算台站接收函数的最大熵谱反褶积方法.地震学报,2003,25(4):382389Wu Q J,Tian X B,ZhangN L,et al.Recei

38、verfunction estimatedby Wienerfiltering.Acta Seismologica Sinica,2003,19(1):414720 ParkJ,Levin V.Receiver functions from multiple2taper spectralcorrelation estimates.Bull.Seismol.Soc.Am.,2000,103(26):89991721 Ligorria J P,Amm on C J.Iterative deconv olution and receiver2function estimation.Bull.Seismol.Soc.Am.,1999,89(1):395400(本文编辑 胡素芳)697地 球 物 理 学 报(ChineseJ.Geophys.)50卷 名师资料总结-精品资料欢迎下载-名师精心整理-第 6 页，共 6 页 -