Micro-fracture porosity seismic rock physics linear inversion method based on dual-dual pore medium

Abstract

The invention discloses a microfracture porosity seismic rock physics linear inversion method, and belongs to the field of oil and gas reservoir exploration. The method comprises the following steps: firstly, constructing an earthquake rock physical forward model by combining a dual-dual pore medium DDP-BR equation as a core and a Zoeppritz equation, establishing a quantitative mapping relation between physical parameters and earthquake data, calculating a Frechet derivative of a forward operator to the physical parameters, and realizing earthquake rock physical linear inversion based on DDP-BR by adopting an LBFGS-MT algorithm. The method comprises the following steps of: firstly, inverting hard hole local porosity, micro-crack local porosity, inclusion content and water saturation by adopting seismic data; and secondly, indirectly estimating the porosity of the hard hole, the porosity of the micro-crack and the total porosity, and avoiding the problem of inversion instability caused by strong coupling among parameters. A synthetic data test verifies that the method has good stability and inversion precision, and microfracture development information of a complex reservoir can be effectively estimated.

Description

A Seismic Rock Physics Linear Inversion Method for Microfracture Porosity Based on Dual-Dual Porosity Media Technical Field

[0001] This invention belongs to the field of oil and gas reservoir exploration and development technology, specifically relating to a seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media. Background Technology

[0002] In the exploration and evaluation of complex oil and gas reservoirs such as carbonate rocks and tight sandstones, seismic rock physics inversion is a core technical means to achieve quantitative prediction of reservoir physical parameters. This technology establishes the correlation between rock elastic parameters and reservoir physical properties by constructing a rock physics model, and then obtains reservoir parameters by combining seismic data inversion, providing key technical support for the formulation of oil and gas reservoir development plans.

[0003] Current seismic rock physics inversion methods mostly employ rock physics models based on the assumption of a homogeneous medium (such as the Gassmann equation). However, complex reservoirs commonly exhibit heterogeneous porous structures such as microfractures and dissolution pores, often accompanied by patchy distributions of multiphase fluids. This makes it difficult for traditional homogeneous medium models to accurately characterize wave propagation patterns and establish reliable mappings between physical properties and elastic parameters, thus limiting the improvement of inversion prediction accuracy. To address the problem of characterizing reservoir heterogeneity, dual-porosity medium models (such as the Biot-Rayleigh, or BR model) have been gradually introduced into this field. Subsequently, the equations have been extended from dual-porosity media saturated with a single fluid to complex cases that simultaneously consider dual-pore structures and partially saturated two-phase fluids, resulting in the double-double-porosity (DDP) structural model. This model can effectively characterize the heterogeneous features of the rock skeleton, significantly overcoming the limitations of traditional homogeneous medium models.

[0004] When conducting seismic prediction studies of reservoir physical parameters based on the BR model, nonlinear inversion methods are typically employed. However, due to the inherently high nonlinearity of the BR model, these methods suffer from low computational efficiency, making them unsuitable for the efficient processing of large-area 3D seismic data. Meanwhile, microfractures, as the primary channels for fluid seepage in complex reservoirs, directly impact the reliability of reserve assessment and the rationality of development scheme design. Traditional models fail to effectively characterize the development features of microfractures, often misinterpreting seismic response anomalies caused by microfractures as changes in fluid type, further reducing the accuracy of reservoir evaluation. Furthermore, while existing linear inversion methods offer high computational efficiency, they often rely on the linearization of simple rock physics models, a process that introduces errors and limits accuracy. Linear methods impose certain requirements on the rock physics model: its mathematical form must be relatively simple to facilitate linearization derivation. However, overly simplistic models often fail to effectively describe the actual formation conditions, especially in complex reservoirs. Therefore, this study investigates a direct linear seismic inversion method for microfracture porosity based on dual-dual-porosity media, employing the central difference method to solve partial derivatives and reduce errors. It can simultaneously take into account the dual heterogeneity of pore structure and fluid distribution, achieving a synergistic improvement in inversion accuracy and computational efficiency, thereby meeting the practical engineering needs of complex reservoir exploration and evaluation. Summary of the Invention

[0005] This invention addresses the problems existing in the prior art by providing a seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media, the steps of which are as follows:

[0006] Step 1: Obtain pre-stack seismic data for the study area and extract seismic wavelets from different angles;

[0007] Step 2: Obtain reservoir logging data and estimate key wellhead parameters, such as hard pore local porosity, using rock physics modeling methods. Microcrack local porosity Inclusion content V inc and water saturation S w Establish an initial parameter model;

[0008] Step 3: Based on the initial model, calculate the elastic parameters, including the longitudinal wave velocity V, using the Biot-Rayleigh model for dual-dual-porosity media. p Shear wave velocity V s and density Based on the calculation results of elastic parameters, the synthetic record is obtained, and the residual between the synthetic seismic data and the actual data is calculated.

[0009] Step 4: Based on the initial model parameters, use the central difference method to calculate the local porosity of hard pores using the seismic rock physics forward modeling operator. Microcrack local porosity Inclusion content V inc and water saturation S w The partial derivatives;

[0010] Step 5: Establish the inversion objective function, calculate the first derivative of the objective function and the model update gradient, and obtain the model update result using the iterative format;

[0011] Step 6: Return to steps 3 to 6 for the next inversion iteration. Use the model parameters to obtain the synthetic record. When the data residuals decrease to the preset range, the iteration stops, and the model property parameter inversion results are output: local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation S w ;

[0012] Step 7: Using the physical property parameter inversion results output in Step 6, and combining them with the porosity mathematical formula, obtain the final physical property parameters: hard pore porosity. Microcrack porosity and total porosity Estimation results.

[0013] Furthermore, in step 1 above, the pre-stack offset data of PP waves in the work area is obtained. Based on the velocity data, the pre-stack offset data is converted into pre-stack angle gathers, with the angle range of the angle gathers being 0°-40°. Seismic wavelet sequences are extracted from the pre-stack angle gathers, and wavelet matrices W(θ) are established at different angles.

[0014] Furthermore, step 2 mentioned above specifically includes:

[0015] Acquire logging data curves at the wellhead, including P-wave velocity V. p Shear wave velocity V s Density ρ, Total porosity Water saturation S w and mud content V sh Rock physics modeling methods are used to estimate the local porosity of hard pores in the well. Microcrack local porosity ϕ 20 Inclusion content V inc ;Will , , and S w Set the parameters to be inverted and establish the inversion target vector z;

[0016] An initial model for the parameters to be inverted is established, based on parameter estimates obtained from rock physics modeling, including local porosity in hard pores. Microcrack local porosity Inclusion content V inc and water saturation First, the Backus averaging algorithm is used to obtain the low-frequency trend of the target parameters. Then, the Kriging interpolation and layer interpretation results are used to obtain the 2D / 3D low-frequency initial model.

[0017] Furthermore, step 3 described above includes the following sub-steps:

[0018] Step 3.1: Based on the dual-dual-porosity media DDP-BR model, the forward modeling of rock physical parameters z to elastic parameters m is as follows:

[0019] ,

[0020] Among them, the elastic parameter m includes the longitudinal wave velocity V p Shear wave velocity V s And density ρ, z includes physical property parameters and auxiliary parameters, among which the physical property parameters involve the local porosity of hard pores. Microcrack local porosity Fluid saturation and inclusion content The parameter S represents the rock physics modeling process based on the DDP-BR model. The specific form of the DDP-BR model is as follows:

[0021] ,

[0022] ,

[0023] ,

[0024] ,

[0025] ,

[0026] ,

[0027] ,

[0028] ,

[0029] Where u is the particle displacement vector of the solid skeleton. and Let ξ represent the first and second derivatives of u, respectively; (1) ξ (2) ξ (3) and ξ (4) U corresponds to the volumetric strain of the four types of pores. (1)U (2) U (3) and U (4) Let ζ be four independent fluid displacement vectors, corresponding to four different porosities. 12 ζ 13 and ζ 24 For the introduced fluid increment; , , and Represents absolute porosity; local porosity in hard pores and microcracks is respectively expressed as... and Indicated; κ1 and κ2 represent the permeability of the matrix framework and inclusion framework, respectively, R 12 R is another skeleton radius embedded in the principal phase skeleton. 13 R 24 These are the radii of the principal phase skeleton and the inclusion skeleton, respectively. and ρ represents the density and viscosity of the matrix fluid. 00 ρ 01 ρ 02 ρ 03 ρ 04 ρ 11 ρ 22 ρ 33 ρ 44 There are nine density coefficients, and A, N, Q1, Q2, Q3, Q4, R1, R2, R3, and R4 are ten stiffness coefficients, all of which are intermediate values.

[0030] Step 3.2: Solve the DDP-BR equations using the plane wave analysis method to obtain the P-wave velocity V. p With transverse wave velocity V s :

[0031] , ,

[0032] Where ω is the angular frequency, μ m Let ρ be the shear modulus of the rock skeleton, and ρ be the density of the rock in saturated fluid. Indicates taking the real part;

[0033] Step 3.3: Construct a forward model of seismic rock physics using the exact Zoeppritz equation:

[0034] ,

[0035] Where d represents pre-stack seismic data, W is the wavelet matrix, r is the reflection coefficient calculated by the exact Zoeppritz equation, and G is the forward process of seismic data calculated by the exact Zoeppritz equation.

[0036] Step 3.4: Calculate the reflection coefficient sequence r, as shown in the following expression:

[0037] ,

[0038] The matrices A and b are expressed as follows:

[0039] ,

[0040] ,

[0041] Among them, V P1 V S1 ρ1 represents the longitudinal wave velocity, transverse wave velocity, and density V of the medium in the upper half-space. P2 V S2 , The values ​​represent the longitudinal wave velocity, transverse wave velocity, and density of the medium in the lower half-space, and i, i', j, and j' represent the longitudinal wave incident angle, longitudinal wave transmission angle, converted transverse wave reflection angle, and converted transverse wave transmission angle, respectively.

[0042] Step 3.5: Based on convolution theory, using the wavelet matrix W and the reflection coefficient sequence r, calculate the synthetic seismic data d. syn :

[0043] ,

[0044] Where G(z) represents the forward modeling of seismic rock physics, and z is the parameter vector to be inverted, including the local porosity of hard pores. Microcrack local porosity Inclusion content and water saturation ; Forward modeling composite data, For wavelet matrix,

[0045] Calculate the actual P-wave dataset d obs and forward simulation records residual :

[0046] .

[0047] Furthermore, in step 3.5 above, d syn It has the following matrix form:

[0048] ,

[0049] In the formula, Angle of incidence The corresponding synthetic data is as follows:

[0050] ,

[0051] Among them, t 1…… t M Representing the tth 1…… t M Each depth time sample point, from the incident angle Wavelet matrices are constructed by extracting wavelets from actual seismic traces. ; Angle of incidence The reflection coefficient sequence obtained by time calculation is used to obtain:

[0052] .

[0053] Furthermore, in step 4 above, based on the parameter values ​​of the initial logging model, the partial derivatives of the seismic rock physics forward model are calculated using the central difference method; Let G(z) be the derivative of the forward operator with respect to the model parameter z, and take the following form:

[0054] ,

[0055] in, Indicates the angle of incidence θ i The corresponding partial derivative matrix, where = 1, 2; This represents the Na-th incident angle.

[0056] Reflection coefficient affects the local porosity of hard pores The partial derivatives are calculated by the following formula:

[0057] ,

[0058] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0059] ,

[0060] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t kLocal porosity of hard pores The partial derivatives, Local porosity of hard pores Tiny variables, Representing a time depth of t j time value;

[0061] Reflection coefficient affects the local porosity of microcracks The partial derivatives are calculated by the following formula:

[0062] ,

[0063] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0064] ,

[0065] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Local porosity of microcracks The partial derivatives, Local porosity of microcracks Tiny variables, Representing a time depth of t j time value.

[0066] Reflectance coefficient affects inclusion content V inc The partial derivatives are calculated by the following formula:

[0067] ,

[0068] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to V inc The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0069] ,

[0070] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Inclusion content V at the locationinc The partial derivatives, V, the content of inclusions inc Tiny variables, Representing a time depth of t j V at time inc value.

[0071] Reflectance coefficient for water saturation S w The partial derivatives are calculated by the following formula:

[0072] ,

[0073] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to S w The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0074] ,

[0075] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Water saturation S at the location w The partial derivatives, Water saturation S w Tiny variables, Representing a time depth of t j S at time w value.

[0076] Furthermore, step 5 mentioned above includes the following sub-steps:

[0077] Step 5.1: Establish the inversion objective function L(z), the specific expression of which is:

[0078] ,

[0079] Among them, C z Let be the covariance matrix of the model parameters, u be the expectation of the model parameters, and λ be the regularization parameter;

[0080] Step 5.2: Set the negative gradient as the update direction of the model. According to the LBFGS-MT algorithm, the iterative update formula for the model parameters is:

[0081] ,

[0082] ,

[0083] Among them, vk v represents the update amount in the k-th iteration. k-1 This represents the update amount in the (k-1)th iteration. Z is the momentum coefficient. k Let z be the model parameter in the k-th iteration. k-1 Let κ be the model parameters in the (k-1)th iteration. k It is the step size of the k-th iteration.

[0084] The specific expression for the first derivative J(z) of the objective function with respect to the model parameters z is as follows:

[0085] ,

[0086] Step 5.3: Let H(z) represent the pseudo-Hessian matrix of the objective function. Using the Gauss-Newton optimization algorithm, the pseudo-Hessian matrix takes the form:

[0087] ,

[0088] Where k is the number of iterations. ε k and S k All of these are intermediate variables, and their specific forms are as follows:

[0089] , , ,

[0090] Where I represents the identity matrix, y k It is also an intermediate variable, in the following form:

[0091] ,

[0092] Among them, J(z) k ) and J(z k+1 ) represent the Jacobian matrix values ​​for the k-th and (k+1)-th iterations, respectively;

[0093] Step 5.4: Set the initial pseudo-Hessian matrix to the following form:

[0094] ,

[0095] The pseudo-Hessian matrix H in the (k+1)th iteration k+1 (z), in the following specific form:

[0096] ,

[0097] Where M is a constant.

[0098] Furthermore, in step 5.2 above, the covariance matrix C of the model parameters zThe specific form is as follows:

[0099] ,

[0100] In the formula, σ is an M×M dimension matrix.

[0101] Furthermore, in step 7 above, calculating the physical property parameters using the parameter inversion results from step 6 includes:

[0102] The parameter inversion results in step 6 include the local porosity of hard pores. Local porosity of microcracks Inclusion content V inc and water saturation S w The final physical property parameters can be calculated using mathematical relationships. The specific mathematical expression is as follows:

[0103] , ,

[0104] ,

[0105] in, For hard pore porosity, For microcrack porosity, Represents total porosity.

[0106] Compared with the prior art, the beneficial technical effects of the present invention using the above technical solution are as follows:

[0107] Previous reservoir property parameter inversion methods employed rock physics theories assuming homogeneous media, such as the Gassmann model. These methods fail to account for the heterogeneity caused by complex pore structures or patchy fluids in complex reservoirs, limiting the accuracy of reservoir predictions. This invention introduces the Biot-Rayleigh model, a dual-porosity medium model, into reservoir property parameter inversion to improve the simulation accuracy of rock physics forward modeling. Furthermore, studies based on the BR equation often employ nonlinear inversion methods; however, due to the inherently high nonlinearity of the BR model, these methods have low computational efficiency and are unsuitable for the efficient processing requirements of large-area 3D seismic data. To improve efficiency while maintaining accuracy, this invention investigates a direct linear seismic inversion method for microfracture porosity based on dual-porosity media, using the central difference method to solve partial derivatives and reduce errors. Compared to conventional homogeneous medium property parameter inversion, this invention offers higher inversion accuracy and is more suitable for predicting complex heterogeneous reservoirs. Attached Figure Description

[0108] Figure 1 is a flowchart of the present invention.

[0109] Figure 2 is a schematic diagram of well logging model data converted to seismic scale.

[0110] Figure 3 is a schematic diagram of the inversion results based on the Gassmann inversion method.

[0111] Figure 4 shows the results of the first step inversion and indirect estimation obtained based on the DDP-BR inversion method; in the figure, (a) is the result of the first step inversion. , V inc and S w Results; (b) are indirectly estimated based on the results of the first step inversion. , as well as result.

[0112] Figure 5 shows the first-step inversion results and indirect estimation results based on the DDP-BR inversion method under different signal-to-noise ratios;

[0113] In the figures, (a) and (b) show the first-step inversion result and the second-step indirect estimation result for a signal-to-noise ratio of 50, respectively; (c) and (d) show the first-step inversion result and the second-step indirect estimation result for a signal-to-noise ratio of 10, respectively; and (e) and (f) show the first-step inversion result and the second-step indirect estimation result for a signal-to-noise ratio of 3, respectively. Detailed Implementation

[0114] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0115] In this invention, various aspects of the invention are described with reference to the accompanying drawings, in which numerous illustrative embodiments are shown. Embodiments of the invention are not limited to those depicted in the drawings. It should be understood that the invention is implemented through any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0116] As shown in Figure 1, this invention provides a seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media, the steps of which are as follows:

[0117] Step 1: Obtain pre-stack seismic data for the study area and extract seismic wavelets from different angles;

[0118] Step 2: Obtain reservoir logging data and estimate key wellhead parameters, such as hard pore local porosity, using rock physics modeling methods. Microcrack local porosity Inclusion content V incand water saturation S w Establish an initial parameter model;

[0119] Step 3: Based on the initial model, calculate the elastic parameters, including the longitudinal wave velocity V, using the Biot-Rayleigh model for dual-dual-porosity media. p Shear wave velocity V s and density Based on the calculation results of elastic parameters, the synthetic record is obtained, and the residual between the synthetic seismic data and the actual data is calculated.

[0120] Step 4: Based on the initial model parameters, use the central difference method to calculate the local porosity of hard pores using the seismic rock physics forward modeling operator. Microcrack local porosity Inclusion content V inc and water saturation S w The partial derivatives;

[0121] Step 5: Establish the inversion objective function, calculate the first derivative of the objective function and the model update gradient, and obtain the model update result using the iterative format;

[0122] Step 6: Return to steps 3 to 6 for the next inversion iteration. Use the model parameters to obtain the synthetic record. When the data residuals decrease to the preset range, the iteration stops, and the model parameter results are output: local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation S w ;

[0123] Step 7: Using the physical property parameter inversion results output in Step 6, and combining them with the porosity mathematical formula, obtain the final physical property parameters: hard pore porosity. Microcrack porosity and total porosity Estimation results.

[0124] In a preferred embodiment of the present invention, step 1 specifically involves: obtaining the pre-stack offset gather of the PP data of the work area, and converting the pre-stack offset data into a pre-stack angle gather based on the 3D velocity data;

[0125] Seismic wavelet sequences were extracted from pre-stack trace sets to establish wavelet matrices at different angles. It is an augmented matrix:

[0126] ,

[0127] in, Representative angle is The wavelet matrix at time is composed of wavelet sequences:

[0128] .

[0129] As a preferred embodiment of the present invention, as shown in Figure 2, the logging data curve at the wellhead is obtained, including the longitudinal wave velocity V. p Shear wave velocity V s ,density Total porosity Water saturation S w and mud content V sh Rock physics modeling methods are used to realize the local porosity of hard pores at the wellhead. Microcrack local porosity and inclusion content V inc Estimating physical property parameters provides prior information for subsequent inversion. This includes estimating the local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation S w The parameters to be inverted are set as follows, and the inversion target parameter vector z is:

[0130] ,

[0131] in, , … Representing the first to the Mth time samples value, , … Representing the first to the Mth time samples Value, (V) inc )1, (V inc )2…(V inc ) M V represents the first to the Mth time samples. inc Value, (S) w )1, (S w )2…(S w ) M S represents the first to the Mth time samples. w value.

[0132] An initial model for the parameters to be inverted is established, based on parameter estimates obtained from rock physics modeling, including local porosity in hard pores. 、 Local porosity of microcracks and inclusion content V incFirst, the Backus averaging algorithm is used to obtain the low-frequency trend of the target parameters. Then, the Kriging interpolation and layer interpretation results are used to obtain the 2D / 3D low-frequency initial model.

[0133] In a preferred embodiment of the present invention, step 3 includes the following sub-steps:

[0134] Step 3.1: Based on the Biot-Rayleigh (DDP-BR) model for dual-dual-porosity media, the forward modeling of rock physical parameters z to elastic parameters m is as follows:

[0135] ,

[0136] Among them, the elastic parameter m includes the longitudinal wave velocity V p Shear wave velocity V s And density ρ, z includes physical property parameters and auxiliary parameters, among which the physical property parameters involve the local porosity of hard pores. Microcrack local porosity Fluid saturation and inclusion content V inc The parameter S represents the rock physics modeling process based on the DDP-BR model. The specific form of the DDP-BR model is as follows:

[0137] ,

[0138] ,

[0139] ,

[0140] ,

[0141] ,

[0142] ,

[0143] ,

[0144] ,

[0145] Where u is the particle displacement vector of the solid skeleton. and Let ξ represent the first and second derivatives of u, respectively; (1) ξ (2) ξ (3) and ξ (4) U corresponds to the volumetric strain of the four types of pores. (1) U (2) U (3) and U(4) Let ζ be four independent fluid displacement vectors, corresponding to four different porosities. 12 ζ 13 and ζ 24 For the introduced fluid increment; , , and Represents absolute porosity; local porosity in hard pores and microcracks is respectively expressed as... and Indicated; κ1 and κ2 represent the permeability of the matrix framework and inclusion framework, respectively, R 12 R is another skeleton radius embedded in the principal phase skeleton. 13 R 24 The radii (R) of the principal phase skeleton and the inclusion skeleton, respectively. 24 < <R 12 ), and ρ represents the density and viscosity of the matrix fluid. 00 ρ 01 ρ 02 ρ 03 ρ 04 ρ 11 ρ 22 ρ 33 ρ 44 There are nine density coefficients, and A, N, Q1, Q2, Q3, Q4, R1, R2, R3, and R4 are ten stiffness coefficients, all of which are intermediate values.

[0146] Step 3.2: Solve the DDP-BR equations using the plane wave analysis method to obtain the P-wave velocity V. p With transverse wave velocity V s :

[0147] , ,

[0148] Where ω is the angular frequency, μ m Let ρ be the shear modulus of the rock skeleton, and ρ be the density of the rock in saturated fluid. Indicates taking the real part;

[0149] Step 3.3: Construct a forward model of seismic rock physics using the exact Zoeppritz equation:

[0150] ,

[0151] Where d represents pre-stack seismic data, W is the wavelet matrix, r is the reflection coefficient calculated by the exact Zoeppritz equation, and G is the forward process of seismic data calculated by the exact Zoeppritz equation.

[0152] Step 3.4: Calculate the reflection coefficient sequence r, as shown in the following expression:

[0153] ,

[0154] The matrices A and b are expressed as follows:

[0155] ,

[0156] ,

[0157] Among them, V P1 V S1 ρ1 represents the longitudinal wave velocity, transverse wave velocity, and density V of the medium in the upper half-space. P2 V S2 ρ2 represents the longitudinal wave velocity, transverse wave velocity, and density of the medium in the lower half-space, and i, i', j, and j' represent the longitudinal wave incident angle, longitudinal wave transmission angle, converted transverse wave reflection angle, and converted transverse wave transmission angle, respectively.

[0158] Step 3.5: Based on convolution theory, using the wavelet matrix W and the reflection coefficient sequence r, calculate the synthetic seismic data d. syn :

[0159] ,

[0160] Where G(z) represents the forward modeling of seismic rock physics, and z is the parameter vector to be inverted, including the local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation S w ; Forward modeling composite data, For wavelet matrix,

[0161] Calculate the actual P-wave dataset d obs and forward simulation records residual :

[0162] .

[0163] As a preferred embodiment of the present invention, the partial derivatives of the seismic rock physics forward model are calculated using the central difference method based on the parameter values ​​of the initial well logging model. Let G(z) be the derivative of the forward operator with respect to the model parameter z, and take the following form:

[0164] ,

[0165] Indicates the angle of incidence θ i (in = 1, 2,) corresponding partial derivative matrices, This represents the reflection coefficient at the Na-th incident angle for the local porosity of the hard pore. The partial derivatives can be calculated using the following formula:

[0166] ,

[0167] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0168] ,

[0169] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Local porosity of hard pores The partial derivatives, Local porosity of hard pores Tiny variables, Representing a time depth of t j time value.

[0170] Reflection coefficient affects the local porosity of microcracks The partial derivatives can be calculated using the following formula:

[0171] ,

[0172] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0173] ,

[0174] in, This represents the angle of incidence being θ. i The time depth is tj Time depth is t k Local porosity of microcracks The partial derivatives, Local porosity of microcracks Tiny variables, Representing a time depth of t j time value.

[0175] Reflectance coefficient affects inclusion content V inc The partial derivatives can be calculated using the following formula:

[0176] ,

[0177] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to V inc The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0178] ,

[0179] in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Inclusion content V at the location inc The partial derivatives, V, the content of inclusions inc Tiny variables, Representing a time depth of t j V at time inc value.

[0180] Reflectance coefficient for water saturation S w The partial derivatives can be calculated using the following formula:

[0181] ,

[0182] in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to S w The partial derivatives are calculated using the central difference method. The partial derivatives are given by the following formula:

[0183] ,

[0184] in, This represents the angle of incidence being θ. i The time depth is t jTime depth is t k Water saturation S at the location w The partial derivatives, Water saturation S w Tiny variables, Representing a time depth of t j S at time w value.

[0185] In a preferred embodiment of the present invention, an inversion objective function is established, the first derivative of the objective function and the model update gradient are calculated, and the model update result is obtained using an iterative format; the inversion objective function L(z) is established, and its specific expression is as follows:

[0186] ,

[0187] Where u is the expectation of the model parameters, λ is the regularization parameter, and C z The covariance matrix of the model parameters has the following specific form:

[0188] ,

[0189] In the formula, σ is an M×M dimension matrix, with... For example, it has the following form:

[0190] ,

[0191] Assuming the negative gradient as the model update direction, and based on the LBFGS-MT (LBFGS combined with momentum techniques) algorithm, the iterative update formula for the model parameters is:

[0192] ,

[0193]

[0194] Among them, v k v represents the update amount in the k-th iteration. k-1 This represents the update amount in the (k-1)th iteration. Z is the momentum coefficient. k Let z be the model parameter in the k-th iteration. k-1 Let κ be the model parameters in the (k-1)th iteration. k It is the step size of the k-th iteration.

[0195] The specific expression for the first derivative J(z) of the objective function with respect to the model parameters z is as follows:

[0196] ,

[0197] H(z) represents the pseudo-Hessian matrix of the objective function. Using the Gauss-Newton optimization algorithm, the pseudo-Hessian matrix takes the following form:

[0198] ,

[0199] Where k is the number of iterations. ε k and S k All of these are intermediate variables, and their specific forms are as follows:

[0200] , , ,

[0201] Where I represents the identity matrix, y k It is also an intermediate variable, in the following form:

[0202] .

[0203] The initial pseudo-Hessian matrix is ​​set to the following form:

[0204] ,

[0205] The pseudo-Hessian matrix for the (k+1)th iteration is as follows:

[0206] .

[0207] In the formula, M is a constant.

[0208] In a preferred embodiment of the present invention, the final estimated physical property parameters are obtained by combining the results of the first step of physical property parameter inversion with the mathematical relationship of porosity. The specific mathematical expression is as follows:

[0209] , ,

[0210] .

[0211] in, For hard pore porosity, For microcrack porosity, Represents total porosity.

[0212] As shown in Figure 3, previous seismic rock physics inversion methods mostly adopted the Gassmann equation based on the homogeneous medium assumption, which could not account for the heterogeneity of pore structure and fluids. Furthermore, nonlinear inversion methods have low computational efficiency and are difficult to apply to actual 3D seismic data. Figures 4-5 show the first-step inversion results and indirect estimation results obtained based on the DDP-BR inversion method. Figure 4(a) shows the results of the first-step inversion. , V inc and S w As a result, (b) shows the indirect estimation based on the inversion results of the first step. , as well as Results. Figure 5 shows the first-step inversion results and indirect estimation results of the DDP-BR inversion method under different signal-to-noise ratios; (a) and (b) show the first-step inversion results and the second-step indirect estimation results with an signal-to-noise ratio of 50, respectively; (c) and (d) show the first-step inversion results and the second-step indirect estimation results with an signal-to-noise ratio of 10, respectively; and (e) and (f) show the first-step inversion results and the second-step indirect estimation results with an signal-to-noise ratio of 3, respectively.

[0213] In summary, the seismic rock physics linear inversion method for microfracture porosity based on the dual-dual-porosity Biot-Rayleigh (DDP-BR) model of this invention aims to solve the problems of traditional inversion methods, such as difficulty in adapting to the heterogeneity of complex reservoirs (e.g., carbonate rocks, tight sandstones), low computational efficiency, and poor result stability. This method first constructs a seismic rock physics forward model based on the DDP-BR equations and the exact Zoeppritz equations, simultaneously describing the dual heterogeneity of pore structure and fluid distribution, and establishing a quantitative mapping relationship between rock physical parameters and pre-stack seismic data. Second, addressing the difficulty of derivative calculation caused by the high nonlinearity of the DDP-BR model, the central difference method is used to calculate the Fréchet derivatives of the forward modeling operators with respect to each rock physics parameter, and the LBFGS-MT algorithm is combined to achieve seismic rock physics linear inversion based on the DDP-BR model. A step-by-step inversion strategy is proposed: the first step involves inverting the local porosity of hard pores, the local porosity of microcracks, the inclusion content, and the water saturation; the second step uses the inversion results from the first step to indirectly estimate the porosity of hard pores, the porosity of microcracks, and the total porosity, thus avoiding instability caused by strong coupling between parameters. Synthetic data testing verifies the feasibility of the proposed inversion method.

[0214] While the present invention has been described above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A linear inversion method for seismic rock physics of microfracture porosity based on dual-dual-pore media, characterized in that, The steps are as follows: Step 1: Obtain pre-stack seismic data for the study area and extract seismic wavelets at different angles; Step 2: Obtain reservoir logging data and estimate the key wellhead parameter, hard pore local porosity, using rock physics modeling methods. Microcrack local porosity Inclusion content V inc and water saturation S w Step 3: Based on the initial model, calculate the elastic parameters, including the longitudinal wave velocity V, using the Biot-Rayleigh model for dual-dual-porosity media. p Shear wave velocity V s and density Step 4: Based on the elastic parameter calculation results, obtain the synthetic record and calculate the residual between the synthetic seismic data and the actual data; based on the initial model parameters, use the central difference method to calculate the local porosity of hard pores of the seismic rock physics forward modeling operator. Microcrack local porosity Inclusion content V inc and water saturation S w Partial derivatives; Step 5: Establish the inversion objective function, calculate the first derivative of the objective function and the model update gradient, and obtain the model update result using the iterative format; Step 6: Return to execute steps 3 to 6 for the next inversion iteration, obtain the synthetic record using the model parameters, and stop the iteration when the data residuals drop to the preset range, outputting the model physical property parameter inversion result: local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation S w Step 7: Using the physical property parameter inversion results output in Step 6, combined with the porosity mathematical relationship, obtain the final physical property parameters: hard pore porosity. Microcrack porosity and total porosity Estimation results.

2. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 1, characterized in that, In step 1, pre-stack offset data of PP waves in the work area are obtained. Based on the velocity data, the pre-stack offset data is converted into pre-stack angle gathers with an angle range of 0°-40°. Seismic wavelet sequences are extracted from the pre-stack angle gathers, and wavelet matrices W(θ) are established at different angles.

3. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 1, characterized in that, Step 2 specifically involves acquiring the logging data curves at the wellhead, including the P-wave velocity V. p Shear wave velocity V s Density ρ, Total porosity Water saturation S w and mud content V sh Rock physics modeling methods are used to estimate the local porosity of hard pores in the well. Microcrack local porosity Inclusion content V inc ;Will 、 、 and S w Set the parameters to be inverted and establish the inversion target vector z; establish an initial model for the parameters to be inverted, based on the parameter results estimated by rock physics modeling, including the local porosity of hard pores. Microcrack local porosity Inclusion content V inc and water saturation First, the Backus averaging algorithm is used to obtain the low-frequency trend of the target parameters. Then, the Kriging interpolation and layer interpretation results are used to obtain the 2D / 3D low-frequency initial model.

4. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 1, characterized in that, Step 3 includes the following sub-steps: Step 3.1, Based on the dual-dual porosity media DDP-BR model, the rock physical forward modeling from rock physical parameters z to elastic parameters m is expressed as follows: The elastic parameter m includes the longitudinal wave velocity V. p Shear wave velocity V s And density ρ, z includes physical property parameters and auxiliary parameters, among which the physical property parameters involve the local porosity of hard pores. Microcrack local porosity Fluid saturation and inclusion content The parameter S represents the rock physics modeling process based on the DDP-BR model. The specific form of the DDP-BR model is as follows: ， ， ， ， ， ， ， Where u is the particle displacement vector of the solid skeleton. and Let ξ represent the first and second derivatives of u, respectively; (1) ξ (2) ξ (3) and ξ (4) U corresponds to the volumetric strain of the four types of pores. (1) U (2) U (3) and U (4) Let ζ be four independent fluid displacement vectors, corresponding to four different porosities. 12 ζ 13 and ζ 24 For the introduced fluid increment; 、 、 and Represents absolute porosity; local porosity in hard pores and microcracks is respectively expressed as... and Indicated; κ1 and κ2 represent the permeability of the matrix framework and inclusion framework, respectively, R 12 R is another skeleton radius embedded in the principal phase skeleton. 13 R 24 These are the radii of the principal phase skeleton and the inclusion skeleton, respectively. and ρ represents the density and viscosity of the matrix fluid. 00 ρ 01 ρ 02 ρ 03 ρ 04 ρ 11 ρ 22 ρ 33 ρ 44 There are nine density coefficients, and A, N, Q1, Q2, Q3, Q4, R1, R2, R3, and R4 are ten stiffness coefficients, all of which are intermediate values; Step 3.2: Solve the DDP-BR equation using the plane wave analysis method to obtain the longitudinal wave velocity V. p With transverse wave velocity V s : ， Where ω is the angular frequency, μ m Let ρ be the shear modulus of the rock skeleton, and ρ be the density of the rock in saturated fluid. Indicates taking the real part; Step 3.3: Construct a forward model of seismic rock physics by combining the exact Zoeppritz equation: Where d is the pre-stack seismic data, W is the wavelet matrix, r is the reflection coefficient calculated by the exact Zoeppritz equation, and G is the forward process of the seismic data calculated by the exact Zoeppritz equation; Step 3.4: Calculate the reflection coefficient sequence r, the expression is as follows: The matrices A and b are expressed as follows: ， , where V P1 V S1 ρ1 represents the longitudinal wave velocity, transverse wave velocity, and density V of the medium in the upper half-space. P2 V S2 ρ2 represents the P-wave velocity, S-wave velocity, and density of the lower half-space medium, and i, i', j, and j' represent the P-wave incident angle, P-wave transmission angle, converted S-wave reflection angle, and converted S-wave transmission angle, respectively; Step 3.5: Based on convolution theory, using the wavelet matrix W and the reflection coefficient sequence r, calculate the synthetic seismic data d. syn : Where G(z) represents the forward modeling of seismic rock physics, and z is the parameter vector to be inverted, including the local porosity of hard pores. Microcrack local porosity Inclusion content and water saturation ； Forward modeling composite data, Given the wavelet matrix, calculate the actual P-wave dataset d. obs and forward simulation records residual ： 。 5. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 4, characterized in that, In step 3.5, d syn It has the following matrix form: In the formula, Angle of incidence The corresponding synthetic data is as follows: , where t 1…… t M Representing the tth 1…… t M Each depth-time sample point, from the angle of incidence. Wavelet matrices are constructed by extracting wavelets from actual seismic traces. ； Angle of incidence The reflection coefficient sequence obtained by time calculation is used to obtain: 。 6. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 5, characterized in that, In step 4, based on the parameter values ​​of the initial well logging model, the partial derivatives of the seismic rock physics forward modeling operator are calculated using the central difference method; Let G(z) be the derivative of the forward operator with respect to the model parameter z, and take the following form: ,in, Indicates the angle of incidence θ i The corresponding partial derivative matrix, where = 1, 2； Represents the Na-th incident angle; the reflection coefficient affects the local porosity of the hard pore. The partial derivatives are calculated by the following formula: ,in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. Partial derivatives, the specific formula is as follows: ,in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Local porosity of hard pores The partial derivatives, Local porosity of hard pores Tiny variables, Representing a time depth of t j time Value; Reflection coefficient for local porosity of microcracks The partial derivatives are calculated by the following formula: ,in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient is related to The partial derivatives are calculated using the central difference method. Partial derivatives, the specific formula is as follows: ,in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Local porosity of microcracks The partial derivatives, Local porosity of microcracks Tiny variables, Representing a time depth of t j time Value; Reflectance coefficient as a function of inclusion content V inc The partial derivatives are calculated by the following formula: ,in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to V inc The partial derivatives are calculated using the central difference method. Partial derivatives, the specific formula is as follows: ,in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Inclusion content V at the location inc The partial derivatives, V, the content of inclusions inc Tiny variables, Representing a time depth of t j V at time inc Value; Reflectance coefficient for water saturation S w The partial derivatives are calculated by the following formula: ,in, Indicates the angle of incidence as θ i At that time, in all time series, the reflection coefficient with respect to S w The partial derivatives are calculated using the central difference method. Partial derivatives, the specific formula is as follows: ,in, This represents the angle of incidence being θ. i The time depth is t j Time depth is t k Water saturation S at the location w The partial derivatives, Water saturation S w Tiny variables, Representing a time depth of t j S at time w value.

7. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 6, characterized in that, Step 5 includes the following sub-steps: Step 5.1, establish the inversion objective function L(z), the specific expression of which is: , where C z Let be the covariance matrix of the model parameters, u be the expectation of the model parameters, and λ be the regularization parameter; Step 5.2: Set the negative gradient as the update direction of the model. According to the LBFGS-MT algorithm, the iterative update formula for the model parameters is: ， , where v k v represents the update amount in the k-th iteration. k-1 This represents the update amount in the (k-1)th iteration. Z is the momentum coefficient. k Let z be the model parameter in the k-th iteration. k-1 Let κ be the model parameters in the (k-1)th iteration. k It is the step size of the k-th iteration; the specific expression of the first derivative J(z) of the objective function with respect to the model parameters z is as follows: Step 5.3: Let H(z) represent the pseudo-Hessian matrix of the objective function. Using the Gauss-Newton optimization algorithm, the pseudo-Hessian matrix takes the form: Where k is the number of iterations, ε k and S k All of these are intermediate variables, and their specific forms are as follows: ， ， Where I represents the identity matrix, y k It is also an intermediate variable, in the following form: , where J(z) k ) and J(z k+1 Let ) represent the Jacobian matrix values ​​for the k-th and (k+1)-th iterations, respectively; Step 5.4: Set the initial pseudo-Hessian matrix to the following form: The pseudo-Hessian matrix H in the (k+1)th iteration k+1 (z), in the following specific form: , where M is a constant.

8. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 7, characterized in that, In step 5.2, the covariance matrix C of the model parameters z The specific form is as follows: In the formula, σ is an M×M dimension matrix.

9. The seismic rock physics linear inversion method for microfracture porosity based on dual-dual-pore media according to claim 1, characterized in that, In step 7, the calculation of physical property parameters using the parameter inversion results from step 6 includes: the parameter inversion results from step 6 include the local porosity of hard pores. Local porosity of microcracks Inclusion content V inc and water saturation S w The final physical property parameters can be calculated using mathematical relationships. The specific mathematical expression is as follows: , ， ,in, For hard pore porosity, For microcrack porosity, Represents total porosity.

Images