CN103135132A - 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
- CN103135132A CN103135132A CN2013100137866A CN201310013786A CN103135132A CN 103135132 A CN103135132 A CN 103135132A CN 2013100137866 A CN2013100137866 A CN 2013100137866A CN 201310013786 A CN201310013786 A CN 201310013786A CN 103135132 A CN103135132 A CN 103135132A
- Authority
- CN
- China
- Prior art keywords
- frequency
- gpu
- field
- wave 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.)
- Granted
Links
Images
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 wave shape method of inversion, relate in particular to by CPU/GPU and work in coordination with the Full wave shape inversion method that parallel computation is carried out on hybrid domain, belong to computer realm.
Background technology
The difficulty of China western part and the exploitation of overseas original oil zone is increasing at present, and the degree of depth is constantly deepened, and is very large to high-precision geophysical imaging and inversion method demand; And high-precision formation method finally also will rely on the structure of high precision velocity model.Inversion method, amplitude inversion method and Full wave shape inverting three class methods when the difference that the inversion method of rate pattern uses according to inverting information is divided into walking.Due to dynamics and the kinematics information of Full wave shape inversion method comprehensive utilization pre-stack seismic wave field, can rebuild accurately underground velocity field, become the important directions of exploration geophysics research both at home and abroad at present.
Existing Full wave shape inverting (Full Waveform Inversion, abbreviation FWI) method comprises the time domain Full wave shape inverting based on the Generalized Least Square inversion theory that Tarantola proposed in 1984, Pratt subsequently, the people such as Shin have been developed respectively the waveform inversion of frequency field and Laplace domain on its basis.The Full wave shape inverting directly utilizes seismologic record inverting underground medium parameter, utilize fully the shape informations such as reflection, refraction, under certain condition, can the Simultaneous Inversion velocity of longitudinal wave, the relative variation of density and shear wave velocity, thereby obtaining the higher rate pattern of precision, is the important directions of at present domestic and international exploration geophysics research.Although the Full wave shape inverting has many technical advantages, due to Full wave shape Inversion Calculation amount huge (being about the 30-100 of reverse-time migration doubly), be difficult to realize 3D processing, cause this method can't adapt to industrial needs.
In existing research, be divided three classes from account form for the Full wave shape inversion method: the first kind is that (Tarantola 1984 for the time domain waveform inversion, 1986), the method is utilized Whole frequency band data compute gradient field, cause it easily to be absorbed in local extremum because the extreme of objective function is non-linear, thereby inverting is not restrained; Equations of The Second Kind is that (Pratt 1998 for the frequency field waveform inversion, 1999), the method is just done in frequency field and is drilled, pursue afterwards (group) inverting frequently from the low frequency to the high frequency in frequency field, the relative high frequency data, low-frequency data and rate pattern have better linear relationship, and the method can effectively be avoided the local minimum problem in the time domain refutation process; The method of just drilling in frequency field is divided into two kinds of direct method (as: LU decompose etc.) and process of iteration, because process of iteration efficient is low, frequency field is just being drilled direct methods such as adopting the LU decomposition more and is being found the solution, the advantage of this method is after impedance matrix is once decomposed, result after decomposition can be used for the calculating of all shot points, counting yield is high, this kind method also exists the drawback that is difficult to overcome, under three-dimensional case, LU decomposition EMS memory occupation amount and calculated amount are all very huge, cause the three-dimensional waveform inverting to implement unrealistic.The 3rd class is Laplace domain waveform inversion (Shin2008,2009), the method is on the basis of frequency field waveform inversion, inverted parameters has been increased decay factor, thus inverting be not only from the low frequency to the high frequency, also from the shallow-layer to the deep layer, thereby reduced the dependence of inversion result to initial model, but the method also faces the memory bottleneck the same with the frequency field method, and 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 research and the application work of waveform inversion, but are generally limited to two dimension, are difficult to use in three-dimensional.In general, the 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 of domestic waveform inversion is used and basically is in the starting stage, and is relatively backward, therefore, studies practical waveform inversion method the technological gap of filling up China is had great importance.
the US Patent No. 11/756 of 2007, 384, " System and method for 3D frequency domain waveform inversion based on 3D time-domain forward modeling ", in the disclosed technical scheme of this invention, provide a kind of waveform inversion method: at first utilize DFT to be transformed into frequency field the big gun numeric field data, then utilize time domain propagation operator and DFT calculated rate territory focus main story wave field, calculate the residual error wave field in frequency field, and further utilize inverse discrete Fourier transform (Inverse DiscreteFourier Transform, abbreviation IDFT) ask for the residual error wave field of time domain, afterwards, with the time domain propagation operator, the residual error wave field is carried out anti-pass, and utilize equally DFT to extract frequency field residual error anti-pass wave field, at last, ask for the gradient fields of current inverting frequency.
Compare with the early stage numerical procedure that adopts in this area, although this technical scheme has obtained certain progress, its process step relative complex is carried out efficient and is under some influence; And in hybrid domain FWI, the parallel granularity of main module consuming time (as DFT, time domain propagation operator etc.) is very thin, only adopts the CPU computing, fails farthest to make program parallelization, and counting yield can be greatly affected.Refutation strategy of the present invention is different from it, has lacked a DFT and an IDFT on the basis of this inverting flow process; And further, adopt the collaborative concurrent operation of GPU/CPU, and farthest make program parallelization, increased significantly the speed-up ratio with respect to the CPU program.
Summary of the invention
The invention discloses the hybrid domain Full wave shape inversion method of the collaborative parallel computation of a kind of CPU/GPU, can carry out the waveform inversion data of two and three dimensions and process.This method combines the advantage of time domain waveform inversion and frequency field waveform inversion, simultaneously working in coordination with computing by means of GPU and CPU has overcome the time domain waveform inversion method and causes objective function easily to be absorbed in the defective of local minimum due to the refutation strategy of Whole frequency band data, and avoided the huge EMS memory occupation amount of frequency field and Laplace domain waveform inversion, the problem that 3D processing is difficult to carry out.
By technical scheme of the present invention, realized following technique effect: (1) is processed for the data of time domain waveform inversion and is improved, not only guarantee the convergence of inverting, and the method for relative frequency territory and Laplace domain, the EMS memory occupation amount greatly reduced; (2) adopt the collaborative parallel speed-up computation method of GPU/CPU, at utmost make the calculation step of hybrid domain waveform inversion obtain parallel processing, thereby significantly improve counting yield, save computing time in a large number.
in the present invention, the terms such as CPU, GPU have been used, and the words such as function, DFT have been used, for the 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, claim again video card, general-purpose computations ability for CPU, but GPU is absorbed in image, the three-dimensional computations abilities such as figure, has the parallel computing unit with computing power far above the CPU number in the GPU of unit, even the computing unit that the low side graphic process unit on market has (claiming again stream handle) is also hundreds of, because the function of GPU is the image of millions of pixels of compound display on screen, namely have simultaneously millions of tasks and need parallel processing, therefore but GPU namely is designed to a lot of tasks of parallel processing in circuit design.Relevant to CPU, GPU is, CPU obtains data and realizes data transmission by internal memory, and the data transmission of GPU is to pass through video memory; Wherein, what DFT represented is a kind of mathematical processing methods: Discrete Fourier Transform, namely discrete Fourier transformation, be continuous fourier transform discrete form all on time domain and frequency domain, the sampling of time-domain signal is transformed to sampling at discrete time Fourier transform (DTFT) frequency domain.
For reaching technique effect of the present invention, the present invention is achieved through the following technical solutions:
The hybrid domain Full wave shape inversion method of the collaborative parallel computation of CPU/GPU comprises the steps: 1, and just drilling of waveform inversion partly is put into time domain; 2, adopt discrete Fourier transformation to extract the frequency field wave field, and forward frequency field to and carry out inverting.Wherein, above-mentioned steps 1, the 2nd is completed calculating in GPU.
Concrete, above-mentioned steps is realized by following inversion method:
The first step is chosen optimum inverting frequency, and is read initial model, inverted parameters and other inverting desired parameters; For those skilled in the art, these parameters are just determined before inverting begins;
Second step carries out group of frequencies circulation or single-frequency point circulation according to the inverting frequency that sets, and recycle design adopts iterative loop;
The 3rd step entered the shot point circulation, read single big gun data, determines the imaging pore diameter range of single big gun, and corresponding velocity field is transferred in the GPU video memory by internal memory;
The 4th step, with Spatial higher order time domain propagation operator, focus is carried out main story, eliminate simultaneously boundary echo; In each propagation constantly of main story, record the earth's surface and accept wave field, utilize simultaneously DFT to extract the frequency field main story wave field of corresponding inverting frequency, namely each wave field snapshot is constantly calculated this moment wave field to the contribution of current frequency field wave field, these contributions are added up;
The 5th step, in time domain, wave field and actual big gun record is accepted on the earth's surface done poorly, calculate simultaneously target function value corresponding to this inverting frequency;
The 6th step, as focus, utilize the time domain propagation operator to carry out anti-pass to it with the residual error wave field, utilize DFT to extract the frequency field wave field of corresponding inverting frequency, all time domain wave fields contributions to current frequency field wave field constantly namely add up;
In the 7th step, calculate the gradient fields that is recorded in correspondence under current frequency or group of frequencies when forward gun;
The 8th step to next big gun collection, repeated for the 3rd step to the 7th step, the gradient fields of this common-shot-gather under current frequency that add up simultaneously, until all shot point circulations are complete, the maximal value of then asking in this gradient fields is used for choosing of Optimal Step Size;
In the 9th step, utilize the Parabolic Fit method to ask for Optimal Step Size;
In the tenth step, with Optimal Step Size and gradient fields renewal speed field, this velocity field is continued to leave in the GPU video memory, be used for next iteration, repeat above-mentioned steps two to nine, until satisfy stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency or group of frequencies;
The 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 the GPU video memory, and in write store, completes refutation process;
Wherein, the process in described the 4th step to the tenth step is walked abreast by GPU and accelerates computing.
By said process, method of the present invention has significantly improved the efficient of Full wave shape inverting, is in particular in:
The time domain propagation operator adopts Spatial higher order, time Finite Difference Scheme of Second Order, wave field is propagated and gradient fields asks for etc. that in process, each net point is separate, all net points can parallel computation, it is the value that each thread only calculates one or more net points, GPU parallel computation core is many, can carry out expeditiously the fine grained parallel computing, make method of the present invention be achieved the acceleration of hybrid domain waveform inversion is processed.
On the other hand, for the characteristics of hybrid domain algorithm, the present invention adopts the thought of matrix burst.Because high-order limited difference algorithm, the calculating of each net point wave field value needs to use the data of a plurality of net points on every side, it is very high that internal memory reads redundance, (share memory) reuses internal storage data by the GPU shared drive, can obviously reduce internal memory and read redundance (read redundancy), so method of the present invention uses the GPU shared drive effectively to promote speed-up ratio; Simultaneously, compare with internal memory, share memory is at the GPU chip internal, have very high bandwidth, its speed has very high data transmission efficiency close to the speed of register with respect to internal memory, existing multiplexing factually by the shared drive logarithm, significantly improve inverting and carried out efficient.
Wherein, eliminating boundary echo in described the 4th step can be accomplished in several ways, and these methods are that widely adopt this area, include but not limited to adopt optimum matching layer condition, spongy absorbing boundary condition, RANDOM BOUNDARY condition, absorption etc.Preferably, the present invention adopts optimum matching layer condition, and namely the PML boundary condition is eliminated boundary echo.
In the present invention, the formula of asking for gradient fields in described the 7th step is
Wherein, P
fBe focus main story wave field, P
bBe residual error anti-pass wave field, v is the speed of underground medium, x
sAnd x
rBe respectively the horizontal level of focus and wave detector, ω is the frequency of current inverting, and * represents variable is got conjugation.
In the present invention, the 9th step built para-curve and can adopt on mathematics existing Parabolic Fit algorithm to realize Parabolic Fit, the following method of preferred employing builds para-curve, and ask for minimum point: with 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 α
0The target function value of=0 o'clock is afterwards by iterative search α
1And α
2Value, its target function value F (α) satisfies 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, with its minimum point as Optimal Step Size.
In the hybrid domain Full wave shape method of inversion of the present invention, the 4th step can also comprise the steps: in the circulation of the time of wave field main story, utilize the DFT operator to ask for each moment to the contribution of frequency field wave field, until that all moment all propagate is complete, obtain frequency field main story wave field corresponding to current inverting frequency or group of frequencies.
Similar to the above, the 6th step can also comprise the steps: with the residual error wave field as the focus anti-pass, in this time circulation, utilize the DFT operator to ask for each residual error anti-pass wave field contribution to the frequency field wave field constantly, until that all moment all propagate is complete, obtain frequency field anti-pass wave field corresponding to current inverting frequency or group of frequencies.
method of the present invention, CPU and GPU are collaborative computings, be not that all computings are all completed at GPU, concrete, collaborative account form both is to call GPU in CPU to carry out Spatial higher order finite difference GPU module, PML absorbing boundary GPU module, the GPU module that wave field extracts is accepted on the earth's surface, discrete Fourier during the wave field main story changes the GPU module, the GPU module that the time domain residual error is asked for, the GPU module that target function value is asked for, the GPU module that the residual error wave field loads, discrete Fourier transformation GPU module during the anti-pass of residual error wave field, the GPU module that gradient fields is asked for, ask for the peaked GPU module of gradient fields, the GPU module that Optimal Step Size is asked for, the GPU module of renewal speed model.
Above-mentioned module, easily understood by those skilled in the art, that the program of carrying out described function (can adopt any programming language to realize, normally C or C++) or the packaged routine package of a plurality of program, by utilizing GPU to research and develop manufacturer, for example the packaged function library that provides of NVIDIA or AMD is called and is revised and realizes, CPU calls these modules and makes instruction in the GPU execution module to realize corresponding function.
By above-mentioned technological improvement, Full wave shape inversion method of the present invention can carry out the hybrid domain waveform inversion of two and three dimensions, namely just does in time domain and drills, and utilizes afterwards DFT to forward frequency field to and does inverting.Method and technology scheme of the present invention has solved the convergence problem of conventional time domain method effectively 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 can't satisfy; Simultaneously, relative forefathers' method, method calculation procedure of the present invention has reduced a DFT and an IDFT.
Description of drawings
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 the 4th step of technical scheme 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 during technical scheme the 9th goes on foot, the Parabolic Fit method is asked for Optimal Step Size.
Fig. 4,5,6 is respectively disclosed initial velocity model in embodiment, true velocity Modelling and inversion rate pattern.
Wherein background color is that the picture frame of grey is to adopt GPU to realize the calculating link of accelerating.
Embodiment
Below in conjunction with accompanying drawing, embodiments of the present invention are described in detail, these accompanying drawings have been expressed a kind of specific implementation of the present invention as a part of the present invention.The present invention allows multiple multi-form embodiment as can be known on the basis of understanding accompanying drawing of the present invention and preferred embodiment and connotation of the present invention, therefore in not departing from the scope of the present invention, also can there be other embodiment to be used and to build, and may makes other optional change.
As shown in Figure 1, hybrid domain Full wave shape inversion method of the present invention has and comprises:
The first step is chosen optimum inverting frequency, and is read initial model, inverted parameters and other parameters.
Second step carries out frequency (group) circulation according to the inverting frequency that sets; Iterative loop in frequency (group) circulation.
The 3rd step entered shot point circulation, and reads single big gun data, was specified to the picture pore diameter range, and with velocity field again the CPU internal memory be transferred in the GPU video memory.
The 4th step, as shown in Fig. 2 a, with the time domain propagation operator, focus is carried out main story, record the earth's surface and accept record, utilize simultaneously DFT to extract the frequency field main story wave field of corresponding inverting frequency.Each wave field snapshot is constantly calculated this moment wave field to the contribution of frequency wave field, and add up it.Because each spatial point is separate, so being fit to the GPU fine grained parallel very much, the realization of DFT calculates.
The 5th step, in time domain, record and actual big gun are accepted in the earth's surface record and do poorly, calculate simultaneously target function value corresponding to this inverting frequency.
The 6th step, as shown in Fig. 2 b, take the residual error wave field as focus, utilize the time domain propagation operator to carry out anti-pass to it, utilize equally DFT to extract the frequency field wave field of corresponding inverting frequency, all time domain wave fields contributions to the frequency wave field constantly add up.
In the 7th step, utilize formula:
Calculating is worked as forward gun and is recorded in corresponding gradient fields under current frequency (group), and asks for the maximal value of gradient fields.
The 8th step repeated for the 3rd step to the 7th step to next big gun collection, and the gradient fields of this common-shot-gather under current frequency that add up simultaneously is until all shot point circulations are complete.
In the 9th step, utilize 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, by the value of iterative search α 1 and α 2, its target function value F (α) satisfies the condition of F (a1)＜F (α 0) and F (α 2)＞F (α 1) afterwards; Further, construct a para-curve with these three groups of points, and ask for minimum point, with this as Optimal Step Size.
In the tenth step, come the renewal speed field with Optimal Step Size and gradient fields, and this velocity field is continued to leave in the GPU video memory, be used for next iteration, until satisfy stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency (group), namely repeat ten steps of second step to the.
The 11 step, complete when the equal inverting of all frequencies, the velocity field of final updated is transferred to CPU by the GPU video memory, and in write store.At this moment, inverting is complete.
As shown in Fig. 4-Fig. 6, the present embodiment adopts the alphabetical model of CAS (Chinese Academy of Sciences writes a Chinese character in simplified form) 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 that the shot point number is 65 big guns, every big gun 400 roads, and record length is 1.6s, and the time sampling interval is 1ms, and shot interval is 30m.
Inverted parameters is set to the inverting frequency and chooses 23 frequencies take 1.5Hz as the interval from 2Hz to 35Hz, and inverting circulates according to group of frequencies, and each group of frequencies is chosen 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) are respectively the target function values of the k time and the k+1 time.
It is the CPU of i7-2600 that the computer that is used for inverting in the present embodiment has adopted model, and model is the GPU of GTX680.Can find out by accompanying drawing, adopt method of the present invention, inversion result is portrayed the true velocity model well.
concrete implementation is as follows: at first big gun data are read into calculator memory from hard disk, and determine the corresponding pore size of this big gun, then with the speed this aperture in and work as the forward gun data from the CPU memory copying to the GPU video memory, call successively afterwards tstep_kernel (high-order limited difference GPU executive routine in CPU, being used for carrying out ripple in time domain propagates), pml_kernel (PML absorbing boundary GPU executive routine, be used for absorbing the reflection from the border), (the GPU program that wave field extracts is accepted on the earth's surface to record_kernel, be used for extracting corresponding earth's surface from the wave field snapshot and accept wave field), (discrete Fourier changes the GPU executive routine to DFT_kernel, be used for asking for current time time domain wave field to the contribution of current frequency field wave field), residual_kernel (the GPU program that the time domain residual error is asked for, be used for asking for the difference of prediction time domain wave field and true wave field), obj_kernel (the GPU program that target function value is asked for, be used for calculating the target function value under the present speed model), addresi_kernel (the GPU executive routine that the residual error wave field loads, be used for when the anti-pass of residual error wave field, the residual error wave field is loaded with the focus form), tstep_kernel and pml_kernel (being used for propagating the residual error wave field), DFT2_kernel (being used for calculating current time domain wave field to the contribution of current frequency field wave field), grad_kernel (the GPU program that gradient fields is asked for, be used 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, comprise a plurality of kernel program modules in this module, be used for asking for Optimal Step Size with the method for Parabolic Fit), updatevel_kernel (the GPU program of renewal speed model, utilize gradient fields and Optimal Step Size to come the renewal speed field).After all GPU routine calls were complete, a complete iterative process just was through with, and enters into afterwards next iteration, until satisfy stopping criterion for iteration or reach maximum iteration time.Iteration namely enters next inverting group of frequencies, until all inverting frequencies are finished after finishing.
In this process, GPU has born most calculated amount, thereby counting yield is significantly promoted.
In above-mentioned concrete enforcement, used a plurality of function names, as pml_kernel etc., all can realize by the function that calls and revise in the CUDA function library that NVIDIA company (being the GTX680 developer) openly provides.
And the experiment that above-mentioned CAS model carries out on the CPU of simple i7-2600 is shown, CPU can't complete this inverting flow process within the time identical with the inverting of above-mentioned collaborative calculating execution Full wave shape, can not get result; Even in fact adopt above-mentioned same algorithm to use merely CPU to carry out the Full wave shape inverting, inverting can't be completed in the time can accept, 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 241s consuming time only this shows that method of the present invention has significantly promoted counting yield, has good speed-up ratio.
Above-mentionedly comprise that the embodiment of accompanying drawing is only example of the present invention, and only be to provide 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 disclosed content in claim and instructions, and is not limited to given embodiment.
Claims (8)
1.CPU/GPU the hybrid domain Full wave shape inversion method of collaborative parallel computation is characterized in that described method comprises the steps: 1, and just drilling of waveform inversion partly is put into time domain; 2, adopt discrete Fourier transformation to extract the frequency field wave field, and forward frequency field to and carry out inverting; Wherein, above-mentioned steps 1, the 2nd is completed calculating in GPU.
2. hybrid domain Full wave shape inversion method according to claim 1 is characterized in that described step specifically comprises: the first step, and choose optimum inverting frequency, and read initial model, inverted parameters and other inverting desired parameters; Second step carries out group of frequencies circulation or single-frequency point circulation according to the inverting frequency that sets, and recycle design adopts iterative loop;
The 3rd step entered the shot point circulation, read single big gun data, determines the imaging pore diameter range of single big gun, and corresponding velocity field is transferred in the GPU video memory by internal memory;
The 4th step, with Spatial higher order time domain propagation operator, focus is carried out main story, eliminate simultaneously boundary echo; In each propagation constantly of main story, record the earth's surface and accept wave field, utilize simultaneously DFT to extract the frequency field main story wave field of corresponding inverting frequency, namely each wave field snapshot is constantly calculated this moment wave field to the contribution of current frequency field wave field, these contributions are added up;
The 5th step, in time domain, wave field and actual big gun record is accepted on the earth's surface done poorly, calculate simultaneously target function value corresponding to this inverting frequency;
The 6th step, as focus, utilize the time domain propagation operator to carry out anti-pass to it with the residual error wave field, utilize DFT to extract the frequency field wave field of corresponding inverting frequency, all time domain wave fields contributions to current frequency field wave field constantly namely add up;
In the 7th step, calculate the gradient fields that is recorded in correspondence under current frequency or group of frequencies when forward gun;
The 8th step to next big gun collection, repeated for the 3rd step to the 7th step, and the cumulative gradient fields of this common-shot-gather under current frequency, until all shot point circulations are complete, then ask for the maximal value in this gradient fields simultaneously, is used for choosing of Optimal Step Size;
In the 9th step, utilize the Parabolic Fit method to ask for Optimal Step Size;
In the tenth step, with Optimal Step Size and gradient fields renewal speed field, this velocity field is continued to leave in the GPU video memory, be used for next iteration, repeat above-mentioned steps two to nine, until satisfy stopping criterion for iteration or reach maximum iteration time, then carry out the inverting of next frequency or group of frequencies;
The 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 the GPU video memory, and in write store, completes refutation process;
Wherein, the process in described the 4th step to the tenth step is walked abreast by GPU and accelerates computing.
3. hybrid domain Full wave shape inversion method according to claim 2, is characterized in that described the 4th step employing optimum matching layer, and namely the PML boundary condition is eliminated boundary echo.
4. hybrid domain Full wave shape inversion method according to claim 1, is characterized in that the formula of asking for gradient fields in described the 7th step is
Wherein, P
fBe focus main story wave field, P
bBe residual error anti-pass wave field, v is the speed of underground medium, x
sAnd x
rBe respectively the horizontal level of focus and wave detector, ω is the frequency of current inverting, and * represents variable is got conjugation.
5. the hybrid domain Full wave shape method of inversion according to claim 1, it is characterized in that described the 9th step utilizes following method to build para-curve, and ask for minimum point: with 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 α
0The target function value of=0 o'clock is afterwards by iterative search α
1And α
2Value, its target function value F (α) satisfies 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, with its minimum point as Optimal Step Size.
6. the described hybrid domain Full wave shape of any one method of inversion according to claim 1-5, it is characterized in that described the 4th step also comprises the steps: in the circulation of the time of wave field main story, utilize the DFT operator to ask for each moment to the contribution of frequency field wave field, until that all moment all propagate is complete, obtain frequency field main story wave field corresponding to current inverting frequency or group of frequencies.
7. the described hybrid domain Full wave shape of any one method of inversion according to claim 1-5, it is characterized in that described the 6th step also comprises the steps: with the residual error wave field as the focus anti-pass, in this time circulation, utilize the DFT operator to ask for each residual error anti-pass wave field contribution to the frequency field wave field constantly, until that all moment all propagate is complete, obtain frequency field anti-pass wave field corresponding to current inverting frequency or group of frequencies.
8. the hybrid domain Full wave shape method of inversion according to claim 1, it is characterized in that when the collaborative parallel computation of CPU/GPU, collaborative account form both is to call GPU in CPU to carry out Spatial higher order finite difference GPU module, PML absorbing boundary GPU module, the GPU module that wave field extracts is accepted on the earth's surface, discrete Fourier during the wave field main story changes the GPU module, the GPU module that the time domain residual error is asked for, the GPU module that target function value is asked for, the GPU module that the residual error wave field loads, discrete Fourier transformation GPU module during the anti-pass of residual error wave field, the GPU module that gradient fields is asked for, ask for the peaked GPU module of gradient fields, the GPU module that Optimal Step Size is asked for, the GPU 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 true CN103135132A (en) | 2013-06-05 |
| CN103135132B 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 (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570090A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Method for extracting full-waveform inversion noise filtering operator and performing noise filtering through full-waveform inversion noise filtering operator |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105005072A (en) * | 2015-06-02 | 2015-10-28 | 中国科学院地质与地球物理研究所 | PML boundary three-dimensional seismic wave propagation simulation method utilizing CUDA |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | Efficient time domain full waveform inversion method |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod 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 |
| CN106294273A (en) * | 2015-06-05 | 2017-01-04 | 中国石油化工股份有限公司 | The communication means of a kind of CPU and coprocessor and system |
| CN106950596A (en) * | 2017-04-11 | 2017-07-14 | 中国石油大学(华东) | A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate |
| CN107656307A (en) * | 2016-07-26 | 2018-02-02 | 中国石油化工股份有限公司 | Full waveform inversion computational methods and system |
| CN107656306A (en) * | 2016-07-26 | 2018-02-02 | 中国石油化工股份有限公司 | Full waveform inversion parallel calculating method and system |
| CN107728199A (en) * | 2017-09-22 | 2018-02-23 | 中国地质大学（北京） | Based on the parallel multi -components anisotropy pre-stack time migration accelerated methods of more GPU |
| CN107918145A (en) * | 2016-10-10 | 2018-04-17 | 中国石油化工股份有限公司 | The parallelization processing method and system of earthquake big gun energy |
| CN108594299A (en) * | 2018-02-28 | 2018-09-28 | 中国科学院地质与地球物理研究所 | High ferro intelligent early-warning method, apparatus and system |
| CN109313619A (en) * | 2016-06-06 | 2019-02-05 | 高通股份有限公司 | Power and the Memory Controller of performance aware ballot mechanism |
| CN109523022A (en) * | 2018-11-13 | 2019-03-26 | Oppo广东移动通信有限公司 | Terminal data processing method, apparatus and terminal |
| CN109633742A (en) * | 2019-01-08 | 2019-04-16 | 中国科学院地质与地球物理研究所 | A kind of full waveform inversion method and device |
| CN110888158A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Full waveform inversion method based on RTM constraint |
| CN110895347A (en) * | 2019-12-09 | 2020-03-20 | 中国科学院地质与地球物理研究所 | Seismic wave forward modeling method and system |
| 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 |
| CN112578431A (en) * | 2019-09-27 | 2021-03-30 | 中国石油化工股份有限公司 | Optimal storage method and system for full waveform inversion wave field in limited state |
| CN113504566A (en) * | 2021-06-01 | 2021-10-15 | 南方海洋科学与工程广东省实验室(湛江) | Seismic inversion method, system, device and medium based on wave equation |
| CN114814843A (en) * | 2022-06-27 | 2022-07-29 | 之江实验室 | 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 |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105891888B (en) * | 2016-03-28 | 2017-03-08 | 吉林大学 | Multiple domain divides multiple dimensioned full waveform inversion method parallel |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040186667A1 (en) * | 2003-03-18 | 2004-09-23 | Lee Ki Ha | Source-independent full waveform inversion of seismic data |
| US20070203673A1 (en) * | 2005-11-04 | 2007-08-30 | Sherrill Francis G | 3d pre-stack full waveform inversion |
| US20070282535A1 (en) * | 2006-05-31 | 2007-12-06 | Bp Corporation North America Inc. | System and method for 3d frequency domain waveform inversion based on 3d time-domain forward modeling |
| 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 |
-
2013
- 2013-01-15 CN CN201310013786.6A patent/CN103135132B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040186667A1 (en) * | 2003-03-18 | 2004-09-23 | Lee Ki Ha | Source-independent full waveform inversion of seismic data |
| US20070203673A1 (en) * | 2005-11-04 | 2007-08-30 | Sherrill Francis G | 3d pre-stack full waveform inversion |
| US20070282535A1 (en) * | 2006-05-31 | 2007-12-06 | Bp Corporation North America Inc. | System and method for 3d frequency domain waveform inversion based on 3d time-domain forward modeling |
| 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》 * |
Cited By (38)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570090A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Method for extracting full-waveform inversion noise filtering operator and performing noise filtering through full-waveform inversion noise filtering operator |
| 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 |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | 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 |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod 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 |
| CN105005072A (en) * | 2015-06-02 | 2015-10-28 | 中国科学院地质与地球物理研究所 | PML boundary three-dimensional seismic wave propagation simulation method utilizing CUDA |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN106294273A (en) * | 2015-06-05 | 2017-01-04 | 中国石油化工股份有限公司 | The communication means of a kind of CPU and coprocessor and system |
| CN106294273B (en) * | 2015-06-05 | 2020-01-10 | 中国石油化工股份有限公司 | Communication method and system of CPU and coprocessor |
| CN109313619A (en) * | 2016-06-06 | 2019-02-05 | 高通股份有限公司 | Power and the Memory Controller of performance aware ballot mechanism |
| CN107656307A (en) * | 2016-07-26 | 2018-02-02 | 中国石油化工股份有限公司 | Full waveform inversion computational methods and system |
| CN107656306A (en) * | 2016-07-26 | 2018-02-02 | 中国石油化工股份有限公司 | Full waveform inversion parallel calculating method and system |
| CN107656306B (en) * | 2016-07-26 | 2019-02-19 | 中国石油化工股份有限公司 | Full waveform inversion parallel calculating method and system |
| CN107918145A (en) * | 2016-10-10 | 2018-04-17 | 中国石油化工股份有限公司 | The parallelization processing method and system of earthquake big gun energy |
| CN106950596A (en) * | 2017-04-11 | 2017-07-14 | 中国石油大学(华东) | A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate |
| 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 |
| CN107728199A (en) * | 2017-09-22 | 2018-02-23 | 中国地质大学（北京） | Based on the parallel multi -components anisotropy pre-stack time migration accelerated methods of more GPU |
| CN108594299A (en) * | 2018-02-28 | 2018-09-28 | 中国科学院地质与地球物理研究所 | High ferro intelligent early-warning method, apparatus and system |
| CN110888158A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Full waveform inversion method based on RTM constraint |
| CN110888158B (en) * | 2018-09-10 | 2021-12-31 | 中国石油化工股份有限公司 | Full waveform inversion method based on RTM constraint |
| CN109523022A (en) * | 2018-11-13 | 2019-03-26 | Oppo广东移动通信有限公司 | Terminal data processing method, apparatus and terminal |
| CN109523022B (en) * | 2018-11-13 | 2022-04-05 | Oppo广东移动通信有限公司 | Terminal data processing method and device and terminal |
| WO2020143645A1 (en) * | 2019-01-08 | 2020-07-16 | 中国科学院地质与地球物理研究所 | Full waveform inversion method and apparatus |
| CN109633742A (en) * | 2019-01-08 | 2019-04-16 | 中国科学院地质与地球物理研究所 | A kind of 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 |
| CN112578431A (en) * | 2019-09-27 | 2021-03-30 | 中国石油化工股份有限公司 | Optimal storage method and system for full waveform inversion wave field in limited state |
| 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 |
| CN113504566A (en) * | 2021-06-01 | 2021-10-15 | 南方海洋科学与工程广东省实验室(湛江) | Seismic inversion method, system, device and medium based on wave equation |
| CN113504566B (en) * | 2021-06-01 | 2024-04-30 | 南方海洋科学与工程广东省实验室(湛江) | Wave equation-based seismic inversion method, system, device and medium |
| CN114814843A (en) * | 2022-06-27 | 2022-07-29 | 之江实验室 | 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 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN103135132B (en) | 2015-07-01 |
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 | |
| CN107783190B (en) | A kind of least square reverse-time migration gradient updating method | |
| CN106842320B (en) | The parallel 3-D seismics wave field generation method of GPU and system | |
| Abdelkhalek et al. | Fast seismic modeling and reverse time migration on a GPU cluster | |
| CN105319581B (en) | A kind of efficient time-domain full waveform inversion method | |
| CN105137486A (en) | Elastic wave reverse-time migration imaging method and apparatus in anisotropic media | |
| CN104570082B (en) | Extraction method for full waveform inversion gradient operator based on green function characterization | |
| CN107894618B (en) | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm | |
| CN103105623B (en) | Data waveform processing method in seismic exploration | |
| WO2023087451A1 (en) | Observation data self-encoding-based multi-scale unsupervised seismic wave velocity inversion method | |
| Poul et al. | Efficient time-domain deconvolution of seismic ground motions using the equivalent-linear method for soil-structure interaction analyses | |
| CN103970960A (en) | Grid-free Galerkin method structural topology optimization method based on GPU parallel acceleration | |
| CN105005072B (en) | The PML border dimensionally seismic wave propagating mode utilizing CUDA intends method | |
| CN108387933A (en) | A kind of method, apparatus and system of definitely interval quality factors | |
| CN105891888A (en) | Multi-domain frequency-division parallel multi-scale full-waveform inversion method | |
| CN107561585A (en) | A kind of multinuclear multi-node parallel 3-D seismics wave field generation method and system | |
| CN104965222B (en) | Three-dimensional longitudinal wave impedance full-waveform inversion method and apparatus | |
| CN105911584B (en) | Implicit staggered-grid finite difference elastic wave numerical simulation method and device | |
| CN110879412A (en) | Underground transverse wave velocity inversion method, device, computing equipment and storage medium | |
| Xue et al. | An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging | |
| CN103698810A (en) | Hybrid network minimum travel time ray tracing tomography method | |
| CN109490948B (en) | Seismic acoustic wave equation vector parallel computing method | |
| US20160034612A1 (en) | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing | |
| CN107526104A (en) | Fracture medium seismic wave field method for numerical simulation based on multimachine multinuclear | |
| CN114002742B (en) | Euler Gaussian beam offset imaging method and device |
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 |