Well three-dimensional earthquake full-waveform inversion velocity modeling method

Abstract

The invention provides an in-well three-dimensional seismic full-waveform inversion velocity modeling method, which comprises the following steps of: 1, acquiring an in-well excitation three-dimensional seismic near-well initial velocity model, and reading seismic data; step 2, calculating a simulated seismic record; step 3, performing data regularization; 4, selecting a target function; step 5, calculating an accompanying seismic source; and step 6, calculating an update gradient to update the speed model. According to the in-well three-dimensional seismic full-waveform inversion velocity modeling method, the characteristics of in-well excitation three-dimensional seismic data are fully considered, wave components in the data are fully utilized, the problem that different wave components are separated and inversed respectively in a traditional method is avoided, and a full-wave field inversion method for in-well excitation three-dimensional seismic data is realized.

Description

Technical Field

[0001] The invention relates to the fields of engineering seismic exploration and medical imaging, and in particular to a method for velocity modeling of three-dimensional seismic full waveform inversion in a well. Background Art

[0002] As oil exploration and development moves toward areas with complex geological conditions, the difficulty and risk of oil drilling continue to increase, placing higher demands on improving the success rate of oil drilling and development and reducing drilling risks. Safety, accuracy, and efficiency have become the ultimate goals of oil drilling.

[0003] The seismic profiles produced by borehole induced 3D seismic are similar to those produced by conventional seismic profiling (VSP), but borehole induced 3D seismic data are collected more quickly and economically. In conventional VSP, the seismic source is placed on the surface and several geophones are placed in the borehole. In borehole induced 3D seismic, the source is placed in the borehole and the survey line is fixed on the surface. The borehole induced 3D seismic data are similar to VSP data obtained by placing the source at many different locations on the surface. The difficulty in using borehole induced 3D seismic technology is how to design a downhole source with enough energy to pass through the formation to the surface without damaging the wellbore. Compared with surface seismic, the advantages of both VSP and borehole induced 3D seismic are that their energy travels a short distance when passing through the formation and passes through the loose near-surface low-velocity zone only once, so the signal attenuation is small, especially the high-frequency components of the data. Borehole induced 3D seismic data can provide physical parameters such as velocity, deconvolution factor, attenuation coefficient, seismic wavelet, reflection coefficient, etc. for other seismic data processing, and can also help surface seismic identify multiple waves and indicate the source and propagation process of multiple waves. At the same time, the information it provides, such as velocity and velocity ratio, wavelength and wavelength ratio, frequency and frequency ratio, amplitude and amplitude ratio, and waveform structure, can be used to study concealed reservoirs and sandstone bodies around the wellbore, detect oil and gas, and delineate the scope of oil storage.

[0004] At present, oil and gas development has gradually shifted from conventional oil and gas reservoirs to complex oil and gas reservoirs. The corresponding requirements for geophysical exploration have also gradually increased. The refined imaging of oil and gas structures has become the main trend of geophysical exploration. The development of high-precision imaging and inversion methods that meet the corresponding requirements has become an urgent problem to be solved. Reverse Time Migration (RTM) and Full Waveform Inversion (FWI) are currently the main methods for imaging underground fine geological structures. Although certain practical application effects have been achieved in offshore data and actual land seismic data, there are still many theoretical and practical challenges. The most significant problem is the high degree of dependence on the initial model. When the initial model and the true model are very different, FWI is prone to frequency jump problems, resulting in inversion instability, and RTM imaging results will also be distorted. The establishment of an accurate initial model is the key to the success of FWI and RTM. The corresponding velocity modeling method, especially the research on low-frequency model construction methods, has become a hot topic in the current research of seismic exploration methods. The 3D seismic data from borehole excitation contains relatively abundant transmission and reflection wave components, and compared with ground seismic, the signal-to-noise ratio is higher, which is very beneficial for high-precision and high-resolution imaging near the well. FWI, as a medium- and high-resolution imaging method commonly used in geophysics, is very suitable for processing 3D seismic data from borehole excitation. However, due to insufficient illumination of the observation system, there are currently few studies on the application of FWI to 3D seismic data from borehole excitation.

[0005] In the Chinese patent application with application number: CN201310133893.2, a frequency domain full waveform inversion seismic velocity modeling method is involved, which includes the following steps: 1) obtaining the original seismic shot record, source wavelet information and the initial model used for inversion; 2) analyzing the information obtained in step 1), and determining the basic inversion parameters and the full waveform inversion framework from low frequency to high frequency based on the forward algorithm and the optimization algorithm; 3) calculating and obtaining the most suitable frequency domain forward and inversion model grid for different frequencies; 4) using the principal component analysis method to compress the data dimension involved in the inversion during low-frequency inversion; 5) judging whether the projection matrix dimension corresponding to different frequencies meets the threshold conversion standard, and when the conversion standard is met, proceeding to the next step, otherwise returning to step 4); 6) introducing the source encoding method, and suppressing the crosstalk noise by using the random phase encoding method; 7) judging whether the iteration cutoff condition is met, and if the iteration cutoff condition is met, proceeding to the next step; otherwise returning to step 6); 8) if the inversion of all frequencies is not completed, returning to step 3), until the inversion of all frequencies is completed, the final velocity model is obtained, and the velocity model is output.

[0006] In the Chinese patent application with application number: CN202010141175.X, a full-waveform velocity modeling inversion method based on geological structure model constraints is involved. The invention inputs seismic data inversion parameters, initial velocity model, actual observed shot records, and structural interpretation data, performs forward simulation according to the initial velocity model to obtain forward shot records and forward wave fields, constructs geological structure inversion operators using geological structure models, constrains the regularization terms of conventional FWI inversion models, obtains full-waveform velocity inversion models based on geological structure model constraints, obtains the residuals between all forward shot records and their corresponding actual observed shot records and back-propagates, calculates the single shot gradient and total velocity gradient of each shot record in combination with the forward wave field, sets the test step size and the number of cycles, updates the step size using the actual observed shot records, completes the iterative update of the velocity field, and obtains the velocity inversion result. The invention uses the geological model as a constraint condition, effectively enhances the inversion effect, improves the velocity inversion accuracy, and is of great significance in the exploration of complex fault block areas.

[0007] In the Chinese patent application with application number: CN201510295941.7, a method for full waveform inversion of seismic waves based on a least squares gradient updated velocity model is involved, which includes the following steps: 1) obtaining an acoustic wave equation of a seismic wave field in the time domain; 2) constructing an initial velocity model, setting the number of velocity model update iterations N and the minimum allowable error value ε; 3) constructing a wave field error vector of observed wave field data and calculated wave field data; 4) constructing an objective function; 5) calculating the objective function to obtain a standard equation for full waveform inversion of seismic waves; 6) introducing a velocity model update gradient direction gk and an update step length α; 7) solving the velocity model update gradient direction gk using the least squares method; 8) interpolating the update step length α; 9) updating the velocity model to obtain: mk=mk-1+αgk; when |αgk|<ε or the number of velocity model updates reaches the number of velocity model update iterations N, the velocity model update ends; otherwise, proceeding to step 3). The invention can quickly complete the velocity model update and is widely used in the full waveform inversion method of seismic waves.

[0008] In the Chinese patent application with application number: CN201910367245.0, a full waveform inversion method based on the fast quasi-Newton method is involved, which belongs to the technical field of seismic exploration velocity modeling, and aims to improve the convergence speed of the quasi-Newton method, thereby improving the computational efficiency of seismic full waveform inversion and reducing the amount of calculation. The method includes: defining a set of weighted coefficient sequences for constructing a parameter model iteration formula; introducing an intermediate iteration parameter model, calculating the quasi-Newton iteration direction and optimal step size based on the parameter model; using the weighted coefficient, iteration direction and optimal step size to update the original parameter model and the intermediate iteration parameter model. The invention calculates a parameter model that is closer to the optimal solution, and is used to calculate the next iteration direction and optimal step size, thereby improving the convergence speed of the objective function. In each iterative calculation, the inventive method does not increase the amount of calculation too much, but only a small amount of multiplication and addition operations with low complexity, which effectively reduces the amount of calculation of full waveform inversion.

[0009] The above existing technologies are all significantly different from the present invention and fail to solve the technical problem we want to solve. For this reason, we have invented a new in-well 3D seismic full waveform inversion velocity modeling method. Summary of the invention

[0010] The purpose of the present invention is to provide a method for in-well three-dimensional seismic full waveform inversion velocity modeling using full wave field information for velocity model inversion.

[0011] The purpose of the present invention can be achieved by the following technical measures: a three-dimensional seismic full waveform inversion velocity modeling method in a well, the three-dimensional seismic full waveform inversion velocity modeling method in a well comprising:

[0012] Step 1, obtaining an initial velocity model of the wellside of a three-dimensional seismic excitation in the well and reading seismic data;

[0013] Step 2, calculating simulated earthquake records;

[0014] Step 3: regularize the data;

[0015] Step 4, select the objective function;

[0016] Step 5, calculating the associated earthquake source;

[0017] Step 6: Calculate the update gradient to update the velocity model.

[0018] The purpose of the present invention can also be achieved by the following technical measures:

[0019] In step 1, the initial velocity model is set to a relatively simple gradient model or a relatively smooth velocity model obtained by picking up actual data, where the actual data is obtained by a tomographic method.

[0020] In step 1, during actual data processing, the observation data is obtained through actual field collection. At this time, file reading is required, and the general data format is SEGY data body.

[0021] In step 1, the SEGY format includes three parts, namely, an EBCDIC file header, a binary file header and a seismic track. Each seismic track contains track header information and seismic track data. The track header information stores the line number, track number, number of sampling points, geodetic coordinates, shot point coordinates and detection point coordinates corresponding to the seismic track.

[0022] In step 2, the constructed initial velocity model is used to perform acoustic wave finite difference forward modeling to obtain simulated seismic records. The acoustic wave equation used here is a three-dimensional constant density acoustic wave equation:

[0023]

[0024] Where u represents the seismic wave field, v is the propagation speed of the seismic wave, (x, y, z) constitutes a two-dimensional Cartesian coordinate system, and t is time; the above equation is discretized using finite differences, and the wave field u is obtained through approximate calculation.

[0025] In step 3, the oscillating seismic data are converted into probability distribution to satisfy the non-negative and equal conditions of the optimal transport theory. The observed data and the simulated records need to be processed separately for data regularization. Here, a linear weighted method is used.

[0026] In step 3, the linear weighting formula is as follows:

[0027]

[0028]

[0029] where d obs is the observation record, d cal is the simulated data, a and b are linear weighting operators, g and f are the probability distribution functions of the regularized observed data and simulated data, respectively, to ensure the non-negative and equal conditions of the optimal transport theory.

[0030] In step 4, considering that the probability distribution functions f and g are functions of x, according to the optimal transportation theory, the objective function of the two distribution functions f and g can be constructed as follows:

[0031]

[0032] Where f and g are the probability distribution functions calculated from the simulated data and observed data of x, respectively. F and G are the cumulative functions of the probability distribution functions f and g, and their expressions are:

[0033] F(x)=∫f(x)dx

[0034] G(x)=∫g(x)dx, (4)

[0035] Where M represents the transportation strategy that maps distribution f to distribution g.

[0036] The in-well 3D seismic full waveform inversion velocity modeling method further includes, after step 4, judging whether the data residual meets the requirements, setting a threshold, and outputting the velocity inversion result when the residual calculated by formula 3 is less than the threshold; when the residual is greater than the threshold, executing step 5.

[0037] In step 5, the corresponding companion source can be derived by the conjugate state method as follows:

[0038]

[0039] in

[0040]

[0041] Where δs represents the conjugate source, L is an upper triangular matrix, and diag is a diagonal matrix; for each y value, the interpolation search method can be used to find the x value corresponding to the monotonically increasing function G(x).

[0042] In step 6, based on the constructed initial model and the obtained conjugate source calculation gradient, the velocity model is continuously updated through an iterative algorithm until the data error meets the requirements.

[0043] The core of the in-well 3D seismic full waveform inversion velocity modeling method of the present invention is the joint inversion of direct waves and reflected waves. The traditional method is to separate the direct waves and reflected waves and invert them separately. Since it is difficult to separate the wave components of in-well excited 3D seismic data, in order to make full use of the wave components in the in-well excited 3D seismic data and at the same time ensure the stability of the inversion to the greatest extent, the present invention has developed a full-wavefield full-waveform inversion method for in-well excited 3D seismic data based on the optimal transmission distance theory, and performed high-precision and high-resolution imaging near the well.

[0044] This method can make full use of the wave components in the seismic data and effectively weaken the frequency skipping problem that occurs in the full waveform inversion when the initial model is poor. A more effective initial velocity model is established, and the high-frequency components of the model are gradually inverted. This method can be applied to medical imaging. Since the three-dimensional seismic data excited in the well contains rich transmission wave components, the problem that the full waveform inversion is more sensitive to the initial velocity model can be effectively weakened. At the same time, the present invention introduces the optimal transport theory to weaken the dependence of the full waveform inversion on the initial model, which can ensure the stability of the inversion process. The present invention fully considers the characteristics of the three-dimensional seismic data excited in the well, makes full use of the wave components in the data, avoids the problem of separating different wave components and inverting them separately in the traditional method, and realizes a full wave field inversion method for the three-dimensional seismic data excited in the well. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a specific implementation flow chart of the method for modeling the full waveform velocity of three-dimensional seismic in-well according to the present invention;

[0046] Figure 2 It is a schematic diagram of a true velocity model in a specific embodiment of the present invention;

[0047] Figure 3 It is a schematic diagram of an initial velocity model in a specific embodiment of the present invention;

[0048] Figure 4 In a specific embodiment of the present invention, Figure 2 Schematic diagram of the slice at the position y = 1.0 km;

[0049] Figure 5 In a specific embodiment of the present invention, Figure 3 Schematic diagram of the slice at the position y = 1.0 km;

[0050] Figure 6 In a specific embodiment of the present invention, Figure 2 Schematic diagram of the final inversion result at the position of y = 1.0 km;

[0051] Figure 7 In a specific embodiment of the present invention, Figure 3 Comparison of velocity curves of the true velocity model, initial velocity model and final inversion result at the position x=1.5km.

[0052] Figure 8 It is a schematic diagram of a true velocity model in a specific embodiment of the present invention;

[0053] Fig. 9 It is a schematic diagram of an initial velocity model in a specific embodiment of the present invention;

[0054] Fig.10In a specific embodiment of the present invention, Figure 2 Schematic diagram of the slice at the position y=2.2km;

[0055] Fig.11 In a specific embodiment of the present invention, Figure 3 Schematic diagram of the slice at the position y=2.2km;

[0056] Fig.12 In a specific embodiment of the present invention, Figure 2 Schematic diagram of the final inversion result at the position of y = 2.2 km;

[0057] Fig.13 In a specific embodiment of the present invention, Figure 3 Comparison of velocity curves of the true velocity model, initial velocity model and final inversion result at x=1.7km;

[0058] Fig.14 It is a schematic diagram of a true velocity model in a specific embodiment of the present invention;

[0059] Fig.15 It is a schematic diagram of an initial velocity model in a specific embodiment of the present invention;

[0060] Fig.16 In a specific embodiment of the present invention, Figure 2 Schematic diagram of the slice at the position of y=2.25km;

[0061] Fig.17 In a specific embodiment of the present invention, Figure 3 Schematic diagram of the slice at the position of y=2.25km;

[0062] Fig.18 In a specific embodiment of the present invention, Figure 2 Schematic diagram of the final inversion result at the position of y = 2.25 km;

[0063] Fig.19 In a specific embodiment of the present invention, Figure 3 Comparison of velocity curves of the true velocity model, initial velocity model and final inversion result at the position x=1.8km. DETAILED DESCRIPTION

[0064] In order to make the above and other purposes, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings.

[0065] like Figure 1 As shown, Figure 1The following is a flow chart of the three-dimensional full waveform inversion velocity modeling method in a well based on the optimal transportation distance of the present invention. The basic process is as follows:

[0066] Step 101: Initial model construction

[0067] The initial model is generally set to a relatively simple gradient model or a relatively smooth velocity model obtained by picking up actual data. The actual data is often obtained through tomography methods.

[0068] Step 102, read observation data

[0069] In actual data processing, the observation data is obtained through actual field collection. At this time, file reading is required, and the general data format is SEGY data body. The SEGY format consists of three parts, namely EBCDIC file header, binary file header and seismic trace. Each seismic trace contains trace header information and seismic trace data; the trace header information stores the line number, trace number, number of sampling points, geodetic coordinates, shot point coordinates and detection point coordinates corresponding to the seismic trace, and the data is read through the program.

[0070] Step 103, simulate earthquake record calculation

[0071] The initial velocity model constructed above is used to perform acoustic wave finite difference forward modeling to obtain simulated seismic records. The acoustic wave equation used here is a three-dimensional constant density acoustic wave equation:

[0072]

[0073] Where u represents the seismic wave field, v is the propagation velocity of the seismic wave, (x, y, z) constitutes a two-dimensional Cartesian coordinate system, and t is time. The above equation is discretized using finite differences, and the wave field u is obtained through approximate calculation.

[0074] Step 104: Data regularization

[0075] To convert the oscillating seismic data into a probability distribution so that it meets the non-negative and equal conditions of the optimal transport theory, it is necessary to perform data regularization on the observed data and the simulated records respectively. Here, the linear weighted method is used. The specific formula is as follows:

[0076]

[0077] where d obs is the observation record, d cal is the simulated data, a and b are linear weighting operators, g and f are the probability distribution functions of the regularized observed data and simulated data, respectively, to ensure the non-negative and equal conditions of the optimal transport theory.

[0078] Step 105: Selection of objective function

[0079] There are many kinds of objective functions based on optimal transport theory (Engquist and Froese, 2014; Yang et al., 2017; Métivier et al., 2018). Here we mainly discuss the 2-Wasserstein distance. Considering that the probability distribution functions f and g are functions of x, according to the optimal transport theory, the objective function of the two distribution functions f and g can be constructed as follows:

[0080]

[0081] Where f and g are the probability distribution functions calculated from the simulated data and observed data of x, respectively. F and G are the cumulative functions of the probability distribution functions f and g, and their expressions are:

[0082]

[0083] The distribution satisfies M represents a transportation strategy that maps distribution f to distribution g.

[0084] Step 106: Determine whether the data residual meets the requirements

[0085] A threshold is set, and the velocity inversion result is output when the residual calculated by Formula 3 is less than the threshold. When the residual is greater than the threshold, step 107 is executed.

[0086] Step 107, calculate the associated earthquake source

[0087] The corresponding companion source can be derived by the conjugate state method as follows:

[0088]

[0089] in

[0090]

[0091] Where δs represents the conjugate source and L is an upper triangular matrix. For each y value, the interpolation search method can be used to find the x value corresponding to the monotonically increasing function G(x).

[0092] Step 108, calculate the update gradient to update the velocity model

[0093] Based on the constructed initial model and the obtained conjugate source calculation gradient, the velocity model is continuously updated through an iterative algorithm until the data error meets the requirements.

[0094] like Figure 1In the process shown, the steps for establishing a high-precision velocity model are: obtain the initial velocity model of the wellside of the three-dimensional seismic excited in the well through the steps shown in step 101; read the seismic data using the method shown in step 102; obtain the seismic record using the method shown in step 103; regularize the data using the method shown in step 104 to ensure the non-negative and equal conditions of the optimal transport theory; obtain the objective function using the method shown in step 105; determine whether the data residual meets the requirements using the method shown in step 106; calculate the accompanying earthquake source using the method shown in step 107; and update the model using the calculation update gradient shown in step 108.

[0095] The present invention does not need to separate the wave components during the calculation process, and the observation and simulation records are input as a whole, converted into distributions, and the distances between the distributions are obtained.

[0096] The following are several specific embodiments of the present invention.

[0097] Example 1

[0098] In a specific embodiment 1 of the present invention, the algorithm is tested based on an actual geological model, which is a three-dimensional thin layer model ( Figure 2 ), the model size is 3000*320*242, 7 wells are set in the model (indicated by asterisks in the figure), the shot spacing is 50m, 10 shots per well, a total of 77,440 (320*242) detectors (full coverage), and the trace spacing is 10m. Initial model ( Figure 3 ) is set as a constant gradient model, where the first layer velocity is the same as the true velocity. The full wave field data is obtained by acoustic forward modeling using the true model. The final inversion result can be obtained by inputting the full wave field record and the initial model for inversion. For the convenience of comparison, the true velocity model at the position x = 1 km is given here ( Figure 4 ), initial model( Figure 5 ) and the inversion results ( Figure 6 ) slice, we can see that the inversion result is close to the real model, and the extracted velocity curve ( Figure 7 ) It can also be seen that this method can effectively restore underground structures and has good stability.

[0099] Example 2

[0100] In a specific embodiment 2 of the present invention, the algorithm is tested based on an actual geological model, which is a three-dimensional thin layer model ( Figure 8 ), the model size is 201*338*676, 2 wells are set in the model (indicated by asterisks in the figure), the shot spacing is 10m, 50 shots in each well, a total of 228488 (338*676) detectors (full coverage), and the trace spacing is 10m. Initial model ( Fig. 9 ) is set as a constant gradient model. The full wave field data is obtained by forward acoustic modeling using the true model. The full wave field record and the initial model are input for inversion to obtain the final inversion result. For the convenience of comparison, the true velocity model at x = 1.7 km is given here ( Fig.10 )、initial model( Fig.11 ) and the inversion results ( Fig.12 ) slice, we can see that the inversion result is close to the real model, and the extracted velocity curve ( Fig.13 ) It can also be seen that this method can effectively restore underground structures and has good stability.

[0101] Example 3

[0102] In a specific embodiment 3 of the present invention, the algorithm is tested based on an actual geological model, which is a three-dimensional thin layer model ( Fig.14 ), the model size is 187*401*401, 2 wells are set in the model (indicated by asterisks in the figure), the shot spacing is 10m, 50 shots in each well, a total of 160801 (401*401) detectors (full coverage), and the trace spacing is 10m. Initial model ( Fig.15 ) is set as a constant gradient model. The full wave field data is obtained by forward acoustic modeling using the true model. The final inversion result can be obtained by inputting the full wave field record and the initial model for inversion. For the convenience of comparison, the true velocity model at the position x = 1.8 km is given here ( Fig.16 )、initial model( Fig.17 ) and the inversion results ( Fig.18 ) slice, we can see that the inversion result is close to the real model, and the extracted velocity curve ( Fig.19 ) It can also be seen that this method can effectively restore underground structures and has good stability.

[0103] The method for inverting the direct wave and reflected wave full wave field velocity based on the three-dimensional seismic data excited in the well in the present invention fully takes into account the characteristics that the three-dimensional seismic data excited in the well has a high signal-to-noise ratio and the wave components of the transmitted wave and the reflected wave are difficult to separate. At the same time, the degree of dependence on the initial model is reduced, and the stability of the full waveform inversion is greatly improved. In the process of inversion, only the observation record and the initial velocity model need to be input, and iterative inversion can be performed to obtain the final high-resolution inversion result. Its objective function is the distance based on the 2-Wasserstein metric. Before inversion, the seismic data needs to be converted into a distributed form, and then the corresponding conjugate source is calculated according to the formula, and the conjugate source is back-transmitted by the wave field, and cross-correlated with the forward wave field, and the corresponding gradient can be obtained. Then, the velocity model is updated iteratively by the conjugate gradient or quasi-Newton method, and the final inversion result is obtained. Although the method does not require high accuracy of the initial model, a better initial velocity model and constraint conditions are also conducive to the convergence speed of the global optimization method. Although the formulas involved in the present invention are all two-dimensional, they can also be applied to three-dimensional seismic data processing.

[0104] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0105] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A method for velocity modeling of 3D seismic full waveform inversion in a well, characterized in that: The well 3D seismic full waveform inversion velocity modeling method includes: Step 1, obtaining an initial velocity model of the wellside of a three-dimensional seismic excitation in the well and reading seismic data; Step 2, calculating simulated earthquake records; Step 3: regularize the data; Step 4, select the objective function; Step 5, calculating the associated earthquake source; Step 6: Calculate the update gradient to update the velocity model.

2. The method for three-dimensional seismic full waveform inversion velocity modeling in a well according to claim 1, characterized in that: In step 1, the initial velocity model is set to a relatively simple gradient model or a relatively smooth velocity model obtained by picking up actual data, where the actual data is obtained by a tomographic method.

3. The method for three-dimensional seismic full waveform inversion velocity modeling in a well according to claim 2, characterized in that: In step 1, during actual data processing, the observation data is obtained through actual field collection. At this time, file reading is required, and the general data format is SEGY data body.

4. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 3, characterized in that: In step 1, the SEGY format includes three parts, namely, the EBCDIC file header, the binary file header and the seismic channel. Each seismic channel contains the channel header information and the seismic channel data. The channel header information stores the line number, channel number, number of sampling points, geodetic coordinates, shot point coordinates and detection point coordinates corresponding to the seismic channel.

5. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 1, characterized in that: In step 2, the constructed initial velocity model is used to perform acoustic wave finite difference forward modeling to obtain simulated seismic records. The acoustic wave equation used here is a three-dimensional constant density acoustic wave equation: Where u represents the seismic wave field, v is the propagation speed of the seismic wave, (x, y, z) constitutes a two-dimensional Cartesian coordinate system, and t is time; the above equation is discretized using finite differences, and the wave field u is obtained through approximate calculation.

6. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 1, characterized in that: In step 3, the oscillating seismic data are converted into probability distribution to satisfy the non-negative and equal conditions of the optimal transport theory. The observed data and the simulated records need to be processed separately for data regularization. Here, a linear weighted method is used.

7. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 6, characterized in that: In step 3, the linear weighting formula is as follows: where d obs is the observation record, d cal is the simulated data, a and b are linear weighting operators, g and f are the probability distribution functions of the regularized observed data and simulated data, respectively, to ensure the non-negative and equal conditions of the optimal transport theory.

8. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 7, characterized in that: In step 4, considering that the probability distribution functions f and g are functions of x, according to the optimal transportation theory, the objective function of the two distribution functions f and g can be constructed as follows: Where f and g are the probability distribution functions calculated from the simulated data and observed data of x, respectively. F and G are the cumulative functions of the probability distribution functions f and g, and their expressions are: F(x)=∫f(x)dx G(x)=∫g(x)dx, (4) The distribution satisfies M represents a transportation strategy that maps distribution f to distribution g.

9. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 8, characterized in that: The in-well 3D seismic full waveform inversion velocity modeling method further includes, after step 4, judging whether the data residual meets the requirements, setting a threshold, and outputting the velocity inversion result when the residual calculated by formula 3 is less than the threshold; when the residual is greater than the threshold, executing step 5.

10. The method for in-well 3D seismic full waveform inversion velocity modeling according to claim 9, characterized in that: In step 5, the corresponding companion source can be derived by the conjugate state method as follows: in Where δs represents the conjugate source, L is an upper triangular matrix, and diag is a diagonal matrix; for each y value, the interpolation search method can be used to find the x value corresponding to the monotonically increasing function G(x).

11. The method for three-dimensional seismic full waveform inversion velocity modeling in a well according to claim 1, characterized in that: In step 6, based on the constructed initial model and the obtained conjugate source calculation gradient, the velocity model is continuously updated through an iterative algorithm until the data error meets the requirements.

Images