CN105319581B - A kind of efficient time-domain full waveform inversion method - Google Patents

A kind of efficient time-domain full waveform inversion method Download PDF

Info

Publication number
CN105319581B
CN105319581B CN201410373265.6A CN201410373265A CN105319581B CN 105319581 B CN105319581 B CN 105319581B CN 201410373265 A CN201410373265 A CN 201410373265A CN 105319581 B CN105319581 B CN 105319581B
Authority
CN
China
Prior art keywords
gpu
source
random
wave field
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201410373265.6A
Other languages
Chinese (zh)
Other versions
CN105319581A (en
Inventor
王杰
王立歆
方伍宝
胡光辉
孙晶梅
王振宇
尹力
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201410373265.6A priority Critical patent/CN105319581B/en
Publication of CN105319581A publication Critical patent/CN105319581A/en
Application granted granted Critical
Publication of CN105319581B publication Critical patent/CN105319581B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention provides a kind of efficient time-domain full waveform inversion method, belongs to oil seismic exploration data efficient rate, low storage, high precision velocity modeling field.The first step, using source polarity is random, source position is random, number random three randomizing scheme in source is randomly formed the super big gun collection of different frequency bands to original observation earthquake record；Second step, using more GPU parallel computings acceleration time domain single order forms based on CUDA ACOUSTIC WAVE EQUATION system and Time Migration of Elastic Wave Equation system high-order staggered-mesh finite difference simulate, obtain GPU simulation the super big gun collection of the i-th freq frequency bands；CPU and GPU coordination techniques are used in forward simulation, are propagated with the different focus of more GPU nuclear mockups, the data boundary and last moment all wave fields of source wavefield are stored in calculator memory.

Description

A kind of efficient time-domain full waveform inversion method
Technical field
The invention belongs to oil seismic exploration data efficient rate, low storage, high precision velocity modeling field, and in particular to one The efficient time-domain full waveform inversion method of kind, for huge amount of calculation and mass memory needed for the full wavefield velocity modeling of time-domain Problem, encoded using focus, more GPU combine parallel computing strategy and no data I/O efficiency frontier storage strategy it is complete Into efficient, high-resolution velocity modeling.
Background technology
Using pre-stack seismic wave number according to kinematics and dynamic information, full waveform inversion, which has, to be disclosed under complex geological condition Geological structure, the ability of rock parameter, it is considered as velocity modeling precision highest technology at present.But three-dimensional full-wave shape Also by a variety of limitations, two of which Major Difficulties are that huge amount of calculation and amount of storage do not obtain preferably always for the realization of inverting Ground solves.
Because the correct implementation of full waveform inversion needs hundreds and thousands of secondary forward simulations, thus the efficiency of forward simulation is to restrict Where the essence of full waveform inversion development.To two dimensional inversion problem, using the frequency domain Wave equation forward modeling of direct solution Technology efficiency when handling more focus simulations is higher, but to three-dimensional problem, even if direct using the higher nested subdivision of efficiency Solution, the computation complexity of frequency domain forward modeling are also up to nsnwO(n4Logn), required internal memory O (n5), its amount of calculation and storage Huge, current computer condition is measured to be difficult to realize.To D integral pin-fin tube, the time computation complexity of time finite element method forward modeling For O (n4), required internal memory complexity is O (n3), relative to the forward modeling of three-dimensional frequency domain, time-domain forward modeling has that amount of calculation is small, internal memory Take low advantage.Mainly there are two kinds of thinkings to the Study on Acceleration of forward simulation at present：First, by the acceleration of computer hardware with Improve the computational efficiency (such as the GPU parallel computations based on CUDA platforms) of inverting, two, raising algorithm computational efficiency in itself, More source inversion technical thoughts based on random phase encoding are preferable, and its basic thought is all single source data superpositions Handled into single or multiple super source datas, and to super source data.Due to full waveform inversion computational efficiency with The number approximation of focus is directly proportional, thus the computational efficiency of full waveform inversion can be greatly improved using super focus.
Another problem that Conventional Time domain full waveform inversion faces is to need mass data to input and store, magnanimity number According to input and output (I/O) higher challenge is not only proposed to computer hardware, and had a strong impact on inverting calculating effect Rate.Even if using GPU parallel computings to when needing the Conventional Time domain full waveform inversion of output data to accelerate, improved efficiency has Limit.Moussa (2009) shows that 91% run time shows for CPU internal memories and GPU when being studied using GPUPU reverse-time migration The data transfer deposited, this is due to that the source wavefield at each moment need to be stored on hard disk, caused by limiting GPU parallel efficiency.
The content of the invention
The problem of present invention is directed to amount of calculation needed for time-domain full waveform inversion and huge amount of storage, there is provided a kind of efficient Time-domain full waveform inversion method, a kind of efficient processing means are provided for the high precision velocity modeling of large-scale data.
The present invention is achieved by the following technical solutions：
The first step, using source polarity is random, source position is random, number random three randomizing scheme in source is to original observation earthquake note Record is randomly formed the super big gun collection of different frequency bands；
Second step, the ACOUSTIC WAVE EQUATION system using more GPU parallel computings acceleration time domain single order forms based on CUDA The high-order staggered-mesh finite difference of system and Time Migration of Elastic Wave Equation system is simulated, and the i-th freq frequency bands for obtaining GPU simulations surpass Big gun collection；CPU and GPU coordination techniques are used in forward simulation, are propagated with the different focus of more GPU nuclear mockups, source wavefield Data boundary and last moment all wave fields be stored in calculator memory；The handle when calculating the source wavefield of forward-propagating The wave field value and all wave field values of maximum moment for the Particle Vibration Velocity component that each moment interior zone has a common boundary with PML regions are deposited Storage is in GPU video memorys；
3rd step, calculating wave field residual error and cost functional value (subtract each other to obtain wave field with the super big gun collection and original super big gun collection of simulation Residual error is superimposed cost functional value of summing to obtain), then judge whether to meet end condition, if it is, representing the Full wave shape of the frequency range Inverting iteration terminates, and using the model of inverting as the initial velocity model of next frequency band full waveform inversion, is transferred to the 6th step, If it is not, then it is transferred to the 4th step；
4th step, wave detector anti-pass wave field is calculated using more GPU parallel computings, while use effective stored boundary plan Slightly reconstruct source wavefield with more GPU parallel computings and calculate gradient of the cost functional relative to model parameter；
5th step, using gradient class iterative algorithm renewal speed model, super big gun is calculated to the rate pattern after renewal and remembered Record and do residual error with the super big gun record of the first step；
6th step, judge whether to meet frequency end condition, if it is, the 7th step is transferred to, if it is not, then returning to first Step；
7th step, output "current" model.
What the first step was realized in：
Input original single shot record and original Ricker wavelet parallel using MPI, original single-shot is remembered using Wiener filter Record and the filtering of original Ricker wavelet, obtain the earthquake single shot record and seismic wavelet of different frequency bands；
It is sub to all earthquakes using function, random function and sign function that random seed is produced in C language built-in function Ripple assigns a number between [- 1,1], that is, assigns each seismic wavelet different colors, according to the positive and negative of seismic wavelet color Property, the corresponding polarity of seismic wavelet is assigned, this process is that source polarity is random；
The color for changing the seismic wavelet of source polarity is sorted using Heap algorithm or quick sort, color value is gradual Increase, and gone out at random according to total super source wavelet number and total single seismic wavelet number using the random function in C language The number of single seismic wavelet in each super focus, this process are that source number is random；
According to each super focus single seismic wavelet number to after sequence change source polarity single seismic wavelet Superposition, the position for each seismic wavelet being now superimposed is also corresponding random, and this process is that source position is random, and it is random thus to complete three Scheme forms the super focus of different frequency bands；
The super big gun collection that three randomizing schemes form different frequency bands is completed to earthquake single shot record in the same way.
Used Wiener filter is：
FWienerIt is Wiener filter, WOriginalIt is primary signal, WTargetFor desired signal, ω is angular frequency, ε Choose 10-6。
The second step includes：
By obtaining initial velocity model after true model large scale Gaussian smoothing；
It is specific as follows using propagation of the high-order staggering mesh finite-difference simulation wave equation in underground medium：
The whole domain of two dimension is considered as a kernel function when being calculated using GPU, the whole domain of two dimension More sub-regions are divided into, be 16 × 16 grid nodes per the sizes of sub-regions, often sub-regions are mapped to CUDA's Thread block, each mesh point is independently updated using the individual process of thread block, the parallel thread in each time step uses GPU All nodal values are updated simultaneously, are calculated so as to complete the wave field of whole region；
For three-dimensional data, three-dimensional data is arranged in two-dimensional slice according to the direction most tieed up slowly, then according to two-dimemsional number According to processing identical mode, three-dimensional nodes are mapped to corresponding process, complete the forward simulation of whole region.
End condition in 3rd step is as follows：
Wherein, dkRepresent the cost functional value of the wave field residual error of kth time iteration, dk-1Represent that the wave field of -1 iteration of kth is residual The cost functional value of difference, ε are the amount of a very little.
4th step includes：
The anti-pass wave field gone out using adjoint equation calculating wave detector；
The particle vibration speed that each moment interior zone and PML regions are had a common boundary when calculating the source wavefield of forward-propagating Spend component wave field value and all wave field values of maximum moment be stored in GPU video memorys, forward-propagating the maximum moment the 3rd Data storage in step is delivered to CPU internal memories by GPU video memorys；During inverse time anti-pass detection residual error wave field, CPU is stored Wave field copies GPU video memorys to, and the residual error of the super big gun collection calculated according to the 3rd step is based on Huygen's principle to the approximate weight of source wavefield Accelerated when structure, reconstruct source wavefield and calculating detection wave field using more GPU parallel computings.
Compared with prior art, the beneficial effects of the invention are as follows：Fluctuated using more GPU parallel computings acceleration time domains The numerical simulation of equation system, the simulation number of focus is reduced using the focus coding techniques of three randomizing schemes, makes time-domain complete The efficiency of waveform inversion greatly promotes, and using GPU calculate caused by the data input output of full waveform inversion for tradition The low problem of efficiency, the present invention use effective stored boundary strategy, the input and output of data are at utmost reduced, so that entirely The computational efficiency of waveform inversion is further lifted, and the final correct implementation for time-domain full waveform inversion provides a kind of height Imitate feasible technical scheme.
Brief description of the drawings
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：The acceleration effect that monokaryon GPU carries out time finite element method forward simulation to three dimensional elasticity wave equation contrasts Figure
Fig. 4-1：True Marmousi2 models
Fig. 4-2：Initial velocity of longitudinal wave model
Fig. 4-3：ACOUSTIC WAVE EQUATION 5.0Hz Multi-scale inversion results, cross-talk noise are serious
Fig. 4-4：ACOUSTIC WAVE EQUATION 15.0Hz Multi-scale inversion result cross-talk noises are suppressed
Fig. 5-1：Initial Overthurst models (Lame constants lambda)
Fig. 5-2：Initial Overthurst models (Lame constants mu)
Fig. 5-3：Time Migration of Elastic Wave Equation 5.0Hz Multi-scale inversion results (Lame constants lambda)
Fig. 5-4：Time Migration of Elastic Wave Equation 5.0Hz Multi-scale inversion results (Lame constants mu)
Fig. 5-5：Time Migration of Elastic Wave Equation 10.0Hz Multi-scale inversion results (Lame constants lambda)
Fig. 5-6：Time Migration of Elastic Wave Equation 10.0Hz Multi-scale inversion results (Lame constants mu)
Fig. 5-7：Time Migration of Elastic Wave Equation 15.0Hz Multi-scale inversion results (Lame constants lambda)
Fig. 5-8：Time Migration of Elastic Wave Equation 15.0Hz Multi-scale inversion results (Lame constants mu)
Embodiment
The present invention is described in further detail below in conjunction with the accompanying drawings：
The present invention is a kind of based on focus phase code and more GPU parallel computings and no data I/O efficient time domain Wave equation full waveform inversion method, a kind of efficient processing means are provided for the high precision velocity modeling of large-scale data.
As shown in figure 1, the inventive method includes：
The first step, using source polarity is random, source position is random, number random three randomizing scheme in source is to original observation earthquake note Record is randomly formed the super big gun collection of different frequency bands；
It is specific as follows：
Input original single shot record and original Ricker wavelet parallel using MPI, original single-shot is remembered using Wiener filter Record and the filtering of original Ricker wavelet, obtain the earthquake single shot record and seismic wavelet of different frequency bands, used Wiener filter For：
Wherein FWienerIt is Wiener filter, WOriginalIt is primary signal, WTargetFor desired signal, ω is angular frequency Rate, ε are the values of a very little, and in order to avoid the generation of wild effect, the present invention chooses 10-6.Using being produced in C language built-in function Srand () function, random function rand () and the sign function of raw random seed all seismic wavelets are assigned [- 1,1] it Between a number, that is, assign each seismic wavelet different colors, according to the positive negativity of seismic wavelet color, assign seismic wavelet Corresponding polarity, it is random now to complete source polarity；
The color for changing the seismic wavelet of source polarity is sorted using Heap algorithm or quick sort, color value is gradual Increase, and using the random function rand () in C language according to total super source wavelet number and total single seismic wavelet number Go out the number of single seismic wavelet in each super focus at random, this process is that source number is random；
Hereafter according to each super focus single seismic wavelet number to after sequence change source polarity single earthquake Wavelet is superimposed, and the position for each seismic wavelet being now superimposed is also corresponding random, and this process is that source position is random, thus completes three Randomizing scheme forms the super focus of different frequency bands, and the super big gun of different frequency bands is formed to earthquake single shot record with same implementation Collection.
Second step, using more GPU parallel computings acceleration time domain single order forms based on CUDA ACOUSTIC WAVE EQUATION and The high-order staggered-mesh finite difference simulation of Time Migration of Elastic Wave Equation system.
Second step is the artificial synthesized super big gun collection of forward simulation and the super big gun collection obtained using the first step and the super big gun collection simulated Residual error is done, wherein the wavelet used is super source wavelet, initial model after true model large scale Gaussian smoothing by obtaining.
Forward problem is the basis of inversion problem, the ACOUSTIC WAVE EQUATION system of single order form that the present invention uses for：
It is write as operational form：
Lv=f
Wherein：
The Time Migration of Elastic Wave Equation system for the single order form that the present invention uses for：
It is write as operational form：
Lv=f
Wherein：
Propagation of the present invention using high-order staggering mesh finite-difference simulation wave equation in underground medium, outside wave field When pushing away, the wave field of each grid node used in renewal the node previous moment wave field value and surroundings nodes before half when The wave field value at quarter, the problem have good concurrency, when being calculated using GPU the whole domain (first step of two dimension Super big gun collection is obtained, second step is to simulate super big gun collection with GPU, it is therefore an objective to seeks both residual errors) it is considered as a kernel function, two dimension Whole domain is divided into more sub-regions, and the size per sub-regions is 16 × 16 grid nodes, is mapped per sub-regions To CUDA thread block, each mesh point is independently updated using the individual process of thread block, and GPU is used in each time step In parallel thread simultaneously update all nodal values, so as to complete the wave field of whole region calculate.For three-dimensional problem, three-dimensional Data are arranged in two-dimensional slice according to the direction (i.e. the y-axis direction of cartesian coordinate system) most tieed up slowly, the processing based on two-dimensional problems Model split process block and process, three-dimensional nodes are mapped to corresponding process, complete the forward simulation of whole region.Due to each The numerical simulation of individual focus (the corresponding node of a focus) is separate, thus CPU/ can be used in forward simulation GPU coordination techniques, propagated with the different focus of more GPU nuclear mockups.
3rd step, using (second step) during GPU speed-up computation source wavefields the data boundary of source wavefield and last One moment all wave field is stored in calculator memory (i.e. " border wave field is deposited in forward modeling (CPU) " in Fig. 1)；
4th step, wave detector anti-pass wave field is calculated using more GPU parallel computings, while use effective stored boundary plan Slightly (efficiency frontier storage strategy is a kind of method for reconstructing source wavefield, and the method needs to use the data of storage, when conventional Between in the full waveform inversion of domain reconstruct source wavefield use CPU, it is of the invention to use GPU.) and the reconstruct shake of more GPU parallel computings Source wave field (i.e. using the border wave field at each moment stored during normal propagation as boundary condition, storage last Then all wave fields at moment solve the wave equation of backpropagation as primary condition, this is asking for Equations of Mathematical Physics Solution problem, i.e. restructural wave field is solved.) and calculate cost functional (gradient can regard focus as relative to the gradient of model parameter Wave field and the cross-correlation with wave field, after two wave fields solve, gradient also can accordingly solve.)；
The present invention is using effective stored boundary strategy reconstruct source wavefield, when calculating the source wavefield of forward-propagating each Wave field value (the border of storage for the Particle Vibration Velocity component that individual moment interior zone has a common boundary with PML (perfectly matched layer) region The number of plies is the half of spatial accuracy difference, when reconstructing wave field using low order differential to the excessive thought of higher difference) and Maximum moment all wave field values are stored in GPU video memorys that (this storage is to be needed when calculating using GPU in the GPU that applies Deposit, after the wave field at maximum moment is obtained, it is necessary to the data transfer needed in GPU internal memories into CPU internal memories, can regard as The pilot process of 3rd step), above-mentioned data storage is delivered to CPU internal memories by GPU video memorys at the maximum moment of forward-propagating. During inverse time anti-pass detection residual error wave field, CPU storage wave fields are copied to GPU video memorys, the super big gun calculated according to second step The residual error of collection is based on Huygen's principle to can be used when source wavefield approximate reconstruction, reconstruct source wavefield and calculating detection wave field More GPU parallel computings accelerate.The process, it is only necessary to the internal memory of video memory and CPU to GPU does data transfer twice, and not Need to data storage and reading, need data storage and the reverse-time migration read and full waveform inversion compared to tradition, it is only more This process is reconstructed to interior zone source wavefield, the process, which can use, to be not added with CPML boundary conditions (equation in front is all Using ACOUSTIC WAVE EQUATION and equations for elastic waves based on CPML (convolution perfectly matched layer -), to suppress the back wave on border, and reconstruct During source wavefield, using the ACOUSTIC WAVE EQUATION and equations for elastic waves for being not added with CPML) ACOUSTIC WAVE EQUATION and Time Migration of Elastic Wave Equation mould Intend, its computational efficiency greatly promotes when using more GPU.After reconstructing wave field, with gradient can be solved with wave field cross-correlation, ladder Degree can use optimized algorithm to complete the renewal of rate pattern iteration, i.e. full waveform inversion after solving.
The anti-pass wave field (i.e. adjoint wave field in Fig. 1) that the present invention is gone out using adjoint equation calculating wave detector, to single order speed Degree-stress ACOUSTIC WAVE EQUATION, its expression formula are：
LTλ=g
Wherein：
To one-order velocity-stress Time Migration of Elastic Wave Equation, the adjoint equation used for：
LTλ=g
Wherein：
The ACOUSTIC WAVE EQUATION cost functional that the present invention uses on the gradient of model parameter expression formula for：
(sound wave)
Wherein NS represents to use the number of the super focus obtained after three randomizing schemes.
The Time Migration of Elastic Wave Equation cost functional that the present invention uses on model parameter pressure gradient expression formula for：
5th step, efficient time domain fluctuation side is realized based on three randomizing schemes, multi-grid method and Wiener filter technique The multiple dimensioned full waveform inversion of journey.
Using gradient class iterative algorithm renewal speed model, super big gun record is calculated the rate pattern after renewal and with the The super big gun record of one step does residual error, when residual error meets end condition：
When the frequency range full waveform inversion iteration terminate, and the model of inverting as next frequency Initial model with full waveform inversion, second step is repeated to the 5th step, when the frequency end condition for meeting to give is：Given is anti- Final multiple dimensioned full waveform inversion flow, most termination of the output "current" model as full waveform inversion are completed when drilling frequency band number Fruit.
During using the super focus full waveform inversion encoded based on focus, there can be cross-talk noise in inversion result, and The introducing of super focus makes the nonlinear degree of inversion problem further deepen, and to suppress random cross-talk noise, reduces inverting Non-linear property, the present invention are (i.e. anti-in the 5th step using the multiple dimensioned time-domain full waveform inversion algorithm based on multi-grid method Drill successively from low-frequency range to the thought of high frequency inverting.), inverting progressively becomes more meticulous from big grid, gradually fine from Anomalies of Backgrounds body To objective body.
Fig. 2 is when video card uses C2060, using single GPU card to the sound wave forward modeling mould under the different mesh scales of two dimension Intend the speed-up ratio comparison diagram accelerated, rhombus broken line represents the calculating time using CPU simulation Acoustic Wave Propagations, square folding Line represents the calculating time using GPU simulation Acoustic Wave Propagations, and triangle broken line represents speed-up ratio.When using block more, speed-up ratio into Linear relationship increase, Fig. 3 be when video card use C2060, using single GPU card to the sound wave under the different mesh scales of three-dimensional just Drill the speed-up ratio comparison diagram that simulation is accelerated.Fig. 4-1 to Fig. 4-4 this to be adopted using the example to Marmousi2 models The velocity profile obtained with acoustic full waveform inverting.Fig. 4-1 is accurate Marmousi2 rate patterns, and Fig. 4-2 is by Gauss Initial velocity model after large scale is smooth, Fig. 4-3 is the inversion result that the super big gun collection inverting that dominant frequency is 5Hz obtains, from result In it can be seen that model in there is stronger cross-talk noise, using Fig. 4-3 result as initial model, and with this to high frequency band Super big gun collection carries out full waveform inversion, and Fig. 4-4 is the super big gun collection inverting gained rate pattern that dominant frequency is 15Hz, can from model Go out, cross-talk noise is preferably suppressed.Fig. 5-1 to Fig. 5-8 is complete using elastic wave to Overthrust models using the example The Lame constants section that waveform inversion Simultaneous Inversion obtains, wherein Fig. 5-1,5-2 are initial Lame constants, Fig. 5-3, Fig. 5-4 It is to use super big gun collection inverting acquired results of the dominant frequency for 5.0Hz, the result is anti-to 10.0Hz super big gun collection as initial model Drill, Fig. 5-5, Fig. 5-6 are 10.0Hz super big gun collection inverting acquired results, using super big gun collection of the result as initial model to 15Hz Inverting is carried out, Fig. 5-7, Fig. 5-8 are its corresponding inversion result.
Huge amount of calculation and amount of storage are always an important factor for restricting the development of 3D Time Migration of Elastic Wave Equation full waveform inversion. For the problem of time-domain wave equation full waveform inversion is computationally intensive and amount of storage is big, the present invention is respectively adopted based on CUDA's More GPU parallel computings and efficiency frontier storage strategy are solved, and combine more focus technologies of three random phase encodings The acceleration of time-domain Time Migration of Elastic Wave Equation full waveform inversion is handled, provides and a set of efficiently may be used for the accurate implementation of full waveform inversion Capable technical scheme.
Above-mentioned technical proposal is one embodiment of the present invention, for those skilled in the art, at this On the basis of disclosure of the invention application process and principle, it is easy to make various types of improvement or deformation, be not limited solely to this Invent the method described by above-mentioned embodiment, therefore previously described mode is simply preferable, and and without limitation The meaning of property.

Claims (5)

  1. A kind of 1. efficient time-domain full waveform inversion method, it is characterised in that：Methods described includes：
    The first step, using source polarity is random, source position is random, random three randomizing scheme of source number to original observation earthquake record with Machine forms the super big gun collection of different frequency bands；
    Second step, using more GPU parallel computings acceleration time domain single order forms based on CUDA ACOUSTIC WAVE EQUATION system and The high-order staggered-mesh finite difference simulation of Time Migration of Elastic Wave Equation system, obtain the GPU simulation super big gun collection of the i-th freq frequency bands； CPU and GPU coordination techniques are used in forward simulation, are propagated with the different focus of more GPU nuclear mockups, the border of source wavefield Data and last moment all wave fields are stored in calculator memory；Calculate forward-propagating source wavefield when it is each when The wave field value and all wave field values of maximum moment for carving the Particle Vibration Velocity component that interior zone has a common boundary with PML regions are stored in In GPU video memorys；
    3rd step, wave field residual error and cost functional value are calculated, then judge whether to meet end condition, if it is, representing the frequency The full waveform inversion iteration of band terminates, and the initial velocity model using the model of inverting as next frequency band full waveform inversion, The 6th step is transferred to, if it is not, then being transferred to the 4th step；
    4th step, calculate wave detector anti-pass wave field using more GPU parallel computings, at the same using effective stored boundary strategy and More GPU parallel computings reconstruct source wavefield and calculate gradient of the cost functional relative to model parameter；
    5th step, using gradient class iterative algorithm renewal speed model, be then back to second step；
    6th step, judge whether to meet frequency end condition, if it is, the 7th step is transferred to, if it is not, then returning to the first step；
    7th step, output "current" model.
  2. 2. efficient time-domain full waveform inversion method according to claim 1, it is characterised in that：The first step is this What sample was realized：
    Input original single shot record and original Ricker wavelet parallel using MPI, using Wiener filter to original single shot record and Original Ricker wavelet filtering, obtains the earthquake single shot record and seismic wavelet of different frequency bands；
    All seismic wavelets are assigned using function, random function and sign function that random seed is produced in C language built-in function A number between [- 1,1] is given, that is, assigns each seismic wavelet different colors, according to the positive negativity of seismic wavelet color, is assigned The corresponding polarity of seismic wavelet is given, this process is that source polarity is random；
    The color for changing the seismic wavelet of source polarity is sorted using Heap algorithm or quick sort, color value gradually increases Greatly, and using the random function in C language according to total super source wavelet number and total single seismic wavelet number go out at random often The number of single seismic wavelet in individual super focus, this process are that source number is random；
    Number according to the single seismic wavelet of each super focus is superimposed to the single seismic wavelet for changing source polarity after sequence, The position for each seismic wavelet being now superimposed is also corresponding random, and this process is that source position is random, thus completes three randomizing schemes Form the super focus of different frequency bands；
    The super big gun collection that three randomizing schemes form different frequency bands is completed to earthquake single shot record in the same way.
  3. 3. efficient time-domain full waveform inversion method according to claim 2, it is characterised in that：The second step bag Include：
    By obtaining initial velocity model after true model large scale Gaussian smoothing；
    It is specific as follows using propagation of the high-order staggering mesh finite-difference simulation wave equation in underground medium：
    The whole domain of two dimension is considered as a kernel function when being calculated using GPU, the whole domain of two dimension is split Into more sub-regions, the size per sub-regions is 16 × 16 grid nodes, and CUDA thread is mapped to per sub-regions Block, each mesh point is independently updated using the individual process of thread block, and the parallel thread in each time step uses GPU is simultaneously All nodal values are updated, are calculated so as to complete the wave field of whole region；
    For three-dimensional data, three-dimensional data is arranged in two-dimensional slice according to the direction most tieed up slowly, then according to at 2-D data Identical mode is managed, three-dimensional nodes are mapped to corresponding process, complete the forward simulation of whole region.
  4. 4. efficient time-domain full waveform inversion method according to claim 3, it is characterised in that：In 3rd step End condition is as follows：
    <mrow> <mo>|</mo> <mfrac> <mrow> <msub> <mi>&amp;delta;d</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>&amp;delta;d</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> <mrow> <msub> <mi>&amp;delta;d</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </mfrac> <mo>|</mo> <mo>&amp;le;</mo> <mi>&amp;epsiv;</mi> </mrow>
    Wherein, dkRepresent the cost functional value of the wave field residual error of kth time iteration, dk-1Represent the wave field residual error of -1 iteration of kth Cost functional value, δ dkRepresent that the record of the same frequency band that the record of kth time iterative modeling and the first step obtain in the frequency band is poor, δ dk-1Represent that the record of the same frequency band that the record of -1 iterative modeling of kth and the first step obtain in the frequency band is poor, ε be one very Small amount.
  5. 5. efficient time-domain full waveform inversion method according to claim 4, it is characterised in that：The 4th step bag Include：
    The anti-pass wave field gone out using adjoint equation calculating wave detector；
    The Particle Vibration Velocity point that each moment interior zone and PML regions are had a common boundary when calculating the source wavefield of forward-propagating The wave field value of amount and all wave field values of maximum moment are stored in GPU video memorys, forward-propagating the maximum moment in the 3rd step Data storage CPU internal memories are delivered to by GPU video memorys；During inverse time anti-pass detection residual error wave field, CPU is stored wave field Copy GPU video memorys to, the residual error of the super big gun collection calculated according to the 3rd step is based on Huygen's principle to source wavefield approximate reconstruction, weight Accelerated when structure source wavefield and calculating detection wave field using more GPU parallel computings.
CN201410373265.6A 2014-07-31 2014-07-31 A kind of efficient time-domain full waveform inversion method Active CN105319581B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410373265.6A CN105319581B (en) 2014-07-31 2014-07-31 A kind of efficient time-domain full waveform inversion method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410373265.6A CN105319581B (en) 2014-07-31 2014-07-31 A kind of efficient time-domain full waveform inversion method

Publications (2)

Publication Number Publication Date
CN105319581A CN105319581A (en) 2016-02-10
CN105319581B true CN105319581B (en) 2018-01-16

Family

ID=55247384

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410373265.6A Active CN105319581B (en) 2014-07-31 2014-07-31 A kind of efficient time-domain full waveform inversion method

Country Status (1)

Country Link
CN (1) CN105319581B (en)

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105842730B (en) * 2016-03-16 2017-03-15 中国海洋大学 A kind of record residual difference acquiring method based on local normalization least square
CN107656306B (en) * 2016-07-26 2019-02-19 中国石油化工股份有限公司 Full waveform inversion parallel calculating method and system
CN106908835B (en) * 2017-03-01 2018-06-08 吉林大学 Band limit Green's function filters multiple dimensioned full waveform inversion method
CN106950596B (en) * 2017-04-11 2019-02-26 中国石油大学(华东) A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate
CN107422379B (en) * 2017-07-27 2019-01-11 中国海洋石油集团有限公司 Multiple dimensioned seismic full-field shape inversion method based on local auto-adaptive convexification method
CN109507726A (en) * 2017-09-15 2019-03-22 中国石油化工股份有限公司 The inversion method and system of time-domain elastic wave multi-parameter Full wave shape
CN108594299B (en) * 2018-02-28 2019-12-31 中国科学院地质与地球物理研究所 Intelligent early warning method, device and system for high-speed rail
CN109541691B (en) * 2018-12-01 2020-05-19 中国石油大学(华东) Seismic velocity inversion method
CN112099936A (en) * 2019-06-17 2020-12-18 中国石油天然气集团有限公司 Heterogeneous parallel computing implementation method and device for three-dimensional acoustic wave NPML algorithm
CN110441816B (en) * 2019-09-20 2020-06-02 中国科学院测量与地球物理研究所 Non-crosstalk multi-seismic-source full-waveform inversion method and device independent of wavelets
CN112578431B (en) * 2019-09-27 2024-04-09 中国石油化工股份有限公司 Method and system for storing full waveform inversion wave field optimization in finite state
CN111208568B (en) * 2020-01-16 2021-04-09 中国科学院地质与地球物理研究所 Time domain multi-scale full waveform inversion method and system
US11281825B2 (en) * 2020-06-30 2022-03-22 China Petroleum & Chemical Corporation Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models
CN114460640B (en) * 2020-11-09 2024-06-25 中国石油天然气集团有限公司 Full waveform inversion method and device for finite difference analog elastic wave
CN112765802B (en) * 2021-01-13 2022-11-29 陕西师范大学 Method for evolving water wave waveform based on high-order water wave model
CN117849879A (en) * 2023-12-20 2024-04-09 京全品质能源技术(北京)有限公司 Wave equation inversion method and device, electronic equipment and storage medium

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2481270A (en) * 2011-01-26 2011-12-21 Nikhil Shah Method of, and apparatus for, full waveform inversion modelling
CN103091711A (en) * 2013-01-24 2013-05-08 中国石油天然气集团公司 Method and device for full-wave-shape inversion
CN103135132A (en) * 2013-01-15 2013-06-05 中国科学院地质与地球物理研究所 Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing
CN103207409A (en) * 2013-04-17 2013-07-17 中国海洋石油总公司 Frequency domain full-waveform inversion seismic velocity modeling method

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101182838B1 (en) * 2010-08-24 2012-09-14 서울대학교산학협력단 Method and Apparatus for Frequency domain Reverse Time Migration with Source Estimation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2481270A (en) * 2011-01-26 2011-12-21 Nikhil Shah Method of, and apparatus for, full waveform inversion modelling
CN103135132A (en) * 2013-01-15 2013-06-05 中国科学院地质与地球物理研究所 Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing
CN103091711A (en) * 2013-01-24 2013-05-08 中国石油天然气集团公司 Method and device for full-wave-shape inversion
CN103207409A (en) * 2013-04-17 2013-07-17 中国海洋石油总公司 Frequency domain full-waveform inversion seismic velocity modeling method

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A new scheme for elastic full waveform inversion based on velocity-stress wave equations in time domain;Jie Wang et al.;《SEG Las Vegas 2012 Annual Meeting》;20121231;第1-5页 *
A strategy for source wavefield reconstruction in reverse time migration;Bo-Feng et al.;《SEG San Antonio 2011 Annual Meeting》;20111231;第3164-3168页 *
全波形反演研究现状及发展趋势;杨勤勇等;《石油物探》;20140131;第53卷(第1期);第77-83页 *
地震叠前逆时偏移的有效边界存储策略;王保利等;《地球物理学报》;20120731;第55卷(第7期);第2412-2421页 *
基于L-BFGS算法和同时激发震源的频率多尺度全波形反演;张生强等;《吉林大学学报（地球科学版）》;20130531;第43卷(第3期);第1004-1012页 *
时间域声波全波形反演及GPU加速;苏超等;《中国地球物理2012》;20121231;第486页 *

Also Published As

Publication number Publication date
CN105319581A (en) 2016-02-10

Similar Documents

Publication Publication Date Title
CN105319581B (en) A kind of efficient time-domain full waveform inversion method
CN108710153B (en) Wave number domain method for magnetic full tensor gradient inversion underground three-dimensional magnetic distribution
Métivier et al. A graph space optimal transport distance as a generalization of Lp distances: application to a seismic imaging inverse problem
CN103713315B (en) A kind of seismic anisotropy parameter full waveform inversion method and device
CN108037531B (en) A kind of seismic inversion method and system based on the full variational regularization of broad sense
CN107783190A (en) A kind of least square reverse-time migration gradient updating method
CN107894618B (en) A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm
CN105137486A (en) Elastic wave reverse-time migration imaging method and apparatus in anisotropic media
WO2023087451A1 (en) Observation data self-encoding-based multi-scale unsupervised seismic wave velocity inversion method
CN106526674A (en) Three-dimensional full waveform inversion energy weighted gradient preprocessing method
CN107765308B (en) Reconstruct low-frequency data frequency domain full waveform inversion method based on convolution thought Yu accurate focus
CN104965222B (en) Three-dimensional longitudinal wave impedance full-waveform inversion method and device
Golubev et al. Compact grid-characteristic scheme for the acoustic system with the piece-wise constant coefficients
CN111580163B (en) Full waveform inversion method and system based on non-monotonic search technology
Li et al. Solving seismic wave equations on variable velocity models with Fourier neural operator
CN105911584A (en) Implicit staggered-grid finite difference elastic wave numerical simulation method and device
Li et al. An immersed boundary method with iterative symmetric interpolation for irregular surface topography in seismic wavefield modelling
Gao et al. An Efficient Multiscale Finite‐Element Method for Frequency‐Domain Seismic Wave Propagation
CN109239776B (en) Seismic wave propagation forward modeling method and device
Zou et al. Deep neural Helmholtz operators for 3D elastic wave propagation and inversion
Wang et al. Acoustic reverse time migration and perfectly matched layer in boundary-conforming grids by elliptic method
Zhou et al. An iterative factored topography-dependent eikonal solver for anisotropic media
CN105974471B (en) A kind of quick forward modelling method of the more GPU of seismic data based on asynchronous flow
Wang et al. A fourth order accuracy summation-by-parts finite difference scheme for acoustic reverse time migration in boundary-conforming grids
CN115373022B (en) Amplitude phase correction-based elastic wave field Helmholtz decomposition method

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant