主动源海底地震仪横波模型的自动搜索与非唯一性分析*

张浩宇, 丘学林, 黄海波

1. 中国科学院边缘海与大洋地质重点实验室, 中国科学院南海海洋研究所, 广东 广州 510301;

2. 中国科学院南海生态环境工程创新研究院, 广东 广州 510301;

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

4. 中国科学院大学, 北京 100049;

5. 中国石油化工股份有限公司石油物探技术研究院, 江苏 南京 211103

主动源广角反射、折射地震探测是对沉积层、地壳和上地幔顶部成像的有力方法, 这一方法促进了构造演化的研究。主动源广角反射、折射地震探测很好地填补了反射地震勘探和被动源地震研究之间的分辨率、有效探测深度的空白, 能够以相对较高的分辨率来研究地壳尺度的构造问题。在海底地震仪的帮助下, 主动源广角反射、折射地震探测拓展进入海洋地球物理研究, 并且在海底构造研究中发挥重要作用。与海上传统的多道反射地震探测相比, 主动源OBS(Ocean Bottom Seismometer)探测系统由于接收系统和震源系统的分离, 能够接收到大角度的地震波, 特别是超临界角的折射波, 它们直接携带了海底构造的速度信息(Rodríguez et al,2008)。这些大角度、远偏移距(可达200km)的地震波在地壳深部甚至上地幔顶部传播了较远的距离,这对于深部结构的成像提供了宝贵信息。此外, OBS探测系统能够记录到来自海底不同波阻抗界面的转换横波(S 波), 这也是传统的海上多道地震探测所难以实现的, 这一特点使得主动源地震海底S 波结构研究成为可能。

横波速度(Vs)结合纵波速度(Vp), 可以获取海底结构的纵横波速比(Vp/Vs)。波速比是反映地下结构岩性和组成成分的良好指标(Watanabe, 1993)。室内测试的结果表明, 波速比对于压强和温度的变化不敏感(Christensen, 1996; Kern et al, 1999; Chevrot et al, 2000), 而对于石英和长石的含量呈现较高的依赖(Christensen et al, 1975; Kern, 1982), 同时, 流体的存在也对波速比产生影响(Nur et al, 1969; Spencer et al, 1976)。因此, 相比于单独的Vp,Vp/Vs能够提供对矿物成分更为精确的约束, 这一特性有助于更好地进行地壳属性研究、共轭陆缘对比研究, 以及确认上地幔蛇纹石化和地下结构中流体的存在, 甚至有助于预测埋藏于海底沉积物中的天然气水合物(张洁 等, 2018)。

基于射线追踪的震相走时正演模拟是广角反射、折射地震数据的主要分析方式, 尤其是对于转换S 波数据分析(Mjelde, 1992; Kodaira et al, 1996;Digranes et al, 2003; Zhao et al, 2010; Wei et al, 2015;Hou et al, 2019)。这类正演模拟以试错法的形式, 通过对初始速度模型的不断调整, 来取得理论走时与观测走时的良好拟合。这一方式能够很好地帮助识别和充分利用震相, 并且在模拟的过程中, 模型每次改动都是明确的。但是, 这一模拟方式也存在着一些问题, 比如工作量较大, 难以找到最优解, 无法对模型的非唯一性进行估计等问题。本文将对这些问题进行分析并尝试提出相应的解决方案。

1 主动源OBS 转换S 波模拟方法的现状

1.1 转换S 波模拟方法

主动源OBS 广角反射、折射地震数据的应用目前主要集中于震相走时部分, 而这一部分主要是通过走时的拟合来进行。对于这类地震走时数据, 要获取速度结构, 一般有两种途径。一种是地震走时的反演, 即通过对观测到的地震相走时来求解速度模型。这是一个典型的反问题, 目前求解此类问题有多款程序可供选择, 比如TOMO2D(Korenaga et al, 2000)、FAST(Zelt et al, 1998)、JIVE3D(Hobro et al,2003)等。以TOMO2D 程序为例, 速度结构的求取归结为如下线性方程组的求解:

式中,d为地震走时残差,G为Frechet 偏导矩阵,m为参数化的速度模型,m1为一次迭代更新后的速度模型, Δm为未知的速度模型扰动向量。通过迭代地求解式(1)的Δm, 并将速度模型m更新, 当m的地震走时射线追踪正演结果的偏差满足预设标准(阈值)时, 反演停止, 当前的速度模型则为最终的反演结果。在实际使用反演程序前, 往往需进行必要的参数测试, 主要涉及到阻尼、平滑因子等。在程序反演过程中, 不再需要人为干预。

另一种速度结构模拟是通过试错法(Trial-anderror)正演的形式获取, 在初始模型的基础上, 根据模型的射线追踪、走时拟合情况, 遵循“由浅至深”的原则, 对速度模型进行手动调整。迭代地进行上述过程, 就可以将初始模型逐步调整为满足走时拟合要求的最终正演模型。

对于转换S 波的速度结构模拟而言, 目前主要采用的是试错法正演方式(Zhao et al, 2010; Wei et al,2015, 2017; Hou et al, 2019)。一般来说, 转换S 波的波型转换界面往往是不唯一的, 现有的反演程序难以处理多个转换界面同时存在的问题, 而试错法正演能够较好地识别和利用多种来自不同地下界面的转换震相。以RAYINVR 中常用的正演模块为例, 其结合P 波速度结构的参数化方式(层-块形式)(Zelt et al, 1992), 以P 波速度结构为基础, 通过赋予速度模型中每个块体的泊松比(或波速比)来定义S 波结构(图1a), 同时速度层之间的界面被作为潜在的波型转换界面, 可实现对传播路径的分段(图1b)。转换界面的运动学可行性可通过射线追踪程序来验证,其动力学可行性可通过求解Zoepprritz 方程来定量化分析(张浩宇 等, 2021)。在试错法调整过程中, 可以根据需要对块体中的泊松比进行调整, 必要的时候还可以调整转换界面的设置。在实际的OBS 转换S 波数据模拟过程中, 由于初始模型与最终模型的差距往往较大, 而且走时拟合要兼顾多台、深浅震相, 手动调整是一个费时、重复的过程。然而, 这些繁琐的手动调整并不一定能够确保获得速度结构的最优解。

图1 转换S 波模拟速度模型参数化方式和射线路径示意图a. 速度模型参数化方式, ri 代表该块体单元内的波速比, Vp 为纵波速度, Vs 为横波速度, x 为水平坐标, z 为垂直坐标, Vp1～Vp4 分别是块体4 个顶点处的纵波速度; b. 转换S 波射线路径, 其中黑色粗实线为P 波路径, 黑色粗虚线为S 波路径, 红色虚线圈定了图a 所示块体的位置, 蓝色实线为海面, 棕色实线为海底面, 黑色细实线为地层界面Fig. 1 Diagram for velocity model parameterization and ray paths

1.2 转换S 波模拟结果的评价

模型的评价方式分为间接评价和直接评价(Zelt,1999)。间接评价包含模型分辨率和不确定度的衡量,而直接评价的方法尝试获取所有能够以可接受的拟合程度满足观测的备选模型。间接方式通常较为简单和快捷, 而直接方式往往较为费时。但是, 间接方式无法像直接方式那样能够充分地解决模型的非唯一性问题。直接法能够帮助我们更好地了解模型的非唯一性边界, 并且建立相应模型参数的权衡关系(trade-off)。

目前, 对于转换S 波模拟结果的评价主要是根据模型的走时残差(root mean square, rms)、χ2和射线覆盖情况(Zhao et al, 2010; Wei et al, 2015; Hou et al, 2019)。rms 是速度模型对应的理论走时和观测走时之间的均方根误差,χ2是以震相拾取不确定度归一化的走时残差。这两个值均可用于衡量理论走时和观测走时的拟合程度,χ2越接近于1, 说明模型能更好地从走时拟合的角度反映地下结构。射线覆盖是根据射线追踪的结果对模型每个网格内射线经过的次数进行统计得到的。射线覆盖次数与追踪到的理论震相数呈正相关, 射线覆盖次数越多,说明模型受震相约束得越完备。射线覆盖次数在模型内分布的均匀程度也会影响模型的可靠性。通过展示具体的射线路径, 并将理论和观测走时并列展示来直接对比两者的偏差, 可以反映速度模型重现观测走时的能力。上述方式均属于间接评价,能够提供最终模型可靠性的信息, 但是无法对模型的非唯一性进行估计。也就是说, 这些方式都无法展示观测数据所要求的模型结构或者说哪些范围的模型能够满足观测数据。尽管最终选定的模型可以提供对于观测数据的较好拟合, 但是模型空间的欠采样和观测误差的存在会导致解的非唯一性问题, 因此需要探究可行的方法来对模型的非唯一性进行分析。

2 主动源OBS 转换S 波模拟新方法

2.1 基于模型解空间和目标函数的转换S 波模拟技术

针对上述问题, 本文提出新的方法来实现主动源转换S 波数据模拟。转换S 波模拟得到的S 波速度结构建立在P 波速度结构的基础之上, P 波速度结构已经约束了速度和深度节点的分布, 因而速度界面和层-块分布也可以视为固定的。试错法正演过程中的模型调整, 实质上是对波速比或泊松比模型的调整。在震相对应的转换模型确定的情况下, 只需要调整速度模型内块体单元的波速比或泊松比, 这为构建S 波速度结构的解空间提供了较大便利。模型解空间定义了模型搜索的范围, 该空间包含了一系列具有相同层-块设置的波速比模型, 层-块单元中的波速比按照一定的间隔分布。理论上, 波速比可以指定任意正值, 但是这会导致无限大的解空间。根据研究区的地质条件限制, 模型的解空间得以缩小为有限大小。假定模型中第k层的波速比值(Vp/Vs)介于rkmin及rkmax之间, 并规定取值步长, 将k层Vp/Vs的集合记为Rk, 而所有层(共n层)Rk的集合则定义为模型的解空间。这里的“集合”相对于数学定义的集合, 有两点不同: 一是可以有重复元素;二是元素存在排列顺序。模型解空间的建立为有序搜索、遍历模型创造了条件。对于取值步长而言, 由于泊松比随着波速比的变化并不是均匀的, 在不同的波速比范围内, 泊松比的变化“快慢”不同(图2)。在常见的海底沉积层波速比范围(>1.80)内, 泊松比相对于波速比的变化率(一阶导数)较低, 变化较慢,尤其是在波速比>3.00 的范围内; 而在常见的地壳内波速比范围(1.70～1.80)内, 泊松比相对于波速比的变化率较大, 变化较快, 此时波速比取值步长应尽量缩小, 以实现对泊松比的充分采样。

从规模庞大的解空间内挑选出合适的模型是模型评价的主要目的。模型的质量可定义为其重现观测走时的准确程度, 模型质量往往受数据质量和模型参数化方式的影响(Loureiro et al, 2016)。在给定的速度模型参数化形式下, 可能存在一系列能够较好地拟合观测走时的模型。模型评价的任务就是对模型追踪震相和拟合走时的能力进行衡量, 并将能够尽可能多地追踪震相同时保证较低走时残差的模型挑选出来。如何设计目标函数是一个关键问题。在设计目标函数时, 需要对多个参数进行组合, 达到同时兼顾走时拟合残差和震相追踪数的效果。

走时拟合的效果可以通过观测走时和理论走时之间的误差来表示, 常用的衡量指标有均方根误差(rms)及其归一化χ2。走时残差的χ2是根据观测走时不确定度归一化得到的, 这一指标使得不同模型的拟合效果对比成为可能(Zelt, 1999)。每组震相的走时残差与观测走时的不确定度相当时,χ2等于1。追踪震相的能力可以用模型能够追踪到的震相数占拾取震相数的比例来衡量。值得注意的是, 速度模型的局部不均匀结构, 可能会削弱其追踪震相的能力。理想的速度模型(或波速比模型)应该可以追踪到大部分的观测走时, 并且以接近1 的χ2实现观测走时和理论走时的拟合。用于大量模型评价的目标函数在不同的模型和观测之间应当是可以比较的,因此追踪的震相数占识别的震相数的比例和归一化的走时残差χ2都应被纳入目标函数中。基于以上原则, 通过将这两个参数加权组合, 本文设计了如下目标函数来进行射线追踪和走时拟合效果的衡量:

式中,Nr为射线追踪到的震相数;Np为观测(拾取)的震相数;w用于调节走时偏差在目标函数中的贡献, 取值范围是0～1, 其值越大表明走时残差越受重视。该目标函数具备下列特性: 1) 值域为[0, 1];2) 随着Nr的增大和χ2向1 靠近而单调递增; 3) 当Nr等于Np并且χ2等于1 时, 目标函数f达到最大值; 4) 对于拟合不足或过拟合的情况, 比如χ2等于2 或者0.5, 目标函数都给出相同的评价。当w取值不同时, 目标函数随着Nr/Np和χ2的变化情况不同, 但整体上都满足上述特征(图3)。使用此目标函数对解空间内的模型进行评价, 当目标函数取得极大值时, 其对应的波速比模型即被认为是最佳模型, 最佳模型能够较好地兼顾震相拟合数量和走时残差。

图3 不同走时残差贡献率(w)下的目标函数分布图Nr 为射线追踪到的震相数, Np 为观测的震相数, χ2 为归一化的走时残差Fig. 3 Objective functions under different contributions (w) from travel-time misfit

通过模型解空间的构建和目标函数的设计, 模型自动搜索和走时拟合效果的比较得以完备而准确地定义。为了将这一方法用程序来实现, 本研究编写了相应的软件, 用于后续的数据分析。

2.2 模型非唯一性分析的实现

速度模型的分析在速度结构模拟中是至关重要的, 其影响着后续地质解释的有效性。速度模型的可靠性和非唯一性是速度模型分析的重要方面。通过观察射线追踪、观测走时和理论走时的对比以及射线覆盖次数的分析, 可以对模型的可靠性进行判断, 但是这些分析方法难以估计模型的非唯一性。而对于模型非唯一性的估计又是建立在获取大量可满足观测的速度模型基础上, 采用试错法手动调整的速度结构模拟方法会耗费巨大的工作量, 且难以保证搜寻到所有符合条件的模型。

本文提出的新模拟方法恰恰可以解决这一问题,实现对模型非唯一性的估计。通过模型解空间的构建, 可以确保包含所有符合区域地质背景的速度模型。通过对解空间内模型的射线追踪, 能够获取各个模型的走时拟合情况, 基于拟合结果进而获取对于模型非唯一性的估计, 即能够以可接受的误差拟合观测走时的模型分布情况。与此同时, 还可以获取到走时拟合情况随波速比的变化情况, 从而有助于反映该模型参数的分辨能力。

3 应用实例

南海北部陆缘位于华南大陆和南海海盆之间,经历了完整的张裂-海底扩张过程, 并且遭受的后期构造运动破坏较少, 因此是研究张裂-破裂过程的理想区域(Nissen et al, 1995)。西沙地块位于南海北部陆缘西段, 是南海陆缘的重要构造单元, 也是研究南海北部陆缘构造演化的重要切入点。西沙地块伴随着南海陆缘的伸展从华南陆缘裂离至现今位置(Qiu et al, 2001), 因此西沙地块与陆缘伸展的作用关系是研究南海陆缘演化的重要参考。西沙地块的深部地壳结构记录了上述信息, 采集于西沙地块的OBS2013-3 深地震测线为揭示西沙地块的深部地壳结构提供了重要资料(图4)。该测线投放了15 台OBS,投放间隔为12km 或18km, 全部成功回收。测线上的OBS 均装配了四分量检波器, 四分量包含一个垂直分量, 两个正交的水平分量以及一个水听器分量。基于该测线的OBS 数据, 本文进行了转换S 波模拟新技术的验证工作。通过前期的处理, 我们获取了OBS 径向和切向分量的地震记录(张浩宇 等,2021)。根据发生波型转换时地震波的传播方向, 转换S 波震相被分为两类, 一类是在向下传播路径上发生的P-S 转换(PSS 型), 另一类是在向上传播路径上发生的P-S 转换(PPS 型)(Kodaira et al, 1996)。该测线共有13 个台站在径向分量地震记录上识别出了有效的转换S 波震相, 其中包含3263 个PPS 型震相和2297 个PSS 型震相。转换S 波模拟的初始模型主要参考了P 波模拟的结果(郭晓然 等, 2016),同时也参考了临近测线的P 波速度模型(Huang et al,2019, 2021)。在P 波速度模型的基础上, 赋予模型中海水层、沉积层1、沉积层2、上地壳、下地壳、上地幔的泊松比分别为0.500、0.437、0.437、0.257、0.257 和0.276。将泊松比转换为波速比, 结合原有的P 波速度, 可进一步定义初始的S 波速度结构。初始的S 波速度模型仅有“层”的泊松比设置, 未包含“块”的泊松比设置。

图4 西沙海域海底地形与OBS2013-3 测线位置示意图图中黑色直线为放炮测线; 红色圆点为OBS 位置Fig. 4 Map showing the location of OBS2013-3 and submarine topography in the Xisha waters

3.1 单台站沉积层波速比结构模拟

张浩宇等(2021)在OBS2013-3 测线上识别并拾取了PPgSb、PPmSb 等PPS 型震相, 以及PmS、SmS 和PSgSw 等PSS 型震相。转换S 波模拟中, 沉积层的S 波结构主要由PPS 型震相约束。以OBS11台站为例, 本研究在该台站的径向分量上拾取了一系列PPS 型震相, 其中包含来自地壳内部的折射震相和来自Moho 面的反射震相。以1.8～4.8 为波速比, 以0.1 为波速比步长, 建立了双层的沉积层波速比结构解空间, 解空间内包含961 个模型,使用自编软件对该解空间进行搜索, 结果如图 5所示。目标函数指示最佳模型沉积层1 的波速比为4.0, 沉积层2 的波速比为2.4, 该模型的射线追踪和走时拟合情况见图6。模拟结果反映了随着深度增加, 压实作用增强, 波速比降低, 且两层沉积层波速比的较大差异指示了沉积间断的存在(Ma et al, 2011)。

图5 沉积层波速比解空间内的模型对OBS11 台站PPS型震相的走时拟合情况图中等值线为走时残差(χ2), 阴影深浅代表了实际追踪震相的丢失率大小, 空白表示该区域对应的模型未能正常完成射线追踪Fig. 5 Models in the solution space of Vp/Vs in the sedimentary layers and their travel-time fitting of PPS phases in OBS11

图6 OBS11 台站最佳沉积层波速比模型对PPS 型震相的射线追踪和走时拟合情况a. 震相射线追踪情况, 虚线段为S 波路径, 实线段为P 波路径; b. 走时拟合情况, 彩色为观测震相, 黑色为理论折合走时, 折合速度为6.0km.s–1Fig. 6 The ray tracing and travel-time fitting of PPS phases from OBS11 with the optimum Vp/Vs model of sedimentary layers

3.2 地壳波速比结构模拟

地壳波速比结构是进行转换S 波速度结构模拟的重要目标。利用地壳内以S 波形式传播的PSS 型震相(PmS、SmS 和PSgSw), 可对地壳内部的S 波速度结构进行研究。拾取测线上所有台站的PSS 型震相, 并且建立地壳内部的波速比模型解空间。结合前人研究和前期测试, 本文将模型中地壳波速比的取值范围限定在1.70～1.80 之间。因为在这一区间内,泊松比随着波速比的变化较快, 相应的岩性指示也变化较快, 故有必要在这一区间采用更为精细的波速比步长。以1.70～1.80 的波速比范围和0.005的波速比步长, 构建上、下地壳的波速比解空间,这一解空间包含了441 个模型。同样使用自编程序对解空间内的模型进行搜索, 得到了如图 7 所示的射线追踪和走时拟合统计结果。目标函数指示最佳模型的上地壳波速比为1.715, 下地壳波速比为1.735, 该模型的射线追踪和走时拟合效果如图8 所示。

图7 地壳波速比解空间内的模型对OBS2013-3 测线上全部台站PSS 型震相的走时拟合情况图中等值线为走时残差(χ2), 阴影深浅代表了实际追踪震相的丢失率大小, 空白表示该区域对应的模型未能正常完成射线追踪Fig. 7 Models in the solution space of crustal Vp/Vs and their travel-time fitting of PSS phases of all stations

图8 最佳地壳波速比模型对PSS 型震相的射线追踪和走时拟合情况。a. 震相射线追踪情况, 虚线段为S 波路径, 实线段为P 波路径; b. 走时拟合情况, 彩色线为观测震相, 黑色线为理论走时, 折合速度为3.5km.s–1Fig. 8 The ray tracing and travel-time fitting of PSS phases from all stations with the optimum crustal Vp/Vs model

3.3 模型非唯一性分析

由于观测走时数据的数量是有限的, 同时观测走时本身也存在误差, 而对真实的、连续的地下结构进行全面描述所需要的参数是无限的, 故地震成像结果是多解的, 因此需要充分评估解的非唯一性,并分析其非唯一性范围。在前述的转换横波速度结构模拟中, 我们获取了设定的模型解空间内全体模型的射线追踪和走时拟合情况。从走时拟合的角度,可以获取一系列满足条件的模型, 这些模型的边界就可以作为模型的非唯一性范围。

利用单个台站的PPS 型震相约束OBS11 台站下方的沉积层波速比结构, 若以χ2=1.2 作为较好的拟合程度, 合理的沉积层速度比结构分布范围将以约0.1 的沉积层2 波速比宽度呈条带状展布于整个模型空间; 若以χ2=1.0 作为较好的拟合程度, 则观测走时所要求的波速比模型收敛至最佳模型一处。这表明OBS11 台站的PPS 震相对沉积层2 波速比的分辨能力更强。

利用全部台站的PSS 型震相约束测线下方地壳的波速比结构, 以χ2=1.6 为依据, 本研究获取了图7中长条状的波速比范围; 以χ2=1.4 为依据, 模型的非唯一性范围缩小为3 个独立的区域, 目标函数揭示的最佳模型落在其中最大的独立区域内。

单台地震仪的PPS 震相, 其射线路径较为单一、集中, 主要分布在台站正下方的沉积层内, 上、下两层沉积层波速比的trade-off 效应较为显著, 当χ2要求不严格时, 存在贯穿模型解空间的非唯一性范围。而PSS 震相的射线路径广泛分布在地壳内,trade-off 效应不明显, 有利于圈定较小的非唯一性范围。这说明射线路径的多样化能够降低模型的非唯一性。

4 结论

通过对主动源OBS 转换S 波模拟技术所存在的问题进行分析, 本文提出了转换S 波模拟的新方法,并在实际的OBS 测线上进行了验证。

1) 通过模型解空间的构建和目标函数的设计,建立了模型自动搜索和评价的基础。这一方法包含了对大量模型的系统性搜索, 并以目标函数对模型的射线追踪和走时拟合能力进行比较。结合编程实现, 该方法可以从预设的解空间中找到最佳模型,达到自动模拟的效果, 从而减少了研究者的工作量,提升了数据分析的自动化程度。

2) 该方法能够得到以一定程度拟合观测走时的模型集合, 实现了对模型的直接评价, 为模型的非唯一性分析提供了可行方案。

3) 将该方法应用于OBS2013-3 测线上, 分别测试了单个台站的PPS 震相和全部台站的PSS 震相,获取了最佳波速比结构, 并对相应模型的非唯一性进行了分析。两种情况下, 模型非唯一性范围的差异都指示了射线路径的多样性能够降低模型的非唯一性。

本文研究提升了主动源OBS 转换S 波数据模拟的效率和可靠性, 有助于该类研究在更多测线上的推广使用, 获取更为丰富的地震学和岩石学信息,为研究海底结构和物质组成提供有力支撑。