High-precision prediction method for few-well-area interlayer type shale oil enrichment layer parameters

Abstract

The invention discloses a high-precision prediction method for parameters of a few-well-area interlayer type shale oil enrichment layer. The method mainly solves the problem that the prediction accuracy of the dessert layer is difficult to meet the dessert prediction requirement of the interlayer type shale oil under the conditions of thin shale single layer, fast lithologic combination transverse change and sparse drilling. The method comprises the following steps: S1, obtaining a sub-offset partial superposition body for inversion; s2, obtaining an elastic parameter body of deterministic inversion; s3, obtaining a final three-dimensional lithofacies probability body and an elastic parameter body of three-dimensional lithofacies trend constraint geostatistics inversion; s4, converting the elastic parameter body inverted by the three-dimensional lithofacies trend constraint geostatistics into an enrichment layer parameter body by using a multiple regression analysis method; and S5, according to an enrichment layer classification evaluation standard, carrying out enrichment layer classification prediction on the obtained enrichment layer parameter body. According to the method, the prediction accuracy of the dessert layer under the conditions that the shale single layer is thin, the lithology combination transverse change is fast and drilling wells are sparse can be improved, and the demand of interlayer type shale oil dessert prediction is met.

Description

Technical field:

[0001] The present invention relates to the technical field of unconventional oil and gas exploration, and in particular to a method for high-precision prediction of parameters of interlayer shale oil enrichment layers in areas with few wells. Background technology:

[0002] In the shale gas exploration and development areas with great potential, the shales are generally widely distributed and have good reservoir properties. However, the lithology in such areas changes rapidly laterally, the single layers are thin, and the shale is interbedded with sandstone and mudstone, resulting in poor lateral continuity of the seismic reflection axis as a whole. In addition, there are few wells drilled in the study area, and the resolution of deterministic inversion is insufficient. The conventional geostatistical inversion prediction results are prone to well circle phenomenon, and often have large errors, low resolution, poor lateral prediction, and it is difficult to effectively identify interlayer thin shale oil (such as Figure 2 b). Parameter prediction relies only on a single parameter and has low accuracy. The particularity and complexity of interlayer thin shale reservoirs in areas with few wells place higher demands on seismic inversion prediction. It is necessary to take into account both the vertical resolution and the lateral extension trend of the reservoir to accurately predict various elastic parameters, thereby improving the prediction accuracy of the sweet spot layer to meet the needs of interlayer shale oil sweet spot prediction. Summary of the invention:

[0003] The present invention aims to solve the problem that the prediction accuracy of the sweet spot layer in the case of thin shale single layer, rapid lateral change of lithologic combination and sparse drilling in the background technology is difficult to meet the prediction demand of interlayer shale oil sweet spot, and provide a high-precision prediction method for parameters of interlayer shale oil rich layer in a well-sparse area. The high-precision prediction method for parameters of interlayer shale oil rich layer in a well-sparse area can improve the prediction accuracy of the sweet spot layer in the case of thin shale single layer, rapid lateral change of lithologic combination and sparse drilling, and meet the prediction demand of interlayer shale oil sweet spot.

[0004] The present invention solves the problem through the following technical solution: the method for high-precision prediction of parameters of interlayer shale oil-rich layers in a few-well area comprises the following steps:

[0005] S1. After optimizing the pre-stack CRP gathers in the work area, the actual gathers with smaller deviation from the AVO characteristics of the underground medium are obtained and partially stacked to obtain the offset partial stacking volume used for inversion;

[0006] S2. Perform well-seismic calibration on the obtained offset stack, extract multi-well sub-waves, establish a well interpolation low-frequency model, and then perform pre-stack deterministic inversion to obtain an elastic parameter body of deterministic inversion;

[0007] S3. Use the elastic parameter body obtained by deterministic inversion to perform lithofacies probability analysis, and combine it with the well lithology curve for comprehensive analysis to obtain the final three-dimensional lithofacies probability body; use the obtained three-dimensional lithofacies probability body to constrain geostatistical inversion, and obtain the elastic parameter body of three-dimensional lithofacies trend constraining geostatistical inversion;

[0008] S4. Using multiple regression analysis, the elastic parameter volume obtained by geostatistical inversion of three-dimensional lithofacies trend constraints is converted into an enrichment layer parameter volume;

[0009] S5. According to the enrichment layer classification evaluation standard, the obtained enrichment layer parameter body is subjected to enrichment layer classification prediction.

[0010] Furthermore, step S1 performs gather optimization processing on the pre-stack CRP gathers of the work area by using a set of technical processes such as Radon transform denoising, spectrum balancing, and elimination of fitting amplitude compensation to perform gather optimization processing;

[0011] The specific method of the optimization process is:

[0012] The parabola Radon transform filtering technology is used to separate the effective signal from the noise in the transform domain, and then the inverse transform is used to restore the seismic gathers; the separation degree between the primary wave and the noise is determined by the moving average time difference between the primary wave and the noise; after processing, the multiple waves are eliminated, the noise is suppressed, and it is closer to the forward gathers;

[0013] The spectrum balancing technique is used to perform time-frequency analysis on the gathers based on Gabor transform. The near-mid-offset gathers or small-angle gathers that match the AVO characteristics of the well are used as standard gathers. The average time-frequency spectrum of the reference gathers is used as the expectation to compensate for the high-frequency energy of the far-offset gathers.

[0014] The technique of eliminating fitting amplitude compensation is adopted to select the trace with correct AVO law as the standard trace. The data points of amplitude variation with angle of the standard trace are fitted by least square method. The overdetermined equation is solved to obtain the intercept and gradient, and the relationship between amplitude variation with angle is obtained. The relationship between amplitude variation with angle is then extrapolated to the near-offset and far-offset trace gathers to achieve the effect close to the real AVO law and improve the inversion result of P-wave and S-wave velocity ratio.

[0015] Furthermore, the method of performing prestack deterministic inversion to obtain the elastic parameter volume in step S2 includes:

[0016] 201) Perform well seismic calibration on the wells in the work area and extract multi-well sub-waves corresponding to each part of the stacked body;

[0017] 202) Using the stratigraphic framework and well curves, a low-frequency model of well difference of P-wave impedance and P-S wave velocity ratio is established for prestack deterministic inversion by well interpolation method;

[0018] 203) Using the low-frequency model of well difference, pre-stack simultaneous inversion is performed on the stacked bodies at different offsets to generate elastic parameter bodies of P-wave impedance and P-S wave velocity ratio;

[0019] 204) The generated elastic parameter bodies of P-wave impedance and P-S wave velocity ratio are used to determine the approximate distribution range of the reservoir and provide parameters such as lateral range and sand-to-ground ratio in the variation function for geostatistical inversion.

[0020] Furthermore, S3 uses the elastic parameter body obtained by deterministic inversion to perform lithofacies probability analysis, and combines the well lithology curve for comprehensive analysis to obtain the final three-dimensional lithofacies probability body, including:

[0021] 301) The elastic parameter curves of the wellbore P-wave impedance, P-wave velocity ratio and lithology curves are intersected and analyzed to establish the response range and quantitative probability density function of the elastic parameter body obtained by deterministic inversion corresponding to different lithologies;

[0022] 302) Based on the Bayesian discriminant principle, lithofacies probability analysis (FFP) is performed, and the elastic parameter body is transformed into a three-dimensional lithofacies probability body of each major lithofacies type through cloud transformation using the probability density function.

[0023] Further, step S3 uses the three-dimensional lithofacies probability body to constrain geostatistical inversion to obtain a pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body, comprising the following steps:

[0024] 401) Obtaining respective prior probability density functions based on the input three-dimensional lithofacies probability volume and seismic data statistics;

[0025] 402) fusing the above-mentioned individual probability density functions through Bayesian reasoning to obtain a posterior probability distribution function of reservoir parameters under known seismic, well logging and geological information;

[0026] 403) The Markov chain-Monte Carlo inversion algorithm is used to obtain the statistically correct random sample point distribution according to the above-mentioned posterior probability distribution function, and the realizations of each reservoir parameter to be determined are calculated to obtain the pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body.

[0027] Furthermore, in step S4, a method for converting the elastic parameter body obtained by geostatistical inversion of three-dimensional lithofacies trend constraints into an enrichment layer parameter body by using a multivariate regression analysis method includes:

[0028] The multivariate regression analysis method is used to fit the curves of wellbore elastic parameters and enrichment layer parameters into a linear regression relationship.

[0029] The linear regression relationship formula between the wellbore elastic parameters and the enriched layer parameters is applied to the pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body to convert it into an enriched layer parameter body.

[0030] Furthermore, the relationship between the elastic parameter and the enrichment layer parameter can be expressed as follows:

[0031] I p =a*TOC+b*POR+c*V sh (1)

[0032] I s =a′*TOC+b′*POR+c′*V sh (2)

[0033] λρ=a″*TOC+b″*POR+c″*V sh (3)

[0034] Where: Ip is the longitudinal wave impedance; Is is the shear wave impedance; λρ is the product of the compression modulus and density (sensitive to fluid); TOC is the total organic carbon; POR is the porosity; Vsh is the mud content; a, a′, a″, b, b′, b″, c, c′, c″ are all constants.

[0035] Furthermore, the above-hole elastic parameters include longitudinal wave impedance, shear wave impedance, and λρ; the above-mentioned enriched layer parameters include porosity, TOC, and mud content.

[0036] Further, the elastic parameters affecting TOC include Ip for longitudinal wave impedance, Is for transverse wave impedance; λρ for the product of compression modulus and density, so the elastic parameters of TOC are respectively related to the longitudinal wave impedance, transverse wave impedance, and the product of compression modulus and density;

[0037] The elastic parameter regression relationship of TOC is expressed as follows:

[0038] TOC=f*I p +g*I s +h*λρ (5)

[0039] Wherein: d, e, f, g, h all represent constant values.

[0040] Compared with the above background technology, the present invention has the following beneficial effects:

[0041] The high-precision prediction method for parameters of interlayer-type shale oil-rich layers in a few-well area of ​​the present invention has the following beneficial effects: after conventional processing, the prestack CRP gathers still have problems such as low signal-to-noise ratio, low resolution, poor consistency of seismic trace frequencies at different offsets, weak near-far trace energy, and poor lateral continuity. These problems lead to a large deviation between the AVO (Amplitude variation with offset) characteristics of the actual gathers and the AVO characteristics of the underground medium. In order to meet the requirements of prestack elastic parameter inversion and restore the AVO characteristics of the gathers, it is necessary to perform a series of optimization treatments on the prestack CRP gathers (common reflection point gathers). The gathers after optimization are compared to obtain certain improvements in signal-to-noise ratio, resolution, and AVO fidelity, and are partially superimposed to obtain partial superpositions for subsequent inversion. Compared with the gathers before optimization, the wavelets extracted from each superposition after gather optimization have improved wavelet symmetry, more stable wavelet phases within the seismic frequency band, and fewer wells are eliminated when extracting multi-well comprehensive wavelets.

[0042] The lithofacies trend constrained inversion uses deterministic inversion results to perform lithofacies probability analysis, and deterministic inversion is completely in accordance with earthquakes. Therefore, this method takes into account the lateral predictability of deterministic inversion, while retaining the high vertical resolution of geostatistical inversion, combining the advantages of the two. The defect of conventional geostatistical methods is that in areas with few wells, it is difficult to give the lithofacies ratio of the statistical model, the inversion results are greatly affected by the stratigraphic framework, and for reservoirs with fast lateral changes, the inversion effect has poor correspondence with the distribution and extension trend of the actual reservoir. Therefore, relative to conventional geostatistical methods, the lithofacies trend constrained inversion used in the present invention effectively solves the problem of high-precision inversion of interlayer thin shale in areas with few wells. And the actual work area data show that the results of elastic parameter inversion are well consistent with the well curve, and the thin layer resolution is increased from 8m to 3m. Compared with the pre-stack deterministic inversion results, the resolution is higher, and it is generally consistent in reflecting the distribution of sand bodies, the scale of geological bodies, and the proportion of lithology. Compared with the results of conventional prestack geostatistical inversion, the resolution is comparable, and the inversion effect corresponds better to the distribution and extension trend of the actual reservoir. It can be seen that the three-dimensional lithofacies trend-constrained prestack geostatistical inversion is more suitable for situations where drilling is sparse and lithology changes rapidly.

[0043] The multivariate regression analysis method used in the method of the present invention fits the curve of the wellbore elastic parameters and the enrichment layer parameters into a linear regression relationship, and applies the relationship formula to the high-precision prestack elastic parameter body obtained by the prestack geostatistical inversion of the previous three-dimensional lithofacies trend constraint, and converts it into an enrichment layer parameter body. Multivariate regression uses a method of jointly predicting multiple influencing elastic parameters, which is more accurate than a single parameter.

[0044] In summary, this method mainly uses three-dimensional lithofacies trend constrained geostatistics on the basis of rock physics analysis to obtain high-precision prestack inversion results that are closer to reality, and then obtain elastic parameter inversion results with higher accuracy. This method avoids the shortcomings of strong randomness of lateral prediction of conventional prestack geostatistical inversion in areas with few wells and low vertical resolution of conventional prestack deterministic inversion, and innovatively proposes a strategy to improve inversion accuracy by using lithofacies constraints, and further improve the prediction accuracy of enrichment layer parameters by using multivariate linear regression parameter prediction methods. The method of the present invention was verified by taking the parameter prediction of interlayer-type thin shale oil enrichment layers of Lianggaoshan Formation in Sichuan and Chongqing PC exploration area as an example. Based on the three-dimensional lithofacies trend constrained geostatistical inversion technology, taking into account both vertical resolution and lateral predictability, the sweet spot layer recognition ability was improved from 8 meters to 3 meters, and the prediction compliance rate of thin layers above 3 meters reached 81%. Then, the enrichment layer parameter prediction method of multivariate regression analysis was applied, and the prediction compliance rate of the main control parameters of the enrichment layer reached 76.5%, supporting the deployment of 8 horizontal wells of Jurassic shale oil. Description of the drawings:

[0045] Figure 1 A comparison diagram of the P-wave impedance inversion profile of the PA1 well 3D lithofacies trend constrained prestack geostatistical inversion (a) and the conventional prestack deterministic inversion (b) in the embodiment of the present invention;

[0046] Figure 2 A comparative diagram of the P-wave impedance inversion profile of the PA1 well 3D lithofacies trend constrained prestack geostatistical inversion (a) and the conventional prestack geostatistical inversion (b) according to an embodiment of the present invention;

[0047] Figure 3 This is a flow chart of the method for high-precision prediction of parameters of interlayer-type shale oil-rich layers in areas with few wells in the present invention;

[0048] Figure 4 This is a flow chart of a three-dimensional lithofacies trend-constrained prestack geostatistical inversion technique according to an embodiment of the present invention;

[0049] Figure 5 This is a probability volume profile of mudstone, sandstone and shale lithofacies in the direction of the PA1 well trajectory according to an embodiment of the present invention;

[0050] Figure 6 It is a comparison diagram of the TOC cross-plot predicted by the single parameter of longitudinal wave impedance and the multiple regression of elastic parameters in the embodiment of the present invention;

[0051] Figure 7 This is a TOC prediction profile of elastic parameters multivariate regression according to an embodiment of the present invention;

[0052] Figure 8 It is a comparison diagram of the cross-plot of porosity predicted by single parameter of longitudinal wave impedance and multiple regression of elastic parameters in the embodiment of the present invention;

[0053] Fig. 9 This is a cross-sectional diagram of porosity prediction by multiple regression of elastic parameters according to an embodiment of the present invention;

[0054] Fig.10 Prediction profile of brittleness index for multiple regression of elastic parameters;

[0055] Fig.11 This is a predicted cross-section of the enriched layer in the direction of the PA1 well trajectory according to an embodiment of the present invention;

[0056] Fig.12 This is a plan view of the thickness of the type I and II enriched layers in sections 1, 2, and 3 of the embodiment of the present invention. Specific implementation method:

[0057] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0058] The following is a detailed description of the high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells from several aspects. 1. Lithofacies Fluid Probability Analysis (FFP)

[0059] Lithofacies fluid probability analysis is to analyze the intersection of the well elastic parameter curve (P-wave impedance, P-wave velocity ratio) and the lithology curve, and establish the response range and quantitative probability density function of the elastic parameter body obtained by deterministic inversion corresponding to different lithologies. Then, the elastic parameter body is intersected in three dimensions through cloud transformation (describing the correlation between the two parameters), and the P-wave impedance and P-wave velocity ratio body obtained by prestack deterministic inversion is transformed into the lithofacies probability body of each major lithofacies type through Bayesian lithofacies classification and discrimination, that is, the spatial distribution probability of each type of lithofacies.

[0060] Divide the lithology required for geostatistical inversion, that is, define discrete attributes. The lithology types collected on the well are simplified into several categories through crossplot analysis, and these major lithofacies types can be distinguished by elastic attributes (P-wave impedance and Vp / Vs), and the lithology types that cannot be distinguished are merged into one.

[0061] 2. Lithofacies trend constrained inversion

[0062] The lithofacies probability body obtained by lithofacies probability analysis is added as a three-dimensional constraint to the inversion process, and a reliable reservoir spatial distribution is obtained through the Markov chain-Monte Carlo algorithm. The prior information input here includes: the low-frequency model constructed by the stratigraphic framework, the three-dimensional lithofacies probability body, probability density function, and variogram obtained in the previous three steps. Compared with conventional geostatistical methods, the three-dimensional lithofacies probability body is used here to replace the conventional zero-dimensional lithofacies ratio.

[0063] (1) Variogram analysis

[0064] The variogram describes the change of the spatial distribution characteristics of a certain attribute with distance. It is a function of distance, describing the structure and characteristic scale of horizontal and vertical geological characteristics, and it has a constraining effect on the scale of the geological body. According to the actual situation, the vertical range and horizontal range of the inversion are given by the lithofacies. The vertical range is based on the thickness of a certain lithofacies, and the horizontal range is combined with the geological knowledge of the area and the extracted seismic attributes. According to the deterministic inversion results, the continuity and ductility of a lithofacies in the plane are quantitatively determined.

[0065] (2) Determination of probability density function

[0066] The probability density function (PDF) describes the probability distribution of a certain attribute in space, that is, the possibility of the distribution of rock physical parameters corresponding to a specific lithology. It includes the PDF of discrete attributes and the PDF of continuous attributes. The PDF of discrete attributes is to define the lithology ratio. The PDF of discrete attributes is tested by referring to the maximum likelihood lithofacies body histogram of the previous lithofacies probability analysis (FFP). The PDF of continuous attributes is to fit and edit the two-dimensional PDF composed of the intersection diagram of the continuous attributes (impedance, longitudinal and transverse wave velocity ratio or porosity, etc.) on the well. You can refer to the rock physical quantity version of the area and use polygons to clearly identify different lithologies.

[0067] The impedance body obtained by inversion without facies separation often has the following problems: large-scale geological bodies are often discontinuous, so it is impossible to obtain a large set of continuous geological bodies, and the elastic properties obtained by inversion tend to be the average value in the intersection diagram. The joint inversion of elastic properties and lithofacies can effectively solve this problem. The defect of conventional geostatistical methods is that it is difficult to give the lithofacies ratio of the statistical model in areas with few wells, and the inversion results are greatly affected by the stratigraphic framework. For reservoirs with rapid lateral changes, the inversion effect has poor correspondence with the distribution and extension trend of the actual reservoir. Since the lithofacies trend constrained inversion uses deterministic inversion results to perform lithofacies probability analysis, and deterministic inversion is completely in accordance with seismic data, this method takes into account the lateral predictability of deterministic inversion while maintaining the high vertical resolution of geostatistical inversion, combining the advantages of the two and effectively solving the problem of high-precision inversion of interlayer thin shales in areas with few wells.

[0068] 3. Multiple regression analysis

[0069] The multivariate regression analysis method fits the curves of the wellbore elastic parameters (P-wave impedance, S-wave impedance, λρ) and the enriched layer parameters (porosity, TOC, shale content) into a linear regression relationship, and applies the relationship formula to the high-precision prestack elastic parameter body obtained by the prestack geostatistical inversion of the previous three-dimensional lithofacies trend constraint, and converts it into an enriched layer parameter body. For example, the multivariate regression method uses a method of jointly predicting multiple elastic parameters that affect TOC, which is more accurate than a single parameter (only using P-wave impedance).

[0070] Based on the above method, the present invention provides a method for high-precision prediction of parameters of interlayer shale oil-rich layers in a well-polluted area, comprising the following steps:

[0071] S1. After optimizing the pre-stack CRP gathers in the work area, the actual gathers with smaller deviation from the AVO characteristics of the underground medium are obtained and partially stacked to obtain the offset partial stacking volume used for inversion;

[0072] The method for optimizing the pre-stack CRP gathers in the work area is to optimize the gathers by using a set of technical processes including Radon transform denoising, spectrum balancing, and elimination of fitting amplitude compensation. The specific process is:

[0073] The parabola Radon transform filtering technology is used to separate the effective signal from the noise in the transform domain, and then the inverse transform is used to restore the seismic gathers; the separation degree between the primary wave and the noise is determined by the moving average time difference between the primary wave and the noise; after processing, the multiple waves are eliminated, the noise is suppressed, and it is closer to the forward gathers;

[0074] The spectrum balancing technique is used to perform time-frequency analysis on the gathers based on Gabor transform. The near-mid-offset gathers or small-angle gathers that match the AVO characteristics of the well are used as standard gathers. The average time-frequency spectrum of the reference gathers is used as the expectation to compensate for the high-frequency energy of the far-offset gathers.

[0075] The technique of eliminating fitting amplitude compensation is adopted to select the trace with correct AVO law as the standard trace. The data points of amplitude variation with angle of the standard trace are fitted by least square method. The overdetermined equation is solved to obtain the intercept and gradient, and the relationship between amplitude variation with angle is obtained. The relationship between amplitude variation with angle is then extrapolated to the near-offset and far-offset trace gathers to achieve the effect close to the real AVO law and improve the inversion result of P-wave and S-wave velocity ratio.

[0076] S2. Perform well-seismic calibration on the obtained offset stack, extract multi-well sub-waves, establish a well interpolation low-frequency model, and then perform pre-stack deterministic inversion to obtain an elastic parameter body of deterministic inversion; the specific methods include:

[0077] 201) Perform well seismic calibration on the wells in the work area and extract multi-well sub-waves corresponding to each part of the stacked body;

[0078] 202) Using the stratigraphic framework and well curves, a low-frequency model of well difference of P-wave impedance and P-S wave velocity ratio is established for prestack deterministic inversion by well interpolation method;

[0079] 203) Using the low-frequency model of well difference, pre-stack simultaneous inversion is performed on the stacked bodies at different offsets to generate elastic parameter bodies of P-wave impedance and P-S wave velocity ratio;

[0080] 204) The generated elastic parameter bodies of P-wave impedance and P-S wave velocity ratio are used to determine the approximate distribution range of the reservoir and provide parameters such as lateral range and sand-to-ground ratio in the variation function for geostatistical inversion.

[0081] S3, such as Figure 4 As shown, the elastic parameter body obtained by deterministic inversion is used for lithofacies probability analysis, and the final three-dimensional lithofacies probability body is obtained by comprehensive analysis of the well lithology curve; the obtained three-dimensional lithofacies probability body is used to constrain geostatistical inversion, and the elastic parameter body and lithofacies body of the three-dimensional lithofacies trend constrained geostatistical inversion are obtained;

[0082] The method of using the elastic parameter body obtained by deterministic inversion to perform lithofacies probability analysis and combining the well lithology curve comprehensive analysis to obtain the final three-dimensional lithofacies probability body includes the following steps:

[0083] 301) The elastic parameter curves of the wellbore P-wave impedance, P-wave velocity ratio and lithology curves are intersected and analyzed to establish the response range and quantitative probability density function of the elastic parameter body obtained by deterministic inversion corresponding to different lithologies;

[0084] 302) Based on the Bayesian discriminant principle, lithofacies probability analysis (FFP) is performed, and the elastic parameter body is transformed into a three-dimensional lithofacies probability body of each major lithofacies type through cloud transformation (describing the correlation between two parameters) using probability density function.

[0085] The method of using a three-dimensional lithofacies probability body to constrain geostatistical inversion to obtain a pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body comprises the following steps:

[0086] 401) Based on the input three-dimensional lithofacies probability body and seismic data statistics, the respective prior probability density functions are obtained; that is, the logging data, lithofacies probability body, seismic data and variogram are represented by the prior probability density function;

[0087] 402) fusing the above-mentioned individual probability density functions through Bayesian reasoning to obtain a posterior probability distribution function of reservoir parameters under known seismic, well logging and geological information;

[0088] 403) The Markov chain-Monte Carlo inversion algorithm is used to obtain the statistically correct random sample point distribution according to the above-mentioned posterior probability distribution function, and the realizations of each reservoir parameter to be determined are calculated to obtain the pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body.

[0089] S4. Use the multivariate regression analysis method to convert the elastic parameter body obtained by the three-dimensional lithofacies trend constraint geostatistical inversion into the enrichment layer parameter body. The specific method is as follows:

[0090] The multivariate regression analysis method is used to fit the curves of well elastic parameters (P-wave impedance, S-wave impedance, λρ) and enrichment layer parameters (porosity, TOC, shale content, oil saturation, brittleness) into a linear regression relationship.

[0091] The linear regression relationship formula between the wellbore elastic parameters and the enriched layer parameters is applied to the pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body to convert it into an enriched layer parameter body.

[0092] The relationship between the elastic parameter and the enrichment layer parameter can be expressed as follows:

[0093] I p =a*TOC+b*POR+c*V sh (1)

[0094] I s =a′*TOC+b′*POR+c′*V sh (2)

[0095] λρ=a″*TOC+b″*POR+c″*V sh (3)

[0096] Wherein: elastic parameters: Ip is longitudinal wave impedance; Is is shear wave impedance; λρ is the product of compression modulus and density (sensitive to fluid); enrichment layer parameters: TOC is total organic carbon; POR is porosity; Vsh is mud content; a, a′, a″, b, b′, b″, c, c′, c″ are all constant values.

[0097] The elastic parameters affecting TOC include Ip longitudinal wave impedance, Is transverse wave impedance, and λρ compression modulus and density product. Therefore, the elastic parameters of TOC are respectively related to the longitudinal wave impedance, transverse wave impedance, and compression modulus and density product.

[0098] The elastic parameter regression relationship of TOC is expressed as follows:

[0099] TOC=f*I p +g*I s +h*λρ (5)

[0100] Wherein: d, e, f, g, h all represent constant values.

[0101] S5. According to the enrichment layer classification evaluation standard, the obtained enrichment layer parameter body is subjected to enrichment layer classification prediction.

[0102] Example 1

[0103] like Figure 3 As shown, the implementation process of the method of the present invention is explained below by taking the prediction of enrichment layer parameters of the Lianggaoshan Formation interlayer shale reservoir in the Sichuan-Chongqing PC exploration area as an example.

[0104] 1. Research background

[0105] The Yilong-Pingchang block is located in the Jurassic foreland basin of Dabashan. The Lianggaoshan Formation is located in the foreland of the basin. The source rocks are well developed, the strata of the Middle and Upper Jurassic are thick, the source rocks are buried deep, and the degree of thermal evolution of organic matter is high. There is only one exploration well, PA1, in the area. The Lianggaoshan Formation is the last large-scale lake transgression in the Jurassic, and the strata thickness is generally between 150-250m. The facies types include far-shore sand beaches, near-shore sand bars, estuary bars and deep lake phases. Among them, the deep lake phase black, gray-black mudstone and shale have well-developed lamellar structure, occasionally intercalated with silty sand bands, and the well logging curve shows high gamma and low resistivity characteristics. The fan delta front sand body at the bottom of the Liangshang section is developed and covers the shale of the Liangxia section, forming a good reservoir combination. The plane distribution in the Pingchang area is stable, the gas content is good, and the resource potential is great. It is the key target of the next exploration of the oil field. Deep lake mudstone is also an important target for shale oil and gas exploration. Pores are the main reservoir space, and the pore types are diverse. Organic pores are developed in lacustrine shale, and microcracks contribute to the total pores. In addition to organic pores, a large number of inorganic pores are also developed. The clay mineral content is high, and clay mineral pores are mainly developed, as well as mineral intercrystalline pores, grain edge pores and dissolution pores. The physical properties of Lianggaoshan shale are relatively good. The average porosity of mud shale is 4.3%, which is higher than that of tight sandstone. The porosity of the target layer shale of PA1 well is 7%. The post-pressure test of PA1 well has achieved good results, confirming that the oil and gas exploration potential of Lianggaoshan shale in Pingchang area is huge. However, because there is only one well PA1 in the area, the shale single layer is thin, and the lateral lithology changes rapidly, conventional reservoir prediction methods are difficult to meet the requirements of high-precision reservoir prediction. Therefore, it is very necessary to study the seismic prediction of interlayer thin shale reservoirs in the Lianggaoshan sparse well area.

[0106] 2. Optimization of superimposed CRP gathers

[0107] After conventional processing, the prestack CRP gathers input from the PC work area in central Sichuan still have problems such as low signal-to-noise ratio, low resolution, poor frequency consistency of seismic traces at different offsets, weak energy of near and far traces, and poor lateral continuity. These problems lead to a large deviation between the AVO characteristics of the actual gathers and the AVO characteristics of the underground medium. In order to meet the requirements of prestack elastic parameter inversion and restore the AVO characteristics of the gathers, a series of optimization treatments are required for the prestack CRP gathers. To this end, the prestack gathers of the Lianggaoshan Formation in the PC work area are optimized using a set of technical processes including Radon denoising, spectral balancing technology, elimination fitting amplitude compensation technology, and alpha filtering technology. The optimized gathers have been improved in terms of signal-to-noise ratio, resolution, and AVO fidelity. Five partial stacks (near offset, near-mid offset, mid offset, mid-far offset, and far offset) are obtained by partial stacking for subsequent inversion. Compared with the gathers before optimization, the wavelets extracted from the five partial stacks after gather optimization have improved wavelet symmetry, more stable wavelet phases within the seismic frequency band, and fewer wells are excluded when extracting multi-well integrated wavelets.

[0108] 3. Prestack deterministic inversion

[0109] Well seismic calibration is performed on the wells in the work area, and multi-well sub-waves are extracted for each part of the stack. The layer fault and the interval velocity after velocity correction are input, and the low-frequency model of P-wave impedance and P-S wave velocity ratio is established by well interpolation method.

[0110] Before the geostatistical inversion, the pre-stack simultaneous inversion is performed. Under the constraint of the logging data, the stack of five different offsets is inverted simultaneously to generate the P-wave impedance ( Figure 1 b) Data bodies such as the ratio of longitudinal and transverse wave velocities determine the approximate distribution range of the reservoir and provide parameters such as lateral range and sand-to-ground ratio for subsequent geostatistical inversion.

[0111] 4. 3D lithofacies trend constrained prestack geostatistical inversion

[0112] The prestack deterministic inversion results are used to perform lithofacies probability analysis (FFP) based on the Bayesian discriminant principle to obtain a three-dimensional lithofacies probability volume ( Figure 5 ), the probability body is used as a constraint on the lithofacies trend to obtain the prior probability density function. The prior probability density function, together with the seismic, wavelet, and well curve data, is input into the Markov chain-Monte Carlo inversion algorithm to obtain a reasonable reservoir lithology and physical property spatial distribution data body ( Figure 1 a), namely, the prestack inversion elastic parameter volume. In the inversion process, the blind well test method is used to evaluate the reliability of the selection of geostatistical parameters, the predictability of the geostatistical inversion method, and the stability of the prediction method.

[0113] The results of elastic parameter inversion in the PC area are in good agreement with the well curves, and the thin layer resolution is improved from 8m to 3m. Compared with the prestack deterministic inversion results, the resolution is higher and is generally consistent in reflecting the distribution of sand bodies, geological body scale, lithology ratio, etc. (see Figure 1 ). By extracting the plane attributes of various reservoir parameters, it is found that the lateral changes in reservoir parameters are also consistent with the geological laws and sedimentary facies of the area. Compared with the pre-stack conventional geostatistical inversion results, the resolution is comparable, and the inversion effect corresponds better to the distribution and extension trend of the actual reservoir (see Figure 2 ). Therefore, the 3D lithofacies trend-constrained prestack geostatistical inversion is more suitable for situations where wells are sparse and lithology changes rapidly.

[0114] 5. Multiple regression analysis to predict enrichment layer parameters

[0115] Changes in enrichment layer parameters such as TOC, porosity, and clay content will affect multiple elastic parameters such as longitudinal wave impedance and shear wave impedance. The Daqing Sichuan-Chongqing exploration area is in the early stages of exploration, and drilling is sparse. If neural network methods are used to fit the relationship between the two, overfitting is likely to occur, resulting in incorrect regression relationships and lack of predictability. Therefore, the prediction of enrichment layer parameters in areas with few wells uses a linear regression relationship, and the relationship between elastic parameters and enrichment layer parameters can be expressed as follows:

[0116] I p =a*TOC+b*POR+c*V sh (1)

[0117] I s =a′*TOC+b′*POR+c′*V sh (2)

[0118] λρ=a″*TOC+b″*POR+c″*V sh (3)

[0119] Where: Elastic parameters: Ip is the longitudinal wave impedance; Is is the transverse wave impedance; λρ is the product of compression modulus and density (sensitive to fluid);

[0120] Enrichment layer parameters: TOC is total organic carbon; POR is porosity; Vsh is mud content;

[0121] a, a′, a″, b, b′, b″, c, c′, c″ are all constant values.

[0122] Taking TOC prediction as an example, the conventional approach is to use the longitudinal wave impedance I p This single parameter is used to fit the regression relationship:

[0123] TOC=d*I p +e (4)

[0124] Where: d and e are constant values.

[0125] Longitudinal wave impedance I p It is not only related to TOC, but the above formula only uses a single parameter, longitudinal wave impedance, to predict TOC, so the fitting accuracy is bound to be low. For this reason, we use a method of jointly predicting multiple elastic parameters that affect TOC to improve the prediction accuracy. The regression relationship is as follows:

[0126] TOC=f*I p +g*I s +h*λρ (5)

[0127] Wherein: d, e, f, g, h all represent constant values.

[0128] The TOC predictions of the multi-parameter and single-parameter methods were compared with the measured curves on the well. Figure 6 ), the correlation coefficient between the predicted TOC and the measured TOC by the multi-parameter prediction method increased from 0.72 to 0.78. From the comparison of TOC prediction profiles ( Figure 7 ), the TOC predicted by the multivariate regression of elastic parameters in the elliptical frame is in better agreement with the well curve, and it can be seen that the latter has a higher fitting accuracy. Based on this principle, a multivariate regression relationship is established between prestack elastic parameters (P-wave impedance, S-wave impedance, λρ) and enriched layer parameters (porosity, TOC, brittleness) on the well, and the high-precision prestack elastic parameter body obtained by the prestack geostatistical inversion constrained by the previous three-dimensional lithofacies trend is converted into an enriched layer parameter body using the relationship formula, as follows:

[0129] (1) Total organic carbon content (TOC): Based on the high-precision prestack elastic parameter inversion volume, the longitudinal wave impedance (I p ), shear wave impedance (I s ), λρ and TOC multivariate regression relationship, to achieve TOC prediction combined with well seismic, compared with the traditional longitudinal wave impedance conversion TOC method, with higher accuracy ( Figure 7 ). The formula is as follows:

[0130] I S =I P / (V P / V S ) (6)

[0131]

[0132] TOC=6.6620-0.00175388*I P / 1000+1.60689e -006 *I s +6.73462e -008 *λρ (8)

[0133] In the formula: Vp / Vs is the ratio of longitudinal and transverse wave velocities.

[0134] (2) Porosity (POR): Based on the high-precision prestack elastic parameter inversion volume, the multivariate regression relationship between P-wave impedance, S-wave impedance, λρ and porosity is used to achieve the porosity prediction by combining well-seismic data. Compared with the conventional method, the correlation coefficient between prediction and measurement is improved from 0.71 to 0.79 ( Figure 8 ), the porosity predicted by the multivariate regression of elastic parameters at the black arrow on the predicted section is more consistent with that on the well ( Fig. 9 ), with higher precision. The formula is as follows:

[0135] POR=11.4748-0.00130691*I S +0.000156502*I P -1.11131e -008 *λρ (9)

[0136] (3) Brittleness index (BRT): The high-precision prestack elastic parameter inversion data volume obtained by inversion uses the elastic parameter method to first calculate the Young's modulus (E) and Poisson's ratio (σ), and then calculate the brittleness index ( Fig.10 ). The formula is as follows:

[0137] σ=[(V P / V S ) 2 -2] / [2(V P / V S ) 2 -2] (10)

[0138]

[0139] BRT σ =(0.4-σ) / (0.14-0.15)*100 (12)

[0140]

[0141] BRT=(BRT e +BRT σ ) / 2 (14)

[0142] Where: σ is Poisson's ratio; Vp / Vs is the ratio of longitudinal and transverse wave velocities; E is Young's modulus; BRT is brittleness; BRT e — Brittleness calculated from Young's modulus; BRT σ — Brittleness of Poisson’s ratio calculation;

[0143] The overall shale brittleness index is above 40%, and fracturing can easily form a network of fractures, realizing the seismic prediction of the brittleness index under well-seismic constraints.

[0144] (4) Mud content (Vsh): The multivariate regression relationship between longitudinal wave impedance, transverse wave impedance, λρ and mud content is used to obtain the formula for mud content as follows:

[0145] Vsh=0.965175-1.47635e -05 *I p -5.29588e -05 *I / -1.36267e -09 *λρ (15)

[0146] 6. Quantitative prediction of enrichment layer

[0147] Table 1 is the evaluation standard table of the shale enrichment layer of Lianggaoshan Formation in the Sichuan-Chongqing exploration area. According to the classification and evaluation standard of the enrichment layer, the distribution prediction of various types of shale oil and gas enrichment layers in Lianggaoshan Formation in the Pingchang 3D area is realized by using the predicted enrichment layer main control parameters (see Fig.11 and Fig.12 ), the well compliance rate of enriched layers above 3 meters is as high as 76.5%, which can provide effective support for well site deployment.

[0148] Table 1

[0149]

[0150] The method of the present invention mainly obtains the lithofacies probability body after lithofacies probability analysis based on deterministic inversion, thereby constraining the prestack geostatistical inversion technology to obtain an inversion body that is more in line with geological characteristics, achieving the purpose of taking into account the vertical resolution accuracy and lateral extension trend of the reservoir. Then, the multi-parameter regression method is used to predict the geological parameters of the enriched layer such as total porosity and total organic carbon content to meet the needs of predicting interlayer shale oil sweet spots in the Sichuan and Chongqing areas with few wells.

Claims

1. A high-precision prediction method for parameters of interlayer shale oil-rich layers in a few-well area, characterized by: The following steps are involved: S1. After optimizing the pre-stack CRP gathers in the work area, the actual gathers with smaller deviation from the AVO characteristics of the underground medium are obtained and partially stacked to obtain the offset partial stacking volume used for inversion; S2. Perform well-seismic calibration on the obtained offset stack, extract multi-well sub-waves, establish a well interpolation low-frequency model, and then perform pre-stack deterministic inversion to obtain an elastic parameter body of deterministic inversion; S3. Use the elastic parameter body obtained by deterministic inversion to perform lithofacies probability analysis, and combine it with the well lithology curve for comprehensive analysis to obtain the final three-dimensional lithofacies probability body; use the obtained three-dimensional lithofacies probability body to constrain geostatistical inversion, and obtain the elastic parameter body of three-dimensional lithofacies trend constraining geostatistical inversion; S4. Using multiple regression analysis, the elastic parameter volume obtained by geostatistical inversion of three-dimensional lithofacies trend constraints is converted into an enrichment layer parameter volume; S5. According to the enrichment layer classification evaluation standard, the obtained enrichment layer parameter body is subjected to enrichment layer classification prediction.

2. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 1 is characterized by: Step S1: Optimize the pre-stack CRP gathers in the work area by using Radon transform denoising, spectrum balance, and elimination fitting amplitude compensation as a set of technical processes to optimize the gathers; The specific method of the optimization process is: The parabola Radon transform filtering technology is used to separate the effective signal from the noise in the transform domain, and then the inverse transform is used to restore the seismic gathers; the separation degree between the primary wave and the noise is determined by the moving average time difference between the primary wave and the noise; after processing, the multiple waves are eliminated, the noise is suppressed, and it is closer to the forward gathers; The spectrum balancing technique is used to perform time-frequency analysis on the gathers based on Gabor transform. The near-mid-offset gathers or small-angle gathers that match the AVO characteristics of the well are used as standard gathers. The average time-frequency spectrum of the reference gathers is used as the expectation to compensate for the high-frequency energy of the far-offset gathers. The technique of eliminating fitting amplitude compensation is adopted to select the trace with correct AVO law as the standard trace. The data points of amplitude variation with angle of the standard trace are fitted by least square method. The overdetermined equation is solved to obtain the intercept and gradient, and the relationship between amplitude variation with angle is obtained. The relationship between amplitude variation with angle is then extrapolated to the near-offset and far-offset trace gathers to achieve the effect close to the real AVO law and improve the inversion result of P-wave and S-wave velocity ratio.

3. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 1 is characterized by: Step S2 is a method for performing prestack deterministic inversion to obtain an elastic parameter volume, comprising: 201) Perform well seismic calibration on the wells in the work area and extract multi-well sub-waves corresponding to each part of the stacked body; 202) Using the stratigraphic framework and well curves, a well interpolation method is used to establish a low-frequency model of well difference of P-wave impedance and P-S wave velocity ratio for prestack deterministic inversion; 203) Using the low-frequency model of well difference, pre-stack simultaneous inversion is performed on the stacked bodies at different offsets to generate elastic parameter bodies of P-wave impedance and P-S wave velocity ratio; 204) Using the generated elastic parameter bodies of P-wave impedance and P-S wave velocity ratio, the approximate distribution range of the reservoir is determined, and the parameters such as lateral range and sand-to-ground ratio in the variation function are provided for geostatistical inversion.

4. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 1 is characterized by: S3 uses the elastic parameter body obtained by deterministic inversion to perform lithofacies probability analysis, and combines the well lithology curve for comprehensive analysis to obtain the final three-dimensional lithofacies probability body method, including: 301) Intersection analysis of elastic parameter curves of P-wave impedance, P-wave velocity ratio and lithology curves on the well, and establishment of response range and quantitative probability density function of elastic parameter bodies obtained by deterministic inversion corresponding to different lithologies; 302) Based on the Bayesian discriminant principle, the lithofacies probability analysis is performed, and the elastic parameter body is transformed into a three-dimensional lithofacies probability body of each major lithofacies type through cloud transformation using the probability density function.

5. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 1 is characterized by: Step S3 uses a three-dimensional lithofacies probability body to constrain geostatistical inversion, and obtains a pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body, comprising the following steps: 401) Obtaining respective prior probability density functions based on the input three-dimensional lithofacies probability volume and seismic data statistics; 402) fusing the above-mentioned individual probability density functions through Bayesian reasoning to obtain a posterior probability distribution function of reservoir parameters under known seismic, well logging and geological information; 403) Using the Markov chain-Monte Carlo inversion algorithm, a statistically correct random sample point distribution is obtained according to the above-mentioned posterior probability distribution function, and each realization of each reservoir parameter to be determined is calculated to obtain a pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body.

6. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 5 is characterized by: Step S4 is a method for converting the elastic parameter body obtained by geostatistical inversion of three-dimensional lithofacies trend constraints into an enrichment layer parameter body by using a multivariate regression analysis method, comprising: The multivariate regression analysis method is used to fit the curves of wellbore elastic parameters and enrichment layer parameters into a linear regression relationship. The linear regression relationship formula between the wellbore elastic parameters and the enriched layer parameters is applied to the pre-stack three-dimensional lithofacies trend constrained geostatistical inversion elastic parameter body to convert it into an enriched layer parameter body.

7. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 6 is characterized by: The relationship between the elastic parameter and the enrichment layer parameter can be expressed as follows: Yo p =a*TOC+b*FOR+c*V sh (1) Yo s =a′*TOC+b′*FOR+c′*V sh (2) λρ=a″*TOC+b″*POR+c″*V sh (3) Where: Ip is the longitudinal wave impedance; Is is the shear wave impedance; λρ is the product of the compression modulus and density (sensitive to fluid); TOC is the total organic carbon; POR is the porosity; Vsh is the mud content; a, a′, a″, b, b′, b″, c, c′, c″ are all constants.

8. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 6 is characterized by: The above-hole elastic parameters include longitudinal wave impedance, shear wave impedance and λρ.

9. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 6 is characterized by: The enrichment layer parameters include porosity, TOC and mud content.

10. The high-precision prediction method for parameters of interlayer shale oil-rich layers in areas with few wells according to claim 9 is characterized by: The elastic parameters affecting TOC include Ip for longitudinal wave impedance, Is for transverse wave impedance; λρ for the product of compression modulus and density, so the elastic parameters of TOC are respectively related to the longitudinal wave impedance, transverse wave impedance, and the product of compression modulus and density; The elastic parameter regression relationship of TOC is expressed as follows: TOC=f*I p +g*I s +h*λρ(5)wherein d, e, f, g, and h all represent constant values.

Images