A method, apparatus and device for determining fracture density

Abstract

The application relates to a crack density determination method, device and equipment, and relates to the technical field of oil and gas field development, and comprises the following steps: acquiring fast shear wave post-stack seismic data and slow shear wave post-stack seismic data of a target area; performing shear wave impedance inversion on the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data to obtain fast shear wave wave impedance data and slow shear wave wave impedance data; calculating shear wave splitting parameter data according to the fast shear wave wave impedance data and the slow shear wave wave impedance data; the shear wave splitting parameter data is used for representing the azimuthal anisotropy strength of the formation medium of the target area; and determining crack density data of the target area according to the shear wave splitting parameter data.

Description

Technical Field

[0001] The embodiments in this specification relate to the field of oil and gas field development technology, specifically to a method, apparatus, and equipment for determining fracture density. Background Technology

[0002] Fracture systems play a crucial role in the accumulation and development of oil and gas, especially in unconventional reservoirs such as shale oil and gas, where fracture prediction has become a core component of reservoir sweet spot evaluation. Compared to P-waves, shear waves are more sensitive to the anisotropic response induced by fractures. When passing through fractured formations, shear waves split, generating fast and slow shear waves.

[0003] Seismic inversion is a crucial technique for extracting formation physical parameters from seismic data. Its fundamental principle is to infer the physical properties of the subsurface medium from the seismic response received at the surface, thereby enabling reservoir prediction and fluid identification. Post-stack inversion is a commonly used technique, which can be broadly categorized by implementation method into trace integral inversion, recursive inversion, and model-constrained inversion, each with its specific applicable conditions and scope. Conventional acoustic impedance inversion typically involves first inverting the reflection coefficient from post-stack seismic data, and then calculating the formation acoustic impedance using a recursive approach. However, this method suffers from several major problems: first, it is highly dependent on the initial acoustic impedance model; inaccurate assignment of the acoustic impedance of the first layer will lead to systematic deviations in the inversion results of all subsequent layers; second, there is an error accumulation effect during trace integration, meaning that the inversion error of the reflection coefficient of a certain layer will propagate and affect the acoustic impedance calculation results of all layers below it.

[0004] The accumulated deviations and errors in the aforementioned wave impedance inversion will directly affect the accuracy and reliability of subsequent crack density prediction. Summary of the Invention

[0005] The purpose of the embodiments in this specification is to provide a method, apparatus, and device for determining crack density, so as to overcome the problems of poor accuracy and reliability of crack density prediction in existing methods.

[0006] To address the aforementioned technical problems, this specification provides a method for determining crack density, comprising: Acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area; Shear wave impedance inversion was performed on the fast shear wave post-stacked seismic data and the slow shear wave post-stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data. Based on the fast shear wave impedance data and the slow shear wave impedance data, shear wave splitting parameter data are calculated; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area. Based on the shear wave splitting parameter data, the crack density data of the target region is determined.

[0007] Further, the step of performing shear wave impedance inversion on the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data to obtain fast shear wave impedance data and slow shear wave impedance data includes: Based on the well logging data and seismic interpretation stratigraphic data of the target area, initial wave impedance models for fast shear waves and slow shear waves are constructed. Using the initial wave impedance model as a priori constraint, the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data are inverted respectively to obtain the fast shear wave impedance data and the slow shear wave impedance data.

[0008] Further, the inversion of the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data using the initial wave impedance model as a priori constraint to obtain the fast shear wave impedance data and the slow shear wave impedance data includes: Construct an inversion objective function; the inversion objective function includes a likelihood term based on fitting post-stack seismic data and a regularization term based on the initial wave impedance model; Solve for the solution that minimizes the inversion objective function to obtain the fast shear wave impedance data and the slow shear wave impedance data.

[0009] Furthermore, the construction of the inversion objective function includes: Construct the following inversion objective function: ; In the formula, Let m be the inversion objective function, m be the parameter matrix to be inverted, d be the post-stack impedance seismic data, and G be the forward modeling operator. For regularization terms, This represents the weighting coefficient of the regularization term, which is proportional to the noise level. This represents the initial wave impedance model.

[0010] Further, the calculation of shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data includes: Based on the fast shear wave impedance data and the slow shear wave impedance data, the shear wave splitting parameter data are calculated using the following first mapping model: ; In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the data for shear wave splitting parameters.

[0011] Furthermore, the method also includes: A second mapping model is constructed between the fast shear wave impedance, the slow shear wave impedance, and the anisotropic parameters of the HTI medium. Construct a third mapping model between the anisotropic parameters of HTI medium and the shear wave splitting parameters; Based on the second and third mapping models, the first mapping model under the preset anisotropic conditions is constructed; the preset anisotropic conditions are used to characterize that the shear wave splitting parameter is less than a preset threshold.

[0012] Further, determining the crack density data of the target region based on the shear wave splitting parameter data includes: Based on the fourth mapping model, the shear wave splitting parameter data is used as the crack density data of the target region; the fourth mapping model is used to characterize that the shear wave splitting parameter is equivalent to the crack density under a preset anisotropic condition; the preset anisotropic condition is used to characterize that the shear wave splitting parameter is less than a preset threshold.

[0013] Furthermore, the method also includes: A fifth mapping model between shear wave splitting parameters and crack density is constructed; Based on the fifth mapping model, the fourth mapping model under the preset anisotropic conditions is determined.

[0014] Furthermore, embodiments of this specification provide a crack density determination device, comprising: The acquisition module is used to acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area. The inversion module is used to perform shear wave impedance inversion on the fast shear wave stacked seismic data and the slow shear wave stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data. The calculation module is used to calculate shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area. The determination module is used to determine the crack density data of the target region based on the shear wave splitting parameter data.

[0015] Furthermore, embodiments of this specification provide a computer device, including: Memory, used to store computer programs; A processor is used to execute the computer program to implement the above-described crack density determination method.

[0016] As can be seen from the technical solutions provided in the embodiments of this specification above, these embodiments directly utilize the high sensitivity of shear waves to fracture anisotropy by acquiring and processing fast and slow shear wave post-stack seismic data in parallel. This provides an effective data foundation for subsequent quantitative analysis and overcomes the inaccuracy problem of indirectly inferring anisotropy solely from P-wave data. Secondly, by performing independent impedance inversion on the two sets of data, seismic amplitude information is transformed into more stable impedance parameters less affected by fluid flow, while fully preserving the impedance difference information caused by fractures. Next, by directly comparing the fast and slow wave impedance data to calculate the shear wave splitting parameters, a direct and quantitative extraction of the azimuth anisotropy intensity of the formation is achieved. Based on these shear wave splitting parameters, fracture density is determined, establishing a clear and quantitative conversion path from original seismic observations to geological engineering parameters. This enables spatial quantitative prediction of fracture density, and the results can be directly used for reservoir sweet spot evaluation and development decisions, significantly improving the accuracy and reliability of the prediction results. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments or prior art of this specification, the accompanying drawings used in the description of the embodiments or prior art will be briefly introduced below.

[0018] Figure 1 This is a flowchart of a crack density determination method provided in the embodiments of this specification; Figure 2 This is a schematic diagram of the fast and slow shear wave impedance inversion process in a crack density determination method provided in the embodiments of this specification; Figure 3 This is a schematic cross-sectional view of a fast shear wave after stacking, provided in a specific embodiment of the embodiments of this specification; Figure 4 This is a schematic cross-sectional view of a slow shear wave stacked in a specific embodiment provided in this specification. Figure 5 This is a schematic diagram of well vibration calibration provided in a specific embodiment of the embodiments of this specification; Figure 6 This is a schematic diagram of an initial fast shear wave impedance model constructed using well logging data and stratigraphic information in a specific embodiment provided in this specification. Figure 7 This is a schematic diagram of an initial slow shear wave impedance model constructed using well logging data and stratigraphic information in a specific embodiment provided in this specification. Figure 8 This is a schematic diagram of the fast shear wave impedance obtained by inversion in a specific embodiment provided in this specification. Figure 9 This is a schematic diagram of the slow shear wave impedance obtained by inversion in a specific embodiment provided in this specification. Figure 10 This is a schematic diagram of the crack density determined based on fast shear wave impedance and slow shear wave impedance in a specific embodiment provided in this specification. Figure 11 This is a schematic diagram of the structural composition of a crack density determination device provided in the embodiments of this specification; Figure 12 This is a schematic diagram of the structural composition of a computer device provided in the embodiments of this specification. Detailed Implementation

[0019] 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.

[0020] It should be noted that the terms "first," "second," etc., used in this specification, claims, and the foregoing drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, apparatus, product, or device that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0021] This specification provides an embodiment of a method for determining crack density, referring to... Figure 1 As shown, the specific implementation may include the following steps: S101: Acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area.

[0022] In some embodiments, the target area can be a specific subsurface space selected for fracture system studies and sweet spot prediction during oil and gas exploration and development. It is the final application object and spatial carrier for all data processing and inversion calculations.

[0023] In some embodiments, the target area may be defined by a three-dimensional geographic coordinate range (longitude, latitude) and a vertical depth or time range (such as a two-way trip from the top to the bottom of a stratum). Its core geological components are potential reservoir segments with possible fractures, such as shale layers, tight sandstone layers, etc.

[0024] In some embodiments, fast shear wave post-stack seismic data can be basic geophysical observation data used to invert fast shear wave impedance and then calculate anisotropy parameters. It records the amplitude information of shear wave components in the seismic wave field whose polarization direction is parallel to the fracture strike (or the direction of the maximum horizontal principal stress of the stratum) after being stacked.

[0025] In some embodiments, fast shear wave post-stack seismic data can be a three-dimensional / two-dimensional matrix containing information in two dimensions: (a) spatial location information, such as common center point (CDP) gather number, line number, etc.; and (b) amplitude-time information, i.e., the amplitude value of each record at different sampling points at different times, which reflects the spatial variation of the reflection intensity of fast shear waves by the subsurface medium.

[0026] In some embodiments, slow shear wave post-stack seismic data can be paired observation data used to invert slow shear wave impedance and to jointly calculate wave velocity differences (i.e., anisotropy) with fast shear wave data. It records the post-stack amplitude response of the shear wave component with polarization perpendicular to the fracture strike.

[0027] In some embodiments, the data structure and dimensions of slow shear wave post-stack seismic data are exactly the same as those of fast shear wave post-stack seismic data, and they are strictly aligned in space and time (i.e., each CDP point and each time sampling point corresponds one-to-one), ensuring comparability between the two. Their amplitude differences directly contain information about the anisotropy caused by cracks.

[0028] In some embodiments, post-stack seismic data can be seismic data that has undergone a series of conventional seismic data processing procedures (such as static correction, dynamic correction, stacking, etc.) and can relatively intuitively reflect the information of subsurface reflection interfaces. Compared with pre-stack data, it has a higher signal-to-noise ratio and is the most commonly used input data for wave impedance inversion.

[0029] In some embodiments, post-stack seismic data can be the result of key processing steps such as azimuth rotation, dynamic correction, and horizontal stacking of raw multi-component seismic acquisition data. Azimuth rotation involves rotating the shear wave components of the original record to a coordinate system aligned with the fracture azimuth, thereby separating the fast and slow components. Horizontal stacking enhances the effective signal and suppresses random noise by superimposing data at different offsets.

[0030] In some embodiments, fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area can be acquired.

[0031] In some embodiments, data files labeled as fast shear wave post-stack data volumes and slow shear wave post-stack data volumes can be identified and read from the multi-component seismic data processing results library of the target work area. These data may be outputs from previous anisotropic processing procedures (including azimuth anisotropic velocity analysis, shear wave splitting correction, and component separation).

[0032] In some embodiments, it can be checked whether the two data volumes are completely consistent in terms of observation system (e.g., survey line number, CDP point range, time sampling length, sampling rate). Spatial alignment is verified to ensure that the seismic traces of the same spatial location point in different data volumes strictly correspond. At the same time, the signal-to-noise ratio of the data is initially evaluated to confirm that it meets the basic requirements of the inversion method for input data.

[0033] In some embodiments, depending on the input requirements of subsequent inversion, the data may need to undergo preprocessing such as format conversion, redefinition of trace head information, or amplitude calibration to ensure that the data can be correctly read and interpreted, resulting in two sets of spatially and temporally fully registered and qualified fast and slow shear wave post-stack seismic data.

[0034] Fast and slow shear wave data are themselves direct seismic responses to strata azimuth anisotropy. Obtaining this pair of data means directly acquiring key differences reflecting the existence and orientation of fractures, avoiding the uncertainties and multiple solutions that might be introduced by indirectly inferring anisotropy from P-wave data. By acquiring professionally processed and separated fast and slow wave data, it is ensured that the input data physically correspond strictly to the responses of the same subsurface point in different polarization directions. This strict pairing relationship is the fundamental guarantee for the accurate extraction of wave impedance differences (i.e., anisotropy intensity) in subsequent joint inversion.

[0035] S102: Perform shear wave impedance inversion on the fast shear wave stacked seismic data and the slow shear wave stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data.

[0036] In some embodiments, shear wave impedance inversion can be a geophysical technique specifically designed to quantitatively estimate the elastic parameter of subsurface shear wave impedance from shear wave seismic data. Its purpose is to transform observable seismic amplitude data into a profile or volume of subsurface physical properties that directly characterizes the mechanical properties of rocks (the product of density and shear wave velocity).

[0037] In some embodiments, shear wave impedance inversion includes: constructing the inversion objective function, establishing an initial model, optimizing the solution algorithm, and iterative updating and termination determination. Essentially, it solves an inverse problem in a geophysical sense.

[0038] In some embodiments, fast shear wave impedance data can be inverted data describing the impedance distribution of the subsurface medium along the fracture direction (fast wave direction). It is one of the fundamental inputs for calculating anisotropy parameters, and its spatial variation mainly reflects background information on lithology, porosity, and fluid variations.

[0039] In some embodiments, fast shear wave impedance data can be a three-dimensional data volume matching the spatial extent of the input seismic data. Each data point represents the fast shear wave impedance value (Z_S_fast) at a spatial location (X, Y, T, or Depth), calculated as the product of formation density (ρ) and fast shear wave velocity (V_S_fast): Z_S_fast = ρ V_S_fast.

[0040] In some embodiments, slow shear wave impedance data can be inverted impedance distribution data describing the vertical fracture orientation (slow wave direction) of the subsurface medium. Its relative difference from fast shear wave impedance data is the direct source for extracting shear wave splitting parameters (i.e., anisotropic intensity).

[0041] In some embodiments, the slow shear wave impedance data can be a three-dimensional data volume spatially aligned with the fast shear wave impedance data, where each data point represents the slow shear wave impedance value (Z_S_slow), calculated using the formula: Z_S_slow = ρ V_S_slow.

[0042] In some embodiments, the initial wave impedance model can provide a reasonable, geologically meaningful starting point and constraint framework for the solution process in model-based inversion. Its role is to provide low-frequency trend information missing from seismic data and guide the inversion towards convergence in a direction consistent with geological laws, avoiding non-uniqueness and instability of the solution.

[0043] In some embodiments, the initial wave impedance model can be established by combining well logging data (which provides high-resolution vertical information) with seismically interpreted stratigraphic levels (which provide a lateral structural framework). Specifically, it can be a smooth, broadband impedance data volume formed by extrapolating and interpolating the wave impedance curves obtained from well logging at well points in three-dimensional space based on the seismically interpreted stratigraphic levels and formation contact relationships.

[0044] In some embodiments, step S102 may specifically include: performing shear wave impedance inversion on the fast shear wave stacked seismic data and the slow shear wave stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data.

[0045] In some embodiments, a complete wave impedance inversion process can be run independently for fast and slow shear wave post-stack seismic data. Each inversion uses the same inversion algorithm (e.g., model-based recursive inversion, sparse pulse inversion, or iterative inversion within a Bayesian framework), but the input data differs. The forward modeling process can be based on a one-dimensional convolution model: synthetic seismic trace = seismic wavelet reflection coefficient sequence, where the reflection coefficients are calculated from the wave impedance. The inversion algorithm continuously adjusts the wave impedance model to minimize the error between the synthesized fast (or slow) shear wave seismic trace and the actually observed fast (or slow) shear wave seismic trace, ultimately outputting the corresponding wave impedance data volume.

[0046] By performing independent but parallel inversions of fast and slow wave data, it is ensured that the wave impedance difference information caused by the crack is completely preserved and encapsulated in two independent wave impedance data volumes, laying the foundation for direct calculation of the difference (anisotropy).

[0047] In some embodiments, step S102 may further include: constructing initial wave impedance models for fast shear waves and slow shear waves based on well logging data and seismic interpretation stratigraphic data of the target area; and using the initial wave impedance models as prior constraints, inverting the fast shear wave post-stack seismic data and slow shear wave post-stack seismic data respectively to obtain the fast shear wave impedance data and slow shear wave impedance data.

[0048] In some embodiments, seismic wavelets for fast and slow wave data inversion can be statistically extracted from well-side seismic trace data and well logging composite records, respectively. These wavelets can be used to construct a forward model (convolution model).

[0049] In some embodiments, an initial shear wave impedance model can be established. Specifically: using logging data from multiple wells within the target area, the fast and slow shear wave impedance curves for each well are calculated; time-depth calibration is completed by matching the synthetic records with the seismic traces near the wells; based on the stratigraphic and fault framework provided by seismic interpretation, geostatistical methods (such as Kriging interpolation) or other interpolation methods are used to extrapolate the impedance information of the well points to the entire three-dimensional space under strata-controlled constraints, generating initial impedance models for fast and slow shear waves respectively. This initial impedance model integrates high-frequency details from well logging with the tectonic morphology of the seismic event.

[0050] In some embodiments, model-based wave impedance inversion can be performed. Specifically, inversions are performed independently but in parallel on fast and slow shear wave data volumes: An inversion objective function is constructed, comprising two terms: a data residual term (minimizing the difference between the forward-modeled composite data and the actual seismic data) and a model constraint term (ensuring the inversion results do not deviate too far from the initial model, i.e., prior constraints). Iterative solution: An optimization algorithm such as the Gauss-Newton method is used for iterative solution. In each iteration, the gradient of the objective function with respect to the model parameters is calculated, and the model parameters (i.e., wave impedance) are updated until the data residual is less than a preset threshold or the maximum number of iterations is reached. This process, using the initial wave impedance model as the starting point for iteration and regularization constraints, significantly improves the efficiency and stability of the solution.

[0051] In some embodiments, after the iteration is completed, the finally converged three-dimensional data volumes of fast shear wave impedance and slow shear wave impedance are output respectively. These two data volumes correspond completely in space, laying the foundation for the next step of calculating their differences (i.e., anisotropy).

[0052] An inversion method based on an initial model was employed to supplement reliable low-frequency information, solving the inherent problem of missing low-frequency components in seismic data and giving the inverted wave impedance results absolute physical meaning. The initial model-based inversion method significantly constrains the solution space, avoiding the inversion from getting trapped in local optima or unphysical solutions, thus improving the lateral continuity and geological rationality of the inversion results. Furthermore, through separate but parallel inversion processes, strictly spatially registered fast and slow wave impedance data were obtained. This strategy of separate inversion and joint use ensures that the independent calculation errors of the two impedance data volumes are minimized. This allows for the maximal preservation and highlighting of the true anisotropic signal caused by fractures when subsequently calculating their relative differences (i.e., shear wave splitting parameters), suppressing shared background noise, and directly improving the sensitivity and accuracy of the final fracture density prediction.

[0053] S103: Calculate the shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area.

[0054] In some embodiments, shear wave splitting parameter data can be geophysical parameters used to quantitatively characterize the intensity of azimuth anisotropy in the subsurface medium. Their magnitude directly reflects the degree of drastic variation in seismic shear wave propagation velocity with azimuth caused by directional fractures, bedding, or stress fields.

[0055] In some embodiments, the shear wave splitting parameter data can be a three-dimensional data volume with a grid that perfectly matches the input wave impedance data volume. At each grid point (e.g., a CDP point, a time sampling point), the data value is a dimensionless scalar. The physical definition of this scalar can be based on velocity, i.e., the relative difference between the velocities of fast and slow shear waves.

[0056] In some embodiments, the azimuth anisotropy intensity can be the degree to which the shear wave propagation characteristics depend on the propagation direction or polarization direction. A higher intensity value indicates a more significant azimuthal difference in the physical properties of the medium.

[0057] In some embodiments, azimuth anisotropy intensity can be an abstract concept, and its specific numerical representation is the shear wave splitting parameter.

[0058] In some embodiments, the first mapping model can be a specific mathematical formula or algorithm for directly calculating shear wave splitting parameters from the impedance of fast and slow shear waves. It maps the difference in wave impedance to a measure of anisotropy intensity.

[0059] In some embodiments, the first mapping model can be derived based on rock physics theory and approximate derivation.

[0060] In some embodiments, the weak anisotropy condition can be a theoretical premise and assumption for the validity and applicability of the first mapping model. It limits the physical applicability of the first mapping model, that is, the anisotropy of the strata to be predicted cannot be too strong, thereby ensuring the accuracy of the simplified model.

[0061] In some embodiments, the weak anisotropy condition can be a qualitative or semi-quantitative technical condition. Quantitatively, it can be expressed as the value of the shear wave splitting parameter being much less than 1. A preset threshold is used to quantify this condition; for example, a threshold of 0.15 is set. When the expected value of the shear wave splitting parameter in the computational region is generally lower than this threshold, the weak anisotropy condition is considered to be met, and the use of the first mapping model is appropriate and highly accurate.

[0062] In some embodiments, step S103 may specifically include: calculating shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area.

[0063] In some embodiments, parallel computation can be performed for each sampling point in the three-dimensional data volume of the target region. The input is the fast shear wave impedance value (Z_fast) and the slow shear wave impedance value (Z_slow) of that point. Using a preset calculation formula or algorithm (which implicitly contains the physical relationship between the impedance difference and the anisotropy intensity), the output is the shear wave splitting parameter value (γ) of that point. After traversing all points, a shear wave splitting parameter data volume with the same dimension as the impedance data volume is generated. This process utilizes the impedance difference to directly substitute for and reflect the velocity difference, because the impedance (Z=Vρ) contains velocity (V) information, and under the assumption of a gradual vertical change, its relative difference approximates the relative velocity difference.

[0064] In some embodiments, step S103 may further include: calculating shear wave splitting parameter data using a first mapping model based on the fast shear wave impedance data and the slow shear wave impedance data; the first mapping model is used to characterize the mapping relationship between the shear wave splitting parameter and the fast shear wave impedance and the slow shear wave impedance under weak anisotropy conditions; the preset anisotropy conditions are used to characterize that the shear wave splitting parameter is less than a preset threshold.

[0065] In some embodiments, based on geological knowledge of the target area (e.g., fracture density is usually not extremely high) and prior knowledge, it can be determined that the area meets the weak anisotropy condition. If so, the pre-inverted fast and slow shear wave impedance data volumes (Z_S_fast, Z_S_slow) are loaded. For each data point, Z_S_fast(i,j,k) and Z_S_slow(i,j,k) are substituted into the specific formula of the first mapping model for calculation. All calculated γ(i,j,k) are arranged in their original spatial order to generate the final shear wave splitting parameter data volume.

[0066] In some embodiments, during or after the calculation, the data volume can be scanned based on a preset threshold (e.g., 0.15). If the calculation results for most regions exceed the threshold, the system can issue a warning, indicating that the "weak anisotropy" assumption may not be applicable and that the data needs to be reviewed or a more complex model should be adopted.

[0067] By explicitly stating that the model is used under weak anisotropy conditions, which is usually met for most fractured reservoirs (where the fracture density is not extremely high), the simplified model can approximate the theoretical true value with extremely high accuracy, ensuring the reliability of the extracted anisotropy intensity information.

[0068] S104: Determine the crack density data of the target region based on the shear wave splitting parameter data.

[0069] In some embodiments, fracture density data can be used to quantitatively describe the degree of fracture development within a target area, serving as a geological basis for evaluating reservoir effectiveness, selecting sweet spots, and optimizing drilling and fracturing strategies. Its numerical value directly indicates the number of fractures per unit volume or unit area.

[0070] In some embodiments, the fracture density data can be a three-dimensional spatial data volume (or two-dimensional data on a two-dimensional profile) with the same spatial grid as the input seismic and intermediate parameter data. The data value at each grid point is a dimensionless scalar that represents the areal or linear density of the fracture. The larger the value, the more developed the fractures at that location.

[0071] In some embodiments, the fourth mapping model can be a specific mathematical relationship or algorithmic rule for directly converting shear wave splitting parameters (γ) into fracture density (e). It establishes a quantitative bridge between geophysical observation responses and geological properties.

[0072] In some embodiments, the fourth mapping model can be derived based on rock physics theory, presenting a direct equivalent or linear proportional relationship. For example, it can be expressed as: e = γ / k, where k is a proportionality coefficient related to properties such as the rock background Poisson's ratio. Under weak anisotropy conditions (Poisson's ratio in the range of 0.2-0.4), the coefficient k can be approximated as 1, thus simplifying the model to e ≈ γ, that is, directly using the value of the shear wave splitting parameter as the value of the fracture density.

[0073] In some embodiments, step S104 may specifically include: determining the crack density data of the target region based on the shear wave splitting parameter data.

[0074] In some embodiments, the shear wave splitting parameter data volume (γ(x,y,t)) of the entire target area calculated in step S103 can be converted or calibrated. Specifically, based on the rock physics principle that there is a quantitative correlation between shear wave splitting parameters and fracture density, a preset conversion relationship (which can be a formula, an empirical chart, or a calibration function) is used to map the γ value at each spatial location to a corresponding fracture density e value. After traversing all data points, a fracture density data volume (e(x,y,t)) corresponding to a spatial location is generated. The establishment of this conversion relationship can rely on prior knowledge of the work area, core experimental calibration, or correction of well logging interpretation results.

[0075] In some embodiments, step S104 may further include: using the shear wave splitting parameter data as the crack density data of the target region based on a fourth mapping model; the fourth mapping model is used to characterize the shear wave splitting parameter as equivalent to crack density under a preset anisotropic condition; the preset anisotropic condition is used to characterize the shear wave splitting parameter as being less than a preset threshold.

[0076] In some embodiments, the input shear wave splitting parameter data volume (γ) can be checked to confirm whether its numerical range meets the weak anisotropy condition (i.e., whether it is generally below a preset threshold, such as 0.15). If it meets the condition, proceed to the next step; if it does not meet the condition, a more complex nonlinear mapping model can be used or a warning can be issued.

[0077] In some embodiments, after confirming the applicable conditions, a fourth mapping model is used for point-by-point calculation. In the most preferred embodiment, this model is a simple equivalence relation: e(i,j,k) = γ(i,j,k). This means that each value in the shear wave splitting parameter data volume is directly assigned to the corresponding position in the fracture density data volume. If a linear model with coefficients, e = γ / k, is used, the applicable k value for the work area (from rock physics experiments or well logging calibration) can be determined first, and then the calculation can be performed.

[0078] In some embodiments, after calculations are completed for all points, a final crack density data volume is generated. This data volume spatially reveals the sweet spots and non-sweet spots of crack development.

[0079] When using the simplified model e ≈ γ, no complex calculations are required; only data assignment or copying is needed, maximizing computational efficiency and making it suitable for real-time or near-real-time processing of massive amounts of 3D data. Furthermore, this simplified model is clearly applicable to weak anisotropy conditions, which aligns with the situation in most real fractured reservoirs, ensuring the scientific validity and reliability of the final fracture density prediction results and avoiding the introduction of unnecessary complexity and parameter uncertainties.

[0080] In some embodiments, step S102 above, using the initial wave impedance model as a priori constraint, involves inverting the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data respectively to obtain the fast shear wave impedance data and the slow shear wave impedance data. Specifically, this may include: constructing an inversion objective function; the inversion objective function includes a likelihood term based on fitting the post-stack seismic data and a regularization term based on the initial wave impedance model; and solving for the solution that minimizes the inversion objective function to obtain the fast shear wave impedance data and the slow shear wave impedance data.

[0081] In some embodiments, prior knowledge about the subsurface model (initial model) can be combined with observed seismic data in a probabilistic manner to derive the posterior probability distribution of the computational model parameters, thereby obtaining the most probable inversion solution.

[0082] In some embodiments, the inversion objective function can be a mathematical expression that transforms the problem of maximizing the posterior probability into a minimization problem that is easier to compute numerically, and it defines what the optimal solution is.

[0083] In some embodiments, the inversion objective function consists of a weighted sum of two terms: a likelihood term, which measures the fit difference between the synthesized theoretical seismic data and the actual observed post-stack seismic data during the inversion process. Its purpose is to force the inversion results to reasonably explain the observed seismic responses. The likelihood term can be represented by the L2 norm (sum of squares) of the residuals between the actual data and the forward modeling data. A regularization term can be used to introduce prior constraints and stabilize the inversion process. Its purpose is to penalize excessive deviations between the inversion results and the initial prior model, and to add smoothness constraints on model space variations to suppress severe oscillations in the solution caused by data noise or ill-posedness. The regularization term can also be in L2 norm form, comprising two parts: the first part ensures that the solution does not deviate too far from the prior model; the second part (gradient term) imposes smoothness constraints by penalizing the gradient in the model space.

[0084] In some embodiments, each iteration may include two steps: Local linearization: At the current model, a first-order Taylor expansion of the forward operator is performed to approximate the objective function as a quadratic function of the model update. Solving the linear system: By solving the resulting normal equation, the optimal model update for this iteration is calculated, thereby updating the model.

[0085] In some embodiments, the values ​​of the fast (or slow) shear wave impedance to be inverted at each sampling point in three-dimensional space can be arranged into a column vector. The corresponding fast (or slow) shear wave post-stack seismic data can also be arranged into a column vector.

[0086] In some embodiments, a forward modeling operator can be defined to describe the process of synthesizing seismic traces from a wave impedance model.

[0087] In some embodiments, a likelihood term can be constructed based on the assumption of a Gaussian distribution of data errors, and a regularization term can be constructed based on the prior Gaussian distribution assumption of model parameters. Based on the likelihood and regularization terms, a complete inversion objective function can be constructed.

[0088] In some embodiments, an iteration counter can be set to iteratively solve for the solution that minimizes the inversion objective function. The iteration termination condition is checked. If satisfied, the solution is output as the final inverted wave impedance data; otherwise, the iteration continues.

[0089] By organically combining seismic observation data (likelihood term) with prior geological logging knowledge (regularization term) in a probabilistic sense, the prior model in the regularization term explicitly introduces low-frequency trends and geological frameworks, while smoothing constraints suppress high-frequency noise. This synergistic constraint fundamentally suppresses the multiple solutions problem commonly encountered in inversion, making the solution more stable and more consistent with geological laws.

[0090] In some embodiments, the above-mentioned construction of the inversion objective function may specifically include: constructing the following inversion objective function: In the formula, For the inversion objective function, The parameter matrix to be inverted, For post-stack impedance seismic data, Forward operands; For regularization terms, This represents the weighting coefficient of the regularization term, which is proportional to the noise level. This represents the initial wave impedance model.

[0091] In some embodiments, the parameter matrix to be inverted It can be the direct target object of the inversion solution, that is, the numerical set of fast (or slow) shear wave impedance at all discrete sampling points in the three-dimensional space of the target region.

[0092] In some embodiments, the parameter matrix to be inverted can be a column vector or a matrix. For example, if the target region has N time sampling points and M spatial channels, then... It is a column vector of dimension (NM, 1), where each element represents the wave impedance value at a specific spatial location and time depth.

[0093] In some embodiments, post-stack impedance seismic data It can be the observation data that needs to be fitted in the inversion, that is, the input fast (or slow) shear wave post-stack seismic amplitude data.

[0094] In some embodiments, the post-stack impedance seismic data can be the same as the parameter matrix to be inverted. This is a column vector corresponding to the spatial and temporal sampling. Its data represents the seismic amplitude values ​​at a specific time sampling point. It is related to the parameter matrix to be inverted. Through forward modeling operators Related to each other.

[0095] In some embodiments, the forward operator It can be a mathematical function or operation used to describe how to extract the wave impedance model (the parameter matrix to be inverted) from the underground. It synthesizes the theoretically observable seismic data. It serves as a bridge connecting model space and data space.

[0096] In some embodiments, the forward operator It can be represented as: In the formula, W is the wavelet matrix, D is the difference operator matrix, and F is the coefficient matrix; their matrix forms are as follows: , ,as well as .

[0097] In some embodiments, in post-stack inversion based on a one-dimensional convolution model, the forward operator can be a nonlinear operator.

[0098] In some embodiments, the regularization term weight coefficient Used to control the relative weights or balance between the likelihood term and the regularization term in the objective function.

[0099] In some embodiments, the regularization term weight coefficients can be set or adaptively adjusted based on the signal-to-noise ratio (SNR) of the data. The higher the noise level, the more... A larger value gives the prior model stronger constraints, suppressing inversion instability caused by noise; conversely, for high signal-to-noise ratio data, The value can be set to a smaller value, allowing the data itself to play a greater role.

[0100] In some embodiments, a data fitting term (likelihood term) is defined: This is the observation data vector (post-stack impedance seismic data). ) and forward-modeled composite data vector The sum of the squares of the L2 norm of the differences (i.e., the sum of the squares of the differences between corresponding elements). It directly quantifies the current model (the parameter matrix to be inverted). This refers to the interpretability of actual seismic data. Minimizing this means pursuing the best match between the inversion results and seismic observations.

[0101] In some embodiments, model constraint terms (prior terms / regularization terms) are defined: This term is the squared L2 norm of the difference between the current model and the prior initial model, multiplied by the weight coefficient. This is the most reasonable initial conjecture based on well logging and seismic interpretation. Minimizing this means limiting the inversion results from deviating excessively from the geologically reasonable initial understanding.

[0102] In some embodiments, the weighted average is used to form the objective function: This function clarifies the dual criteria a good wave impedance model must meet: it must fit the data and be close to the prior. The weighting coefficients act as a mediator, their magnitude determining which criterion to favor when the two criteria conflict.

[0103] This method clearly transforms the geological inversion problem into an optimization problem with a well-defined mathematical expression, enabling subsequent automatic solutions using mature numerical optimization algorithms. This achieves objectivity and automation of the inversion process, avoiding the subjective arbitrariness of human interpretation. The introduction of a regularization term mathematically adds an anchor and buffer to the inversion solution. Even with poor data quality or an ill-posed problem, this mechanism effectively prevents meaningless oscillations or divergences in the solution, significantly improving the robustness and practicality of the inversion algorithm and overcoming the shortcomings of conventional methods, such as sensitivity to initial values ​​and error accumulation.

[0104] In some embodiments, the above-described solution that minimizes the inversion objective function to obtain the fast shear wave impedance data and slow shear wave impedance data may specifically include: calculating the partial derivative of the inversion objective function with respect to the parameter matrix to be inverted, thereby obtaining the update amount of the parameter matrix to be inverted. In the formula, The update amount of the parameter matrix to be inverted. The identity matrix is ​​used; the solution that minimizes the inversion objective function is obtained through iterative solving to obtain the fast shear wave impedance data and the slow shear wave impedance data. In the formula, For the number of iterations, For the first Shear wave impedance data from the next iteration For the first Shear wave impedance data from the next iteration.

[0105] In some embodiments, in the current iteration model At this point, the nonlinear forward modeling operator will be used. Perform a first-order Taylor expansion and update the model size. Find the partial derivative and set it to zero to find the minimum point.

[0106] In some embodiments, the update quantity is obtained by taking partial derivatives. A system of linear equations.

[0107] In some embodiments, initialization: setting =0, = . No. Next iteration: Based on the current model Calculate the residuals and Jacobian matrix. Then solve the linear system to obtain the update variables. Model update: Apply the update formula This yields a new generation model. Convergence criterion: Check if the update amount of the parameter matrix to be inverted or the change in the Jacobian matrix is ​​less than a preset threshold, or if the maximum number of iterations has been reached. If convergence has not occurred, continue iterating; if convergence has occurred, output the final result. As the wave impedance data obtained from the inversion.

[0108] The model is updated by iteratively solving the update equation directly derived mathematically from the objective function. The update direction in each iteration is one of the directions that causes the objective function to decrease the fastest at the current point (which is the optimal direction under the local quadratic approximation), ensuring a fast convergence speed for the algorithm.

[0109] In some embodiments, the calculation of shear wave splitting parameter data using a first mapping model based on the fast shear wave impedance data and the slow shear wave impedance data in step S103 above may specifically include: calculating the shear wave splitting parameter data using the following first mapping model based on the fast shear wave impedance data and the slow shear wave impedance data: In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the data for shear wave splitting parameters.

[0110] In some embodiments, the first mapping model can provide an explicit mathematical relationship that uniquely maps parallel pairs of fast and slow shear wave impedance data to transverse wave splitting parameters. Its purpose is to eliminate ambiguity and enable direct, deterministic calculations from elastic parameters to anisotropic parameters.

[0111] In some embodiments, two spatially fully registered three-dimensional data volumes generated by the inversion in step S102 can be obtained: a fast shear wave impedance data volume and a slow shear wave impedance data volume. It is ensured that for each three-dimensional mesh node (i, j, k) within the target region, there is a definite pair of input values.

[0112] In some embodiments, for each grid node (i, j, k), the corresponding fast shear wave impedance and slow shear wave impedance are substituted into the first mapping model formula to obtain the shear wave splitting parameter corresponding to that grid node. The calculation is performed on all grid nodes (i, j, k), and all calculated shear wave splitting parameters are stored in their original spatial order, thus generating the final three-dimensional data volume of the shear wave splitting parameters.

[0113] Compared to methods that require complex nonlinear inversion to obtain anisotropy parameters, the formula... It is an analytical solution. For each data point, only a single, finite-step arithmetic operation is required to obtain the result, completely avoiding iterative optimization. This makes the computation extremely fast and resource consumption very low, making it particularly suitable for real-time post-processing of massive amounts of 3D seismic inversion results, greatly improving the overall engineering practicality and processing efficiency of the method.

[0114] In some embodiments, the method further includes: constructing a second mapping model of fast shear wave impedance, slow shear wave impedance and HTI medium anisotropy parameters; constructing a third mapping model of HTI medium anisotropy parameters and shear wave splitting parameters; and constructing a first mapping model under preset anisotropy conditions based on the second mapping model and the third mapping model; wherein the preset anisotropy conditions are used to characterize that the shear wave splitting parameters are less than a preset threshold.

[0115] In some embodiments, the second mapping model can be used to characterize the quantitative relationship between the fast and slow shear wave impedances—a pair of parameters directly obtained from seismic data inversion—and the anisotropic parameters of the HTI medium. Its purpose is to connect the observed / inverted elastic parameters with the classical theoretical parameter system describing the inherent anisotropic properties of the medium.

[0116] In some embodiments, the second mapping model can be derived based on elastic wave theory.

[0117] In some embodiments, the HTI medium anisotropy parameter can be a theoretical parameter used to quantitatively describe the elastic properties of a transversely isotropic medium with a horizontal axis of symmetry, serving as a bridge connecting the seismic wave equation and the physical properties of the medium.

[0118] In some embodiments, the HTI medium anisotropy parameter can be a parameter defined for the HTI medium, which directly reflects the difference between the squares of the fast and slow shear wave velocities.

[0119] In some embodiments, the third mapping model can be used to characterize the conversion relationship between the anisotropic parameters of the HTI medium and the shear wave splitting parameters. Its purpose is to bridge the gap between classical theoretical parameters and parameters based on observable seismic splitting phenomena.

[0120] In some embodiments, the third mapping model can be an approximate transformation formula. For weakly anisotropic media, there is an approximate relationship between the shear wave splitting parameters and the anisotropic parameters of the HTI medium: The negative sign originates from the conventional differences in the direction of parameter definition. Its core is the establishment of an equivalent or proportional relationship between the anisotropic parameters of the HTI medium and the shear wave splitting parameters that can be directly estimated from seismic splitting phenomena.

[0121] In some embodiments, based on the constitutive equations of HTI media and plane wave propagation theory, the relationship between fast and slow shear wave velocities and elastic stiffness parameters can be derived. Furthermore, under the assumption of relatively gradual density changes, the velocity relationship is transformed into a wave impedance relationship, resulting in a second mapping model.

[0122] In some embodiments, the HTI medium anisotropy parameter is a static definition based on the medium's elastic stiffness, while the shear wave splitting parameter is a dynamic definition based on the observable difference in fast and slow wave velocities or time differences. Both are physically identical, but their numerical values ​​and definitions differ. Under the premise of weak anisotropy, an approximate relationship between the two can be established by analyzing the reflection coefficient or propagation response of a vertically incident shear wave.

[0123] In some embodiments, the third mapping model is substituted into the second mapping model, and the weak anisotropy condition is applied to obtain the first mapping model.

[0124] In some embodiments, the above-mentioned construction of the second mapping model between the fast shear wave impedance, the slow shear wave impedance, and the anisotropy parameters of the HTI medium may specifically include: constructing the following second mapping model between the fast shear wave impedance, the slow shear wave impedance, and the anisotropy parameters of the HTI medium: In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the anisotropic parameter data for HTI media.

[0125] In some embodiments, the above-mentioned construction of the third mapping model between the anisotropic parameters and shear wave splitting parameters of the HTI medium may specifically include: constructing the following third mapping model between the anisotropic parameters and shear wave splitting parameters of the HTI medium: In the formula, For shear wave splitting parameter data, This is the anisotropic parameter data for HTI media.

[0126] In some embodiments, parameters It is based on the definition of the constitutive equation of the medium, and the parameters This definition is based on seismic wave propagation observations (such as time difference or amplitude difference). By analyzing the reflection coefficient of shear waves at the HTI medium interface under perpendicular or small-angle incidence, or their propagation response in layered media, the functional relationship between the two can be derived. (Formula) This is precisely the manifestation of this relationship. The negative sign indicates that the two changes in opposite directions (this depends on the specific parameter definition conventions).

[0127] In some embodiments, constructing the first mapping model under preset anisotropic conditions based on the second mapping model and the third mapping model may specifically include: determining the third mapping model under the following weak anisotropic conditions: In the formula, For shear wave splitting parameter data, The data pertains to the anisotropy parameters of the HTI medium; wherein the weak anisotropy condition is used to characterize: Based on the second mapping model and the third mapping model under weak anisotropy, the first mapping model under weak anisotropy is constructed as follows: In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the data for shear wave splitting parameters.

[0128] In some embodiments, weak anisotropy conditions can be applied to simplify the third mapping model: due to the presupposition of... , and by It can be known It is also much smaller than 1, thus yielding a highly simplified linear approximation: .

[0129] In some embodiments, a simplified model and a second mapping model can be combined: the simplified relationship Substitute into the second mapping model .Will Replace with Substituting, we get: This is the final simplified formula for connecting the fast and slow wave impedances and the transverse wave splitting parameter under weak anisotropy conditions. This is precisely the formula used for direct calculation in the previous embodiment.

[0130] By introducing the condition of weak anisotropy, the complex nonlinear relationship chain (second and third mapping models) is ultimately simplified into an extremely concise formula. The formula is simple in form, fast in calculation, intuitive in physical meaning, and highly efficient and practical. The derivation clearly shows that the high accuracy of the final simplified formula depends on the premise of weak anisotropy. This is not a technical flaw, but rather a scientific definition of the method's applicability. It guides users to apply the formula in geological areas that meet this condition (most fractured reservoirs), thereby ensuring reliable engineering accuracy in fracture density prediction results.

[0131] In some embodiments, the method further includes: constructing a fifth mapping model of shear wave splitting parameters and crack density; and determining a fourth mapping model under preset anisotropic conditions based on the fifth mapping model.

[0132] In some embodiments, the fifth mapping model can be used to characterize the fundamental quantitative relationship between shear wave splitting parameters and fracture density based on rock physics principles. Its purpose is to provide the original theoretical basis for the conversion from seismic properties to geological parameters.

[0133] In some embodiments, the fifth mapping model can be derived based on a specific fractured rock physics model (such as the Hudson model or the Schoenberg linear slip model) or equivalent medium theory. This model reflects the complex dependence of anisotropic intensity on the degree of fracture development and the background properties of the rock.

[0134] In some embodiments, the fourth mapping model can be a simplified form of the fifth mapping model under specific idealized conditions. Its purpose is to provide an extremely simple direct conversion rule that requires no additional rock physics parameters, thereby greatly improving the efficiency and practicality of the final prediction step.

[0135] In some embodiments, a fractured medium rock physics model suitable for the target reservoir type can be selected or established. For example, the Hudson model can be used for a set of randomly distributed but predominantly oriented dry fractures; for more general cases, the Schoenberg linear slip model can be used.

[0136] In some embodiments, based on the selected model, the relationship between the equivalent elastic stiffness tensor and parameters such as crack density, crack aspect ratio, and background rock matrix elastic modulus (e.g., shear modulus, Poisson's ratio) is derived.

[0137] In some embodiments, parameters characterizing the anisotropy of shear waves can be extracted or calculated from the obtained equivalent elastic stiffness tensor.

[0138] In some embodiments, through the above derivation, an equation can be finally established that expresses the shear wave splitting parameters as a function of fracture density and other rock background and fracture morphology parameters.

[0139] In some embodiments, the fifth mapping model is analyzed. The weak anisotropy condition implies a low crack density, allowing for linear simplification and yielding a simplified relationship. After simplification, the complex fifth mapping model can be reduced to the fourth mapping model.

[0140] In some embodiments, the above-mentioned construction of the fifth mapping model between shear wave splitting parameters and crack density may specifically include: constructing the following fifth mapping model between shear wave splitting parameters and crack density: In the formula, For shear wave splitting parameter data, For the Poisson's ratio data of the formation medium in the target area, This refers to the crack density data for the target area.

[0141] In some embodiments, an equivalent medium model suitable for media containing randomly oriented vertical cracks can be selected, such as Hudson's thin coin-shaped crack model. In this model, the cracks are treated as low-aspect-ratio ellipsoidal inclusions. The elastic flexibility or stiffness increment of the equivalent medium under the first-order approximation is derived.

[0142] In some embodiments, parameters characterizing the anisotropy of the shear wave are calculated from the obtained equivalent elastic tensor. The elastic modulus of the background medium can be expressed as Lamé constant or Young's modulus and Poisson's ratio. Through algebraic operations, the proportionality constant is expressed as a function of Poisson's ratio, yielding... This formula is the concretized fifth mapping model, which clearly shows that when the fracture density is not high (first-order approximation is effective), the shear wave splitting parameter is proportional to the fracture density, and the proportionality coefficient is determined only by the Poisson's ratio of the rock matrix.

[0143] This model abandons purely empirical calibration, providing a precise conversion formula with clear theoretical support that can be achieved with only one key background parameter. This allows crack density prediction to move from qualitative analogy or local empirical calibration to parametric quantitative calculation, significantly improving the method's universality and scientific rigor. (Formula) It is clearly pointed out that, apart from the crack density itself, the Poisson's ratio of the background rock is the most important factor affecting the shear wave splitting response.

[0144] In some embodiments, determining the fourth mapping model under preset anisotropic conditions based on the fifth mapping model may specifically include: determining the fourth mapping model under the following weak anisotropic conditions based on the fifth mapping model: In the formula, For shear wave splitting parameter data, The crack density data for the target region; wherein, the weak anisotropy condition is used to characterize: Under the aforementioned anisotropic conditions: In the formula, This provides Poisson's ratio data for the formation medium in the target area.

[0145] In some embodiments, analyzable functions Its characteristics: Calculate its properties under the condition that Poisson's ratio is weakly anisotropic (at which point...) The range of values ​​is from 0.2 to 0.4. It can be observed that... Based on this, the fifth mapping model Simplified to the fourth mapping model This means that, under weak anisotropy conditions, the numerical values ​​of the shear wave splitting parameters can be directly used as numerical estimates of the crack density.

[0146] Received This relationship simplifies the final conversion process to a direct assignment or scaling operation that requires no additional parameters or complex calculations. This significantly lowers the application threshold and computational costs, enabling the instantaneous generation of crack density data volumes for large-scale 3D work areas and improving production efficiency.

[0147] In some embodiments, after acquiring fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area, shear wave impedance inversion can be performed on the fast shear wave post-stack seismic data and slow shear wave post-stack seismic data. The approximate reflection coefficients of the fast and slow shear waves during the inversion process can be expressed as: ; ; In some embodiments, given that the relationship between the above reflection coefficient equation and the elastic parameters to be inverted is nonlinear, which increases the complexity of the inversion calculation and enhances the instability of the inversion, the equation is linearized. The parametric impedance term in the equation can be linearly expanded as follows: ; ; The reflection coefficient equation can be further reconstructed into the following form: ; ; During the inversion process, fast shear wave impedance data are obtained. With slow shear wave impedance data ;because Therefore, it can be directly from and Further, the crack density data of the target area can be directly determined.

[0148] As can be seen from the technical solutions provided in the embodiments of this specification above, these embodiments directly utilize the high sensitivity of shear waves to fracture anisotropy by acquiring and processing fast and slow shear wave post-stack seismic data in parallel. This provides an effective data foundation for subsequent quantitative analysis and overcomes the inaccuracy problem of indirectly inferring anisotropy solely from P-wave data. Secondly, by performing independent impedance inversion on the two sets of data, seismic amplitude information is transformed into more stable impedance parameters less affected by fluid flow, while fully preserving the impedance difference information caused by fractures. Next, by directly comparing the fast and slow wave impedance data to calculate the shear wave splitting parameters, a direct and quantitative extraction of the azimuth anisotropy intensity of the formation is achieved. Based on these shear wave splitting parameters, fracture density is determined, establishing a clear and quantitative conversion path from original seismic observations to geological engineering parameters. This enables spatial quantitative prediction of fracture density, and the results can be directly used for reservoir sweet spot evaluation and development decisions, significantly improving the accuracy and reliability of the prediction results.

[0149] Reference Figure 2 The embodiments shown in this specification illustrate the impedance inversion process for fast and slow shear waves. In one specific embodiment, the obtained post-stack profiles of fast and slow shear waves are as follows: Figure 3 and Figure 4 As shown in the diagram, the vertical lines represent the locations of the wells. Figure 5 As shown, wellbore calibration is performed. Figure 6 and Figure 7 As shown, this is an initial model constructed using well logging data and stratigraphic information. The vertical lines in the figure represent the well locations. Figure 8 and Figure 9 As shown, the fast shear wave impedance and slow shear wave impedance are obtained through inversion. The vertical lines in the figure represent the well locations. Figure 10 As shown, the fracture density is determined based on the fast shear wave impedance and the slow shear wave impedance. The vertical lines in the figure represent the locations of the wells.

[0150] Based on the above-described method for determining crack density, this specification also provides embodiments of a crack density determination device. For example... Figure 11 As shown, the crack density determination device 1100 may specifically include the following modules: The acquisition module 1101 is used to acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data of the target area; Inversion module 1102 is used to perform shear wave impedance inversion on the fast shear wave stacked seismic data and the slow shear wave stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data. Calculation module 1103 is used to calculate shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area. The determination module 1104 is used to determine the crack density data of the target area based on the shear wave splitting parameter data.

[0151] In some embodiments, the inversion module 1102 described above can be specifically used for: Based on the well logging data and seismic interpretation stratigraphic data of the target area, initial wave impedance models for fast shear waves and slow shear waves are constructed. Using the initial wave impedance model as a priori constraint, the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data are inverted respectively to obtain the fast shear wave impedance data and the slow shear wave impedance data.

[0152] In some embodiments, the inversion module 1102 described above can also be used for: Construct an inversion objective function; the inversion objective function includes a likelihood term based on fitting post-stack seismic data and a regularization term based on the initial wave impedance model; Solve for the solution that minimizes the inversion objective function to obtain the fast shear wave impedance data and the slow shear wave impedance data.

[0153] In some embodiments, the inversion module 1102 described above can also be used for: Construct the following inversion objective function: ; In the formula, Let m be the inversion objective function, m be the parameter matrix to be inverted, d be the post-stack impedance seismic data, and G be the forward modeling operator. For regularization terms, This represents the weighting coefficient of the regularization term, which is proportional to the noise level. This represents the initial wave impedance model.

[0154] In some embodiments, the above-described calculation module 1103 can be specifically used for: Based on the fast shear wave impedance data and the slow shear wave impedance data, the shear wave splitting parameter data are calculated using the following first mapping model: ; In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the data for shear wave splitting parameters.

[0155] In some embodiments, the above-described calculation module 1103 can also be used for: A second mapping model is constructed between the fast shear wave impedance, the slow shear wave impedance, and the anisotropic parameters of the HTI medium. Construct a third mapping model between the anisotropic parameters of HTI medium and the shear wave splitting parameters; Based on the second and third mapping models, the first mapping model under the preset anisotropic conditions is constructed; the preset anisotropic conditions are used to characterize that the shear wave splitting parameter is less than a preset threshold.

[0156] In some embodiments, the determining module 1104 described above can be specifically used for: Based on the fourth mapping model, the shear wave splitting parameter data is used as the crack density data of the target region; the fourth mapping model is used to characterize that the shear wave splitting parameter is equivalent to the crack density under a preset anisotropic condition; the preset anisotropic condition is used to characterize that the shear wave splitting parameter is less than a preset threshold.

[0157] In some embodiments, the determining module 1104 described above can also be used for: A fifth mapping model between shear wave splitting parameters and crack density is constructed; Based on the fifth mapping model, the fourth mapping model under the preset anisotropic conditions is determined.

[0158] As can be seen from the technical solutions provided in the embodiments of this specification above, these embodiments directly utilize the high sensitivity of shear waves to fracture anisotropy by acquiring and processing fast and slow shear wave post-stack seismic data in parallel. This provides an effective data foundation for subsequent quantitative analysis and overcomes the inaccuracy problem of indirectly inferring anisotropy solely from P-wave data. Secondly, by performing independent impedance inversion on the two sets of data, seismic amplitude information is transformed into more stable impedance parameters less affected by fluid flow, while fully preserving the impedance difference information caused by fractures. Next, by directly comparing the fast and slow wave impedance data to calculate the shear wave splitting parameters, a direct and quantitative extraction of the azimuth anisotropy intensity of the formation is achieved. Based on these shear wave splitting parameters, fracture density is determined, establishing a clear and quantitative conversion path from original seismic observations to geological engineering parameters. This enables spatial quantitative prediction of fracture density, and the results can be directly used for reservoir sweet spot evaluation and development decisions, significantly improving the accuracy and reliability of the prediction results.

[0159] This specification also provides a computer device for determining fracture density, including a processor and a memory for storing processor-executable instructions. Specifically, the processor can perform the following tasks according to the instructions: acquiring fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for a target area; performing shear wave impedance inversion on the fast and slow shear wave post-stack seismic data to obtain fast and slow shear wave impedance data; calculating shear wave splitting parameter data based on the fast and slow shear wave impedance data; the shear wave splitting parameter data being used to characterize the azimuth anisotropy intensity of the formation medium in the target area; and determining the fracture density data of the target area based on the shear wave splitting parameter data.

[0160] To execute the above instructions more accurately, please refer to... Figure 12 As shown in the embodiments of this specification, another specific computer device 1200 is also provided, wherein the computer device 1200 includes a network communication port 1201, a processor 1202 and a memory 1203, and the above structures are connected by internal cables so that the various structures can perform specific data interaction.

[0161] The processor 1202 is specifically configured to: acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data of the target area; perform shear wave impedance inversion on the fast shear wave post-stack seismic data and slow shear wave post-stack seismic data to obtain fast shear wave impedance data and slow shear wave impedance data; calculate shear wave splitting parameter data based on the fast shear wave impedance data and slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area; and determine the fracture density data of the target area based on the shear wave splitting parameter data.

[0162] The memory 1203 can be used to store the corresponding instruction program.

[0163] In this embodiment, the network communication port 1201 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.

[0164] In this embodiment, the processor 1202 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.

[0165] In this embodiment, the memory 1203 includes volatile memory and non-volatile memory. The memory 1203 can include multiple layers. In digital systems, anything that can store binary data can be a memory; in integrated circuits, 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.

[0166] Furthermore, embodiments of this specification also provide a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described... Figure 1 Instructions for determining crack density as shown.

[0167] It should be understood that in the various embodiments of this specification, the sequence number of each process does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this specification.

[0168] It should also be understood that, in the embodiments of this specification, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Additionally, the character " / " in this specification generally indicates that the preceding and following related objects have an "or" relationship.

[0169] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0170] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0171] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0172] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational tasks to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The task is a function specified in one or more boxes.

[0173] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for determining crack density, characterized in that, include: Acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area; Shear wave impedance inversion was performed on the fast shear wave post-stacked seismic data and the slow shear wave post-stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data. Based on the fast shear wave impedance data and the slow shear wave impedance data, shear wave splitting parameter data are calculated; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area. Based on the shear wave splitting parameter data, the crack density data of the target region is determined.

2. The method according to claim 1, characterized in that, The step of performing shear wave impedance inversion on the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data to obtain fast shear wave impedance data and slow shear wave impedance data includes: Based on the well logging data and seismic interpretation stratigraphic data of the target area, initial wave impedance models for fast shear waves and slow shear waves are constructed. Using the initial wave impedance model as a priori constraint, the fast shear wave post-stack seismic data and the slow shear wave post-stack seismic data are inverted respectively to obtain the fast shear wave impedance data and the slow shear wave impedance data.

3. The method according to claim 2, characterized in that, The process involves inverting the fast shear wave impedance data and the slow shear wave impedance data using the initial impedance model as a priori constraint, respectively, to obtain the fast shear wave impedance data and the slow shear wave impedance data, including: Construct an inversion objective function; the inversion objective function includes a likelihood term based on fitting post-stack seismic data and a regularization term based on the initial wave impedance model; Solve for the solution that minimizes the inversion objective function to obtain the fast shear wave impedance data and the slow shear wave impedance data.

4. The method according to claim 3, characterized in that, Construct the inversion objective function, including: Construct the following inversion objective function: ； In the formula, Let m be the inversion objective function, m be the parameter matrix to be inverted, d be the post-stack impedance seismic data, and G be the forward modeling operator. For regularization terms, This represents the weighting coefficient of the regularization term, which is proportional to the noise level. This represents the initial wave impedance model.

5. The method according to claim 1, characterized in that, The calculation of shear wave splitting parameters based on the fast shear wave impedance data and the slow shear wave impedance data includes: Based on the fast shear wave impedance data and the slow shear wave impedance data, the shear wave splitting parameter data are calculated using the following first mapping model: ； In the formula, For fast shear wave impedance data, This is the impedance data for slow shear waves. This is the data for shear wave splitting parameters.

6. The method according to claim 5, characterized in that, The method further includes: A second mapping model is constructed between the fast shear wave impedance, the slow shear wave impedance, and the anisotropic parameters of the HTI medium. Construct a third mapping model between the anisotropic parameters of HTI medium and the shear wave splitting parameters; Based on the second and third mapping models, the first mapping model under the preset anisotropic conditions is constructed; the preset anisotropic conditions are used to characterize that the shear wave splitting parameter is less than a preset threshold.

7. The method according to claim 1, characterized in that, The step of determining the crack density data of the target region based on the shear wave splitting parameter data includes: Based on the fourth mapping model, the shear wave splitting parameter data is used as the crack density data of the target region; the fourth mapping model is used to characterize that the shear wave splitting parameter is equivalent to the crack density under a preset anisotropic condition; the preset anisotropic condition is used to characterize that the shear wave splitting parameter is less than a preset threshold.

8. The method according to claim 7, characterized in that, The method further includes: A fifth mapping model between shear wave splitting parameters and crack density is constructed; Based on the fifth mapping model, the fourth mapping model under the preset anisotropic conditions is determined.

9. A device for determining crack density, characterized in that, include: The acquisition module is used to acquire fast shear wave post-stack seismic data and slow shear wave post-stack seismic data for the target area. The inversion module is used to perform shear wave impedance inversion on the fast shear wave stacked seismic data and the slow shear wave stacked seismic data to obtain fast shear wave impedance data and slow shear wave impedance data. The calculation module is used to calculate shear wave splitting parameter data based on the fast shear wave impedance data and the slow shear wave impedance data; the shear wave splitting parameter data is used to characterize the azimuth anisotropy intensity of the formation medium in the target area. The determination module is used to determine the crack density data of the target region based on the shear wave splitting parameter data.

10. A computer device, characterized in that, include: Memory, used to store computer programs; A processor for executing the computer program to implement the method of any one of claims 1-8.

Images