CN104950332B - A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient - Google Patents
A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient Download PDFInfo
- Publication number
- CN104950332B CN104950332B CN201510342764.3A CN201510342764A CN104950332B CN 104950332 B CN104950332 B CN 104950332B CN 201510342764 A CN201510342764 A CN 201510342764A CN 104950332 B CN104950332 B CN 104950332B
- Authority
- CN
- China
- Prior art keywords
- layer
- wave
- stress
- formula
- displacement
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Reverberation, Karaoke And Other Acoustics (AREA)
Abstract
The invention discloses a kind of calculating side of Elastic Multi-layer medium midplane wave reflection coefficient with laxative remedy, it is incident from medium n+1 to purpose series of strata that step following (1) sets plane harmonization, reflected P-wave and shear wave are then produced in n+1 media, transmitted P-wave and shear wave, each layer medium are produced in medium 1 can be write as scalar potential and vector potential form；(2) displacement, stress and relation existing for bit shift are determined；(3) scalar potential in step (1) and vector potential are brought into and obtained in step (2) in (n+1)th layer and the 1st layer of displacement and stress relation formula；(4) according in (n+1)th layer obtained in (3) and the 1st layer of displacement and stress relation formula obtain the displacement comprising each layer coefficient of elasticity, stress transfer matrix, it is solved, can be reflected and transmission coefficient.The present invention has taken into full account that the more subwaves of elastic layer system and conversion involve the influence of thickness, frequency to reflectance factor, is suitable for thin layer AVO analysis modes, and computational efficiency is high.
Description
Technical field
The present invention relates to a kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient, belong to thin layer AVO simulations
Technical field.
Background technology
It is two major classes that simulation thin layer AVO, which is divided to, (1) time-domain convolution model, in the case of a certain special angle, to time-domain
Sampled point calculates reflectance factor one by one, carries out convolution using wavelet and reflectance factor, forms the angle and correspond to seismic channel, this method
Although calculating fast, formula is simple, easily realizes, it can be difficult to embodying thin-bed effect, simulation precision is relatively low.(2) using elasticity
The Integral Solution finite difference calculus simulation thin layer wave field of wave equation, can simulate thin-bed effect, but computational efficiency is relatively low, is built
Mould has a great influence, and industrial production is difficult to practical.
Conventional Zoeppritz equations calculate reflectance factor and are all based on single interface, no gauge variation, although can also count
AVO features are calculated, but are unable to impact effect of the Direct Analysis thickness to each frequency component.
The content of the invention
In view of the deficienciess of the prior art, it is an object of the present invention to provide a kind of Elastic Multi-layer medium midplane wave reflection system
Several computational methods, take into full account that the more subwaves of elastic layer system and conversion involve the influence of thickness, frequency to reflectance factor, suitably
In thin layer AVO analysis modes, computational efficiency is high, and simulation precision is high.
To achieve these goals, the present invention is to realize by the following technical solutions：
A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient of the present invention, specifically include following step
Suddenly：
(1) plane harmonization P is seteFrom medium n+1 to purpose series of strata withIncidence, then it is vertical that reflection is produced in n+1 media
Ripple and reflection wave, its reflectance factor are respectively PrAnd Sr, transmitted P-wave and transmitted shear wave, its transmission coefficient are produced in medium 1
Respectively PtAnd St, each layer medium can be write as scalar potential and vector potential form；
(2) under two-dimensional case, the relational expression between displacement, stress and bit shift is determined；
(3) scalar potential is related in step (1) and vector potential is brought into step (2), can obtain in (n+1)th layer and
1st layer of displacement and stress relation formula；
(4) according in (n+1)th layer obtained in (3) and the 1st layer of displacement and stress relation formula, can obtain including each layer
The displacement of coefficient of elasticity, stress transfer matrix, are then solved to it, you can obtain reflectance factor, transmission coefficient.
Step (1) specifically includes following steps：
If interlayer includes n-1 level course, sequence numbering is 2,3 ... n, bottom medium are 1 layer；In each layer compressional wave and
Transverse wave speed respectively use subscripting α, β represent, under be designated as sequence number, take x coordinate axle to be coincided with n-th layer bottom interface；If one
Individual ripple incides interlayer top surface from n+1 layers direction, and note compressional wave incidence wave is Pe, now, there is a longitudinal wave reflection ripple P in n+1 layersr,
Transverse wave reflection ripple Sr, there is compressional wave transmitted wave P in 1 layertWith shear wave transmitted wave StIf compressional wave and shear wave are to each layer incidence angle
Then bit shift total in n+1 layer media is：
Wherein：
σ=ω/ckn+1=ω/αn+1, Kn+1=ω/βn+1
Wherein, ω is frequency, φeFor incident longitudinal wave bit function, φrReflected P-wave bit function,Enter for (n+1) layer
Penetrate the amplitude of compressional wave, e is constant, j is imaginary unit, z be direction,It is direction for the amplitude of reflected P-wave, x, t is when being
Between, ψrRepresent reflection wave bit function,For the amplitude of reflection wave；
C is apparent velocity of the ripple along interface direction, is all equal to each ripple on interface according to Snell's law
's；
Compressional wave and shear wave bit function expression formula are respectively in n-layer medium：
Wherein,For incident longitudinal wave amplitude,For reflected P-wave amplitude, ψeFor the bit function of incident shear wave,To enter
Penetrate shear wave amplitude,For reflection wave amplitude, knThe horizontal wave number of n-th layer compressional wave, KnThe horizontal wave number of n-th layer shear wave；
Then transmitted wave bit function expression formula is respectively in 1 layer of medium：
Wherein,
φtFor transmitted P-wave bit function,For transmitted P-wave amplitude, ψtFor transmitted shear wave bit function,For transmitted shear wave
Amplitude, k1For the horizontal wave number of first layer compressional wave, K1For the horizontal wave number of first layer shear wave.
In step (2), in two-dimensional case, the relational expression of displacement, stress and bit shift is：
Wherein, u represents displacement, σzzRepresent principal strain along z-axis, τzxShear stress, z are represented as coordinate direction, λ and μ Lames
Constant.
Step (3) specifically includes following steps：
If n-layer thickness is h, displacement component and the components of stress in n-layer are calculated, takes its z=h value, is on n-th layer top surface
Displacement and components of stress value, be designated as u(n)、ω(n)、Bring formula (2-1-2), (2-1-3), (2-1-4) into (2-
1-6) formula, matrix can be obtained：Wherein p=dnH, Q=snH, dnIt is thickness, s for n-th layer compressional wave vertical wavenumber, hnFor n-th layer
Shear wave vertical wavenumber；
Wherein
Wherein, s is the 1st layer of shear wave vertical wavenumber and d is the 1st layer of compressional wave vertical wavenumber, λnAnd μnIt is normal for the Lame of n-th layer
Number
The value of fetch bit shifting and the components of stress at z=0, displacement component and the components of stress on n-layer bottom surface can be obtained, according to point
Interfacial displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer top surface, is designated as u(n-1)、ω(n-1)、 As z=0, p=0, Q=0, can be obtained by matrix (2-1-8)
(n-1) displacement on layer top surface is represented by with components of stress value：
The pass between the displacement component and the components of stress of (n) layer and (n-1) layer top surface can be set up by above-mentioned formula
System, therefore, solving equation group (2-1-10) can have：
Wherein (bij) representInverse matrix：
Wherein, K is the horizontal wave number of shear wave
Formula (2-1-11) substitution (2-1-7) can be obtained
Take mark (aij) representing matrix product (Bij)(bij), 4 × 4 rank square formations are similarly, each element is：
a11=2sin2γcosp-cos2γcosQ
a12=j (tan θ cos2 γ sinp+sin2 γ sinQ)
a22=cos2 γ cosp+2sin2γcosQ
a31=2j ω ρ β sin γ (cosp-cosQ) cos2 γ
a33=cos2 γ cosp+2sin2γcosQ
a34=2j ρ β2(cos2γtanθsinp-sin2γsinQ)
a44=2sin2γcosp+cos2γcosQ (2-1-13)
Wherein, γ=is, θ=idIncidence angle of the ripple in layer is represented respectively, so have sin γ=σ/K, sin θ=σ/k,
α, β are p-and s-wave velocity, and ρ is Media density, established in formula (2-1-12) (n) layer and (n-1) layer top surface displacement component and
Relation between the components of stress, calculate (aij) matrix element, the parameter values of (n) layer will be used in formula (2-1-13), therefore,
Coefficient matrix represents to have with corner brace (n) in formula (2-1-12)：
The recurrence formula of the displacement component and the components of stress on each interface of interlayer has been obtained, has met relational expression (2-1-14)
The displacement for being equivalent to meet on interlayer interface and the stress condition of continuity.
Step (4) specifically includes following steps：
For the transmission coefficient for determining interlayer top surface reflectance factor and being traveled below through interlayer in interlayer bottom surface, according to public affairs
The displacement component and the relation of the components of stress that formula (2-1-14) is established on (n) layer top surface and (1) layer top surface：
(2-1-15) is deformed
Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), can be obtained after solution
Obtain the reflectance factor and transmission coefficient of interlayer；
In view of AijIt is plural number, makes A11=R11+iI11, A12=R12+iI12, A13=R13+iI13...
Wherein, ρn+1(n+1)th layer of density,(n+1)th layer of the horizontal wave number of shear wave
Wherein B is referred to as dematrix, and its each element is：
B11=B22=-i σ, B12=isn+1
B21=idn+1,
Wherein, Dn1For plural number；
An1Imaginary part, Dn2For plural An2Imaginary part, d1For thickness, Dn3For plural An3Imaginary part, ρ1For the close of first layer
Degree,Represent horizontal wave number, the C of first layer shear waven4For plural An4Imaginary part, s1First layer shear wave vertical wavenumber, Cn1It is multiple
Number An1Imaginary part, Cn2For plural An2Imaginary part
Solve dematrix expression formula formula (2-1-17) and can obtain longitudinal wave reflection system of the plane harmonization in elastic layer system
Number V, transverse wave reflection coefficient V/, shear wave transmission coefficient W/, compressional wave transmission coefficient W.
Reflectance factor spectral theory of the invention based on elastic layer system, displacement and stress recurrence formula using elastic layer system,
For solid multi-layer medium, P-wave And S reflection, transmission coefficient formula are deduced, this method can calculate different incidence angles
The Reflection Coefficient of Planar Wave of degree, on petroleum exploration, geologic thin (mutual) layer stratum knot can be represented with elastic layer system model
Structure, therefore the method for the present invention can be used for calculating thin (mutual) layer AVO (Amplitude versus offset), can solve
The influence of more subwaves, converted wave to thin layer AVO, while can be with caused by quantitative simulation different reservoir thickness, different frequency change
AVO changes, while can be used for explanation of seismic road set information and be based on well zoeppritz equation forward modelings road with routine business software
Collect unmatched problem；The present invention can be as an important supplement of current oil exploration industry commercial software.
Brief description of the drawings
Fig. 1 is reflection and the transmission model of layered medium；
Fig. 2 is the reflection and transmission of single thin layer；
Fig. 3 is the reservoir AVO characteristic simulations under different reservoir depth information；
Fig. 4 is the reservoir AVO characteristic simulations in the case of different incident wave frequency rates；
Fig. 5 is that single-layer medium is consistent with traditional commerce software.
Embodiment
To be easy to understand the technical means, the inventive features, the objects and the advantages of the present invention, with reference to
Embodiment, the present invention is expanded on further.
Conventional Zoeppritz equations calculate reflectance factor and are all based on single interface, no gauge variation, although can also count
AVO features are calculated, but are unable to impact effect of the Direct Analysis thickness to each frequency component.The plane wave reflection system of Elastic Multi-layer medium
Number has taken into full account that the more subwaves of elastic layer system and conversion involve the influence of thickness, frequency to reflectance factor, therefore is more suitable for grinding
Study carefully thin layer, the AVO analyses of thin interbed and forward simulation.Therefore the plane wave that the present invention is described in detail in Elastic Multi-layer medium is anti-
Coefficient and forward modelling method are penetrated, and single-layer model has been carried out to derive checking.
Reflection Coefficient of Planar Wave computational methods in Elastic Multi-layer medium are as follows：
An interlayer among two semo-infinite elastic fluids.The interlayer is made up of multiple horizontal levels.Provided with one
Plane wave then produces an interlayer back wave and propagated in upper dielectric, while ripple also will from upper dielectric incidence interlayer top surface
Through interlayer, a transmitted wave in bottom Propagation is produced.Calculate the reflectance factor and transmission coefficient of interlayer.
As shown in figure 1, interlayer is made up of n-1 level course, sequence numbering is 2,3 ... n；Bottom medium is 1 layer.Respectively
Compressional wave and transverse wave speed use the α of subscripting, β to represent respectively in layer, under be designated as sequence number.Take x coordinate axle and n-th layer bottom interface phase
Overlap.There is a P-SV wave systems system to incide interlayer top surface from n+1 layers direction, note compressional wave incidence wave is Pe, now, in n+1 layers
In have a longitudinal wave reflection ripple Pr, transverse wave reflection ripple Sr；There is compressional wave transmitted wave P in 1 layertWith shear wave transmitted wave St.If compressional wave and horizontal stroke
Ripple is to each layer incidence angle
Then bit shift total in n+1 layer media is：
Wherein：σ=ω/
ckn+1=ω/αn+1, Kn+1=ω/βn+1
C is apparent velocity of the ripple along interface direction, is all equal to each ripple on interface according to Snell's law
's.Compressional wave and shear wave bit function expression formula are respectively in n-layer medium：
Then transmitted wave bit function expression formula is respectively in 1 layer of medium：
Wherein
In two-dimensional case, the relational expression of displacement, stress and bit shift is：
If n-layer thickness is h, displacement component and the components of stress in n-layer are calculated, takes its z=h value, is on n-th layer top surface
Displacement and components of stress value, be designated as u(n)、ω(n)、Bring formula (2-1-2), (2-1-3), (2-1-4) into (2-
1-6) formula, matrix can be obtained：Wherein p=dnH, Q=snH,
Wherein
On the one hand, fetch bit shifting and the components of stress can obtain the displacement component and stress on n-layer bottom surface in the value at z=0 for order
Component.According to interface displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer top surface, is designated as u(n-1)、
ω(n-1)、As z=0, p=0, Q=0, can be obtained by matrix (2-1-8)
(n-1) displacement on layer top surface is represented by with components of stress value：
The pass between the displacement component and the components of stress of (n) layer and (n-1) layer top surface can be set up by above-mentioned formula
System.Therefore, solving equation group (2-1-10) can have：
Wherein (bij) representInverse matrix：
Formula (2-1-11) substitution (2-1-7) can be obtained
Take mark (aij) representing matrix product (Bij)(bij), 4 × 4 rank square formations are similarly, each element is：a11=2sin2γ
cosp-cos2γcosQ
a12=j (tan θ cos2 γ sinp+sin2 γ sinQ)
a22=cos2 γ cosp+2sin2γcosQ
a31=2j ω ρ β sin γ (cosp-cosQ) cos2 γ
a33=cos2 γ cosp+2sin2γcosQ
a34=2j ρ β2(cos2γtanθsinp-sin2γsinQ)
a44=2sin2γcosp+cos2γcosQ (2-1-13)
Wherein γ=is, θ=idIncidence angle of the ripple in layer is represented respectively, so have sin γ=σ/K, sin θ=σ/k.Its
Its parameter, such as α, β are p-and s-wave velocity, and ρ is Media density, and ω is frequency.(n) layer and (n- are established in formula (2-1-12)
1) relation between layer top surface displacement component and the components of stress, (a is calculatedij) matrix element, (n) will be used in formula (2-1-13)
The parameter values of layer.Therefore, coefficient matrix represents to have with corner brace (n) in formula (2-1-12)：
The recurrence formula of the displacement component and the components of stress on each interface of interlayer is obtained.Meet relational expression (2-1-14)
The displacement for being equivalent to meet on interlayer interface and the stress condition of continuity.
For the transmission coefficient for determining interlayer top surface reflectance factor and being traveled below through interlayer in interlayer bottom surface, according to public affairs
The displacement component and the relation of the components of stress that formula (2-1-14) is established on (n) layer top surface and (1) layer top surface：
(2-1-15) is deformed
Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), can be obtained after solution
Obtain the reflectance factor and transmission coefficient of interlayer.
In view of AijIt is plural number, makes A11=R11+iI11, A12=R12+iI12, A13=R13+iI13...
Wherein B is referred to as dematrix, and its each element is：
B11=B22=-i σ, B12=isn+1
B21=idn+1,
Solve dematrix expression formula formula (2-1-17) and can obtain longitudinal wave reflection system of the plane harmonization in elastic layer system
Number V, transverse wave reflection coefficient V/, shear wave transmission coefficient W/, compressional wave transmission coefficient W.
Fig. 3 is achievement of the longitudinal wave reflection coefficient corresponding to single thin film model with angle change, and in figure, abscissa is angle
Degree, 60 degree are incremented to from 0 degree, and ordinate is longitudinal wave reflection coefficient, is changed between 5 meters to 30 meters of thickness of thin layer, can be with from figure
Find out that thin layer AVO forms are similar corresponding to different reservoir thickness, but the intercept of four curves is different, the gradients of four curves is not yet
Together, this requires that we judge to consider the thickness change of thin layer during the AVO of thin layer.
Fig. 4 is achievement of the longitudinal wave reflection coefficient corresponding to single thin film model with angle change, and in figure, abscissa is angle
Degree, 45 degree is incremented to from 0 degree, ordinate is longitudinal wave reflection coefficient, and frequency, can from figure by changing between 10hz to 90hz
It is similar to go out thin layer AVO forms corresponding to different frequency, 10 curves are all monotonic increases, but each not phase of gradient of 10 curves
Together, this requires that we judge to consider the dominant frequency change of incidence wave during the AVO of thin layer, and especially same set of thin reservoir changes in buried depth
In the case of larger.
To examine the correctness of above-mentioned theory, by the reflection system in the case of the single thin layer of theory deduction vertical incidence above
Number.
If have among two half infinite mediums a thickness be h solid interlayer, calculate the interlayer reflectance factor and thoroughly
Penetrate coefficient.As shown in Fig. 2 being sandwiched between a thin layer in 1 and 3 two half infinite medium, its thickness is h.There is a compressional wave from Jie
The direction of matter 3 impinges perpendicularly on thin layer top surface, will produce a back wave p propagated in medium 3rWith propagate in medium 1
Transmitted wave pt, according to the coordinate system shown in figure, each ripple expression formula is：
There is bit function to the 3rd medium：
φ3=φe+φr (2-2-2)
There is bit function to the 1st medium：
φ1=φt (2-2-3)
The boundary condition that they should meet, according to (2-1-15), it is
In the case of vertical incidence, m=0, τzx=0.By formula (2-2-2), formula (2-2-3) brings formula (2- into
2-4) formula, have on z=0：
Can have on z=h：
Boundary condition equation group is：
Solution of equations is exactly required reflectance factor and transmission coefficient, and they are：
Wherein TPPTransmission coefficient takes the displacement amplitude of transmitted wave and incidence wave ratio.By the parameter of thin-layered medium 2 and incidence angle number
According to θ=0, γ=0, bringing formula (2-1-13) into, can have：
Formula (2-2-8) is substituted into formula (2-2-6), formula (2-2-7), can be obtained
As h=0, reflectance factor and transmission coefficient during compressional wave vertical incidence during an interface can be obtained：
Its result is it is clear that correctly.
Fig. 5 is situation of the reflectance factor under the single interface conditions calculated using zoppritz equations with angle change, is schemed
In, abscissa is angle, and for scope from 0 degree to 60 degree, ordinate is reflectance factor, change of entirely successively decreasing.In order to verify the present invention
Correctness, we use the computational methods of the present invention, thin intermediate thickness are taken as 0 meter, then model regression of the invention is single boundary
Surface model, reflectance factor calculate with angle change situation using the program of the present invention, by comparing, method of the invention with
Zoeppritz equations are identicals for both single INTERFACE MODELs, and this also demonstrates the correctness of the present invention.
The general principle and principal character and advantages of the present invention of the present invention has been shown and described above.The technology of the industry
Personnel are it should be appreciated that the present invention is not limited to the above embodiments, and the simply explanation described in above-described embodiment and specification is originally
The principle of invention, without departing from the spirit and scope of the present invention, various changes and modifications of the present invention are possible, these changes
Change and improvement all fall within the protetion scope of the claimed invention.The claimed scope of the invention by appended claims and its
Equivalent thereof.
Claims (4)
1. a kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient, it is characterised in that specifically include following
Step：
(1) plane harmonization P is seteFrom medium n+1 to purpose series of strata withIncidence, then in n+1 media produce reflected P-wave and
Reflection wave, its reflectance factor are respectively PrAnd Sr, transmitted P-wave and transmitted shear wave, its transmission coefficient difference are produced in medium 1
For PtAnd St, each layer medium can be write as scalar potential and vector potential form；
(2) under two-dimensional case, the relational expression between displacement, stress and bit shift is determined；
(3) scalar potential is related in step (1) and vector potential is brought into step (2), can obtain in (n+1)th layer and the 1st
The displacement of layer and stress relation formula；
(4) according in (n+1)th layer obtained in (3) and the 1st layer of displacement and stress relation formula, can obtain comprising each layer elasticity
The displacement of coefficient, stress transfer matrix, are then solved to it, you can obtain reflectance factor, transmission coefficient；Step (1) has
Body comprises the following steps：
If interlayer includes n-1 level course, sequence numbering is 2,3 ... n, bottom medium are 1 layer；Compressional wave and shear wave in each layer
Velocity of wave respectively use subscripting α, β represent, under be designated as sequence number, take x coordinate axle to be coincided with n-th layer bottom interface；An if ripple
Interlayer top surface is incided from n+1 layers direction, note compressional wave incidence wave is Pe, now, there is a longitudinal wave reflection ripple P in n+1 layersr, shear wave
Back wave Sr, there is compressional wave transmitted wave P in 1 layertWith shear wave transmitted wave StIf compressional wave and shear wave are to each layer incidence angle
Then bit shift total in n+1 layer media is：
Wherein：
σ=ω/c,kn+1=ω/αn+1, Kn+1=ω/βn+1
Wherein, ω is frequency, φn+1For (n+1)th layer of total compressional wave bit function, φeFor incident longitudinal wave bit function, φrReflected P-wave
Bit function,Amplitude, e for (n+1) layer incident longitudinal wave are constant, j is imaginary unit, z be direction,For reflection
Amplitude, the x of compressional wave are direction, t is time, ψn+1For (n+1)th layer of total shear wave bit function, ψrRepresent reflection wave bit function,For the amplitude of reflection wave,For (n+1)th layer of shear wave incidence angle；dn+1For (n+1)th layer of compressional wave vertical wavenumber, sn+1For n-th
+ 1 layer of shear wave vertical wavenumber, kn+1For the horizontal wave number of (n+1)th layer of compressional wave, Kn+1For the horizontal wave number of (n+1)th layer of shear wave；
C is apparent velocity of the ripple along interface direction, is all equal to each ripple on interface according to Snell's law；
Compressional wave and shear wave bit function expression formula are respectively in n-layer medium：
Wherein, φnFor the total compressional wave bit function of n-th layer, ψnFor the total shear wave bit function of n-th layer,For incident longitudinal wave amplitude,For reflected P-wave amplitude, ψeFor the bit function of incident shear wave,For incident shear wave amplitude,For reflection wave amplitude, kn
The horizontal wave number of n-th layer compressional wave, KnThe horizontal wave number of n-th layer shear wave；dnFor n-th layer compressional wave vertical wavenumber, snFor n-th layer shear wave
Vertical wavenumber；
Then transmitted wave bit function expression formula is respectively in 1 layer of medium：
Wherein,
φ1For the 1st layer of total compressional wave bit function, ψ1For the 1st layer of total shear wave bit function,
φtFor transmitted P-wave bit function,For transmitted P-wave amplitude, ψtFor transmitted shear wave bit function,Shaken for transmitted shear wave
Width, k1For the horizontal wave number of first layer compressional wave, K1For the horizontal wave number of first layer shear wave, d1For the 1st layer of compressional wave vertical wavenumber, s1For the 1st
Layer shear wave vertical wavenumber.
2. the computational methods of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 1, it is characterised in that step
Suddenly in (2), in two-dimensional case, the relational expression of displacement, stress and bit shift is：
Wherein, u represents displacement, σzzRepresent principal strain along z-axis, τzxShear stress, z are represented as coordinate direction, λ and μ Lame constants.
3. the computational methods of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 2, it is characterised in that step
Suddenly (3) specifically include following steps：
If n-layer thickness is h, displacement component and the components of stress in n-layer are calculated, takes its z=h value, is the position on n-th layer top surface
Shifting and components of stress value, are designated as u(n)、ω(n)、Formula (2-1-2), (2-1-3), (2-1-4) are brought into (2-1-6)
Formula, matrix can be obtained：Wherein p=dnH, Q=snH, dnIt is thickness, s for n-th layer compressional wave vertical wavenumber, hnFor n-th layer shear wave
Vertical wavenumber；
Wherein
Wherein, s is the 1st layer of shear wave vertical wavenumber, and d is the 1st layer of compressional wave vertical wavenumber, λnAnd μnFor the Lame constants of n-th layer；μn-1
For (n-1)th layer of Lame constants；
The value of fetch bit shifting and the components of stress at z=0, displacement component and the components of stress on n-layer bottom surface can be obtained, according to interface
Displacement and the stress condition of continuity, it should be equal with the analog value of (n-1) layer top surface, is designated as u(n-1)、ω(n-1)、
As z=0, p=0, Q=0, can be obtained by matrix (2-1-8)
(n-1) displacement on layer top surface is represented by with components of stress value：
The relation between the displacement component and the components of stress of (n) layer and (n-1) layer top surface can be set up by above-mentioned formula, is
This, solving equation group (2-1-10) can have：
Wherein (bij) representInverse matrix, wherein i, j is subscript, and i value is respectively 1,2,3,4；J value difference
For 1,2,3,4：
Wherein, K is the horizontal wave number of shear wave
Formula (2-1-11) substitution (2-1-7) can be obtained
Take mark (aij) representing matrix product (Bij)(bij), 4 × 4 rank square formations are similarly, wherein i, j is subscript, and i value is divided
Wei 1,2,3,4；J value is respectively 1,2,3,4：Each element is：
a11=2sin2γcos p-cos 2γcos Q
a12=j (tan θ cos2 γ sin p+sin2 γ sin Q)
a22=cos2 γ cos p+2sin2γcos Q
a31=2j ω ρ β sin γ (cos p-cos Q) cos2 γ
a33=cos2 γ cos p+2sin2γcos Q
a34=2j ρ β2(cos2γtanθsin p-sin2γsin Q)
a44=2sin2γcos p+cos2γcos Q (2-1-13)
Wherein, γ=is, θ=idIncidence angle of the ripple in layer is represented respectively, so there is sin γ=σ/K, sin θ=σ/k, α, β are
P-and s-wave velocity, ρ are Media density, and (n) layer and (n-1) layer top surface displacement component and stress are established in formula (2-1-12)
Relation between component, calculate (aij) matrix element, the parameter values of (n) layer will be used in formula (2-1-13), therefore, in formula
Coefficient matrix represents to have with corner brace (n) in (2-1-12)：
The recurrence formula of the displacement component and the components of stress on each interface of interlayer has been obtained, has met relational expression (2-1-14) equivalence
In meeting the displacement on interlayer interface and the stress condition of continuity.
4. the computational methods of Elastic Multi-layer medium midplane wave reflection coefficient according to claim 3, it is characterised in that step
Suddenly (4) specifically include following steps：
For the transmission coefficient for determining interlayer top surface reflectance factor and being traveled below through interlayer in interlayer bottom surface, according to formula
The displacement component and the relation of the components of stress that (2-1-14) is established on (n) layer top surface and (1) layer top surface：
(2-1-15) is deformed
Formula (2-1-1), formula (2-1-5) are substituted into formula (2-1-16) and obtain matrix equation (2-1-17), can be pressed from both sides after solution
The reflectance factor and transmission coefficient of layer；
In view of AijIt is plural number, makes A11=R11+iI11, A12=R12+iI12, A13=R13+iI13...
Wherein, ρn+1(n+1)th layer of density,(n+1)th layer of the horizontal wave number of shear wave
Wherein B is referred to as dematrix, and its each element is：
B11=B22=-i σ, B12=isn+1
B21=idn+1,
Wherein, Dn1For plural An1Imaginary part, Dn2For plural An2Imaginary part, d1For first layer compressional wave vertical wavenumber, Cn3For plural An3
Imaginary part, ρ1For the density of first layer, K1 2Represent horizontal wave number, the C of first layer shear waven4For plural An4Imaginary part, s1First layer
Shear wave vertical wavenumber, Cn1For plural An1Imaginary part, Cn2For plural An2Imaginary part, longitudinal wave reflection coefficient V, transverse wave reflection coefficient V/、
Shear wave transmission coefficient W/, compressional wave transmission coefficient W；
Solve dematrix expression formula formula (2-1-17) can obtain plane harmonization elastic layer system longitudinal wave reflection coefficient V,
Transverse wave reflection coefficient V/, shear wave transmission coefficient W/, compressional wave transmission coefficient W.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510342764.3A CN104950332B (en) | 2015-06-18 | 2015-06-18 | A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510342764.3A CN104950332B (en) | 2015-06-18 | 2015-06-18 | A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104950332A CN104950332A (en) | 2015-09-30 |
| CN104950332B true CN104950332B (en) | 2018-02-02 |
Family
ID=54165128
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510342764.3A Expired - Fee Related CN104950332B (en) | 2015-06-18 | 2015-06-18 | A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104950332B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105629301B (en) * | 2016-03-29 | 2018-02-09 | 中国地质大学（北京） | Thin layer elastic wave reflex coefficient fast solution method |
| CN106054247B (en) * | 2016-05-25 | 2020-09-29 | 中国石油集团东方地球物理勘探有限责任公司 | Method for calculating high-precision reflection coefficient based on converted wave seismic data |
| CN106772578B (en) * | 2016-12-07 | 2018-11-09 | 中国矿业大学(北京) | A kind of method and apparatus of synthetic seismogram |
| CN109324343A (en) * | 2017-08-01 | 2019-02-12 | 中国石油化工股份有限公司 | A kind of analogy method and system of thin layer displacement multi-wave seismic wave field |
| CN108196301B (en) * | 2017-12-04 | 2020-07-10 | 中国石油天然气股份有限公司 | Method and device for acquiring gather with amplitude changing along with offset |
| CN111830566B (en) * | 2020-06-12 | 2022-02-11 | 中国海洋大学 | Parameter matching virtual reflection suppression method and marine seismic exploration system |
| CN112130207B (en) * | 2020-09-25 | 2021-07-20 | 中国科学院武汉岩土力学研究所 | Method for calculating underground vibration from ground vibration based on spherical charging condition |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103149586A (en) * | 2013-02-04 | 2013-06-12 | 西安交通大学 | Tilted layered viscoelasticity dielectric medium wave field forward modelling method |
| CN104570072A (en) * | 2013-10-16 | 2015-04-29 | 中国石油化工股份有限公司 | Method for modeling reflection coefficient of spherical PP wave in viscoelastic medium |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140324354A1 (en) * | 2013-04-29 | 2014-10-30 | King Fahd University Of Petroleum And Minerals | Transmission coefficient method for avo seismic analysis |
-
2015
- 2015-06-18 CN CN201510342764.3A patent/CN104950332B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103149586A (en) * | 2013-02-04 | 2013-06-12 | 西安交通大学 | Tilted layered viscoelasticity dielectric medium wave field forward modelling method |
| CN104570072A (en) * | 2013-10-16 | 2015-04-29 | 中国石油化工股份有限公司 | Method for modeling reflection coefficient of spherical PP wave in viscoelastic medium |
Non-Patent Citations (4)
| Title |
|---|
| 利用井中偶极声源远场辐射特性的远探测测井;唐晓明 等;《地球物理学报》;20120831;第55卷(第8期);第2798-2807页 * |
| 多分量地震勘探技术新进展;刘秀娟 等;《西部探矿工程》;20071231(第1期);第123-124,127页 * |
| 弹性层系的反射系数正演;梁立锋 等;《物探与化探》;20070430;第31卷(第2期);第174-177页 * |
| 波在多层弹性介质中一些特性的研究;王育生 等;《地震工程与工程振动》;19831231;第3卷(第4期);第69-82页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104950332A (en) | 2015-09-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104950332B (en) | A kind of computational methods of Elastic Multi-layer medium midplane wave reflection coefficient | |
| CN104407378B (en) | Method and device for inversing anisotropy parameters | |
| CN104267429B (en) | The method and device of stressor layer definitely | |
| Tariq et al. | A new artificial intelligence based empirical correlation to predict sonic travel time | |
| CN106368691B (en) | Three-dimensional abnormal pore pressure prediction method based on rock physics seismic information | |
| Sallarès et al. | Relative contribution of temperature and salinity to ocean acoustic reflectivity | |
| CN102426390A (en) | Method for determining reserve volume of nonhomogeneous sandstone reservoir | |
| CN104237937B (en) | Pre-stack seismic inversion method and system thereof | |
| CN102156297B (en) | Fluid substitution method based on sandstone reservoir post-stack seismic data | |
| CN104502997A (en) | Method for using fracture density curve to forecast fracture density body | |
| CN102253415A (en) | Method for establishing earthquake response mode based on fracture equivalent medium model | |
| CN106353797A (en) | High-precision earthquake forward modeling method | |
| CN105911584B (en) | Implicit staggered-grid finite difference elastic wave numerical simulation method and device | |
| CN103513277A (en) | Earthquake stratum fracture crack density retrieval method and system | |
| Sun et al. | 2-D poroelastic wave modelling with a topographic free surface by the curvilinear grid finite-difference method | |
| Omidyeganeh et al. | Large eddy simulation of interacting barchan dunes in a steady, unidirectional flow | |
| CN106291689A (en) | A kind of extract the processing method of geological data frequency dispersion attribute, device and prognoses system | |
| CN106368687A (en) | Shale reservoir brittleness evaluating method | |
| Ge et al. | An efficient approach for simulating wave propagation with the boundary element method in multilayered media with irregular interfaces | |
| Golubev et al. | Simulation of seismic responses from fractured MARMOUSI2 model | |
| Ouadfeul et al. | Shale gas reservoirs characterization using neural network | |
| CN107576985A (en) | A kind of method and apparatus of seismic inversion | |
| CN102768367B (en) | Two-phase medium amplitude versus offset (AVO) forward modeling method based on triple constraints | |
| CN102288993A (en) | Fluid replacing method of pre-stack earthquake data based on sandstone oil reservoir | |
| Ye et al. | Simulation of field injection experiments in heterogeneous unsaturated media using cokriging and artificial neural network |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180202 Termination date: 20200618 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |