Detailed description of the invention
As it is shown in figure 1, elastic vector reverse-time migration method, step is as follows:
Step 1: according to dielectric model, given source wavelet, it is achieved wave field just pushes out, utilizes decoupling extension equation method to obtain the longitudinal and transverse wave field of decoupling, and concrete grammar is as follows:
(1), according to dielectric model set in advance, load source wavelet, utilize elastic wave propagation operator, it is achieved the forward time continuation of seismic wave field, build source vector wave field:
Utilize following one-order velocity-stress equations for elastic waves
Building source wavefield and just pushing out staggering mesh finite-difference operator, wherein ρ is density, and C is elastic fluid stiffness matrix, v=(vx,vy,vz)TRepresenting Particle Vibration Velocity vector field, superscript notation " T " represents transposition, vxRepresent Particle Vibration Velocity vector field component in the x-direction, vyRepresent the Particle Vibration Velocity vector field component along y, vzRepresent the Particle Vibration Velocity vector field component along z, τ=(σxx,σyy,σzz,τyz,τxz,τxy)TIt is stress tensor, σxx、σyyAnd σzzIt is direct stress, τyz、τxzAnd τxyIt is shearing stress,Represent Particle Vibration Velocity vector field derivative on time orientation,Representing stress tensor derivative on time orientation, L is differential matrix,
Wherein, lx, lyAnd lzRepresenting the derivative along x, y and z direction respectively, in isotropic medium, stiffness matrix C is expressed as
Wherein λ and μ is Lame Coefficient.Above-mentioned wave equation is carried out discretization and obtains following elastic wave continuation operator:
Wherein, τSRepresent the discrete stress field of source wavefield, vSRepresenting the discrete Particle Vibration Velocity field of source wavefield, η is border absorptance, absorptance η=0 in target area, absorptance η=200 (0.5-0.5cos (π r/R)), r=1,2 in absorption region, border, ..., R, R is the thickness of absorbed layer, π represents that pi, Δ t are time sampling interval, and n Δ t represents whole time point, (n+1/2) Δ t is half timing node, n=1,2, ..., N, T0=N Δ t represents that total earthquake record receives duration, DfAnd DbRepresenting high-order staggering mesh finite-difference matrix operator respectively, expression is
With
WhereinWithRepresent forwardly and rearwardly Difference Schemes with Staggered in the x-direction,WithRepresent forwardly and rearwardly Difference Schemes with Staggered in the y-direction, whereinWithRepresenting forwardly and rearwardly Difference Schemes with Staggered in the z-direction, concrete form is as follows:
Wherein, Δ x, Δ y and the Δ z respectively sampling interval along x, y and z direction,Being 2M rank staggering mesh finite-difference coefficients, (i, j k) are space networks lattice point (i, j, smooth function k) to f.
(2), in source wavefield forward time continuation process, the P-wave And S data of the source wavefield of decoupling extension equation method acquisition decoupling are utilized:
Wave field extrapolation process utilizes decoupling extension equation method obtain the P-wave And S field data of decoupling, first passes through formula (3) and calculate source wavefield compressional wave stress wave field:
WhereinSource wavefield compressional wave stress field,Represent compressional wave stress wave field derivative on time orientation,Represent divergence operator, obtain the Particle Vibration Velocity field of focus compressional wave again through compressional wave stress field:
WhereinRepresent source wavefield longitudinal-wave particle vibration velocity field,The Particle Vibration Velocity field of expression source wavefield compressional wave derivative on time orientation,Represent gradient operator,The Particle Vibration Velocity vector field of expression source wavefield compressional wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of source wavefield compressional wave,Represent the Particle Vibration Velocity vector field component along z of source wavefield compressional wave, finally from the total wave field v of focusSIn deduct source wavefield longitudinal wave fieldObtain source wavefield shear wave field, namely
WhereinRepresent the Particle Vibration Velocity field of source wavefield shear wave,The Particle Vibration Velocity vector field of expression source wavefield shear wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of source wavefield shear wave,Represent the Particle Vibration Velocity vector field component along z of source wavefield shear wave.Outside description word describes, can also embodying the present invention in Figure of description intuitively, accompanying drawing 3-1 is the focus longitudinal wave field that obtains after the elastic wave field separation that rate pattern 5.8 kms in the x direction to the 12.4 km regions that 2-1 and Fig. 2 with reference to the accompanying drawings-2 is given are calculated component in the x-direction;Accompanying drawing 3-2 is the focus longitudinal wave field that obtains after the elastic wave field separation that rate pattern 5.8 kms in the x direction to the 12.4 km regions that 2-1 and Fig. 2 with reference to the accompanying drawings-2 is given are calculated component in the z-direction;Accompanying drawing 3-3 is the focus shear wave field that obtains after the elastic wave field that rate pattern 5.8 kms in the x direction to the 12.4 km regions that 2-1 and Fig. 2 with reference to the accompanying drawings-2 is given are calculated separates component in the x-direction, and accompanying drawing 3-4 is the focus shear wave field that obtains after the elastic wave field separation that rate pattern 5.8 kms in the x direction to the 12.4 km regions that 2-1 and Fig. 2-2 with reference to the accompanying drawings is given are calculated component in the z-direction.From accompanying drawing 3-1 to accompanying drawing 3-4 it can clearly be seen that wave field is separated comes.
Step 2: according to dielectric model, using multi-component seismic data as boundary value condition, it is achieved the inverse time extrapolation of seismic wave field, utilizes decoupling extension equation method to obtain the longitudinal and transverse wave field of decoupling, and concrete grammar is as follows:
(1), according to dielectric model set in advance, using the multi-component seismic data of surface seismic records as boundary value condition, utilize elastic wave propagation operator, it is achieved the inverse time continuation of seismic wave field, build detection vector wave field:
Multi-component earthquake data is used as boundary value condition, is loaded in following inverse time extrapolation equation
Wherein, τRRepresent the discrete stress field of detection wave field, vRRepresenting the discrete Particle Vibration Velocity field of detection wave field, η is border absorptance, absorptance η=0 in target area, absorptance η=200 (0.5-0.5cos (π r/R)), r=1,2 in absorption region, border, ..., R, R is the thickness of absorbed layer, π represents that pi, Δ t are time sampling interval, and n Δ t represents whole time point, (n+1/2) Δ t is half timing node, n=1,2, ..., N, T0=N Δ t represents that total earthquake record receives duration, DfAnd DbRepresenting high-order staggering mesh finite-difference matrix operator respectively, expression is
With
WhereinWithRepresent forwardly and rearwardly Difference Schemes with Staggered in the x-direction,WithRepresent forwardly and rearwardly Difference Schemes with Staggered in the y-direction, whereinWithRepresenting forwardly and rearwardly Difference Schemes with Staggered in the z-direction, concrete form is as follows:
Wherein, Δ x, Δ y and the Δ z respectively sampling interval along x, y and z direction,Being 2M rank staggering mesh finite-difference coefficients, (i, j k) are space networks lattice point (i, j, smooth function k) to f.
(2), in detection wave field inverse time continuation process, the P-wave And S data of the detection wave field of decoupling extension equation method acquisition decoupling are utilized:
Wave field extrapolation process utilizes decoupling extension equation method obtain the P-wave And S field data of decoupling, first passes through formula (7) and calculate detection wave field compressional wave stress field:
WhereinDetection wave field compressional wave stress field,Represent detection wave field compressional wave stress field derivative on time orientation,Represent divergence operator, obtain the Particle Vibration Velocity field of detection wave field compressional wave again through detection wave field compressional wave stress field:
WhereinRepresent detection wave field longitudinal-wave particle vibration velocity field,The Particle Vibration Velocity field of expression detection wave field compressional wave derivative on time orientation,Represent gradient operator,The Particle Vibration Velocity vector field of expression detection wave field compressional wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of detection wave field compressional wave,Represent the Particle Vibration Velocity vector field component along z of detection wave field compressional wave, finally from the total wave field v of detection wave fieldRIn deduct detection wave field longitudinal wave fieldObtain detection wave field shear wave field, namely
WhereinRepresent the Particle Vibration Velocity field of detection wave field shear wave,The Particle Vibration Velocity vector field of expression detection wave field shear wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of detection wave field shear wave,Represent the Particle Vibration Velocity vector field component along z of detection wave field shear wave.Accompanying drawing 4-1 is detection wave field by solving vector longitudinal wave field that decoupling extension equation method obtains component in the x-direction;The vector longitudinal wave field that accompanying drawing 4-2 is detection wave field to be obtained by decoupling extension equation method component in the z-direction;The vector shear wave field that accompanying drawing 4-3 is detection wave field to be obtained by decoupling extension equation method component in the x-direction, the vector shear wave field that accompanying drawing 4-4 is detection wave field to be obtained by decoupling extension equation method component in the z-direction.From accompanying drawing 4-1 to accompanying drawing 4-4 it can clearly be seen that wave field is separated comes.
Step 3: utilizing the long-pending cross-correlation image-forming condition of vector wave site of source wavefield illumination compensation to carry out imaging, it is thus achieved that the imaging results of scalar, concrete grammar is as follows:
The detection wave field decoupling P-wave And S that source wavefield decoupling P-wave And S that step 1 obtains and step 2 obtain is utilized to carry out following
Vector dot product cross-correlation imaging, and carry out focus illumination compensation, obtain multi-component seismic data elasticity vector reverse-time migration scalar imaging results.Wherein IPPRepresent PP imaging results, IPSIt is PS imaging results, ISPIt is SP imaging results, ISSIt is SS imaging results, t express time, T0Represent that earthquake record receives duration,Represent source wavefield longitudinal-wave particle vibration velocity field,The Particle Vibration Velocity vector field of expression source wavefield compressional wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of source wavefield compressional wave,Represent the Particle Vibration Velocity vector field component along z of source wavefield compressional wave,Represent detection wave field longitudinal-wave particle vibration velocity field,The Particle Vibration Velocity vector field of expression detection wave field compressional wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of detection wave field compressional wave,Represent the Particle Vibration Velocity vector field component along z of detection wave field compressional wave,Represent the Particle Vibration Velocity field of source wavefield shear wave,The Particle Vibration Velocity vector field of expression source wavefield shear wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of source wavefield shear wave,Represent the Particle Vibration Velocity vector field component along z of source wavefield shear wave,Represent the Particle Vibration Velocity field of detection wave field shear wave,The Particle Vibration Velocity vector field of expression detection wave field shear wave component in the x-direction,Represent the Particle Vibration Velocity vector field component along y of detection wave field shear wave,Representing the Particle Vibration Velocity vector field component along z of detection wave field shear wave, symbol represents vector dot product.Outside description word describes, it is also possible to embody the present invention in Figure of description intuitively, Fig. 5-1 is the PP imaging section that the elastic wave cross-correlation image-forming condition utilizing the vector dot product based on illumination compensation to build obtains;Fig. 5-2 is the PS imaging section that the elastic wave cross-correlation image-forming condition utilizing the vector dot product based on illumination compensation to build obtains;Fig. 5-3 is the SP imaging section that the elastic wave cross-correlation image-forming condition utilizing the vector dot product based on illumination compensation to build obtains;Fig. 5-4 is the SS imaging section that the elastic wave cross-correlation image-forming condition utilizing the vector dot product based on illumination compensation to build obtains.
Step 4: the imaging results obtained is carried out low frequency noise compacting, and specific implementation method is as follows:
Utilize following Laplce's filtering algorithm
It is filtered imaging section processing, obtains final imaging results.WhereinRepresent Laplace operator,WithIt is by the filtered PP of Laplce, PS, SP and SS imaging section respectively.Fig. 6-1 is final superposition the PP imaging section carrying out Laplce's filtering.Fig. 6-2 is final superposition the PS imaging section carrying out Laplce's filtering, and Fig. 6-3 is final superposition the SP imaging section carrying out Laplce's filtering, and Fig. 6-4 is final superposition the SS imaging section carrying out Laplce's filtering.