Three-dimensional anisotropic medium time-space domain earthquake forward modeling method and device, medium and product

Abstract

The invention discloses a three-dimensional anisotropic medium time-space domain seismic forward modeling method and device, a medium and a product, and relates to the technical field of geophysical exploration forward modeling simulation, and the method comprises the steps: presetting an improved time high-order difference operator based on a preset improved space high-order difference operator; respectively discretizing a time derivative and a space derivative in the three-dimensional VTI dielectric acoustic wave equation; based on the three-dimensional VTI dielectric acoustic wave equation, simplifying the time discrete function and the space discrete function to obtain a simplified time-space domain frequency dispersion relational expression; solving the simplified time-space domain dispersion relational expression based on a target optimization algorithm to obtain an optimized time item high-order difference coefficient and an optimized space item high-order difference coefficient; and based on the optimized time item high-order difference coefficient, the optimized space item high-order difference coefficient, the time discrete function and the space discrete function, performing extrapolation solution on the three-dimensional VTI dielectric acoustic wave equation. According to the invention, the numerical precision in anisotropic wave field simulation can be effectively improved.

Description

Technical Field

[0001] The present application relates to the technical field of geophysical prospecting forward modeling, and in particular to a three-dimensional anisotropic medium time-space domain seismic forward modeling method, equipment, medium and product. Background Art

[0002] With the growing demand for oil and gas resources, the exploration and development of new oil and gas fields has become increasingly important. At the same time, in the exploration and development of old oil fields, higher requirements are placed on the accuracy and efficiency of geophysical exploration technology. Traditional seismic exploration technology is mainly based on the assumption that the medium is isotropic, so there is a lot of room for improvement in accuracy. With the increasing complexity of exploration targets and the continuous iteration and updating of exploration technology, it is of great significance to develop new theories, methods and technologies for seismic exploration based on anisotropy theory or even more complex medium characteristics. Seismic data processing plays a connecting role in seismic exploration. Among them, seismic migration imaging and waveform inversion can effectively help understand underground structures or structural characteristics. Imaging and inversion methods and technologies rely on high-precision wave field extension strategies. Therefore, for imaging and inversion technologies based on anisotropy theory, it is necessary to develop high-precision seismic wave field numerical simulation methods and technologies.

[0003] Up to now, scholars have developed various types of numerical algorithms for numerical solutions of wave equations, such as pseudospectral method, finite element method and finite difference method. Among them, the finite difference method is widely favored because of its simplicity, high efficiency and easy implementation. For the wave equation of anisotropic media, the high-order finite difference format (such as regular grid, staggered grid or rotated staggered grid) can be directly used to discretize its spatial partial derivatives, which can obtain higher simulation accuracy in the spatial domain. In terms of time derivative discretization, the second-order difference format is currently mainly used to solve it. The low-precision discrete format is prone to strong time numerical dispersion and even numerical instability. Since the development of high-order finite difference methods, spatial high-order difference discretization has been widely used. However, the existing time high-order difference format is mainly designed for isotropic medium wave equations. Due to its simpler form and fewer parameters involved, this type of numerical solution cannot be directly extended to anisotropic media. The Lax-Wendroff time discretization scheme can be extended to anisotropic media, but its computational efficiency is low. Other high-precision spectral methods involve Fourier transform, so the computational efficiency is difficult to guarantee, especially for three-dimensional wave equations. By introducing paraxial grid nodes into axial nodes, the improved difference template constructed can obtain high-order accuracy in both time and space domains in numerical simulations of isotropic wave equations.

[0004] In general, in the high-precision numerical simulation of anisotropic wave equations, how to construct a reasonable and effective improved differential template to solve it numerically without increasing additional computing consumption has become an urgent problem to be solved. Summary of the invention

[0005] The purpose of this application is to provide a three-dimensional anisotropic medium time-space domain seismic forward modeling method, equipment, medium and product, which can effectively improve the numerical accuracy in anisotropic wave field simulation.

[0006] To achieve the above objectives, this application provides the following solutions:

[0007] In a first aspect, the present application provides a three-dimensional anisotropic medium time-space domain seismic forward modeling method, comprising:

[0008] Obtain the three-dimensional VTI medium acoustic wave equation;

[0009] Based on a preset improved time high-order difference operator, the time derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a time discrete function; based on a preset improved space high-order difference operator, the space derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a space discrete function; the time discrete function contains a time term high-order difference coefficient, and the space discrete function contains a space term high-order difference coefficient;

[0010] Based on the three-dimensional VTI medium acoustic wave equation, simplifying the time discrete function and the space discrete function to obtain a simplified time-space domain dispersion relation;

[0011] Based on the target optimization algorithm, the simplified spatiotemporal dispersion relation is solved to obtain the optimized time term high-order differential coefficient and the optimized space term high-order differential coefficient;

[0012] Based on the optimized high-order differential coefficients of the time term, the optimized high-order differential coefficients of the space term, the time discrete function and the space discrete function, the three-dimensional VTI medium acoustic wave equation is extrapolated and solved;

[0013] Among them, the preset improved time high-order difference operator is:

[0014]

[0015] The preset improved spatial high-order difference operator is:

[0016]

[0017] in, is the time partial derivative of the pressure wave, is the pressure wave at (0,0,0) at time 1, is the pressure wave at the position (0,0,0) at time 0, Δt is the time sampling interval, v pz is the vertical P-wave velocity, h is the spatial step length, N is half the length of the time difference operator, c m,n,0 and c m,n,l is the high-order differential coefficient of the time term, ε and δ are Thomsen anisotropy parameters, e1, e2, e3, f1, f2, f3 are all mathematical operators, is the time partial derivative of the auxiliary wave, is the auxiliary wave at the position (0,0,0) at time 1, is the auxiliary wave at the position (0,0,0) at time 0;

[0018] The pressure waves at time 0 are (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l).

[0019] They are the auxiliary waves at (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,-l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l). is the x-axis spatial partial derivative of the pressure wave, is the y-axis spatial partial derivative of the pressure wave, is the z-axis spatial partial derivative of the auxiliary wave, These are the pressure waves at (m,0,0), (-m,0,0), (0,m,0), (0,-m,0), and (0,0,-m) at time 0, respectively. is the auxiliary wave at the position (0,0,m) at time 0, m, n and l are all serial numbers, M is half the length of the spatial term difference operator, c 0,0,0 and cm,0,0 is the high-order difference coefficient of the spatial term.

[0020] In a second aspect, the present application provides a computer device, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement any one of the three-dimensional anisotropic medium time-space domain seismic forward modeling methods described above.

[0021] In a third aspect, the present application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements any one of the above-mentioned three-dimensional anisotropic medium time-space domain seismic forward modeling methods.

[0022] In a fourth aspect, the present application provides a computer program product, including a computer program, which, when executed by a processor, implements any one of the above-mentioned three-dimensional anisotropic medium time-space domain seismic forward modeling methods.

[0023] According to the specific embodiments provided in this application, this application discloses the following technical effects:

[0024] The present application provides a three-dimensional anisotropic medium time-space domain seismic forward modeling method, equipment, medium and product, and provides a preset improved time high-order difference operator and a preset improved space high-order difference operator to discretize the time derivative and space derivative in the three-dimensional VTI medium acoustic wave equation respectively, thereby adding additional non-axial grid nodes on the basis of the traditional differential template, which can effectively improve the time and space simulation accuracy; then simplifying the process, solving it through the target optimization algorithm, and obtaining optimized high-order difference coefficients for subsequent extrapolation solution of the three-dimensional VTI medium acoustic wave equation, thereby ensuring high-precision extrapolation of the three-dimensional VTI acoustic wave field. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0026] Figure 1 This is an application environment diagram of a three-dimensional anisotropic medium time-space domain seismic forward modeling method in one embodiment of the present application;

[0027] Figure 2 A schematic flow chart of a method for forward modeling a three-dimensional anisotropic medium in time and space domain seismic data according to an embodiment of the present application;

[0028] FIG3 is a schematic diagram of phase velocity surfaces along different directions calculated by different difference schemes in Model I provided in an embodiment of the present application; wherein, Figure 3(a)-Figure 3(c) This is a schematic diagram of model I calculated using the traditional difference method C-Method (12,2). Figure 3(d)-Figure 3(f) This is a schematic diagram of model I calculated according to the method N-Method (12,4) of this application. Figure 3(g)-Figure 3(i) This is a schematic diagram of model I calculated according to the method N-Method (12,6) of this application;

[0029] FIG4 is a schematic diagram of phase velocity surfaces along different directions calculated by different difference schemes in Model II provided in an embodiment of the present application; wherein, Figure 4(a)-Figure 4(c) This is a schematic diagram of Model II calculated using the traditional difference method C-Method (12,2). Figure 4(d)-Figure 4(f) This is a schematic diagram of Model II calculated according to the method N-Method (12,4) of this application. Figure 4(g)-Figure 4(i) This is a schematic diagram of Model II calculated according to the method N-Method (12,6) of this application;

[0030] FIG5 is a wavefield snapshot calculated by the traditional differential method and the improved method in the three-dimensional uniform model at 0.25 s provided in an embodiment of the present application; wherein FIG5(a) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.25 s according to the traditional differential method C-Method (12, 2) and Δt = 1.25 s, FIG5(b) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.25 s according to the traditional differential method C-Method (12, 2) and Δt = 0.625 s, and FIG5(c) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.25 s according to the method N-Method (12, 4) of the present application and Δt = 1.25 s;

[0031] FIG6 is a wavefield snapshot calculated by the traditional differential method and the improved method in the three-dimensional uniform model at 0.5 s provided in an embodiment of the present application; wherein FIG6(a) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.5 s according to the traditional differential method C-Method (12, 2) and Δt = 1.25 s, FIG6(b) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.5 s according to the traditional differential method C-Method (12, 2) and Δt = 0.625 s, and FIG6(c) is a wavefield snapshot calculated by the three-dimensional uniform model at 0.5 s according to the method N-Method (12, 4) of the present application and Δt = 1.25 s;

[0032] FIG. 7 is a modified three-dimensional Marmousi VTI anisotropic medium model provided in an embodiment of the present application; wherein FIG. 7( a ) is a velocity model, and FIG. 7( b ) and FIG. 7( c ) are anisotropic parameter models;

[0033] FIG8 is a snapshot of the instantaneous wave field calculated by the traditional differential method and the method of the present application for a complex model using different time steps provided in an embodiment of the present application; wherein FIG8(a) is a snapshot of the instantaneous wave field calculated by the traditional differential method C-Method (12, 2) and Δt = 0.75 s for the complex model, and FIG8(b) is a snapshot of the instantaneous wave field calculated by the method of the present application N-Method (12, 4) and Δt = 1.0 s for the complex model. DETAILED DESCRIPTION

[0034] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0035] Anisotropy is widely present in the earth's medium, and the development of seismic exploration methods and technologies based on anisotropy theory is crucial to the development of high-precision geophysical technology. Seismic imaging and inversion technology play an important role in the processing of seismic data. Both of these technical methods rely on accurate and efficient numerical simulation of seismic wave fields. However, the existing time-space high-precision numerical algorithms are mainly developed for isotropic media and are difficult to be directly extended to anisotropic media. In view of this, the present application develops a time high-order-space high-order time-space domain finite difference numerical solution based on an improved differential template for the three-dimensional VTI (transverse isotropy with a vertical axis of symmetry, i.e. VTI, vertical symmetry axis transverse isotropy) medium wave equation, which can effectively improve the numerical accuracy in anisotropic wave field simulation.

[0036] In order to make the above-mentioned objects, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0037] The three-dimensional anisotropic medium time-space domain seismic forward modeling method provided in the embodiment of the present application can be applied to Figure 1In the application environment shown. Among them, the terminal 102 communicates with the server 104 through the network. The data storage system can store the data that the server 104 needs to process. The data storage system can be set up separately, or it can be integrated on the server 104, or it can be placed on the cloud or other servers. The terminal 102 can send the three-dimensional uniform VTI medium model to be processed to the server 104. After the server 104 receives the three-dimensional uniform VTI medium model to be processed, for the three-dimensional uniform VTI medium model to be processed, the server 104 performs an improved high-order discrete format derivation for the three-dimensional uniform VTI medium model, simplifies the derivation of the dispersion relationship and the calculation of the differential coefficients, and then realizes a high-precision numerical solution of the three-dimensional VTI acoustic wave equation. The server 104 can feedback the obtained solution results to the terminal 102. In addition, in some embodiments, the three-dimensional anisotropic medium time-space domain seismic forward modeling method can also be implemented independently by the server 104 or the terminal 102. For example, the terminal 102 can directly perform three-dimensional anisotropic medium time-space domain seismic forward modeling on the three-dimensional uniform VTI medium model, or the server 104 can obtain the three-dimensional uniform VTI medium model from the data storage system and perform three-dimensional anisotropic medium time-space domain seismic forward modeling.

[0038] In an exemplary embodiment, Figure 2 As shown, a three-dimensional anisotropic medium time-space domain seismic forward modeling method is provided. The method is executed by a computer device, and can be executed by a computer device such as a terminal or a server alone, or by a terminal and a server together. In the embodiment of the present application, the method is applied to Figure 1 The server 104 in the example is used for explanation, and the steps include the following steps 201 to 205.

[0039] in:

[0040] Step 201, obtain the three-dimensional VTI medium acoustic wave equation. Specifically, VTI medium is a commonly used earth anisotropic medium model in the field of seismic exploration. In the three-dimensional VTI medium, the second-order form wave equation characterized by the pressure component, that is, the three-dimensional VTI medium acoustic wave equation can be expressed as:

[0041]

[0042] Among them, P(x, y, z, t) is the pressure wave field and Q(x, y, z, t) is the auxiliary wave field.

[0043] When the finite difference method is used to numerically solve the three-dimensional VTI medium acoustic wave equation, the traditional scheme usually uses a high-order difference format to calculate the spatial derivative. Solve and apply the second-order difference format to the time derivative However, relevant numerical experiments show that when the time sampling interval is large, the traditional time second-order difference format is prone to cause strong time dispersion and even propagation instability. In addition, the existing time high-order or other advanced difference methods are mainly designed for simple isotropic wave equations and cannot be directly extended to anisotropic wave field simulations. In view of this, the present application provides a time-space high-order finite difference numerical solution in the following step 202 for the three-dimensional anisotropic second-order wave equation.

[0044] Step 202: based on a preset improved time high-order difference operator, the time derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a time discrete function; based on a preset improved space high-order difference operator, the space derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a space discrete function; the time discrete function contains high-order difference coefficients of the time term, and the space discrete function contains high-order difference coefficients of the space term.

[0045] Specifically, relevant research results show that adding additional non-axial grid nodes on the basis of traditional differential templates can effectively improve the temporal and spatial simulation accuracy. Based on this idea, this application constructs a preset improved temporal high-order differential operator and a preset improved spatial high-order differential operator.

[0046] Among them, the preset improved time high-order difference operator is:

[0047]

[0048]

[0049] The preset improved spatial high-order difference operator is:

[0050]

[0051] in, is the time partial derivative of the pressure wave, is the pressure wave at (0,0,0) at time 1, is the pressure wave at the position (0,0,0) at time 0, Δt is the time sampling interval, v pz is the vertical P-wave velocity, h is the spatial step length, N is half the length of the time difference operator, c m,n,0 and c m,n,l is the high-order differential coefficient of the time term, ε and δ are Thomsen anisotropy parameters, e1, e2, e3, f1, f2, f3 are all mathematical operators without clear physical meanings. is the time partial derivative of the auxiliary wave, is the auxiliary wave at the position (0,0,0) at time 1, is the auxiliary wave at the position (0,0,0) at time 0;

[0052] The pressure waves at time 0 are (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l).

[0053] They are the auxiliary waves at (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,-l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l). is the x-axis spatial partial derivative of the pressure wave, is the y-axis spatial partial derivative of the pressure wave, is the z-axis spatial partial derivative of the auxiliary wave, These are the pressure waves at (m,0,0), (-m,0,0), (0,m,0), (0,-m,0), and (0,0,-m) at time 0, respectively. is the auxiliary wave at the position (0,0,m) at time 0, m, n and l are all serial numbers, M is half the length of the spatial term difference operator, c 0,0,0 and c m,0,0 is the high-order difference coefficient of the spatial term.

[0054] In a practical application, the present application also expresses the wave field components in the three-dimensional VTI medium acoustic wave equation based on the characteristics of the plane wave itself as:

[0055]

[0056] in, k x = kcosθcosφ, k y= kcosθsinφ, k z =ksinθ, θ is the plane wave propagation angle, φ is the plane wave azimuth, where the plane wave is obtained by decomposing the pressure wave or the auxiliary wave; that is, the pressure wave field P and the auxiliary wave field Q mentioned above can both be decomposed into plane waves, or can be expressed by plane waves; k x ,k y ,k z are the spatial wave numbers in different directions, and k is the total wave number.

[0057] Step 203, based on the three-dimensional VTI medium acoustic wave equation, simplify the time discrete function and the space discrete function to obtain a simplified time-space domain dispersion relation. Specifically, substitute the equations corresponding to the preset improved time high-order difference operator and the preset improved space high-order difference operator into the three-dimensional VTI medium acoustic wave equation, and after arranging and simplifying, obtain the original time-space domain dispersion relation of the mutual coupling of anisotropic parameters and spatial wave numbers, as shown below:

[0058] G 2 -[(1+2ε)(A+B)+C]r 2 G+(2ε-2δ)(A+B)Cr 4 =0.

[0059] Where G = 2cos(ωΔt)-2, r = v pz Δt / h.

[0060]

[0061] It should be noted that equation G 2 -[(1+2ε)(A+B)+C]r 2 G+(2ε-2δ)(A+B)Cr 4 =0 can be regarded as the variable G is a quadratic form, according to the root-finding formula:

[0062] Δ=[(1+2ε) 2 (A+B) 2 +C 2 +(2-4ε+8δ)(A+B)C]r 4 .

[0063] The anisotropy parameter and the wave number term in the above equation are coupled to each other and are difficult to solve directly. Therefore, in order to facilitate subsequent calculations, this application simplifies it along the fixed propagation direction, setting A=B=C, and the above equation can be simplified to Δ≈(9+8ε+16ε 2 +16δ)A 2 r 4 , substitute this simplified formula into equation G 2-[(1+2ε)(A+B)+C]r 2 G+(2ε-2δ)(A+B)Cr 4 =0, we can get:

[0064]

[0065] Further combined with G = 2cos(ωΔt)-2, the simplified time-space domain dispersion relationship can be obtained:

[0066] 2cos(ωΔt)-2=R 2 (A+B+C),

[0067] Where ω is the angular frequency.

[0068] Step 204, based on the target optimization algorithm, solve the simplified time-space domain dispersion relation to obtain the optimized time term high-order difference coefficient and the optimized space term high-order difference coefficient.

[0069] In another exemplary embodiment of the present application, in order to achieve a more efficient and accurate solution, the specific implementation steps of step 204 are as follows:

[0070] (1) The trigonometric functions at both ends of the simplified spatiotemporal dispersion relation are expanded using Taylor series. In this step, the simplified spatiotemporal dispersion relation is firstly organized into the following equation:

[0071] R -2 [cos(ωΔt)-1]≈F.

[0072]

[0073] Among them, a 0,0,0 、a m,0,0 、a m,n,0 and a m,n,l are all differential coefficients, and each differential coefficient is c 0,0,0 、c m,0,0 、c m,n,0 and c m,n,l The linear combination of

[0074] The above equation can be further organized into the following function formula:

[0075]

[0076] Among them, β=kh, and the range of β is from 0 to π.

[0077] In order to obtain the differential coefficients in the above formula, the differential coefficients in the expanded function formula are organized into the following preset function equations:

[0078]

[0079] Where M1 is the number of differential coefficients, g m is the difference coefficient a 0,0,0 、a m,0,0 、a m,n,0 and a m,n,l In different parameter ranges, χ m The specific form of (β,θ,φ) is as follows:

[0080] When m≤M,

[0081] When M<m≤M+int(N 2 / 4), when m1=m2,

[0082]

[0083] When M<m≤M+int(N 2 / 4), when m1≠m2,

[0084]

[0085] When M+int(N 2 / 4)<m≤M1,m3=m4=m5,

[0086]

[0087] When M+int(N 2 / 4)<m≤M1,m3=m4≠m5,m4=m5≠m3,m3=m5≠m4,

[0088]

[0089] When M+int(N 2 / 4)<m≤M1,m3≠m4≠m5,

[0090]

[0091] Among them, m1, m2, m3, m4, and m5 are all mathematical operators and have no clear physical meanings.

[0092] Considering that the Taylor expansion differential coefficients can only maintain high accuracy in the low-frequency / low-wavenumber region, in order to obtain higher simulation accuracy in the entire wavenumber / frequency range, the target optimization algorithm is used here to solve it.

[0093] (2) With the goal of minimizing the error between the two ends of the preset function equation, an optimization objective function is constructed; the optimization objective function is:

[0094]

[0095] Wherein, E is the error value. It should be noted that in practical applications, the error value is artificially given according to the required data accuracy, such as E is given as 0.0001 or 0.00001, etc.; b is the upper limit value of β.

[0096] (3) Solving the optimization objective function to obtain the high-order differential coefficients of the optimized time term and the high-order differential coefficients of the optimized space term; specifically, using the least squares optimization method to solve the optimization objective function to obtain the high-order differential coefficients of the optimized time term and the high-order differential coefficients of the optimized space term.

[0097] In a specific application example, the least squares algorithm can be used to solve the following equation to obtain the optimized differential coefficients (including the optimized high-order differential coefficients of the time term and the optimized high-order differential coefficients of the space term):

[0098]

[0099] By jointly solving the preset function equation, the optimization objective function and the above equation, it is possible to obtain the optimized differential coefficients for the time-space domain dispersion relationship of the three-dimensional VTI acoustic wave equation using the least squares method.

[0100] Step 205, based on the optimized high-order differential coefficients of the time term, the optimized high-order differential coefficients of the space term, the time discrete function and the space discrete function, the three-dimensional VTI medium acoustic wave equation is extrapolated and solved.

[0101] In addition, in order to verify the feasibility of this application and the accuracy of the data, this application also uses the phase velocity relative error analysis method to test the numerical accuracy of different methods. The traditional method is abbreviated as C-Method (2M, 2), and the improved method in this application is abbreviated as N-Method (2M, 2N). Based on the plane simple harmonic wave theory and equation (i.e., the simplified time-space domain dispersion relationship), the error function is defined as follows:

[0102]

[0103] Among them, κ represents the relative error of the phase velocity. When it is closer to 1, the dispersion error is smaller, and vice versa.

[0104] Two sets of three-dimensional uniform VTI medium models (including model I and model II) are used to calculate and analyze the phase velocity accuracy of the traditional differential method and the differential method in this application. The anisotropy parameters of model I are ε = 0.24, δ = 0.06, and the anisotropy parameters of model II are ε = 0.12, δ = -0.10, and the remaining parameters are v pz=3000m / s, Δt=1.2ms, h=10m. In different difference methods, the length of the spatial domain difference operator is 8, and the length of the time domain difference operator is 2, 4 and 6 respectively. Figures 3 and 4 show the phase velocity surfaces along different directions calculated by different difference schemes in two sets of anisotropic models. Figure 3(a)-Figure 3(i) , Figure 4(a)-Figure 4(i) It can be seen that under the same parameter conditions, the traditional difference method has strong phase velocity dispersion errors along different propagation directions. In contrast, the optimized high-order difference method in this application can significantly suppress the dispersion error and obtain higher numerical accuracy. In addition, in the optimized difference scheme, as the time operator N increases, the numerical accuracy can be further improved.

[0105] This application also uses a three-dimensional uniform VTI medium model to compare the simulation accuracy of the traditional differential method and the method of this application. The dimensions of the three-dimensional model are 200×200×200, and the size of the cubic discrete grid unit is 10m×10m×10m. A Ricker wavelet with a main frequency of 40Hz is selected as the source function and placed in the center of the model to generate vibration. The background velocity value is 3000m / s, and the elliptical anisotropy parameters are ε=0.15 and δ=0.15. The spatial domain differential operator length of the two methods is 12, and the time domain differential operator length is different. The remaining parameters are given in the legend. Figures 5 and 6 show snapshots of the wave fields calculated by the traditional differential method and the improved method in the three-dimensional uniform model at different times. From Figure 5(a)-Figure 5(c) , Figure 6(a)-Figure 6(c) It can be observed that: when the time step is large, the calculation results of the traditional differential method will produce strong time numerical dispersion interference at the wavefront position, as indicated by the arrow. When the time step is reduced, the numerical dispersion is suppressed. In contrast, when the time step is large, the method of the present application can produce high-precision numerical simulation results that are not affected by time dispersion. In particular, when the recording duration is consistent, the traditional differential method must use a smaller time step to ensure a higher time simulation accuracy, thereby increasing the number of iterations and slightly increasing the amount of calculation. In contrast, the method of the present application can use a larger time step to obtain high-precision simulation results, effectively balancing calculation accuracy and efficiency.

[0106] This application also uses a complex model to test the accuracy of different methods. Figure 7 shows a modified three-dimensional Marmousi VTI anisotropic medium model, including a velocity model and an anisotropic parameter model. The size of the three-dimensional model is 2500m×2500m×1800m. The source uses Ricker wavelets and is placed in the shallow surface of the model to generate vibrations. The length of the spatial domain differential operator is 12. A mixed absorbing boundary condition with a thickness of 10 layers is used to suppress artificial truncated boundary reflections. Figure 8 shows snapshots of the instantaneous wave field calculated by the traditional differential method and the method of this application using different time steps for a complex model (i.e., the modified three-dimensional Marmousi VTI anisotropic medium model mentioned above). It should be noted that when the time step is 0.8ms, the traditional differential method will experience propagation instability. Therefore, in this example, the time step of the traditional differential method is 0.75ms, and the time step of the method of this application is 1.0ms. Observation Figure 8(a)-Figure 8(b) It can be seen that in order to ensure the stability of wave field extrapolation, the traditional differential method needs to adopt a smaller time step, which inevitably increases the calculation time. In contrast, the method of this application takes into account the propagation characteristics of waves in both the time domain and the space domain, and can effectively improve the accuracy of numerical simulation while ensuring the calculation efficiency. Therefore, on the whole, the improved time high-order-space high-order finite difference numerical solution in this application can provide an effective wave field continuation tool for subsequent high-precision seismic processing of anisotropic media.

[0107] In summary, in order to develop an accurate and efficient numerical simulation algorithm for seismic wave fields in VTI medium offset imaging, waveform inversion and seismic data processing, this application proposes a time high-order-space high-order time-space domain finite difference numerical solution method for the three-dimensional VTI medium acoustic wave equation. First, the time derivative and space derivative in the equation are discretely solved by using an improved high-order difference format; the improved discrete format is sorted and simplified by combining plane wave theory and mathematical simplification, and a simplified anisotropic time-space domain dispersion relationship is derived; the simplified dispersion relationship is solved by the least squares algorithm to obtain the high-order difference coefficients in the improved difference template; finally, based on the optimization of high-order difference coefficients and the improved discrete format, the high-precision extrapolation of the three-dimensional VTI acoustic wave field can be achieved. Through numerical analysis and model examples, it can be seen that the improved high-order difference scheme in this application can significantly improve the time accuracy and space accuracy of anisotropic wave field simulation. The technical solution in the present invention can provide an effective numerical calculation tool for the development of high-precision seismic exploration technology. The accurate and efficient time-space domain high-order time-space high-order finite difference numerical solution developed for three-dimensional VTI media provided by the present invention can provide an effective theoretical basis for subsequent VTI medium wave field simulation, migration imaging and waveform inversion.

[0108] In an exemplary embodiment, a computer device is provided, which may be a server or a terminal. The computer device includes a processor, a memory, an input / output interface (I / O for short) and a communication interface. The processor, the memory and the input / output interface are connected via a system bus, and the communication interface is connected to the system bus via the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used for three-dimensional VTI medium acoustic wave equations. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal via a network connection. When the computer program is executed by the processor, a three-dimensional anisotropic medium time-space domain seismic forward modeling method is implemented.

[0109] In an exemplary embodiment, a computer device is further provided, including a memory and a processor, wherein a computer program is stored in the memory, and the processor implements the steps in the above-mentioned method embodiments when executing the computer program.

[0110] In an exemplary embodiment, a computer-readable storage medium is provided, storing a computer program, and when the computer program is executed by a processor, the steps in the above method embodiments are implemented.

[0111] In an exemplary embodiment, a computer program product is provided, including a computer program, and when the computer program is executed by a processor, the steps in the above method embodiments are implemented.

[0112] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant regulations.

[0113] Those of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to the memory, database or other medium used in the embodiments provided in the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM may be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0114] The database involved in each embodiment provided in this application may include at least one of a relational database and a non-relational database. The non-relational database may include a distributed database based on blockchain, etc., but is not limited thereto. The processor involved in each embodiment provided in this application may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but is not limited thereto.

[0115] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0116] This article uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, according to the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A three-dimensional anisotropic medium time-space domain seismic forward modeling method, characterized in that: The three-dimensional anisotropic medium time-space domain seismic forward modeling method comprises: Obtain the three-dimensional VTI medium acoustic wave equation; Based on a preset improved time high-order difference operator, the time derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a time discrete function; based on a preset improved space high-order difference operator, the space derivative in the three-dimensional VTI medium acoustic wave equation is discretized to obtain a space discrete function; the time discrete function contains a time term high-order difference coefficient, and the space discrete function contains a space term high-order difference coefficient; Based on the three-dimensional VTI medium acoustic wave equation, simplifying the time discrete function and the space discrete function to obtain a simplified time-space domain dispersion relation; Based on the target optimization algorithm, the simplified spatiotemporal dispersion relation is solved to obtain the optimized time term high-order differential coefficient and the optimized space term high-order differential coefficient; Based on the optimized high-order differential coefficients of the time term, the optimized high-order differential coefficients of the space term, the time discrete function and the space discrete function, the three-dimensional VTI medium acoustic wave equation is extrapolated and solved; Among them, the preset improved time high-order difference operator is: The preset improved spatial high-order difference operator is: in, is the time partial derivative of the pressure wave, is the pressure wave at (0,0,0) at time 1, is the pressure wave at the position (0,0,0) at time 0, Δt is the time sampling interval, v pz is the vertical P-wave velocity, h is the spatial step length, N is half the length of the time difference operator, c m,n,0 and c m,n,l is the high-order differential coefficient of the time term, ε and δ are Thomsen anisotropy parameters, e1, e2, e3, f1, f2, f3 are all mathematical operators, is the time partial derivative of the auxiliary wave, is the auxiliary wave at the position (0,0,0) at time 1, is the auxiliary wave at the position (0,0,0) at time 0; The pressure waves at time 0 are (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l). They are the auxiliary waves at (m,n,0), (-m,n,0), (m,-n,0), (-m,-n,0), (m,0,n), (-m,0,n), (m,0,-n), (-m,0,-n), (0,m,n), (0,-m,n), (0,m,-n), (0,-m,-n), (m,n,l), (m,-n,l), (m,n,-l), (m,-n,-l), (-m,n,l), (-m,-n,-l), (-m,n,-l), (-m,-n,-l), (-m,-n,-l), (-m,-n,-l), and (-m,-n,-l). is the x-axis spatial partial derivative of the pressure wave, is the y-axis spatial partial derivative of the pressure wave, is the z-axis spatial partial derivative of the auxiliary wave, These are the pressure waves at (m,0,0), (-m,0,0), (0,m,0), (0,-m,0), and (0,0,-m) at time 0, respectively. is the auxiliary wave at the position (0,0,m) at time 0, m, n and l are all serial numbers, M is half the length of the spatial term difference operator, c 0,0,0 and c m,0,0 is the high-order difference coefficient of the spatial term.

2. The three-dimensional anisotropic medium time-space domain seismic forward modeling method according to claim 1 is characterized in that: The three-dimensional VTI medium acoustic wave equation is: Among them, P(x, y, z, t) is the pressure wave field and Q(x, y, z, t) is the auxiliary wave field.

3. The three-dimensional anisotropic medium time-space domain seismic forward modeling method according to claim 1, characterized in that: The simplified spatiotemporal dispersion relation is: 2cos(ωΔt)-2=R 2 (A+B+C), Where r = v pz Δt / h, ω is the angular frequency, k x , k y , k z are the spatial wave numbers in different directions respectively.

4. The three-dimensional anisotropic medium time-space domain seismic forward modeling method according to claim 3 is characterized in that: The target-based optimization algorithm solves the simplified spatiotemporal dispersion relation to obtain the optimized time term high-order differential coefficient and the optimized space term high-order differential coefficient, specifically including: The Taylor series is used to expand the trigonometric functions at both ends of the simplified spatiotemporal dispersion relation, and the function formula is as follows: Where β = kh, a 0,0,0 、a m,0,0 、a m,n,0 and a m,n,l are all differential coefficients, and each differential coefficient is c 0,0,0 、c m,0,0 、c m,n,0 and c m,n,l Linear combinations of Arrange the differential coefficients in the expanded function formula into the following preset function equation: Where M1 is the number of differential coefficients, g m is the difference coefficient a 0,0,0 、a m,0,0 、a m,n,0 and a m,n,l θ is the propagation angle of the plane wave, φ is the azimuth of the plane wave, where the plane wave is obtained by decomposing the pressure wave or the auxiliary wave; within different parameter ranges, χ m The specific form of (β,θ,φ) is as follows: When m≤M, When M<m≤M+int(N 2 / 4), when m1=m2, When M<m≤M+int(N 2 / 4), when m1≠m2, When M+int(N 2 / 4)<m≤M1,m3=m4=m5, When M+int(N 2 / 4)<m≤M1,m3=m4≠m5,m4=m5≠m3,m3=m5≠m4, When M+int(N 2 / 4)<m≤M1,m3≠m4≠m5, Taking the minimum error between the two ends of the preset function equation as the goal, constructing an optimization objective function; Solving the optimization objective function to obtain the optimized time term high-order differential coefficient and the optimized space term high-order differential coefficient; Among them, m1, m2, m3, m4, and m5 are all mathematical operators.

5. The three-dimensional anisotropic medium time-space domain seismic forward modeling method according to claim 4, characterized in that: The optimization objective function is: Among them, E is the error value and b is the upper limit value of β.

6. The three-dimensional anisotropic medium time-space domain seismic forward modeling method according to claim 4, characterized in that: Solving the optimization objective function to obtain the optimized time term high-order differential coefficient and the optimized space term high-order differential coefficient specifically includes: solving the optimization objective function by using the least squares optimization method.

7. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the three-dimensional anisotropic medium time-space domain seismic forward modeling method described in any one of claims 1-6.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the three-dimensional anisotropic medium time-space domain seismic forward modeling method described in any one of claims 1-6 is implemented.

9. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the three-dimensional anisotropic medium time-space domain seismic forward modeling method described in any one of claims 1-6 is implemented.

Images