CN112684505A - 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
CN112684505A
CN112684505A CN202011418564.9A CN202011418564A CN112684505A CN 112684505 A CN112684505 A CN 112684505A CN 202011418564 A CN202011418564 A CN 202011418564A CN 112684505 A CN112684505 A CN 112684505A
Authority
CN
China
Prior art keywords
longitudinal
transverse
velocity
wave
wave velocity
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
Application number
CN202011418564.9A
Other languages
Chinese (zh)
Other versions
CN112684505B (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 discloses an elastic wave full waveform inversion method adopting a direct envelope sensitive kernel function, which specifically comprises the following steps: s1, setting a shot point and a demodulator probe, and acquiring a transverse component and a longitudinal component of the observation data; s2, constructing a rectangular grid geological model, and setting a space sampling interval, a time sampling interval dt and maximum sampling time Nt of forward modeling; s3, setting longitudinal wave velocity model vp(x) And transverse wave velocity model vs(x) And an initial longitudinal wave velocity model vp0(x) And an initial shear velocity model vs0(x) Given an objective function J to be optimized in the direct envelope inversion stagee(vp(x),vs(x) ); s4, aiming at an objective function Je(vp(x),vs(x) Optimized to obtain the long wavelength component v of the longitudinal wave velocity modelpl(x) And the long wavelength component v of the transverse wave velocity modelsl(x) (ii) a S5, a new objective function J (v) to be optimized in the next stage is givenp(x),vs(x) ); s6, pair J (v)p(x),vs(x) For example) to perform an optimization,obtaining fine structure v of longitudinal wave velocity modelpf(x) And fine structure v of shear wave velocity modelsf(x) In that respect The method can use a new Frechet derivative 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
Petroleum and natural gas are important strategic resources related to national economic development and national safety, and with the development of the petroleum industry and the continuous deepening of exploration and development, the exploration difficulty gradually increases from a structural exploration stage to a lithological exploration stage. To comply with this requirement, full waveform inversion has rapidly developed and is a hot spot in the geophysical field today. Full waveform inversion is an effective method for stratum parameter estimation, can provide a high-resolution and high-precision velocity model for seismic imaging, and has the capability of simultaneously inverting various different parameters.
The traditional full waveform inversion is based on weak scattering theory assumption, so that the method has strong dependence on an initial model. If the initial model differs too much from the real model, it is prone to fall into local extrema. To overcome this problem, some scholars choose to employ a global optimization approach. However, global optimization requires searching in the whole model space, and finally determining an optimal model parameter, which is quite computationally expensive. Although there have been many improvements in algorithms in recent years, the huge amount of computation is still a major problem hindering the development of global optimization, so that this type of method cannot be widely applied to high-dimensional problems. The center of gravity of the full waveform inversion is still a local optimization method. In this field, researchers 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. The methods can alleviate the local extremum problem to a certain extent by using different wave field information to construct an objective function so as to reduce the nonlinearity degree of inversion. However, these methods are still based on the assumption of weak scattering, and lose their efficacy when strong scattering perturbations, such as large-scale salt dome models, are included in the model parameters. For large-scale salt dome models, the industry typically employs imaging, boundary picking, and salt dome filling strategies to construct the salt dome, however this technique requires a great deal of manual intervention. In recent years, new techniques have been proposed to estimate the velocity parameters inside the salt dome. If a total variation method is introduced to carry out constraint in the self-adaptive full waveform inversion and an optimal transmission distance is defined in the inversion, good effects are achieved.
The method aims at the inversion of the acoustic wave equation, the actual stratum is a complex elastic medium, and the acoustic wave equation can only simulate the longitudinal wave component and cannot reflect the transverse wave information. The elastic wave equation can be used for more accurately simulating the wave field propagation condition in the stratum, so that elastic wave inversion is necessary, particularly the elastic wave inversion of a strong scattering medium.
Disclosure of Invention
The invention aims to provide an elastic wave full waveform inversion method adopting a direct envelope sensitive kernel function, which can carry out energy propagation by using a new Frechet derivative in a direct envelope and realize simultaneous inversion of longitudinal wave velocity and transverse wave velocity 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 a shot point and a wave detection point, and acquiring transverse components of observation data
Figure BDA0002821151500000021
And a longitudinal component
Figure BDA0002821151500000022
Wherein t represents a time variable; x is the number ofr,xsRespectively representing the positions of the demodulator probe and the shot point;
step 2, constructing a rectangular grid geological model, setting the number Nx of transverse grid points and the number Nz of longitudinal grid points of the model, and setting a space sampling interval, a time sampling interval dt and maximum sampling time Nt of forward simulation, wherein the space sampling interval comprises a transverse sampling interval dx and a longitudinal sampling interval dz;
step 3, setting a longitudinal wave velocity model vp(x) And transverse wave velocity model vs(x) And an initial longitudinal wave velocity model vp0(x) And initial transverse wave velocityDegree model vs0(x) Given an objective function J to be optimized in the direct envelope inversion stagee(vp(x),vs(x))；
Step 4, carrying out direct envelope inversion and iterative optimization algorithm on the target function Je(vp(x),vs(x) Optimized to obtain the long wavelength component v of the longitudinal wave velocity modelpl(x) And the long wavelength component v of the transverse wave velocity modelsl(x)；
Step 5, a new objective function J (v) to be optimized in the next stage is givenp(x),vs(x))；
Step 6, aiming at the target function J (v)p(x),vs(x) ) is optimized to obtain a fine structure v of the longitudinal wave velocity modelpf(x) And fine structure v of shear wave velocity modelsf(x)。
The present invention is also characterized in that,
in step 3, an objective function J is sete(vp(x),vs(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 reaction mixture is,
Figure BDA0002821151500000032
and
Figure BDA0002821151500000033
respectively representing the transverse components of the observed wavefield
Figure BDA0002821151500000034
And calculating the wavefield transverse component
Figure BDA0002821151500000035
The envelope of the (c) signal is,
Figure BDA0002821151500000036
and
Figure BDA0002821151500000037
respectively representing the longitudinal components of the observed wavefield
Figure BDA0002821151500000038
And calculating the wavefield longitudinal component
Figure BDA0002821151500000039
Wherein the transverse component and the longitudinal component of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time;
Figure BDA00028211515000000310
indicating that the summation operation is performed with respect to the shot point and the demodulator probe.
In step 3, the transverse component of the wave field is calculated
Figure BDA00028211515000000311
And a longitudinal component
Figure BDA00028211515000000312
Can be obtained by the following process:
firstly, solving an elastic wave equation:
Figure BDA0002821151500000041
in the formula (2), ρ (x) is the density of the medium at the spatial position x, and is set as a constant; λ (x) and μ (x) are values of the Lame parameter at spatial position x; f. ofx(t) and fz(t) are a transverse seismic source function and a longitudinal seismic source function, respectively; delta (x-x)s) As a function of the pulse, when x ═ xsWhen the function value is 1, otherwise, the function value is 0, which indicates that the seismic source function is positioned at xsAt least one of (1) and (b); lame parameters lambda (x) and mu (x) and longitudinal wave velocity vp(x) And transverse wave velocity vs(x) The relationship of (a) to (b) is as follows:
Figure BDA0002821151500000042
u in formula (2)x(t,x；xs) And uz(t,x；xs) Respectively indicate when the shot point is located at xsIn time of processing, the transverse component and longitudinal component of wave field at arbitrary spatial position x are made x ═ xrThe transverse component of the calculated wavefield may be obtained
Figure BDA0002821151500000043
And a longitudinal component
Figure BDA0002821151500000044
Namely:
Figure BDA0002821151500000045
in step 3, the envelope is obtained by hilbert transform, defined as follows:
Figure BDA0002821151500000046
Figure BDA0002821151500000051
Figure BDA0002821151500000052
Figure BDA0002821151500000053
in the formulas (7) and (8), H { } represents the hilbert transform, where the parameter τ is a temporary variable introduced during the hilbert transform and represents a time offset.
The specific implementation of step 4 is as follows:
direct envelope inversion aims at the strong scattering problem, and calculates a new Frechet derivative directly based on energy scattering:
Figure BDA0002821151500000054
in the formula (9), the reaction mixture is,
Figure BDA0002821151500000055
and
Figure BDA0002821151500000056
respectively represent the objective function Je(vp(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function Je(vp(x),vs(x) Gradient with respect to shear wave velocity; f. ofeRepresenting the envelope incidental source, equal to the envelope residual;
Figure BDA0002821151500000057
and
Figure BDA0002821151500000058
the envelope frechet derivatives with respect to the compressional and shear wave velocities, respectively, have the following form:
Figure BDA0002821151500000059
formula (10)
Figure BDA0002821151500000061
And
Figure BDA0002821151500000062
virtual source function of envelope, G, with respect to the velocity of longitudinal and transverse waves, respectivelyeGreen's function which is envelope propagation; after the velocity gradients of the longitudinal wave and the transverse wave are calculated by using the equations (9) and (10), the velocity is iteratively updated by using the following criteria:
Figure BDA0002821151500000063
in the formula (11), the reaction mixture is,
Figure BDA0002821151500000064
and
Figure BDA0002821151500000065
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number;
Figure BDA0002821151500000066
and
Figure BDA0002821151500000067
calculating a longitudinal wave velocity gradient and a transverse wave velocity gradient according to an equation (9) and an equation (10) in the nth step iteration process;
Figure BDA0002821151500000068
and
Figure BDA0002821151500000069
the longitudinal wave velocity and the transverse wave velocity of the step (n + 1) after updating;
calculating the velocity gradient of longitudinal wave and the velocity gradient of transverse wave by the method, iterating, and obtaining the long wavelength component v of the velocity model of longitudinal wave after iteration convergencepl(x) And the long wavelength component v of the transverse wave velocity modelsl(x)。
In step 5, v is obtainedpl(x) And vsl(x) Then, the next stage objective function J (v)p(x),vs(x) Defined as the two-norm of the observed wavefield and the calculated wavefield residual, of the form:
Figure BDA00028211515000000610
in step 6, the objective function J (v) is processedp(x),vs(x) The velocity model gradient equation at this stage of the optimization is as follows:
Figure BDA00028211515000000611
in the formula (13), ζvp(x) And ζvs(x) Respectively represent the objective function J (v)p(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function J (v)p(x),vs(x) Gradient with respect to shear wave velocity; f represents the wavefield satellite, equal to the wavefield residual; fpAnd FpThe derivatives of the wavefield frechet with respect to compressional and shear wave velocities, respectively, have the following form:
Figure BDA0002821151500000071
formula (14) QpAnd QsThe wave field virtual source functions related to longitudinal wave velocity and transverse wave velocity respectively, and G is a Green function of wave field propagation; after the compressional velocity gradient and the shear velocity gradient are calculated by the equations (13) and (14), the compressional velocity and the shear velocity are iteratively updated by the same method as the equation (11), that is:
Figure BDA0002821151500000072
in the formula (15), the reaction mixture is,
Figure BDA0002821151500000073
and
Figure BDA0002821151500000074
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number; (ζ)vp(x))nAnd (ζ)vs(x))nCalculating the longitudinal wave velocity gradient and the transverse wave velocity gradient according to the formula (13) and the formula (14) in the nth step iteration process;
Figure BDA0002821151500000075
and
Figure BDA0002821151500000076
the longitudinal wave velocity and the transverse wave velocity of the step (n + 1) after updating;
through the method, the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated and iterated, and the fine structure v of the longitudinal wave velocity model can be obtained after iteration convergencepf(x) And fine structure v of shear wave velocity modelsf(x)。
The invention has the beneficial effects that: firstly, calculating envelopes of an observation wave field and a calculation wave field, and establishing an envelope objective function; secondly, starting from a linear initial model, obtaining a long-wavelength component of a longitudinal wave velocity and a long-wavelength component of a transverse wave velocity of the salt dome by using a direct envelope inversion method; thirdly, establishing a wave field target function according to the calculated wave field and the observed wave field; and finally, based on the long wavelength velocity component, obtaining a fine velocity structure by utilizing the traditional elastic wave full waveform inversion.
Compared with a common inversion method, the method is more suitable for inversion of the parameters of the strong scattering medium, and can simultaneously and accurately estimate the longitudinal wave velocity and the transverse wave velocity through an elastic wave equation. The reason why the method has the above-described advantages is that: firstly, because a direct envelope inversion method is used, a new Frechet derivative is adopted for energy propagation, and weak scattering approximation of the traditional inversion is avoided; secondly, because the elastic wave equation is used, multi-parameter simultaneous inversion can be realized.
Drawings
FIG. 1 is a flow chart of the method of the present invention;
FIG. 2 is a model of true longitudinal wave velocity in the target zone;
FIG. 3 is a model of true shear wave velocity in the target zone;
FIG. 4 is a linear initial compressional velocity model;
FIG. 5 is a linear initial shear velocity model;
FIG. 6 is a longitudinal velocity model obtained by inversion using a conventional method;
FIG. 7 is a shear velocity model obtained by inversion using conventional methods;
FIG. 8 shows the long wavelength component of longitudinal wave velocity inverted by the method of the present invention;
FIG. 9 shows the long wavelength component of transverse wave velocity obtained by inversion using the method proposed in the present invention;
FIG. 10 is a fine compressional velocity model obtained by inversion using the method proposed by the present invention;
FIG. 11 is a fine shear velocity model obtained by inversion using the method proposed by the present invention.
Detailed Description
The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.
The invention adopts an elastic wave full waveform inversion method of a direct envelope sensitive kernel function, as shown in figure 1, the specific steps are respectively as follows:
step 1, setting a shot point and a demodulator probe, and acquiring original seismic data to obtain a transverse component of observation data
Figure BDA0002821151500000091
And a longitudinal component
Figure BDA0002821151500000092
Wherein t represents a time variable, xr,xsThe positions of the demodulator probe and the shot point are respectively, and the collected seismic data are subjected to conventional preprocessing, including dynamic and static correction, noise elimination and the like.
And 2, constructing a rectangular grid geological model, setting the number Nx of transverse grid points and the number Nz of longitudinal grid points of the model according to technical indexes and requirements of a work area, and setting a forward simulation space sampling interval (comprising a transverse sampling interval dx and a longitudinal sampling interval dz), a time sampling interval dt and maximum sampling time Nt. In practice, relevant parameters can be determined according to the effective frequency band range of seismic data, the seismic recording time length and a seismic data observation system; the criterion for selecting the grid parameters is to make it possible to not only satisfy the stability condition but also effectively suppress the numerical dispersion when performing finite difference forward modeling based on the grid.
Step 3, setting a longitudinal wave velocity model vp(x) And transverse wave velocity model vs(x) And initial longitudinal wave velocity modeForm vp0(x) And an initial shear velocity model vs0(x) Given an objective function J to be optimized in the direct envelope inversion stagee(vp(x),vs(x) ); longitudinal wave velocity model vp(x) And an initial model v of the shear velocity modelp0(x) And vs0(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 transverse uniform longitudinal linear increase; the objective function defines the matching degree between the observation data and the calculation data, and the objective function is optimized to update the model parameters. In the invention, a target function J is set at the direct envelope inversion stagee(vp(x),vs(x) Is a two-norm of the residual between the observed data envelope and the calculated data envelope, having the form:
Figure BDA0002821151500000093
in the formula (1), the reaction mixture is,
Figure BDA0002821151500000094
and
Figure BDA0002821151500000095
respectively representing the transverse components of the observed wavefield
Figure BDA0002821151500000096
And calculating the wavefield transverse component
Figure BDA0002821151500000097
The envelope of the (c) signal is,
Figure BDA0002821151500000098
and
Figure BDA0002821151500000099
respectively representing the longitudinal components of the observed wavefield
Figure BDA00028211515000000910
And meterComputing longitudinal components of wavefields
Figure BDA00028211515000000911
Wherein the transverse component and the longitudinal component of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time;
Figure BDA00028211515000000912
indicating that the summation operation is performed with respect to the shot point and the demodulator probe.
In step 3, the calculation of the transverse component of the wavefield
Figure BDA0002821151500000101
And a longitudinal component
Figure BDA0002821151500000102
This can be obtained by first solving the elastic wave equation as follows:
Figure BDA0002821151500000103
in the formula (2), ρ (x) is the density of the medium at the spatial position x, and is set as a constant; λ (x) and μ (x) are values of the Lame parameter at spatial position x; f. ofx(t) and fz(t) are a transverse seismic source function and a longitudinal seismic source function, respectively; delta (x-x)s) As a function of the pulse, when x ═ xsWhen the function value is 1, otherwise, the function value is 0, which indicates that the seismic source function is positioned at xsAt least one of (1) and (b); lame parameters lambda (x) and mu (x) and longitudinal wave velocity vp(x) And transverse wave velocity vs(x) The relationship of (a) to (b) is as follows:
Figure BDA0002821151500000104
u in formula (2)x(t,x；xs) And uz(t,x；xs) Respectively indicate when the shot point is located at xsIn time of processing, the transverse component and longitudinal component of wave field at arbitrary spatial position x are made x ═ xrThe transverse component of the calculated wavefield may be obtained
Figure BDA0002821151500000105
And a longitudinal component
Figure BDA0002821151500000106
Namely:
Figure BDA0002821151500000107
in step (3), the envelope may be obtained by hilbert transform, and is defined as follows:
Figure BDA0002821151500000111
Figure BDA0002821151500000112
Figure BDA0002821151500000113
Figure BDA0002821151500000114
in the formulas (7) and (8), H { } represents the hilbert transform, where the parameter τ is a temporary variable introduced during the hilbert transform and represents a time offset.
The above objective function Je(vp(x),vs(x) Envelope information of the signal) is used because the envelope of the signal contains abundant low-frequency components, which helps to solve the model long-wavelength components.
Step 4, carrying out direct envelope inversion and iterative optimization algorithm on the target function Je(vp(x),vs(x) Optimized to obtain the long wavelength component v of the longitudinal wave velocity modelpl(x) And the long wavelength component v of the transverse wave velocity modelsl(x) (ii) a Tradition ofThe inversion method is based on weak scattering approximation and uses a chain rule to calculate the derivative of an objective function about model parameters, and the direct envelope inversion in the invention is directed to the problem of strong scattering and directly calculates a new Frechet derivative based on energy scattering:
Figure BDA0002821151500000121
in the formula (9), the reaction mixture is,
Figure BDA0002821151500000122
and
Figure BDA0002821151500000123
respectively represent the objective function Je(vp(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function Je(vp(x),vs(x) Gradient with respect to shear wave velocity; f. ofeRepresenting the envelope incidental source, equal to the envelope residual;
Figure BDA0002821151500000124
and
Figure BDA0002821151500000125
the envelope frechet derivatives with respect to the compressional and shear wave velocities, respectively, have the following form:
Figure BDA0002821151500000126
formula (10)
Figure BDA0002821151500000127
And
Figure BDA0002821151500000128
respectively, the imaginary source function of the envelope, G, with respect to the velocity of the longitudinal and transverse waveseGreen's function of envelope propagation. After the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated by the equations (9) and (10), the velocity is iterated by the following criteriaUpdating:
Figure BDA0002821151500000129
in the formula (11), the reaction mixture is,
Figure BDA00028211515000001210
and
Figure BDA00028211515000001211
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number;
Figure BDA00028211515000001212
and
Figure BDA00028211515000001213
calculating a longitudinal wave velocity gradient and a transverse wave velocity gradient according to an equation (9) and an equation (10) in the nth step iteration process;
Figure BDA00028211515000001214
and
Figure BDA00028211515000001215
the longitudinal wave velocity and the transverse wave velocity of the (n + 1) th step after updating.
The longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated and iterated through the method, and the long wavelength component v of the longitudinal wave velocity model can be obtained after iterative convergencepl(x) And the long wavelength component v of the transverse wave velocity modelsl(x)。
Step 5, a new objective function J (v) to be optimized in the next stage is givenp(x),vs(x) ); to obtain vpl(x) And vsl(x) Then, the next stage objective function J (v)p(x),vs(x) Defined as the two-norm of the observed wavefield and the calculated wavefield residual, of the form:
Figure BDA0002821151500000131
step 6, aiming at the target function J (v)p(x),vs(x) ) is optimized to obtain a fine structure v of the longitudinal wave velocity modelpf(x) And fine structure v of shear wave velocity modelsf(x) (ii) a And (3) a fine structure inversion stage, namely performing iterative optimization on model parameters according to an objective function defined by a formula (12). For the objective function J (v)p(x),vs(x) The velocity model gradient equation at this stage of the optimization is as follows:
Figure BDA0002821151500000132
in the formula (13), ζvp(x) And ζvs(x) Respectively represent the objective function J (v)p(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function J (v)p(x),vs(x) Gradient with respect to shear wave velocity; f represents the wavefield satellite, equal to the wavefield residual; fpAnd FpThe derivatives of the wavefield frechet with respect to compressional and shear wave velocities, respectively, have the following form:
Figure BDA0002821151500000133
formula (14) QpAnd QsThe wavefield virtual source functions for compressional and shear wave velocities, respectively, and G is the green's function of wavefield propagation. After the compressional velocity gradient and the shear velocity gradient are calculated by the equations (13) and (14), the compressional velocity and the shear velocity are iteratively updated by the same method as the equation (11), that is:
Figure BDA0002821151500000141
in the formula (15), the reaction mixture is,
Figure BDA0002821151500000142
and
Figure BDA0002821151500000143
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number; (ζ)vp(x))nAnd (ζ)vs(x))nCalculating the longitudinal wave velocity gradient and the transverse wave velocity gradient according to the formula (13) and the formula (14) in the nth step iteration process;
Figure BDA0002821151500000144
and
Figure BDA0002821151500000145
the longitudinal wave velocity and the transverse wave velocity of the step (n + 1) after updating;
through the method, the longitudinal wave velocity gradient and the transverse wave velocity gradient are calculated and iterated, and the fine structure v of the longitudinal wave velocity model can be obtained after iteration convergencepf(x) And fine structure v of shear wave velocity modelsf(x)。
Model example
The specific implementation process of the invention is applied to a certain strong scattering salt dome geological model. The model has a transverse direction of 16100 m and a longitudinal direction of 3750 m. The model contains a large range of high velocity salt hills, while the background medium exhibits a low velocity distribution, which is typical of strong scattering properties. The compressional velocity model and shear velocity model of the target zone are shown in figures 2 and 3, respectively.
The initial compressional velocity model and the initial shear velocity model were set as linear models as shown in FIGS. 4 and 5, respectively, in which the compressional velocity model gradually linearly increased from 1530m/s to 3930m/s from shallow to deep, and the shear velocity model gradually linearly increased from 950m/s to 2000kg/m from shallow to deep3。
The longitudinal wave velocity model and the transverse wave velocity model obtained by the conventional elastic wave full waveform inversion method are respectively shown in fig. 6 and 7. It can be seen that the traditional 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 show the long wavelength component of the longitudinal wave velocity model and the long wavelength component of the transverse wave velocity model of the target region obtained by performing the inversion through steps 1 to 4 of the method proposed by the present invention. The accurate acquisition of the long wavelength component, i.e. the large-scale background model, is the key to the whole parametric inversion in the inversion process, since it ensures that the iteration can converge to the correct globally optimal solution and the subsequent further inversion of the fine structure. As can be seen from the results shown in fig. 8 and 9, the method proposed in the present invention gives a good profile of the salt dome as well as the velocity profile.
Based on the long wavelength velocity model obtained as described above, the fine structures of the longitudinal and transverse velocity models 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 result is close to the real model, which illustrates 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:
step 1, setting a shot point and a wave detection point, and acquiring transverse components of observation data
Figure FDA0002821151490000011
And a longitudinal component
Figure FDA0002821151490000012
Wherein t represents a time variable; x is the number ofr,xsRespectively representing the positions of the demodulator probe and the shot point;
step 2, constructing a rectangular grid geological model, setting the number Nx of transverse grid points and the number Nz of longitudinal grid points of the model, and setting a space sampling interval, a time sampling interval dt and maximum sampling time Nt of forward simulation, wherein the space sampling interval comprises a transverse sampling interval dx and a longitudinal sampling interval dz;
step 3, setting a longitudinal wave velocity model vp(x) And transverse wave velocity model vs(x) And an initial longitudinal wave velocity model vp0(x) And an initial shear velocity model vs0(x) Given an objective function J to be optimized in the direct envelope inversion stagee(vp(x),vs(x))；
Step 4, carrying out direct envelope inversion and iterative optimization algorithm on the target function Je(vp(x),vs(x) Optimized to obtain the long wavelength component v of the longitudinal wave velocity modelpl(x) And the long wavelength component v of the transverse wave velocity modelsl(x)；
Step 5, a new objective function J (v) to be optimized in the next stage is givenp(x),vs(x))；
Step 6, aiming at the target function J (v)p(x),vs(x) ) is optimized to obtain a fine structure v of the longitudinal wave velocity modelpf(x) And fine structure v of shear wave velocity modelsf(x)。
2. The method of claim 1, wherein in step 3, a target function J is sete(vp(x),vs(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 reaction mixture is,
Figure FDA0002821151490000021
and
Figure FDA0002821151490000022
respectively representing the transverse components of the observed wavefield
Figure FDA0002821151490000023
And calculating the wavefield transverse component
Figure FDA0002821151490000024
The envelope of the (c) signal is,
Figure FDA0002821151490000025
and
Figure FDA0002821151490000026
respectively representing the longitudinal components of the observed wavefield
Figure FDA0002821151490000027
And calculating the wavefield longitudinal component
Figure FDA0002821151490000028
Wherein the transverse component and the longitudinal component of the computed wavefield are computed by finite difference forward modeling; t is the maximum observation time;
Figure FDA0002821151490000029
indicating that the summation operation is performed with respect to the shot point and the demodulator probe.
3. The method of claim 2, wherein in step 3, the transverse component of the wavefield is calculated
Figure FDA00028211514900000210
And a longitudinal component
Figure FDA00028211514900000211
Can be obtained by the following process:
firstly, solving an elastic wave equation:
Figure FDA00028211514900000212
in the formula (2), ρ (x) is the density of the medium at the spatial position x, and is set as a constant; λ (x) and μ (x) are values of the Lame parameter at spatial position x; f. ofx(t) and fz(t) is a transverse seismic source function and a longitudinal seismic source function respectivelyCounting; delta (x-x)s) As a function of the pulse, when x ═ xsWhen the function value is 1, otherwise, the function value is 0, which indicates that the seismic source function is positioned at xsAt least one of (1) and (b); lame parameters lambda (x) and mu (x) and longitudinal wave velocity vp(x) And transverse wave velocity vs(x) The relationship of (a) to (b) is as follows:
Figure FDA00028211514900000213
u in formula (2)x(t,x；xs) And uz(t,x；xs) Respectively indicate when the shot point is located at xsIn time of processing, the transverse component and longitudinal component of wave field at arbitrary spatial position x are made x ═ xrThe transverse component of the calculated wavefield may be obtained
Figure FDA0002821151490000031
And a longitudinal component
Figure FDA0002821151490000032
Namely:
Figure FDA0002821151490000033
4. the method of claim 3, wherein in step 3, the envelope is obtained by a Hilbert transform and is defined as follows:
Figure FDA0002821151490000034
Figure FDA0002821151490000035
Figure FDA0002821151490000036
Figure FDA0002821151490000037
in the formulas (7) and (8), H { } represents the hilbert transform, where the parameter τ is a temporary variable introduced during the hilbert transform and represents a time offset.
5. The method for inverting the full waveform of an elastic wave using a direct envelope sensitivity kernel as claimed in claim 4, wherein the specific implementation of step 4 is as follows:
direct envelope inversion aims at strong scattering problem, and new energy scattering direct-based calculation method
Figure FDA0002821151490000038
Derivative:
Figure FDA0002821151490000041
in the formula (9), the reaction mixture is,
Figure FDA0002821151490000042
and
Figure FDA0002821151490000043
respectively represent the objective function Je(vp(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function Je(vp(x),vs(x) Gradient with respect to shear wave velocity; f. ofeRepresenting the envelope incidental source, equal to the envelope residual;
Figure FDA0002821151490000044
and
Figure FDA0002821151490000045
envelopes relating to longitudinal and transverse wave velocities, respectively
Figure FDA0002821151490000046
Derivative, having the form:
Figure FDA0002821151490000047
formula (10)
Figure FDA0002821151490000048
And
Figure FDA0002821151490000049
respectively, the imaginary source function of the envelope, G, with respect to the velocity of the longitudinal and transverse waveseGreen's function which is envelope propagation; after the velocity gradients of the longitudinal wave and the transverse wave are calculated by using the equations (9) and (10), the velocity is iteratively updated by using the following criteria:
Figure FDA00028211514900000410
in the formula (11), the reaction mixture is,
Figure FDA00028211514900000411
and
Figure FDA00028211514900000412
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number;
Figure FDA00028211514900000413
and
Figure FDA00028211514900000414
according to the formula (9) andcalculating a longitudinal wave velocity gradient and a transverse wave velocity gradient by the formula (10);
Figure FDA00028211514900000415
and
Figure FDA00028211514900000416
the longitudinal wave velocity and the transverse wave velocity of the step (n + 1) after updating;
calculating the velocity gradient of longitudinal wave and the velocity gradient of transverse wave by the method, iterating, and obtaining the long wavelength component v of the velocity model of longitudinal wave after iteration convergencepl(x) And the long wavelength component v of the transverse wave velocity modelsl(x)。
6. The method of claim 5, wherein in step 5, v is obtainedpl(x) And vsl(x) Then, the next stage objective function J (v)p(x),vs(x) Defined as the two-norm of the observed wavefield and the calculated wavefield residual, of the form:
Figure FDA0002821151490000051
7. the method of claim 6, wherein in step 6, the target function J (v) is selectedp(x),vs(x) The velocity model gradient equation at this stage of the optimization is as follows:
Figure FDA0002821151490000052
in the formula (13), ζvp(x) And ζvs(x) Respectively represent the objective function J (v)p(x),vs(x) Gradient with respect to longitudinal wave velocity and objective function J (v)p(x),vs(x) About a cross barA gradient of wave velocity; f represents the wavefield satellite, equal to the wavefield residual; fpAnd FpThe derivatives of the wavefield frechet with respect to compressional and shear wave velocities, respectively, have the following form:
Figure FDA0002821151490000053
formula (14) QpAnd QsThe wave field virtual source functions related to longitudinal wave velocity and transverse wave velocity respectively, and G is a Green function of wave field propagation; after the compressional velocity gradient and the shear velocity gradient are calculated by the equations (13) and (14), the compressional velocity and the shear velocity are iteratively updated by the same method as the equation (11), that is:
Figure FDA0002821151490000061
in the formula (15), the reaction mixture is,
Figure FDA0002821151490000062
and
Figure FDA0002821151490000063
respectively the longitudinal wave speed and the transverse wave speed of the nth step iteration; alpha is alphanThe iteration step length of the nth step is a normal number; (ζ)vp(x))nAnd (ζ)vs(x))nCalculating the longitudinal wave velocity gradient and the transverse wave velocity gradient according to the formula (13) and the formula (14) in the nth step iteration process;
Figure FDA0002821151490000064
and
Figure FDA0002821151490000065
the longitudinal wave velocity and the transverse wave velocity of the step (n + 1) after updating;
calculating the velocity gradient of longitudinal wave and the velocity gradient of transverse wave by the method, iterating, and repeatingAfter convergence, a fine structure v of a longitudinal wave velocity model can be obtainedpf(x) And fine structure v of shear wave velocity modelsf(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 true CN112684505A (en) 2021-04-20
CN112684505B 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 (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180017690A1 (en) * 2016-07-13 2018-01-18 Sirui Tan Joint Full Wavefield Inversion of P-Wave Velocity and Attenuation Using an Efficient First Order Optimization
CN110007340A (en) * 2019-02-01 2019-07-12 西安理工大学 Salt dome speed density estimation method based on the direct envelope inverting of angle domain

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180017690A1 (en) * 2016-07-13 2018-01-18 Sirui Tan Joint Full Wavefield Inversion of P-Wave Velocity and Attenuation Using an Efficient First Order Optimization
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;: "二维时间空间域弹性波全波形反演-理论模型数值试验", 地球物理学进展, no. 02 *

Also Published As

Publication number Publication date
CN112684505B (en) 2023-06-30

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
CN106526674B (en) Three-dimensional full waveform inversion energy weighting gradient preprocessing method
CN109188519B (en) System and method for inverting longitudinal and transverse wave speeds of elastic waves under polar coordinates
CN110187382B (en) Traveling time inversion method for wave equation of reverse wave and reflected wave
WO2011139413A1 (en) Artifact reduction in method of iterative inversion of geophysical data
US20130258810A1 (en) Method and System for Tomographic Inversion
CN113740901B (en) Land seismic data full-waveform inversion method and device based on complex undulating surface
CN111025387B (en) Pre-stack earthquake multi-parameter inversion method for shale reservoir
CN113031068B (en) Reflection coefficient accurate base tracking prestack seismic inversion method
CN111045076A (en) Multi-mode Rayleigh wave frequency dispersion curve parallel joint inversion method
CN111025388B (en) Multi-wave combined prestack waveform inversion method
Ba et al. The dynamic stiffness matrix method for seismograms synthesis for layered transversely isotropic half-space
CN110888159A (en) Elastic wave full waveform inversion method based on angle decomposition and wave field separation
CN111665556B (en) Stratum acoustic wave propagation velocity model construction method
NO20190489A1 (en) Seismic modeling
CN111999764B (en) Method for constructing least square reverse time migration under salt based on time-frequency domain objective function
CN111505714B (en) Elastic wave direct envelope inversion method based on rock physical constraint
CN114624766B (en) Elastic wave least square reverse time migration gradient solving method based on traveling wave separation
CN108680957B (en) Local cross-correlation time-frequency domain Phase-retrieval method based on weighting
CN112684505A (en) Elastic wave full waveform inversion method adopting direct envelope sensitive kernel function
CN110161561A (en) A kind of controllable layer position sublevel interbed multiple analogy method in oil and gas reservoir
CN111175822B (en) Strong scattering medium inversion method for improving direct envelope inversion and disturbance decomposition
CN112305595B (en) Method for analyzing geologic body structure based on refraction wave and storage medium

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