A two-step method for inversion of anisotropic parameters of VTI media

Abstract

The present invention relates to a two-step method for inverting anisotropic parameters of VTI media, which belongs to the field of seismic velocity modeling. Vp , Vs. , r The gradients of the three model parameters are optimized using the conjugate gradient optimization algorithm and then brought into iterative inversion to obtain Vp , Vs. , r Three model parameters. Vp , Vs. , r These three model parameters are used as the real model parameters, and only e and d Get their gradients, and then use the Adam gradient algorithm to optimize their gradients and bring them into iterative inversion until the accuracy requirements are met, and finally get e and d The method of the present invention reduces the crosstalk effect of simultaneous inversion of multiple parameters, and selects different gradient optimization algorithms according to the characteristics of the parameters during the inversion process, thereby improving the accuracy of gradient optimization in a targeted manner and significantly improving the effect of full waveform inversion of anisotropic parameters of VTI media.

Description

Technical Field

[0001] The invention belongs to the field of seismic velocity modeling, and in particular relates to a VTI medium anisotropic parameter inversion method based on a two-step method. Background Art

[0002] The full waveform inversion method is a pre-stack data imaging and velocity modeling algorithm. It estimates underground physical parameters by minimizing the difference between simulated data and observed data through optimization iteration. It can fully utilize the travel time, amplitude, and phase information in the seismic data to reconstruct the velocity, density and other parameter models of the underground medium and obtain a high-quality inversion model. However, most conventional full waveform inversion methods are based on isotropic media such as acoustic waves or elastic waves for inversion, while the real underground media often have anisotropic characteristics. The full waveform inversion under the assumption of isotropic media does not conform to the actual situation. Therefore, the inversion of anisotropic parameters based on vertical transverse isotropy (VTI) media, which is more in line with the actual geological conditions, has attracted widespread attention in the industry and has become a research hotspot in the field of geophysics.

[0003] In the inversion of anisotropic parameters of VTI media, the current inversion method is to use the same gradient optimization algorithm to simultaneously invert the model parameters Vp, Vs, ρ and the anisotropic parameters ε and δ. However, the model parameters Vp, Vs, ρ and the anisotropic parameters ε and δ have a very large difference in magnitude, often more than 1000 times. During simultaneous inversion, the anisotropic parameters ε and δ are greatly affected by the crosstalk of the model parameters Vp, Vs, ρ, and are updated very slowly and easily fall into local extremes, resulting in a poor inversion effect. Summary of the invention

[0004] The present invention provides a two-step method for inverting anisotropic parameters of VTI media. The method reduces the crosstalk effect of simultaneous inversion of multiple parameters, and during the inversion process, different gradient optimization algorithms are selected according to the characteristics of the parameters, so as to improve the accuracy of gradient optimization in a targeted manner, thereby significantly improving the effect of full waveform inversion of anisotropic parameters of VTI media.

[0005] The present invention is achieved through the following technical solutions:

[0006] A method for inverting anisotropic parameters of VTI media based on a two-step method, the method specifically includes the following steps: the first step is to invert the gradients of the three model parameters Vp, Vs, and ρ in the full waveform inversion of the VTI medium, and then use the conjugate gradient optimization algorithm to optimize the gradients of Vp, Vs, and ρ, and use the optimized gradients to iteratively invert until the accuracy requirements are met (generally, 100 inversions are sufficient), and the three model parameters Vp, Vs, and ρ are obtained. The second step is to use the three model parameters Vp, Vs, and ρ obtained by the previous step as the true model parameters, and only invert the two anisotropic parameters ε and δ in the subsequent inversion process to obtain the gradients of ε and δ, and then use the Adam gradient algorithm in deep learning to optimize the gradients of ε and δ, and use the optimized gradients to iteratively invert until the accuracy requirements are met (generally, 100 inversions are sufficient), and finally obtain the two anisotropic parameters ε and δ.

[0007] The present invention also provides a specific VTI medium anisotropic parameter inversion method based on a two-step method, the method specifically comprising the following steps:

[0008] (1) Given the initial model parameters Vp(x, z), Vs(x, z), ρ(x, z) and anisotropy parameters ε(x, z) and δ(x, z), where x and z represent spatial position coordinates, x=1, 2, 3, …, Nx, z=1, 2, 3, …, Nz, and Nx and Nz represent the total number of model horizontal and vertical grid points, respectively. The shot points are evenly placed on the top of the model, with a total of Nx shots. The receiver points are also evenly placed on the top of the model, with a total of Nx 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 the signal simultaneously, and the actual seismic record Record of Nx shots is obtained. Based on the actual seismic record Record obtained, full waveform inversion is performed using the adjoint state method. During the inversion process, only Vp(x, z), Vs(x, z) and ρ(x, z) are inverted, and ε(x, z) and δ(x, z) remain unchanged. In each inversion, the gradient profile is obtained: gVp k (x, z), gVs k (x, z) and gρ k (x, z), where k is the number of inversions, k=1, 2, 3, …, K, and K is the total number of inversions;

[0009] (2) Using the conjugate gradient optimization algorithm formula (1), gVp k (x, z), gVs k (x, z) and gρ k The three gradient profiles (x, z) are optimized to obtain the optimized gradient GVp k (x, z), GVsk (x, z) and Gρ k (x, z), g in formula (1) k represents the gradient before processing, G k represents the processed gradient, T represents the transposition of the array, and the formula (1) is: ;

[0010] (3) Use the optimized gradient to update the model parameters Vp(x, z), Vs(x, z), ρ(x, z), and then repeat steps (2) and (3) until the inversion result meets the accuracy requirements and obtain the final inversion model parameter Vp real (x, z), Vs real (x, z),ρ real (x, z);

[0011] (4) Vp real (x, z), Vs real (x, z),ρ real (x, z) is brought in as a fixed model, and full waveform inversion is performed using the adjoint state method. At this time, only the anisotropic parameters ε(x, z) and δ(x, z) are inverted during the inversion process, and the gradient profile is obtained: gε k (x, z) and gδ k (x, z), where k is the number of inversions, k=1, 2, 3, …, K, and K is the total number of inversions;

[0012] (5) Using the Adam optimization algorithm formula (2) in deep learning, gε k (x, z) and gδ k The (x, z) gradient profile is optimized to obtain the optimized gradient Gε k (x, z) and Gδ k (x, z), g in formula (2) k represents the gradient before processing, G k represents the processed gradient, and Respectively represent the hyperparameters that control the average exponential decay rate of the gradient term and the gradient square term, and are generally taken as , , s0=0, r0=0, represents the multiplication of corresponding elements, is a very small constant, generally taken as 1×10 -10 , the formula (2) is: ;

[0013] (6) Use the optimized gradient to update the anisotropic parameters ε(x, z) and δ(x, z), and then repeat steps (4) to (6) until the inversion result meets the accuracy requirement (generally, 100 inversions are enough to meet the accuracy requirement), and obtain the final anisotropic parameter ε real (x, z) and δ real Inversion result of (x, z).

[0014] The beneficial effects of the present invention compared with the prior art are as follows: the method of the present invention inverts the anisotropic parameters ε, δ and the model parameters Vp, Vs, ρ in two steps, thereby reducing the crosstalk effect of simultaneous inversion of multiple parameters, and selecting different gradient optimization algorithms according to the characteristics of the parameters during the inversion process, thereby specifically improving the accuracy of gradient optimization, accelerating the inversion convergence speed, saving time, and significantly improving the effect of full waveform inversion of anisotropic parameters of VTI media. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 This is the flow chart of inversion of anisotropic parameters of conventional VTI media;

[0016] Figure 2 The flowchart of the inversion method of anisotropic parameters of VTI media based on the two-step method;

[0017] Figure 3 is the real model diagram of Marmousi, where (a) is Vp The real model diagram, (b) is Vs The real model diagram, (c) ρ The real model diagram, (d) is ε The real model diagram 、 (e) δ The real model diagram;

[0018] Figure 4 is the initial model diagram, where (a) is Vp The initial model diagram of Vs The initial model diagram, (c) is ρ The initial model diagram, (d) is ε Initial model diagram of 、 (e) δ The initial model diagram of

[0019] Figure 5 is the model parameter gradient profile of the tenth iterative inversion, where (a) is Vp The gradient profile of Vs The gradient profile of ρ Gradient profile of

[0020] Figure 6is the gradient profile of the model parameters after optimization based on the conjugate gradient method, where (a) is Vp The optimized gradient profile, (b) is Vs The optimized gradient profile (c) is ρ Optimized gradient profile;

[0021] Figure 7 is the model parameter profile obtained by the final inversion, where (a) is Vp The final inversion profile is (b) Vs The final inversion profile (c) is ρ The final inverted profile;

[0022] Figure 8 is the anisotropic parameter gradient profile of the tenth iterative inversion, where (a) is ε The cross-sectional view of ρ sectional view of

[0023] Fig. 9 is the gradient profile of anisotropic parameters after optimization based on the Adam algorithm, where (a) is ε The cross-sectional view of ρ sectional view of

[0024] Fig.10 The inversion results of the anisotropic parameters of VTI media obtained by conventional methods, where (a) is ε The inversion result diagram of ρ The inversion result diagram of

[0025] Fig.11 The inversion results of the anisotropic parameters of VTI media obtained based on the two-step method, where (a) is ε The inversion result diagram of ρ The inversion result diagram. DETAILED DESCRIPTION

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

[0027] The method of the present invention is used in the inversion of anisotropic parameters of conventional VTI media (see the specific process for details). Figure 1 Based on the above, the anisotropic parameters ε, δ and the model parameters Vp, Vs, ρ are inverted in two steps, and different gradient optimization algorithms are selected according to the characteristics of the parameters during the inversion process (see the specific process for details). Figure 2). Experiments have shown that this method reduces the crosstalk effect of simultaneous inversion of multiple parameters, improves the accuracy of gradient optimization in a targeted manner, thereby accelerating the inversion convergence speed and significantly improving the effect of full waveform inversion of anisotropic parameters of VTI media.

[0028] Example 1

[0029] 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 4610m, the vertical depth is 1730m, and the grid step size in the x and z directions of the model is 10m.

[0030] Based on this model, a total of 461 shots were fired, with the shot points evenly distributed at intervals of 10 m directly above the model. 461 detection points were set, and the detection points were also evenly distributed at intervals of 10 m directly above the model. The depths of the shot points and the detection points were both 0 m. The shot points were fired in sequence from left to right at the top of the model. Each time they were fired, all the detection points received the shot at the same time, and actual seismic records of a total of 461 shots were obtained.

[0031] The specific implementation mode of the present invention is described in detail below. Figure 2 As shown, the specific steps are as follows:

[0032] (1) The Marmousi real model with a horizontal length of 4610 m and a vertical depth of 1730 m (e.g. Figure 3 The model is meshed with a mesh step size of 10 m in the x and z directions. The initial velocity model Vp(x, z), Vs(x, z), ρ(x, z) and anisotropic parameters ε(x, z) and δ(x, z) are obtained by smoothing (as shown in Figure 4 As shown in the figure, x and z represent the grid position coordinates (x=1, 2, 3, …, 461, z=1, 2, 3, …, 173); based on the actual seismic record Record, the adjoint state method is used to perform full waveform inversion. During the inversion process, only Vp(x, z), Vs(x, z) and ρ(x, z) are inverted, and ε(x, z) and δ(x, z) remain unchanged. In each inversion, the gradient profile is obtained: gVp k (x, z), gVs k (x, z) and gρ k (x, z) (take the model parameter gradient profile of the tenth iteration inversion as Figure 5 As shown), where k is the number of inversions, k=1, 2, 3, …, 100;

[0033] (2) Using the conjugate gradient optimization algorithm formula (1), gVp k (x, z), gVs k(x, z) and gρ k The three gradient profiles (x, z) are optimized to obtain the optimized gradient profile GVp k (x, z), GVs k (x, z) and Gρ k (x,z) (The results after optimization based on the conjugate gradient method are as follows Figure 6 As shown), g in formula (1) k represents the gradient before processing, G k represents the processed gradient, T represents the transposition of the array, and the formula (1) is: ;

[0034] (3) Use the optimized gradient to update the model parameters Vp(x, z), Vs(x, z), ρ(x, z), and then repeat the above steps for 100 times to obtain the final inverted model parameter Vp real (x, z), Vs real (x, z),ρ real (x, z) Figure 7 shown);

[0035] (4) Vp real (x, z), Vs real (x, z),ρ real (x, z) is brought in as a fixed model, and full waveform inversion is performed using the adjoint state method. At this time, only the anisotropic parameters ε(x, z) and δ(x, z) are inverted during the inversion process, and the gradient profile is obtained: gε k (x, z) and gδ k (x, z) (the anisotropic parameter gradient profile of the tenth iteration inversion is as follows Figure 8 As shown), where k is the number of inversions, k=1, 2, 3, …, 100;

[0036] (5) Using the Adam optimization algorithm formula (2) in deep learning, gε k (x, z) and gδ k The (x, z) gradient profile is optimized to obtain the optimized gradient Gε k (x, z) and Gδ k (x, z) (The anisotropic parameter gradient profile after optimization based on the Adam algorithm is as follows Fig. 9 As shown), g in formula (2) k represents the gradient before processing, G k represents the processed gradient, and Respectively represent the hyperparameters that control the average exponential decay rate of the gradient term and the gradient square term, and are generally taken as , , s0=0, r0=0, represents the multiplication of corresponding elements, is a very small constant, generally taken as 1×10 -10 , the formula (2) is: ;

[0037] (6) Use the optimized gradient to update the anisotropic parameters ε(x, z) and δ(x, z), and then repeat the above steps for 100 times to obtain the final anisotropic parameter ε real (x, z) and δ real Inversion result of (x, z).

[0038] In order to demonstrate the beneficial effects of the method of the present invention, here we compare it with the conventional anisotropic parameter inversion (process as shown in Figure 1 The results of anisotropy parameters ε(x, z) and δ(x, z) obtained by θ(x, z) are compared with those obtained by θ(x, z). Fig.10 The inversion results of ε(x, z) and δ(x, z) obtained by conventional anisotropic parameter inversion are: Fig.11 The inversion results of ε(x, z) and δ(x, z) obtained by the two-step method for inversion of anisotropic parameters of VTI media. By comparing the two figures, it can be found that the inversion results obtained by the method of this embodiment have clearer structural description and stronger continuity of the event axis, especially in the medium and deep areas, and its accuracy has been greatly improved, which proves that this method provides a better solution for the inversion of anisotropic parameters of VTI media.

Claims

1. A two-step method for inversion of anisotropic parameters of VTI media, characterized in that: The method specifically comprises the following steps: (1) Given the initial P-wave velocity model parameters Vp(x,z), S-wave velocity model parameters Vs(x,z), density model parameters ρ(x,z), and anisotropy parameters ε(x,z) and δ(x,z), where x and z represent spatial position coordinates, x = 1, 2, 3, ..., Nx, z = 1, 2, 3, ..., Nz, and Nx and Nz represent the total number of model transverse and longitudinal grid points, respectively, the shot points are evenly placed on the top of the model, and the receiver points are also evenly placed on the top of the model. 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 the signal at the same time, and the actual seismic record Record is obtained. Based on the actual seismic record Record, full waveform inversion is performed using the adjoint state method. During the inversion process, only Vp(x,z), Vs(x,z) and ρ(x,z) are inverted, and ε(x,z) and δ(x,z) remain unchanged. In each inversion, the P-wave velocity gradient profile gVp is obtained. k (x,z), shear wave velocity gradient profile gVs k (x,z) and the density gradient profile gρ k (x, z), where k is the number of inversions, k = 1, 2, 3, ..., K, K is the total number of inversions; (2) Using the conjugate gradient optimization algorithm formula (1), gVp k (x,z),gVs k (x,z) and gρ k The three gradient profiles (x, z) are optimized to obtain the optimized gradient GVp k (x,z), GVs k (x,z) and Gρ k (x,z), g in formula (1) k represents the gradient before processing, G k represents the processed gradient, T represents the transposition of the array, and the formula (1) is: (3) Use the optimized gradient to update the model parameters Vp(x,z), Vs(x,z), ρ(x,z), and then repeat steps (2) and (3) until the inversion result meets the accuracy requirements, and obtain the final inverted P-wave velocity model parameter Vp real (x,z), shear wave velocity model parameter Vs real (x,z), density model parameter ρ real (x,z); (4) Vp real (x,z),Vs real (x,z),ρ real (x, z) is brought in as a fixed model, and full waveform inversion is performed using the adjoint state method. At this time, only the anisotropic parameters ε(x, z) and δ(x, z) are inverted during the inversion process to obtain the ε(x, z) gradient profile gε k The gradient profile gδ of (x,z) and δ(x,z) k (x, z), where k is the number of inversions, k = 1, 2, 3, ..., K, K is the total number of inversions; (5) Using the Adam optimization algorithm formula (2) in deep learning, gε k (x,z) and gδ k The (x,z) gradient profile is optimized to obtain the optimized gradient Gε k (x,z) and Gδ k (x,z), g' in formula (2) k represents the gradient before treatment, G' k represents the processed gradient, β1 and β2 represent the hyperparameters controlling the average exponential decay rate of the gradient term and the gradient square term, ⊙ represents the multiplication of the corresponding elements, and e is a very small constant. The formula (2) is: (6) Use the optimized gradient to update the anisotropic parameters ε(x, z) and δ(x, z), and then repeat steps (4)-(6) until the inversion result meets the accuracy requirement, and obtain the final anisotropic parameter ε real (x,z) and δ real The inversion result of (x,z).

Images