CN111505713A - Pre-stack seismic inversion method based on multi-point geological statistics - Google Patents
Pre-stack seismic inversion method based on multi-point geological statistics Download PDFInfo
- Publication number
- CN111505713A CN111505713A CN202010068941.4A CN202010068941A CN111505713A CN 111505713 A CN111505713 A CN 111505713A CN 202010068941 A CN202010068941 A CN 202010068941A CN 111505713 A CN111505713 A CN 111505713A
- Authority
- CN
- China
- Prior art keywords
- seismic
- data
- inversion
- estimated
- grid
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 30
- 239000011435 rock Substances 0.000 claims abstract description 20
- 238000012549 training Methods 0.000 claims abstract description 15
- 208000035126 Facies Diseases 0.000 claims description 14
- 238000004088 simulation Methods 0.000 claims description 14
- 238000005070 sampling Methods 0.000 claims description 13
- 238000005315 distribution function Methods 0.000 claims description 9
- 239000002131 composite material Substances 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 230000001186 cumulative effect Effects 0.000 claims description 6
- 230000000877 morphologic effect Effects 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 6
- 238000012804 iterative process Methods 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 238000012512 characterization method Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000010587 phase diagram Methods 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
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/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/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
- 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/624—Reservoir parameters
-
- 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
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
The invention discloses a prestack seismic inversion method based on multipoint geological statistics, which comprises the following steps of ①, investigating rock phase and physical attribute statistical relationship, ②, establishing a work area grid and a training image conforming to a work area, ③, distributing measured well data to the nearest grid node, ④, assigning a work area grid initial attribute value, ⑤, performing iterative inversion and ⑥, and calculating seismic record matching rate.
Description
Technical Field
The invention relates to the technical field of oil and gas exploration and development, in particular to a prestack seismic inversion method based on multi-point geological statistics.
Background
At the middle and later stages of oil field exploration and development, a method based on two-point geostatistics inversion becomes a mainstream seismic reservoir inversion technology.
In the traditional two-point geostatistical inversion (Deutsch, 1992; Hass and Dubrule, 1994; Yanze 2012, etc.), the reservoir spatial correlation is characterized mainly by a variogram. However, the two-point geological statistical method only considers the correlation between two spatial points, and is difficult to be competent for characterization of a reservoir with a complex morphology (Yi Yan tree, etc., 2009, 2011; Yangbeigen, etc., 2013, 2014).
In order to better characterize the complex change characteristics of the reservoir morphology, the joint distribution of a plurality of spatial points needs to be considered. Gonzalezet al (2008) for the first time attempted to apply multi-point geostatistics to reservoir inversion, indicating a good prospect for the application of multi-point geostatistics in inversion,
however, the influence of the result of the prior iterative inversion on the later iterative prediction cannot be reflected in the later iterative process, and the inversion randomness and the calculation load are increased.
Later multipoint inversion studies are increasing (Cheolkyun Jeong, 2017; Liuxing, 2018; Yi Yan Tree, 2018, etc.). Although multipoint inversion is greatly improved in the aspect of reservoir morphology characterization compared with two-point inversion, most of the existing multipoint geostatistical inversion is post-stack inversion, and seismic data (such as transverse wave information) are not fully utilized, and the prediction accuracy of the pre-stack inversion to the reservoir is higher than that of the post-stack inversion, so that the pre-stack inversion becomes a consensus.
Therefore, a method for improving the reservoir prediction accuracy is urgently needed to solve the above problems.
Disclosure of Invention
The invention aims to overcome the defects of the background technology and provides a prestack seismic inversion method based on multi-point geological statistics.
The purpose of the invention is implemented by the following technical scheme: the prestack seismic inversion method based on the multipoint geological statistics is characterized by comprising the following steps: it comprises the following steps:
① investigating rock phase and physical attribute statistical relationship
Acquiring elastic parameters of different rock types, such as density, longitudinal wave velocity and transverse wave velocity, according to well data interpretation results; building cumulative distribution maps of different rock elastic parameters; the distribution function of the parameters is used as the basis for sampling the elastic parameters of the rock under the later-stage multi-point geostatistical phased modeling;
② creating a work area grid and training images conforming to the work area
Selecting a proper grid size according to the actual work area range, carrying out grid division on the work area, and establishing a grid model; obtaining rock facies morphological characteristics and a sedimentary model of a research area through basic geological research, wherein the commonly used obtaining method comprises unconditional modeling based on a target, modeling based on a sedimentary process, seismic data, geological personnel sketching geological map data and the like;
③ assigning measured well data to nearest neighbor grid nodes
The measured well data comprises lithofacies, density and longitudinal and transverse wave information, and the data are distributed to the nearest grid nodes as hard data; when a plurality of well data exist in one grid, selecting a data point nearest to the center of the grid as hard data;
④ assigning initial attribute value of work area grid
In order to ensure smooth calculation of the seismic channel synthetic record in the iterative process, initial attribute values including a density value, a longitudinal wave velocity value and a transverse wave velocity value are required to be given to a point to be estimated in a work area; referring to the distribution function of the elastic parameters of different lithofacies of the work area counted in the step 1, if points to be estimated are all given average density values and average longitudinal and transverse wave velocity values;
⑤ iterative inversion
Each iterative inversion needs to traverse all grids, so that the inversion result of the work area is composed of the inversion results of a plurality of positions to be estimated.
⑥ calculating seismic record matching rate
Calculating the seismic matching rate of the current iterative inversion, and comparing the seismic matching rate with a manually set threshold value; when the iterative inversion result is smaller than the seismic record matching rate, repeating the step 4 and the step 5, and performing the next iterative inversion; and when the iterative inversion result is greater than or equal to the seismic record matching rate, ending the whole inversion process to obtain a multi-point geostatistical prestack seismic inversion model.
In the above technical solution, the iterative inversion includes the following sub-steps in step ⑤;
step 5.1: selecting a data template of appropriate size
Determining the shape and size of the data template according to the morphological characteristics of the sedimentary facies;
step 5.2: determining a pseudo-random path for accessing a trellis to be estimated
Firstly, generating a pseudo-random number in each grid random on a top surface grid, and then searching the number of surrounding condition points at each point to be estimated through a data template; the numerical value of each point to be estimated is equal to the sum of the pseudo random number and the number of the condition points, all the points to be estimated are sorted from large to small, and the grid with a large value is simulated preferentially, so that the simulation path of the whole work area is obtained;
step 5.3: data event scanning training image obtaining candidate data pattern
When a point to be estimated is simulated, a data event needs to be determined; randomly scanning a training image by using a data event to obtain a plurality of completely matched data patterns; selecting a plurality of data modes as candidate data modes according to the sequence;
step 5.4: inner loop selects the attribute closest to the original composite record
Sampling elastic parameters of all points to be estimated in the candidate data mode according to the elastic parameter (density, longitudinal wave velocity and transverse wave velocity) distribution functions of different lithofacies obtained in the step ①;
calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of objective functions of a plurality of synthetic seismic channels and actual seismic records, assigning rock facies, velocities and density values contained in corresponding candidate data modes to corresponding simulation areas, and completing inversion of a grid to be estimated;
and sampling the elastic parameters of all points to be estimated in the candidate data mode again. Calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of target functions of a plurality of synthetic seismic traces and actual seismic records, comparing the target function with the target function calculated last time, and selecting the minimum error as an optimal inversion result; repeating the step to perform a plurality of samples, and selecting the smallest sample as an optimal inversion result;
sampling other candidate data modes, repeating the step 5.4, similarly, comparing the target function of each time with the target function calculated last time, and selecting the minimum value as the optimal inversion result;
and (5.1) repeating the steps of 5.1-5.4 until all grids are traversed, namely all the positions to be estimated are simulated, and realizing one-time iterative inversion.
In the technical scheme, in the substep 5.1, the shape of the data template is an ellipsoid or a cuboid, and the cuboid is a lattice cuboid of 5 × 5 × 7.
In the above technical scheme: in substep 5.4;
calculating the sum of target functions of a plurality of synthetic seismic channels and actual seismic records, and selecting the minimum error as an optimal inversion result; the sum of the objective functions of the synthetic seismic traces and the actual seismic records may be calculated using the following formula:
wherein L represents the seismic trace in which the plurality of candidate data patterns are located, l represents the seismic trace in which the several candidate data patterns are located, and d(l) obsRepresenting actual seismic trace data; d(l) synRepresenting a seismic synthetic record calculated from a convolution formula by means of the estimated elastic parameters and reflection coefficients; psi is the sum of the objective functions of the synthetic seismic traces and the actual seismic records;
the seismic synthetic record is calculated according to a convolution formula through the estimated elastic parameters and the estimated reflection coefficients, and the convolution formula is expressed by the following formula:
dsyn＝w*R (2)
wherein w represents a seismic wavelet; r represents a seismic reflection coefficient; dsynA composite record representing the seismic traces;
the calculated seismic reflection coefficient can be approximated by the equation of Levopritz (Zoeppritz) proposed by Aki & Richard:
in the formula Rpp(θ) represents a longitudinal wave reflection coefficient at an adjacent horizon interface position in the target work area; vp is the longitudinal wave velocity; vs is the longitudinal wave velocity; ρ is the density.
In the above technical scheme: in substep 5.4; in the inversion, the more well data known to be contained in the data template, the earlier the access order, typically a gradual transition from a well data rich region to a well data poor region, and finally to a well free region.
In the above technical scheme: in substep 5.3; the data event is a term which can show the correlation among multiple points in space in multi-point geology and indicates lithofacies data in a data template range; the lithofacies data comprises known well data and lithofacies data of simulated points and the geometrical shape (data configuration) of the data space;
in the above technical solution, in step ⑥, when the seismic matching rate is greater than or equal to 80%, the actual requirement can be met, the inversion is finished, and the simulation result is output.
The invention has the following advantages: 1. according to the method, multiple seismic records are compared when data events are selected in multi-point modeling, and an integral assignment means is adopted when the data events are assigned to a grid, so that the inversion prediction precision is improved.
2. In the iterative process, the optimal synthetic seismic record obtained from the current simulation result is compared with the synthetic seismic record from the previous iterative simulation result to determine whether to replace the previous data pattern (depositional facies) with the currently selected data pattern (depositional facies), thereby achieving iterative updating of the depositional facies model.
Because each updating is to complete the updating of the sedimentary facies and the seismic attributes at the same time, and the requirement of being closer to the actual seismic record is met, namely the error is the minimum, the reservoir inversion consistency is better ensured, the reservoir inversion precision is effectively improved, and a method and a technical guarantee are provided for accurately predicting the reservoir parameters in the oil and gas exploration and development.
3. According to the method, a part of data modes are selected to simulate the grid to be estimated, so that the operation efficiency is improved, and the continuity and accuracy of inversion results are improved.
Drawings
FIG. 1 is a process diagram of example 1 of the present invention;
FIG. 2 is well condition data for example 1 of the present invention;
FIG. 3 is a density cumulative distribution chart of example 1 of the present invention;
FIG. 4 is a longitudinal wave cumulative distribution diagram in example 1 of the present invention;
FIG. 5 is a distribution diagram of the cumulative transverse wave in example 1 of the present invention;
FIG. 6 is a training image according to example 1 of the present invention;
FIG. 7 shows data templates in example 1 of the present invention;
FIG. 8 is a data event of embodiment 1 of the present invention;
FIG. 9 shows data patterns in example 1 of the present invention;
FIG. 10 is a partial view of a pseudo random number according to embodiment 1 of the present invention;
FIG. 11 is a phase diagram of a first iterative inversion in accordance with example 1 of the present invention;
FIG. 12 is a density map of a first iterative inversion according to example 1 of the present invention;
FIG. 13 is a longitudinal velocity map of the first iterative inversion in accordance with example 1 of the present invention;
FIG. 14 is a first iterative inversion shear velocity plot of example 1 in accordance with the present invention;
FIG. 15 is a first iterative inversion synthetic seismic record of example 1 of the present invention.
FIG. 16 is a sixth iterative inversion phase diagram according to example 1 of the present invention;
FIG. 17 is a sixth iterative inversion density plot of example 1 in accordance with the present invention;
FIG. 18 is a sixth iteration inversion compressional velocity plot of example 1 in accordance with the present invention;
FIG. 19 is a sixth iteration inversion shear velocity plot of example 1 in accordance with the present invention;
FIG. 20 is a sixth iterative inversion synthetic seismic record of example 1 of the present invention.
Detailed Description
The embodiments of the present invention will be described in detail below with reference to the drawings, but the embodiments are not intended to limit the present invention and are merely examples. The invention is described in further detail by way of example and specific embodiments.
Referring to FIGS. 1-18: the prestack seismic inversion method based on the multipoint geological statistics comprises the following steps:
① investigating rock phase and physical attribute statistical relationship
Acquiring elastic parameters of different rock types, such as density, longitudinal wave velocity and transverse wave velocity, according to well data interpretation results; building cumulative distribution maps of different rock elastic parameters; the distribution function of the parameters is used as the basis for sampling the elastic parameters of the rock under the later-stage multi-point geostatistical phased modeling;
② creating a work area grid and training images conforming to the work area
Selecting a proper grid size according to the actual work area range, carrying out grid division on the work area, and establishing a grid model; obtaining rock facies morphological characteristics and a sedimentary model of a research area through basic geological research, wherein the commonly used obtaining method comprises unconditional modeling based on a target, modeling based on a sedimentary process, seismic data, geological personnel sketching geological map data and the like;
③ assigning measured well data to nearest neighbor grid nodes
The measured well data comprises lithofacies, density and longitudinal and transverse wave information, and the data are distributed to the nearest grid nodes as hard data; when a plurality of well data exist in one grid, selecting a data point nearest to the center of the grid as hard data;
④ assigning initial attribute value of work area grid
In order to ensure smooth calculation of the seismic channel synthetic record in the iterative process, initial attribute values including a density value, a longitudinal wave velocity value and a transverse wave velocity value are required to be given to a point to be estimated in a work area; referring to the distribution function of the elastic parameters of different lithofacies of the work area counted in the step 1, if points to be estimated are all given average density values and average longitudinal and transverse wave velocity values;
⑤ iterative inversion
Each iterative inversion needs to traverse all grids, so that the inversion result of the work area is composed of the inversion results of a plurality of positions to be estimated.
⑥ calculating seismic record matching rate
Calculating the seismic matching rate of the current iterative inversion, and comparing the seismic matching rate with a manually set threshold value; when the iterative inversion result is smaller than the seismic record matching rate, repeating the step 4 and the step 5, and performing the next iterative inversion; and when the iterative inversion result is greater than or equal to the seismic record matching rate, ending the whole inversion process to obtain a multi-point geostatistical prestack seismic inversion model.
The iterative inversion, in step ⑤, includes the following substeps;
step 5.1: selecting a data template of appropriate size
Determining the shape and size of the data template according to the morphological characteristics of the sedimentary facies;
step 5.2: determining a pseudo-random path for accessing a trellis to be estimated
Firstly, generating a pseudo-random number in each grid random on a top surface grid, and then searching the number of surrounding condition points at each point to be estimated through a data template; the numerical value of each point to be estimated is equal to the sum of the pseudo random number and the number of the condition points, all the points to be estimated are sorted from large to small, and the grid with a large value is simulated preferentially, so that the simulation path of the whole work area is obtained;
step 5.3: data event scanning training image obtaining candidate data pattern
When a point to be estimated is simulated, a data event needs to be determined; randomly scanning a training image by using a data event to obtain a plurality of completely matched data patterns; selecting a plurality of data modes as candidate data modes according to the sequence;
step 5.4: inner loop selects the attribute closest to the original composite record
Sampling elastic parameters of all points to be estimated in the candidate data mode according to the elastic parameter (density, longitudinal wave velocity and transverse wave velocity) distribution functions of different lithofacies obtained in the step ①;
calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of objective functions of a plurality of synthetic seismic channels and actual seismic records, assigning rock facies, velocities and density values contained in corresponding candidate data modes to corresponding simulation areas, and completing inversion of a grid to be estimated;
and sampling the elastic parameters of all points to be estimated in the candidate data mode again. Calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of target functions of a plurality of synthetic seismic traces and actual seismic records, comparing the target function with the target function calculated last time, and selecting the minimum error as an optimal inversion result; repeating the step to perform a plurality of samples, and selecting the smallest sample as an optimal inversion result;
sampling other candidate data modes, repeating the step 5.4, similarly, comparing the target function of each time with the target function calculated last time, and selecting the minimum value as the optimal inversion result;
and (5.1) repeating the steps of 5.1-5.4 until all grids are traversed, namely all the positions to be estimated are simulated, and realizing one-time iterative inversion.
In the substep 5.1, the shape of the data template is an ellipsoid or a cuboid, and the cuboid is a lattice cuboid of 5 × 5 × 7.
In substep 5.4; calculating the sum of target functions of a plurality of synthetic seismic channels and actual seismic records, and selecting the minimum error as an optimal inversion result; the sum of the objective functions of the synthetic seismic traces and the actual seismic records may be calculated using the following formula:
wherein L represents the seismic trace in which the plurality of candidate data patterns are located, l represents the seismic trace in which the several candidate data patterns are located, and d(l) obsRepresenting actual seismic trace data; d(l) synRepresenting a seismic synthetic record calculated from a convolution formula by means of the estimated elastic parameters and reflection coefficients; psi is the sum of the objective functions of the synthetic seismic traces and the actual seismic records;
the seismic synthetic record is calculated according to a convolution formula through the estimated elastic parameters and the estimated reflection coefficients, and the convolution formula is expressed by the following formula:
dsyn＝w*R (2)
wherein w represents a seismic wavelet; r represents a seismic reflection coefficient; dsynA composite record representing the seismic traces;
the calculated seismic reflection coefficient can be approximated by the equation of Levopritz (Zoeppritz) proposed by Aki & Richard:
in the formula Rpp(θ) represents a longitudinal wave reflection coefficient at an adjacent horizon interface position in the target work area; vp is the longitudinal wave velocity; vs is the longitudinal wave velocity; ρ is the density.
In substep 5.4; in the inversion, the more well data known to be contained in the data template, the earlier the access order, typically a gradual transition from a well data rich region to a well data poor region, and finally to a well free region.
In substep 5.3; the data event is a term which can show the correlation among multiple points in space in multi-point geology and indicates lithofacies data in a data template range; the lithofacies data comprises known well data and lithofacies data of simulated points and the geometrical shape (data configuration) of the data space;
in step ⑥, when the seismic matching rate is greater than or equal to 80%, the actual requirement can be met, the inversion is finished, and the simulation result is output.
Example 1
Fig. 1 is a diagram of implementation steps of the seismic reservoir inversion method in this embodiment, specifically:
(1) statistical relationship between investigation rock phase and physical property
The sedimentary facies of the work area comprise mudstone and sandstone, and an accumulation distribution map of the sand mudstone is established. As shown in FIG. 2, FIG. 3 and FIG. 4, the density distribution range of mudstone is 2.17-2.33g/cm3The longitudinal wave velocity distribution range is 3930-4272m/s, and the transverse wave velocity distribution range is 2226-2496 m/s; the sandstone density distribution range is 2.27-2.42g/cm3The longitudinal velocity distribution range is 4041-4528m/s, the transverse velocity distribution range is 2300-2561m/s, both the two are in accordance with Gaussian distribution, and the elastic properties of the sand-shale are overlapped and crossed.
(2) Establishing work area grids and training images conforming to work areas
Total length of work area 2000m, total thickness 100m, model mesh divided into 200 × 1 × 45, width of each mesh 10 m, thickness of each mesh about 2.2m, training image as shown in fig. 5 was obtained according to geological features of the work area.
(3) Assigning measured well data to nearest grid nodes
As shown in fig. 6, the phase data of the live log is distributed into the established work area grid. Similarly, statistical density and compressional and shear data in the well data are also assigned to the work area grid.
(4) Giving initial attribute values to the simulation work area
Giving an initial attribute value to unknown points in a work area, and the density of the unknown points is 2.4g/cm3The longitudinal wave velocity was 4200m/s, and the transverse wave velocity was 2600 m/s.
(5) Determining a pseudo-random path for accessing a trellis to be estimated
Firstly, randomly generating a string of pseudo random numbers for all grids to be estimated in a work area, then determining the number of condition points around a data sample plate with each unknown point as a center, then calculating the value of each unknown point, adding random numbers to calculate according to the number of the condition points, and finally sequencing the condition points from large to small to obtain the simulation path of the work area. Taking the grid of the sample size range of the data around the grid points (46, 1, 8) as an example, the first column is the well data and does not participate in the assignment and sorting of the random numbers, the integers of the 2 nd and 3 rd columns from left to right are 7, which indicates that each sample range of the data around the grid contains 7 well data condition points, and similarly, the integer part of the rightmost column is 0, which indicates that no condition point is contained in each sample range of the data around the grid.
(6) Selecting a data template of appropriate size
A data template of 5 × 1 × 7 as shown in fig. 7 was selected.
(7) Simulation of random scanning training image of data event to obtain candidate data
Take the first grid point (46, 1, 8) as an example, wherein the grid point (46, 1, 8) corresponds to the center point of the data template. The data event shown in fig. 8 is obtained by using the data template, the training image is scanned to obtain a series of matched data patterns, 50 data patterns are extracted as the data patterns to be selected, and fig. 9 shows 2 data patterns.
(8) Inner loop selects the attribute closest to the original composite record
Extracting 25 attributes aiming at a mode kernel of each data mode to be selected according to the established elastic parameter distribution, using standard Ricker wavelets with a main frequency of 25HZ, a total time length of 100ms and a sampling rate of 2ms as seismic wavelets, calculating a reflection coefficient according to a formula (1), and calculating seismic records by using the reflection coefficient and seismic wavelet convolution as actual seismic records. The synthetic seismic record is compared to the original seismic record and an optimal one of the modes is selected for retention with its corresponding attribute. And turning to the next grid node until all grid nodes are visited. FIGS. 9, 10, 11, 12 and 13 are the phase, density, longitudinal, shear and synthetic seismogram plots, respectively, for the first iterative inversion.
(8) And establishing a pseudo-random path and starting subsequent iterative inversion. The phase, density, compressional, shear and synthetic seismic traces of the 6 th iterative inversion are shown in figures 14, 15, 16, 17 and 18.
(9) Outputting inversion results
The above-mentioned parts not described in detail are prior art.
Claims (7)
1. The prestack seismic inversion method based on the multipoint geological statistics is characterized by comprising the following steps: it comprises the following steps:
① investigating rock phase and physical attribute statistical relationship
Acquiring elastic parameters of different rock types, such as density, longitudinal wave velocity and transverse wave velocity, according to well data interpretation results; building cumulative distribution maps of different rock elastic parameters; the distribution function of the parameters is used as the basis for sampling the elastic parameters of the rock under the later-stage multi-point geostatistical phased modeling;
② creating a work area grid and training images conforming to the work area
Selecting a proper grid size according to the actual work area range, carrying out grid division on the work area, and establishing a grid model; obtaining rock facies morphological characteristics and a sedimentary model of a research area through basic geological research, wherein the commonly used obtaining method comprises unconditional modeling based on a target, modeling based on a sedimentary process, seismic data, geological personnel sketching geological map data and the like;
③ assigning measured well data to nearest neighbor grid nodes
The measured well data comprises lithofacies, density and longitudinal and transverse wave information, and the data are distributed to the nearest grid nodes as hard data; when a plurality of well data exist in one grid, selecting a data point nearest to the center of the grid as hard data;
④ assigning initial attribute value of work area grid
In order to ensure smooth calculation of the seismic channel synthetic record in the iterative process, initial attribute values including a density value, a longitudinal wave velocity value and a transverse wave velocity value are required to be given to a point to be estimated in a work area; referring to the distribution function of the elastic parameters of different lithofacies of the work area counted in the step 1, if points to be estimated are all given average density values and average longitudinal and transverse wave velocity values;
⑤ iterative inversion
Each iterative inversion needs to traverse all grids, so that the inversion result of the work area is composed of the inversion results of a plurality of positions to be estimated.
⑥ calculating seismic record matching rate
Calculating the seismic matching rate of the current iterative inversion, and comparing the seismic matching rate with a manually set threshold value; when the iterative inversion result is smaller than the seismic record matching rate, repeating the step 4 and the step 5, and performing the next iterative inversion; and when the iterative inversion result is greater than or equal to the seismic record matching rate, ending the whole inversion process to obtain a multi-point geostatistical prestack seismic inversion model.
2. The method of pre-stack seismic inversion based on multi-point geostatistics of claim 1, characterized in that said iterative inversion comprises the following substeps in step ⑤;
step 5.1: selecting a data template of appropriate size
Determining the shape and size of the data template according to the morphological characteristics of the sedimentary facies;
step 5.2: determining a pseudo-random path for accessing a trellis to be estimated
Firstly, generating a pseudo-random number in each grid random on a top surface grid, and then searching the number of surrounding condition points at each point to be estimated through a data template; the numerical value of each point to be estimated is equal to the sum of the pseudo random number and the number of the condition points, all the points to be estimated are sorted from large to small, and the grid with a large value is simulated preferentially, so that the simulation path of the whole work area is obtained;
step 5.3: data event scanning training image obtaining candidate data pattern
When a point to be estimated is simulated, a data event needs to be determined, namely a data template contains condition point information; randomly scanning a training image by using a data event to obtain a plurality of completely matched data patterns; selecting a plurality of data modes as candidate data modes according to the sequence;
step 5.4: inner loop selects the attribute closest to the original composite record
Sampling elastic parameters of all points to be estimated in the candidate data mode according to the elastic parameter (density, longitudinal wave velocity and transverse wave velocity) distribution functions of different lithofacies obtained in the step ①;
calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of objective functions of a plurality of synthetic seismic channels and actual seismic records, assigning rock facies, velocities and density values contained in corresponding candidate data modes to corresponding simulation areas, and completing inversion of a grid to be estimated;
and sampling the elastic parameters of all points to be estimated in the candidate data mode again. Calculating the reflection coefficients of the seismic channels according to the elastic parameter data of the seismic channels where the candidate data patterns are located, and performing convolution on the reflection coefficients and the seismic wavelets to obtain a plurality of synthetic seismic channels; calculating the sum of target functions of a plurality of synthetic seismic traces and actual seismic records, comparing the target function with the target function calculated last time, and selecting the minimum error as an optimal inversion result; repeating the step to perform a plurality of samples, and selecting the smallest sample as an optimal inversion result;
sampling other candidate data modes, repeating the step 5.4, similarly, comparing the target function of each time with the target function calculated last time, and selecting the minimum value as the optimal inversion result;
and (5.1) repeating the steps of 5.1-5.4 until all grids are traversed, namely all the positions to be estimated are simulated, and realizing one-time iterative inversion.
3. The multi-point geostatistical-based prestack seismic inversion method of claim 1 or 2, characterized in that in substep 5.1, the data template is in the shape of an ellipsoid or a cuboid, and the cuboid is a lattice cuboid of 5 × 5 × 7.
4. The multi-point geostatistical-based prestack seismic inversion method of claim 3, characterized in that: in substep 5.4;
calculating the sum of target functions of a plurality of synthetic seismic channels and actual seismic records, and selecting the minimum error as an optimal inversion result; the sum of the objective functions of the synthetic seismic traces and the actual seismic records may be calculated using the following formula:
wherein L represents the seismic trace in which the plurality of candidate data patterns are located, l represents the seismic trace in which the several candidate data patterns are located, and d(l) obsRepresenting actual seismic trace data; d(l) synRepresenting a seismic synthetic record calculated from a convolution formula by means of the estimated elastic parameters and reflection coefficients; psi is the sum of the objective functions of the synthetic seismic traces and the actual seismic records;
the seismic synthetic record is calculated according to a convolution formula through the estimated elastic parameters and the estimated reflection coefficients, and the convolution formula is expressed by the following formula:
dsyn＝w*R (2)
wherein w represents a seismic wavelet; r represents a seismic reflection coefficient; dsynA composite record representing the seismic traces;
the calculated seismic reflection coefficient can be approximated by the equation of Levopritz (Zoeppritz) proposed by Aki & Richard:
in the formula Rpp(θ) represents a longitudinal wave reflection coefficient at an adjacent horizon interface position in the target work area; vp is the longitudinal wave velocity; vs is the longitudinal wave velocity; ρ is the density.
5. The multi-point geostatistical-based prestack seismic inversion method of claim 4, characterized in that: in substep 5.4; in the inversion, the more well data known to be contained in the data template, the earlier the access order, typically a gradual transition from a well data rich region to a well data poor region, and finally to a well free region.
6. The multi-point geostatistical-based prestack seismic inversion method of claim 4 or 5, characterized in that: in substep 5.3; the data event is a term which can show the correlation among multiple points in space in multi-point geology and indicates lithofacies data in a data template range; the lithofacies data comprises known well data and lithofacies data of simulated points and the geometrical shape (data configuration) of the data space;
7. the pre-stack seismic inversion method based on the multi-point geological statistics is characterized in that in step ⑥, when the seismic matching rate is greater than or equal to 80%, the actual requirements can be met, inversion is finished, and a simulation result is output.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010068941.4A CN111505713B (en) | 2020-01-21 | 2020-01-21 | Pre-stack seismic inversion method based on multi-point geological statistics |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010068941.4A CN111505713B (en) | 2020-01-21 | 2020-01-21 | Pre-stack seismic inversion method based on multi-point geological statistics |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111505713A true CN111505713A (en) | 2020-08-07 |
| CN111505713B CN111505713B (en) | 2021-05-07 |
Family
ID=71868921
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010068941.4A Active CN111505713B (en) | 2020-01-21 | 2020-01-21 | Pre-stack seismic inversion method based on multi-point geological statistics |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111505713B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112014882A (en) * | 2020-09-01 | 2020-12-01 | 电子科技大学 | Pre-stack seismic reflection pattern analysis method based on tensor discrimination dictionary |
| CN115963548A (en) * | 2023-01-16 | 2023-04-14 | 中国矿业大学 | Mine microseismic P-wave arrival time picking model construction method based on inverse deductive learning |
| US11686870B2 (en) | 2021-09-29 | 2023-06-27 | China University Of Petroleum (East China) | Interpretive-guided velocity modeling seismic imaging method and system, medium and device |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090164182A1 (en) * | 2007-12-21 | 2009-06-25 | Schlumberger Technology Corporation | Multipoint geostatistics method using branch runlength compression and local grid transformation |
| US20130262051A1 (en) * | 2012-04-03 | 2013-10-03 | Brandon David Plost | Ordered multipoint geostatistics simulation using non-symmetric search mask |
| CN104297787A (en) * | 2014-10-17 | 2015-01-21 | 中国石油天然气股份有限公司 | Method and device for processing three-dimensional lithofacies data of fluvial-facies hypotonic compact sandstone reservoir |
| CN104850682A (en) * | 2015-04-17 | 2015-08-19 | 长江大学 | Multiple-point geostatistics modeling method based on position |
| CN107238862A (en) * | 2016-03-29 | 2017-10-10 | 中国石油化工股份有限公司 | Reflectance factor method of estimation and device based on Bayes's inverting framework |
| CN108645994A (en) * | 2018-04-25 | 2018-10-12 | 中国石油大学(北京) | A kind of geology stochastic inversion methods and device based on Multiple-Point Geostatistics |
| CN108931811A (en) * | 2018-05-17 | 2018-12-04 | 长江大学 | Seismic Reservoir inversion method based on multiple spot geological statistics |
| CN109116421A (en) * | 2018-10-23 | 2019-01-01 | 中海石油(中国)有限公司 | A kind of braid deltas reservoir statistics inverted parameters determine method |
| CN110031896A (en) * | 2019-04-08 | 2019-07-19 | 中国石油天然气集团有限公司 | Earthquake stochastic inversion methods and device based on Multiple-Point Geostatistics prior information |
| CN110031895A (en) * | 2019-03-11 | 2019-07-19 | 西安科技大学 | A kind of Multiple-Point Geostatistics stochastic inversion methods and device based on image stitching |
-
2020
- 2020-01-21 CN CN202010068941.4A patent/CN111505713B/en active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090164182A1 (en) * | 2007-12-21 | 2009-06-25 | Schlumberger Technology Corporation | Multipoint geostatistics method using branch runlength compression and local grid transformation |
| US20130262051A1 (en) * | 2012-04-03 | 2013-10-03 | Brandon David Plost | Ordered multipoint geostatistics simulation using non-symmetric search mask |
| CN104297787A (en) * | 2014-10-17 | 2015-01-21 | 中国石油天然气股份有限公司 | Method and device for processing three-dimensional lithofacies data of fluvial-facies hypotonic compact sandstone reservoir |
| CN104850682A (en) * | 2015-04-17 | 2015-08-19 | 长江大学 | Multiple-point geostatistics modeling method based on position |
| CN107238862A (en) * | 2016-03-29 | 2017-10-10 | 中国石油化工股份有限公司 | Reflectance factor method of estimation and device based on Bayes's inverting framework |
| CN108645994A (en) * | 2018-04-25 | 2018-10-12 | 中国石油大学(北京) | A kind of geology stochastic inversion methods and device based on Multiple-Point Geostatistics |
| CN108931811A (en) * | 2018-05-17 | 2018-12-04 | 长江大学 | Seismic Reservoir inversion method based on multiple spot geological statistics |
| CN109116421A (en) * | 2018-10-23 | 2019-01-01 | 中海石油(中国)有限公司 | A kind of braid deltas reservoir statistics inverted parameters determine method |
| CN110031895A (en) * | 2019-03-11 | 2019-07-19 | 西安科技大学 | A kind of Multiple-Point Geostatistics stochastic inversion methods and device based on image stitching |
| CN110031896A (en) * | 2019-04-08 | 2019-07-19 | 中国石油天然气集团有限公司 | Earthquake stochastic inversion methods and device based on Multiple-Point Geostatistics prior information |
Non-Patent Citations (2)
| Title |
|---|
| 尹艳树 等: "一种基于多点地质统计的储层反演新方法", 《长江大学学报（自科版）》 * |
| 王文鑫 等: "多点地质统计学中训练图像优选方法及其在地质建模中的应用", 《石油勘探与开发》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112014882A (en) * | 2020-09-01 | 2020-12-01 | 电子科技大学 | Pre-stack seismic reflection pattern analysis method based on tensor discrimination dictionary |
| CN112014882B (en) * | 2020-09-01 | 2022-02-15 | 电子科技大学 | Pre-stack seismic reflection pattern analysis method based on tensor discrimination dictionary |
| US11686870B2 (en) | 2021-09-29 | 2023-06-27 | China University Of Petroleum (East China) | Interpretive-guided velocity modeling seismic imaging method and system, medium and device |
| CN115963548A (en) * | 2023-01-16 | 2023-04-14 | 中国矿业大学 | Mine microseismic P-wave arrival time picking model construction method based on inverse deductive learning |
| CN115963548B (en) * | 2023-01-16 | 2024-01-23 | 中国矿业大学 | Mine microseismic P wave arrival time pickup model construction method based on deduction learning |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111505713B (en) | 2021-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10976459B2 (en) | Final statics calculation for automated near surface analysis | |
| CN111273348B (en) | Multipoint geostatistical prestack inversion method based on updated probability ratio constant theory | |
| US11346970B2 (en) | Automatic quality control of seismic travel time | |
| US11016214B2 (en) | Dolomite reservoir prediction method and system based on well and seismic combination, and storage medium | |
| US10067253B2 (en) | Method for determining sedimentary facies using 3D seismic data | |
| CN110031896B (en) | Seismic random inversion method and device based on multi-point geostatistics prior information | |
| US10386511B2 (en) | Seismic survey design using full wavefield inversion | |
| EP0891562B1 (en) | 3-d geologic modelling | |
| US9279905B2 (en) | Method for extracting a thumbnail image from a training image so as to constrain the multipoint geostatistical modeling of the subsoil | |
| CN111505713B (en) | Pre-stack seismic inversion method based on multi-point geological statistics | |
| CN111596364B (en) | Seismic sediment microphase combination analysis method based on high-precision sequence stratum grillwork | |
| CN102176052B (en) | Hierarchical sequence analysis method oriented to generation of three-dimensional hierarchical grids | |
| WO2012139082A1 (en) | Event selection in the image domain | |
| CN108931811B (en) | Seismic Reservoir inversion method based on multiple spot geological statistics | |
| CN108645994A (en) | A kind of geology stochastic inversion methods and device based on Multiple-Point Geostatistics | |
| US6324478B1 (en) | Second-and higher-order traveltimes for seismic imaging | |
| CN112505754B (en) | Method for collaborative partitioning sedimentary microfacies by well-seismic based on high-precision sequence grid model | |
| US11397273B2 (en) | Full waveform inversion in the midpoint-offset domain | |
| CN110954956B (en) | Method for evaluating acquisition trace of observation system and computer-readable storage medium | |
| CN112180440A (en) | AVO characteristic analysis-based prestack stochastic inversion method and system | |
| CN116520407A (en) | Method for predicting thickness of thin sand body based on twice random forest method by utilizing earthquake multi-attribute | |
| WO2022156900A1 (en) | A method of and apparatus for determining a multiple well seismic-to-well tie | |
| CN112147701A (en) | Seismic waveform driven high-resolution seismic inversion method | |
| Schaub et al. | Acoustic impedance estimation after prestack depth migration | |
| Maoshan et al. | Generating high-precision seismic stratigraphic cubes by dip-spreading |
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 |