Triangular mesh tomography velocity inversion optimization method based on node data

Abstract

The application provides a kind of triangulation tomographic velocity inversion optimization method based on node data, comprising: step 1, obtaining imaging profile and angle domain common imaging gather;Step 2, picking horizon on the imaging profile;Step 3, triangulation section of velocity field;Step 4, residual curvature analysis is carried out to angle gather, and depth residual is obtained, which is converted into travel time information;Step 5, ray tracing is carried out to triangulation velocity field, and ray path is obtained;Step 6, using triangulation ray path and travel time information, tomographic inversion is carried out to obtain slowness update, and this iteration is completed;Step 7, Gaussian beam prestack depth migration is carried out to the updated velocity field, and angle gather and imaging profile are obtained.The triangulation tomographic velocity inversion optimization method based on node data reduces the ill-conditioned nature of inversion, improves the calculation efficiency, improves the accuracy of velocity inversion of node data in large offset and deep structure, and realizes velocity field modeling.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic data processing for oil and gas exploration, and particularly to a triangular mesh tomographic velocity inversion optimization method based on node data. BACKGROUND

[0002] At present, the main methods for solving the problem of seismic data velocity inversion include migration velocity analysis, tomographic velocity inversion and full waveform inversion method.

[0003] The conventional migration velocity analysis method obtains the corresponding residual error by imaging gather residual curvature analysis to update the velocity, and takes the best imaging gather flattening and migration imaging quality as the criterion. The full waveform inversion is to take the best fitting of model data and observed data as the criterion, and to realize velocity recovery through seismic full waveform information. The tomographic inversion method is realized in the imaging domain by using the curvature of common imaging point gather, and uses travel time or waveform residual to update the velocity. This method has high precision and is relatively easy to implement.

[0004] Theoretically, the waveform inversion has the highest precision, but there are still many limitations in practical application, such as source wavelet, initial velocity precision and other problems. In the case of difficulty in practical application of the waveform inversion method, the tomographic velocity inversion method has lower requirements for data than the waveform inversion, and the precision and efficiency can meet the production demand, so it is the most widely used velocity inversion method at present.

[0005] At present, the commonly used method for tomographic velocity inversion is rectangular grid travel time tomography. However, when modeling the node data velocity, due to the poor flexibility of rectangular grid, it is difficult to adjust the number of grids in the large offset area according to the signal-to-noise ratio, so that the velocity tomographic precision in the large offset area is low. Secondly, the target area of large offset data is deep structure, and the number of grids is large under rectangular grid division, which leads to strong inversion ill-condition and low calculation efficiency. The signal-to-noise ratio of data in large offset area is low, and the conventional rectangular grid travel time tomography method is difficult to pick up high-precision travel time information in prestack data, which leads to low inversion precision.

[0006] In the Chinese patent application with the application number CN201810855447.5, a kind of imaging field stereotomography velocity inversion method is involved, comprising: step 1, the initial velocity field of input is made prestack depth migration, and depth migration profile is obtained, and angle domain common imaging point gather is extracted;Step 2, in angle domain common imaging point gather, pick up the residual curvature and depth residual at the corresponding position, and convert into travel time residual;Step 3, calculate stereotomography data space;Step 4, ray tracing is carried out in the current velocity model, and the travel time information and inversion kernel function in the model are obtained;Step 5, using the obtained stereotomography kernel function and the residual of real data space and model data space, establish inversion equation set, and the update of model is calculated, and this iteration is completed.The imaging field stereotomography velocity inversion method can obtain more accurate stereotomography data space, improve the precision and stability of inversion, and has better utilization value in the practice of stereotomography inversion in the future.

[0007] In the Chinese patent application with the application number CN201710495784.3, a kind of microseismic positioning and tomography method is involved, comprising steps 1, read initial model and inversion parameter;Step 2, the shortest path algorithm based on interface element is used to calculate the theoretical travel time and ray path of microseismic event;Step 3, when the underground medium structure and microseismic event parameters are simultaneously inverted, the first order partial derivative of time about the to-be-inverted parameter is calculated;Step 4, the tomographic equation set of simultaneously inverting underground medium structure and microseismic event parameters is constructed;Step 5, the tomographic equation set is solved by conjugate gradient method, and the update of microseismic source parameter and velocity structure is obtained;Step 6, update and constrain velocity structure model and source parameter;Step 7, iteration termination judgment.The invention can simultaneously invert source parameter and velocity structure model using ground and well observation travel time data for horizontal layered medium, where source parameter refers to source location and time of origin, and velocity structure model refers to interval velocity and layer interface position.

[0008] In the Chinese patent application with the application number CN202011102625.0, a triangular net Fresnel zone travel time difference tomographic inversion method is involved. The tomographic inversion method used by the present application is not easily affected by environmental factors, has small working difficulty, and has low construction cost. Unlike the traditional tomographic inversion which uses a rectangular grid to subdivide the model, the present application can better fit the irregular bridge model at the edge, reducing the error of forward modeling, so that the subsequent inversion process can also achieve good results. The present application uses the trace closest to the shot point as the reference trace, calculates the difference between the theoretical first arrival time and the actual measured first arrival time, and uses this as the basis to calculate the true first arrival time of all the receiver points. Using this method, the first arrival time can be accurately picked up, which is beneficial to subsequent inversion. The Fresnel zone technology used by the present application avoids the problem of requiring a large number of geophone devices to ensure accuracy in traditional tomographic inversion, so that the tomographic inversion method can still achieve good bridge detection results in relatively small projects such as bridge detection.

[0009] In the Chinese patent application with the application number CN201210015474.4, a multi-scale regular grid tomographic inversion static correction method in geophysical exploration static correction is involved. Seismic data is collected, a rectangular grid is sub-divided and the model velocity is initialized. The ray path and travel time of each shot point and receiver point pair are obtained by forward modeling. The difference between the actual picked first arrival time and the ray forward travel time is calculated. Multi-scale tomographic inversion is performed to update the model velocity. The iteration is repeated until the inverted velocity field is stable, and the grid tomographic inversion static correction is completed. The present application uses a multi-scale model for sub-division, which can better invert the velocity value at the position where the ray is sparse, reduces human error, and improves the accuracy of velocity model inversion and the effect of static correction.

[0010] The above prior art has great differences from the present application and cannot solve the technical problems we want to solve. Therefore, we have invented a new triangular grid tomographic velocity inversion optimization method based on node data. SUMMARY

[0011] The purpose of the present application is to provide a triangular grid tomographic velocity inversion optimization method based on node data which utilizes the flexibility of triangular grid subdivision, adjusts the number of grids according to the signal-to-noise ratio at large offset distances and the complexity of deep structures, and reduces the ill-conditioned nature of inversion.

[0012] The purpose of the present application can be achieved by the following technical measures: a triangular grid tomographic velocity inversion optimization method based on node data, which includes:

[0013] Step 1, obtaining an imaging profile and an angle domain common imaging point gather;

[0014] Step 2, picking up horizons on the imaging profile;

[0015] Step 3, triangulation mesh is divided for the velocity field;

[0016] Step 4, residual curvature analysis is carried out on the angle gather to obtain a depth residual, which is converted into travel time information;

[0017] Step 5, ray tracing is carried out on the triangulation mesh velocity field to obtain a ray path;

[0018] Step 6, the triangulation mesh ray path and travel time information are used to carry out tomographic inversion to obtain a slowness update amount, and the iteration is completed;

[0019] Step 7, the updated velocity field is subjected to Gaussian beam prestack depth migration to obtain an angle gather and an imaging profile.

[0020] The object of the application can also be achieved through the following technical measures:

[0021] In step 1, based on the input rectangular mesh initial velocity field, prestack depth migration is carried out on node data to obtain an imaging profile and an angle domain common imaging point gather.

[0022] In step 1, the Gaussian beam prestack depth migration method is selected to obtain a Gaussian beam prestack depth migration imaging profile and an angle gather of the node data.

[0023] In step 1, first, the kinematic ray tracing method is used to obtain a Gaussian beam at a shot point and a receiver point, then forward and reverse wave field continuation is carried out on the shot point and the Gaussian beam center, respectively, imaging is carried out according to the cross-correlation imaging formula, and then in the wave field continuation process, the angle information of the incident wave field and the reflected wave field at different imaging points is obtained by using the travel time information at the center position of the Gaussian beam and the node information around the center position.

[0024] In step 1, the formula for imaging by the cross-correlation imaging formula is:

[0025]

[0026] Wherein, ω is an angular frequency, x, y and z represent spatial positions, subscript s represents a shot point, subscript r represents a receiver point, subscript-free x represents an imaging position, u is a continuation wave field, G * represents a Green function, and I pre is an imaging value.

[0027] In step 1, in the ray center coordinate system, the travel time at a point A outside the ray can be expressed by the travel time of a point B vertically intersecting the ray:

[0028]

[0029] where M(B) represents the second-order partial derivative of the travel time at point B, n is the normal distance from point A to the ray, i.e. the coordinate value in the ray center coordinate system with B as the origin, and T represents the travel time;

[0030] Taking the derivative of both sides of equation (2) with respect to x, we have:

[0031]

[0032] Let l x ,l z represent the component values of the tangential unit vector at point B in the x and z directions of the rectangular coordinate system, respectively; then we have:

[0033] n = (x-x B )l z -(z-z B )l x (equation 4)

[0034] Substituting the above equation into equation (3) and taking the derivative of both sides of the equation with respect to z, we can obtain the propagation angle θ of the Gaussian beam, i.e. the angle between the ray and the positive direction of the Z axis:

[0035]

[0036] where After obtaining the propagation angles of the Gaussian beams at the source point and the receiving point by using the above formula, the opening angle in the migration process can be obtained, and then the imaging values are arranged according to the angles to complete the extraction of the angle gathers.

[0037] In step 4, the angle gather is obtained according to the Gaussian beam prestack depth migration method by using the travel time tomography method based on the angle gather, and the conversion relationship between the depth residual and the travel time residual is established to provide the travel time information for the tomographic inversion equation.

[0038] In step 4, the conversion relationship between the depth residual and the travel time residual is:

[0039]

[0040] where Δz is the imaging point depth error caused by the inaccuracy of the parameter field; β is the dip angle of the stratum near the imaging point, θ is the seismic wave incidence angle, L1+L2 is the ray path change, and Δv represents the velocity disturbance between the initial model and the true model.

[0041] In step 5, the tomographic velocity inversion is performed by grid division on the geological model, and the travel time information and the ray path in the grid are obtained by using ray tracing. The ray tracing is performed on the velocity field of the triangular grid division by using the slowness square gradient variable.

[0042] In step 5, the Hamilton-Jacobi operator of the eikonal equation in the medium of slowness square is:

[0043]

[0044] In the above formula, p i is the slowness vector, is the function relationship of the n-th power of slowness with respect to the spatial position change; according to the eigenvalue method, the ray tracing equation can be obtained:

[0045]

[0046] In the above formula, p1 and p2 respectively represent the components of the ray parameter in the x and z directions of the spatial position; when the square of the slowness in the medium changes linearly with the spatial position, the slowness square expression can be rewritten as:

[0047]

[0048] In the above formula, s 00 is the square slowness of the reference point, s x ' and s z ' respectively represent the gradients of the square slowness in the x and z directions, and the slowness vector expression can be obtained by combining the above formula:

[0049]

[0050] By combining the above formula, the analytical solution of the ray coordinate and the ray travel time required for tomographic inversion can be obtained:

[0051]

[0052] Based on the above triangular grid ray tracing theory, triangular grid ray tracing is performed on the basis of triangular grid division of the velocity to obtain the travel time and ray path information.

[0053] In step 6, the Radon transform is the basic theory of travel time tomography, and its basic principle is to obtain the underground medium velocity information by back-projecting the seismic data along the ray direction to construct the underground parameter field; among them, the ray-based travel time tomography is to integrate the slowness in the travel time propagation direction to obtain the corresponding ray travel time:

[0054]

[0055] In the formula: s is the reciprocal of the velocity, and l represents the integral path, and the upper and lower limits of the integral are the receiver and the source; the model slowness change will cause the travel time to be disturbed, and the travel time residual vector Δt and the slowness update Δs are introduced:

[0056] LΔs=Δt (Formula 13)

[0057] The relationship between the ray travel time variation and the slowness disturbance can be established according to the above formula, wherein L represents the ray path length;

[0058] The travel time residual and the ray path are brought into the tomography equation, and the velocity update is obtained by solving, and the inversion velocity field is constructed by repeated iteration.

[0059] In step 7, the updated velocity field is subjected to Gaussian beam prestack depth migration to obtain an angle gather and an imaging profile, and according to the angle gather flattening degree, the accuracy of the picked horizon and the velocity precision requirement, it is judged whether to perform the next iteration, if the iteration is continued, it returns to step 1, and the process is repeated, otherwise, the loop is exited, and the final velocity field is obtained.

[0060] The triangular mesh tomography velocity inversion optimization method based on node data in the application is aimed at the problems existing in the traditional rectangular mesh tomography inversion when the node data velocity modeling is performed, and is based on the angle domain imaging gather (angle gather), and high-precision travel time information is obtained through the conversion relationship between the travel time residual and the depth residual. Meanwhile, the flexibility of the triangular mesh is utilized, the form and quantity of the mesh are adjusted according to the signal-to-noise ratio at large offset and the fluctuation degree of the deep structure, the ill-conditioned property of the inversion is reduced, the calculation efficiency is improved, the accuracy of the velocity inversion of the node data at large offset and deep structure is improved, and the velocity field modeling is realized. The triangular mesh tomography velocity inversion optimization method based on node data has advantages that other methods do not have, and the specific advantages and characteristics are shown in the following aspects:

[0061] (1) Flexibility of the method. The flexibility of the triangular mesh division is utilized, the quantity of the mesh is adjusted according to the signal-to-noise ratio of different offset ranges and the complexity of the geological structure, and the ill-conditioned property of the inversion is reduced.

[0062] (2) Efficiency of the method. The mesh boundary is coincided with the stratum boundary under the triangular mesh division, the stratum form can be inverted with high precision, the interface judgment calculation in the ray tracing process is avoided, and the calculation efficiency and the calculation precision are improved.

[0063] (3) Strong adaptability of the method. The triangular mesh ray tracing method has higher processing precision compared with the rectangular mesh ray tracing because the ray coverage times are higher and the coverage is more uniform in the large offset and deep area.

[0064] (4) Simple operation and easy implementation. The method has simple steps and is convenient, and does not need complicated parameter adjustment and test. BRIEF DESCRIPTION OF DRAWINGS

[0065] Figure 1 It is a schematic view of the Gaussian beam migration ray center coordinate system of the application;

[0066] Figure 2 A time-keeping residual conversion relationship diagram of the present application;

[0067] Figure 3 A flow chart of a specific embodiment of the present application of the triangular mesh tomography velocity inversion optimization method based on node data;

[0068] Figure 4 A model triangular mesh ray tracing diagram of a specific embodiment 1 of the present application;

[0069] Figure 5 A model rectangular mesh ray tracing diagram of a specific embodiment 1 of the present application;

[0070] Figure 6 A real velocity model diagram of a specific embodiment 2 of the present application;

[0071] Figure 7 A triangular mesh dissection diagram of an initial velocity model of a specific embodiment 2 of the present application;

[0072] Figure 8 A migrated profile diagram of an initial velocity field of a specific embodiment 2 of the present application;

[0073] Figure 9 An angle gather diagram of an initial velocity field of a specific embodiment 2 of the present application;

[0074] Figure 10 A triangular mesh tomography inversion velocity field diagram of a specific embodiment 2 of the present application;

[0075] Figure 11 An inversion migrated profile diagram of a specific embodiment 2 of the present application;

[0076] Figure 12 An inversion angle gather diagram of a specific embodiment 2 of the present application;

[0077] Figure 13 A rectangular mesh tomography velocity field diagram of a specific embodiment 3 of the present application;

[0078] Figure 14 A rectangular mesh tomography velocity migrated profile diagram of a specific embodiment 3 of the present application;

[0079] Figure 15 A rectangular mesh tomography velocity angle gather diagram of a specific embodiment 3 of the present application;

[0080] Figure 16 A triangular mesh tomography velocity field diagram of a specific embodiment 3 of the present application;

[0081] Figure 17This is a schematic diagram of the offset profile after triangular mesh tomography in a specific embodiment 3 of the present invention;

[0082] Figure 18 This is a schematic diagram of the angle gather after triangular mesh tomography, which is a specific embodiment 3 of the present invention. Detailed Implementation

[0083] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0084] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0085] Large offset data suffers from low signal-to-noise ratio and low travel time acquisition accuracy. To address this issue, this invention leverages the advantages of angle gathers—few artifacts and high resolution—in tomographic inversion, obtaining high-precision travel time information based on the conversion relationship between travel time residuals and depth residuals. Furthermore, it utilizes the advantage of triangular mesh ray tracing's uniform coverage in deep structures and at large offsets to achieve high-precision velocity modeling for large offset data.

[0086] like Figure 3 As shown, Figure 3 This is a flowchart of the node-data-based triangular mesh tomography velocity inversion optimization method of the present invention. The node-data-based triangular mesh tomography velocity inversion optimization method includes:

[0087] (101) Based on the initial velocity field of the input rectangular grid, perform pre-stack depth migration on the node data to obtain the imaging profile and the common imaging point gather (angle gather, ADCIG) in the angle domain.

[0088] This invention acquires Gaussian beam pre-stack depth migration imaging profiles and corner gathers from node data. Migration imaging is fundamental to tomographic velocity inversion, providing imaging results and imaging gathers to obtain stratigraphic constraints and travel time residuals. This invention selects the Gaussian beam pre-stack depth migration method for imaging and acquiring corner gathers. The Gaussian beam migration method has good imaging performance for initial velocity fields with large velocity errors, and the acquired corner gathers are more accurate than other types of imaging gathers.

[0089] Firstly, the Gauss beam at the shot point and receiver point is obtained by using the kinematic ray tracing method, then the forward and backward continuation of wave field is carried out at the shot point and the beam center respectively, and the imaging is carried out according to the cross-correlation imaging formula:

[0090]

[0091] Where ω is the angular frequency, x, y, z represent the spatial position, subscript s represents the source point, subscript r represents the receiving point, x without subscript represents the imaging position, u is the continuation wave field, G * represents the Green function, I pre is the imaging value.

[0092] Then, the angle information of the incident wave field and the reflected wave field at different imaging points can be obtained by using the travel time information at the beam center position and the node information around the center position to calculate the angle information of the incident wave field and the reflected wave field at different imaging points during the wave field continuation.

[0093] In the ray center coordinate system ( Figure 1 ), the travel time at a point A outside the ray can be expressed by the travel time of the point B which is perpendicular to the ray:

[0094]

[0095] Where M(B) represents the second-order partial derivative of the travel time at point B, n is the normal distance of point A from the ray, i.e. the coordinate value in the ray center coordinate system with B as the origin, and T represents the travel time.

[0096] Taking the derivative of both sides of equation (2) with respect to x gives:

[0097]

[0098] Let l x ,l z represent the component values of the tangential unit vector at point B in the x and z directions of the rectangular coordinate system respectively. Then we have:

[0099] n=(x-x B )l z -(z-z B )l x (Formula 4)

[0100] Substituting the above formula into equation (3) and taking the derivative of both sides with respect to z, the propagation angle θ of the Gauss beam (the angle between the ray and the positive direction of the Z axis) is obtained as:

[0101]

[0102] Where After the propagation angles of the Gaussian beams at the source point and the receiving point are obtained by using the above formula, the opening angle in the migration process can be solved, and then the imaging values are arranged according to the angles to complete the extraction of the angle gathers.

[0103] (102) picking up the horizon on the imaging profile.

[0104] (103) under the constraint of the horizon, performing triangular mesh division on the initial velocity field to generate an initial triangular mesh velocity field.

[0105] (104) performing residual curvature analysis on the angle gather to obtain a depth residual, which is converted into travel time information.

[0106] The application adopts the travel time tomography method based on the angle gather, acquires the angle gather according to the above Gaussian beam prestack depth migration method, and provides travel time information for the tomographic inversion equation by establishing the conversion relationship between the depth residual and the travel time residual, as shown in the following formula: Figure 2

[0107] In the figure, Δz is the depth error of the imaging point caused by the inaccuracy of the parameter field. β is the dip angle of the stratum near the point, θ is the seismic wave incidence angle, and L1+L2 is the ray path change. Then, the conversion relationship between the depth residual and the travel time residual can be obtained from the figure as follows:

[0108]

[0109] Where Δv represents the velocity disturbance between the initial model and the real model.

[0110] (105) performing triangular mesh ray tracing on the initial triangular mesh velocity field to solve the ray path corresponding to the angle gather.

[0111] The tomographic velocity inversion performs mesh division on the geological model, and uses ray tracing to solve the travel time information and the ray path in the mesh. For the velocity field of the triangular mesh division, the ray tracing is performed by using the slowness square gradient variable, the expression structure is simple and accurate, there is an analytical solution in a single mesh, and the accuracy is higher than that of the traditional single-step method. The Hamilton-Jacobi operator of the eikonal equation in the slowness square medium is:

[0112]

[0113] In the above formula, p i is the slowness vector, is a function relationship of the n-th power of the slowness with respect to the spatial position. According to the eigenvalue method, the ray tracing equation can be obtained as follows:

[0114]

[0115] ​In the above formula, p1, p2 represent the components of the ray parameter in the x, z direction of the spatial position. When the square of the slowness in the medium varies linearly with the spatial position, the square of the slowness can be rewritten as:

[0116]

[0117] In the formula, s 00 is the square slowness of the reference point, s x ' and s z ' represent the gradients of the square slowness in the x and z directions, and the expression of the slowness vector can be obtained by combining the above formula:

[0118]

[0119] By solving the above formula, the analytical solution of the ray coordinates and the ray travel time required for tomographic inversion can be obtained as:

[0120]

[0121] Based on the above triangular mesh ray tracing theory, triangular mesh ray tracing is performed on the basis of triangular mesh division of the velocity to obtain the travel time and ray path information.

[0122] (106)Using the triangular mesh ray path and travel time information, the slowness update is obtained by tomographic inversion to complete this iteration.

[0123] Radon transformation is the basic theory of travel time tomography, and its basic principle is to obtain the underground medium velocity information by back-projecting the seismic data along the ray direction to construct the underground parameter field. Among them, the ray-based travel time tomography is to integrate the slowness in the travel direction to obtain the corresponding ray travel time:

[0124]

[0125] In the formula: s is the reciprocal of the velocity, and l represents the integral path, and the upper and lower limits of the integral are the receiver and the source. The model slowness change will cause the travel time to be disturbed, and the travel time residual vector Δt and the slowness update Δs are introduced:

[0126] LΔs = Δt (Formula 13)

[0127] According to the above formula, the relationship between the ray travel time change and the slowness disturbance can be established. In the formula, L represents the length of the ray path.

[0128] The travel time residual and the ray path are brought into the tomographic equation, and the velocity update is obtained by solving, and the inversion velocity field is constructed by repeated iteration.

[0129] (107)The updated velocity field is subjected to Gaussian beam pre-stack depth migration to obtain an angle gather and an imaging profile, and whether to perform the next iteration is judged according to the angle gather flattening degree, the accuracy of the picked horizon and the velocity accuracy requirement, such as returning to step 101 to repeat the process if the iteration is continued, or otherwise jumping out of the loop to obtain the final velocity field.

[0130] The following are several specific embodiments of the application

[0131] Embodiment 1

[0132] The triangular mesh and rectangular mesh ray tracing designed by the application are compared and analyzed by taking a complex equivalent model as an example.

[0133] The triangular mesh and rectangular mesh are respectively divided by using an equal complexity parameter field, and 45 rays are uniformly shot downward at a depth of 4 km on the surface with a ray interval of 2. Figure 4 Figure 5 ):

[0134] In terms of the accuracy of boundary description, the arc interface of the triangular mesh is completely consistent with the shape of the stratum interface, and compared with the parameter field rectangular mesh division method, the triangular mesh has higher processing accuracy for stratum interfaces with large relief.

[0135] In terms of interface reflection and transmission interface processing accuracy, the grid boundary in the triangular mesh division method is the parameter field interface, while the parameter interface in the rectangular mesh division method is usually inside the grid, so the calculation of whether there is an interface in the grid can be omitted in the ray tracing process, improving the operation efficiency. In addition, the triangular mesh ray tracing eliminates the smoothing preprocessing required due to the velocity jump in the rectangular mesh division, reducing the processing error. The triangular mesh ray tracing can specify the reflection surface of the ray, and the detailed description of the specific horizon in the tomography process is more accurate.

[0136] According to the relief complexity of the velocity interface, the corresponding triangular mesh control points are selected, and the number of mesh division of different structural forms is adjusted. The selection of control points in the application introduces a stratum dip angle constraint, increases the number of grids in complex relief structure areas, and reduces the number of grids in relatively flat and simple structural forms. Compared with the rectangular mesh division method, the triangular mesh division method can reduce the number of mesh division, reduce the number of traversed grids in the ray tracing process, effectively reduce the calculation amount and storage amount, and correspondingly improve the forward and inverse accuracy.

[0137] Therefore, the triangular mesh ray tracing of the application is superior to the conventional rectangular mesh ray tracing.

[0138] Embodiment 2

[0139] The application takes an analog data embodiment, applies the method to data for angle gather-based triangular mesh velocity modeling to verify the effect of the method, and a specific flow chart is shown in the following Figure 3 . The steps of this embodiment are as shown in the inversion strategy steps above:

[0140] The analog data size is 801*450km( Figure 6 ), the true velocity field is transversely and longitudinally smoothed as an initial velocity field( Figure 7 ), and the initial velocity field is used for migration imaging. It can be known from Figure 8 that, due to the parameter field error, the diffraction wave does not converge, the horizon false image appears, the resolution is reduced, the angle gather is obviously curved upwards( Figure 9 ), and other problems.

[0141] The initial velocity field is subjected to triangular mesh travel time tomographic inversion to obtain a tomographic velocity field( Figure 10 ) and perform migration imaging. Compared with the initial parameter field migration profile, the inversion parameter field migration profile( Figure 11 ) has improved resolution, enhanced deep energy, completely suppressed diffraction wave, and completely flattened angle gather( Figure 12 ), and finally a high-precision inversion velocity field is constructed. Compared with the rectangular mesh, for deep complex structures, the triangular mesh inversion result is more consistent with the true model, and the sawtooth and burr problems caused by mesh partitioning are improved, and the inversion result is closer to the true stratum shape. It is proved that the strategy proposed in the application is adaptive to large offset data.

[0142] Embodiment 3

[0143] In order to further verify the effectiveness and processing effect of the method of the application, actual work area data are subjected to rectangular mesh travel time tomography and triangular mesh velocity tomography based on node data proposed in the application. The steps of this embodiment are as shown in the inversion strategy steps above:

[0144] The root mean square velocity field is obtained by using the prestack velocity analysis method, the initial velocity field is converted according to the Dix formula, the conventional rectangular mesh travel time tomography is used for processing, the rectangular mesh tomographic velocity field( Figure 13 ) is obtained, and prestack depth migration is performed, the tomographic migration profile( Figure 14 ) and the angle gather( Figure 15 ) are obtained. The initial velocity field is subjected to triangular mesh partitioning, and the triangular mesh travel time tomographic inversion process proposed in the application is used to obtain the triangular mesh tomographic inversion velocity field( Figure 16 ), prestack depth migration is performed, the triangular mesh tomographic migration profile( Figure 17 ) and the angle gather( Figure 18By comparing the migration profile and the angle gather of the two kinds of tomography methods, it can be concluded that the velocity field constructed by the triangular grid velocity tomography method has higher imaging resolution and better horizon continuity at large migration distance, and the angle gather is completely flattened, and the imaging effect of the tomographic velocity field of the present application is superior to the conventional travel time tomography method in resolution, phase axis continuity and clear depiction of structural form.

[0145] The actual data test results show that the method has good effect on the velocity inversion of the exploration area of node data, has strong adaptability to large offset low signal-to-noise ratio data, can accurately invert the deep structural form, and obtain more accurate parameter field and high-precision imaging results.

[0146] Finally, it should be pointed out that the above description is only the preferred embodiment of the present application and is not intended to limit the present application, although the present application has been described in detail with reference to the foregoing embodiments, and for those skilled in the art, the technical solutions recorded in the foregoing embodiments can still be modified, or some technical features can be replaced by equivalents. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the present application shall be included in the protection scope of the present application.

[0147] In addition to the technical features described in the specification, they are known to those skilled in the art.

Claims

1. A method for optimizing velocity inversion of tomography based on node data of triangular mesh, characterized in that, The node data-based triangular mesh tomographic velocity inversion optimization method comprises: Step 1, obtaining an imaging profile and an angle domain common imaging point gather; Step 2, picking up a horizon on the imaging profile; Step 3, performing triangular mesh dissection on a velocity field; Step 4, performing residual curvature analysis on an angle gather to obtain a depth residual and convert the depth residual into travel time information; Step 5, performing ray tracing on the triangular mesh velocity field to obtain a ray path; Step 6, using the triangular mesh ray path and the travel time information to perform tomographic inversion to obtain a slowness update amount and complete this iteration; Step 7, performing Gaussian beam prestack depth migration on the updated velocity field to obtain an angle gather and an imaging profile; In step 1, based on an input rectangular mesh initial velocity field, prestack depth migration is performed on node data to obtain an imaging profile and an angle domain common imaging point gather; a Gaussian beam prestack depth migration method is selected to obtain a Gaussian beam prestack depth migration imaging profile and an angle gather of the node data; first, Gaussian beams at shot points and geophones are obtained by using a kinematic ray tracing method, then forward and reverse wave field continuation is performed on the shot points and the Gaussian beam centers, respectively, imaging is performed according to a cross-correlation imaging formula, and then in the wave field continuation process, angle information of incident wave fields and reflected wave fields at different imaging points is obtained by using travel time information at the center positions of the Gaussian beams and node information around the center positions to obtain an angle gather; In step 4, a travel time tomography method based on the angle gather is adopted, the angle gather is obtained according to the Gaussian beam prestack depth migration method, and a conversion relationship between a depth residual and a travel time residual is established to provide travel time information for a tomographic inversion equation.

2. The node data based triangular mesh tomographic velocity inversion optimization method of claim 1, wherein, In step 1, a formula for imaging by the cross-correlation imaging formula is: where ω is the angular frequency, x, y, z represent spatial locations, the subscript sou represents the source point, and the subscript r represents the receiver point, represents the imaging location, u is the seismic wavefield, and G * represents the Green's function, and I pre is the imaging value.

3. The node data based triangular mesh tomographic velocity inversion optimization method of claim 2, wherein, In step 1, in a ray center coordinate system, the travel time at a point A outside a ray is expressed by using a point B vertically intersecting the ray: In the formula, M(B) represents a second-order partial derivative of the travel time at the point B, n is a normal distance of the point A from the ray, that is, a coordinate value in the ray center coordinate system with the point B as the origin, and t represents the travel time; the derivative of both sides of formula (2) with respect to x is obtained as follows: Let l x , z respectively represent the component values of the tangential unit vector at point B in the x and z directions of the rectangular coordinate system; then we have: n = (x - x B ) / (x - x z ) (Formula 3) B ) / (z - z x ) (Formula 4) The above formula is substituted into formula (3) and the derivative of both sides of the equation with respect to z is obtained, and the propagation angle θ of the Gaussian beam, that is, the angle between the ray and the positive direction of the Z axis, is obtained as follows: wherein After the propagation angles of Gaussian beams at the source point and the receiver point are obtained by using the above formula, the opening angle in the migration process is calculated, and then the imaging values are arranged according to the angles to complete the extraction of the angle gathers.

4. The node data based triangular mesh tomographic velocity inversion optimization method of claim 3, wherein, In step 4, the conversion relationship between the depth residual and the travel time residual is: In the formula, Δz is an imaging point depth error caused by inaccuracy of a parameter field; β is a dip angle of a stratum near the imaging point; Δv represents a velocity disturbance between an initial model and a real model; and Δt represents a travel time residual.

5. The node data based triangular mesh tomographic velocity inversion optimization method of claim 1, wherein, In step 5, tomographic velocity inversion is performed by dissectioning a geological model in a mesh and using ray tracing to obtain travel time information and a ray path in the mesh, and for a triangular mesh dissection velocity field, ray tracing is performed by using a slowness square gradient variable.

6. The node data based triangular mesh tomographic velocity inversion optimization method of claim 5, wherein, In step 5, a Hamilton-Jacobi operator of an eikonal equation in a slowness square medium is: p in the above formula i is the slowness vector, is the function relationship of the k-th power of the slowness with respect to the spatial position change; according to the eigenvalue method, the ray tracing equation is obtained: In the formula, when the square of the slowness in the medium varies linearly with a spatial position, the slowness square expression is rewritten as: where s 00 is the square slowness of the reference point, s x and s z ' represent the gradients of the square slowness in the x and z directions, respectively, and the slowness vector expression is obtained by combining the above equations: Analytical solutions of ray coordinates and ray travel times required by tomographic inversion are obtained by solving the above formula as follows: Based on the above-mentioned triangular mesh ray tracing theory, the triangular mesh ray tracing is carried out on the basis of the triangular mesh division of the velocity to obtain the travel time and ray path information.

7. The node data based triangular mesh tomographic velocity inversion optimization method of claim 1, wherein, In step 6, the Radon transform is the basic theory of travel time tomography, and the basic principle is to obtain the underground medium velocity information by back-projecting the seismic data along the ray direction to construct the underground parameter field; wherein the ray-based travel time tomography is to integrate the slowness in the travel time propagation direction to obtain the corresponding ray travel time: In the formula: s m Let l be the reciprocal of velocity, i.e., the slowness, and l represent the integration path. The upper and lower limits of integration are the receiver point r and the source point sou, respectively. Changes in the model's slowness will cause disturbances during travel, so the slowness update Δs m have: LΔs m = Δt (Equation 13) According to the above formula, the relationship between the ray travel time variation and the slowness perturbation is established; wherein L represents the ray path length; The travel time residual and the ray path are brought into the tomography equation, and the velocity update is obtained by solving, and the inversion velocity field is constructed by repeated iteration.

8. The node-data-based triangular mesh tomographic velocity inversion optimization method of claim 1, wherein, In step 7, the updated velocity field is subjected to Gaussian beam prestack depth migration to obtain the angle gather and the imaging profile, and according to the angle gather flattening degree, the accuracy of the picked horizon and the velocity precision requirement, it is judged whether to continue the next iteration, if the iteration is continued, it returns to step 1 and the process is repeated; otherwise, the loop is exited and the final velocity field is obtained.

Images