苏晓康, 赵明辉, 李子正, 袁野, 王星月, 程锦辉, 张佳政

1. 中国科学院边缘海与大洋地质重点实验室(南海海洋研究所), 南海生态环境工程创新研究院, 广东 广州 510301;

2. 南方海洋科学与工程广东省实验室(广州), 广东 广州 511458;

3. 中国科学院大学, 北京 100049

南海是西太平洋的一个大型边缘海, 形成于自晚白垩世起始的渐进式裂谷作用(Taylor et al, 1983;Franke et al, 2011); 位置特殊, 夹在欧亚、印澳、菲律宾海三大板块之间, 构造环境十分复杂, 在其发展演化的一个较短的时间里经历了近乎完整的威尔逊旋回, 包括张裂、破裂、海底扩张、俯冲消减(Ludwig et al, 1979; 李家彪 等, 2004; 李春峰 等,2009); 具有规模小、年龄新的特性, 条件较为优越,是我国走向深海大洋研究的重要突破口(汪品先,2009)。

洋陆转换带(Continent–Ocean Transition, COT)是位于减薄陆壳和洋壳之间的一个过渡地带, 是伸展背景下大陆岩石圈与大洋岩石圈相互作用的特殊区域(郭玲莉 等, 2017), 蕴藏着岩石圈张裂、破裂及海底扩张过程的大量信息(Whitmarsh et al, 2001);其深部地壳结构、物质组成及构造属性是探究南海形成演化、解决北部陆缘性质的关键。被动大陆边缘被划分为贫岩浆型(非火山型)和富岩浆型(火山型)。其中贫岩浆型(非火山型)被动大陆边缘的主要特征有发育拆离断层和向海倾斜的旋转断块、地壳发生强烈拉张减薄、洋陆转换带较宽(宽度可超过200km)且剥露蛇纹石化地幔(Dean et al, 2000;Pérez-Gussinyé et al, 2006; Lavier et al, 2006; Reston,2007, 2009); 富岩浆型(火山型)被动大陆边缘的主要特征有洋陆转换带发育向海倾斜反射层(Seaward Dipping Reflector sequence, SDRs)、存在下地壳高速层(Geoffroy, 2005; Schnabel et al, 2008; Geoffroy et al, 2015)。

在前人的一些研究中, 南海北部陆缘不仅拆离断层发育, 而且在其洋陆转换带处并未发现富岩浆型边缘特有的向海倾斜反射层, 同时发育大量海倾的断层系统(Taylor et al, 1983; Zhu et al, 2012; Gao et al, 2016)。因此, 有人认为南海北部陆缘属于贫岩浆型(非火山型)大陆边缘(Yan et al, 2001; Ding et al,2012)。但是, 南海北部陆缘也表现出富岩浆型(火山型)大陆边缘的一些特征, 比如张裂后岩浆活动活跃(Gao et al, 2015; Sun et al, 2020), 在其洋陆转换带处存在下地壳高速层(Wang et al, 2006; Zhao et al, 2010;Liu et al, 2018; Zhu et al, 2018; Wan et al, 2019); 平均伸展程度比伊比利亚非火山型陆缘要小得多(Clift et al, 2001); 从张裂阶段到海底扩张阶段比较迅速(Larsen et al, 2018a), 洋陆转换带比较狭窄, 仅宽15～25km (Ding et al, 2020)等。另外, Sun等(2019)认为南海北部陆缘的岩浆活动具有阶段性, 先期表现出贫岩浆的特征而后期表现出富岩浆的特征。Gao等(2015)通过多道地震数据和重力模拟研究发现下地壳高速层在南海北部陆缘的分布具有特殊性,仅分布于东段。所以也有人认为南海北部陆缘是介于贫岩浆型(非火山型)与富岩浆型(火山型)之间的中间型(Clift et al, 2001; Zhu et al, 2012; Larsen et al,2018a)。

南海北部陆缘兼具贫岩浆型(非火山型)与富岩浆型(火山型)陆缘的特征, 具有特殊性与复杂性。为了更加深入地了解南海, 揭示其从张裂到破裂的演化机制, 国际大洋发现计划(International Ocean Discovery Program, IODP)在南海北部陆缘开展的大洋钻探IODP367/368/368X航次, 从陆向洋共在7个站位(图1)进行了钻探, 其中站位U1499钻遇同裂谷或前裂谷的砾岩, 但未钻遇基底(图2a); U1502钻到基底, 但钻遇的是强烈热液蚀变的玄武岩(图2a),具有洋中脊型玄武岩的成分组成(Larsen et al,2018a), 并没有钻遇蛇纹石化地幔(Sun et al, 2019)。任金锋等(2018)通过古地貌重建认为, U1499钻遇的砾石层发育于小型断陷陆相湖盆内, 站位U1499处的地壳属于陆壳或洋陆转换带。Larsen等(2018a)由U1499钻遇的砾岩有较大比例的上地壳物质, 部分是沉积成因, 认为其下应属于陆壳范围。那么, 从站位U1499到U1502, 在如此窄的距离(40km)如何完成了从陆壳到玄武岩洋壳的转变？另外, 站位U1502处钻遇的玄武岩被强烈热液蚀变, 其成因是什么？板块破裂早期的局部岩浆活动？以上的钻探结果表明, 钻探并不能全面揭示南海北部陆缘洋陆转换带COT的地壳性质、张裂–破裂机制与动力学过程, 需要地球物理探测获取更多深部的证据(Jian et al, 2018)。

图1 南海北部陆缘OBS2018–H2测线位置与水深地形图图a中的黑色方框表示图b的位置; 图b中黑色粗实线为OBS2018–H2测线; 黑色箭头代表放炮方向; 黄色五角星代表IODP钻探站位; 红色实心圆代表OBS台站; 橙色实心矩形、圆形和三角形分别代表磁异常C11、C10、C9 [据Briais等(1993)]; 水深数据来自GEBECO_2019, http://www.geobco.net。审图号: GS(2016)1609Fig. 1 Location of OBS2018–H2 and the bathymetry map of the northern continental margin of the South China Sea. The black square in Figure a indicates the location of Figure b. The thick solid black line indicates OBS2018–H2. The black arrow indicates the shooting direction. The yellow pentagrams indicate the IODP drilling sites. The red solid circles indicate the OBS sites. Orange solid rectangles, circles and triangles indicate magnetic anomalies C11, C10 and C9, respectively (after Briais et al, 1993). Bathymetric data are from GEBECO_2019 (http://www.geobco.net)

图2 IODP站位U1499和U1502的钻探结果a. 站位U1499和U1502的岩性柱状图[据Larsen等(2018a)修改]; b. 站位U1502的速度间隔模型[据Larsen等(2018b)修改], 空心小正方形和小三角形代表纵波测量仪PWC (P-wave caliper)数据, 空心小圆代表测井数据Fig. 2 Drilling results of IODP drilling sites U1499 and U1502. (a) Stratigraphic columns of U1499 and U1502 (modified after Larsen et al, 2018a). (b) The velocity interval model of U1502 (modified after Larsen et al, 2018b). Small squares and triangles indicate PWC (P-wave caliper) data, and small circles indicate logging data

2018年实施的OBS2018–H2测线穿过了IODP站位U1499和U1502 (图1), 深部地球物理探测与钻探相结合, 深浅优势互补, 可望解决从站位U1499到U1502如何在短距离实现洋陆转换及其深部结构特征如何变化的科学问题。本文以OBS2018–H2测线数据为依托, 详细叙述了前期数据处理过程、初始模型的建立、震相的初步识别, 并对地壳结构进行了初步分析。

1 海上数据采集

2018年5—7月, 通过搭载中国科学院南海海洋研究所的“实验2”号科考船, 国家自然科学基金委“南海北部地球物理共享航次”顺利完成了南海北部洋陆转换带COT区域大规模三维深地震探测实验(杨富东 等, 2020)。实验旨在获得该区域的深部速度结构, 与IODP结果共同揭示南海的形成演化历史。

本次实验使用的是中国科学院地质与地球物理研究所研制的便携式短周期海底地震仪(Ocean Bottom Seismometer, OBS)(郝天珧 等, 2011), 其主要构成部件包括内置三份量检波器、外置的深海水听器、采集电路板、工作状态监控模块和内置Wifi模块。震源系统为中国科学院南海海洋研究所的4支大容量Bolt枪组成的气枪阵列, 总容量达6000in3(98.32L), 激发的气枪信号的主频为4～8Hz, 具有衰减较慢、传播距离较远的特征(赵明辉 等, 2008)。

本文重点讨论OBS2018–H2测线, 它是三维地震探测实验中的一条近东西向的主测线, 共投放OBS 10台, 采样率为100Hz, 最终成功回收10台,回收率达到100%。测线长度约179km, 穿过了IODP站位U1499和U1502, 激发的有效炮数为849炮, 放炮时间间隔为90s, 放炮期间船速约为5节。另外,气枪放炮的同时, 随船还进行了沿测线的多道地震数据的采集。该多道地震数据共24道, 道间距为12.5m, 多道拖缆长度为300m, 采样时间间隔为2ms, 但其中20道数据损坏, 只余4道数据可用。

2 数据处理

海上采集的数据包括Hypack导航数据、计时器数据、OBS原始数据及多道地震(Multi-channel seismic reflection, MCS)数据。这些原始数据不能直接使用, 需要经过一系列的处理, 包括数据格式转换、炮点和OBS位置校正、多道地震数据处理等。

2.1 数据转换处理

2.1.1 OBS2018-H2测线导航数据的整理转换

首先提取原始导航数据(Hypack文件和计时器文件)中记录的气枪每次激发的炮号、精确的放炮时间(GMT时间)、经纬度坐标、航向及航速等信息, 然后根据标准格式输出标准导航文件UKOOA, 用于后期按炮对OBS数据进行裁截处理(赵明辉 等,2004)。

2.1.2 OBS记录数据的格式转换

首先对OBS原始数据进行解编处理, 转化为SAC格式(张佳政 等, 2018), 得到4个分量SAC格式数据(分别为垂直、水平1、水平2和水听器)。国产OBS工作时把数据分段存储在定时生成的每个文件里, 所以还需要把放炮时间内的所有SAC数据进行合并。然后, 使用地震分析软件SAC (William et al, 1991)对得到的有连续时间序列的SAC格式数据进行分析, 查找气枪信号的绝对时间。最后, 从前面得到的标准炮点导航文件UKOOA中读取放炮时间,并截取SAC格式数据中的有效信号段转换成国际通用的SEGY格式。SEGY文件可使用SU软件(Stockwell,1999)进行可视化处理以检查数据质量(图3)。

图3 OBS数据处理流程图Fig. 3 Flow chart of OBS data processing

2.2 炮点和OBS位置校正

由于波浪、海流及船给的水平初速度等因素的影响, OBS投放后在海水中并不是以类似自由落体的形式垂直下落, 而是会从投放位置运动漂移数百甚至上千米(敖威 等, 2010; 张莉 等, 2013)。另外,由于海水表层流及海风的作用, 船体无法按设计的测线以直线行驶, 从而造成放炮轨迹是一条曲线。所以, 需要进行炮点和OBS位置校正。

2.2.1 炮点位置校正

在放炮过程中, 差分全球定位系统(Differential Global Position System, DGPS)授时由Haypack记录到的炮点位置实际应该是DGPS的位置而并非气枪震源所在位置。产生的偏离包括“实验2”号科考船上DGPS点到船尾的硬连接(17m)和由船尾到水中气枪阵列的软连接(35m)。科考船行进过程中,船体受风浪影响会左右摇摆曲线前行, 软连接部分无法随船体实时转向, 从而造成DGPS与船尾的连线和船尾与气枪震源的连线并不平行, 而是形成一定的角度(杨富东 等, 2020)。对于以上问题,杨富东等(2020)提出了利用连续3炮的平均航向来校正炮点位置, 可以最大程度贴合实际放炮轨迹。当炮点为x时, 将其与前一炮(x–1)和后一炮(x+1)分别连接, 得到两个连线航向, 将两个航向平均即为炮点x的航向。假设该平均航向为θ, 校正前的炮点位置为A点(x0,y0), 校正后的炮点位置为B点(x1,y1), 则校正公式如下:

所有炮点校正到实际位置后, 通过Matlab程序拟合为一条直线, 以该直线作为后期速度模拟的测线, 然后将所有炮点均投影到直线上得到新的炮点位置。校正前后炮点位置对比如图4所示, 以50炮为间隔显示。由于OBS2018–H2测线由东向西放炮,船向西行驶, 炮点位置由船上DGPS记录点校正到船后的气枪阵列处, 因此, 图4中可以看到校正后的炮点沿测线向东移动。

图4 炮点和OBS位置校正前后对比图炮点对应的数字为其序号, 每隔50炮突出显示一次; OBS对应的为其代号; 炮点校正后的OBS2018–H2测线: y=0.1243x+4.0068Fig. 4 Shot point locations and OBS locations before and after relocation. The number corresponding to the shot point is its serial number. Highlight every 50 shots. OBS corresponds to its code name. OBS2018–H2 after shot point relocation:y=0.1243x+4.0068

2.2.2 OBS位置校正

OBS实际海底落点与船上记录的投放点有一定的偏离, 如果把投放点位置当作接收点的位置, 后期拾取的震相会有比较大的偏差, 所以需要对OBS的位置进行校正。直达水波震相(Pw)是气枪信号激发后在海水中传播, 直接由OBS检波器记录, 其炮点到OBS的时-距信息最简单可靠, 常用来校正OBS位置。

蒙特卡洛法是最常用的校正方法, 曾应用于多次OBS位置校正工作(敖威 等, 2010; 张莉 等,2013; Du et al, 2018; 杨富东 等, 2020)。具体校正步骤如图3所示。所有台站的OBS位置校正结果见表1, 校正前后的位置见图4。

表1 OBS2018–H2测线10台OBS位置校正结果表Tab. 1 Relocation results of 10 OBS in OBS2018–H2

位置校正的重要判断标准是直达水波在校正后是否左右对称(敖威 等, 2010; 张莉 等, 2013; Du et al, 2018; 杨富东 等, 2020)。以OBS2018–H2测线中OBS37为例。从图5中可以看到, 校正前后, 直达水波形态从左右不对称变为对称形, 说明校正后的OBS位置更接近其实际落点, 校正效果良好。图4展示了OBS校正后的台站位置与其投放点之间的关系, 10个台站漂移距离均在400m范围内(表1),远小于在西南印度洋中脊实验区的OBS位置校正结果(敖威 等, 2010), 其最大偏移达3590m。通过两次实验过程推测, OBS2018–H2测线实验实施过程中的良好海况, 可能是其OBS偏移距离远小于西南印度洋实验的重要因素之一。另外, 位置校正结果中, 走时残差均方根(Root Mean Square, RMS)值均在5ms以内, 校正精度可以为后期速度结构模拟提供足够的保障。

图5 OBS37台站直达水波震相Pw校正前(a)和校正后(b)对比折合速度为6.0km·s-1, 折合时间=绝对走时−偏移距/6Fig. 5 Contrast of direct water phases (Pw) before and after correction of OBS37. Reduced velocity is 6.0 km·s-1, and reduced travel time = travel time − offset/6

2.3 多道地震数据处理

OBS2018–H2测线在进行作业的同时, 随船采集了沿测线的多道反射地震数据。多道反射地震数据的初始文件是SEGY格式记录的, 使用SU软件包进行处理。多道地震的数据是分段独立记录的, 其初始文件炮号范围不能与目标测线对应一致, 所以第一步先进行裁截与拼接。首先读取SEGY文件将其转成SU文件, 然后经过计时器班报和UKOOA文件对比, 剔除无效炮, 按照对应的炮号进行裁截与拼接。第二步, 把4道叠加起来, 得到叠加后的H2–stack.su文件,之后把偏移距offset加入其道头。第三步, 对SU文件进行滤波、增益等处理后绘制出多道地震记录剖面图(图6a)。剖面的层位解释如图6b, 中间的黄色部分为多道剖面LW3, 来自Sun等(2016), 覆盖了OBS2018–H2的一部分区域。因此, 我们将二者结合起来, 共同为建立初始模型提供依据。

图6 多道反射地震剖面及其层位解释a. 随船采集的原始多道反射地震剖面; b. 多道反射地震剖面层位解释, 黄色实心圆代表OBS台站, 黑色竖直实线代表IODP钻探站位U1499和U1502, 粉红色、黑色和褐色虚线分别代表海底(Seafloor)、基底(Tg)、莫霍面(Moho), 绿色、蓝色、紫色虚线分别代表不整合T30、T60、T70, 红色短线段代表断层Fig. 6 Multi-channel reflection seismic profile and its horizon interpretation. (a) Original multi-channel reflection seismic profile acquired by Shiyan 2. (b) Horizon interpretation of multi-channel reflection seismic profile. The yellow dots indicate the OBS sites. The black vertical lines indicate the IODP drilling sites of U1499 and U1502. The pink, black and brown dashed lines indicate Seafloor, Tg and Moho, respectively. The green, blue and purple dashed lines indicate unconformity T30, T60 and T70, respectively. The red short lines indicate faults

3 初始速度模型的建立

多道地震(MCS)数据含有丰富而可靠的浅部信息, 可以很好地限制海底和基底的形态。因此, 这里以多道为基础建立初始速度模型。

多道上获取的层位信息是一个时间域的地质模型, 需要通过时深转换把它转成深度域的地质模型。首先, 根据经验设定, 水层速度取均值1.5km·s-1,海底面(Seafloor)下速度取1.6～1.7km·s-1, 基底(Tg)上速度取3.6～4.2km·s-1, 中间层位间隔取平均。然后, 参考IODP钻探站位U1502的速度间隔模型(图2b), 这里取Seafloor界面下速度为1.7km·s-1, Tg界面上速度为4.0km·s-1, 具体各层初始速度设定见图7。再后, 通过时深转换计算把时间域模型的各界面双程走时转换为深度。转换计算公式:Di=Di-1+（vi上+vi-1上）（t i-ti-1）/4, 各层参数对应关系见图7。具体计算顺序是自上而下, 逐层累加。例如, 界面T30的深度D2即为海底面深度D1加上海底面—界面T30的这个地层单元的厚度, 而地层单元的厚度可通过其平均速度与走时计算得出。之后, 手动设定上地壳(4.5～6.4km·s-1)、下地壳(6.4～7.0km·s-1)、上地幔(8.0～8.2km·s-1)及地壳厚度为10km。最后, 经过整理即可得到OBS2018-H2测线的初始速度模型(图8)。

图7 时深转换计算各参数对应图Fig. 7 Corresponding parameters of time-depth conversion calculation

图8 OBS2018–H2测线的初始速度模型黄色实心圆代表OBS台站; 沉积层内的曲线代表图6中的沉积界面T30、T60、T70; 虚线代表6.4km·s-1速度等值面Fig. 8 Initial velocity model of OBS2018-H2. The yellow dots indicate the OBS sites. The curves within the sediment represent the sedimentary interfaces T30, T60, and T70 in Fig. 6. Dashed line indicates 6.4 km·s-1 velocity iso-surface

4 震相初步识别与地壳结构初步分析

在完成了炮点和OBS位置校正之后, 重新进行数据转换生成新的SEGY数据。然后利用地震处理软件SU进行可视化。这里以OBS36、OBS37、OBS09和OBS42台站为例, 具体操作包括数据读入、自动增益控制、带通滤波, 最终得到综合地震记录剖面(赵明辉 等, 2004; 张佳政 等, 2018)。根据震相对比原则及各震相记录特征, 4个台站的综合地震剖面都能清晰识别出多组清晰的P波震相(图9), 如Pw、Pg、PmP 和Pn震相。其中, Pw是直达水波震相, 呈左右对称、倒置的双曲线形态, 常用于OBS位置校正。Pg是地壳内部折射震相, 是能量最强的一组震相, 可连续追踪, 在OBS36左半支6～30km和右半支6～36km、OBS37左半支6～42km和右半支6～40km、OBS09左半支6～45km和右半支6～45km以及OBS42左半支6～58km和右半支6～32km范围内都能连续追踪到该震相。PmP是莫霍面反射震相,与Pw形态相似, 呈左右对称的双曲线形态, 由于速度更大, 形态更为平缓, 在一定范围内滞后于Pg震相, 在OBS36右半支22～40km、OBS37左半支20～42km、OBS09左半支32～52km和右半支28～45km以及OBS42右半支17～20km范围内能够识别该震相。Pn是来自上地幔顶部的折射震相, 视速度比Pg大且往往与Pg相交, OBS09没有看到明显的Pn震相, 在OBS36右半支36～65km、OBS37左半支42～58km和右半支40～65km、OBS42右半支22～37km可以追踪到该震相。

测线西侧OBS36台站, 位于IODP站位U1499附近, 其下方的海底地形较为平坦(图1), 但两侧的Pg震相特征不同(图9a), 左半支在偏移距17km时,走时突然增加, 而右半支走时增加缓慢, 主要是由于两侧基底起伏情况不同所引起(图6)。根据Pg和Pn震相交点的偏移距可以初步估算地壳厚度, 如OBS36右半支Pg和Pn震相交点的偏移距为36km,取地壳的平均速度为6km·s-1, 则可以估算OBS36台站下方的地壳厚度, 为36/6=6km。OBS37台站位于IODP站位U1499和U1502之间, 其左半支和右半支的Pg和Pn震相交点的偏移距分别为40km和42km, 初步估算地壳厚度为7km (图9b), 与OBS36台站估算的地壳厚度大致相当。

图9 OBS36 (a)、OBS37 (b)、OBS09 (c)和OBS42 (d)的垂直分量综合地震记录剖面折合速度为6.0km·s-1, 折合时间=绝对走时−偏移距/6; Pw代表直达水波震相, Pg代表地壳内部折射震相, PmP代表莫霍面反射震相, Pn代表上地幔顶部折射震相Fig. 9 Comprehensive seismic record profiles of OBS36 (a), OBS37 (b), OBS09 (c) and OBS42 (d). Reduced velocity is 6.0 km·s-1, and reduced travel time = travel time − offset/6. Pw indicates direct water wave phases, Pg indicates crustal refractive phases, PmP indicates reflective phases from Moho, Pn indicates refractive phases from uppermost mantle

测线东侧OBS09和OBS42台站位于石申海山附近(图1), 石申海山由岩浆活动形成(任金锋 等,2018), 导致OBS09和OBS42右半支的震相都表现出弱震相的特征, 破碎模糊(图9c, 9d)。OBS09左半支的Pg震相在6～20km都为比较平坦的形态, 但在20km处突然向上隆起, 走时减少, 主要是其西侧一个基底隆起造成(图6)。OBS42右半支的PmP在17km就出现, Pg和Pn震相在偏移距22km处相交,由此估算台站之下的地壳厚度仅3～4km。

从多道地震剖面(图6)来看, 基底的双程走时在测线的西部和中部起伏较小而在东部的起伏相对较大, 东部地壳受到更强烈改造; 黄色部分多道地震剖面LW3莫霍面的双程走时与基底相比起伏相对较小; 测线西侧OBS36与OBS37台站下的基底和莫霍面之间的双程走时间隔为2～3s, 取地壳的平均速度为6km·s-1计算所得地壳厚度为6～9km, 与根据震相所估算的地壳厚度6～7km大致符合。

5 结论与展望

本文以OBS2018–H2测线数据为依托, 阐述了测线前期的数据处理流程, 进行了震相初步识别,并对地壳结构进行了初步分析, 得到以下几点认识:

1) 炮点与OBS位置校正效果良好。炮点位置校正消除了DGPS到枪阵中心的偏离; 直达水波形态由OBS位置校正前的左右不对称变为校正后的对称形, 校正后的OBS位置更接近于其实际落点;所有OBS的校正精度均在5ms以内, 以平均海水声速1.5km·s-1计算, 位置校正误差在8m以内, 可为后期速度结构模拟提供足够保障。

2) 基底在测线的西部和中部起伏较小而在东部的起伏相对较大; 由测线西侧OBS36、OBS37台站的震相特征初步估算台站下方地壳厚度约为6～7km, 这个厚度与由多道地震剖面LW3的双程走时估算的厚度6～9km大致相当。

3) 此次实验数据质量良好, 在10台OBS的综合地震记录剖面中均识别出多组清晰的P波走时震相, 如: Pw、Pg、PmP和Pn。下一步将以初始速度结构模型为基础, 综合测线上所有OBS台站的震相信息, 开展速度结构的正/反演模拟, 获取最终的速度结构模型, 可望揭示南海北部洋陆转换带的地壳结构特征与张裂-破裂过程。