Walkaway-VSP longitudinal wave velocity modeling method based on stratigraphic structure information constraint

Abstract

The invention discloses a Walkaway-VSP longitudinal wave velocity modeling method based on stratigraphic structure information constraint, is suitable for the field of oil-gas exploration, and aims to solve the problems that in the prior art, the tectonic boundary depicting capacity is insufficient, and seismic section prior information is not fully utilized in a complex tectonic area. The method comprises the following steps: firstly, correcting a ground earthquake initial velocity model by combining high-precision longitudinal wave velocity information provided by a zero offset VSP, improving the accuracy of the initial velocity model, extracting a stratum boundary from an earthquake offset profile through an edge detection algorithm, and accurately obtaining a stratigraphic dip angle by combining anisotropic diffusion filtering optimization structure tensor calculation. In the inversion process, stratum edge and dip angle information is introduced into an inversion objective function as regularization constraints, the change direction of the velocity model gradient is limited, and the continuity and consistency of the model in the construction direction are enhanced. According to the method, the structure prior information in the seismic migration profile is fully utilized, and the speed modeling precision under the complex geological condition is remarkably improved.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas exploration, and specifically relates to a Walkaway-VSP P-wave velocity modeling method based on stratigraphic structure information constraints. Background Technology

[0002] Velocity inversion is crucial in pre-stack depth migration imaging, providing a solid foundation for oil and gas development, reservoir prediction, and geological interpretation. As exploration targets advance into deeper and more complex geological regions, higher accuracy requirements are being placed on velocity inversion. However, the traditional Walkaway-VSP first-arrival inversion method, based on the assumption of a homogeneous medium, struggles to capture abrupt velocity changes in heterogeneous media. While smoothing regularization is introduced to improve stability, it is prone to over-smoothing key structural information, failing to adequately preserve fault and structural boundaries. Therefore, the accuracy of traditional methods in complex geological environments cannot meet practical needs, making the research of high-precision inversion methods for complex structural geological environments a key focus in this field.

[0003] To address these issues, researchers have proposed several improvement methods. Clapp et al. (2008) incorporated geological information into reflection tomography and improved the geological rationality and convergence speed of the velocity model by using an asymmetric regularization operator to constrain the propagation along the dip angle of the strata. Wang et al. (2010) improved the vertical accuracy of the velocity model in complex surface areas by fusing small refraction and micrologging data. Dai et al. (2020) significantly improved the accuracy of shallow velocity modeling in complex mountainous areas by combining micrologging constraints and anisotropic multi-directional inversion. In addition, Lu et al. (2023) proposed a Walkaway-VSP first-arrival travel time inversion method based on velocity structure constraints, which integrates geological structure, stratigraphic sequence, and smoothing regularization terms to improve the reliability and accuracy of the inversion model.

[0004] Existing research results indicate that the velocity inversion method still faces the following problems in practical applications:

[0005] 1. Existing Walkaway-VSP inversion methods fail to effectively integrate the high-precision P-wave velocity information provided by zero-bias VSP, affecting the accuracy of the velocity model.

[0006] 2. Although existing methods have made some progress in integrating geological and well logging information, they are still insufficient in extracting and utilizing the stratigraphic structural features contained in seismic migration profiles, which restricts the high-precision characterization of complex geological regions.

[0007] 3. Traditional methods generally rely on isotropic smoothing regularization to improve inversion stability. Although this can suppress the influence of noise, it often leads to the excessive smoothing of key geological structures, which limits the imaging capability of velocity inversion in complex geological environments. Summary of the Invention

[0008] To address the shortcomings of traditional velocity inversion methods in complex geological regions, such as insufficient ability to characterize structural boundaries and inadequate utilization of prior information, this invention proposes a Walkaway-VSP velocity inversion method based on stratigraphic structural feature constraints. This method first combines high-precision P-wave velocity information provided by zero-biased VSP to correct the initial ground seismic velocity model, improving its accuracy. Stratigraphic boundaries are extracted from seismic migration profiles using an edge detection algorithm, and the structural tensor calculation is optimized using anisotropic diffusion filtering to accurately obtain the stratigraphic dip angle. During the inversion process, stratigraphic edges and dip angle information are introduced as regularization constraints into the inversion objective function, effectively limiting the direction of gradient change in the velocity model and enhancing the continuity and consistency of the model along structural directions. This method not only fully utilizes prior structural information from seismic migration profiles but also significantly improves the accuracy of velocity modeling under complex geological conditions, achieving higher accuracy than traditional methods.

[0009] The technical solution of the present invention includes the following steps:

[0010] S1. Correct the initial ground seismic model using high-precision P-wave velocity information from zero-biased VSP.

[0011] S2. Structural tensor and edge detection algorithm to extract stratigraphic dip angle and edge information from seismic profiles;

[0012] S3. Stratigraphic structure information is introduced as a regularization constraint into the inversion objective function to guide the iterative update of the velocity model.

[0013] Step S1 includes the following sub-steps:

[0014] S11. On the constructed two-dimensional velocity grid model, the well velocity curve is first calculated based on the zero offset VSP data, and all depth velocity information is completed by spline interpolation. Then, the initial surface velocity curve at the well location is extracted, Fourier spectrum analysis is performed, and a bandpass filter is designed according to the main frequency band to filter the zero offset VSP velocity, thereby suppressing high-frequency noise and providing a basis for subsequent model correction.

[0015] S12. Interpolate the velocity ratio within the same tectonic unit to expand it into a two-dimensional matrix covering the entire velocity model range. Multiply this matrix point-by-point with the original ground seismic velocity model to obtain the updated velocity model. Finally, calculate the slowness matrix on the updated model to provide a more accurate initial velocity model for inversion.

[0016] Step S2 includes the following steps:

[0017] S21. The anisotropic diffusion filter is used to optimize the calculation of the structure tensor. The structure tensor is constructed by image gradient and eigenvalue decomposition is performed. Feature vectors are extracted to determine the direction of stratigraphic extension and the dip angle of the stratigraphic layer is calculated accordingly, so as to achieve an accurate characterization of the geological structure trend.

[0018] S22. First, the image is normalized to unify the grayscale range. Then, two-dimensional Gaussian filtering is used to suppress noise and retain the main structure. Finally, the Canny algorithm is used to extract the stratigraphic boundary. Combined with gradient direction non-maximum suppression and double threshold strategy, the accuracy and continuity of edge detection are effectively improved.

[0019] Step S3 includes the following steps:

[0020] S31. This method divides the neighborhood set according to the dip direction and edge state of each grid point in the seismic image. For non-boundary points, a larger weight is assigned along the direction of stratum extension, and a smaller weight is assigned perpendicular to the direction of stratum extension. Other directions are assigned basic weights to maintain basic connectivity. For boundary points, a larger weight is assigned if the neighboring point is also a boundary point, and a smaller weight is assigned otherwise. Finally, a regularization matrix is ​​constructed so that the weight of the central grid is the negative of the sum of the weights of the neighboring grids, thereby achieving the continuity constraint of the velocity model in the direction of stratum strike and the effective preservation of structural abrupt changes.

[0021] S32. This method establishes a joint inversion objective function that includes the first arrival travel time difference, initial model constraints, velocity model coarsening term, and stratigraphic structure regularization term. It solves the objective function through algebraic reconstruction method, taking into account both data fitting and geological rationality, and achieves high-precision velocity model update.

[0022] Through verification using a synthesis rate model, this invention demonstrates significant advantages in the following aspects:

[0023] This method uses the high-precision P-wave velocity information provided by zero-biased VSP to correct the ground seismic initial velocity model, significantly improving the accuracy of the initial velocity model.

[0024] By using structural tensors and the Canny edge detection algorithm to extract stratum dip angle and boundary information from seismic profiles, the ability to identify complex structural features such as stratum dip and interlayer interfaces in seismic migration profiles is enhanced.

[0025] By constructing a stratigraphic structure-oriented regularization term, the model is subjected to stronger smoothing constraints in the strike direction and weaker constraints in the vertical direction, effectively avoiding the loss of structural features caused by excessive smoothing in traditional methods. By combining boundary state and dip direction information, the smoothing weights among neighboring grids are dynamically adjusted to adapt to the continuity and abrupt changes in stratigraphic structure, thereby enhancing the geological rationality of the inversion model.

[0026] Compared with traditional methods, this invention significantly improves the accuracy, stability and resolution of velocity model inversion, providing accurate and reliable technical support for velocity modeling under complex geological conditions. Attached Figure Description

[0027] Figure 1 Flowchart of the present invention;

[0028] Figure 2 Walkaway-VSP forward modeling observation system;

[0029] Figure 3 Non-boundary mesh constraints;

[0030] Figure 4 Boundary mesh constraints;

[0031] Figure 5 Marmousi velocity model and observation system setup;

[0032] Figure 6 Marmousi initial velocity model;

[0033] Figure 7 Comparison of the actual VSP velocity in the well with the initial velocity profile obtained from the ground seismic inversion;

[0034] Figure 8 Comparison of the spatial frequency amplitude spectrum of the real model and the initial model in the depth direction;

[0035] Figure 9 The result after applying a bandpass filter to the actual velocity profile of the VSP;

[0036] Figure 10 Marmousi offset profile constructs the boundary;

[0037] Figure 11 Marmousi offset profile extraction of dip angle;

[0038] Figure 12 Inversion results: (a) Inversion results with only spatial smoothing regularization added; (b) Inversion results with stratigraphic structure constraints added.

[0039] Figure 13 Comparison of velocity curves at different distances: (a) comparison of velocity curves at Distance = 3040m; (b) comparison of velocity curves at Distance = 3760m; (c) comparison of velocity curves at Distance = 4200m. Detailed Implementation

[0040] To facilitate understanding of the technical content of this invention by those skilled in the art, the following description, in conjunction with the accompanying drawings, further illustrates the invention.

[0041] like Figure 1 As shown, the method of the present invention includes the following steps:

[0042] 1) Use zero-bias VSP high-precision P-wave velocity information to correct the initial ground seismic model;

[0043] 11) On the constructed two-dimensional velocity grid model, the well velocity curve is first calculated based on the zero offset VSP data, and all depth velocity information is completed by spline interpolation. Then, the initial surface velocity curve at the well location is extracted, Fourier spectrum analysis is performed, and a bandpass filter is designed according to the main frequency band to filter the well velocity, thereby suppressing high-frequency noise and providing a basis for subsequent model correction.

[0044] Establish a grid for the two-dimensional velocity model. Based on the observation system, set up an M×N grid model with S shot points, T geophones, and the well located at the midpoint of the shot points.

[0045] Obtaining the surface seismic velocity model. Based on the established velocity grid model, the velocity model v of the subsurface medium is extracted from seismic data obtained from oil fields and geophysical service companies. surface ,v surface It is an M-row N-column matrix.

[0046] Pick the initial arrival time of the VSP data. Based on the zero offset VSP data of the work area, pick the actual initial arrival time t. vsp .

[0047] Calculate the medium velocity V in the well. vsp The seismic wave propagation velocity between adjacent geophones in the well is calculated based on VSP data, using the distance between adjacent geophones and the first arrival time difference. vsp It is a column vector of R rows.

[0048]

[0049] The velocity in the well is interpolated to M rows, and spline interpolation is used to establish V. vsp The relationship with depth, V vsp The interpolation value is M rows and column vectors.

[0050] Fourier transform is performed on the velocity signal extracted from the well by the initial velocity model to calculate its spectrum; based on the main frequency bands of energy concentration in the spectrum, a Butterworth bandpass filter is designed.

[0051] Using a Butterworth bandpass filter for zero-bias VSP velocity Perform filtering

[0052] Calculation speed ratio. The speed of VSP data calculation. The column of wells included in the surface seismic velocity model The corresponding velocities are compared to obtain a column vector α with M rows of velocity ratios, which is used to correct the ground seismic velocity model.

[0053]

[0054] 12) Interpolate the velocity ratio within the same tectonic unit to expand it into a two-dimensional matrix covering the entire velocity model range. Multiply this matrix point-by-point with the original ground seismic velocity model to obtain the updated velocity model. Finally, calculate the slowness matrix on the updated model to provide a high-resolution initial velocity model for velocity inversion.

[0055] Obtain the tectonic interpretation model. Obtain the tectonic interpretation model of the surface earthquake. In the same geological structure, the interpolation velocity ratio α is an M-row N-column matrix, which is multiplied with the surface earthquake velocity model to update the velocity model.

[0056] V0=αv surface

[0057] Calculate the slowness matrix. Slowness is the reciprocal of velocity.

[0058]

[0059] 2) The structural tensor and edge detection algorithm extracts the dip angle and edge information of the strata in the seismic profile.

[0060] 21) The anisotropic diffusion filter is used to optimize the calculation of the structure tensor. The structure tensor is constructed by image gradient and eigenvalue decomposition is performed. The feature vector is extracted to determine the direction of stratum extension and the dip angle of stratum is calculated accordingly, so as to achieve accurate characterization of the geological structure trend.

[0061] The computational structure tensor introduces anisotropic diffusion filter G. M The formula for calculating the structure tensor based on anisotropic diffusion filters is as follows:

[0062]

[0063] The structure tensor matrix S corresponding to a certain grid point σ It can be decomposed into the product of eigenvalues ​​and eigenvectors. Its decomposition form is:

[0064] S σ =λ u uu T +λ v vv T

[0065] The unit eigenvector v corresponding to a smaller eigenvalue can reflect the strike information of the geological structure at that grid point. The dip angle θ can be calculated by the following formula:

[0066]

[0067] 22) Seismic profile stratigraphic edge calculation. First, the seismic profile image data is normalized to a fixed grayscale range, a two-dimensional Gaussian filter is introduced to smooth the image, and the Canny algorithm is used for edge detection to obtain stratigraphic edge information.

[0068] Seismic profile image data is normalized to a fixed grayscale range to eliminate the impact of amplitude differences in the original data on gradient calculation.

[0069] A two-dimensional Gaussian filter is introduced to smooth the image, suppressing high-frequency noise while preserving the main structural features.

[0070] The Canny algorithm is used for edge detection to process seismic migration profiles.

[0071] The local first-order gradient (G) of the image is obtained by taking the derivative of the image along the x and z directions using the Sobel operator. x G z M(x,z) represents the intensity of the grayscale change at the current point; θ(x,z) is the direction of this change, which is subsequently used for directional nonmaximum suppression, as shown in the following formula:

[0072]

[0073]

[0074]

[0075]

[0076] To obtain high-precision edge lines, the gradient magnitude map is refined, retaining only the local maxima along the gradient direction. For each point, the gradient direction is checked against neighboring pixels to see if it is the maximum value; if not, the pixel is set to zero. This process sharpens the edges and avoids blurry edges.

[0077] The Canny algorithm employs a dual-threshold strategy: pixels above the high threshold are marked as "strong edges"; pixels below the low threshold are discarded; and pixels between the high and low thresholds, if connected to strong edges, are retained as "weak edges".

[0078] 3) Stratigraphic structure information is introduced as a regularization constraint into the inversion objective function to guide the iterative update of the velocity model.

[0079] 31) Stratigraphic Structure Regularization Term Design. In the stratigraphic structure regularization matrix R, the velocity model gradient variation relationship between adjacent grid points is defined. Differential weight allocation reflects geological characteristics such as stratigraphic continuity and fault abrupt changes, thereby suppressing high-frequency oscillations in ill-conditioned problems and guiding the model to update in a geologically reasonable direction. To ensure that the regularization direction is consistent with the natural extension direction of the strata, this method dynamically adjusts the connection strength between pixels by analyzing the dip angle difference and edge state between each grid and its eight neighboring grids in the migration profile.

[0080] For non-boundary points, the connection strategy is first determined based on the main extension direction within the neighborhood of the grid point, such as... Figure 3 As shown, the neighborhood of the central grid consists of grids numbered 1 to 8. Assuming the arrow of the central grid represents the direction of its v-vector, and the characteristic direction of the central grid is 230°, then the extension direction of the central grid is a line segment passing through the center and along the characteristic direction. Set D contains the gray grids 2, 3, 6, and 7 in the diagram. The perpendicular direction of the extension direction, the u-vector, is a dashed line passing through the center and along the characteristic direction. Set E contains the white grids 1, 3, 5, and 8 in the diagram.

[0081] Based on the dip direction of the central grid, two neighborhood sets are defined: the extended neighborhood set D contains gray grids extending along the dip direction; the perpendicular set E contains white grids perpendicular to the dip direction. Stronger constraints are imposed along the stratigraphic strike, while weaker constraints are imposed perpendicular to the stratigraphic strike. Specifically, if the direction vector v of a neighborhood point belongs to set D and its dip angle difference with the central point satisfies Δθ∈[0,Δθ], then... max If a neighboring point belongs to set E, it is assigned a larger weight 'a'; if the neighboring point belongs to set E, it is assigned a smaller weight 'b'; other neighboring grids are assigned a weight 'c' by default, where a > c > b. Its mathematical expression is:

[0082]

[0083] For boundary points, the abrupt changes in fault and stratigraphic boundaries are emphasized. Weights are determined based on whether both the current point and its neighboring points are edge points. If both are edge points, a weight of 'a' is assigned; otherwise, a weight of 'c' is assigned, where 'a' > 'c'. The mathematical expression is:

[0084]

[0085] exist Figure 4In this example, the neighborhood of the central grid consists of grids numbered 1 through 8. Assume that grids numbered 3 and 7 are also boundary points within this neighborhood, while the remaining grids are non-boundary points. In this example, the gray grids numbered 3 and 7, being boundary points along with the central grid, have a strong correlation and are therefore assigned a weight 'a' to maintain strong connectivity. In contrast, grids numbered 1, 2, 4, 5, 6, and 8 are not on the boundary and are therefore assigned only a weight 'c' to maintain basic connectivity.

[0086] Finally, in the stratigraphic structure regularization matrix, the weight W of each non-central grid is... ii This indicates the smoothness of the non-center grid relative to the center grid, while the weight of the center grid is set to the sum of the negatives of the weights of all non-center grids.

[0087] W ii =-∑ j≠i W ij

[0088] By assigning different weights to edge points and non-edge points, key structural features are preserved while avoiding excessive smoothing. Simultaneously, the regularization matrix is ​​dynamically adjusted to match the stratigraphic extension direction, enhancing geological plausibility. Differential weight allocation for local grid characteristics ensures smoothness along the stratigraphic strike, allows abrupt changes perpendicular to the stratigraphic direction, and guarantees that the regularization direction always aligns with the stratigraphic structure extension. This achieves the optimal balance between detailed stratigraphic characterization and overall smoothness under complex geological conditions.

[0089] 32) This method establishes a joint inversion objective function that includes the first arrival travel time difference, initial model constraints, velocity model coarsening term and stratigraphic structure regularization term, and solves it through algebraic reconstruction method.

[0090] Forward modeling of first arrival and ray path calculation. For example... Figure 2 The initial velocity model V0 was forward modeled using the shortest path ray tracing algorithm to obtain the sparse matrix A1 of the first arrival wave ray paths. An inversion objective function was constructed and solved. The objective function for grid tomographic velocity inversion was established, including the first arrival difference, coarsening terms of the velocity model, and stratigraphic structure terms, and solved using an algebraic reconstruction method to ensure model accuracy and optimization effectiveness.

[0091]

[0092] The first term is based on the travel time difference of the first arrival wave, the second term is the velocity model coarsening term, and the third term is the stratigraphic structure constraint term.

[0093] when At that time, we can obtain:

[0094]

[0095] To verify the applicability of this method in complex velocity models, numerical experiments were conducted using the classic Marmousi velocity model. To adapt to the shortest travel time ray tracing algorithm, the model was appropriately smoothed. Figure 5 The Marmousi velocity model was demonstrated, with a grid size of 379×121, each grid cell measuring 20×20 meters. The model has a horizontal length of 7580 meters and a total depth of 2420 meters. 124 excitation points were horizontally arranged on the ground surface, spaced 90 meters apart; and 111 receiver points were vertically arranged at a distance of 3780 meters, spaced 70 meters apart.

[0096] To simulate the low-frequency velocity field obtained from ground seismic processing, the Marmousi model was downsampled and then reconstructed to the actual model size using two-dimensional interpolation, such as... Figure 6 As shown

[0097] Figure 7 This paper presents a comparison between the actual wellbore velocity (VSP) and the initial velocity profile obtained from ground seismic inversion. Although the initial model is relatively consistent with the actual model in terms of overall velocity trends, it lacks the ability to express detailed variations, has limited spatial resolution, and struggles to accurately characterize local structural features. Therefore, it is necessary to further extract the main frequency components to enhance the model's structural expressiveness.

[0098] Therefore, we first performed a fast Fourier transform analysis on the velocity profile in the well to obtain its spatial spectral characteristics in the depth direction. Figure 8 The results show a comparison of the spatial frequency amplitude spectra of the real model and the initial model. It can be seen that the real model contains more high-frequency components, indicating that it has richer stratigraphic details; while the spectral energy of the initial model is mainly concentrated in the low-frequency band, and the structure is relatively smooth.

[0099] Based on the frequency bands where the spectral energy is mainly distributed, a Butterworth-type bandpass filter was designed to highlight key geological structural features and effectively suppress high-frequency noise. Subsequently, this filter was applied to a real VSP model, and its main frequency components were extracted through mean removal and bidirectional filtering.

[0100] Figure 9 The cross-section of the filtered true velocity model is shown. It is evident that the filtered velocity curve effectively removes high-frequency disturbances while preserving the main velocity gradients and formation structure. This result more realistically reflects the macroscopic variation trends of the subsurface medium and can serve as a valid reference for inverting the initial model.

[0101] Stratigraphic boundary information is extracted from the offset profile. Figure 10The diagram shows the tectonic boundaries (black lines) extracted from the offset profile of the Marmousi model using the Canny algorithm. These black lines delineate the edges of the stratigraphic structure and reflect the general outline of the subsurface strata.

[0102] Stratigraphic dip information is extracted from the offset profile. Figure 11 The study demonstrates that the structural dip angle extracted from the offset profile of the Marmousi model using the anisotropic diffusion filter structure tensor algorithm accurately reflects the dip angle information of the strata, and the dip angle remains consistent within the same structure, thus updating the model along the dip angle direction.

[0103] This paper employs a joint first-arrival inversion method based on low-frequency background field and stratigraphic structure constraints for calculation. In this section, the weight W of non-boundary points... ij The boundary point weights are set as follows: a = 5, b = 0, c = 1; W ij In the given settings, a = 5 and c = 1. The inversion result is as follows: Figure 12 As shown in b, the inversion results with only the addition of spatial smoothing regularization were selected as the control group (e.g., Figure 12 (As shown in the image). The comparison reveals that the joint first-arrival inversion method significantly outperforms the method using only spatial smoothing and regularization in terms of velocity model accuracy. Specifically, the joint inversion method more clearly captures the stratigraphic contours, and the model's structural orientation is more consistent with the actual velocity model, conforming to actual geological patterns. Especially in complex tectonic regions, the joint inversion based on low-frequency background fields and stratigraphic structure constraints significantly improves the stability and accuracy of the velocity model, effectively avoiding the over-smoothing problem in traditional methods, thus more accurately reflecting subsurface geological characteristics.

[0104] In analysis Figure 12 When analyzing the velocity inversion results, velocity curves at different distances were compared (see...). Figure 13 The effect of the inversion method. Comparison of the velocity curve at 3040 meters ( Figure 13 (a) At depths approaching 1000 and 2000 meters, the inversion results deviate significantly from the actual velocity curves, failing to align well. This deviation may be due to the following reasons: First, the resolution of ground seismic migration images is limited, making it difficult to accurately capture complex geological structures; second, the uneven distribution of rays leads to insufficient data coverage in some areas; finally, the underdeterminacy of the inversion problem makes it difficult for the velocity model to accurately fit the actual data in certain depth regions. For the velocity curve comparison at 3760 meters (… Figure 13At this location, the inversion results of the zero-offset VSP (vertical seismic profile) are quite ideal. The high number of ray coverages results in velocity curve shapes and values ​​that closely resemble the actual model, demonstrating high accuracy. The good inversion results here are mainly due to the relatively uniform ray distribution and the accurate travel time of the zero-offset first arrival wave, which makes the inversion process more stable. For the velocity curve comparison at 4200 meters (…),… Figure 13 The inversion results in the peri-well area are quite close to the actual velocity curves, but deviations still exist in some depth regions. This deviation may stem from insufficient ray coverage in deeper areas, resulting in sparse data in certain depth segments during the inversion process, thus affecting the accuracy of the inversion. Nevertheless, the overall inversion results are reasonable and can capture the main structure of the formation quite well.

[0105] Overall, the formation structure regularization inversion method based on low-frequency background field constraints significantly improves the stability and accuracy of velocity inversion results, especially in the peri-well region. Compared with traditional methods using only spatial smoothing regularization, the joint inversion method shows significant improvements in both the shape and numerical values ​​of the velocity curves. However, in deeper regions, the inversion accuracy is still somewhat affected by the uneven distribution of rays, thus requiring further optimization of the inversion method to improve the model's accuracy.

[0106] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. A Walkaway-VSP P-wave velocity modeling method based on stratigraphic structure information constraints, characterized in that, Includes the following steps: S1. Correct the initial ground seismic model using high-precision P-wave velocity information from zero-biased VSP. S2. Structural tensor and edge detection algorithm to extract stratigraphic dip angle and edge information from seismic profiles; S3. Stratigraphic structure information is introduced as a regularization constraint into the inversion objective function to guide the iterative update of the velocity model.

2. The method according to claim 1, characterized in that, Step S1 includes the following steps: S11. On the constructed two-dimensional velocity grid model, the well velocity curve is first calculated based on the zero offset VSP data, and all depth velocity information is completed by spline interpolation. Then, the initial surface velocity curve at the well location is extracted, Fourier spectrum analysis is performed, and a bandpass filter is designed according to the main frequency band to filter the zero offset VSP velocity, thereby suppressing high-frequency noise and providing a basis for subsequent model correction. S12. The velocity ratio is interpolated within the same structural unit to expand into a two-dimensional matrix covering the entire velocity model range. This matrix is ​​then multiplied point-by-point with the original ground seismic velocity model to obtain the updated velocity model. Finally, the slowness matrix is ​​calculated on the updated model to provide a more accurate initial velocity model for inversion.

3. The method according to claim 1, characterized in that, Step S2 includes the following steps: S21. Optimize the calculation of the structure tensor using an anisotropic diffusion filter, construct the structure tensor through image gradient and perform eigenvalue decomposition, extract feature vectors to determine the direction of stratigraphic extension, and calculate the dip angle of the stratigraphic layers accordingly to achieve an accurate characterization of the geological structure trend. S22. First, the image is normalized to unify the grayscale range. Then, noise is suppressed by two-dimensional Gaussian filtering. Finally, the geological boundary is extracted. By combining non-maximum suppression in the gradient direction and a double threshold strategy, the accuracy and continuity of edge detection are effectively improved.

4. The method according to claim 1, characterized in that, Step S3 includes the following steps: S31. This method divides the neighborhood set according to the dip direction and edge state of each grid point in the seismic image. For non-boundary points, the weight along the stratum extension direction is greater than the weight perpendicular to the stratum extension direction, and other directions are assigned basic weights to maintain basic connectivity. For boundary points, the weight of neighboring points that are also boundary points is greater than the weight of neighboring points that are not boundary points. Finally, a regularization matrix is ​​constructed so that the weight of the central grid is the negative of the sum of the weights of the neighborhood grids, thereby realizing the continuity constraint of the velocity model in the stratum strike direction and the effective preservation of structural abrupt changes. S32. This method establishes a joint inversion objective function that includes the first arrival travel time difference, initial model constraints, velocity model coarsening term, and stratigraphic structure regularization term. It solves the objective function through algebraic reconstruction method, taking into account both data fitting and geological rationality, and achieves high-precision velocity model update.

5. The method according to claim 1, characterized in that, By utilizing the high-precision P-wave velocity information provided by the zero-bias VSP, the initial velocity model of the ground earthquake can be corrected, significantly improving the accuracy of the initial velocity model.

6. The method according to claim 1, characterized in that, By using structural tensor and edge detection algorithms to extract stratum dip angle and boundary information from seismic profiles, the ability to identify stratum dip and interlayer interfaces in seismic migration profiles is enhanced.

7. The method according to claim 1, characterized in that, By constructing a stratigraphic structure-oriented regularization term, the model is subjected to stronger smoothing constraints in the stratigraphic strike direction and weaker constraints in the vertical direction.

Images