Triangular mesh chromatography velocity inversion optimization method based on node data

Abstract

The invention provides a triangular mesh tomography velocity inversion optimization method based on node data. The triangular mesh tomography velocity inversion optimization method comprises the following steps: step 1, acquiring an imaging profile and angle domain common imaging point gather; step 2, picking up a horizon on the imaging profile; 3, performing triangular mesh generation on the velocity field; 4, performing residual curvature analysis on the angle gather to obtain a depth residual error, and converting the depth residual error into travel time information; 5, performing ray tracing on the triangular mesh velocity field, and obtaining a ray path; step 6, utilizing the triangular mesh ray path and travel time information to carry out tomographic inversion to obtain slowness update quantity, and completing the iteration; and 7, performing Gaussian beam pre-stack depth migration on the updated velocity field to obtain an angle gather and an imaging profile. According to the triangular mesh tomography velocity inversion optimization method based on the node data, the ill-conditioned property of inversion is reduced, the calculation efficiency is improved, the precision of velocity inversion constructed by the node data in a large offset distance and a deep layer is improved, and velocity field modeling is realized.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration seismic data processing technology, and in particular to a triangular grid tomography velocity inversion optimization method based on node data. Background Art

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

[0003] Conventional migration velocity analysis method obtains the corresponding residual error through residual curvature analysis of imaging gathers to update velocity, and takes imaging gather flattening and the best migration imaging quality as the judgment criteria. Full waveform inversion uses the best fit between model data and observation data as the judgment criteria, and realizes velocity recovery through seismic full waveform information. The tomographic inversion method uses the curvature of the common imaging point gather in the imaging domain, and uses travel time or waveform residual to update velocity. This method has high accuracy and is relatively easy to implement.

[0004] Theoretically, waveform inversion has the highest accuracy, but there are still many limitations in its practical application, such as source wavelets, initial velocity accuracy, etc. In the case of the difficulty in the practical application of waveform inversion methods, the tomographic velocity inversion method does not have as high data requirements as waveform inversion. Its accuracy and efficiency can meet production needs, and it is currently the most widely used velocity inversion method.

[0005] At present, the commonly used method for tomographic velocity inversion is rectangular grid travel time tomography. However, when modeling the velocity of node data, due to the poor flexibility of rectangular grids, it is difficult to adjust the number of grids according to the signal-to-noise ratio in the large offset area, resulting in low velocity tomography accuracy at large offsets; secondly, the target area of ​​large offset data is a deep structure, and the number of grids under rectangular grid division is large, resulting in strong inversion pathology and reduced calculation efficiency; the signal-to-noise ratio of data at large offsets is low, and conventional rectangular grid travel time tomography methods are difficult to pick up high-precision travel time information in pre-stack data, resulting in low inversion accuracy.

[0006] In the Chinese patent application with application number: CN201810855447.5, a method for inversion of stereo tomography velocity in the imaging domain is involved, including: step 1, pre-stack depth migration of the input initial velocity field to obtain the depth migration profile, and extract the angle domain common imaging point gather; step 2, in the angle domain common imaging point gather, pick up the residual curvature and depth residual at the corresponding position, and convert them into travel time residual; step 3, calculate the stereo tomography data space; step 4, perform ray tracing in the current velocity model to obtain the travel time information and inversion kernel function in the model; step 5, use the obtained stereo tomography kernel function and the residual between the real data space and the data space in the model to establish the inversion equation group, and calculate the update amount of the model to complete this iteration. This imaging domain stereo tomography velocity inversion method can obtain a more accurate stereo tomography data space, improve the accuracy and stability of the inversion, and has better utilization value in the future practice of stereo tomography inversion.

[0007] In the Chinese patent application with application number: CN201710495784.3, a microseismic positioning and tomography method is involved, including step 1, reading in the initial model and inversion parameters; step 2, using the shortest path algorithm based on the interface element to calculate the theoretical travel time and ray path of the microseismic event; step 3, calculating the first-order partial derivative of the parameter to be inverted when the underground medium structure and the microseismic event parameters are inverted simultaneously; step 4, constructing a tomographic equation group for simultaneously inverting the underground medium structure and the microseismic event parameters; step 5, solving the tomographic equation group by the conjugate gradient method to obtain the updated amount of the microseismic source parameters and the velocity structure; step 6, updating and constraining the velocity structure model and the source parameters; step 7, iterative termination judgment. This invention is aimed at horizontal layered media, and can simultaneously invert the source parameters and the velocity structure model by using the travel time data observed on the ground and in the well, where the source parameters refer to the source position and the time of earthquake occurrence, and the velocity structure model refers to the layer velocity and the layer interface position.

[0008] In the Chinese patent application with application number: CN202011102625.0, a triangulated Fresnel zone time difference tomographic inversion method is involved. The tomographic inversion method used in the present invention is not easily affected by environmental factors, has low working difficulty and low construction cost; and unlike the traditional tomographic inversion that uses rectangular grids to divide the model, the present invention can better fit the bridge model with irregular edges, reduce the error of forward modeling, and enable the subsequent inversion process to achieve good results. The invention uses the track closest to the shot point as the reference track, calculates the difference between its theoretical first arrival time and the actual measured first arrival time, and uses this as a reference to deduce the true first arrival time of all detection points. This method can be used to pick up the first arrival travel time more accurately, which is conducive to subsequent inversion. The Fresnel zone technology used in the invention avoids the problem that traditional tomographic inversion requires a large number of detection equipment to ensure accuracy, 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 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, rectangular grid division is performed and model velocity is initialized, the ray path and travel time of each excitation point and receiving point pair are obtained by forward modeling, the difference between the actual first arrival time picked and the ray forward modeling travel time is calculated, multi-scale tomographic inversion is performed, the model velocity is updated, and it is repeated until the inverted velocity field is stable, and the grid tomographic inversion static correction is completed. The invention uses a multi-scale model for division to better invert the velocity value at the sparse position of the ray, reduce human errors, and improve the forward and inversion accuracy of the velocity model and the static correction effect.

[0010] The above existing technologies are all significantly different from the present invention and fail to solve the technical problem we want to solve. For this reason, we have invented a new triangular grid tomography velocity inversion optimization method based on node data. Summary of the invention

[0011] The purpose of the present invention is to provide a triangular mesh tomography velocity inversion optimization method based on node data which utilizes the flexibility of triangular mesh subdivision, adjusts the number of meshes according to the signal-to-noise ratio at large offsets and the complexity of deep structures, and reduces inversion pathology; and optimizes the triangular mesh ray tracing method.

[0012] The object of the present invention can be achieved by the following technical measures: a triangular grid tomography velocity inversion optimization method based on node data, the triangular grid tomography velocity inversion optimization method based on node data comprises:

[0013] Step 1, obtaining imaging sections and angle domain common imaging point gathers;

[0014] Step 2, picking the layer on the imaging section;

[0015] Step 3, triangulate the velocity field;

[0016] Step 4, perform residual curvature analysis on the diagonal gathers to obtain depth residuals and convert them into travel time information;

[0017] Step 5, ray tracing is performed on the velocity field of the triangular mesh to obtain the ray path;

[0018] Step 6: Using the triangular mesh ray path and travel time information, perform tomographic inversion to obtain the slowness update amount and complete this iteration;

[0019] Step 7: Perform Gaussian beam prestack depth migration on the updated velocity field to obtain angle gathers and imaging profiles.

[0020] The purpose of the present invention can also be achieved by the following technical measures:

[0021] In step 1, based on the input rectangular grid initial velocity field, the node data is prestacked and depth migrated to obtain the imaging section and the angle domain common imaging point gather.

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

[0023] In step 1, the kinematic ray tracing method is first used to obtain the Gaussian beam at the shot point and the detection point, and then the wave field of the shot point and the center of the Gaussian beam is forwardly extended and reversely extended respectively, and imaging is performed according to the cross-correlation imaging formula. Then, in the process of wave field extension, the travel time information at the center of the Gaussian beam and the node information around the center are used to calculate the angle information of the incident wave field and the reflected wave field at different imaging points to obtain the angle gather.

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

[0025]

[0026] Where ω is the angular frequency, x, y, and z represent spatial positions, subscript s represents the source point, subscript r represents the receiving point, x without a subscript represents the imaging position, u is the extended wave field, and G * represents the Green's function, I pre is the 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 at point B that intersects the ray perpendicularly through the point:

[0028]

[0029] 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, that is, the coordinate value in the ray center coordinate system with point 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 get:

[0031]

[0032] Let l x ,l z Respectively represent the component values ​​of the tangent unit vector at point B in the x and z directions in the rectangular coordinate system; then:

[0033] n=(xx B ) z -(zz B ) x (Formula 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, that is, the angle between the ray and the positive direction of the Z axis:

[0035]

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

[0037] In step 4, the travel time tomography method based on angle gathers is adopted to obtain angle gathers according to the Gaussian beam prestack depth migration method, and the travel time information is provided for the tomographic inversion equation by establishing a conversion relationship between the depth residual and the travel time residual.

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

[0039]

[0040] Among them, Δz is the depth error of the imaging point caused by inaccurate parameter field; β is the inclination of the stratum near the imaging point, θ is the incident angle of the seismic wave, L1+L2 is the change in the ray path, 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 meshing the geological model and using ray tracing to obtain the travel time information and ray paths within the grid. The slowness square gradient variable is used to perform ray tracing on the velocity field of the triangular mesh.

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

[0043]

[0044] In the above formula, p i is the slowness vector, is the functional relationship between the nth power of slowness and the change of spatial position; according to the eigenvalue method, the ray tracing equation can be obtained:

[0045]

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

[0047]

[0048] In the formula, s 00 is the squared slowness of the reference point, s x ′ and s z ′ represents the gradient of the square slowness in the x and z directions respectively. Combining the above formula, the slowness vector expression can be obtained:

[0049]

[0050] Combining the above formulas for solution, the ray coordinates and ray travel time analytical solutions required for tomographic inversion can be obtained as follows:

[0051]

[0052] Using 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.

[0053] In step 6, Radon variation is the basic theory of travel time tomography. Its basic principle is to obtain underground medium velocity information by back-projecting seismic data along the ray direction to construct underground parameter fields. Among them, ray travel time tomography integrates the slowness in the travel time propagation direction to obtain the corresponding ray travel time:

[0054]

[0055] Where: s is the inverse of velocity, l represents the integral path, the upper and lower limits of the integral are the detection point and the source respectively; the change of model slowness will cause disturbance in travel time, and the introduction of travel time residual vector Δt and slowness update Δs is:

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

[0057] According to the above formula, the relationship between the change in ray travel time and the slowness disturbance can be established; where L represents the length of the ray path;

[0058] The travel time residual and ray path are introduced into the tomographic equation, the velocity update is obtained by solving it, and the inversion velocity field is constructed by repeated iterations.

[0059] In step 7, Gaussian beam prestack depth migration is performed on the updated velocity field to obtain angle gathers and imaging profiles. Based on the degree of flattening of the angle gathers, the accuracy of the picked layers, and the velocity accuracy requirements, it is determined whether to perform the next iteration. If the iteration is continued, return to step 1 and repeat the process; otherwise, exit this loop to obtain the final velocity field.

[0060] The triangular grid tomography velocity inversion optimization method based on node data in the present invention aims at the problems existing in the traditional rectangular grid tomography inversion when modeling the node data velocity. Based on the angle domain imaging gather (angle gather), high-precision travel time information is obtained through the conversion relationship between travel time residual and depth residual. At the same time, the flexibility of the triangular grid is utilized to adjust the shape and number of the grid according to the signal-to-noise ratio at large offset and the degree of undulation of deep structure, reduce the pathological nature of the inversion, improve the computational efficiency, improve the accuracy of the velocity inversion of the node data at large offset and deep structure, and realize the modeling of the velocity field. The triangular grid tomography velocity inversion optimization method based on node data has advantages that other methods do not have, and its specific advantages and characteristics are manifested in the following aspects:

[0061] (1) Flexibility of the method. By utilizing the flexibility of triangular mesh generation, the number of meshes can be adjusted according to the signal-to-noise ratio of different offset ranges and the complexity of geological structures to reduce inversion pathology.

[0062] (2) High efficiency of the method. The mesh boundary and the stratum boundary coincide with each other under triangular meshing, which can invert the stratum morphology with high precision, avoid the interface judgment calculation in the ray tracing process, and improve the calculation efficiency and accuracy.

[0063] (3) Strong adaptability of the method. The triangular grid ray tracing method has a higher number of ray coverage times in large offset and deep areas, and the coverage is more uniform. Compared with rectangular grid ray tracing, the processing accuracy is higher. Therefore, the present invention has stronger adaptability to terrain structure in seismic exploration.

[0064] (4) The operation is simple and easy to implement. The method is simple and convenient, and does not require tedious parameter adjustment tests. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 A schematic diagram of the center coordinate system of the Gaussian beam offset ray of the present invention;

[0066] Figure 2 It is a schematic diagram of the travel time residual conversion relationship of the present invention;

[0067] Figure 3 It is a flow chart of a specific embodiment of the triangular grid tomography velocity inversion optimization method based on node data of the present invention;

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

[0069] Figure 5 A schematic diagram of ray tracing of a rectangular grid of a model according to a specific embodiment 1 of the present invention;

[0070] Figure 6 is a schematic diagram of a true velocity model of a specific embodiment 2 of the present invention;

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

[0072] Figure 8 is a schematic diagram of an initial velocity field migration profile according to a specific embodiment 2 of the present invention;

[0073] Fig. 9 is a schematic diagram of an angle gather of an initial velocity field according to a specific embodiment 2 of the present invention;

[0074] Fig.10 It is a schematic diagram of a velocity field inverted by triangular grid tomography according to a specific embodiment 2 of the present invention;

[0075] Fig.11 is a schematic diagram of an inversion migration profile of a specific embodiment 2 of the present invention;

[0076] Fig.12 is a schematic diagram of an inversion angle gather according to a specific embodiment 2 of the present invention;

[0077] Fig.13 is a schematic diagram of a rectangular grid tomographic velocity field according to a specific embodiment 3 of the present invention;

[0078] Fig.14 It is a schematic diagram of a rectangular grid tomography velocity migration profile according to a specific embodiment 3 of the present invention;

[0079] Fig.15 It is a schematic diagram of a rectangular grid tomography velocity angle gather according to a specific embodiment 3 of the present invention;

[0080] Fig.16 is a schematic diagram of a triangular grid tomographic velocity field according to a specific embodiment 3 of the present invention;

[0081] Fig.17It is a schematic diagram of a cross section after the migration of the triangular grid tomography according to a specific embodiment 3 of the present invention;

[0082] Fig.18 It is a schematic diagram of the angle gather after triangular grid tomography according to a specific embodiment 3 of the present invention. DETAILED DESCRIPTION

[0083] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used in the present invention have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

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

[0085] The signal-to-noise ratio of large offset data is low, and the travel time picking accuracy is low. To address this problem, the present invention utilizes the advantages of fewer false images and high resolution of angle gathers in tomographic inversion, and obtains high-precision travel time information based on the conversion relationship between travel time residuals and depth residuals. The advantages of uniform coverage of triangular mesh ray tracing in deep structures and large offsets are utilized to achieve high-precision velocity modeling of large offset data.

[0086] like Figure 3 As shown, Figure 3 The flowchart of the triangular grid tomography velocity inversion optimization method based on node data of the present invention is as follows. The triangular grid tomography velocity inversion optimization method based on node data includes:

[0087] (101) Based on the input rectangular grid initial velocity field, the node data is prestacked and depth migrated to obtain the imaging section and the angle domain common imaging point gather (angular gather, ADCIG).

[0088] Obtain Gaussian beam prestack depth migration imaging profiles and angle gathers of node data. Migration imaging is the basis of tomographic velocity inversion, providing imaging results and imaging gathers for tomographic inversion to obtain layer constraints and travel time residuals. The present invention selects the Gaussian beam prestack depth migration method for imaging and obtaining angle gathers. The Gaussian beam migration method has a better imaging effect for the initial velocity field with a large velocity error, and the obtained angle gathers are more accurate than other types of imaging gathers.

[0089] First, the kinematic ray tracing method is used to obtain the Gaussian beam at the shot point and the receiver point. Then, the wave field of the shot point and the beam center is extended forward and backward respectively, and imaging is performed according to the cross-correlation imaging formula:

[0090]

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

[0092] Then, in the process of wave field continuation, the travel time information at the center of the Gaussian beam and the node information around the center can be used to calculate the angle information of the incident wave field and the reflected wave field at different imaging points to obtain the angle gather.

[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 at point B that intersects the ray perpendicularly through the point:

[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, that is, the coordinate value in the ray center coordinate system with point B as the origin, and T represents the travel time.

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

[0097]

[0098] Let l x ,l z Respectively represent the x and z components of the tangent unit vector at point B in the rectangular coordinate system. Then:

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

[0100] 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 (the angle between the ray and the positive direction of the Z axis):

[0101]

[0102] in After obtaining the propagation angles of the Gaussian beam at the source point and the receiving point using the above formula, the angle during the migration process can be obtained, and then the imaging values ​​are arranged according to the angle to complete the extraction of the angle gather.

[0103] (102) Pick a layer on the imaging section.

[0104] (103) Under the layer constraint, the initial velocity field is triangulated to generate an initial triangulated mesh velocity field.

[0105] (104) Perform residual curvature analysis on the diagonal gathers to obtain depth residuals, which are converted into travel time information.

[0106] The present invention adopts a travel time tomography method based on angle gathers. According to the Gaussian beam prestack depth migration method, angle gathers are obtained, and travel time information is provided for the tomography inversion equation by establishing a conversion relationship between depth residuals and travel time residuals, such as Figure 2 As shown:

[0107] In the figure, Δz is the depth error of the imaging point caused by inaccurate parameter field. β is the inclination of the stratum near the point, θ is the incident angle of the seismic wave, and L1+L2 is the change in the ray path. From the figure, the conversion relationship between the depth residual and the travel time residual is:

[0108]

[0109] Where Δv represents the velocity perturbation between the initial model and the true model.

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

[0111] The tomographic velocity inversion is performed by meshing the geological model and using ray tracing to obtain the travel time information and ray path in the grid. Ray tracing is performed on the velocity field of the triangular mesh using the slowness square gradient variable. The expression structure is simple and accurate, and there is an analytical solution in a single grid. 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 the functional relationship between the nth power of slowness and the change of spatial position. According to the eigenvalue method, the ray tracing equation can be obtained:

[0114]

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

[0116]

[0117] In the formula, s 00 is the squared slowness of the reference point, s x ′ and s z ′ represents the gradient of the square slowness in the x and z directions respectively. Combining the above formula, the slowness vector expression can be obtained:

[0118]

[0119] Combining the above formulas for solution, the ray coordinates and ray travel time analytical solutions required for tomographic inversion can be obtained as follows:

[0120]

[0121] Using 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, tomographic inversion is performed to obtain the slowness update and complete this iteration.

[0123] Radon variation is the basic theory of travel time tomography. Its basic principle is to construct underground parameter fields by back-projecting seismic data along the ray direction to obtain underground medium velocity information. Among them, ray travel time tomography integrates the slowness in the travel time propagation direction to obtain the corresponding ray travel time:

[0124]

[0125] Where: s is the inverse of velocity, l represents the integral path, and the upper and lower limits of the integral are the detection point and the source, respectively. The change in model slowness will cause disturbances in travel time. The introduction of travel time residual vector Δt and slowness update Δs is:

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

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

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

[0129] (107) Perform Gaussian beam prestack depth migration on the updated velocity field to obtain angle gathers and imaging profiles. Determine whether to perform the next iteration based on the degree of flattening of the angle gathers, the accuracy of the picked layers, and the velocity accuracy requirements. If the iteration is to be continued, return to step 101 and repeat the process; otherwise, exit the loop to obtain the final velocity field.

[0130] The following are several specific embodiments of the present invention.

[0131] Example 1

[0132] The present invention uses a complex equivalent model as an embodiment to compare and analyze the designed triangular mesh and rectangular mesh ray tracing.

[0133] Using the equivalent complex parameter field, triangular mesh and rectangular mesh were divided respectively, and the excitation was carried out at 4 km above the surface, and 45 rays were evenly shot downward with a ray interval of 2. It can be seen from the ray tracing diagram that ( Figure 4 Figure 5 ):

[0134] In terms of boundary depiction accuracy: the arc interface of the triangular grid is completely consistent with the shape of the stratum interface. Compared with the parameter field rectangular grid division method, it has higher processing accuracy for stratum interfaces with greater undulations.

[0135] Regarding the accuracy of interface processing such as interface reflection and transmission: In the triangular meshing method, the mesh boundary is the interface of the parameter field, while the rectangular meshing parameter interface is usually inside the mesh. Therefore, in the ray tracing process, the calculation of whether there is an interface in the mesh can be omitted, which improves the operation efficiency. In addition, triangular mesh ray tracing omits the smoothing preprocessing required for the problem of velocity mutation in the mesh based on rectangular meshing, reducing the processing error. Triangular mesh ray tracing can specify the reflection surface of the ray, which can more accurately depict the details of the specific layer of the tomography process.

[0136] According to the complexity of the fluctuation of the velocity interface, the corresponding triangular grid control points are selected, and the number of grid divisions for different structural forms is adjusted. The present invention introduces formation dip constraints to the problem of selecting control points, densifies the number of grids in complex undulating structural areas, and reduces the number of grids in relatively flat and simple structural forms. Compared with the rectangular grid division method, the number of grid divisions can be reduced, so that the number of grids traversed during ray tracing is reduced, which effectively reduces the amount of calculation and storage, and correspondingly improves the forward modeling and inversion accuracy.

[0137] Therefore, the triangular grid ray tracing of the present invention is superior to the conventional rectangular grid ray tracing.

[0138] Example 2

[0139] The present invention uses a simulated data embodiment to apply the method to the data to perform triangular grid velocity modeling based on angle gathers to verify the effect of the method. The specific flow chart is shown in Figure 3 The steps of this embodiment are the same as the inversion strategy steps above:

[0140] The simulation data size is 801*450km( Figure 6 ), the real velocity field is smoothed horizontally and vertically as the initial velocity field ( Figure 7 ), using the initial velocity field for migration imaging. Figure 8 It can be seen that due to the parameter field error, the diffraction wave does not converge, the layer illusion appears, the resolution is reduced, and the angle gather has an obvious upward bend ( Fig. 9 ) and other issues.

[0141] Perform triangular grid travel time tomographic inversion on the initial velocity field to obtain the tomographic velocity field ( Fig.10 ) and perform migration imaging. Compare the initial parameter field migration profile and invert the parameter field migration profile ( Fig.11 ) The medium resolution is improved, the deep energy is enhanced, the diffraction wave is completely suppressed, and the angle gather is completely flattened ( Fig.12 ), and finally construct a high-precision inversion velocity field. Compared with rectangular grids, for deep complex structures, the inversion results of triangular grids are more consistent with the actual model situation, improving the problems of jagged and burred shapes caused by grid division, and the inversion results are closer to the actual formation morphology. It proves the adaptability of the strategy proposed in this invention to large offset data.

[0142] Example 3

[0143] In order to further verify the effectiveness and processing effect of the method of the present invention, rectangular grid travel time tomography and triangular grid velocity tomography based on node data proposed by the present invention are respectively performed on the actual work area data. The steps of this embodiment are the same as the inversion strategy steps above:

[0144] The root mean square velocity field is obtained by prestack velocity analysis method, and the initial velocity field is obtained by conversion according to Dix formula. The rectangular grid tomography velocity field is processed by conventional rectangular grid travel time tomography ( Fig.13 ) and perform prestack depth migration. In the obtained tomographic migration section ( Fig.14 ) and angle gathers ( Fig.15 ). The initial velocity field is triangulated and brought into the triangulated mesh travel time tomographic inversion process proposed in the present invention to obtain the triangulated mesh tomographic inversion velocity field ( Fig.16 ), perform prestack depth migration, and obtain the triangular grid level migration profile ( Fig.17 ) and angle gathers ( Fig.18). By comparing the offset profiles and angle gathers of the two tomography methods, it can be concluded that the velocity field constructed by the triangular grid velocity tomography method of the present invention has a higher imaging resolution at large offset distances and better layer continuity; the angle gathers are completely flattened, and for deep structures, the imaging effect of the tomography velocity field of the present invention is superior to the conventional travel time tomography method in terms of resolution, event continuity and clear depiction of structural morphology.

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

[0146] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0147] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A triangular grid tomographic velocity inversion optimization method based on node data, characterized in that: The triangular grid tomography velocity inversion optimization method based on node data includes: Step 1, obtaining imaging sections and angle domain common imaging point gathers; Step 2, picking the layer on the imaging section; Step 3, triangulate the velocity field; Step 4, perform residual curvature analysis on the diagonal gathers to obtain depth residuals and convert them into travel time information; Step 5, ray tracing is performed on the velocity field of the triangular mesh to obtain the ray path; Step 6: Using the triangular mesh ray path and travel time information, perform tomographic inversion to obtain the slowness update amount and complete this iteration; Step 7: Perform Gaussian beam prestack depth migration on the updated velocity field to obtain angle gathers and imaging profiles.

2. The triangular grid tomography velocity inversion optimization method based on node data according to claim 1 is characterized in that: In step 1, based on the input rectangular grid initial velocity field, the node data is prestacked and depth migrated to obtain the imaging section and the angle domain common imaging point gather.

3. The triangular grid tomography velocity inversion optimization method based on node data according to claim 2 is characterized in that: In step 1, the Gaussian beam prestack depth migration method is selected to obtain Gaussian beam prestack depth migration imaging sections and angle gathers of node data.

4. The triangular grid tomography velocity inversion optimization method based on node data according to claim 3 is characterized in that: In step 1, the kinematic ray tracing method is first used to obtain the Gaussian beam at the shot point and the detection point, and then the wave field of the shot point and the center of the Gaussian beam is forwardly extended and reversely extended respectively, and imaging is performed according to the cross-correlation imaging formula. Then, in the process of wave field extension, the travel time information at the center of the Gaussian beam and the node information around the center are used to calculate the angle information of the incident wave field and the reflected wave field at different imaging points to obtain the angle gather.

5. The triangular grid tomography velocity inversion optimization method based on node data according to claim 4 is characterized in that: In step 1, the formula for imaging by the cross-correlation imaging formula is: Where ω is the angular frequency, x, y, and z represent spatial positions, subscript s represents the source point, subscript r represents the receiving point, x without a subscript represents the imaging position, u is the extended wave field, and G * represents the Green's function, I pre is the imaging value.

6. The triangular grid tomography velocity inversion optimization method based on node data according to claim 5, characterized in that: 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 at point B that intersects the ray perpendicularly through the point: 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, that is, the coordinate value in the ray center coordinate system with point B as the origin, and T represents the travel time; Taking the derivative of both sides of equation (2) with respect to x, we get: Let l x ,l z Respectively represent the component values ​​of the tangent unit vector at point B in the x and z directions in the rectangular coordinate system; then: n=(x - x B )l z -(z - z B )l x (Equation 4) 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, that is, the angle between the ray and the positive direction of the Z axis: in After obtaining the propagation angles of the Gaussian beam at the source point and the receiving point using the above formula, the opening angle during the migration process can be obtained, and then the imaging values ​​are arranged according to the angle to complete the extraction of the angle gather.

7. The triangular grid tomography velocity inversion optimization method based on node data according to claim 1, characterized in that: In step 4, the travel time tomography method based on angle gathers is adopted to obtain angle gathers according to the Gaussian beam prestack depth migration method, and the travel time information is provided for the tomographic inversion equation by establishing a conversion relationship between the depth residual and the travel time residual.

8. The triangular grid tomography velocity inversion optimization method based on node data according to claim 7, characterized in that: In step 4, the conversion relationship between depth residual and travel time residual is: Among them, Δz is the depth error of the imaging point caused by inaccurate parameter field; β is the inclination of the stratum near the imaging point, θ is the incident angle of the seismic wave, L1+L2 is the change in the ray path, and Δv represents the velocity disturbance between the initial model and the true model.

9. The triangular grid tomography velocity inversion optimization method based on node data according to claim 1, characterized in that: In step 5, the tomographic velocity inversion is performed by meshing the geological model and using ray tracing to obtain the travel time information and ray paths within the grid. The slowness square gradient variable is used to perform ray tracing on the velocity field of the triangular mesh.

10. The triangular grid tomography velocity inversion optimization method based on node data according to claim 9, characterized in that: In step 5, the Hamilton-Jacobi operator of the eikonal equation in the slowness square medium is: In the above formula, p i is the slowness vector, is the functional relationship between the nth power of slowness and the change of spatial position; according to the eigenvalue method, the ray tracing equation can be obtained: In the above formula, p1 and p2 represent the components of the ray parameters in the x and z directions in the spatial position respectively; when the square of the slowness in the medium changes linearly with the spatial position, the square of the slowness expression can be rewritten as: In the formula, s 00 is the squared slowness of the reference point, s′ x and s′ z Respectively represent the gradient of the square slowness in the x and z directions. Combining the above formula, the slowness vector expression can be obtained: Combining the above formulas for solution, the ray coordinates and ray travel time analytical solutions required for tomographic inversion can be obtained as follows: Using 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.

11. The triangular grid tomography velocity inversion optimization method based on node data according to claim 1, characterized in that: In step 6, Radon variation is the basic theory of travel time tomography. Its basic principle is to obtain underground medium velocity information by back-projecting seismic data along the ray direction to construct underground parameter fields. Among them, ray travel time tomography integrates the slowness in the travel time propagation direction to obtain the corresponding ray travel time: Where: s is the inverse of velocity, l represents the integral path, the upper and lower limits of the integral are the detection point and the source respectively; the change of model slowness will cause disturbance in travel time, and the introduction of travel time residual vector Δt and slowness update Δs is: LΔs=Δt (Equation 13) According to the above formula, the relationship between the change in ray travel time and the slowness disturbance can be established; where L represents the length of the ray path; The travel time residual and ray path are introduced into the tomographic equation, the velocity update is obtained by solving it, and the inversion velocity field is constructed by repeated iterations.

12. The triangular grid tomography velocity inversion optimization method based on node data according to claim 1, characterized in that: In step 7, Gaussian beam prestack depth migration is performed on the updated velocity field to obtain angle gathers and imaging profiles. Based on the degree of flattening of the angle gathers, the accuracy of the picked layers, and the velocity accuracy requirements, it is determined whether to perform the next iteration. If the iteration is continued, return to step 1 and repeat the process; otherwise, exit this loop to obtain the final velocity field.

Images