CN110703318B - Direct inversion method of pre-stack seismic data - Google Patents

Direct inversion method of pre-stack seismic data Download PDF

Info

Publication number
CN110703318B
CN110703318B CN201910904075.5A CN201910904075A CN110703318B CN 110703318 B CN110703318 B CN 110703318B CN 201910904075 A CN201910904075 A CN 201910904075A CN 110703318 B CN110703318 B CN 110703318B
Authority
CN
China
Prior art keywords
data
mapping
inversion
observation
observation 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.)
Active
Application number
CN201910904075.5A
Other languages
Chinese (zh)
Other versions
CN110703318A (en
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.)
First Institute of Oceanography MNR
Original Assignee
First Institute of Oceanography MNR
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 First Institute of Oceanography MNR filed Critical First Institute of Oceanography MNR
Priority to CN201910904075.5A priority Critical patent/CN110703318B/en
Publication of CN110703318A publication Critical patent/CN110703318A/en
Application granted granted Critical
Publication of CN110703318B publication Critical patent/CN110703318B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • 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/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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

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

According to the method for directly inverting the pre-stack seismic data, the mapping from the observation data to the model parameters is constructed to be used as the inversion operator, and the constructed inversion operator is used for directly inverting the pre-stack seismic data, so that the calculation efficiency is improved, the calculation amount is reduced, and meanwhile, the problems of dependence on the initial model and reduction of the convergence of the inversion result are effectively avoided. The method comprises the following implementation steps: (1) acquiring a limited group of observation data and model parameter data which is less than m and corresponds to the limited group of observation data, wherein D is used as sample data; (2) sorting the obtained sample data; (3) constructing a mapping from observed data to model parameters
Figure DDA0002212744870000011
(4) And directly inverting the data to be inverted by utilizing the constructed mapping relation.

Description

Direct inversion method of pre-stack seismic data
Technical Field
The invention relates to a method for inverting prestack seismic data, in particular to a method for directly inverting after constructing a mapping relation from the prestack seismic data to model parameters by using sample data, and belongs to the field of seismic exploration.
Background
The seismic data prestack inversion is an effective means for solving the physical property parameters of the underground medium based on the change relation of the reflection coefficient along with the incident angle. The prestack inversion technology takes the Zoeppritz (1919) equation and the approximate equation thereof (Bortfeld, 1961; Aki and Richards, 1980; Shuey,1985) as the theoretical basis, has been proposed since the 80 th century (Ostrander,1982,1984), and is continuously developed and plays an important role in the field of oil and gas resource exploration.
In the prestack seismic data, the variation of the seismic reflection coefficient with the incidence angle includes various physical parameters of the stratum, such as: longitudinal wave velocity, shear wave velocity, density, and the like, and combinations thereof. In the pre-stack seismic inversion, observation data are pre-stack gather data which come from the same reflection point and contain reflection coefficient change information along with incident angles, parameters to be inverted are model parameters containing longitudinal wave velocity, transverse wave velocity, density and the like, and the relation is established between the observation data D and a model parameter vector m through a forward operator f:
D＝f(m) (1)
in the above equation (1), the positive operator f is a mapping from the model parameter vector m to the observation data D, and is specifically Zoeppritz (1919) equation in the prestack inversion.
However, the forward operator f is only a mapping from m to D, and the geophysical inversion is a process of solving the model parameter m from the known observation data D, and the process cannot be directly realized by the equation (1). Moreover, the forward modeling process shown in equation (1) is often a multidimensional, highly nonlinear mapping relationship, and a mapping from the observation data D to the model parameter m cannot be directly obtained by inverting f:
m＝f-1(D) (2)
therefore, in the inversion problem, an objective function is constructed through the formula (1), and the optimal model parameters are solved as the inversion result by using a forward iteration optimization method. The procedure is shown in the attached FIG. 1, first, initial model parameters m are given, and simulation data D is obtained by forward calculation (formula (1))synRecording actual observation data as DobsThen the least squares objective function can be written as:
Figure BDA0002212744850000011
determining an objective function
Figure BDA0002212744850000012
Whether or not the minimum value is reached, if
Figure BDA0002212744850000013
If the minimum value is not reached, further adjusting the model parameter m, and obtaining new simulation data D through forward calculation of the adjusted msyn(ii) a Judgment of
Figure BDA0002212744850000014
Whether the minimum value is reached or not, and if the minimum value is reached, outputting m at the moment as an inversion result m'; if the minimum value is not reached, continuing to adjust m, and continuing to iteratively calculate until the minimum value is reached
Figure BDA0002212744850000021
Reaching a minimum value.
The complexity of the forward operator determines the difficulty of inversion, and the optimization solving process based on forward iteration is often large in calculated amount and low in calculation efficiency, and has the problems of large dependence on an initial model and convergence. The exact Zoeppritz equation is also not physically intuitive due to its mathematically complex form.
To overcome the mathematical complexity and physical non-intuitiveness of the Zoeppritz equation, many scholars approximate the equation from different perspectives. The approximation equation, while overcoming the deficiencies of the exact equation, reduces the accuracy of the calculation.
As summarized above, although the precision of the iterative inversion method based on the precise equation is high, the iterative inversion method is large in calculation amount, low in calculation efficiency, large in dependence on an initial model, and has the problem of convergence. While the linear inversion technology based on the approximate equation has high calculation efficiency, the inversion result has low precision.
In view of this, the present patent application is specifically proposed.
Disclosure of Invention
The direct inversion method of the prestack seismic data is designed to solve the problems in the prior art, and the direct inversion method is provided.
The direct inversion method of the pre-stack seismic data mainly comprises the following implementation steps:
(1) acquiring a limited group of observation data and model parameter data which is less than m and corresponds to the limited group of observation data, wherein D is used as sample data;
(2) sorting the obtained sample data;
(3) constructing a mapping from observed data to model parameters
Figure BDA0002212744850000022
(4) And directly inverting the data to be inverted by utilizing the constructed mapping relation.
Further, in the step (1), each set of the primitive data < m, D > needs to satisfy the mapping relationship from the model parameter m to the observation data D within a certain error range, and the sample data can be obtained by theoretical forward calculation or can be obtained from actual observation data. The observation data D is prestack gather data containing the variation of the seismic reflection coefficient with the incidence angle.
In the step (2), the sample data is sorted by adopting data cleaning operations including outlier rejection, noise suppression and repeated sample deletion.
The step (3) specifically includes:
(3-1) regarding the obtained multiple groups of sample data as discrete samples on the mapping F from the observation data to the model parameters to be constructed, and constructing the mapping F from the observation data to the model parameters in the sample data coverage range by utilizing an interpolation or regression method.
(3-2) mapping in step (3)
Figure BDA0002212744850000031
The specific expression of (A) is as follows:
Figure BDA0002212744850000032
wherein,
Figure BDA0002212744850000033
for any pre-stack seismic data to be inverted in the coverage range of a plurality of groups of obtained sample data < m and D,
Figure BDA0002212744850000034
for the inversion result corresponding to the pre-stack data, wiRepresenting their corresponding sub-functions for the weight coefficients
Figure BDA0002212744850000035
The proportion in the mapping.
(3-3) weight coefficient w in step (3-2)iSum subfunction giSamples obtained as described aboveThis data is solved:
Figure BDA0002212744850000036
in the formula, | m-F (D) | represents that the sample data is less than m, D > is the distance between the model parameter m and the observation data D after inversion by mapping F, sigma-represents that the obtained finite group of samples are summed,
Figure BDA0002212744850000037
represents w when the objective function takes the minimum valueiAnd gi。
In summary, compared with the iterative inversion method in the prior art, the direct inversion method has the following advantages:
1. the constructed mapping from the observation data to the model parameters is directly utilized to carry out direct inversion, and an iterative process adopted by solving a target function in the traditional inversion process is not needed, so that the problems of initial values, convergence and the like can be effectively avoided or reduced.
2. The constructed inversion operator is stored for future use, not only for one specific inversion task, thus greatly improving the computational efficiency.
3. On the basis of improving the calculation efficiency, the accuracy of the inversion result cannot be reduced.
Drawings
FIG. 1 is a schematic flow chart of an iterative inversion method adopted in the prior art;
FIG. 2 is a schematic flow diagram of a direct inversion method described herein;
FIG. 3 is a schematic view of the transverse distribution of elastic parameters;
FIG. 4 is a graph of reflection coefficient as a function of angle of incidence;
FIG. 5 is a schematic diagram of the results of a direct inversion;
FIG. 6 is a graph of error comparison of direct inversion and conventional inversion method results;
Detailed Description
Example 1 the present invention will now be described in further detail with reference to the accompanying drawings, which are intended to represent only some, but not all, embodiments of the present application. All other embodiments, which can be obtained by persons skilled in the art without any inventive work, based on the embodiments of the present invention, belong to the scope of the present invention.
As shown in fig. 2, the direct inversion method of pre-stack seismic data mainly includes the following implementation steps:
(1) acquiring a limited group of observation data and model parameter data which is less than m and corresponds to the limited group of observation data, wherein D is used as sample data;
the data of each group is less than m, D is more than the requirement of meeting the mapping relation between the model parameter m and the observation data D within a certain error range, and the sample data can be directly obtained from actual observation data or can be calculated by theoretical forward modeling. In this example, the sample data is obtained by theoretical forward modeling calculation, that is, observation data is obtained by model parameter calculation using the Zoeppritz equation shown in formula (1). Model parameter ml＝[αiβjρk]In this example, the longitudinal wave velocity α of the sea bottomiSea bottom transverse wave velocity betajAnd seafloor density ρkRespectively has a value range of 1530-alphai≤1570、180≤βj≤220、1630≤ρkThe value interval is less than or equal to 1670 and is 2.5. From mlForward calculation to obtain observation data DlWherein the observation data DlThe sea bottom reflection coefficient in the range of 1 deg. to 55 deg. is calculated at intervals of 1 deg..
(2) Sorting the obtained sample data;
the sample data is sorted by adopting data cleaning operations of abnormal value elimination, noise suppression, repeated sample deletion and the like. As can be seen from the above forward calculations, i is 1,2, …,17, j is 1,2, …,17, k is 1,2, …, 17; then the model parameter mlThere are a total of 17 × 17 × 17 — 4913 combinations. We have obtained 4913 sets of non-repetitive sample data altogether<ml,Dl>，l＝1,2,…,4913。
(3) Constructing a mapping from observed data to model parameters
Figure BDA0002212744850000041
(3-1) the obtained multiple groups of sample data < ml,DlRegarding as discrete sampling on the mapping F from the observation data to the model parameters to be constructed, and constructing the mapping F from the observation data to the model parameters in the sample data coverage range by utilizing an interpolation or regression method.
(3-2) mapping to be constructed
Figure BDA0002212744850000042
An equivalent variant written as equation (4) is:
Figure BDA0002212744850000043
in the above formula (6)
Figure BDA0002212744850000044
w＝[w1 w2 … wi … wP]TAnd T denotes transposition.
(3-3) the weight coefficient w in the step (3-2) is obtained by using the formula (5)iSum subfunction giSubstituting all sample data into < ml,DlMore than get:
M＝G·w (7)
in the above formula (7)
Figure BDA0002212744850000045
Selecting
Figure BDA0002212744850000046
To obtain w ═ G-1·M。
To this end, we have obtained the mapping to be constructed
Figure BDA0002212744850000047
(i.e., equation (4)) that the construction of the mapping F from the observation data D to the model parameters m is completed.
(4) And directly inverting the data to be inverted by utilizing the constructed mapping relation.
In order to facilitate the analysis of the accuracy of the inversion result, the data to be inverted used in the embodiment is theoretical forward data, and the longitudinal wave velocity alpha of the sediment at the bottom of the sea is assumed2Transverse wave velocity beta2And density ρ2Distribution along a two-dimensional seismic line, as shown in FIG. 3:
specifically, the elastic parameter distribution is given by equation (8):
Figure BDA0002212744850000051
in the above formula, n is 1,2, …,200 is CRP gather number.
The sea bottom reflection coefficient in the range of 1 ° to 55 ° was calculated at intervals of 1 ° for each reflection point (each CRP), and random errors were added to the obtained reflection coefficients to simulate data containing noise in practice, as shown in fig. 4, which shows the reflection coefficients at CRP numbers 25, 75, and 150, respectively, as a function of the incident angle.
Thus, a two-dimensional line theoretical simulation data comprising 200 CRP gathers was obtained.
The theoretical simulation data used here only replaces the actual observation data, which is convenient for comparing inversion errors; when the actual data is used for inversion, the simulation synthesis of the seismic data is not needed to be carried out firstly, but the actual data is directly obtained and inverted.
Taking the observation data as the mapping from the observation data to the model parameters constructed in step (3)
Figure BDA0002212744850000052
The obtained model parameter values are the direct inversion results, and the inversion results are shown in fig. 5.
In fig. 5, a thick solid line represents model parameters, and a thin solid line represents inversion results when observed data contains noise; the left side shows the results obtained with the direct inversion method proposed by the present invention, and the right side is the results obtained with the conventional inversion. Obviously, under the influence of noise, the inversion result jumps up and down with a given model value as a center, and the error of the transverse wave velocity is the largest among three parameters of the inversion, and the longitudinal wave velocity and the density are the second order, which is consistent with the traditional inversion result.
As shown in fig. 6, a comparison between the inversion error of the direct inversion method proposed in the present application and the inversion error of the conventional inversion method is shown. In the figure, "+" represents the inversion error of the direct inversion method provided by the invention, and a solid line represents the inversion error obtained by the traditional inversion method. The two errors are obviously basically consistent, and the direct inversion method provided by the invention can achieve the inversion accuracy consistent with that of the traditional inversion method.
In summary, the method provided by the present invention has the advantages described above, and the accuracy of the inversion result is not reduced, so that the method provided by the present invention is feasible.

Claims (4)

1. A direct inversion method of pre-stack seismic data is characterized by comprising the following steps: comprises the following implementation steps of the method,
(1) acquiring a limited group of observation data and model parameter data which is less than m and corresponds to the limited group of observation data, wherein D is used as sample data;
(2) sorting the obtained sample data;
(3) constructing a mapping from observed data to model parameters
Figure FDA0002940199720000011
(4) Directly inverting the data to be inverted by utilizing the constructed mapping relation;
the mapping
Figure FDA0002940199720000012
The specific expression of (A) is as follows:
Figure FDA0002940199720000013
wherein,
Figure FDA0002940199720000014
for any pre-stack seismic data to be inverted in the coverage range of a plurality of groups of obtained sample data < m and D,
Figure FDA0002940199720000015
for the inversion result corresponding to the pre-stack data, wiRepresenting their corresponding sub-functions for the weight coefficients
Figure FDA0002940199720000016
The proportion in the mapping;
weight coefficient wiSum subfunction giThe calculation of (1):
Figure FDA0002940199720000017
in the formula, | m-F (D) | represents that the sample data is less than m, D > is the distance between the model parameter m and the observation data D after inversion by mapping F, sigma-represents that the obtained finite group of samples are summed,
Figure FDA0002940199720000018
represents w when the objective function takes the minimum valueiAnd gi。
2. The method of direct inversion of prestack seismic data according to claim 1, characterized by: in the step (1), each group of sample data is less than m, D > needs to satisfy the mapping relationship from the model parameter m to the observation data D within a certain error range, the sample data can be obtained by theoretical forward calculation or obtained from actual observation data, and the observation data D is pre-stack gather data containing the variation relationship of the seismic reflection coefficient with the incident angle.
3. The method of direct inversion of prestack seismic data according to claim 1, characterized by: in the step (2), the sample data is sorted by adopting data cleaning operations including outlier elimination, noise suppression and repeated sample deletion.
4. The method of direct inversion of prestack seismic data according to claim 1, characterized by: in the step (3), the obtained multiple groups of sample data are regarded as discrete samples on the mapping F from the observation data to the model parameters to be constructed, and the mapping F from the observation data to the model parameters in the sample data coverage range is constructed by utilizing an interpolation or regression method.
CN201910904075.5A 2019-09-24 2019-09-24 Direct inversion method of pre-stack seismic data Active CN110703318B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910904075.5A CN110703318B (en) 2019-09-24 2019-09-24 Direct inversion method of pre-stack seismic data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910904075.5A CN110703318B (en) 2019-09-24 2019-09-24 Direct inversion method of pre-stack seismic data

Publications (2)

Publication Number Publication Date
CN110703318A CN110703318A (en) 2020-01-17
CN110703318B true CN110703318B (en) 2021-06-11

Family

ID=69195894

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910904075.5A Active CN110703318B (en) 2019-09-24 2019-09-24 Direct inversion method of pre-stack seismic data

Country Status (1)

Country Link
CN (1) CN110703318B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113569493B (en) * 2021-09-26 2021-12-17 自然资源部第一海洋研究所 Domain decomposition-based physical driving deep learning inversion method

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103097914A (en) * 2010-04-06 2013-05-08 道达尔公司 A process for characterising the evolution of a reservoir
CN107092029A (en) * 2017-04-26 2017-08-25 中国石油大学(北京) A kind of seismic inversion method and device
CN108072892A (en) * 2016-11-09 2018-05-25 中国石油化工股份有限公司 A kind of geological structure constraint chromatography conversion method of automation
CN108139499A (en) * 2015-10-02 2018-06-08 埃克森美孚上游研究公司 The full wave field inversion of Q- compensation

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8923093B2 (en) * 2009-08-25 2014-12-30 Westerngeco L.L.C. Determining the quality of a seismic inversion
US10670751B2 (en) * 2015-03-27 2020-06-02 Cgg Services Sas Full waveform inversion method for seismic data processing using preserved amplitude reverse time migration

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103097914A (en) * 2010-04-06 2013-05-08 道达尔公司 A process for characterising the evolution of a reservoir
CN108139499A (en) * 2015-10-02 2018-06-08 埃克森美孚上游研究公司 The full wave field inversion of Q- compensation
CN108072892A (en) * 2016-11-09 2018-05-25 中国石油化工股份有限公司 A kind of geological structure constraint chromatography conversion method of automation
CN107092029A (en) * 2017-04-26 2017-08-25 中国石油大学(北京) A kind of seismic inversion method and device

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
非高斯AVO三参数反演算法及其应用研究;刘洋;《中国博士学位论文全文数据库 基础科学辑》;20151115(第11期);第24页 *

Also Published As

Publication number Publication date
CN110703318A (en) 2020-01-17

Similar Documents

Publication Publication Date Title
CN106526674B (en) Three-dimensional full waveform inversion energy weighting gradient preprocessing method
US8406081B2 (en) Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
CN110780351B (en) Longitudinal wave and converted wave prestack joint inversion method and system
AU2012220584A1 (en) Sensitivity kernel-based migration velocity analysis in 3D anisotropic media
CN111025387B (en) Pre-stack earthquake multi-parameter inversion method for shale reservoir
CN109521469B (en) Regularization inversion method for elastic parameters of submarine sediments
US9952341B2 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
CN112698390B (en) Pre-stack seismic inversion method and device
CN105629299A (en) Travel-time table and angle table acquisition method for angle domain prestack depth migration and imaging method
IL92132A (en) Homeomorphical imaging method of analyzing the structure of a medium
WO2021127382A1 (en) Full waveform inversion in the midpoint-offset domain
CN111025388B (en) Multi-wave combined prestack waveform inversion method
CN110703318B (en) Direct inversion method of pre-stack seismic data
CN111175821B (en) Step-by-step inversion method for anisotropic parameters of VTI medium
CN107976714B (en) A kind of channel set calculation method of complicated earth surface classification spatial distance weighting
CN111538084B (en) Method and system for converting OVT domain data into azimuth angle domain imaging gather
CN117388944A (en) Multi-physical parameter inversion method of geologic model constraint
TANG et al. Multiple Ray Tracing Within 3‐D Layered Media with the Shortest Path Method
CN110208858B (en) Method and system for directly estimating &#39;sweet spot&#39; probability based on pre-stack inversion
CN112083492A (en) Full-path compensation primary wave and multiple wave combined imaging method under deep sea environment
CN112526605B (en) Method for simulating and exploring natural gas hydrate by adopting seismic numerical value
CN115184986B (en) Global envelope cross-correlation full waveform inversion method independent of seismic source
NL2031350B1 (en) Method and system for determining geosteering irregular observation system
CN118244355B (en) Reflection wave travel time inversion method based on reconstructed observation seismic data
Purba et al. Improving the accuracy of the expanded anisotropic eikonal equation at larger offsets using Levin T-transformation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant