Axial strain rate field simulation method of distributed optical fiber DAS-VSP based on strain tensor projection under non-vertical well layout condition

Abstract

The invention discloses an axial strain rate field simulation method of a distributed optical fiber DAS-VSP based on strain tensor projection under a non-vertical well layout condition, and belongs to the field of exploration. Aiming at the problem that a traditional vertical well simulation method is difficult to accurately simulate wave field propagation characteristics in a curved well trajectory, high-precision wave field simulation under a non-vertical well type condition is realized by establishing a parameterized well trajectory model based on a curvilinear function and combining a first-order speed-stress elastic wave equation and a finite difference algorithm. A tangential direction projection algorithm is introduced to convert a velocity field into optical fiber axial strain rate response, and a dynamic scale distance segmentation method is adopted to truly simulate the spatial sampling characteristics of an actual DAS system. By constructing PML absorption boundary conditions and a stability control mechanism, the numerical stability of highly-deviated well section wave field simulation is ensured. According to the method, a reliable numerical experiment method is provided for DAS-VSP data interpretation and forward modeling under the non-straight well type condition, and important technical support is provided for non-straight well optical fiber earthquake monitoring in oil-gas exploration.

Description

Technical Field

[0001] This invention relates to the field of exploration technology, and in particular to a method for simulating the elastic field response of a distributed optical fiber DAS-VSP in a well under non-vertical well layout conditions. Background Technology

[0002] In the detailed exploration of deep and complex targets in my country, VSP (Vertical Seismic Profile) utilizes sensors embedded deep in underground wells to acquire seismic signals, achieving high-resolution imaging of subsurface structures at a relatively low cost. VSP is receiving increasing attention due to its rich wavefield information and high-quality data performance. Fiber optic distributed acoustic sensing (DAS), with its lower cost, denser spatial sampling at the meter level, wider bandwidth, extensive receiver array coverage over tens of kilometers, and the potential for permanent installation / monitoring without interfering with other activities within the well, enables the rapid application of DAS technology to compensate for or replace some traditional seismic detector measurements.

[0003] DAS-VSP technology refers to using optical fibers to replace in-well seismic detector arrays for measuring seismic signals. Seismic detectors record the particle vibration velocity caused by seismic waves, while the DAS demodulation unit measures the phase change of backscattered laser light along fiber impurities, expressed as the average axial strain rate of the fiber along the gauge length. Current research on seismic response simulation methods for DAS-VSP mainly focuses on vertical wells, and a relatively systematic theoretical framework and numerical implementation methods for wavefield simulation have been established. However, with the deepening of oil and gas exploration into complex structures and unconventional reservoirs, the application of non-vertical wells in actual exploration is becoming increasingly widespread. Existing simulation systems under vertical well conditions struggle to accurately characterize the impact of wellbore bending on the fiber strain rate response, and existing methods lack a complete understanding of the interaction mechanism between elastic wave propagation and fiber sensing in curved well sections. This limitation restricts the forward modeling of DAS-VSP in complex exploration scenarios, while research on fiber optic seismic response simulation methods under non-vertical well deployment conditions remains relatively limited. There is an urgent need to develop a more comprehensive seismic forward modeling method for in-well distributed fiber optic DAS-VSP under non-vertical well deployment conditions. Summary of the Invention

[0004] This invention addresses the problem of simulating distributed fiber optic DAS-VSP data in non-vertical well environments. It proposes a simulation method for the axial strain rate field of distributed fiber optic DAS-VSP based on strain tensor projection under non-vertical well deployment conditions, so as to achieve accurate simulation of wave field propagation characteristics in curved well trajectories and overcome the shortcomings of existing DAS-VSP wave field response simulation methods that are limited to vertical well models.

[0005] The technical solution adopted in this invention is a method for simulating the axial strain rate field of distributed optical fiber DAS-VSP based on strain tensor projection under non-vertical well deployment conditions, which includes the following steps:

[0006] Step 1: Setting up and storing non-vertical well grid coordinates;

[0007] Based on the spatial range of the velocity model and the exploration target, a continuous mathematical function describing the non-vertical well trajectory is established to obtain a well trajectory parameterization model to obtain a continuous well trajectory. The well trajectory parameterization model uses the spatial coordinates of the trajectory start and end points as the boundary conditions of the continuous well trajectory, and adjusts the degree of well curvature and morphological characteristics through control functions.

[0008] The continuous well trajectory is converted into a discrete form and the grid coordinates are calculated: uniform sampling is performed along the parametric curve to generate a dense set of location points; then the physical coordinates corresponding to each location point in the location point set are mapped to a predefined computational grid system, and the coordinate values ​​of each location point are rounded (e.g., rounded to the nearest integer) to convert them into integer grid indices, thus obtaining non-vertical well grid points.

[0009] Step 2: Calculate and store the tangent parameters at non-vertical well grid points;

[0010] Traverse each non-vertical well grid point and calculate the tangent direction parameters at each grid point;

[0011] The tangent direction parameters at each grid point are converted into cosine and sine function values ​​to construct a direction parameter matrix. Each element of the direction parameter matrix is ​​used to characterize the local well trajectory direction characteristics of a non-vertical well grid point.

[0012] Step 3: Perform numerical simulation of elastic wave field based on the direction parameter matrix, and calculate the axial strain rate of each non-vertical well grid point, that is, calculate the strain rate component of each grid point along the well trajectory axis.

[0013] A single non-vertical well grid point is used to characterize a geophone location. The axial strain rate time series of the geophone location is used to construct discrete geophone strain rate time history data in the form of gathers to complete the forward modeling output data.

[0014] Step 4: Perform gauge-length averaging on the strain rate time history data of the discrete detector to obtain the DAS-VSP data simulated wave field and output it.

[0015] Furthermore, in step 1, a parametric model of the well trajectory is established using quadratic Bézier curves:

[0016]

[0017] in, As the starting point of the trajectory, The endpoint of the trajectory, As control points, These are the parameterized variables of the Bézier curve, used to characterize the normalization progress of the curve from the starting point to the ending point. Control points By calculating the midpoint between the starting and ending points of the trajectory And determined by offset along the vertical direction: , For curvature intensity parameters, This is the normalized normal vector. That is, the control function is... .

[0018] Furthermore, in order to eliminate possible overlapping grid points, step 1 eliminates possible overlapping grid points through uniqueness verification.

[0019] Furthermore, in step 1, the non-vertical well grid coordinates are stored in MAT file format.

[0020] Furthermore, in step 2, the tangent direction parameters at each grid point are calculated using a numerical differential method: the coordinate changes between adjacent grid points are calculated using the central difference method; the boundary points are processed using forward / backward difference; and then all coordinate changes are normalized to obtain the unit tangent vector at each grid point.

[0021] Furthermore, in step 2, the calculation of the tangent direction parameters at each grid point using the numerical differentiation method specifically includes:

[0022] Convert grid point coordinates to actual physical coordinates Then, based on the central difference method, the two coordinate components are... and Find the first derivatives respectively: , Where i represents the coordinate index of the current grid point;

[0023] Unit tangent vector of the current grid point , For actual physical coordinates The angle between the tangent and the horizontal axis.

[0024] In this invention, the non-straight well grid is the coordinate in a two-dimensional vertical profile, where x represents the horizontal direction and z represents the vertical direction.

[0025] Furthermore, in step 3, the axial strain rate of each non-vertical well grid point is calculated using the staggered grid finite difference method, specifically including:

[0026] Obtain the three-component velocity field at the grid node: , , ;

[0027] The partial derivatives of the three-component velocity field in the x, y, and z directions are obtained through spatial difference operations: , , , , , , , , ;

[0028] Calculate the axial strain rate of non-vertical well grid points in three-dimensional space. ;

[0029]

[0030] in, The angle between the optical fiber and the horizontal plane.

[0031] Furthermore, the axial strain rate of the non-vertical well grid points in step 3 Two-dimensional spatial form: .

[0032] Furthermore, in step 3, the physical quantities recorded at each detector location include: velocity field. and and axial strain rate .

[0033] Furthermore, step 4, the gauge-length averaging process for the discrete detector strain rate time history data includes:

[0034] Gauge length segmentation: The geophone point set is dynamically segmented along the well trajectory. During segmentation, the distance between adjacent geophones is accumulated until it reaches or exceeds the preset gauge length for the first time (to ensure that the total length of each segment is approximately equal to the gauge length requirement). The resulting segments are used as gauge length segments. The number of geophones and their position indexes contained in each marked segment are recorded.

[0035] Determination of the center representative point: For each gauge length segment, calculate the cumulative distance distribution of each detector within it, and select the detector closest to the geometric center of the segment as the center representative point of the current gauge length segment, so as to use the coordinates of the center representative point as the reference for the position of the gauge length segment.

[0036] Strain rate averaging calculation: The strain rate time history data of all detector points contained in each gauge length segment are arithmetically averaged to generate the composite strain rate response of that gauge length segment.

[0037] The output DAS-VSP data simulates the wavefield including the coordinates of the center representative point of the gauge length segment and the synthetic strain rate response (i.e., the average strain rate), thereby accurately simulating the strain rate measurement results of the actual DAS system in gauge length units.

[0038] The technical solution provided by this invention brings at least the following beneficial effects:

[0039] The method proposed in this invention achieves accurate simulation of DAS-VSP data under non-vertical well conditions by establishing a coupled computational framework of a parametric model of the non-vertical well trajectory and numerical simulation of the elastic wave field. Through data simulation and analysis, the method of this invention has the following significant advantages:

[0040] (1) Non-vertical well modeling and wave field coupling

[0041] This invention establishes a complete simulation system for the elastic field response of non-vertical wells using DAS-VSP. It describes the curved well trajectory by designing a curve function, combining the geometric parameters of the well trajectory with numerical simulation of the elastic wave field. The finite difference method (such as the 8th-order finite difference method) is used to solve the velocity-stress elastic wave equation. A tangent projection algorithm is introduced to convert the velocity field into an optical fiber axial strain rate response, enabling accurate calculation of the distributed optical fiber strain rate response and solving the problem that traditional vertical well simulation methods cannot handle the influence of well inclination changes.

[0042] (2) High-precision simulation of multi-physics coupling

[0043] By constructing the physical transformation relationship between the velocity field and the strain rate field, precise coupling between the elastic wave field and the fiber optic sensing characteristics is achieved. Dynamic gauge length segmentation is employed, and the gauge length segment averaging method is introduced to realistically simulate the spatial sampling characteristics of an actual DAS system. Attached Figure Description

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

[0045] Figure 1 This is a schematic diagram of the processing procedure of the method proposed in the embodiment of the present invention;

[0046] Figure 2 This is a schematic diagram showing the velocity model and the location of the seismic source / detector.

[0047] Figure 3 The image shows the original horizontal component of the velocity field and the shot gather record after adding AGC gain; where (3a) is the actual vertical depth and (3b) is the measured vertical depth.

[0048] Figure 4 The image shows the original vertical component of the velocity field and the shot gather record after adding AGC gain; where (4a) is the actual vertical depth and (4b) is the measured vertical depth.

[0049] Figure 5 The original strain rate and the shot gather record after adding AGC gain are shown for a gauge length of 33m; where (5a) is the actual vertical depth and (5b) is the measured vertical depth.

[0050] Figure 6 A snapshot of the wave field at 0.9s;

[0051] Figure 7 The original strain rate and the shot gather record after adding AGC gain are shown for a gauge length of 21m; where (7a) is the actual vertical depth and (7b) is the measured vertical depth.

[0052] Figure 8 The original strain rate and the shot gather record after adding AGC gain are shown for a gauge length of 73m; where (8a) is the actual vertical depth and (8b) is the measured vertical depth. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be described in detail and completely 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 them. Generally, the components of the embodiments of the present invention described and shown in the accompanying drawings can be arranged and designed using different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely represents selected embodiments of the present invention.

[0054] Existing studies mostly rely on numerical simulations based on the assumption of vertical wells, failing to fully consider the impact of complex well configurations such as deviated and horizontal wells commonly found in actual exploration on wave field propagation and fiber optic sensing characteristics. Addressing the current limitation of DAS-VSP wave field response simulation methods to vertical well models, which struggle to accurately simulate wave field propagation characteristics in curved well trajectories, this invention provides an axial strain rate field simulation method for distributed fiber optic DAS-VSP based on strain tensor projection under non-vertical well deployment conditions. This method establishes a parameterized well trajectory model based on curve functions, combines it with the first-order velocity-stress elastic wave equation (isotropic medium), and employs an explicit finite difference method with staggered grids for numerical solution, establishing a two-dimensional velocity model for non-vertical well DAS-VSP wave field response simulation framework. While retaining the computational efficiency advantages of traditional vertical well simulation methods, this method introduces well trajectory inclination parameters to accurately simulate the coupling process of elastic wave propagation and fiber optic strain rate response in deviated well environments. Specifically, by improving the discretization scheme of the detector positions within the grid, the detector grid array can better adapt to the geometric characteristics of non-vertical wells. Simultaneously, the fiber strain response calculation module is optimized to more accurately reflect the actual response characteristics of distributed optical fibers to elastic wave fields under non-vertical well conditions. This invention's method, specifically optimized for non-vertical well scenarios, effectively fills a significant gap in the practical application of current DAS-VSP simulation technology. This method provides a more reliable numerical experimental platform for the interpretation and processing of DAS-VSP data under complex well conditions such as deviated wells, and has practical value in promoting the application of distributed optical fiber sensing technology in non-vertical well seismic exploration. This invention's method ensures the numerical stability of wavefield simulation in high-angle well sections by constructing a PML (Perfectly Matched Layer) absorbing boundary condition and a stability control mechanism. Experimental results show that this invention's method can accurately characterize the impact of well inclination changes on optical fiber-received seismic signals, especially effectively simulating polarization variation characteristics in curved well sections. This invention provides a reliable numerical experimental method for the interpretation and forward modeling of DAS-VSP data under non-vertical well conditions, and provides important technical support for fiber optic seismic monitoring of non-vertical wells in oil and gas exploration.

[0055] In one embodiment, the axial strain rate field simulation method of distributed optical fiber DAS-VSP based on strain tensor projection provided by the present invention under non-straight well layout conditions mainly includes four steps: (1) layout and storage of non-straight well grid coordinates; (2) calculation and storage of tangent parameters at each grid point of the non-straight well; (3) forward modeling output data; (4) processing of synthetic data.

[0056] Specifically, step (1) involves pre-designing the shape and spatial distribution of the well trajectory in the velocity model and establishing a geometric description of the fiber optic sensing array using mathematical modeling methods. This step realizes the transformation from the theoretical well trajectory to discrete grid coordinates usable for numerical calculation, providing geometric input conditions for subsequent elastic wave field simulation. This process includes the following specific implementation steps:

[0057] (1-a) Parametric Modeling of Well Trajectory: Based on the spatial range of the velocity model and the exploration target, a continuous mathematical function describing the non-vertical well trajectory is established. First, the spatial coordinates of the trajectory's starting and ending points are determined as boundary conditions. Then, the degree of well curvature and morphological characteristics are adjusted by introducing control functions. The modeling process needs to ensure the first-order geometric continuity of the trajectory, making the well path transition smoothly, while controlling the maximum curvature to not exceed the engineering allowable range. The final parametric equation should accurately reflect the geometric characteristics of the designed well type and provide a mathematical basis for subsequent discretization processing.

[0058] (1-b) Gridded Coordinate Generation and Storage: The process of converting continuous well trajectories into discrete computational grid coordinates first involves uniform sampling along the parametric curve to generate a dense set of location points. These physical coordinates are then mapped to a predefined computational grid system, and the coordinate values ​​of each point are rounded to integer grid indices. To ensure all coordinates are within the valid computational domain, boundary points are constrained, and uniqueness checks are used to eliminate any overlapping grid points. The final generated detector location coordinates are stored in standard MAT file format.

[0059] Step (2) specifically involves calculating and storing the tangent direction parameters at each point on the well trajectory based on the established discrete grid coordinates of the well trajectory. This step provides the necessary geometric direction information for subsequent elastic wave field simulation, ensuring that the well-direction strain rate at each discrete point on the optical fiber can be accurately obtained. This process consists of the following two steps, and the specific flow is shown below:

[0060] (2-a) Tangent direction parameter calculation: Based on the discretized well trajectory grid coordinates, the tangent vector at each grid point is calculated using the numerical differentiation method. First, the coordinate change between adjacent grid points is obtained through the central difference method, and then the obtained vector is normalized to obtain the unit tangent vector. Special treatment of boundary points needs to be considered during the calculation process to ensure that the tangent direction is continuous and consistent throughout the entire trajectory.

[0061] (2-b) Direction Parameter Storage: The calculated tangent direction parameters are converted into cosine and sine function values ​​to construct a direction parameter matrix. This matrix corresponds one-to-one with the position coordinate matrix, fully describing the local well trajectory direction characteristics at each grid point. The final direction parameters are output and stored in a file, containing complete tangent direction information.

[0062] Step (3) specifically involves performing an elastic wave field numerical simulation and outputting strain rate response data based on the well trajectory grid coordinates and tangent direction parameters generated in the previous two steps. This process includes the following specific implementation steps:

[0063] (3-a) Wavefield simulation input preparation: Input the well trajectory grid coordinate file and tangent direction parameter file as initial conditions into the forward modeling program. The program automatically establishes a spatial discrete system that matches the velocity model, ensuring that the positions of each detector are accurately mapped to the computational grid nodes.

[0064] (3-b) Strain rate field calculation: Based on the velocity-stress elastic wave equation, the wave field propagation is solved using the staggered grid finite difference method. During the calculation, the three-component velocity field data at each grid node are acquired in real time. , and And the partial derivatives of the three-component velocity field in the x, y, and z directions are obtained through spatial difference operations. , , , , , , , , Based on the principles of elasticity, the strain rate component along the well trajectory axis is calculated by utilizing the conversion relationship between the velocity field and the strain rate field, combined with the tangent direction parameters at each point.

[0065] (3-c) Data output and storage: The axial strain rate time series at each detector location (each grid coordinate corresponds to one detector location) is organized in the form of gathers and output as a standardized data file.

[0066] Step (4) specifically involves gauge-length averaging of the discrete detector strain rate data output from the forward modeling to simulate the data acquisition characteristics of a real distributed fiber optic sensing system. This step realizes the conversion from discrete-point strain response to gauge-length averaged data that conforms to the characteristics of a DAS system. This process specifically includes the following implementation steps:

[0067] (4-a) Gauge length segmentation: Based on the preset gauge length, the geophone point set is dynamically segmented along the well trajectory. During segmentation, the distance between adjacent geophones is accumulated until it reaches or exceeds the gauge length for the first time, ensuring that the total length of each segment is approximately equal to the gauge length requirement. The segment obtained at this time is taken as the gauge length segment, and the number of geophones contained in each gauge length segment and their location index are recorded.

[0068] (4-b) Determination of the center representative point: For each gauge length segment, calculate the cumulative distance distribution of each detector within it, select the detector closest to the geometric center of the segment as the center representative point of the current gauge length segment, and use the coordinates of the center representative point as the reference for the position of the gauge length segment.

[0069] (4-c) Strain rate averaging calculation: Based on the segmented results, the strain rate time history data of all detector points contained in each gauge length segment are arithmetically averaged to generate the synthetic strain rate response of that gauge length segment. The final output is a dataset containing the coordinates of the center point of each gauge length segment and the average strain rate, accurately simulating the strain rate measurement results of the actual DAS system in gauge length units.

[0070] In one embodiment, the implementation process of the axial strain rate field simulation method based on strain tensor projection of distributed fiber optic DAS-VSP under non-vertical well layout conditions provided in this embodiment is as follows: Figure 1 As shown, it includes:

[0071] Step S1: Based on the spatial range of the velocity model and the exploration target, implement parameterized modeling of the well trajectory;

[0072] Step S2: Based on uniform sampling along the parameter curve, position mapping and rounding of coordinate values ​​of each sampling point are performed to generate non-straight well grid points, which are then stored as well trajectory grid coordinate files.

[0073] Step S3: Calculate the tangent parameters at the non-vertical well grid points and save the tangent direction parameter file;

[0074] Step S4: Using the non-vertical well grid points and their tangent parameters as input data for wavefield simulation, the well trajectory grid coordinate file and tangent direction parameter file are input into the forward modeling program as initial conditions. This program automatically establishes a spatial discrete system that matches the velocity model, ensuring that each detector position is accurately mapped to the computational grid nodes. Then, based on the velocity-stress elastic wave equation, the staggered grid finite difference method is used to solve the wavefield propagation, obtaining the axial strain rate of each non-vertical well grid point (i.e., detector position). The axial strain rate time series of the detector positions is then used to construct discrete detector strain rate time history data in gather form, completing the forward modeling data output and storage.

[0075] Step S5: Perform gauge-length averaging on the strain rate time history data of the discrete detector to obtain the DAS-VSP data simulated wave field and output it.

[0076] In one embodiment, step S1, the parameterized modeling of the well trajectory, specifically involves:

[0077] like Figure 2As shown, the velocity model used in this embodiment is 1126 (rows) × 1751 (columns) in size, with a total of 6 layers. A high-speed layer with a diagonal interface is included in the middle in the trend of velocity increasing from top to bottom. The velocity of the model is set from 3500m / s to 6768m / s, and the mesh size is 4m.

[0078] To facilitate construction, this embodiment is pre-designed so that the well trajectory presents a smooth arc from the upper left to the lower right, starting from the first layer and ending in the high-velocity layer. Considering the oblique interface characteristics of the high-velocity layer and the wave propagation characteristics, in order to make the wave propagate towards the detector after reflection at the oblique interface and to allow the detector to receive as much information as possible, the entire optical fiber is placed to the left of the velocity model, and the source is placed to the upper right of the optical fiber.

[0079] The start and end positions of the detector are set to (2000m, 500m) and (4500m, 2740m) respectively. The first position is the horizontal coordinate, and the second position is the vertical coordinate, corresponding to discrete grid coordinates of (500, 125) and (1125, 685) respectively. Figure 2 As shown.

[0080] A parameterized model of the well trajectory is established using quadratic Bézier curves, and its mathematical expression is as follows:

[0081] (1)

[0082] Among them, the starting point = (500, 125) grid units, endpoint = (1125, 685) grid units, The control point is t, which is a parameterized variable of the Bézier curve, representing the normalized progress of the curve from the start point to the end point. t=0 corresponds to the start point of the curve, and t=1 corresponds to the end point of the curve. By calculating the midpoint between the start and end points And determined by offset along the vertical direction,

[0083] (2)

[0084] The offset is determined by the curvature intensity parameter. control, The recommended value range is between 0.1 and 1. A larger value results in a more curved curve, while a value of 0 makes the curve a straight line. Ultimately...

[0085] (3)

[0086] in, This is the normalized normal vector.

[0087] In one embodiment, step S2, the generation and storage of gridded coordinates, specifically involves:

[0088] 180 points were uniformly sampled along the Bézier curve, and the discrete grid coordinates were obtained after rounding. Since the model size is 1126 (rows) × 1751 (columns), the coordinate range is constrained. [0,1750], [0,1125] grid units, and a uniqueness check ensures that the coordinates of each point are not repeated. The final 180 detector position coordinates are stored in MAT file format, and the data is organized as an N×2 integer matrix (N is 180, which is the number of detectors), where each row represents the (x,z) grid coordinates of a detector (x is the horizontal coordinate, z is the vertical coordinate). In this embodiment, the curvature intensity parameter is adjusted. The curvature of the wellbore trajectory can be flexibly controlled when When the value is 0, it degenerates into a straight well. When the radius is 0.5, it exhibits a typical medium-curvature non-vertical well trajectory.

[0089] In one embodiment, the calculation and storage of the tangent parameters at each grid point of the non-vertical well in step S3 specifically involves:

[0090] Step S3-a, Calculation of tangent direction parameters:

[0091] This embodiment uses the central difference method to calculate the tangent direction of each detector position. The specific calculation process is as follows: First, the grid coordinates are converted into actual physical coordinates (multiplied by the grid spacing of 4 meters), and the first derivatives are calculated with respect to the x and z coordinates respectively:

[0092] , (4)

[0093] The boundary points are processed using forward / backward difference. and This represents the actual physical coordinates (in meters) of the i-th detector, obtained by multiplying the grid coordinates by the grid spacing (4 meters). The calculated tangent vector is then normalized to obtain the unit direction vector.

[0094] (5)

[0095] in, is the angle between the tangent at the detector and the horizontal axis, and i is the detector number, which ranges from 1 to 180 in this embodiment. This embodiment specifically handles the differential calculation at the trajectory endpoints to ensure that all 180 detectors obtain physically reasonable tangent directions.

[0096] Step S3-b, Direction Parameter Storage: The calculated tangent direction parameters are stored in the form of an N×2 floating-point matrix, where each row of the matrix contains the cosine of one detector. sin Data. In this embodiment, the data is stored in MAT file format. The file data structure is designed to support direct interface with subsequent forward modeling modules, where the direction parameters and detector position coordinates strictly correspond one-to-one to ensure data consistency.

[0097] In one embodiment, step S4 includes:

[0098] Step S4-a, Wavefield Simulation Input Preparation:

[0099] In this embodiment, the input parameter configuration for forward simulation includes settings in both the time domain and the spatial domain. In the time domain, 12,500 time sampling points are used (i.e., the number of sampling points = 12,500), and the time step is set to 0.0002 seconds, satisfying the CFL (Courant-Friedrichs-Lewy) stability condition.

[0100] (6)

[0101] in, The stability coefficients of the difference scheme are given. This represents the maximum velocity value in the velocity model. For time step, grid spacing rice, , Grid spacing in the horizontal and depth directions respectively; the seismic source uses a dominant frequency of 30Hz (…). Ricker wavelet:

[0102] (7)

[0103] In this embodiment, the amplitude Set to 10, t is the time series, representing the time sampling points of the Ricker wavelet (earthquake source signal). In the spatial domain, the grid spacing... The computational domain covers 1126 × 1751 grid points. The epicenter location is precisely located at grid coordinates (1250, 10), corresponding to the actual physical location (5000 meters, 40 meters). Velocity model data is loaded from the preprocessed file. Fiber optic discrete point coordinate data and tangent parameter data for each point are loaded from the files saved in steps S2 and S3, respectively.

[0104] Step S4-a, Strain rate field calculation:

[0105] This embodiment uses the complete first-order velocity-stress elastic wave equation for forward modeling:

[0106] (8)

[0107] The constitutive relation is:

[0108] (9)

[0109] in, Let be the components of the particle's velocity in the x and z directions. For the normal stresses in the x and z directions, Let be the tangential stress in the xz plane, and t be the time series. Regarding the medium parameter settings, the shear wave velocity... = P-wave velocity / 1.8, density field Assuming a uniform distribution, the elasticity matrix coefficients are calculated as follows:

[0110] , , (10)

[0111] Spatial partial derivatives are calculated using an 8th-order precision finite difference scheme:

[0112] (11)

[0113] This embodiment uses 8th-order precision, with order m=8. The corresponding difference coefficients are given, f is the objective parameter for the difference, and i and j represent the horizontal and vertical coordinates of f, respectively. The boundary treatment uses a PML absorbing boundary with an attenuation coefficient of:

[0114] (12)

[0115] in, To calculate the distance from grid points in the absorption layer to the interface within the absorption layer, The thickness of the absorption layer, This represents the theoretical reflection coefficient (in this embodiment, the value is 10). -6 ).

[0116] The strain rate field is calculated through velocity field transformation. For a vertical well, the relationship between the velocity field and the strain rate field is as follows:

[0117] (13)

[0118] in, Let be the strain rate. Essentially, it is the partial derivative of the velocity component along the fiber direction with respect to the fiber direction. Therefore, for a well with an angle, in three-dimensional space:

[0119] (14)

[0120] in, This represents the velocity component along the fiber direction (the projected velocity). The angle between the optical fiber and the horizontal plane (pitch angle). This is the azimuth angle of the optical fiber in the horizontal plane (the angle between it and the x-axis). Let be the strain rate after rotation (velocity gradient along the fiber direction), and dl be the small displacement increment along the fiber direction. In two-dimensional space:

[0121] (15)

[0122] in, The angle between the tangent at a point on the optical fiber and the horizontal axis represents the local direction of the optical fiber, which is why the tangent parameters at each point need to be calculated in advance.

[0123] Step S4-c, Data Output and Storage:

[0124] The data output system design in this embodiment takes into account the needs of subsequent processing and analysis. The main output includes two types of data: wavefield snapshots and gather records. Wavefield snapshots are generated every 50 seconds. Each snapshot is generated once and simultaneously displays the horizontal velocity component. and vertical velocity component The instantaneous wave field is displayed and overlaid with a velocity model and detector location distribution, facilitating observation of wave propagation and boundary absorption effects. The gather record contains 1250 time sampling points (sampling interval 10). Complete wavefield data, where each time point contains three physical quantities at 180 detector locations: , and (All are 1250×180 arrays). This embodiment uses the MAT file format for data storage. This information is stored together with the wavefield snapshot to ensure data integrity and traceability, facilitating subsequent processing and analysis. , The effect after visualization and adding gain is shown in the image below. Figure 3 , Figure 4 As shown.

[0125] In one embodiment, step S5 includes:

[0126] Step S5-a, Gauge length segmentation:

[0127] This embodiment uses a dynamic accumulation method to segment the fiber optic gauge length. The specific processing flow is as follows: First, the detector grid coordinates are converted to actual physical coordinates (multiplied by the grid spacing of 4 meters). Then, starting from the first detector, the spacing between adjacent detectors is accumulated along the well trajectory direction. When the accumulated length reaches the preset gauge length (33 meters in this example), a segment is completed. The segmentation algorithm calculates the spacing between adjacent points using the following formula:

[0128] (16)

[0129] in, This represents the distance between adjacent detector points. The cumulative length is calculated using the following formula:

[0130] (17)

[0131] Here, the subscript k is the segment identifier. The above process ensures that the actual length of each segment is as close as possible to, but not less than, 33 meters, while maintaining continuity between segments (the starting point of the next segment is the ending point of the previous segment). In the Python implementation, this process can be achieved through a while loop and conditional statements, ultimately generating a segment list containing all detector indices.

[0132] Step S5-b, Determine the center representative point:

[0133] For each gauge length segment, the center representative point is determined by the following steps: First, the cumulative distance distribution of each detector within the segment is calculated:

[0134] (18)

[0135] in, Let be the cumulative distance from the j-th detector in this segment to the start of the segment, where j is the number of detectors in this segment. This represents the total number of detectors. Then, calculate the distance from the midpoint of each segment to the starting point of each segment:

[0136] (19)

[0137] in, This represents the total length of the current segment. This is the distance from the midpoint of the segment to the starting point of the segment. Finally, the detector with the cumulative distance closest to the midpoint is selected as the center representative point.

[0138] (20)

[0139] in, This provides the index of the central representative point within the entire sequence of 180 detectors. During MATLAB processing, segmentation information is obtained by loading the `geophone_segments.mat` file. The `segment_counts` array records the number of detectors in each segment, while the `segment_centers` array stores the coordinates of the center point of each segment. The depth coordinates of the central representative point are converted to physical coordinates by multiplying by the grid spacing.

[0140] (twenty one)

[0141] Step S5-c, Strain rate averaging calculation:

[0142] The strain rate data within each gauge length segment are averaged over time using the following formula:

[0143] (twenty two)

[0144] in, Let be the average strain rate of the k-th segment at time t. Let be the original strain rate of the i-th detector at time t. This represents the set of detector indices contained in the k-th segment. This represents the number of detectors in the segment. In the implementation, all segments are iterated through, and the mean function is called at each time point to calculate the segment average, ultimately generating an average strain rate matrix with dimensions of 1250 × the number of segments N.

[0145] like Figure 5 As shown, this is the average strain rate data output with a gauge length of 33m. In addition to the graph with the actual vertical depth of the detector as the abscissa, this embodiment also adds a graph with the measurement depth along the well trajectory (the well's starting point is zero depth) as the abscissa, facilitating multi-angle data observation. Comparative observation. Figure 3 , Figure 4 The initial arrival wave of the detector at a depth of 500m to 1000m at 0.9s, supplemented by... Figure 6 The wavefield snapshot at this moment shows that this segment of the first arrival wave is dominated by the horizontal component, with a very weak vertical component. Since the trajectory of this segment of the well is nearly vertical, the first arrival wave is approximately perpendicular to the incident direction for this segment. Further observation... Figure 5 The energy in this segment of the strain rate data is very weak, which is consistent with the characteristic of optical fibers that are sensitive to signals along the axial direction but insensitive to signals in the vertical direction.

[0146] Besides the effect of the incident angle on the fiber optic signal response, the effect of different gauge lengths on the fiber optic signal response can also be observed. For example... Figure 7 and 8The gauge lengths were changed to 21m and 73m respectively, and the corresponding number of segments was increased and decreased accordingly. (Comparison) Figure 5 It can be clearly observed that an excessively large gauge length will reduce the resolution of the data, which is consistent with the actual characteristics of optical fiber.

[0147] In this embodiment of the invention, a complete simulation system for the elastic field response of a non-vertical well DAS-VSP is constructed by introducing a non-vertical well geometric model and establishing the physical transformation relationship between the velocity field and the strain rate field. This method uses the velocity-stress elastic wave equation to describe wave field propagation in an isotropic medium and employs an 8-2 staggered grid explicit finite difference algorithm for numerical solution, focusing on solving the key problem of fiber strain rate response calculation under non-vertical well conditions. During implementation, the quantitative relationship between velocity field components and the axial strain rate of the distributed fiber is utilized to accurately calculate the average strain rate within the gauge length along the well, thus completely simulating the response process of the actual DAS system to the elastic wave field. Numerical experiments show that this method can effectively simulate the elastic field response of a distributed fiber DAS-VSP in a well under non-vertical well deployment conditions, providing a reliable numerical simulation tool for the interpretation of DAS-VSP data in complex structural areas.

[0148] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

[0149] The above descriptions are merely some embodiments of the present invention. Those skilled in the art can make various modifications and improvements without departing from the inventive concept of the present invention, and these all fall within the scope of protection of the present invention.

Claims

1. A method for simulating the axial strain rate field of a distributed fiber optic DAS-VSP based on strain tensor projection under non-vertical well layout conditions, characterized in that, Includes the following steps: Step 1: Setting up and storing non-vertical well grid coordinates; Based on the spatial range of the velocity model and the exploration target, a continuous mathematical function describing the non-vertical well trajectory is established to obtain a well trajectory parameterization model to obtain a continuous well trajectory. The well trajectory parameterization model uses the spatial coordinates of the trajectory start and end points as the boundary conditions of the continuous well trajectory, and adjusts the degree of well curvature and morphological characteristics through control functions. The continuous well trajectory is converted into a discrete form and the grid coordinates are calculated: uniform sampling is performed along the parametric curve to generate a dense set of location points; then the physical coordinates corresponding to each location point in the location point set are mapped to a predefined computational grid system, and the coordinate values ​​of each location point are rounded down to convert them into integer grid indices to obtain non-vertical well grid points; Step 2: Calculate and store the tangent parameters at non-vertical well grid points; Traverse each non-vertical well grid point and calculate the tangent direction parameters at each grid point; The tangent direction parameters at each grid point are converted into cosine and sine function values ​​to construct a direction parameter matrix. Each element of the direction parameter matrix is ​​used to characterize the local well trajectory direction characteristics of a non-vertical well grid point. Step 3: Perform numerical simulation of the elastic wave field based on the direction parameter matrix to calculate the axial strain rate of each non-vertical well grid point; A single non-vertical well grid point is used to characterize a geophone location, and the axial strain rate time series of the geophone location is used to construct discrete geophone strain rate time history data in the form of gathers; Step 4: Perform gauge-length averaging on the strain rate time history data of the discrete detector to obtain the DAS-VSP data simulated wave field and output it.

2. The method as described in claim 1, characterized in that, In step 1, a parametric model of the well trajectory is established using quadratic Bézier curves: ； in, As the starting point of the trajectory, The endpoint of the trajectory, As control points, These are the parameterized variables of the Bézier curve, used to characterize the normalization progress of the curve from the starting point to the ending point. Control points By calculating the midpoint between the starting and ending points of the trajectory And determined by offset along the vertical direction: , For curvature intensity parameters, This is the normalized normal vector.

3. The method as described in claim 1, characterized in that, Step 1 also includes: using uniqueness verification to eliminate overlapping grid points.

4. The method as described in claim 1, characterized in that, In step 1, the coordinates of the non-vertical well grid are stored in MAT file format.

5. The method as described in claim 1, characterized in that, In step 2, the tangent direction parameters at each grid point are calculated using the numerical differential method: the coordinate changes between adjacent grid points are calculated using the central difference method; the boundary points are processed using forward / backward difference; and all coordinate changes are then normalized to obtain the unit tangent vector for each grid point.

6. The method as described in claim 5, characterized in that, Step 2, specifically calculating the tangent direction parameters at each grid point using the numerical differentiation method, includes: Convert grid point coordinates to actual physical coordinates Then, based on the central difference method, the two coordinate components are... and Find the first derivatives respectively: , Where i represents the coordinate index of the current grid point; Unit tangent vector of the current grid point , For actual physical coordinates The angle between the tangent and the horizontal axis.

7. The method as described in claim 1, characterized in that, In step 3, the axial strain rate of each non-vertical well grid point is calculated using the staggered grid finite difference method, specifically including: Obtain the three-component velocity field at the grid node: , , ; The partial derivatives of the three-component velocity field in the x, y, and z directions are obtained through spatial difference operations: , , , , , , , , ; Calculate the axial strain rate of non-vertical well grid points in three-dimensional space. ; ； in, The angle between the optical fiber and the horizontal plane.

8. The method as described in claim 7, characterized in that, In step 3, the axial strain rate of the non-vertical well grid points in three-dimensional space is... Replace with two-dimensional space format: .

9. The method as described in claim 7 or 8, characterized in that, In step 3, the physical quantities recorded at each detector position include: velocity field. and and axial strain rate .

10. The method as described in claim 1, characterized in that, Step 4, the gauge-length averaging of the discrete detector strain rate time history data, includes: Gauge length segmentation: The geophone point set is dynamically segmented along the well trajectory. During segmentation, the distance between adjacent geophones is accumulated until it reaches or exceeds the preset gauge length for the first time. The resulting segment is taken as the gauge length segment. The number of geophones and their position indexes contained in each marked segment are recorded. Determination of the center representative point: For each gauge length segment, calculate the cumulative distance distribution of each detector within it, and select the detector closest to the geometric center of the segment as the center representative point of the current gauge length segment; Strain rate averaging calculation: The strain rate time history data of all detector points contained in each gauge length segment are arithmetically averaged to generate the composite strain rate response of that gauge length segment. The output DAS-VSP data simulates the wavefield, including the coordinates of the center representative point of the gauge length segment and the synthetic strain rate response.

Images