Speed and Q wave equation inversion method based on local domain travel time

Abstract

The invention discloses a speed and Q wave equation inversion method based on local domain travel time, and relates to the technical field of oil exploration, and the method comprises the steps: carrying out the viscoacoustic forward modeling to obtain seismic data; calculating travel time difference and obtaining a time difference accompanying source; calculating a back propagation wave field and gradient; and updating the speed model and the Q model. The method is based on wave equation travel time inversion, travel time information is utilized, the cycle jump problem of full waveform inversion under the condition of large time difference is avoided, and dependence on an initial model is reduced; the DTW algorithm is used for calculating the local domain travel time difference, the corresponding time difference can be accurately calculated, the situation that a cross-correlation algorithm tends to a large-amplitude event under the condition of large time difference is avoided, and the inversion precision is improved; according to the method, the speed and the Q can be updated at the same time by using travel time information, a large amount of Fourier transform calculation in a frequency domain method is avoided, and the inversion efficiency is greatly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of oil exploration, and particularly relates to a velocity and Q wave equation inversion method based on local domain travel time. BACKGROUND

[0002] Real strata exist viscosity, which will affect the propagation of seismic waves, resulting in amplitude attenuation and phase distortion. This attenuation characteristic is usually quantified by a quality factor Q, and a low Q value means that the energy loss or attenuation of each periodic wave is large, and if these attenuations are not compensated, it will lead to reduced imaging resolution and inaccurate imaging position, especially in oil and gas reservoir areas. In order to obtain accurate reservoir information, absorption attenuation compensation must be carried out to obtain high-quality imaging results, and the premise of high-resolution imaging is to obtain high-resolution velocity and Q models.

[0003] At present, the method for obtaining velocity model and Q model is mainly based on multi-parameter full waveform inversion (Full-waveform inversion, FWI), and a series of inversion strategies are developed on this basis. However, due to the inaccuracy of the initial velocity, the residual objective function has a highly nonlinear relationship with the velocity, and in the case of poor initial velocity, there is a large travel time difference between the observed seismic record and the simulated seismic record, which is easy to cause the cycle skipping problem, making the inversion fall into a local minimum value. In addition, since velocity and Q are a set of strongly correlated parameters, they significantly influence each other in the inversion process, and any Q tomography method is highly dependent on the accuracy of the velocity model. The mainstream industrial workflow adopts a separate inversion method, but the correlation between velocity and Q will result in insufficient accuracy of the separate inversion result, which brings difficulties to oil and gas exploration.

[0004] Travel time is linearly related to long-wavelength velocity, so a preliminary velocity model is usually obtained by using tomographic inversion of travel time information, which mainly includes ray-based methods and wave equation-based methods. Due to the limitations of ray-based methods such as high-frequency assumption, wave equation travel time inversion can not only cope with complex wave field conditions, but also directly calculate the travel time difference from the seismic record without additional travel time calculation, greatly improving the inversion efficiency. However, in order to obtain more accurate time difference values, the cross-correlation method needs to be windowed for calculation, which will make the time difference value deviate to the large amplitude phase axis in the time window, resulting in reduced resolution of the inversion result.

[0005] How to obtain high-resolution velocity and Q models is one of the problems to be solved in oil exploration. SUMMARY

[0006] To solve the above technical problems, the application discloses a velocity and Q wave equation inversion method based on local domain travel time, which is based on wave equation travel time inversion, combines a time difference target function and a SLS-based viscoacoustic wave equation through a connection function, derives a time difference accompanying source, and obtains a gradient about velocity and Q.

[0007] To achieve the above object, the application adopts the following technical scheme:

[0008] A velocity and Q wave equation inversion method based on local domain travel time comprises the following steps:

[0009] Step one, viscoacoustic forward modeling to obtain seismic data

[0010] Based on the SLS viscoacoustic wave equation, Q display is introduced into the wave equation, a real velocity model and a real Q model are used to forward modeling to obtain observed seismic records, an initial velocity model and an initial Q model are used to forward modeling to obtain the first iteration of simulated seismic records, and the positive wave field part in the gradient is saved;

[0011] Step two, calculating travel time difference and obtaining time difference accompanying source

[0012] The local domain travel time difference between the observed seismic records and the simulated seismic records is calculated by using the DTW algorithm, and the time difference accompanying source is calculated by combining the derived formula;

[0013] Step three, calculating reverse wave field and gradient

[0014] The accompanying source is brought into the accompanying equation to obtain the reverse wave field, and the velocity and Q gradient is obtained by combining the positive wave field;

[0015] Step four, updating velocity model and Q model

[0016] The velocity and Q are updated by the gradient, and the new velocity and Q model are used to forward modeling to obtain new simulated seismic records, and the final velocity model and Q model are obtained by repeating steps two to four.

[0017] Optionally, in step one, based on the SLS viscoacoustic wave equation, Q display is introduced into the wave equation, and the formula is as follows:

[0018] ;

[0019] ;

[0020] ;

[0021] wherein, represents pressure field, represents velocity vector, represents divergence, represents density, represents memory variable related to quality factor Q, represents source function at, , , , and are parameters related to velocity and Q respectively, and their relationship formulas are as follows:

[0022] ;

[0023] ;

[0024] ;

[0025] ;

[0026] wherein, represents frequency.

[0027] Optionally, in step two, based on the adjoint state method, the corresponding adjoint equation is derived, and the formula is as follows:

[0028] ;

[0029] ;

[0030] ;

[0031] ;

[0032] wherein, represents time difference adjoint source, , , , corresponds to , , , in wave equation one by one.

[0033] Optionally, in step two, the DTW algorithm is used instead of the traditional cross-correlation method to calculate the travel time difference between the observed seismic record and the simulated seismic record, and the time difference value is corresponded to each point in the seismic record, so as to accurately calculate the local domain travel time difference.

[0034] Optionally, in step three, the gradient of the time difference objective function with respect to velocity and Q is derived respectively, and the formula is as follows:

[0035] ;

[0036] ;

[0037] ;

[0038] wherein, represents a target function, represents a travel time difference between an observed seismic record and a simulated seismic record, represents a gradient of accelerates the update of Q, represents a velocity gradient, represents a bulk modulus, represents a divergence of a velocity vector, is a forward wave field part obtained and saved in step one through a visco-acoustic wave equation forward process, and is a backward wave field part obtained in step three through calculation in an adjoint equation, and the gradient is calculated simultaneously in the adjoint equation calculation process.

[0039] In the process of deriving the gradient, the form of the travel time adjoint source is obtained, and the formula is as follows:

[0040] ;

[0041] wherein, represents a derivative of the observed seismic record with respect to time after adding the travel time.

[0042] The present application has the beneficial effects that, the present application connects the travel time information and the visco-acoustic wave equation through the wave equation travel time inversion, based on the linear correlation between the velocity and the long wavelength velocity, simultaneously inverts the velocity and the Q model in the case of a poor initial model, and uses the DTW algorithm to calculate the travel time difference between the observed seismic record and the simulated seismic record, compared with the cross-correlation algorithm, more accurate travel time values can be obtained, and the inversion precision is improved; finally, the simulation experiment using the overthrust model proves the accuracy of the method of the present application.

[0043] The method of the present application simultaneously performs the velocity and Q tomographic inversion using the travel time information, avoids the dependence of the full waveform multi-parameter inversion on the initial model, simultaneously inverts the Q model using the travel time information, avoids the large amount of Fourier calculation in the traditional frequency domain method, reduces the calculation cost, and in addition, the DTW algorithm is used to calculate the travel time difference, more accurate travel time values can be obtained, and the inversion precision is improved. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1A flowchart of a velocity and Q wave equation inversion method based on local domain travel time according to an embodiment of the present application;

[0045] Figure 2A A real velocity model according to an embodiment of the present application;

[0046] Figure 2B A real Q model according to an embodiment of the present application;

[0047] Figure 3A An initial inversion velocity model according to an embodiment of the present application;

[0048] Figure 3B An initial Q model according to an embodiment of the present application;

[0049] Figure 4 Forward modeling seismic records according to an embodiment of the present application;

[0050] Figure 5A A travel time calculation result according to an embodiment of the present application;

[0051] Figure 5B A travel time accompanying source according to an embodiment of the present application;

[0052] Figure 6A A final inversion velocity result according to an embodiment of the present application;

[0053] Figure 6B A final inversion Q result according to an embodiment of the present application;

[0054] Figure 7A An inversion velocity and real velocity trace comparison according to an embodiment of the present application;

[0055] Figure 7B An inversion Q and real Q trace comparison according to an embodiment of the present application. DETAILED DESCRIPTION

[0056] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application but not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of the present application. Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed application, but only represents selected embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of the present application.

[0057] A velocity and Q wave equation inversion method based on local domain travel time, such as Figure 1 As shown, the following steps are included:

[0058] Step 1: Obtain seismic data using viscoacoustic forward modeling

[0059] In wave equation inversion, in order to simulate the viscosity of the actual formation, it is necessary to introduce the Q parameter into the wave equation. Different theories have different forms of viscoacoustic wave equations. Based on the viscoacoustic wave equation of the SLS model, this paper explicitly introduces the quality factor Q into the wave equation. The formula is as follows:

[0060] ;

[0061] ;

[0062] ;

[0063] Where, represents the pressure field, represents the velocity vector, represents the divergence, represents density, represents the memory variable associated with the quality factor Q, Representatives in The source function at 、 、 and They are parameters related to speed and Q, and their relationship formulas are as follows:

[0064] ;

[0065] ;

[0066] ;

[0067] ;

[0068] Where, Represents frequency.

[0069] Figure 2A is the real nappe velocity model, Figure 2B To correspond to the Q model, the red area in the figure represents the low Q value area, which is actually the gas reservoir area with strong attenuation. According to the velocity model and Q model, the observed seismic records are simulated by the above viscoacoustic wave equation, as shown in the following figure: Figure 4 shown.

[0070] Figure 3AThe initial velocity model is obtained by smoothing the real model on a large scale. It is linear in the longitudinal direction and simulates the initial velocity model obtained in actual production. Figure 3B is the initial Q model, which is a uniform model and simulates the unknown gas reservoir area in actual production.

[0071] Step 2: Calculate the travel time difference and obtain the accompanying source of the time difference

[0072] The DTW algorithm is used to calculate the local domain travel time difference between the observed seismic records and the simulated seismic records, and the derived formula is used to calculate the time difference accompanying source; the key to wave equation travel time inversion lies in the calculation of travel time difference. In order to obtain accurate travel time difference, windowing is generally used for cross-correlation calculation. However, when the time difference is large, the time difference value will be biased towards the large-amplitude phase axis in the time window, resulting in a decrease in inversion accuracy.

[0073] In order to obtain accurate travel time difference, the present invention adopts DTW algorithm to calculate, such as Figure 5A As shown, each point in the data domain corresponds to its time difference value, rather than calculating a time difference for each time window in the windowing calculation, and obtaining the time difference accompanying source ( Figure 5B ), its form is obtained from the gradient derivation process, the formula is as follows:

[0074] ;

[0075] Where, represents the time difference companion source, represents the time derivative of the observed seismic record after adding the time difference, Represents the travel time difference between the observed earthquake record and the simulated earthquake record.

[0076] Step 3: Calculate the back propagation field and gradient

[0077] The traditional wave equation travel time inversion method is only applied to the acoustic wave equation, but both velocity and Q will affect the travel time. Therefore, it is applied to the viscoacoustic wave equation and the gradient of the time difference objective function with respect to velocity and Q is derived respectively. The formula is as follows:

[0078] ;

[0079] ;

[0080] ;

[0081] Where, represents the objective function, represents the travel time difference between the observed earthquake record and the simulated earthquake record, represent The gradient of Accelerate the update of Q, representing the velocity gradient, representing the bulk modulus, representing the divergence of the velocity vector, is the forward wavefield part obtained and saved in step one through the process of forward modeling of viscoacoustic wave equation, and is the backward wavefield part obtained in step three through the calculation of the adjoint equation, and the gradient is calculated simultaneously in the process of adjoint equation calculation.

[0082] Step four, updating the velocity model and Q model

[0083] The velocity and Q are updated by the gradient, and the iteration is repeated to obtain the final velocity model (v) and Q model (Q). Figure 6A Through the comparison of the extracted channels (v, Q), it can be seen that the inverted velocity curve is highly consistent with the true velocity curve, and the inverted Q model curve correctly reflects the position of the attenuation area. Figure 6B Through this method of inversion, a high-precision velocity model and Q model are obtained in the initial model which is a linear velocity model (v) and a uniform Q model (Q), providing an initial model for the application of full waveform inversion in actual production. Figure 7A Figure 7B The method of the application is based on the travel time inversion theory of wave equation, and a velocity and Q wave equation inversion method based on local domain travel time is developed. Figure 3A Figure 3B Through the viscoacoustic wave equation under the SLS model, the gradient of the travel time objective function to the velocity and Q is derived, the travel time adjoint source is obtained, so that the velocity and Q model can be inverted simultaneously under the condition of poor initial model.

[0084] The method of the application is based on the travel time inversion theory of wave equation, and a velocity and Q wave equation inversion method based on local domain travel time is developed.

[0085] Of course, the above description is not a limitation of the application, and the application is not limited to the above examples.Changes, modifications, additions or substitutions made by those skilled in the art within the scope of the application should also be within the protection scope of the application.​

Claims

1. A velocity and Q wave equation inversion method based on local domain travel time, characterized in that: The following steps are involved: Step 1: Obtain seismic data using viscoacoustic forward modeling Based on the viscoacoustic wave equation of the standard linear solid attenuation model, the Q display is introduced into the wave equation. Through the finite difference method, the observed seismic records are forward modeled using the true velocity model and the true Q model. The simulated seismic records of the first iteration are forward modeled using the initial velocity model and the initial Q model, and the forward wave field part of the gradient is saved. Step 2: Calculate the travel time difference and obtain the accompanying source of the time difference The DTW algorithm is used to calculate the local domain travel time difference between the observed earthquake record and the simulated earthquake record, and the derived formula is used to calculate the time difference companion source; Step 3: Calculate the back propagation field and gradient Substitute the companion source into the companion equation to obtain the reverse wave field, and combine it with the forward wave field to obtain the gradient of velocity and Q; Step 4: Update velocity model and Q model The velocity and Q are updated by gradient correspondence, and the new velocity and Q models are used to forward model new simulated earthquake records. Steps 2 to 4 are repeated to obtain the final velocity model and Q model.

2. The velocity and Q wave equation inversion method based on local domain travel time according to claim 1, characterized in that: In step 1, based on the viscoacoustic wave equation of the standard linear solid attenuation model, the Q display is introduced into the wave equation, and the formula is as follows: ； ； ； Where, represents the pressure field, represents the velocity vector, represents the divergence, represents density, represents the memory variable associated with the quality factor Q, Representatives in The source function at 、 、 and They are parameters related to speed and Q, and their relationship formulas are as follows: ； ； ； ； Where, Represents frequency.

3. The velocity and Q wave equation inversion method based on local domain travel time according to claim 2, characterized in that: In step 2, based on the adjoint state method, the corresponding adjoint equation is derived, and the formula is as follows: ； ； ； ； Where, represents the time difference companion source, 、 、 、 and the wave equation 、 、 、 One to one correspondence.

4. The velocity and Q wave equation inversion method based on local domain travel time according to claim 3, characterized in that: In step 3, the gradients of the time difference objective function with respect to velocity and Q are derived respectively, and the formulas are as follows: ； ； ； Where, represents the objective function, represents the travel time difference between the observed earthquake record and the simulated earthquake record, represent The gradient of Accelerate the update of Q, represents the velocity gradient, represents the bulk modulus, Represents the divergence of the velocity vector, which is the forward wave field part obtained and saved in the forward modeling process of the viscoacoustic wave equation in step 1. and It is the part of the back propagation field obtained by calculating the adjoint equation in step 3, and the gradient is calculated simultaneously during the calculation of the adjoint equation; In the process of deriving the gradient, the form of the time difference companion source is obtained, and the formula is as follows: ； Where, Represents the time derivative of the observed seismic record after adding the time difference.

Images