A “double sweet spot” seismic inversion method for unconventional fractured shale reservoirs

Abstract

The present invention discloses a "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs, including establishing a fractured gas-bearing shale reservoir rock physics model; defining a fracture gas indicator factor based on the fractured gas-bearing shale reservoir rock physics model; GFI ; Constructing a new brittleness index BI 6; Finally, the forward equation is obtained, and the inversion algorithm is used to solve it to obtain X, where X includes the fracture gas indicator factor GFI and crack density e . This application combines both anisotropy and isotropy to effectively carry out the study of shale gas "double sweet spots" seismic inversion methods. It can effectively predict the parameters of fracture-type shale gas "double sweet spots" with arbitrary fracture inclination angles through OVT domain five-dimensional seismic inversion technology, and then more accurately and reliably predict the spatial distribution of shale reservoir "double sweet spots", complete the comprehensive evaluation of the compressibility of unconventional shale reservoirs, and finally achieve the goal of correctly guiding hydraulic fracturing to increase shale gas production.

Description

Technical Field

[0001] The invention discloses a "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs, and belongs to the technical field of oil and gas field development engineering. Background Art

[0002] Shale gas sweet spots have strong anisotropy, thin reservoir thickness, and shales often have low porosity and ultra-low permeability. Hydraulic fracturing of shale reservoirs is an important way to increase gas production. Therefore, predicting the sweet spots of shale reservoirs has become an important issue. Shale gas "double sweet spots" mainly include engineering sweet spots (such as brittleness index, fracture density) and geological sweet spots (such as gas saturation). The prediction results of shale gas double sweet spots are an important basis for judging how to carry out hydraulic fracturing. Through the implementation of hydraulic fracturing, artificial fractures can be effectively connected with natural fractures to form a larger shale gas storage and transfer space. Therefore, for unconventional fractured gas-bearing shale reservoirs, how to accurately and reliably predict the shale gas "double sweet spots" is a key issue to increase shale gas production. Summary of the invention

[0003] The purpose of the present invention is to provide a "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs, so as to solve the problem that shale gas "double sweet spots" cannot be predicted for unconventional fractured gas-bearing shale reservoirs.

[0004] The embodiment of the present application is implemented as follows: a "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs includes:

[0005] Establish rock physics model of fractured gas-bearing shale reservoir;

[0006] Based on the rock physics model of fractured gas-bearing shale reservoirs, the fracture gas indicator factor is defined GFI ;

[0007] Constructing a new fragility index BI 6;

[0008] Obtain a new isotropic reflection coefficient equation including a new brittleness index:

[0009] ;

[0010] ;

[0011] ;

[0012] ;

[0013] In the formula, Indicates the brittleness index of the upper and lower strata BI The average value of 6, Indicates the Young's modulus of the upper and lower strata E The average value of Indicates the rock density of the upper and lower strata The average value of Indicates the brittleness index of the lower and upper strata BI The difference between 6, Represents the Young's modulus of the lower and upper strata E The difference between Indicates the rock density of the lower and upper strata The difference between To include the new brittleness index BI 6 isotropic longitudinal wave reflection coefficient, E is Young's modulus, is the rock density, , , Respectively represent the incident angle The relevant reflection coefficients are: is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio;

[0014] The new anisotropic TTI total reflectivity equation consisting of the fracture density and fracture gas indicator factor is obtained:

[0015] ;

[0016] ;

[0017] ;

[0018] ;

[0019] In the formula, is the total reflection coefficient of the TTI medium, is the isotropic longitudinal wave reflection coefficient including the new brittleness index, is the anisotropic longitudinal wave reflection coefficient including the fracture gas indicator factor, is the incident angle, is the azimuth, is the crack inclination angle, and Respectively represent the incident angle , azimuth and crack inclination The relevant reflection coefficients are: is the square of the ratio of the shear wave velocity to the longitudinal wave velocity, Indicates the fracture density of the lower and upper strata e The difference between Fracture gas indicator factor representing the lower and upper formations GFI The difference between

[0020] Obtain the anisotropic perturbation reflection coefficient represented by the fracture parameters :

[0021] ;

[0022] ;

[0023] ;

[0024] In the formula, is the crack inclination angle, is the azimuth, is the incident angle, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the gas indicator factor of the fractures in the lower and upper formations GFI The difference between is the fracture density of the lower and upper strata e The difference between

[0025] Assuming that the seismic data is synthesized by the anisotropic forward operator and the anisotropic reflection coefficient term, the forward equation is obtained:

[0026] ;

[0027] In the formula, represents the azimuth amplitude difference seismic data, Anisotropic forward operator, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the inversion unknown parameter matrix;

[0028] By solving the inversion algorithm, we can obtain X, where X includes the fracture gas indicator factor. GFI and crack density e .

[0029] Optionally, in some embodiments of the present application, the method for establishing a rock physics model of a fractured gas-bearing shale reservoir comprises:

[0030] The effective elastic modulus of isotropic medium mixture is calculated using the VRH average value;

[0031] Taking clay as the background, organic matter was added using the SCA model, and pores containing immobile water were added using the KT model;

[0032] The Backus averaging method was used to obtain the rock mixture with layered equivalent background medium;

[0033] Use the DEM model to add empty pore space in the rock mixture to obtain a dry rock skeleton;

[0034] A set of parallel inclined fractures were added using a linear sliding model;

[0035] The Brown-Korringa equation is used to add the gas-water mixed equivalent fluid calculated by the Wood equation into the empty pores and fracture spaces to form saturated rocks and obtain the rock physics model of fractured gas-bearing shale reservoirs.

[0036] Optionally, in some embodiments of the present application, the isotropic medium mixture includes quartz, calcite, dolomite and pyrite; and / or

[0037] The effective elastic modulus includes the equivalent bulk modulus and shear modulus of the isotropic medium mixture; and / or

[0038] The rock mixture includes a brittle mineral mixture and an organic matter mixture.

[0039] Optionally, in some embodiments of the present application, the brittle mineral mixture includes quartz, calcite, dolomite and pyrite; and / or

[0040] The organic matter mixture includes clay and organic matter.

[0041] Optionally, in some embodiments of the present application, ;

[0042] In the formula, GFI The crack gas indicator factor proposed by the present invention; e is the crack density, g b is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio, π is a constant, k fluid is the bulk modulus of the fluid, is the shear wave modulus, is the crack aspect ratio.

[0043] Optionally, in some embodiments of the present application,

[0044] ;

[0045] ;

[0046] ;

[0047] ;

[0048] ;

[0049] ;

[0050] ;

[0051] ;

[0052] ;

[0053] ;

[0054] ;

[0055] In the formula, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, M is the number of incident angles, is the incident angle, Indicates i Angle of incidence, is the crack inclination, parameter N is the number of longitudinal sampling points, is a single-parameter wavelet matrix, is the multi-parameter pre-stack seismic inversion wavelet matrix, represents the inversion unknown parameter matrix, T represents the transpose of the matrix, represents the azimuth amplitude difference seismic data, B is the coefficient matrix, express The coefficient matrix represented by the reflection coefficient term coefficients, express The coefficient matrix represented by the reflection coefficient term coefficients, D 0 is a single-parameter difference matrix, D is the difference matrix of prestack inversion, Represents the matrix composed of the fracture gas indicator factors of N sampling points, represents the matrix composed of the crack density of N sampling points, subscript sc represents the length of the wavelet, The symbol represents a diagonal matrix and T represents the transpose of the matrix.

[0056] Optionally, in some embodiments of the present application, the inversion algorithm comprises:

[0057] Use low frequency constraints and sparse constraints Construct the target inversion function:

[0058] ;

[0059] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, represents the initial parameter model, represents the inversion unknown parameter matrix, is the coefficient of the constraint term of the initial model; is the coefficient of the sparse constraint, J ( X ) represents the inversion objective function with respect to the variable X The minimum value of

[0060] Using Lagrange multipliers and the dual term The target inversion function is transformed into:

[0061] ;

[0062] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the prestack inversion coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, is the Lagrange multiplier, is the dual term, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, J ( X , H, Q) represents the minimum value of the inversion objective function with respect to variables X, H, and Q;

[0063] Solve the above formula using the multiplier alternating direction method algorithm to obtain the inversion iterative update formula:

[0064] ;

[0065] ;

[0066] ;

[0067] ;

[0068] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, superscript ( k ) indicates the k The parameter iteration result, superscript ( k +1) indicates the k +1 parameter iteration result, Indicates the maximum value, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, Represents variables The symbol function of .

[0069] Optionally, in some embodiments of the present application, Lagrange multipliers are used and the dual term The nonlinear optimization problem of the target inversion function is transformed into a linear optimization problem, that is, into:

[0070] .

[0071] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0072] 1. This invention proposes a new brittleness index expression BI 6. It is highly sensitive to the volume content of brittle minerals such as quartz, calcite, dolomite and pyrite, and to porosity, gas saturation and organic matter content. It has more advantages than traditional brittleness index and can better predict the brittleness distribution of rocks in shale gas reservoirs.

[0073] 2. This invention proposes a new crack gas indicator factor GFI , the validity of the indicator factor was verified by the constructed rock physics modeling method.

[0074] 3. The present invention derives a new longitudinal wave reflection coefficient equation, mainly based on the inclined crack density e , a new brittleness index to characterize rock brittleness BI 6 and crack gas indicator factor GFI With the core parameters as the core parameters, a direct seismic inversion method for double sweet spot parameters of fractured gas-bearing shale reservoirs is realized, which can provide advanced technical support for the development of shale gas. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 The rock physics modeling process of fractured gas shale in the present invention;

[0076] Figure 2 The results of rock physics modeling of fractured gas-bearing shale in the present invention;

[0077] Figure 3 For the present invention GFI Mechanism of action;

[0078] Figure 4 is the sensitivity of different brittleness indexes of the present invention to brittle minerals;

[0079] Figure 5 is the sensitivity of different brittleness indexes to porosity of the present invention;

[0080] Figure 6 is the sensitivity of the brittleness index to the organic matter content at different gas saturations of the present invention;

[0081] Figure 7 It is the contribution characteristics of different reflection coefficient terms to the longitudinal wave reflection coefficient of the present invention;

[0082] Figure 8 The well logging imaging data and core data related to the target layer section of Well A of the present invention;

[0083] Fig. 9 The seismic stacking data of different azimuths and different incident angles of the present invention;

[0084] Fig.10The azimuth amplitude difference seismic data of different incident angles of the present invention;

[0085] Fig.11 This is the prediction result of the double sweet spot profile through Well A of the present invention. DETAILED DESCRIPTION

[0086] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0087] The technical solution of this application is as follows:

[0088] The present application embodiment provides a "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs, including:

[0089] S01. Establish a rock physics model for fractured gas-bearing shale reservoirs;

[0090] S02. Define fracture gas indicator factors based on the rock physics model of fractured gas-bearing shale reservoirs GFI ;

[0091] S03. Constructing a new brittleness index BI 6;

[0092] S04. Obtain a new isotropic reflection coefficient equation including a new brittleness index:

[0093] ;

[0094] ;

[0095] ;

[0096] ;

[0097] In the formula, Indicates the brittleness index of the upper and lower strata BI The average value of 6, Indicates the Young's modulus of the upper and lower strata E The average value of Indicates the rock density of the upper and lower strata The average value of Indicates the brittleness index of the lower and upper strata BI The difference between 6, Represents the Young's modulus of the lower and upper strata E The difference between Indicates the rock density of the lower and upper strata The difference between To include the new brittleness index BI 6 isotropic longitudinal wave reflection coefficient, E is Young's modulus, is the rock density, , , Respectively represent the incident angle The relevant reflection coefficients are: is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio;

[0098] S05. Obtain a new anisotropic TTI total reflection coefficient equation consisting of fracture density and fracture gas indicator factor:

[0099] ;

[0100] ;

[0101] ;

[0102] ;

[0103] In the formula, is the total reflection coefficient of the TTI medium, is the isotropic longitudinal wave reflection coefficient including the new brittleness index, is the anisotropic longitudinal wave reflection coefficient including the fracture gas indicator factor, is the incident angle, is the azimuth, is the crack inclination angle, and Respectively represent the incident angle , azimuth and crack inclination The relevant reflection coefficients are: is the square of the ratio of the shear wave velocity to the longitudinal wave velocity, Indicates the fracture density of the lower and upper strata e The difference between Fracture gas indicator factor representing the lower and upper formations GFI The difference between

[0104] S06. Obtain anisotropic perturbation reflection coefficient represented by crack parameters :

[0105] ;

[0106] ;

[0107] ;

[0108] In the formula, is the crack inclination angle, is the azimuth, is the incident angle, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the gas indicator factor of the fractures in the lower and upper formations GFI The difference between is the fracture density of the lower and upper strata e The difference between

[0109] S07. Assuming that the seismic data is synthesized by the anisotropic forward operator and the anisotropic reflection coefficient term, the forward equation is obtained:

[0110] ;

[0111] In the formula, represents the azimuth amplitude difference seismic data, Anisotropic forward operator, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the inversion unknown parameter matrix;

[0112] S08. Solve by inversion algorithm to obtain X, where X includes the fracture gas indicator factor GFI and crack density e。

[0113] In the technical solution of the present application, for the prediction of fractured gas-bearing shale reservoirs with arbitrary fracture inclination angles, the present application combines both anisotropy and isotropy to effectively carry out the research on the "double sweet spots" seismic inversion method of shale gas, and realizes the "double sweet spots" seismic identification method of unconventional fractured shale gas reservoirs based on TTI media. It can effectively predict the "double sweet spots" parameters of fractured shale gas with arbitrary fracture inclination angles through OVT domain five-dimensional seismic inversion technology, and then more accurately and reliably predict the spatial distribution of the "double sweet spots" of shale reservoirs, complete the comprehensive evaluation of the compressibility of unconventional shale reservoirs, and finally achieve the goal of correctly guiding hydraulic fracturing to increase shale gas production.

[0114] In the S01:

[0115] See also Figure 1 The method for establishing a rock physics model of a fractured gas-bearing shale reservoir comprises:

[0116] S011. Calculate the effective elastic modulus of isotropic medium mixture using the VRH average value;

[0117] S012, using clay as the background, adding organic matter using the SCA model, and adding pores containing immobile water using the KT model;

[0118] S013, rock mixture with layered equivalent background medium obtained by Backus averaging method (the equivalent stiffness matrix is ​​calculated as layered equivalent background medium by Backus averaging method);

[0119] S014. Add empty pore space in the rock mixture using the DEM model to obtain a dry rock skeleton;

[0120] S015, adding a set of parallel inclined cracks using a linear sliding model;

[0121] S016. The gas-water mixed equivalent fluid calculated by the Wood equation is added into the empty pores and fracture spaces using the Brown-Korringa equation to form saturated rock and obtain a rock physics model of fractured gas-bearing shale reservoirs.

[0122] In said S011:

[0123] In some embodiments, the isotropic medium mixture includes quartz, calcite, dolomite, and pyrite.

[0124] In some embodiments, the effective elastic modulus includes an equivalent bulk modulus and a shear modulus of the isotropic medium mixture.

[0125] In said S012:

[0126] In some embodiments, the immobile water is bound water.

[0127] In said S013:

[0128] In some embodiments, the rock mixture includes a brittle mineral mixture and an organic matter mixture.

[0129] Furthermore, the brittle mineral mixture includes quartz, calcite, dolomite and pyrite.

[0130] Furthermore, the organic mixture includes clay and organic matter.

[0131] In said S016:

[0132] It can be understood that the Brown-Korringa equation is a theory describing the equivalent modulus of anisotropic rock skeletons when saturated with fluid. This equation is used to add the gas-water mixed equivalent fluid calculated by the Wood equation into the empty pores and fracture spaces to form saturated rocks, thereby estimating the elastic parameters and fracture parameters of fractured gas-bearing shale reservoirs.

[0133] Figure 2 The relevant logging curves and rock physics prediction parameter curves of Well A are shown. Figure 2 (a) Mineral composition MCC; (c) P-wave velocity VP; (d) S-wave velocity VS; (e) P-wave velocity prediction relative error REP; (f) Porosity PrT; (g) Water saturation Sw; (h) Fracture normal weakness ; (i) Crack tangential weakness ; (j) TOC content, T1, T2, T3 represent well stratification. Figure 2 It can be seen that based on Figure 1 The rock physics modeling method in Figure 2 The predicted values ​​of P-wave velocity Vp and S-wave velocity Vs in (b) and (c) (dashed lines) are in good agreement with their well logging values ​​(solid lines), which effectively confirms the accuracy of the predicted parameters of the rock physics modeling method proposed in this application. Figure 2 (h) The normal crack weakness and tangential crack weakness in (i) are reliable.

[0134] In this application, rock physics modeling of fractured gas-bearing shale reservoirs is an important bridge connecting physical parameters (porosity, mineral composition, water saturation) and elastic parameters (P-wave velocity, S-wave velocity, Young's modulus), and is an important technical support for predicting fracture parameters. Refined rock physics modeling can effectively optimize and simulate the rock characteristics of fractured gas-bearing shale reservoirs and obtain more accurate and effective rock physics analysis conclusions.

[0135] In the S02:

[0136] In some embodiments, ;

[0137] In the formula, GFI The crack gas indicator factor proposed by the present invention; e is the crack density, g b is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio, π is a constant, k fluid is the bulk modulus of the fluid, is the shear wave modulus, is the crack aspect ratio.

[0138] It can be understood that the crack gas indicator factor GFI Used to characterize the gas content in fractures.

[0139] It can be understood that based on the method of establishing a rock physics model of fractured gas-bearing shale reservoirs, a fracture gas indicator factor is defined to characterize the gas content of fractures. GFI , as shown in the formula, covers the important parameters of fracture density, fluid bulk modulus and fracture aspect ratio, which can characterize the fracture gas content and serve as an effective indicator of shale gas geological sweet spots. Figure 3 In the figure, (a) is at different oil saturation So; (b) is at different gas saturation Sg; (c) is at different fracture aspect ratios Situation. Figure 3 visible, GFI Very sensitive to fracture gas, fracture aspect ratio and fracture density have the least influence, and it is used as the geological sweet spot fracture gas indicator factor in this application GFI , which can effectively predict the fracture gas.

[0140] In said S03:

[0141] It is understandable that shale gas reservoirs usually have the characteristics of low porosity and ultra-low permeability. For actual shale reservoirs, the common brittleness index (see Table 1) BI 1- BI 5) It may not be able to characterize the rock brittleness of shale well, because the brittleness index based on geophysical elastic parameters will be affected to varying degrees by the brittle mineral components, porosity, water saturation and organic matter content of the rock. For example, the presence of porosity and organic matter will reduce the Young's modulus. Therefore, based on the established rock physics model of fractured gas-bearing shale reservoirs, this application defines a new brittleness index that can better characterize the brittleness of shale. Experiments show that it is a new brittleness index that is sensitive to brittle mineral components, especially porosity, gas saturation and organic matter content. BI 6. More accurate brittleness prediction results can be obtained, providing more reliable engineering sweet spot support for shale gas hydraulic fracturing.

[0142] Further, from the brittle mineral components (see Figure 4 ), porosity (see Figure 5 ), organic matter content (see Figure 6 ) demonstrated that the new brittleness index BI6 showed positive correlation sensitivity characteristics and could better characterize the rock brittleness of shale reservoirs. Figure 4 In the figure: (a) is quartz; (b) is calcite; (c) is dolomite; (d) is the sensitivity value weighted by component content; Figure 5Middle: (a) is quartz 30%; (b) is quartz 40%; (c) is quartz 50%; Figure 6 In: (a) BI 3; (b) BI 5; (c) BI 6. Against the traditional brittleness index BI 1- BI 5 and the brittleness index proposed by the present invention BI 6, respectively verify their changing rules with the volume content of quartz, calcite and dolomite. Figure 4 (a) It can be seen that when the quartz volume content increases, the six brittleness indices increase in a positive correlation, and BI 5 and BI 6The change is greater than other BI 1- BI 4 greater; by Figure 4 (b) It can be seen that when the volume content of calcite increases, the six brittleness indices all increase in a positive correlation, and at this time, BI2 changes the most; Figure 4 (c) It can be seen that when the dolomite content increases, the six brittleness indices all increase in positive correlation, and BI2 changes the most. However, in practical applications, it is necessary to consider the contribution of the actual content of quartz, calcite, and dolomite to the brittleness index. Figure 4 In (d), the volume contents of quartz, calcite and dolomite were weighted to obtain the brittle mineral sensitivity coefficient values ​​of the six brittleness indices, and the preliminary optimization was carried out. BI 3. BI 5. BI 6. In order to continue to verify the BI 6 Ratio BI 3 and BI 5 advantages, in Figure 5 (a), (b) and (c) discuss the cases where the quartz volume content is 30%, 40% and 50%, respectively. As the rock porosity increases, BI 3 shows a negative correlation decrease, which does not meet the requirements, while BI 5 and BI 6 shows a positive correlation increase, and BI 6The increase is even greater, indicating that BI 6 Ratio BI 5 is more sensitive to rock porosity and has more advantages. BI 3. BI 5 and BI 6. The characteristics of the change of organic matter content in the case of saturated water and saturated oil. Figure 6 (a) It can be seen that BI 3 In the case of saturated water and saturated oil, with the increase of organic matter content, the decrease is negatively correlated, indicating that it is not sensitive to organic matter content and there is a deficiency; Figure 6 (b) It can be seen that BI 5 In the case of saturated water and saturated oil, with the increase of organic matter content, there is a positive correlation increase, indicating that there is a certain sensitivity to organic matter content; Figure 6 (c) It can be seen that BI 6 In the case of saturated water and saturated oil, with the increase of organic matter content, the increase is positively correlated, and the change range is Figure 6 (b) BI 5 is greater, indicating BI 6 Ratio BI 5. More sensitive to organic matter content.

[0143] Table 1 Traditional brittleness index and new brittleness index

[0144]

[0145] In Table 1: represents the rock density, E is Young's modulus, is Poisson's ratio, E b is Young's modulus brittleness, is the Poisson's ratio brittleness, and the parameter subscripts min and max represent the minimum and maximum values ​​of the parameters, respectively. E min represents the minimum value in Young's modulus, E max represents the maximum value in Young's modulus, represents the minimum value in Young's modulus, represents the maximum value in Young's modulus, M is the longitudinal wave modulus, It is the first Latin American coefficient.

[0146] Specifically: Rickman, R., Mullen, M., Petre, E., et al., 2008. APractical Use of Shale Petrophysics for Stimulation Design Optimization: AllShale Plays Are Not Clones of the Barnett Shale. SPE Annual Technical Conference&Exhibition Society of Petroleum Engineers.doi.org / 10.2118 / 115258-MS.

[0147] Sharma and Chopra, etc. Specifically: Sharma, RK, Chopra, S., 2012. Newattribute for determination of lithology and brittleness. SEG. Tech. Program.Expand. Abstr.doi.org / 10.1190 / segam2012-1389.1.

[0148] Guo et al. Specifically: Guo, ZQ, Chapman, M., Li, XY, 2012. A shale rockphysics model and its application in the prediction of brittleness index, mineralalogy, and porosity of the Barnett Shale. SEG. Tech. Program. Expand.Abstr. doi.org / 10.1190 / segam2012-0777.1.

[0149] Chen et al. are specifically: Chen, JJ, Zhang, GZ, Chen, HZ, et al., 2014. The Construction of Shale Rock Physics Effective Model and Prediction of RockBrittleness. SEG. Tech. Program. Expand. Abstr. doi.org / 10.1190 / segam2014-0716.1.

[0150] In said S04:

[0151] It can be understood that for the isotropic background reflection coefficient term, a quantitative relationship between the new brittleness index and the original elastic parameters (P, S wave velocities and density) is established, and a new isotropic reflection coefficient equation containing the engineering sweet spot (new brittleness index) is derived.

[0152] In said S05:

[0153] It can be understood that for the anisotropic TTI reflection coefficient fracture disturbance term, a quantitative relationship between the engineering sweet spot fracture density, the geological sweet spot fracture gas indicator factor and the fracture weakness is established, and finally a new anisotropic TTI total reflection coefficient equation consisting of the double sweet spot fracture density and the fracture gas indicator factor is derived.

[0154] In said S06:

[0155] It is understandable that due to the weak anisotropy, it is necessary to introduce azimuth differences to eliminate the isotropic background component and reduce the number of unknown parameters in one-time synchronous inversion, which helps to improve the stability and accuracy of the anisotropic parameter inversion problem. Based on the analysis of the contribution of the above parameters to the PP wave reflection coefficient, the sweet spot parameters can be effectively decoupled (see Figure 7 ), Figure 7 middle, is the crack inclination angle, is the azimuth, is the incident angle, from Figure 7 It can be seen that the reflection coefficient terms of different elastic parameters and the reflection coefficient terms of crack parameters have different characteristics with the change of azimuth and incident angle. Figure 7 (a), (b) and (c) correspond to Young's modulus respectively E , brittleness index BI 6 and the reflection coefficient term of rock density ρ do not change with azimuth and do not have azimuthal anisotropy characteristics; Figure 7 (d)-(f) correspond to the crack inclination angles v The crack gas indicator factor at 0°, 45°, and 90° GFI The reflection coefficient term varies with the azimuth and the incident angle, and has azimuthal anisotropy. The greater the crack inclination, the greater the azimuthal anisotropy intensity. Figure 7 (g)-(h) correspond to the crack inclination angles v Crack density at 0°, 45°, and 90° e The reflection coefficient term varies with the azimuth and the angle of incidence, and also has azimuthal anisotropy. The greater the fracture inclination, the greater the azimuthal anisotropy intensity. Therefore, these differences can be used to perform seismic inversion on these parameters. At the same time, this application optimizes the two-step inversion strategy to improve the accuracy and stability of the inversion. The isotropic parameter component is eliminated based on the azimuthal amplitude difference seismic data between the fracture development direction and its orthogonal direction, thereby constructing an anisotropic perturbation reflection coefficient represented by the fracture parameters. .

[0156] It can be understood that the fracture parameters include the fracture gas indicator factor GFI and crack density e .

[0157] In said S07:

[0158] It can be understood that in order to invert the relevant fracture parameters using the azimuth amplitude difference seismic data, it is assumed that the seismic data is synthesized by the anisotropic forward operator and the anisotropic perturbation P-wave reflection coefficient term, and the forward equation is obtained, where:

[0159] ;

[0160] ;

[0161] ;

[0162] ;

[0163] ;

[0164] ;

[0165] ;

[0166] ;

[0167] ;

[0168] ;

[0169] ;

[0170] In the formula, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, M is the number of incident angles, is the incident angle, Indicates i Angle of incidence, is the crack inclination, parameter N is the number of longitudinal sampling points, is a single-parameter wavelet matrix, is the multi-parameter pre-stack seismic inversion wavelet matrix, represents the inversion unknown parameter matrix, T represents the transpose of the matrix, represents the azimuth amplitude difference seismic data, B is the coefficient matrix, express The coefficient matrix represented by the reflection coefficient term coefficients, express The coefficient matrix represented by the reflection coefficient term coefficients, D 0 is a single-parameter difference matrix, D is the difference matrix of prestack inversion, Represents the matrix composed of the fracture gas indicator factors of N sampling points, represents the matrix composed of the crack density of N sampling points, subscript sc represents the length of the wavelet, The symbol represents a diagonal matrix and T represents the transpose of the matrix.

[0171] In said S08:

[0172] The inversion algorithm comprises:

[0173] S081, use low frequency constraints and sparse constraints Construct the target inversion function:

[0174] ;

[0175] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, represents the initial parameter model, represents the inversion unknown parameter matrix, is the coefficient of the constraint term of the initial model; is the coefficient of the sparse constraint, J ( X ) represents the inversion objective function with respect to the variable X The minimum value of

[0176] S082. Using Lagrange multipliers and the dual term Convert the formula in step S081 into:

[0177] ;

[0178] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the prestack inversion coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, is the Lagrange multiplier, is the dual term, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, J ( X , H, Q ) represents the minimum value of the inversion objective function with respect to variables X, H, and Q.

[0179] S083, using the alternating direction method of multipliers (ADMM) algorithm to solve the formula in step S082 to obtain the inversion iterative update formula:

[0180] ;

[0181] ;

[0182] ;

[0183] ;

[0184] In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, superscript ( k ) indicates the k The parameter iteration result, superscript ( k +1) indicates the k +1 parameter iteration result, Indicates the maximum value, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, Represents variables The symbol function of .

[0185] It can be understood that by iterating the formula in step S083, the unknown parameter X within the allowable error range is finally solved. At this time, X includes the crack gas indicator factor GFI and crack density e .

[0186] Understandably, represents the initial parameter model, It represents the inversion unknown parameter matrix. The initial parameter model is obtained by interpolating and extrapolating the logging curve and then performing Gaussian low-pass filtering.

[0187] It can also be understood that using Lagrange multipliers and the dual term The nonlinear optimization problem of the formula in step S081 is converted into a linear optimization problem, that is, converted into:

[0188] ;

[0189] It can also be understood that the longitudinal wave reflection coefficient part of the isotropic parameter is first calculated Then the same inversion algorithm can be used to solve the new brittleness index BI 6 and Young's modulus E and rock density .

[0190] Application Examples

[0191] To verify the effectiveness of this method in predicting fracture development, fracture gas content and shale brittleness distribution in fractured gas-bearing shale reservoirs, the actual working area is located in a block in the Sichuan Basin of China. Figure 4 The well logging data of Well A and the prediction results of fracture density and brittleness on the well are shown. The T2 to T3 sections of Well A are fractured high-quality gas-bearing shale reservoirs with vertical depths of 2483 m to 2516 m. The corresponding time at the bottom of the reservoir is 1.78 s. The reservoir has low porosity and ultra-low permeability characteristics. FMI logging imaging (see Figure 8 ), Figure 8 The FMI borehole imaging in (a) shows the development of inclined fractures in the reservoir. Figure 8 (b) The development of fractures was observed in the core sampling of the well. Figure 8 (c) shows that the crack development direction is southeast-northwest. Figure 8 (d) shows that the fracture dip angle is unevenly distributed, mainly concentrated at 70°. The high organic matter content indicates that the shale in this section has a high hydrocarbon generation capacity, a high brittleness index, and high fracture density and fracture gas indicator factor values. Double sweet spot prediction is carried out for any survey line profile passing through Well A. First, two sets of azimuth gathers are extracted along the fracture development direction and orthogonal to the direction, and then superimposed at three different central incidence angles to improve the signal-to-noise ratio of seismic data (see Fig. 9 ), Fig. 9 (a), (b) and (c) show the position 1: φ obs1 The incident angles are obtained at 45°. θ The seismic stacking data are 5°, 15° and 25° respectively; Fig. 9 (d), (e) and (f) show the position 2: φ obs2 The incident angles are obtained at 135°. θ Seismic stacking data of 5°, 15° and 25° are used to make subtractions of seismic data at different azimuths with the same incident angle to obtain three sets of azimuth amplitude difference seismic data with different incident angles (see Fig.10 ), Fig.10 (a) shows the azimuthal amplitude difference seismic data with an incident angle of 5°; Fig.10 (b) shows the azimuthal amplitude difference seismic data with an incident angle of 15°; Fig.10 (c) shows the azimuth amplitude difference seismic data with an incident angle of 25°. In the first step, the azimuth amplitude difference seismic data is used to predict the anisotropic sweet spot parameters, and the profile inversion results of the fracture density and fracture gas indicator factor through Well A are obtained. In the second step, the residual seismic data with the anisotropic part eliminated is used to predict the isotropic sweet spot parameters, and the profile inversion results of the new high sensitivity brittleness index are obtained (see Fig.11 ), Fig.11 (a) The brittleness index predicted by the method proposed in the present invention is shown BI 6 Sweet spot parameters, it can be seen that the brittleness index detected high abnormal values ​​in the reservoir section, and the shale brittleness is good; Fig.11 (b) shows the fracture effective gas indicator factor predicted by the method proposed in this invention GFI The sweet spot parameters show that the reservoir section exhibits high value anomalies, which is consistent with the logging curve characteristics of Well A, indicating that the reservoir may contain a large amount of shale gas; Fig.11 (c) shows the fracture density sweet spot parameters predicted by the method proposed in the present invention, which also show high value anomalies in the reservoir section, indicating that fractures are developed in the reservoir, providing a certain storage space for oil and gas. It can be seen that the reservoir section (the bold position of the black well trajectory) is predicted to have developed fractures, high brittleness, and unevenly distributed shale gas.

[0192] In the above application example:

[0193] Figure 8 Middle: (a) FMI wellbore imaging, (b) target layer core, (c) fracture orientation rose diagram, (d) fracture dip statistical histogram;

[0194] Fig. 9 middle: is the azimuth, is the angle of incidence;

[0195] Fig.10 middle: is the azimuth, is the angle of incidence;

[0196] Fig.11 Middle: (a) brittleness index, (b) crack gas indicator factor, (c) crack density.

[0197] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs, characterized in that: include: Establish rock physics model of fractured gas-bearing shale reservoir; Based on the rock physics model of fractured gas-bearing shale reservoirs, the fracture gas indicator factor is defined GFI ; Constructing a new fragility index BI 6; Obtain a new isotropic reflection coefficient equation including a new brittleness index: ； ； ； ； In the formula, Indicates the brittleness index of the upper and lower strata BI The average value of 6, Indicates the Young's modulus of the upper and lower strata E The average value of Indicates the rock density of the upper and lower strata The average value of Indicates the brittleness index of the lower and upper strata BI The difference between 6, Represents the Young's modulus of the lower and upper strata E The difference between Indicates the rock density of the lower and upper strata The difference between To include the new brittleness index BI 6 isotropic longitudinal wave reflection coefficient, E is Young's modulus, is the rock density, , , Respectively represent the incident angle The relevant reflection coefficients are: is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio; The new anisotropic TTI total reflectivity equation consisting of the fracture density and fracture gas indicator factor is obtained: ； ； ； ； In the formula, is the total reflection coefficient of the TTI medium, is the isotropic longitudinal wave reflection coefficient including the new brittleness index, is the anisotropic longitudinal wave reflection coefficient including the fracture gas indicator factor, is the incident angle, is the azimuth, is the crack inclination angle, and Respectively represent the incident angle , azimuth and crack inclination The relevant reflection coefficients are: is the square of the ratio of the shear wave velocity to the longitudinal wave velocity, Indicates the fracture density of the lower and upper strata e The difference between Fracture gas indicator factor representing the lower and upper formations GFI The difference between Obtain the anisotropic perturbation reflection coefficient represented by the fracture parameters : ； ； ； In the formula, is the crack inclination angle, is the azimuth, is the incident angle, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the gas indicator factor of the fractures in the lower and upper formations GFI The difference between is the fracture density of the lower and upper strata e The difference between Assuming that the seismic data is synthesized by the anisotropic forward operator and the anisotropic reflection coefficient term, the forward equation is obtained: ； In the formula, represents the azimuth amplitude difference seismic data, Anisotropic forward operator, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the inversion unknown parameter matrix; By solving the inversion algorithm, we can obtain X, where X includes the fracture gas indicator factor. GFI and crack density e .

2. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 1, characterized in that: The method for establishing a rock physics model of a fractured gas-bearing shale reservoir comprises: The effective elastic modulus of isotropic medium mixture is calculated using the VRH average value; Taking clay as the background, organic matter was added using the SCA model, and pores containing immobile water were added using the KT model; The Backus averaging method was used to obtain the rock mixture with layered equivalent background medium; Use the DEM model to add empty pore space in the rock mixture to obtain a dry rock skeleton; A set of parallel inclined fractures were added using a linear sliding model; The Brown-Korringa equation is used to add the gas-water mixed equivalent fluid calculated by the Wood equation into the empty pores and fracture spaces to form saturated rocks and obtain the rock physics model of fractured gas-bearing shale reservoirs.

3. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 2, characterized in that: The isotropic medium mixture comprises quartz, calcite, dolomite and pyrite; and / or The effective elastic modulus includes the equivalent bulk modulus and shear modulus of the isotropic medium mixture; and / or The rock mixture includes a brittle mineral mixture and an organic matter mixture.

4. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 3, characterized in that: The brittle mineral mixture includes quartz, calcite, dolomite and pyrite; and / or The organic matter mixture includes clay and organic matter.

5. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 1, characterized in that: ； In the formula, GFI The crack gas indicator factor proposed by the present invention; e is the crack density, g b is the shear wave velocity V S and the longitudinal wave velocity V P The square of the ratio, π is a constant, k fluid is the bulk modulus of the fluid, is the shear wave modulus, is the crack aspect ratio.

6. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 1, characterized in that: ； ； ； ； ； ； ； ；； ； ； ； In the formula, are two azimuth angles selected based on the fracture development orientation ( ) of the amplitude difference seismic data, To eliminate the anisotropic disturbance of the longitudinal wave reflection coefficient after the isotropic background, M is the number of incident angles, is the incident angle, Indicates i Angle of incidence, is the crack inclination, parameter N is the number of longitudinal sampling points, is a single-parameter wavelet matrix, is the multi-parameter pre-stack seismic inversion wavelet matrix, represents the inversion unknown parameter matrix, T represents the transpose of the matrix, represents the azimuth amplitude difference seismic data, B is the coefficient matrix, express The coefficient matrix represented by the reflection coefficient term coefficients, express The coefficient matrix represented by the reflection coefficient term coefficients, D 0 is a single-parameter difference matrix, D is the difference matrix of prestack inversion, Represents the matrix composed of the fracture gas indicator factors of N sampling points, represents the matrix composed of the crack density of N sampling points, subscript sc represents the length of the wavelet, The symbol represents a diagonal matrix and T represents the transpose of the matrix.

7. The "double sweet spot" seismic inversion method for unconventional fractured shale reservoirs according to claim 1, characterized in that: The inversion algorithm comprises: Use low frequency constraints and sparse constraints Construct the target inversion function: ； In the formula, and They represent L2 norm and L P norm, Represents the mathematical minimum symbol, represents the initial parameter model, represents the inversion unknown parameter matrix, is the coefficient of the constraint term of the initial model; is the coefficient of the sparse constraint, J ( X ) represents the inversion objective function with respect to the variable X The minimum value of Using Lagrange multipliers and the dual term The target inversion function is transformed into: ； In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the prestack inversion coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, is the Lagrange multiplier, is the dual term, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, J ( X , H, Q ) represents the minimum value of the inversion objective function with respect to variables X, H, and Q; Solve the above formula using the multiplier alternating direction method algorithm to obtain the inversion iterative update formula: ； ； ； ； In the formula, and denote the L2 norm and L P norm, Represents the mathematical minimum symbol, is the multi-parameter pre-stack seismic inversion wavelet matrix, B is the coefficient matrix, D is the difference matrix of prestack inversion, represents the initial parameter model, represents the inversion unknown parameter matrix variable, are two azimuth angles selected based on the fracture development orientation ( and ) of the amplitude difference seismic data, superscript ( k ) indicates the k The parameter iteration result, superscript ( k +1) indicates the k +1 parameter iteration result, Indicates the maximum value, represents the inversion unknown parameter matrix, is the low frequency constraint coefficient, is the coefficient of the sparse constraint, is the dual term coefficient, Represents variables The symbol function of .

Images