Joint inversion method for velocity and Q value of seismic viscous sound waves in well

Abstract

The invention provides a well seismic viscous sound wave speed and Q value joint inversion method, which comprises the following steps of: 1, forward modeling a real speed and Q model to obtain an observation record, and constructing an initial speed and Q model; 2, performing forward modeling by utilizing a viscous sound wave equation to obtain a simulation observation record; 3, comparing the simulation record with the observation record, solving a data residual error based on the optimal transportation distance, and calculating a corresponding accompanying seismic source; 4, carrying out back propagation on the residual error by utilizing a conjugate state equation to obtain a corresponding gradient, and updating the model; and step 5, carrying out loop iteration on the process, exiting the loop when an iteration condition is met, and outputting a corresponding inversion result. According to the in-well seismic viscous sound wave speed and Q value joint inversion method, a speed and Q two-parameter inversion full-waveform inversion method which aims at in-well three-dimensional seismic data and takes an optimal transportation distance as a distance measurement method is developed.

Description

Technical Field

[0001] The invention relates to the fields of engineering seismic exploration and medical imaging, and in particular to a joint inversion method of well seismic viscous acoustic wave velocity and Q value. Background Art

[0002] Seismic exploration is currently the most commonly used and effective exploration method in the field of oil and gas resource exploration. With the continuous development of oil and gas exploration, complex areas and deep oil and gas exploration have gradually been valued, so higher requirements have been put forward for seismic technology, especially the rapid high-resolution imaging method of deep complex geological structures has received more attention. The reverse time migration imaging method based on the two-way wave theory can effectively overcome the steep angle limitation of traditional methods, and can accurately locate seismic waves of any propagation path. It is more suitable for high-precision imaging of complex structures such as high-steep-angle rock masses and fault block depressions.

[0003] Absorption attenuation is an inherent property of the medium, which has an important impact on the propagation of seismic waves, causing the energy of the seismic wave field to attenuate and the velocity to disperse. Energy attenuation is manifested as a decrease in the amplitude of seismic waves and a weakening of the reflection signal of the deep medium, resulting in poor imaging quality and low resolution of deep strata; velocity dispersion will cause changes in the travel of seismic waves and distortion of the phase, resulting in deviations in the formation interface and structural position on the imaging profile, affecting the interpretation of seismic data. On the other hand, when the reservoir contains gas, its absorption attenuation effect is enhanced. If the absorption attenuation effect is not compensated, it may lead to the difficulty of imaging the reservoir boundary and the failure of the reservoir prediction method AVO (Amplitude versus offset) inversion, which will have an adverse effect on the interpretation of structure and lithology. Therefore, it is very important to study the propagation law and response characteristics of seismic waves in viscous acoustic wave media and to compensate and correct the seismic data accordingly. Generally, traditional compensation methods directly act on seismic records, while attenuation compensation based on wave equations fully considers the propagation path of seismic waves, and the compensation is more in line with the real physical principles. Therefore, the development of high-resolution Q-compensated imaging methods for viscous acoustic wave media has become a hot topic in current research.

[0004] In actual geophysical exploration for oil and gas, the absorption and attenuation of the near-surface has a serious impact on seismic waves, so the accurate reconstruction of the Q and velocity models near the surface is of great significance for high-resolution and high-precision imaging of deep geological structures. The conventional seismic profiling method (VSP) is an effective means of high-precision imaging of the shallow surface due to the short propagation distance of seismic waves and less signal attenuation. Due to the observation method of ground blasting and receiving in the well, the lighting is relatively limited and the construction restrictions are large, resulting in high construction costs. The seismic profiles produced by reverse vertical seismic profiling (RVSP) are similar to those produced by conventional seismic profiling (VSP), but RVSP data are collected faster and more economically. The advantages of implementing RVSP observations are mainly the following: setting up detectors on the ground is less affected by surface conditions and the relationship between workers and farmers, and a large number of detectors can be deployed, which can minimize the loss of facet data volume, improve the quality of profile imaging, and form a data volume similar to ground seismic; the number of RVSP excitation points is much smaller than that of three-dimensional VSP excitation points, which improves production timeliness, reduces production costs, and reduces the impact of uneven energy generated by different excitation points on data; the number of excitation points in the well is less affected by construction conditions, and the number of excitation points is greatly increased compared to the number of VSP detectors, so the number of coverage times is greatly enhanced; ground-based detectors can use area combinations to eliminate the impact of ground noise to the maximum extent. Based on the above analysis of the advantages of RVSP, it can be seen that the development of RVSP technology has relatively important practical significance.

[0005] The full waveform inversion method is considered to be the most effective method for high-resolution imaging. This method is highly dependent on the initial model and is prone to cycle skipping, causing the inversion to fall into a local minimum. In order to weaken the dependence of full waveform inversion on the initial model, scholars have developed corresponding time domain, frequency domain and Laplace domain methods, and have achieved good application results. In recent years, the full waveform inversion method based on the optimal transmission distance objective function proposed by the Seismic Scope research group can effectively avoid cycle skipping. The robustness of its inversion algorithm has been confirmed by the inversion of synthetic data and actual data. The process of deriving the gradient of this inversion method is more complicated than that of the traditional method, and the calculation of the optimal transmission distance objective function itself is a global optimization process. Compared with the traditional full waveform inversion method, its calculation amount is significantly increased. An important development direction of full waveform inversion is multi-parameter inversion. With the development of numerical simulation technology of wave equations in complex media, multi-parameter full waveform inversion methods have become possible. In the multi-parameter full waveform inversion method, the simultaneous inversion of velocity and Q is a very difficult problem, which stems from two reasons. First, as mentioned above, most of the existing time-domain viscoelastic wave equations do not conform to the constant Q model commonly used in the field of seismic exploration, so they cannot accurately describe the attenuation and dispersion of seismic waves. Secondly, velocity and Q scatter seismic waves in the same way, that is, the seismic response caused by velocity disturbances and Q disturbances is similar. This strong coupling effect causes the update of velocity and Q to interfere with each other during the inversion process, making it difficult to obtain ideal inversion results.

[0006] The wave field components of 3D seismic data in wells are mainly transmission waves and reflection waves. In order to weaken the influence of frequency jumps, traditional full waveform usually separates transmission waves and reflection waves and inverts them separately. For 3D seismic data in wells, wave field separation is more difficult.

[0007] In the Chinese patent application with application number: CN201710994606.5, a method for acoustic impedance inversion with Q absorption attenuation compensation is involved. Impedance inversion is one of the important methods for seismic oil and gas resource exploration. The present invention integrates the formation Q absorption attenuation dispersion factor into the convolution model, and directly inverts the high-resolution impedance results from the non-steady-state seismic recording data of absorption attenuation through an iterative inversion method. The inventive method avoids the error accumulation effect existing in the traditional recursive process and the numerical instability problem existing in the traditional anti-Q filtering method. That is, the Q value compensation technology of the invention is unconditionally stable, and it can compensate for higher frequency components than traditional technologies, further improving the resolution of seismic exploration data.

[0008] In the Chinese patent application with application number: CN201910367448.X, a full waveform inversion method based on a preconditioned conjugate gradient acceleration algorithm is involved, which belongs to the field of seismic exploration inversion imaging technology. It aims to address the defects of slow convergence speed and low deep inversion accuracy of conjugate gradient full waveform inversion methods. By combining the advantages of gradient preprocessing operators and Nesterov accelerated gradient descent methods, the deep inversion accuracy of conjugate gradient full waveform inversion methods is improved while the convergence speed is effectively improved, thereby reducing the computational complexity of full waveform inversion, and providing technical support for the practical application of full waveform inversion.

[0009] In the Chinese patent application with application number: CN201610020556.6, a method for inverting seismic reflection coefficients based on total variation minimization constraints is involved. The post-stack seismic data is inverted by each time window. For a single time window, the post-stack seismic data and the pre-extracted seismic wavelet in the time window are first Fourier transformed, and then the frequency domain expression of the reflection coefficient is obtained. The time domain reflection coefficient is then Fourier transformed, and the real part and the imaginary part are extracted on the basis of the parity decomposition of the reflection coefficient. The corresponding solution equation is constructed, and the minimum total variation constraint is adopted on the basis of the traditional conjugate gradient algorithm to solve the equation, and the odd and even reflection coefficients are obtained and the original reflection coefficient in the time window is reconstructed; all time windows are inverted in turn to obtain the reflection coefficient of a single seismic record. All seismic traces are inverted in turn to obtain the reflection coefficient of each seismic record. The method of the invention can effectively invert the seismic reflection coefficient, which is beneficial to improve the resolution of seismic data and improve the accuracy of reservoir prediction.

[0010] In the Chinese patent application with application number: CN202011178090.5, a method for forward modeling of acoustic seismic data in viscous media is involved. By separating the amplitude attenuation term from the phase dispersion term, the invention can conveniently realize viscoacoustic reverse time migration imaging compared to the traditional forward modeling of the viscoacoustic wave equation in which the attenuation term and the dispersion term are coupled together. Compared with traditional calculation methods, it can greatly improve the calculation efficiency and reduce the calculation memory.

[0011] The above existing technologies are greatly different from the present invention and fail to solve the technical problem we want to solve. Therefore, we have invented a new joint inversion method for well seismic viscous acoustic wave velocity and Q value. Summary of the invention

[0012] The purpose of the present invention is to provide a viscous acoustic full wave field velocity and Q model joint inversion method based on optimal transport distance, which can effectively weaken the frequency jump problem that occurs in full waveform inversion when the velocity and Q initial models are poor, and enhance the stability of inversion.

[0013] The object of the present invention can be achieved by the following technical measures: a joint inversion method of well seismic viscous acoustic wave velocity and Q value, the joint inversion method of well seismic viscous acoustic wave velocity and Q value comprises:

[0014] Step 1, forward model the real velocity and Q model to obtain the observation record and construct the initial velocity and Q model;

[0015] Step 2, using the viscous acoustic wave equation to perform forward modeling to obtain simulated observation records;

[0016] Step 3, compare the simulated records with the observed records, obtain the data residual based on the optimal transportation distance, and calculate the corresponding accompanying earthquake source;

[0017] Step 4: Back propagate the remaining residual using the conjugate state equation to obtain the corresponding gradient and update the model;

[0018] Step 5, iterate the above process, exit the loop when the iteration condition is met, and output the corresponding inversion result.

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

[0020] In step 1, the initial velocity and Q models are relatively simple gradient models or relatively smooth velocity models obtained by picking up actual data, and the actual data are obtained by tomography.

[0021] 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 data format is SEGY data body.

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

[0023] In step 2, the initial velocity and Q model are used for forward modeling to obtain the corresponding observation records and simulation records. The observation system adopts the in-well excitation seismic observation method, that is, the seismic observation is excited in the well and the data is received on the ground.

[0024] In step 2, for the simulation data, the fractional-order constant-Q viscous acoustic wave equation is used for forward modeling. The specific formula is as follows:

[0025]

[0026] Where u is the wave field, c0 is the wave propagation speed in the medium, Q is the quality factor, λ=(ω d / ω0) 2γ ,ωd =2πf d is the main frequency of the source wavelet, ω0 is the reference angular frequency, the first two terms of the equation mainly control the phase distortion, and the third term controls the amplitude attenuation; the finite difference algorithm is used to differentiate the above equation, where the left side of the equation is the second-order derivative of the wave field with respect to time, which can be approximated by the second-order central difference operator, and the first-order partial derivative of time on the right side is approximated by the backward difference operator, and the fractional-order Laplace operator is calculated by FFT. When the above wave equation is used for numerical simulation, three FFT and three IFFT operations are required in each time step.

[0027] The method for joint inversion of well seismic viscous acoustic wave velocity and Q value also includes, after step 2, regularizing the data and selecting an objective function.

[0028] When regularizing the data, the oscillating seismic data is converted into a probability distribution to satisfy 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, using a linear weighted method. The specific formula is as follows:

[0029] g=a*d obs +b

[0030] f=a*d cal +b. (2)

[0031] where d obs and d cal are the observed data and simulated data, respectively, 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 transportation theory

[0032] When selecting the objective function, the source δs is calculated by the partial derivative of the objective function to the data. The objective function here is defined according to the optimal transportation distance. The objective function formula of the optimal transportation theory is J OT as follows

[0033]

[0034] 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:

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

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

[0037] The distribution satisfies f: g: M represents the transportation strategy that maps the probability distribution function f to the probability distribution function g. In step 3, it is determined whether the data residual meets the requirements, and a threshold is set. When the residual calculated by formula 3 is less than the threshold, the inversion results of velocity and Q are output; when the residual is greater than the threshold, the associated earthquake source is calculated.

[0038] In step 3, when calculating the accompanying earthquake source, the conjugate state method is used to perform back propagation on the accompanying earthquake source to obtain an updated velocity model. According to formula (3), the expression of the accompanying earthquake source can be derived as follows:

[0039]

[0040] Where δs represents the associated earthquake source, L is an upper triangular matrix, and diag is a diagonal matrix; F and G are cumulative functions of the probability distribution functions f and g; for each y value, the interpolation search method is used to find the x value corresponding to the monotonically increasing function G(x); for It can be further expressed as:

[0041]

[0042] In step 4, the gradients of the velocity model and the Q model are calculated. Based on the above wave equation, the non-attenuated adjoint equation and gradient expression of the full waveform inversion can be derived using the conjugate state method as follows:

[0043]

[0044]

[0045] Where c0 and Q are the reference speeds; f0 and f d are frequency and reference frequency respectively, and (-▽ 2 ) 0.505 is a constant fractional Laplace operator; dt is the time difference; u is the calculated wave field; ω d is the angular frequency; λ is the accompanying wave field; μ can be expressed as:

[0046]

[0047] The adjoint wave field λ satisfies the following wave equation:

[0048]

[0049] Where δs is the associated earthquake source, They can be expressed as:

[0050]

[0051]

[0052]

[0053] Where ω0 is the reference angular frequency.

[0054] In order to avoid the instability of the growing adjoint operator, a non-decaying adjoint operator is used. By removing the last term in the equation, we can obtain:

[0055]

[0056] The companion wave field obtained by this equation can maintain the same phase, but the amplitude is different.

[0057] In step 4, the velocity model and Q model are updated, and the nonlinear conjugate gradient algorithm is selected. The iterative update formula is: k+1 =m k +α k Δm k , where m is the model parameter, α is the iteration step, and Δm is the gradient obtained in step 108; for the inversion strategy problem in multi-parameter inversion, a simultaneous inversion strategy is adopted, that is, the gradients of the velocity and Q models are calculated simultaneously, and the velocity and Q models are updated at the same time.

[0058] The present invention firstly obtains the observation record by forward modeling with the real velocity and Q model, constructs the initial velocity and Q model, and obtains the simulated observation record by forward modeling with the viscous acoustic wave equation; compares the simulated record with the observation record, obtains the data residual based on the optimal transport distance, and calculates the corresponding accompanying earthquake source; back-propagates the remaining residual by using the conjugate state equation, obtains the corresponding gradient, and updates the model; iterates the above process in a loop, exits the loop when the iteration condition is met, and outputs the corresponding inversion result. The present invention makes full use of the effective wave component of the three-dimensional seismic in the well, provides high-resolution imaging results around the well location, and develops a full waveform inversion method for the three-dimensional seismic data in the well, which uses the optimal transport distance as the distance measurement method for the velocity and Q dual parameter inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a flowchart of a specific embodiment of the method for joint inversion of viscous acoustic wave velocity and Q value in a well based on the optimal transportation distance of the present invention;

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

[0061] Figure 3 It is a schematic diagram of a real Q model in a specific embodiment 1 of the present invention;

[0062] Figure 4It is a schematic diagram of the initial velocity model in the specific embodiment 1 of the present invention;

[0063] Figure 5 It is a schematic diagram of the initial Q model in the specific embodiment 1 of the present invention;

[0064] Figure 6 It is a schematic diagram of seismic data corresponding to the real model in the first specific embodiment of the present invention;

[0065] Figure 7 It is a schematic diagram of seismic data corresponding to the initial model in the first specific embodiment of the present invention;

[0066] Figure 8 It is a schematic diagram of the velocity model inverted in the specific embodiment 1 of the present invention;

[0067] Fig. 9 It is a schematic diagram of the Q model inverted in the specific embodiment 1 of the present invention;

[0068] Fig.10 It is a schematic diagram of seismic data obtained by forward calculation using an inversion model in a specific embodiment 1 of the present invention;

[0069] Fig.11 The corresponding extracted in the specific embodiment 1 of the present invention Fig.10 Schematic diagram of the comparison curve of seismic data calculated by the initial model, true model and inversion model at the red box in the middle;

[0070] Fig.12 It is a schematic diagram of comparison curves of the initial Q model, the real Q model and the inverted Q model at the position x=1.5 km extracted in the first specific embodiment of the present invention;

[0071] Fig.13 It is a schematic diagram of comparison curves of the initial velocity model, the true velocity model and the inversion velocity model at the position x=1.5 km extracted in the first specific embodiment of the present invention;

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

[0073] Fig.15 It is a schematic diagram of a real Q model in the second specific embodiment of the present invention;

[0074] Fig.16 It is a schematic diagram of the initial velocity model in the second specific embodiment of the present invention;

[0075] Fig.17 It is a schematic diagram of the initial Q model in the second specific embodiment of the present invention;

[0076] Fig.18It is a schematic diagram of the velocity model inverted in the specific embodiment 2 of the present invention;

[0077] Fig.19 It is a schematic diagram of the Q model inverted in the specific embodiment 2 of the present invention;

[0078] Fig. 20 It is a schematic diagram of comparison curves of the initial Q model, the real Q model and the inverted Q model at the position x=3.3 km in a specific embodiment of the present invention;

[0079] Fig.21 It is a schematic diagram of comparison curves of the initial velocity model, the true velocity model and the inversion velocity model at the position x=3.3 km in a specific embodiment of the present invention. DETAILED DESCRIPTION

[0080] In order to make the above and other purposes, features and advantages of the present invention more obvious and easy to understand, the following is a preferred embodiment, and a detailed description is given in conjunction with the accompanying drawings. This patent relates to the velocity and Q model inversion of three-dimensional viscous acoustic waves. For the convenience of description, a two-dimensional example is given here. The three-dimensional implementation method is the same as the two-dimensional one, except that the amount of calculation increases due to the addition of one dimension.

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

[0082] Example 1

[0083] like Figure 1 As shown, Figure 1 This is a flow chart of the viscous acoustic wave velocity and Q joint inversion method based on the optimal transport distance of the present invention. The basic process is as follows:

[0084] Step 101: Initial model construction

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

[0086] Step 102, read observation data

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

[0088] Step 103, obtaining simulation data

[0089] The initial velocity and Q model are used for forward modeling (Formula 1) to obtain the corresponding observation records and simulation records. Here, the observation system adopts the well-excited seismic observation method, that is, the well is excited and the ground receives data. For the simulation data, the fractional-order constant Q viscous acoustic wave equation is used for forward modeling. The specific formula is as follows:

[0090]

[0091] Where u is the wave field, c0 is the wave propagation speed in the medium, Q is the quality factor, λ=(ω d / ω0) 2γ ,ω d =2πf d is the main frequency of the source wavelet, ω0 is the reference angular frequency, the first two terms of the equation mainly control the phase distortion, and the third term controls the amplitude attenuation. The above equation can be derived using the finite difference algorithm, where the left side of the equation is the second-order derivative of the wave field with respect to time, which can be approximated using the second-order central difference operator, and the first-order partial derivative of time on the right side is approximated using the backward difference operator, and the fractional Laplace operator is calculated using FFT. When using the above wave equation for numerical simulation, three FFT and three IFFT operations are required in each time step.

[0092] Step 104: Data regularization

[0093] 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:

[0094]

[0095] Step 105: Selection of objective function

[0096] The associated earthquake source δs can be calculated by the partial derivative of the objective function with respect to the data. The objective function here is defined according to the optimal transportation distance, and the formula is as follows:

[0097]

[0098] 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:

[0099]

[0100] The distribution satisfies f: g: M represents the transportation strategy that maps the probability distribution function f to the probability distribution function g.

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

[0102] A threshold is set, and when the residual calculated by formula 3 is less than the threshold, the inversion results of velocity and Q are output. When the residual is greater than the threshold, step 107 is executed.

[0103] Step 107, calculate the associated earthquake source

[0104] By using the conjugate state method, the accompanying source is back-propagated to obtain an updated velocity model. According to formula (3), the accompanying source expression can be derived as follows:

[0105]

[0106] in

[0107]

[0108] Where δs represents the associated 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).

[0109] Step 108, calculate the gradients of the velocity model and the Q model

[0110] According to the above wave equation, using the conjugate state method, the non-attenuated adjoint equation and gradient expression of full waveform inversion can be derived as follows:

[0111]

[0112]

[0113] The adjoint wave field λ satisfies the following wave equation:

[0114]

[0115] in

[0116]

[0117]

[0118]

[0119] In order to avoid numerical instability, some scholars directly use the forward equation as the adjoint operator, which means that the adjoint wave field is attenuated, which will reduce the resolution of the gradient. This adjoint operator is called the attenuated adjoint operator. In order to avoid the instability of the growing adjoint operator, a non-attenuated adjoint operator is used here. By removing the last term in the equation, we can get:

[0120]

[0121] Where δs is the companion source, and the companion wavefield obtained can maintain the same phase, but the amplitude is different. Compared with the attenuated adjoint operator, the gradient constructed by the non-attenuated operator has higher resolution because the companion wavefield does not have amplitude attenuation. In addition, since it contains fewer Laplace operator terms, the non-attenuated adjoint operator also improves computational efficiency.

[0122] Step 109, update the velocity model and Q model

[0123] There are two key issues involved here: the selection of iterative algorithms and the strategy of multi-parameter inversion. The nonlinear conjugate gradient algorithm is used here. For general inversion iterative problems, the nonlinear conjugate gradient algorithm can already meet the corresponding requirements. The iterative update formula is: k+1 =m k +α k Δm k , where m is the model parameter, α is the iteration step, and Δm is the gradient obtained in step 108. For the inversion strategy problem in multi-parameter inversion, a simultaneous inversion strategy is adopted here, that is, the gradients of the velocity and Q models are calculated simultaneously, and the velocity and Q models are updated at the same time.

[0124] Example 2

[0125] In the specific embodiment 2 of the present invention, the real speed and Q model are as follows: Figure 1 As shown in the figure, the initial model is generally set to a relatively simple gradient model or a relatively smooth velocity model obtained by picking up actual data, respectively using the real ( Figure 2 ) and the initial velocity ( Figure 4 )、Q model( Figure 3 , Figure 5 ) is forward modeled (Formula 1) to obtain the corresponding observation records and simulation records. The observation system here adopts the well-excitation three-dimensional seismic observation method, that is, excitation in the well and receiving data on the ground. The well is located at x = 1km, the shot spacing is 10m, there are 30 shots in total, and a total of 320 geophones are laid on the surface (full coverage), the track spacing is 10m, and a 20Hz Ricker wavelet is used as the source. The observation and simulation records of the 10th shot are as follows. Figure 6 , Figure 7As shown, the final inverted velocity and Q model are Figure 8 , Fig. 9 As shown. It can be seen that the inversion result is close to the real result, and the simulated record obtained by re-forward modeling using the final inversion result is close to the real record ( Fig.10 , Fig.11 ), Q model extracted at 1.5km position ( Fig.12 ),speed( Fig.13 )The curve fitting effect is better.

[0126] Example 3

[0127] In the specific embodiment 3 of the present invention, the real speed and Q model are as follows: Fig.14 As shown in the figure, the initial model is generally set to a relatively simple gradient model or a relatively smooth velocity model obtained by picking up actual data, respectively using the real ( Fig.14 ) and the initial velocity ( Fig.16 ), the real Q model ( Fig.15 ) and the initial Q model ( Fig.17 ) is forward modeled (Formula 1) to obtain the corresponding observation records and simulation records. The observation system here adopts the well-excitation three-dimensional seismic observation method, that is, the well is excited in the well and the ground receives data. The well is located at x = 3km, the shot spacing is 10m, and there are 40 shots in total. A total of 663 geophones are laid on the surface (full coverage), the track spacing is 10m, and 15Hz Ricker wavelet is used as the source. The final inverted velocity and Q model are as follows Fig.18 , Fig.19 As shown. It can be seen that the inversion results are close to the actual results. The Q model ( Fig. 20 ),speed( Fig.21 )The curve fitting effect is better.

[0128] The viscous acoustic wave velocity and Q full waveform inversion method based on the optimal transport theory in the present invention achieves the purpose of joint inversion of velocity and Q models by reducing the distance between the observed record and the simulated record in the theoretical sense of the optimal transport distance. In the application example, when the initial model and the real model differ greatly, a good application effect can still be achieved, the nonlinearity of the full waveform inversion is weakened, and the stability of the inversion is enhanced. Although the method does not require high accuracy of the initial model, a good initial velocity model and constraint conditions are also conducive to the convergence speed of the global optimization method.

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

[0130] 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 joint inversion method for borehole seismic viscous acoustic wave velocity and Q value, characterized in that: The joint inversion method of seismic viscous acoustic wave velocity and Q value in the well includes: Step 1, forward model the real velocity and Q model to obtain the observation record and construct the initial velocity and Q model; Step 2, using the viscous acoustic wave equation to perform forward modeling to obtain simulated observation records; Step 3, compare the simulated records with the observed records, obtain the data residual based on the optimal transportation distance, and calculate the corresponding accompanying earthquake source; Step 4: back propagate the remaining residual using the conjugate state equation to obtain the corresponding gradient and update the model; Step 5, iterate the above process, exit the loop when the iteration condition is met, and output the corresponding inversion result.

2. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 1, characterized in that: In step 1, the initial velocity and Q models are relatively simple gradient models or relatively smooth velocity models obtained by picking up actual data, and the actual data are obtained by tomography.

3. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value 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 data format is SEGY data body.

4. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value 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 joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 1, characterized in that: In step 2, the initial velocity and Q model are used for forward modeling to obtain the corresponding observation records and simulation records. The observation system adopts the in-well excitation seismic observation method, that is, the seismic observation is excited in the well and the data is received on the ground.

6. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 5, characterized in that: In step 2, for the simulation data, the fractional-order constant-Q viscous acoustic wave equation is used for forward modeling. The specific formula is as follows: Where u is the wave field, c0 is the wave propagation speed in the medium, Q is the quality factor, λ=(ω d / ω0) 2γ ,ω d =2πf d is the angular frequency of the source wavelet, ω0 is the reference angular frequency, the first two terms of the equation mainly control the phase distortion, and the third term controls the amplitude attenuation; the finite difference algorithm is used to differentiate the above equation, where the left side of the equation is the second-order derivative of the wave field with respect to time, which can be approximated by the second-order central difference operator, and the first-order partial derivative of time on the right side is approximated by the backward difference operator, and the fractional-order Laplace operator is calculated by FFT. When the above wave equation is used for numerical simulation, three FFT and three IFFT operations are required in each time step.

7. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 1, characterized in that: The method for joint inversion of well seismic viscous acoustic wave velocity and Q value also includes, after step 2, regularizing the data and selecting an objective function.

8. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 7, characterized in that: When regularizing the data, the oscillating seismic data is converted into a probability distribution to satisfy 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, using a linear weighted method. The specific formula is as follows: g=a*d obs +b f=a*d cal +b. (2) where d obs and d cal are the observed data and simulated data, respectively, 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 transportation theory.

9. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 7, characterized in that: When selecting the objective function, the source δs is calculated by the partial derivative of the objective function to the data. The objective function here is defined according to the optimal transportation distance. The objective function formula of the optimal transportation theory is J OT 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 f: g: M represents the transportation strategy that maps the probability distribution function f to the probability distribution function g.

10. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 9, characterized in that: In step 3, it is determined whether the data residual meets the requirements and a threshold is set. When the residual calculated by formula 3 is less than the threshold, the inversion results of velocity and Q are output; when the residual is greater than the threshold, the associated earthquake source is calculated.

11. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 10, characterized in that: In step 3, when calculating the accompanying earthquake source, the conjugate state method is used to perform back propagation on the accompanying earthquake source to obtain an updated velocity model. According to formula (3), the expression of the accompanying earthquake source can be derived as follows: Where δs represents the associated earthquake source, L is an upper triangular matrix, and diag is a diagonal matrix; F and G are cumulative functions of the probability distribution functions f and g; for each y value, the interpolation search method is used to find the x value corresponding to the monotonically increasing function G(x); for It can be further expressed as:

12. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 11, characterized in that: In step 4, the gradients of the velocity model and the Q model are calculated. Based on the above wave equation, the non-attenuated adjoint equation and gradient expression of the full waveform inversion can be derived using the conjugate state method as follows: Where c0 and Q are the reference speeds; and is a constant fractional Laplace operator; dt is the time difference; u is the calculated wave field; ω d is the angular frequency; λ is the accompanying wave field, and μ can be expressed as: where f d and f o are the frequency and the reference frequency respectively; The adjoint wave field λ satisfies the following wave equation: Where δs is the associated earthquake source, They can be expressed as: Where ω0 is the reference angular frequency; In order to avoid the instability of the growing adjoint operator, a non-decaying adjoint operator is used. By removing the last term in the equation, we can obtain: The companion wave field obtained by this equation can maintain the same phase, but the amplitude is different.

13. The method for joint inversion of borehole seismic viscous acoustic wave velocity and Q value according to claim 12, characterized in that: In step 4, the velocity model and Q model are updated, and the nonlinear conjugate gradient algorithm is selected. The iterative update formula is: k+1 =m k +α k Δm k , where m is the model parameter, α is the iteration step, and Δm is the gradient obtained in step 108; for the inversion strategy problem in multi-parameter inversion, a simultaneous inversion strategy is adopted, that is, the gradients of the velocity and Q models are calculated simultaneously, and the velocity and Q models are updated at the same time.

Images