Shale facies combination characterization and analysis method based on time domain division of spectral analysis

Abstract

This invention discloses a method for characterizing and analyzing shale lithofacies assemblages based on time-domain spectral analysis. The specific steps include: (1) correcting the natural gamma curve of a vertical well to true vertical depth, obtaining the corrected natural gamma curve; (2) performing time-series analysis on the natural gamma curve to obtain the time-domain eccentricity cyclic curve and sedimentation rate values ​​at different depths; (3) converting the depth-domain shale lithofacies assemblages into time-domain shale lithofacies assemblages based on the sedimentation rate values, and classifying the time-domain shale lithofacies assemblages based on the time-domain eccentricity cyclic curve to obtain the classification results of the time-domain shale lithofacies; (4) obtaining and analyzing the classification results of the time-domain shale lithofacies from multiple vertical wells to obtain the lithofacies assemblages variation trend within the region. This invention provides a detailed characterization and analysis of lithofacies assemblages under a time-domain framework, clearly defining the lateral lithofacies assemblages variation trend within the region.

Description

Technical Field

[0001] This invention relates to the field of petroleum and natural gas geology, and in particular to a method for characterizing and analyzing shale lithofacies assemblages based on time-domain division of spectral analysis. Background Technology

[0002] The study of shale lithofacies is fundamental and a crucial component of shale geological characteristics. Currently, the classification of shale lithofacies is becoming increasingly refined, facilitating the identification of favorable lithofacies assemblages and their distribution, and providing important geological evidence for shale gas exploration and evaluation. However, existing technologies primarily analyze shale lithofacies within a depth framework, neglecting the detailed characterization of shale lithofacies assemblages within a time framework. This fails to adequately reflect the trends in regional shale lithofacies assemblages, hindering the analysis of some of the reasons for vertical differences in regional lithofacies assemblages, resulting in incomplete analysis and unsatisfactory results.

[0003] Furthermore, extensive follow-up studies of shale work areas have revealed differences in lithofacies assemblages within the shale formations, and these assemblages exhibit certain trends of change based on various factors. Traditional analyses of shale lithofacies under deep frameworks are insufficient to explore these trends or clarify the reasons for these changes, thus failing to provide detailed and reliable evidence for more refined 3D modeling of the work areas. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies that mainly focus on the analysis of lithofacies under a deep framework and cannot clearly identify the trend of changes in lateral shale lithofacies assemblages. This invention provides a method for characterizing and analyzing shale lithofacies assemblages based on time domain division of spectral analysis. It provides a detailed characterization of lithofacies assemblages under a time framework and can explore the trend changes in shale lithofacies assemblages.

[0005] To achieve the above-mentioned objectives, the present invention provides the following technical solution:

[0006] A method for characterizing and analyzing shale lithofacies assemblages based on time-domain spectral analysis includes the following steps:

[0007] S1, Based on well logging curve data and electrical imaging data, perform true vertical depth correction on the natural gamma curve of the vertical well to obtain the true vertical depth corrected natural gamma curve;

[0008] S2, perform time series analysis on the natural gamma curve to obtain the eccentricity cycle curve in the depth domain, convert the eccentricity cycle curve in the depth domain into the eccentricity cycle curve in the time domain, and obtain the deposition rate values ​​at different depths based on the eccentricity cycle curve in the time domain.

[0009] S3, obtain the classification results of shale facies in the depth domain, combine the sedimentation rate value to transform the shale facies assemblage in the depth domain into a shale facies assemblage in the time domain, and divide the shale facies assemblage in the time domain based on the eccentricity cyclic curve in the time domain to obtain the classification results of shale facies in the time domain.

[0010] S4. Obtain and analyze the classification results of the time-domain shale facies of multiple vertical wells to obtain the trend of facies combination changes in the region.

[0011] As a preferred embodiment of the present invention, the shale lithofacies assemblage characterization and analysis method based on time-domain division of spectral analysis, step S1 specifically includes the following steps:

[0012] S11, Obtain logging curve data and electrical imaging data from the vertical well, and perform depth correction on the logging curve data and electrical imaging data;

[0013] S12, Based on the depth-corrected electrical imaging logging data, obtain dynamic and static electrical imaging images, and pick up and record the stratigraphic structure features within the electrical imaging images;

[0014] S13, collect the formation bedding structure features and the logging curve data of the depth-corrected natural gamma data, wellbore inclination and azimuth data into a dataset;

[0015] S14. Based on the dataset, perform true vertical depth correction on the natural gamma curve of the vertical well to eliminate the influence of well inclination and formation dip on the shape of the natural gamma curve, and obtain the natural gamma curve after true vertical depth correction.

[0016] As a preferred embodiment of the present invention, the method for characterizing and analyzing shale facies assemblages based on time-domain spectral analysis, wherein the method for picking up and recording the stratigraphic bedding structure of dynamic electrical imaging is as follows: at least one point of stratigraphic bedding structure is picked up and recorded within every 1-meter depth of the dynamic electrical imaging image.

[0017] As a preferred embodiment of the present invention, the method for characterizing and analyzing shale lithofacies assemblages based on time-domain division of spectral analysis, specifically includes the following steps in S2:

[0018] S21, Perform detrending analysis on the natural gamma curve after true vertical depth correction to obtain the trend line of the natural gamma curve after true vertical depth correction;

[0019] S22, perform spectral analysis on the trend line to obtain a spectrum diagram, obtain high confidence frequencies based on the spectrum diagram, and convert the spectrum diagram to a logarithmic scale thickness spectrum diagram to obtain the deposition cycle thickness corresponding to the high confidence frequencies;

[0020] S23. By conducting research and analysis, the stable eccentricity period of the formation in the area where the vertical well is located is obtained. The ratio of the eccentricity period is compared with the ratio of the depositional cycle thickness corresponding to the peak power of the trend line to obtain the depositional cycle thickness corresponding to the eccentricity period.

[0021] S24, calculate the Gaussian filter window length based on the deposition cycle thickness corresponding to the eccentricity period, and perform Gaussian filtering on the natural gamma curve after true vertical depth correction based on the window length to obtain the eccentricity cycle curve in the depth domain.

[0022] S25, the eccentricity cycle curve in the depth domain is processed according to the tuning method in the study of cyclic stratigraphy to obtain the astronomical time scale. Based on the astronomical time scale, the eccentricity cycle curve in the depth domain is transformed into the eccentricity cycle curve in the time domain. At the same time, the sedimentation rate values ​​at different depths are calculated and recorded during the transformation process.

[0023] As a preferred embodiment of the present invention, the method for characterizing and analyzing shale facies assemblages based on time-domain spectral analysis includes the following specific steps for calculating the deposition rate value: by calculating the ratio of the top-to-bottom difference of the eccentricity cycle curve in the depth domain of the vertical well to the top-to-bottom difference of the eccentricity cycle curve in the time domain, the ratio is obtained as the deposition rate value.

[0024] As a preferred embodiment of the present invention, the method for characterizing and analyzing shale lithofacies assemblages based on time-domain division of spectral analysis, specifically includes the following steps in S3:

[0025] S31, obtain Lithoscanner data from the vertical well, perform depth correction on the Lithoscanner data, summarize and classify the depth-corrected Lithoscanner data. The Lithoscanner data contains minerals such as quartz, feldspar, calcite, dolomite and clay minerals. Based on the content of the minerals and the differences in the content of minerals in different shale lithofacies types in the three-end-member classification scheme, the lithofacies of the shale are divided to obtain the depth domain shale lithofacies classification results.

[0026] S32, calculate the ratio of the depth data in the depth domain shale facies classification result to the deposition rate under the depth data, and obtain the time data of each shale facies in the depth domain shale facies classification result, and reassemble the time domain shale facies assemblage from the time data;

[0027] S33, the lithofacies of the time domain are divided based on the variation period of the eccentricity cycle curve of the time domain to obtain the secondary lithofacies assemblage;

[0028] S34, the secondary lithofacies assemblages divided by the eccentricity cyclic curves in the same time domain are combined into a primary lithofacies assemblage, and the primary lithofacies assemblages and the secondary lithofacies assemblages together construct a shale lithofacies assemblage in a single vertical well time domain.

[0029] As a preferred embodiment of the present invention, the shale lithofacies assemblage characterization and analysis method based on spectral analysis time domain division includes the following specific steps in S33: dividing the shale lithofacies assemblage in the time domain into independent units of the specified number of cycles, i.e., secondary lithofacies assemblages, according to the number of cycles of the eccentricity cyclic curve in the time domain, and naming the secondary lithofacies assemblage based on the lithofacies with the largest proportion in the secondary lithofacies assemblage.

[0030] As a preferred embodiment of the present invention, the shale lithofacies assemblage characterization and analysis method based on spectral analysis time domain division, step S4 specifically includes: performing steps S1-S3 on multiple vertical wells located at different locations in a region to obtain the time domain shale lithofacies classification results of the multiple vertical wells, judging the lithofacies assemblage change trend among the multiple vertical wells based on the time domain shale lithofacies classification results of the multiple vertical wells, and then combining the geographical location of the multiple vertical wells to obtain the lithofacies assemblage change trend within the region.

[0031] As a preferred embodiment of the present invention, the method for characterizing and analyzing shale lithofacies combinations based on time-domain spectral analysis includes the following specific steps for determining the lithofacies combination change trend among multiple vertical wells based on the time-domain shale lithofacies combinations: obtaining the proportion of the same secondary lithofacies combination in the multiple vertical wells based on the time-domain shale lithofacies combinations, and determining the change in the proportion, which is the lithofacies combination change trend among the multiple vertical wells.

[0032] Based on the same concept, a shale lithofacies assemblage characterization and analysis device based on spectral analysis time domain division is also proposed, including at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, which are executed by the at least one processor to enable the at least one processor to execute the shale lithofacies assemblage characterization and analysis method based on spectral analysis time domain division described above.

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

[0034] This invention transforms depth-domain data into time-domain data based on time-series analysis, and finally divides shale lithofacies assemblages based on the time domain, achieving a fine characterization of shale lithofacies assemblages under a time framework, which well reflects the trend of shale lithofacies assemblages changing within the region. Attached image description:

[0035] Figure 1 This is a schematic diagram of the shale lithofacies assemblage characterization and analysis method based on time domain division of spectral analysis in Example 1;

[0036] Figure 2 This is a detrended graph of the natural gamma curve after true vertical depth correction in Example 1.

[0037] Figure 3 The MTM spectrum obtained from the trend line spectrum analysis in Example 1;

[0038] Figure 4 This is a diagram of the three-terminal classification scheme in Example 1;

[0039] Figure 5 This is a classification result diagram of the shale lithofacies assemblage in the depth domain of well 1A in Example 1;

[0040] Figure 6 This is a classification result diagram of the time-domain shale lithofacies assemblage obtained from the X405 eccentricity cycle curve in the time domain of well 1A in Example 1A;

[0041] Figure 7 This is a combination diagram of the classification results of time-domain shale lithofacies assemblages obtained from the time-domain eccentricity cyclic curve X405 of wells A, B, C and D in Example 1. Detailed Implementation

[0042] The present invention will be further described in detail below with reference to experimental examples and specific embodiments. However, this should not be construed as limiting the scope of the above-mentioned subject matter of the present invention to the following embodiments; all technologies implemented based on the content of the present invention fall within the scope of the present invention.

[0043] Example 1

[0044] like Figure 1 The diagram shows a flowchart of a shale lithofacies assemblage characterization and analysis method based on time-domain partitioning of spectral analysis. The specific steps are as follows:

[0045] S1, Based on well logging curve data and electrical imaging data, perform true vertical depth correction on the natural gamma curve of the vertical well to obtain the true vertical depth corrected natural gamma curve;

[0046] Specifically, S1 includes the following steps:

[0047] S11, Obtain logging curve data and electrical imaging data from the vertical well, and perform depth correction on the logging curve data and electrical imaging data;

[0048] Specifically, well logging data, electrical imaging data, Lithoscanner data, and XRD data are obtained from the vertical well. The natural gamma curve of the conventional well logging data is used as the standard curve. The natural gamma curves of the electrical imaging data and Lithoscanner data are moved to have the same characteristics at the same depth as the standard curve. Then, the XRD data points are projected onto the Lithoscanner data curve. The XRD data depth is moved until it has the same trend as the Lithoscanner data, thereby determining the accuracy of the Lithoscanner data and completing the depth correction of the conventional well logging data and the electrical imaging data.

[0049] S12, Based on the depth-corrected electrical imaging logging data, obtain dynamic and static electrical imaging images, and pick up and record the stratigraphic structure features within the electrical imaging images;

[0050] Specifically, since the electro-imaging processing results also contain multiple types of data, such as low-resistivity fractures, high-resistivity fractures, and induced fractures, it is necessary to extract the bedding fractures separately to avoid the influence of other non-bedding signals. The bedding fractures contain stratigraphic bedding structure features.

[0051] Furthermore, to reduce errors in picking and recording stratigraphic bedding features in dynamic electro-imaging images and avoid error accumulation that leads to unsatisfactory subsequent processing results, this embodiment picks and records the stratigraphic bedding structure at least one point within every 1-meter depth of the dynamic electro-imaging image. This picking and recording lays a good foundation for the subsequent true vertical depth correction of the natural gamma curve of the vertical well, thereby making the results obtained from subsequent time series analysis more accurate and improving the lithofacies classification accuracy of this embodiment.

[0052] S13, collect the formation bedding structure features and the logging curve data of the depth-corrected natural gamma data, wellbore inclination and azimuth data into a dataset;

[0053] S14. Based on the dataset, perform true vertical depth correction on the natural gamma curve of the vertical well to eliminate the influence of well inclination and formation dip on the shape of the natural gamma curve, and obtain the natural gamma curve after true vertical depth correction.

[0054] S2, perform time series analysis on the natural gamma curve to obtain the eccentricity cycle curve in the depth domain, convert the eccentricity cycle curve in the depth domain into the eccentricity cycle curve in the time domain, and obtain the deposition rate values ​​at different depths based on the eccentricity cycle curve in the time domain.

[0055] S21, Perform detrending analysis on the natural gamma curve after true vertical depth correction to obtain the trend line of the natural gamma curve after true vertical depth correction;

[0056] Specifically, four detrending methods—LOWESS locally weighted regression scatter smoothing, rLOWESS robust locally weighted regression scatter smoothing, LOESS locally weighted scatter regression, and rLOESS robust locally weighted scatter regression—were used to perform detrending analysis on the natural gamma curve of the single vertical well after true vertical depth correction, resulting in the following... Figure 2 The results show that the LOESS curve obtained by using the LOESS local weighted scatter regression method is closest to the natural gamma curve of a single vertical well after true vertical depth correction. Therefore, the LOESS curve is used as the trend line of the natural gamma curve of a single vertical well after true vertical depth correction for subsequent processing.

[0057] S22, perform spectral analysis on the trend line to obtain a spectrum diagram, obtain high confidence frequencies based on the spectrum diagram, and convert the spectrum diagram to a logarithmic scale thickness spectrum diagram to obtain the deposition cycle thickness corresponding to the high confidence frequencies;

[0058] Specifically, the LOESS curve is processed using the Fast Fourier Transform provided by Matlab to obtain the following results: Figure 3 The MTM spectrum shown records high-confidence frequencies with a confidence level exceeding 95%, and the high-confidence deposition cycle thickness corresponding to these high-confidence frequencies is obtained through logarithmic scaling. Figure 3 The deposition cycle thicknesses corresponding to the peak power frequencies with a confidence level exceeding 95% are 11.40m, 3.65m, 3.28m, 2.99m, 2.74m, 2.27m, 2.12m, 1.69m, and 1.62m, respectively.

[0059] S23. By conducting research and analysis, the stable eccentricity period of the formation in the area where the vertical well is located is obtained. The ratio of the eccentricity period is compared with the ratio of the depositional cycle thickness corresponding to the peak power of the trend line to obtain the depositional cycle thickness corresponding to the eccentricity period.

[0060] Specifically, the stable eccentricity period of the strata was obtained through surveys, such as the stable long eccentricity period of 405 kyr and the short eccentricity periods of 125 kyr and 95 kyr in the Wufeng Formation-Longyi Formation. The ratio of the sedimentary cycle thickness of 11.40:3.65:2.74 is closest to the ratio of the eccentricity period of 405:125:95 = 4.26:1.32:1. Therefore, 11.40 m, 3.65 m, and 2.74 m are determined to be the sedimentary cycle thicknesses of the 405 kyr long eccentricity period, the 125 kyr short eccentricity period, and the 95 kyr short eccentricity period, respectively, and are denoted as S405, S125, and S95. In this embodiment, S405 is selected for subsequent research.

[0061] S24, calculate the Gaussian filter window length based on the deposition cycle thickness corresponding to the eccentricity period, and perform Gaussian filtering on the natural gamma curve after true vertical depth correction based on the window length to obtain the eccentricity cycle curve in the depth domain.

[0062] Specifically, the reciprocal of S405 is calculated using numerical values, and the Gaussian filter window length is set based on 80% and 120% of the reciprocal, denoted as C405. Gaussian filtering is applied to the trend line within the C405 window length to obtain the Gaussian-filtered curve D405, which is the long eccentricity cyclic curve in the depth domain.

[0063] Furthermore, after calculation, C405 is 0.070175-0.105263, therefore the Gaussian filter window length corresponding to the 405kyr long eccentricity period is 0.070175-0.105263.

[0064] S25, the eccentricity cycle curve in the depth domain is processed according to the tuning method in the study of cyclic stratigraphy to obtain the astronomical time scale. Based on the astronomical time scale, the eccentricity cycle curve in the depth domain is transformed into the eccentricity cycle curve in the time domain. At the same time, the sedimentation rate values ​​at different depths are calculated and recorded during the transformation process.

[0065] Specifically, based on the tuning method in cyclic stratigraphy research, an astronomical time scale corresponding to a period of 405kyr is established based on the D405 curve, and the long eccentricity cyclic curve in the depth domain is transformed into a long eccentricity cyclic curve in the time domain using the scale, denoted as X405.

[0066] Furthermore, the calculation of the deposition rate value specifically includes: setting the top and bottom difference of the long eccentricity cycle curve X405 of the time domain for different formation depths in well A as T, and using the formation thickness H / deposition time T = deposition rate value V, the deposition rate value V is obtained.

[0067] S3, obtain the classification results of shale facies in the depth domain, combine the sedimentation rate value to transform the shale facies assemblage in the depth domain into a shale facies assemblage in the time domain, and divide the shale facies assemblage in the time domain based on the eccentricity cyclic curve in the time domain to obtain the classification results of shale facies in the time domain.

[0068] Specifically, the steps in S3 include:

[0069] S31, obtain Lithoscanner data from the vertical well, perform depth correction on the Lithoscanner data, summarize and classify the depth-corrected Lithoscanner data. The Lithoscanner data contains minerals such as quartz, feldspar, calcite, dolomite and clay minerals. Based on the content of the minerals and the differences in the content of minerals in different shale lithofacies types in the three-end-member classification scheme, the lithofacies of the shale are divided to obtain the depth domain shale lithofacies classification results.

[0070] Specifically, the Lithoscanner data is first depth-corrected, as described in S11. The depth-corrected Lithoscanner special logging data is then summarized and categorized. The Lithoscanner data contains minerals including quartz, feldspar, calcite, dolomite, and clay minerals. Specifically, silica = quartz + feldspar, and carbonate minerals = calcite + dolomite. The mineral composition is determined based on the content of silica, carbonate, and clay minerals, as follows: Figure 4 The three-terminal-member classification scheme shown divides shale lithofacies. For example, if the mineral content at a certain depth shows calcareous <25%, clayous <25%, and silicatic >75%, then the strata at that depth are identified as siliceous rocks. If calcareous <25%, clayous <25%, and 50% < silicatic <75%, then it is identified as S-2 siliceous shale. After the classification is completed, the following results are obtained: Figure 5 The depth domain shale facies classification results shown are based primarily on Lithoscanner logging data, which is highly reliable, accurate, and continuous, facilitating detailed modeling and description of the region.

[0071] S32, calculate the ratio of the depth data in the depth domain shale facies classification result to the deposition rate under the depth data, and obtain the time data of each shale facies in the depth domain shale facies classification result, and reassemble the time domain shale facies assemblage from the time data;

[0072] Specifically, by dividing the thickness of each shale facies assemblage in the depth domain shale facies classification results by the corresponding sedimentation rate value V at that depth, the temporal data of that shale facies assemblage can be obtained. The temporal data is then used to obtain the result from the curve X405. Figure 6 The time-domain shale facies of well A is shown.

[0073] S33, the lithofacies of the time domain are divided based on the variation period of the eccentricity cycle curve of the time domain to obtain the secondary lithofacies assemblage;

[0074] Specifically, such as Figure 6The X405 curve shown has 10 periods. Based on these 10 periods, the lithofacies are divided into 10 independent units, i.e., second-order lithofacies assemblages. It can be clearly seen that the lithofacies within each second-order lithofacies assemblage are relatively stable. Generally, a certain lithofacies in a second-order lithofacies assemblage has an absolute dominance, such as... Figure 6 The main M-3 type secondary lithofacies assemblage shown is mainly composed of M-3 mixed shale, accounting for more than 93%.

[0075] S34, the secondary lithofacies assemblages divided by the eccentricity cyclic curves in the same time domain are combined into a primary lithofacies assemblage. The primary lithofacies assemblages and the secondary lithofacies assemblages together construct a shale lithofacies assemblage in a single vertical well time domain. The multi-level lithofacies assemblage division scheme not only makes the division results more refined, but also makes it more conducive to judging the trend of shale lithofacies changes in the region.

[0076] Specifically, Figure 6 The 10 secondary lithofacies assemblages shown form a primary lithofacies assemblage unit corresponding to the Wufeng Formation-Longyi Formation. Analysis of the primary lithofacies assemblages reveals that, from bottom to top, they contain mixed, main M-2, main M-3, and main CM-1 secondary lithofacies assemblages. Figure 5 It can be seen that the superposition relationship of the secondary lithofacies assemblages within the primary lithofacies assemblages is that the mixed secondary lithofacies assemblages are located at the bottom, followed by the development of the main M-2 type, then two main M-3 types, and then one main CM-1 type, two main M-2 types, and three main CM-1 types are superimposed on top of it.

[0077] S4. Obtain and analyze the classification results of the time-domain shale facies of multiple vertical wells to obtain the trend of facies combination changes in the region.

[0078] Specifically, S4 includes the following steps:

[0079] Perform steps S1-S3 on the other three vertical wells B, C, and D in region X from west to east, and obtain the following results: Figure 7 The time-domain shale lithofacies classification results of the four vertical wells A, B, C, and D shown are as follows: Figure 7 The secondary lithofacies assemblages of the four vertical wells shown differ in their vertical stacking relationships. Specifically, the bottom of all four wells has a mixed type, while wells A, B, and D in the next cycle all have a main M-2 type assemblages, but well C does not have this secondary lithofacies assemblages. In the next cycle, wells A, B, and C all have two main M-3 type secondary lithofacies assemblages, while well D has only one main M-3 type secondary lithofacies assemblages.

[0080] Furthermore, the proportions of each shale type within the secondary lithofacies assemblages corresponding to different wells in the four vertical wells A, B, C, and D are also different. The lateral comparability can be further refined by analyzing the proportions of each lithofacies within the secondary lithofacies assemblages. For example, after analyzing the lithofacies proportions of the main M-3 type lithofacies assemblages, it can be found that the proportion of M-3 mixed shale within the M-3 type lithofacies assemblages from well A to well C continuously decreases. Therefore, the decrease in the number of main M-3 type assemblages from well A to well D is also reflected in the number of secondary lithofacies assemblages.

[0081] Since the geographical locations of wells A to D are generally oriented east-west, with well A located in the western part of region X and well D located in the eastern part of region X, the analysis of the single-well time-domain shale lithofacies classification results of the above four vertical wells yields the lithofacies combination variation trend among the four wells. This leads to the lithofacies combination variation trend in region X, including the length of the combined depositional period and the corresponding lithofacies combination type and thickness variation trend.

[0082] The method involved in this invention can be effectively applied not only in shale, but also in sandstone and mudstone formations with frequent interbedded layers, which is of great significance for subsequent tight sandstone development and marine-continental shale oil and gas reservoir analysis.

[0083] 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 within the protection scope of the present invention.

Claims

1. A shale facies combination characterization and analysis method based on time domain division of spectral analysis, characterized in that, The method comprises the following steps: S1, true vertical depth correction is performed on a natural gamma ray curve of a vertical well based on logging curve data and electrical imaging data, and a true vertical depth corrected natural gamma ray curve is obtained; S2, time series analysis is performed on the natural gamma ray curve, an eccentricity cycle curve in the depth domain is obtained, the eccentricity cycle curve in the depth domain is converted into an eccentricity cycle curve in the time domain, and a sedimentation rate value at different depths is obtained based on the eccentricity cycle curve in the time domain; S3, a classification result of shale facies in the depth domain is obtained, the shale facies combination in the depth domain is converted into a shale facies combination in the time domain in combination with the sedimentation rate value, the shale facies combination in the time domain is divided based on the eccentricity cycle curve in the time domain, and a classification result of shale facies in the time domain is obtained; S4, the classification results of shale facies in the time domain of multiple vertical wells are obtained and analyzed, and a facies combination variation trend in the region is obtained; The specific steps of S3 comprise: S31, Lithoscanner data is obtained from the vertical well, depth correction is performed on the Lithoscanner data, the depth corrected Lithoscanner data is classified and summarized, the Lithoscanner data contains mineral substances of quartz, feldspar, calcite, dolomite and clay minerals, the facies of shale is divided according to the content of the mineral substances and the content difference of the mineral substances of different shale facies types in a three-terminal classification scheme, and a classification result of shale facies in the depth domain is obtained; S32, a ratio of depth data in the classification result of shale facies in the depth domain to the sedimentation rate at the depth data is calculated, the ratio is time data of each shale facies in the classification result of shale facies in the depth domain, and a shale facies combination in the time domain is re-composed from the time data; S33, the time domain facies is divided based on a change period of the eccentricity cycle curve in the time domain, and a secondary facies combination is obtained; S34, the secondary facies combination divided by the same eccentricity cycle curve in the time domain is composed into a primary facies combination, and the primary facies combination and the secondary facies combination jointly construct a shale facies combination in the time domain of a single vertical well.

2. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 1, characterized in that, Step S1 specifically comprises the following steps: S11, logging curve data and electrical imaging data are obtained from the vertical well, and depth correction is performed on the logging curve data and the electrical imaging data; S12, dynamic and static electrical imaging images are obtained based on the depth corrected electrical imaging logging data, and a formation layering structure feature in the electrical imaging images is picked up and recorded; S13, the formation layering structure feature, the depth corrected natural gamma data and wellbore deviation and azimuth data logging curve data are collected into a data set; S14, true vertical depth correction is performed on a natural gamma curve of a vertical well based on the data set, the influence of well deviation and formation dip on the shape of the natural gamma curve is eliminated, and a true vertical depth corrected natural gamma curve is obtained.

3. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 2, characterized in that, The method for picking up and recording the formation layer structure of the dynamic electro-imaging is that at least one point of the formation layer structure is picked up and recorded in each 1-meter depth of the dynamic electro-imaging image.

4. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 1, characterized in that, The S2 specific steps include: S21, a detrend analysis is performed on the true vertical depth corrected natural gamma ray curve to obtain a trend line of the true vertical depth corrected natural gamma ray curve; S22, a frequency spectrum analysis is performed on the trend line to obtain a frequency spectrum graph, and a high confidence frequency is obtained based on the frequency spectrum graph, and the high confidence frequency is converted to a log scale sedimentary cycle thickness to obtain the sedimentary cycle thickness corresponding to the high confidence frequency; S23, a stable eccentricity period of the formation in the area where the vertical well is located is obtained through research and analysis, and a ratio of the eccentricity period is compared with a ratio of the sedimentary cycle thickness corresponding to the power peak value of the trend line to obtain the sedimentary cycle thickness corresponding to the eccentricity period; S24, a Gaussian filter window length is calculated based on the sedimentary cycle thickness corresponding to the eccentricity period, and a Gaussian filter is performed on the true vertical depth corrected natural gamma ray curve based on the window length to obtain a depth domain eccentricity cycle curve; S25, the depth domain eccentricity cycle curve is processed according to a tuning method in cyclostratigraphy research to obtain an astronomical time scale, the depth domain eccentricity cycle curve is converted to a time domain eccentricity cycle curve based on the astronomical time scale, and a sedimentation rate value at different depths in the conversion process is calculated and recorded.

5. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 4, characterized in that, The specific steps for calculating the sedimentation rate value include: calculating a ratio of a top-bottom difference value of the depth domain eccentricity cycle curve of the vertical well to a top-bottom difference value of the time domain eccentricity cycle curve to obtain the ratio as the sedimentation rate value.

6. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 1, characterized in that, The S33 specific steps include: According to the number of periods of the change period of the time domain eccentricity cycle curve, the shale facies combination in the time domain is divided into the number of independent units, i.e., secondary facies combinations, and the secondary facies combination is named according to the facies with the largest proportion in the secondary facies combination.

7. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 1, characterized in that, The step S4 specifically includes: performing the steps S1-S3 on a plurality of vertical wells located at different positions in an area to obtain classification results of the time domain shale facies of the plurality of vertical wells, judging a facies combination change trend between the plurality of vertical wells according to the classification results of the time domain shale facies of the plurality of vertical wells, and obtaining a facies combination change trend in the area in combination with the geographical positions of the plurality of vertical wells.

8. The shale facies combination characterization and analysis method based on spectrum analysis time domain division according to claim 7, characterized in that, The specific steps for judging the facies combination change trend between the plurality of vertical wells according to the time domain shale facies combination of the plurality of vertical wells include: obtaining a proportion of a same secondary facies combination of the plurality of vertical wells according to the time domain shale facies combination of the plurality of vertical wells, judging a change condition of the proportion, and obtaining the facies combination change trend between the plurality of vertical wells.

9. The shale facies combination characterization and analysis device based on spectrum analysis time domain division, characterized in that, The method comprises the following steps: acquiring a shale rock spectrum; performing time domain division on the shale rock spectrum; and performing shale facies combination characterization and analysis based on the time domain division. The method comprises the following steps: acquiring a shale rock spectrum; performing time domain division on the shale rock spectrum; and performing shale facies combination characterization and analysis based on the time domain division.

Images