CN104597488B - Optimum design method of finite difference template of non-equiangular long-grid wave equation - Google Patents

Abstract

The invention discloses an optimum design method of a finite difference template of a non-equiangular long-grid wave equation. The method comprises the following steps: according to a time sampling interval and a spatial sampling interval, mesh subdivision is executed to the simulating region of an actual geologic model; according to the maximum appointed permissible error and the wave number range, the corresponding operator lengths of different acoustic velocities are calculated; based on the finite difference method of the least square optimization time space domain and the corresponding operator length, the optimized finite difference coefficient of the grid is acquired; the acquired optimized finite difference coefficient is brought into the wave equation with a difference scheme to perform forward modeling of the wave equation. According to the method, the finite difference method of the time space domain only applied to square and cube grids is expanded into the rectangular or rectangular parallelepiped grids, thus a special simulation precision requirement and a requirement for saving a calculated amount are met in actual production.

Description

Non-equilateral long-grid wave equation finite difference template optimization design method

Technical Field

The invention relates to the technical field of seismic wave forward numerical simulation, in particular to a finite difference template optimization design method of a non-equilateral long-grid wave equation.

Background

The seismic wave forward numerical simulation technology is a process of propagating waves in a complex geological model (including isotropic media, anisotropic media, Biot multiphase anisotropic media, random hole media and the like) under the condition that the complex geological model is known, and receiving the waves by detectors arranged on the earth surface or underground through multiple transmission, reflection and scattering of underground geological structures by using a numerical calculation method. The accurate wave equation numerical solution is utilized to simulate the seismic response of the underground complex geological structure, and a more scientific mathematical physical basis is provided for the research of seismic wave propagation mechanism, special processing method of seismic data, the explanation of complex stratum and other aspects. In recent years, wave equation numerical simulation methods have been widely used in reverse time migration and full waveform inversion.

There are several methods of wave equation forward modeling, the more common being: finite difference methods, pseudo-spectral methods, finite element methods, boundary element methods, spectral element methods, and the like. The finite difference method is widely used because of small calculation amount and high calculation efficiency, and can be suitable for a more complex speed model. The finite difference method can be classified into: explicit finite differences and implicit finite differences; regular grid finite differences, staggered grid finite differences, and rotated staggered grid finite differences. In the finite difference method, the difference coefficient may be obtained by taylor series expansion or an optimization method, and corresponds to a finite difference based on taylor series expansion and a finite difference based on optimization, respectively. In the conventional finite difference method, the difference coefficient is obtained by minimizing the dispersion relation of the spatial domain. In recent years, a time-space domain finite difference method has appeared, which has higher simulation accuracy and better stability by minimizing the dispersion relation between a time domain and a space domain to obtain a difference coefficient.

The current finite difference method of the time-space domain requires equal space sampling intervals in all directions, namely, the model needs to be divided into square or cubic grids. In actual production, in order to meet specific accuracy requirements or to save computation, we often need to subdivide the model into rectangular or rectangular parallelepiped grids (usually, the grid spacing in the depth direction is different from the horizontal direction). The time-space domain finite difference method suitable for square and cube grids cannot meet specific precision requirements or cannot save calculated amount.

Disclosure of Invention

The embodiment of the invention provides a non-equilateral long grid wave equation finite difference template optimization design method, which expands a time-space domain finite difference method only suitable for square and cube grids into rectangular or cuboid grids to meet the specific simulation precision requirement and the requirement of saving calculated amount in actual production, and comprises the following steps:

performing mesh generation on a simulation area of the actual geological model according to the time sampling interval and the space sampling interval;

calculating corresponding operator lengths for different acoustic wave speeds according to the given maximum allowable error and wave number range;

obtaining an optimized finite difference coefficient of the grid based on a time-space domain finite difference method of least square optimization and the corresponding operator length;

substituting the obtained optimized finite difference coefficient into a wave equation in a difference format to carry out forward simulation of the wave equation;

the differential mesh generation is carried out on the simulation area of the actual geological model according to the time sampling interval and the space sampling interval, and the differential mesh generation comprises the following steps: dividing a simulation area of the actual geological model into a cuboid grid;

based on the least square optimization time-space domain finite difference method and the corresponding operator length, the method obtains the optimization finite difference coefficient of the grid, and comprises the following steps: obtaining an optimized finite difference coefficient of the rectangular grid;

the optimized finite difference coefficient of the cuboid grid is obtained and calculated according to the following formula:

wherein,

where b is the wavenumber, M is the operator length, a_mIn order to optimize the finite difference coefficient, theta is the included angle between the plane wave propagation direction and the horizontal plane, and theta ∈ [0, pi ]]Phi is the azimuth angle of plane wave propagation, phi ∈ [0,2 pi ]]；V is acoustic velocity, tau is time sampling interval, h is sampling interval in x and y directions, c is delta z/h, c is parameter variable, delta z is sampling interval in z direction, β is kh, k is parameter variable, β is wave number range, β∈ [0, b]M is a parameter variable, M is an integer, M ∈ [1, M]N is a parameter variable, n is an integer, n ∈ [1, M ]]。

In one embodiment, the maximum error corresponding to the optimized finite difference coefficient of the rectangular parallelepiped grid satisfies the following constraint condition:

ξ_1max<η；

wherein, ξ_1maxThe maximum error corresponding to the optimized finite difference coefficient of the rectangular grid, and η is the maximum allowable error.

In one embodiment, the maximum error ξ corresponding to the optimized finite difference coefficient of the cuboid mesh_1maxCalculated according to the following formula:

wherein,

in one embodiment, substituting the obtained optimized finite difference coefficients into the wave equation in the difference format to perform wave equation forward simulation includes: substituting the obtained optimized finite difference coefficient of the cuboid grid into a three-dimensional sound wave equation of a difference format to carry out forward simulation of the three-dimensional sound wave equation; the three-dimensional sound wave equation of the difference format is as follows:

wherein P is sound pressure.

The embodiment of the invention provides a finite difference template optimization design method of a non-equilateral long mesh wave equation, which comprises the following steps:

performing mesh generation on a simulation area of the actual geological model according to the time sampling interval and the space sampling interval;

calculating corresponding operator lengths for different acoustic wave speeds according to the given maximum allowable error and wave number range;

obtaining an optimized finite difference coefficient of the grid based on a time-space domain finite difference method of least square optimization and the corresponding operator length;

substituting the obtained optimized finite difference coefficient into a wave equation in a difference format to carry out forward simulation of the wave equation;

the differential mesh generation is carried out on the simulation area of the actual geological model according to the time sampling interval and the space sampling interval, and the differential mesh generation comprises the following steps: dividing a simulation area of the actual geological model into rectangular grids;

based on the least square optimization time-space domain finite difference method and the corresponding operator length, the method obtains the optimization finite difference coefficient of the grid, and comprises the following steps: obtaining an optimized finite difference coefficient of the rectangular grid;

the optimized finite difference coefficient of the rectangular grid is obtained and calculated according to the following formula:

wherein,

where b is the wavenumber, M is the operator length, a_mIn order to optimize the finite difference coefficient, theta is the included angle between the plane wave propagation direction and the horizontal plane, and theta ∈ [0, pi ]]；V is acoustic velocity, tau is time sampling interval, h is x direction sampling interval, c is delta z/h, c is parameter variable, delta z is z direction sampling interval, β is kh, k is parameter variable, β is wave number range, β∈ [0, b]M is a parameter variable, M is an integer, M ∈ [1, M]N is a parameter variable, n is an integer, n ∈ [1, M ]]。

In one embodiment, the maximum error corresponding to the optimized finite difference coefficient of the rectangular grid satisfies the following constraint condition:

ξ_2max<η；

wherein, ξ_2maxMaximum error corresponding to the optimized finite difference coefficient for the rectangular grid, η is the maximum allowable error.

In one embodiment, the maximum error ξ corresponding to the optimized finite difference coefficients of the rectangular grid_2maxCalculated according to the following formula:

wherein,

in one embodiment, the substituting the obtained optimized finite difference coefficient into the wave equation in the difference format to perform the wave equation forward modeling further includes: substituting the obtained optimized finite difference coefficient of the rectangular grid into a two-dimensional sound wave equation in a difference format to carry out forward simulation of the two-dimensional sound wave equation; the two-dimensional acoustic wave equation of the difference format is as follows:

wherein P is sound pressure.

In the embodiment of the invention, the method for optimally designing the non-equilateral long-mesh wave equation finite difference template expands the time-space domain finite difference method only suitable for square and cubic meshes into rectangular or cuboid meshes, and meets the specific simulation precision requirement and the requirement of saving calculated amount in actual production.

Drawings

The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this application, illustrate embodiment(s) of the invention and together with the description serve to explain the principles of the invention. In the drawings:

FIG. 1 is a flow chart of a finite difference template optimization design method of a non-equilateral long mesh wave equation provided in an embodiment of the present invention;

FIG. 2 is a diagram illustrating the variation of an error curve with propagation direction obtained by a conventional finite difference method according to an embodiment of the present invention;

FIG. 3 is a diagram of the variation of an error curve along with the propagation direction obtained by the finite difference template optimization design method of the non-equilateral long mesh wave equation provided by the embodiment of the present invention;

FIG. 4 is a graph of the variation law of an error curve with operator length obtained by a conventional finite difference method according to an embodiment of the present invention;

FIG. 5 is a diagram of the variation law of an error curve with operator length obtained by a finite difference template optimization design method of a non-equilateral long mesh wave equation provided by the embodiment of the present invention;

fig. 6 is a comparison graph of wave field snapshots at 1.0s of a uniform acoustic wave medium, obtained by a conventional finite difference method and a non-equilateral long mesh wave equation finite difference template optimization design method, respectively, according to an embodiment of the present invention;

FIG. 7 is a comparison of waveforms at the location of the dashed line in FIG. 6;

fig. 8 is a comparison diagram of wave field snapshots at 1.0s of a Marmousi model respectively obtained by a conventional finite difference method and a non-equilateral long mesh wave equation finite difference template optimization design method according to an embodiment of the present invention;

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail with reference to the following embodiments and accompanying drawings. The exemplary embodiments and descriptions of the present invention are provided to explain the present invention, but not to limit the present invention.

The inventor finds that the current time-space domain finite difference method requires equal space sampling intervals in all directions, namely, a model needs to be divided into square or cubic grids, but the method cannot meet specific precision requirements or can not save calculation amount. If the model is divided into rectangular or cuboid grids, specific precision requirements can be met or calculated amount can be saved. Based on the method, the invention provides a finite difference template optimization design method of a non-equilateral long mesh wave equation.

Fig. 1 is a flowchart of a finite difference template optimization design method for a non-equilateral long mesh wave equation, as shown in fig. 1, the method includes:

step 101: performing mesh generation on a simulation area of the actual geological model according to the time sampling interval and the space sampling interval;

step 102: calculating corresponding operator lengths for different acoustic wave speeds according to the given maximum allowable error and wave number range;

step 103: obtaining an optimized finite difference coefficient of the grid based on a time-space domain finite difference method of least square optimization and the corresponding operator length;

step 104: and substituting the obtained optimized finite difference coefficient into a wave equation in a difference format to carry out forward simulation of the wave equation.

In specific implementation, the simulation region (or the calculation region) of the actual geological model needs to be divided into rectangular meshes or rectangular grids (usually, the grid spacing in the depth direction is different from the horizontal direction).

The three-dimensional wave equation of sound waves is as follows:

in the formula, P represents a sound pressure, and V represents a sound wave velocity.

Taking the case that the grid distance in the z direction is different from the x and y directions as an example, the specific optimized finite difference format of the rectangular parallelepiped grid is as follows:

wherein τ is a time sampling interval; h is sampling interval in x and y directions; c is delta z/h, c is a parameter variable, and delta z is a sampling interval in the z direction; a is_mFor the optimized finite difference coefficient, M represents the operator length, M is a parameter variable, M is an integer, M ∈ [1, M]。

To simplify notation, we willIs abbreviated asAccording to plane wave theory, we order:

wherein,

wherein i, l, j and n are parameter variables, n is an integer, and n belongs to [1, M ]; k is a parameter variable; theta belongs to [0, pi ] is the included angle between the plane wave propagation direction and the horizontal plane, and phi belongs to [0,2 pi ] is the azimuth angle of the plane wave propagation. Substituting equations (3) and (4) into equation (2), and appropriately simplifying, can obtain:

wherein,then according to the constraint condition;

the following can be obtained:

wherein, beta is kh, beta is wave number range, beta belongs to [0, b ]. The finite difference algorithm is optimized by solving a set of difference coefficients to minimize the error at both ends of equation (7). According to the least squares theory, this problem can be transformed to solve the following system of equations:

wherein,

by solving the system of linear equations (8), we can obtain the optimized difference coefficients of the cuboid grids. And substituting the difference coefficient into the equation (2) to solve the three-dimensional sound wave equation.

In order to ensure the effectiveness of the method of the present invention, the following constraint conditions need to be added in the previous optimization process, that is, the maximum error corresponding to the optimized finite difference coefficient of the rectangular parallelepiped grid satisfies the following constraint conditions:

ξ_1max<η (10)

we calculate the maximum error for the optimized difference coefficient by:

wherein,

where η is the maximum allowable error.

ξ_1maxRelating only to b and M when M is fixed, ξ_1maxDetermined by b, ξ_1maxTherefore, if the initial b does not satisfy equation (10), we should go through progressively decreasing b up to ξ_1max<η to obtain optimum b, additionally, when b is fixed, ξ_1maxRelated to M only, ξ_1maxTherefore, if the initial M does not satisfy equation (10), we should increase M gradually until ξ_1max<η to obtain the best M.

For a two-dimensional rectangular grid, taking the case that the grid distance in the z direction is different from that in the x direction as an example, a specific optimized finite difference format is as follows:

wherein P represents sound pressure and V represents sound wave velocity; wherein τ is a time sampling interval; h is a sampling interval in the x direction; c is delta z/h, c is a parameter variable, and delta z is a sampling interval in the z direction; a is_mFor the optimized finite difference coefficient, M represents the operator length, M is a parameter variable, M is an integer, M ∈ [1, M]。

We can obtain the optimized difference coefficients for the rectangular grid by solving the following linear equations:

wherein,

by solving the system of linear equations (15), we can obtain the optimized difference coefficients for the rectangular grid. And substituting the difference coefficient into equation (13) to solve the two-dimensional sound wave equation.

Similarly, in order to ensure the effectiveness of the method of the present invention, the following constraint condition needs to be added in the foregoing optimization process, that is, the maximum error corresponding to the optimized finite difference coefficient of the rectangular grid satisfies the following constraint condition:

ξ_2max<η (16)

we calculate the maximum error for the optimized difference coefficient by:

wherein,

where η is the maximum allowable error.

The advantages of the present invention are illustrated below by way of example with a uniform acoustic medium. The longitudinal wave velocity is 1500M/s, the time sampling interval is 1ms, the sampling interval in the x direction is 10M, the sampling interval in the z direction is 5M, and the operator length M belongs to [2, 10 ].

The method (namely the finite difference template optimization design method of the non-equilateral long mesh wave equation) and the traditional finite difference method are adopted to obtain the change rule of the error curve along with the propagation direction, which is shown in the figures 2 and 3 respectively, and the rule η is 10^-2.5Approximately equals to 0.003, the method of the invention and the traditional finite difference method are adopted to obtain the change rule of the error curve along with the operator length, as shown in fig. 4 and 5 respectively, as can be seen from fig. 2 to 5, compared with the taylor-based time-space domain finite difference method (traditional finite difference method), the method of the invention has wider effective frequency band and smaller numerical dispersion, and the regulation η is 10^-2.5And is approximately equal to 0.003, a wave field snapshot comparison diagram of a uniform acoustic wave medium is obtained by adopting a traditional finite difference method and the method disclosed by the invention, and is shown in figure 6. It is composed ofIn fig. 6, the left graph is a wave field snapshot of the uniform acoustic wave medium obtained by the conventional finite difference method, and the right graph is a wave field snapshot of the uniform acoustic wave medium obtained by the method of the present invention. Fig. 7 is a comparison graph of waveforms at the positions of the dotted lines in fig. 6. As can be seen from FIGS. 6 and 7, the simulation accuracy of the method of the present invention is higher when the operator lengths are the same.

To further illustrate the effectiveness of the present invention, we performed forward modeling with a Marmousi model. The time sampling interval is 1ms, the x-direction sampling interval is 12m, and the z-direction sampling interval is 9 m. The operator length M of the conventional finite difference method is 9. The finite difference of the time-space domain adopts an algorithm of changing the length of an operator, operators with different lengths are selected according to different speeds, and M belongs to [2, 9 ]. Fig. 8 is a comparison diagram of wave field snapshots of the Marmousi model at the time of 1.0s, which are respectively obtained by the conventional finite difference method and the method of the present invention, provided in the embodiment of the present invention, the upper graph is a wave field snapshot of the Marmousi model at the time of 1.0s, which is obtained by the conventional finite difference method, and the lower graph is a wave field snapshot of the Marmousi model at the time of 1.0s, which is obtained by the method of the present invention; as can be seen from FIG. 8, the method of the present invention is suitable for forward simulation of complex media, and the simulation effect is significantly better than that of the conventional finite difference method.

In summary, the method of the present invention has the following advantages: 1. the numerical dispersion of the middle and high frequency bands is reduced. 2. The method is more practical based on the time-space domain dispersion relation, and has higher simulation precision. 3. The method is suitable for rectangular and cuboid grids, and can meet specific precision requirements in actual production and save the calculation amount.

The above description is only a preferred embodiment of the present invention, and is not intended to limit the present invention, and various modifications and changes may be made to the embodiment of the present invention by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (8)

1. A finite difference template optimization design method of a non-equilateral long mesh wave equation is characterized by comprising the following steps:

performing mesh generation on a simulation area of the actual geological model according to the time sampling interval and the space sampling interval;

calculating corresponding operator lengths for different acoustic wave speeds according to the given maximum allowable error and wave number range;

obtaining an optimized finite difference coefficient of the grid based on a time-space domain finite difference method of least square optimization and the corresponding operator length;

substituting the obtained optimized finite difference coefficient into a wave equation in a difference format to carry out forward simulation of the wave equation;

the differential mesh generation is carried out on the simulation area of the actual geological model according to the time sampling interval and the space sampling interval, and the differential mesh generation comprises the following steps: dividing a simulation area of the actual geological model into a cuboid grid;

based on the least square optimization time-space domain finite difference method and the corresponding operator length, the method obtains the optimization finite difference coefficient of the grid, and comprises the following steps: obtaining an optimized finite difference coefficient of the rectangular grid;

the optimized finite difference coefficient of the cuboid grid is obtained and calculated according to the following formula:

wherein,

where b is the wavenumber, M is the operator length, a_mIn order to optimize the finite difference coefficient, theta is the included angle between the plane wave propagation direction and the horizontal plane, and theta ∈ [0, pi ]]Phi is the azimuth angle of plane wave propagation, phi ∈ [0,2 pi ]]；V is acoustic velocity, tau is time sampling interval, h is sampling interval in x and y directions, c is delta z/h, c is parameter variable, delta z is sampling interval in z direction, β is kh, k is parameter variable, β is wave number range, β∈ [0, b]M is a parameter variable, M is an integer, M ∈ [1, M]N is a parameter variable, n is an integer, n ∈ [1, M ]]。

2. The method of claim 1, wherein the maximum error corresponding to the optimized finite difference coefficients of the rectangular parallelepiped lattice satisfies the following constraint:

ξ_1max<η；

wherein, ξ_1maxThe maximum error corresponding to the optimized finite difference coefficient of the rectangular grid, and η is the maximum allowable error.

3. The method of claim 2, wherein the optimized finite difference coefficients of the cuboid mesh correspond to a maximum error of ξ_1maxCalculated according to the following formula:

${\mathrm{\&xi;}}_{1max}=\underset{\mathrm{\&beta;}\&Element;\&lsqb;0,b\&rsqb;,\mathrm{\&theta;}\&Element;\&lsqb;0,2\mathrm{\&pi;}\&rsqb;}{max}\left|{\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

wherein,

${\mathrm{\&xi;}}_1(\mathrm{\&beta;},\mathrm{\&theta;})=\frac{2}{r\mathrm{\&beta;}}\arcsin \sqrt{r^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\left[\begin{array}{c}si{n}^2\left(\frac{m\mathrm{\&beta;}}{2}cos\mathrm{\&theta;}cos\mathrm{\&phi;}\right)-si{n}^2\left(\frac{m\mathrm{\&beta;}}{2}cos\mathrm{\&theta;}sin\mathrm{\&phi;}\right)\\ -\frac{1}{c^2}{\sin }^2\left(\frac{mc\mathrm{\&beta;}}{2}sin\mathrm{\&theta;}\right)\end{array}\right]}-1.$

4. the method of claim 3, wherein substituting the obtained optimized finite difference coefficients into the wave equation in the difference format for wave equation forward simulation comprises: substituting the obtained optimized finite difference coefficient of the cuboid grid into a three-dimensional sound wave equation of a difference format to carry out forward simulation of the three-dimensional sound wave equation; the three-dimensional sound wave equation of the difference format is as follows:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,y,z}^{t-1}-2P_{x,y,z}^t+P_{x,y,z}^{t+1})=\frac{1}{{\left(ch\right)}^2}\&lsqb;a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,y,z-m}^t+P_{x,y,z+m}^t)\&rsqb;\\ +\frac{1}{h^2}\&lsqb;2a_0{P}_{x,y,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,y,z}^t+P_{x+m,y,z}^t+P_{x,y-m,z}^t+P_{x,y+m,z}^t)\&rsqb;\end{array};$

wherein P is sound pressure.

5. A finite difference template optimization design method of a non-equilateral long mesh wave equation is characterized by comprising the following steps:

performing mesh generation on a simulation area of the actual geological model according to the time sampling interval and the space sampling interval;

calculating corresponding operator lengths for different acoustic wave speeds according to the given maximum allowable error and wave number range;

obtaining an optimized finite difference coefficient of the grid based on a time-space domain finite difference method of least square optimization and the corresponding operator length;

substituting the obtained optimized finite difference coefficient into a wave equation in a difference format to carry out forward simulation of the wave equation;

the differential mesh generation is carried out on the simulation area of the actual geological model according to the time sampling interval and the space sampling interval, and the differential mesh generation comprises the following steps: dividing a simulation area of the actual geological model into rectangular grids;

based on the least square optimization time-space domain finite difference method and the corresponding operator length, the method obtains the optimization finite difference coefficient of the grid, and comprises the following steps: obtaining an optimized finite difference coefficient of the rectangular grid;

the optimized finite difference coefficient of the rectangular grid is obtained and calculated according to the following formula:

wherein,

where b is the wavenumber, M is the operator length, a_mIn order to optimize the finite difference coefficient, theta is the included angle between the plane wave propagation direction and the horizontal plane, and theta ∈ [0, pi ]]；V is acoustic velocity, tau is time sampling interval, h is x direction sampling interval, c is delta z/h, c is parameter variable, delta z is z direction sampling interval, β is kh, k is parameter variable, β is wave number range, β∈ [0, b]M is a parameter variable, M is an integer, M ∈ [1, M]N is a parameter variable, n is an integer, n ∈ [1, M ]]。

6. The method of claim 5, wherein the maximum error corresponding to the optimized finite difference coefficients of the rectangular grid satisfies the following constraint:

ξ_2max<η；

wherein, ξ_2maxOptimization for rectangular gridThe maximum error corresponding to the finite difference coefficient, η is the maximum allowable error.

7. The method of claim 6, wherein the optimized finite difference coefficients of the rectangular grid correspond to a maximum error of ξ_2maxCalculated according to the following formula:

${\mathrm{\&xi;}}_{2max}=\underset{\mathrm{\&beta;}\&Element;\&lsqb;0,b\&rsqb;,\mathrm{\&theta;}\&Element;\&lsqb;0,2\mathrm{\&pi;}\&rsqb;}{max}\left|{\mathrm{\&xi;}}_2(\mathrm{\&beta;},\mathrm{\&theta;})\right|;$

wherein,

8. the method of claim 7, wherein the fitting the obtained optimized finite difference coefficients into the wave equation in the difference format for wave equation forward simulation further comprises: substituting the obtained optimized finite difference coefficient of the rectangular grid into a two-dimensional sound wave equation in a difference format to carry out forward simulation of the two-dimensional sound wave equation; the two-dimensional acoustic wave equation of the difference format is as follows:

$\begin{array}{c}\frac{1}{V^2{\mathrm{\&tau;}}^2}(P_{x,z}^{t-1}-2P_{x,z}^t+P_{x,z}^{t+1})=\frac{1}{{\left(ch\right)}^2}\&lsqb;a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x,z-m}^t+P_{x,z+m}^t)\&rsqb;\\ +\frac{1}{h^2}\&lsqb;a_0{P}_{x,z}^t+\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{x-m,z}^t+P_{x+m,z}^t)\&rsqb;\end{array};$

wherein P is sound pressure.