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

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

Info

Publication number
WO2020029015A1
WO2020029015A1 PCT/CN2018/098980 CN2018098980W WO2020029015A1 WO 2020029015 A1 WO2020029015 A1 WO 2020029015A1 CN 2018098980 W CN2018098980 W CN 2018098980W WO 2020029015 A1 WO2020029015 A1 WO 2020029015A1
Authority
WO
WIPO (PCT)
Prior art keywords
surface wave
dispersion curve
dispersion
function
data
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.)
Ceased
Application number
PCT/CN2018/098980
Other languages
English (en)
French (fr)
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.)
Southern University of Science and Technology
Original Assignee
Southern 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 Southern University of Science and Technology filed Critical Southern University of Science and Technology
Priority to US17/051,247 priority Critical patent/US11137511B2/en
Priority to PCT/CN2018/098980 priority patent/WO2020029015A1/zh
Priority to CN201880000953.4A priority patent/CN111164462B/zh
Publication of WO2020029015A1 publication Critical patent/WO2020029015A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

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

Definitions

  • FIG. 3 is a schematic diagram of a surface wave exploration device according to an embodiment of the present application.
  • the term “if” can be construed as “when” or “once” or “in response to a determination” or “in response to a detection” depending on the context .
  • the phrase “if determined” or “if [the described condition or event] is detected” can be interpreted, depending on the context, to mean “once determined” or “in response to the determination” or “once [the condition or event described ] “Or” In response to [Description of condition or event] detected ".
  • the surface wave is a kind of seismic wave, which mainly propagates on the ground surface and has the largest energy.
  • the wave velocity is about 3.8 km / s, which is lower than the body wave and is often recorded last.
  • Surface waves are actually secondary waves derived from body waves on the surface. The propagation of surface waves is more complicated, which can cause both the ups and downs of the ground surface and the ground surface to perform lateral shearing, among which the shearing motion destroys the building most strongly.
  • Surface waves include Rayleigh waves, Love waves, hydraulic waves, Stone Geb waves and so on.
  • studying low-frequency Rayleigh wave dispersion in natural seismic waves can solve deep geological structural problems; studying higher-frequency Rayleigh waves excited by artificial sources can solve shallow geology such as engineering survey, site and foundation treatment evaluation, obstacle and cavity detection, etc. problem. Therefore, Rayleigh waves are preferably used in the embodiments of the present application.
  • step S102 a dispersion energy map is calculated based on the vector wave number transformation algorithm and calculated according to the surface wave data.
  • a dispersion curve is extracted from the dispersion energy map, and the dispersion curve includes a fundamental-order surface wave dispersion curve and a higher-order surface wave dispersion curve.
  • the actual received background noise data consists of waves generated by various vibrations, including not only surface waves but also body waves.
  • the surface wave will be dispersed in a non-uniform medium, that is, the surface wave is composed of modes with different phase velocities.
  • a large number of studies have proved that the higher-order surface wave part of the dispersion curve plays a key role in the analysis of formation structure inversion.
  • the calculated dispersion energy map can effectively separate surface waves with different velocities (i.e., capable of separating fundamental-order and high-order surface waves) and body waves. Components.
  • the Rayleigh wave energy in the frequency range interval corresponding to the buried depth of the interlayer steps from the fundamental order to the first-order or higher-order modal steps, resulting in the appearance of the fundamental-order and higher-order surface wave dispersion curves.
  • the imaging quality may be worse in actual data. Therefore, in the dispersive energy map, the distribution of each modal energy is also closely related to the stratigraphic structure.
  • the frequency intervals are classified according to the energy distribution relationship of each modal in different frequency intervals in the dispersion energy map, so as to quickly establish a simple layered stratum model as an initial model for subsequent accurate inversion .
  • an inversion algorithm such as a simulated annealing algorithm or a genetic algorithm may be used to invert the dispersion curve to obtain formation information and / or velocity information of the vibration wave.
  • formation depth information and velocity profile can be obtained to achieve the exploration of the stratum structure; stratum information can include formation depth, S-wave velocity, density, and P-wave velocity.
  • the step of establishing an initial stratum model based on the fundamental-order surface wave dispersion curve and the higher-order surface wave dispersion curve includes:
  • the initial stratigraphic model is established according to the correspondence between the classified frequency interval and the stratum.
  • Figure 5a is the FC dispersion energy diagram, and the dotted line in the figure is the theoretical surface wave dispersion curve;
  • Figure 5b) is the discrete dispersion points obtained based on the dispersion energy, and according to the distribution characteristics, The frequency interval is divided into 4 categories;
  • Figure 5c) is the theoretical FC dispersion energy map obtained by the Green's function kernel function.
  • the surface wave data transmitted from the seismic source is collected by a detector at a preset station, and the collecting device can be arbitrarily arranged, which reduces the requirements for the layout site and improves the adaptability of the surface wave exploration site.
  • the initial stratigraphic model is established by using the extracted dispersion curve, which reduces the calculation time of the inversion operation. Then, based on the initial stratigraphic model, joint inversion of the base-order surface wave dispersion curve and the higher-order surface wave dispersion curve are performed.
  • the order surface wave dispersion information is added to the formation inversion operation, which reduces the uncertainty of the inversion operation and improves the robustness of the inversion operation.
  • FIG. 2 is a schematic diagram of an implementation flow of calculating a dispersion energy map in an artificial source surface wave exploration method according to an embodiment of the present application.
  • step S102 in the embodiment shown in FIG. 1 may include the following steps:
  • step S201 a seismic source time function is extracted from the surface wave data.
  • Step S202 Calculate a gun offset between the seismic source and the geophone, and calculate a Green function between the preset station and the seismic source according to the gun offset.
  • Step S203 Calculate a dispersion energy map according to the source time function and the Green function.
  • step S203 may further include:
  • Convolution processing is performed on the source time function and the Green function to obtain vertical component data of the surface wave data at the preset station in the time domain.
  • U is the vertical classification data in the frequency domain
  • ⁇ is the angular frequency
  • ⁇ 2 ⁇ f
  • f is the frequency
  • F 0 is the source time function in the frequency domain
  • G is the Green's function in the frequency domain.
  • g ( ⁇ , k) is the kernel function
  • J 0 (kr) is a kind of zero-order Bessel function.
  • the performing a vector wave number transformation on the vertical component data in the frequency domain to obtain a dispersion energy map further includes:
  • the intermediate calculation formula is converted into a final calculation formula.
  • A F 0 ( ⁇ ) is a constant.
  • the vector wave number transformation is performed on the vertical component data U in the frequency domain to obtain
  • the kernel function g ( ⁇ , k) has the following characteristics: the value of g ( ⁇ , k) is inversely proportional to the value of the duration function that determines the surface wave dispersion characteristics, that is,
  • VWTM Vector Wavenumber Transform Method
  • the dispersion energy map can be calculated.
  • An embodiment of the present application proposes a vector wave number transformation algorithm, and based on this vector wave number transformation algorithm, calculates a dispersion energy map from surface wave data.
  • the dispersion energy map obtained from the above method can effectively separate surface waves of different speeds. , That is, the fundamental-order surface wave dispersion curve and the higher-order surface wave dispersion curve.
  • the preset stations need not be placed according to certain rules, and can be placed arbitrarily, which improves the adaptability of surface wave exploration sites.
  • FIG. 3 is a schematic diagram of a surface wave exploration device provided in an embodiment of the present application. For convenience of explanation, only a part related to the embodiment of the present application is shown.
  • the surface wave exploration device shown in FIG. 3 may be a software unit, a hardware unit, or a combination of hardware and software that is built into an existing terminal device, or may be integrated into the terminal device as an independent pendant, or may be used as A separate terminal device exists.
  • the surface wave exploration device 3 includes:
  • the collecting unit 31 is configured to collect surface wave data transmitted from the seismic source through a detector at a preset station.
  • the calculation unit 32 is configured to calculate a dispersion energy map based on the vector wave number transformation algorithm and calculate according to the surface wave data.
  • the extraction unit 33 is configured to extract a dispersion curve from the dispersion energy map, where the dispersion curve includes a fundamental-order surface wave dispersion curve and a higher-order surface wave dispersion curve.
  • An inversion unit 34 is configured to establish an initial stratigraphic model according to the fundamental-order surface wave dispersion curve and the higher-order surface wave dispersion curve, and to analyze the fundamental-order surface wave dispersion curve and The high-order surface wave dispersion curve is subjected to joint inversion to obtain inversion data of the stratum structure.
  • the calculation unit 32 includes:
  • An extraction subunit is used to extract a seismic source time function from the surface wave data.
  • a first calculation subunit is configured to calculate a gun offset between the seismic source and the geophone, and calculate a Green's function between the preset station and the seismic source according to the gun offset.
  • a second calculation subunit is configured to calculate a dispersion energy map according to the source time function and the Green function.
  • the second calculation subunit includes:
  • a convolution module is configured to perform convolution processing on the source time function and the Green function to obtain vertical component data of the surface wave data at the preset station in the time domain.
  • the first transformation module is configured to perform Fourier transform on the vertical component data in the time domain to obtain vertical component data in the frequency domain.
  • a result module configured to perform vector wave number transformation on the vertical component data in the frequency domain to obtain a dispersion energy map.
  • the result module includes:
  • a transformation sub-module is configured to perform vector wave number transformation on the vertical component data in the frequency domain to obtain an intermediate calculation formula after performing vector wave number transformation on the vertical component data in the frequency domain.
  • a transformation submodule is used to transform the intermediate calculation formula into a final calculation formula based on the orthogonal property of the Bezier function.
  • a calculation sub-module is configured to calculate a dispersion energy map based on the final calculation formula.
  • the vertical component data in the time domain is:
  • u z (x, t) represents the vertical component data
  • f 0 represents the source time function
  • g zz represents the Green's function
  • the vertical component data in the frequency domain is:
  • U vertical classification data in the frequency domain
  • I the distance between two observation stations
  • ⁇ the angular frequency
  • ⁇ 2 ⁇ f
  • f the frequency
  • F 0 the source time function in the frequency domain
  • G the Green's function in the frequency domain.
  • g ( ⁇ , k) is the kernel function
  • J 0 (kr) is a kind of zero-order Bessel function.
  • the intermediate calculation formula is:
  • R du is the reflection coefficient matrix of the downstream wave
  • Rud is the reflection coefficient matrix of the upstream wave
  • I is the identity matrix
  • det ⁇ is the matrix determinant.
  • the inversion unit 24 includes:
  • a classification subunit is configured to classify the frequency interval according to an energy distribution of a surface wave mode in each frequency interval in the basic-order surface wave dispersion curve and the higher-order surface wave dispersion curve.
  • a subunit is established, and is configured to establish the initial stratum model according to the correspondence between the classified frequency interval and the stratum.
  • FIG. 4 is a schematic diagram of a terminal device according to an embodiment of the present application.
  • the terminal device 4 of this embodiment includes a processor 40, a memory 41, and a computer program 42 stored in the memory 41 and executable on the processor 40.
  • the processor 40 executes the computer program 42, the steps in the embodiments of the artificial source surface wave exploration method described above are implemented, for example, steps S101 to S104 shown in FIG. 1.
  • the processor 40 executes the computer program 42
  • the functions of the modules / units in the foregoing device embodiments are implemented, for example, the functions of modules 31 to 34 shown in FIG. 3.
  • the computer program 42 may be divided into one or more modules / units, and the one or more modules / units are stored in the memory 41 and executed by the processor 40 to complete This application.
  • the one or more modules / units may be a series of computer program instruction segments capable of performing specific functions, and the instruction segments are used to describe an execution process of the computer program 42 in the terminal device 4.
  • the computer program 42 can be divided into an acquisition unit, a calculation unit, an extraction unit, and an inversion unit. The specific functions of each unit are as follows:
  • the collecting unit 31 is configured to collect surface wave data transmitted from the seismic source through a detector at a preset station.
  • the calculation unit 32 is configured to calculate a dispersion energy map based on the vector wave number transformation algorithm and calculate according to the surface wave data.
  • the extraction unit 33 is configured to extract a dispersion curve from the dispersion energy map, where the dispersion curve includes a fundamental-order surface wave dispersion curve and a higher-order surface wave dispersion curve.
  • An inversion unit 34 is configured to establish an initial stratigraphic model according to the fundamental-order surface wave dispersion curve and the higher-order surface wave dispersion curve, and to analyze the fundamental-order surface wave dispersion curve and The high-order surface wave dispersion curve is subjected to joint inversion to obtain inversion data of the stratum structure.
  • the calculation unit 32 includes:
  • An extraction subunit is used to extract a seismic source time function from the surface wave data.
  • a first calculation subunit is configured to calculate a gun offset between the seismic source and the geophone, and calculate a Green's function between the preset station and the seismic source according to the gun offset.
  • a second calculation subunit is configured to calculate a dispersion energy map according to the source time function and the Green function.
  • the second calculation subunit includes:
  • a convolution module is configured to perform convolution processing on the source time function and the Green function to obtain vertical component data of the surface wave data at the preset station in the time domain.
  • the first transformation module is configured to perform Fourier transform on the vertical component data in the time domain to obtain vertical component data in the frequency domain.
  • a result module configured to perform vector wave number transformation on the vertical component data in the frequency domain to obtain a dispersion energy map.
  • the result module includes:
  • a second transforming sub-module is configured to perform vector wave number transformation on the vertical component data in the frequency domain to obtain an intermediate calculation formula after performing vector wave number transformation on the vertical component data in the frequency domain.
  • a transformation submodule is used to transform the intermediate calculation formula into a final calculation formula based on the orthogonal property of the Bezier function.
  • a calculation sub-module is configured to calculate a dispersion energy map based on the final calculation formula.
  • the vertical component data in the time domain is:
  • u z (x, t) represents the vertical component data
  • f 0 represents the source time function
  • g zz represents the Green's function
  • the vertical component data in the frequency domain is:
  • U vertical classification data in the frequency domain
  • I the distance between two observation stations
  • ⁇ the angular frequency
  • ⁇ 2 ⁇ f
  • f the frequency
  • F 0 the source time function in the frequency domain
  • G the Green's function in the frequency domain.
  • g ( ⁇ , k) is the kernel function
  • J 0 (kr) is a kind of zero-order Bessel function.
  • the intermediate calculation formula is:
  • R du is the reflection coefficient matrix of the downstream wave
  • Rud is the reflection coefficient matrix of the upstream wave
  • I is the identity matrix
  • det ⁇ is the matrix determinant.
  • the inversion unit 24 includes:
  • a classification subunit is configured to classify the frequency interval according to an energy distribution of a surface wave mode in each frequency interval in the basic-order surface wave dispersion curve and the higher-order surface wave dispersion curve.
  • a subunit is established, and is configured to establish the initial stratum model according to the correspondence between the classified frequency interval and the stratum.
  • the terminal device 4 may be a computing device such as a desktop computer, a notebook, a palmtop computer, and a cloud server.
  • the terminal device may include, but is not limited to, a processor 40 and a memory 41.
  • FIG. 4 is only an example of the terminal device 4 and does not constitute a limitation on the terminal device 4. It may include more or fewer components than shown in the figure, or combine some components or different components
  • the terminal device may further include an input / output device, a network access device, a bus, and the like.
  • the so-called processor 40 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application specific integrated circuits (ASICs), Ready-made programmable gate array (Field-Programmable Gate Array, FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc.
  • a general-purpose processor may be a microprocessor or the processor may be any conventional processor or the like.
  • the memory 41 may be an internal storage unit of the terminal device 4, such as a hard disk or a memory of the terminal device 4.
  • the memory 41 may also be an external storage device of the terminal device 4, such as a plug-in hard disk, a Smart Media Card (SMC), and a Secure Digital (SD) provided on the terminal device 4. Card, flash card, etc.
  • the memory 41 may include both an internal storage unit of the terminal device 4 and an external storage device.
  • the memory 41 is configured to store the computer program and other programs and data required by the terminal device.
  • the memory 41 may also be used to temporarily store data that has been or will be output.
  • the disclosed apparatus / terminal device and method may be implemented in other ways.
  • the device / terminal device embodiments described above are only schematic.
  • the division of the modules or units is only a logical function division.
  • components can be combined or integrated into another system, or some features can be ignored or not implemented.
  • the displayed or discussed mutual coupling or direct coupling or communication connection may be indirect coupling or communication connection through some interfaces, devices or units, which may be electrical, mechanical or other forms.
  • the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, may be located in one place, or may be distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the objective of the solution of this embodiment.
  • the functional units in the embodiments of the present application may be integrated into one processing unit, or each of the units may exist separately physically, or two or more units may be integrated into one unit.
  • the above integrated unit may be implemented in the form of hardware or in the form of software functional unit.
  • the integrated module / unit When the integrated module / unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, this application implements all or part of the processes in the method of the above embodiment, and can also be completed by a computer program instructing related hardware.
  • the computer program can be stored in a computer-readable storage medium.
  • the computer When the program is executed by a processor, the steps of the foregoing method embodiments can be implemented.
  • the computer program includes computer program code, and the computer program code may be in a source code form, an object code form, an executable file, or some intermediate form.
  • the computer-readable medium may include: any entity or device capable of carrying the computer program code, a recording medium, a U disk, a mobile hard disk, a magnetic disk, an optical disk, a computer memory, a read-only memory (ROM, Read-Only Memory) , Random Access Memory (RAM, Random Access Memory), electric carrier signals, telecommunication signals, and software distribution media.
  • ROM Read-Only Memory
  • RAM Random Access Memory
  • electric carrier signals telecommunication signals
  • software distribution media any entity or device capable of carrying the computer program code
  • a recording medium a U disk, a mobile hard disk, a magnetic disk, an optical disk, a computer memory, a read-only memory (ROM, Read-Only Memory) , Random Access Memory (RAM, Random Access Memory), electric carrier signals, telecommunication signals, and software distribution media.

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是本申请实施例提供的终端设备的示意图；
图5是本申请实施例提供的提取到的F-C频散谱(a)、频率区间分类(b)和格林函数核函数得到的理论F-C频散谱(c)的示意图；
图6是本申请实施例提供的深度域的频散谱和地层模型的示意图。
本发明的实施方式
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本申请实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本申请。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本申请的描述。
应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当理解，在此本申请说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本申请。如在本申请说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
还应当进一步理解，在本申请说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
如在本说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
为了说明本申请所述的技术方案，下面通过具体实施例来进行说明。
图1是本申请实施例提供的人工源面波勘探方法的实现流程示意图，如图所示，所述方法可以包括以下步骤：
步骤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波速度等。
在本申请实施例中，所述根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，包括：
根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类。
根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
如图5所示，图5a)为F-C频散能量图，图中点线为理论面波频散曲线；图5b)为根据频散能量提取得到的离散的频散点，并根据分布特性将频率区间分为4类；图5c)为格林函数核函数得到的理论F-C频散能量图。
将频率-速度域的频散点，根据半波长理论转换到深度-速度域，见图6。图6a)深度域的频散能量图，图6b)为地层模型。可以看到，频散能量图上1号和3号点线上的点为基阶频散曲线上的点，2号和4号点线上的点为高阶频散曲线上的点。将地层模型与得到的深度域的频散曲线进行对比，可以看到位移20-40m埋深的第三层(低速层)与深度-速度剖面中的4号频散点的分布基本一致；位移0-10m的第一层与深度-速度剖面中的2号频散点的分布也基本一致。这样我们可以看到高阶频散曲线的分布与地层确实存在一一对应的关系。由此证明，通过在频率域将频散点进行分类，对地层进行分层，建立初始建模的思路是正确的。
本申请实施例通过预设台站处的检波器采集从震源传播过来的面波数据，可以任意布置采集装置，降低了对布置场地的要求，提高了面波勘探场地的适应性。通过利用提取出的频散曲线建立初始地层模型，降低了反演运算的计算时间；之后根据初始地层模型对基阶面波频散曲线和高阶面波频散曲线进行联合反演，将高阶面波频散的信息加入到地层的反演运算中，从而降低了反演运算的不确定性，提高了反演运算的鲁棒性。
图2是本申请实施例提供的人工源面波勘探方法中计算频散能量图的实现流程示意图。如图所示，图1所示实施例中步骤S102可以包括以下步骤：
步骤S201，从所述面波数据中提取震源时间函数。
在实际应用中，可以有多种方法提取震源时间函数。例如，可以直接将震源时间函数近似为雷克子波；还可以从多道地震记录中根据相关性提取震源时间函数。只要能够提取出震源时间函数即可，不对提取方法做具体限定。
步骤S202，计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
其中，炮检距为震源与检波器之间的距离。因为检波器安装于预设台站处，所以炮检距也可看作是预设台站与震源之间的距离。在计算式中，可以用g zz表示预设台站与震源之间的格林函数。
步骤S203，根据所述震源时间函数和所述格林函数计算得到频散能量图。
在本申请实施例中，步骤S203还可以包括：
将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
其中，所述时间域的垂直分量数据为：
u z(x,t)＝f 0(t)*g zz(x,t)
其中，u z(x,t)表示所述垂直分量数据，f 0表示所述震源时间函数，g zz表示所述格林函数。
对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
其中，所述频率域的垂直分量数据为：
U(r,ω)＝F 0(ω)G(r,ω)
式中，U为所述频率域的垂直分类数据，
Figure PCTCN2018098980-appb-000001
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F 0为频率域的震源时间函数，G为频率域的格林函数，
Figure PCTCN2018098980-appb-000002
g(ω,k)是核函数，
Figure PCTCN2018098980-appb-000003
为波数，J 0(kr)是一类零阶贝塞尔函数。
对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
在本申请实施例中，所述对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图，还包括：
对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
基于所述最终计算式，计算得到频散能量图。
在实际应用中，当震源时间函数为雷克子波时，频率域中F 0(ω)为纯实数函数，即可得预设台站接收到的面波数据的谱函数近似为格林函数的虚部，二者仅在幅值上有差异：
U(r,ω)＝A·{G(r,ω)}
其中，A＝F 0(ω)为常量。对频率域的垂直分量数据U进行矢量波数变换，可得
Figure PCTCN2018098980-appb-000004
将
Figure PCTCN2018098980-appb-000005
带入上式，可得中间计算式为：
Figure PCTCN2018098980-appb-000006
根据贝塞尔函数的正交性
Figure PCTCN2018098980-appb-000007
可将上式化简得到最终计算式为：
Figure PCTCN2018098980-appb-000008
这里核函数g(ω,k)具有如下特点：g(ω,k)的值反比于确定面波频散特性的久期函数值，即
Figure PCTCN2018098980-appb-000009
其中，
Figure PCTCN2018098980-appb-000010
R du为下行波的反射系数矩阵，R ud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。当k＝k n(ω)(n＝1,2,3,...)是核函数的g(ω,k)的极点时，核函数的值趋于无穷大。利用这一性质，提出了矢量波数变换法(Vector Wavenumber Transform Method，VWTM)提取频散曲线。
基于上述最终计算式，即可计算得到频散能量图。
本申请实施例提出了一种矢量波数变换算法，并基于这种矢量波数变换算法根据面波数据计算频散能量图，从上述方法得到的频散能量图中能够有效分离出不同速度的面波，即基阶面波频散曲线和高阶面波频散曲线。另外，采用矢量波数变换算法，预设台站无须按一定规则进行摆放，可以任意摆放，提高了面波勘探场地的适应性。
应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。
图3是本申请实施例提供的面波勘探装置的示意图，为了便于说明，仅示出与本申请 实施例相关的部分。
图3所示的面波勘探装置可以是内置于现有的终端设备内的软件单元、硬件单元、或软硬结合的单元，也可以作为独立的挂件集成到所述终端设备中，还可以作为独立的终端设备存在。
所述面波勘探装置3包括：
采集单元31，用于通过预设台站处的检波器采集从震源传播过来的面波数据。
计算单元32，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图。
提取单元33，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线。
反演单元34，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
可选的，所述计算单元32包括：
提取子单元，用于从所述面波数据中提取震源时间函数。
第一计算子单元，用于计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
第二计算子单元，用于根据所述震源时间函数和所述格林函数计算得到频散能量图。
可选的，所述第二计算子单元包括：
卷积模块，用于将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
第一变换模块，用于对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
结果模块，用于对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
可选的，所述结果模块包括：
变换子模块，用于在对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图之后，对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
转化子模块，用于基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
计算子模块，用于基于所述最终计算式，计算得到频散能量图。
其中，所述时间域的垂直分量数据为：
u z(x,t)＝f 0(t)*g zz(x,t)
其中，u z(x,t)表示所述垂直分量数据，f 0表示所述震源时间函数，g zz表示所述格林函 数。
所述频率域的垂直分量数据为：
U(r,ω)＝F 0(ω)G(r,ω)
其中，U为所述频率域的垂直分类数据，
Figure PCTCN2018098980-appb-000011
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F 0为频率域的震源时间函数，G为频率域的格林函数，
Figure PCTCN2018098980-appb-000012
g(ω,k)是核函数，
Figure PCTCN2018098980-appb-000013
为波数，J 0(kr)是一类零阶贝塞尔函数。
所述中间计算式为：
Figure PCTCN2018098980-appb-000014
所述最终计算式为：
Figure PCTCN2018098980-appb-000015
其中，
Figure PCTCN2018098980-appb-000016
R du为下行波的反射系数矩阵，R ud为上行波的反射系数矩阵，I为单位矩阵，det{}为矩阵行列式。
可选的，所述反演单元24包括：
分类子单元，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线中各个频率区间的面波模态的能量分布，对所述频率区间进行分类。
建立子单元，用于根据分类后的频率区间与地层的对应关系建立所述初始地层模型。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
图4是本申请实施例提供的终端设备的示意图。如图4所示，该实施例的终端设备4包括：处理器40、存储器41以及存储在所述存储器41中并可在所述处理器40上运行的计算机程序42。所述处理器40执行所述计算机程序42时实现上述各个人工源面波勘探方法实施例中的步骤，例如图1所示的步骤S101至S104。或者，所述处理器40执行所述计算机程序42时实现上述各装置实施例中各模块/单元的功能，例如图3所示模块31至34的功能。
示例性的，所述计算机程序42可以被分割成一个或多个模块/单元，所述一个或者多个模块/单元被存储在所述存储器41中，并由所述处理器40执行，以完成本申请。所述一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述所述计算机程序42在所述终端设备4中的执行过程。例如，所述计算机程序42可以被分割成采集单元、计算单元、提取单元、反演单元，各单元具体功能如下：
采集单元31，用于通过预设台站处的检波器采集从震源传播过来的面波数据。
计算单元32，用于基于矢量波数变换算法，并根据所述面波数据计算得到频散能量图。
提取单元33，用于从所述频散能量图中提取频散曲线，所述频散曲线包括基阶面波频散曲线和高阶面波频散曲线。
反演单元34，用于根据所述基阶面波频散曲线和所述高阶面波频散曲线建立初始地层模型，并根据所述初始地层模型对所述基阶面波频散曲线和所述高阶面波频散曲线进行联合反演，得到地层结构的反演数据。
可选的，所述计算单元32包括：
提取子单元，用于从所述面波数据中提取震源时间函数。
第一计算子单元，用于计算所述震源与所述检波器之间的炮检距，并根据所述炮检距计算所述预设台站与所述震源之间的格林函数。
第二计算子单元，用于根据所述震源时间函数和所述格林函数计算得到频散能量图。
可选的，所述第二计算子单元包括：
卷积模块，用于将所述震源时间函数与所述格林函数进行卷积处理，得到所述预设台站处的面波数据在时间域的垂直分量数据。
第一变换模块，用于对所述时间域的垂直分量数据进行傅里叶变换得到频率域的垂直分量数据。
结果模块，用于对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图。
可选的，所述结果模块包括：
第二变换子模块，用于在对所述频率域的垂直分量数据进行矢量波数变换得到频散能量图之后，对所述频率域的垂直分量数据进行矢量波数变换得到中间计算式。
转化子模块，用于基于贝塞尔函数的正交性质，将所述中间计算式转化为最终计算式。
计算子模块，用于基于所述最终计算式，计算得到频散能量图。
其中，所述时间域的垂直分量数据为：
u z(x,t)＝f 0(t)*g zz(x,t)
其中，u z(x,t)表示所述垂直分量数据，f 0表示所述震源时间函数，g zz表示所述格林函数。
所述频率域的垂直分量数据为：
U(r,ω)＝F 0(ω)G(r,ω)
其中，U为所述频率域的垂直分类数据，
Figure PCTCN2018098980-appb-000017
为两个观测台站之间的距离，ω为角频率，ω＝2πf，f为频率，F 0为频率域的震源时间函数，G为频率域的格林函数，
Figure PCTCN2018098980-appb-000018
g(ω,k)是核函数，
Figure PCTCN2018098980-appb-000019
为波数，J 0(kr)是一类零阶贝塞尔函数。
所述中间计算式为：
Figure PCTCN2018098980-appb-000020
所述最终计算式为：
Figure PCTCN2018098980-appb-000021
其中，
Figure PCTCN2018098980-appb-000022
R du为下行波的反射系数矩阵，R ud为上行波的反射系数矩阵，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 (10)

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

Priority Applications (3)

Application Number Priority Date Filing Date Title
US17/051,247 US11137511B2 (en) 2018-08-06 2018-08-06 Active source surface wave prospecting method, surface wave exploration device and computer-readable storage medium
PCT/CN2018/098980 WO2020029015A1 (zh) 2018-08-06 2018-08-06 一种人工源面波勘探方法、面波勘探装置及终端设备
CN201880000953.4A CN111164462B (zh) 2018-08-06 2018-08-06 一种人工源面波勘探方法、面波勘探装置及终端设备

Applications Claiming Priority (1)

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

Publications (1)

Publication Number Publication Date
WO2020029015A1 true WO2020029015A1 (zh) 2020-02-13

Family

ID=69413934

Family Applications (1)

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

Country Status (3)

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

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111723491A (zh) * 2020-06-29 2020-09-29 重庆大学 一种基于任意分层土壤格林函数的接地参数获取方法
CN114994755A (zh) * 2022-05-25 2022-09-02 平安煤炭开采工程技术研究院有限责任公司 一种无损高效的矸石山地质结构调查方法
CN115236738A (zh) * 2022-06-27 2022-10-25 中国地质科学院地球物理地球化学勘查研究所 一种基于汉克尔变换的勒夫波频散提取方法
CN116027411A (zh) * 2021-10-25 2023-04-28 中国石油化工股份有限公司 一种数据驱动的面波预测方法、装置及电子设备和系统
CN119720604A (zh) * 2025-02-27 2025-03-28 华东交通大学 一种非线性超表面瑞利波频散关系计算方法及系统
CN119846714A (zh) * 2025-01-25 2025-04-18 中国矿业大学 一种边坡裂隙被动源震波成像方法
CN120370399A (zh) * 2025-06-26 2025-07-25 湖北省水文地质工程地质勘察院有限公司 基于震源的面波体波综合勘探方法及勘探系统

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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 南方科技大学 一种面波智能反演勘探方法、存储介质及终端设备
CN114201794B (zh) * 2021-11-03 2025-05-13 中国水利水电第八工程局有限公司 一种基于bim的异形土石方工程量测算方法及装置
CN114114401B (zh) * 2021-12-01 2024-04-23 北京华晖探测科技股份有限公司 利用轴对称探头在地面激发的sh波进行浅层勘探的方法
CN116299689A (zh) * 2021-12-21 2023-06-23 中国石油天然气集团有限公司 构建近地表速度模型的方法及装置
CN115903016B (zh) * 2022-07-15 2025-05-16 青岛黄海学院 道路平整度低情况下的车辆震源信号面波频散提取方法
CN115270451B (zh) * 2022-07-22 2023-06-06 鹏城实验室 一种获取多介质模型的漏能振型的方法、介质及终端
CN115407397B (zh) * 2022-09-01 2025-03-18 青岛地质工程勘察院(青岛地质勘查开发局) 一种瑞雷波频散曲线有监督学习反演方法及系统
CN118033728B (zh) * 2022-11-14 2025-09-09 中国石油天然气集团有限公司 基于混合源面波的速度结构模型的构建方法
CN115755176B (zh) * 2022-11-22 2023-06-13 南方科技大学 利用频率汉克尔变换分离波场的面波勘探方法及相关装置
CN115877451A (zh) * 2022-12-09 2023-03-31 中铁第四勘察设计院集团有限公司 一种多激发点瞬态面波勘探方法、设备及存储介质
CN116660974B (zh) * 2023-04-11 2024-02-23 中国地震局地球物理研究所 一种基于结构耦合约束的体波和面波三维联合反演方法
CN116482759B (zh) * 2023-04-25 2025-09-23 上海城勘信息科技有限公司 一种可模态分离的高分辨率频率-波数法
CN117518265B (zh) * 2023-11-01 2024-10-18 大连理工大学 层状介质场地下基于频散特性的多模态面波波场反演方法
CN118068422B (zh) * 2024-02-20 2025-07-29 陕西小保当矿业有限公司 一种体波槽波的波场拟合方法及设备
CN118483743B (zh) * 2024-05-06 2025-05-27 南方科技大学 一种基于应变场提取面波频散曲线的方法
CN118884531B (zh) * 2024-09-26 2025-01-21 鹏城实验室 提取高频面波群速度方法、装置、设备、存储介质及产品
CN119224853B (zh) * 2024-11-28 2025-03-04 中国海洋大学 基于Bi-LSTM网络频散曲线反演方法、介质和设备
CN120044595B (zh) * 2025-02-24 2025-07-22 北京市勘察设计研究院有限公司 一种城市道路塌陷风险评价方法及系统
CN119936995B (zh) * 2025-03-04 2025-12-16 安徽省地质调查院(安徽省地质科学研究所) 一种被动源和主动源瑞雷面波联合勘探方法及系统
CN120468936B (zh) * 2025-05-23 2025-12-09 甘肃省地震局（中国地震局兰州地震研究所） 地震波形分类检测方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN104730579A (zh) * 2013-12-18 2015-06-24 中国石油化工股份有限公司 一种基于表层横波速度反演的纵横波联合静校正方法
WO2016187252A1 (en) * 2015-05-20 2016-11-24 Conocophillips Company Surface wave tomography using sparse data acquisition
CN106772575A (zh) * 2016-11-28 2017-05-31 安徽理工大学 一种基于折射波与面波联合反演剩余煤层厚度的方法
CN106950599A (zh) * 2017-05-08 2017-07-14 北京瑞威工程检测有限公司 一种隧道基底密实性检测系统、检测方法及存储介质
US20180100947A1 (en) * 2016-10-06 2018-04-12 The Curators Of The University Of Missouri Spectral analysis of surface waves to detect subsurface voids

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2337591B (en) * 1998-05-20 2000-07-12 Geco As Adaptive seismic noise and interference attenuation method
AU2014221254B2 (en) * 2007-06-29 2016-04-14 Geco Technology B.V. Estimating and Using Slowness Vector Attributes in Connection with a Multi-Component Seismic Gather
AU2009282330B2 (en) * 2008-08-11 2013-10-10 Exxonmobil Upstream Research Company Estimation of soil properties using waveforms of seismic surface waves
US20140236487A1 (en) * 2013-02-21 2014-08-21 Westerngeco L.L.C. Methods and computing systems for processing seismic data
US10605937B2 (en) * 2016-05-26 2020-03-31 Cgg Services Sas Device and method for smart picking surface waves dispersion curves
CN108064348B (zh) * 2017-10-12 2020-05-05 南方科技大学 一种基于两点射线追踪的地震走时层析反演方法
CN107561589B (zh) * 2017-10-25 2019-04-30 中国石油化工股份有限公司 一种近地表横波层速度模型建立方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104730579A (zh) * 2013-12-18 2015-06-24 中国石油化工股份有限公司 一种基于表层横波速度反演的纵横波联合静校正方法
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
WO2016187252A1 (en) * 2015-05-20 2016-11-24 Conocophillips Company Surface wave tomography using sparse data acquisition
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 安徽理工大学 一种基于折射波与面波联合反演剩余煤层厚度的方法
CN106950599A (zh) * 2017-05-08 2017-07-14 北京瑞威工程检测有限公司 一种隧道基底密实性检测系统、检测方法及存储介质

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111723491A (zh) * 2020-06-29 2020-09-29 重庆大学 一种基于任意分层土壤格林函数的接地参数获取方法
CN111723491B (zh) * 2020-06-29 2024-03-19 重庆大学 一种基于任意分层土壤格林函数的接地参数获取方法
CN116027411A (zh) * 2021-10-25 2023-04-28 中国石油化工股份有限公司 一种数据驱动的面波预测方法、装置及电子设备和系统
CN114994755A (zh) * 2022-05-25 2022-09-02 平安煤炭开采工程技术研究院有限责任公司 一种无损高效的矸石山地质结构调查方法
CN115236738A (zh) * 2022-06-27 2022-10-25 中国地质科学院地球物理地球化学勘查研究所 一种基于汉克尔变换的勒夫波频散提取方法
CN119846714A (zh) * 2025-01-25 2025-04-18 中国矿业大学 一种边坡裂隙被动源震波成像方法
CN119720604A (zh) * 2025-02-27 2025-03-28 华东交通大学 一种非线性超表面瑞利波频散关系计算方法及系统
CN119720604B (zh) * 2025-02-27 2025-05-13 华东交通大学 一种非线性超表面瑞利波频散关系计算方法及系统
CN120370399A (zh) * 2025-06-26 2025-07-25 湖北省水文地质工程地质勘察院有限公司 基于震源的面波体波综合勘探方法及勘探系统
CN120370399B (zh) * 2025-06-26 2025-08-22 湖北省水文地质工程地质勘察院有限公司 基于震源的面波体波综合勘探方法及勘探系统

Also Published As

Publication number Publication date
US20210141113A1 (en) 2021-05-13
CN111164462B (zh) 2022-05-06
US11137511B2 (en) 2021-10-05
CN111164462A (zh) 2020-05-15

Similar Documents

Publication Publication Date Title
WO2020029015A1 (zh) 一种人工源面波勘探方法、面波勘探装置及终端设备
CN109923440B (zh) 面波勘探方法及终端设备
CN104678435A (zh) 一种提取Rayleigh面波频散曲线的方法
CN104570066B (zh) 地震反演低频模型的构建方法
WO2017024523A1 (zh) 一种射线弹性参数的反演方法
CN111458749A (zh) 应用于被动源地震勘探的面波与体波的分离方法及系统
CN107894613A (zh) 弹性波矢量成像方法、装置、存储介质及设备
CN103926623B (zh) 一种压制逆时偏移低频噪音的方法
WO2021174434A1 (zh) 一种震电波场联合提取瑞雷波频散特征的面波勘探方法
CN114740528B (zh) 一种超微分拉普拉斯块约束的叠前多波联合反演方法
CN113075747A (zh) 储层裂缝发育区域的预测方法及装置
CN112231974B (zh) 基于深度学习的tbm破岩震源地震波场特征恢复方法及系统
CN102590856A (zh) 基于小波频谱分析的位场异常分离方法
Mehra et al. Acoustic pulse propagation in an urban environment using a three-dimensional numerical simulation
CN113238280A (zh) 一种基于格林函数的地震监测方法
CN107861156B (zh) 绕射波的提取方法及装置
CN114578445A (zh) 基于重力资料确定断裂位置的方法及装置
CN111965707A (zh) 一种复杂构造含逆掩断裂的地震反演储层预测方法
CN116520418A (zh) 一种弹性波角度域共成像点道集高效提取方法
CN113433588B (zh) 一种分偏移距扫描叠加的近地表速度分析方法
CN118884531B (zh) 提取高频面波群速度方法、装置、设备、存储介质及产品
CN105589099B (zh) 一种盲源地震波场的多边形带通滤波方法
CN104280774A (zh) 一种单频地震散射噪声的定量分析方法
CN109856672B (zh) 基于深度波数谱的瞬变波包提取方法、存储介质与终端
CN114488346B (zh) 一种地下空间异常体探测方法、装置、设备及介质

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 18929491

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18929491

Country of ref document: EP

Kind code of ref document: A1