CN114721044B - 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统 - Google Patents

一种多频率接收函数和振幅比联合反演地壳结构的方法及系统 Download PDF

Info

Publication number
CN114721044B
CN114721044B CN202210424327.6A CN202210424327A CN114721044B CN 114721044 B CN114721044 B CN 114721044B CN 202210424327 A CN202210424327 A CN 202210424327A CN 114721044 B CN114721044 B CN 114721044B
Authority
CN
China
Prior art keywords
data
amplitude ratio
longitudinal wave
function
frequency receiving
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
CN202210424327.6A
Other languages
English (en)
Other versions
CN114721044A (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 CN202210424327.6A priority Critical patent/CN114721044B/zh
Publication of CN114721044A publication Critical patent/CN114721044A/zh
Application granted granted Critical
Publication of CN114721044B publication Critical patent/CN114721044B/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/282—Application of seismic models, synthetic seismograms
    • 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
    • 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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/66—Subsurface modeling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

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

Abstract

本发明公开了一种多频率接收函数和振幅比联合反演地壳结构的方法及系统，该包括：获取待反演的多频率的接收函数和纵波振幅比观测数据；对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加；随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；建立多频率的接收函数和纵波振幅比观测数据和理论数据之间的目标函数，求取目标函数的全局最优解，确定最优地壳结构模型。本发明联合多频率接收函数和纵波振幅比进行反演，可以消除单频率数据集反演地壳结构的非唯一性问题，实现对地壳间断面和速度结构的双重约束，提高地壳结构的成像分辨率。

Description

一种多频率接收函数和振幅比联合反演地壳结构的方法及 系统
技术领域
本发明涉及地震数据处理技术领域，特别地涉及一种多频率接收函数和振幅比联合反演地壳结构的方法及系统。
背景技术
接收函数是基于等效震源假定，利用反褶积技术从远震三分量纵波波形数据中分离出的一种地震波形数据，其可应用于反演地下介质结构信息。由于接收函数反演方法具有较高的垂向分辨率，不仅可以获取台站下方的横波速度结构，而且对于速度间断面具有很好的约束。这一方法在研究地壳结构方面具有独特的优势，但是接收函数反演仍然存在高度的非线性和非唯一性问题。由于不同频率的接收函数波形振幅存在差异，为了更好的约束初始地壳结构模型，有学者提出利用不同频率接收函数迭代反演方法对地震台站下方地壳结构进行了研究，该方法首先利用低频率接收函数反演获取大尺度地下结构模型，作为下次高频率接收函数的输入模型，这样逐次迭代反演得到最终结果。由于线性反演方法本身存在着局限性，虽然这种递进式的反演策略可以约束不同尺度的速度结构，在一定程度上为接收函数反演提供了较为可靠的初始速度模型，但由于每一种尺度都采取了单一频率，且忽略了地震波形中的振幅信息，反演的地壳结构仍然存在不唯一性。
发明内容
有鉴于此，本发明提出一种多频率接收函数和振幅比联合反演地壳结构的方法及系统，以解决单一频率接收函数对地下结构约束能力较低的问题。
本发明第一方面提供一种多频率接收函数和振幅比联合反演地壳结构的方法，该方法包括：获取待反演的多频率的接收函数和纵波振幅比观测数据；对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加；随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；建立多频率的接收函数和纵波振幅比观测数据与接收函数和纵波振幅比理论数据之间的目标函数，求取目标函数的全局最优解，确定最优地壳结构模型。
进一步的，所述获取待反演的多频率的接收函数和纵波振幅比观测数据的步骤包括：获取三分量地震波形数据；利用时间域迭代反褶积方法对所述三分量地震波形数据进行处理，提取待反演的多频率的接收函数观测数据和多频率的纵波振幅比观测数据。
进一步的，所述对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加的步骤包括：采用互相关方法对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选，将挑选后的多频率的接收函数和纵波振幅比观测数据按射线参数进行分组叠加。
进一步的，所述随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据的步骤包括：获取最优深度-波速比参数，设定反演搜索空间；在设定的反演搜索空间内，随机生成若干地壳结构模型，并正演多频率的接收函数和纵波振幅比理论数据。
进一步的，所述地壳结构模型的生成方法包括：对多频率的纵向振幅比观测数据进行深度敏感核分析，得到最优时间窗口内的纵波振幅比数据；利用多频率的接收函数观测数据和最优时间窗口内的纵波振幅比数据进行反演，得到若干地壳结构模型。
进一步的，所述对多频率的纵向振幅比观测数据进行深度敏感核分析的步骤包括：选取某一秒时的纵波振幅比观测数据，分别对纵波、横波、密度和波速比进行振幅比深度敏感核分析；选取不同时间窗口的纵向振幅比观测数据对横波速度结构进行深度敏感核分析；确定最优时间窗口内的纵波振幅比数据。
进一步的，所述目标函数的建立方法包括：将每个地壳结构模型的多频率的接收函数和纵向振幅比观测数据与接收函数和纵向振幅比理论数据进行归一化，建立每个地壳结构模型对应的目标函数；采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，确定最优地壳结构模型。
进一步的，所述采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，确定最优地壳结构模型的步骤包括：采用差异演化算法对目标函数进行迭代求取，计算目标函数对应的拟合残差；判断拟合残差是否收敛，若拟合残差收敛，则比较所有目标函数对应的拟合残差值，获取拟合残差最小值，将拟合残差最小值作为目标函数的全局最优解；根据目标函数的全局最优解，将全局最优解的目标函数对应的地壳结构模型选作为最优地壳结构模型。
本发明第二方面提供一种多频率接收函数和振幅比联合反演地壳结构的系统，该系统包括：数据获取模块，用于获取待反演的多频率的接收函数和纵波振幅比观测数据；数据处理模块，用于对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加；模型建立模块，随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；模型筛选模块，用于建立多频率的接收函数和纵波振幅比观测数据与接收函数和纵波振幅比观测数据理论数据之间的目标函数，求取目标函数的全局最优解，确定最优地壳结构模型。
上述的多频率接收函数和振幅比联合反演地壳结构的方法及系统，利用时间域迭代反褶积技术提取多频率的接收函数和纵波振幅比观测数据；在给出反演搜索空间内，利用差异演化算法随机生成若干地壳结构模型；正演模拟每个地壳结构模型的不同频率的接收函数和纵波振幅比理论数据；计算每个模型中，同一频率的接收函数和纵波振幅比观测数据，和正演模拟的接收函数和纵波振幅比理论数据之间的拟合残差；如果拟合残差收敛，找到的拟合残差最小且收敛的全局最优解对应的最优地壳结构模型，即是最后反演的最佳地壳结构，可以消除单频率数据集反演地壳结构的非唯一性问题，实现对地壳间断面和速度结构的双重约束，提高地壳结构的成像分辨率。
附图说明
为了说明而非限制的目的，现在将根据本发明的优选实施例、特别是参考附图来描述本发明，其中：
图1是本发明一实施例提出的多频率接收函数和振幅比联合反演地壳结构的方法的流程图；
图2是依照本发明实施例的高斯因子(频率大小)对接收函数产生的影响示意图；
图3是依照本发明实施例的新引入纵波振幅比深度敏感核分析的示意图；
图4是依照本发明实施例的印度尼西亚一个地震台站记录的三分量地震波形数据提取的多频率接收函数和纵波振幅比数据的示意图；
图5是依照本发明实施例的印度尼西亚一个地震台站下方地壳结构的联合反演结果的示意图；
图6是本发明另一实施例提出的多频率接收函数和振幅比联合反演地壳结构的系统的结构示意图。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施例对本发明进行详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
在下面的描述中阐述了很多具体细节以便于充分理解本发明，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
图1是本发明一实施例提出的多频率接收函数和振幅比联合反演地壳结构的方法的流程图。该多频率接收函数和振幅比联合反演地壳结构的方法利用时间域迭代反褶积技术提取多频率的接收函数和纵波振幅比观测数据；在给出反演搜索空间内，利用差异演化算法随机生成若干地壳结构模型；正演模拟每个地壳结构模型的不同频率的接收函数和纵波振幅比理论数据；计算每个模型中，同一频率的实际观测接收函数和纵波振幅比数据，和正演模拟的接收函数和纵波振幅比理论数据之间的拟合残差；如果拟合残差收敛，找到的拟合残差最小且收敛的全局最优解对应的最优地壳结构模型，即是最后反演的最佳地壳结构。
请参阅图1，该多频率接收函数和振幅比联合反演地壳结构的方法包括以下步骤：
S100，获取待反演的多频率的接收函数和纵波振幅比观测数据。
在本实施例中，步骤S100的具体实现方式如下：
S101，获取三分量地震波形数据。
在本实施例中，步骤101中，获取三分量地震波形数据的具体实现方式为：
S101-1，获取原始地震波形数据。
S101-2，对原始地震波形数据进行预处理，包括去均值、去倾斜、带通滤波等。
S101-3，根据中国地震台网发布的全球地震目录，对预处理后的三分量地震波形数据进行截取，从而得到三分量地震波形数据。
本实施例在反演地壳结构之前，获取了预处理后的三分量地震波形数据，以便后续提取待反演的多频率的接收函数观测数据和多频率的纵波振幅比观测数据。
S102，提取待反演的多频率的接收函数观测数据和多频率的纵波振幅比观测数据。
本实施例利用时间域迭代反褶积技术，对三分量地震波形数据进行处理，提取待反演的不同频率的接收函数观测数据和不同频率的纵波振幅比观测数据。
地震波纵波振幅比数据反映了其在介质中传播的动力学特征，引入振幅信息约束地壳结构，对于从地震波场理论研究地壳结构具有重要意义。由于地震波形的复杂性以及数据信噪比的影响，从原始地震波形记录中直接识别地震体波并计算其纵波振幅比相对困难。通过理论推导发现，可以利用接收函数代替原始波形记录来测量视入射角，获取纵波振幅比信息，即基于接收函数的纵波振幅比法。由于该方法是利用反褶积技术去除震源时间函数、传播路径以及仪器响应的影响，仅保留与地震台站下方地质结构相关的信息，这种提取振幅比的方法，不仅减少了原始波形的复杂性而且信噪比也得到了明显提高。通过反褶积技术提取地震纵波振幅比，联合多频率接收函数和纵波振幅比进行反演，可以消除单频率接收函数反演地壳结构的非唯一性问题，实现对地壳间断面和速度结构的双重约束，提高地壳结构的成像分辨率。
准确的地壳速度结构对岩石岩性组分研究至关重要，然而仅利用单一频率数据集不能够反演地壳的精细结构，尤其是速度间断面和梯度信息。有鉴于此，本实施例首先模拟了不同频率的接收函数数据，分析其对不同地壳结构的响应特征。图2中，右侧图示的接收函数由左侧不同速度结构合成而来，其中每个速度结构合成高斯因子分别为1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5和5.0的接收函数。对于存在莫霍面过渡带的结构，随着接收函数频率的增加，接收函数的震相衰减很快(a)；对于存在低速带的结构，随着频率的增加，壳内的震相更加明显。所有合成理论接收函数时，选取的射线参数为0.06秒/千米，对应于60度震中距。通过该模拟可以发现，低频数据对深部大尺度结构敏感，而高频数据对浅部精细结构敏感，联合反演高低频数据即可获取浅部精细的速度结构，又兼顾了深部速度结构的成像。
准确的提取新引入的多频率纵波振幅比数据。理论推导和实际数据表明，平面纵波从均匀半空间入射到自由表面，实际入射角ip，视入射角
Figure BDA0003607885260000061
射线参数(水平慢度)p，纵波速度VP、横波速度VS之间关系如下：
Figure BDA0003607885260000062
根据上式可知，横波速度是关于纵波慢度与视入射角的函数，而视入射角可以根据实际地震波形记录直接进行测量。由于地震波形的复杂性以及数据信噪比的影响，从原始地震波形记录中直接识别地震体波并计算其水平-垂直比值(H/V)比较困难。本实施例利用接收函数代替原始波形记录来测量视入射角：
Figure BDA0003607885260000063
上式中，RRF和ZRF分别表示径向和垂向接收函数，T是测量采用的时间窗口长度的二分之一，该方法称为接收函数径向-垂向积分比值法(RFHV)，即基于远震接收函数的纵波振幅比法。
由于实际观测到的地震子波并不是脉冲函数，是具有一定宽度的复杂时间序列。对于一定宽度的地震子波信号，则有：
Figure BDA0003607885260000071
Figure BDA0003607885260000072
上式中，
Figure BDA0003607885260000073
表示纵波的地震信号；
Figure BDA0003607885260000074
表示横波的地震信号；RP与ZP和RS与ZS分别表示纵波和横波在径向与垂向分量上地震台站记录的频率域形式；f表示高斯截止频率，由高斯因子α控制(f＝α/π)；F-1表示反傅里叶变换；T为低通平滑因子；ω表示角频率；
Figure BDA0003607885260000075
为公式中RP/ZP和ZS/RS为1时的积分值。根据上式可知，RP/ZP和ZS/RS的计算实际上是一个反褶积过程，以此消除震源时间函数、仪器响应和传播路径的影响。
频率依赖性的控制参数。本实施例将引入高斯低通滤波
Figure BDA0003607885260000076
中的高斯因子α来控制多个频率范围。与面波频散类似，地震体波振幅比是关于频率的函数，低频数据对深部大尺度结构敏感，而高频数据对浅部精细结构敏感。因此联合多种频率的纵波振幅比数据可以有效降低反演结果的非唯一性，提高成像结果分辨率。
S200，对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加。
在本实施例中，采用互相关方法对提取的多频率的接收函数观测数据和纵波振幅比观测数据进行挑选，将挑选后的多频率的接收函数观测数据和纵波振幅比观测数据按射线参数进行分组叠加。这样反演出的结果是不同组内接收函数和纵波振幅比的平均效应，提升了多频率的接收函数和振幅比挑选的效率和质量。
本实施例通过时间域迭代反褶积技术提取多频率的接收函数和纵波振幅比观测数据，利用互相关方法对多频率的接收函数和纵波振幅比观测数据进行挑选，提升了多频率接收函数挑选的效率和质量。
S300，利用挑选及分组叠加后的多频率的接收函数和纵波振幅比观测数据，采用差异演化算法反演若干地壳结构模型，正演模拟多频率的接收函数和理论纵向振幅比理论数据。
在本实施例中，步骤S300的具体实现方式如下：
S301，给定反演搜索空间。
为了约束反演所用的搜索空间，利用深度-波速比叠加方法获取初始地壳结构模型的平均深度和波速比最佳取值范围，然后设定反演搜索空间，并在整个不断的测试中不断调整反演搜索空间，获取最优地壳结构模型拟合效果。
本实施例在反演搜索空间中施加深度-波速比参数约束，使得反演搜索空间得以实时更新，能够促进反演过程的收敛效率和提升最终反演结果质量。
给定反演搜索空间后，将界面的深度值和速度值作为模型参数向量，利用多频率的接收函数和纵波振幅比观测数据通过差异演化算法生成若干地壳结构模型。
S302，利用多频率的接收函数和纵波振幅比观测数据生成若干地壳结构模型。
在本实施例中，步骤S302中利用多频率的接收函数和纵波振幅比观测数据生成若干地壳结构模型的具体实现方式如下：
S3021，对多频率的纵向振幅比观测数据进行深度敏感核分析，得到最优时间窗口内的纵波振幅比数据。
首先，选取某一秒时的纵波振幅比观测数据，分别对纵波、横波、密度和波速比开展振幅比深度敏感核分析；同时，选取不同时间窗口的纵向振幅比观测数据对横波速度结构进行深度敏感核分析。
为了确认新引入的纵波振幅比对地壳结构具有约束作用，本实施例对理论数据进行深度敏感核分析，有助于实际数据联合反演时选取合适的参数。选取6秒时的纵波振幅比数据，分别对纵波、横波、密度和波速比开展振幅比深度敏感核分析，结果表明(如图3中(A)所示)，地壳模型深度范围以内纵波振幅比对横波速度敏感，因此可以用以反演地壳横波速度结构；同时，选取不同时间窗口的纵向振幅比数据(例如，1秒、6秒、10秒和20秒)对横波速度结构进行深度敏感核分析，结果表明(如图3中(B)所示)，10秒以内的振幅比数据对地壳横波速度都具有较高的敏感性，因此本实施例采用10秒以内的纵波振幅比数据参与联合反演工作。
S3013，利用多频率的接收函数观测数据和最优时间窗口内的纵波振幅比数据进行反演，得到若干地壳结构模型。
S303，正演模拟多频率的接收函数理论数据和纵向振幅比理论数据。
本实施例可以在无地下介质模型先验信息情况下，进行地壳结构的参数建模。
S400，建立多频率的接收函数和纵向振幅比观测数据与接收函数和纵向振幅比理论数据之间的目标函数，求取目标函数的全局最优解，获取最优地壳结构模型。
在本实施例中，步骤S400中建立目标函数的具体实现方式如下：
S401，将每个地壳结构模型的多频率接收函数和纵向振幅比观测数据与多频率的接收函数和纵向振幅比理论数据进行归一化，建立每个地壳结构模型对应的目标函数。
本实施例用目标函数来定量计算多频率的接收函数和纵向振幅比理论数据与多频率的接收函数和振幅比纵向观测数据之间的拟合残差。目标函数δ表达式如下：
Figure BDA0003607885260000091
其中，δ表示联合反演中，理论数据与实际观测数据的失配值(目标函数)；m表示数据类型，可以根据实际观测数据类型选择(例如，m＝1表示多频率接收函数数据，m＝2表示多频率纵波振幅比数据)；Wm表示数据类型m的权重系数；WS为光滑因子权重；Nm表示每种数据类型用于反演的数目；i表示第i个数据；j表示数据的第j个采样点；NPTim表示m类型数据第i个数据的采样点数目；
Figure BDA0003607885260000092
表示m类型数据第i个数据第j个采样点的观测值
Figure BDA0003607885260000093
表示m类型数据第i个数据第j个采样点的理论值σim(j)表示m类型数据第i个数据第j个采样点的测量误差。
在目标函数迭代求取中，为了避免异常的速度跳跃，本发明引入了光滑因子：
Figure BDA0003607885260000094
式中，Nlarer是地壳结构模型的层数，Vj是第j层的速度，Vj+1是第j+1层的速度，Vmean是地壳平均速度，n表示第n层。
图4展示了本实施例所使用的印度尼西亚一个地震台站的多频率接收函数和纵波振幅比数据，一共使用了9个不同高斯因子(频率大小)的接收函数和纵波振幅比数据，横轴表示反演的地壳横波速度结构，纵轴表示深度(千米)，models表示联合反演中绘制的模型数目图例。
S402，采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，寻找最优地壳结构模型。
在本实施中，步骤S403中采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，寻找最优地壳结构模型的具体实现方式如下：
采用差异演化算法对目标函数进行迭代求取，计算每个地壳结构模型中同一频率的接收函数和纵向振幅比实际观测数据，和正演模拟的接收函数和纵向振幅比理论数据之间的拟合残差；
判断拟合残差是否收敛，若拟合残差收敛，则比较所有目标函数对应的拟合残差值，获取拟合残差最小值，将拟合残差最小值作为目标函数的全局最优解；若拟合残差不收敛，重复步骤S100-S400，直至拟合残差收敛。
根据目标函数的全局最优解，将全局最优解的目标函数对应的地壳结构模型选作为最优地壳结构模型。
图5展示了本实施例的印度尼西亚一个地震台站的联合反演最优地壳结构模型，右侧曲线表示联合反演的最优地壳结构，左侧曲线表示联合反演的地壳波速比信息，黑白相间的背景区域表示联合反演的搜索范围，该成像结构清晰的刻画了地壳结构在3公里左右存在一个小的低速层，同时对地壳结构24至28公里莫霍面过渡带也进行了很好的成像。
S500，利用最优地壳结构模型对待反演的地壳结构进行成像。
上述的多频率接收函数和振幅比联合反演地壳结构的方法，利用多频率的接收函数数据和纵向振幅比数据联合反演地壳结构，新引入多频率纵向振幅比数据对地壳结构进行约束，可以兼顾大尺度宏观结构与小尺度细节构造，最终获得地下介质高精度地壳结构。
上述的多频率接收函数和振幅比联合反演地壳结构的方法，利用反褶积方法获取待反演的多频率的接收函数和纵波振幅比观测数据，给定反演搜索空间，利用差异演化算法随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据，将观测数据和理论数据归一化到同一个目标函数，基于差异演化算法求取目标函数的全局最优收敛解，获取最优地壳结构模型，对地下介质结构进行精细成像。
上述的多频率接收函数和振幅比联合反演地壳结构的方法，联合多频率的接收函数和纵向振幅比数据进行反演，可以消除单频率数据集反演地壳结构的非唯一性问题，实现对地壳间断面和速度结构的双重约束，提高地壳结构的成像分辨率。
图6是本发明另一实施例提供的多频率接收函数和振幅比联合反演地壳结构的系统的结构示意图。
在本实施例中，所述多频率接收函数和振幅比联合反演地壳结构的系统100可以应用于计算机装置中，所述多频率接收函数和振幅比联合反演地壳结构的系统100可以包括多个由程序代码段所组成的功能模块。所述多频率接收函数和振幅比联合反演地壳结构的系统中的各个程序段的程序代码可以存储于计算机装置的存储器中，并由所述计算机装置的至少一个处理器所执行，以实现(详见图1描述)多频率接收函数和振幅比联合反演地壳结构功能。
本实施例中，所述多频率接收函数和振幅比联合反演地壳结构的系统100根据其所执行的功能，可以被划分为多个功能模块。所述功能模块可以包括：
数据获取模块101，用于获取待反演的多频率的接收函数和纵波振幅比观测数据；
数据处理模块102，用于对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加；
模型建立模块103，随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；
模型筛选模块104，用于建立多频率的接收函数和纵波振幅比观测数据与接收函数和纵波振幅比理论数据之间的目标函数，求取目标函数的全局最优解，获取最优地壳结构模型。
本发明所称的模块是指一种能够被至少一个处理器所执行并且能够完成固定功能的一系列计算机程序段，其存储在存储器中。在本实施例中，关于各模块的具体功能不再赘述。
上述具体实施方式，并不构成对本发明保护范围的限制。本领域技术人员应该明白的是，取决于设计要求和其他因素，可以发生各种各样的修改、组合、子组合和替代。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明保护范围之内。

Claims (6)

1.一种多频率接收函数和振幅比联合反演地壳结构的方法，其特征在于，包括：
获取待反演的多频率的接收函数和纵波振幅比观测数据；
对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加；
随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；
所述随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据的步骤包括：
获取最优深度-波速比参数，设定反演搜索空间；
在设定的反演搜索空间内，随机生成若干地壳结构模型，并正演多频率的接收函数和纵向振幅比理论数据；
其中，所述地壳结构模型的生成方法包括：
对多频率的纵向振幅比观测数据进行深度敏感核分析，得到最优时间窗口内的纵波振幅比数据；
利用多频率的接收函数观测数据和最优时间窗口内的纵波振幅比数据进行反演，得到若干地壳结构模型；
建立多频率的接收函数和纵波振幅比观测数据与接收函数和纵波振幅比理论数据之间的目标函数，求取目标函数的全局最优解，确定最优地壳结构模型；所述目标函数的建立方法包括：
将每个地壳结构模型的多频率的接收函数和纵向振幅比观测数据与接收函数和纵向振幅比理论数据进行归一化，建立每个地壳结构模型对应的目标函数；
采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，确定最优地壳结构模型；其中，所述目标函数的表达式为：
Figure FDA0003909996680000011
其中，δ表示目标函数；m表示数据类型；Wm表示数据类型m的权重系数；Ws为光滑因子权重；Nm表示每种数据类型用于反演的数目；i表示第i个数据；j表示数据的第j个采样点；NPTim表示m类型数据第i个数据的采样点数目；
Figure FDA0003909996680000021
表示m类型数据第i个数据第j个采样点的观测值；
Figure FDA0003909996680000022
表示m类型数据第i个数据第j个采样点的理论值；σim(j)表示m类型数据第i个数据第j个采样点的测量误差；Smoothness表示光滑因子。
2.根据权利要求1所述的多频率接收函数和振幅比联合反演地壳结构的方法，其特征在于，所述获取待反演的多频率的接收函数和纵波振幅比观测数据的步骤包括：
获取三分量地震波形数据；
利用时间域迭代反褶积方法对所述三分量地震波形数据进行处理，提取待反演的多频率的接收函数和纵波振幅比观测数据。
3.根据权利要求1所述的多频率接收函数和振幅比联合反演地壳结构的方法，其特征在于，所述对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选并分组叠加的步骤包括：
采用互相关方法对待反演的多频率的接收函数和纵波振幅比观测数据进行挑选，将挑选后的多频率的接收函数和纵波振幅比观测数据按射线参数进行分组叠加。
4.根据权利要求1所述的多频率接收函数和振幅比联合反演地壳结构的方法，其特征在于，所述对多频率的纵向振幅比观测数据进行深度敏感核分析的步骤包括：
选取某一秒时的纵波振幅比观测数据，分别对纵波、横波、密度和波速比进行振幅比深度敏感核分析；
选取不同时间窗口的纵向振幅比观测数据对横波速度结构进行深度敏感核分析；
确定最优时间窗口内的纵波振幅比数据。
5.根据权利要求1所述的多频率接收函数和振幅比联合反演地壳结构的方法，其特征在于，所述采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，确定最优地壳结构模型的步骤包括：
采用差异演化算法对目标函数进行迭代求取，计算目标函数对应的拟合残差；
判断拟合残差是否收敛，若拟合残差收敛，则比较所有目标函数对应的拟合残差值，获取拟合残差最小值，将拟合残差最小值作为目标函数的全局最优解；
根据目标函数的全局最优解，将全局最优解的目标函数对应的地壳结构模型选作为最优地壳结构模型。
6.一种多频率接收函数和振幅比联合反演地壳结构的系统，其特征在于，包括：
数据获取模块，用于获取待反演的多频率的接收函数和纵波振幅比观测数据；
数据处理模块，用于对待反演的多频率的接收函数和纵波振幅比数据进行挑选并分组叠加；
模型建立模块，用于随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据；所述随机生成若干地壳结构模型，正演多频率的接收函数和纵波振幅比理论数据的步骤包括：
获取最优深度-波速比参数，设定反演搜索空间；
在设定的反演搜索空间内，随机生成若干地壳结构模型，并正演多频率的接收函数和纵向振幅比理论数据；
其中，所述地壳结构模型的生成方法包括：
对多频率的纵向振幅比观测数据进行深度敏感核分析，得到最优时间窗口内的纵波振幅比数据；
利用多频率的接收函数观测数据和最优时间窗口内的纵波振幅比数据进行反演，得到若干地壳结构模型；
模型筛选模块，用于建立多频率的接收函数和纵波振幅比观测数据与接收函数和纵波振幅比理论数据之间的目标函数，求取目标函数的全局最优解，确定最优地壳结构模型；
所述目标函数的建立方法包括：
将每个地壳结构模型的多频率的接收函数和纵向振幅比观测数据与接收函数和纵向振幅比理论数据进行归一化，建立每个地壳结构模型对应的目标函数；
采用差异演化算法对所有目标函数进行迭代反演，求取目标函数极小值，获取目标函数的全局最优解，确定最优地壳结构模型；其中，所述目标函数的表达式为：
Figure FDA0003909996680000041
其中，δ表示目标函数；m表示数据类型；Wm表示数据类型m的权重系数；Ws为光滑因子权重；Nm表示每种数据类型用于反演的数目；i表示第i个数据；j表示数据的第j个采样点；NPTim表示m类型数据第i个数据的采样点数目；
Figure FDA0003909996680000042
表示m类型数据第i个数据第j个采样点的观测值；
Figure FDA0003909996680000043
表示m类型数据第i个数据第j个采样点的理论值；σim(j)表示m类型数据第i个数据第j个采样点的测量误差；Smoothness表示光滑因子。
CN202210424327.6A 2022-04-21 2022-04-21 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统 Active CN114721044B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210424327.6A CN114721044B (zh) 2022-04-21 2022-04-21 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210424327.6A CN114721044B (zh) 2022-04-21 2022-04-21 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统

Publications (2)

Publication Number Publication Date
CN114721044A CN114721044A (zh) 2022-07-08
CN114721044B true CN114721044B (zh) 2023-03-10

Family

ID=82245770

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210424327.6A Active CN114721044B (zh) 2022-04-21 2022-04-21 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统

Country Status (1)

Country Link
CN (1) CN114721044B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118295017B (zh) * 2024-05-10 2024-10-01 成都理工大学 基于密集台阵多频接收函数的莫霍面形态反演方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106886047A (zh) * 2017-02-28 2017-06-23 中国地质大学（北京） 一种接收函数和重力联合反演地壳厚度和波速比的方法
CN107203002B (zh) * 2017-06-12 2019-05-24 中国科学院地质与地球物理研究所 反演速度模型的建立方法和地下结构的像的获得方法
CN108873103A (zh) * 2018-09-14 2018-11-23 吉林大学 一种结构约束的二维重力梯度和大地电磁联合反演方法
CN110221344B (zh) * 2019-06-17 2020-08-28 中国地质大学（北京） 一种地壳三维密度结构的地震全波形与重力联合反演方法
CN110888159B (zh) * 2019-11-15 2021-08-06 西安理工大学 基于角度分解与波场分离的弹性波全波形反演方法
CN112698390B (zh) * 2020-11-11 2022-12-02 中国石油天然气股份有限公司 叠前地震反演方法及装置
CN113740915B (zh) * 2021-07-28 2024-04-19 中国人民武装警察部队黄金第五支队 一种球坐标系重力和接收函数联合反演地壳结构参数的方法
CN113671570B (zh) * 2021-08-23 2022-04-19 湖南工商大学 一种地震面波走时和重力异常联合反演方法与系统

Also Published As

Publication number Publication date
CN114721044A (zh) 2022-07-08

Similar Documents

Publication Publication Date Title
CN110058303B (zh) 声波各向异性逆时偏移混合方法
CN107065013B (zh) 一种地震尺度下的层速度确定方法及装置
US7460437B2 (en) Seismic data processing method and system for migration of seismic signals incorporating azimuthal variations in the velocity
EA032186B1 (ru) Сейсмическая адаптивная фокусировка
EA001509B1 (ru) Способ обработки сигналов сейсмических данных
CN109655918B (zh) 地面浅井微地震监测观测台站位置确定方法及系统
CN113391351B (zh) 一种基于被动源地震波场分析提取矿集区结构的方法
CN112285778B (zh) 一种粘声TTI介质中纯qP波的逆时偏移成像方法
CN112285767A (zh) 海底地震仪四分量海洋面波多阶频散能量成像装置及方法
JP2000512385A (ja) 伝搬する波動場のサンプリング及び復元
CN114721044B (zh) 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统
KR20170009609A (ko) 반복적 직접 파형 역산 및 완전 파형 역산을 이용한 탄성파 영상화 장치 및 방법
CN111665556B (zh) 地层声波传播速度模型构建方法
CN114415234A (zh) 基于主动源面波频散和h/v确定浅地表横波速度的方法
CN110579799B (zh) 一种等旅行时间间隔的地震采集观测方法及系统
Li et al. Gaussian beam imaging of fractures near the wellbore using sonic logging tools after removing dispersive borehole waves
CN115598704A (zh) 一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质
CN111158053B (zh) 裂缝预测方法及装置
CN111665546B (zh) 用于可燃冰探测的声学参数获取方法
CN111665550A (zh) 地下介质密度信息反演方法
CN111665549A (zh) 地层声波衰减因子反演方法
CN111665551B (zh) 用于桥梁基底探测的声学参数获取方法
CN116125535B (zh) 三维vsp成像的方法及装置
CN113075734B (zh) 一种基于信噪比约束的剩余曲率谱计算方法及装置
CN112462428B (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