Anisotropy-based shale reservoir plastic method geostress correction method

Abstract

The present application relates to the technical field of geophysical ground stress correction, in particular to a shale reservoir plastic method ground stress correction method based on anisotropy. The present application firstly acquires modal components of fast shear wave signals and slow shear wave signals; further, in each signal, according to the shape similarity between different wave crests in the frequency domain of each modal component, and combining the frequency domain difference between different modal components, an interference factor is acquired; further, the modal components are weighted and reconstructed based on the interference factor and band-pass filtered; further, the cross-correlation method is used to accurately extract the arrival time of fast and slow shear waves to calculate anisotropy parameters; finally, based on lithology data and anisotropy parameters, a ground stress calculation model is corrected, and the corrected minimum horizontal principal stress value is obtained, effectively solving the stress calculation deviation problem caused by anisotropy in complex strata, and providing a reliable basis for accurate drilling trajectory design, optimized fracturing parameters and efficient reservoir development.

Description

Technical Field

[0001] This invention relates to the field of geophysical stress correction technology, specifically to a stress correction method for anisotropic shale reservoirs based on the plasticity method. Background Technology

[0002] Accurate geostress information is a crucial foundation for achieving low-carbon extraction in the oil industry. In some areas, the lithology changes rapidly vertically, with alternating clastic rocks and shale, and shale reservoirs exhibiting thin interbedded structures, high clay content, and significant anisotropy. This results in a different acoustic logging response mechanism in horizontal wells compared to vertical wells, making them more susceptible to dispersion. Furthermore, the correlation of slow shear wave (STC) curves within the well section is poor, making extraction difficult. Therefore, anisotropy correction of acoustic logging curves is necessary to obtain data that truly reflects the formation properties.

[0003] However, the underground conditions are quite complex, and different geological influencing factors are not independent entities. This leads to differences in the stress state at different points in the formation. In actual processing, the traditional Huang pressure method is affected by the complex structure of shale reservoirs, which affects the judgment of shear wave signal time difference. At the same time, the strong anisotropy also affects stress calibration, causing the logging curve to deviate from the true value. This results in large errors in the calculation of pore pressure and geostress, making it difficult to provide a reliable basis for accurate drilling trajectory design, optimization of fracturing parameters, and efficient reservoir utilization. Summary of the Invention

[0004] To address the technical problem that existing shale reservoir in-situ stress analysis is affected by complex shale reservoir structures and strong anisotropy, leading to large in-situ stress errors, the present invention aims to provide an anisotropic shale reservoir plasticity-based in-situ stress correction method. The specific technical solution adopted is as follows: Acquire lithological data and well logging data; perform signal decomposition on the array acoustic data contained in the well logging data to obtain fast shear wave signals and slow shear wave signals; Extract the modal components of each signal; within each signal, based on the morphological similarity between different peaks in the frequency domain of each modal component and the frequency domain differences between different modal components, obtain the interference factor of each modal component; based on the interference factor, reconstruct the modal components by weighting and bandpass filtering; use the cross-correlation method to analyze the reconstructed signal to obtain the arrival times of the fast and slow shear waves, and obtain the anisotropy parameters; Based on the lithological data and the anisotropic parameters, the geostress calculation model is modified to obtain the corrected minimum horizontal principal stress value.

[0005] Furthermore, the method for obtaining the interference factor includes: The spectrum of each modal component is obtained and the peaks are extracted. The peak value, kurtosis and half-peak width of each peak are used to form a first feature vector. The anti-interference coefficient is obtained based on the similarity of the first feature vectors between all the peaks of each modal component. Within the same signal, based on 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 respectively, and in conjunction with the anti-interference coefficient, the interference factor of each modal component is obtained.

[0006] Furthermore, the method for obtaining the anti-interference coefficient includes: The anti-interference coefficient is obtained based on the overall characteristics of the cosine similarity of the first feature vectors of all the peak pairs of each modal component.

[0007] Furthermore, the method for obtaining the interference factor of each of the modal components includes: Each modal component is taken as a target component, and each other modal components within the same signal of the target component are taken as comparison components; the peak value, the kurtosis, and the half-peak width are taken as comparison dimensions; the mean and variance of the eigenvalues ​​of the comparison dimensions within each modal component are used to form a second feature vector; The Euclidean distance between the target component and the comparison component in the second feature vector of the comparison dimension is used as the numerator, the sum of the anti-interference coefficients of the two components is mapped by an exponential function with the natural constant e as the base as the denominator, and the ratio of the fractions is used as the interference factor. The interference factor is obtained by fusing the target component with the interference sub-factors of all the contrast components across all the contrast dimensions.

[0008] Furthermore, the method for weighted reconstruction of the modal components based on the interference factor includes: Within each signal, the negative correlation mapping result of the interference factor is normalized by Softmax and used as the reconstruction weight. The modal components are then reconstructed based on the reconstruction weight.

[0009] Furthermore, the method for obtaining the minimum horizontal principal stress value includes: The anisotropy of resistivity and sonic transit time in horizontal wells is corrected based on the aforementioned anisotropy parameters. The formation pore pressure is obtained by processing the acoustic transit time of overlying strata pressure, normal compaction pressure, formation water hydrostatic column pressure, and target stratum pressure after anisotropic parameter correction using 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: This invention first acquires fast and slow shear wave signals to provide a foundation for subsequent well logging response characteristic analysis. Further, within each signal, based on the morphological similarity between different peaks in the frequency domain of each modal component and the frequency domain differences between different modal components, the interference status of the signal components is analyzed, and an interference factor is obtained to characterize the degree of influence of dispersion and superposition interference on each modal component. Further, based on the interference factor, the modal components are weighted and reconstructed, and bandpass filtered. The cross-correlation method is then used to accurately obtain the arrival times of the fast and slow shear waves of the reconstructed signal at the corresponding receivers, improving the accuracy of anisotropy parameter calculation and obtaining the anisotropy parameters to provide a more reliable basis for final correction. Finally, based on lithological data and anisotropy parameters, the geostress calculation model is corrected to obtain the corrected minimum horizontal principal stress value, providing technical support for large-scale exploration and development of strata. This scheme performs denoising and weighted reconstruction on the array acoustic signal, accurately extracts the arrival times of fast and slow shear waves to calculate anisotropic parameters, and corrects the in-situ stress by combining lithological data, thereby achieving accurate acquisition of the minimum horizontal principal stress and effectively solving the problem of stress calculation deviation caused by anisotropy in complex strata. Attached Figure Description

[0015] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 A flowchart illustrating a method for in-situ stress correction in anisotropic shale reservoirs using the plasticity method, as provided in an embodiment of the present invention; Figure 2 A comparison diagram of acoustic time difference provided in one embodiment of the present invention; Figure 3 A flowchart illustrating a method for obtaining an interference factor according to an embodiment of the present invention; Figure 4 An embodiment of the present invention provides a correlation diagram between clay content and transverse wave anisotropy ratio. Detailed Implementation

[0017] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of an anisotropic shale reservoir plasticity-based in-situ stress correction method proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, 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 this invention pertains.

[0019] The following description, in conjunction with the accompanying drawings, details a specific scheme for an anisotropic shale reservoir plasticity method for in-situ stress correction provided by the present invention.

[0020] Please see Figure 1 The diagram illustrates a flowchart of an anisotropic shale reservoir plasticity method for in-situ stress correction, provided by an embodiment of the present invention, specifically including: Step S1: Acquire lithological data and well logging data; perform signal decomposition on the array acoustic data contained in the well logging data to obtain fast shear wave signals and slow shear wave signals.

[0021] In one embodiment of the present invention, array acoustic wave data, lithological data, geostress pore pressure experimental data, and logging curve data of horizontal and vertical wells are first collected. Data acquisition is an existing technology and will not be described in detail here. The collected data provides a basis for subsequent logging response characteristic analysis.

[0022] Statistical analysis of the logging response characteristics of vertical wells and horizontal wells with shale target layers leads to the conclusion that the sonic transit time of vertical wells in shale formations is: The sonic transit time of horizontal well shale is The sonic transit time in horizontal wells is less affected by anisotropy. The resistivity of shale in vertical wells is... The resistivity of shale in horizontal wells is The resistivity of horizontal wells is more affected by anisotropy, and the logging curves of natural gamma and neutron are consistent between horizontal and vertical wells.

[0023] Please see Figure 2 It shows a comparison diagram of acoustic time difference provided by an embodiment of the present invention. Figure 2 It includes the acoustic transit time of one vertical well, and the acoustic transit time of three horizontal wells: horizontal well 1, horizontal well 2, and horizontal well 3. The vertical axis corresponds to the acoustic transit time. Figure 2 It can be seen that the acoustic time difference of horizontal wells is less affected by anisotropy.

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

[0025] Preferably, in one embodiment of the present invention, a four-component rotation technique (existing technology) is used to process the collected array acoustic wave data to obtain the corresponding fast shear wave signal and slow shear wave signal.

[0026] In anisotropic shale reservoir in-situ stress correction, the original fast and slow shear wave signals obtained by decomposition contain a mixture of information such as propagation time, amplitude, and frequency, which cannot be directly substituted into the physical model for calculation as quantitative parameters. Therefore, it is necessary to extract the fast shear wave transit time, which quantitatively characterizes the anisotropy of the formation. and slow transverse wave time difference The specific calculation process is as follows: Since shear waves in shale horizontal wells split into fast and slow shear waves, decomposing and recombining the original data yields dipole acoustic array data with four components. Taking a sound source in the X direction as an example, a shear wave propagating in the X direction is emitted by the sound source. When the shear wave propagates into the formation, due to the anisotropy of the formation, the bending wave in its orthogonal direction will split, i.e., split into fast and slow shear waves. The corresponding polarization directions are the fast and slow principal axes of the formation, respectively. The shear wave splitting equations for the fast and slow shear waves generated by the shear wave splitting are: ; In the formula, This represents the fast transverse wave signal corresponding to the propagation of the sound source in the X direction. This represents the slow transverse wave signal corresponding to the propagation of the sound source in the X direction; Represents the source function of the sound wave signal. Indicates the rotation angle.

[0027] Furthermore, after being generated in the formation, the fast and slow shear waves propagate at slow speeds s1 and s2, respectively. After traveling a certain distance (from the source to the receiver), they are received by receivers in the X and Y directions, respectively. Projecting the fast and slow waves back to the X and Y directions, respectively, the relationship between the four-component data of the cross-dipole acoustic logging and the fast and slow main curved waves can be simplified to a four-component data equation: ; In the formula, This represents the time-domain signal received by the receiver in the X direction after the sound source emits sound in the X direction; This represents the time-domain signal received by the receiver in the Y direction after the sound source emits sound in the X direction; This represents the time-domain signal received by the receiver in the Y direction after the sound source emits its signal in the Y direction. This represents the time-domain signal received by the receiver in the X direction after the sound source emits sound in the Y direction; This indicates the slowness corresponding to a fast transverse wave; This indicates the slowness of a slow transverse wave. Indicates the distance from the sound source to the receiver; Indicates the rotation angle; This represents the source function of the sound wave signal.

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

[0029] The four-component data equation is further converted into matrix form by rotating it by an angle θ to diagonalize the matrix; the two diagonal elements represent the fast and slow main waves. The resulting core matrix equation is: ; By further utilizing matrix rotation to transform the core matrix equation, the relationship between the four-component data and the main wave and the anisotropic main directions is derived. The specific relationship equation is as follows: ; Furthermore, the transverse wave in the X direction is precisely decomposed into fast transverse wave signals and slow transverse wave signals, and the arrival time of the signal wave at the corresponding receiver is obtained based on the decomposed fast and slow transverse wave signals. and Ultimately, it depends on the distance between the sound source and the receiver. The fast shear wave time difference was calculated separately for the arrival time of the decomposed signal wave at the receiver. and slow transverse wave time difference The formula for calculating time difference is: ; Step S2: Extract the modal components of each signal; within each signal, based on the morphological similarity between different peaks in the frequency domain of each modal component and the frequency domain differences between different modal components, obtain the interference factor of each modal component; based on the interference factor, reconstruct the modal components by weighting and bandpass filtering; use the cross-correlation method to analyze the reconstructed signal to obtain the arrival times of the fast and slow shear waves, and obtain the anisotropy parameters.

[0030] In the process of correcting the in-situ stress of shale reservoirs using the plastic method based on anisotropy, the accuracy of the anisotropy parameter calculation directly affects the accuracy of the in-situ stress correction calculation. In actual processing, fast and slow shear wave signals are affected by wellbore interference, instrument noise, and signal dispersion, which in turn affects the large deviation of the arrival time of the fast and slow shear waves at the corresponding receivers, ultimately affecting the accuracy of the anisotropy parameters.

[0031] Therefore, it is necessary to analyze the interference of signal components, reconstruct the signal and filter it to remove noise, accurately obtain the arrival times of the fast and slow shear waves at the corresponding receivers, improve the accuracy of anisotropic parameter calculation, and provide a more reliable basis for calibration.

[0032] Considering that shear waves will disperse during propagation around the wellbore in actual acquisition, and that instrument noise superposition will affect the signal variation of shear waves in different frequency bands, resulting in jitter and spike interference, the modal components of each signal are extracted first. The modal components correspond to the signals in different frequency bands after signal decomposition. Jitter and spike interference affect the similarity between wave peaks in the frequency domain, which is also reflected in the frequency domain differences of different modal components. Therefore, for each signal, based on the morphological similarity between different wave peaks in the frequency domain of each modal component, combined with the frequency domain differences between different modal components, the interference factor of each modal component is obtained. This characterizes the degree to which each modal component is affected by dispersion and superposition interference, preparing for subsequent weighted reconstruction.

[0033] Preferably, in one embodiment of the present invention, please refer to Figure 3 The diagram illustrates a flowchart of a method for obtaining an interference factor according to an embodiment of the present invention, specifically including: Step S201: Obtain the spectrum of each modal component and extract the peaks. Construct the first feature vector by the peak value, kurtosis and half-peak width of each peak. Obtain the anti-interference coefficient based on the similarity of the first feature vectors among all peaks of each modal component.

[0034] First, the spectrum of each modal component is obtained and the peaks are extracted. Considering that the peak value of each peak reflects the energy intensity of the frequency component, the kurtosis reflects the sharpness of the peak shape, and the half-peak width reflects the distribution range or bandwidth of the peak on the frequency axis, the three reflect the morphological characteristics of the peak from different perspectives. Therefore, the first feature vector is constructed to characterize the morphological characteristics of each peak, providing a basis for subsequent comparison of the morphological similarity between peaks. Considering that the greater the similarity of the first eigenvectors between the peaks within the same modal component, the smaller the impact of jitter and spike interference, the stronger the anti-interference ability of the corresponding modal component, and the larger the anti-interference coefficient.

[0035] As an example, the morphological similarity between different peaks is represented by the cosine similarity between vectors. The anti-interference coefficient is obtained based on the overall characteristics of the cosine similarity of the first eigenvectors of all peak pairs for each modal component. Here, a peak pair is a pair consisting of any two peaks of the same modal component.

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

[0037] It should be noted that the spectrum can be obtained by fast Fourier transform, and taking each local maximum point as the peak point corresponding to each wave peak is a commonly used technique. In other embodiments of the present invention, the implementer can also obtain the anti-interference coefficient by weighted summation of the mean, median and mode of the cosine similarity, such as weights of 0.6, 0.2 and 0.2.

[0038] Step S202: Within the same signal, based on the differences in the distribution characteristics of peak value, kurtosis and half-peak width between each modal component and other modal components, and combined with the anti-interference coefficient, the interference factor of each modal component is obtained.

[0039] Considering the distribution characteristics of modal components in peak value, kurtosis, and half-width at half-maximum, which reflect the stability and purity of signal energy distribution within the frequency band, the degree of "abnormality" or "deviation" of a modal component relative to the overall average characteristics of the signal can be expressed by the differences between the distribution characteristics of each modal component and other modal components in peak value, kurtosis, and half-width at half-maximum within the same signal. The anti-interference coefficient reflects the inherent resistance of the modal component to interference, and thus the interference factor is obtained.

[0040] As an example, each modal component is taken as the target component, and each other modal components within the same signal of the target component are taken as the comparison components; peak value, kurtosis, and half-width are taken as the comparison dimensions; and comparisons are made one by one. A modal component with minimal interference should exhibit concentrated energy (high and prominent peaks), stable shape (moderate and consistent kurtosis), and pure composition (narrow and consistent half-peak width). In contrast, a modal component severely affected by dispersion and noise superposition will have dispersed energy and distorted or broadened peaks, causing these distribution characteristics (such as mean and variance) to differ significantly from other normal components. Therefore, the mean and variance of the eigenvalues ​​of the contrast dimension within each modal component are used to construct a second eigenvector.

[0041] Among them, the feature value is the data value of the data corresponding to the comparison dimension, such as the data value of kurtosis.

[0042] Considering that the larger the Euclidean distance between the second eigenvectors, the more significant the morphological difference between the two modal components in this feature dimension (such as peak distribution), the greater the frequency domain difference, and the larger the interference factor; the larger the sum of the anti-interference coefficients of the two corresponding components, the stronger the resistance to interference in the frequency band where the two modal components are located, and the smaller the interference factor. At the same time, considering that the anti-interference coefficient may be negative, which may affect the relevant logical relationship, a mapping adjustment is performed. Specifically, the Euclidean distance between the target component and the contrast component in the second eigenvector of the contrast dimension is used as the numerator, the sum of the anti-interference coefficients of the two components is mapped by an exponential function with the natural constant e as the base as the denominator, and the ratio of the fractions is used as the interference factor.

[0043] Among them, the exponential function It can map real numbers to the positive number interval, and the denominator will not be negative or zero; is the independent variable.

[0044] The interference sub-factor characterizes the relative intensity of the morphological difference between two components caused by systematic interference. Therefore, the interference factor is obtained by fusing the interference sub-factors of the target component and all contrast components across all contrast dimensions.

[0045] Specifically, the sum of all interfering factors of the target component and other modal components across all contrast dimensions is used as the interfering factor of the target component.

[0046] It should be noted that the analysis process is the same for each signal; only one example is described here, and will not be repeated.

[0047] To effectively reduce the impact of signal quality on the arrival times of fast and slow shear waves in shale reservoirs, the modal components are weighted and reconstructed based on the interference factor and then bandpass filtered to provide a more reliable basis for obtaining anisotropic parameters. The cross-correlation method is then used to analyze the reconstructed signal to obtain the arrival times of fast and slow shear waves and to obtain anisotropic parameters, providing a more reliable basis for the final correction.

[0048] Preferably, in one embodiment of the present invention, considering that the larger the interference factor, the greater the interference of the corresponding modal component to the signal, and the smaller the weight ratio during reconstruction, the interference factor is first negatively correlated within each signal, and the result after Softmax normalization of the negative correlation mapping result is used as the reconstruction weight, and the modal component is reconstructed based on the reconstruction weight.

[0049] As an example, adopt The function is negatively correlated and reconstructed using the weighted envelope average reconstruction method.

[0050] Since the main frequency of the shear wave signal is usually 5~15kHz (shale reservoir), the Buterworth bandpass filter algorithm is used to process the reconstructed signal, with a stop frequency of 4kHz and 16kHz, to further filter out the influence of low-frequency wellbore interference (<4kHz) and high-frequency environmental noise (>16kHz) on signal quality.

[0051] The arrival times of the fast and slow shear waves are further determined 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. Similarly, 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 it takes for a sound wave to travel from its source to the specific receiver (on its first arrival).

[0053] After analyzing the reconstructed signal and obtaining the arrival time of the shear wave, the time difference is substituted into the time difference calculation formula to obtain the accurate fast shear wave time difference. and slow transverse wave time difference Then, the ratio of the slow shear wave time difference to the fast shear wave time difference is used as the anisotropy parameter.

[0054] After obtaining the anisotropic parameters, the process also includes: establishing a model relating the anisotropic parameters to the clay content by combining lithological data.

[0055] Specifically, experimental data on clay content were collected and repositioned. Correlation analysis was performed between clay content and anisotropy parameters calculated from fast and slow shear wave time differences. Clay content was then used to systematically characterize anisotropy parameters, and a relationship model between anisotropy parameters and clay content was established. The calculated clay content-anisotropy relationship is as follows: ; In the formula, Here, represents the anisotropy parameter; a, b, and c are the fitting coefficients (regression coefficients). , , e is the natural constant; This represents the clay content.

[0056] As an example, please see Figure 4 It shows a correlation diagram between clay content and shear wave anisotropy ratio provided by an embodiment of the present invention. Figure 4 The horizontal axis represents clay content, and the vertical axis represents the shear wave anisotropy ratio, i.e., the anisotropy parameter. The data points in the figure are experimental data points, and the curve is the fitted relationship curve between clay content and anisotropy; the corresponding formula is: ; In another embodiment of the invention, normalization is performed first, followed by negative correlation processing. As an example, within each signal, the interference factor is used as input, and Softmax is used for normalization. The processing result is... , Indicates the first The characteristic coefficients of each modal component under the influence of dispersion and superposition interference; the calculation formula for the reconstruction weight includes: ; In the formula, This represents the reconstruction weight for reconstructing the x-th modal component. Indicates the first The characteristic coefficients of each modal component under the influence of dispersion and superposition interference; n represents the total number of modal components.

[0057] In other embodiments of the present invention, a weighted superposition reconstruction method can also be used for reconstruction, which, along with the weighted envelope average reconstruction method and the Buterworth bandpass filtering algorithm, are all well-known techniques; a, b, and c can be: , , a, b, and c can also be: , , I will not go into details.

[0058] Step S3: Based on lithological data and anisotropic parameters, correct the geostress calculation model to obtain the corrected minimum horizontal principal stress value.

[0059] After obtaining accurate anisotropy parameters, the calculation of geostress can be corrected. Geostress calculation requires the combination of lithological data, so the geostress calculation model is finally corrected based on lithological data and anisotropy parameters to obtain the corrected minimum horizontal principal stress value. This provides a reliable basis for accurate drilling trajectory design, optimization of fracturing parameters and efficient reservoir utilization, reduces the blindness of exploration, helps to reduce ineffective drilling, avoid repeated operations, and reduce construction energy consumption, thereby reducing the energy intensity and carbon emission intensity per unit of production from the source.

[0060] Preferably, in one embodiment of the present invention, the correction process includes: First, the anisotropy of resistivity and sonic transit time in horizontal wells is corrected based on anisotropic parameters to eliminate the logging response differences caused by formation bedding. As an example, the relationship between resistivity of vertical and horizontal wells is established based on the anisotropic index curve, the resistivity anisotropy coefficient is obtained, and a 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: ; In the formula, Pore ​​pressure, in MPa; This represents the pressure of the overlying strata, expressed in MPa. This represents the hydrostatic pressure of the formation water column, in MPa. The acoustic transit time at normal compaction pressure is the measured acoustic transit time, expressed in μs / ft. The acoustic transit time of the target layer pressure is corrected using anisotropy parameters, with units of μs / ft; r is the regional Eaton power function, with a default value of r=5.

[0064] The target segment is the current analysis segment. The correction formula is: ; In the formula, For P-wave time difference, corresponding The original value; It is based on The constructed correction factor function is used to eliminate the propagation time error caused by anisotropy, thereby ensuring that the calculation results of pore pressure are more accurate and reliable.

[0065] The correction factor function is generally a proportional or logarithmic relationship, obtained by fitting actual data. Specifically, the actual sonic transit time of vertical wells in the target area and the measured sonic transit time of horizontal wells, as well as the corresponding anisotropy parameters, are collected. The function form and corresponding fitting constant are obtained by fitting the data using statistical regression methods such as linear and exponential regression. After error testing and optimization (selecting the method with the smallest fitting residual), the corresponding correction factor function is obtained. Function fitting is a well-known technique in the art and will not be described in detail here.

[0066] After obtaining the formation pore pressure, a mechanical relationship model between the vertical principal stress and the minimum horizontal principal stress is established using viscoplastic stress relaxation theory to characterize the evolution of stress over time. The mechanical relationship model includes: ; In the formula, is the time-varying differential stress, measured in MPa, reflecting the stress relaxation effect; B is a measure of the rock's elastic flexibility; n is the stress relaxation index, describing the deformation trend over time. It is the total strain. The time of burial depth is expressed in hundreds of millions of years.

[0067] Applying this constitutive law to reservoirs, we can further obtain: ; In the formula, The vertical principal stress is expressed in MPa. The minimum horizontal principal stress is expressed in MPa. It is Young's modulus, corrected for by anisotropy parameters, and is expressed in MPa. It reflects the elastic response characteristics of the formation.

[0068] Furthermore, based on the Young's modulus, formation pore pressure, and formation burial time calculated from the sonic transit time of the corrected horizontal well, and combined with the mechanical relationship model, the minimum horizontal principal stress value was calculated.

[0069] The formulas for calculating the minimum horizontal principal stress include: ; Where Y is the regional coefficient, used to quantify the influence of regional geological background on the minimum horizontal principal stress; Z is the correction coefficient, used to correct systematic errors in the formula derivation process.

[0070] In one embodiment of the present invention, Y takes the value of 0.8 to 1 and Z takes the value of -1.2 to -0.8; as an example, Y takes the value of 0.9 and Z takes the value of -1; in other embodiments of the present invention, the implementer may adjust them as appropriate.

[0071] The correction process fully considers the anisotropic characteristics and stress relaxation effect of shale reservoirs, and realizes dynamic correction and accurate inversion of horizontal well stress parameters, thereby significantly improving the accuracy and reliability of stress prediction.

[0072] It should be noted that the lithological data includes the formation characteristics data required for calculation, such as clay content, formation burial time, and formation water hydrostatic pressure; the models and algorithms used in the calibration process are well-known technologies and will not be described in detail here.

[0073] In another embodiment of the present invention, after correcting the resistivity and sonic transit time of the horizontal well, the method further includes verifying the characterization effect of the correction by comparing the logging curve data of the same formation and lithology of the vertical well.

[0074] After experimental verification, it was proven that the corrected formation logging characteristics can be obtained through anisotropic curves. The anisotropic parameters have a high degree of agreement with the experimental data, proving that anisotropic curves are feasible.

[0075] In summary, to address the technical problem of significant errors in existing shale reservoir in-situ stress analysis due to the influence of complex shale reservoir structures and strong anisotropy, this invention provides an anisotropic shale reservoir plasticity-based in-situ stress correction method. This invention first acquires fast and slow shear wave signals and extracts the modal components of each signal. Further, within each signal, based on the morphological similarity between different peaks in the frequency domain of each modal component and the frequency domain differences between different modal components, an interference factor is obtained. Further, based on the interference factor, the modal components are weighted and reconstructed, and bandpass filtered. Further, the arrival times of the reconstructed fast and slow shear waves are obtained using the cross-correlation method, and anisotropy parameters are acquired. Finally, based on lithological data and anisotropy parameters, the in-situ stress calculation model is corrected to obtain the corrected minimum horizontal principal stress value. This scheme performs modal decomposition and interference-free weighted reconstruction on the array acoustic signal, and uses the cross-correlation method to accurately extract the arrival times of fast and slow shear waves to calculate anisotropic parameters. Combined with lithological data, the geostress model is corrected to achieve accurate acquisition of the minimum horizontal principal stress. This effectively solves the problem of stress calculation deviation caused by anisotropy in complex formations, and provides a reliable basis for accurate drilling trajectory design, optimization of fracturing parameters, and efficient reservoir utilization.

[0076] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0077] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for in-situ stress correction in anisotropic shale reservoirs using the plasticity method, characterized in that, The method includes: Acquire lithological data and well logging data; perform signal decomposition on the array acoustic data contained in the well logging data to obtain fast shear wave signals and slow shear wave signals; Extract the modal components of each signal; within each signal, based on the morphological similarity between different peaks in the frequency domain of each modal component and the frequency domain differences between different modal components, obtain the interference factor of each modal component; based on the interference factor, reconstruct the modal components by weighting and bandpass filtering; use the cross-correlation method to analyze the reconstructed signal to obtain the arrival times of the fast and slow shear waves, and obtain the anisotropy parameters; Based on the lithological data and the anisotropic parameters, the geostress calculation model is modified to obtain the corrected minimum horizontal principal stress value.

2. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, The method for obtaining the interference factor includes: The spectrum of each modal component is obtained and the peaks are extracted. The peak value, kurtosis and half-peak width of each peak are used to form a first feature vector. The anti-interference coefficient is obtained based on the similarity of the first feature vectors between all the peaks of each modal component. Within the same signal, based on 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 respectively, and in conjunction with the anti-interference coefficient, the interference factor of each modal component is obtained.

3. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 2, characterized in that, The method for obtaining the anti-interference coefficient includes: The anti-interference coefficient is obtained based on the overall characteristics of the cosine similarity of the first feature vectors of all the peak pairs of each modal component.

4. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 2, characterized in that, The method for obtaining the interference factor of each of the modal components includes: Each modal component is taken as a target component, and each other modal components within the same signal of the target component are taken as comparison components; the peak value, the kurtosis, and the half-peak width are taken as comparison dimensions; the mean and variance of the eigenvalues ​​of the comparison dimensions within each modal component are used to form a second feature vector; The Euclidean distance between the target component and the comparison component in the second feature vector of the comparison dimension is used as the numerator, the sum of the anti-interference coefficients of the two components is mapped by an exponential function with the natural constant e as the base as the denominator, and the ratio of the fractions is used as the interference factor. The interference factor is obtained by fusing the target component with the interference sub-factors of all the contrast components across all the contrast dimensions.

5. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, The method for weighted reconstruction of the modal components based on the interference factor includes: Within each signal, the negative correlation mapping result of the interference factor is normalized by Softmax and used as the reconstruction weight. The modal components are then reconstructed based on the reconstruction weight.

6. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, The method for obtaining the minimum horizontal principal stress value includes: The anisotropy of resistivity and sonic transit time in horizontal wells is corrected based on the aforementioned anisotropy parameters. The formation pore pressure is obtained by processing the acoustic transit time of overlying strata pressure, normal compaction pressure, formation water hydrostatic column pressure, and target stratum pressure after anisotropic parameter correction using 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.

7. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 6, characterized in that, The method for correcting the anisotropy of resistivity and acoustic time difference 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.

8. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, 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.

9. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, The cutoff frequencies for bandpass filtering are 4kHz and 16kHz.

10. The method for correcting in-situ stress in anisotropic shale reservoirs using the plasticity method according to claim 1, characterized in that, 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.

Images