CN103135132B - Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing - Google Patents
Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing Download PDFInfo
- Publication number
- CN103135132B CN103135132B CN201310013786.6A CN201310013786A CN103135132B CN 103135132 B CN103135132 B CN 103135132B CN 201310013786 A CN201310013786 A CN 201310013786A CN 103135132 B CN103135132 B CN 103135132B
- Authority
- CN
- China
- Prior art keywords
- frequency
- wave field
- gpu
- field
- inverting
- 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.)
- Expired - Fee Related
Links
Landscapes
- Complex Calculations (AREA)
Abstract
The invention discloses a hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing. Compared with a traditional method, the method can be adopted to conduct the CPU/GPU synergetic parallel computing so as to obviously improve computing efficiency. The method enables a forward part of full wave form inversion to be placed to a time domain, namely forward is conducted in the time domain, the forward is converted to be conducted as the inversion in a frequency domain by the discrete fourier transform (DFT) being utilized, namely the DFT is adopted to pick wave field components corresponding to inversion frequency, and the inversion is conducted in the frequency domain from low frequency to high frequency. The method effectively resolves the problem of astringency of a standard time domain method, and avoids the problem that internal memory of wave form inversion of a standard three-dimensional frequency domain and a Laplace domain cannot be satisfied. The method is few in step, reduces computing expenses, and due to the fact the method can be adopted to conduct the CPU/GPU synergetic parallel computing, largely improves speed-up ratio of the computing method.
Description
Technical field
The present invention relates to a kind of full waveform inversion method, particularly relate to and work in coordination with by CPU/GPU the full waveform inversion method that parallel computation performs in hybrid domain, belong to computer realm.
Background technology
Current China difficulty that the is western and exploitation of overseas original oil zone is increasing, and the degree of depth is constantly deepened, to high-precision geophysical imaging and inversion method demand very large; And high-precision formation method finally also will rely on the structure of high precision velocity model.The inversion method when inversion method of rate pattern is divided into walking according to the difference that inverting information uses, amplitude inversion method and full waveform inversion three class methods.Due to dynamics and the kinematics information of full waveform inversion method comprehensive utilization pre-stack seismic wave field, subsurface velocities field can be rebuild accurately, become the important directions of exploration geophysics research both at home and abroad at present.
Existing full waveform inversion (Full Waveform Inversion, be called for short FWI) method comprise Tarantola 1984 propose the time domain full waveform inversion based on Generalized Least Square inversion theory, the people such as Pratt, Shin has developed the waveform inversion of frequency field and Laplace domain on its basis respectively subsequently.Full waveform inversion directly utilizes seismologic record inverting underground medium parameter, utilize the shape informations such as reflection, refraction fully, under certain condition, can Simultaneous Inversion velocity of longitudinal wave, density and shear wave velocity relative change, thus obtain the higher rate pattern of precision, be the important directions of exploration geophysics research both at home and abroad at present.Although full waveform inversion has many technical advantages, due to full waveform inversion calculated amount huge (being about the 30-100 of reverse-time migration doubly), be difficult to realize 3D processing, cause this method cannot adapt to industrial needs.
In existing research, full waveform inversion method is divided three classes from account form: the first kind is time domain waveform inversion (Tarantola 1984,1986), the method utilizes Whole frequency band data compute gradient field, extreme due to objective function is non-linear causes it to be easily absorbed in local extremum, thus inverting is not restrained; Equations of The Second Kind is frequency field waveform inversion (Pratt 1998,1999), the method is done in frequency field and is just drilled, carry out by (group) inverting frequently from low to high in frequency field afterwards, relative high frequency data, low-frequency data and rate pattern have better linear relationship, and the method effectively can avoid the local minimum problem in time domain refutation process; The method just drilled in frequency field is divided into direct method (as: LU decomposes) and process of iteration two kinds, because process of iteration efficiency is low, frequency field is just being drilled the direct methods such as LU decomposition that adopt more and is being solved, the advantage of this method is after once decomposing impedance matrix, result after decomposition may be used for the calculating of all shot points, counting yield is high, this kind of method also also exists the drawback being difficult to overcome, under three-dimensional case, LU decompose EMS memory occupation amount and calculated amount all very huge, cause three-dimensional waveform inverting to implement unrealistic.3rd class is Laplace domain waveform inversion (Shin2008,2009), the method is on the basis of frequency field waveform inversion, decay factor is added to inverted parameters, so inverting is not only from low to high, also from shallow-layer to deep layer, thus reduce the dependence of inversion result to initial model, but the method also faces the memory bottleneck the same with frequency field method, namely just drilling part at it needs LU to decompose, and causes three-dimensional computations to realize difficulty.
At present at home, along with the increase of exploration complexity, many experts and scholars have also carried out the Study and appliance work of waveform inversion, but are generally limited to two dimension, are difficult to use in three-dimensional.In general, waveform inversion technology was the hot issue of academia's research especially in recent years from being born from it always; What the research of waveform inversion was carried out in the world is more, and the research of especially external industry member links is all carried out in anxiety.On the contrary, the research application of domestic waveform inversion is in the starting stage, relatively backwardness substantially, therefore, studies practical waveform inversion method and has great importance to the technological gap filling up China.
The US Patent No. 11/756 of 2007, 384, " System and method for 3D frequencydomain waveform inversion based on 3D time-domain forward modeling ", in technical scheme disclosed in this invention, provide a kind of waveform inversion method: first utilize DFT to be transformed into frequency field big gun numeric field data, then time domain propagation operator and DFT calculated rate territory focus main story wave field is utilized, residual error wave field is calculated in frequency field, and utilize inverse discrete Fourier transform (Inverse DiscreteFourier Transform further, be called for short IDFT) ask for the residual error wave field of time domain, afterwards, with time domain propagation operator, anti-pass is carried out to residual error wave field, and utilize DFT to extract frequency field residual error anti-pass wave field equally, finally, ask for the gradient fields of current inverting frequency.
Compared with the numerical procedure adopted in early days with this area, although this technical scheme achieves certain progress, its process step relative complex, execution efficiency is under some influence; And the parallel granularity of module consuming time (as DFT, time domain propagation operator etc.) main in hybrid domain FWI is very thin, only adopts CPU computing, fail farthest to make program parallelization, counting yield can be greatly affected.Refutation strategy of the present invention is different from it, and the basis of this inverting flow process has been lacked a DFT and time IDFT; And further, adopt GPU/CPU to work in coordination with concurrent operation, farthest make program parallelization, significantly increase the speed-up ratio relative to CPU program.
Summary of the invention
The invention discloses the hybrid domain full waveform inversion method that a kind of CPU/GPU works in coordination with parallel computation, the waveform inversion data processing of two and three dimensions can be carried out.This method combines the advantage of time domain waveform inversion and frequency field waveform inversion, simultaneously work in coordination with computing by means of GPU and CPU to overcome time domain waveform inversion method due to the refutation strategy of Whole frequency band data and cause objective function to be easily absorbed in the defect of local minimum, and avoid frequency field and the huge EMS memory occupation amount of Laplace domain waveform inversion, the problem that 3D processing is difficult to carry out.
By technical scheme of the present invention, achieve following technique effect: (1) is improved for the data processing of time domain waveform inversion, not only ensure the convergence of inverting, and the method for relative frequency territory and Laplace domain, greatly reduce EMS memory occupation amount; (2) adopt GPU/CPU to work in coordination with parallel accelerate computing method, at utmost make the calculation step of hybrid domain waveform inversion obtain parallel processing, thus significantly improve counting yield, save computing time in a large number.
In the present invention, employ the terms such as CPU, GPU, and employ the word such as function, DFT, for computer realm technician, these terms are apparent, CPU represents central processing unit, is a kind of electronic component with general-purpose computations ability, and this common series products is produced by Intel or AMD, GPU representative of graphics processor, also known as video card, general-purpose computations ability for CPU, but GPU is absorbed in image, the three-dimensional computations abilities such as figure, there is the parallel computing unit with computing power far above CPU number in unit GPU, even if the computing unit (also known as stream handle) that the low side graphic process unit on market has is also hundreds of, function due to GPU is the image of compound display millions of pixels on screen, namely have millions of tasks simultaneously and need parallel processing, therefore be namely designed in circuit design can a lot of task of parallel processing for GPU.Relevant to CPU, GPU, CPU obtains data and realizes data transmission by internal memory, and the data transmission of GPU is then pass through video memory; Wherein, what DFT represented is a kind of mathematical processing methods: Discrete Fourier Transform, namely discrete Fourier transformation, being the form all discrete in time domain and frequency domain of continuous fourier transform, is the sampling at discrete time Fourier transform (DTFT) frequency domain by the sampling transformation of time-domain signal.
For reaching technique effect of the present invention, the present invention is achieved through the following technical solutions:
CPU/GPU works in coordination with the hybrid domain full waveform inversion method of parallel computation, comprises the steps: 1, and the part of just drilling of waveform inversion is put into time domain; 2, adopt discrete Fourier transformation to extract frequency field wave field, and forward frequency field to and carry out inverting.Wherein, above-mentioned steps 1,2 completes calculating in GPU.
Concrete, above-mentioned steps is realized by following inversion method:
The first step, chooses optimum inverting frequency, and reads initial model, inverted parameters and other inverting desired parameters; For those skilled in the art, these parameters are just determined before inverting starts;
Second step, carry out group of frequencies circulation or single-frequency point circulation according to the inverting frequency set, recycle design adopts iterative loop;
3rd step, enters shot point circulation, reads single big gun data, determine the imaging pore diameter range of single big gun, and be transferred in GPU video memory by corresponding velocity field by internal memory;
4th step, carries out main story with Spatial higher order time domain propagation operator to focus, eliminates boundary echo simultaneously; Within each propagation moment of main story, record earth's surface accepts wave field, utilize DFT to extract the frequency field main story wave field of corresponding inverting frequency simultaneously, namely this moment wave field is calculated to the contribution of ongoing frequency territory wave field to the wave field snapshot in each moment, these contributions are added up;
5th step, in time domain, accepts wave field to earth's surface and actual big gun record does difference, calculates the target function value that this inverting frequency is corresponding simultaneously;
6th step, using residual error wave field as focus, utilizes time domain propagation operator to carry out anti-pass to it, utilizes DFT to extract the frequency field wave field of corresponding inverting frequency, and namely cumulative all moment time domain wave fields are to the contribution of ongoing frequency territory wave field;
7th step, calculates gradient fields corresponding under forward gun is recorded in ongoing frequency or group of frequencies;
8th step, to next big gun collection, repeat the 3rd step to the 7th step, the cumulative gradient fields of this common-shot-gather under ongoing frequency, until the circulation of all shot points is complete, then asks for maximal value the choosing for Optimal Step Size in this gradient fields simultaneously;
9th step, utilizes curve-parabola-fitting method to ask for Optimal Step Size;
Tenth step, with Optimal Step Size and gradient fields renewal speed field, continues to leave in GPU video memory by this velocity field, for next iteration, repeat above-mentioned steps two to nine, until meet stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency or group of frequencies;
11 step, when all frequencies or the equal inverting of group of frequencies complete, the velocity field of final updated is transferred to internal memory by GPU video memory, and writes in storer, complete refutation process;
Wherein, described 4th step carries out parallel accelerate computing to the process of the tenth step by GPU.
By said process, method of the present invention significantly improves the efficiency of full waveform inversion, is in particular in:
Time domain propagation operator adopts Spatial higher order, time Finite Difference Scheme of Second Order, wave field is propagated and gradient fields each net point in process such as to ask for be separate, all net points can parallel computation, namely each thread only calculates the value of one or more net point, GPU parallel computation core is many, fine grained parallel computing can be carried out expeditiously, make method of the present invention be achieved acceleration process to hybrid domain waveform inversion.
On the other hand, for the feature of hybrid domain algorithm, the present invention adopts the thought of matrix burst.Because high-order finite difference method algorithm, the calculating of each net point wave field value needs to use the data of multiple net point around, it is very high that internal memory reads redundance, internal storage data is reused by GPU shared drive (share memory), obviously can reduce internal memory and read redundance (read redundancy), use GPU shared drive effectively to improve speed-up ratio with the inventive method; Meanwhile, compared with internal memory, share memory is at GPU chip internal, have very high bandwidth, its speed, close to the speed of register, has very high data transmission efficiency relative to internal memory, existing multiplexing factually by shared drive logarithm, significantly improve inverting execution efficiency.
Wherein, eliminating boundary echo in described 4th step can be accomplished in several ways, and these methods are that this area widely adopts, and includes but not limited to adopt perfectly matched layer condition, sponge absorbing boundary condition, random boundary condition, absorption etc.Preferably, the present invention adopts perfectly matched layer condition, and namely PML boundary condition eliminates boundary echo.
In the present invention, the formula asking for gradient fields in described 7th step is
Wherein, P
ffor focus main story wave field, P
bfor residual error anti-pass wave field, v is the speed of underground medium, x
sand x
rbe the horizontal level of focus and wave detector respectively, ω is the frequency of current inverting, and * represents and gets conjugation to variable.
In the present invention, 9th step builds para-curve can adopt mathematically existing Parabolic Fit algorithm realization Parabolic Fit, the following method of preferred employing builds para-curve, and ask for minimum point: by the value of the longitudinal axis objective definition function of parabolic coordinates, transverse axis is step-length, choose three step-lengths, and corresponding target function value; Known α
0target function value when=0, afterwards by iterative search α
1and α
2value, its target function value F (α) meets F (α
1) < F (α
0) and F (α
2) > F (α
1) condition; With (α
0, F (α
0)), (α
1, F (α
1)), (α
2, F (α
2)) three groups of points determine a para-curve, using its minimum point as Optimal Step Size.
In hybrid domain full waveform inversion method of the present invention, 4th step can also comprise the steps: in the time circulation of wave field main story, DFT operator is utilized to ask for the contribution of each moment to frequency field wave field, until all moment all propagate complete, obtain current inverting frequency or frequency field main story wave field corresponding to group of frequencies.
Similar to the above, 6th step can also comprise the steps: using residual error wave field as focus anti-pass, in the circulation of this time, DFT operator is utilized to ask for each moment residual error anti-pass wave field to the contribution of frequency field wave field, until all moment all propagate complete, obtain current inverting frequency or frequency field anti-pass wave field corresponding to group of frequencies.
Method of the present invention, CPU and GPU is collaborative computing, and not all computing all completes at GPU, concrete, the cooperated computing mode of the two is call GPU in CPU to perform Spatial higher order finite difference CPU module, PML absorbing boundary CPU module, earth's surface accepts the CPU module that wave field extracts, discrete Fourier change CPU module during wave field main story, the CPU module that time domain residual error is asked for, the CPU module that target function value is asked for, the CPU module that residual error wave field loads, discrete Fourier transformation CPU module during the anti-pass of residual error wave field, the CPU module that gradient fields is asked for, ask for the CPU module of gradient fields maximal value, the CPU module that Optimal Step Size is asked for, the CPU module of renewal speed model.
Above-mentioned module, easily by those skilled in the art are understood, that the program performing described function (can adopt any programming language to realize, normally C or C++) or the packaged routine package of multiple program, manufacturer is researched and developed by utilizing GPU, the packaged function library that such as NVIDIA or AMD provide is called and is revised to realize, and CPU calls these modules and makes the instruction in GPU execution module to realize corresponding function.
By above-mentioned technological improvement, full waveform inversion method of the present invention can carry out the hybrid domain waveform inversion of two and three dimensions, namely does in time domain and just drills, and utilizes DFT to forward frequency field to afterwards and does inverting.Method and technology scheme of the present invention efficiently solves the convergence problem of Conventional Time territory method based on the collaborative parallel computation of GPU and CPU, and has avoided the problem that conventional three-dimensional frequency field and Laplace domain waveform inversion EMS memory occupation amount cannot meet; Meanwhile, the method for relative forefathers, method calculation procedure of the present invention decreases a DFT and time IDFT.
Accompanying drawing explanation
Fig. 1 is the process flow diagram of technical solution of the present invention.
Fig. 2 is the process flow diagram of twice DFT in technical scheme the 4th step and the 6th step: (a) DFT extracts frequency domain focus main story wave field; B () DFT extracts frequency-domain residual anti-pass wave field.
Fig. 3 is the function curve schematic diagram that in technical scheme the 9th step, curve-parabola-fitting method asks for Optimal Step Size.
Fig. 4,5,6 is initial velocity model, true velocity model and inversion speed model disclosed in embodiment respectively.
Wherein background color is the picture frame of grey is the calculating link adopting GPU to achieve acceleration.
Embodiment
Be described in detail embodiments of the present invention below in conjunction with accompanying drawing, these accompanying drawings, as a part of the present invention, indicate a kind of specific implementation of the present invention.On the basis understanding accompanying drawing of the present invention and preferred embodiment and connotation of the present invention, known the present invention allows multiple multi-form embodiment, therefore without departing from the scope of the present invention, also other embodiment can be had to be used and build, and may be made other optionally change.
As shown in Figure 1, hybrid domain full waveform inversion method of the present invention, has and comprises:
The first step, chooses optimum inverting frequency, and reads initial model, inverted parameters and other parameters.
Second step, carries out frequency (group) circulation according to the inverting frequency set; Iterative loop in frequency (group) circulation.
3rd step, enters shot point circulation, and reads single big gun data, is specified to picture pore diameter range, and by velocity field again CPU internal memory be transferred in GPU video memory.
4th step, as shown in Figure 2 a, carries out main story with time domain propagation operator to focus, and record earth's surface accepts record, utilizes DFT to extract the frequency field main story wave field of corresponding inverting frequency simultaneously.This moment wave field is calculated to the contribution of frequency wave field to the wave field snapshot in each moment, and adds up it.Because each spatial point is separate, calculate so the realization of DFT is applicable to GPU fine grained parallel very much.
5th step, in time domain, does difference to earth's surface acceptance record and actual big gun record, calculates the target function value that this inverting frequency is corresponding simultaneously.
6th step, as shown in Figure 2 b, with residual error wave field for focus, utilizes time domain propagation operator to carry out anti-pass to it, utilizes DFT to extract the frequency field wave field of corresponding inverting frequency equally, and cumulative all moment time domain wave fields are to the contribution of frequency wave field.
7th step, utilizes formula:
Calculate corresponding gradient fields under forward gun is recorded in ongoing frequency (group), and ask for the maximal value of gradient fields.
8th step, repeats the 3rd step to the 7th step to next big gun collection, simultaneously the cumulative gradient fields of this common-shot-gather under ongoing frequency, until the circulation of all shot points is complete.
9th step, utilizes the Parabolic Fit in Fig. 3 to ask for Optimal Step Size.Wherein the longitudinal axis is the value of objective function, and transverse axis is step-length, chooses three step-lengths and corresponding target function value; Target function value during known α 0=0, afterwards by the value of iterative search α 1 and α 2, its target function value F (α) meets F (a1) < F (α 0) and the condition of F (α 2) > F (α 1); Further, construct a para-curve with these three groups of points, and ask for minimum point, in this, as Optimal Step Size.
Tenth step, renewal speed field is come by Optimal Step Size and gradient fields, and this velocity field is continued to leave in GPU video memory, for next iteration, until meet stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency (group), namely repeat second step to the tenth step.
11 step, when the equal inverting of all frequencies is complete, is transferred to CPU by the velocity field of final updated by GPU video memory, and writes in storer.Now, inverting is complete.
As shown in Fig. 4-Fig. 6, the present embodiment adopts CAS (Chinese Academy of Sciences writes a Chinese character in simplified form) mode letters to verify the validity of hybrid domain method.
Selected model size is 2000m × 500m, and sampling interval is 5m in length and breadth, and the data of selected big gun are shot point number is 65 big guns, Mei Bao 400 road, and record length is 1.6s, and time sampling interval is 1ms, and shot interval is 30m.
It is that 23 frequencies are chosen at interval with 1.5Hz that inverted parameters is set to inverting frequency from 2Hz to 35Hz, and inverting circulates according to group of frequencies, and each group of frequencies chooses two frequencies.
Stopping criterion for iteration is:
(E(m
k+1)-E(m
k))/E(m
k+1)＜0.001
Wherein, E (mk) and E (mk+1) is the target function value of kth time and kth+1 time respectively.
Have employed for the computer of inverting the CPU that model is i7-2600 in the present embodiment, model is the GPU of GTX680.Can be found out by accompanying drawing, adopt method of the present invention, inversion result portrays true velocity model well.
Concrete implementation is as follows: first big gun data be read into calculator memory from hard disk, and determine the pore size corresponding to this big gun, then by the speed in this aperture with when forward gun data are from CPU memory copying to GPU video memory, tstep_kernel (high-order finite difference method GPU executive routine is called afterwards successively in CPU, for carrying out wave traveling in time domain), pml_kernel (PML absorbing boundary GPU executive routine, for absorbing the reflection from border), (earth's surface accepts the GPU program that wave field extracts to record_kernel, wave field is accepted) for extracting corresponding earth's surface from wave field snapshot, DFT_kernel (discrete Fourier change GPU executive routine, for asking for current time time domain wave field to the contribution of ongoing frequency territory wave field), residual_kernel (the GPU program that time domain residual error is asked for, for asking for the difference of prediction time domain wavefield and true wavefield), obj_kernel (the GPU program that target function value is asked for, for calculating the target function value under present speed model), addresi_kernel (the GPU executive routine that residual error wave field loads, for when the anti-pass of residual error wave field, residual error wave field is loaded with focus form), tstep_kernel and pml_kernel (for propagating residual error wave field), DFT2_kernel (for calculating the contribution of current time zone wave field to ongoing frequency territory wave field), grad_kernel (the GPU program that gradient fields is asked for, for calculating the gradient fields that current input big gun set pair is answered), step_cuda (the GPU program that step-length is asked for, multiple kernel program module is comprised in this module, Optimal Step Size is asked for) for the method with Parabolic Fit, updatevel_kernel (the GPU program of renewal speed model, gradient fields and Optimal Step Size is utilized to come renewal speed field).All GPU program call complete after, a complete iterative process just finishes, and enters into next iteration afterwards, until meet stopping criterion for iteration or reach maximum iteration time.After iteration terminates, namely enter next inverting group of frequencies, until all inverting frequencies are finished.
In this process, GPU assume responsibility for most calculated amount, thus counting yield is significantly promoted.
In above-mentioned concrete enforcement, employ multiple function name, as pml_kernel etc., all realize by the function calling and revise in CUDA function library that NVIDIA company (i.e. GTX680 developer) openly provides.
And the experiment display that above-mentioned CAS model is carried out on the CPU of simple i7-2600, CPU cannot complete this inverting flow process performing in the full waveform inversion identical time with above-mentioned cooperated computing, can not get result; Even if in fact adopt above-mentioned same algorithm to use merely CPU to carry out full waveform inversion, inverting cannot can accept to complete in the time, based on above-mentioned hardware platform, adopt identical algorithm, the simple calculating process single iteration based on CPU needs 5034s consuming time, and method single iteration of the present invention only 241s consuming time, method of the present invention significantly improves counting yield as can be seen here, has good speed-up ratio.
The above-mentioned embodiment comprising accompanying drawing is only example of the present invention, and only there is provided for diagram of the present invention.The embodiment of other variation is apparent for technician in the art.Therefore, protection scope of the present invention is as the criterion by content disclosed in claim and instructions, and is not limited to given embodiment.
Claims (7)
1.CPU/GPU works in coordination with the hybrid domain full waveform inversion method of parallel computation, it is characterized in that described method comprises the steps: 1, and the part of just drilling of waveform inversion is put into time domain; 2, adopt discrete Fourier transformation to extract frequency field wave field, and forward frequency field to and carry out inverting; Wherein, above-mentioned steps 1,2 completes calculating in GPU; Described step specifically comprises:
The first step, chooses optimum inverting frequency, and reads initial model, inverted parameters and other inverting desired parameters;
Second step, carry out group of frequencies circulation or single-frequency point circulation according to the inverting frequency set, recycle design adopts iterative loop;
3rd step, enters shot point circulation, reads single big gun data, determine the imaging pore diameter range of single big gun, and be transferred in GPU video memory by corresponding velocity field by internal memory;
4th step, carries out main story with Spatial higher order time domain propagation operator to focus, eliminates boundary echo simultaneously; Within each propagation moment of main story, record earth's surface accepts wave field, utilize DFT to extract the frequency field main story wave field of corresponding inverting frequency simultaneously, namely this moment wave field is calculated to the contribution of ongoing frequency territory wave field to the wave field snapshot in each moment, these contributions are added up;
5th step, in time domain, accepts wave field to earth's surface and actual big gun record does difference, calculates the target function value that this inverting frequency is corresponding simultaneously;
6th step, using residual error wave field as focus, utilizes time domain propagation operator to carry out anti-pass to it, utilizes DFT to extract the frequency field wave field of corresponding inverting frequency, and namely cumulative all moment time domain wave fields are to the contribution of ongoing frequency territory wave field;
7th step, calculates gradient fields corresponding under forward gun is recorded in ongoing frequency or group of frequencies;
8th step, to next big gun collection, repeat the 3rd step to the 7th step, the cumulative gradient fields of this common-shot-gather under ongoing frequency, until the circulation of all shot points is complete, then asks for the maximal value in this gradient fields, choosing for Optimal Step Size simultaneously;
9th step, utilizes curve-parabola-fitting method to ask for Optimal Step Size;
Tenth step, with Optimal Step Size and gradient fields renewal speed field, continues to leave in GPU video memory by this velocity field, for next iteration, repeat above-mentioned steps two to nine, until meet stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency or group of frequencies;
11 step, when all frequencies or the equal inverting of group of frequencies complete, the velocity field of final updated is transferred to internal memory by GPU video memory, and writes in storer, complete refutation process;
Wherein, described 4th step carries out parallel accelerate computing to the process of the tenth step by GPU.
2. hybrid domain full waveform inversion method according to claim 1, it is characterized in that described 4th step adopts perfectly matched layer, namely PML boundary condition eliminates boundary echo.
3. hybrid domain full waveform inversion method according to claim 1, is characterized in that the formula asking for gradient fields in described 7th step is
Wherein, P
ffor focus main story wave field, P
bfor residual error anti-pass wave field, v is the speed of underground medium, x
sand x
rbe the horizontal level of focus and wave detector respectively, ω is the frequency of current inverting, and * represents and gets conjugation to variable.
4. hybrid domain full waveform inversion method according to claim 1, it is characterized in that described 9th step profit builds para-curve with the following method, and ask for minimum point: by the value of the longitudinal axis objective definition function of parabolic coordinates, transverse axis is step-length, choose three step-lengths, and corresponding target function value; Known α
0target function value when=0, afterwards by iterative search α
1and α
2value, its target function value F (α) meets F (α
1) < F (α
0) and F (α
2) > F (α
1) condition; With (α
0, F (α
0)), (α
1, F (α
1)), (α
2, F (α
2)) three groups of points determine a para-curve, using its minimum point as Optimal Step Size.
5. the hybrid domain full waveform inversion method according to any one of claim 1-4, it is characterized in that described 4th step also comprises the steps: in the time circulation of wave field main story, DFT operator is utilized to ask for the contribution of each moment to frequency field wave field, until all moment all propagate complete, obtain current inverting frequency or frequency field main story wave field corresponding to group of frequencies.
6. the hybrid domain full waveform inversion method according to any one of claim 1-4, it is characterized in that described 6th step also comprises the steps: using residual error wave field as focus anti-pass, in the circulation of this time, DFT operator is utilized to ask for each moment residual error anti-pass wave field to the contribution of frequency field wave field, until all moment all propagate complete, obtain current inverting frequency or frequency field anti-pass wave field corresponding to group of frequencies.
7. hybrid domain full waveform inversion method according to claim 1, it is characterized in that when CPU/GPU works in coordination with parallel computation, the cooperated computing mode of the two is call GPU in CPU to perform Spatial higher order finite difference CPU module, PML absorbing boundary CPU module, earth's surface accepts the CPU module that wave field extracts, discrete Fourier change CPU module during wave field main story, the CPU module that time domain residual error is asked for, the CPU module that target function value is asked for, the CPU module that residual error wave field loads, discrete Fourier transformation CPU module during the anti-pass of residual error wave field, the CPU module that gradient fields is asked for, ask for the CPU module of gradient fields maximal value, the CPU module that Optimal Step Size is asked for, the CPU module of renewal speed model.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310013786.6A CN103135132B (en) | 2013-01-15 | 2013-01-15 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310013786.6A CN103135132B (en) | 2013-01-15 | 2013-01-15 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN103135132A CN103135132A (en) | 2013-06-05 |
| CN103135132B true CN103135132B (en) | 2015-07-01 |
Family
ID=48495211
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201310013786.6A Expired - Fee Related CN103135132B (en) | 2013-01-15 | 2013-01-15 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN103135132B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105891888A (en) * | 2016-03-28 | 2016-08-24 | 吉林大学 | Multi-domain frequency-division parallel multi-scale full-waveform inversion method |
Families Citing this family (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570090B (en) * | 2013-10-29 | 2017-07-28 | 中国石油化工股份有限公司 | The extraction of full waveform inversion noise filter operator and the method filtered using its noise |
| CN105319581B (en) * | 2014-07-31 | 2018-01-16 | 中国石油化工股份有限公司 | A kind of efficient time-domain full waveform inversion method |
| CN105445798B (en) * | 2014-08-21 | 2018-08-07 | 中国石油化工股份有限公司 | A kind of full waveform inversion method and system based on gradient processing |
| CN106199692A (en) * | 2015-05-30 | 2016-12-07 | 中国石油化工股份有限公司 | A kind of wave equation inverse migration method based on GPU |
| CN105005076B (en) * | 2015-06-02 | 2017-05-03 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105005072B (en) * | 2015-06-02 | 2016-08-17 | 中国科学院地质与地球物理研究所 | The PML border dimensionally seismic wave propagating mode utilizing CUDA intends method |
| CN106294273B (en) * | 2015-06-05 | 2020-01-10 | 中国石油化工股份有限公司 | Communication method and system of CPU and coprocessor |
| US10331195B2 (en) * | 2016-06-06 | 2019-06-25 | Qualcomm Incorporated | Power and performance aware memory-controller voting mechanism |
| CN107656306B (en) * | 2016-07-26 | 2019-02-19 | 中国石油化工股份有限公司 | Full waveform inversion parallel calculating method and system |
| CN107656307B (en) * | 2016-07-26 | 2019-02-19 | 中国石油化工股份有限公司 | Full waveform inversion calculation method and system |
| CN107918145B (en) * | 2016-10-10 | 2020-09-15 | 中国石油化工股份有限公司 | Parallelization processing method and system for seismic gun energy |
| CN106950596B (en) * | 2017-04-11 | 2019-02-26 | 中国石油大学(华东) | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
| CN107728199B (en) * | 2017-09-22 | 2019-05-31 | 中国地质大学（北京） | Based on the parallel multi -components anisotropy pre-stack time migration accelerated method of more GPU |
| CN108594299B (en) * | 2018-02-28 | 2019-12-31 | 中国科学院地质与地球物理研究所 | Intelligent early warning method, device and system for high-speed rail |
| CN110888158B (en) * | 2018-09-10 | 2021-12-31 | 中国石油化工股份有限公司 | Full waveform inversion method based on RTM constraint |
| CN109523022B (en) * | 2018-11-13 | 2022-04-05 | Oppo广东移动通信有限公司 | Terminal data processing method and device and terminal |
| CN109633742B (en) * | 2019-01-08 | 2020-03-27 | 中国科学院地质与地球物理研究所 | Full waveform inversion method and device |
| CN111638551A (en) * | 2019-03-01 | 2020-09-08 | 中国石油天然气集团有限公司 | Seismic first-motion wave travel time chromatography method and device |
| CN112099936A (en) * | 2019-06-17 | 2020-12-18 | 中国石油天然气集团有限公司 | Heterogeneous parallel computing implementation method and device for three-dimensional acoustic wave NPML algorithm |
| CN112578431B (en) * | 2019-09-27 | 2024-04-09 | 中国石油化工股份有限公司 | Method and system for storing full waveform inversion wave field optimization in finite state |
| CN110895347A (en) * | 2019-12-09 | 2020-03-20 | 中国科学院地质与地球物理研究所 | Seismic wave forward modeling method and system |
| CN113504566B (en) * | 2021-06-01 | 2024-04-30 | 南方海洋科学与工程广东省实验室(湛江) | Wave equation-based seismic inversion method, system, device and medium |
| CN114814843B (en) * | 2022-06-27 | 2022-10-04 | 之江实验室 | Earth surface deformation inversion method and system based on CPU-GPU heterogeneous parallel |
| WO2024087002A1 (en) * | 2022-10-25 | 2024-05-02 | Saudi Arabian Oil Company | Methods and systems for determining attenuated traveltime using parallel processing |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101630018A (en) * | 2008-07-16 | 2010-01-20 | 中国石油天然气股份有限公司 | Method for processing seismic exploration data in process of controlling full acoustic wave equation inversion |
| CN101688926A (en) * | 2007-06-26 | 2010-03-31 | 申昌秀 | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6999880B2 (en) * | 2003-03-18 | 2006-02-14 | The Regents Of The University Of California | Source-independent full waveform inversion of seismic data |
| AU2006235820B2 (en) * | 2005-11-04 | 2008-10-23 | Westerngeco Seismic Holdings Limited | 3D pre-stack full waveform inversion |
| US7725266B2 (en) * | 2006-05-31 | 2010-05-25 | Bp Corporation North America Inc. | System and method for 3D frequency domain waveform inversion based on 3D time-domain forward modeling |
-
2013
- 2013-01-15 CN CN201310013786.6A patent/CN103135132B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101688926A (en) * | 2007-06-26 | 2010-03-31 | 申昌秀 | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging |
| CN101630018A (en) * | 2008-07-16 | 2010-01-20 | 中国石油天然气股份有限公司 | Method for processing seismic exploration data in process of controlling full acoustic wave equation inversion |
Non-Patent Citations (1)
| Title |
|---|
| 苏超,等.时间域声波全波形反演及GPU加速.《中国地球物理2012》.2012,第486页. * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105891888A (en) * | 2016-03-28 | 2016-08-24 | 吉林大学 | Multi-domain frequency-division parallel multi-scale full-waveform inversion method |
| CN105891888B (en) * | 2016-03-28 | 2017-03-08 | 吉林大学 | Multiple domain divides multiple dimensioned full waveform inversion method parallel |
Also Published As
| Publication number | Publication date |
|---|---|
| CN103135132A (en) | 2013-06-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN103135132B (en) | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing | |
| Huthwaite | Accelerated finite element elastodynamic simulations using the GPU | |
| Wu et al. | A procedure for 3D simulation of seismic wave propagation considering source‐path‐site effects: Theory, verification and application | |
| Weiss et al. | Solving 3D anisotropic elastic wave equations on parallel GPU devices | |
| Abdelkhalek et al. | Fast seismic modeling and reverse time migration on a GPU cluster | |
| CN106842320B (en) | The parallel 3-D seismics wave field generation method of GPU and system | |
| CN104570082B (en) | Extraction method for full waveform inversion gradient operator based on green function characterization | |
| CN103105623B (en) | Data waveform processing method in seismic exploration | |
| CN107783190A (en) | A kind of least square reverse-time migration gradient updating method | |
| Poul et al. | Efficient time-domain deconvolution of seismic ground motions using the equivalent-linear method for soil-structure interaction analyses | |
| CN107894618B (en) | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm | |
| CN105005072B (en) | The PML border dimensionally seismic wave propagating mode utilizing CUDA intends method | |
| CN105891888A (en) | Multi-domain frequency-division parallel multi-scale full-waveform inversion method | |
| CN108387933A (en) | A kind of method, apparatus and system of definitely interval quality factors | |
| CN107561585A (en) | A kind of multinuclear multi-node parallel 3-D seismics wave field generation method and system | |
| CN105911584B (en) | Implicit staggered-grid finite difference elastic wave numerical simulation method and device | |
| Yuan et al. | FUNWAVE‐GPU: Multiple‐GPU acceleration of a Boussinesq‐type wave model | |
| CN104965222B (en) | Three-dimensional longitudinal wave impedance full-waveform inversion method and apparatus | |
| Xue et al. | An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging | |
| CN102830431B (en) | Self-adaption interpolating method for real ground-surface ray tracking | |
| CN103698810A (en) | Hybrid network minimum travel time ray tracing tomography method | |
| Chabot et al. | A high-order discontinuous Galerkin method for 1D wave propagation in a nonlinear heterogeneous medium | |
| CN109490948A (en) | Seismoacoustics wave equation vector parallel calculating method | |
| Xie et al. | GPU acceleration of time gating based reverse time migration using the pseudospectral time-domain algorithm | |
| Wang et al. | An optimized parallelized SGFD modeling scheme for 3D seismic wave propagation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20150701 Termination date: 20190115 |