摘要：随着地震学的发展以及观测资料的积累，从20世纪90年代开始，人们开始直接利用地震波形成像的方法来研究地球内部界面，主要是借鉴工业界发展的一些成熟的主动源方法。最开始前人将工业界CMP（共中心点）叠加方法引入到接收函数中，后来朱露培教授及其合作者在2000年提出了CCP（共转换点）叠加的接收函数方法并得到国际同行的广泛使用。报告主要介绍了他们将工业界使用的逆时偏移成像方法引入到远震接收函数的研究进展，阐述了将逆时偏移成像引入到远震接收函数中的基本算法思想，通过一系列合成数据的理论测试验证了该方法的正确性，最后将这一新方法应用到一条长300 km台间距为5-10 km的测线数据中。通过与传统CCP成像结果进行对比，报告指出针对目前的天然地震观测数据情况，CCP成像结果已经差不多是最好的了，但是将来随着更为密集的台阵观测资料的积累，这种逆时偏移的远震接收函数方法可以发挥更大的优势，建议开展天然地震观测研究时，尽量布设密集的(台间距小于5 km)观测台阵，以获得更为精细的成像结果。

　　我今天的报告主要给大家介绍一下我们近期将相移和插值（PSPI, Phase Shift Plus Interpolation）方法引入到远震接收函数的研究。首先我会对PSPI方法进行一下介绍，然后为了大家能够看到这个方法的有效性，我们设计了一系列的合成数据进行理论测试，接着将这种方法与我之前发展的共转换点叠加（CCP, Common Conversion Point）的接收函数方法进行一个比较，看一下它们两个各自的优缺点，最后我把这个方法应用到了我在美国的一个线性密集台阵数据中，得到了相应的成像结果并进行讨论。 　　  

　　一、科学问题　　  

　　首先是科学问题及其重要性。利用接收函数对地球内部的界面进行成像是地震学研究的一个重要部分，我们知道地球内部包含有多个间断面，例如莫霍面、410 km、660 km界面等，这些内部界面蕴含了很多关于地球构造的信息，它们的深度位置、横向变化能够帮助我们在认识地球内部，例如物质成分、温度、相变以及一些动力学过程等提供相关的依据。以前对界面的成像主要是通过主动源（active source）的方法，工业界用得较多的就是近垂直反射，以及天然地震和人工测深研究用到的宽角反射/折射，这些早期研究地球界面的传统方法主要是利用各种震相的到时来确定界面的位置。 

　　随着地震学的发展以及观测资料的积累，大概是从20世纪90年代开始，人们开始直接利用地震波形成像的方法来研究地球内部界面，主要是借鉴工业界发展的一些成熟的主动源方法进行成像。现在天然地震中对界面成像用得最多的可能就是接收函数方法了。这里以远震接收函数为例，一个P波从下面传播到界面上，一部分能量转换成S波，然后可以通过记录到的P波和S波的相对到时来确定这个界面的位置。相对于主动源方法，接收函数主要有两个明显的优势：一是在主动源方法中，地表的震源激发的地震波先向下传播再反射回来，两次穿过介质，如果浅部结构比较复杂的话，接收到的信号噪音就比较大，而接收函数的地震波在地球内部传播，到达地表台站只需要一次穿过浅部介质，因此噪音相对较低；同时，地表主动源的穿透深度有限，因此能够探测到的界面也相对较浅，而接收函数方法因为使用的是天然地震源，地震波从下面传播上来，不存在穿透深度的问题。 　　  

　　二、历史背景 

　　最早利用接收函数做成像工作的，我认为是Dueker和Sheehan^[1]。他们在1997年的文章将工业界共中心点叠加方法（CMP, Common Mid-Point）引入到接收函数中，其思想很简单，先通过一个背景速度模型将倾斜入射的接收函数动校正到垂直入射，然后将一个小范围内的接收函数进行叠加以提高信噪比，接着再用一个速度模型将时间域的信号转换到深度域，在不同的点做完这样的处理后就可以拼成一个最终的图像。借鉴CMP方法这样处理有几个问题，首先是做动校正时需要一个区域速度模型，在计算中会带入误差，而且它也没有进行偏移，直接从时间域转换到了深度域。 

　　我1998年博士毕业后到南加州做博后时，他们在南加州布设了一条一百多公里的短周期测线，大概有一个月的远震事件记录，我拿到数据后想确定莫霍面的形状，但是短周期数据噪音很大，当时我就想，与其利用CMP这种叠加方法还不如直接进行叠前偏移，也就是利用一个背景模型直接将每个接收函数转换到空间域，转换完成后再在空间域进行叠加，我把这个方法称为CCP(Common Conversion Point)^[2]。 

　　CCP与CMP一个明显的区别就是不需要再进行时间域到空间域的转换，因为第一步做射线计算的时候就已经考虑了，它最大的好处是对每个深度都是聚焦的，而CMP方法只对某一个深度聚焦。CCP方法也需要一个1D的模型，我现在的code（代码）可以在每个台都用一个不同的1D模型，相当于一个准三维的模型。以上是对这两种方法的比较，我之所以在这个地方花点时间介绍，主要是我发现有人发表的文章中实际上用的是CMP，但把它叫做了CCP，实际上这是一个错误。我们当时的文章发表后，发现欧洲的学者Kosarev^[3]在1999年提出了一个类似的方法，但他们没有把其方法叫做CCP，CCP是当时我为了和CMP区别而取的名字。因为使用1D模型很难处理复杂介质的情况，因此后来又有人将工业界中使用的三维偏移成像方法引入到了接收函数中，例如广义射线方法（GRT, Generalized Ray Theory），它的缺点是不能处理介质变化太剧烈的情况。另外一类就是贵所陈凌提出的基于波动方程的方法^[4]，这种方法对复杂介质也能进行较好的处理。 

　　我在美国有时会去石油公司做consulting（技术顾问）,可以学到很多工业界的方法。我在给一个公司做consulting时就学到了逆时偏移这个方法。逆时偏移的原理其实很简单也很巧妙，其涉及到两个波场，一个是震源的波场，一个是接收点的波场，首先把震源的波场按正的时间传播到各个深度上，可以把每个深度上的波场求出来（可以是时间域也可以是频率域），接下来把记录到的接收点的波场进行时间反转，相当于逆时地从地表传回到空间中每个点上，然后把两个波场做相关，得到最后的图像。但是这种方法的计算量非常大，波场正传和逆传都要用有限差分计算，我当时去做consulting时他们就是想提高计算效率，有一种近似的方法可以做到，就是我前面提到的PSPI，这种方法也很巧妙，同样也是工业界发展起来的。 

　　在一维模型的情况下，不同深度的波场在波数域里有一个相位差，其中一个波场乘以这个相位差就可以得到第二个波场，但如果不是一维均匀模型，也就是在横向上有速度变化的话，就不能直接这样处理了，为此工业界上发展了一种叫做插值的方法，因为你不知道用哪个速度好，所以就用一系列的速度把它从一个深度传到另一个深度，这个速度不一定是真实的速度，做完以后再把它变到空间域，在每个点找到它对应的速度，然后用这个速度值找到最近的两个参考速度，然后插值得到最后的速度。工业界上他们并没有去证明这个方法，他们觉得如果能用就是好的。这是我当时给他们做的一个例子（图1），这是一个很复杂的速度模型，地表有一系列的炮点和接收点，然后用PSPI方法去做成像，我们可以看到成像结果里，这些界面的位置都能跟已知模型很好地对上，由此可以看出这种插值方法可以有效解决三维复杂介质中波传播的问题。 

图1  对一个复杂的速度模型（上）进行PSPI逆时偏移成像（下）　  

　　三、接收函数的逆时偏移　　  

　　PSPI方法的优势概括来说，就是提升了逆时偏移的计算效率，接下来我们就把PSPI用到接收函数上面。在接收函数中唯一要改变的就是不在是地表的震源了，而是从地底下传播上来的平面波，我们把平面波从下往上传上去，然后接收函数的波场逆时地从上往下传下来。这个过程有个问题，就是做完反褶积算的接收函数已经没有了走时的信息，它的P波都是在零时刻，所以直接用这个波场去传肯定不满足声波方程了，我当时跟学生讨论以后，发现这个问题其实也好解决，我们把接收函数的每个点都做个时移，把它移到相应的到时，这样操作之后这两个波场就全部match（匹配）了。接收函数里面除了界面上的一次转换波外，还有多次转换波，例如P波可以传到地表，反射回去再在界面上发生转换，所以用这种方法也可以对多次波进行成像。　　  

　　四、合成数据的理论测试　　  

　　下面就是来测试这个方法的正确性了。我们设计了一系列的理论模型，第一个是呈正弦形态变化的莫霍面（横向变化周期为50 km），地壳一个速度，地幔一个速度，在地表设置台间距为1 km 的100个接收点，接下来利用这个模型计算远震P波接收函数，从图像上我们可以看到PS转换波振幅有变化，这是因为振幅取决于平面波的入射方向和界面形成的夹角大小，近垂直入射的话振幅是很小的，另外仔细看的话，我们能观察到S波造成的的衍射，射线理论当然不会计算出这样的东西，但是用波场算的话，因为有限频的原因，所以会出现衍射的信号。 

　　同样的数据我们用CCP也做了一个成像结果，可以看出这两种方法的结果差别，CCP做出来的界面形态基本上也是对的，差别就是CCP做出来的深度比实际的要浅一点，这是因为CCP用的是每个台下面的一维速度模型，这个速度模型并不能反映射线经过的真实速度结构，另外CCP也不能解决这种衍射的情况，即使在比较光滑的速度变化情况下，可以看出逆时偏移方法要比CCP结果好。接下来我把莫霍面的形态变化变得更剧烈一点（横向变化周期降为20km），然后采用逆时偏移和CCP分别得到各自的成像结果，这个时候可以看到CCP的结果把界面就搞错了（界面变化剧烈的时候它成像的位置都是错的），而逆时偏的成像结果还是比较干净的。 

　　第三个测试模型是莫霍面出现陡变的情况，这种构造在真实地球介质中也是可能存在的，算出来的接收函数在这种陡变情况下有很多衍射的信号，CCP的结果在远离陡变带的界面结果都是对的，但在陡变区域附近搞错了。我曾见过有人写文章将衍射信号认为是俯冲的信息，这是没考虑可能出现的衍射情况，其实工业界对这个早就注意了，就是成像中所谓的ghost现象。前面都是用的二维模型，接下来用了一个简单的含有一个倾斜界面的三维模型，射线可以从各个方向传过来，前面二维模型中射线只从左边或者右边来，这个模型是用来测试对三维模型的有效性，从结果可以看到成像出来的界面都是对的，说明我们的code是没有问题的。　　  

　　五、实际数据处理　　  

　　最后我们把这个方法用到了实际数据的处理上，美国USArray项目观测2014年到2016年经过Wabash Valley时我们同时布设了一条测线，前人在该区域做的深部构造工作很少，我们主要想看看这个区域的地震带和它南边的地震带在结构上有些什么样的差异。测线长300 km，内侧台间距5 km，外侧是10 km，一共用了40多套宽频带地震仪，观测了两年多。逆时偏方法对数据要求很高，对一个地震事件，它要求最好这条线上所有台站都能记录到，台站最好还是等间隔的。我们一共找到了40多个地震，台站分布并不是等间距的，这里我把每个接收函数画到5 km宽，有些地方有空缺，没办法所以只好插值。我有个去年毕业的博士生，发展了一种叫做非线性伸缩插值的方法，给这样一个波场，他可以把它插成等间距的。背景模型是一个二维的，来自于背景噪音和接收函数联合反演得到的S波速度结果，因为要构造一个三维的模型，因此我们假定在垂直方向上速度结构没有变化。这是最后的成像结果（图2），上面是用逆时偏移得到的成像结果，下面是CCP的成像结果。　　 

图2  逆时偏移远震接收函数成像（上）和远震接收函数CCP叠加成像（下）对比

　　如果我要说上面逆时偏移得到的结果更好，大家可能不一定同意，但是首先我认为逆时偏移成像结果的细节要比CCP多，CCP因为要做水平光滑，很多细节都被抹掉了。总体来看，就我们目前这个数据水平，实际上我觉得如果构造变化不是很剧烈的话，CCP应该是最好的了，尽管这是二十年前我用的了，但是我们是为了将来，将来我们有了密集台阵，可以发挥这个方法的优势。 

　　关于台间距对成像结果的影响，我也做了一些测试和推导。刚才是1 km台间距做理论测试，现在改成2 km台间距，Moho变化还是50 km，2 km台间距的也能成像出来，没有问题。对复杂的构造，2 km台间距也大体能成出来。到了8 km台间距基本就不行了，变化缓的还能大体上看出来形态，变化快的基本上没有了。我在文章里面讨论了这个问题，成像结果是与你的空间采样有关的，我做了公式推算，我的推算结果就是在1Hz情况下，5 km以下台间距的成像结果都还可以的，到了8 km结果就已经不好了，所以我的结论就是布台要在5 km以内，10 km以上的只能做CCP，要用波场的话就不行了。 

　　最后是关于多次波的成像讨论，给的模型中有一个盆地，带一个地壳界面，理论接收函数中可以看到多次波的存在，如果你的速度模型是正确的话，逆时偏移的结果都能找到正确的界面深度，假如你知道这个盆地了，然后你进行成像这没问题，但如果你不知道上面有个低速的盆地存在，直接使用一个平均的速度去成像的话，底下两个界面都错了，变深了，这就是盆地的影响，它会把界面拉深，这做成像都是知道的。　　  

　　六、结论　　  

　　关于今天报告的结论是，我们发展了一个有效的远震接收函数三维逆时偏移方法，逆时偏移这两年在天然地震中有不少新的工作发表，我们追踪了那些文献，这些文献基本都没用在接收函数上，主要是应用在远震波形上。同时这个方法可以有效提升计算效率。合成数据的理论测试表明这个方法是正确的，一个优点是它能够处理横向介质变化较大造成的衍射现象，而且对界面的倾角没有限定。最后我的建议是，被动震源研究布设的台站最好把台间距控制在5 km以内，也最好是面状的，没有条件的话也最好让这条测线与构造线垂直，如果你觉得这个地区有很强的横向变化的话，5 km可能还不够，那就再加密，越密越好，最后谢谢大家！　　  

　　【参考文献】　　  

　　[1] Dueker K G, Sheehan A F. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track[J]. Journal of Geophysical Research: Solid Earth, 1997, 102(B4): 8313-8327..（原文链接） 

　　[2] Zhu L, Kanamori H. Moho depth variation in southern California from teleseismic receiver functions[J]. Journal of Geophysical Research: Solid Earth, 2000, 105(B2): 2969-2980.（原文链接） 

　　[3] Kosarev G, Kind R, Sobolev S V, et al. Seismic evidence for a detached Indian lithospheric mantle beneath Tibet[J]. Science, 1999, 283(5406): 1306-1309.（原文链接） 

　　[4] Chen L, Wen L, Zheng T. A wave equation migration method for receiver function imaging: 1. Theory[J]. Journal of Geophysical Research: Solid Earth, 2005, 110(B11).DOI:10.1029/2005JB003665.（原文链接）