WO2017015954A1 - 一种地震信号处理方法、装置和系统 - Google Patents

一种地震信号处理方法、装置和系统 Download PDF

Info

Publication number
WO2017015954A1
WO2017015954A1 PCT/CN2015/085576 CN2015085576W WO2017015954A1 WO 2017015954 A1 WO2017015954 A1 WO 2017015954A1 CN 2015085576 W CN2015085576 W CN 2015085576W WO 2017015954 A1 WO2017015954 A1 WO 2017015954A1
Authority
WO
WIPO (PCT)
Prior art keywords
seismic signal
reflected wave
offset
wave arrival
arrival time
Prior art date
Application number
PCT/CN2015/085576
Other languages
English (en)
French (fr)
Inventor
张金海
Original Assignee
中国科学院地质与地球物理研究所
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 中国科学院地质与地球物理研究所 filed Critical 中国科学院地质与地球物理研究所
Priority to PCT/CN2015/085576 priority Critical patent/WO2017015954A1/zh
Priority to CN201580083576.1A priority patent/CN108260359A/zh
Priority to US15/533,247 priority patent/US20170336523A1/en
Publication of WO2017015954A1 publication Critical patent/WO2017015954A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/40—Transforming data representation
    • G01V2210/41—Arrival times, e.g. of P or S wave or first break
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/52—Move-out correction
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/58—Media-related
    • G01V2210/586—Anisotropic media
    • 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/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time
    • 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/62—Physical property of subsurface
    • G01V2210/626—Physical property of subsurface with anisotropy

Definitions

  • the present invention relates to the field of seismic exploration, and in particular, to a seismic signal processing method, apparatus and system.
  • Seismic exploration is the primary means of finding oil and gas.
  • the seismic signal includes an artificial source seismic signal and a reflected seismic signal from a subsurface target
  • the processed seismic signal refers to a reflected seismic signal from the underground received by the detector.
  • VTI anisotropic medium refers to a transversely isotropic medium with a vertical axis of symmetry, which is the most commonly accepted case of various anisotropic models, and has a good process with the deposition of underground media and the formation of shale.
  • anisotropic parameter refers to the parameter describing the different propagation speeds of seismic waves along different directions of the medium
  • vertical propagation velocity refers to the propagation velocity of seismic waves along the vertical axis of symmetry of the medium.
  • Anisotropic parameters and vertical propagation velocity can provide reliable constraints for the lithology analysis of underground rocks, and thus analyze the structural information and even seismic attribute information of underground targets.
  • the extracted anisotropic parameters also provide the necessary anisotropy initial model for migration imaging and inversion.
  • the attribute parameters of the subsurface medium such as the vertical propagation velocity v 0 , the anisotropic parameters ⁇ and ⁇ , are extracted by the dynamic correction technique.
  • Dynamic correction refers to the process of eliminating the difference between the travel time (or time) of the seismic wave and the time t 0 of the shot position.
  • the time (or time) here refers to the time that the seismic wave passes from the source to the observation point.
  • Existing motion correction methods such as the dynamic correction method based on the Dix formula, Siliqi (2001) directly utilizes the method of non-hyperbolic approximation, and Ursin and Stovas (2006) use the method of continuous fractionation to extract anisotropic parameters. There are defects in the precision and only applicable to weak anisotropic media. Therefore, the prior art seismic signal processing method cannot simultaneously deal with the case of long offset and strong anisotropy, and cannot be processed under high precision. A case with either one of long offset and strong anisotropy.
  • embodiments of the present invention are expected to provide a seismic signal processing method, apparatus and system capable of processing seismic signals with long offset and strong anisotropy, and having a long bias In the case of either shifting and strong anisotropy, the accuracy of the method for processing seismic signals is higher.
  • Embodiments of the present invention provide a seismic signal processing method, including:
  • a non-hyperbolic correction formula based on Pade approximation is constructed according to the offset of the seismic signal reflected by the sampling point and the corresponding reflected wave arrival time;
  • the vertical propagation velocity and the anisotropy parameter of the reflected seismic signal are extracted according to the constructed Pade approximation-based non-hyperbolic correction formula.
  • the Pade approximation based on the correspondence between the normalized offset and the normalized reflected wave arrival time includes:
  • x is the normalized offset
  • ⁇ (x) is the normalized reflected wave arrival
  • n is the order of the Pade approximation
  • P k and Q k are the kth order undetermined coefficients
  • the extracting the vertical propagation velocity and the anisotropy parameter of the reflected seismic signal according to the Pade approximation-based non-hyperbolic correction formula includes:
  • the vertical propagation velocity and the anisotropy parameter corresponding to the seismic signal are obtained according to the actual offset distance and the corresponding reflected wave arrival.
  • the method when the seismic signal includes a plurality of geological layers, the method further includes:
  • the seismic signal is subjected to a layer peeling treatment in an anisotropic velocity analysis of the multilayer.
  • the embodiment of the invention further provides a seismic signal processing device, which comprises: an information acquisition module, a formula construction module, and a parameter extraction module; wherein
  • An information acquiring module configured to acquire an offset of a seismic signal reflected by the sampling point and a corresponding reflected wave arrival time
  • a formula building module is configured to construct a non-hyperbolic correction formula based on the Pade approximation according to the offset of the seismic signal reflected by the sampling point and the corresponding reflected wave arrival time;
  • a parameter extraction module configured to extract a vertical propagation velocity and an anisotropy parameter of the reflected seismic signal according to the constructed Pade approximation-based non-hyperbolic correction formula.
  • the formula building module includes:
  • a normalization unit configured to normalize the offset of the acquired seismic signal and the reflected wave arrival time
  • Pade-type building unit for constructing Pade approximations of their correspondence based on normalized offset and normalized reflected wave arrival.
  • the parameter extraction module includes:
  • a formula scanning unit configured to obtain a normalized reflected wave arrival time ⁇ (x) by the constructed Pade approximation-based non-hyperbolic correction formula
  • the result obtaining unit is configured to obtain a vertical propagation speed and an anisotropy parameter corresponding to the seismic signal according to the converted actual offset distance and the corresponding reflected wave arrival time.
  • the seismic signal processing device further includes:
  • a layer stripping module for layer stripping of seismic signals in an anisotropic velocity analysis of multiple layers.
  • the present invention also provides a seismic signal processing system including any one of the above seismic signal processing devices.
  • the present invention processes the seismic signal by using the non-hyperbolic correction formula based on the Pade approximation, and can process the underground medium with strong anisotropy and long offset. Situation; and in the case of one of strong anisotropy and long offset, the accuracy of the processing method of the existing seismic signal is higher.
  • FIG. 1 is a schematic diagram of a model of a VTI medium according to an embodiment of the present invention.
  • FIG. 2 is a schematic flowchart of a seismic signal processing method according to an embodiment of the present invention.
  • FIG. 3 is a schematic diagram showing a comparison between a reflected wave arrival time curve and a ray tracing result obtained by the Pade approximation-based motion correction method and the existing motion correction method according to an embodiment of the present invention
  • Figure 4 is a schematic diagram showing the distribution of Thomsen anisotropic parameters measured in nature
  • FIG. 5 is a seismic record overlay diagram of ray tracing and finite difference simulation according to an embodiment of the present invention
  • FIG. 6 is a schematic diagram of propagation time errors of a Pade approximation based motion correction method and an existing motion correction method according to an embodiment of the present invention
  • FIG. 7 is a schematic diagram of a result of performing an Alkhalifah method scan on a seismic data according to an embodiment of the present invention.
  • FIG. 8 is a schematic diagram of a result of scanning a seismic data using a Pade[7,7] approximation method according to an embodiment of the present invention
  • FIG. 9 is a schematic structural diagram of a seismic signal processing apparatus according to an embodiment of the present invention.
  • the model of the VTI medium involved is as shown in FIG. 1.
  • the distance X of the shot point from the detector is the offset distance
  • the model thickness is D
  • the propagation speed of the seismic wave in the model is v 0 .
  • ⁇ and ⁇ are Thomsen anisotropic parameters.
  • the processed seismic signal refers to an electrical signal obtained by the detector receiving the reflected seismic wave from the underground target through the detector.
  • the other three existing methods for extracting anisotropic parameters are respectively based on the dynamic correction method of the Dix formula, and Siliqi (2001) directly utilizes the method of non-hyperbolic approximation, Ursin and Stovas. (2006) The method of using the continuous expansion method.
  • the dynamic correction method based on the Dix formula uses the following non-hyperbolic approximation formula to calculate the square of the reflected wave arrival time of different offsets:
  • X is the actual offset (ie the distance from the shot to the detection point)
  • ⁇ nmo is moveout velocity (normal moveout (NMO) velocity)
  • t 0 is a reflected wave to zero-offset case
  • [eta] non ellipsometric parameters It is defined as follows (Alkhalifah and Tsvankin, 1995):
  • ⁇ and ⁇ are Thomsen anisotropic parameters. When they are not zero, they will cause anisotropy in the propagation of seismic waves, and the larger their absolute values, the more obvious the anisotropic effect of geological medium; Siliqi( 2001) Directly using the non-hyperbolic approximation to obtain the arrival time of the reflected wave, the form is as follows:
  • Ursin and Stovas (2006) use the continuous expansion to obtain the square of the reflected wave arrival time, in the form of:
  • FIG. 2 is a schematic flow chart of a seismic signal processing method provided by the present invention. As shown in FIG. 2, the method includes:
  • Step 201 Obtain an offset distance of the seismic signal reflected by the sampling point and a corresponding reflected wave arrival time
  • the reflected seismic signal received by the detector is used to obtain the offset of the reflected seismic signal at different sampling points and the corresponding reflected wave arrival time;
  • Step 202 Construct a non-hyperbolic correction formula based on Pade approximation according to the offset distance of the seismic signal and the arrival time of the reflected wave according to the acquired sampling point;
  • step 202 includes the following steps:
  • Step A normalizing the offset of the acquired seismic signal and the reflected wave arrival time
  • Step B constructing a Pade approximation of the corresponding relationship based on the normalized offset and the normalized reflected wave arrival time;
  • Step C determining the undetermined coefficients P k and Q k ;
  • Step 203 Extract a vertical propagation velocity and an anisotropy parameter of the reflected seismic signal according to the constructed Pade approximation-based non-hyperbolic correction formula
  • step 203 includes the following steps:
  • Step a obtaining a normalized reflected wave arrival time ⁇ (x) (time-walking parameter) by the constructed non-hyperbolic correction formula based on the Pade approximation;
  • the non-hyperbolic correction formula can calculate the normalized reflected wave corresponding to the zero offset to the time ⁇ ;
  • Figure 3 shows the method based on Pade approximation in a specific model.
  • Siliqi (2001) directly uses non-hyperbolic approximation.
  • Figure 4 the ray tracing (dashed line) and the finite difference simulation of the seismic record overlay are shown.
  • the model used is a VTI medium model with a model thickness of 500 m.
  • the shot point and the detection point are both on the surface.
  • Figure 5 is a measured Thomsen anisotropic parameter distribution map.
  • the anisotropic parameters corresponding to the two "pentagrams" in the figure are the two sets of anisotropic parameters used in the comparison with the existing methods in the embodiment of the present invention.
  • Fig. 6 gives the reflected wave arrival time error curves of various approximation methods compared with ray tracing on the basis of Fig. 3, and more clearly demonstrates the wider anisotropy parameter using the Pade approximation method. The advantage of higher accuracy over previous methods over distributed and wider offset ranges.
  • step c the vertical propagation velocity and the anisotropy parameter corresponding to the seismic signal are obtained according to the converted actual offset distance and the corresponding reflected wave arrival time.
  • the method of determining the anisotropic parameter according to the reflected wave arrival time ⁇ belongs to the conventional method in the art, and therefore will not be redundantly described herein.
  • Fig. 7 is a view showing the result of performing the existing Alkhalifah method scan on the seismic data derived from the model, wherein the left column is the seismic data, the middle column is the scan result of the equivalent anisotropy parameter ⁇ , and the right column is The result of the velocity scan;
  • FIG. 8 is a schematic diagram showing the result of scanning the seismic data derived from the model by Pade[7,7], the left column is the seismic data, and the middle column is the equivalent anisotropic parameter ⁇ . The result of the scan, the right column is the speed scan result.
  • the seismic data derived from the model here is simulated by the finite difference method according to the elastic wave equation.
  • Fig. 7 the result obtained by the existing Alkhalifah method is 0.39 with a deviation of 22%; as shown in Fig. 8, the result of approximating Pade [7, 7] is 0.5, and the deviation is 0%. It can be seen that the ⁇ precision obtained by Pade[7,7] approximation is higher than that obtained by the Alkhalifah method.
  • the method further includes:
  • the seismic signal is subjected to a layer peeling treatment in an anisotropic velocity analysis of the multilayer.
  • multi-layer anisotropy velocity analysis can be performed, and the layer of the current layer is subjected to velocity independent layer peeling, and the multilayer layer is The reflection record is converted into a single-layer reflection record, and then the velocity analysis process of the single-layer model is repeated until all the reflections are processed, and a multilayer anisotropic velocity analysis is realized.
  • FIG. 9 is a schematic structural diagram of a seismic signal processing apparatus according to an embodiment of the present invention. As shown in FIG. 9, the seismic signal processing apparatus includes:
  • the information acquiring module 901 is configured to acquire an offset of the reflected seismic signal of the sampling point and a corresponding reflected wave arrival time;
  • the formula construction module 902 is configured to reflect the offset of the seismic signal and the corresponding reflection according to the sampling point When the arrival is reached, a non-hyperbolic correction formula based on Pade approximation is constructed;
  • the parameter extraction module 903 is located between the information acquisition module and the formula construction module, and is configured to extract the vertical propagation velocity and the anisotropy parameter of the reflected seismic signal according to the constructed Pade approximation-based non-hyperbolic correction formula.
  • the formula construction module 902 includes:
  • a normalization unit configured to normalize the offset of the acquired seismic signal and the reflected wave arrival time
  • Pade-type building unit for constructing Pade approximations of their correspondence based on normalized offset and normalized reflected wave arrival.
  • the parameter extraction module 903 includes:
  • a formula scanning unit configured to obtain a normalized reflected wave arrival time ⁇ (x) by the constructed Pade approximation-based non-hyperbolic correction formula
  • the result obtaining unit is configured to obtain a vertical propagation speed and an anisotropy parameter corresponding to the seismic signal according to the converted actual offset distance and the corresponding reflected wave arrival time.
  • the processing device for the seismic signal further includes:
  • a layer stripping module for layer stripping of seismic signals in an anisotropic velocity analysis of multiple layers.
  • the velocity independent layer peeling is performed on the recording of the current layer, the multi-layer reflection record is converted into a single layer reflection record, and then the single layer model analysis process is repeated until All of the reflections are processed, thereby enabling anisotropic velocity analysis of the multilayer.
  • the embodiment of the invention further provides a seismic signal processing system, which includes any one of the above seismic signal processing devices.
  • the present invention processes the seismic signal by using the non-hyperbolic correction formula based on the Pade approximation, and can process the underground medium with strong anisotropy and long offset. Situation; and in the case of one of strong anisotropy and long offset, the accuracy of the processing method of the existing seismic signal is higher.
  • embodiments of the present invention can be provided as a method, system, or computer program product. Accordingly, the present invention can take the form of a hardware embodiment, a software embodiment, or a combination of software and hardware. Moreover, the invention can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage and optical storage, etc.) including computer usable program code.
  • the computer program instructions can also be stored in a computer readable memory that can direct a computer or other programmable data processing device to operate in a particular manner, such that the instructions stored in the computer readable memory produce an article of manufacture comprising the instruction device.
  • the apparatus implements the functions specified in one or more blocks of a flow or a flow and/or block diagram of the flowchart.
  • These computer program instructions can also be loaded onto a computer or other programmable data processing device such that a series of operational steps are performed on a computer or other programmable device to produce computer-implemented processing for execution on a computer or other programmable device.
  • the instructions provide steps for implementing the functions specified in one or more of the flow or in a block or blocks of a flow diagram.

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)
  • Materials For Photolithography (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种地震信号的处理方法、装置和系统，其中，所述方法包括：获取采样点反射地震信号的偏移距和对应的反射波到时；根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式；根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。

Description

一种地震信号处理方法、装置和系统 技术领域
本发明涉及地震勘探领域，尤其涉及一种地震信号处理方法、装置和系统。
背景技术
地震勘探是寻找石油和天然气的最主要手段。人们可以在地表钻孔中埋置炸药，通过引爆炸药激发一个近于脉冲的地震信号(即激发人工震源)，然后利用埋置在地表的检波器接收来自地下目标的反射地震信号。最后，通过对所接收的地震信号的处理，获得地下目标的结构信息甚至地震属性信息，从而为后续的钻井布置提供重要参考。这里，所述地震信号包括人工震源地震信号和来自地下目标的反射地震信号，所处理的地震信号是指检波器接收的来自地下的反射地震信号。
当前，基于地下介质为VTI各向异性介质的假设，从地震信号提取地质层的各向异性参数和垂向传播速度等，是地震信号处理中的重要环节。所谓VTI各向异性介质是指具有垂直对称轴的横向各向同性介质，它是各种各向异性模型中最为普遍接受的情形，与地下介质的沉积、页岩的形成等过程具有很好的对应性；而各向异性参数是指描述地震波沿介质不同方向的传播速度不同的参数；垂向传播速度是指地震波沿介质垂直对称轴的传播速度。各向异性参数和垂向传播速度能够为地下岩石的岩性分析给出可靠的约束，从而据此分析地下目标的结构信息甚至地震属性信息。同时，所提取的各向异性参数也为偏移成像和反演提供了必要的各向异性初始模型。
地下介质的属性参数如垂向传播速度v0、各向异性参数ε和δ，是通过动校正技术提取的。动校正是指消除地震波的走时(或到时)较炮点位置的到时t0之间差异的过程。这里走时(或到时)是指地震波从震源传到观测点所经过的时间。现有的动校正方法，例如基于Dix公式的动校正方法，Siliqi(2001)直接利用非双曲近似的方法，以及Ursin和Stovas(2006)利用连分式展开的方法，在提取各向异性参数方面存在精度不高和只适用于弱各向异性介质的缺陷，因此现有技术的地震信号处理方法不仅无法同时处理具有长偏移距和强各向异性的情形，也无法在高精度下处理具有长偏移距和强各向异性二者之一的情形。
发明内容
为解决现有存在的技术问题，本发明实施例期望提供一种地震信号处理方法、装置和系统，能在具有长偏移距和强各向异性的情形下处理地震信号，并且在具有长偏移距和强各向异性二者之一的情形下比现有处理地震信号的方法精度更高。
本发明实施例的技术方案是这样实现的：
本发明实施例提供了一种地震信号处理方法，该方法包括：
获取采样点反射地震信号的偏移距和对应的反射波到时；
根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式；
根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。
上述方案中，所述根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式包括：
对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式。
上述方案中，所述基于归一化偏移距和归一化反射波到时的对应关系的Pade近似式包括：
Figure PCTCN2015085576-appb-000001
其中，x为归一化偏移距，τ(x)为归一化反射波到时，n为Pade近似式的阶数，Pk与Qk为第k阶待定系数；
上述方案中，所述根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数包括：
由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)；
利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时；
根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
上述方案中，当地震信号中包含多个地质层时，所述方法还包括：
在多层的各向异性速度分析中对地震信号进行层剥离处理。
本发明实施例还提供一种地震信号处理装置，该装置包括：信息获取模块、公式构建模块、参数提取模块；其中，
信息获取模块，用于获取采样点反射地震信号的偏移距和对应的反射波到时；
公式构建模块，用于根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式；
参数提取模块，用于根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。
上述方案中，所述公式构建模块包括：
归一化单元，用于对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
Pade式构建单元，用于基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式。
上述方案中，所述参数提取模块包括：
公式扫描单元，用于由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)；
参数换算单元，用于利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时；
结果求取单元，用于根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
上述方案中，所述地震信号处理装置还包括：
层剥离模块，用于在多层的各向异性速度分析中对地震信号进行层剥离处理。
本发明还提供一种地震信号的处理系统，该系统中包括上述任意一种地震信号处理装置。
本发明实施例提供地震信号处理方法、装置和系统的优势是：本发明利用基于Pade近似的非双曲动校正公式处理地震信号，可以处理地下介质同时具有强各向异性和长偏移距的情形；并且在有强各向异性和长偏移距二者之一的情形时，比现有地震信号的处理方法精度更高。
附图说明
图1为本发明实施例提供的所涉及的VTI介质的模型示意图；
图2为本发明实施例提供的地震信号处理方法的流程示意图；
图3为本发明实施例提供的基于Pade近似的动校正方法与现有的动校正方法得到的反射波到时曲线与射线追踪结果的对比示意图；
图4为自然界实测的Thomsen各向异性参数分布示意图；
图5为本发明实施例提供的射线追踪与有限差分模拟的地震记录叠合图；
图6为本发明实施例提供的基于Pade近似的动校正方法与现有的动校正方法的传播时间误差示意图；
图7为本发明实施例提供的对一地震数据进行Alkhalifah方法扫描的结果示意图；
图8为本发明实施例提供的对一地震数据进行利用Pade[7,7]近似的方法扫描的结果示意图；
图9为本发明实施例提供的地震信号处理装置的结构示意图。
具体实施方式
为了更清楚地说明本发明实施例和技术方案，下面将结合附图及实施例对 本发明的技术方案进行更详细的说明，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明的实施例，本领域普通技术人员在不付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明实施例中，所涉及的VTI介质的模型如图1所示，图中炮点距检波器的距离X为偏移距，模型厚度为D，地震波在模型中的传播速度为v0，ε和δ为Thomsen各向异性参数。
在本发明实施例中，所处理的地震信号，是指检波器接收的来自地下目标的反射地震波经检波器转化而成的电信号。
在本发明实施例中，所比较的另外三种现有的提取各向异性参数的方法，分别是基于Dix公式的动校正方法，Siliqi(2001)直接利用非双曲近似的方法，Ursin和Stovas(2006)利用连分式展开的方法。其中基于Dix公式的动校正方法采用如下的非双曲近似公式来计算不同偏移距的反射波到时的平方：
Figure PCTCN2015085576-appb-000002
这里τ(x)＝t(X)/t0是归一化的反射波到时，X是实际的偏移距(即炮点到检波点的距离)，x＝X/(t0υnmo)是归一化的偏移距(Stovas，2006)，υnmo是动校正速度(normal moveout(NMO)velocity)，t0是零偏移距情形的反射波到时，η为非椭圆参数，其定义如下(Alkhalifah and Tsvankin，1995)：
Figure PCTCN2015085576-appb-000003
其中ε和δ是Thomsen各向异性参数，它们不为零时就会导致地震波的传播出现各向异性现象，而且它们的绝对值越大，表示地质介质的各向异性效应就越明显；Siliqi(2001)直接利用非双曲近似获得反射波的到时，形式如下：
Figure PCTCN2015085576-appb-000004
Ursin和Stovas(2006)利用连分式展开获得反射波到时的平方，形式如下：
Figure PCTCN2015085576-appb-000005
如上几类方法是业界使用最为广泛的动校正的代表性方法，但是，这些方法有如下缺陷：不适用于当地质介质具有强各向异性(|ε|≥0.2和/或|δ|≥0.2)和并且偏移距较大(偏移距/目标层深度>2)的情形，具体而言，只适用于各向异性参数|ε|＜0.2和|δ|＜0.2，和/或偏移距/目标层深度<2时的情形，超出上述范围则导致动校正公式误差过大而无法使用。
图2为本发明提供的地震信号处理方法的流程示意图，如图2所示，该方法包括：
步骤201，获取采样点反射地震信号的偏移距和对应的反射波到时；
具体的，通过检波器接收的反射地震信号，获取不同采样点反射地震信号的偏移距和对应的反射波到时；
步骤202，根据获取的采样点反射地震信号的偏移距和反射波到时，构建基于Pade近似的非双曲动校正公式；
具体的，步骤202包括如下步骤：
步骤A，对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
具体的，定义归一化偏移距为x，由x＝X/t0υnmo给出；定义归一化反射波到时为τ(x)，由τ(x)＝t(X)/t0给出；
步骤B，基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式；
具体的，Pade近似式的一般表达形式为：
Figure PCTCN2015085576-appb-000006
则将归一化偏移距x作为上式中的特征值，归一化反射波到时τ(x)作为x的函数，利用对角Padé近似构建非双曲动校正公式：
Figure PCTCN2015085576-appb-000007
这里，K＝L＝n且n为任一大于0的整数，Pk与Qk为待定系数。
步骤C，确定所述的待定系数Pk与Qk；
具体的，确定待定系数Pk与Qk的过程如下：
对已知的精确时距函数f(x)做泰勒展开：
Figure PCTCN2015085576-appb-000008
当满足f(x)-RLM(x)＝O(xL+M+1)，即有理展开RLM(x)同原函数f(x)的误差是xL+M+1的高阶小量(x＜1)时，RLM(x)即为f(x)的Pade近似。将f(x)-RLM(x)＝O(xL+M+1)两端与RLM(x)的分母相乘，并比较方程两端同阶次xk的系数，就可以得到系数Pk(k＝0,2,…,L)和Qk(k＝0,2,…,M)。
特别的，在一个实施例中，以n＝4的Pade[4,4]近似式来确定归一化偏移距x和归一化反射波到时τ(x)的对应关系，即
Figure PCTCN2015085576-appb-000009
利用上述方法求取的系数Pk和Qk(k＝1,2,3,4)分别为：
Figure PCTCN2015085576-appb-000010
Figure PCTCN2015085576-appb-000011
Figure PCTCN2015085576-appb-000012
Figure PCTCN2015085576-appb-000013
Figure PCTCN2015085576-appb-000014
Figure PCTCN2015085576-appb-000015
Figure PCTCN2015085576-appb-000016
Figure PCTCN2015085576-appb-000017
在另一个实施例中，以n＝7的Pade[7,7]近似式来确定，归一化偏移距x和归一化反射波到时τ(x)的对应关系，即
Figure PCTCN2015085576-appb-000018
则，求取的系数Pk和Qk(k＝1,2,3,4,5,6,7)分别为：
P1＝(1928934+177875271×η+7837827462
×η2+219632829372×η3+4393725401232×η4+66762558665424×η5+800282483624416
×η6+7757101152212224
×η7+61816601652878720×η8+409500878760915968×η9+2270674620269779456
×η10+10577274688104635392×η11+41421467806172823552
×η12+136091600042563719168×η13+373312046919000924160
×η14+848038798280195538944×η15+1576067219935627804672
×η16+2354425649632364789760×η17+2755347906978845556736
×η18+2429849593769511354368×η19+1516096215071553748992
×η20+595558663830221357056×η21+110526504681837953024×η22)/(275562+23048793
×η+922612194×η2+23507343720×η3+427711600368×η4+5909128384752
×η5+64339900368352×η6+565512119070400×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560
×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
P2＝(5786802+575589969×η+27402488892×η2+830945491596×η3+18017393437944×η4+297266
111100144×η5+3876968360705088×η6+40985450564014144
×η7+357252032884224384×η8+2597717830723195392×η9+15879010903089374208
×η10+81969396932994297856×η11+358017626105963333632
×η12+1322357168238845911040×η13+4118100751204746092544
×η14+10752767982215247659008×η15+23337340412356598693888
×η16+41574115402550069166080×η7+59708407067590095798272
×η18+67362589572056487559168×η19+57416164905232919166976
×η20+34711631653490820907008×η21+13244594083153838080000
×η22+2393573928540571697152×η23)/(275562+23048793
×η+922612194×η2+23507343720×η3+427711600368×η4+5909128384752
×η5+64339900368352×η6+565512119070400×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560
×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
P3＝(9644670+1016528535×η+51423263424×η2+1661413259772×η3+38484900883368
×η4+680152787928000×η5+9528602090114432×η6+108526847831855296
×η7+1022520875817401856×η8+8066274270183933696×η9+53717060689940852224
×η10+303577662557579200512×η11+1459995142346725052416×η12+5978678247368511492096
×η13+20814233677401535979520×η14+61375505707148677595136
×η15+152343562780817632362496×η16+315415772152951498211328
×η17+537719364162729998548992×η18+741188745099913124642816
×η19+804734216362744583553024×η20+661814793824070652657664
×η21+387035816399908829659136×η22+143219393455221827960832
×η23+25168401484774230720512×η24)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752
×η5+64339900368352×η6+565512119070400×η7+4076204157164160×η8+24339343255939840
×η9+121096428605260288×η10+503183373567916032×η11+1744689526345007104
×η12+5027600355175223296×η13+11950937819425972224×η14+23161421974674882560
×η15+35970414631180271616×η16+43640591992815747072×η17+39793720202662510592
×η18+25609553049671958528×η19+10350908157396516864×η20+1971601553789812736×η21)
P4＝(9644670+1061012115×η+56187769770×η2+1906082406360
×η3+46500030064440×η4+868173241337712×η5+12889555987110176
×η6+156091894087908096×η7+1569122021197143680×η8+13256308026663413760
×η9+94931570607418164224×η10+579580149141234107392×η11+3027009907748502867968
×η12+13543073431265094049792×η13+51881658354936538816512
×η14+169777709930246048923648×η15+472523780210153518759936
×η16+1111087418606925566705664×η17+2186403098161881572966400
×η18+3553307500792885418131456×η19+4682291860203479095050240
×η20+4873064867923789925580800×η21+3851662186402023565426688
×η22+2170557943353311710150656×η23+776096351861666830352384
×η24+132168601394774970204160×η25)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104
×η12+5027600355175223296×η13+11950937819425972224×η14+23161421974674882560
×η15+35970414631180271616×η16+43640591992815747072×η17+39793720202662510592
×η18+25609553049671958528×η19+10350908157396516864×η20+1971601553789812736×η21)
P5＝(5786802+655660413×η+35856061830×η2+1259630524044
×η3+31916132469504×η4+620788649928336×η5+9632260563079296
×η6+122305923730049536×η7+1293561923916277632×η8+11539439285652457728
×η9+87595635509731739136×η10+569272882492753199104×η11+3179636656647032514560
×η12+15293848493477868236800×η13+63367358384507785412608
×η14+225861235067735815340032×η15+690448514286351392735232
×η16+1801417024902636767477760
×η17+3983110792356283753693184×η18+7390979589931870290182144
×η19+11355991893219311783247872×η20+14182024077702696639922176
×η21+14022149691826641035067392×η22+10554855242750478780465152
×η23+5679097991222977795981312×η24+1944192338541529440190464
×η25+318006564232006006210560×η26)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560×η15+35970414631180271616
×η16+43640591992815747072×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
P6＝(1928934+222358851×η+12398103000×η2+445102784100
×η3+11554334143416×η4+230871063265776×η5+3690547741618528×η6+48423553021297536
×η7+530898122265830400×η8+4925505742734286592
×η9+39019858824218919424×η10+265608094280186059776×η11+1559949387124244523008
×η12+7923397797552678592512×η13+34832234201656515723264
×η14+132440615814478859452416×η15+434609242879999474434048
×η16+1226339352391968944357376×η17+2959331825491941912346624
×η18+6061451039088885284732928×η19+10431257171576328759017472
×η20+14876750084992650232987648×η21+17255659150501004091326464
×η22+15852942901944734233657344×η23+11092112127981678147665920
×η24+5550124461870301593993216×η25+1768179670704684621889536
×η26+269466871890791760920576×η27)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560×η15+35970414631180271616
×η16+43640591992815747072×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
P7＝(275562+31945509×η+1793821140×η2+64961412588×η3+1704160627776
×η4+34483134232080×η5+559500187781824×η6+7470118654010240×η7+83561375615403648
×η8+793207525502117632×η9+6447977398614715904×η10+45172677920381583360×η
11+273887607136658833408×η12+1440713955626014019584×η13+6581074491526837641216
×η14+26093407726900327022592
×η15+89638299455484154675200×η16+265943700421102279262208
×η17+678189291691971688529920×η18+1476802093336346384400384
×η19+2721935498188550444679168×η20+4196921622869561497878528
×η21+5329335799228928878444544×η22+5455238404086497720401920
×η23+4367657736695455202410496×η24+2615332663994183525072896
×η25+1089922493261880949211136×η26+277010504006563182673920
×η27+31241190228217552699392×η28)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560×η15+35970414631180271616
×η16+43640591992815747072×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
Q1＝(2(826686+77413239×η+3457607634×η2+98062742826×η3+1983006900432
×η4+30426715140336×η5+367971291628032×η6+3595794516570912×η7+28870198747857280
×η8+192580767752488064×η9+1074789095832259584×η10+5037045657268359680×η
11+19838389139913908224×η12+65531999843694247936
×η13+180680554549787475968×η14+412438688152760328192×η15+770048402652223766528
×η16+1155392528819774521344×η17+1357777093388091523072
×η18+1202120020359919697920×η19+752872653457078616064×η20+296793531138215772160
×η21+55263252340918976512×η22))/(275562+23048793
×η+922612194×η2+23507343720×η3+427711600368×η4+5909128384752
×η5+64339900368352×η6+565512119070400×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296×η13+11950937819425972224
×η14+23161421974674882560×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528×η19+10350908157396516864
×η20+1971601553789812736×η21)
Q2＝(4133430+421314615×η+20533371210×η2+636665230332×η3+14098394324520
×η4+237268104020208×η5+3152844034218528×η6+33922541331609024
×η7+300642659626650624×η8+2220708703532547584×η9+13778111397936734720
×η10+72137498475668099072×η11+319347214573271349248×η12+1194782547604147429376
×η13+3766794842815521587200×η14+9951792481548578947072
×η15+21843566451001500925952×η16+39335271174172880666624
×η17+57080134064799544246272×η18+65037936971741973184512
×η19+55961638704418105851904×η20+34138746407529182396416
×η21+13138010781579579752448×η22+2393573928540571697152×η23)/(275562+23048793×η
+922612194×η2+23507343720×η3+427711600368×η4+5909128384752×η5+64339900368352
×η6+565512119070400×η7+4076204157164160×η8+24339343255939840
×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296×η13+11950937819425972224
×η14+23161421974674882560×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528×η19+10350908157396516864
×η20+1971601553789812736×η21)
Q3＝(4(1377810+149492385×η+7787535210×η2+259114162518×η3+6180167874096
×η4+112419800046036×η5+1620128530120832×η6+18969150581351168
×η7+183589572803618464×η8+1486526951974910528×η9+10152919806654055040
×η10+58801263872239585280×η11+289584146628032046592×η12+1213430419197128672256
×η13+4319643839733575213056×η14+13015550590974573481984
×η15+32990004246161177829376×η16+69702704174107263369216
×η17+121185468513412439408640×η18+170244660489351100891136
×η19+188263078497808715612160×η20+157589880394364732047360
×η21+93739206410471443791872×η22+35255803329349863604224
×η23+6292100371193557680128×η24))/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560×η15+35970414631180271616
×η16+43640591992815747072×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
Q4＝(4133430+468553815×η+25605680310×η2+897615594648×η3+22656390826320
×η4+438108906241008×η5+6742228840223712×η6+84681107416547328
×η7+883198773982473728×η8+7742645820390831616×η9+57536031426832966656
×η10+364455198712506151936×η11+1974417941419785115648
×η12+9159988806870836191232×η13+36372615088882275975168
×η14+123319839835275326144512×η15+355432883417909228699648
×η16+865041016551896630034432×η17+1760876157015194696155136
×η18+2958545607405671325630464×η19+4027783529552811685576704
×η20+4327745007067682018689024×η21+3528725301548983557029888
×η22+2049606239540481072562176×η23+754593826739385828114432
×η24+132168601394774970204160×η25)/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160×η8+24339343255939840×η9+121096428605260288
×η10+503183373567916032×η11+1744689526345007104×η12+5027600355175223296
×η13+11950937819425972224×η14+23161421974674882560×η15+35970414631180271616
×η16+43640591992815747072×η17+39793720202662510592×η18+25609553049671958528
×η19+10350908157396516864×η20+1971601553789812736×η21)
Q5＝(2(826686+96308919×η+5427009882×η2+196894550034×η3+5164414595280
×η4+104238873016848×η5+1682460628070880×η6+22275247273303488
×η7+246203045442158592×η8+2299954242925295488×η9+18316966430868216832
×η10+125095616610351866368×η11+735309014792072728576×η12+3726590593337458358272
×η13+16285963937114949967872×η14+61280124231354545930240×η
15+197902014621189397643264×η16+545788559575343499411456×η
17+1276194617518571555127296×η18+2505058528646059114889216×η
19+4072316943659743137759232×η20+5381115380052287009849344×η
21+5628762556804139084414976×η22+4481042550617360469852160×η
23+2548584529843927975460864×η24+921493653485101950959616×η
25+159003282116003003105280×η26))/(275562+23048793×η+922612194×η2+23507343720
×η3+427711600368×η4+5909128384752×η5+64339900368352×η6+565512119070400
×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296×η13+11950937819425972224
×η14+23161421974674882560×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528×η19+10350908157396516864
×η20+1971601553789812736×η21)
Q6＝(275562+32496633×η+1857712158×η2+68549054868×η3+1834083452952
×η4+37891455487632×η5+628466456245984×η6+8589030869889280×η7+98494511987346304
×η8+960056884707777280×η9+8027700648044889088×η10+57951806860999945216×η
11+362686086208153231360×η12+1972364022839890993152×η13+9326922280871311966208
×η14+38321172966570882416640
×η15+136490783498625906245632×η16+419830781828325790646272
×η17+1109016177366253329317888×η18+2496753260464358062555136
×η19+4741556136527862118744064×η20+7491432251681512087355392
×η21+9662752109219610376536064×η22+9907682471796831223283712
×η23+7763656208371979326062592×η24+4364326686884151485792256
×η25+1566341751558105880068096×η26+269466871890791760920576
×η27)/(275562+23048793×η+922612194×η2+23507343720×η3+427711600368
×η4+5909128384752×η5+64339900368352×η6+565512119070400×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104×η12+5027600355175223296×η13+11950937819425972224
×η14+23161421974674882560×η15+35970414631180271616×η16+43640591992815747072
×η17+39793720202662510592×η18+25609553049671958528×η19+10350908157396516864
×η20+1971601553789812736×η21)
Q7＝(128×η7(86093442+6841239993×η+259789618440×η2+6270533376300
×η3+107931841249680×η4+1408828996702656×η5+14476512654173952
×η6+119978955650199744×η7+815076528301108224×η8+4587220895934600704
×η9+21525572079369882624×η10+84473360948626243584×η11+277205625021341642752
×η12+758245030780603625472×η13+1717357845011212746752
×η14+3186134285268034453504×η15+4763395695841300316160
×η16+5600076296653792083968×η17+4986276874779718320128
×η18+3161381112437777891328×η19+1271871501522830360576
×η20+244071798657949630464×η21))/(275562+23048793
×η+922612194×η2+23507343720×η3+427711600368×η4+5909128384752
×η5+64339900368352×η6+565512119070400×η7+4076204157164160
×η8+24339343255939840×η9+121096428605260288×η10+503183373567916032
×η11+1744689526345007104
×η12+5027600355175223296×η13+11950937819425972224×η14+23161421974674882560
×η15+35970414631180271616×η16+43640591992815747072×η17+39793720202662510592
×η18+25609553049671958528×η19+10350908157396516864×η20+1971601553789812736
×η21)。
步骤203，根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数；
具体的，步骤203包括如下步骤：
步骤a，由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)(走时参数)；
具体的，根据不同采样点获取的反射波到时，利用所构建的基于Pade近似 的非双曲动校正公式，即可计算零偏移距所对应的归一化反射波到时τ；
为了证明利用Pade近似的精度比前人的方法精度更高，图3给出了在一个具体模型中基于Pade近似的方法，基于Dix公式的动校正方法，Siliqi(2001)直接利用非双曲近似的方法，以及Ursin和Stovas(2006)利用连分式展开的方法得到的走时，其中“Analytical”是利用射线追踪的结果，它用来作为参考，考察上述近似的精度。在图4中，展示了所述的射线追踪(虚线)与有限差分模拟的地震记录叠合图，从图中可以看到，射线追踪的结果同有限差分模拟得到的结果吻合很好，从而可以作为上述近似的精度的参考标准。从图3可见，Pade[4,4]和Pade[7,7]方法都具有很高的精度，尤其是Pade[7,7]的精度在列出的两种情形下都较前人的方法具有明显的提升。
在图3和图4中，所采用的模型为VTI介质模型，模型厚度为500m，炮点和检波点均位于地表，模型速度为2000m/s，对应的零偏移距到时为t0＝0.5s。图5为实测到的Thomsen各向异性参数分布图，图中两个“五角星”处所对应的各向异性参数即为本发明实施例中与已有方法比较时采用的两组各向异性参数，具体为：(a)ε＝0.3,δ＝-0.1,υnmo＝1.7889km/s,η＝0.5和(b)ε＝0.4,δ＝-0.3,υnmo＝1.2649km/s,η＝1.75两组参数，两组数据中各向异性参数均为强各向异性参数。
进一步的，图6在图3的基础上给出了各种近似方法同射线追踪相比较的反射波到时误差曲线，更为清楚的展示了利用Pade近似的方法在更宽广的各向异性参数分布和更广泛的偏移距范围内较以往方法具有更高精度这一优点。
步骤b，利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时。
步骤c，根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
具体的，根据反射波到时τ确定各向异性参数的方法，属于本领域的常规方法故在此不做多余的阐述。
为了证明采用Pade近似获取的各向异性参数精度更高，对相同的由模型导 出的地震数据分别用已有的Alkhalifah方法和利用Pade[7,7]近似的方法进行扫描，分别获取各向异性参数及速度，如图7和图8所示。图7表示对所述的由模型导出的地震数据进行已有的Alkhalifah方法扫描的结果示意图，图中左边一列为地震数据，中间一列为等效各向异性参数η的扫描结果，右侧一列为速度扫描结果；图8表示对所述的由模型导出的地震数据进行利用Pade[7,7]近似的方法扫描的结果示意图，左边一列为地震数据，中间一列为等效各向异性参数η的扫描结果，右侧一列为速度扫描结果。
这里由模型导出的地震数据是根据弹性波方程利用有限差分法模拟得到的。所采用的模型参数如下：v0＝2000m/s，ε＝0.3，δ＝-0.1，η＝(ε-δ)/(1+2δ)＝0.5。如图7所示，已有的Alkhalifah方法得到的结果η为0.39，偏差为22％；如图8所示，利用Pade[7,7]近似的结果η为0.5，偏差为0％。由此可见，利用Pade[7,7]近似得到的η精度比利用Alkhalifah方法得到的精度更高。
进一步的，求取各向异性参数后，就可以获得地下目标的结构信息甚至地震属性信息。
进一步的，当地震信号中包含多个地质层的地震信号时，所述方法还包括：
在多层的各向异性速度分析中对地震信号进行层剥离处理。
具体的，为实现多层各向异性速度分析，通过对地震信号进行层剥离处理，就可以进行多层各向异性速度分析，通过对当前层的记录进行速度独立的层剥离，将多层的反射记录转换为单层的反射记录，然后重复单层模型的速度分析过程，直到将所有的反射全部处理完，就实现了多层的各向异性速度分析。
图9是本发明实施例提供的地震信号处理装置的组成结构示意图，如图9所示，该地震信号处理装置包括：
信息获取模块901，用于获取采样点反射地震信号的偏移距和对应的反射波到时；
公式构建模块902，用于根据采样点反射地震信号的偏移距和对应的反射 波到时，构建基于Pade近似的非双曲动校正公式；
参数提取模块903，位于信息获取模块与公式构建模块之间，用于根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。
进一步的，上述地震信号的处理装置中，所述公式构建模块902包括：
归一化单元，用于对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
Pade式构建单元，用于基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式。
进一步的，上述地震信号的处理装置中，所述参数提取模块903包括：
公式扫描单元，用于由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)；
参数换算单元，用于利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时；
结果求取单元，用于根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
上述方案中，所述地震信号的处理装置还包括：
层剥离模块，用于在多层的各向异性速度分析中对地震信号进行层剥离处理。
具体的，基于Pade近似的非双曲动校正公式，对当前层的记录进行速度独立的层剥离，将多层的反射记录转换为单层的反射记录，然后重复单层模型的分析过程，直到将所有的反射全部处理完，由此就能实现对多层的各向异性速度分析。
本发明实施例还提供一种地震信号处理系统，该系统中包括上述任意一种地震信号处理装置。
本发明实施例提供地震信号处理方法、装置和系统的优势是：本发明利用基于Pade近似的非双曲动校正公式处理地震信号，可以处理地下介质同时具有强各向异性和长偏移距的情形；并且在有强各向异性和长偏移距二者之一的情形时，比现有地震信号的处理方法精度更高。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (10)

  1. 一种地震信号处理方法，其特征在于，所述方法包括：
    获取采样点反射地震信号的偏移距和对应的反射波到时；
    根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式；
    根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。
  2. 根据权利要求1所述的地震信号处理方法，其特征在于，所述根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式包括：
    对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
    基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式。
  3. 根据权利要求2所述的地震信号处理方法，其特征在于，所述基于归一化偏移距和归一化反射波到时的对应关系的Pade近似式包括：
    Figure PCTCN2015085576-appb-100001
    其中，x为归一化偏移距，τ(x)为归一化反射波到时，n为Pade近似式的阶数，Pk与Qk为第k阶待定系数；
  4. 根据权利要求1所述的地震信号处理方法，其特征在于，所述根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数包括：
    由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)；
    利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时；
    根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
  5. 根据权利要求1所述的地震信号处理方法，其特征在于，当地震信号中 包含多个地质层的地震信号时，所述方法还包括：
    在多层的各向异性速度分析中对地震信号进行层剥离处理。
  6. 一种地震信号处理装置，其特征在于，所述装置包括：信息获取模块、公式构建模块、参数提取模块；其中，
    信息获取模块，用于获取采样点反射地震信号的偏移距和对应的反射波到时；
    公式构建模块，用于根据采样点反射地震信号的偏移距和对应的反射波到时，构建基于Pade近似的非双曲动校正公式；
    参数提取模块，用于根据构建的基于Pade近似的非双曲动校正公式提取所述反射地震信号的垂向传播速度和各向异性参数。
  7. 根据权利要求6所述的地震信号处理装置，其特征在于，所述公式构建模块包括：
    归一化单元，用于对获取的采样点反射地震信号的偏移距和反射波到时进行归一化处理；
    Pade式构建单元，用于基于归一化偏移距和归一化反射波到时，构建它们对应关系的Pade近似式。
  8. 根据权利要求6所述的地震信号处理装置，其特征在于，所述参数提取模块包括：
    公式扫描单元，用于由所构建的基于Pade近似的非双曲动校正公式获取归一化反射波到时τ(x)；
    参数换算单元，用于利用归一化偏移距x的定义x＝X/t0υnmo和归一化反射波到时τ(x)的定义τ(x)＝t(X)/t0，换算地震信号实际的偏移距和对应的反射波到时；
    结果求取单元，用于根据换算出的实际的偏移距和对应的反射波到时，求取地震信号对应的垂向传播速度和各向异性参数。
  9. 根据权利要求6所述的地震信号处理装置，其特征在于，所述地震信号的处理装置还包括：
    层剥离模块，用于对地震信号进行层剥离处理。
  10. 一种地震信号的处理系统，其特征在于，该系统中包括根据权利要求5至9任一项所述的地震信号处理装置。
PCT/CN2015/085576 2015-07-30 2015-07-30 一种地震信号处理方法、装置和系统 WO2017015954A1 (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/CN2015/085576 WO2017015954A1 (zh) 2015-07-30 2015-07-30 一种地震信号处理方法、装置和系统
CN201580083576.1A CN108260359A (zh) 2015-07-30 2015-07-30 一种地震信号处理方法、装置和系统
US15/533,247 US20170336523A1 (en) 2015-07-30 2015-07-30 Seismic signal processing method, apparatus and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2015/085576 WO2017015954A1 (zh) 2015-07-30 2015-07-30 一种地震信号处理方法、装置和系统

Publications (1)

Publication Number Publication Date
WO2017015954A1 true WO2017015954A1 (zh) 2017-02-02

Family

ID=57883963

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2015/085576 WO2017015954A1 (zh) 2015-07-30 2015-07-30 一种地震信号处理方法、装置和系统

Country Status (3)

Country Link
US (1) US20170336523A1 (zh)
CN (1) CN108260359A (zh)
WO (1) WO2017015954A1 (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109490962A (zh) * 2018-11-30 2019-03-19 中国海洋大学 一种浅地层剖面非零偏移距的消除方法
CN111624656B (zh) * 2019-02-27 2023-02-10 中国石油天然气集团有限公司 基于近偏移距地震信号的节点时钟漂移校正方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100054082A1 (en) * 2008-08-29 2010-03-04 Acceleware Corp. Reverse-time depth migration with reduced memory requirements
CN102053261A (zh) * 2009-10-29 2011-05-11 中国石油化工股份有限公司 一种地震数据处理方法
CN102073064A (zh) * 2009-11-25 2011-05-25 中国石油天然气集团公司 一种利用相位信息提高速度谱分辨率的方法
CN103257363A (zh) * 2013-05-21 2013-08-21 西南石油大学 一种探测地下裂缝性储层中裂缝倾角的方法
CN104678439A (zh) * 2013-11-27 2015-06-03 中国石油化工股份有限公司 一种多波连续介质速度分析方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7343245B2 (en) * 2005-06-09 2008-03-11 Eseis, Inc. Method of processing seismic data to extract and portray AVO information
CN103105623B (zh) * 2012-12-13 2013-08-21 东北石油大学 一种地震勘探中的数据波形处理方法
CN104155691B (zh) * 2014-08-27 2017-12-26 中国石油天然气集团公司 转换波各向异性速度分析方法及装置
CN104635264B (zh) * 2014-08-28 2017-03-08 中国石油天然气股份有限公司 叠前地震数据的处理方法及设备
CN104407378B (zh) * 2014-11-25 2017-05-10 中国石油天然气股份有限公司 一种各向异性参数反演方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100054082A1 (en) * 2008-08-29 2010-03-04 Acceleware Corp. Reverse-time depth migration with reduced memory requirements
CN102053261A (zh) * 2009-10-29 2011-05-11 中国石油化工股份有限公司 一种地震数据处理方法
CN102073064A (zh) * 2009-11-25 2011-05-25 中国石油天然气集团公司 一种利用相位信息提高速度谱分辨率的方法
CN103257363A (zh) * 2013-05-21 2013-08-21 西南石油大学 一种探测地下裂缝性储层中裂缝倾角的方法
CN104678439A (zh) * 2013-11-27 2015-06-03 中国石油化工股份有限公司 一种多波连续介质速度分析方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
FU, QIANG ET AL.: "Velocity Analysis and Dynamic Correction of P-Wave Based on VTI Medium Theory", COMPUTING TECHNIQUES GEOPHYSICAL AND GEOCHEMICAL EXPLORATION, vol. 30, no. 1, 31 January 2008 (2008-01-31), pages 10 - 16, ISSN: 1001-1749 *
TIAN, CHUNZHI ET AL.: "Study on Seismic Wave Propagation for Deep Layer Seismic Exploration", PROGRESS IN GEOPHYSICS, vol. 17, no. 4, 31 December 2002 (2002-12-31), pages 564 - 574, ISSN: 1004-2903 *

Also Published As

Publication number Publication date
US20170336523A1 (en) 2017-11-23
CN108260359A (zh) 2018-07-06

Similar Documents

Publication Publication Date Title
CA2542577C (en) Radial profiling of slowness: methods and apparatus for near-wellbore alteration estimation
CN106405651B (zh) 一种基于测井匹配的全波形反演初始速度模型构建方法
Ji et al. Deep neural network-based permittivity inversions for ground penetrating radar data
CN103492910B (zh) 时域中的同步小波提取和反卷积
Haldorsen et al. Locating microseismic sources using migration-based deconvolution
RU2631407C1 (ru) Способ и устройство для обработки сейсмических сигналов
CN103852785B (zh) 地层各向异性的评价方法
Duret et al. Near-surface velocity modeling using a combined inversion of surface and refracted P-waves
US6980482B2 (en) Seismic survey system
WO2017015954A1 (zh) 一种地震信号处理方法、装置和系统
Guo et al. Becoming effective velocity-model builders and depth imagers, Part 2—The basics of velocity-model building, examples and discussions
Zhu et al. A multi-channel approach for automatic microseismic event association using RANSAC-based arrival time event clustering (RATEC)
US10126450B2 (en) Black hole boundary conditions
CN116755148A (zh) 正交各向异性介质多方位反射波走时反演方法
CN108693560A (zh) 一种基于互相关道的散射波成像方法及系统
Mildner et al. Target-oriented velocity analysis using Marchenko-redatumed data
Khadhraoui et al. Real-time detection and localization of microseismic events
CN109143398B (zh) 一种自动网格层析深度域速度的建模方法
GB2503640A (en) Quality Assurance in a Full Waveform Inversion Process
Hogarth et al. Controlled-source velocity calibration for real-time downhole microseismic monitoring
CN111474580B (zh) 一种基于炮检距矢量片的方位角道集提取方法和系统
US10338248B2 (en) Method for processing seismic images
CN106338760B (zh) 误差补偿的起伏地表偏移方法
Baziw et al. Deriving interval velocities from downhole seismic data
CN108387937B (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: 15899306

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: 15899306

Country of ref document: EP

Kind code of ref document: A1