A wave equation inversion method for constructing background velocity fields in deep and medium depths

Abstract

The present application relates to the field of exploration geophysics and discloses a wave equation inversion method for constructing a medium-deep background velocity field, including S1. inputting observation shot record data and an initial velocity model, S2. forward modeling the initial velocity field to calculate the forward wave field, and back-propagating the observation data to calculate the back wave field; cross-correlating the forward wave field and the back wave field according to imaging conditions to calculate the offset imaging result, S3. simulating a single reflection wave data based on the offset imaging result and the initial velocity field back-migration calculation, obtaining the data residual using the optimal transmission distance, constructing a companion source and substituting it into the reflection wave inversion, and calculating the gradient under the optimal transmission metric, S4. iteratively updating the background velocity field based on the calculated gradient, and performing full waveform inversion using the background velocity field information to obtain the final high-precision inversion result. By more accurately extracting the wave field residual based on the optimal transmission distance strategy, background field gradient information that is more consistent with the underground structure is constructed.

Description

Technical Field

[0001] The present invention relates to the technical field of exploration geophysics, and in particular to a wave equation inversion method for constructing a medium-deep background velocity field. Background Art

[0002] With the continuous deepening of exploration and development, oil and gas exploration is gradually tending towards deeper and more complex geological targets. How to improve the accuracy of seismic inversion imaging in complex geological structures at medium and deep depths has become a new challenge. The full waveform inversion method is a high-precision seismic exploration method based on a local optimization algorithm. It uses the full wavefield information of seismic waves to accurately model the velocity field representing the geological structure.

[0003] Full waveform inversion mainly relies on latent wave information for inversion, but the limited migration aperture and inaccurate initial velocity make it easy for full waveform inversion to fail to obtain the global optimal solution due to cycle skipping, making it difficult to accurately characterize medium-deep structural information. Although the reflection wave inversion method based on the gradient construction of the first reflection wave information provides good medium-deep background velocity information for full waveform inversion, the data residual measurement method of conventional reflection wave inversion based on the point-to-point least squares distance easily ignores the phase difference between wave fields and causes cycle skipping problems. For example, an initial velocity field with too large a difference in travel time information will make it impossible to correctly extract the wave field difference in the background velocity inversion, which will cause the inverted background velocity field to fall into a local extreme value, and it is impossible to obtain high-precision inversion results and migration images. Summary of the Invention

[0004] In response to the shortcomings of the existing technology, the present invention provides a wave equation inversion method for constructing a medium-deep background velocity field, which solves the problem that the existing technology relies on latent wave information for inversion, the limited offset aperture and the inaccurate initial velocity make it easy for full waveform inversion to fail to obtain a global optimal solution due to frequency skipping.

[0005] To achieve the above objectives, the present invention is implemented through the following technical solutions: a wave equation inversion method for constructing a mid-deep background velocity field, comprising:

[0006] S1. Input observation gun record data and initial velocity model;

[0007] S2. Calculate the forward wavefield by forward modeling the initial velocity field, and calculate the reverse wavefield by reverse propagation of the observation data. Cross-correlate the forward and reverse wavefields based on the imaging conditions to calculate the migration imaging results.

[0008] S3. Based on the migration imaging results and the initial velocity field demigration calculation, simulate the primary reflection wave data, use the optimal transmission distance to obtain the data residual, construct the companion source and substitute it into the reflection wave inversion, and calculate the gradient under the optimal transmission metric;

[0009] S4. Based on the calculated gradient, the background velocity field is iteratively updated, and full waveform inversion is performed using the background velocity field information to obtain the final high-precision inversion result.

[0010] Preferably, in step S1, the observation gun recording data and the initial velocity model are input, and the preset parameters include the number of model transverse sampling points n x , number of longitudinal sampling points n z , horizontal sampling interval d x , longitudinal sampling interval d z , number of time sampling points n t , time sampling interval d t , the main frequency of the earthquake source f.

[0011] Preferably, in step S2, the forward wavefield is calculated by combining preset parameters using the following formula:

[0012]

[0013] Where v is the velocity field, u s is the forward wave field, t is the time index, x and z are the horizontal and vertical grid indices respectively, and s is the source term;

[0014] The back propagation field is calculated using the following formula:

[0015]

[0016] Among them, u r is the back propagation wave field, d obs It is the observation gun record data;

[0017] Then, the forward wavefield and the reverse wavefield obtained above are cross-correlated using the following formula to calculate the migration imaging result:

[0018]

[0019] Where ref is the migration imaging result.

[0020] Preferably, in step S3, the simulated primary reflection wave data is de-migrated and calculated using the following formula:

[0021]

[0022] Where v(x) is the velocity, is the Laplace operator, u ref (x, t) is the simulated primary reflection wave field, u s (x, t) is the forward wave field of the initial velocity field forward modeling in claim 3; the simulated primary reflection wave field u is recorded according to the preset parameters. ref (x,t), denoted as Using the optimal transmission distance per channel, the objective function is calculated using the following formula:

[0023]

[0024] in, is the objective function under the optimal transmission distance, To observe the primary reflection wave data, are the probability distributions of observed and simulated primary reflection wave data, t is the time index of the observed data, and the accompanying source under the optimal transmission distance is calculated by the following formula:

[0025]

[0026] in, is the adjoint source, U is an upper triangular matrix with one non-zero element, diag is the diagonalization operation, and t′ is the time index of the corresponding simulation data when the optimal transmission metric is satisfied. The reflected wave inversion gradient under the optimal transmission distance is calculated by the following formula:

[0027]

[0028] Among them, g OTRWI is the gradient, u0, are the forward and reverse background wave fields, u p 、 They are the forward and reverse disturbance wave fields respectively.

[0029] Preferably, in S3, when using the optimal transmission distance to measure data, the data needs to be converted into a probability distribution form, and the data conversion is achieved by the following formula:

[0030]

[0031] in, Represents x s Cannon x r Earthquake data from the channel, is the regularized seismic data to satisfy the non-negativity, For the corresponding x s Cannon x r Probability distribution of seismic data.

[0032] Preferably, in said S4, when constructing the background velocity field gradient, the forward and reverse background wave fields u0, Solved by the following formula:

[0033]

[0034] Forward and reverse disturbance wave field u p and Solved by the following formula:

[0035]

[0036] Substituting the two calculated background wave fields and two disturbance wave fields into the gradient formula can obtain the result of the background gradient.

[0037] Preferably, in S4, the background velocity field is iteratively updated by the following formula:

[0038] v k+1 =v k +αg OTRWI k

[0039] Among them, k is the iteration index, v is the updated background velocity field, α is the optimal iteration step size, g OTRWI is the background velocity field gradient.

[0040] Preferably, the conjugate gradient method is used to optimize the gradient, and the parabola fitting method is used to calculate the optimal iterative step size.

[0041] Preferably, the calculated final background velocity field inversion result is substituted into a conventional full waveform inversion to calculate the final high-precision velocity profile.

[0042] Preferably, the objective function and the accompanying source are measured according to the optimal transmission distance, and a channel-by-channel calculation is implemented in a 2D case.

[0043] The present invention provides a wave equation inversion method for constructing a background velocity field in a medium-deep layer. It has the following beneficial effects:

[0044] 1. The present invention more accurately extracts wave field residuals based on the optimal transmission distance strategy and constructs background field gradient information that is more consistent with the underground structure. This can avoid the false updates on the wave path of traditional methods, making the final inverted background velocity field more consistent with the structural characteristics, thereby providing a good initial velocity field for full waveform inversion, greatly reducing the frequency jump problem of full waveform inversion, improving the accuracy and stability of inversion, and further improving the inversion resolution and accuracy of medium-deep geological targets, ultimately improving the inversion quality in complex geological conditions at medium and deep depths.

[0045] 2. The present invention adopts the conjugate gradient method to optimize the gradient during the inversion process and uses the parabola fitting method to calculate the optimal iterative step size. These optimization methods can significantly improve the computational efficiency and the convergence speed of the inversion process. Especially when processing large-scale, high-resolution seismic data, these methods can effectively reduce computing time and resource consumption, and improve the overall efficiency of the inversion.

[0046] 3. The present invention uses the optimal transmission distance to measure the difference between observed data and simulated data, constructs a new objective function that focuses more on data wavefront matching, and can more accurately evaluate the similarities and differences of data, thereby more accurately updating the velocity model during the inversion process. Ultimately, the present invention can calculate high-resolution velocity field inversion results, especially in the inversion of medium-deep geological targets, which can significantly improve the inversion accuracy and resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 A schematic flow chart of a wave equation inversion method for constructing a mid-deep background velocity field according to an embodiment of the present invention;

[0048] Figure 2 It is the true velocity model diagram of the present invention;

[0049] Figure 3 It is the initial velocity model diagram of the present invention;

[0050] Figure 4 This is a diagram showing the calculation results of reverse time migration imaging of the initial velocity model of the present invention;

[0051] Figure 5 This is a diagram showing the calculation results of the companion source under the optimal transmission distance of the present invention;

[0052] Figure 6 This is a graph showing the gradient calculation results of the reflection wave inversion for the optimal transmission distance of the present invention;

[0053] Figure 7 This is a diagram showing the calculation results of the background velocity field for the optimal transmission distance reflection wave inversion of the present invention;

[0054] Figure 8 This is a diagram showing the final velocity field calculation results of the full waveform inversion using the inverted background velocity field as the initial model. DETAILED DESCRIPTION

[0055] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the present specification. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0056] Example:

[0057] Please see the attached Figure 1 -Attached Figure 8 The embodiment of the present invention provides a wave equation inversion method for constructing a mid-deep background velocity field, comprising the following steps:

[0058] S1, input observation gun record data and initial velocity model;

[0059] The preset parameters include the number of horizontal sampling points n of the model x , number of longitudinal sampling points n z , horizontal sampling interval d x , longitudinal sampling interval d z , number of time sampling points n t , time sampling interval d t , the main frequency f of the earthquake source, what needs to be prepared is the actual observed seismic shot recording data and an initial velocity model for starting the inversion process. These data and models are generated in the forward calculation and compared with the real velocity model when the final result is confirmed to verify the effect of the inversion. Parameter settings, including the number of horizontal and vertical sampling points and the sampling interval, are the basis for building these models and ensure the spatial resolution of the model. The number of time sampling points and the time sampling interval affect the time resolution of the model. The choice of the main frequency of the earthquake source is also very important for the accuracy of the simulation. It determines the frequency range of the simulated waveform. In the embodiment of the present application, the observation shot recording data is obtained by forward calculation of the real velocity field model. The real velocity field is only used to generate the observation shot recording data and compare the results. It does not participate in the inversion calculation. The input real velocity field and background velocity field are as follows. Figure 1 As shown, the model has 670 horizontal sampling points, 260 vertical sampling points, and both the horizontal and vertical sampling intervals are 10 m.

[0060] S2: Forward modeling of the initial velocity field calculates the forward wavefield, and inverse modeling of the observed data calculates the inverse wavefield. Forward modeling is based on the initial velocity model, simulating the propagation of seismic waves in the underground medium by solving the wave equation to calculate the forward wavefield. Inverse modeling is based on the actual observed seismic data, reversing the wave equation to calculate the inverse wavefield. The forward and inverse wavefields are then cross-correlated to calculate the migration imaging result. Based on the imaging conditions, the forward and inverse wavefields are cross-correlated to calculate the migration imaging result.

[0061] Combined with the preset parameters, the forward wave field is calculated using formula (1):

[0062]

[0063] Where v is the velocity field, u s is the forward wave field, t is the time index, x and z are the horizontal and vertical grid indices respectively, and s is the source term;

[0064] The back propagation field is calculated by formula (2):

[0065]

[0066] Where u ris the back propagation wave field, d obs It is the observation gun record data;

[0067] Then, use formula (3) to perform cross-correlation on the forward wavefield and the reverse wavefield obtained above to calculate the migration imaging result:

[0068]

[0069] Where ref is the migration imaging result, and the calculated migration imaging result is as follows: Figure 4 shown.

[0070] S3, based on the migration imaging results and the initial velocity field, simulates primary reflection wave data through back-migration calculation. The optimal transmission distance is used to obtain data residuals, and a companion source is constructed and substituted into the reflection wave inversion to calculate the gradient under the optimal transmission metric. The migration imaging results are used together with the initial velocity model for back-migration calculation to simulate primary reflection wave data. The optimal transmission distance is used to measure the residual between the observed data and the simulated data. A new objective function is constructed. This objective function focuses on data wavefront matching rather than traditional waveform matching, which helps to avoid spurious updates on the wave path. By using the companion source under the optimal transmission distance to calculate the gradient, it can more accurately guide the update of the velocity model and improve the stability of the inversion.

[0071] The simulated primary reflection wave data is calculated by de-migration using formula (4):

[0072]

[0073] Where v(x) is the velocity, is the Laplace operator, u ref (x, t) is the simulated primary reflection wave field, u s (x, t) is the forward wave field of the initial velocity field in formula (1); the simulated primary reflection wave field u is recorded according to the preset parameters. ref (x,t), denoted as Using the optimal transmission distance per channel, the objective function is calculated using formula (5):

[0074]

[0075] Where, is the objective function under the optimal transmission distance, To observe the primary reflection wave data, are the probability distributions of observed and simulated primary reflection wave data, t is the time index of the observed data, and the companion source under the optimal transmission distance is calculated by formula (6):

[0076]

[0077] Where, is the adjoint source, U is an upper triangular matrix with one non-zero element, diag is a diagonalization operation, t′ is the time index of the corresponding simulation data when the optimal transmission metric is satisfied, that is, the probability distribution of the simulation data corresponding to the index at time t′ is equal to the probability distribution of the observation data corresponding to the index at time t. The binary search method is used here. Here, the objective function and the adjoint source are measured according to the optimal transmission distance. The calculation is performed channel by channel in the 2D case. The adjoint source calculation results are as follows: Figure 5 As shown, the reflected wave inversion gradient under the optimal transmission distance is calculated by formula (7):

[0078]

[0079] Where g OTRWI is the gradient, u0, are the forward and reverse background wave fields, u p 、 They are the forward and reverse disturbance wave fields respectively.

[0080] When using the optimal transmission distance to measure data, the data needs to be converted into a probability distribution form, and the data conversion is achieved through formulas (8) and (9):

[0081]

[0082] Where, Represents x s Cannon x r Earthquake data from the channel, is the regularized seismic data to satisfy the non-negativity, For the corresponding x s Cannon x r Probability distribution of seismic data.

[0083] Among them, when formula (7) constructs the background velocity field gradient, the forward and reverse background wave fields u0, Solve by formula (10) and (11):

[0084]

[0085]

[0086] Forward and reverse disturbance wave field u p and Solve by formula (12) and (13):

[0087]

[0088] Substituting the calculated two background wave fields and two disturbance wave fields into the gradient formula (7) can obtain the background gradient result. The calculation result of the background gradient is as follows: Figure 6 shown.

[0089] S4. Finally, the calculated gradient is used to update the background velocity field through an iterative algorithm (such as the conjugate gradient method). In this process, the optimal iterative step size strategy (such as parabola fitting) is adopted to ensure that the updated velocity, direction, and amplitude can minimize the objective function. The updated background velocity field can provide a more accurate starting point for the subsequent full waveform inversion, thereby obtaining a higher-precision underground medium velocity distribution image. Specifically, according to the calculated gradient, the background velocity field is iteratively updated, and the background velocity field information is used for full waveform inversion to obtain the final high-precision inversion result.

[0090] The background velocity field is iteratively updated using formula (14):

[0091] v k+1 =v k +αg OTRWI k (14)

[0092] Where k is the iteration index, v is the updated background velocity field, α is the optimal iteration step size, g OTRWI is the background velocity field gradient. During iteration, the conjugate gradient method is used to optimize the gradient, and the parabola fitting method is used to calculate the optimal iteration step size. The background velocity field finally calculated is as follows: Figure 7 shown.

[0093] Finally, the calculated final background velocity field inversion result is substituted into the conventional full waveform inversion to calculate the final high-precision velocity profile, such as Figure 8 The final results show that the method of the present invention can calculate relatively accurate high-resolution velocity field inversion results.

[0094] The method of the present invention can avoid the false updates on the wave path of traditional methods, make the background velocity field of the final inversion more consistent with the structural characteristics, thereby providing a good initial velocity field for full waveform inversion, greatly reducing the cycle jump problem of full waveform inversion, improving the accuracy and stability of inversion, and further improving the inversion resolution and accuracy of medium-deep geological targets, and ultimately improving the inversion quality in complex geological conditions of medium-high depths.

[0095] The full waveform inversion imaging performed in the above embodiment is 2D, which can be expanded to 3D with the following changes: the number of transverse sampling points, the number of longitudinal sampling points, the transverse sampling interval, and the longitudinal sampling interval are changed to parameters in three directions, for example, 670 transverse sampling points and 260 longitudinal sampling points, and the number of vertical sampling points is added, for example, 300. The transverse sampling interval is 10 m, the longitudinal sampling interval is 10 m, and the vertical sampling interval is 10 m. The calculation of the forward wavefield and the backward wavefield needs to be expanded to three dimensions, and the cross-correlation operation needs to be performed in three-dimensional space. The calculation of data residuals and accompanying sources needs to be performed on a channel-by-channel basis in three-dimensional space.

[0096] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A wave equation inversion method for constructing the background velocity field in the deep sea. comprising, characterized in that, S1. Input observation gun record data and initial velocity model; S2. Calculate the forward wavefield by forward modeling the initial velocity field, and calculate the reverse wavefield by reverse propagation of the observation data. Cross-correlate the forward and reverse wavefields based on the imaging conditions to calculate the migration imaging results. S3. Calculate and simulate primary reflection wave data based on the migration imaging results and the initial velocity field de-migration. Obtain data residuals using the optimal transmission distance. Construct a companion source and substitute it into the reflection wave inversion. Calculate the gradient under the optimal transmission metric. Specifically, de-migrate and calculate the simulated primary reflection wave data using the following formula: Where v(x) is the velocity, is the Laplace operator, u ref (x, t) is the simulated primary reflection wave field, u s (x, t) is the forward wave field of the initial velocity field; the simulated primary reflection wave field u is recorded according to the preset parameters. ref (x,t), denoted as Using the optimal transmission distance per channel, the objective function is calculated using the following formula: in, is the objective function under the optimal transmission distance, To observe the primary reflection wave data, are the probability distributions of observed and simulated primary reflection wave data, t is the time index, x is the horizontal grid index, and the accompanying source under the optimal transmission distance is calculated by the following formula: in, is the adjoint source, U is an upper triangular matrix with one non-zero element, diag is the diagonalization operation, and t′ is the time index of the corresponding simulation data when the optimal transmission metric is satisfied. The reflected wave inversion gradient under the optimal transmission distance is calculated by the following formula: Among them, g OTRWI is the gradient, u0, are the forward and reverse background wave fields, u p 、 They are the forward and reverse disturbance wave fields respectively; S4. Based on the calculated gradient, the background velocity field is iteratively updated, and full waveform inversion is performed using the background velocity field information to obtain the final high-precision inversion result.

2. A wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In S1, the observation gun record data and the initial velocity model are input, including: Observation gun recorded data, initial velocity model and preset parameters, the preset parameters include the number of model lateral sampling points n x , number of longitudinal sampling points n z , horizontal sampling interval d x , longitudinal sampling interval d z , number of time sampling points n t , time sampling interval d t , the main frequency of the earthquake source f.

3. The wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In S2, the forward wave field is calculated by combining the preset parameters using the following formula: Among them, v is the background velocity field, u s is the forward wave field, t is the time index, x and z are the horizontal and vertical grid indices respectively, and s is the source term; The back propagation field is calculated using the following formula: Among them, u r is the back propagation wave field, d obs It is the observation gun record data; Then, the forward wavefield and the reverse wavefield obtained above are cross-correlated using the following formula to calculate the migration imaging result: Where ref is the migration imaging result.

4. The wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In S3, when using the optimal transmission distance to measure data, the data needs to be converted into a probability distribution form, and the data conversion is achieved through the following formula: in, Represents x s Cannon x r Earthquake data from the channel, is the regularized seismic data to satisfy the non-negativity, For the corresponding x s Cannon x r Probability distribution of seismic data.

5. The wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In the above S4, when constructing the background velocity field gradient, the forward and reverse background wave fields u0, Solved by the following formula: Forward disturbance wave field u p and back-propagation disturbance wave field Solved by the following formula: Substituting the two calculated background wave fields and two disturbance wave fields into the gradient formula can obtain the result of the background gradient.

6. The wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In S4, the background velocity field is iteratively updated by the following formula: in k+1 =in k +αg OTRWI k Among them, k is the iteration index, v is the updated background velocity field, α is the optimal iteration step size, g OTRWI is the background velocity field gradient.

7. The method according to any one of claims 1 to 6, characterized in that: The conjugate gradient method is used to optimize the gradient, and the parabola fitting method is used to calculate the optimal iteration step size.

8. The method according to any one of claims 1 to 6, characterized in that: The calculated final background velocity field inversion result is substituted into the conventional full waveform inversion to calculate the final high-precision velocity profile.

9. The method according to claim 6, characterized in that According to the optimal transmission distance measurement objective function and the accompanying source, the channel-by-channel calculation is implemented in the 2D case.

Images