Shale thermal maturity earthquake prediction method based on v-Ro model

Abstract

The invention discloses a shale thermal maturity earthquake prediction method based on a v-Ro model, and the method specifically comprises the following steps: 1, collecting logging and earthquake velocity spectrum data, constructing a Morlet-BP neural network model, reconstructing sound wave velocity data containing lithologic information, and obtaining a low-frequency earthquake velocity through low-frequency filtering; 2, medium-high frequency seismic velocities are obtained through well-seismic combination inversion, arithmetic summation is carried out on the low-frequency seismic velocities and the high-frequency seismic velocities in the step 1, and seismic horizon velocities capable of reflecting the real situation of a shale section stratum are obtained; and step 3, converting the seismic velocity of the shale section in the step 2 into vitrinite reflectivity, and realizing plane continuous prediction of the thermal maturity of the shale. The problem that in the prior art, the plane prediction precision of the thermal maturity of the shale is low is solved.

Description

Technical Field

[0001] The invention belongs to the technical field of shale maturity prediction and relates to a shale thermal maturity seismic prediction method based on a v-Ro model. Background Art

[0002] Organic pores are the main storage space for shale gas and an important part of the natural seepage channels of shale. They contribute up to more than 50% to the total porosity of shale, and thermal maturity is an important parameter that restricts the organic pores in shale. The development types of organic pores at different evolutionary levels vary greatly, but in the late stage of overmaturity, carbonization of organic matter will destroy shale pores and increase the risk of shale gas exploration. Therefore, a comprehensive and detailed evaluation of thermal maturity is of great significance to the development of shale oil and gas.

[0003] In actual exploration, due to factors such as core mud contamination, limited sampling, and high test and analysis costs, there is less test and analysis data, resulting in insufficient understanding of the planar thermal evolution of shale, which reduces the success rate of shale oil and gas drilling. At present, there are three main methods for evaluating shale thermal maturity, namely geochemical methods, numerical simulation methods, and geophysical methods.

[0004] The more popular geochemical methods include the TTI method, the vitrinite reflectance method and the rock pyrolysis method. The advantages are that they are based on first-hand data and have high prediction accuracy, but they cannot achieve continuous evaluation of source rocks. The numerical simulation method mainly simulates the temperature and Ro parameters of the target layer by comprehensively considering the burial history, thermal history, erosion amount and other parameters of a single well. The single-point accuracy is relatively high, but it is limited to the simulation of a certain point and has no specific understanding of the thermal maturity evaluation of the entire set of shales. The geophysical method mainly extracts the attributes of the target layer to establish the correlation between the attributes and thermal maturity, and then predicts the thermal maturity. It can basically perform continuous evaluation, but the accuracy is not high. Summary of the invention

[0005] The purpose of the present invention is to provide a shale thermal maturity seismic prediction method based on the v-Ro model, which solves the problem of low accuracy of shale thermal maturity plane prediction in the prior art.

[0006] The technical solution adopted by the present invention is a shale thermal maturity seismic prediction method based on the v-Ro model, which is specifically implemented according to the following steps:

[0007] Step 1: Collect well logging and seismic velocity spectrum data, build a Morlet-BP neural network model, reconstruct the acoustic velocity data containing lithology information, and then obtain low-frequency seismic velocity through low-frequency filtering;

[0008] Step 2: obtain medium-high frequency seismic velocity by combining well-seismic inversion, perform arithmetic addition of the low-frequency seismic velocity and the high-frequency seismic velocity in step 1 to obtain the seismic layer velocity that can reflect the true situation of the shale section;

[0009] Step 3: Convert the seismic velocity of the shale section in step 2 into vitrinite reflectance to achieve a planar continuous prediction of shale thermal maturity.

[0010] The present invention is also characterized in that:

[0011] The well logging and seismic velocity spectrum data collected in step 1 are used to establish a Morlet-BP neural network model for learning and training to obtain low-frequency seismic velocity;

[0012] Step 1.1, use Morlet wavelet function to replace Sigmoid excitation function in network neurons to form a radial basis function network;

[0013] Step 1.2, select a lithology curve that is highly correlated with the acoustic wave velocity, use the Morlet-BP neural network model for learning and training, and obtain a reconstructed acoustic wave velocity curve;

[0014] Step 1.3: Using the acoustic wave velocity reconstructed in step 1.2 as output and the seismic layer velocity converted from the stacked velocity spectrum as input, a three-layer network structure is constructed for learning and training, and the seismic layer velocity obtained is subjected to low-frequency filtering with the critical frequency as the boundary to produce a low-frequency velocity model;

[0015] Step 1.3 For horizontal layered media, the stacking velocity is equal to the root mean square velocity. The continuous velocity spectrum data is converted into the root mean square velocity, which is then converted into the layer velocity after inclination correction.

[0016] When the stratum has a certain dip angle and the overburden is a uniform medium, the stacking velocity is converted to the root mean square velocity:

[0017]

[0018] The root mean square velocity is converted to seismic layer velocity using the DIX formula:

[0019]

[0020] Where V d is the superposition velocity of the reflection interface; V R is the root mean square speed (RMS); Δt 0 is the travel time difference of the same reflection interface between two adjacent velocity spectra; L is the horizontal distance between two adjacent velocity spectra; V int is the layer velocity of layer n; V r,n and V r,n-1are the root mean square velocities of the nth and n-1th layers respectively; t 0,n and t 0,n-1 The two-way reflection times of the n layer and the n-1 layer are respectively used. Taking the logging acoustic wave velocity as the target, deep learning of the seismic layer velocity is performed to obtain the low-frequency seismic velocity.

[0021] Step 2 is as follows:

[0022] Step 2.1, determine the critical frequency of the medium-high frequency velocity and the low-frequency velocity, and determine the lower limit of the medium-high frequency of the seismic data as the low cutoff frequency by comprehensively analyzing the frequency spectrum of the seismic trace near the well, that is, the upper limit of the low-frequency velocity as the high cutoff frequency;

[0023] Step 2.2: Perform arithmetic addition of the low-frequency seismic velocity obtained in step 1 and the medium-high frequency seismic velocity data to obtain high-precision seismic velocity that can reflect the actual stratigraphic conditions, and identify the shale and sandstone in the main shale section.

[0024] Step 2.1.1, select the acoustic time difference curve with the whole well section, convert it into acoustic velocity data, and fit the logging acoustic velocity spectrum curve through spectrum analysis;

[0025] Step 2.1.2, extract the frequency spectrum curve of the seismic trace near the well, perform fitting, and determine the finite bandwidth matching operator of the well data and the seismic data;

[0026] Step 2.1.3, merge the velocities according to the critical frequency, and use the operator and the seismic trace convolution operation to perform inversion to achieve the purpose of improving the low-frequency and high-frequency components; the medium-high frequency velocity obtained by the Colored inversion method is then arithmetically added with the low-frequency velocity obtained in step 1;

[0027] Step 3 is as follows:

[0028] Step 3.1, integrating the thermal evolution mechanism of organic matter and the porosity evolution mechanism, respectively obtaining the power function relationship between the porosity and thermal maturity of mud shale and between the velocity and thermal maturity of mud shale;

[0029] Step 3.2: Based on the geophysical and geochemical data of different study areas, a power function quantitative relationship between shale velocity and thermal maturity is established. Using the inverted velocity of the shale section, the vitrinite reflectance of the shale is calculated to predict the planar distribution of shale Ro.

[0030] Step 3.1 is as follows:

[0031] Step 3.1.1: Porosity gradually decreases with increasing burial depth. There is an exponential relationship between shale porosity and shale burial depth:

[0032] φ=φ o e-ch (3)

[0033] Where φ is the porosity of shale, %; φ o is the initial porosity, %; h is the burial depth of the formation, m; c is a constant;

[0034] The porosity is expressed in terms of the overlying formation pressure:

[0035] φ=f(σ) (4)

[0036] Where φ is the porosity, %; σ is the overlying formation pressure, MPa;

[0037] According to the Maxwell curve of the viscoelastic creep body, there is a linear relationship between creep velocity and time:

[0038]

[0039] where ε is the strain, dimensionless; σ o is stress, MPa; E elastic modulus, MPa; η is viscosity coefficient, kg / ms; t is time, Ma;

[0040] In the case of small deformation, the strain basically reflects the change of porosity as:

[0041] ε≈Δφ(6)

[0042] If Δφ is used to represent the change in porosity, then:

[0043]

[0044] When the compaction stage of shale is divided into a series of creep sub-stages, the porosity at a given time is expressed as:

[0045]

[0046] where φ o is the porosity at initial deposition, %; σ(t) is the overlying formation pressure varying with time, MPa;

[0047] In a continuous sedimentation basin, the logarithm of the vitrinite reflectance is linearly related to the burial depth:

[0048] LqCy o =Ah+B (9)

[0049] Among them, R o is the vitrinite reflectance, %; h is the burial depth, m; A and B are constants;

[0050] In a certain part of the basin, the geothermal gradient changes little, and the relationship between ground temperature and burial depth is expressed as a linear relationship:

[0051] T=To +Gh (10)

[0052] Where T o is the surface temperature, °C; h is the burial depth, m; G is the geothermal gradient, °C / m;

[0053] The overlying stratum pressure is proportional to the burial depth:

[0054] σ=a+bh (11)

[0055] a and b are coefficients, and σ is the overlying formation pressure, MPa.

[0056] Step 3.1.2: Porosity and TTI are basically functions of burial depth and time, but the function forms are different. Substituting the linear relationship between vitrinite reflectance and burial depth into the exponential relationship of porosity evolution, we can get:

[0057]

[0058] Since there is a linear relationship between porosity and velocity, simplifying it, we can get:

[0059] v=aR o b (13)

[0060] Where a and b are coefficients; v is the velocity of shale, m / s; R o is the vitrinite reflectance of shale, %.

[0061] The beneficial effect of the present invention is that the seismic prediction method of shale thermal maturity based on the v-Ro model considers extracting the low-frequency components in the logging acoustic wave velocity, and then supplementing the seismic medium-high frequency velocity data by means of inversion, so that the newly formed velocity body has both good lateral resolution and certain vertical recognition ability, and can distinguish mud shale from sandstone; then the power function relationship between mud shale velocity and thermal maturity is derived by integrating the organic matter thermal evolution mechanism and the porosity evolution mechanism; then the seismic velocity of the mud shale section is converted into Ro according to the quantitative formula, and the plane continuous evaluation of shale thermal maturity is realized. By predicting the distribution of shale layers and estimating thermal maturity, it is helpful to find shale "sweet spots", avoid over-maturity blocks, and reduce the risk of shale exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a flow chart of the seismic prediction method of shale thermal maturity based on the v-Ro model of the present invention;

[0063] FIG. 2( a ) is a spectrum analysis diagram of the well acoustic velocity in Example 2 of the present invention;

[0064] FIG2( b ) is a fitting curve diagram of the well acoustic wave velocity in Example 2 of the present invention;

[0065] FIG2(c) is a spectrum analysis diagram of a wellside seismic trace in Example 2 of the present invention;

[0066] FIG2(d) is a frequency spectrum fitting curve diagram of the well-side seismic trace in Example 2 of the present invention;

[0067] FIG2(e) is a graph showing the frequency spectrum matching between the well and the seismic trace in Example 2 of the present invention;

[0068] FIG2(f) is a time domain Colored inversion matching operator diagram in Example 2 of the present invention;

[0069] Figure 3 This is a comprehensive analysis diagram of the relative velocity spectrum of the well bypass in Example 3 of the present invention;

[0070] Figure 4 (a) is a superposition diagram of low frequency and relative velocity components in Example 3 of the present invention;

[0071] Figure 4 (b) is a partial frequency missing superposition diagram in Example 3 of the present invention;

[0072] Figure 4 (c) is a partial frequency repetition superposition diagram in Example 3 of the present invention;

[0073] Figure 5 is a relative velocity profile of Line 1040 in Example 4 of the present invention;

[0074] FIG6( a ) is a spectrum analysis diagram of the original seismic profile in Example 4 of the present invention;

[0075] FIG6( b ) is a spectrum analysis diagram of the Colored inversion seismic profile in Example 4 of the present invention;

[0076] FIG. 7( a ) is a graph showing the AC-GR correlation analysis in Example 4 of the present invention;

[0077] FIG7( b ) is a graph showing the AC-DEN correlation analysis in Example 4 of the present invention;

[0078] Figure 8 This is a precision analysis diagram of the acoustic wave time difference curve reconstructed by the Morlet-BP neural network in Example 4 of the present invention;

[0079] Fig. 9 is an AC-DEN correlation analysis diagram in Example 4 of the present invention;

[0080] Fig.10 is an AC-DEN correlation analysis diagram in Example 4 of the present invention;

[0081] FIG. 11( a ) is a cross-plot of shale maturity and velocity in the Upper Paleozoic Carboniferous-Permian system in the Majiatan area in Example 4 of the present invention;

[0082] FIG. 11( b ) is a cross-plot of the Cambrian-Ordovician system of the Paleozoic in the Majiatan area in Example 4 of the present invention in terms of shale maturity and velocity;

[0083] Fig.12 is a thermal maturity contour map of the Uralik Formation shale in the Majiatan area in Example 4 of the present invention;

[0084] FIG. 13( a ) is an intersection diagram of the measured and seismic predicted Ro values ​​of the Wulalike Formation shale in the Majiatan area in Example 4 of the present invention;

[0085] FIG. 13( b ) is an intersection diagram of the measured Ro of the Cambrian-Ordovician system in the Shatan 1 well in the Majiatan area and the seismic predicted Ro in Example 4 of the present invention. DETAILED DESCRIPTION

[0086] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0087] Example 1

[0088] The shale thermal maturity seismic prediction method based on the v-Ro model of the present invention is specifically implemented according to the following steps:

[0089] Step 1: Collect well logging and seismic velocity spectrum data, build a Morlet-BP neural network model, reconstruct the acoustic velocity data containing lithology information, and then obtain low-frequency seismic velocity through low-frequency filtering;

[0090] Step 2: obtain medium-high frequency seismic velocity by combining well-seismic inversion, perform arithmetic addition of the low-frequency seismic velocity and the high-frequency seismic velocity in step 1 to obtain the seismic layer velocity that can reflect the true situation of the shale section;

[0091] Step 3: Convert the seismic velocity of the shale section in step 2 into vitrinite reflectance to achieve a planar continuous prediction of shale thermal maturity.

[0092] Example 2

[0093] The present invention is based on the v-Ro model of shale thermal maturity seismic prediction method, the process is as follows Figure 1 As shown, the specific implementation steps are as follows:

[0094] Step 1: Collect well logging and seismic velocity spectrum data, build a Morlet-BP neural network model, reconstruct the acoustic velocity data containing lithology information, and then obtain low-frequency seismic velocity through low-frequency filtering;

[0095] The well logging and seismic velocity spectrum data collected in step 1 are used to establish a Morlet-BP neural network model for learning and training to obtain low-frequency seismic velocity;

[0096] Step 1.1, use Morlet wavelet function to replace Sigmoid excitation function in network neurons to form a radial basis function network;

[0097] Step 1.2, select a lithology curve that is highly correlated with the acoustic wave velocity, use the Morlet-BP neural network model for learning and training, and obtain a reconstructed acoustic wave velocity curve;

[0098] Step 1.3: Using the acoustic wave velocity reconstructed in step 1.2 as output and the seismic layer velocity converted from the stacked velocity spectrum as input, a three-layer network structure is constructed for learning and training, and the seismic layer velocity obtained is subjected to low-frequency filtering with the critical frequency as the boundary to produce a low-frequency velocity model;

[0099] Step 1.3 For horizontal layered media, the stacking velocity is equal to the root mean square velocity. The continuous velocity spectrum data is converted into the root mean square velocity, which is then converted into the layer velocity after inclination correction.

[0100] When the stratum has a certain dip angle and the overburden is a uniform medium, the stacking velocity is converted to the root mean square velocity:

[0101]

[0102] The root mean square velocity is converted to seismic layer velocity using the DIX formula:

[0103]

[0104] Where V d is the superposition velocity of the reflection interface; V R is the root mean square speed (RMS); Δt 0 is the travel time difference of the same reflection interface between two adjacent velocity spectra; L is the horizontal distance between two adjacent velocity spectra; V int is the layer velocity of layer n; V r,n and V r,n-1 are the root mean square velocities of the nth and n-1th layers respectively; t 0,n and t 0,n-1 The two-way reflection times of the n layer and the n-1 layer are respectively used. Taking the logging acoustic wave velocity as the target, deep learning of the seismic layer velocity is performed to obtain the low-frequency seismic velocity.

[0105] Step 2: obtain medium-high frequency seismic velocity by combining well-seismic inversion, perform arithmetic addition of the low-frequency seismic velocity and the medium-high frequency seismic velocity in step 1, and obtain the seismic layer velocity that can reflect the true situation of the shale section formation;

[0106] Step 2 is as follows:

[0107] Step 2.1, determine the critical frequency of the medium-high frequency velocity and the low-frequency velocity, and determine the lower limit of the medium-high frequency of the seismic data as the low cutoff frequency by comprehensively analyzing the frequency spectrum of the seismic trace near the well, that is, the upper limit of the low-frequency velocity as the high cutoff frequency;

[0108] Step 2.1.1, select the acoustic time difference curve with the whole well section, convert it into acoustic velocity data, and fit the logging acoustic velocity spectrum curve through spectrum analysis, such as Figure 2(a)-Figure 2(b) As shown;

[0109] Step 2.1.2, extract the spectrum curve of the seismic trace near the well, as shown in Figure 2(c), and perform fitting, as shown in Figure 2(d); determine the finite bandwidth matching operator of the well data and seismic data, such as Figure 2(e)-Figure 2(f) As shown;

[0110] Step 2.1.3, merge the velocities according to the critical frequency, and use the operator and the seismic trace convolution operation to perform inversion to achieve the purpose of improving the low-frequency and high-frequency components; the medium-high frequency velocity obtained by the Colored inversion method is then arithmetically added with the low-frequency velocity obtained in step 1;

[0111] Step 2.2: Perform arithmetic addition of the low-frequency seismic velocity obtained in step 1 and the medium-high frequency seismic velocity data to obtain a seismic velocity that can reflect the actual stratigraphic conditions, and identify the shale and sandstone in the main shale section.

[0112] Step 3: Convert the seismic velocity of the shale section in step 2 into vitrinite reflectance to achieve a planar continuous prediction of shale thermal maturity.

[0113] Step 3 is as follows:

[0114] Step 3.1, integrating the thermal evolution mechanism of organic matter and the porosity evolution mechanism, respectively obtaining the power function relationship between the porosity and thermal maturity of mud shale and between the velocity and thermal maturity of mud shale;

[0115] Step 3.2: Based on the geophysical and geochemical data of different study areas, a power function quantitative relationship between shale velocity and thermal maturity is established. Using the inverted velocity of the shale section, the vitrinite reflectance of the shale is calculated to predict the planar distribution of shale Ro.

[0116] Example 3

[0117] The shale thermal maturity seismic prediction method based on the v-Ro model of the present invention is specifically implemented according to the following steps:

[0118] Step 1: Collect well acoustic wave time difference data and seismic velocity spectrum data, extract the low-frequency data components in the logging acoustic wave velocity, and establish a Morlet-BP neural network model for learning and training in order to construct a better trend background velocity and obtain high-precision low-frequency seismic velocity.

[0119] This model is based on the basic principles and applications of the traditional BP algorithm, using Morlet wavelet function to replace the Sigmoid excitation function in the network neurons to form a radial basis function network.

[0120] A lithology curve with a high correlation with the acoustic wave velocity is selected, and a neural network model is used for learning and training to obtain a reconstructed acoustic wave velocity curve, which can not only reflect the law of seismic layer velocity changes, but also has a certain lithology identification ability.

[0121] With the reconstructed acoustic wave velocity as output and the seismic layer velocity converted from the stacked velocity spectrum as input, a three-layer network structure is constructed for learning and training. The seismic layer velocity obtained is subjected to low-frequency filtering with the critical frequency as the boundary to produce a low-frequency velocity model.

[0122] For horizontal layered media, stacking velocity is equal to root mean square velocity; when converting continuous velocity spectrum data into root mean square velocity, the conversion of layer velocity must first consider the formation dip. When the formation has a certain dip and the overburden is a uniform medium, the stacking velocity is converted to root mean square velocity according to the following formula:

[0123]

[0124] The root mean square velocity is converted into seismic layer velocity by the DIX formula:

[0125]

[0126] Where V d is the superposition velocity of the reflection interface; V R is the root mean square speed (RMS); Δt 0 is the travel time difference of the same reflection interface between two adjacent velocity spectra; L is the horizontal distance between two adjacent velocity spectra; V int is the layer velocity of layer n; V r,n and V r,n-1 are the root mean square velocities of the nth and n-1th layers respectively; t 0,n and t 0,n-1 are the round-trip reflection times of layer n and layer n-1 respectively.

[0127] The accuracy of seismic layer velocity converted from stacked velocity spectrum is very low. Taking logging acoustic wave velocity as the target, deep learning of layer velocity is performed to obtain velocity with higher accuracy.

[0128] Step 2: Use inversion to add the low-frequency data in step 1 to the medium-high frequency velocity data of the seismic data, distinguish the shale and sandstone, and compensate the low-frequency velocity information into the medium-high frequency velocity. In order to ensure the effectiveness and integrity of the velocity compensation, the newly formed velocity body has both good lateral resolution and certain vertical recognition ability.

[0129] Velocity filtering is performed at the critical frequencies of low-frequency velocity and medium-high frequency velocity respectively. By comprehensively analyzing the spectrum of the seismic trace near the well, the lower limit of the medium-high frequency of the seismic data is determined as the low cutoff frequency, which is also the upper limit of the low-frequency velocity as the high cutoff frequency, such as Figure 3 As shown in the figure, the acoustic time difference curves with the whole well section are screened and converted into acoustic velocity data. Through spectrum analysis, the logging acoustic velocity spectrum curves are fitted and merged into high-precision seismic layer velocity.

[0130] Extract the frequency spectrum curve of the seismic trace near the well, perform fitting, and determine the limited bandwidth matching operator of the well data and seismic data; when merging the velocity, it is necessary to pay attention to strictly merge according to the critical frequency; Figure 4 As shown in (a), the information merging method in this embodiment is as follows: Figure 4 (b) and Figure 4 The merging method shown in (c) will result in missing information and repeated superposition, which will inevitably affect the final result.

[0131] The operator and seismic trace convolution operation are used for inversion to improve the low-frequency and high-frequency components. The medium-high frequency velocity and low-frequency velocity obtained by Colored inversion are effectively combined to obtain the seismic velocity that can reflect the real formation situation. The seismic true velocity is used for lithology identification to identify the main shale section.

[0132] The Colored inversion method is used to obtain medium-high frequency seismic velocities without wavelet extraction, so it is not affected by wavelet extraction; the calculation is simple and fast, requiring only one convolution operator; it is faithful to earthquakes, and performs band-limited inversion within the effective bandwidth of seismic data, without excessive human interference, and has high stability; it is globally optimized and highly reliable. Therefore, Colored inversion can be applied to more complex areas, such as special work areas such as complex lithology areas, complex structural areas, and low exploration areas.

[0133] Step 3: Integrate the thermal evolution mechanism of organic matter and the porosity evolution mechanism to derive the power function relationship between shale velocity and thermal maturity; convert the seismic velocity of the shale section into Ro to achieve planar continuous evaluation of shale thermal maturity, and provide a method flow for shale layer distribution prediction and thermal maturity estimation.

[0134] Combining the porosity evolution mechanism and the organic matter thermal evolution mechanism, the power function relationship between the shale porosity and thermal maturity is derived. Athy (1930) first pointed out that there is an exponential relationship between porosity and shale burial depth based on the fact that porosity gradually decreases with increasing burial depth:

[0135] φ=φ o e -ch (3)

[0136] Where φ is the porosity of shale, %; φ o is the initial porosity, %; h is the burial depth of the formation, m; c is a constant.

[0137] The evolution of porosity is only related to the pressure of the overlying formation under normal compaction conditions. When the burial depth of shale is constant, the matrix and pore fluid jointly bear the overlying formation pressure, and the abnormal pressure generated by pore fluid drainage is suppressed. The change of shale porosity is affected by the above two factors. The higher the pressure, the higher the porosity. Since the fundamental cause of particle strain and fluid pressure is the overlying formation pressure, the porosity can be expressed by the overlying formation pressure:

[0138] φ=f(σ) (4)

[0139] Where φ is the porosity, %; σ is the overlying formation pressure, MPa.

[0140] When considering the time factor, the degree of rock compaction under the same pressure varies with the duration of overloading. The mechanical properties of the formation are similar to those of rheological bodies, showing creep characteristics. Under certain conditions of the overburden, the medium is gradually compacted and the rock porosity gradually decreases over time;

[0141] According to the Maxwell curve of viscoelastic creep body, there is a linear relationship between creep velocity and time:

[0142]

[0143] where ε is the strain, dimensionless; σ o is stress, MPa; E is elastic modulus, MPa; η is viscosity coefficient, kg / ms; t is time, Ma.

[0144] In the case of small deformation, the strain basically reflects the change in porosity:

[0145] ε≈Δφ (6)

[0146] If Δφ is used to represent the change in porosity, then:

[0147]

[0148] When the compaction stage of shale is divided into a series of creep sub-stages, the porosity at a given time can be expressed as:

[0149]

[0150] where φ o is the porosity at the time of initial deposition, %; σ(t) is the overlying formation pressure that changes with time, MPa; Obviously, when the shale compaction process is observed as the sum of a series of creep processes during the increase of burial depth, the porosity is affected by the overlying formation pressure and time.

[0151] Many scholars have conducted a comprehensive study on the relationship between the reflectivity of the vitrinite and the burial depth (Dow, 1977; Tissot et al., 1987). The reflectivity of the vitrinite increases with the increase of burial depth (temperature), and is linear in the semi-logarithmic coordinate, which has been proved by experiments and practice. In a continuous sedimentation basin, the logarithm of the reflectivity of the vitrinite is linearly related to the burial depth:

[0152] LqCy o =Ah+B (9)

[0153] Where R o is the vitrinite reflectance, %; h is the burial depth, m; A and B are constants.

[0154] From the above derivation, it can be concluded that there is a linear relationship between porosity and reaction time.

[0155] Porosity is affected by overburden pressure, and TTI is controlled by geotemperature.

[0156] In a certain part of the basin, the geothermal gradient changes little, and the relationship between ground temperature and burial depth is expressed as a linear relationship:

[0157] T=T o +Gh (10)

[0158] Where T o is the surface temperature, ℃; h is the burial depth, m; G is the geothermal gradient, ℃ / m.

[0159] The overlying stratum pressure is proportional to the burial depth:

[0160] σ=a+bh (11)

[0161] a and b are coefficients, and σ is the overlying formation pressure, MPa.

[0162] Porosity and TTI are basically functions of burial depth and time, but the function forms are different. Substituting the linear relationship between vitrinite reflectance and burial depth into the exponential relationship of porosity evolution, we can get:

[0163]

[0164] Since there is a linear relationship between porosity and velocity, simplifying it, we can get:

[0165] v=aR o b (13)

[0166] Where a and b are coefficients; v is the velocity of shale, m / s; R o is the vitrinite reflectance of shale, %.

[0167] According to the geophysical and geochemical data of different study areas, a quantitative relationship between shale velocity and thermal maturity is established, and the vitrinite reflectance of shale is calculated using the inverted velocity of the shale section.

[0168] The present invention analyzes the spectrum of seismic data and finds that seismic data usually contains medium-high frequency information but lacks low-frequency information. The spectrum of well logging data contains both low-frequency information and medium-high frequency information, that is, well logging data has a relatively full-segment spectrum information, so its vertical resolution is relatively high. If the low-frequency part of the well logging data is compensated into the seismic data during seismic inversion, the resolution of the seismic data can be improved, and the reconstructed seismic data can be used in actual work, which may produce better results.

[0169] Example 4

[0170] The western margin of the Ordos Basin refers to a large fault-fold belt that stretches from the Inner Mongolia Autonomous Region in the north to Longxian County in Shaanxi Province in the south. It stretches more than 600 km from north to south and is 50 to 100 km wide from east to west. For a long time, the western margin of the basin has been the focus of oilfield exploration. It is the strategic successor area for the next step of oil and gas exploration and the main block for increasing reserves and production. Many exploration wells have obtained low-yield gas flow in the Uralik Formation on the western margin. Zhongping 1 well and Zhongping 4 well have obtained 6.42×10 4 m 3 / d and 4×10 4 m 3 / d of industrial gas flow, and breakthrough progress has been made in shale gas exploration on the western edge, indicating that although the shale in this study area has a low organic matter abundance content, it has a strong hydrocarbon generation potential and is an effective marine shale.

[0171] The Majiatan area is a transition zone between the southern and northern parts of the western margin. Under the joint influence and superposition of the Liupan thrust tectonic system in the south and the Hetao-Yinchuan extensional tectonic system in the north, it has developed structural styles such as fault propagation folds, double structures, imbricate structures, structural wedges and back-thrust fault combinations (thrust structures), forming oil and gas trap types such as fault anticlines, fault noses and fault blocks. Therefore, shale gas development has great technical difficulties. Accurate prediction of the distribution, development, and evolution of shale is helpful to find "sweet spots" and reduce drilling risks. Among them, shale thermal maturity is a very important parameter. The pore types and development of shales with different evolution degrees are quite different. If the evolution degree is too low, it is not conducive to the formation of organic pores, and the physical properties are poor; if the evolution degree is too high, the carbonization of organic matter will destroy the shale pores and deteriorate the physical properties. Therefore, a comprehensive and detailed evaluation of shale thermal maturity is helpful to the study of the occurrence mechanism of shale oil and gas, and can also provide a geological basis for the development of shale oil and gas and reduce drilling risks.

[0172] By converting the well logging acoustic time difference data in the study area into acoustic wave velocity data, matching it with the seismic data, forming a matching operator, and performing Colored inversion, the relative velocity, that is, the medium-high frequency velocity, is obtained, such as Figure 5 As shown in Figure 6, negative values ​​indicate shale, positive values ​​indicate sandstone, and the relative velocity is basically distributed around -800 to 1400 m / s, representing the relative size of the seismic layer velocity. Comparing the relative velocity profile and the original seismic profile, the main frequency of the original seismic profile is 20 Hz, while the main frequency of the information obtained by Colored inversion is 9 Hz, as shown in Figure 6 (a) and Figure 6 (b), indicating that Colored inversion has significantly lowered the main frequency of the seismic data. The information obtained by Colored inversion has increased the components of 0 to 5 Hz and 60 to 80 Hz, indicating that Colored inversion has significantly supplemented the low-frequency and high-frequency components in the seismic data, which is an effective improvement of the seismic data.

[0173] Comprehensive analysis shows that, as shown in Figure 7(a) and Figure 7(b), the acoustic time difference data and the natural gamma and density data have a high correlation. The AC-GR correlation analysis and AC-DEN correlation analysis graphs use these two curves to reconstruct the acoustic time difference curve, and the correlation coefficient can reach above 0.75.

[0174] According to the constructed Morlet-BP neural network model, the natural gamma and density curves are used as input, and the acoustic time difference curve is used as output to reconstruct the acoustic time difference data in order to enhance the ability of the acoustic time difference curve to identify lithology. The calculation results show that if Figure 8As shown in the figure, the Morlet-BP neural network model has a strong prediction ability with a compliance rate of more than 80%. The reconstructed acoustic wave velocity data is used to establish a low-frequency velocity model. The seismic layer velocity is used as the input of the network again, and the reconstructed acoustic wave time difference data is used as the output. The seismic layer velocity converted using the DIX formula is reconstructed, and then under the constraint of the isochronous stratigraphic framework, the critical frequency is used as the boundary, as shown in the figure. Fig. 9 As shown in the figure, the reconstructed seismic layer velocity is low-frequency filtered to obtain the low-frequency velocity. The low-frequency velocity is also called the background trend velocity, which can reflect the velocity change trend of the stratum. From the low-frequency profile, the background velocity of the study area is basically distributed between 1000 and 7000 m / s. The deeper the seismic layer, the greater the velocity, and the lateral direction basically conforms to the change trend of the structure.

[0175] The relative velocity and low-frequency velocity models are combined to obtain the absolute seismic layer velocity that can reflect the true situation of the stratum. The absolute seismic layer velocity reflects the true change of the seismic layer velocity. The seismic layer velocity in the study area is mainly distributed between 1000 and 7000 m / s. On the whole, the greater the depth, the greater the velocity. The lateral seismic layer velocity changes with the changes in structural undulations and sedimentary facies, such as Fig.10 As shown in Figure 11(a), it is relatively easy to identify the distribution of favorable sand bodies. Based on the measured Ro data and shale acoustic velocity data in the study area, the power function relationship between the two is fitted. Since the compaction laws of the Upper Paleozoic and Lower Paleozoic in the study area are quite different, it is necessary to fit the models separately. As shown in Figure 11(a), the fitting model of the Carboniferous-Permian system in the Upper Paleozoic is φ=13.403x -1.164 , R 2 =0.7055; as shown in Figure 11(b), the fitting model of the Cambrian-Ordovician system of the Lower Paleozoic is φ=7.418x -4.454 , R 2 =0.7215. When making predictions, it is necessary to stratify the systems and use different models for prediction.

[0176] Based on the quantitative relationship between the velocity and thermal maturity of the Uralik Formation shale, the Ro plane distribution of the shale is predicted, such as Fig.12 As shown. According to the prediction results, the overall shale maturity in the Majiatan area is relatively high, generally reaching between 1.5% and 1.9%, which is a high-maturity shale. From the perspective of change trend, there are two centers with higher maturity in the center of the study area, which gradually decreases to the surrounding areas, generally reaching 1.7% to 2.0%, which is a high-maturity shale, and locally reaching the over-maturity stage. From the perspective of mining, the overall maturity is relatively high, and the organic pores are in the growth stage. The local area in the center is in the over-maturity stage, and there may be a situation where the organic carbonization loses pores.

[0177] In order to verify the accuracy of the prediction, the seismic predicted Ro value and the measured Ro value in the exploration well were correlated. The seismic prediction results and error analysis of the thermal maturity of the Uralik Formation shale in the Majiatan area are shown in Figure 13(a). It is found that the prediction coincidence rate of the Uralik Formation shale Ro in the Zhong 2, Zhong 5, Kushen 1 and Maji 1 wells is 90%, as shown in Figure 13(b). The prediction coincidence rate of the Shatan 1 well is 89.6%, both of which can achieve high accuracy.

[0178] The present invention is based on the v-Ro model of the seismic prediction method of shale thermal maturity, which supplements the low-frequency information in the well logging and velocity spectrum into the seismic data to complete the spectrum of the seismic data, making it closer to the spectrum of the well logging data, thereby ensuring the lateral resolution of the seismic data while enhancing the vertical resolution of the seismic data; comprehensively derives the power function relationship between v-Ro based on the organic matter thermal evolution mechanism and the porosity evolution mechanism theory; combines geophysical inversion and geological theory derivation to preliminarily predict the thermal maturity of shale in complex tectonic areas, and achieves good results.

Claims

1. Seismic prediction method of shale thermal maturity based on v-Ro model, It is characterized in that Follow the steps below to implement it: Step 1: Collect well logging and seismic velocity spectrum data, build a Morlet-BP neural network model, reconstruct the acoustic velocity data containing lithology information, and then obtain low-frequency seismic velocity through low-frequency filtering; Step 2: obtain medium-high frequency seismic velocity by combining well-seismic inversion, perform arithmetic addition of the low-frequency seismic velocity and the high-frequency seismic velocity in step 1 to obtain the seismic layer velocity that can reflect the true situation of the shale section; Step 3: Convert the seismic velocity of the shale section in step 2 into vitrinite reflectance to achieve a planar continuous prediction of shale thermal maturity.

2. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 1, It is characterized in that The logging and seismic velocity spectrum data collected in step 1 are used to establish a Morlet-BP neural network model for learning and training to obtain low-frequency seismic velocity; Step 1.1, use Morlet wavelet function to replace Sigmoid excitation function in network neurons to form a radial basis function network; Step 1.2, select a lithology curve that is highly correlated with the acoustic wave velocity, use the Morlet-BP neural network model for learning and training, and obtain a reconstructed acoustic wave velocity curve; Step 1.3: Take the acoustic wave velocity reconstructed in step 1.2 as output and the seismic layer velocity converted from the stacked velocity spectrum as input, build a three-layer network structure for learning and training, and perform low-frequency filtering on the obtained seismic layer velocity with the critical frequency as the boundary to produce a low-frequency velocity model.

3. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 2, It is characterized in that In step 1.3, for horizontal layered media, the stacking velocity is equal to the root mean square velocity, and the continuous velocity spectrum data is converted into the root mean square velocity, which is then converted into the layer velocity after the inclination correction; When the stratum has a certain dip angle and the overburden is a uniform medium, the stacking velocity is converted to the root mean square velocity: The root mean square velocity is converted to seismic layer velocity using the DIX formula: Where V d is the superposition velocity of the reflection interface; V R is the root mean square speed (RMS); Δt 0 is the travel time difference of the same reflection interface between two adjacent velocity spectra; L is the horizontal distance between two adjacent velocity spectra; V int is the layer velocity of layer n; V r,n and V r,n-1 are the root mean square velocities of the nth and n-1th layers respectively; t 0,n and t 0,n-1 The two-way reflection times of the n layer and the n-1 layer are respectively used. Taking the logging acoustic wave velocity as the target, deep learning of the seismic layer velocity is performed to obtain a velocity with higher accuracy.

4. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 2, It is characterized in that The step 2 is specifically as follows: Step 2.1, determine the critical frequency of the medium-high frequency velocity and the low-frequency velocity, and determine the lower limit of the medium-high frequency of the seismic data as the low cutoff frequency by comprehensively analyzing the frequency spectrum of the seismic trace near the well, that is, the upper limit of the low-frequency velocity as the high cutoff frequency; Step 2.2: Perform arithmetic addition of the low-frequency seismic velocity obtained in step 1 and the medium-high frequency seismic velocity data to obtain a seismic velocity that can reflect the actual stratigraphic conditions, and identify the shale and sandstone in the main shale section.

5. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 4, It is characterized in that The step 2.1 is specifically as follows: Step 2.1.1, select the acoustic time difference curve with the whole well section, convert it into acoustic velocity data, and fit the logging acoustic velocity spectrum curve through spectrum analysis; Step 2.1.2, extract the frequency spectrum curve of the seismic trace near the well, perform fitting, and determine the finite bandwidth matching operator of the well data and the seismic data; Step 2.1.3: Merge the velocities according to the critical frequency, and use the operator and the seismic trace convolution operation to perform inversion to achieve the purpose of improving the low-frequency and high-frequency components; the medium-high frequency velocities obtained by the Colored inversion method are then arithmetically added with the low-frequency velocities obtained in step 1.

6. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 5, It is characterized in that The step 3 is specifically as follows: Step 3.1, integrating the thermal evolution mechanism of organic matter and the porosity evolution mechanism, respectively obtaining the power function relationship between the porosity and thermal maturity of shale and between the velocity and thermal maturity of shale; Step 3.2: Based on the geophysical and geochemical data of different study areas, a power function quantitative relationship between shale velocity and thermal maturity is established. Using the inverted velocity of the shale section, the vitrinite reflectance of the shale is calculated to predict the planar distribution of shale Ro.

7. The shale thermal maturity seismic prediction method based on the v-Ro model according to claim 6, It is characterized in that The step 3.1 is specifically as follows: Step 3.1.1: Porosity decreases gradually with increasing burial depth. There is an exponential relationship between shale porosity and shale burial depth: f=f o e -ch (3) Where φ is the porosity of shale, %; φ o is the initial porosity, %; h is the burial depth of the formation, m; c is a constant; The porosity is expressed in terms of the overlying formation pressure: φ=f(σ) (4) Where φ is the porosity, %; σ is the overlying formation pressure, MPa; According to the Maxwell curve of the viscoelastic creep body, there is a linear relationship between creep velocity and time: where ε is the strain, dimensionless; σ o is stress, MPa; E is elastic modulus, MPa; η is the viscosity coefficient, kg / ms; t is the time, Ma; In the case of small deformation, the strain basically reflects the change of porosity as: ε≈Δφ(6) If Δφ is used to represent the change in porosity, then: When the compaction stage of shale is divided into a series of creep sub-stages, the porosity at a given time is expressed as: where φ o is the porosity at initial deposition, %; σ(t) is the overlying formation pressure varying with time, MPa; In a continuous sedimentation basin, the logarithm of the vitrinite reflectance is linearly related to the burial depth: lnR o =Ah+B (9) Among them, R o is the vitrinite reflectance, %; h is the burial depth, m; A and B are constants; In a certain part of the basin, the geothermal gradient changes little, and the relationship between ground temperature and burial depth is expressed as a linear relationship: T=T o +Gh (10) Where T o is the surface temperature, °C; h is the burial depth, m; G is the geothermal gradient, °C / m; The overlying stratum pressure is proportional to the burial depth: σ=a+bh (11) a and b are coefficients, σ ​​is the overlying formation pressure, MPa; Step 3.1.2: Porosity and TTI are basically functions of burial depth and time, but the function forms are different; substituting the linear relationship between vitrinite reflectance and burial depth into the exponential relationship of porosity evolution, we can get: Since there is a linear relationship between porosity and velocity, simplifying it, we can get: v=aR o b (13) Where a and b are coefficients; v is the velocity of shale, m / s; R o is the vitrinite reflectance of shale, %.

Images