Earthquake velocity inversion method

Abstract

The invention relates to a seismic velocity inversion method, and belongs to the field of oil and gas geophysical prospecting engineering. Estimating the maximum amplitude value of each shot according to the shot point position of the obtained initial observation seismic data, determining the combined time shift parameter of the iteration process according to the maximum time shift parameter of the combined superposition number and the total number of iterations, and performing delay superposition on the single shots by using the combined time shift parameter to generate an iterative combined single shot gather; performing delay superposition on the seismic wavelets by using the combined time-shifting parameters to generate a combined seismic source, performing forward modeling by using a combined seismic source initial speed model and an observation system to obtain a new simulated single shot, determining a proportionality coefficient according to the amplitude of the new simulated single shot and the amplitude of initial observation data, updating the combined seismic source according to the amplitude obtained after forward modeling of the proportionality coefficient, and obtaining the seismic source according to the updated combined seismic source. And the combined seismic source is updated, the initial velocity field is subjected to forward modeling again to obtain a forward modeling combined single shot, the forward modeling combined single shot and the combined single shot gather are subjected to cross-correlation gradient calculation to obtain a final velocity field, the energy difference of each shot is fully considered, and the stability and accuracy of inversion are improved.

Description

Technical Field

[0001] The invention relates to a seismic velocity inversion method and belongs to the field of oil and gas geophysical exploration engineering. Background Art

[0002] In the current exploration field, the inversion of seismic velocity field is of great significance for the migration imaging and interpretation of seismic data. At present, the more accurate velocity modeling method is the full waveform inversion method, which uses information such as amplitude, waveform, and frequency to finely characterize the underground velocity model. Full waveform inversion is an iterative velocity inversion process. The amplitude and seismic wavelet of seismic data will have a relatively large impact on the inversion results of seismic data. However, due to the influence of many factors such as acquisition cost, bad tracks, noise, and terrain, the energy difference between the actual data shots is relatively large, the extracted seismic wavelets are not accurate, and the initial velocity model is quite different from the actual underground model. Therefore, the inversion process is unstable, which causes certain troubles for the velocity inversion of seismic data. It 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 method and device for seismic velocity inversion. The method includes: combining the maximum time shift, the total number of iterations and the number of combined superpositions to obtain the combined time shift of the kth iteration; combining the polarity coded supergun and the combined time shift of the kth iteration, and the weighting coefficient corresponding to each time shift to obtain the combined coded supergun, using the combined time shift of the kth iteration and the weighting coefficient corresponding to each time shift, the seismic wavelet and the polarity coding matrix to obtain the combined source of the kth iteration, and using and the velocity field of the k-1th iteration to obtain the forward supergun Dcal,k(v) of the kth iteration; using the forward supergun Dcal,k(v) and to obtain the inversion gradient field of the kth iteration; combining the update step size of the kth iteration and the inversion gradient field and the velocity field of the k-1th iteration to obtain the velocity field of the kth iteration, and iterating cyclically to obtain the final velocity field. This method mainly considers the acceleration effect, but does not consider the influence of the large energy difference between shots. When the energy difference between shots is large, the extracted seismic wavelet is not accurate and the inversion process is unstable, which leads to low accuracy of the inverted velocity field of the seismic data. Summary of the invention

[0004] The purpose of the present invention is to provide a seismic velocity inversion method to solve the problem of improving the stability and accuracy of the seismic data inversion velocity field when the shot-to-shot energy difference of the seismic data is large.

[0005] To achieve the above object, the solution of the present invention includes:

[0006] A seismic velocity inversion method of the present invention comprises the following steps:

[0007] 1) Obtaining initial observed seismic data, initial velocity model and observation system, and performing spectrum analysis on the initial observed seismic data to obtain the main frequency of the initial observed seismic data;

[0008] 2) De-noise the initial observed seismic data, and estimate the maximum seismic amplitude of each shot based on the shot point location of the de-noised observed seismic data;

[0009] 3) Determine the combined time shift parameter of each iteration process according to the number of combined superpositions, the maximum time shift parameter and the total number of iterations, and use the combined time shift parameter to delay superposition of the single shot to generate the combined single shot collection of each iteration;

[0010] 4) Determine the main frequency of the seismic wavelet according to the main frequency of the initial observed seismic data, and use the combined time-shift parameters to delay and superimpose the seismic wavelet to generate a combined source;

[0011] 5) Using the combined source, initial velocity model and observation system to perform forward modeling to obtain a new simulated single shot, determining a proportionality coefficient based on the maximum amplitude of the simulated single shot obtained after forward modeling and the maximum amplitude of the observed seismic data after denoising in each iteration, and using the proportionality coefficient, the maximum amplitude after forward modeling and the combined source to obtain a new combined source;

[0012] 6) forward modeling again using the new combined source and initial velocity field to obtain the forward modeled combined single shot of each iteration, performing differential backpropagation between the forward modeled combined single shot and the combined shot set in step 3), cross-correlating to obtain the inversion gradient field of the kth iteration, calculating the optimal step length based on the inversion gradient, and updating the velocity field;

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

[0014] Beneficial effects: The seismic velocity inversion method of the present invention first performs time-shift combined superposition on seismic wavelets and initial observed seismic data, generates a combined single-shot set of each iteration, uses combined time-shift parameters to delay and superpose the seismic wavelets to generate a combined source, uses the obtained combined source and the acquired initial velocity model and observation system to perform forward modeling to obtain a new simulated single shot, determines the amplitude proportional coefficient according to the new simulated single shot amplitude maximum value and the amplitude maximum value obtained from the initial data, updates the combined source by using the proportional coefficient and the forward modeled amplitude maximum value, forward models the initial velocity field again according to the updated combined source, obtains the final forward modeled combined shot, and then obtains the cross-correlation gradient. Compared with the traditional multi-scale inversion method, the pre-processing step of low-pass filtering the observed seismic data is avoided, and the gradient can be obtained, realizing the multi-scale inversion method, improving the accuracy of the seismic waveform inversion velocity field, avoiding the truncation effect of direct frequency division filtering of seismic data, and improving the accuracy of the seismic waveform inversion velocity field. At the same time, the present invention fully considers the energy difference of each shot and improves the stability and balance of the inversion.

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

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

[0017] 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, K is the total number of iterations, C max To set the number of combined superpositions, the number of combined superpositions is 3.

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

[0019]

[0020] In the formula, is the combined source of the nth shot in the kth iteration after update, is the combined source of the nth shot of the kth iteration without updating the combined source, λ k is the kth amplitude proportional coefficient, s n (t-τ i ) is the source of the nth shot source wavelet passing through the i-th time shift at time t, τ i is the time shift of the ith time, Amp(ds n ) is the maximum amplitude of the trace where the shot point is located.

[0021] Furthermore, the proportional coefficient corresponding to each iterative combination of earthquake 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 is the kth amplitude proportional coefficient, Amp(d cal,sn ,k) is the maximum amplitude value of the forward simulation of the position of the n-th shot point in the k-th iteration, Amp(ds n ) is the maximum amplitude of the trace where the nth shot is located, abs(max(d cal,sn，k (t))) is the maximum forward simulation amplitude of the position of the nth shot point in the kth iteration, d cal,sn，k (t) is the forward simulation shot trace set of the shot point position of the nth shot in the kth iteration.

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

[0025]

[0026] In the formula, It is the absolute value of the maximum value of the nth gun point from time 0 to T.

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

[0028]

[0029] Where, d n (t-τ i ) is the time at which the gun set passes through τ at the nth gun at time t i Shot record of time shift, τ i is the time shift of the ith time, C max To set the number of combinations to be superimposed.

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

[0031]

[0032] Where L is the finite difference forward operator and v is the initial velocity field.

[0033] Furthermore, the time shift of the delay is less than 2 times the inverse of the main frequency of the seismic data. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

[0037] Figure 4 Schematic diagram showing single-source observation seismic data used in an embodiment of the present invention;

[0038] Figure 5 It is a display diagram of combined source observation seismic data used in the embodiment of the present invention;

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

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

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

[0042] Fig. 9 The inverted velocity field is obtained by the method of the present invention in the embodiment of the present invention. DETAILED DESCRIPTION

[0043] The present invention will be further described in detail below in conjunction with the accompanying drawings.

[0044] Seismic velocity inversion method embodiment:

[0045] like Figure 1 A seismic velocity inversion method shown in the figure uses a combined time-shift parameter to delay and superimpose seismic wavelets to generate a combined source, uses the obtained combined source and the acquired initial velocity model and observation system to perform forward modeling to obtain a new simulated single shot, determines the amplitude proportionality coefficient based on the new simulated single shot amplitude maximum value and the amplitude maximum value obtained from the initial data, updates the combined source using the proportionality coefficient and the forward modeled maximum amplitude, and forward models the initial velocity field again based on the updated combined source to obtain the final forward modeled combined shot, then obtains the cross-correlation gradient, updates the velocity field, and obtains the final velocity field. The specific implementation is as follows:

[0046] 1) Obtain initial seismic observation data, initial velocity model and observation system, and perform spectrum analysis on the observation data to obtain the main frequency of the observation data;

[0047] The seismic observation data is single-source seismic data. The initial velocity model in this embodiment is selected as follows: Figure 2 The international standard Marmous model shown in the figure, for the initial velocity model, the seismic wavelet used for forward modeling in this embodiment adopts the Ricker wavelet, and uses finite difference forward modeling to obtain multiple single shot records, and uses multiple single shot data as observation data. In order to test the effectiveness of the method, this embodiment selects energy amplitude differences of 1-50, and selects the shot points as follows Figure 4 Single shot data for the left position shown.

[0048] Perform spectrum analysis on the initial observed seismic data to obtain the main frequency of the initial observed seismic data. Figure 4 The data shown in the figure is analyzed by Fourier transform to obtain the main frequency information of the data. If the main frequency is f, 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 point location.

[0050] In order to avoid the influence of the noise of the initial seismic observation data on the shot energy, the seismic observation data needs to be denoised in this embodiment. The maximum amplitude of the seismic energy of each shot is estimated for the shot point location channel, and the maximum amplitude of the channel where the shot point is located is determined by statistics for each shot:

[0051]

[0052] Where Amp(s n ) is the maximum amplitude value of the nth shot point, is the absolute value of the maximum value of the nth gun point from time 0 to T.

[0053] 3) De-noising the initial seismic observation data, determining the combined time-shift parameters of each iteration process according to the number of combined stacking, the maximum time-shift parameter and the total number of iterations, and using the combined time-shift parameters to delay the single shot and generate the combined single shot gun set of each iteration.

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

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

[0056] 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, K is the total number of iterations, C maxIn this embodiment, the maximum time shift is selected to be less than the inverse of twice the main frequency of the seismic data, that is, 1 / (2f). The number of combined stacking can be set according to the accuracy of the actual seismic velocity inversion results. In order to achieve a better inversion effect, the number of combined stacking C max The value is selected as 3, (int(k / 10)+1) is the rounding operation of dividing by 10 in the kth iteration, and the purpose of dividing (int(k / 10)+1) by 10 is to maintain the stability of one operation of 10 iterations.

[0057] In this embodiment, the initial observed seismic data is used to synthesize the combined shot collection, and the combined time shift is used to update the combined single shot in the combined shot collection, wherein the combined single shot gun collection The calculation formula is:

[0058]

[0059] Where, d n (t-τ i ) is the time at which the gun set passes through τ at the nth gun at time t i Shot record of time shift, τ i is the time shift of the ith time (wherein the time shift τ in this embodiment is i is the combined time shift parameter), C max To set the number of combinations to be superimposed. is the combined single shot of the nth shot in the kth iteration, such as Figure 5 As shown in the figure, the event axis of the combined single-shot gun set becomes significantly coarser and the main frequency of the seismic data is significantly reduced, thereby enhancing the stability of the inversion process.

[0060] 4) Determine the main frequency of the seismic wavelet according to the main frequency of the initial observed seismic data, and use the combined time-shift parameters to delay and superimpose the seismic wavelet to generate a combined earthquake source.

[0061] The main frequency of the initial seismic observation data is used to obtain the seismic wavelet, wherein the seismic wavelet used in this embodiment is the Ricker wavelet, and in this embodiment, the wavelet main frequency is 15 Hz. Using the Ricker wavelet with a main frequency of 15 Hz, the combined source is obtained by combining and superimposing according to the following formula:

[0062]

[0063] In the formula, is the kth generation iterative synthesis combined source, s n (t-τ i ) indicates that the source wavelet of the nth shot passes through τ at time t i The source of the time shift, τ iis the time shift of the ith time (the time shift of this embodiment is less than 2 times the inverse of the main frequency of the seismic data). Figure 7 As shown in the figure, the main frequency of the combined source calculated by this formula is Figure 6 The main frequency shown is significantly lower than that shown in Figure 1, which can increase the stability of the inversion results.

[0064] 5) Using the combined source, initial velocity model and observation system to obtain a new simulated single shot, determine the proportionality coefficient 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, and use the proportionality coefficient, the maximum amplitude after forward modeling and the combined source to obtain a new combined source.

[0065] The simulated single shot is updated by using the combined source and initial velocity field and the observation system, and the amplitude proportional coefficient is calculated. Assuming the amplitude value of the combined source forward modeling of the nth shot is sn, the amplitude proportional coefficient of 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] Among them, λ k is the kth amplitude proportional coefficient, Amp(d cal,sn ,k) is the maximum amplitude value of the forward simulation of the position of the n-th shot point in the k-th iteration, Amp(ds n ) is the maximum amplitude of the trace where the nth shot is located, abs(max(d cal,sn，k (t))) is the maximum forward simulation amplitude of the position of the nth shot point in the kth iteration, d cal,sn，k (t) is the forward simulation shot trace set of the shot point position of the nth shot in the kth iteration.

[0068] Then the updated combined source of the kth iteration and the nth shot is:

[0069]

[0070] in, is the combined source of the nth shot in the kth iteration after update, is the combined source of the nth shot of the kth iteration without updating the combined source, λ k is the kth amplitude proportional coefficient, s n (t-τ i ) is the source of the nth shot source wavelet passing through the i-th time shift at time t, τ iis the time shift of the ith time, Amp(ds n ) is the maximum amplitude of the trace where the shot point is located.

[0071] Specifically, taking the amplitude value s1 of the first shot combination source forward modeling as an example, the amplitude proportional coefficient of the kth iteration is:

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

[0073] Among them, Amp(ds1) is the maximum amplitude of the trace where the first shot point is located, Amp(ds1) is the maximum amplitude of the trace where the first shot point is located, abs(max(d cal,sn，k (t))) is the maximum forward simulation amplitude of the first shot position in the kth iteration, d cal,sn，k (t) is the forward simulation shot trace set of the first shot point position in the kth iteration.

[0074] Combined source with updated polarity coding The velocity field model v obtained in the k-1th iteration k-1 Using finite difference to generate the kth and nth forward combined single shots (new forward shots) For example, when k = 1, the velocity field model v k-1 is the initial velocity field model, such as Figure 3 As shown, the forward combined single gun (new forward gun) for:

[0075]

[0076] Where L is the finite difference forward operator and v is the initial velocity field.

[0077] 6) Using the new combined source and initial velocity field, forward modeling is performed again to obtain the forward modeling combined single shot of each iteration. Differential backpropagation is performed between the forward modeling combined single shot and the combined shot set in step 3), and the inversion gradient field of the kth iteration is obtained by cross-correlation. The optimal step size is calculated based on the inversion gradient, and the velocity field is updated.

[0078] The process of inverting the velocity field is to minimize an objective function through continuous iteration. The objective function here is the forward single shot and the combined shot set obtained in step 3). The specific formula of the objective function for the difference is:

[0079]

[0080] Among them, E(v) is the data residual second norm corresponding to the velocity field v, Forward combined single gun (new forward gun) It is the sum of the shots from the 1st to the Nsth shots, and the total number of shots is Ns.

[0081] The forward combined single shot and the combined shot set in step 3) are used for differential back propagation, and the inversion gradient is obtained by cross-correlation, wherein the gradient is normalized and superimposed to obtain the final gradient. The calculation formula of the inversion gradient is:

[0082]

[0083] In the formula, g n,k is the inverted velocity field of the kth iteration, δ is the derivative operator, is the second-order partial derivative, L is the finite difference forward operator, represents the complex conjugate, i.e. the conjugate device, <,> is the inner product operation, where is the second-order derivative of the forward modeling data with respect to time t, The forward and reverse propagation of the data residuals is It is the forward lighting term and also the gradient normalization term.

[0084] The optimal step length is calculated using the inversion gradient to update the velocity field. The step length can be calculated by using a linear search method, a parabolic fitting method, or a parabolic interpolation method, so that the k-th iteration inversion velocity field corresponds to the smallest inversion gradient error. In this embodiment, the parabolic interpolation method is used to obtain the update step length α of the k-th iteration. k , use the inverted gradient to update the velocity field:

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

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

[0087] When the error is less than the preset value, the velocity field model obtained is the final velocity field model. During the inversion process, the maximum time shift gradually decreases, and the main frequency of the combined coded super gun gradually increases, achieving the purpose of multi-scale inversion and avoiding the truncation effect of direct frequency division filtering of seismic data. Figure 2 When inverting the Marmousi model shown, use Figure 3 After repeated iterations, the initial velocity field model shown in Fig. 9 The inversion results shown are consistent with Figure 8 Compared with the inversion results without using the method of this embodiment, it can be seen that Fig. 9It is better for inverting background fields and complex structural layers, and can invert the velocity field more accurately, eliminating the influence of instability in the inversion process.

Claims

1. A seismic velocity inversion method, characterized in that: The following steps are involved: 1) Obtaining initial observed seismic data, initial velocity model and observation system, and performing spectrum analysis on the initial observed seismic data to obtain the main frequency of the initial observed seismic data; 2) De-noise the initial observed seismic data, and estimate the maximum seismic amplitude of each shot based on the shot point location of the de-noised observed seismic data; 3) Determine the combined time shift parameter of each iteration process according to the number of combined superpositions, the maximum time shift parameter and the total number of iterations, and use the combined time shift parameter to delay superposition of the single shot to generate the combined single shot collection of each iteration; 4) Determine the main frequency of the seismic wavelet according to the main frequency of the initial observed seismic data, and use the combined time-shift parameters to delay and superimpose the seismic wavelet to generate a combined source; 5) Using the combined source, initial velocity model and observation system to perform forward modeling to obtain a new simulated single shot, determining a proportionality coefficient based on the maximum amplitude of the simulated single shot obtained after forward modeling and the maximum amplitude of the observed seismic data after denoising in each iteration, and using the proportionality coefficient, the maximum amplitude after forward modeling and the combined source to obtain a new combined source; 6) Forward modeling the new combined source, initial velocity field and observation system again to obtain the forward modeled combined single shot of each iteration, perform differential backpropagation between the forward modeled combined single shot and the combined shot set in step 3), and obtain the inversion gradient field of the kth iteration by cross correlation, calculate the optimal step size according to the inversion gradient, and update the velocity field; 7) Repeat steps 2) to 6) until the number of iterations is met or the velocity field is less than a preset value, and then output the final velocity field.

2. The seismic velocity inversion method according to claim 1, characterized in that: The calculation formula used for the combined time shift parameter in 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, K is the total number of iterations, C max To set the number of combinations to be superimposed.

3. The seismic velocity inversion method according to claim 1, characterized in that: The calculation formula used for the new combined source in step 5) is: In the formula, is the combined source of the nth shot in the kth iteration after update, is the combined source of the nth shot of the kth iteration without updating the combined source, λ k is the kth amplitude proportional coefficient, s n (t-τ i ) is the source of the nth shot source wavelet passing through the i-th time shift at time t, τ i is the time shift of the ith time, Amp(ds n ) is the maximum amplitude of the trace where the shot point is located.

4. The seismic velocity inversion method according to claim 3, characterized in that: The proportional coefficient corresponding to each iterative combination of earthquake sources in 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 is the kth amplitude proportional coefficient, Amp(d cal,sn ,k) is the maximum amplitude value of the forward simulation of the position of the n-th shot point in the k-th iteration, Amp(ds n ) is the maximum amplitude of the trace where the nth shot is located, abs(max(d cal,sn，k (t))) is the maximum forward simulation amplitude of the position of the nth shot point in the kth iteration, d cal,sn，k (t) is the forward simulation shot trace set of the shot point position of the nth shot in the kth iteration.

5. The seismic velocity inversion method according to claim 4, characterized in that: The maximum amplitude of the trace where the nth shot point is located is: In the formula, It is the absolute value of the maximum value of the nth gun point from time 0 to T.

6. The seismic velocity inversion method according to claim 1, characterized in that: In step 3), the single-gun gun collection is combined The calculation formula is: Where, d n (t-τ i ) is the time at which the gun set passes through τ at the nth gun at time t i Shot record of time shift, τ i is the time shift of the ith time, C max To set the number of combinations to be superimposed.

7. The seismic velocity inversion method according to claim 3, characterized in that: The forward combined single shot in step 6) The calculation formula is: Where L is the finite difference forward operator and v is the initial velocity field.

8. The seismic velocity inversion method according to claim 3, characterized in that: The time shift of the delay is less than 2 times the reciprocal of the main frequency of the seismic data.

Images