High-precision three-dimensional first-motion wave travel time chromatography velocity inversion method and device

Abstract

The invention discloses a high-precision three-dimensional first-motion wave travel time chromatography velocity inversion method and device, and relates to the field of seismic exploration. A multi-template windward difference algorithm and an undulating surface narrow-band continuation strategy are adopted to determine three-dimensional global forward first-motion wave propagation travel time; ray tracing is carried out according to the three-dimensional global forward first-motion wave propagation travel time, and a first-motion wave propagation path and a Frechet differential quotient matrix are obtained; establishing a three-dimensional initial inversion velocity field constrained by micro-logging information; constructing a tomographic inversion equation based on the apparent slowness constraint and the regularization constraint according to the travel time residual error and the Frechet differential quotient matrix; and solving the tomographic inversion equation to obtain a tomographic inversion slowness solution, and converting the slowness solution into a velocity solution. The problem of non-uniform travel time delay of seismic wave propagation of a complex near-surface structure can be solved, a more accurate superficial layer velocity structure can be obtained, and the quality of seismic data processing and imaging inversion can be improved.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration, and in particular to a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method and device. Background Art

[0002] With the continuous expansion of the oil and gas exploration field, areas dominated by complex surface structures such as piedmont belts, Gobi, and deserts have dramatic surface fluctuations, large lithology changes, extremely thick or eroded weathering layers, and very dramatic lateral velocity changes, which are the difficulties faced by oil and gas exploration and development. Such a complex near-surface structure will cause uneven travel time delays in seismic wave propagation. If its velocity structure cannot be accurately described, it will greatly interfere with subsequent static correction work, mid-deep reflection wave tomography inversion, and migration imaging results. On the one hand, it will lead to incorrect static correction results, and on the other hand, it will cause a large deviation between the inversion ray path and the forward ray path. Therefore, the velocity inversion result will significantly deteriorate the depth imaging effect and even produce imaging illusions. Summary of the invention

[0003] In order to solve the above problems, the inventors have made the present invention, and through specific implementation methods, provide a high-precision three-dimensional first-arrival wave traveltime tomographic velocity inversion method and device.

[0004] In a first aspect, an embodiment of the present invention provides a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method, comprising the following steps:

[0005] The multi-template upwind difference algorithm and the narrow-band extension strategy for undulating terrain are used to determine the propagation travel time of the first arrival wave in the three-dimensional global forward modeling.

[0006] According to the three-dimensional global forward modeling of the first arrival wave propagation travel time, ray tracing is performed, and the first arrival wave propagation path and Frechet derivative matrix are obtained by combining the forward modeling travel time of the inverted initial velocity field;

[0007] Establishing the 3D initial inversion velocity field constrained by micro-logging information;

[0008] According to the travel time residual and Frechet derivative matrix, the tomographic inversion equation based on apparent slowness constraint and regularization constraint is constructed.

[0009] Solve the tomographic inversion equation to obtain a tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

[0010] Specifically, picking up the first arrival wave includes the following steps:

[0011] Using machine learning algorithms, the first arrival waves are picked from the gathers in the work area. When seismic data generated by both explosive sources and controllable vibrators exist at the same time, the first arrival waves are picked by combining the wave peaks picked from the controllable vibrator seismic data and the take-off positions picked from the explosive sources.

[0012] Specifically, the high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method further includes the following steps:

[0013] According to the effective distribution range of rays, the maximum offset distance of seismic data used in first-arrival wave tomography inversion is determined;

[0014] The tomographic velocity inversion is performed using the seismic data whose offset does not exceed the maximum value.

[0015] Specifically, a multi-template upwind difference algorithm and a narrow-band extension strategy for undulating terrain are used to determine the propagation travel time of the three-dimensional global forward first arrival wave, including the following steps:

[0016] The multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid nodes of the initial inversion velocity field;

[0017] A three-dimensional narrow-band extension strategy is adopted to simulate the seismic wavefront expansion process on the undulating surface and determine the travel time of the three-dimensional global grid nodes.

[0018] Specifically, a multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid nodes of the initial inversion velocity field, including the following steps:

[0019] A multi-template fast advancing algorithm is used. With each grid node as the center, six upwind difference templates are established to cover the travel time calculation of 26 adjacent grid nodes around each grid node.

[0020] The travel time of each grid node in the three-dimensional local initial inversion velocity field is determined by the following model:

[0021]

[0022]

[0023]

[0024]

[0025] Among them, x is the current local grid node of the initial inversion velocity field in three dimensions, x v It is the grid point with the smallest travel time in the diagonal direction of the cube surrounded by the six upwind differential templates with the current grid node as the center; v represents the direction, and U1, U2, and U3 are respectively along and The directional derivative of are the unit vectors along the directions of l1l2, p1p2, and q1q2 respectively; l1l2, p1p2, and q1q2 are the intersection points of the grid; T vis the travel time of the corresponding grid point position, in ms; T1, T2 and T3 are the minimum travel times of different template arm directions of each template, in ms; represents the travel time at the grid node (i,j,k), are the travel times at grid nodes (i-1,j,k), (i-2,j,k), (i+1,j,k), (i+2,j,k), (i,j+1,k-1), (i,j+2,k-2), (i,j-1,k+1), (i,j-2,k+2), (i,j+1,k+1), (i,j+2,k+2), (i,j-1,k-1), and (i,j-2,k-2).

[0026] Specifically, a three-dimensional narrow-band extension strategy is adopted to simulate the expansion process of seismic wavefronts on the undulating surface and determine the travel time of the three-dimensional global grid nodes of the initial inversion velocity field, including the following steps:

[0027] According to the travel time of the grid nodes before the expansion of the seismic wavefront on the undulating surface, the grid nodes are divided into four categories, including completion points, far-away points, narrow-band points and points above the surface. The completion points correspond to the grid nodes whose travel time has been calculated and cannot be modified, the far-away points correspond to the grid nodes whose travel time has not been calculated, the narrow-band points correspond to the grid nodes whose travel time has been calculated and can be updated, and the points above the surface correspond to the grid nodes located on the undulating surface and do not participate in the wavefront expansion.

[0028] A three-dimensional narrow-band extension strategy is adopted to simulate the expansion process of seismic wavefront on an undulating surface, determine the travel time of far-away points, update the travel time of narrow-band points, and thus determine the travel time of three-dimensional global grid nodes.

[0029] Specifically, according to the three-dimensional global forward modeling first arrival wave propagation travel time, ray tracing is performed, and the first arrival wave propagation path and Frechet derivative matrix are obtained by combining the forward modeling travel time of the inverted initial velocity field, including the following steps:

[0030] According to the derivative of the slowness of the three-dimensional global grid node travel time of the initial inversion velocity field, a sensitive kernel function of the tomographic velocity inversion is constructed, and the global grid node travel time is back-projected to perform ray tracing;

[0031] Under isotropic medium conditions, taking the detection point position as the initial point, ray tracing is performed in the reverse direction along the travel time gradient according to the preset step size until the earthquake excitation point is reached. The ray path is obtained as follows:

[0032]

[0033] Where α is the ray tracing step size; is the travel time gradient, in ms, which can be obtained through the MSFM algorithm; x nis the current position of the ray tracing; x n+1 is the position of the next ray tracing point;

[0034] The differential formula is used to calculate the travel time gradient in the x, y, and z directions respectively. The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node. According to the travel time gradient of the global grid node, the corresponding Frechet derivative matrix is ​​determined.

[0035] Specifically, the differential formula is used to calculate the travel time gradients in the x, y, and z directions respectively. The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node, including the following steps:

[0036] Use the travel time gradient weights of multiple grid nodes to obtain the x at the current ray tracing position n The gradient value of:

[0037]

[0038]

[0039]

[0040]

[0041]

[0042] Among them, h is the calculation step length;

[0043] Set the termination threshold of ray tracing. When ||x is satisfied n+1 -x n When ||<ε·h, ray tracing is terminated. ε is a constant and is set according to different accuracy requirements.

[0044] Specifically, establishing a three-dimensional initial inversion velocity field constrained by micro-logging information includes the following steps:

[0045] The initial inversion velocity field is interpolated along the layers to obtain a three-dimensional initial inversion velocity field containing micro-logging information in the work area.

[0046] Specifically, according to the travel time residual and the Frechet derivative matrix, a tomographic inversion equation based on the apparent slowness constraint and the regularization constraint is constructed, including the following steps:

[0047] Construct the objective function of the apparent slowness constraint and the regularization constraint,

[0048] Φ(m)=(1-ω)S1+ωS2+τS3

[0049] Among them, Φ(m) is the objective function, m is the slowness, ω is the weight factor between the average slowness mismatch and the apparent slowness mismatch, and S1 is represents the normalization of the true travel time d, and l represents the ray path length; represents the normalization of the forward travel time G(m), S2 is represents the apparent slowness, which is the gradient of the travel time and offset curve, D x is the differential operator of travel time with respect to distance, τ is the smoothing weight factor, S3 is Rm, and R is the regularization operator;

[0050] The steady-state equation that minimizes the objective function is iteratively updated to obtain:

[0051]

[0052]

[0053]

[0054] m k+1 =m k +Δm k (k=1,2,3...,N,)

[0055] Among them, for the kth slowness m k , A k is the average slowness kernel function, B k is the apparent slowness kernel function, ε k I is the variable damping term, Δm k is the iteration step size, and N represents the maximum number of slownesses.

[0056] In a second aspect, an embodiment of the present invention provides a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion device, comprising:

[0057] The global travel time determination module is used to determine the three-dimensional global forward modeling first arrival wave propagation travel time by adopting a multi-template upwind difference algorithm and a narrow-band extension strategy for undulating terrain;

[0058] The tomographic inversion model building module is used to perform ray tracing based on the three-dimensional global forward modeling first arrival wave propagation travel time, and combine the forward modeling travel time of the inverted initial velocity field to obtain the first arrival wave propagation path and the Frechet differential matrix; establish a three-dimensional initial inversion velocity field constrained by micro-logging information; and construct a tomographic inversion equation based on apparent slowness constraints and regularization constraints based on the travel time residual and the Frechet differential matrix.

[0059] The velocity solution determination module is used to solve the tomographic inversion equation, obtain the tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

[0060] Based on the same inventive concept, an embodiment of the present invention provides a computer storage medium, in which computer executable instructions are stored. When the computer executable instructions are executed, the above-mentioned high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method is implemented.

[0061] The beneficial effects of the above technical solution provided by the embodiment of the present invention include at least:

[0062] The present invention can overcome the problem of uneven travel time delay of seismic wave propagation in complex near-surface structures, and obtain a more accurate shallow surface velocity structure through tomographic inversion of first-arrival wave travel times under complex undulating surface conditions, which can improve the quality of seismic data processing and imaging inversion, and can effectively improve the velocity inversion accuracy within 1000m underground, eliminating abnormal "bull's eye" illusions and boundary artifact illusions, and the well point inversion velocity curve is highly consistent with the actual logging data.

[0063] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the written description, claims, and drawings.

[0064] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0066] Figure 1 This is a flow chart of a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method in an embodiment of the present invention;

[0067] Figure 2 This is a rendering of the intelligent first arrival picking of the Fusha 4 wire harness in the Talimu Industrial Zone in an embodiment of the present invention;

[0068] Figure 3 6 Cartesian domain templates of the three-dimensional multi-template rapid marching algorithm (MSFM) in the embodiment of the present invention;

[0069] Figure 4 It is a schematic diagram of the calculation state of the three-dimensional grid node when traveling under the narrow-band extension strategy of the three-dimensional undulating surface in an embodiment of the present invention;

[0070] Figure 5 A comparison diagram of an inversion initial velocity field established by interpolating along the layer using micro-logging information and an inversion initial velocity field established by not using micro-logging information in an embodiment of the present invention;

[0071] Figure 6 A comparison diagram of the velocity field obtained by conventional inversion and the velocity field obtained by the method used in the present invention for the Fusha 4 line bundle in the Talimu industrial area in an embodiment of the present invention;

[0072] Figure 7 A comparison diagram of the ray path distribution diagram of the Fusha 4 line bundle in the Talimu industrial area in an embodiment of the present invention obtained by conventional inversion and the ray path distribution diagram obtained by the method used in the present invention;

[0073] Figure 8 A comparison diagram of the velocity field obtained by conventional inversion in a block of Yingmai in Tarim Industrial Zone in an embodiment of the present invention and the velocity field obtained by the method used in the present invention;

[0074] Fig. 9 It is a comparison diagram of the ray path distribution diagram obtained by conventional inversion of three-dimensional complex data of a certain block in Yingmai, Tarim Industrial Zone in an embodiment of the present invention and the ray path distribution diagram obtained by the method used in the present invention;

[0075] Fig.10 A comparison diagram of formation velocities obtained by inversion before and after optimization in an embodiment of the present invention;

[0076] Fig.11 4 is a block diagram of a high-precision three-dimensional first-arrival wave traveltime tomographic velocity inversion device in an embodiment of the present invention. DETAILED DESCRIPTION

[0077] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0078] In order to solve the problems existing in the prior art, an embodiment of the present invention provides a high-precision three-dimensional first-arrival wave traveltime tomographic velocity inversion method and device.

[0079] The embodiment of the present invention provides a high-precision three-dimensional first arrival wave travel time tomographic velocity inversion method, the process of which is as follows: Figure 1 As shown, the following steps are included:

[0080] Step S1: Use the multi-template upwind difference algorithm and the narrow-band extension strategy on the undulating surface to determine the three-dimensional global forward modeling first arrival wave propagation travel time. Travel time is also called travel time.

[0081] In some specific embodiments, picking up the first arrival wave includes the following steps: using a machine learning algorithm to pick up the first arrival wave from the gather of the work area, and when there are seismic data generated by both the explosive source and the vibrator, combining the wave crest picked up from the vibrator seismic data and the take-off position picked up from the explosive source, picking up the first arrival wave. The receiving record of a single seismic detector is called a seismic trace, and the collection of multiple seismic traces is referred to as a gather. Figure 2 The figure shows the effect of intelligent first arrival picking of Fusha 4 line in Tarim Industrial Area, where the horizontal axis is the track number and the vertical axis is the time in ms. The red line marks the first arrival of earthquake picked up based on machine learning method.

[0082] In some specific embodiments, the high-precision three-dimensional first-arrival wave traveltime tomographic velocity inversion method also includes the following steps: determining the maximum offset distance of the seismic data used for the first-arrival wave tomographic inversion based on the effective distribution range of the rays; and performing tomographic velocity inversion using the seismic data whose offset distance does not exceed the maximum value.

[0083] In some specific embodiments, a multi-template upwind difference algorithm and a narrow-band extension strategy for an undulating surface are used to determine the travel time of the three-dimensional global forward first-arrival wave propagation, including the following steps: in a gradient field, a multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid nodes of the initial inversion velocity field; a three-dimensional narrow-band extension strategy is used to simulate the expansion process of the seismic wavefront on an undulating surface to determine the travel time of the three-dimensional global grid nodes.

[0084] In some specific embodiments, a multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid node of the initial inversion velocity field, including the following steps: using a multi-template rapid advancement algorithm, with each grid node as the center, 6 upwind difference templates are established to cover the travel time calculation of 26 adjacent grid nodes around each grid node, and 6 three-dimensional multi-template rapid advancement algorithm (MSFM) Cartesian domain template reference Figure 3 As shown. The following model is used to determine the travel time of each grid node in the three-dimensional local initial inversion velocity field.

[0085]

[0086]

[0087]

[0088]

[0089] Among them, x is the current local grid node of the initial inversion velocity field in three dimensions, x vIt is the grid point with the smallest travel time in the diagonal direction of the cube surrounded by the six upwind differential templates with the current grid node as the center; v represents the direction, and U1, U2, and U3 are respectively along and The directional derivative of are the unit vectors along the directions of l1l2, p1p2, and q1q2 respectively; l1l2, p1p2, and q1q2 are the intersection points of the grid; T v is the travel time of the corresponding grid point position, in ms; T1, T2 and T3 are the minimum travel times of different template arm directions of each template, in ms; represents the travel time at the grid node (i,j,k), are the travel times at grid nodes (i-1,j,k), (i-2,j,k), (i+1,j,k), (i+2,j,k), (i,j+1,k-1), (i,j+2,k-2), (i,j-1,k+1), (i,j-2,k+2), (i,j+1,k+1), (i,j+2,k+2), (i,j-1,k-1), and (i,j-2,k-2).

[0090] In some specific embodiments, a three-dimensional narrowband extension strategy is adopted to simulate the expansion process of the earthquake wavefront on the undulating surface and determine the travel time of the three-dimensional global grid nodes, including the following steps: according to the travel time of the grid nodes before the expansion of the earthquake wavefront on the undulating surface, the grid nodes are divided into four categories, including completion points, far-away points, narrowband points and points above the surface, the completion points correspond to the grid nodes whose travel time has been calculated and cannot be modified, the far-away points correspond to the grid nodes whose travel time has not been calculated, the narrowband points correspond to the grid nodes whose travel time has been calculated and can be updated, and the points above the surface correspond to the grid nodes located above the undulating surface and do not participate in the wavefront expansion; the three-dimensional narrowband extension strategy is adopted to simulate the expansion process of the earthquake wavefront on the undulating surface, determine the travel time of the far-away points, update the travel time of the narrowband points, and thus determine the travel time of the three-dimensional global grid nodes of the initial inversion velocity field. Reference Figure 4 As shown in the figure, the calculation status of the travel time of the three-dimensional grid node under the three-dimensional undulating surface narrowband extension strategy. Among them, the black dots are completed points, the travel time has been calculated and cannot be modified; the white dots are far away points, the travel time has not been calculated; the gray dots are narrowband points, the travel time has been calculated and can still be updated.

[0091] Step S2: According to the three-dimensional global forward modeling first arrival wave propagation travel time, ray tracing is performed, and the first arrival wave propagation path and Frechet derivative matrix are obtained by combining the forward modeling travel time of the inverted initial velocity field.

[0092] In some specific embodiments, ray tracing is performed according to the three-dimensional global forward modeling first arrival wave propagation travel time to obtain the first arrival wave propagation path and the Frechet derivative matrix, including the following steps: constructing a sensitive kernel function for tomographic velocity inversion according to the derivative of the three-dimensional global grid node travel time of the initial inversion velocity field to the slowness, back-projecting the global grid node travel time, and performing ray tracing; under isotropic medium conditions, taking the detection point position as the initial point, ray tracing is performed in the reverse direction along the travel time gradient direction according to a preset step size until the seismic excitation point is reached, and the ray path is obtained as follows,

[0093]

[0094] Where α is the ray tracing step size; is the travel time gradient, in ms, which can be obtained through the MSFM algorithm; x n is the current position of the ray tracing; x n+1 is the position of the next ray tracing point; the differential formula is used to calculate the travel time gradient in the x, y, and z directions respectively The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node. According to the travel time gradient of the global grid node, the corresponding Frechet derivative matrix is ​​determined.

[0095] In some specific embodiments, the differential formula is used to calculate the travel time gradient in the x, y, and z directions respectively. The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node, including the following steps: using the travel time gradient weights of multiple grid nodes to obtain the current ray tracing position x n The gradient value of:

[0096]

[0097]

[0098]

[0099]

[0100]

[0101] Where h is the calculation step size. Set the termination threshold of ray tracing. When ||x n+1 -x n When ||<ε·h, ray tracing is terminated. ε is a constant and is set according to different accuracy requirements.

[0102] Step S3: Establish a three-dimensional initial inversion velocity field constrained by micro-logging information; construct a tomographic inversion equation based on apparent slowness constraints and regularization constraints according to travel time residuals and Frechet derivative matrices.

[0103] In some specific embodiments, establishing a three-dimensional initial inversion velocity field constrained by micro-logging information includes the following steps: interpolating the initial inversion velocity field along the layer to obtain a three-dimensional initial inversion velocity field containing micro-logging information in the work area.

[0104] In some specific embodiments, a tomographic inversion equation based on apparent slowness constraints and regularization constraints is constructed according to the travel time residual and the Frechet derivative matrix, comprising the following steps:

[0105] Construct the objective function of the apparent slowness constraint and the regularization constraint,

[0106] Φ(m)=(1-ω)S1+ωS2+τS3

[0107] Among them, Φ(m) is the objective function, m is the slowness, ω is the weight factor between the average slowness mismatch and the apparent slowness mismatch, and S1 is represents the normalization of the true travel time d, and l represents the ray path length; represents the normalization of the forward travel time G(m), S2 is represents the apparent slowness, which is the gradient of the travel time and offset curve, D x is the differential operator of travel time with respect to distance, τ is the smoothing weight factor, S3 is Rm, and R is the regularization operator;

[0108] The steady-state equation that minimizes the objective function is iteratively updated to obtain:

[0109]

[0110]

[0111]

[0112] m k+1 =m k +Δm k (k=1,2,3...,N,,)

[0113] Among them, for the kth slowness m k , A k is the average slowness kernel function, B k is the apparent slowness kernel function, ε k I is the variable damping term, Δm k is the iteration step size, and N represents the maximum number of slownesses.

[0114] refer to Figure 5 As shown in the figure, the initial velocity field (corresponding to Figure 5 Figure b) on the right shows the initial velocity field of the inversion without using micro-logging information (corresponding to Figure 5 Comparison with the left figure a) shows that the initial velocity field of the inversion containing micro-logging information is slightly different from the conventional initial velocity field in the shallow layer. In the subsequent tomographic inversion process, the two initial velocity fields will be used for tomographic inversion respectively, and then the effects will be compared.

[0115] Step S4: Solve the tomographic inversion equation to obtain a tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

[0116] refer to Figure 6 As shown, the velocity field obtained by conventional inversion of the Fusha 4 line bundle in the Tarim Industrial Zone is compared with the velocity field obtained by the method used in the present invention. Figure 6 In the left figure a, the conventional first arrival tomography inversion algorithm and the optimized inversion algorithm use the same grid size of 15m×15m×10m. The conventional first arrival tomography inversion results of this block have a serious "bull's eye" phenomenon, including many high-speed anomalies, and the low and falling speed zones and the high-speed layer top interface are not clear. Figure 6 In Figure b on the right, the optimized first-arrival wave tomography inversion results basically eliminate the "bull's eye" phenomenon, and have a good effect on depicting the low and decreasing velocity zones under undulating surface conditions. The speed of the low-velocity zone is about 800m / s, and the speed of the decreasing velocity zone is about 1300m / s. The top interface of the high-speed layer is clear, and the speed is about 2000m / s, which are consistent with the surface survey results.

[0117] In a specific embodiment, a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method comprises the following steps:

[0118] Step S11: Use the intelligent first arrival wave picking method based on machine learning to perform quality control on all the picked first arrival travel times and remove abnormal first arrival picking values, such as 0 value.

[0119] Step S12, using a multi-template upwind difference algorithm and an undulating surface narrowband extension strategy to calculate the travel time of the three-dimensional global first-arrival wave propagation based on the initial inversion velocity field.

[0120] Step S13, back-project the three-dimensional global forward modeled first arrival wave propagation travel time obtained in the previous step, complete ray tracing, and combine the forward modeled travel time of the inverted initial velocity field to obtain the first arrival wave propagation path and the Frechet derivative matrix.

[0121] Step S14: establishing a three-dimensional initial inversion velocity field constrained by micro-logging information.

[0122] Step S15: constructing a tomographic inversion equation based on apparent slowness constraints and regularization constraints on the basis of the travel time residual and the Frechet derivative matrix.

[0123] Step S16, solving the tomographic inversion equation to obtain a tomographic inversion slowness solution, which is further converted into a velocity solution.

[0124] In terms of improving computing efficiency, distributed memory programming MPI (a programming model based on message passing) is adopted, combined with MSFM technology and narrowband extension strategy to calculate the travel time of the forward model and Runge-Kutta algorithm to complete reverse ray tracing; in terms of improving the accuracy of tomographic inversion, on the basis of ensuring the efficiency and accuracy of first-arrival picking, a set of multi-information constrained first-arrival wave tomographic inversion optimization strategies is proposed for the core links in the first-arrival tomographic inversion process, that is, increasing the apparent slowness and regularization to constrain the tomographic inversion objective function, reducing the morbidity and multiple solutions of the objective equation, and using micro-logging information to constrain the initial velocity field of tomographic inversion.

[0125] Specifically, a first-arrival picking method based on machine learning is first used to improve the efficiency and accuracy of first-arrival picking, thereby ensuring the accuracy of the prior information of tomographic inversion during the travel of the first-arrival wave. At the same time, when processing seismic data generated by explosive sources and controllable seismic sources, in order to address the problem of obvious phase differences between the two methods, the present invention accurately picks the first arrival by combining the wave peak of the controllable seismic source data with the starting position of the explosive source.

[0126] The larger the offset distance of the seismic data used for first-arrival wave tomographic inversion, the richer the stratigraphic travel time information contained, and the deeper the inversion depth. However, if the offset distance is too large, on the one hand, the ray path will have a shadow area as the inversion depth increases, causing distortion of the tomographic inversion result and reducing the credibility of the inversion result; on the other hand, if the offset distance is too large, the amount of first arrivals contained in the seismic data will increase, and the probability of first arrival mispickup will increase. If the accuracy of first arrival picking cannot be guaranteed, the tomographic inversion result will be affected from the source. In response to this problem, the present invention only retains seismic data with an offset distance of nearly 3000m for tomographic inversion, which can ensure that the effective distribution range of the rays is about 1000m while ensuring the inversion accuracy of the shallow and middle layer velocity as much as possible.

[0127] In the forward modeling traveltime calculation of the three-dimensional first-arrival wave, the multi-template fast marching algorithm (MSFM) is adopted to establish 6 upwind difference templates, which can cover the traveltime calculation of 26 adjacent grid nodes around the grid, effectively improving the calculation accuracy of the local forward modeling traveltime.

[0128] Take the template F2 in the second-order difference format as an example ( Figure 3 ), U1, U2, U3 can be approximated as:

[0129]

[0130] Where: x v are the grid points with smaller travel time in the diagonal direction; U1, U2, and U3 are the grid points along and The directional derivative of are the unit vectors along the directions of l1l2, p1p2, q1q2 respectively; l1l2, p1p2, q1q2 are the grid intersection points; T v is the travel time of the corresponding grid point position, in ms; T1, T2 and T3 are the minimum travel times of different template arm directions of each template, in ms, which can be expressed as

[0131]

[0132]

[0133]

[0134] Since the multi-template fast advancing algorithm calculates the travel time of the grid node through multiple Cartesian domain templates, it takes into account the travel time information of 26 adjacent nodes around the grid node. Compared with the single template, it can more effectively improve the accuracy of the travel time calculation in the diagonal direction.

[0135] Only local grid node travel time calculation can be realized through the multi-template upwind difference algorithm. The present invention adopts three-dimensional narrow-band extension technology to simulate the seismic wavefront expansion process to realize global travel time calculation: the grid nodes can be divided into four categories according to whether the node completes the travel time calculation during wavefront expansion: ① completed points, the travel time has been calculated and cannot be modified; ② far away points, the travel time has not been calculated; ③ narrow-band points, the travel time has been calculated and can still be updated; ④ points above the surface, located on the undulating surface, do not participate in wavefront expansion.

[0136] Ray tracing is required in the process of travel time tomographic inversion, that is, the differential of travel time with respect to model slowness is used to construct the sensitive kernel function of tomographic inversion. Therefore, after calculating the global travel time through the multi-template rapid marching algorithm (MSFM), the ray tracing can be completed by back-projecting the global travel time calculation results. Under isotropic medium conditions, the propagation direction of seismic waves is the direction of the wavefront normal, that is, the direction of the travel time gradient. Therefore, starting with the forward global travel time obtained by MSFM, taking the position of the detection point as the initial point, ray tracing is performed in the reverse direction along the travel time gradient direction at a certain step size until it reaches the excitation point to obtain the ray path.

[0137]

[0138] Where: α is the step size of ray tracing; is the travel time gradient, in ms, which can be obtained through the MSFM algorithm; xn is the current position of the ray tracing; x n+1 is the position of the next ray tracing point.

[0139] refer to Figure 7 As shown in the figure, the ray path distribution diagram of the Fusha 4 beam in the Tarim area is compared with the ray path distribution diagram obtained by the method used in the present invention. The ray distribution of the conventional tomographic inversion result is messy and unreasonable ( Figure 7 a), the ray path distribution of the optimized inversion result is more reasonable and uniform ( Figure 7 b). From the comparison of ray density, it can be seen that after optimization, the ray distribution has obvious velocity anomalies near the depth of 500m and obvious refraction layers at the depth of 1000m. There is basically no ray path distribution at a depth greater than 1000m, and the inversion accuracy is not high. This is related to the fact that we only use seismic data with an offset distance of less than 3000m.

[0140] Here we take the x direction as an example (the calculation method for the y and z directions is the same), and use the difference formula to calculate the travel time gradient The first-order component of:

[0141]

[0142] The second-order component,

[0143] The third-order component,

[0144] Where: represents the travel time gradient in the x-direction of the coordinate axis, in ms; Δx is the grid spacing in the x-direction, in m; t i,j,k Represents the travel time at the grid node (i,j,k) in ms.

[0145] The travel time gradients of all grid nodes can be obtained through the above formula. In order to further improve the calculation accuracy of the global travel time gradient, the Runge-Kutta algorithm is introduced. By using the Taylor expansion of the approximate ray equation, the travel time gradients of multiple grid nodes are weighted to obtain the x at the current ray tracing position. n The gradient value of:

[0146]

[0147] Where: h is the calculation step length; K1, K2, K3, K4 can be expressed as

[0148]

[0149]

[0150]

[0151]

[0152] By combining equations (7) and (8), the travel time gradient calculation can be completed and the termination threshold of ray tracing can be set, that is, the ray tracing path length is much smaller than the ray tracing step length, satisfying

[0153] ||x n+1 -x n ||<ε·h (9)

[0154] Then the ray tracing ends, where ε is a constant, which can take different values ​​according to different accuracy requirements. Combining the advantages of different order travel time gradient calculation formulas, Runge-Kutta algorithm and flexible adjustment of step size h, this method can meet the requirements of complex surface three-dimensional actual data tomography inversion, and has the advantages of high computational efficiency and stability.

[0155] The present invention uses a variety of information to constrain the target equation when constructing the tomographic inversion target equation. The travel time tomographic inversion gradually updates the slowness model by fitting the true travel time d and the forward travel time G(m). The inversion process can be expressed as:

[0156]

[0157] Where: m is slowness, s / m; Φ(m) is the objective function; It represents the standardization of the true first arrival travel time, in ms; It represents the normalized forward travel time of the model, in ms; l represents the ray path length, in m.

[0158] For 3D travel time tomographic inversion involving massive travel time data, the model grids are large and the insufficient number of rays will lead to underdetermined tomographic inversion equations, that is, the number of equations is less than the number of unknowns, which aggravates the ill-conditioned nature of the inversion problem. Prior information needs to be added to constrain the inversion solution space to speed up the iterative convergence speed, so as to obtain an approximate solution (maximum likelihood solution) that meets the needs of practical applications. To address this problem, the apparent slowness information is added to the objective function to constrain the average slowness of the slowness model inversion. It helps to improve the accuracy of shallow velocity, and inverting the apparent slowness helps to restore the deep velocity structure. Therefore, by minimizing the mismatch between the average slowness and the apparent slowness, the inversion multi-solution is reduced. After introducing the apparent slowness information, the objective function Φ(m) can be expressed as:

[0159]

[0160] or

[0161] Where: D x is the differential operator of travel time with respect to distance; is the gradient of the travel time and offset curve (also called apparent slowness); ω is the weight factor between the average slowness mismatch and the apparent slowness mismatch.

[0162] After introducing the visual slowness information constraint, the present invention adds a regularization information constraint. On the basis of the visual slowness information constraint, Tikhonov regularization is applied to constrain the smoothness of the model, that is, the minimization objective function Φ(m) of formula (12) can be rewritten as:

[0163]

[0164] Where: τ is the smoothing weight factor; C l is an operator that scales the travel time corresponding to the ray length l, and R is a regularization operator (such as a derivative operator, a smoothing operator).

[0165] The Gauss-Newton (GN) method is used to construct the steady-state equations related to the minimization objective function, and the conjugate gradient (CG) algorithm is used to complete each iterative update to obtain

[0166]

[0167]

[0168]

[0169] m k+1 =m k +Δm k (k=1,2,3...,N,) (17)

[0170] Where: A k is the average slowness kernel function; B k is the apparent slowness kernel function; ε k I is the variable damping term.

[0171] Among them, when establishing the initial inversion velocity field, the present invention adopts the method of interpolation along the layer to obtain the inversion initial arrival velocity field containing the micro-logging information in the work area through interpolation. The use of the initial inversion velocity field constrained by micro-logging helps to accelerate the convergence of the objective function and improve the accuracy of the tomographic inversion results.

[0172] refer to Figure 8 The figure shows the comparison of the velocity field obtained by conventional inversion and the velocity field obtained by the method used in the present invention in a block of Yingmai in the Tarim Industrial Zone. Figure 8In the left figure, the conventional inversion velocity field has velocity anomalies in the longitudinal slices, and due to the uneven distribution of rays and the small depth, the boundary layer shadow phenomenon is severe; while in the tomographic inversion velocity field after multi-information constraint optimization, Figure 8 In the right figure, due to the regularization factor, the velocity is smoother, the "bull's eye" abnormal velocity block is eliminated, and the boundary shadow phenomenon is significantly improved.

[0173] refer to Fig. 9 The figure shows the comparison between the ray path distribution diagram obtained by conventional inversion of the three-dimensional complex data of a certain block in Yingmai, Tarim Industrial Zone and the ray path distribution diagram obtained by the method used in the present invention. From the ray density field before optimization (corresponding to the left figure) and after optimization (corresponding to the right figure), it can be seen that the ray distribution inverted after optimization is more uniform than the inverted ray distribution of the conventional algorithm, and the ray density at the same position is greater, which reduces the sparsity of the tomographic inversion equation, and thus the velocity inversion result is more credible.

[0174] refer to Fig.10 The figure shows the formation velocity obtained by inversion before and after quantitative comparison optimization. The velocity curves are extracted at the same position of the two tomographic inversion results and superimposed and compared with the well velocity curve. Fig.10 It can be seen from the above that the inversion velocity before and after optimization is slightly different at a depth of 200-800 m. When the depth is greater than 800 m, the optimized velocity curve is more consistent with the logging velocity curve, which means that the optimized inversion result is more accurate and closer to the actual formation velocity.

[0175] In the above method of this embodiment, the first arrival picking method under complex surface conditions and the first arrival picking of massive seismic data are picked by machine learning methods; the apparent slowness information and regularization are used to constrain the inversion solution for the tomographic inversion target equation, and the initial velocity field of the inversion is constrained by micro-logging prior information to improve the accuracy of tomographic inversion; a variety of gradient solution methods and line search methods are used for the tomographic inversion equation itself to make the tomographic equation iterative convergence as much as possible. Finally, through the accurate picking of seismic first arrivals and the research and improvement of forward and inverse algorithms, high-precision three-dimensional first arrival wave travel time tomographic velocity inversion with multi-information constraints is performed to obtain the accurate velocity structure of seismic waves. It can overcome the problem of uneven travel time delay of seismic wave propagation in complex near-surface structures, and obtain a more accurate shallow surface velocity structure through tomographic inversion of first-arrival wave travel time under complex undulating surface conditions. It can improve the quality of seismic data processing and imaging inversion, and effectively improve the velocity inversion accuracy within 1000m underground, eliminating abnormal "bull's eye" illusions and boundary artifact illusions, and the well point inversion velocity curve is highly consistent with the actual logging data.

[0176] Those skilled in the art can change the above sequence without departing from the protection scope of the present disclosure.

[0177] Another embodiment of the present invention provides a high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion device, such as Fig.11 As shown, including:

[0178] The global travel time determination module uses a multi-template upwind difference algorithm and a narrow-band extension strategy for undulating terrain to determine the three-dimensional global forward modeling first-arrival wave propagation travel time;

[0179] The tomographic inversion model building module is used to perform ray tracing based on the travel time of the first arrival wave propagation in the three-dimensional global forward modeling to obtain the first arrival wave propagation path and the Frechet differential matrix; establish a three-dimensional initial inversion velocity field constrained by micro-logging information; and construct a tomographic inversion equation based on apparent slowness constraints and regularization constraints based on the travel time residual and the Frechet differential matrix.

[0180] The velocity solution determination module is used to solve the tomographic inversion equation, obtain the tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

[0181] Regarding the device in the above embodiment, the specific manner in which each module performs operations has been described in detail in the embodiment of the method, and will not be elaborated here.

[0182] In this embodiment, the problem of uneven travel time delay of seismic wave propagation in complex near-surface structures can be overcome, and a more accurate shallow surface velocity structure can be obtained through tomographic inversion of the first-arrival wave travel time under complex undulating surface conditions, which can improve the quality of seismic data processing and imaging inversion, and can effectively improve the velocity inversion accuracy within 1000m underground, eliminating abnormal "bull's eye" illusions and boundary artifact illusions, and the well point inversion velocity curve is highly consistent with the actual logging data.

[0183] Based on the same inventive concept, an embodiment of the present invention provides a computer storage medium, in which computer executable instructions are stored. When the computer executable instructions are executed, the above-mentioned high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method is implemented.

[0184] Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention shall still fall within the patent scope of the present invention.

Claims

1. A high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method, characterized in that: The following steps are involved: The multi-template upwind difference algorithm and the narrow-band extension strategy for undulating terrain are used to determine the propagation travel time of the first arrival wave in the three-dimensional global forward modeling. According to the three-dimensional global forward modeling of the first arrival wave propagation travel time, ray tracing is performed, and the first arrival wave propagation path and Frechet derivative matrix are obtained by combining the forward modeling travel time of the inverted initial velocity field; Establishing the 3D initial inversion velocity field constrained by micro-logging information; According to the travel time residual and Frechet derivative matrix, the tomographic inversion equation based on apparent slowness constraint and regularization constraint is constructed. Solve the tomographic inversion equation to obtain a tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

2. The method according to claim 1, characterized in that Picking the first arrival wave includes the following steps: Using machine learning algorithms, the first arrival waves are picked from the gathers in the work area. When seismic data generated by both explosive sources and controllable vibrators exist at the same time, the first arrival waves are picked by combining the wave peaks picked from the controllable vibrator seismic data and the take-off positions picked from the explosive sources.

3. The method according to claim 1, characterized in that The high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method further comprises the following steps: According to the effective distribution range of rays, the maximum offset distance of seismic data used in first-arrival wave tomography inversion is determined; The tomographic velocity inversion is performed using the seismic data whose offset does not exceed the maximum value.

4. The method according to claim 1, characterized in that The multi-template upwind difference algorithm and the narrow-band extension strategy for undulating terrain are used to determine the travel time of the three-dimensional global forward first arrival wave propagation, including the following steps: The multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid nodes of the initial inversion velocity field; A three-dimensional narrow-band extension strategy is adopted to simulate the seismic wavefront expansion process on the undulating surface and determine the travel time of the three-dimensional global grid nodes.

5. The method according to claim 4, characterized in that The multi-template upwind difference algorithm is used to determine the travel time of the three-dimensional local grid nodes of the initial inversion velocity field, including the following steps: A multi-template fast advancing algorithm is used. With each grid node as the center, six upwind difference templates are established to cover the travel time calculation of 26 adjacent grid nodes around each grid node. The travel time of each grid node in the three-dimensional local initial inversion velocity field is determined by the following model: Among them, x is the current local grid node of the initial inversion velocity field in three dimensions, x v It is the grid point with the smallest travel time in the diagonal direction of the cube surrounded by the six upwind differential templates with the current grid node as the center; v represents the direction, and U1, U2, and U3 are respectively along and The directional derivative of are the unit vectors along the directions of l1l2, p1p2, and q1q2 respectively; l1l2, p1p2, and q1q2 are the intersection points of the grid; T v is the travel time of the corresponding grid point position; T1, T2 and T3 are the minimum travel times of different template arm directions of each template respectively; represents the travel time at the grid node (i,j,k), are the travel times at grid nodes (i-1,j,k), (i-2,j,k), (i+1,j,k), (i+2,j,k), (i,j+1,k-1), (i,j+2,k-2), (i,j-1,k+1), (i,j-2,k+2), (i,j+1,k+1), (i,j+2,k+2), (i,j-1,,k-1), and (i,j-2,k-2).

6. The method according to claim 4, characterized in that The three-dimensional narrow-band extension strategy is used to simulate the expansion process of seismic wavefronts on the undulating surface and determine the travel time of the three-dimensional global grid nodes, including the following steps: According to the travel time of the grid nodes before the expansion of the seismic wavefront on the undulating surface, the grid nodes are divided into four categories, including completion points, far-away points, narrow-band points and points above the surface. The completion points correspond to the grid nodes whose travel time has been calculated and cannot be modified, the far-away points correspond to the grid nodes whose travel time has not been calculated, the narrow-band points correspond to the grid nodes whose travel time has been calculated and can be updated, and the points above the surface correspond to the grid nodes located on the undulating surface and do not participate in the wavefront expansion. A three-dimensional narrow-band extension strategy is adopted to simulate the expansion process of seismic wavefront on an undulating surface, determine the travel time of far-away points, update the travel time of narrow-band points, and thus determine the travel time of three-dimensional global grid nodes.

7. The method according to claim 1, characterized in that According to the three-dimensional global forward modeling first arrival wave propagation travel time, ray tracing is performed, and combined with the inverted initial velocity field forward modeling travel time, the first arrival wave propagation path and Frechet derivative matrix are obtained, including the following steps: According to the derivative of the slowness of the three-dimensional global grid node travel time of the initial inversion velocity field, a sensitive kernel function of the tomographic velocity inversion is constructed, and the global grid node travel time is back-projected to perform ray tracing; Under isotropic medium conditions, taking the detection point as the initial point, ray tracing is performed in the reverse direction along the travel time gradient according to the preset step size until the earthquake excitation point is reached. The ray path is obtained as follows: Where α is the ray tracing step size; is the travel time gradient, which can be obtained through the MSFM algorithm; x n is the current position of the ray tracing; x n+1 is the position of the next ray tracing point; The differential formula is used to calculate the travel time gradient in the x, y, and z directions respectively. The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node. According to the travel time gradient of the global grid node, the corresponding Frechet derivative matrix is ​​determined.

8. The method according to claim 7, characterized in that The differential formula is used to calculate the travel time gradient in the x, y, and z directions respectively. The first-order, second-order, and third-order components are used to obtain the travel time gradient of the global grid node, including the following steps: Use the travel time gradient weights of multiple grid nodes to obtain the x at the current ray tracing position n The gradient value of: Among them, h is the calculation step length; Set the termination threshold of ray tracing. When ||x is satisfied n+1 -x n When ||<ε·h, ray tracing is terminated. ε is a constant and is set according to different accuracy requirements.

9. The method according to claim 1, characterized in that Establishing the three-dimensional initial inversion velocity field constrained by micro-logging information includes the following steps: The initial inversion velocity field is interpolated along the layers to obtain a three-dimensional initial inversion velocity field containing micro-logging information in the work area.

10. The method according to claim 1, characterized in that According to the travel time residual and the Frechet derivative matrix, a tomographic inversion equation based on apparent slowness constraint and regularization constraint is constructed, which includes the following steps: Construct the objective function of the apparent slowness constraint and the regularization constraint, Φ(m)=(1-ω)S1+ωS2+τS3 Among them, Φ(m) is the objective function, m is the slowness, ω is the weight factor between the average slowness mismatch and the apparent slowness mismatch, and S1 is represents the normalization of the true travel time d, and l represents the ray path length; represents the normalization of the forward travel time G(m), S2 is represents the apparent slowness, which is the gradient of the travel time and offset curve, D x is the differential operator of travel time with respect to distance, τ is the smoothing weight factor, S3 is Rm, and R is the regularization operator; The steady-state equation that minimizes the objective function is iteratively updated to obtain: m k+1 =m k +Δm k (k=1,2,3...,N,) Among them, for the kth slowness m k , A k is the average slowness kernel function, B k is the apparent slowness kernel function, ε k I is the variable damping term, Δm k is the iteration step size, and N represents the maximum number of slownesses.

11. A high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion device, characterized in that: include: The global travel time determination module is used to determine the three-dimensional global forward modeling first arrival wave propagation travel time through a multi-template upwind difference algorithm and a narrow-band extension strategy for undulating terrain; The tomographic inversion model building module is used to perform ray tracing based on the three-dimensional global forward modeling first arrival wave propagation travel time, and combine the forward modeling travel time of the inverted initial velocity field to obtain the first arrival wave propagation path and the Frechet differential matrix; establish a three-dimensional initial inversion velocity field constrained by micro-logging information; and construct a tomographic inversion equation based on apparent slowness constraints and regularization constraints based on the travel time residual and the Frechet differential matrix. The velocity solution determination module is used to solve the tomographic inversion equation, obtain the tomographic inversion slowness solution, and convert the slowness solution into a velocity solution.

12. A computer storage medium, characterized in that: The computer storage medium stores computer executable instructions, which, when executed, implement the high-precision three-dimensional first-arrival wave travel time tomographic velocity inversion method described in any one of claims 1 to 10.

Images