CN110221344B - 一种地壳三维密度结构的地震全波形与重力联合反演方法 - Google Patents

Abstract

本发明提供了一种地壳三维密度结构的地震全波形与重力联合反演方法，涉及区域地壳三维结构成像技术领域，能够充分发挥重力与地震各自优势，降低重力方法反演多解性，有效提高密度成像分辨率；该方法获取宽频地震波形数据和布格重力异常数据；构建地壳三维Vp、Vs速度和密度初始模型；利用布格重力异常数据和宽频地震波形数据进行地震全波形与重力联合反演，得到地壳三维密度精细结构模型。本发明提供的技术方案适用于地壳三维结构成像的过程中。

Description

一种地壳三维密度结构的地震全波形与重力联合反演方法

【技术领域】

本发明涉及区域地壳三维结构成像技术领域，尤其涉及一种地壳三维密度结构的地震全波形与重力联合反演方法。

【背景技术】

重力方法是深部结构与区域构造研究的重要方法之一，容易实现区域高分辨率覆盖，水平分辨率高。重力反演是获取壳幔密度结构的主要手段，但重力反演存在固有的多解性且深度分辨率低，需要引入先验信息约束以降低多解性、提高深度分辨率。地震层析成像是获取壳幔速度结构的主要手段，深度分辨率高。因此，发展地震与重力联合反演技术，实现地震与重力相互约束、优势互补，降低单一方法反演多解性，提高成像分辨率，成为壳幔结构反演的研究热点之一。

当前国际上地震与重力联合反演技术较多，反演壳幔结构的主要有体波走时与重力联合反演、面波与重力联合反演、体波走时、面波与重力联合反演等，反演形式有序贯反演和同步反演两类。然而，上述联合反演技术大都需要已知速度与密度的关联式作速度和密度的互换，而此关联式在很多地区往往未知且呈非均匀性，制约了其实际应用与推广。

全波形反演利用地震波场运动学和动力学信息重建地下复杂结构、揭示复杂地质构造细节，有效提高速度成像分辨率，为区域深部结构与构造精细成像等提供有力支撑，已成为当前地震层析成像发展的新趋势之一。天然地震全波形反演方法主要有伴随场法(Tromp et al.,2005)和散射积分法(Zhao et al.,2005)两种，后者计算需求相对小、效率高，适用于区域壳幔结构成像，但数据存储需求大。近年来，国内外学者通过全波形反演获得了美国南加州、西华盛顿州和中国东北的地壳三维速度精细结构，通过全波形反演获得了欧洲大陆、美国西北部和东亚的壳幔三维速度结构，这些成果分辨率比常规地震成像技术有明显提高，充分展现全波形反演的优越性和应用前景。

全波形反演虽然可以有效提高速度成像分辨率但并没有很好提高密度成像分辨率，主要由于地震波场对密度的敏感度低，需要引入重力约束改善密度成像。因此，发展全波形与重力联合反演技术，在全波形反演提高速度成像分辨率的同时有效提高密度成像分辨率，这具有重要的理论意义和应用价值。

因此，有必要设计一种新的地壳三维密度结构的地震全波形与重力联合反演方法，以解决或减轻上述一个或多个问题。

【发明内容】

有鉴于此，本发明提供了一种地壳三维密度结构的地震全波形与重力联合反演方法，能够利用研究区内地震台站宽频地震波形数据和区域高精度布格重力异常数据，通过对其进行全波形与重力联合反演，从而获得地壳三维密度精细结构，进而达到对研究区内的地壳构造及深部成矿机制、孕震机制等进行推断解释的目的。

一方面，本发明提供一种地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述方法的步骤包括：

S1、获取宽频地震波形数据和布格重力异常数据；

S2、筛选出震相清晰度和信噪比满足要求的宽频地震波形数据；

S3、根据先验信息构建地壳三维Vp、Vs速度和密度初始模型；

S4、利用布格重力异常数据和S2中筛选出的宽频地震波形数据进行地震全波形与重力联合反演，得到地壳三维密度精细结构模型。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述方法的步骤还包括：

S5、对S4中的联合反演的结果进行质量评价；

S6、根据所述地壳三维密度精细结构对地壳构造及深部成矿机制、孕震机制进行推断和解释。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述S3的具体步骤包括：

S31、构建起伏地表的地壳模型的三维网格架构；

S32、收集先验信息；所述先验信息包括地表地质、物性资料以及已探测得到的壳幔结构和界面结构；

S33、根据收集的先验信息构建地壳三维速度结构初始模型；

S34、根据波速比关系实现Vp、Vs速度结构互换，再根据速度和密度经验式，换算出地壳三维Vp、Vs速度和密度初始模型。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述S4的具体步骤包括：

S41、对S3中获得的地壳三维Vp、Vs速度和密度初始模型进行地震全波形模拟和重力模拟，得到理论地震波形数据、理论重力异常数据以及对应的地震波场三维核函数和重力场三维核函数；

S42、将实际地震波形数据、实际布格重力异常数据分别与理论地震波形数据和理论重力异常数据相减，得到相应的震相走时残差δd和重力异常残差δg；

S43、将计算得到的地震波场三维核函数、重力场三维核函数以及震相走时残差δd和重力异常残差δg，进行地震全波形与重力联合反演，得到地壳三维Vp、Vs速度和密度扰动模型；

S44、更新初始地壳三维模型，并重复S41～S43，直至震相走时残差δd和重力异常残差δg满足误差限值，得到地壳三维密度精细结构模型。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述地震全波形模拟和所述地震波场三维核函数采用空间四阶、时间二阶的交错网格有限差分算法进行模拟和计算。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述密度初始模型的重力模拟和重力场三维核函数采用直立长方体空间域解析式进行计算。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述重力模拟和重力场三维核函数的具体计算公式为：

其中，

G为万有引力常数；

x₀、y₀,、z₀为测点坐标；

x、y、z为直立长方体内任意一点的坐标；

(x₁,x₂)，(y₁,y₂)，(z₁,z₂)分别为直立长方体延x、y、z方向的坐标范围。

如上所述的方面和任一可能的实现方式，进一步提供一种实现方式，所述S43中地震全波形与重力联合反演的公式为：

sL_ρ·δρ＝0,λ_ρI·δρ＝0

其中：

δV_p、δV_s、δρ分别为Vp、Vs速度和密度模型的扰动量；

分别为Vp、Vs速度和密度的波场三维Fréchet核函数；

为密度的重力场三维核函数；

δd为地震数据函数定义的地震观测数据；

δg为实测布格重力异常Δg_B和模型模拟重力异常Δg_cal的残差，δg＝Δg_B-Δg_cal；

w^s为全波形拟合的权系数；

为重力异常拟合的权系数；

L_ρ分别为Vp、Vs速度和密度模型的光滑函数；

s为光滑系数；

I为单位矩阵；

λ_ρ分别为Vp、Vs速度和密度模型的正则化约束系数。

与现有技术相比，本发明可以获得包括以下技术效果：采用地震全波形与重力联合反演技术，充分发挥重力与地震各自优势，降低重力方法反演多解性，有效提高密度成像分辨率；反演过程不需要作速度和密度的互换，适用于众多地区的深部三维密度结构成像，更具实用性。

当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。

【附图说明】

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是本发明一个实施例提供的地壳三维密度结构的地震全波形与重力联合反演方法流程图；

图2为美国南加州研究区地形图、地震台站分布及2000～2018年部分历史地震震中分布(方框-研究区，三角形-固定地震台站，圆圈-近震震中)。

图3为美国南加州研究区实测布格重力异常(a)与联合反演模型理论重力异常(b)。

图4为美国南加州研究区地壳三维密度初始模型(a)与联合反演结果模型(b)。

图5为美国南加州研究区联合反演密度模型剖面图与构造解释：(a)地质简图与剖面位置(虚线)，(b)剖面A密度模型、岩性与断裂推断(虚线)及历史地震震中(黑点)叠合，(c)剖面B密度模型、岩性与断裂推断(虚线)及历史地震震中(黑点)叠合；EF-Elsinore断裂带，SJF-San Jacinto断裂带，SAF-San Andreas断裂带。

【具体实施方式】

为了更好的理解本发明的技术方案，下面结合附图对本发明实施例进行详细描述。

应当明确，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

在本发明实施例中使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本发明。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。

一种地壳三维密度结构的地震全波形与重力联合反演方法，所述方法首先收集研究区区域高精度布格重力异常数据和地震台站宽频地震波形数据；筛选出震相相对清晰、信噪比较高的宽频地震波形数据；结合研究区已知先验信息，构建地壳三维Vp、Vs速度和密度初始模型；对地震波形数据和布格重力异常数据进行地震全波形与重力联合反演，进而获取地壳三维密度结构模型；从而对研究区内的地壳构造或深部成矿机制、孕震机制等进行推断解释。

本发明的地壳三维密度结构的地震全波形与重力联合反演方法，如图1所示。所述方法包括如下步骤：

(1)收集研究区区域高精度布格重力异常数据和地震台站宽频地震波形数据；

(2)筛选出震相相对清晰、信噪比较高的宽频地震波形数据；

(3)研究区地壳三维Vp、Vs速度和密度初始模型取自Shaw et al(2015)的地壳Vp、Vs速度和密度模型；

(4)利用区域高精度布格重力异常数据和筛选后的波形数据进行地震全波形与重力联合反演，得到地壳三维密度精细结构模型；

(5)对上述联合反演结果进行质量评价；

(6)对研究区内的地壳构造或深部成矿机制、孕震机制等进行推断解释；例如，地壳内部花岗岩的展布、深大断裂展布、地震孕震机制等；

其中，步骤(3)中，构建研究区地壳三维Vp、Vs速度和密度初始模型的具体过程为：

(31)设计研究区地壳模型的三维网格框架，横向网格间距1km，深度范围0～48km；模型深度步长从近地表的0.2km递增到上地幔顶部的2km；采用基于曲形网格的三维模型剖分方法，构建起伏地表的地壳模型的三维网格架构；

(32)从Shaw et al(2015)的地壳Vp、Vs速度和密度模型插值获得研究区地壳Vp、Vs速度和密度初始模型。

步骤(4)中，地震全波形与重力联合反演的流程具体为：

(41)对研究区内地壳三维Vp、Vs速度和密度初始模型进行地震全波形模拟和重力模拟，得到理论地震波形数据和理论重力异常数据，同时得到对应的地震波场三维核函数

和重力场三维核函数

其中，地震全波形模拟和重力模拟方法具体为：

地震全波形模拟和三维核函数

计算采用空间四阶、时间二阶的交错网格有限差分算法精确计算，具体计算公式见文献(Zhao et al.,2005)；密度模型的重力模拟及重力三维核函数计算根据常规的直立长方体空间域解析式高精度计算而得，重力模拟和核函数计算公式分别如下：

其中，G为万有引力常数，G＝6.67×10^-11N·m²/kg²，σ为直立长方体剩余密度，单位为kg/m³，测点坐标为(x₀,y₀,z₀)，直立长方体内任意一点坐标为(x,y,z)，且有

直立长方体延x、y、z方向的坐标范围分别为(x₁,x₂)，(y₁,y₂)，(z₁,z₂)；

由于数值模拟、核函数计算和联合反演计算量大、数据存储需求大，采用集群并行快速计算和大容量数据存储设备高速存取；

(42)将实际地震波形数据(筛选后)、实际高精度布格重力异常数据分别与理论计算的地震波形数据和理论重力异常数据相减，得到相应的震相走时残差δd和重力异常残差δg；

(43)将计算得到的地震波场三维核函数

重力场三维核函数

以及震相走时残差δd和重力异常残差δg，进行地震全波形与重力联合反演，得到地壳三维Vp、Vs速度和密度扰动模型；

地震全波形与重力联合反演公式为：

sL_ρ·δρ＝0,λ_ρI·δρ＝0

其中，δV_p、δV_s、δρ分别为Vp、Vs速度和密度模型的扰动量；

分别为Vp、Vs速度和密度的波场三维Fréchet核函数(敏感度矩阵)，

为密度的重力场三维核函数(敏感度矩阵)；δd为广义地震数据函数(GSDF)定义的地震观测数据，它由实际波形数据和模型数值模拟数据震相走时残差描述；δg为实测布格重力异常Δg_B和模型模拟重力异常Δg_cal的残差，即δg＝Δg_B-Δg_cal；w^s是全波形拟合的权系数，

是重力异常拟合的权系数；

L_ρ为Vp、Vs速度和密度模型的光滑函数，s为光滑系数；I为单位矩阵，

λ_ρ为Vp、Vs速度和密度模型的正则化约束系数；

λ_ρ和s采用折中曲线(L-curve)评价机制确定；

模型正则化和光滑化约束的引入降低联合反演不确定性、提高收敛性和模型的地质意义，采用共轭梯度法或高斯-牛顿迭代法迭代求解上式联合反演方程，最终得到地壳三维密度扰动模型；

(44)更新初始地壳三维模型，并重复步骤(41)～(43)，直至震相走时残差δd和重力异常残差δg满足误差限值，最终得到地壳三维精细密度结构模型。

步骤(5)中，质量评价方法具体为：

质量评价采用多渠道评价机制，比如，利用研究区其它若干地震事件对结果模型作波场数值模拟，再与观测波形数据作对比和评价；分析结果模型的重力拟合残差或方差等。

在本实施例所述方法中，采用地震全波形与重力联合反演技术，使重力与地震相互约束、优势互补，反演不需要速度与密度互换，降低重力方法反演多解性，有效提高密度成像分辨率，有利于更精确地获得研究区地壳三维密度精细结构。

本发明以美国南加州研究区地壳三维密度结构重建为例，在美国南加州研究区用上述方法进行地壳三维密度结构重建：

步骤一：收集美国南加州研究区地震台站宽频地震波形数据和高精度布格重力异常数据；

步骤二：筛选出震相相对清晰、信噪比较高的波形数据；

步骤三：根据美国南加州研究区先验信息构建地壳三维Vp、Vs速度和密度初始模型；

步骤四：利用筛选后的波形数据和高精度布格重力异常数据进行地震全波形与重力联合反演，得到美国南加州研究区地壳三维密度精细模型；

步骤五：对上述联合反演结果进行质量评价；

步骤六：对研究区内的断裂、岩浆构造与孕震机制进行推断解释。

联合反演结果的密度成像分辨率得到显著提高，速度成像分辨率也得到进一步提高。依据联合反演结果取得了重要的新认识：(1)研究区地表大型花岗岩延深可达中地壳；(2)三条北西向深断裂带位于密度高低过渡带上且近陡直；(3)历史地震多发生在中、上地壳的密度高低过渡带和低密度带中。

本发明相较于现有技术的有益效果如下：

1)本发明采用地震全波形与重力联合反演技术，使重力与地震相互约束，充分发挥地震方法的高垂向分辨率、重力方法的高水平分辨率和区域高分辨率覆盖优势，降低重力方法反演多解性，有效提高密度成像分辨率；

2)本发明的全波形与重力联合反演技术不需求研究区的速度与密度关联式，反演过程不需要作速度和密度的互换，适用于众多地区的深部三维结构成像，更具实用性；

3)本发明采用前人的基于散射积分法的地震全波形反演技术(Zhao et al.,2005)，散射积分法虽然需求较大的数据存储能力，但计算需求小、效率高，适用于区域地壳三维结构成像，而且当今数据存储能力和集群并行计算技术较为发达，促使大数据的散射积分法全波形反演可行；

4)本发明采用前人的基于曲形网格的三维模型剖分方法(Zhang et al.,2012)构建起伏地表的三维反演模型，在此基础上获得的反演结果更贴近实际地质意义。

以上对本申请实施例所提供的一种地壳三维密度结构的地震全波形与重力联合反演方法，进行了详细介绍。以上实施例的说明只是用于帮助理解本申请的方法及其核心思想；同时，对于本领域的一般技术人员，依据本申请的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本申请的限制。

如在说明书及权利要求书当中使用了某些词汇来指称特定组件。本领域技术人员应可理解，硬件制造商可能会用不同名词来称呼同一个组件。本说明书及权利要求书并不以名称的差异来作为区分组件的方式，而是以组件在功能上的差异来作为区分的准则。如在通篇说明书及权利要求书当中所提及的“包含”、“包括”为一开放式用语，故应解释成“包含/包括但不限定于”。“大致”是指在可接收的误差范围内，本领域技术人员能够在一定误差范围内解决所述技术问题，基本达到所述技术效果。说明书后续描述为实施本申请的较佳实施方式，然所述描述乃以说明本申请的一般原则为目的，并非用以限定本申请的范围。本申请的保护范围当视所附权利要求书所界定者为准。

还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的商品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种商品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的商品或者系统中还存在另外的相同要素。

应当理解，本文中使用的术语“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。

上述说明示出并描述了本申请的若干优选实施例，但如前所述，应当理解本申请并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述申请构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本申请的精神和范围，则都应在本申请所附权利要求书的保护范围内。

Claims (6)

1.一种地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述方法的步骤包括：

S1、获取宽频地震波形数据和布格重力异常数据；

S2、筛选出震相清晰度和信噪比满足要求的宽频地震波形数据；

S3、根据先验信息构建地壳三维Vp、Vs速度和密度初始模型；

S4、利用布格重力异常数据和S2中筛选出的宽频地震波形数据进行地震全波形与重力联合反演，得到地壳三维密度精细结构模型；

所述S4的具体步骤包括：

S41、对S3中获得的地壳三维Vp、Vs速度和密度初始模型进行地震全波形模拟和重力模拟，得到理论地震波形数据、理论重力异常数据以及对应的地震波场三维核函数和重力场三维核函数；

S42、将实际地震波形数据、实际布格重力异常数据分别与理论地震波形数据和理论重力异常数据相减，得到相应的震相走时残差δd和重力异常残差δg；

S43、将计算得到的地震波场三维核函数、重力场三维核函数以及震相走时残差δd和重力异常残差δg，进行地震全波形与重力联合反演，得到地壳三维Vp、Vs速度和密度扰动模型；

S44、更新初始地壳三维模型，并重复S41～S43，直至震相走时残差δd和重力异常残差δg满足误差限值，得到地壳三维密度精细结构模型；

所述S43中地震全波形与重力联合反演的公式为：

sL_ρ·δρ＝0,λ_ρI·δρ＝0

其中：

δV_p、δV_s、δρ分别为Vp、Vs速度和密度模型的扰动量；

分别为Vp、Vs速度和密度的波场三维Fréchet核函数；

为密度的重力场三维核函数；

δd为地震数据函数定义的地震观测数据；

δg为实测布格重力异常△g_B和模型模拟重力异常△g_cal的残差，δg＝△g_B-△g_cal；

w^s为全波形拟合的权系数；

为重力异常拟合的权系数；

L_ρ分别为Vp、Vs速度和密度模型的光滑函数；

s为光滑系数；

I为单位矩阵；

λ_ρ分别为Vp、Vs速度和密度模型的正则化约束系数。

2.根据权利要求1所述的地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述方法的步骤还包括：

S5、对S4中的联合反演的结果进行质量评价；

S6、根据所述地壳三维密度精细结构对地壳构造及深部成矿机制、孕震机制进行推断和解释。

3.根据权利要求1所述的地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述S3的具体步骤包括：

S31、构建起伏地表的地壳模型的三维网格架构；

S32、收集先验信息；所述先验信息包括地表地质、物性资料以及已探测得到的壳幔结构和界面结构；

S33、根据收集的先验信息构建地壳三维速度结构初始模型；

S34、根据波速比关系实现Vp、Vs速度结构互换，再根据速度和密度经验式，换算出地壳三维Vp、Vs速度和密度初始模型。

4.根据权利要求1所述的地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述地震全波形模拟和所述地震波场三维核函数采用空间四阶、时间二阶的交错网格有限差分算法进行模拟和计算。

5.根据权利要求1所述的地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述密度初始模型的重力模拟和重力场三维核函数采用直立长方体空间域解析式进行计算。

6.根据权利要求5所述的地壳三维密度结构的地震全波形与重力联合反演方法，其特征在于，所述重力模拟和重力场三维核函数的具体计算公式为：

其中，

G为万有引力常数；

x₀、y₀,、z₀为测点坐标；

x、y、z为直立长方体内任意一点的坐标；

(x₁,x₂)，(y₁,y₂)，(z₁,z₂)分别为直立长方体延x、y、z方向的坐标范围。

Images