CN115201913A - Least squares reverse time migration imaging method, system and storage medium based on meshless finite difference method - Google Patents

Abstract

The invention relates to the technical field of geophysical exploration, in particular to a least square reverse time migration imaging method and system based on a non-grid finite difference method and a storage medium. According to the method provided by the invention, the grid-free distribution generated according to the target work area underground model is used, the grid-free finite difference method reverse time migration is used for obtaining the initial imaging result, and the least square reverse time migration imaging method is carried out by using the migration and reverse migration operator based on the grid-free finite difference method, so that the accurate description and the efficient imaging of the complex earth surface and the underground irregular structure are realized, and the underground complex geological structure can be finely described under any complex earth surface condition. Meanwhile, the method effectively reduces algorithm calculation and memory requirements, improves imaging precision and can provide more accurate imaging results for seismic exploration.

Description

Least squares inverse time migration imaging method and system based on meshless finite difference method and storage media

technical field

The invention relates to the technical field of geophysical exploration, in particular to a least squares reverse time migration imaging method, system and storage medium based on a gridless finite difference method.

Background technique

With the increasing demand for seismic exploration, the current industrial requirements for imaging technology are also increasing. Among the many migration imaging technologies, the reverse time migration technology based on the two-way wave equation has extremely high imaging accuracy, can describe the spatial position and geometric structure of the underground reflection interface of any ray path, and can accurately explain and implement for the structure. Well location provides reliable information. However, the conventional reverse time migration method has inherent shortcomings such as strong low-frequency interference, uneven illumination and poor amplitude preservation. The least squares reverse time migration method is a wave equation inversion method based on reverse time migration and taking reflection coefficients as parameters, which can improve various defects and deficiencies in reverse time migration and effectively improve seismic wave migration imaging. accuracy.

However, traditional least squares reverse time migration imaging methods mostly use finite difference method with regular grid to simulate seismic waves. The grid cells are usually rectangular, which makes it difficult to accurately describe the location of curved interfaces or irregular model boundaries, and it is also difficult to accurately describe the location of curved interfaces or irregular model boundaries. Regions achieve local refinement. This leads to scattering noise at the irregular boundary when numerically calculating the subsurface wave field, which seriously affects the accuracy of the wave field calculation. In order to solve the problem of any shape abnormal body or irregular terrain, many irregular mesh finite difference methods have been proposed. Although these methods have certain geometric flexibility and can solve the model interface or boundary adaptability problems, they often require complex meshes. The meshing or meshing process increases the cost of numerical simulation. In addition, the finite element and spectral element methods can also solve such problems to a certain extent, but their requirements for calculation and memory are much higher than those of the finite difference method. Applicability of the method is poor. As a new numerical calculation method, the gridless finite difference method can effectively solve the problems existing in the above calculation methods, and has been initially applied in the field of seismic exploration. For the least squares inverse time migration imaging method, the meshless finite difference method can be effectively applied.

At present, the focus of domestic exploration has shifted from the eastern plain area to the western plateau and mountainous area, while most of the western exploration areas are faced with complex surface and underground structures. In order to cope with the ever-increasing demand for exploration, and for the problems of difficulty in accurately describing irregular surface and subsurface interfaces, excessive computation and memory requirements in the current traditional least squares reverse time migration method, there is an urgent need for a method that can accurately describe the An efficient imaging method for complex surface and underground irregular structures, providing high-quality imaging results for geological censuses and energy exploration.

SUMMARY OF THE INVENTION

In view of this, the first object of the present invention is to provide a least squares inverse time migration imaging method based on the meshless finite difference method, through the meshless distribution generated according to the underground model of the target work area, and using The migration and demigration operators of the grid finite difference method enable accurate description and efficient imaging of complex surface and subsurface irregular structures.

Based on the same inventive concept, the second object of the present invention is to provide a least squares inverse time migration imaging system based on the meshless finite difference method.

Based on the same inventive concept, the third object of the present invention is to provide a storage medium.

The first object of the present invention can be realized by adopting the following technical solutions:

A least squares inverse time migration imaging method based on meshless finite difference method, comprising the following steps:

Observing the seismic records, and obtaining the underground velocity model of the target work area according to the seismic records obtained by observation;

According to the underground velocity model of the target work area, a gridless node distribution suitable for the underground velocity model is generated, and the underground velocity model is discretized;

According to the observed seismic records and the discretized underground velocity model, the reverse time migration imaging based on the meshless finite difference method is applied to obtain the initial imaging results;

According to the initial imaging results or the updated imaging results, the forward operator based on the gridless finite difference method is applied to calculate the simulated seismic records;

Calculate the record residuals between the simulated seismic records and the observed seismic records;

According to the recorded residuals, apply the migration operator of the meshless finite difference method to obtain the updated gradient of the imaging results corresponding to the recorded residuals, and then obtain the updated imaging results through the linear inversion solver;

It is judged whether the preset termination condition is met, and if so, the final imaging result is output; otherwise, the recording residuals are continuously calculated and the imaging results are updated, and iteration is performed.

Further, in the underground velocity model of the target work area, the discrete degree of nodes is positively correlated with the magnitude of the velocity field.

Further, generating a network-free node distribution adapted to the underground velocity model according to the underground velocity model of the target work area, including the following steps:

According to the underground velocity model of the target work area, determine the size of the computational domain, and determine the mapping relationship between the particle radius of the gridless distribution and the size of the velocity field;

The node distribution is randomly generated at the bottom of the underground velocity model of the target work area, the density of the node distribution is given according to the preset value, and these node positions are defined as potential node positions;

Select a potential node position with the lowest position as the new valid node position;

Take the new valid node as the center of the circle, and use the particle radius at the valid node as the radius of the circle to determine a circular area, and delete all nodes in the circular area except the valid node;

Mark one nearest potential node from the left and right of the center of the circle, determine the direction of the line connecting the center of the circle and the two potential nodes, and place several new potential nodes equidistantly on the arc between the two directions. ;

Continue to select a potential node position with the lowest position as a new valid node position, and repeat the above steps until the nodes fill the entire computational domain.

Further, according to the observed seismic records and the discretized underground velocity model, the reverse time migration imaging based on the meshless finite difference method is applied to obtain the initial imaging results, including the following steps:

Use radial basis functions to construct linear equations, solve the constructed linear equations, and obtain meshless difference coefficients;

The wave equation propagation operator is constructed according to the gridless difference coefficients;

According to the wave equation propagation operator, the reverse time migration imaging based on the meshless finite difference method is applied to obtain the initial imaging result.

Further, the linear equation system constructed by using radial basis function, the specific form is:

表示拉普拉斯算子，c_i为待求的无网格差分系数。Among them, φ(||XX _i ||) represents the radial basis function, which is only related to the coordinate position, and the subscript i=0, 1, 2...n represents the ith node in the differential template, where n represents the node in the differential template number, X represents the spatial coordinates of the gridless nodes, ||·|| represents the two-norm,

represents the Laplacian operator, and ci is the _gridless difference coefficient to be obtained.

Further, an additional basis function is added to the linear equation system constructed by using the radial basis function, so that the static error can be eliminated and the convergence speed of the error can reach the highest order of the additional polynomial.

Further, the additional basis function added to the linear equation system constructed by using radial basis functions is Taylor monomial, and the specific form of the obtained linear equation system is:

表示X-X₀时的径向基函数。Among them, c _n+1 , c _n+2 and c _n+3 are auxiliary coefficients, A is the radial basis function matrix, c is the difference coefficient matrix, and Δx _i is the horizontal coordinate difference between the node X _i and the node X ₀ , _{Δzi is the vertical coordinate difference between node X i and} node X ₀ ,

represents the radial basis function at XX ₀ .

Further, constructing a wave equation propagation operator according to the grid-free difference coefficients includes the following steps:

对应的差分系数，进而求解波动方程中的空间导数项，具体形式为：Find the Laplacian

The corresponding difference coefficient, and then solve the spatial derivative term in the wave equation, the specific form is:

where v is the velocity, p is the seismic wave field, t is the time, c _i is the difference coefficient, and p _i is the wave field value at the ith node in the difference template;

The second-order central difference scheme is used to estimate the time derivative term in the wave equation, and the specific form is:

Among them, p is the seismic wave field, p ₀ is the background wave field, t is the time, the superscript -1, 0 and 1 represent the previous moment, the current moment and the next moment, respectively, Δt is the time step.

The second object of the present invention can be realized by adopting the following technical solutions:

A least squares reverse time migration imaging system based on gridless finite difference method, comprising:

The velocity model unit is used to observe the seismic data and generate the underground velocity model of the target work area according to the observed seismic data;

The gridless node generation unit is used for generating a gridless node distribution suitable for the underground velocity model of the target working area according to the underground velocity model of the target working area, and discretizing the underground velocity model;

The initial imaging unit is used to obtain the initial imaging results by applying reverse time migration imaging based on the gridless finite difference method according to the seismic records obtained by observation and the discretized underground velocity model;

Least squares migration imaging unit for imaging using the principle of least squares migration imaging, including:

According to the initial imaging results or the updated imaging results, the forward operator based on the gridless finite difference method is applied to calculate the simulated seismic records;

Calculate the record residuals between the simulated seismic records and the observed seismic records;

According to the recorded residuals, apply the migration operator of the meshless finite difference method to obtain the updated gradient of the imaging results, and then obtain the updated imaging results through the linear inversion solver;

It is judged whether the preset termination condition is met, and if so, the final imaging result is output; otherwise, the recording residuals are continuously calculated and the imaging results are updated, and iteration is performed.

The third object of the present invention can be realized by adopting the following technical solutions:

A storage medium stores a program, and when the program is executed by a processor, realizes the above-mentioned least squares inverse time migration imaging method based on the meshless finite difference method.

The present invention has the following beneficial effects with respect to the prior art:

(1) The meshless node generation algorithm in the least squares inverse time migration imaging method based on the meshless finite difference method provided by the present invention can generate curved interfaces and irregular boundaries suitable for any model, and the node distribution density can be Based on the gridless node distribution preset by the model speed, the generation process is simple and efficient, and the algorithm can be guaranteed to have extremely high geometric flexibility without high computational cost, and can effectively reduce the memory requirement of the algorithm.

(2) The calculation method of meshless finite difference coefficient and wave equation operator construction method ensures the accuracy of wave field calculation under meshless node distribution. "The influence of scattering noise introduced by the approximation provides a convenient condition for the implementation of the high-precision least squares inverse time migration imaging algorithm, and helps to speed up the convergence efficiency of the imaging algorithm.

(3) Generally speaking, the method provided by the present invention can accurately describe the complex underground geological structure under any complex surface conditions. Compared with the traditional method, the method can effectively reduce the algorithm calculation and memory requirements while improving the imaging performance. It can provide more accurate imaging results for seismic exploration.

Description of drawings

FIG. 1 is a schematic flowchart of a least squares inverse time migration imaging method based on a meshless finite difference method provided in Embodiment 1 of the present invention.

FIG. 2 is a schematic diagram comparing the regular grid distribution and the gridless distribution provided in Embodiment 1 of the present invention.

FIG. 3 is a schematic diagram of a fast generation algorithm for gridless node distribution provided in 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 Schematic.

FIG. 4 is a schematic diagram of an underground velocity model of a target work area with three sets of strata in Embodiment 1 of the present invention.

FIG. 5 is a schematic diagram of the number of nodes required for different node discrete manners when performing migration imaging on the three-layer velocity model in Embodiment 1 of the present invention.

6 is a schematic diagram of an imaging result obtained by performing meshless finite difference least squares reverse time migration on 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 traditional regular grid-based finite difference least squares inverse time migration on the three-layer velocity model in Embodiment 1 of the present invention.

8 is a schematic diagram of an underground velocity model of a target construction area with a complex surface in Embodiment 1 of the present invention.

9 is a schematic diagram of an imaging result obtained by performing conventional reverse time migration based on gridless finite difference on the complex surface velocity model in Embodiment 1 of the present invention.

10 is a schematic diagram of an imaging result obtained by performing gridless finite difference least squares inverse time migration on the complex surface velocity model in Embodiment 1 of the present invention.

FIG. 11 is a schematic diagram of a least squares inverse time migration imaging system based on a meshless finite difference method in Embodiment 2 of the present invention.

Detailed ways

In order to make the purposes, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in 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. Obviously, the described embodiments It is a part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative work are protected by the present invention. scope.

Example 1:

As shown in FIG. 1 , this embodiment provides a least squares inverse time migration imaging method based on a meshless finite difference method, including the following steps:

S100, observe the seismic records, and obtain the underground velocity model of the target work area according to the observed seismic records, including the following steps:

S110. Observing seismic records. In this embodiment, the observed seismic records are the excitation of artificial explosive sources (explosives in the onshore work area, and air guns in the offshore work area), and active source seismic data received by a high-density geophone array. In order to improve the quality of seismic imaging and the effect of record reconstruction, the observed seismic data can be denoised (removal of random noise, surface wave noise, radio interference, adjacent shot interference, linear noise, etc.), removal (removal of refracted wave components in seismic records) , deconvolution, etc.

S120, obtaining an underground velocity model of the target work area according to the seismic records obtained by observation. In this embodiment, according to the seismic records obtained by observation, a relatively accurate underground velocity model of the target work area is obtained by velocity modeling means such as tomography or full waveform inversion.

S200 , as shown in FIG. 3 , according to the underground velocity model of the target work area, a gridless node distribution suitable for the underground velocity model is generated, and the underground velocity model is discretized.

In this step, the meshless node distribution is generated to discretize the underground velocity model. After the meshless node distribution is generated, the underground wave field will be numerically solved discretely according to these nodes. This step includes the following steps:

S210. Determine the size of the calculation domain according to the underground velocity model of the target work area, and determine the mapping relationship between the particle radius of the gridless distribution and the size of the velocity field; The larger the value, the larger the corresponding particle radius;

S220, as shown in subgraph (a) of FIG. 3, randomly generate node distribution at the bottom of the underground velocity model of the target work area, the density of the node distribution is given according to the preset value, and define these node positions as potential node positions (PDP);

S230, selecting a potential node position with the lowest position as a new valid node position;

S240. As shown in subgraphs (b) and (c) of FIG. 3, take the new effective node as the center of the circle, and use the particle radius at the effective node as the radius of the circle to determine a circular area, delete the inner part of the circular area All nodes except valid nodes;

S250. As shown in the sub-graph (d) of FIG. 3, mark a nearest potential node from the left and right of the circle center, determine the direction of the line connecting the circle center and the two potential nodes, and locate the two potential nodes between the two directions. 5 new potential nodes are placed equidistantly on the arc of ;

S260. Continue to select a potential node position with the lowest position as a new valid node position, and repeat steps S240-S260 until the nodes fill the entire computational domain.

S300, according to the seismic records obtained by observation and the discretized underground velocity model, apply reverse time migration imaging based on the gridless finite difference method to obtain initial imaging results, including the following steps:

S310 , constructing a linear equation system by using a radial basis function, and solving the constructed linear equation system to obtain a meshless difference coefficient.

In this embodiment, the specific form of the linear equation system constructed by using the radial basis function is shown in the following formula:

表示拉普拉斯算子，c_i为待求的无网格差分系数。The above system of linear equations can be abbreviated as Ac=LΦ. Among them, φ(||XX _i ||) represents the radial basis function, which is a function only related to the coordinate position, and the subscript i=0,1,2...n represents the ith node in the difference template, where n represents the difference The number of nodes in the template, X represents the spatial coordinates of the meshless nodes generated in step S200, ||·|| represents the two-norm,

represents the Laplacian operator, and ci is the _gridless difference coefficient to be obtained.

Additional basis functions (generally polynomials) can be added to the above linear equation system constructed by using radial basis functions, so that the static error can be eliminated and the error convergence speed can reach the highest order of the additional polynomial. In this embodiment, the Taylor monomial is added to the above linear equation system to obtain the linear equation system shown in the following formula:

表示X-X₀时的径向基函数。Among them, c _n+1 , c _n+2 and c _n+3 are auxiliary coefficients, A is the radial basis function matrix, c is the difference coefficient matrix, and Δx _i is the horizontal coordinate difference between the node X _i and the node X ₀ , _{Δzi is the vertical coordinate difference between node X i and} node X ₀ ,

represents the radial basis function at XX ₀ .

S320, constructing a wave equation propagation operator according to the gridless difference coefficient, including the following steps:

对应的差分系数，进而求解波动方程中的空间导数项，具体形式如下式所示：S321. Obtain Laplace operator

The corresponding difference coefficient, and then solve the spatial derivative term in the wave equation, the specific form is shown in the following formula:

where v is the velocity, p is the seismic wave field, t is the time, c _i is the difference coefficient, and p _i is the wave field value at the ith node.

S322, the second-order central difference format is used to estimate the time derivative term in the wave equation, and the specific form is shown in the following formula:

Among them, p is the seismic wave field, t is the time, the superscript -1, 0 and 1 represent the previous moment, the current moment and the next moment, respectively, p ₀ is the background wave field, Δt is the time step.

S330, according to the wave equation propagation operator, applying the reverse time migration imaging based on the gridless finite difference method to obtain the initial imaging result, including the following steps:

S331. Input the seismic records obtained by the observation through the wave equation propagation operator, and calculate the source wave field and the detection wave field. The specific form is shown in the following formula:

where p _s and _pr are the source wave field and receiver wave field, respectively, x _s and x _r are the spatial positions of the source and receiver points, respectively, S is the source function, and d _obs is the observed seismic record.

S332, applying imaging conditions to the seismic source and the detector wave field. In this embodiment, the applied imaging conditions are cross-correlation imaging conditions, and the specific form is shown in the following formula:

I(x)= _∫∫ps (x,xs, _{t)pr(x,xr , t} ) _dtdxs

where I represents the imaging result.

S400 , according to the initial imaging result or the updated imaging result, apply the forward operator based on the meshless finite difference method to calculate and obtain the simulated seismic record, and the specific form is shown in the following formula:

where p ₀ and δp are the background wave field and disturbance wave field, respectively, and d is the simulated seismic record.

S500. Calculate the recording residual between the simulated seismic record and the observed seismic record. In this embodiment, the simulated seismic record obtained in step S400 is compared with the observed seismic record in step S110 to obtain the record residual. .

S600. Apply the offset operator of the gridless finite difference method according to the recorded residual to obtain the updated gradient of the imaging result corresponding to the recorded residual, and then obtain the updated imaging result through the linear inversion solver, including the following steps :

S610. Input the recorded residual, and use the offset operator based on the gridless finite difference method to calculate the update gradient of the imaging result, and the specific form is shown in the following formula:

Where g is the update gradient of the imaging result, λ is the accompanying wave field of the background wave field p ₀ , and v ₀ is the background velocity field.

S620. Obtain the updated imaging result through an optimization method, such as the steepest descent method, the conjugate gradient method, and the like.

S700. Determine whether the preset termination condition is met, and if so, output the final imaging result; otherwise, continue to calculate the recording residual and update the imaging result, and perform iteration.

In this embodiment, the iterative method is to input the new imaging result obtained in step S620 as the imaging result of step S400, and perform steps S400-S700 again.

In this embodiment, the preset termination condition is: satisfying any one of the preset conditions, including:

The first preset condition, the objective function related to the residual is less than the set threshold;

The second preset condition is that the number of algorithm iterations is greater than the set threshold.

Figure 2 is a schematic diagram of the comparison between the regular grid distribution and the gridless distribution. It can be seen that for the irregular interface (curved solid line) in the computational domain, when the regular grid distribution is discrete, there will be a "stepped" approximation (black dotted line). ), and through the meshless distribution discretization, the node distribution can fit the interface well, and can effectively avoid the scattering noise introduced by the "staircase" approximation in the numerical calculation of the subsurface wave field. Figure 3 is a schematic diagram of a fast generation algorithm for gridless node distribution. According to this algorithm, the underground velocity model can be rapidly discretized through gridless nodes.

Figure 4 shows the underground velocity model of the target work area with three sets of strata adopted in this embodiment, and its size is 3km×2km. The number of nodes required, it can be seen that without grid, only 2/3 of the number of nodes in the regular grid can be discretized to discretize the velocity model, which means that when the imaging algorithm is calculated, the calculation and memory requirements required by the algorithm are also corresponding. reduce.

FIG. 6 and FIG. 7 are respectively the imaging results obtained by using the gridless finite difference-based least squares inverse time migration imaging method and the conventional regular grid finite difference-based least squares inverse time migration imaging method in this embodiment. . It can be seen that both can describe the subsurface structure and provide amplitude information similar to the true reflectivity model. However, due to the step-like approximation of the regular grid, the continuity of the events in the imaging results in Fig. 7 is poor, which may be more serious and even affect the algorithm convergence when imaging more complex subsurface structures. The imaging results in Figure 5 can almost perfectly describe the subsurface structure. By comparing the imaging results in FIG. 6 and FIG. 7 , it can be found that the least squares inverse time migration imaging method based on gridless finite difference in this embodiment can provide more accurate imaging methods on the premise of effectively reducing algorithm calculation and memory requirements. Subsurface imaging results.

FIG. 8 is an underground velocity model of a target construction area with a complex surface adopted in this embodiment, and its size is 4.16km×2.49km. FIG. 9 and FIG. 10 are respectively the imaging results obtained by using the conventional inverse time migration imaging method based on gridless finite difference and the least squares inverse time migration imaging method in the embodiment of the present invention. It can be seen that the imaging accuracy at the surface is very high in the imaging results. Compared with the imaging results in Figure 9, the imaging results in Figure 10 have significantly higher resolution, which can provide more reliable reflectivity for structural and lithologic interpretation. Model. This also verifies the adaptability of the gridless finite difference-based least squares reverse time migration imaging method in the embodiment of the present invention to complex surface conditions, and demonstrates that the method provides high-precision imaging for seismic exploration under complex geological conditions ability to result.

Therefore, applying the meshless node generation algorithm in the least squares inverse time migration imaging method based on the meshless finite difference method provided in this embodiment can generate curved interfaces and irregular boundaries suitable for any model, and the node distribution density can be based on The gridless 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 spending a high computational cost, and the memory requirements of the algorithm can be effectively reduced; at the same time, this embodiment The gridless finite difference coefficient calculation and wave equation operator construction method adopted ensures the accuracy of wavefield calculation under gridless node distribution. At the same time, because the wavefield simulation in this method is not subject to the "staircase" at the irregular boundary The influence of scattering noise introduced by the approximation provides a convenient condition for the implementation of the high-precision least squares inverse time migration imaging algorithm, and helps to speed up the convergence efficiency of the imaging algorithm; in general, the method provided in this embodiment can be used in any complex surface Compared with the traditional method, this method can effectively reduce the calculation and memory requirements of the algorithm, and at the same time improve the imaging accuracy, and can provide more accurate imaging results for seismic exploration.

Example 2:

As shown in FIG. 11 , this embodiment provides a least squares inverse time migration imaging system based on the gridless finite difference method, including:

The velocity model unit 1101 is used to observe the seismic data, and generate an underground velocity model of the target work area according to the observed seismic data;

The meshless node generation unit 1102 generates a meshless node distribution suitable for the underground velocity model according to the underground velocity model of the target work area, and discretizes the underground velocity model;

The initial imaging unit 1103, according to the seismic records obtained by observation and the discretized underground velocity model, applies reverse time migration imaging based on the gridless finite difference method to obtain initial imaging results;

The least squares migration imaging unit 1104 is configured to perform imaging using the least squares migration imaging principle, including:

According to the imaging results, the forward operator based on the meshless finite difference method is applied to obtain the simulated seismic records;

Calculate the record residuals between the simulated seismic records and the observed seismic records;

According to the residual, apply the migration operator of the meshless finite difference method to obtain the updated gradient of the imaging result, and then obtain the updated imaging result through the linear inversion solver;

It is judged whether the preset termination condition is met, and if so, the imaging result is output; otherwise, the residual error and the imaging result are continued to be calculated, and the iteration is performed.

Example 3:

This embodiment provides a storage medium, where the storage medium is a computer-readable storage medium and stores a computer program. When the computer program is executed by a processor, the gridless finite difference method based on the above-mentioned Embodiment 1 is implemented. The least squares inverse time migration imaging method is as follows:

Observing the seismic records, and obtaining the underground velocity model of the target work area according to the seismic records obtained by observation;

According to the underground velocity model of the target work area, a gridless node distribution suitable for the underground velocity model is generated, and the underground velocity model is discretized;

According to the observed seismic records and the discretized underground velocity model, the reverse time migration imaging based on the meshless finite difference method is applied to obtain the initial imaging results;

According to the initial imaging results or the updated imaging results, the forward operator based on the gridless finite difference method is applied to calculate the simulated seismic records;

Calculate the record residuals between the simulated seismic records and the observed seismic records;

According to the recorded residuals, apply the migration operator of the meshless finite difference method to obtain the updated gradient of the imaging results corresponding to the recorded residuals, and then obtain the updated imaging results through the linear inversion solver;

It is judged whether the preset termination condition is met, and if so, the final imaging result is output; otherwise, the recording residuals are continuously calculated and the imaging results are updated, and iteration is performed.

It should be noted that the computer-readable storage medium in this embodiment may be a computer-readable signal medium or a computer-readable storage medium, or any combination of the above two. The computer-readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus or device, or a combination of any of the above. More specific examples of computer readable storage media may include, but are not limited to, electrical connections with one or more wires, portable computer disks, hard disks, random access memory (RAM), read only memory (ROM), erasable Programmable read only memory (EPROM or flash memory), fiber optics, portable compact disk read only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.

In this embodiment, the computer-readable storage medium can be any tangible medium that contains or stores a program that can be used by or in conjunction with an instruction execution system, apparatus, or device. In this embodiment, however, the computer-readable signal medium may include a data signal in baseband or propagated as part of a carrier wave, carrying a computer-readable program therein. Such propagated data signals may take a variety of forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination of the foregoing. A computer-readable signal medium can also be any computer-readable storage medium, other than a computer-readable storage medium, that can be sent, propagated, or transmitted for use by or in connection with the instruction execution system, apparatus, or device. program. A computer program embodied on a computer-readable storage medium may be transmitted using any suitable medium including, but not limited to, electrical wire, optical fiber cable, RF (radio frequency), etc., or any suitable combination of the foregoing.

The above-mentioned computer-readable storage medium can write a computer program for executing the present embodiment in one or more programming languages or a combination thereof, the above-mentioned programming languages include object-oriented programming languages—such as Java, Python, C++, Also included are conventional procedural programming languages - such as C or similar programming languages. 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. Where a remote computer is involved, the remote computer may be connected to the user's computer through any kind of network, including a local area network (LAN) or wide area network (WAN), or may be connected to an external computer (eg, using an Internet service provider to connect through the Internet) ).

Obviously, the above-mentioned embodiments are only a part of the embodiments of the present invention, rather than all the embodiments, and the present invention is not limited to the details of the above-mentioned embodiments, any suitable changes or modifications made by those of ordinary skill in the art, All are regarded as not departing from the patent scope of the present invention.

Claims (10)

1. A least square reverse time migration imaging method based on a meshless finite difference method is characterized by comprising the following steps:

observing the seismic records, and acquiring an underground velocity model of the target work area according to the seismic records obtained by observation;

generating non-grid node distribution suitable for an underground speed model according to the underground speed model of the target work area, and discretizing the underground speed model;

according to the seismic record obtained by observation and the discretized underground velocity model, applying reverse time migration imaging based on a non-grid finite difference method to obtain an initial imaging result;

calculating to obtain a simulated seismic record by applying a forward operator based on a non-grid finite difference method according to the initial imaging result or the updated imaging result;

calculating a recording residual error between the simulated seismic record and the seismic record obtained by observation;

according to the recorded residual error, applying the recorded residual error to an offset operator of a non-grid finite difference method to obtain an updated gradient of an imaging result corresponding to the recorded residual error, and obtaining an updated imaging result through a linear inversion solver;

judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuously calculating and recording residual errors and updating the imaging result, and performing iteration.

2. The grid-free finite difference method-based least square reverse time migration imaging method according to claim 1, wherein in the target work area underground velocity model, the discrete degree of the nodes is positively correlated with the size of the velocity field.

3. The least square reverse time migration imaging method based on the meshless finite difference method as claimed in claim 1, wherein a meshless node distribution adapted to the underground velocity model is generated according to the target work area underground 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 non-grid 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 out the density of the node distribution according to a preset value, and defining the positions of the nodes 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 the new effective node as the center of a circle and the particle radius 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 respectively, determining the direction of a connecting line of the circle center and the two potential nodes, and placing a plurality of new potential nodes on an arc clamped by the two directions at equal intervals;

and continuously selecting a potential node position with the lowest position as a new effective node position, and repeating the steps until the nodes fill the whole calculation domain.

4. The method of claim 1, wherein a reverse time migration imaging based on the meshless finite difference method is applied according to the observed seismic record and the discretized subsurface velocity model to obtain an initial imaging result, comprising the steps of:

constructing a linear equation set by using the radial basis function, and solving the constructed linear equation set to obtain a non-grid differential coefficient;

constructing a wave equation propagation operator according to the non-grid difference coefficient;

and acquiring the initial imaging result by applying reverse time migration imaging based on a non-grid finite difference method according to a wave equation propagation operator.

5. The grid-free finite difference method-based least square reverse time migration imaging method according to claim 4, wherein the linear equation set constructed by using the radial basis function is in a specific form:

wherein, phi (| | X-X) _i I |) represents the radial basis function, only in relation to the coordinate position, subscript i =0,1,2 \8230nrepresents the ith node in the difference template, where n represents the number of nodes in the difference template, X represents the space coordinate without grid nodes, | | | · | represents the two-norm,

represents the Laplace operator, c _i The non-grid difference coefficient to be solved.

6. The method of claim 5, wherein an additional basis function is added to the linear system of equations constructed using radial basis functions to eliminate static errors and to increase the convergence rate of errors to the highest order of the additional polynomial.

7. The method of claim 6, wherein the additional basis functions added to the linear system of equations constructed using radial basis functions are Taylor monomials, and the linear system of equations obtained by the method has the specific form:

wherein, c _n+1 、c _n+2 And c _n+3 For auxiliary coefficients, A is a matrix of radial basis functions, c is a matrix of difference coefficients, Δ x _i Is X _i Node and X ₀ Difference in horizontal coordinate of node, Δ z _i Is X _i Node and X ₀ Vertical coordinate difference of node, phi (| X-X) ₀ ||)|x ₀ Represents X-X ₀ Radial basis function of time.

8. The grid-free finite difference method-based least square reverse time migration imaging method according to claim 4, wherein a wave equation propagation operator is constructed according to the grid-free difference coefficient, comprising the following steps:

calculating laplacian

And solving a spatial derivative term in the wave equation according to the corresponding difference coefficient, wherein the specific form is as follows:

where v is velocity, p is the seismic wavefield, t is time, c _i Is a difference coefficient, p _i The wave field value at the ith node in the differential template;

estimating a time derivative term in the wave equation by adopting a second-order central difference format, wherein the specific form is as follows:

where p is the seismic wavefield and p is ₀ For the background wavefield, t is time, and the superscripts-1, 0, and 1 represent the previous time, the current time, and the next time, respectively.

9. A least square reverse time migration imaging system based on a non-grid finite difference method is characterized by comprising the following steps:

the velocity model unit is used for observing seismic data and generating an underground velocity model of the target work area according to the observed seismic data;

the non-grid node generation unit is used for generating non-grid node distribution which is suitable for the target work area underground speed model according to the target work area underground speed model and discretizing the underground speed model;

the initial imaging unit is used for applying reverse time migration imaging based on a non-grid finite difference method according to the seismic record obtained by observation and the discretized underground velocity model to obtain an initial imaging result;

a least squares offset imaging unit for imaging using the least squares offset imaging principle, comprising:

calculating to obtain a simulated seismic record by applying a forward operator based on a non-grid finite difference method according to the initial imaging result or the updated imaging result;

calculating a recording residual error between the simulated seismic record and the seismic record obtained by observation;

according to the recorded residual error, applying the recorded residual error to an offset operator of a non-grid finite difference method to obtain an updated gradient of an imaging result corresponding to the recorded residual error, and obtaining an updated imaging result through a linear inversion solver;

judging whether a preset termination condition is met, and if so, outputting a final imaging result; otherwise, continuously calculating and recording residual errors and updating the imaging result, and performing iteration.

10. A storage medium storing a program which, when executed by a processor, realizes the least-squares reverse-time migration imaging method based on the meshless finite difference method according to any one of claims 1 to 8.

Images