A method and device for determining fracture density from pure shear reflection amplitude domain information

Abstract

This specification provides a method and apparatus for determining fracture density based on pure shear wave reflection amplitude domain information. The method involves acquiring first and second shear wave data, as well as logging data, for the target area; preprocessing the first and second shear wave data and logging data to determine the first and second shear wave amplitude gathers, seismic wavelets, and stratigraphic interpretation information; simultaneously inverting the first and second shear wave amplitude gathers, seismic wavelets, and stratigraphic interpretation information using a pre-defined pre-stack joint inversion objective function to determine the shear wave splitting parameters; the pre-defined pre-stack joint inversion objective function is constructed based on a linearized reconstruction reflection coefficient equation; and determining the fracture density of the target area based on the shear wave splitting parameters and the pre-defined mapping relationship between the shear wave splitting parameters and fracture density. This method effectively solves the problem of weak capture and low resolution of fracture details by time-difference attributes in existing technologies.

Description

Technical Field

[0001] This manual belongs to the field of oil and gas geological exploration technology, and in particular relates to a method and device for determining fracture density based on pure shear wave reflection amplitude domain information. Background Technology

[0002] Fracture density prediction is an important task in oil and gas reservoir evaluation. However, existing technologies mainly rely on the propagation time difference of fast and slow shear waves for prediction, which has the problems of weak ability to capture fracture details and low resolution, making it difficult to achieve precise quantitative prediction of fracture density.

[0003] There is currently no effective solution to the above problems. Summary of the Invention

[0004] This specification provides a method and apparatus for determining crack density based on pure transverse wave reflection amplitude domain information, which solves the technical problem that existing technologies mainly rely on the propagation time difference of fast and slow transverse waves for prediction, resulting in weak ability to capture crack details and low resolution.

[0005] This specification provides a method for determining crack density based on pure transverse wave reflection amplitude domain information, including: Acquire first shear wave data, second shear wave data, and logging data for the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively; Data preprocessing is performed on the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information; Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information are synchronously inverted to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; The crack density of the target region is determined based on the shear wave splitting parameters, the preset linear mapping relationship between the shear wave splitting parameters and the crack density.

[0006] In one embodiment, the data preprocessing of the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and stratigraphic interpretation information includes: Pre-stack processing is performed on the first shear wave data and the second shear wave data to determine the first shear wave amplitude gather and the second shear wave amplitude gather; The well logging data is used to perform well-seismic calibration on the first and second shear wave amplitude gathers to determine seismic wavelets and layer interpretation information.

[0007] In one embodiment, the step of simultaneously inverting the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters includes: Based on the logging data and the stratigraphic interpretation information, an initial background parameter volume containing information outside the seismic frequency band is determined; Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume are synchronously inverted to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region. The shear wave splitting parameters are determined based on the first shear wave velocity, the second shear wave velocity, and the density.

[0008] In one embodiment, determining the initial background parameter volume containing out-of-seismic-band information based on the logging data and the stratigraphic interpretation information includes: Using the P-wave velocity curve, S-wave velocity curve, and density curve in the logging data, spatial interpolation is performed within the structural framework constrained by the stratigraphic interpretation information to construct a spatially continuous initial parameter field. By performing frequency filtering on the initial parameter field, low-frequency trend components and high-frequency detail components outside the effective seismic frequency band are extracted. The low-frequency trend component and the high-frequency detail component are spectrally reconstructed to determine the initial background parameter volume. In one embodiment, the simultaneous inversion of the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume using a preset pre-stack joint inversion objective function to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region includes: Using the linearized and reconstructed reflection coefficient equation and the seismic wavelet, a mapping operator is constructed; wherein, the mapping operator is used to characterize the mapping relationship between subsurface elastic parameters and seismic response; The data residual term is constructed based on the mapping operator, and the prior constraint term is constructed using the initial background parameter volume. By combining the data residual term and the prior constraint term, the preset pre-stack joint inversion objective function is obtained. Using the preset pre-stack joint inversion objective function, the first shear wave amplitude gather and the second shear wave amplitude gather are iteratively solved to determine the parameter update amount; The initial background parameter volume is corrected using the parameter update amount until the preset convergence condition is met, thereby obtaining the first and second shear wave velocities of the target region.

[0009] In one embodiment, determining the shear wave splitting parameters based on the first shear wave velocity, the second shear wave velocity, and the density includes: The square ratio term is determined based on the ratio of the first shear wave velocity to the second shear wave velocity; Subtracting the preset constant correction term from the square ratio term yields the intermediate difference term; The intermediate difference term is scaled using a preset scaling factor to determine the shear wave splitting parameters of the target region.

[0010] In one embodiment, the method further includes: Based on the crack density, determine the spatial distribution characteristics of crack development intensity within the target area; Based on the spatial distribution characteristics, the exploration potential of the target area is determined.

[0011] This specification provides a crack density determination device based on pure transverse wave reflection amplitude domain information, comprising: The data acquisition module is used to acquire first shear wave data, second shear wave data, and logging data of the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively. The data processing module is used to preprocess the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information. The parameter determination module is used to simultaneously invert the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; The density determination module is used to determine the crack density of the target area based on the shear wave splitting parameters and the preset linear mapping relationship between the shear wave splitting parameters and crack density.

[0012] This specification also provides a computer-readable storage medium storing computer instructions that, when executed, implement a method for determining crack density based on pure transverse wave reflection amplitude domain information.

[0013] Based on the fracture density determination method provided in this specification using pure shear wave reflection amplitude domain information, the following steps are taken: First shear wave data, second shear wave data, and logging data are acquired for the target area. The first and second shear wave data are fast and slow pre-stack shear wave data acquired from two orthogonal observation azimuths, respectively. Data preprocessing is performed on the first, second, and logging data to determine the first and second shear wave amplitude gathers, seismic wavelets, and stratigraphic interpretation information. A preset pre-stack joint inversion objective function is used to simultaneously invert the first, second, and seismic wave amplitude gathers, the seismic wavelets, and the stratigraphic interpretation information to determine the shear wave splitting parameters. The preset pre-stack joint inversion objective function is constructed based on a linearized reconstruction reflection coefficient equation. The fracture density of the target area is determined according to the linear mapping relationship between the shear wave splitting parameters, the preset shear wave splitting parameters, and the fracture density. In this way, by extracting the first and second shear wave amplitude gathers from pre-stack data of fast and slow shear waves from two orthogonal observation azimuths, and performing synchronous inversion using a pre-set pre-stack joint inversion objective function, the physical basis of crack prediction is shifted from the traditional propagation time difference attribute to the amplitude dynamic attribute, which is more sensitive to anisotropy. This effectively solves the problem of weak crack detail capture and low resolution of the time difference attribute in existing technologies. By directly inverting to determine the shear wave splitting parameters, the prediction resolution of crack density is significantly improved, achieving more precise and refined quantitative analysis of crack density in the target area. Attached Figure Description

[0014] To more clearly illustrate the embodiments of this specification, the accompanying drawings used in the embodiments will be briefly introduced below. The drawings described below are only some embodiments recorded in this specification. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 This is a flowchart illustrating a method for determining crack density based on pure transverse wave reflection amplitude domain information, provided in one embodiment of this specification. Figure 2 This is a schematic diagram of the electronic device structure provided in one embodiment of this specification; Figure 3 This is a schematic diagram of the structural composition of a crack density determination device for pure transverse wave reflection amplitude domain information, provided in one embodiment of this specification. Figure 4 This is a schematic diagram of fast and slow shear wave excitation and reception under two special observation azimuths, provided as an embodiment of this specification; Figure 5This is a schematic diagram of pre-stack joint inversion of fast and slow transverse wave reflection amplitudes under two special observation azimuths, provided by one embodiment of this specification. Detailed Implementation

[0016] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this specification, and not all embodiments. Based on the embodiments in this specification, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this specification.

[0017] Fracture prediction is a key task in evaluating sweet spots in hydrocarbon reservoirs. In horizontally symmetric, transversely isotropic (HTI) media induced by directional vertical fractures, shear waves (SWs) exhibit superior potential for detailed characterization compared to PWs due to their high sensitivity to anisotropy and the fast and slow SW characteristics generated by splitting. Addressing the limitations of traditional SW splitting time-of-flight analysis techniques, which suffer from low resolution and difficulty in achieving fine characterization, this invention leverages the development of high-quality pure SW seismic acquisition technology. It fully utilizes the physical characteristic that SWs do not split in two specific observation azimuths: the isotropic plane (along the fracture strike) and the symmetry plane (perpendicular to the fracture strike). By comparing and analyzing the differences in reflection coefficients of fast and slow SWs dominated by SW splitting parameters in these azimuths, this invention constructs a pre-stack joint inversion method based on four-component fast and slow SW data. This method achieves high-precision estimation of SW splitting parameters and, combined with rock physics mapping relationships, successfully achieves refined quantitative prediction of fracture density, providing higher-resolution technical support for the description of complex fractured reservoirs.

[0018] To address the root cause of the aforementioned problems, this specification extracts first and second shear wave amplitude gathers using pre-stack data of fast and slow shear waves from two orthogonal observation azimuths. Simultaneous inversion is then performed using a pre-defined pre-stack joint inversion objective function. This shifts the physical basis of crack prediction from the traditional propagation time difference attribute to the amplitude dynamic attribute, which is more sensitive to anisotropy. This effectively solves the problem of weak crack detail capture and low resolution of existing technologies using time difference attributes. By directly inverting and determining shear wave splitting parameters, the prediction resolution of crack density is significantly improved, achieving more precise and refined quantitative analysis of crack density in the target area.

[0019] See Figure 1 As shown in the embodiments of this specification, a method for determining crack density based on pure transverse wave reflection amplitude domain information is provided, wherein the method is specifically applied to the server side. In specific implementation, the method may include the following: S101: Acquire the first shear wave data, the second shear wave data, and logging data of the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively. S102: Perform data preprocessing on the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information; S103: Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information are synchronously inverted to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; S104: Determine the crack density of the target region based on the shear wave splitting parameters, the preset linear mapping relationship between the shear wave splitting parameters and the crack density.

[0020] The aforementioned first and second shear wave data refer to two sets of seismic wave information with orthogonal polarization characteristics collected for anisotropic underground media. They correspond to the fast shear wave component propagating along the fracture direction and the slow shear wave component propagating perpendicular to the fracture direction, respectively, and are the basic observation data for quantitatively characterizing the shear wave splitting phenomenon.

[0021] The aforementioned logging data can refer to records of underground rock physical parameters obtained through drilling engineering, including sonic transit time, density, imaging logging, and logging evaluation data.

[0022] The two orthogonal observation azimuths mentioned above can refer to two mutually perpendicular azimuth intervals preset according to regional geological stress or fracture orientation during seismic data processing. By extracting azimuth information within these intervals, fast and slow shear wave seismic information with different degrees of anisotropy influence can be separated.

[0023] The aforementioned first shear wave amplitude gather and the aforementioned second shear wave amplitude gather can refer to the set of reflection amplitudes obtained after fidelity processing, dynamic correction and azimuth extraction of the original seismic record. They reflect the variation characteristics of seismic wave intensity under different offsets (angles) and are the direct input carriers for performing pre-stack synchronous inversion.

[0024] The aforementioned seismic wavelet can refer to the source waveform with specific frequency, amplitude, and phase characteristics extracted from seismic wave records, which serves as an important component of the forward modeling operator during the inversion process.

[0025] The aforementioned stratigraphic interpretation information may refer to stratigraphic interface information determined using the continuity characteristics of seismic profiles, which is used to spatially constrain the scope of inversion calculations.

[0026] The aforementioned pre-stack joint inversion objective function can refer to a mathematical functional used to measure the degree of matching between theoretically synthesized seismic data and measured amplitude gathers. It is usually constructed based on a Bayesian theoretical framework, integrating residual terms and constraint terms to guide the inversion algorithm to converge toward the optimal shear wave velocity parameter value.

[0027] The aforementioned transverse wave splitting parameters can refer to physical quantities that reflect the velocity differences and time delays generated when shear waves propagate in anisotropic media.

[0028] The aforementioned pre-defined linear mapping relationship can refer to a pre-determined conversion rule that characterizes the relationship between anisotropic physical fields and geological engineering indicators. It uses linear mapping coefficients to directly map shear wave splitting parameters to fracture density, thus realizing the quantitative conversion from physical parameters to evaluation indicators.

[0029] The linearized reflection coefficient equation described above can be based on HTI medium theory. Under two special orthogonal orientations, namely the isotropic surface and the symmetry axis surface, the conventional reflection coefficient approximation formula is linearized and the parameters are reconstructed. This simplifies the formula to a mapping form that includes only three core variables: fast shear wave velocity, slow shear wave velocity, and density. This significantly enhances the stability and accuracy of extracting crack density-sensitive parameters while reducing the dimensionality of the inversion parameters.

[0030] In some embodiments, acquiring the first shear wave data, the second shear wave data, and the logging data of the target area may specifically include: Using the acquired imaging logging data and geological data of the target area, the fracture development orientation of the strata within the target area is determined; Based on the development direction of the crack, a first observation azimuth and a second observation azimuth are determined, wherein the first observation azimuth is parallel to the development direction of the crack, and the second observation azimuth is perpendicular to the development direction of the crack. The first shear wave data is obtained by extracting the original seismic record of the target area by azimuth using the first observation azimuth. The original seismic record is extracted by azimuth using the second observation azimuth to obtain the second shear wave data.

[0031] It should be noted that conventional shear wave prediction methods usually use all-around seismic records or random azimuth gathers without crack orientation correction for inversion, and rely heavily on P-wave data to indirectly derive shear wave characteristics. This results in the inability to effectively separate the anisotropic information of fast and slow shear waves, and the prediction results often have the defects of strong multiple solutions, low spatial resolution, and inability to quantitatively characterize crack strength.

[0032] This manual pre-determines the fracture development orientation using imaging logging or geological data, and constructs mutually orthogonal observation azimuths accordingly. This enables the physical separation of fast and slow shear wave dynamics during the seismic record extraction stage. This data acquisition method, guided by prior geological information, not only eliminates interference from azimuth stray signals but also provides independent observations with clear physical meaning for subsequent synchronous inversion. It significantly improves the extraction accuracy of shear wave splitting parameters and lays a reliable data foundation for the quantitative and high-precision characterization of fracture density.

[0033] It should be noted that when parallel vertical fissures exist underground, seismic waves undergo shear wave splitting. In this case, the seismic waves propagating along the fissure direction (parallel) travel the fastest and are called fast shear waves; while the seismic waves propagating perpendicular to the fissure direction are most obstructed by the fissure and travel the slowest, and are called slow shear waves. The term "orthogonal" in this instruction manual refers to one orientation being locked "parallel to the fissure direction" and the other orientation being locked "perpendicular to the fissure direction."

[0034] By employing two orthogonal observation azimuths for data extraction, the fast and slow shear wave signals were completely decoupled and purified at the physical level, effectively eliminating signal interference caused by azimuth aliasing in conventional all-around seismic records. This directional extraction method can obtain a pure wavefield that reflects the true differences in stratigraphic anisotropy, providing an independent and reliable observational basis for subsequent inversion and calculation of high-precision shear wave splitting parameters, thereby elevating the stratigraphic fracture development characteristics from fuzzy qualitative estimation to highly reliable quantitative characterization.

[0035] In some embodiments, the linearized reconstructed reflection coefficient equation can be obtained according to the following steps, which may include, in specific implementation: First, we introduce approximate reflection coefficient equations for fast and slow shear waves, characterized by specific anisotropic parameters, under the observation azimuth of the fracture surface. These initial equations establish a physical relationship between the reflection coefficients and formation elastic parameters (including shear wave impedance, shear wave velocity, and density). However, since the parameters to be inverted in these equations often exist in complex nonlinear coupling forms, directly using them as the basis for inversion would significantly increase computational complexity, and the inversion process would be susceptible to noise interference, leading to distorted results.

[0036] In the linearization expansion stage, to address the computational challenges posed by nonlinearity, linear approximations were applied to the parameter impedance, velocity, and density terms in the aforementioned equations. Specifically, the logarithmic difference approximation method was used to rewrite the originally proportional elastic parameter perturbation terms in the equations as differences in logarithmic physical quantities between adjacent stratum sampling points. Through this mathematical transformation, the originally nonlinear reflection coefficient equation was converted into a linear algebraic form concerning the logarithmic elastic parameter perturbation, thus providing theoretical support for achieving fast and stable linear inversion.

[0037] In the parameter reconstruction and simplification stage, to further improve the stability of the inversion and reduce the dimensionality of the parameters to be solved, this embodiment reconstructs the linearized equations based on the intrinsic correlation of physical quantities. Utilizing the additive mapping relationship of shear wave impedance, shear wave velocity, and density in the logarithmic domain, the originally independent logarithmic impedance perturbation terms in the equations are split and merged, uniformly transforming them into a form characterized by the three core physical quantities: fast shear wave velocity, slow shear wave velocity, and density. Through this reconstruction step, the dimensionality reduction and decoupling of the parameters to be inverted are successfully achieved, making the final constructed equations not only possess clear physical meaning but also easier to extract key anisotropic features from seismic data.

[0038] After the above linearization derivation and parameter reconstruction, the specific form of the linearized and reconstructed reflection coefficient equation is as follows:

[0039] in, Angle of incidence To achieve faster transverse wave speed, For slow shear wave velocity, For density.

[0040] It should be noted that existing techniques (such as the classic Rüger approximation formula) generally employ "linearization" methods, but their essence is "mathematical approximation" rather than the "physical reconstruction" of this application. Traditional linearized equations typically retain P-wave velocity, S-wave velocity, density, and three anisotropic parameters, resulting in high dimensionality of the parameters to be determined (usually 5 to 6), strong coupling interference between the parameters, and poor inversion stability.

[0041] This application differs in that it utilizes the physical characteristics of two special observation orientations—the "symmetry axis plane" and the "isotropic plane"—and reconstructs the equations into a simplified form containing only the first shear wave velocity, the second shear wave velocity, and density by substituting the impedance terms. This is not only a mathematical simplification but also a "parameter decoupling" based on physical constraints, which is currently key to achieving high-precision shear wave inversion in the industry. This eliminates insensitive P-wave terms and redundant anisotropic terms, removing "parameter crosstalk" and making the inversion results more unique and stable. Since the equations directly establish a mapping between seismic amplitude and fast and slow shear wave velocities, the cumbersome anisotropic parameter conversion process in traditional methods is avoided, significantly reducing accumulated errors and thus achieving a refined characterization of fracture density. Three-parameter inversion has a much lower matrix computational load than six-parameter inversion, making it more suitable for industrial applications in large-scale work areas.

[0042] Furthermore, the linearized reconstruction reflection coefficient equation in this application is not a general mathematical transformation, but an innovative physical model deeply bound to a specific orthogonal observation azimuth (symmetry axis and isotropic surface). Without this specific observation azimuth, the shear wave will undergo complex splitting, and the three-parameter equation will fail due to the instability of its physical premises. Therefore, the combination of a specific acquisition layout and the reconstruction equation makes a significant technical contribution by decoupling key parameters and significantly improving prediction accuracy, possessing outstanding substantive characteristics.

[0043] In some embodiments, two sets of comparative schemes are used to process synthetic seismic gathers containing random noise. The first set uses conventional linear approximation equations containing P-wave velocity and multiple anisotropic parameters from existing technologies for six-parameter synchronous inversion. The results show that under noise interference, the predicted values ​​of S-wave velocity and anisotropic parameters exhibit significant oscillations and large deviations from well logging data, indicating that high-dimensional parameter coupling enhances the instability of the inversion. The second set uses the reconstruction equations proposed in this application, which only include fast S-wave velocity (first S-wave velocity), slow S-wave velocity (second S-wave velocity), and density, for three-parameter joint synchronous inversion. Under the same noise level, the inverted fast and slow S-wave velocity curves maintain extremely high consistency with the well logging results, and the computation time is only one-third of the first set of schemes. These experimental results strongly demonstrate that by reconstructing and decoupling the reflection coefficient equations, this application not only significantly reduces the ambiguity and computational cost of inversion calculations but also achieves unexpected technological advancements in noise resistance and fracture density prediction accuracy.

[0044] In some embodiments, determining the crack density of the target region based on the shear wave splitting parameters, a preset linear mapping relationship between the shear wave splitting parameters and crack density may include, in practice: The inverted shear wave splitting parameters are substituted into a preset linear mapping relationship between shear wave splitting parameters and crack density for numerical mapping processing to obtain the crack density of the target region. The linear mapping relationship between the preset shear wave splitting parameter and the fracture density is a linear proportional function with the shear wave splitting parameter as the independent variable and the fracture density as the dependent variable, and the linear proportional function has a preset linear mapping coefficient. The linear mapping coefficient is a constant predetermined based on the equivalent correlation characteristics between the physical anisotropy intensity of the strata and the fracture development density in the target area, and is used to convert the shear wave splitting parameter into the fracture density in terms of magnitude.

[0045] In some embodiments, the linear mapping relationship between the preset shear wave splitting parameters and the crack density can be obtained by the following method, which may include: First, the shear wave splitting parameters are estimated using the fast and slow shear wave synchronous inversion method as described above. To achieve a quantitative assessment of subsurface fracture development, a mapping relationship needs to be established between these splitting parameters and geologically significant fracture density. A classic rock physics model is introduced, establishing the shear wave splitting parameters as a function of fracture density, Poisson's ratio, and proportionality coefficient.

[0046] Considering the statistical regularity of the physical properties of underground reservoir rocks, a constraint analysis was performed on the Poisson's ratio term in the model. Based on actual well logging data and geological experience, the Poisson's ratio value of rocks in the formation is usually stable within a specific range. Sensitivity analysis within this range revealed that the coefficient term composed of Poisson's ratio has extremely high stability within this range, and its calculation result always approaches the unit value of 1.

[0047] Based on the aforementioned physical understanding that coefficient terms approach 1, the complex rock physics equations were linearly simplified. By ignoring the influence of coefficients with minimal fluctuations, the original proportional relationships were simplified to approximately equivalent relationships. This processing step shows that the inverted shear wave splitting parameters are numerically equivalent to fracture density. Through the establishment of this linear mapping relationship, a direct conversion from abstract seismic anisotropy parameters to intuitive geological fracture density parameters was achieved, greatly simplifying the calculation process for subsequent quantitative prediction of reservoir fractures.

[0048] Specifically, the linear mapping relationship between the preset shear wave splitting parameters and the crack density can be determined according to the following formula:

[0049] in, The shear wave splitting parameters are... Poisson's ratio, For linear mapping coefficients, The crack density is given.

[0050] The Poisson's ratio of rocks in the strata is generally between 0.2 and 0.4, corresponding to the coefficient term 8(1- in the above equation). υ ) / 3(2- υ Since )≈1, the relationship between the shear wave splitting parameter and the crack density can be approximately expressed as:

[0051] The above equation shows that the transverse wave splitting parameter is equivalent to the crack density.

[0052] In some embodiments, the data preprocessing of the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and stratigraphic interpretation information may further include the following in specific implementations: S1: Perform pre-stack processing on the first shear wave data and the second shear wave data to determine the first shear wave amplitude gather and the second shear wave amplitude gather; S2: Use the logging data to perform well-seismic calibration on the first shear wave amplitude gather and the second shear wave amplitude gather to determine the seismic wavelet and layer interpretation information.

[0053] Specifically, the first and second shear wave data are first acquired using multi-wave seismic exploration equipment. To address clutter interference and energy attenuation in the raw data, this embodiment employs a pre-stack fidelity processing workflow. Specifically, the first and second shear wave data undergo surface consistency amplitude compensation, pre-stack random noise attenuation, and spherical diffusion compensation sequentially to eliminate the influence of near-surface conditions and propagation paths on seismic amplitude. Subsequently, precise dynamic correction (NMO) and gather removal processes are used to remove stretching distortion caused by large incident angles, ultimately forming high-fidelity, high-signal-to-noise ratio first and second shear wave amplitude gathers, providing a highly reliable foundation for subsequent splitting parameter inversion.

[0054] We fully utilize well logging data (such as sonic transit time, density, and shear wave velocity) within the study area as a "geological benchmark." First, we establish accurate time-depth conversion relationships using well logging curves through first-arrival matching or synthetic recording, ensuring a high degree of agreement between seismic and geological strata in both the depth and time domains. Based on this, we employ multi-well constrained statistical or deterministic extraction methods to extract representative seismic wavelets from the first and second shear wave amplitude collections. These seismic wavelets not only contain the temporal characteristics of the source but also reflect the absorption and compensation properties of the subsurface medium, making them a core element in constructing the inversion-forward modeling operator.

[0055] By combining the layering information provided by the logging data, manual or semi-automatic stratigraphic tracking is performed on the preprocessed seismic profile to determine the stratigraphic interpretation information of the target reservoir section. This stratigraphic interpretation information not only defines the time window range for inversion but also provides spatial lateral continuity constraints for the inversion process. Through the above preprocessing steps, discrete logging point information is organically combined with continuous seismic surface information, outputting four core parameters that meet the inversion requirements: high-quality two-component shear wave amplitude gathers, accurate seismic wavelets, and stratigraphic interpretation information that reflects the geological structural framework, thereby significantly reducing the ambiguity of the inversion process.

[0056] In some embodiments, the method of simultaneously inverting the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters may further include the following: S1: Based on the logging data and the stratigraphic interpretation information, determine the initial background parameter volume containing information outside the seismic frequency band; S2: Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume are synchronously inverted to determine the first shear wave velocity, the second shear wave velocity, and the density of the target area. S3: Determine the shear wave splitting parameters based on the first shear wave velocity, the second shear wave velocity, and the density.

[0057] Specifically, firstly, using the collected well logging data and stratigraphic interpretation information of the target area, a spatial low-frequency background framework is established by converting the well logging curves from the depth domain to the time domain and combining the stratigraphic interpretation results. Based on this, spatial interpolation and frequency band fusion processing are performed on the well point data to obtain an initial background parameter volume covering low-frequency trends and high-frequency characteristics outside the effective seismic frequency band, providing stable initial values ​​and background references for subsequent inversion iterations.

[0058] Subsequently, the obtained initial background parameter volume, the first and second shear wave amplitude gathers acquired from different orthogonal observation azimuths, and the seismic wavelet are substituted into the pre-stack joint inversion objective function constructed based on a Bayesian framework. This objective function establishes a mapping relationship between the model and the data through a linearized reconstructed reflection coefficient equation, decoupling the parameters to be determined into a three-parameter form containing only the first shear wave velocity (fast shear wave velocity), the second shear wave velocity (slow shear wave velocity), and density. The Gauss-Newton iterative algorithm is used to solve the objective function. During the inversion process, the residuals between the synthetic records and the actual gathers are calculated in real time, and iterative corrections are performed in conjunction with the model's prior constraints, thereby obtaining high-resolution inversion results for the fast shear wave velocity, slow shear wave velocity, and density of the target area.

[0059] The obtained fast S&W velocity, slow S&W velocity, and density inversion results are directly used for quantitative estimation of S&W splitting parameters. Through a pre-defined S&W splitting parameter calculation logic, the ratio difference between fast and slow S&W velocities is used to obtain S&W splitting parameters characterizing the anisotropy intensity of the formation. Because a decoupled reconstruction equation is used for simultaneous three-parameter inversion in the inversion stage, the crosstalk effects of P-wave parameters and redundant anisotropy parameters on the inversion process are effectively avoided. This results in more stable fast and slow S&W velocity results, significantly improving the ability of S&W splitting parameters to resolve micro-fractures, and laying a precise physical foundation for subsequent quantitative prediction of fracture density.

[0060] By introducing an initial background parameter volume containing information outside the seismic frequency band, the deficiency of missing low-frequency components in seismic data is effectively compensated, significantly enhancing the convergence and uniqueness of the solution in the inversion process. At the same time, by adopting a pre-stack joint synchronous inversion strategy, the first shear wave velocity, the second shear wave velocity, and the density are extracted in a coordinated manner. This not only improves the calculation accuracy and spatial resolution of shear wave splitting parameters, but also more realistically portrays the anisotropic characteristics of the subsurface medium, providing a more scientific and reliable technical means for describing complex fractured-vuggy reservoirs.

[0061] In some embodiments, the method for determining an initial background parameter volume containing out-of-seismic frequency band information based on the logging data and the stratigraphic interpretation information may further include the following: S1: Using the P-wave velocity curve, S-wave velocity curve and density curve in the logging data, spatial interpolation is performed within the structural framework constrained by the stratigraphic interpretation information to construct a spatially continuous initial parameter field. S2: By performing frequency filtering on the initial parameter field, low-frequency trend components and high-frequency detail components outside the effective seismic frequency band are extracted; S3: Perform spectral reconstruction on the low-frequency trend component and the high-frequency detail component to determine the initial background parameter volume.

[0062] Specifically, the P-wave velocity, S-wave velocity, and density curves of multiple key wells within the target area are first extracted. Using the structural stratigraphic planes determined by the stratigraphic interpretation information as spatial constraints, a three-dimensional structural mesh model is established. Variational function analysis is employed, under the constraints of the structural framework, to perform three-dimensional kriging interpolation or inverse distance weighted interpolation on the discrete well logging curves, thereby expanding the point-like logging information into a spatially continuous initial parameter field. Because this process strictly adheres to the formation attitude and fault morphology, it ensures the continuity of geological regularities in the horizontal direction and retains the high sampling rate characteristics of the logging data in the vertical direction.

[0063] Considering that seismic exploration data typically exhibits a "bandpass" characteristic due to seismic wave attenuation and limitations in excitation and reception conditions, i.e., it lacks extremely low-frequency and high-frequency components, multi-scale frequency filtering is applied to the initial parameter field. A low-pass filter is used to extract low-frequency trend components representing the sedimentary background and crustal evolution trend, while a high-pass filter is used to extract high-frequency detail components exceeding the seismic resolution limit from high-resolution well logging information. This step uses mathematical filtering techniques to decouple the full-band information in the well logging data according to its relationship with the seismic frequency band.

[0064] The extracted low-frequency trend components and high-frequency detail components are linearly superimposed and their spectroscopic reconstructions are performed to determine the final initial background parameter volume. During the reconstruction process, the cutoff frequency and suppression slope of the filter are adjusted to make the reconstructed parameter volume exhibit a "concave" characteristic in its spectral distribution, leaving space for inversion within the effective seismic frequency band. This constructed initial background parameter volume contains both a low-frequency "skeleton" capable of stabilizing the inversion convergence process and high-frequency "details" that enhance stratigraphic characterization, laying a solid foundation for subsequent wideband, high-resolution synchronous inversion.

[0065] By spatially interpolating and recombining the spectral decomposition of logging data, an initial background parameter volume containing information outside the seismic frequency band was constructed, which has the following significant advantages: Enhanced inversion stability: The introduced low-frequency trend component provides correct large-scale constraints for nonlinear inversion, effectively overcoming the multiple solutions in seismic inversion and ensuring the accuracy of the shear wave splitting parameters in absolute terms; Improved vertical resolution: By integrating high-frequency detail features from logging into the initial model through spectral recombination, the shortcomings of seismic records in thin-layer identification are compensated, enabling the finally determined first and second shear wave velocity fields to more delicately characterize the microscopic heterogeneity within the reservoir; Strong geological consistency: The interpolation method based on the stratigraphic constraint structural framework ensures a high degree of consistency between the background model and the actual subsurface structural morphology, avoiding the layer-crossing phenomenon or artificial artifacts caused by simple mathematical interpolation.

[0066] In some embodiments, the method utilizes a preset pre-stack joint inversion objective function to simultaneously invert the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region. In specific implementations, the method may further include the following: S1: Using the linearized and reconstructed reflection coefficient equation and the seismic wavelet, a mapping operator is constructed; wherein, the mapping operator is used to characterize the mapping relationship between subsurface elastic parameters and seismic response; S2: Construct a data residual term based on the mapping operator, and construct a prior constraint term using the initial background parameter volume. By combining the data residual term and the prior constraint term, obtain the preset pre-stack joint inversion objective function. S3: Using the preset pre-stack joint inversion objective function, iteratively solve the first shear wave amplitude gather and the second shear wave amplitude gather to determine the parameter update amount; S4: The initial background parameter volume is corrected using the parameter update amount until the preset convergence condition is met, and the first shear wave velocity and the second shear wave velocity of the target region are obtained.

[0067] Specifically, the linearized and reconstructed reflection coefficient equation determined in the preceding steps is first combined with the extracted seismic wavelet. More specifically, the seismic wavelet is used as a convolution kernel, and spatial convolution is performed with the coefficients of the linearized and reconstructed reflection coefficient equation to construct a three-dimensional matrix mapping operator. This mapping operator physically characterizes the linear mapping relationship between subsurface elastic parameter disturbances (such as changes in velocity and density) and the seismic response observed at the surface (first and second shear wave amplitude gathers). Through this operator, complex subsurface geological features can be projected onto the seismic wavefield space, providing a theoretical prediction basis for subsequent residual calculations.

[0068] The pre-stack joint inversion objective function is constructed by combining two key terms. First, the predicted seismic gathers are calculated based on the mapping operator, and the difference between these gathers and the actually observed first and second shear wave amplitude gathers is calculated to construct a data residual term, which measures the degree of matching between the inversion results and the original seismic data. Second, a priori constraint term (regularization term) is constructed using an initial background parameter volume containing out-of-band information to constrain the spatial continuity of the inversion solution and introduce low-frequency geological trends. By setting weight coefficients, the data residual term and the prior constraint term are weighted and combined to form the final pre-set pre-stack joint inversion objective function, thereby effectively reducing the ambiguity of the inversion while ensuring the fidelity of the seismic data.

[0069] The algorithm employs gradient descent, conjugate gradient, or quasi-Newton methods (such as L-BFGS) to optimize the pre-stack joint inversion objective function. In each iteration, the algorithm first calculates the gradient direction under the current model to determine the parameter update amounts, i.e., the correction values ​​for fast and slow shear wave velocities and densities. Subsequently, the initial background parameter volume is progressively corrected using these parameter update amounts. After each correction, the objective function value is recalculated, and it is determined whether it meets the pre-defined convergence conditions (such as the sum of squared residuals reaching a pre-defined threshold or the improvement between two iterations being less than a given value).

[0070] Once the iterative process meets the convergence condition, the computation stops and the final corrected elastic parameter model is output. At this point, the values ​​in the model represent the first shear wave velocity (fast shear wave velocity) and the second shear wave velocity (slow shear wave velocity) that are highly geologically representative within the target region. Because this process is completed synchronously within a unified objective function framework, it ensures a high degree of coordination between the first and second shear wave velocities in terms of spatial location and dynamic characteristics, laying a reliable foundation for subsequent high-precision calculation of shear wave splitting parameters.

[0071] Based on the above logic, the mathematical expression of the preset pre-stack joint inversion objective function used in this embodiment is as follows:

[0072] in, The pre-stack joint inversion objective function is set as follows: Let be the vector of model parameters to be solved. For the observed data vector, For mapping operators, This is the initial background parameter volume. For regularization parameters, This is the matrix transpose operator.

[0073] In the formula, It can be represented as

[0074] in, It is the seismic wavelet matrix. For difference operator matrices, It is a coefficient matrix.

[0075] To obtain the model parameter values ​​through inversion, the functional J is partially derived with respect to the model parameter m, and... J / When m = 0, the update amount of the model parameters can be obtained as follows:

[0076] in, This represents the identity matrix. The model parameters can be obtained iteratively using the Gauss-Newton method.

[0077] Here, k represents the number of iterations, and its selection is highly dependent on factors such as model complexity and the signal-to-noise ratio of the data. In actual data inversion, a suitable number of iterations can be set by comparing the inversion results of the well bypass with the actual logging data.

[0078] In some embodiments, the step of correcting the initial background parameter volume using the parameter update amount until a preset convergence condition is met, thereby obtaining the first and second shear wave velocities of the target region, may specifically include: In the first iteration, the initial background parameter volume is used as the current elastic parameter model to be inverted; The current elastic parameter model to be inverted is summed and corrected using the parameter update amount to obtain the updated elastic parameter model to be inverted. Based on the mapping operator, the predicted earthquake record corresponding to the updated elastic parameter model to be inverted is calculated, and the data residual functional value between the predicted earthquake record and the first and second shear wave amplitude gathers is calculated. Determine whether the residual functional value of the data is less than a preset error threshold, or whether the current iteration number has reached the preset maximum iteration number; If any convergence condition is met, the iteration stops, the current updated elastic parameter model to be inverted is determined as the final inversion model, and the corresponding first and second shear wave velocity components are extracted from it. If the convergence condition is not met, the current updated elastic parameter model to be inverted is used as the starting model for the next iteration, and the process returns to the step of determining the parameter update amount.

[0079] Specifically, the initial state of the model to be inverted is first established. During the initial iteration (i.e., iteration counter $k=1$), the initial background parameter volume, which contains low-frequency information outside the seismic frequency band and was constructed in the preceding steps, is read and used as the current elastic parameter model to be inverted. This initial background parameter volume provides the basic velocity framework and density trend of the subsurface medium for inversion, ensuring that the subsequent nonlinear search process is conducted within a geophysically meaningful solution space, thereby effectively avoiding the common multi-solution trap in seismic inversion.

[0080] The current elastic parameter model to be inverted is supplemented and corrected using the currently determined parameter update amounts. Specifically, the update vectors reflecting formation velocity and density perturbation characteristics are superimposed point-by-point onto the original model according to the spatial grid points, thereby generating the updated elastic parameter model to be inverted. Subsequently, a forward modeling simulation is performed on the updated model using a pre-constructed mapping operator (composed of a linearized reconstruction of the reflection coefficient equation and seismic wavelet convolution) to synthesize the predicted seismic record. By calculating the energy difference between this predicted seismic record and the first and second shear wave amplitude gathers observed in actual seismic data, a data residual functional is constructed to quantitatively evaluate the fitting accuracy of the current model to the actual subsurface geological response.

[0081] The algorithm monitors the dynamic changes of the functional values ​​of the data residuals in real time. Two convergence criteria are preset: one is the error threshold criterion, which determines whether the current data residuals have decreased to a preset minimum allowable error range; the other is the scale criterion, which determines whether the current iteration count has reached a preset maximum iteration count (e.g., 50 or 100 iterations). If either of these conditions is met, the inversion is considered to have reached convergence, and iteration stops. At this point, the parameter model obtained from the last update is determined as the final inversion model, and the first shear wave velocity component (fast shear wave velocity) and the second shear wave velocity component (slow shear wave velocity) corresponding to each sampling point are precisely extracted using index mapping or parameter extraction operators.

[0082] If the current iteration result does not yet meet the convergence condition, the program automatically enters the closed-loop control process. The program uses the currently updated elastic parameter model to be inverted as the starting input model for the next iteration and returns to the step of determining the parameter update amount, calculating the new search gradient and iteration step size. Through this continuous "prediction-difference-correction" cycle, the elastic parameter model, driven by seismic data, gradually approximates the actual physical parameter field underground, ultimately achieving a high-precision quantitative characterization of the first and second shear wave velocities in the target area.

[0083] In some embodiments, the method for determining the shear wave splitting parameters based on the first shear wave velocity, the second shear wave velocity, and the density may further include the following: S1: Determine the square ratio term based on the ratio of the first shear wave velocity to the second shear wave velocity; S2: Subtract the preset constant correction term from the square ratio term to obtain the intermediate difference term; S3: Scale the intermediate difference term using a preset scaling factor to determine the shear wave splitting parameters of the target region.

[0084] Specifically, the first shear wave velocity (fast shear wave velocity) and the second shear wave velocity (slow shear wave velocity) obtained from the aforementioned synchronous inversion step are received. The quotient of the first shear wave velocity and the second shear wave velocity is calculated point by point for each sampling point. Then, this ratio is squared to obtain a squared ratio term reflecting the energy distribution of the velocity difference in the two orthogonal directions. Since the changes in fast and slow shear wave velocities are directly affected by the development of formation fractures and the intensity of anisotropy, this squared ratio term can sensitively capture the splitting characteristics of shear waves within the geological parameter space, laying the foundation for subsequent quantitative transformation.

[0085] Considering that the velocities of the two components of a seismic wave should theoretically remain consistent when propagating in an isotropic background medium, a preset constant correction term is introduced to subtract the background value from the square ratio term. Specifically, the constant correction term (e.g., an offset constant set according to the properties of the background medium) is subtracted from the square ratio term to obtain the intermediate difference term. The physical significance of this step is to remove the interference from the stratigraphic background and highlight the velocity asymmetry caused by fractures through differential processing, making the calculation results more focused on the effective anisotropic disturbances.

[0086] The intermediate difference term is linearly scaled using a preset scaling factor to determine the shear wave splitting parameters of the target area. This scaling factor is typically pre-calibrated based on rock physics experimental data or empirical models within the work area, and is used to map the mathematical domain difference to the order of magnitude of physically meaningful shear wave splitting parameters. Through this scaling process, the finally determined shear wave splitting parameters can characterize the development intensity and anisotropy of fractures within the reservoir in a standardized numerical form, providing crucial reference data for subsequent oil and gas reservoir exploration and evaluation.

[0087] The method of determining shear wave splitting parameters through a combination of square ratio calculation, constant correction, and scaling has the following significant advantages: Improved parameter sensitivity: The square ratio processing method amplifies the small velocity differences between fast and slow shear waves, enhancing the ability to identify weak anisotropic features; Strong computational stability: By introducing a constant correction term to subtract background interference, computational redundancy caused by isotropic formations is effectively eliminated, ensuring that the splitting parameters only reflect the effective changes caused by fractures; Wide application adaptability: The introduction of scaling factors allows the algorithm to be flexibly adjusted according to logging data from different work areas, ensuring a high degree of consistency between the inversion results and actual geological conditions, providing a more accurate and reliable technical means for the quantitative description of complex reservoirs.

[0088] In some embodiments, the method may further include the following: S1: Determine the spatial distribution characteristics of crack development intensity within the target area based on the crack density; S2: Determine the exploration potential of the target area based on the spatial distribution characteristics.

[0089] Specifically, the fracture density data of the entire work area obtained from the aforementioned inversion steps are processed into a three-dimensional spatial grid. Using geological modeling software and the structural morphology determined by stratigraphic interpretation, a planar distribution map of fracture density at the target stratigraphic level is extracted. Through spatial trend analysis, high-value clusters and low-value areas of fracture development are identified, and the heterogeneity of fracture development intensity is determined by combining the distribution characteristics of the fault system. Furthermore, the dominant orientation and connectivity characteristics of fracture development are statistically analyzed, and a three-dimensional spatial distribution characteristic model of fracture development intensity within the target area is constructed using spatial kriging interpolation or stochastic simulation methods. This allows for the quantitative and visual representation of fracture density and geometric distribution in three-dimensional space.

[0090] Taking into account the spatial distribution characteristics and reservoir hydrocarbon-bearing patterns, a multi-parameter fusion potential evaluation system is established. Specifically, high-fracture density areas (i.e., areas with high development intensity) are overlaid with reservoir properties, structural high points, and known well gas / oil production for analysis. Based on the regression relationship between fracture development intensity and reservoir productivity, "fracture development sweet spots" are identified as key exploration targets. By comprehensively scoring the fracture connectivity and structural location of different areas, the exploration potential level of each evaluation unit within the target area is determined, and the final exploration potential result map is output. The exploration potential results directly guide subsequent well placement plans and horizontal well fracturing designs, ensuring efficient drilling in areas with the most developed fractures and the greatest potential.

[0091] As can be seen from the above, the embodiment of this specification provides a method for determining fracture density based on pure shear wave reflection amplitude domain information. This method acquires first shear wave data, second shear wave data, and logging data for a target area. The first and second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data, respectively, acquired from two orthogonal observation azimuths. Data preprocessing is performed on the first and second shear wave data and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and stratigraphic interpretation information. A preset pre-stack joint inversion objective function is used to simultaneously invert the first and second shear wave amplitude gathers, the seismic wavelet, and the stratigraphic interpretation information to determine the shear wave splitting parameters. The preset pre-stack joint inversion objective function is constructed based on a linearized reconstruction reflection coefficient equation. The fracture density of the target area is determined according to the shear wave splitting parameters, the preset shear wave splitting parameters, and the linear mapping relationship between the fracture density and the fracture parameters. In this way, by extracting the first and second shear wave amplitude gathers from pre-stack data of fast and slow shear waves from two orthogonal observation azimuths, and performing synchronous inversion using a pre-set pre-stack joint inversion objective function, the physical basis of crack prediction is shifted from the traditional propagation time difference attribute to the amplitude dynamic attribute, which is more sensitive to anisotropy. This effectively solves the problem of weak crack detail capture and low resolution of the time difference attribute in existing technologies. By directly inverting to determine the shear wave splitting parameters, the prediction resolution of crack density is significantly improved, achieving more precise and refined quantitative analysis of crack density in the target area.

[0092] See Figure 2 As shown in the embodiments of this specification, a specific electronic device is also provided, wherein the electronic device includes a network communication port 201, a processor 202 and a memory 203, and the above structures are connected by internal cables so that the various structures can perform specific data interaction.

[0093] Specifically, the network communication port 201 can be used to acquire first shear wave data, second shear wave data, and logging data of the target area; wherein the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively.

[0094] The processor 202 is specifically used to preprocess the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information; using a preset pre-stack joint inversion objective function, it performs synchronous inversion on the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; and according to the linear mapping relationship between the shear wave splitting parameters, the preset shear wave splitting parameters, and the fracture density, it determines the fracture density of the target area.

[0095] The memory 203 can be used to store the corresponding instruction program.

[0096] Based on the above method, the relevant structural performance of electronic devices can be effectively utilized to improve the data processing speed of electronic devices and efficiently realize a method for determining crack density in the pure transverse wave reflection amplitude domain.

[0097] In this embodiment, the network communication port 201 can be a virtual port bound to different communication protocols, thereby enabling the sending or receiving of different data. For example, the network communication port can be a port responsible for web data communication, a port responsible for FTP data communication, or a port responsible for email data communication. Furthermore, the network communication port can also be a physical communication interface or communication chip. For example, it can be a wireless mobile network communication chip, such as GSM or CDMA; it can also be a Wi-Fi chip; or it can be a Bluetooth chip.

[0098] In this embodiment, the processor 202 can be implemented in any suitable manner. For example, the processor can take the form of a microprocessor or processor and a computer-readable medium storing computer-readable program code (e.g., software or firmware) executable by the (micro)processor, logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers, and embedded microcontrollers, etc. This specification is not limiting.

[0099] In this embodiment, the memory 203 may include a hierarchy. In a digital system, anything that can store binary data can be a memory. In an integrated circuit, a circuit with storage function but no physical form is also called a memory, such as RAM, FIFO, etc. In a system, a storage device with a physical form is also called a memory, such as a memory stick, TF card, etc.

[0100] This specification also provides a computer-readable storage medium based on the above-described method for determining fracture density using pure shear wave reflection amplitude domain information. The method acquires first shear wave data, second shear wave data, and logging data for a target area. The first and second shear wave data are fast and slow pre-stack shear wave data acquired from two orthogonal observation azimuths, respectively. Data preprocessing is performed on the first, second, and logging data to determine the first and second shear wave amplitude gathers, seismic wavelets, and stratigraphic interpretation information. A preset pre-stack joint inversion objective function is used to simultaneously invert the first and second shear wave amplitude gathers, the seismic wavelets, and the stratigraphic interpretation information to determine shear wave splitting parameters. The preset pre-stack joint inversion objective function is constructed based on a linearized reconstruction reflection coefficient equation. The fracture density of the target area is determined according to the linear mapping relationship between the shear wave splitting parameters, the preset shear wave splitting parameters, and the fracture density.

[0101] In this embodiment, the storage medium includes, but is not limited to, Random Access Memory (RAM), Read-Only Memory (ROM), Cache, Hard Disk Drive (HDD), or Memory Card. The memory can be used to store computer program instructions. The network communication unit can be an interface configured according to standards specified in the communication protocol for network connection communication.

[0102] In this embodiment, the specific functions and effects implemented by the program instructions stored in the computer-readable storage medium can be explained in comparison with other embodiments, and will not be repeated here.

[0103] See Figure 3 At the software level, this specification also provides a crack density determination device based on pure transverse wave reflection amplitude domain information. This device may specifically include the following structural modules: The data acquisition module 301 is used to acquire first shear wave data, second shear wave data and logging data of the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data respectively acquired based on two orthogonal observation azimuths. Data processing module 302 is used to perform data preprocessing on the first shear wave data, the second shear wave data and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet and the stratigraphic interpretation information; The parameter determination module 303 is used to simultaneously invert the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; The density determination module 304 is used to determine the crack density of the target area based on the shear wave splitting parameters and the preset linear mapping relationship between the shear wave splitting parameters and the crack density.

[0104] In some embodiments, the data processing module 302, in its specific implementation, performs pre-stack processing on the first shear wave data and the second shear wave data to determine the first shear wave amplitude gather and the second shear wave amplitude gather; and uses the logging data to perform well-seismic calibration on the first shear wave amplitude gather and the second shear wave amplitude gather to determine seismic wavelet and layer interpretation information.

[0105] In some embodiments, the parameter determination module 303, specifically implemented as follows: An initial background parameter volume determination module is used to determine an initial background parameter volume containing information outside the seismic frequency band based on the logging data and the stratigraphic interpretation information; using a preset pre-stack joint inversion objective function, synchronously inverting the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region; and a shear wave splitting parameter determination module is used to determine shear wave splitting parameters based on the first shear wave velocity, the second shear wave velocity, and the density.

[0106] In some embodiments, the initial background parameter volume determination module, in its specific implementation, utilizes the P-wave velocity curve, S-wave velocity curve, and density curve in the logging data to perform spatial interpolation within the structural framework constrained by the stratigraphic interpretation information, constructing a spatially continuous initial parameter field; by performing frequency filtering on the initial parameter field, low-frequency trend components and high-frequency detail components outside the effective seismic frequency band are extracted; the low-frequency trend components and the high-frequency detail components are then spectrally reconstructed to determine the initial background parameter volume.

[0107] In some embodiments, the parameter determination module 303, in its specific implementation, constructs a mapping operator using the linearized reconstructed reflection coefficient equation and the seismic wavelet; wherein, the mapping operator is used to characterize the mapping relationship between subsurface elastic parameters and seismic response; a data residual term is constructed based on the mapping operator, and a priori constraint term is constructed using the initial background parameter volume; by combining the data residual term and the prior constraint term, the preset pre-stack joint inversion objective function is obtained; the preset pre-stack joint inversion objective function is used to iteratively solve the first shear wave amplitude gather and the second shear wave amplitude gather to determine the parameter update amount; the parameter update amount is used to correct the initial background parameter volume until the preset convergence condition is met, thereby obtaining the first shear wave velocity and the second shear wave velocity of the target region.

[0108] In some embodiments, the above-mentioned shear wave splitting parameter determination module specifically implements the following steps: determining a square ratio term based on the ratio of the first shear wave velocity to the second shear wave velocity; subtracting a preset constant correction term from the square ratio term to obtain an intermediate difference term; and scaling the intermediate difference term using a preset scaling factor to determine the shear wave splitting parameters of the target region.

[0109] In some embodiments, the above-described apparatus further includes: determining the spatial distribution characteristics of fracture development intensity within the target area based on the fracture density; and determining the exploration potential of the target area based on the spatial distribution characteristics.

[0110] It should be noted that the units, devices, or modules described in the above embodiments can be implemented by computer chips or physical entities, or by products with certain functions. For ease of description, the above devices are described by dividing them into various modules according to their functions. Of course, in implementing this specification, the functions of each module can be implemented in the same software and / or hardware, or modules that implement the same function can be implemented by a combination of sub-modules or sub-units, etc. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and there may be other division methods in actual implementation. For example, units or components can be combined or integrated into another system, or some features can be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection between the devices or units shown or discussed can be through some interfaces, and the indirect coupling or communication connection between devices or units can be electrical, mechanical, or other forms.

[0111] As can be seen from the above, the fracture density determination device based on the pure shear wave reflection amplitude domain information provided in the embodiments of this specification acquires first shear wave data, second shear wave data, and logging data of the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively; the first shear wave data, the second shear wave data, and the logging data are preprocessed to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information; using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information are synchronously inverted to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; the fracture density of the target area is determined according to the shear wave splitting parameters, the preset shear wave splitting parameters, and the linear mapping relationship between fracture density.

[0112] In a specific scenario example, the crack density determination method and apparatus based on pure shear wave reflection amplitude domain information provided in this specification can be applied to address the technical problem that existing technologies mainly rely on the propagation time difference of fast and slow shear waves for prediction, resulting in weak ability to capture crack details and low resolution. The specific implementation process may include the following:

[0113] Verification was conducted using measured data from a certain region. To address the missing slow shear wave velocity and shear wave splitting parameters in the original logging curve of the target well (hereinafter referred to as Well A) within the work area, this scheme introduced a nine-component zero-source-distance vertical seismic profile (VSP) signal from a neighboring well (hereinafter referred to as Well B).

[0114] Based on the geological background of the target strata in the work area being gently undulating and fracture development mainly controlled by the compaction of the overlying strata, the fracture density in this area was determined to have spatial consistency. Therefore, the shear wave splitting parameters obtained from well B were used as constraints for the model construction of well A. Using this test model, combined with a Ricker wavelet with a dominant frequency of 30 Hz, fast and slow shear wave pre-stack angle gathers were synthesized at two observation azimuths: the isotropic plane and the symmetry axis plane. The reflection coefficient was calculated using an exact equation, and the incident angles were selected at 10°, 20°, and 30° to provide basic data for subsequent inversion.

[0115] Considering the strong nonlinearity between the original reflection coefficient equation and the elastic parameters to be inverted, which significantly increases the complexity and instability of the inversion calculation, the reflection coefficient equation is first linearized. Specifically, the impedance, velocity, and density terms in the equation are linearized and expanded, rewritten as the difference form of physical quantities in the logarithmic domain.

[0116] Building upon this, to reduce the dimensionality of the inversion parameters, the reflection coefficient equation was further reconstructed. Utilizing the physical correlation between impedance, velocity, and density in the logarithmic domain, the equation was rewritten as a representation driven only by three core parameters: fast shear wave velocity, slow shear wave velocity, and density. This reconstruction strategy not only simplifies the computational process but also effectively improves the robustness of the inversion results by reducing coupling variables.

[0117] When constructing the pre-stack multi-parameter synchronous inversion strategy, an objective function is established within a Bayesian framework. This function seeks the optimal solution for the model parameter vector by integrating the data residual term and the prior constraint term. The model-to-data mapping operator integrates incident angle information, the statistical wavelet matrix extracted from actual data, the difference operator matrix, and the coefficient matrix containing weighting coefficients.

[0118] In the specific solution process, since the inverse matrix of the objective function is difficult to obtain directly, the Gauss-Newton method is used for iterative optimization. The model update is calculated and the initial background parameter volume is gradually corrected until the preset convergence criterion is met. During this process, the number of iterations is set according to the degree of agreement between the wellbore inversion results and the actual logging data, thereby accurately separating the first shear wave velocity (fast shear wave) and the second shear wave velocity (slow shear wave) in the target area from the complex background noise.

[0119] After obtaining high-precision fast and slow shear wave velocity fields, quantitative prediction of anisotropy parameters is performed. According to rock physics theory, when the Poisson's ratio of the rocks in the work area is between 0.2 and 0.4, the corresponding conversion factor term approaches the unit value of 1.

[0120] Based on this physical discovery, a linear mapping relationship is established between the inverted shear wave splitting parameters and the fracture density. Under this specific constraint, the shear wave splitting parameters are numerically equivalent to the fracture density. Through this equivalent transformation, this scheme directly converts seismic dynamic parameters into geological indicators reflecting the intensity of underground fracture development, achieving a quantitative characterization of fracture distribution.

[0121] By complementing data from adjacent wells, the limitations of missing logging information in single wells were overcome. Logarithmic domain linearization and parameter reconstruction techniques were used to significantly reduce nonlinear interference in the reflection coefficient equation, improving the convergence efficiency and uniqueness of the solution in the inversion process. Under the Bayesian framework, mapping operators and iterative correction strategies were introduced to achieve synchronous and high-precision extraction of fast and slow shear wave velocities. Finally, by utilizing linear equivalence relations under rock physics constraints, a direct and efficient conversion from seismic observation data to fracture density was achieved, providing a more scientific and reliable decision-making basis for sweet spot prediction of complex fractured oil and gas reservoirs.

[0122] In some embodiments, a fracture density prediction method based on pre-stack multi-parameter synchronous inversion of four-component fast and slow shear waves is provided. This method first selects measured logging data from well A in a certain area, and uses the precise Zoeppritz equation combined with a Ricker wavelet with a dominant frequency of 30Hz to synthesize fast and slow shear wave angular datasets with different signal-to-noise ratios under two specific observation azimuths. Time-domain matching correction from slow shear wave to fast shear wave is then performed. Subsequently, a joint pre-stack inversion algorithm of fast and slow shear waves is used to synchronously estimate the shear wave splitting parameters reflecting the anisotropic characteristics of the formation. Finally, based on the quantitative mapping relationship between the shear wave splitting parameters and fracture density, the inversion results are used to achieve a fine quantitative prediction of the fracture density in the target area.

[0123] In some embodiments, see Figure 4 As shown, this involves the excitation and reception process of four-component fast and slow shear waves under two special observation azimuths. In the top view shown in the figure: Line 1 is arranged parallel to the fracture strike, and Line 2 is arranged perpendicular to the fracture strike; at each excitation point (Source), shear wave components parallel to the fracture strike are excited respectively. ) and the transverse wave component perpendicular to the crack direction ( The system records the corresponding horizontal shear wave (SH) and vertical shear wave (SV) signals at each receiver point. This specific acquisition layout ensures that the observation azimuth strictly corresponds to the azimuth axis of the crack, thereby enabling the direct capture of the dynamic characteristics differences of fast and slow shear waves on different symmetry components. This provides a high-fidelity four-component seismic data foundation for subsequent synchronous inversion based on linearized reconstruction of the reflection coefficient equation.

[0124] In some embodiments, see Figure 5 As shown, firstly, three sets of seismic gather data recorded under a specific acquisition layout were acquired. The first set is isotropic fast shear wave data, whose excitation and reception components are both vertical shear waves (SV); the second set is isotropic slow shear wave data, whose excitation and reception components are both horizontal shear waves (SH); and the third set is symmetry axis fast shear wave data, whose excitation and reception components are also horizontal shear waves (SH).

[0125] The three sets of seismic gathers, encompassing different amplitude and phase characteristics, are used as input and substituted into the pre-defined pre-stack joint inversion framework. During the synchronous inversion process, the algorithm directly solves for and outputs three key physical parameters of the target region using the linearized reconstruction form of the reflection coefficient equation: fast shear wave velocity. Slow transverse wave velocity and formation density .

[0126] After obtaining the high-precision velocity field, quantitative characterization is performed using a pre-defined calculation model for shear wave splitting parameters. Specifically, based on the inverted ratio of fast to slow shear wave velocities, the shear wave splitting parameters are determined using the following mapping equation. :

[0127] This step enables the direct conversion from seismic dynamic response to anisotropic physical parameters, effectively characterizing the fracture development intensity of underground reservoirs.

[0128] While this specification provides the steps of operation for the methods described in the embodiments or flowcharts, more or fewer steps may be included based on conventional or non-inventive means. The order of steps listed in the embodiments is merely one possible order of execution among many steps and does not represent the only possible order. In actual device or client product execution, the methods shown in the embodiments or drawings may be executed sequentially or in parallel (e.g., in a parallel processor or multi-threaded processing environment, or even a distributed data processing environment). The terms "comprising," "including," or any other variations thereof are intended to cover a non-exclusive inclusion, such that a process, method, product, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, product, or apparatus. Without further limitations, the presence of other identical or equivalent elements in a process, method, product, or apparatus that includes said elements is not excluded. The terms "first," "second," etc., are used to denote names and do not indicate any particular order.

[0129] Those skilled in the art will also know that, besides implementing the controller using purely computer-readable program code, the same functions can be achieved by logically programming the method steps, making the controller function as logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers (PLCs), and embedded microcontrollers. Therefore, such a controller can be considered a hardware component, and the devices within it used to implement various functions can also be considered structures within that hardware component. Alternatively, the devices used to implement various functions can be considered as both software modules implementing the method and structures within a hardware component.

[0130] As can be seen from the above description of the embodiments, those skilled in the art can clearly understand that this specification can be implemented by means of software plus necessary general-purpose hardware platforms. Based on this understanding, the technical solutions of this specification can essentially be embodied in the form of a software product. This computer software product can be stored in a storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, mobile terminal, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments of this specification.

[0131] Although this specification has been described by way of examples, those skilled in the art will recognize that many variations and modifications are possible without departing from the spirit of this specification, and it is intended that the appended claims cover such variations and modifications without departing from the spirit of this specification.

Claims

1. A method for determining crack density based on pure transverse wave reflection amplitude domain information, characterized in that, include: Acquire first shear wave data, second shear wave data, and logging data for the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively; Data preprocessing is performed on the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information; Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information are synchronously inverted to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; The crack density of the target region is determined based on the shear wave splitting parameters, the preset mapping relationship between the shear wave splitting parameters and the crack density.

2. The method according to claim 1, characterized in that, The data preprocessing of the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and stratigraphic interpretation information includes: Pre-stack processing is performed on the first shear wave data and the second shear wave data to determine the first shear wave amplitude gather and the second shear wave amplitude gather; The well logging data is used to perform well-seismic calibration on the first and second shear wave amplitude gathers to determine seismic wavelets and layer interpretation information.

3. The method according to claim 2, characterized in that, The method involves simultaneously inverting the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters, including: Based on the logging data and the stratigraphic interpretation information, an initial background parameter volume containing information outside the seismic frequency band is determined; Using a preset pre-stack joint inversion objective function, the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume are synchronously inverted to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region. The shear wave splitting parameters are determined based on the first shear wave velocity, the second shear wave velocity, and the density.

4. The method according to claim 3, characterized in that, The step of determining the initial background parameter volume containing out-of-seismic frequency band information based on the logging data and the stratigraphic interpretation information includes: Using the P-wave velocity curve, S-wave velocity curve, and density curve in the logging data, spatial interpolation is performed within the structural framework constrained by the stratigraphic interpretation information to construct a spatially continuous initial parameter field. By performing frequency filtering on the initial parameter field, low-frequency trend components and high-frequency detail components outside the effective seismic frequency band are extracted. The low-frequency trend component and the high-frequency detail component are spectrally reconstructed to determine the initial background parameter volume.

5. The method according to claim 3, characterized in that, The method of simultaneously inverting the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the initial background parameter volume using a preset pre-stack joint inversion objective function to determine the first shear wave velocity, the second shear wave velocity, and the density of the target region includes: Using the linearized and reconstructed reflection coefficient equation and the seismic wavelet, a mapping operator is constructed; wherein, the mapping operator is used to characterize the mapping relationship between subsurface elastic parameters and seismic response; The data residual term is constructed based on the mapping operator, and the prior constraint term is constructed using the initial background parameter volume. By combining the data residual term and the prior constraint term, the preset pre-stack joint inversion objective function is obtained. Using the preset pre-stack joint inversion objective function, the first shear wave amplitude gather and the second shear wave amplitude gather are iteratively solved to determine the parameter update amount; The initial background parameter volume is corrected using the parameter update amount until the preset convergence condition is met, thereby obtaining the first and second shear wave velocities of the target region.

6. The method according to claim 3, characterized in that, The step of determining the shear wave splitting parameters based on the first shear wave velocity, the second shear wave velocity, and the density includes: The square ratio term is determined based on the ratio of the first shear wave velocity to the second shear wave velocity; Subtracting the preset constant correction term from the square ratio term yields the intermediate difference term; The intermediate difference term is scaled using a preset scaling factor to determine the shear wave splitting parameters of the target region.

7. The method according to claim 1, characterized in that, The method further includes: Based on the crack density, determine the spatial distribution characteristics of crack development intensity within the target area; Based on the spatial distribution characteristics, the exploration potential of the target area is determined.

8. A device for determining crack density based on pure transverse wave reflection amplitude domain information, characterized in that, include: The data acquisition module is used to acquire first shear wave data, second shear wave data, and logging data of the target area; wherein, the first shear wave data and the second shear wave data are fast shear wave pre-stack data and slow shear wave pre-stack data acquired based on two orthogonal observation azimuths, respectively. The data processing module is used to preprocess the first shear wave data, the second shear wave data, and the logging data to determine the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the stratigraphic interpretation information. The parameter determination module is used to simultaneously invert the first shear wave amplitude gather, the second shear wave amplitude gather, the seismic wavelet, and the layer interpretation information using a preset pre-stack joint inversion objective function to determine the shear wave splitting parameters; wherein, the preset pre-stack joint inversion objective function is constructed based on the reflection coefficient equation of linearized reconstruction; The density determination module is used to determine the crack density of the target area based on the shear wave splitting parameters and the preset linear mapping relationship between the shear wave splitting parameters and crack density.

9. An electronic device, characterized in that, It includes a processor and a memory for storing processor-executable instructions, wherein the processor, when executing the instructions, implements the steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, It stores computer instructions that, when executed by a processor, implement the steps of the method according to any one of claims 1 to 7.

Images