Variable grid constrained Gaussian beam target imaging method

Abstract

The invention belongs to the technical field of seismic exploration, and relates to a variable grid constrained Gaussian beam target imaging method, which comprises the following steps of: firstly, acquiring seismic data of a target, observation data of an observation system and offset parameters of Gaussian beam offset for a set imaging target; then establishing a migration velocity field model, and setting values of an average sampling rate and a highest sampling rate; performing first global mesh generation according to a mesh generation formula on the basis of the speed difference in the vertical direction of the migration speed field model; carrying out second mesh generation according to the characteristics of different target imaging bodies, and refining a target layer mesh; and finally, a reverse continuation wave field is constructed by adopting an imaging point uplink ray tracing technology, and target-oriented high-precision imaging is realized. On the basis of Gaussian beam migration imaging, the imaging precision is further improved.

Description

Technical Field

[0001] This invention belongs to the field of seismic exploration technology and relates to a Gaussian beam target imaging method with variable grid constraints. Background Technology

[0002] Migration imaging is a seismic data processing technique that uses geophysical theory to backpropagate the wavefield of seismic records observed at the Earth's surface. After eliminating the propagation effects of seismic waves, it yields images of the subsurface structure. For decades, migration imaging has been widely used in geophysical exploration, geological resource investigation, and engineering fields such as geotechnical engineering, bridges, and tunnels. As an indispensable and crucial step in the seismic data processing workflow, migration imaging's accuracy not only plays a vital role in seismic data interpretation but also provides significant guidance for oil and gas exploration and development.

[0003] However, migration imaging is a complex mathematical and physical problem that requires a comprehensive consideration of computational accuracy and efficiency. Simply pursuing high accuracy or solely focusing on high efficiency is not advisable. Currently, migration imaging methods are mainly divided into two types: ray-based and wave equation-based. Ray-based methods primarily employ the ray integral method, which has higher computational efficiency, but its computational accuracy is far inferior to the differential method commonly used in wave equation-based methods. The differential method requires more memory and computation time, resulting in relatively lower computational efficiency. Therefore, considering practical considerations, ray-based migration, specifically ray integral migration, remains the most widely used imaging method in actual production.

[0004] Currently, with the continuous deepening of oil and gas resource exploration, the target of seismic exploration has gradually shifted from large-scale structural oil and gas reservoirs to medium- and small-scale complex lithological oil and gas reservoirs. This places higher demands on the accuracy of exploration data, and conventional seismic exploration techniques are no longer sufficient to meet the current needs of oil and gas resource exploration. Therefore, how to maintain the high efficiency of ray migration while improving imaging accuracy has become a research focus for both academia and industry.

[0005] Currently, the most studied ray migration imaging methods are Kirchhoff migration and Gaussian beam migration. Kirchhoff migration suffers from singularities in the caustic region and the absence of a wavefield in the shadow region, which degrades imaging quality. Gaussian beam migration addresses these issues by using a complex exponential term with a Gaussian envelope decay relative to the central ray. Furthermore, by incorporating a side-axis ray compared to asymptotic rays, it considers more propagation information and is not limited by multipath propagation problems. Therefore, with its higher imaging accuracy and comparable efficiency, Gaussian beam migration is often considered an optimized alternative to Kirchhoff migration. Maintaining high imaging efficiency while further improving imaging accuracy remains a crucial research direction. Summary of the Invention

[0006] The purpose of this invention is to provide a Gaussian beam target imaging method with variable grid constraints. Based on Gaussian beam migration imaging, it further improves imaging accuracy and promotes the progress of research on high-efficiency and high-precision imaging methods in current seismic exploration technology.

[0007] The technical solution adopted in this invention is a Gaussian beam target imaging method with variable mesh constraints, comprising the following steps: S1: For the set imaging target, acquire the target's seismic data, the observation data of the observation system, and the migration parameters of the Gaussian beam migration; S2: Based on the data and parameters collected in S1, establish the offset velocity field model and set the average sampling rate. and highest sampling rate The value; S3: Based on the vertical velocity difference of the offset velocity field model, the first global mesh of the imaging domain is performed according to the meshing formula; S4: Perform a second mesh subdivision of the imaging domain based on the characteristics of different target images to refine the mesh of the imaging target layer; S5: Employs up-ray tracing technology for imaging points to construct a reverse extended wavefield, achieving high-precision imaging of the target.

[0008] The invention is further characterized in that: the offset parameter in S1 includes the lateral sampling points of the offset velocity field. and longitudinal sampling points Spatial sampling interval Time sampling interval Number of time sampling points , main frequency Reference frequency highest frequency , number of offset shots in earthquake records and the angle of ray emission and The parameters.

[0009] The first global mesh generation of the imaging domain in S3 includes the following steps: S3.1: Perform a global scan analysis on the initial offset velocity field and calculate the average velocity. and minimum speed Simultaneously, based on the velocity change gradient, velocity layers are divided and the minimum velocity arrangement of each layer is calculated. ; S3.2: Based on the calculation results of S3.1 and the average sampling rate set in S2. and highest sampling rate The first global mesh partitioning of the imaging domain is performed using the following formula 1. (1) in, This indicates the grid scale arrangement.

[0010] The method for dividing velocity layers based on velocity change gradients in S3.1 is as follows: Velocity layers are defined by the rate of change of velocity between adjacent sampling points; If the rate of change of velocity at a sampling point changes drastically, that is, the rate of change of velocity at that sampling point is greater than 20% of the average rate of change of velocity at the two adjacent sampling points, then this point is considered to be the interface of a layer.

[0011] The method for performing the second mesh subdivision of the imaging domain in S4 is as follows: According to each layer The secondary grid size is determined to optimize the grid scale of the target layer. The grid size is then adjusted based on the Gaussian beam migration imaging results from the first grid division to ensure fine imaging. One method for adjusting the mesh size is to adjust the speed of a certain layer. When the velocity of a layer is simultaneously greater than or simultaneously less than the velocity of its adjacent layers, the imaging domain of that layer is subdivided into grids, adjusted from 5x5 to 2x2 or 1x1.

[0012] The present invention is further characterized in that S5 includes the following steps: S5.1: Represent the Green's function using Gaussian beam integrals; S5.2: Representing the uplink and downlink wave fields using Green's function; S5.3: Perform correlation imaging of the uplink and downlink wave fields under the constraint of a Gaussian window; S5.4: Based on the grid size of the secondary subdivision, interpolate the Green's function from the shot point and receiver point at the imaging point calculated in S5.2, so that there is a corresponding Green's function at each grid point of the secondary subdivision. Repeat the uplink and downlink wavefield of S5.3 under the constraint of Gaussian window to obtain the Gaussian beam migration imaging result with variable grid constraint.

[0013] The beneficial effects of this invention are as follows: This invention provides an imaging method that combines high efficiency and high precision, improving the overall efficiency and target layer accuracy of beam migration imaging. Based on Gaussian beam migration imaging, this invention employs a variable mesh partitioning method for the first global mesh partitioning of the imaging domain, followed by a second mesh partitioning tailored to the characteristics of different target images, refining the target layer mesh. Finally, it utilizes up-ray tracing technology at imaging points to construct a reverse-extended wavefield, achieving high-precision target-oriented imaging. The imaging method proposed in this invention further improves imaging accuracy, promotes advancements in research on high-efficiency, high-precision imaging methods in current seismic exploration technology, better meets the needs of current oil and gas resource exploration, and has higher practical application value. Attached Figure Description

[0014] Figure 1 This is a flowchart of the present invention; Figure 2 This is a schematic diagram of the three-dimensional ray center coordinate system in step 5 of the present invention; Figure 3 This is a schematic diagram of the upward ray tracing imaging in step 5 of the present invention; Figure 4 The velocity field and the offset results of the three methods are shown in Example 3 of this invention; Figure 5 The offset velocity field and mesh subdivision curves are shown in Example 3 of this invention; Figure 6 The memory usage and computation time are illustrated in Example 3 of this invention. Figure 7 The velocity field offset is an example of the velocity field in Embodiment 4 of this invention; Figure 8 The offset results are for the two methods exemplified in Embodiment 4 of this invention; Figure 9 The images show the imaging results of the two methods at the target location in Example 4 of this invention. Detailed Implementation

[0015] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0016] Example 1: like Figure 1 As shown, the variable mesh-constrained Gaussian beam target imaging method includes the following steps: S1: For the set imaging target, acquire the target's seismic data, the observation data of the observation system, and the migration parameters of the Gaussian beam migration; The offset parameters include the lateral sampling points of the offset velocity field. and longitudinal sampling points Spatial sampling interval Time sampling interval Number of time sampling points , main frequency Reference frequency highest frequency , number of offset shots in earthquake records and the angle of ray emission and The parameters.

[0017] S2: Based on the data and parameters collected in S1, establish the offset velocity field model and set the average sampling rate. and highest sampling rate The value is used for coarse mesh generation of the imaging domain in subsequent steps; S3: Based on the vertical velocity difference of the offset velocity field model, the first global mesh of the imaging domain is performed according to the mesh generation formula, which includes the following steps; S3.1: Perform a global scan analysis on the initial offset velocity field and calculate the average velocity. and minimum speed Simultaneously, based on the velocity change gradient, velocity layers are divided and the minimum velocity arrangement of each layer is calculated. ; The method for dividing velocity layers is to divide velocity layers by the velocity change rate of adjacent sampling points; if the velocity change rate of a sampling point changes drastically, that is, the velocity change rate at that sampling point is greater than 20% of the average velocity change rate of the two adjacent sampling points, then this point is considered to be the interface of a layer. S3.2: Based on the calculation results of S3.1 and the average sampling rate set in S2. and highest sampling rate The first global mesh partitioning of the imaging domain is performed according to the following formula 1. (1) in, Indicates grid-scale arrangement; S4: Perform a second mesh subdivision of the imaging domain based on the characteristics of different target images to refine the mesh of the imaging target layer. The specific method is as follows: According to each layer The secondary grid size is determined to optimize the grid scale of the target layer. The grid size is then adjusted based on the Gaussian beam migration imaging results from the first grid division to ensure fine imaging. One method for adjusting the mesh size is to adjust the speed of a certain layer. When the velocities of the two layers are simultaneously greater than or simultaneously less than the velocities of the adjacent layers, the imaging domain of that layer is subdivided into grids, adjusted from 5x5 to 2x2 or 1x1. S5: Employing uplink ray tracing technology at imaging points, a reverse-extended wavefield is constructed to achieve high-precision target-oriented imaging. This step is based on the principle of Gaussian beam shift and specifically includes the following steps: S5.1: Represent the Green's function using Gaussian beam integrals; S5.2: Representing the uplink and downlink wave fields using Green's function; S5.3: Perform correlation imaging of the uplink and downlink wave fields under the constraint of a Gaussian window; S5.4: Based on the grid size of the secondary subdivision, interpolate the Green's function from the shot point and receiver point at the imaging point calculated in S5.2, so that there is a corresponding Green's function at each grid point of the secondary subdivision. Repeat the uplink and downlink wave fields of S5.3 under the constraint of Gaussian window to perform correlation imaging and obtain the Gaussian beam migration imaging result with variable grid constraint.

[0018] In this embodiment, the highest sampling rate The setting is to limit the occurrence of extremely coarse grids, in order to ensure stable imaging in high-speed regions.

[0019] In this embodiment, the first global mesh partitioning of the imaging domain adopts a variable mesh partitioning method, the principle of which is as follows: Currently, conventional beam shifting in the imaging domain employs coarse-grid algorithms. Existing methods for coarse-grid imaging domain partitioning involve taking five sampling points in each of the x, y, and z directions of the velocity field. While this significantly improves the efficiency of beam scanning imaging, it also introduces discontinuous, jagged breaks and imaging artifacts, particularly noticeable in regions with rapid velocity changes. Using a globally fine-grid algorithm, while avoiding these distortions, results in a geometric increase in computational cost and data caching.

[0020] On the other hand, considering the correlation between mesh size and background velocity, the time step is constant during ray tracing based on the Runge-Kutta method. In the low-velocity region, two adjacent points on the same ray are very close, and it is necessary to divide them with a fine mesh to reduce the amount of information contained in a single mesh. However, in the high-velocity region, a coarse mesh is sufficient, and mesh refinement would be wasteful. Therefore, the idea of ​​variable meshing can be perfectly integrated with the process of ray tracing.

[0021] In this embodiment, the principle of performing the second mesh subdivision of the imaging domain is as follows: Based on the initial global mesh generation, and considering the potential for drastically changing velocity boundaries, anomalous velocity volumes, and the need for higher-precision imaging of the target layer, a second-level mesh generation was implemented. This second-level mesh generation involves adjusting the grid scale arrangement after the initial mesh generation. Local mesh refinement is performed. Based on the abnormal velocity fluctuations, a secondary mesh size is determined to optimize the mesh scale of the target layer. The mesh size is then adjusted in conjunction with the migration results to ensure fine imaging.

[0022] Example 2: Based on Example 1, the specific details of step S5 in this variable mesh constrained Gaussian beam target imaging method are as follows: S5.1: The Green's function is expressed using Gaussian beam integrals, as follows: The Green's function is essentially the wave field generated by a point source. Due to the independent imaging characteristic of the Gaussian beam method, when solving for the wave field effect at a certain point, it can be decomposed into a superposition of the effects produced by a series of Gaussian beams at that location. Here, a three-dimensional ray center coordinate system is introduced (e.g., Figure 2 The expression for the Gaussian beam (as shown) is given by the following formula 2: (2) in, These are the three base coordinates of the imaging point in the three-dimensional ray center coordinate system, and among them... , is along the coordinate axis in the coordinate system of the ray center. The two components; It is the curve length of the imaging point from the origin. These are the coordinate components of the initial point. The velocity at that point along both coordinate axes The second derivative matrix; It is a solution to the three-dimensional dynamics ray tracing equations. , among them This is another solution to the above system of equations. Indicates travel information. The Gaussian beam form solution representing the seismic wave field; In the Gaussian beam method, the Green's function is expressed by the superposition integral of a series of Gaussian beams emitted from the source point in different directions, as shown in Equation 3 below: (3) in, It is a Green's function expression; These are ray parameter vectors representing the imaging points in three directions; It is a wave field represented by a Gaussian beam; Substituting the above expression for the Gaussian beam into Formula 3, we obtain the following Formula 4: (4) in, It is a vector of ray parameters at the surface receiving point, representing three directions; The Green's function represents the contribution of the image emitted from the receiving point on the ground to the imaging point; The coordinates representing the position of the bundle center; Here, we introduce The exponent term is used to transform the Green's function from a point... Extending to a range from its beam center point, this is because the Gaussian beam has a certain width, and the detector point... It can receive energy from different beams, which can be seen as a weighted superposition of them; S5.2: The up and down wave fields are represented using Green's function, as follows: The scalar wave equation can be expressed in the form of the following formula 5: (5) in, Indicates underground imaging points, This represents the wavefield value at the underground imaging point; It is the angular frequency; This indicates the velocity at the imaging point.

[0023] There are many solutions to the scalar wave equation of the form Equation 5. The solution to this equation is the wave field at the imaging point, which can be approximated by the Green's function. If we use the Rayleigh integral formula, we can obtain the following equation, Equation 6: (6) in, Indicates underground imaging points The location comes from the artillery position. Contributed seismic wave field; The coordinates of the receiving point on the Earth's surface are represented in the formula. These are the components in the three directions of the receiving point; Represents the imaging point Source from receiving point The Green's function representation of the contributed seismic wavefield; This represents the earthquake records received at the Earth's surface. In Equation 6, the partial differential form of the equation can be replaced by the following Equation 7: (7) Similarly, the contribution from the receiving point at the imaging point can also be represented by the Green's function. Nearby wave field As shown in the following formula 8: (8) in, Represents the imaging point The location comes from the artillery position. The Green's function representation of the contributed seismic wavefield; S5.3: Correlation imaging of the uplink and downlink wave fields under the constraint of a Gaussian window is performed as follows: In the equation shown in Formula 4, it is clear that when the receiver point Distance beam center point At greater distances, errors are more likely to occur. To prevent the accumulation of errors, a series of overlapping Gaussian window functions can be introduced into the surface observation system to impose certain restrictions. Since the Gaussian function has the property that the integral sum of its continuous form is one, and the superposition sum of its discrete form is close to one, it can be expressed as the following formula 9: (9) in, Represents angular frequency. It is a reference frequency. This represents the initial beam width of the Gaussian beam. The bundle center spacing can be represented by the following formula 10: (10) As can be seen from the equation represented by Formula 8, the beam center spacing is related to the angular frequency, the reference frequency, and the initial beam width. However, experiments have shown that if the beam center spacing calculated using the above equation is added to the discrete form of the Gaussian function, the relative error can be guaranteed to be within one percent. It can be said that the approximation effect is quite good. Therefore, if the Gaussian window function is introduced into the equation represented by Formula 4, it will not only have no effect on the points near the center of the bundle, but also give the weights of points far from the center of the bundle close to zero, that is, we can ignore its influence. This is the constraint of the Gaussian window form. In imaging, we introduce the cross-correlation imaging condition, namely the complex conjugate product integral of the upgoing and downgoing wave fields, and substitute the wave field representation in the form of the Green's function into it (e.g. Figure 3 As shown), we can obtain the following formula 11: (11) in, This is the final pre-stack offset imaging value. Represents complex conjugation; Substituting the Green's function, expressed as a Gaussian beam, into Equation 11 above, we obtain the most basic Gaussian beam pre-stack migration imaging formula, as shown in Equation 12 below: (12) in, Indicates the coordinates of the imaging point. Indicates the coordinates of the receiver point. Represents the coordinates of the source point. The ray parameter vector representing the imaging point, and This is the final pre-stack imaging value; S5.4: Based on the mesh size of the secondary subdivision, calculate the Green's function at the imaging point obtained in S5.2, which originates from the shot point and the receiver point. and Interpolation is performed so that each grid point in the secondary subdivision has a corresponding Green's function. The up and down wave fields of S5.3 are repeated under the constraint of Gaussian window to perform correlation imaging, and Gaussian beam migration imaging results with variable grid constraint are obtained.

[0024] The working principle of this embodiment is the same as that of Embodiment 1, and will not be repeated here.

[0025] To illustrate the method of the present invention more specifically, based on Examples 1 and 2, the following two sets of imaging targets from Examples 3 and 4 will be used as examples for further explanation.

[0026] Example 3: The computational efficiency advantage of the method of this invention was tested by selecting a simple depression-shaped imaging target, such as... Figure 4 As shown, where Figure 4 (a) shows the velocity field of the depression-shaped imaging target. Figure 4 (b) is the result of fine-grid Gaussian beam migration, with a grid sampling rate of 2, meaning that 4 points make up one grid; Figure 4 (c) is the result of coarse-grid Gaussian beam migration, with a grid sampling rate of 30, that is, 900 points make up one grid; Figure 4 (d) is the result of the variable grid Gaussian beam offset of the present invention.

[0027] By scanning a smooth velocity field, such as Figure 5 As shown in (a), the minimum velocity arrangement can be obtained; then, according to Formula 1, the mesh is divided, as follows: Figure 5 As shown in (b), finally, a variable mesh model field is applied for migration imaging, as follows: Figure 4 As shown in (d).

[0028] Simple comparison Figure 4 The imaging results of the three methods show that the fine-mesh algorithm and the variable-mesh algorithm of this invention produce similar results, while the coarse-mesh algorithm is slightly less accurate, exhibiting localized artifacts and jagged discontinuities in the reflective layer. Furthermore, in terms of computation time and memory usage, ... Figure 6As shown, the fine-mesh algorithm has a much higher memory footprint and computation time than the variable-mesh and coarse-mesh algorithms, making it extremely expensive. While the variable-mesh algorithm is more computationally expensive than the coarse-mesh algorithm, its cost is still within a reasonable range and acceptable. Therefore, considering both computational accuracy and cost, the method of this invention is more competitive and offers high cost-effectiveness.

[0029] The working principle of this embodiment is the same as that of Embodiment 1, and will not be repeated here.

[0030] Example 4: To further demonstrate the advantages of this invention in target imaging, a complex Sigsbee-type imaging target was selected for testing, such as... Figure 7 As shown, where Figure 7 (a) is the velocity field model of the imaging target, with a grid size of 1601×601 and a grid spacing of 10m in both the horizontal and vertical directions. The strata of the imaging target are undulating, with a huge irregular high-speed salt body in the central region and two rows of small-scale diffraction points in the middle and deep layers. Whether diffraction target identification and high-precision imaging under salt can be achieved is the key to this test. Figure 7 (b) is Figure 7 (a) The result after Gaussian smoothing, as the offset velocity field.

[0031] The forward simulation dataset consists of 200 shots, with 801 receivers per shot, a receiver spacing of 20m, a time sampling interval of 2ms, and a record length of 4s. The final migration results are as follows: Figure 8 As shown, where, Figure 8 (a) is the result of conventional Gaussian beam migration. Figure 8 (b) shows the result of the variable grid Gaussian beam migration used in this invention. In this process, we... Figure 8 (b) During the variable-grid Gaussian beam migration process, three regions were subjected to secondary fine-grid refinement. These regions contain obvious features such as small-scale diffractors or steeply dipping strata. The results are as follows: Figure 9 As shown, where, Figure 9 (a), 9(c), and 9(e) are the results of the conventional Gaussian beam method. Figure 9 (b), (d), and (f) are the results of the method of the present invention. A comparison clearly shows that the method of the present invention is superior in both the identification of small-scale diffractors and the imaging of steeply tilted structures.

[0032] The working principle of this embodiment is the same as that of Embodiment 1, and will not be repeated here.

Claims

1. A Gaussian beam target imaging method with variable mesh constraints, characterized in that, Includes the following steps: S1: For the set imaging target, acquire the target's seismic data, the observation data of the observation system, and the migration parameters of the Gaussian beam migration; S2: Based on the data and parameters collected in S1, establish the offset velocity field model and set the average sampling rate. and highest sampling rate The value; S3: Based on the vertical velocity difference of the offset velocity field model, the first global mesh of the imaging domain is performed according to the meshing formula; S4: Perform a second mesh subdivision of the imaging domain based on the characteristics of different target images to refine the mesh of the imaging target layer; S5: Employs up-ray tracing technology for imaging points to construct a reverse extended wavefield, achieving high-precision imaging of the target.

2. The variable mesh-constrained Gaussian beam target imaging method according to claim 1, characterized in that, The offset parameters include the lateral sampling points of the offset velocity field. and longitudinal sampling points Spatial sampling interval Time sampling interval Number of time sampling points , main frequency Reference frequency highest frequency , number of offset shots in earthquake records and the angle of ray emission and The parameters.

3. The Gaussian beam target imaging method with variable mesh constraints according to claim 1, characterized in that, The first global mesh generation of the imaging domain in S3 includes the following steps: S3.1: Perform a global scan analysis on the initial offset velocity field and calculate the average velocity. and minimum speed Simultaneously, based on the velocity change gradient, velocity layers are divided and the minimum velocity arrangement of each layer is calculated. ; S3.2: Based on the calculation results of S3.1 and the average sampling rate set in S2. and highest sampling rate The first global mesh partitioning of the imaging domain is performed using the following formula 1. （1） in, This indicates the grid scale arrangement.

4. The variable mesh-constrained Gaussian beam target imaging method according to claim 3, characterized in that, The method for dividing velocity layers based on velocity change gradients in S3.1 is as follows: Velocity layers are defined by the rate of change of velocity between adjacent sampling points; If the rate of change of velocity at a sampling point changes drastically, that is, the rate of change of velocity at that sampling point is greater than 20% of the average rate of change of velocity at the two adjacent sampling points, then this point is considered to be the interface of a layer.

5. The Gaussian beam target imaging method with variable mesh constraints according to claim 1, characterized in that, The method for performing the second mesh subdivision of the imaging domain in S4 is as follows: According to each layer The secondary grid size is determined to optimize the grid scale of the target layer. The grid size is then adjusted based on the Gaussian beam migration imaging results from the first grid division to ensure fine imaging.

6. The Gaussian beam target imaging method with variable mesh constraints according to claim 5, characterized in that, The method for adjusting the mesh size during the second mesh generation process is as follows: When the speed of a certain layer When the velocity of a layer is simultaneously greater than or simultaneously less than the velocity of its adjacent layers, the imaging domain of that layer is subdivided into grids, adjusted from 5x5 to 2x2 or 1x1.

7. The Gaussian beam target imaging method with variable mesh constraints according to claim 1, characterized in that, S5 includes the following steps: S5.1: Represent the Green's function using Gaussian beam integrals; S5.2: Representing the uplink and downlink wave fields using Green's function; S5.3: Perform correlation imaging of the uplink and downlink wave fields under the constraint of a Gaussian window; S5.4: Based on the grid size of the secondary subdivision, interpolate the Green's function from the shot point and receiver point at the imaging point calculated in S5.2, so that there is a corresponding Green's function at each grid point of the secondary subdivision. Repeat the uplink and downlink wave fields of S5.3 under the constraint of Gaussian window to perform correlation imaging and obtain the Gaussian beam migration imaging result with variable grid constraint.

Images