CN105319581A - Efficient time domain full waveform inversion method - Google Patents

Abstract

The present invention provides an efficient time domain full waveform inversion method, belonging to the field of efficient, low storage, high precision and high speed modeling of petroleum seismography data. The efficient time domain full waveform inversion method comprises the following steps: (1) employing original seismic observation records to randomly form supper-shot gathers with different frequency bands through adoption of a source polarity randomizing scheme, a source position randomizing scheme and a source number randomizing scheme; (2) accelerating the forward modeling of high order staggered mesh finite differences of a sound wave equation system and an elastic wave system of the one-form basis of the time domain through adoption of a multi-GPU parallel computing technique on the basis of CUDA, and obtaining a supper-shot gather with the ifreq frequency band simulated by the GPU; and (3) simulating propagations with different seismic sources through multi-GPU nucleus and through adoption of the collaborative technology of the CPU and the GPU at the forward modeling, and storing boundary data of a seismic source wave field and all the wave fields at the last moment to a computer memory.

Description

A kind of time domain full waveform inversion method efficiently

Technical field

The invention belongs to oil seismic exploration data efficient rate, low storage, high precision velocity modeling field, be specifically related to a kind of time domain full waveform inversion method efficiently, for huge calculated amount and mass memory problem needed for the modeling of time domain full wavefield velocity, adopt focus coding, parallel computing strategy that many GPU combine and the countless efficiency frontier storage policy according to I/O completes efficiently, high resolving power velocity modeling.

Background technology

Utilize pre-stack seismic wave datum kinematics and dynamic information, full waveform inversion has the ability of the tectonic structure disclosed under complex geological condition, lithologic parameter, and it is considered to the highest technology of velocity modeling precision at present.But the realization of three-dimensional full-wave shape inverting is also subject to all restrictions, wherein two Major Difficulties are that huge calculated amount and memory space are solved always preferably.

Because the correct enforcement of full waveform inversion needs hundreds and thousands of forward simulations, thus the efficiency of forward simulation is the essential place of restriction full waveform inversion development.To two dimensional inversion problem, adopt the frequency field Wave equation forward modeling technology of direct solution efficiency when processing the simulation of many focus higher, but to three-dimensional problem, even if adopt the nested subdivision direct solution that efficiency is higher, the computation complexity that frequency field is just being drilled also will reach n _sn _wo (n ⁴logn), required memory O (n ⁵), its calculated amount and memory space huge, current computing machine condition be difficult to realize.To D integral pin-fin tube, the Time Calculation complexity that time finite element method is just being drilled is O (n ⁴), required memory complexity is O (n ³), just drill relative to three-dimensional frequency field, time domain is just drilling that to have calculated amount little, the advantages such as EMS memory occupation is low.At present two kinds of thinkings are mainly contained to the Study on Acceleration of forward simulation: one, by the acceleration of computer hardware to improve the counting yield (the GPU parallel computation etc. as based on CUDA platform) of inverting, two, the counting yield of algorithm itself is improved, many source inversions technical thought based on random phase encoding is better, its basic thought is that all single source data are superimposed as single or multiple super source data, and processes super source data.Because the counting yield of full waveform inversion and the number of focus are approximated to direct ratio, thus adopt super focus greatly can improve the counting yield of full waveform inversion.

The another one difficult problem that Conventional Time territory full waveform inversion faces needs mass data to input and stores, and the input and output (I/O) of mass data not only propose higher challenge to computer hardware, and have had a strong impact on the counting yield of inverting.Even if when adopting GPU parallel computing to accelerate the Conventional Time territory full waveform inversion that need export data, improved efficiency is limited.Moussa (2009) adopts GPUPU to show the data transmission of the working time of 91% for CPU internal memory and GPU video memory when studying reverse-time migration, and this is because the source wavefield in each moment need be stored on hard disk, caused by the parallel efficiency limiting GPU.

Summary of the invention

The present invention is directed to calculated amount needed for time domain full waveform inversion and the huge problem of memory space, provide a kind of time domain full waveform inversion method efficiently, the high precision velocity modeling for large-scale data provides a kind of process means efficiently.

The present invention is achieved by the following technical solutions:

The first step, adopt that random, the source position of source polarity is random, number random three randomizing schemes in source are randomly formed different frequency bands to original observation seismologic record and surpass big gun collection;

Second step, adopt high-order staggered-mesh finite difference simulation based on the ACOUSTIC WAVE EQUATION system of many GPU parallel computing acceleration time territory single order form of CUDA and Time Migration of Elastic Wave Equation system, the i-th freq frequency band obtaining GPU simulation surpasses big gun collection; Adopt CPU and GPU coordination technique when forward simulation, the focus different with many GPU nuclear mockup is propagated, and the data boundary of source wavefield and last moment all wave fields are stored in calculator memory; Wave field value and all wave field value of maximum moment of the Particle Vibration Velocity component each moment interior zone and PML region had a common boundary when calculating the source wavefield of forward-propagating are stored in GPU video memory;

3rd step, calculating wave field residual error and cost functional value (subtracting each other wave field residual error superposes cost functional value of suing for peace to obtain with the super big gun collection of simulating and original super big gun collection), then judge whether to meet end condition, if, represent that the full waveform inversion iteration of this frequency range terminates, and using the initial velocity model of the model of inverting as next frequency band full waveform inversion, proceed to the 6th step, if not, then proceed to the 4th step;

4th step, adopt many GPU parallel computing to calculate wave detector anti-pass wave field, adopt effective stored boundary strategy and many GPU parallel computing reconstruct source wavefield simultaneously and calculate the gradient of cost functional relative to model parameter;

5th step, adopt gradient class iterative algorithm renewal speed model, super big gun record is calculated to the rate pattern after upgrading and and the super big gun record of the first step do residual error;

6th step, judge whether to meet frequency end condition, if so, then proceed to the 7th step, if not, then return the first step;

7th step, output "current" model.

The described first step is achieved in that

Adopt the MPI original single shot record of parallel input and original Ricker wavelet, utilize S filter to original single shot record and original Ricker wavelet filtering, obtain earthquake single shot record and the seismic wavelet of different frequency bands;

Function, random function and the sign function producing random seed in C language built-in function is utilized to give [-1 to all seismic wavelets, 1] number between, namely the color that each seismic wavelet is different is given, according to the positive negativity of seismic wavelet color, give seismic wavelet corresponding polarity, this process is that source polarity is random;

Heap algorithm or the color of quick sort to the seismic wavelet changing source polarity is adopted to sort, color value increases gradually, and adopting the random function in C language to go out the number of single seismic wavelet in each super focus at random according to total super source wavelet number and total single seismic wavelet number, this process is that source number is random;

According to the number of the single seismic wavelet of each super focus, the single seismic wavelet changing source polarity after sequence is superposed, the position of each seismic wavelet now superposed is also at random corresponding, this process is that source position is random, completes the super focus that three randomizing schemes form different frequency bands thus;

In the same way the super big gun collection that three randomizing schemes form different frequency bands is completed to earthquake single shot record.

The S filter adopted is:

$F_{\mathrm{Wiener}}\left(\mathrm{\ω}\right)=\frac{W_{\mathrm{Original}}^{*}\left(\mathrm{\ω}\right)W_{T\arg \mathrm{et}}\left(\mathrm{\ω}\right)}{{\left|W_{\mathrm{Original}}\left(\mathrm{\ω}\right)\right|}^2+\mathrm{\ϵ}}$

F _wieners filter, W _originaloriginal signal, W _targetfor expecting the signal obtained, ω is angular frequency, and ε chooses 10 ^-6.

Described second step comprises:

Initial velocity model is obtained by after true model large scale Gaussian smoothing;

Adopt the propagation of high-order staggering mesh finite-difference simulation wave equation in underground medium, specific as follows:

When adopting GPU to calculate, the whole domain of two dimension is considered as a kernel function, the whole domain of two dimension is divided into multiple subregion, the size of every sub regions is 16 × 16 grid nodes, every sub regions is mapped to the thread block of a CUDA, each net point adopts the individual process of thread block independently updated, adopt the parallel thread in GPU to upgrade all nodal values at each time step, thus the wave field completing whole region calculate simultaneously;

For three-dimensional data, three-dimensional data is arranged in two-dimentional sheet according to the direction of tieing up the most slowly, then according to the mode identical with 2-D data process, three-dimensional nodes is mapped to corresponding process, completes the forward simulation in whole region.

End condition in described 3rd step is as follows:

$\left|\frac{{\mathrm{\δd}}_k-\mathrm{\δ}{d}_{k-1}}{\mathrm{\δ}{d}_{k-1}}\right|\≤\mathrm{\ϵ}$

Wherein, d _krepresent the cost functional value of the wave field residual error of kth time iteration, d _k-1represent the cost functional value of the wave field residual error of kth-1 iteration, ε is a very little amount.

Described 4th step comprises:

Adopt the anti-pass wave field that adjoint equation calculating wave detector goes out;

Wave field value and all wave field value of maximum moment of the Particle Vibration Velocity component each moment interior zone and PML region had a common boundary when calculating the source wavefield of forward-propagating are stored in GPU video memory, in the maximum moment of forward-propagating, the storage data in the 3rd step are delivered to CPU internal memory by GPU video memory; In the process of inverse time anti-pass detection residual error wave field, CPU is stored wave field and copies GPU video memory to, the residual error of super big gun collection calculated according to the 3rd step based on Huygens' principle to source wavefield approximate reconstruction, reconstruct source wavefield and all adopt many GPU parallel computing to accelerate when calculating detection wave field.

Compared with prior art, the invention has the beneficial effects as follows: the numerical simulation utilizing many GPU parallel computing acceleration time territory wave equation system, the focus coding techniques of three randomizing schemes is adopted to reduce the simulation number of focus, the efficiency of time domain full waveform inversion is promoted greatly, the low problem of counting yield that the data input and output utilizing GPU to carry out full waveform inversion for tradition cause, the present invention adopts effective stored boundary strategy, at utmost reduce the input and output of data, thus the counting yield of full waveform inversion is promoted further, the final correct enforcement for time domain full waveform inversion provides a kind of efficiently feasible technical scheme.

Accompanying drawing explanation

Fig. 1: the step block diagram of the inventive method

Fig. 2: monokaryon GPU carries out the acceleration effect comparison diagram of finite-difference forward modeling to two-dimension elastic wave equation

Fig. 3: monokaryon GPU carries out the acceleration effect comparison diagram of time finite element method forward simulation to three dimensional elasticity wave equation

Fig. 4-1: true Marmousi2 model

Fig. 4-2: initial velocity of longitudinal wave model

Fig. 4-3: ACOUSTIC WAVE EQUATION 5.0Hz Multi-scale inversion result, cross-talk noise is serious

Fig. 4-4: ACOUSTIC WAVE EQUATION 15.0Hz Multi-scale inversion result cross-talk noise is suppressed

Fig. 5-1: initial Overthurst model (Lame's constant lambda)

Fig. 5-2: initial Overthurst model (Lame's constant mu)

Fig. 5-3: Time Migration of Elastic Wave Equation 5.0Hz Multi-scale inversion result (Lame's constant lambda)

Fig. 5-4: Time Migration of Elastic Wave Equation 5.0Hz Multi-scale inversion result (Lame's constant mu)

Fig. 5-5: Time Migration of Elastic Wave Equation 10.0Hz Multi-scale inversion result (Lame's constant lambda)

Fig. 5-6: Time Migration of Elastic Wave Equation 10.0Hz Multi-scale inversion result (Lame's constant mu)

Fig. 5-7: Time Migration of Elastic Wave Equation 15.0Hz Multi-scale inversion result (Lame's constant lambda)

Fig. 5-8: Time Migration of Elastic Wave Equation 15.0Hz Multi-scale inversion result (Lame's constant mu)

Embodiment

Below in conjunction with accompanying drawing, the present invention is described in further detail:

The present invention is a kind of based on focus phase encoding and many GPU parallel computing and the countless efficient time territory wave equation full waveform inversion method according to I/O, and the high precision velocity modeling for large-scale data provides a kind of process means efficiently.

As shown in Figure 1, the inventive method comprises:

The first step, adopt that random, the source position of source polarity is random, number random three randomizing schemes in source are randomly formed different frequency bands to original observation seismologic record and surpass big gun collection;

Specific as follows:

Adopt the MPI original single shot record of parallel input and original Ricker wavelet, utilize S filter to original single shot record and original Ricker wavelet filtering, obtain earthquake single shot record and the seismic wavelet of different frequency bands, the S filter adopted is:

$F_{\mathrm{Wiener}}\left(\mathrm{\ω}\right)=\frac{W_{\mathrm{Original}}^{*}\left(\mathrm{\ω}\right)W_{T\arg \mathrm{et}}\left(\mathrm{\ω}\right)}{{\left|W_{\mathrm{Original}}\left(\mathrm{\ω}\right)\right|}^2+\mathrm{\ϵ}}$

Wherein F _wieners filter, W _originaloriginal signal, W _targetfor expecting the signal obtained, ω is angular frequency, and ε is a very little value, and in order to avoid the generation of wild effect, the present invention chooses 10 ^-6.Srand () function, random function rand () and the sign function producing random seed in C language built-in function is utilized to give [-1 to all seismic wavelets, 1] number between, namely the color that each seismic wavelet is different is given, according to the positive negativity of seismic wavelet color, give seismic wavelet corresponding polarity, now complete source polarity random;

Heap algorithm or the color of quick sort to the seismic wavelet changing source polarity is adopted to sort, color value increases gradually, and adopting the random function rand () in C language to go out the number of single seismic wavelet in each super focus at random according to total super source wavelet number and total single seismic wavelet number, this process is that source number is random;

After this according to the number of the single seismic wavelet of each super focus, the single seismic wavelet changing source polarity after sequence is superposed, the position of each seismic wavelet now superposed is also at random corresponding, this process is that source position is random, complete three randomizing schemes thus and form the super focus of different frequency bands, by same implementation, different frequency bands is formed to earthquake single shot record and surpass big gun collection.

Second step, adopt high-order staggered-mesh finite difference simulation based on the ACOUSTIC WAVE EQUATION of many GPU parallel computing acceleration time territory single order form of CUDA and Time Migration of Elastic Wave Equation system.

Second step is that forward simulation Prof. Du Yucang surpasses big gun collection and the super big gun collection of the super big gun collection adopting the first step to obtain and simulation does residual error, and the wavelet wherein adopted is super source wavelet, and initial model obtains by after true model large scale Gaussian smoothing.

Forward problem is the basis of inverse problem, and the ACOUSTIC WAVE EQUATION system of the single order form that the present invention adopts is:

$\left\{\begin{array}{c}\mathrm{\ρ}\left(x\right)\frac{\∂v_{x_i}(x,t)}{\∂t}=\frac{\∂P(x,t)}{\∂x_i}+f_{x_i},x_i=x,y,z\\ \frac{\∂P(x,t)}{\∂t}=\mathrm{\κ}\left(x\right)\underset{x_i}{\mathrm{\Σ}}\frac{\∂v_{x_i}(x,t)}{\∂x_i}\end{array}\right.$

It is write as operational form:

Lv＝f

Wherein: $L=A{\∂}_x+B{\∂}_y+C{\∂}_z+D{\∂}_t,v={(v_x,v_y,y_z,P)}^T,f={(f_x,f_y,f_z,0)}^T$

$A=\left[\begin{array}{cccc}0& 0& 0& -1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ -\mathrm{\κ}& 0& 0& 0\end{array}\right],B=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& -1\\ 0& 0& 0& 0\\ 0& -\mathrm{\κ}& 0& 0\end{array}\right],C=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -1\\ 0& 0& -\mathrm{\κ}& 0\end{array}\right],D\left[\begin{array}{cccc}\mathrm{\ρ}& 0& 0& 0\\ 0& \mathrm{\ρ}& 0& 0\\ 0& 0& \mathrm{\ρ}& 0\\ 0& 0& 0& 1\end{array}\right]$

The Time Migration of Elastic Wave Equation system of the single order form that the present invention adopts is:

$\left\{\begin{array}{c}\mathrm{\ρ}\left(x\right){\stackrel{\·}{v}}_x(x,t)={\mathrm{\σ}}_{\mathrm{xx},x}(x,t)+{\mathrm{\σ}}_{\mathrm{xy},x}(x,t)+{\mathrm{\σ}}_{\mathrm{xz},z}(x,t)\\ \mathrm{\ρ}\left(x\right){\stackrel{\·}{v}}_y(x,t)={\mathrm{\σ}}_{\mathrm{xy},x}(x,t)+{\mathrm{\σ}}_{\mathrm{yy},y}(x,t)+{\mathrm{\σ}}_{\mathrm{yz},z}(x,t)\\ \mathrm{\ρ}\left(x\right){\stackrel{\·}{v}}_z(x,t)={\mathrm{\σ}}_{\mathrm{xz},x}(x,t)+{\mathrm{\σ}}_{\mathrm{yz},y}(x,t)+{\mathrm{\σ}}_{\mathrm{zz},z}(x,t)\end{array}\right.$

$\left[\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xx}}(x,t)=(\mathrm{\λ}+2\mathrm{\μ})\left(x\right)v_{x,x}(x,t)+\mathrm{\λ}(x,t)v_{y,y}(x,t)+\mathrm{\λ}\left(x\right)v_{z,z}(x,t)+f_x(x,t)\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yy}}(x,t)=\mathrm{\λ}\left(x\right)v_{x,x}\left(x\right)+(\mathrm{\λ}+2\mathrm{\μ})(x,t)v_{y,y}(x,t)+\mathrm{\λ}\left(x\right)v_{z,z}(x,t)+f_y(x,t)\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{zz}}(x,t)=\mathrm{\λ}\left(x\right)v_{x,x}\left(x\right)+\mathrm{\λ}(x,t)v_{y,y}(x,t)+(\mathrm{\λ}+2\mathrm{\μ})\left(x\right)v_{z,z}(x,t)+f_z(x,t)\end{array}\right]$

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xy}}(x,t)=\mathrm{\μ}\left(x\right)(v_{y,x}(x,t)+v_{x,y}(x,t))\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{xz}}(x,t)=\mathrm{\μ}\left(x\right)(v_{z,x}(x,t)+v_{x,z}(x,t))\\ {\stackrel{\·}{\mathrm{\σ}}}_{\mathrm{yz}}(x,t)=\mathrm{\μ}\left(x\right)(v_{z,y}(x,t)+v_{y,z}(x,t))\end{array}\right.$

It is write as operational form:

Lv＝f

Wherein:

$L=A{\∂}_x+B{\∂}_y+C{\∂}_z+D{\∂}_t,v={(v_x,v_y,v_z,{\mathrm{\σ}}_{\mathrm{xx}},{\mathrm{\σ}}_{\mathrm{yy}},{\mathrm{\σ}}_{\mathrm{zz}},{\mathrm{\σ}}_{\mathrm{xy}},{\mathrm{\σ}}_{\mathrm{xz}},{\mathrm{\σ}}_{\mathrm{yz}})}^T,f={(\mathrm{0,0,0},f_x,f_y,f_z,\mathrm{0,0,0})}^T$

$A=\left[\begin{array}{ccccccccc}0& 0& 0& -1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -1& 0\\ -\mathrm{\λ}-2\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0& 0\\ -\mathrm{\λ}& 0& 0& 0& 0& 0& 0& 0& 0\\ -\mathrm{\λ}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -\mathrm{\μ}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right],B=\left[\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -1\\ 0& -\mathrm{\λ}& 0& 0& 0& 0& 0& 0& 0\\ 0& -\mathrm{\λ}-2\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0\\ 0& -\mathrm{\λ}& 0& 0& 0& 0& 0& 0& 0\\ -\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -\mathrm{\μ}& 0& 0& 0& 0& 0& 0\end{array}\right]$

$C=\left[\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -1\\ 0& 0& 0& 0& 0& -1& 0& 0& 0\\ 0& 0& -\mathrm{\λ}& 0& 0& 0& 0& 0& 0\\ 0& 0& -\mathrm{\λ}& 0& 0& 0& 0& 0& 0\\ 0& 0& -\mathrm{\λ}-2\mathrm{\μ}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -\mathrm{\μ}& 0& 0& 0& 0& 0& 0& 0\end{array}\right],D=\left[\begin{array}{ccccccccc}\mathrm{\ρ}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \mathrm{\ρ}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \mathrm{\ρ}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

The present invention adopts high-order staggering mesh finite-difference to simulate the propagation of wave equation in underground medium, when wave field extrapolation, the wave field of each grid node uses the wave field value in half moment before the wave field value in this node previous moment and surroundings nodes when upgrading, this problem has good concurrency, when adopting GPU to calculate, the whole domain of two dimension, (first step obtains super big gun collection, second step simulates super big gun collection with GPU, object asks both residual errors) be considered as a kernel function, the whole domain of two dimension is divided into multiple subregion, the size of every sub regions is 16 × 16 grid nodes, every sub regions is mapped to the thread block of a CUDA, each net point adopts the individual process of thread block independently updated, adopt the parallel thread in GPU to upgrade all nodal values at each time step simultaneously, thus the wave field completing whole region calculates.For three-dimensional problem, three-dimensional data is arranged in two-dimentional sheet according to the direction of tieing up the most slowly (i.e. the y-axis direction of cartesian coordinate system), processing mode based on two-dimensional problems divides process block and process, three-dimensional nodes is mapped to corresponding process, completes the forward simulation in whole region.Because the numerical simulation of each focus (a corresponding node of focus) is separate, thus can adopt CPU/GPU coordination technique when forward simulation, the focus different with many GPU nuclear mockup is propagated.

3rd step, when utilizing GPU speed-up computation source wavefield (second step) data boundary of source wavefield and last moment all wave fields are stored in calculator memory (" just drill (CPU) and deposit border wave field " in Fig. 1);

4th step, many GPU parallel computing is adopted to calculate wave detector anti-pass wave field, adopt effective stored boundary strategy (efficiency frontier storage policy is a kind of way reconstructing source wavefield simultaneously, this way needs the data using storage, reconstruct source wavefield in conventional time domain full waveform inversion and adopt CPU, the present invention adopts GPU.) and many GPU parallel computing reconstruct source wavefield (namely using the border wave field in each moment stored in normal propagation process as boundary condition, all wave fields in last moment stored are as starting condition, then the wave equation of backpropagation is separated, this is the Solve problems of Equations of Mathematical Physics, has separated i.e. restructural wave field.) and (gradient can regard source wavefield as and with the cross-correlation of wave field, after two wave fields solve, gradient is also corresponding to be solved relative to the gradient of model parameter to calculate cost functional.)；

The present invention adopts effective stored boundary strategy to reconstruct source wavefield, (the border number of plies of storage is the half of spatial accuracy difference to the wave field value of the Particle Vibration Velocity component each moment interior zone and PML (perfectly matched layer) region had a common boundary when calculating the source wavefield of forward-propagating, the low order differential thought excessive to higher difference is adopted when reconstructing wave field) and all wave field value of maximum moment are stored in GPU video memory, and (this stores is the GPU internal memory needing application in employing GPU calculates, when after the wave field obtaining the maximum moment, the data needed in GPU internal memory are needed to be delivered in CPU internal memory, the pilot process of the 3rd step can be regarded as), in the maximum moment of forward-propagating, above-mentioned storage data are delivered to CPU internal memory by GPU video memory.In the process of inverse time anti-pass detection residual error wave field, CPU is stored wave field and copies GPU video memory to, the residual error of super big gun collection calculated according to second step based on Huygens' principle to source wavefield approximate reconstruction, reconstruct source wavefield and many GPU parallel computing all can be adopted to accelerate when calculating detection wave field.This process, only need to do two secondary data transmission to the video memory of GPU and the internal memory of CPU, and do not need to store data and read, the reverse-time migration needing data to store compared to tradition and to read and full waveform inversion, only many this process is reconstructed to interior zone source wavefield, this process can adopt and not add CPML boundary condition (equation in front is all adopt based on the ACOUSTIC WAVE EQUATION of CPML (convolution perfectly matched layer-) and equations for elastic waves, to suppress the reflection wave on border, and when reconstructing source wavefield, adopt be the ACOUSTIC WAVE EQUATION and the equations for elastic waves that do not add CPML) ACOUSTIC WAVE EQUATION and Time Migration of Elastic Wave Equation simulation, when adopting many GPU, its counting yield promotes greatly.After reconstruct wave field, can gradient be solved with wave field cross-correlation, optimized algorithm can be adopted after gradient solves to complete rate pattern iteration and upgrade, be i.e. full waveform inversion.

The anti-pass wave field (the adjoint wave field namely in Fig. 1) that the present invention adopts adjoint equation calculating wave detector to go out, to one-order velocity-stress ACOUSTIC WAVE EQUATION, its expression formula is:

L ^Tλ＝g

Wherein: $L^T=A^T{\∂}_x+B^T{\∂}_y+C^T{\∂}_z+D^T{\∂}_t,\mathrm{\λ}={({\mathrm{\λ}}_x,{\mathrm{\λ}}_y,{\mathrm{\λ}}_z,Q)}^T,g={(g_x,g_y,g_z,0)}^T$

To one-order velocity-stress Time Migration of Elastic Wave Equation, the adjoint equation of employing is:

L ^Tλ＝g

Wherein:

$L^T=A^T{\∂}_x+B^T{\∂}_y+C^T{\∂}_z+D^T{\∂}_t,\mathrm{\λ}={({\mathrm{\λ}}_x,{\mathrm{\λ}}_y,{\mathrm{\λ}}_z,Q_{\mathrm{xx}},Q_{\mathrm{yy}},Q_{\mathrm{zz}},Q_{\mathrm{xy}},Q_{\mathrm{xz}},Q_{\mathrm{yz}})}^T,g={(g_x,g_y,g_z,\mathrm{0,0,0,0,0,0})}^T$

The ACOUSTIC WAVE EQUATION cost functional that the present invention adopts about the expression formula of the gradient of model parameter is:

$\mathrm{grad}\left(v_p\right)=-2\mathrm{\ρ}{v}_p\underset{\mathrm{is}=1}{\overset{\mathrm{NS}}{\mathrm{\Σ}}}{\∫}_0^T\underset{x_i}{\mathrm{\Σ}}Q\frac{\∂v_{x_i}}{\∂x_i}\mathrm{dt},x_i=x,y,z$ (sound wave)

The wherein number of super focus that obtains after representing employing three randomizing scheme of NS.

The Time Migration of Elastic Wave Equation cost functional that the present invention adopts about the gradient expression formula of model parameter is:

$\mathrm{grad}\left(\mathrm{\λ}\right)=\underset{\mathrm{is}=1}{\overset{\mathrm{NS}}{\mathrm{\Σ}}}{\∫}_0^T\underset{x_i}{\mathrm{\Σ}}{Q}_{x_i{x}_i}\underset{x_i}{\mathrm{\Σ}}\frac{\∂v_{x_i}}{\∂x_i}\mathrm{dt}$

$\mathrm{grad}\left(\mathrm{\μ}\right)=\underset{\mathrm{is}=1}{\overset{\mathrm{ns}}{\mathrm{\Σ}}}{\∫}_0^T\underset{x_i}{\mathrm{\Σ}}{Q}_{x_i{x}_i}\frac{\∂v_{x_i}}{\∂x_i}+\frac{1}{2}\underset{x_i{y}_i}{\mathrm{\Σ}}{Q}_{x_i{y}_i}(\frac{\∂v_{x_i}}{\∂y_i}+\frac{\∂v_{y_i}}{\∂x_i})\mathrm{dt},x_i=x,y,z;y_i=z,y,z$

5th step, realize the multiple dimensioned full waveform inversion of efficient time territory wave equation based on three randomizing schemes, multi-grid method and Wiener filter technique.

Adopt gradient class iterative algorithm renewal speed model, super big gun record is calculated to the rate pattern after upgrading and and the super big gun record of the first step do residual error, when residual error meets end condition:

time this frequency range full waveform inversion iteration terminate, and the initial model of the model of inverting as next frequency band full waveform inversion, repeat second step to the 5th step, when meeting given frequency end condition namely: complete final multiple dimensioned full waveform inversion flow process during given inverting frequency band number, export the net result of "current" model as full waveform inversion.

When adopting the super focus full waveform inversion based on focus coding, cross-talk noise can be there is in inversion result, and the introducing of super focus makes the nonlinear degree of inverse problem deepen further, for suppressing random cross-talk noise, reduce the non-linearity of inverting, (inverting namely in the 5th step is successively from low-frequency range to the thought of high frequency inverting for the multiple dimensioned time domain full waveform inversion algorithm that the present invention adopts based on multi-grid method.), inverting progressively becomes more meticulous from macrolattice, is fine to objective body gradually from Anomalies of Backgrounds body.

Fig. 2 is for when video card adopts C2060, adopt the speed-up ratio comparison diagram that single GPU card accelerates the Acoustic Forward Modeling under the different mesh scale of two dimension, rhombus broken line represents the computing time adopting CPU simulated sound wave traveling, square broken line represents the computing time adopting GPU simulated sound wave traveling, and triangle broken line represents speed-up ratio.When employing blocks more, the linear increase of speed-up ratio, Fig. 3, for when video card adopts C2060, adopts the speed-up ratio comparison diagram that single GPU card accelerates the Acoustic Forward Modeling under the different mesh scale of three-dimensional.Fig. 4-1 to Fig. 4-4 this be the velocity profile adopting this example to adopt acoustic full waveform inverting to obtain to Marmousi2 model.Fig. 4-1 is Marmousi2 rate pattern accurately, Fig. 4-2 is the initial velocity model after Gauss's large scale is level and smooth, the inversion result that the super big gun collection inverting of Fig. 4-3 to be dominant frequency be 5Hz obtains, can find out in model from result and there is stronger cross-talk noise, using the result of Fig. 4-3 as initial model, and with this, big gun collection is surpassed to high frequency band and carry out full waveform inversion, the super big gun collection inverting gained rate pattern of Fig. 4-4 to be dominant frequency be 15Hz, as can be seen from model, cross-talk noise is better suppressed.Fig. 5-1 to Fig. 5-8 is the Lame's constant sections adopting this example to adopt elastic wave full waveform inversion Simultaneous Inversion to obtain to Overthrust model, wherein Fig. 5-1,5-2 is initial Lame's constant, Fig. 5-3, Fig. 5-4 adopts dominant frequency to be the super big gun collection inverting acquired results of 5.0Hz, using this result as the super big gun collection inverting of initial model to 10.0Hz, Fig. 5-5, Fig. 5-6 is super big gun collection inverting acquired results of 10.0Hz, this result is carried out inverting as the super big gun collection of initial model to 15Hz, Fig. 5-7, Fig. 5-8 is its corresponding inversion result.

Huge calculated amount and memory space are the key factors of restriction 3D Time Migration of Elastic Wave Equation full waveform inversion development always.The problem large for time domain wave equation full waveform inversion calculated amount and memory space is large, the present invention adopts respectively and is solved based on many GPU parallel computing of CUDA and efficiency frontier storage policy, and in conjunction with many focus technology of three random phase encodings, process is accelerated to time domain Time Migration of Elastic Wave Equation full waveform inversion, the accurate enforcement for full waveform inversion provides a set of efficiently feasible technical scheme.

Technique scheme is one embodiment of the present invention, for those skilled in the art, on the basis that the invention discloses application process and principle, be easy to make various types of improvement or distortion, and the method be not limited only to described by the above-mentioned embodiment of the present invention, therefore previously described mode is just preferred, and does not have restrictive meaning.

Claims (5)

1. an efficient time domain full waveform inversion method, is characterized in that: described method comprises:

The first step, adopt that random, the source position of source polarity is random, number random three randomizing schemes in source are randomly formed different frequency bands to original observation seismologic record and surpass big gun collection;

Second step, adopt high-order staggered-mesh finite difference simulation based on the ACOUSTIC WAVE EQUATION system of many GPU parallel computing acceleration time territory single order form of CUDA and Time Migration of Elastic Wave Equation system, obtain GPU and simulate the i-th freq frequency band and surpass big gun collection; Adopt CPU and GPU coordination technique when forward simulation, the focus different with many GPU nuclear mockup is propagated, and the data boundary of source wavefield and last moment all wave fields are stored in calculator memory; Wave field value and all wave field value of maximum moment of the Particle Vibration Velocity component each moment interior zone and PML region had a common boundary when calculating the source wavefield of forward-propagating are stored in GPU video memory;

3rd step, calculating wave field residual error and cost functional value, then judge whether to meet end condition, if, represent that the full waveform inversion iteration of this frequency range terminates, and using the initial velocity model of the model of inverting as next frequency band full waveform inversion, proceed to the 6th step, if not, then proceed to the 4th step;

4th step, adopt many GPU parallel computing to calculate wave detector anti-pass wave field, adopt effective stored boundary strategy and many GPU parallel computing reconstruct source wavefield simultaneously and calculate the gradient of cost functional relative to model parameter;

5th step, employing gradient class iterative algorithm renewal speed model, then return second step;

6th step, judge whether to meet frequency end condition, if so, then proceed to the 7th step, if not, then return the first step;

7th step, output "current" model.

2. efficient time domain full waveform inversion method according to claim 1, is characterized in that: the described first step is achieved in that

Adopt the MPI original single shot record of parallel input and original Ricker wavelet, utilize S filter to original single shot record and original Ricker wavelet filtering, obtain earthquake single shot record and the seismic wavelet of different frequency bands;

Function, random function and the sign function producing random seed in C language built-in function is utilized to give [-1 to all seismic wavelets, 1] number between, namely the color that each seismic wavelet is different is given, according to the positive negativity of seismic wavelet color, give seismic wavelet corresponding polarity, this process is that source polarity is random;

Heap algorithm or the color of quick sort to the seismic wavelet changing source polarity is adopted to sort, color value increases gradually, and adopting the random function in C language to go out the number of single seismic wavelet in each super focus at random according to total super source wavelet number and total single seismic wavelet number, this process is that source number is random;

According to the number of the single seismic wavelet of each super focus, the single seismic wavelet changing source polarity after sequence is superposed, the position of each seismic wavelet now superposed is also at random corresponding, this process is that source position is random, completes the super focus that three randomizing schemes form different frequency bands thus;

In the same way the super big gun collection that three randomizing schemes form different frequency bands is completed to earthquake single shot record.

3. efficient time domain full waveform inversion method according to claim 2, is characterized in that: described second step comprises:

Initial velocity model is obtained by after true model large scale Gaussian smoothing;

Adopt the propagation of high-order staggering mesh finite-difference simulation wave equation in underground medium, specific as follows:

When adopting GPU to calculate, the whole domain of two dimension is considered as a kernel function, the whole domain of two dimension is divided into multiple subregion, the size of every sub regions is 16 × 16 grid nodes, every sub regions is mapped to the thread block of a CUDA, each net point adopts the individual process of thread block independently updated, adopt the parallel thread in GPU to upgrade all nodal values at each time step, thus the wave field completing whole region calculate simultaneously;

For three-dimensional data, three-dimensional data is arranged in two-dimentional sheet according to the direction of tieing up the most slowly, then according to the mode identical with 2-D data process, three-dimensional nodes is mapped to corresponding process, completes the forward simulation in whole region.

4. efficient time domain full waveform inversion method according to claim 3, is characterized in that: the end condition in described 3rd step is as follows:

$\left|\frac{{\mathrm{\δd}}_k-\mathrm{\δ}{d}_{k-1}}{\mathrm{\δ}{d}_{k-1}}\right|\≤\mathrm{\ϵ}$

Wherein, d _krepresent the cost functional value of the wave field residual error of kth time iteration, d _k-1represent the cost functional value of the wave field residual error of kth-1 iteration, ε is a very little amount.

5. efficient time domain full waveform inversion method according to claim 4, is characterized in that: described 4th step comprises:

Adopt the anti-pass wave field that adjoint equation calculating wave detector goes out;

Wave field value and all wave field value of maximum moment of the Particle Vibration Velocity component each moment interior zone and PML region had a common boundary when calculating the source wavefield of forward-propagating are stored in GPU video memory, in the maximum moment of forward-propagating, the storage data in the 3rd step are delivered to CPU internal memory by GPU video memory; In the process of inverse time anti-pass detection residual error wave field, CPU is stored wave field and copies GPU video memory to, the residual error of super big gun collection calculated according to the 3rd step based on Huygens' principle to source wavefield approximate reconstruction, reconstruct source wavefield and all adopt many GPU parallel computing to accelerate when calculating detection wave field.