A 3D seismic velocity inversion method based on sparse constraints

Abstract

The present invention relates to a three-dimensional seismic velocity inversion method based on sparse constraints. The present invention establishes an objective function for sparse constrained full waveform inversion based on polarity coding, obtains the difference between coded observation seismic data and coded forward seismic data, cross-correlates the back-propagation wave field with the coded forward seismic wave field, obtains a gradient field for full waveform inversion, obtains a suitable iteration step size, performs velocity update, then uses fast dictionary learning to learn the velocity field to obtain a suitable dictionary set, performs sparse constraint denoising on the seismic data to obtain a sparse constrained velocity field, weights the sparse constrained velocity field and the original velocity field to obtain a new velocity field, and outputs the final seismic velocity field after satisfying the number of iterations. Because the synthesized polarity coding observation records are far less than the observation data, the amount of calculation is reduced. The synthesized sparse constrained inversion velocity field has better noise resistance for the seismic data, thereby improving the accuracy of the seismic waveform inversion velocity field.

Description

Technical Field

[0001] The present invention belongs to the field of oil and gas geophysical exploration engineering, and in particular relates to a three-dimensional seismic velocity inversion method based on sparse constraints. Background Art

[0002] In current exploration, the inversion of three-dimensional seismic velocity fields is crucial for migration imaging and subsequent geological interpretation in seismic data processing. Full waveform inversion accurately inverts the subsurface seismic velocity field by utilizing information such as travel time, phase, and amplitude in seismic data. However, full waveform inversion of seismic velocity fields presents numerous challenges, including susceptibility to noise interference, high computational complexity, and inversion stability. Due to factors such as acquisition costs, noise, topography, and bad tracks, actual seismic data often contain noise and missing earthquake data, which poses challenges for velocity inversion. Furthermore, the computational complexity of 3D seismic data inversion is enormous. Therefore, developing optimization methods for full waveform inversion is crucial.

[0003] The Chinese invention patent with the authorization number CN105353405B discloses a full waveform inversion method and system. The method includes the following steps: inversion based on onshore seismic data to obtain the background velocity field of the seismic profile; obtaining the well velocity vector in the seismic profile based on the logging data of known wells; using the well velocity vector to interpolate the background velocity field to obtain an initial velocity model coupled with the logging data in the low-frequency range; forward calculation based on the initial velocity model to obtain a velocity perturbation model, updating the initial velocity model according to the perturbation model, and determining the inversion model in the absence of low-frequency information in the background velocity field. The present invention uses the logging data of the wells to constrain the initial velocity model, and uses the rich high-frequency information and complete low-frequency components of the logging data to supplement the limited bandwidth of the seismic data, and obtains an accurate inversion final model in the subsequent iterative calculation. It mainly uses the logging constraint to enhance the low-frequency component.

[0004] Chinese invention patent application publication number CN111505714A discloses a rock physics-constrained direct envelope inversion method for elastic waves. This method uses wavefield mode decomposition and direct envelope inversion to obtain the P-wave velocity structure of strongly scattering media, and then calculates the S-wave velocity structure of the medium based on rock physics relationships. First, the elastic wavefield undergoes wavefield mode decomposition to obtain the P-wave field. The P-wave velocity gradient is then obtained by cross-correlating the forward P-wave and accompanying P-wave fields. Second, an update to the S-wave velocity is calculated based on rock physics relationships, which yields the large-scale P- and S-wave velocity structures of the strongly scattering medium. Finally, full waveform inversion is performed using the direct envelope inversion results as the initial model to obtain a high-precision P- and S-wave velocity model for the strongly scattering medium. This method improves the decoupling of P- and S-wave velocities in strongly scattering media, resulting in a high-quality P-wave velocity structure. By applying rock physics constraints during the inversion, the S-wave velocity structure of the strongly scattering medium is obtained. This method primarily utilizes rock physics constraints, requires knowledge of rock physics properties, and primarily operates on S-waves.

[0005] Therefore, the use of coded super guns can greatly improve computational efficiency, but it will also bring about the problem of crosstalk noise. At the same time, the noise of field observation data will also affect the inversion results. It is very necessary to develop a fast seismic velocity inversion method suitable for seismic data containing noise. Summary of the Invention

[0006] The purpose of the present invention is to provide a three-dimensional seismic velocity inversion method based on sparse constraints, which has good adaptability to seismic data containing more noise and improves the stability and computational efficiency of the seismic data inversion velocity field.

[0007] In order to achieve the above objectives, the technical solution adopted by the present invention is:

[0008] A 3D seismic velocity inversion method based on sparse constraints includes the following steps:

[0009] S1, obtain seismic observation data, initial velocity field, source wavelet, polar coding matrix and observation system;

[0010] S2, synthesizes super-coded observation data using seismic observation data, polar coding matrix and observation system;

[0011] S3, synthesizing a coded source wavelet using the source wavelet and the polarity coding matrix, and performing forward modeling using the coded seismic wavelet and the initial velocity field to obtain the coded seismic forward wavefield and coded forward seismic data at each moment;

[0012] S4, establishing the objective function of full waveform inversion based on polarity-coded sparse constrained T distribution, and obtaining the difference between synthetic super-coded observation data and coded forward seismic data;

[0013] S5, performs wavefield backpropagation on the seismic data difference obtained in S4, and cross-correlates it with the coded seismic forward wavefield in S3 to obtain the gradient field of full waveform inversion;

[0014] S6, using the gradient field of the full waveform inversion in S5, calculates the iteration step size and performs velocity update;

[0015] S7, performing dictionary learning on the velocity field to obtain a dictionary set;

[0016] S8, using the dictionary set to perform sparse constraint denoising on the seismic data, obtain a sparse constraint velocity field, and perform velocity update;

[0017] S9, repeat S2 to S8, and output the final seismic velocity field when the number of iterations meets the requirements.

[0018] The three-dimensional seismic velocity inversion method of the present invention reduces the amount of calculation because the synthesized polarity-coded observation records are far less than the observation data. The inversion velocity field with synthesized sparse constraints has better noise resistance for seismic data and improves the accuracy of the seismic waveform inversion velocity field.

[0019] Preferably, in S2, the super-encoded observation data is synthesized using formula (1):

[0020]

[0021] In formula (1), d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the synthetic super-encoded observation data for k iterations.

[0022] Preferably, in S3, the coded source wavelet is synthesized using formula (2):

[0023]

[0024] In formula (2), s n (x s ,y s ,t) is the shot point position (x s ,y s )’s earthquake source wavelet at the tth moment of the nth shot, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the coded source wavelet of k iterations, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=ys ,δ(x|x s ,y|y s )=1, otherwise it is 0.

[0025] Further preferably, in S3, the forward wave equation applied is:

[0026]

[0027] In formula (3), v is the medium velocity of the model, f is the source term of the earthquake, which is constructed using formula (2), u represents the seismic wave field, x, y represent the lateral coordinates and z represents the depth coordinate, and t represents time.

[0028] Preferably, in S4, the objective function is constructed according to formula (4):

[0029]

[0030] Among them, s and ε are the degree of freedom parameter and scale factor respectively; Represents the forward simulation of the velocity model v(x,y,z), where is the source term used, s n (x s ,y s ,t) is the shot point position (x s ,y s ) of the earthquake source wavelet at the tth moment of the nth shot, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=y s ,δ(x|x s ,y|y s )=1, otherwise 0; represents the synthetic coded observation data item, d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the encoding polarity of the nth shot in the kth iteration, Ns represents the total number of observations; λ is the regularization parameter, S represents the sparse constraint on v(x, y, z), and the fast dictionary learning denoising constraint is used, E(v) is the objective function of the inversion, Ne is the total number of encoded super-observations, and j is the sequence number of the encoded super-observation data.

[0031] Preferably, in S5, the gradient is updated according to formula (5):

[0032]

[0033] Among them, v krepresents the inverted velocity field of the kth iteration, and δ is the derivative operator.

[0034] Preferably, in S6, the current background velocity field is iterated by calculating the optimal iteration step size.

[0035] Further preferably, in S6, the speed is updated according to formula (6):

[0036] v k =v k-1 -α k g k (6)

[0037] Among them, v k represents the inverted velocity field of the kth iteration, g k is the kth weighted gradient field, 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.

[0038] Preferably, in S7, dictionary learning is performed according to formula (7):

[0039]

[0040] Among them, G is the dictionary set to be trained, v k is the velocity field obtained for the kth iteration, S is the sparse representation term, is the sparse coefficient corresponding to the k-th velocity field, and λ is the regularization parameter.

[0041] Preferably, in S8, the speed is updated using formula (8):

[0042]

[0043] in, is the new updated velocity field, S is the sparse representation term, is the sparse coefficient corresponding to the k-th velocity field, Indicates that only , ε is the retention coefficient threshold. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 Schematic diagram of a flow chart of an embodiment of the present invention;

[0045] Figure 2 Schematic diagram of three-dimensional seismic data of a single source in an embodiment of the present invention;

[0046] Figure 3 This is a standard depression model diagram in an embodiment of the present invention;

[0047] Figure 4Schematic diagram of the initial velocity field according to an embodiment of the present invention;

[0048] Figure 5 is the polar coding matrix in the embodiment of the present invention;

[0049] Figure 6 is the target velocity field obtained by using the method of the embodiment of the present invention;

[0050] Figure 7 To obtain the velocity field using the conventional full waveform inversion method. DETAILED DESCRIPTION

[0051] The present invention first synthesizes super-encoded observation data using input seismic observation data, a polar coding matrix, and an observation system. Then, a coded source wavelet is synthesized using a seismic wavelet and the polar coding matrix. A forward simulation is performed on the coded seismic wavelet and the initial velocity field to obtain coded forward seismic data. According to the objective function of sparse constrained full waveform inversion based on polar coding, the difference between the coded observation seismic data and the coded forward seismic data is obtained. The backpropagation wavefield is cross-correlated with the coded seismic forward wavefield to obtain a gradient field for full waveform inversion. A suitable iteration step size is obtained to perform velocity update. Then, the velocity field is learned using fast dictionary learning to obtain a suitable dictionary set. Sparse constrained denoising is performed on the seismic data to obtain a sparse constrained velocity field. The sparse constrained velocity field and the original velocity field are weighted to obtain a new velocity field. The final seismic velocity field is outputted after the number of iterations meets the requirements. Since the synthesized polar coding observation records are far less than the observation data, the computational complexity is reduced. The synthesized sparse constrained inversion velocity field has better noise resistance for the seismic data, thereby improving the accuracy of the seismic waveform inversion velocity field.

[0052] The implementation process of the present invention will be further described below with reference to specific embodiments.

[0053] Example 1

[0054] The 3D seismic velocity inversion method based on sparse constraints in this embodiment has a workflow diagram as shown in FIG. Figure 1 As shown, the following steps are included:

[0055] S1, obtain observation data, initial velocity field, source wavelet, polar coding matrix and observation system.

[0056] In actual oil and gas geophysical exploration, the observation data obtained from seismic exploration is often 3D seismic data. In this example, the 3D seismic data of a single source is shown in the following figure: Figure 2 The indicator observation system is a full-receiver observation system.

[0057] This example uses the standard depression model Figure 3 To test the effectiveness of the method, the initial velocity field is given as Figure 4 The gradient field is shown.

[0058] The source wavelet adopts the Ricker wavelet commonly used in seismic exploration, with a main frequency of 15 Hz and an amplitude of 1.0.

[0059] The observed seismic data is filtered using frequency division filtering or Wiener filtering. The present invention adopts Wiener filtering and performs Wiener filtering with Ricker wavelets having a main frequency ranging from 3 Hz to 15 Hz as target wavelets to obtain inverted observed seismic data.

[0060] The polarity coding matrix is a random polarity coding matrix, such as Figure 5 As shown, black represents positive 1 and white represents negative 1.

[0061] S2, uses seismic observation data with polarity coding matrix and observation system to synthesize super-coded observation data.

[0062] The traditional method uses a single-shot stacking calculation method, but the computational complexity of 3D seismic data is huge and inefficient. This example uses polar coding, which requires the observation system to encode and synthesize super-encoded observation data. Inversion is an iterative process, so the formula for encoding the polar coding matrix and observation data input by S1 is:

[0063]

[0064] Where, d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the coding polarity of the n-th shot in the k-th iteration, which comes from the polarity coding matrix of S1. Ns represents the total number of shots of the observation data. The computational effort of the synthetic super-encoded observation data for k iterations is equivalent to that of one shot of seismic data.

[0065] S3, synthesize the coded source wavelet using the seismic wavelet and the polarity coding matrix, and perform forward modeling using the coded source wavelet and the initial velocity field to obtain the coded forward wavefield and coded forward seismic data at each moment;

[0066] Corresponding to S2, the seismic wavelet must be polarity-encoded to obtain the coded source wavelet. The encoding matrix is consistent with the encoding matrix used in S2 and can be obtained in the following way:

[0067]

[0068] Where s n (x s ,y s ,t) is the shot point position (x s ,y s) of the earthquake source wavelet at the tth moment of the nth shot, generally the Ricker wavelet is selected, and the case of the present invention is the 15 Hz Ricker wavelet, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the coded source wavelet of k iterations, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=y s ,δ(x|x s ,y|y s )=1, otherwise it is 0.

[0069] use Figure 4 The initial velocity field and the synthesized coded seismic wavelet are shown Finite difference forward modeling can obtain the seismic wave field at each moment, store the seismic wave field, and then use the S2 synthetic super-encoded observation data consistent observation system to obtain coded forward seismic data. The calculation process is equivalent to the calculation amount of a single shot. The forward wave equation used in this example is:

[0070]

[0071] Where v is the medium velocity of the model, f is the source term of the earthquake. This case is constructed using Equation (2), u represents the seismic wave field, x and y represent the lateral coordinates, z represents the depth coordinate, and t represents time.

[0072] S4, establish the objective function of full waveform inversion based on polarity-coded sparse constrained T distribution, and calculate the difference between the synthetic super-coded observation data obtained in S2 and the coded forward seismic data.

[0073] The process of inverting the velocity field is to minimize a target functional through continuous iteration. First, a target functional is given, and then the difference is obtained according to the target functional. The present invention adopts a sparse constrained T distribution target function based on dictionary learning in the inversion, which is expressed as:

[0074]

[0075] Among them, s and ε are the degree of freedom parameter and scale factor, respectively. Their values are flexible and related to the noise of the observation data. In this case, when s = 1.5 and ε = 1.0, it is also suitable for most cases. n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the synthetic super-encoded observation data for k iterations. Represents the forward simulation of the velocity model v(x,y,z), where is the source term used, s n (x s ,y s ,t) is the shot point position (x s ,y s ) of the earthquake source wavelet at the tth moment of the nth shot, generally the Ricker wavelet is selected, and the case of the present invention is the 15 Hz Ricker wavelet, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observed data, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=y s ,δ(x|x s ,y|y s )=1, otherwise it is 0. represents the synthetic coded observation data item, d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the encoding polarity of the nth shot at the kth iteration, Ns represents the total number of observations. λ is the regularization parameter, S represents the sparse constraint on v(x, y, z), where fast dictionary learning denoising is used, E(v) is the objective function of the inversion, Ne is the total number of encoded super-observations, and can be set to 1, 5, 10, etc. A larger value increases the computational complexity. In this case, 1 is selected. j is the number of encoded super-observations.

[0076] Since full waveform inversion is an iterative computational process, based on the initial velocity field and the dominant frequency from 3 Hz to 15 Hz, 10 iterations are performed at each dominant frequency, for a total of 130 iterations. Using a T-distributed objective function effectively suppresses noise in the observation data. Polarity encoding significantly reduces computational effort, but introduces crosstalk noise. Using sparse constraints effectively suppresses crosstalk noise introduced by coded synthetic super-observation data and mitigates the impact of observation data noise. Fast dictionary learning significantly improves denoising.

[0077] S5, performs wave field back propagation on the seismic data difference obtained in S4, and cross-correlates it with the coded seismic forward wave field in S3 to obtain the gradient field of full waveform inversion.

[0078] The gradient is determined by the steepest descent method or the conjugate gradient method according to the data residual. The present invention uses the steepest descent method to update the gradient. It is necessary to derive the corresponding gradient formula according to the formula of S4. Specifically, the data residual is back-propagated to obtain the initial updated gradient value g for the kth time. k, the gradient update formula is as follows:

[0079]

[0080] Among them, v k represents the inverted velocity field of the kth iteration, and δ is the derivative operator.

[0081] S6, using the gradient field of the full waveform inversion in S5, finds the appropriate iteration step size and performs velocity update;

[0082] The parabolic interpolation method is used to calculate the optimal step size and iterate the current background velocity field. The parabolic fitting method or linear search method can be used to obtain the step size so that the error in the formula (4) corresponding to the k-th iteration inversion velocity field is minimized. In this paper, the parabolic interpolation method is used to obtain the update step size α of the k-th iteration. k , and then use the gradient in S5 to update the velocity field:

[0083] v k =v k-1 -α k g k (6)

[0084] Among them, v k represents the inverted velocity field of the kth iteration, g k is the kth weighted gradient field, 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.

[0085] S7, using fast dictionary learning to learn the velocity field and obtain a suitable dictionary set;

[0086] Because the objective function in S4 contains sparse constraints, the velocity field in S6 needs to be denoised using sparse constraints. Dictionary learning is performed on the velocity field obtained in S6. Since it is a three-dimensional model, fast dictionary learning is particularly important. Tightly constrained dictionary methods can recover n-dimensional prestack seismic data under varying signal-to-noise ratios. Compared with standard Fourier and directional transform reconstruction methods, tightly constrained dictionary methods are more likely to preserve subtle features. This paper uses a tightly constrained dictionary learning construction method to perform dictionary learning.

[0087]

[0088] Among them, G is the dictionary set to be trained, v k The velocity field is obtained by the kth iteration of formula (6), S is the sparse representation term, is the sparse coefficient corresponding to the k-th velocity field, and λ is the regularization parameter. Because this dictionary is trained based on the k-th iterative velocity field, it can more sparsely represent the current velocity field, so the denoising effect will be better.

[0089] S8, using a suitable dictionary set to perform sparse constraint denoising on the velocity field in S6, obtaining a sparse constrained velocity field, and performing velocity update;

[0090] By learning the tight constraint dictionary to perform constrained denoising on the velocity field, the iterative formula corresponding to the objective function of S3 is modified as follows:

[0091]

[0092] in, is the new updated velocity field, S is the sparse representation term, is the sparse coefficient corresponding to the k-th velocity field, Indicates that only , ε is the retention coefficient threshold, which is generally selected according to the noise level. The greater the noise, the larger the value. In this case, 1% of the maximum value is taken.

[0093] S9, repeat S2 to S8, and output the final seismic velocity field when the number of iterations meets the requirements.

[0094] As mentioned above, S2 to S8 are repeated until the number of iterations meets the preset number, and when the last iteration is completed, the target velocity field is output.

[0095] The target velocity field obtained by the above method is as follows: Figure 6 As shown, it can be seen that the depression structure is well inverted.

[0096] If the initial velocity field ( Figure 4 ), the conventional coded full waveform inversion method is performed from 2 Hz to 15 Hz using Wiener filtering, without using the sparse constraint method, and the objective function is the conventional two-norm objective functional, and the velocity field is obtained as follows Figure 7 , it can be seen that there is a certain amount of noise. The above embodiment makes the velocity field of the depression structure inversion clearer. Figure 3 The comparison shows that the velocity field inversion of the present invention is more accurate for depressions.

Claims

1. A 3D seismic velocity inversion method based on sparse constraints, characterized in that: The following steps are involved: S1, obtain seismic observation data, initial velocity field, source wavelet, polar coding matrix and observation system; S2, synthesizes super-coded observation data using seismic observation data, polar coding matrix and observation system; S3, synthesizing a coded source wavelet using the source wavelet and the polarity coding matrix, and performing forward modeling using the coded source wavelet and the initial velocity field to obtain the coded seismic forward wavefield and coded forward seismic data at each moment; S4, establishing the objective function of full waveform inversion based on polarity-coded sparse constrained T distribution, and obtaining the difference between synthetic super-coded observation data and coded forward seismic data; S5, performs wavefield backpropagation on the seismic data difference obtained in S4, and cross-correlates it with the coded seismic forward wavefield in S3 to obtain the gradient field of full waveform inversion; S6, using the gradient field of the full waveform inversion in S5, calculates the iteration step size and performs velocity update; S7, performing dictionary learning on the velocity field to obtain a dictionary set; S8, using the dictionary set to perform sparse constraint denoising on the seismic data, obtain a sparse constraint velocity field, and perform velocity update; S9, repeat S2 to S8, and output the final seismic velocity field when the number of iterations meets the requirements.

2. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1, characterized in that: In S2, the super-encoded observation data is synthesized using formula (1): In formula (1), d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the synthetic super-encoded observation data for k iterations.

3. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1, wherein: In S3, the coded source wavelet is synthesized using formula (2): In formula (2), s n (x s ,y s ,t) is the shot point position (x s ,y s )’s earthquake source wavelet at the tth moment of the nth shot, is the coding polarity of the nth shot in the kth iteration, Ns represents the total number of shots of the observation data, is the coded source wavelet of k iterations, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=y s ,δ(x|x s ,y|y s )=1, otherwise it is 0.

4. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 3, wherein: In S3, the forward wave equation applied is: In formula (3), v is the medium velocity of the model, f is the source term of the earthquake, which is constructed using formula (2), u represents the seismic wave field, x, y represent the lateral coordinates and z represents the depth coordinate, and t represents time.

5. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1, wherein: In S4, the objective function is constructed according to formula (4): Among them, s and ε are the degree of freedom parameter and scale factor respectively; Represents the forward simulation of the velocity model v(x,y,z), where is the source term used, s n (x s ,y s ,t) is the shot point position (x s ,y s ) of the earthquake source wavelet at the tth moment of the nth shot, δ(x|x s ,y|y s ) is the pulse function, when x=x s ,y=y s ,δ(x|x s ,y|y s )=1, otherwise 0; represents the synthetic coded observation data item, d n (x,y,t) is the seismic data shot record of the tth time of the nth shot at the (x,y) coordinates, is the encoding polarity of the nth shot in the kth iteration, Ns represents the total number of observations; λ is the regularization parameter, S represents the sparse constraint on v(x, y, z), and the fast dictionary learning denoising constraint is used, E(v) is the objective function of the inversion, Ne is the total number of encoded super-observations, and j is the sequence number of the encoded super-observation data.

6. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1 or 2, characterized in that: In S5, the gradient is updated according to formula (5): Among them, g k is the kth weighted gradient field, E(V k ) is the objective function of the kth inversion, v k represents the inverted velocity field of the kth iteration, and δ is the derivative operator.

7. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1, characterized in that: In S6, the current background velocity field is iterated using the calculated optimal iteration step size.

8. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 7, characterized in that: In S6, the speed is updated according to formula (6): v k =v k-1 -α k g k (6) Among them, v k represents the inverted velocity field of the kth iteration, α k is the update step size of the kth iteration, g k is the kth weighted gradient field, 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.

9. The sparse constraint-based three-dimensional seismic velocity inversion method according to claim 1, wherein: In S7, dictionary learning is performed according to formula (7): Among them, G is the dictionary set to be trained, v k The velocity field is obtained for the kth iteration, where (x, y, z) represents the lateral coordinates, z represents the depth coordinate, and S is the sparse representation term. is the sparse coefficient corresponding to the k-th velocity field, and λ is the regularization parameter.

10. The sparse constraint-based three-dimensional seismic velocity inversion method according to any one of claims 1 to 5, characterized in that: In S8, the speed is updated using formula (8): in, is the new updated velocity field, S is the sparse representation term, is the sparse coefficient corresponding to the k-th velocity field, Indicates that only , ε is the retention coefficient threshold.

Images