Elastic wave full waveform inversion method based on optimal transmission function

Abstract

The invention relates to an elastic wave full waveform inversion method based on an optimal transmission function, and belongs to the field of seismic velocity modeling. The method comprises the following steps: firstly, establishing an observation system, selecting an excitation source with a certain frequency and a speed model, and receiving an actual seismic record and a simulated seismic record through the excitation source; secondly, introducing an optimal transmission function as a target function of inversion, and establishing a target functional formula; a timing wave field and a reverse time wave field are obtained through forward modeling and reverse time continuation, derivation is carried out on a target functional formula to obtain a gradient calculation formula, the gradient calculation formula is substituted into the timing wave field and the reverse time wave field, and the gradient is obtained. And optimizing the obtained gradient by using an Adam gradient optimization algorithm, and substituting the optimized gradient into iterative inversion to update the model. And repeating the steps until an inversion result meeting the precision requirement is obtained. Model experiments show that compared with a traditional L2 objective function, the method has the advantages that the event continuity of the inversion model is enhanced, the deep structure is clearer, and the overall precision is greatly improved.

Description

Technical Field

[0001] The invention belongs to the field of seismic velocity modeling, and in particular relates to an elastic wave full waveform inversion method based on an optimal transmission function. Background Art

[0002] Elastic Full Waveform Inversion (EFWI) is a high-resolution geophysical imaging technology that can form a detailed geological model underground by utilizing the propagation characteristics of seismic waves excited by the earthquake source. EFWI compares the observed seismic data with the data generated by the forward simulation, and continuously adjusts the model parameters through an optimization iterative process, so that the simulated data gradually approaches the observed data, thereby inverting the physical parameters of the underground structure. This method can accurately characterize the underground structure and provide important information support for oil and gas exploration and geological disaster prediction.

[0003] In the traditional elastic wave full waveform inversion process, the L2 norm is usually used as the objective function to measure the difference between the observed data and the simulated data. However, due to the influence of nonlinearity, multipath and complex geological bodies in the process of seismic wave propagation, the traditional method is prone to fall into local minima during the inversion process, and it is difficult to obtain the global optimal solution. The model update process is easily affected by the local minima, resulting in error accumulation, and the final generated geological model deviates greatly from the actual situation. Summary of the invention

[0004] The present invention provides an elastic wave full waveform inversion method based on the optimal transmission function to solve the above problems. The method first establishes an observation system and obtains seismic records, selects the optimal transmission function as the inversion target functional, derives the target functional, obtains the gradient of the optimal transmission, and finally adopts the gradient optimization algorithm to iteratively update the model, and finally obtains a high-precision elastic wave full waveform inversion result. Compared with the traditional L2 objective function, the continuity of the event axis of the inversion model is enhanced, the deep structure is clearer, and the overall accuracy is greatly improved.

[0005] The present invention is achieved through the following technical solutions: The invention discloses an elastic wave full waveform inversion method based on an optimal transmission function, which specifically comprises the following steps: Step 1. Seismic record acquisition: firstly, an observation system is established and an exciting source of a certain frequency and a velocity model are selected, and actual seismic records and simulated seismic records are received at a detection point by exciting the source; Step 2. Selection of objective function: an optimal transmission function in the mathematical field is introduced as the objective function of inversion, and a target functional formula is established based on actual seismic records and simulated seismic records by considering the overall shape and manifold structure of data distribution; Step 3. Wave field acquisition: forward wave field and reverse wave field are obtained by forward simulation and reverse time extension; Step 4. Gradient acquisition: the target functional formula is derived to obtain a gradient calculation formula, and the gradient is obtained by bringing the forward wave field and reverse wave field into the formula; Step 5. Model update: the required gradient is optimized by using the Adam gradient optimization algorithm and then brought into the iterative inversion to update the model; Steps 3-5 are repeated until an inversion result that meets the accuracy requirements is obtained.

[0006] The present invention also provides a specific elastic wave full waveform inversion method based on an optimal transfer function, the method specifically comprising the following steps: Step 1. Seismic record acquisition: A linear model is selected as the initial model, and then a velocity model is given as the true model. The model contains model parameters Vp(x, z), Vs(x, z), ρ(x, z), x and z represent the coordinates in space, x=[1, 2, 3, …, Nx], z=[1, 2, 3, …, Nz], Nx and Nz represent the number of vertical and horizontal grid points of the velocity model respectively. The shot points are evenly placed on the top of the model, with a total of Nx shots. The detection points are also evenly placed on the top of the model, with a total of Nx detection points. The shot points are excited from left to right at the top of the model. Each time they are excited, all the detection points receive at the same time, and the actual seismic records with a total of Nx shots are obtained. Simulated earthquake records of Nx guns ; Step 2. Selection of objective function: Introduce the optimal transfer function in the mathematical field as the objective function of inversion, based on actual seismic records Simulated earthquake records Establish the target functional (1), where To simulate earthquake records, For actual earthquake records, , for , The probability density function obtained after the preprocessing operator is: is the objective function of elastic wave inversion; simulated data is measured by Wasserstein distance and The degree of matching; is the probability density function and of Suboptimal transmission distance; Indicates the number of receiving points; Indicates the receiving point index; is the wave field data of longitudinal wave velocity, Wave field data of shear wave velocity, Wavefield data of density distribution; is the reverse time wave field distribution of the longitudinal wave velocity, is the reverse time wave field distribution of shear wave velocity, is the reverse time wave field distribution of density distribution; the formula (1): ; Step 3. Wave field calculation: Based on the seismic wavelet and initial velocity model, the elastic wave equation (2) is used for forward modeling, where is the displacement, is the stress tensor, is the source term, is the medium density, For time, represents the differential operator, Represents the divergence of the stress tensor; the positive wave field extension starts from the earthquake source and calculates different parameters The forward wave field under the condition of , , Wavefield data; using actual earthquake records and simulated earthquake records The difference between the two is taken as the earthquake source, and the elastic wave equation (2) is used to perform reverse time wave field extension, and the following are obtained: Reverse time wave field distribution of parameters , and , the formula (2) is: ; Step 4. Gradient calculation: Derivative the optimal transfer function objective functional formula (1) to obtain the gradient calculation formula as shown in formula (3). Substitute the forward and reverse wave fields obtained in step (3) to calculate the model gradient, where the parameters are are model parameters, including The value of; the formula (3) is: ; Step 5. Update model parameters: Update model parameters using gradients according to the gradient optimization algorithm The optimization algorithm is the Adam gradient optimization algorithm, which performs gradient optimization through formula (4) to update the model, where k is the number of iterations. is the initial learning rate, and are the first-order and second-order momentum error correction terms, respectively, is the momentum-corrected gradient; the gradient obtained by the gradient formula in step (4) is , and then calculate each The first moment estimate of In the formula, is a decay factor less than 1, is the current iteration number; use formula (6) to obtain an unbiased estimate through bias correction ; is the value after the deviation correction of the second-order moment estimate in the Adam optimization algorithm, which represents the variance estimate of the current gradient; by , use formula (7) to calculate the first-order displacement average of its second-order moment, is the attenuation factor of the second-order moment, is the current iteration number; since There is a certain deviation at the beginning, which can be corrected using formula (8) ; is a smoothing term, and its value is 0.000001 to prevent the denominator from being 0; the equations (4) to (8) are respectively: ; ; ; ; ; Iterative calculation is repeated by performing (3) to (5) to obtain the inversion result of the optimal transfer function method for elastic wave full waveform inversion.

[0007] The beneficial effects of the present invention compared with the prior art are as follows: As a robust distance measurement method, the optimal transfer function can effectively characterize the distribution difference between the observed data and the simulated data. This method jointly optimizes the distribution of P-wave and S-wave velocities and densities to reduce the optimal distance between the observed data and the simulated data. Compared with the traditional L2 norm method, it can better handle data offset and shape differences, making the inversion results less dependent on the initial model, while reducing the local minimum traps in the inversion process and improving the stability and accuracy of the inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 Elastic wave equation inversion flow chart; Figure 2 Elastic wave full waveform inversion calculation process for optimal transfer function; Figure 3The Marmousi model diagram, where (a) is the true model of longitudinal wave velocity, and (b) is the true model of shear wave velocity; Figure 4 It is a linear initial velocity model diagram, where (a) is the original model of longitudinal wave velocity, and (b) is the original model of shear wave velocity; Figure 5 is the earthquake sub-wave map; Figure 6 The elastic wave full waveform inversion result diagram of the traditional L2 objective function method, where (a) the true model of P-wave velocity, (b) the true model of S-wave velocity, (c) the original model of P-wave velocity, (d) the original model of S-wave velocity, (e) the inversion result of P-wave velocity, and (f) the inversion result of S-wave velocity; Figure 7 The results of elastic wave full waveform inversion using the optimal transfer function method, including (a) the true model of P-wave velocity, (b) the true model of S-wave velocity, (c) the original model of P-wave velocity, (d) the original model of S-wave velocity, (e) the inversion result of P-wave velocity, and (f) the inversion result of S-wave velocity; Figure 8 Comparison of the extracted traces of elastic wave full waveform inversion using the traditional L2 objective function method; Fig. 9 Comparison of elastic wave full waveform inversion and extraction gathers of elastic waves using the optimal transfer function method. DETAILED DESCRIPTION

[0009] The technical solution of the present invention is further explained below through embodiments in combination with the accompanying drawings, but the protection scope of the present invention is not limited in any form by the embodiments.

[0010] In the conventional elastic wave full waveform inversion, the present invention selects the optimal transmission function as the inversion target functional, derives the target functional, obtains the gradient of the optimal transmission, and finally adopts the gradient optimization algorithm to iteratively update the model. Experiments show that the results of elastic wave full waveform inversion based on the optimal transmission function are more accurate than the traditional elastic wave full waveform inversion based on the L2 norm.

[0011] Example 1 The present invention takes the Marmousi model as an example (such as Figure 3 The model is illustrated in Figure 2. The horizontal length of the model is 9200m, the vertical depth is 3400m, and the grid step size in the x and z directions of the model is 20m.

[0012] Based on this model, a total of 153 blastings were conducted, with the blasting points evenly distributed above the model at intervals of 60 m. 153 detection points were set, with the detection points also evenly distributed above the model at intervals of 60 m. The blasting point depth was 20 m, and the detection point depth was 500 m.

[0013] The specific implementation of the present invention is described in detail below. Figure 2 .

[0014] Step 1. Seismic record acquisition: The horizontal layered model with a horizontal length of 9200m and a vertical depth of 3440m is gridded. The grid step size in the x and z directions of the model is 20m, and a true velocity model with a grid size of 461×173 is obtained. The model contains model parameters Vp(x, z), Vs(x, z), ρ(x, z), x and z are expressed as spatial coordinates, x=[1, 2, 3, …, 461], z=[1, 2, 3, …, 173]. The shot points are evenly placed on the top of the model, with a total of 153 shots. The receiver points are also evenly placed on the top of the model, with a total of 154 receiver points. The shot points are excited from left to right on the top of the model. Each time they are excited, all the receiver points receive them simultaneously, and a total of 153 shots are obtained. Actual seismic records. Simulated earthquake records of 153 guns , select a linear model as the initial velocity model, such as Figure 4 , given the seismic wavelet W(t), which is a 10 Hz Ricker wavelet, such as Figure 5 .

[0015] Step 2. Selection of objective function: Introduce the optimal transfer function in the mathematical field as the objective function of inversion, based on actual seismic records Simulated earthquake records Establish the target functional (1), where To simulate earthquake records, For actual earthquake records, , for , The probability density function obtained after the preprocessing operator is: ; Step 2. Wave field calculation: Based on the seismic wavelet and initial velocity model, the elastic wave equation (2) is used for forward modeling, where is the displacement, is the stress tensor, is the source term. The positive wave field extension starts from the source and calculates different parameters The forward wave field under the condition of , , Wavefield data; using actual earthquake records and simulated earthquake records The difference between the two is taken as the earthquake source, and the elastic wave equation (2) is used to perform reverse time wave field extension to obtain Reverse time wave field distribution of parameters , and , the formula (2) is: ; Step 3. Gradient calculation: Derivative the optimal transfer function objective functional formula (1) to obtain the gradient calculation formula as shown in formula (3). Substitute the forward and reverse wave fields obtained in step (3) to calculate the model gradient, where the parameters are are model parameters, including The value of. The formula (3) is: ; Step 4. Update model parameters. Update model parameters using gradients according to the gradient optimization algorithm The optimization algorithm used here is the Adam gradient optimization algorithm, which uses formula (4) to perform gradient optimization to update the model, where k is the number of iterations. is the initial learning rate. The appropriate learning rate value can be selected based on multiple experiments. When the objective function is the optimal transfer function, The value range is 20~25. and are the first-order and second-order momentum error correction terms, respectively, is the momentum-corrected gradient. Let the gradient obtained by the gradient formula in step (4) be , and then calculate each The first moment estimate of In the formula, is a decay factor less than 1, usually around 0.9. is the current iteration number. It is an exponential moving average estimate of the first-order momentum, which will produce bias in the initial stage, so the unbiased estimate is obtained by using formula (6) through bias correction. ,and Can improve the stability and convergence speed of gradient descent. is the value after the deviation correction of the second-order moment estimate in the Adam optimization algorithm, which represents the variance estimate of the current gradient. , use formula (7) to calculate the first-order displacement average of its second-order moment, is the attenuation factor of the second-order moment, usually taken as 0.999, is the current iteration number. There is a certain deviation at the beginning, which can be corrected using formula (8) . is a smoothing term, which is generally taken as 0.000001 to prevent the denominator from being 0. The formulas (4) to (8) are respectively: ; ; ; ; ; Step 5. Iterative calculation. Repeat (3) to (5) to obtain the inversion result of the optimal transfer function method for elastic wave full waveform inversion.

[0016] In order to demonstrate the beneficial effects of the method of the present invention, the elastic wave full waveform inversion results of the traditional 𝐿2 objective function are compared here. Figure 6 is the elastic wave full waveform inversion result of the traditional 𝐿2 objective function, Figure 7 The elastic wave full waveform inversion result of the optimal transfer function method. By comparing the two figures, we can find that Figure 7 The continuity of the event axis of the inversion model has been enhanced, the deep structure has become clearer, and the overall accuracy has been greatly improved. Figure 8 Comparison of the traditional 𝐿2 objective function of elastic wave full waveform inversion with extracted traces. Fig. 9 Comparison of elastic wave full waveform inversion gathers using the optimal transfer function method shows that: Fig. 9 The inversion results are more consistent with the real model and the inversion effect is better.

Claims

1. A full waveform inversion method for elastic waves based on an optimal transfer function, characterized in that: The specific steps of the method are as follows: Step 1. Seismic record acquisition: firstly, an observation system is established and an exciting source of a certain frequency and a velocity model are selected. By exciting the source, actual seismic records and simulated seismic records are received at the detection point; Step 2. Selection of objective function: an optimal transmission function in the mathematical field is introduced as the objective function of inversion, and the overall shape and manifold structure of the data distribution are considered, and a target functional formula is established based on actual seismic records and simulated seismic records; Step 3. Wave field acquisition: forward modeling and reverse time extension are used to obtain the forward wave field and reverse time wave field; Step 4. Gradient acquisition: the target functional formula is differentiated to obtain the gradient calculation formula, which is brought into the forward wave field and reverse time wave field to obtain the gradient; Step 5. Model update: the required gradient is optimized by using the Adam gradient optimization algorithm and then brought into the iterative inversion to update the model; the above steps 3-5 are repeated until an inversion result that meets the accuracy requirements is obtained.

2. A full waveform inversion method for elastic waves based on an optimal transfer function, characterized in that: The method specifically comprises the following steps: Step 1. Seismic record acquisition: A linear model is selected as the initial model, and then a velocity model is given as the true model. The model contains model parameters Vp(x, z), Vs(x, z), ρ(x, z), x and z represent the coordinates of the space, x=[1,2, 3, …, Nx], z=[1, 2, 3, …, Nz], Nx and Nz represent the number of vertical and horizontal grid points of the velocity model respectively. The shot points are evenly placed on the top of the model, with a total of Nx shots. The detection points are also evenly placed on the top of the model, with a total of Nx detection points. The shot points are excited from left to right at the top of the model. Each time they are excited, all the detection points receive at the same time, and the actual seismic records with a total of Nx shots are obtained. Simulated earthquake records of Nx guns ; Step 2. Selection of objective function: Introduce the optimal transfer function in the mathematical field as the objective function of inversion, based on actual seismic records Simulated earthquake records Establish the target functional (1), where To simulate earthquake records, For actual earthquake records, , for , The probability density function obtained after the preprocessing operator is: is the objective function of elastic wave inversion; simulated data is measured by Wasserstein distance and The degree of matching; is the probability density function and of Suboptimal transmission distance; Indicates the number of receiving points; Indicates the receiving point index; is the wave field data of longitudinal wave velocity, Wave field data of shear wave velocity, Wavefield data of density distribution; is the reverse time wave field distribution of the longitudinal wave velocity, is the reverse time wave field distribution of shear wave velocity, is the reverse time wave field distribution of density distribution; the formula (1): ; Step 3. Wave field calculation: Based on the seismic wavelet and initial velocity model, the elastic wave equation (2) is used for forward modeling, where is the displacement, is the stress tensor, is the source term, is the medium density, For time, represents the differential operator, Represents the divergence of the stress tensor; the positive wave field extension starts from the earthquake source and calculates different parameters The forward wave field under the condition of , , Wavefield data; using actual earthquake records and simulated earthquake records The difference between the two is taken as the earthquake source, and the elastic wave equation (2) is used to perform reverse time wave field extension, and the following are obtained: Reverse time wave field distribution of parameters , and , the formula (2) is: ; Step 4. Gradient calculation: Derivative the optimal transfer function objective functional formula (1) to obtain the gradient calculation formula as shown in formula (3). Substitute the forward and reverse wave fields obtained in step (3) to calculate the model gradient, where the parameters are are model parameters, including The value of; the formula (3) is: ； Step 5. Update model parameters: Update model parameters using gradients according to the gradient optimization algorithm The optimization algorithm is the Adam gradient optimization algorithm, which performs gradient optimization through formula (4) to update the model, where k is the number of iterations. is the initial learning rate, and are the first-order and second-order momentum error correction terms, respectively, is the momentum-corrected gradient; the gradient obtained by the gradient formula in step (4) is , and then calculate each The first moment estimate of In the formula, is a decay factor less than 1, is the current iteration number; use formula (6) to obtain an unbiased estimate through bias correction ; is the value after the deviation correction of the second-order moment estimate in the Adam optimization algorithm, which represents the variance estimate of the current gradient; by , use formula (7) to calculate the first-order displacement average of its second-order moment, is the attenuation factor of the second-order moment, is the current iteration number; since There is a certain deviation at the beginning, which can be corrected using formula (8) ; is a smoothing term, and its value is 0.000001 to prevent the denominator from being 0; the equations (4) to (8) are respectively: ; ; ; ; ; Iterative calculation is performed, and steps (3) to (5) are repeated to obtain the inversion result of the optimal transfer function method for elastic wave full waveform inversion.

Images