A seismic velocity inversion method

Abstract

The present application relates to a kind of seismic velocity inversion method, belong to oil and gas geophysical engineering field.By the amplitude maximum of each shot estimated from the initial observed seismic data shot point position obtained, the combined time shift parameter of iteration process is determined by the maximum time shift parameter of combined stack and total iteration number, the combined single shot gather of iteration is generated by using combined time shift parameter to single shot delay stack, the combined seismic source is generated by using combined time shift parameter to seismic wavelet delay stack, the new simulated single shot is obtained by using combined seismic source initial velocity model and observation system forward, the proportion coefficient is determined by the amplitude of new simulated single shot and the amplitude of initial observation data, the amplitude of combined seismic source after forward is updated by proportion coefficient, the forward combined single shot is obtained by updating combined seismic source, initial velocity field and again forward, the final velocity field is obtained by the cross-correlation gradient of forward combined single shot and combined single shot gather, the energy difference of each shot is fully considered, and the stability and accuracy of inversion are improved.

Description

Technical Field

[0001] This invention relates to a seismic velocity inversion method, belonging to the field of oil and gas geophysical exploration engineering. Background Technology

[0002] In current exploration, seismic velocity field inversion is crucial for seismic data migration imaging and interpretation. Currently, the most accurate velocity modeling method is the full waveform inversion method, which comprehensively utilizes information such as amplitude, waveform, and frequency to finely characterize the subsurface velocity model. Full waveform inversion is an iterative velocity inversion process, and the amplitude and wavelet of the seismic data significantly influence the inversion results. However, due to factors such as acquisition costs, poor paths, noise, and topography, the energy differences between actual acquired shots are significant, leading to inaccurate extracted wavelets and a large discrepancy between the initial velocity model and the actual subsurface model. This instability in the inversion process poses a challenge to seismic data velocity inversion and is one of the key issues limiting the practical application of full waveform inversion.

[0003] The invention patent document with publication number CN114063160A discloses a seismic velocity inversion method and apparatus. The method includes: obtaining the combined time shift of the k-th iteration by combining the maximum time shift, the total number of iterations, and the number of combined superpositions; obtaining the combined coded supershot by combining the polarity-coded supershot of the k-th iteration and the combined time shift, and the weighting coefficients corresponding to each time shift; obtaining the combined source of the k-th iteration using the combined time shift of the k-th iteration, the weighting coefficients corresponding to each time shift, the seismic wavelet, and the polarity-coded matrix; obtaining the forward supershot Dcal,k(v) of the k-th iteration using the velocity field of the (k-1)-th iteration; obtaining the inversion gradient field of the k-th iteration using the forward supershot Dcal,k(v); obtaining the velocity field of the k-th iteration by combining the update step size of the k-th iteration, the inversion gradient field, and the velocity field of the (k-1)-th iteration; and obtaining the final velocity field through iterative iteration. This method mainly considers the acceleration effect, but does not take into account the impact of large energy differences between shots. When the energy difference between shots is large, the extracted seismic wavelets are not accurate, and the inversion process is unstable, which leads to low accuracy of the inversion velocity field of the obtained seismic data. Summary of the Invention

[0004] The purpose of this invention is to provide a seismic velocity inversion method to address the problem of improving the stability and accuracy of the velocity field inverted from seismic data when there is a large difference in shot energy between seismic data.

[0005] To achieve the above objectives, the present invention includes:

[0006] The present invention provides a seismic velocity inversion method, comprising the following steps:

[0007] 1) Acquire initial seismic observation data, initial velocity model, and observation system, and perform spectral analysis on the initial seismic observation data to obtain the dominant frequency of the initial seismic observation data;

[0008] 2) Denoise the initial seismic observation data and estimate the maximum seismic amplitude of each shot based on the shot point location traces of the denoised seismic observation data;

[0009] 3) Determine the combined time-shift parameters for each iteration based on the number of combined superpositions, the maximum time-shift parameter, and the total number of iterations, and use the combined time-shift parameters to perform time-delay superposition on single shots to generate the combined single-shot set for each iteration;

[0010] 4) Determine the dominant frequency of the seismic wavelet based on the dominant frequency of the initial seismic observation data, and use the combined time shift parameters to delay and superimpose the seismic wavelet to generate a combined source.

[0011] 5) Use the combined seismic source, initial velocity model and observation system to perform forward modeling to obtain a new simulated single shot. Determine the scaling factor based on the maximum amplitude of the simulated single shot obtained after forward modeling and the maximum amplitude obtained from the denoised observed seismic data in each iteration. Use the scaling factor, the maximum amplitude after forward modeling and the combined seismic source to obtain a new combined seismic source.

[0012] 6) Use the new combined seismic source and initial velocity field to perform forward modeling again to obtain the forward modeled combined shot for each iteration. Based on the difference and backpropagation of the forward modeled combined shot and the combined shot set in step 3), cross-correlate to obtain the inversion gradient field for the kth iteration. Calculate the optimal step size based on the inversion gradient and update the velocity field.

[0013] 7) Repeat steps 2) to 6) until the number of iterations or the velocity field is less than the preset value, then output the final velocity field.

[0014] Beneficial Effects: The seismic velocity inversion method of this invention first performs time-shifted combination superposition of seismic wavelets and initial observed seismic data to generate combined single-shot gathers for each iteration. The seismic wavelets are then delayed and superimposed using combined time-shift parameters to generate combined sources. Forward modeling is performed using the obtained combined sources, the acquired initial velocity model, and the observation system to obtain new simulated single shots. An amplitude scaling factor is determined based on the maximum amplitude of the new simulated single shot and the maximum amplitude obtained from the initial data. Using the scaling factor and the maximum amplitude after forward modeling, the combined sources are updated. Based on the updated combined sources, the initial velocity field is forward modeled again to obtain the final forward-modeled combined shot. Then, cross-correlation gradients are calculated. Compared with traditional multi-scale inversion methods, this method avoids the preprocessing step of low-pass filtering on the observed seismic data, thus enabling gradient calculation and realizing multi-scale inversion. This improves the accuracy of velocity field inversion from seismic waveforms, avoids the truncation effect of direct frequency-division filtering of seismic data, and enhances the stability and balance of the inversion. Furthermore, this invention fully considers the energy differences of each shot, improving the stability and balance of the inversion.

[0015] Furthermore, the calculation formula used for combining the time-shift parameters in step 3) is as follows:

[0016] τ i,k =2i*(K-10*(int(k / 10)+1)*τ Cmax,K / (C max -1) / K

[0017] In the formula, τ i,k Let τ be the time shift parameter for the k-th iteration and the i-th iteration. Cmax,K Let K be the maximum time shift, K be the total number of iterations, and C be the maximum time shift. max To set the number of combinations, the number of combinations is 3.

[0018] Furthermore, the calculation formula used for the new combined seismic source in step 5) is as follows:

[0019]

[0020] In the formula, For the combined seismic source of the nth shot in the kth iteration after the update, For the combined source of the nth shot in the kth iteration where the combined source has not been updated, λ k s is the amplitude proportionality coefficient for the kth order. n (t-τ i Let τ be the source wavelet of the nth shot source at time t, undergoing the i-th time shift. i Let Amp(ds) be the time shift for the i-th time iteration. n C represents the maximum amplitude of the path where the shot point is located. max To set the number of combinations and stacks.

[0021] Furthermore, the scaling factor corresponding to each iterative combination of seismic sources in step 5) is:

[0022] λ k =Amp(d cal,sn ,k) / Amp(ds n ) = abs(max(d cal,sn，k (t))) / Amp(ds n )

[0023] In the formula, λ k Let Amp(d) be the amplitude scaling factor for the kth oscillation. cal,sn ,k) represents the maximum amplitude value of the forward simulation of the nth shot point position in the kth iteration, Amp(ds) n ) represents the maximum amplitude of the path where the nth shot point is located, abs(max(d cal,sn，k (t))) represents the maximum value of the forward simulation amplitude at the position of the nth shot in the kth iteration, d cal,sn，k (t) is the forward simulation of the shot trajectory set for the position of the nth shot in the kth iteration.

[0024] Furthermore, the maximum amplitude of the path where the nth shot point is located is:

[0025]

[0026] In the formula, Let be the absolute value of the maximum value of the nth shot point from time 0 to T.

[0027] Furthermore, in step 3), the combined single-gun assembly... The calculation formula is:

[0028]

[0029] In the formula, d n (t-τ i () represents the distance τ from the nth shot at time t. i Shot record of time shift, τ i Let C be the time shift amount for the i-th time. max To set the number of combinations and stacks.

[0030] Furthermore, in step 6), the forward modeling combined single-shot... The calculation formula is:

[0031]

[0032] In the formula, L is the finite difference forward modeling operator, and v is the initial velocity field.

[0033] Furthermore, the time shift of the delay is less than twice the reciprocal of the dominant frequency of the seismic data. Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the method flow in an embodiment of the present invention;

[0035] Figure 2 This is a velocity field model diagram in an embodiment of the present invention;

[0036] Figure 3 This is a schematic diagram of the initial velocity field in an embodiment of the present invention;

[0037] Figure 4 This is a schematic diagram illustrating the single-source seismic data observed in this embodiment of the invention.

[0038] Figure 5 This is a diagram illustrating the combined source observation seismic data used in this embodiment of the invention.

[0039] Figure 6 This is a schematic diagram of a seismic wavelet in an embodiment of the present invention;

[0040] Figure 7 This is a schematic diagram of the combined seismic source in an embodiment of the present invention;

[0041] Figure 8 In this embodiment of the invention, the inverted velocity field is obtained without considering energy differences;

[0042] Figure 9 The inverted velocity field is obtained by the method of the present invention in this embodiment of the invention. Detailed Implementation

[0043] The present invention will now be described in further detail with reference to the accompanying drawings.

[0044] Example of seismic velocity inversion method:

[0045] like Figure 1 The method for seismic velocity inversion shown involves delaying and superimposing seismic wavelets using combined time-shift parameters to generate combined seismic sources. Forward modeling is then performed using the obtained combined seismic sources, the acquired initial velocity model, and the observation system to obtain new simulated single-shot data. An amplitude scaling factor is determined based on the maximum amplitude of the new simulated single-shot data and the maximum amplitude obtained from the initial data. Using the scaling factor and the maximum amplitude after forward modeling, the combined seismic sources are updated. Based on the updated combined seismic sources, the initial velocity field is forward modeled again to obtain the final forward-modeled combined shot data. Finally, cross-correlation gradients are calculated to update the velocity field, resulting in the final velocity field. The specific implementation is as follows:

[0046] 1) Acquire initial seismic observation data, initial velocity model, and observation system, and perform spectral analysis on the observation data to obtain the dominant frequency of the observation data;

[0047] The earthquake observation data are single-source earthquake data. The initial velocity model selected in this embodiment is as follows: Figure 2 The international standard Marmous model shown in this embodiment uses the Ricker wavelet for forward modeling of the initial velocity model, and obtains multiple single-shot records using finite-difference forward modeling. These multiple single-shot data are then used as observation data. To test the effectiveness of the method, this embodiment selects energy amplitude differences of 1-50, and selects shot points as shown... Figure 4 The data for a single gun is shown on the left side.

[0048] Spectral analysis was performed on the initial seismic observation data to obtain its dominant frequency. For the initial seismic observation data (such as...) Figure 4 The data shown is analyzed using Fourier transform to obtain the main frequency information. Let the main frequency be f, then the maximum time shift is calculated to be 1 / (2f). In this embodiment, the main frequency is 15 Hz, and the time interval is 0.001 s, so the maximum time shift is 16 calculation points.

[0049] 2) Estimate the maximum seismic amplitude of each shot based on the shot location.

[0050] To avoid the impact of noise in the initial seismic observation data on shot energy, this embodiment requires denoising the seismic observation data. The maximum seismic energy amplitude for each shot is estimated for the trace at the shot point location. The maximum amplitude for each shot trace is statistically determined as follows:

[0051]

[0052] In the formula, Amp(s) n ) represents the maximum amplitude value at the point of attack for the nth shot. It is the absolute value of the maximum value of the nth shot point from time 0 to T.

[0053] 3) Denoise the initial seismic observation data, determine the combined time-shift parameters for each iteration based on the number of combined stacks, the maximum time-shift parameter, and the total number of iterations, and use the combined time-shift parameters to perform time-delayed stacking of single shots to generate combined single-shot sets for each iteration.

[0054] The specific formula for combining the time shift parameters is as follows:

[0055] τ k,i =2i*(K-10*(int(k / 10)+1)*τ Cmax,K / (C max -1) / K

[0056] In the formula, τ i,k Let τ be the time shift parameter for the k-th iteration and the i-th iteration. Cmax,KLet K be the maximum time shift, K be the total number of iterations, and C be the maximum time shift. max To determine the number of stacking operations, in this embodiment, the maximum time shift is selected to be less than the reciprocal of twice the dominant frequency of the seismic data, i.e., 1 / (2f). The number of stacking operations can be set according to the accuracy requirements of the actual seismic velocity inversion results. To achieve better inversion results, the number of stacking operations, C, is set. max The value is chosen as 3. (int(k / 10)+1) is the integer division operation of the k-th iteration by 10. The purpose of (int(k / 10)+1) dividing by 10 is to maintain the stability of one operation in 10 iterations.

[0057] In this embodiment, combined shot sets are synthesized using initial seismic observation data, and the combined individual shots in the combined shot sets are updated using combined time shifts. The calculation formula is:

[0058]

[0059] In the formula, d n (t-τ i () represents the distance τ from the nth shot at time t. i Shot record of time shift, τ i Let τ be the time shift amount for the i-th time (where τ is the time shift amount in this embodiment). i That is, the combined time-shift parameters), C max To set the number of combinations and stacks. For the combined single-shot method of the nth shot in the kth iteration, such as... Figure 5 As shown, the phase axis of the combined single-shot gather becomes significantly thicker, and the dominant frequency of the seismic data is significantly reduced, thereby enhancing the stability of the inversion process.

[0060] 4) Determine the dominant frequency of the seismic wavelet based on the dominant frequency of the initial seismic observation data, and use the combined time-shift parameters to delay and superimpose the seismic wavelet to generate a combined seismic source.

[0061] The seismic wavelet is obtained using the dominant frequency of the initial seismic observation data. In this embodiment, the seismic wavelet used is the Ricker wavelet, and its dominant frequency is 15 Hz. Using the Ricker wavelet with a dominant frequency of 15 Hz, the combined seismic source is obtained by combining and superimposing the waves using the following formula.

[0062]

[0063] In the formula, For the combined seismic source synthesized in the kth generation iteration, s n (t-τ i ) indicates that the nth shot source wavelet passes through τ at time t. i The source of time shift, τ iThis represents the time shift of the i-th time (in this embodiment, the time shift is less than twice the reciprocal of the dominant frequency of the seismic data). In this embodiment, as... Figure 7 As shown, the dominant frequency of the combined seismic source obtained after calculation using this formula is... Figure 6 The frequency shown is significantly lower than that of the previous method, which can increase the stability of the inversion results.

[0064] 5) Using the combined seismic source, the initial velocity model and the new simulated single shot obtained by forward modeling of the observation system, the scaling factor is determined based on the maximum amplitude of the simulated single shot obtained after forward modeling in each iteration and the maximum amplitude obtained from the initial observed seismic data. Using the scaling factor, the maximum amplitude after forward modeling and the combined seismic source, a new combined seismic source is obtained.

[0065] The simulated single shot is generated using the combined source and initial velocity field, along with forward modeling updates from the observation system, and the amplitude scaling factor is calculated. Let the amplitude value of the nth shot combined source forward model be sn, then the amplitude scaling factor for the kth iteration is:

[0066] λ k =Amp(d cal,sn ,k) / Amp(ds n ) = abs(max(d cal,sn，k (t))) / Amp(ds n )

[0067] Where, λ k Let Amp(d) be the amplitude scaling factor for the kth oscillation. cal,sn ,k) represents the maximum amplitude value of the forward simulation of the nth shot point position in the kth iteration, Amp(ds) n ) represents the maximum amplitude of the path where the nth shot point is located, abs(max(d cal,sn，k (t))) represents the maximum value of the forward simulation amplitude at the position of the nth shot in the kth iteration, d cal,sn，k (t) is the forward simulation of the shot trajectory set for the position of the nth shot in the kth iteration.

[0068] The combined seismic source for the updated k-th iteration and n-th shot is:

[0069]

[0070] in, For the combined seismic source of the nth shot in the kth iteration after the update, For the combined source of the nth shot in the kth iteration where the combined source has not been updated, λ k s is the amplitude proportionality coefficient for the kth order. n (t-τ i Let τ be the source wavelet of the nth shot source at time t, undergoing the i-th time shift. iLet Amp(ds) be the time shift for the i-th time iteration. n ) represents the maximum amplitude of the path where the shot point is located.

[0071] Specifically, taking the amplitude value s1 of the first shot combined source forward modeling as an example, the amplitude scaling factor for the kth iteration is:

[0072] λ k =Amp(d cal,s1 ,k) / Amp(ds1)=abs(max(d cal,s1，k (t))) / Amp(ds1)

[0073] Where Amp(ds1) is the maximum amplitude of the trace where the first shot point is located, and abs(max(d cal,sn，k (t))) represents the maximum value of the forward simulation amplitude at the position of the first shot in the k-th iteration, d cal,sn，k (t) is the forward simulation of the shot trajectory set for the first shot in the k-th iteration.

[0074] Using the updated polarity coding to combine seismic sources The velocity field model v obtained in the (k-1)th iteration k-1 Using finite difference forward modeling, generate the forward modeling combination of the k-th and n-th shots (the new forward modeling shot). For example, when k=1, the velocity field model v k-1 For the initial velocity field model, such as Figure 3 As shown, the forward modeling combination single-shot (new forward modeling shot) is then performed. for:

[0075]

[0076] In the formula, L is the finite difference forward modeling operator, and v is the initial velocity field.

[0077] 6) Using the new combined source and initial velocity field, perform forward modeling again to obtain the forward modeled combined shot for each iteration. Based on the difference and backpropagation of the forward modeled combined shot and the combined shot set in step 3), cross-correlate to obtain the inversion gradient field for the kth iteration. Calculate the optimal step size based on the inversion gradient and update the velocity field.

[0078] The process of inverting the velocity field involves iteratively minimizing an objective function. This objective function is the combined shot set obtained from the forward modeling of a single shot and step 3). The specific formula for the difference objective function is as follows:

[0079]

[0080] Where E(v) is the L2 norm of the data residuals corresponding to the velocity field v. For the main simulation combination single cannon (new main simulation cannon) The summation of shots from the 1st to the Nsth shot, with a total number of shots of Ns.

[0081] The inversion gradient is obtained by performing a differential inverse propagation between the combined single-shot model and the combined shot set in step 3), and cross-correlating the results. The gradient is then normalized and superimposed to obtain the final gradient. The formula for calculating the inversion gradient is as follows:

[0082]

[0083] In the formula, g n,k Let δ be the inverted velocity field of the k-th iteration, and let δ be the derivative operator. Let L be the second-order partial derivative, and L be the finite-difference forward operand. The complex conjugate, or conjugate mechanism, is represented by <,>, where <,> represents the inner product operation. For the second derivative of the forward modeling data with respect to time t, Forward and backward propagation of data residuals This is the forward illumination term, and also the gradient normalization term.

[0084] The optimal step size is calculated using the inversion gradient to update the velocity field. Methods for calculating the step size include linear search, parabolic fitting, and parabolic interpolation, aiming to minimize the inversion gradient error corresponding to the velocity field in the k-th iteration. In this embodiment, parabolic interpolation is used to determine the update step size α for the k-th iteration. k Update the velocity field using the inverted gradient:

[0085] v k =v k-1 -α k g k

[0086] 7) Repeat steps 2) to 6) until the number of iterations or the velocity field is less than the preset value, then output the final velocity field.

[0087] When the error is less than the preset value, the obtained velocity field model is the final velocity field model. During the inversion process, as the maximum time shift gradually decreases, the dominant frequency of the combined coded supershot gradually increases, achieving the goal of multi-scale inversion and avoiding the truncation effect of direct frequency division filtering of seismic data. Figure 2 When performing the inversion of the Marmousi model shown, the following method is used: Figure 3 The initial velocity field model shown was iterated repeatedly. After the 100th iteration, the result was as follows: Figure 9 The inversion results shown are consistent with Figure 8 The comparison of inversion results not obtained using the method of this embodiment shows that... Figure 9It performs better inversion of background fields and complex tectonic layers, and can accurately invert velocity fields, eliminating the influence of instability during the inversion process.

Claims

1. A seismic velocity inversion method characterized by, The method comprises the following steps: 1) obtaining initial observed seismic data, an initial velocity model and an observation system, and performing spectrum analysis on the initial observed seismic data to obtain a main frequency of the initial observed seismic data; 2) denoising the initial observed seismic data, and estimating a maximum seismic amplitude of each shot according to a shot location trace of the denoised observed seismic data; 3) determining a combined time shift parameter of each iteration process according to a combined stacking number, a maximum time shift parameter and a total iteration number, and generating a combined single shot gather of each iteration by time delay stacking of a single shot using the combined time shift parameter; 4) determining a main frequency of a seismic wavelet according to the main frequency of the initial observed seismic data, and generating a combined seismic source by time delay stacking of the seismic wavelet using the combined time shift parameter; 5) performing forward calculation using the combined seismic source, the initial velocity model and the observation system to obtain a new simulated single shot, determining a proportional coefficient according to a maximum amplitude of the simulated single shot obtained after the forward calculation in each iteration process and the maximum amplitude obtained from the denoised observed seismic data, and obtaining a new combined seismic source using the proportional coefficient, the maximum amplitude after the forward calculation and the combined seismic source; 6) performing forward calculation again on the new combined seismic source, the initial velocity field and the observation system to obtain a forward calculation combined single shot of each iteration, performing reverse transmission by difference between the forward calculation combined single shot and the combined shot gather in the step 3), and obtaining an inversion gradient field of the kth iteration number by cross-correlation, calculating an optimal step length according to the inversion gradient, and updating the velocity field; 7) repeating the steps 2) to 6) until the iteration number is satisfied or the velocity field is less than a preset value, and outputting a final velocity field.

2. The seismic velocity inversion method of claim 1, wherein, The calculation formula of the combined time shift parameter in the step 3) is: τ i,k = 2i * (K - 10 * (int(k / 10) + 1) * τ Cmax,K / (C max - 1) / K In the formula, τ i,k is the time shift parameter of the i-th iteration of the k-th iteration, τ Cmax,K is the maximum time shift amount, K is the total number of iterations, C max is the set number of combination superposition.

3. The seismic velocity inversion method of claim 1, wherein, The calculation formula of the new combined seismic source in the step 5) is: In the formula, For the combined seismic source of the nth shot in the kth iteration after the update, For the combined source of the nth shot in the kth iteration where the combined source has not been updated, λ k s is the amplitude proportionality coefficient for the kth order. n (t-τ i Let τ be the source wavelet of the nth shot source at time t, undergoing the i-th time shift. i Let Amp(ds) be the time shift for the i-th time. n C represents the maximum amplitude of the path where the shot point is located. max To set the number of combinations and stacks.

4. The seismic velocity inversion method of claim 3, wherein, The proportional coefficient corresponding to the combined seismic source of each iteration in the step 5) is: λ k =Amp(d cal,sn ,k) / Amp(ds n )=abs(max(d cal,sn，k (t))) / Amp(ds n ) In the formula, λ k Let Amp(d) be the amplitude scaling factor for the kth oscillation. cal,sn ,k) represents the maximum amplitude value of the forward simulation of the nth shot point position in the kth iteration, Amp(ds) n ) represents the maximum amplitude of the path where the nth shot point is located, abs(max(d cal,sn，k (t))) represents the maximum value of the forward simulation amplitude at the position of the nth shot in the kth iteration, d cal,sn，k (t) is the forward simulation of the shot trajectory set for the position of the nth shot in the kth iteration.

5. The seismic velocity inversion method of claim 4, wherein, The maximum amplitude of the n-th shot location trace is: In the formula, is the absolute value of the maximum of the n-th shot point from time 0 to T.

6. The seismic velocity inversion method of claim 1, wherein, The step 3) combines the shot gathers The calculation formula is: In the formula, d n (t-τ i () represents the distance τ from the nth shot at time t. i Shot record of time shift, τ i Let C be the time shift amount for the i-th time. max To set the number of combinations and stacks.

7. The seismic velocity inversion method of claim 3, wherein, The forward combination single shot in the step 6) The calculation formula is: In the formula, L is a finite difference forward calculation operator, and v is the initial velocity field.

8. The seismic velocity inversion method of claim 3, wherein, The time shift amount of the time delay is less than 2 times the reciprocal of the main frequency of the seismic data.

Images