Pre-stack seismic inversion method and device, electron and storage medium
By using the pre-stack seismic inversion method, and combining lithofacies training maps and probability functions of elastic parameter distribution with forward seismic records and pre-stack time migration angle gathers, the problem of insufficient accuracy in reservoir morphology characterization of the post-stack inversion method is solved, achieving higher inversion accuracy and continuity.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA PETROLEUM & CHEMICAL CORP
- Filing Date
- 2024-11-12
- Publication Date
- 2026-05-12
AI Technical Summary
Existing geostatistical post-stack inversion methods are insufficient in reservoir morphology characterization and accuracy, making it difficult to meet the needs of oil and gas field exploration and development.
The pre-stack seismic inversion method is adopted. By acquiring the pre-stack time migration angle gather of the reservoir, and based on the probability function of the lithofacies training map and elastic parameter distribution, the lithofacies and elastic parameters of the grid points are determined by random path. Combined with forward seismic records and pre-stack time migration angle gathers, iterative calculation is performed to improve the accuracy of the inversion results.
It improves the accuracy and continuity of inversion results, enabling better characterization of reservoir morphology and meeting the needs of oil and gas field exploration and development.
Smart Images
Figure CN122018001A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of oil and gas geophysics, and in particular to a pre-stack seismic inversion method, apparatus, electronic device and storage medium. Background Technology
[0002] Oil and gas field exploration and development primarily focus on reservoirs, and reservoir modeling and stochastic inversion methods have always been hot topics in industry and academia. Geostatistical inversion is an inversion method that combines geological modeling and seismic inversion based on well logging, seismic, and geological data. Existing geostatistical inversion methods are mostly based on traditional geostatistical post-stack inversions, and the reservoir morphology characterization and accuracy of the inversion results need improvement. Summary of the Invention
[0003] To address the aforementioned issues, this application provides a pre-stack seismic inversion method, apparatus, electronic device, and storage medium, which can improve the accuracy of the inversion results obtained.
[0004] This application provides a pre-stack seismic inversion method, including:
[0005] Obtain pre-stack time-off angle gathers of reservoirs in the target region;
[0006] Each grid point in the reservoir is determined sequentially based on a random path, and the lithofacies of the grid points is determined based on the lithofacies training map of the target area.
[0007] The target probability function is determined based on the correspondence between the described lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions;
[0008] The elastic parameters of each time sampling point of the grid points are obtained by random sampling based on the target probability function.
[0009] The forward-modeled seismic record for each time sampling point is determined based on the elastic parameters of each time sampling point.
[0010] The elastic parameters of each time sampling point of the grid are determined based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
[0011] In some embodiments, determining the elastic parameters of each time sampling point of the grid points based on the forward-modeled seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point includes:
[0012] The error value of the objective function is determined based on the forward model seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point, wherein the objective function is the least squares sum of the forward model seismic record and the pre-stack time migration angle gather;
[0013] If the error value is less than the error threshold, the elastic parameter for each time sampling point is output.
[0014] In some embodiments, the method further includes:
[0015] If the error value is greater than the error threshold, the lithofacies of the grid points are updated, and iterative calculations are performed to obtain the updated forward seismic record for each time sampling point;
[0016] The error value of the objective function is determined based on the updated forward model seismic records for each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point;
[0017] If the error value is less than the error threshold or the number of iterations reaches the iteration count threshold, the elastic parameters corresponding to the forward model seismic records of each updated time sampling point are determined as the elastic parameters of each time sampling point of the grid points.
[0018] In some embodiments, determining the lithofacies of the grid points based on the lithofacies training map of the target region includes:
[0019] The similarity between data events in the lithofacies training map and data events at the grid points was calculated using a direct sampling multi-point geostatistical simulation algorithm.
[0020] The lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold are determined as the lithofacies of the grid points.
[0021] In some embodiments, determining the forward-modeled seismic record for each time sampling point based on the elastic parameters of each time sampling point includes:
[0022] The reflection coefficient for each time sampling point at different incident angles is determined based on the elastic parameters at each time sampling point.
[0023] The forward seismic record for each time sampling point is obtained by convolving the reflection coefficients at different incident angles corresponding to each time sampling point with the wavelet.
[0024] In some embodiments, determining the reflection coefficient corresponding to each time sampling point based on the elastic parameter at each time point includes:
[0025] The elastic parameters at each time sampling point are input into the calculation model to determine the reflection coefficient corresponding to each time sampling point. The calculation model includes:
[0026] R(θ) = A + B sin 2 θ;
[0027] in, vs, Δvp, Δvs, and Δρ are the average values of the P-wave velocity, S-wave velocity, and density on both sides of the interface, respectively, and the differences between the P-wave velocity, S-wave velocity, and density on both sides of the interface.
[0028] In some embodiments, the method further includes:
[0029] Acquire well logging data for the target area;
[0030] Based on the well logging data, a correspondence is established between probability functions of different lithofacies and elastic parameter distributions.
[0031] This application provides a pre-stack seismic inversion device, comprising:
[0032] The acquisition module is used to acquire pre-stack time offset angle gathers of reservoirs in the target region;
[0033] The first determining module is used to sequentially determine each grid point in the reservoir based on a random path, and to determine the lithofacies of the grid points based on the lithofacies training map of the target area.
[0034] The second determining module is used to determine the target probability function based on the correspondence between the lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions;
[0035] The module is used to obtain the elastic parameters of each time sampling point of the grid points by random sampling based on the target probability function;
[0036] The third determination module is used to determine the forward modeling seismic record for each time sampling point based on the elastic parameters of each time sampling point.
[0037] The fourth determination module is used to determine the elastic parameters of each time sampling point of the grid points based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
[0038] This application provides an electronic device, including a memory and a processor. The memory stores a computer program, which, when executed by the processor, performs any of the pre-stack seismic inversion methods described above.
[0039] This application provides a storage medium storing a computer program that can be executed by one or more processors and can be used to implement the pre-stack seismic inversion method described in any of the above claims.
[0040] This application provides a pre-stack seismic inversion method, apparatus, electronic device, and storage medium. The method involves: acquiring pre-stack time-migrating angle gathers of reservoirs in a target region; sequentially determining grid points in the reservoir based on random paths; determining the lithofacies of the grid points based on lithofacies training maps of the target region; determining a target probability function based on the correspondence between the lithofacies and pre-established probability functions of different lithofacies and elastic parameter distributions; randomly sampling to obtain the elastic parameters of each time sampling point of the grid points based on the target probability function; determining the forward seismic record of each time sampling point based on the elastic parameters of each time sampling point; and determining the elastic parameters of each time sampling point of the grid points based on the forward seismic record of each time sampling point and the pre-stack time-migrating angle gathers corresponding to each time sampling point. This method improves the accuracy of the inversion results. Attached Figure Description
[0041] The present application will be described in more detail below based on embodiments and with reference to the accompanying drawings.
[0042] Figure 1 A schematic diagram illustrating the implementation process of a pre-stack seismic inversion method provided in this application embodiment;
[0043] Figure 2 A schematic diagram illustrating the implementation process of another pre-stack seismic inversion method provided in this application embodiment;
[0044] Figure 3 This is a schematic diagram of the composition structure of the electronic device provided in the embodiments of this application.
[0045] In the accompanying drawings, the same parts are referred to by the same reference numerals, and the drawings are not drawn to scale. Detailed Implementation
[0046] To make the objectives, technical solutions, and advantages of this application clearer, the application will be further described in detail below with reference to the accompanying drawings. The described embodiments should not be regarded as limitations on this application. All other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.
[0047] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.
[0048] If the application documents contain similar descriptions such as "first, second, third", the following explanation shall be added: In the following description, the terms "first, second, third" are used only to distinguish similar objects and do not represent a specific order of objects. It is understood that "first, second, third" may be interchanged in a specific order or sequence where permitted, so that the embodiments of this application described herein can be implemented in an order other than that illustrated or described herein.
[0049] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used herein is for the purpose of describing embodiments of this application only and is not intended to limit this application.
[0050] To address the problems existing in related technologies, this application provides a pre-stack seismic inversion method. This method is applied to an electronic device, such as a computer or mobile terminal. The functions implemented by the pre-stack seismic inversion method provided in this application can be achieved by the processor of the electronic device calling program code, which can be stored in a computer storage medium.
[0051] Example 1
[0052] This application provides a pre-stack seismic inversion method. Figure 1 This is a schematic diagram illustrating the implementation process of a pre-stack seismic inversion method provided in an embodiment of this application, as shown below. Figure 1 As shown, it includes:
[0053] Step S101: Obtain the pre-stack time offset angle gather of the reservoir in the target region.
[0054] In this embodiment of the application, the pre-stack time migration angle gather of the reservoir in the target area can be obtained by inputting through an input device, which can be a seismic trace measurement device, a storage device, etc.
[0055] In some embodiments, pre-stack time-off angle gathers stored in the target region can be obtained from the Internet.
[0056] Step S102: Determine each grid point in the reservoir sequentially based on a random path, and determine the lithofacies of the grid points based on the lithofacies training map of the target area.
[0057] In this embodiment of the application, a random access path can be set, and each grid point of the reservoir can be determined sequentially based on the random access path.
[0058] In this embodiment of the application, a grid template can be set to divide the reservoir into grids. The grid template can be a 3×3 or 5×5 square, cross, etc.
[0059] In this embodiment of the application, a lithofacies training map of the target area can be pre-established.
[0060] In this embodiment of the application, a lithofacies training map of the target area can be established by combining basic geological research, well logging lithofacies analysis, and seismic attribute analysis with modern sedimentary models.
[0061] In this embodiment of the application, determining the lithofacies of the grid points based on the lithofacies training map of the target area can be achieved through the following steps:
[0062] The similarity between data events in the lithofacies training map and data events at the grid points was calculated using a direct sampling multi-point geostatistical simulation algorithm.
[0063] The lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold are determined as the lithofacies of the grid points.
[0064] In this embodiment, the similarity can be represented by Euclidean distance. In some embodiments, the similarity between data events in the lithofacies training map and data events at the grid points can be calculated using the cosine formula.
[0065] In this embodiment of the application, a lithology classification probability database based on the training image can be established by statistically analyzing the probability that the center point of the grid is a certain lithology, given that the lithology type of the surrounding grid points is determined, based on a set grid template (which can be a 3×3 or 5×5 square and cross shape).
[0066] p(C i |C s ) = p i,s ;
[0067] The similarity is based on the data events formed by the condition points near the point to be estimated and the data patterns in the training image. When it is found that the data events in the training image are exactly the same as the condition data events near the point to be estimated, or the similarity distance between them is less than a certain threshold, the data pattern value can be assigned to the area to be estimated centered on the point to be estimated. That is, the lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold is determined as the lithofacies of the grid point.
[0068] In this embodiment of the application, after determining the lithofacies of a grid point, the lithofacies of the next grid point is determined along a random path, and the above method is repeated until every grid point is determined.
[0069] Step S103: Determine the target probability function based on the correspondence between the lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions.
[0070] In this embodiment, the correspondence between the probability functions of different rock phases and elastic parameter distributions can be expressed as: p(m i |C j ), where m i (i = 1, 2, 3) represent the longitudinal wave velocity V, respectively. p Shear wave velocity V s and density ρ. C j (j=1,2,3) represent different lithofacies.
[0071] In this embodiment of the application, the probability density function of the elastic parameter distribution of different lithofacies can be statistically analyzed based on well logging data.
[0072] In this embodiment of the application, once the lithofacies is determined, the target probability function can be determined based on the correspondence.
[0073] Step S104: Based on the target probability function, random sampling is performed to obtain the elastic parameters of each time sampling point of the grid points.
[0074] In this embodiment of the application, the elastic parameters may include: longitudinal wave velocity, transverse wave velocity, and density.
[0075] Step S105: Determine the forward modeling seismic record for each time sampling point based on the elastic parameters of each time sampling point.
[0076] In this embodiment of the application, step S105 can be implemented through the following steps:
[0077] Step S1051: Determine the reflection coefficient of each time sampling point at different incident angles based on the elastic parameters of each time sampling point.
[0078] In this embodiment of the application, determining the reflection coefficient corresponding to each time sampling point includes:
[0079] The elastic parameters at each time sampling point are input into the calculation model to determine the reflection coefficient corresponding to each time sampling point. The calculation model includes:
[0080] R(θ) = A + B sin 2 θ;
[0081] in, vs, Δvp, Δvs, and Δρ are the average values of the P-wave velocity, S-wave velocity, and density on both sides of the interface, respectively, and the differences between the P-wave velocity, S-wave velocity, and density on both sides of the interface.
[0082] That is, the average longitudinal wave velocity on both sides of the interface is Δvp, the average transverse wave velocity on both sides of the interface is Δvs, and the average density on both sides of the interface is... The difference between the longitudinal wave velocities on both sides of the interface is Δvp, the difference between the transverse wave velocities on both sides of the interface is Δvs, and the difference between the densities on both sides of the interface is Δρ.
[0083] In this embodiment of the application, when the incident angle is less than 30°, the approximate formula for the longitudinal wave reflection coefficient is the calculation model described above.
[0084] Step S1052: The reflection coefficients of different incident angles corresponding to each time sampling point are convolved with the wavelet to obtain the forward seismic record of each time sampling point.
[0085] In this embodiment of the application, each time sampling point can be traversed to obtain the forward modeling seismic record of each time sampling point.
[0086] Step S106: Determine the elastic parameters of each time sampling point of the grid points based on the forward modeling seismic records of each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point.
[0087] In this embodiment of the application, step S106 can be implemented through the following steps:
[0088] Step S1061: Determine the error value of the objective function based on the forward model seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point, wherein the objective function is the least squares sum of the forward model seismic record and the pre-stack time migration angle gather.
[0089] Step S1062: If the error value is less than the error threshold, output the elastic parameter for each time sampling point.
[0090] In this embodiment of the application, the elastic parameters output for each time sampling point may include: the current lithofacies C i and longitudinal wave velocity V p Shear wave velocity V s and density ρ.
[0091] Step S1063: If the error value is greater than the error threshold, update the lithofacies of the grid points and perform iterative calculations to obtain the updated forward seismic record for each time sampling point;
[0092] Step S1064: Determine the error value of the objective function based on the updated forward model seismic records for each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point;
[0093] Step S1065: If the error value is less than the error threshold or the number of iterations reaches the iteration count threshold, the elastic parameters corresponding to the forward model seismic records of each updated time sampling point are determined as the elastic parameters of each time sampling point of the grid points.
[0094] In this embodiment of the application, each grid point can be traversed to output the current lithofacies C of all grid points. i and longitudinal wave velocity V p Shear wave velocity V s and density ρ.
[0095] This application provides a pre-stack seismic inversion method, which involves: acquiring pre-stack time-migrating angle gathers of reservoirs in a target region; sequentially determining grid points in the reservoir based on random paths; determining the lithofacies of the grid points based on lithofacies training maps of the target region; determining a target probability function based on the correspondence between the lithofacies and pre-established probability functions of different lithofacies and elastic parameter distributions; randomly sampling to obtain the elastic parameters of each time sampling point of the grid points based on the target probability function; determining the forward seismic record of each time sampling point based on the elastic parameters of each time sampling point; and determining the elastic parameters of each time sampling point of the grid points based on the forward seismic record of each time sampling point and the pre-stack time-migrating angle gathers corresponding to each time sampling point. This method can improve the accuracy of the inversion results and ensure the continuity of the inversion results.
[0096] Example 2
[0097] Based on the foregoing embodiments, this application provides another pre-stack seismic inversion method, which directly combines a multi-point geostatistical algorithm with a Markov chain Monte Carlo multi-angle pre-stack simultaneous inversion method to improve the continuity and accuracy of the inversion results.
[0098] First, we will introduce the multi-point geostatistical simulation method:
[0099] Since Guardiano & Srivastava first proposed the concepts of training images and multi-point geostatistics in 1993, multi-point geostatistics has not been widely applied. In training images, based on a set grid template (which can be a 3×3 or 5×5 square or cross shape), the probability of the grid center point being a certain lithology is statistically calculated, given that the lithology types of surrounding grid points are determined, thus establishing a lithology classification probability database based on the training images.
[0100] p(C i |C s ) = p i,s ;
[0101] Both the Direct Sampling algorithm and traditional multi-point geostatistical methods employ a pixel-based sequential approach. Compared to traditional SNESIM and SIMPAT multi-point modeling, the Direct Sampling algorithm differs in its assignment of values to the points to be estimated. It bases its approach on the similarity between the data events formed by conditional points near the point to be estimated and the data patterns in the training image. When it is found that the data events in the training image are exactly the same as the conditional data events near the point to be estimated, or that the similarity distance between them is less than a certain threshold, the data pattern value can be assigned to the area to be estimated centered on the point to be estimated. Then, the next point to be estimated is simulated according to a specified path until the entire grid of the work area has been traversed. The specific steps are as follows:
[0102] (1) Establish training images that conform to the geological understanding of the work area;
[0103] (2) Determine a random access path for accessing the grid points to be estimated, and then randomly access the grid points of the simulated work area according to the random path.
[0104] (3) Determine the conditional data events surrounding the grid point to be estimated;
[0105] (4) Calculate the similarity between the data events of the grid point to be estimated and the data events of the training image, and assign the data event with the highest similarity in the training image to the grid point to be estimated;
[0106] (5) Visit the next grid node along a random path, repeating steps (3) and (4) until all grid nodes have been visited.
[0107] Calculation of the reflection coefficient of longitudinal waves at any incident angle:
[0108] When the incident angle is less than 30°, the approximate formula for the longitudinal wave reflection coefficient is:
[0109] R(θ) = A + B sin 2 θ;
[0110] in, and Δvp,vs, Δρ represents the average and difference of the longitudinal wave velocity, transverse wave velocity, and density on both sides of the interface.
[0111] Figure 2 This is a flowchart illustrating a pre-stack seismic inversion method provided in an embodiment of this application, as shown below. Figure 2 As shown, the steps are as follows:
[0112] 1. Input the pre-stack time offset angle gather S o .
[0113] 2. Based on well logging data, statistically analyze the probability density function p(m) of the elastic parameter distribution for different lithofacies.i |C j m i (i = 1, 2, 3) represent the longitudinal wave velocity V, respectively. p Shear wave velocity V s and density ρ. C j (j=1,2,3) represent different lithofacies.
[0114] 3. Establish training images of the lithofacies that are consistent with the work area.
[0115] By combining basic geological research, well logging lithofacies analysis, and seismic attribute analysis with modern sedimentary models, lithofacies training images that conform to the work area are established.
[0116] 4. Randomly access CDP grid points sequentially using a random path. Based on the lithofacies distribution characteristics around each time sampling point of the CDP grid point, match the training image library. Use a direct sampling multi-point geostatistical simulation algorithm to calculate the distance between data events in the training images and data events in the CDP grid points. If the distance is less than a given threshold, directly assign this point from the training image to the point to be estimated, obtaining the lithofacies C of the point to be estimated. i .
[0117] 5. In lithofacies C i Under constraints, based on the probability density function p(m) of the elastic parameter distribution of different lithofacies according to well logging statistics. i |C j Random sampling simulates the P-wave velocity V at the point to be estimated. p Shear wave velocity V s and density ρ.
[0118] 6. Repeat step 5 until every time sampling point of the current CDP grid point has been traversed.
[0119] 7. Substitute the elastic parameters of each time sampling point into the approximate formula for the P-wave reflection coefficient (same as the calculation model in the above embodiment), perform forward modeling of the reflection coefficient R at different incident angles, and convolve it with the wavelet to obtain the forward modeled seismic record S. M .
[0120] 8. Calculate the forward modeled seismic record S M Compared with actual earthquake records S o The objective function is the least squares sum of errors. When the error is less than a set threshold, or the number of iterations reaches a set number, the current lithofacies C is output. i and longitudinal wave velocity V p Shear wave velocity V s And density ρ; if it cannot exceed the set threshold value, or the number of iterations is less than the set number, resample and update the lithofacies C. i Proceed to step 5.
[0121] 9. Continue until all time sampling points of all CDP grid points in space are traversed, and the final lithofacies C is output. i and longitudinal wave velocity V p Shear wave velocity V s And density ρ inversion results.
[0122] Based on the foregoing embodiments, this application provides a specific application example, taking a channel reservoir model as an example. The channel reservoir model includes structural and stratigraphic features. The structural features of the channel reservoir are not obvious; it is a lithological trap. The top depth is 1000m, the surface depth is 1045m, and the channel reservoir is a two-dimensional regular stratigraphic model with 200×1×45 elements. The grid size is 25m in the horizontal direction (x and y) and 1m in the vertical direction (z). The reservoir scale is 5000m east-west and 45m thick, constituting a single-layer river system. In the channel model, the channel sand bodies are basically isolated, with no obvious contact between channels. The width of the channel sand bodies is 150–300m, and the thickness ranges from 5–15m. The entire reservoir consists of two facies: sandstone facies and mudstone facies.
[0123] Rock physical properties include density, P-wave velocity, and S-wave velocity. It is assumed that the reservoir's rock physical properties follow a multi-Gaussian distribution, the parameter distribution of each phase follows a Gaussian distribution, and there is overlap between the rock physical properties of different phases. The density distribution of sandstone is approximately 2.3–2.6 g / cm³. 3 The longitudinal wave velocity ranges approximately from 4200 to 4600 m / s, and the transverse wave velocity ranges approximately from 2400 to 2700 m / s; the density of the mudstone ranges approximately from 2.15 to 2.35 g / cm³. 3Between these values, the P-wave velocity ranges approximately from 3900 to 4300 m / s, and the S-wave velocity ranges approximately from 2200 to 2500 m / s. Seven virtual wells were extracted from the channel reservoir model as conditional data. In the initial iteration, the reservoir lithology model showed that the channel sand bodies had a flat-top, convex-bottom morphology on the cross-section, and the width-to-thickness ratio of the sand bodies was similar to the theoretical model. However, the individual channel sand bodies were generally small in size, and the sand bodies were randomly distributed with complex contact relationships between them, including isolated single sand bodies and composite sand bodies in lateral, oblique, and vertical contact. The synthetic seismic record makes it difficult to identify the reservoir as a channel reservoir; it resembles more of a long, continuous, thin-layered sand body. After the third iteration, the main distribution locations of the channel sand bodies became clearer, with only scattered sand bodies in the background sediment. The morphology and distribution of the channel could also be seen in the synthetic seismic records from the red high-value and blue low-value seismic data. However, the vertical resolution of the synthetic records remained low. Although the width of the channel could be estimated from the synthetic records, the low resolution of the seismic data itself prevented a relatively accurate estimation of the thickness of the channel sand bodies. For example, within the gray elliptical box, the sand bodies appeared isolated in the reservoir model, while the corresponding seismic data might show a stacked relationship. The third iteration was able to reflect the reservoir distribution quite well. The results of the sixth iteration showed a more pronounced flat-top, convex-bottom morphology of the channel sand bodies, better matching the designed geological model and reflecting the distribution characteristics of the sand bodies more effectively. In the first and second iterations, the seismic traces showed good correlation only in the near-wellbore region, while the correlation was poor in the far-well region (such as near traces 30, 81, and 130), and even negatively correlated trace sets appeared. During the iteration process, as the number of iterations increased, the correlation coefficient of each seismic trace tended to 1, indicating that the correlation between the inversion results and the theoretical model was gradually increasing. After the sixth iteration, the average correlation coefficient of the seismic traces reached 0.92, indicating that the inversion results after 6 iterations had a strong correlation and good inversion effect.
[0124] Example 3
[0125] Based on the foregoing embodiments, this application provides a pre-stack seismic inversion device. The various modules and units included in the device can be implemented by a processor in a computer device; of course, they can also be implemented by specific logic circuits. In the implementation process, the processor can be a central processing unit (CPU), a microprocessor unit (MPU), a digital signal processor (DSP), or a field programmable gate array (FPGA), etc.
[0126] This application provides a pre-stack seismic inversion device, which includes:
[0127] The acquisition module is used to acquire pre-stack time offset angle gathers of reservoirs in the target region;
[0128] The first determining module is used to sequentially determine each grid point in the reservoir based on a random path, and to determine the lithofacies of the grid points based on the lithofacies training map of the target area.
[0129] The second determining module is used to determine the target probability function based on the correspondence between the lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions;
[0130] The module is used to obtain the elastic parameters of each time sampling point of the grid points by random sampling based on the target probability function;
[0131] The third determination module is used to determine the forward modeling seismic record for each time sampling point based on the elastic parameters of each time sampling point.
[0132] The fourth determination module is used to determine the elastic parameters of each time sampling point of the grid points based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
[0133] In some embodiments, determining the elastic parameters of each time sampling point of the grid points based on the forward-modeled seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point includes:
[0134] The error value of the objective function is determined based on the forward model seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point, wherein the objective function is the least squares sum of the forward model seismic record and the pre-stack time migration angle gather;
[0135] If the error value is less than the error threshold, the elastic parameter for each time sampling point is output.
[0136] In some embodiments, the method further includes:
[0137] If the error value is greater than the error threshold, the lithofacies of the grid points are updated, and iterative calculations are performed to obtain the updated forward seismic record for each time sampling point;
[0138] The error value of the objective function is determined based on the updated forward model seismic records for each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point;
[0139] If the error value is less than the error threshold or the number of iterations reaches the iteration count threshold, the elastic parameters corresponding to the forward model seismic records of each updated time sampling point are determined as the elastic parameters of each time sampling point of the grid points.
[0140] In some embodiments, determining the lithofacies of the grid points based on the lithofacies training map of the target region includes:
[0141] The similarity between data events in the lithofacies training map and data events at the grid points was calculated using a direct sampling multi-point geostatistical simulation algorithm.
[0142] The lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold are determined as the lithofacies of the grid points.
[0143] In some embodiments, determining the forward-modeled seismic record for each time sampling point based on the elastic parameters of each time sampling point includes:
[0144] The reflection coefficient for each time sampling point at different incident angles is determined based on the elastic parameters at each time sampling point.
[0145] The forward seismic record for each time sampling point is obtained by convolving the reflection coefficients at different incident angles corresponding to each time sampling point with the wavelet.
[0146] In some embodiments, determining the reflection coefficient corresponding to each time sampling point based on the elastic parameter at each time point includes:
[0147] The elastic parameters at each time sampling point are input into the calculation model to determine the reflection coefficient corresponding to each time sampling point. The calculation model includes:
[0148] R(θ) = A + B sin 2 θ;
[0149] in, vs, Δvp, Δvs, and Δρ are the average values of the P-wave velocity, S-wave velocity, and density on both sides of the interface, respectively, and the differences between the P-wave velocity, S-wave velocity, and density on both sides of the interface.
[0150] In some embodiments, the method further includes:
[0151] Acquire well logging data for the target area;
[0152] Based on the well logging data, a correspondence is established between probability functions of different lithofacies and elastic parameter distributions.
[0153] Example 4
[0154] This application provides an electronic device; Figure 3 This is a schematic diagram of the composition structure of the electronic device provided in the embodiments of this application, such as... Figure 3 As shown, the electronic device 700 includes: a processor 701, at least one communication bus 702, a user interface 703, at least one external communication interface 704, and a memory 705. The communication bus 702 is configured to enable communication between these components. The user interface 703 may include a display screen, and the external communication interface 704 may include standard wired and wireless interfaces. The processor 701 is configured to execute a program of a pre-stack seismic inversion method stored in the memory to implement the steps in the pre-stack seismic inversion method provided in the above embodiment.
[0155] The pre-stack seismic inversion method includes:
[0156] This application provides a pre-stack seismic inversion method, including:
[0157] Obtain pre-stack time-off angle gathers of reservoirs in the target region;
[0158] Each grid point in the reservoir is determined sequentially based on a random path, and the lithofacies of the grid points is determined based on the lithofacies training map of the target area.
[0159] The target probability function is determined based on the correspondence between the described lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions;
[0160] The elastic parameters of each time sampling point of the grid points are obtained by random sampling based on the target probability function.
[0161] The forward-modeled seismic record for each time sampling point is determined based on the elastic parameters of each time sampling point.
[0162] The elastic parameters of each time sampling point of the grid are determined based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
[0163] In some embodiments, determining the elastic parameters of each time sampling point of the grid points based on the forward-modeled seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point includes:
[0164] The error value of the objective function is determined based on the forward model seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point, wherein the objective function is the least squares sum of the forward model seismic record and the pre-stack time migration angle gather;
[0165] If the error value is less than the error threshold, the elastic parameter for each time sampling point is output.
[0166] In some embodiments, the method further includes:
[0167] If the error value is greater than the error threshold, the lithofacies of the grid points are updated, and iterative calculations are performed to obtain the updated forward seismic record for each time sampling point;
[0168] The error value of the objective function is determined based on the updated forward model seismic records for each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point;
[0169] If the error value is less than the error threshold or the number of iterations reaches the iteration count threshold, the elastic parameters corresponding to the forward model seismic records of each updated time sampling point are determined as the elastic parameters of each time sampling point of the grid points.
[0170] In some embodiments, determining the lithofacies of the grid points based on the lithofacies training map of the target region includes:
[0171] The similarity between data events in the lithofacies training map and data events at the grid points was calculated using a direct sampling multi-point geostatistical simulation algorithm.
[0172] The lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold are determined as the lithofacies of the grid points.
[0173] In some embodiments, determining the forward-modeled seismic record for each time sampling point based on the elastic parameters of each time sampling point includes:
[0174] The reflection coefficient for each time sampling point at different incident angles is determined based on the elastic parameters at each time sampling point.
[0175] The forward seismic record for each time sampling point is obtained by convolving the reflection coefficients at different incident angles corresponding to each time sampling point with the wavelet.
[0176] In some embodiments, determining the reflection coefficient corresponding to each time sampling point based on the elastic parameter at each time point includes:
[0177] The elastic parameters at each time sampling point are input into the calculation model to determine the reflection coefficient corresponding to each time sampling point. The calculation model includes:
[0178] R(θ) = A + Bsin 2 θ;
[0179] in, vs, Δvp, Δvs, and Δρ are the average values of the P-wave velocity, S-wave velocity, and density on both sides of the interface, respectively, and the differences between the P-wave velocity, S-wave velocity, and density on both sides of the interface.
[0180] In some embodiments, the method further includes:
[0181] Acquire well logging data for the target area;
[0182] Based on the well logging data, a correspondence is established between probability functions of different lithofacies and elastic parameter distributions.
[0183] Example 5
[0184] In this application embodiment, if the above-described pre-stack seismic inversion method is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application embodiment, in essence, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), magnetic disks, or optical disks. Thus, this application embodiment is not limited to any specific hardware and software combination.
[0185] Accordingly, this application provides a storage medium storing a computer program thereon, characterized in that the computer program, when executed by a processor, implements the steps in the pre-stack seismic inversion method provided in the above embodiments.
[0186] The descriptions of the above embodiments of the electronic devices and storage media are similar to those of the above method embodiments, and have similar beneficial effects. For technical details not disclosed in the embodiments of the computer devices and storage media of this application, please refer to the descriptions of the method embodiments of this application for understanding.
[0187] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. It should be understood that in the various embodiments of this application, the sequence numbers of the above-described processes do not imply a sequential order of execution; the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application. The sequence numbers of the above-described embodiments are merely descriptive and do not represent the superiority or inferiority of the embodiments.
[0188] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.
[0189] In the several embodiments provided in this application, it should be understood that the disclosed devices and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods, such as: multiple units or components can be combined, or integrated into another system, or some features can be ignored or not executed. In addition, the coupling, direct coupling, or communication connection between the various components shown or discussed can be through some interfaces, and the indirect coupling or communication connection between devices or units can be electrical, mechanical, or other forms.
[0190] The units described above as separate components may or may not be physically separate. The components shown as units may or may not be physical units. They may be located in one place or distributed across multiple network units. Some or all of the units may be selected to achieve the purpose of this embodiment according to actual needs.
[0191] In addition, each functional unit in the various embodiments of this application can be integrated into one processing unit, or each unit can be a separate unit, or two or more units can be integrated into one unit; the integrated unit can be implemented in hardware or in the form of hardware plus software functional units.
[0192] Those skilled in the art will understand that all or part of the steps of the above method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it performs the steps of the above method embodiments. The aforementioned storage medium includes various media that can store program code, such as mobile storage devices, read-only memory (ROM), magnetic disks, or optical disks.
[0193] Alternatively, if the integrated units described above are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this application, or the parts that contribute to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a controller to execute all or part of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, ROMs, magnetic disks, or optical disks.
[0194] The above description is merely an embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.
Claims
1. A pre-stack seismic inversion method, characterized in that, include: Obtain pre-stack time-off angle gathers of reservoirs in the target region; Each grid point in the reservoir is determined sequentially based on a random path, and the lithofacies of the grid points is determined based on the lithofacies training map of the target area. The target probability function is determined based on the correspondence between the described lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions; The elastic parameters of each time sampling point of the grid points are obtained by random sampling based on the target probability function. The forward-modeled seismic record for each time sampling point is determined based on the elastic parameters of each time sampling point. The elastic parameters of each time sampling point of the grid are determined based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
2. The method according to claim 1, characterized in that, The determination of the elastic parameters of each time sampling point of the grid points based on the forward modeled seismic records of each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point includes: The error value of the objective function is determined based on the forward model seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point, wherein the objective function is the least squares sum of the forward model seismic record and the pre-stack time migration angle gather; If the error value is less than the error threshold, the elastic parameter for each time sampling point is output.
3. The method according to claim 2, characterized in that, The method further includes: If the error value is greater than the error threshold, the lithofacies of the grid points are updated, and iterative calculations are performed to obtain the updated forward seismic record for each time sampling point; The error value of the objective function is determined based on the updated forward model seismic records for each time sampling point and the pre-stack time migration angle gathers corresponding to each time sampling point; If the error value is less than the error threshold or the number of iterations reaches the iteration count threshold, the elastic parameters corresponding to the forward model seismic records of each updated time sampling point are determined as the elastic parameters of each time sampling point of the grid points.
4. The method according to claim 1, characterized in that, The determination of the lithofacies of the grid points based on the lithofacies training map of the target region includes: The similarity between data events in the lithofacies training map and data events at the grid points was calculated using a direct sampling multi-point geostatistical simulation algorithm. The lithofacies corresponding to the lithofacies training image with a similarity greater than the similarity threshold are determined as the lithofacies of the grid points.
5. The method according to claim 1, characterized in that, The determination of the forward-modeled seismic record for each time sampling point based on the elastic parameters of each time sampling point includes: The reflection coefficient for each time sampling point at different incident angles is determined based on the elastic parameters at each time sampling point. The forward seismic record for each time sampling point is obtained by convolving the reflection coefficients at different incident angles corresponding to each time sampling point with the wavelet.
6. The method according to claim 5, characterized in that, The determination of the reflection coefficient corresponding to each time sampling point based on the elastic parameters at each time point includes: The elastic parameters at each time sampling point are input into the calculation model to determine the reflection coefficient corresponding to each time sampling point. The calculation model includes: R(θ)=A+B sin 2 I; in, vs, Δvp, Δvs, and Δρ are the average values of the P-wave velocity, S-wave velocity, and density on both sides of the interface, respectively, and the differences between the P-wave velocity, S-wave velocity, and density on both sides of the interface.
7. The method according to claim 1, characterized in that, The method further includes: Acquire well logging data for the target area; Based on the well logging data, a correspondence is established between probability functions of different lithofacies and elastic parameter distributions.
8. A pre-stack seismic inversion device, characterized in that, include: The acquisition module is used to acquire pre-stack time offset angle gathers of reservoirs in the target region; The first determining module is used to sequentially determine each grid point in the reservoir based on a random path, and to determine the lithofacies of the grid points based on the lithofacies training map of the target area. The second determining module is used to determine the target probability function based on the correspondence between the lithofacies and the pre-established probability functions of different lithofacies and elastic parameter distributions; The module is used to obtain the elastic parameters of each time sampling point of the grid points by random sampling based on the target probability function; The third determination module is used to determine the forward modeling seismic record for each time sampling point based on the elastic parameters of each time sampling point. The fourth determination module is used to determine the elastic parameters of each time sampling point of the grid points based on the forward modeling seismic record of each time sampling point and the pre-stack time migration angle gather corresponding to each time sampling point.
9. An electronic device, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program that, when executed by the processor, performs the pre-stack seismic inversion method as described in any one of claims 1 to 7.
10. A storage medium, characterized in that, The computer program stored in the storage medium can be executed by one or more processors and can be used to implement the pre-stack seismic inversion method as described in any one of claims 1 to 7.