CN115201913B - Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method - Google Patents
Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method Download PDFInfo
- Publication number
- CN115201913B CN115201913B CN202210889256.7A CN202210889256A CN115201913B CN 115201913 B CN115201913 B CN 115201913B CN 202210889256 A CN202210889256 A CN 202210889256A CN 115201913 B CN115201913 B CN 115201913B
- Authority
- CN
- China
- Prior art keywords
- imaging
- node
- grid
- finite difference
- free
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
- G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Geophysics (AREA)
- Pure & Applied Mathematics (AREA)
- Acoustics & Sound (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Mathematical Optimization (AREA)
- Geology (AREA)
- Environmental & Geological Engineering (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Computer Graphics (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to the technical field of geophysical exploration, in particular to a least square reverse time migration imaging method, a system and a storage medium based on a gridless finite difference method. According to the method provided by the invention, through the grid-free distribution generated according to the underground model of the target work area, the initial imaging result is obtained by using the grid-free finite difference method for reverse time migration, and the least square reverse time migration imaging method is performed by using the migration and reverse migration operators based on the grid-free finite difference method, so that the accurate description and the high-efficiency imaging of the irregular structures of the complex earth surface and underground are realized, and the underground complex geological structure can be finely depicted under any complex earth surface condition. Meanwhile, the method effectively reduces the algorithm calculation and memory requirements, improves the imaging precision, and can provide more accurate imaging results for seismic exploration.
Description
Technical Field
The invention relates to the technical field of geophysical exploration, in particular to a least square reverse time migration imaging method, a system and a storage medium based on a gridless finite difference method.
Background
With the increasing demand for seismic exploration, the current industry is increasingly demanding in terms of imaging technology. Among a plurality of offset imaging technologies, the inverse time offset technology based on the double-pass wave equation has extremely high imaging precision, can describe the space position and the geometric structure of the underground reflection interface of any ray path, and can provide reliable data for construction fine interpretation and implementation of well positions. However, the conventional reverse time shift method has the inherent disadvantages of strong low-frequency interference, uneven illumination, poor amplitude preservation and the like. The least square reverse time migration method is an inversion method of wave equation based on reverse time migration and taking reflection coefficient as parameter, can improve various defects and shortages in reverse time migration, and effectively improves the precision of seismic wave migration imaging.
However, the conventional least square inverse time migration imaging method mostly adopts a finite difference method of a regular grid to perform seismic wave numerical simulation, grid units of the method are usually rectangular, positions of a bending interface or an irregular model boundary are difficult to accurately describe, and local refinement is difficult to achieve in a key area. This results in scattered noise at irregular boundaries when numerically computing the subsurface wavefield, severely affecting the accuracy of the wavefield computation. In order to solve the problems of abnormal bodies with any shape or irregular terrains, a plurality of irregular grids finite difference methods are proposed, and although the methods have certain geometric flexibility, the problems of model interface or boundary adaptability can be solved, complex grid subdivision or grid mapping processes are often required, so that the numerical simulation cost is increased. In addition, the finite element method and the spectral element method can solve the problems to a certain extent, but the demand of the finite element method and the spectral element method on the self is far higher than that of the finite difference method on the calculation amount and the calculation memory, and the applicability of the least square reverse time offset method with larger calculation amount is poor. The grid-free finite difference method is used as an emerging numerical value calculation method, can effectively solve the problems of the calculation method, and is preliminarily applied to the field of seismic exploration. For the least square inverse time shift imaging method, the gridless finite difference method can be effectively applied.
At present, the domestic exploration focus is changed from the eastern plain area to the western plateau area and the mountain area, and most of the western exploration areas are faced with complex surface and underground structures. In order to cope with increasing exploration demands and solve the problems that the conventional least square reverse time migration method is difficult to accurately describe irregular surface and underground interfaces and excessive in calculation amount and memory demands, a high-efficiency imaging method capable of accurately describing complex surface and underground irregular structures is needed, and high-quality imaging results are provided for geological screening and energy exploration.
Disclosure of Invention
In view of the above, a first object of the present invention is to provide a least squares inverse time migration imaging method based on a gridless finite difference method, which achieves accurate description and efficient imaging of complex surface and subsurface irregular structures by generating gridless distributions according to a target work area subsurface model, and using migration and inverse migration operators based on the gridless finite difference method.
Based on the same inventive concept, a second object of the present invention is to provide a least squares inverse time offset imaging system based on a gridless finite difference method.
Based on the same inventive concept, a third object of the present invention is to provide a storage medium.
The first object of the present invention can be achieved by adopting the following technical scheme:
a least squares reverse time shift imaging method based on a gridless finite difference method, comprising the steps of:
observing the seismic record, and acquiring a target work area underground speed model according to the observed seismic record;
generating mesh-free node distribution suitable for an underground speed model according to the underground speed model of a target work area, and discretizing the underground speed model;
applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic records and the discretized underground speed model to obtain an initial imaging result;
according to the initial imaging result or the updated imaging result, a positive algorithm based on a grid-free finite difference method is applied to calculate to obtain a simulated seismic record;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the recorded residual errors, an offset operator applied to a grid-free finite difference method is used for obtaining updated gradients of imaging results corresponding to the recorded residual errors, and updated imaging results are obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating.
Further, in the underground speed model of the target work area, the degree of node discretization is positively related to the size of the speed field.
Further, generating a network-free node distribution adapted to the subsurface velocity model from the target work area subsurface velocity model, comprising the steps of:
determining the size of a calculation domain according to the underground speed model of the target work area, and determining the mapping relation between the radius of the grid-free distributed particles and the size of a speed field;
randomly generating node distribution at the bottom of the underground speed model of the target work area, giving the density of the node distribution according to a preset value, and defining the node positions as potential node positions;
selecting a potential node position with the lowest position as a new effective node position;
determining a circular area by taking a new effective node as a circle center and taking the radius of the particle at the effective node as the radius of the circle, and deleting all nodes except the effective node in the circular area;
marking a nearest potential node from the left side and the right side of the circle center, determining the connecting direction of the circle center and the two potential nodes, and placing a plurality of new potential nodes on the arc clamped in the two directions at equal intervals;
and continuing to select a potential node position with the lowest position as a new effective node position, and repeating the steps until the node fills the whole calculation domain.
Further, according to the observed seismic record and the discretized underground velocity model, inverse time migration imaging based on a gridless finite difference method is applied to obtain an initial imaging result, and the method comprises the following steps:
constructing a linear equation set by using a radial basis function, and solving the constructed linear equation set to obtain a gridless differential coefficient;
constructing a wave equation propagation operator according to the gridless differential coefficient;
and according to the wave equation propagation operator, applying reverse time migration imaging based on a gridless finite difference method, and acquiring the initial imaging result.
Further, a linear equation set constructed by using the radial basis function is specifically formed as follows:
wherein phi (||X-X) i I) represents a radial basis function, only related to coordinate positions, subscripts i=0, 1,2 … n represent the i-th node in the differential template, where n represents the number of nodes in the differential template, X represents the spatial coordinates of mesh-free nodes, i·i represents a two-norm,
representing the Laplace operator, c i Is the gridless differential coefficient to be solved.
Further, adding additional basis functions to the linear equation set constructed by the radial basis functions eliminates static errors and makes the convergence rate of the errors reach the highest order of the additional added polynomials.
Further, the additional basis function added in the linear equation set constructed by using the radial basis function is Taylor single equation, and the specific form of the obtained linear equation set is as follows:
wherein c n+1 、c n+2 And c n+3 For auxiliary coefficients, A is a radial basis function matrix, c is a differential coefficient matrix, Δx i Is X i Node and X 0 Difference in horizontal coordinates of nodes, Δz i Is X i Node and X 0 The difference in the vertical direction coordinates of the nodes,
represents X-X 0 Radial basis function at that time.
Further, according to the gridless differential coefficient, constructing a wave equation propagation operator, which comprises the following steps:
solving for Laplace operator
The corresponding differential coefficient is used for further solving a spatial derivative term in the wave equation, and the specific form is as follows:
where v is velocity, p is seismic wavefield, t is time, c i Is a difference coefficient, p i A wavefield value at an ith node in the differential template;
estimating a time derivative term in the wave equation by adopting a second-order center difference format, wherein the specific form is as follows:
wherein p is the seismic wavefield, p 0 For the background wavefield, t is time, superscripts-1, 0, and 1 represent the previous, current, and next time, respectively, and Δt is the time step.
The second object of the present invention can be achieved by adopting the following technical scheme:
a grid-less finite difference method based least squares reverse time offset imaging system comprising:
the speed model unit is used for observing the seismic data and generating a target work area underground speed model according to the observed seismic data;
the grid-free node generating unit is used for generating grid-free node distribution suitable for the underground speed model of the target work area according to the underground speed model of the target work area, and discretizing the underground speed model;
the initial imaging unit is used for obtaining an initial imaging result by applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic record and the discretized underground speed model;
a least squares offset imaging unit for imaging using a least squares offset imaging principle, comprising:
according to the initial imaging result or the updated imaging result, a positive algorithm based on a grid-free finite difference method is applied to calculate to obtain a simulated seismic record;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the recorded residual error, an offset operator applied to a grid-free finite difference method is used for obtaining an updated gradient of an imaging result, and then the updated imaging result is obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating.
The third object of the present invention can be achieved by adopting the following technical scheme:
a storage medium storing a program which, when executed by a processor, implements the above-described least squares inverse time-shift imaging method based on a gridless finite difference method.
Compared with the prior art, the invention has the following beneficial effects:
(1) The grid-free node generation algorithm in the least square inverse time migration imaging method based on the grid-free finite difference method can generate a curved interface and an irregular boundary which are suitable for any model, the node distribution density can be based on grid-free node distribution preset by the model speed, the generation process is simple and efficient, the algorithm can be guaranteed to have extremely high geometric flexibility without high calculation cost, and the memory requirement of the algorithm can be effectively reduced.
(2) The grid-free finite difference coefficient solving and wave equation operator constructing method ensures the accuracy of wave field calculation under grid-free node distribution, and meanwhile, the wave field simulation is not influenced by scattering noise approximately introduced by 'step-like' at irregular boundaries in the method, so that convenience is provided for implementation of a high-accuracy least square reverse time migration imaging algorithm, and convergence efficiency of the imaging algorithm is improved.
(3) In general, the method provided by the invention can finely depict underground complex geological structures under any complex earth surface condition, and compared with the traditional method, the method has the advantages that the imaging precision is improved while the algorithm calculation and the memory requirement are effectively reduced, and more accurate imaging results can be provided for seismic exploration.
Drawings
Fig. 1 is a flow chart of a least squares inverse time shift imaging method based on the gridless finite difference method according to embodiment 1 of the present invention.
Fig. 2 is a comparative schematic diagram of a regular grid distribution and a no grid distribution provided in example 1 of the present invention.
Fig. 3 is a schematic diagram of a grid-free node distribution rapid generation algorithm according to embodiment 1 of the present invention, wherein (a) is a schematic diagram of step S220, (b) and (c) are schematic diagrams of step S240, and (d) is a schematic diagram of step S250.
FIG. 4 is a schematic representation of a model of subsurface velocities for a target work area having three sets of formations according to example 1 of the present invention.
FIG. 5 is a schematic diagram of the number of nodes required for different node discrete modes in offset imaging of the three-layer velocity model in example 1 of the present invention.
Fig. 6 is a schematic diagram of imaging results obtained based on a mesh-free finite difference least squares reverse time offset for the three-layer velocity model in embodiment 1 of the present invention.
Fig. 7 is a schematic diagram of an imaging result obtained by performing conventional rule-grid-based finite difference least squares reverse time offset acquisition on the three-layer velocity model in embodiment 1 of the present invention.
FIG. 8 is a schematic representation of a model of subsurface velocities for a target work area having a complex surface in accordance with example 1 of the present invention.
FIG. 9 is a schematic diagram of imaging results obtained based on a conventional reverse time migration without grid finite difference for the complex surface velocity model in embodiment 1 of the present invention.
FIG. 10 is a schematic diagram of imaging results obtained based on grid-free finite difference least squares reverse time migration for the complex surface velocity model in example 1 of the present invention.
Fig. 11 is a schematic diagram of a least squares inverse time offset imaging system based on the gridless finite difference method in embodiment 2 of the present invention.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are some embodiments of the present invention, but not all embodiments, and all other embodiments obtained by persons of ordinary skill in the art without making any inventive effort based on the embodiments of the present invention are within the scope of protection of the present invention.
Example 1:
as shown in fig. 1, the present embodiment provides a least squares reverse time shift imaging method based on a gridless finite difference method, which includes the steps of:
s100, observing the seismic record, and acquiring a target work area underground speed model according to the observed seismic record, wherein the method comprises the following steps of:
s110, observing the seismic record. In this embodiment, the observed seismic record is active source seismic data received by a high-density geophone array, which is excited by an artificial explosive source (explosive in the land work area and air gun in the marine work area). In order to improve the seismic imaging quality and the record transformation effect, operations such as denoising (removing random noise, surface wave noise, radio interference, adjacent shot interference, linear noise and the like), cutting (cutting off refractive wave components in seismic records), deconvolution and the like can be performed on the observed seismic data.
S120, acquiring a target work area underground speed model according to the observed seismic records. In this embodiment, according to the seismic record obtained by observation, a relatively accurate underground velocity model of the target work area is obtained by means of a tomography or full waveform inversion constant velocity modeling means.
S200, as shown in fig. 3, generating mesh-free node distribution suitable for the underground speed model according to the underground speed model of the target work area, and discretizing the underground speed model.
The step of generating the gridless node distribution is to discrete the subsurface velocity model, and after generating the gridless node distribution, the subsurface wavefield will be solved according to discrete values of the nodes. The method comprises the following steps:
s210, determining the size of a calculation domain according to an underground speed model of a target work area, and determining the mapping relation between the radius of grid-free distributed particles and the size of a speed field; in this embodiment, the radius of the particles distributed without grid is positively correlated with the size of the velocity field, i.e. the larger the velocity, the larger the corresponding particle radius;
s220, as shown in a sub-graph (a) of FIG. 3, randomly generating node distribution at the bottom of the underground speed model of the target work area, wherein the density of the node distribution is given according to a preset value, and defining the node positions as potential node positions (PDP);
s230, selecting a potential node position with the lowest position as a new effective node position;
s240, as shown in sub-graphs (b) and (c) of FIG. 3, determining a circular area by taking a new effective node as a circle center and taking the radius of the particle at the effective node as the radius of the circle, and deleting all nodes except the effective node in the circular area;
s250, as shown in a sub-graph (d) of FIG. 3, marking a nearest potential node from the left side and the right side of the circle center, determining the connecting direction of the circle center and the two potential nodes, and placing 5 new potential nodes on an arc clamped in the two directions at equal intervals;
and S260, continuing to select a potential node position with the lowest position as a new effective node position, and repeating the steps S240-S260 until the node fills the whole calculation domain.
S300, applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic record and the discretized underground speed model to obtain an initial imaging result, wherein the method comprises the following steps of:
s310, constructing a linear equation set by utilizing a radial basis function, and solving the constructed linear equation set to obtain the gridless differential coefficient.
In this embodiment, a system of linear equations constructed using radial basis functions is shown in the following formula:
the above linear system of equations may be abbreviated as ac=lΦ. Wherein phi (||X-X) i I) represents a radial basis function, which is a function related to the coordinate position only, the subscript i=0, 1,2 … n represents the i-th node in the differential template, where n represents the number of nodes in the differential template, X represents the spatial coordinates of the mesh-free node generated in step S200, i·i represents a two-norm,
representing the Laplace operator, c i Is the gridless differential coefficient to be solved.
Additional basis functions (typically polynomials) may be added to the above-described system of linear equations constructed using radial basis functions to eliminate static errors and to bring the convergence rate of the errors to the highest order of the additional added polynomials. In this embodiment, a Taylor single equation is added to the above-mentioned linear equation set to obtain a linear equation set shown in the following equation:
wherein c n+1 、c n+2 And c n+3 For auxiliary coefficients, A is a radial basis function matrix, c is a differential coefficient matrix, Δx i Is X i Node and X 0 Difference in horizontal coordinates of nodes, Δz i Is X i Node and X 0 The difference in the vertical direction coordinates of the nodes,
represents X-X 0 Radial basis function at that time.
S320, constructing a wave equation propagation operator according to the gridless differential coefficient, wherein the wave equation propagation operator comprises the following steps of:
s321, solving Laplacian operator
The corresponding differential coefficient is used for further solving a spatial derivative term in the wave equation, and the specific form is as follows:
where v is velocity, p is seismic wavefield, t is time, c i Is a difference coefficient, p i Is the wavefield value at the ith node.
S322, estimating a time derivative term in the wave equation by adopting a second-order center differential format, wherein the specific form is as follows:
wherein p is the seismic wavefield, t is the time, the superscripts-1, 0 and 1 represent the previous, current and next times, p 0 For the background wavefield, Δt is the time step.
S330, applying reverse time migration imaging based on a gridless finite difference method according to a wave equation propagation operator to obtain an initial imaging result, wherein the method comprises the following steps of:
s331, inputting the observed seismic records through the wave equation propagation operator, and calculating a seismic source wave field and a detection wave field, wherein the specific form is shown as follows:
wherein p is s And p r Respectively a source wave field and a detection wave field, x s And x r The spatial positions of the source and the detector respectively, S is the function of the source, d obs To observe seismic records。
S332, applying imaging conditions to the seismic source and the detected wave field, wherein in the embodiment, the applied imaging conditions are cross-correlation imaging conditions, and the specific form is as follows:
I(x)＝∫∫p s (x，x s ，t)p r (x，x r ，t)dtdx s
wherein I represents the imaging result.
S400, according to an initial imaging result or an updated imaging result, a positive algorithm based on a grid-free finite difference method is applied, and a simulated seismic record is obtained through calculation, wherein the specific form is as follows:
wherein p is 0 And δp is the background wavefield and the disturbance wavefield, respectively, and d is the simulated seismic record.
S500, calculating a record residual between the simulated seismic record and the observed seismic record, wherein in the embodiment, the simulated seismic record obtained in the step S400 is differenced with the seismic record observed in the step S110 to obtain the record residual.
S600, according to the recorded residual, the method is applied to an offset operator without grid finite difference method, an updated gradient of an imaging result corresponding to the recorded residual is obtained, and then the updated imaging result is obtained through a linear inversion solver, and the method comprises the following steps:
s610, inputting a record residual, calculating an imaging result update gradient by using an offset operator based on a gridless finite difference method, wherein the specific form is as follows:
wherein g is the imaging result updated gradient, lambda is the background wavefield p 0 V of the concomitant wave field of (2) 0 Is the background velocity field.
S620, acquiring an updated imaging result through an optimization method, such as a steepest descent method, a conjugate gradient method and the like.
S700, judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating.
In this embodiment, the iteration method is as follows: the new imaging result acquired in step S620 is input as the imaging result of step S400, and steps S400 to S700 are performed again.
In this embodiment, the preset termination conditions are: any one of the preset conditions is satisfied, including:
the method comprises the steps that under a first preset condition, an objective function related to residual errors is smaller than a set threshold value;
and the algorithm iteration times are larger than a set threshold value under a second preset condition.
Fig. 2 is a schematic diagram comparing regular grid distribution and non-grid distribution, and it can be seen that for an irregular interface (curved solid line) in the calculation domain, when the regular grid distribution is discrete, a "step-like" approximation (black dotted line) appears, and when the non-grid distribution is discrete, the node distribution can be well attached to the interface, and scattering noise introduced when the "step-like" approximation is used for numerically calculating the subsurface wave field can be effectively avoided. FIG. 3 is a schematic diagram of a grid-free node distribution rapid generation algorithm according to which a subsurface velocity model may be rapidly discretized by grid-free nodes.
Fig. 4 is a model of the underground speed of the target work area with three strata, which is adopted in the embodiment, and has a size of 3km×2km, fig. 5 shows the number of nodes required for respectively performing regular grid dispersion and no grid dispersion on the three-layer speed model, and it can be seen that the no grid can be used for dispersing the speed model only by using the number of nodes of 2/3 of the regular grid dispersion, which means that when the imaging algorithm calculation is performed, the calculation and the memory requirement required by the algorithm are correspondingly reduced.
Fig. 6 and 7 are imaging results obtained by using a least squares inverse time offset imaging method based on a mesh-free finite difference and a conventional least squares inverse time offset imaging method based on a regular mesh finite difference, respectively, in the present embodiment. It can be seen that both describe subsurface structures and provide amplitude information similar to a true reflectivity model. However, due to the stepped approximation of the regular grid, the in-phase axis continuity in the imaging results of fig. 7 is poor, which may be more severe and even affect algorithm convergence when imaging more complex subsurface formations. While the imaging results of fig. 5 can almost perfectly describe the morphology of the subsurface structure. By comparing the imaging results of fig. 6 and fig. 7, it can be found that the least square inverse time offset imaging method based on the mesh-free finite difference in the present embodiment can provide a more accurate underground imaging result on the premise of effectively reducing the algorithm calculation and the memory requirement.
FIG. 8 is a model of the subsurface speed of a target work area with a complex surface, with a size of 4.16km 2.49km, as used in this embodiment. Fig. 9 and 10 are imaging results obtained using a conventional inverse time-shift imaging method and a least squares inverse time-shift imaging method based on a gridless finite difference, respectively, in an embodiment of the present invention. It can be seen that the imaging accuracy at the surface is high in the imaging results, and the imaging results of fig. 10 are significantly higher resolution than the imaging results of fig. 9, which can provide a more reliable reflectivity model for construction and lithology interpretation. The adaptability of the least square inverse time migration imaging method based on the grid-free finite difference to complex earth surface conditions in the embodiment of the invention is also verified, and the capability of providing high-precision imaging results for seismic exploration under complex geological conditions by the method is demonstrated.
Therefore, the grid-free node generation algorithm in the least square inverse time migration imaging method based on the grid-free finite difference method can generate grid-free node distribution which is adaptive to a bending interface and an irregular boundary of any model and can be preset based on the model speed, the generation process is simple and efficient, the algorithm can be guaranteed to have extremely high geometric flexibility without high calculation cost, and the memory requirement of the algorithm can be effectively reduced; meanwhile, the gridless finite difference coefficient solving and wave equation operator constructing method adopted by the embodiment ensures the accuracy of wave field calculation under gridless node distribution, and meanwhile, the wave field simulation is not influenced by scattering noise approximately introduced by the step-shaped at the irregular boundary in the method, so that a convenient condition is provided for implementing a high-accuracy least square reverse time migration imaging algorithm, and the convergence efficiency of the imaging algorithm is improved; in general, the method provided by the embodiment can finely depict underground complex geological structures under any complex earth surface condition, and compared with the traditional method, the method has the advantages that the imaging precision is improved while the algorithm calculation and the memory requirements are effectively reduced, and more accurate imaging results can be provided for seismic exploration.
Example 2:
as shown in fig. 11, the present embodiment provides a least squares inverse time shift imaging system based on a gridless finite difference method, including:
a velocity model unit 1101, configured to observe seismic data, and generate a target work area underground velocity model according to the observed seismic data;
the mesh-free node generating unit 1102 generates mesh-free node distribution adapted to an underground speed model according to the underground speed model of the target work area, and discretizes the underground speed model;
an initial imaging unit 1103, which applies reverse time migration imaging based on a grid-free finite difference method according to the observed seismic record and the discretized underground velocity model to obtain an initial imaging result;
a least squares offset imaging unit 1104 for imaging using the least squares offset imaging principle, comprising:
according to the imaging result, a positive algorithm based on a grid-free finite difference method is applied, and a simulated seismic record is obtained through calculation;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the residual error, an offset operator without grid finite difference method is applied to obtain an updated gradient of an imaging result, and then the updated imaging result is obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting an imaging result; otherwise, continuing to calculate the residual error and the imaging result, and iterating.
Example 3:
the present embodiment provides a storage medium, which is a computer readable storage medium storing a computer program, where the computer program is executed by a processor to implement the least squares inverse time offset imaging method based on the gridless finite difference method of the above embodiment 1, specifically as follows:
observing the seismic record, and acquiring a target work area underground speed model according to the observed seismic record;
generating mesh-free node distribution suitable for an underground speed model according to the underground speed model of a target work area, and discretizing the underground speed model;
applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic records and the discretized underground speed model to obtain an initial imaging result;
according to the initial imaging result or the updated imaging result, a positive algorithm based on a grid-free finite difference method is applied to calculate to obtain a simulated seismic record;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the recorded residual errors, an offset operator applied to a grid-free finite difference method is used for obtaining updated gradients of imaging results corresponding to the recorded residual errors, and updated imaging results are obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating.
The computer readable storage medium of the present embodiment may be a computer readable signal medium or a computer readable storage medium, or any combination of the two. The computer readable storage medium can be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or a combination of any of the foregoing. More specific examples of the computer-readable storage medium may include, but are not limited to: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a Random Access Memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.
In this embodiment, a computer-readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device. In the present embodiment, however, the computer-readable signal medium may include a data signal propagated in baseband or as part of a carrier wave, with a computer-readable program embodied therein. Such a propagated data signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination of the foregoing. A computer readable signal medium may also be any computer readable storage medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device. A computer program embodied on a computer readable storage medium may be transmitted using any appropriate medium, including but not limited to: electrical wires, fiber optic cables, RF (radio frequency), and the like, or any suitable combination of the foregoing.
The computer readable storage medium may be written in one or more programming languages, including an object oriented programming language such as Java, python, C ++ and conventional procedural programming languages, such as the C-language or similar programming languages, or combinations thereof for performing the present embodiments. The program may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the case of a remote computer, the remote computer may be connected to the user's computer through any kind of network, including a Local Area Network (LAN) or a Wide Area Network (WAN), or may be connected to an external computer (for example, through the Internet using an Internet service provider).
It is apparent that the above-described embodiments are only some embodiments of the present invention, but not all embodiments, and the present invention is not limited to the details of the above-described embodiments, and any appropriate changes or modifications made by those skilled in the art will be deemed to be within the scope of the present invention.
Claims (8)
1. A least squares reverse time shift imaging method based on a gridless finite difference method, comprising the steps of:
observing the seismic record, and acquiring a target work area underground speed model according to the observed seismic record;
generating mesh-free node distribution suitable for an underground speed model according to the underground speed model of a target work area, and discretizing the underground speed model;
applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic records and the discretized underground speed model to obtain an initial imaging result;
according to the initial imaging result or the updated imaging result, a positive algorithm based on a grid-free finite difference method is applied to calculate to obtain a simulated seismic record;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the recorded residual errors, the method is applied to an offset operator based on a grid-free finite difference method, an updated gradient of an imaging result corresponding to the recorded residual errors is obtained, and then the updated imaging result is obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating;
generating a gridless node distribution adapted to the subsurface velocity model from the target work area subsurface velocity model, comprising the steps of:
determining the size of a calculation domain according to the underground speed model of the target work area, and determining the mapping relation between the radius of the grid-free distributed particles and the size of a speed field;
randomly generating node distribution at the bottom of the underground speed model of the target work area, giving the density of the node distribution according to a preset value, and defining the node positions as potential node positions;
selecting a potential node position with the lowest position as a new effective node position;
determining a circular area by taking a new effective node as a circle center and taking the radius of the particle at the effective node as the radius of the circle, and deleting all nodes except the effective node in the circular area;
marking a nearest potential node from the left side and the right side of the circle center, determining the connecting direction of the circle center and the two potential nodes, and placing a plurality of new potential nodes on the arc clamped in the two directions at equal intervals;
continuously selecting a potential node position with the lowest position as a new effective node position, and repeating the steps until the node fills the whole calculation domain;
according to the observed seismic records and the discretized underground speed model, reverse time migration imaging based on a grid-free finite difference method is applied to obtain an initial imaging result, and the method comprises the following steps of:
constructing a linear equation set by using a radial basis function, and solving the constructed linear equation set to obtain a gridless differential coefficient;
constructing a wave equation propagation operator according to the gridless differential coefficient;
and according to the wave equation propagation operator, applying reverse time migration imaging based on a gridless finite difference method, and acquiring the initial imaging result.
2. The grid-free finite difference method-based least squares reverse time migration imaging method of claim 1, wherein the degree of node discretization in the subsurface velocity model of the target work area is positively correlated with the size of the velocity field.
3. The grid-free finite difference method-based least squares inverse time offset imaging method of claim 1, wherein the system of linear equations constructed using radial basis functions is in the specific form of:
wherein phi (|X-X) i Ii) denotes a radial basis function, which is related to the coordinate position only, the subscript i=0, 1,2 … n denotes the i-th node in the differential template, where n denotes the number of nodes in the differential template, X denotes the spatial coordinates of the mesh-free node, i·i denotes the two norms,
representing the Laplace operator, c i Is the gridless differential coefficient to be solved.
4. A method of grid-less finite difference based least squares reverse time migration imaging according to claim 3, wherein additional basis functions are added to the system of linear equations constructed using radial basis functions to eliminate static errors and to bring the convergence rate of the errors to the highest order of the additional additive polynomial.
5. The method of grid-free finite difference method-based least squares inverse time shift imaging of claim 4, wherein the additional basis function added to the system of linear equations constructed using radial basis functions is Taylor's single term, and the specific form of the system of linear equations obtained is:
wherein c n+1 、c n+2 And c n+3 For auxiliary coefficients, A is a radial basis function matrix, c is a differential coefficient matrix, Δx i Is X i Node and X 0 Difference in horizontal coordinates of nodes, Δz i Is X i Node and X 0 Nodes ofThe difference in the coordinates in the vertical direction,
represents X-X 0 Radial basis function at that time.
6. The method for least squares inverse time offset imaging based on the gridless finite difference method according to claim 1, wherein the wave equation propagation operator is constructed according to gridless differential coefficients, comprising the steps of:
solving for Laplace operator
The corresponding differential coefficient is used for further solving a spatial derivative term in the wave equation, and the specific form is as follows:
where v is velocity, p is seismic wavefield, t is time, c i Is a difference coefficient, p i A wavefield value at an ith node in the differential template;
estimating a time derivative term in the wave equation by adopting a second-order center difference format, wherein the specific form is as follows:
wherein p is the seismic wavefield, p 0 For the background wavefield, t is time, and superscripts-1, 0, and 1 represent the previous, current, and next times, respectively.
7. A grid-less finite difference method-based least squares inverse time offset imaging system, comprising:
the speed model unit is used for observing the seismic data and generating a target work area underground speed model according to the observed seismic data;
the grid-free node generating unit is used for generating grid-free node distribution suitable for the underground speed model of the target work area according to the underground speed model of the target work area, and discretizing the underground speed model;
the initial imaging unit is used for obtaining an initial imaging result by applying reverse time migration imaging based on a grid-free finite difference method according to the observed seismic record and the discretized underground speed model;
a least squares offset imaging unit for imaging using a least squares offset imaging principle, comprising:
according to the initial imaging result or the updated imaging result, a positive algorithm based on a grid-free finite difference method is applied to calculate to obtain a simulated seismic record;
calculating a record residual error between the simulated seismic record and the observed seismic record;
according to the recorded residual errors, the method is applied to an offset operator based on a grid-free finite difference method, an updated gradient of an imaging result corresponding to the recorded residual errors is obtained, and then the updated imaging result is obtained through a linear inversion solver;
judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuing to calculate the record residual error and update the imaging result, and iterating;
generating a gridless node distribution adapted to the subsurface velocity model from the target work area subsurface velocity model, comprising the steps of:
determining the size of a calculation domain according to the underground speed model of the target work area, and determining the mapping relation between the radius of the grid-free distributed particles and the size of a speed field;
randomly generating node distribution at the bottom of the underground speed model of the target work area, giving the density of the node distribution according to a preset value, and defining the node positions as potential node positions;
selecting a potential node position with the lowest position as a new effective node position;
determining a circular area by taking a new effective node as a circle center and taking the radius of the particle at the effective node as the radius of the circle, and deleting all nodes except the effective node in the circular area;
marking a nearest potential node from the left side and the right side of the circle center, determining the connecting direction of the circle center and the two potential nodes, and placing a plurality of new potential nodes on the arc clamped in the two directions at equal intervals;
continuously selecting a potential node position with the lowest position as a new effective node position, and repeating the steps until the node fills the whole calculation domain;
according to the observed seismic records and the discretized underground speed model, reverse time migration imaging based on a grid-free finite difference method is applied to obtain an initial imaging result, and the method comprises the following steps of:
constructing a linear equation set by using a radial basis function, and solving the constructed linear equation set to obtain a gridless differential coefficient;
constructing a wave equation propagation operator according to the gridless differential coefficient;
and according to the wave equation propagation operator, applying reverse time migration imaging based on a gridless finite difference method, and acquiring the initial imaging result.
8. A storage medium storing a program which, when executed by a processor, implements the grid-less finite difference method-based least squares inverse time offset imaging method of any one of claims 1 to 6.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210889256.7A CN115201913B (en) | 2022-07-27 | 2022-07-27 | Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210889256.7A CN115201913B (en) | 2022-07-27 | 2022-07-27 | Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115201913A CN115201913A (en) | 2022-10-18 |
| CN115201913B true CN115201913B (en) | 2023-05-12 |
Family
ID=83584493
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210889256.7A Active CN115201913B (en) | 2022-07-27 | 2022-07-27 | Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115201913B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116774292B (en) * | 2023-08-22 | 2023-11-24 | 浙江大学海南研究院 | Seismic wave travel time determining method, system, electronic equipment and storage medium |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101915938A (en) * | 2010-07-05 | 2010-12-15 | 中国科学院地质与地球物理研究所 | Offset imaging method and device for converted waves |
| CN104216012A (en) * | 2014-08-27 | 2014-12-17 | 中国科学院地质与地球物理研究所 | Three-dimensional Born-Kirchhoff variable-step interpolation imaging method |
| CN105717547A (en) * | 2015-12-22 | 2016-06-29 | 吉林大学 | Anisotropy medium magnetotelluric meshless value simulation method |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
| CN109946741A (en) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | Pure qP wave least square reverse-time migration imaging method in a kind of TTI medium |
| CN110133713A (en) * | 2019-04-24 | 2019-08-16 | 中国石油大学(华东) | A kind of the multiple wave least square reverse-time migration imaging method and system of full propagation path attenuation compensation |
| CN110824558A (en) * | 2019-11-20 | 2020-02-21 | 中国石油大学(华东) | Seismic wave numerical simulation method |
| CN110888166A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Least square offset imaging method and device based on L-BFGS algorithm |
| CN112384937A (en) * | 2018-05-12 | 2021-02-19 | 地质探索系统公司 | Seismic data interpretation system |
| CN112649859A (en) * | 2019-10-12 | 2021-04-13 | 中国石油化工股份有限公司 | Seismic wave speed self-adaptive grid-free field node establishment method and system |
| CN114357831A (en) * | 2021-12-29 | 2022-04-15 | 中海石油(中国)有限公司 | Non-grid generalized finite difference forward modeling method, device, storage medium and equipment |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10088596B2 (en) * | 2013-03-15 | 2018-10-02 | Schlumberger Technology Corporation | Meshless representation of a geologic environment |
| BR112015028562A2 (en) * | 2013-05-15 | 2017-07-25 | Logined Bv | method, system, and one or more computer readable storage media |
| US20160131800A1 (en) * | 2014-11-07 | 2016-05-12 | Schlumberger Technology Corporation | Modeling fluid-conducting fractures in reservoir simulation grids |
| US10295683B2 (en) * | 2016-01-05 | 2019-05-21 | Schlumberger Technology Corporation | Amplitude inversion on partitioned depth image gathers using point spread functions |
-
2022
- 2022-07-27 CN CN202210889256.7A patent/CN115201913B/en active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101915938A (en) * | 2010-07-05 | 2010-12-15 | 中国科学院地质与地球物理研究所 | Offset imaging method and device for converted waves |
| CN104216012A (en) * | 2014-08-27 | 2014-12-17 | 中国科学院地质与地球物理研究所 | Three-dimensional Born-Kirchhoff variable-step interpolation imaging method |
| CN105717547A (en) * | 2015-12-22 | 2016-06-29 | 吉林大学 | Anisotropy medium magnetotelluric meshless value simulation method |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
| CN112384937A (en) * | 2018-05-12 | 2021-02-19 | 地质探索系统公司 | Seismic data interpretation system |
| CN110888166A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Least square offset imaging method and device based on L-BFGS algorithm |
| CN109946741A (en) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | Pure qP wave least square reverse-time migration imaging method in a kind of TTI medium |
| CN110133713A (en) * | 2019-04-24 | 2019-08-16 | 中国石油大学(华东) | A kind of the multiple wave least square reverse-time migration imaging method and system of full propagation path attenuation compensation |
| CN112649859A (en) * | 2019-10-12 | 2021-04-13 | 中国石油化工股份有限公司 | Seismic wave speed self-adaptive grid-free field node establishment method and system |
| CN110824558A (en) * | 2019-11-20 | 2020-02-21 | 中国石油大学(华东) | Seismic wave numerical simulation method |
| CN114357831A (en) * | 2021-12-29 | 2022-04-15 | 中海石油(中国)有限公司 | Non-grid generalized finite difference forward modeling method, device, storage medium and equipment |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115201913A (en) | 2022-10-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Rawlinson et al. | Seismic ray tracing and wavefront tracking in laterally heterogeneous media | |
| US5539704A (en) | Bayesian sequential Gaussian simulation of lithology with non-linear data | |
| Martinez et al. | Denoising of gravity gradient data using an equivalent source technique | |
| CN105353405B (en) | A kind of full waveform inversion method and system | |
| EA032186B1 (en) | Seismic adaptive focusing | |
| CN115201913B (en) | Least square reverse time migration imaging method, system and storage medium based on gridless finite difference method | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| US11119233B2 (en) | Method for estimating elastic parameters of subsoil | |
| AU739128B2 (en) | A method of seismic processing, and in particular a 3D seismic prospection method implementing seismic data migration | |
| CN115166827B (en) | Least square offset imaging method, device and storage medium based on deconvolution imaging condition | |
| CN114167511B (en) | Bit field data rapid inversion method based on continuous expansion downward continuation | |
| CN113608262B (en) | Seismic data processing method and device for calculating rotation component by using translation component | |
| CN112630840B (en) | Random inversion method based on statistical characteristic parameters and processor | |
| CN115586575A (en) | 3D passive source pseudo-shot gather computing method, system, equipment and storage medium | |
| CN111665549A (en) | Inversion method of stratum acoustic wave attenuation factor | |
| CN111665550A (en) | Underground medium density information inversion method | |
| CN115951402B (en) | Offset imaging method, system and medium based on vector reflectivity forward operator | |
| CN113219534B (en) | Pre-stack depth migration speed quality control method and device, medium and electronic equipment | |
| Taillandier et al. | Refraction traveltime tomography based on adjoint state techniques | |
| CN114002742B (en) | Euler Gaussian beam offset imaging method and device | |
| CN113219533B (en) | Pre-stack time migration velocity modeling method and device, medium and electronic equipment | |
| US11397271B2 (en) | Systems and methods to correct seismic data for undesired effects caused by near surface velocity variations by correcting for long-wavelength statics | |
| WO2024051834A1 (en) | Full-waveform inversion method and device, and storage medium | |
| CN111665547A (en) | Inversion method for formation acoustic wave impedance information | |
| CN110858000B (en) | Seismic data reconstruction method and device, computer equipment and storage medium |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |