A method and device for predicting sandstone fractures based on VSP technology

Abstract

The application provides a sandstone fracture prediction method and device based on VSP technology. The method obtains three-dimensional three-component vertical seismic profile VSP data containing downgoing P-wave and upgoing converted S-wave around a CO2 injection well; separates the downgoing P-wave and upgoing converted S-wave from the VSP data; obtains a background P-wave velocity field of a well area based on the first arrival time of the downgoing P-wave by using a travel time tomography inversion method; and determines the main strike of the fracture based on the travel time difference between the fast S-wave and the slow S-wave in the upgoing converted S-wave at different azimuths. The method realizes high-resolution three-dimensional prediction of the fracture density and the main strike of the tight sandstone reservoir.

Description

Technical Field

[0001] This application relates to the field of CO2 geological storage technology, and more specifically, to a method and apparatus for predicting sandstone fractures based on VSP technology. Background Technology

[0002] CO2 geological storage is one of the key technologies for mitigating global warming, with deep saline aquifers and sandstone formations in depleted oil and gas reservoirs being important storage targets. Tight sandstone reservoirs generally exhibit low porosity and low permeability, and the development of natural fractures within them directly affects CO2 injection capacity, migration pathways, and long-term storage safety. Accurately predicting the spatial distribution, density, and main orientation of natural formation fractures is a crucial prerequisite for CO2 injection site selection, well placement, and injection scheme optimization.

[0003] Currently, the industry primarily relies on conventional 3D ground seismic P-wave exploration methods to predict fractures in tight sandstone. These methods include: indirectly inferring fracture zones by analyzing geometric properties such as coherence and curvature volumes from post-stack seismic data; or using pre-stack P-wave azimuth anisotropy inversion to estimate fracture orientation and density by observing differences in P-wave travel time or amplitude at different azimuth angles. In addition, some studies have attempted to utilize the principle of shear wave splitting, but conventional ground seismic methods struggle to effectively excite and receive pure shear waves. Summary of the Invention

[0004] The purpose of this application is to provide a sandstone fracture prediction method and device based on VSP technology, which uses the travel time information of downlink P-waves and uplink converted S-waves to achieve high-resolution three-dimensional prediction of fractures and orientation in tight sandstone reservoirs, effectively overcoming the limitation of conventional seismic methods in identifying micro-fractures.

[0005] Firstly, a method for predicting sandstone fractures based on VSP technology is provided, which may include: Acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well; Downward P-waves and upward converted S-waves are separated from the VSP data; based on the first arrival time of the downward P-waves, the background P-wave velocity field in the well-side region is obtained using the travel-time tomography inversion method. The main orientation of the crack is determined based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

[0006] In one possible implementation, three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves is acquired around the CO2 injection well, including: A three-dimensional, three-component VSP observation system consisting of multiple concentrically arranged geophones cascaded together is deployed around the CO2 injection well, and excitation source points are distributed at equal angles on each concentric circle to collect the VSP data.

[0007] In one possible implementation, wavefield separation is performed from VSP data to separate the downlink P-wave and the uplink converted S-wave, including: A combination of median filtering and polarization filtering was used to separate the downlink longitudinal wave and the uplink converted transverse wave from the VSP data.

[0008] In one possible implementation, based on the first arrival time of the downlink P-wave, the background P-wave velocity field in the well-side region is obtained using travel-time tomography inversion, including: The first arrival times of the downlink P-wave at different receiving depths were collected as observation travel time data; For each ray path between the excitation point and the receiver point, a set of linear or nonlinear equations is established between the observed travel time data and the longitudinal wave velocity within the spatial grid cell based on the ray path along the straight or curved ray. The system of equations is solved by the least squares method to obtain the P-wave velocity of each spatial grid cell, thus forming the background P-wave velocity field.

[0009] In one possible implementation, the main orientation of the crack is determined based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles, including: Based on the propagation time difference between fast and slow shear waves in the upward converted shear waves at different azimuth angles, a functional relationship between the crack density along the seismic wave ray path and the propagation time difference is established. Based on this functional relationship, the fracture density of each spatial grid cell in the well-side region is calculated, and a three-dimensional fracture density distribution volume is generated. Based on the three-dimensional crack density distribution and propagation time difference, the main orientation of the crack is determined.

[0010] In one possible implementation, based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles, a functional relationship between the crack density along the seismic wave ray path and the propagation time difference is established, including: The propagation time difference corresponding to different azimuth angles at each receiving point is determined as the integral of the product of the crack density and the shear wave splitting coefficient along the ray path. The shear wave splitting coefficient is calculated based on the background P-wave velocity field and a preset background shear wave velocity field; the ray path is determined by ray tracing with the background P-wave velocity field as input; the crack density is an unknown parameter within the spatial grid cell to be inverted, and this parameter is used as one of the integrands in the integral expression.

[0011] In one possible implementation, the main orientation of the cracks is determined based on the three-dimensional crack density distribution and propagation time difference, including: For each spatial grid cell in the well-side area, the propagation time difference between the fast and slow shear waves at different azimuth angles is extracted to form a discrete propagation time difference sequence. The discrete sequence is fitted with an ellipse or a sine curve to obtain the extreme azimuth angle of the propagation time difference as a function of the azimuth angle. The azimuth angle corresponding to the maximum propagation time difference is determined as the main direction of the crack at that spatial grid cell.

[0012] Secondly, a sandstone fracture prediction device based on VSP technology is provided, the device may include: The acquisition unit is used to acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well. The separation unit is used to separate the downlink P-wave and uplink converted S-wave from the VSP data. The acquisition unit is also used to obtain the background P-wave velocity field in the well-side area based on the first arrival time of the downflow P-wave using the travel-time tomography inversion method. The determination unit is used to determine the main orientation of the crack based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

[0013] Thirdly, an electronic device is provided, which includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; When a processor executes a program stored in memory, it implements any of the steps described in the first aspect above.

[0014] Fourthly, a computer-readable storage medium is provided, wherein a computer program is stored therein, and when executed by a processor, the computer program implements the steps of any of the methods described in the first aspect above.

[0015] This application provides a sandstone fracture prediction method and apparatus based on VSP technology. The method acquires three-dimensional, three-component vertical seismic profile (VSP) data containing downflow P-waves and upflow converted S-waves around a CO2 injection well. Downflow P-waves and upflow converted S-waves are separated from the VSP data through wavefield separation. Based on the first arrival time of the downflow P-waves, the background P-wave velocity field in the well-side region is obtained using travel-time tomography inversion. The main orientation of the fracture is determined based on the propagation time difference between fast and slow S-waves in the upflow converted S-waves at different azimuth angles. This method, by jointly utilizing the travel-time information of downflow P-waves and upflow converted S-waves, achieves high-resolution three-dimensional prediction of fracture density and main orientation in tight sandstone reservoirs, effectively overcoming the limitation of conventional seismic methods in identifying microfractures. Attached Figure Description

[0016] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used in the embodiments of this application will be briefly introduced below. It should be understood that the following drawings only show some embodiments of this application and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 A schematic flowchart illustrating a sandstone fracture prediction method based on VSP technology provided in this application embodiment; Figure 2 A schematic diagram of a sandstone fracture prediction device based on VSP technology provided in this application embodiment; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application. Unless otherwise defined, the technical or scientific terms used in this application should have the ordinary meaning understood by those skilled in the art. The words "first," "second," and similar terms used in this application do not indicate any order, quantity, or importance, but are only used to distinguish different components. The words "comprising" or "including," etc., mean that the element or object preceding the word covers the element or object listed after the word and its equivalents, but do not exclude other elements or objects. The words "connected," "coupled," or "connected," etc., are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. "Up," "down," "left," "right," etc., 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.

[0019] For ease of understanding, the terms used in the embodiments of this application are explained below: 1. Downward P-wave: A P-wave that propagates directly downwards from a ground source without reflection to the geophone in the well. Its propagation direction is consistent with the direction of particle vibration. In this application, the first arrival time is used for travel-time tomography inversion to obtain the background P-wave velocity field.

[0020] 2. First arrival time: The travel time of a seismic wave from the excitation point to a receiving point (such as a downhole geophone) for the first time. In VSP data, the downlink P-wave usually arrives first, and its first arrival time is the moment corresponding to the first obvious jump or amplitude peak in the seismic record.

[0021] 3. Upward-converted S-wave: During the propagation of a seismic wave, when a P-wave encounters an interface or crack, it undergoes a wave mode conversion, transforming into a S-wave (transverse wave) and propagating upwards to reach the detector. Its propagation direction is upward, and the direction of particle vibration is perpendicular to the propagation direction.

[0022] 4. Fast and Slow Shear Waves: When a shear wave passes through an anisotropic medium (such as a reservoir containing oriented fractures), it splits into two waves with mutually perpendicular polarization directions and different propagation speeds. The faster wave is called a fast shear wave, with its polarization direction parallel to the fracture strike; the slower wave is called a slow shear wave, with its polarization direction perpendicular to the fracture strike. In this application, the propagation time difference (splitting time difference) between the two is directly related to the fracture density, and the azimuth of the largest time difference indicates the main fracture strike.

[0023] 5. Travel-time tomography inversion: A geophysical inversion method. By utilizing the travel time (travel time data) of seismic waves from the excitation point to the receiving point, and establishing a mathematical relationship between travel time and the velocity of the subsurface medium, the velocity distribution at various points in space can be deduced.

[0024] 6. Observation travel time data: The actual propagation time of seismic waves (here, downward P-waves) from the ground excitation point to the geophone at a specific depth in the well, obtained from the actual raw seismic records by manual or automatic algorithms.

[0025] 7. Background P-wave velocity field: The distribution of P-wave velocity in three-dimensional space of the tight sandstone matrix and its pore fluids, without considering the influence of fractures or with zero fracture density. In this application, it serves as the background model for subsequent fracture density inversion and is used to calculate ray paths and shear wave splitting coefficients.

[0026] 8. Shear wave splitting coefficient: A physical quantity characterizing the contribution of anisotropic intensity in the medium to the shear wave splitting time difference. Its magnitude depends on the background shear wave velocity and the flexibility characteristics of the crack.

[0027] 9. Main direction of cracks: The preferred direction of crack extension on the horizontal projection plane (usually expressed as azimuth, 0~180°). For a group of parallel or nearly parallel cracks, the main direction is their average direction.

[0028] 10. Polarization Filtering: A filtering method that uses the polarization characteristics of seismic wave particle motion (e.g., the linear polarization direction of P-waves is consistent with the propagation direction) to separate different wave types. In this application: the descending P-wave and the ascending converted S-wave are separated based on the difference between the polarization angle and the theoretical polarization angle of the P-wave.

[0029] This application focuses on fracture prediction in tight sandstone reservoirs where CO2 injection wells are located. First, a three-dimensional, three-component vertical seismic profile (VSP) observation system is deployed around the CO2 injection well. The system consists of multiple cascaded geophones arranged concentrically, each a three-component geophone capable of simultaneously recording the vertical component and two orthogonal horizontal components of the seismic wavefield. Excitation source points (such as explosive or controlled seismic sources) are distributed at equal angles along each concentric circle to obtain three-dimensional, three-component VSP data containing downlink P-waves and uplink converted S-waves.

[0030] Before deploying the observation system, the following steps are also included: Based on well logging data (such as sonic transit time logging, density logging, and natural gamma logging) and prior geological information (such as the direction of regional tectonic stress field, fault distribution, and core fracture observation statistics) in the area where the CO2 injection well is located (tight sandstone reservoir section), the dominant orientation and density range of fractures are estimated. Based on this, the number of concentric circles, the number of detectors per circle, and the azimuth interval of the excitation source points in the three-dimensional three-component VSP observation system are designed. Specifically, areas with higher estimated fracture density correspond to more concentric circles and smaller detector spacing. For example, if logging and core data indicate that the fracture density is expected to reach 0.05 to 0.10 (unit: number of fractures per meter) within a 500m radius around the well, and the dominant azimuth is 60° east of north, then three concentric circles are set with radii of 200m, 400m, and 600m respectively. Within a ±30° interval in the direction perpendicular to the dominant azimuth (i.e., 150° east of north), the azimuth interval of the excitation source point in each layer is set to 15°, while the interval in the interval parallel to the dominant azimuth is set to 30°, thus achieving non-uniform densification.

[0031] This design significantly improves the sampling accuracy of extreme values ​​of transverse wave splitting time difference as a function of azimuth angle, enhances the accuracy of determining the main direction of the crack, and reduces the cost of redundant data acquisition in insensitive azimuths.

[0032] The preferred embodiments of this application are described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit this application. Furthermore, the embodiments and features in the embodiments of this application can be combined with each other without conflict.

[0033] Figure 1 This is a schematic flowchart illustrating a sandstone fracture prediction method based on VSP technology, provided as an embodiment of this application. Figure 1 As shown, the method may include: Step S110: Obtain three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well.

[0034] In practice, a three-dimensional three-component VSP observation system consisting of multiple concentrically arranged detectors is deployed around the CO2 injection well, and excitation source points are distributed at equal angles on each concentric circle to collect the VSP data.

[0035] A typical implementation is as follows: Three concentric circles are set on the ground with radii R1=200m, R2=400m, and R3=600m, centered on the injection wellhead. Twenty-four geophone cascade points are evenly distributed on each concentric circle (the actual number of cascades can be adjusted according to cost). Each geophone cascade includes three seismic geophones in orthogonal directions (vertical component Z, radial component R, and tangential component T). Excitation source points are arranged at equal angular intervals (e.g., 15°) along each concentric circle, totaling 24 source points. After each source point is excited, geophones at all depths downhole (e.g., depth range 500m~2500m, spaced 10m apart) simultaneously record seismic waveforms, with a sampling rate of 1ms and a recording duration of 2s. The acquired raw data are three-component seismic gathers for each source-receiver pair.

[0036] Step S120: Separate the downlink longitudinal wave and the uplink converted transverse wave from the VSP data.

[0037] In practice, a combination of median filtering and polarization filtering is used to separate the downlink longitudinal wave and the uplink converted transverse wave from the VSP data, and to separate the fast transverse wave and the slow transverse wave from the uplink converted transverse wave.

[0038] In a preferred embodiment, the three-component records in the VSP data are first rotated to a radial-tangential-vertical coordinate system. Then, a mean-mode filter is applied in the radial-vertical plane to remove horizontal in-phase axes, enhancing the up-wave energy. Next, polarization angle filtering is applied to the radial component, separating the down-current P-wave and the up-converted S-wave based on the difference between the polarization angle and the theoretical polarization angle of the P-wave. For the separated up-converted S-wave, covariance matrix eigenvalue decomposition is further applied to the tangential and radial components to obtain the direction and energy of the fast and slow S-waves. This process specifically includes: First, the three-component records in the VSP data are rotated. The original data is typically recorded in a geographic coordinate system (East, North, Vertical), and needs to be rotated to a radial-tangential-vertical coordinate system (RTZ). The rotation angle is determined by the azimuth angle between the source point and the detector: the horizontal component is projected onto the line connecting the source point and the detector to obtain the radial component R, and the vertical component T is obtained. The rotation formula is:

[0039] Where R is the radial component, T is the tangential component, E is the east component, N is the north component, and θ is the azimuth angle of the source point relative to the detector (calculated clockwise from true north). E, N, and Z are the east, north, and vertical components, respectively, which can be obtained directly from the original three-component record; θ can be calculated from the geographical coordinates of the source point and the detector.

[0040] Then, apply mean filtering in the radial-vertical plane (RZ plane). Select a time window length (e.g., 5 sampling points) and slide it along the time direction, taking the median of the R and Z components respectively. This can effectively filter out horizontal in-phase axes (e.g., shear waves, surface waves) and enhance the energy of rising waves (especially converted shear waves).

[0041] Then, polarization angle filtering is applied to the radial component R. The instantaneous polarization angle at each sampling point is calculated:

[0042] in, Let Z(t) be the instantaneous polarization angle at time t, Z(t) be the vertical component seismic record, and R(t) be the radial component seismic record. Z(t) and R(t) can be directly read from the median-filtered data. The theoretical polarization angle of the P-wave is close to 45° in the radial to vertical plane (when Poisson's ratio is 0.25), while the polarization angle of the S-wave is close to 0° or 90°. Set an angle threshold range (e.g., 25° to 65°), retaining signals within this range as downlink P-waves, and signals outside the range as uplink converted S-wave candidates.

[0043] For the separated upward-converting shear waves, covariance matrix eigendecomposition is further applied to the tangential component T and the radial component R to obtain the direction and energy of the fast and slow shear waves. Specifically, a covariance matrix is ​​constructed for the three-component data within each time window:

[0044] Where R is the radial component, T is the tangential component, and Z is the perpendicular component. · This represents the average within the time window. R, T, and Z can be obtained from the data after the coordinate rotation and filtering described above. Regarding the covariance matrix... Eigenvalue decomposition is performed, and the eigenvector corresponding to the largest eigenvalue indicates the principal direction of particle motion. For shear waves, the direction of particle motion is approximately perpendicular to the propagation direction. By rotating, the energy of one component is concentrated in the radial direction R (fast shear wave), and the other component is concentrated in the tangential direction T (slow shear wave), thus achieving fast and slow shear wave separation. This method avoids the problems of high computational cost and low directional resolution of conventional rotating scans, and enables the separation accuracy of fast and slow shear waves to meet the detection requirements of weak anisotropy in dense sandstone. The joint filtering strategy combines the separation advantages of median filtering in the apparent velocity domain and polarization filtering in the polarization domain, effectively overcoming the limitations of single filtering methods in complex wavefield environments and improving the fidelity of wavefield separation.

[0045] Step S130: Based on the first arrival time of the downlink P-wave, the background P-wave velocity field in the well-side area is obtained using the travel-time tomography inversion method.

[0046] In practice, the arrival times of the downlink P-waves at different receiving depths are picked up as observation travel time data. For each ray path between the excitation point and the receiving point, a set of linear or nonlinear equations between the observation travel time data and the P-wave velocity in the spatial grid cell is established based on the ray path along the straight or curved ray. The P-wave velocity of each spatial grid cell is obtained by solving the set of equations using the least squares method, thus forming the background P-wave velocity field.

[0047] Specifically, the area around the well to be predicted is discretized into three-dimensional spatial grid cells, where the grid size is, for example, taken as... The P-wave velocity within each grid cell is considered constant. For the first The excitation point and the first Ray path between receiving points .

[0048] (1) Straight ray hypothesis Under the straight-ray assumption, seismic waves propagate along a straight line from the excitation point to the receiver point. The data is equal to the path length of the ray through each grid divided by the sum of the P-wave velocities of that grid cell:

[0049] in, The data represents the travel time from the i-th excitation point to the j-th receiving point; K is the total number of three-dimensional mesh cells. The length of the line segment through which the straight ray between the i-th excitation point and the j-th receiver point passes through the k-th grid cell is obtained by analytical calculation based on the spatial coordinates of the excitation point and the receiver point and the grid geometric boundary. Let be the P-wave velocity in the k-th grid cell, and be the unknown parameter to be solved. This equation is about... The linear equations can be solved using the damped least squares method.

[0050] (2) Bending ray hypothesis Under the assumption of a curved ray, the ray path The velocity field itself follows the longitudinal wave. The distribution varies according to Fermat's principle (minimum travel time path). In this case, the observed travel time data is represented as a slow field along the actual ray path. Line integral:

[0051] in, To observe travel time data; For dependent on the longitudinal wave velocity field The curved ray path from excitation point i to receiving point j is obtained by solving the ray tracing equation based on the velocity field model. Let r be the slow-degree field, representing the propagation time per unit length at spatial location r, and ; Let be the longitudinal wave velocity at spatial position r; dl be the arc length infinitesimal element along the ray path; and r be the three-dimensional spatial position vector. Due to the ray path With velocity field The equations are mutually coupled and form a system of nonlinear integral equations, which need to be solved by iterative algorithms such as the nonlinear conjugate gradient method.

[0052] In a preferred embodiment, considering the dramatic lateral variations in the velocity field of tight sandstone, a single inversion method is prone to getting trapped in local extrema. This embodiment employs a multi-step, multi-scale tomographic inversion strategy: The first step utilizes only the initial arrival time of the downlink P-wave at zero offset (very small distance between the source and the wellhead) or small offset (e.g., less than 200m). In the inversion method, the straight ray assumption is used for ray path tracing, and the observed travel time data is compared with the propagation time. A system of linear equations was established and solved using the damped least squares method. This step aims to obtain the low wavenumber components (i.e., macroscopic velocity background) of the background P-wave velocity field in the well-side region, reflecting the overall trend of formation velocity variation with depth and large-scale structural features, thus providing a stable initial model for subsequent fine inversion.

[0053] The second step, based on the macroscopic velocity model obtained in the first step, further introduces the first arrival time of the downflow P-wave at mid-to-long-range offsets (e.g., 200m~600m). In terms of inversion methods, the curved ray assumption is used for ray path tracing, and the observed travel time data is expressed as an integral along the slow-degree line of the curved ray, forming a nonlinear equation system, which is then solved iteratively using the nonlinear conjugate gradient method. This step uses the low-wavenumber velocity model from the first step as the initial iterative model to invert and obtain the high-wavenumber details (i.e., intra-layer heterogeneity) of the background P-wave velocity field. This allows for the characterization of lateral velocity fluctuations within tight sandstone reservoirs caused by lithological variations, local tight zones, or fracture development.

[0054] Furthermore, spatial smoothing constraints (such as the Laplace smoothing operator) and logging P-wave velocity curves can be added as soft constraints in each inversion step to prevent overfitting. This multi-step, multi-scale strategy solves the problem that single-tomographic inversion is prone to getting trapped in local extrema when the background velocity field changes drastically laterally in tight sandstone, providing a more accurate ray path and background model for subsequent fracture density inversion.

[0055] Step S140: Determine the main orientation of the crack based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

[0056] First, based on the propagation time difference between fast and slow shear waves in the ascending converted shear waves at different azimuth angles, a functional relationship between crack density and propagation time difference along the seismic wave ray path is established. Based on this functional relationship, the crack density of each spatial grid cell in the well-side region is calculated through inversion, generating a three-dimensional crack density distribution. Specifically, the well-side region to be predicted is discretized into three-dimensional spatial grid cells.

[0057] Specifically, the area around the well to be predicted is first discretized into three-dimensional spatial grid cells. Let the propagation time difference between the fast and slow shear waves of the upward-converted shear wave emitted from the source point at the j-th receiving point at the n-th azimuth angle be denoted as the shear wave splitting time difference. The integral expression for this time difference along the ray path is constructed as follows: The propagation time difference corresponding to different azimuth angles at each receiving point is determined as the integral of the product of the crack density and the shear wave splitting coefficient along the ray path; wherein, the shear wave splitting coefficient is calculated based on the background P-wave velocity field and the preset background shear wave velocity field; the ray path is determined by ray tracing with the background P-wave velocity field as input; the crack density is an unknown parameter in the spatial grid cell to be inverted, which is used as one of the integrands in the integral expression.

[0058] For example, suppose that at the j-th receiving point, an upward-converting shear wave is emitted from the source point at the n-th azimuth angle, and the propagation time difference (i.e., shear wave splitting time difference) between its fast and slow shear waves is denoted as... The integral expression for this time difference along the ray path is: ; in: This is the ray path originating from the receiving point and following the upward converted shear wave propagation path (from the conversion point to the receiving point). This path passes through the background P-wave velocity field and the preset background shear wave velocity field (determined by empirical formulas). (or well logging shear wave data) is used as input for ray tracing determination. For spatial location The crack density at a given location (an unknown parameter to be inverted) is defined as the total area of ​​cracks per unit volume or the number of cracks per unit length, as needed. This is the shear wave splitting coefficient, whose value depends on the background velocity and the crack's compliance characteristics. It can be calculated based on the background P-wave velocity field and the background S-wave velocity field, for example... ,in Background shear wave velocity, and denoted as the normal compliance and tangential compliance of the fracture (which can be estimated using a rock physics model), and dl is the arc length infinitesimal element along the ray path.

[0059] To solve the above integral equation, it is discretized into a system of linear equations. Suppose there are M receiver-source pairs (i.e., M propagation time difference observations at different azimuth angles) and N grid cells to be inverted. Then, for the p-th observation:

[0060] in, Let be the propagation time difference vector between the fast and slow shear waves at different azimuth angles at the p-th receiving point. Let be the path length of the ray from the p-th receiving point through the q-th grid cell. Let be the shear wave splitting coefficient within the q-th grid cell. Let be... Then the system of equations is Where A is an M×N matrix, and the element in the p-th row and q-th column is... ; Let be the crack density to be determined, which is an N×1 column vector storing the crack density of all mesh elements to be inverted. It is an M×1 column vector that stores the propagation time difference of fast and slow shear waves measured at all M receiving points with different azimuth angles.

[0061] The system of equations can be solved using either algebraic reconstruction techniques or the conjugate gradient method. Taking algebraic reconstruction techniques as an example, the iterative format is as follows: ; in, and These are the crack densities of the q-th mesh element at the (k+1)-th and k-th iterations, respectively. The coefficient matrix elements are denoted by μ, which is the relaxation factor (0 < μ < 2). The iteration continues until the residual is less than a preset threshold. Finally, the crack density of each spatial grid cell is obtained, forming a three-dimensional crack density distribution.

[0062] As an extension of the preferred logic, in actual calculations, the ray path is related to the background velocity field, which has already been obtained. However, changes in crack density can cause local velocity variations, thus affecting the ray path. Therefore, an iterative re-forward modeling strategy can be adopted: after updating the crack density in each inversion iteration, the shear wave splitting coefficient α and the ray path L are updated according to the new crack density, matrix A is recalculated, and the inversion is performed again until the time difference residuals converge. This strategy ensures the consistency between the inversion results and the physical process of seismic wave propagation, preventing the accumulation of errors from fixed-path inversion.

[0063] Subsequently, the main orientation of the cracks was determined based on the three-dimensional crack density distribution and propagation time difference.

[0064] In practice, for each spatial grid cell in the well-side area, the propagation time difference between the fast and slow shear waves at different azimuth angles is extracted, forming a discrete azimuth-propagation time difference sequence. Elliptic or sine curve fitting is performed on the discrete sequence to obtain the extreme azimuth angles of the propagation time difference as a function of azimuth angle. The azimuth angle corresponding to the maximum propagation time difference is determined as the main fracture orientation at that spatial grid cell. Three-dimensional spatial interpolation is then performed on the main fracture orientations of all spatial grid cells to output a three-dimensional fracture orientation volume.

[0065] Specifically: For the j-th grid cell, the different azimuth angles corresponding to that grid cell can be obtained from the inversion process. Propagation time difference (k=1,…,K) The time difference data itself is the integral value along the ray path. To assign it to the grid cells, a ray density-weighted average or the time difference corresponding to the ray passing through the grid is used. A more robust approach is to utilize the established ray path and time difference observations to simultaneously obtain the contribution weight of each grid cell to the time difference in each azimuth angle when solving for the crack density, thereby back-projecting the time difference onto the grid cells.

[0066] After obtaining the azimuth-time difference sequence of the grid, ellipse fitting is used: assuming ,in, Azimuth The corresponding propagation time difference; A is the average value, B is the amplitude. This represents the extreme azimuth angle. The parameters A, B, and C are fitted using least squares. ,but This is the azimuth angle at which the propagation time difference is maximum, which is also the main direction of the crack. A sine curve can also be used for fitting. Where a, b, and c are fitting coefficients, and the extreme value azimuth angle is... The azimuth angle corresponding to the maximum propagation time difference is determined as the main crack orientation at that spatial grid cell.

[0067] Furthermore, the above fitting is performed on all grid cells to obtain the principal orientation value (range 0~180°) for each grid cell. Since the grid is discrete, there may be abrupt changes between adjacent grid cells. Therefore, inverse distance weighted interpolation or Kriging interpolation is used to interpolate the orientation values ​​at grid points to the entire space, forming a smooth 3D fracture orientation volume. This orientation volume, together with the fracture density volume, constitutes the final result of fracture prediction and can be used for CO2 injection scheme optimization and leakage risk assessment. Inverse distance weighted interpolation, based on the distance decay principle, assigns higher weights to nearby grid points, is computationally simple, and can better preserve local features; Kriging interpolation considers the spatial autocorrelation of spatial variables, providing unbiased optimal estimates and is suitable for scenarios where the fracture orientation has spatial continuity. The discrete grid orientation values ​​are not used directly because they are significantly affected by observation noise and inversion errors, and abrupt changes between adjacent grid cells may occur, failing to truly reflect the gradual spatial changes in fracture orientation. Interpolation smoothing can effectively suppress high-frequency noise and generate a continuously transitioning three-dimensional orientation volume, which, together with the crack density volume, constitutes the final result of crack prediction, thus providing more intuitive guidance for the optimization of CO2 injection schemes and leakage risk assessment.

[0068] Corresponding to the above method, this application also provides a sandstone fracture prediction device based on VSP technology, such as... Figure 2 As shown, the device includes: Acquisition unit 210 is used to acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well. Separation unit 220 is used to separate the downlink longitudinal wave and the uplink converted transverse wave from the VSP data; The acquisition unit 210 is also used to obtain the background P-wave velocity field in the well-side area based on the first arrival time of the downflow P-wave using the travel-time tomography inversion method. The determination unit 230 is used to determine the main orientation of the crack based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

[0069] The functions of each functional unit of the sandstone fracture prediction device of VSP provided in the above embodiments of this application can be realized through the above methods and steps. Therefore, the specific working process and beneficial effects of each unit in the sandstone fracture prediction device of VSP provided in the embodiments of this application will not be repeated here.

[0070] This application also provides an electronic device, such as... Figure 3 As shown, it includes a processor 310, a communication interface 320, a memory 330, and a communication bus 340, wherein the processor 310, the communication interface 320, and the memory 330 communicate with each other through the communication bus 340.

[0071] Memory 330 is used to store computer programs; When the processor 310 executes the program stored in the memory 330, it performs the following steps: Acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well; Downward P-waves and upward converted S-waves are separated from the VSP data; based on the first arrival time of the downward P-waves, the background P-wave velocity field in the well-side region is obtained using the travel-time tomography inversion method. The main orientation of the crack is determined based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

[0072] The communication bus mentioned above can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. This communication bus can be divided into address bus, data bus, control bus, etc. For ease of illustration, only one thick line is used to represent it in the diagram, but this does not mean that there is only one bus or one type of bus.

[0073] The communication interface is used for communication between the aforementioned electronic devices and other devices.

[0074] The memory may include random access memory (RAM) or non-volatile memory (NVM), such as at least one disk storage device. Optionally, the memory may also be at least one storage device located remotely from the aforementioned processor.

[0075] The processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0076] The implementation methods and beneficial effects of the various components of the electronic device in the above embodiments for solving the problem can be found in [reference needed]. Figure 1 The steps in the illustrated embodiments are used to implement the electronic device. Therefore, the specific working process and beneficial effects of the electronic device provided in this application will not be repeated here.

[0077] In another embodiment provided in this application, a computer-readable storage medium is also provided, which stores instructions that, when executed on a computer, cause the computer to perform the sandstone fracture prediction method of VSP as described in any of the above embodiments.

[0078] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to execute the sandstone fracture prediction method of VSP as described in any of the above embodiments.

[0079] Those skilled in the art will understand that the embodiments in this application can be provided as methods, systems, or computer program products. Therefore, the embodiments in this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the embodiments in this application can take the form of a computer program product implemented 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.

[0080] This application describes embodiments of methods, apparatus (systems), and computer program products according to embodiments of this application with reference to flowchart illustrations and / or block diagrams. 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. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

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

[0082] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps 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 steps of the function specified in one or more boxes.

[0083] Although preferred embodiments have been described in this application, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of the embodiments of this application.

[0084] Obviously, those skilled in the art can make various modifications and variations to the embodiments of this application without departing from the spirit and scope of the embodiments of this application. Therefore, if these modifications and variations to the embodiments of this application fall within the scope of the claims in this application and their equivalents, then this application also intends to include these modifications and variations.

Claims

1. A method for predicting sandstone fractures based on VSP technology, characterized in that, The method includes: Acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well; Wave field separation from VSP data reveals downlink longitudinal waves and uplink converted transverse waves; Based on the initial arrival time of the downflow P-wave, the background P-wave velocity field in the well-side region is obtained using the travel-time tomography inversion method. The main orientation of the crack is determined based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

2. The method as described in claim 1, characterized in that, Acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well, including: A three-dimensional, three-component VSP observation system consisting of multiple concentrically arranged geophones cascaded together is deployed around the CO2 injection well, and excitation source points are distributed at equal angles on each concentric circle to collect the VSP data.

3. The method as described in claim 1, characterized in that, Wavefield separation from VSP data reveals descending P-waves and ascending converted S-waves, including: A combination of median filtering and polarization filtering was used to separate the downlink longitudinal wave and the uplink converted transverse wave from the VSP data.

4. The method as described in claim 1, characterized in that, Based on the initial arrival time of the downflow P-wave, the background P-wave velocity field in the well-side region is obtained using the travel-time tomography inversion method, including: The first arrival times of the downlink P-wave at different receiving depths were collected as observation travel time data; For each ray path between the excitation point and the receiver point, a set of linear or nonlinear equations is established between the observed travel time data and the longitudinal wave velocity within the spatial grid cell based on the ray path along the straight or curved ray. The system of equations is solved by the least squares method to obtain the P-wave velocity of each spatial grid cell, thus forming the background P-wave velocity field.

5. The method as described in claim 1, characterized in that, Based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles, the main orientation of the crack is determined, including: Based on the propagation time difference between fast and slow shear waves in the upward converted shear waves at different azimuth angles, a functional relationship between the crack density along the seismic wave ray path and the propagation time difference is established. Based on this functional relationship, the fracture density of each spatial grid cell in the well-side region is calculated, and a three-dimensional fracture density distribution volume is generated. Based on the three-dimensional crack density distribution and propagation time difference, the main orientation of the crack is determined.

6. The method as described in claim 5, characterized in that, Based on the propagation time difference between fast and slow shear waves in the upward converted shear waves at different azimuth angles, a functional relationship between the crack density along the seismic wave ray path and the propagation time difference is established, including: The propagation time difference corresponding to different azimuth angles at each receiving point is determined as the integral of the product of the crack density and the shear wave splitting coefficient along the ray path. The shear wave splitting coefficient is calculated based on the background P-wave velocity field and a preset background shear wave velocity field; the ray path is determined by ray tracing with the background P-wave velocity field as input; the crack density is an unknown parameter within the spatial grid cell to be inverted, and this parameter is used as one of the integrands in the integral expression.

7. The method as described in claim 5, characterized in that, Based on the three-dimensional crack density distribution and propagation time difference, the main orientation of the cracks is determined, including: For each spatial grid cell in the well-side area, the propagation time difference between the fast and slow shear waves at different azimuth angles is extracted to form a discrete propagation time difference sequence. The discrete sequence is fitted with an ellipse or a sine curve to obtain the extreme azimuth angle of the propagation time difference as a function of the azimuth angle. The azimuth angle corresponding to the maximum propagation time difference is determined as the main direction of the crack at that spatial grid cell.

8. A sandstone fracture prediction device based on VSP technology, characterized in that, The device includes: The acquisition unit is used to acquire three-dimensional three-component vertical seismic profile (VSP) data containing downlink P-waves and uplink converted S-waves around the CO2 injection well. The separation unit is used to separate the downlink P-wave and uplink converted S-wave from the VSP data. The acquisition unit is also used to obtain the background P-wave velocity field in the well-side area based on the first arrival time of the downflow P-wave using the travel-time tomography inversion method. The determination unit is used to determine the main orientation of the crack based on the propagation time difference between the fast and slow shear waves in the upward converted shear waves at different azimuth angles.

9. An electronic device, characterized in that, The electronic device includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; A processor, when executing a program stored in memory, implements the method of any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the method described in any one of claims 1-7.

Images