Full waveform inversion method, device, medium and equipment based on self-adjoint VTI equation

Abstract

The present invention relates to a full waveform inversion method, device, medium and equipment based on a self-accompanied VTI equation. The method comprises the following steps: splitting a VTI acoustic wave equation and a VTI elastic wave equation to obtain split VTI acoustic wave equation and VTI elastic wave equation; constructing auxiliary variables p and q and introducing them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain a self-accompanied VTI acoustic wave equation and a self-accompanied VTI elastic wave equation; solving the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate a forward wave field u and predicted data d. pre ; Using the objective function and prediction data d pre Calculate the residual d res Solve the self-adjoint VTI acoustic wave equation and the self-adjoint VTI elastic wave equation, and back-propagate the residual d res Calculate the accompanying wave field u * Based on the gradient formula of the objective function to the model parameters, the forward wave field u and the accompanying wave field u are used to calculate the * Calculate the gradient calculation step size α and update the model m through the gradient and iterative update formula to determine whether the objective function converges.

Description

Technical Field

[0001] The present invention relates to a full waveform inversion method, device, medium and equipment based on a self-adjoint VTI equation, and belongs to the technical field of geological exploration. Background Art

[0002] Full-waveform inversion (FWI) of vertically symmetric transversely isotropic (VTI) media can effectively model the subsurface in regions with VTI anisotropy. FWI typically achieves high-precision inversion of the subsurface by iteratively updating a good initial model using a local optimization algorithm. This iterative update requires the gradient of the objective function as the update direction, and the calculation of the gradient requires accurate numerical simulation of both the forward and adjoint wavefields.

[0003] Currently, the VTI elastic wave equation can be used to simulate the forward propagation of the VTI elastic wave field. To simulate the VTI acoustic wave field, the shear wave velocity along the symmetry axis in the elastic wave equation is usually set to 0 (Duveneck et al., 2008), thereby obtaining the VTI acoustic anisotropic wave equation and accurately realizing the forward propagation of the VTI acoustic wave field ( Figure 1 (a)). In the process of realizing full waveform inversion of VTI anisotropic media, the adjoint state method is usually used to obtain the gradient of the objective function with respect to the model parameters. The adjoint state method requires the derivation of the adjoint wave equation and the adjoint wave field. However, the VTI elastic wave equation and the VTI acoustic wave equation that are currently commonly used are not self-symmetric. Therefore, there are obvious problems with the adjoint wave equation obtained directly by transposing the wave equation operator. First, although the elastic wave adjoint wave equation of VTI elastic media can accurately realize the numerical simulation of the adjoint wave field, the free boundary conditions of the adjoint wave equation of the direct transposed wave equation operator are difficult to implement effectively. The main reason is that the physical meaning of the adjoint variables of the adjoint wave equation after transposition is unclear. Therefore, in the current references, most scholars only consider the absorbing boundary in the wave field simulation and do not implement the free boundary (Ren and Liu, 2016). Secondly, there are two problems with the VTI acoustic wave adjoint state equation derived by direct transposition: 1. When the adjoint equation does not consider the influence of the model parameters changing with space, the wave field simulated by the adjoint wave equation will have simulation instability problems at the mutation position of the anisotropic parameters (such as Figure 1 (b)); 2. When the adjoint equation considers the influence of model parameters changing with space, the adjoint wave field will have the problem of weak amplitude value in deep part ( Figure 1 (d)). Summary of the Invention

[0004] To address the above-mentioned technical problems, the present invention provides a full waveform inversion method, apparatus, medium, and equipment based on a self-adjoint VTI equation. The adjoint equation derived by this method does not suffer from the problems of unstable adjoint wavefield simulation and weak deep amplitude values. Secondly, the adjoint wavefield of this method can accurately realize the surface free boundary condition. Finally, this method can use the same set of codes to realize the numerical simulation of the forward wavefield and the adjoint wavefield, thereby reducing the complexity of programming.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A full waveform inversion method based on the self-adjoint VTI equation includes:

[0007] The normal stress τ containing the anisotropy parameter ε in the VTI acoustic wave equation and the VTI elastic wave equation xx Perform splitting to obtain split VTI acoustic wave equation and VTI elastic wave equation;

[0008] Construct auxiliary variables p and q and substitute them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation;

[0009] Solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate the forward wave field u and the predicted data d pre ;

[0010] Using the objective function and prediction data d pre Calculate the residual d res ;

[0011] Solve the self-adjoint VTI acoustic wave equation and the self-adjoint VTI elastic wave equation, and backpropagate the residual d res Calculate the accompanying wave field u * ;

[0012] Based on the gradient formula of the objective function to the model parameters, the forward wave field u and the accompanying wave field u are used to calculate the gradient of the model parameters. * Calculating gradients

[0013] Calculate the step size α and pass the gradient Update the model m with the iterative update formula, determine whether the objective function converges, output the result if the objective function converges, and perform the next iteration if the objective function does not converge until the objective function converges.

[0014] In the full waveform inversion method based on the self-adjoint VTI equation, preferably, the specific formula of the split VTI acoustic wave equation is as follows:

[0015]

[0016] Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz Represents the normal stress in the horizontal and vertical directions respectively; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent the density, P-wave velocity and anisotropy parameter, respectively.

[0017] In the full waveform inversion method based on the self-adjoint VTI equation, preferably, the specific formula of the split VTI elastic wave equation is as follows:

[0018]

[0019] Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz represent the normal stress in the horizontal and vertical directions respectively; τ xz represents the shear stress in the shear direction; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter, respectively.

[0020] The full waveform inversion method based on the self-accompanied VTI equation is preferably formulated as follows:

[0021]

[0022] Among them, p and q represent auxiliary variables; τ xx2 represents the partial normal stress in the horizontal direction after splitting; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent density, P-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

[0023] The full waveform inversion method based on the self-adjoint VTI equation is preferably a self-adjoint VTI elastic wave equation having the following specific formula:

[0024]

[0025] Among them, p and q represent auxiliary variables; τ xx2 represents the partial normal stress in the horizontal direction after splitting; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

[0026] In the full waveform inversion method based on the self-adjoint VTI equation, preferably, the gradient formula of the objective function with respect to the model parameters is as follows:

[0027]

[0028] Where, “T” represents the transpose of the matrix; m represents the model parameters; u represents the forward wave field; A represents the forward operator; u * represents the accompanying wave field; x represents the spatial coordinate of a point underground.

[0029] In the full waveform inversion method based on the self-adjoint VTI equation, preferably, the objective function is formulated as follows:

[0030]

[0031] Among them, m represents the model parameter, d pre Represents the predicted data generated by model m, d obs Represents observation data collected in the field.

[0032] A second aspect of the present invention provides a full waveform inversion device based on a self-adjoint VTI equation, comprising:

[0033] The first processing unit is used to process the normal stress τ containing the anisotropy parameter ε in the VTI acoustic wave equation and the VTI elastic wave equation. xx Perform splitting to obtain split VTI acoustic wave equation and VTI elastic wave equation;

[0034] The second processing unit is used to construct auxiliary variables p and q and bring them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain a self-accompanied VTI acoustic wave equation and a self-accompanied VTI elastic wave equation;

[0035] The third processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate the forward wave field u and the predicted data d pre ;

[0036] The fourth processing unit is used to use the objective function and the prediction data d pre Calculate the residual d res ;

[0037] The fifth processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation, and back-propagate the residual d res Calculate the accompanying wave field u * ;

[0038] The sixth processing unit is used to calculate the gradient formula of the model parameters based on the objective function, using the forward wave field u and the accompanying wave field u * Calculating gradients

[0039] The seventh processing unit is used to calculate the step size α and pass the gradient Update the model m with the iterative update formula, determine whether the objective function converges, output the result if the objective function converges, and perform the next iteration if the objective function does not converge until the objective function converges.

[0040] A third aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of any one of the above-mentioned full waveform inversion methods based on the self-adjoint VTI equations.

[0041] A fourth aspect of the present invention provides a computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the steps of any one of the above-mentioned full waveform inversion methods based on the self-adjoint VTI equation are implemented.

[0042] The present invention has the following advantages due to the adoption of the above technical solution:

[0043] 1. In the full waveform inversion based on the self-adjoint VTI equation, the present invention firstly normalizes the stress τ xx The VTI equation is self-adjointly implemented by splitting the wavefield and then constructing auxiliary variables p and q. This technology can reduce the programming complexity of full waveform inversion of VTI media and solve the problems of unstable adjoint wavefield simulation and weak deep amplitude values ​​in the existing technology. In addition, the forward wavefield and the adjoint wavefield can use the same free boundary implementation scheme.

[0044] 2. The present invention uses the normal stress τ of the VTI acoustic wave equation to xxThe self-adjoint VTI acoustic wave equation is realized by splitting and constructing auxiliary variables p and q. The full waveform inversion based on the self-adjoint VTI acoustic wave equation can effectively avoid the simulation instability of the adjoint wave field and the weak amplitude value in the deep.

[0045] 3. The present invention uses the normal stress τ of the VTI elastic wave equation to calculate the xx The self-adjoint VTI elastic wave equation is realized by splitting and constructing auxiliary variables p and q. The full waveform inversion based on the self-adjoint VTI elastic wave equation can effectively avoid the problem that the adjoint wave field cannot effectively implement the free boundary. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 The forward and companion wave fields corresponding to the Marmousi2 model provided in one embodiment of the present invention at 1700 milliseconds; Figures (a) and (c) are the forward wave fields calculated by the VTI acoustic wave equation (Formula (2)); Figure (b) is the companion wave field calculated by the companion equation without considering the spatial variation of model parameters (Formula (8)); Figure (d) is the companion wave field calculated by the companion equation considering the spatial variation of model parameters (Formula (9)); Figure (e) is the forward wave field calculated by the self-accompanying VTI acoustic wave equation (Formula (15)) proposed by the present invention, and Figure (f) is the companion wave field calculated by the self-accompanying VTI acoustic wave equation (Formula (15)) proposed by the present invention. The earthquake source location is 6.5 km, the depth is 10 m, and the surface is a free boundary;

[0047] Figure 2 are the real model and initial model of Marmousi2; among them, Figures (a), (c), (e), (g) and (i) are the real P-wave velocity model, S-wave velocity model, density model and anisotropic parameter model respectively, and Figures (b), (d), (f), (h) and (j) are the corresponding initial models respectively;

[0048] Figure 3 The VTI acoustic wave full waveform inversion results of the Marmousi2 model, wherein Figure (a) is the P-wave velocity inverted by the prior art; Figure (b) is the P-wave velocity inverted by the present invention;

[0049] Figure 4 Figure 1 is the full waveform inversion result of VTI elastic waves of the Marmousi2 model, where (a) and (c) are the P-wave velocity and S-wave velocity inverted by the prior art, respectively; (b) and (d) are the P-wave velocity and S-wave velocity inverted by the present invention, respectively.

[0050] Figure 5 This is a flow chart of the full waveform inversion equation based on the self-adjoint VTI equation. DETAILED DESCRIPTION

[0051] To make the objectives, technical solutions, and advantages of the present invention more clear, the technical solutions of the present invention are described clearly and completely below. Obviously, the embodiments described are only some of the embodiments of the present invention, not all of them. All other embodiments derived by ordinary persons in this field based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0052] Unless otherwise defined, the technical or scientific terms used in the present invention shall have the usual meanings understood by persons of ordinary skill in the field to which the present invention belongs. The words "first", "second", "third", "fourth" and similar terms used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. Words such as "include" or "comprise" mean that the elements or objects preceding the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Words such as "connect" or "connected" are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect.

[0053] For ease of description, spatially relative terms may be used herein to describe the relationship of one element or feature relative to another element or feature as shown in the figures, such as "inside," "outside," "inner side," "outer side," "lower," "upper," etc. Such spatially relative terms are intended to encompass different orientations of the device in use or operation in addition to the orientation depicted in the figures.

[0054] Existing full waveform inversion usually uses the least squares objective function to fit the surface observation records and simulated prediction records to invert the underground medium parameters. Its objective function can be expressed as:

[0055]

[0056] Among them, m represents the model parameter, d pre Represents the predicted data generated by model m, d obs Represents observation data collected in the field.

[0057] For full waveform inversion of acoustic waves in VTI media, the objective function (1) is constrained by the VTI acoustic wave equation:

[0058]

[0059] Among them, τ xx and τ zz Represents the normal stress in the horizontal and vertical directions respectively; v x and v zRepresent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent density, P-wave velocity and anisotropy parameters respectively; s represents the source, and is loaded in the form of a pure P-wave source. The above formula (2) gives the two-dimensional VTI acoustic wave equation, and there is also a three-dimensional VTI acoustic wave equation. The two are similar in form, and there is a variable τ along the y-axis yy and v y Its specific form will not be described here.

[0060] For full waveform inversion of elastic waves in VTI media, the objective function (1) is constrained by the VTI elastic wave equation:

[0061]

[0062] Among them, ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter, respectively; f is a newly defined variable, which can be accessed through v p and v s Calculated; τ xx and τ zz represent the normal stress in the horizontal and vertical directions respectively; τ xz represents the shear stress in the shear direction; v x and v z Represent the particle vibration velocity in the horizontal and vertical directions respectively; s represents the earthquake source, and is loaded in the form of a pure longitudinal wave source. The same as the acoustic wave equation, the above formula (3) gives the two-dimensional VTI elastic wave equation, and there is also a three-dimensional VTI elastic wave equation. The two are similar in form, and there is a variable τ related to the y-axis yy , τ xy , τ yz and v y Its specific form will not be described here.

[0063] Formulas (2) and (3) can both be expressed in matrix form as:

[0064] Au(x,t)=s(x s ,t) (4)

[0065] Where u represents the forward wave field, which includes multiple components, specifically defined in formulas (2) and (3), x represents the spatial coordinates of a point underground; s represents the earthquake source, x s represents the spatial coordinates of the earthquake source, and t represents the time of seismic wave propagation. A represents the forward operator. For the VTI acoustic wave equation shown in formula (2), its forward operator A can be specifically expressed as:

[0066]

[0067] For the VTI elastic wave equation shown in formula (3), its forward operator A can be specifically expressed as:

[0068]

[0069] The forward wave field in full waveform inversion of VTI media can be numerically simulated using formula (4). The companion wave field requires the corresponding companion equation for numerical simulation. By deriving it through the companion state method, the companion equation corresponding to formula (4) can be expressed in matrix form as follows:

[0070] A T u * (x,t)=d res (x r ,t) (7)

[0071] Among them, u * Represents the accompanying wave field, including and Four components, A T represents the adjoint operator of the adjoint wavefield, which is the transpose of formulas (5) and (6); d res represents the data residual as the accompanying source, x r Represents the spatial coordinates of the detector.

[0072] From formula (7), we can get that without considering the spatial variation of model parameters, the adjoint equation of formula (2) can be expressed as:

[0073]

[0074] in, and Represent the four components of the accompanying wave field; δp represents the residual data of the pressure as the accompanying source. Similarly, considering the model parameters changing with space, the accompanying equation of formula (2) can be expressed as:

[0075]

[0076] The main difference between formula (8) and (9) lies in whether the corresponding spatial derivatives of the model parameters are obtained during the solution of the accompanying equations.

[0077] Through formulas (1), (4) and (7), the gradient of the objective function with respect to the model parameters can be expressed as:

[0078]

[0079] Where, “T” represents the transpose of the matrix; m represents the model parameters; u represents the forward wave field; A represents the forward operator; u *represents the accompanying wave field; x represents the spatial coordinate of a point underground.

[0080] The gradient and local optimization algorithm shown in formula (10) can be used to iteratively update the initial model parameters, thereby achieving the inversion of the subsurface model parameters, that is, to obtain the specific values ​​of the subsurface model parameters, including the P-wave and S-wave velocities and anisotropy parameters. The specific iterative update formula is as follows:

[0081]

[0082] Where m i represents the model parameters of the i-th iteration, m i+1 represents the updated gradient.

[0083] From the forward operators shown in formulas (5) and (6), it can be seen that operator A is not self-symmetric. The forward operator A and its adjoint operator A T are not equal, so different codes are required to realize the numerical simulation of the forward wave field and the adjoint wave field respectively during the implementation of full waveform inversion. At the same time, the adjoint equation derived from the full waveform inversion based on the VTI acoustic wave equation (Formula 2) has two problems: 1. When A does not consider the spatial variation of model parameters, the adjoint equation will simulate instability at the location where the anisotropic parameters have a sudden change, such as Figure 1 (b) ; 2. When A T When considering the spatial variation of model parameters, the accompanying wave field will have the problem of weak amplitude values ​​in the deep part, such as Figure 1 (d) In addition, the adjoint equation derived from the full waveform inversion of the VTI elastic wave equation (Equation 3) is difficult to implement accurately for free boundaries.

[0084] In view of the above shortcomings of the existing full waveform inversion of VTI media, the present invention proposes a full waveform inversion method based on the self-adjoint VTI equation. This method can accurately realize the stable and correct numerical simulation of the forward wavefield and the adjoint wavefield in the full waveform inversion of VTI media under the condition of free surface. Figure 1 (e) and 1(f)). At the same time, due to the self-adjoint nature of the forward operator, i.e., its self-symmetry, the forward operator and the adjoint operator are the same. Therefore, the numerical simulation of the forward wave field and the adjoint wave field in programming only requires one set of code to execute, which can reduce the complexity of the code implementation.

[0085] like Figure 1 As shown, the full waveform inversion method based on the self-adjoint VTI equation provided by the present invention includes the following specific steps:

[0086] In the present invention, in order to realize the self-adjoint of the forward operator A, the normal stress τ containing the anisotropy parameter ε in formulas (2) and (3) is firstly xxThe VTI acoustic wave equation shown in formula (2) can be expressed as follows after splitting:

[0087]

[0088] Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz Represents the normal stress in the horizontal and vertical directions respectively; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent the density, P-wave velocity and anisotropy parameter, respectively.

[0089] The VTI elastic wave equation shown in formula (3) can be expressed as follows after decomposition:

[0090]

[0091] Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz represent the normal stress in the horizontal and vertical directions respectively; τ xz represents the shear stress in the shear direction; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter, respectively.

[0092] Then the forward operator A is made self-adjoint by constructing auxiliary variables p and q:

[0093]

[0094] and

[0095]

[0096] Substituting formulas (14) and (15) into formula (12) yields the self-adjoint VTI acoustic wave equation:

[0097]

[0098] Among them, p and q represent auxiliary variables; τ xx2 represents the partial normal stress in the horizontal direction after splitting; v x and v zRepresent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent density, P-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

[0099] Substituting formulas (14) and (15) into formula (13) yields the self-adjoint VTI elastic wave equation:

[0100]

[0101] Among them, p and q represent auxiliary variables; τ xx2 represents the partial normal stress in the horizontal direction after splitting; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

[0102] The present invention also comprises the steps of:

[0103] The forward wave field u and the predicted data d are calculated by solving the self-accompanying VTI acoustic wave (Formula (16)) or elastic wave (Formula (17)) pre ;

[0104] By using the objective function (Formula (1)) and the predicted data d pre Calculate the residual d res ;

[0105] By backpropagating the residual d from the self-adjoint VTI acoustic wave (Formula (16)) and elastic wave (Formula (17)) res Calculate the accompanying wave field u *

[0106] By using formula (10) to use the forward wave field u and the accompanying wave field u * Calculating gradients

[0107] Calculate the step size α and pass the gradient Update the model m using formula (11) and determine whether the objective function converges. If the objective function converges, output the result. If the objective function does not converge, perform the next iteration until the objective function converges.

[0108] The self-accompanied VTI acoustic and elastic wave equations shown in formulas (16) and (17) have forward operators A equal to their transpose A T Therefore, the full waveform inversion based on formula (16) and formula (17) has the same forward operator A and adjoint operator A T In the implementation of full waveform inversion, the same set of codes can be used to realize the numerical simulation of the forward wavefield and the companion wavefield, and the companion wavefield can use the same free boundary implementation scheme as the forward wavefield.

[0109] The following describes in detail through specific examples how the present invention can effectively achieve full waveform inversion of acoustic and elastic waves in VTI media.

[0110] The example used in this invention is the Marmousi2 model. The real longitudinal wave velocity model is as follows: Figure 2 As shown in (a), the real shear wave velocity model is as follows Figure 2 As shown in (c), the real density model is Figure 2 As shown in (e), the true anisotropic parameter ε model is as follows Figure 2 As shown in (g), the true anisotropy parameter δ is Figure 2 (i) The initial velocity models required for full waveform inversion are as follows: Figure 1 (b), 1(d) and 1(f).

[0111] The present invention only shows the implementation process of full waveform inversion of acoustic waves in VTI media. Figure 2 The initial P-wave velocity model shown in (b) is updated, the density is generated using the Gardner formula, and the anisotropy parameters are not updated during the inversion process using the initial model. The final P-wave velocity inversion result obtained by the existing technology is as follows: Figure 3 (a) shows the inversion results of the longitudinal wave velocity of the present invention. Figure 3 By comparison, it can be seen that the prior art suffers from the problem of weak deep amplitude values ​​in the accompanying wave field, and the deep update effect of the final inversion result is not as good as that of the present invention. The inversion result of the present invention is also closer to the true P-wave velocity model.

[0112] The present invention is to implement the elastic wave full waveform inversion process of VTI medium. Figure 2 (b) and Figure 2 The initial P-wave velocity and S-wave velocity models shown in (d) are updated, the density is generated using the Gardner formula, and the anisotropy parameters are not updated during the inversion process using the initial model. The final P-wave velocity and S-wave velocity inversion results obtained by the existing technology are as follows: Figure 4 (a) and 4 (c), the inversion results of the longitudinal wave velocity and the shear wave velocity obtained by the present invention are as follows Figure 4(b) and 4(d). Figure 4 Comparison between (a) and 4(c) shows that the adjoint wave field of the prior art cannot effectively implement the free boundary, and the final inversion result contains a lot of noise. Figure 4 (b) and 4(d) can accurately implement the free boundary because of their accompanying wave fields. Therefore, the final inversion results are closer to the true velocity model than the existing technology.

[0113] In the full waveform inversion based on the self-adjoint VTI equation, the present invention firstly adjusts the normal stress τ xx The VTI equation is self-adjointly implemented by splitting the wavefield and then constructing auxiliary variables p and q. This technology can reduce the programming complexity of full waveform inversion of VTI media and solve the problems of unstable adjoint wavefield simulation and weak deep amplitude values ​​in the existing technology. In addition, the forward wavefield and the adjoint wavefield can use the same free boundary implementation scheme.

[0114] The method of the present invention has the following beneficial effects: first, the adjoint equation derived by the method does not have the problems of unstable adjoint wavefield simulation and weak deep amplitude value; second, the adjoint wavefield of the method can accurately realize the free boundary; finally, the method can use the same set of codes to realize the numerical simulation of the forward wavefield and the adjoint wavefield, thereby reducing the complexity of programming.

[0115] In the technical solution of the present invention, the objective function of full waveform inversion is an objective function based on least squares fitting. This objective function can also be other publicly available objective functions, such as the optimal transmission objective function, the normalized objective function, etc. However, these other objective functions still need to be combined with the self-adjoint equations of the present invention to achieve stable and efficient inversion. In addition, in order to simplify the derivation, the present invention only describes the derivation of the self-adjoint VTI acoustic wave equation and elastic wave equation in the two-dimensional case. Self-adjointness can also be achieved in the three-dimensional case through the above process.

[0116] A second aspect of the present invention provides a full waveform inversion device based on a self-adjoint VTI equation, comprising:

[0117] The first processing unit is used to process the normal stress τ containing the anisotropy parameter ε in the VTI acoustic wave equation and the VTI elastic wave equation. xx Perform splitting to obtain split VTI acoustic wave equation and VTI elastic wave equation;

[0118] The second processing unit is used to construct auxiliary variables p and q and bring them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain a self-accompanied VTI acoustic wave equation and a self-accompanied VTI elastic wave equation;

[0119] The third processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate the forward wave field u and the predicted data d pre ;

[0120] The fourth processing unit is used to use the objective function and the prediction data d pre Calculate the residual d res ;

[0121] The fifth processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation, and back-propagate the residual d res Calculate the accompanying wave field u * ;

[0122] The sixth processing unit is used to calculate the gradient formula of the model parameters based on the objective function, using the forward wave field u and the accompanying wave field u * Calculating gradients

[0123] The seventh processing unit is used to calculate the step size α and pass the gradient Update the model m with the iterative update formula, determine whether the objective function converges, output the result if the objective function converges, and perform the next iteration if the objective function does not converge until the objective function converges.

[0124] A third aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of any one of the above-mentioned full waveform inversion methods based on the self-adjoint VTI equations.

[0125] A fourth aspect of the present invention provides a computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the steps of any one of the above-mentioned full waveform inversion methods based on the self-adjoint VTI equation are implemented.

[0126] The present invention is described in terms of flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to specific embodiments. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as a combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0127] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0128] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0129] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A full waveform inversion method based on the self-adjoint VTI equation, characterized in that: include: The normal stress τ in the VTI acoustic wave equation and the VTI elastic wave equation xx Split into τ xx1 and τ xx2 , and obtain the split VTI acoustic wave equation and VTI elastic wave equation; Construct auxiliary variables p and q and substitute them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation; Solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate the forward wave field u and the predicted data d pre ; Using the objective function and prediction data d pre Calculate the residual d res ; Solve the self-adjoint VTI acoustic wave equation and the self-adjoint VTI elastic wave equation, and backpropagate the residual d res Calculate the accompanying wave field u * ; Based on the gradient formula of the objective function to the model parameters, the forward wave field u and the accompanying wave field u are used to calculate the gradient of the model parameters. * Calculating gradients Calculate the step size α and pass the gradient Update the model m with the iterative update formula, determine whether the objective function converges, output the result if the objective function converges, and perform the next iteration if the objective function does not converge until the objective function converges.

2. The full waveform inversion method based on the self-adjoint VTI equation according to claim 1 is characterized in that: The specific formula of the split VTI acoustic wave equation is as follows: Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz Represents the normal stress in the horizontal and vertical directions respectively; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent the density, P-wave velocity and anisotropy parameter, respectively.

3. The full waveform inversion method based on the self-adjoint VTI equation according to claim 2 is characterized in that: The specific formula of the split VTI elastic wave equation is as follows: Among them, τ xx =τ xx1 +τ xx2 , τ xx and τ zz represent the normal stress in the horizontal and vertical directions respectively; τ xz represents the shear stress in the shear direction; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter, respectively.

4. The full waveform inversion method based on the self-adjoint VTI equation according to claim 3 is characterized in that: The specific formula of the self-adjoint VTI acoustic wave equation is as follows: Among them, p and q represent auxiliary variables; τ xx2 represents the normal stress in the horizontal direction after splitting; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , ε and δ represent density, P-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

5. The full waveform inversion method based on the self-adjoint VTI equation according to claim 4 is characterized in that: The specific formula of the self-adjoint VTI elastic wave equation is as follows: Among them, p and q represent auxiliary variables; τ xx2 represents the normal stress in the horizontal direction after splitting; v x and v z Represent the particle vibration speed in the horizontal and vertical directions respectively; ρ, v p , v s , ε and δ represent density, P-wave velocity, S-wave velocity and anisotropy parameter respectively; s represents earthquake source; and They represent the horizontal spatial first-order derivative operator, the vertical spatial first-order derivative operator, and the time first-order derivative operator respectively.

6. The full waveform inversion method based on the self-adjoint VTI equation according to claim 5 is characterized in that: The gradient formula of the objective function with respect to the model parameters is as follows: Where "T" represents the transpose of the matrix; m represents the model parameters; u represents the forward wave field, A represents the forward operator; u * represents the accompanying wave field; x represents the spatial coordinate of a point underground.

7. The full waveform inversion method based on the self-adjoint VTI equation according to claim 6 is characterized in that: The objective function formula is as follows: Among them, m represents the model parameter, d pre Represents the predicted data generated by model m, d obs Represents observation data collected in the field.

8. A full waveform inversion device based on the self-adjoint VTI equation, characterized in that: include: The first processing unit is used to process the normal stress τ in the VTI acoustic wave equation and the VTI elastic wave equation. xx Split into τ xx1 and τ xx2 , and obtain the split VTI acoustic wave equation and VTI elastic wave equation; The second processing unit is used to construct auxiliary variables p and q and substitute them into the split VTI acoustic wave equation and VTI elastic wave equation to obtain a self-accompanied VTI acoustic wave equation and a self-accompanied VTI elastic wave equation; The third processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation to calculate the forward wave field u and the predicted data d pre ; The fourth processing unit is used to use the objective function and the prediction data d pre Calculate the residual d res ; The fifth processing unit is used to solve the self-accompanied VTI acoustic wave equation and the self-accompanied VTI elastic wave equation, and back-propagate the residual d res Calculate the accompanying wave field u * ; The sixth processing unit is used to calculate the gradient formula of the model parameters based on the objective function, using the forward wave field u and the accompanying wave field u * Calculating gradients The seventh processing unit is used to calculate the step size α and pass the gradient Update the model m with the iterative update formula, determine whether the objective function converges, output the result if the objective function converges, and perform the next iteration if the objective function does not converge until the objective function converges.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the full waveform inversion method based on the self-adjoint VTI equation according to any one of claims 1 to 7 are implemented.

10. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the steps of the full waveform inversion method based on the self-adjoint VTI equation according to any one of claims 1 to 7 are implemented.

Images