CN106570201A - Numerical simulation method and system for viscoelastic medium - Google Patents

Abstract

The invention discloses a numerical simulation method and system for a viscoelastic medium. The method comprises the steps: obtaining time steps under the stability conditions in a plurality of speed regions based on a viscoelastic medium model; carrying out the numerical simulation based on the time steps of the plurality of speed regions and a fluctuation equation, and obtaining wave fields of the plurality of speed regions; carrying out the boundary processing based on the time steps of the plurality of speed regions, and obtaining the wave field of the boundary of the plurality of speed regions; and obtaining the wave field of the whole field based on the seismic wave fields of the plurality of speed regions and the seismic wave field of the boundary.

Description

Visco-elastic medium method for numerical simulation and system

Technical field

The present invention relates to field of seismic exploration, more particularly, to a kind of visco-elastic medium method for numerical simulation and System.

Background technology

In field of seismic exploration, earthquake numerical simulation is seismic prospecting and seismographic important foundation, while Be understand complex dielectricss in seimic wave propagation rule important tool, its effect through whole earthquake-capturing, In processing and explaining.With going deep into for Seismic exploration and development, conventional elastic fluid theory is difficult to meet actual Medium demand.When seismic wave is propagated in actual formation, energy and phase place all change, and directly affect ground The resolution of shake data, actual formation show certain viscous-elastic behaviour, therefore by visco-elastic medium Seismic wave carries out numerical simulation, the decay characteristics come during research and analysis seimic wave propagation, to practically The raising of shake data resolution is very meaningful.But ensure that visco-elastic medium numerical simulation is stably a urgent need The problem (i.e. the determination of time step) of solution, is related to the success or failure of simulation algorithm, is also related to simulation precision Height.Characterize the viscosity property of real medium presently mainly with viscoelastic solid model, most widely used, Most popular is Kelvin-Fo Gete (Kelvin-Voigt) solid model.

Inventor has found, in the prior art, usually on the basis of elastic stability condition, by tentative calculation Determining the time step of visco-elastic medium numerical simulation, its shortcoming is that the parameter for being given can not be optimal to method, And then extra operand is brought, cause the inefficient of simulation, or even time of occurrence step-length to can not meet fortune Steady-state conditionss are calculated, partway self-interrupting situation is run.Therefore, it is necessary to develop a kind of efficient viscoelastic Dielectric value analogy method and system.

The information for being disclosed in disclosure background section is merely intended to deepen general background technology of this disclosure Understanding, and be not construed as recognizing or imply in any form that the information structure has been people in the art Prior art well known to member.

The content of the invention

The present disclosure proposes a kind of visco-elastic medium method for numerical simulation and system, which can be by friction speed In region, reasonable selection different time step-length carries out visco-elastic medium numerical simulation, realizes efficient visco-elastic medium Numerical simulation.

According to the one side of the disclosure, it is proposed that a kind of visco-elastic medium method for numerical simulation, methods described can be with Including：Based on viscoelastic solid model, the time step under the stability condition of multiple velocity bands is obtained； Time step and wave equation based on multiple velocity bands carries out numerical simulation, obtains multiple velocity bands Wave field；Time step based on multiple velocity bands carries out boundary process, obtains the boundary of multiple velocity bands The wave field at place；And the seismic wave field of the seismic wave field and intersection based on multiple velocity bands, obtain whole The wave field in region.

According to another aspect of the present disclosure, it is proposed that a kind of visco-elastic medium numerical simulation system, the system can To include：For based on viscoelastic solid model, obtaining the time under the stability condition of multiple velocity bands The unit of step-length；Numerical simulation is carried out for the time step and wave equation based on multiple velocity bands, is obtained Obtain the unit of the wave field of multiple velocity bands；Intersection is carried out for the time step based on multiple velocity bands Reason, obtains the unit of the wave field of the intersection of multiple velocity bands；And for based on multiple velocity bands The seismic wave field of seismic wave field and intersection, obtains the unit of the wave field of whole region.

Disclosed method and device have other characteristics and advantage, and these characteristics and advantage are from being expressly incorporated herein In accompanying drawing and subsequent specific embodiment in will be apparent, or by the accompanying drawing being incorporated herein With stated in subsequent specific embodiment in detail, these the drawings and specific embodiments be provided commonly for explain this Disclosed certain principles.

Description of the drawings

Disclosure exemplary embodiment is described in more detail by combining accompanying drawing, the disclosure it is above-mentioned with And other purposes, feature and advantage will be apparent from, wherein, in disclosure exemplary embodiment, Identical reference number typically represents same parts.

The flow chart of the step of Fig. 1 shows the visco-elastic medium method for numerical simulation according to the disclosure.

Fig. 2 shows the schematic diagram of Kelvin-Fo Gete (Kelvin-Voigt) viscoelastic solid model.

Fig. 3 a and 3b show each of the visco-elastic medium method for numerical simulation of the one embodiment according to the disclosure Plant the schematic diagram of wavefield component and physical parameter.

Fig. 4 shows the schematic diagram of the intersection in the friction speed region of the one embodiment according to the disclosure.

Fig. 5 a-5c show the schematic diagram for carrying out visco-elastic medium simulation respectively with three time steps.

Fig. 6 shows the schematic diagram of the rate pattern of the one embodiment according to the disclosure.

Fig. 7 a and 7b show the different time step-length analog result and phase of the one embodiment according to the disclosure With the schematic diagram of the comparison of time step analog result.

Specific embodiment

Preferred embodiment of the present disclosure is more fully described below with reference to accompanying drawings.Although this is shown in accompanying drawing Disclosed preferred embodiment, however, it is to be appreciated that may be realized in various forms the disclosure and should not be by here The embodiment of elaboration is limited.Conversely, these embodiments are provided so that the disclosure is more thorough and complete, And the scope of the present disclosure intactly can be conveyed to those skilled in the art.

Embodiment 1

The flow chart of the step of Fig. 1 shows the visco-elastic medium method for numerical simulation according to the disclosure.

In this embodiment, can be included according to the visco-elastic medium method for numerical simulation of the disclosure：Step 101, Based on viscoelastic solid model, the time step under the stability condition of multiple velocity bands is obtained；Step 102, Time step and wave equation based on multiple velocity bands carries out numerical simulation, obtains multiple velocity bands Wave field；Step 103, the time step based on multiple velocity bands carry out boundary process, obtain multiple speed areas The wave field of the intersection in domain；And step 104, the ground of seismic wave field and intersection based on multiple velocity bands Seismic wave field, obtains the wave field of whole region.

The embodiment is based on viscoelastic solid model, and under the stability condition in friction speed region, reasonable selection is not Same time step carries out visco-elastic medium numerical simulation, realizes efficient visco-elastic medium numerical simulation.

The following detailed description of the concrete steps of the visco-elastic medium method for numerical simulation according to the disclosure.

Obtain the time step under stability condition

In one example, viscoelastic solid model can be based on, obtains the stability bar in multiple velocity bands Time step under part.

In one example, the viscoelastic solid model can be Kelvin-Fo Gete (Kelvin-Voigt) viscoelastic Solid model.It is to be understood by those skilled in the art that the viscoelastic solid model of the disclosure is not limited to this, Various viscoelastic solid models known in the art can be adopted.

Fig. 2 shows the schematic diagram of Kelvin-Fo Gete (Kelvin-Voigt) viscoelastic solid model.Such as Fig. 2 Shown, Kelvin-Voigt viscoelastic solids model can (coefficient of elasticity be M by an elastomer₁) and antivibrator (coefficient of viscosity is M₂) compose in parallel, the viscous-elastic behaviour of medium can be characterized after being constructed such that.

In one example, can based on the Kelvin-Voigt viscoelastic solid models set up, according to The composition of Kelvin-Voigt viscoelastic solid models, derives its constitutive equation.Kelvin-Voigt has been illustrated below The constitutive equation of viscoelastic solid model：

Wherein, P can represent the stress of viscoelastic solid, σ₁For can be with elastic stress, σ₂Can answer for viscous Power, t can be with express times, and ε can be the displacement (strain) that spring and antivibrator are produced.As shown in Figure 2, The displacement that spring and antivibrator are produced is equal.

It is then possible to set up stability equations.

Relation between strain stress and particle displacement (u, v, w) can be expressed as：

Wherein, x, y, z are mutually perpendicular to, and can represent three directions of particle displacement respectively.

Formula (2) is substituted into formula (1), can be obtained：

Secondary partial derivative is asked to t simultaneously to formula (3) the right and left, can be obtained：

For Na Weier (Navier) equation of sound wave can be expressed as：

Wherein, ρ can represent the density of viscoelastic solid.

Formula (5) is substituted in formula (4), can be obtained：

Stress P in formula (6) is taken into spatial fourier transform, the expression formula that can be obtained in wave-number domain is：

Wherein,Can be the spatial fourier transform of P, k can be wave number.By in formula (7) when Between derivative difference approximation, can obtain after arrangement：

Wherein(n-1)th, n can be respectively, the n+1 momentValue, Δ t can be between express time samplings Every (that is, time step).

It is then possible to makeThen have：

Because,Wherein v can be the speed of seismic wave, and ω can To represent circular frequency, Q can be quality factor.

Such that it is able to obtain：

Obviously, stability condition can be the mould of the eigenvalue λ of state-transistion matrix A less than 1.

Wherein,

Thus, it is possible to the stability condition for obtaining Kelvin-Voigt models is：

Wherein, v ' can represent visco-elastic medium seismic wave propagation speed, and have ω=2 π f, f are frequency.In the case where space difference accuracy is 2N, if the grid on tri- directions of x, y, z Step-length is respectively Δ x, and Δ y, Δ z, wavenumber components are respectively k_x,k_y,k_z.Can draw：

Wherein a_lCan be the difference coefficient of 2N precision, the span of l can be [- N, N].So, can To have：

Formula (11) is substituted into formula (10), obtaining steady-state conditionss expression formula can be：

In one example, in the case of the Difference Calculation of two-dimensional wave equation, viscoelastic solid model is stablized Property condition is：

Wherein, Δ t express times step-length, h represent numerical simulation mesh spacing, and v' represents visco-elastic medium seismic wave Speed, ω represent circular frequency, and Q represents quality factor, and N represents difference accuracy, a_lFor the difference of 2N precision Coefficient, the span of l is [- N, N].In realistic simulation, two-dimensional wave equation can be generally carried out Difference Calculation, as Δ x=Δ z=h, stability condition can be reduced to formula (13).

As set forth above, it is possible to obtain the stability condition formula of viscoelastic solid model.Meeting friction speed area In the case of the stability condition in domain, the time step in friction speed region can be manually set.

Obtain the wave field of multiple velocity bands

In one example, Numerical-Mode can be carried out based on the time step of multiple velocity bands and wave equation Intend, obtain the wave field of multiple velocity bands.

Specifically, in inhomogeneous isotropic elastic medium, two dimension glues the single order stress-rate equation shape of ACOUSTIC WAVE EQUATION Formula can be expressed as follows shown in formula (13).

Wherein, x, z can represent horizontal direction and vertical direction, v respectively_x、v_zCan represent in x, z directions Speed, u can represent density, v with express time, ρ with representation uniaxial stress, t_pCan represent Velocity of longitudinal wave, v_p' visco-elastic medium velocity of longitudinal wave can be represented.And can haveWherein ω is circle Frequency, Q_pFor compressional wave quality factor.

Fig. 3 a and 3b show each of the visco-elastic medium method for numerical simulation of the one embodiment according to the disclosure Plant the schematic diagram of wavefield component and physical parameter.Various wavefield components and physical parameter are according to Fig. 3 a and 3b point Cloth, in conjunction with the difference thinking of staggered-mesh, can obtain 2N ranks space difference accuracy, second-order time difference Precision staggered-mesh high-order finite difference method form, then difference scheme such as formula (15)-(17) of formula (14) It is shown.

Wherein,

Wherein, u^t(i, j) can represent wave field of the mesh point (i, j) in time t, u^t+(i, j) represents that mesh point (i, j) exists The wave field of the later moment in time of time t, u^t-(i, j) can represent ripple of the mesh point (i, j) in the previous moment of time t .X direction velocity component of the mesh point (i, j) in time t can be represented,Can represent Mesh point (i⁺, j) in the x directions velocity component of time t.Mesh point (i, j) can be represented time t's Z directions velocity component,Represent in mesh point (i, j⁺) in the z directions velocity component of time t.Δ t is represented Time step, Δ x, Δ z represent x, the mesh spacing in z directions, v_x、v_zMesh point is represented in x, the speed in z directions, ρ represents density, v_pRepresent velocity of longitudinal wave,X directions forward difference is represented,X directions backward difference is represented,Z directions forward difference is represented,Z directions backward difference is represented,Represent to time t forward difference, Represent to time t backward difference.

In formula (15)-(17), can meet：

During formula (18) and (19) show formula (15)Solution Process.Wherein, N is difference accuracy；For the difference weight coefficient of correspondence m rank precision, the value model of m Enclose for 1,2 ..., N.

According to simulated time t, the maximum displacement of whole seimic wave propagation during each mesh point can be counted, And retain the maximum displacement of each mesh point.After numerical simulation terminates, can obtain obtaining multiple velocity bands Wave field.

Obtain the wave field of the intersection of multiple velocity bands

In one example, boundary process can be carried out based on the time step of multiple velocity bands, is obtained many The wave field of the intersection of individual velocity band.

In one example, the time step based on multiple velocity bands carries out boundary process can include：Protect Stay the wave field in transitional region.

In the visco-elastic medium method for numerical simulation according to the disclosure, for multiple different velocity bands are adopted Multiple different time steps, it is therefore desirable to which boundary process is carried out to the intersection of multiple velocity bands, with Solve the problems, such as the data transfer of friction speed region intersection.

Detailed process can be with as described below：

Fig. 4 shows the schematic diagram of the intersection in the friction speed region of the one embodiment according to the disclosure. It is assumed that the time step Δ t of a determination, starts on time point n Δ t.First, intermediate zone The region of change must be retained (gray shade part).Then, in zone 1, time step is with 5 Δ t's Time step is spaced to calculate, and time stage (n-3) Δ t and n Δ t is merged.Time step quilt in zone 1 Complete.

Then, the time integral on the time step interval of 1 Δ t, 2 Δ t, 3 Δ t and 4 Δ t, with time stage n Δs t and (n-1) Δ t, time stage n Δs t and (n-2) Δ t, time stage n Δs t and (n-3) Δ t and time stage n Δs t and (n-4) Δ t are in mistake Cross and calculated in band.Here us are given regional value in transit time (the empty center of circle, the real center of circle, rectangle and Rhombus), the wave field value of these time steps is needed in region 2.Because the value of mistake is beyond zone boundary Issuable, the time step in region 2 is calculated and is spaced 1 Δ t until time stage with time step (n+5) Δ t is reached (time step in region 2 is the 1/5 of region 1).

The numerical value of intermediate zone zone line is needed with little need for extra computing capability, only extra storage It is allocated in order to the change of holding area is in the interlude level (gray shade in Fig. 3 near zone boundary Region), and transition region is very narrow, therefore its internal memory for being consumed and increased amount of calculation can be neglected substantially Slightly.

Additionally, if staggered-mesh, intermediate zone width needs to double, because each time step provides two Individual first order spatial derivative.

As set forth above, it is possible to obtain the wave field of the intersection of multiple velocity bands.

Obtain the wave field of whole region

In one example, can based on the seismic wave field of multiple velocity bands and the seismic wave field of intersection, Obtain the wave field of whole region.

Using example

For ease of understanding the scheme and its effect of embodiment of the disclosure, concrete application example given below.This Art personnel should be understood that the example only for the purposes of understanding the present invention, and its any detail is not anticipated The present invention is being limited by any way.

Example 1：Stability condition quantification function analysis

Be given below a two-dimensional example carry out visco-elastic medium stability condition quantitatively determine to calculate effect shadow Ring analysis.

The model of the example is homogeneous half space model, and speed is 3000m/s, and model size is 1000m*1000m, space interval Δ x=Δ z=5m, quality factor q=5, the difference coefficient point of 6 rank spatial accuracies It is not：a₀=-2.7222222, a_-1=a₁=1.5000000, a_-2=a₂=-0.15000000, a_-3=a₃=0.011111111.

Based on the parameter for providing, the time step for meeting stable simulation computing is calculated by formula (13) Dt=0.422ms (that is, Δ t=0.422ms).Dt=0.38ms computings can be quantitatively selected, relatively be managed By value, amount of calculation less, can guarantee that and calculate successfully.Two time steps dt=0.45ms manually can be set And dt=0.28ms, visco-elastic medium simulation is carried out with three time steps.

Fig. 5 a-5c show the schematic diagram for carrying out visco-elastic medium simulation respectively with three time steps, wherein, Fig. 5 a are dt=0.45ms operation results, and Fig. 5 b are dt=0.38ms operation results, and Fig. 5 c are dt=0.28ms Operation result.Table 1 shows the corresponding analysis result of different time step-length.As illustrated in figs. 5 a-5 c, and join 1 is shown in Table, if it can be found that computing is gone using dt=0.45ms, fail certainly, affect working effect, if Computing is gone using dt=0.28ms, can also be smoothly completed, but extra operand will be brought, waste a large amount of Time.

Table 1

Sampling time (dt/ms)	Cycle-index	It is stable
			0.45	2222	It is no
0.38	2632	It is
			0.28	3571	It is

Example 2：Variable time step improves simulation trial efficiency analysiss

Fig. 6 shows the schematic diagram of the rate pattern of the one embodiment according to the disclosure.

The model of the example is a horizontal layer rate pattern, and per layer of velocity amplitude is as shown in fig. 6, model area Domain grid number is 500*500, and sizing grid is 5m, and shot point coordinate is for (250,4), demarcation line is 750m Place's (ground floor).

Due to being simulated to identical model, identical algorithm calculates the used time on different configuration of computer May be different, the memory requirements of Moving grids difference algorithm has close relationship with the selection of model again, now just with As a example by stratified model shown in Fig. 6, to save as 3.49G-Byte, microcomputers of the CPU for PIV-3.0GHz interior Based on the trial result is simulated on calculation machine, contrast is analyzed.

Fig. 7 a and 7b show the different time step-length analog result and phase of the one embodiment according to the disclosure With the schematic diagram of the comparison of time step analog result.Table 2 shows the relative analyses of Fig. 7 a and 7b.

By the calculating of steady-state conditionss for the calculating time sampling interval that this model maximal rate meets is at least 1ms, the calculating time sampling interval that minimum speed meets are at least 2ms, and both sampling time intervals are all to work as Model space mesh spacing in calculating process is constant.When calculating process is all to use 1ms, stability is met Condition, calculates the used time for 4380ms to this region, and analog result is as shown in Figure 7b；When model is in calculating process Middle sampling variable time step, ground floor low velocity layer (LVL) time interval 2ms, other regions 1ms, this also meets Steady-state conditionss, calculate the used time for 3780ms to this region, and analog result as shown in Figure 7a, save by the time that calculates 13.70%.When whole model calculates time sampling interval all uses 2ms, operation time can be more saved, but It is now to be unsatisfactory for stability condition so that calculate failure (labor is shown in Table 2).It can thus be seen that If model low velocity layer (LVL) velocity amplitude is lower, region is bigger, then variable time step method carries out sound wave simulation energy It is enough that more calculating times are saved in the case where stability condition is met.

Table 2

It will be understood by those skilled in the art that the purpose of the description of embodiment of this disclosure is only for example above Property ground explanation embodiment of the disclosure beneficial effect, be not intended to by embodiment of the disclosure be limited to Any example for going out.

Embodiment 2

In accordance with an embodiment of the present disclosure, there is provided a kind of visco-elastic medium numerical simulation system, the system can be with Including：For based on viscoelastic solid model, obtaining the time step under the stability condition of multiple velocity bands Long unit；Numerical simulation is carried out for the time step and wave equation based on multiple velocity bands, is obtained The unit of the wave field of multiple velocity bands；Boundary process is carried out for the time step based on multiple velocity bands, Obtain the unit of the wave field of the intersection of multiple velocity bands；And for the earthquake based on multiple velocity bands The seismic wave field of wave field and intersection, obtains the unit of the wave field of whole region.

The embodiment is based on viscoelastic solid model, and under the stability condition in friction speed region, reasonable selection is not Same time step carries out visco-elastic medium numerical simulation, realizes efficient visco-elastic medium numerical simulation.

In one example, the viscoelastic solid model can be the extremely viscous bullet solid models of Kelvin-Fo Ge.

In one example, in the case of the Difference Calculation of two-dimensional wave equation, the viscoelastic solid model Stability condition can be：

Wherein, Δ t express times step-length, h represent numerical simulation mesh spacing, and v' represents visco-elastic medium seismic wave Speed, ω represent circular frequency, and Q represents quality factor, and N represents difference accuracy, a^lFor the difference of 2N precision Coefficient, the span of l is [- N, N].

In one example, the time step and wave equation based on multiple velocity bands carry out numerical simulation can To include：The staggered-mesh numerical simulation of mesh point (i, j) is realized by equation below：

Wherein,

Wherein, u^t(i, j) represents wave field of the mesh point (i, j) in time t, u^t+(i, j) represents mesh point (i, j) in the time The wave field of the later moment in time of t, u^t-(i, j) represents wave field of the mesh point (i, j) in the previous moment of time t； X direction velocity component of the mesh point (i, j) in time t is represented,Represent in mesh point (i⁺, j) in time t X directions velocity component；Z direction velocity component of the mesh point (i, j) in time t is represented,Table Show in mesh point (i, j⁺) in the z directions velocity component of time t；Δ t express time step-lengths, Δ x, Δ z represent x, z side To mesh spacing, v_x、v_zMesh point is represented in x, the speed in z directions, ρ represent density, v_pRepresent compressional wave Speed, v_p' visco-elastic medium velocity of longitudinal wave is represented,X directions forward difference is represented,Represent x directions backward Difference,Z directions forward difference is represented,Z directions backward difference is represented,It is poor forward to time t to represent Point,Represent to time t backward difference.

In one example, the time step based on multiple velocity bands carries out boundary process can include：Protect Stay the wave field in transitional region.

It will be understood by those skilled in the art that the purpose of the description of embodiment of this disclosure is only for example above Property ground explanation embodiment of the disclosure beneficial effect, be not intended to by embodiment of the disclosure be limited to Any example for going out.

The disclosure can be system, method and/or computer program.Computer program can include Computer-readable recording medium, containing the computer for being used to make processor realize various aspects of the disclosure Readable program instructions.

Computer-readable recording medium can be to keep and store the instruction used by instruction execution equipment Tangible device.Computer-readable recording medium can for example be-- but be not limited to-- storage device electric, Magnetic storage apparatus, light storage device, electromagnetism storage device, semiconductor memory apparatus or above-mentioned any conjunction Suitable combination.The more specifically example (non exhaustive list) of computer-readable recording medium includes：It is portable Formula computer disks, hard disk, random access memory (RAM), read only memory (ROM), erasable type can Program read-only memory (EPROM or flash memory), static RAM (SRAM), portable pressure Contracting disk read only memory (CD-ROM), digital versatile disc (DVD), memory stick, floppy disk, machinery are compiled Decoding apparatus, the punch card for being for example stored thereon with instruction or groove internal projection structure and above-mentioned any conjunction Suitable combination.Computer-readable recording medium used herein above is not construed as instantaneous signal itself, such as The electromagnetic wave of radio wave or other Free propagations, the electromagnetic wave propagated by waveguide or other transmission mediums (for example, by the light pulse of fiber optic cables) or the signal of telecommunication by wire transfer.

Computer-readable program instructions as described herein can download to each from computer-readable recording medium Calculating/processing equipment, or downloaded by network, such as the Internet, LAN, wide area network and/or wireless network To outer computer or External memory equipment.Network can include copper transmission cable, fiber-optic transfer, wireless biography Defeated, router, fire wall, switch, gateway computer and/or Edge Server.Each calculates/processes and sets Adapter or network interface in standby receives computer-readable program instructions from network, and forwards the meter Calculation machine readable program instructions, for being stored in the computer-readable recording medium in each calculating/processing equipment.

Computer program instructions for performing disclosure operation can be assembly instruction, instruction set architecture (ISA) instruction, machine instruction, machine-dependent instructions, microcode, firmware instructions, condition setup data, Or the source code write with the combination in any of one or more programming language or object code, the programming language Speech includes OO programming language-Smalltalk, C++ etc., and the process type programming language of routine Speech-such as " C " language or similar programming language.Computer-readable program instructions fully can be counted in user Perform on calculation machine, partly perform on the user computer, performing as an independent software kit, part Part performs or completely on remote computer or server on the remote computer on the user computer Perform.In the situation of remote computer is related to, remote computer can be by the network-bag of any kind LAN (LAN) or wide area network (WAN)-be connected to subscriber computer are included, or, it may be connected to outside meter Calculation machine (such as using ISP come by Internet connection).In certain embodiments, pass through Using the status information of computer-readable program instructions come personalized customization electronic circuit, such as FPGA Circuit, field programmable gate array (FPGA) or programmable logic array (PLA), the electronic circuit can To perform computer-readable program instructions, so as to realize various aspects of the disclosure.

Referring herein to the stream of the method according to the embodiment of the present disclosure, device (system) and computer program Journey figure and/or block diagram describe various aspects of the disclosure.It should be appreciated that each side of flow chart and/or block diagram The combination of each square frame in frame and flow chart and/or block diagram, can be realized by computer-readable program instructions.

These computer-readable program instructions can be supplied to general purpose computer, special-purpose computer or other can compile The processor of journey data processing equipment, so as to produce a kind of machine so that these instructions are by computer Or other programmable data processing units computing device when, generate in flowchart and/or block diagram The device of function/action specified in one or more square frames.Can also be these computer-readable program instructions Storage in a computer-readable storage medium, these instruct cause computer, programmable data processing unit and/ Or other equipment is worked in a specific way, so as to the computer-readable medium of the instruction that is stored with then includes one Manufacture, which includes flowchart and/or function/action specified in one or more square frames in block diagram The instruction of various aspects.

Can also computer-readable program instructions be loaded into computer, other programmable data processing units, Or on miscellaneous equipment so that one is performed on computer, other programmable data processing units or miscellaneous equipment Series of operative steps, to produce computer implemented process, so that in computer, other programmable numbers According to one or more sides in the instruction flowchart and/or block diagram performed in processing meanss or miscellaneous equipment Function/action specified in frame.

Flow chart and block diagram in accompanying drawing shows system, method and the meter of the multiple embodiments according to the disclosure The architectural framework in the cards of calculation machine program product, function and operation.At this point, flow chart or block diagram In each square frame can represent the part of module, program segment or an instruction, the module, program segment Or a part for instruction is used for the executable instruction of the logic function that realization specifies comprising one or more.Having In a little realizations as replacement, the function of being marked in square frame can also be with suitable different from what is marked in accompanying drawing Sequence occurs.For example, two continuous square frames can essentially be performed substantially in parallel, and they sometimes can also Perform in the opposite order, this is depending on involved function.It is also noted that block diagram and/or flow chart In each square frame and block diagram and/or the square frame in flow chart combination, can use perform regulation function Or the special hardware based system of action is realizing, or specialized hardware can be used with computer instruction Combine to realize.

It is described above the presently disclosed embodiments, described above is exemplary, and non-exclusive, And it is also not necessarily limited to disclosed each embodiment.In the scope and spirit without departing from illustrated each embodiment In the case of, many modifications and changes will be apparent from for those skilled in the art. The selection of term used herein, it is intended to best explain the principle of each embodiment, practical application or to market In technology improvement, or make the art other those of ordinary skill be understood that it is disclosed herein each Embodiment.

Claims (10)

1. a kind of visco-elastic medium method for numerical simulation, including：

Based on viscoelastic solid model, the time step under the stability condition of multiple velocity bands is obtained；

Time step and wave equation based on multiple velocity bands carries out numerical simulation, obtains multiple speed areas The wave field in domain；

Time step based on multiple velocity bands carries out boundary process, obtains the intersection of multiple velocity bands Wave field；And

The seismic wave field of seismic wave field and intersection based on multiple velocity bands, obtains the wave field of whole region.

2. visco-elastic medium method for numerical simulation according to claim 1, wherein, the viscoelastic solid mould Type is the extremely viscous bullet solid models of Kelvin-Fo Ge.

3. visco-elastic medium method for numerical simulation according to claim 2, wherein, in two-dimensional wave equation Difference Calculation in the case of, the stability condition of the viscoelastic solid model is：

$\mathrm{\Δ}t\≤\frac{h}{v^{\′}}\sqrt{\frac{2}{\underset{l=-N}{\overset{N}{\Σ}}\left|a_l\right|}+\frac{v^{\′2}}{Q^2{\mathrm{\ω}}^2{h}^2}}-\frac{1}{Q\mathrm{\ω}}$

Wherein, Δ t express times step-length, h represent numerical simulation mesh spacing, and v' represents visco-elastic medium seismic wave Speed, ω represent circular frequency, and Q represents quality factor, and N represents difference accuracy, a_lFor the difference of 2N precision Coefficient, the span of l is [- N, N].

4. visco-elastic medium method for numerical simulation according to claim 1, wherein, based on multiple speed areas The time step and wave equation in domain carries out numerical simulation to be included：The friendship of mesh point (i, j) is realized by equation below Wrong grid values simulation：

$\begin{array}{c}{u}^{t+}(i,j)=u^{t-}(i,j)-\frac{{\mathrm{\Δt\ρv}}_p^2}{\mathrm{\Δ}x}\left\{L_x^{-}\[v_x^t\left(i^{+},j\right)\]\right\}-\frac{{\mathrm{\Δt\ρv}}_p^2}{\mathrm{\Δ}z}\left\{L_z^{-}\[v_z^t\left(i,j^{+}\right)\]\right\}\\ -\frac{{\mathrm{\Δt\ρv}}_p^{\′2}}{\mathrm{\Δ}x}\left\{L_t^{-}\[v_x^t\left(i^{+},j\right)\]\right\}-\frac{{\mathrm{\Δt\ρv}}_p^{\′2}}{\mathrm{\Δ}z}\left\{L_t^{-}\[v_z^t\left(i,j^{+}\right)\]\right\}\end{array}$

$v_x^t(i^{+},j)=v_x^{t-1}(i^{+},j)-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{L_x^{+}\[u^{t-}(i,j)\]\}$

$v_z^t(i,j^{+})=v_z^{t-1}(i,j^{+})-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{L_z^{+}\[u^{t-}(i,j)\]\}$

Wherein, $t+=t+\frac{1}{2}\mathrm{\Δ}t,t-=t-\frac{1}{2}\mathrm{\Δ}t,t-1=t-\mathrm{\Δ}t,i^{+}=i+\frac{1}{2}\mathrm{\Δ}x,i^{-}=i-\frac{1}{2}\mathrm{\Δ}x,j^{+}=j+\frac{1}{2}\mathrm{\Δ}z,$ $j^{-}=j-\frac{1}{2}\mathrm{\Δ}z,$

Wherein, u^t(i, j) represents wave field of the mesh point (i, j) in time t, u^t+(i, j) represents mesh point (i, j) in the time The wave field of the later moment in time of t, u^t-(i, j) represents wave field of the mesh point (i, j) in the previous moment of time t； X direction velocity component of the mesh point (i, j) in time t is represented,Represent in mesh point (i⁺, j) in time t X directions velocity component；Z direction velocity component of the mesh point (i, j) in time t is represented,Table Show in mesh point (i, j⁺) in the z directions velocity component of time t；Δ t express time step-lengths, Δ x, Δ z represent x, z side To mesh spacing, v_x、v_zMesh point is represented in x, the speed in z directions, ρ represent density, v_pRepresent compressional wave Speed, v'_pVisco-elastic medium velocity of longitudinal wave is represented,X directions forward difference is represented,Represent x directions backward Difference,Z directions forward difference is represented,Z directions backward difference is represented,It is poor forward to time t to represent Point,Represent to time t backward difference.

5. visco-elastic medium method for numerical simulation according to claim 1, wherein, based on multiple speed areas The time step in domain carries out boundary process to be included：Retain the wave field in transitional region.

6. a kind of visco-elastic medium numerical simulation system, including：

For based on viscoelastic solid model, obtaining the time step under the stability condition of multiple velocity bands Unit；

Numerical simulation is carried out for the time step and wave equation based on multiple velocity bands, multiple speed are obtained The unit of the wave field in degree region；

Boundary process is carried out for the time step based on multiple velocity bands, the friendship of multiple velocity bands is obtained The unit of the wave field at boundary；And

For the seismic wave field of seismic wave field and intersection based on multiple velocity bands, whole region is obtained The unit of wave field.

7. visco-elastic medium numerical simulation system according to claim 6, wherein, the viscoelastic solid mould Type is the extremely viscous bullet solid models of Kelvin-Fo Ge.

8. visco-elastic medium numerical simulation system according to claim 7, wherein, in two-dimensional wave equation Difference Calculation in the case of, the stability condition of the viscoelastic solid model is：

$\mathrm{\Δ}t\≤\frac{h}{v^{\′}}\sqrt{\frac{2}{\underset{l=-N}{\overset{N}{\Σ}}\left|a_l\right|}+\frac{v^{\′2}}{Q^2{\mathrm{\ω}}^2{h}^2}}-\frac{1}{Q\mathrm{\ω}}$

Wherein, Δ t express times step-length, h represent numerical simulation mesh spacing, and v' represents visco-elastic medium seismic wave Speed, ω represent circular frequency, and Q represents quality factor, and N represents difference accuracy, a_lFor the difference of 2N precision Coefficient, the span of l is [- N, N].

9. visco-elastic medium numerical simulation system according to claim 6, wherein, based on multiple speed areas The time step and wave equation in domain carries out numerical simulation to be included：The friendship of mesh point (i, j) is realized by equation below Wrong grid values simulation：

$\begin{array}{c}{u}^{t+}(i,j)=u^{t-}(i,j)-\frac{{\mathrm{\Δt\ρv}}_p^2}{\mathrm{\Δ}x}\left\{L_x^{-}\[v_x^t\left(i^{+},j\right)\]\right\}-\frac{{\mathrm{\Δt\ρv}}_p^2}{\mathrm{\Δ}z}\left\{L_z^{-}\[v_z^t\left(i,j^{+}\right)\]\right\}\\ -\frac{{\mathrm{\Δt\ρv}}_p^{\′2}}{\mathrm{\Δ}x}\left\{L_t^{-}\[v_x^t\left(i^{+},j\right)\]\right\}-\frac{{\mathrm{\Δt\ρv}}_p^{\′2}}{\mathrm{\Δ}z}\left\{L_t^{-}\[v_z^t\left(i,j^{+}\right)\]\right\}\end{array}$

$v_x^t(i^{+},j)=v_x^{t-1}(i^{+},j)-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x\mathrm{\ρ}}\{L_x^{+}\[u^{t-}(i,j)\]\}$

$v_z^t(i,j^{+})=v_z^{t-1}(i,j^{+})-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}z\mathrm{\ρ}}\{L_z^{+}\[u^{t-}(i,j)\]\}$

Wherein, $t+=t+\frac{1}{2}\mathrm{\Δ}t,t-=t-\frac{1}{2}\mathrm{\Δ}t,t-1=t-\mathrm{\Δ}t,i^{+}=i+\frac{1}{2}\mathrm{\Δ}x,i^{-}=i-\frac{1}{2}\mathrm{\Δ}x,j^{+}=j+\frac{1}{2}\mathrm{\Δ}z,$ $j^{-}=j-\frac{1}{2}\mathrm{\Δ}z,$

Wherein, u^t(i, j) represents wave field of the mesh point (i, j) in time t, u^t+(i, j) represents mesh point (i, j) in the time The wave field of the later moment in time of t, u^t-(i, j) represents wave field of the mesh point (i, j) in the previous moment of time t； X direction velocity component of the mesh point (i, j) in time t is represented,Represent in mesh point (i⁺, j) in time t X directions velocity component；Z direction velocity component of the mesh point (i, j) in time t is represented,Table Show in mesh point (i, j⁺) in the z directions velocity component of time t；Δ t express time step-lengths, Δ x, Δ z represent x, z side To mesh spacing, v_x、v_zMesh point is represented in x, the speed in z directions, ρ represent density, v_pRepresent compressional wave Speed, v'_pVisco-elastic medium velocity of longitudinal wave is represented,X directions forward difference is represented,Represent x directions backward Difference,Z directions forward difference is represented,Z directions backward difference is represented,It is poor forward to time t to represent Point,Represent to time t backward difference.

10. visco-elastic medium numerical simulation system according to claim 6, wherein, based on multiple speed The time step in region carries out boundary process to be included：Retain the wave field in transitional region.