Prestack brittleness index inversion method and system based on Zoeppritz equation

Abstract

The invention provides a pre-stack brittleness index inversion method based on a Zoeppritz equation. The inversion method comprises the steps of S1, constructing a relationship between a brittleness index and a seismic record; s2, constructing an objective function about a brittleness index based on a least square principle; s3, establishing an iterative inversion equation for deducing the brittleness index of the reservoir through a Gaussian-Newton iteration method in combination with the objective function; and S4, giving an initial model and the iterative inversion equation, and carrying out iterative operation to obtain a final reservoir brittleness index. The ill-conditioned inversion problem is solved through various constraints, and the stability and accuracy of brittleness index estimation are improved.

Description

Technical Field

[0001] The invention relates to the field of oil and gas geophysical exploration, and in particular to a prestack brittleness index inversion method and system based on the Zoeppritz equation. Background Art

[0002] Conventional prestack brittleness index inversion is usually based on the inversion of the Zoeppritz equation to obtain the longitudinal and transverse wave velocities or Young's modulus and Poisson's ratio, and then indirectly solves the brittleness index through methods such as the Rickman brittleness formula. However, the conventional prestack brittleness index inversion is divided into two or even three steps, which will introduce certain cumulative errors and cause inaccurate results.

[0003] There are problems such as cumulative errors in the brittleness index inversion method used in the prior art. Summary of the invention

[0004] In view of the above problems, the present invention is proposed to provide a prestack brittleness index inversion method and system based on the Zoeppritz equation that overcomes the above problems or at least partially solves the above problems.

[0005] According to one aspect of the present invention, a prestack brittleness index inversion method based on the Zoeppritz equation is provided, the inversion method comprising:

[0006] Step S1, constructing the relationship between the brittleness index and the earthquake record;

[0007] Step S2, constructing an objective function of the brittleness index based on the least squares principle;

[0008] Step S3, establishing an iterative inversion equation for deriving the reservoir brittleness index using a Gauss-Newton iteration method in combination with the objective function;

[0009] Step S4: Given the initial model and the iterative inversion equation, an iterative operation is performed to obtain a final reservoir brittleness index.

[0010] Optionally, the step of constructing a relationship between the brittleness index and the seismic record specifically includes:

[0011] Based on the intersection analysis of brittleness index and elastic parameters, the relationship between brittleness index and elastic parameters is obtained, and then based on the Zoeppritz equation, the relationship between brittleness index and seismic records is constructed.

[0012] Optionally, the constraint items of the objective function specifically include: a Cauchy distribution constraint item and a low-frequency model constraint item.

[0013] Optionally, the relationship between the brittleness index and the elastic parameter obtained by the intersection analysis of the brittleness index and the elastic parameter is ρ=φ(BI,V P,V S ).

[0015] Optionally, constructing the relationship between the brittleness index and the seismic record based on the Zoeppritz equation specifically includes:

[0016] Introducing the relationship between the brittleness index and the elastic parameter into the Zoeppritz equation; obtaining the mathematical relationship between the longitudinal wave reflection coefficient of the interface and the longitudinal and transverse wave velocities and brittleness index on both sides of the interface at a certain incident angle;

[0017]

[0018]

[0019]

[0020] Where p = sinθ1 / V P1 is the ray parameter, θ1 is the incident angle, BI, V P 、V S , ρ represent the brittleness index, longitudinal wave velocity, shear wave velocity, and density of the medium, respectively. Subscript 1 represents the upper medium, and subscript 2 represents the lower medium.

[0021] Optionally, the step of constructing a relationship between the brittleness index and the seismic record based on the Zoeppritz equation further includes:

[0022] When there are n reflection interfaces, the P-wave reflection coefficient is expressed in the form of reservoir brittleness index and elastic parameters: P =f(θ,BI1,...,BI n+1 ,V P1 ,...,V Pn+1 ,V S1 ,...,V Sn+1 );

[0023] The form of the vector function is r(θ)=f(θ,m);

[0024] in f is a forward operator, which represents the nonlinear function of the reflection coefficient sequence r with respect to the reservoir brittleness index and elastic parameter m;

[0025] Considering that the seismic record is the convolution of the wavelet and the reflection coefficient, and considering that each gather has multiple records corresponding to different incident angles, the equation d=g(m) for synthesizing the seismic angle gather record using the reservoir brittleness index and elastic parameters is obtained;

[0026] Wherein, d represents the seismic angle gather record, g is the nonlinear forward operator, and represents the function of the angle gather record d with respect to the vector m composed of the reservoir brittleness index and elastic parameters.

[0027] Optionally, the constraint term of the objective function is expressed as:

[0028]

[0029] in, It represents the objective function of the reservoir brittleness index m, which consists of three parts, which represent the data error term, the Cauchy prior constraint term, and the low-frequency model constraint term from left to right;

[0030] In the data error term, the superscript T is the transposition symbol; d represents the seismic observation record; m represents the reservoir brittleness index; g is the nonlinear forward operator, which represents the relationship between the seismic record and the reservoir brittleness index; W d represents the weighting matrix;

[0031] In the Cauchy prior constraint, is the weight coefficient; m c represents the relative change of the reservoir brittleness index m, which conforms to the three-variable Cauchy distribution; Φ i =(D i ) T Ψ -1 D i ; Ψ is a 3×3 scale matrix; D i is a 3×3n matrix; i represents the i-th seismic record sampling point, 1≤i≤n, and n represents the number of seismic record sampling points;

[0032] In the low-frequency model constraint term, ε m 2 Indicates the weight coefficient of this item; m r represents the low-frequency model; W m represents the weighting matrix.

[0033] Optionally, the step S3, establishing an iterative inversion equation for deriving the reservoir brittleness index by using the Gauss-Newton iteration method in combination with the objective function, specifically comprises:

[0034] Step S31: Calculate according to the objective function expression The first-order derivative of m gives the gradient of the objective function

[0035] Step S32: Calculate according to the objective function gradient expression The first-order derivative of m, ignoring the higher-order derivatives in the result, obtains the Hessian matrix of the objective function

[0036] Step S33: Assume that the initial model is m(0) , respectively substitute it into the gradient of the objective function and the Hessian matrix expression to calculate the objective function in m (0) The gradient at and the Hessian matrix The iterative equation of reservoir brittleness index m is obtained.

[0037] Optionally, the gradient of the objective function The expression is:

[0038]

[0039] in

[0040] Optionally, the Hessian matrix of the objective function The expression is

[0041]

[0042] Optionally, the iterative equation of the reservoir brittleness index m is:

[0043]

[0044] Optionally, the step S4, performing iterative calculation on a given initial model and the iterative inversion equation to obtain a final reservoir brittleness index specifically includes:

[0045] Given an initial model m (0) , substitute into the iterative equation of reservoir brittleness index to calculate m;

[0046] Determine the calculated result m and the initial model m (0) Whether the difference between is less than the set threshold or reaches the set maximum number of iterations, if so, m is the final result, otherwise m is set as the initial model for cyclic iteration.

[0047] The present invention also provides a prestack brittleness index inversion system based on the Zoeppritz equation, which is applied to the prestack brittleness index inversion method based on the Zoeppritz equation described above. The inversion system comprises:

[0048] Obtain parameter module to build the relationship between brittleness index and earthquake records;

[0049] An objective function building module, used to build an objective function about the brittleness index based on the least squares principle;

[0050] An inversion method module, used for establishing an iterative inversion equation for deriving a reservoir brittleness index using a Gauss-Newton iteration method in combination with the objective function;

[0051] The output module is used to perform iterative calculations based on the initial model and the iterative inversion equation to obtain the final reservoir brittleness index.

[0052] The present invention provides a prestack brittleness index inversion method and system based on the Zoeppritz equation, the inversion method comprising: step S1, constructing the relationship between the brittleness index and the seismic record; step S2, constructing the objective function of the brittleness index based on the least squares principle; step S3, establishing the iterative inversion equation for deriving the reservoir brittleness index by the Gauss-Newton iteration method in combination with the objective function; step S4, performing iterative calculations given the initial model and the iterative inversion equation to obtain the final reservoir brittleness index. The ill-conditioned inversion problem is solved by multiple constraints to improve the stability and accuracy of the brittleness index estimation.

[0053] The above description is only an overview of the technical solution of the present invention. In order to more clearly understand the technical means of the present invention, it can be implemented according to the contents of the specification. In order to make the above and other purposes, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are listed below. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0055] Figure 1 A flowchart of a prestack brittleness index inversion method based on the Zoeppritz equation provided in an embodiment of the present invention;

[0056] Figure 2 P-wave velocity profile of the layered model provided by the embodiment of the present invention

[0057] Figure 3 The inversion result section of the layered model data provided by the embodiment of the present invention

[0058] Figure 4 An actual well data and inversion result of a certain work area provided by the embodiment of the present invention

[0059] Figure 5 A structural diagram of a prestack brittleness index inversion system based on the Zoeppritz equation provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0060] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0061] The terms "comprises" and "having" and any variations thereof in the description embodiments, claims and drawings of the present invention are intended to cover non-exclusive inclusions, for example, including a series of steps or units.

[0062] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments.

[0063] like Figure 1 As shown, the present invention provides a direct inversion method for prestack brittleness index based on the Zoeppritz equation, and the method specifically comprises the following steps:

[0064] Step 1: Based on the Zoeppritz equation, the relationship between the brittleness index and the elastic parameters, the relationship between the brittleness index and the earthquake record is directly constructed;

[0065] Step 2: Construct the objective function of the brittleness index based on the least squares principle, and add the Cauchy distribution constraint term and the low-frequency model constraint term;

[0066] Step 3: Use the Gauss-Newton iteration method to derive the iterative inversion equation of the reservoir brittleness index;

[0067] Step 4: Given the initial model, perform iterative calculations. If the difference between the calculated result and the initial model is small enough or reaches the set maximum number of iterations, stop the calculation and obtain the final reservoir brittleness index.

[0068] Embodiment 1:

[0069] First, a reservoir brittleness index theoretical model section is designed, and then the reservoir brittleness index is forward modeled to synthesize the prestack angle gather record. Finally, the synthetic record is used to invert the reservoir brittleness index through the method of the present invention, and the inversion result is compared with the theoretical model, such as Figure 2 As shown, Figure 2 a is the accurate brittleness index model, Figure 2 b is the inversion result, Figure 2 c is the error statistics of the inversion, and the average relative error of the brittleness index inversion is 0.95%. Take the single-channel data at the position of 500m, and the inversion results are as follows: Figure 3 As shown in a, the solid line represents the true result and the dotted line represents the inversion result. Figure 3b is the brittleness index inversion error curve, and its average relative error is 0.8%; the comparison between the seismic records synthesized by the inversion result and the original seismic records is as follows: Figure 3 c~ Figure 3 e, where Figure 3 c represents the original earthquake record, Figure 3 d is the seismic record synthesized by the inversion result, Figure 3 e is the difference between the original seismic record and the seismic record synthesized by the inversion result. It can be seen that the seismic record synthesized by the inversion result has a good agreement with the actual observation record. After inversion error evaluation and synthetic seismic record comparison, it can be seen that the present invention can accurately and stably invert the reservoir brittleness index.

[0070] Embodiment 2:

[0071] The seismic pre-stack angle gather data near a well in a certain work area is used to invert the underground reservoir brittleness index by the method of the present invention, such as Figure 4 As shown in a, the solid line represents the brittleness index curve obtained from the logging data, and the dotted line represents the brittleness index inversion result. The inversion result is basically consistent with the brittleness index curve obtained from the logging data. Figure 4 b is the brittleness index inversion error curve, and its average relative error is 1.68%; the inversion results are used to synthesize seismic records and compared with the original seismic records, such as Figure 4 c~ Figure 4 e, where Figure 4 c represents the original earthquake record, Figure 4 d is the seismic record synthesized by the inversion result, Figure 4 e is the difference between the original seismic record and the seismic record synthesized by the inversion result. It can be seen that the seismic record synthesized by the inversion result has a good agreement with the actual observation record, indicating that the present invention has a good effect in the brittleness index inversion in the actual work area.

[0072] Figure 5 The system structure diagram of the present invention includes the following modules:

[0073] The parameter acquisition module L1 is used to construct the relationship between the seismic record and the brittleness index based on the data of the well data, the wavelet data, the rock physics model and the convolution model;

[0074] The objective function construction module L2 is used to construct the objective function of brittleness index inversion based on the relationship between seismic records and brittleness index constructed based on Bayesian theory and other constraints;

[0075] The inversion method module L3 is used to perform inversion and solve the objective function calculated by the objective function construction module, and construct an inversion equation for the brittleness index. The inversion method module of this invention uses the Gauss-Newton method as the inversion method, and the inversion method can also be adjusted according to the situation;

[0076] The output module L4 substitutes the initial model obtained by the parameter acquisition module L1 into the inversion equation of the brittleness index established by the inversion method module L3, performs iterative calculations, and determines whether the convergence conditions are met; when the convergence conditions are met or the set maximum number of iterations is reached, the current result is output as the final result; otherwise, the current value is used as the initial model to continue to the next iteration until the convergence conditions are met or the set maximum number of iterations is reached.

[0077] Beneficial effects: The present invention provides a direct inversion method for prestack brittleness index based on the Zoeppritz equation, and directly inverts the prestack seismic record to obtain elastic parameters such as reservoir brittleness index, P-wave velocity, S-wave velocity, etc. The present invention obtains the relationship between the brittleness index and the elastic parameters based on the intersection analysis of the brittleness index and the elastic parameters, and then constructs the relationship between the brittleness index and the seismic record based on the Zoeppritz equation, so that the result is more accurate and the theory is more perfect; multiple constraints are used to solve the ill-conditioned inversion problem, and improve the stability and accuracy of the brittleness index estimation.

[0078] The above specific implementation methods further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above are only specific implementation methods of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A prestack brittleness index inversion method based on the Zoeppritz equation, characterized in that: The inversion method comprises: Step S1, constructing the relationship between the brittleness index and the earthquake record; Step S2, constructing an objective function of the brittleness index based on the least squares principle; Step S3, establishing an iterative inversion equation for deriving the reservoir brittleness index using a Gauss-Newton iteration method in combination with the objective function; Step S4: Given the initial model and the iterative inversion equation, an iterative operation is performed to obtain a final reservoir brittleness index.

2. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The construction of the relationship between the brittleness index and the seismic record specifically includes: Based on the intersection analysis of brittleness index and elastic parameters, the relationship between brittleness index and elastic parameters is obtained, and then based on the Zoeppritz equation, the relationship between brittleness index and seismic records is constructed.

3. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The constraint items of the objective function specifically include: a Cauchy distribution constraint item and a low-frequency model constraint item.

4. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 2, characterized in that: The relationship between the brittleness index and the elastic parameter obtained based on the intersection analysis of the brittleness index and the elastic parameter is ρ=φ(BI,V P ,V S ).

5. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 2, characterized in that: The relationship between the brittleness index and the seismic record based on the Zoeppritz equation specifically includes: Introducing the relationship between the brittleness index and the elastic parameter into the Zoeppritz equation; Obtain the mathematical relationship between the longitudinal wave reflection coefficient of the interface and the longitudinal and transverse wave velocities and brittleness index on both sides of the interface at a certain incident angle; Where p = sinθ1 / V P1 is the ray parameter, θ1 is the incident angle, BI, V P 、V S , ρ represent the brittleness index, longitudinal wave velocity, shear wave velocity, and density of the medium, respectively. Subscript 1 represents the upper medium, and subscript 2 represents the lower medium.

6. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The relationship between the brittleness index and the seismic record based on the Zoeppritz equation also includes: When there are n reflection interfaces, the P-wave reflection coefficient is expressed in the form of reservoir brittleness index and elastic parameters: P =f(θ,BI1,...,BI n+1 ,V P1 ,...,V Pn+1 ,V S1 ,...,V Sn+1 ); The form of the vector function is r(θ)=f(θ,m); in f is a forward operator, which represents the nonlinear function of the reflection coefficient sequence r with respect to the reservoir brittleness index and elastic parameter m; Considering that the seismic record is the convolution of the wavelet and the reflection coefficient, and considering that each gather has multiple records corresponding to different incident angles, the equation d=g(m) for synthesizing the seismic angle gather record using the reservoir brittleness index and elastic parameters is obtained; Wherein, d represents the seismic angle gather record, g is the nonlinear forward operator, which represents the function of the angle gather record d with respect to the vector m composed of the reservoir brittleness index and elastic parameters.

7. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The constraint expression of the objective function is: in, It represents the objective function about the reservoir brittleness index m, which consists of three parts, which represent the data error term, the Cauchy prior constraint term, and the low-frequency model constraint term from left to right; In the data error term, the superscript T is the transposition symbol; d represents the seismic observation record; m represents the reservoir brittleness index; g is the nonlinear forward operator, which represents the relationship between the seismic record and the reservoir brittleness index; W d represents the weighting matrix; In the Cauchy prior constraint, is the weight coefficient; m c represents the relative change of the reservoir brittleness index m, which conforms to the three-variable Cauchy distribution; Φ i =(D i ) T Ψ -1 D i ; Ψ is a 3×3 scale matrix; D i is a 3×3n matrix; i represents the i-th seismic record sampling point, 1≤i≤n, and n represents the number of seismic record sampling points; In the low-frequency model constraints, Indicates the weight coefficient of this item; m r represents the low-frequency model; W m represents the weighting matrix.

8. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The step S3, establishing an iterative inversion equation for deriving the reservoir brittleness index by using the Gauss-Newton iteration method in combination with the objective function, specifically comprises: Step S31: Calculate according to the objective function expression The first-order derivative of m gives the gradient of the objective function Step S32: Calculate according to the objective function gradient expression The first-order derivative of m, ignoring the higher-order derivatives in the result, obtains the Hessian matrix of the objective function Step S33: Assume that the initial model is m (0) , respectively substitute it into the gradient of the objective function and the Hessian matrix expression to calculate the objective function in m (0) The gradient at and the Hessian matrix The iterative equation of reservoir brittleness index m is obtained.

9. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 8, characterized in that: The gradient of the objective function The expression is: in 10. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 8, characterized in that: The Hessian matrix of the objective function The expression is 11. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 8, characterized in that: The iterative equation of the reservoir brittleness index m is:

12. The prestack brittleness index inversion method based on the Zoeppritz equation according to claim 1, characterized in that: The step S4, performing iterative calculations on a given initial model and the iterative inversion equation to obtain a final reservoir brittleness index specifically includes: Given an initial model m (0) , substitute into the iterative equation of reservoir brittleness index to calculate m; Determine the calculated result m and the initial model m (0) Whether the difference between is less than the set threshold or reaches the set maximum number of iterations, if so, m is the final result, otherwise m is set as the initial model for cyclic iteration.

13. A prestack brittleness index inversion system based on the Zoeppritz equation, applied to the prestack brittleness index inversion method based on the Zoeppritz equation as claimed in any one of claims 1 to 12, characterized in that: The inversion system comprises: Obtain parameter module to build the relationship between brittleness index and earthquake records; An objective function building module, used to build an objective function about the brittleness index based on the least squares principle; An inversion method module, used for establishing an iterative inversion equation for deriving a reservoir brittleness index using a Gauss-Newton iteration method in combination with the objective function; The output module is used to perform iterative calculations based on the initial model and the iterative inversion equation to obtain the final reservoir brittleness index.

Images