CN111164462B - 一种人工源面波勘探方法、面波勘探装置及终端设备 - Google Patents

一种人工源面波勘探方法、面波勘探装置及终端设备 Download PDF

Info

Publication number
CN111164462B
CN111164462B CN201880000953.4A CN201880000953A CN111164462B CN 111164462 B CN111164462 B CN 111164462B CN 201880000953 A CN201880000953 A CN 201880000953A CN 111164462 B CN111164462 B CN 111164462B
Authority
CN
China
Prior art keywords
surface wave
frequency
dispersion curve
frequency dispersion
function
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
CN201880000953.4A
Other languages
English (en)
Other versions
CN111164462A (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.)
Southwest University of Science and Technology
Original Assignee
Southwest University of Science and 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 Southwest University of Science and Technology filed Critical Southwest University of Science and Technology
Publication of CN111164462A publication Critical patent/CN111164462A/zh
Application granted granted Critical
Publication of CN111164462B publication Critical patent/CN111164462B/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/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • 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/284—Application of the shear wave component and/or several components of the seismic signal
    • 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
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • 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/32—Transforming one recording into another or one representation into another
    • 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/67—Wave propagation modeling
    • G01V2210/675—Wave equation; Green's functions

Landscapes

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

Abstract

一种人工源面波勘探方法，适用于地质勘探技术领域，该方法包括：通过预设台站处的检波器采集从震源传播过来的面波数据；基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图；从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线；根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。通过上述方法，有效提高了面波勘探结果的准确率。还提供一种面波勘探装置及终端设备。

Description

一种人工源面波勘探方法、面波勘探装置及终端设备
技术领域
本申请涉及地质勘探技术领域，尤其涉及一种人工源面波勘探方法、面波勘探装置及终端设备。
背景技术
面波，是地震波中一种特殊类型的波，是纵波和横波在震源区域内各界面处经过复杂的反射、透射后相互干涉叠加而成的。面波在传播过程中携带了大量的地层信息，呈现出频散特征，也能够间接反映出层状介质本身所固有的一些特征。因此，通常利用主动源面波对地层结构进行勘探。
但是目前利用主动源面波对地层结构进行勘探，需要观测系统和震源成线性排列，而且要求检波器等间距放置。在城市复杂区域经常无法达到上述施工条件，即使可以勉强施工，也无法得到高分辨率的面波勘探图像，进而无法得到准确的地层结构的勘探结果。
发明内容
有鉴于此，本申请实施例提供了一种人工源面波勘探方法、面波勘探装置及终端设备，以解决现有技术中面波勘探结果不准确的问题。
本申请实施例的第一方面提供了一种人工源面波勘探方法，包括：
通过预设台站处的检波器采集从震源传播过来的面波数据；
基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图；
从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线；
根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
本申请实施例的第二方面提供了一种面波勘探装置，包括：
采集单元，用于通过预设台站处的检波器采集从震源传播过来的面波数据；
计算单元，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图；
提取单元，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线；
反演单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
本申请实施例的第三方面提供了一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本申请实施例第一方面提供的所述方法的步骤。
本申请实施例的第四方面提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被一个或多个处理器执行时实现本申请实施例第一方面提供的所述方法的步骤。
本申请实施例通过预设台站处的检波器采集从震源传播过来的面波数据，可以任意布置采集装置，降低了对布置场地的要求，提高了面波勘探场地的适应性。通过利用提取出的频散曲线建立初始地层模型，降低了反演运算的计算时间；之后根据初始地层模型对基阶面波频散曲线和高阶面波频散曲线进行联合反演，将高阶面波频散的信息加入到地层的反演运算中，从而降低了反演运算的不确定性，提高了反演运算的鲁棒性。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是本申请实施例提供的人工源面波勘探方法中计算频散能量图的实现流程示意图；
图2是本申请实施例提供的面波勘探装置的示意图；
图3是本申请实施例提供的终端设备的示意图；
图4是本申请实施例提供的提取到的F-C频散谱(a)、频率区间分类(b)和格林函数核函数得到的理论F-C频散谱(c)的示意图；
图5是本申请实施例提供的深度域的频散谱和地层模型的示意图。
具体实施方式
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本申请实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本申请。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本申请的描述。
应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当理解，在此本申请说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本申请。如在本申请说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
还应当进一步理解，在本申请说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
如在本说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
为了说明本申请所述的技术方案，下面通过具体实施例来进行说明。
所述方法可以包括以下步骤：
步骤S101，通过预设台站处的检波器采集从震源传播过来的面波数据。
在实际应用中，采集从震源传播过来的面波数据所用的装置包括但不限于检波器，例如，可以采用多道有线连接的工程地震仪，或者独立无线连接的地震仪。优选地，检波器可以为主频不高于4hz的宽频带检波器，采集带宽越快越有利于各种频率的面波的采集。检波器的个数大于等于预设个数，例如检波器的数量大于等于12个。检波器的采样频率应满足勘探目的，工程勘探采样率一般不低于200hz。另外，预设台站是人为预先设定的，每个预设台站处放置有一台检波器。
其中，面波是地震波的一种，主要在地表传播，能量最大，波速约为3.8千米/秒，低于体波，往往最后被记录到。面波实际上是体波在地表衍生而成的次生波。面波的传播较为复杂，既可以引起地表上下的起伏，也可以是地表做横向的剪切，其中剪切运动对建筑物的破坏最为强烈。面波包括瑞雷波、拉夫波、水力波、斯通利尔波等。而经研究发现Rayleigh波(瑞雷波)在层状介质中相速度随频率改变而改变，呈现明显的频散特性。水平层状介质中的Rayleigh波实际上是纵波和横波在震源区域内各界面处经过复杂的反射、透射后相互干涉叠加而成。它携带了各层介质的P波速度、S波速度、密度等参数信息，且速度主要取决于层状介质中S波速度的分布。Rayleigh波在传播过程中能量和速度的变化特征携带了大量地下地层的信息，呈现出的频散特征，也间接反映了层状介质本身所固有的一些特征。由此研究天然地震波中的低频Rayleigh波频散可以解决深部地质构造问题；研究人工震源激发的较高频率的Rayleigh波可以解决工程勘察、场地和地基处理评价、障碍物和空洞探测等浅层地质问题。因此，本申请实施例中优选地采用瑞雷波。
步骤S102，基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图。
采用矢量波数变换算法，预设台站处的检波器可以不必按照预设规则(现有技术中要求检波器和震源线性排列且等间距放置)进行放置，这样降低了施工难度，增加了施工场地的适应性。
步骤S102的具体实施步骤可参见图2实施例中的描述。
步骤S103，从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线。
实际接收的背景噪声数据由各种震动产生的波组成，不仅包含面波，也包含了体波。而且面波在非均匀介质中会发生频散现象，即面波由不同相速度的模态组成。经大量研究证明，频散曲线中的高阶面波部分在地层结构反演分析中起着关键作用。利用本申请实施例中的矢量波数变换算法，通过计算得到的频散能量图就能有效分离出由不同速度的面波(即能够分离出基阶面波和高阶面波)和体波的组分。
步骤S104，根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
目前工程应用中的人工源面波勘探方法只在频散能量图中按能量极大值，手动或自动连接频散曲线，根据频散曲线中的“之”字型特征来反演地层深度和厚度。上述反演方法，必须对高阶面波的阶数有个一个准确的判断。但当地层存在低速层或高速层时，不但瑞雷波各个模态的能量分布发生变化，各个模态的速度随频率的变化也会发生改变，从而经常会产生“模式接吻(mode kissing)”现象，这样就会对高阶模态频散曲线的判断带来很大困难。而且当在高频范围内在水平层状的地层模型中存在软弱夹层时，高阶面波比基阶面波具有更大的能量，这就意味着在一定频率范围内通过目前的方法是无法得到基阶面波的，而仅能得到高阶面波。在实际勘探中，真实地层并不是理想的水平层状各向同性的结构，从而导致瑞雷波(Rayleigh)在频散谱波高阶模态的成像质量通常不高，以上这些因素都制约了利用高阶频散曲线进行反演。
经研究发现当地层存在低速或高速夹层时，夹层埋深对应的频率范围区间内Rayleigh波能量从基阶向一阶或更高阶模态阶跃，从而导致基阶和高阶面波频散曲线出现只在某一频率范围内连续，在实际数据中成像质量可能更差。因此在频散能量图中，各个模态能量的分布也与地层结构有着密切的联系。在本实施例中，根据在频散能量图中各个不同频率区间各模态的能量分布的关系，对频率区间进行分类，从而迅速建立简单的层状地层模型，作为后续精确反演的初始模型。
另外，在本实施例中，可以采用模拟退火算法、遗传算法等反演算法，对频散曲线进行反演，得到地层信息和/或振动波的速度信息。例如，可以得到地层深度信息和速度剖面，从而实现对地层结构的勘探；地层信息可以包括地层埋深、S波速度、密度、P波速度等。
在本申请实施例中，所述根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，包括：
根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类。
根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
如图4所示，图4a)为F-C频散能量图，图中点线为理论面波频散曲线；图4b)为根据频散能量提取得到的离散的频散点，并根据分布特性将频率区间分为4类；图4c)为格林函数核函数得到的理论F-C频散能量图。
将频率-速度域的频散点，根据半波长理论转换到深度-速度域，见图5。图5a)深度域的频散能量图，图5b)为地层模型。可以看到，频散能量图上1号和3号点线上的点为基阶频散曲线上的点，2号和4号点线上的点为高阶频散曲线上的点。将地层模型与得到的深度域的频散曲线进行对比，可以看到位移20-40m埋深的第三层(低速层)与深度-速度剖面中的4号频散点的分布基本一致；位移0-10m的第一层与深度-速度剖面中的2号频散点的分布也基本一致。这样我们可以看到高阶频散曲线的分布与地层确实存在一一对应的关系。由此证明，通过在频率域将频散点进行分类，对地层进行分层，建立初始建模的思路是正确的。
本申请实施例通过预设台站处的检波器采集从震源传播过来的面波数据，可以任意布置采集装置，降低了对布置场地的要求，提高了面波勘探场地的适应性。通过利用提取出的频散曲线建立初始地层模型，降低了反演运算的计算时间；之后根据初始地层模型对基阶面波频散曲线和高阶面波频散曲线进行联合反演，将高阶面波频散的信息加入到地层的反演运算中，从而降低了反演运算的不确定性，提高了反演运算的鲁棒性。
图1是本申请实施例提供的人工源面波勘探方法中计算频散能量图的实现流程示意图。如图所示，步骤S102可以包括以下步骤：
步骤S201，从所述面波数据中提取震源时间函数。
在实际应用中，可以有多种方法提取震源时间函数。例如，可以直接将震源时间函数近似为雷克子波；还可以从多道地震记录中根据相关性提取震源时间函数。只要能够提取出震源时间函数即可，不对提取方法做具体限定。
步骤S202，计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
其中，炮检距为震源与检波器之间的距离。因为检波器安装于预设台站处，所以炮检距也可看作是预设台站与震源之间的距离。在计算式中，可以用gzz表示预设台站与震源之间的格林函数。
步骤S203，根据所述震源时间函数和所述格林函数计算得到频散能量图。
在本申请实施例中，步骤S203还可以包括：
将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
其中，所述时间域的垂直分量数据为：
uz(x,t)＝f0(t)*gzz(x,t)
其中，uz(x,t)表示所述垂直分量数据，f0表示所述震源时间函数，gzz表示所述格林函数。
对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
其中，所述频率域的垂直分量数据为：
U(r,ω)＝F0(ω)G(r,ω)
式中，U为所述频率域的垂直分类数据，
Figure GDA0003467162410000081
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F0为频率域的震源时间函数，G为频率域的格林函数，
Figure GDA0003467162410000082
g(ω,k)是核函数，
Figure GDA0003467162410000083
为波数，J0(kr)是一类零阶贝塞尔函数。
对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
在本申请实施例中，所述对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图，还包括：
对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
基于所述最终计算式，计算得到频散能量图。
在实际应用中，当震源时间函数为雷克子波时，频率域中F0(ω)为纯实数函数，即可得预设台站接收到的面波数据的谱函数近似为格林函数的虚部，二者仅在幅值上有差异：
U(r,ω)＝A·{G(r,ω)}
其中，A＝F0(ω)为常量。对频率域的垂直分量数据U进行矢量波数变换，可得
Figure GDA0003467162410000091
将
Figure GDA0003467162410000092
带入上式，可得中间计算式为：
Figure GDA0003467162410000093
根据贝塞尔函数的正交性
Figure GDA0003467162410000094
可将上式化简得到最终计算式为：
Figure GDA0003467162410000095
这里核函数g(ω,k)具有如下特点：g(ω,k)的值反比于确定面波频散特性的久期函数值，即
Figure GDA0003467162410000096
其中，
Figure GDA0003467162410000097
Rdu为下行波的反射系数矩阵，Rud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。当k＝kn(ω)(n＝1,2,3,...)是核函数的g(ω,k)的极点时，核函数的值趋于无穷大。利用这一性质，提出了矢量波数变换法(VectorWavenumberTransformMethod，VWTM)提取频散曲线。
基于上述最终计算式，即可计算得到频散能量图。
本申请实施例提出了一种矢量波数变换算法，并基于这种矢量波数变换算法根据面波数据计算频散能量图，从上述方法得到的频散能量图中能够有效分离出不同速度的面波，即基阶面波频散曲线和高阶面波频散曲线。另外，采用矢量波数变换算法，预设台站无须按一定规则进行摆放，可以任意摆放，提高了面波勘探场地的适应性。
应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。
图2是本申请实施例提供的面波勘探装置的示意图，为了便于说明，仅示出与本申请实施例相关的部分。
图2所示的面波勘探装置可以是内置于现有的终端设备内的软件单元、硬件单元、或软硬结合的单元，也可以作为独立的挂件集成到所述终端设备中，还可以作为独立的终端设备存在。
所述面波勘探装置3包括：
采集单元31，用于通过预设台站处的检波器采集从震源传播过来的面波数据。
计算单元32，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图。
提取单元33，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线。
反演单元34，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
可选的，所述计算单元32包括：
提取子单元，用于从所述面波数据中提取震源时间函数。
第一计算子单元，用于计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
第二计算子单元，用于根据所述震源时间函数和所述格林函数计算得到频散能量图。
可选的，所述第二计算子单元包括：
卷积模块，用于将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
第一变换模块，用于对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
结果模块，用于对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
可选的，所述结果模块包括：
变换子模块，用于在对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图之后，对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
转化子模块，用于基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
计算子模块，用于基于所述最终计算式，计算得到频散能量图。
其中，所述时间域的垂直分量数据为：
uz(x,t)＝f0(t)*gzz(x,t)
其中，uz(x,t)表示所述垂直分量数据，f0表示所述震源时间函数，gzz表示所述格林函数。
所述频率域的垂直分量数据为：
U(r,ω)＝F0(ω)G(r,ω)
其中，U为所述频率域的垂直分类数据，
Figure GDA0003467162410000111
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F0为频率域的震源时间函数，G为频率域的格林函数，
Figure GDA0003467162410000112
g(ω,k)是核函数，
Figure GDA0003467162410000113
为波数，J0(kr)是一类零阶贝塞尔函数。
所述中间计算式为：
Figure GDA0003467162410000114
所述最终计算式为：
Figure GDA0003467162410000115
其中，
Figure GDA0003467162410000116
Rdu为下行波的反射系数矩阵，Rud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。
可选的，所述反演单元24包括：
分类子单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类。
建立子单元，用于根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
图3是本申请实施例提供的终端设备的示意图。如图3所示，该实施例的终端设备4包括：处理器40、存储器41以及存储在所述存储器41中并可在所述处理器40上运行的计算机程序42。所述处理器40执行所述计算机程序42时实现上述各个人工源面波勘探方法实施例中的步骤，例如图1所示的步骤S101至S104。或者，所述处理器40执行所述计算机程序42时实现上述各装置实施例中各模块/单元的功能，例如图3所示模块31至34的功能。
示例性的，所述计算机程序42可以被分割成一个或多个模块/单元，所述一个或者多个模块/单元被存储在所述存储器41中，并由所述处理器40执行，以完成本申请。所述一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述所述计算机程序42在所述终端设备4中的执行过程。例如，所述计算机程序42可以被分割成采集单元、计算单元、提取单元、反演单元，各单元具体功能如下：
采集单元，用于通过预设台站处的检波器采集从震源传播过来的面波数据。
计算单元，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图。
提取单元，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线。
反演单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
可选的，所述计算单元包括：
提取子单元，用于从所述面波数据中提取震源时间函数。
第一计算子单元，用于计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
第二计算子单元，用于根据所述震源时间函数和所述格林函数计算得到频散能量图。
可选的，所述第二计算子单元包括：
卷积模块，用于将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
第一变换模块，用于对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
结果模块，用于对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
可选的，所述结果模块包括：
第二变换子模块，用于在对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图之后，对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
转化子模块，用于基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
计算子模块，用于基于所述最终计算式，计算得到频散能量图。
其中，所述时间域的垂直分量数据为：
uz(x,t)＝f0(t)*gzz(x,t)
其中，uz(x,t)表示所述垂直分量数据，f0表示所述震源时间函数，gzz表示所述格林函数。
所述频率域的垂直分量数据为：
U(r,ω)＝F0(ω)G(r,ω)
其中，U为所述频率域的垂直分类数据，
Figure GDA0003467162410000141
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F0为频率域的震源时间函数，G为频率域的格林函数，
Figure GDA0003467162410000142
g(ω,k)是核函数，
Figure GDA0003467162410000143
为波数，J0(kr)是一类零阶贝塞尔函数。
所述中间计算式为：
Figure GDA0003467162410000144
所述最终计算式为：
Figure GDA0003467162410000145
其中，
Figure GDA0003467162410000146
Rdu为下行波的反射系数矩阵，Rud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。
可选的，所述反演单元24包括：
分类子单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类。
建立子单元，用于根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
所述终端设备4可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述终端设备可包括，但不仅限于，处理器40、存储器41。本领域技术人员可以理解，图4仅仅是终端设备4的示例，并不构成对终端设备4的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述终端设备还可以包括输入输出设备、网络接入设备、总线等。
所称处理器40可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
所述存储器41可以是所述终端设备4的内部存储单元，例如终端设备4的硬盘或内存。所述存储器41也可以是所述终端设备4的外部存储设备，例如所述终端设备4上配备的插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)等。进一步地，所述存储器41还可以既包括所述终端设备4的内部存储单元也包括外部存储设备。所述存储器41用于存储所述计算机程序以及所述终端设备所需的其他程序和数据。所述存储器41还可以用于暂时地存储已经输出或者将要输出的数据。
在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述或记载的部分，可以参见其它实施例的相关描述。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
在本申请所提供的实施例中，应该理解到，所揭露的装置/终端设备和方法，可以通过其它的方式实现。例如，以上所描述的装置/终端设备实施例仅仅是示意性的，例如，所述模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性，机械或其它的形式。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
所述集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本申请实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括是电载波信号和电信信号。
以上所述实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。

Claims (8)

1.一种人工源面波勘探方法，其特征在于，包括：
通过预设台站处的检波器采集从震源传播过来的面波数据；
基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图；
从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线；
根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据；
所述基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图，包括：
从所述面波数据中提取震源时间函数；
计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数；
根据所述震源时间函数和所述格林函数计算得到频散能量图。
2.如权利要求1所述的人工源面波勘探方法，其特征在于，所述根据所述震源时间函数和所述格林函数计算得到频散能量图，包括：
将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据；
对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据；
对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
3.如权利要求2所述的人工源面波勘探方法，其特征在于，所述对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图，还包括：
对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式；
基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式；
基于所述最终计算式，计算得到频散能量图。
4.如权利要求3所述的人工源面波勘探方法，其特征在于，所述时间域的垂直分量数据为：
uz(x,t)＝f(t)*gzz(x,t)
其中，uz(x,t)表示所述垂直分量数据，f表示所述震源时间函数，gzz表示所述格林函数；
所述频率域的垂直分量数据为：
U(r,ω)＝F(ω)G(r,ω)
其中，U为所述频率域的垂直分类数据，
Figure FDA0003488196190000021
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F为频率域的震源时间函数，G为频率域的格林函数，
Figure FDA0003488196190000022
g(ω,k)是核函数，
Figure FDA0003488196190000023
为波数，J0(kr)是一类零阶贝塞尔函数；
所述中间计算式为：
Figure FDA0003488196190000024
所述最终计算式为：
Figure FDA0003488196190000025
其中，
Figure FDA0003488196190000026
Rdu为下行波的反射系数矩阵，Rud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。
5.如权利要求1所述的人工源面波勘探方法，其特征在于，所述根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，包括：
根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类；
根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
6.一种面波勘探装置，其特征在于，包括：
采集单元，用于通过预设台站处的检波器采集从震源传播过来的面波数据；
计算单元，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图；
提取单元，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线；
反演单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据；
所述计算单元包括：
提取模块，用于从所述面波数据中提取震源时间函数；
第一计算模块，用于计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数；
第二计算模块，用于根据所述震源时间函数和所述格林函数计算得到频散能量图。
7.一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至5任一项所述方法的步骤。
8.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至5任一项所述方法的步骤。
CN201880000953.4A 2018-08-06 2018-08-06 一种人工源面波勘探方法、面波勘探装置及终端设备 Active CN111164462B (zh)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2018/098980 WO2020029015A1 (zh) 2018-08-06 2018-08-06 一种人工源面波勘探方法、面波勘探装置及终端设备

Publications (2)

Publication Number Publication Date
CN111164462A CN111164462A (zh) 2020-05-15
CN111164462B true CN111164462B (zh) 2022-05-06

Family

ID=69413934

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201880000953.4A Active CN111164462B (zh) 2018-08-06 2018-08-06 一种人工源面波勘探方法、面波勘探装置及终端设备

Country Status (3)

Country Link
US (1) US11137511B2 (zh)
CN (1) CN111164462B (zh)
WO (1) WO2020029015A1 (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111723491B (zh) * 2020-06-29 2024-03-19 重庆大学 一种基于任意分层土壤格林函数的接地参数获取方法
CN112505749B (zh) * 2020-10-19 2024-04-26 中国地质调查局南京地质调查中心(华东地质科技创新中心) 一种基于线形台阵多次覆盖的微动数据采集方法
CN112307422A (zh) * 2020-10-30 2021-02-02 天津光电通信技术有限公司 一种低信噪比下信号时频分析方法、装置及设备
CN112379403B (zh) * 2020-12-14 2024-01-16 北京华晖探测科技股份有限公司 一种地下采空区的探测方法及系统
CN112861721B (zh) * 2021-02-09 2024-05-07 南方科技大学 一种自动提取背景噪声频散曲线的方法及装置
CN113189641B (zh) * 2021-03-25 2024-01-19 西安石油大学 一种两道多模式瑞利波地下探测系统及方法
CN113126146A (zh) * 2021-04-15 2021-07-16 北京市水利规划设计研究院 针对峡谷地区复杂地质的探测方法、控制装置和存储介质
CN113642232B (zh) * 2021-07-22 2024-01-12 南方科技大学 一种面波智能反演勘探方法、存储介质及终端设备
CN114114401B (zh) * 2021-12-01 2024-04-23 北京华晖探测科技股份有限公司 利用轴对称探头在地面激发的sh波进行浅层勘探的方法
CN115270451B (zh) * 2022-07-22 2023-06-06 鹏城实验室 一种获取多介质模型的漏能振型的方法、介质及终端
CN115755176B (zh) * 2022-11-22 2023-06-13 南方科技大学 利用频率汉克尔变换分离波场的面波勘探方法及相关装置
CN116660974B (zh) * 2023-04-11 2024-02-23 中国地震局地球物理研究所 一种基于结构耦合约束的体波和面波三维联合反演方法
CN117518265A (zh) * 2023-11-01 2024-02-06 大连理工大学 层状介质场地下基于频散特性的多模态面波波场反演方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1306621A (zh) * 1998-05-20 2001-08-01 施鲁博格控股有限公司 自适应地震噪声和干扰衰减方法
AU2014221254A1 (en) * 2007-06-29 2014-09-25 Geco Technology B.V. Estimating and Using Slowness Vector Attributes in Connection with a Multi-Component Seismic Gather
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN106950599A (zh) * 2017-05-08 2017-07-14 北京瑞威工程检测有限公司 一种隧道基底密实性检测系统、检测方法及存储介质
CN107561589A (zh) * 2017-10-25 2018-01-09 中国石油化工股份有限公司 一种近地表横波层速度模型建立方法
CN108064348A (zh) * 2017-10-12 2018-05-22 南方科技大学 一种基于两点射线追踪的地震走时层析反演方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2503035C2 (ru) * 2008-08-11 2013-12-27 Эксонмобил Апстрим Рисерч Компани Оценивание свойств почвы с использованием волновых сигналов сейсмических поверхностных волн
US20140236487A1 (en) * 2013-02-21 2014-08-21 Westerngeco L.L.C. Methods and computing systems for processing seismic data
CN104730579B (zh) * 2013-12-18 2018-02-13 中国石油化工股份有限公司 一种基于表层横波速度反演的纵横波联合静校正方法
US10338242B2 (en) * 2015-05-20 2019-07-02 Conocophillips Company Surface wave tomography using sparse data acquisition
US10605937B2 (en) * 2016-05-26 2020-03-31 Cgg Services Sas Device and method for smart picking surface waves dispersion curves
US20180100947A1 (en) 2016-10-06 2018-04-12 The Curators Of The University Of Missouri Spectral analysis of surface waves to detect subsurface voids
CN106772575A (zh) 2016-11-28 2017-05-31 安徽理工大学 一种基于折射波与面波联合反演剩余煤层厚度的方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1306621A (zh) * 1998-05-20 2001-08-01 施鲁博格控股有限公司 自适应地震噪声和干扰衰减方法
AU2014221254A1 (en) * 2007-06-29 2014-09-25 Geco Technology B.V. Estimating and Using Slowness Vector Attributes in Connection with a Multi-Component Seismic Gather
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN106950599A (zh) * 2017-05-08 2017-07-14 北京瑞威工程检测有限公司 一种隧道基底密实性检测系统、检测方法及存储介质
CN108064348A (zh) * 2017-10-12 2018-05-22 南方科技大学 一种基于两点射线追踪的地震走时层析反演方法
CN107561589A (zh) * 2017-10-25 2018-01-09 中国石油化工股份有限公司 一种近地表横波层速度模型建立方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基阶与高阶瑞利波联合反演研究;罗银河 等;《地球物理学报》;20080131;第51卷(第1期);第242-249页 *
近岸带三组成波耦合作用的观测与分析;陈子燊 等;《热带海洋学报》;20031130;第22卷(第6期);第46-53页 *

Also Published As

Publication number Publication date
US11137511B2 (en) 2021-10-05
WO2020029015A1 (zh) 2020-02-13
CN111164462A (zh) 2020-05-15
US20210141113A1 (en) 2021-05-13

Similar Documents

Publication Publication Date Title
CN111164462B (zh) 一种人工源面波勘探方法、面波勘探装置及终端设备
CN109923440B (zh) 面波勘探方法及终端设备
CN103926623B (zh) 一种压制逆时偏移低频噪音的方法
CN104678435A (zh) 一种提取Rayleigh面波频散曲线的方法
KR20120019055A (ko) 탄성 매질에서의 주파수 영역 역시간 구조보정을 이용한 지하구조의 영상화 장치 및 방법
CN104375188A (zh) 一种地震波透射衰减补偿方法及装置
US20150309200A1 (en) A method for processing acoustic waveforms
Delph et al. Constraining crustal properties using receiver functions and the autocorrelation of earthquake‐generated body waves
Petrovic et al. Joint deconvolution of building and downhole strong‐motion recordings: Evidence for the seismic wavefield being radiated back into the shallow geological layers
CN105629299A (zh) 角度域叠前深度偏移的走时、角度表获取方法及成像方法
CN111458749B (zh) 应用于被动源地震勘探的面波与体波的分离方法及系统
CN114428324A (zh) 叠前高角度快速傅里叶变换地震成像方法、系统、设备
Ullah et al. Soft sediment thickness and shear-wave velocity estimation from the H/V technique up to the bedrock at meteorite impact crater site, Sao Paulo city, Brazil
CN113534243B (zh) 一种被动源Marchenko成像方法及系统
CN106257309A (zh) 叠后地震数据体处理方法及装置
WO2021174434A1 (zh) 一种震电波场联合提取瑞雷波频散特征的面波勘探方法
CN113126161B (zh) 一种基于冲击波预测溶洞埋深和尺寸的方法及系统
Suzuki et al. New method to estimate bedrock shape of small-scale basin using modal properties of sediment
CN111123361A (zh) 垂直地震剖面地震数据规则化重建方法及装置、存储介质
CN115170428A (zh) 一种声波远探测成像图的降噪方法
CN114200522B (zh) 深度域地震子波提取方法和装置、存储介质及电子设备
CN113721296A (zh) 远震数据处理方法及装置
CN112649848B (zh) 利用波动方程求解地震波阻抗的方法和装置
US20170336523A1 (en) Seismic signal processing method, apparatus and system
Kalil et al. The efficiency of horizontal to vertical spectral ratio technique for buried monuments delineation (case study, Saqqara (Zoser) pyramid, Egypt)

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
EE01 Entry into force of recordation of patent licensing contract
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20200515

Assignee: SHANDONG University

Assignor: SOUTH University OF SCIENCE AND TECHNOLOGY OF CHINA

Contract record no.: X2022980022879

Denomination of invention: An artificial source surface wave exploration method, surface wave exploration device and terminal equipment

Granted publication date: 20220506

License type: Common License

Record date: 20221122