Regularized seismic full waveform inversion method and system

Abstract

This application belongs to the field of seismic exploration technology, and provides a regularized seismic full waveform inversion method and system, including: establishing a travel time tomography velocity model as the initial velocity model of the velocity parameter model for inversion according to the travel time information of seismic observation records; extracting frequency band-divided seismic data from single-shot seismic records; obtaining a synthetic seismic forward modeling dataset according to the velocity parameter model and the actual seismic observation system; constructing a regularized objective function, solving the regular gradient of the regularized objective function and the regular Hessian matrix of the regularized objective function to obtain the update amount of the velocity parameter model, and updating the velocity parameter model for the calculation of the next frequency; the velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency band-divided seismic data is the inversion result, which improves the inversion accuracy and convergence speed.

Description

Technical Field

[0001] This application belongs to the technical field of seismic exploration, and more specifically, it is a regularization seismic full waveform inversion method and system. Background Art

[0002] Seismic full waveform inversion (FWI) is an important seismic exploration technique aimed at inversing the physical parameters of the subsurface medium by minimizing the difference between observed seismic data and simulated seismic data. Due to its ability to provide high-resolution subsurface structure imaging, FWI has important application values in fields such as oil and gas exploration, crustal dynamics research, and earthquake disaster prediction. However, FWI is a highly nonlinear and ill-posed inverse problem, mainly facing the following challenges:

[0003] FWI is a typical nonlinear optimization problem and is prone to falling into local minima, especially when the accuracy of the initial model is low. In addition, the low-frequency information in seismic data is usually insufficient, resulting in difficulties in recovering long-wavelength structures and exacerbating the model convergence difficulty.

[0004] The computational cost of FWI is huge, 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-like algorithms difficult to apply in practice, while gradient-based methods (such as the conjugate gradient method) have a slow convergence speed and require multiple iterations. Summary of the Invention

[0005] The first aspect of the embodiments of this application provides a regularization seismic full waveform inversion method, which solves the problems that the full waveform inversion method is limited by the strong nonlinearity commonly existing in gradient-based inversion methods, resulting in inaccurate inversion results and difficulty in adapting to the large computational requirements of actual data.

[0006] The first aspect of the embodiments of this application provides a regularization seismic full waveform inversion system.

[0007] This application is implemented as follows:

[0008] The first aspect of the embodiments of this application provides a regularization seismic full waveform inversion method, including:

[0009] Based on the travel time information of seismic observation records, establish a travel time tomography velocity model as the initial velocity model of the velocity parameter model for inversion;

[0010] Extract frequency-band seismic data from single-shot seismic records;

[0011] Based on the velocity parameter model and the actual seismic observation system, obtain a seismic forward simulation synthetic data set;

[0012] Construct a regularization objective function, which includes a data fitting term and a model fitting term. The data fitting term is constructed based on the residual between the synthetic dataset obtained from seismic forward modeling and the actual observed seismic dataset. The model fitting term uses a minimum support velocity parameter model fitting operator; solve the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function to obtain the update amount of the velocity parameter model, and update the velocity parameter model for the calculation of the next frequency;

[0013] The velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency-divided seismic data is the inversion result.

[0014] Further, obtaining the synthetic dataset of seismic forward modeling according to the velocity parameter model and the actual seismic observation system includes:

[0015] Divide the calculation medium into a background medium and a perturbed medium, and divide the seismic wave field generated by the calculation medium into a background wave field and a perturbed wave field;

[0016] Generate a source function using the main frequency of the wavelet, simulate the seismic wave field through the source function, and calculate the background wave field through the background medium Green's function;

[0017] Use the integral calculation 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 of action to obtain the perturbed wave field;

[0018] Superimpose the background wave field and the perturbed wave field to obtain the synthetic dataset of seismic forward modeling.

[0019] Further, set a regularization parameter for the model fitting term, and update the regularization parameter in each iteration. The initial value of the regularization parameter is set as the ratio between the data fitting term and the model fitting term under the initial velocity model.

[0020] Further, the data fitting term uses the second norm of the residual to measure the fitting degree.

[0021] Further, the regularization objective function is expressed as: ,

[0022] is the synthetic dataset of seismic forward modeling, is the actual observed seismic data, represents the updated velocity parameter model, represents the initial velocity parameter model at the current frequency, represents the regularization objective function, It represents a stabilization parameter, which is used to measure the difference between the updated velocity parameter model and the previous velocity parameter model and adjust the update direction of the model.

[0023] Furthermore, obtaining the update amount of the velocity parameter model by solving the regular gradient of the regularized objective function and the Hessian matrix of the regularized objective function includes:

[0024] Obtaining the gradient of the data fitting term through the data fitting term, and obtaining the Hessian matrix of the data fitting term through the conjugate transpose of the sensitivity kernel of any source-detector pair of the actual observed seismic data set with respect to the velocity parameter model and the sensitivity kernel;

[0025] Obtaining the gradient and Hessian matrix of the model fitting term through the model fitting term;

[0026] Obtaining the regular gradient of the regularized objective function by superimposing the gradient of the data fitting term and the gradient of the model fitting term; obtaining the regular Hessian matrix of the regularized objective function by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term;

[0027] Obtaining the update amount of the velocity parameter model by using the regular gradient and the regular Hessian matrix.

[0028] 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 acting on the residuals between the actual observed seismic data set and the synthetic seismic forward modeling 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.

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

[0030] In the second aspect of the embodiments of the present application, a regularized seismic full waveform inversion system includes: an initial velocity model construction module, which is used to establish a travel time tomography velocity model as the initial velocity model of the velocity parameter model for inversion according to the travel time information of seismic observation records;

[0031] A frequency band seismic data extraction module, which is used to extract frequency band seismic data from single-shot seismic records;

[0032] A simulation data synthesis module, which is used to obtain a synthetic seismic forward modeling data set according to the velocity parameter model and the actual seismic observation system;

[0033] An update module constructs a regularization objective function, which includes a data fitting term and a model fitting term. The data fitting term is constructed based on the residual between the synthetic dataset obtained from seismic forward modeling and the actual observed seismic dataset. The model fitting term uses a minimum support velocity parameter model fitting operator. The update amount of the velocity parameter model is obtained by solving the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function, and the velocity parameter model is updated for the calculation of the next frequency.

[0034] An inversion module uses the velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency-divided seismic data for full waveform inversion of the seismic data.

[0035] One or more of the above technical solutions in this application have at least the following beneficial effects: The method of the embodiment of this application improves the inversion accuracy and convergence speed; the constructed regularization 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. Description of the Drawings

[0036] Figure 1 It is a flowchart of the seismic full waveform inversion method provided by the embodiment of this application;

[0037] Figure 2 It is a flowchart of the method for obtaining the synthetic dataset of seismic forward modeling provided by the embodiment of this application;

[0038] Figure 3 It is the test velocity model provided by the embodiment of this application;

[0039] Figure 4 It is the initial velocity model provided by the embodiment of this application;

[0040] Figure 5 It is the (a) real part and (b) imaginary part of the frequency domain wave field at the first frequency provided by the embodiment of this application;

[0041] Figure 6 It is the velocity parameter model at the first frequency provided by the embodiment of this application;

[0042] Figure 7 It is the (a) real part and (b) imaginary part of the frequency domain wave field at the last frequency provided by the embodiment of this application;

[0043] Figure 8 It is the velocity parameter model at the last frequency provided by the embodiment of this application, that is, the final inversion result;

[0044] Figure 9Velocity profile comparison between the initial velocity model provided by the embodiments of the present application and the inversion result of the true velocity model; wherein, (a) is the longitudinal profile and (b) is the transverse profile.

[0045] Figure 10 Convergence of the objective function at different frequencies provided by the embodiments of the present application.

[0046] Figure 11 Block diagram of the seismic full waveform inversion system provided by the embodiments of the present application. Detailed implementation manners

[0047] In order to make the objectives, technical solutions and advantages of the present application more clear and understandable, the present application will be further described in detail below with reference to the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0048] Refer to Figure 1 As shown in the flowchart of the seismic full waveform inversion method, a regularization seismic full waveform inversion method provided by the embodiments of the present application is based on preprocessed LBFGS optimization (Limited-memory Broyden–Fletcher–Goldfarb–Shanno, a limited-memory quasi-Newton method), and includes:

[0049] S101 Establish a travel time tomography velocity model as the initial velocity model of the velocity parameter model for inversion according to the travel time information of the seismic observation records; refer to Figure 4 As the initial velocity model in an application scenario.

[0050] It should be noted that the execution subject of the seismic full waveform inversion method provided by the embodiments of the present application can be a server or a computer device, such as a smart phone, a tablet computer, a notebook computer, etc.

[0051] Seismic observation records refer to seismic waveform data recorded by a seismic observation system (including seismographs, geophones, etc.). These seismic waveform data reflect the propagation of seismic waves in the underground medium, including reflected waves, refracted waves, direct waves, etc. Usually, it includes multiple channels (i.e., data recorded by multiple geophones), forming a seismic trace gather. These seismic trace gathers can be common midpoint (CMP) gathers, common shot gathers, etc.

[0052] The travel time tomography velocity model is to obtain the propagation time (travel time) of seismic waves from the excitation point to the receiving point through seismic observation records, and use the travel time data, combined with the geological model and mathematical algorithms, to invert the velocity distribution of the underground medium.

[0053] The travel-time tomography velocity model established from the travel-time information recorded in seismic observations is used as the initial velocity model for inversion. The initial velocity model is the starting point of inversion, and the inversion process starts from the initial velocity model and is gradually optimized.

[0054] In the initial iteration, the velocity parameter model is the initial velocity model, and the velocity parameter model is updated after each iteration.

[0055] S102 extracts frequency-band seismic data from single-shot seismic records, and the frequency range of the frequency-band seismic data is 4 - 25 Hz;

[0056] The seismic observation record contains multiple single-shot seismic records. Frequency-band seismic data is extracted from each single-shot seismic record. Since the seismic observation record is a time-domain signal, it is necessary to convert the time-domain signal into a frequency-domain signal. The Fourier transform can be used to perform the Fourier transform on the single-shot seismic record, thereby intercepting the 4 - 25 Hz frequency-band seismic data. The range of 4 - 25 Hz 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 range of 4 - 25 Hz.

[0057] S103 obtains a synthetic seismic forward modeling dataset according to the velocity parameter model and the actual seismic observation system;

[0058] It can be understood that in the initial calculation, the velocity parameter model refers to the initial velocity model. After updating the velocity parameter model based on the calculation result of the previous time, the updated velocity parameter model is obtained. Therefore, the velocity parameter model here is the latest velocity parameter model.

[0059] The actual seismic observation system is a seismic observation system that includes the actual source location and geophone location. The calculated data obtained by simulating data between the actual seismic observation system and the velocity parameter model is the synthetic seismic forward modeling dataset.

[0060] S104 constructs a regularization objective function. The regularization objective function includes a data fitting term and a model fitting term. The data fitting term is constructed from the residuals between the synthetic seismic forward modeling dataset and the actual observed seismic dataset, and the model fitting term uses the minimum support velocity parameter model fitting operator; the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function are solved to obtain the update amount of the velocity parameter model, and the velocity parameter model is updated for the calculation of the next frequency;

[0061] The actual observed seismic dataset here is the observed seismic dataset actually collected by the corresponding actual seismic observation system.

[0062] The minimum support model fitting operator is used to evaluate the change rate of model updates, which can avoid premature model updates of high wavenumber components, better fit the sharp edges of the model, and dynamically adjust regularization parameters or other control parameters in combination with an adaptive L-curve-like strategy to improve the stability and robustness of the fitting process.

[0063] S105 obtains a velocity parameter model by traversing all integers within the frequency range of the frequency-divided seismic data for full waveform inversion of seismic data.

[0064] In one embodiment, referring to Figure 2 the method flowchart for obtaining a synthetic seismic forward modeling dataset shown in

[0065] S201 divides the computational medium into a background medium and a perturbed medium, and divides the seismic wavefield generated by the computational medium into a background wavefield and a perturbed wavefield, that is, divides the seismic wavefield in the computational region into the background wavefield generated by the background medium and the perturbed wavefield generated by the perturbed medium. In an application scenario, for example, as Figure 3 shown in the test velocity model, the computational medium is divided into a background medium and a perturbed medium, and the seismic wavefield generated by the computational medium is divided into a background wavefield and a perturbed wavefield;

[0066] S202 generates a source function using the dominant frequency of the wavelet, simulates the seismic wavefield through the source function, and calculates the background wavefield through the background medium Green's function;

[0067] S203 performs integral calculation using the interaction between the background medium Green's function between two points in the computational region and the velocity parameter model of the target point within the scope of action to obtain the perturbed wavefield;

[0068] S204 adds the background wavefield and the perturbed wavefield together to obtain the synthetic seismic forward modeling dataset.

[0069] Specifically, the computational region refers to the computational space used for simulation calculations, and the shape of the computational region is not limited. For example, the computational region can be set as a rectangle and equally spaced grid division can be performed. By establishing a coordinate system, the position vector of the geophone coordinates and the point coordinate position vector of the source are set according to the actual seismic observation system. The computational medium is located within the computational region.

[0070] Set the dominant frequency of the source wavelet. The dominant frequency of the wavelet refers to the frequency value with the largest amplitude in the amplitude spectrum of the seismic wavelet. It reflects the frequency range where the energy of the seismic wavelet is most concentrated. After setting the dominant frequency of the wavelet, the source function is obtained, which is used to simulate the seismic wavefield and is calculated based on the frequency-domain acoustic wave Helmholtz equation:

[0071] ,

[0072] wherein, represents the spatial derivative operator, is an arbitrary point in the calculation region, is the wave number, is the dominant frequency of the wavelet. The calculation region is divided into a background medium and a perturbed medium, and the seismic wave field can be decomposed into a background wave field and a perturbed wave field , that is: , As the perturbed wave field, the magnitude is and the background medium Green's function between and the velocity parameter model at point within the scope of action D of the integral calculation result of the interaction, is an arbitrary point in the calculation region except , and the calculation formula is;

[0073] ,

[0074] wherein, represents the seismic wave field at point and the background medium Green's function between is solved by the zero-order Hankel function of the first kind (representing the propagation characteristics of cylindrical waves):

[0075] .

[0076] wherein, , is the zero-order Hankel function of the first kind, represents the background medium wave number, represents the imaginary number.

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

[0078] The background wave field is calculated through the background medium Green's function , wherein is the source function generated by the dominant frequency of the wavelet. The source function is used to describe the process of energy release at the earthquake source when an earthquake occurs. By selecting an appropriate dominant frequency of the wavelet, a seismic wavelet with specific frequency characteristics can be generated. In this embodiment, the dominant frequency of the wavelet is selected to obtain the source function representing the seismic wavelet.

[0079] The integral equation for obtaining the synthetic seismic forward modeling dataset by superimposing the background wavefield and the perturbed wavefield. The integral equation of the synthetic seismic forward modeling dataset belongs to the frequency-domain particle displacement field and can be expressed as:

[0080] , and the synthetic seismic forward modeling dataset under different dominant frequencies of the wavelet and velocity parameter models can be obtained through the integral equation of the synthetic seismic forward modeling dataset.

[0081] In one embodiment, the integral equation of the synthetic seismic forward modeling dataset can be represented in operator form. To accelerate the operation speed, the Krylov subspace iterative solver, i.e., the generalized minimum residual method, is used to solve the integral equation of the synthetic seismic forward modeling dataset represented in operator form to obtain the synthetic seismic forward modeling dataset. In an application scenario, refer to Figure 5 the frequency-domain wavefield at the first frequency of Figure 5 where (a) in Figure 5 is the real part, Figure 7 and (b) in Figure 7 is the imaginary part, showing a low-frequency dataset of 4 Hz; refer to Figure 7 the frequency-domain wavefield at the last frequency of

[0082] where (a) in

[0083] is the real part,

[0084] and (b) in

[0085] is the imaginary part, which is a high-frequency dataset of 20 Hz.

[0086] The regular gradient of the regularization 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 regularization objective function is obtained by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term.

[0087] The update of the velocity parameter model is obtained by solving the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function through the quasi - Newton method. This process is also a process of using the LBFGS method to solve for the velocity parameter model that minimizes the regularization objective function by using the information of the regular gradient and the regular Hessian matrix.

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

[0089] The regularization objective function of the velocity parameter model established above includes a data fitting term and a model fitting term, which is expressed as: , where is the data fitting term, is the model fitting term, is the regularization parameter, is the velocity parameter model after the previous update, is the update of the velocity parameter model. At the -th iteration, the update of the velocity parameter model can be obtained by solving the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function: represents the updated velocity parameter model. The update of the velocity parameter model can be calculated by the following formula:

[0090] ,

[0091] represents the velocity parameter model at the -th iteration.

[0092] The velocity parameter model is continuously updated and iterated. When the regularization objective function reaches the set threshold, the inversion result at the first frequency is obtained. This result is used as the input for the second frequency, and a new round of iteration for the next frequency is carried out 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, that is, the velocity parameter model (4 Hz in the example) is as Figure 6As shown, it can be observed that at the first frequency, the general outline of the salt dome is already clearly visible. After obtaining the inversion result of the first frequency, it is input as the initial model for the next frequency for a new round of iteration 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 is the velocity parameter model as Figure 8 shown. According to Figure 8 it can be observed that the deep region of the salt dome and the edge of the salt dome are both accurately restored.

[0093] In one embodiment, the regularization objective function includes a data fitting term and a model fitting term. Among them, the data fitting term uses the two-norm (i.e., Euclidean distance) of the residual between the actual observed seismic data set and the synthetic seismic forward modeling data set to measure the fitting degree. 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 a stabilization parameter is introduced to enhance stability. The regularization objective function is expanded as: ,

[0094] is the synthetic seismic forward modeling data set, is the actual observed seismic data, and are the same and represent the updated velocity parameter model, and are the same and represent the previously updated velocity parameter model.

[0095] To solve for the updated velocity parameter model, the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function (i.e., the first-order derivative matrix and the second-order derivative matrix) need to be calculated. Here, it can be calculated separately for the model fitting term part and the data fitting term part.

[0096] The regular gradient of the data fitting term is the sum of the sensitivity kernels of the actual observed seismic data with respect to any source-receiver pair of the velocity parameter model acting on the residual between the actual observed seismic data set and the synthetic seismic forward modeling data set at all source 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.

[0097] Since the model fitting term and the Hessian matrix of the model fitting term have no direct connection with the source wavefield data and are only related to the rate of change of the velocity parameter model, the first-order derivative and the second-order derivative of the model fitting term can be directly obtained by numerical methods.

[0098] Therefore, the solution formulas for the regular gradient of the regularized objective function and the regular Hessian matrix of the regularized objective function are as follows:

[0099] ,

[0100] ;

[0101] ;

[0102] where and are the coordinate positions of the geophone and the source point respectively, is the sensitivity kernel of any source-receiver pair of the actual observed seismic data with respect to the velocity parameter model, is the conjugate transpose operator.

[0103] Based on the update formula of the velocity parameter model, the update amount of the velocity parameter model can be solved, is the residual between the actual observed seismic data set and the synthetic data set of seismic forward modeling at the position of any source-receiver pair.

[0104] However, directly inverting the Hessian matrix places a huge burden on the computer memory, making it impractical to solve large-scale problems. To solve this problem, in actual operations, it is only calculated once at the beginning, and then based on the LBFGS framework, an iterative method is used to update the Hessian matrix. At the same time, the Sherman–Morrison formula is used to implement the inversion operation of the Hessian matrix, directly avoiding the direct inversion operation and storage of large matrices.

[0105] Furthermore, the sensitivity kernel needs to be solved. The physical role of the sensitivity kernel is to quantify the relationship between the update amount of the velocity parameter model and the perturbed wave field in the case where the perturbation of the velocity parameter model approaches infinitesimal, and can be defined by the mathematical expression:

[0106] ,

[0107] where is the perturbed wave field, which has the same meaning as and defines the difference between the seismic wave field and the background wave field in any source-receiver pair.

[0108] In one embodiment, the regularization parameter can be updated using an adaptive L-curve-like strategy. The update method after the th iteration is , where:

[0109] ,

[0110] Initial value of the regularization parameter is set to the ratio between the data fitting term and the model fitting term of the regularization objective function under the initial velocity model, is the regularization parameter in the -th iteration, is the regularization parameter after the -th iteration.

[0111] To more intuitively compare the improvement of the method proposed in the embodiments of the present application on the initial velocity model and the accuracy of the final inversion result, velocity profiles located at the middle positions in the horizontal and vertical directions of the velocity parameter model can be extracted for comparison, as shown in Figure 9 the comparison of velocity profiles between the inversion results of the initial velocity model and the true velocity model shown; among them, Figure 9 (a) in is the vertical profile, Figure 9 (b) in is the horizontal profile. It can be seen from the vertical profile and the horizontal profile that in the case of low accuracy of the initial velocity model, whether in the shallow layer or the deep layer, the method of the present application can still well fit the seismic observation data.

[0112] The convergence of the regularization objective function at different frequencies during the inversion process is as shown in Figure 10 . It can be seen from the curves between the mean square objective function and the number of iterations at frequencies 4 Hz, 10 Hz, and 24 Hz that the method of the embodiments of the present application can achieve full waveform inversion of a complex velocity parameter model within a limited number of times.

[0113] The regularization seismic full waveform inversion system provided by the embodiments of the present application will be described below. The seismic full waveform inversion system described below can be correspondingly referred to the regularization seismic full waveform inversion method described above.

[0114] Referring to Figure 11 the structural block diagram of the seismic full waveform inversion system provided by the embodiments of the present application shown, a regularization seismic full waveform inversion system provided by the embodiments of the present application includes:

[0115] An initial velocity model construction module, configured to establish a travel time tomography velocity model as the initial velocity model of the velocity parameter model for inversion according to the travel time information of seismic observation records;

[0116] A frequency band seismic data extraction module, configured to extract frequency band seismic data from single-shot seismic records;

[0117] A simulated data synthesis module, configured to obtain a seismic forward simulation synthesis data set according to the velocity parameter model and the actual seismic observation system;

[0118] An update module constructs a regularization objective function, which includes a data fitting term and a model fitting term. The data fitting term is constructed based on the residual between the synthetic data set obtained by seismic forward modeling and the actual observed seismic data set. The model fitting term uses a minimum support velocity parameter model fitting operator. Solve the regular gradient of the regularization objective function and the regular Hessian matrix of the regularization objective function to obtain the update amount of the velocity parameter model, and update the velocity parameter model for the calculation of the next frequency.

[0119] An inversion module uses the velocity parameter model obtained after traversing all integer frequencies within the frequency range of the frequency-divided seismic data for full waveform inversion of seismic data.

[0120] In one embodiment, the simulation data synthesis module is also used to divide the seismic wave field in the calculation area into a background wave field generated by the background medium and a perturbation wave field generated by the perturbed medium.

[0121] Generate a source function using the main frequency of the wavelet, simulate the seismic wave field through the source function, and calculate the perturbation wave field rate through the Green's function of the background medium.

[0122] Calculate the perturbation wave field using the integral of 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 of action.

[0123] Superimpose the background wave field and the perturbation wave field to obtain the synthetic data set of seismic forward modeling.

[0124] In one embodiment, the obtaining of the update amount of the velocity parameter model by solving the regular gradient of the regularization objective function and the Hessian matrix of the regularization objective function includes:

[0125] Obtain the regular gradient of the data fitting term through the data fitting term, and obtain the Hessian matrix of the data fitting term through 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.

[0126] Obtain the gradient of the model fitting term and the Hessian matrix of the model fitting term through the model fitting term.

[0127] Obtain the regular gradient of the regularization objective function by superimposing the gradient of the data fitting term and the gradient of the model fitting term; obtain the regular Hessian matrix of the regularization objective function by superimposing the Hessian matrix of the data fitting term and the Hessian matrix of the model fitting term.

[0128] The update amount of the velocity parameter model is obtained by using the regularized gradient and the regularized Hessian matrix. This process is also a process of solving the velocity parameter model that minimizes the regularized objective function by using the LBFGS method and the information of the regularized gradient and the regularized Hessian matrix.

[0129] In one embodiment, the gradient of the data fitting term for the update module is the accumulation of the sensitivity kernels of the actual observed seismic data with respect to any source-receiver pair of the velocity parameter model acting on the residuals of the actual observed seismic data set and the seismic forward simulation synthetic data set at all source and receiver 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.

[0130] The update amount of the velocity parameter model is better fitted by the seismic full waveform inversion system, enhancing the recovery ability of the high wavenumber components and improving the inversion effect and stability.

[0131] The above are only the preferred embodiments of the present application and are not intended to limit the present application. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present application shall be included within 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 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, 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; 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 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.

6. A regularized seismic full waveform inversion method according to claim 5, 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.

7. A regularized seismic full waveform inversion method according to claim 5, 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.

8. 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 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, 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; 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