A seismic velocity inversion method based on adaptive weighting

Abstract

The present invention relates to a seismic velocity inversion method based on adaptive weighting, comprising the following steps: 1) acquiring initial seismic data, which includes an initial velocity field, and taking the initial velocity field as the last iterative velocity field; 2) calculating a current iterative gradient field using the last iterative velocity field; 3) calculating a depth weighting coefficient of a particle according to a velocity of a particle in the last iterative velocity field and a maximum value of the velocity of the particle in the last iterative velocity field; 4) calculating an adaptive weighting coefficient of the particle according to the depth weighting coefficient of the particle and a set number of iterations; 5) multiplying the adaptive weighting coefficient of the particle by the gradient of the particle in the current iterative gradient field to obtain the current weighted gradient field; 6) using the current weighted gradient field to update the iteration step size for velocity update to obtain the current iterative velocity field, taking the current iterative velocity field as the last iterative velocity field, repeating steps 2) to 6), and taking the last iterative velocity field as the final inverted velocity field.

Description

Technical Field

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

[0002] In current exploration, seismic velocity field inversion plays a crucial role in migration imaging and reservoir prediction in seismic data processing. Full waveform inversion (FWI) fully utilizes the full wavefield information (travel time, amplitude, phase, etc.) in seismic data to accurately invert subsurface model parameters (velocity, density, etc.). However, FWI involves significant computational complexity, suffers from low inversion stability, and requires low-frequency information in the seismic data.

[0003] The Chinese invention patent document with the authorization announcement number CN108680957B discloses a weighted local cross-correlation time-frequency domain phase inversion method, which introduces the time-frequency domain phase information of seismic data into the cross-correlation objective function, and is called a weighted local cross-correlation time-frequency domain phase inversion method. The introduction of phase information alleviates the dependence of full waveform inversion on the initial velocity model. At the same time, adding a weight factor to the time-frequency domain objective function can greatly enhance the anti-noise ability and the stability of the inversion. In the weighted local cross-correlation time-frequency domain phase inversion method in the low-frequency band, a good initial velocity model can be obtained, which is then used in the weighted local cross-correlation time-frequency domain phase inversion method in the high-frequency band, and finally a high-resolution inversion result can be obtained.

[0004] However, due to the multi-solution nature of inversion, during the inversion process, due to the influence of shallow noise, the inversion of deep layers cannot achieve good results, which has a certain impact on the velocity inversion accuracy of seismic data. Summary of the Invention

[0005] The purpose of the present invention is to provide a seismic velocity inversion method based on adaptive weighting to solve the problem of low seismic data velocity inversion accuracy caused by shallow noise affecting deep inversion effect.

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

[0007] A seismic velocity inversion method based on adaptive weighting of the present invention comprises the following steps:

[0008] 1) Acquire initial seismic data, including an initial velocity field, which includes the velocity of each particle, and use the initial velocity field as the last iterative velocity field;

[0009] 2) Using the velocity field of the previous iteration, the gradient field of the current iteration is calculated. The gradient field includes the gradient of each particle.

[0010] 3) Calculate the depth weighting coefficient corresponding to each particle based on the velocity of each particle in the last iterative velocity field and the maximum value of the particle velocity in the last iterative velocity field;

[0011] 4) Calculate the adaptive weighting coefficient corresponding to each particle based on the depth weighting coefficient corresponding to each particle, the number of iterations at that time, and the set total number of iterations;

[0012] 5) Multiply the adaptive weighting coefficient corresponding to each particle by the gradient of the particle in the current iterative gradient field to calculate the current weighted gradient field;

[0013] 6) Use the current weighted gradient field to update the iteration step size and perform velocity update to obtain the current iteration velocity field. Use the current iteration velocity field as the previous iteration velocity field. Repeat steps 2) to 6) based on the total number of iterations, and use the final iteration velocity field as the final seismic velocity field obtained by inversion.

[0014] The present invention obtains an initial velocity field from initial seismic data, uses the initial velocity field as the first iterative velocity field, calculates a second iterative gradient field based on the first iterative velocity field, and calculates a depth weighting coefficient for each particle based on the maximum value of the particle velocity in the first iterative velocity field and the velocity of each particle in the first iterative velocity field. In addition to changing with depth, the iterative velocity field also changes with the number of iterations. The depth weighting coefficient is combined with the number of iterations to calculate an adaptive weighting coefficient for each particle. The adaptive weighting coefficient is multiplied by the second iterative gradient field to obtain a second weighted gradient field. The second weighted gradient field is used to calculate the iteration step size and perform velocity update to obtain a second iterative velocity field. The second iterative velocity field is used as the third iterative velocity field for repeated iterative calculations, and the last iterative velocity field is used as the final velocity field obtained by inversion.

[0015] The velocity field finally inverted by the present invention can better invert deep complex structures, and the velocity field inversion of depressions is more accurate.

[0016] Furthermore, in step 3), the depth weighting coefficient corresponding to each particle is calculated using the following formula:

[0017]

[0018] Among them, Cm max，k-1 is the maximum value of the particle velocity in the k-1th iteration velocity field, α k (x, z) is the depth weighting coefficient of the particle (x, z), m k-1(x, z) is the velocity of the particle (x, z) in the k-1th iteration velocity field, x represents the lateral coordinate, and z represents the depth coordinate.

[0019] Furthermore, in step 4), the adaptive weighting coefficient corresponding to each particle is calculated using the following formula:

[0020]

[0021] Among them, N iter is the number of iterations set, k is the kth iteration, α k (x, z) is the depth weighting coefficient of the particle (x, z), β k (x, z) is the adaptive weighting coefficient of the particle (x, z), x represents the lateral coordinate, and z represents the depth coordinate.

[0022] Furthermore, in step 6), the iteration step size is obtained using a parabolic fitting method or a linear search method.

[0023] Furthermore, in step 6), the speed is updated using the following formula:

[0024]

[0025] Among them, v k represents the k-th iteration velocity field, represents the kth weighted gradient field, v k-1 is the k-1th iteration velocity field, β k Represents the updated iteration step size.

[0026] Furthermore, the initial earthquake data also includes source wavelets, an observation system, and observed earthquake data. The method used in step 2) to calculate the gradient field of the current iteration using the velocity field of the previous iteration includes the following steps:

[0027] ①Use the seismic wavelet and the last iterative velocity field to perform forward simulation to obtain the forward wave field and forward seismic data at each moment;

[0028] ② According to the target functional of full waveform inversion, the difference between the observed seismic data and the forward seismic data is obtained;

[0029] ③ The difference obtained in step ② is back-propagated and cross-correlated with the forward wave field obtained in step ① to obtain the gradient field of the current iteration.

[0030] Furthermore, the initial seismic data is filtered, and the filtering method adopted is frequency division filtering or Wiener filtering.

[0031] Furthermore, in step ①, the velocity field and seismic wavelet of the last iteration are subjected to finite difference forward modeling using the following formula to obtain forward seismic data:

[0032]

[0033] Among them, m represents the particle velocity, f represents the source term, u represents the seismic wave field, x represents the lateral coordinate, z represents the depth coordinate, and t represents time.

[0034] Furthermore, in step ③, the gradient of each particle is determined using the steepest descent method or the conjugate gradient method, thereby obtaining the gradient field of the current iteration.

[0035] Furthermore, in step ②, the target functional used is a two-norm target functional. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 1 is a flow chart of the seismic velocity inversion method based on adaptive weighting of the present invention;

[0037] Figure 2 This is a single shot record chart used in the embodiment of the method of the present invention;

[0038] Figure 3 This is the standard velocity field diagram given by the Marmousi model in the present invention;

[0039] Figure 4 Schematic diagram of the initial velocity field of the seismic velocity inversion method of the present invention;

[0040] Figure 5 is a gradient field map updated by the method of the present invention;

[0041] Figure 6 is a gradient field map updated using prior art methods;

[0042] Figure 7 It is a seismic velocity field map obtained by using the inversion method of the present invention;

[0043] Figure 8 It is a velocity field map obtained by conventional inversion using the standard velocity field. DETAILED DESCRIPTION

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

[0045] Example:

[0046] The present invention provides a seismic velocity inversion method. First, a forward modeling simulation is performed using input shot records, seismic wavelets and an initial velocity field to obtain a forward wave field and forward seismic data at each moment. According to a target functional of full waveform inversion, a difference between the observed seismic data and the forward seismic data is obtained. Wave field backpropagation is performed using the seismic data difference, and cross-correlation is performed with the seismic forward wave field to obtain a gradient field of full waveform inversion. Then, a depth weighting coefficient matrix is obtained from the seismic velocity field, and an adaptive weighting coefficient is obtained according to the number of iterations to obtain an adaptive weighting matrix. Finally, the gradient field is processed using the adaptive weighting matrix to obtain a new weighted gradient field. The weighted gradient field is used to update the velocity, thereby achieving effective updating of the deep velocity field and improving the accuracy of the seismic waveform inversion velocity field.

[0047] like Figure 1 The flowchart of the seismic velocity inversion method of the present invention is shown, which specifically includes the following steps:

[0048] 1) Obtain observed seismic data, initial velocity field and source wavelet.

[0049] The observation data in actual oil and gas seismic exploration is two-dimensional seismic data, such as Figure 2 The figure shows a schematic diagram of two-dimensional seismic data of a single source provided by an embodiment of the present invention. The observation system used in the present invention is a full-receiver observation system. The initial velocity field is given by a gradient model. Figure 3 The Marmousi model shown in the figure is used to test the effectiveness of the test method. The initial velocity field is given as Figure 4 As shown in the figure, the source wavelet utilizes the Ricker wavelet commonly used in seismic exploration. The Ricker wavelet selected in this invention has a dominant frequency of 15 Hz and an amplitude of 1.0. Wiener filtering is performed on the acquired seismic data using Ricker wavelets with dominant frequencies ranging from 3 Hz to 15 Hz as the target wavelet to obtain more accurate seismic data. In addition to this embodiment, frequency division filtering can also be used for filtering.

[0050] 2) Use seismic wavelets and initial velocity fields to perform forward simulation to obtain the forward wave field and forward data at each moment.

[0051] First use Figure 4 The initial velocity field and the given seismic wavelet are used for finite difference forward modeling to obtain the seismic wave field at each moment. The seismic wave field is stored and then the forward modeling data of the earthquake is obtained using a full-receiver observation system. The forward wave equation used in the present invention is:

[0052]

[0053] Among them, m represents the particle velocity, f represents the source term, u represents the seismic wave field, x represents the lateral coordinate, z represents the depth coordinate, and t represents time.

[0054] 3) According to the target functional of full waveform inversion, the difference between the observed seismic data and the forward modeled seismic data is obtained.

[0055] The process of inverting the velocity field is to use the target functional, which is continuously minimized through iteration, to obtain the difference between the observed seismic data in step 1) and the forward seismic data in step 2). Specifically, a target functional is first given, and then the difference is obtained based on the target functional. This embodiment uses the two-norm target functional commonly used in inversion, expressed as:

[0056]

[0057] Among them, u obs and u cal They represent the field observation earthquake record and forward simulation record at time t, respectively, and the earthquake source point is x s , the receiving point is x r , m is the velocity parameter model, E(m) is the data residual two-norm corresponding to the velocity field m, which is the target functional used in the present invention.

[0058] 4) The seismic data difference obtained in step 3) is subjected to wavefield backpropagation and cross-correlated with the seismic forward wavefield in step 3) to obtain the gradient field of full waveform inversion (i.e., the gradient field of the current iteration).

[0059] In this embodiment, the gradient is updated using the steepest descent method based on the seismic data difference. The data difference is back-propagated to obtain the Kth updated gradient value g k , the gradient update formula is as follows:

[0060]

[0061] Where mk represents the inverted velocity field of the Kth iteration, and δ is the derivative operator.

[0062] In addition to this embodiment, as another embodiment, the conjugate gradient method can also be used to determine the gradient.

[0063] 5) Use the input velocity field (i.e. the velocity field of the last iteration) to obtain the depth weighting coefficient matrix.

[0064] The existing technology uses the gradient formula in step 4) to update the velocity. The key of the present invention lies in the processing of the gradient. It is necessary to obtain a depth weighting coefficient matrix to achieve better inversion of the deep part. The present invention uses the law of velocity field change with depth to construct the depth weighting coefficient, and it can adaptively change with the number of iterations k:

[0065]

[0066] Among them, Cm max，k-1 is the maximum value of the particle velocity in the k-1th iteration velocity field, α k (x, z) is the depth weighting coefficient of the particle (x, z), m k-1 (x, z) represents the particle velocity of the velocity field obtained by k-1 iterations, x represents the lateral coordinate, and z represents the depth coordinate.

[0067] It should be noted that α k (x, z) is the depth weighting coefficient for the particle (x, z). Once the depth weighting coefficient corresponding to each particle is obtained, the depth weighting coefficient matrix can be obtained. In addition, if this is the first iteration, the input velocity field is the initial velocity field. If this is not the first iteration, the input velocity field is the velocity field obtained in the previous iteration.

[0068] 6) Obtain an adaptive weighting coefficient matrix through the number of iterations, and apply it to the depth weighting coefficient matrix in step 5) to obtain an adaptive weighting matrix.

[0069] The process of iteratively obtaining the velocity field varies not only with depth but also with the number of iterations, requiring the shallow layers to be updated first before the deep layers. The formula used to obtain the adaptive weighting coefficient using the depth weighting coefficient matrix in step 5) is:

[0070]

[0071] Among them, N iter is the total number of iterations, k is the kth iteration, α k (x, z) is the depth weighting coefficient of the particle (x, z) in step 5), β k (x, z) is the adaptive weighting coefficient corresponding to the particle (x, z). From the formula, it can be seen that the depth weighting coefficient increases with the increase of the number of iterations.

[0072] It should be noted that β k (x, z) is the adaptive weighting coefficient corresponding to the mass point (x, z). Once the adaptive weighting coefficient corresponding to each mass point is obtained, the adaptive weighting coefficient matrix can be obtained.

[0073] 7) Apply the adaptive weighting matrix obtained in step 6) to the gradient field in step 4) to obtain a new weighted gradient field (i.e., the current weighted gradient field). The specific formula is:

[0074]

[0075] Among them, β k (x, z) is the adaptive weighting matrix obtained in step 6), g k is the k-th gradient field obtained in step 4), Represents the k-th weighted gradient field.

[0076] like Figure 5 The gradient field updated by the present invention is as follows Figure 6 The gradient field of the present invention is updated from the prior art. The energy of the gradient field of the present invention is more prominent in the deep part than the traditional method, and the gradient field is more balanced.

[0077] 8) Use the weighted gradient field from step 7) to find the appropriate iterative step size and perform velocity update.

[0078] Use the parabolic interpolation method to obtain the iterative step size γ for the kth iteration update k , and then use the weighted gradient field in step 7) to update the velocity field to obtain the velocity field of the current iteration:

[0079]

[0080] Among them, v k represents the inverted velocity field of the kth iteration, is the kth weighted gradient field in step 7), v k-1 is the inverted velocity field of the k-1th iteration. If it is the first iteration, k = 1, then v k-1 is the initial velocity field.

[0081] 9) Output the final seismic velocity field after satisfying the total number of iterations.

[0082] Repeat steps 2) to 8) until the set number of iterations is met. When the last iteration is completed, the target velocity field is output.

[0083] After the first iteration, a velocity field is updated and the updated velocity field is used as the next initial velocity field to calculate the gradient field. The velocity field is iterated based on the gradient field. According to the initial velocity field and the main frequency from 3Hz to 15Hz, each main frequency is iterated 10 times, and the total number of iterations is 120 times. Until the last iteration is completed, the target velocity field is output.

[0084] The present invention uses Wiener filtering to perform filtering inversion from 3Hz to 15Hz, completes full waveform inversion and iterates, and uses the target velocity field obtained by the present invention to perform further inversion to obtain the inversion results as shown below: Figure 7 As shown in Figure 2, it can be seen that the deep complex structure is well inverted. Using the existing technology of gradient update, such as Figure 8 As shown, it can be seen that the velocity field of deep structure inversion is not accurate, and the standard Figure 3 In contrast, the velocity field inversion of the present invention is more accurate in depressions.

Claims

1. A seismic velocity inversion method based on adaptive weighting, characterized in that: The steps include: 1) Acquire initial seismic data, including an initial velocity field, which includes the velocity of each particle, and use the initial velocity field as the last iterative velocity field; 2) Using the velocity field of the previous iteration, the gradient field of the current iteration is calculated. The gradient field includes the gradient of each particle. 3) Calculate the depth weighting coefficient corresponding to each particle based on the velocity of each particle in the last iterative velocity field and the maximum value of the particle velocity in the last iterative velocity field; 4) Calculate the adaptive weighting coefficient corresponding to each particle based on the depth weighting coefficient corresponding to each particle, the number of iterations at that time, and the set total number of iterations; 5) Multiply the adaptive weighting coefficient corresponding to each particle by the gradient of the particle in the current iterative gradient field to calculate the current weighted gradient field; 6) Use the current weighted gradient field to update the iteration step size and perform velocity update to obtain the current iteration velocity field. Use the current iteration velocity field as the previous iteration velocity field. Repeat steps 2) to 6) based on the total number of iterations, and use the final iteration velocity field as the final seismic velocity field obtained by inversion.

2. The adaptive weighted seismic velocity inversion method according to claim 1, characterized in that: In step 3), the depth weighting coefficient corresponding to each particle is calculated using the following formula: Among them, Cm max，k-1 is the maximum value of the particle velocity in the k-1th iteration velocity field, α k (x, z) is the depth weighting coefficient of the particle (x, z), m k-1 (x, z) is the velocity of the particle (x, z) in the k-1th iteration velocity field, x represents the lateral coordinate, and z represents the depth coordinate.

3. The adaptive weighted seismic velocity inversion method according to claim 1, characterized in that: In step 4), the adaptive weighting coefficient corresponding to each particle is calculated using the following formula: Among them, N iter is the total number of iterations set, k is the kth iteration, α k (x, z) is the depth weighting coefficient of the particle (x, z), β k (x, z) is the adaptive weighting coefficient of the particle (x, z), x represents the lateral coordinate, and z represents the depth coordinate.

4. The adaptive weighted seismic velocity inversion method according to claim 1, characterized in that: In step 6), the iterative step size is obtained using a parabolic fitting method or a linear search method.

5. The adaptive weighted seismic velocity inversion method according to claim 1, characterized in that: In step 6), the speed is updated using the following formula: Among them, v k represents the k-th iteration velocity field, represents the kth weighted gradient field, v k-1 is the k-1th iteration velocity field, γ k Represents the updated iteration step size.

6. The adaptive weighted seismic velocity inversion method according to claim 1, characterized in that: The initial seismic data also includes seismic wavelets, an observation system, and observed seismic data. The method used in step 2) to calculate the gradient field of the current iteration using the velocity field of the previous iteration includes the following steps: ①Use the seismic wavelet and the last iterative velocity field to perform forward simulation to obtain the forward wave field and forward seismic data at each moment; ② According to the target functional of full waveform inversion, the difference between the observed seismic data and the forward seismic data is obtained; ③ The difference obtained in step ② is back-propagated and cross-correlated with the forward wave field obtained in step ① to obtain the gradient field of the current iteration.

7. The method for seismic velocity inversion based on adaptive weighting according to claim 1 or 6, characterized in that: The initial seismic data is also filtered, and the filtering method used is frequency division filtering or Wiener filtering.

8. The adaptive weighted seismic velocity inversion method according to claim 6, characterized in that: In step ①, the finite difference forward modeling is performed on the velocity field and seismic wavelet of the last iteration using the following formula to obtain the forward seismic data: Among them, m represents the particle velocity, f represents the source term, u represents the seismic wave field, x represents the lateral coordinate, z represents the depth coordinate, and t represents time.

9. The adaptive weighted seismic velocity inversion method according to claim 6, characterized in that: In step ③, the steepest descent method or conjugate gradient method is used to determine the gradient of each particle, and then the gradient field of the current iteration is obtained.

10. The method for seismic velocity inversion based on adaptive weighting according to claim 6, characterized in that: In step ②, the target functional used is the two-norm target functional.

Images