Wave equation inversion method for constructing middle-deep part background velocity field

Abstract

The invention relates to the field of exploration geophysics, and discloses a wave equation inversion method for constructing a middle-deep part background velocity field, which comprises the steps of S1, inputting observation gun record data and an initial velocity model, S2, performing forward modeling on the initial velocity field to calculate a forward propagation wave field, and performing back propagation on observation data to calculate a back propagation wave field; s3, calculating and simulating primary reflected wave data according to the migration imaging result and the initial velocity field back migration, solving a data residual error by using an optimal transmission distance, constructing an adjoint source, substituting the adjoint source into reflected wave inversion, calculating a gradient under the optimal transmission measure, and calculating a gradient under the optimal transmission measure; and S4, iteratively updating the background velocity field according to the calculated gradient, and performing full-waveform inversion by using the background velocity field information to obtain a final high-precision inversion result. Wave field residual errors are extracted more accurately based on an optimal transmission distance strategy, and background field gradient information more conforming to an underground structure is constructed.

Description

Technical Field

[0001] The 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 has gradually tended 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 that characterizes the geological structure.

[0003] Full waveform inversion mainly relies on latent wave information for inversion, but the limited offset 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 mid-deep structural information. Although the reflection wave inversion method that constructs a gradient based on the information of a single reflection wave provides good mid-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 for the background velocity inversion to correctly extract the wave field difference, which will cause the inverted background velocity field to fall into a local extreme value, and no high-precision inversion result and offset image can be obtained. Summary of the invention

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

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

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

[0007] S2. Forward modeling of the initial velocity field calculates the forward wave field, and reverse modeling of the observation data calculates the reverse wave field; cross-correlating the forward wave field and the reverse wave field according to the imaging conditions, and calculating the migration imaging result;

[0008] S3. Simulate the reflection wave data based on the migration imaging results and the initial velocity field back-migration calculation, 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. According to the calculated gradient, the background velocity field is iteratively updated, and the background velocity field information is used to perform full waveform inversion 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 lateral sampling points n of the model. x , the number of longitudinal sampling points n z , horizontal sampling interval d x , longitudinal sampling interval d z , the 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 wave field 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 wave field is calculated by the following formula:

[0015]

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

[0017] Then, the forward wave field and the reverse wave field 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 first 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 by 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, respectively, t is the time index of the observed data, and the companion source under the optimal transmission distance is calculated by the following formula:

[0025]

[0026] in, is the companion source, U is an upper triangular matrix with one non-zero element, diag is a 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 calculated two background wave fields and two disturbance wave fields into the gradient formula, the background gradient result can be obtained.

[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 iteration step size.

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

[0042] Preferably, the objective function and the companion 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 medium-deep background velocity field. It has the following beneficial effects:

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

[0045] 2. The present invention uses the conjugate gradient method to optimize the gradient during the inversion process and uses the parabola fitting method to calculate the optimal iteration 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 measures the difference between observed data and simulated data by using the optimal transmission distance, 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 accuracy and resolution of the inversion. 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 of the reverse time migration imaging calculation results 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 of the optimal transmission distance reflection wave inversion of the present invention;

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

[0055] The following will be combined with the drawings in the specification of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0056] Example:

[0057] Please refer to 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 lateral sampling points n of the model x , the 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 are the actual observed seismic shot record 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 setting, including the number of horizontal and vertical sampling points and the sampling interval, is the basis for building these models and ensures 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 record data is obtained by forward calculation of the real velocity field model. The real velocity field is only used to generate the observation shot record 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 the horizontal and vertical sampling intervals are both 10 m.

[0060] S2, forward calculation of the initial velocity field, reverse calculation of the reverse wave field by observation data; forward calculation is based on the initial velocity model, and the propagation process of seismic waves in the underground medium is simulated by solving the wave equation to calculate the forward wave field, while inversion calculation is based on the actual observed seismic data, and the reverse wave field is calculated by inversely analyzing the wave equation. Then, the result of the migration imaging can be calculated by cross-correlating the forward wave field and the reverse wave field, and the forward wave field and the reverse wave field are cross-correlated according to the imaging conditions 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] In the formula, u ris the anti-transmission wave field, d obs It is the observation gun record data;

[0067] Then, the forward wave field and the reverse wave field obtained above are cross-correlated using formula (3) to calculate the migration imaging result:

[0068]

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

[0070] S3, simulate the reflection wave data according to the migration imaging result and the initial velocity field anti-migration calculation, use the optimal transmission distance to obtain the data residual, construct the companion source and substitute it into the reflection wave inversion, calculate the gradient under the optimal transmission metric, the migration imaging result and the initial velocity model are used for anti-migration calculation, so as to simulate the reflection wave data, use the optimal transmission distance to measure the residual between the observed data and the simulated data, and construct a new objective function, which pays more attention to the matching of the data wavefront rather than the traditional waveform matching, which helps to avoid false updates on the wave path. By using the companion source under the optimal transmission distance to calculate the gradient, the update of the velocity model can be more accurately guided, and the stability of the inversion can be improved;

[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 first 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 by formula (5):

[0074]

[0075] In the formula, 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, respectively, 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] In the formula, is the companion source, U is an upper triangular matrix with one non-zero element, diag is a diagonalization operation, t′ is the time index corresponding to the simulated data when the optimal transmission metric is satisfied, that is, the probability distribution of the simulated data corresponding to the index at time t′ is equal to the probability distribution of the observed data corresponding to the index at time t. The binary search is used here. Here, according to the optimal transmission distance metric objective function and the companion source, the calculation is implemented in the 2D case. The calculation results of the companion source are as follows: Figure 5 As shown, the reflected wave inversion gradient under the optimal transmission distance is calculated by formula (7):

[0078]

[0079] In the formula, 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] In the formula, Represents x s Cannon x r The earthquake data of 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 iteration step 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 by 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, and g is 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. 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 As shown in the final result, it can be seen 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 update of the traditional method on the wave path, 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 weakening 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 under 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. Specifically, there are 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 into parameters in three directions respectively, for example, the number of transverse sampling points is 670, the number of longitudinal sampling points is 260, the number of vertical sampling points is added: for example, 300, the transverse sampling interval is 10m, the longitudinal sampling interval is 10m, and the vertical sampling interval is: for example, 10m. The calculation of the forward wave field and the reverse wave field 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 track by track in three-dimensional space.

[0096] Although 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 the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present 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 and medium depths, comprising, characterized in that, S1. Input observation gun record data and initial velocity model; S2. Forward modeling of the initial velocity field calculates the forward wave field, and reverse modeling of the observation data calculates the reverse wave field; cross-correlating the forward wave field and the reverse wave field according to the imaging conditions, and calculating the migration imaging result; S3. Simulate the reflection wave data based on the migration imaging results and the initial velocity field back-migration calculation, 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; S4. According to the calculated gradient, the background velocity field is iteratively updated, and the background velocity field information is used to perform full waveform inversion 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 , the 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. A 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: 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; The back propagation wave field is calculated by the following formula: Among them, u r is the anti-transmission wave field, d obs It is the observation gun record data; Then, the forward wave field and the reverse wave field obtained above are cross-correlated using the following formula to calculate the migration imaging result: Where ref is the migration imaging result.

4. A wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In S3, the simulated primary reflection wave data is de-migrated by the following formula: Where v(x) is the velocity, is the Laplace operator, u ref (x, t) is the simulated first reflection wave field, u s (x, t) is the forward wave field of the initial velocity field; the simulated 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 by 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, respectively, t is the time index of the observed data, and the companion source under the optimal transmission distance is calculated by the following formula: in, is the companion source, U is an upper triangular matrix with one non-zero element, diag is a 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.

5. A 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.

6. A wave equation inversion method for constructing a mid-deep background velocity field according to claim 1, characterized in that: In 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-propagating disturbance wave field Solved by the following formula: Substituting the calculated two background wave fields and two disturbance wave fields into the gradient formula, the background gradient result can be obtained.

7. A 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.

8. The method according to any one of claims 1 to 7, 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.

9. The method according to any one of claims 1 to 7, 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.

10. The method according to claim 1, 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