接收函数波动方程叠后偏移方法

1、接收函数波动方程叠后偏移方法介绍陈凌，张耀阳中国科学院地质与地球物理研究所接收函数波动方程叠后偏移方法是我们借鉴勘探地震学中发展成熟的反射波偏移成像技术，2005年以来发展的一种新的天然地震接收函数偏移成像方法（Chenetal.,2005a;2005b）。该方法包括两个基本步骤：时间域CCP叠加和深度域波场反向延拓。与反射地震学中处理反射波资料的CMP叠加类似，CCP叠加用以提高资料的信噪比。波场反向延拓则是一个将产生Ps转换波（图1d）的转换波源偏移至其真实位置的过程。在CCP叠加中，用一个1D参考模型将来自单个地震事件的接收函数按垂直入射情况做时间校正后，根据其转换点位置进行组合叠加。所

2、得到的CCP叠加道集可以近似看作是零偏移距资料，即转换波从速度间断面的转换点处垂直向上传播至地表而被记录到的波场。对于复杂结构，只采用CCP叠加将会给成像结果带来假像，而偏移对于转换界面的正确成像是必不可少的（RybergandWeber,2000;PoppeliersandPavlis,2003a）。我们采用的波场偏移原理与反射地震学中普遍采用的爆炸反射面原理（Claerbout,1985;SheriffandGeldart,1995）类似。在爆炸反射面模型中：每一个反射界面都被当作是布满了密度正比于界面反射强度的爆炸源。将这些爆炸源同时在时间t=0激发，地震波以实际速度的二分之一向上传播以

3、满足双向传播走时（图1a,b）。这样在地表上接收的波场可以模拟由CMP叠加得到的零偏移距（零炮点-接收点距离）道集。倒转上述过程，即将地表观测的零偏移距波场反向传播至t=0时刻，就可以得到一个使所有反射体正确归位的深度偏移成像结果。在接收函数偏移中，通过CCP叠加产生的零偏移距接收函数道集可以用来模拟地表记录到的将所有地下转换波源同时在t=0时刻激发得到的波场，而源的强度正比于Ps的转换系数（图1c）。这种模型称为爆炸转换波源模型，是由Pavlis（2003）首先提出并应用于接收函数成像的。将接收函数的CCP叠加道集反向延拓至整个空间并在t=0时刻成像，就能使所有转换波偏移至产生转换波的相应地

4、下速度间断面的真实位置处。Ryberg和Weber（2000）提出了同样的接收函数偏移方案。我们所提出的偏移方法与Ryberg和Weber方法的本质区别在于采用的波场延拓技术不同。在他们的研究中采用基于射线的Kirchhoff高频近似方法，而我们则应用波动方程波场传播技术。1.CCP叠力口传统的接收函数CCP叠加过程包括接收函数的走时校正和组合叠加(DuekerandSheehan,1997,Zhu,2000;Aietal.,2003)。这种数据处理过程与反射地震学中的CMP叠加类似，但两者之间仍存在着显著差别。首先，在接收函数处理过程中，时间校正是为了补偿由不同震中距入射P波造成的P波转换成

5、S波的走时差。因此，通常是以一个1D参考模型计算得到的参考震中距的走时作为标准，对所有接收函数进行时间校正。参考震中距一般选取为所有资料的平均震中距，以使时间校正量达到最小。由于接收函数的典型震中距是30。至U90。，参考震中距也应落在该范围内，因此必然对应于P波倾斜入射的情况(入射平面波的水平慢度p丰0)。另一方面，在反射资料分析中，反射波的走时均按垂直入射(p=0)进行校正。如此得到的CMP叠加道集才能近似作为叠后偏移所需要的零偏移距资料。为了应用叠后偏移对接收函数进行处理，定义相对于假设零偏移距情况(震中距为180,p=0)的Ps转换波走时校正量为source(d)Psconversio

6、nandmultiplereflectionsTps(P)=Tps(p)-Tps(0)其中接收函数中PS转换波走时Tps是基于展平地球模型中平面波入射假定得到的直达P波与转换S波的走时差(Gurrolaetal.,1994;DuekerandSheehan,1997)0Tps(P)=(、Vs1z)-p2一,VP2(z)-p2)dz(2)D上式中D,Vp,Vs分别为Ps转换深度，随深度变化的P波速度和S波速度。在走时校正之后，根据各个深度处的转换点位置将所有接收函数重新组合到各自的共转换点单元中，以达到最佳聚焦的目的。对共转换点单元中的接收函数分别进行叠加而构成的CCP道集近似相当于P波垂直入射

7、的观测资料，因而将作为下一步叠后深度偏移的初始波场。2.基于波动方程的波场反向延拓为了由接收函数的CCP叠加道集对产生Ps转换波的地下间断面真实成像，需要采用一个适当的波场延拓方法，从地表的时间域波场(观测)出发推导出空间域t=0时刻的波场(像)。目前已经发展的接收函数偏移方法大多以射线近似作为波场延拓的理论基础(RybergandWeber,2000;Sheehanetal.,2000;PoppeliersandPavlis,2003a,2003b;Niuetal.,2005)。我们则采用波动方程的单向波算子分解，并基于扰动理论，在频率-波数域进行波场反向延拓。为了便于表达，我们只考虑2D模

8、型。下面的公式推导可以非常直接地推广到3D情况。频率-空间域的2D标量波动方程可表示为rc2c2e2以；u(x,zf)=0(3)_cxczv(x,z)其中圆频率，V(x,z)为介质波速，u(x,zp)表示频率-空间域的波场。对上式采用不同的波场分解和近似技术将得到不同类型的波场传播算子和延拓方法。这里我们采用相位屏算子进行接收函数的波场延拓。相位屏算子根据扰动理论，介质的速度场可以分解为背景速度场V0和相应的扰动场,V(x,z)v(x,z)=Vo(z)v(x,z)将公式带入公式(3),得到一2一2ccI+.2.2xz其中k0/v0(z)为背景波数，b(x,z)=v0(z)/v(x,z)为扰动函

9、数。分解公式(5)将产生两个单向波波动方程，一个表示前向传播，另一个表示后向传播。控制前向传播的方程为u(x,z,).2.2F2/、=ikob2u(x,zj)(6).:z：x基于不同的考虑，上式右端的平方根算子可以展开为扰动项(DeHoopetal.,2000)或Pad赤歹U(XieandWu,1998;Jinetal.,2002)。波动方程相应地表示为cu(x,z,)铲=ikb|1+Qk0,b,-u(x,z,6)也;；打)一上式中函数Q弋表平方根算子的高阶展开。另一方面，波场可以分解为不同波数的平面波u(x,z,)=dkxUkx,z,eikxx(8)其中U(kx,z,co)eikxx表示振幅

10、为U(kx,z,)的平面波，kx为相应的横向波数。对公式沿x方向进行傅里叶变换，得到-U(kv,z,)iI：；vx=ikz-k0(FTx)Y(kxk,b)U(kx,z,)(9)z_v其中kz=Jk；-k2为背景垂直波数。FTxU表示沿x方向的傅里叶变换，表示波数域卷积，而Y(kx,k0,b)为kQ沿x方向的傅里叶变换。在速度小扰动(6v(x,z)v(x,z)和小角度近似(/心)条件下，可以忽略上式中的Y。进一步考虑屏近似，即在波场从深度乙传播到zi+1的薄板中，如果薄板足够薄,则可以忽略扰动函数b(x,z)在z方向上的变化。于是波场简化为(Stoffaetal.,1990;XieandWu,1

11、998)U(kx,z+,)=exp(ikz(zi+-z)FTxexp-i, v(x,zi)v(x,z)%(Z+-zj u(x,4 (10)再通过一次傅里叶反变换，就得到了空间域的波场u（x,z中Q）。在横向扰动较大情况下，对公式（9）中的函数Y给出不同近似表示就构成了不同类型的广义屏传播算子（XieandWu,1998;DeHoopetal.,2000）。公式（10）表明，每一深度处的相位屏波场延拓过程包含了两步：计及横向速度扰动的空间域相位校正exp产叱1z）附基于相移算子expQkzd-z）的波数域背景介质v（x,z|）I中的自由传播。通过迭代执行这两个步骤，波场就由地表延拓到整个空间。3

12、 .偏移速度上述波场延拓的实现需要一个速度模型。如公式（2）所示，接U函数中Ps转换波CCP00 dtpsM0)M逐-.1(11)的走时与反射资料不同，是直达P波与转换S波的走时差。在时间校正后的叠加道集中，该走时差变为其中V(z)Vp(z) Vs(z)Vsz)-Vpz)Vp(z)-Vs(z)(12)为折合波速。将上式拓展到速度存在横向变化的情况，我们可以定义任意介质中的折合速度为：Vp(x,z)Vs(x,z)V(x,z)=(13)Vp(x,z)-Vs(x,z)V（x,z）即为波场反向延拓公式（10）中的偏移速度。4 .成像条件在叠后偏移过程中，采用上述波动方程传播算子将地表观测波场逐个频率地

13、反向延拓至整个空间。根据爆炸转换波源模型，将偏移后的波场对频率进行叠加得到的1=0时刻的波场即为最终的像I(x,z)=u(x,z,t=0)=u(x,z,)(14)5 .讨论在上述波动方程叠后深度偏移方法中，虽然接收函数的CCP叠加隐含了水平层状结构的假设，但横向速度变化可以包含在偏移速度模型中。因此，叠后偏移过程能够计及结构横向不均匀性所带来的传播效应。相位屏传播算子已被证明在处理速度横向变化平缓情况时具有足够的精度，甚至当速度反差达到40%寸仍能较准确地对结构成像(Stoffaetal.,1990)。在上面介绍的偏移框架下，其它波动方程传播算子，如基于波场小波分解的小波束域传播算子(陈凌，2

14、002;陈凌等，2004;Chenetal.,2006a)、各种广义屏算子等(XieandWu,1998;DeHoopetal.,2000)可以很便捷地插入到偏移程序中替代这里所采用的相位屏传播算子。因此，可以根据所要研究的结构对象选取最适合的波场传播算子，以达到最佳的偏移成像精度和效率。止匕外，当结构特别复杂，叠后偏移无法满足成像精度时，需要采用叠前偏移的成像方法；而当2D假设不再成立时，就必须采用3D波场延拓。基于上述偏移方法，叠前偏移和3D偏移都可以参照反射地震学中的相应过程来实现。另外，上面介绍的偏移过程虽然是对Ps转换波成像，但其原理和方法同样适用于其它类型的地震震相，如SpW换波，

15、Pp和Ss震相，以及各种地表多次反射波。我们已经将Ps转换波波动方程偏移方法进行了拓展，使其同样适用于地表多次波PpPs(Chenetal.,2006c)和SpW换波震相的偏移(Chenetal.,2008)。我们已分别或联合应用上述不同地震波震相偏移技术，开展了华北、东北及日本俯冲带地区壳幔间断面结构的成像研究(如Chen,2009;ChenandAi,2009;Chenetal.,2008;2009;2014;王炳瑜等,2013)。值得提出的是，由于波场延拓是在频率域完成的，因此，不论是采用相位屏算子、小波束算子，还是广义屏算子等，都可以十分方便地获得不同频率成分的偏移图像，即直接对相应频

16、率的波场进行叠加。而在基于叠加或时间域的接收函数成像方法中，当考虑资料的不同频率成分时，需要对资料进行滤波处理。与这些方法相比，频率域的波动方程偏移技术自然成为研究频率(或尺度)相关结构特性的有效工具(如Chenetal.,2005b;2006b;2014;Chen,2009)。参考文献博士论文.北陈凌.小波束域波场的分解、传播及其在地震偏移成像中的应用京:中国地震局地球物理研究所,2002。陈凌，吴如山，王伟君，2004，基于Gabor-Daubechies小波束叠前深度偏移的角度域共成像道集，地球物理学报，47(5)，876-885。王炳瑜，陈凌，艾印双，何玉梅，华北克拉通东北部及邻区地壳

17、和地幔转换带厚度研究，地球物理学报，56(11)，60-68，2013。Ai,Y.,T.Zheng,W.Xu,Y.He,andD.Dong,Acomplex660kmdiscontinuitybeneathnortheastChina,EarthPlanet.Sci.Lett.,212,63-71,2003.Claerbout,J.F.,ImagingtheEarthsInterior,BlackwellSci.,Malden,Mass.,1985.Chen,L.,L.X.Wen,T.Zheng,AWaveEquationMigrationMethodforReceiverFunctionIm

18、aging:1.Theory,J.Geophys.Res.,110,B11309,doi:10.1029/2005JB003665，2005a.Chen,L.,L.X.Wen,T.Zheng,AWaveEquationMigrationMethodforReceiverFunctionImaging:2.ApplicationtotheJapansubductionzone,J.Geophys.Res.,110,B11310,doi:10.1029/2005JB003666,2005b.Chen,L.,R.S.Wu,Y.Chen,Target-orientedbeamletmigrationb

19、asedonGabor-Daubechiesframedecomposition,Geophysics,71(2),S37-S52,2006a.Chen,L.,T.Zheng,W.Xu,AThinnedLithosphericImageoftheTanluFaultZone,EasternChina:ConstructedfromWaveEquationBasedReceiverFunctionMigration,J.Geophys.Res.,111,B09312,doi:10.1029/2005JB003974,2006b.Chen,L.,T.Zheng,W.Xu,Receiverfunct

20、ionmigrationimageofthedeepstructureintheBohaiBayBasin,easternChina,Geophys.Res.Lett.,33,L20307,doi:10.1029/2006GL027593,2006c.Chen,L.,T.Wang,L.Zhao,T.Zheng,DistinctLateralVariationofLithosphericthicknessintheNortheasternNorthChinaCraton,EarthPlanet.Sci.Lett.,267,56-68,2008.Chen,L.,C.Cheng,Z.G.Wei,Se

21、ismicevidenceforsignificantlateralvariationsinlithosphericthicknessbeneaththecentralandwesternNorthChinaCraton,EarthPlanet.Sci.Lett.,286,171-183,2009.Chen,L.,Y.S.Ai,DiscontinuityStructureoftheMantleTransitionZonebeneaththeNorthChinaCratonfromReceiverFunctionMigration,J.Geophys.Res.,114,B06307,doi:10

22、.1029/2008JB006221,2009.Chen,L.,LithosphericstructurevariationsbetweentheeasternandcentralNorthChinaCratonfromS-andP-receiverfunctionmigration,Phys.Earth.Planet.Inter.,173,216-227,2009.Chen,L.,M.M.Jiang,J.H.Yang,Z.G.Wei,C.Z.Liu,Y.Ling,PresenceofanintralithosphericdiscontinuityinthecentralandwesternN

23、orthChinaCraton:Implicationsfordestructionofthecraton,Geology,42(3),223-226,2014.DeHoop,M.,J.Rousseau,andR.S.Wu,Generalizationofthephase-screenapproximationforthescatteringofacousticwaves,WaveMotion,31,43-70,2000.Dueker,K.G.andA.F.Sheehan,MantlediscontinuitystructurefrommidpointstacksofconvertedPtoS

24、wavesacrosstheYellowstonehotspottrack,J.Geophys.Res.,102,8313-8327,1997.Gurrola,H.,J.B.Minster,andT.Owens,Theuseofvelocityspectrumforstackingreceiverfunctionsandimaginguppermantlediscontinuities,Geophys.J.Int.,117,427-440,1994.Jin,S.,C.C.Mosher,andR.S.Wu,Offset-domainpseudoscreenprestackdepthmigrati

25、on,Geophysics,67,1895-1902,2002.Niu,F.,A.Levander,S.Ham,M.Obayashi,MappingthesubductingPacificslabbeneathsouthwestJapanwithHi-netreceiverfunctions,EarthPlanet.Sci.Lett.,239,9-17,2005.Pavlis,G.L.,ImagingtheEarthwithpassiveseismicarrays,LeadingEdge,22,224-331,2003.Poppeliers,C.andG.L.Pavlis,Three-dimensional,prestack,planewavemigrationofteleseismicP-to-Sconvertedphases:1.Theory,J.Geophys.Res.,108(B2),2112,doi:10.1029/2001JB000216,2003a.Poppeliers,C.andG.L.Pavlis,Three-dimensional,prestack,planewavemigrationofteleseismicP-to-Sconvertedphases:2.Stackingmultipleevents,J.Geophy