CN108333628A - Elastic wave least square reverse-time migration method based on regularization constraint - Google Patents
Elastic wave least square reverse-time migration method based on regularization constraint Download PDFInfo
- Publication number
- CN108333628A CN108333628A CN201810042691.XA CN201810042691A CN108333628A CN 108333628 A CN108333628 A CN 108333628A CN 201810042691 A CN201810042691 A CN 201810042691A CN 108333628 A CN108333628 A CN 108333628A
- Authority
- CN
- China
- Prior art keywords
- wave
- gradient
- equation
- elastic wave
- time migration
- 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
- 238000000034 method Methods 0.000 title claims abstract description 71
- 238000013508 migration Methods 0.000 title claims abstract description 71
- 230000005012 migration Effects 0.000 title claims abstract description 70
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 238000002939 conjugate gradient method Methods 0.000 claims abstract description 5
- 238000002310 reflectometry Methods 0.000 claims abstract description 5
- 239000012530 fluid Substances 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 238000012986 modification Methods 0.000 claims description 5
- 230000004048 modification Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 238000003672 processing method Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims description 3
- 238000013316 zoning Methods 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 abstract description 22
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000013461 design Methods 0.000 abstract description 2
- 230000001976 improved effect Effects 0.000 description 5
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000005303 weighing Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000003989 dielectric material Substances 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses the elastic wave least square reverse-time migration methods based on regularization constraint.Design new object function；Derive the elastic wave inverse migration operator and reflectance factor gradient formula under fresh target function；Calculate the gradient of reflectance factor；Gradient is handled using conjugate gradient method or quasi-Newton method inversion algorithm；Iteration step length is sought using curve-parabola-fitting method；Reflectivity model is updated, until meeting the condition of convergence.The beneficial effects of the invention are as follows improve imaging resolution and stability by using new full variational regularization constraints policy.
Description
Technical field
The invention belongs to seismic wave migration and imaging techniques fields, are related to more points of geophysics field (especially seismic prospecting)
Imaging precision and resolution ratio are improved in the offset of amount data.
Background technology
The most basic purpose of seism processing is migration imaging, and the quality of image quality directly determines interface location
The signal-to-noise ratio of accuracy, the height of resolution ratio and section.Migration processing can make tilted interface playback, diffracted wave convergence, carry
High lateral resolution.After decades of development, migration technology has gone to prestack from poststack, and depth has been developed to from time-domain
Domain.Classify by Method And Principle, migration technology is divided into：Kirchhoff offsets, F-K offsets, one-way wave wave equation migration and inverse time
It deviates (reverse time migration, RTM).Compared with other methods, reverse-time migration as a kind of prestack, Depth Domain,
Round trip wave wave equation migration method, it is assumed that condition is minimum, precision highest, adapts to arbitrary velocity distribution, aclinal limitation, becomes
The most commonly used imaging method under the conditions of complex geological structure.But conventional reverse-time migration utilizes the adjoint of wave field forward-propagating operator
Operator replaces the inverse of it, is still inaccurate.In addition, by acquisition aperture, underground lighting, data itself (with it is sex-limited, do not advise
Then, noisy) etc. factors influence, in the imaging section of reverse-time migration there are acquisition footprint, the problems such as resolution ratio is low, amplitude unbalance.
To further increase imaging precision, there is least square reverse-time migration.This method passes through minimization simulated reflections
The error of wave and observation back wave seeks best reflectance factor.Due to establishing specific object between reflectance factor and earthquake record
Reason relationship (inverse migration operator), least square reverse-time migration method, which has, preferably protects width ability.The introducing of inverting thought also makes
Obtaining least square reverse-time migration has higher resolution ratio and less migration noise.Least square reverse-time migration is constantly sent out
Exhibition, has been applied in sound wave and viscous sound wave media imaging.But it is suitable for complex dielectrics (elasticity, viscoplasticity, anisotropy
Deng) least square reverse-time migration method it is relatively fewer.Compared with acoustic wave methodogy, the reverse-time migration of elastic wave least square can obtain
Accurate longitudinal wave reflection coefficient and transverse wave reflection coefficient are obtained, can preferably identify lithology and fluid, prediction geological disaster etc., tool
Have broad application prospects.At the same time, longitudinal wave, shear wave, converted wave etc. are a variety of involved in the reverse-time migration of elastic wave least square
A variety of model parameters such as wave mode and Lame Coefficient, speed, density, impedance, the crosstalk phenomenon between wave field and parameter are tight
Weight.Therefore, existing elastic wave least square reverse-time migration method often hardly results in satisfied result.The resolution of migrated section
Rate is low, it is serious to protect prismatic error, crosstalk noise.In addition, as a kind of Multi-parameters conversion, elastic wave least square reverse-time migration it is steady
It is qualitative also to be improved.
Invention content
The purpose of the present invention is to provide the elastic wave least square reverse-time migration methods based on regularization constraint.
The technical solution adopted in the present invention is to follow the steps below：
(1) new object function is designed；
Fresh target function includes two：Simulated reflections wave and the difference and regularization term for observing back wave, the contribution of the two
It is adjusted by regularization coefficient；
(2) the elastic wave inverse migration operator and reflectance factor gradient formula under fresh target function are derived；
It is public to the gradient of reflectance factor that elastic wave adjoint equation/inverse migration operator and object function are derived based on adjoint method
Formula；
(3) gradient of reflectance factor is calculated；
It specifically includes：Source wavefield forward-propagating；Back wave residual error backpropagation；Forward and reverse wave field correlation obtains often
Advise gradient；Along with regularization term is to the gradient of reflectance factor；
(4) conjugate gradient method or quasi-Newton method inversion algorithm is used to handle gradient；
(5) iteration step length is sought using curve-parabola-fitting method；
(6) reflectivity model is updated, until meeting the condition of convergence.
Further, elastic wave least square reverse-time migration process, object function are constrained using TV regularizations in step (1)
For：
Wherein：T is maximum time, and H is zoning, d1、d2And d3For regularization coefficient, β1、β2And β3For stability because
Son, Δ vxWith Δ vzFor simulated reflections wave horizontal component and vertical component,WithFor observation back wave horizontal component and
Vertical component,
Wherein, λ and μ is Lame constants, and ρ is density, and equation (2) indicates the opposite variation of model parameter, dimensionless, can be with
For weighing the size of reflectance factor, elastic wave least square reverse-time migration is exactly to seek optimal Rρ、RλAnd RμProcess, new mesh
Scalar functions include two parts：Simulated reflections wave and the difference and regularization term for observing back wave, the contribution of the two pass through regularization
Coefficient d1、d2And d3It adjusts.
Further, the elastic wave inverse migration operator and reflectance factor gradient formula under fresh target function are derived in step (2)
Method is as follows：
Elastic wave velocity-stress equation is：
Wherein, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is stress vector；
In elastic fluid, for background model parameters [λ, μ, ρ], background wave field [vx,vz,τxx,τzz,τxz] pass through solution side
Journey obtains, and when there are model disturbance [Δ λ, Δ μ, Δ ρ], wave field knots modification is [Δ vx,Δvz,Δτxx,Δτzz,Δτxz],
And meet：
Abbreviation is simultaneously ignored high-order small quantity and is obtained：
To given parameter perturbation [Δ λ, Δ μ, Δ ρ], solves equation and obtain back wave [Δ vx,Δvz,Δτxx,Δτzz,
Δτxz], the inverse migration process as in elastic fluid, in least square reverse-time migration, context parameter [λ, μ, ρ] is constant, the back of the body
Scape wave field is also constant, and the power of back wave is directly determined by parameter perturbation item；
By Rρ、RλAnd RμEquation is substituted into obtain：
Only consider the difference item of simulated reflections wave and observation back wave in object function：
Wherein, reflected wave field [the Δ v of simulationx,Δvz,Δτxx,Δτzz,Δτxz], it is solved using method of Lagrange multipliers
The constrained optimization problem, cost functional become：
Wherein,For Lagrange multiplier function,
Integration by parts obtains：
Wherein,
It enablesCorresponding adjoint equation is obtained, form is as follows：
Object function is about the gradient formula of parameter perturbation：
In the case of TV regularizations, reflectance factor gradient formula becomes：
Further, in step (3)
A. it solves equation (3) and equation (6) obtains reflected wave field [Δ vx,Δvz,Δτxx,Δτzz,Δτxz]T, initial strip
Part is：
[vx(x,z,0),vz(x,z,0),τxx(x,z,0),τzz(x,z,0),τxz(x,z,0)]T=0,
[Δvx(x,z,0),Δvz(x,z,0),Δτxx(x,z,0),Δτzz(x,z,0),Δτxz(x,z,0)]T=0 (15)
B. it solves adjoint equation (12) and obtains backward extension wave fieldFinal value condition is：
C. the gradient by equation (14) calculating target function about reflectance factor.
Further, processing method is as follows in step (4)：
Using L-BFGS methods：
Wherein, HkFor the approximate matrix that Hessian matrix is inverse, H is directly calculatedkLarger calculation amount is needed, passes through several groups here
Column vector carrys out approximate Hk。
Further, to seek iteration step length method in step (5) as follows：
Iteration step length is sought using Parabolic Fit
Wherein, α1And α2To sound out step-length, J1And J2For corresponding target function value, J0For the object function of current iteration
Value calculates J1And J2Need four times additional forward modeling operations；
Then the optimum stepsize of current iteration is：
Further, reflectance factor is updated by following formula in step (6)：
Wherein, mkAnd mk+1The respectively model parameter of current iteration and next iteration：
The beneficial effects of the invention are as follows differentiated by using new full variation (TV) regularization constraint strategy to improve imaging
Rate and stability.The purpose of invention is in order to improve the imaging precision of multi-component seismic data, for subsequent explanation and inverting work
Make to provide reliable migrated section.
Description of the drawings
The flow chart of elastic wave least square reverse-time migrations of the Fig. 1 based on regularization constraint；
Fig. 2 groove models；
The imaging results of Fig. 3 groove model difference offset methods；
Fig. 4 Marmousi models；
The imaging results of Fig. 5 Marmousi model difference offset methods.
Specific implementation mode
The present invention is described in detail With reference to embodiment.
As shown in Figure 1, implement the flow chart of the elastic wave least square reverse-time migration based on regularization constraint for the present invention,
It specifically includes：
(1) new object function is designed.
Fresh target function includes two：Simulated reflections wave and the difference and regularization term for observing back wave, the contribution of the two
It is adjusted by regularization coefficient.
(2) the elastic wave inverse migration operator and reflectance factor gradient formula under fresh target function are derived.
Elastic wave adjoint equation/inverse migration operator and object function pair are derived based on adjoint method (Adjoint method)
The gradient formula of reflectance factor.
(3) gradient of reflectance factor is calculated.
It specifically includes：Source wavefield forward-propagating；Back wave residual error backpropagation；Forward and reverse wave field correlation obtains often
Advise gradient；Along with regularization term is to the gradient of reflectance factor.
(4) gradient is handled using suitable inversion algorithm.
Fore condition processing is carried out to gradient using conjugate gradient method or quasi-Newton method (such as L_BFGS).
(5) iteration step length is sought.
Iteration step length is sought using curve-parabola-fitting method.
(6) reflectivity model is updated, until meeting the condition of convergence.
The new object function method of design is as follows in step (1)：
Regularization Strategy can improve the precision and stability of inverting.Tikhonov regularizations are smooth by applying to model
It constrains (such as derivative), precision can be reduced while improving stability, and full variational regularization method is believed with better high frequency
Cease recovery capability.To obtain high-precision speed and density reflectance factor simultaneously, bullet is constrained using TV regularizations in the present invention
Property wave least square reverse-time migration process.Object function becomes：
Wherein：T is maximum time, and H is zoning, d1、d2And d3For regularization coefficient, β1、β2And β3For stability because
Son.ΔvxWith Δ vzFor simulated reflections wave horizontal component and vertical component,WithFor observation back wave horizontal component and
Vertical component,
Wherein, λ and μ is Lame constants, and ρ is density.Equation (2) indicates the opposite variation of model parameter, dimensionless, can be with
For weighing the size of reflectance factor.Elastic wave least square reverse-time migration is exactly to seek optimal Rρ、RλAnd RμProcess.
Fresh target function (equation 15) includes two parts：The difference (first two) and just of simulated reflections wave and observation back wave
Then change item (latter three), the contribution of the two passes through regularization coefficient (d1、d2And d3) adjust.
Elastic wave inverse migration operator and reflectance factor gradient formula method in step (2) under derivation fresh target function are such as
Under：
Elastic wave velocity-stress equation is：
Wherein, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is stress vector.
In elastic fluid, for background model parameters [λ, μ, ρ], background wave field [vx,vz,τxx,τzz,τxz] can be by asking
3 are solved equation to obtain.When there are model disturbance [Δ λ, Δ μ, Δ ρ], wave field knots modification is [Δ vx,Δvz,Δτxx,Δτzz,Δ
τxz], and meet：
Equation (4a-4e) subtracts each other with equation (3), and abbreviation is simultaneously ignored high-order small quantity and can be obtained：
To given parameter perturbation [Δ λ, Δ μ, Δ ρ], back wave [Δ can be obtained by solving equation (3) and equation (5)
vx,Δvz,Δτxx,Δτzz,Δτxz], the inverse migration process as in elastic fluid.In least square reverse-time migration, background
Parameter [λ, μ, ρ] is constant, and background wave field is also constant, and the power of back wave is directly determined by parameter perturbation item.
By Rρ、RλAnd RμEquation (3) is substituted into obtain：
Only consider the difference item of simulated reflections wave and observation back wave in object function (equation 1)：
Wherein, reflected wave field [the Δ v of simulationx,Δvz,Δτxx,Δτzz,Δτxz] it must satisfy equation (6).Using glug
Bright day multiplier method solves the constrained optimization problem.Cost functional becomes：
Wherein,For Lagrange multiplier function,
Integration by parts equation (8) can obtain：
Wherein,
It enablesCorresponding adjoint equation is obtained, form is as follows：
Object function is about the gradient formula of parameter perturbation：
Regularization constraint item does not influence adjoint equation (remaining as equation 12), but influences gradient formula.TV regularization situations
Under, reflectance factor gradient formula becomes：
The gradient method that reflectance factor is calculated in step (3) is as follows：
A. it solves equation (3) and equation (6) obtains reflected wave field [Δ vx,Δvz,Δτxx,Δτzz,Δτxz]T, initial strip
Part is：
[vx(x,z,0),vz(x,z,0),τxx(x,z,0),τzz(x,z,0),τxz(x,z,0)]T=0,
[Δvx(x,z,0),Δvz(x,z,0),Δτxx(x,z,0),Δτzz(x,z,0),Δτxz(x,z,0)]T=0.
(15)
B. it solves adjoint equation (12) and obtains backward extension wave fieldFinal value condition is：
C. the gradient by equation (14) calculating target function about reflectance factor.
It is as follows to gradient progress processing method using suitable inversion algorithm in step (4)：
The precision and convergence rate of inverting can be improved by being pre-processed to gradient.Common method have conjugate gradient method,
Quasi-Newton method and Newton method etc..Inversion accuracy and computational efficiency in order to balance use L-BFGS methods in of the invention：
Wherein, HkFor the approximate matrix that Hessian matrix is inverse.Directly calculate HkLarger calculation amount is needed, passes through several groups here
Column vector carrys out approximate Hk.L-BFGS inversion algorithm concrete implementations step can refer to the related books and document optimized, this
In repeat no more.
It is as follows that step 5 seeks iteration step length method：
(5) iteration step length is sought using Parabolic Fit in the present invention.
Wherein, α1And α2To sound out step-length, J1And J2For corresponding target function value, J0For the object function of current iteration
Value.Calculate J1And J2Need four times additional forward modeling operations.
Then the optimum stepsize of current iteration is：
It is as follows that step (6) updates reflectivity model method：
Based on above step, reflectance factor is updated by following formula：
Wherein, mkAnd mk+1The respectively model parameter of current iteration and next iteration：
Step (3)-(6) are repeated, are stopped until meeting the condition of convergence (such as residual error is less than 1e-5 or iterations are less than 30)
Only iteration exports final elastic reflectance factor.
The invention adopts the above technical scheme, which has the following advantages：1. it is inverse to improve elastic wave least square
The imaging precision and resolution ratio of hour offset.2. the stability of elastic wave least square reverse-time migration can be improved.3. can weaken
Crosstalk effect between model parameter improves the imaging precision of weak responsive parameter.
The precision of the elastic wave least square reverse-time migration method proposed in the present invention is analyzed below by several examples
And stability.
As shown in Fig. 2, illustrating the advantage of the present invention by taking a groove model as an example first.Time step 1ms, between space
Every 10m, shot point (20 big gun) and geophone station are uniformly distributed in earth's surface.Focus is the Ricker wavelet of 15Hz, is added on direct stress.Fig. 3
For the migration result of different offset methods.The PP pictures (a) and PS pictures (b) of conventional reverse-time migration.Conventional least square reverse-time migration
RλAs (c), RμAs (d) and RρAs (e).The R of least square reverse-time migration based on regularization constraintλAs (f), RμAs (g) and
RρAs (h).The elastic wave least square reverse-time migration method based on regularization newly proposed can carry out subsurface structure accurate
Imaging, the resolution ratio higher of migration result, the continuity of lineups are more preferable.In addition, can to obtain preferable density anti-for new method
Coefficient section is penetrated, and the density reflectance factor section of conventional least square reverse-time migration method is poor.
Complicated Marmousi models (as shown in Figure 4) are used to test the offset method newly proposed below.Time
Step-length 1ms, space interval 10m, shot point (26 big gun) and geophone station are uniformly distributed in earth's surface.Focus is the Ricker wavelet of 15Hz, is added
On direct stress.Fig. 5 provides the imaging results of Marmousi model difference offset methods.Marmousi model difference offset methods
Imaging results.The PP pictures (a) and PS pictures (b) of conventional reverse-time migration.The R of conventional least square reverse-time migrationλAs (c) and RμPicture
(d).The R of least square reverse-time migration based on regularization constraintλAs (e) and RμAs (f).As seen from the figure, the least square inverse time is inclined
Shifting method is higher than the imaging precision of conventional reverse-time migration method.The elastic wave based on regularization constraint proposed in the present invention is minimum
Two to multiply reverse-time migration method imaging effect best.
The present invention is a kind of new seismic wave offset imaging method, can greatly improve multi component signal imaging precision and
Stability；The crosstalk effect between different parameters can be effectively suppressed, the imaging of weak responsive parameter (such as density reflectance factor) is improved
Effect；Reliable reflectance factor can be provided for following explanations in seismic prospecting and inverting, and then improve the identification of lithology and oil gas
Precision.
The above is only the better embodiment to the present invention, not makees limit in any form to the present invention
System, every any simple modification that embodiment of above is made according to the technical essence of the invention, equivalent variations and modification,
Belong in the range of technical solution of the present invention.
Claims (7)
1. the elastic wave least square reverse-time migration method based on regularization constraint, it is characterised in that follow the steps below：
(1) new object function is designed；
Fresh target function includes two：Simulated reflections wave and the difference and regularization term for observing back wave, the contribution of the two pass through
Regularization coefficient is adjusted；
(2) the elastic wave inverse migration operator and reflectance factor gradient formula under fresh target function are derived；
The gradient formula of elastic wave adjoint equation/inverse migration operator and object function to reflectance factor is derived based on adjoint method；
(3) gradient of reflectance factor is calculated；
It specifically includes：Source wavefield forward-propagating；Back wave residual error backpropagation；Forward and reverse wave field correlation obtains conventional ladder
Degree；Along with regularization term is to the gradient of reflectance factor；
(4) conjugate gradient method or quasi-Newton method inversion algorithm is used to handle gradient；
(5) iteration step length is sought using curve-parabola-fitting method；
(6) reflectivity model is updated, until meeting the condition of convergence.
2. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
Elastic wave least square reverse-time migration process is constrained using TV regularizations in the step (1), object function is：
Wherein：T is maximum time, and H is zoning, d1、d2And d3For regularization coefficient, β1、β2And β3For stability factor, Δ
vxWith Δ vzFor simulated reflections wave horizontal component and vertical component, Δ vx obsWith Δ vz obsFor observation back wave horizontal component and hang down
Straight component,
Wherein, λ and μ is Lame constants, and ρ is density, and equation (2) indicates the opposite variation of model parameter, dimensionless, can be used for
The size of reflectance factor is weighed, elastic wave least square reverse-time migration is exactly to seek optimal Rρ、RλAnd RμProcess, fresh target letter
Number includes two parts：Simulated reflections wave and the difference and regularization term for observing back wave, the contribution of the two pass through regularization coefficient
d1、d2And d3It adjusts.
3. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
Elastic wave inverse migration operator and reflectance factor gradient formula method in the step (2) under derivation fresh target function is as follows：
Elastic wave velocity-stress equation is：
Wherein, (vx,vz) it is Particle Vibration Velocity vector, (τxx,τzz,τxz) it is stress vector；
In elastic fluid, for background model parameters [λ, μ, ρ], background wave field [vx,vz,τxx,τzz,τxz] obtained by solving equation
It arrives, when there are model disturbance [Δ λ, Δ μ, Δ ρ], wave field knots modification is [Δ vx,Δvz,Δτxx,Δτzz,Δτxz], and it is full
Foot：
Abbreviation is simultaneously ignored high-order small quantity and is obtained：
To given parameter perturbation [Δ λ, Δ μ, Δ ρ], solves equation and obtain back wave [Δ vx,Δvz,Δτxx,Δτzz,Δ
τxz], the inverse migration process as in elastic fluid, in least square reverse-time migration, context parameter [λ, μ, ρ] is constant, background
Wave field is also constant, and the power of back wave is directly determined by parameter perturbation item；
By Rρ、RλAnd RμEquation is substituted into obtain：
Only consider the difference item of simulated reflections wave and observation back wave in object function：
Wherein, reflected wave field [the Δ v of simulationx,Δvz,Δτxx,Δτzz,Δτxz], this is solved about using method of Lagrange multipliers
Beam optimization problem, cost functional become：
Wherein,For Lagrange multiplier function,
Integration by parts obtains：
Wherein,
Enable θ J/ θ [Δ vx,Δvz,Δτxx,Δτzz,Δτxz]T=0 obtains corresponding adjoint equation, and form is as follows：
Object function is about the gradient formula of parameter perturbation：
In the case of TV regularizations, reflectance factor gradient formula becomes：
4. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
In the step (3)
A. it solves equation (3) and equation (6) obtains reflected wave field [Δ vx,Δvz,Δτxx,Δτzz,Δτxz]T, primary condition is：
B. it solves adjoint equation (12) and obtains backward extension wave fieldFinal value condition is：
C. the gradient by equation (14) calculating target function about reflectance factor.
5. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
Processing method is as follows in the step (4)：
Using L-BFGS methods：
Wherein, HkFor the approximate matrix that Hessian matrix is inverse, H is directly calculatedkNeed larger calculation amount, here by several groups of row to
Amount carrys out approximate Hk。
6. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
It is as follows that iteration step length method is sought in the step (5)：
Iteration step length is sought using Parabolic Fit
Wherein, α1And α2To sound out step-length, J1And J2For corresponding target function value, J0For the target function value of current iteration, meter
Calculate J1And J2Need four times additional forward modeling operations；
Then the optimum stepsize of current iteration is：
7. according to the elastic wave least square reverse-time migration method based on regularization constraint described in claim 1, it is characterised in that：
Reflectance factor is updated by following formula in the step (6)：
Wherein, mkAnd mk+1The respectively model parameter of current iteration and next iteration：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810042691.XA CN108333628B (en) | 2018-01-17 | 2018-01-17 | Elastic wave least square reverse-time migration method based on regularization constraint |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810042691.XA CN108333628B (en) | 2018-01-17 | 2018-01-17 | Elastic wave least square reverse-time migration method based on regularization constraint |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108333628A true CN108333628A (en) | 2018-07-27 |
| CN108333628B CN108333628B (en) | 2019-09-03 |
Family
ID=62925121
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810042691.XA Active CN108333628B (en) | 2018-01-17 | 2018-01-17 | Elastic wave least square reverse-time migration method based on regularization constraint |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108333628B (en) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108845355A (en) * | 2018-09-26 | 2018-11-20 | 中国矿业大学(北京) | Seismic migration imaging method and device |
| CN109471053A (en) * | 2018-10-18 | 2019-03-15 | 电子科技大学 | A kind of dielectric property iteration imaging method based on double constraints |
| CN109521469A (en) * | 2018-11-16 | 2019-03-26 | 国家海洋局第海洋研究所 | A kind of regularization inversion method of bottom sediment elastic parameter |
| CN109946741A (en) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | Pure qP wave least square reverse-time migration imaging method in a kind of TTI medium |
| CN110879415A (en) * | 2018-09-06 | 2020-03-13 | 中国石油化工股份有限公司 | Sticky sound reverse time migration method and system based on wave field decomposition |
| CN110888166A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Least square offset imaging method and device based on L-BFGS algorithm |
| CN111596346A (en) * | 2019-02-20 | 2020-08-28 | 中国石油天然气集团有限公司 | Elastic wave velocity inversion method and device |
| CN112130199A (en) * | 2020-07-31 | 2020-12-25 | 西安工程大学 | Optimized least square reverse time migration imaging method |
| CN113126149A (en) * | 2018-12-28 | 2021-07-16 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
| CN113552631A (en) * | 2021-08-16 | 2021-10-26 | 中煤科工集团西安研究院有限公司 | Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals |
| US11733413B2 (en) | 2021-04-30 | 2023-08-22 | Saudi Arabian Oil Company | Method and system for super resolution least-squares reverse time migration |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2796899A2 (en) * | 2013-04-23 | 2014-10-29 | CGG Services SA | Seismic data processing and apparatus |
| CN104216011A (en) * | 2013-06-05 | 2014-12-17 | 上海青凤致远地球物理地质勘探科技有限公司 | Reverse time migration method of stable qP wave in TTI (tilted transversely isotropic) media |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
-
2018
- 2018-01-17 CN CN201810042691.XA patent/CN108333628B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2796899A2 (en) * | 2013-04-23 | 2014-10-29 | CGG Services SA | Seismic data processing and apparatus |
| CN104216011A (en) * | 2013-06-05 | 2014-12-17 | 上海青凤致远地球物理地质勘探科技有限公司 | Reverse time migration method of stable qP wave in TTI (tilted transversely isotropic) media |
| CN106970416A (en) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | Elastic wave least square reverse-time migration system and method based on wave field separation |
Non-Patent Citations (3)
| Title |
|---|
| ZHIMING REN ET AL.: "Elastic full-waveform inversion using the second-generation wavelet and an adaptive-operator-length scheme", 《GEOPHYSICS》 * |
| ZHIMING REN ET AL.: "Least-squares reverse time migration in elastic media", 《GEOPHYSICAL JOURNAL INTERNATIONAL》 * |
| 王华忠 等: "最小二乘叠前深度偏移成像理论与方法", 《石油物探》 * |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110879415A (en) * | 2018-09-06 | 2020-03-13 | 中国石油化工股份有限公司 | Sticky sound reverse time migration method and system based on wave field decomposition |
| CN110888166A (en) * | 2018-09-10 | 2020-03-17 | 中国石油化工股份有限公司 | Least square offset imaging method and device based on L-BFGS algorithm |
| CN108845355A (en) * | 2018-09-26 | 2018-11-20 | 中国矿业大学(北京) | Seismic migration imaging method and device |
| CN109471053A (en) * | 2018-10-18 | 2019-03-15 | 电子科技大学 | A kind of dielectric property iteration imaging method based on double constraints |
| CN109521469A (en) * | 2018-11-16 | 2019-03-26 | 国家海洋局第海洋研究所 | A kind of regularization inversion method of bottom sediment elastic parameter |
| CN113126149A (en) * | 2018-12-28 | 2021-07-16 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
| CN113126149B (en) * | 2018-12-28 | 2024-04-09 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
| CN111596346A (en) * | 2019-02-20 | 2020-08-28 | 中国石油天然气集团有限公司 | Elastic wave velocity inversion method and device |
| CN109946741A (en) * | 2019-03-29 | 2019-06-28 | 中国石油大学(华东) | Pure qP wave least square reverse-time migration imaging method in a kind of TTI medium |
| CN112130199A (en) * | 2020-07-31 | 2020-12-25 | 西安工程大学 | Optimized least square reverse time migration imaging method |
| US11733413B2 (en) | 2021-04-30 | 2023-08-22 | Saudi Arabian Oil Company | Method and system for super resolution least-squares reverse time migration |
| CN113552631A (en) * | 2021-08-16 | 2021-10-26 | 中煤科工集团西安研究院有限公司 | Time-frequency double-domain regularization sparse deconvolution method and device for narrow-band signals |
| CN113552631B (en) * | 2021-08-16 | 2023-11-03 | 中煤科工集团西安研究院有限公司 | Time-frequency double-domain regularized sparse deconvolution method and device for narrowband signals |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108333628B (en) | 2019-09-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108333628B (en) | Elastic wave least square reverse-time migration method based on regularization constraint | |
| Liu et al. | Migration velocity analysis: Theory and an iterative algorithm | |
| US6826484B2 (en) | 3D prestack time migration method | |
| US10088588B2 (en) | Device and method for stable least-squares reverse time migration | |
| US20090257308A1 (en) | Migration velocity analysis methods | |
| EP3094992B1 (en) | Velocity model building for seismic data processing using pp-ps tomography with co-depthing constraint | |
| CN108873066A (en) | Elastic fluid fluctuates equation back wave Travel Time Inversion method | |
| Wang et al. | Inversion of seismic refraction and reflection data for building long-wavelength velocity models | |
| CN102841375A (en) | Method for tomography velocity inversion based on angle domain common imaging gathers under complicated condition | |
| CN107817526B (en) | Prestack seismic gather segmented amplitude energy compensation method and system | |
| US20100118652A1 (en) | Determination of depth moveout and of residual radii of curvature in the common angle domain | |
| EP2548052B1 (en) | System and method of 3d salt flank vsp imaging with transmitted waves | |
| WO2011124532A1 (en) | A process for characterising the evolution of a reservoir | |
| CN113740901B (en) | Land seismic data full-waveform inversion method and device based on complex undulating surface | |
| EP3067718A1 (en) | Boundary layer tomography method and device | |
| CN110187382A (en) | A kind of diving Wave and back wave wave equation Travel Time Inversion method | |
| US20180059276A1 (en) | System and method for focusing seismic images | |
| Hole et al. | Interface inversion using broadside seismic refraction data and three‐dimensional travel time calculations | |
| Liu | Migration velocity analysis | |
| DUAN et al. | Computation of density perturbation based point spread function and simulation of migration image in 3D depth domain | |
| Gibson Jr et al. | Modeling and velocity analysis with a wavefront-construction algorithm for anisotropic media | |
| Fei et al. | 3D common-reflection-point-based seismic migration velocity analysis | |
| Xiong et al. | An improved constrained velocity inversion algorithm for geological structures | |
| Barnes et al. | Diving wave tomography: A robust method for velocity estimation in a foothills geological context | |
| Kim et al. | Migration velocity analysis with the Kirchhoff integral |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |