A wave equation full-wave q tomography method

Abstract

The application provides a wave equation full-wave Q tomography method, which comprises the following steps: performing FFT transformation on existing seismic records to determine the frequency of a Ricker wavelet; performing forward simulation using the frequency of the Ricker wavelet determined in step 1 to obtain seismic correlation data; converting the seismic correlation data to a local domain and performing FFT on the data in the local domain to obtain the frequency of a single seismic event; obtaining the peak frequency shift of each seismic event in the local domain according to the frequency of the single seismic event; obtaining a companion source according to the peak frequency shift; performing wave field continuation calculation according to the companion source to obtain a back-propagating wave field; obtaining a gradient using the back-propagating wave field; and performing iterative updating using a conjugate gradient method according to the gradient to finally obtain an updated Q model. The method provided by the application can be simulated to obtain a relatively accurate Q model.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic exploration technology, and particularly relates to a full-wave Q tomography method based on wave equation. BACKGROUND

[0002] The earth is inelastic, which distorts the amplitude and phase of the propagating seismic wave. The attenuation of seismic wave can be quantified by a quality factor Q, which indicates the phase shift and amplitude loss as a function of the frequency content and travel distance of the propagating wave. A lower Q value means a greater loss of energy or a greater attenuation of each periodic wave, which will lead to a lower imaging resolution and an inaccurate imaging position if not compensated. Therefore, how to accurately solve these influences is crucial to the quality of imaging.

[0003] Q compensation migration algorithm can eliminate these unnecessary disturbances to obtain a higher resolution migration image. However, all Q migration algorithms need a relatively accurate Q model. The traditional wave equation Q tomography can also invert the Q model by eliminating the peak frequency difference between the observed early wave and the artificially synthesized early wave, but can only tomographically invert the Q model near the surface or in a large range, and cannot accurately estimate the Q value of each stratum. Using full waveform inversion can also obtain a higher precision Q model, but the computational efficiency is too low. SUMMARY

[0004] In view of the above problems, the present application aims to provide a full-wave Q tomography method based on wave equation, which can more accurately calculate the underground Q value and is suitable for early wave and full wave respectively.

[0005] To achieve the above-mentioned purpose, the present application comprises the following steps:

[0006] A full-wave Q tomography method based on wave equation comprises the following steps:

[0007] Step 1: According to the existing seismic record, perform FFT transformation to obtain the main frequency and determine the frequency of the Ricker wavelet;

[0008] Step 2: According to the velocity model and the initial Q model for tomography, perform forward modeling using the Ricker wavelet frequency determined in step 1 to obtain seismic correlation data;

[0009] Step 3: Convert the seismic correlation data to a local domain and perform FFT on the data in the local domain to obtain the frequency of a single seismic event;

[0010] Step 4: According to the frequency of the single seismic event, obtain the peak frequency shift of each seismic event in the local domain;

[0011] Step 5: Obtain the ghost source according to the peak frequency shift.

[0012] Step 6: According to the accompanying source, the accurate wave field continuation calculation is carried out to obtain the back propagation wave field q(x s ,t,x r ) and s(x s ,t,x r );

[0013] Step 7: The gradient is obtained by using the back propagation wave field through the cross correlation of the forward propagation wave field and the back propagation wave field;

[0014] Step 8: The conjugate gradient method is used according to the gradient to carry out iterative updating, and finally the updated Q model is obtained.

[0015] Further, the step 2 of the wave equation full wave Q tomography method as described above comprises: using the Ricker wavelet frequency and formula (1), formula (2) to carry out forward simulation:

[0016]

[0017] Wherein, P is the pressure field data, v = {v x ,v z} represents the particle velocity vector, r p represents the memory variable, K represents the bulk modulus of the medium, S(x s ,t) represents that the source is located at x = x s , the stress strain relaxation parameters τ, τ ε and τ σ are related to Q and angular frequency ω, and their relationship is as follows:

[0018] The relationship is as follows:

[0019]

[0020] Further, the step 3 of the wave equation full wave Q tomography method as described above comprises:

[0021] Comprise:

[0022] According to formula (3), the seismic related data is converted to the local domain and the data in the local domain is subjected to FFT to obtain the frequency of a single seismic event;

[0023] S(t,k)=s(t)w(t-k),S f (f,t)=FFT(S(k,t)), (3)

[0024] Wherein, S and s represent seismic records in the local domain and the data domain, w represents a Gaussian time window, S f represents the local frequency, which contains the frequency of each seismic event.

[0025] Further, as described above, the wave equation based full-wave Q tomography method, step 4

[0026] comprises: calculating the peak frequency shift of each seismic event in the local domain by using formula (4) and formula (5):

[0027]

[0028]

[0029] wherein, is the local spectrum of the observed data, is the local spectrum of the predicted data, f is frequency, f s is the weight coefficient, the weight coefficient is the two norm of the cross-correlation product f s representing the minimum f s is the peak frequency shift; C(f

[0030] Further, as described above, the wave equation based full-wave Q tomography method, step 5 comprises: using formula (6) to accompany the source:

[0031] ΔP(x s ,t,x r )=P(x s ,t,x r )*Δf(x s ,t,x r )*C. (6)

[0032] wherein, ΔP is the accompanied source, X s represents the position of the source, X r represents the position of the receiver, t represents time; C represents the cross-correlation of the spectrum of the observed data and the spectrum of the predicted data in the local domain.

[0033] Further, as described above, the wave equation based full-wave Q tomography method, step 6

[0034] comprises: obtaining the back-propagating wave field q(x s ,t,x r ) and s(x s ,t,x r ) by accurate wave field continuation calculation by formula (7):

[0035]

[0036] wherein ΔP is the accompanied source is the source term, q and u are the accompanied vectors of P and v, s is the accompanied vector of r p .

[0037] Further, the step 7 of the full-wave Q tomography method based on wave equation as described above comprises:

[0038] The gradient is calculated by using the back-propagation wave field and formula (8);

[0039]

[0040] where, v(x s ,t,x r ) is the particle velocity vector in the forward process, q(x s ,t,x r ) is the particle velocity vector in the back-propagation process, and s(x s ,t,x r ) is the memory variable in the back-propagation process.

[0041] The present application has the following advantages by adopting the above technical solution:

[0042] (1) In deep underground, if there is a gas layer or a strong attenuation area, the conventional early-wave-based wave equation Q tomography cannot be simulated, while the full-wave-based Q tomography method proposed in the present application can be simulated and a relatively accurate Q model can be obtained;

[0043] (2) A new method for calculating the peak frequency in a local domain is proposed, which can avoid crosstalk noise caused by surrounding seismic events;

[0044] (3) By using all seismic events, a high-resolution Q model can be obtained, and the tomographic Q model can be used as a background model for any migration compensation imaging. The method provided in the present application can be widely applied in the field of migration imaging in strong attenuation areas such as gas reservoirs. BRIEF DESCRIPTION OF DRAWINGS

[0045] Figure 1 is a flow chart of the full-wave Q tomography method based on wave equation of the present application;

[0046] Figure 2 is a schematic diagram of a method for calculating the peak frequency shift of a seismic event in a local domain, wherein (2a) is two simple seismic records each containing different peak frequency shifts, and (2b) is the peak frequency shift calculated by using the method of the present application;

[0047] Figure 3 is a model diagram for proving the correctness of the method of the present application; wherein (3a) is a velocity model, (3b) is a real quality factor Q model, and (3c) is a tomographic initial model;

[0048] Figure 4Comparison of the results of the tomography method of the present application and conventional wave equation Q tomography: (4a) is conventional wave equation Q tomography, (4b) is the tomography method of the present application;

[0049] Figure 5 Comparison of the results of the tomography method of the present application and conventional tomography in a single channel;

[0050] Figure 6 Four different reverse-time migration results: (6a) acoustic reverse-time migration; (6b) acoustic reverse-time migration with attenuation; (6c) viscoacoustic reverse-time migration using a real Q model for compensation; (6d) Q model obtained by the method of the present application for acoustic reverse-time migration. DETAILED DESCRIPTION

[0051] In order to make the purpose, technical solutions and advantages of the present application clearer, the technical solutions in the present application will be described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0052] The wave equation-based Q model modeling method provided by the present application first converts seismic data to a local domain by using a sliding Gaussian window to reduce crosstalk noise between different seismic events, so as to realize the goal of calculating the frequency of a single seismic event. Then, the cross-correlation algorithm of the frequency spectrum of the observed data and the synthetic data is used to calculate the peak frequency shift of each seismic event in the local domain. The gradient is obtained by calculating the companion source and performing reverse continuation, so that a higher-precision Q model can be obtained by using all seismic events in a seismic record. Finally, a simulation experiment is performed using a BP gas model to prove that a higher-resolution Q model can be obtained, and an imaging image is obtained by using a Q migration compensation algorithm, thereby further proving the accuracy of the Q tomography imaging method proposed by the present application. As shown in FIG. 1, the method specifically comprises the following steps: Figure 1

[0053] Step 1: According to the existing seismic record, perform FFT transformation on the seismic record, calculate the main frequency, and determine the frequency of the Ricker wavelet used in the present application;

[0054] Specifically, the frequency of the Ricker wavelet is ultimately calculated by performing FFT transformation on the existing seismic record, which can make the tomography process converge faster.

[0055] Step 2: According to the velocity model and the initial Q model for tomography, use the Ricker wavelet frequency determined in step 1 and formulas (1) and (2) to perform forward modeling,

[0056]

[0057] ​​acquiring seismic related data;

[0058]

[0059] where P is pressure field data, v = {v x ,v z} represents particle velocity vector, r p represents memory variable, K represents bulk modulus of medium, S(x s ,t) represents source located at x = x s .

[0060] Stress-strain relaxation parameters τ, τ ε and τ σ are related to Q and angular frequency ω, and their relationship is shown in the following formula.

[0061]

[0062] Specifically, the purpose of the step of the present application is to acquire seismic records without attenuation, and then acquire the difference between the seismic records and the real seismic records, so as to make the tomography result more accurate.

[0063] Step 3: According to formula (3), the seismic related data is converted to a local domain and the data therein is subjected to FFT to acquire the frequency of a single seismic event;

[0064] S(t,k) = s(t)w(t-k), S f (f,t) = FFT(S(k,t)), (3)

[0065] where S and s represent seismic records in a local domain and a data domain, w represents a Gaussian time window, S f represents local frequency (including the frequency of each seismic event);

[0066] Specifically, the purpose of the step is to acquire the frequency of a single seismic event, and compared with the traditional FFT for acquiring the main frequency, the method provided by the present application can avoid the crosstalk noise caused by adjacent seismic events in the seismic records.

[0067] Step 4: According to the frequency of the single seismic event, the peak frequency shift of each seismic event in the local domain is calculated by using formula (4) and formula (5);

[0068]

[0069]

[0070] where is the local spectrum of the observed data, is the local spectrum of the predicted data, f is frequency, and fs is the weight coefficient, and the two-norm of the cross-correlation product f represents the peak frequency shift s is the peak frequency shift; C(f s ,t) represents the cross-correlation of the frequency spectrum of the observed data and the predicted data in the local domain;

[0071] Specifically, the purpose of this step is to obtain the peak frequency shift of a single seismic event, so as to make the tomography result more accurate.

[0072] Step 5: Obtain the companion source according to the peak frequency shift using formula (6);

[0073] ΔP(x s ,t,x r )=P(x s ,t,x r )*Δf(x s ,t,x r )*C. (6)

[0074] where ΔP is the companion source, X s represents the position of the seismic source, X r represents the position of the receiver, and t represents time; C represents the cross-correlation of the frequency spectrum of the observed data and the predicted data in the local domain;

[0075] Specifically, the purpose of this step is to obtain the companion source for constructing the counter-propagating wave field, so that the computational efficiency can be improved using the companion state method.

[0076] Step 6: According to the companion source, accurate wave field continuation calculation is performed through formula (7) to obtain the counter-propagating wave fields q(x s ,t,x r ) and s(x s ,t,x r );

[0077]

[0078] where ΔP is the companion source is the source term, q and u are the companion vectors of P and v, and s is the companion vector of r p ;

[0079] Step 7: By cross-correlating the forward-propagating wave field and the counter-propagating wave field, the gradient is obtained using the counter-propagating wave field and formula (8);

[0080]

[0081] where v(x s ,t,x r ) is the particle velocity vector in the forward process, q(xs ,t,x r ) is the particle velocity vector during the backward propagation process, s(x s ,t,x r ) is the memory variable during the backward propagation process;

[0082] Specifically, the purpose of the step is to obtain the gradient for updating the Q model; so that the result converges faster.

[0083] Step 8: iteratively update using the conjugate gradient method according to the gradient, and finally obtain the updated Q model.

[0084] Embodiment:

[0085] The full-wave Q tomography method based on wave equation provided by the embodiment of the application comprises the following steps:

[0086] 1) As shown in the model designed, we use a real BP gas model to prove the accuracy and effectiveness of the method proposed, and the velocity model is as shown in Figure a. There is a strong attenuation area on the top of the model, and the Q value is about 20.0, as shown in Figure b. The grid numbers in the x and z directions are 398 and 161, and the spatial sampling is 10m. Here, we use a Ricker wavelet to obtain the propagation wave field, and the frequency is 12Hz. The length of the seismic record and the time sampling interval are 2.8s and 0.8ms respectively. There are a total of 130 shots, and the first shot is located at (60m, 10m), and the shot interval is 30m. We set the initial Q value to 1000.0, as shown in Figure c. We set a Q value of 20 on the top to remove the direct wave. Figure 3 Figure 3 Figure 3 Figure 3 2) By using formula (1) and (2) for forward modeling, we can obtain two kinds of seismic records, one is obtained by the real Q model, and the other is calculated by the given initial Q model, and then use formula 3, 4 and 5 to calculate the peak frequency shift.

[0087] a is two seismic records representing observed data and predicted data respectively, and each data contains two seismic events. Figure 2 b is the frequency curve and the frequency shift curve, wherein the solid line represents the frequency of the predicted data, the dotted line represents the frequency of the observed data, and the dashed line represents the frequency difference. Figure 2 Figure 2

[0088] ​​​​​3) Put the result of step 2) into equation 6 to calculate the accompanying source, i.e. ΔP. The accompanying source obtained by the method of the present application contains more deep information than the traditional wave equation Q tomography. The reverse wave field is obtained by reverse continuation of the accompanying source, and then the gradient is calculated by equation (8).

[0089] 4) According to the result of step 3), the final Q model is obtained by several iterations. The result obtained by the method of the present application is better than the result obtained by the traditional method. The single channel curve is extracted for comparison (Fig. 3c). Figure 5

[0090] 5) The model obtained in step 4) is used to obtain the compensated reverse time migration result by using the compensation migration imaging technology. The present application provides four imaging results to verify the accuracy of the technology proposed by the present application.

[0091] 6) First Figure 6 a shows the acoustic reverse time migration result, which is used as a reference model, Figure 6 b is the acoustic reverse time migration with attenuation, Figure 6 c is the visco-acoustic reverse time migration compensated by using the real Q model, Figure 6 d is the visco-acoustic reverse time migration using the Q model obtained by the technology of the present application. It can be obviously seen that the amplitude anomaly caused by attenuation after compensation using the Q model obtained by our method is corrected, which is the same as the effect after compensation using the real Q model. Therefore, our technology is accurate and feasible.

[0092] Figure 4 Comparison of the results of the tomography method of the present application and the conventional wave equation Q tomography: (4a) is the conventional wave equation Q tomography, and (4b) is the tomography method of the present application. Through Figure 4 It can be seen that the tomography result obtained by our method accurately represents the strong attenuation area existing in the model, while the conventional method only tomographs the shallow part. For the deep part, our method can more clearly depict the underground profile, which is not achieved by the conventional method.

[0093] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, but not to limit them. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement to some technical features. The modification or replacement does not make the essence of the corresponding technical solution deviate from the spirit and scope of the technical solutions of the embodiments of the present application.​

Claims

1. A full-wave method based on the wave equation Q Chromatographic method, characterized in that, Includes the following steps: Step 1: Based on the existing seismic records, perform an FFT transformation on them to obtain the dominant frequency and determine the frequency of the Ricker wavelet; Step 2: Initial tomography based on the velocity model Q The model is used to perform forward modeling using the Rack wavelet frequency determined in step 1 to obtain earthquake-related data; Step 3: Convert the earthquake-related data into a local domain and perform FFT on the data to obtain the frequency of individual earthquake events; Step 4: Calculate the peak frequency shift of each seismic event in the local domain based on the frequency of the individual seismic events; Step 5: Obtain the accompanying source based on the peak frequency shift; Step 6: Based on the accompanying source, perform accurate wavefield extension calculations to obtain the reverse propagating wavefield; Step 7: Calculate the gradient using the anti-propagation wavefield by cross-correlation between the forward and reverse propagation wavefields; Step 8: Iteratively update the gradient using the conjugate gradient method to obtain the updated value. Q Model.

2. The full-wave equation based method according to claim 1 Q Chromatographic method, characterized in that, Step 2, acquiring earthquake-related data, includes: performing forward modeling using the Ricker wavelet frequency and formulas (1) and (2): （1） in, For pressure field data, Represents the particle velocity vector. Represents a memory variable. Indicates the bulk modulus of the medium. Indicates that the epicenter is located at Stress-strain relaxation parameters and Is with Q The relationship between them is as follows: （2） Where t is the seismic wave propagation time, This refers to the density of underground strata.

3. The full-wave equation based method according to claim 1 Q Chromatographic method, characterized in that, Step 3 includes: According to formula (3), the earthquake-related data is converted into a local domain and FFT is performed on the data therein to obtain the frequency of a single earthquake event; （3） in, S and s Represents earthquake records in the local domain and data domain. w Indicates a Gaussian time window. This represents the local frequency, containing the frequency of each seismic event. and These are wave number and seismic wave propagation time, respectively.

4. The full-wave equation based method according to claim 1 Q Chromatographic method, characterized in that, Step 4 includes: using formulas (4) and (5) to calculate the peak frequency shift of each seismic event in the local domain. : （4） （5） in, For the local spectrum of the observed data, To predict the local spectrum of the data, For frequency, The L2 norm of the product of the weight coefficients and the cross-correlation coefficients is given by the weight coefficients. The minimum represents It refers to peak frequency shift; C(f s ,t) This represents the cross-correlation between the spectrum of observed data and the spectrum of predicted data in the local domain. This represents the frequency shift obtained using equation (5).

5. The full-wave equation based method according to claim 1 Q Chromatographic method, characterized in that, Step 5 includes: using formula (6) to accompany the source: （6） in, As an accompanying source, X s Represents the location of the epicenter 、X r Represents the location of the detector 、t C represents time; C represents the cross-correlation between the spectrum of the observed data and the spectrum of the predicted data in the local domain. The frequency shift obtained using equation (5) represents the frequency shift. It represents a pressure field.

6. The full-wave equation based method according to claim 1 Q Chromatographic method, characterized in that, Step 6 includes: obtaining the reverse propagating wavefield by performing accurate wavefield extension calculations using formula (7). and : ； ； （7） in For the propagation time of seismic waves, and x represents the density and velocity of the subsurface strata, respectively. s and These represent the position coordinates of the shot point and the receiver point, respectively. Bulk modulus and For the attenuation parameters of underground strata, As an accompanying source, For the focal term, and for and The adjoint vector, for The adjoint vector, For pressure field data, Represents the particle velocity vector. This represents a memory variable.

7. The full-wave equation-based wave according to claim 1 Q Chromatographic method, characterized in that, Step 7 includes: The gradient is obtained using the reverse propagation wave field and formula (8); （8） in, Let be the objective function. and For the attenuation parameters of underground strata, For bulk modulus, For the propagation time of seismic waves, This represents the particle velocity vector during the forward modeling process. This represents the particle velocity vector during backward propagation. This is a memory variable used in the backpropagation process.

Images