CN116774281B - 一种地震面波与重力同步联合反演方法与系统 - Google Patents

Abstract

本发明提供了一种地震面波与重力同步联合反演方法与系统，包括：根据先验信息建立三维速度和密度初始模型；对三维速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；根据理论观测数据计算地震面波走时残差和重力异常残差；构建地震核函数和重力核函数，并基于残差数据建立目标函数；对目标函数进行求解得到最终的速度和密度模型。本发明通过重构地震核函数，使得关于密度的地震核函数得以重现，摆脱了传统地震面波与重力联合反演过程中对V_s速度‑密度关系式的依赖，同时引入重力核函数以提高密度的反演精度，最终实现了地震面波与重力的同步联合反演，为后续地球物理和地质解释提供可靠的速度和密度模型支撑。

Description

一种地震面波与重力同步联合反演方法与系统

技术领域

本发明涉及地震面波和重力数据处理技术领域，特别是涉及一种地震面波与重力同步联合反演方法与系统。

背景技术

单一地震面波层析成像方法和单一重力反演方法在许多领域已经有了成功应用。然而，由于受物性参数敏感度的限制，单一地球物理反演方法往往只能在模型空间囊括局部信息，容易造成单一反演方法缺乏约束，导致反演精度的降低。例如，重力数据对横向上的密度变化具有较高敏感度，但径向上的分辨率往往有所不足，而地震面波数据虽然有较高的径向分辨率，却很难为近距观测提供足够的速度信息。因此，将地震面波与重力数据联合反演，可以使得两种单一的地球物理反演方法在一定程度上形成优势互补，从而提高反演结果的可靠性。

当地震面波和重力联合反演时，需要将各单一地球物理数据所对应的较为敏感的物性参数之间进行联合，目前，国内外地震面波和重力联合反演主要采用的联合方式为物性耦合。基于该联合方式的常规地震面波与重力联合反演方法通常是将关于密度的重力核函数通过V_s速度与密度经验关系式转化为关于V_s速度的重力核函数，可在一定程度上提高V_s速度的反演精度。然而，这种利用V_s速度-密度关系式在各单一反演方法之间进行物性参数核函数转化的方法，其实质仍是对单一物性参数进行反演，既难以在缺乏V_s速度-密度关系式的区域进行推广，也无法完全摆脱联合反演过程中某单一物性参数起主导作用的现状，还制约了地震面波与重力同步联合反演的实现。

发明内容

为解决上述问题，本发明实施例的目的在于提供一种地震面波与重力同步联合反演方法与系统。

为实现上述目的，本发明提供了如下方案：

一种地震面波与重力同步联合反演方法，包括：

步骤1：获取反演区域的先验信息；

步骤2：根据所述先验信息建立三维V_s速度和密度初始模型；

步骤3：对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

步骤4：根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

步骤5：构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

步骤6：对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

优选的，所述步骤4：根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差，包括：

采用公式：

δt＝t^obs-t^cal

δg＝g_B+g₀-g_cal

计算地震面波走时残差δt和重力异常残差δg；其中，t^obs为地震面波频散观测走时，t^cal为地震面波理论走时，g_B为布格重力异常，g₀为常密度模型正演得到的重力背景场，g_cal为变密度模型正演得到的理论重力场。

优选的，在所述步骤5中，通过公式：

对地震核函数和重力核函数进行构建；其中，δt为地震面波走时残差，ω为均衡因子，为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>为/>核函数转换系数，/>为/> 核函数转换系数，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，δg为重力异常残差，/>为采用球坐标系二阶泰勒级数展开公式计算得到的重力核函数。

优选的，所述球坐标系二阶泰勒级数展开公式为：

其中，G为万有引力常数，i表示第i个观测点，j表示模型网格第j个Tesseroid单元体，N表示模型网格Tesseroid单元体总数，Δr、Δλ分别为各个Tesseroid单元体在径向、纬向、经向的大小，L₀₀₀为模型网格的零阶泰勒系数，L₂₀₀、L₀₂₀、L₀₀₂分别为模型网格在径向、纬向、经向上的二阶泰勒系数。

优选的，所述步骤5：构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数，包括：

采用公式：

构建地震面波与重力同步联合反演目标函数；其中，Ws为地震面波数据质量有关的权重系数，为重力相对权重系数；ω为均衡因子，/>为V_s速度的地震敏感核函数，为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>为密度的重力核函数，分别由地震面波频散对V_s、V_p、ρ的一阶偏导数得来；/>为/>核函数转换系数，/>为/>核函数转换系数，δt为地震面波走时残差，δg为重力异常残差，/>L_ρ分别为V_s速度和密度模型的平滑函数，I为单位矩阵，s为平滑系数，/>λ_ρ分别为V_s速度和密度模型的阻尼因子，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，L为矩阵的L范数。

本发明还提供了一种地震面波与重力同步联合反演系统，包括：

先验信息获取模块，用于获取反演区域的先验信息；

初始模型构建模块，用于根据所述先验信息建立三维V_s速度和密度初始模型；

正演模块，用于对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

残差计算模块，用于根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

目标函数构建模块，用于构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

目标函数求解模块，用于对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现上述的一种地震面波与重力同步联合反演方法中的步骤。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种基于地震核函数重构的面波与重力同步联合反演方法及系统，与现有技术相比，本发明通过重构地震核函数，使得关于密度的地震核函数得以重现，摆脱了传统地震面波与重力联合反演过程中对V_s速度-密度关系式的依赖，同时引入重力核函数以提高密度的反演精度，最终实现了地震面波与重力的同步联合反演，可以同步得到精度较高的速度和密度模型，为后续地球物理和地质解释提供可靠的V_s速度和密度模型支撑。

附图说明

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

图1是本发明提供的一种地震面波与重力同步联合反演方法流程图；

图2为依照本发明实施例的检测板测试结果中的V_s速度模型切片图，(a)-(d)为真实模型，(e)-(h)为单一面波反演结果，(i)-(l)为同步联合反演结果；

图3为依照本发明实施例的检测板测试结果中的密度模型切片图；(a)-(d)为真实模型，(e)-(h)为单一重力反演结果，(i)-(l)为单一面波反演结果，(m)-(p)为同步联合反演结果；

图4为依照本发明实施例的检测板测试结果中的图2V_s速度模型A-A'剖面对比图；(a)为真实模型，(b)为单一面波反演结果，(c)为同步联合反演结果；

图5为依照本发明实施例的检测板测试结果中的图3密度模型B-B'剖面对比图；(a)为真实模型，(b)为单一重力反演结果，(c)为单一面波反演结果，(d)为同步联合反演结果。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为实现上述目的，本发明提供了如下方案：

请参阅图1，一种地震面波与重力同步联合反演方法，包括：

步骤1：获取反演区域的先验信息；

步骤2：根据所述先验信息建立三维V_s速度和密度初始模型；

进一步地，在步骤2中，首先需要根据研究区地震面波频散的数据量选择合适的三维模型网格，为了便于计算，V_s速度和密度可选择相同的模型网格。根据地球物理反演理论，V_s速度和密度模型网格点为需要求取的模型扰动未知量，为了得到较好的反演效果，两套模型网格点总数应小于地震面波频散的数据量；其次，在确定V_s速度和密度模型网格后，结合已有的地质、地球物理等资料，选择合适的V_s速度初始模型；最后，根据速度-密度关系式(Brocher，2005)将V_s速度初始模型转化成密度模型，再结合已有的密度信息对其进行微调，得到了密度初始模型。

步骤3：对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

步骤4：根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

需要说明的是，本发明的步骤4，包括：

步骤4.1：采用公式：

δt＝t^obs-t^cal

计算地震面波走时残差δt；其中，t^obs为地震面波频散观测走时，t^cal为地震面波理论走时。t^obs由地震面波频散观测数据通过大圆路径公式计算得到。t^cal分为两步计算得到，首先，采用Dunkin矩阵计算法正演得到频散数据，该方法由Dunkin(1965)基于Thomson(1950)和Haskell(1953)发展而来；其次，利用快速行进算法(Fast Marching Method，简称FMM)计算得到面波理论走时。

步骤4.2：采用公式：

δg＝g_B+g₀-g_cal

计算重力异常残差；其中，g_B为布格重力异常，g₀为由常密度模型正演得到的重力背景场，g_cal为变密度模型正演得到的理论重力场。

本发明采用球坐标系二阶泰勒级数展开公式计算径向重力场：

其中，G为万有引力常数，ρ为Tesseroid单元体密度，Δr、Δλ分别为Tesseroid单元体在径向、纬向、经向的大小，L₀₀₀为零阶泰勒系数，L₂₀₀、L₀₂₀、L₀₀₂分别为径向、纬向、经向二阶泰勒系数。

步骤5：构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

在本发明的步骤5中，包括：

步骤5.1：重新构建地震核函数；

与传统地震面波与重力联合反演采用关于V_s速度的单一核函数不同的是，本发明采用公式：

对地震核函数进行了重构，使得关于密度的地震核函数得以重现；其中，ω为均衡因子，由公式求得；/>为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>分别由地震面波频散对V_s、V_p、ρ的一阶偏导数得来；/>为/>核函数转换系数，/>为/>核函数转换系数，分别由V_p-V_s、V_p-ρ关系式一阶倒数得来；δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动。

步骤5.2：构建重力核函数；

与传统联合反演中重力核函数计算不同的是，本发明采用球坐标系二阶泰勒级数展开公式：

计算关于密度的重力核函数；其中，G为万有引力常数，i表示第i个观测点，j表示模型网格第j个Tesseroid单元体，N表示模型网格Tesseroid单元体总数，Δr、Δλ分别为各Tesseroid单元体在径向、纬向、经向的大小，L₀₀₀为模型网格的零阶泰勒系数，L₂₀₀、L₀₂₀、L₀₀₂分别为模型网格在径向、纬向、经向上的二阶泰勒系数，δρ为待反演的密度模型扰动，δg为重力异常残差，/>为采用球坐标系二阶泰勒级数展开公式计算得到的重力核函数。

步骤5.3：构建目标函数；

与传统地震面波与重力联合反演不同的是，本发明建立了地震面波与重力同步联合反演架构，并采用公式：

构建地震面波与重力同步联合反演目标函数；其中，W^s为地震面波数据质量有关的权重系数，如果需要计算方便，可令W^s＝I，为重力相对权重系数；ω为均衡因子，由公式/>求得；/>为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，为密度的地震敏感核函数，/>为密度的重力核函数，/>分别由地震面波频散对V_s、V_p、ρ的一阶偏导数得来；/>为/>核函数转换系数，/>为/>核函数转换系数，/> 分别由V_p-V_s、V_p-ρ关系式一阶倒数得来；δt为地震面波走时残差，δg为重力异常残差，/>L_ρ分别为V_s速度和密度模型的平滑函数，I为单位矩阵，s为平滑系数，/>λ_ρ分别为V_s速度和密度模型的阻尼因子，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，L为矩阵的L范数。

步骤6：对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

需要说明的是，本发明的步骤6，包括：

步骤一：采用LSMR算法(Least Squares MR算法简称)求解V_s速度和密度模型扰动；

LSMR算法是在LSQR算法(Least Squares QR算法简称)基础上引入了迭代优化算法。

步骤二：选取参数s，λ_ρ，/>

采用公式：

对参数s，λ_ρ，/>进行选取，实现同步联合反演结果最优化；其中，s为平滑系数，/>λ_ρ分别为V_s速度、密度模型的阻尼因子，/>为重力相对权重；/>为归一化后的地震面波走时残差模量，R_s为面波走时残差，min(R_s)为R_s最小值，max(R_s)为R_s最大值；/>为归一化后的重力异常残差模量，R_g重力异常残差，min(R_g)为R_g最小值，max(R_g)为R_g最大值；/>为归一化后的模型模量；M_s为单一地震面波层析成像的模量和，M_u为同步联合反演的模量和，先使/>通过逐一调整s，/>λ_ρ(顺序可以不固定)使得M_s最小时，地震面波层析成像达到最优化，再调整相对重力权重/>使得M_u取最小值，则同步联合反演结果达到最优化。

当确定s，λ_ρ，/>后，地震面波与重力同步联合反演得到了V_s速度和密度模型扰动量，采用公式：

修正V_s速度和密度初始模型后，判断是否满足设定的迭代次数上限或者误差限制，若未满足条件，则将修正后的V_s速度和密度模型作为下次迭代的初始模型，返回步骤3直到满足迭代条件后输出最终V_s速度和密度模型。

下面本发明结合具体的实施例对本发明的地震面波与重力同步联合反演过程做进一步的说明：

图2和图3给出了青藏高原东北缘的检测板测试结果，图2为V_s速度结果，图3为密度结果。该模型试验中，V_s速度和密度的模型网格大小为0.6°×0.6°，初始V_s速度和密度扰动为±10％，检测板网格大小为1.8°×1.8°，深度范围为0-60km，检测板模型正演得到的面波频散和重力异常作为实测数据用于同步联合反演。

从测试结果可以看出，本发明的地震面波与重力同步联合反演方法可以同步获取分辨率较高的V_s速度和密度模型；由于面波频散数据对V_s速度较为敏感，在单一地震面波层析成像和同步联合反演结果中，V_s速度反演精度会略高于密度的反演精度；对比单一地震面波层析成像、单一重力反演和同步联合反演的密度结果可以发现，同步联合反演方法对各单一反演方法的密度结果分辨率具有明显改善效果；同步联合反演对单一重力反演的各深度密度分辨率都具有明显的改善效果；同步联合反演对单一地震面波层析成像的浅层密度分辨率的改善效果较深部更为明显，这也与重力反演的特征相吻合。

图4和图5分别为图2(a)所示的A-A′剖面和图3(a)所示的B-B′剖面。从图4和图5可以发现，在单一地震面波层析成像和同步联合反演结果中，V_s速度反演精度略高于密度的反演精度；同步联合反演对单一地震面波层析成像和单一重力反演的密度结果分辨率都具有明显的改善效果。

根据本发明提供的具体实施例，本发明公开了以下有益效果：

1)本发明实现了一种基于地震核函数重构的面波与重力同步联合反演方法。

2)本发明摆脱了传统地震面波与重力联合反演过程中对V_s速度-密度关系式的依赖。

3)本发明重构了地震核函数，并引入重力核函数，实现了地震面波与重力的同步联合反演。

4)本发明采用了迭代优化算法，实现了同步联合反演的结果最优化。

5)利用本发明的地震面波与重力同步联合反演方法，可以同步获取精度较高的三维V_s速度和密度模型。

本发明还提供了一种地震面波与重力同步联合反演系统，包括：

先验信息获取模块，用于获取反演区域的先验信息；

初始模型构建模块，用于根据所述先验信息建立三维V_s速度和密度初始模型；

正演模块，用于对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

残差计算模块，用于根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

目标函数构建模块，用于构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

目标函数求解模块，用于对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现上述的一种地震面波与重力同步联合反演方法中的步骤。

与现有技术相比，本发明通过对地震核函数的重构，并引入重力核函数，利用迭代优化算法，实现了地震面波走时数据与布格重力异常数据的同步联合反演，可以同步获取精度较高的三维V_s速度和密度模型，同步提高了单一地震面波层析成像和单一重力反演的精度，可为地质和地球物理解释提供高分辨率速度和密度模型。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上，本说明书内容不应理解为对本发明的限制。

Claims (5)

1.一种地震面波与重力同步联合反演方法，其特征在于，包括：

步骤1：获取反演区域的先验信息；

步骤2：根据所述先验信息建立三维V_s速度和密度初始模型；

步骤3：对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

步骤4：根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

所述步骤4：根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差，包括：

采用公式：

δt＝t^obs-t^cal

δg＝g_B+g₀-g_cal

计算地震面波走时残差δt和重力异常残差δg；其中，t^obs为地震面波频散观测走时，t^cal为地震面波理论走时，g_B为布格重力异常，g₀为常密度模型正演得到的重力背景场，g_cal为变密度模型正演得到的理论重力场；

步骤5：构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

在所述步骤5中，通过公式：

对地震核函数和重力核函数进行构建；其中，δt为地震面波走时残差，ω为均衡因子，为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>为/>核函数转换系数，/>为/>核函数转换系数，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，δg为重力异常残差，/>为采用球坐标系二阶泰勒级数展开公式计算得到的重力核函数；

步骤6：对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

2.根据权利要求1所述的一种地震面波与重力同步联合反演方法，其特征在于，所述球坐标系二阶泰勒级数展开公式为：

其中，G为万有引力常数，i表示第i个观测点，j表示模型网格第j个Tesseroid单元体，N表示模型网格Tesseroid单元体总数，Δr、Δλ分别为各个Tesseroid单元体在径向、纬向、经向的大小，L₀₀₀为模型网格的零阶泰勒系数，L₂₀₀、L₀₂₀、L₀₀₂分别为模型网格在径向、纬向、经向上的二阶泰勒系数。

3.根据权利要求2所述的一种地震面波与重力同步联合反演方法，其特征在于，所述步骤5：构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数，包括：

采用公式：

构建地震面波与重力同步联合反演目标函数；其中，W^s为地震面波数据质量有关的权重系数，为重力相对权重系数；ω为均衡因子，/>为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>为密度的重力核函数，分别由地震面波频散对V_s、V_p、ρ的一阶偏导数得来；/>为/>核函数转换系数，/>为/>核函数转换系数，δt为地震面波走时残差，δg为重力异常残差，L_ρ分别为V_s速度和密度模型的平滑函数，I为单位矩阵，s为平滑系数，/>λ_ρ分别为V_s速度和密度模型的阻尼因子，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，L为矩阵的L范数。

4.一种地震面波与重力同步联合反演系统，其特征在于，包括：

先验信息获取模块，用于获取反演区域的先验信息；

初始模型构建模块，用于根据所述先验信息建立三维V_s速度和密度初始模型；

正演模块，用于对所述三维V_s速度和密度初始模型进行正演，得到地震面波频散和重力异常的理论观测数据；

残差计算模块，用于根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差；

所述根据反演区域的地震面波频散数据、布格重力异常数据和所述理论观测数据计算地震面波走时残差和重力异常残差，包括：

采用公式：

δt＝t^obs-t^cal

δg＝g_B+g₀-g_cal

计算地震面波走时残差δt和重力异常残差δg；其中，t^obs为地震面波频散观测走时，t^cal为地震面波理论走时，g_B为布格重力异常，g₀为常密度模型正演得到的重力背景场，g_cal为变密度模型正演得到的理论重力场；

目标函数构建模块，用于构建地震核函数和重力核函数，并基于所述地震面波走时残差和所述重力异常残差建立地震面波与重力同步联合反演目标函数；

在所述目标函数构建模块中，通过公式：

对地震核函数和重力核函数进行构建；其中，δt为地震面波走时残差，ω为均衡因子，为V_s速度的地震敏感核函数，/>为V_p速度的地震敏感核函数，/>为密度的地震敏感核函数，/>为/>核函数转换系数，/>为/>核函数转换系数，δV_s为待反演的V_s速度模型扰动，δρ为待反演的密度模型扰动，δg为重力异常残差，/>为采用球坐标系二阶泰勒级数展开公式计算得到的重力核函数；

目标函数求解模块，用于对所述地震面波与重力同步联合反演目标函数进行迭代求解得到最终的V_s速度和密度模型。

5.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1-3中任一项所述的一种地震面波与重力同步联合反演方法中的步骤。