A method for full waveform inversion of elastic waves based on stochastic gradient sampling

Abstract

This invention relates to a method for elastic wave full waveform inversion based on stochastic gradient sampling, comprising the following steps: Step 1) preprocessing the original seismic data before inversion; Step 2) performing Fourier transform on the preprocessed seismic data to obtain the true seismic record in the frequency domain; Step 3) establishing and initializing the elastic parameter model, setting the inversion parameters of the elastic wave full waveform and defining the observation system; Step 4) using the stochastic gradient sampling elastic wave full waveform inversion method, updating the elastic parameter model based on the true seismic record in the frequency domain to obtain the final elastic parameter model, thus completing the inversion. Compared with existing technologies, this invention can effectively perform elastic wave full waveform inversion even with poor initial elastic parameter models or a lack of low-frequency seismic data, and to a certain extent reduces the adverse effects of period jump phenomena on the inversion results. It also has advantages such as simple algorithm implementation and no additional memory consumption or computational load.

Description

Technical Field

[0001] This invention relates to the field of velocity modeling in exploration seismology, and in particular to a method for full waveform inversion of elastic waves based on stochastic gradient sampling. Background Technology

[0002] Elastic wave full-waveform inversion can utilize multi-component elastic wave seismic data to obtain subsurface elastic parameters. With the rapid development of computer technology and seismic acquisition techniques in recent years, the significant computational and memory consumption issues associated with elastic wave full-waveform inversion have been largely resolved, leading to its successful application to large amounts of real-world data. However, while elastic wave full-waveform inversion yields rich information on subsurface physical parameters, it also faces significant challenges. Because full-waveform inversion is a highly nonlinear and nonconvex optimization problem, it is prone to getting trapped in local extrema. Lower shear wave velocities exacerbate the nonlinearity of the elastic wave full-waveform inversion problem; compared to P-wave velocities, inverting shear wave velocities requires constructing a more accurate initial model to avoid getting trapped in local extrema. Furthermore, the lack of low-frequency seismic data is even more detrimental to reconstructing the shear wave velocity model, making the "period jump" phenomenon more pronounced. On the other hand, due to the complex propagation mechanism of elastic waves in the subsurface medium, phenomena such as diffraction and converted waves occur, and parameter coupling effects can also have a very negative impact on the inversion results.

[0003] To address the parameter coupling effect in the full waveform inversion of elastic waves, numerous scholars have proposed a series of solutions. These include multi-step inversion strategies to gradually recover P-wave velocity (impedance), S-wave velocity (impedance), and density; appropriate selection of parameterization methods; and P / S wavefield separation. These measures can effectively eliminate aliasing interference from P-wave velocity disturbances in the inversion results, achieving parameter decoupling. Xu and McMechan proposed a multi-step inversion method to constrain the updates of parameters with varying strengths, which can also partially suppress the parameter coupling effect.

[0004] However, research on the "period jump" problem is still limited. While multi-scale inversion strategies (from low-frequency data to high-frequency data) can prioritize obtaining background information and avoid premature recovery of high-wavenumber information, they do not fundamentally solve the problem of missing low-frequency data. Other methods, such as using travel-time tomography to provide a better initial background model, or using seismic waveform envelope and phase information for inversion, as well as using subsurface scattering angle filtering to obtain parameter background and optimal transport distance elastic wave inversion, can also address the "period jump" problem to some extent, improve inversion results, and obtain parameter models with higher accuracy and resolution.

[0005] However, all of the above methods change the objective functional of the original elastic wave full waveform inversion, increasing the algorithm complexity. Furthermore, some algorithms involve a large amount of additional computation and memory consumption during implementation, hindering their practical application. Summary of the Invention

[0006] The purpose of this invention is to provide a method for elastic wave full waveform inversion based on stochastic gradient sampling, which can effectively perform elastic wave full waveform inversion under poor initial elastic parameter models or lack of low-frequency seismic data, reduce the adverse effects of period jump phenomenon on the inversion results, simplify algorithm complexity, and improve calculation speed.

[0007] The objective of this invention can be achieved through the following technical solutions:

[0008] A method for full waveform inversion of elastic waves based on stochastic gradient sampling includes the following steps:

[0009] Step 1) Perform preprocessing on the raw seismic data before inversion;

[0010] Step 2) Perform Fourier transform on the preprocessed seismic data to obtain the true seismic record in the frequency domain;

[0011] Step 3) Establish and initialize the elastic parameter model, set the inversion parameters of the full elastic wave waveform and define the observation system;

[0012] Step 4) The full waveform inversion method of elastic waves using stochastic gradient sampling is adopted. The elastic parameter model is updated based on the real seismic records in the frequency domain to obtain the final elastic parameter model and complete the inversion.

[0013] The elastic parameters include longitudinal wave velocity, transverse wave velocity, and density.

[0014] The inversion parameters include the inversion frequency.

[0015] Step 4) includes the following steps:

[0016] Step 41) Execute the frequency group loop;

[0017] Step 42) Within each frequency group loop, execute an iterative loop, performing a forward propagation elastic wavefield simulation at the shot point for each shot under the elastic parameter model updated in the previous iteration. Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records. With the simulated wave field {u(x,z,ω)=(u x ,u z ) T} n ,(n=1,2,3…Ns), where g represents the spatial position of the detector, ω is the angular frequency, and Ns is the total number of shots;

[0018] Step 43) From the frequency domain real seismic record d obs (g,ω) and simulated seismic record d cal(g,ω) Calculate the residual seismic record res(g,ω)=(res x ,res z ) T The residuals between real and simulated seismic records under the L2 norm are calculated as the objective function, where... This is the frequency domain real seismic record obtained in step 2);

[0019] Step 44) Determine whether the inversion meets the convergence or stopping iteration conditions based on the objective function. If it does, exit the current frequency group loop and return to step 41) to enter the next frequency group loop until the final inversion result is obtained; if it does not meet the conditions, proceed to step 45).

[0020] Step 45) For each shot, the residual record is used as the accompanying source for the elastic wave field simulation. The elastic wave backpropagation wave field simulation at the receiver is performed. The Fourier transform of the simulated surface seismic record is then used to obtain the frequency domain simulated backpropagation seismic record.

[0021] Step 46) Calculate the gradient of the objective function with respect to each elastic parameter using a calculation method based on stochastic gradient sampling;

[0022] Step 47) Update the elastic parameter model using the conjugate gradient method;

[0023] Step 48) Return to step 42) and proceed to the next iteration.

[0024] Step 42) specifically includes the following steps:

[0025] Step 421) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary, and perform a shot-point end elastic wave propagation simulation for each shot, satisfying the elastic wave displacement stress equation in an isotropic medium:

[0026]

[0027] Among them, u x with u z The X and Z components of the particle displacement; σ xx , σ xz , σ zx , σ zz λ is the stress component; μ is the Lamé coefficient, λ = (α 2 -2β 2 )ρ,μ=β 2 ρ, α, β, and ρ represent the longitudinal wave velocity, transverse wave velocity, and density, respectively; f x with f z The X and Z components of the body wave source;

[0028] Step 422) Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records {d}. cal (g,ω)} n With the simulated wave field {u(x,z,ω)} n ,(n=1,2,3…Ns).

[0029] Step 43) specifically includes the following steps:

[0030] Step 431) Calculate the residual seismic record:

[0031] res(g,ω)=d cal (g,ω)-d obs (g,ω)

[0032] Where, d cal (g,ω),d obs (g,ω) represent the frequency domain simulated seismic record and the real seismic record received by the nth shot at spatial location g, respectively;

[0033] Step 432) Calculate the residuals between the real and simulated seismic records under the L2 norm as the objective function:

[0034]

[0035] Among them, res x,n (g,ω),res z,n (g,ω) represent the X and Z components of the residual seismic record of the nth shot at spatial position g, respectively; Nr is the total number of receiver points for the nth shot; For conjugate.

[0036] Step 45) specifically includes the following steps:

[0037] Step 451) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary. For each shot, use the residual record as the accompanying source for the elastic wave field simulation, and perform elastic wave back propagation simulation at the receiver end, satisfying:

[0038]

[0039] Among them, u x with u z The X and Z components of the particle displacement; σ xx , σ xz , σ zx , σ zz λ is the stress component; μ is the Lamé coefficient, λ = (α 2 -2β 2 )ρ,μ=β 2ρ, α, β, and ρ represent the longitudinal wave velocity, transverse wave velocity, and density, respectively; res x with res z The X and Z components of the accompanying seismic source are the X and Z components of the residual seismic record.

[0040] Step 452) Perform a Fourier transform on the simulated surface seismic record to obtain the frequency domain simulated backpropagation seismic record {u * (x,z,ω)} n ,(n=1,2,3…Ns).

[0041] Step 46) specifically includes the following steps:

[0042] Step 461) Calculate the spatial derivative of the simulated seismic wave field using the fourth-order spatial difference. in To simulate seismic wave fields;

[0043] Step 462) Define the range of random spatial movement:

[0044] h αx ∈[-k1λ α ,k1λ α ]

[0045] h αz ∈[-k2λ α ,k2λ α ]

[0046] h βx ∈[-k3λ β ,k3λ β ]

[0047] h βz ∈[-k4λ β ,k4λ β ]

[0048] λ α ,λ β Given the longitudinal wave wavelength and the transverse wave wavelength, satisfying:

[0049]

[0050] Where α and β are the P-wave velocity and S-wave velocity, respectively, and f0 = ω / 2π is the inversion frequency; k i (i = 1, 2, 3, 4 are proportionality coefficients;)

[0051] Step 463) Randomly select h αx ,h αz ,h βx ,h βzRandom spatial shift of the spatial derivative of the simulated seismic wave field:

[0052]

[0053]

[0054] Step 464) Calculate the gradient of the elastic wave objective function with respect to the P-wave velocity and S-wave velocity using the seismic wave field after random spatial movement:

[0055]

[0056]

[0057] Among them, i=x,z,j=x,z,l=x,z,m=x,z,δ ij It is a pulse function;

[0058] Step 465) Apply gradient g RSS (α) and g RSS (β) Perform Gaussian smoothing.

[0059] Step 47) specifically includes the following steps:

[0060] Step 471) Calculate the conjugate gradient:

[0061]

[0062] y k (v)=g k (v)-g k-1 (v)

[0063]

[0064]

[0065] Where k is the iteration number, v = {α, β}, g k (v)={g k (α),g k (β)} represents the gradient of the longitudinal wave velocity and the transverse wave velocity during the k-th iteration;

[0066] Step 472) Determine the update step size:

[0067]

[0068] Where ε is the scaling factor for updating the model;

[0069] Step 473) Update the elasticity parameter model:

[0070] v k+1 =vk +λ k (v)z k (v)

[0071] The preprocessing of the raw seismic data before inversion includes noise reduction filtering.

[0072] Compared with the prior art, the present invention has the following beneficial effects:

[0073] (1) Reduce the influence of period jump phenomenon: The present invention can effectively perform full waveform inversion of elastic waves under poor initial elastic parameter model or lack of low frequency seismic data, reduce the adverse effects of period jump phenomenon on inversion results, and improve inversion accuracy.

[0074] (2) No additional memory consumption and computation: This invention uses the random spatial movement of a single reference model wave field to approximate the calculation of several random sampled model wave fields. There is no need to calculate the gradients of multiple random models to expand the search space, nor is there a need to store the gradients of multiple random sampled models. This improves the calculation speed and reduces memory consumption.

[0075] (3) Simple algorithm: Under the framework of traditional elastic wave full waveform inversion program, this invention only needs to add a random spatial movement module of wave field spatial derivative, and use the wave field spatial derivative after random spatial movement to calculate the gradient. Other contents remain basically unchanged, and the algorithm is simple and easy to implement. Attached Figure Description

[0076] Figure 1 This is a flowchart of the method of the present invention;

[0077] Figure 2 The model is the actual elastic parameter model in Example 1, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0078] Figure 3 The first-generation update direction (negative gradient) for conventional elastic wave full waveform inversion in Example 1 is shown, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0079] Figure 4 The results are the conventional full waveform inversion results of elastic waves in Example 1, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0080] Figure 5 The first-generation update direction (negative gradient) for the stochastic gradient elastic wave full waveform inversion proposed in Example 1 is shown in this invention, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0081] Figure 6The results of the full waveform inversion of the stochastic gradient elastic wave proposed in this invention in Example 1 are shown, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0082] Figure 7 The results of the full waveform inversion of stochastic gradient elastic waves and conventional elastic waves proposed in this invention in Example 1 are shown, where (a) represents the longitudinal wave velocity and (b) represents the transverse wave velocity.

[0083] Figure 8 This is the actual shear wave velocity model in Example 2;

[0084] Figure 9 The first-generation update direction (negative gradient) for the conventional elastic wave full waveform inversion in Example 2;

[0085] Figure 10 This is the conventional full waveform inversion result of elastic waves in Example 2;

[0086] Figure 11 The first-generation update direction (negative gradient) for the stochastic gradient elastic wave full waveform inversion proposed in Example 2 of this invention;

[0087] Figure 12 The result of the full waveform inversion of the stochastic gradient elastic wave proposed in this invention in Example 2;

[0088] Figure 13 The result of the random gradient elastic wave full waveform inversion—conventional elastic wave full waveform inversion proposed in this invention in Example 2 is shown. Detailed Implementation

[0089] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. These embodiments are based on the technical solution of the present invention and provide detailed implementation methods and specific operating procedures. However, the scope of protection of the present invention is not limited to the following embodiments.

[0090] In some complex underground structures with large velocity variations, especially strong low-velocity anomalies, complex wave phenomena, including complex converted waves, can occur. Furthermore, the velocity of shear waves is lower than that of P-waves, which exacerbates the nonlinearity of elastic wave full waveform inversion and increases the difficulty of the inversion. In actual exploration areas, due to the absorption and attenuation of low-frequency seismic data and the limitation of technology, detectors cannot receive extremely low-frequency signals. The lack of low-frequency seismic data easily leads to periodic jumps in nonlinear inversion, causing it to fall into local extrema and severely affecting the accuracy of the inversion results. Conventional elastic wave full waveform inversion is difficult to adapt to these problems. Therefore, overcoming periodic jumps and avoiding falling into local extrema in full waveform inversion is an urgent problem to be solved. This invention proposes an elastic wave full waveform inversion method based on stochastic gradient sampling, such as... Figure 1As shown, this method can alleviate the periodic jump phenomenon to a certain extent and improve its adverse effects. It has more advantages than conventional elastic wave full waveform inversion. The specific steps include:

[0091] Step 1) Perform preprocessing on the raw seismic data before inversion, including but not limited to noise reduction filtering.

[0092] Step 2) Perform Fourier transform on the preprocessed seismic data to obtain the true frequency domain seismic record. Where g represents the spatial position of the detector, ω is the angular frequency, and Ns is the total number of shots.

[0093] Step 3) Establish and initialize the elastic parameter model, set the inversion parameters of the full elastic wave waveform and define the observation system. The elastic parameters include the longitudinal wave velocity, transverse wave velocity and density, and the inversion parameters include the inversion frequency.

[0094] Step 4) The full waveform inversion method of elastic waves using stochastic gradient sampling is adopted. The elastic parameter model is updated based on the real seismic records in the frequency domain to obtain the final elastic parameter model and complete the inversion.

[0095] Step 41) Execute the frequency group loop;

[0096] Step 42) Within each frequency group loop, execute an iterative loop, performing a forward propagation elastic wavefield simulation at the shot point for each shot under the elastic parameter model updated in the previous iteration. Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records. With the simulated wave field {u(x,z,ω)=(u x ,u z ) T} n (n = 1, 2, 3…Ns),

[0097] Step 43) From the frequency domain real seismic record d obs (g,ω) and simulated seismic record d cal (g,ω) Calculate the residual seismic record res(g,ω)=(res x ,res z ) T The residuals between real and simulated seismic records under the L2 norm are calculated as the objective function.

[0098] Step 44) Determine whether the inversion meets the convergence or stopping iteration conditions based on the objective function. If it does, exit the current frequency group loop and return to step 41) to enter the next frequency group loop until the final inversion result is obtained; if it does not meet the conditions, proceed to step 45).

[0099] Step 45) For each shot, the residual record is used as the accompanying source for the elastic wave field simulation. The elastic wave backpropagation wave field simulation at the receiver is performed. The Fourier transform of the simulated surface seismic record is then used to obtain the frequency domain simulated backpropagation seismic record.

[0100] Step 46) Calculate the gradient of the objective function with respect to each elastic parameter using a calculation method based on stochastic gradient sampling;

[0101] Step 47) Update the elastic parameter model using the conjugate gradient method;

[0102] Step 48) Return to step 42) and proceed to the next iteration.

[0103] Step 42) specifically includes the following steps:

[0104] Step 421) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary, and perform a shot-point end elastic wave propagation simulation for each shot, satisfying the elastic wave displacement stress equation in an isotropic medium:

[0105]

[0106] Among them, u x with u z The X and Z components of the particle displacement; σ xx , σ xz , σ zx , σ zz λ is the stress component; μ is the Lamé coefficient, λ = (α 2 -2β 2 )ρ,μ=β 2 ρ, α, β, and ρ represent the longitudinal wave velocity, transverse wave velocity, and density, respectively; f x with f z The X and Z components of the body wave source;

[0107] Step 422) Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records {d}. cal (g,ω)} b With the simulated wave field {u(x,z,ω)} n ,(n=1,2,3…Ns).

[0108] Step 43) specifically includes the following steps:

[0109] Step 431) Calculate the residual seismic record:

[0110] res(g,ω)=d cal (g,ω)-d obs (g,ω)

[0111] Where, d cal (g,ω),d obs (g,ω) represent the frequency domain simulated seismic record and the real seismic record received by the nth shot at spatial location g, respectively;

[0112] Step 432) Calculate the residuals between the real and simulated seismic records under the L2 norm as the objective function:

[0113]

[0114] Among them, res x,n (g,ω),res z,n (g,ω) represent the X and Z components of the residual seismic record of the nth shot at spatial position g, respectively; Nr is the total number of receiver points for the nth shot; For conjugate.

[0115] Step 45) specifically includes the following steps:

[0116] Step 451) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary. For each shot, use the residual record as the accompanying source for the elastic wave field simulation, and perform elastic wave back propagation simulation at the receiver end, satisfying:

[0117]

[0118] Among them, u x with u z The X and Z components of the particle displacement; σ xx , σ xz , σ zx , σ zz λ is the stress component; μ is the Lamé coefficient, λ = (α 2 -2β 2 )ρ,μ=β 2 ρ, α, β, and ρ represent the longitudinal wave velocity, transverse wave velocity, and density, respectively; res x with res z The X and Z components of the accompanying seismic source are the X and Z components of the residual seismic record.

[0119] Step 452) Perform a Fourier transform on the simulated surface seismic record to obtain the frequency domain simulated backpropagation seismic record {u * (x,z,ω)} n ,(n=1,2,3…Ns).

[0120] Step 46) specifically includes the following steps:

[0121] Step 461) Calculate the spatial derivative of the simulated seismic wave field using the fourth-order spatial difference. in To simulate seismic wave fields;

[0122] Step 462) Define the range of random spatial movement:

[0123] h αx ∈[-k1λ α ,k1λ α ]

[0124] h αz ∈[-k2λ α ,k2λ α ]

[0125] h βx ∈[-k3λ β ,k3λ β ]

[0126] h βz ∈[-k4λ β ,k4λ β ]

[0127] λ α ,λ β Given the longitudinal wave wavelength and the transverse wave wavelength, satisfying:

[0128]

[0129] Where α and β are the P-wave velocity and S-wave velocity, respectively, and f0 = ω / 2π is the inversion frequency; k i (i = 1, 2, 3, 4 are proportionality coefficients;)

[0130] Step 463) Randomly select h αx ,h αz ,h βx ,h βz Random spatial shift of the spatial derivative of the simulated seismic wave field:

[0131]

[0132]

[0133] Step 464) Calculate the gradient of the elastic wave objective function with respect to the P-wave velocity and S-wave velocity using the seismic wave field after random spatial movement:

[0134]

[0135]

[0136] Among them, i=x,z,j=x,z,l=x,z,m=x,z,δ ij It is a pulse function;

[0137] Step 465) Apply gradient g RSS (α) and g RSS (β) Perform 11-point Gaussian window smoothing.

[0138] Step 47) specifically includes the following steps:

[0139] Step 471) Calculate the conjugate gradient:

[0140]

[0141] y k (v)=g k (v)-g k-1 (v)

[0142]

[0143]

[0144] Where k is the iteration number, v = {α, β}, g k (v)={g k (α),g k (β)} represents the gradient of the longitudinal wave velocity and the transverse wave velocity during the k-th iteration;

[0145] Step 472) Determine the update step size:

[0146]

[0147] Wherein, ε is the scaling factor for updating the model, and in this embodiment, ε is taken as 1%;

[0148] Step 473) Update the elasticity parameter model:

[0149] v k+1 =v k +λ k (v)z k (v)

[0150] Example 1

[0151] This embodiment uses the above-mentioned elastic wave full waveform inversion method based on stochastic gradient sampling as the real model, employing a heterogeneous model of low-velocity P- and S-wave velocity anomalies in a two-dimensional Gaussian sphere as the true model (e.g., Figure 2As shown in the figure, the density is 2000 kg / m³ and remains constant. The model has 401 × 201 grids with a grid spacing of 10 m × 10 m. The maximum and minimum P-wave velocities are 3000 m / s and 2100 m / s, respectively, and the maximum and minimum S-wave velocities are 2400 m / s and 1600 m / s, respectively. A two-dimensional isotropic elastic wave forward modeling simulation was performed on this model, with a total of 48 shots. Shot points were uniformly distributed on the surface with a shot spacing of 80 m, and the first shot was at a horizontal position of 80 m. Receiver points were uniformly distributed at a depth of 2000 m underground at 10 m intervals. The excitation source was a 10 Hz dominant frequency Ricker wavelet (containing both P-wave and S-wave energy), and the seismic record reception length was 2.0 seconds with an interval of 1 ms. Elastic wave full waveform inversion was performed using 5-7 Hz frequency domain data. The inversion uses a uniform velocity model as the initial model, with P-wave and S-wave velocities of 3000 m / s and 2400 m / s, respectively. To avoid surface model updates caused by shot point excitation, the actual P-wave and S-wave velocities are maintained and kept constant above a depth of 50 meters in the subsurface. After applying the elastic wave full waveform inversion based on stochastic gradient sampling (RSS-EFWI) proposed in this invention, the inversion results are used as model initialization parameters, and conventional elastic wave full waveform inversion (RSS-EFWI-CEFWI) is performed. Simultaneously, conventional elastic wave full waveform inversion (CEFWI) is also applied to highlight the effectiveness and superiority of this invention through comparison.

[0152] The first-generation model update directions (negative gradients) calculated using CEFWI and RSS-EFWI are as follows: Figure 3 and 5 As shown. According to the elastic parameters of the real model, near the velocity anomaly centers on the left and right sides, updates should logically be made in the negative velocity direction. However, CEFWI updates positively at these locations, which is unexpected and indicates that the inversion is heading in the wrong direction, resulting in periodic jumps. Conversely, RSS-EFWI updates negatively at these locations, which is expected and makes the inverted model closer to the real model. Figure 4 The CEFWI inversion results exhibited periodic jumps, falling into local extrema. The results showed strong velocity discontinuities, deviating significantly from the actual model, making further seismic interpretation difficult. Figure 6 The displayed RSS-EFWI results better match the background characteristics of the real model, providing a better initial model reference for full waveform inversion. The final RSS-EFWI-CEFWI inversion results are as follows: Figure 7As shown, both the P-wave velocity and the S-wave velocity are very close to the actual model, and no period jump phenomenon is observed, resulting in a high-precision, high-resolution inversion model. This verifies that, compared with conventional elastic wave full-waveform inversion, the present invention is better suited to poor initial models and missing low-frequency data, and to a certain extent improves the adverse effects of period jump phenomenon on the inversion results.

[0153] Example 2

[0154] This embodiment applies the stochastic gradient sampling elastic wave full waveform inversion method proposed in this invention to a reflection observation system, specifically to a horizontally layered medium model (such as...). Figure 8 In the transverse wave velocity inversion shown (as indicated), the longitudinal wave velocity is 3000 m / s, and the density is 2000 kg / m³, both remaining constant. The model has 301 × 151 grids with a grid spacing of 10 m × 10 m. The maximum and minimum transverse wave velocities are 2000 m / s and 1200 m / s, respectively. Two-dimensional isotropic elastic wave forward modeling was performed on this model, with a total of 59 shots. Shot points were uniformly distributed on the surface, with a shot spacing of 50 m, and the first shot was at a horizontal position of 50 m. Receiver points were uniformly distributed on the surface at 10-meter intervals. The excitation source was a 10 Hz dominant frequency Ricker wavelet (containing both longitudinal and transverse wave energy), and the seismic record reception length was 2.5 seconds with a 1 ms interval. Elastic wave full waveform inversion was performed using frequency domain data from 5–22.5 Hz. A uniform transverse wave velocity model with a magnitude of 2000 m / s was used as the initial model for the inversion. To avoid abnormal surface model updates caused by shot point activation, the true shear wave velocity (1500 m / s) was maintained and kept constant at a depth of 50 meters or more in the shallow subsurface. RSS-EFWI was applied to the low-frequency band (5-7.5 Hz frequency domain data), and then the inversion results were used as the initial model to perform RSS-EFWI-CEFWI on the full frequency domain data (5-22.5 Hz). Simultaneously, conventional elastic wave full waveform inversion (CEFWI) was also applied (5-7.5 Hz frequency domain data).

[0155] The first-generation model update directions (negative gradients) calculated using CEFWI and RSS-EFWI are as follows: Figure 9 and 11 As shown, while CEFWI updates in the negative velocity direction at shallow depths, several segments of positive updates occur at deeper depths. These gradient fluctuations and oscillations have a very detrimental effect on subsequent inversion. Especially above 800 meters in depth, the real model exhibits strong low-velocity anomalies, which are not reflected in the traditional CEFWI gradient. In contrast, RSS-EFWI updates in the negative direction, as expected, and its update range is deeper than CEFWI, aligning with the low-velocity anomaly range of the real model. Figure 10The displayed CEFWI inversion results do not conform to the true model in terms of either velocity values ​​or layer interface depth, indicating that it has fallen into a local extremum. Figure 12 The RSS-EFWI results shown can, to some extent, invert the approximate location of the horizontal layer, and also reflect the corresponding high and low velocity characteristics to some extent, making it a good initial model for subsequent full waveform inversion. The final RSS-EFWI-CEFWI inversion results are as follows: Figure 13 As shown, the differences between the shallow and middle layers and the actual model are minimal in terms of layer interface depth and layer velocity. Due to the lack of deep layer reflection data, there are only some errors in the numerical values ​​of the velocities in the middle and deep layers. However, the reflection of the layer interface is clearly characterized and its location is accurate. Therefore, for the shear wave velocity inversion work of the reflection observation system, this invention also performs well in the case of poor initial models and missing low-frequency data, and can obtain high-precision and high-resolution elastic parameter models.

[0156] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A method for full waveform inversion of elastic waves based on stochastic gradient sampling, characterized in that, Includes the following steps: Step 1) Perform preprocessing on the raw seismic data before inversion; Step 2) Perform Fourier transform on the preprocessed seismic data to obtain the true seismic record in the frequency domain; Step 3) Establish and initialize the elastic parameter model, set the inversion parameters of the full elastic wave waveform and define the observation system; Step 4) The elastic parameter model is updated based on the real seismic records in the frequency domain using the stochastic gradient sampling elastic wave full waveform inversion method to obtain the final elastic parameter model and complete the inversion. Step 4) includes the following steps: Step 41) Execute the frequency group loop; Step 42) Within each frequency group loop, execute an iterative loop, performing a forward propagation elastic wavefield simulation at the shot point for each shot under the elastic parameter model updated in the previous iteration. Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records. With simulated wave field ,in, Indicates the spatial position of the detector. Angular frequency, This represents the total number of guns; Step 43) From the real seismic records in the frequency domain Compared with simulated earthquake records Calculate residual seismic records The residuals between real and simulated seismic records under the L2 norm are calculated as the objective function, where... This is the frequency domain real seismic record obtained in step 2); Step 44) Determine whether the inversion meets the convergence or stopping iteration conditions based on the objective function. If it does, exit the current frequency group loop and return to step 41) to enter the next frequency group loop until the final inversion result is obtained; if it does not meet the conditions, proceed to step 45). Step 45) For each shot, using the residual record as the accompanying source for the elastic wave field simulation, perform the receiver-end elastic wave backpropagation wave field simulation, and perform a Fourier transform on the simulated surface seismic record to obtain the frequency domain simulated backpropagation seismic record. ; Step 46) Calculate the gradient of the objective function with respect to each elastic parameter using a calculation method based on stochastic gradient sampling; Step 47) Update the elastic parameter model using the conjugate gradient method; Step 47) specifically includes the following steps: Step 471) Calculate the conjugate gradient: in, For the number of iterations, , For the first The gradient of P-wave velocity and S-wave velocity in the next iteration process; Step 472) Determine the update step size: in, To update the scaling factor of the model; Step 473) Update the elasticity parameter model: ; Step 48) Return to step 42) and proceed to the next iteration.

2. The elastic wave full waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The elastic parameters include longitudinal wave velocity, transverse wave velocity, and density.

3. The elastic wave full waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, The inversion parameters include the inversion frequency.

4. The method for full waveform inversion of elastic waves based on stochastic gradient sampling according to claim 1, characterized in that, Step 42) specifically includes the following steps: Step 421) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary, and perform a shot-end elastic wave propagation simulation for each shot, satisfying the elastic wave displacement stress equation in an isotropic medium: in, and These are the X and Z components of the particle displacement; , , These are stress components; The Lamé coefficient is given. , These are the longitudinal wave velocity, the transverse wave velocity, and the density; and The X and Z components of the body wave source; Step 422) Perform Fourier transforms on the simulated surface seismic records and wavefields to obtain the frequency domain simulated seismic records. With simulated wave field .

5. The elastic wave full waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, Step 43) specifically includes the following steps: Step 431) Calculate the residual seismic record: in, The first Cannon in spatial position The received frequency domain simulated seismic records and real seismic records; Step 432) Calculate the residuals between the real and simulated seismic records under the L2 norm as the objective function: in, The first Cannon in spatial position The residual seismic records of X and Z components; For the first Total number of receiver points for the shot; For conjugate.

6. The method for full waveform inversion of elastic waves based on stochastic gradient sampling according to claim 1, characterized in that, Step 45) specifically includes the following steps: Step 451) Using the fourth-order spatial and second-order temporal finite difference method, implement CPML absorbing boundary conditions at the boundary. For each shot, use the residual record as the accompanying source for the elastic wave field simulation, and perform elastic wave back propagation simulation at the receiver end, satisfying: in, and These are the X and Z components of the particle displacement; , , These are stress components; The Lamé coefficient is given. , These are the longitudinal wave velocity, the transverse wave velocity, and the density; and The X and Z components of the accompanying seismic source are the X and Z components of the residual seismic record. Step 452) Perform a Fourier transform on the simulated surface seismic records to obtain the frequency domain simulated back-propagation seismic records. .

7. The elastic wave full waveform inversion method based on stochastic gradient sampling according to claim 1, characterized in that, Step 46) specifically includes the following steps: Step 461) Calculate the spatial derivative of the simulated seismic wave field using the fourth-order spatial difference. ,in To simulate seismic wave fields; Step 462) Define the range of random spatial movement: Given the longitudinal wave wavelength and the transverse wave wavelength, satisfying: in, These are the longitudinal wave velocity and the transverse wave velocity, respectively. For inversion frequency; This is the proportionality coefficient; Step 463) Random selection Random spatial shift of the spatial derivative of the simulated seismic wavefield: Step 464) Calculate the gradient of the elastic wave objective function with respect to the P-wave velocity and S-wave velocity using the seismic wave field after random spatial shift: in, , , , , It is a pulse function; Step 465) Gradient and Perform Gaussian smoothing.

8. The method for full waveform inversion of elastic waves based on stochastic gradient sampling according to claim 1, characterized in that, The preprocessing of the raw seismic data before inversion includes noise reduction filtering.

Images