Method and apparatus for crack prediction based on seismic data

Abstract

This application provides a method and apparatus for fracture prediction based on seismic data. The method includes: dividing seismic data into multiple sector regions according to the shot point and detector azimuth angle; determining the azimuth domain common imaging gathers of converted waves in each sector region; performing fast and slow shear wave separation on the azimuth domain common imaging gathers; determining the fast shear wave gathers, slow shear wave gathers, and the angle between the fracture orientation and the main survey line direction of the detector in each formation fracture; determining the fracture development direction of the target formation based on the fast and slow shear wave gathers; converting the azimuth domain common imaging gathers corresponding to the fast and slow shear waves into snail gathers according to the azimuth and offset of the sector region, and correcting the time difference of each channel in the snail gathers; determining the fracture density based on the time difference, and determining the fracture development index based on the density.

Description

Technical Field

[0001] This application relates to the field of geophysical technology, and in particular to a method and apparatus for predicting cracks based on seismic data. Background Technology

[0002] The presence of natural fracture systems in subsurface media plays a crucial role in controlling the formation, storage capacity, and fluid migration patterns of oil and gas reservoirs. Therefore, accurately characterizing key properties such as fracture orientation and density is essential.

[0003] Shear wave splitting is a core phenomenon and analytical method for characterizing the anisotropy of subsurface media, especially fracture characteristics. Shear wave splitting refers to the phenomenon where, when a shear wave passes obliquely across a fracture surface, it splits into a fast shear wave and a slow shear wave, with polarization directions parallel and perpendicular to the fracture plane, respectively. The geophone receives a mixture of partial projections of the fast and slow shear waves. If these two waves are not separated, the accuracy of the analysis of azimuth anisotropy variations will be affected.

[0004] PS waves simultaneously contain P-wave azimuth anisotropy and shear wave splitting information, providing richer seismic responses for fracture identification in complex strata. However, due to limitations in the coverage inhomogeneity of actual observation systems and data signal-to-noise ratio, predicting fracture parameters using pre-stack PS wave data still faces challenges. Summary of the Invention

[0005] In view of this, the purpose of this application is to propose a crack prediction method and apparatus based on seismic data that overcomes or at least partially solves the above problems.

[0006] To achieve the above objectives, a first aspect of this application provides a crack prediction method based on seismic data, characterized in that it includes:

[0007] The seismic data is divided into multiple sector regions according to the azimuth of the shot point and the geophone, and the azimuth domain common imaging gather of the converted wave in each sector region is determined.

[0008] Fast and slow shear wave separation is performed on the azimuth domain common imaging gather to determine the fast shear wave gather, the slow shear wave gather, and the angle between the fracture orientation and the main survey line direction of the geophone in each formation fracture.

[0009] Based on the fast shear wave gather and the slow shear wave gather, the development direction of fractures in the target formation is determined;

[0010] The azimuth domain co-imaging gathers corresponding to the fast and slow shear waves are converted into snail gathers according to the azimuth and offset of the sector region, and the time difference of each channel in the snail gathers is corrected.

[0011] Based on the time difference, the density of the cracks is determined, and the development index of the cracks is determined based on the density.

[0012] Optionally, fast and slow shear wave separation is performed on the azimuth domain common imaging gather to determine the fast shear wave gather and the slow shear wave gather, including:

[0013] Using the main measurement line direction and the connecting line direction of the detector as the X-axis and Y-axis of the propagation coordinate system, the propagation vector composed of the fast shear wave and the slow shear wave is determined, so that the propagation vector is transformed into the coordinate system composed of the polarization direction of the fast shear wave and the polarization direction of the slow shear wave.

[0014] Determine the reference coordinate system in which the detector is located, and rotate the propagation vector composed of the fast shear wave and the slow shear wave into the reference coordinate system to determine the seismic records after multiple shear wave separations of different strata fractures in the reference coordinate system.

[0015] The earthquake record is represented as follows:

[0016]

[0017] Among them, PSV(t) i ) and PSH(t i ) represent the seismic records of multiple shear waves received by the fracture in the i-th stratum in the X and Y components of the reference coordinate system after separation, t represents the propagation time, T represents the transpose of the rotation matrix, and Λ i Let d(t1) represent the time delay of the fracture in the i-th formation, d(t1) represent the propagation vector received on the detector, and R(θ) represent the Alford rotation matrix.

[0018] Optionally, the azimuth domain common imaging gather is subjected to fast and slow shear wave separation to determine the fast shear wave gather, the slow shear wave gather, and the angle between the fracture orientation in each formation fracture and the main survey line direction of the geophone. This also includes:

[0019] The azimuth domain common imaging gathers for each formation fracture, from shallow to deep, are determined in the direction of the main survey line and the direction of the connecting line.

[0020] The azimuth domain common imaging gathers in the main survey line direction and the connecting line direction are rotated to the formation fractures to determine the fast shear wave and slow shear wave corresponding to each formation fracture.

[0021] The scanning angle is used to rotate at different offset distances on the azimuth domain common imaging gather. In response to the maximum energy ratio of the fast shear wave and the slow shear wave at the scanning angle, the scanning angle is determined to be the relative angle between the dominant development azimuth of the current formation fracture and the fracture orientation of the current formation fracture at the current offset distance.

[0022] Based on the relative angle, the angle between the direction of the crack and the main measuring line direction of the detector is determined.

[0023] Optionally, the offset distance includes near offset distance and mid-to-far offset distance;

[0024] Based on the fast shear wave gather and the slow shear wave gather, the development direction of fractures in the target formation is determined, including:

[0025] The average value of the circular offset is determined to be the scanning angle;

[0026] Determine a first weight corresponding to the mid-to-long offset distance and a second weight corresponding to the near offset distance; wherein the first weight is greater than the second weight.

[0027] The fracture development direction of the target stratum is determined based on the first weight, the second weight, and the circular mean.

[0028] Optionally, the azimuth domain co-imaging gathers corresponding to the fast and slow shear waves are converted into snail gathers according to the azimuth and offset of the fan-shaped region, and the time difference between the fast and slow shear wave gathers in each azimuth is corrected, including:

[0029] The fast shear wave gather and the slow shear wave gather are converted into snail gathers according to the orientation and offset of the sector region;

[0030] Using the near offset of the snail track set as a reference track, a standard reference overlay profile model is generated.

[0031] The remaining traces in the snail trace set are taken as traces to be corrected. The traces to be corrected are cross-correlated with the reference overlay profile model to determine the time shift corresponding to the largest correlation.

[0032] The time difference of the channel to be corrected is determined based on the time shift.

[0033] Optionally, based on the time difference, the density of the cracks is determined, and the crack development index is determined according to the density, including:

[0034] Based on the time difference, determine the magnitude of the time difference between the fast and slow shear waves, and use the sliding time window Z-score to determine the local mean and standard deviation of the residuals of the time difference magnitude within the sliding time window;

[0035] The crack development index is determined based on the time difference, the local mean, and the standard deviation.

[0036] Optionally, the circular mean of the offset distance is represented as:

[0037]

[0038] Where j = 1, 2, ..., N i N represents the sampling points within the time window, and θ i,j This indicates all scan angles within the time window.

[0039] Optionally, the fracture development direction of the target stratum is represented as follows:

[0040]

[0041] Where M represents the maximum offset distance, ω i These represent the weighting coefficients for different offset distances.

[0042] Optionally, the time difference of the channel to be corrected is expressed as:

[0043]

[0044] Among them, s j (t) represents the signal of the j-th near offset channel, m represents the number of channels stacked, and S m (t) represents the reference superimposed profile model trace signal, s i (t) represents the signal of the i-th channel to be corrected. and t1 and t2 represent the average values ​​of the traces in the reference superimposed profile model and the traces to be corrected, respectively, and [t1, t2] represent the time window for cross-correlation calculation.

[0045] A second aspect of this application provides a crack prediction device based on seismic data, comprising:

[0046] The imaging processing module is used to divide the seismic data into multiple sector regions according to the azimuth of the shot point and the detector, and to determine the azimuth domain common imaging gather of the converted wave in each sector region.

[0047] The wave field separation module is used to separate the fast and slow shear waves in the azimuth domain common imaging gather, and to determine the fast shear wave gather, the slow shear wave gather, and the angle between the fracture orientation and the main survey line direction of the geophone in each formation fracture.

[0048] The parameter inversion module is used to determine the fracture development direction of the target formation based on the fast shear wave gather and the slow shear wave gather.

[0049] The time difference correction module is used to convert the azimuth domain co-imaging gathers corresponding to the fast shear wave and slow shear wave into snail gathers according to the azimuth and offset of the fan-shaped region, and to correct the time difference of each channel in the snail gathers.

[0050] A crack development module is used to determine the density of the cracks based on the time difference, and to determine the crack development index based on the density.

[0051] As can be seen from the above, the crack prediction method and device based on seismic data provided in this application adopts a sub-azimuth processing flow based on the azimuth of the shot receiver. Especially for three-dimensional seismic data with uneven azimuth coverage, it shows better applicability and stability, and makes the most of the data information acquired by sub-azimuth.

[0052] By analyzing the reflected amplitude energy and travel time characteristics of converted waves in different azimuths, we constrain the dominant azimuth of fracture development from shallow to deep. Secondly, we perform layer-by-layer Alford rotation on the azimuth domain common imaging gather (ACIG) to achieve layer-by-layer separation of fast and slow shear waves. Then, considering that the mid-to-far offset is more sensitive to fracture response and that the angle has periodic characteristics, we assign appropriate weighting factors to different offsets and use weighted circular mean to calculate the main fracture direction of the target layer.

[0053] A Z-score normalization method based on a sliding time window is introduced to characterize the spatial variation trend of fracture density, highlighting the local changes in anomalous splitting time difference, thereby enhancing the ability to identify fracture enrichment zones. Compared with the traditional method of characterizing fracture density using the time difference of fast and slow shear waves, the Z-score method has stronger resistance to background interference, and is particularly stable in regions with strong tectonic undulations.

[0054] The above description is merely an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention and to implement it in accordance with the contents of the specification, and in order to make the above and other objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention are described below. Attached Figure Description

[0055] To more clearly illustrate the technical solutions in this application or related technologies, the drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0056] Figure 1 This is a flowchart of the crack prediction method 100 based on seismic data according to an embodiment of this application;

[0057] Figure 2 This is a schematic diagram of a multiple shear wave separation phenomenon according to an embodiment of this application;

[0058] Figure 3 for Figure 2 Schematic diagram of the travel time of PS1 and PS2 waves with azimuth in the second layer of the crack;

[0059] Figure 4 for Figure 2Schematic diagram of the travel time of PS1 and PS2 waves with azimuth in the third layer of the fracture;

[0060] Figure 5a This is a schematic diagram of a reference coordinate system according to an embodiment of this application;

[0061] Figure 5b This is a schematic diagram of the polarization direction for horizontal X-component rotation according to an embodiment of this application;

[0062] Figure 6a This is a schematic diagram showing the two cracks with an included angle of 0° in an embodiment of this application;

[0063] Figure 6b This is a schematic diagram showing a 10° angle between two cracks in an embodiment of this application.

[0064] Figure 6c This is a schematic diagram showing a 20° angle between two cracks in an embodiment of this application.

[0065] Figure 6d This is a schematic diagram showing a 30° angle between two cracks in an embodiment of this application.

[0066] Figure 6e This is a schematic diagram showing the angle between the two cracks in an embodiment of this application, which is 40°.

[0067] Figure 6f This is a schematic diagram showing a 50° angle between two cracks in an embodiment of this application.

[0068] Figure 7 This is a schematic diagram of the azimuth sector division of the shot point-detector in an embodiment of this application;

[0069] Figure 8 This is the azimuth domain common imaging trace set for the R and T components in this embodiment of the application;

[0070] Figure 9a for Figure 2 The time window where the second layer of cracks is located;

[0071] Figure 9b for Figure 2 The time window in which the third layer crack is located;

[0072] Figure 10 These are the fast shear wave gathers and slow shear wave gathers obtained using the Alford rotation method in the embodiments of this application;

[0073] Figure 11 This is a schematic diagram of a snail gating system according to an embodiment of this application;

[0074] Figure 12 This is a schematic diagram of the fast shear wave gather and slow shear wave gather before and after time difference correction for fast and slow shear waves in an embodiment of this application.

[0075] Figure 13 This is a schematic diagram of the coverage number and azimuth distribution of a shot point-detector according to an embodiment of this application;

[0076] Figure 14 This is a schematic diagram comparing the effects of fast and slow shear wave separation before and after on the azimuth domain common imaging gather in an embodiment of this application.

[0077] Figure 15a This is a crack principal direction rose diagram obtained from the dominant orientation of crack development according to an embodiment of this application;

[0078] Figure 15b This is another crack principal direction rose diagram obtained from the crack development dominance orientation according to an embodiment of this application;

[0079] Figure 16 This is a schematic diagram of the time difference correction for fast and slow shear waves in an embodiment of this application.

[0080] Figure 17 This is a schematic diagram of the R and T components of a converted wave according to an embodiment of this application;

[0081] Figure 18 for Figure 17 A schematic diagram of the imaging area within the blue box;

[0082] Figure 19 This is a schematic diagram illustrating crack density characterization according to an embodiment of this application;

[0083] Figure 20 This is another schematic diagram illustrating crack density characterization according to an embodiment of this application;

[0084] Figure 21 This is a schematic diagram of a crack prediction device based on seismic data according to an embodiment of this application;

[0085] Figure 22 This is a schematic diagram of an electronic device according to an embodiment of this application. Detailed Implementation

[0086] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.

[0087] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this application should have the ordinary meaning understood by one of ordinary skill in the art to which this application pertains. The terms "first," "second," and similar terms used in the embodiments of this application do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are only used to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0088] As described in the background section above, shear wave splitting is a core phenomenon and analytical method for characterizing the anisotropy of subsurface media, especially fracture characteristics. Shear wave splitting refers to the phenomenon where, when a shear wave passes obliquely across a fracture surface, it splits into a fast shear wave and a slow shear wave, with polarization directions parallel and perpendicular to the fracture plane, respectively. The detector receives a mixture of partial projections of the fast and slow shear waves. If these two waves are not separated, the accuracy of the analysis of azimuth anisotropy variations will be affected.

[0089] Oriented vertical fractures in the formation induce significant azimuthal anisotropy in the medium. The ray path of the converted wave incorporates both P-waves and S-waves, including not only the azimuthal anisotropy information of the P-wave, but more importantly, the shear wave splitting phenomenon provides richer seismic responses for fracture parameter inversion, which is more conducive to improving the accuracy of fracture identification and the reliability of quantitative characterization. However, related techniques for shear wave splitting analysis of PS-waves are still limited to post-stack domain processing, failing to fully utilize pre-stack information.

[0090] When using 3D seismic data for pre-stack fracture prediction, dividing the data into azimuth sectors is a common strategy for extracting azimuth anisotropy information. Two main technical approaches exist: OVT domain processing and azimuth-based migration processing based on shot-receiver azimuth. OVT processing, by merging seismic traces with the same migration vector (including offset and azimuth) into OVT gathers, theoretically provides complete azimuth information for wide-azimuth seismic data. However, OVT domain processing faces numerous challenges in practical applications: ideal OVT domain data acquisition is costly and difficult to fully realize; the irregular spatial distribution of shot and receiver points in conventional observation systems often leads to severely uneven OVT domain coverage, inconsistent gather quality, and the introduction of noise problems such as acquisition footprints and spatial aliasing, typically requiring complex regularization techniques for compensation.

[0091] Based on this, refer to Figure 1As shown, this application embodiment provides a crack prediction method 100 based on seismic data, including: S101, dividing the seismic data into multiple sector regions according to the shot point and the azimuth angle of the detector, and determining the azimuth domain common imaging gather of the converted wave in each sector region.

[0092] In some exemplary embodiments, the seismic data in this application is mainly OBN (Ocean Bottom Node) seismic data. OBN seismic data is a technology used for marine seismic exploration. It is usually deployed on the seabed and can independently acquire and record seismic signals.

[0093] In some exemplary embodiments, the converted wave is a PS wave.

[0094] First, the directional vertical fractures in the strata induce significant azimuth anisotropy in the medium, resulting in a clear azimuth dependence in the propagation characteristics of seismic waves. The PP wave travel time characteristics based on azimuth angle can be expressed as:

[0095]

[0096] Where T and T0 represent the travel time of the PP wave at offset x and zero offset, respectively, α represents the azimuth angle of the CMP line, and V nmo η represents the NMO velocity, and η represents the non-elliptic coefficient.

[0097] It should be noted that CMP stands for Common Center Point Gather, and NMO stands for Normal Time Difference Correction, which is used to eliminate the difference in seismic wave arrival time caused by different shot-receiver distances.

[0098] Furthermore, the travel time equation for PS waves in VIT (Vertical Transverse Isotropy) media is as follows:

[0099]

[0100] Among them, T P0 V represents the travel time of a downlink P-wave with zero offset in a single trip. p and η P T represents the downlink P-wave NMO velocity and non-elliptic coefficients, respectively. S0 V represents the travel time of an upward S-wave with zero offset in a single trip. S and η S These represent the upward S-wave NMO velocity and non-elliptic coefficients, respectively.

[0101] like Figure 2As shown, when converted waves propagate in strata containing multiple vertical fractures, the propagation characteristics of fast and slow shear waves become extremely complex due to multiple shear wave splitting effects, generating PS1 (fast shear wave) and PS2 (slow shear wave) waves with different propagation speeds and polarization directions. The records received by the geophone may contain superposition and interference of multiple sets of fast and slow shear waves from deeper layers, which increases the difficulty of shear wave splitting analysis. Based on the shot-receiver azimuth, the seismic data is divided into multiple sector regions at certain intervals, and pre-stack time migration is performed on each azimuth sector region to generate the azimuth-domain common-image gather (ACIG) of PS waves.

[0102] A horizontal layered model was further designed, with a model size of 1000m×800m×2000m. The first and fourth layers are isotropic media (counting from top to bottom), with true north as zero degrees, counterclockwise as positive direction, and clockwise as negative direction. The crack direction in the second layer is 50°, and the crack direction in the third layer is 168°. The model parameters are shown in Table 1.

[0103] Table 1

[0104]

[0105] Figure 3 The travel times of PS1 and PS2 in the second fractured medium are shown to exhibit a clear elliptical shape in polar coordinates. The PS1 wave has the fastest velocity along the fracture strike (red dashed line), corresponding to the minor axis of the blue ellipse, and the slowest velocity perpendicular to the fracture strike (major axis). Conversely, the PS2 wave has the slowest velocity along the fracture strike (major axis of the green ellipse) and the fastest velocity perpendicular to the fracture strike (minor axis). Furthermore, the time delay Δt between PS1 and PS2 reaches its maximum value in the fracture strike direction.

[0106] exist Figure 4 In the middle layer, the travel time of the third layer is affected by the overlying fractured layer, causing the travel time patterns and time delay Δt of PS1 and PS2 to no longer align with the fracture strike (red dashed line). This indicates the superposition of multiple anisotropies. The fractures in the overlying strata alter the time delay and polarization direction, obscuring the fracture strike of the target layer. Single-layer elliptic fitting becomes insufficient to solve the problem, thus requiring layer-by-layer anisotropy correction. To address this issue, subsequent steps will introduce how to separate the fractured shear wave modes through layer peeling Alford rotation and determine the rotation angle to infer the fracture development direction of each layer.

[0107] In step S102, fast and slow shear wave separation is performed on the common imaging gather in the azimuth domain to determine the fast shear wave gather, the slow shear wave gather, and the angle between the direction of the crack in each bottom layer crack and the main measurement line direction of the detector.

[0108] Specifically, step S102 includes determining the propagation vector composed of the fast shear wave and the slow shear wave, using the main measurement line direction and the connecting line direction of the detector as the X-axis and Y-axis of the propagation coordinate system, so that the propagation vector is transformed into a coordinate system composed of the polarization directions of the fast shear wave and the slow shear wave. The direction of the detector's X component is the Inline direction, and the direction of the Y component is the Xline direction. The reference coordinate system and the detector's coordinate system (propagation coordinate system) are the same coordinate system.

[0109] Specifically, using the inline (main survey line direction) and Xline (connection line direction) of the geophone as the X and Y axes of the coordinate system (the geophone's coordinate system), the converted wave undergoes shear wave splitting when it obliquely crosses a vertically fractured formation. Let S(t) be the initial waveform of the converted wave, and θ represent the angle between the polarization direction of the fast shear wave and the X-axis. Then, in this coordinate system, the propagation vector S after initial polarization... P (t) transformed into the propagation coordinate system is expressed as:

[0110]

[0111] During propagation, a time difference arises between the fast and slow shear waves due to their velocity difference, which can be described by the diagonal time delay Λ(t), where the diagonal elements of the diagonal time delay are the delays of the fast and slow shear waves, respectively. Correspondingly, the propagation vector d composed of the fast and slow shear waves... P Represented as:

[0112] d P (t)=Λ(t)*S P (t), (4)

[0113] Determine the reference coordinate system in which the detector is located, and rotate the propagation vector composed of fast and slow shear waves into the reference coordinate system to determine the seismic records after multiple shear wave separations of fractures in different strata under the reference coordinate system.

[0114] To obtain the signals recorded in the detector coordinate system, i.e., the reference coordinate system where the detector is located (horizontal X and Y components), namely the PSV wave and PSH wave, respectively. The PSV wave represents the signal of the horizontal X component in the reference coordinate system, that is, the mixed signal of the converted wave projected onto the horizontal X component, and the PSH wave represents the signal of the horizontal Y component in the reference coordinate system, that is, the mixed signal of the converted wave projected onto the horizontal Y component. Next, the wavelength after time delay needs to be rotated from the propagation coordinate system back to the reference coordinate system. Using the Alford rotation matrix R(θ), the propagation vector d(t) finally received at the detector can be expressed as:

[0115] d(t)=R T d P (t)=RT (θ)Λ(t)*S P (t), (5)

[0116] When a converted wave propagates in a medium containing multiple vertical cracks, it undergoes multiple transverse wave splitting. Taking a two-layer crack as an example, neglecting amplitude attenuation, in time window t1, the fast and slow transverse waves received by the horizontal X and Y components from the first crack are:

[0117]

[0118] Where θ1 represents the angle between the orientation of the first fracture layer and the direction of the main survey line (X-axis), and Λ1 represents the time delay of the first fracture medium. When the reflected shear wave from the second fracture medium passes through the upper fracture medium, a secondary shear wave split will occur. During this process, the reflected wave energy in the second fracture medium can be considered equivalent to a secondary seismic source. In time window t2, the PSV and PSH waves received on the horizontal X and Y components can be expressed as:

[0119]

[0120] Where θ2 represents the angle between the direction of the second-layer fracture and the main survey line (X-axis), and Λ2 represents the medium delay of the second fracture. Similarly, the seismic records received on the horizontal X and Y components after multiple shear wave separations can be represented as follows:

[0121]

[0122] Among them, PSV(t) i ) and PSH(t i ) represent the seismic records of multiple shear waves received by the fracture in the i-th stratum in the X and Y components of the reference coordinate system after separation, t represents the propagation time, T represents the transpose of the rotation matrix, and Λ i Let represent the time delay of the fractured medium in the i-th formation, d(t1) represent the propagation vector received on the detector, and R(θ) represent the Alford rotation matrix.

[0123] Furthermore, for formations containing multiple fractures, layer stripping is required to separate fast and slow shear waves. The time window range for shear wave splitting is selected. From shallow to deep, the azimuth domain common imaging gathers for each formation fracture are determined in both the main survey line direction and the shear line direction.

[0124] During field observations, observations are generally taken along the survey line direction X and perpendicular to the survey line direction Y. For subsequent processing, the X and Y components first need to be rotated to the direction of the shot-detector line (R) and the direction orthogonal to it (T). Here, the azimuth α formed by the shot pointing towards the detector and the X-axis is used to rotate the X and Y components to the R and T components, as shown below. Figure 5a and Figure 5bAs shown.

[0125] The azimuth domain common imaging gathers in the main survey line direction and the connecting line direction are rotated to the formation fractures to determine the fast and slow shear waves corresponding to each formation fracture.

[0126] Specifically, from shallow to deep strata, Alford rotation separation of fast and slow shear waves is performed on the ACIG gathers of the R component (along the shot-to-detector line) and T component (perpendicular to the shot-to-detector line) to separate the fast and slow shear waves. Taking a model containing two fractured media as an example, firstly, based on geological data and the travel time and reflection energy characteristics of PS waves on the ACIG, the range where shear wave splitting occurs is defined. The first fractured layer on the ACIG gathers of the R and T components is taken as the target layer, and the fast shear wave PS1 of the first fractured layer is... (1) Slow transverse wave PS2 (1) In ACIG sets, this can be represented as:

[0127]

[0128] The smaller the relative angle between the azimuth α and the orientation θ1 of the first layer of fracture, the earlier the in-phase axis of the R component arrives, and the greater the energy ratio of the in-phase axes of R and T. At this time, α is the dominant azimuth for the development of the first layer of fracture, and it is used as prior information.

[0129] Rotation is performed at different offsets on the azimuth domain common imaging gather using a scanning angle. In response to the maximum energy ratio of fast shear wave and slow shear wave at the scanning angle, the scanning angle is determined to be the relative angle between the dominant azimuth of the current formation fracture and the fracture orientation of the current formation fracture at the current offset.

[0130] At different offsets on the ACIG, a series of scanning angles κ are used for rotation. When the scanning angle κ satisfies the maximum energy ratio of fast shear waves and slow shear waves, the scanning angle at this time is the relative angle between the dominant azimuth and the strike of the first-layer fracture at the current offset.

[0131]

[0132] Based on the relative angle, determine the angle between the crack direction and the main measuring line direction of the detector.

[0133] Therefore, for ACIG in different orientations, the angle between the crack orientation and the main survey line direction is: θ1=α+κ (1) offset Alford rotation was used to rotate the ACIG of the R and T components to the fast and slow shear waves of the first fracture layer.

[0134] After correcting for the shear wave splitting effect of the first fracture layer, the second fracture layer on the ACIG gather was selected as the target layer, and the fast shear wave PS1 to be corrected was... (2)’ and slow transverse wave PS2 (2)’ It can be represented as:

[0135]

[0136] The fast transverse wave PS1 to be corrected at this time (2)’ and slow transverse wave PS2 (2)’ This reflects the shear wave splitting information of the current fractured strata. Referring to the wavefield separation method of the first fractured strata, the azimuth with the earliest arrival time of the in-phase axis on the R component and the highest R and T energies is selected as the preferred method for the development of the second fracture layer. The angle between the travel time of the second fracture layer and the preferred azimuth is scanned at each offset. When the energy ratio of fast shear wave to slow shear wave is the largest, the angle between the azimuth and strike of the second fracture layer and the direction of the main survey line is: θ2=α+κ (2) offset Therefore, the parameters of each crack layer must be obtained by analyzing and correcting the shear wave splitting effect of each layer from top to bottom in order to obtain the true and accurate rotation angle.

[0137] To accurately characterize the propagation and superposition effects of PS waves (converted waves) in multi-layered fractured media, the angle difference between adjacent fractured strata must be considered. When the difference is small (less than 20°), in Figures 6a-6b In the 1000-1200 ms window, the X component of the second-layer crack shows a distinct single peak, and the vector tip plot shows a dominant polarization direction, reflecting a nearly uniform shear wave splitting similar to a single-layer effect. However, in Figure 6c , Figure 6d , Figure 6e and Figure 6f In the analysis, as the angle difference increases (greater than 20°), the X component exhibits multiple peaks and troughs within the same time window, while the vector tip plot reveals two main polarization directions, indicating the interference and superposition of shear waves from layers with different fracture orientations. This suggests that the shear wave splitting and accumulation effects become more complex, thus requiring layer stripping to accurately determine the fracture orientation of each layer. Figures 6a-6f In Figure 6, the left column shows the simulated X-component (black curve) and Y-component (red curve) records with 10% random noise added, and the right column shows the vector plots of the X-component and Y-component for a time window of 1000-1200 ms. In Figure 6(a), Δθ = 0°; in (b), Δθ = 10°; in (c), Δθ = 20°; in (d), Δθ = 30°; in (e), Δθ = 40°; and in (f), Δθ = 50°.

[0138] Subsequently, in S103, the development direction of fractures in the target stratum was determined based on the fast shear wave gather and the slow shear wave gather.

[0139] Since the orientation of fractures in the formation changes with depth from shallow to deep, and each offset within the time window will generate one or more scanning angles representing the fracture response, it is necessary to extract the main orientation of the overall fracture development from these angles.

[0140] Let α be the dominant orientation of fracture development in the i-th stratum. i Then, the principal direction θi of the crack in the i-th layer satisfies θi∈[α]. i -Δφ,α i +Δφ], where Δφ is the half-width of the azimuth interval. The azimuth interval is as follows: Figure 2 The various sector-shaped regions shown.

[0141] Within the time window of the fracture layer, at each offset... i (i = 1, 2, ..., M), where M is the maximum offset distance, and a set of scan angle data {κ} is obtained. i,j}, j = 1, 2, ..., N i N represents the sampling points within the time window.

[0142] In some embodiments, factors such as uneven actual signal-to-noise ratio and coverage times can lead to abnormal scanning angles, when the scanning angle κ... i,j If the value is greater than 2Δφ, the data is considered an outlier and is removed.

[0143] In some embodiments, step S103 includes:

[0144] The circular homogeneous element with the offset distance is used as the scanning angle.

[0145] Angular data is periodic, making conventional arithmetic averaging inapplicable. For a set of scan angle data {κ} within the time window containing the fracture in the i-th formation, i,j For each offset after removing outliers, i The average value of the circle at this offset distance is calculated as the scanning angle.

[0146] In some embodiments, the circular mean of the offset distance is represented as:

[0147]

[0148] Where j = 1, 2, ..., N i N represents the sampling points within the time window, and θ i,j This indicates all scan angles within the time window.

[0149] Its circular standard deviation is:

[0150]

[0151] Where R represents the length of the mean vector, a larger R indicates that the data is more concentrated, σ θi Smaller; conversely, when R is smaller, the data is more dispersed, σ θi Relatively large.

[0152] Determine the first weight corresponding to the mid-to-long offset distance and the second weight corresponding to the near offset distance; wherein the first weight is greater than the second weight.

[0153] In PS-wave (converted wave) PSTM (PreStack Time Migration), the accuracy of the conversion point calculation has a significant impact on imaging. This conversion point is the point where the P-wave converts to the S-wave.

[0154] At near offsets, the conversion point is close to the detector, causing the ray path to tilt and the in-phase axis to become unfocused. At far offsets, the conversion point is closer to the reflection point, the in-phase axis is focused, and the reflected energy is stronger. Regarding the changes in ray path geometry and focusing effect of PS waves at different offsets, this application's embodiments assign larger weighting factors (first weight) to mid-to-far offsets and offsets with high signal-to-noise ratios, and smaller weighting factors (second weight) to near offsets and offsets with low signal-to-noise ratios.

[0155] The direction of fracture development in the target stratum is determined based on the first weight, the second weight, and the circular mean.

[0156] The principal direction of fractures in this layer is calculated using a weighted circular mean. Multiplying the circular mean of all offsets by their respective weights, the fracture development direction of the target stratum is then expressed as:

[0157]

[0158] Where M represents the maximum offset distance, ω i These represent the weighting coefficients for different offset distances.

[0159] Based on the model parameters in Table 1, forward modeling was performed using the reflectivity method with 10% random noise added to verify the effectiveness of the method. The ground seismic observation system consisted of 357 shots with a detector spacing of 25 meters and a shot point spacing of 50 meters. Figure 7 At the center of the black triangle marked (500, 400), the azimuth angle of the shot point to the detector is divided into six sectors. Figure 8 The PS-wave azimuth gathers (ACIGs) of the R and T components obtained at this imaging location using pre-stack time migration (PSTM) are shown.

[0160] Figure 8In the middle (a), at approximately 820ms, indicated by the red arrow at azimuth 2 in the first row of R component, the PS wave reflection energy is relatively strong and the time is shortest, while the PS wave reflection energy at azimuth 5 is relatively strong and the time is longest. Figure 8 In (b), the T component of the second row shows the opposite trend: the PS wave reflection energy at azimuth 2 and azimuth 5 is relatively weak, and the arrival time at azimuth 5 is longer. Around 1220 ms (as shown by the red arrow), due to the influence of the overlying fracture layer, the shear waves in the R and T components undergo secondary splitting. The interference between the PS waves forms a distinct in-phase axis, exhibiting multiple in-phase axes in different azimuths. The waveform broadens, and two adjacent peaks can be observed in some azimuths. The split shear waves of the first layer are referred to as PS1. (1) Wave and PS2 (1) The transverse waves that split twice in the second layer are called PS. 11 (2) Wave, PS 12 (2) Wave, PS 21 (2) Wave and PS 22 (2) Wave.

[0161] Targeting the first crack layer, a scan was performed within a time window of 820-920ms. Figure 9a The angle between the principal direction of azimuth 2 and the strike of the first-layer crack was determined to be approximately 4.4°. Based on this, the R and T components of azimuth 2 were rotated, successfully eliminating the in-phase axis of the PS wave near 820 ms. Figure 10 This indicates that the crack location was accurately estimated and that the leaked PS1 (1) Energy is redirected to the PS1 wave. After correcting for shear wave splitting in the first layer, the target is the second fracture layer. Scan ( Figure 9b The angle between the main azimuth of azimuth 6 and the direction of the second layer crack was determined to be approximately 3.5°. The azimuth of the second layer was then rotated to the final PS1 and PS2 waves, thereby effectively separating the PS1 and PS2 waves. Figure 10This paper presents the PS1 and PS2 ACIG gathers obtained after applying the Alford rotation method with layer stripping. This method statistically separates the fast and slow shear waves, enhances the PS1 energy and arrival time, and weakens the PS2 energy by eliminating fast wave interference, thus making the overall fast and slow shear wave time differences clearer. Although secondary splitting information may be difficult to extract due to wavefield interference and energy attenuation in multi-layered fractures, this method significantly improves the PS1 / PS2 separation effect. It is worth noting that when the wavefield passes through multiple or even more layers of fractured media, interference and superposition of the wavefield in different azimuths, especially the rapid energy attenuation of the secondary wavefield generated by secondary splitting, makes it difficult to effectively separate independent wavefield components, thus making it difficult to stably extract secondary splitting information. Nevertheless, this method significantly improves the PS1 / PS2 separation effect, thus making the overall fast and slow shear wave time differences clearer. However, since the travel times still differ in different azimuths, further correction of the remaining time differences is still needed.

[0162] Subsequently, in step S104, the azimuth domain co-imaging gathers corresponding to the fast and slow shear waves are converted into snail gathers according to the azimuth and offset of the fan-shaped region, and the time difference of each channel in the snail gathers is corrected.

[0163] Specifically, step S104 includes:

[0164] The fast shear wave gather and the slow shear wave gather are converted into snail gathers according to the azimuth and offset of the sector region.

[0165] like Figure 11 As shown, the PS1 and PS2 wave gathers after separating the fast and slow shear waves are converted into snail gathers according to their azimuth and offset. The remaining time difference is calculated using the normalized cross-correlation method. After correcting the time difference in each direction, the PS1 and PS2 wave profiles are obtained by superimposing them.

[0166] Using the near offset of the snail track set as a reference track, a standard reference overlay profile model is generated.

[0167] Each layer of cracks has its own dominant crack development direction, and the near offset exhibits weak anisotropy. Therefore, the near offset of the snail gather can be used as a reference trace to generate a standard reference stacked profile model:

[0168]

[0169] The remaining traces in the snail trace set are taken as traces to be corrected. The traces to be corrected are cross-correlated with the reference superimposed profile model to determine the time shift corresponding to the largest correlation.

[0170] In some exemplary implementations, the time shift refers to the distance the signal is shifted along the time axis.

[0171] For each arrival to be corrected, cross-correlation is performed with the reference overlay profile model to find the time shift corresponding to the largest correlation, thereby estimating and correcting the remaining time difference of each arrival to be corrected as follows:

[0172]

[0173] Based on the time shift, determine the time difference of the track to be corrected.

[0174] Determine the time difference Δt to be corrected i for:

[0175] Δt i =argmaxC i (τ),

[0176] Right now

[0177] Among them, s j (t) represents the signal of the j-th near offset channel, m represents the number of channels stacked, and S m (t) represents the reference superimposed profile model trace signal, s i (t) represents the signal of the i-th channel to be corrected. and t1 and t2 represent the average values ​​of the traces in the reference superimposed profile model and the traces to be corrected, respectively, and [t1, t2] represent the time window for cross-correlation calculation.

[0178] Then in step S105, the density of the cracks is determined based on the time difference, and the development index of the cracks is determined based on the density.

[0179] Specifically, step S105 includes:

[0180] Based on the time difference, the magnitude of the time difference between the fast and slow shear waves is determined. Using the sliding time window Z-score, the local mean and standard deviation of the residuals of the time difference magnitude within the sliding time window are determined.

[0181] The time difference between fast and slow shear waves is related to the degree of fracture development, and is usually expressed as the time difference between fast and slow shear waves. This method characterizes the density of fractures. However, it is easily affected by background time variations caused by tectonic factors such as faults and stratigraphic dips, which can lead to the fracture response being masked or misjudged.

[0182] Therefore, the crack density is characterized by Z-score normalization using a sliding time window, and the time difference δt is adjusted accordingly. i Perform dimensionless processing.

[0183] The crack development index is determined based on the time difference, local mean, and standard deviation.

[0184]

[0185] Where i represents the CMP point, and These represent the mean and standard deviation of the residuals within the window during sliding, respectively.

[0186] In some embodiments, a low degree of crack development is determined in response to a crack development index less than 0; that is, e < 0.

[0187] If the crack development index is greater than 0, then the crack development level of the current area is determined to be high. That is, e > 0.

[0188] When the crack development index approaches 0, it is determined that cracks exist in the current region and are evenly distributed. That is, when e≈0.

[0189] like Figure 12 As shown in (a) and (b), the PS1 and PS2 waves, separated by azimuth, are converted into a snail gather using ACIG. The first layer at approximately 430 ms is an isotropic medium, with PS waves exhibiting consistent orientations and straight phase axes in all directions. The second layer at approximately 820 ms is a PS1 wave... (1) Wave( Figure 11 (a) and PS2 (1) Wave( Figure 11 In (b) the in-phase axis exhibits wavy jitter. At approximately 1220 ms, the third-layer PS wave undergoes secondary splitting, resulting in attenuation of wave field energy. The PS wave in the snail channel of PS1 is concentrated in the PS... 11 (2) Wave and PS 21 (2) Wave interference forms a strong reflection axis ( Figure 12 In (a), the PS2 wave in the snail track concentration of PS 12 (2) Wave and PS 22 (2) The large time difference in wave arrival causes the wave field to lose its in-phase axis characteristics. Figure 12 (b)

[0190] After correcting for time differences in all directions, Figure 12 In (c)-(d), PS1 is located at the second layer. (1) Wave and PS2 (1) The wave's phase axis is essentially flat, and the third layer also forms a statistically strong phase axis. The flattened phase axis does not strictly correspond to the PS. 11 (2 Wave, PS 21 (2) Wave, PS 12 (2) Wave and PS 22 (2)Instead of representing the actual arrival time of the waves, it forms an in-phase axis that effectively characterizes the arrival times of PS1 and PS2 waves, which can effectively characterize the overall travel time of the PS1 and PS2 wave fields in the third layer. The model is horizontally layered and does not need to consider the background time difference caused by the correction structure. The fast and slow shear wave time differences of the second and third layers picked on the superimposed PS1 and PS2 wave profiles are 23.97 ms and 28.14 ms, respectively. The crack distribution of the horizontally layered model is uniform, and the overall Z-score distribution will be concentrated near 0.

[0191] In some exemplary embodiments, the fracture prediction method based on OBN seismic data from this application is applied to the fast and slow shear wave separation and fracture prediction of multi-component seismic data in a public area of ​​a certain sea area. The public area seismic data was acquired using a conventional three-dimensional towed cable method, and its shot-receiver coverage times and azimuth distribution are as follows: Figure 13 As shown in the figure, the coverage area is concentrated in the azimuth range of 0°-30° and 180°-210°. A azimuth-based processing strategy was adopted, with true north as 0° and counterclockwise as the positive direction. First, the entire azimuth range (0°-360°) was divided into twelve initial sector regions at 30° intervals. Then, the diagonal sector regions were merged, integrating the data into six azimuth datasets. Next, PSTM processing was performed on each of these six azimuth datasets to generate the azimuth domain common imaging gather (ACIG).

[0192] Figure 14 Images (a)-(b) show the R and T components before the separation of fast and slow shear waves. Figure 14 Figures (c)-(d) show the fast and slow shear wave components after separation of the fast and slow shear waves. Figure 14 (a)-(b) show the ACIG of the R and T components at imaging point 2000 in six azimuth directions. The reflected energy of the wavefield is basically uniform in the shallow layer, indicating that the shallow layer has a low degree of fracture development and is close to an isotropic medium. The blue dashed box at approximately 3400-3600 ms represents the time window of fracture development in group W. Within this window, the R component (blue arrow) in azimuth sector 6 shows the strongest energy, while azimuth sector 3 shows the weakest energy. After PS wave splitting correction, Figure 14 In the blue box in (c), the reflected energy of the PS1 wave is enhanced. Figure 14 The reduced wavefield reflection energy within the blue box in (d) indicates that the leaked PS1 wave in the T component has refocused into the PS1 wave. Further analysis of the fracture's principal orientation at the dominant azimuth reveals that the reflection energy of the seismic traces at the mid-to-far offset is stronger. With a weighting factor set to 0.6 and the others to 0.4, the principal orientation of the fracture in group W is approximately 158°. Figure 15aAfter correcting for the shear wave splitting effect of the overlying strata, the target layer was set at approximately 4000 ms, indicated by the blue dashed line. Within this window, the R component of azimuth sector 2 showed the strongest energy (blue arrow), while azimuth sector 5 showed the weakest energy (blue arrow). The T component showed the opposite trend. After correction, Figure 14 In the middle (c), the amplitude energy of the PS1 wave, indicated by the blue dashed line, is enhanced. Figure 14 The PS1 wave leaking in (d) was separated, and the amplitude energy of the PS1 wave was weakened, indicating that the fast and slow shear waves were effectively separated. Further statistical analysis of the main fracture direction in the dominant azimuth yielded that the main fracture azimuth in group L was approximately 51°. Figure 15b ).in Figure 15a The blue box in the middle indicates the main direction of the cracks in group W. Figure 15b The blue dashed line indicates the main direction of the cracks in group L.

[0193] like Figure 16 Images (a) and (b) show the snail gathers of the fast and slow shear waves before correction, respectively. The ACIGs of PS1 and PS2 waves are converted to snail gathers. The snail gathers show a significant azimuth change in the travel time of the in-phase axis, indicating the presence of anisotropy caused by a horizontally transverse isotropy (HTI) medium. Figure 16 In (c)-(d), after the residual time difference is corrected by the snail gather of PS1 and PS2 waves, the periodic jitter of the in-phase axis is improved, and the in-phase axis indicated by the red arrow is more focused. Figure 17 Images (a)-(b) show PS-wave imaging profiles after superimposing snail gathers of the R and T components. Green represents the fourth segment of the W group, and blue represents the L group. Green lines indicate faults, and blue boxes indicate fracture interpretation areas. This stratum underwent a Paleogene rifting phase, depositing terrestrial clastic rocks including the L and W groups. The fault system on the right side is complex, dominated by extensional structures, with local development of Y-shaped extensional-strike-slip structures and compressional structures such as folds and truncation unconformities. Continued fault activity is often accompanied by fracture development; the blue boxes indicate predicted fracture areas. After fast / slow shear wave separation and residual time difference correction, Figure 16 The energy of the PS1 wave in (c) is enhanced. Figure 16 In (d), the overall energy of the PS2 wave is weakened relative to the T component, indicating that the projection of the PS1 wave on the T component is effectively separated.

[0194] Figure 18 Focus on Figure 17 The imaging area within the blue box. Before shear wave splitting correction ( Figure 18In (a)-(b)), the phase axes of the PS1 wave (marked by the pink dashed line) and the PS2 wave (marked by the blue dashed line) both have partial projections onto the radial (R) and tangential (T) components, and the phase axis of the PS2 wave exhibits opposite polarities in the R and T components. After transverse wave splitting correction ( Figure 18 In (c)-(d), the phase axes of fast and slow shear waves become clearer and more continuous, the polarity reversal phenomenon is eliminated, and travel time differences can be extracted for crack density estimation. Figure 19 (a) and Figure 20 Figure (a) shows the variation trends of the travel time of the PS1 wave (green curve) and PS2 wave (blue curve) with Xline in groups W and L, respectively. In group W, the travel time curve fluctuates significantly overall, especially in the Xline region of 2050-2150, where the curve fluctuates violently and the difference between the travel times of PS1 and PS2 increases significantly, corresponding to the region of enhanced anisotropy. Figure 19 As shown in (b), the fracture density also exhibits two peaks in this Xline region; while in areas with smaller travel time differences, the fracture density is negative, indicating that the fracture development in this area is relatively low. Figure 20 In Figure (a), the travel time curves of PS1 and PS2 waves in group L show an overall upward trend, while the fracture density fluctuates periodically, with values ​​ranging from -2 to 2. Local peaks correspond to areas with high fracture development, while troughs indicate areas with lower fracture development.

[0195] Directional vertical fractures in strata induce significant azimuth anisotropy in the medium. The ray paths of converted waves incorporate both P-waves and S-waves, containing not only P-wave azimuth anisotropy information but, more importantly, shear wave splitting provides richer seismic responses for fracture parameter inversion, thus improving fracture identification accuracy and the reliability of quantitative characterization. When using 3D seismic data for pre-stack fracture prediction, dividing the data into azimuth sectors is a common strategy for extracting azimuth anisotropy information. Two main technical approaches exist: OVT domain processing and azimuth migration based on shot-receiver azimuth. OVT processing, by merging seismic traces with the same migration vector (including offset and azimuth), theoretically provides complete azimuth information for wide-azimuth seismic data. However, OVT domain processing faces numerous challenges in practical applications: ideal OVT domain data acquisition is costly and difficult to fully achieve; the irregular spatial distribution of shot and receiver points in conventional observation systems often leads to severely uneven OVT domain coverage, inconsistent gather quality, and the introduction of noise issues such as acquisition footprints and spatial aliasing, typically requiring complex regularization techniques for compensation. In contrast, the embodiment of this application employs a sub-azimuth processing flow based on shot and receiver point azimuth, which demonstrates better applicability and stability, especially for conventional 3D seismic data with highly uneven azimuth coverage. Based on the expected characteristics of geological targets and fracture development in the work area, the azimuth angles of shot and receiver points from 0° to 360° can be divided into several sector-shaped regions, thereby maximizing the utilization of sub-azimuth acquired data.

[0196] When the orientation of underground fractures changes with depth, accurate separation of fast and slow shear waves and analysis of the superposition effects of anisotropy are crucial. This application first analyzes the reflected amplitude energy and travel time characteristics of converted waves in different orientations to constrain the dominant orientation of fracture development from shallow to deep. Second, layer-by-layer Alford rotation is performed on the azimuth domain common imaging gather (ACIG) to achieve layer-by-layer separation of fast and slow shear waves. Then, considering that mid-to-far offsets are more sensitive to fracture response and that the angle has periodic characteristics, appropriate weighting factors are assigned to different offsets, and the principal orientation of fractures in the target layer is calculated using a weighted circular mean. Furthermore, to address the potential impact of background tectonic changes, this application introduces a Z-score normalization method based on a sliding time window to characterize the spatial variation trend of fracture density, highlighting local changes in anomalous splitting time differences, thereby enhancing the ability to identify fracture-rich zones. Compared with the traditional method of characterizing fracture density using the time difference of fast and slow shear waves, the Z-score method has stronger resistance to background interference, especially in areas with strong tectonic undulations.

[0197] This application proposes a fracture prediction method based on seismic data, particularly for OBN seismic data, demonstrating its potential for predicting fracture parameters in multi-layered fracture systems. The Alford rotation method, involving layer stripping, ensures the accuracy of fracture orientation prediction. After correcting for the residual time difference in the snail gather, the fast and slow SWW time differences are picked up more robustly on the stacked profiles. By comprehensively utilizing P-wave and SWW azimuth attributes as well as multi-azimuth migration gathers, a relatively robust estimation method is provided for the problem of fracture parameter estimation in complex strata.

[0198] This application presents a fracture prediction method based on seismic data. First, the azimuth anisotropy of PS waves and the travel time response characteristics of shear wave splitting in multi-layered fractured media are clarified. Based on this, a correction process for shear wave splitting effects is constructed, including: an azimuth-based segmentation strategy for finite azimuth coverage data, fast and slow shear wave separation based on azimuth domain co-imaging gathers, and a time difference correction method based on snail gathers. Regarding fracture parameter prediction, in pre-stack gathers, the principal fracture direction of the target layer is calculated using a weighted circular mean, combining the signal-to-noise ratio at different offsets and the periodic distribution characteristics of angle data. In post-stack profiles, a Z-score normalization method based on a sliding time window is introduced to quantitatively characterize the degree of fracture density development.

[0199] Meanwhile, the embodiments of this application effectively correct for the shear wave splitting effect. This separates the leaked PS1 wave in the T component, compensates for the imaging energy of the radial component, thereby improving the imaging quality of the PS wave profile and enabling stable prediction of fracture parameters at multiple fracture layers. Both model and actual data tests demonstrate that this method can effectively correct for the shear wave splitting effect, enhancing the reliability and resolution of fracture orientation and density identification.

[0200] It should be noted that the method in this embodiment can be executed by a single device, such as a computer or server. The method can also be applied in a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method in this embodiment, and the multiple devices will interact with each other to complete the method described.

[0201] It should be noted that the above description describes some embodiments of this application. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in a different order than that shown in the above embodiments and still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0202] Based on the same technical concept, corresponding to any of the above embodiments, this application also provides a crack prediction device based on seismic data.

[0203] refer to Figure 21 The crack prediction device based on seismic data includes:

[0204] The imaging processing module 201 is used to divide the seismic data into multiple sector regions according to the azimuth of the shot point and the detector, and to determine the azimuth domain common imaging gather of the converted wave in each sector region.

[0205] The wave field separation module 202 is used to perform fast and slow shear wave separation on the azimuth domain common imaging gather, and to determine the fast shear wave gather, the slow shear wave gather, and the angle between the fracture orientation and the main survey line direction of the geophone in each formation fracture.

[0206] The parameter inversion module 203 is used to determine the fracture development direction of the target formation based on the fast shear wave gather and the slow shear wave gather.

[0207] The time difference correction module 204 is used to convert the azimuth domain co-imaging gathers corresponding to the fast shear wave and the slow shear wave into snail gathers according to the azimuth and offset of the fan-shaped region, and to correct the time difference of each channel in the snail gathers.

[0208] Crack development module 205 is used to determine the density of the cracks based on the time difference, and to determine the crack development index based on the density.

[0209] For ease of description, the above devices are described in terms of function, divided into various modules. Of course, in implementing this application, the functions of each module can be implemented in one or more software and / or hardware.

[0210] The apparatus described above is used to implement the corresponding crack prediction method based on seismic data in any of the foregoing embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0211] Based on the same technical concept, corresponding to the methods of any of the above embodiments, this application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the crack prediction method based on seismic data as described in any of the above embodiments.

[0212] Figure 22This embodiment illustrates a more specific hardware structure of an electronic device, which may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.

[0213] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0214] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.

[0215] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.

[0216] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).

[0217] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.

[0218] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.

[0219] The electronic devices described above are used to implement the corresponding crack prediction methods based on seismic data in any of the foregoing embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0220] Based on the same technical concept, corresponding to the methods of any of the above embodiments, this application also provides a non-transitory computer-readable storage medium storing computer instructions for causing the computer to execute the crack prediction method based on seismic data as described in any of the above embodiments.

[0221] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0222] The computer instructions stored in the storage medium of the above embodiments are used to cause the computer to execute the crack prediction method based on seismic data as described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0223] Based on the same inventive concept, corresponding to the crack prediction method based on seismic data described in any of the above embodiments, this disclosure also provides a computer program product, which includes computer program instructions. In some embodiments, the computer program instructions can be executed by one or more processors of a computer to cause the computer and / or the processor to perform the crack prediction method based on seismic data. Corresponding to the execution entity for each step in each embodiment of the crack prediction method based on seismic data, the processor executing the corresponding step can belong to the corresponding execution entity.

[0224] The computer program product of the above embodiments is used to enable the computer and / or the processor to execute the crack prediction method based on seismic data as described in any of the above embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0225] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this application (including the claims) is limited to these examples; within the framework of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this application as described above, which are not provided in the details for the sake of brevity.

[0226] Additionally, to simplify the description and discussion, and to avoid obscuring the embodiments of this application, the well-known power / ground connections to integrated circuit (IC) chips and other components may or may not be shown in the provided drawings. Furthermore, the apparatus may be shown in block diagram form to avoid obscuring the embodiments of this application, and this also takes into account the fact that the details of the implementation of these block diagram apparatuses are highly dependent on the platform on which the embodiments of this application will be implemented (i.e., these details should be fully understood by those skilled in the art). While specific details (e.g., circuits) have been set forth to describe exemplary embodiments of this application, it will be apparent to those skilled in the art that the embodiments of this application can be implemented without these specific details or with variations thereof. Therefore, these descriptions should be considered illustrative rather than restrictive.

[0227] Although this application has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.

[0228] The embodiments of this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the embodiments of this application should be included within the protection scope of this application.

Claims

1. A crack prediction method based on seismic data, characterized in that, include: The seismic data is divided into multiple sector regions according to the azimuth of the shot point and the geophone, and the azimuth domain common imaging gather of the converted wave in each sector region is determined. The azimuth domain common imaging gather is subjected to fast and slow shear wave separation to determine the fast shear wave gather, the slow shear wave gather, and the angle between the fracture orientation and the main survey line direction of the geophone in each formation fracture. This includes: rotating the azimuth domain common imaging gather at different offsets using a scanning angle; in response to the scanning angle satisfying the maximum energy ratio of the fast and slow shear waves, the scanning angle is determined to be the relative angle between the dominant azimuth of the current formation fracture and the fracture orientation at the current offset; based on the relative angle, the angle between the fracture orientation and the main survey line direction of the geophone is determined. Determining the fracture development direction of the target formation based on the fast shear wave gather and the slow shear wave gather includes: the offset distance includes near offset distance and mid-to-far offset distance; determining the circular mean of the offset distance as the scanning angle; determining a first weight corresponding to the mid-to-far offset distance and a second weight corresponding to the near offset distance; wherein the first weight is greater than the second weight; determining the fracture development direction of the target formation based on the first weight, the second weight, and the circular mean; the circular mean of the offset distance is expressed as: ,in, , For sampling points within the time window, This indicates all scan angles within the time window; The azimuth domain co-imaging gathers corresponding to the fast and slow shear waves are converted into snail gathers according to the azimuth and offset of the fan-shaped region, and the time difference of each channel in the snail gathers is corrected. This includes: using the near offset of the snail gathers as a reference channel to generate a standard reference stacking profile model; using the remaining channels in the snail gathers as channels to be corrected; performing cross-correlation calculations between the channels to be corrected and the reference stacking profile model to determine the time shift corresponding to the maximum correlation; and determining the time difference of the channels to be corrected based on the time shift. Based on the time difference, the density of the cracks is determined, and the development index of the cracks is determined according to the density, including: determining the magnitude of the time difference between the fast shear wave and the slow shear wave based on the time difference; using a sliding time window Z-score, determining the local mean and standard deviation of the residuals of the time difference magnitude within the sliding time window; and determining the crack development index based on the time difference magnitude, the local mean, and the standard deviation.

2. The method according to claim 1, characterized in that, Performing fast and slow shear wave separation on the azimuth domain common imaging gather to determine the fast and slow shear wave gathers includes: Using the main measurement line direction and the connecting line direction of the detector as the X-axis and Y-axis of the propagation coordinate system, the propagation vector composed of the fast shear wave and the slow shear wave is determined, so that the propagation vector is transformed into the coordinate system composed of the polarization direction of the fast shear wave and the polarization direction of the slow shear wave. Determine the reference coordinate system in which the detector is located, and rotate the propagation vector composed of the fast shear wave and the slow shear wave into the reference coordinate system to determine the seismic records after multiple shear wave separations of different strata fractures in the reference coordinate system. The earthquake record is represented as follows: ， in, and These represent the seismic records of multiple shear waves received by the fracture in the i-th stratum in the X and Y components of the reference coordinate system, after separation. This represents the angle between the direction of the crack in the i-th layer and the direction of the main survey line. This represents the transpose of a rotation matrix. This represents the time delay of the fracture in the i-th formation. This represents the propagation vector received at the detector.

3. The method according to claim 2, characterized in that, The azimuth domain common imaging gather is subjected to fast and slow shear wave separation to determine the fast and slow shear wave gathers, as well as the angle between the fracture orientation in each formation fracture and the main survey line direction of the geophone. This also includes: The azimuth domain common imaging gathers for each formation fracture, from shallow to deep, are determined in the direction of the main survey line and the direction of the connecting line. The azimuth domain common imaging gathers in the main survey line direction and the connecting line direction are rotated to the formation fractures to determine the fast shear wave and slow shear wave corresponding to each formation fracture.

4. A crack prediction device based on seismic data, characterized in that, include: The imaging processing module is used to divide the seismic data into multiple sector regions according to the azimuth of the shot point and the detector, and to determine the azimuth domain common imaging gather of the converted wave in each sector region. The wavefield separation module is used to separate fast and slow shear waves in the azimuth domain common imaging gather, and to determine the angle between the fast shear wave gather, the slow shear wave gather, and the fracture orientation of each formation fracture and the main survey line direction of the geophone. This includes: rotating the azimuth domain common imaging gather at different offsets using a scanning angle; in response to the scanning angle satisfying the maximum energy ratio of the fast and slow shear waves, the scanning angle is determined to be the relative angle between the dominant azimuth of the current formation fracture and the fracture orientation of the current formation fracture at the current offset; and determining the angle between the fracture orientation and the main survey line direction of the geophone based on the relative angle. The parameter inversion module is used to determine the fracture development direction of the target formation based on the fast shear wave gather and the slow shear wave gather, including: the offset distance includes near offset distance and mid-to-far offset distance; the circular mean of the offset distance is determined as the scanning angle; a first weight corresponding to the mid-to-far offset distance and a second weight corresponding to the near offset distance are determined; wherein the first weight is greater than the second weight; the fracture development direction of the target formation is determined based on the first weight, the second weight, and the circular mean; the circular mean of the offset distance is expressed as: ,in, , For sampling points within the time window, This indicates all scan angles within the time window; The time difference correction module is used to convert the azimuth domain co-imaging gathers corresponding to the fast and slow shear waves into snail gathers according to the azimuth and offset of the fan-shaped region, and to correct the time difference of each channel in the snail gathers. This includes: using the near offset of the snail gathers as a reference channel to generate a standard reference stacking profile model; using the remaining channels in the snail gathers as channels to be corrected; performing cross-correlation calculations between the channels to be corrected and the reference stacking profile model to determine the time shift corresponding to the maximum correlation; and determining the time difference of the channels to be corrected based on the time shift. A crack development module is used to determine the density of cracks based on the time difference and to determine the crack development index based on the density. This includes: determining the time difference between the fast and slow shear waves based on the time difference; using a sliding time window Z-score to determine the local mean and standard deviation of the residuals of the time difference within the sliding time window; and determining the crack development index based on the time difference, the local mean, and the standard deviation.

Images