基于球面波的叠前密度反演方法及装置

摘要

本发明实施例提供了一种基于球面波的叠前密度反演方法及装置，其中，该方法包括：获取入射角为单角度的地震资料，地震资料包括观测地震资料和理论地震资料，地震资料包括低频段的频点；计算观测地震资料的球面波反射系数和理论地震资料的球面波反射系数，低频段的频点的球面波反射系数随角频率变化而变化；采用非常快速模拟退火算法，对观测地震资料的球面波反射系数和理论地震资料的球面波反射系数进行频率域的叠前密度反演。该方案利用低频频变信息弥补了传统叠前密度反演信息量不足的缺点，并利用球面波反射系数的频变效应是弹性参数的客观反映这一事实，使得基于低频频变信息的叠前密度反演只需要入射角为单角度的地震资料即可实现。

权利要求

1.一种基于球面波的叠前密度反演方法，其特征在于，包括：

获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；

计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；

采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演；

通过以下公式计算所述理论地震资料中各频点的球面波反射系数：

其中， 是所述理论地震资料中各频点的球面波反射系数，

是平面波反射系数，平面波对应的入射角是θ，

x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。

2.如权利要求1所述的基于球面波的叠前密度反演方法，其特征在于，所述地震资料是入射角小于30度的地震资料，所述地震资料是界面深度小于500米的地震资料。

3.如权利要求1或2所述的基于球面波的叠前密度反演方法，其特征在于，采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演，包括：采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数同时进行密度比和速度比反演。

4.如权利要求1或2所述的基于球面波的叠前密度反演方法，其特征在于，在采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演之前，还包括：通过以下目标函数公式确定与密度同时进行反演的弹性参数：

其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

5.一种基于球面波的叠前密度反演装置，其特征在于，包括：

获取模块，用于获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；

计算模块，用于计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；

反演模块，用于采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演；

所述计算模块具体用于，通过以下公式计算所述理论地震资料中各频点的球面波反射系数：其中， 是所述理论地震资料中各频点的球面波反射系数，

是平面波反射系数，平面波对应的入射角是θ，

x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。

6.如权利要求5所述的基于球面波的叠前密度反演装置，其特征在于，所述地震资料是入射角小于30度的地震资料，所述地震资料是界面深度小于500米的地震资料。

7.如权利要求5或6所述的基于球面波的叠前密度反演装置，其特征在于，所述反演模块具体用于，采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数同时进行密度比和速度比反演。

8.如权利要求5或6所述的基于球面波的叠前密度反演装置，其特征在于，还包括：反演确定模块，用于在采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演之前，通过以下目标函数公式确定与密度同时进行反演的弹性参数：其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

说明书

基于球面波的叠前密度反演方法及装置

技术领域

[0001] 本发明涉及石油、天然气能源的地球物理勘探技术领域，特别涉及一种基于球面波的叠前密度反演方法及装置。

背景技术

[0002] 地震勘探领域中，密度是区分地下岩石物性的重要依据，并且通过密度和速度导出的如泊松比、流体因子等物性参数也是地下储层岩性、孔隙性和流体成分的重要指示工具，对储层表征和油藏开发与检测等环节非常重要。但现阶段的基于平面波理论的AVO(AVA)反演方法，由于能够利用的信息量有限，一般难以获取独立的密度估计，而多采用两参数耦合反演策略，如反演截距梯度和纵横波阻抗等，并采用多种正则化方法以增强反演的稳定性和可靠性。此外，反演密度信息是基于平面波AVO(AVA)理论反演方法的瓶颈，基于平面波AVO(AVA)理论反演方法一定要应用多个角度的地震资料，并且还要用到大角度地震资料，然而，大角度地震资料的获取和处理需要格外用心，稍有不慎，地震资料的质量就会影响反演的效果。

发明内容

[0003] 本发明实施例提供了一种基于球面波的叠前密度反演方法，以解决现有技术中基于平面波理论的AVO的密度参数反演方法一定要应用多个角度的地震资料、并且还要用到大角度地震资料的技术问题。该方法包括：获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演。

[0004] 在一个实施例中，所述地震资料是入射角小于30度的地震资料，所述地震资料是界面深度小于500米的地震资料。

[0005] 在一个实施例中，通过以下公式计算所述理论地震资料中各频点的球面波反射系数：

[0006]

[0007] 其中， 是所述理论地震资料中各频点的球面波反射系数，是平面波反射系数，平面波对应的入射角是
θ，x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。

[0008] 在一个实施例中，采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演，包括：采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数同时进行密度比和速度比反演。

[0009] 在一个实施例中，在采用非常快速模拟退火寻优方法，对所述球面波反射系数进行频率域的叠前密度反演之前，还包括：

[0010] 通过以下目标函数公式确定与密度同时进行反演的弹性参数：

[0011]

[0012] 其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

[0013] 本发明实施例还提供了一种基于球面波的叠前密度反演装置，以解决现有技术中基于平面波理论反演方法一定要应用多个角度的地震资料、并且还要用到大角度地震资料的技术问题。该装置包括：获取模块，用于获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；计算模块，用于计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；反演模块，用于采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演。

[0014] 在一个实施例中，所述地震资料是入射角小于30度的地震资料，所述地震资料是界面深度小于500米的地震资料。

[0015] 在一个实施例中，所述计算模块具体用于，通过以下公式计算所述理论地震资料中各频点的球面波反射系数：

[0016]

[0017] 其中， 是所述理论地震资料中各频点的球面波反射系数，是平面波反射系数，平面波对应的入射角是
θ，x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。

[0018] 在一个实施例中，所述反演模块具体用于，采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数同时进行密度比和速度比反演。

[0019] 在一个实施例中，还包括：反演确定模块，用于在采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演之前，通过以下目标函数公式确定与密度同时进行反演的弹性参数：

[0020]

[0021] 其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

[0022] 在本发明实施例中，通过获取入射角为单角度的地震资料，该地震资料包括低频段的频点，进而计算观测地震资料的球面波反射系数和理论地震资料的球面波反射系数，即获得了包括低频段频点的球面波反射系数，最后，采用非常快速模拟退火算法，对计算得到的观测地震资料的球面波反射系数和理论地震资料的球面波反射系数进行频率域的叠前密度反演，实现了基于低频段频点的球面波反射系数的叠前密度反演。由于球面波理论能够更精确地反映地下地震波的传播规律，低频段的频点的球面波反射系数随角频率变化而变化，即引入了低频频变信息，可以丰富反演过程中可采用的信息量，利用低频频变信息弥补了传统叠前密度反演信息量不足的缺点，并利用球面波反射系数的频变效应(即随角频率变化而变化)是弹性参数的客观反映这一事实，使得基于低频频变信息的叠前密度反演只需要入射角为单角度(小角度)的地震资料即可实现，且对入射角的角度大小并无特别要求。同时，采用非常快速模拟退火寻优方法进行反演，以求取准确的反演结果，减小反演结果的不确定性。

附图说明

[0023] 此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

[0024] 图1是本发明实施例提供的一种基于球面波的叠前密度反演方法的流程图；

[0025] 图2是本发明实施例提供的一种两层介质模型示意图；

[0026] 图3是本发明实施例提供的一种平面波反射系数曲线的示意图；

[0027] 图4是本发明实施例提供的一种球面波反射系数曲线的示意图；

[0028] 图5是本发明实施例提供的一种反演参数多次反演的结果示意图；

[0029] 图6是本发明实施例提供的一种反演参数多次反演的目标函数收敛曲线示意图；

[0030] 图7是本发明实施例提供的一种基于球面波的叠前密度反演装置的结构框图。

具体实施方式

[0031] 为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

[0032] 在油气勘探程度日益提高的情况下，利用现有资料实现对已开发或即将开发的油气资源的最大化利用是地球物理方法需要解决的关键问题，解决这一问题需要更准确模拟地下特征，更有针对性地选取方法和利用更多地下地震波提供的资源。现有技术中基于平面波理论的AVO(Amplitude Variation With Offset，振幅随偏移距的变化)(非垂直入射振幅关系，简称为AVA)反演方法，要应用多个角度的地震资料，并且还要用到大角度地震资料。但是，发明人发现，在AVO(AVA)分析和反演领域，球面波理论较平面波理论具有两个主要优势，一：球面波反射系数可以准确描述近或超临界角(大角度)处的反射波振幅和相位信息；二：球面波反射系数在低频段反射系数随频率变化的特征明显。考虑到球面波反射系数随频率变化是弹性参数的客观反映，并且对所用资料的角度没有严格限制，避开了易受动校正拉伸、调谐效应影响的大角度资料，发明人提出利用球面波反射系数的频变特征进行小角度资料频率域弹性参数反演，具体的方法即上述基于球面波的叠前密度反演方法。

[0033] 在本发明实施例中，提供了一种基于球面波的叠前密度反演方法，如图1所示，该方法包括：

[0034] 步骤101：获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；

[0035] 步骤102：计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；

[0036] 步骤103：采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演。

[0037] 由图1所示的流程可知，在本发明实施例中，通过获取入射角为单角度的地震资料，该地震资料包括低频段的频点，进而计算观测地震资料的球面波反射系数和理论地震资料的球面波反射系数，即获得了包括低频段频点的球面波反射系数，最后，采用非常快速模拟退火算法，对计算得到的观测地震资料的球面波反射系数和理论地震资料的球面波反射系数进行频率域的叠前密度反演，实现了基于低频段频点的球面波反射系数的叠前密度反演。由于球面波理论能够更精确地反映地下地震波的传播规律，低频段的频点的球面波反射系数随角频率变化而变化，即引入了低频频变信息，可以丰富反演过程中可采用的信息量，利用低频频变信息弥补了传统叠前密度反演信息量不足的缺点，并利用球面波反射系数的频变效应(即随角频率变化而变化)是弹性参数的客观反映这一事实，使得基于低频频变信息的叠前密度反演只需要入射角为单角度(小角度)的地震资料即可实现，且对入射角的角度大小并无特别要求。同时，采用非常快速模拟退火寻优方法进行反演，以求取准确的反演结果，减小反演结果的不确定性。

[0038] 具体实施时，为了实现计算理论地震资料中各频点的球面波反射系数，在本实施例中，通过以下公式计算所述理论地震资料中各频点的球面波反射系数：

[0039]

[0040] 其中， 是所述理论地震资料中各频点的球面波反射系数，是平面波反射系数，平面波对应的入射角是
θ，x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。如图2所示，我们假设界面深度已知并且激发点和接收点的位置已知，对于两均匀半空间声学介质模型，接触平面处的球面波反射系数的表达式如公式(1)所示。为了方便后续分析，我们可以将球面波反射系数写成入射角和角频率的函数，如公式(2)所示：

[0041]

[0042] 是采用地震资料的入射角的角度，可以通过界面深度、接收点和激发点的位置进行计算，ω是角频率。ρ＝ρ2/ρ1是密度比，c是积分路径，

[0043] 显然，球面波反射系数与上、下层的密度比和速度参数有关，并且从公式(2)中可以看出，球面波反射系数具有随角频率变化而变化的特征。

[0044] 具体实施时，通过上述公式(1)或(2)计算得到理论地震资料的球面波反射系数后，我们可以根据球面波反射系数的频变特征反演上下介质的弹性参数，并且我们可以利用的信息将极大丰富。图3说明了这一问题，在图3中示出的是传统基于平面波反射系数AVO(AVA)分析和反演所利用的信息点，可以看出平面波反射系数曲线随入射角变化，利用平面波反射系数发展出的AVO(AVA)分析和反演技术已经获得了极大的成功，并在工业界广泛应用。而图4是球面波反射系数随角频率和入射角变化的示意图，可见，相同状态下，传统基于平面波反射系数的AVO(AVA)方法利用了3个点的信息(如图3所示，)，而基于球面波反射系数进行反演可以利用9个点的信息量(如图4所示)，即使利用某一入射角 我们还是可以利用不同角频率的信息获取界面上、下两侧的弹性参数信息，基于球面波反射系数发展的反演方法也将成为有力的工具。

[0045] 具体实施时，现阶段利用球面波理论的反演均基于两弹性/声学介质，已知上层或下层的某些参数，反演上下介质三参数(纵横波速度和密度)中未知的参数，还是没有利用低频频变信息，并且参数选择随意，未有理论和数据的支撑。在本申请本实施例中，可以根据频变球面波反射系数(即球面波反射系数随角频率变化)的表达式，先通过可行性分析确定欲反演的弹性参数组合，并在反演之前给出反演关键参数的选取策略(例如，哪些参数和密度同时进行反演可以减小反演结果的不确定性)，最终实现有效的密度反演。例如，在采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行叠前密度反演之前，还包括：

[0046] 通过以下目标函数公式确定与密度同时进行反演的弹性参数：

[0047]

[0048] 其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

[0049] 具体的，采用的模型如图2所示，在不同的选定频率下，将利用公式(1)计算获得的理论地震资料的球面波反射系数与正演获得的观测地震资料的球面波反射系数进行匹配，将每个频点处的理论地震资料的球面波反射系数结果与观测地震资料的球面波反射系数的残差的2范数作为目标函数，具体如公式(3)所示：

[0050]

[0051] 在子波已知的前提下，可以通过傅里叶变换获得观测地震资料中各个频点的球面波反射系数记录AVAOBS(ωn)，利用公式(1)可以计算得到理论地震资料中各个频点的球面波反射系数记录AVATHEO(ωn)。

[0052] 在确定目标函数后，可以先考察不同反演参数组合对目标函数的影响，这样可以确定一组反演结果不确定性最小的弹性参数组合，即反演最可靠的弹性参数组合。具体的做法是对不同弹性参数组合的扰动试验。m代表同时反演的弹性参数的集合，针对每一组观测地震数据来说应该有一组真解m与之对应，反演是通过匹配观测地震数据和理论计算地震数据以求取m'，如果m'使得F(m)达到极小值，m'即为对真解m的估计，即获得了参数的反演结果。但是不同的m'会产生不同的F(m)，如果F(m)的形态很平缓，平缓可以理解成当F(m)小于某一小值З时，m'的可接受的变化范围很大，反演方法很难获得好的估计。因此，通过测试不同弹性参数组合的m'的F(m)的形态来确定反演弹性参数组合。在一定范围内扰动(m'＝mtrue*(1±20％))，代入公式(1)，计算理论地震资料的正演球面波反射系数记录AVATHEO(ωn)，然后代入公式(3)中获得目标函数值F(m)。

[0053] 例如，当我们扰动上、下层速度时，反演的目标函数是一个“峡谷”的特征，目标函数的等值线呈条带分布，说明组成目标函数的参数组合的线性相关程度比较高，在一定的扰动范围(例如，20％)内，未出现局部极小点。而且目标函数并没有收敛成理想的圆或者是椭圆，说明同时反演上、下层速度的不确定性很高，而这种不确定性是上、下层速度同时反演产生的，因此，应采用上、下层速度比而非上、下层速度作为参数代表。而选取速度比和密度比同时进行反演时，目标函数形态良好，目标函数逐渐收敛为椭圆，因此，在进行反演时，应该采用速度比与密度比作为弹性参数组合同时进行反演，可以进一步减小反演结果的不确定性。采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行叠前密度反演，包括：采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行密度比和速度比同时反演。

[0054] 具体实施时，上述基于球面波的叠前密度反演方法基于频率域匹配，因此需要选择一定频点的数据进行反演，从原理上讲，利用低频频变特征反演弹性参数，需要有足够的低频成分，我们通过目标函数(上述公式(3)所示的目标函数)的形态也可以直观的证实这一点，如果选取15-40Hz频点进行目标函数的绘制，可以观察到目标函数的形态很差，极小值点不收敛，呈发散状；如果选取40Hz单频进行目标函数的绘制，可以观察到目标函数值几乎没有变化，无法有效反应参数变化，也就无法进行参数反演，也证明了低频信息对频变反演的重要性。加之，现阶段低频信息的重要意义被广泛认知，地震采集的资料中低频成分的信息越来越多，质量越来越好，这都有利于利用低频信息特征的反演方法的研究。

[0055] 然而球面波反演利用的球面波反射系数的计算比较复杂，耗时较多，如何正确的选择地震资料的频率点作为输入资料，使得反演更加稳定也是值得研究的课题。鉴于低频重要并且需要一定的中高频进行速度反演的控制，在本实施例中，上述地震资料包括低频段的频点、中频段的频点以及高频段的频点，具体的，可以选取低频段(低频段频点一般为2-12Hz)2Hz采样，中频段(中频段频点一般为12-37Hz)5Hz采样，高频段(高频段频点一般为
37-87Hz)10Hz采样的策略，充分利用低频信息，并合理抽稀中高频信息，即采用[2 4 6 8 
10 12 17 22 27 32 37 47 57 67 87]Hz的频点组合作为上述地震资料。

[0056] 具体实施时，在确定上述地震资料的频点组合后，还可以对不同入射角的角度、界面深度的反演目标函数进行分析。具体的，可以选择不同入射角、采用相同频率点代入公式(3)的反演目标函数，例如，选择界面深度为150m的一组频率点，界面深度为500m的一组频率点，入射角为10度的一组频率点，入射角为30度的一组频率点，入射角为50度的一组频率点。根据各组频率点对应的目标函数可以看出，入射角角度越大界面深度越浅，目标函数的形态越好。其中，界面深度在150m，入射角角度为50度时，目标函数的极小值最为收敛，可见界面深度较小时，入射角角度较大时反演效果会更好。由于球面波的频变特征，可以选取小到中角度的地震资料进行密度反演，且这些资料的信噪比高，有益于反演的进行，因此，对于基于球面波的叠前密度反演方法，上述地震资料是入射角可以小于30度的地震资料，相对于基于平面波反射系数的反演方法可以采用小于30度的地震资料，而基于平面波反射系数的反演方法却不能实现，上述地震资料是界面深度可以小于500米的地震资料。

[0057] 此外，为了比较平面波反射系数与球面波反射系数的差异，利用平面波反射系数公式作为目标函数(公式3)的正演公式去计算AVATHEO(ωn)，可以获得界面深度h＝150m，入射角角度为50度时的目标函数形态图。根据目标函数形态图可见，目标函数的最小值点与模型的真实值之间是存在差异的，这说明平面波理论无法描述频变现象，无法获得良好的反演效果(即目标函数最小值没有收敛到真实值)。

[0058] 具体实施时，上述基于球面波的叠前密度反演方法采用非线性全局优化反演算法，发明人进行了多次试验，验证反演的可行性和准确性。具体的，本专利中采用的是非常快速模拟退火(Very Fast Simulated Annealing，简称为VFSA)寻优方法，首先输入选定频点(即上述地震资料的各频点)处的观测球面波频变反射系数AVAOBS(ωn)作为观测数据，并给定任一组参数的初始值m0＝(α2/α1,ρ)0∈[mmin,mmax]，m*为参数的真实值，mmin和mmax是参数的上下限。设置数据匹配误差限ε和最大迭代次数itmax，计算系统初始能量，即目标函数对于第i次迭代，利用随机扰动算子Pi可以更新参数为：mi＝mi-1+Δmi且mi∈[mmin,mmax]，Δmi＝Pi×(mmax-mmin)，令ξ∈(0,1)为满足均匀|2ξ-1|
分布的随机数，则Pi＝Tisign(ξ-0.5)[(1+1/Ti) -1]，然后合成新参数对应的频变反射系数理论计算资料AVATHEO(ωn,m0)，更新系统能量：
根据概率接受准则判断是否接受本次
更新的参数值，如果接受且Ei＜{Ek}min，k＝0,1,2,3,…,i-1，则将当前接受的参数更新值存入到最优解mopt记录器中。如果Ei＞ε且i＜itmax，则降低系统温度，继续迭代；否则输出记录器中的最优解mopt即参数的估计。通过模拟自然界金属退火结晶的过程(即上述非常快速模拟退火寻优方法)，采用的搜索方法更符合自然现象，可以获得全局最优解，确保了密度和速度信息反演的稳定性。并且可以通过输入不同的初始模型获得多个反演结果，以进行反演结果不确定性分析，并证实了提出的基于球面波低频频变信息的反演方法的准确性和可行性。

[0059] 具体的，以入射角为20度的地震资料为例时，可以得出基于球面波反射系数的非常快速模拟退火算法多次反演的反演结果(如图5所示)和目标函数的收敛情况(如图6所示)。多次实现是通过随机更改非常快速模拟退火算法最初给定的初始值来实现的。如图5所示，可知多次反演结果投影点集中在真实值投影点周围，误差(err)小于2％；如图6所示，可知目标函数的收敛曲线体现了反演收敛较快在2500次左右就可以达到收敛条件，收敛曲线收敛快速正确，说明了反演的可行性和准确性。

[0060] 如果我们将反演结果投影到目标函数图并将之放大显示，可知目标函数收敛很好，并且反演获得的密度和速度的估计非常稳定。使得可以利用单一角度地震数据反演获得了密度和速度的同时反演，这是传统AVO(AVA)反演无法实现的，说明了利用传统技术中忽略的低频频变信息进行密度反演是可行的也是有效的。

[0061] 基于同一发明构思，本发明实施例中还提供了一种基于球面波的叠前密度反演装置，如下面的实施例所述。由于基于球面波的叠前密度反演装置解决问题的原理与基于球面波的叠前密度反演方法相似，因此基于球面波的叠前密度反演装置的实施可以参见基于球面波的叠前密度反演方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

[0062] 图7是本发明实施例的基于球面波的叠前密度反演装置的一种结构框图，如图5所示，包括：获取模块701、计算模块702和反演模块703，下面对该结构进行说明。

[0063] 获取模块701，用于获取入射角为单角度的地震资料，所述地震资料包括通过直接观测获取的观测地震资料和通过计算得到的理论地震资料，所述地震资料包括低频段的频点；

[0064] 计算模块702，用于计算所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数，其中，低频段的频点的球面波反射系数随角频率变化而变化；

[0065] 反演模块703，用于采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演。

[0066] 在一个实施例中，所述地震资料是入射角小于30度的地震资料，所述地震资料是界面深度小于500米的地震资料。

[0067] 在一个实施例中，所述计算模块具体用于，通过以下公式计算所述理论地震资料中各频点的球面波反射系数：

[0068]

[0069] 其中， 是所述理论地震资料中各频点的球面波反射系数，是平面波反射系数，平面波对应的入射角是
θ，x＝cos(θ)，ω是角频率，J0是零阶贝塞尔函数，r是偏移距，h是激发源到反射平面的距离，zr是接收点到界面的距离，α1是上层的纵波速度，α2是下层的纵波速度，ρ1是上层的密度，ρ2是下层的密度，i是虚数单位。

[0070] 在一个实施例中，所述反演模块具体用于，采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数同时进行密度比和速度比反演。

[0071] 在一个实施例中，还包括：反演确定模块，用于在采用非常快速模拟退火算法，对所述观测地震资料的球面波反射系数和所述理论地震资料的球面波反射系数进行频率域的叠前密度反演之前，通过以下目标函数公式确定与密度同时进行反演的弹性参数：

[0072]

[0073] 其中，其中，F(m)表示所述观测地震资料与所述理论地震资料的匹配程度，N表示所选的参与计算的频点个数，m表示同时进行反演的弹性参数的集合，AVAOBS(ωn)是通过傅里叶变换获得的所述观测地震资料中各个频点的球面波反射系数记录，AVATHEO(ωn)是计算得到的所述理论地震资料中各个频点的球面波反射系数记录，n是所选频点的编号，使F(m)的数值小于预设小值З的弹性参数的变化范围的大小与反演结果的不确定性成正比，确定在F(m)小于预设小值З时，确定所有弹性参数中变换范围最小的弹性参数与密度同时进行反演。

[0074] 在本发明实施例中，通过获取入射角为单角度的地震资料，该地震资料包括低频段的频点，进而计算观测地震资料的球面波反射系数和理论地震资料的球面波反射系数，即获得了包括低频段频点的球面波反射系数，最后，采用非常快速模拟退火算法，对计算得到的观测地震资料的球面波反射系数和理论地震资料的球面波反射系数进行频率域的叠前密度反演，实现了基于低频段频点的球面波反射系数的叠前密度反演。由于球面波理论能够更精确地反映地下地震波的传播规律，低频段的频点的球面波反射系数随角频率变化而变化，即引入了低频频变信息，可以丰富反演过程中可采用的信息量，利用低频频变信息弥补了传统叠前密度反演信息量不足的缺点，并利用球面波反射系数的频变效应(即随角频率变化而变化)是弹性参数的客观反映这一事实，使得基于低频频变信息的叠前密度反演只需要入射角为单角度(小角度)的地震资料即可实现，且对入射角的角度大小并无特别要求。同时，采用非常快速模拟退火寻优方法进行反演，以求取准确的反演结果，减小反演结果的不确定性。

[0075] 显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

[0076] 以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。