Shale reservoir plasticity method crustal stress correction method based on anisotropy

Abstract

The invention relates to the technical field of geophysical crustal stress correction, in particular to a shale reservoir plastic method crustal stress correction method based on anisotropy. The method comprises the following steps: firstly, acquiring modal components of a fast shear wave signal and a slow shear wave signal; further obtaining interference factors in each signal according to the morphological similarity between different wave crests in the frequency domain of each modal component in combination with the frequency domain difference between different modal components; further performing weighted reconstruction and band-pass filtering on the modal component based on the interference factor; further accurately extracting the arrival time of the fast and slow transverse waves by using a cross-correlation method so as to calculate anisotropic parameters; and finally, the crustal stress calculation model is corrected based on the lithology data and the anisotropy parameters, a corrected minimum horizontal principal stress value is obtained, the problem of stress calculation deviation caused by anisotropy in a complex stratum is effectively solved, and a reliable basis is provided for precise drilling trajectory design, fracturing parameter optimization and efficient reservoir utilization.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical ground stress correction, and particularly relates to a plastic method for ground stress correction of anisotropic shale reservoirs. BACKGROUND

[0002] Accurate ground stress information is an important basis for realizing low-carbon exploitation in the oil field. Due to the rapid lithology change in the vertical direction of the stratum in some areas, the clastic rock and the shale are alternately developed, the shale reservoir thin interbedded structure is developed, the clay content is high, and the anisotropy is significant; the response mechanism of the sonic logging curve in the horizontal well is different from that in the vertical well, is greatly affected by the frequency dispersion, and the correlation of the slow shear wave STC in the well section is poor, which is difficult to extract, so the anisotropy correction of the sonic logging curve is needed, so as to obtain the data reflecting the true properties of the stratum.

[0003] However, the underground situation is more complex, and different geological influencing factors are not independent individuals, which leads to the difference of the stress state of each point of the stratum, the traditional Huang pressure method in the actual processing process, due to the complex structure of the shale reservoir, the influence of the time difference of the shear wave signal, and the strong anisotropy also affects the stress calibration, which leads to the deviation of the logging curve from the true value, causes the large error of the pore pressure and the ground stress calculation, and is difficult to provide reliable basis for the accurate drilling trajectory design, the optimization of the fracturing parameters and the efficient use of the reservoir. SUMMARY

[0004] In order to solve the technical problem that the analysis of the ground stress of the existing shale reservoir is affected by the complex structure of the shale reservoir and the strong anisotropy, and the ground stress error is large, the purpose of the present application is to provide a plastic method for ground stress correction of anisotropic shale reservoirs, and the technical scheme is as follows: Obtain lithology data and logging curve data; perform signal decomposition on the array sonic data contained in the logging curve data to obtain fast shear wave signals and slow shear wave signals; Extract the modal components of each signal; in each signal, according to the morphological similarity between different wave peaks in the frequency domain of each modal component, and combining the frequency domain difference between different modal components, obtain the interference factor of each modal component; weight and reconstruct the modal components based on the interference factor and perform band-pass filtering; use the cross-correlation method to analyze the reconstructed signals to obtain the arrival time of the fast shear wave and the slow shear wave, and obtain the anisotropy parameter; Based on the lithology data and the anisotropy parameter, the ground stress calculation model is corrected to obtain the corrected minimum horizontal principal stress value.

[0005] Further, the interference factor acquisition method comprises: acquire a spectrum of each of the modal components and extract a wave peak, and construct a first feature vector with a peak value, a kurtosis and a half-peak width corresponding to each of the wave peaks; acquire an anti-interference coefficient according to a similarity of the first feature vectors between all of the wave peaks of each of the modal components; acquire an interference factor of each of the modal components according to a difference between the distribution characteristics of the peak value, the kurtosis and the half-peak width of each of the modal components and other modal components respectively in the same signal, in combination with the anti-interference coefficient.

[0006] Further, the method for acquiring the anti-interference coefficient comprises: acquire an anti-interference coefficient according to an overall feature of a cosine similarity of the first feature vectors of all of the wave peak pairs of each of the modal components.

[0007] Further, the method for acquiring the interference factor of each of the modal components comprises: take the modal components one by one as target components, take other modal components in the same signal of the target components one by one as comparison components, take the peak value, the kurtosis and the half-peak width one by one as comparison dimensions, and construct a second feature vector with a mean value and a variance of characteristic values of the comparison dimensions in each of the modal components; take an Euclidean distance of the second feature vectors of the comparison dimensions of the target component and the comparison component as a numerator, take a mapped value of a sum value of the anti-interference coefficients of the two components mapped by an exponential function with a natural constant e as a denominator, and take a ratio of the numerator and the denominator as an interference sub-factor; fuse the interference sub-factors of the target component and all of the comparison components in all of the comparison dimensions to obtain an interference factor.

[0008] Further, the method for weighting and reconstructing the modal components based on the interference factor comprises: respectively in each signal, take a result of a negative correlation mapping of the interference factor after Softmax normalization as a reconstruction weight, and weight and reconstruct the modal components based on the reconstruction weight.

[0009] Further, the method for acquiring the minimum horizontal principal stress value comprises: correct anisotropy of resistivity and acoustic travel time of the horizontal well based on the anisotropy parameter; acquire a formation pore pressure by processing acoustic travel time of overburden pressure, normal compaction pressure, static column pressure of formation water and acoustic travel time of pressure of the target layer section corrected by the anisotropy parameter through the Eaton method; Based on the viscoplastic stress relaxation constitutive relation, a mechanical relationship model between the vertical principal stress and the minimum horizontal principal stress is established; based on the Young's modulus calculated from the corrected acoustic transit time, the formation pore pressure, and the formation burial time, combined with the mechanical relationship model, the minimum horizontal principal stress value is calculated.

[0010] Furthermore, the method for correcting the anisotropy of resistivity and acoustic transit time based on the anisotropy parameters includes: The relationship between the resistivity of vertical wells and horizontal wells is established based on the anisotropy index curve, the resistivity anisotropy coefficient is obtained, and the horizontal well resistivity correction formula is obtained. The acoustic transit time of the horizontal well is corrected based on the anisotropic parameters, and the corrected acoustic transit time of the horizontal well is the slow shear wave transit time.

[0011] Furthermore, after obtaining the anisotropy parameters, the method further includes: establishing a relationship model between the anisotropy parameters and the clay content based on the lithological data.

[0012] Furthermore, the cutoff frequencies for bandpass filtering are 4kHz and 16kHz.

[0013] Furthermore, the method for obtaining the arrival time includes: Cross-correlation analysis was performed between the reconstructed fast shear wave signal and the standard fast shear wave signal in the receiver, and the time corresponding to the maximum cross-correlation coefficient was taken as the arrival time of the fast shear wave; cross-correlation analysis was also performed between the reconstructed slow shear wave signal and the standard slow shear wave signal in the receiver, and the time corresponding to the maximum cross-correlation coefficient was taken as the arrival time of the slow shear wave.

[0014] The present invention has the following beneficial effects: The present application firstly acquires fast transverse wave signals and slow transverse wave signals to provide a basis for subsequent well logging response characteristic analysis; further, in each signal, according to the morphological similarity between different wave peaks in the frequency domain of each modal component, combined with the frequency domain difference between different modal components, the interference situation of the signal component is analyzed, the interference factor is obtained, and the influence degree of each modal component in frequency dispersion and superimposed interference is represented; further, based on the interference factor, the modal components are weighted and reconstructed and band-pass filtered, and the mutual correlation method is used to accurately obtain the time when the reconstructed signals of fast transverse wave and slow transverse wave arrive at the corresponding receivers, the accuracy of anisotropy parameter calculation is improved, and the anisotropy parameter is obtained, which provides a more reliable basis for final correction; finally, based on the lithology data and the anisotropy parameter, the geostress calculation model is corrected, and the corrected minimum horizontal principal stress value is obtained, which provides technical support for large-scale exploration and development of the formation. The scheme denoises and weights the array acoustic wave signals, accurately extracts the arrival time of fast and slow transverse waves to calculate the anisotropy parameter, corrects the geostress in combination with the lithology data, realizes accurate acquisition of the minimum horizontal principal stress, and effectively solves the stress calculation deviation problem caused by anisotropy in complex formations. BRIEF DESCRIPTION OF DRAWINGS

[0015] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, and the advantages thereof, a brief introduction will be given to the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor based on these drawings.

[0016] Figure 1 A flow chart of a shale reservoir plastic method geostress correction method based on anisotropy provided by an embodiment of the present application; Figure 2 A comparison chart of acoustic travel time provided by an embodiment of the present application; Figure 3 A flow chart of an interference factor acquisition method provided by an embodiment of the present application; Figure 4 A correlation chart of clay content and transverse wave anisotropy ratio provided by an embodiment of the present application. DETAILED DESCRIPTION

[0017] In order to further illustrate the technical means and effects taken by the present application to achieve the predetermined inventive objectives, the following describes in detail the specific implementation, structure, features and effects of a shale reservoir plastic method anisotropic geostress correction method according to the present application, in combination with the accompanying drawings and preferred embodiments. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. In addition, the specific features, structures or characteristics in one or more embodiments can be combined in any suitable form.

[0018] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the present application belongs.

[0019] The following specifically describes a specific scheme of a shale reservoir plastic method anisotropic geostress correction method according to the present application in combination with the accompanying drawings.

[0020] Please refer to Figure 1 , which shows a flowchart of a shale reservoir plastic method anisotropic geostress correction method according to one embodiment of the present application, specifically including: Step S1: Obtain lithology data and logging curve data; perform signal decomposition on array acoustic wave data contained in the logging curve data to obtain fast shear wave signals and slow shear wave signals.

[0021] In one embodiment of the present application, array acoustic wave data, lithology data, geostress pore pressure experimental data, and logging curve data of horizontal wells and vertical wells are first collected. The data collection is prior art and will not be described again. The collected data provides a basis for subsequent logging response characteristic analysis.

[0022] Statistical analysis of the logging curve response characteristics of vertical wells and horizontal wells in shale target layers shows that the acoustic travel time of vertical well shale is , the acoustic travel time of horizontal well shale is , the acoustic travel time of horizontal well is affected by anisotropy and is smaller. The resistivity of vertical well shale is , the resistivity of horizontal well shale is , the resistivity of horizontal well is affected by anisotropy and is larger. The natural gamma and neutron logging curves of horizontal wells and vertical wells are consistent.

[0023] Please refer to Figure 2 , which shows a comparison chart of acoustic travel time according to one embodiment of the present application, Figure 2 which contains one vertical well acoustic travel time and three horizontal well acoustic travel times of horizontal well 1, horizontal well 2 and horizontal well 3, the vertical axis corresponds to the acoustic travel time, and Figure 2 From the comparison chart, it can be seen that the acoustic travel time of horizontal well is affected by anisotropy and is smaller.

[0024] Since the shear wave in the shale horizontal well splits into fast and slow shear waves, signal decomposition is performed on the array sonic data contained in the logging curve data to obtain fast and slow shear wave signals.

[0025] Preferably, in one embodiment of the present application, the collected array sonic data is processed by using a four-component rotation technique (prior art) to obtain corresponding fast and slow shear wave signals.

[0026] In the anisotropic shale reservoir geostress correction, the decomposed fast and slow shear wave original signals contain a mixture of propagation time, amplitude, frequency and other information, which cannot be directly used as a quantitative parameter to be substituted into a physical model for calculation, so it is necessary to extract the fast and slow shear wave interval transit times which quantitatively characterize the anisotropy of the formation. The specific calculation process is as follows: Since the shear wave in the shale horizontal well splits into fast and slow shear waves, four-component dipole sonic array data are obtained after decomposition and recombination of the original data; taking the X-direction sound source as an example, the shear wave propagating in the X-direction is emitted by the sound source, and when the shear wave propagates into the formation, due to the anisotropy of the formation, the bending wave in the orthogonal direction will split, i.e. split into fast and slow shear waves, and the corresponding polarization directions are the fast and slow principal axis directions of the formation. The shear wave splitting equation of the fast and slow shear waves generated by the shear wave splitting is: ; In the formula, represents the fast shear wave signal corresponding to the X-direction sound source propagation, represents the slow shear wave signal corresponding to the X-direction sound source propagation; represents the sound wave signal source function, represents the rotation angle.

[0027] Further, after the fast and slow shear waves are generated in the formation, they propagate at slowness s1 and s2 respectively, and after a distance (from the sound source to the receiver), they are received by the X and Y direction receivers respectively. After projecting the fast and slow waves back to the X and Y directions, the relationship between the four-component data of the cross-dipole sonic logging and the fast and slow principal bending waves can be simplified as the four-component data equation: ; In the formula, represents the time-domain signal received by the X-direction receiver after the X-direction sound source is emitted; represents the time-domain signal received by the Y-direction receiver after the X-direction sound source is emitted; represents the time-domain signal received by the Y-direction receiver after the Y-direction sound source is emitted; represents the time domain signal received by the X-direction receiver after the Y-direction sound source emits; represents the slowness corresponding to the fast transverse wave; represents the slowness corresponding to the slow transverse wave; represents the distance from the sound source to the receiver; represents the rotation angle; represents the sound wave signal source function.

[0028] wherein the X and Y directions refer to two mutually perpendicular radial directions in the instrument coordinate system.

[0029] Further, the four-component data equation is converted into a matrix form, and the matrix is diagonalized by rotating an angle θ; wherein the two diagonal elements are the fast and slow main waves. The core matrix equation converted is: ; Further, the four-component data equation is converted into a matrix form, and the matrix is diagonalized by rotating an angle θ; wherein the two diagonal elements are the fast and slow main waves. The core matrix equation converted is: ; Further, the transverse wave in the X direction is accurately decomposed into fast transverse wave signal and slow transverse wave signal, and the time when the signal wave reaches the corresponding receiver is obtained based on the decomposed fast transverse wave signal and slow transverse wave signal and , and finally the distance between the sound source and the receiver and the time when the decomposed signal wave reaches the receiver are calculated respectively to obtain the fast transverse wave time difference and the slow transverse wave time difference , and the time difference calculation formula is: ; Step S2: Extracting the modal components of each signal; in each signal, the interference factor of each modal component is obtained according to the morphological similarity between different wave peaks in the frequency domain of each modal component, combined with the frequency domain difference between different modal components; the modal components are weighted and reconstructed based on the interference factor and band-pass filtered; the arrival time of the fast transverse wave and the slow transverse wave is obtained by using the cross-correlation method to analyze the reconstructed signal, and the anisotropy parameters are obtained.

[0030] In the process of correcting the plastic method of shale reservoir based on anisotropy, the accuracy of the anisotropy parameter calculation directly affects the accuracy of the stress correction calculation, and in the actual processing process, the fast transverse wave signal and the slow transverse wave signal are disturbed by wellbore interference, instrument noise and signal dispersion, etc., which further affects the deviation of the arrival time of the fast transverse wave and the slow transverse wave at the corresponding receiver, and finally affects the accuracy of the anisotropy parameters.

[0031] Therefore, it is necessary to analyze the interference condition of the signal components, reconstruct the signal, filter and denoise, accurately obtain the arrival time of the fast and slow transverse waves at the corresponding receivers, and improve the accuracy of the anisotropy parameter calculation to provide a more reliable basis for correction.

[0032] Considering that the transverse wave will be dispersed during the propagation around the wellbore in the actual acquisition process, and the superposition of instrument noise will affect the signal change of the transverse wave at different frequency bands, jitter and peak interference will occur, the modal components of each signal are first extracted, and the modal components correspond to the signals at different frequency bands after signal decomposition; the jitter and peak interference will affect the similarity between the wave peaks in the frequency domain, and also appear in the frequency domain difference of different modal components, so in each signal, according to the shape similarity between different wave peaks in the frequency domain of each modal component, and combining the frequency domain difference between different modal components, the interference factor of each modal component is obtained, which represents the affected degree of each modal component in the dispersion and superimposed interference, and prepares for subsequent weighted reconstruction.

[0033] Preferably, in an embodiment of the present application, please refer to Figure 3 which shows a flow chart of an interference factor acquisition method provided by an embodiment of the present application, and specifically includes: Step S201, obtaining the frequency spectrum of each modal component and extracting the wave peak, and constructing the first feature vector composed of the peak value, kurtosis and half peak width corresponding to each wave peak; according to the similarity of the first feature vectors between all wave peaks of each modal component, the anti-interference coefficient is obtained.

[0034] Firstly, the frequency spectrum of each modal component is obtained and the wave peak is extracted, considering that the peak value corresponding to each wave peak reflects the energy intensity of the frequency component, the kurtosis reflects the sharpness of the wave peak shape, and the half peak width reflects the distribution range or bandwidth of the wave peak on the frequency axis, and the three reflect the morphological characteristics of the wave peak from different angles, so the first feature vector is constructed to represent the morphological characteristics of each wave peak, and to provide a basis for subsequent comparison of the morphological similarity between the wave peaks. Considering that the greater the similarity of the first feature vectors between the wave peaks in the same modal component, the smaller the jitter and peak interference, the stronger the anti-interference ability of the corresponding modal component, and the greater the anti-interference coefficient.

[0035] As an example, the morphological similarity between different wave peaks is represented by the cosine similarity between vectors, and the anti-interference coefficient is obtained according to the overall characteristics of the cosine similarity of the first feature vectors of all wave peak pairs of each modal component. Wherein, the wave peak pair is a pair of two wave peaks in the same modal component.

[0036] Specifically, the average value of all cosine similarities corresponding to each modal component is taken as the anti-interference coefficient of each modal component, and the average value represents the overall characteristics of the cosine similarity.

[0037] It should be noted that the spectrum diagram can be obtained by fast Fourier transform, and each local maximum point is taken as the peak point corresponding to each wave peak, which is a common technical means. In other embodiments of the present application, the implementer can also weight and sum the average value, median and mode of the cosine similarity to obtain the anti-interference coefficient, such as weights of 0.6, 0.2 and 0.2.

[0038] In step S202, according to the difference between the distribution characteristics of peak value, kurtosis and half-peak width of each modal component and other modal components in the same signal, and combining the anti-interference coefficient, the interference factor of each modal component is obtained.

[0039] Considering that the distribution characteristics of modal components in peak value, kurtosis and half-peak width reflect the stability and purity of signal energy distribution in the frequency band, according to the difference between the distribution characteristics of peak value, kurtosis and half-peak width of each modal component and other modal components in the same signal, the "abnormal" or "deviation" degree of the modal component relative to the overall average characteristics of the signal can be shown. The anti-interference coefficient reflects the inherent resistance of the modal component to interference, so the interference factor is obtained.

[0040] As an example, each modal component is taken as a target component, and each other modal component in the same signal of the target component is taken as a comparison component; peak value, kurtosis and half-peak width are taken as comparison dimensions; and comparison is made one by one. A modal component with small interference should have concentrated energy (high and prominent peak value), stable shape (moderate and consistent kurtosis value) and pure composition (narrow and consistent half-peak width); and a modal component with serious interference of frequency dispersion and noise superposition will have dispersed energy and distorted or widened peak shape, resulting in significant differences in these distribution characteristics (such as mean and variance) from other normal components, so the mean and variance of the characteristic values of the comparison dimensions in each modal component form a second feature vector.

[0041] Wherein, the characteristic value is the data value of the data corresponding to the comparison dimension, such as the data value of kurtosis.

[0042] The greater the Euclidean distance between the second feature vectors is, the more significant the modal difference of two modal components on the feature dimension (such as peak distribution) is, the greater the frequency domain difference is, and the greater the interference sub-factor is; the sum of the anti-interference coefficients of the two components is greater, which indicates that the frequency band where the two modal components are located is stronger in resisting interference, and the interference sub-factor is smaller; meanwhile, the anti-interference coefficient can be negative, which can affect the related logical relationship, so mapping adjustment is performed. Specifically, the Euclidean distance of the second feature vector of the target component and the contrast component in the contrast dimension is taken as the numerator, the sum of the anti-interference coefficients of the two components is mapped by an exponential function with a natural constant e as the base to obtain a mapping value as the denominator, and the ratio of the numerator and the denominator is taken as the interference sub-factor.

[0043] The exponential function is The denominator cannot be negative or zero. The exponential function is

[0044] The interference sub-factor represents the relative strength of the modal difference between the two components caused by systematic interference, so the sum of the interference sub-factors of the target component and all the contrast components in all the contrast dimensions is obtained as the interference factor.

[0045] Specifically, the sum of the interference sub-factors of the target component and other modal components in all the contrast dimensions is taken as the interference factor of the target component.

[0046] It should be noted that the analysis process of each signal is consistent, and only one example is described here, and the description is not repeated.

[0047] To effectively reduce the influence of signal quality in shale reservoirs on the arrival time of fast and slow transverse wave signals, the modal components are weighted and reconstructed based on the interference factor and band-pass filtered to provide a more reliable basis for obtaining anisotropy parameters; the cross-correlation method is used to analyze the reconstructed signals to obtain the arrival time of fast and slow transverse waves and anisotropy parameters, which provide a more reliable basis for final correction.

[0048] Preferably, in an embodiment of the present application, the greater the interference factor is, the greater the interference of the corresponding modal component to the signal is, so the weight ratio should be smaller during reconstruction, therefore, in each signal, the negative correlation of the interference factor is first mapped, and the result of the negative correlation mapping is normalized by the Softmax function, and the result is taken as the reconstruction weight, and the modal components are weighted and reconstructed based on the reconstruction weight.

[0049] As an example, the function is used for negative correlation mapping, and the weighted envelope average reconstruction method is used for weighted reconstruction.

[0050] Since the main frequency of the shear wave signal is usually 5-15 kHz (shale reservoir), the reconstructed signal is processed by using a Buterworth band-pass filtering algorithm, wherein the cut-off frequency is 4 kHz and 16 kHz, thereby further filtering out the influence of low-frequency borehole interference (<4 kHz) and high-frequency environmental noise (>16 kHz) on the signal quality.

[0051] Further, the arrival time of the fast shear wave and the slow shear wave is judged by using the cross-correlation method. Specifically, the reconstructed fast shear wave signal is cross-correlated with the standard fast shear wave signal in the receiver, and the time corresponding to the maximum cross-correlation coefficient is taken as the arrival time of the fast shear wave. The reconstructed slow shear wave signal is cross-correlated with the standard slow shear wave signal in the receiver, and the time corresponding to the maximum cross-correlation coefficient is taken as the arrival time of the slow shear wave.

[0052] The arrival time is the time taken by the sound wave to propagate from the sound source to the specific receiver (first arrival).

[0053] After analyzing the reconstructed signal to obtain the arrival time of the shear wave, the accurate fast shear wave time difference and the slow shear wave time difference are obtained by substituting the time difference calculation formula, and the ratio of the slow shear wave time difference to the fast shear wave time difference is taken as the anisotropy parameter.

[0054] After obtaining the anisotropy parameter, the relationship model between the anisotropy parameter and the clay content is established by combining the lithology data.

[0055] Specifically, the clay content experimental data is collected, and the correlation analysis between the clay content and the anisotropy parameter calculated by the fast and slow shear wave time difference is performed, and then the clay content is used to systematically represent the anisotropy parameter, and the relationship model between the anisotropy parameter and the clay content is established, and the clay content-anisotropy relationship formula calculated and analyzed is: ; In the formula, is the anisotropy parameter; a, b, and c are fitting coefficients (regression coefficients), , , ; e is a natural constant; is the clay content.

[0056] As an example, refer to Figure 4 which shows a correlation graph of the clay content and the shear wave anisotropy ratio provided by an embodiment of the present application, Figure 4 wherein the horizontal axis is the clay content, the vertical axis is the shear wave anisotropy ratio, i.e., the anisotropy parameter, the data points in the graph are experimental data points, and the curve is the fitting clay content-anisotropy relationship curve; and the corresponding relationship formula is: ; In another embodiment of the present application, normalization is performed first, and then negative correlation processing is performed. As an example, in each signal, an interference factor is taken as an input, and normalization processing is performed by using Softmax, and the processing result is , indicates the feature coefficient of the xth modal component under the influence of frequency dispersion and superimposed interference; the calculation formula of the reconstruction weight includes: ; In the formula, indicates the reconstruction weight of the xth modal component for reconstruction, indicates the feature coefficient of the xth modal component under the influence of frequency dispersion and superimposed interference; n indicates the total number of modal components.

[0057] In other embodiments of the present application, weighted superposition reconstruction can also be used for reconstruction, which is a known technology together with weighted envelope average reconstruction and Butterworth band-pass filtering algorithm; a, b and c can be: , , ; a, b and c can also be: , , , and will not be described again.

[0058] Step S3: correcting the geostress calculation model based on the lithology data and the anisotropy parameter to obtain a corrected minimum horizontal principal stress value.

[0059] After obtaining the accurate anisotropy parameter, the calculation of the geostress can be corrected, and the geostress calculation needs to be combined with the lithology data, so finally the geostress calculation model is corrected based on the lithology data and the anisotropy parameter to obtain a corrected minimum horizontal principal stress value, thereby providing a reliable basis for accurate drilling trajectory design, optimized fracturing parameters and efficient reservoir development, reducing the blindness of exploration, helping to reduce invalid drilling, avoid repeated operations, reduce construction energy consumption, and reduce the energy intensity and carbon emission intensity per unit of production from the source.

[0060] Preferably, in an embodiment of the present application, the correction process includes: First, the anisotropy of the resistivity and acoustic travel time of the horizontal well is corrected based on the anisotropy parameter to eliminate the difference in logging response caused by the development of the formation bedding. As an example, the relationship between the resistivity of the vertical well and the horizontal well is established based on the anisotropy index curve, the resistivity anisotropy coefficient is obtained, and the horizontal well resistivity correction formula is established. ; In the formula,​​ Vertical resistivity, in units of ; Horizontal resistivity, in units of ; Here, is the anisotropy parameter; d is the anisotropy coefficient, ranging from 1.5 to 3. In one embodiment of the present invention, the anisotropy is finally determined by statistical analysis of the back-calculation results of multiple wells in the target work area. .

[0061] The horizontal well resistivity correction process can effectively reflect the true response characteristics of resistivity along different directions, providing a basis for subsequent anisotropic correction of acoustic parameters.

[0062] because (Fast shear wave time difference) corresponds to the acoustic wave propagation characteristics "parallel to the bedding direction" (horizontal time difference) in horizontal well logging; while (Slow transverse wave transit time) corresponds to the acoustic wave propagation characteristic "perpendicular to the bedding direction" (vertical transit time, i.e., the true acoustic transit time of the formation). The "measured value" of acoustic transit time logging in horizontal wells (denoted as...) Essentially, it is the time difference of a fast transverse wave propagating in the horizontal direction. Because the wellbore trajectory of a horizontal well is parallel to the formation plane, the sonic logging tool mainly receives fast shear wave signals propagating in the horizontal direction; therefore, the measured time difference... = ; To obtain the slow shear wave transit time corresponding to the true formation stress state, correction based on anisotropic parameters is required; the correction formula for the sonic transit time of horizontal wells is as follows: ; Therefore, the acoustic time difference is corrected based on anisotropic parameters. During the correction process, the measured time difference... After correction This refers to the slow transverse wave time difference during vertical propagation, which reflects the actual sound wave propagation time in the strata.

[0063] Based on this, the Eaton method is used to solve for the formation pore pressure. The Eaton method comprehensively considers the sonic transit time of overburden pressure, normal compaction pressure, formation water hydrostatic pressure, and the sonic transit time of the target section pressure after anisotropic parameter correction to obtain the formation pore pressure, thus obtaining pore pressure results that are more consistent with the actual formation conditions. By introducing anisotropic parameters and corrected horizontal well sonic transit time and resistivity curves, and comprehensively considering the influence of uncompaction using the Eaton method, the formation pore pressure is calculated by combining overlying strata pressure, target section compaction pressure, normal compaction pressure, and sonic transit time parameters. The formula for calculating formation pore pressure includes: ; wherein, is the pore pressure, in MPa; is the overburden pressure, in MPa; is the formation water static column pressure, in MPa; is the normal compaction pressure acoustic travel time, i.e. the measured acoustic travel time of the acoustic wave, in μs / ft; is the target interval pressure acoustic travel time, corrected by anisotropy parameters, in μs / ft; r is the regional Eaton power function, and by default r = 5.

[0064] wherein, the target interval is the current analysis interval, the correction formula is: ; wherein, is the P-wave travel time, corresponding to the original value of ; is a correction factor function constructed according to , which is used to eliminate the propagation time error caused by anisotropy, so as to ensure that the calculation result of the pore pressure is more accurate and reliable.

[0065] wherein, the correction factor function is generally in proportion or logarithmic relationship, and is obtained by fitting actual data. Specifically, the real acoustic travel time of a straight well in a target area, the measured acoustic travel time of a horizontal well and the corresponding anisotropy parameters are collected, and the function form and the corresponding fitting constant are obtained by linear and exponential statistical regression methods respectively. After error testing and optimization (selecting the method with the smallest fitting residual), the corresponding correction factor function is obtained. The function fitting is a technology well known to those skilled in the art, and will not be described in detail.

[0066] After obtaining the formation pore pressure, a mechanical relationship model of the vertical principal stress and the minimum horizontal principal stress is established to represent the evolution law of the stress with time, according to the visco-plastic stress relaxation theory. The mechanical relationship model includes: ; wherein, is the differential stress changing with time, in MPa, reflecting the stress relaxation effect; B is a measure of the rock elastic compliance; n is the stress relaxation index, describing the deformation trend changing with time; is the total strain, is the formation burial time, in billion years.

[0067] Applying this constitutive law to the reservoir, the following can be further obtained: ; wherein, is the vertical principal stress, unit: MPa; is the minimum horizontal principal stress, unit: MPa; is the Young's modulus, which is corrected by the anisotropy parameter, unit: MPa, and reflects the elastic response characteristics of the formation.

[0068] Further, the minimum horizontal principal stress value is calculated according to the Young's modulus calculated by the corrected acoustic travel time of the horizontal well, the formation pore pressure and the formation depth time, and in combination with a mechanical relationship model.

[0069] The calculation formula of the minimum horizontal principal stress includes: Wherein, Y is a regional coefficient, which is used to quantify the influence degree of the regional geological background on the minimum horizontal principal stress; and Z is a correction coefficient, which is used to correct the systematic error in the formula derivation process.

[0070] In an embodiment of the present application, Y takes a value of 0.8-1, and Z takes a value of-1.2-0.8; as an example, Y takes 0.9, and Z takes-1; in other embodiments of the present application, the implementer can adjust by himself.

[0071] The correction process fully considers the anisotropy characteristics and stress relaxation effect of the shale reservoir, realizes the dynamic correction and accurate inversion of the horizontal well ground stress parameters, and thus significantly improves the accuracy and reliability of the ground stress prediction.

[0072] It should be noted that the lithology data includes the clay content, formation depth time and formation water static column pressure and other calculated formation characteristic data; the model and algorithm used in the correction process are known technologies, and will not be described here.

[0073] In another embodiment of the present application, after the resistivity and acoustic travel time of the horizontal well are corrected, the corrected representation effect is also verified by comparing the same formation and same lithology logging curve data of the vertical well.

[0074] After the experiment verification, it is proved that the corrected real formation logging characteristics can be obtained through the anisotropy curve, the anisotropy parameter has a high coincidence rate with the experimental data, and it is proved that the anisotropy curve has feasibility.

[0075] ​In summary, in order to solve the technical problem that the analysis of the existing shale reservoir ground stress is affected by the complex shale reservoir structure and strong anisotropy, resulting in a large error of the ground stress, the present application provides a shale reservoir plasticity method ground stress correction method based on anisotropy. The present application first acquires the fast shear wave signal and the slow shear wave signal, and extracts the modal components of each signal; further, in each signal, according to the morphological similarity between different wave peaks in the frequency domain of each modal component, combined with the frequency domain difference between different modal components, an interference factor is obtained; further, based on the interference factor, the modal components are weighted and reconstructed and band-pass filtered; further, the cross-correlation method is used to obtain the arrival time of the reconstructed fast shear wave and slow shear wave, and the anisotropy parameter is obtained; finally, based on the lithology data and the anisotropy parameter, the ground stress calculation model is corrected, and the corrected minimum horizontal principal stress value is obtained. The present scheme carries out modal decomposition and de-interference weighted reconstruction on the array acoustic wave signal, accurately extracts the arrival time of the fast and slow shear waves by using the cross-correlation method to calculate the anisotropy parameter; in combination with the lithology data, the ground stress model is corrected, the accurate minimum horizontal principal stress is obtained, and the stress calculation deviation problem caused by anisotropy in the complex formation is effectively solved, which provides a reliable basis for accurate drilling trajectory design, optimized fracturing parameters and efficient reservoir development.

[0076] It should be noted that the above-mentioned sequence of the embodiments of the present application is only for description, and does not represent the advantages and disadvantages of the embodiments. The processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multi-task processing and parallel processing are also possible or can be advantageous.

[0077] Each embodiment in the specification is described in a progressive manner, and the same and similar parts between each embodiment can be referred to each other, and each embodiment mainly describes the difference from other embodiments.

Claims

1. A method for anisotropic-based shale reservoir plastic method geostress correction, characterized in that, The method comprises: obtaining lithology data and logging curve data; performing signal decomposition on array acoustic data contained in the logging curve data to obtain fast shear wave signals and slow shear wave signals; extracting modal components of each signal; in each signal, obtaining an interference factor of each modal component according to the morphological similarity between different wave peaks in the frequency domain of each modal component, and combining the frequency domain differences between different modal components; weighting and reconstructing the modal components based on the interference factor and performing band-pass filtering; obtaining the arrival time of fast shear waves and slow shear waves by using cross-correlation method to analyze the reconstructed signals, and obtaining anisotropy parameters; correcting a ground stress calculation model based on the lithology data and the anisotropy parameters to obtain a corrected minimum horizontal principal stress value.

2. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for obtaining the interference factor comprises: obtaining a spectrum of each modal component and extracting wave peaks, and constructing a first feature vector by using a peak value, kurtosis and half-peak width corresponding to each wave peak; obtaining an anti-interference coefficient according to the similarity of the first feature vectors between all wave peaks of each modal component; in the same signal, obtaining an interference factor of each modal component according to the differences between the distribution characteristics of the peak value, the kurtosis and the half-peak width of each modal component and other modal components, and combining the anti-interference coefficient.

3. The method of claim 2, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for obtaining the anti-interference coefficient comprises: obtaining an anti-interference coefficient according to the overall characteristics of the cosine similarity of the first feature vectors of all wave peak pairs of each modal component.

4. The method of claim 2, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for obtaining the interference factor of each modal component comprises: taking each modal component as a target component one by one, taking other modal components in the same signal of the target component as comparison components one by one, taking the peak value, the kurtosis and the half-peak width as comparison dimensions one by one, and constructing a second feature vector by using the mean and variance of the characteristic values of the comparison dimensions in each modal component; taking the Euclidean distance of the second feature vectors of the target component and the comparison components in the comparison dimensions as a numerator, taking the sum of the anti-interference coefficients of the two components mapped by an exponential function with a natural constant e as a denominator, and taking the ratio of the numerator and the denominator as an interference sub-factor; fusing the interference sub-factors of the target component and all comparison components in all comparison dimensions to obtain an interference factor.

5. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for weighting and reconstructing the modal components based on the interference factor comprises: in each signal, taking the result of the negative correlation mapping of the interference factor after Softmax normalization as a reconstruction weight, and weighting and reconstructing the modal components based on the reconstruction weight.

6. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for obtaining the minimum horizontal principal stress value comprises: correcting the anisotropy of resistivity and acoustic travel time of a horizontal well based on the anisotropy parameters; obtaining the formation pore pressure by processing the acoustic travel time of the overlying rock pressure, the normal compaction pressure, the static column pressure of the formation water and the acoustic travel time of the pressure of the target layer section corrected by the anisotropy parameters through the Eaton method. Based on the visco-plastic stress relaxation constitutive relation, a mechanical relationship model of the vertical principal stress and the minimum horizontal principal stress is established; and according to the Young's modulus calculated from the corrected acoustic travel time, the formation pore pressure and the formation burial depth time, and in combination with the mechanical relationship model, the minimum horizontal principal stress value is calculated.

7. The method of claim 6, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for correcting the anisotropy of the resistivity and the acoustic travel time based on the anisotropic parameter comprises: A relationship between the vertical well and the horizontal well resistivity is established based on the anisotropic index curve, the resistivity anisotropy coefficient is obtained, and a horizontal well resistivity correction formula is obtained; The acoustic travel time of the horizontal well is corrected based on the anisotropic parameter, and the corrected acoustic travel time of the horizontal well is the slow shear wave travel time.

8. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. After the anisotropic parameter is obtained, a relationship model of the anisotropic parameter and the clay content is established in combination with the lithology data.

9. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The cut-off frequency of the band-pass filtering is 4KHz and 16kHz.

10. The method of claim 1, wherein the method is based on anisotropic shale reservoir plastic method of stress correction. The method for obtaining the arrival time comprises: The reconstructed fast shear wave signal and the standard fast shear wave signal in the receiver are cross-correlated, the time corresponding to the maximum cross-correlation coefficient is taken as the arrival time of the fast shear wave; and the reconstructed slow shear wave signal and the standard slow shear wave signal in the receiver are cross-correlated, the time corresponding to the maximum cross-correlation coefficient is taken as the arrival time of the slow shear wave.

Images