Pre-stack seismic inversion method for biot coefficient and reservoir parameter based on rock physics model

Abstract

The present application belongs to the technical field of oil and gas reservoir geophysical exploration, and particularly relates to a Biot coefficient and reservoir parameter pre-stack seismic inversion method based on a rock physics model. First, a rock physics model capable of representing the influence of complex pore-fracture structure and fluid properties is constructed, and the model parameters are calibrated, so as to calculate the well continuous curve of the Biot coefficient and the fluid bulk modulus. Then, the PP wave reflection coefficient expression is derived, and a random nonlinear pre-stack seismic inversion method based on the Bayesian framework is constructed. The seismic waveform indication method is used to construct the seismic prior model of the parameters to be inverted, and the joint inversion is carried out by using the partial stack seismic data, so as to obtain the spatial distribution of the target parameters. The present application can comprehensively depict the reservoir characteristics from the aspects of reservoir effective stress sensitivity, pore development degree and pore fluid properties, so as to provide important geophysical basis for reservoir evaluation and oil and gas exploration and development.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas reservoir geophysical exploration technology, and particularly relates to a pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model. Background Technology

[0002] In the elastic theory of porous media, the Biot coefficient is used to describe the influence of pore fluid pressure on the stress response of the rock skeleton. It is an important parameter connecting total stress, effective stress, and pore pressure, and has significant application value in fracturing pressure prediction, hydraulic fracturing simulation, reservoir deformation analysis, wellbore stability evaluation, and formation pressure prediction. Therefore, accurately obtaining the Biot coefficient is of great significance for oil and gas reservoir evaluation and development.

[0003] Currently, the Biot coefficient is mainly obtained through laboratory measurements and calculations based on well logging data. Laboratory measurements are complex and costly, and the sample representativeness is limited. Calculations based on well logging data typically rely on empirical formulas or calculations using elastic modulus based on rock-physical relationships. However, these empirical formulas are often affected by lithological conditions and critical porosity assumptions, introducing uncertainty in multi-lithological or highly heterogeneous reservoirs. Furthermore, these methods are usually limited to the wellpoint scale and fail to reflect the spatial variation characteristics of the Biot coefficient at the reservoir regional scale. In practical applications, the Biot coefficient is sometimes simplified to a constant, which can lead to significant errors under complex reservoir conditions.

[0004] Seismic rock physics inversion is an important technique that combines seismic observation data with rock physics parameters. By introducing rock physics models, the relationship between seismic response and rock physics parameters can be established, enabling quantitative prediction of reservoir parameters. However, existing seismic inversion methods mostly focus on predicting elastic parameters or single reservoir parameters. There is a lack of research on the joint inversion of Biot coefficients, which reflect the coupling relationship between the rock skeleton and pore fluids, as well as related parameters characterizing pore development and fluid properties. Summary of the Invention

[0005] In view of this, the present invention aims to provide a pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model, so as to achieve joint seismic inversion of Biot coefficients, porosity and fluid bulk modulus, thereby providing a reliable geophysical basis for reservoir evaluation and oil and gas exploration and development.

[0006] To achieve the above objectives, the technical solution of the present invention is implemented as follows:

[0007] A pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model includes the following steps:

[0008] Step 1: Based on well logging mineral composition, pore-fracture structure and pore fluid properties, construct a rock physics model of the volcanic reservoir in the study area according to rock physics theory to describe the elastic response characteristics of the rock;

[0009] Step 2: Develop a well-based rock physics inversion method driven by the aforementioned rock physics model to determine well mineral composition and total rock porosity. Gas saturation As input parameters, the measured P-wave velocity and S-wave velocity in the well are fitted with the predicted P-wave velocity and S-wave velocity by the rock physics model to obtain the fracture porosity. and the aspect ratio of the crack The curve is used to calibrate the rock physics model;

[0010] Step 3: Based on the calibrated rock physics model, calculate the equivalent bulk modulus of the rock solid matrix. The bulk modulus of the skeleton of dry rocks And based on the definition of the Biot coefficient ,calculate The continuous curve is obtained simultaneously, and the fluid bulk modulus is acquired. The curve;

[0011] Step 4: Further calculate the combined parameters by combining well logging data. and ,in , , The density of the rock, The shear modulus of the rock is used, and the parameters to be inverted are constructed using the seismic waveform indication method. , , , and Earthquake prior models;

[0012] Step 5: Based on the simplified Gassmann fluid term, , , An approximate model of amplitude variation with offset distance is introduced to construct a theoretical expression for the reflection coefficient of PP waves;

[0013] Step 6: Introduce a stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework, and construct the pre-stack seismic inversion process based on the theoretical expression of PP wave reflection coefficient and the earthquake prior model;

[0014] Step 7: Under the constraints of the prior seismic model, input partially stacked seismic data, and use the pre-stack seismic inversion workflow constructed in Step 6 to perform the inversion. , , A joint inversion was performed to obtain the spatial distribution results of the three.

[0015] In step 1, when constructing the rock physics model, the mineral composition of the rock is determined based on the well logging data of the study area. The equivalent bulk modulus of the solid matrix and the shear modulus of the rock are calculated based on the Hashin-Shtrikman limit theory. The self-consistent approximation theory is used to introduce fluid-saturated pores and fractures to calculate the equivalent bulk modulus and shear modulus of the saturated rock, thereby obtaining the P-wave velocity and S-wave velocity predicted by the rock physics model.

[0016] In step 2, the objective function for well rock physics inversion is: The expression is: ;

[0017] In the formula, The P-wave velocity predicted by the rock physics model. The transverse wave velocity predicted by the rock physics model. The measured P-wave velocity in the well. This represents the transverse wave velocity measured in the well.

[0018] In step 3, based on the rock physics model calibrated in step 2, the equivalent bulk modulus of the rock solid matrix is ​​calculated using well logging mineral components and the Hashin-Shtrikman limit theory. Based on the calibrated rock physics model, the skeletal bulk modulus of dry rock was calculated by setting the fluid bulk modulus to 0. Simultaneously, the Wood equation was used to calculate the fluid bulk modulus. .

[0019] In step 5, the simplified Gassmann fluid term The expression is .

[0020] In step 5, the theoretical expression for the reflection coefficient of the PP wave is constructed based on the incident angle. As variables, by , , , and The linear weighting is used to construct the weighting coefficients, which are respectively: , , , , ;

[0021] Specifically as follows:

[0022] ;

[0023] in , , , ,as well as Representing the Biot coefficients on both sides of the reflective interface Total porosity of rocks Fluid bulk modulus ,parameter and parameters The difference.

[0024] Step 7 involves comparing and analyzing the inverted parameter profiles with the corresponding well logging curves to verify the rationality of the inversion results.

[0025] A pre-stack seismic inversion system for Biot coefficients and reservoir parameters based on a rock physics model, used to implement the aforementioned pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model, comprising:

[0026] The model building module is used to construct a rock physics model of the volcanic reservoir in the study area based on well logging mineral composition, pore-fracture structure and pore fluid properties, and rock physics theory.

[0027] The model calibration module is used to develop a well rock physics inversion method driven by the aforementioned rock physics model, for logging mineral composition and total rock porosity. Gas saturation Using known inputs, fracture porosity is obtained by fitting the measured P-wave and S-wave velocities in the well with the P-wave and S-wave velocities predicted by the rock physics model. and the aspect ratio of the crack The curve is used to calibrate the rock physics model;

[0028] The parameter calculation module is used to calculate the equivalent bulk modulus of the solid matrix of rock based on the calibrated rock physics model. The bulk modulus of the skeleton of dry rocks And based on the definition of the Biot coefficient ,calculate The continuous curve is obtained simultaneously, and the fluid bulk modulus is acquired. The curve;

[0029] The earthquake prior model construction module is used to calculate combined parameters by combining well logging data. and ,in , , The density of the rock, The shear modulus of the rock is used, and the parameters to be inverted are constructed using the seismic waveform indication method. , , , and Earthquake prior models;

[0030] A module for constructing theoretical expressions for the PP wave reflection coefficient, used to construct expressions based on simplified Gassmann fluid terms. , , An approximate model of amplitude variation with offset distance is introduced to construct a theoretical expression for the reflection coefficient of PP waves;

[0031] The pre-stack seismic inversion workflow establishment module is used to introduce the stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework and construct the pre-stack seismic inversion workflow based on the theoretical expression of PP wave reflection coefficient and the earthquake prior model.

[0032] The computation module is used to perform calculations on partially stacked seismic data, under the constraints of a prior seismic model, using a constructed pre-stack seismic inversion process, to... , , A joint inversion was performed to obtain the spatial distribution results of the three.

[0033] Compared with the prior art, the beneficial effects of the present invention are:

[0034] This invention proposes a pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model. By constructing a reservoir rock physics model for the study area and calibrating the model parameters using model-driven in-well rock physics inversion, the Biot coefficients can be calculated. and fluid bulk modulus The continuous curve in the well. Based on this, establish a system containing... , , The parameterized expression for the PP wave reflection coefficient is derived and introduced into a stochastic nonlinear pre-stack seismic inversion framework to achieve joint seismic inversion of the aforementioned parameters. This method combines rock physical relationships with pre-stack seismic data, creating clear physical constraints between the inversion parameters, thereby improving the stability and geological rationality of the inversion results. Simultaneously, by introducing a seismic waveform indication method to construct a priori models of the parameters to be inverted, the heterogeneity and spatial variation characteristics of the reservoir can be better characterized. Through... , , Quantitative prediction of spatial distribution can provide reliable geophysical basis for reservoir evaluation and oil and gas exploration and development. Attached Figure Description

[0035] Figure 1 This is a flowchart of the method of the present invention.

[0036] Figure 2 Flowchart for rock physics modeling.

[0037] Figure 3 This is a schematic diagram of the logging data and rock physics inversion results for well A1; where (a) represents the P-wave velocity. (a) includes the logging P-wave velocity and the model-predicted P-wave velocity; (b) is the shear wave velocity. (c) Porosity, including well logging shear wave velocity and model-predicted shear wave velocity; With crack porosity (d) represents the aspect ratio of the crack. (e) represents the bulk modulus. The equivalent bulk modulus of the rock solid matrix The bulk modulus of the skeleton of dry rocks ;(f) represents the Biot coefficient. (g) represents the fluid bulk modulus. ;(h) represents the rock mineral composition.

[0038] Figure 4 This is a schematic diagram of the logging data and rock physics inversion results for well A2; where (a) represents the P-wave velocity. (a) includes the logging P-wave velocity and the model-predicted P-wave velocity; (b) is the shear wave velocity. (c) Porosity, including well logging shear wave velocity and model-predicted shear wave velocity; With crack porosity (d) represents the aspect ratio of the crack. (e) represents the bulk modulus. Including the bulk modulus of the rock solid matrix The bulk modulus of the skeleton of dry rocks ;(f) represents the Biot coefficient. (g) represents the fluid bulk modulus. ;(h) represents the rock mineral composition.

[0039] Figure 5 This is a comparison chart of well logging data and pre-stack seismic inversion results for well A2 under noise-free conditions; where (a) represents the Biot coefficient. Comparison chart; (b) shows the total porosity of the rock. Comparison chart; (c) shows the fluid bulk modulus. Comparison chart; (d) represents parameters. Comparison chart; (e) represents parameters. Comparison charts; in (a)-(e), the blue curve represents well logging data, the black curve represents the prior model, and the red curve represents the inversion result; (f) shows the seismic record synthesized based on Equation 11 under noise-free conditions, where black represents the seismic record synthesized using the continuous curve in the well using the parameters to be inverted, and red represents the seismic record synthesized based on the inversion result; (g) shows the difference between the two types of synthesized seismic records mentioned above.

[0040] Figure 6 This is a comparison chart of well logging data and pre-stack seismic inversion results for well A2 under noisy conditions; where (a) represents the Biot coefficient. Comparison chart; (b) shows the total porosity of the rock. Comparison chart; (c) shows the fluid bulk modulus. Comparison chart; (d) represents parameters. Comparison chart; (e) represents parameters. Comparison charts; in (a)-(e), the blue curve represents well logging data, the black curve represents the prior model, and the red curve represents the inversion result; (f) shows the seismic record synthesized based on Equation 11 under noise-free conditions, where black represents the seismic record synthesized using the continuous curve in the well using the parameters to be inverted, and red represents the seismic record synthesized based on the inversion result; (g) shows the difference between the two types of synthesized seismic records mentioned above.

[0041] Figure 7 The images show the stacked profiles of seismic data from wells A1 and A2 at different incident angles. (a) is the stacked profile at an incident angle of 5°; (b) is the stacked profile at an incident angle of 10°; (c) is the stacked profile at an incident angle of 15°; (d) is the stacked profile at an incident angle of 20°; and (e) is a partial stacked profile at an incident angle of 25°.

[0042] Figure 8 This is a well-connected profile diagram of the prior model of the parameters to be inverted in wells A1 and A2, established based on the SWI method, where (a) is... (a) is the prior model profile; (b) is (c) is the prior model profile; The prior model profile; (d) is The prior model profile; (e) is A cross-sectional view of the prior model.

[0043] Figure 9 The diagram shows the well profile results obtained based on pre-stack seismic inversion, where (a) is... (a) is a cross-sectional view; (b) is (c) is a cross-sectional view; (d) is a cross-sectional view; (e) is a cross-sectional view; Cross-sectional view. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not constitute a limitation thereof.

[0045] In this embodiment of the invention, a pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model is proposed. First, a rock physics model of the volcanic reservoir in the study area is constructed to describe the influence of rock mineral composition, pore structure, and pore fluid properties on rock elastic parameters. Based on this, a model-driven in-well rock physics inversion method is developed. By fitting measured P-wave velocities and S-wave velocities, fracture-related parameters in the rock physics model are constrained, thereby calibrating the rock physics model. Subsequently, the equivalent bulk modulus of the rock solid matrix is ​​calculated based on the calibrated rock physics model. and the bulk modulus of the skeleton of dry rocks And calculate the Biot coefficient according to the definition of the Biot coefficient. The continuous curve in the well was obtained, and the fluid bulk modulus was calculated in conjunction with the fluid properties. Based on this, establish a system that includes , , The parameterized expression for the PP wave reflection coefficient was derived and introduced into the stochastic nonlinear pre-stack seismic inversion framework. Furthermore, combining well logging data and seismic waveform information, the target parameters were constructed using a seismic waveform indication method. , , and parameters and A priori seismic model was developed. Under the constraints of the priori model, pre-stack seismic inversion was performed using partially angle-stacked seismic data to obtain the seismic data for the study area. , , The spatial distribution of the reservoir parameters is used to achieve quantitative prediction of the Biot coefficient and reservoir parameters in seismic data.

[0046] The invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0047] Example:

[0048] A pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model is described below. Figure 1 As shown, the method includes the following steps:

[0049] Step 1: Taking into account the logging mineral composition, pore-fracture structure and pore fluid properties, a rock physics model of the volcanic reservoir in the study area is constructed based on rock physics theory to describe the elastic response characteristics of the rock.

[0050] Step 2: Develop a well-based rock physics inversion method driven by the aforementioned rock physics model to determine well mineral composition and total rock porosity. Gas saturation Inverting fracture porosity as a known input and the aspect ratio of the crack Fracture porosity is obtained by fitting the measured P-wave and S-wave velocities in the well with the P-wave and S-wave velocities predicted by the rock physics model. and the aspect ratio of the crack The parameter curves are used to calibrate the rock physics model;

[0051] Step 3: Based on the calibrated rock physics model, calculate the equivalent bulk modulus of the rock solid matrix. The bulk modulus of the skeleton of dry rocks And based on the definition of the Biot coefficient Calculation well The continuous curve is obtained simultaneously, and the fluid bulk modulus is acquired. The curve;

[0052] Step 4: Based on the parameters obtained in Step 3, further calculate the combined parameters by combining them with well logging data. and ,in , , The density of the rock, The shear modulus of the rock was determined, and the parameters to be inverted were constructed using the Seismic Waveform Indication (SWI) method. , , , and Earthquake prior models;

[0053] Step 5: Based on the simplified Gassmann fluid term, , , An approximate model of amplitude variation with offset distance (AVO) is introduced to construct a theoretical expression for the reflection coefficient of PP waves;

[0054] Step 6: Introduce the stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework, and construct the pre-stack seismic inversion process based on the theoretical expression of PP wave reflection coefficient established in Step 5 and the seismic prior model established in Step 4.

[0055] Step 7: Under the constraints of the prior seismic model, input partially stacked seismic data, and use the pre-stack seismic inversion workflow constructed in Step 6 to perform the inversion. , , A joint inversion was performed to obtain the spatial distribution results of the three.

[0056] Furthermore, in step 1, when constructing the rock physics model, the mineral composition of the rock is determined based on the well logging data of the study area. The equivalent bulk modulus of the solid matrix and the shear modulus of the rock are calculated based on the Hashin-Shtrikman limit theory. The self-consistent approximation (SCA) theory is used to introduce fluid-saturated pores and fractures to calculate the equivalent bulk modulus and shear modulus of the saturated rock, thereby obtaining the P-wave velocity and S-wave velocity predicted by the rock physics model.

[0057] Furthermore, in step 2, the objective function for well rock physics inversion is: The expression is as follows:

[0058] ;

[0059] In the formula, The P-wave velocity predicted by the rock physics model. The transverse wave velocity predicted by the rock physics model. The measured P-wave velocity in the well. This represents the transverse wave velocity measured in the well.

[0060] Continuously adjust the unknown and The curve, each time it is adjusted, will have a new... , and fixed known inputs (mineral composition, , The data is fed into the model together to calculate a new predicted velocity curve. This process continues until the predicted velocity curve and the measured velocity curve fit each other well enough. , The curve is the result of the inversion; the curve obtained from the inversion... , The curve, as a known parameter, is substituted back into the original rock physics model to obtain the calibrated rock physics model.

[0061] Furthermore, in step 3, based on the rock physics model calibrated in step 2, the equivalent bulk modulus of the rock solid matrix is ​​calculated using the well logging mineral composition curves and the Hashin-Shtrikman limit theory. After inputting the fracture porosity and fracture aspect ratio obtained in step 2 into the rock physics model, the skeletal bulk modulus of the dry rock is calculated by setting the fluid bulk modulus of the pores and fractures to 0. Simultaneously, the Wood equation was used to calculate the fluid bulk modulus. .

[0062] Furthermore, in step 4, the Seismic Waveform Indication (SWI) method can construct a priori model that conforms to the geological characteristics of complex reservoirs by integrating well logging data and seismic waveform information.

[0063] Furthermore, in step 5, the simplified Gassmann fluid term The expression is .

[0064] Furthermore, in step 5, the theoretical expression for the PP wave reflection coefficient is constructed. With the angle of incidence As variables, by , , , and The linear weighting is used to construct the weighting coefficients, which are respectively: , , , , All are related to the angle of incidence. Correlation coefficient;

[0065] Specifically as follows:

[0066] ;

[0067] in , , , ,as well as Representing the Biot coefficients on both sides of the reflective interface Total porosity of rocks Fluid bulk modulus ,parameter and parameters The difference.

[0068] Furthermore, in step 7, the rationality of the inversion results is verified by comparing and analyzing the parameter profiles obtained from the inversion with the corresponding logging curves in the well.

[0069] The specific implementation process is as follows:

[0070] Based on the mineral composition of well logging data from the study area, the mineral composition of the rocks was determined, including felsic materials, clay, potassium feldspar, pyrite, and carbonates. Subsequently, the equivalent bulk modulus of the rock matrix was calculated using the Hashin-Shtrikman limit theory. and shear modulus The specific calculation formula is as follows:

[0071] Formula 1: ; ;

[0072] ; ;

[0073] In the formula, These represent the upper and lower limits of the effective bulk modulus of the rock solid matrix, respectively. These represent the upper and lower limits of the shear modulus of the solid rock matrix, respectively. for The average of the upper and lower limits. for The average of the upper and lower limits; and These represent the maximum and minimum bulk moduli of each mineral in the solid matrix of the rock, respectively. and These represent the maximum and minimum shear moduli of each mineral in the solid matrix of the rock, respectively. , , For operators, the calculation method is as follows:

[0074] Formula 2: ;

[0075] ;

[0076] ;

[0077] In the formula, Refers to the variable within the parentheses corresponding to the operator; Referential variable or ; Referential variable or ; and Representing the first solid matrix of the rock respectively The bulk modulus and shear modulus of each component; This indicates a weighted average calculation.

[0078] Volcanic reservoir rocks typically exhibit a porous-fracture structure. To accurately characterize their elastic response, the Self-Compatible Approximation (SCA) theory is used to introduce saturated fluid into the porous-fracture structure of the rock's solid matrix. The specific calculation formula is as follows:

[0079] Formula 3: ;

[0080] ;

[0081] In the formula, subscripts Indicates the first in the rock Inclusions; In rocks Volume fraction of the included substances; Indicates the first in the rock The bulk modulus of the inclusions; Indicates the first in the rock The shear modulus of the inclusions; and These are the equivalent bulk modulus and shear modulus of saturated rock, respectively. and Each of the components (the first) For details on the calculation method of geometric factors related to the inclusions, please refer to Mavko et al. (2009). This refers to the amount of matter contained within the rock.

[0082] Fluid bulk modulus in pores and cracks The calculation is performed using the Wood equation, as shown in the following formula:

[0083] Formula 4: ;

[0084] In the formula, and These are the bulk moduli of gas and water, respectively. This represents the water saturation level of the rock pores.

[0085] Based on the equivalent bulk modulus and shear modulus of saturated rock and Calculate the P-wave velocity of the rock and transverse wave velocity The formula is as follows:

[0086] Formula 5: ;

[0087] Formula 6: ;

[0088] in This represents the density of the rock.

[0089] A rock physics model of the volcanic reservoir in the study area was constructed using the aforementioned rock physics theories, and a model-driven method for inverting well rock physics parameters was developed. Mineral composition, total porosity, and gas saturation logging curves were input into the rock physics model. By fitting the measured P-wave and S-wave velocities in the well with the P-wave and S-wave velocities predicted by the rock physics model, the fracture parameters (fracture porosity) in the well were quantitatively calculated. The aspect ratio of the crack The curve. The objective function for inverting the rock physical parameters in the well. for:

[0090] Formula 7: ;

[0091] By solving the above objective function, the optimal crack porosity is obtained. Compared with the aspect ratio of cracks At this point, the model's predicted velocity and the measured velocity are optimally matched, thus achieving the calibration of the rock physics model and ensuring the equivalent bulk modulus of the rock solid matrix required for subsequent steps. , bulk modulus of rock skeleton Fluid bulk modulus within pores and Biot coefficient It is true and reliable.

[0092] Based on the above inversion results, a calibrated rock physics model applicable to the reservoirs in the study area was obtained.

[0093] Based on this, recalculate and According to the definition of the Biot coefficient ,calculate Continuous curve in the well; fluid bulk modulus The curve is calculated using Equation 4; combined with well logging data, further information including target parameters can be obtained. , , and combined parameters and The wellbore continuous curves of all parameters to be inverted, including , The parameter curves provide reliable input for the construction of the prior model and lay the foundation for subsequent pre-stack seismic inversion.

[0094] Russell et al. derived the equation including the Gassmann fluid term. shear modulus and density The approximate expression for AVO is shown in Equation 8 below:

[0095] ;

[0096] in, , , These represent the differences in the Gassmann fluid term, shear modulus, and density between the media above and below the reflecting interface, respectively. , , Representing the Gassmann fluid terms respectively shear modulus and density The average value of the upper and lower media on both sides of the interface; The angle of incidence of the earthquake; and These represent the ratio of longitudinal wave velocity to transverse wave velocity for saturated rock and dry rock, respectively.

[0097] Gassmann fluid terms The expression is as follows:

[0098] Formula 9: ;

[0099] in, This refers to the total porosity of the rock, which is the same as the porosity mentioned above. The equivalent bulk modulus of the solid rock matrix. The bulk modulus of the dry rock skeleton. This represents the bulk modulus of the fluid within the pores. Based on existing rock physics theories, the Gassmann fluid term... The approximate expression is:

[0100] Formula 10: ;

[0101] Applying Equation 10 to Equation 8, a new theoretical expression for the reflection coefficient of PP waves is constructed, with the incident angle... As variables, by , , , and The linear weighting is used to construct the weighting coefficients, which are respectively: , , , , All are related to the angle of incidence. The relevant coefficients are as follows:

[0102] Formula 11: ;

[0103] in, and The combined parameters introduced in the parameterized PP wave reflection coefficient formula are defined as follows: , ; , , , , Representing the Biot coefficients on both sides of the reflective interface Total porosity of rocks Fluid bulk modulus ,parameter and parameters The difference.

[0104] coefficient , , , , The expression is as follows:

[0105] Formula 12:

[0106] ;

[0107] ;

[0108] ;

[0109] ;

[0110] ;

[0111] Based on the derived theoretical formula for reflection coefficient, a stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework is introduced to construct a pre-stack seismic inversion workflow and invert the target parameters. , and .

[0112] Equation 13 below can be reformulated as a linear combination of the relative changes of the parameters to be inverted:

[0113] ;

[0114] in, , , , and , respectively represent , , , and The relative change. Assume the number of sampling points for each seismic trace is... Multiple incident angles ( Reflection coefficient sequence and seismic wavelet matrix under the condition of incident angle (number of incident angles) By performing convolution, the forward model of the seismic response is constructed as follows:

[0115] Formula 14: ;

[0116] In the formula, Indicates the first The observed seismic records corresponding to each incident angle are: Column vector; , , , , All diagonal weighted matrix, , , , , Each represents a corresponding parameter ( , , , and The reflection coefficient of PP waves constitutes vector.

[0117] The earthquake forward model shown in Equation 14 can be simplified as follows:

[0118] Formula 15: ;

[0119] in:

[0120] Formula 16: ;

[0121] Formula 17: ;

[0122] Formula 18: ;

[0123] Among them, superscript This represents the transpose of a matrix.

[0124] Building upon this, a nonlinear stochastic pre-stack seismic inversion method based on Bayesian theory is introduced, utilizing multi-incident angle seismic observation data. Stable estimation model parameter vector According to Bayes' theorem, the posterior probability distribution of the model parameters... From the prior distribution of model parameters Likelihood function of seismic data They are jointly determined, and their relationship can be expressed as the product of the two:

[0125] Formula 19:

[0126] In the formula, This is the normalization factor.

[0127] For matrices with orthogonal operators For a linear seismic forward model, assuming that the prior distribution of model parameters and the distribution of observation errors both follow a Gaussian distribution, the posterior probability density function of the model parameters is expressed as: Equation 20:

[0128] ;

[0129] in, The normalization constant is and Let these represent the prior distribution of model parameters and the covariance matrix of seismic data, respectively. This is the prior mean. The objective function is defined by the model constraint terms. Data mismatch Composition, expressed as follows:

[0130] Equation 21: ;

[0131] in,

[0132] Equation 22: ;

[0133] Equation 23: ;

[0134] This invention employs a stochastic inversion framework based on an ensemble smoother with multiple data assimilation (ES-MDA). By minimizing the objective function, it obtains the maximum a posteriori (MAP) solution for the model parameters, which corresponds to the posterior mean. The expression is as follows:

[0135] Formula 24: ;

[0136] in:

[0137] Formula 25: ;

[0138] in, For the forward calculus operator matrix, the corresponding posterior covariance matrix of the model parameters is... It can be represented as:

[0139] Equation 26: ;

[0140] in, The analytical Kalman gain matrix is ​​expressed as follows:

[0141] Equation 27: ;

[0142] In the multiple data assimilation set smoothing method, each set member is updated in the following way:

[0143] Equation 28: , ;

[0144] in, Represents the size of the set, superscript and The ensemble gain matrix corresponds to the parameter update state and the prior state, respectively. Determine as follows:

[0145] Equation 29: , ;

[0146] in, Represents prior model parameters With corresponding synthetic data The covariance matrix between them Represents the covariance matrix of the synthesized data; This represents an observation data vector containing random perturbations.

[0147] Through the above stochastic nonlinear inversion process based on Bayesian theory, the model parameters are solved. Quantitative inversion can be performed from pre-stack seismic data. , , Considering the impact of prior models on the stability and reliability of pre-stack seismic inversion results, the Seismic Waveform Indication (SWI) method is introduced to construct a prior model that conforms to the reservoir geological characteristics of the study area, providing credible and physically consistent prior constraints for the proposed pre-stack seismic inversion process.

[0148] Based on a systematic analysis of the reservoir geological background, core data, and well logging data in the study area, and considering the complex mineral composition of the reservoir, the development characteristics of pores and fractures in the rocks, and the complexity of pore fluid properties, a petrophysical model suitable for the volcanic reservoirs in the study area is constructed to reasonably characterize their elastic response features. The modeling process is as follows: Figure 2 As shown.

[0149] Based on the established rock physics model, a model-driven method for inverting well rock physics parameters was developed. This method was applied to well logging data in the study area to calibrate the established rock physics model. The inversion results for wells A1 and A2 are shown below. Figure 3 and Figure 4 As shown in the figure. The results indicate that the measured P-wave velocity and S-wave velocity in the well show good consistency with the model predictions in both overall trend and local details.

[0150] In well A1, the correlation coefficients between the measured and predicted P-wave and S-wave velocities were 0.8965 and 0.9712, respectively, with corresponding root mean square errors of 0.1391 and 0.0469. In well A2, the correlation coefficients were 0.8442 and 0.9130, respectively, with corresponding root mean square errors of 0.1710 and 0.0809. These results demonstrate that the constructed rock physics model can effectively characterize the elastic response features of the reservoir rocks in the study area, verifying the rationality and effectiveness of the established rock physics model and the developed model-driven well rock physics parameter inversion method.

[0151] Based on the above inversion results, the equivalent bulk modulus of the rock matrix was further calculated. With skeleton bulk modulus And based on the Biot coefficient The definition formula was used to calculate the well continuity curve; simultaneously, the pore fluid bulk modulus was calculated according to Equation 4. This provides a reliable parameter basis for the construction of prior models in the subsequent earthquake inversion process.

[0152] To evaluate the accuracy and robustness of the theoretical expression for the PP wave reflection coefficient constructed in this invention and the pre-stack seismic inversion process developed based on it in terms of quantitative inversion target parameters, the corresponding PP wave reflection coefficient sequence was calculated using Equation 11 based on the logging data of well A2, and the reflection coefficient was convolved with the Ricker wavelet with a dominant frequency of 25Hz to generate a synthetic angle gather for synthetic data inversion testing.

[0153] Synthetic seismic records generated under noise-free conditions and their corresponding inversion results are as follows: Figure 5 As shown in the figure, the results indicate that the inversion results (red) and well logging parameter curves (blue) show high consistency in both overall trends and local details, demonstrating that the method of this invention can accurately recover the target parameters under noise-free conditions. Furthermore, comparing the seismic record synthesized based on well logging parameters (black) with the seismic record reconstructed based on the inversion results (red), the two show good agreement, with smaller residual gather amplitudes, further validating the accuracy and reliability of the proposed inversion framework.

[0154] To further verify the robustness of the method under actual noise conditions, Gaussian random noise with a signal-to-noise ratio of 5 was added to the synthetic seismic data, and the inversion experiment was repeated under the same inversion procedure. The results are as follows. Figure 6As shown, although the introduction of noise affects the inversion accuracy to some extent, the overall inversion results of each parameter still maintain good consistency with the well logging parameters; at the same time, the pre-stack seismic record reconstructed from the inversion parameters still has high consistency with the noisy synthetic seismic data, and its residuals are within an acceptable range. The above results indicate that the method of the present invention still has good stability and noise resistance under the condition of noise.

[0155] Figure 7 Partial angle-stacked seismic profiles under different incident angle conditions are presented as input data for pre-stack seismic inversion.

[0156] A priori model of the parameters to be inverted was constructed using the seismic waveform indication simulation method. By fusing post-stack seismic data with well logging parameters from wells A1 and A2, a priori model conforming to the reservoir geological characteristics of the study area was established and applied to the pre-stack seismic inversion process. Figure 8 As shown, the constructed prior model can well characterize the low-frequency variation characteristics of logging parameters at the well location, indicating that the SWI method can provide stable and reasonable prior constraints under the constraint of seismic waveform information.

[0157] Under the constraints of the prior model, the pre-stack seismic inversion process proposed in this invention is used to jointly invert multi-incident angle partially stacked seismic data to obtain the spatial distribution of the target parameters within the study area, such as... Figure 9 As shown in the figure. By comparing and verifying the in-well logging parameter curves with the inversion results, it is demonstrated that the method of the present invention has high stability and reliability in the inversion of parameters in complex reservoirs.

[0158] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0159] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model, characterized in that, Includes the following steps: Step 1: Based on well logging mineral composition, pore-fracture structure and pore fluid properties, construct a rock physics model of the volcanic reservoir in the study area according to rock physics theory to describe the elastic response characteristics of the rock; Step 2: Develop a well-based rock physics inversion method driven by the aforementioned rock physics model to determine well mineral composition and total rock porosity. Gas saturation As input parameters, the measured P-wave velocity and S-wave velocity in the well are fitted with the predicted P-wave velocity and S-wave velocity by the rock physics model to obtain the fracture porosity. and the aspect ratio of the crack The curve is used to calibrate the rock physics model; Step 3: Based on the calibrated rock physics model, calculate the equivalent bulk modulus of the rock solid matrix. The bulk modulus of the skeleton of dry rocks And based on the definition of the Biot coefficient ,calculate The continuous curve is obtained simultaneously, and the fluid bulk modulus is acquired. The curve; Step 4: Further calculate the combined parameters by combining well logging data. and ,in , , The density of the rock, The shear modulus of the rock is used, and the parameters to be inverted are constructed using the seismic waveform indication method. , , , and Earthquake prior models; Step 5: Based on the simplified Gassmann fluid term, , , An approximate model of amplitude variation with offset distance is introduced to construct a theoretical expression for the reflection coefficient of PP waves; Step 6: Introduce a stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework, and construct the pre-stack seismic inversion process based on the theoretical expression of PP wave reflection coefficient and the earthquake prior model; Step 7: Under the constraints of the prior seismic model, input partially stacked seismic data, and use the pre-stack seismic inversion workflow constructed in Step 6 to perform the inversion. , , A joint inversion was performed to obtain the spatial distribution results of the three.

2. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, In step 1, when constructing the rock physics model, the mineral composition of the rock is determined based on the well logging data of the study area. The equivalent bulk modulus of the solid matrix and the shear modulus of the rock are calculated based on the Hashin-Shtrikman limit theory. The self-consistent approximation theory is used to introduce fluid-saturated pores and fractures to calculate the equivalent bulk modulus and shear modulus of the saturated rock, thereby obtaining the P-wave velocity and S-wave velocity predicted by the rock physics model.

3. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, In step 2, the objective function for well rock physics inversion is: The expression is: In the formula, The P-wave velocity predicted by the rock physics model. The transverse wave velocity predicted by the rock physics model. The measured P-wave velocity in the well. This represents the transverse wave velocity measured in the well.

4. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, In step 3, based on the rock physics model calibrated in step 2, the equivalent bulk modulus of the rock solid matrix is ​​calculated using well logging mineral components and the Hashin-Shtrikman limit theory. Based on the calibrated rock physics model, the skeletal bulk modulus of dry rock was calculated by setting the fluid bulk modulus to 0. Simultaneously, the Wood equation was used to calculate the fluid bulk modulus. .

5. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, In step 5, the simplified Gassmann fluid term The expression is .

6. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, In step 5, the theoretical expression for the reflection coefficient of the PP wave is constructed based on the incident angle. As variables, by , , , and The linear weighting is used to construct the weighting coefficients, which are respectively: , , , , The details are as follows: ； in , , , as well as Representing the Biot coefficients on both sides of the reflective interface Total porosity of rocks Fluid bulk modulus ,parameter and parameters The difference.

7. The pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model according to claim 1, characterized in that, Step 7 involves comparing and analyzing the inverted parameter profiles with the corresponding well logging curves to verify the rationality of the inversion results.

8. A pre-stack seismic inversion system for Biot coefficients and reservoir parameters based on a rock physics model, used to implement the pre-stack seismic inversion method for Biot coefficients and reservoir parameters based on a rock physics model as described in any one of claims 1-7, characterized in that, include: The model building module is used to construct a rock physics model of the volcanic reservoir in the study area based on well logging mineral composition, pore-fracture structure and pore fluid properties, and rock physics theory. The model calibration module is used to develop a well rock physics inversion method driven by the aforementioned rock physics model, for logging mineral composition and total rock porosity. Gas saturation Using known inputs, fracture porosity is obtained by fitting the measured P-wave and S-wave velocities in the well with the P-wave and S-wave velocities predicted by the rock physics model. and the aspect ratio of the crack The curve is used to calibrate the rock physics model; The parameter calculation module is used to calculate the equivalent bulk modulus of the solid matrix of rock based on the calibrated rock physics model. The bulk modulus of the skeleton of dry rocks And based on the definition of the Biot coefficient ,calculate The continuous curve is obtained simultaneously, and the fluid bulk modulus is acquired. The curve; The earthquake prior model construction module is used to calculate combined parameters by combining well logging data. and ,in , , The density of the rock, The shear modulus of the rock is used, and the parameters to be inverted are constructed using the seismic waveform indication method. , , , and Earthquake prior models; A module for constructing theoretical expressions for the PP wave reflection coefficient, used to construct expressions based on simplified Gassmann fluid terms. , , An approximate model of amplitude variation with offset distance is introduced to construct a theoretical expression for the reflection coefficient of PP waves; The pre-stack seismic inversion workflow establishment module is used to introduce the stochastic nonlinear pre-stack seismic inversion method under the Bayesian framework and construct the pre-stack seismic inversion workflow based on the theoretical expression of PP wave reflection coefficient and the earthquake prior model. The computation module is used to perform calculations on partially stacked seismic data, under the constraints of a prior seismic model, using a constructed pre-stack seismic inversion process, to... , , A joint inversion was performed to obtain the spatial distribution results of the three.

Images