CN112684505B - Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function - Google Patents

Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function Download PDF

Info

Publication number
CN112684505B
CN112684505B CN202011418564.9A CN202011418564A CN112684505B CN 112684505 B CN112684505 B CN 112684505B CN 202011418564 A CN202011418564 A CN 202011418564A CN 112684505 B CN112684505 B CN 112684505B
Authority
CN
China
Prior art keywords
wave velocity
longitudinal
transverse
envelope
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011418564.9A
Other languages
Chinese (zh)
Other versions
CN112684505A (en
Inventor
罗静蕊
卢金波
张爽
王璐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Technology
Original Assignee
Xian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Technology filed Critical Xian University of Technology
Priority to CN202011418564.9A priority Critical patent/CN112684505B/en
Publication of CN112684505A publication Critical patent/CN112684505A/en
Application granted granted Critical
Publication of CN112684505B publication Critical patent/CN112684505B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Geophysics And Detection Of Objects (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

The invention disclosesThe elastic wave full waveform inversion method adopting the direct envelope sensitive kernel function comprises the following steps: s1, setting shot points and wave detection points, and acquiring transverse components and longitudinal components of observed data; s2, constructing a rectangular grid geological model, and setting a space sampling interval, a time sampling interval dt and a maximum sampling time Nt of forward modeling; s3, giving a longitudinal wave speed model v p (x) And transverse wave velocity model v s (x) And an initial longitudinal wave velocity model v p0 (x) And an initial shear wave velocity model v s0 (x) Given an objective function J to be optimized in the direct envelope inversion stage e (v p (x),v s (x) A) is provided; s4, aiming at an objective function J e (v p (x),v s (x) Optimized to obtain long wavelength component v of longitudinal wave velocity model pl (x) And long wavelength component v of transverse wave velocity model sl (x) The method comprises the steps of carrying out a first treatment on the surface of the S5, giving a new objective function J (v p (x),v s (x) A) is provided; s6, pair J (v) p (x),v s (x) Optimizing to obtain fine structure v of longitudinal wave velocity model pf (x) And fine structure v of transverse wave velocity model sf (x) A. The invention relates to a method for producing a fibre-reinforced plastic composite The method of the invention can use the new derivative of Fre chet in the elastic wave direct envelope to carry out energy propagation.

Description

Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function
Technical Field
The invention belongs to the technical field of geophysical exploration, and particularly relates to an elastic wave full waveform inversion method adopting a direct envelope sensitive kernel function.
Background
Oil and natural gas are important strategic resources for national economic development and national security, and as the development of the oil industry and the exploration and development continue to go deep, the exploration difficulty gradually increases from the construction exploration stage to the lithologic exploration stage. To meet this requirement, full waveform inversion is rapidly evolving and is a research hotspot in the geophysical world today. Full waveform inversion is an effective method for estimating formation parameters, and not only can provide a high-resolution high-precision velocity model for seismic imaging, but also has the capability of simultaneously inverting various different parameters.
Traditional full waveform inversion is based on weak scattering theory assumptions, resulting in a very strong dependence of the method on the initial model. If the initial model differs too much from the real model, local extrema are easily trapped. To overcome this problem, some scholars choose to employ a global optimization approach. However, global optimization requires searching through the whole model space, and finally determining an optimal model parameter, which is quite computationally intensive. Although algorithms have been improved in recent years, huge computational effort is still a major problem that hampers the development of global optimization, so that this type of approach cannot be widely used in high-dimensional problems. The center of gravity of the full waveform inversion remains the local optimization method. In this field, scholars have studied and developed various methods such as a multi-scale inversion method, a tomographic full waveform inversion method, an adaptive full waveform inversion method, and the like. These methods can alleviate the local extremum problem to some extent by using different wavefield information to construct the objective function to reduce the degree of nonlinearity of the inversion. However, these methods are still based on weak scattering assumptions, which lose their efficacy when strong scattering perturbations, such as a large scale salt dome model, are included in the model parameters. For large scale salt dome models, the industry typically uses a strategy of imaging, picking up boundaries, and filling the salt dome to construct the salt dome, however this technique requires significant manual intervention. In recent years, new techniques have been proposed to estimate the speed parameters inside the salt dome. The full variation method is introduced to carry out constraint in the self-adaptive full waveform inversion, and the optimal transmission distance is defined in the inversion, so that good effects are achieved.
The method aims at the inversion of the acoustic wave equation, the actual stratum is a complex elastic medium, and only longitudinal wave components can be simulated by using the acoustic wave equation, so that transverse wave information cannot be reflected. The wave field propagation condition in the stratum can be more accurately simulated by adopting an elastic wave equation, so that elastic wave inversion, particularly elastic wave inversion of a strong scattering medium, is necessary.
Disclosure of Invention
The invention aims to provide an elastic wave full waveform inversion method adopting a direct envelope sensitive kernel function, which can use a new Frechet derivative in the direct envelope to conduct energy propagation and realize simultaneous inversion of longitudinal wave speed and transverse wave speed aiming at an elastic medium.
The technical scheme adopted by the invention is that an elastic wave full waveform inversion method of a direct envelope sensitive kernel function is adopted, and the method is implemented according to the following steps:
step 1, setting shot points and wave detection points, and acquiring transverse components of observed data
Figure BDA0002821151500000021
And longitudinal component
Figure BDA0002821151500000022
Wherein t represents a time variable; x is x r ,x s Respectively representing the positions of the detector points and the shot points;
step 2, constructing a rectangular grid geological model, setting the number Nx of transverse grids and the number Nz of longitudinal grids of the model, and setting a space sampling interval, a time sampling interval dt and a maximum sampling time Nt of forward modeling, wherein the space sampling interval comprises a transverse sampling interval dx and a longitudinal sampling interval dz;
step 3, giving a longitudinal wave velocity model v p (x) And transverse wave velocity model v s (x) And an initial longitudinal wave velocity model v p0 (x) And an initial shear wave velocity model v s0 (x) Given an objective function J to be optimized in the direct envelope inversion stage e (v p (x),v s (x))；
Step 4, the objective function J is subjected to direct envelope inversion and iterative optimization algorithm e (v p (x),v s (x) Optimized to obtain long wavelength component v of longitudinal wave velocity model pl (x) And long wavelength component v of transverse wave velocity model sl (x)；
Step 5, giving a new objective function J (v p (x),v s (x))；
Step 6, for the objective function J (v p (x),v s (x) Optimizing to obtain fine structure v of longitudinal wave velocity model pf (x) And fine structure v of transverse wave velocity model sf (x)。
The present invention is also characterized in that,
in step 3, an objective function J is set e (v p (x),v s (x) Is a two-norm of the residual between the observed data envelope and the calculated data envelope, having the form:
Figure BDA0002821151500000031
in the formula (1), the components are as follows,
Figure BDA0002821151500000032
and->
Figure BDA0002821151500000033
Representing the transverse components of the observed wavefield, respectively>
Figure BDA0002821151500000034
And calculating the transverse wave field component>
Figure BDA0002821151500000035
Envelope of->
Figure BDA0002821151500000036
And->
Figure BDA0002821151500000037
Representing the longitudinal components of the observed wavefield, respectively->
Figure BDA0002821151500000038
And calculating the wave field longitudinal component +.>
Figure BDA0002821151500000039
Wherein the transverse and longitudinal components of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time; />
Figure BDA00028211515000000310
Representing the summation operation with respect to the shot and the detector.
In step 3, the meterCalculating transverse components of wave fields
Figure BDA00028211515000000311
And longitudinal component->
Figure BDA00028211515000000312
Can be obtained by the following steps:
first, solving an elastic wave equation:
Figure BDA0002821151500000041
in the formula (2), ρ (x) is the medium density at the spatial position x, and is set to be constant; λ (x) and μ (x) are the pull Mei Can values at spatial position x; f (f) x (t) and f z (t) a transverse source function and a longitudinal source function, respectively; delta (x-x) s ) As a pulse function, when x=x s When the function value is 1, otherwise, the function value is 0, which means that the source function is located at x s A place; pull Mei Canshu lambda (x) and mu (x) vs. longitudinal wave velocity v p (x) And transverse wave velocity v s (x) The relationship of (2) is as follows:
Figure BDA0002821151500000042
u in formula (2) x (t,x；x s ) And u z (t,x；x s ) Respectively when the shot point is positioned at x s At any spatial position x, the transverse and longitudinal components of the wavefield are such that x=x r Can obtain the transverse component of the calculated wave field
Figure BDA0002821151500000043
And longitudinal component
Figure BDA0002821151500000044
Namely:
Figure BDA0002821151500000045
in step 3, the envelope is obtained by hilbert transformation, and is defined as follows:
Figure BDA0002821151500000046
Figure BDA0002821151500000051
Figure BDA0002821151500000052
Figure BDA0002821151500000053
in the formulas (7) and (8), H { } represents a hilbert transform in which the parameter τ is a temporary variable introduced in performing the hilbert transform and represents a time offset.
The specific implementation mode of the step 4 is as follows:
direct envelope inversion calculates a new derivative of friechet directly based on energy scattering for the strong scattering problem:
Figure BDA0002821151500000054
in the formula (9), the amino acid sequence of the compound,
Figure BDA0002821151500000055
and->
Figure BDA0002821151500000056
Respectively represent the objective functions J e (v p (x),v s (x) Gradient and objective function J with respect to longitudinal wave velocity e (v p (x),v s (x) A gradient with respect to shear wave velocity; f (f) e Representing an envelope concomitant source, equal to an envelope residual;
Figure BDA0002821151500000057
and->
Figure BDA0002821151500000058
The envelope friechet derivatives, with respect to the longitudinal and transversal wave velocities, respectively, have the following form:
Figure BDA0002821151500000059
(10)
Figure BDA0002821151500000061
And->
Figure BDA0002821151500000062
Envelope virtual source function related to longitudinal wave velocity and transverse wave velocity, G e Green's function for envelope propagation; after the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by using the formula (9) and the formula (10), the velocity is iteratively updated by using the following criteria:
Figure BDA0002821151500000063
in the formula (11), the amino acid sequence of the compound,
Figure BDA0002821151500000064
and->
Figure BDA0002821151500000065
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; />
Figure BDA0002821151500000066
And->
Figure BDA0002821151500000067
For the longitudinal wave velocity step calculated according to the formulas (9) and (10) in the n-th iterationDegree and shear wave velocity gradients; />
Figure BDA0002821151500000068
And->
Figure BDA0002821151500000069
The longitudinal wave speed and the transverse wave speed of the n+1th step after updating;
the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by the method and iterated, and the long wavelength component v of the longitudinal wave velocity model is obtained after the iteration converges pl (x) And long wavelength component v of transverse wave velocity model sl (x)。
In step 5, v is obtained pl (x) And v sl (x) Thereafter, the next-stage objective function J (v p (x),v s (x) Defined as the two norms of the observed wavefield and the calculated wavefield residual, in the form:
Figure BDA00028211515000000610
。
in step 6, the objective function J (v p (x),v s (x) The velocity model gradient formula at this stage of optimization is as follows:
Figure BDA00028211515000000611
zeta in the formula (13) vp (x) And zeta vs (x) Respectively represent the objective functions J (v p (x),v s (x) Gradient and objective function J (v) with respect to longitudinal wave velocity p (x),v s (x) A gradient with respect to shear wave velocity; f represents the wavefield-concomitant source, equal to the wavefield residual; f (F) p And F p The wave field friechet derivatives, with respect to longitudinal and transverse wave velocities, respectively, have the following form:
Figure BDA0002821151500000071
q (14) p And Q s A virtual source function of the wave field with respect to longitudinal wave velocity and transverse wave velocity, respectively, G being a green's function of wave field propagation; after calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient by using the formulas (13) and (14), the longitudinal wave velocity and the transverse wave velocity are iteratively updated by the same method as the formula (11), that is:
Figure BDA0002821151500000072
in the formula (15), the amino acid sequence of the compound,
Figure BDA0002821151500000073
and->
Figure BDA0002821151500000074
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; (ζ) vp (x)) n And (ζ) vs (x)) n The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated according to the formula (13) and the formula (14) in the n-th iteration process; />
Figure BDA0002821151500000075
And->
Figure BDA0002821151500000076
The longitudinal wave speed and the transverse wave speed of the n+1th step after updating;
the method is used for calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient and iterating, and the fine structure v of the longitudinal wave velocity model can be obtained after the iteration converges pf (x) And fine structure v of transverse wave velocity model sf (x)。
The beneficial effects of the invention are as follows: firstly, calculating an observed wave field and an envelope of the calculated wave field, and establishing an envelope objective function; secondly, starting from a linear initial model, obtaining a salt dome longitudinal wave velocity long wavelength component and a transverse wave velocity long wavelength component by using a direct envelope inversion method; thirdly, establishing a wave field objective function according to the calculated wave field and the observed wave field; finally, based on the long wavelength velocity component, a fine velocity structure is obtained by utilizing the conventional elastic wave full waveform inversion.
Compared with the common inversion method, the method is more suitable for inversion of the strong scattering medium parameters, and can accurately estimate the longitudinal wave speed and the transverse wave speed through an elastic wave equation. The reason for the above advantages is that: firstly, because a direct envelope inversion method is used, a new Fre chet derivative is adopted for energy propagation, and weak scattering approximation of the traditional inversion is avoided; second, because an elastic wave equation is used, simultaneous inversion of multiple parameters can be achieved.
Drawings
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2 is a true longitudinal wave velocity model of the target zone;
FIG. 3 is a real transverse wave velocity model of the target zone;
FIG. 4 is a linear initial longitudinal wave velocity model;
FIG. 5 is a linear initial shear wave velocity model;
FIG. 6 is a longitudinal wave velocity model obtained by inversion using a conventional method;
FIG. 7 is a shear wave velocity model obtained using inversion of conventional methods;
FIG. 8 is a longitudinal velocity long wavelength component inverted using the method of the present invention;
FIG. 9 is a plot of the long wavelength components of shear wave velocity inverted using the method of the present invention;
FIG. 10 is a fine longitudinal velocity model inverted using the method of the present invention;
FIG. 11 is a model of the velocity of a fine shear wave inverted using the method of the present invention.
Detailed Description
The invention will be described in detail below with reference to the drawings and the detailed description.
The invention adopts an elastic wave full waveform inversion method of a direct envelope sensitive kernel function, as shown in figure 1, and comprises the following specific steps:
step 1, setting a gunThe point and the wave detection point acquire the original seismic data to obtain the transverse component of the observed data
Figure BDA0002821151500000091
And longitudinal component->
Figure BDA0002821151500000092
Wherein t represents a time variable, x r ,x s And the positions of the detection points and the shot points are respectively used for carrying out conventional preprocessing on the acquired seismic data, including dynamic and static correction, noise elimination and other processing.
And 2, constructing a rectangular grid geological model, setting the number Nx of transverse grids and the number Nz of longitudinal grids of the model according to technical indexes and work area requirements, and setting the space sampling interval (comprising the transverse sampling interval dx and the longitudinal sampling interval dz), the time sampling interval dt and the maximum sampling time Nt of forward modeling. In practice, relevant parameters can be determined according to the effective frequency band range of the seismic data, the time length of the seismic recording and the seismic data observation system; the criteria for selecting the grid parameters are such that when finite difference forward modeling is performed based on the grid, not only stability conditions are met but also numerical dispersion is effectively suppressed.
Step 3, giving a longitudinal wave velocity model v p (x) And transverse wave velocity model v s (x) And an initial longitudinal wave velocity model v p0 (x) And an initial shear wave velocity model v s0 (x) Given an objective function J to be optimized in the direct envelope inversion stage e (v p (x),v s (x) A) is provided; longitudinal wave velocity model v p (x) And an initial model v of the transverse wave velocity model p0 (x) And v s0 (x) The method is a rough estimation of a real longitudinal wave velocity model and a real transverse wave velocity model, and can be generally set as a linear initial model with uniform transverse longitudinal linear increase; the objective function defines the degree of matching between the observed data and the calculated data, and is optimized to achieve updating of the model parameters. In the invention, an objective function J is set in a direct envelope inversion stage e (v p (x),v s (x) Is the second norm of the residual between the observed data envelope and the calculated data envelope,has the following form:
Figure BDA0002821151500000093
in the formula (1), the components are as follows,
Figure BDA0002821151500000094
and->
Figure BDA0002821151500000095
Representing the transverse components of the observed wavefield, respectively>
Figure BDA0002821151500000096
And calculating the transverse wave field component>
Figure BDA0002821151500000097
Envelope of->
Figure BDA0002821151500000098
And->
Figure BDA0002821151500000099
Representing the longitudinal components of the observed wavefield, respectively->
Figure BDA00028211515000000910
And calculating the wave field longitudinal component +.>
Figure BDA00028211515000000911
Wherein the transverse and longitudinal components of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time; />
Figure BDA00028211515000000912
Representing the summation operation with respect to the shot and the detector.
In step 3, the transverse component of the wavefield is calculated
Figure BDA0002821151500000101
And longitudinal component->
Figure BDA0002821151500000102
The method can be obtained by the following steps of firstly solving an elastic wave equation:
Figure BDA0002821151500000103
in the formula (2), ρ (x) is the medium density at the spatial position x, and is set to be constant; λ (x) and μ (x) are the pull Mei Can values at spatial position x; f (f) x (t) and f z (t) a transverse source function and a longitudinal source function, respectively; delta (x-x) s ) As a pulse function, when x=x s When the function value is 1, otherwise, the function value is 0, which means that the source function is located at x s A place; pull Mei Canshu lambda (x) and mu (x) vs. longitudinal wave velocity v p (x) And transverse wave velocity v s (x) The relationship of (2) is as follows:
Figure BDA0002821151500000104
u in formula (2) x (t,x；x s ) And u z (t,x；x s ) Respectively when the shot point is positioned at x s At any spatial position x, the transverse and longitudinal components of the wavefield are such that x=x r Can obtain the transverse component of the calculated wave field
Figure BDA0002821151500000105
And longitudinal component
Figure BDA0002821151500000106
Namely:
Figure BDA0002821151500000107
in step (3), the envelope may be obtained by hilbert transform, defined as follows:
Figure BDA0002821151500000111
Figure BDA0002821151500000112
Figure BDA0002821151500000113
Figure BDA0002821151500000114
in the formulas (7) and (8), H { } represents a hilbert transform in which the parameter τ is a temporary variable introduced in performing the hilbert transform and represents a time offset.
The objective function J e (v p (x),v s (x) Envelope information of the signal is used because the signal envelope contains abundant low-frequency components, which facilitates the solution of long-wavelength components of the model.
Step 4, the objective function J is subjected to direct envelope inversion and iterative optimization algorithm e (v p (x),v s (x) Optimized to obtain long wavelength component v of longitudinal wave velocity model pl (x) And long wavelength component v of transverse wave velocity model sl (x) The method comprises the steps of carrying out a first treatment on the surface of the The traditional inversion method is based on weak scattering approximation and uses a chain rule to calculate the derivative of an objective function with respect to model parameters, and the direct envelope inversion in the invention is used for directly calculating a new Frenchet derivative based on energy scattering aiming at the strong scattering problem:
Figure BDA0002821151500000121
in the formula (9), the amino acid sequence of the compound,
Figure BDA0002821151500000122
and->
Figure BDA0002821151500000123
Respectively represent the objective functions J e (v p (x),v s (x) Gradient and objective function J with respect to longitudinal wave velocity e (v p (x),v s (x) A gradient with respect to shear wave velocity; f (f) e Representing an envelope concomitant source, equal to an envelope residual;
Figure BDA0002821151500000124
and->
Figure BDA0002821151500000125
The envelope friechet derivatives, with respect to the longitudinal and transversal wave velocities, respectively, have the following form:
Figure BDA0002821151500000126
(10)
Figure BDA0002821151500000127
And->
Figure BDA0002821151500000128
Envelope virtual source functions for longitudinal wave velocity and transverse wave velocity, respectively, G e Is a green's function of envelope propagation. After the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by using the formula (9) and the formula (10), the velocity is iteratively updated by using the following criteria:
Figure BDA0002821151500000129
in the formula (11), the amino acid sequence of the compound,
Figure BDA00028211515000001210
and->
Figure BDA00028211515000001211
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; />
Figure BDA00028211515000001212
And->
Figure BDA00028211515000001213
The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated according to the formula (9) and the formula (10) in the n-th iteration process; />
Figure BDA00028211515000001214
And->
Figure BDA00028211515000001215
The longitudinal wave velocity and the transverse wave velocity of the n+1th step after updating.
The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by the method and iterated, and the long wavelength component v of the longitudinal wave velocity model can be obtained after the iteration converges pl (x) And long wavelength component v of transverse wave velocity model sl (x)。
Step 5, giving a new objective function J (v p (x),v s (x) A) is provided; obtain v pl (x) And v sl (x) Thereafter, the next-stage objective function J (v p (x),v s (x) Defined as the two norms of the observed wavefield and the calculated wavefield residual, in the form:
Figure BDA0002821151500000131
step 6, for the objective function J (v p (x),v s (x) Optimizing to obtain fine structure v of longitudinal wave velocity model pf (x) And fine structure v of transverse wave velocity model sf (x) The method comprises the steps of carrying out a first treatment on the surface of the And in the fine structure inversion stage, the model parameters are subjected to iterative optimization according to the objective function defined by the formula (12). For the objective function J (v p (x),v s (x) The velocity model gradient formula at this stage of optimization is as follows:
Figure BDA0002821151500000132
zeta in the formula (13) vp (x) And zeta vs (x) Respectively represent the objective functions J (v p (x),v s (x) Gradient and objective function J (v) with respect to longitudinal wave velocity p (x),v s (x) A gradient with respect to shear wave velocity; f represents the wavefield-concomitant source, equal to the wavefield residual; f (F) p And F p The wave field friechet derivatives, with respect to longitudinal and transverse wave velocities, respectively, have the following form:
Figure BDA0002821151500000133
q (14) p And Q s G is the green's function of wave field propagation, which is the virtual source function of the wave field with respect to longitudinal and transverse wave velocities, respectively. After calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient by using the formulas (13) and (14), the longitudinal wave velocity and the transverse wave velocity are iteratively updated by the same method as the formula (11), that is:
Figure BDA0002821151500000141
in the formula (15), the amino acid sequence of the compound,
Figure BDA0002821151500000142
and->
Figure BDA0002821151500000143
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; (ζ) vp (x)) n And (ζ) vs (x)) n The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated according to the formula (13) and the formula (14) in the n-th iteration process; />
Figure BDA0002821151500000144
And->
Figure BDA0002821151500000145
The longitudinal wave speed and the transverse wave speed of the n+1th step after updating;
the method is used for calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient and iterating, and the fine structure v of the longitudinal wave velocity model can be obtained after the iteration converges pf (x) And fine structure v of transverse wave velocity model sf (x)。
Model calculation example
The implementation process of the invention is applied to a strong scattering salt dome geologic model. The model is 16100 meters in transverse direction and 3750 meters in longitudinal direction. The model contains a large range of high-speed salt domes, while the background medium exhibits a low-speed distribution, representing typical strong scattering properties. The longitudinal wave velocity model and the transverse wave velocity model of the target zone are shown in fig. 2 and 3, respectively.
The initial longitudinal wave velocity model and the initial transverse wave velocity model are set as linear models, as shown in fig. 4 and 5, respectively, in which the longitudinal wave velocity model gradually increases linearly from shallow to deep 1530m/s to 3930m/s, and the transverse wave velocity model gradually increases linearly from shallow to deep 950m/s to 2000kg/m 3 。
The longitudinal wave velocity model and the transverse wave velocity model obtained by the conventional elastic wave full waveform inversion method are shown in fig. 6 and fig. 7, respectively. It can be seen that the conventional method cannot accurately invert the salt dome structure because the initial velocity model is too far from the true velocity model.
Fig. 8 and 9 are long wavelength components of a longitudinal wave velocity model and a transverse wave velocity model of a target area obtained by inversion in steps 1 to 4 of the method according to the present invention. The accurate acquisition of long wavelength components, i.e. large scale background models, is critical for the whole parametric inversion during inversion, since it ensures that the iteration can converge to the correct globally optimal solution and further inversion of the subsequent fine structure. From the results shown in fig. 8 and 9, it can be seen that the method proposed in the present invention gives a very good profile of the salt dome and the velocity profile.
Based on the long wavelength velocity model obtained as described above, the fine structures of the longitudinal wave velocity model and the transverse wave velocity model obtained by inversion using the method proposed in the present invention are shown in fig. 10 and 11. It can be seen that the inversion results are close to the real model, illustrating the effectiveness of the method proposed in the present invention.

Claims (7)

1. The elastic wave full waveform inversion method adopting the direct envelope sensitive kernel function is characterized by comprising the following steps of:
step 1, setting shot points and wave detection points, and acquiring transverse components of observed data
Figure FDA0002821151490000011
And longitudinal component
Figure FDA0002821151490000012
Wherein t represents a time variable; x is x r ,x s Respectively representing the positions of the detector points and the shot points;
step 2, constructing a rectangular grid geological model, setting the number Nx of transverse grids and the number Nz of longitudinal grids of the model, and setting a space sampling interval, a time sampling interval dt and a maximum sampling time Nt of forward modeling, wherein the space sampling interval comprises a transverse sampling interval dx and a longitudinal sampling interval dz;
step 3, giving a longitudinal wave velocity model v p (x) And transverse wave velocity model v s (x) And an initial longitudinal wave velocity model v p0 (x) And an initial shear wave velocity model v s0 (x) Given an objective function J to be optimized in the direct envelope inversion stage e (v p (x),v s (x))；
Step 4, the objective function J is subjected to direct envelope inversion and iterative optimization algorithm e (v p (x),v s (x) Optimized to obtain long wavelength component v of longitudinal wave velocity model pl (x) And long wavelength component v of transverse wave velocity model sl (x)；
Step 5, giving a new objective function J (v p (x),v s (x))；
Step 6, for the objective function J (v p (x),v s (x) Optimizing to obtain fine structure v of longitudinal wave velocity model pf (x) And fine structure v of transverse wave velocity model sf (x)。
2. The method for full waveform inversion of elastic waves using direct envelope sensitivity kernel function as claimed in claim 1, wherein in step 3, an objective function J is set e (v p (x),v s (x) Is a two-norm of the residual between the observed data envelope and the calculated data envelope, having the form:
Figure FDA0002821151490000013
in the formula (1), the components are as follows,
Figure FDA0002821151490000021
and->
Figure FDA0002821151490000022
Representing the transverse components of the observed wavefield, respectively>
Figure FDA0002821151490000023
And calculating the transverse wave field component>
Figure FDA0002821151490000024
Envelope of->
Figure FDA0002821151490000025
And->
Figure FDA0002821151490000026
Representing the longitudinal components of the observed wavefield, respectively->
Figure FDA0002821151490000027
And calculating the wave field longitudinal component +.>
Figure FDA0002821151490000028
Wherein the transverse and longitudinal components of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time; />
Figure FDA0002821151490000029
Representing the summation operation with respect to the shot and the detector.
3. The elastic wave full waveform inversion method using direct envelope sensitivity kernel function as claimed in claim 2, wherein in step 3, the transverse component of the wave field is calculated
Figure FDA00028211514900000210
And longitudinal component->
Figure FDA00028211514900000211
Can be obtained by the following steps:
first, solving an elastic wave equation:
Figure FDA00028211514900000212
in the formula (2), ρ (x) is the medium density at the spatial position x, and is set to be constant; λ (x) and μ (x) are the pull Mei Can values at spatial position x; f (f) x (t) and f z (t) a transverse source function and a longitudinal source function, respectively; delta (x-x) s ) As a pulse function, when x=x s When the function value is 1, otherwise, the function value is 0, which means that the source function is located at x s A place; pull Mei Canshu lambda (x) and mu (x) vs. longitudinal wave velocity v p (x) And transverse wave velocity v s (x) The relationship of (2) is as follows:
Figure FDA00028211514900000213
u in formula (2) x (t,x；x s ) And u z (t,x；x s ) Respectively when the shot point is positioned at x s At any spatial position x, the transverse and longitudinal components of the wavefield are such that x=x r Can obtain the transverse component of the calculated wave field
Figure FDA0002821151490000031
And longitudinal component
Figure FDA0002821151490000032
Namely:
Figure FDA0002821151490000033
4. a method of full waveform inversion of elastic waves using a direct envelope sensitive kernel function as claimed in claim 3 wherein in step 3 the envelope is obtained by hilbert transform defined as follows:
Figure FDA0002821151490000034
Figure FDA0002821151490000035
Figure FDA0002821151490000036
Figure FDA0002821151490000037
in the formulas (7) and (8), H { } represents a hilbert transform in which the parameter τ is a temporary variable introduced in performing the hilbert transform and represents a time offset.
5. The method of full waveform inversion of elastic waves using a direct envelope sensitivity kernel of claim 4 wherein the embodiment of step 4 is as follows:
direct envelope inversion is used for calculating new directly based on energy scattering aiming at strong scattering problem
Figure FDA0002821151490000038
Derivative:
Figure FDA0002821151490000041
in the formula (9), the amino acid sequence of the compound,
Figure FDA0002821151490000042
and->
Figure FDA0002821151490000043
Respectively represent the objective functions J e (v p (x),v s (x) Gradient and objective function J with respect to longitudinal wave velocity e (v p (x),v s (x) A gradient with respect to shear wave velocity; f (f) e Representing an envelope concomitant source, equal to an envelope residual; />
Figure FDA0002821151490000044
And->
Figure FDA0002821151490000045
Envelope about longitudinal wave velocity and transverse wave velocity, respectively +.>
Figure FDA0002821151490000046
Derivative, having the form:
Figure FDA0002821151490000047
(10)
Figure FDA0002821151490000048
And->
Figure FDA0002821151490000049
Envelope virtual source functions for longitudinal wave velocity and transverse wave velocity, respectively, G e Green's function for envelope propagation; after the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by using the formula (9) and the formula (10), the velocity is iteratively updated by using the following criteria:
Figure FDA00028211514900000410
in the formula (11), the amino acid sequence of the compound,
Figure FDA00028211514900000411
and->
Figure FDA00028211514900000412
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; />
Figure FDA00028211514900000413
And->
Figure FDA00028211514900000414
The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated according to the formula (9) and the formula (10) in the n-th iteration process; />
Figure FDA00028211514900000415
And->
Figure FDA00028211514900000416
The longitudinal wave speed and the transverse wave speed of the n+1th step after updating;
calculating longitudinal wave velocity gradient by the methodAnd transverse wave velocity gradient and iterating to obtain long wavelength component v of longitudinal wave velocity model after iteration convergence pl (x) And long wavelength component v of transverse wave velocity model sl (x)。
6. The method for full waveform inversion of elastic waves using a direct envelope sensitivity kernel of claim 5, wherein in step 5, v is obtained pl (x) And v sl (x) Thereafter, the next-stage objective function J (v p (x),v s (x) Defined as the two norms of the observed wavefield and the calculated wavefield residual, in the form:
Figure FDA0002821151490000051
。
7. the method of full waveform inversion of elastic waves using direct envelope sensitivity kernel of claim 6 wherein in step 6, the objective function J (v p (x),v s (x) The velocity model gradient formula at this stage of optimization is as follows:
Figure FDA0002821151490000052
zeta in the formula (13) vp (x) And zeta vs (x) Respectively represent the objective functions J (v p (x),v s (x) Gradient and objective function J (v) with respect to longitudinal wave velocity p (x),v s (x) A gradient with respect to shear wave velocity; f represents the wavefield-concomitant source, equal to the wavefield residual; f (F) p And F p The wave field friechet derivatives, with respect to longitudinal and transverse wave velocities, respectively, have the following form:
Figure FDA0002821151490000053
q (14) p And Q s Respectively with respect to longitudinal wave velocity and transverse wave velocityA virtual source function of the wave field, G is a green function of wave field propagation; after calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient by using the formulas (13) and (14), the longitudinal wave velocity and the transverse wave velocity are iteratively updated by the same method as the formula (11), that is:
Figure FDA0002821151490000061
in the formula (15), the amino acid sequence of the compound,
Figure FDA0002821151490000062
and->
Figure FDA0002821151490000063
The longitudinal wave speed and the transverse wave speed of the n-th iteration are respectively; alpha n The iteration step length of the nth step is a positive constant; (ζ) vp (x)) n And (ζ) vs (x)) n The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated according to the formula (13) and the formula (14) in the n-th iteration process; />
Figure FDA0002821151490000064
And->
Figure FDA0002821151490000065
The longitudinal wave speed and the transverse wave speed of the n+1th step after updating;
the method is used for calculating the longitudinal wave velocity gradient and the transverse wave velocity gradient and iterating, and the fine structure v of the longitudinal wave velocity model can be obtained after the iteration converges pf (x) And fine structure v of transverse wave velocity model sf (x)。
CN202011418564.9A 2020-12-07 2020-12-07 Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function Active CN112684505B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011418564.9A CN112684505B (en) 2020-12-07 2020-12-07 Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011418564.9A CN112684505B (en) 2020-12-07 2020-12-07 Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function

Publications (2)

Publication Number Publication Date
CN112684505A CN112684505A (en) 2021-04-20
CN112684505B true CN112684505B (en) 2023-06-30

Family

ID=75447504

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011418564.9A Active CN112684505B (en) 2020-12-07 2020-12-07 Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function

Country Status (1)

Country Link
CN (1) CN112684505B (en)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007340A (en) * 2019-02-01 2019-07-12 西安理工大学 Salt dome speed density estimation method based on the direct envelope inverting of angle domain

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10459096B2 (en) * 2016-07-13 2019-10-29 Exxonmobil Upstream Research Company Joint full wavefield inversion of P-wave velocity and attenuation using an efficient first order optimization

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007340A (en) * 2019-02-01 2019-07-12 西安理工大学 Salt dome speed density estimation method based on the direct envelope inverting of angle domain

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
二维时间空间域弹性波全波形反演-理论模型数值试验;廖建平;刘和秀;戴世鑫;赵延林;Andrew Hursthouse;;地球物理学进展(第02期);全文 *

Also Published As

Publication number Publication date
CN112684505A (en) 2021-04-20

Similar Documents

Publication Publication Date Title
CN108873066B (en) Elastic medium wave equation reflected wave travel time inversion method
CN110007340B (en) Salt dome velocity density estimation method based on angle domain direct envelope inversion
CN111239819B (en) Direct envelope inversion method with polarity based on seismic channel attribute analysis
CN113341465B (en) Directional anisotropic medium ground stress prediction method, device, medium and equipment
CN110187382B (en) Traveling time inversion method for wave equation of reverse wave and reflected wave
CN111025387B (en) Pre-stack earthquake multi-parameter inversion method for shale reservoir
WO1994029750A1 (en) Method of determining earth elastic parameters in anisotropic media
CN109459789B (en) Time-domain full waveform inversion method based on amplitude decaying and linear interpolation
CN110888159B (en) Elastic wave full waveform inversion method based on angle decomposition and wave field separation
US10132945B2 (en) Method for obtaining estimates of a model parameter so as to characterise the evolution of a subsurface volume
CN113031068A (en) Reflection coefficient accurate base tracking prestack seismic inversion method
Waheed et al. An iterative fast sweeping based eikonal solver for tilted orthorhombic media
CN111999764B (en) Method for constructing least square reverse time migration under salt based on time-frequency domain objective function
CN111025388B (en) Multi-wave combined prestack waveform inversion method
CN112684505B (en) Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function
Warner et al. Full-elastic AVA extraction using acoustic FWI
Zhou et al. New numerical schemes for time-domain seismic wave modeling in viscoelastic media
CN111208568B (en) Time domain multi-scale full waveform inversion method and system
CN109901221B (en) Seismic data anisotropy modeling method based on dynamic correction velocity parameter
Saraswat et al. Simultaneous stochastic inversion of prestack seismic data using hybrid evolutionary algorithm
Li et al. Small-scale Fracture Detection Via Anisotropic Bayesian Ant-tracking Colony Optimization Driven By Azimuthal Seismic Data
CN112099079B (en) Adaptive frequency division series reflectivity inversion method and system
Zhou et al. Anisotropic model building with well control
CN114624766B (en) Elastic wave least square reverse time migration gradient solving method based on traveling wave separation
CN115016002A (en) Direct envelope inversion method including angular domain illumination compensation

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