CN105277980A - High-precision spatial and temporal arbitrary multiple variable grid finite difference forward modeling method - Google Patents

Abstract

The invention provides a high-precision spatial and temporal arbitrary multiple variable grid finite difference forward modeling method, which comprises the steps of: step 1, establishing a forward modeling speed model of an underground medium; step 2, subjecting the forward modeling speed model to two-dimensional grid discretization, and subjecting an acoustic wave field in the forward modeling speed model to two-dimensional grid discretization, wherein the acoustic wave field is located on grid nodes; step 3, subjecting a perfectly matched layer boundary condition to grid discretization; and step 4, conducting time-domain finite difference forward modeling simulation through an acoustic wave equation, wherein time sampling step lengths are variable step lengths at grids of different sizes. The high-precision spatial and temporal arbitrary multiple variable grid finite difference forward modeling method can effectively suppress forward modeling frequency dispersion, significantly increases signal-to-noise ratio, and greatly improves precision and efficiency of forward modeling simulation.

Description

High-precision spatial and time any multiple variable grid finite difference forward modeling method

Technical field

The present invention relates to seismic forward modeling analogue technique, particularly relate to a kind of high-precision spatial and time any multiple variable grid finite difference forward modeling method.

Background technology

Seismic forward modeling analogue technique is an important research field of domestic and international geophysics educational circles, and based on the computing method of wave theory, there is application distant view most, because it can the kinematics and dynamics feature of retentively seismic wave, the object of accurate analog PROPAGATION CHARACTERISTICS OF SEISMIC can be reached.Wherein, method of finite difference is one of the most frequently used seismic forward modeling method for numerical simulation.It has concept and lists computation scheme simple, to the advantage that uniform dielectric, nonhomogeneous media, anisotropic medium are all suitable for, therefore in widespread attention in earthquake educational circles.

The obvious shortcoming of method of finite difference is difficult to overcome frequency dispersion effect, solve frequency dispersion effect, necessary computations grid, but this can cause the increase of calculated amount, efficiency declines, and due to the solution on net point can only be provided, simulation precision is subject to the restriction of sizing grid etc., in addition when earth's surface big rise and fall or geologic structure complexity, because traditional finite-difference algorithm have employed fixing mesh spacing, the method validity is in some instances restricted.

(" Journal of Engineering Geophysics " Vol.4 such as Li Shengjun, No.3,2007,207-212 page) realize carrying out Moving grids in one direction, now Moving grids does not need to carry out interpolation calculation, ensure that arithmetic accuracy to a certain extent, but method still there is significant limitation, Moving grids calculating can only be carried out in whole model scope.In order to improve dirigibility and the simulation precision of method, Zhao Haibo (" Science Bulletin " Vol.52, No.12,2007,1387-1395 page) have employed Taylor and optimizes Moving grids optimized algorithm, improves operator precision; Grid change direction is extended to two dimension by Zhu Shengwang etc. (" geophysical prospecting for oil " Vol.42, No.6,2007,634-639 page), and have employed Spatial higher order difference, improves dirigibility and the simulation precision of space Moving grids.

But the thinking of above-mentioned variable grid only carries out variable process on space lattice, the efficiency improving simulation is had some limitations.Because require to adopt less time step with the stability of satisfied calculating at employing space, the region small grid discrete way that media variations is violent.And if use traditional step-length unified time, the problem in medium graded region with regard to over-sampling in life period, significantly can reduce simulation precision.Time Problem On Step-size Change remains unsolved.We have invented a kind of new high-precision spatial and time any multiple variable grid finite difference forward modeling method for this reason, solve above technical matters.

Summary of the invention

The object of this invention is to provide one room and time variable grid is combined, realize high-precision spatial and time any multiple variable grid finite difference forward modeling method of sound wave.

Object of the present invention realizes by following technical measures: high-precision spatial and time any multiple variable grid finite difference forward modeling method, this high-precision spatial and time any multiple variable grid finite difference forward modeling method comprise: step 1, that sets up underground medium just drills rate pattern, macrolattice and small grid is divided according to model area, and according to sizing grid, the forward simulation time-sampling step-length of different spaces net region is set; Step 2, align and drill rate pattern and carry out two-dimensional grid discretize, align the acoustic wavefield of drilling in rate pattern and carry out two-dimensional grid discretize, this acoustic wavefield is positioned on grid node; Step 3, to perfectly matched layer boundary condition grid discretization; And step 4, carry out Fdtd Method Forward simulation by Acoustic Wave-equation, wherein time-sampling step-length is variable step at the large small grid of difference.

Object of the present invention also realizes by following technical measures:

In step 2, acoustic wavefield two-dimensional grid discretize formula is as follows:

$\begin{array}{c}{u}_{(i,k)}^{n+1}={2u}_{i,k}^n-u_{i,k}^{n-1}+\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δt}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\\ +\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δz}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\end{array},$

Wherein, u is acoustic wavefield; I and k is two-dimensional grid sequence number; V is speed; Δ t is time-sampling step-length; Δ x is x direction sizing grid; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice.

3. high-precision spatial according to claim 1 and time any multiple variable grid finite difference forward modeling method, is characterized in that, in step 2, the acoustic wavefield two-dimensional grid discretize for the excessive grid in space is as follows:

To zone of transition each point in macrolattice, difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array},$

In small grid, zone of transition each point difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array},$

Wherein, u is acoustic wavefield; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice; N is space lattice multiple.

In step 3, perfectly matched layer boundary condition grid discretization formula is as follows:

$\begin{array}{c}{u}^{n+1}=u_1^{n+1}+u_2^{n+1}\\ =\frac{1}{1+d\left(x\right)\mathrm{\Δt}}[2-d^2\left(x\right){\mathrm{\Δt}}^2)u_1^n+(d\left(x\right)\mathrm{\Δt}-1)u_1^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δx}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\right)]\\ =\frac{1}{1+d\left(z\right)\mathrm{\Δt}}[(2-d^2\left(z\right){\mathrm{\Δt}}^2)u_2^n+(d\left(z\right)\mathrm{\Δt}-1)u_2^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δz}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\right)]\end{array}$

Wherein, u is acoustic wavefield; I and k is two-dimensional grid sequence number; Δ t is time-sampling step-length; Δ x is x direction sizing grid; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice, and d (x) and d (z) are respectively the perfectly matched layer border attenuation term relevant with x and z directional derivative.

In step 4, when changing time-sampling step-length, a transitional zone is set in large time-sampling step-length district, utilize large time-sampling step-length to obtain the wave field value of subsequent time in transitional zone region, employing interpolation method obtains each moment wave field value needed for little time sampling step length calculating.

In step 4, time-sampling step-length is any multiple, and time-sampling step-length interpolation method carries out high-order nonlinear interpolation to multiple spot time-sampling step-length for using Lagrangian simultaneously.

High-precision spatial in the present invention and time any multiple variable grid finite difference forward modeling method, it is sound wave finite difference forward modeling method, can spatially any multiple Moving grids, any multiple variable step time on, fully adapts to the forward simulation calculating of all kinds of complicated underground condition.Such as, weathering layer model, fracture medium model, organic reef model and solution cavity model etc.Small grid is used to complex structure region in forward model, other regions use macrolattice, thus calculate territory, model division, the encryption time slice circulation when calculating small grid region simultaneously, namely adopt large time step in macrolattice region, small grid region adopts little time step.Contrast with conventional finite difference method, counting yield and computational accuracy can be significantly improved.The present invention and other space Moving grids finite difference method contrast, and can just drill variable step on the time in refined net region, and the forward simulation fully adapting to all kinds of complicated underground condition calculates.What the present invention effectively can suppress that weathering layer causes just drills frequency dispersion, improves the signal to noise ratio (S/N ratio) of the entire profile.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of high-precision spatial and time any multiple variable grid finite difference forward modeling method;

Fig. 2 is multistage variable grid model subdivision schematic diagram;

Fig. 3 is zone of transition difference scheme schematic diagram in macrolattice;

Fig. 4 is zone of transition difference scheme schematic diagram in small grid;

Fig. 5 is that completely permutation absorbing boundary loads schematic diagram;

Fig. 6 is variable time step zone of transition processing mode schematic diagram;

Fig. 7 is uniform dielectric model wave field snapshot;

Fig. 8 is fracture medium model schematic;

Fig. 9 is fracture medium model single shot record;

Figure 10 is organic reef model schematic;

Figure 11 is organic reef model single shot record;

Figure 12 is as the criterion northern mountain front regional model schematic diagram;

Figure 13 is as the criterion band model Moving grids region before the Beishan Mountain;

Figure 14 is as the criterion band model regular grid single shot record before the Beishan Mountain;

Figure 15 is as the criterion band model Moving grids single shot record before the Beishan Mountain.

Embodiment

For making above and other object of the present invention, feature and advantage can become apparent, cited below particularly go out preferred embodiment, and coordinate institute's accompanying drawings, be described in detail below.

Fig. 1 is the process flow diagram of high-precision spatial of the present invention and time any multiple variable grid finite difference forward modeling method.Its specific implementation flow process is as follows:

In step 101, that sets up underground medium just drills rate pattern, divides macrolattice and small grid according to model area.And according to sizing grid, the forward simulation time-sampling step-length of different spaces net region is set.

As shown in Figure 2, the present invention realizes multistage variable grid model subdivision.The inside and outside two-layer small grid in the region needing in model to utilize small grid to calculate is carried out nested subdivision, reaches the final simulation precision needed by twice Moving grids.The multiple of inside and outside two-layer Moving grids can be all arbitrary constant.Nested multistage variable mesh generation mode can realize the mesh spacing change of high power by the Moving grids algorithm of low multiple.Require high for computational accuracy, the situation that complex geologic body yardstick is little, this method is more flexible, and simultaneously relative to conventional Moving grids algorithm, this method also corresponding saving can calculate internal memory.Flow process enters into step 102.

In step 102, carry out two-dimensional grid discretize to forward model, carry out two-dimensional grid discretize to the acoustic wavefield in forward model, described acoustic wavefield is positioned on grid node.Wherein sound wave difference discrete formula is:

$\begin{array}{c}{u}_{(i,k)}^{n+1}={2u}_{i,k}^n-u_{i,k}^{n-1}+\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δt}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\\ +\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δz}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\end{array},$

Wherein, u is acoustic wavefield; I and k is two-dimensional grid sequence number; V is speed; Δ t is time-sampling step-length; Δ x is x direction sizing grid; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice.

Difference scheme for the excessive grid in space needs to revise, and to zone of transition each point in macrolattice, difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array},$

Wherein, N is space lattice multiple.

Fig. 3 is zone of transition difference scheme schematic diagram in macrolattice.Zone of transition in macrolattice region 5 points, its Difference Calculation get a little in small grid region and non-sequential get a little, but maintain the form about central point, still use macroreticular difference scheme and difference coefficient.

In small grid, zone of transition each point difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array}.$

Wherein, N is space lattice multiple.

Fig. 4 is zone of transition difference scheme schematic diagram in small grid.Zone of transition in small grid region 5 points, when Difference Calculation needs to use macrolattice point, use point symmetrical with it to participate in calculating in small grid.Now step value there occurs change, needs to revise difference scheme.Flow process enters into step 103.

In step 103, to perfectly matched layer boundary condition grid discretization, perfectly matched layer boundary condition grid discretization formula is as follows:

$\begin{array}{c}{u}^{n+1}=u_1^{n+1}+u_2^{n+1}\\ =\frac{1}{1+d\left(x\right)\mathrm{\Δt}}[2-d^2\left(x\right){\mathrm{\Δt}}^2)u_1^n+(d\left(x\right)\mathrm{\Δt}-1)u_1^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δx}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\right)]\\ =\frac{1}{1+d\left(z\right)\mathrm{\Δt}}[(2-d^2\left(z\right){\mathrm{\Δt}}^2)u_2^n+(d\left(z\right)\mathrm{\Δt}-1)u_2^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δz}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\right)]\end{array}$

Wherein, d (x) and d (z) are respectively the perfectly matched layer border attenuation term relevant with x and z directional derivative.

The structure of perfectly matched layer as shown in Figure 5.Perfectly matched layer border is positioned at outside model zoning, and this border is divided into several zones of different, and the effect difference absorbed wave field in each region, therefore decay factor value is different.In perfectly matched layer absorption region representated by B1 and B2, only decay to the ripple propagated along the z-axis direction in this region; In PML absorption region representated by B3 and B4, decay to the ripple propagated along the x-axis direction in this region; In the perfectly matched layer absorption region representated by C1, C2, C3 and C4, angle point is all decayed to the ripple propagated with z-axis direction along the x-axis direction.Flow process enters into step 104.

In step 104, carry out Fdtd Method Forward simulation by Acoustic Wave-equation, wherein time-sampling step-length is variable step at the large small grid of difference.

When spatial sampling interval changes, in order to ensure the stability of time domain explicit difference method, necessary regulation time sampling interval, makes minimum space sampling step length meet stability condition.If whole model all adopts unified time sampling interval, the time over-sampling in large space step-length region will inevitably be caused.In order to improve the simulation precision of algorithm, need the encryption time slice circulation when calculating small grid region, namely adopt large time step in macrolattice region, small grid region adopts little time step.

Change time step need in large time step district, transitional zone is set, utilize large time step to obtain the wave field value of subsequent time in transitional zone region, and adopt interpolating method obtain little time step calculate needed for each moment wave field value.Fig. 6 is variable time step transitional zone process schematic diagram, and in figure, the point of cross star mark is for needing the point of interpolation.Because time difference operator precision is less than space difference operator precision, comparatively space is low for the precision susceptibility on time dimension, and the error that therefore interpolation is introduced can not cause significantly artificial reflection.

Time step can be any multiple.Time step interpolation method uses Lagrangian to carry out high-order nonlinear interpolation to multiple spot time step simultaneously.Flow process terminates.

Application embodiments of the invention 1 are uniform dielectric model, use the availability of uniform dielectric model measurement variable grid algorithm, namely verify that variable grid algorithm can not produce reflection at Moving grids place.

Homogeneous model size is 1km × 1km, and velocity of longitudinal wave is 5000m/s, and macrolattice spatial mesh size is 10m, and time step is 1ms, tests 40 times of spatially-variable trellis algorithm.Small grid room and time step-length is respectively 0.25m, 0.025ms.Respectively Moving grids district is placed in model central authorities and edge of model.Focus is positioned at model center, and source wavelet all uses the Ricker wavelet of dominant frequency 50Hz.

Fig. 7 is uniform dielectric model snapshot, and net region of attempting to change, a left side is placed in model central authorities, and net region of attempting to change, the right side is placed in edge of model.Test result display variable grid algorithm can not produce artificial interference at Moving grids intersection, and namely this algorithm can ensure simulation precision.

Application embodiments of the invention 2 are fracture medium model, and Fig. 8 is fractured model schematic diagram, and wherein scheming (a) is fractured model schematic diagram, and figure (b) is macrolattice sampled result, and figure (c) is small grid sampled result.The size of model is 100m × 100m, be that 50m place has a bed interface in the degree of depth, upper interval velocity is 4500m/s, lower interval velocity is 6000m/s, and interface central authorities have a height to be the small-sized projection of 3m, and at protruding interface, place exists many cracks, fracture width is 5mm, density is 9/meter, and position angle is 90 °, and in crack, the velocity of longitudinal wave of filled media is 2000m/s.To this model usage space high power variable grid algorithm, space macrolattice sampling step length is 1m, 0.08ms, and small grid step-length is 2.5mm, 0.002ms.Use the conventional finite difference algorithm that sampling step length is 1m, 0.08ms to calculate same model, the Comparative result of the result obtained and Moving grids algorithm simultaneously.The source wavelet of two kinds of method uses is Ricker wavelet, and wavelet dominant frequency is 40Hz.

Two kinds of method fracture models just drill result as shown in Figure 9, wherein left figure is the single shot record that variable grid algorithm obtains, and right figure is the single shot record that regular grid algorithm obtains.The sampling point in crack is not reflected in model after regular grid is discrete, therefore the big gun record obtained in regular grid analogy method only has the reflection of the reflection wave of horizontal interface and small-sized Raised key axis, the obvious epirelief in reflection wave top of Raised key axis, the reflected energy of its lower interface is slightly weakened by Raised key axis reflected wave effects.Moving grids analogy method employs the sampling interval of millimeter magnitude in crack district, ensure that the model of each crack after discrete all has sampling point.

Fractured model confirms that Moving grids algorithm can reach millimetre-sized simulation precision.Run in needed for this program and save as about 1G, if accurate analog crack and all use small grid to sample to whole model, then save as about 2T in needing, current computing machine is difficult to meet the requirement calculated, also considerable during the machine expended.Moving grids algorithm improves simulation precision to a great extent, when having saved internal memory and machine.

Application embodiments of the invention 3 are organic reef model, and Figure 10 is organic reef model schematic, and wherein upper figure is macrolattice sampled result, and figure below is small grid sampled result.Model size is 5km × 2km.Organic reef is passed by tomography, edge is unsmooth, and organic reef inside is dispersed with the solution cavity that multiple radius is less than 5m, use conventional finite difference algorithm and 40 times of variable grid algorithms to calculate this model simultaneously, wherein conventional finite difference method spatial sampling interval is 10m, time sampling interval is 0.5ms, Moving grids algorithm space macrolattice sampling interval is 10m, small grid sampling interval is 0.25m, time sampling interval is respectively 0.5ms and 0.0125ms, and Moving grids scope is mainly for organic reef and small-sized solution cavity wherein.Lateral extent is 3.2km ~ 4.2km, and depth range is 500m ~ 900m.The source wavelet of two kinds of method uses is the Ricker wavelet that dominant frequency is 40Hz.

Figure 11 shows the analog result that two kinds of distinct methods obtain, and left figure is 40 times of Moving grids single shot records, and right figure is conventional macrolattice single shot record.Bottom organic reef and edge has tomography to distribute, and cause bioherm reflection wave along with very multibreak diffracted wave, conventional finite difference method and Moving grids method all can obtain the response of each bed interface and the complicated diffracted wave of organic reef breakpoint.Organic reef inside has more small-sized solution cavity to grow, and conventional finite difference method uses macrolattice to sample to model, does not reflect the sampling point of small-sized solution cavity, the big gun record therefore obtained at conventional finite difference method does not have the response of the inner solution cavity of organic reef.Moving grids method uses small grid to carry out meticulous sampling to organic reef inside, feature the feature of solution cavity comparatively accurately, the big gun record of Moving grids method can be seen, under organic reef top reflective ripple, there is a large amount of mixed and disorderly diffracted waves, its energy is less than the reflection of aspect.

Application embodiments of the invention 4 are as the criterion northern mountain front regional model, and as shown in figure 12, model size is 2.1km × 1.2km to accurate northern mountain front regional model schematic diagram, and model surface layer bit rate is relatively low, and velocity of longitudinal wave is 2472m/s.And model deep layer velocity of longitudinal wave can reach 6500m/s, velocity contrast is large, and complex geologic conditions is changeable, can cause certain influence to Forward modelling result.Moving grids weathering layer region being carried out to five times calculates, and at horizontal 12000m, within the scope of degree of depth 1000m, uses the fine grid blocks of 5m to carry out analog computation to model.Model resampling result as shown in figure 13.Model trace spacing is 25m, and Temporal sampling is 1.5ms, and focus is combination p-wave source, and source wavelet is the Ricker wavelet of dominant frequency 20Hz, and record length is 6s.Shot point is placed on the left of model.

Figure 14 and Figure 15 is band model regular grid single shot record and Moving grids single shot record before the accurate Beishan Mountain respectively.Can be found out by analog result, regular grid is just drilled, due to the existence of weathering layer, make frequency dispersion in analog result quite serious, serious interference is caused to analog result, simultaneously due to the existence of the numerical solidification of strong energy, model lower floor reflected energy is weakened relatively, is unfavorable for the further Treatment Analysis of geological data.And limited Moving grids calculates the frequency dispersion effectively can suppressed low velocity layer (LVL) and cause, thus improve the signal to noise ratio (S/N ratio) of the entire profile.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, within the spirit and principles in the present invention all, any amendment done, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (6)

1. high-precision spatial and time any multiple variable grid finite difference forward modeling method, is characterized in that, this high-precision spatial and time any multiple variable grid finite difference forward modeling method comprise:

Step 1, that sets up underground medium just drills rate pattern, divides macrolattice and small grid, and according to sizing grid, arrange the forward simulation time-sampling step-length of different spaces net region according to model area;

Step 2, align and drill rate pattern and carry out two-dimensional grid discretize, align the acoustic wavefield of drilling in rate pattern and carry out two-dimensional grid discretize, this acoustic wavefield is positioned on grid node;

Step 3, to perfectly matched layer boundary condition grid discretization; And

Step 4, carries out Fdtd Method Forward simulation by Acoustic Wave-equation, and wherein time-sampling step-length is variable step at the large small grid of difference.

2. high-precision spatial according to claim 1 and time any multiple variable grid finite difference forward modeling method, it is characterized in that, in step 2, acoustic wavefield two-dimensional grid discretize formula is as follows:

$\begin{array}{c}{u}_{(i,k)}^{n+1}={2u}_{i,k}^n-u_{i,k}^{n-1}+\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δt}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\\ +\frac{v^2{\left(\mathrm{\Δt}\right)}^2}{{\mathrm{\Δz}}^2}\underset{m=1}{\overset{\frac{M}{2}}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\end{array},$

Wherein, u is acoustic wavefield; I and k is two-dimensional grid sequence number; V is speed; Δ t is time-sampling step-length; Δ x is x direction sizing grid; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice.

3. high-precision spatial according to claim 1 and time any multiple variable grid finite difference forward modeling method, is characterized in that, in step 2, the acoustic wavefield two-dimensional grid discretize for the excessive grid in space is as follows:

To zone of transition each point in macrolattice, difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array},$

In small grid, zone of transition each point difference scheme is:

$\begin{array}{c}\frac{{\∂}^{2}u}{{\∂z}^2}=\frac{1}{{\left(\mathrm{\Δz}\right)}^2}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{m\Δz})+u(x,z-\mathrm{m\Δz})-2u(x,z)]\\ +\frac{1}{{\left(\mathrm{N\Δz}\right)}^2}\underset{m=1+n}{\overset{M/2}{\mathrm{\Σ}}}{c}_m^{\′}[u(x,z+\mathrm{mN\Δz})+u(x,z-\mathrm{mN\Δz})-2u(x,z)]\end{array},$

Wherein, u is acoustic wavefield; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice; N is space lattice multiple.

4. high-precision spatial according to claim 1 and time any multiple variable grid finite difference forward modeling method, it is characterized in that, in step 3, perfectly matched layer boundary condition grid discretization formula is as follows:

$\begin{array}{c}{u}^{n+1}=u_1^{n+1}+u_2^{n+1}\\ =\frac{1}{1+d\left(x\right)\mathrm{\Δt}}[2-d^2\left(x\right){\mathrm{\Δt}}^2)u_1^n+(d\left(x\right)\mathrm{\Δt}-1)u_1^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δx}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i+m,k}^n+u_{i-m,k}^n-{2u}_{i,k}^n)\right)]\\ =\frac{1}{1+d\left(z\right)\mathrm{\Δt}}[(2-d^2\left(z\right){\mathrm{\Δt}}^2)u_2^n+(d\left(z\right)\mathrm{\Δt}-1)u_2^{n-1}+v^2\frac{{\mathrm{\Δt}}^2}{{\mathrm{\Δz}}^2}\left(\underset{m=1}{\overset{M/2}{\mathrm{\Σ}}}{c}_m(u_{i,k+m}^n+u_{i,k-m}^n-{2u}_{i,k}^n)\right)]\end{array}$

Wherein, u is acoustic wavefield; I and k is two-dimensional grid sequence number; Δ t is time-sampling step-length; Δ x is x direction sizing grid; Δ z is z direction sizing grid; M is M rank, space precision; c _mit is difference coefficient; N is time slice, and d (x) and d (z) are respectively the perfectly matched layer border attenuation term relevant with x and z directional derivative.

5. high-precision spatial according to claim 1 and time any multiple variable grid finite difference forward modeling method, it is characterized in that, in step 4, when changing time step, a transitional zone is set in large time-sampling step-length district, utilize large time-sampling step-length to obtain the wave field value of subsequent time in transitional zone region, employing interpolation method obtains each moment wave field value needed for little time sampling step length calculating.

6. high-precision spatial according to claim 5 and time any multiple variable grid finite difference forward modeling method, it is characterized in that, in step 4, time-sampling step-length is any multiple, and time-sampling step-length interpolation method carries out high-order nonlinear interpolation to multiple spot time-sampling step-length for using Lagrangian simultaneously.