CN113126151A - Elastic reflection wave travel time inversion method based on pure wave continuation equation - Google Patents

Elastic reflection wave travel time inversion method based on pure wave continuation equation Download PDF

Info

Publication number
CN113126151A
CN113126151A CN202110260212.3A CN202110260212A CN113126151A CN 113126151 A CN113126151 A CN 113126151A CN 202110260212 A CN202110260212 A CN 202110260212A CN 113126151 A CN113126151 A CN 113126151A
Authority
CN
China
Prior art keywords
wave
pure
background
equation
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
CN202110260212.3A
Other languages
Chinese (zh)
Other versions
CN113126151B (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.)
Oceanographic Instrumentation Research Institute Shandong Academy of Sciences
Original Assignee
Oceanographic Instrumentation Research Institute Shandong Academy of Sciences
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 Oceanographic Instrumentation Research Institute Shandong Academy of Sciences filed Critical Oceanographic Instrumentation Research Institute Shandong Academy of Sciences
Priority to CN202110260212.3A priority Critical patent/CN113126151B/en
Publication of CN113126151A publication Critical patent/CN113126151A/en
Application granted granted Critical
Publication of CN113126151B publication Critical patent/CN113126151B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/675—Wave equation; Green's functions

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses an elastic reflection wave travel time inversion method based on a pure wave continuation equation, which comprises the following steps: (1) omitting a cross term from the approximate elastic wave decoupling continuation equation to obtain a pure wave continuation equation, and calculating a background wave field based on the pure wave continuation equation; (2) calculating a reverse time migration result by using the obtained background wave field; (3) solving a disturbance pure wave continuation equation based on the pure wave continuation equation and the reverse time migration result, and calculating a disturbance wave field by using the disturbance pure wave continuation equation; (4) respectively calculating longitudinal and transverse wave gradients based on the background wave field obtained in the step (1) and the disturbance wave field obtained in the step (3); (5) and calculating the step size according to the longitudinal and transverse wave gradients, updating the longitudinal and transverse wave speeds, and outputting the updated longitudinal and transverse wave speeds as a final inversion result when the iteration value converges to the global minimum value. The method disclosed by the invention can improve the inversion calculation efficiency and the inversion precision and provide technical support for elastic wave seismic exploration.

Description

Elastic reflection wave travel time inversion method based on pure wave continuation equation
Technical Field
The invention belongs to the field of elastic wave seismic exploration, and particularly relates to an elastic reflection wave travel time inversion method based on a pure wave continuation equation.
Background
The elastic reflection wave travel time inversion is based on an elastic wave theory, the simulated seismic reflection data and the observed seismic reflection data are matched, the velocity component of the low and medium wave number is recovered along the propagation path of the elastic wave, an accurate initial velocity model can be provided for elastic full waveform inversion, elastic reverse time migration and the like, and the method is one of the leading edge technologies in the field of elastic wave seismic exploration at present. However, an offset/anti-offset process needs to be introduced to calculate the gradient sensitive kernel function of the inversion during the elastic reflected wave travel, and the elastic wave fluctuation equation is solved for 6 times in one iteration process, so that the step of practical production and application is restricted by huge calculation amount and massive storage requirements.
In addition, due to the simultaneous existence of multi-wave modes in the elastic wave field, serious high wave number crosstalk exists in the inversion of elastic reflection waves during traveling, which is mainly reflected in a sensitive kernel function of transverse wave gradient. Elastic wave decoupling continuation can effectively weaken the coupling effect of elastic parameter inversion, but transverse wave stress wave field decoupling of reverse continuation still can introduce high wave number noise, and the precision of transverse wave velocity inversion is reduced.
Disclosure of Invention
In order to solve the technical problems, the invention provides an elastic reflection wave travel time inversion method based on a pure wave continuation equation so as to achieve the purpose of improving inversion calculation efficiency and inversion accuracy.
In order to achieve the purpose, the technical scheme of the invention is as follows:
an elastic reflection wave travel time inversion method based on a pure wave continuation equation comprises the following steps:
(1) omitting a cross term from the approximate elastic wave decoupling continuation equation to obtain a pure wave continuation equation, and calculating a background wave field based on the pure wave continuation equation;
(2) calculating a reverse time migration result by using the obtained background wave field;
(3) solving a disturbance pure wave continuation equation based on the pure wave continuation equation and the reverse time migration result, and calculating a disturbance wave field by using the disturbance pure wave continuation equation;
(4) respectively calculating longitudinal and transverse wave gradients based on the background wave field obtained in the step (1) and the disturbance wave field obtained in the step (3);
(5) and calculating the step size according to the longitudinal and transverse wave gradients, updating the longitudinal and transverse wave speeds, and outputting the updated longitudinal and transverse wave speeds as a final inversion result when the iteration value converges to the global minimum value.
In the above scheme, the step (1) is specifically as follows: the pure wave continuation equation comprises a pure P wave continuation equation and a pure S wave continuation equation, a forward-propagating P wave vibration velocity background wave field and a backward-propagating P wave vibration velocity background wave field are calculated based on the pure P wave continuation equation respectively, a backward-propagating S wave vibration velocity background wave field is calculated based on the pure S wave continuation equation, a backward-propagating P wave stress background wave field is calculated based on the pure P wave continuation equation, a backward-propagating S wave stress background wave field is calculated based on the pure S wave continuation equation, and numerical simulation calculation of the background wave field is implemented by high-order staggered grid finite difference.
In a further technical scheme, a pure P wave continuation equation is as follows:
Figure BDA0002969635470000021
wherein the subscript "0" represents the background, x and z represent the spatial horizontal and vertical directions, respectively, t represents time, vPx0And vPz0Representing the horizontal and vertical components, tau, of the P-wave vibration velocity background wave field, respectivelyPxx0And τPzz0Representing the P-wave stress background wavefield, P, in xx and zz directions, respectively0And VP0Respectively representing background density and background longitudinal wave velocity;
the pure S wave continuation equation is:
Figure BDA0002969635470000022
wherein v isSx0And vSz0Representing the horizontal and vertical components, tau, of the S-wave vibration velocity background wave field, respectivelySxz0And τSzx0S-wave stress background wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity.
In the above scheme, the step (2) is specifically as follows: and respectively calculating a PP wave reverse time migration result by cross-correlation of the forward-transmitted P wave vibration velocity background wave field and the backward-transmitted P wave vibration velocity background wave field, and calculating a PS wave reverse time migration result by using the forward-transmitted P wave vibration velocity background wave field and the backward-transmitted S wave vibration velocity background wave field.
In the above scheme, the step (3) is specifically as follows: establishing a Born forward equation, namely a disturbed pure wave continuation equation, by taking the interaction of the background wave field and the velocity disturbance as a seismic source, wherein the disturbed pure wave continuation equation comprises a disturbed pure P wave continuation equation and a disturbed pure S wave continuation equation; wherein, the speed disturbance is replaced by a reverse time migration result;
and calculating a forward-transmitted P wave vibration velocity disturbance wave field and a backward-transmitted P wave stress disturbance wave field based on a disturbed pure P wave continuation equation, and calculating a forward-transmitted S wave vibration velocity disturbance wave field based on a disturbed pure S wave continuation equation.
In the above scheme, the pure P wave continuation equation of the disturbance is:
Figure BDA0002969635470000031
where the symbol "δ" represents a disturbance, ρ0Represents the background density, δ VPRepresenting longitudinal wave velocity disturbances, δ VSRepresenting shear wave velocity perturbations, vPx0And vPz0Respectively representing the horizontal and vertical components, V, of the P-wave vibration velocity background wave fieldP0Representing background longitudinal wave velocity, δ vPxAnd δ vPzRespectively representing the horizontal and vertical components, delta tau, of the P-wave vibration velocity disturbance wavefieldPxxAnd δ τPzzRespectively representing P wave stress disturbance wave fields in xx and zz directions; in practice, δ V in equation (3)PAnd δ VSExpressed by PP wave reverse time migration result;
the pure S wave continuation equation of the disturbance is as follows:
Figure BDA0002969635470000032
in practice, δ V in equation (4)SExpressed as a result of reverse time shift of the PS wave, δ vSxAnd δ vSzRespectively representing the horizontal and vertical components, delta tau, of the S-wave vibration velocity disturbance wavefieldSxzAnd δ τSzxS-wave stress disturbance wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity.
In the scheme, the specific method in the step (4) is as follows: reading the background wave field calculated in the step (1) and the disturbance wave field calculated in the step (3), and respectively solving the velocity gradient of the longitudinal wave
Figure BDA0002969635470000033
And transverse wave velocity gradient
Figure BDA0002969635470000034
Maintaining background density ρ throughout the calculation0Is a constant;
Figure BDA0002969635470000041
wherein, χPHexix-SRespectively representing the objective functions, V, of the matched P-wave data and the matched S-wave dataP0Representing background longitudinal wave velocity, vPx0And vPz0Respectively representing the horizontal and vertical components, deltav, of the forward propagating P-wave vibration velocity background wavefieldPxAnd δ vPzRespectively representing the horizontal component and the vertical component of the forward propagating P-wave vibration velocity disturbance wave field,
Figure BDA0002969635470000042
and
Figure BDA0002969635470000043
representing the xx and zz components of the backward propagating P-wave stress perturbation wavefield,
Figure BDA0002969635470000044
and
Figure BDA0002969635470000045
representing the xx and zz components, V, respectively, of the backward propagating P-wave stress background wavefieldS0Representing the background shear wave velocity, δ vSxAnd δ vSzRespectively representing the horizontal component and the vertical component of the forward propagating S-wave vibration velocity perturbation wave field,
Figure BDA0002969635470000046
and
Figure BDA0002969635470000047
representing the xz and zx components, respectively, of the back-propagating S-wave stress background wavefield.
In the scheme, in the step (5), the update amount of the longitudinal wave velocity of each iteration is ensured to be between 20m/s and 100m/s, and the update amount of the transverse wave velocity is ensured to be between 10m/s and 60 m/s.
Through the technical scheme, the elastic reflection wave travel time inversion method based on the pure wave continuation equation has the following beneficial effects:
according to the invention, the corresponding pure wave continuation equation is obtained by utilizing the approximate elastic wave decoupling continuation equation, and in the implementation process of the elastic reflection wave traveling time inversion, all wave fields are calculated based on the pure wave continuation equation, so that the calculation amount and the storage amount of the elastic reflection wave traveling time inversion can be effectively reduced, and the calculation efficiency is improved by more than 3 times. In addition, longitudinal wave components exist in the shear wave stress wave field calculated by elastic wave decoupling continuation, high wave number noise is brought to shear wave gradient calculation, the longitudinal wave components in the shear wave stress wave field can be effectively eliminated by adopting a pure wave continuation equation, and inversion accuracy is improved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below.
FIG. 1 is a schematic flow chart of an implementation of an elastic reflection wave travel time inversion method based on a pure wave continuation equation;
FIGS. 2a-2d are single layer velocity models, wherein FIG. 2a is a true compressional velocity model; FIG. 2b is a true shear velocity model; FIG. 2c is an initial compressional velocity model; FIG. 2d is an initial shear velocity model
3a-3f are comparisons of the backward shear wave stress wavefield calculated based on the elastic wave decoupling prolongation equation and the pure wave prolongation equation, where FIGS. 3a, 3b, and 3c are xx, zz, and xz components, respectively, of the backward shear wave stress wavefield calculated based on the elastic wave decoupling prolongation equation; 3d, 3e and 3f are xx, zz and xz components, respectively, of an inverse shear wave stress wavefield calculated based on a pure wave prolongation equation;
4a-4b are comparisons of shear wave gradients calculated based on an elastic wave decoupling prolongation equation and a pure wave prolongation equation, where FIG. 4a is the shear wave gradient calculated based on the elastic wave decoupling prolongation equation; FIG. 4b is a shear wave gradient calculated based on the pure wave continuation equation;
FIGS. 5a-5d are Sigsbee2A velocity models, where FIG. 5a is a true longitudinal velocity model; FIG. 5b is a true shear velocity model; FIG. 5c is an initial compressional velocity model; FIG. 5d is an initial shear velocity model;
6a-6b are imaging results of initial velocity, wherein FIG. 6a is PP wave imaging results; FIG. 6b shows the result of PS wave imaging;
FIGS. 7a-7b are final velocity inversion results, where FIG. 7a is compressional velocity; FIG. 7b is the shear wave velocity;
8a-8f are pseudo-well plot comparisons of true velocity, initial velocity, and inversion velocity, where FIGS. 8a, 8b, and 8c are pseudo-well plot comparisons of compressional velocity at lateral distances 1200m, 1500m, and 2200m, respectively; 8d, 8e and 8f are cross-wave velocity pseudo-well curve comparisons at lateral distances 1200m, 1500m and 2200m, respectively;
FIGS. 9a-9b are imaging results of inversion velocities, where FIG. 9a is a PP wave imaging result; fig. 9b shows the PS wave imaging result.
Detailed Description
The technical solution in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.
The invention provides an elastic reflection wave travel time inversion method based on a pure wave continuation equation, as shown in figure 1, a specific technical scheme is explained through model testing:
first, input data
Inputting seismic wavelet data of each shot point; observing seismic data, including P-wave seismic data and S-wave seismic data with longitudinal and transverse wave vector separation; in the longitudinal and transverse wave initial velocity model, the velocity value of the initial velocity generally increases from the shallow layer to the deep layer.
Second, calculating background wave field based on pure wave continuation equation
And (3) utilizing an approximate elastic wave decoupling continuation equation, omitting a cross term and obtaining a pure wave continuation equation comprising a pure P wave continuation equation and a pure S wave continuation equation.
The pure P wave continuation equation is as follows:
Figure BDA0002969635470000051
wherein the subscript "0" represents the background, x and z represent the spatial horizontal and vertical directions, respectively, t represents time, vPx0And vPz0Representing the horizontal and vertical components, tau, of the P-wave vibration velocity background wave field, respectivelyPxx0And τPzz0Representing the P-wave stress background wavefield, P, in xx and zz directions, respectively0And VP0Respectively representing background density and background longitudinal wave velocity;
the pure S wave continuation equation is:
Figure BDA0002969635470000061
wherein，vSx0And vSz0Representing the horizontal and vertical components, tau, of the S-wave vibration velocity background wave field, respectivelySxz0And τSzx0S-wave stress background wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity.
And respectively calculating a forward-propagating P-wave vibration velocity background wave field and a backward-propagating P-wave vibration velocity background wave field based on a pure P-wave continuation equation (1)), and calculating a backward-propagating S-wave vibration velocity background wave field based on a pure S-wave continuation equation (2)), wherein the numerical simulation calculation of the wave fields is implemented by adopting high-order staggered grid finite difference.
A single-layer velocity model (fig. 2a to 2d) was established for numerical simulation testing. Fig. 3a-3f show the comparison of the backward shear wave stress wave field calculated based on the elastic wave decoupling continuation equation and the pure wave continuation equation, and it can be seen that the interference of the compressional wave component exists in the backward shear wave stress wave field calculated based on the elastic wave decoupling continuation equation, and the interference of the compressional wave component can be eliminated by the shear wave stress wave field calculated based on the pure wave continuation equation. It can be seen from fig. 4a that the interference of the longitudinal wave components generates high wave number noise (arrows in fig. 4 a) in the calculation of the shear wave gradient, and it can be seen from fig. 4b that the shear wave gradient calculated based on the pure wave continuation equation can avoid the influence of such high wave number noise.
Thirdly, calculating the reverse time migration result
And utilizing the background wave field obtained in the second step, respectively utilizing the P wave vibration velocity wave field propagated in the forward direction and the P wave vibration velocity wave field propagated in the backward direction to perform cross-correlation calculation on the PP wave reverse time migration result, and utilizing the P wave vibration velocity wave field propagated in the forward direction and the S wave vibration velocity wave field propagated in the backward direction to perform cross-correlation calculation on the PS wave reverse time migration result.
Fig. 5-9 use the sigbee 2A model for numerical testing based on elastic reflection wave travel-time inversion. Fig. 5a to 5d show the initial velocity and the real velocity of the Sigsbee2A model. Fig. 6a and 6b show the reverse time migration result of the initial velocity, where fig. 6a is the PP wave reverse time migration result and fig. 6b is the PS wave reverse time migration result, it can be seen that the reverse time migration results of the PP wave and the PS wave are not effectively returned due to the inaccuracy of the low wave number velocity component in the initial velocity, and particularly the "drawing arc" phenomenon occurs at the position of the mid-deep layer diffraction body.
Fourthly, calculating a disturbance wave field based on a pure wave continuation equation and a reverse time migration result
And constructing a Born forward equation, namely a pure wave continuation equation of the disturbance by taking the interaction of the background wave field and the velocity disturbance as a seismic source, wherein the velocity disturbance is replaced by a reverse time migration result of the third step.
The pure P-wave continuation equation of the disturbance is:
Figure BDA0002969635470000071
where the symbol "δ" represents a disturbance, ρ0Represents the background density, δ VPRepresenting longitudinal wave velocity disturbances, δ VSRepresenting shear wave velocity perturbations, vPx0And vPz0Respectively representing the horizontal and vertical components, V, of the P-wave vibration velocity background wave fieldP0Representing background longitudinal wave velocity, δ vPxAnd δ vPzRespectively representing the horizontal and vertical components, delta tau, of the P-wave vibration velocity disturbance wavefieldPxxAnd δ τPzzRepresenting the P-wave stress perturbation wavefield in xx and zz directions, respectively. In practice, δ V in equation (3)PAnd δ VSThe result of the reverse time shift of the PP wave is shown in FIG. 6 a.
The pure S wave continuation equation of the disturbance is as follows:
Figure BDA0002969635470000072
wherein, δ vSxAnd δ vSzRespectively representing the horizontal and vertical components, delta tau, of the S-wave vibration velocity disturbance wavefieldSxzAnd δ τSzxS-wave stress disturbance wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity. In practice, δ V in equation (4)SThe result of the reverse time shift of the PS wave is shown in FIG. 6 b.
Respectively calculating a backward-propagating P wave stress background wave field based on a pure P wave continuation equation (1)), and calculating a forward-propagating P wave vibration velocity disturbance wave field and a backward-propagating P wave stress disturbance wave field based on a disturbed pure P wave continuation equation (3)); and calculating a backward-propagating S-wave stress background wave field based on a pure S-wave continuation equation (2)), and calculating a forward-propagating S-wave vibration velocity disturbance wave field based on a disturbed pure S-wave continuation equation (4)).
The comparison of wavefield time calculated based on the elastic wave decoupling prolongation equation and the pure wave prolongation equation is compared in tables 1 and 2, and the input speed is shown in fig. 2.
TABLE 1 comparison of calculated time of longitudinal wave gradients
Figure BDA0002969635470000081
Wherein u isP0A P wave background wave field which is forward propagating; delta uPA P wave disturbance wave field which is propagated in the forward direction;
Figure BDA00029696354700000812
is a back-propagating P-wave background wavefield;
Figure BDA00029696354700000813
the wavefield is perturbed for the back-propagating P-waves.
TABLE 2 calculated time comparison of shear wave gradients
Figure BDA0002969635470000082
Wherein u isP0A P wave background wave field which is forward propagating; delta uSAn S-wave disturbance wave field which is propagated in the forward direction;
Figure BDA00029696354700000814
is a back-propagating S-wave background wavefield.
The comparison shows that the time for calculating the wave field based on the pure wave continuation equation is about 1/3 for calculating the wave field time based on the elastic wave decoupling continuation equation, which shows that the invention can effectively improve the calculation efficiency of the elastic reflection wave traveling time inversion.
The fifth step, calculating the gradient of longitudinal and transverse waves
The wave fields calculated in the second step and the fourth step are read, and the velocity gradients of the longitudinal waves are respectively calculated according to equation (5)
Figure BDA0002969635470000083
And transverse wave velocity gradient
Figure BDA0002969635470000084
Maintaining background density ρ throughout the calculation0Is a constant.
Figure BDA0002969635470000085
Wherein, χPHexix-SRespectively representing the objective functions, V, of the matched P-wave data and the matched S-wave dataP0Representing background longitudinal wave velocity, vPx0And vPz0Respectively representing the horizontal and vertical components, deltav, of the forward propagating P-wave vibration velocity background wavefieldPxAnd δ vPzRespectively representing the horizontal component and the vertical component of the forward propagating P-wave vibration velocity disturbance wave field,
Figure BDA0002969635470000086
and
Figure BDA0002969635470000087
representing the xx and zz components of the backward propagating P-wave stress perturbation wavefield,
Figure BDA0002969635470000088
and
Figure BDA0002969635470000089
representing the xx and zz components, V, respectively, of the backward propagating P-wave stress background wavefieldS0Representing the background shear wave velocity, δ vSxAnd δ vSzRespectively representing forward propagation S-wave vibration velocity disturbance wave fieldThe horizontal component and the vertical component of (a),
Figure BDA00029696354700000810
and
Figure BDA00029696354700000811
representing the xz and zx components, respectively, of the back-propagating S-wave stress background wavefield.
Wherein the target function χPHexix-SThe method comprises the following specific steps: recording forward-propagating vibration velocity disturbance wave field data at the surface position as simulated seismic data, respectively matching the simulated P-wave data with the observed P-wave seismic data to establish an objective function equation (6), and matching the simulated S-wave data with the observed S-wave seismic data to establish an objective function equation (8).
Figure BDA0002969635470000091
Wherein, χPRepresenting an objective function established with P-wave data, | · | | non-woven phosphor2Denotes the L2 norm,. DELTA.tauPRepresenting the travel moveout between the simulated P-wave data and the observed P-wave seismic data. Delta tauPThe specific calculation uses equation (7):
Figure BDA0002969635470000092
wherein the content of the first and second substances,
Figure BDA0002969635470000093
and
Figure BDA0002969635470000094
respectively representing simulated P-wave data and observed P-wave seismic data. When the cross-correlation value CPWhen (τ) is maximized, τ and Δ τPAre equal.
Figure BDA0002969635470000095
Wherein, χSRepresenting an objective function, Δ τ, built from S-wave dataSRepresenting the travel moveout between the simulated S-wave data and the observed S-wave seismic data. Delta tauSThe specific calculation uses equation (9):
Figure BDA0002969635470000096
wherein the content of the first and second substances,
Figure BDA0002969635470000097
and
Figure BDA0002969635470000098
respectively representing simulated S-wave data and observed S-wave seismic data. When the cross-correlation value CSWhen (τ) is maximized, τ and Δ τSAre equal.
And finally, calculating the longitudinal and transverse wave gradients of the single cannons through the first step to the fifth step, and then accumulating the gradients of all the cannons to finish the calculation of the longitudinal and transverse wave gradients.
Sixth, updating the longitudinal and transverse wave velocities
And calculating the step size according to the longitudinal wave and transverse wave gradients, and ensuring that the update quantity of the longitudinal wave velocity of each iteration is between 20 and 100m/s and the update quantity of the transverse wave velocity is between 10 and 60 m/s.
Seventh, determining convergence
Calculating the iteration | | delta tauP||2+||ΔτS||2Is compared with the value of the last iteration, and if the value is reduced, the calculation of the second step is continued.
Eighth step, obtaining inversion result of longitudinal and transverse wave velocity
When | | | Δ τP||2+||ΔτS||2If the value of (2) is in a non-descending state for 5 times continuously, the inversion method is considered to be converged to a global minimum value, and at the moment, the output longitudinal and transverse wave velocity data is the final inversion result.
Fig. 7 shows the inversion result of the converged longitudinal and transverse wave velocities, and it can be seen that the low wave number component in the velocity model is mainly updated in the elastic reflected wave travel time inversion based on the pure wave continuation equation, the pseudo-well curves of the real velocity, the initial velocity and the inversion velocity at different transverse distances are respectively extracted and compared with fig. 8, and the updated velocity pseudo-well curve is more fit with the real velocity, which indicates that the method can effectively recover the background velocity component with large scale. The results of reverse time migration imaging performed at the inversion speed are shown in fig. 9, and a comparison of the results in fig. 6 shows that after the low wave number component in the speed is accurately recovered, the reverse time migration results of the PP wave and the PS wave are improved to a higher degree, the imaging results are more focused, the arc drawing phenomenon of the diffraction body at the middle and deep layer position is eliminated, and the diffraction body is well migrated and returned.
The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims (8)

1. An elastic reflection wave travel time inversion method based on a pure wave continuation equation is characterized by comprising the following steps:
(1) omitting a cross term from the approximate elastic wave decoupling continuation equation to obtain a pure wave continuation equation, and calculating a background wave field based on the pure wave continuation equation;
(2) calculating a reverse time migration result by using the obtained background wave field;
(3) solving a disturbance pure wave continuation equation based on the pure wave continuation equation and the reverse time migration result, and calculating a disturbance wave field by using the disturbance pure wave continuation equation;
(4) respectively calculating longitudinal and transverse wave gradients based on the background wave field obtained in the step (1) and the disturbance wave field obtained in the step (3);
(5) and calculating the step size according to the longitudinal and transverse wave gradients, updating the longitudinal and transverse wave speeds, and outputting the updated longitudinal and transverse wave speeds as a final inversion result when the iteration value converges to the global minimum value.
2. The elastic reflection wave travel-time inversion method based on the pure wave continuation equation as claimed in claim 1, wherein the step (1) is as follows: the pure wave continuation equation comprises a pure P wave continuation equation and a pure S wave continuation equation, a forward-propagating P wave vibration velocity background wave field and a backward-propagating P wave vibration velocity background wave field are calculated based on the pure P wave continuation equation respectively, a backward-propagating S wave vibration velocity background wave field is calculated based on the pure S wave continuation equation, a backward-propagating P wave stress background wave field is calculated based on the pure P wave continuation equation, a backward-propagating S wave stress background wave field is calculated based on the pure S wave continuation equation, and numerical simulation calculation of the background wave field is implemented by high-order staggered grid finite difference.
3. The elastic reflection wave travel time inversion method based on the pure wave continuation equation as claimed in claim 2, wherein the pure P wave continuation equation is:
Figure FDA0002969635460000011
wherein the subscript "0" represents the background, x and z represent the spatial horizontal and vertical directions, respectively, t represents time, vPx0And vPz0Representing the horizontal and vertical components, tau, of the P-wave vibration velocity background wave field, respectivelyPxx0And τPzz0Representing the P-wave stress background wavefield, P, in xx and zz directions, respectively0And VP0Respectively representing background density and background longitudinal wave velocity;
the pure S wave continuation equation is:
Figure FDA0002969635460000021
wherein v isSx0And vSz0Representing the horizontal and vertical components, tau, of the S-wave vibration velocity background wave field, respectivelySxz0And τSzx0S-wave stress background wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity.
4. The elastic reflection wave travel-time inversion method based on the pure wave continuation equation as claimed in claim 2, wherein the step (2) is specifically as follows: and respectively calculating a PP wave reverse time migration result by cross-correlation of the forward-transmitted P wave vibration velocity background wave field and the backward-transmitted P wave vibration velocity background wave field, and calculating a PS wave reverse time migration result by using the forward-transmitted P wave vibration velocity background wave field and the backward-transmitted S wave vibration velocity background wave field.
5. The elastic reflection wave travel-time inversion method based on the pure wave continuation equation is characterized in that the step (3) is as follows: establishing a Born forward equation, namely a disturbed pure wave continuation equation, by taking the interaction of the background wave field and the velocity disturbance as a seismic source, wherein the disturbed pure wave continuation equation comprises a disturbed pure P wave continuation equation and a disturbed pure S wave continuation equation; wherein, the speed disturbance is replaced by a reverse time migration result;
and calculating a forward-transmitted P wave vibration velocity disturbance wave field and a backward-transmitted P wave stress disturbance wave field based on a disturbed pure P wave continuation equation, and calculating a forward-transmitted S wave vibration velocity disturbance wave field based on a disturbed pure S wave continuation equation.
6. The elastic reflection wave travel time inversion method based on the pure wave continuation equation as claimed in claim 5, wherein the disturbed pure P wave continuation equation is:
Figure FDA0002969635460000022
where the symbol "δ" represents a disturbance, ρ0Represents the background density, δ VPRepresenting longitudinal wave velocity disturbances, δ VSRepresenting shear wave velocity perturbations, vPx0And vPz0Respectively representing the horizontal and vertical components, V, of the P-wave vibration velocity background wave fieldP0Representing background longitudinal wave velocity, δ vPxAnd δ vPzRespectively representing the horizontal and vertical components, delta tau, of the P-wave vibration velocity disturbance wavefieldPxxAnd δ τPzzRespectively representing P wave stress disturbance wave fields in xx and zz directions; in practice, δ V in equation (3)PAnd δ VSExpressed by PP wave reverse time migration result;
the pure S wave continuation equation of the disturbance is as follows:
Figure FDA0002969635460000031
in practice, δ V in equation (4)SExpressed as a result of reverse time shift of the PS wave, δ vSxAnd δ vSzRespectively representing the horizontal and vertical components, delta tau, of the S-wave vibration velocity disturbance wavefieldSxzAnd δ τSzxS-wave stress disturbance wavefields, V, representing xz and zx directions, respectivelyS0Representing the background shear wave velocity.
7. The elastic reflection wave travel time inversion method based on the pure wave continuation equation as claimed in claim 1, wherein the specific method in the step (4) is as follows: reading the background wave field calculated in the step (1) and the disturbance wave field calculated in the step (3), and respectively solving the velocity gradient of the longitudinal wave
Figure FDA0002969635460000032
And transverse wave velocity gradient
Figure FDA0002969635460000033
Maintaining background density ρ throughout the calculation0Is a constant;
Figure FDA0002969635460000034
wherein, χPHexix-SRespectively representing the objective functions, V, of the matched P-wave data and the matched S-wave dataP0Representing background longitudinal wave velocity, vPx0And vPz0Respectively representing the horizontal and vertical components, deltav, of the forward propagating P-wave vibration velocity background wavefieldPxAnd δ vPzRespectively representing the horizontal component and the vertical component of the forward propagating P-wave vibration velocity disturbance wave field,
Figure FDA0002969635460000035
and
Figure FDA0002969635460000036
representing the xx and zz components of the backward propagating P-wave stress perturbation wavefield,
Figure FDA0002969635460000037
and
Figure FDA0002969635460000038
representing the xx and zz components, V, respectively, of the backward propagating P-wave stress background wavefieldS0Representing the background shear wave velocity, δ vSxAnd δ vSzRespectively representing the horizontal component and the vertical component of the forward propagating S-wave vibration velocity perturbation wave field,
Figure FDA0002969635460000039
and
Figure FDA00029696354600000310
representing the xz and zx components, respectively, of the back-propagating S-wave stress background wavefield.
8. The elastic reflection wave travel time inversion method based on the pure wave continuation equation is characterized in that in the step (5), the update amount of the velocity of the longitudinal wave in each iteration is ensured to be between 20m/s and 100m/s, and the update amount of the velocity of the transverse wave is ensured to be between 10m/s and 60 m/s.
CN202110260212.3A 2021-03-10 2021-03-10 Elastic reflection wave travel time inversion method based on pure wave continuation equation Active CN113126151B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110260212.3A CN113126151B (en) 2021-03-10 2021-03-10 Elastic reflection wave travel time inversion method based on pure wave continuation equation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110260212.3A CN113126151B (en) 2021-03-10 2021-03-10 Elastic reflection wave travel time inversion method based on pure wave continuation equation

Publications (2)

Publication Number Publication Date
CN113126151A true CN113126151A (en) 2021-07-16
CN113126151B CN113126151B (en) 2022-06-07

Family

ID=76772859

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110260212.3A Active CN113126151B (en) 2021-03-10 2021-03-10 Elastic reflection wave travel time inversion method based on pure wave continuation equation

Country Status (1)

Country Link
CN (1) CN113126151B (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101358867A (en) * 2008-09-03 2009-02-04 山东省科学院海洋仪器仪表研究所 Ocean water level real-time monitoring system
US20090257308A1 (en) * 2008-04-11 2009-10-15 Dimitri Bevc Migration velocity analysis methods
CN103293553A (en) * 2013-04-17 2013-09-11 中国海洋石油总公司 Continuation and correction method for boundary element of earthquake data collected through upper cables and lower cables in complex seabed
CN108873066A (en) * 2018-06-26 2018-11-23 中国石油大学(华东) Elastic fluid fluctuates equation back wave Travel Time Inversion method
CN111983703A (en) * 2020-07-24 2020-11-24 中国石油天然气集团有限公司 Method, system and device for imaging fluid through interwell electromagnetic measurement

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090257308A1 (en) * 2008-04-11 2009-10-15 Dimitri Bevc Migration velocity analysis methods
CN101358867A (en) * 2008-09-03 2009-02-04 山东省科学院海洋仪器仪表研究所 Ocean water level real-time monitoring system
CN103293553A (en) * 2013-04-17 2013-09-11 中国海洋石油总公司 Continuation and correction method for boundary element of earthquake data collected through upper cables and lower cables in complex seabed
CN108873066A (en) * 2018-06-26 2018-11-23 中国石油大学(华东) Elastic fluid fluctuates equation back wave Travel Time Inversion method
CN111983703A (en) * 2020-07-24 2020-11-24 中国石油天然气集团有限公司 Method, system and device for imaging fluid through interwell electromagnetic measurement

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
张晓语，等: "基于能量密度的自解耦互相关成像条件", 《地球物理学报》 *

Also Published As

Publication number Publication date
CN113126151B (en) 2022-06-07

Similar Documents

Publication Publication Date Title
CN108873066B (en) Elastic medium wave equation reflected wave travel time inversion method
KR102020759B1 (en) Q-compensated full wave field reversal
Etgen et al. Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial
US20080192572A1 (en) Method for processing borehole seismic data
AU2015383134B2 (en) Multistage full wavefield inversion process that generates a multiple free data set
CN110007340B (en) Salt dome velocity density estimation method based on angle domain direct envelope inversion
WO2017162731A1 (en) Method of operating a data-processing system for the simulation of the acoustic wave propagation in the transversely isotropic media comprising an hydrocarbon reservoir
Shabelansky et al. Converted-wave seismic imaging: Amplitude-balancing source-independent imaging conditions
CN110187382A (en) A kind of diving Wave and back wave wave equation Travel Time Inversion method
CN109946742A (en) The pure rolling land qP shakes digital simulation method in a kind of TTI medium
Datta et al. Full-waveform inversion of salt models using shape optimization and simulated annealing
CN114460646A (en) Wave field excitation approximation-based reflection wave travel time inversion method
NO20190489A1 (en) Seismic modeling
CN107340537A (en) A kind of method of P-SV converted waves prestack reverse-time depth migration
Jia et al. Superwide-angle one-way wave propagator and its application in imaging steep salt flanks
CN113126151B (en) Elastic reflection wave travel time inversion method based on pure wave continuation equation
Jin et al. 2D multiscale non‐linear velocity inversion
Lin et al. Time-domain wavefield reconstruction inversion solutions in the weighted full waveform inversion form
WO2017189180A1 (en) Fwi with areal and point sources
Jiang et al. Full waveform inversion based on inversion network reparameterized velocity
Moura et al. Progressive matching optimisation method for FWI
Saraswat et al. Simultaneous stochastic inversion of prestack seismic data using hybrid evolutionary algorithm
Shin et al. Laplace-domain full waveform inversion using irregular finite elements for complex foothill environments
Iversen Velocity rays for heterogeneous anisotropic media: Theory and implementation
CN114624766B (en) Elastic wave least square reverse time migration gradient solving method based on traveling wave separation

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