中国东北地区地幔转换带接收函数三维Kirchhoff偏移成像研究

2023-01-10朱敏吴庆举宁杰远张瑞青

地球物理学报 2023年1期

朱敏，吴庆举，宁杰远，张瑞青

1 中国地震局地球物理研究所，北京 100081 2 北京大学理论与应用地球物理研究所，北京 100871 3 中国地震局地震观测与地球物理成像重点实验室，北京 100081

0 引言

中国东北地区位于兴蒙造山带的东端，南北分别与华北克拉通和西伯利亚克拉通相接，东临西太平洋俯冲带，自古生代以来，先后经历古亚洲洋和蒙古—鄂霍茨克洋闭合的影响以及太平洋板块俯冲的多期构造作用(林强等, 1998; 吴福元和曹林, 1999)，这些构造作用对东北地区深部结构演化产生了重要影响.

东北地区广泛发育了一系列的新生代板内火山，其中比较活跃的是长白山火山和五大连池火山.前人利用多种地球物理方法对火山区的深部结构展开研究(Ai et al., 2003; Lei et al., 2013; 张风雪等, 2014; 潘佳铁等, 2014; Zhang et al., 2014; 强正阳和吴庆举, 2015; Li et al., 2016; Liu et al., 2017)，研究表明这些板内火山活动的形成与太平洋板块的俯冲作用存在重要的联系，然而这些火山作用的起源和形成机制仍然是一个存在争议的问题(Zhao, 2021).大地幔楔模型的研究认为板片的俯冲和后撤驱动的地幔楔中的地幔对流是这些板内火山的深部热源(Lei and Zhao, 2005; Duan et al., 2009; Tian et al., 2016; 田有等, 2019)；板片间隙模型的研究认为地幔转换带中停滞的俯冲板片存在空缺，由俯冲板片携带的热物质或者过渡带之下的热物质从这个空缺上涌为长白山火山提供热源(Tang et al., 2014; Liu et al., 2015; Guo et al., 2018; Tao et al., 2018)；最近的数值模拟研究表明，东北亚的板内火山活动可以用俯冲太平洋板块与含水MTZ的相互作用来解释(Yang and Faccenda, 2020).

对于中国东北地区地幔转换带结构的研究，有助于我们研究地幔转换带内俯冲板片的结构及其俯冲过程，进一步探究东北地区火山活动的深部成因.地幔转换带上界面(简记为410-km)和地幔转换带下界面(简记为660-km)分别对应不同矿物的相变面，其中410-km是橄榄石到瓦兹利石的相变面(Ringwood, 1968)，对应的克拉伯龙斜率(Clapeyron slope)为正，而660-km是林伍德石分解为钙钛矿和铁方镁石的相变面(Ito and Takahashi, 1989)，对应的克拉伯龙斜率为负.根据这一相反变化的克拉伯龙斜率性质，若MTZ附近存在高温物质(热流等)，410-km会变深，660-km会变浅，MTZ的厚度会减薄；若MTZ附近存在低温异常例如残存的俯冲板片或岩石圈拆沉物质，410-km会抬升，660-km会下沉，MTZ的厚度会增厚，因此MTZ的厚度是地幔转换带附近温度和物质组分的重要表征参数.

接收函数是研究速度间断面的有效手段，在全球深部结构研究中发挥着重要的作用.在接收函数成像方法中比较常用的是共转换点(Common Conversion Point, CCP)叠加方法(Dueker and Sheehan, 1997; Zhu, 2000)，这个方法通常基于一维速度模型，并且依赖水平界面的假设，对于速度横向不均匀性明显和存在大角度倾斜界面的地区，难以准确成像.借鉴反射地震学中的偏移思想，基于射线理论和波动方程的接收函数偏移成像方法也逐渐发展起来，并广泛应用于不同地区的间断面结构成像中.Bostock 和 Rondenay(1999)通过Born近似的引入实现了接收函数的波动方程偏移；Levander 等(2005)发展了Kirchhoff积分的散射波波动方程偏移方法；Chen 等(2005)提出了基于波动方法的叠后偏移方法；此外，Shang 等(2012)尝试将逆时偏移方法(Revers-time Migration, RTM)引入接收函数偏移成像，并得到了一定的发展(Burdick et al., 2014; Shang et al., 2017; Li et al., 2018)；Liu和Levander(2013)在Rondenay等(2001)基础上提出了三维广义Radon变换(Generalized Radon Transform, GRT)方法并验证其在火山口及圆顶状界面成像中的有效性.波动方程叠加、RTM、GRT等方法对于复杂界面可以恢复，但是计算成本较高，并且对于数据的要求也较高，大多运用在区域台阵或线性阵列研究中.我们仍然有必要发展新的数据处理手段，以挖掘现有数据的潜力.

近年来，基于散射核的接收函数偏移成像方法也得到进一步发展.Cheng等(2016)基于前人Kirchhoff偏移思想，结合Fast Marching Method(FMM)三维走时场计算方法，提出了接收函数Kirchhoff 3-D叠前偏移成像方法，通过理论测试验证了该方法对高角度倾斜界面的成像有效性，并应用到北美Cascadia俯冲板片的研究中(Cheng et al., 2017).Hansen 和Schmandt(2017)利用散射核的计算实现P波接收函数和S波接收函数的三维偏移成像.在此方法基础上，Zhang和Schmandt(2019)通过引入慢度窗口，减少了由于数据密度不够造成的成像假象.Millet等(2019)进一步发展了联合多种散射震相的接收函数Kirchhoff偏移成像方法.

本文拟在前人提出的接收函数Kirchhoff偏移成像方法基础上，综合考虑多种叠加因素，进一步实现接收函数三维Kirchhoff偏移成像方法，并利用中国东北地区丰富的接收函数资料，得到中国东北地区(116°E—134°E，38°N—52°N)地幔转换带间断面的起伏形态，进而为研究区板内火山的深部成因提供约束.

1 研究数据与研究方法

1.1 研究数据

本研究中用到的接收函数由Zhang等(2016)提供，主要由三部分组成：研究区东北近乎平行的两条线性台阵(Linear arrays)由中国地震局地球物理研究所布设于2009年6月至2011年8月，包括121个宽频带地震台站，台间距约为20 km；几乎同时期中美国际合作布设的NECESSArray台阵由127个宽频带流动地震台站组成，平均台间距约为80 km；穿插其中的124个固定台(CEA stations)的记录时间是2009年5月至2013年5月(Zheng et al., 2010).整体上，研究区布设共计372个地震台站，记录到约1200个地震事件，为我们研究中国东北地区(116°E—136°E, 38°N—52°N)上地幔间断面结构成像提供了良好的数据基础(图1).

图1 研究区示意图及地震台站和地震事件分布情况

接收函数由时间域迭代反褶积方法(Ligorría and Ammon, 1999)提取得到，由于东北地区地幔转换带结构复杂，为了对小尺度异常进行成像，我们选取高斯系数2.5，对应的截止频率约1.2 Hz.对收集到的接收函数进一步筛选，选出0时刻附近(约0.5 s)振幅积分为正、且该振幅为整条接收函数的最大振幅、后续震相清晰的6万多条接收函数用于偏移成像研究(图2).

图2 P波径向接收函数实例

1.2 研究方法

1.2.1 Kirchhoff叠前偏移成像

基于Kirchhoff积分的偏移成像方法最早广泛应用在石油行业中，由Levander等(2005)引入天然地震成像研究中.该方法的基本思路是通过给定背景速度模型，把地震波能量反传到地震波等时双曲线上.对于Ps接收函数而言，研究区三维空间的所有网格点都可以当作潜在的P-to-S转换点，通过给定的转换波与直达波的走时差把接收函数的能量偏移到对应的深度双曲线上，进而对每一条观测记录的能量进行加权叠加得到地下界面成像结果.图3给出了Kirchhoff散射波场成像的基本框架.

图3 散射波场示意图(改自Millet et al., 2019)

根据Rondenay(2009)总结的转换波偏移成像的基本思路，偏移方程可以表示为

fPS(x)=∑i∑jW(i,j,k)·RF(i,j,TPS(i,j,x)),

(1)

其中，fPS(x)是成像格点x处的散射势能，i和j分别代表第j个台站和第i个地震事件，RF表示接收函数，本研究中只用到P波径向接收函数.W表示相应的叠加因子，包含几何扩散、散射花样、偏振极性等信息，与事件和台站的相对位置、散射角有关，将在后文中详细展开.TPS(i,j,x)=TP(i,x)+TS(j,x)-Td(i,j)表示走时差，其中TP(i,x)是第i个地震到散射点x的P波走时，TS(j,x)表示从转换点x到第j个台站的S波走时，而Td(i,j)表示参考速度模型下第i个震源到第j个台站的直达P波的走时.

三维走时场的计算我们利用基于快速行进法(Fast Marching Method, FMM)的FM3D程序(De Kool et al., 2006)实现，这个方法的核心思想是利用由节点组成的波前窄带模拟曲面的演化，通过向后差分获得程函方程的弱解，与基于斯奈尔定律的传统射线追踪方法相比，可以得到空间中所有节点的完整走时场，且可以同时计算多个接收点的走时，计算效率较高.在实际计算过程中，我们只需要计算由震源出发、台站接收的正传P波走时场和由台站出发(以台站为虚拟震源)、散射点接收的反传S波走时场，计算量由地震事件个数、台站个数以及研究区三维网格划分尺度决定.

权重因子进一步可以写为

W(i,j,k)=G(j,x)·εPs(i,j,x)·δPs(x),

(2)

其中G(j,x)是散射点x到第j个台站的几何扩散因子，在三维参考坐标下，G(j,x)=1/d(j,x)，d(j,x)近似为散射点与台站之间的直线距离.εPs(i,j,x)为转换波的散射花样，根据Rondenay(2009)的推导，可以表示为

(3)

综合以上的叠加因子，研究区下方每个散射格点的Ps转换波散射势能可以表示为

(4)

1.2.2 理论测试

前人为了验证此方法对于倾斜界面的分辨能力进行了诸多理论测试(Cheng et al., 2016; Hansen and Schmandt, 2017; Millet et al., 2019; Zhang and Schmandt, 2019)，测试结果表明相较于CCP共转换点叠加方法而言，基于Kirchhoff叠前偏移的接收函数成像方法可以考虑速度的横向不均匀性，对倾斜角度高达60°的界面都有很好的恢复能力，因此对于复杂形态的界面有更好的分辨能力.为了验证接收函数三维Kirchhoff偏移成像方法在中国东北复杂俯冲构造下的成像可能性，我们设计了二维板片俯冲停滞模型进行测试.

采用的参考速度模型是AK135模型，其中高速俯冲板片从地幔过渡带之上以30°倾角俯冲进入地幔过渡带，并在660-km界面上方水平停滞，其下界面下沉至710 km深度.板片厚度为100 km，板片的P波、S波速度以及密度异常均为+5%左右(图5)，与前人在这一地区速度成像得到的速度异常幅度相当(Tao et al., 2018).根据东北地区实际数据的台站台间距(20～80 km)，设计台间距为0.3弧度，共计101个台站线性排列在地表，以台站中心为基准，在台站两侧，50°～90°震中距10°等间距分布共10个远震事件(图4a).本测试主要对板片俯冲结构进行成像，为了成像结果更加清晰直观，设置地震事件均为深源地震，震源深度为550 km，以保证pP，sP波等深度震相至少比P660s波晚到20 s，减少震相的混杂对成像结果的干扰.

利用二维谱元法(SPECFEM2D)(Komatitsch et al., 2001)合成理论地震图，最小分辨周期为1 s左右，数据采样间隔为0.025 s，长度为1200 s，采用中心频率为0.3 Hz的Ricker子波作为震源时间函数(图4b).选取震中距30°～90°的合成数据经坐标系旋转，时间域迭代提取高斯系数为2.0的890条理论P波接收函数(图6).

图4 震源与台站的分布情况及震源时间函数

图5 研究区速度模型示意图

图6 正演计算的接收函数实例

基于图5的参考速度模型光滑后计算三维走时场并进行偏移叠加成像，得到如图7所示的成像结果.结果显示，Kirchhoff三维偏移成像方法可以对地幔转换带上下界面(410 km和660 km)进行清晰的成像，并且对于以30°倾角俯冲进入地幔转换带并水平停滞在地幔转换带中的俯冲板片，其上下界面的形态也有较好的成像效果.

图7 研究区200～798 km偏移成像结果(以整个剖面最大绝对幅值做归一化)

2 成像结果

前人在中国东北地区地壳上地幔开展了一系列的速度成像研究(Zhao, 2004; Huang and Zhao, 2006; Tang et al., 2014; 张风雪等, 2014; 潘佳铁等, 2014; Shen et al., 2016; Tao et al., 2018; 田有等, 2019)，结果显示，中国东北地区存在明显的速度横向不均匀性，主要表现在长白山下方显示明显低速异常，松辽盆地表现为显著高速异常，速度异常可达±6%，因而在利用接收函数进行偏移成像时要考虑研究区三维速度横向不均匀性的影响.本研究采用Tao等(2018)利用全波形反演的方法得到的三维速度模型FWEA18模型，这个模型同时利用了体波及面波的走时和振幅信息，是东亚地区较为精细的地壳上地幔三维P波以及S波速度模型，模型的横向分辨率为50～80 km，纵向分辨率为20～50 km.

以FWEA18模型为参考，基于三维Kirchhoff偏移成像方法对研究区(115°E—136°E，36°N—52°N)采用0.1°×0.1°×3 km的网格划分，计算走时场，偏移叠加得到研究区深部间断面结构，由于台站的台间距高达80 km，而地表至两倍台间距深度成像结果会受空间假频较大影响(Rondenay et al., 2005)，因此我们主要讨论地幔转换带的结构.图8给出了研究区410-km界面以及660-km界面转换波射线穿透点的密度，在研究区大部分地区射线密度都有很好的覆盖.

图8 研究区410-km界面和660-km界面的穿透点密度

图9展示了41°N至49°N每一度间隔共九条东西向剖面的成像结果，每条剖面以最大绝对振幅进行归一化处理，颜色深浅只表征散射势能的相对大小.结果显示，在410 km和660 km深度附近存在两条连续性较好、正极性(速度由小到大变化)的强信号，分别代表410-km和660-km间断面.从各剖面结果中可以看出，410-km界面整体起伏幅度较小，平均深度在410 km左右，而660-km界面深度变化较大，并且以45°N为界呈现明显的南北差异，在41°N—45°N，660-km界面下沉明显，最大下沉在43°N剖面中130°E附近，深度可达711 km，结合剖面中深源地震的分布，此处的剧烈下沉可能与俯冲板片与660-km间断面在此处直接交汇有关.在45°N以北，660-km界面的下沉幅度减小，范围向西扩展至125°E附近，呈现整体性下沉的特点.Wang等(2020) 的接收函数剖面结果显示，在410 km与660 km深度之间还存在一组极性相反的震相，推测是过渡带内停滞的俯冲板片的上、下界面.本文剖面结果中也可以观测到地幔转换带中存在一些不连续的界面，但是这些界面的埋深及展布形态都与之存在一定差异。因此，这些不连续界面是否是俯冲板片的上、下界面仍值得进一步深入研究.

图9 研究区东西方向41°N—49°N共9条垂直剖面250 km至790 km散射势能以每条剖面上最大绝对幅值做归一化进行绘图.图中绿色虚线是410-km和660-km界面的参考位置，深度分别是410 km和675 km.

根据各剖面410-km和660-km界面最大散射势能深度的拾取得到了研究区410-km以及660-km界面的空间分布，并根据这两个界面的深度差得到了地幔转换带厚度的变化(图10a,10c,10e).410-km界面的深度范围为391～430 km，平均深度为411 km.在45°N两侧差异明显，45°N以北主要表现在松辽盆地北侧五大连池火山群附近大范围抬升，抬升幅度达10 km；45°N以南表现为几处明显的下沉，集中在长白山火山附近以及松辽盆地以南，下沉幅度达到15～20 km.660-km界面的平均深度为675 km，深度范围为651～711 km，长白山火山群附近显示大范围明显下沉，下沉超过25 km，沿北偏西近30°方向呈现大范围不连续的下沉，下沉幅度稍小，约15 km范围内，此外，在研究区东南角也呈现小范围10 km左右的下沉.研究区MTZ厚度从250 km到313 km变化，与660-km的结果呈现相似的特点.在后文中，我们以410 km和675 km作为410-km界面和660-km界面的参考深度，并以265 km作为MTZ的参考厚度对研究区火山群附近界面的起伏以及地幔转换带的厚度分布进行讨论.

图10 研究区的410-km和660-km界面深度以及地幔转换带厚度变化

结合图10中410-km界面、660-km界面起伏和MTZ厚度的分布以及图11中这两个界面深度以及MTZ厚度之间相关关系可以看出，410-km界面和660-km界面的深度呈较弱的正相关；MTZ的厚度与410-km界面深度呈现明显负相关，与660-km界面深度呈明显的正相关，受这两个界面深度变化的共同影响.

图11 410-km界面深度、660-km界面深度、MTZ厚度的相关关系蓝色圆圈是两两对应散点，紫色实线是线性拟合关系. (a)660-km界面深度与410-km界面深度的相关关系； (b) MTZ厚度与410-km界面深度的相关关系； (c) MTZ厚度与660-km界面深度的相关关系.

3 讨论

以三维速度模型FWEA18模型为参考模型，我们得到了中国东北地区地幔转换带界面的起伏变化，为了进一步分析三维速度横向不均匀性的影响，我们利用全球一维模型IASP91(Kennett and Engdahl, 1991)采用相同的处理流程进行成像，并与基于三维速度模型的成像结果比较.

图10中左右两列分别为基于FWEA18模型的成像结果(图10a,10c,10e)以及基于IASP91模型的成像结果(图10b,10d,10f).结果显示，基于这两个模型的成像结果主要特征相似，在细节上存在一定差异.410-km界面的成像结果，相比于一维模型，基于FWEA18模型的成像结果没有显示45°N左右松辽盆地西侧明显的抬升，并且长白山火山周围的下沉异常范围收缩到长白山以北、龙岗火山附近，在长白山正下方没有明显异常，此外在松辽盆地以南呈现了局部的下沉.660-km界面的成像结果中，基于FWEA18模型的结果显示诺敏河火山以北的异常不明显，五大连池和长白山火山附近660-km下沉的范围收缩，呈现不连续的特征，长白山火山以西的C异常范围更大，抬升更明显.总体来说，FWEA18模型的结果异常范围更小，主要火山区的异常不连续，长白山火山以西的异常更加明显.对于结果的差异，可以归因于三维速度对走时差的校正在一定程度上修正了各界面的深度，考虑速度横向不均匀性、基于三维速度模型的成像结果能够更加准确的恢复东北地区地幔过渡带结构，在后文中主要对基于FWEA18模型的成像结果进行讨论.

410-km间断面和660-km间断面的起伏形态是地幔转换带物质组分、热结构、含水量等因素的综合表现.根据这两个速度密度间断面深度的分布以及地幔转换带厚度的变化，我们对研究区深部结构进行讨论，并尝试对主要火山区的深部成因作出解释.

以675 km作为660-km界面的参考深度，该界面表现为大范围的下沉(图10中A、B异常)和局部的隆起(图10中C异常).界面下沉最明显的是位于长白山火山群附近的A异常，异常幅度超过25 km，相同位置410-km界面表现为小范围的下沉，对应的MTZ的A异常也显示明显下沉，前人接收函数研究(Ai et al., 2003; Liu et al., 2015; Zhang et al., 2016)也显示相似的结果，并解释这一异常与较冷的西太平洋俯冲板片有关.C异常位于镜泊湖火山以西200 km位置，660-km界面显示近15 km的抬升，与Liu等(2015)CCP成像得到的结果一致，并且在FWEA18模型中同样深度相近位置，S波也表现为明显的低速异常，暗示该处存在高温异常.结合Tang等(2014)提出的板片空缺模型，我们认为俯冲停滞的板片可能由于重力失稳发生撕裂，形成热物质上涌的通道，热物质上涌并与地幔楔流汇聚至410-km界面下方，造成410-km界面的下沉，由于板块后撤的影响使得上升热流向东倾斜，为长白山火山群提供深部热源.

在C异常以北存在呈长条状南北分布的异常B，表现为660-km界面15 km左右的下沉，较A异常幅度稍小，可能是俯冲板片的残存造成的低温异常，在410-km界面中表现为五大连池以及诺敏河火山群附近的明显抬升，与Zhang等(2016)的结果类似，结合前人岩石圈上地幔的速度成像结果，松辽盆地下方存在高速异常，可能是由于松辽盆地下方岩石圈发生拆沉(Li et al., 2012; 潘佳铁等, 2014; 张风雪等，2014)，较冷的岩石圈拆沉物质下沉至410-km附近，造成410-km界面的抬升，同时岩石圈拆沉造成的软流圈局部热流上涌为五大连池火山、诺敏河火山提供深部热源.Guo 等(2018)联合体波和面波数据反演得到的速度结构也显示松辽盆地周围阿尔山火山以及阿巴嘎火山下方的低速异常深度较浅，不超过200 km，认为这些火山的的岩浆来源于局部的软流圈热物质上涌，我们的观测支持这一观点.

除这几处主要异常外，松辽盆地以南还存在不连续的660-km以及410-km的下沉，结合长白山火山以南660-km界面没有明显的下沉，我们认为西太平洋俯冲板片的主体可能已经进入下地幔，残存的俯冲板片造成这些局部幅值较小的间断面异常.