Regularized seismic full-waveform inversion method and system

Abstract

The invention belongs to the technical field of seismic exploration, and relates to a regularization seismic full waveform inversion method and system, and the method comprises the steps: building a travel time chromatography velocity model as an initial velocity model of an inversion velocity parameter model according to the travel time information of a seismological observation record; extracting sub-band seismic data from the single-shot seismic record; acquiring a seismic forward modeling synthetic data set according to the velocity parameter model and an actual seismic observation system; constructing a regularization objective function, solving a regularization gradient of the regularization objective function and a regularization Hessian matrix of the regularization objective function to obtain an update quantity of the speed parameter model, and updating the speed parameter model for calculation of a next frequency; and a speed parameter model obtained after traversing all integer frequencies in the frequency range of the frequency-band-divided seismic data is an inversion result, so that the inversion precision and the convergence speed are improved.

Description

Technical Field

[0001] The present application belongs to the field of seismic exploration technology, and specifically relates to a regularized seismic full waveform inversion method and system. Background Art

[0002] Full Waveform Inversion (FWI) is an important seismic exploration technology that aims to invert the physical parameters of the underground medium by minimizing the difference between observed seismic data and simulated seismic data. Due to its ability to provide high-resolution underground structure imaging, FWI has important application value in oil and gas exploration, crustal dynamics research, and earthquake disaster prediction. However, FWI is a highly nonlinear and ill-posed inverse problem, and faces the following challenges: FWI is a typical nonlinear optimization problem that is prone to fall into local minima, especially when the initial model accuracy is low. In addition, the low-frequency information in seismic data is usually insufficient, which makes it difficult to recover the long-wavelength structure, exacerbating the difficulty of model convergence.

[0003] FWI is computationally intensive, especially in three-dimensional complex media, where a large number of forward and adjoint problems need to be solved. The high-dimensional computational requirements of the Hessian matrix (especially in three-dimensional scenarios) make traditional Newton-type algorithms difficult to apply in practice, while gradient-type methods (such as the conjugate gradient method) converge slowly and require multiple iterations. Summary of the invention

[0004] The first aspect of the embodiments of the present application provides a regularized seismic full waveform inversion method to solve the problem that the full waveform inversion method is limited by the strong nonlinearity commonly existing in the gradient inversion method, resulting in inaccurate inversion results and difficulty in adapting to the large demand for actual data calculation.

[0005] A first aspect of an embodiment of the present application provides a regularized seismic full waveform inversion system.

[0006] This application is implemented in this way. The first aspect of the embodiment of the present application provides a method for normalized seismic full waveform inversion, comprising: According to the travel time information recorded by seismic observation, a travel time tomographic velocity model is established as the initial velocity model of the inverted velocity parameter model; Extract frequency band seismic data from single shot seismic records; According to the velocity parameter model and the actual seismic observation system, a synthetic data set for earthquake forward modeling is obtained; A regularized objective function is constructed, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum supported velocity parameter model fitting operator; a regularized gradient of the regularized objective function and a regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency; The velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency-band seismic data is the inversion result.

[0007] Furthermore, obtaining a synthetic data set for earthquake forward modeling based on the velocity parameter model and the actual earthquake observation system includes: The calculation medium is divided into background medium and disturbance medium, and the seismic wave field generated by the calculation medium is divided into background wave field and disturbance wave field; Generate source function using wavelet dominant frequency, simulate seismic wave field through source function, and calculate background wave field through background medium Green function; The disturbance wave field is obtained by integrating the interaction between the Green's function of the background medium between two points in the calculation area and the velocity parameter model of the target point within the scope. The background wave field and the disturbance wave field are superimposed to obtain a synthetic data set for seismic forward modeling.

[0008] Furthermore, a regularization parameter is set for the model fitting item, and the regularization parameter is updated in each iteration, and the initial value of the regularization parameter is set to the ratio between the data fitting item and the model fitting item under the initial velocity model.

[0009] Furthermore, the data fitting term uses the second norm of the residual to measure the goodness of fit.

[0010] Furthermore, the regularized objective function is expressed as: , Synthetic dataset for earthquake forward modeling, To actually observe earthquake data, represents the updated velocity parameter model, represents the initial speed parameter model at the current frequency, represents the regularized objective function, It represents the stabilization parameter, which is used to measure the difference between the updated speed parameter model and the previous speed parameter model to adjust the update direction of the model.

[0011] Furthermore, the updating amount of the velocity parameter model obtained by solving the regularized gradient of the regularized objective function and the Hessian matrix of the regularized objective function includes: The gradient of the data fitting term is obtained through the data fitting term, and the Hessian matrix of the data fitting term is obtained through the conjugate transposition of the sensitivity kernel and the sensitivity kernel of any source-receiver pair of the actual observed seismic data set with respect to the velocity parameter model; The gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained through the model fitting term; The regularized gradient of the regularized objective function is obtained by superimposing the gradient of the data fitting term and the gradient of the model fitting term; the regularized Hessian matrix of the regularized objective function is obtained by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term; The update amount of the velocity parameter model is obtained using the regularized gradient and regularized Hessian matrix.

[0012] Furthermore, the gradient of the data fitting term is the accumulation of the sensitivity kernel of any source-detector pair of the actual observed seismic data with respect to the velocity parameter model on the residuals of the actual observed seismic data set and the seismic forward simulation synthetic data set at all source and detector positions; the Hessian matrix of the data fitting term is calculated as the product of the conjugate transpose of the sensitivity kernel and the sensitivity kernel itself.

[0013] Furthermore, the gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained by taking the first-order derivative and the second-order derivative of the model fitting term, respectively.

[0014] In a second aspect of the embodiment of the present application, a regularized seismic full waveform inversion system includes: an initial velocity model construction module for establishing a travel time tomographic velocity model as an initial velocity model of an inverted velocity parameter model according to travel time information recorded by seismic observations; A frequency band seismic data extraction module, used to extract frequency band seismic data from single shot seismic records; A simulation data synthesis module is used to obtain a synthetic data set of earthquake forward simulation based on a velocity parameter model and an actual earthquake observation system; An updating module is provided to construct a regularized objective function, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of a seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum support velocity parameter model fitting operator; the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency; The inversion module uses a velocity parameter model obtained by traversing all integer frequencies within the frequency range of the frequency-divided-band seismic data for full waveform inversion of the seismic data.

[0015] The above one or more technical solutions in the present application have at least the following beneficial effects: the method of the embodiment of the present application improves the inversion accuracy and convergence speed; the constructed regularized objective function better fits the update amount of the velocity parameter model, enhances the recovery ability of high wavenumber components, and improves the inversion effect and stability. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A flowchart of a seismic full waveform inversion method provided in an embodiment of the present application; Figure 2 A flow chart of a method for obtaining a synthetic data set for earthquake forward modeling provided in an embodiment of the present application; Figure 3 A test speed model provided for an embodiment of the present application; Figure 4 An initial velocity model provided in an embodiment of the present application; Figure 5 (a) real part and (b) imaginary part of the frequency domain wave field at the first frequency provided by an embodiment of the present application; Figure 6 A speed parameter model at a first frequency provided in an embodiment of the present application; Figure 7 (a) real part and (b) imaginary part of the frequency domain wave field at the last frequency provided by an embodiment of the present application; Figure 8 The velocity parameter model at the last frequency provided in the embodiment of the present application, i.e., the final inversion result; Fig. 9 The velocity profile comparison between the inversion results of the initial velocity model and the true velocity model provided in the embodiment of the present application; wherein (a) is a longitudinal profile, and (b) is a transverse profile; Fig.10 The convergence of the objective function at different frequencies provided in the embodiments of the present application; Fig.11 A structural block diagram of a seismic full waveform inversion system provided in an embodiment of the present application. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical solution and advantages of the present application more clear, the present application is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described here are only used to explain the present application and are not used to limit the present application.

[0018] See also Figure 1The flowchart of the seismic full waveform inversion method shown in the figure, a regularized seismic full waveform inversion method provided in the embodiment of the present application, based on the preprocessing LBFGS optimization (Limited-memory Broyden–Fletcher–Goldfarb–Shanno, a limited memory quasi-Newton method), includes: S101 establishes a travel time tomographic velocity model based on the travel time information recorded by the seismic observation as the initial velocity model of the inverted velocity parameter model; see Figure 4 is the initial velocity model in an application scenario.

[0019] It should be noted that the executor of the seismic full waveform inversion method provided in the embodiment of the present application may be a server or a computer device, such as a smart phone, a tablet computer, a laptop computer, etc.

[0020] Seismic observation records refer to seismic waveform data recorded by seismic observation systems (including seismographs, geophones, etc.). These seismic waveform data reflect the propagation of seismic waves in underground media, including reflected waves, refracted waves, direct waves, etc. Usually, multiple channels (i.e., data recorded by multiple geophones) are included to form a seismic gather. These seismic gathers can be common center point (CMP) gathers, common shot point gathers, etc.

[0021] The travel-time tomography velocity model obtains the propagation time (travel time) of seismic waves from the excitation point to the receiving point through seismic observation records, and uses the travel-time data in combination with geological models and mathematical algorithms to invert the velocity distribution of the underground medium.

[0022] The time-to-tomography velocity model established by the travel time information recorded by seismic observation is used as the initial velocity model for inversion. The initial velocity model is the starting point of the inversion, and the inversion process is gradually optimized starting from the initial velocity model.

[0023] The speed parameter model is the initial speed model at the initial iteration, and the speed parameter model is updated after each iteration.

[0024] S102 extracts frequency-band seismic data from the single-shot seismic record, wherein the frequency range of the frequency-band seismic data is 4-25 Hz; The seismic observation record contains multiple single-shot seismic records. The frequency band seismic data is extracted from each single-shot seismic record. Since the seismic observation record is a time domain signal, the time domain signal needs to be converted into a frequency domain signal. The single-shot seismic record can be Fourier transformed by Fourier transform to extract the 4-25Hz frequency band seismic data. The 4-25Hz range is defined by empirical values. It can be understood that the required information is included in this range, but it is not limited to the 4-25Hz range.

[0025] S103, obtaining a synthetic data set of earthquake forward simulation according to the velocity parameter model and the actual earthquake observation system; It can be understood that in the initial calculation, the speed parameter model refers to the initial speed model. After the speed parameter model is updated by the previous calculation result, the updated speed parameter model is obtained. Therefore, the speed parameter model here is the latest speed parameter model.

[0026] The actual seismic observation system is a seismic observation system that includes the actual earthquake source location and the detector location. The calculated data obtained by simulating the data of the actual seismic observation system and the velocity parameter model is the seismic forward simulation synthetic data set.

[0027] S104 constructs a regularized objective function, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between the synthetic data set of the seismic forward simulation and the actual observed seismic data set, and the model fitting term adopts a minimum support velocity parameter model fitting operator; solves the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function to obtain an update amount of the velocity parameter model, and updates the velocity parameter model for calculation of the next frequency; The actual observed seismic data set here refers to the observed seismic data set actually collected by the corresponding actual seismic observation system.

[0028] The minimum support model fitting operator is used to evaluate the rate of change of model update, which can avoid the model from updating high-wavenumber components too early, while better fitting the sharp edges of the model, and dynamically adjust the regularization parameters or other control parameters in combination with the adaptive L-curve strategy to improve the stability and robustness of the fitting process.

[0029] S105 traverses all integer frequencies within the frequency range of the sub-band seismic data to obtain a velocity parameter model for full waveform inversion of the seismic data.

[0030] In one embodiment, see Figure 2 The flowchart of the method for obtaining a synthetic data set for seismic forward modeling is shown in S103, wherein the synthetic data set for seismic forward modeling is obtained according to the velocity parameter model and the actual seismic observation system, including: S201 divides the calculation medium into background medium and disturbance medium, and divides the seismic wave field generated by the calculation medium into background wave field and disturbance wave field, that is, the seismic wave field of the calculation area is divided into the background wave field generated by the background medium and the disturbance wave field generated by the disturbance medium. In an application scenario, for example, Figure 3 The test velocity model shown is used as a calculation medium, the calculation medium is divided into a background medium and a disturbance medium, and the seismic wave field generated by the calculation medium is divided into a background wave field and a disturbance wave field; S202 generates a source function using the wavelet main frequency, simulates the seismic wave field through the source function, and calculates the background wave field through the background medium Green's function; S203 uses the Green's function of the background medium between two points in the calculation area and the velocity parameter model of the target point to perform integral calculation within the scope to obtain the disturbance wave field; S204 adds the background wave field and the disturbance wave field to obtain a synthetic data set for seismic forward modeling.

[0031] Specifically, the calculation area refers to the calculation space used for simulation calculation, and the shape of the calculation area is not limited. For example, the calculation area can be set to a rectangle and divided into equally spaced grids. By establishing a coordinate system, the position vector of the detector coordinates and the point coordinate position vector of the source are set according to the actual seismic observation system. The calculation medium is in the calculation area.

[0032] Set the main frequency of the wavelet of the earthquake source. The main frequency of the wavelet refers to the frequency value with the largest amplitude in the amplitude spectrum of the earthquake wavelet. It reflects the frequency range in which the energy of the earthquake wavelet is most concentrated. After setting the main frequency of the wavelet, the earthquake source function is obtained. , used to simulate seismic wave fields , calculated based on the frequency domain acoustic wave Helmholtz equation: , in, represents the spatial derivative operator, is any point in the computational domain, is the wave number, The calculation area is divided into background medium and disturbance medium, and the seismic wave field can be Background wave field and disturbance wave field ,Right now: , As a disturbance wave field, the size is and Green's function of the background medium and Velocity parameter model of a point The integral calculation results of the interaction within the scope D are: The calculation area is divided by For any point outside the , in, express The seismic wave field at the point, and Green's function of the background medium Solved using the zero-order Hankel function of the first kind (representing the propagation characteristics of cylindrical waves): .

[0033] in, , is the zero-order Hankel function of the first kind, represents the background medium wave number, Represents an imaginary number.

[0034] The seismic wave field in the frequency domain under a given velocity parameter model can be the superposition of the background wave field and the disturbance wave field, that is, the seismic wave field data.

[0035] The background wave field is calculated by the background medium Green's function ,in The source function is generated for the wavelet main frequency. The source function is used to describe the process of energy release at the source when an earthquake occurs. By selecting a suitable wavelet main frequency, an earthquake wavelet with a specific frequency characteristic can be generated. In this embodiment, the wavelet main frequency is selected to obtain a source function representing the earthquake wavelet.

[0036] The background wave field and the disturbance wave field are superimposed to obtain the integral equation of the synthetic data set of seismic forward simulation. The integral equation of the synthetic data set of seismic forward simulation belongs to the particle displacement field in the frequency domain and can be expressed as: , through the integral equation of the seismic forward simulation synthetic data set, the seismic forward simulation synthetic data set under different wavelet main frequencies and velocity parameter models can be obtained.

[0037] In one embodiment, the integral equation of the seismic forward simulation synthetic data set can be expressed in the form of an operator, and in order to speed up the operation, the Krylov subspace iterative solver, that is, the generalized minimum residual method, is used to solve the integral equation of the seismic forward simulation synthetic data set expressed in the form of an operator to obtain the seismic forward simulation synthetic data set. In an application scenario, see Figure 5 The frequency domain wave field at the first frequency of Figure 5 (a) in the equation is the real part, Figure 5 (b) is the imaginary part, showing a 4 Hz low-frequency data set; see Figure 7 The frequency domain wave field at the last frequency, Figure 7 (a) in the equation is the real part, Figure 7 (b) in the figure is the imaginary part, which is a 20Hz high frequency data set.

[0038] In one embodiment, S104 constructs a regularized objective function, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of a seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum support velocity parameter model fitting operator; the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency.

[0039] The purpose of calculating the residual between the synthetic data set of earthquake forward simulation and the actual observed earthquake data set is to make the synthetic data set of earthquake forward simulation close to the actual observed earthquake data set. In this process, the velocity parameter model is continuously updated to obtain an updated velocity parameter model.

[0040] The data fitting term is constructed through the residual, the gradient of the data fitting term is obtained through the data fitting term, and the Hessian matrix of the data fitting term is obtained through the sensitivity kernel and the conjugate transpose of the sensitivity kernel of any source-receiver pair of the actual observed seismic data set with respect to the velocity parameter model; The gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained through the model fitting term; The regularized gradient of the regularized objective function is obtained by superimposing the gradient of the data fitting term and the gradient of the model fitting term; the Hessian matrix of the regularized objective function is obtained by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term; The update of the speed parameter model is obtained by solving the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function through the quasi-Newton method. This process is also the process of solving the speed parameter model when the regularized objective function is minimized by using the regularized gradient and regularized Hessian matrix information through the LBFGS method.

[0041] The update amount of the speed parameter model is calculated by an optimization algorithm. In one embodiment, a quasi-Newton method is used to find a local optimal solution using the first-order derivative and the second-order derivative (corresponding to the regularized gradient and the regularized Hessian matrix) of the regularized objective function.

[0042] The regularized objective function of the velocity parameter model established above includes data fitting terms and model fitting terms, which can be expressed as: ,in is the data fitting term, is the model fitting term, is the regularization parameter, is the speed parameter model after the previous update, is the update amount of the speed parameter model. In the iteration, the update amount of the velocity parameter model can be obtained by solving the regularized gradient of the regularized objective function And the regularized Hessian matrix of the regularized objective function get: represents the updated speed parameter model. The update amount of the speed parameter model can be calculated by the following formula: , Indicates The velocity parameter model at the iteration.

[0043] The velocity parameter model is continuously updated and iterated. When the regularized objective function reaches the set threshold, the inversion result at the first frequency is obtained. This result is used as the input of the second frequency, and a new round of iteration is performed at the next frequency until the iteration of the set maximum frequency is completed. The output velocity parameter model is the final calculation result. The inversion result at the first frequency is the velocity parameter model (4Hz in this example). Figure 6 As shown in the figure, it can be observed that the general outline of the salt dome is clearly visible at the first frequency. After obtaining the inversion result of the first frequency, it is input into the next frequency as the initial model, and a new round of iteration is performed until the inversion model is obtained in the maximum frequency iteration, which is the final inversion result. In a specific application scenario, the inversion result, i.e., the velocity parameter model, is as follows: Figure 8 As shown. Figure 8 It can be observed that the deep areas and edges of the salt dome have been restored more accurately.

[0044] In one embodiment, the regularization objective function includes a data fitting term and a model fitting term. The data fitting term uses the second norm (i.e., the Euclidean distance) of the residual between the actual observed seismic data set and the seismic forward simulation synthetic data set to measure the degree of fit. The model fitting term uses the minimum support velocity parameter model fitting operator to measure the difference between the updated velocity parameter model and the previous velocity parameter model to adjust the update direction of the model, and introduces a stabilization parameter To enhance stability, the regularized objective function is expanded as: , Synthetic dataset for earthquake forward modeling, To actually observe earthquake data, and The same represents the updated velocity parameter model, and It also represents the speed parameter model after the previous update.

[0045] In order to solve the updated velocity parameter model, it is necessary to calculate the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function (i.e., the first-order derivative matrix and the second-order derivative matrix). This can be divided into the model fitting part and the data fitting part for separate calculation.

[0046] The canonical gradient of the data fitting term is the accumulation of the sensitivity kernel of the actual observed seismic data with respect to any source-receiver pair of the velocity parameter model on the residuals of the actual observed seismic data set and the synthetic data set of seismic forward simulation at all sources and receiver positions. The Hessian matrix of the data fitting term uses the first-order Gauss-Newton approximation of the Hessian and is calculated as the product of the conjugate transpose of the sensitivity kernel and the sensitivity kernel itself.

[0047] The model fitting items and their Hessian matrices have no direct connection with the source wavefield data and are only related to the rate of change of the velocity parameter model. Therefore, the first-order and second-order derivatives of the model fitting items can be directly obtained by numerical methods.

[0048] Therefore, the solution formulas for calculating the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function are: , ; ; in as well as are the coordinate positions of the geophone and the source point respectively, is the sensitivity kernel of any source-receiver pair of actual observed seismic data with respect to the velocity parameter model, is the conjugate transpose operator.

[0049] Based on the update formula of the speed parameter model, the update amount of the speed parameter model can be solved. is the residual between the actual observed earthquake dataset and the synthetic dataset of earthquake forward simulation at any source-receiver pair position.

[0050] However, the direct inversion of the Hessian matrix will place a huge burden on computer memory, making it impractical to solve large-scale problems. To solve this problem, in actual operations, the calculation is only performed once at the beginning, and then the Hessian matrix is ​​updated iteratively based on the LBFGS framework. At the same time, the Sherman–Morrison formula is used to implement the inversion of the Hessian matrix, directly avoiding the direct inversion and storage of large matrices.

[0051] Furthermore, the sensitivity kernel needs to be solved. The physical function of the sensitivity kernel is to quantify the update amount of the velocity parameter model when the velocity parameter model perturbation is close to infinitesimal. The relationship between the disturbance wave field can be defined as follows: , in, is the disturbance wave field, and The meaning is the same, defining the difference between the seismic wavefield and the background wavefield in any source-receiver pair.

[0052] In one embodiment, the regularization parameter can be updated using an adaptive L-curve-like strategy. After iterations, the update method is ,in: , Initial value of the regularization parameter is set as the ratio between the data fitting term and the model fitting term of the regularized objective function under the initial velocity model, For the The regularization parameter in the iteration, For the The regularization parameter after iterations.

[0053] In order to more intuitively compare the improvement of the initial velocity model by the method proposed in the embodiment of the present application and the accuracy of the final inversion result, the velocity profiles at the middle positions of the velocity parameter model in the horizontal and vertical directions can be extracted for comparison, such as Fig. 9 The velocity profile comparison between the initial velocity model and the inversion results of the true velocity model shown in FIG. Fig. 9 (a) is the longitudinal section. Fig. 9 (b) in the figure is a transverse section. From the longitudinal section and the transverse section, it can be seen that when the initial velocity model has a low accuracy, the method of the present application can still fit the seismic observation data well, whether in shallow or deep layers.

[0054] The convergence of the regularized objective function at different frequencies during the inversion process is shown in Fig.10 As shown, through the curves between the mean square objective function and the number of iterations at frequencies of 4 Hz, 10 Hz and 24 Hz, it can be seen that the method of the embodiment of the present application can achieve full waveform inversion of a complex velocity parameter model within a limited number of times.

[0055] The regularized seismic full waveform inversion system provided in an embodiment of the present application is described below. The seismic full waveform inversion system described below and the regularized seismic full waveform inversion method described above can be referenced to each other.

[0056] See also Fig.11 The structural block diagram of the seismic full waveform inversion system provided in the embodiment of the present application is shown. A regularized seismic full waveform inversion system provided in the embodiment of the present application includes: An initial velocity model building module is used to build a travel time tomographic velocity model as an initial velocity model of the inverted velocity parameter model according to the travel time information recorded by the seismic observation; A frequency band seismic data extraction module, used to extract frequency band seismic data from single shot seismic records; A simulation data synthesis module is used to obtain a synthetic data set of earthquake forward simulation based on a velocity parameter model and an actual earthquake observation system; An updating module is provided to construct a regularized objective function, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of a seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum support velocity parameter model fitting operator; the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency; The inversion module uses a velocity parameter model obtained by traversing all integer frequencies within the frequency range of the frequency-divided-band seismic data for full waveform inversion of the seismic data.

[0057] In one embodiment, the simulation data synthesis module is further used to divide the seismic wave field of the calculation area into a background wave field generated by the background medium and a disturbance wave field generated by the disturbance medium; Generate source function using wavelet dominant frequency, simulate seismic wave field through source function, and calculate disturbance wave field rate through background medium Green function; The disturbance wave field is calculated by using the integral of the interaction between the background medium Green's function between two points in the calculation area and the velocity parameter model of the target point within the scope; The background wave field and the disturbance wave field are superimposed to obtain a synthetic data set for seismic forward modeling.

[0058] In one embodiment, solving the regularized gradient of the regularized objective function and the Hessian matrix of the regularized objective function to obtain the update amount of the velocity parameter model includes: The canonical gradient of the data fitting term is obtained through the data fitting term, and the Hessian matrix of the data fitting term is obtained through the sensitivity kernel of any source-receiver pair of the actual observed seismic data set with respect to the velocity parameter model and the conjugate transpose of the sensitivity kernel; The gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained through the model fitting term; The regularized gradient of the regularized objective function is obtained by superimposing the gradient of the data fitting term and the gradient of the model fitting term; the regularized Hessian matrix of the regularized objective function is obtained by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term; The update amount of the speed parameter model is obtained by using the regularized gradient and regularized Hessian matrix. This process is also the process of solving the speed parameter model when the regularized objective function is minimized by using the LBFGS method and the regularized gradient and regularized Hessian matrix information.

[0059] In one embodiment, the update module is used to calculate the gradient of the data fitting term, which is the accumulation of the sensitivity kernel of any source-detector pair of the actual observed seismic data with respect to the velocity parameter model on the residuals of the actual observed seismic data set and the seismic forward simulation synthetic data set at all source and detector positions; the gradient Hessian matrix of the data fitting term is calculated as the product of the conjugate transpose of the sensitivity kernel and the sensitivity kernel itself.

[0060] The seismic full waveform inversion system is used to better fit the update amount of the velocity parameter model, enhance the recovery ability of high wavenumber components, and improve the inversion effect and stability.

[0061] The above description is only a preferred embodiment of the present application and is not intended to limit the present application. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A regularized seismic full waveform inversion method, characterized in that: The method includes: According to the travel time information recorded by seismic observation, a travel time tomographic velocity model is established as the initial velocity model of the inverted velocity parameter model; Extract frequency band seismic data from single shot seismic records; According to the velocity parameter model and the actual seismic observation system, a synthetic data set for earthquake forward modeling is obtained; A regularized objective function is constructed, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum supported velocity parameter model fitting operator; a regularized gradient of the regularized objective function and a regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency; The velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency-band seismic data is the inversion result.

2. A regularized seismic full waveform inversion method according to claim 1, characterized in that: The method of obtaining a synthetic data set for earthquake forward modeling based on the velocity parameter model and the actual earthquake observation system includes: The calculation medium is divided into background medium and disturbance medium, and the seismic wave field generated by the calculation medium is divided into background wave field and disturbance wave field; Generate source function using wavelet dominant frequency, simulate seismic wave field through source function, and calculate background wave field through background medium Green function; The disturbance wave field is obtained by integrating the interaction between the Green's function of the background medium between two points in the calculation area and the velocity parameter model of the target point within the scope. The background wave field and the disturbance wave field are superimposed to obtain a synthetic data set for seismic forward modeling.

3. A regularized seismic full waveform inversion method according to claim 1, characterized in that: A regularization parameter is set for the model fitting item, and the regularization parameter is updated in each iteration. The initial value of the regularization parameter is set to the ratio between the data fitting item and the model fitting item under the initial velocity model.

4. A regularized seismic full waveform inversion method according to claim 1, characterized in that: The data fitting item uses the second norm of the residual to measure the goodness of fit.

5. A regularized seismic full waveform inversion method according to claim 1, characterized in that: The regularized objective function is expressed as: , Synthetic dataset for earthquake forward modeling, To actually observe earthquake data, represents the updated velocity parameter model, represents the initial speed parameter model at the current frequency, represents the regularized objective function, It represents the stabilization parameter, which is used to measure the difference between the updated speed parameter model and the previous speed parameter model to adjust the update direction of the model.

6. A regularized seismic full waveform inversion method according to claim 1, characterized in that: The updating amount of the velocity parameter model obtained by solving the regularized gradient of the regularized objective function and the Hessian matrix of the regularized objective function includes: The gradient of the data fitting term is obtained through the data fitting term, and the Hessian matrix of the data fitting term is obtained through the conjugate transposition of the sensitivity kernel and the sensitivity kernel of any source-receiver pair of the actual observed seismic data set with respect to the velocity parameter model; The gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained through the model fitting term; The regularized gradient of the regularized objective function is obtained by superimposing the gradient of the data fitting term and the gradient of the model fitting term; the regularized Hessian matrix of the regularized objective function is obtained by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term; The update amount of the velocity parameter model is obtained using the regularized gradient and regularized Hessian matrix.

7. A regularized seismic full waveform inversion method according to claim 6, characterized in that: The gradient of the data fitting term is the accumulation of the sensitivity kernel of any source-detector pair of the actual observed seismic data with respect to the velocity parameter model on the residuals of the actual observed seismic data set and the seismic forward simulation synthetic data set at all source and detector positions; the Hessian matrix of the data fitting term is calculated as the product of the conjugate transpose of the sensitivity kernel and the sensitivity kernel itself.

8. A regularized seismic full waveform inversion method according to claim 6, characterized in that: The gradient of the model fitting term and the Hessian matrix of the model fitting term are obtained by taking the first-order derivative and the second-order derivative of the model fitting term, respectively.

9. A regularized seismic full waveform inversion system, characterized in that: include: An initial velocity model building module is used to build a travel time tomographic velocity model as an initial velocity model of the inverted velocity parameter model according to the travel time information recorded by the seismic observation; A frequency band seismic data extraction module, used to extract frequency band seismic data from single shot seismic records; A simulation data synthesis module is used to obtain a synthetic data set of earthquake forward simulation based on a velocity parameter model and an actual earthquake observation system; An updating module is provided to construct a regularized objective function, wherein the regularized objective function includes a data fitting term and a model fitting term, wherein the data fitting term is constructed by the residual between a synthetic data set of a seismic forward simulation and an actual observed seismic data set, and the model fitting term adopts a minimum support velocity parameter model fitting operator; the regularized gradient of the regularized objective function and the regularized Hessian matrix of the regularized objective function are solved to obtain an update amount of the velocity parameter model, and the velocity parameter model is updated for calculation of the next frequency; The inversion module uses a velocity parameter model obtained by traversing all integer frequencies within the frequency range of the frequency-divided-band seismic data for full waveform inversion of the seismic data.

Images