Anisotropic pre-stack inversion method and device, equipment and storage medium

Abstract

The embodiment of the invention provides a thin reservoir geostatistical anisotropy pre-stack inversion method, device and equipment and a storage medium. The method comprises the following steps: obtaining a positive operator in inversion by using an anisotropy fast reflectivity method; a multiple order in the positive operator is controlled; a full-band interpolation model is established in combination with logging data, and an inverted low-frequency initial model is obtained through smoothing processing; an inversion objective function based on a geostatistical method is constructed, the resolution of an inversion result is improved, and full wave field information is fully utilized to predict reservoir parameters; and determining an optimal elastic parameter model according to the inversion objective function and the PP wave inversion residual error. According to the method, the resolution of the inversion result is higher, theoretical support can be provided for prediction of the tight sandstone thin reservoir, and the tight sandstone reservoir can be better identified.

Description

Technical Field

[0001] This application relates to the field of reservoir geophysics, specifically to a method, apparatus, equipment, and storage medium for pre-stack inversion of anisotropic thin reservoir geostatistics. Background Technology

[0002] Tight sandstone thin reservoirs with horizontal fractures can be equivalent to thin interbedded media with transversely isotropic (VTI) media. Their internal structures are diverse, with numerous combinations and complex time-frequency response characteristics, leading to complex converted waves. This renders formulas based on the single-interface assumption, such as the Zoeppritz equation and the Rüger formula, inapplicable for exploring the intrinsic response mechanisms and time-frequency variations of this type of reservoir.

[0003] To overcome the inapplicability of the theoretical formula under the single-interface assumption, Booth and Crampin (1983), based on Kennett's (1974) recursive algorithm, improved the method of Fuchs and Muller (1971) under the assumption of weak anisotropy, extending the reflectivity method to anisotropic media. Fryer and Frazer (1984) modified Kennett's (1974) recursive algorithm, extending the reflectivity method to anisotropic layered media. Subsequently, Fryer and Frazer (1987) gave an analytical expression for the elastic dynamics of media with horizontal symmetry planes, making the formula applicable to layered media with multiple sets of intersecting vertical cracks, but not to media with inclined cracks.

[0004] To address this problem, Frazer and Fryer (1989) developed two numerical methods for the general triclinic case. Mallick and Frazer (1988) rearranged the source and receiver response formulas in Kennett's (1983) reflectivity method, extending the reflectivity method to VSP and implementing a vector algorithm to accelerate computation. This formula is applicable to both isotropic and anisotropic media. Mallick and Frazer (1990) combined the numerical techniques of Frazer and Fryer (1989) with Mallick and Frazer's (1988) efficient Kennett reflectivity method, providing a reflectivity method applicable to arbitrary anisotropic layered media. Carcione (2001, 2007) combined the recursive method in Kennett's (1983) reflectivity method to provide a wavefield simulation method describing horizontally multilayered viscoelastic anisotropic media. Deng et al. (2018) used the anisotropic reflectance method (Fryer and Frazer, 1987) to perform full-wavefield seismic forward modeling of thin-layered shale reservoirs and calculated the AVA response of anisotropic shale oil reservoirs.

[0005] In terms of inversion, several scholars have introduced anisotropic propagation matrix theory into pre-stack inversion methods for VTI media (Kamath and Tsvankin, 2013; Padhi and Mallick, 2014). Chang and McMechan (2009) conducted a feasibility study of FWI for layered media containing VTI media and successfully completed the inversion of model data. Kamath and Tsvankin (2013) established a depth-domain FWI method based on anisotropic reflectivity and applied it to thin interlayered model media containing VTI media. Under the assumption of local one-dimensionality, Padhi and Mallick (2013) established a PWI method for simple models of isotropic and VTI media and successfully performed inversion using multi-component seismic data. Padhi and Mallick (2014) extended the PWI method for VTI media to complex models and applied it to actual data containing thin layers. Li and Mallick (2015) extended VTI inversion to azimuth anisotropy inversion and applied it to real-world data with four components and two azimuths, including thin layers. Although several scholars have published related articles, pre-stack waveform inversion based on anisotropic propagation matrices still requires further research. The main reason is that most published applications to real-world data are limited to a few CMP gathers or within a small time window (Li and Mallick, 2015). The main factors restricting the development of this type of inversion technique are the complexity of the formulas themselves, their higher degree of nonlinearity, and the lack of rapid inversion similar to that in isotropic media (e.g., Yang and Lu, 2020b; Li et al., 2021).

[0006] Due to limitations of deterministic inversion methods, the inversion results based on the aforementioned methods are insufficient for distinguishing thin reservoirs. Geostatistical inversion is a stochastic inversion method that combines geostatistics with seismic inversion. Based on geostatistics, using well logging data as conditional data, and seismic data as constraints, stochastic inversion methods have a significant advantage in resolution compared to conventional deterministic inversion.

[0007] The stochastic inversion method was first proposed by Haas and Dubrule at the First Break conference in 1994. This method significantly improves the vertical resolution of the inversion and captures frequency information beyond the seismic bandwidth. Contreras et al. (2005) successfully applied stochastic inversion to pre-stack seismic data by combining pre-stack seismic data with geostatistical simulation. In 2010, Cordua et al. incorporated Bayesian theory to integrate the prior information described by the geostatistical algorithm, thus obtaining a good posterior probability distribution of the model parameters and realizing AVO stochastic inversion. Gan Lideng et al. (2011) pointed out that stochastic inversion has greater advantages in seismic inversion of thin reservoirs, with less uncertainty compared to deterministic inversion. Huang Handong et al. (2011) combined Bayesian theory with Gaussian prior information to realize the combination of wave impedance inversion algorithm and multichannel inversion, and applied it to thin layer prediction. However, the above methods are all based on the isotropic assumption and have not fully considered the anisotropic characteristics of the reservoir.

[0008] There are relatively few studies on anisotropic inversion within a statistical framework. This is because anisotropic media require consideration of numerous physical parameters, including porosity, saturation, fracture density, fracture morphology, and even mineral composition. Inverting so many unknown physical parameters under seismic data constraints is not only computationally intensive but also results in highly uncertain inversion results. Therefore, it is necessary to use assumptions or experience to reduce the number of unknown parameters. Ran et al. (2013) proposed using statistical rock physics to obtain anisotropic parameters, adding constraints to pre-stack seismic anisotropic parameter inversion to reduce its ambiguity. Zhang et al. (2018) proposed a statistical rock physics inversion method for anisotropic shale reservoirs, obtaining the posterior probability density function of reservoir parameters within a Bayesian framework, and deriving the optimal estimates of the parameters and a quantitative description of their uncertainties.

[0009] Although stochastic inversion can improve the resolution of inversion results and offers higher certainty in the inversion of thin reservoirs (thin interbedded layers) (Sams and Saussus, 2008), no scholars have yet performed statistical stochastic inversion based on anisotropic propagation matrices. The main reason is that geostatistical stochastic inversion itself is computationally intensive, and the complex formulas further reduce its efficiency. Therefore, there is an urgent need to design a method, apparatus, and equipment for thin reservoir geostatistical anisotropic pre-stack inversion that can achieve anisotropic pre-stack geostatistical stochastic inversion and more accurately predict the development of thin interbedded layers in tight sandstone. Summary of the Invention

[0010] In view of this, this application provides a thin reservoir geostatistical anisotropic pre-stack inversion method, apparatus, equipment and storage medium. Based on the reflectance method, thin-layer VTI medium simulation theory and geostatistical inversion theory, a high-precision, high-stability and high-efficiency thin reservoir geostatistical anisotropic pre-stack inversion method is constructed, providing theoretical support for the prediction of tight sandstone thin reservoirs, solving some of the problems in the exploration of tight sandstone thin reservoirs, and providing guidance for finding high-quality tight gas reservoir development areas.

[0011] In a first aspect, embodiments of this application provide a pre-stack inversion method for thin reservoir geostatistical anisotropy, including:

[0012] The reflection coefficients were determined for both the VTI single-layer model and the thin interlayer dielectric model with VTI single-layer sandwiched between them.

[0013] Based on geostatistical methods, the reflection coefficient is used to perform pre-stack anisotropic inversion of tight clastic thin reservoirs.

[0014] In one possible implementation, the formula for calculating the reflection coefficient is:

[0015]

[0016] Among them, V 04 V 01 These are the first and fourth elements of array v0, respectively.

[0017] In one possible implementation, v0 is a 6-element array defined in the slowness-frequency domain (p, ω), and the expression for v0 is:

[0018] v0=[Δ -R PS Δ -R SS Δ R PP Δ R SP Δ detRΔ] T

[0019] Where Δ represents the proportionality coefficient, R PS R SS R PP R SP denoted by and detR respectively, representing the total reflection coefficients for the corresponding wave types.

[0020] In one possible implementation, the pre-stack anisotropic inversion of tight clastic thin reservoirs using the reflection coefficient based on geostatistical methods includes:

[0021] Acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times;

[0022] Define a random path that accesses all grid nodes in the model space;

[0023] Within the neighborhood of a node, well logging data and simulated point data are selected as class A data d. obsA Select angle gather data as class B data d obsB Constructing hybrid observation data d obsA+B ;

[0024] Based on the mixed observation data d obsA+B By combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance

[0025] The random realization value m of the current node is obtained by uniformly sampling the local Gaussian cumulative distribution function. sim (j) and add the random implementation value to the known point data;

[0026] The calculations are performed sequentially on all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result;

[0027] Regenerate random paths that can access all grid nodes in the model space, and perform anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementations of the model parameters m. sim Posterior expectation and posterior covariance

[0028] In one possible implementation, the observation data d are mixed. obsA+B The calculation formula is as follows:

[0029] d obsA+B =G A+B m A+B +n A+B

[0030] Among them, G A+B The forward operator represents a mixture of data of type A and type B, m. A+B n represents the model parameters for mixed data of type A and type B. A+B This represents the error term for mixed data of type A and type B.

[0031] In one possible implementation, the posterior expectation and posterior covariance The calculation formula is as follows:

[0032]

[0033] Where, m priorA+B C represents the initial pre-stack elasticity parameter. mA+B Let G represent the model covariance matrix. A+B The forward operator representing mixed data of type A and type B, d obsA+B Let G represent the mixed observation data, and let G represent the nonlinear operator that maps the model parameter vector m to the observation data vector d.

[0034] In one possible implementation, the model covariance matrix C mA+B and data covariance matrix C dA+B The calculation formula is as follows:

[0035]

[0036] Among them, C mAA and C mBB C represents the model covariance matrix for class A data and class B data, respectively; mAB Represents the cross-covariance matrix of the two types of data; C represents the transpose of the cross-covariance matrix of the two classes of data; dBB The data covariance matrix representing class B data.

[0037] Secondly, embodiments of this application provide a thin reservoir geostatistical anisotropy pre-stack inversion device, comprising:

[0038] The reflection coefficient determination module is used to determine the reflection coefficient for VTI single-layer models and thin interlayer media models with VTI single-layers.

[0039] A pre-stack inversion module is used to perform anisotropic pre-stack inversion of tight clastic thin reservoirs based on geostatistical methods and utilizing the reflection coefficient; wherein, the pre-stack inversion module includes:

[0040] The data acquisition module is used to acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times;

[0041] The path definition module is used to define a random path that accesses all grid nodes in the model space;

[0042] The data construction module is used to select well logging data and simulated point data as class A data within the neighborhood of a node. obsA Select angle gather data as class B data d obsB Constructing hybrid observation data d obsA+B ;

[0043] The node calculation module is used to calculate the mixed observation data d. obsA+BBy combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance

[0044] The random implementation value acquisition module is used to uniformly sample the local Gaussian cumulative distribution function to obtain the random implementation value m of the current node. sim (j) and add the random implementation value to the known point data;

[0045] The result generation module is used to sequentially calculate all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result. The path definition module 502 regenerates a random path that can access all grid nodes in the model space, and performs anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementation results m of the model parameters. sim Posterior expectation and posterior covariance

[0046] Thirdly, embodiments of this application provide an electronic device, including:

[0047] processor;

[0048] Memory;

[0049] And a computer program, wherein the computer program is stored in the memory, the computer program including instructions that, when executed by the processor, cause the electronic device to perform the method described in any one of the first aspects.

[0050] Fourthly, embodiments of this application provide a computer-readable storage medium including a stored program, wherein, when the program is executed, it controls the device where the computer-readable storage medium is located to perform the method described in any one of the first aspects.

[0051] The present invention adopts the above technical solution and has the following advantages:

[0052] (1) Tight sandstone reservoirs with low-angle fractures are approximated as VTI media. Anisotropic parameters are calculated based on pseudo-anisotropic theory. Then, the calculation method for single-interface reflection coefficient and the recursive algorithm in the total reflection coefficient are improved to obtain a reflectivity method formula applicable to anisotropic media. In addition, to speed up the calculation, the previously derived reflectivity method formula is reorganized by analogy with the fast reflectivity method proposed by Phinney (1987), and vectorization is realized to improve the calculation speed.

[0053] (2) The pre-stack anisotropic inversion process based on geostatistics can obtain higher resolution inversion results and is applicable to tight sandstone reservoirs with thin reservoir development. Attached Figure Description

[0054] To more clearly illustrate the technical solutions of the embodiments of this application, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0055] Figure 1 A schematic flowchart of the thin reservoir geostatistical anisotropy pre-stack inversion method provided in the embodiments of this application;

[0056] Figure 2 The target region PP wave angle gather obtained by the thin reservoir geostatistical anisotropic pre-stack inversion method provided in the embodiments of this application;

[0057] Figure 3 The inversion results obtained by the thin reservoir geostatistical anisotropic pre-stack inversion method provided in this application embodiment are shown in black, representing the actual logging curves, blue the initial model, and red the inversion results. Figure 3 In the figure, (a) represents the longitudinal wave velocity. Figure 3 In the diagram, (b) represents the transverse wave velocity. Figure 3 In this context, (c) represents density. Figure 3 In this context, (d) represents the anisotropy parameter. Figure 3 In this context, (e) represents the anisotropy parameter;

[0058] Figures 4A-4E The inversion results obtained by the thin reservoir geostatistical anisotropic pre-stack inversion method provided in this application embodiment are connected to well profiles, where the red curve is the well logging curve. Figure 4A The inversion results representing the P-wave velocity are connected to the well profile. Figure 4B The inversion results representing the shear wave velocity are connected to the well profile. Figure 4C The density inversion results are linked to the well profile. Figure 4D This represents the well profile along with the inversion results of anisotropy parameters. Figure 4E The inversion results of anisotropy parameters are represented by well profiles.

[0059] Figure 5 A structural block diagram of the thin reservoir geostatistical anisotropy pre-stack inversion device provided in the embodiments of this application;

[0060] Figure 6 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0061] To better understand the technical solution of this application, the embodiments of this application will be described in detail below with reference to the accompanying drawings.

[0062] It should be understood that the described embodiments are merely some, not all, of the embodiments in this application. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

[0063] The terminology used in the embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. The singular forms “a,” “the,” and “the” used in the embodiments of this application and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise.

[0064] It should be understood that the term "and / or" used in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0065] In recent years, tight sandstone gas has become a major growth point for proven natural gas reserves in my country, with huge resource potential and broad exploration space. Exploration studies of tight sandstone gas in the Xujiahe Formation in the Chongzhou area of ​​the Sichuan Basin have revealed that these high-quality tight sandstone reservoirs are generally very thin and have developed horizontal fractures, exhibiting strong seismic anisotropy. This results in diverse internal structural characteristics, numerous combinations, and complex time-frequency response features, leading to the formation of complex converted waves. Consequently, formulas based on the single-interface assumption, such as the Zoeppritz equation and the Rüger formula, are no longer applicable for exploring the intrinsic response mechanisms and time-frequency variation patterns of this type of reservoir.

[0066] To address the aforementioned issues, this application provides a geostatistical pre-stack anisotropic inversion method, apparatus, and equipment for thin reservoirs. By equating tight sandstone thin reservoirs with horizontal fractures to thin interbedded media containing transversely isotropic media, it studies analytical solutions for reflection coefficients applicable to this type of medium, achieving rapid full-wavelength simulation. Based on this, it conducts research on geostatistical pre-stack stochastic inversion algorithms, fully utilizing full-wavelength information to predict reservoir parameters. This research will provide theoretical support for the prediction of tight sandstone thin reservoirs, solve some of the challenges in exploring tight sandstone thin reservoirs, and provide guidance for finding high-quality tight gas reservoir development areas. A detailed description is provided below with reference to the accompanying drawings.

[0067] See Figure 1This is a schematic flowchart of the pre-stack inversion method for thin reservoir geostatistical anisotropy provided in this application embodiment. Figure 1 As shown, it mainly includes the following steps:

[0068] Step S1: Determine the reflection coefficient for the VT1 single-layer model and the thin interlayer medium model with VT1 sandwiched between them.

[0069] For the VTl single-layer medium model, based on the relationship between displacement and stress at the two interfaces of the single-layer, and by analogy with the reflectivity method, we study the solution of the PP and PS wave reflection and transmission coefficient spectra of anisotropic single-layer media under arbitrary incident angles. For the thin interbedded medium model with VTI single-layer, which can be regarded as an approximation of a tight sandstone reservoir with low-angle fractures, it can be decomposed into thin interbedded layers of isotropic single-layer and VTl single-layer. The reflection coefficient of the anisotropic single-layer studied above is extended to the thin interbedded layers to analyze its overall response characteristics.

[0070] In the VT1 single-thin-layer dielectric model, the reflection coefficient formula for the frequency domain PP wave is:

[0071]

[0072] Among them, V 04 V 01 These are the first and fourth elements of array v0, which is a 6-element array defined in the slowness-frequency domain (p, ω), where p is the slowness and ω is the angular frequency. The expression for v0 is:

[0073] v0=[Δ -R PS Δ -R SS Δ R PP Δ R SP Δ detRΔ] T (2)

[0074] Where Δ represents the proportionality coefficient, R PS R SS R PP R SP denoted by and detR respectively, representing the total reflection coefficients for the corresponding wave types.

[0075] v0 can be calculated recursively:

[0076] V0 = Q0Q1…Q n …Q N-1 V N (3)

[0077] Among them, v N As the initial term, Q n Let be the propagation matrix of stratum n, where N is the total number of strata.

[0078] V N =[1 0 0 0 0 0] T (4)

[0079] Q n The calculation formula is:

[0080] Q n =E n F n (5)

[0081] Among them, E n F is the layer-crossing matrix, which describes the phase shift of a plane wave as it propagates through stratum n. n E is the interface matrix. n and F n It can be written specifically as:

[0082]

[0083] Where i is the imaginary unit, ω is the frequency, and d n Let be the thickness of the nth layer of the formation. Let be the vertical slowness of the P-wave in the nth formation medium. Let be the vertical slowness of the shear wave in the nth formation medium. Let T be the upward interface propagation matrix of the nth formation medium. n+1 The interface propagation matrix of the nth formation medium downwards.

[0084] In isotropic and VTI media And E n The form remains a matrix (2×2), that is, the form remains unchanged, so E in VTI medium n The form remains the same:

[0085]

[0086] F n Let F be the interface matrix. n It can be written specifically as:

[0087]

[0088] in, Let be the vertical slowness of the P-wave in the nth formation medium. Let S be the vertical slowness of the shear wave in the nth formation medium, and its specific expression is as follows:

[0089]

[0090] in,

[0091]

[0092] Among them, c ij Let ρ be the stiffness tensor, ρ be the medium density, and p be the horizontal slowness.

[0093] Substituting the Thomsen anisotropy parameters and The expression is obtained as follows:

[0094]

[0095]

[0096] Where p is the horizontal slowness, v p0 The longitudinal wave velocity in the vertical direction (axis of symmetry), v s0 ε is the transverse wave velocity in the vertical direction (axis of symmetry), δ and ε are the Thomsen anisotropy parameters, and θ is the incident angle.

[0097] Step S2: Based on geostatistical methods, the reflection coefficient is used to perform pre-stack anisotropic inversion of tight clastic thin reservoirs.

[0098] For tight sandstone thin reservoirs with low-angle fractures, a set of Kriging equations constrained by mixed observation data is constructed based on the reflection coefficients of thin interbedded media with VTI single thin layers derived above. This is used to perform anisotropic pre-stack inversion based on geostatistics, thereby estimating the elastic parameters of the target area (P-wave velocity α, S-wave velocity βn, medium density ρ, and anisotropic parameters δ and ε).

[0099] Typically, the forward modeling process can be represented as:

[0100] d=Ψ*R+n=G(m)+n (18)

[0101] Where d is the observation data vector, Ψ is the seismic wavelet matrix, R is the reflection coefficient, m is the model parameter vector, G is the nonlinear operator mapping m to d, G(m) is the result of forward modeling using the composite matrix method, and n is the error between the observation data vector and the forward modeling result.

[0102] The prior distribution of the model parameters is as follows:

[0103] P(m)=cρ(m)L(m) (19)

[0104] Where L(m) is the regularization term of the model parameters, and ρ(m) is the prior probability density function. Assuming that the noise is independent and follows a Gaussian distribution, the posterior distribution of the model parameters can be obtained using Bayes' theorem:

[0105]

[0106] Where P(d|m) is the likelihood function of the data; Con is the constant coefficient; m is the model parameter vector; m prior It is a 3n-dimensional column vector, composed of a smooth background model of P-wave velocities and densities, m prior =[lnα0,lnβ0,lnρ] o [lnε0, lnδ0] T C m For 3k m ×3k m The prior model covariance matrix of dimension:

[0107] C m (p, q) = C0v(h) (21)

[0108] Where h is the distance between any two points p and q in the model space, v(h) is the spatial structure covariance matrix, which mainly includes Gaussian model, spherical model, exponential model, and nested structures of various models, and C0 is the time-invariant covariance matrix.

[0109] The expression for the regularization term L(m) of the model parameters is:

[0110]

[0111] Where c2 is a constant coefficient; d obs For the observed data; G is the nonlinear operator mapping m to d; m prior It is a 3n-dimensional column vector, composed of a smooth background model of P-wave velocities and densities, m prior =[lnα0, lnβ0, lnρ0, lnε0, lnδ0] T C d is the seismic covariance matrix, estimated by combining well logging synthetic records and actual well bypass seismic data; m is the model parameter vector.

[0112] Substituting equations (20) and (22) into equation (19), we obtain the posterior solution space for the whole probability, where the posterior expectation is:

[0113]

[0114] And the posterior covariance is

[0115]

[0116] Within the theoretical framework of geostatistics, formulas (23) and (24) are simple Kriging estimates under the constraint of observational data, respectively. As Kriging expected, Let Kriging variance be denoted as .

[0117] In geostatistics, the basic idea of ​​sequential Gaussian simulation is to use simulated points as conditional data, together with the original known points, to participate in the calculation of the next point within a given neighborhood, until all nodes in the random path have been simulated. Borrowing from this idea, in sequential inversion, for ease of solution, the observed data within the neighborhood is divided into two categories, here assumed to be category A and category B. Category A data consists of direct measurements of model parameters, including well logging data and previously simulated points, also known as "hard data" or "point data"; Category B data consists of linearly averaged observed values ​​of model parameters, i.e., measurement data obtained through linear forward modeling operators, also known as "soft data" or "volume data," referring here to pre-stack seismic angle gather data.

[0118] Combining the two types of observation data, formula (18) can be rewritten as:

[0119]

[0120] Or, to simplify:

[0121] d obsA+B =G A+B m A+B +n A+B (26)

[0122] Where d obsA G A m A and n A These represent the observed data, forward modeling operators, model parameters, and error terms for class A data, respectively; d obsB G B m B and n B These represent the observed data, forward modeling operator, model parameters, and error term for type B data, respectively; d obsA+B G A+B m A+B and n A+B These represent the joint parameters corresponding to the two types of data, respectively. It is important to note here that this is related to the seismic forward modeling operator G. B Unlike the forward modeling operator G for directly observed type A data. A It is a simple identity matrix that does not depend on any physical laws.

[0123] Similar to the posterior solutions of formulas (23) and (24) for single data types, Bayesian inversion is performed on the joint data in formula (26) to obtain the posterior expectation under the joint constraints of well logging data (type A) and pre-stack seismic data (type B):

[0124]

[0125] and the posterior covariance is

[0126]

[0127] in, This represents the posterior expectation of well logging data (Type A) and pre-stack seismic data (Type B); Denotes the posterior covariance; m priorA+B This represents the initial pre-stack elastic parameters, which are (3k A k B )-dimensional column vector, k A and k B These represent the number of sampling points for well logging data (Category A) and pre-stack seismic data (Category B) within the neighborhood, respectively; C mA+B Indicates (3k) A ,k,k B )-dimensional model covariance matrix; G represents the transpose of the seismic forward calculus operator matrix; A+B Represents the seismic forward modeling operator matrix; C dA+B Indicates (3k) A ,k,k B )-dimensional data covariance matrix; d obsA+B This represents mixed observation data, which is (3k A k B A 1)-dimensional column vector; G represents the nonlinear operator that maps the model parameter vector m to the observation data vector d.

[0128] C mA+B and C dA+B The expression is:

[0129]

[0130] Among them, C mAA and C mBB C represents the model covariance matrix for well logging data (Type A) and pre-stack seismic data (Type B), respectively; mAB Represents the cross-covariance matrix of the two types of data; C represents the transpose of the cross-covariance matrix of the two classes of data; dBB This represents the data covariance matrix of pre-stack seismic data (Type B). The data covariance matrix C in the above formula... dA+B In this context, it is assumed that the directly observed logging data (Type A) are accurate values, and that the errors of the two types of data are uncorrelated. Therefore, the covariance matrix of the logging data (Type A) and the cross-covariance matrix of the two types of data are both 0, as shown in formula (29) 29. dA+B As shown.

[0131] The specific process for performing pre-stack anisotropic inversion of tight clastic thin reservoirs using the reflection coefficient based on geostatistical methods is as follows:

[0132] Acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times;

[0133] Define a random path that accesses all grid nodes in the model space;

[0134] Within the neighborhood of a node, well logging data and simulated point data are selected as class A data d. obsA Select angle gather data as class B data d obsB Constructing hybrid observation data d obsA+B ;

[0135] Based on the mixed observation data d obsA+B By combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance

[0136] The random realization value m of the current node is obtained by uniformly sampling the local Gaussian cumulative distribution function. sim (j) and add the random implementation value to the known point data;

[0137] The calculations are performed sequentially on all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result;

[0138] Regenerate random paths that can access all grid nodes in the model space, and perform anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementations of the model parameters m. sim Posterior expectation and posterior covariance

[0139] In one embodiment, the pre-stack anisotropic inversion of tight clastic thin reservoirs using the reflection coefficient based on geostatistical methods can be implemented through the following steps:

[0140] Step S1: Obtain well logging data, pre-stack angle gathers, seismic wavelet, initial smooth background model parameters, and random inversion number H;

[0141] Step S2: Define a random path that accesses all grid nodes in the model space, with Q nodes;

[0142] Step S3: Within the neighborhood of the currently accessed node, select well logging data and simulated point data as Class A data d. obsASelect angle gather data as class B data d obsB Constructing hybrid observation data d obsA+B ;

[0143] Step S4: Solve the Kriging equations under the constraints of the mixed observation data to obtain the local Kriging expectation of the current node. With variance

[0144] Step S5: Uniformly sample the local Gaussian cumulative distribution function to obtain the random realization value m of the current node. sim (j) and add the random implementation value to the known point data;

[0145] Step S6: Repeat steps S3-S5 until Q nodes in the random path have been calculated, and obtain one inversion result;

[0146] Step S7: Repeat steps S2-S6 until the set number of random inversions H is reached, to obtain multiple realization results m of the model parameters. sim Posterior expectation Posterior covariance

[0147] Taking the Xinchang tectonic belt in western Sichuan as an example, the target layer lithology is a combination of tight sandstone and mudstone interbedded, with the underground top generally buried at a depth of more than 4600m. The reservoir has relatively well-developed low-angle fractures, thin single-layer sand bodies, and a relatively flat stratigraphic structure, which can be equivalent to VTI medium.

[0148] See Figure 2 The PP wave angle gather for the target region is obtained by the thin reservoir geostatistical anisotropic pre-stack inversion method provided in this application embodiment. See also Figure 3 The figures represent the wellbore inversion results obtained using the thin reservoir geostatistical anisotropic pre-stack inversion method provided in this application embodiment. The black curve represents the actual logging curve, the blue curve represents the initial model, and the red curve represents the actual value. Figure 3 In the figure, (a) represents the longitudinal wave velocity. Figure 3 In the diagram, (b) represents the transverse wave velocity. Figure 3 In this context, (c) represents density. Figure 3 In this context, (d) represents the anisotropy parameter. Figure 3 In the figure, (e) represents the anisotropy parameter. It can be seen that the inversion results have high resolution and can basically reflect the changing trend of the curve, and also have significant identification and differentiation of thin layers. In addition, the inversion here requires that the initial model can basically reflect the changing trend of the strata.

[0149] The thin reservoir geostatistical anisotropic pre-stack inversion method disclosed in this application was used to perform inversion calculations on the entire work area.

[0150] See Figures 4A-4E The image shows the well profile obtained from the pre-stack inversion method for thin reservoir geostatistical anisotropy provided in this application embodiment, where the red curve represents the well logging curve. Figure 4A The inversion results representing the P-wave velocity are connected to the well profile. Figure 4B The inversion results representing the shear wave velocity are connected to the well profile. Figure 4C The density inversion results are linked to the well profile. Figure 4D This represents the well profile along with the inversion results of anisotropy parameters. Figure 4E This represents the well profile obtained from the inversion of anisotropic parameters. It's important to note that the basic assumption of the reflectivity method is that the data is local one-dimensional data or horizontal stratigraphic data. Combined with pre-stack time migration, this method is applicable to areas with simple geological conditions. In areas with complex geological conditions, pre-stack time-depth migration or reverse time migration is required. The data for this work area has undergone pre-stack migration processing, therefore it can be considered a horizontal stratigraphic profile. The inversion results show that they have high resolution and accurately reflect the trend of stratigraphic changes, which is very helpful in identifying effective tight sandstone reservoirs.

[0151] Corresponding to the above embodiments, this application also provides a thin reservoir geostatistical anisotropy pre-stack inversion device.

[0152] See Figure 5 This is a structural block diagram of a pre-stack inversion device for thin reservoir geostatistical anisotropy provided in an embodiment of this application. Figure 5 As shown, it mainly includes the following modules:

[0153] The reflection coefficient determination module 501 is used to determine the reflection coefficient for the VTI single thin layer model and the thin interlayer medium model with VTI single thin layer sandwiched in between.

[0154] The pre-stack inversion module 502 is used to perform anisotropic pre-stack inversion of tight clastic thin reservoirs based on geostatistical methods and using the reflection coefficient.

[0155] The pre-stack inversion module 502 includes the following modules:

[0156] The data acquisition module is used to acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times;

[0157] The path definition module is used to define a random path that accesses all grid nodes in the model space;

[0158] The data construction module is used to select well logging data and simulated point data as class A data within the neighborhood of a node. obsA Select angle gather data as class B data dobsB Constructing hybrid observation data d obsA+B ;

[0159] The node calculation module is used to calculate the mixed observation data d. obsA+B By combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance

[0160] The random implementation value acquisition module is used to uniformly sample the local Gaussian cumulative distribution function to obtain the random implementation value m of the current node. sim (j) and add the random implementation value to the known point data;

[0161] The result generation module is used to sequentially calculate all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result. The path definition module 502 regenerates a random path that can access all grid nodes in the model space, and performs anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementation results m of the model parameters. sim Posterior expectation and posterior covariance

[0162] It should be noted that the specific content involved in the embodiments of this application can be found in the description of the above method embodiments, and will not be repeated here for the sake of brevity.

[0163] Corresponding to the above embodiments, this application also provides an electronic device.

[0164] See Figure 6 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Figure 6 As shown, the electronic device 600 may include a processor 601, a memory 602, and a communication unit 603. These components communicate via one or more buses. Those skilled in the art will understand that the electronic device structure shown in the figures does not constitute a limitation on the embodiments of this application. It may be a bus topology or a star topology, and may include more or fewer components than shown, or combine certain components, or have different component arrangements.

[0165] The communication unit 603 is used to establish a communication channel, thereby enabling the electronic device to communicate with other devices.

[0166] The processor 601 serves as the control center of the electronic device, connecting various parts of the device via various interfaces and lines. It executes software programs and / or modules stored in the memory 602, and calls data stored in the memory to perform various functions and / or process data. The processor can be composed of integrated circuits (ICs), such as a single packaged IC or multiple packaged ICs with the same or different functions connected together. For example, the processor 601 may consist only of a central processing unit (CPU). In this embodiment, the CPU may have a single processing core or include multiple processing cores.

[0167] Memory 602 is used to store the execution instructions of processor 601. Memory 602 can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic storage, flash memory, magnetic disk or optical disk.

[0168] When the execution instructions in memory 602 are executed by processor 601, the electronic device 600 is able to perform some or all of the steps in the above method embodiments.

[0169] Corresponding to the above embodiments, this application also provides a computer-readable storage medium, wherein the computer-readable storage medium may store a program, wherein when the program runs, it can control the device where the computer-readable storage medium is located to execute some or all of the steps in the above method embodiments. Specifically, the computer-readable storage medium may be a magnetic disk, an optical disk, read-only memory (ROM), or random access memory (RAM), etc.

[0170] Corresponding to the above embodiments, this application also provides a computer program product containing executable instructions that, when executed on a computer, cause the computer to perform some or all of the steps in the above method embodiments.

[0171] In this application embodiment, "at least one" refers to one or more, and "more than one" refers to two or more. "And / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent the existence of A alone, the simultaneous existence of A and B, or the existence of B alone. A and B can be singular or plural. The character " / " generally indicates that the preceding and following related objects are in an "or" relationship. "At least one of the following" and similar expressions refer to any combination of these items, including any combination of single or plural items. For example, at least one of a, b, and c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be single or multiple.

[0172] Those skilled in the art will recognize that the units and algorithm steps described in the embodiments disclosed herein can be implemented using electronic hardware, computer software, or a combination of electronic hardware and software. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0173] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0174] In the several embodiments provided in this application, any function, if implemented as a software functional unit and sold or used as an independent product, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a part of the technical solution, 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 steps 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), random access memory (RAM), magnetic disks, or optical disks.

[0175] The above description is merely a specific embodiment of this application. 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 protection scope of this application. The protection scope of this application should be determined by the protection scope of the claims.

Claims

1. A pre-stack inversion method for geostatistical anisotropy of thin reservoirs, characterized in that, include: The reflection coefficients were determined for both the VTI single-layer model and the thin interlayer dielectric model with VTI single-layer sandwiched between them. Based on geostatistical methods, the reflection coefficient is used to perform pre-stack anisotropic inversion of tight clastic thin reservoirs.

2. The method according to claim 1, characterized in that, The formula for calculating the reflection coefficient is: Among them, v 04 v 01 These are the first and fourth elements of array v0, respectively.

3. The method according to claim 2, characterized in that, v0 is a 6-element array defined in the slow-frequency domain (p, ω). The expression for v0 is: v0=[△ -R PS △-R SS △ R PP △ R SP △ detR△] T Where △ represents the proportionality coefficient, R PS R SS R PP R SP denoted by and detR respectively, representing the total reflection coefficients for the corresponding wave types.

4. The method according to claim 1, characterized in that, The method based on geostatistics, utilizing the reflection coefficient, performs pre-stack anisotropic inversion of tight clastic thin reservoirs, including: Acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times; Define a random path that accesses all grid nodes in the model space; Within the neighborhood of a node, well logging data and simulated point data are selected as class A data d. obsA Select angle gather data as type B data d obsB Constructing hybrid observation data d obsA+B ; Based on the mixed observation data d obsA+B By combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance The random realization value m of the current node is obtained by uniformly sampling the local Gaussian cumulative distribution function. sim (j) and add the random implementation value to the known point data; The calculations are performed sequentially on all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result; Regenerate random paths that can access all grid nodes in the model space, and perform anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementations of the model parameters m. sim posterior expectation and posterior covariance 5. The method according to claim 4, characterized in that, Mixed observation data d obsA+B The calculation formula is as follows: d obsA+B =G A+B m A+B +n A+B Among them, G A+B The forward operator represents a mixture of data of type A and type B, m. A+B n represents the model parameters for mixed data of type A and type B. A+B This represents the error term for mixed data of type A and type B.

6. The method according to claim 5, characterized in that, The posterior expectation and posterior covariance The calculation formula is as follows: Where, m priorA+B C represents the initial pre-stack elasticity parameter. mA+B Let G represent the model covariance matrix. A+B The forward operator representing mixed data of type A and type B, d obsA+B Let G represent the mixed observation data, and let G represent the nonlinear operator that maps the model parameter vector m to the observation data vector d.

7. The method according to claim 6, characterized in that, Model covariance matrix C mA+B and data covariance matrix C dA+B The calculation formula is as follows: Among them, C mAA and C mBB C represents the model covariance matrix for class A data and class B data, respectively; mAB Represents the cross-covariance matrix of the two types of data; C represents the transpose of the cross-covariance matrix of the two classes of data; dBB The data covariance matrix representing class B data.

8. A thin reservoir geostatistical anisotropic pre-stack inversion device, characterized in that, include: The reflection coefficient determination module is used to determine the reflection coefficient for VTI single-layer models and thin interlayer media models with VTI single-layers. A pre-stack inversion module is used to perform anisotropic pre-stack inversion of tight clastic thin reservoirs based on geostatistical methods and utilizing the reflection coefficient; wherein, the pre-stack inversion module includes: The data acquisition module is used to acquire well logging data, pre-stack angle gathers, seismic wavelets, initial smooth background model parameters, and random inversion times; The path definition module is used to define a random path that accesses all grid nodes in the model space; The data construction module is used to select well logging data and simulated point data as class A data within the neighborhood of a node. obsA Select angle gather data as type B data d obsB Constructing hybrid observation data d obsA+B ; The node calculation module is used to calculate the mixed observation data d. obsA+B By combining well logging data, pre-stack angle gathers, seismic wavelets, and initial smooth background model parameters, the local kriging expectation of the current node is calculated. With variance The random implementation value acquisition module is used to uniformly sample the local Gaussian cumulative distribution function to obtain the random implementation value m of the current node. sim (j) and add the random implementation value to the known point data; The result generation module is used to sequentially calculate all nodes along the random path to obtain a single implementation of the anisotropic pre-stack inversion result. The path definition module 502 regenerates a random path that can access all grid nodes in the model space, and performs anisotropic pre-stack inversion again until the set number of random inversions is reached, obtaining multiple implementation results m of the model parameters. sim posterior expectation and posterior covariance 9. An electronic device, characterized in that, include: processor; Memory; And a computer program, wherein the computer program is stored in the memory, the computer program including instructions that, when executed by the processor, cause the electronic device to perform the method of any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored program, wherein, when the program is executed, it controls the device on which the computer-readable storage medium is located to perform the method according to any one of claims 1 to 7.

Images