Earthquake velocity inversion method

Abstract

The invention relates to a seismic velocity inversion method, and belongs to the technical field of seismic exploration. According to the seismic velocity inversion method, seismic data are divided into different frequency bands, and then full-waveform inversion is carried out by using a low-frequency band to obtain a velocity field; fusing the velocity field with the well velocity field, and performing grid chromatography inversion to obtain a well constraint velocity field; and finally, carrying out grid chromatography inversion and full waveform inversion on the well constraint velocity field to obtain a target velocity field. The seismic velocity inversion method can solve the problem of cycle jump caused by low-frequency loss, has good adaptability to data with low-frequency loss in seismic data, and can improve the stability and accuracy of the seismic data inversion velocity field.

Description

TECHNICAL FIELD

[0001] The present application relates to a seismic velocity inversion method, belonging to the technical field of seismic exploration. BACKGROUND

[0002] In the current oil and gas seismic exploration field, the accurate construction of the velocity field is the basis for subsequent migration imaging and geological interpretation. Full waveform inversion can accurately invert the underground seismic velocity field by comprehensively utilizing the information in the seismic data. However, full waveform inversion has the phenomenon of cycle skipping and is easily affected by local extreme values. Full waveform inversion of seismic velocity field has many problems such as being easily disturbed by noise and poor inversion stability. Therefore, it is necessary to perform segmented inversion on seismic data. However, due to the frequent absence of low-frequency data in the collected data, it causes certain disturbance to the velocity inversion of seismic data. Due to the multi-solution nature of inversion, in the inversion process, due to the absence of low-frequency data and the inaccuracy of the background field, it causes certain disturbance to the velocity inversion of seismic data. Therefore, it is necessary to accurately construct the background velocity field.

[0003] For example, Chinese patent document CN108680957B discloses a weighted local cross-correlation time-frequency domain phase inversion method. The inversion method disclosed in the patent document is to introduce the time-frequency domain phase information of seismic data into the cross-correlation objective function, which 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, the addition of a weight factor in the time-frequency domain objective function can well enhance the anti-noise ability and the stability of 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 can then be 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. In the test of simultaneous absence of low-frequency components and strong Gaussian background noise, the weighted local cross-correlation time-frequency domain phase inversion method has strong anti-noise ability and does not depend on the initial velocity model, etc. However, the extraction of phase information causes complex operation.

[0004] The patent document CN110888158A discloses a full waveform inversion method based on RTM constraint. The method disclosed in the patent document establishes a velocity model by using a full waveform inversion gradient decomposition method and a RTM constrained full waveform inversion method, and performs flat gather processing on the velocity model by using a time difference correction method to obtain a high-precision migration velocity field. The method disclosed in the patent document fully utilizes low-frequency information and deep reflection wave energy by using the full waveform inversion gradient decomposition technology and the RTM constrained full waveform inversion technology, and improves the precision and stability of the inversion result. The method disclosed in the patent document introduces an iterative inversion technology of time difference correction, reduces the travel time error caused by velocity abnormal bodies by using the time difference correction technology, improves the applicability of the velocity model to migration imaging, and improves the accuracy of inversion by using the RTM constraint. However, the method needs additional reverse-time migration calculation and processing, which causes complex operation and cannot fully utilize logging and other geological information.

[0005] Therefore, there is an urgent need to provide a seismic velocity inversion method which is simple to operate and has high stability and accuracy. SUMMARY

[0006] The present application aims to provide a seismic velocity inversion method which can solve the problem of complex operation of the current seismic velocity inversion method.

[0007] To achieve the above object, the technical scheme of the seismic velocity inversion method of the present application is as follows:

[0008] A seismic velocity inversion method comprises the following steps:

[0009] S1, performing full waveform inversion on an initial velocity field by using low-frequency seismic observation data to obtain a velocity field;

[0010] S2, fusing the velocity field obtained in step S1 with a well velocity field to obtain a well-constrained velocity field;

[0011] S3, performing grid tomography inversion on the well-constrained velocity field to obtain a tomography inversion velocity field;

[0012] S4, performing full waveform inversion on the tomography inversion velocity field;

[0013] S5, repeatedly operating steps S1 to S4 for multiple times, and the frequency of the low-frequency data is increased in turn in the repeated operation, and the initial velocity field used in each operation is the velocity field obtained at the end of the previous operation.

[0014] The seismic velocity inversion method of the present application divides the seismic data into different frequency bands, then uses the low frequency band to obtain the velocity field by full waveform inversion; then fuses the velocity field with the well velocity field to obtain the well-constrained velocity field by grid tomography inversion; finally, the well-constrained velocity field is subjected to grid tomography inversion and full waveform inversion to obtain the target velocity field. The seismic velocity inversion method of the present application can solve the problem of cycle skipping caused by low frequency missing, has good adaptability to data with low frequency missing in seismic data, and can improve the stability and accuracy of the inversion velocity field of seismic data.

[0015] Preferably, the initial velocity field is obtained by grid tomography inversion of the seismic observation data and the original velocity field.

[0016] Preferably, the frequency of the low frequency band seismic observation data in step S1 is 5Hz, and the number of times of repeating the steps S1 to S4 is 3, and the frequency of the low frequency band seismic observation data in each repeated operation is 7Hz, 9Hz and 10Hz in turn.

[0017] Preferably, the method for full waveform inversion of the initial velocity field using the low frequency band data comprises the following steps: first, determining the low frequency band seismic data and the seismic wavelet, then using the seismic wavelet and the initial velocity field to perform forward simulation to obtain the seismic forward wave field and the forward seismic data at each time, then according to the target functional of the full waveform inversion, calculating the difference between the seismic observation data and the forward seismic data, then performing wave field back propagation on the calculated difference, then performing cross-correlation with the seismic forward wave field to obtain the gradient field of the full waveform inversion, and finally, using the gradient field to calculate the appropriate iteration step length to perform velocity update.

[0018] Preferably, the forward wave equation used for forward simulation using the seismic wavelet and the initial velocity field is as shown in the following formula:

[0019]

[0020] In the formula, m represents the particle velocity, f represents the source term, u is the seismic wave field, x and z represent the coordinates, and t represents the time.

[0021] Preferably, the target functional of the full waveform inversion is the two-norm target functional.

[0022] Preferably, the iteration step length is obtained by parabolic interpolation calculation.

[0023] Preferably, the calculation formula of the well-constrained velocity field is as shown in the following formula:

[0024] v rh =α w v well +(1-α w )v fwi

[0025] wherein v rh is the well constrained velocity field, v well is the well velocity field, v fwi is the velocity field obtained in step S2, a w is the well constraint weighting coefficient.

[0026] Preferably, the well constraint weighting coefficient is determined according to the following formula:

[0027]

[0028] wherein H is the depth of the well, h is the depth of the velocity field, f max is the maximum frequency of the full waveform inversion, f min is the minimum frequency of the full waveform inversion, and a is the initial coefficient of the well constraint. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1 is a flowchart of the seismic velocity inversion method according to an embodiment of the present application;

[0030] Figure 2 is a schematic diagram of a single shot record according to an embodiment of the present application;

[0031] Figure 3 is a schematic diagram of an original velocity field according to an embodiment of the present application;

[0032] Figure 4 is a schematic diagram of a background velocity field according to an embodiment of the present application;

[0033] Figure 5 is a schematic diagram of two-dimensional seismic data (standard model, Marmousi model) according to an embodiment of the present application;

[0034] Figure 6 is a schematic diagram of the velocity field inversion result obtained according to an embodiment of the present application;

[0035] Figure 7 is a schematic diagram of the velocity field inversion result obtained by using a conventional inversion method according to an embodiment of the present application. DETAILED DESCRIPTION

[0036] The seismic velocity inversion method according to the present application is an inventive creation. The seismic velocity inversion method according to the present application comprises the following steps:

[0037] S1, obtaining an initial velocity field;

[0038] S2, performing frequency division on seismic data, performing full waveform inversion on the initial velocity field by using low frequency band data to obtain a velocity field;

[0039] S3, fusing the velocity field obtained in step S2 with a well velocity field to obtain a well constrained velocity field.

[0040] S4, grid tomography inversion is performed on the well-constrained velocity field to obtain a tomography inversion velocity field;

[0041] S5, full waveform inversion is performed on the tomography inversion velocity field.

[0042] The seismic velocity inversion method can solve the problem of cycle skipping caused by low frequency missing, has good adaptability to data with missing low frequencies in seismic data, and can improve the stability and accuracy of the inversion velocity field of seismic data.

[0043] In some preferred embodiments, the initial velocity field is obtained by grid tomography inversion of observed seismic data and an original velocity field.

[0044] In some preferred embodiments, the frequency of the frequency-divided seismic data is 5 Hz, 7 Hz, 9 Hz or 10 Hz.

[0045] In some preferred embodiments, the method for performing full waveform inversion on the initial velocity field using low frequency data comprises the following steps: first, determining low frequency seismic data and a seismic wavelet, then performing forward modeling on the initial velocity field using the seismic wavelet to obtain seismic forward wave field and forward seismic data at each time, then calculating the difference between the observed seismic data and the forward seismic data according to the target functional of the full waveform inversion, then performing wave field back propagation on the calculated difference, then performing cross-correlation with the seismic forward wave field to obtain the gradient field of the full waveform inversion, and finally, calculating a suitable iteration step length using the gradient field to update the velocity.

[0046] In some preferred embodiments, the forward wave equation used for performing forward modeling on the initial velocity field using the seismic wavelet is as shown in the following formula:

[0047]

[0048] In the formula, m represents the particle velocity, f represents the source term, u is the seismic wave field, x and z represent coordinates, and t represents time.

[0049] In some preferred embodiments, the target functional of the full waveform inversion is a two-norm target functional.

[0050] In some preferred embodiments, the iteration step length is obtained by parabolic interpolation calculation.

[0051] In some preferred embodiments, the calculation formula of the well-constrained velocity field is as shown in the following formula:

[0052] v rh =α w v well +(1-α w )v fwi

[0053] In the formula, v rh For the well-constrained velocity field, v well For the well velocity field, v fwi For the velocity field obtained in step S2, α w The weighting coefficients for well constraints.

[0054] In some preferred embodiments, the weighting coefficients for well constraints are determined using the following formula:

[0055]

[0056] In the formula, H is the depth of the well, h is the depth of the velocity field, and f max f is the maximum frequency for full waveform inversion. min The minimum frequency for full waveform inversion is τ, and the initial coefficient of the well constraint is 0.1 to 0.5.

[0057] The technical solution of the present invention will be further described below with reference to specific embodiments.

[0058] Example

[0059] The seismic velocity inversion method in this embodiment, such as Figure 1 As shown, the specific steps include:

[0060] S1, using seismic observation data and the original velocity field, performs grid tomography inversion to obtain the initial velocity field.

[0061] The seismic observation data in this embodiment is two-dimensional seismic data, and the single-shot record is as follows: Figure 2 As shown, using the original velocity field (such as...) Figure 3 As shown, a gather is generated, then the seismic attribute volume required for mesh tomography is picked out using the gather, and then the residual velocity volume is generated using the seismic attribute volume. The mesh tomography is iterated 200 times to solve the problem, and then the background velocity field is smoothed to obtain the result (as shown). Figure 4 (As shown), the remaining velocity body and the background velocity field are then algebraically added together to obtain the initial velocity field.

[0062] S2 divides the seismic data into different frequency bands, and then uses the low-frequency seismic data to perform full waveform inversion to obtain the velocity field.

[0063] S21, Acquire low-frequency seismic data and source wavelet (seismic wavelet).

[0064] The schematic diagram of the two-dimensional seismic observation data (standard model, Marmousi model) used in this embodiment is shown below. Figure 5 As shown, the effectiveness of the test method is verified using the Marmousi model; the two-dimensional seismic observation data is filtered to obtain low-frequency seismic observation data with a frequency of 5Hz. In this embodiment, the source wavelet (seismic wavelet) selected is the Rick wavelet commonly used in seismic exploration, and the dominant frequency is equal to the frequency of the low-frequency seismic observation data, which is 5Hz.

[0065] S22, forward modeling is performed using seismic wavelets and the initial velocity field to obtain the forward modeling wavefield and forward modeling seismic data at each time step.

[0066] Using the initial velocity field obtained in step S1 and the seismic wavelet selected in step S21, a finite-difference forward modeling simulation was performed to obtain seismic forward modeling data. The results are as follows: Figure 5 As shown; the forward wave equation used in this embodiment is shown in Equation 1:

[0067]

[0068] In Equation 1, m represents the particle velocity, f represents the source term, u is the seismic wave field, x and z represent coordinates, and t represents time.

[0069] S23, Based on the target functional obtained from the full waveform inversion, calculate the difference between the low-frequency seismic observation data with a frequency of 5Hz obtained in step S21 and the forward-modeled seismic data obtained in step S22.

[0070] The process of inverting the velocity field involves iteratively minimizing a target functional (by continuously iterating to minimize the target functional). First, a target functional is given, and then the difference is calculated based on the target functional. In this embodiment, the L2 norm target functional, which is commonly used in inversion, is used, as shown in Equation 2:

[0071]

[0072] In Equation 2, m is the velocity parameter model, representing the particle velocity; E(m) is the L2 norm of the data residual corresponding to the velocity field (particle velocity) m; and u is the objective functional in this embodiment. obs and u ca l represents the field observation seismic record at time t ( Figure 4 The single-shot record shown is compared with the forward modeling record (forward modeling seismic data obtained in step S22), x s x is the epicenter. r This is the receiving point.

[0073] S24, wave field back propagation is performed on the difference between the low frequency observed seismic data with a frequency of 5 Hz obtained in step S23 and the forward seismic data, and then cross correlation is performed on the seismic forward wave field obtained in step S22, to obtain a gradient field of full waveform inversion

[0074] The process of obtaining the gradient field of full waveform inversion through cross correlation is the process of velocity field updating. According to the data difference, the gradient is determined by using the steepest descent method or the conjugate gradient method. In this embodiment, the steepest descent method is used for gradient updating, and the velocity field model is iterated based on the gradient. According to the original velocity field (such as Figure 3 shown), the main frequency 5 Hz is iterated for 10 times, and the update gradient value g k of the kth iteration is obtained by back propagation on the data difference. The gradient update formula is shown in formula 3:

[0075]

[0076] In formula 3, m k represents the inversion velocity field of the kth iteration, E(m k ) is the data residual norm corresponding to the inversion velocity field m k of the kth iteration, δ is the derivative operator, m is the particle velocity, and the main frequency of the Ricker wavelet used in the kth iteration is the main frequency of the Ricker wavelet in step S21.

[0077] S25, using the gradient field to obtain a suitable iteration step length for velocity updating

[0078] The optimal step length is calculated by using parabolic interpolation, and the current background velocity field is iterated for 10 times, so that the error (i.e. the data residual norm in step S23) corresponding to the inversion velocity field of the kth iteration in step S23 is minimized, and the velocity field is output.

[0079] S3, fusing the velocity field obtained in step S2 with the well velocity field to obtain a well-constrained velocity field

[0080] First, the well curve is extrapolated and interpolated along the layer to obtain a well velocity field; then in the velocity inversion updating process, the full waveform inversion velocity field obtained in step S2 and the well velocity field are weighted and fused to obtain a well-constrained velocity field; wherein the calculation formula of the well-constrained velocity field is shown in formula 4:

[0081] v rh = α w v well + (1- α w ) v fwi Formula 4

[0082] In formula 4, v rh is the well-constrained velocity field, v well is the well velocity field, v fwi is the velocity field obtained in step S2, and αw The well-constrained weighting coefficient matrix is determined according to the principle that the well velocity information is reasonably utilized, the macro effect of the well velocity information on tomography is fully played, and the fine effect of full waveform inversion is also played. When the depth is greater than the well velocity, the well constraint is not affected. When the background velocity inversion is dominant, the well constraint is strengthened, and the lower the frequency is, the greater the contribution to the background field inversion is. In the embodiment, the determination method of the well-constrained weighting coefficient matrix is shown in formula 5:

[0083]

[0084] In formula 5, H is the depth of the well, h is the depth of the velocity field, f max is the maximum frequency of full waveform inversion, which is 10 Hz in the embodiment (in the embodiment, considering the influence of the velocity field on imaging, it is found that the maximum frequency of full waveform inversion is most suitable to be set as 10 Hz), f min is the minimum frequency of full waveform inversion, which is 5 Hz in the embodiment, and τ is the initial coefficient of the well constraint, which is 0.2 in the embodiment (the initial coefficient of the well constraint is generally 0.1-0.5, and experimental results show that the effect is best when the value is 0.2, that is, the absolute value error of the velocity is smallest). As can be seen from formula 5, when the depth is greater than the well velocity, the well-constrained weighting coefficient is 0, and when the depth is less than the well depth, the well-constrained weighting coefficient is τ*(f max -f) / (f max -f min ).

[0085] S4, using the well-constrained velocity field obtained in step S3 to perform grid tomography inversion using the gather, to obtain a tomography inversion velocity field;

[0086] The gather is generated using the well-constrained velocity field obtained in step S3, the gather is not flat at this time, then the gather is used to pick up the grid tomography seismic attribute volume, then the seismic attribute volume is used to generate a residual velocity volume, and the iteration is performed 200 times to solve, and the inversion velocity field is output. The gather is flattened to obtain the tomography inversion velocity field.

[0087] S5, performing full waveform inversion on the tomography inversion velocity field obtained in step S4 to obtain a velocity field

[0088] The step specifically repeats steps S21-S25 once to perform full waveform inversion using the tomography inversion velocity field obtained in step S4 to obtain a new velocity field.

[0089] In the specific operation, first, the source wavelet (i.e. the seismic wavelet, the Ricker wavelet with the main frequency of 5 Hz) is forward modeled with the tomographic inversion velocity field obtained in step S4 to obtain the seismic forward wave field and the forward seismic data at each time (the seismic forward wave field and the forward seismic data obtained by forward modeling based on the tomographic inversion velocity field obtained in step S4); then, the difference between the low-frequency seismic observation data obtained in step S21 and the forward seismic data (the forward seismic data obtained by forward modeling based on the tomographic inversion velocity field obtained in step S4) is calculated according to the target functional of the full waveform inversion; then, the obtained difference is wave field back propagated, and then cross correlated with the seismic forward wave field (the seismic forward wave field obtained by forward modeling based on the tomographic inversion velocity field obtained in step S4) to obtain the gradient field of the full waveform inversion, and then the suitable iteration step is calculated by using the gradient field, the velocity is updated, and the new velocity field is obtained.

[0090] S6, increase the data frequency band, repeat steps S2 to S5, and output the inversion velocity field of the target

[0091] In this step, in the specific operation, the main frequency band of the seismic data is increased from 5 Hz to 7 Hz, and then sequentially increased to 9 Hz and 10 Hz, and then the new seismic data is obtained, and then the steps S2 to S5 are repeated to obtain the target velocity field. The results are shown in Figure 6

[0092] The specific operation steps are as follows:

[0093] ①First, the two-dimensional seismic observation data is filtered to obtain low-frequency seismic observation data with a frequency of 7 Hz, and then the seismic wavelet (the Ricker wavelet with the main frequency of 7 Hz) is forward modeled with the velocity field obtained in step S5 to obtain the seismic forward wave field and the forward seismic data at each time, and then the difference between the low-frequency seismic observation data with the frequency of 7 Hz and the forward seismic data is calculated according to the target functional of the full waveform inversion, and then the difference is wave field back propagated, and then cross correlated with the seismic forward wave field to obtain the gradient field of the full waveform inversion, and then the suitable iteration step is calculated by using the gradient field, the velocity is updated, and then the well-constrained velocity field is obtained by fusing the obtained velocity field and the well velocity field, and then the well-constrained velocity field is grid tomographic inverted by using the gather to obtain the tomographic inversion velocity field, and then the tomographic inversion velocity field is full waveform inverted to obtain the velocity field;

[0094] ​②Then the two-dimensional seismic observation data is filtered to obtain low-frequency seismic observation data with a frequency of 9Hz, then forward modeling is carried out by using a seismic wavelet (a Ricker wavelet with a main frequency of 9Hz) and the velocity field obtained in step ① (the velocity field obtained by using low-frequency seismic observation data with a frequency of 7Hz), to obtain seismic forward wave field and forward seismic data at each time, then the difference between the low-frequency seismic observation data with a frequency of 9Hz and the forward seismic data is obtained according to the objective function of full waveform inversion, then the difference is wave field back propagation, then cross correlation is carried out with the seismic forward wave field to obtain the gradient field of full waveform inversion, then the appropriate iteration step is obtained by using the gradient field, velocity updating is carried out, then the obtained velocity field is fused with the well velocity field to obtain the well-constrained velocity field, then the well-constrained velocity field is used for grid tomography inversion by using the gather to obtain the tomography inversion velocity field, then the tomography inversion velocity field is used for full waveform inversion to obtain the velocity field;

[0095] ③Finally, the two-dimensional seismic observation data is filtered to obtain low-frequency seismic observation data with a frequency of 10Hz, then forward modeling is carried out by using a seismic wavelet (a Ricker wavelet with a main frequency of 10Hz) and the velocity field obtained in step ② (the velocity field obtained by using low-frequency seismic observation data with a frequency of 9Hz), to obtain seismic forward wave field and forward seismic data at each time, then the difference between the low-frequency seismic observation data with a frequency of 10Hz and the forward seismic data is obtained according to the objective function of full waveform inversion, then the difference is wave field back propagation, then cross correlation is carried out with the seismic forward wave field to obtain the gradient field of full waveform inversion, then the appropriate iteration step is obtained by using the gradient field, velocity updating is carried out, then the obtained velocity field is fused with the well velocity field to obtain the well-constrained velocity field, then the well-constrained velocity field is used for grid tomography inversion by using the gather to obtain the tomography inversion velocity field, then the tomography inversion velocity field is used for full waveform inversion to obtain the target velocity field.

[0096] In the inversion of the Marmousi model in step S2, Figure 5 the original velocity field model Figure 3 is used to obtain the velocity field inversion result through repeated iterations, as shown in Figure 6 , and it can be known from Figure 6 that the seismic velocity inversion method has better inversion effect on complex structure horizons, and can accurately invert the velocity field, and eliminates the unstable influence in the conventional inversion method.

[0097] In order to compare the seismic velocity inversion method with the conventional inversion method, the original velocity field is processed by using the conventional inversion method, that is, first, grid tomography inversion is carried out by using the seismic observation data and the original velocity field to obtain the initial velocity field, then full waveform inversion (the same as the full waveform inversion in step S2) is carried out, and the result is shown in Figure 7 , and it can be known fromFigure 7 It can be known that the stability of the conventional inversion method is poor. Compared with the conventional method, the sum of absolute value errors of the velocity values (predicted velocity values) of the test points obtained by using the method of the application (the sum of absolute value errors is equal to the sum of velocity absolute value errors of each test point, and the velocity absolute value error of the test point = |predicted velocity value - actual velocity value| / actual velocity value*100%) is reduced from 25% to 5%; wherein, the test model selected in the method of the application and the conventional method is the same, and the test points are all data points in the test model.

Claims

1. A seismic velocity inversion method, characterized in that, Includes the following steps: S1, the initial velocity field is obtained by full waveform inversion using low-frequency seismic observation data; S2, fuse the velocity field obtained in step S1 with the well velocity field to obtain the well constrained velocity field; S3, perform grid tomographic inversion on the well-constrained velocity field to obtain the tomographic inversion velocity field; S4, perform full waveform inversion of the tomographically inverted velocity field; S5, repeat steps S1 to S4 multiple times. During repeated operations, the frequency of the low-frequency data increases sequentially. The initial velocity field used in a certain operation is the velocity field obtained after the previous operation ends.

2. The seismic velocity inversion method as described in claim 1, characterized in that, The initial velocity field was obtained by grid tomography inversion of seismic observation data and the original velocity field.

3. The seismic velocity inversion method as described in claim 1, characterized in that, The frequency of the low-frequency seismic observation data in step S1 is 5Hz. Steps S1 to S4 are repeated 3 times, and the frequencies of the low-frequency seismic observation data in each repetition are 7Hz, 9Hz, and 10Hz respectively.

4. The seismic velocity inversion method as described in any one of claims 1-3, characterized in that, The method for full waveform inversion of the initial velocity field using low-frequency data includes the following steps: First, determine the low-frequency seismic data and seismic wavelet. Then, perform forward modeling using the seismic wavelet and the initial velocity field to obtain the forward modeling wavefield and forward modeling seismic data at each time step. Next, calculate the difference between the seismic observation data and the forward modeling seismic data based on the objective functional of the full waveform inversion. Then, backpropagate the calculated difference to the wavefield and cross-correlate it with the forward modeling seismic wavefield to obtain the gradient field of the full waveform inversion. Finally, use the gradient field to determine an appropriate iteration step size for velocity update.

5. The seismic velocity inversion method as described in claim 4, characterized in that, The forward wave modeling equation used for forward modeling using seismic wavelets and the initial velocity field is shown below: In the formula, m represents the particle velocity, f represents the source term, u is the seismic wave field, x and z represent coordinates, and t represents time.

6. The seismic velocity inversion method as described in claim 4, characterized in that, The objective functional for full waveform inversion is a norm 2 objective functional.

7. The seismic velocity inversion method as described in claim 4, characterized in that, The iteration step size is obtained by parabolic interpolation.

8. The seismic velocity inversion method as described in any one of claims 1-3, characterized in that, The formula for calculating the well-constrained velocity field is as follows: v rh =α w v well +(1-α w )v fwi In the formula, v rh For the well-constrained velocity field, v well For the well velocity field, v fwi For the velocity field obtained in step S2, α w The weighting coefficients for well constraints.

9. The seismic velocity inversion method as described in claim 8, characterized in that, The method for determining the weighting coefficients of well constraints is shown in the following formula: In the formula, H is the depth of the well, h is the depth of the velocity field, and f max f is the maximum frequency for full waveform inversion. min τ is the minimum frequency for full waveform inversion, and τ is the initial coefficient of the well constraint.

Images