CN104977605B - A kind of high-resolution multiple dimensioned wave equation inversion method - Google Patents

A kind of high-resolution multiple dimensioned wave equation inversion method Download PDF

Info

Publication number
CN104977605B
CN104977605B CN201410129114.6A CN201410129114A CN104977605B CN 104977605 B CN104977605 B CN 104977605B CN 201410129114 A CN201410129114 A CN 201410129114A CN 104977605 B CN104977605 B CN 104977605B
Authority
CN
China
Prior art keywords
mrow
msub
mfrac
msup
rho
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201410129114.6A
Other languages
Chinese (zh)
Other versions
CN104977605A (en
Inventor
宋建勇
秦臻
李劲松
姚逢昌
魏超
马晓宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Natural Gas Co Ltd
Original Assignee
China Petroleum and Natural Gas Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Natural Gas Co Ltd filed Critical China Petroleum and Natural Gas Co Ltd
Priority to CN201410129114.6A priority Critical patent/CN104977605B/en
Publication of CN104977605A publication Critical patent/CN104977605A/en
Application granted granted Critical
Publication of CN104977605B publication Critical patent/CN104977605B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

A kind of high-resolution multiple dimensioned wave equation inversion method, the invention provides a kind of multiple dimensioned wave equation inversion method, and it utilizes earthquake big gun collection, and based on Acoustic Wave-equation, inverting is carried out to given initial volume modulus and density.Refutation process is carried out in different frequency, space scale, and enters row constraint to destination layer using prior model, and the inverting section of bulk modulus and density is obtained by iteration, the well-posedness of inverting is enhanced, improves inverting iterative convergence speed.

Description

A kind of high-resolution multiple dimensioned wave equation inversion method
Technical field
The present invention relates to technical field of geological exploration, is concretely a kind of high-resolution multiple dimensioned wave equation inversion Method.
Background technology
Reservoir properties characteristic parameter controls the spatial distribution of oil gas and the reserves estimation of oil gas, directly affects oil-gas mining The design and recovery ratio of scheme, it is always an emphasis of oil-gas exploration and development research.With oilfield prospecting developing at this stage The increase of difficulty and the development of computing technique, in recent years, reservoir prediction and fluid detection geophysical method and the research of technology Make remarkable progress.Seismic inversion is one of reservoir prediction and fluid detection technology series.
Seismic inversion is the geological data using field inspection, and geologic objective is asked for by mathematical method and computing technique Physical parameter.Round trip Acoustic Wave-equation can simulate the seismic waveform in acoustic medium, with traditional inverting（Wave resistance anti-reflective Drill, elastic impedance inverting and AVO invertings）Compare, the equation can handle arbitrarily complicated non-uniform dielectric situation, improve inverting essence True property.
Acoustic Wave-equation inverting is faced with two large problems：First, it is bottleneck that actual industrial production faces that computational efficiency is low Sex chromosome mosaicism.It is the convergence speed for improving the computational efficiency for calculating wave equation simulation or improving inverting to solve this key to the issue Degree.Tarantola etc. (1988) m- Simulation of depth seismic wave field, Pratt etc. when proposing（1998）The frequency space domain of proposition Simulation single-frequency wave field.Both approaches reduce the time required for an inverting iteration.In addition, efficient parallel computing Calculating speed greatly improved.2nd, resolution of inversion is relatively low prevents wave equation inversion from being applied to reservoir inversion very well.Earthquake The dominant frequency of data is mainly 30Hz-50Hz, if only utilizing seismic data inversion, final inversion result resolution ratio is also than relatively low.
The content of the invention
In order to solve problem of the prior art, the efficiency of seismic data inversion is improved, there is provided a kind of high-resolution more Yardstick wave equation inversion method, can Fast Convergent, improve the efficiency of inverting.
The embodiments of the invention provide a kind of high-resolution multiple dimensioned wave equation inversion method, including：
Step 101, shot gather data is gathered；
Step 102, stratigraphic horizon data H (x) is inputted；
Step 103, restricted model m is establishedp(x, z)=(Kp,ρp), wherein, the model of input is the Reservoir Section of required inverting Bulk modulus KpAnd density pp；
Step 104, it is a 2-D data P to remember each shot gather datao(t, x), wherein, x represents distance, and t represents sampling Time point, if big gun integrates number as ns, to shot gather data Po(t, x) carries out multi-resolution decomposition, obtains the ground under different dominant frequency yardsticks Shake big gun collection Po(t,x,ωi), ωi=i Δ ω, multi-resolution decomposition formula are represented by：
Wherein ω represents the centre frequency of Gaussian function；σ represents the width of Gaussian function, and multiple dimensioned lower earthquake big gun collection can table It is shown as
Po(t,x,ωi)=FFT-1[FFT(P(t,x))·G(ωi)]
To the shot gather data P of each yardsticko(t,x,ωi), following step 105 is repeated to step 111；
Step 105, if the input m of initial model0(x, z)=(K0,ρ0), K0For the initial volume of required inverting Reservoir Section Modulus, ρ0For the initial density of required inverting Reservoir Section；
Step 106, for frequencies omegai, using dominant frequency as ωiZero phase Ricker wavelet, according to Acoustic Wave-equation and absorption Boundary condition calculating simulation shot gather data Ps(t,x,ωi)；
Step 107, the earthquake shot gather data P received is calculatedo(t,x,ωi) with simulating shot gather data Ps(t,x,ωi) Residual error d (t, x), is designated as
D (t, x)=Po-Ps（4）
Wherein, PoFor Po(t,x,ωi), PsFor Ps(t,x,ωi)；
Step 108, with-d (- t, x) for focus, the focus item of Acoustic Wave-equation is replaced, simulation residual error records reverse Shot gather data Ps'；
Step 109, the ladder changed according to the above-mentioned shot gather data being calculated and reverse shot gather data construction inverse model Degree：
E (p) is shot gather data energy error, above-mentioned formula（5）The conjugation modification amount of inverse model parameter be represented by：
Wherein, x is the position of receiving point, and t is time variable, xsFor the position of focus, ▽ represents gradient operator, and T is big gun Collection record time span, K, ρ are inverse model parameter,For inverse model parameter modification amount, iωRepresent i-th Individual frequency, jsRepresent jth big gun；Concentrate all big gun repeat steps 106 whole to step 109, accumulation calculating in the big gun under the frequency Model modification amount；
Step 110, all model modification amounts calculated using step 109, generate new inverse model,
The new inverse model is：
mn+1=mn-αngn（8）
Wherein, m=(K, ρ) is inverse model, and n is iterations, αnFor nth iteration step-length, gn=(δ K, δ ρ) is n-th The inverse model variable gradient of secondary iteration；
Step 111, energy error is calculated according to the following formula,
E is energy error in formula, Δ dtFor Δ d transposition, Δ d is the residual error in above-mentioned steps 108, CDFor data association side Poor matrix；
As E ＜ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is anti-as next time in this frequency range The initial model drilled, into repeat step 105；Wherein iωFor i-th of frequency, nωFor frequency sampling number；
As E ＞ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is as next frequency range inverting Initial model, into repeat step 104；
As E ＜ δ, iω=nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is as final inverting knot Fruit, wherein, δ is arbitrarily small real number.
One of a kind of high-resolution multiple dimensioned wave equation inversion method according to embodiments of the present invention is further Aspect, denoising is carried out to the shot gather data collected in the step 101, removes environmental noise and system noise.
Another of a kind of high-resolution multiple dimensioned wave equation inversion method according to embodiments of the present invention enters one The aspect of step, in step 108, the Acoustic Wave-equation is：
Wherein, s (x0,z0, t) and it is focus item, P is stress vector, vxAnd vzRespectively simulate seismic profile Ps(t,x, ωi), timely m- space wave field P (t, x, z).
Another of a kind of high-resolution multiple dimensioned wave equation inversion method according to embodiments of the present invention enters one The aspect of step, in the step 108, using improved nondividing PML absorbing boundary conditions, i.e.,
Wherein Ωxx、Ωzz、Ψxx、Ψzz、Txx、TzzFor the middle auxiliary variable of introducing, subsidiary equation corresponding to them is：
Above-mentioned equation group（2）With（3）The equation of absorbing boundary condition is formed,
The equation of simultaneous Acoustic Wave-equation and absorbing boundary condition obtains big gun collection wave field Ps。
Pass through the Multi-scale inversion of the embodiment of the present invention so that wave equation inversion is stable, convergence is quick.Theoretical model is surveyed Test result indicates the validity of this method.Exemplary application is China's Sichuan Gas Fields real data inverting, by 5 iteration Afterwards, it is preferable to obtain result.
Brief description of the drawings
Read the detailed description to embodiment in conjunction with the following drawings, features described above and advantage of the invention, and extra Feature and advantage, it will clearer.
Fig. 1 show a kind of flow chart of high-resolution multiple dimensioned wave equation inversion method of the embodiment of the present invention；
Fig. 2 show the schematic diagram of the big gun in the big gun of field of embodiment of the present invention shot gather data 1000；
Fig. 3 a- Fig. 3 d are the schematic diagram of 4 different frequency ranges in geological data of the embodiment of the present invention；
Fig. 4 a are the schematic diagram of speed initial model of the embodiment of the present invention；
Fig. 4 b are the schematic diagram of density initial model of the embodiment of the present invention；
Fig. 5 a to Fig. 5 d are the increment schematic diagram for the bulk modulus that the embodiment of the present invention obtains according to four different frequencies；
Fig. 6 a to Fig. 6 d are the bulk modulus inversion result schematic diagram that four different frequencies of the embodiment of the present invention obtain；
Fig. 7 show each frequency range of the embodiment of the present invention and passes through the final inversion result schematic diagram that 5 iteration obtain.
Embodiment
Following description can make any those skilled in the art using the present invention.Provided in specific embodiment and application Description information it is merely illustrative.The various extensions and combination of embodiment as described herein are aobvious for those skilled in the art And be clear to, without departing from the spirit and scope of the present invention, the rule that the present invention defines may apply to it In his embodiment and application.Therefore, the present invention is not limited solely to shown embodiment, and the present invention covers and principle illustrated herein and spy Levy consistent maximum magnitude.
It is as shown in Figure 1 a kind of flow chart of high-resolution multiple dimensioned wave equation inversion method of the embodiment of the present invention.
Including step 101, shot gather data is gathered, and denoising is carried out to the shot gather data, environmental noise is removed and is System noise.
Step 102, stratigraphic horizon data H (x) is inputted.In the inventive solutions with inverting of the prior art not Together, the data inputted in the present invention are depth numeric field data.
Step 103, restricted model m is establishedp(x, z)=(Kp,ρp), wherein, the model of input is the Reservoir Section of required inverting Bulk modulus KpAnd density pp；The model be space in parameter, full waveform inversion in the prior art（Full Waveform Inversion）In do not need prior model.
Step 104, it is a 2-D data P to remember each shot gather datao(t, x), wherein, x represents distance, and t represents sampling Time point, if big gun integrates number as ns。
To big gun collection Po(t, x) carries out multi-resolution decomposition, obtains the earthquake big gun collection P under different dominant frequency yardstickso(t,x,ωi), ωi=i Δs ω.Multi-resolution decomposition formula is represented by：
Wherein ω represents the centre frequency of Gaussian function；σ represents the width of Gaussian function, and σ is bigger, and Gaussian function is wider. Multiple dimensioned lower earthquake big gun collection is represented by
Po(t,x,ωi)=FFT-1[FFT(P(t,x))·G(ωi)]
To the shot gather data P of each yardsticko(t,x,ωi), following step 105 is repeated to step 111；
Step 105, if the input m of initial model0(x, z)=(K0,ρ0), K0For the initial volume of required inverting Reservoir Section Modulus, ρ0For the initial density of required inverting Reservoir Section；The two models are the coarse smooth of the condition foundation known to Depth domain model.If inversion technique just updates this model by different iterationses.
Step 106, for frequencies omegai, using dominant frequency as ωiZero phase Ricker wavelet, according to Acoustic Wave-equation and absorption Boundary condition calculating simulation shot gather data Ps(t,x,ωi)。
Step 107, the earthquake shot gather data P received is calculatedo(t,x,ωi) with simulating shot gather data Ps(t,x,ωi) Residual error d (t, x), is designated as
D (t, x)=Po-Ps
Wherein, PoFor Po(t,x,ωi), PsFor Ps(t,x,ωi)。
Step 108, with-d (- t, x) for focus, the focus item of Acoustic Wave-equation is replaced, simulation residual error records reverse Shot gather data Ps'。
Step 109, the ladder changed according to the above-mentioned shot gather data being calculated and reverse shot gather data construction inverse model Degree：
E (p) is shot gather data energy error, above-mentioned formula（5）The conjugation modification amount of inverse model parameter be represented by：
Wherein, x is the position of receiving point, and t is time variable, xsFor the position of focus, ▽ represents gradient operator, and T is big gun Collection record time span, K, ρ are inverse model parameter,For inverse model parameter modification amount, iωRepresent i-th Individual frequency, jsRepresent jth big gun.
All big gun repeat steps 106 are concentrated to big gun to step 109, the whole model modification amount of accumulation calculating；
Step 110, all model modification amounts calculated using step 109, generate new inverse model,
The new inverse model is：
mn+1=mn-αngn（8）
Wherein, m=(K, ρ) is inverse model, and n is iterations, αnFor nth iteration step-length, gn=(δ K, δ ρ) is n-th The inverse model variable gradient of secondary iteration；
Step 111, energy error is calculated according to the following formula,
E is energy error in formula, Δ dtFor Δ d transposition, Δ d is the residual error in above-mentioned steps 108, CDFor data association side Poor matrix；
As E ＜ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is anti-as next time in this frequency range The initial model drilled, into repeat step 105；Wherein, iωFor i-th of frequency, nωFor frequency sampling number.
As E ＞ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is as next frequency range inverting Initial model, into repeat step 104；
As E ＜ δ, iω=nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is as final inverting knot Fruit, wherein, δ is arbitrarily small real number, such as can be 10-6Or 10-5Predetermined value.
As one embodiment of the present of invention, the Acoustic Wave-equation in the step 108 can be：
Wherein, s (x0,z0, t) and it is focus item, P is stress vector, vxAnd vzRespectively simulate seismic profile Ps(t,x, ωi), timely m- space wave field P (t, x, z).
Using improved nondividing PML absorbing boundary conditions, i.e.,
Wherein Ωxx、Ωzz、Ψxx、Ψzz、Txx、TzzFor the middle auxiliary variable of introducing, can all be filled in the starting stage Zero, subsidiary equations corresponding to them are：
More than simultaneous two formula（1）、（2）With（3）The solution of PML wave equation problems is can obtain, obtains solving Ps, as big gun collection Wave field Ps。
It is quick by Multi-scale inversion wave equation inversion stabilization, convergence by embodiments of the invention.Theoretical mould Type test result indicates the validity of this method.
The schematic diagram of a big gun in the big gun of field of embodiment of the present invention shot gather data 1000 is illustrated in figure 2, as the defeated of inverting Enter；
Step 201, actual seismic shot gather data is gathered（Fig. 2）, static correction, terrestrial reference uniformity amplitude are carried out to shot gather data Compensation and prestack denoising, obtain the higher earthquake shot gather data of signal to noise ratio, the input as inverting.
Step 202, stratigraphic horizon H (x) is inputted.The layer position of this example is the layer position of Depth Domain, and unconventional prestack inversion Required Depth Domain layer position.
Step 203, restricted model m is establishedp(x, z)=(Kp,ρp), in this example, utilize well-log information and layer position interpolation of data Go out prior model（Bulk modulus Kp, density pp）, as constraint：Included in prior model known as the true model of underground Information, the data of each inverting are controlled without departing from the numerical value certain limit of prior model with these information, so that instead Drill result and obtain more accurately result.
Step 204, multi-resolution decomposition is carried out according to Gaussian function to geological data, specific decomposable process may be referred to Fig. 1 Middle embodiment corresponding steps, decompose 4 different frequency ranges（Yardstick）（Fig. 3 a- Fig. 3 d）, each frequency range shot gather data is repeated Step 205 is to step 211.
Step 205, more accurately macroscopic velocity initial velocity model is obtained with tomographic inversion, as shown in fig. 4 a, wherein, Tomographic inversion is a kind of conventional velocity modeling method, is not repeated herein, and it is initial to obtain density by GAEDENAR formula Model, as shown in Figure 4 b；By model and density, initial volume modulus model is calculated.
Step 206, for each frequencies omegai, using dominant frequency as ωiZero phase Ricker wavelet, it is public according to Acoustic Wave-equation Formula（1）And the equation formulations of absorbing boundary condition（2）With（3）Calculating simulation shot gather data Ps(t,x,ωi)。
Step 207, the earthquake shot gather data P received is calculatedo(t,x,ωi) with simulating shot gather data Ps(t,x,ωi) Residual error d (t, x).
Step 208, with-d (- t, x) for focus, the focus item of Acoustic Wave-equation is replaced, simulation residual error records reverse Shot gather data Ps'；
Step 209, the ladder changed according to the above-mentioned shot gather data being calculated and reverse shot gather data construction inverse model DegreeAll big gun repeat steps 206 are concentrated to big gun to step 209, the whole model modification amount of accumulation calculating.Fig. 5 a The increment of the bulk modulus obtained to Fig. 5 d four different frequencies of expression.
Step 210, all model modification amounts calculated using step 209, generate new inverse model（Fig. 6 a- Fig. 6 d）. Fig. 6 a to Fig. 6 d represent the bulk modulus inversion result that four different frequencies obtain.
Step 211, with four frequency bands, each frequency band iteration five times（Other integer can also be set）, such as Fig. 7 institutes Show.
It is initial to what is given based on Acoustic Wave-equation using earthquake big gun collection by the method for the embodiments of the present invention Bulk modulus and density carry out inverting.Refutation process is carried out in different frequency, space scale, and using prior model to target Layer enters row constraint, and the inverting section of bulk modulus and density is obtained by iteration, the well-posedness of inverting is enhanced, improves inverting Iterative convergence speed.The present invention passes through Multi-scale inversion so that wave equation inversion is stable, convergence is quick.Theoretical model is tested As a result the validity of this method is indicated.Exemplary application is China's Sichuan Gas Fields real data inverting, after 5 iteration, It is preferable to obtain result.
The present invention can realize in any suitable form, including hardware, software, firmware or their any combination.This Invention can according to circumstances selectively part be realized, for example, software performing in one or more data processors and Digital signal processor.The element and component of this paper each embodiment can be fitted with any physically, functionally, in logic When mode realize.In fact, One function can in separate unit, in one group of unit or be used as other functional units A part realize.Therefore, the system and method can both be realized in separate unit, can also be physically and functionally It is distributed between different unit and processor.
Technical staff in the related art will recognize that embodiments of the invention have many possible modifications and group Close, although form is slightly different, still using identical fundamental mechanism and method.For purposes of explanation, it is described above to reference to Several specific embodiments.However, above-mentioned illustrative discussion is not intended to the precise forms that exhaustive or limitation is invented herein.Before Shown in text, many modifications and variations are possible.Selected and described embodiment, to explain the principle and in fact of the present invention Border is applied, to enable repairing for application-specific of the those skilled in the art best using of the invention and each embodiment Change, deform.

Claims (4)

  1. A kind of 1. high-resolution multiple dimensioned wave equation inversion method, it is characterised in that including：
    Step 101, shot gather data is gathered；
    Step 102, stratigraphic horizon data H (x) is inputted；
    Step 103, restricted model m is establishedp(x, z)=(Kp,ρp), wherein, the model of input is the body of the Reservoir Section of required inverting Product module amount KpAnd density pp；
    Step 104, it is a 2-D data P to remember each shot gather datao(t, x), wherein, x and z represent space coordinates, when t is represented Between variable, if big gun integrates number as ns, to shot gather data Po(t, x) carries out multi-resolution decomposition, obtains the ground under different dominant frequency yardsticks Shake big gun collection Po(t,x,ωi), ωi=i Δ ω, multi-resolution decomposition formula are represented by：
    <mrow> <mi>G</mi> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msqrt> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> </msqrt> <mi>&amp;sigma;</mi> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mfrac> <msup> <mi>&amp;omega;</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> </mrow> </msup> </mrow>
    Wherein, i represents imaginary number, and ω represents the centre frequency of Gaussian function, and Δ ω represents frequency sampling rate；σ represents Gaussian function Width, multiple dimensioned lower earthquake big gun collection is represented by
    Po(t,x,ωi)=FFT-1[FFT(P(t,x))·G(ωi)]
    To the shot gather data P of each yardsticko(t,x,ωi), following step 105 is repeated to step 111；
    Step 105, if the input m of initial model0(x, z)=(K0,ρ0), K0For the initial volume modulus of required inverting Reservoir Section, ρ0For the initial density of required inverting Reservoir Section；
    Step 106, for frequencies omegai, using dominant frequency as ωiZero phase Ricker wavelet, according to Acoustic Wave-equation and absorbing boundary Condition calculating simulation shot gather data Ps(t,x,ωi)；
    Step 107, the earthquake shot gather data P received is calculatedo(t,x,ωi) with simulating shot gather data Ps(t,x,ωi) residual error d (t, x), it is designated as
    D (t, x)=Po-Ps (4)
    Wherein, PoFor Po(t,x,ωi), PsFor Ps(t,x,ωi)；
    Step 108, with-d (- t, x) for focus, the focus item of Acoustic Wave-equation, the reverse big gun collection of simulation residual error record are replaced Data Ps'；
    Step 109, the gradient changed according to the above-mentioned shot gather data being calculated and reverse shot gather data construction inverse model：
    <mrow> <mi>g</mi> <mo>=</mo> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mi>E</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&lt;</mo> <mi>&amp;delta;</mi> <mi>&amp;rho;</mi> <mo>,</mo> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;CenterDot;</mo> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <msup> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;prime;</mo> </msup> <mo>&gt;</mo> <mo>+</mo> <mo>&lt;</mo> <mi>&amp;delta;</mi> <mi>K</mi> <mo>,</mo> <mo>&amp;dtri;</mo> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;CenterDot;</mo> <mo>&amp;dtri;</mo> <msup> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;prime;</mo> </msup> <mo>&gt;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
    E (p) is shot gather data energy error, and the conjugation modification amount of the inverse model parameter of above-mentioned formula (5) is represented by：
    <mrow> <msub> <mi>&amp;delta;&amp;rho;</mi> <mrow> <msub> <mi>i</mi> <mi>&amp;omega;</mi> </msub> <mo>,</mo> <msub> <mi>j</mi> <mi>s</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msup> <mi>&amp;rho;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msubsup> <mo>&amp;Integral;</mo> <mn>0</mn> <mi>T</mi> </msubsup> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <msub> <mi>P</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>;</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <msup> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>;</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>d</mi> <mi>t</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
    <mrow> <msub> <mi>&amp;delta;K</mi> <mrow> <msub> <mi>i</mi> <mi>&amp;omega;</mi> </msub> <mo>,</mo> <msub> <mi>j</mi> <mi>s</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msup> <mi>K</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msubsup> <mo>&amp;Integral;</mo> <mn>0</mn> <mi>T</mi> </msubsup> <mo>&amp;dtri;</mo> <msub> <mi>P</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>;</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mo>&amp;dtri;</mo> <msup> <msub> <mi>P</mi> <mi>s</mi> </msub> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>;</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>d</mi> <mi>t</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
    Wherein, r is the position of receiving point, and t is time variable, xsFor the position of focus,Gradient operator is represented, T records for big gun collection Time span, K, ρ are inverse model parameter,For inverse model parameter modification amount, iωI-th of frequency is represented, jsRepresent jth big gun；All big gun repeat steps 106 are concentrated to the big gun under the frequency to step 109, the whole model modification of accumulation calculating Amount；
    Step 110, all model modification amounts calculated using step 109, generate new inverse model,
    The new inverse model is：
    mn+1=mn-αngn (8)
    Wherein, m=(K, ρ) is inverse model, and n is iterations, αnFor nth iteration step-length, gn=(δ K, δ ρ) changes for n-th The inverse model variable gradient in generation；
    Step 111, energy error is calculated according to the following formula,
    <mrow> <mi>E</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mi>&amp;Delta;d</mi> <mi>t</mi> </msup> <msubsup> <mi>C</mi> <mi>D</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mi>&amp;Delta;</mi> <mi>d</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
    E is energy error in formula, Δ dtFor Δ d transposition, Δ d is the residual error in above-mentioned steps 107, CDFor data covariance square Battle array；
    As E ＜ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is as next inverting in this frequency range Initial model, into repeat step 105；Wherein iωFor i-th of frequency, nωFor frequency sampling number；
    As E ＞ δ, iω＜ nωWhen, inversion result is exported, i.e. inverse model new m=(K, ρ) is initial as next frequency range inverting Model, into repeat step 104；
    As E ＜ δ, iω=nωWhen, export inversion result, i.e. inverse model new m=(K, ρ) as final inversion result, its In, δ is arbitrarily small real number.
  2. 2. the high-resolution multiple dimensioned wave equation inversion method of one kind according to claim 1, it is characterised in that described Denoising is carried out to the shot gather data collected in step 101, removes environmental noise and system noise.
  3. 3. the high-resolution multiple dimensioned wave equation inversion method of one kind according to claim 1, it is characterised in that step In 108, the Acoustic Wave-equation is：
    <mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;rho;</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>&amp;rho;</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>+</mo> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>,</mo> <msub> <mi>z</mi> <mn>0</mn> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
    Wherein, s (x0,z0, t) and it is focus item, P is stress vector, vxAnd vzRespectively speed v is the component on x, z direction, according to Formula (1) simulation seismic profile Ps(t,x,ωi), timely m- space wave field P (t, x, z).
  4. 4. the high-resolution multiple dimensioned wave equation inversion method of one kind according to claim 1, it is characterised in that described In step 108, using improved nondividing PML absorbing boundary conditions, i.e.,
    <mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;rho;</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>-</mo> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>&amp;rho;</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>-</mo> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>P</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>-</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> <mo>-</mo> <msup> <mi>&amp;rho;&amp;upsi;</mi> <mn>2</mn> </msup> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
    Wherein Ωxx、Ωzz、Ψxx、Ψzz、Txx、TzzFor the middle auxiliary variable of introducing, subsidiary equation corresponding to them is：
    <mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>&amp;rho;</mi> <mi>x</mi> </msub> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;rho;</mi> <mi>x</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>T</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>&amp;rho;</mi> <mi>z</mi> </msub> <msub> <mi>&amp;Omega;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;rho;</mi> <mi>z</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>T</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>&amp;rho;</mi> <mi>x</mi> </msub> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;rho;</mi> <mi>x</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>V</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>&amp;rho;</mi> <mi>z</mi> </msub> <msub> <mi>&amp;Psi;</mi> <mrow> <mi>z</mi> <mi>z</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;rho;</mi> <mi>z</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>V</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
    Above-mentioned equation group (2) and (3) form the equation of absorbing boundary condition, wherein, ρx、ρzX is represented, the component on z directions；
    The equation of simultaneous Acoustic Wave-equation and absorbing boundary condition obtains big gun collection wave field Ps。
CN201410129114.6A 2014-04-01 2014-04-01 A kind of high-resolution multiple dimensioned wave equation inversion method Active CN104977605B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410129114.6A CN104977605B (en) 2014-04-01 2014-04-01 A kind of high-resolution multiple dimensioned wave equation inversion method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410129114.6A CN104977605B (en) 2014-04-01 2014-04-01 A kind of high-resolution multiple dimensioned wave equation inversion method

Publications (2)

Publication Number Publication Date
CN104977605A CN104977605A (en) 2015-10-14
CN104977605B true CN104977605B (en) 2018-01-02

Family

ID=54274280

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410129114.6A Active CN104977605B (en) 2014-04-01 2014-04-01 A kind of high-resolution multiple dimensioned wave equation inversion method

Country Status (1)

Country Link
CN (1) CN104977605B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101369024A (en) * 2008-09-09 2009-02-18 中国石油天然气集团公司 Earthquake wave equation generation method and system
CN101630018A (en) * 2008-07-16 2010-01-20 中国石油天然气股份有限公司 Method for processing seismic exploration data in process of controlling full acoustic wave equation inversion
CN102385066A (en) * 2010-09-06 2012-03-21 中国石油天然气股份有限公司 Pre-stack seismic quantitative imaging method
CN103592685A (en) * 2013-10-22 2014-02-19 中国石油天然气股份有限公司 Method and device for removing wave equation simulation direct wave from full-waveform inversion

Patent Citations (4)

* Cited by examiner, † Cited by third party
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
CN101369024A (en) * 2008-09-09 2009-02-18 中国石油天然气集团公司 Earthquake wave equation generation method and system
CN102385066A (en) * 2010-09-06 2012-03-21 中国石油天然气股份有限公司 Pre-stack seismic quantitative imaging method
CN103592685A (en) * 2013-10-22 2014-02-19 中国石油天然气股份有限公司 Method and device for removing wave equation simulation direct wave from full-waveform inversion

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Multiscale seismic waveform inversion;Carey Bunks et al.;《GEOPHYSICS》;19951031;第60卷(第5期);第1457－1473页 *
基于声学全波形反演的油气藏地震成像方法;石玉梅，等;《地球物理学报》;20140228;第57卷(第2期);第607－617页 *
频率域反射波全波形速度反演;成景旺，等;《地球科学－中国地质大学学报》;20130330;第38卷(第2期);第391－396页 *
高分辨率的频率空间域声波全波形速度反演－理论模型;廖建平，等;《地球物理学进展》;20111031;第26卷(第5期);第1691页 *

Also Published As

Publication number Publication date
CN104977605A (en) 2015-10-14

Similar Documents

Publication Publication Date Title
CN104407378B (en) Method and device for inversing anisotropy parameters
CN105445791B (en) A kind of formation pore pressure Forecasting Methodology based on a variety of seismic properties
CN101251604B (en) Method for analyzing and NMO correcting two parameters transformation wave speed
CN102466816B (en) Inversion method for stratum elasticity constant parameter of pre-stack seismic data
CN102393532B (en) Seismic signal inversion method
CN101630018B (en) Method for processing seismic exploration data in process of controlling full acoustic wave equation inversion
CN107894612B (en) A kind of the sound impedance inversion method and system of Q attenuation by absorption compensation
CN104570082B (en) Extraction method for full waveform inversion gradient operator based on green function characterization
Dadashpour et al. A derivative-free approach for the estimation of porosity and permeability using time-lapse seismic and production data
CN101201409B (en) Method for revising earthquake data phase
CN107329171A (en) Depth Domain reservoir seismic inversion method and device
CN104237937B (en) Pre-stack seismic inversion method and system thereof
CN103713315A (en) Seismic anisotropy parameter full waveform inversion method and device
CN103499835A (en) Method for inverting near-surface velocity model by utilizing preliminary waveforms
CN103257361A (en) Petroleum-gas prediction method and system based on Zoeppritz equation approximate expression
CN101630013A (en) Method for inverting Poisson ratio parameters of pre-stack seismic data
CN110501744A (en) Hydrocarbon source rock organic carbon geophysics quantitative forecasting technique, device, equipment and storage medium
CN106353797A (en) High-precision earthquake forward modeling method
CN105319589A (en) Full-automatic three-dimensional tomography inversion method using a local event slope
CN104155693A (en) Angle gather seismic response numerical computation method of reservoir fluid fluidity
CN106033124A (en) Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization
CN106483559A (en) A kind of construction method of subsurface velocity model
CN111007567A (en) Sand shale thin interbed prediction method and system based on seismic waveform inversion
CN107390270A (en) A kind of AVA analysis methods based on elastic wave reverse-time migration ADCIGs
CN104597489A (en) Seismic source wavelet optimal setting 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
GR01 Patent grant
GR01 Patent grant