A seismic velocity inversion method

Abstract

This invention relates to a seismic velocity inversion method, belonging to the field of seismic exploration technology. The seismic velocity inversion method of this invention divides seismic data into different frequency bands, then uses the low-frequency band for full waveform inversion to obtain the velocity field; next, the velocity field is fused with the well velocity field, and a grid tomographic inversion is performed to obtain the well-constrained velocity field; finally, the well-constrained velocity field is subjected to grid tomographic inversion and full waveform inversion to obtain the target velocity field. This seismic velocity inversion method can solve the cycle skipping problem caused by missing low-frequency data, has good adaptability to low-frequency missing data in seismic data, and can improve the stability and accuracy of the velocity field inverted from seismic data.

Description

Technical Field

[0001] This invention relates to a seismic velocity inversion method, belonging to the field of seismic exploration technology. Background Technology

[0002] In the current field of oil and gas seismic exploration, the accurate construction of the velocity field is the foundation for subsequent migration imaging and geological interpretation. Full waveform inversion, by comprehensively utilizing information from seismic data, can accurately invert the subsurface seismic velocity field. However, full waveform inversion suffers from cycle skipping and is easily affected by local extrema. Full waveform inversion of the seismic velocity field also presents numerous problems, such as susceptibility to noise interference and poor inversion stability. Therefore, it is necessary to perform segmented inversion of the seismic data. However, the acquisition data often lacks low-frequency data, causing some difficulties for the velocity inversion. Due to the ambiguity of the inversion process, the lack of low-frequency data, coupled with an inaccurate background field, further complicates the velocity inversion. Therefore, it is essential 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. This method incorporates time-frequency domain phase information from seismic data into the cross-correlation objective function, hence the name "weighted local cross-correlation time-frequency domain phase inversion method." Introducing phase information alleviates the dependence of full-waveform inversion on the initial velocity model. Furthermore, adding weighting factors to the time-frequency domain objective function significantly enhances noise resistance and inversion stability. In the low-frequency band weighted local cross-correlation time-frequency domain phase inversion method, a good initial velocity model can be obtained, which is then used in the high-frequency band weighted local cross-correlation time-frequency domain phase inversion method, ultimately yielding high-resolution inversion results. Tests with both missing low-frequency components and strong Gaussian background noise demonstrate that the weighted local cross-correlation time-frequency domain phase inversion method has advantages such as strong noise resistance and independence from the initial velocity model; however, the need to extract phase information complicates the operation.

[0004] Chinese patent document CN110888158A discloses a full-waveform inversion method based on RTM constraints. This method utilizes a full-waveform inversion gradient decomposition method and an RTM-constrained full-waveform inversion method to establish a velocity model. A time-difference correction method is then used to flatten the gathers in the velocity model, resulting in a high-precision migrated velocity field. This method, through full-waveform inversion gradient decomposition and RTM-constrained full-waveform inversion techniques, fully utilizes low-frequency information and deep reflected wave energy, improving the accuracy and stability of the inversion results. The method also introduces an iterative inversion technique with time-difference correction to reduce travel time errors caused by velocity anomalies, improving the applicability of the velocity model to migration imaging. While RTM constraints enhance the accuracy of the inversion, the method requires additional reverse-time migration calculations, leading to operational complexity and limiting the full utilization of geological information such as well logging data.

[0005] Therefore, there is an urgent need to provide a seismic velocity inversion method that is easy to operate and has high stability and accuracy. Summary of the Invention

[0006] The purpose of this invention is to provide a seismic velocity inversion method that can solve the problem of complex operation in current seismic velocity inversion methods.

[0007] To achieve the above objectives, the technical solution adopted by the seismic velocity inversion method of the present invention is as follows:

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

[0009] S1, the initial velocity field is obtained by full waveform inversion using low-frequency seismic observation data;

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

[0011] S3, perform grid tomographic inversion on the well-constrained velocity field to obtain the tomographic inversion velocity field;

[0012] S4, perform full waveform inversion of the tomographically inverted velocity field;

[0013] 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.

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

[0015] Preferably, the initial velocity field is obtained by grid tomography inversion of 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 steps S1 to S4 are repeated is 3 times, with the frequency of the low-frequency band seismic observation data being 7Hz, 9Hz, and 10Hz respectively during each repetition.

[0017] Preferably, the method for full waveform inversion of the initial velocity field using low-frequency band data includes the following steps: First, determine the low-frequency band seismic data and seismic wavelet. Then, perform forward modeling using the seismic wavelet and the initial velocity field to obtain the seismic 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, perform wavefield backpropagation on the calculated difference and cross-correlate it with the seismic forward modeling 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.

[0018] Preferably, the forward wave equation used for forward modeling using seismic wavelets and the initial velocity field is as follows:

[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 coordinates, and t represents time.

[0021] Preferably, the objective functional for full waveform inversion is a L2 objective functional.

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

[0023] Preferably, the formula for calculating the well-constrained velocity field is as follows:

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

[0025] 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.

[0026] Preferably, the weighting coefficients for well constraints are determined using the following formula:

[0027]

[0028] 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. Attached Figure Description

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

[0030] Figure 2 This is a schematic diagram of single-shot recording in an embodiment of the present invention;

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

[0032] Figure 4 This is a schematic diagram of the background velocity field in an embodiment of the present invention;

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

[0034] Figure 6 This is a schematic diagram of the velocity field inversion results obtained in an embodiment of the present invention;

[0035] Figure 7 This is a schematic diagram of the velocity field inversion results obtained using the conventional inversion method in this invention. Detailed Implementation

[0036] The seismic velocity inversion method of this invention is a pioneering invention. The seismic velocity inversion method of this invention includes the following steps:

[0037] S1, obtain the initial velocity field;

[0038] S2, the seismic data is divided into frequencies, and the initial velocity field is inverted using the low-frequency data to obtain the velocity field;

[0039] S3, fuse the velocity field obtained in step S2 with the well velocity field to obtain the well constrained velocity field;

[0040] S4. Perform mesh tomographic inversion on the well-constrained velocity field to obtain the tomographic inversion velocity field;

[0041] S5 performs a full waveform inversion of the tomographically inverted velocity field.

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

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

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

[0045] In some preferred embodiments, 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 the source 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 observed seismic data and the forward modeling seismic data based on the objective functional of the full waveform inversion. Then, perform wavefield backpropagation on the calculated difference 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.

[0046] In some preferred embodiments, the forward wave modeling equation used for forward modeling with seismic wavelets and initial velocity fields is shown below:

[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 for full waveform inversion is a L2-norm target functional.

[0050] In some preferred embodiments, the iteration step size is obtained using parabolic interpolation.

[0051] In some preferred embodiments, the formula for calculating the well-constrained velocity field is as follows:

[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, the difference between the 5Hz low-frequency observed seismic data obtained in step S23 and the forward modeled seismic data is backpropagated using wavefield inversion, and then cross-correlated with the seismic forward modeled wavefield obtained in step S22 to obtain the gradient field of the full waveform inversion.

[0074] The process of obtaining the gradient field from the full waveform inversion through cross-correlation is the process of updating the velocity field. Based on the data difference, the gradient is determined using either the steepest descent method or the conjugate gradient method. This embodiment uses the steepest descent method for gradient updating, and iterates the velocity field model based on the gradient, according to the original velocity field (e.g., ...). Figure 3 As shown), the main frequency is 5Hz and the iteration is performed 10 times. The data difference is backpropagated to obtain the update gradient value g of the kth iteration. k The gradient update formula is shown in Equation 3:

[0075]

[0076] In Equation 3, m k Let E(m) represent the inversion velocity field of the k-th iteration. k Let m be the inversion velocity field of the kth iteration. k The corresponding data residual L2 norm, δ is the differentiation operator, m is the particle velocity, and the dominant frequency of the Ricker wavelet used in the k-th iteration is the dominant frequency of the Ricker wavelet in step S21.

[0077] S25, use the gradient field to determine a suitable iteration step size and perform velocity update.

[0078] The optimal step size is calculated using parabolic interpolation. The current background velocity field is iterated 10 times to minimize the error in step S23 (i.e., the L2 norm of the data residual in step S23) corresponding to the velocity field in the k-th iteration, and the velocity field is output.

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

[0080] First, the well velocity field is obtained by extrapolating the well curve along the layer; then, during the velocity inversion update process, the full waveform inversion velocity field obtained in step S2 is weighted and fused with the well velocity field to obtain the well-constrained velocity field; the calculation formula for the well-constrained velocity field is shown in Equation 4:

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

[0082] In Equation 4, 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 This is the weighting coefficient matrix for well constraints. The principle for selecting the weighting coefficients for well constraints is to make reasonable use of well velocity information, give full play to the macroscopic role of well velocity information in tomography, and also give full play to the fine role of full waveform inversion. When the depth is greater than the well velocity, it is not affected by the well velocity. When the background velocity inversion is dominant, the well constraint should be strengthened, and the lower the frequency, the greater the contribution to the background field inversion. In this embodiment, the method for determining the weighting coefficient matrix of well constraints is shown in Equation 5:

[0083]

[0084] In Equation 5, H is the depth of the well, h is the depth of the velocity field, and f max The maximum frequency for full waveform inversion is 10Hz in this embodiment (in this embodiment, considering the influence of the velocity field on imaging, setting the maximum frequency for full waveform inversion to 10Hz is found to be most suitable), f min The minimum frequency for full waveform inversion is 5Hz in this embodiment. τ is the initial coefficient of the well constraint, which is 0.2 in this embodiment (the initial coefficient of the well constraint is generally between 0.1 and 0.5; experimental results show that a value of 0.2 yields the best results, i.e., the smallest absolute velocity error). As shown in Equation 5, when the depth is greater than the well velocity, the weighting coefficient of the well constraint is 0; when the depth is less than the well depth, the weighting coefficient of the well constraint is τ*(f). max -f) / (f max -f min ).

[0085] S4. The well-constrained velocity field obtained in step S3 is used to perform grid tomographic inversion using gathers to obtain the tomographic inversion velocity field.

[0086] The well-constrained velocity field obtained in step S3 is used to generate a gather. At this time, the gather is not flat. Then, the seismic attribute volume of the mesh tomography is picked up using the gather. The residual velocity volume is then generated using the seismic attribute volume. The solution is iterated 200 times and the inverted velocity field is output. The gather is flattened and the tomographic inverted velocity field is obtained.

[0087] S5, perform full waveform inversion on the tomographic inversion velocity field obtained in step S4 to obtain the velocity field.

[0088] Specifically, this step involves using the tomographic inversion velocity field obtained in step S4, and repeating steps S21 to S25 once to perform full waveform inversion, thereby obtaining a new velocity field.

[0089] In specific operation, firstly, the source wavelet (i.e., the seismic wavelet, with a dominant frequency of 5Hz, the Ricker wavelet) is subjected to forward modeling with the tomographic inversion velocity field obtained in step S4 to obtain the forward modeling wavefield and forward modeling seismic data at each time step (the forward modeling wavefield and forward modeling seismic data obtained from the forward modeling simulation based on the tomographic inversion velocity field obtained in step S4); then, according to the objective functional of the full waveform inversion, the difference between the low-frequency seismic observation data obtained in step S21 and the forward modeling seismic data (the forward modeling seismic data obtained from the forward modeling simulation based on the tomographic inversion velocity field obtained in step S4) is calculated; then, the obtained difference is backpropagated with the wavefield, and then cross-correlated with the forward modeling wavefield (the forward modeling wavefield obtained from the forward modeling simulation based on the tomographic inversion velocity field obtained in step S4) to obtain the gradient field of the full waveform inversion; then, the gradient field is used to determine the appropriate iteration step size, and the velocity is updated to obtain the new velocity field.

[0090] S6, increase the data frequency band, repeat steps S2 to S5, and output the target's inverted velocity field.

[0091] In this step, the main frequency band of the seismic data is increased from 5Hz to 7Hz, and then sequentially to 9Hz and 10Hz to obtain new seismic data. Steps S2-S5 are then repeated to obtain the target velocity field. The results are as follows: Figure 6 As shown.

[0092] The specific operating 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 7Hz. Then, forward modeling is performed using the seismic wavelet (the dominant frequency of the Ricker wavelet is equal to 7Hz) and the velocity field obtained in step S5 to obtain the seismic forward modeling wavefield and forward modeling seismic data at each time step. Then, according to the objective functional of the full waveform inversion, the difference between the low-frequency seismic observation data with a frequency of 7Hz and the forward modeling seismic data is calculated. Then, the difference is backpropagated through the wavefield and cross-correlated with the seismic forward modeling wavefield to obtain the gradient field of the full waveform inversion. Then, the gradient field is used to determine the appropriate iteration step size for velocity update. 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 tomographic inversion using gathers to obtain the tomographic inversion velocity field. Finally, the tomographic inversion velocity field is used for full waveform inversion to obtain the velocity field.

[0094] ② The two-dimensional seismic observation data is then filtered to obtain low-frequency seismic observation data with a frequency of 9Hz. Then, forward modeling is performed using the seismic wavelet (the Ricker wavelet with a dominant frequency of 9Hz) and the velocity field obtained in step ① (the velocity field obtained from the low-frequency seismic observation data with a frequency of 7Hz) to obtain the seismic forward modeling wavefield and forward modeling seismic data at each time step. Then, based on the target functional of the full waveform inversion, the difference between the low-frequency seismic observation data with a frequency of 9Hz and the forward modeling seismic data is calculated. The difference is then backpropagated through the wavefield and cross-correlated with the seismic forward modeling wavefield to obtain the gradient field of the full waveform inversion. Then, the gradient field is used to determine the appropriate iteration step size for velocity update. The obtained velocity field is then fused with the well velocity field to obtain the well-constrained velocity field. Then, the well-constrained velocity field is used for grid tomographic inversion using gathers to obtain the tomographic inversion velocity field. Finally, the tomographic 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 performed using the seismic wavelet (the Ricker wavelet with a dominant frequency of 10Hz) and the velocity field obtained in step ② (the velocity field obtained from the low-frequency seismic observation data with a frequency of 9Hz). The forward modeling wavefield and forward modeling seismic data at each time moment are obtained. Then, according to the target functional of the full waveform inversion, the difference between the low-frequency seismic observation data with a frequency of 10Hz and the forward modeling seismic data is calculated. The difference is then backpropagated through the wavefield and cross-correlated with the forward modeling wavefield to obtain the gradient field of the full waveform inversion. The gradient field is then used to determine the appropriate iteration step size for velocity update. The obtained velocity field is then fused with the well velocity field to obtain the well-constrained velocity field. The well-constrained velocity field is then subjected to grid tomography inversion using gathers to obtain the tomography inverted velocity field. Finally, the tomography inverted velocity field is subjected to full waveform inversion to obtain the target velocity field.

[0096] In step S2, the Marmousi model ( Figure 5 In the inversion of ), the original velocity field model is used. Figure 3 After repeated iterations, the velocity field inversion results obtained are as follows: Figure 6 As shown, by Figure 6 It can be seen that the seismic velocity inversion method of the present invention has a better inversion effect on complex structural layers, can invert the velocity field more accurately, and eliminates the instability effects present in conventional inversion methods.

[0097] To compare the seismic velocity inversion method of this invention with conventional inversion methods, the original velocity field is processed using a conventional inversion method. Specifically, a grid tomography inversion is first performed using seismic observation data and the original velocity field to obtain the initial velocity field. Then, a full waveform inversion is performed (same as the full waveform inversion in step S2). The result is as follows: Figure 7 As shown. By Figure 7 It is known that conventional inversion methods have poor stability. Compared with conventional methods, the sum of the absolute errors between the velocity values ​​(predicted velocity values) and actual velocity values ​​obtained by using the method of this invention (the sum of absolute errors is equal to the sum of the absolute errors of the velocity values ​​of each test point, where the absolute error of the velocity at a 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 this invention is the same as that in the conventional method, 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, repeating operation steps S1 to S4 multiple times, with the frequency of low-frequency data increasing sequentially during repeated operations, and the initial velocity field used in a certain operation being the velocity field obtained after the previous operation ended. The formula for calculating the well-constrained velocity field is as follows: In the formula, For well-constrained velocity field, For the well velocity field, The velocity field obtained in step S1, The weighting coefficients for well constraints; 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, and h is the depth of the velocity field. The maximum frequency for full waveform inversion. The minimum frequency for full waveform inversion. represents the initial coefficients for well constraints.

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.

Images