CN113671570B - 一种地震面波走时和重力异常联合反演方法与系统 - Google Patents

一种地震面波走时和重力异常联合反演方法与系统 Download PDF

Info

Publication number
CN113671570B
CN113671570B CN202110966090.XA CN202110966090A CN113671570B CN 113671570 B CN113671570 B CN 113671570B CN 202110966090 A CN202110966090 A CN 202110966090A CN 113671570 B CN113671570 B CN 113671570B
Authority
CN
China
Prior art keywords
inversion
data
model
wave
velocity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110966090.XA
Other languages
English (en)
Other versions
CN113671570A (zh
Inventor
李雪垒
许辉
王元
郝天珧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan University of Technology
Original Assignee
Hunan University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan University of Technology filed Critical Hunan University of Technology
Priority to CN202110966090.XA priority Critical patent/CN113671570B/zh
Publication of CN113671570A publication Critical patent/CN113671570A/zh
Application granted granted Critical
Publication of CN113671570B publication Critical patent/CN113671570B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02—Details
    • G01V7/06—Analysis or interpretation of gravimetric records

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种地震面波走时和重力异常联合反演方法与系统，本发明提供的一种地震面波走时和重力异常联合反演方法包括：首先对待反演的地震面波频散数据和布格重力异常数据进行正演得到理论观测数据并计算敏感核；然后根据理论观测数据和敏感核进行反演得到线性化反演模型；基于岩石的速度‑密度关系对线性化反演模型进行修正得到修正后的线性化反演模型；最后根据修正后的线性化反演模型建立联合反演目标函数；利用反演目标函数的解对待反演的地震面波进行反演成像。本发明通过利用地震面波频散数据和布格重力异常数据对地震面波进行联合反演成像，不仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性。

Description

一种地震面波走时和重力异常联合反演方法与系统
技术领域
本发明涉及地震数据处理技术领域，特别是涉及一种地震面波走时和重力异常联合反演方法与系统。
背景技术
面波层析成像从数据来源上，可以分为地震面波和背景噪声面波层析成像，这主要是依据面波是从一个真实的地震震源处激发，还是从一个虚拟的台站处激发来考虑的。对于后一种情况，其中的“面波”更准确来说是一种在特定条件下面波占据主要成分的叠加信号。在射线类方法中，又可以按照是否将反演群相速度图作为中间步骤分为两步法和一步法。同时，由于面波本质上是一个(粘)弹性介质的本征值问题，它在低频下也常常通过地球自由震荡的方式表现出来，因而可以用来研究地球的大尺度结构。尽管分类比较繁多，但这也从侧面说明了目前已经对面波开展了多手段的研究工作。
重力异常数据反映了地球真实的密度分布相对于参考椭球体的扰动。考虑到重力大小随距离平方衰减的性质，它通常被认为反映了地壳浅部和莫霍面附近的密度异常。因此，重力异常数据通常用来研究莫霍面起伏或者地壳内的密度变化。目前已经有很多文章利用重力异常来反演莫霍面结构，反演浅地表和地壳内密度结构等。重力反演研究的主要思路是首先利用特定的场分离方法从观测重力异常中提取需要研究的异常体对应的重力异常信号，之后再利用反演方法获取该异常体的密度结构。
由于地震面波具有频散特征，各个周期的面波传播速度不同。由于这种频散效应是由地下介质唯一决定的，因此利用频散信息可以约束地下速度结构信息。重力异常数据反映了地球真实的密度分布相对于参考椭球体的扰动，考虑到重力大小随距离平方衰减的性质，其反映了浅层地质结构的密度异常。由于每一种类型的数据仅仅对某些区域或者某些尺度的物理量敏感，不能完全地刻画地下介质的精细结构，而联合反演可以通过取长补短的方式弥补这一缺陷。通过利用岩石物理实验中的速度-密度经验关系，联合地震面波频散和重力异常数据可以改善单一数据速度结构成像的分辨率，提高反演结果的可靠性。
发明内容
本发明的目的是提供一种地震面波走时和重力异常联合反演方法与系统，以解决单一数据集对地下结构约束能力较低的问题。
为实现上述目的，本发明提供了如下方案：
一种地震面波走时和重力异常联合反演方法，包括：
获取待反演的地震面波频散数据和布格重力异常数据；
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
计算敏感核；
根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
根据修正后的线性化反演模型建立联合反演目标函数；
对所述反演目标函数进行求解得到最优解；
利用所述最优解对待反演的地震面波进行反演成像。
优选的，所述对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据，包括：
采用模型：
F(m)＝F(α,β,ρ)＝d
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；其中，m表示真实模型，d表示理论观测数据，F表示表示真实模型到观测数据的映射，α表示P波速度，β表示S波速度，ρ表示介质密度。
优选的，所述线性化反演模型，包括：
Figure BDA0003223988470000031
其中，ω表示面波频率，
Figure BDA0003223988470000032
表示在面波频率为ω下的第i个观测面波走时，ti(ω)表示第i个理论面波走时，δ0ti(ω)表示第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δcj(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，Kij(ω)是相应的敏感核。
优选的，所述基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型，包括：
获取岩石的速度-密度关系；
将根据所述速度-密度关系将S波速度结构转换成密度和P波速度公式；所述密度和P波速度公式为：
α＝Fα(β),ρ＝Fρ(β)
其中，α表示P波速度，β表示S波速度，ρ表示密度，Fα表示S波速度模型到P波速度的映射，Fρ表示S波速度到密度的映射；
根据密度和P波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型。
优选的，所述修正后的线性化反演模型为：
Figure BDA0003223988470000033
d/dβ是对S波速度的全偏导数，ti(ω)表示第i个理论面波走时，δti(ω)表示修正后的第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的S波速度差值。
优选的，所述根据修正后的线性化反演模型建立联合反演目标函数，包括：
采用公式：
Figure BDA0003223988470000041
建立联合反演目标函数；其中，Gs表示面波数据理论算子，Gg表示重力数据理论算子，σ是二阶吉洪诺夫正则化系数，S是二阶导数矩阵，γ是阻尼系数，I是单位矩阵，Ws是面波走时权重矩阵，Wg是重力数据的权重矩阵，
Figure BDA0003223988470000042
Figure BDA0003223988470000043
其中p是一个0-1之间的常数，Ns是面波的数据个数，Ng是重力的数据个数，
Figure BDA0003223988470000044
表示面波数据中理论模型与真实模型之间的标准差，
Figure BDA0003223988470000045
表示重力数据中理论模型与真实模型之间的标准差。
本发明还提供了一种地震面波走时和重力异常联合反演系统，包括：
数据集获取模块，用于获取待反演的地震面波频散数据和布格重力异常数据；
正演模块，用于对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
计算模块，用于计算敏感核；
反演模块，用于根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
线性化反演模型修正模块，用于基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
联合反演目标函数构建模块，用于根据修正后的线性化反演模型建立联合反演目标函数；
求解模块，用于对所述反演目标函数进行求解得到最优解；
反演成像模块，用于利用所述最优解对待反演的地震面波进行反演成像。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供了一种地震面波走时和重力异常联合反演方法与系统，本发明提供的一种地震面波走时和重力异常联合反演方法包括：首先对待反演的地震面波频散数据和布格重力异常数据进行正演得到理论观测数据并计算敏感核；然后根据理论观测数据和敏感核进行反演得到线性化反演模型；基于岩石的速度-密度关系对线性化反演模型进行修正得到修正后的线性化反演模型；最后根据修正后的线性化反演模型建立联合反演目标函数；对反演目标函数进行求解得到最优解；利用最优解对待反演的地震面波进行反演成像。本发明通过利用地震面波频散数据和布格重力异常数据对地震面波进行联合反演成像，不仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是依照本发明实施例的速度模型的3次均匀B样条插值结果；
图2是依照本发明实施例的地震面波走时和重力异常联合反演方法流程图；
图3是依照本发明实施例的检测板测试分析图，(a-d)为真实模型，(e-h)为面波单独反演结果，(i-l)为联合反演结果；
图4是依照本发明实施例的异常体测试分析图，(a-b)为真实模型，(c-d)为面波单独反演结果，(e-f)为联合反演结果；
图5是依照本发明实施例的4条剖面上的速度反演结果，第一列为真实模型，第二列为面波单独反演结果，第三列为联合反演结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种地震面波走时和重力异常联合反演方法与系统，以解决单一数据集对地下结构约束能力较低的问题。
为实现上述目的，本发明提供了如下方案：
一种地震面波走时和重力异常联合反演方法，包括：
步骤1:获取待反演的地震面波频散数据和布格重力异常数据；
步骤2:对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
其中，步骤2具体包括：
采用了高斯-勒让德数值积分方法计算重力场的径向分量，然后基于径向分量采用模型：
F(m)＝F(α,β,ρ)＝d
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；其中，m表示真实模型，d表示理论观测数据，F表示表示真实模型到观测数据的映射，α表示P波速度，β表示S波速度，ρ表示介质密度。
步骤3:计算敏感核；
步骤4:根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；在本发明中，所述线性化反演模型，包括：
Figure BDA0003223988470000061
其中，ω表示面波频率，
Figure BDA0003223988470000062
表示在面波频率为ω下的第i个观测面波走时，ti(ω)表示第i个理论面波走时，δ0ti(ω)表示第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δcj(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，Kij(ω)是相应的敏感核。
步骤5:基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
其中，步骤5具体包括：
获取岩石的速度-密度关系；
将根据所述速度-密度关系将S波速度结构转换成密度和P波速度公式；所述密度和P波速度公式为：
α＝Fα(β),ρ＝Fρ(β)
其中，α表示P波速度，β表示S波速度，ρ表示密度，Fα表示S波速度模型到P波速度的映射，Fρ表示S波速度到密度的映射；
根据密度和P波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型；所述修正后的线性化反演模型为：
Figure BDA0003223988470000071
d/dβ是对S波速度的全偏导数，ti(ω)表示第i个理论面波走时，δti(ω)表示修正后的第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的S波速度差值。
步骤6:根据修正后的线性化反演模型建立联合反演目标函数；
其中，步骤6具体包括：
采用公式：
Figure BDA0003223988470000072
建立联合反演目标函数；其中，σ是二阶吉洪诺夫正则化系数，S是二阶导数矩阵，γ是阻尼系数，I是单位矩阵，Ws是面波走时权重矩阵，Wg是重力数据的权重矩阵，
Figure BDA0003223988470000073
Figure BDA0003223988470000074
其中p是一个0-1之间的常数，Ns是面波的数据个数，Ng是重力的数据个数，
Figure BDA0003223988470000075
表示面波数据中理论模型与真实模型之间的标准差，
Figure BDA0003223988470000076
表示重力数据中理论模型与真实模型之间的标准差。
步骤7:对所述反演目标函数进行求解得到最优解；
步骤8:利用所述最优解对待反演的地震面波进行反演成像。
下面结合具体的实施例对本发明提供的一种地震面波走时和重力异常联合反演方法进行进一步的说明：
每一种类型的数据仅仅对某些区域或者某些尺度的物理量敏感，它们不能完全地刻画地下介质的精细结构。而本发明中的联合反演可以通过取长补短的方式弥补这一缺陷。因此，本发明需要用两套数据来构建一个联合反演问题，包括：模型参数化：通过一系列参数作为模型空间的“基矢量”，使得模型空间内每一点都能够用这些参数表示(即给定插值方式),同时给定这些参数的初始值。一般来说，表征一个模型通常包括模型类的参数化方法和节点类的参数化方法。他们的区别是在前者中这些参数表征的是模型基函数的系数，而在后者中表示的是网格节点上的模型物理量的具体数值。当然在某些情况下两种表述是完全一致的。正演计算：通过这些参数计算出该模型对应的理论观测数据，在本发明的方法中，在正演建模时考虑了面波射线的几何形态以及地球的曲率对面波和重力数据的影响。
请一并参阅图1-5，为使本发明的目的和技术方案更加的便于理解，下面将参照附图更为详细的说明本发明的实施例。
本发明采用了高斯-勒让德数值积分方法计算重力场的径向分量：
Figure BDA0003223988470000081
其中，(r,θ,φ)表示球坐标系的三个方向，对应径向分量、水平分量和纵向分量，i,j,k表示每个方向对应的积分阶数，gr是在接收点(r0,θ0,φ0)处的重力场的径向分量，ωijk表示权重系数，NI,J,K表示使用的高斯-勒让德积分的最高阶数。
由于本发明在面波走时正演中使用了快速多极算法方法，该方法要求利用球面坐标系的均匀网格节点来进行有限差分型的计算，因此需要对这种均匀分布的节点网格选择一个比较好的插值方法。对于均匀分布的网格，一个比较好的参数化方法是利用3次均匀B样条来对速度模型插值。图1给出了经三次均匀B样条插值后的一个速度模型的分布特征，可以看到有间断的模型在参数化后具有了很好的平滑特征。
同时，为了更准确地计算面波走时，在模型参数化过程中有一些小技巧。第一，由于实际问题使用的初始模型的尺度可能不足以对速度模型进行足够精确的采样，这通常会导致走时场计算不准确。因此，在做正演计算时，首先需要构建一套计算网格。这套网格是用模型网格在更密的网格节点上插值得到的。
另一个问题是震源附近的速度结构对走时的影响非常大，而且震源位置并不一定就位于网格节点上。这样会造成初始几个临近节点中的走时计算不准确，从而影响最终的结果。因此先要利用模型网格将震源附近的网格局部加密，然后将这个局部网格上的计算出来的走时插值到计算网格上。最终，为了执行一次正演计算，本发明在模型网格节点之外，还需要另外构建两套用于计算的网格。
对于面波一维模型的正演，本发明仍然采用节点式的模型，但是在深度方向让网格的间距可以自由变化以适应速度模型的跳变。在局部化假设中，为了利用汤姆森-哈斯克尔矩阵方法计算频散曲线，需要将网格节点模型转换为层状模型：这里本发明假设垂直方向的速度线性变化，在两个网格节点之间插入几个均匀层，层间的弹性参数用其中点位置的值来代替。在这样的假设下，就可以利用三维节点模型，首先计算出频散图，之后用快速多极算法方法计算各个周期的面波走时。
同时，为了计算的方便，本发明让重力正演计算的网格与面波正演的网格位于同一套网格节点上。最终就可以利用这套三维网格节点上的模型参数来合成理论观测数据,正演问题也可以进一步表示为：
F(m)＝F(α,β,ρ)＝d#(2.2)
其中，m表示真实模型，d表示理论观测数据,F表示真实模型到观测数据的映射，α,β,ρ表示真实模型的三个模型参数，分别表示P波速度，S波速度，介质密度。
从之前的论述可以看出，面波的正演分为两个步骤。因此，面波走时的反演问题也可以相应地分成两步:首先要利用面波走时反演每个周期的相速度、群速度频散图，之后再用一系列的一维反演来获得地下三维结构。
第一步线性化反演可以写成：
Figure BDA0003223988470000101
其中，ω表示面波频率，
Figure BDA0003223988470000102
是该周期(频率)下的第i个观测面波走时，ti(ω)是该周期下的第i个理论面波走时，δ0ti(ω)表示第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δcj(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，Kij(ω)是敏感核函数(模型参数的广义导数)。
由于面波频散对S波的速度结构更敏感，可以用岩石物理中的经验关系将S波速度结构转换成密度和P波速度公式：
α＝Fα(β),ρ＝Fρ(β)#(2.4)
其中，α表示P波速度，β表示S波速度，ρ表示密度，Fα表示S波速度模型到P波速度的映射，Fρ表示S波速度到密度的映射。
传统的两步法利用(2.3)式来反演频散图，之后再用(2.4)式来反演S波速度。但第一步反演可能会将误差传递到第二步反演中。为了克服上述问题，可以将(2.3)-(2.4)式结合起来，得到：
Figure BDA0003223988470000103
其中，d/dβ是对S波速度的全偏导数，ti(ω)是该周期下的第i个理论面波走时，δti(ω)表示第i个观测面波走时和理论面波走时之间的残差，cj(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的S波速度差值。
需要说明的是，在当前的局部化假设下，每一个走时只依赖于两个台站的之间的面波射线路径之下的弹性参数。因此，面波走时敏感核也是稀疏矩阵。利用(2.4)式，可以把密度的扰动量转换为对S波速度的扰动量：
Figure BDA0003223988470000111
其中，β表示S波速度，ρ表示密度，d/dβ是对S波速度的全偏导数，δβ表示真实模型和理论模型之间的S波速度差值。
之后利用(2.5)式，可以将每一次迭代过程中面波和重力数据组装到一个稀疏线性系统中：
Figure BDA0003223988470000112
其中，Gs表示面波数据理论算子，Gg表示重力数据理论算子，δβi表示真实模型和理论模型之间的S波速度差值，δt表示面波走时残差，δg表示重力异常残差。
上式的解可以看作一个优化问题：
Figure BDA0003223988470000113
其中，σ和S分别是二阶吉洪诺夫正则化系数和二阶导数矩阵，γ是阻尼系数，I是单位矩阵，Ws和Wg分别是面波走时和重力数据的(对角)权重矩阵：
Figure BDA0003223988470000114
Figure BDA0003223988470000115
其中，p是一个0-1之间的常数，需要根据实际问题的特点来调整，Ns，g分别是面波和重力的数据个数，σs，g是两套数据集中理论模型与真实模型之间的标准差。(2.8)式可以在2范数下通过最小二乘算法高效地求解。之后，可以通过：
βi+1＝βi+δβi#(2.10)
来更新模型，直到误差降到合理的范围内。
图3示出了中国四川盆地地区检测板测试分析图，利用初始模型加上±10％的速度扰动构成的，在水平方向上的尺度大约是1.5°×1.5°，在垂直方向上，考虑到两种数据敏感核的有效范围，异常体仅放置在60km以内，利用检测板模型正演得到所有台站对走时和所有观测点的布格重力异常后，在正演的面波数据中加入标准差1.4s的高斯噪声，重力异常数据中加入10mGal的高斯噪声作为观测数据，之后利用和实际数据同样的方法进行反演。
从该结果可以看出可以看出，在60km之内的各个深度范围内，面波单独反演总体上恢复了速度结构，但是在射线路径之外的分辨率有所下降。而联合反演在射线路径内部得到了和单独反演类似的结果，同时在浅层射线路径之外的模型分辨率也有了显著提高。
为了研究联合反演方法相比单独反演在减少反演多解性上的优势，本发明进行了一个异常体合成测试。本发明在一个一维模型上的不同深度的位置加入了4个异常值为±0.4km/s的矩形异常体(图5)。之后用了和检测板测试同样的方法来恢复模型，反演结果如图4-5所示。
图4和图5给出了反演结果在两个深度上的水平剖面和在4条垂直剖面上的反演结果。可以看出，面波单独反演恢复了模型的主要特征，但同时也引入了很多假异常。这些异常可能是由于对带噪声的走时数据过拟合引起的。而联合反演的结果明显更加“干净”，更接近真实模型。
本发明实施例提供了一种地震面波走时和重力异常联合反演方法。其中，方法包括：获取待反演的地震面波频散数据和布格重力异常数据；基于程函方程的直接面波层析成像和球坐标系下的自适应高斯-勒让德数值积分来进行两套数据集相应正演计算和反演敏感核；利用岩石物理实验中的速度-密度经验关系，将面波直接反演方法以及球坐标系下的高斯-勒让德重力积分方法相结合；对观测到的面波走时数据和布格重力异常直接建模进行反演成像。在本发明的方法中，在正演建模时考虑了面波射线的几何形态以及地球的曲率对面波和重力数据的影响。理论合成测试表明，和面波单独反演相比，尤其是在重力有较强约束能力的区域内，本发明能获得更合理的结果。最后，将该方法在中国四川盆地区域进行测试分析，实施例验证了新方法反演地下速度结构的可靠性。从反演误差上看，联合反演的速度模型相对于单独反演大大减小了重力数据的残差，说明引入重力数据可以提高联合反演可靠性，反演更为准确的地下速度结构信息，可为深部矿产资源勘查和地震灾害预警提供先验的速度模型。
本发明还提供了一种地震面波走时和重力异常联合反演系统，包括：
数据集获取模块，用于获取待反演的地震面波频散数据和布格重力异常数据；
正演模块，用于对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
计算模块，用于计算敏感核；
反演模块，用于根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
线性化反演模型修正模块，用于基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
联合反演目标函数构建模块，用于根据修正后的线性化反演模型建立联合反演目标函数；
求解模块，用于对所述反演目标函数进行求解得到最优解；
反演成像模块，用于利用所述最优解对待反演的地震面波进行反演成像。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供了一种地震面波走时和重力异常联合反演方法与系统，本发明提供的一种地震面波走时和重力异常联合反演方法包括：首先对待反演的地震面波频散数据和布格重力异常数据进行正演得到理论观测数据并计算敏感核；然后根据理论观测数据和敏感核进行反演得到线性化反演模型；基于岩石的速度-密度关系对线性化反演模型进行修正得到修正后的线性化反演模型；最后根据修正后的线性化反演模型建立联合反演目标函数；对反演目标函数进行求解得到最优解；利用最优解对待反演的地震面波进行反演成像。本发明通过利用地震面波频散数据和布格重力异常数据对地震面波进行联合反演成像，不仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上，本说明书内容不应理解为对本发明的限制。

Claims (4)

1.一种地震面波走时和重力异常联合反演方法，其特征在于，包括：
获取待反演的地震面波频散数据和布格重力异常数据；
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
计算敏感核；
根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
所述基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型，包括：
获取岩石的速度-密度关系；
将根据所述速度-密度关系将S波速度结构转换成密度和P波速度公式；所述密度和P波速度公式为：
Figure DEST_PATH_IMAGE002
其中，
Figure DEST_PATH_IMAGE004
表示P波速度，
Figure DEST_PATH_IMAGE006
表示S波速度，
Figure DEST_PATH_IMAGE008
表示介质密度，
Figure DEST_PATH_IMAGE010
表示S波速度模型到P波速度的映射，
Figure DEST_PATH_IMAGE012
表示S波速度到密度的映射；
根据密度和P波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型；
所述修正后的线性化反演模型为：
Figure DEST_PATH_IMAGE014
Figure DEST_PATH_IMAGE016
是对S波速度的全偏导数，
Figure DEST_PATH_IMAGE018
表示第
Figure DEST_PATH_IMAGE020
个理论面波走时，
Figure DEST_PATH_IMAGE022
表示修正后的第
Figure 789135DEST_PATH_IMAGE020
个观测面波走时和理论面波走时之间的残差，
Figure DEST_PATH_IMAGE024
是在自由表面第
Figure DEST_PATH_IMAGE026
个网格节点处的频散速度值，
Figure DEST_PATH_IMAGE028
表示真实模型和理论模型之间的S波速度差值；
所述根据修正后的线性化反演模型建立联合反演目标函数，包括：
采用公式：
Figure DEST_PATH_IMAGE030
建立联合反演目标函数；其中，
Figure DEST_PATH_IMAGE032
表示面波数据理论算子，
Figure DEST_PATH_IMAGE034
表示重力数据理论算子，
Figure DEST_PATH_IMAGE036
是二阶吉洪诺夫正则化系数，
Figure DEST_PATH_IMAGE038
是二阶导数矩阵，
Figure DEST_PATH_IMAGE040
是阻尼系数，
Figure DEST_PATH_IMAGE042
是单位矩阵，
Figure DEST_PATH_IMAGE044
是面波走时权重矩阵，
Figure DEST_PATH_IMAGE046
是重力数据的权重矩阵，
Figure DEST_PATH_IMAGE048
，
Figure DEST_PATH_IMAGE050
，其中
Figure DEST_PATH_IMAGE052
是一个0-1之间的常数，
Figure DEST_PATH_IMAGE054
是面波的数据个数，
Figure DEST_PATH_IMAGE056
是重力的数据个数，
Figure DEST_PATH_IMAGE058
表示面波数据中理论模型与真实模型之间的标准差，
Figure DEST_PATH_IMAGE060
表示重力数据中理论模型与真实模型之间的标准差，
Figure DEST_PATH_IMAGE062
表示面波走时残差，
Figure DEST_PATH_IMAGE064
表示重力异常残差；
根据修正后的线性化反演模型建立联合反演目标函数；
对所述反演目标函数进行求解得到最优解；
利用所述最优解对待反演的地震面波进行反演成像。
2.根据权利要求1所述的一种地震面波走时和重力异常联合反演方法，其特征在于，所述对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据，包括：
采用模型：
Figure DEST_PATH_IMAGE066
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；其中，
Figure DEST_PATH_IMAGE068
表示真实模型，
Figure DEST_PATH_IMAGE070
表示理论观测数据，
Figure DEST_PATH_IMAGE072
表示真实模型到观测数据的映射，
Figure DEST_PATH_IMAGE074
表示P波速度，
Figure DEST_PATH_IMAGE076
表示S波速度，
Figure DEST_PATH_IMAGE078
表示介质密度。
3.根据权利要求2所述的一种地震面波走时和重力异常联合反演方法，其特征在于，所述线性化反演模型，包括：
Figure DEST_PATH_IMAGE080
其中，
Figure DEST_PATH_IMAGE082
表示面波频率，
Figure DEST_PATH_IMAGE084
表示在面波频率为
Figure 614265DEST_PATH_IMAGE082
下的第
Figure 987478DEST_PATH_IMAGE020
个观测面波走时，
Figure 212048DEST_PATH_IMAGE018
表示第
Figure 495262DEST_PATH_IMAGE020
个理论面波走时，
Figure DEST_PATH_IMAGE086
表示第
Figure 973516DEST_PATH_IMAGE020
个观测面波走时和理论面波走时之间的残差，
Figure 443812DEST_PATH_IMAGE024
是在自由表面第
Figure 236187DEST_PATH_IMAGE026
个网格节点处的频散速度值，
Figure DEST_PATH_IMAGE088
表示自由表面第
Figure 952338DEST_PATH_IMAGE026
个网格节点处观测频散速度值和理论频散速度值之间的差值，
Figure DEST_PATH_IMAGE090
是相应的敏感核。
4.一种地震面波走时和重力异常联合反演系统，其特征在于，包括：
数据集获取模块，用于获取待反演的地震面波频散数据和布格重力异常数据；
正演模块，用于对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
计算模块，用于计算敏感核；
反演模块，用于根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
线性化反演模型修正模块，用于基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；所述基于岩石的速度-密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型，包括：
获取岩石的速度-密度关系；
将根据所述速度-密度关系将S波速度结构转换成密度和P波速度公式；所述密度和P波速度公式为：
Figure DEST_PATH_IMAGE002A
其中，
Figure 663811DEST_PATH_IMAGE004
表示P波速度，
Figure 247501DEST_PATH_IMAGE006
表示S波速度，
Figure 453355DEST_PATH_IMAGE008
表示介质密度，
Figure 304636DEST_PATH_IMAGE010
表示S波速度模型到P波速度的映射，
Figure 62376DEST_PATH_IMAGE012
表示S波速度到密度的映射；
根据密度和P波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型；
所述修正后的线性化反演模型为：
Figure DEST_PATH_IMAGE014A
Figure 927171DEST_PATH_IMAGE016
是对S波速度的全偏导数，
Figure 671136DEST_PATH_IMAGE018
表示第
Figure 642503DEST_PATH_IMAGE020
个理论面波走时，
Figure 305565DEST_PATH_IMAGE022
表示修正后的第
Figure 503329DEST_PATH_IMAGE020
个观测面波走时和理论面波走时之间的残差，
Figure 910039DEST_PATH_IMAGE024
是在自由表面第
Figure 237378DEST_PATH_IMAGE026
个网格节点处的频散速度值，
Figure 212287DEST_PATH_IMAGE028
表示真实模型和理论模型之间的S波速度差值；
所述根据修正后的线性化反演模型建立联合反演目标函数，包括：
采用公式：
Figure DEST_PATH_IMAGE030A
建立联合反演目标函数；其中，
Figure 553138DEST_PATH_IMAGE032
表示面波数据理论算子，
Figure 763540DEST_PATH_IMAGE034
表示重力数据理论算子，
Figure 948314DEST_PATH_IMAGE036
是二阶吉洪诺夫正则化系数，
Figure 94125DEST_PATH_IMAGE038
是二阶导数矩阵，
Figure 125535DEST_PATH_IMAGE040
是阻尼系数，
Figure 874048DEST_PATH_IMAGE042
是单位矩阵，
Figure 143355DEST_PATH_IMAGE044
是面波走时权重矩阵，
Figure 555007DEST_PATH_IMAGE046
是重力数据的权重矩阵，
Figure 73713DEST_PATH_IMAGE048
，
Figure 501284DEST_PATH_IMAGE050
，其中
Figure 890677DEST_PATH_IMAGE052
是一个0-1之间的常数，
Figure 768503DEST_PATH_IMAGE054
是面波的数据个数，
Figure 384292DEST_PATH_IMAGE056
是重力的数据个数，
Figure 238722DEST_PATH_IMAGE058
表示面波数据中理论模型与真实模型之间的标准差，
Figure 748201DEST_PATH_IMAGE060
表示重力数据中理论模型与真实模型之间的标准差，
Figure 406715DEST_PATH_IMAGE062
表示面波走时残差，
Figure 634434DEST_PATH_IMAGE064
表示重力异常残差；
联合反演目标函数构建模块，用于根据修正后的线性化反演模型建立联合反演目标函数；
求解模块，用于对所述反演目标函数进行求解得到最优解；
反演成像模块，用于利用所述最优解对待反演的地震面波进行反演成像。
CN202110966090.XA 2021-08-23 2021-08-23 一种地震面波走时和重力异常联合反演方法与系统 Active CN113671570B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110966090.XA CN113671570B (zh) 2021-08-23 2021-08-23 一种地震面波走时和重力异常联合反演方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110966090.XA CN113671570B (zh) 2021-08-23 2021-08-23 一种地震面波走时和重力异常联合反演方法与系统

Publications (2)

Publication Number Publication Date
CN113671570A CN113671570A (zh) 2021-11-19
CN113671570B true CN113671570B (zh) 2022-04-19

Family

ID=78544899

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110966090.XA Active CN113671570B (zh) 2021-08-23 2021-08-23 一种地震面波走时和重力异常联合反演方法与系统

Country Status (1)

Country Link
CN (1) CN113671570B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114721044B (zh) * 2022-04-21 2023-03-10 湖南工商大学 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统
CN116774281B (zh) * 2023-06-29 2024-01-30 中国地质大学（北京） 一种地震面波与重力同步联合反演方法与系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113189641A (zh) * 2021-03-25 2021-07-30 西安石油大学 一种两道多模式瑞利波地下探测系统及方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013012470A1 (en) * 2011-07-21 2013-01-24 Exxonmobil Upstream Research Company Adaptive weighting of geophysical data types in joint inversion
WO2013052035A1 (en) * 2011-10-04 2013-04-11 Westerngeco, L.L.C. Methods and systems for multiple-domain inversion of collected data
CA2904008C (en) * 2013-03-15 2020-10-27 Schlumberger Canada Limited Methods of characterizing earth formations using physiochemical model
CN104237970B (zh) * 2014-09-23 2017-07-07 中国石油天然气集团公司 地震电磁联合勘探系统及其数据采集装置和数据采集方法
CN105549106B (zh) * 2016-01-07 2018-06-08 中国科学院地质与地球物理研究所 一种重力多界面反演方法
CN110221344B (zh) * 2019-06-17 2020-08-28 中国地质大学（北京） 一种地壳三维密度结构的地震全波形与重力联合反演方法
CN111045076A (zh) * 2019-12-10 2020-04-21 核工业北京地质研究院 多模式瑞雷波频散曲线并行联合反演方法
CN111366987A (zh) * 2020-04-21 2020-07-03 中油奥博(成都)科技有限公司 地面地震微重力联合测量系统及数据采集处理方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113189641A (zh) * 2021-03-25 2021-07-30 西安石油大学 一种两道多模式瑞利波地下探测系统及方法

Also Published As

Publication number Publication date
CN113671570A (zh) 2021-11-19

Similar Documents

Publication Publication Date Title
Hole Nonlinear high‐resolution three‐dimensional seismic travel time tomography
EP2567063B1 (en) Artifact reduction in method of iterative inversion of geophysical data
CN113671570B (zh) 一种地震面波走时和重力异常联合反演方法与系统
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
Lelièvre et al. Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the fast marching method
MX2011003850A (es) Estimado de señal de dominio de imagen a interferencia.
CN115508908A (zh) 一种地震面波走时和重力异常联合反演方法与系统
CN109636912A (zh) 应用于三维声呐图像重构的四面体剖分有限元插值方法
Hu et al. Ray-illumination compensation for adjoint-state first-arrival traveltime tomography
Guo et al. Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique
CN113109883B (zh) 基于等参变换全球离散网格球坐标下卫星重力场正演方法
Liang et al. Scattering pattern analysis and true‐amplitude generalized Radon transform migration for acoustic transversely isotropic media with a vertical axis of symmetry
CN112596106B (zh) 一种球坐标系下重震联合反演密度界面分布的方法
CN111665556A (zh) 地层声波传播速度模型构建方法
CN110824558B (zh) 一种地震波数值模拟方法
CN111736208A (zh) 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质
CN111158059A (zh) 基于三次b样条函数的重力反演方法
CN116088048A (zh) 包含不确定性分析的各向异性介质多参数全波形反演方法
CN116466402A (zh) 一种基于地质信息和电磁数据联合驱动的电磁反演方法
Wang et al. Rock fracture monitoring based on high-precision microseismic event location using 3D multiscale waveform inversion
CN115373022A (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
CN113267830A (zh) 基于非结构网格下二维重力梯度与地震数据联合反演方法
CN113608262A (zh) 利用平动分量计算旋转分量的地震数据处理方法及装置
Li et al. A dual-layer equivalent-source method for deriving gravity field vector and gravity tensor components from observed gravity data
CN116660974B (zh) 一种基于结构耦合约束的体波和面波三维联合反演方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant